Commit 131b81a8 authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

introduced square error test

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 810
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 8d62c805
c this is <gremlin.f>
c------------------------------------------------------------------------------
c $Id: gremlin.f,v 1.10 2000-06-02 19:42:00 thof Exp $
c $Id: gremlin.f,v 1.11 2002-04-19 15:20:33 forbrig Exp $
c
c 05/12/97 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -77,13 +77,15 @@ c 25/05/00 V4.10 changed dat_fcamp to calculate weighted least-squares
c 02/06/00 V4.11 hat to correct follow reporting schemem in
c mod_showparcor
c 02/06/00 V4.12 added mweight condition in res_opt
c 19/04/02 V4.13 support square error ense-test by command swense
c
c==============================================================================
c
program gremlin
c
character*79 version
parameter(version='GREMLIN V4.12 GREens Matrix Linearized INversion')
parameter(version=
& 'GREMLIN V4.13 GREens Matrix Linearized INversion')
c common blocks
include 'libs/glq_dim.inc'
include 'libs/glq_para.inc'
......@@ -494,7 +496,14 @@ c
elseif (arg(1:6).eq.'wense ') then
read(arg(7:80), *, err=98, end=98) ival,ival2,rval,rval2,bn
if (inv_mat()) then
call res_wresmod(ival, ival2, rval, rval2, bn)
call res_wresmod(ival, ival2, rval, rval2, bn, .true.)
else
print *,'calculation of partial derivatives failed'
endif
elseif (arg(1:6).eq.'swense ') then
read(arg(7:80), *, err=98, end=98) ival,ival2,rval,rval2,bn
if (inv_mat()) then
call res_wresmod(ival, ival2, rval, rval2, bn, .false.)
else
print *,'calculation of partial derivatives failed'
endif
......
c this is <gremlin_help.f>
c------------------------------------------------------------------------------
c $Id: gremlin_help.f,v 1.2 2000-05-24 18:41:34 thof Exp $
c $Id: gremlin_help.f,v 1.3 2002-04-19 15:20:34 forbrig Exp $
c
c 25/03/98 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -12,6 +12,7 @@ c 14/01/99 V1.1 removed link to old gremlin code
c 04/03/99 V1.2 file/mtt explanation
c 24/05/00 V1.3 - introduced reso option wgpd
c - introduced reso option wense
c 19/04/02 V1.4 introduced command swense
c
c==============================================================================
c
......@@ -109,6 +110,12 @@ c
print *,' wense f,t,nu,rms,name'
print *,' just like ''sense'' but puts rating for each'
print *,' parameter to a file with base ''name'' '
print *,' swense f,t,nu,err,name'
print *,' just like ''wense'' but checks square'
print *,' error rather than rms error'
print *,' ''err'' is the allowed increase of the'
print *,' square error relative to the misfit'
print *,' of the reference synthetics'
call gremhelp_subtit('tests')
print *,' orth calculate scalar products and normalized scalar'
print *,' products of vectors of partial derivatives'
......
c this is <res_opt.f>
c------------------------------------------------------------------------------
c $Id: res_opt.f,v 1.3 2000-06-02 20:31:20 thof Exp $
c $Id: res_opt.f,v 1.4 2002-04-19 15:20:36 forbrig Exp $
c
c 08/04/98 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -14,11 +14,14 @@ c 17/04/98 V1.1 normalize stabilization factor
c 02/06/00 V2.0 introduced mweightcondition
c changed algorithm completely
c V2.1 take real part of lq_dssd everywhere!
c 19/04/02 V3.0 reduced to core calculations
c
c==============================================================================
cS
c
logical function res_opt(mindex, nu, finalrms, x2ref)
logical function res_opt(mindex, nu, deltapar, deltaerror,
& deltamisfit, deltanofiterror,
& parselect)
c
c rate the total resolution for a model parameter by changing this
c parameter and optimizing with the rest of free parameters keeping
......@@ -26,13 +29,24 @@ c the resulting relative rms error increase at finalrms
c
c lq_dssd should be ready calculated (using inv_part, inv_dss and inv_dssd)
c
c mindex: index of model parameter to check
c nu: stabilization factor
c finalrms: rms-error to find parameter changes for
c x2ref: X2 error value to be calculated for reference model synthetics
c input:
c mindex: index of model parameter to check
c nu: stabilization factor
c
c output:
c deltapar: variation of parameter mindex
c deltaerror: change of error square sum due to change of parameter
c mindex and optimization of all other parameters
c deltamisfit: change of misfit due to (no damping term) du to change
c of parameter mindex and optimization of all other
c parameters
c deltanofiterror: change of error square sum only due to change of
c parameter mindex (no check for trade-off)
c parselect: flags indicating whether a parameter took part in
c optimization
c
integer mindex
real nu, finalrms, x2ref
real nu, deltapar, deltaerror, deltamisfit, deltanofiterror
c
include 'glq_dim.inc'
include 'glq_para.inc'
......@@ -43,20 +57,14 @@ c
c
cE
integer im,jm,imm,jmm,nmm, info, nrhs, ldb
integer msec, mpar, mpol
logical parselect(glqm_mano)
character*20 parname
character*1 uplo
logical result
double precision rmserror, nofiterror, mweightmax, errstep
real normnu, modelfactor
double precision mweightmax
real normnu, errstep, penaltystep
c
if (verb_subaction) print *,'ENTER res_opt(',mindex,',',nu,
& ',',finalrms,')'
if (verb_subaction) print *,'ENTER res_opt(',mindex,',',nu,')'
c
if (verb_subaction) print *,
& 'NOTICE (res_opt): rate resolution for parameter ',mindex,
& ' and nu ',nu,' and final rms ',finalrms
result=.true.
c find maximum search range value
mweightmax=0.d0
......@@ -64,6 +72,12 @@ c find maximum search range value
mweightmax=max(mweightmax,mweight(im))
enddo
c find selected parameters
c
c we set up a reduced system of linear equations. all parameters not matching
c the condition below will be omitted (although they are selected for
c inversion). This is necessary due to the fact, that only full sections (all
c polynomial orders) may be activated for inversion. If we wanr to check and
c optimize only for the mean value, we have to exclude the other orders.
do im=1,mod_n
if ((mweight(im)*mweightcondition).gt.mweightmax) then
parselect(im)=.true.
......@@ -74,6 +88,22 @@ c find selected parameters
c normalize stabilization factor
normnu=nu/float(mod_n)
c calculate dssd+nuww and right hand side
c
c this sets up the system of linear equations
c
c sum_(n != K) R_ln dm_n = - R_lK Dm_K
c
c where dm_n are the parameter optimizations, Dm_K is a specified change in
c the tested parameter (K=mindex), which is one search range here. And
c
c R_ln = Re( M_ln )
c
c where
c
c M = D*S*S*D + normnu * W*W
c
c where D are the partial derivatives, S is the matrix of data tolerance
c weights and W is the matrix of reziprocal search ranges.
imm=0
do im=1,mod_n
if ((parselect(im)).and.(im.ne.mindex)) then
......@@ -94,6 +124,8 @@ c calculate "right hand side"
enddo
endif
enddo
c we changed the test-paremeter (mindex) by one search range:
deltapar=mweight(mindex)
c remember reduced system size
nmm=imm
c
......@@ -129,84 +161,60 @@ c
c calculate linear rms error estimate for that model parameter change
c order of IF/ELSEIF matters!
c a parameter may be mindex and not be selected - however even if mindex is
c selected will will not be contained in the reduced model vector lq_mestim
rmserror=0.d0
c selected it will not be contained in the reduced model vector lq_mestim
c
deltaerror=0.
deltamisift=0.
deltanofiterror=0.
imm=0
do im=1,mod_n
jmm=0
if (im.eq.mindex) then
do jm=1,mod_n
if (jm.eq.mindex) then
errstep=mweight(im)*mweight(jm)*real(lq_dssd(im,jm))
rmserror=rmserror+errstep
nofiterror=errstep
penaltystep=normnu
deltaerror=deltaerror+errstep+penaltystep
deltanofiterror=deltanofiterror+errstep+penaltystep
deltamisfit=deltamisfit+errstep
elseif (parselect(jm)) then
jmm=jmm+1
errstep=mweight(im)*lq_mestim(jmm)*real(lq_dssd(im,jm))
rmserror=rmserror+errstep
deltaerror=deltaerror+errstep
deltanofiterror=deltanofiterror+errstep
deltamisfit=deltamisfit+errstep
endif
enddo
elseif (parselect(im)) then
imm=imm+1
jmm=0
do jm=1,mod_n
if (jm.eq.mindex) then
errstep=lq_mestim(imm)*mweight(mindex)*real(lq_dssd(im,jm))
rmserror=rmserror+errstep
deltaerror=deltaerror+errstep
deltanofiterror=deltanofiterror+errstep
deltamisfit=deltamisfit+errstep
elseif (parselect(jm)) then
jmm=jmm+1
errstep=lq_mestim(imm)*lq_mestim(jmm)*real(lq_dssd(im,jm))
rmserror=rmserror+errstep
errstep=lq_mestim(imm)*lq_mestim(jmm)*
& real(lq_dssd(im,jm))
if (jmm.eq.imm) then
penaltystep=normnu*lq_mestim(imm)**2/(mweight(im)**2)
else
penaltystep=0.
endif
deltaerror=deltaerror+errstep+penaltystep
deltanofiterror=deltanofiterror+errstep+penaltystep
deltamisfit=deltamisfit+errstep
endif
enddo
endif
enddo
rmserror=sqrt(rmserror)
nofiterror=sqrt(nofiterror)
modelfactor=sqrt((finalrms+1.)**2-1.)*sqrt(x2ref)/rmserror
c
c expand solution vector
imm=0
do im=1,mod_n
if (im.eq.mindex) then
mdelta(im)=mweight(im)*modelfactor
else
if ((mweight(im)*mweightcondition).gt.mweightmax) then
imm=imm+1
mdelta(im)=lq_mestim(imm)*modelfactor
else
mdelta(im)=0.d0
endif
endif
enddo
c
c print result
call mod_identify(mindex, msec, mpol, mpar, parname)
print *,'RESULTS (res_opt):'
print *,' testing model parameter ',mindex
print 50,mpol-1,parname(1:index(parname,' ')),msec
print 51,mweight(mindex),sqrt(1.+(nofiterror**2)/x2ref)-1.
print 52,sqrt(1.+(rmserror**2)/x2ref)-1.
print 53,finalrms
do im=1,mod_n
call mod_identify(im, msec, mpol, mpar, parname)
print 54,im,mpol-1,parname,msec,mdelta(im)
enddo
endif
c
res_opt=result
c
if (verb_subaction) print *,'LEAVE res_opt (result=',result,')'
c
return
50 format(3x,'that is ord. ',i1,' of ',a,' in sec. ',i2)
51 format(3x,'changing just this by ',g9.2,
& ' leads to relative rms-error increase of ',g12.5)
52 format(3x,' optimizing others in addition ',
& ' leads to relative rms-error increase of ',g12.5)
53 format(3x,'parameter changes to by applied',
& ' for a relative rms-error increase of ',g12.5,':')
54 format(5x,i2,': ord. ',i1,' of ',a14,' in sec. ',i2,
& ' to be changed by ',g12.5)
end
c
c ----- END OF res_opt.f -----
c this is <res_optrms.f>
cS
c ----------------------------------------------------------------------------
c ($Id: res_optrms.f,v 1.1 2002-04-16 18:04:13 forbrig Exp $)
c ($Id: res_optrms.f,v 1.2 2002-04-19 15:20:36 forbrig Exp $)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
......@@ -11,5 +12,104 @@ c 16/04/2002 V1.0 Thomas Forbriger
c
c ============================================================================
c
logical function res_optrms(mindex, nu, finalrms, x2ref)
c
integer mindex
real nu, finalrms, x2ref
c
c input:
c nu: damping weight factor
c finalrms: the relative increase in rms-error allowed
c (relative to the misfit of the reference synthetics)
c x2ref: the misfit for the reference synthetics
c
c output:
c on output the mdelta array in the glq_model common block contains
c appropriate changes for the test results
c
c return .true. if executed correctly
c
c
include 'glq_dim.inc'
include 'glq_para.inc'
include 'glq_model.inc'
include 'glq_data.inc'
include 'glq_inv.inc'
include 'glq_verbose.inc'
c
c local variables
real deltapar, deltaerror, deltamisfit, deltanofiterror
real rmserror, nofiterror, modelfactor
integer im,jm,imm,jmm,msec,mpar,mpol
logical parselect(glqm_mano)
logical res_opt
character*20 parname
c
cE
c
if (verb_subaction) print *,'ENTER res_optrms(',mindex,',',nu,
& ',',finalrms,')'
c
if (verb_subaction) print *,
& 'NOTICE (res_optrms): rate resolution for parameter ',mindex,
& ' and nu ',nu,' and final rms ',finalrms
c calculate raw test
if (res_opt(mindex, nu, deltapar, deltaerror,
& deltamisfit, deltanofiterror,
& parselect)) then
c
result=.true.
c
rmserror=sqrt(deltaerror)
nofiterror=sqrt(deltanofiterror)
modelfactor=sqrt((finalrms+1.)**2-1.)*sqrt(x2ref)/rmserror
c
c expand solution vector
imm=0
do im=1,mod_n
if (im.eq.mindex) then
mdelta(im)=deltapar*modelfactor
else
if (parselect(im)) then
imm=imm+1
mdelta(im)=lq_mestim(imm)*modelfactor
else
mdelta(im)=0.d0
endif
endif
enddo
c
c print result
call mod_identify(mindex, msec, mpol, mpar, parname)
print *,'RESULTS (res_optrms):'
print *,' testing model parameter ',mindex
print 50,mpol-1,parname(1:index(parname,' ')),msec
print 51,deltapar,sqrt(1.+(nofiterror**2)/x2ref)-1.
print 52,sqrt(1.+(rmserror**2)/x2ref)-1.
print 53,finalrms
do im=1,mod_n
call mod_identify(im, msec, mpol, mpar, parname)
print 54,im,mpol-1,parname,msec,mdelta(im)
enddo
endif
else
result=.false.
endif
c
res_opt=result
c
if (verb_subaction) print *,'LEAVE res_optrms (result=',result,')'
c
return
50 format(3x,'that is ord. ',i1,' of ',a,' in sec. ',i2)
51 format(3x,'changing just this by ',g9.2,
& ' leads to relative rms-error increase of ',g12.5)
52 format(3x,' optimizing others in addition ',
& ' leads to relative rms-error increase of ',g12.5)
53 format(3x,'parameter changes to be applied',
& ' for a relative rms-error increase of ',g12.5,':')
54 format(5x,i2,': ord. ',i1,' of ',a14,' in sec. ',i2,
& ' to be changed by ',g12.5)
end
c
c ----- END OF res_optrms.f -----
c this is <res_optsqr.f>
c ----------------------------------------------------------------------------
c ($Id: res_optsqr.f,v 1.1 2002-04-16 18:04:13 forbrig Exp $)
c ($Id: res_optsqr.f,v 1.2 2002-04-19 15:20:37 forbrig Exp $)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
......@@ -11,5 +11,104 @@ c 16/04/2002 V1.0 Thomas Forbriger
c
c ============================================================================
c
logical function res_optsqr(mindex, nu, finaldeltaerr)
c
integer mindex
real nu, finaldeltaerr
c
c input:
c nu: damping weight factor
c finaldeltaerr: the relative increase in error suqare sum allowed
c (relative to the misfit of the reference synthetics)
c
c output:
c on output the mdelta array in the glq_model common block contains
c appropriate changes for the test results
c
c return .true. if executed correctly
c
c
include 'glq_dim.inc'
include 'glq_para.inc'
include 'glq_model.inc'
include 'glq_data.inc'
include 'glq_inv.inc'
include 'glq_verbose.inc'
c
c local variables
real deltapar, deltaerror, deltamisfit, deltanofiterror
real rmserror, nofiterror, modelfactor
integer im,jm,imm,jmm,msec,mpar,mpol
logical parselect(glqm_mano)
logical res_opt
character*20 parname
c
cE
c
if (verb_subaction) print *,'ENTER res_optsqr(',mindex,',',nu,
& ',',finaldeltaerr,')'
c
if (verb_subaction) print *,
& 'NOTICE (res_optsqr): rate resolution for parameter ',mindex,
& ' and nu ',nu,' and final error increase ',finaldeltaerr
c calculate raw test
if (res_opt(mindex, nu, deltapar, deltaerror,
& deltamisfit, deltanofiterror,
& parselect)) then
c
result=.true.
c
modelfactor=sqrt(x2ref*finaldeltaerr/deltaerror)
c
c expand solution vector
imm=0
do im=1,mod_n
if (im.eq.mindex) then
mdelta(im)=mweight(im)*modelfactor
else
if (parselect(im)) then
imm=imm+1
mdelta(im)=lq_mestim(imm)*modelfactor
else
mdelta(im)=0.d0
endif
endif
enddo
c
c print result
call mod_identify(mindex, msec, mpol, mpar, parname)
print *,'RESULTS (res_optsqr):'
print *,' testing model parameter ',mindex
print 50,mpol-1,parname(1:index(parname,' ')),msec
print 51,deltapar,deltanofiterror/x2ref
print 52,deltaerror/x2ref
print 52,deltamisfit/x2ref
print 53,finaldeltaerr
do im=1,mod_n
call mod_identify(im, msec, mpol, mpar, parname)
print 54,im,mpol-1,parname,msec,mdelta(im)
enddo
endif
else
result=.false.
endif
c
res_opt=result
c
if (verb_subaction) print *,'LEAVE res_optsqr (result=',result,')'
c
return
50 format(3x,'that is ord. ',i1,' of ',a,' in sec. ',i2)
51 format(3x,'changing just this by ',g9.2,
& ' leads to relative square-error increase of ',g12.5)
52 format(3x,' optimizing others in addition ',
& ' leads to relative square-error increase of ',g12.5)
52 format(3x,' optimizing others in addition ',
& ' leads to relative misfit increase of ',g12.5)
53 format(3x,'parameter changes to be applied',
& ' for a relative square-error increase of ',g12.5,':')
54 format(5x,i2,': ord. ',i1,' of ',a14,' in sec. ',i2,
& ' to be changed by ',g12.5)
end
c
c ----- END OF res_optsqr.f -----
......@@ -2,7 +2,7 @@ c this is <res_wresmod.f>
c------------------------------------------------------------------------------
cS
c ($Source: /home/tforb/svnbuild/cvssource/CVS/thof/src/green/gremlin1/libs/res_wresmod.f,v $)
c ($Id: res_wresmod.f,v 1.3 2002-04-16 18:04:13 forbrig Exp $)
c ($Id: res_wresmod.f,v 1.4 2002-04-19 15:20:37 forbrig Exp $)
c
c 24/05/2000 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -15,19 +15,23 @@ c 16/04/2002 V1.2 now calls res_optrms
c
c==============================================================================
c
subroutine res_wresmod(mimin, mimax, nu, finalrms, basename)
subroutine res_wresmod(mimin, mimax, nu, finalrms, basename,
& dormstest)
c
c call res_optrms for given partial derivatives and precalculated
c matrices for model parameters mimin to mimax and write resulting
c ensemble of models - each to a seperate file
c
c nu: stabilization factor
c finalrms: rms error increase to find parameter changes for
c basename: base for filenames
c nu: stabilization factor
c finalrms: rms error increase to find parameter changes for
c basename: base for filenames
c dormstest: checks for rms error if true (else finalrms gives the
c relative increase in the square-error)
c
integer mimin,mimax
real nu, finalrms
character*(*) basename
logical dormstest
c
include 'glq_dim.inc'
include 'glq_model.inc'
......@@ -36,18 +40,21 @@ c
cE
c declare local variables
logical res_optrms
logical res_optsqr
logical result
integer mi, i
character*120 parname, title, filename
integer mpol, mpar, msec