diff --git a/DESCRIPTION b/DESCRIPTION index beafe8a..aa9f4ba 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,23 +1,29 @@ Package: lfe -Version: 2.9-1 -Date: 2023-04-03 +Version: 3.0-0 Title: Linear Group Fixed Effects -Authors@R: c(person("Simen", "Gaure", email="Simen.Gaure@frisch.uio.no", role=c("aut"), - comment=c(ORCID="https://orcid.org/0000-0001-7251-8747")), - person("Grant","McDermott", email="grantmcd@uoregon.edu", role="ctb"), - person(given = "Mauricio", - family = "Vargas", - role = "ctb", - email = "mavargas11@uc.cl", - comment = c(ORCID = "0000-0003-1017-7574")), - person(given = "Paulo", - family = "Alcazar", - role = "ctb"), - person("Karl", "Dunkle Werner", role="ctb"), - person("Matthieu","Stigler", email = "Matthieu.Stigler@gmail.com", role= c("ctb", "cre"), - comment =c(ORCID="0000-0002-6802-4290")), - person("Daniel","Lüdecke",email="mail@danielluedecke.de",role="ctb")) -Copyright: 2011-2019, Simen Gaure +Authors@R: c( + person(given = "Simen", + family = "Gaure", + role = "aut", + comment = c(ORCID="https://orcid.org/0000-0001-7251-8747")), + person(given = "Grant", + family = "McDermott", + role = "ctb"), + person(given = "Mauricio", + family = "Vargas Sepulveda", + role = c("ctb", "cre"), + email = "m.sepulveda@mail.utoronto.ca", + comment = c(ORCID = "0000-0003-1017-7574")), + person(given = "Karl", + family = "Dunkle Werner", + role ="ctb"), + person(given = "Matthieu", + family = "Stigler", + role = "ctb", + comment = c(ORCID="0000-0002-6802-4290")), + person(given = "Daniel", + family = "Lüdecke", + role ="ctb")) Depends: R (>= 2.15.2), Matrix (>= 1.1-2) Imports: Formula, xtable, compiler, utils, methods, sandwich, parallel Suggests: @@ -36,13 +42,13 @@ Description: Transforms away factors with many levels prior to doing an OLS. Useful for estimating linear models with multiple group fixed effects, and for estimating linear models which uses factors with many levels as pure control variables. See Gaure (2013) Includes support for instrumental variables, conditional F statistics for weak instruments, - robust and multi-way clustered standard errors, as well as limited mobility bias correction (Gaure 2014 ). - WARNING: This package is NOT under active development anymore, no further improvements are to be expected, and the package is at risk of being removed from CRAN. -License: Artistic-2.0 + robust and multi-way clustered standard errors, as well as limited mobility bias correction (Gaure 2014 ). + Since version 3.0, it provides dedicated functions to estimate Poisson models. +License: MIT + file LICENSE Classification/JEL: C13, C23, C60 Classification/MSC: 62J05, 65F10, 65F50 -URL: https://github.com/MatthieuStigler/lfe -BugReports: https://github.com/MatthieuStigler/lfe/issues +URL: https://github.com/r-econometrics/lfe +BugReports: https://github.com/r-econometrics/lfe/issues Encoding: UTF-8 RoxygenNote: 7.2.3 Roxygen: list(markdown = TRUE) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..a1e92e2 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +YEAR: 2023 +COPYRIGHT HOLDER: lfe authors diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..289f91e --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,21 @@ +# MIT License + +Copyright (c) 2023 lfe authors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/R/feglm.R b/R/feglm.R index fb3d24f..dbe36b3 100644 --- a/R/feglm.R +++ b/R/feglm.R @@ -10,6 +10,8 @@ #' @param cluster optional variable to group by and compute sandwich-type #' robust standard errors. Should be a formula of the form `~x_j` or #' an object that be coerced to a formula. +#' @param pseudo_rsq logical value to return a a pseudo-R2 based on Kendall's +#' correlation. #' @param tol tolerance value for GLM convergence criteria. #' @importFrom Formula Formula #' @importFrom Matrix Diagonal @@ -21,6 +23,7 @@ fepois <- function(formula, data, subset = NULL, robust = TRUE, cluster = NULL, + pseudo_rsq = FALSE, tol = 1e-10) { if (!is.null(subset)) { data <- data[subset, ] } @@ -41,6 +44,10 @@ fepois <- function(formula, data, vardep <- all.vars(formula(formula, lhs = 1, rhs = 0)) vardep <- data[, vardep, drop = TRUE] + if (isTRUE(pseudo_rsq)) { + vardep2 <- vardep + } + if (min(vardep) < 0) { stop("y should be greater or equals to zero.") } @@ -181,6 +188,10 @@ fepois <- function(formula, data, names(reg$fitted.values) <- rownames(data) + if (isTRUE(pseudo_rsq)) { + reg$pseudo_rsq <- cor(vardep2, reg$fitted.values, method = "kendall")^2 + } + class(reg) <- "fepois" return(reg) } diff --git a/dev/test-ppml.R b/dev/test-ppml.R new file mode 100644 index 0000000..ef51a03 --- /dev/null +++ b/dev/test-ppml.R @@ -0,0 +1,73 @@ +load_all() + +library(dplyr) + +ch1_application1 <- tradepolicy::agtpa_applications %>% + select(exporter, importer, pair_id, year, trade, dist, cntg, lang, clny) %>% + filter(year %in% seq(1986, 2006, 4)) + +ch1_application1 <- ch1_application1 %>% + mutate( + log_trade = log(trade), + log_dist = log(dist) + ) + +ch1_application1 <- ch1_application1 %>% + # Create Yit + group_by(exporter, year) %>% + mutate( + y = sum(trade), + log_y = log(y) + ) %>% + + # Create Eit + group_by(importer, year) %>% + mutate( + e = sum(trade), + log_e = log(e) + ) + +ch1_application1 <- ch1_application1 %>% + # Replicate total_e + group_by(exporter, year) %>% + mutate(total_e = sum(e)) %>% + group_by(year) %>% + mutate(total_e = max(total_e)) %>% + + # Replicate rem_exp + group_by(exporter, year) %>% + mutate( + remoteness_exp = sum(dist * total_e / e), + log_remoteness_exp = log(remoteness_exp) + ) %>% + + # Replicate total_y + group_by(importer, year) %>% + mutate(total_y = sum(y)) %>% + group_by(year) %>% + mutate(total_y = max(total_y)) %>% + + # Replicate rem_imp + group_by(importer, year) %>% + mutate( + remoteness_imp = sum(dist / (y / total_y)), + log_remoteness_imp = log(remoteness_imp) + ) + +ch1_application1 <- ch1_application1 %>% + # This merges the columns exporter/importer with year + mutate( + exp_year = paste0(exporter, year), + imp_year = paste0(importer, year) + ) + +ch1_application1 <- ch1_application1 %>% + filter(exporter != importer) + +fit_ppml <- lfe::fepois(trade ~ log_dist + cntg + lang + clny | exp_year + imp_year, + data = ch1_application1, + cluster = ~pair_id, + pseudo_rsq = T +) + +summary(fit_ppml) diff --git a/man/fepois.Rd b/man/fepois.Rd index 93abe1a..a2fb8d7 100644 --- a/man/fepois.Rd +++ b/man/fepois.Rd @@ -11,6 +11,7 @@ fepois( subset = NULL, robust = TRUE, cluster = NULL, + pseudo_rsq = FALSE, tol = 1e-10 ) } @@ -37,6 +38,9 @@ used in the fitting process.} robust standard errors. Should be a formula of the form \code{~x_j} or an object that be coerced to a formula.} +\item{pseudo_rsq}{logical value to return a a pseudo-R2 based on Kendall's +correlation.} + \item{tol}{tolerance value for GLM convergence criteria.} } \description{ diff --git a/src/crowsum.c b/src/crowsum.c index d630818..3826a26 100644 --- a/src/crowsum.c +++ b/src/crowsum.c @@ -2,10 +2,8 @@ // #ifdef _OPENMP // #include // #endif -SEXP crowsum(SEXP Rmat, SEXP Rfactor, SEXP Rmean) -{ - if (!IS_NUMERIC(Rmat)) - error("Only numeric matrices accepted"); +SEXP crowsum(SEXP Rmat, SEXP Rfactor, SEXP Rmean) { + if (!IS_NUMERIC(Rmat)) error("Only numeric matrices accepted"); if (!isInteger(Rfactor) && !isFactor(Rfactor)) error("Only factor or integer vector accepted"); R_xlen_t len = xlength(Rmat); @@ -19,13 +17,10 @@ SEXP crowsum(SEXP Rmat, SEXP Rfactor, SEXP Rmean) int *table = NULL; // int nthr = INTEGER(threads)[0]; mat = REAL(Rmat); - if (isMatrix(Rmat)) - { + if (isMatrix(Rmat)) { cols = ncols(Rmat); rows = nrows(Rmat); - } - else - { + } else { cols = 1; rows = len; } @@ -34,21 +29,15 @@ SEXP crowsum(SEXP Rmat, SEXP Rfactor, SEXP Rmean) nlev = nlevels(Rfactor); - for (int i = 0; i < rows; i++) - { - if (f[i] < 1 || ISNA(f[i])) - error("Missing levels not supported"); - if (f[i] > nlev) - error("Level for %d is %d, too large %d", i, f[i], nlev); + for (int i = 0; i < rows; i++) { + if (f[i] < 1 || ISNA(f[i])) error("Missing levels not supported"); + if (f[i] > nlev) error("Level for %d is %d, too large %d", i, f[i], nlev); } - if (mean) - { + if (mean) { table = (int *)R_alloc(nlev, sizeof(int)); - for (int i = 0; i < nlev; i++) - table[i] = 0; - for (int i = 0; i < rows; i++) - table[f[i] - 1]++; + for (int i = 0; i < nlev; i++) table[i] = 0; + for (int i = 0; i < rows; i++) table[f[i] - 1]++; } // Allocate resultant matrix @@ -58,8 +47,7 @@ SEXP crowsum(SEXP Rmat, SEXP Rfactor, SEXP Rmean) SEXP rdn = GET_DIMNAMES(Rmat); PROTECT(dn = allocVector(VECSXP, 2)); SET_VECTOR_ELT(dn, 0, GET_LEVELS(Rfactor)); - if (!isNull(rdn)) - SET_VECTOR_ELT(dn, 1, VECTOR_ELT(rdn, 1)); + if (!isNull(rdn)) SET_VECTOR_ELT(dn, 1, VECTOR_ELT(rdn, 1)); SET_DIMNAMES(res, dn); UNPROTECT(1); @@ -67,20 +55,16 @@ SEXP crowsum(SEXP Rmat, SEXP Rfactor, SEXP Rmean) // Now, run through column by column, summing the levels // #pragma omp parallel for num_threads(nthr) memset(mres, 0, nlev * cols * sizeof(double)); - mres--; // factor levels are 1-based - for (int k = 0; k < cols; k++, mres += nlev) - { - for (int i = 0; i < rows; i++) - { + mres--; // factor levels are 1-based + for (int k = 0; k < cols; k++, mres += nlev) { + for (int i = 0; i < rows; i++) { mres[f[i]] += *mat++; } } - if (mean) - { + if (mean) { mres = REAL(res); for (int k = 0; k < cols; k++, mres += nlev) - for (int i = 0; i < nlev; i++) - mres[i] /= table[i]; + for (int i = 0; i < nlev; i++) mres[i] /= table[i]; } UNPROTECT(1); return (res); diff --git a/src/demean.c b/src/demean.c index b04ec5c..0e4e12e 100644 --- a/src/demean.c +++ b/src/demean.c @@ -1,6 +1,6 @@ #include "lfe.h" /* Need snprintf */ -#include +#include /* Argument block for threads. @@ -42,51 +42,48 @@ extern SEXP df_string; Centre on all factors in succession. Vector v. In place. */ - -static R_INLINE void centre(double *v, int N, - FACTOR *factors[], int e, double *means, - double *weights) { - +static R_INLINE void centre(double *v, int N, FACTOR *factors[], int e, + double *means, double *weights) { const int hw = (weights != NULL); - for(int i=0; i < e; i++) { + for (int i = 0; i < e; i++) { FACTOR *f = factors[i]; const int *gp = f->group; const int hx = (f->x != NULL); /* compute means */ - memset(means,0,sizeof(double)* f->nlevels); - for(int j = 0; j < N; j++) { - double w = hw ? (hx ? f->x[j]*weights[j] : weights[j]) : (hx ? f->x[j] : 1.0); - if(gp[j] > 0) means[gp[j]-1] += v[j]*w; + memset(means, 0, sizeof(double) * f->nlevels); + for (int j = 0; j < N; j++) { + double w = + hw ? (hx ? f->x[j] * weights[j] : weights[j]) : (hx ? f->x[j] : 1.0); + if (gp[j] > 0) means[gp[j] - 1] += v[j] * w; } #ifdef USEOMP - #pragma omp parallel for +#pragma omp parallel for #endif - for(int j = 0; j < f->nlevels; j++) { + for (int j = 0; j < f->nlevels; j++) { means[j] *= f->invgpsize[j]; } /* subtract means */ #ifdef USEOMP - #pragma omp parallel for +#pragma omp parallel for #endif - for(int j = 0; j < N; j++) { - double w = hw ? (hx ? f->x[j]*weights[j] : weights[j]) : (hx ? f->x[j] : 1.0); - if(gp[j] > 0) v[j] -= means[gp[j]-1]*w; + for (int j = 0; j < N; j++) { + double w = + hw ? (hx ? f->x[j] * weights[j] : weights[j]) : (hx ? f->x[j] : 1.0); + if (gp[j] > 0) v[j] -= means[gp[j] - 1] * w; } } } - /* Method of alternating projections. In place in vectors pointed to by vec. */ -static int demean(int N, double *vec, double *weights,int *scale, - FACTOR *factors[],int e,double eps, double *means, - double *prev, double *prev2, - int *stop,int vecnum, LOCK_T lock, const int gkacc, - const int quiet) { +static int demean(int N, double *vec, double *weights, int *scale, + FACTOR *factors[], int e, double eps, double *means, + double *prev, double *prev2, int *stop, int vecnum, + LOCK_T lock, const int gkacc, const int quiet) { double delta; double norm2 = 0.0; int i; @@ -103,146 +100,157 @@ static int demean(int N, double *vec, double *weights,int *scale, // for(int i = 0; i < N; i++) { // if(isnan(vec[i])) vec[i]=0.0; // } - if(weights != NULL && scale[0]) for(int i = 0; i < N; i++) vec[i] = vec[i]*weights[i]; + if (weights != NULL && scale[0]) + for (int i = 0; i < N; i++) vec[i] = vec[i] * weights[i]; centre(vec, N, factors, e, means, weights); - if(e <= 1 || factors[0]->oneiter) { - if(weights != NULL && scale[1]) for(int i = 0; i < N; i++) vec[i] = vec[i]/weights[i]; - return(1); + if (e <= 1 || factors[0]->oneiter) { + if (weights != NULL && scale[1]) + for (int i = 0; i < N; i++) vec[i] = vec[i] / weights[i]; + return (1); } // Initialize - memcpy(prev2, vec, N*sizeof(double)); + memcpy(prev2, vec, N * sizeof(double)); norm2 = 1.0; - for(i = 0; i < N; i++) norm2 += vec[i]*vec[i]; + for (i = 0; i < N; i++) norm2 += vec[i] * vec[i]; norm2 = sqrt(norm2); - delta = prevdelta = 2*norm2; - neps = norm2*eps; + delta = prevdelta = 2 * norm2; + neps = norm2 * eps; iter = 0; last = time(NULL); lastiter = 0; - double target=0; - target = 1e-4*neps; - double prevdp=0.0; + double target = 0; + target = 1e-4 * neps; + double prevdp = 0.0; do { iter++; - if(gkacc) memcpy(prev, vec, N*sizeof(double)); - centre(vec,N,factors,e,means, weights); + if (gkacc) memcpy(prev, vec, N * sizeof(double)); + centre(vec, N, factors, e, means, weights); - if(gkacc) { + if (gkacc) { // Try the Gearhart-Koshy acceleration, AS 3.2 - // The result vector is an affine combination of the previous x and the centred vector Px - // tPx + (1-t) x - // where t = (x, x-Px)/|x - Px|^2 (t can be > 1). - // It's quite much faster for some datasets (factor 2-3), perhaps more. - // Though this seems to be related to lower precision. - double ip = 0, nm=0, nm2=0; - for(int i = 0; i < N; i++) { - nm += prev[i]*prev[i]; - ip += prev[i]*(prev[i]-vec[i]); - nm2 += (prev[i]-vec[i])*(prev[i] - vec[i]); + // The result vector is an affine combination of the previous x and the + // centred vector Px tPx + (1-t) x where t = (x, x-Px)/|x - Px|^2 (t can + // be > 1). It's quite much faster for some datasets (factor 2-3), perhaps + // more. Though this seems to be related to lower precision. + double ip = 0, nm = 0, nm2 = 0; + for (int i = 0; i < N; i++) { + nm += prev[i] * prev[i]; + ip += prev[i] * (prev[i] - vec[i]); + nm2 += (prev[i] - vec[i]) * (prev[i] - vec[i]); } - if(nm2 > 1e-18*nm) { - double t = ip/nm2; - - // By Lemma 3.9, t >= 0.5. If we're below zero we're numerically out of whack - // Or does it mean we have converged? - if(t < 0.49) { - char buf[256]; - if(nm2 > 1e-15*nm) { - snprintf(buf,sizeof(buf),"Demeaning of vec %d failed after %d iterations, t=%.1e, nm2=%.1enm\n", - vecnum,iter,t,nm2/nm); - pushmsg(buf,lock); - okconv = 0; - } - break; - } - - for(int i = 0; i < N; i++) { - vec[i] = t*vec[i] + (1.0-t)*prev[i]; - } + if (nm2 > 1e-18 * nm) { + double t = ip / nm2; + + // By Lemma 3.9, t >= 0.5. If we're below zero we're numerically out of + // whack Or does it mean we have converged? + if (t < 0.49) { + char buf[256]; + if (nm2 > 1e-15 * nm) { + snprintf(buf, sizeof(buf), + "Demeaning of vec %d failed after %d iterations, t=%.1e, " + "nm2=%.1enm\n", + vecnum, iter, t, nm2 / nm); + pushmsg(buf, lock); + okconv = 0; + } + break; + } + + for (int i = 0; i < N; i++) { + vec[i] = t * vec[i] + (1.0 - t) * prev[i]; + } } } delta = neps = 0.0; - for(int i = 0; i < N; i++) { - delta += (prev2[i]-vec[i])*(prev2[i]-vec[i]); - neps += vec[i]*vec[i]; + for (int i = 0; i < N; i++) { + delta += (prev2[i] - vec[i]) * (prev2[i] - vec[i]); + neps += vec[i] * vec[i]; } - memcpy(prev2,vec,N*sizeof(double)); + memcpy(prev2, vec, N * sizeof(double)); delta = sqrt(delta); - neps = sqrt(1.0+neps)*eps; - c = delta/prevdelta; - /* c is the convergence rate per iteration, at infinity they add up to 1/(1-c). - This is true for e==2, but also in the ballpark for e>2 we guess. - Only fail if it's after some rounds. It seems a bit unstable in the - beginning. This is to be expected due to Deutsch and Hundal's result + neps = sqrt(1.0 + neps) * eps; + c = delta / prevdelta; + /* c is the convergence rate per iteration, at infinity they add up to + 1/(1-c). This is true for e==2, but also in the ballpark for e>2 we + guess. Only fail if it's after some rounds. It seems a bit unstable in + the beginning. This is to be expected due to Deutsch and Hundal's result If we're doing acceleration, we test above for t < 0.5 for failure */ - if(!gkacc && c >= 1.0 && iter > 100) { + if (!gkacc && c >= 1.0 && iter > 100) { char buf[256]; - snprintf(buf,sizeof(buf),"Demeaning of vec %d failed after %d iterations, c-1=%.0e, delta=%.1e\n", - vecnum,iter,c-1.0,delta); - pushmsg(buf,lock); + snprintf(buf, sizeof(buf), + "Demeaning of vec %d failed after %d iterations, c-1=%.0e, " + "delta=%.1e\n", + vecnum, iter, c - 1.0, delta); + pushmsg(buf, lock); okconv = 0; break; } - if(c < 1.0) target = (1.0-c)*neps; - // target should not be smaller than 3 bits of precision in each coordinate - if(target < N*1e-15) target = N*1e-15; + if (c < 1.0) target = (1.0 - c) * neps; + // target should not be smaller than 3 bits of precision in each coordinate + if (target < N * 1e-15) target = N * 1e-15; now = time(NULL); // should we report? - if(quiet > 0 && now - last >= quiet && delta > 0.0) { + if (quiet > 0 && now - last >= quiet && delta > 0.0) { int reqiter; double eta; char buf[256]; - double targ; - if(prevdp != 0.0) { - // use an average c - c = pow(delta/prevdp, 1.0/(iter-lastiter)); - targ = (1.0-c)*neps; - } else targ = target; - - reqiter = log(targ/delta)/log(c); - eta = 1.0*(now-last)*reqiter/(iter-lastiter); - if(eta < 0) eta = 0; - - time_t arriv = now + (time_t) (eta+1.0); + double targ; + if (prevdp != 0.0) { + // use an average c + c = pow(delta / prevdp, 1.0 / (iter - lastiter)); + targ = (1.0 - c) * neps; + } else + targ = target; + + reqiter = log(targ / delta) / log(c); + eta = 1.0 * (now - last) * reqiter / (iter - lastiter); + if (eta < 0) eta = 0; + + time_t arriv = now + (time_t)(eta + 1.0); char timbuf[50]; struct tm tmarriv; #ifdef WIN localtime_s(&tmarriv, &arriv); #else - localtime_r(&arriv,&tmarriv); + localtime_r(&arriv, &tmarriv); #endif strftime(timbuf, sizeof(timbuf), "%c", &tmarriv); - snprintf(buf,sizeof(buf),"...centering vec %d i:%d c:%.1e d:%.1e(t:%.1e) ETA:%s\n", - vecnum,iter,1.0-c,delta,target,timbuf); - pushmsg(buf,lock); + snprintf(buf, sizeof(buf), + "...centering vec %d i:%d c:%.1e d:%.1e(t:%.1e) ETA:%s\n", + vecnum, iter, 1.0 - c, delta, target, timbuf); + pushmsg(buf, lock); lastiter = iter; prevdp = delta; last = now; } - - prevdelta = delta; - + + prevdelta = delta; + #ifdef NOTHREADS R_CheckUserInterrupt(); #else int stp; - LOCK(lock); stp = *stop; UNLOCK(lock); - if(stp != 0) {okconv = 0; break;} + LOCK(lock); + stp = *stop; + UNLOCK(lock); + if (stp != 0) { + okconv = 0; + break; + } #endif - - } while(delta > target); - if(weights != NULL && scale[1]) for(int i = 0; i < N; i++) vec[i] = vec[i]/weights[i]; - return(okconv); -} - + } while (delta > target); + if (weights != NULL && scale[1]) + for (int i = 0; i < N; i++) vec[i] = vec[i] / weights[i]; + return (okconv); +} /* Thread routine for demeaning. @@ -251,8 +259,7 @@ static int demean(int N, double *vec, double *weights,int *scale, by a mutex. */ - -/* The thread routine +/* The thread routine Each thread centres an entire vector. */ #ifdef WIN @@ -260,7 +267,7 @@ DWORD WINAPI demeanlist_thr(LPVOID varg) { #else static void *demeanlist_thr(void *varg) { #endif - PTARG *arg = (PTARG*) varg; + PTARG *arg = (PTARG *)varg; double *means; double *tmp1, *tmp2; int vecnum; @@ -273,42 +280,43 @@ static void *demeanlist_thr(void *varg) { means = arg->means[myid]; tmp1 = arg->tmp1[myid]; tmp2 = arg->tmp2[myid]; - while(1) { + while (1) { time_t now; /* Find next vector */ LOCK(arg->lock); vecnum = (arg->nowdoing)++; UNLOCK(arg->lock); - if(vecnum >= arg->K) break; + if (vecnum >= arg->K) break; #ifdef HAVE_THREADNAME #ifndef USEOMP char thrname[16]; char buf[256]; - snprintf(buf,256,"Ct %d/%d",vecnum+1, arg->K); + snprintf(buf, 256, "Ct %d/%d", vecnum + 1, arg->K); buf[15] = 0; - memcpy(thrname,buf,16); + memcpy(thrname, buf, 16); STNAME(thrname); #endif #endif // If in place, res is pointing to v, otherwise we copy v to res - if(arg->res[vecnum] != arg->v[vecnum]) memcpy(arg->res[vecnum],arg->v[vecnum],arg->N*sizeof(double)); - okconv = demean(arg->N,arg->res[vecnum],arg->weights,arg->scale, - arg->factors,arg->e,arg->eps,means, - tmp1,tmp2, - &arg->stop,vecnum+1,arg->lock,arg->gkacc, arg->quiet); + if (arg->res[vecnum] != arg->v[vecnum]) + memcpy(arg->res[vecnum], arg->v[vecnum], arg->N * sizeof(double)); + okconv = demean(arg->N, arg->res[vecnum], arg->weights, arg->scale, + arg->factors, arg->e, arg->eps, means, tmp1, tmp2, + &arg->stop, vecnum + 1, arg->lock, arg->gkacc, arg->quiet); now = time(NULL); LOCK(arg->lock); - if(!okconv) { + if (!okconv) { arg->badconv++; } (arg->done)++; - if(arg->quiet > 0 && now > arg->last + arg->quiet && arg->K > 1) { + if (arg->quiet > 0 && now > arg->last + arg->quiet && arg->K > 1) { char buf[256]; - snprintf(buf,sizeof(buf),"...finished centering %d of %d vectors in %d seconds\n", - arg->done,arg->K,(int)(now-arg->start)); + snprintf(buf, sizeof(buf), + "...finished centering %d of %d vectors in %d seconds\n", + arg->done, arg->K, (int)(now - arg->start)); arg->last = now; - UNLOCK(arg->lock); - pushmsg(buf,arg->lock); + UNLOCK(arg->lock); + pushmsg(buf, arg->lock); } else { UNLOCK(arg->lock); } @@ -321,7 +329,7 @@ static void *demeanlist_thr(void *varg) { Last leaver turns out the light.*/ arg->running--; #ifdef HAVE_SEM - if(arg->running == 0) sem_post(&arg->finished); + if (arg->running == 0) sem_post(&arg->finished); #endif UNLOCK(arg->lock); #endif @@ -329,11 +337,10 @@ static void *demeanlist_thr(void *varg) { return 0; } - /* main C entry-point for demeaning. Called from the R-entrypoint */ -static int demeanlist(double **vp, int N, int K, double **res, double *weights,int *scale, - FACTOR *factors[], int e, double eps,int cores, - int quiet, const int gkacc) { +static int demeanlist(double **vp, int N, int K, double **res, double *weights, + int *scale, FACTOR *factors[], int e, double eps, + int cores, int quiet, const int gkacc) { PTARG arg; int numthr = 1; int thr; @@ -355,39 +362,40 @@ static int demeanlist(double **vp, int N, int K, double **res, double *weights,i #endif numthr = cores; - if(numthr > K) numthr = K; - if(numthr < 1) numthr = 1; + if (numthr > K) numthr = K; + if (numthr < 1) numthr = 1; #ifdef WIN - lock = CreateMutex(NULL,FALSE,NULL); - if(lock == NULL) { + lock = CreateMutex(NULL, FALSE, NULL); + if (lock == NULL) { error("Couldn't create mutex (error=%d)", GetLastError()); } - threads = (HANDLE*) R_alloc(numthr,sizeof(HANDLE)); - threadids = (DWORD*) R_alloc(numthr,sizeof(DWORD)); + threads = (HANDLE *)R_alloc(numthr, sizeof(HANDLE)); + threadids = (DWORD *)R_alloc(numthr, sizeof(DWORD)); arg.lock = lock; #else - threads = (pthread_t*) R_alloc(numthr,sizeof(pthread_t)); + threads = (pthread_t *)R_alloc(numthr, sizeof(pthread_t)); arg.lock = &lock; #ifdef HAVE_SEM - if(sem_init(&arg.finished,0,0) != 0) error("sem_init failed, errno=%d",errno); + if (sem_init(&arg.finished, 0, 0) != 0) + error("sem_init failed, errno=%d", errno); #endif arg.running = numthr; #endif #endif maxlev = 0; - for(i = 0; i < e; i++) - if(maxlev < factors[i]->nlevels) maxlev = factors[i]->nlevels; - - arg.means = (double **) R_alloc(numthr,sizeof(double*)); - arg.tmp1 = (double **) R_alloc(numthr,sizeof(double*)); - arg.tmp2 = (double **) R_alloc(numthr,sizeof(double*)); - for(thr = 0; thr < numthr; thr++) { - arg.means[thr] = (double*) R_alloc(maxlev,sizeof(double)); - if(e > 1) { - arg.tmp1[thr] = (double*) R_alloc(N,sizeof(double)); - arg.tmp2[thr] = (double*) R_alloc(N,sizeof(double)); + for (i = 0; i < e; i++) + if (maxlev < factors[i]->nlevels) maxlev = factors[i]->nlevels; + + arg.means = (double **)R_alloc(numthr, sizeof(double *)); + arg.tmp1 = (double **)R_alloc(numthr, sizeof(double *)); + arg.tmp2 = (double **)R_alloc(numthr, sizeof(double *)); + for (thr = 0; thr < numthr; thr++) { + arg.means[thr] = (double *)R_alloc(maxlev, sizeof(double)); + if (e > 1) { + arg.tmp1[thr] = (double *)R_alloc(N, sizeof(double)); + arg.tmp2[thr] = (double *)R_alloc(N, sizeof(double)); } } @@ -411,52 +419,53 @@ static int demeanlist(double **vp, int N, int K, double **res, double *weights,i arg.gkacc = gkacc; #ifdef NOTHREADS - demeanlist_thr((void*)&arg); + demeanlist_thr((void *)&arg); #else /* Do it in separate threads */ - for(thr = 0; thr < numthr; thr++) { + for (thr = 0; thr < numthr; thr++) { #ifdef WIN - threads[thr] = CreateThread(NULL,0,demeanlist_thr,&arg,0,&threadids[thr]); - if(0 == threads[thr]) error("Failed to create thread"); + threads[thr] = + CreateThread(NULL, 0, demeanlist_thr, &arg, 0, &threadids[thr]); + if (0 == threads[thr]) error("Failed to create thread"); #else - int stat = pthread_create(&threads[thr],NULL, demeanlist_thr,&arg); - if(0 != stat) error("Failed to create thread, stat=%d",stat); + int stat = pthread_create(&threads[thr], NULL, demeanlist_thr, &arg); + if (0 != stat) error("Failed to create thread, stat=%d", stat); #endif } /* wait for completion */ /* We want to check for interrupts regularly, and set a stop flag */ - while(1) { + while (1) { printmsg(arg.lock); - if(arg.stop == 0 && checkInterrupt()) { + if (arg.stop == 0 && checkInterrupt()) { REprintf("...stopping centering threads...\n"); - arg.stop=1; + arg.stop = 1; } #ifdef WIN - if(WaitForMultipleObjects(numthr,threads,TRUE,3000) != WAIT_TIMEOUT) { - for(thr = 0; thr < numthr; thr++) { - CloseHandle(threads[thr]); + if (WaitForMultipleObjects(numthr, threads, TRUE, 3000) != WAIT_TIMEOUT) { + for (thr = 0; thr < numthr; thr++) { + CloseHandle(threads[thr]); } break; } #else #ifndef HAVE_SEM - struct timespec atmo = {0,50000000}; + struct timespec atmo = {0, 50000000}; /* Kludge when no timedwait, i.e. MacOSX */ int done; - if(arg.stop == 0) nanosleep(&atmo,NULL); + if (arg.stop == 0) nanosleep(&atmo, NULL); LOCK(arg.lock); done = (arg.running == 0); UNLOCK(arg.lock); - if(arg.stop == 1 || done) { + if (arg.stop == 1 || done) { #else - struct timespec tmo = {time(NULL)+3,0}; - if(arg.stop == 1 || sem_timedwait(&arg.finished,&tmo) == 0) { + struct timespec tmo = {time(NULL) + 3, 0}; + if (arg.stop == 1 || sem_timedwait(&arg.finished, &tmo) == 0) { #endif - for(thr = 0; thr < numthr; thr++) { - (void)pthread_join(threads[thr], NULL); + for (thr = 0; thr < numthr; thr++) { + (void)pthread_join(threads[thr], NULL); } break; } @@ -464,7 +473,7 @@ static int demeanlist(double **vp, int N, int K, double **res, double *weights,i } #endif - /* print remaining messages */ + /* print remaining messages */ printmsg(arg.lock); /* Release synchronization gear */ #ifdef WIN @@ -476,31 +485,32 @@ static int demeanlist(double **vp, int N, int K, double **res, double *weights,i pthread_mutex_destroy(arg.lock); #endif + if (arg.stop == 1) error("centering interrupted"); + if (quiet > 0 && arg.start != arg.last && K > 1) + REprintf("...%d vectors centred in %d seconds\n", K, + (int)(time(NULL) - arg.start)); - if(arg.stop == 1) error("centering interrupted"); - if(quiet > 0 && arg.start != arg.last && K > 1) - REprintf("...%d vectors centred in %d seconds\n",K,(int)(time(NULL)-arg.start)); - - return(arg.badconv); + return (arg.badconv); } - SEXP inplace(SEXP x) { - setAttrib(x,install("LFE_inplace"),R_QuoteSymbol); // Just set it to something non-NULL + setAttrib(x, install("LFE_inplace"), + R_QuoteSymbol); // Just set it to something non-NULL return x; } SEXP outplace(SEXP x) { - setAttrib(x,install("LFE_inplace"),R_NilValue); // Just set it to something non-NULL + setAttrib(x, install("LFE_inplace"), + R_NilValue); // Just set it to something non-NULL return x; } int isinplace(SEXP x) { - return getAttrib(x,install("LFE_inplace")) != R_NilValue; + return getAttrib(x, install("LFE_inplace")) != R_NilValue; } /* Now, for the R-interface. We only export a demeanlist function - It should be called from R with + It should be called from R with demeanlist(matrix,list(f1,f2,...),double eps) and preferably update matrix inplace to save copying memory unless that's completely un-R'ish. Hmm, it is. @@ -510,9 +520,9 @@ int isinplace(SEXP x) { and call our demeanlist and return */ -SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, - SEXP scores, SEXP quiet, SEXP gkacc, SEXP Rmeans, - SEXP Rweights, SEXP Rscale, SEXP attrs) { +SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, SEXP scores, + SEXP quiet, SEXP gkacc, SEXP Rmeans, SEXP Rweights, + SEXP Rscale, SEXP attrs) { int numvec; int numfac; int cnt; @@ -524,43 +534,44 @@ SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, FACTOR **factors; int cores; int i; - int N=0; + int N = 0; int icpt; int *vicpt; int icptlen; SEXP badconv; int listlen; - int domeans=0; - int protectcount=0; + int domeans = 0; + int protectcount = 0; int scale[2]; int wraplist = 0; int inplacelist = 0; - + // Find the length of the data // We are never called with length(flist) == 0 - if(!isNewList(flist)) error("flist is not a list"); - if(LENGTH(flist) == 0) { + if (!isNewList(flist)) error("flist is not a list"); + if (LENGTH(flist) == 0) { warning("demeanlist called with length(fl)==0, internal error?"); N = 0; } else { - N = LENGTH(VECTOR_ELT(flist,0)); + N = LENGTH(VECTOR_ELT(flist, 0)); } - if(isinplace(vlist) || (!isNewList(vlist) && NO_REFERENCES(vlist))) { + if (isinplace(vlist) || (!isNewList(vlist) && NO_REFERENCES(vlist))) { inplacelist = 1; outplace(vlist); } - if(LENGTH(Rscale) == 1) { + if (LENGTH(Rscale) == 1) { scale[0] = LOGICAL(AS_LOGICAL(Rscale))[0]; scale[1] = scale[0]; - } else if(LENGTH(Rscale) > 1) { + } else if (LENGTH(Rscale) > 1) { scale[0] = LOGICAL(AS_LOGICAL(Rscale))[0]; scale[1] = LOGICAL(AS_LOGICAL(Rscale))[1]; } else { error("scale must have length > 0"); } domeans = LOGICAL(AS_LOGICAL(Rmeans))[0]; - vicpt = INTEGER(PROTECT(AS_INTEGER(Ricpt)));protectcount++; + vicpt = INTEGER(PROTECT(AS_INTEGER(Ricpt))); + protectcount++; icptlen = LENGTH(Ricpt); icpt = vicpt[0] - 1; /* convert from 1-based to zero-based */ eps = REAL(AS_NUMERIC(Reps))[0]; @@ -568,105 +579,115 @@ SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, cores = 1; #else cores = INTEGER(AS_INTEGER(scores))[0]; - if(cores < 1) cores = 1; + if (cores < 1) cores = 1; #endif - if(!isNull(Rweights)) { - if(LENGTH(Rweights) != N) error("Length of weights (%d) must equal length of data (%d)", - LENGTH(Rweights),N); - weights = REAL(PROTECT(AS_NUMERIC(Rweights))); protectcount++; + if (!isNull(Rweights)) { + if (LENGTH(Rweights) != N) + error("Length of weights (%d) must equal length of data (%d)", + LENGTH(Rweights), N); + weights = REAL(PROTECT(AS_NUMERIC(Rweights))); + protectcount++; } factors = makefactors(flist, 1, weights); numfac = 0; - for(FACTOR **f = factors; *f != NULL; f++) numfac++; + for (FACTOR **f = factors; *f != NULL; f++) numfac++; - int isdf = inherits(vlist,"data.frame"); - if(isdf) { + int isdf = inherits(vlist, "data.frame"); + if (isdf) { icptlen = 1; icpt = -1; } - if(!isNewList(vlist)) { + if (!isNewList(vlist)) { // Just put it in a list to avoid tests further down. - SEXP tmp = PROTECT(allocVector(VECSXP,1)); + SEXP tmp = PROTECT(allocVector(VECSXP, 1)); protectcount++; - SET_VECTOR_ELT(tmp,0,vlist); + SET_VECTOR_ELT(tmp, 0, vlist); vlist = tmp; wraplist = 1; } - PROTECT(vlist = AS_LIST(vlist)); protectcount++; + PROTECT(vlist = AS_LIST(vlist)); + protectcount++; - listlen = LENGTH(vlist); - if(icptlen != 1 && icptlen != listlen) - error("Length of icpt (%d) should be 1 or the number of arguments (%d)", icptlen, listlen); + listlen = LENGTH(vlist); + if (icptlen != 1 && icptlen != listlen) + error("Length of icpt (%d) should be 1 or the number of arguments (%d)", + icptlen, listlen); - PROTECT(reslist = NEW_LIST(listlen)); protectcount++; + PROTECT(reslist = NEW_LIST(listlen)); + protectcount++; SET_NAMES(reslist, GET_NAMES(vlist)); - if(isdf) { - if(inherits(vlist, "pdata.frame")) { + if (isdf) { + if (inherits(vlist, "pdata.frame")) { DUPLICATE_ATTRIB(reslist, vlist); } else { setAttrib(reslist, R_RowNamesSymbol, getAttrib(vlist, R_RowNamesSymbol)); - classgets(reslist,df_string); + classgets(reslist, df_string); } } /* First, count the number of vectors in total */ numvec = 0; - for(int i = 0; i < listlen; i++) { - SEXP elt = VECTOR_ELT(vlist,i); + for (int i = 0; i < listlen; i++) { + SEXP elt = VECTOR_ELT(vlist, i); /* Each entry in the list is either a vector or a matrix */ - if(isNull(elt)) continue; - if(!isMatrix(elt)) { - if(LENGTH(elt) != N) - error("mtx[%d]: Factors and vectors must have the same length %d != %d",i, N,LENGTH(elt)); + if (isNull(elt)) continue; + if (!isMatrix(elt)) { + if (LENGTH(elt) != N) + error("mtx[%d]: Factors and vectors must have the same length %d != %d", + i, N, LENGTH(elt)); numvec++; } else { - if(nrows(elt) != N) - error("mtx[%d]: Factors and vectors must have the same length %d != %d",i, N,nrows(elt)); + if (nrows(elt) != N) + error("mtx[%d]: Factors and vectors must have the same length %d != %d", + i, N, nrows(elt)); numvec += ncols(elt); - if(icptlen != 1) icpt = vicpt[i]-1; - if(icpt >= 0 && icpt < ncols(elt)) numvec--; + if (icptlen != 1) icpt = vicpt[i] - 1; + if (icpt >= 0 && icpt < ncols(elt)) numvec--; } } /* Allocate pointers to source vectors */ - vectors = (double **)R_alloc(numvec,sizeof(double*)); + vectors = (double **)R_alloc(numvec, sizeof(double *)); /* Allocate pointers to result vectors */ - target = (double**) R_alloc(numvec,sizeof(double*)); + target = (double **)R_alloc(numvec, sizeof(double *)); /* Loop through list again to set up result structure */ cnt = 0; - for(i = 0; i < listlen; i++) { - SEXP elt = VECTOR_ELT(vlist,i); + for (i = 0; i < listlen; i++) { + SEXP elt = VECTOR_ELT(vlist, i); int inplaceelt = 0; - if(isinplace(elt)) { + if (isinplace(elt)) { outplace(elt); inplaceelt = 1; } - if(isNull(elt)) continue; + if (isNull(elt)) continue; - if(!isReal(elt)) { - PROTECT(elt = coerceVector(elt, REALSXP)); protectcount++; + if (!isReal(elt)) { + PROTECT(elt = coerceVector(elt, REALSXP)); + protectcount++; inplaceelt = 1; } - if(!isMatrix(elt)) { + if (!isMatrix(elt)) { /* It's a vector */ SEXP resvec; vectors[cnt] = REAL(elt); - if(!domeans && (inplacelist || inplaceelt) ) { - // in place centering - SET_VECTOR_ELT(reslist,i,elt); - target[cnt] = vectors[cnt]; + if (!domeans && (inplacelist || inplaceelt)) { + // in place centering + SET_VECTOR_ELT(reslist, i, elt); + target[cnt] = vectors[cnt]; } else { - if(LENGTH(elt) != N) error("lengths of vectors (%d) != lengths of factors(%d)",LENGTH(elt),N); - PROTECT(resvec = allocVector(REALSXP,LENGTH(elt))); - SET_VECTOR_ELT(reslist,i,resvec); - UNPROTECT(1); - SET_NAMES(resvec, GET_NAMES(elt)); - target[cnt] = REAL(resvec); + if (LENGTH(elt) != N) + error("lengths of vectors (%d) != lengths of factors(%d)", + LENGTH(elt), N); + PROTECT(resvec = allocVector(REALSXP, LENGTH(elt))); + SET_VECTOR_ELT(reslist, i, resvec); + UNPROTECT(1); + SET_NAMES(resvec, GET_NAMES(elt)); + target[cnt] = REAL(resvec); } cnt++; } else { @@ -675,69 +696,72 @@ SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, int rcols = cols; int j; SEXP mtx; - if(nrows(elt) != N) error("nrow(matrix) (%d) != lengths of factors(%d)",nrows(elt),N); + if (nrows(elt) != N) + error("nrow(matrix) (%d) != lengths of factors(%d)", nrows(elt), N); /* Allocate a matrix */ - if(icptlen != 1) icpt = vicpt[i]-1; - if(icpt >= 0 && icpt < cols) rcols--; - - if(!domeans && (inplacelist || inplaceelt) ) { - // in place centering - SET_VECTOR_ELT(reslist,i,elt); - for(int j=0; j < cols; j++) { - vectors[cnt] = target[cnt] = REAL(elt) + j*(mybigint_t)N; - cnt++; - } + if (icptlen != 1) icpt = vicpt[i] - 1; + if (icpt >= 0 && icpt < cols) rcols--; + + if (!domeans && (inplacelist || inplaceelt)) { + // in place centering + SET_VECTOR_ELT(reslist, i, elt); + for (int j = 0; j < cols; j++) { + vectors[cnt] = target[cnt] = REAL(elt) + j * (mybigint_t)N; + cnt++; + } } else { - PROTECT(mtx = allocMatrix(REALSXP,N,rcols)); - SET_VECTOR_ELT(reslist,i,mtx); - UNPROTECT(1); - // copy colnames - if(cols == rcols) { - SET_DIMNAMES(mtx, GET_DIMNAMES(elt)); - } - /* Set up pointers */ - rcols = 0; - for(j = 0; j < cols; j++) { - if(j == icpt) continue; - vectors[cnt] = REAL(elt) + j * (mybigint_t)N; - target[cnt] = REAL(mtx) + (rcols++) * (mybigint_t)N; - cnt++; - } + PROTECT(mtx = allocMatrix(REALSXP, N, rcols)); + SET_VECTOR_ELT(reslist, i, mtx); + UNPROTECT(1); + // copy colnames + if (cols == rcols) { + SET_DIMNAMES(mtx, GET_DIMNAMES(elt)); + } + /* Set up pointers */ + rcols = 0; + for (j = 0; j < cols; j++) { + if (j == icpt) continue; + vectors[cnt] = REAL(elt) + j * (mybigint_t)N; + target[cnt] = REAL(mtx) + (rcols++) * (mybigint_t)N; + cnt++; + } } } } /* Then do stuff */ - PROTECT(badconv = allocVector(INTSXP,1));protectcount++; - INTEGER(badconv)[0] = demeanlist(vectors,N,numvec,target,weights,scale, - factors,numfac, - eps,cores,INTEGER(AS_INTEGER(quiet))[0], - INTEGER(AS_INTEGER(gkacc))[0]); - - if(INTEGER(badconv)[0] > 0) - warning("%d vectors failed to centre to tolerance %.1le",INTEGER(badconv)[0],eps); + PROTECT(badconv = allocVector(INTSXP, 1)); + protectcount++; + INTEGER(badconv) + [0] = demeanlist(vectors, N, numvec, target, weights, scale, factors, numfac, + eps, cores, INTEGER(AS_INTEGER(quiet))[0], + INTEGER(AS_INTEGER(gkacc))[0]); + if (INTEGER(badconv)[0] > 0) + warning("%d vectors failed to centre to tolerance %.1le", + INTEGER(badconv)[0], eps); - if(domeans) { - for(int i = 0; i < numvec; i++) { + if (domeans) { + for (int i = 0; i < numvec; i++) { double *srcvec = vectors[i]; double *dstvec = target[i]; - for(int j = 0; j < N; j++) { - dstvec[j] = (srcvec[j]-dstvec[j]); + for (int j = 0; j < N; j++) { + dstvec[j] = (srcvec[j] - dstvec[j]); } } } SEXP ret = reslist; - if(wraplist) { - ret = VECTOR_ELT(reslist,0); + if (wraplist) { + ret = VECTOR_ELT(reslist, 0); } - if(INTEGER(badconv)[0] > 0) setAttrib(ret,install("badconv"),badconv); + if (INTEGER(badconv)[0] > 0) setAttrib(ret, install("badconv"), badconv); // Now, there may be more attributes to set - if(!isNull(attrs) && LENGTH(attrs) > 0) { - SEXP nm = PROTECT(GET_NAMES(attrs)); protectcount++; - for(int i = 0; i < LENGTH(attrs); i++) { - setAttrib(ret, installChar(STRING_ELT(nm,i)), VECTOR_ELT(attrs,i)); + if (!isNull(attrs) && LENGTH(attrs) > 0) { + SEXP nm = PROTECT(GET_NAMES(attrs)); + protectcount++; + for (int i = 0; i < LENGTH(attrs); i++) { + setAttrib(ret, installChar(STRING_ELT(nm, i)), VECTOR_ELT(attrs, i)); } } @@ -746,4 +770,3 @@ SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, return ret; } - diff --git a/src/factor.c b/src/factor.c index 659879f..32c7355 100644 --- a/src/factor.c +++ b/src/factor.c @@ -1,6 +1,5 @@ #include "lfe.h" -static void invertfactor(FACTOR *f, int N) -{ +static void invertfactor(FACTOR *f, int N) { int nlev = f->nlevels; int *curoff; int i; @@ -10,23 +9,19 @@ static void invertfactor(FACTOR *f, int N) memset(f->ii, 0, sizeof(int) * (nlev + 1)); /* find sizes of groups */ - for (i = 0; i < N; i++) - { + for (i = 0; i < N; i++) { int gp = f->group[i]; - if (gp < 1) - error("Factors can not have missing levels"); + if (gp < 1) error("Factors can not have missing levels"); f->ii[gp]++; } /* cumulative */ - for (i = 1; i < nlev + 1; i++) - { + for (i = 1; i < nlev + 1; i++) { f->ii[i] += f->ii[i - 1]; } curoff = Calloc(nlev + 1, int); - for (i = 0; i < N; i++) - { + for (i = 0; i < N; i++) { int gp = f->group[i] - 1; f->gpl[f->ii[gp] + curoff[gp]] = i; curoff[gp]++; @@ -34,96 +29,77 @@ static void invertfactor(FACTOR *f, int N) Free(curoff); } -FACTOR **makefactors(SEXP flist, int allowmissing, double *weights) -{ +FACTOR **makefactors(SEXP flist, int allowmissing, double *weights) { FACTOR **factors; int numfac = LENGTH(flist); int N = 0; int oneiter = 0; numfac = 0; - for (int i = 0; i < LENGTH(flist); i++) - { + for (int i = 0; i < LENGTH(flist); i++) { SEXP sf = VECTOR_ELT(flist, i); SEXP xattr = getAttrib(sf, install("x")); - if (isNull(xattr)) - { + if (isNull(xattr)) { numfac++; continue; - } - else if (LENGTH(flist) == 1) - { + } else if (LENGTH(flist) == 1) { SEXP ortho = getAttrib(xattr, install("ortho")); - if (isLogical(ortho)) - oneiter = LOGICAL(ortho)[0]; + if (isLogical(ortho)) oneiter = LOGICAL(ortho)[0]; } - if (!isMatrix(xattr)) - { + if (!isMatrix(xattr)) { numfac++; continue; } numfac += ncols(xattr); } - if (!oneiter) - { + if (!oneiter) { SEXP Roneiter = getAttrib(flist, install("oneiter")); - if (isLogical(Roneiter)) - oneiter = LOGICAL(Roneiter)[0]; + if (isLogical(Roneiter)) oneiter = LOGICAL(Roneiter)[0]; } factors = (FACTOR **)R_alloc(numfac + 1, sizeof(FACTOR *)); factors[numfac] = NULL; int truefac = 0; - for (int i = 0; i < LENGTH(flist); i++) - { + for (int i = 0; i < LENGTH(flist); i++) { int len; FACTOR *f; len = LENGTH(VECTOR_ELT(flist, i)); - if (i == 0) - { + if (i == 0) { N = len; - } - else if (len != N) - { + } else if (len != N) { error("All factors must have the same length %d %d", len, N); } f = factors[truefac++] = (FACTOR *)R_alloc(1, sizeof(FACTOR)); f->group = INTEGER(VECTOR_ELT(flist, i)); f->nlevels = LENGTH(getAttrib(VECTOR_ELT(flist, i), R_LevelsSymbol)); - if (f->nlevels <= 0) - error("factor %d in list has no levels\n", i + 1); + if (f->nlevels <= 0) error("factor %d in list has no levels\n", i + 1); f->oneiter = oneiter; SEXP xattr = getAttrib(VECTOR_ELT(flist, i), install("x")); - if (isNull(xattr)) - { + if (isNull(xattr)) { f->x = NULL; - } - else - { - if (isMatrix(xattr)) - { - if (nrows(xattr) != len) - { - error("Factor interaction terms (%d) must have the same length (%d) as the factor", - LENGTH(xattr), len); + } else { + if (isMatrix(xattr)) { + if (nrows(xattr) != len) { + error( + "Factor interaction terms (%d) must have the same length (%d) as " + "the factor", + LENGTH(xattr), len); } truefac--; - for (int j = 0; j < ncols(xattr); j++) - { + for (int j = 0; j < ncols(xattr); j++) { FACTOR *g = factors[truefac++] = (FACTOR *)R_alloc(1, sizeof(FACTOR)); g->group = f->group; g->nlevels = f->nlevels; g->oneiter = f->oneiter; g->x = &REAL(xattr)[j * (mybigint_t)nrows(xattr)]; } - } - else - { - if (LENGTH(xattr) != len) - { - error("Factor interaction terms (%d) must have the same length (%d) as the factor", - LENGTH(xattr), len); + } else { + if (LENGTH(xattr) != len) { + error( + "Factor interaction terms (%d) must have the same length (%d) as " + "the factor", + LENGTH(xattr), len); } f->x = REAL(xattr); } @@ -135,32 +111,27 @@ FACTOR **makefactors(SEXP flist, int allowmissing, double *weights) I.e. NA entries, how is that handled, I wonder. seems to be negative. Anyway, we fail on them. No, we don't */ - for (int i = 0; i < truefac; i++) - { + for (int i = 0; i < truefac; i++) { FACTOR *f = factors[i]; f->gpsize = (double *)R_alloc(f->nlevels, sizeof(double)); f->invgpsize = (double *)R_alloc(f->nlevels, sizeof(double)); memset(f->gpsize, 0, f->nlevels * sizeof(double)); /* first count it */ - for (int j = 0; j < N; j++) - { + for (int j = 0; j < N; j++) { /* skip entries without a group, do we need that? */ /* if(f->group[j] < 1) error("Factors can't have missing levels"); */ - if (f->group[j] > 0) - { - double w = (f->x == NULL) ? (weights == NULL ? 1.0 : weights[j]) : (weights == NULL ? f->x[j] : f->x[j] * weights[j]); + if (f->group[j] > 0) { + double w = (f->x == NULL) + ? (weights == NULL ? 1.0 : weights[j]) + : (weights == NULL ? f->x[j] : f->x[j] * weights[j]); f->gpsize[f->group[j] - 1] += w * w; - } - else - { - if (!allowmissing) - error("Factors can't have missing levels"); + } else { + if (!allowmissing) error("Factors can't have missing levels"); } } /* then invert it, it's much faster to multiply than to divide */ /* in the iterations */ - for (int j = 0; j < f->nlevels; j++) - { + for (int j = 0; j < f->nlevels; j++) { f->invgpsize[j] = 1.0 / f->gpsize[j]; } } @@ -172,8 +143,7 @@ FACTOR **makefactors(SEXP flist, int allowmissing, double *weights) but we would suffer from a deep call-stack. Hence, we make our own stack and push/pop in a loop */ -static int Components(int **vertices, FACTOR **factors, int K) -{ +static int Components(int **vertices, FACTOR **factors, int K) { int stacklevel = 0; int *stack; int curfac, curvert, curcomp, candvert; @@ -184,8 +154,7 @@ static int Components(int **vertices, FACTOR **factors, int K) /* How big a stack ? */ /* The number of vertices */ - for (i = 0; i < K; i++) - numvert += factors[i]->nlevels; + for (i = 0; i < K; i++) numvert += factors[i]->nlevels; /* Never used in threads, so use R's Calloc */ stack = Calloc(numvert * 4, int); @@ -210,27 +179,24 @@ static int Components(int **vertices, FACTOR **factors, int K) curcomp = 1; candvert = 0; - do - { + do { curvert = candvert; curfac = 0; /* Find the entire component */ - while (1) - { + while (1) { /* At the top here, curfac,curvert is a candidate for marking off as a vertex in curcomp For each iteration: If it's not marked, mark it, push on stack, go to first datapoint - If it's already marked, go to next datapoint for the vertex (incidence matrix) - If data-points are exhausted, go to next factor, start over with datapoints. - If factors are exhausted, pop the stack - If final stack-frame, we're done with component. + If it's already marked, go to next datapoint for the vertex (incidence + matrix) If data-points are exhausted, go to next factor, start over with + datapoints. If factors are exhausted, pop the stack If final stack-frame, + we're done with component. */ - if (vertices[curfac][curvert] == 0) - { + if (vertices[curfac][curvert] == 0) { /* Mark new vertex, find incidence list */ vertices[curfac][curvert] = curcomp; PUSHALL; @@ -238,27 +204,20 @@ static int Components(int **vertices, FACTOR **factors, int K) startfac = curfac; curfac = (startfac + 1) % K; ii = factors[startfac]->ii[startvert]; - } - else - { + } else { /* Been there, try next in group */ ii++; } - if (ii >= factors[startfac]->ii[startvert + 1]) - { + if (ii >= factors[startfac]->ii[startvert + 1]) { /* No more, move to next factor */ curfac = (curfac + 1) % K; - if (curfac == startfac) - { + if (curfac == startfac) { /* This is where we began, pop */ /* No more neighbours, go back to previous */ POPALL; /* Get outta here */ - if (0 == stacklevel) - break; - } - else - { + if (0 == stacklevel) break; + } else { /* start over with data-points */ ii = factors[startfac]->ii[startvert]; } @@ -281,11 +240,10 @@ static int Components(int **vertices, FACTOR **factors, int K) } // Algorithm from: -// "A Note on the Determination of Connectedness in an N-Way Cross Classification" -// D.L. Weeks and D.R. Williams, Technometrics, vol 6 no 3, August 1964 -// There probably exists faster algorithms, this one is quite slow -static void wwcomp(FACTOR *factors[], int numfac, int N, int *newlevels) -{ +// "A Note on the Determination of Connectedness in an N-Way Cross +// Classification" D.L. Weeks and D.R. Williams, Technometrics, vol 6 no 3, +// August 1964 There probably exists faster algorithms, this one is quite slow +static void wwcomp(FACTOR *factors[], int numfac, int N, int *newlevels) { int level = 0; int chead = 0; int *newlist = Calloc(N, int); @@ -294,18 +252,15 @@ static void wwcomp(FACTOR *factors[], int numfac, int N, int *newlevels) // For cache-efficiency, make a numfac x N matrix of the factors int *facmat = Calloc(numfac * N, int); - for (mysize_t i = 0; i < N; i++) - { + for (mysize_t i = 0; i < N; i++) { int *obsp = &facmat[i * numfac]; newlevels[i] = 0; oldlist[i] = i; - for (int j = 0; j < numfac; j++) - { + for (int j = 0; j < numfac; j++) { obsp[j] = factors[j]->group[i]; } } - while (oldstart < N) - { + while (oldstart < N) { int newidx; // set component number for first node in component // find the next we haven't checked, it's oldlist[oldstart] @@ -318,13 +273,12 @@ static void wwcomp(FACTOR *factors[], int numfac, int N, int *newlevels) newidx = 1; // loop over the list of newly added nodes, including the head // Note that we may increase newidx during the loop - for (int i = 0; i < newidx; i++) - { + for (int i = 0; i < newidx; i++) { mysize_t newnode = newlist[i]; int *newp = &facmat[newnode * numfac]; - // search for observations with distance 1 from newnode, mark them with level - for (int jidx = oldstart; jidx < N; jidx++) - { + // search for observations with distance 1 from newnode, mark them with + // level + for (int jidx = oldstart; jidx < N; jidx++) { mysize_t trynode = oldlist[jidx]; int *tryp = &facmat[trynode * numfac]; int dist = 0; @@ -334,8 +288,7 @@ static void wwcomp(FACTOR *factors[], int numfac, int N, int *newlevels) // dist += (factors[fi]->group[newnode] != factors[fi]->group[trynode]); // if close, set its level, add it to the list, move the start node // to the empty place in oldlist. - if (dist < 2) - { + if (dist < 2) { newlevels[trynode] = level; newlist[newidx++] = trynode; oldlist[jidx] = oldlist[oldstart++]; @@ -352,26 +305,22 @@ static void wwcomp(FACTOR *factors[], int numfac, int N, int *newlevels) R entry-point for conncomp. Takes a list of factors as input. */ -SEXP MY_wwcomp(SEXP flist) -{ +SEXP MY_wwcomp(SEXP flist) { int numfac, N; FACTOR **factors; SEXP result; numfac = LENGTH(flist); - if (numfac < 2) - error("At least two factors must be specified"); + if (numfac < 2) error("At least two factors must be specified"); N = LENGTH(VECTOR_ELT(flist, 0)); - for (int i = 0; i < numfac; i++) - { + for (int i = 0; i < numfac; i++) { if (N != LENGTH(VECTOR_ELT(flist, i))) error("Factors must have the same length"); } factors = (FACTOR **)R_alloc(numfac, sizeof(FACTOR *)); - for (int i = 0; i < numfac; i++) - { + for (int i = 0; i < numfac; i++) { FACTOR *f; factors[i] = (FACTOR *)R_alloc(1, sizeof(FACTOR)); @@ -381,26 +330,22 @@ SEXP MY_wwcomp(SEXP flist) PROTECT(result = allocVector(INTSXP, N)); int *fac = INTEGER(result); wwcomp(factors, numfac, N, fac); - // Now it's time to order the levels by decreasing size, so let's compute the sizes + // Now it's time to order the levels by decreasing size, so let's compute the + // sizes int levels = 0; for (int i = 0; i < N; i++) - if (fac[i] > levels) - levels = fac[i]; + if (fac[i] > levels) levels = fac[i]; double *levsize = (double *)R_alloc(levels, sizeof(double)); int *index = (int *)R_alloc(levels, sizeof(int)); - for (int i = 0; i < levels; i++) - { + for (int i = 0; i < levels; i++) { levsize[i] = 0.0; index[i] = i; } - for (int i = 0; i < N; i++) - levsize[fac[i] - 1] = levsize[fac[i] - 1] + 1; + for (int i = 0; i < N; i++) levsize[fac[i] - 1] = levsize[fac[i] - 1] + 1; revsort(levsize, index, levels); int *rindex = (int *)R_alloc(levels, sizeof(int)); - for (int i = 0; i < levels; i++) - rindex[index[i]] = i; - for (int i = 0; i < N; i++) - { + for (int i = 0; i < levels; i++) rindex[index[i]] = i; + for (int i = 0; i < N; i++) { fac[i] = rindex[fac[i] - 1] + 1; } @@ -414,8 +359,7 @@ a factor of the same length with the connection components */ -SEXP MY_conncomp(SEXP flist) -{ +SEXP MY_conncomp(SEXP flist) { int numfac; int i; @@ -431,18 +375,15 @@ SEXP MY_conncomp(SEXP flist) int *idx; numfac = LENGTH(flist); - if (numfac < 2) - error("At least two factors must be specified"); + if (numfac < 2) error("At least two factors must be specified"); N = LENGTH(VECTOR_ELT(flist, 0)); - for (i = 0; i < numfac; i++) - { + for (i = 0; i < numfac; i++) { if (N != LENGTH(VECTOR_ELT(flist, i))) error("Factors must have the same length"); } factors = (FACTOR **)R_alloc(numfac, sizeof(FACTOR *)); PROTECT(flist = AS_LIST(flist)); - for (i = 0; i < numfac; i++) - { + for (i = 0; i < numfac; i++) { FACTOR *f; factors[i] = (FACTOR *)R_alloc(1, sizeof(FACTOR)); @@ -457,8 +398,7 @@ SEXP MY_conncomp(SEXP flist) /* Create vertices */ vertices = (int **)R_alloc(numfac, sizeof(int *)); /* Create arrays for them */ - for (i = 0; i < numfac; i++) - { + for (i = 0; i < numfac; i++) { vertices[i] = (int *)R_alloc(factors[i]->nlevels, sizeof(int)); /* Assign no component to them*/ memset(vertices[i], 0, sizeof(int) * factors[i]->nlevels); @@ -471,8 +411,7 @@ SEXP MY_conncomp(SEXP flist) PROTECT(result = allocVector(INTSXP, N)); resgroup = INTEGER(result); group = factors[0]->group; - for (i = 0; i < N; i++) - { + for (i = 0; i < N; i++) { resgroup[i] = vertices[0][group[i] - 1]; } @@ -484,22 +423,18 @@ SEXP MY_conncomp(SEXP flist) gpsiz = Calloc(comps, double); idx = Calloc(comps, int); - for (i = 0; i < comps; i++) - idx[i] = i; - for (i = 0; i < N; i++) - { + for (i = 0; i < comps; i++) idx[i] = i; + for (i = 0; i < N; i++) { gpsiz[resgroup[i] - 1]++; } revsort(gpsiz, idx, comps); Free(gpsiz); levtrl = Calloc(comps, int); - for (i = 0; i < comps; i++) - levtrl[idx[i]] = i + 1; + for (i = 0; i < comps; i++) levtrl[idx[i]] = i + 1; Free(idx); - for (i = 0; i < N; i++) - { + for (i = 0; i < N; i++) { resgroup[i] = levtrl[resgroup[i] - 1]; } Free(levtrl); diff --git a/src/kaczmarz.c b/src/kaczmarz.c index 1268ce9..0f47801 100644 --- a/src/kaczmarz.c +++ b/src/kaczmarz.c @@ -2,8 +2,7 @@ /* Need snprintf */ #include -typedef struct -{ +typedef struct { int nowdoing; double **source; @@ -31,18 +30,18 @@ typedef struct #endif } KARG; -static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double *x, - double eps, mysize_t *work, int *stop, LOCK_T lock) -{ - +static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, + double *x, double eps, mysize_t *work, int *stop, + LOCK_T lock) { /* The factors define a matrix D, we will solve Dx = R There are N rows in D, each row a_i contains e non-zero entries, one for each factor, the level at that position. We iterate on k, start with x=0, an iteration consists - of adding a multiple of row i (with i = k %% N), the multiplying coefficient - is (R[i] - (a_i,x))/e + of adding a multiple of row i (with i = k %% N), the multiplying + coefficient is (R[i] - (a_i,x))/e - To get better memory locality, we create an (e X N)-matrix with the non-zero indices + To get better memory locality, we create an (e X N)-matrix with the + non-zero indices Now, as an afterthought we have made possible to interact a covariate with each factor, it's in factors[i]->x @@ -69,44 +68,38 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * int hasinteract = 0; for (int i = 0; i < e; i++) - if (NULL != factors[i]->x) - hasinteract = 1; + if (NULL != factors[i]->x) hasinteract = 1; mysize_t workpos = 0; /* Do the doubles first to keep alignment */ double *newR = (double *)work; - mysize_t *indices = (mysize_t *)&work[workpos += N * sizeof(double) / sizeof(mysize_t)]; + mysize_t *indices = + (mysize_t *)&work[workpos += N * sizeof(double) / sizeof(mysize_t)]; mysize_t *perm = (mysize_t *)&work[workpos += e * N]; mysize_t *prev = (mysize_t *)&work[workpos += N]; mysize_t *this = (mysize_t *)&work[workpos += e]; if (!hasinteract) - for (int i = 0; i < e; i++) - prev[i] = 0; + for (int i = 0; i < e; i++) prev[i] = 0; newN = 0; ie = 0; - for (mysize_t i = 0; i < N; i++) - { + for (mysize_t i = 0; i < N; i++) { perm[i] = i; - for (int j = 0; j < e; j++) - { + for (int j = 0; j < e; j++) { this[j] = factors[j]->group[i]; } - if (hasinteract || (memcmp(this, prev, e * sizeof(int)) != 0)) - { + if (hasinteract || (memcmp(this, prev, e * sizeof(int)) != 0)) { int nlev = 0; /* not duplicate, store in indices */ - for (int j = 0; j < e; j++) - { + for (int j = 0; j < e; j++) { indices[ie + j] = this[j] - 1 + nlev; nlev += factors[j]->nlevels; } newR[newN] = R[i]; newN++; ie += e; - if (!hasinteract) - memcpy(prev, this, e * sizeof(int)); + if (!hasinteract) memcpy(prev, this, e * sizeof(int)); } } @@ -124,15 +117,12 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * of the routine, only the speed. */ LOCK(lock); - if (e > 1) - { - for (mysize_t i = newN - 1; i > 0; i--) - { + if (e > 1) { + for (mysize_t i = newN - 1; i > 0; i--) { mysize_t j; /* Pick j between 0 and i inclusive */ j = (mysize_t)floor((i + 1) * unif_rand()); - if (j == i) - continue; + if (j == i) continue; /* exchange newR[i] and newR[j] as well as indices[i*e:i*e+e-1] and indices[j*e:j*e+e-1] */ @@ -142,8 +132,7 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * mysize_t itmp = perm[j]; perm[j] = perm[i]; perm[i] = itmp; - for (mysize_t k = 0; k < e; k++) - { + for (mysize_t k = 0; k < e; k++) { mysize_t itmp; itmp = indices[j * e + k]; indices[j * e + k] = indices[i * e + k]; @@ -155,26 +144,23 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * /* Then, do the Kaczmarz iterations */ norm2 = 0.0; - for (mysize_t i = 0; i < newN; i++) - norm2 += newR[i] * newR[i]; + for (mysize_t i = 0; i < newN; i++) norm2 += newR[i] * newR[i]; norm2 = sqrt(norm2); prevdiff = 2 * norm2; neweps = eps * norm2; - do - { - mysize_t ie = 0; /* equals i*e; integer multiplication is slow, keep track instead */ + do { + mysize_t ie = + 0; /* equals i*e; integer multiplication is slow, keep track instead */ double diff = 0.0; - for (mysize_t i = 0; i < newN; i++, ie += e) - { + for (mysize_t i = 0; i < newN; i++, ie += e) { const mysize_t ip = perm[i]; double upd = 0.0; double ai2 = 0.0; upd = newR[i]; /* Subtract inner product */ - for (int j = 0; j < e; j++) - { + for (int j = 0; j < e; j++) { const mysize_t idx = indices[ie + j]; const double *fx = factors[j]->x; const double w = (fx == NULL) ? 1.0 : fx[ip]; @@ -183,8 +169,7 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * } /* Update */ upd /= ai2; - for (int j = 0; j < e; j++) - { + for (int j = 0; j < e; j++) { const mysize_t idx = indices[ie + j]; const double *fx = factors[j]->x; const double w = (fx == NULL) ? upd : upd * fx[ip]; @@ -197,18 +182,19 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * sd = sqrt(diff); c = sd / prevdiff; prevdiff = sd; - if (c >= 1.0 && iter > 20) - { + if (c >= 1.0 && iter > 20) { char buf[256]; - snprintf(buf, sizeof(buf), "Kaczmarz failed in iter %d*%d, c=1-%.0e, delta=%.1e, eps=%.1e\n", iter, newN, 1.0 - c, sd, neweps); + snprintf( + buf, sizeof(buf), + "Kaczmarz failed in iter %d*%d, c=1-%.0e, delta=%.1e, eps=%.1e\n", + iter, newN, 1.0 - c, sd, neweps); pushmsg(buf, lock); break; } #ifdef NOTHREADS R_CheckUserInterrupt(); #else - if (*stop != 0) - return (0); + if (*stop != 0) return (0); #endif } while (sd >= neweps * (1.0 - c) && neweps > 1e-15); @@ -216,11 +202,9 @@ static double kaczmarz(FACTOR *factors[], int e, mysize_t N, double *R, double * } #ifdef WIN -DWORD WINAPI kaczmarz_thr(LPVOID varg) -{ +DWORD WINAPI kaczmarz_thr(LPVOID varg) { #else -static void *kaczmarz_thr(void *varg) -{ +static void *kaczmarz_thr(void *varg) { #endif KARG *arg = (KARG *)varg; int myid; @@ -229,14 +213,12 @@ static void *kaczmarz_thr(void *varg) LOCK(arg->lock); myid = arg->threadnum++; UNLOCK(arg->lock); - while (1) - { + while (1) { LOCK(arg->lock); vecnum = arg->nowdoing++; UNLOCK(arg->lock); - if (vecnum >= arg->numvec) - break; + if (vecnum >= arg->numvec) break; #ifdef HAVE_THREADNAME char thrname[16]; char buf[256]; @@ -246,17 +228,16 @@ static void *kaczmarz_thr(void *varg) STNAME(thrname); #endif - (void)kaczmarz(arg->factors, arg->e, arg->N, - arg->source[vecnum], arg->target[vecnum], - arg->eps, arg->work[myid], &arg->stop, arg->lock); + (void)kaczmarz(arg->factors, arg->e, arg->N, arg->source[vecnum], + arg->target[vecnum], arg->eps, arg->work[myid], &arg->stop, + arg->lock); } #ifndef NOTHREADS #ifndef WIN LOCK(arg->lock); arg->running--; #ifdef HAVE_SEM - if (arg->running == 0) - sem_post(&arg->finished); + if (arg->running == 0) sem_post(&arg->finished); #endif UNLOCK(arg->lock); #endif @@ -264,9 +245,7 @@ static void *kaczmarz_thr(void *varg) return 0; } -SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) -{ - +SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) { double eps = REAL(Reps)[0]; /* double *R = REAL(RR);*/ FACTOR **factors; @@ -300,8 +279,7 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) /* Set up the message stack */ initmsg(); #endif - if (!isNull(initial)) - { + if (!isNull(initial)) { init = REAL(initial); } PROTECT(flist = AS_LIST(flist)); @@ -310,14 +288,13 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) factors = makefactors(flist, 0, NULL); numfac = 0; - for (FACTOR **f = factors; *f != NULL; f++) - numfac++; + for (FACTOR **f = factors; *f != NULL; f++) numfac++; N = LENGTH(VECTOR_ELT(flist, 0)); - for (int i = 0; i < numfac; i++) - sumlev += factors[i]->nlevels; + for (int i = 0; i < numfac; i++) sumlev += factors[i]->nlevels; if (!isNull(initial) && LENGTH(initial) != sumlev) - error("Initial vector must have length %d, but is %d\n", sumlev, LENGTH(initial)); + error("Initial vector must have length %d, but is %d\n", sumlev, + LENGTH(initial)); /* Then the vectors */ @@ -329,20 +306,18 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) /* First, count the number of vectors in total */ numvec = 0; - for (int i = 0; i < listlen; i++) - { + for (int i = 0; i < listlen; i++) { SEXP elt = VECTOR_ELT(vlist, i); /* Each entry in the list is either a vector or a matrix */ - if (!isMatrix(elt)) - { + if (!isMatrix(elt)) { if (LENGTH(elt) != N) - error("Vector length (%d) must be equal to factor length (%d)", LENGTH(elt), N); + error("Vector length (%d) must be equal to factor length (%d)", + LENGTH(elt), N); numvec++; - } - else - { + } else { if (nrows(elt) != N) - error("Vector length must be equal to factor length %d %d ", nrows(elt), N); + error("Vector length must be equal to factor length %d %d ", nrows(elt), + N); numvec += ncols(elt); } } @@ -353,16 +328,13 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) target = (double **)R_alloc(numvec, sizeof(double *)); /* Loop through list again to set up result structure */ cnt = 0; - for (int i = 0; i < listlen; i++) - { + for (int i = 0; i < listlen; i++) { SEXP elt = VECTOR_ELT(vlist, i); - if (!isReal(elt)) - { + if (!isReal(elt)) { elt = PROTECT(coerceVector(elt, REALSXP)); protectcount++; } - if (!isMatrix(elt)) - { + if (!isMatrix(elt)) { /* It's a vector */ SEXP resvec; vectors[cnt] = REAL(elt); @@ -371,9 +343,7 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) SET_VECTOR_ELT(reslist, i, resvec); UNPROTECT(1); cnt++; - } - else - { + } else { /* It's a matrix */ int cols = ncols(elt); SEXP mtx; @@ -382,8 +352,7 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) SET_VECTOR_ELT(reslist, i, mtx); UNPROTECT(1); /* Set up pointers */ - for (int j = 0; j < cols; j++) - { + for (int j = 0; j < cols; j++) { vectors[cnt] = REAL(elt) + j * N; target[cnt] = REAL(mtx) + j * sumlev; cnt++; @@ -391,23 +360,18 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) } } - for (int cnt = 0; cnt < numvec; cnt++) - { + for (int cnt = 0; cnt < numvec; cnt++) { if (init != 0) - for (int i = 0; i < sumlev; i++) - target[cnt][i] = init[i]; + for (int i = 0; i < sumlev; i++) target[cnt][i] = init[i]; else - for (int i = 0; i < sumlev; i++) - target[cnt][i] = 0.0; + for (int i = 0; i < sumlev; i++) target[cnt][i] = 0.0; } /* set up for threading */ numthr = cores; - if (numthr > numvec) - numthr = numvec; - if (numthr < 1) - numthr = 1; + if (numthr > numvec) numthr = numvec; + if (numthr < 1) numthr = 1; GetRNGstate(); #ifndef NOTHREADS #ifdef WIN @@ -438,53 +402,53 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) arg.stop = 0; arg.work = (mysize_t **)R_alloc(numthr, sizeof(mysize_t *)); - // when allocating the work, we use mysize_t, but parts of it is accessed as double which may be - // larger. So we allocate some more (8*sizeof(mysize_t) more), and adjust the address so it's aligned on a double - // When using it, make sure we do all the doubles first. + // when allocating the work, we use mysize_t, but parts of it is accessed as + // double which may be larger. So we allocate some more (8*sizeof(mysize_t) + // more), and adjust the address so it's aligned on a double When using it, + // make sure we do all the doubles first. #ifdef NOTHREADS - arg.work[0] = (mysize_t *)R_alloc(numfac * N + N * sizeof(double) / sizeof(mysize_t) + N + 2 * numfac + 8, sizeof(mysize_t)); + arg.work[0] = (mysize_t *)R_alloc( + numfac * N + N * sizeof(double) / sizeof(mysize_t) + N + 2 * numfac + 8, + sizeof(mysize_t)); uintptr_t amiss = (uintptr_t)arg.work[0] % sizeof(double); if (amiss != 0) arg.work[0] = (mysize_t *)((uintptr_t)arg.work[0] + sizeof(double) - amiss); kaczmarz_thr((void *)&arg); #else /* Do it in separate threads */ - for (thr = 0; thr < numthr; thr++) - { + for (thr = 0; thr < numthr; thr++) { // allocate some thread-specific storage, we can't use R_alloc in a thread - arg.work[thr] = (mysize_t *)R_alloc(numfac * N + N * sizeof(double) / sizeof(mysize_t) + N + 2 * numfac + 8, sizeof(mysize_t)); + arg.work[thr] = (mysize_t *)R_alloc( + numfac * N + N * sizeof(double) / sizeof(mysize_t) + N + 2 * numfac + 8, + sizeof(mysize_t)); uintptr_t amiss = (uintptr_t)arg.work[thr] % sizeof(double); if (amiss != 0) - arg.work[thr] = (mysize_t *)((uintptr_t)arg.work[thr] + sizeof(double) - amiss); + arg.work[thr] = + (mysize_t *)((uintptr_t)arg.work[thr] + sizeof(double) - amiss); #ifdef WIN - threads[thr] = CreateThread(NULL, 0, kaczmarz_thr, &arg, 0, &threadids[thr]); - if (0 == threads[thr]) - error("Failed to create kaczmarz thread"); + threads[thr] = + CreateThread(NULL, 0, kaczmarz_thr, &arg, 0, &threadids[thr]); + if (0 == threads[thr]) error("Failed to create kaczmarz thread"); #else int stat = pthread_create(&threads[thr], NULL, kaczmarz_thr, &arg); - if (0 != stat) - error("Failed to create kaczmarz thread, stat=%d", stat); + if (0 != stat) error("Failed to create kaczmarz thread, stat=%d", stat); #endif } /* wait for completion */ /* We want to check for interrupts regularly, and set a stop flag */ - while (1) - { + while (1) { printmsg(arg.lock); - if (arg.stop == 0 && checkInterrupt()) - { + if (arg.stop == 0 && checkInterrupt()) { REprintf("...stopping Kaczmarz threads...\n"); arg.stop = 1; } #ifdef WIN - if (WaitForMultipleObjects(numthr, threads, TRUE, 3000) != WAIT_TIMEOUT) - { - for (thr = 0; thr < numthr; thr++) - { + if (WaitForMultipleObjects(numthr, threads, TRUE, 3000) != WAIT_TIMEOUT) { + for (thr = 0; thr < numthr; thr++) { CloseHandle(threads[thr]); } /* Print any remaining messages */ @@ -498,17 +462,13 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) struct timespec atmo = {0, 50000000}; /* Kludge in MacOSX because no timedwait */ - if (arg.stop == 0) - nanosleep(&atmo, NULL); - if (arg.stop == 1 || arg.running == 0) - { + if (arg.stop == 0) nanosleep(&atmo, NULL); + if (arg.stop == 1 || arg.running == 0) { #else struct timespec tmo = {time(NULL) + 3, 0}; - if (arg.stop == 1 || sem_timedwait(&arg.finished, &tmo) == 0) - { + if (arg.stop == 1 || sem_timedwait(&arg.finished, &tmo) == 0) { #endif - for (thr = 0; thr < numthr; thr++) - { + for (thr = 0; thr < numthr; thr++) { (void)pthread_join(threads[thr], NULL); } #ifdef HAVE_SEM @@ -524,8 +484,7 @@ SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores) } #endif PutRNGstate(); - if (arg.stop == 1) - error("Kaczmarz interrupted"); + if (arg.stop == 1) error("Kaczmarz interrupted"); UNPROTECT(protectcount); return (reslist); } diff --git a/src/lfe.c b/src/lfe.c index 3fdd8f0..01374eb 100644 --- a/src/lfe.c +++ b/src/lfe.c @@ -2,13 +2,12 @@ SEXP df_string; int LFE_GLOBAL_THREADS = 1; -SEXP MY_threads(SEXP rt) -{ - if (LENGTH(rt) < 1) - return R_NilValue; +SEXP MY_threads(SEXP rt) { + if (LENGTH(rt) < 1) return R_NilValue; LFE_GLOBAL_THREADS = INTEGER(rt)[0]; return R_NilValue; } + static R_CallMethodDef callMethods[] = { {"conncomp", (DL_FUNC)&MY_conncomp, 1}, {"wwcomp", (DL_FUNC)&MY_wwcomp, 1}, @@ -28,8 +27,7 @@ static R_CallMethodDef callMethods[] = { // {"threads", (DL_FUNC) &MY_threads, 1}, {NULL, NULL, 0}}; -void attribute_visible R_init_lfe(DllInfo *info) -{ +void attribute_visible R_init_lfe(DllInfo *info) { /* register our routines */ (void)R_registerRoutines(info, NULL, callMethods, NULL, NULL); (void)R_useDynamicSymbols(info, FALSE); @@ -38,10 +36,8 @@ void attribute_visible R_init_lfe(DllInfo *info) LFE_GLOBAL_THREADS = 1; } -void attribute_visible R_unload_lfe(DllInfo *info) -{ - if (info != NULL) - { - }; // avoid pedantic warning about unused parameter +void attribute_visible R_unload_lfe(DllInfo *info) { + if (info != NULL) { + }; // avoid pedantic warning about unused parameter (void)R_ReleaseObject(df_string); } diff --git a/src/lfe.h b/src/lfe.h index 532bdaf..d166dcd 100644 --- a/src/lfe.h +++ b/src/lfe.h @@ -16,8 +16,8 @@ #define _GNU_SOURCE /* to find pthread_setname_np */ #endif -#include #include +#include #endif #endif @@ -32,7 +32,8 @@ int pthread_setname_np(pthread_t thread, const char *name); int pthread_setname_np(const char *); #define STNAME(s) pthread_setname_np(s) #elif __FreeBSD__ -// FreeBSD & OpenBSD: function name is slightly different, and has no return value +// FreeBSD & OpenBSD: function name is slightly different, and has no return +// value void pthread_set_name_np(pthread_t tid, const char *name); #define STNAME(s) pthread_set_name_np(pthread_self(), s) #else @@ -47,19 +48,19 @@ void pthread_set_name_np(pthread_t tid, const char *name); #define STNAME(s) #endif -#include -#include -#include -#include -#include -#include -#include #include -#include -#include +#include #include #include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #if defined(R_VERSION) && R_VERSION >= R_Version(3, 0, 0) typedef R_xlen_t mybigint_t; @@ -96,8 +97,7 @@ typedef int mysize_t; #endif /* My internal definition of a factor */ -typedef struct -{ +typedef struct { /* group[i] is the level of observation i */ int *group; /* invgpsize[j] is the 1/(size of level j) */ @@ -121,9 +121,9 @@ void printmsg(LOCK_T lock); SEXP MY_kaczmarz(SEXP flist, SEXP vlist, SEXP Reps, SEXP initial, SEXP Rcores); SEXP MY_wwcomp(SEXP flist); SEXP MY_conncomp(SEXP flist); -SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, - SEXP scores, SEXP quiet, SEXP gkacc, SEXP Rmeans, - SEXP weights, SEXP Rscale, SEXP attrs); +SEXP MY_demeanlist(SEXP vlist, SEXP flist, SEXP Ricpt, SEXP Reps, SEXP scores, + SEXP quiet, SEXP gkacc, SEXP Rmeans, SEXP weights, + SEXP Rscale, SEXP attrs); SEXP MY_scalecols(SEXP mat, SEXP vec); SEXP MY_pdaxpy(SEXP inX, SEXP inY, SEXP inbeta); SEXP MY_piproduct(SEXP inX, SEXP inY); diff --git a/src/utils.c b/src/utils.c index 7b2fe49..a6f3a45 100644 --- a/src/utils.c +++ b/src/utils.c @@ -5,40 +5,31 @@ #define FCONE #endif -SEXP MY_scalecols(SEXP mat, SEXP vec) -{ - if (!isMatrix(mat)) - error("first argument should be a matrix"); +SEXP MY_scalecols(SEXP mat, SEXP vec) { + if (!isMatrix(mat)) error("first argument should be a matrix"); mybigint_t col = ncols(mat), row = nrows(mat); - if (isMatrix(vec)) - { + if (isMatrix(vec)) { if (row != nrows(vec)) error("Rows of matrix should be the same as rows of vector"); double *cmat = REAL(mat); double *cvec = REAL(AS_NUMERIC(vec)); - for (mybigint_t j = 0; j < col; j++) - { + for (mybigint_t j = 0; j < col; j++) { double *cc = &cmat[j * row]; - for (mybigint_t i = 0; i < row; i++) - { + for (mybigint_t i = 0; i < row; i++) { double tmp = 0.0; - for (int k = 0; k < ncols(vec); k++) - tmp += cc[i] * cvec[i + k * row]; + for (int k = 0; k < ncols(vec); k++) tmp += cc[i] * cvec[i + k * row]; cc[i] = tmp; } } - } - else - { + } else { if (row != LENGTH(vec)) - error("length of vector %d is different from number of rows %d", LENGTH(vec), row); + error("length of vector %d is different from number of rows %ld", + LENGTH(vec), row); double *cmat = REAL(mat); double *cvec = REAL(AS_NUMERIC(vec)); - for (mybigint_t j = 0; j < col; j++) - { + for (mybigint_t j = 0; j < col; j++) { double *cc = &cmat[j * row]; - for (mybigint_t i = 0; i < row; i++) - cc[i] *= cvec[i]; + for (mybigint_t i = 0; i < row; i++) cc[i] *= cvec[i]; } } return mat; @@ -53,8 +44,7 @@ SEXP MY_scalecols(SEXP mat, SEXP vec) X + t(beta * t(Y)), so this is a utility to save some costly transposes. used in cgsolve() */ -SEXP MY_pdaxpy(SEXP inX, SEXP inY, SEXP inbeta) -{ +SEXP MY_pdaxpy(SEXP inX, SEXP inY, SEXP inbeta) { mybigint_t col = ncols(inX), row = nrows(inX); if (col != ncols(inY) || row != nrows(inY)) error("X and Y should have the same shape"); @@ -66,14 +56,12 @@ SEXP MY_pdaxpy(SEXP inX, SEXP inY, SEXP inbeta) SEXP res; PROTECT(res = allocMatrix(REALSXP, row, col)); double *pres = REAL(res); - for (mybigint_t j = 0; j < col; j++) - { + for (mybigint_t j = 0; j < col; j++) { double b = beta[j]; double *out = &pres[j * row]; double *xin = &X[j * row]; double *yin = &Y[j * row]; - for (mybigint_t i = 0; i < row; i++) - { + for (mybigint_t i = 0; i < row; i++) { out[i] = xin[i] + b * yin[i]; } } @@ -85,8 +73,7 @@ SEXP MY_pdaxpy(SEXP inX, SEXP inY, SEXP inbeta) compute inner product pairwise of the columns of matrices X and Y This is diag(crossprod(X,Y)). */ -SEXP MY_piproduct(SEXP inX, SEXP inY) -{ +SEXP MY_piproduct(SEXP inX, SEXP inY) { mybigint_t col = ncols(inX), row = nrows(inX); if (col != ncols(inY) || row != nrows(inY)) error("X and Y should have the same shape"); @@ -95,13 +82,11 @@ SEXP MY_piproduct(SEXP inX, SEXP inY) SEXP res; PROTECT(res = allocVector(REALSXP, col)); double *pres = REAL(res); - for (mybigint_t j = 0; j < col; j++) - { + for (mybigint_t j = 0; j < col; j++) { double *xin = &X[j * row]; double *yin = &Y[j * row]; pres[j] = 0.0; - for (mybigint_t i = 0; i < row; i++) - { + for (mybigint_t i = 0; i < row; i++) { pres[j] += xin[i] * yin[i]; } } @@ -110,27 +95,20 @@ SEXP MY_piproduct(SEXP inX, SEXP inY) } // copy-free dimnames<- -SEXP MY_setdimnames(SEXP obj, SEXP nm) -{ - if (!isNull(obj)) - setAttrib(obj, R_DimNamesSymbol, nm); +SEXP MY_setdimnames(SEXP obj, SEXP nm) { + if (!isNull(obj)) setAttrib(obj, R_DimNamesSymbol, nm); return (R_NilValue); } /* Compute and return alpha * bread %*% meat %*% bread */ -SEXP MY_sandwich(SEXP inalpha, SEXP inbread, SEXP inmeat) -{ +SEXP MY_sandwich(SEXP inalpha, SEXP inbread, SEXP inmeat) { double alpha = REAL(AS_NUMERIC(inalpha))[0]; - if (!isMatrix(inbread)) - error("bread must be a matrix"); - if (!isMatrix(inmeat)) - error("bread must be a matrix"); - if (ncols(inbread) != nrows(inbread)) - error("bread must be square matrix"); + if (!isMatrix(inbread)) error("bread must be a matrix"); + if (!isMatrix(inmeat)) error("bread must be a matrix"); + if (ncols(inbread) != nrows(inbread)) error("bread must be square matrix"); - if (ncols(inmeat) != nrows(inmeat)) - error("meat must be square matrix"); + if (ncols(inmeat) != nrows(inmeat)) error("meat must be square matrix"); if (ncols(inmeat) != ncols(inbread)) error("bread and meat must have the same size"); @@ -156,34 +134,28 @@ SEXP MY_sandwich(SEXP inalpha, SEXP inbread, SEXP inmeat) } // copy-free dsyrk C = beta*C + alpha * A' A -SEXP MY_dsyrk(SEXP inbeta, SEXP inC, SEXP inalpha, SEXP inA) -{ +SEXP MY_dsyrk(SEXP inbeta, SEXP inC, SEXP inalpha, SEXP inA) { double beta = REAL(AS_NUMERIC(inbeta))[0]; double alpha = REAL(AS_NUMERIC(inalpha))[0]; - if (!isMatrix(inC)) - error("C must be a matrix"); - if (!isMatrix(inA)) - error("A must be a matrix"); + if (!isMatrix(inC)) error("C must be a matrix"); + if (!isMatrix(inA)) error("A must be a matrix"); - if (ncols(inC) != nrows(inC)) - { + if (ncols(inC) != nrows(inC)) { error("C must be a square matrix, it is %d x %d", nrows(inC), ncols(inC)); } int N = nrows(inC); double *C = REAL(inC); - if (ncols(inA) != ncols(inC)) - { - error("A (%d x %d) must have the same number of columns as C (%d x %d)", nrows(inA), ncols(inA), nrows(inC), nrows(inC)); + if (ncols(inA) != ncols(inC)) { + error("A (%d x %d) must have the same number of columns as C (%d x %d)", + nrows(inA), ncols(inA), nrows(inC), nrows(inC)); } int K = nrows(inA); double *A = REAL(inA); F77_CALL(dsyrk) ("U", "T", &N, &K, &alpha, A, &K, &beta, C, &N FCONE FCONE); // fill in the lower triangular part - for (mybigint_t row = 0; row < N; row++) - { - for (mybigint_t col = 0; col < row; col++) - { + for (mybigint_t row = 0; row < N; row++) { + for (mybigint_t col = 0; col < row; col++) { C[col * N + row] = C[row * N + col]; } } @@ -191,17 +163,14 @@ SEXP MY_dsyrk(SEXP inbeta, SEXP inC, SEXP inalpha, SEXP inA) } // debugging memory copy -SEXP MY_address(SEXP x) -{ +SEXP MY_address(SEXP x) { char chr[30]; snprintf(chr, sizeof(chr), "adr=%p, named=%d", (void *)x, NAMED(x)); return (mkString(chr)); } -SEXP MY_named(SEXP x, SEXP n) -{ - if (isNull(n)) - { +SEXP MY_named(SEXP x, SEXP n) { + if (isNull(n)) { // return NAMED status SEXP res = allocVector(INTSXP, 1); PROTECT(res); @@ -217,19 +186,16 @@ SEXP MY_named(SEXP x, SEXP n) } /* Trickery to check for interrupts when using threads */ -static void chkIntFn(void *dummy) -{ - if (dummy == NULL) - { - }; // avoid pedantic warning about unused parameter +static void chkIntFn(void *dummy) { + if (dummy == NULL) { + }; // avoid pedantic warning about unused parameter R_CheckUserInterrupt(); } /* this will call the above in a top-level context so it won't longjmp-out of context */ -extern int checkInterrupt(void) -{ +extern int checkInterrupt(void) { return (R_ToplevelExec(chkIntFn, NULL) == 0); } @@ -241,46 +207,36 @@ extern int checkInterrupt(void) static char *msgstack[MSGLIM]; static int msgptr; -extern void initmsg(void) -{ - msgptr = 0; -} +extern void initmsg(void) { msgptr = 0; } /* Craft our own strdup, it's not supported everywhere */ -static char *mystrdup(char *s) -{ +static char *mystrdup(char *s) { char *sc = (char *)malloc(strlen(s) + 1); - if (sc != NULL) - strcpy(sc, s); + if (sc != NULL) strcpy(sc, s); return (sc); } -void pushmsg(char *s, LOCK_T lock) -{ +void pushmsg(char *s, LOCK_T lock) { #ifdef NOTHREADS REprintf(s); #else LOCK(lock); - if (msgptr < MSGLIM) - { + if (msgptr < MSGLIM) { msgstack[msgptr++] = mystrdup(s); } UNLOCK(lock); #endif } -void printmsg(LOCK_T lock) -{ +void printmsg(LOCK_T lock) { #ifdef NOTHREADS return; #else char *s; int i; LOCK(lock); - for (i = 0; i < msgptr; i++) - { + for (i = 0; i < msgptr; i++) { s = msgstack[i]; - if (s != NULL) - { - REprintf(s); + if (s != NULL) { + REprintf("%s", s); free(s); } }