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