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