Xisco Fauli
7003eb6fc6
Change-Id: Ia2732cc041367bc113654bc424ea6621f3961e9f Reviewed-on: https://gerrit.libreoffice.org/c/core/+/211724 Tested-by: Jenkins Reviewed-by: Xisco Fauli <xiscofauli@libreoffice.org>
87 lines
2.9 KiB
C++
87 lines
2.9 KiB
C++
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; fill-column: 100 -*- */
|
|
/*
|
|
* This file is part of the LibreOffice project.
|
|
*
|
|
* This Source Code Form is subject to the terms of the Mozilla Public
|
|
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
|
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
*/
|
|
|
|
#pragma once
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <optional>
|
|
#include <span>
|
|
|
|
namespace chart
|
|
{
|
|
/** The Pearson correlation coefficient of two runs of values, paired up by
|
|
position.
|
|
|
|
A position where either side is not a finite number counts for nothing, and
|
|
so does a position past the end of the shorter run. The answer is empty when
|
|
fewer than two positions are left, or when one of the two runs holds the
|
|
same value throughout, because then there is no coefficient to give. The
|
|
answer otherwise lies between minus one and one.
|
|
*/
|
|
inline std::optional<double> calculateCorrelationCoefficient(std::span<const double> aFirst,
|
|
std::span<const double> aSecond)
|
|
{
|
|
const std::size_t nLength = std::min(aFirst.size(), aSecond.size());
|
|
|
|
auto lclIsPair = [&aFirst, &aSecond](std::size_t nIndex) {
|
|
return std::isfinite(aFirst[nIndex]) && std::isfinite(aSecond[nIndex]);
|
|
};
|
|
|
|
std::size_t nCount = 0;
|
|
double fSumFirst = 0.0;
|
|
double fSumSecond = 0.0;
|
|
for (std::size_t nIndex = 0; nIndex < nLength; ++nIndex)
|
|
{
|
|
if (lclIsPair(nIndex))
|
|
{
|
|
++nCount;
|
|
fSumFirst += aFirst[nIndex];
|
|
fSumSecond += aSecond[nIndex];
|
|
}
|
|
}
|
|
|
|
if (nCount < 2)
|
|
return {};
|
|
|
|
const double fMeanFirst = fSumFirst / nCount;
|
|
const double fMeanSecond = fSumSecond / nCount;
|
|
|
|
// Subtracting the means before multiplying keeps the answer accurate for
|
|
// values that are large next to how far they spread.
|
|
double fProductSum = 0.0;
|
|
double fSquareSumFirst = 0.0;
|
|
double fSquareSumSecond = 0.0;
|
|
for (std::size_t nIndex = 0; nIndex < nLength; ++nIndex)
|
|
{
|
|
if (lclIsPair(nIndex))
|
|
{
|
|
const double fDeltaFirst = aFirst[nIndex] - fMeanFirst;
|
|
const double fDeltaSecond = aSecond[nIndex] - fMeanSecond;
|
|
fProductSum += fDeltaFirst * fDeltaSecond;
|
|
fSquareSumFirst += fDeltaFirst * fDeltaFirst;
|
|
fSquareSumSecond += fDeltaSecond * fDeltaSecond;
|
|
}
|
|
}
|
|
|
|
if (fSquareSumFirst <= 0.0 || fSquareSumSecond <= 0.0)
|
|
return {};
|
|
|
|
const double fCoefficient = fProductSum / std::sqrt(fSquareSumFirst * fSquareSumSecond);
|
|
|
|
// Rounding can carry the division a little past one, so hold the answer
|
|
// inside the range a coefficient has.
|
|
return std::clamp(fCoefficient, -1.0, 1.0);
|
|
}
|
|
|
|
} // namespace chart
|
|
|
|
/* vim:set shiftwidth=4 softtabstop=4 expandtab cinoptions=b1,g0,N-s cinkeys+=0=break: */
|