source: trunk/Mars/mcore/zfits.h@ 16840

Last change on this file since 16840 was 16840, checked in by lyard, 12 years ago
bug fix and emplace_back
File size: 18.0 KB
Line 
1/*
2 * zfits.h
3 *
4 * Created on: May 16, 2013
5 * Author: lyard
6 */
7
8#ifndef MARS_ZFITS
9#define MARS_ZFITS
10
11#include <stdexcept>
12
13#include "fits.h"
14#include "huffman.h"
15
16
17#define FACT_RAW 0x0
18#define FACT_SMOOTHING 0x1
19#define FACT_HUFFMAN16 0x2
20
21#define FACT_COL_MAJOR 'C'
22#define FACT_ROW_MAJOR 'R'
23
24
25#ifndef __MARS__
26namespace std
27{
28#endif
29
30class zfits : public fits
31{
32public:
33
34 // Basic constructor
35 zfits(const string& fname, const string& tableName="",
36 bool force=false) : fits(fname, tableName, force),
37 fBuffer(0),
38 fTransposedBuffer(0),
39 fCompressedBuffer(0),
40 fNumTiles(0),
41 fNumRowsPerTile(0),
42 fCurrentRow(-1),
43 fHeapOff(0),
44 fTileSize(0)
45 {
46 InitCompressionReading();
47 }
48
49 // Alternative contstructor
50 zfits(const string& fname, const string& fout, const string& tableName,
51 bool force=false) : fits(fname, fout, tableName, force),
52 fBuffer(0),
53 fTransposedBuffer(0),
54 fCompressedBuffer(0),
55 fNumTiles(0),
56 fNumRowsPerTile(0),
57 fCurrentRow(-1),
58 fHeapOff(0),
59 fTileSize(0)
60 {
61 InitCompressionReading();
62 }
63
64 // Skip the next row
65 bool SkipNextRow()
66 {
67 if (!fTable.isCompressed)
68 return fits::SkipNextRow();
69
70 fRow++;
71 return true;
72 }
73protected:
74
75 // Stage the requested row to internal buffer
76 // Does NOT return data to users
77 virtual void StageRow(size_t row, char* dest)
78 {
79 if (!fTable.isCompressed)
80 {
81 fits::StageRow(row, dest);
82 return;
83 }
84
85 ReadBinaryRow(row, dest);
86 }
87
88private:
89
90 //Structure helper for reading tiles headers
91 typedef struct TileHeader
92 {
93 char id[4];
94 uint32_t numRows;
95 uint64_t size;
96
97 TileHeader() {}
98
99 TileHeader(uint32_t nRows,
100 uint64_t s) : id({'T', 'I', 'L', 'E'}),
101 numRows(nRows),
102 size(s)
103 { };
104 } __attribute__((__packed__)) TileHeader;
105
106 //Structure helper for reading blocks headers and compresion schemes
107 typedef struct BlockHeader
108 {
109 uint64_t size;
110 char ordering;
111 unsigned char numProcs;
112 uint16_t processings[];
113
114 BlockHeader(uint64_t s=0,
115 char o=FACT_ROW_MAJOR,
116 unsigned char n=1) : size(s),
117 ordering(o),
118 numProcs(n)
119 {}
120 } __attribute__((__packed__)) BlockHeader;
121
122 // Do what it takes to initialize the compressed structured
123 void InitCompressionReading()
124 {
125 if (!fTable.isCompressed)
126 return;
127
128 //The constructor may have failed
129 if (!good())
130 return;
131
132 if (fTable.isCompressed)
133 for (auto it=fTable.sortedCols.begin(); it!= fTable.sortedCols.end(); it++)
134 if (it->comp != FACT)
135 {
136#ifdef __EXCEPTIONS
137 throw runtime_error("ERROR: Only the FACT compression scheme is handled by this reader.");
138#else
139 gLog << ___err___ << "ERROR: Only the FACT compression scheme is handled by this reader." << endl;
140 return;
141#endif
142 }
143
144 fColumnOrdering.resize(fTable.sortedCols.size());
145 for (auto it=fColumnOrdering.begin(); it != fColumnOrdering.end(); it++)
146 (*it) = FACT_ROW_MAJOR;
147 //Get compressed specific keywords
148 fNumTiles = fTable.isCompressed ? GetInt("NAXIS2") : 0;
149 fNumRowsPerTile = fTable.isCompressed ? GetInt("ZTILELEN") : 0;
150
151 //give it some space for uncompressing
152 AllocateBuffers();
153
154 //read the file's catalog
155 ReadCatalog();
156
157 //check that heap agrees with head
158 //CheckIfFileIsConsistent();
159 }
160
161 // Copy decompressed data to location requested by user
162 void MoveColumnDataToUserSpace(char* dest, const char* src, const Table::Column& c)
163 {
164 if (!fTable.isCompressed)
165 {
166 fits::MoveColumnDataToUserSpace(dest, src, c);
167 return;
168 }
169
170 memcpy(dest, src, c.num*c.size);
171 }
172
173 vector<char> fBuffer; ///<store the uncompressed rows
174 vector<char> fTransposedBuffer; ///<intermediate buffer to transpose the rows
175 vector<char> fCompressedBuffer; ///<compressed rows
176 vector<char> fColumnOrdering; ///< ordering of the column's rows. Can change from tile to tile.
177
178 size_t fNumTiles; ///< Total number of tiles
179 size_t fNumRowsPerTile; ///< Number of rows per compressed tile
180 size_t fCurrentRow; ///< current row in memory.
181
182 streamoff fHeapOff; ///< offset from the beginning of the file of the binary data
183
184 vector<vector<pair<int64_t, int64_t> > > fCatalog;///< Catalog, i.e. the main table that points to the compressed data.
185 vector<size_t> fTileSize; ///< size in bytes of each compressed tile
186 vector<vector<size_t> > fTileOffsets; ///< offset from start of tile of a given compressed column
187
188 // Get buffer space
189 void AllocateBuffers()
190 {
191 fBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile);
192
193 fTransposedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile);
194 fCompressedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile + fTable.num_cols*(sizeof(BlockHeader)+256)); //use a bit more memory for block headers
195 }
196
197 // Read catalog data. I.e. the address of the compressed data inside the heap
198 void ReadCatalog()
199 {
200 char readBuf[16];
201 fCatalog.resize(fNumTiles);
202
203 const streampos catalogStart = tellg();
204
205 //do the actual reading
206 for (uint32_t i=0;i<fNumTiles;i++)
207 for (uint32_t j=0;j<fTable.num_cols;j++)
208 {
209 read(readBuf, 2*sizeof(int64_t));
210
211 //swap the bytes
212 int64_t tempValues[2] = {0,0};
213 revcpy<8>(reinterpret_cast<char*>(tempValues), readBuf, 2);
214 if (tempValues[0] < 0 || tempValues[1] < 0)
215 {
216 clear(rdstate()|ios::badbit);
217#ifdef __EXCEPTIONS
218 throw runtime_error("Negative value in the catalog");
219#else
220 gLog << ___err___ << "ERROR - negative value in the catalog" << endl;
221 return;
222#endif
223 }
224 //add catalog entry
225 fCatalog[i].emplace_back(tempValues[0], tempValues[1]);
226 }
227
228 //compute the total size of each compressed tile
229 fTileSize.resize(fNumTiles);
230 fTileOffsets.resize(fNumTiles);
231 for (uint32_t i=0;i<fNumTiles;i++)
232 {
233 fTileSize[i] = 0;
234 for (uint32_t j=0;j<fTable.num_cols;j++)
235 {
236 fTileSize[i] += fCatalog[i][j].first;
237 fTileOffsets[i].emplace_back(fCatalog[i][j].second - fCatalog[i][0].second);
238 }
239 }
240 //see if there is a gap before heap data
241 fHeapOff = tellg()+fTable.GetHeapShift();
242
243 if (!fCopy.is_open())
244 return;
245
246 //write catalog and heap gap to target file
247 seekg(catalogStart);
248
249 const size_t catSize = fTable.GetHeapShift() + fTable.total_bytes;
250
251 vector<char> buf(catSize);
252 read(buf.data(), catSize);
253
254 fCopy.write(buf.data(), catSize);
255 if (!fCopy)
256 clear(rdstate()|ios::badbit);
257 }
258
259 //overrides fits.h method with empty one
260 //work is done in ReadBinaryRow because it requires volatile data from ReadBinaryRow
261 virtual void WriteRowToCopyFile(size_t )
262 {
263
264 }
265
266 // Compressed version of the read row
267 bool ReadBinaryRow(const size_t &rowNum, char *bufferToRead)
268 {
269 if (rowNum >= GetNumRows())
270 return false;
271
272 const uint32_t requestedTile = rowNum/fNumRowsPerTile;
273 const uint32_t currentTile = fCurrentRow/fNumRowsPerTile;
274 const size_t previousRow = fCurrentRow;
275
276 fCurrentRow = rowNum;
277
278 //should we read yet another chunk of data ?
279 if (requestedTile != currentTile)
280 {
281 //read yet another chunk from the file
282 const int64_t sizeToRead = fTileSize[requestedTile];
283
284 //skip to the beginning of the tile
285 seekg(fHeapOff+fCatalog[requestedTile][0].second - sizeof(TileHeader));
286 TileHeader tHead;
287 read((char*)(&tHead), sizeof(TileHeader));
288 read(fCompressedBuffer.data(), sizeToRead);
289
290 if (fCurrentRow == previousRow+1 &&
291 fCopy.is_open() &&
292 fCopy.good())
293 {
294 fCopy.write((char*)(&tHead), sizeof(TileHeader));
295 fCopy.write(fCompressedBuffer.data(), sizeToRead);
296 if (!fCopy)
297 clear(rdstate()|ios::badbit);
298 }
299 else
300 if (fCopy.is_open())
301 clear(rdstate()|ios::badbit);
302
303 const uint32_t thisRoundNumRows = (GetNumRows()<fCurrentRow + fNumRowsPerTile) ? GetNumRows()%fNumRowsPerTile : fNumRowsPerTile;
304
305 //uncompress it
306 UncompressBuffer(requestedTile, thisRoundNumRows);
307
308 // pointer to column (source buffer)
309 const char *src = fTransposedBuffer.data();
310
311 uint32_t i=0;
312 for (auto it=fTable.sortedCols.begin(); it!=fTable.sortedCols.end(); it++, i++)
313 {
314 char *buffer = fBuffer.data() + it->offset; // pointer to column (destination buffer)
315
316 switch (fColumnOrdering[i])
317 {
318 case FACT_ROW_MAJOR:
319 // regular, "semi-transposed" copy
320 for (char *dest=buffer; dest<buffer+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row
321 {
322 memcpy(dest, src, it->bytes);
323 src += it->bytes; // next column
324 }
325 break;
326
327 case FACT_COL_MAJOR:
328 // transposed copy
329 for (char *elem=buffer; elem<buffer+it->bytes; elem+=it->size) // element-by-element (arrays)
330 {
331 for (char *dest=elem; dest<elem+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row
332 {
333 memcpy(dest, src, it->size);
334 src += it->size; // next element
335 }
336 }
337 break;
338 default:
339 #ifdef __EXCEPTIONS
340 throw runtime_error("Unkown column ordering scheme");
341 #else
342 gLog << ___err___ << "ERROR - unkown column ordering scheme" << endl;
343 return;
344 #endif
345 break;
346 };
347 }
348 }
349
350 //Data loaded and uncompressed. Copy it to destination
351 memcpy(bufferToRead, fBuffer.data()+fTable.bytes_per_row*(fCurrentRow%fNumRowsPerTile), fTable.bytes_per_row);
352 return good();
353 }
354
355 // Read a bunch of uncompressed data
356 uint32_t UncompressUNCOMPRESSED(char* dest,
357 const char* src,
358 uint32_t numElems,
359 uint32_t sizeOfElems)
360 {
361 memcpy(dest, src, numElems*sizeOfElems);
362 return numElems*sizeOfElems;
363 }
364
365 // Read a bunch of data compressed with the Huffman algorithm
366 uint32_t UncompressHUFFMAN16(char* dest,
367 const char* src,
368 uint32_t numChunks)
369 {
370 vector<uint16_t> uncompressed;
371
372 //read compressed sizes (one per row)
373 const uint32_t* compressedSizes = reinterpret_cast<const uint32_t*>(src);
374 src += sizeof(uint32_t)*numChunks;
375
376 //uncompress the rows, one by one
377 uint32_t sizeWritten = 0;
378 for (uint32_t j=0;j<numChunks;j++)
379 {
380 Huffman::Decode(reinterpret_cast<const unsigned char*>(src), compressedSizes[j], uncompressed);
381
382 memcpy(dest, uncompressed.data(), uncompressed.size()*sizeof(uint16_t));
383
384 sizeWritten += uncompressed.size()*sizeof(uint16_t);
385 dest += uncompressed.size()*sizeof(uint16_t);
386 src += compressedSizes[j];
387 }
388 return sizeWritten;
389 }
390
391 // Apply the inverse transform of the integer smoothing
392 uint32_t UnApplySMOOTHING(int16_t* data,
393 uint32_t numElems)
394 {
395 //un-do the integer smoothing
396 for (uint32_t j=2;j<numElems;j++)
397 data[j] = data[j] + (data[j-1]+data[j-2])/2;
398
399 return numElems*sizeof(uint16_t);
400 }
401
402 // Data has been read from disk. Uncompress it !
403 void UncompressBuffer(const uint32_t &catalogCurrentRow,
404 const uint32_t &thisRoundNumRows)
405 {
406 char *dest = fTransposedBuffer.data();
407
408 //uncompress column by column
409 for (uint32_t i=0; i<fTable.sortedCols.size(); i++)
410 {
411 const fits::Table::Column &col = fTable.sortedCols[i];
412 if (col.num == 0)
413 continue;
414
415 //get the compression flag
416 const int64_t compressedOffset = fTileOffsets[catalogCurrentRow][i];
417
418 BlockHeader* head = reinterpret_cast<BlockHeader*>(&fCompressedBuffer[compressedOffset]);
419
420 fColumnOrdering[i] = head->ordering;
421
422 uint32_t numRows = (head->ordering==FACT_ROW_MAJOR) ? thisRoundNumRows : col.num;
423 uint32_t numCols = (head->ordering==FACT_COL_MAJOR) ? thisRoundNumRows : col.num;
424
425 const char *src = fCompressedBuffer.data()+compressedOffset+sizeof(BlockHeader)+sizeof(uint16_t)*head->numProcs;
426
427 for (int32_t j=head->numProcs-1;j >= 0; j--)
428 {
429 uint32_t sizeWritten=0;
430
431 switch (head->processings[j])
432 {
433 case FACT_RAW:
434 sizeWritten = UncompressUNCOMPRESSED(dest, src, numRows*numCols, col.size);
435 break;
436 case FACT_SMOOTHING:
437 sizeWritten = UnApplySMOOTHING(reinterpret_cast<int16_t*>(dest), numRows*numCols);
438 break;
439 case FACT_HUFFMAN16:
440 sizeWritten = UncompressHUFFMAN16(dest, src, numRows);
441 break;
442 default:
443#ifdef __EXCEPTIONS
444 throw runtime_error("Unknown processing applied to data. Aborting");
445#else
446 gLog << ___err___ << "ERROR - Unknown processing applied to data. Aborting" << endl;
447 return;
448#endif
449 break;
450 }
451 //increment destination counter only when processing done.
452 if (j==0) dest+= sizeWritten;
453 }
454 }
455 }
456
457 void CheckIfFileIsConsistent()
458 {
459 //goto start of heap
460 streamoff whereAreWe = tellg();
461 seekg(fHeapOff);
462
463 //init number of rows to zero
464 uint64_t numRows = 0;
465
466 //get number of columns from header
467 size_t numCols = fTable.num_cols;
468
469 vector<vector<pair<int64_t, int64_t> > > catalog;
470
471 TileHeader tileHead;
472 BlockHeader columnHead;
473 streamoff offsetInHeap = 0;
474 //skip through the heap
475 while (true)
476 {
477 read((char*)(&tileHead), sizeof(TileHeader));
478 //end of file
479 if (!good())
480 break;
481 //padding or corrupt data
482 if (tileHead.id[0] != 'T' || tileHead.id[1] != 'I' ||
483 tileHead.id[2] != 'L' || tileHead.id[3] != 'E')
484 break;
485
486 //a new tile begins here
487 catalog.push_back(vector<pair<int64_t, int64_t> >(0));
488 offsetInHeap += sizeof(TileHeader);
489
490 //skip through the columns
491 for (size_t i=0;i<numCols;i++)
492 {
493 //zero sized column do not have headers. Skip it
494 if (fTable.sortedCols[i].num == 0)
495 {
496 catalog.back().push_back(make_pair(0,0));
497 continue;
498 }
499 //read column header
500 read((char*)(&columnHead), sizeof(BlockHeader));
501 //corrupted tile
502 if (!good())
503 break;
504 catalog.back().emplace_back(make_pair((int64_t)(columnHead.size),offsetInHeap));
505 offsetInHeap += columnHead.size;
506 seekg(fHeapOff+offsetInHeap);
507 }
508
509 //if we ain't good, this means that something went wrong inside the current tile.
510 if (!good())
511 {
512 catalog.pop_back();
513 break;
514 }
515
516 //current tile is complete. Add rows
517 numRows += tileHead.numRows;
518 }
519
520 if (catalog.size() != fCatalog.size() ||
521 numRows != fTable.num_rows)
522 {
523#ifdef __EXCEPTIONS
524 throw runtime_error("Heap data does not agree with header. Aborting");
525#else
526 gLog << ___err___ << "ERROR - Heap data does not agree with header." << endl;
527 return;
528#endif
529 }
530
531 for (uint32_t i=0;i<catalog.size(); i++)
532 for (uint32_t j=0;j<numCols;j++)
533 {
534 if (catalog[i][j].first != fCatalog[i][j].first ||
535 catalog[i][j].second != fCatalog[i][j].second)
536 {
537#ifdef __EXCEPTIONS
538 throw runtime_error("Heap data does not agree with header. Aborting");
539#else
540 gLog << ___err___ << "ERROR - Heap data does not agree with header." << endl;
541 return;
542#endif
543 }
544 }
545 //go back to start of heap
546 seekg(whereAreWe);
547 }
548
549};//class zfits
550
551#ifndef __MARS__
552}; //namespace std
553#endif
554
555#endif
Note: See TracBrowser for help on using the repository browser.