1fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber/*
2fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * Copyright 2012, The Android Open Source Project
3fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber *
4fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * Licensed under the Apache License, Version 2.0 (the "License");
5fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * you may not use this file except in compliance with the License.
6fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * You may obtain a copy of the License at
7fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber *
8fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber *     http://www.apache.org/licenses/LICENSE-2.0
9fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber *
10fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * Unless required by applicable law or agreed to in writing, software
11fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * distributed under the License is distributed on an "AS IS" BASIS,
12fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * See the License for the specific language governing permissions and
14fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber * limitations under the License.
15fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber */
16fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
17fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber//#define LOG_NDEBUG 0
18fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber#define LOG_TAG "LinearRegression"
19fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber#include <utils/Log.h>
20fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
21fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber#include "LinearRegression.h"
22fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
23fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber#include <math.h>
24fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber#include <string.h>
25fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
26fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Hubernamespace android {
27fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
28fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas HuberLinearRegression::LinearRegression(size_t historySize)
29fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    : mHistorySize(historySize),
30fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber      mCount(0),
31fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber      mHistory(new Point[mHistorySize]),
32fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber      mSumX(0.0),
33fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber      mSumY(0.0) {
34fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber}
35fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
36fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas HuberLinearRegression::~LinearRegression() {
37fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    delete[] mHistory;
38fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    mHistory = NULL;
39fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber}
40fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
41fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Hubervoid LinearRegression::addPoint(float x, float y) {
42fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    if (mCount == mHistorySize) {
43fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        const Point &oldest = mHistory[0];
44fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
45fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        mSumX -= oldest.mX;
46fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        mSumY -= oldest.mY;
47fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
48fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        memmove(&mHistory[0], &mHistory[1], (mHistorySize - 1) * sizeof(Point));
49fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        --mCount;
50fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    }
51fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
52fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    Point *newest = &mHistory[mCount++];
53fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    newest->mX = x;
54fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    newest->mY = y;
55fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
56fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    mSumX += x;
57fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    mSumY += y;
58fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber}
59fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
60fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huberbool LinearRegression::approxLine(float *n1, float *n2, float *b) const {
61fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    static const float kEpsilon = 1.0E-4;
62fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
63fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    if (mCount < 2) {
64fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        return false;
65fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    }
66fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
67fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float sumX2 = 0.0f;
68fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float sumY2 = 0.0f;
69fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float sumXY = 0.0f;
70fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
71fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float meanX = mSumX / (float)mCount;
72fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float meanY = mSumY / (float)mCount;
73fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
74fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    for (size_t i = 0; i < mCount; ++i) {
75fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        const Point &p = mHistory[i];
76fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
77fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        float x = p.mX - meanX;
78fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        float y = p.mY - meanY;
79fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
80fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        sumX2 += x * x;
81fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        sumY2 += y * y;
82fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        sumXY += x * y;
83fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    }
84fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
85fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float T = sumX2 + sumY2;
86fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float D = sumX2 * sumY2 - sumXY * sumXY;
87fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float root = sqrt(T * T * 0.25 - D);
88fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
89fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    float L1 = T * 0.5 - root;
90fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
91fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    if (fabs(sumXY) > kEpsilon) {
92fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        *n1 = 1.0;
93fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        *n2 = (2.0 * L1 - sumX2) / sumXY;
94fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
95fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        float mag = sqrt((*n1) * (*n1) + (*n2) * (*n2));
96fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
97fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        *n1 /= mag;
98fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        *n2 /= mag;
99fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    } else {
100fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        *n1 = 0.0;
101fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber        *n2 = 1.0;
102fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    }
103fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
104fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    *b = (*n1) * meanX + (*n2) * meanY;
105fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
106fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber    return true;
107fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber}
108fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
109fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber}  // namespace android
110fbe9d81ff5fbdc5aecdcdd13e4a5d7f019824f96Andreas Huber
111