[cig-commits] [commit] master: make obsres.py tolerant to wrong input directories; output warnings start with # for easy post-processing. (abc4c6a)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/relax

On branch  : master
Link       : https://github.com/geodynamics/relax/compare/392758611720e5bf2bede2fa63db44922d54edda...6b5c2f26e0e2cb532d632dff1e43c559bc030db4

>---------------------------------------------------------------

commit abc4c6a73aa21e5aaddde5f2897844fe94582ce4
Author: Sylvain Barbot <sbarbot at ntu.edu.sg>
Date:   Sun Aug 24 22:25:22 2014 -0700

    make obsres.py tolerant to wrong input directories; output warnings start with # for easy post-processing.


>---------------------------------------------------------------

abc4c6a73aa21e5aaddde5f2897844fe94582ce4
 util/obsres.py | 34 +++++++++++++++++++++++-----------
 1 file changed, 23 insertions(+), 11 deletions(-)

diff --git a/util/obsres.py b/util/obsres.py
index bd54472..a73e2e5 100755
--- a/util/obsres.py
+++ b/util/obsres.py
@@ -16,7 +16,8 @@ Options:
     --ddir=<data/path>            data directory.
  
 Description:
-    obsres.py computes residuals between Relax and GPS data time series. 
+    obsres.py computes the weighted sum of the squares of residuals between Relax and GPS data time 
+    series. 
 
     The --range option is for testing overall scaling of the model predictions in the time domain.
     For --range=-1/1/1, the program evaluates the residuals with the modified time t/tm, with tm 
@@ -86,7 +87,12 @@ def main():
 	for i in xrange(len(args)):
 		wdir=args[i]
 		fname=wdir+'/'+network
-		f=file(fname,'r')
+                try:
+		    f=file(fname,'r')
+                except IOError as e:
+                    print >> sys.stderr, '# obsres.py: could not find '+fname+', skipping.'
+                    continue
+
 		name=np.loadtxt(f,comments='#',unpack=True,dtype='S4',usecols=[3])
 		coverage=np.zeros(len(name))
 		#print name
@@ -97,6 +103,8 @@ def main():
 
 			# initialize norm of residuals
 			norm=0
+			norm0=0
+                        count=0
 			index=0
 			for s in name:
                             # load data (test upper and lower case)
@@ -118,7 +126,7 @@ def main():
                                 with open(fname) as f: pass
                             except IOError as e:
                                 # skipping station s
-                                print >> sys.stderr, 'obsres.py: could not find '+fname+', skipping.'
+                                print >> sys.stderr, '# obsres.py: could not find '+fname+', skipping.'
                                 continue
 			
                             f=file(fname,'r')
@@ -132,11 +140,12 @@ def main():
                                 se=np.atleast_1d(se)
                                 sd=np.atleast_1d(sd)
                             except:
-                                print >> sys.stderr, 'obsres.py: error loading file '+s+'. skipping.'
+                                print >> sys.stderr, '# obsres.py: error loading file '+s+'. skipping.'
                                 coverage[index]=1
                                 index+=1
                                 continue
 
+                            pos=[]
                             pos=np.logical_and(np.logical_and(np.isfinite(nr+er+dr),tr>=t1),tr<=t2)
                             tr=tr[pos]
                             nr=nr[pos]
@@ -145,10 +154,9 @@ def main():
                             sn=sn[pos]
                             se=se[pos]
                             sd=sd[pos]
-                            pos=[]
 
                             if 0==len(tr):
-                                print >> sys.stderr, 'obsres.py: skipping station '+s+' because of insufficient data coverage.'
+                                print >> sys.stderr, '# obsres.py: skipping station '+s+' because of insufficient data coverage.'
                                 continue
 
                             # load model
@@ -163,6 +171,7 @@ def main():
 
                             # time interval
                             tmax=min([max(tm),max(tr)])
+                            pos=[]
                             pos=tr<=tmax
                             coverage[index]=float(sum(pos))/float(len(pos))
                             index+=1
@@ -173,10 +182,9 @@ def main():
                             sn=sn[pos]
                             se=se[pos]
                             sd=sd[pos]
-                            pos=[]
 
                             if 0==len(tr):
-                                print >> sys.stderr, 'obsres.py: skipping station '+s+' because of insufficient model coverage. Try reducing scaling.'
+                                print >> sys.stderr, '# obsres.py: skipping station '+s+' because of insufficient model coverage. Try reducing scaling.'
                                 continue
 
                             if isrelax:
@@ -196,14 +204,18 @@ def main():
                                     em/=scale
                                     dm/=scale
 
-                            dnorm=sum(pow(nr-nm,2)/sn*wn+pow(er-em,2)/se*we+pow(dr-dm,2)/sd*wd)
+                            dnorm0=sum(pow(nr/sn,2)*wn+pow(er/se,2)*we+pow(dr/sd,2)*wd)
+                            dnorm=sum(pow((nr-nm)/sn,2)*wn+pow((er-em)/se,2)*we+pow((dr-dm)/sd,2)*wd)
+                            count+=wn*float(len(nm))+we*float(len(em))+wd*float(len(em))
                             norm+=dnorm
+                            norm0+=dnorm0
                             #print '{0:8.2e}'.format(dnorm)
 
                         if 1==exprange:
-                            print '{0}    {1:8.2e}'.format(wdir.ljust(max(map(len,args))+1),norm)
+                            print '{0}    {1:8.2e}'.format(wdir.ljust(max(map(len,args))+1),norm/count)
                         else:
-                            print '{0}    {1:8.2e}  {2:9.2e} {3:8.2e}'.format(wdir.ljust(max(map(len,args))+1),norm,tscale,sum(coverage)/float(len(coverage)))
+                            print '{0}    {1:8.2e}  {2:9.2e} {3:8.2e}'.format(wdir.ljust(max(map(len,args))+1),norm/count,tscale,sum(coverage)/float(len(coverage)))
+                            #print 1. - norm/norm0
                         sys.stdout.flush()
 
 if __name__ == "__main__":