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