ColumnVectorData.h

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

Generated on Fri Nov 3 17:09:05 2006 for CCfits by  doxygen 1.4.7