Skip to content

Commit e75fd3d

Browse files
jangoreckijangorecki
committed
fix broken automatic git merge
1 parent c410d6d commit e75fd3d

File tree

2 files changed

+0
-659
lines changed

2 files changed

+0
-659
lines changed

‎src/froll.c‎

Lines changed: 0 additions & 387 deletions
Original file line numberDiff line numberDiff line change
@@ -1071,390 +1071,3 @@ void frollprodExact(double *x, uint64_t nx, ans_t *ans, int k, double fill, bool
10711071
}
10721072
}
10731073
}
1074-
1075-
static inline void wmin(const double * restrict x, uint64_t o, int k, double * restrict w, uint64_t *iw, bool narm) {
1076-
if (narm) {
1077-
for (int i=0; i<k; i++) {
1078-
if (x[o+i-k+1] <= w[0]) { // this never true if all x NAs and narm=TRUE
1079-
iw[0] = o+i-k+1;
1080-
w[0] = x[iw[0]];
1081-
}
1082-
}
1083-
} else {
1084-
double ww = R_PosInf;
1085-
uint64_t iww = 0;
1086-
for (int i=0; i<k; i++) {
1087-
uint64_t ii = o+i-k+1;
1088-
if (ISNAN(x[ii])) {
1089-
if (ISNA(x[ii])) {
1090-
error("internal error: frollmin reached untested branch of code, for now use frollmin algo='exact', please provide reproducible example of your frollmin usage to data.table github issue tracker\n"); // # nocov
1091-
//iww = ii; ww = NA_REAL;
1092-
} else if (ISNA(ww)) {
1093-
error("internal error: frollmin reached untested branch of code, for now use frollmin algo='exact', please provide reproducible example of your frollmin usage to data.table github issue tracker\n"); // # nocov
1094-
// do nothing because w > x[i]: NA > NaN
1095-
} else { // no NA in window so NaN >= than any non-NA
1096-
iww = ii; ww = R_NaN;
1097-
}
1098-
} else if (ISNAN(ww)) {
1099-
// w still within the window and is NA or NaN, x[i] is not NA - already checked above, therefore to nothing
1100-
} else if (x[ii] <= ww) {
1101-
iww = ii; ww = x[iww];
1102-
}
1103-
}
1104-
iw[0] = iww;
1105-
w[0] = ww;
1106-
}
1107-
}
1108-
/* fast rolling min - fast
1109-
* see rolling max fast details
1110-
*/
1111-
void frollminFast(double *x, uint64_t nx, ans_t *ans, int k, double fill, bool narm, int hasnf, bool verbose) {
1112-
if (verbose)
1113-
snprintf(end(ans->message[0]), 500, _("%s: running for input length %"PRIu64", window %d, hasnf %d, narm %d\n"), "frollminFast", (uint64_t)nx, k, hasnf, (int)narm);
1114-
if (k == 0) {
1115-
if (verbose)
1116-
snprintf(end(ans->message[0]), 500, _("%s: window width of size 0, returning all +Inf vector\n"), __func__);
1117-
for (uint64_t i=0; i<nx; i++) {
1118-
ans->dbl_v[i] = R_PosInf;
1119-
}
1120-
return;
1121-
}
1122-
double w = R_PosInf; // window min
1123-
uint64_t cmin = 0; // counter of nested loops for verbose
1124-
uint64_t iw = 0; // index of window min
1125-
uint64_t i;
1126-
if (narm || hasnf==-1) {
1127-
for (i=0; i<k-1; i++) { // #loop_counter_not_local_scope_ok
1128-
if (x[i] <= w) { // <= rather than < because we track most recent minimum using iw
1129-
iw = i; w = x[iw];
1130-
}
1131-
ans->dbl_v[i] = fill;
1132-
}
1133-
for (i=k-1; i<nx; i++) {
1134-
if (iw+k <= i) { // min left current window // note that it is still <= same as in max
1135-
iw = i-k; w = R_PosInf;
1136-
wmin(x, i, k, &w, &iw, true); cmin++;
1137-
} else if (x[i] <= w) {
1138-
iw = i; w = x[iw];
1139-
}
1140-
ans->dbl_v[i] = w;
1141-
}
1142-
} else {
1143-
bool truehasnf = hasnf>0;
1144-
for (i=0; i<k-1; i++) { // up to first full window only #loop_counter_not_local_scope_ok
1145-
if (ISNAN(x[i])) {
1146-
truehasnf = true;
1147-
if (ISNA(x[i])) {
1148-
iw = i; w = NA_REAL;
1149-
} else if (ISNA(w)) {
1150-
// do nothing because w > x[i]: NA > NaN
1151-
} else {
1152-
iw = i; w = R_NaN;
1153-
}
1154-
} else if (x[i] <= w) {
1155-
iw = i; w = x[iw];
1156-
}
1157-
ans->dbl_v[i] = fill;
1158-
}
1159-
if (!truehasnf) { // maybe no NAs
1160-
for (; i<nx; i++) {
1161-
if (ISNAN(x[i])) { // Inf properly propagates
1162-
if (verbose)
1163-
ansSetMsg(ans, 0, "%s: non-finite values are present in input, continue with extra care for NFs\n", __func__);
1164-
truehasnf = true;
1165-
break; // does not increment, w stays as from previous iteration
1166-
}
1167-
if (iw+k <= i) { // min left current window // note that it is still <= same as in max
1168-
iw = i-k; w = R_PosInf;
1169-
wmin(x, i, k, &w, &iw, true); cmin++;
1170-
} else if (x[i] <= w) {
1171-
iw = i; w = x[iw];
1172-
}
1173-
ans->dbl_v[i] = w;
1174-
}
1175-
}
1176-
if (truehasnf) {
1177-
for (; i<nx; i++) { // this loop continues from where "maybe no NAs" loop left
1178-
if (ISNAN(x[i])) {
1179-
if (ISNA(x[i])) {
1180-
iw = i; w = NA_REAL;
1181-
} else if (ISNA(w)) {
1182-
// do nothing because w > x[i]: NA > NaN
1183-
} else { // no NA in window so NaN >= than any non-NA
1184-
iw = i; w = R_NaN;
1185-
}
1186-
} else if (iw+k <= i) { // min left current window // note that it is still <= same as in max
1187-
iw = i-k; w = R_PosInf;
1188-
wmin(x, i, k, &w, &iw, false); cmin++;
1189-
} else if (ISNAN(w)) {
1190-
// w still within the window and is NA or NaN, x[i] is not NA - already checked above, therefore do nothing
1191-
} else if (x[i] <= w) {
1192-
iw = i; w = x[iw];
1193-
}
1194-
ans->dbl_v[i] = w;
1195-
}
1196-
}
1197-
}
1198-
if (verbose)
1199-
snprintf(end(ans->message[0]), 500, _("%s: nested window min calculation called %"PRIu64" times\n"), __func__, cmin);
1200-
}
1201-
/* fast rolling min - exact
1202-
* see rolling max exact details
1203-
*/
1204-
void frollminExact(double *x, uint64_t nx, ans_t *ans, int k, double fill, bool narm, int hasnf, bool verbose) {
1205-
if (verbose)
1206-
snprintf(end(ans->message[0]), 500, _("%s: running in parallel for input length %"PRIu64", window %d, hasnf %d, narm %d\n"), "frollminExact", (uint64_t)nx, k, hasnf, (int)narm);
1207-
if (k == 0) {
1208-
if (verbose)
1209-
snprintf(end(ans->message[0]), 500, _("%s: window width of size 0, returning all +Inf vector\n"), __func__);
1210-
for (uint64_t i=0; i<nx; i++) {
1211-
ans->dbl_v[i] = R_PosInf;
1212-
}
1213-
return;
1214-
}
1215-
for (int i=0; i<k-1; i++) {
1216-
ans->dbl_v[i] = fill;
1217-
}
1218-
if (narm || hasnf==-1) { // ignore NAs as > does not propagate
1219-
#pragma omp parallel for num_threads(getDTthreads(nx, true))
1220-
for (uint64_t i=k-1; i<nx; i++) {
1221-
double w = R_PosInf;
1222-
for (int j=-k+1; j<=0; j++) {
1223-
if (x[i+j] < w)
1224-
w = x[i+j];
1225-
}
1226-
ans->dbl_v[i] = w;
1227-
}
1228-
} else {
1229-
bool *isnan = malloc(nx*sizeof(bool)); // isnan lookup - we use it to reduce ISNAN calls in nested loop
1230-
if (!isnan) { // # nocov start
1231-
ansSetMsg(ans, 3, "%s: Unable to allocate memory for isnan", __func__); // raise error
1232-
free(isnan);
1233-
return;
1234-
} // # nocov end
1235-
bool truehasnf = hasnf>0;
1236-
for (uint64_t i=0; i<nx; i++) { // no openmp as this should be very fast
1237-
if (ISNAN(x[i])) {
1238-
truehasnf = true;
1239-
isnan[i] = true;
1240-
} else {
1241-
isnan[i] = false;
1242-
}
1243-
}
1244-
if (!truehasnf) { // not found any NAs
1245-
#pragma omp parallel for num_threads(getDTthreads(nx, true))
1246-
for (uint64_t i=k-1; i<nx; i++) {
1247-
double w = R_PosInf;
1248-
for (int j=-k+1; j<=0; j++) {
1249-
if (x[i+j] < w)
1250-
w = x[i+j];
1251-
}
1252-
ans->dbl_v[i] = w;
1253-
}
1254-
} else { // there are some NAs
1255-
#pragma omp parallel for num_threads(getDTthreads(nx, true))
1256-
for (uint64_t i=k-1; i<nx; i++) {
1257-
double w = R_PosInf;
1258-
if (isnan[i] && ISNA(x[i])) {
1259-
w = NA_REAL;
1260-
} else {
1261-
for (int j=-k+1; j<=0; j++) {
1262-
if (isnan[i+j]) {
1263-
if (ISNA(x[i+j])) {
1264-
w = NA_REAL;
1265-
break; // break because NA > NaN
1266-
} else {
1267-
w = R_NaN; // continue nested loop in case there is NA there
1268-
}
1269-
} else if (x[i+j] < w)
1270-
w = x[i+j];
1271-
}
1272-
}
1273-
ans->dbl_v[i] = w;
1274-
}
1275-
}
1276-
}
1277-
}
1278-
1279-
/* fast rolling prod - fast
1280-
* same as mean fast
1281-
*/
1282-
void frollprodFast(double *x, uint64_t nx, ans_t *ans, int k, double fill, bool narm, int hasnf, bool verbose) {
1283-
if (verbose)
1284-
snprintf(end(ans->message[0]), 500, _("%s: running for input length %"PRIu64", window %d, hasnf %d, narm %d\n"), "frollprodFast", (uint64_t)nx, k, hasnf, (int)narm);
1285-
if (k == 0) {
1286-
if (verbose)
1287-
snprintf(end(ans->message[0]), 500, _("%s: window width of size 0, returning all 1 vector\n"), __func__);
1288-
for (uint64_t i=0; i<nx; i++) {
1289-
ans->dbl_v[i] = 1.0;
1290-
}
1291-
return;
1292-
}
1293-
long double w = 1.0;
1294-
bool truehasnf = hasnf>0;
1295-
if (!truehasnf) {
1296-
int i;
1297-
for (i=0; i<k-1; i++) { // #loop_counter_not_local_scope_ok
1298-
w *= x[i];
1299-
ans->dbl_v[i] = fill;
1300-
}
1301-
w *= x[i];
1302-
ans->dbl_v[i] = (double) w;
1303-
if (R_FINITE((double) w)) {
1304-
for (uint64_t i=k; i<nx; i++) {
1305-
w /= x[i-k];
1306-
w *= x[i];
1307-
ans->dbl_v[i] = (double) w;
1308-
}
1309-
if (!R_FINITE((double) w)) {
1310-
if (hasnf==-1)
1311-
ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__);
1312-
if (verbose)
1313-
ansSetMsg(ans, 0, "%s: non-finite values are present in input, re-running with extra care for NFs\n", __func__);
1314-
w = 1.0; truehasnf = true;
1315-
}
1316-
} else {
1317-
if (hasnf==-1)
1318-
ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__);
1319-
if (verbose)
1320-
ansSetMsg(ans, 0, "%s: non-finite values are present in input, skip non-finite inaware attempt and run with extra care for NFs straighaway\n", __func__);
1321-
w = 1.0; truehasnf = true;
1322-
}
1323-
}
1324-
if (truehasnf) {
1325-
int nc = 0, pinf = 0, ninf = 0; // NA counter within sliding window
1326-
int i; // iterator declared here because it is being used after for loop
1327-
1328-
#undef PROD_WINDOW_STEP_FRONT
1329-
#define PROD_WINDOW_STEP_FRONT \
1330-
if (R_FINITE(x[i])) { \
1331-
w *= x[i]; \
1332-
} else if (ISNAN(x[i])) { \
1333-
nc++; \
1334-
} else if (x[i]==R_PosInf) { \
1335-
pinf++; \
1336-
} else if (x[i]==R_NegInf) { \
1337-
ninf++; \
1338-
}
1339-
#undef PROD_WINDOW_STEP_BACK
1340-
#define PROD_WINDOW_STEP_BACK \
1341-
if (R_FINITE(x[i-k])) { \
1342-
w /= x[i-k]; \
1343-
} else if (ISNAN(x[i-k])) { \
1344-
nc--; \
1345-
} else if (x[i-k]==R_PosInf) { \
1346-
pinf--; \
1347-
} else if (x[i-k]==R_NegInf) { \
1348-
ninf--; \
1349-
}
1350-
#undef PROD_WINDOW_STEP_VALUE
1351-
#define PROD_WINDOW_STEP_VALUE \
1352-
if (nc == 0) { \
1353-
if (pinf == 0 && ninf == 0) { \
1354-
ans->dbl_v[i] = (double) w; \
1355-
} else { \
1356-
ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \
1357-
} \
1358-
} else if (nc == k) { \
1359-
ans->dbl_v[i] = narm ? 1.0 : NA_REAL; \
1360-
} else { \
1361-
if (narm) { \
1362-
if (pinf == 0 && ninf == 0) { \
1363-
ans->dbl_v[i] = (double) w; \
1364-
} else { \
1365-
ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \
1366-
} \
1367-
} else { \
1368-
ans->dbl_v[i] = NA_REAL; \
1369-
} \
1370-
}
1371-
1372-
for (i=0; i<k-1; i++) { // loop over leading observation, all partial window only; #loop_counter_not_local_scope_ok
1373-
PROD_WINDOW_STEP_FRONT
1374-
ans->dbl_v[i] = fill; // partial window fill all
1375-
}
1376-
PROD_WINDOW_STEP_FRONT // i==k-1
1377-
PROD_WINDOW_STEP_VALUE
1378-
for (uint64_t i=k; i<nx; i++) { // loop over obs, complete window, all remaining after partial window
1379-
PROD_WINDOW_STEP_BACK
1380-
PROD_WINDOW_STEP_FRONT
1381-
PROD_WINDOW_STEP_VALUE
1382-
}
1383-
}
1384-
}
1385-
/* fast rolling prod - exact
1386-
* same as mean exact
1387-
*/
1388-
void frollprodExact(double *x, uint64_t nx, ans_t *ans, int k, double fill, bool narm, int hasnf, bool verbose) {
1389-
if (verbose)
1390-
snprintf(end(ans->message[0]), 500, _("%s: running in parallel for input length %"PRIu64", window %d, hasnf %d, narm %d\n"), "frollprodExact", (uint64_t)nx, k, hasnf, (int)narm);
1391-
if (k == 0) {
1392-
if (verbose)
1393-
snprintf(end(ans->message[0]), 500, _("%s: window width of size 0, returning all 1 vector\n"), __func__);
1394-
for (uint64_t i=0; i<nx; i++) {
1395-
ans->dbl_v[i] = 1.0;
1396-
}
1397-
return;
1398-
}
1399-
for (int i=0; i<k-1; i++) {
1400-
ans->dbl_v[i] = fill;
1401-
}
1402-
bool truehasnf = hasnf>0;
1403-
if (!truehasnf || !narm) {
1404-
#pragma omp parallel for num_threads(getDTthreads(nx, true))
1405-
for (uint64_t i=k-1; i<nx; i++) {
1406-
if (narm && truehasnf) {
1407-
continue;
1408-
}
1409-
long double w = 1.0;
1410-
for (int j=-k+1; j<=0; j++) {
1411-
w *= x[i+j];
1412-
}
1413-
if (R_FINITE((double) w)) {
1414-
ans->dbl_v[i] = (double) w;
1415-
} else if (ISNAN((double) w)) {
1416-
if (!narm) {
1417-
ans->dbl_v[i] = (double) w;
1418-
}
1419-
truehasnf = true;
1420-
} else {
1421-
ans->dbl_v[i] = (double) w;
1422-
}
1423-
}
1424-
if (truehasnf) {
1425-
if (hasnf==-1)
1426-
ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__);
1427-
if (verbose) {
1428-
if (narm)
1429-
ansSetMsg(ans, 0, "%s: non-finite values are present in input, re-running with extra care for NFs\n", __func__);
1430-
else
1431-
ansSetMsg(ans, 0, "%s: non-finite values are present in input, na.rm=FALSE and algo='exact' propagates NFs properply, no need to re-run\n", __func__);
1432-
}
1433-
}
1434-
}
1435-
if (truehasnf && narm) {
1436-
#pragma omp parallel for num_threads(getDTthreads(nx, true))
1437-
for (uint64_t i=k-1; i<nx; i++) {
1438-
long double w = 1.0;
1439-
int nc = 0;
1440-
for (int j=-k+1; j<=0; j++) {
1441-
if (ISNAN(x[i+j])) {
1442-
nc++;
1443-
} else {
1444-
w *= x[i+j];
1445-
}
1446-
}
1447-
if (w > DBL_MAX) { // in contrast to mean, here we can overflow long double more than DBL_MAX
1448-
ans->dbl_v[i] = R_PosInf;
1449-
} else if (w < -DBL_MAX) {
1450-
ans->dbl_v[i] = R_NegInf;
1451-
} else {
1452-
if (nc < k) {
1453-
ans->dbl_v[i] = (double) w;
1454-
} else {
1455-
ans->dbl_v[i] = 1.0;
1456-
}
1457-
}
1458-
}
1459-
}
1460-
}

0 commit comments

Comments
 (0)