1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
#include "utility.h"
#include <algorithm>
#include <cmath>
namespace meow{
template<class T>
inline T noEPS(T value, T eps){
T epsp((eps < T(0)) ? -eps : eps);
return ((value < -epsp || value > epsp)
? value
: T(0));
}
template<class T>
inline T normalize(T lower, T upper, T value){
return (value - lower) / (upper - lower);
}
template<class T>
inline T denormalize(T lower, T upper, T _ratio){
return lower + _ratio * (upper - lower);
}
template<class T>
inline T ratioMapping(T l1, T u1, T m1, T l2, T u2){
return denormalize(l2, u2, normalize(l1, u1, m1));
}
template<class T>
inline T inRange(T const& mn, T const& mx, T const& v){
return std::min(mx, std::max(mn, v));
}
template<class T>
inline T squ(T const& x){
return x * x;
}
template<class T>
inline T cub(T const& x){
return x * x * x;
}
template<class T>
inline double average( T const& beg, T const& end, double sigs){
int N = 0;
double av = 0;
for(T it = beg; it != end; it++, N++){
av += *it;
}
av /= N;
double sig = 0;
for(T it = beg; it != end; it++){
sig += (*it - av) * (*it - av);
}
sig = sqrt(sig / N);
double lower = av - sig * sigs, upper = av + sig * sigs;
double ret = 0, retn = 0;
for(T it = beg; it != end; it++){
if(lower <= *it && *it <= upper){
ret += *it;
retn++;
}
}
return ret / retn;
}
template<class T>
inline double average( T const& beg, T const& end, T const& p, double sigs){
int N = 0;
double ps = 0;
for(T it = beg, ip = p; it != end; it++, N++, ip++){
ps += *ip;
}
double av = 0;
for(T it = beg, ip = p; it != end; it++, ip++){
av += *it * *ip / ps;
}
double sig = 0;
for(T it = beg, ip = p; it != end; it++, ip++){
sig += *ip / ps * (*it - av) * (*it - av);
}
sig = sqrt(sig);
double lower = av - sig * sigs, upper = av + sig * sigs;
double ret = 0, retn = 0;
for(T it = beg, ip = p; it != end; it++, ip++){
if(lower <= *it && *it <= upper){
ret += *it * *ip;
retn += *ip;
}
}
if(retn <= 1e-10) return av;
return ret / retn;
}
}
|
