MikkelSchubert / adapterremoval / #36

Committed 22 Jul 2024 09:33AM CUT coverage: 87.26% (-12.7%) from 100.0%

Build # #36

Build Type

push

travis-ci

Committed by

MikkelSchubert

Commit Message

remove duplicate tests

Run Details

2185 of 2504 relevant lines covered (87.26%)

16293.15 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

33.33

/src/mathutils.cpp

/*************************************************************************\
 * AdapterRemoval - cleaning next-generation sequencing reads            *
 *                                                                       *
 * Copyright (C) 2024 by Mikkel Schubert - mikkelsch@gmail.com           *
 *                                                                       *
 * This program is free software: you can redistribute it and/or modify  *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * This program is distributed in the hope that it will be useful,       *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have received a copy of the GNU General Public License     *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>. *
\*************************************************************************/
#include "mathutils.hpp" // declarations
#include "debug.hpp"     // for AR_REQUIRES
#include <array>         // for array
#include <cmath>         // for sqrt
#include <numeric>       // for accumulate

namespace adapterremoval {

namespace {

//! Critical Values of the Student's t Distribution at 0.995 for 0 to 100 df,
//! calculated via `qt(0.995, 1:100)` in R
const std::array<double, 101> STUDENTS_T_CRITICAL_VALUES = {
  NAN,    63.657, 9.9248, 5.8409, 4.6041, 4.0321, 3.7074, 3.4995, 3.3554,
  3.2498, 3.1693, 3.1058, 3.0545, 3.0123, 2.9768, 2.9467, 2.9208, 2.8982,
  2.8784, 2.8609, 2.8453, 2.8314, 2.8188, 2.8073, 2.7969, 2.7874, 2.7787,
  2.7707, 2.7633, 2.7564, 2.7500, 2.7440, 2.7385, 2.7333, 2.7284, 2.7238,
  2.7195, 2.7154, 2.7116, 2.7079, 2.7045, 2.7012, 2.6981, 2.6951, 2.6923,
  2.6896, 2.6870, 2.6846, 2.6822, 2.6800, 2.6778, 2.6757, 2.6737, 2.6718,
  2.6700, 2.6682, 2.6665, 2.6649, 2.6633, 2.6618, 2.6603, 2.6589, 2.6575,
  2.6561, 2.6549, 2.6536, 2.6524, 2.6512, 2.6501, 2.6490, 2.6479, 2.6469,
  2.6459, 2.6449, 2.6439, 2.6430, 2.6421, 2.6412, 2.6403, 2.6395, 2.6387,
  2.6379, 2.6371, 2.6364, 2.6356, 2.6349, 2.6342, 2.6335, 2.6329, 2.6322,
  2.6316, 2.6309, 2.6303, 2.6297, 2.6291, 2.6286, 2.6280, 2.6275, 2.6269,
  2.6264, 2.6259
};

} // namespace

double
arithmetic_mean(const std::vector<uint64_t>& values)
{
  AR_REQUIRE(!values.empty());

  return std::accumulate(values.begin(), values.end(), uint64_t()) /
         static_cast<double>(values.size());
}

double
standard_deviation(const std::vector<uint64_t>& values)
{
  AR_REQUIRE(values.size() > 1);

  double error = 0;
  const auto m = arithmetic_mean(values);
  for (const auto it : values) {
    error += (it - m) * (it - m);
  }

  return std::sqrt(error / (values.size() - 1));
}

double
students_t_critical_value(size_t df)
{
  return STUDENTS_T_CRITICAL_VALUES.at(
    std::min<size_t>(df, STUDENTS_T_CRITICAL_VALUES.size() - 1));
}

bool
grubbs_test_prune(std::vector<uint64_t>& values)
{
  const size_t original_size = values.size();
  bool found_any = false;

  do {
    const double N = values.size();
    const double t = students_t_critical_value(N - 2);
    const auto t2 = t * t;

    const double m = arithmetic_mean(values);
    const double sd = standard_deviation(values);

    const double cutoff = ((N - 1) / std::sqrt(N)) * sqrt(t2 / (N - 2 + t2));

    found_any = false;
    for (size_t i = 0; i < values.size(); ++i) {
      const auto G = (values.at(i) - m) / sd;
      if (G < -cutoff || G > cutoff) {
        values.at(i) = values.back();
        values.pop_back();
        found_any = true;
        break;
      }
    }
  } while (found_any);

  return original_size != values.size();
}

} // namespace adapterremoval

