[cig-commits] r16222 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: matlab src
carltape at geodynamics.org
carltape at geodynamics.org
Wed Feb 3 05:38:46 PST 2010
Author: carltape
Date: 2010-02-03 05:38:45 -0800 (Wed, 03 Feb 2010)
New Revision: 16222
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub4.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_variables.f90
Log:
Modified procedure for combining gradients of different parameter classes (structure, source location, source origin time).
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2010-02-03 01:58:11 UTC (rev 16221)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2010-02-03 13:38:45 UTC (rev 16222)
@@ -43,8 +43,6 @@
%parms = [500 0 3]; % joint inversion: 25 sources, Nfac=3
%parms = [600 0 3]; % structure inversion: 25 sources, Nfac=2
-parms = [1000 0 3];
-
%---------------------------------------------------------
irun0 = parms(1); % irun for m0
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-02-03 01:58:11 UTC (rev 16221)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-02-03 13:38:45 UTC (rev 16222)
@@ -50,8 +50,8 @@
!------------------------
! source and receiver variables (and fake receivers for spectral plots)
- integer :: istf, nsrc, nrec, nrecf, nevent, nevent_dat, nevent_run, nparm_src_run, idata
- double precision :: dnparm_src_run
+ integer :: istf, nsrc, nrec, nrecf, nevent, nevent_dat, nevent_run, idata
+ double precision :: dnevent_run
!integer, dimension(:), allocatable :: sglob, sglob_dat
!integer, dimension(:), allocatable :: sglob_temp, rglob_temp, fglob_temp
@@ -148,19 +148,19 @@
double precision :: cov_fac, cov_src_fac, cov_str_fac
double precision :: afac_target, coverage_str, coverage_src
double precision :: joint_str_src_grad_ratio, jsw_str, jsw_src, joint_total_weight, jfac
- double precision :: fac_str, fac_ts, fac_xs, fac_ys, fac_total
+ double precision :: fac_str, fac_ts, fac_xs, fac_ys, fac_total_m, fac_total_g
double precision :: ugrad_str, ugrad_ts, ugrad_xs, ugrad_ys, ugrad_src
- double precision, dimension(:), allocatable :: cov_imetric, icov_metric, cov_model, icov_model
+ double precision, dimension(:), allocatable :: cov_imetric, cov_model ! icov_metric, icov_model
! conjugate gradient optimization (using line-search)
double precision, dimension(:), allocatable :: pk, gk, g0, p0, m0, gt, mt, mdiff, mprior
- double precision :: beta_val, lam_0_val, lam_t_val, chi_t_val, chi_k_val, mu_val
+ double precision :: beta_val, beta_val_top, beta_val_bot, lam_0_val, lam_t_val, lam_t_val_bot, chi_t_val, chi_k_val, mu_val
double precision :: Pa,Pb,Pc,Pd,qfac,xx1,xx2,yy1,yy2,g1,g2,xmin
! arrays describing the structure and source parameters
- double precision, dimension(:), allocatable :: m0_vec, mt_vec, mtarget ! m0_vec_initial
- double precision, dimension(:), allocatable :: m_src_syn_vec, m_src_syn ! m_src_syn_vec_initial, m_src_syn_initial
- double precision, dimension(:), allocatable :: m_src_dat_vec, m_src_dat, m_src_prior
+ ! OBSOLETE: m_src_syn_vec_initial, m_src_syn_initial, m_src_syn_vec, m_src_dat_vec, m0_vec_initial,
+ double precision, dimension(:), allocatable :: m0_vec, mt_vec, mtarget
+ double precision, dimension(:), allocatable :: m_src_syn, m_src_dat, m_src_prior
double precision :: vel_min, vel_max
integer :: itest, nmod, nmod_src, nmod_str
@@ -191,7 +191,7 @@
integer :: itemp,itemp1,itemp2,itemp3, hmod
integer :: i, j, k, iq, itime, iglob, irec, isrc, ipick, ievent, icomp ! irunz
- integer :: isolver, irun0, irun, iopt, ispec, istep, istep_switch, imod, ipar
+ integer :: isolver, irun0, irun, iopt, ispec, istep, istep_switch, imod, ipar, imnorm
!double precision :: meas_pert, rand1, meas, ppert
@@ -448,12 +448,8 @@
ievent_min = 5 ; ievent_max = ievent_min
nevent_run = ievent_max - ievent_min + 1
- ! number of source parameters
- ! NOTE: we used to normalize the whole source portion of the model vector (3 x 25 = 75),
- ! but now we allow for variable balancing among structure, ts, xy, and ys.
- !nparm_src_run = nevent_run * NVAR_SOURCE
- nparm_src_run = nevent_run
- dnparm_src_run = dble(nparm_src_run)
+ ! number of sources in inversion
+ dnevent_run = dble(nevent_run)
! MAX source perturbations for location (lat-lon) and origin time
! dt0 = dt0_dat - dt0_syn ( < 0 for delayed synthetics)
@@ -1071,14 +1067,16 @@
! all vectors should be this length
nevent = MAX_EVENT
- ! convert from lat-lon to mesh coordinates (in meters) -- also update nevent
+ ! convert from lat-lon to mesh coordinates (in meters)
call mesh_geo(nevent,x_eve_lon0_syn,z_eve_lat0_syn,x_eve0_syn,z_eve0_syn,UTM_PROJECTION_ZONE,ILONLAT2MESH)
endif
! filter target points (utm-mesh) (update nevent)
print *
+ print *, 'nevent: ', nevent
print *, 'events for synthetics:'
call station_filter(nevent, x_eve0_syn, z_eve0_syn, ifilter_eve_syn, SOURCE_GRID_BUFFER)
+ print *, 'nevent: ', nevent
if(nevent < 1) stop 'Must have at least one event'
@@ -1091,7 +1089,7 @@
! origin time for synthetics
! initialize to mprior
- otime_syn(1:nevent) = otime(1:nevent) ! otime_syn is length MAX_EVENT, otime is length nevent
+ otime_syn(1:nevent) = otime(1:nevent) ! otime_syn is length nevent, otime is length MAX_EVENT
if (PERT_SOURCE_T == 1) then
otime_syn(:) = otime_syn(:) + otime_eve_syn_pert(:)
endif
@@ -1720,8 +1718,9 @@
close(18)
! allocate model vector
- allocate(m0(nmod),mt(nmod),mdiff(nmod),mprior(nmod))
- allocate(cov_imetric(nmod),icov_metric(nmod),cov_model(nmod),icov_model(nmod))
+ allocate(m0(nmod),mt(nmod),mdiff(nmod),mprior(nmod),mtarget(nmod))
+ allocate(cov_imetric(nmod),cov_model(nmod))
+ !allocate(icov_metric(nmod),icov_model(nmod))
allocate(gradient(nmod),gradient_data(nmod),gradient_model(nmod))
allocate(m0_vec(nmod),mt_vec(nmod))
allocate(g0(nmod),gt(nmod),gk(nmod),p0(nmod),pk(nmod))
@@ -1733,20 +1732,28 @@
!allocate(m_vel(nmod_str),m0_vel(nmod_str),mt_vel(nmod_str))
+ ! initialize model vectors
m0(:) = 0.0
mt(:) = 0.0
m0_vec(:) = 0.0
mt_vec(:) = 0.0
- mprior(:) = 0.0
+ mprior(:) = 0.0
+ mtarget(:) = 0.0
+ ! initialize gradient vectors
+ g0(:) = 0.0
+ gk(:) = 0.0
+ gt(:) = 0.0
+ p0(:) = 0.0
+ pk(:) = 0.0
+
!------------------
! source-related vectors
! scaling for source model coefficients
allocate(source_gradient(nmod_src))
- allocate(m_src_syn_vec(nmod_src),m_src_syn(nmod_src))
- allocate(m_src_dat_vec(nmod_src),m_src_dat(nmod_src))
- allocate(m_src_prior(nmod_src))
+ allocate(m_src_syn(nmod_src),m_src_dat(nmod_src),m_src_prior(nmod_src))
+ !allocate(m_src_syn_vec(nmod_src),m_src_dat_vec(nmod_src))
!allocate(m_src_syn_vec_initial(nmod_src),m_src_syn_initial(nmod_src))
!allocate(m_scale_src_all(nmod_src))
!allocate(m_src(nmod_src),mt_src(nmod_src),m0_src(nmod_src),m_scale_src_all(nmod_src))
@@ -1754,8 +1761,8 @@
m_src_prior(:) = 0.0
m_src_dat(:) = 0.0
m_src_syn(:) = 0.0
- m_src_dat_vec(:) = 0.0
- m_src_syn_vec(:) = 0.0
+ !m_src_dat_vec(:) = 0.0
+ !m_src_syn_vec(:) = 0.0
!m_src_syn_vec_initial(:) = 0.0
! fill source model vector for SYNTHETICS AND DATA
@@ -1774,13 +1781,13 @@
itemp2 = index_source(2,ievent)
itemp3 = index_source(3,ievent)
- m_src_syn_vec(itemp1) = otime_syn(ievent)
- m_src_syn_vec(itemp2) = x_eve_syn(ievent)
- m_src_syn_vec(itemp3) = z_eve_syn(ievent)
+ m_src_syn(itemp1) = otime_syn(ievent)
+ m_src_syn(itemp2) = x_eve_syn(ievent)
+ m_src_syn(itemp3) = z_eve_syn(ievent)
- m_src_dat_vec(itemp1) = otime_dat(ievent)
- m_src_dat_vec(itemp2) = x_eve_dat(ievent)
- m_src_dat_vec(itemp3) = z_eve_dat(ievent)
+ m_src_dat(itemp1) = otime_dat(ievent)
+ m_src_dat(itemp2) = x_eve_dat(ievent)
+ m_src_dat(itemp3) = z_eve_dat(ievent)
m_src_prior(itemp1) = otime(ievent)
m_src_prior(itemp2) = x_eve0_syn(ievent)
@@ -1798,54 +1805,41 @@
! perturbations are w.r.t. INITIAL SYNTHETICS
!m_src_dat(:) = m_src_dat_vec(:) - m_src_syn_vec_initial(:)
!m_src_syn(:) = m_src_syn_vec(:) - m_src_syn_vec_initial(:)
- m_src_dat(:) = m_src_dat_vec(:)
- m_src_syn(:) = m_src_syn_vec(:)
! create m0 -- should be identical to what is "un-done" in the CG algorithm
temp_local1(:,:,:) = log( beta_syn(:,:,:) / beta0 )
call local2mvec(temp_local1, nmod_src, m_src_syn, nmod, m0)
! create m0_vec
- call local2mvec(beta_syn, nmod_src, m_src_syn_vec, nmod, m0_vec)
+ call local2mvec(beta_syn, nmod_src, m_src_syn, nmod, m0_vec)
!m0_vec_initial(:) = m0_vec(:)
! create mprior
temp_local1(:,:,:) = 0.0
call local2mvec(temp_local1, nmod_src, m_src_prior, nmod, mprior)
- ! write out prior model
- open(88,file=trim(out_dir2)//'prior_model_vector.dat',status='unknown')
- do i = 1,nmod
- write(88,'(1e24.12)') mprior(i)
- enddo
- close(88)
+ ! create mtarget
+ temp_local1(:,:,:) = log( beta_dat(:,:,:) / beta0 )
+ call local2mvec(temp_local1, nmod_src, m_src_dat, nmod, mtarget)
- ! write out initial model vector
- if( .not. READ_IN ) then
- open(88,file=trim(out_dir2)//'initial_model_vector.dat',status='unknown')
+ ! write out initial source vector
+ !if( .not. READ_IN ) then
+ open(89,file=trim(out_dir2)//'initial_structure.dat',status='unknown')
do i = 1,nmod
- write(88,'(2e24.12)') m0_vec(i), m0(i)
+ write(89,'(3e24.12)') mprior(i), m0(i), mtarget(i)
enddo
- close(88)
- endif
+ close(89)
- ! write out initial source vector
- if( .not. READ_IN ) then
open(88,file=trim(out_dir2)//'initial_source.dat',status='unknown')
do i = 1,nmod_src
!write(88,'(5e24.12)') m_scale_src_all(i), m_src_syn_vec_initial(i), m_src_syn_initial(i), m_src_dat_vec(i), m_src_dat(i)
!write(88,'(4e24.12)') m_src_syn_vec_initial(i), m_src_syn_initial(i), m_src_dat_vec(i), m_src_dat(i)
- write(88,'(4e24.12)') m_src_syn_vec(i), m_src_syn(i), m_src_dat_vec(i), m_src_dat(i)
+ !write(88,'(4e24.12)') m_src_syn_vec(i), m_src_syn(i), m_src_dat_vec(i), m_src_dat(i)
+ write(88,'(3e24.12)') m_src_prior(i), m_src_syn(i), m_src_dat(i)
enddo
close(88)
+ !endif
- open(88,file=trim(out_dir2)//'prior_source.dat',status='unknown')
- do i = 1,nmod_src
- write(88,'(1e24.12)') m_src_prior(i)
- enddo
- close(88)
- endif
-
!------------------
! MODEL COVARIANCE MATRIX = inverse metric tensor
! NOTE: This part of the code is crucial for doing joint source-structure inversions,
@@ -1911,36 +1905,42 @@
!------------------------------
- ! DEFAULT: do not incorporate the norms of the initial gradients into the
- ugrad_str = 1.0 ; ugrad_ts = 1.0 ; ugrad_xs = 1.0 ; ugrad_ys = 1.0
+ ! initialize to no weights
+ covm_weight_parts(:) = 0.0
- ! DEFAULT: equal weight to each part
- fac_str = 1.0 ; fac_ts = 1.0 ; fac_xs = 1.0 ; fac_ys = 1.0
- fac_total = 1.0
-
- ! (1) unbalanced initial gradient values (gradient_norm_all_unbalanced.dat) -- 3 x 3 checker model
- ! (2) how much toward the norm should each variable contribute
+ ! factors for balancing the model norm term
+ ! NOTE: we ignore parts of the model norm that do not participate in the inversion
if ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0 ) then ! structure
fac_str = 1.00
-
elseif ( INV_STRUCT_BETA == 0 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 0 ) then ! origin time
fac_ts = 1.00
-
elseif ( INV_STRUCT_BETA == 0 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 1 ) then ! location
fac_xs = 0.50 ; fac_ys = 0.50
-
elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 0 ) then ! structure + origin time
fac_str = 0.50 ; fac_ts = 0.50
-
elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 1 ) then ! structure + location
fac_str = 0.50 ; fac_xs = 0.25 ; fac_ys = 0.25
-
elseif ( INV_STRUCT_BETA == 0 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! source
- !ugrad_ts = 6.538943677 ; ugrad_xs = 3.288688153 ; ugrad_ys = 3.900641302 ! fac_total = 43.0
- !fac_ts = 0.50 ; fac_xs = 0.25 ; fac_ys = 0.25
- !fac_total = (ugrad_ts/fac_ts) + (ugrad_xs/fac_xs) + (ugrad_ys/fac_ys)
+ fac_ts = 0.50 ; fac_xs = 0.25 ; fac_ys = 0.25
+ elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! joint
+ fac_str = 1.0/2.0 ; fac_ts = 1.0/6.0 ; fac_xs = 1.0/6.0 ; fac_ys = 1.0/6.0
+ else
+ stop 'you must invert for something'
+ endif
- elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! joint
+ ! NOTE: normalize to ensure that the total norm is 1
+ covm_weight_parts(1) = fac_str
+ covm_weight_parts(2) = fac_ts
+ covm_weight_parts(3) = fac_xs
+ covm_weight_parts(4) = fac_ys
+ covm_weight_parts(:) = covm_weight_parts(:) / sum(covm_weight_parts)
+ !covm_weight_parts(:) = 1.0 ! OLD VERSION
+
+ ! unbalanced initial gradient values (gradient_norm_all_unbalanced.dat)
+ covg_weight_parts(:) = 1.0
+
+ fac_total_g = 1.0
+ if ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! joint
! NOTE: we obtain these values by looking at the initial unbalanced data gradient (gradient_norm_data_all_unbalanced.dat)
! TEST CASES: 1-25 events, 1-5 events, 5-5 event
!ugrad_str = 0.1296496361d6 ; ugrad_ts = 0.2706999550d4; ugrad_xs = 0.1334372694d4 ; ugrad_ys = 0.1743549898d4 ! Nfac = 2, 25 events
@@ -1950,23 +1950,25 @@
!ugrad_str = 0.1205146534d6 ; ugrad_ts = 0.2926409454d3 ; ugrad_xs = 0.9695606936d3 ; ugrad_ys = 0.7224889563d3 ! Nfac = 3, 1 event
! 1300 - with 0.7,0.1,0.1,0.1 and old source
- !ugrad_str = 0.1191920326d6 ; ugrad_ts = 0.4533353832d3 ; ugrad_xs = 0.6191223030d2 ; ugrad_ys = 0.4387672586d2 ! Gaussians, 1 event
+ ugrad_str = 0.1191920326d6 ; ugrad_ts = 0.4533353832d3 ; ugrad_xs = 0.6191223030d2 ; ugrad_ys = 0.4387672586d2 ! Gaussians, 1 event
! 3200/5000
- ugrad_str = 0.3026988450d6 ; ugrad_ts = 0.8260031407d3 ; ugrad_xs = 0.1187469038d4 ; ugrad_ys = 0.1381078706d4 ! Gaussians, 25 events
+ !ugrad_str = 0.3026988450d6 ; ugrad_ts = 0.8260031407d3 ; ugrad_xs = 0.1187469038d4 ; ugrad_ys = 0.1381078706d4 ! Gaussians, 25 events
! 3300/1700
!ugrad_str = 0.3352965529d6 ; ugrad_ts = 0.8068082292d3 ; ugrad_xs = 0.1209790150d4 ; ugrad_ys = 0.1391331978d4 ! Gaussians, 25 events
! ad hoc: choose balance among the four parts of the gradient
- fac_str = 1.0/2.0 ; fac_ts = 1.0/6.0 ; fac_xs = 1.0/6.0 ; fac_ys = 1.0/6.0
+ !fac_str = 1.0/2.0 ; fac_ts = 1.0/6.0 ; fac_xs = 1.0/6.0 ; fac_ys = 1.0/6.0
!fac_str = 1.0/4.0 ; fac_ts = 1.0/4.0 ; fac_xs = 1.0/4.0 ; fac_ys = 1.0/4.0
- !fac_str = 0.7 ; fac_ts = 0.1 ; fac_xs = 0.1 ; fac_ys = 0.1
+ fac_str = 0.7 ; fac_ts = 0.1 ; fac_xs = 0.1 ; fac_ys = 0.1
!fac_str = 0.85 ; fac_ts = 0.05 ; fac_xs = 0.05 ; fac_ys = 0.05
- fac_total = (ugrad_str/fac_str) + (ugrad_ts/fac_ts) + (ugrad_xs/fac_xs) + (ugrad_ys/fac_ys)
+
+ covg_weight_parts(1) = ugrad_str/fac_str
+ covg_weight_parts(2) = ugrad_ts/fac_ts
+ covg_weight_parts(3) = ugrad_xs/fac_xs
+ covg_weight_parts(4) = ugrad_ys/fac_ys
+ covg_weight_parts(:) = covg_weight_parts(:) / sum(covg_weight_parts)
endif
- ! weight for source, to compare with ugrad_str (hence factor of 2)
- !ugrad_src = 2.0 * (fac_ts*ugrad_ts + fac_xs*ugrad_xs + fac_ys*ugrad_ys)
-
! Re-set any weights that are 0.0 to 1.0; these portions of the covariance matrix
! should not play a role, since the corresponding gradients will always be 0.0.
!if(fac_str < EPS) fac_str = 1.0
@@ -1988,28 +1990,21 @@
coverage_src = 1.0
cov_model(m_inds(1,1):m_inds(1,2)) = ( sigma_beta )**2 / da_local_vec(:) * AREA * coverage_str
- cov_model(m_inds(2,1):m_inds(2,2)) = sigma_ts**2 * dnparm_src_run * coverage_src
- cov_model(m_inds(3,1):m_inds(3,2)) = sigma_xs**2 * dnparm_src_run * coverage_src
- cov_model(m_inds(4,1):m_inds(4,2)) = sigma_zs**2 * dnparm_src_run * coverage_src
+ cov_model(m_inds(2,1):m_inds(2,2)) = sigma_ts**2 * dnevent_run * coverage_src
+ cov_model(m_inds(3,1):m_inds(3,2)) = sigma_xs**2 * dnevent_run * coverage_src
+ cov_model(m_inds(4,1):m_inds(4,2)) = sigma_zs**2 * dnevent_run * coverage_src
! incorporate relative weighting to make the final metric
! STRUCTURE: (fac_str / ugrad_str) * ugrad_str --> fac_str
- cov_imetric(m_inds(1,1):m_inds(1,2)) = cov_model(m_inds(1,1):m_inds(1,2)) * (fac_str / ugrad_str) * fac_total
- cov_imetric(m_inds(2,1):m_inds(2,2)) = cov_model(m_inds(2,1):m_inds(2,2)) * (fac_ts / ugrad_ts) * fac_total
- cov_imetric(m_inds(3,1):m_inds(3,2)) = cov_model(m_inds(3,1):m_inds(3,2)) * (fac_xs / ugrad_xs) * fac_total
- cov_imetric(m_inds(4,1):m_inds(4,2)) = cov_model(m_inds(4,1):m_inds(4,2)) * (fac_ys / ugrad_ys) * fac_total
+ cov_imetric(m_inds(1,1):m_inds(1,2)) = cov_model(m_inds(1,1):m_inds(1,2)) / covg_weight_parts(1)
+ cov_imetric(m_inds(2,1):m_inds(2,2)) = cov_model(m_inds(2,1):m_inds(2,2)) / covg_weight_parts(2)
+ cov_imetric(m_inds(3,1):m_inds(3,2)) = cov_model(m_inds(3,1):m_inds(3,2)) / covg_weight_parts(3)
+ cov_imetric(m_inds(4,1):m_inds(4,2)) = cov_model(m_inds(4,1):m_inds(4,2)) / covg_weight_parts(4)
! METRIC TENSOR: inverse diagponal covariance matrix
- icov_metric = 1.0 / cov_imetric
- icov_model = 1.0 / cov_model ! unbalanced version
+ !icov_metric = 1.0 / cov_imetric
+ !icov_model = 1.0 / cov_model ! unbalanced version
- ! compute norm of the target model
- ! create mtarget -- should be identical to what is shown in the CG algorithm
- allocate(mtarget(nmod))
- mtarget(:) = 0.0
- temp_local1(:,:,:) = log( beta_dat(:,:,:) / beta0 )
- call local2mvec(temp_local1, nmod_src, m_src_dat, nmod, mtarget)
-
! possible stopping criteria based on the target model
! NOTE 1: this depends on what you pick as your initial model (prior mean, or a sample)
! NOTE 2: THIS IS NOT USED
@@ -2030,51 +2025,61 @@
close(19)
open(unit=19,file=trim(out_dir2)//'scaling_values_covm.dat',status='unknown')
- write(19,*) 'fac_str', fac_str
- write(19,*) 'fac_ts', fac_ts
- write(19,*) 'fac_xs', fac_xs
- write(19,*) 'fac_ys', fac_ys
- write(19,*) 'fac_total', fac_total
+ do i=1,NVAR
+ write(19,*) 'covm_weight_parts', covm_weight_parts(i)
+ enddo
+ write(19,*) 'sum(covm_weight_parts)', sum(covm_weight_parts)
+ do i=1,NVAR
+ write(19,*) 'covg_weight_parts', covg_weight_parts(i)
+ enddo
+ write(19,*) 'sum(covg_weight_parts)', sum(covg_weight_parts)
+
+!!$ write(19,*) 'fac_str', fac_str
+!!$ write(19,*) 'fac_ts', fac_ts
+!!$ write(19,*) 'fac_xs', fac_xs
+!!$ write(19,*) 'fac_ys', fac_ys
+!!$ write(19,*) 'fac_total', fac_total
write(19,*) 'ugrad_str', ugrad_str
write(19,*) 'ugrad_ts', ugrad_ts
write(19,*) 'ugrad_xs', ugrad_xs
write(19,*) 'ugrad_ys', ugrad_ys
- write(19,*) 'dnparm_src_run', dnparm_src_run
+ write(19,*) 'dnevent_run', dnevent_run
write(19,*) 'coverage_str', coverage_str
write(19,*) 'coverage_src', coverage_src
- write(19,*) 'cov_imetric_fac_str', (fac_str / ugrad_str) * fac_total
- write(19,*) 'cov_imetric_fac_ts', (fac_ts / ugrad_ts) * fac_total
- write(19,*) 'cov_imetric_fac_xs', (fac_xs / ugrad_xs) * fac_total
- write(19,*) 'cov_imetric_fac_ys', (fac_ys / ugrad_ys) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_str', (fac_str / ugrad_str) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_ts', (fac_ts / ugrad_ts) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_xs', (fac_xs / ugrad_xs) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_ys', (fac_ys / ugrad_ys) * fac_total
close(19)
open(unit=19,file=trim(out_dir2)//'scaling_values.dat',status='unknown')
write(19,*) 'GJI_PAPER = ', GJI_PAPER
write(19,*) 'IRUNZ = ', IRUNZ
write(19,'(6i10)') 1, NLOCAL, NLOCAL+1, nmod_str, nmod_str+1, nmod
- write(19,'(2i10)') nparm_src_run, NLOCAL
+ write(19,'(2i10)') nevent_run, NLOCAL
write(19,'(5e14.6)') sum(da_local(:,:,:)), sum(da_local_vec(:)), sum(da_global(:)), LENGTH*HEIGHT, AREA
write(19,'(3e14.6)') minval(da_local_vec(:)), maxval(da_local_vec(:)), sum(da_local_vec(:))/NLOCAL
write(19,'(3e14.6)') minval(da_global(:)), maxval(da_global(:)), sum(da_global(:))/NGLOB
- write(19,*)
- write(19,*) 'COVERAGE WEIGHTS: '
- write(19,'(2f14.6)') coverage_str, coverage_src
- write(19,*) 'JOINT WEIGHTS (based on unbalanced gradients): '
- write(19,'(4f14.6)') fac_str, fac_ts, fac_xs, fac_ys
- write(19,'(1f14.6)') fac_total
- write(19,'(4e14.6)') ugrad_str, ugrad_ts, ugrad_xs, ugrad_ys
- !write(19,'(1e14.6)') ugrad_src
- write(19,'(4e14.6)') fac_str/ugrad_str, fac_ts/ugrad_ts, fac_xs/ugrad_xs, fac_ys/ugrad_ys
- write(19,'(4e14.6)') (fac_str/ugrad_str)*fac_total, (fac_ts/ugrad_ts)*fac_total, &
- (fac_xs/ugrad_xs)*fac_total, (fac_ys/ugrad_ys)*fac_total
- !write(19,'(a20,5e14.6)') 'JOINT WEIGHTS: ', joint_str_src_grad_ratio, joint_total_weight, jfac, jsw_str, jsw_src
- !write(19,'(a20,3e14.6)') 'STRUCTURE: ', jsw_str, joint_total_weight, jsw_str/joint_total_weight
- !write(19,'(a20,3e14.6)') 'SOURCE: ', jsw_src, joint_total_weight, jsw_src/joint_total_weight
+!!$ write(19,*)
+!!$ write(19,*) 'COVERAGE WEIGHTS: '
+!!$ write(19,'(2f14.6)') coverage_str, coverage_src
+!!$ write(19,*) 'JOINT WEIGHTS (based on unbalanced gradients): '
+!!$ write(19,'(4f14.6)') fac_str, fac_ts, fac_xs, fac_ys
+!!$ write(19,'(1f14.6)') fac_total
+!!$ write(19,'(4e14.6)') ugrad_str, ugrad_ts, ugrad_xs, ugrad_ys
+!!$ !write(19,'(1e14.6)') ugrad_src
+!!$ write(19,'(4e14.6)') fac_str/ugrad_str, fac_ts/ugrad_ts, fac_xs/ugrad_xs, fac_ys/ugrad_ys
+!!$ write(19,'(4e14.6)') (fac_str/ugrad_str)*fac_total, (fac_ts/ugrad_ts)*fac_total, &
+!!$ (fac_xs/ugrad_xs)*fac_total, (fac_ys/ugrad_ys)*fac_total
+!!$ !write(19,'(a20,5e14.6)') 'JOINT WEIGHTS: ', joint_str_src_grad_ratio, joint_total_weight, jfac, jsw_str, jsw_src
+!!$ !write(19,'(a20,3e14.6)') 'STRUCTURE: ', jsw_str, joint_total_weight, jsw_str/joint_total_weight
+!!$ !write(19,'(a20,3e14.6)') 'SOURCE: ', jsw_src, joint_total_weight, jsw_src/joint_total_weight
close(19)
open(unit=19,file=trim(out_dir2)//'cov_model_diagonal.dat',status='unknown')
do i = 1,nmod
- write(19,'(4e20.10)') cov_imetric(i), icov_metric(i), cov_model(i), icov_model(i)
+ write(19,'(4e20.10)') cov_model(i), cov_imetric(i)
+ !write(19,'(4e20.10)') cov_imetric(i), icov_metric(i), cov_model(i), icov_model(i)
enddo
close(19)
@@ -2212,13 +2217,13 @@
endif
! KEY: obtain the structure and source arrays from the model vector:
- ! structure parameters fill the top portion (bulk_syn, beta_syn);
- ! source parameters fill the bottom portion (m_src_syn_vec)
+ ! structure parameters fill the top portion (beta_syn);
+ ! source parameters fill the bottom portion (m_src_syn)
if(itest==0) then ! reference model (current model)
- call mvec2local(nmod, nmod_src, m0_vec, beta_syn, m_src_syn_vec)
+ call mvec2local(nmod, nmod_src, m0_vec, beta_syn, m_src_syn)
elseif(itest==1) then ! test model
- call mvec2local(nmod, nmod_src, mt_vec, beta_syn, m_src_syn_vec)
+ call mvec2local(nmod, nmod_src, mt_vec, beta_syn, m_src_syn)
endif
! for CG algorithm, update kappa_syn and mu_syn
@@ -2234,7 +2239,7 @@
write(filename,'(a,i4.4,a)') trim(out_dir2)//'src_syn_m',hmod,'.dat'
open(unit=18,file=filename,status='unknown')
- m_src_syn_vec(:) = m_src_dat_vec(:) ! initialize to no perturbation from target sources
+ m_src_syn(:) = m_src_dat(:) ! initialize to no perturbation from target sources
do i = 1,nevent
!itemp1 = (i-1)*3 + 1 ! origin time
!itemp2 = (i-1)*3 + 2 ! x position
@@ -2245,14 +2250,14 @@
itemp3 = index_source(3,i)
read(18,*) temp1, temp2, &
- m_src_syn_vec(itemp1), m_src_syn_vec(itemp2), m_src_syn_vec(itemp3), &
+ m_src_syn(itemp1), m_src_syn(itemp2), m_src_syn(itemp3), &
temp3, temp4, temp5
enddo
close(18)
! update the entries in the model vector, then re-compute the norm
!m0(nmod_str+1 : nmod) = m_src_syn_vec(:) - m_src_syn_vec_initial(:)
- m0(nmod_str+1 : nmod) = m_src_syn_vec(:)
+ m0(nmod_str+1 : nmod) = m_src_syn(:)
endif
@@ -2262,60 +2267,91 @@
! NOTE 2: THE FACTOR OF 0.5 IS NOT INCLUDED HERE
! NOTE 3: This is norm-squared (no square root is taken)
- ! model vector, unbalanced Cm (cov_model)
- if(itest==0) then ! reference model (current model)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, m0, cov_model, &
- model_norm, model_norm_struct, model_norm_source, model_norm_parts, mprior)
- elseif(itest==1) then ! test model
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, mt, cov_model, &
- model_norm, model_norm_struct, model_norm_source, model_norm_parts, mprior)
- endif
- filename = trim(out_dir1)//'model_norm_all_unbalanced.dat'
- call write_norm_sq(filename,model_norm_parts)
+ imnorm = 1 ! =1 for m-mprior ; =0 for g
- ! model vector, balanced Cm (cov_imetric)
+ ! model vector, cov_model
if(itest==0) then ! reference model (current model)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, m0, cov_imetric, &
- model_norm, model_norm_struct, model_norm_source, model_norm_parts, mprior)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ m0, mprior, cov_model, model_norm_parts, covm_weight_parts)
elseif(itest==1) then ! test model
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, mt, cov_imetric, &
- model_norm, model_norm_struct, model_norm_source, model_norm_parts, mprior)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ mt, mprior, cov_model, model_norm_parts, covm_weight_parts)
endif
- filename = trim(out_dir1)//'model_norm_all.dat'
+ filename = trim(out_dir1)//'model_norm_all_cov_model.dat'
call write_norm_sq(filename,model_norm_parts)
- ! target model vector, unbalanced Cm (cov_model)
- filename = trim(out_dir1)//'model_norm_target_all_unbalanced.dat'
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, mtarget, cov_model, &
- model_target_norm, model_target_norm_struct, model_target_norm_source, model_target_norm_parts, mprior)
- call write_norm_sq(filename,model_target_norm_parts)
+!!$ ! model vector, cov_imetric
+!!$ if(itest==0) then ! reference model (current model)
+!!$ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, m0, mprior, cov_imetric, &
+!!$ model_norm_parts, covm_weight_parts)
+!!$ elseif(itest==1) then ! test model
+!!$ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, mt, mprior, cov_imetric, &
+!!$ model_norm_parts, covm_weight_parts)
+!!$ endif
+!!$ filename = trim(out_dir1)//'model_norm_all_cov_imetric.dat'
+!!$ call write_norm_sq(filename,model_norm_parts)
- ! target model vector, balanced Cm (cov_imetric)
- filename = trim(out_dir1)//'model_norm_target_all.dat'
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, mtarget, cov_imetric, &
- model_target_norm, model_target_norm_struct, model_target_norm_source, model_target_norm_parts, mprior)
- call write_norm_sq(filename,model_target_norm_parts)
+ ! model_norm is used in the misfit function
+ model_norm_struct = model_norm_parts(1)
+ model_norm_source = sum( model_norm_parts(2:4) )
+ model_norm = sum( model_norm_parts(1:4) )
+!!$ ! write CG vectors to file
+!!$ if(1==1) then
+!!$ open(unit=19,file=trim(out_dir1)//'model_vectors_structure.dat',status='unknown')
+!!$ do i = 1,NLOCAL
+!!$ write(19,'(4e16.6)') mprior(i), m0(i), mt(i), mtarget(i)
+!!$ enddo
+!!$ close(19)
+!!$ open(unit=19,file=trim(out_dir1)//'model_vectors_source.dat',status='unknown')
+!!$ do i = NLOCAL+1,nmod
+!!$ !write(19,'(4e16.6)') mprior(i), m0(i), mt(i), mtarget(i)
+!!$ write(19,'(3e16.6)') mprior(i), cov_model(i), mtarget(i)
+!!$ enddo
+!!$ close(19)
+!!$ endif
+
+ ! target model vector, cov_model
+ ! NOTE: THIS SHOULD NOT CHANGE
+ filename = trim(out_dir1)//'model_norm_target_all_cov_model.dat'
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ mtarget, mprior, cov_model, model_norm_target_parts, covm_weight_parts)
+ call write_norm_sq(filename,model_norm_target_parts)
+
+!!$ ! target model vector, cov_imetric
+!!$ filename = trim(out_dir1)//'model_norm_target_all_cov_imetric.dat'
+!!$ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, mtarget, mprior, cov_imetric, &
+!!$ model_norm_target_parts, covm_weight_parts)
+!!$ call write_norm_sq(filename,model_norm_target_parts)
+
+ model_norm_target_struct = model_norm_target_parts(1)
+ model_norm_target_source = sum( model_norm_target_parts(2:4) )
+ model_norm_target = sum( model_norm_target_parts(1:4) )
+
! compute the difference between the present model and the target model
mdiff(:) = 0.0
if(itest==0) then
- mdiff = mtarget - m0
+ mdiff = mtarget - m0 + mprior
elseif(itest==1) then
- mdiff = mtarget - mt
+ mdiff = mtarget - mt + mprior
endif
- ! diff vector, unbalanced Cm (cov_model)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, mdiff, cov_model, &
- model_diff_norm, model_diff_norm_struct, model_diff_norm_source, model_diff_norm_parts)
- filename = trim(out_dir1)//'model_norm_diff_all_unbalanced.dat'
- call write_norm_sq(filename,model_diff_norm_parts)
+ ! diff vector, cov_model
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ mdiff, mprior, cov_model, model_norm_diff_parts, covm_weight_parts)
+ filename = trim(out_dir1)//'model_norm_diff_all_cov_model.dat'
+ call write_norm_sq(filename,model_norm_diff_parts)
- ! diff vector, balanced Cm (cov_imetric)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, mdiff, cov_imetric, &
- model_diff_norm, model_diff_norm_struct, model_diff_norm_source, model_diff_norm_parts)
- filename = trim(out_dir1)//'model_norm_diff_all.dat'
- call write_norm_sq(filename,model_diff_norm_parts)
+!!$ ! diff vector, cov_imetric
+!!$ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, mdiff, mprior, cov_imetric, &
+!!$ model_norm_diff_parts, covm_weight_parts)
+!!$ filename = trim(out_dir1)//'model_norm_diff_all_cov_imetric.dat'
+!!$ call write_norm_sq(filename,model_norm_diff_parts)
+ model_norm_diff_struct = model_norm_diff_parts(1)
+ model_norm_diff_source = sum( model_norm_diff_parts(2:4) )
+ model_norm_diff = sum( model_norm_diff_parts(1:4) )
+
!stop 'TESTING NORMS'
!--------------------------
@@ -2619,9 +2655,9 @@
allocate(ispec_src(nsrc),xi_src(nsrc),gamma_src(nsrc))
! KEY: source values for current model and event ievent
- origin_time_syn = m_src_syn_vec(itemp1) ! new origin time
- x_src(1) = m_src_syn_vec(itemp2) ! new x position
- z_src(1) = m_src_syn_vec(itemp3) ! new z position
+ origin_time_syn = m_src_syn(itemp1) ! new origin time
+ x_src(1) = m_src_syn(itemp2) ! new x position
+ z_src(1) = m_src_syn(itemp3) ! new z position
! filter target points (utm-mesh) -- update nsrc
call station_filter(nsrc, x_src, z_src, ifilter_src, SOURCE_GRID_BUFFER)
@@ -2702,19 +2738,19 @@
write(91,*) 'istep, imod, ievent, irun : '
write(91,*) istep, imod, ievent, irun
write(91,*) 'SOURCE MODEL (xs, ys, t0) :'
- write(91,'(a12,3a18)') ' ','m_src_syn_vec','m_src_dat_vec','dat - syn'
+ write(91,'(a12,3a18)') ' ','m_src_syn','m_src_dat','dat - syn'
write(91,'(a12,3f18.8)') ' t0, s : ', &
- m_src_syn_vec(itemp1), &
- m_src_dat_vec(itemp1), &
- m_src_dat_vec(itemp1) - m_src_syn_vec(itemp1)
+ m_src_syn(itemp1), &
+ m_src_dat(itemp1), &
+ m_src_dat(itemp1) - m_src_syn(itemp1)
write(91,'(a12,3f18.8)') ' xs, km : ', &
- m_src_syn_vec(itemp2)/1000.0, &
- m_src_dat_vec(itemp2)/1000.0, &
- (m_src_dat_vec(itemp2) - m_src_syn_vec(itemp2))/1000.0
+ m_src_syn(itemp2)/1000.0, &
+ m_src_dat(itemp2)/1000.0, &
+ (m_src_dat(itemp2) - m_src_syn(itemp2))/1000.0
write(91,'(a12,3f18.8)') ' zs, km : ', &
- m_src_syn_vec(itemp3)/1000.0, &
- m_src_dat_vec(itemp3)/1000.0, &
- (m_src_dat_vec(itemp3) - m_src_syn_vec(itemp3))/1000.0
+ m_src_syn(itemp3)/1000.0, &
+ m_src_dat(itemp3)/1000.0, &
+ (m_src_dat(itemp3) - m_src_syn(itemp3))/1000.0
endif
!--------------------------------------
@@ -2883,17 +2919,17 @@
! sources for data
write(19,'(8e20.10)') xtemp, ztemp, &
- m_src_dat_vec(itemp1), m_src_dat_vec(itemp2), m_src_dat_vec(itemp3), &
- m_src_dat_vec(itemp1) - m_src_dat_vec(itemp1), &
- m_src_dat_vec(itemp2) - m_src_dat_vec(itemp2), &
- m_src_dat_vec(itemp3) - m_src_dat_vec(itemp3)
+ m_src_dat(itemp1), m_src_dat(itemp2), m_src_dat(itemp3), &
+ m_src_dat(itemp1) - m_src_dat(itemp1), &
+ m_src_dat(itemp2) - m_src_dat(itemp2), &
+ m_src_dat(itemp3) - m_src_dat(itemp3)
! sources for synthetics
write(20,'(8e20.10)') xtemp, ztemp, &
- m_src_syn_vec(itemp1), m_src_syn_vec(itemp2), m_src_syn_vec(itemp3), &
- m_src_dat_vec(itemp1) - m_src_syn_vec(itemp1), &
- m_src_dat_vec(itemp2) - m_src_syn_vec(itemp2), &
- m_src_dat_vec(itemp3) - m_src_syn_vec(itemp3)
+ m_src_syn(itemp1), m_src_syn(itemp2), m_src_syn(itemp3), &
+ m_src_dat(itemp1) - m_src_syn(itemp1), &
+ m_src_dat(itemp2) - m_src_syn(itemp2), &
+ m_src_dat(itemp3) - m_src_syn(itemp3)
enddo
close(19)
close(20)
@@ -3281,9 +3317,12 @@
if(ISOURCE_LOG) then
write(91,*) 'SUMMED MISFIT FOR ALL EVENTS'
write(91,'(a30,1i16)') ' N-data : ', nmeas_run
- write(91,'(a30,1e16.8)') ' chi-data(m) : ', data_norm
- write(91,'(a30,3e16.8)') ' chi-model(m) : ', model_norm, model_norm_struct, model_norm_source
- write(91,'(a30,3e16.8)') ' chi_model_target(m) : ', model_target_norm, model_target_norm_struct, model_target_norm_source
+ write(91,'(a30,1e16.8)') ' data_norm(m) : ', data_norm
+ write(91,'(a30,3e16.8)') ' model_norm(m) : ', model_norm, model_norm_struct, model_norm_source
+ write(91,'(a30,3e16.8)') ' model_norm_target(m) : ', &
+ model_norm_target, model_norm_target_struct, model_norm_target_source
+ write(91,'(a30,3e16.8)') ' model_norm_diff(m) : ', &
+ model_norm_diff, model_norm_diff_struct, model_norm_diff_source
!write(91,'(a30,1e16.8)') ' chi_model_stop(mtarget) : ', chi_model_stop
write(91,'(a30,1e16.8)') ' chi_data_stop : ', chi_data_stop
write(91,'(a30,1e16.8)') ' chi(m) : ', chi_val
@@ -3328,7 +3367,6 @@
m0(:) = 0.0
temp_local1(:,:,:) = log( beta_syn(:,:,:) / beta0 )
- m_src_syn(:) = m_src_syn_vec(:)
call local2mvec(temp_local1, nmod_src, m_src_syn, nmod, m0)
!!$ ! m_src_syn: w.r.t initial source values (m_src_syn_vec_initial)
@@ -3392,11 +3430,12 @@
call local2mvec(temp_local1, nmod_src, source_gradient, nmod, gradient_data)
! add the gradient term associated with the MODEL NORM in the misfit function
+ ! KEY QUESTION: cov_model or cov_imetric?
if( INCLUDE_MODEL_NORM ) then
if (itest == 0) then
- gradient_model(:) = (m0(:) - mprior(:)) / cov_imetric(:)
+ gradient_model(:) = (m0(:) - mprior(:)) / cov_model(:)
else
- gradient_model(:) = (mt(:) - mprior(:)) / cov_imetric(:)
+ gradient_model(:) = (mt(:) - mprior(:)) / cov_model(:)
endif
endif
@@ -3425,46 +3464,64 @@
enddo
close(19)
- ! compute the overall gradient norm
- ! NOTE THE TRICK: To compute the norm of the gradient, we would MULTIPLY by Cm, rather than divide.
- ! We achieve this by passing in 1/Cm instead of Cm.
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient, icov_metric, &
- gradient_norm, gradient_norm_struct, gradient_norm_source, gradient_norm_parts)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient_data, icov_metric, &
- gradient_data_norm, gradient_data_norm_struct, gradient_data_norm_source, gradient_data_norm_parts)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient_model, icov_metric, &
- gradient_model_norm, gradient_model_norm_struct, gradient_model_norm_source, gradient_model_norm_parts)
+ ! NORMS OF THE GRADIENT (C^-1 is used, not C ; mprior is not used)
+ imnorm = 0
- ! old scaling used in GJI-2007 paper
- !scale_struct_gji = sqrt( norm_source_2 / (norm_bulk_2 + norm_beta_2) )
- scale_struct_gji = sqrt( gradient_norm_source / gradient_norm_struct )
- open(unit=19,file=trim(out_dir1)//'scaling_gji.dat',status='unknown')
- write(19,'(1e20.10)') scale_struct_gji
- close(19)
+ ! compute the balanced gradient norm
+ if (0==1) then
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient, mprior, cov_imetric, gradient_norm_parts)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient_data, mprior, cov_imetric, gradient_norm_data_parts)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient_model, mprior, cov_imetric, gradient_norm_model_parts)
+ else
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient, mprior, cov_model, gradient_norm_parts, covg_weight_parts)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient_data, mprior, cov_model, gradient_norm_data_parts, covg_weight_parts)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient_model, mprior, cov_model, gradient_norm_model_parts, covg_weight_parts)
+ endif
+! gradient_norm, gradient_norm_struct, gradient_norm_source
+! gradient_norm_data, gradient_norm_data_struct, gradient_norm_data_source
+! gradient_norm_model, gradient_norm_model_struct, gradient_norm_model_source
+
+!!$ ! old scaling used in GJI-2007 paper
+!!$ !scale_struct_gji = sqrt( norm_source_2 / (norm_bulk_2 + norm_beta_2) )
+!!$ scale_struct_gji = sqrt( gradient_norm_source / gradient_norm_struct )
+!!$ open(unit=19,file=trim(out_dir1)//'scaling_gji.dat',status='unknown')
+!!$ write(19,'(1e20.10)') scale_struct_gji
+!!$ close(19)
+
! write gradient norms to file
filename1 = trim(out_dir1)//'gradient_norm_all.dat'
filename2 = trim(out_dir1)//'gradient_norm_data_all.dat'
filename3 = trim(out_dir1)//'gradient_norm_model_all.dat'
call write_norm_sq(filename1,gradient_norm_parts)
- call write_norm_sq(filename2,gradient_data_norm_parts)
- call write_norm_sq(filename3,gradient_model_norm_parts)
+ call write_norm_sq(filename2,gradient_norm_data_parts)
+ call write_norm_sq(filename3,gradient_norm_model_parts)
- ! SAME AS ABOVE, but use the unbalanced icov_model instead of icov_metric
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient, icov_model, &
- gradient_norm, gradient_norm_struct, gradient_norm_source, gradient_norm_parts)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient_data, icov_model, &
- gradient_data_norm, gradient_data_norm_struct, gradient_data_norm_source, gradient_data_norm_parts)
- call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient_model, icov_model, &
- gradient_model_norm, gradient_model_norm_struct, gradient_model_norm_source, gradient_model_norm_parts)
+ ! SAME AS ABOVE, but use cov_model instead of cov_imetric
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient, mprior, cov_model, gradient_norm_parts)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient_data, mprior, cov_model, gradient_norm_data_parts)
+ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ gradient_model, mprior, cov_model, gradient_norm_model_parts)
+! gradient_norm, gradient_norm_struct, gradient_norm_source
+! gradient_norm_data, gradient_norm_data_struct, gradient_norm_data_source
+! gradient_norm_model, gradient_norm_model_struct, gradient_norm_model_source
+
! write (unbalanced) gradient norms to file
filename1 = trim(out_dir1)//'gradient_norm_all_unbalanced.dat'
filename2 = trim(out_dir1)//'gradient_norm_data_all_unbalanced.dat'
filename3 = trim(out_dir1)//'gradient_norm_model_all_unbalanced.dat'
call write_norm_sq(filename1,gradient_norm_parts)
- call write_norm_sq(filename2,gradient_data_norm_parts)
- call write_norm_sq(filename3,gradient_model_norm_parts)
+ call write_norm_sq(filename2,gradient_norm_data_parts)
+ call write_norm_sq(filename3,gradient_norm_model_parts)
!!$ ! scaling for joint inversions
!!$ if( istep==0 .and. 0==1 ) then
@@ -3479,7 +3536,7 @@
!!$ ! NORM-SQUARED of the steepest ascent vector : sqrt( ghat C ghat )
!!$ norm_bulk_2 = sum( kbulk_smooth(:,:,:)**2 * da_local(:,:,:) * (m_scale_str(1))**2 * AREA )
!!$ norm_beta_2 = sum( kbeta_smooth(:,:,:)**2 * da_local(:,:,:) * (m_scale_str(2))**2 * AREA )
-!!$ norm_source_2 = sum( ( source_gradient(:) * m_scale_src_all(:) )**2 * dnparm_src_run ) ! m_scale_src_all NOT ALLOCATED
+!!$ norm_source_2 = sum( ( source_gradient(:) * m_scale_src_all(:) )**2 * dnevent_run * NVAR_SOURCE ) ! m_scale_src_all NOT ALLOCATED
!!$
!!$ ! This shows the balance used in the GJI-2007 paper
!!$ scale_struct_gji = sqrt( norm_source_2 / (norm_bulk_2 + norm_beta_2) )
@@ -3498,7 +3555,7 @@
!!$ write(19,*) 'GJI_PAPER = ', GJI_PAPER
!!$ write(19,*) 'IRUNZ = ', IRUNZ
!!$ write(19,'(6i14)') 1, NLOCAL, NLOCAL+1, nmod_str, nmod_str+1, nmod
-!!$ write(19,'(2i14)') nparm_src_run, NLOCAL
+!!$ write(19,'(2i14)') nevent_run, NLOCAL
!!$ write(19,'(5e14.6)') sum(da_local(:,:,:)), sum(da_local_vec(:)), sum(da_global(:)), LENGTH*HEIGHT, AREA
!!$ write(19,'(3e14.6)') minval(da_local_vec(:)), maxval(da_local_vec(:)), sum(da_local_vec(:))/NLOCAL
!!$ write(19,'(3e14.6)') minval(da_global(:)), maxval(da_global(:)), sum(da_global(:))/NGLOB
@@ -3539,7 +3596,9 @@
beta_val = 0.0
p0(:) = 0.0 ! initialize
else
- beta_val = sum((gk(:) - g0(:)) * (cov_imetric(:)*gk(:)) ) / sum(g0(:) * (cov_imetric(:)*g0(:)) )
+ beta_val_top = sum((gk(:) - g0(:)) * (cov_imetric(:)*gk(:)) )
+ beta_val_bot = sum(g0(:) * (cov_imetric(:)*g0(:)) )
+ beta_val = beta_val_top / beta_val_bot
endif
pk(:) = -cov_imetric(:) * gk(:) + beta_val * p0(:)
@@ -3563,7 +3622,8 @@
mu_val = chi_data_stop
!mu_val = 0.0
- lam_t_val = 2.0*(mu_val - chi_k_val) / sum( gk(:) * pk(:) ) ! ghat dot p
+ lam_t_val_bot = sum( gk(:) * pk(:) ) ! ghat dot p
+ lam_t_val = 2.0*(mu_val - chi_k_val) / lam_t_val_bot
! alternative approach (more ad hoc)
!lam_t_val = -2.0*(1.0 - mu_val)*chi_k_val / sum( gk(:) * pk(:) )
@@ -3795,14 +3855,15 @@
! conjugate-gradient variables
deallocate(m0,mt,mdiff,mtarget,mprior)
- deallocate(cov_imetric,icov_metric,cov_model,icov_model)
+ deallocate(cov_imetric,cov_model)
+ !deallocate(icov_metric,icov_model)
deallocate(gradient,gradient_data,gradient_model)
deallocate(gradient_data_all)
deallocate(g0,gt,gk,p0,pk)
deallocate(source_gradient)
deallocate(m0_vec,mt_vec)
- deallocate(m_src_syn_vec,m_src_syn)
- deallocate(m_src_dat_vec,m_src_dat)
+ deallocate(m_src_syn,m_src_dat)
+ !deallocate(m_src_syn_vec,m_src_dat_vec)
!deallocate(m0_vec_initial,m_src_syn_vec_initial,m_src_syn_initial)
!deallocate(m_scale_src_all)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 2010-02-03 01:58:11 UTC (rev 16221)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 2010-02-03 13:38:45 UTC (rev 16222)
@@ -99,6 +99,7 @@
! reference model and target model choice
integer, parameter :: IMODEL_SYN = 3
integer, parameter :: IMODEL_DAT = 3
+ logical, parameter :: M0ISMPRIOR = .false.
!-----------------------------------------------------------------------------------------
! 0 1 2 3
!-----------------------------------------------------------------------------------------
@@ -111,8 +112,8 @@
! smooth the structure model or kernels
! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
@@ -202,12 +203,12 @@
! what to perturb, what to invert
! (For the inverse tests, we only allow perturbations in beta.)
integer, parameter :: PERT_STRUCT_BETA = 1
- integer, parameter :: PERT_SOURCE_T = 0
- integer, parameter :: PERT_SOURCE_X = 0
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
integer, parameter :: INV_STRUCT_BETA = 1
- integer, parameter :: INV_SOURCE_T = 0
- integer, parameter :: INV_SOURCE_X = 0
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90 2010-02-03 01:58:11 UTC (rev 16221)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90 2010-02-03 13:38:45 UTC (rev 16222)
@@ -1120,58 +1120,71 @@
!----------------------------------------------------
- subroutine compute_norm_sq(ievent_min, ievent_max, nevent, &
- index_source, nmod, mvec, cov_model, &
- norm_total, norm_struct, norm_source, norm_parts, &
- mvec_prior)
+ subroutine compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ mvec, mvec_prior, cov_model, norm_parts, norm_parts_weight)
! This computes the norm-squared of a model vector using the model covariance.
! The dimensions of the input model vector are always the same, but this norm
! takes into account the number of events used, which may be less than nevent.
- ! It outputs the overall norm, the source part, and the structure part.
! UPDATE: If mvec_prior is present, then subtract from mvec: (m-mprior)^T Cm (m-mprior)
+ ! NOTE 1: mprior is only used is imnorm = 1
- integer, intent(in) :: ievent_min, ievent_max, nevent, nmod
+ integer, intent(in) :: imnorm, ievent_min, ievent_max, nevent, nmod
integer, dimension(NVAR_SOURCE, nevent), intent(in) ::index_source
- double precision, dimension(nmod), intent(in) :: mvec, cov_model
- double precision, dimension(nmod), intent(in), optional :: mvec_prior
+ double precision, dimension(nmod), intent(in) :: mvec, mvec_prior, cov_model
+ double precision, dimension(NVAR), intent(in), optional :: norm_parts_weight
!double precision, intent(out) :: norm_total, norm_struct, norm_source
- double precision, intent(out) :: norm_total, norm_struct, norm_source
double precision, dimension(NVAR), intent(out) :: norm_parts
- double precision, dimension(nmod) :: mtemp
+ double precision, dimension(nmod) :: mtemp, ctemp
+ double precision, dimension(NVAR) :: npw
integer :: ievent, itemp1, itemp2, itemp3
!----------
-
- norm_parts(:) = 0.0
- mtemp = mvec
- if (present(mvec_prior)) mtemp = mvec - mvec_prior
+ norm_parts(:) = 0.0
+ ctemp(:) = 0.0
+ mtemp(:) = 0.0
+ npw(:) = 1.0
+ ! NOTE 1: if taking the norm of a gradient, use the inverse covariance matrix
+ ! NOTE 2: if taking the norm of a model, use m - mprior
+ if (imnorm == 0) then
+ ctemp(:) = 1.0 / cov_model(:)
+ mtemp(:) = mvec(:)
+ if (present(norm_parts_weight)) npw(:) = 1.0 / norm_parts_weight(:)
+
+ elseif (imnorm == 1) then
+ ctemp(:) = cov_model(:)
+ mtemp(:) = mvec(:) - mvec_prior(:)
+ if (present(norm_parts_weight)) npw(:) = norm_parts_weight(:)
+ else
+ stop 'imnorm must = 0 or 1'
+ endif
+
! structure part of the norm -- BETA only
- norm_parts(1) = sum( mtemp(1 : NLOCAL)**2 / cov_model(1 : NLOCAL) )
+ ! NOTE: division corresponds to inversion of a diagonal covariance matrix
+ norm_parts(1) = sum( mtemp(1 : NLOCAL)**2 / ctemp(1 : NLOCAL) )
! source part of the norm -- only events that you are inverting for
do ievent = ievent_min, ievent_max
- !itemp1 = 2*NLOCAL + (ievent-1)*3 + 1
- !itemp2 = 2*NLOCAL + (ievent-1)*3 + 2
- !itemp3 = 2*NLOCAL + (ievent-1)*3 + 3
-
itemp1 = NLOCAL + index_source(1,ievent)
itemp2 = NLOCAL + index_source(2,ievent)
itemp3 = NLOCAL + index_source(3,ievent)
- norm_parts(2) = norm_parts(2) + mtemp(itemp1)**2 / cov_model(itemp1)
- norm_parts(3) = norm_parts(3) + mtemp(itemp2)**2 / cov_model(itemp2)
- norm_parts(4) = norm_parts(4) + mtemp(itemp3)**2 / cov_model(itemp3)
+ norm_parts(2) = norm_parts(2) + mtemp(itemp1)**2 / ctemp(itemp1)
+ norm_parts(3) = norm_parts(3) + mtemp(itemp2)**2 / ctemp(itemp2)
+ norm_parts(4) = norm_parts(4) + mtemp(itemp3)**2 / ctemp(itemp3)
enddo
- norm_struct = norm_parts(1)
- norm_source = sum( norm_parts(2:4) )
- norm_total = sum( norm_parts(1:4) )
+ !norm_struct = norm_parts(1)
+ !norm_source = sum( norm_parts(2:4) )
+ !norm_total = sum( norm_parts(1:4) )
+ ! weight each part of the norm
+ norm_parts(:) = norm_parts(:) * npw(:)
+
end subroutine compute_norm_sq
!-----------------------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub4.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub4.f90 2010-02-03 01:58:11 UTC (rev 16221)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub4.f90 2010-02-03 13:38:45 UTC (rev 16222)
@@ -380,12 +380,15 @@
! normalization factors
! NOTE sign convention for N (differs from 2005 GJI paper) so that Nnorm >= 0
Nnorm = -DT * sum( syn_displ(:) * syn_accel(:) )
-
Mnorm = DT * sum( syn_displ(:) * syn_displ(:) )
- ! cross-correlation traveltime adjoint source
- ft_bar_t(:) = -syn_veloc(:) / Nnorm ! BD kernel
- ft_t(:) = -( tshift_xc_pert / cov_data(imeasure) ) * ft_bar_t(:) ! misfit kernel (note sign)
+ ! cross-correlation traveltime adjoint source for banana-doughnut kernel
+ ft_bar_t(:) = -syn_veloc(:) / Nnorm
+
+ ! cross-correlation traveltime adjoint source for misfit kernel
+ ! NOTE 1: sign convention
+ ! NOTE 2: weighted by measurement and (diagonal) data covariance matrix
+ ft_t(:) = -( tshift_xc_pert / cov_data(imeasure) ) * ft_bar_t(:)
! cross-correlation amplitude adjoint source
! You have TWO OPTIONS: measure the amplitudes based on DISPLACEMENT or VELOCITY
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_variables.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_variables.f90 2010-02-03 01:58:11 UTC (rev 16221)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_variables.f90 2010-02-03 13:38:45 UTC (rev 16222)
@@ -124,19 +124,20 @@
! misfit function
double precision :: chi_data(MAX_EVENT, MAX_SR, MAX_COMP, MAX_PHASE), chi_val, chi_val_0
double precision :: data_norm, chi_data_stop, chi_model_stop
- double precision :: model_target_norm, model_norm, model_diff_norm
- double precision :: model_target_norm_struct, model_norm_struct, model_diff_norm_struct
- double precision :: model_target_norm_source, model_norm_source, model_diff_norm_source
- double precision :: gradient_norm, gradient_model_norm, gradient_data_norm
- double precision :: gradient_norm_struct, gradient_model_norm_struct, gradient_data_norm_struct
- double precision :: gradient_norm_source, gradient_model_norm_source, gradient_data_norm_source
- double precision, dimension(NVAR) :: model_target_norm_parts, model_norm_parts, model_diff_norm_parts
- double precision, dimension(NVAR) :: gradient_norm_parts, gradient_model_norm_parts, gradient_data_norm_parts
+ double precision :: model_norm_target, model_norm, model_norm_diff
+ double precision :: model_norm_target_struct, model_norm_struct, model_norm_diff_struct
+ double precision :: model_norm_target_source, model_norm_source, model_norm_diff_source
+ double precision :: gradient_norm, gradient_norm_model, gradient_norm_data
+ double precision :: gradient_norm_struct, gradient_norm_model_struct, gradient_norm_data_struct
+ double precision :: gradient_norm_source, gradient_norm_model_source, gradient_norm_data_source
+ double precision, dimension(NVAR) :: covm_weight_parts, covg_weight_parts
+ double precision, dimension(NVAR) :: model_norm_target_parts, model_norm_parts, model_norm_diff_parts
+ double precision, dimension(NVAR) :: gradient_norm_parts, gradient_norm_model_parts, gradient_norm_data_parts
! double precision :: chi_model_norm, chi_model_norm_struct, chi_model_norm_source
! double precision :: chi_model_norm_target, chi_model_norm_target_struct, chi_model_norm_target_source
! double precision :: chi_gradient_norm, chi_gradient_norm_struct, chi_gradient_norm_source
-! double precision :: chi_gradient_data_norm, chi_gradient_data_norm_struct, chi_gradient_data_norm_source
-! double precision :: chi_gradient_model_norm, chi_gradient_model_norm_struct, chi_gradient_model_norm_source
+! double precision :: chi_gradient_norm_data, chi_gradient_norm_data_struct, chi_gradient_norm_data_source
+! double precision :: chi_gradient_norm_model, chi_gradient_norm_model_struct, chi_gradient_norm_model_source
double precision :: var_red_val
! kernels
More information about the CIG-COMMITS
mailing list