ColumnVectorData.h

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

Generated on Thu Jun 28 11:49:08 2007 for CCfits by  doxygen 1.4.7