/* 
 *------------------------------------------------------------------
 * Copyright (c) 2006-2016 Cisco and/or its affiliates.
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at:
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

/* see "Numerical Recipies in C, 2nd ed." p 665 */

#include <stdio.h>
#include <math.h>

/*
 * linreg
 * Linear regression of (xi, yi), returns parameters for least-squares
 * fit y = a + bx.  Also, compute Pearson's R.
 */
void linreg (double *x, double *y, int nitems, double *a, double *b,
             double *minp, double *maxp, double *meanp, double *r)
{
    double sx = 0.0;
    double sy = 0.0;
    double st2 = 0.0;
    double min = y[0];
    double max = 0.0;
    double ss, meanx, meany, t;
    double errx, erry, prodsum, sqerrx, sqerry;
    int i;

    *b = 0.0;
    
    for (i = 0; i < nitems; i++) {
        sx += x[i];
        sy += y[i];
        if (y[i] < min)
            min = y[i];
        if (y[i] > max)
            max = y[i];
    }
    ss = nitems;
    meanx = sx / ss;
    meany = *meanp = sy / ss;
    *minp = min;
    *maxp = max;

    for (i = 0; i < nitems; i++) {
        t = x[i] - meanx;
        st2 += t*t;
        *b += t*y[i];
    }

    *b /= st2;
    *a = (sy-sx*(*b))/ss;

    prodsum = 0.0;
    sqerrx = 0.0;
    sqerry = 0.0;

    /* Compute numerator of Pearson's R */
    for (i = 0; i < nitems; i++) {
        errx = x[i] - meanx;
        erry = y[i] - meany;
        prodsum += errx * erry;
        sqerrx += errx*errx;
        sqerry += erry*erry;
    }

    *r = prodsum / (sqrt(sqerrx)*sqrt(sqerry));
}