SPlus implementation of hellinger routines THIS SITE IS UNDER CONSTRUCTION The following SPlus macros are included HELLDATIN Reads a datafile in the following format line 1: nsect (number of sections) lines 2 and following: one line for each data point, each line has five items: iside (side number), isect (section number), latitude, longitude, standard deviation The output is an Splus list object suitable for use in HELLINGER1 or HELLINGER3. If iside > 3 or isect > nsect, the point is discarded. HELLINGER1 Fits a single spreading ridge. Input is output list object from HELLDATIN, together with an initial guess for axis latitude, longitude, and angle of rotation. There is no grid search here so the initial guess should be good! Output is an Splus list object which contains all the objects of the input list object together with various fit results. MDATAPLOT and MDATAPLOT2 plot the output of HELLDATIN (for simple ridges) in Mercator projection. The sections are distinguished by letter and the side by case. MDATAPLOT reconstructs side 2 to side 1, MDATAPLOT2 plots data in original locations. MDATAPLOT and MDATAPLOT2 produce postscript files. These functions require the Splus maps library which, although not supported, is included with the Splus software. HELLINGER3 Same as HELLINGER1, but for triple junctions BINGHAM Calculations necessary for a plot of the confidence region. Default level is 95%, but can be changed by including an argument alpha. Input is an Splus list object such as that created by HELLINGER1. Output is an Splus list object containing the coordinates of the boundary, the upper and lower data matrix, together with variables icase and ncurve covering the pathological cases allowed in the original BINGHAM program (e.g. confidence region includes north pole). Can be used for triple junctions by setting indside to either c(k1,k2), where k1 and k2 are 1,2, or 3. The result will be a confidence region from side k1 to side k2. Can be used to change coordinates so that latitude and longitude are expressed relative to a new "North pole" with standard coordinates (plat,plong) and a new "Greenwich" with standard coordinates (glat,glong). BINGPLOT plots the output of BINGHAM as a postscript file. In addition it saves the coordinates xlong and ylat of the grid matrices so that the user can easily run the Splus macro contour with changed plot parameters. LONGRHO takes the output from BINGHAM and prepares an Splus list object appending bound.long.rho and bound.lat.rho which are simultaneous 95% confidence regions for long-rho and lat-rho, respectively. RUNNING THE PROGRAMS IN THE Splus ENVIRONMENT: The heavy computation is done with fortran subroutines in the file sub_hellinger.o. These need to be compiled using the UNIX command Splus COMPILE sub_hellinger.f Splus stores everything in a directory called .Data. By default, .Data is placed at the top level of the user's main directory. However if a .Data directory is created in the current working directory, Splus will use it. Since Splus saves all variables and references them in subsequent Splus sessions, it is useful to put different projects in different directories and to create separate .Data directories for each of them. In most installations the UNIX command Splus will run the Splus software. From within Splus 1. Load the FORTRAN utilities using the Splus command dyn.load("sub_hellinger.o") This step must be repeated at the start of each Splus session, but look at the Splus documentation under help(.First) for suggestions on how to automate step 1. 2. Load the Splus macros source("helldatin.splus") together with similar commands for the remaining macro files. This step only needs to be done once for each data file (since Splus stores the macros in the .Data directory.) 3. The following Splus session takes the data in the file "testdat" and loads it into the Splus object testdata, calls hellinger1 to fit a reconstruction from side 1 to side 2 and appends the fit parameters to the Splus object testdata, plots the data in the file "test.ps", calls bingham and puts the results in the Splus object testdata$bingham, plots a 95% confidence region in the file "test2.ps", calls bingham with plat=45, plong=0, glat=0, glong=180 and puts the results in testdata$bingham.rot, calls LONGRHO to get a 95% confidence region for axis longitude and rho relative to the pole at latitude 45, longitude 0, and plots these putting the results in "test3.ps". *******************journalized Splus session************************** S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc. S : Copyright AT&T. Version 3.4 Release 1 for Sun SPARC, SunOS 5.3 : 1996 Working data will be in .Data > dyn.load("sub_hellinger.o") > testdata<-HELLDATIN("test.dat",names=c("Africa","N. Amer.")) > #let's see what is now in testdata > names(testdata) [1] "filename" "names" "nsect" "ndat" "iside" "isect" [7] "xlat" "xlong" "sd" "call.l" "date.l" > testdata$nsect [1] 22 > testdata$call.l HELLDATIN("test.dat", names = c("Africa", "N. Amer.")) > testdata$date.l [1] "Thu Oct 14 16:17:21 EDT 1999" > #Let's try a fit with initial guess axis lat = 75, axis long = -90, rho = -5 > testdata<-HELLINGER1(testdata,alat,75,-90,-5) > names(testdata) [1] "filename" "names" "nsect" "ndat" "iside" "isect" [7] "xlat" "xlong" "sd" "call.l" "date.l" "alat" [13] "along" "rho" "etalat" "etalong" "resi" "ahat" [19] "kappahat" "df" "Xmat" "H11.2" "ier" "call" [25] "date" > #look at all the goodies now in testdata, what is the fit? > testdata$alat [1] -77.5465 > testdata$along [1] -94.10024 > testdata$rho [1] 2.403236 > MDATAPLOT(testdata,"test.ps") Generated postscript file "test.ps". null device 1 Warning messages: The functions and datasets in library section maps are not supported by StatSci. in: library(maps) > #what a bummer; but AT&T Bell Labs wrote the maps library so hopefully it is good > #run bingham and store the results in testdata > testdata$bingham <- BINGHAM(testdata) > #hopefully this is a nice football shaped region (icase = 1) > names(testdata) [1] "filename" "names" "nsect" "ndat" "iside" "isect" [7] "xlat" "xlong" "sd" "call.l" "date.l" "alat" [13] "along" "rho" "etalat" "etalong" "resi" "ahat" [19] "kappahat" "df" "Xmat" "H11.2" "ier" "call" [25] "date" "bingham" > names(testdata$bingham) [1] "filename" "names" "indside" "kappahat" "df" "alat" [7] "along" "rho" "ahat" "H11.2" "icase" "ncurve" [13] "nptbd" "bound" "uppmat" "lowmat" "rotmat" "ier" [19] "call.b" "date.b" > testdata$bingham$icase [1] 1 > #so far so good - now to plot it > BINGPLOT(testdata$bingham,"test2.ps") Generated postscript file "test2.ps". $xlong: [1] -109.00 -108.87 -108.74 -108.61 -108.48 -108.35 -108.22 -108.09 -107.96 [10] -107.83 -107.70 -107.57 -107.44 -107.31 -107.18 -107.05 -106.92 -106.79 [19] -106.66 -106.53 -106.40 -106.27 -106.14 -106.01 -105.88 -105.75 -105.62 [28] -105.49 -105.36 -105.23 -105.10 -104.97 -104.84 -104.71 -104.58 -104.45 [37] -104.32 -104.19 -104.06 -103.93 -103.80 -103.67 -103.54 -103.41 -103.28 [46] -103.15 -103.02 -102.89 -102.76 -102.63 -102.50 -102.37 -102.24 -102.11 [55] -101.98 -101.85 -101.72 -101.59 -101.46 -101.33 -101.20 -101.07 -100.94 [64] -100.81 -100.68 -100.55 -100.42 -100.29 -100.16 -100.03 -99.90 -99.77 [73] -99.64 -99.51 -99.38 -99.25 -99.12 -98.99 -98.86 -98.73 -98.60 [82] -98.47 -98.34 -98.21 -98.08 -97.95 -97.82 -97.69 -97.56 -97.43 [91] -97.30 -97.17 -97.04 -96.91 -96.78 -96.65 -96.52 -96.39 -96.26 [100] -96.13 -96.00 -95.87 -95.74 -95.61 -95.48 -95.35 -95.22 -95.09 [109] -94.96 -94.83 -94.70 -94.57 -94.44 -94.31 -94.18 -94.05 -93.92 [118] -93.79 -93.66 -93.53 -93.40 -93.27 -93.14 -93.01 -92.88 -92.75 [127] -92.62 -92.49 -92.36 -92.23 -92.10 -91.97 -91.84 -91.71 -91.58 [136] -91.45 -91.32 -91.19 -91.06 -90.93 -90.80 -90.67 -90.54 -90.41 [145] -90.28 -90.15 -90.02 -89.89 -89.76 -89.63 -89.50 -89.37 -89.24 [154] -89.11 -88.98 -88.85 -88.72 -88.59 -88.46 -88.33 -88.20 -88.07 [163] -87.94 -87.81 -87.68 -87.55 -87.42 -87.29 -87.16 -87.03 -86.90 [172] -86.77 -86.64 -86.51 -86.38 -86.25 -86.12 -85.99 -85.86 -85.73 [181] -85.60 -85.47 -85.34 -85.21 -85.08 -84.95 -84.82 -84.69 -84.56 [190] -84.43 -84.30 -84.17 -84.04 -83.91 -83.78 -83.65 -83.52 -83.39 [199] -83.26 -83.13 -83.00 $ylat: [1] -80.00 -79.94 -79.88 -79.82 -79.76 -79.70 -79.64 -79.58 -79.52 -79.46 [11] -79.40 -79.34 -79.28 -79.22 -79.16 -79.10 -79.04 -78.98 -78.92 -78.86 [21] -78.80 -78.74 -78.68 -78.62 -78.56 -78.50 -78.44 -78.38 -78.32 -78.26 [31] -78.20 -78.14 -78.08 -78.02 -77.96 -77.90 -77.84 -77.78 -77.72 -77.66 [41] -77.60 -77.54 -77.48 -77.42 -77.36 -77.30 -77.24 -77.18 -77.12 -77.06 [51] -77.00 -76.94 -76.88 -76.82 -76.76 -76.70 -76.64 -76.58 -76.52 -76.46 [61] -76.40 -76.34 -76.28 -76.22 -76.16 -76.10 -76.04 -75.98 -75.92 -75.86 [71] -75.80 -75.74 -75.68 -75.62 -75.56 -75.50 -75.44 -75.38 -75.32 -75.26 [81] -75.20 -75.14 -75.08 -75.02 -74.96 -74.90 -74.84 -74.78 -74.72 -74.66 [91] -74.60 -74.54 -74.48 -74.42 -74.36 -74.30 -74.24 -74.18 -74.12 -74.06 [101] -74.00 > #if we wanted to save this output for future use (such as to redraw the plot): > testdata$bingham$xlong<-.Last.value$xlong > testdata$bingham$ylat<-.Last.value$ylat > names(testdata$bingham) [1] "filename" "names" "indside" "kappahat" "df" "alat" [7] "along" "rho" "ahat" "H11.2" "icase" "ncurve" [13] "nptbd" "bound" "uppmat" "lowmat" "rotmat" "ier" [19] "call.b" "date.b" "xlong" "ylat" > testdata$bingham$ylat [1] -80.00 -79.94 -79.88 -79.82 -79.76 -79.70 -79.64 -79.58 -79.52 -79.46 [11] -79.40 -79.34 -79.28 -79.22 -79.16 -79.10 -79.04 -78.98 -78.92 -78.86 [21] -78.80 -78.74 -78.68 -78.62 -78.56 -78.50 -78.44 -78.38 -78.32 -78.26 [31] -78.20 -78.14 -78.08 -78.02 -77.96 -77.90 -77.84 -77.78 -77.72 -77.66 [41] -77.60 -77.54 -77.48 -77.42 -77.36 -77.30 -77.24 -77.18 -77.12 -77.06 [51] -77.00 -76.94 -76.88 -76.82 -76.76 -76.70 -76.64 -76.58 -76.52 -76.46 [61] -76.40 -76.34 -76.28 -76.22 -76.16 -76.10 -76.04 -75.98 -75.92 -75.86 [71] -75.80 -75.74 -75.68 -75.62 -75.56 -75.50 -75.44 -75.38 -75.32 -75.26 [81] -75.20 -75.14 -75.08 -75.02 -74.96 -74.90 -74.84 -74.78 -74.72 -74.66 [91] -74.60 -74.54 -74.48 -74.42 -74.36 -74.30 -74.24 -74.18 -74.12 -74.06 [101] -74.00 > #let's look at everything relative to a pole at lat 45, long 0 > testdata$bingham.rot<- BINGHAM(testdata,plat=45,plong=0,glat=0,glong=180) > names(testdata) [1] "filename" "names" "nsect" "ndat" "iside" [6] "isect" "xlat" "xlong" "sd" "call.l" [11] "date.l" "alat" "along" "rho" "etalat" [16] "etalong" "resi" "ahat" "kappahat" "df" [21] "Xmat" "H11.2" "ier" "call" "date" [26] "bingham" "bingham.rot" > names(testdata$bingham.rot) [1] "filename" "names" "indside" "kappahat" "df" "alat" [7] "along" "rho" "ahat" "H11.2" "icase" "ncurve" [13] "nptbd" "bound" "uppmat" "lowmat" "rotmat" "ier" [19] "call.b" "date.b" > testdata$bingham.rot$call.b BINGHAM(testdata, plat = 45, plong = 0, glat = 0, glong = 180) > testdata$bingham.rot$icase [1] 1 > #yum yum > testdata$bingham$.rot<-LONGRHO(testdata$bingham.rot) > names(testdata$bingham.rot) [1] "filename" "names" "indside" "kappahat" [5] "df" "alat" "along" "rho" [9] "ahat" "H11.2" "icase" "ncurve" [13] "nptbd" "bound" "uppmat" "lowmat" [17] "rotmat" "ier" "call.b" "date.b" [21] "bound.lat.rho" "bound.long.rho" > #Let's plot this into a postscript file > ps.options(command="lpr ",paper="letter",pointsize=10,horizontal=F) > postscript(file="test3.ps",print.it=F) > par(oma=c(0,0,1.5,0),mfrow=c(2,1)) > matplot(testdata$bingham.rot$bound[,1],testdata$bingham.rot$bound[,2],t ype="l" ) > title("Axis latitude vs. Axis longitude") ,xlab="longitude",ylab="latitude") > matplot(testdata$bingham.rot$bound.long.rho[,1],testdata$bingham.rot$bo und.long.rho[,2],type="l") > title("Axis longitude vs angle of rotation",xlab="longitude",ylab="rho") > mtext(outer=T,"file test.dat, 95% conf. region relative to pole (0,45) and 'Greenwich' (180,0), ", side=3, cex=1.5) > #plot is finished > dev.off() Generated postscript file "test3.ps". null device 1 > #you can look at the plot to see what all these commands do; I will go home > q() ***********************end of journalization************************** 4. You can open up a graphics window using the Splus command X11() In this case you see the results of your plot commands as you enter them. Just eliminate the ps.options, postscript and dev.off, commands from the above session. 5. If you wish to write the results to a file for use with the other programs (in KRCG_GJI.dir) try the Splus command help(write)