/*
* c++ implementation of Lowess weighted regression by
* Peter Glaus http://www.cs.man.ac.uk/~glausp/
*
*
* Based on fortran code by Cleveland downloaded from:
* http://netlib.org/go/lowess.f
* original author:
* wsc@research.bell-labs.com Mon Dec 30 16:55 EST 1985
* W. S. Cleveland
* Bell Laboratories
* Murray Hill NJ 07974
*
* See original documentation below the code for details.
*
*/
#include<algorithm>
#include<cmath>
#include<fstream>
using namespace std;
#include "lowess.h"
#include "common.h"
void lowess(const vector<double> &x, const vector<double> &y, double f, long nsteps, vector<double> &ys){//{{{
vector<double> rw,res;
lowess(x,y,f,nsteps,0.,ys,rw,res);
}//}}}
void lowess(const vector<double> &x, const vector<double> &y, double f, long nsteps, double delta, vector<double> &ys, vector<double> &rw, vector<double>&res){ //{{{
long n=(long)x.size();
bool ok=false;
long nleft,nright, i, j, iter, last, m1, m2, ns;
double cut, cmad, r, d1, d2, c1, c9, alpha, denom;
if((n==0)||((long)y.size()!=n)) return;
ys.resize(n);
rw.resize(n);
res.resize(n);
if(n==1){
ys[0]=y[0];
return;
}
// ns - at least 2, at most n
ns = max(min((long)(f*n),n),(long)2);
for(iter=0;iter<nsteps+1; iter++){
// robustnes iterations
nleft = 0;
nright = ns-1;
// index of last estimated point
last = -1;
// index of current point
i=0;
do{
while(nright<n-1){
// move <nleft,nright> right, while radius decreases
d1 = x[i]-x[nleft];
d2 = x[nright+1] - x[i];
if(d1<=d2)break;
nleft++;
nright++;
}
// fit value at x[i]
lowest(x,y,x[i],ys[i],nleft,nright,res,iter>0,rw,ok);
if(!ok) ys[i]=y[i];
if(last<i-1){
// interpolate skipped points
if(last<0){
warning("Lowess: out of range.\n");
}
denom = x[i] - x[last];
for(j=last+1;j<i;j++){
alpha = (x[j]-x[last])/denom;
ys[j] = alpha * ys[i] + (1.0-alpha)*ys[last];
}
}
last = i;
cut = x[last]+delta;
for(i=last+1;i<n;i++){
if(x[i]>cut)break;
if(x[i]==x[last]){
ys[i]=ys[last];
last=i;
}
}
i=max(last+1,i-1);
}while(last<n-1);
for(i=0;i<n;i++)
res[i] = y[i]-ys[i];
if(iter==nsteps)break ;
for(i=0;i<n;i++)
rw[i]=abs(res[i]);
sort(rw.begin(),rw.end());
m1 = n/2+1;
m2 = n-m1;
m1 --;
cmad = 3.0 *(rw[m1]+rw[m2]);
c9 = .999*cmad;
c1 = .001*cmad;
for(i=0;i<n;i++){
r = abs(res[i]);
if(r<=c1) rw[i]=1;
else if(r>c9) rw[i]=0;
else rw[i] = (1.0-(r/cmad)*(r/cmad))*(1.0-(r/cmad)*(r/cmad));
}
}
}//}}}
void lowest(const vector<double> &x, const vector<double> &y, double xs, double &ys, long nleft, long nright, vector<double> &w, bool userw, vector<double> &rw, bool &ok){//{{{
long n = (long)x.size();
long nrt, j;
double a, b, c, h, r, h1, h9, range;
range = x[n-1]-x[0];
h = max(xs-x[nleft],x[nright]-xs);
h9 = 0.999*h;
h1 = 0.001*h;
// sum of weights
a = 0;
for(j=nleft;j<n;j++){
// compute weights (pick up all ties on right)
w[j]=0.;
r = abs(x[j]-xs);
if(r<=h9){
// small enough for non-zero weight
if(r>h1) w[j] = (1.0-(r/h)*(r/h)*(r/h))*(1.0-(r/h)*(r/h)*(r/h))*(1.0-(r/h)*(r/h)*(r/h));
else w[j] = 1.;
if(userw) w[j] *= rw[j];
a += w[j];
}else if(x[j]>xs) break; // get out at first zero wt on right
}
nrt = j-1;
// rightmost pt (may be greater than nright because of ties)
if(a<=0.) ok = false;
else{
// weighted least squares
ok = true;
// normalize weights
for(j=nleft;j<=nrt;j++)
w[j] /= a;
if(h>0.){
// use linear fit
a = 0.;
for(j=nleft;j<=nrt;j++)
a += w[j]*x[j]; // weighted centre of values
b = xs-a;
c = 0;
for(j=nleft;j<=nrt;j++)
c += w[j]*(x[j]-a)*(x[j]-a);
if(sqrt(c)>0.001*range){
// points are spread enough to compute slope
b /= c;
for(j=nleft;j<=nrt;j++)
w[j] *= (1.0+b*(x[j]-a));
}
}
ys = 0;
for(j=nleft;j<=nrt;j++)
ys += w[j]*y[j];
}
}//}}}
/* {{{ Documentation
* wsc@research.bell-labs.com Mon Dec 30 16:55 EST 1985
* W. S. Cleveland
* Bell Laboratories
* Murray Hill NJ 07974
*
* outline of this file:
* lines 1-72 introduction
* 73-177 documentation for lowess
* 178-238 ratfor version of lowess
* 239-301 documentation for lowest
* 302-350 ratfor version of lowest
* 351-end test driver and fortran version of lowess and lowest
*
* a multivariate version is available by "send dloess from a"
*
* COMPUTER PROGRAMS FOR LOCALLY WEIGHTED REGRESSION
*
* This package consists of two FORTRAN programs for
* smoothing scatterplots by robust locally weighted
* regression, or lowess. The principal routine is LOWESS
* which computes the smoothed values using the method
* described in The Elements of Graphing Data, by William S.
* Cleveland (Wadsworth, 555 Morego Street, Monterey,
* California 93940).
*
* LOWESS calls a support routine, LOWEST, the code for
* which is included. LOWESS also calls a routine SORT, which
* the user must provide.
*
* To reduce the computations, LOWESS requires that the
* arrays X and Y, which are the horizontal and vertical
* coordinates, respectively, of the scatterplot, be such that
* X is sorted from smallest to largest. The user must
* therefore use another sort routine which will sort X and Y
* according to X.
* To summarize the scatterplot, YS, the fitted values,
* should be plotted against X. No graphics routines are
* available in the package and must be supplied by the user.
*
* The FORTRAN code for the routines LOWESS and LOWEST has
* been generated from higher level RATFOR programs
* (B. W. Kernighan, ``RATFOR: A Preprocessor for a Rational
* Fortran,'' Software Practice and Experience, Vol. 5 (1975),
* which are also included.
*
* The following are data and output from LOWESS that can
* be used to check your implementation of the routines. The
* notation (10)v means 10 values of v.
*
*
*
*
* X values:
* 1 2 3 4 5 (10)6 8 10 12 14 50
*
* Y values:
* 18 2 15 6 10 4 16 11 7 3 14 17 20 12 9 13 1 8 5 19
*
*
* YS values with F = .25, NSTEPS = 0, DELTA = 0.0
* 13.659 11.145 8.701 9.722 10.000 (10)11.300 13.000 6.440 5.596
* 5.456 18.998
*
* YS values with F = .25, NSTEPS = 0 , DELTA = 3.0
* 13.659 12.347 11.034 9.722 10.511 (10)11.300 13.000 6.440 5.596
* 5.456 18.998
*
* YS values with F = .25, NSTEPS = 2, DELTA = 0.0
* 14.811 12.115 8.984 9.676 10.000 (10)11.346 13.000 6.734 5.744
* 5.415 18.998
*
*
*
*
* LOWESS
*
*
*
* Calling sequence
*
* CALL LOWESS(X,Y,N,F,NSTEPS,DELTA,YS,RW,RES)
*
* Purpose
*
* LOWESS computes the smooth of a scatterplot of Y against X
* using robust locally weighted regression. Fitted values,
* YS, are computed at each of the values of the horizontal
* axis in X.
*
* Argument description
*
* X = Input; abscissas of the points on the
* scatterplot; the values in X must be ordered
* from smallest to largest.
* Y = Input; ordinates of the points on the
* scatterplot.
* N = Input; dimension of X,Y,YS,RW, and RES.
* F = Input; specifies the amount of smoothing; F is
* the fraction of points used to compute each
* fitted value; as F increases the smoothed values
* become smoother; choosing F in the range .2 to
* .8 usually results in a good fit; if you have no
* idea which value to use, try F = .5.
* NSTEPS = Input; the number of iterations in the robust
* fit; if NSTEPS = 0, the nonrobust fit is
* returned; setting NSTEPS equal to 2 should serve
* most purposes.
* DELTA = input; nonnegative parameter which may be used
* to save computations; if N is less than 100, set
* DELTA equal to 0.0; if N is greater than 100 you
* should find out how DELTA works by reading the
* additional instructions section.
* YS = Output; fitted values; YS(I) is the fitted value
* at X(I); to summarize the scatterplot, YS(I)
* should be plotted against X(I).
* RW = Output; robustness weights; RW(I) is the weight
* given to the point (X(I),Y(I)); if NSTEPS = 0,
* RW is not used.
* RES = Output; residuals; RES(I) = Y(I)-YS(I).
*
*
* Other programs called
*
* LOWEST
* SSORT
*
* Additional instructions
*
* DELTA can be used to save computations. Very roughly the
* algorithm is this: on the initial fit and on each of the
* NSTEPS iterations locally weighted regression fitted values
* are computed at points in X which are spaced, roughly, DELTA
* apart; then the fitted values at the remaining points are
* computed using linear interpolation. The first locally
* weighted regression (l.w.r.) computation is carried out at
* X(1) and the last is carried out at X(N). Suppose the
* l.w.r. computation is carried out at X(I). If X(I+1) is
* greater than or equal to X(I)+DELTA, the next l.w.r.
* computation is carried out at X(I+1). If X(I+1) is less
* than X(I)+DELTA, the next l.w.r. computation is carried out
* at the largest X(J) which is greater than or equal to X(I)
* but is not greater than X(I)+DELTA. Then the fitted values
* for X(K) between X(I) and X(J), if there are any, are
* computed by linear interpolation of the fitted values at
* X(I) and X(J). If N is less than 100 then DELTA can be set
* to 0.0 since the computation time will not be too great.
* For larger N it is typically not necessary to carry out the
* l.w.r. computation for all points, so that much computation
* time can be saved by taking DELTA to be greater than 0.0.
* If DELTA = Range (X)/k then, if the values in X were
* uniformly scattered over the range, the full l.w.r.
* computation would be carried out at approximately k points.
* Taking k to be 50 often works well.
*
* Method
*
* The fitted values are computed by using the nearest neighbor
* routine and robust locally weighted regression of degree 1
* with the tricube weight function. A few additional features
* have been added. Suppose r is FN truncated to an integer.
* Let h be the distance to the r-th nearest neighbor
* from X(I). All points within h of X(I) are used. Thus if
* the r-th nearest neighbor is exactly the same distance as
* other points, more than r points can possibly be used for
* the smooth at X(I). There are two cases where robust
* locally weighted regression of degree 0 is actually used at
* X(I). One case occurs when h is 0.0. The second case
* occurs when the weighted standard error of the X(I) with
* respect to the weights w(j) is less than .001 times the
* range of the X(I), where w(j) is the weight assigned to the
* j-th point of X (the tricube weight times the robustness
* weight) divided by the sum of all of the weights. Finally,
* if the w(j) are all zero for the smooth at X(I), the fitted
* value is taken to be Y(I).
*
*
*
*
* subroutine lowess(x,y,n,f,nsteps,delta,ys,rw,res)
* real x(n),y(n),ys(n),rw(n),res(n)
* logical ok
* if (n<2){ ys(1) = y(1); return }
* ns = max0(min0(ifix(f*float(n)),n),2) # at least two, at most n points
* for(iter=1; iter<=nsteps+1; iter=iter+1){ # robustness iterations
* nleft = 1; nright = ns
* last = 0 # index of prev estimated point
* i = 1 # index of current point
* repeat{
* while(nright<n){
* # move nleft, nright to right if radius decreases
* d1 = x(i)-x(nleft)
* d2 = x(nright+1)-x(i)
* # if d1<=d2 with x(nright+1)==x(nright), lowest fixes
* if (d1<=d2) break
* # radius will not decrease by move right
* nleft = nleft+1
* nright = nright+1
* }
* call lowest(x,y,n,x(i),ys(i),nleft,nright,res,iter>1,rw,ok)
* # fitted value at x(i)
* if (!ok) ys(i) = y(i)
* # all weights zero - copy over value (all rw==0)
* if (last<i-1) { # skipped points -- interpolate
* denom = x(i)-x(last) # non-zero - proof?
* for(j=last+1; j<i; j=j+1){
* alpha = (x(j)-x(last))/denom
* ys(j) = alpha*ys(i)+(1.0-alpha)*ys(last)
* }
* }
* last = i # last point actually estimated
* cut = x(last)+delta # x coord of close points
* for(i=last+1; i<=n; i=i+1){ # find close points
* if (x(i)>cut) break # i one beyond last pt within cut
* if(x(i)==x(last)){ # exact match in x
* ys(i) = ys(last)
* last = i
* }
* }
* i=max0(last+1,i-1)
* # back 1 point so interpolation within delta, but always go forward
* } until(last>=n)
* do i = 1,n # residuals
* res(i) = y(i)-ys(i)
* if (iter>nsteps) break # compute robustness weights except last time
* do i = 1,n
* rw(i) = abs(res(i))
* call sort(rw,n)
* m1 = 1+n/2; m2 = n-m1+1
* cmad = 3.0*(rw(m1)+rw(m2)) # 6 median abs resid
* c9 = .999*cmad; c1 = .001*cmad
* do i = 1,n {
* r = abs(res(i))
* if(r<=c1) rw(i)=1. # near 0, avoid underflow
* else if(r>c9) rw(i)=0. # near 1, avoid underflow
* else rw(i) = (1.0-(r/cmad)**2)**2
* }
* }
* return
* end
*
*
*
*
* LOWEST
*
*
*
* Calling sequence
*
* CALL LOWEST(X,Y,N,XS,YS,NLEFT,NRIGHT,W,USERW,RW,OK)
*
* Purpose
*
* LOWEST is a support routine for LOWESS and ordinarily will
* not be called by the user. The fitted value, YS, is
* computed at the value, XS, of the horizontal axis.
* Robustness weights, RW, can be employed in computing the
* fit.
*
* Argument description
*
*
* X = Input; abscissas of the points on the
* scatterplot; the values in X must be ordered
* from smallest to largest.
* Y = Input; ordinates of the points on the
* scatterplot.
* N = Input; dimension of X,Y,W, and RW.
* XS = Input; value of the horizontal axis at which the
* smooth is computed.
* YS = Output; fitted value at XS.
* NLEFT = Input; index of the first point which should be
* considered in computing the fitted value.
* NRIGHT = Input; index of the last point which should be
* considered in computing the fitted value.
* W = Output; W(I) is the weight for Y(I) used in the
* expression for YS, which is the sum from
* I = NLEFT to NRIGHT of W(I)*Y(I); W(I) is
* defined only at locations NLEFT to NRIGHT.
* USERW = Input; logical variable; if USERW is .TRUE., a
* robust fit is carried out using the weights in
* RW; if USERW is .FALSE., the values in RW are
* not used.
* RW = Input; robustness weights.
* OK = Output; logical variable; if the weights for the
* smooth are all 0.0, the fitted value, YS, is not
* computed and OK is set equal to .FALSE.; if the
* fitted value is computed OK is set equal to
*
*
* Method
*
* The smooth at XS is computed using (robust) locally weighted
* regression of degree 1. The tricube weight function is used
* with h equal to the maximum of XS-X(NLEFT) and X(NRIGHT)-XS.
* Two cases where the program reverts to locally weighted
* regression of degree 0 are described in the documentation
* for LOWESS.
*
*
*
*
* subroutine lowest(x,y,n,xs,ys,nleft,nright,w,userw,rw,ok)
* real x(n),y(n),w(n),rw(n)
* logical userw,ok
* range = x(n)-x(1)
* h = amax1(xs-x(nleft),x(nright)-xs)
* h9 = .999*h
* h1 = .001*h
* a = 0.0 # sum of weights
* for(j=nleft; j<=n; j=j+1){ # compute weights (pick up all ties on right)
* w(j)=0.
* r = abs(x(j)-xs)
* if (r<=h9) { # small enough for non-zero weight
* if (r>h1) w(j) = (1.0-(r/h)**3)**3
* else w(j) = 1.
* if (userw) w(j) = rw(j)*w(j)
* a = a+w(j)
* }
* else if(x(j)>xs)break # get out at first zero wt on right
* }
* nrt=j-1 # rightmost pt (may be greater than nright because of ties)
* if (a<=0.0) ok = FALSE
* else { # weighted least squares
* ok = TRUE
* do j = nleft,nrt
* w(j) = w(j)/a # make sum of w(j) == 1
* if (h>0.) { # use linear fit
* a = 0.0
* do j = nleft,nrt
* a = a+w(j)*x(j) # weighted center of x values
* b = xs-a
* c = 0.0
* do j = nleft,nrt
* c = c+w(j)*(x(j)-a)**2
* if(sqrt(c)>.001*range) {
* # points are spread out enough to compute slope
* b = b/c
* do j = nleft,nrt
* w(j) = w(j)*(1.0+b*(x(j)-a))
* }
* }
* ys = 0.0
* do j = nleft,nrt
* ys = ys+w(j)*y(j)
* }
* return
* end
*
}}}*/
C语音版本
#define FALSE 0
#define TRUE 1
long min (long x, long y)
{
return((x < y) ? x : y);
}
long max (long x, long y)
{
return((x > y) ? x : y);
}
static double pow2(double x) { return(x * x); }
static double pow3(double x) { return(x * x * x); }
double fmax(double x, double y) { return (x > y ? x : y); }
int static compar(const void *aa, const void *bb)
{
const double* a = (double*)aa;
const double* b = (double*)bb;
if (*a < *b) return(-1);
else if (*a > *b) return(1);
else return(0);
}
static void lowest(const double *x, const double *y, size_t n, double xs, double *ys, long nleft, long nright,
double *w, int userw, double *rw, int *ok)
{
double range, h, h1, h9, a, b, c, r;
long j, nrt;
range = x[n - 1] - x[0];
h = fmax(xs - x[nleft], x[nright] - xs);
h9 = .999 * h;
h1 = .001 * h;
/* compute weights (pick up all ties on right) */
a = 0.0; /* sum of weights */
for(j = nleft; j < n; j++)
{
w[j]=0.0;
r = fabs(x[j] - xs);
if (r <= h9)
{ /* small enough for non-zero weight */
if (r > h1) w[j] = pow3(1.0-pow3(r/h));
else w[j] = 1.0;
if (userw) w[j] = rw[j] * w[j];
a += w[j];
}
else if (x[j] > xs) break; /* get out at first zero wt on right */
}
nrt = j - 1; /* rightmost pt (may be greater than nright because of ties) */
if (a <= 0.0) *ok = FALSE;
else { /* weighted least squares */
*ok = TRUE;
/* make sum of w[j] == 1 */
for (j = nleft; j <= nrt; j++) w[j] = w[j] / a;
if (h > 0.0) { /* use linear fit */
/* find weighted center of x values */
for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j];
b = xs - a;
for (j = nleft, c = 0.0; j <= nrt; j++)
c += w[j] * (x[j] - a) * (x[j] - a);
if(sqrt(c) > .001 * range) {
/* points are spread out enough to compute slope */
b = b/c;
for (j = nleft; j <= nrt; j++)
w[j] = w[j] * (1.0 + b*(x[j] - a));
}
}
for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j];
}
}
static void sort(double *x, size_t n)
{
qsort(x, n, sizeof(double), compar);
}
int lowess(const double *x, const double *y, size_t n,double f, size_t nsteps,
double delta, double *ys, double *rw, double *res)
{
int iter, ok;
long i, j, last, m1, m2, nleft, nright, ns;
double d1, d2, denom, alpha, cut, cmad, c9, c1, r;
if (n < 2) { ys[0] = y[0]; return(1); }
ns = max(min((long) (f * n), n), 2); /* at least two, at most n points */
for(iter = 1; iter <= nsteps + 1; iter++){ /* robustness iterations */
nleft = 0; nright = ns - 1;
last = -1; /* index of prev estimated point */
i = 0; /* index of current point */
do {
while(nright < n - 1){
/* move nleft, nright to right if radius decreases */
d1 = x[i] - x[nleft];
d2 = x[nright + 1] - x[i];
/* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
if (d1 <= d2) break;
/* radius will not decrease by move right */
nleft++;
nright++;
}
lowest(x, y,
n, x[i],
&ys[i],
nleft, nright,
res, (iter > 1), rw, &ok);
/* fitted value at x[i] */
if (! ok) ys[i] = y[i];
/* all weights zero - copy over value (all rw==0) */
if (last < i - 1) { /* skipped points -- interpolate */
denom = x[i] - x[last]; /* non-zero - proof? */
for(j = last + 1; j < i; j = j + 1){
alpha = (x[j] - x[last]) / denom;
ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last];
}
}
last = i; /* last point actually estimated */
cut = x[last] + delta; /* x coord of close points */
for(i=last + 1; i < n; i++) { /* find close points */
if (x[i] > cut) break; /* i one beyond last pt within cut */
if(x[i] == x[last]) { /* exact match in x */
ys[i] = ys[last];
last = i;
}
}
i = max(last + 1,i - 1);
/* back 1 point so interpolation within delta, but always go forward */
} while(last < n - 1);
for (i = 0; i < n; i++) /* residuals */
res[i] = y[i] - ys[i];
if (iter > nsteps) break; /* compute robustness weights except last time */
for (i = 0; i < n; i++)
rw[i] = fabs(res[i]);
sort(rw,n);
m1 = 1 + n / 2; m2 = n - m1 + 1;
cmad = 3.0 * (rw[m1] + rw[m2]); /* 6 median abs resid */
c9 = .999 * cmad; c1 = .001 * cmad;
for (i = 0; i < n; i++) {
r = fabs(res[i]);
if(r <= c1) rw[i] = 1.0; /* near 0, avoid underflow */
else if(r > c9) rw[i] = 0.0; /* near 1, avoid underflow */
else rw[i] = pow2(1.0 - pow2(r / cmad));
}
}
return(0);
}
测试main
void test_smooth()
{
const double in[] = {
-55.1047,
-56.3673,
-56.314,
-55.8626,
-56.3733,
-55.8694,
-54.4476,
-56.1106,
-58.3752,
-56.4632,
-57.5418,
-57.259,
-54.999,
-56.5903,
-56.5675,
-57.2702,
-56.5198,
-59.275,
-58.0706,
-55.4509,
-58.1618,
-57.2596,
-55.4555,
-55.7893,
-55.7953,
-55.8368,
-57.3642,
-57.5559,
-56.0514,
-56.8001,
-58.4702,
-55.175,
-53.7347,
-54.4519,
-54.5773,
-56.9627,
-57.1959,
-55.6704,
-55.0481,
-57.9426,
-57.3462,
-55.6331,
-55.629,
-55.3975,
-56.4719,
-58.1078,
-56.1705,
-54.758,
-58.8711,
-57.9153,
-56.4004,
-55.1685,
-55.7985,
-54.3574,
-55.6311,
-55.4626,
-56.6099,
-55.4795,
-54.0479,
-56.069,
-56.2238,
-56.643,
-57.3297,
-57.2569,
-55.3871,
-55.8629,
-55.827,
-55.3129,
-56.7753,
-54.7421,
-53.2383,
-53.1972,
-54.2125,
-55.1294,
-55.3516,
-54.4107,
-56.1692,
-55.6607,
-54.1289,
-56.2226,
-54.9853,
-54.5406,
-55.8668,
-54.4344,
-51.34,
-53.4457,
-55.3933,
-56.4022,
-57.494,
-57.042,
-53.8239,
-54.7248,
-55.1078,
-54.9422,
-57.6964,
-57.2593,
-54.7605,
-56.342,
-57.4363,
-53.8504,
-52.5132,
-54.1004,
-55.4099,
-55.062,
-54.2594,
-52.8564,
-52.14,
-51.7081,
-52.2992,
-52.3724,
-51.8007,
-55.592,
-55.7873,
-53.5016,
-54.1608,
-53.7607,
-52.8233,
-54.0887,
-54.6511,
-54.4701,
-54.7901,
-56.7217,
-55.1668,
-54.6029,
-55.2335,
-53.67,
-53.9694,
-56.1003,
-55.6258,
-55.8682,
-54.7942,
-54.6461,
-56.4503,
-57.3815,
-55.4552,
-56.4655,
-55.3854,
-54.1829,
-53.3174,
-54.3715,
-53.5423,
-54.4937,
-56.8845,
-54.4562,
-53.0783,
-54.8609,
-52.7257,
-53.1482,
-55.1311,
-54.7786,
-54.3794,
-55.2594,
-52.3897,
-50.3529,
-50.6989,
-50.8013,
-50.2843,
-51.4467,
-52.3954,
-51.4057,
-51.6931,
-53.6986,
-52.1103,
-49.9167,
-51.4783,
-53.218,
-53.8505,
-52.805,
-51.589,
-51.8991,
-53.1859,
-50.7663,
-51.6103,
-52.6432,
-50.0238,
-52.5902,
-54.4426,
-51.1014,
-51.337,
-52.8024,
-53.7283,
-53.3006,
-54.6558,
-56.1147,
-53.3179,
-52.8044,
-52.334,
-51.9468,
-51.2366,
-56.9915,
-54.5127,
-53.0841,
-54.4758,
-53.2447,
-54.7147,
-54.3793,
-52.4646,
-52.7251,
-53.2872,
-51.1914,
-52.1654,
-53.3249,
-52.2418,
-50.3992,
-51.7577,
-50.617,
-50.6632,
-54.7326,
-52.635,
-51.3294,
-54.1903,
-53.3178,
-51.235,
-53.1232,
-52.5514,
-51.5221,
-49.9557,
-52.2744,
-53.2325,
-51.3947,
-51.9394,
-52.0155,
-51.9813,
-52.9384,
-51.6752,
-51.4657,
-53.9598,
-56.0678,
-55.6356,
-54.9773,
-52.1095,
-49.7851,
-51.2385,
-52.6269,
-53.3314,
-55.0205,
-55.7239,
-53.4701,
-52.5249,
-53.1064,
-55.33,
-53.1046,
-52.853,
-53.7369,
-54.7797,
-54.5858,
-54.2175,
-53.2216,
-50.8936,
-36.8913,
-35.1439,
-37.0516,
-50.6355,
-52.1987,
-53.0451,
-53.1897,
-52.8646,
-52.6694,
-52.9763,
-53.23,
-54.3301,
-55.2726,
-54.0729,
-51.3799,
-50.7924,
-51.3911,
-51.1238,
-50.1222,
-51.9869,
-51.6642,
-50.5145,
-50.098,
-49.673,
-51.0346,
-48.4426,
-46.532,
-49.659,
-49.8172,
-47.0652,
-48.479,
-50.3125,
-50.6249,
-52.0916,
-49.325,
-46.4799,
-48.3024,
-51.8701,
-48.8237,
-50.4471,
-50.5064,
-49.4765,
-51.0378,
-49.8951,
-50.867,
-51.7528,
-49.7907,
-50.9414,
-50.196,
-50.7166,
-48.2638,
-48.1027,
-49.4721,
-51.5115,
-49.6891,
-49.1679,
-50.4271,
-50.3478,
-50.5238,
-49.163,
-48.1769,
-48.8715,
-51.5548,
-48.3888,
-47.5323,
-50.1061,
-49.1536,
-49.2668,
-49.7307,
-51.1017,
-54.1429,
-50.4325,
-49.1318,
-48.6643,
-50.4365,
-49.6167,
-48.396,
-49.2512,
-50.9879,
-49.5467,
-50.9555,
-53.2533,
-50.8848,
-50.4579,
-50.1226,
-49.8508,
-49.3174,
-50.3957,
-47.6939,
-47.1738,
-49.8836,
-49.2091,
-49.4535,
-53.6354,
-52.5986,
-51.6961,
-52.7263,
-48.666,
-49.2609,
-52.923,
-52.6971,
-51.1352,
-50.2645,
-48.4548,
-48.1124,
-47.8769,
-48.4172,
-46.8714,
-48.2583,
-50.2179,
-48.3703,
-49.6104,
-49.8433,
-46.8929,
-47.5154,
-50.6053,
-51.3085,
-51.9856,
-50.2469,
-47.5982,
-49.3081,
-51.3249,
-48.7754,
-49.3255,
-50.9454,
-49.4825,
-50.0666,
-49.8912,
-48.8952,
-48.1874,
-49.0473,
-48.5272,
-47.6251,
-50.1978,
-51.7487,
-48.6737,
-49.9187,
-50.7832,
-48.591,
-47.5307,
-53.948,
-49.7888,
-47.7759,
-49.4646,
-49.2951,
-48.0835,
-50.9995,
-47.7416,
-47.4029,
-51.7832,
-47.9682,
-46.4668,
-49.6323,
-50.7472,
-48.5351,
-48.8773,
-49.049,
-48.3435,
-50.9687,
-49.1747,
-46.6598,
-48.7942,
-49.6085,
-47.7136,
-46.3717,
-47.9584,
-48.8272,
-47.049,
-48.0484,
-48.5147,
-47.7235,
-48.9985,
-48.3041,
-48.7325,
-52.1014,
-48.6559,
-45.9734,
-47.9816,
-49.9001,
-49.744,
-49.0705,
-49.5089,
-50.6053,
-51.4918,
-50.2007,
-49.1247,
-50.8951,
-53.99,
-50.9961,
-51.2368,
-54.649,
-50.5483,
-49.6662,
-49.7394,
-48.4737,
-50.7644,
-52.0322,
-51.6659,
-48.8891,
-50.7512,
-51.6192,
-49.0519,
-48.6193,
-49.1657,
-49.9413,
-50.675,
-50.7688,
-47.5235,
-48.7736,
-51.2266,
-50.0237,
-48.7437,
-51.924,
-52.801,
-49.2737,
-46.3321,
-48.3736,
-48.1676,
-46.3604,
-48.1548,
-51.7357,
-48.9502,
-48.1962,
-48.5177,
-49.363,
-48.0272,
-45.426,
-48.0143,
-48.0975,
-45.1166,
-46.3444,
-48.1079,
-47.5219,
-47.5311,
-47.5127,
-45.4789,
-47.1243,
-48.7736,
-46.554,
-48.3435,
-50.5012,
-46.9389,
-47.8544,
-48.8036,
-50.0142,
-50.8305,
-51.3919,
-50.5138,
-47.5832,
-47.4375,
-47.9406,
-48.4136,
-47.9574,
-46.4125,
-46.3805,
-49.6796,
-48.7645,
-46.8555,
-48.4917,
-48.5139,
-45.8423,
-47.3436,
-49.4883,
-47.0694,
-46.6695,
-48.0118,
-47.1591,
-49.1952,
-50.3417,
-48.5627,
-47.191,
-48.3246,
-48.2971,
-47.0113,
-48.5018,
-49.5539,
-48.1113,
-48.9499,
-48.0196,
-44.7132,
-46.8052,
-49.3817,
-48.9602,
-51.7208,
-47.954,
-45.4842,
-49.5872,
-47.2174,
-44.829,
-46.5841,
-48.7515,
-47.0837,
-47.6698,
-49.6554,
-49.0501,
-48.0787,
-47.2196,
-48.6574,
-46.4617,
-47.7125,
-47.2467,
-49.2226,
-50.1231,
-47.5745,
-45.2543,
-46.22,
-46.6847,
-45.2459,
-45.3651,
-48.327,
-45.4815,
-45.0398,
-45.7934,
-45.2835,
-44.7138,
-45.9321,
-45.2813,
-47.0122,
-46.104,
-45.367,
-45.9876,
-47.9313,
-48.014,
-45.8913,
-45.8209,
-44.539,
-43.8343,
-47.3054,
-48.6662,
-47.7547,
-46.8956,
-48.9327,
-47.4592,
-48.1918,
-47.9374,
-45.3557,
-45.4929,
-45.3678,
-43.5012,
-43.1875,
-45.2978,
-46.5465,
-49.1348,
-49.1225,
-46.0337,
-46.0285,
-47.0877,
-44.4196,
-44.1308,
-46.6495,
-45.2522,
-44.6563,
-45.8002,
-44.6618,
-44.5779,
-44.4855,
-44.2913,
-43.1099,
-43.9997,
-45.2453,
-45.0873,
-46.2542,
-47.2774,
-45.2886,
-44.5012,
-45.9774,
-44.1381,
-45.4811,
-48.5527,
-47.2975,
-44.8606,
-46.5022,
-46.7482,
-46.028,
-48.0085,
-47.6895,
-45.0938,
-46.893,
-47.9799,
-46.101,
-46.2139,
-48.2228,
-47.5895,
-45.3641,
-45.9702,
-45.4339,
-47.2054,
-50.1577,
-47.9846,
-48.8655,
-48.461,
-48.6598,
-48.6151,
-49.7885,
-47.47,
-49.6225,
-47.6268,
-47.3095,
-48.6953,
-47.3902,
-47.3501,
-46.8645,
-49.2446,
-46.3305,
-44.4973,
-46.9743,
-48.1871,
-48.1099,
-46.9727,
-45.9809,
-47.0705,
-48.4405,
-45.9353,
-45.4907,
-46.1851,
-44.4699,
-47.7013,
-46.1764,
-46.9609,
-46.7163,
-45.2347,
-45.2775,
-45.7384,
-45.4096,
-46.1625,
-43.9906,
-42.2237,
-43.7476,
-44.4925,
-42.0441,
-43.5432,
-44.4195,
-44.1445,
-45.2037,
-43.0051,
-41.7508,
-43.4584,
-45.7652,
-43.0487,
-44.2161,
-45.7777,
-44.9954,
-45.5829,
-45.2463,
-46.069,
-46.7692,
-44.4243,
-45.4386,
-45.8462,
-45.5292,
-42.6122,
-42.72,
-45.0977,
-47.1575,
-41.8999,
-42.8585,
-44.6372,
-45.2421,
-45.8997,
-43.8037,
-43.0463,
-43.7606,
-49.1001,
-43.6061,
-42.4586,
-43.3683,
-43.391,
-44.8932,
-46.0093,
-46.4953,
-46.0914,
-45.0499,
-44.6887,
-44.5065,
-44.6259,
-43.2564,
-41.6824,
-44.2345,
-46.3672,
-44.4248,
-44.7575,
-45.4479,
-46.4974,
-46.4438,
-48.1768,
-47.8746,
-48.1585,
-46.0249,
-44.2233,
-45.0271,
-47.7074,
-45.1241,
-45.602,
-48.034,
-45.7877,
-45.7414,
-48.8029,
-46.1515,
-44.7503,
-46.0245,
-48.1391,
-42.7574,
-45.0704,
-44.4313,
-42.3392,
-42.7127,
-42.63,
-42.1699,
-43.1941,
-44.8777,
-43.7871,
-44.7024,
-43.1134,
-41.052,
-41.5106,
-42.7823,
-43.9133,
-46.6444,
-44.1329,
-42.4254,
-41.9902,
-43.1927,
-42.4992,
-43.2412,
-45.32,
-43.1171,
-41.8517,
-45.7062,
-44.5682,
-43.9681,
-42.4479,
-42.7704,
-45.2259,
-46.2546,
-43.2546,
-42.3056,
-42.9318,
-41.2086,
-40.1974,
-38.86,
-41.7293,
-45.3347,
-43.45,
-43.1411,
-42.5701,
-42.2307,
-42.1508,
-40.1295,
-39.5435,
-44.5262,
-44.3852,
-41.6828,
-42.4575,
-41.5466,
-41.2296,
-40.6723,
-40.7443,
-41.0065,
-44.0477,
-44.4363,
-41.7784,
-43.8286,
-44.3334,
-40.4965,
-40.7184,
-42.1522,
-40.372,
-40.0213,
-42.1974,
-43.9423,
-41.4528,
-41.3604,
-40.0896,
-40.5994,
-43.1084,
-44.5182,
-40.5199,
-43.2448,
-43.9948,
-42.4505,
-43.8804,
-41.658,
-41.7391,
-43.9621,
-42.5052,
-41.941,
-43.953,
-44.4372,
-42.9089,
-41.4266,
-44.5048,
-44.0139,
-43.1635,
-42.9775,
-42.5989,
-44.2387,
-43.6555,
-42.2137,
-40.8761,
-41.2583,
-41.6775,
-39.2895,
-38.7132,
-39.5674,
-41.6928,
-38.4184,
-36.5258,
-36.6337,
-36.45,
-35.4038,
-34.3154,
-33.8786,
-34.3,
-32.2024,
-31.6128,
-29.9677,
-28.604,
-26.7203,
-21.491,
-20.0435,
-20.0388,
-19.4482,
-19.675,
-15.5392,
-15.7468,
-15.7612,
-14.1639,
-17.6221,
-20.3647,
-15.2342,
-14.4646,
-15.9993,
-16.7893,
-20.0614,
-18.6413,
-15.156,
-15.7714,
-18.9955,
-11.9265,
-13.1928,
-17.4033,
-15.9239,
-15.6211,
-15.1626,
-17.2121,
-15.1454,
-14.846,
-17.4043,
-14.9008,
-17.8761,
-13.6937,
-12.7696,
-16.4788,
-17.8299,
-14.1835,
-15.9506,
-15.5977,
-16.4307,
-16.6612,
-17.2893,
-17.1491,
-15.2785,
-14.299,
-11.5101,
-12.0057,
-14.3695,
-12.3953,
-14.5376,
-18.6463,
-17.6179,
-17.2854,
-14.838,
-15.901,
-15.4921,
-16.3385,
-14.7699,
-11.5637,
-13.0916,
-16.2336,
-16.4366,
-15.3438,
-16.4631,
-16.434,
-15.4712,
-14.356,
-16.113,
-16.3932,
-14.813,
-10.4706,
-14.0698,
-18.1615,
-16.0609,
-15.0685,
-18.28,
-15.6291,
-15.2425,
-18.0123,
-16.8178,
-13.9477,
-18.1202,
-20.0955,
-18.5818,
-18.6161,
-17.0272,
-15.1053,
-14.8645,
-14.2902,
-14.2992,
-18.4018,
-18.3067,
-14.4854,
-17.361,
-16.466,
-12.033,
-12.1618,
-18.0388,
-14.6222,
-14.4357,
-16.3651,
-14.5199,
-16.5932,
-18.0895,
-18.2942,
-14.4461,
-14.9826,
-13.7885,
-11.4138,
-14.7876,
-17.2821,
-15.1468,
-14.2192,
-14.4969,
-14.6106,
-19.7936,
-18.5169,
-15.4286,
-17.5611,
-17.1634,
-13.5942,
-13.6943,
-16.9909,
-16.9429,
-16.4109,
-18.8415,
-16.5409,
-15.0603,
-15.8583,
-15.0508,
-14.7035,
-20.1458,
-14.2932,
-11.0877,
-14.5908,
-18.6891,
-16.1547,
-15.6604,
-17.4981,
-15.1965,
-16.2621,
-15.3162,
-15.8825,
-18.2088,
-17.8679,
-12.9174,
-13.0332,
-14.5191,
-15.1663,
-17.884,
-19.5843,
-16.3794,
-15.5378,
-15.8095,
-16.8979,
-16.1351,
-16.554,
-14.7715,
-11.7863,
-15.7083,
-16.0304,
-16.8404,
-19.9122,
-15.3532,
-17.1025,
-16.169,
-18.5948,
-18.5845,
-17.9948,
-20.0339,
-16.0732,
-15.9746,
-18.4749,
-17.5057,
-15.7567,
-18.9827,
-18.5605,
-19.4898,
-19.8482,
-19.9323,
-19.4093,
-20.9551,
-19.3766,
-18.472,
-16.8604,
-15.9094,
-16.6785,
-18.2195,
-20.4397,
-17.7038,
-16.4312,
-20.627,
-20.5637,
-17.8994,
-19.3752,
-15.7811,
-15.8661,
-21.2333,
-17.3369,
-17.6729,
-17.944,
-16.6844,
-15.3104,
-18.4577,
-17.4173,
-14.2345,
-16.6316,
-15.1417,
-12.4047,
-16.2736,
-18.1997,
-13.7877,
-14.0002,
-15.3601,
-15.24,
-16.6319,
-17.412,
-15.6287,
-16.9564,
-18.5734,
-12.4446,
-13.2975,
-16.535,
-16.7369,
-16.7826,
-16.7616,
-17.2345,
-15.6208,
-15.6083,
-14.7413,
-15.5053,
-19.3031,
-14.9355,
-12.4761,
-16.9753,
-19.1808,
-14.589,
-15.7518,
-17.9452,
-17.4412,
-18.4867,
-19.3394,
-15.8374,
-15.3029,
-16.1687,
-13.7927,
-12.9113,
-14.1068,
-14.4755,
-16.5459,
-20.7771,
-18.5065,
-16.205,
-14.6572,
-14.6384,
-15.8938,
-17.6287,
-14.2392,
-12.21,
-14.5071,
-18.1666,
-16.0671,
-14.9861,
-16.4561,
-18.072,
-15.2579,
-13.6526,
-14.9221,
-15.7395,
-13.5946,
-10.9806,
-14.4211,
-17.8838,
-15.2315,
-16.0815,
-17.462,
-16.6821,
-15.6187,
-19.4224,
-14.3781,
-13.0367,
-18.8141,
-19.094,
-17.2449,
-19.502,
-18.4997,
-16.9683,
-18.2979,
-16.559,
-14.2958,
-17.2245,
-16.5611,
-15.965,
-16.2773,
-16.0848,
-13.4482,
-13.628,
-17.3424,
-15.2759,
-15.4864,
-15.5957,
-15.4528,
-17.081,
-18.0077,
-16.725,
-16.0681,
-16.6705,
-14.7472,
-11.1078,
-15.5698,
-18.9295,
-14.5944,
-15.3533,
-14.4935,
-14.1263,
-18.3201,
-16.4641,
-14.0456,
-16.3689,
-15.9044,
-13.9083,
-15.1376,
-18.0881,
-15.4433,
-18.0392,
-19.8115,
-17.2612,
-16.1393,
-14.6346,
-12.639,
-13.3451,
-20.3823,
-15.5033,
-11.6252,
-17.2243,
-19.6212,
-15.889,
-16.2187,
-17.2933,
-15.3106,
-15.2818,
-15.3002,
-15.3482,
-16.0188,
-15.113,
-13.0713,
-12.5087,
-15.7176,
-17.3474,
-16.7773,
-19.4358,
-16.6256,
-15.8939,
-17.3975,
-16.3073,
-15.8907,
-17.3589,
-15.003,
-11.9169,
-16.0895,
-14.965,
-15.6599,
-18.6817,
-17.8752,
-18.2858,
-17.1836,
-18.8071,
-17.1642,
-16.1271,
-16.7472,
-19.1897,
-18.4565,
-20.2236,
-18.0552,
-18.4603,
-20.0452,
-18.7502,
-17.8809,
-22.238,
-29.1337,
-30.655,
-30.9884,
-31.0855,
-26.6656,
-30.7148,
-31.103,
-29.2949,
-26.9712,
-21.8321,
-17.1371,
-19.202,
-20.7523,
-18.956,
-18.5722,
-22.2649,
-17.2581,
-17.3095,
-17.8291,
-16.7546,
-18.3583,
-19.4196,
-17.1012,
-19.0665,
-19.6574,
-18.7763,
-15.6548,
-17.41,
-14.4814,
-12.649,
-16.4868,
-19.6048,
-15.9429,
-14.6092,
-16.2591,
-16.6707,
-16.7363,
-17.5889,
-15.5971,
-17.7072,
-16.9493,
-12.379,
-14.5237,
-17.3588,
-18.0528,
-18.7503,
-17.4407,
-17.773,
-16.1955,
-16.316,
-17.5686,
-16.1528,
-17.5352,
-15.1959,
-12.5491,
-16.2168,
-17.951,
-16.2831,
-16.8028,
-21.7191,
-17.5696,
-18.5079,
-17.8364,
-16.8529,
-17.3568,
-16.4887,
-14.5074,
-13.6032,
-15.0561,
-14.8592,
-17.1061,
-18.2881,
-17.3267,
-17.4816,
-15.2771,
-15.9633,
-17.1671,
-16.8132,
-15.4977,
-13.4716,
-14.6176,
-15.4932,
-15.5906,
-16.0732,
-16.6312,
-19.8594,
-17.7166,
-14.7387,
-16.1423,
-17.7272,
-16.0497,
-12.0825,
-14.8194,
-18.6578,
-16.2647,
-16.952,
-21.209,
-17.2636,
-15.4765,
-21.519,
-18.9253,
-15.6929,
-18.4423,
-22.0175,
-20.0796,
-19.7565,
-18.4389,
-17.4147,
-16.1838,
-15.7579,
-15.7313,
-17.9307,
-17.6456,
-15.472,
-16.2165,
-18.3629,
-13.9219,
-13.3908,
-18.0883,
-15.8902,
-15.8117,
-17.6531,
-15.399,
-17.1985,
-20.5278,
-18.2017,
-16.0071,
-16.6188,
-15.3262,
-12.4708,
-16.0418,
-18.4806,
-15.2803,
-14.1438,
-15.4045,
-17.4351,
-21.2766,
-18.956,
-14.8291,
-15.6915,
-17.9646,
-13.796,
-14.2603,
-17.9241,
-17.1372,
-20.2816,
-18.5349,
-16.7273,
-15.9286,
-17.0292,
-15.6254,
-16.2442,
-21.4171,
-15.2941,
-11.6143,
-16.4831,
-19.1565,
-17.3051,
-16.2204,
-16.7853,
-14.3267,
-16.5206,
-17.4987,
-15.8609,
-16.1717,
-18.7535,
-13.7566,
-13.515,
-16.0579,
-16.2419,
-17.9404,
-20.2519,
-18.8312,
-16.4888,
-16.1905,
-18.4802,
-15.4975,
-16.8566,
-15.1004,
-12.0728,
-16.1559,
-16.5879,
-17.5269,
-17.6385,
-17.0088,
-18.2969,
-17.6654,
-18.2338,
-18.5142,
-18.5781,
-23.0018,
-18.7451,
-18.7929,
-20.3262,
-18.2918,
-20.992,
-20.8563,
-21.0641,
-19.8925,
-20.8524,
-19.6837,
-20.3027,
-20.8755,
-18.7651,
-19.745,
-19.0357,
-22.9735,
-19.6601,
-17.5248,
-22.4799,
-22.2368,
-20.741,
-21.4475,
-19.1367,
-20.0864,
-20.6374,
-17.0022,
-17.0454,
-20.2269,
-20.2276,
-18.1309,
-18.3089,
-19.0382,
-17.7965,
-17.7277,
-19.0534,
-16.6147,
-16.9542,
-15.08,
-13.8843,
-18.0411,
-20.2552,
-17.1335,
-15.8885,
-16.4164,
-19.8657,
-21.2539,
-21.6384,
-17.9591,
-16.4973,
-19.4645,
-13.7794,
-14.5086,
-17.499,
-18.7878,
-22.569,
-18.4554,
-19.112,
-17.7439,
-16.5406,
-19.7502,
-17.2812,
-19.5686,
-16.5935,
-14.6281,
-16.7502,
-19.8096,
-16.3894,
-17.3978,
-17.9746,
-16.3131,
-16.844,
-18.1621,
-16.3371,
-16.2464,
-18.0401,
-14.3329,
-13.5673,
-16.2798,
-15.4207,
-17.4112,
-20.4657,
-19.6506,
-18.5436,
-18.5071,
-18.5784,
-18.1828,
-18.302,
-15.7905,
-14.5212,
-15.8785,
-17.87,
-17.5357,
-17.2979,
-19.5085,
-15.4981,
-16.5418,
-16.4873,
-16.6953,
-18.9113,
-17.337,
-12.8396,
-15.5374,
-18.1843,
-18.6142,
-17.9819,
-20.3665,
-18.874,
-17.8146,
-21.8009,
-19.4206,
-15.1749,
-16.941,
-21.3388,
-19.2361,
-18.5137,
-19.5836,
-17.8294,
-18.0813,
-16.2072,
-15.2342,
-17.3446,
-17.9416,
-17.0639,
-19.0464,
-18.1247,
-14.0078,
-14.1669,
-17.2219,
-16.7482,
-15.523,
-17.009,
-17.2274,
-17.2415,
-20.7098,
-22.7188,
-15.6795,
-16.1834,
-15.076,
-12.8413,
-17.1034,
-21.6745,
-16.3876,
-15.0064,
-15.6645,
-15.5831,
-19.8088,
-18.1582,
-15.6864,
-16.8537,
-16.605,
-14.3675,
-14.9494,
-18.2275,
-17.9394,
-17.9904,
-21.1667,
-18.6842,
-17.2572,
-16.2678,
-15.2376,
-15.7674,
-20.0147,
-15.5518,
-12.6766,
-16.4039,
-15.9753,
-16.8163,
-16.7973,
-17.6159,
-17.4915,
-17.6545,
-17.1689,
-16.1077,
-18.1998,
-18.469,
-13.9715,
-13.9777,
-16.6906,
-16.3311,
-17.1617,
-19.8228,
-17.2765,
-18.7073,
-20.429,
-19.0216,
-16.8083,
-20.0403,
-18.0369,
-14.1958,
-18.8201,
-19.6189,
-18.9291,
-21.8515,
-21.4555,
-19.4145,
-21.1048,
-24.6117,
-30.1201,
-30.2087,
-32.4339,
-34.0687,
-34.3296,
-35.9902,
-35.9791,
-36.9887,
-38.0751,
-37.0478,
-38.1653,
-39.9887,
-41.1763,
-40.8435,
-40.5692,
-42.2582,
-43.4551,
-43.1281,
-41.4753,
-41.3605,
-44.7594,
-48.0518,
-45.8695,
-44.5699,
-44.9164,
-46.7043,
-46.9795,
-45.987,
-46.8277,
-46.9546,
-48.517,
-48.9181,
-48.5357,
-47.1731,
-46.1186,
-46.5777,
-49.208,
-50.2139,
-47.5279,
-47.6892,
-45.6928,
-42.3576,
-45.6602,
-50.0009,
-44.0409,
-44.8954,
-46.0335,
-46.3665,
-47.0359,
-46.5136,
-43.8663,
-45.1994,
-48.9513,
-43.6399,
-44.5852,
-48.9535,
-48.7395,
-50.6469,
-47.9366,
-50.2694,
-53.1197,
-51.219,
-48.9197,
-49.6145,
-49.3445,
-46.2711,
-46.8294,
-48.7328,
-48.689,
-46.8953,
-46.455,
-47.4365,
-49.5511,
-47.2848,
-48.2091,
-48.48,
-49.6225,
-49.5303,
-45.9227,
-45.6378,
-47.6079,
-48.3332,
-51.1252,
-50.0722,
-47.7211,
-46.3724,
-47.9191,
-46.9733,
-46.4178,
-47.9862,
-48.7368,
-47.3079,
-48.573,
-47.8612,
-47.7978,
-47.9014,
-47.465,
-48.1392,
-51.7731,
-50.9328,
-48.3029,
-48.4842,
-51.4898,
-48.723,
-48.3427,
-50.284,
-50.1632,
-49.0694,
-49.4466,
-49.4901,
-51.7125,
-51.408,
-49.0947,
-46.7358,
-48.8949,
-50.5354,
-51.6669,
-48.0811,
-46.9125,
-47.2465,
-48.6119,
-48.0179,
-47.1497,
-49.5535,
-49.6195,
-51.0516,
-50.4786,
-49.7187,
-47.9761,
-49.0253,
-50.6974,
-49.7176,
-48.4694,
-49.1365,
-48.0496,
-48.3849,
-49.6228,
-48.9914,
-46.2694,
-47.6119,
-49.1836,
-50.8548,
-50.1212,
-49.1347,
-47.36,
-49.579,
-47.6367,
-48.0236,
-50.94,
-51.1781,
-50.3151,
-47.9266,
-49.14,
-48.8484,
-47.8906,
-51.3558,
-51.0253,
-51.8225,
-52.3138,
-50.3064,
-50.2402,
-50.6688,
-47.7897,
-48.0829,
-50.0277,
-49.8401,
-46.2139,
-48.474,
-48.5406,
-46.5788,
-46.0183,
-48.064,
-47.0858,
-48.4572,
-51.3744,
-49.071,
-48.7529,
-49.3962,
-47.3005,
-47.7479,
-49.1209,
-45.7782,
-46.6839,
-49.1914,
-48.0888,
-48.3552,
-48.7767,
-50.5968,
-50.9137,
-49.9867,
-47.3503,
-46.2074,
-50.2556,
-49.8635,
-46.9128,
-48.4134,
-50.8621,
-48.2746,
-47.2561,
-48.3277,
-48.462,
-49.5798,
-51.1221,
-49.8699,
-49.8763,
-52.7402,
-51.8437,
-50.9413,
-49.4195,
-48.6313,
-48.7089,
-48.26,
-49.1006,
-51.7097,
-51.2647,
-50.0222,
-46.266,
-47.6016,
-51.4242,
-51.2863,
-50.4255,
-50.777,
-50.6262,
-52.4123,
-51.1417,
-48.6119,
-49.1773,
-49.7257,
-51.0008,
-51.1137,
-49.4485,
-47.3073,
-46.4001,
-49.1763,
-49.6873,
-48.8188,
-47.7588,
-48.8807,
-49.5315,
-50.5,
-49.1801,
-46.205,
-49.0411,
-48.2937,
-47.7029,
-49.864,
-49.8018,
-51.6115,
-50.3661,
-49.6347,
-48.2495,
-51.1527,
-50.0531,
-50.046,
-51.5537,
-52.1198,
-52.4282,
-51.7337,
-50.2507,
-50.2803,
-49.4137,
-49.3984,
-50.0768,
-48.9663,
-49.2549,
-49.3198,
-48.9314,
-51.3619,
-50.9843,
-49.4912,
-50.9405,
-48.5417,
-46.2278,
-48.7857,
-49.9857,
-48.5979,
-49.3343,
-49.8103,
-51.2711,
-50.6175,
-52.9629,
-46.3871,
-47.3597,
-50.0049,
-51.3863,
-50.9092,
-49.2139,
-49.7167,
-48.2345,
-49.7085,
-47.9198,
-47.2056,
-48.779,
-49.5685,
-47.9345,
-48.9486,
-48.4122,
-46.8779,
-48.9406,
-47.3671,
-49.5441,
-50.6269,
-48.1014,
-47.5496,
-50.8901,
-50.5588,
-50.3705,
-49.6213,
-50.9238,
-51.3464,
-51.3144,
-51.2102,
-51.0348,
-48.9954,
-50.0406,
-51.2378,
-48.9083,
-50.871,
-49.1628,
-49.3872,
-47.8336,
-45.5733,
-45.3397,
-48.3252,
-49.8814,
-49.5617,
-52.5514,
-50.8642,
-47.4364,
-46.9441,
-49.2649,
-48.6666,
-48.7245,
-52.2812,
-50.4663,
-48.9467,
-50.4195,
-49.3005,
-48.0456,
-49.7586,
-51.4669,
-49.4233,
-49.3954,
-51.918,
-52.5363,
-53.245,
-50.1197,
-49.5715,
-48.6225,
-51.4814,
-48.6501,
-48.435,
-49.0636,
-46.99,
-47.1601,
-51.3511,
-51.7712,
-50.9188,
-49.9227,
-47.5559,
-47.5108,
-47.5403,
-48.5864,
-48.7457,
-52.6943,
-47.607,
-45.5889,
-46.5265,
-47.079,
-46.4308,
-46.6699,
-46.537,
-46.3403,
-50.2296,
-49.3748,
-48.8292,
-48.3511,
-49.7581,
-48.8874,
-49.5481,
-51.5764,
-53.0174,
-51.1775,
-49.355,
-48.1056,
-50.0939,
-52.3718,
-49.6871,
-48.6252,
-48.5005,
-49.1256,
-47.7784,
-49.501,
-49.6147,
-50.4504,
-52.848,
-52.5858,
-51.5435,
-50.0252,
-50.4699,
-49.7307,
-50.7538,
-52.7033,
-51.7882,
-52.9639,
-53.4579,
-50.473,
-50.4015,
-52.9672,
-52.2223,
-51.081,
-52.4087,
-53.6466,
-53.7278,
-51.9755,
-50.2002,
-48.3145,
-50.6177,
-51.2229,
-48.8809,
-52.0095,
-53.5939,
-50.5413,
-52.9807,
-53.6694,
-53.6789,
-54.1299,
-51.5615,
-51.0408,
-54.0533,
-51.3977,
-49.306,
-51.6357,
-53.7669,
-55.5275,
-54.4678,
-51.1255,
-50.2686,
-48.2089,
-50.3321,
-48.6478,
-47.8024,
-51.0462,
-51.0639,
-49.8515,
-52.1705,
-50.1129,
-49.5335,
-51.425,
-52.8038,
-50.2589,
-50.7481,
-48.5987,
-47.1937,
-48.4121,
-50.5755,
-51.7702,
-50.9743,
-51.1422,
-50.6376,
-51.6163,
-50.8467,
-52.6728,
-50.7194,
-50.9636,
-53.1127,
-53.6698,
-51.9594,
-49.8107,
-48.4892,
-46.9902,
-49.1573,
-49.527,
-49.2974,
-53.1054,
-51.2926,
-50.9592,
-50.8261,
-48.1979,
-48.3261,
-49.4646,
-48.531,
-49.1327,
-50.5197,
-50.4211,
-52.3129,
-52.3579,
-49.2716,
-47.8853,
-48.2628,
-46.8968,
-47.9587,
-50.7874,
-50.925,
-49.1492,
-52.0662,
-51.2057,
-49.733,
-49.6531,
-50.086,
-51.8516,
-52.064,
-49.5581,
-48.727,
-52.7099,
-51.1081,
-48.84,
-51.0247,
-54.5521,
-55.9998,
-51.3357,
-51.0436,
-53.4228,
-51.0122,
-50.6261,
-53.7108,
-54.3617,
-52.8832,
-50.9914,
-48.9337,
-48.592,
-49.4738,
-50.0585,
-50.9642,
-52.1417,
-52.0936,
-52.8795,
-52.6958,
-49.9302,
-47.6595,
-52.1132,
-51.3054,
-50.0328,
-50.9853,
-50.9056,
-52.5706,
-52.3182,
-49.0835,
-48.2247,
-52.7862,
-51.8768,
-49.5388,
-51.7117,
-53.3644,
-52.5562,
-53.1913,
-53.2061,
-52.064,
-51.4112,
-52.2875,
-51.7046,
-53.4237,
-51.6789,
-48.9602,
-49.615,
-54.3065,
-52.7595,
-51.7539,
-52.4489,
-51.0018,
-48.8141,
-49.9922,
-50.1697,
-51.5571,
-55.0326,
-49.2505,
-47.7099,
-50.4369,
-49.4718,
-47.7357,
-49.81,
-49.7595,
-49.5549,
-52.5567,
-50.5915,
-51.071,
-50.3581,
-48.8385,
-47.9719,
-50.531,
-53.754,
-54.5275,
-51.6008,
-50.5125,
-50.2016,
-51.6389,
-52.0413,
-52.9119,
-52.7238,
-51.7382,
-51.692,
-50.3941,
-52.276,
-53.055,
-51.7107,
-51.5907,
-56.1152,
-54.4785,
-53.1961,
-52.6983,
-51.7241,
-55.7689,
-51.4337,
-52.7248,
-53.1615,
-56.9085,
-52.0995,
-54.9093,
-56.176,
-55.7428,
-55.3085,
-55.296,
-55.3061,
-54.5521,
-38.7349,
-36.6404,
-35.6112,
-36.3573,
-38.6825,
-40.9994,
-54.838,
-55.4042,
-55.9963,
-56.5734,
-55.3016,
-52.379,
-54.4772,
-57.5687,
-55.7754,
-56.1422,
-58.4314,
-56.3081,
-58.3597,
-54.8967,
-55.1681,
-54.94,
-56.9699,
-54.5725,
-53.3058,
-55.4576,
-57.0444,
-54.5734,
-55.4882,
-58.8719,
-57.8147,
-57.161,
-55.9259,
-54.1616,
-54.9256,
-57.9064,
-55.5892,
-54.4938,
-56.6013,
-53.3273,
-53.5378,
-56.8796,
-54.122,
-54.2593,
-56.193,
-53.9324,
-55.5943,
-57.948,
-55.9327,
-55.7052,
-54.8921,
-55.0404,
-54.565,
-55.2447,
-53.395,
-54.9005,
-56.3495,
-57.184,
-56.2635,
-57.2707,
-57.4239,
-56.9224,
-56.0371,
-56.0452,
-57.16,
-58.3847,
-58.6806,
-56.2668,
-58.4003,
-56.7104,
-56.2934,
-57.4744,
-56.0358,
-54.3808,
-54.3828,
-54.9765,
-55.1621,
-55.7418,
-55.4705,
-58.0448,
-56.9212,
-57.1734,
-55.65,
-54.5887,
-57.506,
-60.9133,
-57.8444,
-57.463,
-60.6389,
-57.9056,
-56.315,
-57.7275,
-59.5114,
-58.9604,
-57.4103,
-56.0218,
-55.2313,
-55.5317,
-57.8286,
-55.9575,
-57.8296,
-58.639,
-55.4288,
-55.5388,
-54.5277,
-55.7242,
-56.5171,
-57.9004,
-58.5774,
-58.0525,
-57.3392,
-56.6642,
-57.3167,
-57.4948,
-57.6944,
-58.4729,
-56.7523,
-56.0368,
-56.9025,
-56.2476,
-55.9178,
-56.8256,
-57.2025,
-56.0975,
-56.4654,
-55.9343,
-58.3103,
-56.9723,
-56.2622,
-58.8845,
-57.8157,
-58.9594,
-58.1085,
-57.076,
-58.3772,
-61.9133,
-59.8652,
-58.0047,
-60.263,
-58.0358,
-58.024,
-61.0247,
-59.0908,
-58.1108,
-58.3913,
-57.9262,
-58.2836,
-59.0939,
-58.8803,
-58.1376,
-57.91,
-58.8683,
-58.4782,
-59.3043,
-58.0344,
-57.5405,
-57.2232,
-56.4142,
-55.2308,
-58.6367,
-60.2894,
-58.3026,
-57.6476,
-58.7072,
-57.2483,
-57.0139,
-58.3794,
-58.8506,
-58.9216,
-61.2998,
-60.9501,
-60.7788,
-57.7742,
-56.8946,
-56.6763,
-58.6238,
-59.6427,
-58.9964,
-60.5465,
-58.0314,
-57.6388,
-59.3954,
-58.3974,
-57.6562,
-58.0348,
-58.047,
-60.6587,
-60.8311,
-57.935,
-58.2454,
-58.5662,
-58.1951,
-58.0662,
-59.569,
-59.0606,
-58.507,
-57.3905,
-58.0025,
-59.7966,
-58.2268,
-58.0629,
-57.8326,
-57.1805,
-59.3286,
-59.5327,
-59.4585,
-61.4747,
-59.7711,
-58.6101,
-58.3887,
-58.4864,
-58.0447,
-56.8286,
-55.9196,
-56.3552,
-59.9929,
-59.1182,
-60.4975,
-59.3041,
-57.4506,
-56.5174,
-57.1703,
-60.0005,
-59.6632,
-59.0451,
-59.9379,
-57.636,
-56.5565,
-58.4742,
-58.3485,
-55.38,
-59.3836,
-57.7059,
-57.7295,
-57.0076,
-56.6843,
-58.024,
-58.3442,
-58.6589,
-58.0776,
-57.2852,
-59.8936,
-58.368,
-58.8886,
-61.1485,
-57.9197,
-59.1586,
-60.8125,
-57.1896,
-56.3304,
-60.1898,
-56.4633,
};
double *out = (double *)malloc(sizeof(in));
FILE *fp = fopen("smooth.txt","wb");
printf("\n=================================\n");
if(1)
{
//const double f = 0.25;
double f = 0.01;
const size_t nsteps = 3;
const size_t n = sizeof(in)/sizeof(in[0]);
const double delta = 0.3;
double *ys = (double *) malloc(sizeof(double)*n);
double *rw = (double *) malloc(sizeof(double)*n);
double *res = (double *) malloc(sizeof(double)*n);
const double *y = &in[0];
double *x = (double *) malloc(sizeof(double)*n);
double start_freq = -125000000;
double rbw = 100000;
for (int i = 0; i < n; i++)
{
x[i] = start_freq+(rbw*i);
}
lowess((const double*)x, (const double*)y, n,f, nsteps,delta, ys, rw, res);
memcpy(out,in,sizeof(in));
int center_point = 250;
int smooth_count = 10;
for (int i = 0; i < smooth_count; i++)
{
if (abs(abs(in[(center_point-5)+i]) - abs(ys[(center_point-5)+i])) > 3)
{
out[(center_point-5)+i] = ys[(center_point-5)+i];
}
}
center_point = 2500 - 250;
for (int i = 0; i < smooth_count; i++)
{
if (abs(abs(in[(center_point-5)+i]) - abs(ys[(center_point-5)+i])) > 3)
{
out[(center_point-5)+i] = ys[(center_point-5)+i];
}
}
for (int i = 0; i < n; i++)
{
fprintf(fp,"%.4f,%.4f\r\n",in[i],out[i]);
}
fclose(fp);
fp = fopen("smooth_2.csv","wb");
const int column = 50;
double *result[column] = {NULL};
for (int i = 0; i < column; i++)
{
f = 0.001 + (i*0.005);
result[i] = (double *) malloc(sizeof(double)*n);
lowess((const double*)x, (const double*)y, n,f, nsteps,delta,result[i], rw, res);
}
for (int j = 0; j < n; j++)
{
for (int i = 0; i < column; i++)
{
fprintf(fp,"%.4f",result[i][j]);
if (i != column-1)
{
fprintf(fp,",");
}
}
fprintf(fp,"\r\n");
}
for (int i = 0; i < column; i++)
{
free(result[i]);
}
free(x);
free(ys);
free(rw);
free(res);
}
fclose(fp);
printf("create ok\n");
getchar();
}
扫一扫
6738
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
