diff --git a/DESCRIPTION b/DESCRIPTION index 95d5207..4f78c45 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,17 +3,16 @@ 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. -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-09-22 19:46:06 UTC; killick NeedsCompilation: yes - 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/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); - } -} - - - - - - - - - - - - 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)