From 3dc7b3156b7975a0826bbc79a4f3974bc6438a6f Mon Sep 17 00:00:00 2001 From: Huriel Reichel <53496565+hurielreichel@users.noreply.github.com> Date: Sat, 29 Jul 2023 19:00:53 +0200 Subject: [PATCH] Update eyefit.R with an autofit With the extension of the "autofit" parameter, with the eyefit function, it would be possible to automatically fit a variogram without the dynamic window. This is especially interesting for batch processes. Moreover, the automatic assigned values for the variogram parameters are pretty good, so... why not? :) --- R/eyefit.R | 374 ++++++++++++++++++++++++++++------------------------- 1 file changed, 200 insertions(+), 174 deletions(-) diff --git a/R/eyefit.R b/R/eyefit.R index d76fe8c..90d3f13 100644 --- a/R/eyefit.R +++ b/R/eyefit.R @@ -37,9 +37,9 @@ print.eyefit <- function(x, ...){ return(invisible()) } -eyefit <- function(vario, silent=FALSE){ +eyefit <- function(vario, silent=FALSE, autofit=FALSE){ if(!requireNamespace("tcltk", quietly=TRUE)) stop("package tcltk is required to run eyefit()") - + geterrmessage() done <- tclVar(0) @@ -48,7 +48,7 @@ eyefit <- function(vario, silent=FALSE){ ## .eyefit.tmp <<- list() dmax <- max(vario$u) - + kappa1 <- tclVar("0.5") kappa2 <- tclVar("1.0") kernel<- tclVar("exponential") @@ -67,200 +67,226 @@ eyefit <- function(vario, silent=FALSE){ tausq <- as.numeric(tclObj(nugget)) eval(substitute(plot(vario))) - + fit <- get("eyefit.tmp", envir=eyefit.env) lapply(fit, function(x) - lines.variomodel(seq(0, maxdist, l=100), cov.model=x$cov.model, - kappa=x$kappa, cov.pars=x$cov.pars, - nug=x$nug, max.dist=x$max.dist)) - + lines.variomodel(seq(0, maxdist, l=100), cov.model=x$cov.model, + kappa=x$kappa, cov.pars=x$cov.pars, + nug=x$nug, max.dist=x$max.dist)) + if (k == "gneiting.matern" |k == "gencauchy" ) - { - lines.variomodel(x=seq(0,maxdist,l=100), cov.model=k, kappa=c(kp1,kp2), - cov.pars = c(sigmasq,phi),nug=tausq,max.dist=maxdist) - } + { + lines.variomodel(x=seq(0,maxdist,l=100), cov.model=k, kappa=c(kp1,kp2), + cov.pars = c(sigmasq,phi),nug=tausq,max.dist=maxdist) + } else if (k == "powered.exponential" || k =="cauchy" || k == "matern") - { - lines.variomodel(x=seq(0,maxdist,l=100), cov.model=k, kappa=kp1, - cov.pars = c(sigmasq,phi),nug=tausq,max.dist=maxdist) - } + { + lines.variomodel(x=seq(0,maxdist,l=100), cov.model=k, kappa=kp1, + cov.pars = c(sigmasq,phi),nug=tausq,max.dist=maxdist) + } else lines.variomodel(x=seq(0,maxdist,l=100), cov.model=k, cov.pars=c(sigmasq,phi), nug=tausq,max.dist=maxdist) - + } redraw <- function(...) + { + var <- as.character(tclObj(kernel)) + if (var == "gneiting.matern" | var == "gencauchy") { - var <- as.character(tclObj(kernel)) - if (var == "gneiting.matern" | var == "gencauchy") - { - tkconfigure(entry.kappa1, state="normal") - tkconfigure(ts5, state="normal") - tkfocus(entry.kappa1) - tkconfigure(entry.kappa2, state="normal") - tkconfigure(ts6, state="normal") - } - else if (var == "powered.exponential" || var =="cauchy" || var == "matern") - { - tkconfigure(entry.kappa1,state="normal") - tkconfigure(ts5,state="normal") - tkfocus(entry.kappa1) - tkconfigure(entry.kappa2,state="disabled") - tkconfigure(ts6,state="disabled") - } - else - { - tkconfigure(ts5,state="disabled") - tkconfigure(ts6,state="disabled") - tkconfigure(entry.kappa1,state="disabled") - tkconfigure(entry.kappa2,state="disabled") - } - replot() - } - - base <- tktoplevel() - tkwm.title(base, "Eyefit 1.0") - - spec.frm <- tkframe(base,borderwidth=2) - left.frm <- tkframe(spec.frm) - right.frm <- tkframe(spec.frm) - - frame1 <- tkframe(left.frm, relief="groove", borderwidth=2) - tkpack(tklabel(frame1,text="Parameters"),fill="both",side="top") - - entry.mdist <-tkentry(frame1,width="4",textvariable=mdist) - tkpack(ts1<-tkscale(frame1, label="Max. Distance",command=replot, from=0.00, to=dmax, - showvalue=0, variable=mdist, - resolution=0.01, orient="horiz",relief="groove"), - fill="both",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") - tkpack(entry.mdist,fill="none",side="right") - - frame3 <- tkframe(left.frm, relief="groove", borderwidth=2) - tkpack(tklabel(frame3,text="Cov. Parameters"),fill="both",side="top") - - entry.sill <-tkentry(frame3,width="4",textvariable=sill) - tkpack(ts2<-tkscale(frame3, label="Sill (sigmasq):",command=replot, from=0.00, to=2*max(vario$v), - showvalue=0, variable=sill, - resolution=0.01, orient="horiz",relief="groove"), - fill="none",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") - tkpack(entry.sill,side="right") - - frame4 <- tkframe(left.frm, relief="groove", borderwidth=2) - entry.range <-tkentry(frame4,width="4",textvariable=range) - tkpack(ts3<-tkscale(frame4, label="Range (phi):",command=replot, from=0.00, to=2*dmax, - showvalue=1, variable=range, - resolution=0.01, orient="horiz",relief="groove"), - fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") - tkpack(entry.range,side="right") - - frame5 <- tkframe(left.frm, relief="groove", borderwidth=2) - tkpack(tklabel(frame5,text="Nugget"),fill="both",side="top") - entry.nugget <-tkentry(frame5,width="4",textvariable=nugget) - tkpack(ts4<-tkscale(frame5, label="Nugget (tausq):",command=replot, - from=0.00, to=2*max(vario$v), - showvalue=1, variable=nugget, - resolution=0.01, orient="horiz",relief="groove"), - fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") - tkpack(entry.nugget,side="right") - - frame6 <- tkframe(left.frm, relief="groove", borderwidth=2) - tkpack(tklabel(frame6,text="Kappa"),fill="both",side="top") - entry.kappa1<-tkentry(frame6,width="4",textvariable=kappa1,state="disabled") - tkpack(ts5<-tkscale(frame6, label="Kappa 1:",command=replot, - from=0.00, to=10.00, - showvalue=1, variable=kappa1,state="disabled", - resolution=0.01, orient="horiz",relief="groove"), - fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") - tkpack(entry.kappa1,side="right",fill="none") - - - frame7 <- tkframe(left.frm, relief="groove", borderwidth=2) - entry.kappa2 <-tkentry(frame7,width="4",textvariable=kappa2,state="disabled") - tkpack(ts6<-tkscale(frame7, label="Kappa 2:",command=replot, - from=0.00, to=10.00, - showvalue=1, variable=kappa2,state="disabled", - resolution=0.01, orient="horiz",relief="groove"), - fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") - tkpack(entry.kappa2,side="right",fill="none") - - frame2 <- tkframe(right.frm, relief="groove", borderwidth=2) - tkpack(tklabel(frame2, text="Function")) - for ( i in c("cauchy","gencauchy","circular","cubic","exponential","gaussian", - "gneiting","gneiting.matern","linear","matern","power", - "powered.exponential","pure.nugget","spherical","wave") ) { - - tmp <- tkradiobutton(frame2, command=redraw, - text=i, value=i, variable=kernel) - tkpack(tmp, anchor="w") - } - - OnOK <- function(){ - replot() - } - - OnQuit <- function(){ - tclvalue(done)<-2 - } - - OnClear <- function(aux=vario){ - assign("eyefit.tmp", list(), envir=eyefit.env) - plot(aux) + tkconfigure(entry.kappa1, state="normal") + tkconfigure(ts5, state="normal") + tkfocus(entry.kappa1) + tkconfigure(entry.kappa2, state="normal") + tkconfigure(ts6, state="normal") + } + else if (var == "powered.exponential" || var =="cauchy" || var == "matern") + { + tkconfigure(entry.kappa1,state="normal") + tkconfigure(ts5,state="normal") + tkfocus(entry.kappa1) + tkconfigure(entry.kappa2,state="disabled") + tkconfigure(ts6,state="disabled") + } + else + { + tkconfigure(ts5,state="disabled") + tkconfigure(ts6,state="disabled") + tkconfigure(entry.kappa1,state="disabled") + tkconfigure(entry.kappa2,state="disabled") + } + replot() } - OnSave <- function(){ + if (autofit == TRUE){ k <- as.character(tclObj(kernel)) kp1 <- as.numeric(tclObj(kappa1)) - if(k == "gneiting.matern") - kp2 <- as.numeric(tclObj(kappa2)) + if (k == "gneiting.matern") + kp2 <- as.numeric(tclObj(kappa2)) else kp2 <- NULL maxdist <- as.numeric(tclObj(mdist)) - sigmasq <- as.numeric(tclObj(sill)) + sigmasq <- as.numeric(tclObj(sill)) phi <- as.numeric(tclObj(range)) tausq <- as.numeric(tclObj(nugget)) - aux <- list(cov.model=k, cov.pars=c(sigmasq, phi), - nugget=tausq, kappa=c(kp1,kp2), lambda = vario$lambda, - trend = vario$trend, - practicalRange = practicalRange(cov.model=k, - phi = phi, kappa = kp1), - max.dist=maxdist) + aux <- list(cov.model = k, cov.pars = c(sigmasq, phi), + nugget = tausq, kappa = c(kp1, kp2), lambda = vario$lambda, + trend = vario$trend, practicalRange = practicalRange(cov.model = k, + phi = phi, kappa = kp1), max.dist = maxdist) oldClass(aux) <- "variomodel" - assign("eyefit.tmp", c(get("eyefit.tmp", envir=eyefit.env),list(aux)), envir=eyefit.env) - + assign("eyefit.tmp", c(get("eyefit.tmp", envir = eyefit.env), + list(aux)), envir = eyefit.env) replot() - } + if (!silent) { + fit <- get("eyefit.tmp", envir = eyefit.env) + oldClass(fit) <- "eyefit" + return(fit) + } else return(invisible()) + } else{ - tkpack(frame1, frame3, frame4, frame5, frame6, frame7, fill="x") - tkpack(frame2, fill="x") - tkpack(left.frm, right.frm,side="left", anchor="n") - - c.but <- tkbutton(base,text="Clear", - command=function(){OnClear(vario)}) - - q.but <- tkbutton(base, text="Quit", command=OnQuit) - - save.but <- tkbutton(base, text="Save", command=OnSave) - tkpack(spec.frm) - tkpack(q.but, side="right") - tkpack(c.but, side="left") - tkpack(save.but, side="right") - - replot() - - tkbind(entry.kappa1, "", function(){replot()}) - tkbind(entry.kappa2, "", function(){replot()}) - tkbind(base, "", function() tclvalue(done)<-2) - - tkwait.variable(done) + base <- tktoplevel() + tkwm.title(base, "Eyefit 1.0") - tkdestroy(base) - - if (!silent){ - fit <- get("eyefit.tmp", envir=eyefit.env) - oldClass(fit) <- "eyefit" - return(fit) + spec.frm <- tkframe(base,borderwidth=2) + left.frm <- tkframe(spec.frm) + right.frm <- tkframe(spec.frm) + + frame1 <- tkframe(left.frm, relief="groove", borderwidth=2) + tkpack(tklabel(frame1,text="Parameters"),fill="both",side="top") + + entry.mdist <-tkentry(frame1,width="4",textvariable=mdist) + tkpack(ts1<-tkscale(frame1, label="Max. Distance",command=replot, from=0.00, to=dmax, + showvalue=0, variable=mdist, + resolution=0.01, orient="horiz",relief="groove"), + fill="both",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") + tkpack(entry.mdist,fill="none",side="right") + + frame3 <- tkframe(left.frm, relief="groove", borderwidth=2) + tkpack(tklabel(frame3,text="Cov. Parameters"),fill="both",side="top") + + entry.sill <-tkentry(frame3,width="4",textvariable=sill) + tkpack(ts2<-tkscale(frame3, label="Sill (sigmasq):",command=replot, from=0.00, to=2*max(vario$v), + showvalue=0, variable=sill, + resolution=0.01, orient="horiz",relief="groove"), + fill="none",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") + tkpack(entry.sill,side="right") + + frame4 <- tkframe(left.frm, relief="groove", borderwidth=2) + entry.range <-tkentry(frame4,width="4",textvariable=range) + tkpack(ts3<-tkscale(frame4, label="Range (phi):",command=replot, from=0.00, to=2*dmax, + showvalue=1, variable=range, + resolution=0.01, orient="horiz",relief="groove"), + fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") + tkpack(entry.range,side="right") + + frame5 <- tkframe(left.frm, relief="groove", borderwidth=2) + tkpack(tklabel(frame5,text="Nugget"),fill="both",side="top") + entry.nugget <-tkentry(frame5,width="4",textvariable=nugget) + tkpack(ts4<-tkscale(frame5, label="Nugget (tausq):",command=replot, + from=0.00, to=2*max(vario$v), + showvalue=1, variable=nugget, + resolution=0.01, orient="horiz",relief="groove"), + fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") + tkpack(entry.nugget,side="right") + + frame6 <- tkframe(left.frm, relief="groove", borderwidth=2) + tkpack(tklabel(frame6,text="Kappa"),fill="both",side="top") + entry.kappa1<-tkentry(frame6,width="4",textvariable=kappa1,state="disabled") + tkpack(ts5<-tkscale(frame6, label="Kappa 1:",command=replot, + from=0.00, to=10.00, + showvalue=1, variable=kappa1,state="disabled", + resolution=0.01, orient="horiz",relief="groove"), + fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") + tkpack(entry.kappa1,side="right",fill="none") + + + frame7 <- tkframe(left.frm, relief="groove", borderwidth=2) + entry.kappa2 <-tkentry(frame7,width="4",textvariable=kappa2,state="disabled") + tkpack(ts6<-tkscale(frame7, label="Kappa 2:",command=replot, + from=0.00, to=10.00, + showvalue=1, variable=kappa2,state="disabled", + resolution=0.01, orient="horiz",relief="groove"), + fill="x",expand=1,padx=3,pady=2,ipadx=3,ipady=2,side="left") + tkpack(entry.kappa2,side="right",fill="none") + + frame2 <- tkframe(right.frm, relief="groove", borderwidth=2) + tkpack(tklabel(frame2, text="Function")) + for ( i in c("cauchy","gencauchy","circular","cubic","exponential","gaussian", + "gneiting","gneiting.matern","linear","matern","power", + "powered.exponential","pure.nugget","spherical","wave") ) { + + tmp <- tkradiobutton(frame2, command=redraw, + text=i, value=i, variable=kernel) + tkpack(tmp, anchor="w") + } + + OnOK <- function(){ + replot() + } + + OnQuit <- function(){ + tclvalue(done)<-2 + } + + OnClear <- function(aux=vario){ + assign("eyefit.tmp", list(), envir=eyefit.env) + plot(aux) + } + + OnSave <- function(){ + k <- as.character(tclObj(kernel)) + kp1 <- as.numeric(tclObj(kappa1)) + if(k == "gneiting.matern") + kp2 <- as.numeric(tclObj(kappa2)) + else kp2 <- NULL + maxdist <- as.numeric(tclObj(mdist)) + sigmasq <- as.numeric(tclObj(sill)) + phi <- as.numeric(tclObj(range)) + tausq <- as.numeric(tclObj(nugget)) + aux <- list(cov.model=k, cov.pars=c(sigmasq, phi), + nugget=tausq, kappa=c(kp1,kp2), lambda = vario$lambda, + trend = vario$trend, + practicalRange = practicalRange(cov.model=k, + phi = phi, kappa = kp1), + max.dist=maxdist) + oldClass(aux) <- "variomodel" + assign("eyefit.tmp", c(get("eyefit.tmp", envir=eyefit.env),list(aux)), envir=eyefit.env) + + replot() + } + + tkpack(frame1, frame3, frame4, frame5, frame6, frame7, fill="x") + tkpack(frame2, fill="x") + tkpack(left.frm, right.frm,side="left", anchor="n") + + c.but <- tkbutton(base,text="Clear", + command=function(){OnClear(vario)}) + + q.but <- tkbutton(base, text="Quit", command=OnQuit) + + save.but <- tkbutton(base, text="Save", command=OnSave) + tkpack(spec.frm) + tkpack(q.but, side="right") + tkpack(c.but, side="left") + tkpack(save.but, side="right") + + replot() + + tkbind(entry.kappa1, "", function(){replot()}) + tkbind(entry.kappa2, "", function(){replot()}) + tkbind(base, "", function() tclvalue(done)<-2) + + tkwait.variable(done) + + tkdestroy(base) + + if (!silent){ + fit <- get("eyefit.tmp", envir=eyefit.env) + oldClass(fit) <- "eyefit" + return(fit) + } + else + return(invisible()) } - else - return(invisible()) }