Line data Source code
1 : /** \file vectorUtils.hpp
2 : *
3 : * \brief Header for the std::vector utilities
4 : *
5 : * \ingroup gen_math_files
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015-2026 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef math_vectorUtils_hpp
28 : #define math_vectorUtils_hpp
29 :
30 : #include <vector>
31 : #include <algorithm>
32 : #include <functional>
33 : #include <numeric>
34 :
35 : #include "func/gaussian.hpp"
36 :
37 : namespace mx
38 : {
39 : namespace math
40 : {
41 :
42 : /** \ingroup vectorutils
43 : *@{
44 : */
45 :
46 : /// Fill in a vector with a regularly spaced scale
47 : /** Fills in the vector with a 0....N-1 scale. The spacing of the points
48 : * can be changed with the scale parameter, and the starting point can be
49 : * changed with the offset.
50 : *
51 : * Example:
52 : \code
53 : std::vector vec;
54 :
55 : mx::vectorScale(vec, 1000, 0.001, 0.001); // Produces a vector with values 0.001,0.002,.... 1.000
56 : \endcode
57 : *
58 : * \tparam vectorT is a std::vector type.
59 : */
60 : template <typename vectorT>
61 1 : void vectorScale(
62 : vectorT &vec, /**< [out] the vector to fill in, can be
63 : pre-allocated or not*/
64 : size_t N = 0, /**< [in] [optional] if specified > 0, then vec is resize()-ed. Default
65 : is 0, and vec is not resize()-ed.*/
66 : typename vectorT::value_type scale = 0, /**< [in] [optional] if specified !=0, then the points are
67 : spaced by this value. Default spacing is 1.*/
68 : typename vectorT::value_type offset = 0 /**< [in] [optional] if specified !=0, then the starting
69 : point of the scale is this value.*/
70 : )
71 : {
72 1 : if( scale == 0 )
73 0 : scale = 1.0;
74 :
75 1 : if( N > 0 )
76 1 : vec.resize( N );
77 :
78 2 : for( int i = 0; i < vec.size(); ++i )
79 1 : vec[i] = i * scale + offset;
80 1 : }
81 :
82 : /// Return the indices of the vector in sorted order, without altering the vector itself
83 : /** Example:
84 : \code
85 : std::vector<double> x;
86 : // -> fill in x with values
87 :
88 : std::vector<size_t> idx;
89 : idx = mx::sortOrder(x);
90 :
91 : //Print x to stdout in sorted order
92 : for(int i=0; i< x.size(); ++i) std::cout << x[idx[i]] << "\n";
93 : \endcode
94 :
95 : * \tparam memberT is the member type of the vector. Must have < comparison defined.
96 : */
97 : template <typename memberT>
98 0 : std::vector<size_t> vectorSortOrder( std::vector<memberT> const &values /**< [in] the vector to sort */ )
99 : {
100 0 : std::vector<size_t> indices( values.size() );
101 :
102 0 : std::iota( begin( indices ), end( indices ), static_cast<size_t>( 0 ) );
103 :
104 0 : std::sort( begin( indices ), end( indices ), [&]( size_t a, size_t b ) { return values[a] < values[b]; } );
105 :
106 0 : return indices; /// \returns the indices of the vector in sorted order.
107 0 : }
108 :
109 : /// Calculate the sum of a vector.
110 : /**
111 : *
112 : * \returns the sum of vec
113 : *
114 : * \tparam vectorT the std::vector type of vec
115 : *
116 : */
117 : template <typename valueT>
118 : valueT vectorSum( const valueT *vec, ///< [in] the vector
119 : size_t sz ///< [in] the size of the vector
120 : )
121 : {
122 : valueT sum = 0;
123 :
124 : for( size_t i = 0; i < sz; ++i )
125 : {
126 : sum += vec[i];
127 : }
128 :
129 : return sum;
130 : }
131 :
132 : /// Calculate the sum of a vector.
133 : /**
134 : *
135 : * \returns the sum of vec
136 : *
137 : * \tparam vectorT the std::vector type of vec
138 : *
139 : */
140 : template <typename vectorT>
141 : typename vectorT::value_type vectorSum( const vectorT &vec /**< [in] the vector */ )
142 : {
143 : typename vectorT::value_type sum = 0;
144 :
145 : for( size_t i = 0; i < vec.size(); ++i )
146 : {
147 : sum += vec[i];
148 : }
149 :
150 : return sum;
151 : }
152 :
153 : /// Calculate the mean of a vector.
154 : /**
155 : *
156 : * \returns the mean of vec
157 : *
158 : * \tparam vectorT the std::vector type of vec
159 : *
160 : */
161 : template <typename valueT>
162 0 : valueT vectorMean( const valueT *vec, ///< [in] the vector
163 : size_t sz ///< [in] the size of the vector
164 : )
165 : {
166 0 : valueT mean = 0;
167 :
168 0 : for( size_t i = 0; i < sz; ++i )
169 : {
170 0 : mean += vec[i];
171 : }
172 :
173 0 : mean /= sz;
174 :
175 0 : return mean;
176 : }
177 :
178 : /// Calculate the mean of a vector.
179 : /**
180 : *
181 : * \returns the mean of vec
182 : *
183 : * \tparam vectorT the std::vector type of vec
184 : *
185 : */
186 : template <typename vectorT>
187 36 : typename vectorT::value_type vectorMean( const vectorT &vec /**< [in] the vector */ )
188 : {
189 36 : typename vectorT::value_type mean = 0;
190 :
191 280 : for( size_t i = 0; i < vec.size(); ++i )
192 : {
193 244 : mean += vec[i];
194 : }
195 :
196 36 : mean /= vec.size();
197 :
198 36 : return mean;
199 : }
200 :
201 : /// Calculate the weighted mean of a vector.
202 : /**
203 : * \returns the weighted mean of vec
204 : *
205 : * \tparam vectorT the std::vector type of vec and w
206 : *
207 : */
208 : template <typename vectorT>
209 4 : typename vectorT::value_type vectorMean( const vectorT &vec, ///< [in] the vector */
210 : const vectorT &w ///< [in] the weights
211 : )
212 : {
213 4 : typename vectorT::value_type mean = 0, wsum = 0;
214 :
215 16 : for( size_t i = 0; i < vec.size(); ++i )
216 : {
217 12 : mean += w[i] * vec[i];
218 12 : wsum += w[i];
219 : }
220 :
221 4 : mean /= wsum;
222 :
223 4 : return mean;
224 : }
225 :
226 : /// Calculate median of a vector in-place, altering the vector.
227 : /** Returns the center element if vec has an odd number of elements. Returns the mean of the center 2 elements if vec
228 : * has an even number of elements.
229 : *
230 : * \returns the median of vec
231 : *
232 : * \tparam vectorT the std::vector type of vec
233 : *
234 : */
235 : template <typename vectorT>
236 109 : typename vectorT::value_type vectorMedianInPlace( vectorT &vec /**< [in] the vector, is altered by std::nth_element*/ )
237 : {
238 : typename vectorT::value_type med;
239 :
240 109 : int n = 0.5 * vec.size();
241 :
242 218 : std::nth_element( vec.begin(), vec.begin() + n, vec.end() );
243 :
244 109 : med = vec[n];
245 :
246 : // Average two points if even number of points
247 109 : if( vec.size() % 2 == 0 )
248 : {
249 147 : med = 0.5 * ( med + *std::max_element( vec.begin(), vec.begin() + n ) );
250 : }
251 :
252 109 : return med;
253 : }
254 :
255 : /// Calculate median of a vector, leaving the vector unaltered.
256 : /** Returns the center element if vec has an odd number of elements. Returns the mean of the center 2 elements if vec
257 : * has an even number of elements.
258 : *
259 : * \returns the median of vec
260 : *
261 : * \tparam vectorT the std::vector type of vec
262 : *
263 : */
264 : template <typename vectorT>
265 4 : typename vectorT::value_type vectorMedian( const vectorT &vec, ///< [in] the vector for which the median is desired
266 : vectorT *work = 0 /**< [in] [optional] an optional vector to use as
267 : workspace, use to avoid re-allocation*/
268 : )
269 : {
270 : typename vectorT::value_type med;
271 :
272 4 : bool localWork = false;
273 4 : if( work == 0 )
274 : {
275 0 : work = new vectorT;
276 0 : localWork = true;
277 : }
278 :
279 4 : work->resize( vec.size() );
280 :
281 16 : for( int i = 0; i < vec.size(); ++i )
282 : {
283 12 : ( *work )[i] = vec[i];
284 : }
285 :
286 4 : med = vectorMedianInPlace( *work );
287 :
288 4 : if( localWork )
289 : {
290 0 : delete work;
291 : }
292 :
293 4 : return med;
294 : }
295 :
296 : /// Calculate the variance of a vector relative to a supplied mean value.
297 : /**
298 : * \returns the variance of vec w.r.t. mean
299 : *
300 : * \tparam valueT the data type
301 : *
302 : */
303 : template <typename valueT>
304 0 : valueT vectorVariance( const valueT *vec, ///< [in] the vector
305 : size_t sz, ///< [in] the size of the vector
306 : valueT mean ///< [in] the mean value with which to calculate the variance
307 : )
308 : {
309 : valueT var;
310 :
311 0 : var = 0;
312 0 : for( size_t i = 0; i < sz; ++i )
313 : {
314 0 : var += pow( vec[i] - mean, 2 );
315 : }
316 :
317 0 : var /= ( sz - 1 );
318 :
319 0 : return var;
320 : }
321 :
322 : /// Calculate the variance of a vector relative to a supplied mean value.
323 : /**
324 : * \returns the variance of vec w.r.t. mean
325 : *
326 : * \tparam vectorT the std::vector type of vec
327 : *
328 : */
329 : template <typename vectorT>
330 : typename vectorT::value_type
331 235 : vectorVariance( const vectorT &vec, ///< [in] the vector
332 : const typename vectorT::value_type &mean ///< [in] the mean value with which to calculate the variance
333 : )
334 : {
335 : typename vectorT::value_type var;
336 :
337 235 : var = 0;
338 614878 : for( size_t i = 0; i < vec.size(); ++i )
339 : {
340 614643 : var += ( vec[i] - mean ) * ( vec[i] - mean );
341 : }
342 :
343 235 : var /= ( vec.size() - 1 );
344 :
345 235 : return var;
346 : }
347 :
348 : /// Calculate the variance of a vector.
349 : /**
350 : * \returns the variance of vec
351 : *
352 : * \tparam vectorT the std::vector type of vec
353 : *
354 : */
355 : template <typename valueT>
356 0 : valueT vectorVariance( const valueT *vec, ///< [in] the vector
357 : size_t sz ///< [in] the size of the vector
358 : )
359 : {
360 : valueT mean;
361 0 : mean = vectorMean( vec, sz );
362 :
363 0 : return vectorVariance( vec, sz, mean );
364 : }
365 :
366 : /// Calculate the variance of a vector.
367 : /**
368 : * \returns the variance of vec
369 : *
370 : * \tparam vectorT the std::vector type of vec
371 : *
372 : * \overload
373 : */
374 : template <typename vectorT>
375 1 : typename vectorT::value_type vectorVariance( const vectorT &vec /**< [in] the vector */ )
376 : {
377 : typename vectorT::value_type mean;
378 1 : mean = vectorMean( vec );
379 :
380 1 : return vectorVariance( vec, mean );
381 : }
382 :
383 : /// Calculate the sigma-clipped mean of a vector
384 : /** Performas sigma-clipping relative to the median, removing any values with deviation from the median > sigma.
385 : * Continues until either no values are removed, or maxPasses iterations. If maxPasses == 0, then it is ignored.
386 : *
387 : * \returns the sigma clipped mean of vec
388 : *
389 : * \tparam vectorT the std::vector type of vec
390 : * \tparam sigmaT the type of sigma, which is converted to value_type
391 : *
392 : */
393 : template <typename vectorT, typename sigmaT>
394 4 : typename vectorT::value_type vectorSigmaMean( const vectorT &vec, ///< [in] the vector (unaltered)
395 : const vectorT *weights, ///< [in] [optional] the weights (unaltered)
396 : const sigmaT &sigma, /**< [in] the standard deviation threshold to
397 : apply. */
398 : int &maxPasses /**< [in/out] [optional] the maximum number of
399 : sigma-clipping passes. Is set
400 : to actual number of passes on
401 : return. */
402 : )
403 : {
404 4 : vectorT work, wwork;
405 :
406 : typename vectorT::value_type med, var, Vsig, dev;
407 :
408 4 : bool doWeight = false;
409 4 : if( weights )
410 : {
411 2 : if( weights->size() == vec.size() )
412 : {
413 2 : doWeight = true;
414 : }
415 : }
416 :
417 4 : Vsig = sigma * sigma;
418 :
419 4 : med = vectorMedian( vec, &work );
420 4 : var = vectorVariance( work, med );
421 :
422 : int nclip;
423 4 : int passes = 0;
424 :
425 : // If weighting, have to synchronize work with weights since work will be
426 : // partially sorted by median.
427 4 : if( doWeight )
428 : {
429 2 : wwork.resize( vec.size() );
430 8 : for( int i = 0; i < vec.size(); ++i )
431 : {
432 6 : work[i] = vec[i];
433 6 : wwork[i] = ( *weights )[i];
434 : }
435 : }
436 :
437 4 : while( passes < maxPasses || maxPasses == 0 )
438 : {
439 4 : ++passes;
440 :
441 4 : nclip = 0;
442 :
443 16 : for( size_t i = 0; i < work.size(); ++i )
444 : {
445 12 : dev = pow( work[i] - med, 2 ) / var;
446 12 : if( dev > Vsig )
447 : {
448 0 : work.erase( work.begin() + i );
449 :
450 0 : if( doWeight )
451 0 : wwork.erase( wwork.begin() + i );
452 :
453 0 : --i;
454 0 : ++nclip;
455 : }
456 : }
457 :
458 4 : if( nclip == 0 )
459 4 : break;
460 0 : med = vectorMedian( work );
461 0 : var = vectorVariance( work, med );
462 : }
463 :
464 4 : maxPasses = passes;
465 :
466 4 : if( doWeight )
467 : {
468 2 : return vectorMean( work, wwork );
469 : }
470 : else
471 : {
472 2 : return vectorMean( work );
473 : }
474 4 : }
475 :
476 : /**
477 : * \overload
478 : */
479 : template <typename vectorT>
480 : typename vectorT::value_type
481 2 : vectorSigmaMean( const vectorT &vec, ///< [in] the vector (unaltered)
482 : typename vectorT::value_type sigma ///< [in] the standard deviation threshold to apply.
483 : )
484 : {
485 2 : int maxPasses = 0;
486 :
487 4 : return vectorSigmaMean( vec, (vectorT *)0, sigma, maxPasses );
488 : }
489 :
490 : /**
491 : * \overload
492 : */
493 : template <typename vectorT>
494 : typename vectorT::value_type
495 : vectorSigmaMean( const vectorT &vec, ///< [in] the vector (unaltered)
496 : const vectorT &weights, ///< [in] [optional] the weights (unaltered)
497 : typename vectorT::value_type sigma, ///< [in] the standard deviation threshold to apply.
498 : int &maxPasses /**< [in/out] [optional] the maximum number of sigma-clipping
499 : passes. Set to actual number of passes on return.*/
500 : )
501 : {
502 :
503 : return vectorSigmaMean( vec, &weights, sigma, maxPasses );
504 : }
505 :
506 : /**
507 : * \overload
508 : */
509 : template <typename vectorT>
510 : typename vectorT::value_type
511 2 : vectorSigmaMean( const vectorT &vec, ///< [in] the vector (unaltered)
512 : const vectorT &weights, ///< [in] [optional] the weights (unaltered)
513 : typename vectorT::value_type sigma ///< [in] the standard deviation threshold to apply.
514 : )
515 : {
516 2 : int maxPasses = 0;
517 :
518 4 : return vectorSigmaMean( vec, &weights, sigma, maxPasses );
519 : }
520 :
521 : /// Subtract a constant value from a vector
522 : template <typename valueT, typename constT>
523 0 : void vectorSub( valueT *vec, ///< [in/out] the vector, each element will have the constant subtracted from it
524 : size_t sz, ///< [in] the size of the vector
525 : const constT &c ///< [in] the constant to subtract from each element
526 : )
527 : {
528 0 : for( size_t n = 0; n < sz; ++n )
529 0 : vec[n] -= c;
530 0 : }
531 :
532 : /// Subtract a constant value from a vector
533 : template <typename vecT, typename constT>
534 : void vectorSub( vecT &vec, ///< [in/out] the vector, each element will have the constant subtracted from it
535 : const constT &c ///< [in] the constant to subtract from each element
536 : )
537 : {
538 : vectorSub( vec.data(), vec.size(), c );
539 : }
540 :
541 : /// Subtract the mean from a vector
542 : template <typename valueT>
543 0 : void vectorMeanSub( valueT *vec, ///< [in/out] the vector, each element will have the mean subtracted from it
544 : size_t sz ///< [in] the vector size
545 : )
546 : {
547 0 : valueT m = vectorMean( vec, sz );
548 0 : vectorSub( vec, sz, m );
549 0 : }
550 :
551 : /// Subtract the mean from a vector
552 : template <typename vecT>
553 0 : void vectorMeanSub( vecT &vec /**< [in/out] the vector, each element will have the mean subtracted from it*/ )
554 : {
555 0 : vectorMeanSub( vec.data(), vec.size() );
556 0 : }
557 :
558 : /// Subtract the median from a vector
559 : template <typename vecT>
560 : void vectorMedianSub( vecT &vec /**< [in/out] the vector, each element will have the median subtracted from it*/ )
561 : {
562 : typename vecT::value_type m = vectorMedian( vec );
563 : vectorSub( vec, m );
564 : }
565 :
566 : /// Smooth a vector using the mean in a window specified by its full-width
567 : template <typename realT>
568 : int vectorSmoothMean( realT *smVec, ///< [out] the smoothed version of the vector. At least as large as \p vec.
569 : realT *vec, ///< [in] the input vector, unaltered.
570 : size_t vecSize, ///< [in] the size of \p vec
571 : int win /**< [in] the full-width of the smoothing window.
572 : Should be even. 0 results in a slow memcpy.*/
573 : )
574 : {
575 : realT sum;
576 : int n;
577 : for( int i = 0; i < vecSize; ++i )
578 : {
579 : int j = i - 0.5 * win;
580 : if( j < 0 )
581 : {
582 : j = 0;
583 : }
584 :
585 : sum = 0;
586 : n = 0;
587 : while( j <= i + 0.5 * win && j < vecSize )
588 : {
589 : sum += vec[j];
590 : ++n;
591 : ++j;
592 : }
593 :
594 : smVec[i] = sum / n;
595 : }
596 :
597 : return 0;
598 : }
599 :
600 : /// Smooth a vector using the mean in a window specified by its full-width
601 : /** \overload
602 : */
603 : template <typename realT>
604 : int vectorSmoothMean( std::vector<realT> &smVec, ///< [out] the smoothed version of the vector. Will be resize()-ed.
605 : std::vector<realT> &vec, ///< [in] the input vector, unaltered.
606 : int win /**< [in] the full-width of the smoothing window.
607 : Should be even. 0 results in a slow memcpy.*/
608 : )
609 : {
610 : smVec.resize( vec.size() );
611 : return vectorSmoothMean( smVec.data(), vec.data(), vec.size(), win );
612 : }
613 :
614 : /// Smooth a vector using the mean in windows specified by their full-widths
615 : /** You supply a window width for each point. This is useful for, say, logarithmically growing bin sizes in a
616 : * PSD.
617 : *
618 : * \returns 0 on success
619 : * \returns -1 but who are we kidding it we don't check for errors
620 : *
621 : */
622 : template <typename realT>
623 : int vectorSmoothMean( std::vector<realT> &smVec, ///< [out] the smoothed version of the vector
624 : std::vector<realT> &vec, ///< [in] the input vector, unaltered.
625 : std::vector<int> &wins, ///< [in] the full-widths of the smoothing windows, same size as vec.
626 : bool norm = false /**< [in] if true the output will normalized to have the same
627 : integral as the input.*/
628 : )
629 : {
630 : smVec = vec;
631 :
632 : realT sum;
633 : int n;
634 : for( int i = 0; i < vec.size(); ++i )
635 : {
636 : int j = i - 0.5 * wins[i];
637 : if( j < 0 )
638 : {
639 : j = 0;
640 : }
641 :
642 : sum = 0;
643 : n = 0;
644 : while( j <= i + 0.5 * wins[i] && j < vec.size() )
645 : {
646 : sum += vec[j];
647 : ++n;
648 : ++j;
649 : }
650 :
651 : smVec[i] = sum / n;
652 : }
653 :
654 : if( norm )
655 : {
656 : realT sumin = 0, sumout = 0;
657 :
658 : for( int i = 0; i < vec.size(); ++i )
659 : {
660 : sumin += vec[i];
661 : sumout += smVec[i];
662 : }
663 :
664 : for( int i = 0; i < vec.size(); ++i )
665 : {
666 : smVec[i] *= sumin / sumout;
667 : }
668 : }
669 :
670 : return 0;
671 : }
672 :
673 : /// Smooth a vector using the median in a window specified by its full width.
674 : /** For even widths, the window is associated with the higher-index member of the central pair: it contains win/2
675 : * samples before the output sample and one fewer after it. Windows are truncated at the vector boundaries.
676 : *
677 : * \returns 0 on success
678 : * \returns -1 if win is not positive
679 : */
680 : template <typename realT>
681 2 : int vectorSmoothMedian( std::vector<realT> &smVec, /**< [out] the smoothed version of the vector */
682 : std::vector<realT> &vec, /**< [in] the input vector, unaltered */
683 : int win /**< [in] the full width of the smoothing window */
684 : )
685 : {
686 2 : if( win <= 0 )
687 : {
688 1 : return -1;
689 : }
690 :
691 1 : smVec = vec;
692 :
693 1 : const int before = win / 2;
694 1 : const int after = win - before - 1;
695 1 : std::vector<realT> tvec;
696 1 : tvec.reserve( std::min<size_t>( static_cast<size_t>( win ), vec.size() ) );
697 :
698 6 : for( int i = 0; i < vec.size(); ++i )
699 : {
700 5 : const int first = std::max( 0, i - before );
701 5 : const int last = std::min<int>( vec.size(), i + after + 1 );
702 :
703 5 : tvec.clear();
704 21 : for( int j = first; j < last; ++j )
705 : {
706 16 : tvec.push_back( vec[j] );
707 : }
708 :
709 5 : smVec[i] = vectorMedianInPlace( tvec );
710 : }
711 :
712 1 : return 0;
713 1 : }
714 :
715 : /// Smooth a vector using the max in a window specified by its full-width
716 : template <typename realT>
717 : int vectorSmoothMax( std::vector<realT> &smVec, ///< [out] the smoothed version of the vector
718 : std::vector<realT> &vec, ///< [in] the input vector, unaltered.
719 : int win ///< [in] the full-width of the smoothing window
720 : )
721 : {
722 : smVec = vec;
723 :
724 : for( int i = 0; i < vec.size(); ++i )
725 : {
726 : int j = i - 0.5 * win;
727 : if( j < 0 )
728 : {
729 : j = 0;
730 : }
731 :
732 : smVec[i] = vec[j];
733 : ++j;
734 : while( j <= i + 0.5 * win && j < vec.size() )
735 : {
736 : if( vec[j] > smVec[i] )
737 : {
738 : smVec[i] = vec[j];
739 : }
740 : ++j;
741 : }
742 : }
743 :
744 : return 0;
745 : }
746 :
747 : /// Re-bin a vector by summing (or averaging) in bins of size n points.
748 : /**
749 : * \returns 0 on success
750 : * \returns -1 on error
751 : *
752 : * \tparam vectorT is any vector-like type with resize(), size(), and the operator()[].
753 : */
754 : template <typename vectorT>
755 0 : int vectorRebin( vectorT &binv, ///< [out] the re-binned vector. will be resized.
756 : const vectorT &v, ///< [in] the vector to bin.
757 : unsigned n, ///< [in] the size of the bins, in points
758 : bool binMean = false /**< [in] [optional] flag controlling whether sums (false) or means
759 : (true) are calculated.*/
760 : )
761 : {
762 0 : if( n == 0 )
763 : {
764 0 : return -1;
765 : }
766 :
767 0 : binv.resize( v.size() / n );
768 :
769 0 : for( size_t i = 0; i < binv.size(); ++i )
770 : {
771 0 : binv[i] = 0;
772 :
773 : unsigned j;
774 0 : for( j = 0; j < n; ++j )
775 : {
776 0 : if( i * n + j >= v.size() )
777 : {
778 0 : break;
779 : }
780 :
781 0 : binv[i] += v[i * n + j];
782 : }
783 :
784 0 : if( binMean )
785 : {
786 0 : if( j > 0 )
787 : {
788 0 : binv[i] /= j;
789 : }
790 : }
791 : }
792 :
793 0 : return 0;
794 : }
795 :
796 : /// Calculate and accumulate the means of a timeseries in bins of various sizes.
797 : /** Useful mainly to calculate the variance of the mean as a function of sample size.
798 : * The output is a vector of vectors, where each element is a vector which contains the means in the
799 : * unique bins of size corresponding to the same index in the in put binSzs vector.
800 : *
801 : * \returns 0 on success.
802 : */
803 : template <typename vectorT, typename binVectorT>
804 0 : int vectorBinMeans( std::vector<vectorT> &means, /**< [out] the means in each distinct bin. Not cleared, but Will be
805 : resized with new means appended.*/
806 : binVectorT &binSzs, ///< [in] the bin sizes in which to calculate the means
807 : const vectorT &v ///< [in] the input vector to bin .
808 : )
809 : {
810 0 : vectorT binv;
811 :
812 0 : means.resize( binSzs.size() );
813 :
814 0 : for( size_t i = 0; i < binSzs.size(); ++i )
815 : {
816 0 : vectorRebin( binv, v, binSzs[i], true );
817 :
818 0 : means[i].resize( means[i].size() + binv.size() );
819 0 : for( size_t j = 0; j < binv.size(); ++j )
820 : {
821 0 : means[i][means[i].size() - binv.size() + j] = binv[j];
822 : }
823 : }
824 :
825 0 : return 0;
826 0 : }
827 :
828 : /// Convolve (smooth) a vector with a Gaussian.
829 : /**
830 : * \returns 0 on success.
831 : * \returns -1 on error.
832 : *
833 : * \tparam realT is the real floating point type of the data.
834 : */
835 : template <typename realT, typename fwhmT, typename winhwT>
836 : int vectorGaussConvolve( std::vector<realT> &dataOut, ///< [out] The smoothed data vector. Resized.
837 : const std::vector<realT> &dataIn, ///< [in] The data vector to smooth.
838 : const std::vector<realT> &scale, ///< [in] The scale vector used to calculate the kernel.
839 : const fwhmT fwhm, /**< [in] The FWHM of the Gaussian kernel,
840 : same units as scale.*/
841 : const winhwT winhw ///< [in] The window half-width in pixels.
842 : )
843 : {
844 :
845 : if( dataIn.size() != scale.size() )
846 : {
847 : return -1;
848 : }
849 :
850 : realT _fwhm = fwhm;
851 : int _winhw = winhw;
852 :
853 : dataOut.resize( dataIn.size() );
854 :
855 : realT sigma = func::fwhm2sigma<realT>( _fwhm );
856 :
857 : for( int i = 0; i < dataIn.size(); ++i )
858 : {
859 : realT G;
860 : realT sum = 0;
861 : realT norm = 0;
862 : for( int j = i - _winhw; j < i + _winhw; ++j )
863 : {
864 : if( j < 0 )
865 : {
866 : continue;
867 : }
868 : if( j > dataIn.size() - 1 )
869 : {
870 : continue;
871 : }
872 :
873 : G = func::gaussian<realT>( scale[j], 0.0, 1.0, scale[i], sigma );
874 :
875 : sum += G * dataIn[j];
876 : norm += G;
877 : }
878 :
879 : dataOut[i] = sum / norm;
880 : }
881 :
882 : return 0;
883 : }
884 :
885 : /// Calculate a cumulative histogram of a vector.
886 : /** Sorts the vector and sums.
887 : *
888 : * \retval 0 on success, -1 otherwise.
889 : *
890 : * \tparam floatT the floating point type of the vector contens.
891 : */
892 : template <typename floatT>
893 : int vectorCumHist( std::vector<floatT> &svec, ///< [out] Contains the sorted vector.
894 : std::vector<floatT> &sum, ///< [out] Contains the cumulative or running sum of the sorted vector
895 : std::vector<floatT> &vec ///< [in] The vector to sort and sum.
896 : )
897 : {
898 : svec = vec;
899 :
900 : std::sort( svec.begin(), svec.end() );
901 :
902 : sum.resize( svec.size() );
903 :
904 : sum[0] = svec[0];
905 :
906 : for( int i = 1; i < svec.size(); ++i )
907 : {
908 : sum[i] = sum[i - 1] + svec[i];
909 : }
910 :
911 : return 0;
912 : }
913 :
914 : /// Calculate a reverse cumulative histogram of a vector.
915 : /** Reverse-sorts the vector and sums.
916 : *
917 : * \retval 0 on success, -1 otherwise.
918 : *
919 : * \tparam floatT the floating point type of the vector contens.
920 : */
921 : template <typename floatT>
922 : int vectorCumHistReverse(
923 : std::vector<floatT> &svec, ///< [out] Contains the reverse-sorted vector.
924 : std::vector<floatT> &sum, ///< [out] Contains the cumulative or running sum of the reverse-sorted vector
925 : std::vector<floatT> &vec ///< [in] The vector to reverse-sort and sum.
926 : )
927 : {
928 : svec = vec;
929 :
930 : std::sort( svec.begin(), svec.end(), std::greater<floatT>() );
931 :
932 : sum.resize( svec.size() );
933 :
934 : sum[0] = svec[0];
935 :
936 : for( int i = 1; i < svec.size(); ++i )
937 : {
938 : sum[i] = sum[i - 1] + svec[i];
939 : }
940 :
941 : return 0;
942 : }
943 :
944 : ///@}
945 :
946 : } // namespace math
947 : } // namespace mx
948 :
949 : #endif // math_vectorUtils_hpp
|
