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