#!/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")