#!/usr/bin/Rscript ## Copyright Hugh Pumphrey 2012 ## 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 3 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, see ## . ## Gravity model for several simple line sources as per page xx fig xx ## of Lowrie. Each source has 3 parameters: depth Z, position x0 and ## radius R. The density contrast, which we keep the same for all ## sources, is entered by the user: you can really only estimate the ## mass/unit length of the source. ## No version numbers, but changes recorded here with date. ## 23 Nov 2012: Removed some uninformative printout from stdout. ## Added printout of fitted model values library(tcltk) ## Forward Model for line sources. ## The vector b is things that the FM needs but which are not to be ## estimated In this case the density contrast is in b because we can ## not estimate both it and the radii. The measurement vector is calculated ## from the formula on page XX of Lowrie. The K matrix is found by ## differentiating it analytically. formod <- function(xs,x,b,doks=FALSE){ ##print(paste("entering formod.xs has ",length(xs),"x has",length(x),"xs=")) ##print(xs) ns <- (length(xs)-1) %/% 3 ##print(paste("ns=",ns)) ## Big G, the gravitational constant G <- 6.67421e-5 ### SI * 1.0e6, so that y comes out in GU drho <- b[1] x0 <- xs[1:ns] Z <- xs[(ns+1):(2*ns)] R <- xs[(2*ns+1):(3*ns)] dg0 <- xs[3*ns+1] yf <- x * 0; for( isource in 1:(ns) ){ thisr <- R[isource] thisz<-Z[isource] thisx <- x0[isource] thisdg <- 2*pi*G*drho*thisr*abs(thisr) * thisz /(thisz^2 + (thisx-x)^2) ##print(paste("Adding source",thisr,thisz,thisx,"thisdg=")) ##print(thisdg) yf <- yf + thisdg } y <- dg0 + yf if(doks){ k <- array(0,c(length(y),length(xs))) ## Brute force Ks del <- 0.1 for(ipt in 1:length(xs)){ xpt <-xs xpt[ipt] <- xpt[ipt] + del ypt <- formod(xpt,x,b,doks=FALSE)$y k[,ipt] <- (ypt - y) / del } }else{ k <- NA } list(y=y,k=k) } ## Function to fit a set of observations to the model. We use the ## Marquadt-levenberg approach with a simple doubling / halving of the ## ML parameter gamma depending on whether the cost function goes down ## or up. Just using the simpler inverse Hessian approach only works ## if the starting point is _very_ close to the solution. nlls <- function(xs,x,y,b){ xi <- xs oldcost <- 1.0e31 gamma <- 1024 nx <- length(xs) D <- diag(nx) gamfac <- 2 convok <- FALSE for(i in 1:200){ fmr <- formod(xi,x,b,doks=TRUE) yc <- fmr$y k <- fmr$k dy <- y - yc oldcost <- t(dy) %*% dy ##print(k) ktk <- t(k) %*% k ##print(ktk) mli <- try(solve(ktk + gamma * D)) if(class(mli) == "try-error"){ print("Oh dear. Your problem is not well-posed.") ### print(ktk) xhat<-xs * NA break } xhat <- xi + mli %*% t(k) %*% (y - yc) dif <- xhat -xi crit <- t(dif) %*% dif ##print(paste("iteration",i,"convergence",crit,"gamma=",gamma)) ##print(xhat) if(crit < 1.e-3 && gamma < 1/1024){ print(paste("finished at iteration",i)) convok <-TRUE break } fmr <- formod(xhat,x,b,doks=FALSE) ynew <- fmr$y dy <- y - ynew newcost <- t(dy) %*% dy ##print(paste("old cost=",oldcost,"new cost=",newcost)) if( newcost < oldcost*1.1 ){ xi <- xhat gamma <- gamma / gamfac } else { gamma <- gamma * gamfac } } if(!convok) print("Failed to converge") return(list(xhat=xhat,cost=newcost,crit=crit,gamma=gamma,convok=convok)) } ## This is the "main function that is called when you press "run" lscalc <- function(){ ## Make sure that we have a window. But don't destroy and re-make ## if we do have one. print("========================") dvl <- dev.list() ##print("device list is:") ##print(dvl) if(length(dvl)<2 && is.null(dvl[1])) { x11(width=5,height=5,pointsize=12,xpos=-1) } par(mai=c(0.8,0.8,0.1,0.2)) layout(matrix(c(1,2),ncol=1)) ## This is where we get stuff from the Tcl controls Z <- eval(parse(text=paste("c(",tclvalue(tkget(Z.widget)),")" ) )) drho <- eval(parse(text=tclvalue(tkget(drho.widget)))) x0 <- eval(parse(text=paste("c(",tclvalue(tkget(x0.widget)),")" ) )) R <- eval(parse(text=paste("c(",tclvalue(tkget(R.widget)),")" ) )) offset <- eval(parse(text=paste("c(",tclvalue(tkget(offset.widget)),")" ) )) offset <- offset[1] ##print("Z,x0,R") ##print(Z) ##print(x0) ##print(R) ## Guard against the user entering incompatible lengths for Z,x0 nf <- min(c(length(Z), length(x0),length(R))) Z <- Z[1:nf] x0 <- x0[1:nf] R <- R[1:nf] ## Control whether we fit the data dofit <- TRUE ## These are not model parameters, but they need to come from somewhere. ## they are the limits of the range to be considered x1 <- eval(parse(text=tclvalue(tkget(x1.widget)))) x2 <- eval(parse(text=tclvalue(tkget(x2.widget)))) ### How the fuck do we get the values out of a text box ??? thisline <- "foobie bletch" xdat <-NA gdat <-NA lc <-0 for(nlines in 1:200){ ixs1 <- paste(nlines,".0",sep="") ixs2 <- paste(nlines,".500",sep="") inline <- tclvalue(tkget(dat.widget,ixs1,ixs2)) ##print(inline) if(inline == "") break if(length(grep("#",inline)) > 0 ) { ##print("Comment line: ignoring") next } chunx<-strsplit(inline,"[[:space:]=a-z_]") chunx<-chunx[[1]] chunx<-chunx[which(nchar(chunx) > 0)] chunx<-as.numeric(chunx) ##print(chunx) if(length(chunx) < 2) { ##print("Line witth less than 2 elements. Ignoring") next } lc <- lc+1 ##print(chunx[1:2]) xdat[lc] <- chunx[1] gdat[lc] <- chunx[2] } npp <- 201 ## This is just for plotting purposes x <- seq(x1,x2,len=npp) xs <-c(x0,Z,R,offset) dg <- formod(xs,x,drho)$y ##print(dg) plot(x,dg,xlab="Distance / m",ylab="Bouguer anomaly / Grav units", type="l",xaxs="i",ylim=range(c(dg,gdat),na.rm=TRUE)) abline(h=0,lty="dotted",col="gray") if(length(xdat) > 1){ points(xdat,gdat,pch=19,col="blue") points(xdat,formod(xs,xdat,drho)$y) } ## fit the model parameters to the data if(dofit & length(gdat) > 3){ fr <- nlls(xs,xdat,gdat,drho) xhat <- fr$xhat if(!is.na(xhat[1])){ if(fr$convok){ fitcol <-"green" }else{ fitcol <- "red" } ## Fit succeeded. We plot it. points(xdat,formod(xhat,xdat,drho)$y,col=fitcol) lines(x ,formod(xhat,x,drho)$y,col=fitcol) } }else{ xhat <- NA } legend("topleft",c("Entered","Fitted OK","Fitted (conv fail)","Data"), col=c("black","green","red","blue"), lty=c("solid","solid","solid","blank"),pch=c(1,1,1,19),cex=0.7, bg="#ffffff7f") if(!is.na(xhat[1])){ x0r <- xhat[1:nf] Zr<- xhat[(nf+1):(2*nf)] Rr<- xhat[(2*nf+1):(3*nf)] bot <- -( max(c(Z,Zr)) + max(c(R,Rr)) ) print("-----------------------") print("Original horizontal positions") print(xs[1:nf]) print("Original depths") print(xs[(nf+1):(2*nf)]) print("Original radii (-ve ==> -ve density contrast)") print(xs[(2*nf+1):(3*nf)]) print(paste("Original grav offset =",xs[3*nf+1],"GU")) print("Fitted horizontal positions") print(xhat[1:nf]) print("Fitted depths") print(xhat[(nf+1):(2*nf)]) print("Fitted radii (-ve ==> -ve density contrast)") print(xhat[(2*nf+1):(3*nf)]) print(paste("Fitted grav offset =",xhat[3*nf+1],"GU")) } else bot <- -(max(Z)+max(R)) plot(0,0,type="n",xlab="Distance / m", ylab="Depth / m",xlim=range(x), ylim=c(bot,0),xaxs="i",asp=1) ## polygon for ground xl <- -abs(x1*10) xr <- abs(x2*10) polygon(c(xl,xl,xr,xr,xl),5*c(bot,0,0,bot,bot),col="gray80") abline(h=0,lty="dashed") for(i in 1:nf){ if(R[i] < 0 ) fcol <- "skyblue" else fcol <- "pink" drawcircle(x0[i],-Z[i],abs(R[i]),col=fcol,border="black") } if(!is.na(xhat[1])){ for(i in 1:nf){ if(Rr[i] < 0 ) fcol <- "blue" else fcol <- "red" drawcircle(x0r[i],-Zr[i],abs(Rr[i]),col=fcol,border=fcol,density=10) } # points(x0r,-Zr,cex=abs(Rr/max(Rr)),col="red") } legend("topleft",c("Entered (-)","Entered (+)","Fitted (-)","Fitted (+)"), col=c("skyblue","pink","blue","red"),pch=c(19,19,1,1),cex=0.7) dev.copy2pdf(file="gravls.pdf",width=7,height=6) # points(x0,-Z,cex=abs(R/max(R))) } ### --- end of function lscalc that does all the interesting stuff --- drawcircle <- function(xc,yc,r,npts=30,...){ ## Use "polygon" to draw a circle theta <- seq(0,2*pi,len=npts) px <-xc + r*cos(theta) py <- yc + r*sin(theta) polygon(px,py,...) } ### ---Main program. Set up TK window and enter event loop. ------- tt <- tktoplevel() ## Wizardry to prevent R from quitting immediately when run with Rscript, ## but to allow it to do so when the quit button is pressed. done <- tclVar(0) tclvalue(done) <- 0 ## capture destroy (e.g. [X] button on window controls) ## otherwise the tkwait hangs with nowhere to go tkbind(tt, "", function()tclvalue(done)<-2) ## function to quit R without any questions ## This works because we have a tkwait.variable() at the end of the code ## which waits for the Tcl variable "done" to be altered ## If run from R prompt, this does nothing. die <- function(){ print("fault gravity demo killed R") tclvalue(done)<-3 } tktitle(tt) <- "Simple Fault Gravity Model" labeltext <- 33 ## For some weird reason this has to be numeric, even ## though the value is never used. And it has to be a _different_ value to ## any of the other Tcl variables you use, or they become the _same_ ## Tcl variable. Very flaky! tclvalue(labeltext)<-"Fault Gravity Model" ## Define widgets like this ## These ones are frames to contain other widgets contpar.widget <- tkframe(tt,relief="groove",borderwidth=3) datf.widget <- tkframe(tt,relief="groove",borderwidth=3) ## Not sure what to use for this. It is for the entry of data dat.widget <- tktext(datf.widget,width=24) ## No easy way to label these text entry widgets. Option text= does ## nothing and option label= causes an error. So we just put a text ## widget above each one. Zlab.widget <- tklabel(contpar.widget,text="Line source depth / m") Z.widget <- tkentry(contpar.widget) drholab.widget <- tklabel(contpar.widget,text="Density Contrast / (kg/m^3) ") drho.widget <- tkentry(contpar.widget) x0lab.widget <- tklabel(contpar.widget,text="Source position / m") x0.widget <- tkentry(contpar.widget) Rlab.widget <- tklabel(contpar.widget,text="Source radius R / m") R.widget <- tkentry(contpar.widget) offsetlab.widget <- tklabel(contpar.widget,text="Grav anomaly offset") offset.widget <- tkentry(contpar.widget) x1lab.widget <- tklabel(contpar.widget,text="left-hand plot limit") x1.widget <- tkentry(contpar.widget) x2lab.widget <- tklabel(contpar.widget,text="right-hand plot limit") x2.widget <- tkentry(contpar.widget) ## Defaults don't change these tkdelete(x1.widget,0,200) tkinsert(x1.widget,0,"-2000") tkdelete(x2.widget,0,200) tkinsert(x2.widget,0,"3000") ## Bogus data inserted into data widget so we don't have to paste it ## in every time for testing. Prob. want to remove this later tkinsert(dat.widget,"1.0","-2000 5.3513 -1750 5.1155 -1500 4.4533 -1250 2.4417 -1000 -2.3527 -750 3.0104 -500 5.0561 -250 7.2824 0 7.8927 250 8.9243 500 10.829 750 14.023 1000 19.679 1250 29.225 1500 34.941 1750 28.818 2000 20.006 2250 14.576 2500 11.041 2750 9.2719 3000 7.4395" ) ## This one is a button to run the main function run.widget <- tkbutton(contpar.widget, text="Run",command=lscalc) ## This one is a button to exit from the program when run from Rscript. ## It does nothing mutch when run from the command line. die.widget <- tkbutton(contpar.widget, text="quit",command=die) ### Function to set all controls to the default value settodef <- function(){ ## How can we clear the text widget? tkdelete(Z.widget,0,200) tkinsert(Z.widget,0,"200,600") tkdelete(drho.widget,0,200) tkinsert(drho.widget,0,"400") tkdelete(offset.widget,0,200) tkinsert(offset.widget,0,"0") tkdelete(x0.widget,0,200) tkinsert(x0.widget,0,"-800,1200") tkdelete(R.widget,0,200) tkinsert(R.widget,0,"-150, 200") } ## Button to reset controls to default values defaults.widget <- tkbutton(contpar.widget, text="Defaults",command=settodef) ### Widgets all defined. Pack 'em in. settodef() tkpack(run.widget) tkpack(die.widget) tkpack(defaults.widget) tkpack(Zlab.widget) tkpack(Z.widget) tkpack(drholab.widget) tkpack(drho.widget) tkpack(x0lab.widget) tkpack(x0.widget) tkpack(Rlab.widget) tkpack(R.widget) tkpack(offsetlab.widget) tkpack(offset.widget) tkpack(x1lab.widget) tkpack(x1.widget) tkpack(x2lab.widget) tkpack(x2.widget) tkpack(dat.widget) tkpack(contpar.widget,side="left") tkpack(datf.widget,side="right") ## This waits until the Tcl variable "done" changes, either because ## the user closed the Tcl window with its [X] button or because he ## pressed the Tcl "quit" button tkwait.variable(done) if(tclvalue(done)=="2") stop("TK window closed") if(tclvalue(done)=="3") stop("Die button pressed")