CCfits
2.7
|
00001 // Astrophysics Science Division, 00002 // NASA/ Goddard Space Flight Center 00003 // HEASARC 00004 // http://heasarc.gsfc.nasa.gov 00005 // e-mail: ccfits@legacy.gsfc.nasa.gov 00006 // 00007 // Original author: Ben Dorman 00008 00009 #ifndef COLUMNVECTORDATA_H 00010 #define COLUMNVECTORDATA_H 1 00011 #ifdef _MSC_VER 00012 #include "MSconfig.h" 00013 #endif 00014 #include "CCfits.h" 00015 00016 // valarray 00017 #include <valarray> 00018 // vector 00019 #include <vector> 00020 // Column 00021 #include "Column.h" 00022 00023 #ifdef SSTREAM_DEFECT 00024 #include <strstream> 00025 #else 00026 #include <sstream> 00027 #endif 00028 00029 #include <memory> 00030 #include <numeric> 00031 #include <algorithm> 00032 00033 namespace CCfits { 00034 00035 class Table; 00036 00037 } 00038 00039 #include "FITS.h" 00040 #include "FITSUtil.h" 00041 using std::complex; 00042 00043 00044 namespace CCfits { 00045 00046 00047 00048 template <typename T> 00049 class ColumnVectorData : public Column //## Inherits: <unnamed>%38BAD1D4D370 00050 { 00051 00052 public: 00053 ColumnVectorData(const ColumnVectorData< T > &right); 00054 ColumnVectorData (Table* p = 0); 00055 ColumnVectorData (int columnIndex, const string &columnName, ValueType type, const string &format, const string &unit, Table* p, int rpt = 1, long w = 1, const string &comment = ""); 00056 ~ColumnVectorData(); 00057 00058 virtual void readData (long firstrow, long nelements, long firstelem = 1); 00059 virtual ColumnVectorData<T>* clone () const; 00060 virtual void setDimen (); 00061 void setDataLimits (T* limits); 00062 const T minLegalValue () const; 00063 void minLegalValue (T value); 00064 const T maxLegalValue () const; 00065 void maxLegalValue (T value); 00066 const T minDataValue () const; 00067 void minDataValue (T value); 00068 const T maxDataValue () const; 00069 void maxDataValue (T value); 00070 const std::vector<std::valarray<T> >& data () const; 00071 void setData (const std::vector<std::valarray<T> >& value); 00072 const std::valarray<T>& data (int i) const; 00073 void data (int i, const std::valarray<T>& value); 00074 00075 // Additional Public Declarations 00076 friend class Column; 00077 protected: 00078 // Additional Protected Declarations 00079 00080 private: 00081 ColumnVectorData< T > & operator=(const ColumnVectorData< T > &right); 00082 00083 virtual bool compare (const Column &right) const; 00084 void resizeDataObject (const std::vector<std::valarray<T> >& indata, size_t firstRow); 00085 // Reads a specified number of column rows. 00086 // 00087 // There are no default arguments. The function 00088 // Column::read(firstrow,firstelem,nelements) 00089 // is designed for reading the whole column. 00090 virtual void readColumnData (long first, long last, T* nullValue = 0); 00091 virtual std::ostream& put (std::ostream& s) const; 00092 void writeData (const std::valarray<T>& indata, long numRows, long firstRow = 1, T* nullValue = 0); 00093 void writeData (const std::vector<std::valarray<T> >& indata, long firstRow = 1, T* nullValue = 0); 00094 // Reads a specified number of column rows. 00095 // 00096 // There are no default arguments. The function 00097 // Column::read(firstrow,firstelem,nelements) 00098 // is designed for reading the whole column. 00099 virtual void readRow (size_t row, T* nullValue = 0); 00100 // Reads a variable row.. 00101 virtual void readVariableRow (size_t row, T* nullValue = 0); 00102 void readColumnData (long firstrow, long nelements, long firstelem, T* nullValue = 0); 00103 void writeData (const std::valarray<T>& indata, const std::vector<long>& vectorLengths, long firstRow = 1, T* nullValue = 0); 00104 void writeFixedRow (const std::valarray<T>& data, long row, long firstElem = 1, T* nullValue = 0); 00105 void writeFixedArray (T* data, long nElements, long nRows, long firstRow, T* nullValue = 0); 00106 // Insert one or more blank rows into a FITS column. 00107 virtual void insertRows (long first, long number = 1); 00108 virtual void deleteRows (long first, long number = 1); 00109 virtual size_t getStoredDataSize() const; 00110 void doWrite (T* array, long row, long rowSize, long firstElem, T* nullValue); 00111 00112 // Additional Private Declarations 00113 00114 private: //## implementation 00115 // Data Members for Class Attributes 00116 T m_minLegalValue; 00117 T m_maxLegalValue; 00118 T m_minDataValue; 00119 T m_maxDataValue; 00120 00121 // Data Members for Associations 00122 std::vector<std::valarray<T> > m_data; 00123 00124 // Additional Implementation Declarations 00125 00126 }; 00127 00128 // Parameterized Class CCfits::ColumnVectorData 00129 00130 template <typename T> 00131 inline void ColumnVectorData<T>::readData (long firstrow, long nelements, long firstelem) 00132 { 00133 readColumnData(firstrow,nelements,firstelem,static_cast<T*>(0)); 00134 } 00135 00136 template <typename T> 00137 inline const T ColumnVectorData<T>::minLegalValue () const 00138 { 00139 return m_minLegalValue; 00140 } 00141 00142 template <typename T> 00143 inline void ColumnVectorData<T>::minLegalValue (T value) 00144 { 00145 m_minLegalValue = value; 00146 } 00147 00148 template <typename T> 00149 inline const T ColumnVectorData<T>::maxLegalValue () const 00150 { 00151 return m_maxLegalValue; 00152 } 00153 00154 template <typename T> 00155 inline void ColumnVectorData<T>::maxLegalValue (T value) 00156 { 00157 m_maxLegalValue = value; 00158 } 00159 00160 template <typename T> 00161 inline const T ColumnVectorData<T>::minDataValue () const 00162 { 00163 return m_minDataValue; 00164 } 00165 00166 template <typename T> 00167 inline void ColumnVectorData<T>::minDataValue (T value) 00168 { 00169 m_minDataValue = value; 00170 } 00171 00172 template <typename T> 00173 inline const T ColumnVectorData<T>::maxDataValue () const 00174 { 00175 return m_maxDataValue; 00176 } 00177 00178 template <typename T> 00179 inline void ColumnVectorData<T>::maxDataValue (T value) 00180 { 00181 m_maxDataValue = value; 00182 } 00183 00184 template <typename T> 00185 inline const std::vector<std::valarray<T> >& ColumnVectorData<T>::data () const 00186 { 00187 return m_data; 00188 } 00189 00190 template <typename T> 00191 inline void ColumnVectorData<T>::setData (const std::vector<std::valarray<T> >& value) 00192 { 00193 m_data = value; 00194 } 00195 00196 template <typename T> 00197 inline const std::valarray<T>& ColumnVectorData<T>::data (int i) const 00198 { 00199 return m_data[i - 1]; 00200 } 00201 00202 template <typename T> 00203 inline void ColumnVectorData<T>::data (int i, const std::valarray<T>& value) 00204 { 00205 if (m_data[i-1].size() != value.size()) 00206 m_data[i-1].resize(value.size()); 00207 m_data[i - 1] = value; 00208 } 00209 00210 // Parameterized Class CCfits::ColumnVectorData 00211 00212 template <typename T> 00213 ColumnVectorData<T>::ColumnVectorData(const ColumnVectorData<T> &right) 00214 :Column(right), 00215 m_minLegalValue(right.m_minLegalValue), 00216 m_maxLegalValue(right.m_maxLegalValue), 00217 m_minDataValue(right.m_minDataValue), 00218 m_maxDataValue(right.m_maxDataValue), 00219 m_data(right.m_data) 00220 { 00221 } 00222 00223 template <typename T> 00224 ColumnVectorData<T>::ColumnVectorData (Table* p) 00225 : Column(p), 00226 m_minLegalValue(0), 00227 m_maxLegalValue(0), 00228 m_minDataValue(0), 00229 m_maxDataValue(0), 00230 m_data() 00231 { 00232 } 00233 00234 template <typename T> 00235 ColumnVectorData<T>::ColumnVectorData (int columnIndex, const string &columnName, ValueType type, const string &format, const string &unit, Table* p, int rpt, long w, const string &comment) 00236 : Column(columnIndex,columnName,type,format,unit,p,rpt,w,comment), 00237 m_minLegalValue(0), 00238 m_maxLegalValue(0), 00239 m_minDataValue(0), 00240 m_maxDataValue(0), 00241 m_data() 00242 { 00243 } 00244 00245 00246 template <typename T> 00247 ColumnVectorData<T>::~ColumnVectorData() 00248 { 00249 // valarray destructor should do all the work. 00250 } 00251 00252 00253 template <typename T> 00254 bool ColumnVectorData<T>::compare (const Column &right) const 00255 { 00256 if ( !Column::compare(right) ) return false; 00257 const ColumnVectorData<T>& that = static_cast<const ColumnVectorData<T>&>(right); 00258 size_t n = m_data.size(); 00259 // m_data is of type valarray<T>. 00260 if ( that.m_data.size() != n ) return false; 00261 for (size_t i = 0; i < n ; i++) 00262 { 00263 const std::valarray<T>& thisValArray=m_data[i]; 00264 const std::valarray<T>& thatValArray=that.m_data[i]; 00265 size_t nn = thisValArray.size(); 00266 if (thatValArray.size() != nn ) return false; 00267 00268 for (size_t j = 0; j < nn ; j++ ) 00269 { 00270 if (thisValArray[j] != thatValArray[j]) 00271 return false; 00272 } 00273 } 00274 return true; 00275 } 00276 00277 template <typename T> 00278 ColumnVectorData<T>* ColumnVectorData<T>::clone () const 00279 { 00280 return new ColumnVectorData<T>(*this); 00281 } 00282 00283 template <typename T> 00284 void ColumnVectorData<T>::resizeDataObject (const std::vector<std::valarray<T> >& indata, size_t firstRow) 00285 { 00286 // the rows() call is the value before updating. 00287 // the updateRows() call at the end sets the call to return the 00288 // value from the fits pointer - which is changed by writeFixedArray 00289 // or writeFixedRow. 00290 00291 const size_t lastInputRow(indata.size() + firstRow - 1); 00292 const size_t newLastRow = std::max(lastInputRow,static_cast<size_t>(rows())); 00293 00294 // if the write instruction increases the rows, we need to add 00295 // rows to the data member and preserve its current contents. 00296 00297 // rows() >= origNRows since it is the value for entire table, 00298 // not just this column. 00299 const size_t origNRows(m_data.size()); 00300 // This will always be an expansion. vector.resize() doesn't 00301 // invalidate any data on an expansion. 00302 if (newLastRow > origNRows) m_data.resize(newLastRow); 00303 00304 if (varLength()) 00305 { 00306 // The incoming data will determine each row size, thus 00307 // no need to preserve any existing values in the row. 00308 // Each value will eventually be overwritten. 00309 for (size_t iRow = firstRow-1; iRow < lastInputRow; ++iRow) 00310 { 00311 std::valarray<T>& current = m_data[iRow]; 00312 const size_t newSize = indata[iRow - (firstRow-1)].size(); 00313 if (current.size() != newSize) 00314 current.resize(newSize); 00315 } 00316 } 00317 else 00318 { 00319 // All row sizes in m_data should ALWAYS be either repeat(), 00320 // or 0 if they haven't been initialized. This is true regardless 00321 // of the incoming data row size. 00322 00323 // Perform LAZY initialization of m_data rows. Only 00324 // expand a row valarray when it is first needed. 00325 for (size_t iRow = firstRow-1; iRow < lastInputRow; ++iRow) 00326 { 00327 if (m_data[iRow].size() != repeat()) 00328 m_data[iRow].resize(repeat()); 00329 } 00330 } 00331 } 00332 00333 template <typename T> 00334 void ColumnVectorData<T>::setDimen () 00335 { 00336 int status(0); 00337 FITSUtil:: auto_array_ptr<char> dimValue (new char[FLEN_VALUE]); 00338 00339 #ifdef SSTREAM_DEFECT 00340 std::ostrstream key; 00341 #else 00342 std::ostringstream key; 00343 #endif 00344 key << "TDIM" << index(); 00345 00346 #ifdef SSTREAM_DEFECT 00347 fits_read_key_str(fitsPointer(), key.str(), dimValue.get(),0,&status); 00348 #else 00349 fits_read_key_str(fitsPointer(),const_cast<char*>(key.str().c_str()),dimValue.get(),0,&status); 00350 #endif 00351 00352 if (status == 0) 00353 { 00354 dimen(String(dimValue.get())); 00355 } 00356 } 00357 00358 template <typename T> 00359 void ColumnVectorData<T>::readColumnData (long first, long last, T* nullValue) 00360 { 00361 makeHDUCurrent(); 00362 00363 00364 if ( rows() < last ) 00365 { 00366 std::cerr << "CCfits: More data requested than contained in table. "; 00367 std::cerr << "Extracting complete column.\n"; 00368 last = rows(); 00369 } 00370 00371 long nelements = (last - first + 1)*repeat(); 00372 00373 00374 readColumnData(first,nelements,1,nullValue); 00375 if (first <= 1 && last == rows()) isRead(true); 00376 } 00377 00378 template <typename T> 00379 std::ostream& ColumnVectorData<T>::put (std::ostream& s) const 00380 { 00381 // output header information 00382 Column::put(s); 00383 if ( FITS::verboseMode() ) 00384 { 00385 s << " Column Legal limits: ( " << m_minLegalValue << "," << m_maxLegalValue << " )\n" 00386 << " Column Data limits: ( " << m_minDataValue << "," << m_maxDataValue << " )\n"; 00387 } 00388 if (!m_data.empty()) 00389 { 00390 for (size_t j = 0; j < m_data.size(); j++) 00391 { 00392 size_t n = m_data[j].size(); 00393 if ( n ) 00394 { 00395 s << "Row " << j + 1 << " Vector Size " << n << '\n'; 00396 for (size_t k = 0; k < n - 1; k++) 00397 { 00398 s << m_data[j][k] << '\t'; 00399 } 00400 s << m_data[j][n - 1] << '\n'; 00401 } 00402 } 00403 } 00404 00405 return s; 00406 } 00407 00408 template <typename T> 00409 void ColumnVectorData<T>::writeData (const std::valarray<T>& indata, long numRows, long firstRow, T* nullValue) 00410 { 00411 // This version of writeData is called by Column write functions which 00412 // can only write the same number of elements to each row. 00413 // For fixed width columns, this must be equal to the repeat value 00414 // or an exception is thrown. For variable width, it only requires 00415 // that indata.size()/numRows is an int. 00416 00417 // won't do anything if < 0, and will give divide check if == 0. 00418 if (numRows <= 0) throw InvalidNumberOfRows(numRows); 00419 00420 #ifdef SSTREAM_DEFECT 00421 std::ostrstream msgStr; 00422 #else 00423 std::ostringstream msgStr; 00424 #endif 00425 if (indata.size() % static_cast<size_t>(numRows)) 00426 { 00427 msgStr << "To use this write function, input array size" 00428 <<"\n must be exactly divisible by requested num rows: " 00429 << numRows; 00430 throw InsufficientElements(msgStr.str()); 00431 } 00432 const size_t cellsize = indata.size()/static_cast<size_t>(numRows); 00433 00434 if (!varLength() && cellsize != repeat() ) 00435 { 00436 msgStr << "column: " << name() 00437 << "\n input data size: " << indata.size() 00438 << " required: " << numRows*repeat(); 00439 String msg(msgStr.str()); 00440 throw InsufficientElements(msg); 00441 } 00442 00443 std::vector<std::valarray<T> > internalFormat(numRows); 00444 00445 // support writing equal row lengths to variable columns. 00446 00447 for (long j = 0; j < numRows; ++j) 00448 { 00449 internalFormat[j].resize(cellsize); 00450 internalFormat[j] = indata[std::slice(cellsize*j,cellsize,1)]; 00451 } 00452 00453 // change the size of m_data based on the first row to be written 00454 // and on the input data vector sizes. 00455 00456 writeData(internalFormat,firstRow,nullValue); 00457 } 00458 00459 template <typename T> 00460 void ColumnVectorData<T>::writeData (const std::vector<std::valarray<T> >& indata, long firstRow, T* nullValue) 00461 { 00462 // This is called directly by Column's writeArrays functions, and indirectly 00463 // by both categories of write functions, ie. those which allow differing 00464 // lengths per row and those that don't. 00465 const size_t nInputRows(indata.size()); 00466 using std::valarray; 00467 00468 resizeDataObject(indata,firstRow); 00469 // After the above call, can assume all m_data arrays to be written to 00470 // have been properly resized whether we're dealing with fixed or 00471 // variable length. 00472 00473 if (varLength()) 00474 { 00475 // firstRow is 1-based, but all these internal row variables 00476 // will be 0-based. 00477 const size_t endRow = nInputRows + firstRow-1; 00478 for (size_t iRow = firstRow-1; iRow < endRow; ++iRow) 00479 { 00480 m_data[iRow] = indata[iRow - (firstRow-1)]; 00481 // doWrite wants 1-based rows. 00482 doWrite(&m_data[iRow][0], iRow+1, m_data[iRow].size(), 1, nullValue); 00483 } 00484 parent()->updateRows(); 00485 } 00486 else 00487 { 00488 // Check for simplest case of all valarrays of size repeat(). 00489 // If any are greater, throw an error. 00490 const size_t colRepeat = repeat(); 00491 bool allEqualRepeat = true; 00492 for (size_t i=0; i<nInputRows; ++i) 00493 { 00494 const size_t sz = indata[i].size(); 00495 if (sz > colRepeat) 00496 { 00497 #ifdef SSTREAM_DEFECT 00498 std::ostrstream oss; 00499 #else 00500 std::ostringstream oss; 00501 #endif 00502 oss << " vector column length " << colRepeat 00503 <<", input valarray length " << sz; 00504 throw InvalidRowParameter(oss.str()); 00505 } 00506 if (sz < colRepeat) 00507 allEqualRepeat = false; 00508 } 00509 00510 if (allEqualRepeat) 00511 { 00512 // concatenate the valarrays and write. 00513 const size_t nElements (colRepeat*nInputRows); 00514 FITSUtil::CVAarray<T> convert; 00515 FITSUtil::auto_array_ptr<T> pArray(convert(indata)); 00516 T* array = pArray.get(); 00517 00518 // if T is complex, then CVAarray returns a 00519 // C-array of complex objects. But FITS requires an array of complex's 00520 // value_type. 00521 00522 // This writes to the file and also calls updateRows. 00523 writeFixedArray(array,nElements,nInputRows,firstRow,nullValue); 00524 00525 for (size_t j = 0; j < nInputRows ; ++j) 00526 { 00527 const valarray<T>& input = indata[j]; 00528 valarray<T>& current = m_data[j + firstRow - 1]; 00529 // current should be resized by resizeDataObject. 00530 current = input; 00531 } 00532 } 00533 else 00534 { 00535 // Some input arrays have fewer than colRepeat elements. 00536 const size_t endRow = nInputRows + firstRow-1; 00537 for (size_t iRow = firstRow-1; iRow<endRow; ++iRow) 00538 { 00539 // resizeDataObject should already have resized all 00540 // corresponding m_data rows to repeat(). 00541 const valarray<T>& input = indata[iRow-(firstRow-1)]; 00542 writeFixedRow(input, iRow, 1, nullValue); 00543 } 00544 parent()->updateRows(); 00545 } 00546 00547 } // end if !varLength 00548 } 00549 00550 template <typename T> 00551 void ColumnVectorData<T>::readRow (size_t row, T* nullValue) 00552 { 00553 makeHDUCurrent(); 00554 00555 00556 00557 if ( row > static_cast<size_t>(rows()) ) 00558 { 00559 #ifdef SSTREAM_DEFECT 00560 std::ostrstream msg; 00561 #else 00562 std::ostringstream msg; 00563 #endif 00564 msg << " row requested: " << row << " row range: 1 - " << rows(); 00565 #ifdef SSTREAM_DEFECT 00566 msg << std::ends; 00567 #endif 00568 00569 throw Column::InvalidRowNumber(msg.str()); 00570 } 00571 00572 // this is really for documentation purposes. I expect the optimizer will 00573 // remove this redundant definition . 00574 bool variable(type() < 0); 00575 00576 00577 long nelements(repeat()); 00578 00579 if (variable) 00580 { 00581 readVariableRow(row,nullValue); 00582 } 00583 else 00584 { 00585 readColumnData(row,nelements,1,nullValue); 00586 } 00587 } 00588 00589 template <typename T> 00590 void ColumnVectorData<T>::readVariableRow (size_t row, T* nullValue) 00591 { 00592 int status(0); 00593 long offset(0); 00594 long repeat(0); 00595 if (fits_read_descript(fitsPointer(),index(),static_cast<long>(row), 00596 &repeat,&offset,&status)) throw FitsError(status); 00597 readColumnData(row,repeat,1,nullValue); 00598 } 00599 00600 template <typename T> 00601 void ColumnVectorData<T>::readColumnData (long firstrow, long nelements, long firstelem, T* nullValue) 00602 { 00603 int status=0; 00604 00605 FITSUtil::auto_array_ptr<T> pArray(new T[nelements]); 00606 T* array = pArray.get(); 00607 int anynul(0); 00608 00609 00610 00611 if (fits_read_col(fitsPointer(), abs(type()),index(), firstrow, firstelem, 00612 nelements, nullValue, array, &anynul, &status) != 0) 00613 throw FitsError(status); 00614 00615 size_t countRead = 0; 00616 const size_t ONE = 1; 00617 00618 if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows()); 00619 size_t vectorSize(0); 00620 if (!varLength()) 00621 { 00622 00623 vectorSize = std::max(repeat(),ONE); // safety check. 00624 00625 } 00626 else 00627 { 00628 // assume that the user specified the correct length for 00629 // variable columns. This should be ok since readVariableColumns 00630 // uses fits_read_descripts to return this information from the 00631 // fits pointer, and this is passed as nelements here. 00632 vectorSize = nelements; 00633 } 00634 size_t n = nelements; 00635 00636 int i = firstrow; 00637 int ii = i - 1; 00638 while ( countRead < n) 00639 { 00640 std::valarray<T>& current = m_data[ii]; 00641 if (current.size() != vectorSize) current.resize(vectorSize); 00642 int elementsInFirstRow = vectorSize-firstelem + 1; 00643 bool lastRow = ( (nelements - countRead) < vectorSize); 00644 if (lastRow) 00645 { 00646 int elementsInLastRow = nelements - countRead; 00647 std::valarray<T> ttmp(array + vectorSize*(ii-firstrow) + elementsInFirstRow, 00648 elementsInLastRow); 00649 for (int kk = 0; kk < elementsInLastRow; kk++) current[kk] = ttmp[kk]; 00650 countRead += elementsInLastRow; 00651 00652 } 00653 // what to do with complete rows 00654 else 00655 { 00656 if (firstelem == 1 || (firstelem > 1 && i > firstrow) ) 00657 { 00658 std::valarray<T> ttmp(array + vectorSize*(ii - firstrow) + 00659 elementsInFirstRow,vectorSize); 00660 current = ttmp; 00661 ii++; 00662 i++; 00663 countRead += vectorSize; 00664 } 00665 else 00666 { 00667 if (i == firstrow) 00668 { 00669 std::valarray<T> ttmp(array,elementsInFirstRow); 00670 for (size_t kk = firstelem ; kk < vectorSize ; kk++) 00671 current[kk] = ttmp[kk-firstelem]; 00672 countRead += elementsInFirstRow; 00673 i++; 00674 ii++; 00675 } 00676 } 00677 } 00678 } 00679 } 00680 00681 template <typename T> 00682 void ColumnVectorData<T>::writeData (const std::valarray<T>& indata, const std::vector<long>& vectorLengths, long firstRow, T* nullValue) 00683 { 00684 // Called from Column write functions which allow differing lengths 00685 // for each row. 00686 using namespace std; 00687 const size_t N(vectorLengths.size()); 00688 vector<long> sums(N); 00689 // pre-calculate partial sums of vector lengths for use as array offsets. 00690 partial_sum(vectorLengths.begin(),vectorLengths.end(),sums.begin()); 00691 // check that sufficient data have been supplied to carry out the entire operation. 00692 if (indata.size() < static_cast<size_t>(sums[N-1]) ) 00693 { 00694 #ifdef SSTREAM_DEFECT 00695 ostrstream msgStr; 00696 #else 00697 ostringstream msgStr; 00698 #endif 00699 msgStr << " input data size: " << indata.size() << " vector length sum: " << sums[N-1]; 00700 #ifdef SSTREAM_DEFECT 00701 msgStr << std::ends; 00702 #endif 00703 00704 String msg(msgStr.str()); 00705 throw InsufficientElements(msg); 00706 } 00707 00708 vector<valarray<T> > vvArray(N); 00709 long& last = sums[0]; 00710 vvArray[0].resize(last); 00711 for (long jj = 0; jj < last; ++jj) vvArray[0][jj] = indata[jj]; 00712 00713 for (size_t j = 1; j < N; ++j) 00714 { 00715 valarray<T>& __tmp = vvArray[j]; 00716 // these make the code much more readable 00717 long& first = sums[j-1]; 00718 long& jlast = sums[j]; 00719 __tmp.resize(jlast - first); 00720 for (long k = first; k < jlast; ++k) 00721 { 00722 __tmp[k - first] = indata[k]; 00723 } 00724 } 00725 00726 writeData(vvArray,firstRow,nullValue); 00727 } 00728 00729 template <typename T> 00730 void ColumnVectorData<T>::writeFixedRow (const std::valarray<T>& data, long row, long firstElem, T* nullValue) 00731 { 00732 00733 // This is to be called only for FIXED length vector columns. It will 00734 // throw if data.size()+firstElem goes beyond the repeat value. 00735 // If data.size() is less than repeat, it leaves the remaining values 00736 // undisturbed both in the file and in m_data storage. 00737 00738 #ifdef SSTREAM_DEFECT 00739 std::ostrstream msgStr; 00740 #else 00741 std::ostringstream msgStr; 00742 #endif 00743 if (varLength()) 00744 { 00745 msgStr <<"Calling ColumnVectorData::writeFixedRow for a variable length column.\n"; 00746 throw FitsFatal(msgStr.str()); 00747 } 00748 00749 std::valarray<T>& storedRow = m_data[row]; 00750 long inputSize = static_cast<long>(data.size()); 00751 long storedSize(storedRow.size()); 00752 if (storedSize != static_cast<long>(repeat())) 00753 { 00754 msgStr<<"stored array size vs. column width mismatch in ColumnVectorData::writeFixedRow.\n"; 00755 throw FitsFatal(msgStr.str()); 00756 } 00757 00758 if (inputSize + firstElem - 1 > storedSize) 00759 { 00760 msgStr << " requested write " << firstElem << " to " 00761 << firstElem + inputSize - 1 << " exceeds vector length " << repeat(); 00762 throw InvalidRowParameter(msgStr.str()); 00763 } 00764 00765 // CANNOT give a strong exception safety guarantee because writing 00766 // data changes the file. Any corrective action that could be taken 00767 // [e.g. holding initial contents of the row and writing it back after 00768 // an exception is thrown] could in principle throw the same exception 00769 // we are trying to protect from. 00770 00771 // routine does however give the weak guarantee (no resource leaks). 00772 00773 // It's never a good thing to cast away a const, but doWrite calls the 00774 // CFITSIO write functions which take a non-const pointer (though 00775 // it shouldn't actually modify the array), and I'd rather not 00776 // copy the entire valarray just to avoid this problem. 00777 std::valarray<T>& lvData = const_cast<std::valarray<T>&>(data); 00778 T* inPointer = &lvData[0]; 00779 doWrite(inPointer, row+1, inputSize, firstElem, nullValue); 00780 00781 // Writing to disk was successful, now update FITS object and return. 00782 const size_t offset = static_cast<size_t>(firstElem) - 1; 00783 for (size_t iElem=0; iElem < static_cast<size_t>(inputSize); ++iElem) 00784 { 00785 // This doesn't require inPointer's non-constness. It's just 00786 // used here to speed things up a bit. 00787 storedRow[iElem + offset] = inPointer[iElem]; 00788 } 00789 } 00790 00791 template <typename T> 00792 void ColumnVectorData<T>::writeFixedArray (T* data, long nElements, long nRows, long firstRow, T* nullValue) 00793 { 00794 int status(0); 00795 00796 // check for sanity of inputs, then write to file. 00797 // this function writes only complete rows to a table with 00798 // fixed width rows. 00799 00800 00801 if ( nElements < nRows*static_cast<long>(repeat()) ) 00802 { 00803 #ifdef SSTREAM_DEFECT 00804 std::ostrstream msgStr; 00805 #else 00806 std::ostringstream msgStr; 00807 #endif 00808 msgStr << " input array size: " << nElements << " required " << nRows*repeat(); 00809 String msg(msgStr.str()); 00810 00811 throw Column::InsufficientElements(msg); 00812 } 00813 00814 if (nullValue) 00815 { 00816 if (fits_write_colnull(fitsPointer(),abs(type()),index(),firstRow, 00817 1,nElements,data,nullValue,&status)) throw FitsError(status); 00818 } 00819 else 00820 { 00821 if (fits_write_col(fitsPointer(),abs(type()),index(),firstRow, 00822 1,nElements,data,&status)) throw FitsError(status); 00823 } 00824 00825 parent()->updateRows(); 00826 } 00827 00828 template <typename T> 00829 void ColumnVectorData<T>::insertRows (long first, long number) 00830 { 00831 if (first >= 0 && first <= static_cast<long>(m_data.size())) 00832 { 00833 typename std::vector<std::valarray<T> >::iterator in; 00834 if (first !=0) 00835 { 00836 in = m_data.begin()+first; 00837 } 00838 else 00839 { 00840 in = m_data.begin(); 00841 } 00842 00843 // non-throwing operations. 00844 m_data.insert(in,number,std::valarray<T>(T(),0)); 00845 } 00846 } 00847 00848 template <typename T> 00849 void ColumnVectorData<T>::deleteRows (long first, long number) 00850 { 00851 // Don't assume the calling routine (ie. Table's deleteRows) 00852 // knows Column's current m_data size. m_data may still be 00853 // size 0 if no read operations have been performed on Column. 00854 // Therefore perform range checking before erasing. 00855 const long curSize = static_cast<long>(m_data.size()); 00856 if (curSize>0 && first <= curSize) 00857 { 00858 const long last = std::min(curSize, first-1+number); 00859 m_data.erase(m_data.begin()+first-1,m_data.begin()+last); 00860 } 00861 00862 } 00863 00864 template <typename T> 00865 size_t ColumnVectorData<T>::getStoredDataSize() const 00866 { 00867 return m_data.size(); 00868 } 00869 00870 template <typename T> 00871 void ColumnVectorData<T>::setDataLimits (T* limits) 00872 { 00873 m_minLegalValue = limits[0]; 00874 m_maxLegalValue = limits[1]; 00875 m_minDataValue = std::max(limits[2],limits[0]); 00876 m_maxDataValue = std::min(limits[3],limits[1]); 00877 } 00878 00879 template <typename T> 00880 void ColumnVectorData<T>::doWrite (T* array, long row, long rowSize, long firstElem, T* nullValue) 00881 { 00882 int status(0); 00883 // internal functioning of write_colnull forbids its use for writing 00884 // variable width columns. If a nullvalue argument was supplied it will 00885 // be ignored. 00886 if ( !varLength()) 00887 { 00888 if (fits_write_colnull(fitsPointer(),type(),index(),row, firstElem, rowSize, 00889 array, nullValue,&status)) throw FitsError(status); 00890 } 00891 else 00892 { 00893 if (fits_write_col(fitsPointer(),abs(type()),index(),row,firstElem,rowSize, 00894 array,&status)) throw FitsError(status); 00895 00896 } 00897 } 00898 00899 // Additional Declarations 00900 00901 // all functions that operate on complex data that call cfitsio 00902 // need to be specialized. The signature containing complex<T>* objects 00903 // is unfortunate, perhaps, for this purpose, but the user will access 00904 // rw operations through standard library containers. 00905 00906 00907 00908 00909 00910 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT 00911 template <> 00912 inline void ColumnVectorData<complex<float> >::setDataLimits (complex<float>* limits) 00913 { 00914 m_minLegalValue = limits[0]; 00915 m_maxLegalValue = limits[1]; 00916 m_minDataValue = limits[2]; 00917 m_maxDataValue = limits[3]; 00918 } 00919 #else 00920 template <> 00921 void 00922 ColumnVectorData<complex<float> >::setDataLimits (complex<float>* limits); 00923 #endif 00924 00925 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT 00926 template <> 00927 inline void ColumnVectorData<complex<double> >::setDataLimits (complex<double>* limits) 00928 { 00929 m_minLegalValue = limits[0]; 00930 m_maxLegalValue = limits[1]; 00931 m_minDataValue = limits[2]; 00932 m_maxDataValue = limits[3]; 00933 } 00934 #else 00935 template <> 00936 void 00937 ColumnVectorData<complex<double> >::setDataLimits (complex<double>* limits); 00938 #endif 00939 00940 00941 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT 00942 template <> 00943 inline void ColumnVectorData<std::complex<float> >::readColumnData(long firstRow, 00944 long nelements, long firstElem, std::complex<float>* null ) 00945 { 00946 int status=0; 00947 float nulval (0); 00948 FITSUtil::auto_array_ptr<float> pArray(new float[2*nelements]); 00949 float* array = pArray.get(); 00950 int anynul(0); 00951 00952 if (fits_read_col_cmp(fitsPointer(),index(),firstRow, firstElem, 00953 nelements,nulval,array,&anynul,&status) ) throw FitsError(status); 00954 00955 if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows()); 00956 00957 std::valarray<std::complex<float> > readData(nelements); 00958 for (long j = 0; j < nelements; ++j) 00959 { 00960 readData[j] = std::complex<float>(array[2*j],array[2*j+1]); 00961 } 00962 size_t countRead = 0; 00963 const size_t ONE = 1; 00964 00965 if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows()); 00966 size_t vectorSize(0); 00967 if (!varLength()) 00968 { 00969 vectorSize = std::max(repeat(),ONE); // safety check. 00970 } 00971 else 00972 { 00973 // assume that the user specified the correct length for 00974 // variable columns. This should be ok since readVariableColumns 00975 // uses fits_read_descripts to return this information from the 00976 // fits pointer, and this is passed as nelements here. 00977 vectorSize = nelements; 00978 } 00979 size_t n = nelements; 00980 00981 int i = firstRow; 00982 int ii = i - 1; 00983 while ( countRead < n) 00984 { 00985 std::valarray<complex<float> >& current = m_data[ii]; 00986 if (current.size() != vectorSize) current.resize(vectorSize,0.); 00987 int elementsInFirstRow = vectorSize-firstElem + 1; 00988 bool lastRow = ( (nelements - countRead) < vectorSize); 00989 if (lastRow) 00990 { 00991 int elementsInLastRow = nelements - countRead; 00992 std::copy(&readData[countRead],&readData[0]+nelements,¤t[0]); 00993 countRead += elementsInLastRow; 00994 } 00995 // what to do with complete rows. if firstElem == 1 the 00996 else 00997 { 00998 if (firstElem == 1 || (firstElem > 1 && i > firstRow) ) 00999 { 01000 current = readData[std::slice(vectorSize*(ii-firstRow)+ 01001 elementsInFirstRow,vectorSize,1)]; 01002 ++ii; 01003 ++i; 01004 countRead += vectorSize; 01005 } 01006 else 01007 { 01008 if (i == firstRow) 01009 { 01010 std::copy(&readData[0],&readData[0]+elementsInFirstRow, 01011 ¤t[firstElem]); 01012 countRead += elementsInFirstRow; 01013 ++i; 01014 ++ii; 01015 } 01016 } 01017 } 01018 } 01019 } 01020 #else 01021 template <> 01022 void ColumnVectorData<complex<float> >::readColumnData(long firstRow, 01023 long nelements, 01024 long firstElem, complex<float>* null); 01025 #endif 01026 01027 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT 01028 template <> 01029 inline void ColumnVectorData<complex<double> >::readColumnData (long firstRow, 01030 long nelements,long firstElem, 01031 complex<double>* nullValue) 01032 { 01033 01034 // duplicated for each complex type to work around imagined or 01035 // actual compiler deficiencies. 01036 int status=0; 01037 double nulval (0); 01038 FITSUtil::auto_array_ptr<double> pArray(new double[2*nelements]); 01039 double* array = pArray.get(); 01040 int anynul(0); 01041 01042 if (fits_read_col_dblcmp(fitsPointer(),index(),firstRow, firstElem, 01043 nelements,nulval,array,&anynul,&status) ) throw FitsError(status); 01044 01045 if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows()); 01046 01047 std::valarray<std::complex<double> > readData(nelements); 01048 for (long j = 0; j < nelements; ++j) 01049 { 01050 readData[j] = std::complex<double>(array[2*j],array[2*j+1]); 01051 } 01052 size_t countRead = 0; 01053 const size_t ONE = 1; 01054 01055 if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows()); 01056 size_t vectorSize(0); 01057 if (!varLength()) 01058 { 01059 vectorSize = std::max(repeat(),ONE); // safety check. 01060 } 01061 else 01062 { 01063 // assume that the user specified the correct length for 01064 // variable columns. This should be ok since readVariableColumns 01065 // uses fits_read_descripts to return this information from the 01066 // fits pointer, and this is passed as nelements here. 01067 vectorSize = nelements; 01068 } 01069 size_t n = nelements; 01070 01071 int i = firstRow; 01072 int ii = i - 1; 01073 while ( countRead < n) 01074 { 01075 std::valarray<std::complex<double> >& current = m_data[ii]; 01076 if (current.size() != vectorSize) current.resize(vectorSize,0.); 01077 int elementsInFirstRow = vectorSize-firstElem + 1; 01078 bool lastRow = ( (nelements - countRead) < vectorSize); 01079 if (lastRow) 01080 { 01081 int elementsInLastRow = nelements - countRead; 01082 std::copy(&readData[countRead],&readData[0]+nelements,¤t[0]); 01083 countRead += elementsInLastRow; 01084 } 01085 // what to do with complete rows. if firstElem == 1 the 01086 else 01087 { 01088 if (firstElem == 1 || (firstElem > 1 && i > firstRow) ) 01089 { 01090 current = readData[std::slice(vectorSize*(ii-firstRow)+ 01091 elementsInFirstRow,vectorSize,1)]; 01092 ++ii; 01093 ++i; 01094 countRead += vectorSize; 01095 } 01096 else 01097 { 01098 if (i == firstRow) 01099 { 01100 std::copy(&readData[0],&readData[0]+elementsInFirstRow, 01101 ¤t[firstElem]); 01102 countRead += elementsInFirstRow; 01103 ++i; 01104 ++ii; 01105 } 01106 } 01107 } 01108 } 01109 } 01110 #else 01111 template <> 01112 void ColumnVectorData<complex<double> >::readColumnData (long firstRow, 01113 long nelements, 01114 long firstElem, complex<double>* null); 01115 #endif 01116 01117 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT 01118 template <> 01119 inline void ColumnVectorData<complex<float> >::writeFixedArray 01120 (complex<float>* data, long nElements, long nRows, long firstRow, 01121 complex<float>* nullValue) 01122 { 01123 01124 int status(0); 01125 01126 // check for sanity of inputs, then write to file. 01127 // this function writes only complete rows to a table with 01128 // fixed width rows. 01129 01130 01131 if ( nElements < nRows*static_cast<long>(repeat()) ) 01132 { 01133 #ifdef SSTREAM_DEFECT 01134 std::ostrstream msgStr; 01135 #else 01136 std::ostringstream msgStr; 01137 #endif 01138 msgStr << " input array size: " << nElements 01139 << " required " << nRows*repeat(); 01140 #ifdef SSTREAM_DEFECT 01141 msgStr << std::ends; 01142 #endif 01143 01144 01145 String msg(msgStr.str()); 01146 01147 throw Column::InsufficientElements(msg); 01148 } 01149 01150 FITSUtil::auto_array_ptr<float> realData(new float[2*nElements]); 01151 01152 for (int j = 0; j < nElements; ++j) 01153 { 01154 realData[2*j] = data[j].real(); 01155 realData[2*j+1] = data[j].imag(); 01156 } 01157 01158 01159 01160 if (fits_write_col_cmp(fitsPointer(),index(),firstRow, 01161 1,nElements,realData.get(),&status)) throw FitsError(status); 01162 01163 parent()->updateRows(); 01164 } 01165 #else 01166 template <> 01167 void ColumnVectorData<complex<float> >::writeFixedArray 01168 (complex<float>* data, long nElements, long nRows, long firstRow, std::complex<float>* null); 01169 #endif 01170 01171 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT 01172 template <> 01173 inline void ColumnVectorData<complex<double> >::writeFixedArray 01174 (complex<double>* data, long nElements, long nRows, long firstRow, 01175 complex<double>* nullValue) 01176 { 01177 int status(0); 01178 01179 // check for sanity of inputs, then write to file. 01180 // this function writes only complete rows to a table with 01181 // fixed width rows. 01182 01183 01184 if ( nElements < nRows*static_cast<long>(repeat()) ) 01185 { 01186 #ifdef SSTREAM_DEFECT 01187 std::ostrstream msgStr; 01188 #else 01189 std::ostringstream msgStr; 01190 #endif 01191 msgStr << " input array size: " << nElements 01192 << " required " << nRows*repeat(); 01193 #ifdef SSTREAM_DEFECT 01194 msgStr << std::ends; 01195 #endif 01196 01197 String msg(msgStr.str()); 01198 01199 throw Column::InsufficientElements(msg); 01200 } 01201 01202 FITSUtil::auto_array_ptr<double> realData(new double[2*nElements]); 01203 01204 for (int j = 0; j < nElements; ++j) 01205 { 01206 realData[2*j] = data[j].real(); 01207 realData[2*j+1] = data[j].imag(); 01208 } 01209 01210 01211 01212 if (fits_write_col_dblcmp(fitsPointer(),index(),firstRow, 01213 1,nElements,realData.get(),&status)) throw FitsError(status); 01214 01215 parent()->updateRows(); 01216 01217 } 01218 #else 01219 template <> 01220 void ColumnVectorData<complex<double> >::writeFixedArray 01221 (complex<double>* data, long nElements, long nRows, long firstRow, 01222 std::complex<double>* null); 01223 #endif 01224 01225 #ifdef SPEC_TEMPLATE_DECL_DEFECT 01226 template <> 01227 inline void 01228 ColumnVectorData<std::complex<float> >::doWrite 01229 (std::complex<float>* data, long row, long rowSize, long firstElem, std::complex<float>* nullValue ) 01230 { 01231 int status(0); 01232 FITSUtil::auto_array_ptr<float> carray( new float[2*rowSize]); 01233 for ( long j = 0 ; j < rowSize; ++ j) 01234 { 01235 carray[2*j] = data[j].real(); 01236 carray[2*j + 1] = data[j].imag(); 01237 } 01238 if (fits_write_col_cmp(fitsPointer(),index(),row,firstElem,rowSize, 01239 carray.get(),&status)) throw FitsError(status); 01240 } 01241 01242 01243 template <> 01244 inline void 01245 ColumnVectorData<std::complex<double> >::doWrite 01246 (std::complex<double>* data, long row, long rowSize, long firstElem, std::complex<double>* nullValue ) 01247 { 01248 int status(0); 01249 FITSUtil::auto_array_ptr<double> carray( new double[2*rowSize]); 01250 for ( long j = 0 ; j < rowSize; ++ j) 01251 { 01252 carray[2*j] = data[j].real(); 01253 carray[2*j + 1] = data[j].imag(); 01254 } 01255 if (fits_write_col_dblcmp(fitsPointer(),index(),row,firstElem,rowSize, 01256 carray.get(),&status)) throw FitsError(status); 01257 01258 } 01259 01260 #else 01261 template<> 01262 void 01263 ColumnVectorData<complex<float> >::doWrite 01264 ( complex<float>* data, long row, long rowSize, long firstElem, complex<float>* nullValue); 01265 01266 template<> 01267 void 01268 ColumnVectorData<complex<double> >::doWrite 01269 ( complex<double>* data, long row, long rowSize, long firstElem, complex<double>* nullValue ); 01270 #endif 01271 } // namespace CCfits 01272 01273 01274 #endif