rm(list=ls()) library(foreach) library(doParallel) library(RSpectra) library(gsl) # library(lattice) #### code for the nonparametric isotropic model used ### input are: simstimes - (ordered) times of the events ### xcoords - x coordinates of events ### ycoords - y coordinates of events ### itermax - number of iterations to run the em-algorithm for ### horizon - time horizon for events ### backgroundbw - the bandwidth of the two dimensional Gaussian kernel used to estimate the background rate ### nearneigh - the number of nearest neighbours used to select $D_i$ in the triggering function ### maxt - maximum $t$ at which the triggering function can have an effect ### maxr - maximum distance at which the triggering function can have an effect. ### output of the function: a vector which gives the estimated probability the input event was a background event ### the times at which an event is said to be triggered ### the euclidean distance at which an event is said to have been triggered. ### the raw distance in the x-direction at which an event is said to be triggered ### the raw distance in the y-direction at which an event is said to be triggered. EMkdeisotropictrig_fixedbwadaptedbg <- function(simstimes, xcoords, ycoords, itermax=100, horizon, backgroundbw=0.15, nearneigh=15, maxt=50, maxr=1){ iter <- 0 p2 <- matrix(0,length(simstimes), length(simstimes)) ###### matrix with time differences for (k in 1:length(simstimes)){ for (l in 1:k){ if (k > l){ p2[k,l] <- simstimes[k] - simstimes[l] } else{ p2[k,l] <- 0 } } } p3 <- matrix(0,length(simstimes), length(simstimes)) ###### matrix with euclidean distances for (k in 1:length(simstimes)){ for (l in 1:k){ if (k > l){ p3[k,l] <- sqrt((xcoords[k] - xcoords[l])^2 + (ycoords[k] - ycoords[l])^2) } else{ p3[k,l] <- 0 } } } omega=0.1 sigmax=0.3 initp <- matrix(0,length(simstimes), length(simstimes)) for (k in 1:length(simstimes)){ for (l in 1:k){ if (k == l){ initp[k,l] <- 1 } else if ((k > l) ){ initp[k,l] <- exp(-omega*p2[k,l])*exp((-p3[k,l]^2/(2*sigmax^2)))/(2*pi*sigmax^2) } } } dividep <- numeric(length(simstimes)) for (i in 1:length(dividep)){ dividep[i] <- sum(initp[i,]) } p <- matrix(0,length(simstimes), length(simstimes)) for (k in 1:length(simstimes)){ for (l in 1:k){ p[k,l] <- initp[k,l]/dividep[k] } } p_x <- matrix(0,length(simstimes), length(simstimes)) ###### matrix with x differences for (k in 1:length(simstimes)){ for (l in 1:k){ if (k > l){ p_x[k,l] <- xcoords[k] - xcoords[l] } else{ p_x[k,l] <- 0 } } } p_y <- matrix(0,length(simstimes), length(simstimes)) ###### matrix with y differences for (k in 1:length(simstimes)){ for (l in 1:k){ if (k > l){ p_y[k,l] <- ycoords[k] - ycoords[l] } else{ p_y[k,l] <- 0 } } } while((iter < itermax) & (sum(diag(p)) < (0.999*length(simstimes)))){ index <- 1:length(simstimes) sampleindex <- numeric(length(simstimes)) backgroundtimes <- NULL backgroundx <- NULL backgroundy <- NULL trigtimes <- NULL trigdist <- NULL trigxdir <- NULL trigydir <- NULL backnum <- 0 trignum <- 0 for (i in 1:length(simstimes)){ samp <- sample(index, size=1, prob=p[i,]) sampleindex[i] <- samp if (samp == i){ backnum <- backnum + 1 backgroundtimes[backnum] <- simstimes[i] backgroundx[backnum] <- xcoords[i] backgroundy[backnum] <- ycoords[i] } else{ trignum <- trignum + 1 trigtimes[trignum] <- p2[i, samp] trigdist[trignum] <- p3[i, samp] trigxdir[trignum] <- p_x[i, samp] trigydir[trignum] <- p_y[i, samp] } } ###### to find nearest neighbours for triggered times, we scale them to have unit variance trigrate <- trignum/length(simstimes) expecttrig <- (length(simstimes) - sum(diag(p)))/(length(simstimes)) nn2 <- min((nearneigh + 1), trignum) sdtrigtimes <- sd(trigtimes) sdtrigdist <- sd(trigdist) scaledtrigtimes <- trigtimes/sdtrigtimes scaledtrigdist <- trigdist/sdtrigdist scaleddistance <- matrix(0,length(trigtimes), length(trigtimes)) ###### matrix to find scaled distance for nearest neighbour for (k in 1:length(trigtimes)){ for (l in 1:length(trigtimes)){ scaleddistance[k,l] <- sqrt((scaledtrigtimes[k] - scaledtrigtimes[l])^2 + (scaledtrigdist[k] - scaledtrigdist[l])^2 ) } } D_trig <- numeric(length(trigtimes)) for (i in 1:length(D_trig)){ D_trig[i] <- scaleddistance[i,][order(scaleddistance[i,])[nn2]] } normconst <- numeric(length(trigtimes)) for (i in 1:length(normconst)){ normconst[i] <- 0.5*(1 + erf(trigdist[i]/(sqrt(2)*D_trig[i]*sdtrigdist))) } g <- matrix(0,length(simstimes), length(simstimes)) ###### matrix with time differences numCores <- detectCores() registerDoParallel(numCores) gparallel <- foreach(k=1:length(simstimes)) %dopar% { grow <- numeric(length(simstimes)) for (l in 1:k){ if ((k > l) & (p2[k,l] > 0) & (p2[k,l] < maxt) & (abs(p3[k,l]) < maxr) & (abs(p3[k,l])>0)){ grow[l] <- expecttrig*sum(exp(-((p3[k,l]-trigdist)^2)/(2*(sdtrigdist^2)*(D_trig^2)))*(exp(-((p2[k,l]-trigtimes)^2)/(2*(sdtrigtimes^2)*(D_trig^2))) + exp(-((-p2[k,l]-trigtimes)^2)/(2*(sdtrigtimes^2)*(D_trig^2))))/((D_trig^2)*normconst*2*pi*p3[k,l])) } } grow } stopImplicitCluster() for (i in 1:length(simstimes)){ g[i,] <- gparallel[[i]]/(trignum*sdtrigtimes*sdtrigdist*((2*pi))) } rm(gparallel) numCores <- detectCores() registerDoParallel(numCores) bgpar <- foreach(i=1:length(simstimes)) %dopar% { bgrate <- 0 for (j in 1:length(simstimes)){ if ((xcoords[i] != xcoords[j]) | (ycoords[i] != ycoords[j])){ bgrate <- bgrate + p[j,j]*exp(-((xcoords[i]-xcoords[j])^2)/(2*backgroundbw^2) -((ycoords[i]-ycoords[j])^2)/(2*backgroundbw^2))/(backgroundbw^2) } } bgrate } stopImplicitCluster() bgrate <- numeric(length(simstimes)) for (i in 1:length(simstimes)){ bgrate[i] <- bgpar[[i]]/(horizon*2*pi) } rm(bgpar) dividep <- numeric(length(simstimes)) for (i in 1:length(dividep)){ dividep[i] <- sum(g[i,]) } numCores <- detectCores() registerDoParallel(numCores) ppar <- foreach(i=1:length(simstimes)) %dopar% { pnum <- numeric(length(simstimes)) for (j in 1:length(simstimes)){ if (i==j){ pnum[j] <- bgrate[i]/(bgrate[i] + dividep[i]) } else if (i > j){ pnum[j] <- g[i,j]/(bgrate[i] + dividep[i]) } } pnum } stopImplicitCluster() p <- matrix(0, length(simstimes), length(simstimes)) for (i in 1:length(simstimes)){ p[i,] <- ppar[[i]] } expecttrig <- (length(simstimes) - sum(diag(p)))/(length(simstimes)) rm(ppar) iter <- iter + 1 print(iter) print(expecttrig) } newlist <- list(diag(p), trigtimes, trigdist, trigxdir, trigydir) return(newlist) } #### code for the etas model with gaussian spatial components ### input are: simstimes - (ordered) times of the events ### xcoords - x coordinates of events ### ycoords - y coordinates of events ### itermax - number of iterations to run the em-algorithm for ### horizon - time horizon for events ### initguess - values relating to an initial starting guess ### maxt - maximum $t$ at which the triggering function can have an effect ### maxr - maximum distance at which the triggering function can have an effect. ### backgroundbw - the bandwidth of the two dimensional Gaussian kernel used to estimate the background rate ### output of the function: a vector which gives the estimated probability the input event was a background event ### estimated alpha value ### estimated omega value ### estimated sigma value etasgaussianspacefixedbgbw <- function(simstimes, xcoords, ycoords, itermax=100, horizon, initguess=c(0.2,0.2), maxt=90, maxr=1, backgroundbw=0.15){ iter <- 0 timemat <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) timeparallel <- foreach(k=1:length(simstimes)) %dopar% { timerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ timerow[l] <- simstimes[l] - simstimes[k] } timerow } stopImplicitCluster() for (i in 1:length(simstimes)){ timemat[i,] <- timeparallel[[i]] } rm(timeparallel) distancemat <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) distanceparallel <- foreach(k=1:length(simstimes)) %dopar% { distancerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ distancerow[l] <- ((xcoords[k] - xcoords[l])^2 + (ycoords[k] - ycoords[l])^2) } distancerow } stopImplicitCluster() for (i in 1:length(simstimes)){ distancemat[i,] <- distanceparallel[[i]] } omega <- initguess[1] sigma <- initguess[2] inittrigalpha <- 0.3 initbgalpha <- 0.7 inittrigp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) poneparallel <- foreach(k=1:length(simstimes)) %dopar% { ponerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (distancemat[k,l]>0)){ ponerow[l] <- inittrigalpha*omega*exp(-omega*(timemat[k,l]))*exp(-distancemat[k,l]/(2*sigma^2))/(2*pi*sigma^2) } } ponerow } stopImplicitCluster() for (i in 1:length(simstimes)){ inittrigp[i,] <- poneparallel[[i]] } rm(poneparallel) initbgp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwoparallel <- foreach(k=1:length(simstimes)) %dopar% { ptworow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((distancemat[k,l]>0)){ ptworow[l] <- initbgalpha*exp(-distancemat[k,l]/(2*backgroundbw^2))/(2*pi*backgroundbw^2*horizon) } } ptworow } stopImplicitCluster() for (i in 1:length(simstimes)){ initbgp[i,] <- ptwoparallel[[i]] } rm(ptwoparallel) lambdatotal <- numeric(length(simstimes)) for (i in 1:length(lambdatotal)){ lambdatotal[i] <- sum(inittrigp[,i]) + sum(initbgp[,i]) } inittrigpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ponefinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ponefinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (distancemat[k,l]>0)){ ponefinalrow[l] <- inittrigp[k,l]/(lambdatotal[l]) } } ponefinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ inittrigpfinal[i,] <- ponefinalparallel[[i]] } rm(ponefinalparallel) initbgpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwofinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ptwofinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ ptwofinalrow[l] <- initbgp[k,l]/(lambdatotal[l]) } ptwofinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ initbgpfinal[i,] <- ptwofinalparallel[[i]] } rm(ptwofinalparallel) bgptotal <- numeric(length(simstimes)) for (i in 1:length(bgptotal)){ bgptotal[i] <- sum(initbgpfinal[,i]) } omega <- sum(inittrigpfinal)/(sum(inittrigpfinal*timemat)) sigmanewsq <- (sum(inittrigpfinal*distancemat))/(2*sum(inittrigpfinal)) sigma <- sqrt(sigmanewsq) alpha <- sum(inittrigpfinal)/length(simstimes) while(iter < itermax){ trigp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) poneparallel <- foreach(k=1:length(simstimes)) %dopar% { ponerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (timemat[k,l] < maxt) & (distancemat[k,l]>0) & (distancemat[k,l] < (maxr^2))){ ponerow[l] <- alpha*omega*exp(-omega*(timemat[k,l]))*exp(-distancemat[k,l]/(2*sigma^2))/(2*pi*sigma^2) } } ponerow } stopImplicitCluster() for (i in 1:length(simstimes)){ trigp[i,] <- poneparallel[[i]] } rm(poneparallel) bgp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwoparallel <- foreach(k=1:length(simstimes)) %dopar% { ptworow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((distancemat[k,l]>0)){ ptworow[l] <- bgptotal[l]*exp(-distancemat[k,l]/(2*backgroundbw^2))/(2*pi*backgroundbw^2*horizon) ####need to check if it's k or l } } ptworow } stopImplicitCluster() for (i in 1:length(simstimes)){ bgp[i,] <- ptwoparallel[[i]] } rm(ptwoparallel) lambdatotal <- numeric(length(simstimes)) for (i in 1:length(lambdatotal)){ lambdatotal[i] <- sum(trigp[,i]) + sum(bgp[,i]) } trigpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ponefinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ponefinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (distancemat[k,l]>0)){ ponefinalrow[l] <- trigp[k,l]/(lambdatotal[l]) } } ponefinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ trigpfinal[i,] <- ponefinalparallel[[i]] } rm(ponefinalparallel) omega <- sum(trigpfinal)/(sum(trigpfinal*timemat)) sigmanewsq <- (sum(trigpfinal*distancemat))/(2*sum(trigpfinal)) sigma <- sqrt(sigmanewsq) alpha <- sum(trigpfinal)/length(simstimes) iter <- iter + 1 print(c(alpha,omega,sigma)) } results <- list(bgptotal,alpha,omega,sigma) return(results) } #### code for the dr model with gaussian spatial components ### input are: simstimes - (ordered) times of the events ### xcoords - x coordinates of events ### ycoords - y coordinates of events ### itermax - number of iterations to run the em-algorithm for ### horizon - time horizon for events ### initguess - values relating to an initial starting guess ### maxt - maximum $t$ at which the triggering function can have an effect ### maxr - maximum distance at which the triggering function can have an effect. ### backgroundbw - the bandwidth of the two dimensional Gaussian kernel used to estimate the background rate ### output of the function: a vector which gives the estimated probability the input event was a background event ### estimated alpha value ### estimated omega value ### estimated sigma value dcgaussianspacefixedbgbw <- function(simstimes, xcoords, ycoords, itermax=100, horizon, initguess=c(0.2,0.2), maxt=90, maxr=1, backgroundbw=0.15){ iter <- 0 timemat <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) timeparallel <- foreach(k=1:length(simstimes)) %dopar% { timerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ timerow[l] <- simstimes[l] - simstimes[k] } timerow } stopImplicitCluster() for (i in 1:length(simstimes)){ timemat[i,] <- timeparallel[[i]] } rm(timeparallel) distancemat <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) distanceparallel <- foreach(k=1:length(simstimes)) %dopar% { distancerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ distancerow[l] <- ((xcoords[k] - xcoords[l])^2 + (ycoords[k] - ycoords[l])^2) } distancerow } stopImplicitCluster() for (i in 1:length(simstimes)){ distancemat[i,] <- distanceparallel[[i]] } omega <- initguess[1] sigma <- initguess[2] inittrigalpha <- 0.3 initbgalpha <- 0.7 inittrigp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) poneparallel <- foreach(k=1:length(simstimes)) %dopar% { ponerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (distancemat[k,l]>0)){ ponerow[l] <- inittrigalpha*omega^2*timemat[k,l]*exp(-omega*(timemat[k,l]))*exp(-distancemat[k,l]/(2*sigma^2))/(2*pi*sigma^2) } } ponerow } stopImplicitCluster() for (i in 1:length(simstimes)){ inittrigp[i,] <- poneparallel[[i]] } rm(poneparallel) initbgp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwoparallel <- foreach(k=1:length(simstimes)) %dopar% { ptworow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((distancemat[k,l]>0)){ ptworow[l] <- initbgalpha*exp(-distancemat[k,l]/(2*backgroundbw^2))/(2*pi*backgroundbw^2*horizon) } } ptworow } stopImplicitCluster() for (i in 1:length(simstimes)){ initbgp[i,] <- ptwoparallel[[i]] } rm(ptwoparallel) lambdatotal <- numeric(length(simstimes)) for (i in 1:length(lambdatotal)){ lambdatotal[i] <- sum(inittrigp[,i]) + sum(initbgp[,i]) } inittrigpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ponefinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ponefinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (distancemat[k,l]>0)){ ponefinalrow[l] <- inittrigp[k,l]/(lambdatotal[l]) } } ponefinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ inittrigpfinal[i,] <- ponefinalparallel[[i]] } rm(ponefinalparallel) initbgpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwofinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ptwofinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ ptwofinalrow[l] <- initbgp[k,l]/(lambdatotal[l]) } ptwofinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ initbgpfinal[i,] <- ptwofinalparallel[[i]] } rm(ptwofinalparallel) bgptotal <- numeric(length(simstimes)) for (i in 1:length(bgptotal)){ bgptotal[i] <- sum(initbgpfinal[,i]) } omega <- 2*sum(inittrigpfinal)/(sum(inittrigpfinal*timemat)) sigmanewsq <- (sum(inittrigpfinal*distancemat))/(2*sum(inittrigpfinal)) sigma <- sqrt(sigmanewsq) alpha <- sum(inittrigpfinal)/length(simstimes) while(iter < itermax){ trigp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) poneparallel <- foreach(k=1:length(simstimes)) %dopar% { ponerow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (timemat[k,l] < maxt) & (distancemat[k,l]>0) & (distancemat[k,l] < (maxr^2))){ ponerow[l] <- alpha*omega^2*timemat[k,l]*exp(-omega*(timemat[k,l]))*exp(-distancemat[k,l]/(2*sigma^2))/(2*pi*sigma^2) } } ponerow } stopImplicitCluster() for (i in 1:length(simstimes)){ trigp[i,] <- poneparallel[[i]] } rm(poneparallel) bgp <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwoparallel <- foreach(k=1:length(simstimes)) %dopar% { ptworow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((distancemat[k,l]>0)){ ptworow[l] <- bgptotal[l]*exp(-distancemat[k,l]/(2*backgroundbw^2))/(2*pi*backgroundbw^2*horizon) ####need to check if it's k or l } } ptworow } stopImplicitCluster() for (i in 1:length(simstimes)){ bgp[i,] <- ptwoparallel[[i]] } rm(ptwoparallel) lambdatotal <- numeric(length(simstimes)) for (i in 1:length(lambdatotal)){ lambdatotal[i] <- sum(trigp[,i]) + sum(bgp[,i]) } trigpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ponefinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ponefinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ if ((timemat[k,l]>0) & (distancemat[k,l]>0)){ ponefinalrow[l] <- trigp[k,l]/(lambdatotal[l]) } } ponefinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ trigpfinal[i,] <- ponefinalparallel[[i]] } rm(ponefinalparallel) bgpfinal <- matrix(0,length(simstimes), length(simstimes)) numCores <- detectCores() registerDoParallel(numCores) ptwofinalparallel <- foreach(k=1:length(simstimes)) %dopar% { ptwofinalrow <- numeric(length(simstimes)) for (l in 1:length(simstimes)){ ptwofinalrow[l] <- bgp[k,l]/(lambdatotal[l]) } ptwofinalrow } stopImplicitCluster() for (i in 1:length(simstimes)){ bgpfinal[i,] <- ptwofinalparallel[[i]] } rm(ptwofinalparallel) bgptotal <- numeric(length(simstimes)) for (i in 1:length(bgptotal)){ bgptotal[i] <- sum(bgpfinal[,i]) } omega <- 2*sum(trigpfinal)/(sum(trigpfinal*timemat)) sigmanewsq <- (sum(trigpfinal*distancemat))/(2*sum(trigpfinal)) sigma <- sqrt(sigmanewsq) alpha <- sum(trigpfinal)/length(simstimes) iter <- iter + 1 print(c(alpha,omega,sigma,iter)) } results <- list(bgptotal,alpha,omega,sigma) return(results) }