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