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