# R/SAPP.R In SAPP: Statistical Analysis of Point Processes

#### Documented in eptrenlinsimpgraphptspec

```eptren <- function(data, mag = NULL, threshold = 0.0, nparam, nsub, cycle = 0,
tmpfile = NULL, nlmax = 1000, plot = TRUE)
{

# data                   #  point process data
n <- length(data)      #  total number of variable
# mag                    #  magnitude
# threshold              #  threshold magnitude
# nparam                 #  maximum number of parameters
# nsub                   #  number of subdivisions in either (0,t) or (0,cycle)
#  for the numerical integration of intensity.
# cycle                  # > 0 : periodicity to be investigated in days
if (cycle == 0)
nfunct = 1     #  trend  (exponential polynomial trend fitting)
if (cycle > 0)
nfunct = 2     #  cycle  (exponential fourier series fitting)

nn <- 0
xx <- NULL
if (is.null(mag)) {
xx <- data
magnitude <- data
nn <- n
t <- data[nn]
if (nfunct == 2)
for (i in 1:nn)
magnitude[i] <- magnitude[i] - as.integer(data[i]/t) * t
} else {
magnitude <- NULL
for (i in 1:n)
if ((mag[i] >= threshold) && (data[i] >= 0)) {
nn <- nn + 1
xx <- c(xx, data[i])
magnitude <- c(magnitude, mag[i])
}
t <- data[nn]
}

nmax <- nparam
if (nfunct == 2)
nmax <- 2 * nparam - 1
np <- 101
nlm <- nlmax
if (is.null(tmpfile))
nlm <- 0

z <- .Call("EptrenC",
as.double(xx),
as.double(t),
as.integer(nn),
as.integer(nfunct),
as.integer(nparam),
as.integer(nsub),
as.double(cycle),
as.integer(nmax),
as.integer(np),
as.integer(nlm))

xa <- array(z[[1L]], dim = c(nmax, nparam))
param <- list()
for (i in 1:nparam) {
m <- i
if (nfunct == 2)
m <- m * 2 - 1
param[[i]] <- xa[1:m, i]
}

if (plot == TRUE) {
par(mfrow = c(2, 1))
data1 <-matrix(, nrow = nn, ncol = 2)
data1[, 1] <- xx
data1[, 2] <- magnitude
if (is.null(mag))
data1[1:nn, 2] <- rep(1.0, nn)
data2 <- matrix(, nrow = np, ncol = 2)
data2[, 1] <- z[[5L]]
data2[, 2] <- z[[6L]]
if (cycle == 0) {
plot(data2, type = "l", main = "EPTREN-Trend", xlab = "Time",
ylab = "Intensity Rate")
plot(data1, type = "h", main = "", xlab = "Time", ylab = "Magnitude")
} else if (cycle > 0) {
for (i in 1:nn)
data1[i, 1] <- data1[i, 1] - as.integer(data1[i, 1] / cycle) * cycle
plot(data2, type = 'l', main = "EPTREN > EXP(CYCLE)",
xlab = "Superposed Occurrence Times", ylab = "Intensity Rate")
plot(data1, type = "h", main = "", xlab = "Superposed Occurrence Times",
ylab = "Magnitude")
}
par(mfrow = c(1, 1))
}

if (nlm > 0) {
nl <- z[[12L]]
x <- array(z[[7L]], c(nmax, nparam))
g <-array(z[[8L]], c(nmax, nparam))
id <- z[[9L]][1:nl]
ramda <- z[[10L]][1:nl]
ee <- z[[11L]][1:nl]
print.process(nfunct, nparam, x, g, id, ramda, ee, tmpfile)
}

eptren.out <- list(aic = z[[2L]], param = param, aicmin = z[[3L]],
maice.order = z[[4L]], time = z[[5L]], intensity = z[[6L]])
class(eptren.out) <- "eptren"
return(eptren.out)
}

print.eptren <- function(x, ...)
{

nparam <- length(x\$aic)
for (i in 1:nparam) {
message(gettextf("\n AIC\t%f", x\$aic[i]), domain = NA)
message(" parameters", domain = NA)
ii <- length(x\$param[[i]])
pmsg <- NULL
for (j in 1:ii) {
pmsg <- paste(pmsg, gettextf("  %e", x\$param[[i]][j]))
if ((j%%5 == 0) || (j == ii)) {
message(pmsg, domain = NA)
pmsg <- NULL
}
}
}

message(gettextf("\n\n minimum AIC\t%f", x\$aicmin), domain = NA)
message(" parameters", domain = NA)
i <- x\$maice.order
ii <- length(x\$param[[i]])
pmsg <- NULL
for (j in 1:ii) {
pmsg <- paste(pmsg, gettextf("  %e",x\$param[[i]][j]))
if ((j%%5 == 0) || (j == ii)) {
message(pmsg, domain = NA)
pmsg <- NULL
}
}
}

pgraph <- function(data, mag, threshold = 0.0, h, npoint, days, delta = 0.0,
dmax = 0.0, separate.graphics = FALSE)
{
nfunct <- 0
isi <- 0
n1 <- length(data)     #  total number of variable

if (delta == 0.0 || dmax == 0.0) {
dmax <- data[n1] / 4
delta <- dmax / 100
}
td <- data[n1] / delta
kmax <- as.integer(td/16.0 + 2)

zd <- rep(0, n1)
xmg <- rep(0, n1)
xmg1 <- 0.0
xmg2 <- 10.0
nn <- 0
for (i in 1:n1)
if (mag[i]==threshold || mag[i]>threshold)
if (mag[i]==xmg2 || mag[i]<xmg2)
if (mag[i]==xmg1 || mag[i]>xmg1)
if (data[i]== 0.0 || data[i]>0.0) {
nn <- nn + 1
zd[nn] <- data[i]
xmg[nn] <- mag[i]
#        if (xmg[nn]== 0.0 ) xmg[nn]=6.0
}
zd <- zd[1:nn]
xmg <- xmg[1:nn]
nn1 = nn - 1

z <- .Call("PgraphC",
as.integer(nfunct),
as.integer(isi),
as.double(zd),
as.integer(nn),
as.integer(npoint),
as.double(days),
as.double(h),
as.double(delta),
as.double(dmax),
as.integer(kmax) )

kn <- z[[3L]]
k <- z[[16L]]

tn <- data[n1] / as.double(nn)
xtau <- z[[1L]][1:kn]
y <- z[[2L]][1:kn]
xl <- z[[4L]][1:nn1] * tn
xx <- array(z[[5L]], dim = c(nn1, 6)) * tn
xx <- xx[1:nn1, 1:6]
ydev <- z[[6L]][1:nn1]
ui <- z[[7L]][1:nn1]
ncn <- z[[8L]][1:nn1]
sui <- z[[9L]][1:nn1]
xp <- z[[10L]]
xrate <- z[[11L]][1:npoint]
dlt <- z[[12L]]
xtime <- z[[13L]][1:k]
yvar <- array(z[[14L]], dim = c(5, kmax))
yvar <- yvar[1:5, 1:k]
sigma <- z[[15L]][1:k]

plot.pgraph(zd, xmg, h, kn, xtau, y, xl, xx, ydev, ui, ncn, sui, xp, npoint,
dlt, xrate, k, xtime, sigma, yvar, separate.graphics)

pgraph.out <- list(cnum = zd, lintv = xl, tau = xtau, nevent = y,
survivor = xx, deviation = ydev, normal.cnum = ncn,
normal.lintv = ui, success.intv = sui, occur = xrate,
time = xtime, variance = sigma, error = yvar)
return(pgraph.out)
}

plot.pgraph <- function(zd, xmg, h, kn, xtau, y, xl, xx, ydev, ui, ncn, sui, xp,
npoint, dlt, xrate, k, xtime, sigma, yvar,
separate.graphics)
{
nn <- length(zd)
nn1 <- nn-1
rnn <- as.double(nn)
ymax <- nn
xmax <- zd[nn]
nband1 <- 1.35810 * sqrt(rnn)
nband2 <- 1.62762 * sqrt(rnn)

# r.pgCumMT
par(mfrow=c(2, 1))
plot(x = zd, y = c(1:nn), type = "s", xlim = c(0, xmax), ylim = c(0, ymax),
main = "Cumulative Curve & M-T plot",
xlab = "Time", ylab = "Cumulative Number")
points(c(0, xmax), c(0, ymax), lty = 2, type = "l")
points(c(0, xmax), c(0, ymax)+nband1, col = 2, lty = 2, type = "l")
points(c(0, xmax), c(0, ymax)+nband2, col = 2, lty = 2, type = "l")
points(c(0, xmax), c(0, ymax)-nband1, col = 2, lty = 2, type = "l")
points(c(0, xmax), c(0, ymax)-nband2, col = 2, lty = 2, type = "l")

plot(x = zd, y = xmg, type = "h", xlim = c(0, xmax), xlab = "Time",
ylab = "Magnitude")

# r.pgPTnum
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(2, 1))
plot(x = xtau, y = y, type = "l", xlim = c(0, nn),
xlab = "tau = Time x (Total Number of Events) / (Time End)",
ylab = "Number of Events in [tau, tau+h]")
abline(h = seq(-3, 3, 1), lty = 2, col = 2)
title(main = paste("Normalized number of point in [t, t+h], where h=", h, sep = ""))
plot(x = zd, y = xmg, type = "h", xlim = c(0, xmax),
xlab = "Ordinary or Transformed Time", ylab = "Magnitude")

# r.pgSurviv
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(2, 1))
yy <- log10(c(nn1:1))
plot(x = xl, y = yy, type = "p", pch = "+", cex = 0.5,
main="Survivor Function", xlab = "Interval Length",
ylab = "Cumulative Number", axes = F)
for (j in 1:6)
points(xx[, j], yy, type = "l", col = 2)
box()
axis(1)
n9 <- c(1:9)
axis(2, at = log10(c(n9, n9*10, n9*10^2, n9*10^3, n9*10^4, 10^5)),
labels = rep("", 46))
axis(2, at = log10(c(1, 10, 100, 1000, 10000, 10000)),
labels = c(expression(10^0), expression(10^1), expression(10^2),
expression(10^3), expression(10^4), expression(10^5)))

xx <- c(1:nn1)
plot(x = xx, y = ydev, type = "n",
main = "Deviation of Survivor Function from the Poisson",
xlab = "Interval Length", ylab = "Deviations")
points(xx, ydev, type = "p", pch = "+", cex = 0.5)
abline(h = seq(-3, 3, 1), lty = 2, col = 2)
abline(h = 0, lty = 1)

# r.pgInterP
# Inter1
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(2, 1))
xmax <- ui[1]
ymax <- ncn[1]
nband1 <- nband1 / nn1
nband2 <- nband2 / nn1
plot(x = ui, y = ncn, type = "l", xlim = c(0, xmax), ylim = c(0, ymax),
xaxs = "i", yaxs = "i", main = "Interval-Length Distribution",
xlab = "U(i) = exp{-(Normalized Interval Length)}",
ylab = "Normalized Cumulative Number")
points(c(0, xmax), c(0, ymax)+nband1, type = "l", lty=2, col = 2)
points(c(0, xmax), c(0, ymax)+nband2, type = "l", lty=2, col = 2)
points(c(0, xmax), c(0, ymax)-nband1, type = "l", lty=2, col = 2)
points(c(0, xmax), c(0, ymax)-nband2, type = "l", lty=2, col = 2)
points(c(0, xmax), c(0, ymax), type = "l", lty=2, col = 2)

# Inter2
data<- matrix(, ncol = 2, nrow = nn1-1)
data[,1] <- sui[1:(nn1-1)]
data[,2] <- sui[2:nn1]
plot(data, type = "p", pch = "+", xaxs = "i", yaxs = "i",
main="Successive Pair of Intervals", xlab = "U(i)", ylab = "U(i+1)")

# r.pgPalm
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(1, 1))
band <- xp
data <- matrix(, ncol = 2, nrow = npoint)
data[, 1] <- c(1:npoint) * dlt
data[, 2] <- xrate
plot(data, type = "p", pch = "*", main = "Palm Intensity",
xlab = "Elapsed Time After Each Event", ylab = "Occurrence Rate")
points(data, type = "l")
abline(h = band, lty = 2)
abline(h = mean(band))

# r.pgVTC
if (separate.graphics == TRUE)
dev.new()
plot(x = xtime, y = sigma, type = "p", pch = "+", cex = 0.5,
ylim = range(pretty(c(sigma, yvar))), main = "Variance-Time Curve",
xlab = "Time", ylab = "Var{N(0,Time)}")
points(xtime, yvar[1, ], type = "l", col = 2)
for (i in 2:5)
points(xtime, yvar[i, ], type = "l", lty=2, col = 2)
abline(h = 0)
abline(v = 0)
}

ptspec <- function(data, nfre, prdmin, prd, nsmooth = 1, pprd, interval,
plot = TRUE)
{
# data                #  data of events
n <- length(data)   #  number of events of data set
# nfre                #  number of sampling frequencies of spectra
nh1 <- nfre+1
# prdmin              #  the minimum periodicity of the sampling
# prd                 #  a periodicity for calculating the rayleigh probability
# smooth              #  number for smoothing of periodgram
is <- nsmooth
# pprd                #  particular periodcities to be investigated among others
nt <- length(pprd)  #  number of particular periodcities to be investigated
# interval            #  length of observed time interval of events

z <- .Call("PtspecC",
as.double(data),
as.integer(n),
as.double(interval),
as.double(pprd),
as.double(prdmin),
as.double(prd),
as.integer(nfre),
as.integer(nt),
as.integer(is))

wt <- z[[6L]]
ht <- z[[7L]]
w <- z[[8L]]
h <- z[[9L]]
g <- z[[10L]]

if (plot == TRUE) {
par(mfrow = c(2, 1))
pt.spec <- matrix(, nrow = nh1-1, ncol = 2)
pt.spec[, 1] <- w[2:nh1]
pt.spec[, 2] <- h[2:nh1]
ymin <- min(0.0, pt.spec[, 2])

ram <- nh1 / interval
sgm <- log(interval / prdmin / pi) - log(0.05)
sgm1 <- log10(ram * sgm) * 10.e0
sgm <- log(interval /prdmin / pi) - log(0.1)
sgm2 <- log10(ram * sgm) * 10.e0
sgm <- log(interval / prdmin / pi) - log(0.01)
sgm3 <- log10(ram * sgm) * 10.e0

confidence.bar <- c(sgm1, sgm2, sgm3)
vertical.bar <- matrix(, nrow = nt, ncol = 2)
vertical.bar[, 1] <- wt / (2 * pi)
vertical.bar[, 2] <- ht
pt.spec[, 1] <- pt.spec[, 1] / (2 * pi)
plot(pt.spec, pch = 1, cex = 0.5, xlab = "Frequency", ylab = "DB")
points(vertical.bar, pch = "+", cex = 2, col = 2)
points(vertical.bar, lty = 2, type = "h")
abline(h=confidence.bar, lty = 3)

plot(pt.spec[, 1], 10^(pt.spec[, 2] / 10), pch = 1, cex = 0.5,
xlab = "Frequency", ylab = "")
points(vertical.bar[, 1], 10^(vertical.bar[, 2] / 10), pch = "+", cex = 2,
col = 2)
points(vertical.bar[, 1], 10^(vertical.bar[, 2] / 10), lty = 2, type = "h")
abline(h = 10^(confidence.bar / 10), lty = 3)
par(mfrow = c(1, 1))
}

ptspec.out <- list(f = w, db = h, power = g, rayleigh.prob=z[[1L]],
distance = z[[2L]], rwx = z[[3L]], rwy = z[[4]],
phase = z[[5L]], n = n, nfre = nfre, prdmin = prdmin,
nsmooth = nsmooth, interval = interval)

class(ptspec.out) <- "ptspec"
return(ptspec.out)
}

print.ptspec <- function(x, ...)
{
om <- 2 * pi / x\$prdmin
message(gettextf("\n maximum frequency= %e divided into %d points",
om, x\$nfre), domain = NA)

ave <- x\$n / x\$interval
message(gettextf(" average= %e\n", ave), domain = NA)

message(gettextf(" rayleigh probability = %f", x\$rayleigh.prob),
domain = NA)
message(gettextf(" distance = %f", x\$distance), domain = NA)
message(gettextf(" rwx = %f\trwy = %f", x\$rwx, x\$rwy), domain = NA)
message(gettextf(" phase = %f", x\$phase), domain = NA)

nh1 <- length(x\$f)
message("\n frequency\tD.B.\t\tpower", domain = NA)
for(i in 2:nh1)
message(gettextf(" %f\t%f\t%f", x\$f[i], x\$db[i], x\$power[i]), domain = NA)

if (x\$nsmooth == 1) {
xi <- x\$f / (2 * pi)
lt4 <- 4 / x\$interval
message("\n\n list of high powers(>4) with period < t/4", domain = NA)
message(" frequency\tperiods\tpowers", domain = NA)
for (i in 1:nh1)
if (xi[i] >= lt4 && x\$power[i] >= 4) {
prd4 <- 1 / xi[i]
message(gettextf(" %f\t%f\t%f", xi[i],prd4,x\$power[i]), domain = NA)
}
message("\n", domain = NA)
}
}

linsim <- function(data, interval, c, d, ax, ay, at, ptmax)
{
kxx <- length(ax)      # kxx-1 is the order of lgp trf; xx --> xx
kxy <- length(ay)      # kxy-1 is the oeder of lgp trf; yy --> xx
kxz <- length(at)      # kxz-1 is the order of polynomial for xx data
t <- interval          # length of time interval in which events take place
# c,d                    # exponential coefficients in lgp corresponding to xx
# and xy, respectively
# ax,ay,at               # coefficients of the corresponding polynomials
mm <- length(data)
yy <- c(data, t)
# ptmax                  # ptxmax: an upper bound of trend polynomial
kmax <- max(kxx, kxy, kxz)
kmax <- max(kmax, 3)

z <- .Call("LinsimC",
as.integer(kxx),
as.integer(kxy),
as.integer(kxz),
as.double(t),
as.double(c),
as.double(d),
as.double(ax),
as.double(ay),
as.double(at),
as.double(yy),
as.integer(mm),
as.double(ptmax),
as.integer(kmax))

ier <- z[[5L]]
if (ier != 0)
stop("Simulated data length is greater than 2*(original data length)")
err <- z[[4L]]
if (err != 0.)
stop(gettextf("'ptmax' is incorrect. prob = %f", err), domain = NA)

simul <- z[[1L]][1:z[[2L]]]
input <- data[1:z[[3L]]]

linsim.out <- list(in.data=input,sim.data=simul )
return(linsim.out)
}

linlin <- function(external, self.excit, interval, c, d, ax = NULL, ay = NULL,
ac = NULL, at = NULL, opt = 0, tmpfile = NULL, nlmax = 1000)
{
yy <- external       # input data set
xx <- self.excit     # self-exciting data set
t <- interval        # length of time interval in which events take place

x <- c(c,d)
nn <- length(xx)
mm <- length(yy)
if (is.null(ax)) {
kkx <- 0
ax <- 0.0
}  else {
x <- c(x, ax)
kkx <- length(ax)  # kxx-1 is the order of lgp trf; xx --> xx
}
if (is.null(ay)) {
kky <- 0
ay <- 0.0
} else {
x <- c(x, ay)
kky <- length(ay)  # kxy-1 is the oeder of lgp trf; yy --> xx
}
if (is.null(ac)) {
kkc <- 0
ac <- 0.0
} else {
x <- c(x, ac)
kkc <- length(ac)  # kxz-1 is the order of polynomial for xx data
}
if (is.null(at)) {
kkt <- 0
at <- 0.0
} else {
x <- c(x, at)
kkt <- length(at)  # kxz-1 is the order of polynomial for xx data
}
n <- length(x)

# opt        # =0 : minimize the likelihood with fixed exponential coefficient c
# =1 :  not fixed d

kmax <- max(kkx, kky) + 1
kmax <- max(kmax, 3)
nlm <- nlmax
if (is.null(tmpfile))
nlm <- 0

z <- .Call("LinlinC",
as.integer(n),
as.double(x),
as.integer(opt),
as.double(t),
as.integer(nn),
as.integer(mm),
as.double(xx),
as.double(yy),
as.integer(kkx),
as.integer(kky),
as.integer(kmax),
as.integer(kkc),
as.integer(kkt),
as.integer(nlm) )

ier <- z[[16L]]
if (ier == -1)
stop(" subroutine funct : n or kkx or kky kkc kkt error")

# initial estimates
x1 <- z[[1L]]
c1 <- x1[1]
d1 <- x1[2]
ax1 <- NULL
ay1 <- NULL
ac1 <- NULL
at1 <- NULL
if (kkx > 0)
ax1 <- x1[3:(kkx+2)]
if (kky > 0)
ay1 <- x1[(kkx+3):(kkx+kky+2)]
if (kkc > 0)
ac1 <- x1[(kkx+kky+3):(kkx+kky+kkc+2)]
if (kkt > 0)
at1 <- x1[(kkx+kky+kkc+3):n]

# final estimates
x2 <- z[[2L]]
c2 <- x2[1]
d2 <- x2[2]
ax2 <- NULL
ay2 <- NULL
ac2 <- NULL
at2 <- NULL
if (kkx > 0)
ax2 <- x2[3:(kkx+2)]
if (kky > 0)
ay2 <- x2[(kkx+3):(kkx+kky+2)]
if (kkc > 0)
ac2 <- x2[(kkx+kky+3):(kkx+kky+kkc+2)]
if (kkt > 0)
at2 <- x2[(kkx+kky+kkc+3):n]

if (nlm > 0) {
nfunct <- 0
nl <- z[[15L]]
x <- array(z[[10L]], c(n, 5))
g <-array(z[[11L]], c(n, 5))
id <- z[[12L]][1:nl]
ramda <- z[[13L]][1:nl]
ee <- z[[14L]][1:nl]
print.process(nfunct, n, x, g, id, ramda, ee, tmpfile)
}

linlin.out <- list(c1 = c1, d1 = d1, ax1 = ax1, ay1 = ay1, ac1 = ac1,
at1 = at1, c2 = c2, d2 = d2, ax2 = ax2, ay2 = ay2,
ac2 = ac2, at2=at2, aic2 = z[[3L]], ngmle = z[[4L]],
rayleigh.prob = z[[5L]], distance = z[[6L]], rwx = z[[7L]],
rwy = z[[8L]], phase = z[[9L]])

class(linlin.out) <- "linlin"
return(linlin.out)
}

print.linlin <- function(x, ...)
{
message(gettextf("\n rayleigh probability = %f", x\$rayleigh.prob), domain = NA)
message(gettextf(" distance = %f", x\$distance), domain = NA)
message(gettextf(" rwx = %f\t\trwy = %f", x\$rwx, x\$rwy), domain = NA)
message(gettextf(" phase = %f", x\$phase), domain = NA)

kkx <- length(x\$ax1)
kky <- length(x\$ay1)
kkc <- length(x\$ac1)
kkt <- length(x\$at1)

message("\n initial estimates", domain = NA)
message(gettextf(" c, d\t%f\t%f", x\$c1, x\$d1), domain = NA)
if (kkx > 0) {
pmsg <- " ax"
for (i in 1:kkx)
pmsg <- paste(pmsg, gettextf("\t%e", x\$ax1[i]))
message(pmsg, domain = NA)
}
if (kky > 0) {
pmsg <- " ay"
for (i in 1:kky)
pmsg <- paste(pmsg, gettextf("\t%e", x\$ay1[i]))
message(pmsg, domain = NA)
}
if (kkc > 0) {
pmsg <- " ac"
for (i in 1:kkc)
pmsg <- paste(pmsg, gettextf("\t%e", x\$ac1[i]))
message(pmsg, domain = NA)
}
if (kkt > 0) {
pmsg <- " at"
for (i in 1:kkt)
pmsg <- paste(pmsg, gettextf("\t%e", x\$at1[i]))
message(pmsg, domain = NA)
}

message("\n final outputs")
message(gettextf(" c, d\t%f\t%f", x\$c2, x\$d2), domain = NA)
if (kkx > 0) {
pmsg <- " ax"
for (i in 1:kkx)
pmsg <- paste(pmsg, gettextf("\t%e", x\$ax2[i]))
message(pmsg, domain = NA)
}
if (kky > 0) {
pmsg <- " ay"
for (i in 1:kky)
pmsg <- paste(pmsg, gettextf("\t%e", x\$ay2[i]))
message(pmsg, domain = NA)
}
if (kkc > 0) {
pmsg <- " ac"
for (i in 1:kkc)
pmsg <- paste(pmsg, gettextf("\t%e", x\$ac2[i]))
message(pmsg, domain = NA)
}
if (kkt > 0) {
pmsg <- " at"
for (i in 1:kkt)
pmsg <- paste(pmsg, gettextf("\t%e", x\$at2[i]))
message(pmsg, domain = NA)
}

message(gettextf("\n negative max likelihood = %f", x\$ngmle), domain = NA)
message(gettextf(" AIC/2 = %f\n", x\$aic2), domain = NA)

}

simbvh <- function(interval, axx = NULL, axy = NULL, axz = NULL, ayx = NULL,
ayy = NULL, ayz = NULL, c, d, c2, d2, ptxmax, ptymax)
{
t <- interval          # length of time interval in which events take place

# axx(i),axy(i),ayx(i),ayy(i),axz(i),ayz(i)
#    coefficients of the corresponding polynomials

if (is.null(axx)) {
axx <- 0
kxx <- 0
} else {
kxx <- length(axx)     # kxx-1 is the order of lgp trf; xx --> xx
}
if (is.null(axy)) {
axy <- 0
kxy <- 0
} else {
kxy <- length(axy)     # kxy-1 is the oeder of lgp trf; yy --> xx
}
if (is.null(ayx)) {
ayx <- 0
kyx <- 0
} else {
kyx <- length(ayx)     # kyx-1 is the oeder of lgp trf; xx --> yy
}
if (is.null(ayy)) {
ayy <- 0
kyy <- 0
} else {
kyy <- length(ayy)     # kyy-1 is the oeder of lgp trf; yy --> yy
}
if (is.null(axz)) {
axz <- 0
kxz <- 0
} else {
kxz <- length(axz)     # kxz-1 is the order of polynomial for xx data
}
if (is.null(ayz)) {
ayz <- 0
kyz <- 0
} else {
kyz <- length(ayz)     # kyz-1 is the order of polynomial for yy data
}

# c,d,c2,d2              # exponential coefficients in lgp corresponding to
# xx,xy,yx and yy, respectively

# ptxmax                 # upper bounds of trend polynomials corresponding to xz
# ptymax                 # upper bounds of trend polynomials corresponding to yz
kmax <- max(kxx, kxy, kyx, kyy)
kmax <- max(kmax, 3)

nnmax <- 10000
mmmax <- 10000

z <- .Call("SimbvhC",
as.integer(kxx),
as.integer(kxy),
as.integer(kxz),
as.integer(kyx),
as.integer(kyy),
as.integer(kyz),
as.double(t),
as.double(c),
as.double(d),
as.double(c2),
as.double(d2),
as.double(axx),
as.double(axy),
as.double(axz),
as.double(ayx),
as.double(ayy),
as.double(ayz),
as.double(ptxmax),
as.double(ptymax),
as.integer(kmax),
as.integer(nnmax),
as.integer(mmmax))

ier <- z[[6L]]
if (ier != 0)
stop("Simulated data length is greater than 10000.")
err <- z[[5L]]
if (err != 0.)
warning(gettextf("Are ptxmax & ptymax correct? prob=%f", err), domain = NA)

nx <- z[[3L]]
ny <- z[[4L]]
simbvh.out <- list(x=z[[1L]][1:nx], y=z[[2L]][1:ny] )
return(simbvh.out)
}

momori <- function(data, mag = NULL, threshold = 0.0, tstart, tend, parami,
tmpfile = NULL, nlmax = 1000)
{

# data                   #  point process data
n <- length(data)      #  total number of variable
# mag                    #  magnitude
# threshold              #  threshold magnitude
# tstart                 #  the start of the target perio
zts = tstart
zte=tend
# tend                   #  the end of the target period
# parami                 #  initial estimates of parameters
np <- length(parami)
np1 <- np + 1
if (np1 != 5)
stop("Number of parameters is worng.")
parami <- c(parami[1:3], 0, parami[4])

nn <- 0
xx <- NULL
magnitude <- NULL

if (is.null(mag)) {
xx <- data
nn <- n
} else {
for (i in 1:n)
if (mag[i] >= threshold)
if ((data[i] >= zts) && (data[i] <= zte)) {
nn <- nn + 1
xx <- c(xx, data[i])
magnitude <- c(magnitude, mag[i])
}
}

nfunct <- 6
ncount <- 1
nlm <- nlmax
if (is.null(tmpfile))
nlm <- 0

z <- .Call("MomoriC",
as.double(xx),
as.integer(nn),
as.double(parami),
as.integer(np),
as.double(zts),
as.double(zte),
as.integer(ncount),
as.integer(nfunct),
as.integer(nlm) )

#  param <- c(t=z\$x[1], k=z\$x[2], c=z\$x[3], p=x0, cls=z\$x[4])
pa1 <- z[[4L]]
pa2 <- list(t_i = z[[6L]], K = z[[8L]], c = z[[9L]], p = z[[10L]],
cls = z[[11L]])

if (nlm > 0) {
x <- array(z[[2L]], c(np, 2))
g <- array(z[[3L]], c(np, 2))
nl <- z[[17L]]
id <- z[[12L]][1:nl]
ramda <- z[[13L]][1:nl]
xx0 <- array(z[[14L]], c(np, nlmax))
xx1 <- xx0[1:np, 1:nl]
h <- array(z[[15L]], c(np, np, 2))
hf <- array(z[[16L]], c(np, np, 2, 2))
print.process6(id, x, g, h, hf, ramda, xx1, tmpfile)
}

momori.out <- list(param = pa1, ngmle = z[[1L]], aic = z[[5L]], plist = pa2)
return(momori.out)
}

respoi <- function(time, mag, param, zts, tstart, zte, threshold = 0.0,
plot = TRUE)
{
# time                   #  time from the main shock in days
nd <- length(time)
# mag                    #  magnitude

# param                  #  estimates of parameters
np <- length(param) + 1
if (np != 5)
stop("Number of parameters is worng.")
param <- c(param[1:3],0,param[4])

# zts                    #  the start of the precursory period
# zstart                 #  the start of the target period
# zte                    #  the end of the target period

dep <- rep(0, nd)
xp <- rep(0, nd)
yp <- rep(0, nd)

z <- .Call("RespoiC",
as.double(time),
as.double(mag),
as.double(dep),
as.double(xp),
as.double(yp),
as.integer(nd),
as.double(param),
as.integer(np),
as.double(zts),
as.double(zte),
as.double(tstart),
as.double(threshold))

n <- z[[8L]]
cn <- c(1:n) - z[[5L]]
ti <- z[[7L]][1:n]
mag1 <- z[[1L]][1:n]

if (plot == TRUE) {
t <- threshold
level <- min(0, cn[1])
xrange <- c(min(ti), max(ti))
mgrange <- max(cn) / 4
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(max(cn), max(ti)))
plot(xrange, yrange, type = "n", main = "Omori-Utsu Residual",
xlab = "Transformed Time", ylab = "Cumulative Number of Events",
lwd = 1, pty = "s", xaxs = 'r')
points(ti, 1:length(ti)+cn[1]+1, type = 's')
mgmax <- max(mag1 - t + 1)
mag1 <- mag1 - t + 0.5
segments(ti, bottom, ti, mag1/mgmax*mgrange+bottom)
abline(h = bottom)
abline(h = 0)
abline(v = 0, lty = 2)
abline(0,1, lty = 1, col = 'red')
s <- tstart
}

respoi.out <- list(trans.time = ti, cnum = cn)
return(respoi.out)
}

respoi2 <- function(etas, param, zts, tstart, zte, threshold = 0.0, plot = TRUE)
{
if (dim(etas)[2] != 9)
stop("etas is a sequential data with 9 columns-format.")

# param                  #  estimates of parameters
np <- length(param) + 1
if (np != 5)
stop("Number of parameters is worng.")
eparam <- c(param[1:3], 0, param[4])

nd <- dim(etas)[1]
#  xp <- etas[,2]         #  longitude
#  yp <- etas[,3]         #  latitude
#  mag <- etas[,4]        #  magnitude
#  time <- etas[,5]       #  time from the main shock in days
#  dep <- etas[,6]        #  depth
xp <- rep(0, nd)
yp <- rep(0, nd)
mag <- rep(0, nd)
time <- rep(0, nd)
dep <- rep(0, nd)
for (i in 1: nd) {
xp[i] <- etas[i, 2]
yp[i] <- etas[i, 3]
mag[i] <- etas[i, 4]
time[i] <- etas[i, 5]
dep[i] <- etas[i, 6]
}

z <- .Call("RespoiC",
as.double(time),
as.double(mag),
as.double(dep),
as.double(xp),
as.double(yp),
as.integer(nd),
as.double(eparam),
as.integer(np),
as.double(zts),
as.double(zte),
as.double(tstart),
as.double(threshold))

n <- z[[8L]]
cn <- c(1:n) - z[[5L]]
ti <- z[[7L]][1:n]
mag1 <- z[[1L]][1:n]

if (plot == TRUE) {
t <- threshold
level <- min(0, cn[1])
xrange <- c(min(ti), max(ti))
mgrange <- max(cn) / 4
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(max(cn), max(ti)))
plot(xrange, yrange, type = "n", main = "Omori-Utsu Residual",
xlab = "Transformed Time", ylab = "Cumulative Number of Events",
lwd = 1, pty = "s", xaxs = 'r')
points(ti, 1:length(ti)+cn[1]+1, type = 's')
mgmax <- max(mag1 - t + 1)
mag1 <- mag1 - t + 0.5
segments(ti, bottom, ti, mag1/mgmax*mgrange+bottom)
abline(h = bottom)
abline(h = 0)
abline(v = 0, lty = 2)
abline(0, 1, lty = 1, col = 'red')
s <- tstart
}

res <- matrix(, ncol = 7, nrow = n)
res[,1] <- cn
res[,2] <- z[[3L]][1:n]
res[,3] <- z[[4L]][1:n]
res[,4] <- z[[1L]][1:n]
res[,5] <- z[[6L]][1:n]
res[,6] <- z[[2L]][1:n]
res[,7] <- ti
res <- data.frame(res)
names(res) <- c("no.", "longitude", "latitude", "magnitude", "time", "depth",
"trans.time")

respoi2.out <- list(resData = res)
return(respoi2.out)

}

etasap <- function(time, mag, threshold = 0.0, reference = 0.0, parami,
zts = 0.0, tstart, zte, approx = 2, tmpfile = NULL,
nlmax = 1000, plot = TRUE)
{
# time                   #  time from the main shock in days
n <- length(time)
# mag                    #  magnitude
# threshold              #  threshold magnitude
# reference              #  reference magnitude
# parami                 #  initial estimates of parameters
np <- length(parami)
if (np != 5)
stop("Number of parameters is worng.")
# zts                    #  the start of the precursory period
# tstart                 #  the start of the target period
# zte                    #  the end of the target period
# approx                 #  the level for approximation version (1, 2, 4, 8, 16)
#  =0 : the exact version

nn <- 0
xx <- NULL
magnitude <- NULL
amx1 <- threshold
amx2 <- 10.0

for (i in 1:n )
if (mag[i] >= amx1 && mag[i] <= amx2)
if (time[i] >= zts && time[i] <= zte) {
nn <- nn + 1
xx <- c(xx, time[i])
magnitude <- c(magnitude, mag[i])
}

nfunct <- 9
if (approx == 0)
nfunct <- 4
nlm <- nlmax
if (is.null(tmpfile))
nlm <- 0

z <- .Call("EtasapC",
as.double(xx),
as.double(magnitude),
as.integer(nn),
as.double(reference),
as.double(threshold),
as.double(parami),
as.integer(np),
as.double(zts),
as.double(zte),
as.double(tstart),
as.integer(nfunct),
as.integer(approx),
as.integer(nlm) )

x <- z[[2L]]
pa <- list(mu = x[1], K = x[2], c = x[3], alpha = x[4], p = x[5])

if (nlm > 0) {
nl <- z[[8L]]
id <- z[[5L]][1:nl]
g <- z[[3L]]
ee <- z[[6L]][1:nl]
x0 <- array(z[[7L]], c(np, nlmax))
x1 <- x0[1:np, 1:nl]
print.process9(id, x, g, ee, x1, tmpfile)
}

if (plot == TRUE) {
mag1 <- min(magnitude)
ti <- xx
zero <- rep(0, nn)
cn <- 1:nn
xrange <- c(min(xx), max(xx))
mgrange <- max(cn)/4
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(cn))
plot(xrange, yrange, type = "n", main = "Seismicity CUM & M-T plot",
xlab = "Ordinary Time (Days)", ylab = "Cumulative Number of Events",
lwd = 1, pty = "s", xaxs = 'r', yaxs = 'i')
lines(ti, cn, type = 'S')
mgmax <- max(magnitude - mag1 + 1)
magnitude1 <- magnitude - mag1 + 0.5
segments(xx, bottom, xx, magnitude1/mgmax*mgrange+bottom)
abline(h = 0)
abline(h = bottom)
}

etasap.out <- list(ngmle = z[[1L]], aic2 = z[[4L]], param = pa)
return(etasap.out)
}

etarpp <- function(time, mag, threshold = 0.0, reference = 0.0, parami,
zts = 0.0, tstart, zte, ztend = NULL, plot = TRUE)
{
# time                   #  time from the main shock in days
n <- length(time)
# mag                    #  magnitude
# threshold              #  threshold magnitude
# reference              #  reference magnitude
# parami                 #  initial estimates of parameters
np <- length(parami)
if (np != 5)
stop("Number of parameters is worng.")
# zts                    #  the start of the precursory period
# tstart                 #  the start of the target period
# zte                    #  the end of the target period
# ztend                  #  the end of the predictionperiod
if (is.null(ztend))
ztend <- time[n]

mm <- 0
nn <- 0
time0 <- NULL
mag0 <- NULL
time1 <- NULL
mag1 <- NULL
for (i in 1:n)
if (mag[i] >= threshold) {
mm <- mm + 1
time0 <- c(time0, time[i])
mag0 <- c(mag0, mag[i])
if (time[i] >= zts && time[i] <= ztend) {
nn <- nn + 1
time1 <- c(time1, time[i])
mag1 <- c(mag1, mag[i])
}
}

z <- .Call("EtarppC",
as.double(time1),
as.double(mag1),
as.double(reference),
as.integer(nn),
as.double(parami),
as.integer(np),
as.double(zts),
as.double(ztend),
as.double(tstart) )

x <- z[[1L]]
ntstar <- z[[2L]]
cnum <- 1:nn - ntstar

if (plot == TRUE) {
# r.seisetas
## scan("work.etas")
mgmin <- min(mag0)
zero <- rep(0, mm)
cn <- 1:mm
cna <- append(cn, 0, after = 0)
cn1 <- cna[1:mm]
tia <- append(time0, 0, after = 0)
ti1 <- tia[1:mm]
xrange <- c(min(time0), ztend)
mgrange <- max(cn) / 4
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(max(cn), max(x-ntstar)))
mtitle <- paste("ETAS Fit and Prediction\nM>=", threshold, "S=", zts, "T=",
zte, "Tend=", ztend)
plot(xrange, yrange, type = "n", main = mtitle,
xlab = "Ordinary Time (Days)", ylab = "Cumulative Number of Events",
lwd = 1, xlim = c(xrange), pty = "s", xaxs = "r", yaxs = "r")
## scan("work.res")
lines(time1, x + ntstar, type = 'l', col = 'red')
segments(ti1, cn1, time0, cn1)
segments(time0, cn1, time0, cn)
mgmax <- max(mag0 - mgmin + 1)
mag2 <- mag0 - mgmin + 0.5
segments(time0, bottom, time0, mag2/mgmax*mgrange+bottom)
abline(h = 0)
abline(h = bottom)
mark0 <- tstart; abline(v = mark0, lty = 2)
#   mark1 <- ngmle ; abline(v = mark1, lty = 2)
mark2 <- ztend;  abline(v = mark2, lty = 2)
abline(v = tstart, lty = 2)
abline(v = zte, lty = 2)
#    text(max(time0)*0.3, max(cn)*0.95, paste('M>=',mgmin,'S=',zts,'T=',zte,'Tend=',ztend))

# r.retas
## scan("work.res")
#    X11()
level <- min(0, cnum[1])
cn <- 1:mm + cnum[1]
ti <- x
xrange <- c(min(ti), max(ti))
mgrange <- max(cn / 4)
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(max(cn), max(ti)))
mtitle <- paste("ETAS Residual\nM>=", threshold, "S=", zts, "T=", zte,
"Tend=", ztend)
plot(xrange, yrange, type = "n", main = mtitle, xlab = "Transformed Time",
ylab = "Cumulative Number of Events", lwd = 1, pty = "s", xaxs = "r")
points(x, 1:nn+cnum[1]+1, type = 's')
mgmax <- max(mag1 - threshold + 1)
mag3 <- mag1 - threshold + 0.5
segments(ti, bottom, ti, mag3/mgmax*mgrange+bottom)
abline(h = bottom)
abline(h = 0)
abline(v = 0, lty = 2)
abline(0, 1, lty = 1, col = 'red')
timax <- max(ti[time1 <= zte])
abline(v = timax, lty = 2)
#    text(max(ti)*0.4, max(cn)*0.9, paste('M>=',threshold,'S=',zts,'T=',zte,'Tend=',ztend))
}

etarpp.out <- list(trans.time = x, no.tstart=ntstar)
return(etarpp.out)
}

etarpp2 <- function(etas, threshold = 0.0, reference = 0.0, parami, zts= 0.0,
tstart, zte, ztend = NULL, plot = TRUE)
{
n <- dim(etas)[1]
# threshold              #  threshold magnitude
# reference              #  reference magnitude
# parami                 #  initial estimates of parameters
np <- length(parami)
if (np != 5)
stop("Number of parameters is worng.")
# zts                    #  the start of the precursory period
# tstart                 #  the start of the target period
# zte                    #  the end of the target period
# ztend                  #  the end of the predictionperiod
if (is.null(ztend))
ztend <- time[n]

xp <- etas[, 2]        #  longitude
yp <- etas[, 3]        #  latitude
mag <- etas[, 4]       #  magnitude
time <- etas[, 5]      #  time from the main shock in days
dep <- etas[, 6]       #  depth

mm <- 0
nn <- 0
time0 <- NULL
mag0 <- NULL
xp1 <- NULL
yp1 <- NULL
mag1 <- NULL
time1 <- NULL
dep1 <- NULL
for (i in 1:n )
if (mag[i] >= threshold) {
mm <- mm + 1
time0 <- c(time0, time[i])
mag0 <- c(mag0, mag[i])
if (time[i] >= zts && time[i] <= ztend) {
nn <- nn + 1
xp1 <- c(xp1, xp[i])
yp1 <- c(yp1, yp[i])
mag1 <- c(mag1, mag[i])
time1 <- c(time1, time[i])
dep1 <- c(dep1, dep[i])
}
}

x <- rep(0, nn)
ntstar <- 0

z <- .Call("EtarppC",
as.double(time1),
as.double(mag1),
as.double(reference),
as.integer(nn),
as.double(parami),
as.integer(np),
as.double(zts),
as.double(ztend),
as.double(tstart))

x <- z[[1L]]
ntstar <- z[[2L]]
cnum <- 1:nn - ntstar

if (plot == TRUE) {
# r.seisetas
## scan("work.etas")
mgmin <- min(mag0)
zero <- rep(0, mm)
cn <- 1:mm
cna <- append(cn, 0, after = 0)
cn1 <- cna[1:mm]
tia <- append(time0, 0, after = 0)
ti1 <- tia[1:mm]
par(pty = "s", xaxs = "r", yaxs = "r")
xrange <- c(min(time0), ztend)
mgrange <- max(cn) / 4
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(max(cn), max(x-ntstar)))
mtitle <- paste("ETAS Fit and Prediction\nM>=", threshold, "S=", zts, "T=",
zte, "Tend=", ztend)
plot(xrange, yrange, type = "n", main = mtitle,
xlab = "Ordinary Time (Days)",ylab = "Cumulative Number of Events",
lwd = 1, xlim = c(xrange))
## scan("work.res")
lines(time1, x+ntstar, type = 'l', col = 'red')
segments(ti1, cn1, time0, cn1)
segments(time0, cn1, time0, cn)
mgmax <- max(mag0 - mgmin + 1)
mag2 <- mag0 - mgmin + 0.5
segments(time0, bottom, time0, mag2/mgmax*mgrange+bottom)
abline(h = 0)
abline(h = bottom)
mark0 <- tstart; abline(v = mark0, lty = 2)
#   mark1 <- ngmle ; abline(v = mark1, lty = 2)
mark2 <- ztend;  abline(v = mark2, lty = 2)
abline(v = tstart, lty = 2)
abline(v = zte, lty = 2)
#    text(max(time0)*0.3, max(cn)*0.95, paste('M>=',mgmin,'S=',zts,'T=',zte,'Tend=',ztend))

# r.retas
## scan("work.res")
#    X11()
level <- min(0, cnum[1])
cn <- 1:mm+cnum[1]
ti <- x
xrange <- c(min(ti), max(ti))
mgrange <- max(cn / 4)
bottom <- min(cn) - mgrange
yrange <- c(bottom, max(max(cn), max(ti)))
par(pty = "s", xaxs = "r")
mtitle <- paste("ETAS Residual\nM>=", threshold, "S=", zts, "T=", zte,
"Tend=", ztend)
plot(xrange, yrange, type = "n", main = mtitle, xlab = "Transformed Time",
ylab = "Cumulative Number of Events", lwd = 1)
points(x, 1:nn+cnum[1]+1, type = 's')
mgmax <- max(mag1 - threshold + 1)
mag3 <- mag1 - threshold + 0.5
segments(ti, bottom, ti, mag3/mgmax*mgrange+bottom)
abline(h = bottom)
abline(h = 0)
abline(v = 0, lty = 2)
abline(0, 1, lty = 1, col = 'red')
timax <- max(ti[time1 <= zte])
abline(v = timax, lty = 2)
#    text(max(ti)*0.4,max(cn)*0.9,paste('M>=',threshold,'S=',zts,'T=',zte,'Tend=',ztend))
}

res <- matrix(, ncol = 7, nrow = nn)
res[, 1] <- cnum
res[, 2] <- xp1
res[, 3] <- yp1
res[, 4] <- mag1
res[, 5] <- time1
res[, 6] <- dep1
res[, 7] <- x
res <- data.frame(res)
names(res) <- c("no.", "longitude", "latitude", "magnitude", "time", "depth",
"trans.time")

etarpp2.out <- list(resData = res)
return(etarpp2.out)
}

etasim1 <- function(bvalue, nd, threshold = 0.0, reference = 0.0, param)
{
# bvlalue                #  b-value of G-R law
# nd                     #  the number of the simulated events
# threshold              #  threshold magnitude
# reference              #  reference magnitude
# param                  #  initial estimates of parameters
np <- length(param)
if (np != 5)
stop("Number of parameters is worng.")
a <- param[1]
b <- param[2]
c <- param[3]
d <- param[4]
p <- param[5]

ic <- 0
tstart <- 0
xm <- rep(0, nd)
zz <- rep(0, nd)

z <- .Call("EtasimC",
as.integer(ic),
as.double(bvalue),
as.double(tstart),
as.integer(nd),
as.double(threshold),
as.double(reference),
as.double(a),
as.double(b),
as.double(c),
as.double(d),
as.double(p),
as.double(xm),
as.double(zz))

sim <- matrix(, ncol = 9, nrow = nd)
sim[,1] <- 1:nd
sim[,2] <- rep(0, nd)
sim[,3] <- rep(0, nd)
sim[,4] <- z[[1L]]
sim[,5] <- z[[2L]]
sim[,6] <- rep(0, nd)
sim[,7] <- rep(0, nd)
sim[,8] <- rep(0, nd)
sim[,9] <- rep(0, nd)
sim <- data.frame(sim)
names(sim) <- c("no.", "longitude", "latitude", "magnitude", "time", "depth",
"year", "month", "days")

return(etasim = sim)
}

etasim2 <- function(etas, tstart, threshold = 0.0, reference = 0.0, param)
{
n <- dim(etas)[1]
mag <- etas[, 4]       #  magnitude
time <- etas[, 5]      #  time from the main shock in days
# tstart                #
# threshold             #  threshold magnitude
# reference             #  reference magnitude
# param                 #  initial estimates of parameters
np <- length(param)
if (np != 5)
stop("Number of parameters is worng.")
a <- param[1]
b <- param[2]
c <- param[3]
d <- param[4]
p <- param[5]

nd <- 0
mag1 <- NULL
time1 <- NULL
for (i in 1:n)
if (mag[i] >= threshold) {
nd <- nd + 1
mag1 <- c(mag1, mag[i])
time1 <- c(time1, time[i])
}

ic <- 1
bvalue <- 0
xx <- rep(0, nd)
probx <- 0

z <- .Call("EtasimC",
as.integer(ic),
as.double(bvalue),
as.double(tstart),
as.integer(nd),
as.double(threshold),
as.double(reference),
as.double(a),
as.double(b),
as.double(c),
as.double(d),
as.double(p),
as.double(mag1),
as.double(time1))

if (z[[3L]] > 1)
stop(gettextf("prob = %f > 1", z[[3L]]))

sim <- matrix(, ncol = 9, nrow = nd)
sim[,1] <- 1:nd
sim[,2] <- rep(0, nd)
sim[,3] <- rep(0, nd)
sim[,4] <- mag1
sim[,5] <- z[[2L]]
sim[,6] <- rep(0, nd)
sim[,7] <- rep(0, nd)
sim[,8] <- rep(0, nd)
sim[,9] <- rep(0, nd)
sim <- data.frame(sim)
names(sim) <- c("no.", "longitude", "latitude", "magnitude", "time",
"depth", "year", "month", "days")

return(etasim2 = sim)
}

print.process <- function(nfunct, n, x, g, id, ramda, ee, tmpfile)
{
nl <- length(ramda)
if (nfunct == 0) {
np <- n
ist <- 0
} else if (nfunct == 1) {
np <- 1
ist <- 1
} else if (nfunct == 2) {
np <- 1
ist <- 2
}
j2 <- 1

for (i in 1:nl) {
k <- id[i]
if (k > 0 & k < 7)
write(sprintf("lambda = %e     e%i = %e", ramda[i], k, ee[i]), tmpfile,
append = TRUE);
if (k == 330)
write(sprintf("lambda = %e     log likelihood  = %e", ramda[i], ee[i]),
tmpfile, append = TRUE);
if (k == 340)
write(sprintf("\n                          log-likelihood = %e",
ramda[i]), tmpfile, append = TRUE);
if (k == -1) {
write("\n ----- x -----", tmpfile, append = TRUE)
out <- NULL
for (j in 1:np)
out <- paste(out, sprintf("  %e", x[j,j2]))
write(out, tmpfile, append = TRUE)

write("\n *** gradient ***", tmpfile, append = TRUE)
out <- NULL
for (j in 1:np)
out <- paste(out, sprintf("  %e", g[j,j2]**2))
out <- paste(out, "\n")
write(out, tmpfile, append = TRUE)

np <- np+ist
j2 <- j2+1
}
}
write(sprintf("................................ Process steps  1- %d ................................\n", nl),
tmpfile, append = TRUE)
}

print.process6 <- function(id, x, g, h, hf, ramda, xx1, tmpfile)
{
np <- dim(x)[1]
nl <- length(id)
j2 <- 1

for (i in 1:nl) {
k <- id[i]
if (k == 8) {
out <- sprintf(" -ll = %e", ramda[i])
for (j in 1:np)
out <- paste(out, sprintf("\t%e", xx1[j, i]**2))
write(out, tmpfile, append = TRUE)
}
if (k == 330)
write(sprintf(" lambda = %e     -LL =%e   %e   %e", ramda[i],xx1[1,i],xx1[2,i],xx1[3,i]),
tmpfile, append = TRUE)
if (k == 340)
write(sprintf("\n\tinitial (-1)*Log-Likelihood = %e", xx1[1,i]), tmpfile)
if (k == 600) {
write("\n ----- x -----", tmpfile, append = TRUE)
out <- NULL
for (j in 1:np)
out <- paste(out, sprintf("  %e", x[j,j2]))
write(out, tmpfile, append = TRUE)

write("\n *** gradient ***", tmpfile, append = TRUE)
out <- NULL
for (j in 1:np)
out <- paste(out, sprintf("  %e", g[j,j2]))
write(out, tmpfile, append = TRUE)

write("\n *** estimated inverse hessian gradient ***", tmpfile, append = TRUE)
for (j0 in 1:np) {
out <- NULL
for (j1 in 1:np)
out <- paste(out, sprintf("  %e", h[j0,j1,j2]))
write(out, tmpfile, append = TRUE)
}

write("\n *** fisher matrix ***", tmpfile, append = TRUE)
for (j0 in 1:np)  {
out <- NULL
for (j1 in 1:np)
out <- paste(out, sprintf("  %e", hf[j0,j1,j2,1]))
write(out, tmpfile, append = TRUE)
}

write("\n *** inverse fisher ***", tmpfile, append = TRUE)
for (j0 in 1:np)  {
out <- NULL
for (j1 in 1:np)
out <- paste(out, sprintf("  %e", hf[j0,j1,j2,2]))
out <- paste(out, "\n")
write(out, tmpfile, append = TRUE)
}
j2 <- j2 + 1
}
}
write(sprintf("................................ Process steps  1- %d ................................\n", nl),
tmpfile, append = TRUE)
}

print.process9 <- function(id, x, g, ee, x1, tmpfile)
{
np <- dim(x1)[1]
nl <- length(id)

for (i in 1:nl) {
k <- id[i]
if (k == 8) {
out <- sprintf(" -ll = %e", ee[i])
for (j in 1:np)
out <- paste(out, sprintf("\t%e", x1[j,i]**2))
write(out, tmpfile, append = TRUE)
}
if (k == 330)
write(sprintf(" lambda = %e     -LL =%e   %e   %e", ee[i], x1[1,i],
x1[2,i], x1[3,i]), tmpfile, append = TRUE)
if (k == 340)
write(sprintf("\n\tinitial (-1)*Log-Likelihood = %e", x1[1,i]), tmpfile)
if (k == 600) {
write("\n ----- x -----", tmpfile, append = TRUE)
out <- NULL
for (j in 1:np)
out <- paste(out, sprintf("  %e", x[j]))
write(out, tmpfile, append = TRUE)
write("\n\n *** gradient ***", tmpfile, append = TRUE)
out <- NULL
for (j in 1:np)
out <- paste(out, sprintf("  %e", g[j]**2))
out <- paste(out, "\n")
write(out, tmpfile, append = TRUE)
}
}
write(sprintf("................................ Process steps  1- %d ................................\n", nl),
tmpfile, append = TRUE)
}
```

## Try the SAPP package in your browser

Any scripts or data that you put into this service are public.

SAPP documentation built on July 2, 2020, 2:59 a.m.