CCfits  2.5
ColumnVectorData.h
1 // Astrophysics Science Division,
2 // NASA/ Goddard Space Flight Center
3 // HEASARC
4 // http://heasarc.gsfc.nasa.gov
5 // e-mail: ccfits@legacy.gsfc.nasa.gov
6 //
7 // Original author: Ben Dorman
8 
9 #ifndef COLUMNVECTORDATA_H
10 #define COLUMNVECTORDATA_H 1
11 #ifdef _MSC_VER
12 #include "MSconfig.h"
13 #endif
14 #include "CCfits.h"
15 
16 // valarray
17 #include <valarray>
18 // vector
19 #include <vector>
20 // Column
21 #include "Column.h"
22 
23 #ifdef SSTREAM_DEFECT
24 #include <strstream>
25 #else
26 #include <sstream>
27 #endif
28 
29 #include <memory>
30 #include <numeric>
31 #include <algorithm>
32 
33 namespace CCfits {
34 
35  class Table;
36 
37 }
38 
39 #include "FITS.h"
40 #include "FITSUtil.h"
41 using std::complex;
42 
43 
44 namespace CCfits {
45 
46 
47 
48  template <typename T>
49  class ColumnVectorData : public Column //## Inherits: <unnamed>%38BAD1D4D370
50  {
51 
52  public:
53  ColumnVectorData(const ColumnVectorData< T > &right);
54  ColumnVectorData (Table* p = 0);
55  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 = "");
56  ~ColumnVectorData();
57 
58  virtual void readData (long firstrow, long nelements, long firstelem = 1);
59  virtual ColumnVectorData<T>* clone () const;
60  virtual void setDimen ();
61  void setDataLimits (T* limits);
62  const T minLegalValue () const;
63  void minLegalValue (T value);
64  const T maxLegalValue () const;
65  void maxLegalValue (T value);
66  const T minDataValue () const;
67  void minDataValue (T value);
68  const T maxDataValue () const;
69  void maxDataValue (T value);
70  const std::vector<std::valarray<T> >& data () const;
71  void setData (const std::vector<std::valarray<T> >& value);
72  const std::valarray<T>& data (int i) const;
73  void data (int i, const std::valarray<T>& value);
74 
75  // Additional Public Declarations
76  friend class Column;
77  protected:
78  // Additional Protected Declarations
79 
80  private:
81  ColumnVectorData< T > & operator=(const ColumnVectorData< T > &right);
82 
83  virtual bool compare (const Column &right) const;
84  void resizeDataObject (const std::vector<std::valarray<T> >& indata, size_t firstRow);
85  // Reads a specified number of column rows.
86  //
87  // There are no default arguments. The function
88  // Column::read(firstrow,firstelem,nelements)
89  // is designed for reading the whole column.
90  virtual void readColumnData (long first, long last, T* nullValue = 0);
91  virtual std::ostream& put (std::ostream& s) const;
92  void writeData (const std::valarray<T>& indata, long numRows, long firstRow = 1, T* nullValue = 0);
93  void writeData (const std::vector<std::valarray<T> >& indata, long firstRow = 1, T* nullValue = 0);
94  // Reads a specified number of column rows.
95  //
96  // There are no default arguments. The function
97  // Column::read(firstrow,firstelem,nelements)
98  // is designed for reading the whole column.
99  virtual void readRow (size_t row, T* nullValue = 0);
100  // Reads a variable row..
101  virtual void readVariableRow (size_t row, T* nullValue = 0);
102  void readColumnData (long firstrow, long nelements, long firstelem, T* nullValue = 0);
103  void writeData (const std::valarray<T>& indata, const std::vector<long>& vectorLengths, long firstRow = 1, T* nullValue = 0);
104  void writeFixedRow (const std::valarray<T>& data, long row, long firstElem = 1, T* nullValue = 0);
105  void writeFixedArray (T* data, long nElements, long nRows, long firstRow, T* nullValue = 0);
106  // Insert one or more blank rows into a FITS column.
107  virtual void insertRows (long first, long number = 1);
108  virtual void deleteRows (long first, long number = 1);
109  void doWrite (T* array, long row, long rowSize, long firstElem, T* nullValue);
110 
111  // Additional Private Declarations
112 
113  private: //## implementation
114  // Data Members for Class Attributes
115  T m_minLegalValue;
116  T m_maxLegalValue;
117  T m_minDataValue;
118  T m_maxDataValue;
119 
120  // Data Members for Associations
121  std::vector<std::valarray<T> > m_data;
122 
123  // Additional Implementation Declarations
124 
125  };
126 
127  // Parameterized Class CCfits::ColumnVectorData
128 
129  template <typename T>
130  inline void ColumnVectorData<T>::readData (long firstrow, long nelements, long firstelem)
131  {
132  readColumnData(firstrow,nelements,firstelem,static_cast<T*>(0));
133  }
134 
135  template <typename T>
136  inline const T ColumnVectorData<T>::minLegalValue () const
137  {
138  return m_minLegalValue;
139  }
140 
141  template <typename T>
142  inline void ColumnVectorData<T>::minLegalValue (T value)
143  {
144  m_minLegalValue = value;
145  }
146 
147  template <typename T>
148  inline const T ColumnVectorData<T>::maxLegalValue () const
149  {
150  return m_maxLegalValue;
151  }
152 
153  template <typename T>
154  inline void ColumnVectorData<T>::maxLegalValue (T value)
155  {
156  m_maxLegalValue = value;
157  }
158 
159  template <typename T>
160  inline const T ColumnVectorData<T>::minDataValue () const
161  {
162  return m_minDataValue;
163  }
164 
165  template <typename T>
166  inline void ColumnVectorData<T>::minDataValue (T value)
167  {
168  m_minDataValue = value;
169  }
170 
171  template <typename T>
172  inline const T ColumnVectorData<T>::maxDataValue () const
173  {
174  return m_maxDataValue;
175  }
176 
177  template <typename T>
178  inline void ColumnVectorData<T>::maxDataValue (T value)
179  {
180  m_maxDataValue = value;
181  }
182 
183  template <typename T>
184  inline const std::vector<std::valarray<T> >& ColumnVectorData<T>::data () const
185  {
186  return m_data;
187  }
188 
189  template <typename T>
190  inline void ColumnVectorData<T>::setData (const std::vector<std::valarray<T> >& value)
191  {
192  m_data = value;
193  }
194 
195  template <typename T>
196  inline const std::valarray<T>& ColumnVectorData<T>::data (int i) const
197  {
198  return m_data[i - 1];
199  }
200 
201  template <typename T>
202  inline void ColumnVectorData<T>::data (int i, const std::valarray<T>& value)
203  {
204  if (m_data[i-1].size() != value.size())
205  m_data[i-1].resize(value.size());
206  m_data[i - 1] = value;
207  }
208 
209  // Parameterized Class CCfits::ColumnVectorData
210 
211  template <typename T>
212  ColumnVectorData<T>::ColumnVectorData(const ColumnVectorData<T> &right)
213  :Column(right),
214  m_minLegalValue(right.m_minLegalValue),
215  m_maxLegalValue(right.m_maxLegalValue),
216  m_minDataValue(right.m_minDataValue),
217  m_maxDataValue(right.m_maxDataValue),
218  m_data(right.m_data)
219  {
220  }
221 
222  template <typename T>
223  ColumnVectorData<T>::ColumnVectorData (Table* p)
224  : Column(p),
225  m_minLegalValue(0),
226  m_maxLegalValue(0),
227  m_minDataValue(0),
228  m_maxDataValue(0),
229  m_data()
230  {
231  }
232 
233  template <typename T>
234  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)
235  : Column(columnIndex,columnName,type,format,unit,p,rpt,w,comment),
236  m_minLegalValue(0),
237  m_maxLegalValue(0),
238  m_minDataValue(0),
239  m_maxDataValue(0),
240  m_data()
241  {
242  }
243 
244 
245  template <typename T>
246  ColumnVectorData<T>::~ColumnVectorData()
247  {
248  // valarray destructor should do all the work.
249  }
250 
251 
252  template <typename T>
253  bool ColumnVectorData<T>::compare (const Column &right) const
254  {
255  if ( !Column::compare(right) ) return false;
256  const ColumnVectorData<T>& that = static_cast<const ColumnVectorData<T>&>(right);
257  size_t n = m_data.size();
258  // m_data is of type valarray<T>.
259  if ( that.m_data.size() != n ) return false;
260  for (size_t i = 0; i < n ; i++)
261  {
262  const std::valarray<T>& thisValArray=m_data[i];
263  const std::valarray<T>& thatValArray=that.m_data[i];
264  size_t nn = thisValArray.size();
265  if (thatValArray.size() != nn ) return false;
266 
267  for (size_t j = 0; j < nn ; j++ )
268  {
269  if (thisValArray[j] != thatValArray[j])
270  return false;
271  }
272  }
273  return true;
274  }
275 
276  template <typename T>
277  ColumnVectorData<T>* ColumnVectorData<T>::clone () const
278  {
279  return new ColumnVectorData<T>(*this);
280  }
281 
282  template <typename T>
283  void ColumnVectorData<T>::resizeDataObject (const std::vector<std::valarray<T> >& indata, size_t firstRow)
284  {
285  // the rows() call is the value before updating.
286  // the updateRows() call at the end sets the call to return the
287  // value from the fits pointer - which is changed by writeFixedArray
288  // or writeFixedRow.
289 
290  const size_t lastInputRow(indata.size() + firstRow - 1);
291  const size_t newLastRow = std::max(lastInputRow,static_cast<size_t>(rows()));
292 
293  // if the write instruction increases the rows, we need to add
294  // rows to the data member and preserve its current contents.
295 
296  // rows() >= origNRows since it is the value for entire table,
297  // not just this column.
298  const size_t origNRows(m_data.size());
299  // This will always be an expansion. vector.resize() doesn't
300  // invalidate any data on an expansion.
301  if (newLastRow > origNRows) m_data.resize(newLastRow);
302 
303  if (varLength())
304  {
305  // The incoming data will determine each row size, thus
306  // no need to preserve any existing values in the row.
307  // Each value will eventually be overwritten.
308  for (size_t iRow = firstRow-1; iRow < lastInputRow; ++iRow)
309  {
310  std::valarray<T>& current = m_data[iRow];
311  const size_t newSize = indata[iRow - (firstRow-1)].size();
312  if (current.size() != newSize)
313  current.resize(newSize);
314  }
315  }
316  else
317  {
318  // All row sizes in m_data should ALWAYS be either repeat(),
319  // or 0 if they haven't been initialized. This is true regardless
320  // of the incoming data row size.
321 
322  // Perform LAZY initialization of m_data rows. Only
323  // expand a row valarray when it is first needed.
324  for (size_t iRow = firstRow-1; iRow < lastInputRow; ++iRow)
325  {
326  if (m_data[iRow].size() != repeat())
327  m_data[iRow].resize(repeat());
328  }
329  }
330  }
331 
332  template <typename T>
333  void ColumnVectorData<T>::setDimen ()
334  {
335  int status(0);
336  FITSUtil:: auto_array_ptr<char> dimValue (new char[FLEN_VALUE]);
337 
338 #ifdef SSTREAM_DEFECT
339  std::ostrstream key;
340 #else
341  std::ostringstream key;
342 #endif
343  key << "TDIM" << index();
344 
345 #ifdef SSTREAM_DEFECT
346  fits_read_key_str(fitsPointer(), key.str(), dimValue.get(),0,&status);
347 #else
348  fits_read_key_str(fitsPointer(),const_cast<char*>(key.str().c_str()),dimValue.get(),0,&status);
349 #endif
350 
351  if (status == 0)
352  {
353  dimen(String(dimValue.get()));
354  }
355  }
356 
357  template <typename T>
358  void ColumnVectorData<T>::readColumnData (long first, long last, T* nullValue)
359  {
360  makeHDUCurrent();
361 
362 
363  if ( rows() < last )
364  {
365  std::cerr << "CCfits: More data requested than contained in table. ";
366  std::cerr << "Extracting complete column.\n";
367  last = rows();
368  }
369 
370  long nelements = (last - first + 1)*repeat();
371 
372 
373  readColumnData(first,nelements,1,nullValue);
374  if (first <= 1 && last == rows()) isRead(true);
375  }
376 
377  template <typename T>
378  std::ostream& ColumnVectorData<T>::put (std::ostream& s) const
379  {
380  // output header information
381  Column::put(s);
382  if ( FITS::verboseMode() )
383  {
384  s << " Column Legal limits: ( " << m_minLegalValue << "," << m_maxLegalValue << " )\n"
385  << " Column Data limits: ( " << m_minDataValue << "," << m_maxDataValue << " )\n";
386  }
387  if (!m_data.empty())
388  {
389  for (size_t j = 0; j < m_data.size(); j++)
390  {
391  size_t n = m_data[j].size();
392  if ( n )
393  {
394  s << "Row " << j + 1 << " Vector Size " << n << '\n';
395  for (size_t k = 0; k < n - 1; k++)
396  {
397  s << m_data[j][k] << '\t';
398  }
399  s << m_data[j][n - 1] << '\n';
400  }
401  }
402  }
403 
404  return s;
405  }
406 
407  template <typename T>
408  void ColumnVectorData<T>::writeData (const std::valarray<T>& indata, long numRows, long firstRow, T* nullValue)
409  {
410  // This version of writeData is called by Column write functions which
411  // can only write the same number of elements to each row.
412  // For fixed width columns, this must be equal to the repeat value
413  // or an exception is thrown. For variable width, it only requires
414  // that indata.size()/numRows is an int.
415 
416  // won't do anything if < 0, and will give divide check if == 0.
417  if (numRows <= 0) throw InvalidNumberOfRows(numRows);
418 
419 #ifdef SSTREAM_DEFECT
420  std::ostrstream msgStr;
421 #else
422  std::ostringstream msgStr;
423 #endif
424  if (indata.size() % static_cast<size_t>(numRows))
425  {
426  msgStr << "To use this write function, input array size"
427  <<"\n must be exactly divisible by requested num rows: "
428  << numRows;
429  throw InsufficientElements(msgStr.str());
430  }
431  const size_t cellsize = indata.size()/static_cast<size_t>(numRows);
432 
433  if (!varLength() && cellsize != repeat() )
434  {
435  msgStr << "column: " << name()
436  << "\n input data size: " << indata.size()
437  << " required: " << numRows*repeat();
438  String msg(msgStr.str());
439  throw InsufficientElements(msg);
440  }
441 
442  std::vector<std::valarray<T> > internalFormat(numRows);
443 
444  // support writing equal row lengths to variable columns.
445 
446  for (long j = 0; j < numRows; ++j)
447  {
448  internalFormat[j].resize(cellsize);
449  internalFormat[j] = indata[std::slice(cellsize*j,cellsize,1)];
450  }
451 
452  // change the size of m_data based on the first row to be written
453  // and on the input data vector sizes.
454 
455  writeData(internalFormat,firstRow,nullValue);
456  }
457 
458  template <typename T>
459  void ColumnVectorData<T>::writeData (const std::vector<std::valarray<T> >& indata, long firstRow, T* nullValue)
460  {
461  // This is called directly by Column's writeArrays functions, and indirectly
462  // by both categories of write functions, ie. those which allow differing
463  // lengths per row and those that don't.
464  const size_t nInputRows(indata.size());
465  using std::valarray;
466 
467  resizeDataObject(indata,firstRow);
468  // After the above call, can assume all m_data arrays to be written to
469  // have been properly resized whether we're dealing with fixed or
470  // variable length.
471 
472  if (varLength())
473  {
474  // firstRow is 1-based, but all these internal row variables
475  // will be 0-based.
476  const size_t endRow = nInputRows + firstRow-1;
477  for (size_t iRow = firstRow-1; iRow < endRow; ++iRow)
478  {
479  m_data[iRow] = indata[iRow - (firstRow-1)];
480  // doWrite wants 1-based rows.
481  doWrite(&m_data[iRow][0], iRow+1, m_data[iRow].size(), 1, nullValue);
482  }
483  parent()->updateRows();
484  }
485  else
486  {
487  // Check for simplest case of all valarrays of size repeat().
488  // If any are greater, throw an error.
489  const size_t colRepeat = repeat();
490  bool allEqualRepeat = true;
491  for (size_t i=0; i<nInputRows; ++i)
492  {
493  const size_t sz = indata[i].size();
494  if (sz > colRepeat)
495  {
496 #ifdef SSTREAM_DEFECT
497  std::ostrstream oss;
498 #else
499  std::ostringstream oss;
500 #endif
501  oss << " vector column length " << colRepeat
502  <<", input valarray length " << sz;
503  throw InvalidRowParameter(oss.str());
504  }
505  if (sz < colRepeat)
506  allEqualRepeat = false;
507  }
508 
509  if (allEqualRepeat)
510  {
511  // concatenate the valarrays and write.
512  const size_t nElements (colRepeat*nInputRows);
513  FITSUtil::CVAarray<T> convert;
514  FITSUtil::auto_array_ptr<T> pArray(convert(indata));
515  T* array = pArray.get();
516 
517  // if T is complex, then CVAarray returns a
518  // C-array of complex objects. But FITS requires an array of complex's
519  // value_type.
520 
521  // This writes to the file and also calls updateRows.
522  writeFixedArray(array,nElements,nInputRows,firstRow,nullValue);
523 
524  for (size_t j = 0; j < nInputRows ; ++j)
525  {
526  const valarray<T>& input = indata[j];
527  valarray<T>& current = m_data[j + firstRow - 1];
528  // current should be resized by resizeDataObject.
529  current = input;
530  }
531  }
532  else
533  {
534  // Some input arrays have fewer than colRepeat elements.
535  const size_t endRow = nInputRows + firstRow-1;
536  for (size_t iRow = firstRow-1; iRow<endRow; ++iRow)
537  {
538  // resizeDataObject should already have resized all
539  // corresponding m_data rows to repeat().
540  const valarray<T>& input = indata[iRow-(firstRow-1)];
541  writeFixedRow(input, iRow, 1, nullValue);
542  }
543  parent()->updateRows();
544  }
545 
546  } // end if !varLength
547  }
548 
549  template <typename T>
550  void ColumnVectorData<T>::readRow (size_t row, T* nullValue)
551  {
552  makeHDUCurrent();
553 
554 
555 
556  if ( row > static_cast<size_t>(rows()) )
557  {
558 #ifdef SSTREAM_DEFECT
559  std::ostrstream msg;
560 #else
561  std::ostringstream msg;
562 #endif
563  msg << " row requested: " << row << " row range: 1 - " << rows();
564 #ifdef SSTREAM_DEFECT
565  msg << std::ends;
566 #endif
567 
568  throw Column::InvalidRowNumber(msg.str());
569  }
570 
571  // this is really for documentation purposes. I expect the optimizer will
572  // remove this redundant definition .
573  bool variable(type() < 0);
574 
575 
576  long nelements(repeat());
577 
578  if (variable)
579  {
580  readVariableRow(row,nullValue);
581  }
582  else
583  {
584  readColumnData(row,nelements,1,nullValue);
585  }
586  }
587 
588  template <typename T>
589  void ColumnVectorData<T>::readVariableRow (size_t row, T* nullValue)
590  {
591  int status(0);
592  long offset(0);
593  long repeat(0);
594  if (fits_read_descript(fitsPointer(),index(),static_cast<long>(row),
595  &repeat,&offset,&status)) throw FitsError(status);
596  readColumnData(row,repeat,1,nullValue);
597  }
598 
599  template <typename T>
600  void ColumnVectorData<T>::readColumnData (long firstrow, long nelements, long firstelem, T* nullValue)
601  {
602  int status=0;
603 
604  FITSUtil::auto_array_ptr<T> pArray(new T[nelements]);
605  T* array = pArray.get();
606  int anynul(0);
607 
608 
609 
610  if (fits_read_col(fitsPointer(), abs(type()),index(), firstrow, firstelem,
611  nelements, nullValue, array, &anynul, &status) != 0)
612  throw FitsError(status);
613 
614  size_t countRead = 0;
615  const size_t ONE = 1;
616 
617  if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows());
618  size_t vectorSize(0);
619  if (!varLength())
620  {
621 
622  vectorSize = std::max(repeat(),ONE); // safety check.
623 
624  }
625  else
626  {
627  // assume that the user specified the correct length for
628  // variable columns. This should be ok since readVariableColumns
629  // uses fits_read_descripts to return this information from the
630  // fits pointer, and this is passed as nelements here.
631  vectorSize = nelements;
632  }
633  size_t n = nelements;
634 
635  int i = firstrow;
636  int ii = i - 1;
637  while ( countRead < n)
638  {
639  std::valarray<T>& current = m_data[ii];
640  if (current.size() != vectorSize) current.resize(vectorSize);
641  int elementsInFirstRow = vectorSize-firstelem + 1;
642  bool lastRow = ( (nelements - countRead) < vectorSize);
643  if (lastRow)
644  {
645  int elementsInLastRow = nelements - countRead;
646  std::valarray<T> ttmp(array + vectorSize*(ii-firstrow) + elementsInFirstRow,
647  elementsInLastRow);
648  for (int kk = 0; kk < elementsInLastRow; kk++) current[kk] = ttmp[kk];
649  countRead += elementsInLastRow;
650 
651  }
652  // what to do with complete rows
653  else
654  {
655  if (firstelem == 1 || (firstelem > 1 && i > firstrow) )
656  {
657  std::valarray<T> ttmp(array + vectorSize*(ii - firstrow) +
658  elementsInFirstRow,vectorSize);
659  current = ttmp;
660  ii++;
661  i++;
662  countRead += vectorSize;
663  }
664  else
665  {
666  if (i == firstrow)
667  {
668  std::valarray<T> ttmp(array,elementsInFirstRow);
669  for (size_t kk = firstelem ; kk < vectorSize ; kk++)
670  current[kk] = ttmp[kk-firstelem];
671  countRead += elementsInFirstRow;
672  i++;
673  ii++;
674  }
675  }
676  }
677  }
678  }
679 
680  template <typename T>
681  void ColumnVectorData<T>::writeData (const std::valarray<T>& indata, const std::vector<long>& vectorLengths, long firstRow, T* nullValue)
682  {
683  // Called from Column write functions which allow differing lengths
684  // for each row.
685  using namespace std;
686  const size_t N(vectorLengths.size());
687  vector<long> sums(N);
688  // pre-calculate partial sums of vector lengths for use as array offsets.
689  partial_sum(vectorLengths.begin(),vectorLengths.end(),sums.begin());
690  // check that sufficient data have been supplied to carry out the entire operation.
691  if (indata.size() < static_cast<size_t>(sums[N-1]) )
692  {
693 #ifdef SSTREAM_DEFECT
694  ostrstream msgStr;
695 #else
696  ostringstream msgStr;
697 #endif
698  msgStr << " input data size: " << indata.size() << " vector length sum: " << sums[N-1];
699 #ifdef SSTREAM_DEFECT
700  msgStr << std::ends;
701 #endif
702 
703  String msg(msgStr.str());
704  throw InsufficientElements(msg);
705  }
706 
707  vector<valarray<T> > vvArray(N);
708  long& last = sums[0];
709  vvArray[0].resize(last);
710  for (long jj = 0; jj < last; ++jj) vvArray[0][jj] = indata[jj];
711 
712  for (size_t j = 1; j < N; ++j)
713  {
714  valarray<T>& __tmp = vvArray[j];
715  // these make the code much more readable
716  long& first = sums[j-1];
717  long& jlast = sums[j];
718  __tmp.resize(jlast - first);
719  for (long k = first; k < jlast; ++k)
720  {
721  __tmp[k - first] = indata[k];
722  }
723  }
724 
725  writeData(vvArray,firstRow,nullValue);
726  }
727 
728  template <typename T>
729  void ColumnVectorData<T>::writeFixedRow (const std::valarray<T>& data, long row, long firstElem, T* nullValue)
730  {
731 
732  // This is to be called only for FIXED length vector columns. It will
733  // throw if data.size()+firstElem goes beyond the repeat value.
734  // If data.size() is less than repeat, it leaves the remaining values
735  // undisturbed both in the file and in m_data storage.
736 
737 #ifdef SSTREAM_DEFECT
738  std::ostrstream msgStr;
739 #else
740  std::ostringstream msgStr;
741 #endif
742  if (varLength())
743  {
744  msgStr <<"Calling ColumnVectorData::writeFixedRow for a variable length column.\n";
745  throw FitsFatal(msgStr.str());
746  }
747 
748  std::valarray<T>& storedRow = m_data[row];
749  long inputSize = static_cast<long>(data.size());
750  long storedSize(storedRow.size());
751  if (storedSize != static_cast<long>(repeat()))
752  {
753  msgStr<<"stored array size vs. column width mismatch in ColumnVectorData::writeFixedRow.\n";
754  throw FitsFatal(msgStr.str());
755  }
756 
757  if (inputSize + firstElem - 1 > storedSize)
758  {
759  msgStr << " requested write " << firstElem << " to "
760  << firstElem + inputSize - 1 << " exceeds vector length " << repeat();
761  throw InvalidRowParameter(msgStr.str());
762  }
763 
764  // CANNOT give a strong exception safety guarantee because writing
765  // data changes the file. Any corrective action that could be taken
766  // [e.g. holding initial contents of the row and writing it back after
767  // an exception is thrown] could in principle throw the same exception
768  // we are trying to protect from.
769 
770  // routine does however give the weak guarantee (no resource leaks).
771 
772  // It's never a good thing to cast away a const, but doWrite calls the
773  // CFITSIO write functions which take a non-const pointer (though
774  // it shouldn't actually modify the array), and I'd rather not
775  // copy the entire valarray just to avoid this problem.
776  std::valarray<T>& lvData = const_cast<std::valarray<T>&>(data);
777  T* inPointer = &lvData[0];
778  doWrite(inPointer, row+1, inputSize, firstElem, nullValue);
779 
780  // Writing to disk was successful, now update FITS object and return.
781  const size_t offset = static_cast<size_t>(firstElem) - 1;
782  for (size_t iElem=0; iElem < static_cast<size_t>(inputSize); ++iElem)
783  {
784  // This doesn't require inPointer's non-constness. It's just
785  // used here to speed things up a bit.
786  storedRow[iElem + offset] = inPointer[iElem];
787  }
788  }
789 
790  template <typename T>
791  void ColumnVectorData<T>::writeFixedArray (T* data, long nElements, long nRows, long firstRow, T* nullValue)
792  {
793  int status(0);
794 
795  // check for sanity of inputs, then write to file.
796  // this function writes only complete rows to a table with
797  // fixed width rows.
798 
799 
800  if ( nElements < nRows*static_cast<long>(repeat()) )
801  {
802 #ifdef SSTREAM_DEFECT
803  std::ostrstream msgStr;
804 #else
805  std::ostringstream msgStr;
806 #endif
807  msgStr << " input array size: " << nElements << " required " << nRows*repeat();
808  String msg(msgStr.str());
809 
810  throw Column::InsufficientElements(msg);
811  }
812 
813  if (nullValue)
814  {
815  if (fits_write_colnull(fitsPointer(),abs(type()),index(),firstRow,
816  1,nElements,data,nullValue,&status)) throw FitsError(status);
817  }
818  else
819  {
820  if (fits_write_col(fitsPointer(),abs(type()),index(),firstRow,
821  1,nElements,data,&status)) throw FitsError(status);
822  }
823 
824  parent()->updateRows();
825  }
826 
827  template <typename T>
828  void ColumnVectorData<T>::insertRows (long first, long number)
829  {
830  typename std::vector<std::valarray<T> >::iterator in;
831  if (first !=0)
832  {
833  in = m_data.begin()+first;
834  }
835  else
836  {
837  in = m_data.begin();
838  }
839 
840  // non-throwing operations.
841  m_data.insert(in,number,std::valarray<T>(T(),0));
842  }
843 
844  template <typename T>
845  void ColumnVectorData<T>::deleteRows (long first, long number)
846  {
847  // the following is an ugly workaround for a bug in g++ v3.0 that
848  // does not erase vector elements cleanly in this case.
849 
850  long N = static_cast<long>(m_data.size());
851  size_t newSize = static_cast<size_t>(N - number);
852  std::vector<std::valarray<T> > __tmp(newSize);
853 
854  long lastDeleted( number + first - 1 );
855  long firstDeleted(first);
856  long count(0);
857  {
858  for (long j = 1; j <= N; ++j)
859  {
860  if ( (j - firstDeleted)*(lastDeleted - j) >= 0 )
861  { ++count;
862  }
863  else
864  {
865  __tmp[j - 1 - count].resize(m_data[j - 1].size());
866  __tmp[j - 1 - count] = m_data[j - 1];
867  }
868  }
869  }
870 
871  m_data.clear();
872  m_data.resize(newSize);
873  {
874  for (size_t j = 0; j < newSize; ++j)
875  {
876  m_data[j].resize(__tmp[j].size());
877  m_data[j] = __tmp[j];
878  }
879  }
880  }
881 
882  template <typename T>
883  void ColumnVectorData<T>::setDataLimits (T* limits)
884  {
885  m_minLegalValue = limits[0];
886  m_maxLegalValue = limits[1];
887  m_minDataValue = std::max(limits[2],limits[0]);
888  m_maxDataValue = std::min(limits[3],limits[1]);
889  }
890 
891  template <typename T>
892  void ColumnVectorData<T>::doWrite (T* array, long row, long rowSize, long firstElem, T* nullValue)
893  {
894  int status(0);
895  // internal functioning of write_colnull forbids its use for writing
896  // variable width columns. If a nullvalue argument was supplied it will
897  // be ignored.
898  if ( !varLength())
899  {
900  if (fits_write_colnull(fitsPointer(),type(),index(),row, firstElem, rowSize,
901  array, nullValue,&status)) throw FitsError(status);
902  }
903  else
904  {
905  if (fits_write_col(fitsPointer(),abs(type()),index(),row,firstElem,rowSize,
906  array,&status)) throw FitsError(status);
907 
908  }
909  }
910 
911  // Additional Declarations
912 
913  // all functions that operate on complex data that call cfitsio
914  // need to be specialized. The signature containing complex<T>* objects
915  // is unfortunate, perhaps, for this purpose, but the user will access
916  // rw operations through standard library containers.
917 
918 
919 
920 
921 
922 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT
923 template <>
924 inline void ColumnVectorData<complex<float> >::setDataLimits (complex<float>* limits)
925  {
926  m_minLegalValue = limits[0];
927  m_maxLegalValue = limits[1];
928  m_minDataValue = limits[2];
929  m_maxDataValue = limits[3];
930  }
931 #else
932 template <>
933  void
934  ColumnVectorData<complex<float> >::setDataLimits (complex<float>* limits);
935 #endif
936 
937 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT
938 template <>
939 inline void ColumnVectorData<complex<double> >::setDataLimits (complex<double>* limits)
940  {
941  m_minLegalValue = limits[0];
942  m_maxLegalValue = limits[1];
943  m_minDataValue = limits[2];
944  m_maxDataValue = limits[3];
945  }
946 #else
947  template <>
948  void
949  ColumnVectorData<complex<double> >::setDataLimits (complex<double>* limits);
950 #endif
951 
952 
953 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT
954  template <>
955  inline void ColumnVectorData<std::complex<float> >::readColumnData(long firstRow,
956  long nelements, long firstElem, std::complex<float>* null )
957  {
958  int status=0;
959  float nulval (0);
960  FITSUtil::auto_array_ptr<float> pArray(new float[2*nelements]);
961  float* array = pArray.get();
962  int anynul(0);
963 
964  if (fits_read_col_cmp(fitsPointer(),index(),firstRow, firstElem,
965  nelements,nulval,array,&anynul,&status) ) throw FitsError(status);
966 
967  if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows());
968 
969  std::valarray<std::complex<float> > readData(nelements);
970  for (long j = 0; j < nelements; ++j)
971  {
972  readData[j] = std::complex<float>(array[2*j],array[2*j+1]);
973  }
974  size_t countRead = 0;
975  const size_t ONE = 1;
976 
977  if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows());
978  size_t vectorSize(0);
979  if (!varLength())
980  {
981  vectorSize = std::max(repeat(),ONE); // safety check.
982  }
983  else
984  {
985  // assume that the user specified the correct length for
986  // variable columns. This should be ok since readVariableColumns
987  // uses fits_read_descripts to return this information from the
988  // fits pointer, and this is passed as nelements here.
989  vectorSize = nelements;
990  }
991  size_t n = nelements;
992 
993  int i = firstRow;
994  int ii = i - 1;
995  while ( countRead < n)
996  {
997  std::valarray<complex<float> >& current = m_data[ii];
998  if (current.size() != vectorSize) current.resize(vectorSize,0.);
999  int elementsInFirstRow = vectorSize-firstElem + 1;
1000  bool lastRow = ( (nelements - countRead) < vectorSize);
1001  if (lastRow)
1002  {
1003  int elementsInLastRow = nelements - countRead;
1004  std::copy(&readData[countRead],&readData[0]+nelements,&current[0]);
1005  countRead += elementsInLastRow;
1006  }
1007  // what to do with complete rows. if firstElem == 1 the
1008  else
1009  {
1010  if (firstElem == 1 || (firstElem > 1 && i > firstRow) )
1011  {
1012  current = readData[std::slice(vectorSize*(ii-firstRow)+
1013  elementsInFirstRow,vectorSize,1)];
1014  ++ii;
1015  ++i;
1016  countRead += vectorSize;
1017  }
1018  else
1019  {
1020  if (i == firstRow)
1021  {
1022  std::copy(&readData[0],&readData[0]+elementsInFirstRow,
1023  &current[firstElem]);
1024  countRead += elementsInFirstRow;
1025  ++i;
1026  ++ii;
1027  }
1028  }
1029  }
1030  }
1031  }
1032 #else
1033 template <>
1034 void ColumnVectorData<complex<float> >::readColumnData(long firstRow,
1035  long nelements,
1036  long firstElem, complex<float>* null);
1037 #endif
1038 
1039 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT
1040  template <>
1041  inline void ColumnVectorData<complex<double> >::readColumnData (long firstRow,
1042  long nelements,long firstElem,
1043  complex<double>* nullValue)
1044  {
1045 
1046  // duplicated for each complex type to work around imagined or
1047  // actual compiler deficiencies.
1048  int status=0;
1049  double nulval (0);
1050  FITSUtil::auto_array_ptr<double> pArray(new double[2*nelements]);
1051  double* array = pArray.get();
1052  int anynul(0);
1053 
1054  if (fits_read_col_dblcmp(fitsPointer(),index(),firstRow, firstElem,
1055  nelements,nulval,array,&anynul,&status) ) throw FitsError(status);
1056 
1057  if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows());
1058 
1059  std::valarray<std::complex<double> > readData(nelements);
1060  for (long j = 0; j < nelements; ++j)
1061  {
1062  readData[j] = std::complex<double>(array[2*j],array[2*j+1]);
1063  }
1064  size_t countRead = 0;
1065  const size_t ONE = 1;
1066 
1067  if (m_data.size() != static_cast<size_t>(rows())) m_data.resize(rows());
1068  size_t vectorSize(0);
1069  if (!varLength())
1070  {
1071  vectorSize = std::max(repeat(),ONE); // safety check.
1072  }
1073  else
1074  {
1075  // assume that the user specified the correct length for
1076  // variable columns. This should be ok since readVariableColumns
1077  // uses fits_read_descripts to return this information from the
1078  // fits pointer, and this is passed as nelements here.
1079  vectorSize = nelements;
1080  }
1081  size_t n = nelements;
1082 
1083  int i = firstRow;
1084  int ii = i - 1;
1085  while ( countRead < n)
1086  {
1087  std::valarray<std::complex<double> >& current = m_data[ii];
1088  if (current.size() != vectorSize) current.resize(vectorSize,0.);
1089  int elementsInFirstRow = vectorSize-firstElem + 1;
1090  bool lastRow = ( (nelements - countRead) < vectorSize);
1091  if (lastRow)
1092  {
1093  int elementsInLastRow = nelements - countRead;
1094  std::copy(&readData[countRead],&readData[0]+nelements,&current[0]);
1095  countRead += elementsInLastRow;
1096  }
1097  // what to do with complete rows. if firstElem == 1 the
1098  else
1099  {
1100  if (firstElem == 1 || (firstElem > 1 && i > firstRow) )
1101  {
1102  current = readData[std::slice(vectorSize*(ii-firstRow)+
1103  elementsInFirstRow,vectorSize,1)];
1104  ++ii;
1105  ++i;
1106  countRead += vectorSize;
1107  }
1108  else
1109  {
1110  if (i == firstRow)
1111  {
1112  std::copy(&readData[0],&readData[0]+elementsInFirstRow,
1113  &current[firstElem]);
1114  countRead += elementsInFirstRow;
1115  ++i;
1116  ++ii;
1117  }
1118  }
1119  }
1120  }
1121  }
1122 #else
1123 template <>
1124 void ColumnVectorData<complex<double> >::readColumnData (long firstRow,
1125  long nelements,
1126  long firstElem, complex<double>* null);
1127 #endif
1128 
1129 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT
1130  template <>
1131  inline void ColumnVectorData<complex<float> >::writeFixedArray
1132  (complex<float>* data, long nElements, long nRows, long firstRow,
1133  complex<float>* nullValue)
1134  {
1135 
1136  int status(0);
1137 
1138  // check for sanity of inputs, then write to file.
1139  // this function writes only complete rows to a table with
1140  // fixed width rows.
1141 
1142 
1143  if ( nElements < nRows*static_cast<long>(repeat()) )
1144  {
1145 #ifdef SSTREAM_DEFECT
1146  std::ostrstream msgStr;
1147 #else
1148  std::ostringstream msgStr;
1149 #endif
1150  msgStr << " input array size: " << nElements
1151  << " required " << nRows*repeat();
1152 #ifdef SSTREAM_DEFECT
1153  msgStr << std::ends;
1154 #endif
1155 
1156 
1157  String msg(msgStr.str());
1158 
1159  throw Column::InsufficientElements(msg);
1160  }
1161 
1162  FITSUtil::auto_array_ptr<float> realData(new float[2*nElements]);
1163 
1164  for (int j = 0; j < nElements; ++j)
1165  {
1166  realData[2*j] = data[j].real();
1167  realData[2*j+1] = data[j].imag();
1168  }
1169 
1170 
1171 
1172  if (fits_write_col_cmp(fitsPointer(),index(),firstRow,
1173  1,nElements,realData.get(),&status)) throw FitsError(status);
1174 
1175  parent()->updateRows();
1176  }
1177 #else
1178 template <>
1179 void ColumnVectorData<complex<float> >::writeFixedArray
1180  (complex<float>* data, long nElements, long nRows, long firstRow, std::complex<float>* null);
1181 #endif
1182 
1183 #if SPEC_TEMPLATE_IMP_DEFECT || SPEC_TEMPLATE_DECL_DEFECT
1184  template <>
1185  inline void ColumnVectorData<complex<double> >::writeFixedArray
1186  (complex<double>* data, long nElements, long nRows, long firstRow,
1187  complex<double>* nullValue)
1188  {
1189  int status(0);
1190 
1191  // check for sanity of inputs, then write to file.
1192  // this function writes only complete rows to a table with
1193  // fixed width rows.
1194 
1195 
1196  if ( nElements < nRows*static_cast<long>(repeat()) )
1197  {
1198 #ifdef SSTREAM_DEFECT
1199  std::ostrstream msgStr;
1200 #else
1201  std::ostringstream msgStr;
1202 #endif
1203  msgStr << " input array size: " << nElements
1204  << " required " << nRows*repeat();
1205 #ifdef SSTREAM_DEFECT
1206  msgStr << std::ends;
1207 #endif
1208 
1209  String msg(msgStr.str());
1210 
1211  throw Column::InsufficientElements(msg);
1212  }
1213 
1214  FITSUtil::auto_array_ptr<double> realData(new double[2*nElements]);
1215 
1216  for (int j = 0; j < nElements; ++j)
1217  {
1218  realData[2*j] = data[j].real();
1219  realData[2*j+1] = data[j].imag();
1220  }
1221 
1222 
1223 
1224  if (fits_write_col_dblcmp(fitsPointer(),index(),firstRow,
1225  1,nElements,realData.get(),&status)) throw FitsError(status);
1226 
1227  parent()->updateRows();
1228 
1229  }
1230 #else
1231 template <>
1232 void ColumnVectorData<complex<double> >::writeFixedArray
1233  (complex<double>* data, long nElements, long nRows, long firstRow,
1234  std::complex<double>* null);
1235 #endif
1236 
1237 #ifdef SPEC_TEMPLATE_DECL_DEFECT
1238  template <>
1239  inline void
1240  ColumnVectorData<std::complex<float> >::doWrite
1241  (std::complex<float>* data, long row, long rowSize, long firstElem, std::complex<float>* nullValue )
1242  {
1243  int status(0);
1244  FITSUtil::auto_array_ptr<float> carray( new float[2*rowSize]);
1245  for ( long j = 0 ; j < rowSize; ++ j)
1246  {
1247  carray[2*j] = data[j].real();
1248  carray[2*j + 1] = data[j].imag();
1249  }
1250  if (fits_write_col_cmp(fitsPointer(),index(),row,firstElem,rowSize,
1251  carray.get(),&status)) throw FitsError(status);
1252  }
1253 
1254 
1255  template <>
1256  inline void
1257  ColumnVectorData<std::complex<double> >::doWrite
1258  (std::complex<double>* data, long row, long rowSize, long firstElem, std::complex<double>* nullValue )
1259  {
1260  int status(0);
1261  FITSUtil::auto_array_ptr<double> carray( new double[2*rowSize]);
1262  for ( long j = 0 ; j < rowSize; ++ j)
1263  {
1264  carray[2*j] = data[j].real();
1265  carray[2*j + 1] = data[j].imag();
1266  }
1267  if (fits_write_col_dblcmp(fitsPointer(),index(),row,firstElem,rowSize,
1268  carray.get(),&status)) throw FitsError(status);
1269 
1270  }
1271 
1272 #else
1273 template<>
1274 void
1275 ColumnVectorData<complex<float> >::doWrite
1276  ( complex<float>* data, long row, long rowSize, long firstElem, complex<float>* nullValue);
1277 
1278 template<>
1279 void
1280 ColumnVectorData<complex<double> >::doWrite
1281  ( complex<double>* data, long row, long rowSize, long firstElem, complex<double>* nullValue );
1282 #endif
1283 } // namespace CCfits
1284 
1285 
1286 #endif
fitsfile * fitsPointer()
fits pointer corresponding to fits file containing column data.
Definition: Column.cxx:264
int index() const
get the Column index (the n in TTYPEn etc).
Definition: Column.h:1357
const String & format() const
return TFORMn keyword
Definition: Column.h:1519
Table * parent() const
return a pointer to the Table which owns this Column
Definition: Column.cxx:312
int rows() const
return the number of rows in the table.
Definition: Column.cxx:275
void updateRows()
update the number of rows in the table
Definition: Table.cxx:301
void makeHDUCurrent()
make HDU containing this the current HDU of the fits file.
Definition: Column.cxx:270
Definition: MSconfig.h:123
virtual std::ostream & put(std::ostream &s) const
internal implementation of << operator.
Definition: Column.cxx:302
static bool verboseMode()
return verbose setting for library
Definition: FITS.h:891
const String & comment() const
retrieve comment for Column
Definition: Column.h:1514
const String & dimen() const
return TDIMn keyword
Definition: Column.h:1424
bool isRead() const
flag set to true if the entire column data has been read from disk
Definition: Column.h:1367
Namespace enclosing all CCfits classes and globals definitions.
Definition: AsciiTable.cxx:26
ValueType
CCfits value types and their CFITSIO equivalents (in caps)
Definition: CCfits.h:81
Column(const Column &right)
copy constructor, used in copying Columns to standard library containers.
Definition: Column.cxx:171
const String & unit() const
get units of data in Column (TUNITn keyword)
Definition: Column.h:1524
const String & name() const
return name of Column (TTYPEn keyword)
Definition: Column.h:1529
virtual void readData(long firstRow, long nelements, long firstElem=1)=0
Read (or reread) data from the disk into the Column object&#39;s internal arrays.
size_t repeat() const
get the repeat count for the rows
Definition: Column.h:1387
ValueType type() const
returns the data type of the column
Definition: Column.h:1434
bool varLength() const
boolean, set to true if Column has variable length vector rows.
Definition: Column.h:1392