Commit 1f739fc3 authored by thomas.forbriger's avatar thomas.forbriger

[MERGE] (unix): merge with current version from TFSoftware

Merge remote-tracking branch 'TFSoftware'

The files now reflect the recent version in TFSoftware.
parents a577f9ae cc04e7bb
#!/bin/gawk -f
# this is <calexiterextract.awk>
# ----------------------------------------------------------------------------
# $Id: $
#
# Copyright (c) 2014 by Thomas Forbriger (BFO Schiltach)
#
# extract values for all iterations from a calex.out file
#
# ----
# 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 2 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, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# ----
#
# This program reads one calex.out log file and extracts the numerical values
# for the progress obtained in all inversion parameters during iteration. This
# values are provided on the terminal as an ASCII table appropriate for
# gnuplot script sg056iter.gpt or similar. The standard invocation is
#
# calexiterextract.awk calex.out > sg056iter.dat
#
# Invoking the program by
#
# calexiterextract.awk -v MODE="persearchrange" calex.out > sg056iter.dat
#
# provides the inversion progress in realtive parameter change per search
# range. This might be of interest, when checking appropriateness of search
# ranges. Notice the these values cannot be used to calculate the location of
# poles and zeroes!
#
# REVISIONS and CHANGES
# 20/01/2014 V1.0 Thomas Forbriger
#
# ============================================================================
#
# function to check derived final parameters
function abs(a) {
rv=a;
if (a<0) { rv = -1 * a; }
return rv;
}
#
function amax(a,b) {
rv = abs(a);
if (rv < abs(b)) { rv = abs(b); }
return rv;
}
#
function dev(a, b) {
deno=amax(a,b);
if (deno < 1.e-20)
{
rv = 0;
}
else
{
rv=(a - b) / deno;
}
return rv;
}
#
# ----------------------------------------------------------------------------
#
BEGIN {
#
# set to 1 for DEBUG output
DEBUG=0;
#
# set current section to undefined
SEC=-1;
# accepted section index values:
# 1: initial column headers
# 2: initial values
# 3: search ranges
# 4: iteration block
# 5: final column headers
# 6: final values
#
# set current iteration
ITER=0;
#
# number of columns in calex output
NCOL=0;
#
# initial iteration line pattern with non-matching character sequence
ITERPATTERN="^XXXXXXXXXXXX";
}
#
# ============================================================================
# start scanning
# --------------
#
# detect first line of section
# ============================
#
# header keys initial or final
# ----------------------------
/^ iter RMS/ {
if (SEC<0)
{
SEC=1;
NCOL=NF;
# extract column headers of first header line
for (i=1; i<=NF; ++i)
{
COLHEAD[i]=$i;
}
}
else
{
SEC=5;
# last iteration will be reported next, but with full final values
ITER -= 1;
ITERPATTERN=sprintf("^%5d ", ITER);
}
if (DEBUG>1) { print; }
next;
}
#
# initial values
# --------------
/^ 0 / {
SEC=2;
for (i=1; i<=NF; ++i)
{
INIVAL[i]=$i;
}
NINI=NF;
# printf "# ";
for (i=1; i<=NCOL; ++i)
{
printf("%s ",COLHEAD[i]);
}
printf"\n";
next;
}
# search ranges
# -------------
/^ \+- / {
if (DEBUG>1) { print "search range!"; }
if (SEC == 2)
{
NSR=2;
SEARCHRANGE[1]=0;
SEARCHRANGE[2]=0;
# output initial values
printf("%d ",INIVAL[1]);
for (i=2; i<=NCOL; ++i)
{
printf("%f ",INIVAL[i]);
}
printf "\n";
}
for (i=2; i<=NF; ++i)
{
SEARCHRANGE[NSR + i - 1]=$i;
}
NSR += NF - 1;
SEC=3;
ITER=1;
ITERPATTERN=sprintf("^%5d ", ITER);
if (DEBUG) { print "ITERPATTERN: \"" ITERPATTERN "\""; }
next;
}
# iteration
# ---------
$0 ~ ITERPATTERN {
if (SEC != 5)
{
SEC=4;
if (DEBUG>1) { print; print "in iter"; }
for (i=1; i<=NF; ++i)
{
VAL[i]=$i;
}
NVAL=NF;
ITER += 1;
ITERPATTERN=sprintf("^%5d ", ITER);
}
else
{
SEC=6;
if (DEBUG>1) { print; print "in iter"; }
for (i=1; i<=NF; ++i)
{
FINVAL[i]=$i;
}
NFINVAL=NF;
}
}
# operate on values within section
# ================================
# (no specific initial sequence of line)
/^ / {
switch (SEC) {
# initial column headers
case 1:
for (i=1; i<=NF; ++i)
{
COLHEAD[NCOL + i]=$i;
}
NCOL += NF;
break;
# initial values
case 2:
for (i=1; i<=NF; ++i)
{
INIVAL[NINI + i]=$i;
}
NINI += NF;
break;
# search ranges
case 3:
print "ERROR: search range line should not match!";
next;
break;
# iteration block
case 4:
for (i=1; i<=NF; ++i)
{
VAL[NVAL + i]=$i;
}
NVAL += NF;
break;
# final column headers
case 5:
next;
break;
# final values
case 6:
for (i=1; i<=NF; ++i)
{
FINVAL[NFINVAL + i]=$i;
}
NFINVAL += NF;
break;
default:
print "WARNING: missed block type " SEC "!";
next;
}
}
# ----------------------------------------------------------------------------
# check whether all values are read
# should be triggered also in cases where only one line per iteration is
# written
// {
if (SEC == 2)
{
if (NINI == NCOL)
{
if (MODE == "persearchrange")
{
# this might appear silly, but in 'persearchrange' mode the initial values are
# just zero
printf("%d ",0);
printf("%f ",0);
for (i=3; i<=NCOL; ++i)
{
printf("%f ",0);
}
printf "\n";
}
else
{
printf("%d ",0);
printf("%f ",INIVAL[2]);
for (i=3; i<=NCOL; ++i)
{
printf("%f ",INIVAL[i]);
}
printf "\n";
}
}
}
else if (SEC == 4)
{
if (NVAL == NCOL)
{
if (MODE == "persearchrange")
{
printf("%d ",VAL[1]);
printf("%f ",VAL[2]);
for (i=3; i<=NCOL; ++i)
{
printf("%f ",VAL[i]);
}
}
else
{
printf("%d ",VAL[1]);
printf("%f ",VAL[2]);
for (i=3; i<=NCOL; ++i)
{
printf("%f ",VAL[i] * SEARCHRANGE[i] + INIVAL[i]);
}
}
printf "\n";
# increase NVAL by 1 to avoid retriggering of this section
NVAL += 1;
}
}
else if (SEC == 6)
{
if (NFINVAL == NCOL)
{
printf("# %3s %6s %11s %11s %8s\n", "no.", "par",
"value", "deriv. val.", "residual");
printf("# final: %3d %6s %11.7f %11.7f %5.2f%\n", 2,
COLHEAD[2], FINVAL[2], VAL[2],
100*dev(VAL[2], FINVAL[2]));
for (i=3; i<=NCOL; ++i)
{
printf("# final: %3d %6s %11.7f %11.7f %5.2f%\n", i,
COLHEAD[i], FINVAL[i],
VAL[i] * SEARCHRANGE[i] + INIVAL[i],
100*dev(VAL[i] * SEARCHRANGE[i] + INIVAL[i], FINVAL[i]));
}
printf "\n";
NFINVAL = 0;
}
}
}
#
# ============================================================================
END {
if (DEBUG)
{
for (i=1; i<=NCOL; ++i)
{
print "COLHEAD[" i "]: " COLHEAD[i];
}
}
}
# ============================================================================
# Format of the iterations block in the file to be scanned. The hash tag must
# be deleted, i.e. each line starts with a single blank.
#
# -----
#
# misfit from start model:
# writing file rest
#
# iter RMS amp del plp hlp plp
# hlp plp hlp plp hlp
#
# 0 0.051710 1.000000 0.000000 9.099180 0.988060 8.833920
# 0.893550 8.258500 0.703380 7.390980 0.404802
# +- 0.500000 0.010000 1.000000 0.010000 1.000000
# +- 0.010000 1.000000 0.010000 1.000000 0.010000
#
#
# 1 0.047760 -0.039046 0.000052 0.001695 0.000142 0.001593
# 0.000144 0.001332 0.000146 0.000742 0.000144
# 2 0.044152 -0.005157 0.008957 0.286338 0.023581 0.270220
# 0.024013 0.227982 0.024580 0.130224 0.024404
# 3 0.040857 -0.038128 0.009001 0.287762 0.023702 0.271562
# 0.024136 0.229113 0.024706 0.130871 0.024527
#
# .
# .
# .
#
# iter RMS amp del plp hlp plp
# hlp plp hlp plp hlp
# 108 0.003147 0.982706 0.000682 11.099974 0.989818 10.726266
# 0.895326 9.878787 0.705166 8.423401 0.406519
#
# QUAD called 2670 times: converged
# output signal for final model...
# writing file synt
# residual error...
#
# -----
#
# ----- END OF calexiterextract.awk -----
#!/bin/gawk -f
# this is <calexoutextract.awk>
# ----------------------------------------------------------------------------
# $Id: $
#
# Copyright (c) 2013 by Thomas Forbriger (BFO Schiltach)
#
# Extract initial and final model parameters from calex output
#
# ----
# 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 2 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, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# ----
#
# This program reads one or more output log files of calex and passes the
# essential lines to the standard output while skipping most iterations.
# Actually the file is passed to the output including the first INI
# iterations. Further iteration are skipped until the last iteration. The
# last iteration and all following lines are passed to the output. The default
# for INI is 3.
#
# To print all calex results from subdirectories:
#
# calexoutextract.awk */calex.out
#
# To print the calex output including the first 20 iterations:
#
# calexoutextract.awk -v INI=20 */calex.out
#
# If variable PSOUT is set, Postscript output files are produced:
#
# calexoutextract.awk -v PSOUT=1 */calex.out
#
# REVISIONS and CHANGES
# 20/12/2013 V1.0 Thomas Forbriger
#
# ============================================================================
#
BEGIN {
#
# hot indicates that a part of the input file is scanned which has to be
# reported to the output
hot=1;
#
# switch to a2ps output is PSOUT != 0
if (PSOUT)
{
output="a2ps -o calex.out.ps --center-title=calex.out";
}
else
{
output="cat";
}
#
# if INI != 0 set a user selected number of iterations to be reported
if (INI)
{
FIRSTITER=INI+1;
}
else
{
FIRSTITER=4;
}
FIRSTITERLINE=sprintf("^%5d", FIRSTITER);
#
# reset iteration counter
ITER=0;
}
#
# ============================================================================
# main body of awk script (initial settings are done)
# ---------------------------------------------------
#
# in case we start a new file (FNR == 1)
# --------------------------------------
FNR == 1 {
if (PSOUT) {
close (output);
output=sprintf("a2ps -o %s.ps --center-title=%s", FILENAME, FILENAME);
}
print "\nNew file: " FILENAME;
print "--------------------------------------------------";
}
#
# ----------------------------------------------------------------------------
# stop reporting at largest iteration as selected in BEGIN block
# --------------------------------------------------------------
#
$0 ~ FIRSTITERLINE {
if (ITER) { hot=0; }
print "\n .\n .\n .\n" | output;
}
#
# ----------------------------------------------------------------------------
# restart reporting after iteration has finished and RMS is reported
# ------------------------------------------------------------------
#
/^ iter RMS/ {
ITER=1;
if (!hot) { hot=1; }
}
#
# ----------------------------------------------------------------------------
# report any content if hot
# -------------------------
#
// {
if (hot) { print | output; }
}
# ----- END OF calexoutextract.awk -----
......@@ -118,6 +118,10 @@ c routine from SEIFE. Other programs had been updated earlier.
c
c ********** Version for homogeneous-triaxial sensors ********** March 2011
c reads aus3 output files, transforms the signals to UVW, applies CALEX to each
c
c May 2014 (Thomas Forbriger): reset normalized model parameters when
c switching to next component
c
c one, finally transforms parameters back to XYZ (Method of Peter Davis)
c File names 'calex.*' chenged to 'trical.*'
......@@ -126,6 +130,29 @@ c use nas_system
c the 'use' statement is compiler specific, and not required for Linux
parameter (mpar=12,msys=25,ns=480000)
c needs heap 40000 kBytes for ns=480000, 60000 for ns=800000
c
c DANGER: subroutines (in particular subroutine trans) rely on number of
c samples being set to 480000; is this value is changed here only,
c arrays will become mis-aligned!
c
c affected are:
c 801: subroutine congrd(x,q,gnorm,iter,ein,aus,sum,einf,sta)
c 803: parameter (mpar=12,ns=480000)
c 857: subroutine mini(x,q,d,ein,aus,sum,einf,sta)
c 861: parameter (mpar=12,ns=480000)
c 955:c original 'spaghetti' version of subroutine mini (of around 1985)
c 958: parameter (mpar=12,ns=480000)
c 1038: subroutine quadr(x,quad,ein,aus,sum,einf,sta)
c 1049: parameter (mpar=12,ns=480000,msys=25)
c 1487: subroutine rekf1(dt,r,n,ein,sum)
c 1488: parameter(ns=480000)
c 1504: subroutine rekf2(sre,sim,rre,rim,n,ein,sum)
c 1505: parameter(ns=480000)
c 1532: subroutine polytrend(mm,x,n)
c 1534: parameter(ndi=8,ns=480000)
c 1596: subroutine trans(aus3,n)
c 1597: parameter(ns=480000)
c
implicit double precision (a-h,o-z)
character endtxt*9, titel*72,name*3,typ*3,plop*12
character*50 eing,ausx,ausy,ausz
......@@ -203,6 +230,12 @@ c Loop over U, V, W output signals
init=-1
mq=0
print *,'initial parameters: '
do k=1,mpar
x(k)=0.d0
enddo
c print *,(x(k),k=1,mpar)
c corner period and order of the anti-alias filter
call inpar(m,maxit,qac,finac,init,dt,ns1,ns2)
ns1=1
......@@ -664,7 +697,7 @@ c read data in SEIFE format
write(6,*) ' header: ',trim(text)
21 read(7,'(a)') zeile
if(zeile(1:1).eq.'%') goto 21
20 read(zeile,1) nn,iform,dt,tmin,tsec
read(zeile,1) nn,iform,dt,tmin,tsec
1 format(i10,a20,3f10.3)
if(nn.gt.n) then
write(*,*) 'sorry, too many data points. can handle only ', n
......@@ -689,13 +722,13 @@ c read data in ASL format (such as written by Quanterra's Cimarron)
text=zeil(1:72)
write(6,*) 'header: ',trim(text)
read(zeil,'(120a)') zarr
nn=0
do i=1,119
if(zarr(i).ne.b.and.zarr(i+1).eq.b) then
nn=nn+1
ix(nn)=i
endif
enddo
nn=0
do i=1,119
if(zarr(i).ne.b.and.zarr(i+1).eq.b) then
nn=nn+1
ix(nn)=i
endif
enddo
read(zeil(1:ix(1)),'(a)') code1
read(zeil(ix(1)+2:ix(2)),'(a)') code2
read(zeil(ix(2)+2:ix(7)),*) iyear,iday,ithr,itmin,itsec
......@@ -714,7 +747,7 @@ c read data in ASL format (such as written by Quanterra's Cimarron)
endif
dt=1./srate
tmin=tmin+60.*thr
tsec=tsec+t100/100.
c tsec=tsec+t100/100.
c read(7,*,err=25,end=23) (x(j),j=1,n)
c free-format read doesn't always work - compiler bug?
......@@ -961,7 +994,7 @@ c! maximum encountered
c! go on
if(ql.lt.qm.or.qr.lt.qm) goto 2
c! interval too small
11 if(step.lt.eight*finac) goto 8
if(step.lt.eight*finac) goto 8
c! divide interval
step=step/eight
......@@ -1226,7 +1259,7 @@ c if (icount.lt.3) write(6,547) i,name(ipar),x0(i),xx(i),rho(i)
xx(i)=x0(i)
c if (icount.lt.3) write(6,547) i,'***',x0(i),xx(i),rho(i)
endif
547 format(i2,2x,a3,3f10.3)
c 547 format(i2,2x,a3,3f10.3)
546 continue
c icount=icount+1
c if(icount.ge.3) stop
......@@ -1535,7 +1568,7 @@ c remove polynomial trend