1	/*************************************************************************\
2	* AdapterRemoval - cleaning next-generation sequencing reads *
3	* *
4	* Copyright (C) 2024 by Mikkel Schubert - mikkelsch@gmail.com *
5	* *
6	* This program is free software: you can redistribute it and/or modify *
7	* it under the terms of the GNU General Public License as published by *
8	* the Free Software Foundation, either version 3 of the License, or *
9	* (at your option) any later version. *
10	* *
11	* This program is distributed in the hope that it will be useful, *
12	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
13	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14	* GNU General Public License for more details. *
15	* *
16	* You should have received a copy of the GNU General Public License *
17	* along with this program. If not, see <http://www.gnu.org/licenses/>. *
18	\*************************************************************************/
19	#include "mathutils.hpp" // declarations
20	#include "debug.hpp" // for AR_REQUIRES
21	#include <array> // for array
22	#include <cmath> // for sqrt
23	#include <numeric> // for accumulate
24
25	namespace adapterremoval {
26
27	namespace {
28
29	//! Critical Values of the Student's t Distribution at 0.995 for 0 to 100 df,
30	//! calculated via `qt(0.995, 1:100)` in R
31	const std::array<double, 101> STUDENTS_T_CRITICAL_VALUES = {
32	NAN, 63.657, 9.9248, 5.8409, 4.6041, 4.0321, 3.7074, 3.4995, 3.3554,
33	3.2498, 3.1693, 3.1058, 3.0545, 3.0123, 2.9768, 2.9467, 2.9208, 2.8982,
34	2.8784, 2.8609, 2.8453, 2.8314, 2.8188, 2.8073, 2.7969, 2.7874, 2.7787,
35	2.7707, 2.7633, 2.7564, 2.7500, 2.7440, 2.7385, 2.7333, 2.7284, 2.7238,
36	2.7195, 2.7154, 2.7116, 2.7079, 2.7045, 2.7012, 2.6981, 2.6951, 2.6923,
37	2.6896, 2.6870, 2.6846, 2.6822, 2.6800, 2.6778, 2.6757, 2.6737, 2.6718,
38	2.6700, 2.6682, 2.6665, 2.6649, 2.6633, 2.6618, 2.6603, 2.6589, 2.6575,
39	2.6561, 2.6549, 2.6536, 2.6524, 2.6512, 2.6501, 2.6490, 2.6479, 2.6469,
40	2.6459, 2.6449, 2.6439, 2.6430, 2.6421, 2.6412, 2.6403, 2.6395, 2.6387,
41	2.6379, 2.6371, 2.6364, 2.6356, 2.6349, 2.6342, 2.6335, 2.6329, 2.6322,
42	2.6316, 2.6309, 2.6303, 2.6297, 2.6291, 2.6286, 2.6280, 2.6275, 2.6269,
43	2.6264, 2.6259
44	};
45
46	} // namespace
47
48	double
49	arithmetic_mean(const std::vector<uint64_t>& values)	8✔
50	{
51	AR_REQUIRE(!values.empty());	20✔
52
53	return std::accumulate(values.begin(), values.end(), uint64_t()) /	21✔
54	static_cast<double>(values.size());	14✔
55	}
56
57	double
58	standard_deviation(const std::vector<uint64_t>& values)	6✔
59	{
60	AR_REQUIRE(values.size() > 1);	20✔
61
62	double error = 0;	4✔
63	const auto m = arithmetic_mean(values);	4✔
64	for (const auto it : values) {	49✔
65	error += (it - m) * (it - m);	11✔
66	}
67
68	return std::sqrt(error / (values.size() - 1));	8✔
69	}
70
71	double
72	students_t_critical_value(size_t df)	×
73	{
74	return STUDENTS_T_CRITICAL_VALUES.at(	×
75	std::min<size_t>(df, STUDENTS_T_CRITICAL_VALUES.size() - 1));	×
76	}
77
78	bool
79	grubbs_test_prune(std::vector<uint64_t>& values)	×
80	{
81	const size_t original_size = values.size();	×
82	bool found_any = false;	×
83
84	do {	×
85	const double N = values.size();	×
86	const double t = students_t_critical_value(N - 2);	×
87	const auto t2 = t * t;	×
88
89	const double m = arithmetic_mean(values);	×
90	const double sd = standard_deviation(values);	×
91
92	const double cutoff = ((N - 1) / std::sqrt(N)) * sqrt(t2 / (N - 2 + t2));	×
93
94	found_any = false;	×
95	for (size_t i = 0; i < values.size(); ++i) {	×
96	const auto G = (values.at(i) - m) / sd;	×
97	if (G < -cutoff \|\| G > cutoff) {	×
98	values.at(i) = values.back();	×
99	values.pop_back();	×
100	found_any = true;	×
101	break;	×
102	}
103	}
104	} while (found_any);
105
106	return original_size != values.size();	×
107	}
108
109	} // namespace adapterremoval

MikkelSchubert / adapterremoval / #36

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous