Nothing
#
# Copyright 2014 Petr Matousu <pmatousu@more-praha.cz>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA 02110-1301, USA.
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: desc.std makes description metadata
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
desc.std <- function(Description=NA,Title=NA,Sample=NA,Unit=NA,Date=NA,x=NULL) {
# merges lists, if not NA y overwrites x
ml <- function(x, y) {
x.names <- names(x)
y.names <- names(y)
m.names <- sort(unique(c(x.names, y.names)))
sapply(m.names, function(i) {
if (is.list(x[[i]]) & is.list(y[[i]])) ml(x[[i]], y[[i]])
else if (i %in% y.names & !is.na(y[[i]])) y[[i]]
else x[[i]]
}, simplify = FALSE)
}
new <- list(Description=Description,
Title=Title,
Sample=Sample,
Unit=Unit,
Date=Date)
if(is.null(x)) {
ret <- new
} else {
stopifnot(is.std(x),length(x) == 1)
old <- list(Description=x[[1]]$Description,
Title=x[[1]]$Title,
Sample=x[[1]]$Sample,
Unit=x[[1]]$Unit,
Date=x[[1]]$Date)
ret <- ml(old,new)
}
return(lapply(ret,as.character))
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: create std object by hand
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
std <- function(a,r,desc=desc.std(),lmargs=list()) {
# a ... apertures
# r ... mass % retained on sieve
stopifnot(length(a) == length(r),is.list(lmargs),is.list(desc))
tab <- data.frame(a,r)
tab <- stats::na.exclude(tab)
metadata <- desc
names(tab) <- c("a","r")
if (!sum(tab$r) == 100) {
warning("Sum of retained is not 100, recalculating to percents, check results.")
tab$r <- 100*tab$r/sum(tab$r)
}
if (!all(tab$r > 0)) {
warning("Zero entries in retained field, excluding, check results.")
ind <- tab$r > 0
tab <- tab[ind,]
}
if (all(sapply(desc[-1],is.na)))
warning("Insufficient metadata.")
# Ordering
tab <- tab[order(tab$a,decreasing=T),]
# Cumulative Oversize
tab$O <- cumsum(tab$r)
# Based on approach, the last number of tab$O should be 100,
# as WAR of R acurracy, the 100 is forced
# to avoid log(log(100/100.00000000000001)) problem (NaNs produced)
tmp <- rev(tab$O)
if(tmp[1] - 100 < 1e-12) {
tmp[1] <- 100
tab$O <- rev(tmp)
}
# Undersize
tab$U <- 100 - tab$O
tmp <- rev(tab$U)
if(tmp[1] < 1e-12) {
tmp[1] <- 0
tab$U <- rev(tmp)
}
# log x
tab$X <- x2rrx(tab$a)
# log log y
tab$Y <- y2rry(tab$O)
# Linear fit
l <- do.call('lm', c(list(formula=Y ~ X, data=tab[tab$a > 0,]),lmargs))
# RR coefficients
xs <- rrx2x((-1)*l$coefficients[1]/l$coefficients[2])
names(xs) <- "xs"
ex <- l$coefficients[2]
names(ex) <- "ex"
RRcoeff <- list(xs=xs,ex=ex)
names(RRcoeff) <- c("diameter","exponent")
# size features
a_modus <- unname(xs*((ex-1)/ex)^(1/ex))
a_mean <- unname(xs*(log(2))^(1/ex))
o90um <- unname(orr(90,ex,xs))
size <- list(modus=a_modus,mean=a_mean,o90um=o90um)
ret <- list(stdata=tab,lmfit=l,RRcoefficients=RRcoeff,size=size)
ret <- list(c(metadata,ret))
# class
class(ret) <- c("std",class(ret))
return(ret)
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: read std data from CSV
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
read.std <- function(file, sep="\t", dec=".") {
stopifnot(file.exists(file))
# get number of columns
hd <- ncol(utils::read.table(file,header=F,sep=sep,stringsAsFactors=F,skip=7,nrows=1))
if(hd%%2 != 0) stop("Data file must have even number of columns!")
idx <- matrix(seq_len(hd),ncol=2,byrow=T)
ret <- list()
for (i in seq_len(hd/2)) {
cC1 <- rep("NULL",hd)
cC2 <- rep("NULL",hd)
cC1[idx[i,]] <- "character"
cC2[idx[i,]] <- "numeric"
tab <- utils::read.table(file,header=T,sep=sep,dec=dec,stringsAsFactors=F,skip=7,colClasses=cC2)
meta <- utils::read.table(file,header=F,sep=sep,dec=dec,stringsAsFactors=F,skip=0,nrows=5,comment.char="#",colClasses=cC1)
metadata <- meta[,2]
names(metadata) <- meta[,1]
ret[[i]] <- std(tab[,1],tab[,2],desc=as.list(metadata))
}
do.call(c.std,ret)
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: modify std object
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
tweak.std <- function(x,desc=desc.std(x=x),lmargs=as.list(x[[1]]$lmfit$call[-c(1:3)])) {
stopifnot(is.std(x),length(x) == 1,is.list(lmargs),is.list(desc))
# get data
a <- x[[1]]$stdata$a
r <- x[[1]]$stdata$r
# make new object
std(a=a,r=r,desc=desc,lmargs=lmargs)
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: concatenate
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
c.std <- function(...) {
dots <- list(...)
if(all(sapply(dots,is.std))){
dots <- lapply(dots,unclass)
ret <- do.call(c,dots)
class(ret) <- c("std",class(ret))
return(ret)
} else stop("All arguments must have std class.")
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: subset
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
'[.std' <- function(x, i, ...) {
structure(NextMethod("["), class = class(x))
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: transform from usr to axis coordinates and vise versa
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
y2rry <- function(x) {log(log(100/x))}
#' @export
rry2y <- function(x) {100/exp(exp(x))}
#' @export
x2rrx <- function(x) {log(x)}
#' @export
rrx2x <- function(x) {exp(x)}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: plot
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
plot.std <- function(x,
type=c("rr","rrdist"),
n=30, lgd=T,
col=c("#1B9E77", "#D95F02", "#7570B3", "#E7298A",
"#66A61E", "#E6AB02", "#A6761D", "#666666"),
bx=1,
by=5,
lgd.x="topleft",
...)
{
# x ... std oject
# n ... hint on labels density
# rounds to nearest b
mround <- function(x,b=5){
ord <- floor(log(base=10,x))
base <- b*10^(ord-1)
base*round(x/base)
}
# pretty labels
pretty_loglog <- function(x,n,b) {
# x ... loglog numbers ~ stdata$Y
# n ... hint for count
rge <- range(x,finite=T,na.rm=T)
unique(mround(rry2y(seq(rge[1],rge[2],length.out=n)),b=b))
}
pretty_log <- function(x,n,b) {
# x ... log numbers ~ stdata$X
# n ... hint for count
rge <- range(x,finite=T,na.rm=T)
unique(mround(rrx2x(seq(rge[1],rge[2],length.out=n)),b=b))
}
stopifnot(is.std(x), is.character(type))
if (length(x) > length(col)) {
warning("There is more samples than colors in palette! Forcing rainbow and disabling legend.")
col <-grDevices::rainbow(length(x))
lgd <- F
}
dots <- list(...)
if (is.character(type)) {
type <- match.arg(type)
}
# legend entries
lgdtext <- sapply(x, FUN=function(x){with(x,paste(Title,Sample,Unit,Date,sep=" | "))})
bg <- col[seq_len(length(x))]
# point type
mypch <- seq_len(length(x))
lty <- 2
lty.d <- 3
lty.u <- 1
lty.o <- 2
fill <- NA
border <- NA
rr <- expression({
# arguments
if(is.null(dots$xlab)) dots$xlab <- "Diameter of mesh aperture [ um ]"
if(is.null(dots$ylab)) dots$ylab <- "Cumulative Oversize, [ mass % ]"
if(is.null(dots$main)) dots$main <- "Sieve Test - Rosin - Rammler Plot"
if(is.null(dots$main)) dots$main <- "Sieve Test - Rosin - Rammler Plot"
if(is.null(dots$cex.axis)) dots$cex.axis <- 0.7
if(is.null(dots$las)) dots$las <- 1
if(is.null(dots$lwd)) dots$lwd <- 1
# common limits
if (is.null(dots$xlim)) {
limX <- lapply(x, FUN=function(x){x$stdata$X})
limX <- range(unique(sort(do.call(c,limX))),finite=T,na.rm=T)
} else {
limX <- dots$xlim
}
if (is.null(dots$ylim)) {
limY <- lapply(x, FUN=function(x){x$stdata$Y})
limY <- range(unique(sort(do.call(c,limY))),finite=T,na.rm=T)
} else {
limY <- dots$ylim
}
# common labels
# xaxis
xl <- pretty_log(limX,n,bx)
xa <- x2rrx(xl)
# yaxis
yl <- pretty_loglog(limY,n,by)
ya <- y2rry(yl)
# base plot
arg <- c(list(limX, limY, type='n', axes=FALSE), dots)
do.call('plot.default',arg)
box()
axis(1,at=xa,labels=xl,las=2,cex.axis=dots$cex.axis)
axis(2,at=ya,labels=yl,las=2,cex.axis=dots$cex.axis)
abline(col="grey", v=xa, h=ya, lwd=0.5, lty=2)
abline(col="grey", h=0, lwd=1.5, lty=1)
for (i in seq_len(length(x))) {
xi <- x[[i]]
X <- xi$stdata$X
Y <- xi$stdata$Y
arg <- c(list(Y ~ X,bg=col[i],col=col[i],pch=mypch[i]),dots)
do.call('points',arg)
usr <- par('usr')
clip(min(X[is.finite(X)]),max(X),usr[3],usr[4])
abline(xi$lmfit, col=col[i], lty=lty.o, lwd=dots$lwd)
do.call("clip", as.list(par('usr')))
}
})
rrdist <- expression({
# arguments
if(is.null(dots$xlab)) dots$xlab <- "Diameter of mesh aperture [ um ]"
if(is.null(dots$ylab)) dots$ylab <- "CDF & CCDF Value (left), PDF Value (right)"
if(is.null(dots$main)) dots$main <- "Distribution Functions"
if(is.null(dots$log)) dots$log <- ""
if(is.null(dots$cex.axis)) dots$cex.axis <- 0.7
if(is.null(dots$las)) dots$las <- 1
if(is.null(dots$lwd)) dots$lwd <- 1
if(n < 1000 | is.null(n)) n <- 1000
limX <- lapply(x, FUN=function(x){x$stdata$a})
limX <- range(unique(sort(do.call(c,limX))),finite=T,na.rm=T)
limY <- c(0,1)
if (grepl(x=dots$log,pattern="x",ignore.case=T)) limX[1] <- limX[2]/n
# base plot
arg <- c(list(limY ~ limX, type='n', axes=FALSE),dots)
do.call('plot.default',arg)
i <- which(names(dots) == "log")
dots <- dots[-i]
for (i in seq_len(length(x))) {
xi <- x[[i]]
# maximalni velikost castice
xmax <- max(limX)
# rr parameters
ex <- xi$RRcoefficients$exponent
xs <- xi$RRcoefficients$diameter
yp <- xi$stdata$O/100
xp <- xi$stdata$a
# calc
# d curve
xx <- seq(0,xmax,length.out=n)
y_d <- drr(x=xx,ex=ex,xs=xs)
y_d_pl <- y_d/max(pretty(y_d))
# oversize, undersize
y_u <- urr(x=xx,ex=ex,xs=xs)
y_o <- orr(x=xx,ex=ex,xs=xs)
arg <- c(list(y_u ~ xx,type='l',col=col[i],lty=lty.u),dots)
do.call('lines',arg)
arg <- c(list(y_o ~ xx,type='l',col=col[i],lty=lty.o),dots)
do.call('lines',arg)
arg <- c(list(y_d_pl ~ xx,type='l',col=col[i],lty=lty.d),dots)
do.call('lines',arg)
arg <- c(list(y=yp,x=xp,col=col[i],pch=mypch[i]),dots)
do.call('points',arg)
}
grid()
box()
axis(1,cex.axis=dots$cex.axis)
axis(2,cex.axis=dots$cex.axis)
at <- axTicks(side=2)
lab <- at*max(pretty(y_d))
axis(4,cex.axis=dots$cex.axis,labels=lab,at=at)
mypch <- c(mypch,NA,NA,NA)
lgdtext <- c(lgdtext,"CDF (undersize)","CCDF (oversize vs. measured)","PDF (density)")
lty <- c(rep(NA,times=length(x)),lty.u,lty.o,lty.d)
})
cmd <- get(type)
eval(cmd)
if (lgd) {
graphics::legend(x=lgd.x, fill=fill, border=border, pch=mypch, pt.lwd=1, lty=lty, lwd=2 , col=bg, pt.bg=bg, legend=lgdtext,bty='n', cex=1.14*dots$cex.axis)
}
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: summary
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
summary.std <- function(object,...) {
su <- lapply(object,FUN=function(x){
with(x,
{
data.frame(Date,Title,Sample,Unit,Description,
RRxs=RRcoefficients$diameter,
RRex=RRcoefficients$exponent,
RRms=size$modus,
RRmn=size$mean,
RRo90umPPC=size$o90um*100,
stringsAsFactors=F)
})})
su <- do.call(rbind,su)
rownames(su) <- NULL
return(su)
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: test is std?
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
is.std <- function(x) {
nms <- c("Description", "Title", "Sample", "Unit", "Date", "stdata", "lmfit", "RRcoefficients","size")
is_list <- is.list(x)
all_lists <- all(sapply(x,is.list))
is_std <- inherits(x, "std")
all_nms <- all(sapply(x,FUN=function(x){all(names(x) %in% nms)}))
all_has_lm <- all(sapply(x,FUN=function(x){inherits(x$lmfit, "lm")}))
all_has_tab <- all(sapply(x,FUN=function(x){inherits(x$stdata, "data.frame")}))
is_list & all_lists & is_std & all_nms & all_has_lm & all_has_tab
}
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#: RR distribution
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#' @export
drr <- function(x,ex,xs){(ex/xs)*(x/xs)^(ex-1)*exp(- (x/xs)^ex)}
#' @export
urr <- function(x,ex,xs){1-exp(- (x/xs)^ex)}
#' @export
orr <- function(x,ex,xs){exp(- (x/xs)^ex)}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.