From d128914e1f65e642054adca67ccd27869b68021e Mon Sep 17 00:00:00 2001 From: AndrewC1998 Date: Wed, 5 Jun 2019 10:16:22 +0100 Subject: [PATCH 1/4] 05/06/2019 Updated to remove the cpt.reg functionality and then port through to changepoint instead. No further changes should need to be made as it will all be done in changepoint from here. Users should not notice any change. AC --- R/CptReg.R | 199 --------------------- R/envcpt.R | 12 +- src/C_cptreg.c | 474 ------------------------------------------------- 3 files changed, 6 insertions(+), 679 deletions(-) delete mode 100755 R/CptReg.R delete mode 100755 src/C_cptreg.c diff --git a/R/CptReg.R b/R/CptReg.R deleted file mode 100755 index 862e2aa..0000000 --- a/R/CptReg.R +++ /dev/null @@ -1,199 +0,0 @@ -cpt.reg <- function(data, penalty="MBIC", pen.value=0, method="AMOC", dist="Normal", - class=TRUE, param.estimates=TRUE, minseglen=3, shape = 0, tol = 1e-07){ - MBIC=0 - ##Check arguments are valid - if(!is.array(data) || !is.numeric(data)) ##Further checks applied later - stop("Argument 'data' must be a numerical matrix/array.") - if(!is.character(penalty) || length(penalty)>1) - stop("Argument 'penelty' is invalid.") - #value of 'penalty' & 'pen.value' checked within changepoint::penalty_decision - if(!is.character(method) || length(method)>1) - stop("Argument 'method' is invalid.") - if(method!="AMOC" && method != "PELT") ##RESTRICTION IN USE - stop("Invalid method, must be AMOC or PELT.") - if(!is.character(dist) || length(dist)>1) - stop("Argument 'dist' is invalid.") - if(dist != "Normal"){ ##RESTRICTION IN USE - warning(paste0("dist = ",dist," is not supported. Converted to dist='Normal'")) - dist <- "Normal" - } - if(!is.logical(class) || length(class)>1) - stop("Argument 'class' is invalid.") - if(!is.logical(param.estimates) || length(param.estimates)>1) - stop("Argument 'param.estimates' is invalid.") - if(!is.numeric(minseglen) || length(minseglen)>1) - stop("Argument 'minseglen' is invalid.") - if(minseglen <= 0 || minseglen%%1 != 0) ##Further checks applied later - stop("Argument 'minseglen' must be positive integer.") - if(!is.numeric(tol) || length(tol)!=1) - stop("Argument 'tol' is invalid.") - if(tol<0) stop("Argument 'tol' must be positive.") - ##Argument shape is assessed by the command where it is to be used. - - #Single data set? convert to multiple data set array. - if(length(dim(data)) == 2) data <- array(data,dim=c(1,dim(data))) - - ans <- vector("list",dim(data)[1]) - for(i in 1:dim(data)[1]){ ##To each data set. - #Check format for ith data set - datai <- check_data(data[i,,]) - - #Evaluate penalty value - pen.value <- changepoint::penalty_decision(penalty=penalty, - pen.value=pen.value, n=nrow(datai), diffparam=ncol(datai), - asymcheck="cpt.reg", method=method) - if(penalty=="MBIC"){MBIC=1} - - #Check value of minseglen - if(minseglen < (ncol(datai)-1)){ - warning(paste("minseglen is too small, set to:",ncol(datai))) - minsegleni <- ncol(datai) - }else if(nrow(datai) < (2*minseglen)){ - stop("Minimum segment length is too large to include a change in this data.") - }else{ - minsegleni <- minseglen - } - - CPTS <- ChangepointRegression(data=datai, penalty=penalty, - penalty.value=pen.value, method=method, dist=dist, minseglen=minseglen, - shape = shape, tol=tol, cpts.only=class, MBIC=MBIC) - - if(class){ - #Convert to cpt.reg object - ansi <- new("cpt.reg") - data.set(ansi) <- datai - cpttype(ansi) <- "regression" - method(ansi) <- method - distribution(ansi) <- dist - pen.type(ansi) <- penalty - pen.value(ansi) <- pen.value - cpts(ansi) <- CPTS - if(method=="PELT") ncpts.max(ansi) <- Inf - if(param.estimates) ansi = param(ansi) - ans[[i]] <- ansi - }else{ - #store what is returned from calculation - ans[[i]] <- CPTS - } - } - - ##Return result, only first case if there is only a singe data set - if(dim(data)[1]==1){ return(ans[[1]]) }else{ return(ans) } -} - -#################### - -check_data <- function(data, minseglen=3){ - if(!is.array(data) || !is.numeric(data)) - stop("Argument 'data' must be a numerical matrix.") - if(length(dim(data))!=2) - stop("Argument 'data' must be a numerical matrix.") - if(!is.numeric(minseglen) || length(minseglen)>1) - stop("Argument 'minseglen' is invalid.") - - n <- nrow(data) - p <- ncol(data)-1 - if(p==0) stop("Dimension of data is 1, no regressors found.") - if(n1){ - i <- which(intercept)[-1]+1 - warning("Multiple intercepts found. Keeping only first instance.") - data <- data[,-i] - } - return(data) -} - - -ChangepointRegression <- function(data, penalty="MBIC", penalty.value=0, - method="AMOC", dist="Normal", minseglen=3, cpts.only=TRUE, shape=0, MBIC=0, tol=1e-07){ - ##Assume all input arguments are entered correctly - if(!is.logical(cpts.only) && length(cpts.only)>1) - stop("Argument 'cpts.only' is invalid.") - if(method=="AMOC" && dist=="Normal"){ - out <- CptReg_AMOC_Normal(data, penalty, penalty.value, minseglen, shape, MBIC, tol) - }else if(method=="PELT" && dist=="Normal"){ - out <- CptReg_PELT_Normal(data, penalty.value, minseglen, shape,MBIC, tol) - }else{ - stop("Changepoint in regression method not recognised.") - } - if(cpts.only){ - return(sort(out$cpts)) - }else{ - return(out) - } -} - -CptReg_AMOC_Normal <- function(data, penalty="MBIC", penalty.value=0, minseglen=3, - shape=0, MBIC=0, tol=1e-07){ - n <- as.integer(nrow(data)) - p <- as.integer(ncol(data)-1) - if(p<1 || n0=-2logLik with this fixed variance - if(!is.numeric(shape) || length(shape)!=1) - stop("Argument 'shape' is invalid.") - - ##Clean-up on exit - answer=list() - answer[[5]]=1 - on.exit(.C("Free_CptReg_Normal_AMOC",answer[[6]],PACKAGE="EnvCpt")) - - answer <- .C("CptReg_Normal_AMOC", data=as.double(data), n=as.integer(n), - m=as.integer(p+1), pen=as.double(penalty.value), err=0L, - shape=as.double(shape), minseglen=as.integer(minseglen), tol=as.double(tol), - tau=0L, nulllike=vector("double",1), taulike=vector("double",1), - tmplike=vector("double",n), MBIC=as.integer(MBIC),PACKAGE="EnvCpt") - - #Check if error has occured - if(answer$err!=0) stop("C code error:",answer$err,call.=F) - tmp <- c(answer$tau,answer$nulllike,answer$taulike) - #Following line in cpt.mean (RSS) but not in old cpt.reg (-2Loglik) - #if(penalty=="MBIC") tmp[3] = tmp[3] + log(tmp[1]) + log(n-tmp[1]+1) - out <- changepoint::decision(tau = tmp[1], null = tmp[2], alt = tmp[3], - penalty = penalty, n=n, diffparam=p, pen.value = penalty.value) - names(out) <- c("cpts","pen.value") - return(out) #return list of cpts & pen.value -} - -CptReg_PELT_Normal <- function(data, penalty.value=0, minseglen=3, shape=0, - MBIC=0, tol=1e-07){ - n <- as.integer(nrow(data)) - p <- as.integer(ncol(data)-1) - #tol <- 1e-07 #Rank tolerance (see lm.fit) - #shape <- -1 #-1=RSS,0=-2logLik,>0=-2logLik with this fixed variance - if(!is.numeric(shape) || length(shape)!=1) - stop("Argument 'shape' is invalid.") - - #Check if error has occured - answer=list() - answer[[6]]=1 - on.exit(.C("Free_CptReg_Normal_PELT",answer[[6]],PACKAGE="EnvCpt")) - - answer <- .C("CptReg_Normal_PELT", data=as.double(data), n=n, - m=as.integer(p+1), pen=as.double(penalty.value), cpt=vector("integer",n), - err=0L, shape=as.double(shape), minseglen=as.integer(minseglen), - tol=as.double(tol), lastchangelike=vector("double",n+1), - lastchangecpts=vector("integer",n+1), numchangecpts=vector("integer",n+1), - MBIC=as.integer(MBIC), PACKAGE="EnvCpt") - - if(answer$err!=0){ - stop("C code error:",answer$err,call.=F) - } - return(list(lastchangecpts=answer$lastchangecpts, - cpts=sort(answer$cpt[answer$cpt>0]), - lastchangelike=answer$lastchangelike, - ncpts=answer$numchangecpts)) -} - - - -######################################################################## - - - diff --git a/R/envcpt.R b/R/envcpt.R index 325e144..d17a99a 100755 --- a/R/envcpt.R +++ b/R/envcpt.R @@ -55,7 +55,7 @@ envcpt=function(data,models=c("mean","meancpt","meanar1","meanar2","meanar1cpt", # mean+AR1+cpt if(any(models=="meanar1cpt")|any(models==5)){ - meanarcpt.fit=cpt.reg(cbind(data[-1],rep(1,n-1),data[-n]),method="PELT",minseglen=minseglen,...) # default MBIC penalty + meanarcpt.fit=changepoint::cpt.reg(cbind(data[-1],rep(1,n-1),data[-n]),method="PELT",minseglen=minseglen,...) # default MBIC penalty if(ncpts(meanarcpt.fit)==0){ if(any(models=="meanar1")|any(models==3)){meanarcpt.loglik=meanar.loglik} else{meanarcpt.loglik=logLik(meanarcpt.fit)[1]} # function from changepoint @@ -72,7 +72,7 @@ envcpt=function(data,models=c("mean","meancpt","meanar1","meanar2","meanar1cpt", # mean+AR2+cpt if(any(models=="meanar2cpt")|any(models==6)){ - meanar2cpt.fit=cpt.reg(cbind(data[-c(1:2)],rep(1,n-2),data[2:(n-1)],data[1:(n-2)]),method="PELT",minseglen=minseglen,...) # default MBIC penalty + meanar2cpt.fit=changepoint::cpt.reg(cbind(data[-c(1:2)],rep(1,n-2),data[2:(n-1)],data[1:(n-2)]),method="PELT",minseglen=minseglen,...) # default MBIC penalty if(ncpts(meanar2cpt.fit)==0){ if(any(models=="meanar2")|any(models==4)){meanar2cpt.loglik=meanar.loglik} else{meanar2cpt.loglik=logLik(meanar2cpt.fit)[1]} # function from changepoint @@ -98,7 +98,7 @@ envcpt=function(data,models=c("mean","meancpt","meanar1","meanar2","meanar1cpt", # trend+cpt if(any(models=="trendcpt")|any(models==8)){ - trendcpt.fit=cpt.reg(cbind(data,rep(1,n),1:n),method='PELT',minseglen=minseglen,...) # default MBIC penalty + trendcpt.fit=changepoint::cpt.reg(cbind(data,rep(1,n),1:n),method='PELT',minseglen=minseglen,...) # default MBIC penalty if(ncpts(trendcpt.fit)==0){ if(any(models=="trend")|any(models==7)){trendcpt.loglik=trend.loglik} else{trendcpt.loglik=logLik(trendcpt.fit)[1]} # function from changepoint @@ -131,7 +131,7 @@ envcpt=function(data,models=c("mean","meancpt","meanar1","meanar2","meanar1cpt", # trend+AR1+cpt if(any(models=="trendar1cpt")|any(models==11)){ - trendarcpt.fit=cpt.reg(cbind(data[-1],rep(1,n-1),1:(n-1),data[-n]),method="PELT",minseglen=minseglen,...) # default MBIC penalty + trendarcpt.fit=changepoint::cpt.reg(cbind(data[-1],rep(1,n-1),1:(n-1),data[-n]),method="PELT",minseglen=minseglen,...) # default MBIC penalty if(ncpts(trendarcpt.fit)==0){ if(any(models=="trendar1")|any(models==9)){trendarcpt.loglik=trendar.loglik} else{trendarcpt.loglik=logLik(trendarcpt.fit)[1]} # function from changepoint @@ -148,7 +148,7 @@ envcpt=function(data,models=c("mean","meancpt","meanar1","meanar2","meanar1cpt", # trend+AR2+cpt if(any(models=="trendar2cpt")|any(models==12)){ - trendar2cpt.fit=cpt.reg(cbind(data[-c(1:2)],rep(1,n-2),1:(n-2),data[2:(n-1)],data[1:(n-2)]),method="PELT",minseglen=minseglen,...) # default MBIC penalty + trendar2cpt.fit=changepoint::cpt.reg(cbind(data[-c(1:2)],rep(1,n-2),1:(n-2),data[2:(n-1)],data[1:(n-2)]),method="PELT",minseglen=minseglen,...) # default MBIC penalty if(ncpts(trendar2cpt.fit)==0){ if(any(models=="trendar2")|any(models==10)){trendar2cpt.loglik=trendar2.loglik} else{trendar2cpt.loglik=logLik(trendar2cpt.fit)[1]} # function from changepoint @@ -186,4 +186,4 @@ envcpt=function(data,models=c("mean","meancpt","meanar1","meanar2","meanar1cpt", if(verbose==TRUE){close(pb)} class(out)="envcpt" return(out) -} \ No newline at end of file +} diff --git a/src/C_cptreg.c b/src/C_cptreg.c deleted file mode 100755 index 558a664..0000000 --- a/src/C_cptreg.c +++ /dev/null @@ -1,474 +0,0 @@ -#include -#include -#include // RK addition -#include // RK addition -#include // RK addition -#include // RK addition -#include // ST addition -#include // ST addition -#include -#include -#include -#include -#include - -static int *checklist; -static double *tmplike; -static double *Sumstats; - -//To register C functions to be read within R!!! -void R_init_EnvCpt(DllInfo *info){ - R_registerRoutines(info, NULL, NULL, NULL, NULL); - R_useDynamicSymbols(info, TRUE); -} - - -//Calculate summary statistics for the regression quadratic cost function. -void RegQuadCost_SS(double *X, int *n, int *nc, double *SS, int *m){ - //X - data as double vector of length n*p - //n - number or rows of X - //nc - numnber of columns of X - //SS - summary statistics vector of length (n+1)*m - //m - number of rows of SS, m = p*(p+1)/2 - - int pos; //position to fill in SS - int i,j,l; //loop indicies - pos = 0; - //Set first column of SS to zero - for(i = 0; i < *m; i++){ - SS[pos] = 0; - pos++; - } - - //For each new row in X, calculate the unique cross prod of X (lower tri) and - // add to previous summary statistic. - for(i = 0; i < *n; i++){ - for(j = 0; j < *nc; j++){ - for(l = j; l < *nc; l++){ - SS[pos] = SS[pos - *m] + X[j * *n + i]*X[l * *n + i]; - pos++; - } - } - } - - return; -} - -//Evaluate the regression quadratic cost function based on summary statistics -void RegQuadCostFunc(double *SS, int *m, int *n, int *p, int *start, int *end, - double *cost, double *tol, int *error, double *scale, int *MBIC){ - //SS - Summary statistics - //m - number of rows of SS - //n - number of columns of SS - //p - number of regressors, nb m = (p+1)*(p+2)/2 - //start - index at start of segment - //end - index at end of segment - //tol - tolerence for lm - //error - error index - //scale - if>0, cost=-2loglike @ fixed scale/variance - // if=0, cost=-2logLik() - // if<0, cost=RSS - //MBIC - 1 if MBIC penalty used, 0 if not. - - //-- memory allocation -- - //Summary statistics for segment - double *Sumstats = calloc(*m, sizeof(double)); - if(Sumstats == NULL){ - *error = 101; - goto err101; - } - //t(X) %*% X matrix (populated by Sumstats) - double *XX = calloc(*p * *p, sizeof(double)); - if(XX == NULL){ - *error = 102; - goto err102; - } - //t(X) %*% y vector (populated by Sumstats) - double *Xy = calloc(*p, sizeof(double)); - if(Xy == NULL){ - *error = 103; - goto err103; - } - //coefficients reordered according to pivot (see below) - double *beta = calloc(*p, sizeof(double)); - if(beta == NULL){ - *error = 104; - goto err104; - } - //copy of XX (converst to qr on return of dqrls) - double *qr = calloc(*p * *p, sizeof(double)); - if(qr == NULL){ - *error = 105; - goto err105; - } - //returns coefficient estimates of lm - double *coef = calloc(*p, sizeof(double)); - if(coef == NULL){ - *error = 106; - goto err106; - } - //returns auxilary information on qr decomposition - double *qraux = calloc(*p, sizeof(double)); - if(qraux == NULL){ - *error = 107; - goto err107; - } - //working space - double *work = calloc(2 * *p, sizeof(double)); - if(work == NULL){ - *error = 108; - goto err108; - } - //1:p (changes if cols of XX are permuted) - int *pivot = calloc(*p, sizeof(int)); - if(pivot == NULL){ - *error = 109; - goto err109; - } - //copy of Xy (converts to y-Xb) - double *residuals = calloc(*p, sizeof(double)); - if(residuals == NULL){ - *error = 110; - goto err110; - } - //copy of Xy (converts to orthog effects: t(q) %*% Xy) - double *effects = calloc(*p, sizeof(double)); - if(effects == NULL){ - *error = 111; - goto err111; - } - //copy of Xy - double *Xytmp = calloc(*p, sizeof(double)); - if(Xytmp == NULL){ - *error = 112; - goto err112; - } - - int i, j, tri; //indices - int rank = 0; //returns estimated rank of system (XX) - int ny = 1; //number of columns of Xy (here, always 1) - int nn = *end - *start; //number of observations - double RSS; //residual sum of squares - - //Evaluate summary statistics for segment - for(i = 0; i < *m; i++){ - //[yy, Xy, XX(lower.tri)] - Sumstats[i] = SS[*end * *m + i] - SS[*start * *m + i]; - } - - //Populate XX and Xy with values from Sumstats. - //Also, make copies and initialise values for dqrls - tri = 0; - for(i = 0; i < *p; i++){ - Xy[i] = Sumstats[i+1]; - Xytmp[i] = Xy[i]; //copy of Xy - coef[i] = 0.0; //initialise to zero - residuals[i] = Xy[i]; //copy of Xy - effects[i] = Xy[i]; //copy of Xy - pivot[i] = i+1; //index from 1:p - qraux[i] = 0.0; //initialise to zero - work[i] = 0.0; //initialise first half to zero - work[i + *p] = 0.0; //initialise second half to zero - for(j = i; j < *p; j++){ - XX[i * *p + j] = Sumstats[*p + 1 + tri + j]; - qr[i * *p + j] = XX[i * *p + j]; //copy of XX - } - tri += *p - i - 1; - for(j = 0; j < i; j++){ - XX[i * *p + j] = XX[j * *p + i]; //Copy to make XX symmetric - qr[i * *p + j] = XX[i * *p + j]; //Copy of XX - } - } - - //Call to routine used by lm to solve system of linear equations - F77_CALL(dqrls)(qr, p, p, Xytmp, &ny, tol, coef, residuals, effects, &rank , - pivot, qraux, work); - - //Extract coefficients in order - for(i = 0; i < *p; i++){ - beta[pivot[i]-1] = coef[i]; - } - - //Calculate the quadratic cost - RSS = Sumstats[0]; - for(i = 0; i < *p; i++){ - RSS -= 2 * beta[i] * Xy[i]; - for(j = 0; j < *p; j++){ - RSS += beta[i] * beta[j] * XX[i * *p + j]; - } - } - - if(*scale == 0){ - *cost = nn + nn*log(2 * M_PI * RSS) - nn*log(nn); //-2*logLik() - }else if(*scale > 0){ - *cost = nn*log(2 * M_PI * *scale) + (RSS / *scale); //-2LL(sig2=scale) - }else{ - *cost = RSS; //RSS - } - if(*MBIC==1){ - //*cost = RSS+log(nn); // extra log(length segment) for MBIC cost - } - //Free allocated memory - free(Xytmp); - err112: free(effects); - err111: free(residuals); - err110: free(pivot); - err109: free(work); - err108: free(qraux); - err107: free(coef); - err106: free(qr); - err105: free(beta); - err104: free(Xy); - err103: free(XX); - err102: free(Sumstats); - err101: return; - -} - -//Find the minimum case -void min_which(double *data, int *n, double *minval, int *minid){ - //data - values for which to find the minimum - //n - number of items to search - //minval - minimum value - //minid - index from start where to find minimum - int i; - *minid = 0; - *minval = data[*minid]; - for(i = 1; i < *n; i++){ - if(data[i] < *minval){ - *minid = i; - *minval = data[i]; - } - } - return; -} - - -//Determine cpts for Normal regression under method PELT -void CptReg_Normal_PELT(double *data, int *n, int *m, double *pen, int *cptsout, - int *error, double *shape, int *minseglen, double *tol, double *lastchangelike, - int *lastchangecpts, int *numchangecpts, int *MBIC){ - - //data - vectorised matrix of size n x m - //n - number of records - //m - number of data points per record - //pen - penalty value - //cptsout - position of estimated changepoints - //error - Index stating if an error has occured (error=0 if ok) - //shape - known variance (if=0,use MLE, if<0, cost=RSS) - //minseglen - minimum segment length - //lastchangelike - working space (cost value at last changepoint) - //lastchangecpts - working space (index of last changepoint) - //numchangecpts - working space (number of cpts that have occured so far) - //MBIC - 1 if MBIC penalty, 0 if not - - int p = *m - 1; //number or regressors - int np1 = *n + 1; //ncols of summary statistcs array - int size = (*m * (*m + 1)) * 0.5; //nrows of summary statistics array - int nchecklist, nchecktmp; //number of items in the checklist - double minout; //minimum cost value - int tstar, i, j, start; //indicies - int whichout; //index corresponding to the minimum cost value - double segcost; //cost over specified segment - *error = 0; - - //working space: position of last changepoint for all 'active' branches - int *checklist = (int *)calloc(np1, sizeof(int)); - if(checklist==NULL){ - *error = 1; - goto err1; - } - //working space: updated cost wrt latest observation for all 'active' branches - double *tmplike = (double *)calloc(np1, sizeof(double)); - if(tmplike == NULL){ - *error = 2; - goto err2; - } - //Summary statistcs - double *Sumstats = (double *)calloc(np1 * size, sizeof(double)); - if(Sumstats == NULL){ - *error = 3; - goto err3; - } - - //Evaluate the summary statistics - RegQuadCost_SS(data, n, m, Sumstats, &size); - - //Initialise - for(j = 0; j <= *minseglen; j++){ - if(j==0){ - lastchangelike[j] = -*pen; - }else{ - lastchangelike[j] = 0; - } - lastchangecpts[j] = 0; - numchangecpts[j] = 0; - } - //Evaluate cost for second minseglen - for(j = *minseglen+1; j <= (2 * *minseglen); j++){ - start = 0; - RegQuadCostFunc(Sumstats, &size, &np1, &p, &start, &j, lastchangelike + j, - tol, error, shape, MBIC); - if(*error != 0){ - goto err4; - } - lastchangecpts[j] = 0; - numchangecpts[j] = 1; - } - //setup check list; - nchecklist = 2; - checklist[0] = 0; - checklist[1] = *minseglen+1; - - //progress through time series - for(tstar = 2 * *minseglen+1; tstar < np1; tstar++){ - R_CheckUserInterrupt(); //Has interrupred the R session? quits if true. - - if(lastchangelike[tstar] == 0){ //Should always be true! - for(i = 0; i < nchecklist; i++){ - //Evaluate cost&penalty based on - // total cost&pen at last cpt + cost over current segment + penalty - start = checklist[i]; //last point of last segment - RegQuadCostFunc(Sumstats, &size, &np1, &p, &start, &tstar, &segcost, - tol, error, shape, MBIC); - if(*error != 0){ - goto err4; - } - tmplike[i] = lastchangelike[start] + segcost + *pen; - } - //Find and store branch with minimum cost&pen - min_which(tmplike, &nchecklist, &minout, &whichout); - lastchangelike[tstar] = minout; - lastchangecpts[tstar] = checklist[whichout]; - numchangecpts[tstar] = numchangecpts[lastchangecpts[tstar]] + 1; - - //Prune out non-minimum branches - nchecktmp = 0; - for(i = 0; i < nchecklist; i++){ - if(tmplike[i] <= (lastchangelike[tstar] +*pen)){ - checklist[nchecktmp] = checklist[i]; - nchecktmp++; - } - } - nchecklist = nchecktmp; - } - //Add new cpt to checklist - checklist[nchecklist] = tstar - *minseglen; - nchecklist++; - } - - //Extract optimal changepoint set - int ncpts = 0; - int last = *n; - while(last != 0){ - cptsout[ncpts] = last; - last = lastchangecpts[last]; - ncpts++; - } - - //Free allocated memory - err4: free(Sumstats); - err3: free(tmplike); - err2: free(checklist); - err1: return; -} - - -//Free allocated memory in case R session has been interrupred -void Free_CptReg_Normal_PELT(int *error){ - // Error code from CptReg_Normal_PELT C function, non-zero => error - if(*error==0){ - free((void *)checklist); - free((void *)tmplike); - free((void *)Sumstats); - } -} - -//Determine cpts for Normal regression under method AMOC -void CptReg_Normal_AMOC(double *data, int *n, int *m, double *pen, - int *error, double *shape, int *minseglen, double *tol, int *tau, - double *nulllike, double *taulike, double *tmplike, int *MBIC){ - - //data - vectorised matrix of size n x m - //n - number of records - //m - number of data points per record - //pen - penalty value - //error - Index stating if an error has occured (error=0 if ok) - //shape - known variance (if=0,use MLE, if<0, cost=RSS) - //minseglen - minimum segment length - //tau - estimated single changepoint position - //nulllike - estimated cost of no changepoints - //taulike - estimated cost of single changepoint - //tmplike - working memory: store all single changepoint costs - //MBIC - 1 if MBIC penalty, 0 if not - - int p = *m - 1; //number or regressors - int np1 = *n + 1; //ncols of summary statistcs array - int size = (*m * (*m + 1)) * 0.5; //nrows of summary statistics array - int tstar; - double seg1cost, seg2cost; - int zero; - int neval; - *error = 0; - - //Summary statistcs - double *Sumstats = (double *)calloc(np1 * size, sizeof(double)); - if(Sumstats == NULL){ - *error = 1; - goto err1; - } - //Evaluate the summary statistics - RegQuadCost_SS(data, n, m, Sumstats, &size); - - //cost with no changepoints - zero = 0; - RegQuadCostFunc(Sumstats, &size, &np1, &p, &zero, n, nulllike, tol, error, shape, MBIC); - if(*error != 0){ - goto err2; - } - - //cost with at most one changepoint - neval = 0; - for(tstar = *minseglen; tstar <= (*n - *minseglen); tstar++){ //?? - R_CheckUserInterrupt(); //Has interrupred the R session? quits if true. - RegQuadCostFunc(Sumstats, &size, &np1, &p, &zero, &tstar, &seg1cost, - tol, error, shape, MBIC); - if(*error != 0){ - goto err2; - } - RegQuadCostFunc(Sumstats, &size, &np1, &p, &tstar, n, &seg2cost, - tol, error, shape, MBIC); - if(*error != 0){ - goto err2; - } - tmplike[tstar-1] = seg1cost + seg2cost; - neval++; - } - //Which tstar returns the minimum cost - min_which(tmplike + *minseglen - 1, &neval, taulike, tau); - *tau += *minseglen; - - //Free allocated memory - err2: free(Sumstats); - err1: return; -} - -//Free allocated memory in case R session has been interrupred -void Free_CptReg_Normal_AMOC(int *error){ - // Error code from CptReg_Normal_AMOC C function, non-zero => error - if(*error==0){ - free((void *)Sumstats); - } -} - - - - - - - - - - - - From c779f27bac023f07691eda780432e6a20c9627fe Mon Sep 17 00:00:00 2001 From: AndrewC1998 Date: Wed, 5 Jun 2019 10:19:03 +0100 Subject: [PATCH 2/4] 05/06/2019 Namespace amended to match changes. AC --- NAMESPACE | 1 - src/Makevars | 1 - 2 files changed, 2 deletions(-) delete mode 100755 src/Makevars diff --git a/NAMESPACE b/NAMESPACE index e5c15cc..33123ce 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,3 @@ -useDynLib(EnvCpt) import(methods) import(graphics) importFrom(MASS,fitdistr) diff --git a/src/Makevars b/src/Makevars deleted file mode 100755 index 22ebc63..0000000 --- a/src/Makevars +++ /dev/null @@ -1 +0,0 @@ -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) From d6dfca1a8f17655fbaf656088d081146c2658ed2 Mon Sep 17 00:00:00 2001 From: Andrew Connell <37156349+AndrewC1998@users.noreply.github.com> Date: Thu, 4 Jul 2019 16:26:33 +0100 Subject: [PATCH 3/4] Update DESCRIPTION Slight changes to increase changepoint package version dependency. --- DESCRIPTION | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 95d5207..ddb466d 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,17 +3,17 @@ Type: Package Title: Detection of Structural Changes in Climate and Environment Time Series Version: 1.1.1 -Date: 2018-09-13 +Date: 2019-07-04 Authors@R: c(person("Rebecca", "Killick", role=c("aut","cre"),email="r.killick@lancs.ac.uk"), - person("Claudie", "Beaulieu", role="aut"), person("Simon", "Taylor", role="aut"),person("Harjit","Hullait",role="aut")) +person("Claudie", "Beaulieu", role="aut"), person("Simon", "Taylor", role="aut"),person("Harjit","Hullait",role="aut"), person("Andrew", "Connell", role=c("ctbr"))) Maintainer: Rebecca Killick URL: https://github.com/rkillick/EnvCpt/ Description: Tools for automatic model selection and diagnostics for Climate and Environmental data. In particular the envcpt() function does automatic model selection between a variety of trend, changepoint and autocorrelation models. The envcpt() function should be your first port of call. -Depends: R(>= 3.3), changepoint(>= 2.2.2), grDevices, MASS, +Depends: R(>= 3.3), changepoint(>= 2.3.1), grDevices, MASS, methods, stats, utils, zoo Suggests: testthat License: GPL LazyLoad: yes -Packaged: 2018-09-13 13:46:06 UTC; killick +Packaged: 2019-07-04 16:23:06 UTC; killick NeedsCompilation: yes From 87fa6306fe76da5fd14c8b544994f8d9052e64f1 Mon Sep 17 00:00:00 2001 From: AndrewC1998 Date: Thu, 22 Aug 2019 19:47:59 +0100 Subject: [PATCH 4/4] Updated description --- DESCRIPTION | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 95d5207..e17baa4 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,9 +3,9 @@ Type: Package Title: Detection of Structural Changes in Climate and Environment Time Series Version: 1.1.1 -Date: 2018-09-13 +Date: 2019-09-22 Authors@R: c(person("Rebecca", "Killick", role=c("aut","cre"),email="r.killick@lancs.ac.uk"), - person("Claudie", "Beaulieu", role="aut"), person("Simon", "Taylor", role="aut"),person("Harjit","Hullait",role="aut")) + person("Claudie", "Beaulieu", role="aut"), person("Simon", "Taylor", role="aut")) Maintainer: Rebecca Killick URL: https://github.com/rkillick/EnvCpt/ Description: Tools for automatic model selection and diagnostics for Climate and Environmental data. In particular the envcpt() function does automatic model selection between a variety of trend, changepoint and autocorrelation models. The envcpt() function should be your first port of call. @@ -14,6 +14,5 @@ Depends: R(>= 3.3), changepoint(>= 2.2.2), grDevices, MASS, Suggests: testthat License: GPL LazyLoad: yes -Packaged: 2018-09-13 13:46:06 UTC; killick +Packaged: 2019-09-22 19:46:06 UTC; killick NeedsCompilation: yes -