vector.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_VECTOR_H
24 #define O2SCL_VECTOR_H
25 
26 /** \file vector.h
27  \brief Assorted generic vector functions
28 
29  This file contains a set of template functions which can be
30  applied to almost any vector or matrix type which allow element
31  access through <tt>operator[]</tt> (for vectors) or
32  <tt>operator(,)</tt> for matrices. Detailed requirements on the
33  template parameters are given in the functions below.
34 
35  For a general discussion of vectors and matrices in \o2, see the
36  \ref vecmat_section of the User's Guide.
37 
38  For statistics operations not included here, see \ref vec_stats.h
39  in the directory \c src/other . Also related are the matrix output
40  functions, \ref o2scl::matrix_out(), which is defined in \ref
41  columnify.h because they utilize the class \ref o2scl::columnify to
42  format the output.
43 
44  For functions which search for a value in an ordered (either
45  increasing or decreasing) vector, see the class \ref
46  o2scl::search_vec .
47 
48  \future Create a matrix transpose copy function?
49  \future Create matrix swap row and column functions
50 */
51 #include <iostream>
52 #include <cmath>
53 #include <string>
54 #include <fstream>
55 #include <sstream>
56 
57 #include <gsl/gsl_vector.h>
58 #include <gsl/gsl_sys.h>
59 #include <gsl/gsl_matrix.h>
60 
61 #include <o2scl/misc.h>
62 #include <o2scl/uniform_grid.h>
63 
64 namespace o2scl {
65 
66  /** \brief A simple convenience wrapper for GSL vector objects
67 
68  \warning This uses typecasts on externally allocated GSL
69  pointers and is not safe or fully const-correct.
70  */
72  const double *d;
73  size_t sz;
74  public:
75  gsl_vector_wrap(gsl_vector *m) {
76  d=(const double *)m->data;
77  sz=m->size;
78  }
79  double operator[](size_t i) const {
80  return d[i];
81  }
82  size_t size() {
83  return sz;
84  }
85  };
86 
87  /** \brief A simple convenience wrapper for GSL matrix objects
88 
89  \warning This uses typecasts on externally allocated GSL
90  pointers and is not safe or fully const-correct.
91  */
93  const double *d;
94  size_t sz1;
95  size_t sz2;
96  size_t tda;
97  public:
98  gsl_matrix_wrap(gsl_matrix *m) {
99  d=(const double *)m->data;
100  sz1=m->size1;
101  sz2=m->size2;
102  tda=m->tda;
103  }
104  double operator()(size_t i, size_t j) const {
105  return *(d+i*tda+j);
106  }
107  size_t size1() {
108  return sz1;
109  }
110  size_t size2() {
111  return sz2;
112  }
113  };
114 
115  /// \name Copying vectors and matrices
116  //@{
117  /** \brief Simple vector copy
118 
119  Copy \c src to \c dest, resizing \c dest if it is too small
120  to hold <tt>src.size()</tt> elements.
121 
122  This function will work for any classes \c vec_t and
123  \c vec2_t which have suitably defined <tt>operator[]</tt>,
124  <tt>size()</tt>, and <tt>resize()</tt> methods.
125  */
126  template<class vec_t, class vec2_t>
127  void vector_copy(const vec_t &src, vec2_t &dest) {
128  size_t N=src.size();
129  if (dest.size()<N) dest.resize(N);
130  size_t i, m=N%4;
131  for(i=0;i<m;i++) {
132  dest[i]=src[i];
133  }
134  for(i=m;i+3<N;i+=4) {
135  dest[i]=src[i];
136  dest[i+1]=src[i+1];
137  dest[i+2]=src[i+2];
138  dest[i+3]=src[i+3];
139  }
140  return;
141  }
142 
143  /** \brief Simple vector copy of the first N elements
144 
145  Copy the first \c N elements of \c src to \c dest.
146  It is assumed that the memory allocation for \c dest
147  has already been performed.
148 
149  This function will work for any class <tt>vec2_t</tt> which has
150  an operator[] which returns a reference to the corresponding
151  element and class <tt>vec_t</tt> with an operator[] which
152  returns either a reference or the value of the corresponding
153  element.
154  */
155  template<class vec_t, class vec2_t>
156  void vector_copy(size_t N, const vec_t &src, vec2_t &dest) {
157  size_t i, m=N%4;
158  for(i=0;i<m;i++) {
159  dest[i]=src[i];
160  }
161  for(i=m;i+3<N;i+=4) {
162  dest[i]=src[i];
163  dest[i+1]=src[i+1];
164  dest[i+2]=src[i+2];
165  dest[i+3]=src[i+3];
166  }
167  return;
168  }
169 
170  /** \brief Simple matrix copy
171 
172  Copy \c src to \c dest, resizing \c dest if it is too small.
173 
174  This function will work for any classes \c mat_t and
175  \c mat2_t which have suitably defined <tt>operator()</tt>,
176  <tt>size()</tt>, and <tt>resize()</tt> methods.
177  */
178  template<class mat_t, class mat2_t>
179  void matrix_copy(mat_t &src, mat2_t &dest) {
180  size_t m=src.size1();
181  size_t n=src.size2();
182  if (dest.size1()<m || dest.size2()<n) dest.resize(m,n);
183  for(size_t i=0;i<m;i++) {
184  for(size_t j=0;j<n;j++) {
185  dest(i,j)=src(i,j);
186  }
187  }
188  }
189 
190  /** \brief Simple matrix copy of the first \f$ (M,N) \f$
191  matrix elements
192 
193  Copy the first <tt>(M,N)</tt> elements of \c src to \c dest. It
194  is assumed that the memory allocation for \c dest has already
195  been performed.
196 
197  This function will work for any class <tt>vec2_t</tt> which has
198  an operator[][] which returns a reference to the corresponding
199  element and class <tt>vec_t</tt> with an operator[][] which
200  returns either a reference or the value of the corresponding
201  element.
202  */
203  template<class mat_t, class mat2_t>
204  void matrix_copy(size_t M, size_t N, mat_t &src, mat2_t &dest) {
205  for(size_t i=0;i<M;i++) {
206  for(size_t j=0;j<N;j++) {
207  dest(i,j)=src(i,j);
208  }
209  }
210  }
211  //@}
212 
213  /// \name Tranpositions
214  //@{
215  /** \brief Simple transpose
216 
217  Copy the transpose of \c src to \c dest, resizing \c dest if it
218  is too small.
219 
220  This function will work for any classes \c mat_t and
221  \c mat2_t which have suitably defined <tt>operator()</tt>,
222  <tt>size()</tt>, and <tt>resize()</tt> methods.
223  */
224  template<class mat_t, class mat2_t>
225  void matrix_transpose(mat_t &src, mat2_t &dest) {
226  size_t m=src.size1();
227  size_t n=src.size2();
228  if (dest.size1()<n || dest.size2()<m) dest.resize(n,m);
229  for(size_t i=0;i<m;i++) {
230  for(size_t j=0;j<n;j++) {
231  dest(i,j)=src(j,i);
232  }
233  }
234  }
235 
236  /** \brief Simple transpose of the first \f$ (m,n) \f$
237  matrix elements
238 
239  Copy the transpose of the first \c m rows and the first \c cols
240  of the matrix \c src into the matrix \c dest
241 
242  This function will work for any classes \c mat_t and \c mat2_t
243  which has a suitably defined <tt>operator()</tt> method.
244  */
245  template<class mat_t, class mat2_t>
246  void matrix_transpose(size_t m, size_t n, mat_t &src, mat2_t &dest) {
247  for(size_t i=0;i<m;i++) {
248  for(size_t j=0;j<n;j++) {
249  dest(i,j)=src(j,i);
250  }
251  }
252  }
253 
254  /** \brief Simple transpose in-place
255 
256  Transpose the matrix \c src . If the matrix is not square,
257  only the upper-left square part of the matrix will be
258  transposed.
259 
260  This function will work for any classes \c mat_t and
261  \c mat2_t which have suitably defined <tt>operator()</tt>,
262  <tt>size()</tt>, and <tt>resize()</tt> methods.
263  */
264  template<class mat_t, class data_t>
265  void matrix_transpose(mat_t &src) {
266  size_t m=src.size1();
267  size_t n=src.size2();
268  // Choose the smaller of n and m
269  if (m<n) n=m;
270  data_t tmp;
271  for(size_t i=0;i<n;i++) {
272  for(size_t j=0;j<n;j++) {
273  tmp=src(i,j);
274  src(i,j)=src(j,i);
275  src(j,i)=tmp;
276  }
277  }
278  }
279 
280  /** \brief Simple in-place transpose of the first \f$ (m,n) \f$
281  matrix elements
282 
283  Copy the transpose of the first \c m rows and the first \c cols
284  of the matrix \c src into the matrix \c dest
285 
286  This function will work for any classes \c mat_t and \c mat2_t
287  which has a suitably defined <tt>operator()</tt> method.
288  */
289  template<class mat_t, class data_t>
290  void matrix_transpose(size_t m, size_t n, mat_t &src) {
291  data_t tmp;
292  for(size_t i=0;i<m;i++) {
293  for(size_t j=0;j<n;j++) {
294  tmp=src(i,j);
295  src(i,j)=src(j,i);
296  src(j,i)=tmp;
297  }
298  }
299  }
300  //@}
301 
302  /// \name Upper and lower triangular functions
303  //@{
304  /** \brief Simple test that a matrix is lower triangular
305  */
306  template<class mat_t> bool matrix_is_lower(mat_t &src) {
307  size_t m=src.size1();
308  size_t n=src.size2();
309  bool ret=true;
310  for(size_t i=0;i<m;i++) {
311  for(size_t j=i+1;j<n;j++) {
312  if (src(i,j)!=0.0) return false;
313  }
314  }
315  return ret;
316  }
317 
318  /** \brief Simple test that a matrix is upper triangular
319  */
320  template<class mat_t> bool matrix_is_upper(mat_t &src) {
321  size_t m=src.size1();
322  size_t n=src.size2();
323  bool ret=true;
324  for(size_t j=0;j<n;j++) {
325  for(size_t i=j+1;i<m;i++) {
326  if (src(i,j)!=0.0) return false;
327  }
328  }
329  return ret;
330  }
331 
332  /** \brief Make a matrix lower triangular by setting the upper
333  triangular entries to zero
334  */
335  template<class mat_t> void matrix_make_lower(mat_t &src) {
336  size_t m=src.size1();
337  size_t n=src.size2();
338  for(size_t i=0;i<m;i++) {
339  for(size_t j=i+1;j<n;j++) {
340  src(i,j)=0.0;
341  }
342  }
343  return;
344  }
345 
346  /** \brief Make a matrix upper triangular by setting the lower
347  triangular entries to zero
348  */
349  template<class mat_t> void matrix_make_upper(mat_t &src) {
350  size_t m=src.size1();
351  size_t n=src.size2();
352  for(size_t j=0;j<n;j++) {
353  for(size_t i=j+1;i<m;i++) {
354  src(i,j)=0.0;
355  }
356  }
357  return;
358  }
359 
360  /** \brief Simple test that a matrix is lower triangular
361  for the first \c m rows and \c n columns
362  */
363  template<class mat_t> bool matrix_is_lower(size_t m, size_t n,
364  mat_t &src) {
365  bool ret=true;
366  for(size_t i=0;i<m;i++) {
367  for(size_t j=i+1;j<n;j++) {
368  if (src(i,j)!=0.0) return false;
369  }
370  }
371  return ret;
372  }
373 
374  /** \brief Simple test that a matrix is upper triangular
375  for the first \c m rows and \c n columns
376  */
377  template<class mat_t> bool matrix_is_upper(size_t m, size_t n,
378  mat_t &src) {
379  bool ret=true;
380  for(size_t j=0;j<n;j++) {
381  for(size_t i=j+1;i<m;i++) {
382  if (src(i,j)!=0.0) return false;
383  }
384  }
385  return ret;
386  }
387 
388  /** \brief Make the first \c m rows and \c n columns of a matrix
389  lower triangular by setting the upper triangular entries to zero
390  */
391  template<class mat_t> void matrix_make_lower(size_t m, size_t n,
392  mat_t &src) {
393  for(size_t i=0;i<m;i++) {
394  for(size_t j=i+1;j<n;j++) {
395  src(i,j)=0.0;
396  }
397  }
398  return;
399  }
400 
401  /** \brief Make the first \c m rows and \c n columns of a matrix
402  upper triangular by setting the lower triangular entries to zero
403  */
404  template<class mat_t> void matrix_make_upper(size_t m, size_t n,
405  mat_t &src) {
406  for(size_t j=0;j<n;j++) {
407  for(size_t i=j+1;i<m;i++) {
408  src(i,j)=0.0;
409  }
410  }
411  return;
412  }
413  //@}
414 
415  /// \name Swapping parts of vectors and matrices
416  //@{
417  /** \brief Swap the first \c N elements of two vectors
418 
419  This function swaps the elements of \c v1 and \c v2, one element
420  at a time.
421  */
422  template<class vec_t, class vec2_t, class data_t>
423  void vector_swap(size_t N, vec_t &v1, vec2_t &v2) {
424  data_t temp;
425  size_t i, m=N%4;
426  for(i=0;i<m;i++) {
427  temp=v1[i];
428  v1[i]=v2[i];
429  v2[i]=temp;
430  }
431  for(i=m;i+3<N;i+=4) {
432  temp=v1[i];
433  v1[i]=v2[i];
434  v2[i]=temp;
435  temp=v1[i+1];
436  v1[i+1]=v2[i+1];
437  v2[i+1]=temp;
438  temp=v1[i+2];
439  v1[i+2]=v2[i+2];
440  v2[i+2]=temp;
441  temp=v1[i+3];
442  v1[i+3]=v2[i+3];
443  v2[i+3]=temp;
444  }
445  return;
446  }
447 
448  /** \brief Swap all elements in two vectors
449 
450  This function swaps the elements of \c v1 and \c v2, one element
451  at a time.
452 
453  \note It is almost always better to use <tt>std::swap</tt>
454  than this function, which is provided only in cases where
455  one knows one is going to be forced to use a vector type
456  without a properly defined <tt>std::swap</tt> method.
457  */
458  template<class vec_t, class vec2_t, class data_t>
459  void vector_swap(vec_t &v1, vec2_t &v2) {
460  size_t N=v1.size();
461  if (v2.size()<N) N=v2.size();
462  data_t temp;
463  size_t i, m=N%4;
464  for(i=0;i<m;i++) {
465  temp=v1[i];
466  v1[i]=v2[i];
467  v2[i]=temp;
468  }
469  for(i=m;i+3<N;i+=4) {
470  temp=v1[i];
471  v1[i]=v2[i];
472  v2[i]=temp;
473  temp=v1[i+1];
474  v1[i+1]=v2[i+1];
475  v2[i+1]=temp;
476  temp=v1[i+2];
477  v1[i+2]=v2[i+2];
478  v2[i+2]=temp;
479  temp=v1[i+3];
480  v1[i+3]=v2[i+3];
481  v2[i+3]=temp;
482  }
483  return;
484  }
485 
486  /** \brief Swap of of the first N elements of two
487  double-precision vectors
488 
489  This function swaps the elements of \c v1 and \c v2, one element
490  at a time.
491  */
492  template<class vec_t, class vec2_t>
493  void vector_swap_double(size_t N, vec_t &v1, vec2_t &v2) {
494  return vector_swap<vec_t,vec2_t,double>(N,v1,v2);
495  }
496 
497  /** \brief Swap of all the elements in two
498  double-precision vectors
499 
500  This function swaps the elements of \c v1 and \c v2, one element
501  at a time.
502 
503  \note It is almost always better to use <tt>std::swap</tt>
504  than this function, which is provided only in cases where
505  one knows one is going to be forced to use a vector type
506  without a properly defined <tt>std::swap</tt> method.
507  */
508  template<class vec_t, class vec2_t>
509  void vector_swap_double(vec_t &v1, vec2_t &v2) {
510  return vector_swap<vec_t,vec2_t,double>(v1,v2);
511  }
512 
513  /** \brief Swap two elements in a vector
514 
515  This function swaps the element \c i and element \c j of vector
516  \c v1.
517  */
518  template<class vec_t, class data_t>
519  void vector_swap(vec_t &v, size_t i, size_t j) {
520  data_t temp=v[i];
521  v[i]=v[j];
522  v[j]=temp;
523  return;
524  }
525 
526  /** \brief Swap two elements in a double-precision vector
527 
528  This function swaps the element \c i and element \c j of vector
529  \c v1.
530 
531  This function is used in \ref o2scl_linalg::QRPT_decomp() .
532  */
533  template<class vec_t>
534  void vector_swap_double(vec_t &v, size_t i, size_t j) {
535  return vector_swap<vec_t,double>(v,i,j);
536  }
537 
538  /** \brief Swap of the first \f$ (M,N) \f$ elements in two
539  matrices
540 
541  This function swaps the elements of \c v1 and \c v2, one element
542  at a time.
543  */
544  template<class mat_t, class mat2_t, class data_t>
545  void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2) {
546  data_t temp;
547  for(size_t i=0;i<M;i++) {
548  for(size_t j=0;j<N;j++) {
549  temp=v1[i][j];
550  v1[i][j]=v2[i][j];
551  v2[i][j]=temp;
552  }
553  }
554  return;
555  }
556 
557  /** \brief Swap of the first \f$ (M,N) \f$ elements in two
558  double-precision matrices
559 
560  This function swaps the elements of \c m1 and \c m2, one element
561  at a time.
562  */
563  template<class mat_t, class mat2_t, class data_t>
564  void matrix_swap_double(size_t M, size_t N, mat_t &m1, mat2_t &m2) {
565  return matrix_swap<mat_t,mat2_t,double>(M,N,m1,m2);
566  }
567 
568  /** \brief Swap two elements in a matrix
569 
570  This function swaps the element <tt>(i1,j1)</tt> and
571  element <tt>(i2,j2)</tt> of matrix \c m1.
572  */
573  template<class mat_t, class data_t>
574  void matrix_swap(mat_t &m, size_t i1, size_t j1, size_t i2, size_t j2) {
575  data_t temp=m(i1,j1);
576  m(i1,j1)=m(i2,j2);
577  m(i2,j2)=temp;
578  return;
579  }
580 
581  /** \brief Swap two elements in a double-precision matrix
582 
583  This function swaps the element \c i and element \c j of matrix
584  \c v1.
585  */
586  template<class mat_t>
587  void matrix_swap_double(mat_t &m, size_t i1, size_t j1,
588  size_t i2, size_t j2) {
589  return matrix_swap<mat_t,double>(m,i1,j1,i2,j2);
590  }
591 
592  /** \brief Swap the first \c M rows of two columns in a matrix
593 
594  This function swaps the element <tt>(i1,j1)</tt> and
595  element <tt>(i2,j2)</tt> of matrix \c m1.
596  */
597  template<class mat_t, class data_t>
598  void matrix_swap_cols(size_t M, mat_t &m, size_t j1, size_t j2) {
599  data_t temp;
600  for(size_t i=0;i<M;i++) {
601  temp=m(i,j1);
602  m(i,j1)=m(i,j2);
603  m(i,j2)=temp;
604  }
605  return;
606  }
607 
608  /** \brief Swap the first \c M rows of two columns in a
609  double-precision matrix
610 
611  This function swaps the element \c i and element \c j of matrix
612  \c v1.
613  */
614  template<class mat_t>
615  void matrix_swap_cols_double(size_t M, mat_t &m, size_t j1, size_t j2) {
616  return matrix_swap_cols<mat_t,double>(M,m,j1,j2);
617  }
618 
619  /** \brief Swap the first \c N columns of two rows in a
620  matrix
621 
622  This function swaps the element <tt>(i1,j1)</tt> and
623  element <tt>(i2,j2)</tt> of matrix \c m1.
624  */
625  template<class mat_t, class data_t>
626  void matrix_swap_rows(size_t N, mat_t &m, size_t i1, size_t i2) {
627  data_t temp;
628  for(size_t j=0;j<N;j++) {
629  temp=m(i1,j);
630  m(i1,j)=m(i2,j);
631  m(i2,j)=temp;
632  }
633  return;
634  }
635 
636  /** \brief Swap the first \c N columns of two rows in a
637  double-precision matrix
638 
639  This function swaps the element \c i and element \c j of matrix
640  \c v1.
641  */
642  template<class mat_t>
643  void matrix_swap_rows_double(size_t N, mat_t &m, size_t i1, size_t i2) {
644  return matrix_swap_rows<mat_t,double>(N,m,i1,i2);
645  }
646  //@}
647 
648  /// \name Sorting vectors
649  //@{
650  /** \brief Provide a downheap() function for vector_sort()
651  */
652  template<class vec_t, class data_t>
653  void sort_downheap(vec_t &data, size_t n, size_t k) {
654 
655  data_t v=data[k];
656 
657  while (k<=n/2) {
658  size_t j=2*k;
659 
660  if (j<n && data[j] < data[j+1]) j++;
661  if (!(v < data[j])) break;
662  data[k]=data[j];
663  k=j;
664  }
665 
666  data[k]=v;
667 
668  return;
669  }
670 
671  /** \brief Sort a vector (in increasing order)
672 
673  This is a generic sorting template function using a heapsort
674  algorithm. It will work for any types \c data_t and \c vec_t for
675  which
676  - \c data_t has a non-const version of <tt>operator=</tt>
677  - \c data_t has a less than operator to compare elements
678  - <tt>vec_t::operator[]</tt> returns a non-const reference
679  to an object of type \c data_t
680 
681  In particular, it will work with the STL template class
682  <tt>std::vector</tt>, and arrays and pointers of numeric,
683  character, and string objects.
684 
685  For example,
686  \code
687  std::string list[3]={"dog","cat","fox"};
688  vector_sort<std::string[3],std::string>(3,list);
689  \endcode
690 
691  \note With this function template alone, the user cannot avoid
692  explicitly specifying the template types for this function
693  because there is no parameter of type \c data_t, and function
694  templates cannot handle default template types. For this
695  reason, the function template \ref o2scl::vector_sort_double() was
696  also created which provides the convenience of not requiring
697  the user to specify the vector template type.
698 
699  \note This sorting routine is not stable, i.e. equal elements
700  have arbtrary final ordering
701 
702  \note If \c n is zero, this function will do nothing and will
703  not call the error handler.
704 
705  This works similarly to the GSL function <tt>gsl_sort_vector()</tt>.
706  */
707  template<class vec_t, class data_t>
708  void vector_sort(size_t n, vec_t &data) {
709 
710  size_t N;
711  size_t k;
712 
713  if (n==0) return;
714 
715  N=n-1;
716  k=N/2;
717  k++;
718  do {
719  k--;
720  sort_downheap<vec_t,data_t>(data,N,k);
721  } while (k > 0);
722 
723  while (N > 0) {
724  data_t tmp=data[0];
725  data[0]=data[N];
726  data[N]=tmp;
727  N--;
728  sort_downheap<vec_t,data_t>(data,N,0);
729  }
730 
731  return;
732  }
733 
734  /** \brief Provide a downheap() function for vector_sort_index()
735  */
736  template<class vec_t, class vec_size_t>
737  void sort_index_downheap(size_t N, const vec_t &data, vec_size_t &order,
738  size_t k) {
739 
740  const size_t pki = order[k];
741 
742  while (k <= N / 2) {
743  size_t j = 2 * k;
744 
745  if (j < N && data[order[j]] < data[order[j + 1]]) {
746  j++;
747  }
748 
749  // [GSL] Avoid infinite loop if nan
750  if (!(data[pki] < data[order[j]])) {
751  break;
752  }
753 
754  order[k] = order[j];
755 
756  k = j;
757  }
758 
759  order[k] = pki;
760 
761  return;
762  }
763 
764  /** \brief Create a permutation which sorts
765  the first \c n elements of a vector (in increasing order)
766 
767  This function takes a vector \c data and arranges a list of
768  indices in \c order, which give a sorted version of the vector.
769  The value <tt>order[i]</tt> gives the index of entry in in \c
770  data which corresponds to the <tt>i</tt>th value in the sorted
771  vector. The vector \c data is unchanged by this function, and
772  the initial values in \c order are ignored. Before calling this
773  function, \c order must already be allocated as a vector of size
774  \c n.
775 
776  For example, after calling this function, a sorted version the
777  vector can be output with
778  \code
779  size_t n=5;
780  double data[5]={3.1,4.1,5.9,2.6,3.5};
781  permutation order(n);
782  vector_sort_index(n,data,order);
783  for(size_t i=0;i<n;i++) {
784  cout << data[order[i]] << endl;
785  }
786  \endcode
787 
788  To create a permutation which stores as its <tt>i</tt>th element,
789  the index of <tt>data[i]</tt> in the sorted vector, you can
790  invert the permutation created by this function.
791 
792  This is a generic sorting template function. It will work for
793  any types \c vec_t and \c vec_size_t for which
794  - \c vec_t has an <tt>operator[]</tt>, and
795  - \c vec_size_t has an <tt>operator[]</tt> which returns
796  a \c size_t .
797  One possible type for \c vec_size_t is \ref o2scl::permutation.
798 
799  This works similarly to the GSL function <tt>gsl_sort_index()</tt>.
800 
801  */
802  template<class vec_t, class vec_size_t>
803  void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order) {
804  size_t N;
805  size_t i, k;
806 
807  if (n == 0) return;
808 
809  // [GSL] Set permutation to identity
810 
811  for (i = 0 ; i < n ; i++) {
812  order[i] = i;
813  }
814 
815  /* [GSL] We have n_data elements, last element is at 'n_data-1',
816  first at '0' Set N to the last element number.
817  */
818  N = n - 1;
819 
820  k = N / 2;
821  // [GSL] Compensate the first use of 'k--'
822  k++;
823  do {
824  k--;
825  sort_index_downheap<vec_t,vec_size_t>(N,data,order,k);
826  } while (k > 0);
827 
828  while (N > 0) {
829 
830  // [GSL] First swap the elements
831  size_t tmp = order[0];
832  order[0] = order[N];
833  order[N] = tmp;
834 
835  // [GSL] Then process the heap
836  N--;
837 
838  sort_index_downheap<vec_t,vec_size_t>(N,data,order,0);
839  }
840 
841  return;
842  }
843 
844  /** \brief Create a permutation which sorts a vector
845  (in increasing order)
846 
847  This function takes a vector \c data and arranges a list of
848  indices in \c order, which give a sorted version of the vector.
849  The value <tt>order[i]</tt> gives the index of entry in in \c
850  data which corresponds to the <tt>i</tt>th value in the sorted
851  vector. The vector \c data is unchanged by this function, and
852  the initial values in \c order are ignored. Before calling this
853  function, \c order must already be allocated as a vector of size
854  \c n.
855 
856  For example, after calling this function, a sorted version the
857  vector can be output with
858  \code
859  size_t n=5;
860  double data[5]={3.1,4.1,5.9,2.6,3.5};
861  permutation order(n);
862  vector_sort_index(n,data,order);
863  for(size_t i=0;i<n;i++) {
864  cout << data[order[i]] << endl;
865  }
866  \endcode
867 
868  To create a permutation which stores as its <tt>i</tt>th element,
869  the index of <tt>data[i]</tt> in the sorted vector, you can
870  invert the permutation created by this function.
871 
872  This is a generic sorting template function. It will work for
873  any types \c vec_t and \c vec_size_t for which
874  - \c vec_t has an <tt>operator[]</tt>, and
875  - \c vec_size_t has an <tt>operator[]</tt> which returns
876  a \c size_t .
877  One possible type for \c vec_size_t is \ref o2scl::permutation.
878 
879  This works similarly to the GSL function <tt>gsl_sort_index()</tt>.
880 
881  */
882  template<class vec_t, class vec_size_t>
883  void vector_sort_index(const vec_t &data, vec_size_t &order) {
884  vector_sort_index(data.size(),data,order);
885  return;
886  }
887 
888  /** \brief Sort a vector of doubles (in increasing order)
889 
890  This function is just a wrapper for
891  \code
892  vector_sort<vec_t,double>(n,data);
893  \endcode
894  See the documentation of \ref o2scl::vector_sort() for more
895  details.
896  */
897  template<class vec_t>
898  void vector_sort_double(size_t n, vec_t &data) {
899  return vector_sort<vec_t,double>(n,data);
900  }
901  //@}
902 
903  /// \name Smallest or largest subset functions
904  //@{
905  /** \brief Find the k smallest entries of the first \c n elements
906  of a vector
907 
908  Given a vector \c data of size \c n, this function sets the
909  first \c k entries of the vector \c smallest to the k smallest
910  entries from vector \c data in ascending order. The vector \c
911  smallest must be allocated beforehand to hold at least \c k
912  elements.
913 
914  This works similarly to the GSL function <tt>gsl_sort_smallest()</tt>.
915 
916  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
917  \f$ k << N \f$.
918 
919  If \c k is zero, then this function does nothing and
920  returns \ref o2scl::success .
921  */
922  template<class vec_t, class data_t>
923  void vector_smallest(size_t n, vec_t &data, size_t k, vec_t &smallest) {
924  if (k>n) {
925  O2SCL_ERR2("Subset length greater than size in ",
926  "function vector_smallest().",exc_einval);
927  }
928  if (k==0 || n==0) {
929  O2SCL_ERR2("Vector size zero or k zero in ",
930  "function vector_smallest().",exc_einval);
931  }
932 
933  // Take the first element
934  size_t j=1;
935  data_t xbound=data[0];
936  smallest[0]=xbound;
937 
938  // Examine the remaining elements
939  for(size_t i=1;i<n;i++) {
940  data_t xi=data[i];
941  if (j<k) {
942  j++;
943  } else if (xi>=xbound) {
944  continue;
945  }
946  size_t i1;
947  for(i1=j-1;i1>0;i1--) {
948  if (xi>smallest[i1-1]) break;
949  smallest[i1]=smallest[i1-1];
950  }
951  smallest[i1]=xi;
952  xbound=smallest[j-1];
953  }
954  return;
955  }
956 
957  /** \brief Find the k smallest entries of a vector
958  of a vector
959 
960  Given a vector \c data, this function sets the first \c k
961  entries of the vector \c smallest to the k smallest entries from
962  vector \c data in ascending order. The vector \c smallest
963  is resized if necessary to hold at least \c k elements.
964 
965  This works similarly to the GSL function <tt>gsl_sort_smallest()</tt>.
966 
967  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
968  \f$ k << N \f$.
969 
970  If \c k is zero, then this function does nothing and
971  returns \ref o2scl::success .
972  */
973  template<class vec_t, class data_t>
974  void vector_smallest(vec_t &data, size_t k, vec_t &smallest) {
975  size_t n=data.size();
976  if (smallest.size()<k) smallest.resize(k);
977  return vector_smallest(n,data,k,smallest);
978  }
979 
980  /** \brief Find the indexes of the k smallest entries among the
981  first \c n entries of a vector
982 
983  Given a vector \c data, this function sets the first \c k
984  entries of the vector \c smallest equal to the indexes of the k
985  smallest entries from vector \c data in ascending order. The
986  vector \c smallest is resized if necessary to hold at least \c k
987  elements.
988 
989  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
990  \f$ k << N \f$.
991 
992  If \c k is zero or \c n is zero or \f$ k > n\f$, then this
993  function calls the error handler.
994  */
995  template<class vec_t, class data_t, class vec_size_t>
996  void vector_smallest_index(size_t n, const vec_t &data, size_t k,
997  vec_size_t &index) {
998  if (k>n) {
999  O2SCL_ERR2("Subset length greater than size in ",
1000  "function vector_smallest_index().",exc_einval);
1001  }
1002  if (k==0 || n==0) {
1003  O2SCL_ERR2("Vector size zero or k zero in ",
1004  "function vector_smallest_index().",exc_einval);
1005  }
1006 
1007  index.resize(k);
1008 
1009  // [GSL] Take the first element
1010  size_t j=1;
1011  data_t xbound=data[0];
1012  index[0]=0;
1013 
1014  // [GSL] Examine the remaining elements
1015  for(size_t i=1;i<n;i++) {
1016  data_t xi=data[i];
1017  if (j<k) {
1018  j++;
1019  } else if (xi>=xbound) {
1020  continue;
1021  }
1022  size_t i1;
1023  for(i1=j-1;i1>0;i1--) {
1024  if (xi>data[index[i1-1]]) break;
1025  index[i1]=index[i1-1];
1026  }
1027  index[i1]=i;
1028  xbound=data[index[j-1]];
1029  }
1030  return;
1031  }
1032 
1033  /** \brief Find the indexes of the k smallest entries of a vector
1034  */
1035  template<class vec_t, class data_t, class vec_size_t>
1036  void vector_smallest_index(const vec_t &data, size_t k,
1037  vec_size_t &index) {
1038  size_t n=data.size();
1039  if (index.size()<k) index.resize(k);
1040  return o2scl::vector_smallest_index<vec_t,data_t,vec_size_t>
1041  (n,data,k,index);
1042  }
1043 
1044  /** \brief Find the k largest entries of the first \c n elements
1045  of a vector
1046 
1047  Given a vector \c data of size \c n this sets the first \c k
1048  entries of the vector \c largest to the k largest entries from
1049  vector \c data in descending order. The vector \c largest must
1050  be allocated beforehand to hold at least \c k elements.
1051 
1052  This works similarly to the GSL function <tt>gsl_sort_largest()</tt>.
1053 
1054  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
1055  \f$ k << N \f$.
1056 
1057  If \c k is zero, then this function does nothing and
1058  returns \ref o2scl::success .
1059  */
1060  template<class vec_t, class data_t>
1061  void vector_largest(size_t n, vec_t &data, size_t k, vec_t &largest) {
1062  if (k>n) {
1063  O2SCL_ERR2("Subset length greater than size in ",
1064  "function vector_largest().",exc_einval);
1065  }
1066  if (k==0 || n==0) {
1067  O2SCL_ERR2("Vector size zero or k zero in ",
1068  "function vector_largest().",exc_einval);
1069  }
1070 
1071  // Take the first element
1072  size_t j=1;
1073  data_t xbound=data[0];
1074  largest[0]=xbound;
1075 
1076  // Examine the remaining elements
1077  for(size_t i=1;i<n;i++) {
1078  data_t xi=data[i];
1079  if (j<k) {
1080  j++;
1081  } else if (xi<=xbound) {
1082  continue;
1083  }
1084  size_t i1;
1085  for(i1=j-1;i1>0;i1--) {
1086  if (xi<largest[i1-1]) break;
1087  largest[i1]=largest[i1-1];
1088  }
1089  largest[i1]=xi;
1090  xbound=largest[j-1];
1091  }
1092  return;
1093  }
1094 
1095  /** \brief Find the k largest entries of a vector
1096  of a vector
1097 
1098  Given a vector \c data, this function sets the first \c k
1099  entries of the vector \c largest to the k largest entries from
1100  vector \c data in ascending order. The vector \c largest
1101  is resized if necessary to hold at least \c k elements.
1102 
1103  This works similarly to the GSL function <tt>gsl_sort_largest()</tt>.
1104 
1105  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
1106  \f$ k << N \f$.
1107 
1108  If \c k is zero, then this function does nothing and
1109  returns \ref o2scl::success .
1110  */
1111  template<class vec_t, class data_t>
1112  void vector_largest(vec_t &data, size_t k, vec_t &largest) {
1113  size_t n=data.size();
1114  if (largest.size()<k) largest.resize(k);
1115  return vector_largest(n,data,k,largest);
1116  }
1117  //@}
1118 
1119  /// \name Vector minimum and maximum functions
1120  //@{
1121  /** \brief Compute the maximum of the first \c n elements of a vector
1122  */
1123  template<class vec_t, class data_t>
1124  data_t vector_max_value(size_t n, const vec_t &data) {
1125 
1126  if (n==0) {
1127  O2SCL_ERR("Sent size=0 to vector_max_value().",exc_efailed);
1128  }
1129  data_t max=data[0];
1130  for(size_t i=1;i<n;i++) {
1131  if (data[i]>max) {
1132  max=data[i];
1133  }
1134  }
1135  return max;
1136  }
1137 
1138  /** \brief Compute the maximum value of a vector
1139  */
1140  template<class vec_t, class data_t>
1141  data_t vector_max_value(const vec_t &data) {
1142 
1143  size_t n=data.size();
1144  if (n==0) {
1145  O2SCL_ERR("Sent empty vector to vector_max_value().",exc_efailed);
1146  }
1147  data_t max=data[0];
1148  for(size_t i=1;i<n;i++) {
1149  if (data[i]>max) {
1150  max=data[i];
1151  }
1152  }
1153  return max;
1154  }
1155 
1156  /** \brief Compute the index which holds the
1157  maximum of the first \c n elements of a vector
1158  */
1159  template<class vec_t, class data_t>
1160  size_t vector_max_index(size_t n, const vec_t &data) {
1161 
1162  if (n==0) {
1163  O2SCL_ERR("Sent size=0 to vector_max_index().",exc_efailed);
1164  }
1165  data_t max=data[0];
1166  size_t ix=0;
1167  for(size_t i=1;i<n;i++) {
1168  if (data[i]>max) {
1169  max=data[i];
1170  ix=i;
1171  }
1172  }
1173  return ix;
1174  }
1175 
1176  /** \brief Compute the maximum of the first \c n elements of a vector
1177  */
1178  template<class vec_t, class data_t>
1179  void vector_max(size_t n, const vec_t &data, size_t &index,
1180  data_t &val) {
1181 
1182  if (n==0) {
1183  O2SCL_ERR("Sent size=0 to vector_max().",exc_efailed);
1184  }
1185  val=data[0];
1186  index=0;
1187  for(size_t i=1;i<n;i++) {
1188  if (data[i]>val) {
1189  val=data[i];
1190  index=i;
1191  }
1192  }
1193  return;
1194  }
1195 
1196  /** \brief Compute the minimum of the first \c n elements of a vector
1197  */
1198  template<class vec_t, class data_t>
1199  data_t vector_min_value(size_t n, const vec_t &data) {
1200 
1201  if (n==0) {
1202  O2SCL_ERR("Sent size=0 to vector_min_value().",exc_efailed);
1203  }
1204  data_t min=data[0];
1205  for(size_t i=1;i<n;i++) {
1206  if (data[i]<min) {
1207  min=data[i];
1208  }
1209  }
1210  return min;
1211  }
1212 
1213  /** \brief Compute the minimum value in a vector
1214  */
1215  template<class vec_t, class data_t>
1216  data_t vector_min_value(const vec_t &data) {
1217 
1218  size_t n=data.size();
1219  if (n==0) {
1220  O2SCL_ERR("Sent empty vector to vector_min_value().",exc_efailed);
1221  }
1222  data_t min=data[0];
1223  for(size_t i=1;i<n;i++) {
1224  if (data[i]<min) {
1225  min=data[i];
1226  }
1227  }
1228  return min;
1229  }
1230 
1231  /** \brief Compute the index which holds the
1232  minimum of the first \c n elements of a vector
1233  */
1234  template<class vec_t, class data_t>
1235  size_t vector_min_index(size_t n, const vec_t &data) {
1236 
1237  if (n==0) {
1238  O2SCL_ERR("Sent size=0 to vector_min_index().",exc_efailed);
1239  }
1240  data_t min=data[0];
1241  size_t ix=0;
1242  for(size_t i=1;i<n;i++) {
1243  if (data[i]<min) {
1244  min=data[i];
1245  ix=i;
1246  }
1247  }
1248  return ix;
1249  }
1250 
1251  /** \brief Compute the minimum of the first \c n elements of a vector
1252  */
1253  template<class vec_t, class data_t>
1254  void vector_min(size_t n, const vec_t &data, size_t &index,
1255  data_t &val) {
1256 
1257  if (n==0) {
1258  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1259  }
1260  val=data[0];
1261  index=0;
1262  for(size_t i=1;i<n;i++) {
1263  if (data[i]<val) {
1264  val=data[i];
1265  index=i;
1266  }
1267  }
1268  return;
1269  }
1270 
1271  /** \brief Compute the minimum and maximum of the first
1272  \c n elements of a vector
1273  */
1274  template<class vec_t, class data_t>
1275  void vector_minmax_value(size_t n, vec_t &data,
1276  data_t &min, data_t &max) {
1277 
1278  if (n==0) {
1279  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1280  }
1281  min=data[0];
1282  max=min;
1283  for(size_t i=1;i<n;i++) {
1284  if (data[i]<min) {
1285  min=data[i];
1286  }
1287  if (data[i]>max) {
1288  max=data[i];
1289  }
1290  }
1291  return;
1292  }
1293 
1294  /** \brief Compute the minimum and maximum of the first
1295  \c n elements of a vector
1296  */
1297  template<class vec_t, class data_t>
1298  void vector_minmax_index(size_t n, vec_t &data,
1299  size_t &ix_min, size_t &ix_max) {
1300 
1301  if (n==0) {
1302  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1303  }
1304  data_t min=data[0];
1305  data_t max=min;
1306  ix_min=0;
1307  ix_max=0;
1308  for(size_t i=1;i<n;i++) {
1309  if (data[i]<min) {
1310  min=data[i];
1311  ix_min=i;
1312  }
1313  if (data[i]>max) {
1314  max=data[i];
1315  ix_max=i;
1316  }
1317  }
1318  return;
1319  }
1320 
1321  /** \brief Compute the minimum and maximum of the first
1322  \c n elements of a vector
1323  */
1324  template<class vec_t, class data_t>
1325  void vector_minmax(size_t n, vec_t &data,
1326  size_t &ix_min, data_t &min,
1327  size_t &ix_max, data_t &max) {
1328 
1329  if (n==0) {
1330  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1331  }
1332  min=data[0];
1333  max=min;
1334  ix_min=0;
1335  ix_max=0;
1336  for(size_t i=1;i<n;i++) {
1337  if (data[i]<min) {
1338  min=data[i];
1339  ix_min=i;
1340  }
1341  if (data[i]>max) {
1342  max=data[i];
1343  ix_max=i;
1344  }
1345  }
1346  return;
1347  }
1348  //@}
1349 
1350  /// \name Minima and maxima of vectors through quadratic fit
1351  //@{
1352  /** \brief Maximum of vector by quadratic fit
1353  */
1354  template<class vec_t, class data_t>
1355  data_t vector_max_quad(size_t n, const vec_t &data) {
1356  size_t ix=vector_max_index<vec_t,data_t>(n,data);
1357  if (ix==0) {
1358  return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1359  } else if (ix==n-1) {
1360  return quadratic_extremum_y<data_t>
1361  (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1362  }
1363  return quadratic_extremum_y<data_t>
1364  (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1365  }
1366 
1367  /** \brief Maximum of vector by quadratic fit
1368  */
1369  template<class vec_t, class data_t>
1370  data_t vector_max_quad(size_t n, const vec_t &x, const vec_t &y) {
1371  size_t ix=vector_max_index<vec_t,data_t>(n,y);
1372  if (ix==0 || ix==n-1) return y[ix];
1373  return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1374  y[ix-1],y[ix],y[ix+1]);
1375  }
1376 
1377  /** \brief Location of vector maximum by quadratic fit
1378  */
1379  template<class vec_t, class data_t>
1380  data_t vector_max_quad_loc(size_t n, const vec_t &x, const vec_t &y) {
1381  size_t ix=vector_max_index<vec_t,data_t>(n,y);
1382  if (ix==0 || ix==n-1) return y[ix];
1383  return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1384  y[ix-1],y[ix],y[ix+1]);
1385  }
1386 
1387  /** \brief Minimum of vector by quadratic fit
1388  */
1389  template<class vec_t, class data_t>
1390  data_t vector_min_quad(size_t n, const vec_t &data) {
1391  size_t ix=vector_min_index<vec_t,data_t>(n,data);
1392  if (ix==0) {
1393  return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1394  } else if (ix==n-1) {
1395  return quadratic_extremum_y<data_t>
1396  (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1397  }
1398  return quadratic_extremum_y<data_t>
1399  (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1400  }
1401 
1402  /** \brief Minimum of vector by quadratic fit
1403  */
1404  template<class vec_t, class data_t>
1405  data_t vector_min_quad(size_t n, const vec_t &x, const vec_t &y) {
1406  size_t ix=vector_min_index<vec_t,data_t>(n,y);
1407  if (ix==0 || ix==n-1) return y[ix];
1408  return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1409  y[ix-1],y[ix],y[ix+1]);
1410  }
1411 
1412  /** \brief Location of vector minimum by quadratic fit
1413  */
1414  template<class vec_t, class data_t>
1415  data_t vector_min_quad_loc(size_t n, const vec_t &x, const vec_t &y) {
1416  size_t ix=vector_min_index<vec_t,data_t>(n,y);
1417  if (ix==0 || ix==n-1) return y[ix];
1418  return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1419  y[ix-1],y[ix],y[ix+1]);
1420  }
1421  //@}
1422 
1423  /// \name Matrix minimum and maximum functions
1424  //@{
1425  /** \brief Compute the maximum of the lower-left part of a matrix
1426  */
1427  template<class mat_t, class data_t>
1428  data_t matrix_max_value(size_t m, const size_t n, const mat_t &data) {
1429 
1430  if (m==0 || n==0) {
1431  std::string str=((std::string)"Matrix with zero size (")+
1432  o2scl::itos(m)+","+o2scl::itos(n)+") in "+
1433  "matrix_max_value().";
1434  O2SCL_ERR(str.c_str(),exc_einval);
1435  }
1436  data_t max=data(0,0);
1437  for(size_t i=0;i<m;i++) {
1438  for(size_t j=0;j<n;j++) {
1439  if (data(i,j)>max) {
1440  max=data(i,j);
1441  }
1442  }
1443  }
1444  return max;
1445  }
1446 
1447  /** \brief Compute the maximum of a matrix
1448  */
1449  template<class mat_t, class data_t> data_t
1450  matrix_max_value(const mat_t &data) {
1451  size_t m=data.size1();
1452  size_t n=data.size2();
1453  if (m==0 || n==0) {
1454  std::string str=((std::string)"Matrix with zero size (")+
1455  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1456  "matrix_max_value().";
1457  O2SCL_ERR(str.c_str(),exc_einval);
1458  }
1459  data_t max=data(0,0);
1460  for(size_t i=0;i<m;i++) {
1461  for(size_t j=0;j<n;j++) {
1462  if (data(i,j)>max) {
1463  max=data(i,j);
1464  }
1465  }
1466  }
1467  return max;
1468  }
1469 
1470  /** \brief Compute the maximum of a matrix
1471  */
1472  template<class mat_t> double
1473  matrix_max_value_double(const mat_t &data) {
1474  size_t m=data.size1();
1475  size_t n=data.size2();
1476  if (m==0 || n==0) {
1477  std::string str=((std::string)"Matrix with zero size (")+
1478  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1479  "matrix_max_value_double().";
1480  O2SCL_ERR(str.c_str(),exc_einval);
1481  }
1482  double max=data(0,0);
1483  for(size_t i=0;i<m;i++) {
1484  for(size_t j=0;j<n;j++) {
1485  if (data(i,j)>max) {
1486  max=data(i,j);
1487  }
1488  }
1489  }
1490  return max;
1491  }
1492 
1493  /** \brief Compute the maximum of a matrix and return
1494  the indices of the maximum element
1495  */
1496  template<class mat_t, class data_t>
1497  void matrix_max_index(size_t m, size_t n, const mat_t &data,
1498  size_t &i_max, size_t &j_max, data_t &max) {
1499 
1500  if (m==0 || n==0) {
1501  std::string str=((std::string)"Matrix with zero size (")+
1502  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1503  "matrix_max_index().";
1504  O2SCL_ERR(str.c_str(),exc_einval);
1505  }
1506  max=data(0,0);
1507  i_max=0;
1508  j_max=0;
1509  for(size_t i=0;i<m;i++) {
1510  for(size_t j=0;j<n;j++) {
1511  if (data(i,j)>max) {
1512  max=data(i,j);
1513  i_max=i;
1514  j_max=j;
1515  }
1516  }
1517  }
1518  return;
1519  }
1520 
1521  /** \brief Compute the maximum of a matrix and return
1522  the indices of the maximum element
1523  */
1524  template<class mat_t, class data_t>
1525  void matrix_max_index(const mat_t &data,
1526  size_t &i_max, size_t &j_max, data_t &max) {
1527 
1528  size_t m=data.size1();
1529  size_t n=data.size2();
1530  if (m==0 || n==0) {
1531  std::string str=((std::string)"Matrix with zero size (")+
1532  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1533  "matrix_max_index().";
1534  O2SCL_ERR(str.c_str(),exc_einval);
1535  }
1536  max=data(0,0);
1537  i_max=0;
1538  j_max=0;
1539  for(size_t i=0;i<m;i++) {
1540  for(size_t j=0;j<n;j++) {
1541  if (data(i,j)>max) {
1542  max=data(i,j);
1543  i_max=i;
1544  j_max=j;
1545  }
1546  }
1547  }
1548  return;
1549  }
1550 
1551  /** \brief Compute the minimum of a matrix
1552  */
1553  template<class mat_t, class data_t>
1554  data_t matrix_min_value(size_t m, size_t n, const mat_t &data) {
1555 
1556  if (m==0 || n==0) {
1557  std::string str=((std::string)"Matrix with zero size (")+
1558  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1559  "matrix_min_value().";
1560  O2SCL_ERR(str.c_str(),exc_einval);
1561  }
1562  data_t min=data(0,0);
1563  for(size_t i=0;i<m;i++) {
1564  for(size_t j=0;j<n;j++) {
1565  if (data(i,j)<min) {
1566  min=data(i,j);
1567  }
1568  }
1569  }
1570  return min;
1571  }
1572 
1573  /** \brief Compute the minimum of a matrix
1574  */
1575  template<class mat_t, class data_t>
1576  data_t matrix_min_value(const mat_t &data) {
1577 
1578  size_t m=data.size1();
1579  size_t n=data.size2();
1580  if (m==0 || n==0) {
1581  std::string str=((std::string)"Matrix with zero size (")+
1582  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1583  "matrix_min_value().";
1584  O2SCL_ERR(str.c_str(),exc_einval);
1585  }
1586  data_t min=data(0,0);
1587  for(size_t i=0;i<m;i++) {
1588  for(size_t j=0;j<n;j++) {
1589  if (data(i,j)<min) {
1590  min=data(i,j);
1591  }
1592  }
1593  }
1594  return min;
1595  }
1596 
1597  /** \brief Compute the minimum of a matrix
1598  */
1599  template<class mat_t>
1600  double matrix_min_value_double(const mat_t &data) {
1601 
1602  size_t m=data.size1();
1603  size_t n=data.size2();
1604  if (m==0 || n==0) {
1605  std::string str=((std::string)"Matrix with zero size (")+
1606  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1607  "matrix_min_value().";
1608  O2SCL_ERR(str.c_str(),exc_einval);
1609  }
1610  double min=data(0,0);
1611  for(size_t i=0;i<m;i++) {
1612  for(size_t j=0;j<n;j++) {
1613  if (data(i,j)<min) {
1614  min=data(i,j);
1615  }
1616  }
1617  }
1618  return min;
1619  }
1620 
1621  /** \brief Compute the minimum of a matrix and return
1622  the indices of the minimum element
1623  */
1624  template<class mat_t, class data_t>
1625  void matrix_min_index(size_t n, size_t m, const mat_t &data,
1626  size_t &i_min, size_t &j_min, data_t &min) {
1627 
1628  if (m==0 || n==0) {
1629  std::string str=((std::string)"Matrix with zero size (")+
1630  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1631  "matrix_min_index().";
1632  O2SCL_ERR(str.c_str(),exc_einval);
1633  }
1634  min=data(0,0);
1635  i_min=0;
1636  j_min=0;
1637  for(size_t i=0;i<m;i++) {
1638  for(size_t j=0;j<n;j++) {
1639  if (data(i,j)<min) {
1640  min=data(i,j);
1641  i_min=i;
1642  j_min=j;
1643  }
1644  }
1645  }
1646  return;
1647  }
1648 
1649  /** \brief Compute the minimum of a matrix and return
1650  the indices of the minimum element
1651  */
1652  template<class mat_t, class data_t>
1653  void matrix_min_index(const mat_t &data,
1654  size_t &i_min, size_t &j_min, data_t &min) {
1655 
1656  size_t m=data.size1();
1657  size_t n=data.size2();
1658  if (m==0 || n==0) {
1659  std::string str=((std::string)"Matrix with zero size (")+
1660  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1661  "matrix_min_index().";
1662  O2SCL_ERR(str.c_str(),exc_einval);
1663  }
1664  min=data(0,0);
1665  i_min=0;
1666  j_min=0;
1667  for(size_t i=0;i<m;i++) {
1668  for(size_t j=0;j<n;j++) {
1669  if (data(i,j)<min) {
1670  min=data(i,j);
1671  i_min=i;
1672  j_min=j;
1673  }
1674  }
1675  }
1676  return;
1677  }
1678 
1679  /** \brief Compute the minimum and maximum of a matrix
1680  */
1681  template<class mat_t, class data_t>
1682  void matrix_minmax(size_t n, size_t m, const mat_t &data,
1683  data_t &min, data_t &max) {
1684 
1685  if (n==0 || m==0) {
1686  O2SCL_ERR("Sent size=0 to matrix_min().",exc_efailed);
1687  }
1688  min=data(0,0);
1689  max=data(0,0);
1690  for(size_t i=0;i<n;i++) {
1691  for(size_t j=0;j<m;j++) {
1692  if (data(i,j)<min) {
1693  min=data(i,j);
1694  } else if (data(i,j)>max) {
1695  max=data(i,j);
1696  }
1697  }
1698  }
1699  return;
1700  }
1701 
1702  /** \brief Compute the minimum and maximum of a matrix
1703  */
1704  template<class mat_t, class data_t>
1705  void matrix_minmax(const mat_t &data,
1706  data_t &min, data_t &max) {
1707  return matrix_minmax<mat_t,data_t>
1708  (data.size1(),data.size2(),data,min,max);
1709  }
1710 
1711  /** \brief Compute the minimum and maximum of a matrix and
1712  return their locations
1713  */
1714  template<class mat_t, class data_t>
1715  void matrix_minmax_index(size_t n, size_t m, const mat_t &data,
1716  size_t &i_min, size_t &j_min, data_t &min,
1717  size_t &i_max, size_t &j_max, data_t &max) {
1718 
1719  if (n==0 || m==0) {
1720  O2SCL_ERR2("Sent size=0 to function ",
1721  "matrix_minmax_index().",exc_efailed);
1722  }
1723  min=data(0,0);
1724  i_min=0;
1725  j_min=0;
1726  max=data(0,0);
1727  i_max=0;
1728  j_max=0;
1729  for(size_t i=0;i<n;i++) {
1730  for(size_t j=0;j<m;j++) {
1731  if (data(i,j)<min) {
1732  min=data(i,j);
1733  i_min=i;
1734  j_min=j;
1735  } else if (data(i,j)>max) {
1736  max=data(i,j);
1737  i_max=i;
1738  j_max=j;
1739  }
1740  }
1741  }
1742  return;
1743  }
1744 
1745  /** \brief Compute the sum of matrix elements
1746  */
1747  template<class mat_t, class data_t>
1748  data_t matrix_sum(size_t m, size_t n, const mat_t &data) {
1749 
1750  data_t sum=0.0;
1751  for(size_t i=0;i<m;i++) {
1752  for(size_t j=0;j<n;j++) {
1753  sum+=data(i,j);
1754  }
1755  }
1756  return sum;
1757  }
1758 
1759  /** \brief Compute the sum of matrix elements
1760  */
1761  template<class mat_t, class data_t>
1762  data_t matrix_sum(const mat_t &data) {
1763  return matrix_sum<mat_t,data_t>(data.size1(),data.size2(),data);
1764  }
1765  //@}
1766 
1767  /// \name Searching vectors and matrices
1768  //@{
1769  /** \brief Lookup the value \c x0 in the first \c n elements of
1770  vector \c x
1771 
1772  The function finds the element among the first \c n elements of
1773  \c x which is closest to the value \c x0. It ignores all
1774  elements in \c x which are not finite. If the vector is empty,
1775  or if all of the first \c n elements in \c x are not finite,
1776  then the error handler will be called.
1777 
1778  This function works for all classes \c vec_t where an operator[]
1779  is defined which returns a double (either as a value or a
1780  reference).
1781  */
1782  template<class vec_t>
1783  size_t vector_lookup(size_t n, const vec_t &x, double x0) {
1784  if (n==0) {
1785  O2SCL_ERR("Empty vector in function vector_lookup().",
1786  exc_einval);
1787  return 0;
1788  }
1789  size_t row=0, i=0;
1790  while(!std::isfinite(x[i]) && i<n-1) i++;
1791  if (i==n-1) {
1792  O2SCL_ERR2("Entire vector not finite in ",
1793  "function vector_lookup()",exc_einval);
1794  return 0;
1795  }
1796  double best=x[i], bdiff=fabs(x[i]-x0);
1797  for(;i<n;i++) {
1798  if (std::isfinite(x[i]) && fabs(x[i]-x0)<bdiff) {
1799  row=i;
1800  best=x[i];
1801  bdiff=fabs(x[i]-x0);
1802  }
1803  }
1804  return row;
1805  }
1806 
1807  /** \brief Lookup element \c x0 in vector \c x
1808 
1809  This function finds the element in vector \c x which is closest
1810  to \c x0. It ignores all elements in \c x which are not finite.
1811  If the vector is empty, or if all of the
1812  elements in \c x are not finite, then the error handler will be
1813  called.
1814 
1815  This function works for all classes \c vec_t with a
1816  <tt>size()</tt> method and where an operator[] is defined which
1817  returns a double (either as a value or a reference).
1818  */
1819  template<class vec_t>
1820  size_t vector_lookup(const vec_t &x, double x0) {
1821  return vector_lookup(x.size(),x,x0);
1822  }
1823 
1824  /** \brief Lookup an element in the first $(m,n)$ entries in a matrix
1825 
1826  Return the location <tt>(i,j)</tt> of the element closest to
1827  \c x0.
1828  */
1829  template<class mat_t>
1830  void matrix_lookup(size_t m, size_t n, const mat_t &A,
1831  double x0, size_t &i, size_t &j) {
1832  if (m==0 || n==0) {
1833  O2SCL_ERR("Empty matrix in matrix_lookup().",
1834  exc_einval);
1835  }
1836  double dist=0.0;
1837  bool found_one=false;
1838  for(size_t i2=0;i2<m;i2++) {
1839  for(size_t j2=0;j2<n;j2++) {
1840  if (std::isfinite(A(i,j))) {
1841  if (found_one==false) {
1842  dist=fabs(A(i,j)-x0);
1843  found_one=true;
1844  i=i2;
1845  j=j2;
1846  } else {
1847  if (fabs(A(i,j)-x0)<dist) {
1848  dist=fabs(A(i,j)-x0);
1849  i=i2;
1850  j=j2;
1851  }
1852  }
1853  }
1854  }
1855  }
1856  if (found_one==false) {
1857  O2SCL_ERR2("Entire matrix not finite in ",
1858  "function matrix_lookup()",exc_einval);
1859  }
1860  return;
1861  }
1862 
1863  /** \brief Lookup an element in a matrix
1864 
1865  Return the location <tt>(i,j)</tt> of the element closest to
1866  \c x0.
1867  */
1868  template<class mat_t>
1869  void matrix_lookup(const mat_t &A,
1870  double x0, size_t &i, size_t &j) {
1871  matrix_lookup(A.size1(),A.size2(),x0,i,j);
1872  return;
1873  }
1874 
1875  /** \brief Binary search a part of an increasing vector for
1876  <tt>x0</tt>.
1877 
1878  This function performs a binary search of between
1879  <tt>x[lo]</tt> and <tt>x[hi]</tt>. It returns
1880  - \c lo if \c x0 < <tt>x[lo]</tt>
1881  - \c i if <tt>x[i]</tt> <= \c x0 < <tt>x[i+2]</tt>
1882  for \c lo <= \c i < \c hi
1883  - \c hi-1 if \c x0 >= \c <tt>x[hi-1]</tt>
1884 
1885  This function is designed to find the interval containing \c x0,
1886  not the index of the element closest to x0. To perform the
1887  latter operation, you can use \ref vector_lookup().
1888 
1889  The element at <tt>x[hi]</tt> is never referenced by this
1890  function. The parameter \c hi can be either the index of the
1891  last element (e.g. <tt>n-1</tt> for a vector of size <tt>n</tt>
1892  with starting index <tt>0</tt>), or the index of one element
1893  (e.g. <tt>n</tt> for a vector of size <tt>n</tt> and a starting
1894  index <tt>0</tt>) for a depending on whether or not the user
1895  wants to allow the function to return the index of the last
1896  element.
1897 
1898  This function operates in the same way as
1899  <tt>gsl_interp_bsearch()</tt>.
1900 
1901  The operation of this function is undefined if the data is
1902  not strictly monotonic, i.e. if some of the data elements are
1903  equal.
1904 
1905  This function will call the error handler if \c lo is
1906  greater than \c hi.
1907  */
1908  template<class vec_t, class data_t>
1909  size_t vector_bsearch_inc(const data_t x0, const vec_t &x,
1910  size_t lo, size_t hi) {
1911  if (lo>hi) {
1912  O2SCL_ERR2("Low and high indexes backwards in ",
1913  "function vector_bsearch_inc().",exc_einval);
1914  }
1915  while (hi>lo+1) {
1916  size_t i=(hi+lo)/2;
1917  if (x[i]>x0) {
1918  hi=i;
1919  } else {
1920  lo=i;
1921  }
1922  }
1923 
1924  return lo;
1925  }
1926 
1927  /** \brief Binary search a part of an decreasing vector for
1928  <tt>x0</tt>.
1929 
1930  This function performs a binary search of between
1931  <tt>x[lo]</tt> and <tt>x[hi]</tt> (inclusive). It returns
1932  - \c lo if \c x0 > <tt>x[lo]</tt>
1933  - \c i if <tt>x[i]</tt> >= \c x0 > <tt>x[i+1]</tt>
1934  for \c lo <= \c i < \c hi
1935  - \c hi-1 if \c x0 <= \c <tt>x[hi-1]</tt>
1936 
1937  The element at <tt>x[hi]</tt> is never referenced by this
1938  function. The parameter \c hi can be either the index of the
1939  last element (e.g. <tt>n-1</tt> for a vector of size <tt>n</tt>
1940  with starting index <tt>0</tt>), or the index of one element
1941  (e.g. <tt>n</tt> for a vector of size <tt>n</tt> and a starting
1942  index <tt>0</tt>) for a depending on whether or not the user
1943  wants to allow the function to return the index of the last
1944  element.
1945 
1946  The operation of this function is undefined if the data is
1947  not strictly monotonic, i.e. if some of the data elements are
1948  equal.
1949 
1950  This function will call the error handler if \c lo is
1951  greater than \c hi.
1952  */
1953  template<class vec_t, class data_t>
1954  size_t vector_bsearch_dec(const data_t x0, const vec_t &x,
1955  size_t lo, size_t hi) {
1956  if (lo>hi) {
1957  O2SCL_ERR2("Low and high indexes backwards in ",
1958  "function vector_bsearch_dec().",exc_einval);
1959  }
1960  while (hi>lo+1) {
1961  size_t i=(hi+lo)/2;
1962  if (x[i]<x0) {
1963  hi=i;
1964  } else {
1965  lo=i;
1966  }
1967  }
1968 
1969  return lo;
1970  }
1971 
1972  /** \brief Binary search a part of a monotonic vector for
1973  <tt>x0</tt>.
1974 
1975  This wrapper just calls \ref o2scl::vector_bsearch_inc() or
1976  \ref o2scl::vector_bsearch_dec() depending on the ordering
1977  of \c x.
1978  */
1979  template<class vec_t, class data_t>
1980  size_t vector_bsearch(const data_t x0, const vec_t &x,
1981  size_t lo, size_t hi) {
1982  if (x[lo]<x[hi-1]) {
1983  return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
1984  }
1985  return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
1986  }
1987 
1988  /** \brief Binary search a monotonic vector for
1989  <tt>x0</tt>.
1990 
1991  This function calls \ref o2scl::vector_bsearch_inc() or
1992  \ref o2scl::vector_bsearch_dec() depending on the ordering
1993  of \c x.
1994  */
1995  template<class vec_t, class data_t>
1996  size_t vector_bsearch(const data_t x0, const vec_t &x) {
1997  size_t lo=0;
1998  size_t hi=x.size();
1999  if (x[lo]<x[hi-1]) {
2000  return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
2001  }
2002  return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
2003  }
2004  //@}
2005 
2006  /// \name Ordering and finite tests
2007  //@{
2008  /** \brief Test if the first \c n elements of a vector are
2009  monotonic and increasing or decreasing
2010 
2011  If \c n is zero or one, this function will return 0 without
2012  calling the error handler. If all the vector's elements are equal,
2013  this function will return 3. Otherwise, if the vector is not
2014  monotonic, then this function will return 0. Finally, if the
2015  vector is nondecreasing (increasing or equal intervals), this
2016  function will return 1, and if the vector is nonincreasing
2017  (decreasing or equal intervals), this function will return 2.
2018  This function assumes that simple comparison operators have been
2019  defined for the type of each vector element.
2020  */
2021  template<class vec_t>
2022  int vector_is_monotonic(size_t n, vec_t &data) {
2023 
2024  if (n<1) return 0;
2025  if (n<2) {
2026  if (data[0]==data[1]) {
2027  return 3;
2028  } else if (data[0]<data[1]) {
2029  return 1;
2030  } else {
2031  return 2;
2032  }
2033  }
2034 
2035  // Find first non-flat interval
2036  size_t start=0;
2037  bool done=false;
2038  for(size_t i=0;i<n-1 && done==false;i++) {
2039  if (data[i]!=data[i+1]) {
2040  done=true;
2041  } else {
2042  start++;
2043  }
2044  }
2045 
2046  // If all elements in the vector are equal, or if the only
2047  // distinct element is at the end, then return true.
2048  if (done==false) {
2049  return 3;
2050  }
2051 
2052  if (start==n-2) {
2053  if (data[start]<data[start+1]) return 1;
2054  else return 2;
2055  }
2056 
2057  // Determine if the vector is increasing (inc=true) or decreasing
2058  // (inc=false)
2059  bool inc=true;
2060  if (data[start]>data[start+1]) inc=false;
2061 
2062  if (inc) {
2063  for(size_t i=start+1;i<n-1;i++) {
2064  // If there is one decreasing interval, return false
2065  if (data[i]>data[i+1]) return 0;
2066  }
2067  return 1;
2068  }
2069 
2070  // If there is one increasing interval, return false
2071  for(size_t i=start+1;i<n-1;i++) {
2072  if (data[i]<data[i+1]) return 0;
2073  }
2074  return 2;
2075  }
2076 
2077  /** \brief Test if the first \c n elements of a vector are
2078  monotonic and increasing or decreasing
2079 
2080  If \c n is zero or one, this function will return 0 without
2081  calling the error handler. If all the vector's elements are equal,
2082  this function will return 3. Otherwise, if the vector is not
2083  monotonic, then this function will return 0. Finally, if the
2084  vector is nondecreasing (increasing or equal intervals), this
2085  function will return 1, and if the vector is nonincreasing
2086  (decreasing or equal intervals), this function will return 2.
2087  This function assumes that simple comparison operators have been
2088  defined for the type of each vector element.
2089  */
2090  template<class vec_t> int vector_is_monotonic(vec_t &data) {
2091  return vector_is_monotonic(data.size(),data);
2092  }
2093 
2094  /** \brief Test if the first \c n elements of a vector are
2095  strictly monotonic and determine if they are increasing or decreasing
2096 
2097  If \c n is zero this function will return 0 without calling the
2098  error handler. Also, if the vector is not monotonic, this
2099  function will return 0. If the vector is strictly
2100  monotonic, then this function will return 1 if it is
2101  increasing and 2 if it is decreasing.
2102  */
2103  template<class vec_t>
2104  int vector_is_strictly_monotonic(size_t n, vec_t &data) {
2105 
2106  if (n<1) return 0;
2107 
2108  // Determine if the vector is increasing (inc=true) or decreasing
2109  // (inc=false)
2110  bool inc=true;
2111  if (data[0]==data[1]) {
2112  return 0;
2113  } else if (data[0]>data[1]) {
2114  inc=false;
2115  }
2116 
2117  if (inc) {
2118  for(size_t i=1;i<n-1;i++) {
2119  // If there is one nonincreasing interval, return 0
2120  if (data[i]>=data[i+1]) return 0;
2121  }
2122  return 1;
2123  }
2124 
2125  // If there is one increasing interval, return 0
2126  for(size_t i=1;i<n-1;i++) {
2127  if (data[i]<=data[i+1]) return 0;
2128  }
2129  return 2;
2130  }
2131 
2132  /** \brief Test if the first \c n elements of a vector are
2133  strictly monotonic and determine if they are increasing or decreasing
2134 
2135  If \c n is zero this function will return 0 without calling the
2136  error handler. Also, if the vector is not monotonic, this
2137  function will return 0. If the vector is strictly
2138  monotonic, then this function will return 1 if it is
2139  increasing and 2 if it is decreasing.
2140  */
2141  template<class vec_t>
2142  int vector_is_strictly_monotonic(vec_t &data) {
2143  return vector_is_strictly_monotonic(data.size(),data);
2144  }
2145 
2146  /** \brief Test if the first \c n elements of a vector are finite
2147 
2148  If \c n is zero, this will return true without throwing
2149  an exception.
2150 
2151  The corresponding tests for matrix functions are
2152  in clbas_base.h .
2153  */
2154  template<class vec_t>
2155  bool vector_is_finite(size_t n, vec_t &data) {
2156  for(size_t i=0;i<n;i++) {
2157  if (!std::isfinite(data[i])) return false;
2158  }
2159  return true;
2160  }
2161 
2162  /** \brief Test if a vector is finite
2163 
2164  If \c n is zero, this will return true without throwing
2165  an exception.
2166 
2167  The corresponding tests for matrix functions are
2168  in clbas_base.h .
2169  */
2170  template<class vec_t> bool vector_is_finite(vec_t &data) {
2171  return vector_is_finite(data.size(),data);
2172  }
2173  //@}
2174 
2175  /// \name Miscellaneous mathematical functions
2176  //@{
2177  /** \brief Compute the sum of the first \c n elements of a vector
2178 
2179  If \c n is zero, this will return 0 without throwing
2180  an exception.
2181  */
2182  template<class vec_t, class data_t>
2183  data_t vector_sum(size_t n, vec_t &data) {
2184 
2185  data_t sum=0.0;
2186  for(size_t i=0;i<n;i++) {
2187  sum+=data[i];
2188  }
2189  return sum;
2190  }
2191 
2192  /** \brief Compute the sum of all the elements of a vector
2193 
2194  If the vector has zero size, this will return 0 without
2195  calling the error handler.
2196  */
2197  template<class vec_t, class data_t> data_t vector_sum(vec_t &data) {
2198  data_t sum=0.0;
2199  for(size_t i=0;i<data.size();i++) {
2200  sum+=data[i];
2201  }
2202  return sum;
2203  }
2204 
2205  /** \brief Compute the sum of the first \c n elements of a vector
2206  of double-precision numbers
2207 
2208  If \c n is zero, this will return 0 without throwing
2209  an exception.
2210  */
2211  template<class vec_t>double vector_sum_double(size_t n, vec_t &data) {
2212  double sum=0.0;
2213  for(size_t i=0;i<n;i++) {
2214  sum+=data[i];
2215  }
2216  return sum;
2217  }
2218 
2219  /** \brief Compute the sum of all the elements of a vector
2220  of double-precision numbers
2221 
2222  If the vector has zero size, this will return 0 without
2223  calling the error handler.
2224  */
2225  template<class vec_t> double vector_sum_double(vec_t &data) {
2226  double sum=0.0;
2227  for(size_t i=0;i<data.size();i++) {
2228  sum+=data[i];
2229  }
2230  return sum;
2231  }
2232 
2233  /** \brief Compute the norm of the first \c n entries of a
2234  vector of floating-point (single or double precision) numbers
2235 
2236  This function is a more generic version of
2237  \ref o2scl_cblas::dnrm2 .
2238  */
2239  template<class vec_t, class data_t>
2240  data_t vector_norm(size_t n, const vec_t &x) {
2241 
2242  data_t scale = 0.0;
2243  data_t ssq = 1.0;
2244 
2245  if (n <= 0) {
2246  return 0.0;
2247  } else if (n == 1) {
2248  return fabs(x[0]);
2249  }
2250 
2251  for (size_t i = 0; i < n; i++) {
2252  const data_t xx = x[i];
2253 
2254  if (xx != 0.0) {
2255  const data_t ax = fabs(xx);
2256 
2257  if (scale < ax) {
2258  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2259  scale = ax;
2260  } else {
2261  ssq += (ax / scale) * (ax / scale);
2262  }
2263  }
2264 
2265  }
2266 
2267  return scale * sqrt(ssq);
2268  }
2269 
2270  /** \brief Compute the norm of a vector of floating-point
2271  (single or double precision) numbers
2272  */
2273  template<class vec_t, class data_t> data_t vector_norm(const vec_t &x) {
2274  return vector_norm<vec_t,data_t>(x.size(),x);
2275  }
2276 
2277  /** \brief Compute the norm of the first \c n entries of a
2278  vector of double precision numbers
2279 
2280  This function is a more generic version of
2281  \ref o2scl_cblas::dnrm2 .
2282  */
2283  template<class vec_t>
2284  double vector_norm_double(size_t n, const vec_t &x) {
2285 
2286  double scale = 0.0;
2287  double ssq = 1.0;
2288 
2289  if (n <= 0) {
2290  return 0.0;
2291  } else if (n == 1) {
2292  return fabs(x[0]);
2293  }
2294 
2295  for (size_t i = 0; i < n; i++) {
2296  const double xx = x[i];
2297 
2298  if (xx != 0.0) {
2299  const double ax = fabs(xx);
2300 
2301  if (scale < ax) {
2302  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2303  scale = ax;
2304  } else {
2305  ssq += (ax / scale) * (ax / scale);
2306  }
2307  }
2308 
2309  }
2310 
2311  return scale * sqrt(ssq);
2312  }
2313 
2314  /** \brief Compute the norm of a vector of double precision numbers
2315  */
2316  template<class vec_t> double vector_norm_double(const vec_t &x) {
2317  return vector_norm_double<vec_t>(x.size(),x);
2318  }
2319  //@}
2320 
2321  /// \name Other vector and matrix functions
2322  //@{
2323  /** \brief Set the first N entries in a vector to a particular value
2324  */
2325  template<class vec_t, class data_t>
2326  void vector_set_all(size_t N, vec_t &src, data_t val) {
2327  for(size_t i=0;i<N;i++) {
2328  src[i]=val;
2329  }
2330  return;
2331  }
2332 
2333  /** \brief Set all entries in a vector to a particular value
2334  */
2335  template<class vec_t, class data_t>
2336  void vector_set_all(vec_t &src, data_t val) {
2337  o2scl::vector_set_all(src.size(),src,val);
2338  return;
2339  }
2340 
2341  /** \brief Set the first (M,N) entries in a matrix to a particular value
2342  */
2343  template<class mat_t, class data_t>
2344  void matrix_set_all(size_t M, size_t N, mat_t &src, data_t val) {
2345  for(size_t i=0;i<M;i++) {
2346  for(size_t j=0;j<N;j++) {
2347  src(i,j)=val;
2348  }
2349  }
2350  return;
2351  }
2352 
2353  /** \brief Set all entries in a matrix to a particular value
2354  */
2355  template<class mat_t, class data_t>
2356  void matrix_set_all(mat_t &src, data_t val) {
2357  o2scl::matrix_set_all(src.size1(),src.size2(),src,val);
2358  return;
2359  }
2360 
2361  /** \brief From a given vector, create a new vector by removing a
2362  specified element
2363 
2364  This funciton is used in \ref o2scl::interp_krige_optim::qual_fun() .
2365  */
2366  template<class vec_t, class vec2_t>
2367  void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest) {
2368  if (src.size()==0) {
2369  O2SCL_ERR("Empty source vector.",o2scl::exc_einval);
2370  }
2371  if (iout>=src.size()) {
2372  O2SCL_ERR("Requested element beyond end.",o2scl::exc_einval);
2373  }
2374  dest.resize(src.size()-1);
2375  size_t j=0;
2376  for(size_t i=0;i<src.size();i++) {
2377  if (i!=iout) {
2378  dest[j]=src[i];
2379  j++;
2380  }
2381  }
2382  return;
2383  }
2384 
2385  /** \brief "Rotate" a vector so that the kth element is now the beginning
2386 
2387  This is a generic template function which will work for
2388  any types \c data_t and \c vec_t for which
2389  - \c data_t has an <tt>operator=</tt>
2390  - <tt>vec_t::operator[]</tt> returns a reference
2391  to an object of type \c data_t
2392 
2393  This function is used, for example, in \ref o2scl::pinside.
2394 
2395  \note This function is not the same as a Givens rotation,
2396  which is typically referred to in BLAS routines as <tt>drot()</tt>.
2397  */
2398  template<class vec_t, class data_t>
2399  void vector_rotate(size_t n, vec_t &data, size_t k) {
2400 
2401  data_t *tmp=new data_t[n];
2402  for(size_t i=0;i<n;i++) {
2403  tmp[i]=data[(i+k)%n];
2404  }
2405  for(size_t i=0;i<n;i++) {
2406  data[i]=tmp[i];
2407  }
2408  delete[] tmp;
2409 
2410  return;
2411  }
2412 
2413  /** \brief Reverse the first \c n elements of a vector
2414 
2415  If \c n is zero, this function will silently do nothing.
2416  */
2417  template<class vec_t, class data_t>
2418  void vector_reverse(size_t n, vec_t &data) {
2419  data_t tmp;
2420 
2421  for(size_t i=0;i<n/2;i++) {
2422  tmp=data[n-1-i];
2423  data[n-1-i]=data[i];
2424  data[i]=tmp;
2425  }
2426  return;
2427  }
2428 
2429  /** \brief Reverse a vector
2430 
2431  If the <tt>size()</tt> method returns zero, this function will
2432  silently do nothing.
2433  */
2434  template<class vec_t, class data_t>
2435  void vector_reverse(vec_t &data) {
2436  data_t tmp;
2437  size_t n=data.size();
2438 
2439  for(size_t i=0;i<n/2;i++) {
2440  tmp=data[n-1-i];
2441  data[n-1-i]=data[i];
2442  data[i]=tmp;
2443  }
2444  return;
2445  }
2446 
2447  /** \brief Reverse the first n elements in a vector of double
2448  precision numbers
2449 
2450  If \c n is zero, this function will silently do nothing.
2451  */
2452  template<class vec_t>
2453  void vector_reverse_double(size_t n, vec_t &data) {
2454  double tmp;
2455 
2456  for(size_t i=0;i<n/2;i++) {
2457  tmp=data[n-1-i];
2458  data[n-1-i]=data[i];
2459  data[i]=tmp;
2460  }
2461  return;
2462  }
2463 
2464  /** \brief Reverse a vector of double precision numbers
2465 
2466  If the <tt>size()</tt> method returns zero, this function will
2467  silently do nothing.
2468  */
2469  template<class vec_t> void vector_reverse_double(vec_t &data) {
2470  double tmp;
2471  size_t n=data.size();
2472 
2473  for(size_t i=0;i<n/2;i++) {
2474  tmp=data[n-1-i];
2475  data[n-1-i]=data[i];
2476  data[i]=tmp;
2477  }
2478  return;
2479  }
2480 
2481  /** \brief Construct a row of a matrix
2482 
2483  This class template works with combinations of ublas
2484  <tt>matrix</tt> and <tt>matrix_row</tt> objects,
2485  <tt>arma::mat</tt> and <tt>arma::rowvec</tt>, and
2486  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2487 
2488  \note When calling this function with ublas objects, the
2489  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2490  otherwise some compilers will use argument dependent lookup and
2491  get (justifiably) confused with matrix_row in the ublas
2492  namespace.
2493 
2494  \note The template parameters must be explicitly specified
2495  when calling this template function.
2496  */
2497  template<class mat_t, class mat_row_t>
2498  mat_row_t matrix_row(mat_t &M, size_t row) {
2499  return mat_row_t(M,row);
2500  }
2501 
2502  /** \brief Generic object which represents a row of a matrix
2503 
2504  \note This class is experimental.
2505 
2506  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2507  to construct a row of a matrix object of type
2508  \code
2509  std::function<double &(size_t,size_t)>
2510  \endcode
2511  */
2512  template<class mat_t> class matrix_row_gen {
2513 
2514  protected:
2515 
2516  mat_t &m_;
2517 
2518  size_t row_;
2519 
2520  public:
2521 
2522  /// Create a row object from row \c row of matrix \c m
2523  matrix_row_gen(mat_t &m, size_t row) : m_(m), row_(row) {
2524  }
2525 
2526  /// Return a reference to the ith column of the selected row
2527  double &operator[](size_t i) {
2528  return m_(row_,i);
2529  }
2530 
2531  /// Return a const reference to the ith column of the selected row
2532  const double &operator[](size_t i) const {
2533  return m_(row_,i);
2534  }
2535  };
2536 
2537  /** \brief Generic object which represents a row of a matrix
2538 
2539  \note This class is experimental.
2540 
2541  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2542  to construct a row of a matrix object of type
2543  \code
2544  std::function<double &(size_t,size_t)>
2545  \endcode
2546  */
2547  template<class mat_t> class const_matrix_row_gen {
2548 
2549  protected:
2550 
2551  const mat_t &m_;
2552 
2553  size_t row_;
2554 
2555  public:
2556 
2557  /// Create a row object from row \c row of matrix \c m
2558  const_matrix_row_gen(const mat_t &m, size_t row) : m_(m), row_(row) {
2559  }
2560 
2561  /// Return a const reference to the ith column of the selected row
2562  const double &operator[](size_t i) const {
2563  return m_(row_,i);
2564  }
2565  };
2566 
2567  /** \brief A simple matrix view object
2568  */
2569  class matrix_view {
2570 
2571  public:
2572 
2573  /** \brief Return a reference to the element at row \c row
2574  and column \c col
2575  */
2576  const double &operator()(size_t row, size_t col) const;
2577  /** \brief Return the number of rows
2578  */
2579  size_t size1();
2580  /** \brief Return the number of columns
2581  */
2582  size_t size2();
2583 
2584  };
2585 
2586  /** \brief Construct a column of a matrix
2587 
2588  This class template works with combinations of ublas
2589  <tt>matrix</tt> and <tt>matrix_column</tt> objects,
2590  <tt>arma::mat</tt> and <tt>arma::colvec</tt>, and
2591  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2592 
2593  \note When calling this function with ublas objects, the
2594  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2595  otherwise some compilers will use argument dependent lookup and
2596  get (justifiably) confused with matrix_column in the ublas
2597  namespace.
2598 
2599  \note The template parameters must be explicitly specified
2600  when calling this template function.
2601  */
2602  template<class mat_t, class mat_column_t>
2603  mat_column_t matrix_column(mat_t &M, size_t column) {
2604  return mat_column_t(M,column);
2605  }
2606 
2607  /** \brief Generic object which represents a column of a matrix
2608 
2609  \note This class is experimental.
2610 
2611  The only requirement on the type <tt>mat_t</tt> is that
2612  it must have an operator(size_t,size_t) method which
2613  accesses elements in the matrix.
2614 
2615  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2616  to construct a row of a matrix object of type
2617  \code
2618  std::function<double &(size_t,size_t)>
2619  \endcode
2620  */
2621  template<class mat_t> class matrix_column_gen {
2622  protected:
2623  mat_t &m_;
2624  size_t column_;
2625  public:
2626  matrix_column_gen(mat_t &m, size_t column) : m_(m), column_(column) {
2627  }
2628  double &operator[](size_t i) {
2629  return m_(i,column_);
2630  }
2631  const double &operator[](size_t i) const {
2632  return m_(i,column_);
2633  }
2634  };
2635 
2636  /** \brief Generic object which represents a column of a const matrix
2637 
2638  \note This class is experimental.
2639 
2640  The only requirement on the type <tt>mat_t</tt> is that
2641  it must have an operator(size_t,size_t) method which
2642  accesses elements in the matrix.
2643 
2644  This class is used in one of
2645  the \ref o2scl::prob_dens_mdim_gaussian constructors.
2646  */
2647  template<class mat_t> class const_matrix_column_gen {
2648  protected:
2649  const mat_t &m_;
2650  size_t column_;
2651  public:
2652  const_matrix_column_gen(const mat_t &m, size_t column) :
2653  m_(m), column_(column) {
2654  }
2655  const double &operator[](size_t i) const {
2656  return m_(i,column_);
2657  }
2658  };
2659 
2660  /** \brief Output the first \c n elements of a vector to a stream,
2661  \c os
2662 
2663  No trailing space is output after the last element, and an
2664  endline is output only if \c endline is set to \c true. If the
2665  parameter \c n is zero, this function silently does nothing.
2666 
2667  This works with any class <tt>vec_t</tt> which has an operator[]
2668  which returns either the value of or a reference to the ith
2669  element and the element type has its own output operator which
2670  has been defined.
2671  */
2672  template<class vec_t>
2673  void vector_out(std::ostream &os, size_t n, const vec_t &v,
2674  bool endline=false) {
2675 
2676  // This next line is important since n-1 is not well-defined if n=0
2677  if (n==0) return;
2678 
2679  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
2680  os << v[n-1];
2681  if (endline) os << std::endl;
2682  return;
2683  }
2684 
2685  /** \brief Output a vector to a stream
2686 
2687  No trailing space is output after the last element, and an
2688  endline is output only if \c endline is set to \c true. If the
2689  parameter \c n is zero, this function silently does nothing.
2690 
2691  This works with any class <tt>vec_t</tt> which has an operator[]
2692  which returns either the value of or a reference to the ith
2693  element and the element type has its own output operator which
2694  has been defined.
2695  */
2696  template<class vec_t>
2697  void vector_out(std::ostream &os, const vec_t &v, bool endline=false) {
2698 
2699  size_t n=v.size();
2700 
2701  // This next line is important since n-1 is not well-defined if n=0
2702  if (n==0) return;
2703 
2704  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
2705  os << v[n-1];
2706  if (endline) os << std::endl;
2707  return;
2708  }
2709 
2710  /** \brief Fill a vector with a specified grid
2711  */
2712  template<class vec_t, class data_t>
2713  void vector_grid(uniform_grid<data_t> g, vec_t &v) {
2714  g.template vector<vec_t>(v);
2715  return;
2716  }
2717 
2718  /// Set a matrix to unity on the diagonal and zero otherwise
2719  template<class mat_t>
2720  void matrix_set_identity(size_t M, size_t N, mat_t &m) {
2721  for(size_t i=0;i<M;i++) {
2722  for(size_t j=0;j<N;j++) {
2723  if (i==j) m(i,j)=1.0;
2724  else m(i,j)=0.0;
2725  }
2726  }
2727  return;
2728  }
2729 
2730  /// Set a matrix to unity on the diagonal and zero otherwise
2731  template<class mat_t>
2732  void matrix_set_identity(mat_t &m) {
2733  matrix_set_identity(m.size1(),m.size2(),m);
2734  return;
2735  }
2736  //@}
2737 
2738  /// \name Vector range classes and functions
2739  //@{
2740  /** \brief Vector range function for pointers
2741 
2742  \note In this case, the return type is the same as the
2743  type of the first parameter.
2744  */
2745  template<class dat_t> dat_t *vector_range
2746  (dat_t *v, size_t start, size_t last) {
2747  return v+start;
2748  }
2749 
2750  /** \brief Vector range function for const pointers
2751 
2752  \note In this case, the return type is the same as the
2753  type of the first parameter.
2754  */
2755  template<class dat_t> const dat_t *const_vector_range
2756  (const dat_t *v, size_t start, size_t last) {
2757  return v+start;
2758  }
2759 
2760  /** \brief Vector range function template for ublas vectors
2761 
2762  The element with index \c start in the original vector
2763  will become the first argument in the new vector, and
2764  the new vector will have size <tt>last-start</tt> .
2765 
2766  \note In this case, the return type is not the same as the
2767  type of the first parameter.
2768  */
2769  template<class dat_t> boost::numeric::ublas::vector_range
2772  size_t start, size_t last) {
2775  (v,boost::numeric::ublas::range(start,last));
2776  }
2777 
2778  /** \brief Const vector range function template for ublas vectors
2779 
2780  The element with index \c start in the original vector
2781  will become the first argument in the new vector, and
2782  the new vector will have size <tt>last-start</tt> .
2783 
2784  \note In this case, the return type is not the same as the
2785  type of the first parameter.
2786  */
2787  template<class dat_t> const boost::numeric::ublas::vector_range
2790  size_t start, size_t last) {
2793  (v,boost::numeric::ublas::range(start,last));
2794  }
2795 
2796  /** \brief Const vector range function template for const ublas
2797  vectors
2798 
2799  The element with index \c start in the original vector
2800  will become the first argument in the new vector, and
2801  the new vector will have size <tt>last-start</tt> .
2802 
2803  \note In this case, the return type is not the same as the
2804  type of the first parameter.
2805  */
2806  template<class dat_t> const boost::numeric::ublas::vector_range
2809  size_t start, size_t last) {
2812  (v,boost::numeric::ublas::range(start,last));
2813  }
2814 
2815  /** \brief Vector range function template for ublas vector
2816  ranges of ublas vectors
2817 
2818  The element with index \c start in the original vector
2819  will become the first argument in the new vector, and
2820  the new vector will have size <tt>last-start</tt> .
2821 
2822  \note In this case, the return type is not the same as the
2823  type of the first parameter.
2824  */
2825  template<class dat_t>
2829  vector_range
2832  size_t start, size_t last) {
2836  (v,boost::numeric::ublas::range(start,last));
2837  }
2838 
2839  /** \brief Const vector range function template for ublas vector
2840  ranges of ublas vectors
2841 
2842  The element with index \c start in the original vector
2843  will become the first argument in the new vector, and
2844  the new vector will have size <tt>last-start</tt> .
2845 
2846  \note In this case, the return type is not the same as the
2847  type of the first parameter.
2848  */
2849  template<class dat_t>
2856  size_t start, size_t last) {
2860  (v,boost::numeric::ublas::range(start,last));
2861  }
2862 
2863  /** \brief Const vector range function template for const
2864  ublas vector ranges of ublas vectors
2865 
2866  The element with index \c start in the original vector
2867  will become the first argument in the new vector, and
2868  the new vector will have size <tt>last-start</tt> .
2869 
2870  \note In this case, the return type is not the same as the
2871  type of the first parameter.
2872  */
2873  template<class dat_t>
2880  size_t start, size_t last) {
2884  (v,boost::numeric::ublas::range(start,last));
2885  }
2886 
2887  /** \brief Const vector range function template for const
2888  ublas vector ranges of const ublas vectors
2889 
2890  The element with index \c start in the original vector
2891  will become the first argument in the new vector, and
2892  the new vector will have size <tt>last-start</tt> .
2893 
2894  \note In this case, the return type is not the same as the
2895  type of the first parameter.
2896  */
2897  template<class dat_t>
2904  size_t start, size_t last) {
2908  (v,boost::numeric::ublas::range(start,last));
2909  }
2910 
2911  // Forward definition for friendship
2912  template<class vec_t> class const_vector_range_gen;
2913 
2914  /** \brief Experimental vector range object
2915  */
2916  template<class vec_t> class vector_range_gen {
2917 
2918  protected:
2919 
2920  friend class const_vector_range_gen<vec_t>;
2921 
2922  /// A reference to the original vector
2923  vec_t &v_;
2924 
2925  /// The index offset
2926  size_t start_;
2927 
2928  /// The end() iterator
2929  size_t last_;
2930 
2931  public:
2932 
2933  /// Create an object starting with index \c start in vector \c v
2934  vector_range_gen(vec_t &v, size_t start, size_t last) : v_(v),
2935  start_(start), last_(last) {
2936 #if !O2SCL_NO_RANGE_CHECK
2937  if (last<start) {
2938  O2SCL_ERR2("End before beginning in vector_range_gen::",
2939  "vector_range_gen(vec_t,size_t,size_t)",
2941  }
2942 #endif
2943  }
2944 
2945  /// Create an object from a previously constructed range object
2946  vector_range_gen(const vector_range_gen &v2, size_t start,
2947  size_t last) : v_(v2.v_),
2948  start_(start+v2.start_), last_(last+v2.start_) {
2949 #if !O2SCL_NO_RANGE_CHECK
2950  if (last<start) {
2951  O2SCL_ERR2("End before beginning in vector_range_gen::",
2952  "vector_range_gen(vector_range_gen,size_t,size_t)",
2954  }
2955  if (last>v2.last_) {
2956  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2957  "vector_range_gen(vector_range_gen,size_t,size_t)",
2959  }
2960 #endif
2961  }
2962 
2963  /// Return the vector size
2964  size_t size() const {
2965  return last_-start_;
2966  }
2967 
2968  /// Return a reference ith element
2969  double &operator[](size_t i) {
2970 #if !O2SCL_NO_RANGE_CHECK
2971  if (i+start_>=last_) {
2972  O2SCL_ERR("Index out of range in vector_range_gen::operator[].",
2974  }
2975 #endif
2976  return v_[i+start_];
2977  }
2978 
2979  /// Return a const reference ith element
2980  const double &operator[](size_t i) const {
2981 #if !O2SCL_NO_RANGE_CHECK
2982  if (i+start_>=last_) {
2983  O2SCL_ERR2("Index out of range in ",
2984  "vector_range_gen::operator[] const.",o2scl::exc_einval);
2985  }
2986 #endif
2987  return v_[i+start_];
2988  }
2989  };
2990 
2991  /** \brief Experimental const vector range object
2992  */
2993  template<class vec_t> class const_vector_range_gen {
2994 
2995  protected:
2996 
2997  /// A reference to the original vector
2998  const vec_t &v_;
2999 
3000  /// The index offset
3001  size_t start_;
3002 
3003  /// The end() iterator
3004  size_t last_;
3005 
3006  public:
3007 
3008  /// Create an object starting with index \c start in vector \c v
3009  const_vector_range_gen(const vec_t &v, size_t start, size_t last) : v_(v),
3010  start_(start), last_(last) {
3011 #if !O2SCL_NO_RANGE_CHECK
3012  if (last<start) {
3013  O2SCL_ERR2("End before beginning in vector_range_gen::",
3014  "vector_range_gen(vec_t,size_t,size_t)",
3016  }
3017 #endif
3018  }
3019 
3020  /// Create an object from a previously constructed range object
3022  size_t last) : v_(v2.v_),
3023  start_(start+v2.start_), last_(last+v2.start_) {
3024 #if !O2SCL_NO_RANGE_CHECK
3025  if (last<start) {
3026  O2SCL_ERR2("End before beginning in vector_range_gen::",
3027  "vector_range_gen(vector_range_gen,size_t,size_t)",
3029  }
3030  if (last>v2.last_) {
3031  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
3032  "vector_range_gen(vector_range_gen,size_t,size_t)",
3034  }
3035 #endif
3036  }
3037 
3038  /// Create an object from a previously constructed range object
3040  size_t last) : v_(v2.v_),
3041  start_(start+v2.start_), last_(last+v2.start_) {
3042 #if !O2SCL_NO_RANGE_CHECK
3043  if (last<start) {
3044  O2SCL_ERR2("End before beginning in vector_range_gen::",
3045  "vector_range_gen(vector_range_gen,size_t,size_t)",
3047  }
3048  if (last>v2.last_) {
3049  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
3050  "vector_range_gen(vector_range_gen,size_t,size_t)",
3052  }
3053 #endif
3054  }
3055 
3056  /// Return the vector size
3057  size_t size() const {
3058  return last_-start_;
3059  }
3060 
3061  /// Return a const reference ith element
3062  const double &operator[](size_t i) const {
3063 #if !O2SCL_NO_RANGE_CHECK
3064  if (i+start_>=last_) {
3065  O2SCL_ERR2("Index out of range in ",
3066  "vector_range_gen::operator[] const.",o2scl::exc_einval);
3067  }
3068 #endif
3069  return v_[i+start_];
3070  }
3071  };
3072 
3073  /** \brief Create a \ref o2scl::vector_range_gen object
3074  from a <tt>std::vector</tt>
3075  */
3076  template<class data_t> vector_range_gen<std::vector<data_t> >
3077  vector_range(std::vector<data_t> &v, size_t start, size_t last) {
3078  return vector_range_gen<std::vector<data_t> >(v,start,last);
3079  }
3080 
3081  /** \brief Create a \ref o2scl::vector_range_gen object
3082  from a <tt>std::vector</tt>
3083  */
3084  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
3085  const_vector_range(const std::vector<data_t> &v, size_t start,
3086  size_t last) {
3087  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
3088  }
3089 
3090  /** \brief Create a \ref o2scl::vector_range_gen object
3091  from a <tt>std::vector</tt>
3092  */
3093  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
3094  const_vector_range(std::vector<data_t> &v, size_t start,
3095  size_t last) {
3096  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
3097  }
3098 
3099  /** \brief Recursively create a \ref o2scl::vector_range_gen object
3100  from a vector range
3101  */
3102  template<class vec_t> vector_range_gen<vec_t>
3103  vector_range(vector_range_gen<vec_t> &v, size_t start, size_t last) {
3104  return vector_range_gen<vec_t>(v,start,last);
3105  }
3106 
3107  /** \brief Recursively create a const \ref o2scl::vector_range_gen
3108  object from a vector range
3109  */
3110  template<class vec_t> const const_vector_range_gen<vec_t>
3112  size_t start, size_t last) {
3113  return const_vector_range_gen<vec_t>(v,start,last);
3114  }
3115 
3116  /** \brief Recursively create a const \ref o2scl::vector_range_gen
3117  object from a const vector range
3118  */
3119  template<class vec_t> const const_vector_range_gen<vec_t>
3121  size_t start, size_t last) {
3122  return const_vector_range_gen<vec_t>(v,start,last);
3123  }
3124 
3125  /** \brief Recursively create a const \ref o2scl::vector_range_gen
3126  object from a const vector range
3127  */
3128  template<class vec_t> const const_vector_range_gen<vec_t>
3130  size_t start, size_t last) {
3131  return const_vector_range_gen<vec_t>(v,start,last);
3132  }
3133 
3134  /** \brief Vector range function template for <tt>std::vector</tt>
3135 
3136  The element with index \c start in the original vector
3137  will become the first argument in the new vector, and
3138  the new vector will have size <tt>last-start</tt> .
3139 
3140  \note In this case, the return type is the same as the
3141  type of the first parameter.
3142  \note Unlike the ublas and pointer cases, this forces
3143  a copy.
3144  */
3145  template<class dat_t> std::vector<dat_t>
3146  vector_range_copy(const std::vector<dat_t> &v, size_t start, size_t last) {
3147  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3148  }
3149 
3150  /** \brief Const vector range function template for <tt>std::vector</tt>
3151 
3152  The element with index \c start in the original vector
3153  will become the first argument in the new vector, and
3154  the new vector will have size <tt>last-start</tt> .
3155 
3156  \note In this case, the return type is the same as the
3157  type of the first parameter.
3158  \note Unlike the ublas and pointer cases, this forces
3159  a copy.
3160  */
3161  template<class dat_t> const std::vector<dat_t>
3162  vector_range_copy(const std::vector<dat_t> &v, size_t start,
3163  size_t last) {
3164  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3165  }
3166 
3167  //@}
3168 
3169 }
3170 
3171 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
3172 
3173 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
3174 #include <armadillo>
3175 namespace o2scl {
3176 
3177  /// \name Armadillo specializations
3178  //@{
3179  /// Armadillo version of \ref matrix_max()
3180  double matrix_max(const arma::mat &data);
3181 
3182  /// Armadillo version of \ref matrix_min()
3183  double matrix_min(const arma::mat &data);
3184 
3185  /// Armadillo version of \ref matrix_row()
3186  template<> arma::subview_row<double>
3187  matrix_row<arma::mat,arma::subview_row<double> >
3188  (arma::mat &M, size_t row);
3189 
3190  /// Armadillo version of \ref matrix_column()
3191  template<> arma::subview_col<double>
3192  matrix_column<arma::mat,arma::subview_col<double> >
3193  (arma::mat &M, size_t column);
3194  //@}
3195 
3196 }
3197 
3198 #endif
3199 
3200 #if defined (O2SCL_EIGEN) || defined (DOXYGEN)
3201 #include <eigen3/Eigen/Dense>
3202 
3203 namespace o2scl {
3204 
3205  /// \name Eigen specializations
3206  //@{
3207  /// Eigen version of \ref matrix_max()
3208  double matrix_max(const Eigen::MatrixXd &data);
3209 
3210  /// Eigen version of \ref matrix_min()
3211  double matrix_min(const Eigen::MatrixXd &data);
3212 
3213  /// Eigen version of \ref matrix_row()
3214  template<> Eigen::MatrixXd::RowXpr
3215  matrix_row<Eigen::MatrixXd,Eigen::MatrixXd::RowXpr>
3216  (Eigen::MatrixXd &M, size_t row);
3217 
3218  /// Eigen version of \ref matrix_column()
3219  template<> Eigen::MatrixXd::ColXpr
3220  matrix_column<Eigen::MatrixXd,Eigen::MatrixXd::ColXpr>
3221  (Eigen::MatrixXd &M, size_t column);
3222  //@}
3223 
3224 }
3225 
3226 #endif
3227 
3228 #else
3229 
3230 #include <o2scl/vector_special.h>
3231 
3232 // End of "#if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)"
3233 #endif
3234 
3235 #ifdef DOXYGEN
3236 /** \brief Placeholder documentation of some related Boost objects
3237  */
3238 namespace boost {
3239  /** \brief Documentation of Boost::numeric objects
3240  */
3241  namespace numeric {
3242  /** \brief Documentation of uBlas objects
3243  */
3244  namespace ublas {
3245  /** \brief The default vector type from uBlas
3246 
3247  The uBlas types aren't documented here, but the full documentation
3248  is available at
3249  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
3250 
3251  Internally in \o2, this is often typedef'd using
3252  \code
3253  typedef boost::numeric::ublas::vector<double> ubvector;
3254  typedef boost::numeric::ublas::vector<size_t> ubvector_size_t;
3255  typedef boost::numeric::ublas::vector<int> ubvector_int;
3256  \endcode
3257 
3258  This is documented in \ref vector.h .
3259  */
3260  template<class T, class A> class vector {
3261  };
3262  /** \brief The default matrix type from uBlas
3263 
3264  The uBlas types aren't documented here, but the full documentation
3265  is available at
3266  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
3267 
3268  Internally in \o2, this is often typedef'd using
3269  \code
3270  typedef boost::numeric::ublas::matrix<double> ubmatrix;
3271  typedef boost::numeric::ublas::matrix<size_t> ubmatrix_size_t;
3272  typedef boost::numeric::ublas::matrix<int> ubmatrix_int;
3273  \endcode
3274 
3275  This is documented in \ref vector.h .
3276  */
3277  template<class T, class F, class A> class matrix {
3278  };
3279  }
3280  }
3281 }
3282 // End of "#ifdef DOXYGEN"
3283 #endif
3284 
3285 // End of "#ifndef O2SCL_VECTOR_H"
3286 #endif
data_t vector_min_quad(size_t n, const vec_t &data)
Minimum of vector by quadratic fit.
Definition: vector.h:1390
size_t vector_lookup(size_t n, const vec_t &x, double x0)
Lookup the value x0 in the first n elements of vector x.
Definition: vector.h:1783
size_t vector_min_index(size_t n, const vec_t &data)
Compute the index which holds the minimum of the first n elements of a vector.
Definition: vector.h:1235
Experimental vector range object.
Definition: vector.h:2916
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
Definition: vector.h:1199
data_t matrix_sum(size_t m, size_t n, const mat_t &data)
Compute the sum of matrix elements.
Definition: vector.h:1748
void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest)
From a given vector, create a new vector by removing a specified element.
Definition: vector.h:2367
data_t vector_max_quad_loc(size_t n, const vec_t &x, const vec_t &y)
Location of vector maximum by quadratic fit.
Definition: vector.h:1380
data_t vector_min_quad_loc(size_t n, const vec_t &x, const vec_t &y)
Location of vector minimum by quadratic fit.
Definition: vector.h:1415
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
void matrix_swap_rows(size_t N, mat_t &m, size_t i1, size_t i2)
Swap the first N columns of two rows in a matrix.
Definition: vector.h:626
void matrix_make_upper(mat_t &src)
Make a matrix upper triangular by setting the lower triangular entries to zero.
Definition: vector.h:349
bool vector_is_finite(size_t n, vec_t &data)
Test if the first n elements of a vector are finite.
Definition: vector.h:2155
Placeholder documentation of some related Boost objects.
Definition: vector.h:3238
size_t vector_bsearch_dec(const data_t x0, const vec_t &x, size_t lo, size_t hi)
Binary search a part of an decreasing vector for x0.
Definition: vector.h:1954
void matrix_max_index(size_t m, size_t n, const mat_t &data, size_t &i_max, size_t &j_max, data_t &max)
Compute the maximum of a matrix and return the indices of the maximum element.
Definition: vector.h:1497
const_vector_range_gen(const vec_t &v, size_t start, size_t last)
Create an object starting with index start in vector v.
Definition: vector.h:3009
void vector_set_all(size_t N, vec_t &src, data_t val)
Set the first N entries in a vector to a particular value.
Definition: vector.h:2326
Generic object which represents a column of a matrix.
Definition: vector.h:2621
bool matrix_is_upper(mat_t &src)
Simple test that a matrix is upper triangular.
Definition: vector.h:320
bool matrix_is_lower(mat_t &src)
Simple test that a matrix is lower triangular.
Definition: vector.h:306
double & operator[](size_t i)
Return a reference ith element.
Definition: vector.h:2969
size_t vector_bsearch_inc(const data_t x0, const vec_t &x, size_t lo, size_t hi)
Binary search a part of an increasing vector for x0.
Definition: vector.h:1909
A simple convenience wrapper for GSL vector objects.
Definition: vector.h:71
void matrix_set_all(size_t M, size_t N, mat_t &src, data_t val)
Set the first (M,N) entries in a matrix to a particular value.
Definition: vector.h:2344
Generic object which represents a column of a const matrix.
Definition: vector.h:2647
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2562
invalid argument supplied by user
Definition: err_hnd.h:59
Generic object which represents a row of a matrix.
Definition: vector.h:2512
A simple matrix view object.
Definition: vector.h:2569
vector_range_gen< vec_t > vector_range(vector_range_gen< vec_t > &v, size_t start, size_t last)
Recursively create a o2scl::vector_range_gen object from a vector range.
Definition: vector.h:3103
void vector_reverse(size_t n, vec_t &data)
Reverse the first n elements of a vector.
Definition: vector.h:2418
size_t last_
The end() iterator.
Definition: vector.h:2929
void vector_rotate(size_t n, vec_t &data, size_t k)
"Rotate" a vector so that the kth element is now the beginning
Definition: vector.h:2399
void vector_swap_double(size_t N, vec_t &v1, vec2_t &v2)
Swap of of the first N elements of two double-precision vectors.
Definition: vector.h:493
size_t last_
The end() iterator.
Definition: vector.h:3004
double matrix_min_value_double(const mat_t &data)
Compute the minimum of a matrix.
Definition: vector.h:1600
The default matrix type from uBlas.
Definition: vector.h:3277
void vector_min(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the minimum of the first n elements of a vector.
Definition: vector.h:1254
Generic object which represents a row of a matrix.
Definition: vector.h:2547
const_matrix_row_gen(const mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2558
double vector_norm_double(size_t n, const vec_t &x)
Compute the norm of the first n entries of a vector of double precision numbers.
Definition: vector.h:2284
generic failure
Definition: err_hnd.h:61
void matrix_make_lower(mat_t &src)
Make a matrix lower triangular by setting the upper triangular entries to zero.
Definition: vector.h:335
double matrix_max_value_double(const mat_t &data)
Compute the maximum of a matrix.
Definition: vector.h:1473
data_t vector_sum(size_t n, vec_t &data)
Compute the sum of the first n elements of a vector.
Definition: vector.h:2183
void vector_swap(size_t N, vec_t &v1, vec2_t &v2)
Swap the first N elements of two vectors.
Definition: vector.h:423
void vector_largest(size_t n, vec_t &data, size_t k, vec_t &largest)
Find the k largest entries of the first n elements of a vector.
Definition: vector.h:1061
std::vector< dat_t > vector_range_copy(const std::vector< dat_t > &v, size_t start, size_t last)
Vector range function template for std::vector
Definition: vector.h:3146
const_vector_range_gen(const vector_range_gen< vec_t > &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
Definition: vector.h:3039
matrix_row_gen(mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2523
vec_t & v_
A reference to the original vector.
Definition: vector.h:2923
void matrix_swap_double(size_t M, size_t N, mat_t &m1, mat2_t &m2)
Swap of the first elements in two double-precision matrices.
Definition: vector.h:564
void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order)
Create a permutation which sorts the first n elements of a vector (in increasing order) ...
Definition: vector.h:803
void matrix_swap_cols(size_t M, mat_t &m, size_t j1, size_t j2)
Swap the first M rows of two columns in a matrix.
Definition: vector.h:598
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
void matrix_copy(mat_t &src, mat2_t &dest)
Simple matrix copy.
Definition: vector.h:179
size_t size() const
Return the vector size.
Definition: vector.h:3057
void matrix_minmax_index(size_t n, size_t m, const mat_t &data, size_t &i_min, size_t &j_min, data_t &min, size_t &i_max, size_t &j_max, data_t &max)
Compute the minimum and maximum of a matrix and return their locations.
Definition: vector.h:1715
data_t vector_max_value(size_t n, const vec_t &data)
Compute the maximum of the first n elements of a vector.
Definition: vector.h:1124
const vec_t & v_
A reference to the original vector.
Definition: vector.h:2998
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:2746
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
const_vector_range_gen(const const_vector_range_gen &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
Definition: vector.h:3021
void vector_minmax_value(size_t n, vec_t &data, data_t &min, data_t &max)
Compute the minimum and maximum of the first n elements of a vector.
Definition: vector.h:1275
void matrix_set_identity(size_t M, size_t N, mat_t &m)
Set a matrix to unity on the diagonal and zero otherwise.
Definition: vector.h:2720
data_t vector_max_quad(size_t n, const vec_t &data)
Maximum of vector by quadratic fit.
Definition: vector.h:1355
vector_range_gen(const vector_range_gen &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
Definition: vector.h:2946
mat_column_t matrix_column(mat_t &M, size_t column)
Construct a column of a matrix.
Definition: vector.h:2603
size_t start_
The index offset.
Definition: vector.h:2926
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
mat_row_t matrix_row(mat_t &M, size_t row)
Construct a row of a matrix.
Definition: vector.h:2498
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
Definition: vector.h:2673
void vector_smallest(size_t n, vec_t &data, size_t k, vec_t &smallest)
Find the k smallest entries of the first n elements of a vector.
Definition: vector.h:923
void vector_grid(uniform_grid< data_t > g, vec_t &v)
Fill a vector with a specified grid.
Definition: vector.h:2713
void vector_minmax(size_t n, vec_t &data, size_t &ix_min, data_t &min, size_t &ix_max, data_t &max)
Compute the minimum and maximum of the first n elements of a vector.
Definition: vector.h:1325
data_t vector_norm(size_t n, const vec_t &x)
Compute the norm of the first n entries of a vector of floating-point (single or double precision) nu...
Definition: vector.h:2240
A simple convenience wrapper for GSL matrix objects.
Definition: vector.h:92
size_t vector_max_index(size_t n, const vec_t &data)
Compute the index which holds the maximum of the first n elements of a vector.
Definition: vector.h:1160
void vector_max(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the maximum of the first n elements of a vector.
Definition: vector.h:1179
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2532
void vector_minmax_index(size_t n, vec_t &data, size_t &ix_min, size_t &ix_max)
Compute the minimum and maximum of the first n elements of a vector.
Definition: vector.h:1298
void vector_smallest_index(size_t n, const vec_t &data, size_t k, vec_size_t &index)
Find the indexes of the k smallest entries among the first n entries of a vector. ...
Definition: vector.h:996
double matrix_max(const arma::mat &data)
Armadillo version of matrix_max()
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:2980
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:3062
double vector_sum_double(size_t n, vec_t &data)
Compute the sum of the first n elements of a vector of double-precision numbers.
Definition: vector.h:2211
vector_range_gen(vec_t &v, size_t start, size_t last)
Create an object starting with index start in vector v.
Definition: vector.h:2934
A class representing a uniform linear or logarithmic grid.
Definition: uniform_grid.h:38
void matrix_minmax(size_t n, size_t m, const mat_t &data, data_t &min, data_t &max)
Compute the minimum and maximum of a matrix.
Definition: vector.h:1682
Experimental const vector range object.
Definition: vector.h:2912
size_t vector_bsearch(const data_t x0, const vec_t &x, size_t lo, size_t hi)
Binary search a part of a monotonic vector for x0.
Definition: vector.h:1980
void matrix_transpose(mat_t &src, mat2_t &dest)
Simple transpose.
Definition: vector.h:225
void sort_index_downheap(size_t N, const vec_t &data, vec_size_t &order, size_t k)
Provide a downheap() function for vector_sort_index()
Definition: vector.h:737
void vector_reverse_double(size_t n, vec_t &data)
Reverse the first n elements in a vector of double precision numbers.
Definition: vector.h:2453
int vector_is_monotonic(size_t n, vec_t &data)
Test if the first n elements of a vector are monotonic and increasing or decreasing.
Definition: vector.h:2022
int vector_is_strictly_monotonic(size_t n, vec_t &data)
Test if the first n elements of a vector are strictly monotonic and determine if they are increasing ...
Definition: vector.h:2104
void matrix_swap_cols_double(size_t M, mat_t &m, size_t j1, size_t j2)
Swap the first M rows of two columns in a double-precision matrix.
Definition: vector.h:615
void matrix_min_index(size_t n, size_t m, const mat_t &data, size_t &i_min, size_t &j_min, data_t &min)
Compute the minimum of a matrix and return the indices of the minimum element.
Definition: vector.h:1625
double matrix_min(const arma::mat &data)
Armadillo version of matrix_min()
void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2)
Swap of the first elements in two matrices.
Definition: vector.h:545
void matrix_swap_rows_double(size_t N, mat_t &m, size_t i1, size_t i2)
Swap the first N columns of two rows in a double-precision matrix.
Definition: vector.h:643
data_t matrix_max_value(size_t m, const size_t n, const mat_t &data)
Compute the maximum of the lower-left part of a matrix.
Definition: vector.h:1428
size_t start_
The index offset.
Definition: vector.h:3001
std::string itos(int x)
Convert an integer to a string.
std::string szttos(size_t x)
Convert a size_t to a string.
const dat_t * const_vector_range(const dat_t *v, size_t start, size_t last)
Vector range function for const pointers.
Definition: vector.h:2756
void matrix_lookup(size_t m, size_t n, const mat_t &A, double x0, size_t &i, size_t &j)
Lookup an element in the first $(m,n)$ entries in a matrix.
Definition: vector.h:1830
double & operator[](size_t i)
Return a reference to the ith column of the selected row.
Definition: vector.h:2527
data_t matrix_min_value(size_t m, size_t n, const mat_t &data)
Compute the minimum of a matrix.
Definition: vector.h:1554
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
Definition: vector.h:898
The default vector type from uBlas.
Definition: vector.h:3260
size_t size() const
Return the vector size.
Definition: vector.h:2964
void vector_sort(size_t n, vec_t &data)
Sort a vector (in increasing order)
Definition: vector.h:708
void sort_downheap(vec_t &data, size_t n, size_t k)
Provide a downheap() function for vector_sort()
Definition: vector.h:653

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).