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