Nothing
R/drsub.R
In kzfs: Multi-Scale Motions Separation with Kolmogorov-Zurbenko Periodogram Signals
# --------------------------------------------------------------------------
# @title
# Internal Function Used In Wave Parameter Estimation
#
# @description
# Function \code{commcore} is designed to find the common elements
# in a given list, which is the pool of all potential parameters. The
# common elements should appear in most data groups of paired directional
# periodograms.
#
# \code{kzpdr.rec} will search all variables in R global environment
# for available records of directional periodograms, and return the
# specified part as output.
#
# \code{kzpdr.pairs} returns a list of all available pairs of sampling
# directions and the counts number for the spikes on each direction.
#
# \code{kzpdr.prj} is used for calculating the projected parameters
# based on directional periodograms.
#
# Function \code{ckGV} refers to a given record set of directional
# periodogram and checks duplicated directions. The returned value is a
# list of sampling directions without those duplications.
#
# \code{shade} decides if a spectral signals would sometime disappear
# in the given sampling directions. The return value is a vector of 0/1.
# 0 means the signal should be seen in all given directions. 0 means not.
# 2 and larger may happen when there are very closed sampling angles.
#
# \code{tight} is used to get the minimum difference in a group of
# potential directions or frequencies.
#
# \code{tolerance} is used to calculate the tolerance of the estimations.
#
# Function \code{closure} returns the closure of the nearest neighbors.
#
# \code{distill} excludes the duplicate supports of a cluster. Usually,
# duplications are supports for the same parameters but from the same pairs,
# grouping under the given tolerance condition.
#
# @param pool List for the potential parameters.
# @param STEP The tolerance interval for a data element.
# Defaults to 1.
# @param DIFF The allowed absents times for a suggested
# common element. Default value is 1.
# @param digit Significant digits. Defaults to 0.
# @rdname kzprj
# @export
# @keywords internal
# ---------------------------------------------------------------------------
commcore <- function(pool, STEP=1, DIFF=1, digit=0) {
u <- as.vector(unlist(pool))
flg <- names(unlist(pool))
flg <- gsub(")[0-9]+",")",flg)
a <- round(u, digit)
o <- order(u)
v <- data.frame(v=a,u=u,d=DIFF,flg=flg)
v <- v[o,]
while (TRUE) {
v$upp <- c(999,diff(v[,1]))
v$dwn <- abs(c(rev(diff(rev(v[,1]))),999))
v$x <- ifelse((v$upp == 0 | v$dwn == 0), 1, 0)
v[v$upp > STEP,]$upp <- 0
v[v$dwn > STEP,]$dwn <- 0
v$delta <- round(v$dwn - v$upp,digit)
v$cc <- c(0, v[1:(dim(v)[1]-1), "delta"])
v$cc <- round(sign(v$cc)*(abs(v$delta)),digit)
v$cc <- ifelse(v$delta<0,v$cc,0)
v$grp <- v[,1] + v$delta + v$cc
v$v <- v$grp
v <- v[order(v$v),]
if (all((v$delta + v$cc)==0)) break
}
v$one <- 1
t <- aggregate(v[,c("one","u")], by=list(v$flg,v$grp), FUN=sum)
names(t)[1:2] <- c("flg", "grp")
t$u <- t$u / t$one
t$one <- 1;
w <- aggregate(t[,c("one","u")], by=list(t$grp), FUN=median)
w$nm <- round(w$u/w$one,digit)
w2 <- aggregate(t[,c("one","u")], by=list(t$grp), FUN=sum)
DIF <- aggregate(v$d, by=list(v$grp), FUN=max)
condition.1 <- (w2$one >= length(pool)-max(DIF[,2])) | (w2$one > 6)
condition.2 <- (w2$one >= length(pool)-DIF[,2]) | (w2$one > 5)
sgstn <- unique(w[which(condition.1),"nm"])
sgstn2 <- unique(w[which(condition.2),"nm"])
cc <- STEP; tol <- STEP
itv <- setdiff(sgstn,sgstn2)
if (length(itv)>0) {
igv <- w[w$nm %in% itv,1]
tmp <- abs(sapply(igv, FUN="-", v$v))
tmp[round(tmp,8)==0] <- NA
cc <- max(apply(tmp, FUN="min", MARGIN=2, na.rm=TRUE))
flag <- matrix(rep(cc,2),ncol=2)
for (i in 1:length(cc)) {
flag[i,1] <- which(tmp[,i]==cc[i], arr.ind =TRUE)[1]
flag[i,2] <- i
}
cc <- ifelse(is.na(tmp[flag]),cc,tmp[flag])
if (cc >= 3*STEP) { tol <- cc ; cc <- STEP ; sgstn <- sgstn2 }
}
nearest <- function(x, Y, cc=1) {
na.omit(ifelse(abs(x - Y) > cc, NA, Y))[1]
}
if (length(sgstn)>0) {
v$ok <- sapply(v$u, FUN=nearest, sgstn, cc)
a[o] <- unlist(v$ok)
er <- aggregate(na.omit(abs(v$u - v$ok)),
by=list(v[!is.na(v$ok),]$ok), FUN=max)
er <- round(er$x, digits=nchar(v$ok[1])-1)
} else {
er <- NULL
}
return(list(suggest=sgstn, grp=unlist(a), tolerance=tol, error=er))
}
# -----------------------------------------------------------------------------
# Calculating of the projected parameters based on directional periodogram
#
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
kzpdr.proj <- function(rec = ls(1)) {
rec <- kzpdr.rec(rec)[,-c(3,4)]
fok <- kzpdr.pairs(rec)
ag1 <- fok$ag1
ag2 <- fok$ag2
s1 <- rec[round(rec$ang*180/pi,2) %in% ag1,]
s2 <- rec[round(rec$ang*180/pi,2) %in% ag2,]
s1$ag1 <- round(s1$ang*180/pi,2)
s1$ag2 <- s1$ag1 + 90
s2$ag2 <- round(s2$ang*180/pi,2)
s2$ag1 <- s2$ag2 - 90
s1$direction <- paste("(",s1$ag1,", ",s1$ag2,")",sep="")
s2$direction <- paste("(",s2$ag1,", ",s2$ag2,")",sep="")
names(s1) <- c("dir","f1","agp.1","g1","g2")
names(s2) <- c("dir","f2","agp.2","g2","g1")
s1$fp <- sapply(s1$f1/cos(s1$agp.1),min, 100)
s2$fp <- sapply(s2$f2/cos(s2$agp.2),min, 100)
s3 <- s1[s1$f1>0.29 & s1$fp>0.5 & s1$fp<=1,]
s4 <- s2[s2$f2>0.29 & s2$fp>0.5 & s2$fp<=1,]
# s3 <- s3[s3$g2>=45 & s3$g2<=90,]
# s4 <- s4[s4$g2>=45 & s4$g2<=90,]
s3$f1 <- 1 - s3$f1
s4$f2 <- 1 - s4$f2
s1 <- rbind(s1, s3)
s2 <- rbind(s2, s4)
cmb <- merge(s1[,1:5], s2[,1:5])
cmb$f0 <- sqrt(cmb$f1^2 + cmb$f2^2)
cmb$b1 <- cmb$g2>45 & cmb$g2<=90
cmb$sf <- ifelse(cmb$b1, sin(cmb$agp.2), -sin(cmb$agp.1))*cmb$f0
cmb$beta <- atan(cmb$f2/cmb$f1) # + cmb$agp.1
cmb$apha <- atan(cmb$f1/cmb$f2) # + cmb$agp.2
cmb$apha -> cmb$agp.2
cmb$beta -> cmb$agp.1
cmb$apha <- round(cmb$apha*180/pi, 0)
cmb$beta <- round(cmb$beta*180/pi, 0)
cmb <- cmb[,-c(4,6,8,9)]
s1 <- cmb[,c("dir", "g1", "g2", "sf", "agp.2", "apha")]
s2 <- cmb[,c("dir", "g1", "g2", "sf", "agp.1", "beta")]
names(s2) <- c("dir", "g1", "g2", "sf", "agp.2", "apha")
s1$formula <- "f2/f1"
s2$formula <- "f1/f2"
cmbf <- rbind(s1, s2)
cmbf$b2 <- ifelse(cmbf$formula=="f1/f2", cmbf$g1, cmbf$g2)
cmbf$apha <- cmbf$b2 + cmbf$apha
cmbf$agp.2 <- (cmbf$b2*pi/180) + cmbf$agp.2
s3 <- cmbf
s3$apha <- (s3$b2 - s3$apha)
s3$agp.2 <- (s3$b2*pi/180 - s3$agp.2)
cmbf <- rbind(cmbf, s3)
cmbf$apha <- ifelse(cmbf$apha > 90, cmbf$apha-180, cmbf$apha)
cmbf$apha <- ifelse(cmbf$apha < -90, (cmbf$apha)%%180, cmbf$apha)
cmbf$agp.2 <- ifelse(cmbf$agp.2 > pi/2, cmbf$agp.2-pi, cmbf$agp.2)
cmbf$agp.2 <- ifelse(cmbf$agp.2 < -pi/2, (cmbf$agp.2)%%pi, cmbf$agp.2)
cmbf <- cmbf[order(cmbf$dir,cmbf$sf,cmbf$agp.2),]
cmbf <- unique(cmbf[,1:6])
cmbf$sd <- shade(cmbf$apha, cmbf$g1)
return(cmbf)
}
# -----------------------------------------------------------------------------
# Find available records of directional periodograms in R environment
#
# @param rec Data frame or list of the outputs from function \code{kzpdr}.
# It includes the marked spike frequencies and sampling directions.
# Default is for searching from all variables in the R environment.
# @param md5 MD5 code for the specified data array. The returned data will be
# the records for the array with the same MD5 code. The default is
# to return the first record variable in alphabet order.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
kzpdr.rec <- function(rec = ls(1), md5="") {
if (missing(rec)) {
nm <- ls(1)
} else if (is.character(rec) & is.vector(rec)) {
nm <- rec
}
if (exists("nm")) {
oknm <- c("direction","freq","spg","max")
nc <- rep(FALSE,length(nm))
for (i in (1:length(nm))) {
rec <- get0(nm[i])
nc[i] <- is.list(rec) & !is.data.frame(rec)
nc[i] <- nc[i] & ("MD5" %in% names(rec))
nc[i] <- nc[i] & ("rec" %in% names(rec))
}
for (i in which(nc)) {
rec <- get0(nm[i])
nc[i] <- nc[i] & all(oknm %in% names(rec$rec))
if (nchar(md5)==32) {
nc[i] <- nc[i] & (rec$MD5 == md5)
}
}
if (length(which(nc)) > 1) {
cat("\nMultiple kzpdr objects. Use the first one: ")
cat(nm[which(nc)[1]],"\n\n")
} else if (length(nc)==0) {
return(NULL)
}
nc <- which(nc)
rec <- get0(nm[nc[1]])
md5 <- as.vector(rec$MD5)
# if (length(md5)>0) cat("\n",md5,"\n\n")
}
if (length(rec)==0) return(rec)
if ("MD5" %in% names(rec)) rec <- rec$rec
rec$ang <- as.numeric(gsub(enc2utf8("\xB0"), "", rec$direction))*pi/180
rec$ang <- ifelse(rec$ang<=(-pi/2), rec$ang%%pi, rec$ang)
rec[,1] <- paste(round(rec$ang*180/pi,2), enc2utf8("\xB0"), sep="")
rec <- unique(rec)
return(rec)
}
# -----------------------------------------------------------------------------
# Get all the pairs of directional periodograms
# and their counts for the spikes on each direction
#
# Returned value is a list of counts and direction pairs.
#
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
kzpdr.pairs <- function(rec = ls(1)) {
rec <- kzpdr.rec(rec)
fc <- aggregate(rec$freq, by=list(rec$direction),FUN=length)
names(fc) <- c("direction","count")
fc <- fc[order(as.numeric(gsub(enc2utf8("\xB0"), "", fc[,1]))),]
fg <- split(fc[,1], fc[,2])
rc <- as.numeric(gsub(enc2utf8("\xB0"), "", as.vector(fc[,1])))
rc -> tmpv
rc <- rev(rc[rc > 0])
rc <- rc[(rc - 90) %in% tmpv] # evaluate in pairs
ag2 <- rc; ag1 <- ag2 - 90;
return(list(fc=fc, fg=fg, rc=rc, ag1=ag1, ag2=ag2))
}
# -----------------------------------------------------------------
# Check and delete duplicated sampling directions
#
# @param angle Vector of sampling directions in radian.
# @rdname kzprj
# @export
# -----------------------------------------------------------------
ckGV <- function(angle, rec) {
if (length(rec)>0) {
relc <- rec
aglc <- unique(relc$direction)
okag <- paste(round((180/pi)*angle,2),enc2utf8("\xB0"),sep="")
okat <- which(okag %in% aglc)
if (length(okat)>0) angle <- angle[-okat]
}
return(angle)
}
# ------------------------------------------------------------------------------
# Check if spectral signals would sometimes disappear in given directions
#
# @param alpha Vector of checked directions in degree.
# @param gs Vector of directions for all available periodograms.
# @param cc Tolerance of wave parameters. Defaults to 3.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
shade <- function(alpha, gs, cc=3){
neighbors <- function(x, Y) {
if (sum(abs(Y-x)<= cc) == 0) return(0)
length(unique(Y[abs(Y-x)<= cc]))
}
if (length(alpha)==0) return(NA)
df <- data.frame(alpha=alpha, g1=gs, g2=90+gs)
df$c1 <- sapply(df$alpha, FUN=neighbors, df$g1)
df$c2 <- sapply(df$alpha, FUN=neighbors, df$g2)
df$c1 <- df$c1 + df$c2
df$c2 <- (abs(df$alpha - df$g1)<=cc) + (abs(df$alpha - df$g2)<=cc)
df$d1 <- df$c1 - df$c2
return(df$d1)
}
# ------------------------------------------------------------------------------
# Get the minimum gap of potential directions or frequencies
#
# Used to find the tolerance of the estimations.
# The argument \code{alpha} and \code{gs} could be frequencies for this function.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
tight <- function(alpha, gs, cc=3){
sv <- split(alpha, gs)
sv <- lapply(sv,sort)
svd <- lapply(sv,diff)
ok <- unlist(svd)
ok <- ok[round(ok,8)!=0]
ok <- ok[ok<cc*2]
if (length(ok)>0) {
mstp <- max(ok, na.rm = TRUE)
return(mstp*1.1)
} else {
svd <- diff(sort(alpha))
ok <- svd[round(svd,8)!=0]
ok <- ok[ok<cc*2]
if (length(ok)==0) return(ok)
Mstp <- max(ok, na.rm = TRUE)
return(Mstp*1.1)
}
}
# ------------------------------------------------------------------------------
# Calculate the tolerance of the estimations.
# @param tm2 Data frame. Internal data table for potential parameters.
# @param scale The scale of gridding data.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
tolerance <- function(tm2, scale, tol=1) {
rg.f <- split(tm2$freq, f=tm2$grp)
if (length(rg.f)>0) {
t.f <- max(sapply(rg.f, FUN="max") - sapply(rg.f, FUN="min"))
} else {
t.f <- scale[1]*tol
}
rg.d <- split(tm2$direction, f=tm2$grp)
if (length(rg.d)>0) {
t.d <- max(sapply(rg.d, FUN="max") - sapply(rg.d, FUN="min"))
} else {
t.d <- scale[2]*tol
}
return(c(t.f, t.d))
}
# ------------------------------------------------------------------------------
# Find the closure of nearest neighbors.
# @param mdc Integer. Specified approximate size of the closure.
# @param sc Numeric. Specified tolerance level of the cluster.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
closure <- function(tm2, mdc, sc) {
distance <- function(id, V, mdc) {
nm <- names(V)
nm <- nm[nm!="id"]
cpt <- V[V$id==id,nm]
sd <- cpt$sd
vx <- abs(V[,nm[1]] - cpt[,nm[1]])/0.5
vy <- abs(V[,nm[2]] - cpt[,nm[2]])/180
vz <- sqrt(vx^2 + vy^2)
vid <- V[order(vz)[1:max((mdc-sd),2)],"id"]
return(list(mdist=max(vz[vid]), mid=paste(sort(vid))))
}
nearest.neighbors <- function(id, V, mdc) {
nm <- names(V)
nm <- nm[nm!="id"]
vid <- as.numeric(unlist(V[V$id==id,"grpd"]))
v <- V[V$id %in% vid, ]
sd <- v[v$id==id, "sd"]
vx <- (paste(v$grpd) == paste(v[v$id==id,"grpd"]))
len <- sum(vx)
if (sum(vx) >= (mdc - sd)) {
group <- max(vid)
freq <- median(v$freq)
degree <- median(v$direction)
} else {
group <- 0
freq <- NA
degree <- NA
}
return(list(len=len, freq=freq, degree=degree, group=group))
}
dist <- sapply(tm2[,c("id")], tm2, FUN=distance, mdc)
tm2$grpd <- dist[2*(1:dim(tm2)[1])]
nbr <- sapply(tm2$id, tm2, FUN=nearest.neighbors, mdc)
tm2$nbr <- unlist(nbr["len",])
tm2$grp <- unlist(nbr["group",])
tm22 <- tm2[tm2$grp>0,]
rg.f <- split(tm22$freq, f=tm22$grp)
t.f <- lapply(lapply(rg.f, FUN="range"),FUN="diff")
rg.d <- split(tm22$direction, f=tm22$grp)
t.d <- lapply(lapply(rg.d, FUN="range"), FUN="diff")
rule <- (t.f < 2*sc[1]) & (t.d < 2*sc[2])
tm2$closure <- tm2$grp > 0
tm2$closest <- tm2$grp %in% names(rule[rule])
return(tm2)
}
# ------------------------------------------------------------------------------
# Exclude the duplicate supports in a cluster
#
# Duplication are supports for one potential parameter but from the same pair of
# directional periodogram, grouping under the given tolerance condition.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
distill <- function(tm2, mdc) {
tm2$nbr <- as.numeric(tm2$nbr)
tm2$grp <- as.numeric(tm2$grp)
tm2 <- tm2[order(tm2$nbr, tm2$grp, decreasing = TRUE),]
tmpv <- split(rownames(tm2),f=tm2$id)
tmpv <- as.vector(sapply(tmpv, FUN=function(m) m[1]))
tm2 <- tm2[rownames(tm2) %in% tmpv, 1:7]
# tmpv <- aggregate(cbind(tm2[,1:2]), by=list(tm2$grp,tm2$pair), FUN=mean)
# names(tmpv)[1:2] <- c("grp","pair")
# tmpw <- aggregate(cbind(tm2[,c(3,5)]), by=list(tm2$grp,tm2$pair), FUN=min)
# names(tmpw)[1:2] <- c("grp","pair")
# tmpv <- merge(tmpv, tmpw)
# tmpw <- aggregate(cbind(tm2[,c(3,6)]), by=list(tm2$grp,tm2$pair), FUN=sum)
# names(tmpw)[1:2] <- c("grp","pair")
# tmpv <- merge(tmpv, tmpw[,c(1,2,4)])
# tm2 <- tmpv[,names(tm2)]
tm2 <- tm2[order(tm2$grp),]
tmpv <- split(tm2$pair,f=tm2$grp)
tmpw <- sapply(sapply(tmpv,unique),FUN=function(m) !is.na(m))
tmpw <- as.vector(sapply(tmpw, sum))
tmpv <- as.vector(sapply(sapply(tmpv,FUN=function(m) !is.na(m)), sum))
tm2$nbr <- rep(tmpw, tmpv)
tm2 <- tm2[tm2$nbr >= mdc - tm2$sd, ]
if (dim(tm2)[1]==0) return(tm2)
tm2 <- tm2[order(tm2$grp, tm2$pair),]
tmpv <- split(tm2$pair, f=paste(tm2$pair, tm2$grp))
tmpv <- sapply(tmpv, FUN=function(m) !is.na(m))
tmpv <- lapply(lapply(tmpv, sum), ">", 1)
tmpv <- names(tmpv[unlist(tmpv)])
tm2$tw <- paste(tm2$pair, tm2$grp) %in% tmpv
tm2 <- tm2[!(tm2$tw & tm2$nbr==1), ]
tmpv <- split(tm2[!tm2$tw,]$freq, f=tm2[!tm2$tw,]$grp)
tmpw <- as.vector(sapply(tmpv, median))
tmpv <- split(tm2$freq, f=tm2$grp)
tmpv <- as.vector(sapply(lapply(tmpv, '>', -1),sum))
tm2$medf <- rep(tmpw, tmpv)
tmpw <- split(tm2[!tm2$tw,]$direction, f=tm2[!tm2$tw,]$grp)
tmpw <- as.vector(sapply(tmpw, median))
tm2$medd <- rep(tmpw, tmpv)
tm2$dif.f <- abs(tm2$freq - tm2$medf)
tm2$dif.d <- abs(tm2$direction - tm2$medd)
tm2 <- tm2[order(tm2$grp, tm2$pair, tm2$dif.f, tm2$dif.d),]
tm2$g1 <- substring(gsub(", [0-9]+)","",tm2$pair),2)
tm2$g2 <- gsub("[(+-]+[0-9]+, ","",tm2$pair)
tm2$g2 <- gsub("[)]","",tm2$g2)
tm2$g1 <- as.numeric(tm2$g1)
tm2$g2 <- as.numeric(tm2$g2)
tm2$dif.gf <- c(999, diff(tm2$grp))
tm2$dif.g1 <- c(999, diff(tm2$g1))
tm2$dif.g2 <- c(999, diff(tm2$g2))
tm2$c1 <- tm2$dif.gf==0 & tm2$dif.g1==0 & tm2$dif.g2==0
tm2 <- tm2[!(tm2$c1),1:10]
tm2 <- tm2[,-8]
return(tm2)
}
Try the kzfs package in your browser
Any scripts or data that you put into this service are public.
kzfs documentation built on June 2, 2019, 5:04 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# --------------------------------------------------------------------------
# @title
# Internal Function Used In Wave Parameter Estimation
#
# @description
# Function \code{commcore} is designed to find the common elements
# in a given list, which is the pool of all potential parameters. The
# common elements should appear in most data groups of paired directional
# periodograms.
#
# \code{kzpdr.rec} will search all variables in R global environment
# for available records of directional periodograms, and return the
# specified part as output.
#
# \code{kzpdr.pairs} returns a list of all available pairs of sampling
# directions and the counts number for the spikes on each direction.
#
# \code{kzpdr.prj} is used for calculating the projected parameters
# based on directional periodograms.
#
# Function \code{ckGV} refers to a given record set of directional
# periodogram and checks duplicated directions. The returned value is a
# list of sampling directions without those duplications.
#
# \code{shade} decides if a spectral signals would sometime disappear
# in the given sampling directions. The return value is a vector of 0/1.
# 0 means the signal should be seen in all given directions. 0 means not.
# 2 and larger may happen when there are very closed sampling angles.
#
# \code{tight} is used to get the minimum difference in a group of
# potential directions or frequencies.
#
# \code{tolerance} is used to calculate the tolerance of the estimations.
#
# Function \code{closure} returns the closure of the nearest neighbors.
#
# \code{distill} excludes the duplicate supports of a cluster. Usually,
# duplications are supports for the same parameters but from the same pairs,
# grouping under the given tolerance condition.
#
# @param pool List for the potential parameters.
# @param STEP The tolerance interval for a data element.
# Defaults to 1.
# @param DIFF The allowed absents times for a suggested
# common element. Default value is 1.
# @param digit Significant digits. Defaults to 0.
# @rdname kzprj
# @export
# @keywords internal
# ---------------------------------------------------------------------------
commcore <- function(pool, STEP=1, DIFF=1, digit=0) {
u <- as.vector(unlist(pool))
flg <- names(unlist(pool))
flg <- gsub(")[0-9]+",")",flg)
a <- round(u, digit)
o <- order(u)
v <- data.frame(v=a,u=u,d=DIFF,flg=flg)
v <- v[o,]
while (TRUE) {
v$upp <- c(999,diff(v[,1]))
v$dwn <- abs(c(rev(diff(rev(v[,1]))),999))
v$x <- ifelse((v$upp == 0 | v$dwn == 0), 1, 0)
v[v$upp > STEP,]$upp <- 0
v[v$dwn > STEP,]$dwn <- 0
v$delta <- round(v$dwn - v$upp,digit)
v$cc <- c(0, v[1:(dim(v)[1]-1), "delta"])
v$cc <- round(sign(v$cc)*(abs(v$delta)),digit)
v$cc <- ifelse(v$delta<0,v$cc,0)
v$grp <- v[,1] + v$delta + v$cc
v$v <- v$grp
v <- v[order(v$v),]
if (all((v$delta + v$cc)==0)) break
}
v$one <- 1
t <- aggregate(v[,c("one","u")], by=list(v$flg,v$grp), FUN=sum)
names(t)[1:2] <- c("flg", "grp")
t$u <- t$u / t$one
t$one <- 1;
w <- aggregate(t[,c("one","u")], by=list(t$grp), FUN=median)
w$nm <- round(w$u/w$one,digit)
w2 <- aggregate(t[,c("one","u")], by=list(t$grp), FUN=sum)
DIF <- aggregate(v$d, by=list(v$grp), FUN=max)
condition.1 <- (w2$one >= length(pool)-max(DIF[,2])) | (w2$one > 6)
condition.2 <- (w2$one >= length(pool)-DIF[,2]) | (w2$one > 5)
sgstn <- unique(w[which(condition.1),"nm"])
sgstn2 <- unique(w[which(condition.2),"nm"])
cc <- STEP; tol <- STEP
itv <- setdiff(sgstn,sgstn2)
if (length(itv)>0) {
igv <- w[w$nm %in% itv,1]
tmp <- abs(sapply(igv, FUN="-", v$v))
tmp[round(tmp,8)==0] <- NA
cc <- max(apply(tmp, FUN="min", MARGIN=2, na.rm=TRUE))
flag <- matrix(rep(cc,2),ncol=2)
for (i in 1:length(cc)) {
flag[i,1] <- which(tmp[,i]==cc[i], arr.ind =TRUE)[1]
flag[i,2] <- i
}
cc <- ifelse(is.na(tmp[flag]),cc,tmp[flag])
if (cc >= 3*STEP) { tol <- cc ; cc <- STEP ; sgstn <- sgstn2 }
}
nearest <- function(x, Y, cc=1) {
na.omit(ifelse(abs(x - Y) > cc, NA, Y))[1]
}
if (length(sgstn)>0) {
v$ok <- sapply(v$u, FUN=nearest, sgstn, cc)
a[o] <- unlist(v$ok)
er <- aggregate(na.omit(abs(v$u - v$ok)),
by=list(v[!is.na(v$ok),]$ok), FUN=max)
er <- round(er$x, digits=nchar(v$ok[1])-1)
} else {
er <- NULL
}
return(list(suggest=sgstn, grp=unlist(a), tolerance=tol, error=er))
}
# -----------------------------------------------------------------------------
# Calculating of the projected parameters based on directional periodogram
#
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
kzpdr.proj <- function(rec = ls(1)) {
rec <- kzpdr.rec(rec)[,-c(3,4)]
fok <- kzpdr.pairs(rec)
ag1 <- fok$ag1
ag2 <- fok$ag2
s1 <- rec[round(rec$ang*180/pi,2) %in% ag1,]
s2 <- rec[round(rec$ang*180/pi,2) %in% ag2,]
s1$ag1 <- round(s1$ang*180/pi,2)
s1$ag2 <- s1$ag1 + 90
s2$ag2 <- round(s2$ang*180/pi,2)
s2$ag1 <- s2$ag2 - 90
s1$direction <- paste("(",s1$ag1,", ",s1$ag2,")",sep="")
s2$direction <- paste("(",s2$ag1,", ",s2$ag2,")",sep="")
names(s1) <- c("dir","f1","agp.1","g1","g2")
names(s2) <- c("dir","f2","agp.2","g2","g1")
s1$fp <- sapply(s1$f1/cos(s1$agp.1),min, 100)
s2$fp <- sapply(s2$f2/cos(s2$agp.2),min, 100)
s3 <- s1[s1$f1>0.29 & s1$fp>0.5 & s1$fp<=1,]
s4 <- s2[s2$f2>0.29 & s2$fp>0.5 & s2$fp<=1,]
# s3 <- s3[s3$g2>=45 & s3$g2<=90,]
# s4 <- s4[s4$g2>=45 & s4$g2<=90,]
s3$f1 <- 1 - s3$f1
s4$f2 <- 1 - s4$f2
s1 <- rbind(s1, s3)
s2 <- rbind(s2, s4)
cmb <- merge(s1[,1:5], s2[,1:5])
cmb$f0 <- sqrt(cmb$f1^2 + cmb$f2^2)
cmb$b1 <- cmb$g2>45 & cmb$g2<=90
cmb$sf <- ifelse(cmb$b1, sin(cmb$agp.2), -sin(cmb$agp.1))*cmb$f0
cmb$beta <- atan(cmb$f2/cmb$f1) # + cmb$agp.1
cmb$apha <- atan(cmb$f1/cmb$f2) # + cmb$agp.2
cmb$apha -> cmb$agp.2
cmb$beta -> cmb$agp.1
cmb$apha <- round(cmb$apha*180/pi, 0)
cmb$beta <- round(cmb$beta*180/pi, 0)
cmb <- cmb[,-c(4,6,8,9)]
s1 <- cmb[,c("dir", "g1", "g2", "sf", "agp.2", "apha")]
s2 <- cmb[,c("dir", "g1", "g2", "sf", "agp.1", "beta")]
names(s2) <- c("dir", "g1", "g2", "sf", "agp.2", "apha")
s1$formula <- "f2/f1"
s2$formula <- "f1/f2"
cmbf <- rbind(s1, s2)
cmbf$b2 <- ifelse(cmbf$formula=="f1/f2", cmbf$g1, cmbf$g2)
cmbf$apha <- cmbf$b2 + cmbf$apha
cmbf$agp.2 <- (cmbf$b2*pi/180) + cmbf$agp.2
s3 <- cmbf
s3$apha <- (s3$b2 - s3$apha)
s3$agp.2 <- (s3$b2*pi/180 - s3$agp.2)
cmbf <- rbind(cmbf, s3)
cmbf$apha <- ifelse(cmbf$apha > 90, cmbf$apha-180, cmbf$apha)
cmbf$apha <- ifelse(cmbf$apha < -90, (cmbf$apha)%%180, cmbf$apha)
cmbf$agp.2 <- ifelse(cmbf$agp.2 > pi/2, cmbf$agp.2-pi, cmbf$agp.2)
cmbf$agp.2 <- ifelse(cmbf$agp.2 < -pi/2, (cmbf$agp.2)%%pi, cmbf$agp.2)
cmbf <- cmbf[order(cmbf$dir,cmbf$sf,cmbf$agp.2),]
cmbf <- unique(cmbf[,1:6])
cmbf$sd <- shade(cmbf$apha, cmbf$g1)
return(cmbf)
}
# -----------------------------------------------------------------------------
# Find available records of directional periodograms in R environment
#
# @param rec Data frame or list of the outputs from function \code{kzpdr}.
# It includes the marked spike frequencies and sampling directions.
# Default is for searching from all variables in the R environment.
# @param md5 MD5 code for the specified data array. The returned data will be
# the records for the array with the same MD5 code. The default is
# to return the first record variable in alphabet order.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
kzpdr.rec <- function(rec = ls(1), md5="") {
if (missing(rec)) {
nm <- ls(1)
} else if (is.character(rec) & is.vector(rec)) {
nm <- rec
}
if (exists("nm")) {
oknm <- c("direction","freq","spg","max")
nc <- rep(FALSE,length(nm))
for (i in (1:length(nm))) {
rec <- get0(nm[i])
nc[i] <- is.list(rec) & !is.data.frame(rec)
nc[i] <- nc[i] & ("MD5" %in% names(rec))
nc[i] <- nc[i] & ("rec" %in% names(rec))
}
for (i in which(nc)) {
rec <- get0(nm[i])
nc[i] <- nc[i] & all(oknm %in% names(rec$rec))
if (nchar(md5)==32) {
nc[i] <- nc[i] & (rec$MD5 == md5)
}
}
if (length(which(nc)) > 1) {
cat("\nMultiple kzpdr objects. Use the first one: ")
cat(nm[which(nc)[1]],"\n\n")
} else if (length(nc)==0) {
return(NULL)
}
nc <- which(nc)
rec <- get0(nm[nc[1]])
md5 <- as.vector(rec$MD5)
# if (length(md5)>0) cat("\n",md5,"\n\n")
}
if (length(rec)==0) return(rec)
if ("MD5" %in% names(rec)) rec <- rec$rec
rec$ang <- as.numeric(gsub(enc2utf8("\xB0"), "", rec$direction))*pi/180
rec$ang <- ifelse(rec$ang<=(-pi/2), rec$ang%%pi, rec$ang)
rec[,1] <- paste(round(rec$ang*180/pi,2), enc2utf8("\xB0"), sep="")
rec <- unique(rec)
return(rec)
}
# -----------------------------------------------------------------------------
# Get all the pairs of directional periodograms
# and their counts for the spikes on each direction
#
# Returned value is a list of counts and direction pairs.
#
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
kzpdr.pairs <- function(rec = ls(1)) {
rec <- kzpdr.rec(rec)
fc <- aggregate(rec$freq, by=list(rec$direction),FUN=length)
names(fc) <- c("direction","count")
fc <- fc[order(as.numeric(gsub(enc2utf8("\xB0"), "", fc[,1]))),]
fg <- split(fc[,1], fc[,2])
rc <- as.numeric(gsub(enc2utf8("\xB0"), "", as.vector(fc[,1])))
rc -> tmpv
rc <- rev(rc[rc > 0])
rc <- rc[(rc - 90) %in% tmpv] # evaluate in pairs
ag2 <- rc; ag1 <- ag2 - 90;
return(list(fc=fc, fg=fg, rc=rc, ag1=ag1, ag2=ag2))
}
# -----------------------------------------------------------------
# Check and delete duplicated sampling directions
#
# @param angle Vector of sampling directions in radian.
# @rdname kzprj
# @export
# -----------------------------------------------------------------
ckGV <- function(angle, rec) {
if (length(rec)>0) {
relc <- rec
aglc <- unique(relc$direction)
okag <- paste(round((180/pi)*angle,2),enc2utf8("\xB0"),sep="")
okat <- which(okag %in% aglc)
if (length(okat)>0) angle <- angle[-okat]
}
return(angle)
}
# ------------------------------------------------------------------------------
# Check if spectral signals would sometimes disappear in given directions
#
# @param alpha Vector of checked directions in degree.
# @param gs Vector of directions for all available periodograms.
# @param cc Tolerance of wave parameters. Defaults to 3.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
shade <- function(alpha, gs, cc=3){
neighbors <- function(x, Y) {
if (sum(abs(Y-x)<= cc) == 0) return(0)
length(unique(Y[abs(Y-x)<= cc]))
}
if (length(alpha)==0) return(NA)
df <- data.frame(alpha=alpha, g1=gs, g2=90+gs)
df$c1 <- sapply(df$alpha, FUN=neighbors, df$g1)
df$c2 <- sapply(df$alpha, FUN=neighbors, df$g2)
df$c1 <- df$c1 + df$c2
df$c2 <- (abs(df$alpha - df$g1)<=cc) + (abs(df$alpha - df$g2)<=cc)
df$d1 <- df$c1 - df$c2
return(df$d1)
}
# ------------------------------------------------------------------------------
# Get the minimum gap of potential directions or frequencies
#
# Used to find the tolerance of the estimations.
# The argument \code{alpha} and \code{gs} could be frequencies for this function.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
tight <- function(alpha, gs, cc=3){
sv <- split(alpha, gs)
sv <- lapply(sv,sort)
svd <- lapply(sv,diff)
ok <- unlist(svd)
ok <- ok[round(ok,8)!=0]
ok <- ok[ok<cc*2]
if (length(ok)>0) {
mstp <- max(ok, na.rm = TRUE)
return(mstp*1.1)
} else {
svd <- diff(sort(alpha))
ok <- svd[round(svd,8)!=0]
ok <- ok[ok<cc*2]
if (length(ok)==0) return(ok)
Mstp <- max(ok, na.rm = TRUE)
return(Mstp*1.1)
}
}
# ------------------------------------------------------------------------------
# Calculate the tolerance of the estimations.
# @param tm2 Data frame. Internal data table for potential parameters.
# @param scale The scale of gridding data.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
tolerance <- function(tm2, scale, tol=1) {
rg.f <- split(tm2$freq, f=tm2$grp)
if (length(rg.f)>0) {
t.f <- max(sapply(rg.f, FUN="max") - sapply(rg.f, FUN="min"))
} else {
t.f <- scale[1]*tol
}
rg.d <- split(tm2$direction, f=tm2$grp)
if (length(rg.d)>0) {
t.d <- max(sapply(rg.d, FUN="max") - sapply(rg.d, FUN="min"))
} else {
t.d <- scale[2]*tol
}
return(c(t.f, t.d))
}
# ------------------------------------------------------------------------------
# Find the closure of nearest neighbors.
# @param mdc Integer. Specified approximate size of the closure.
# @param sc Numeric. Specified tolerance level of the cluster.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
closure <- function(tm2, mdc, sc) {
distance <- function(id, V, mdc) {
nm <- names(V)
nm <- nm[nm!="id"]
cpt <- V[V$id==id,nm]
sd <- cpt$sd
vx <- abs(V[,nm[1]] - cpt[,nm[1]])/0.5
vy <- abs(V[,nm[2]] - cpt[,nm[2]])/180
vz <- sqrt(vx^2 + vy^2)
vid <- V[order(vz)[1:max((mdc-sd),2)],"id"]
return(list(mdist=max(vz[vid]), mid=paste(sort(vid))))
}
nearest.neighbors <- function(id, V, mdc) {
nm <- names(V)
nm <- nm[nm!="id"]
vid <- as.numeric(unlist(V[V$id==id,"grpd"]))
v <- V[V$id %in% vid, ]
sd <- v[v$id==id, "sd"]
vx <- (paste(v$grpd) == paste(v[v$id==id,"grpd"]))
len <- sum(vx)
if (sum(vx) >= (mdc - sd)) {
group <- max(vid)
freq <- median(v$freq)
degree <- median(v$direction)
} else {
group <- 0
freq <- NA
degree <- NA
}
return(list(len=len, freq=freq, degree=degree, group=group))
}
dist <- sapply(tm2[,c("id")], tm2, FUN=distance, mdc)
tm2$grpd <- dist[2*(1:dim(tm2)[1])]
nbr <- sapply(tm2$id, tm2, FUN=nearest.neighbors, mdc)
tm2$nbr <- unlist(nbr["len",])
tm2$grp <- unlist(nbr["group",])
tm22 <- tm2[tm2$grp>0,]
rg.f <- split(tm22$freq, f=tm22$grp)
t.f <- lapply(lapply(rg.f, FUN="range"),FUN="diff")
rg.d <- split(tm22$direction, f=tm22$grp)
t.d <- lapply(lapply(rg.d, FUN="range"), FUN="diff")
rule <- (t.f < 2*sc[1]) & (t.d < 2*sc[2])
tm2$closure <- tm2$grp > 0
tm2$closest <- tm2$grp %in% names(rule[rule])
return(tm2)
}
# ------------------------------------------------------------------------------
# Exclude the duplicate supports in a cluster
#
# Duplication are supports for one potential parameter but from the same pair of
# directional periodogram, grouping under the given tolerance condition.
# @rdname kzprj
# @export
# -----------------------------------------------------------------------------
distill <- function(tm2, mdc) {
tm2$nbr <- as.numeric(tm2$nbr)
tm2$grp <- as.numeric(tm2$grp)
tm2 <- tm2[order(tm2$nbr, tm2$grp, decreasing = TRUE),]
tmpv <- split(rownames(tm2),f=tm2$id)
tmpv <- as.vector(sapply(tmpv, FUN=function(m) m[1]))
tm2 <- tm2[rownames(tm2) %in% tmpv, 1:7]
# tmpv <- aggregate(cbind(tm2[,1:2]), by=list(tm2$grp,tm2$pair), FUN=mean)
# names(tmpv)[1:2] <- c("grp","pair")
# tmpw <- aggregate(cbind(tm2[,c(3,5)]), by=list(tm2$grp,tm2$pair), FUN=min)
# names(tmpw)[1:2] <- c("grp","pair")
# tmpv <- merge(tmpv, tmpw)
# tmpw <- aggregate(cbind(tm2[,c(3,6)]), by=list(tm2$grp,tm2$pair), FUN=sum)
# names(tmpw)[1:2] <- c("grp","pair")
# tmpv <- merge(tmpv, tmpw[,c(1,2,4)])
# tm2 <- tmpv[,names(tm2)]
tm2 <- tm2[order(tm2$grp),]
tmpv <- split(tm2$pair,f=tm2$grp)
tmpw <- sapply(sapply(tmpv,unique),FUN=function(m) !is.na(m))
tmpw <- as.vector(sapply(tmpw, sum))
tmpv <- as.vector(sapply(sapply(tmpv,FUN=function(m) !is.na(m)), sum))
tm2$nbr <- rep(tmpw, tmpv)
tm2 <- tm2[tm2$nbr >= mdc - tm2$sd, ]
if (dim(tm2)[1]==0) return(tm2)
tm2 <- tm2[order(tm2$grp, tm2$pair),]
tmpv <- split(tm2$pair, f=paste(tm2$pair, tm2$grp))
tmpv <- sapply(tmpv, FUN=function(m) !is.na(m))
tmpv <- lapply(lapply(tmpv, sum), ">", 1)
tmpv <- names(tmpv[unlist(tmpv)])
tm2$tw <- paste(tm2$pair, tm2$grp) %in% tmpv
tm2 <- tm2[!(tm2$tw & tm2$nbr==1), ]
tmpv <- split(tm2[!tm2$tw,]$freq, f=tm2[!tm2$tw,]$grp)
tmpw <- as.vector(sapply(tmpv, median))
tmpv <- split(tm2$freq, f=tm2$grp)
tmpv <- as.vector(sapply(lapply(tmpv, '>', -1),sum))
tm2$medf <- rep(tmpw, tmpv)
tmpw <- split(tm2[!tm2$tw,]$direction, f=tm2[!tm2$tw,]$grp)
tmpw <- as.vector(sapply(tmpw, median))
tm2$medd <- rep(tmpw, tmpv)
tm2$dif.f <- abs(tm2$freq - tm2$medf)
tm2$dif.d <- abs(tm2$direction - tm2$medd)
tm2 <- tm2[order(tm2$grp, tm2$pair, tm2$dif.f, tm2$dif.d),]
tm2$g1 <- substring(gsub(", [0-9]+)","",tm2$pair),2)
tm2$g2 <- gsub("[(+-]+[0-9]+, ","",tm2$pair)
tm2$g2 <- gsub("[)]","",tm2$g2)
tm2$g1 <- as.numeric(tm2$g1)
tm2$g2 <- as.numeric(tm2$g2)
tm2$dif.gf <- c(999, diff(tm2$grp))
tm2$dif.g1 <- c(999, diff(tm2$g1))
tm2$dif.g2 <- c(999, diff(tm2$g2))
tm2$c1 <- tm2$dif.gf==0 & tm2$dif.g1==0 & tm2$dif.g2==0
tm2 <- tm2[!(tm2$c1),1:10]
tm2 <- tm2[,-8]
return(tm2)
}
Try the kzfs package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.