[cig-commits] r16202 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: . INPUT PLOTTING matlab src
carltape at geodynamics.org
carltape at geodynamics.org
Sun Jan 31 13:45:04 PST 2010
Author: carltape
Date: 2010-01-31 13:45:03 -0800 (Sun, 31 Jan 2010)
New Revision: 16202
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_initial.dat
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_sigmas.dat
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_target.dat
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_sigmas.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_run.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.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
Log:
Added an option for a third model in the synthetic 2D inversions: mprior. The most general inversion uses a prior model, and initial model, and a target model. The parameter M0ISMPRIOR controls whether the initial model is set to be the prior model.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_initial.dat
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_initial.dat (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_initial.dat 2010-01-31 21:45:03 UTC (rev 16202)
@@ -0,0 +1,25 @@
+ -3.536835037163e-01 -2.386685777941e+02 -1.924548663627e+03
+ -1.881576167538e-01 -1.042879043462e+02 -5.339549004841e+01
+ -2.277842775899e-02 4.890166921559e+02 1.486799867792e+03
+ 1.904860781068e-01 1.466610751464e+03 8.225874607935e+02
+ 2.977253191260e-01 1.463891872924e+03 -3.323701514803e+03
+ 1.760292362496e-01 6.509873245940e+02 -1.205978299711e+03
+ -2.354483040599e-01 -8.672916261062e+02 -2.633583938356e+02
+ -4.053812296059e-02 1.238918609589e+03 -1.274503992146e+03
+ 1.063281279690e-01 -6.524093722175e+02 -9.133714195961e+02
+ -1.881159312061e-01 -3.038375050267e+03 4.390025628386e+02
+ 1.590961434494e-01 5.320958961441e+02 -2.129866153020e+02
+ 5.486779578383e-03 -1.132837148733e+03 -3.192951815800e+03
+ 9.821927425374e-02 4.533744206902e+02 2.514603479788e+03
+ -1.325931037854e-01 -1.338958420396e+03 1.894042185975e+03
+ 2.205036528620e-01 -1.175893596315e+03 3.235857062970e+02
+ 2.891277920173e-01 5.709607091382e+02 -2.497226892390e+02
+ 1.737994592790e-01 -3.992399541996e+02 1.334166581139e+03
+ -2.282989794563e-01 2.027992537581e+02 -6.286970292329e+02
+ -1.104160300014e-01 3.932107415916e+03 -4.006905958110e+02
+ 4.753646700002e-02 1.792802574399e+02 -3.911767104531e+02
+ -1.260375289651e-01 -1.175519002512e+03 1.078917271940e+02
+ 1.493269447788e-01 -1.156491247969e+03 -8.815532967918e+02
+ 2.505222349116e-01 2.125261150226e+03 5.271405158113e+02
+ -3.683170932055e-01 -2.859674289660e+02 -7.793544591443e+02
+ 1.577051140744e-01 -1.009426030469e+03 -1.201704593132e+03
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_sigmas.dat
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_sigmas.dat (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_sigmas.dat 2010-01-31 21:45:03 UTC (rev 16202)
@@ -0,0 +1 @@
+ 2.000000000000e-01 1.400000000000e+03 1.400000000000e+03
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_target.dat
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_target.dat (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_txy_pert_target.dat 2010-01-31 21:45:03 UTC (rev 16202)
@@ -0,0 +1,25 @@
+ 1.160775060720e-01 -9.944076592739e+02 1.681483130184e+03
+ 2.295204481617e-01 9.887647914473e+01 1.524297747240e+03
+ 2.241758767645e-01 8.756730581052e+02 -2.728471278381e+03
+ -1.887103681581e-01 2.551573352352e+01 -2.300826221248e+02
+ 7.975610515441e-02 -1.599978278510e+03 -6.502537573573e+02
+ 2.334750818602e-02 8.696364180818e+02 -1.849726529968e+03
+ -3.367222291166e-01 7.062839556562e+02 1.903053424246e+03
+ 3.840391963698e-01 -6.328389730655e+02 1.390155644706e+03
+ -1.146043946295e-01 -2.023797889493e+01 -1.587390555107e+03
+ -3.018745708645e-01 1.660838687452e+03 2.682720430261e+03
+ -1.552197664068e-02 -2.460625018794e+02 -3.664534071485e+02
+ -4.296284196720e-02 -1.396073553711e+03 -9.247054599965e+02
+ 1.087083767281e-01 1.781633767059e+03 9.802272276703e+02
+ 3.220444508241e-01 -8.502726610638e+02 1.699657077216e+02
+ -2.714351891995e-01 4.314503525290e+02 -4.404834761421e+02
+ -9.758934299250e-03 1.031070499578e+03 1.853675270537e+03
+ 3.232490476556e-01 1.898577379972e+03 -8.337524624502e+02
+ -9.857980617259e-02 -1.092169238099e+03 -2.250918218853e+03
+ 1.391451317585e-01 1.265008843569e+03 2.388349722488e+02
+ -3.378974504494e-01 -6.134957109442e+02 -2.013929505267e+02
+ 9.849649966840e-02 -1.742447956667e+03 1.483855848018e+03
+ -1.534172803983e-01 1.251823429703e+03 -2.215510718741e+03
+ 2.079603848032e-01 2.274332387072e+02 -6.243670342737e+02
+ 9.135382207105e-03 4.130605389406e+03 5.108737665895e+02
+ 6.454447257253e-02 1.005642040108e+03 -1.074550681503e+03
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2010-01-31 21:45:03 UTC (rev 16202)
@@ -204,13 +204,29 @@
#$lam_lab = sprintf("\@~\154\@~ = %.0f km",$lam/1000);
#print "\n$per_lab, $c0_lab, $lam_lab \n";
-open(IN,"$dir0/reference_values.dat"); @lines = <IN>; ($alpha0,$beta0) = split(" ",$lines[0]);
-open(IN,"$dir0/reference_period.dat"); $per = <IN>; chomp($per);
+$fileR1 = "$dir0/reference_values.dat";
+$fileR2 = "$dir0/reference_period.dat";
+$fileR3 = "$dir0/sigma_values.dat";
+if (not -f $fileR1) {die("Check if fileR1 $fileR1 exist or not\n");}
+if (not -f $fileR2) {die("Check if fileR2 $fileR2 exist or not\n");}
+if (not -f $fileR3) {die("Check if fileR3 $fileR3 exist or not\n");}
+
+open(IN,$fileR1); @lines = <IN>; ($alpha0,$beta0) = split(" ",$lines[0]);
+open(IN,$fileR2); $per = <IN>; chomp($per);
$lam = $beta0*$per;
$per_lab = sprintf("T = %3.3f s",$per);
$beta0_lab = sprintf("beta0 = %3.3f km/s",$beta0/1000);
$lam_lab = sprintf("\@~\154\@~ = %.0f km",$lam/1000);
print "\n$per_lab, $beta0_lab, $lam_lab\n";
+
+open(IN,$fileR3); @lines = <IN>;
+($sigma_B,undef) = split(" ",$lines[0]);
+($sigma_ts,undef) = split(" ",$lines[1]);
+($sigma_xs,undef) = split(" ",$lines[2]);
+($sigma_zs,undef) = split(" ",$lines[3]);
+print "\n SIGMAS: $sigma_B, $sigma_ts, $sigma_xs, $sigma_zs\n";
+#die("TESTING");
+
#-------------------------
# write plotting scripts
@@ -530,7 +546,7 @@
$B = "$B0".$Bopts[1];
# make colorpoint file (-Z for continuous; -I to invert)
-$cmax = 1;
+$cmax = 2*$sigma_ts;
$dc = $cmax/40;
$T2 = "-T-$cmax/$cmax/$dc";
$cptfile = "color_src.cpt";
@@ -599,16 +615,16 @@
print CSH "psxy $dots_file $J $R $dot_info -C$cptfile -K -O -V >> $psfile \n";
if ($i==0) {
- # depth scale
+ # origin time scale
#$Dx = 0.25; $Dy = -0.2;
$Dx = 0.25; $Dy = 3.0;
-
- $Dscale_otime = "-D$Dx/$Dy/0.8/0.10h";;
- $Bscale_otime = "-B1f0.25:\"Origin time error (s)\": -E5p";
+ $Dscale_otime = "-D$Dx/$Dy/0.8/0.10h";
+ $Bscale_otime = sprintf("-B%.2ff0.1:\"Origin time error (s)\": -E8p",2*$sigma_ts);
print CSH "psscale -C$cptfile $Dscale_otime $Bscale_otime -K -O -V >> $psfile\n";
# convert source errors to dot size
- @misloc_dots = (0,2,4);
+ #$sigma_xs = 2.0;
+ @misloc_dots = (0,$sigma_xs/1000,2*$sigma_xs/1000);
$ndot = @misloc_dots;
$xlon0 = -117; $dlon = 1.0;
for ($j = 0; $j < $ndot; $j = $j+1) {
@@ -632,9 +648,9 @@
# source error scale -- text
print CSH "pstext $J $R -N -K -V -O $origin_box >> $psfile <<EOF
- $xlon[2] $yp2 $fsize2 0 $fontno CM 4
- $xlon[1] $yp2 $fsize2 0 $fontno CM 2
- $xlon[0] $yp2 $fsize2 0 $fontno CM 0
+ $xlon[2] $yp2 $fsize2 0 $fontno CM $misloc_dots[2]
+ $xlon[1] $yp2 $fsize2 0 $fontno CM $misloc_dots[1]
+ $xlon[0] $yp2 $fsize2 0 $fontno CM $misloc_dots[0]
$xlon[1] $yp1 $fsize2 0 $fontno CM Source
$xlon[1] $yp0 $fsize2 0 $fontno CM mislocation (km)
EOF\n";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2010-01-31 21:45:03 UTC (rev 16202)
@@ -76,13 +76,9 @@
Example 5: 25 sources, joint inversion using CG algorithm
Example 6: 25 sources, structure inversion using CG algorithm
Example 7: 25 sources, structure inversion using subspace method
-------------
-Using 2 Gaussian fields:
-Example 8: 1 source, structure inversion
-Example 9: 1 source, source inversion
-Example 10: 1 source, joint inversion
-Example 11: 5 sources, joint inversion
-Example 12: 25 sources, joint inversion
+----------
+2 Gaussian fields
+Example 8: 25 sources, joint inversion, new source perturbations
To get a sense for the output from these examples, open the composite PDF files
/PLOTTING/FIGURES/Example1_run_0100.pdf
@@ -100,6 +96,8 @@
OUTPUT/run_0000/chi.dat 100.1837219688
OUTPUT/run_0001/chi.dat 75.2034179848
OUTPUT/run_0002/chi.dat 31.7804791149
+ OUTPUT/run_0003/chi.dat 23.5724627945
+ OUTPUT/run_0004/chi.dat 10.7315686377
------------------
EXAMPLE 1 (run_0100)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m 2010-01-31 21:45:03 UTC (rev 16202)
@@ -1,13 +1,11 @@
%
% wave2d_cg_figs.m
-% CARL TAPE, 21-Nov-2008
-% printed xxx
+% Carl Tape, 21-Nov-2008
%
% This program is called following wave2d_cg_poly.m.
%
% calls xxx
% called by xxx
-%
%irun0 = 140;
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-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2010-01-31 21:45:03 UTC (rev 16202)
@@ -1,7 +1,6 @@
%
% wave2d_cg_poly.m
-% CARL TAPE, 21-Nov-2008
-% printed xxx
+% Carl Tape, 21-Nov-2008
%
% This program takes the output from wave2d.f90 and plots the polynomials
% used in the line-search conjugate gradient algorithm within the code.
@@ -17,7 +16,6 @@
% linefit.m
% axes_expand.m, fontsize.m, colors.m, wysiwyg.m
% called by xxx
-%
format long
format compact
@@ -45,6 +43,8 @@
%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/matlab/wave2d_cg_run.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_run.m 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_run.m 2010-01-31 21:45:03 UTC (rev 16202)
@@ -36,14 +36,14 @@
if ~exist(dirrand,'dir'), error(['dirrand ' dirrand ' does not exist']); end
% KEY
-parms = [1500 5];
+parms = [2100 0];
%parms = [1500 4];
% icg = 1: emulate wave2d.f90 algorithm
% icg = 2: emulate wave2d.f90 algorithm, but use 32 x 32 cells
% icg = 3: CG algorithm with full covariance matrix and 32 x 32 cells
icg = 1;
-iwrite = 1;
+iwrite = 0;
iplotmod = 1;
iplotker = 0;
sticg = sprintf('%2.2i',icg);
@@ -197,6 +197,9 @@
%write(88,'(2e24.12)') m0_vec(i), m0(i)
[m0_vec_initial, m0_initial] = textread([dir0 'initial_model_vector.dat'],'%f%f');
+% PRIOR model vector
+mprior = textread([dir0 'prior_model_vector.dat'],'%f');
+
% CURRENT model vector
%write(19,'(4e16.6)') m0(i), mt(i), m0_vec(i), mt_vec(i)
[m0, mt0, m0_vec, mt_vec0] = textread([idir1 'cg_model_vectors.dat'],'%f%f%f%f');
@@ -221,7 +224,8 @@
summed_ker = sum(kernels_all,2);
% source gradient
-source_gradient = load([idir1 'source_gradient.dat']);
+%source_gradient = load([idir1 'source_gradient.dat']);
+source_gradient = zeros(nmod_src,1);
% data covariance
% nmeas_run = nevent_run * nrec * NCOMP
@@ -280,7 +284,7 @@
% icg = 1 ==> emulate CG variables in wave2d.f90
if icg==1
- if 0==1
+ if 1==1
cov_model = zeros(nmod,1);
cov_model(indB) = sigma_beta^2 ./ da_local_vec * Atot * coverage_str;
cov_model(indT) = sigma_ts^2 * dnparm_src_run * coverage_src;
@@ -360,17 +364,19 @@
% check chi model values
model_norm_parts = zeros(4,1);
- model_norm_parts(1) = sum( m0(indB).^2 ./ cov_imetric(indB) );
- model_norm_parts(2) = sum( m0(indT).^2 ./ cov_imetric(indT) );
- model_norm_parts(3) = sum( m0(indX).^2 ./ cov_imetric(indX) );
- model_norm_parts(4) = sum( m0(indY).^2 ./ cov_imetric(indY) );
- model_norm = sum( m0.^2 ./ cov_imetric );
+ model_norm_parts(1) = sum( (m0(indB)-mprior(indB)).^2 ./ cov_imetric(indB) );
+ model_norm_parts(2) = sum( (m0(indT)-mprior(indT)).^2 ./ cov_imetric(indT) );
+ model_norm_parts(3) = sum( (m0(indX)-mprior(indX)).^2 ./ cov_imetric(indX) );
+ model_norm_parts(4) = sum( (m0(indY)-mprior(indY)).^2 ./ cov_imetric(indY) );
+ model_norm = sum( (m0-mprior).^2 ./ cov_imetric );
% check
model_norm0, model_norm, sum(model_norm_parts)
% check chi model values
chi_k_val = 0.5*(data_norm + model_norm)
+ break
+
%---------------------------------------------------------
% emulate CG algorithm in wave2d.f90 -- THIS ASSUMES A DIAGONAL COVARIANCE MATRIX
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_sigmas.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_sigmas.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_sigmas.m 2010-01-31 21:45:03 UTC (rev 16202)
@@ -0,0 +1,476 @@
+%
+% wave2d_sigmas.m
+% CARL TAPE, 20-Jan-2010
+% printed xxx
+%
+% Produces Gaussian errors using the Matlab function randn.
+% randn(1,N) generates an N x 1 column vector, the values having a
+% standard deviation of 1.
+%
+% calls plot_histo.m
+% called by xxx
+%
+
+format short
+format compact
+close all
+clear
+clc
+
+igauss_diff = 0;
+isigma_source = 1;
+isigma_struct = 0;
+isigma_data = 0;
+
+iwrite = 0;
+
+%odir = '/net/denali/home2/carltape/wave2d/2d_adjoint/INPUT/';
+odir = '/home/carltape/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work/INPUT/';
+
+% labels for source parameters
+xlabs = {'origin time (s)','x position (m)','y position (m)'};
+xlabs2 = {'ts (s)','xs (m)','ys (m)'};
+
+NLOCAL = 40000;
+NPARM_SOURCE = 3;
+
+%---------------------------------------------------------
+% demonstrate the analytical expression for a linear combination of two Gaussians
+
+if igauss_diff == 1
+
+ % input parameters
+ N = 1000;
+ sig1 = 2; mu1 = 4;
+ sig2 = 6; mu2 = -20;
+ A = 0.5; B = -3;
+
+ g1 = mu1 + sig1 * randn(N,1);
+ g2 = mu2 + sig2 * randn(N,1);
+ dmat = zeros(N,N);
+ for ii=1:N
+ for jj=1:N
+ % linear sum of every pair of samples
+ dmat(ii,jj) = A*g1(ii) + B*g2(jj);
+ end
+ end
+ g1g2 = dmat(:);
+
+ % KEY: analytical expressions for combined mu and sig
+ sigs = [sig1 sig2 sqrt(A^2*sig1^2 + B^2*sig2^2)];
+ mus = [mu1 mu2 A*mu1+B*mu2];
+
+ figure; nr=3; nc=1; ylims = [0 0.15];
+ subplot(nr,nc,1); sig = sigs(1); mu = mus(1);
+ Nbin = plot_histo(g1,mu + [-4*sig : sig/4 : 4*sig]);
+ grid on; xlabel(sprintf('g1: mean %.2f sigma %.2f',mu,sig)); ylim(ylims);
+
+ subplot(nr,nc,2); sig = sigs(2); mu = mus(2);
+ Nbin = plot_histo(g2,mu + [-4*sig : sig/4 : 4*sig]);
+ grid on; xlabel(sprintf('g2: mean %.2f sigma %.2f',mu,sig)); ylim(ylims);
+
+ subplot(nr,nc,3); sig = sigs(3); mu = mus(3);
+ Nbin = plot_histo(g1g2,mu + [-4*sig : sig/4 : 4*sig]);
+ grid on; xlabel(sprintf('A*g1 + B*g2: mean %.2f sigma %.2f (A = %.2f, B = %.2f)',mu,sig,A,B)); ylim(ylims);
+ orient tall, wysiwyg;
+
+end
+
+%---------------------------------------------------------
+
+if isigma_source == 1
+ ifig = 0;
+
+ %N = 10000; RUNS = 1;
+ N = 25; RUNS = 1000;
+
+ % KEY COMMANDS
+ %sigs = [0.5 2000 2000]; % sigmas to approximate the GJI2007 perturbations
+ sigs = [0.2 1400 1400];
+
+ % the STD for the difference (or sum) of two samples
+ sigssum = sqrt( 2*sigs.^2 )
+
+ cov_model = (N*NPARM_SOURCE) * repmat(sigs.^2,N,1);
+
+ % generate RUNS samples
+ ts_samples = sigs(1)*randn(N,RUNS);
+ xs_samples = sigs(2)*randn(N,RUNS);
+ ys_samples = sigs(3)*randn(N,RUNS);
+
+ % compute norms
+ norm2_ts = sum( ts_samples.^2 / (N*NPARM_SOURCE*sigs(1)^2), 1)';
+ norm2_xs = sum( xs_samples.^2 / (N*NPARM_SOURCE*sigs(2)^2), 1)';
+ norm2_ys = sum( ys_samples.^2 / (N*NPARM_SOURCE*sigs(3)^2), 1)';
+ norm2_all = [norm2_ts norm2_xs norm2_ys];
+ norm2_tot = sum(norm2_all,2);
+
+ figure; nr=2; nc=2;
+ ncen = 1/3;
+ for ii=1:3
+ subplot(nr,nc,ii);
+ Nbin = plot_histo(norm2_all(:,ii),[0:ncen/8:2*ncen]);
+ grid on; xlabel(['Norm of ' xlabs{ii} ' for each sample']);
+ ylim([0 0.3]);
+ end
+ ncen = 1; subplot(nr,nc,4);
+ Nbin = plot_histo(norm2_tot,[0:ncen/8:2*ncen]);
+ grid on; xlabel('Norm of each sample');
+ ylim([0 0.3]);
+ orient tall, wysiwyg;
+
+ % pick a subset of samples whose overall norms are close to 1.0 and
+ % whose and norm-parts are close to 0.33
+ i1 = find( abs( log(norm2_tot / 1.0 )) < 0.1 );
+ i2 = find( abs( log(norm2_ts / (1/3) )) < 0.1 );
+ i3 = find( abs( log(norm2_xs / (1/3) )) < 0.1 );
+ i4 = find( abs( log(norm2_ys / (1/3) )) < 0.1 );
+ iset = mintersect(i1,i2,i3,i4);
+ nset = length(iset)
+
+ if ifig==1
+ for kk=1:nset
+ jj = iset(kk);
+ Ymat = [ts_samples(:,jj) xs_samples(:,jj) ys_samples(:,jj)];
+ Nvec = [norm2_all(jj,:) norm2_tot(jj)];
+
+ figure; nr=3; nc=1;
+ for ii = 1:3
+ sig = sigs(ii);
+ yvec = Ymat(:,ii);
+ subplot(nr,nc,ii);
+ Nbin = plot_histo(yvec,[-5*sig : sig/4 : 5*sig]);
+ grid on; xlabel(['Perturbation in ' xlabs{ii}]); ylim([0 0.2]);
+ title(sprintf('%s -- sigma = %.3f -- Norm of %i parameters = %.4f',...
+ xlabs{ii},sig,N,Nvec(ii)));
+ end
+ orient tall, wysiwyg;
+ end
+ end
+
+ % from the set, pick two samples whose DIFFERENCE is not an outlier
+ % --> the DIFFERENCE represents the net perturbation needed
+ n = 0;
+ Imat = zeros(n,3);
+ Nmat = zeros(n,3);
+ for kk=1:nset
+ ik = iset(kk);
+ Kmat = [ts_samples(:,ik) xs_samples(:,ik) ys_samples(:,ik)];
+ for jj=kk+1:nset
+ ij = iset(jj);
+ Jmat = [ts_samples(:,ij) xs_samples(:,ij) ys_samples(:,ij)];
+ Rmat = Kmat - Jmat;
+
+ n = n+1;
+ % NOTE siggssum instead of sigs
+ Nmat(n,1) = sum( Rmat(:,1).^2 / (N*NPARM_SOURCE*sigssum(1)^2) );
+ Nmat(n,2) = sum( Rmat(:,2).^2 / (N*NPARM_SOURCE*sigssum(2)^2) );
+ Nmat(n,3) = sum( Rmat(:,3).^2 / (N*NPARM_SOURCE*sigssum(3)^2) );
+ Imat(n,:) = [n ik ij];
+ end
+ end
+ Nmat_tot = sum(Nmat, 2);
+ %figure; plot(Nmat_tot,'.')
+
+ % pick a pair whose norm of the DIFFERENCE (and its parts) is not an outlier
+ i1 = find( abs( log(Nmat_tot / 1.0 )) < 0.1 );
+ i2 = find( abs( log(Nmat(:,1) / (1/3) )) < 0.1 );
+ i3 = find( abs( log(Nmat(:,2) / (1/3) )) < 0.1 );
+ i4 = find( abs( log(Nmat(:,3) / (1/3) )) < 0.1 );
+ ipairs = mintersect(i1,i2,i3,i4);
+ npair = length(ipairs);
+ %Imat(ipairs,:)
+
+ % plot the pairs
+ for ii=1:npair
+ inds = Imat(ipairs(ii),:);
+ i1 = inds(2);
+ i2 = inds(3);
+
+ g1 = [ts_samples(:,i1) xs_samples(:,i1) ys_samples(:,i1)];
+ g2 = [ts_samples(:,i2) xs_samples(:,i2) ys_samples(:,i2)];
+ g1g2 = g2 - g1;
+
+ % event 5
+ d5 = sqrt( g1g2(5,2)^2 + g1g2(5,3)^2 );
+ t5 = g1g2(5,1);
+ disp(sprintf('%6i%6i%6.2f s%10.1f m',i1,i2,t5,d5));
+
+ if ifig == 1
+ figure; nr=3; nc=3;
+ for kk=1:3
+ subplot(nr,nc,kk); sig = sigs(kk);
+ Nbin = plot_histo(g1(:,kk),[-4*sig : sig/4 : 4*sig]);
+ title(sprintf('sample %i (%.2f)',i1,max(abs(g1(:,kk)))/sig)); xlabel(xlabs2{kk});
+ subplot(nr,nc,3+kk); sig = sigs(kk);
+ Nbin = plot_histo(g2(:,kk),[-4*sig : sig/4 : 4*sig]);
+ title(sprintf('sample %i (%.2f)',i2,max(abs(g2(:,kk)))/sig)); xlabel(xlabs2{kk});
+
+ subplot(nr,nc,6+kk); sig = sigssum(kk);
+ Nbin = plot_histo(g1g2(:,kk),[-4*sig : sig/4 : 4*sig]);
+ title(sprintf('sample %i - %i (%.2f)',i2,i1,max(abs(g1g2(:,kk)))/sig)); xlabel(xlabs2{kk});
+
+ end
+ fontsize(8), orient landscape, wysiwyg;
+ end
+ end
+
+% nr=3; nc=1;
+% for kk = 1:RUNS
+% kk
+%
+% if ifig==1, figure; end
+% Ymat = zeros(N,3);
+% for ii = 1:3 % loop over three parameters
+% sig = sigs(ii);
+% yvec = sig*randn(N,1);
+% Ymat(:,ii) = yvec;
+%
+% if ifig==1
+% subplot(nr,nc,ii);
+% Nbin = plot_histo(yvec,[-5*sig : sig/4 : 5*sig]);
+% grid on; xlabel(['Perturbation in ' xlabs{ii}]);
+% ylim([0 0.2]);
+% end
+% end
+% end
+%
+% for kk = 1:NRUNS
+%
+% % pick the sample based on
+% % (1) the ith event perturbation
+% % (2) the overall norm (close to 1)
+% % (3) the STD of the norms of each parameter (low)
+%
+% % ith event
+% %irec = 5; dloc_min = 4000; dotime_min = 0.8;
+% irec = 5; dloc_min = 4000; dotime_min = 0.2;
+% misloc = sqrt( Ymat(irec,2)^2 + Ymat(irec,3)^2 );
+% otime = Ymat(irec,1);
+%
+% % overall norm
+% norms = sum( Ymat .* Ymat ./ cov_model);
+%
+% if and( and( misloc > dloc_min, abs(otime) > dotime_min ), ...
+% and( abs( sum(norms) - 1 ) <= 0.05, std(norms) <= 0.02) )
+% disp(' Event 5:');
+% disp([' Mislocation (m) : ' num2str(misloc) ]);
+% disp([' Origin time perturbation (s) : ' num2str(otime) ]);
+%
+% figure;
+% for ii = 1:3
+% sig = sigs(ii); yvec =Ymat(:,ii);
+% subplot(nr,nc,ii);
+% Nbin = plot_histo(yvec,[-5*sig : sig/4 : 5*sig]);
+% grid on; xlabel(['Perturbation in ' xlabs{ii}]); ylim([0 0.2]);
+% title(sprintf('%s -- sigma = %.3f -- Norm of %i parameters = %.4f',...
+% xlabs{ii},sig,N,norms(ii)));
+% end
+% orient tall, wysiwyg;
+% break
+% end
+% if ifig==1, orient tall, wysiwyg; end
+% end
+
+ if iwrite==1
+ disp('Are you sure you want to write two samples to files?');
+ i1 = input(' Enter index of initial model: ');
+ i2 = input(' Enter index of target model: ');
+ Y1 = [ts_samples(:,i1) xs_samples(:,i1) ys_samples(:,i1)];
+ Y2 = [ts_samples(:,i2) xs_samples(:,i2) ys_samples(:,i2)];
+ YD = Y2 - Y1;
+
+ % check norms
+ norm1 = [ sum( Y1(:,1).^2 / (N*NPARM_SOURCE*sigs(1)^2) )
+ sum( Y1(:,2).^2 / (N*NPARM_SOURCE*sigs(2)^2) )
+ sum( Y1(:,3).^2 / (N*NPARM_SOURCE*sigs(3)^2) ) ]
+ sum(norm1)
+ norm2 = [ sum( Y2(:,1).^2 / (N*NPARM_SOURCE*sigs(1)^2) )
+ sum( Y2(:,2).^2 / (N*NPARM_SOURCE*sigs(2)^2) )
+ sum( Y2(:,3).^2 / (N*NPARM_SOURCE*sigs(3)^2) ) ]
+ sum(norm1)
+ normD = [ sum( YD(:,1).^2 / (N*NPARM_SOURCE*sigssum(1)^2) )
+ sum( YD(:,2).^2 / (N*NPARM_SOURCE*sigssum(2)^2) )
+ sum( YD(:,3).^2 / (N*NPARM_SOURCE*sigssum(3)^2) ) ]
+ sum(normD)
+
+ ofile = [odir 'events_txy_pert_initial.dat'];
+ fid = fopen(ofile,'w');
+ for ii = 1:N
+ fprintf(fid,'%20.12e%20.12e%20.12e\n', Y1(ii,:));
+ end
+ fclose(fid);
+
+ ofile = [odir 'events_txy_pert_target.dat'];
+ fid = fopen(ofile,'w');
+ for ii = 1:N
+ fprintf(fid,'%20.12e%20.12e%20.12e\n', Y2(ii,:));
+ end
+ fclose(fid);
+
+ ofile = [odir 'events_txy_pert_sigmas.dat'];
+ fid = fopen(ofile,'w');
+ fprintf(fid,'%20.12e%20.12e%20.12e\n',sigs);
+ fclose(fid);
+ end
+
+ %xmin = ans(1); xmax = ans(2);
+ %xpt = linspace(xmin,xmax,100);
+ %ypt = 1/sqrt(2*sig^2) * exp(-xpt.^2 / (2*sig^2) )
+end
+
+%---------------------------------------------------------
+% see what the sine-like checkerboard maps look like as a distribution
+% --> What is the sigma value that characterizes the distribution?
+
+if isigma_struct == 1
+
+ idir = '/net/denali/scratch1/carltape/OUTPUT_2/run_9100/'; % 7100
+ temp = load([idir 'structure_dat.dat']); beta = temp(:,7);
+ N = length(beta);
+
+ sig0 = 0.1;
+ sfac = [1 2 3 4 5];
+
+ for kk = 1:length(sfac);
+ sig = sig0 / sfac(kk);
+
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1);
+ Nbin = plot_histo(beta,[-5*sig : sig/4 : 5*sig]);
+ grid on; xlabel(['Perturbation in ' xlabs{1}]);
+ title('checkerboard with +/- 0.10');
+ ylim([0 0.2]); %xlim([-0.5 0.5]);
+
+ subplot(nr,nc,2);
+ Nbin = plot_histo(sig*randn(N,1),[-5*sig : sig/4 : 5*sig]);
+ grid on; xlabel(['Perturbation in ' xlabs{1}]);
+ title(['sigma = ' sprintf('%.3f',sig)]);
+ ylim([0 0.2]); %xlim([-0.5 0.5]);
+ orient tall, wysiwyg
+ end
+
+ %---------------------
+
+ % load the jacobian for constructing the "event gradient" from the event kernel
+ lmesh_all = load([idir 'local_mesh.dat']);
+ Ai = lmesh_all(:,9); % local dA area element
+ Atot = sum(Ai); % total area
+
+ % load the model covariance
+ cov_model = load([idir 'cov_imetric_diagonal.dat']);
+ cov_str = cov_model(1:NLOCAL);
+
+ sig0 = 0.05
+ beta_const = sig0 * ones(N,1);
+ beta_gaus = sig0 * randn(N,1);
+
+ % compute norms
+ Mnorm0_gaus = sum( beta_gaus.^2 ./ cov_str )
+ Mnorm0 = sum( beta.^2 ./ cov_str )
+ Mnorm0_const = sum( beta_const.^2 ./ cov_str )
+
+ %---------------------
+
+ figure; nr=3; nc=1;
+ bins = [-5*sig0 : sig0/4 : 5*sig0];
+
+ subplot(nr,nc,1);
+ Nbin = plot_histo(sig0*randn(N,1),bins);
+ grid on; xlabel(['Perturbation in ' xlabs{1}]);
+ title(sprintf('Gaussian with sigma = %.3f; norm (m^T C^{-1} m) = %.3f',...
+ sig0,Mnorm0_gaus)); ylim([0 0.2]);
+
+ subplot(nr,nc,2);
+ Nbin = plot_histo(beta,bins);
+ grid on; xlabel(['Perturbation in ' xlabs{1}]);
+ title(sprintf('Checkerboard with amplitude 0.100; norm (m^T C^{-1} m) = %.3f',...
+ sig0,Mnorm0)); ylim([0 0.2]);
+
+ subplot(nr,nc,3);
+ Nbin = plot_histo(beta_const,bins);
+ grid on; xlabel(['Perturbation in ' xlabs{1}]);
+ title(sprintf('Constant function with values = %.3f; norm (m^T C^{-1} m) = %.3f',...
+ sig0,Mnorm0_const)); ylim([0 0.2]);
+
+ orient tall, wysiwyg
+
+ %---------------------------
+
+ sigvec = 10.^linspace(-2,0,121); % test sigma values
+ for kk = 1:length(sigvec);
+ sig = sigvec(kk);
+ Mnorm(kk) = (1/Atot) * sum(beta.^2 / sig^2 .* Ai);
+ Mnorm_const(kk) = (1/Atot) * sum(beta_const.^2 / sig^2 .* Ai);
+ end
+
+ figure; nr=2; nc=2;
+ st1 = ['checkerboard (PERT = 0.100)'];
+ st2 = ['uniform (PERT = ' sprintf('%.3f',sig0) ')'];
+
+ subplot(nr,nc,1);
+ plot(sigvec, Mnorm,'b.', sigvec, Mnorm_const,'r.'); grid on;
+ legend(st1,st2); xlabel(' sigma value'); ylabel(' Norm');
+
+ subplot(nr,nc,2);
+ loglog(sigvec, Mnorm,'b.', sigvec, Mnorm_const,'r.'); grid on;
+ legend(st1,st2); xlabel(' sigma value'); ylabel(' Norm');
+
+ subplot(nr,nc,3);
+ plot(sigvec, Mnorm,'b.', sigvec, Mnorm_const,'r.'); grid on;
+ legend(st1,st2); xlabel(' sigma value'); ylabel(' Norm');
+ axis([0.03 0.12 0 2]);
+
+ orient tall, wysiwyg
+end
+
+%---------------------------------------------------------
+% see what the sine-like checkerboard maps look like as a distribution
+% --> What is the sigma value that characterizes the distribution?
+
+if isigma_data == 1
+
+ %-------------------------------
+ % load traveltime measurments from 2D code
+
+ idir = '/net/denali/scratch1/carltape/OUTPUT_2/run_7000/';
+ meas_all = load([idir 'measure_vec.dat']);
+ dT_all = meas_all(:,2);
+ Nmeas = length(dT_all);
+
+ sig = 1.7;
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1); Nbin = plot_histo(dT_all,[-5*sig : sig/4 : 5*sig]); title('Actual data'); grid on;
+ subplot(nr,nc,2); Nbin = plot_histo(sig*randn(Nmeas,1),[-5*sig : sig/4 : 5*sig]); grid on;
+ title(['Gaussian distribution with sigma = ' sprintf('%.3f',sig)]);
+ orient tall, wysiwyg
+
+ %-------------------------------
+ % create errors to add to data
+
+ sig = 0.1; stag = '0p1';
+ N = 10000;
+
+ vals = sig*randn(N,1);
+
+ figure; Nbin = plot_histo(vals,[-5*sig : sig/4 : 5*sig]);
+ grid on; xlabel('Perturbation in measurement');
+ title(['sigma = ' sprintf('%.3f',sig) ', std = ' sprintf('%.6f',std(vals)) ]);
+
+ if iwrite==1
+ ofile = [odir 'sigma_' stag '_pert.dat'];
+ fid = fopen(ofile,'w');
+ for ii = 1:length(vals)
+ fprintf(fid,'%20.12e\n', vals(ii));
+ end
+ fclose(fid);
+
+ % check
+ sigt = load(ofile);
+ figure; Nbin = plot_histo(sigt,[-5*sig : sig/4 : 5*sig]);
+ grid on; xlabel('Perturbation in measurement');
+ end
+end
+
+%=============================================================
+
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-01-31 21:45:03 UTC (rev 16202)
@@ -130,6 +130,8 @@
! socal earthquake location
double precision, dimension(MAX_EVENT) :: socal_events_lon, socal_events_lat, socal_events_mag
+ double precision, dimension(MAX_EVENT) :: x_eve_syn_pert, z_eve_syn_pert, otime_eve_syn_pert, otime
+ double precision, dimension(MAX_EVENT) :: x_eve_dat_pert, z_eve_dat_pert, otime_eve_dat_pert
integer, dimension(MAX_EVENT) :: socal_events_lab
integer nmeas, nmeas_run, ievent_min, ievent_max, ifirst_event
@@ -151,14 +153,14 @@
double precision, dimension(:), allocatable :: cov_imetric, icov_metric, cov_model, icov_model
! conjugate gradient optimization (using line-search)
- double precision, dimension(:), allocatable :: pk, gk, g0, p0, m0, gt, mt, mdiff ! m0_prior
+ 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 :: 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, m0_vec_initial, mtarget
- 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
+ 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
double precision :: vel_min, vel_max
integer :: itest, nmod, nmod_src, nmod_str
@@ -236,6 +238,16 @@
if(ISRC_SPACE==1) FINITE_SOURCE = 0
if(ISRC_SPACE/=1) FINITE_SOURCE = 1
+ origin_time = tshift ! reference origin time (s)
+ origin_time_dat = origin_time
+ origin_time_syn = origin_time
+ otime(:) = origin_time ! origin time for prior model
+
+ ! STANDARD DEVIATIONS of source perturbations (see wave2d_sigmas.m)
+ ! NOTE: THESE ARE OVER-WRITTEN FOR ISURFACE=1
+ src_pert_dist = 2000.0 ! source spatial perturbation (m)
+ src_pert_time = 0.5 ! source temporal perturbation (s)
+
!--------------------------------------
! load or specify events, depending on whether you are computing
! membrane surface waves (ISURFACE==1) or body waves
@@ -263,6 +275,80 @@
enddo
!----------
+ ! read in perturbations for initial and target models
+
+ otime_eve_syn_pert(:) = 0.0 ; x_eve_syn_pert(:) = 0.0 ; z_eve_syn_pert(:) = 0.0
+ otime_eve_dat_pert(:) = 0.0 ; x_eve_dat_pert(:) = 0.0 ; z_eve_dat_pert(:) = 0.0
+
+ if(0==1) then ! NEW perturbations
+
+ ! perturbations generated from wave2d_sigmas.m
+ open(19,file='INPUT/events_txy_pert_initial.dat',status='old')
+ do ievent = 1,nevent
+ read(19,*) temp1,temp2,temp3
+ otime_eve_syn_pert(ievent) = temp1
+ x_eve_syn_pert(ievent) = temp2
+ z_eve_syn_pert(ievent) = temp3
+ enddo
+ close(19)
+
+ ! perturbations generated from wave2d_sigmas.m
+ open(20,file='INPUT/events_txy_pert_target.dat',status='old')
+ do ievent = 1,nevent
+ read(20,*) temp1,temp2,temp3
+ otime_eve_dat_pert(ievent) = temp1
+ x_eve_dat_pert(ievent) = temp2
+ z_eve_dat_pert(ievent) = temp3
+ enddo
+ close(20)
+
+ ! read sigma values used to generate the source parameter perturbations
+ ! NOTE: these over-write what is defined above
+ open(21,file='INPUT/events_txy_pert_sigmas.dat',status='old')
+ read(21,*) src_pert_time,src_pert_dist,temp1
+ close(21)
+
+ else ! OLD perturbations
+
+ ! perturbations generated from gji_src_inversion.m
+ ! NOTE: OLD ordering scheme for source indexing
+ open(19,file='INPUT/events_xyt_pert.dat',status='unknown')
+ do ievent = 1,nevent
+ read(19,*) temp1,temp2,temp3
+ !if(ievent==5) then
+ ! temp1 = -1.599978278510d3
+ ! temp2 = -6.502537573573d2
+ ! temp3 = 7.975610515441d-2
+ !endif
+ x_eve_dat_pert(ievent) = temp1
+ z_eve_dat_pert(ievent) = temp2
+ otime_eve_dat_pert(ievent) = temp3
+ enddo
+ close(19)
+
+ open(20,file='INPUT/events_xyt_pert_sigmas.dat',status='unknown')
+ read(20,*) temp1,src_pert_dist,src_pert_time
+ close(20)
+
+ endif
+
+ ! no perturbations for synthetics --> m0 = mprior
+ if (M0ISMPRIOR) then
+ otime_eve_syn_pert(:) = 0.0
+ x_eve_syn_pert(:) = 0.0
+ z_eve_syn_pert(:) = 0.0
+ endif
+
+!!$ ! check
+!!$ do ievent = 1,nevent
+!!$ write(*,'(a20,i6,3f12.3)') 'Initial model :',ievent, &
+!!$ otime_eve_syn_pert(ievent), x_eve_syn_pert(ievent), z_eve_syn_pert(ievent)
+!!$ write(*,'(a20,i6,3f12.3)') 'Target model :',ievent, &
+!!$ otime_eve_dat_pert(ievent), x_eve_dat_pert(ievent), z_eve_dat_pert(ievent)
+!!$ enddo
+!!$ stop 'testing perturbations'
+
+ !----------
! load socal coast and shelf points
! (This is useful for some ocean microseismi time-rerversal experiments.)
@@ -293,6 +379,7 @@
!write(*,'(4f16.6)') (coast_lon(i), coast_lat(i), coast_x(i)/1000.0, coast_z(i)/1000.0, i = 1,ncoast)
else
+
!----------
! events for data
@@ -373,15 +460,6 @@
!src_pert_dist = 5000.0 ! source spatial perturbation (m)
!src_pert_time = 1.0 ! source temporal perturbation (s)
- ! STANDARD DEVIATIONS of source perturbations (see wave2d_sigmas.m)
- ! NOTE: These should be READ IN below
- src_pert_dist = 2000.0 ! source spatial perturbation (m)
- src_pert_time = 0.5 ! source temporal perturbation (s)
-
- origin_time = tshift ! reference origin time (s)
- origin_time_dat = origin_time
- origin_time_syn = origin_time
-
! reference directory for the run
!irunz = 500
@@ -1012,13 +1090,22 @@
allocate(otime_syn(nevent))
! origin time for synthetics
- otime_syn(:) = origin_time
+ ! initialize to mprior
+ otime_syn(1:nevent) = otime(1:nevent) ! otime_syn is length MAX_EVENT, otime is length nevent
+ if (PERT_SOURCE_T == 1) then
+ otime_syn(:) = otime_syn(:) + otime_eve_syn_pert(:)
+ endif
! assign initial synthetic events to (x_eve_syn, z_eve_syn)
+ ! initialize to mprior
do i = 1,nevent
x_eve_syn(i) = x_eve0_syn(ifilter_eve_syn(i))
z_eve_syn(i) = z_eve0_syn(ifilter_eve_syn(i))
enddo
+ if (PERT_SOURCE_X == 1) then
+ x_eve_syn(:) = x_eve_syn(:) + x_eve_syn_pert(:)
+ z_eve_syn(:) = z_eve_syn(:) + z_eve_syn_pert(:)
+ endif
!write(*,'(2e24.16)') (x_eve_syn(i), z_eve_syn(i), i = 1,nevent)
@@ -1037,7 +1124,7 @@
! display target events and final events -- and the distance between (in meters)
if(ISURFACE==1) then
- ! convert from mesh to lat-lon
+ ! convert from mesh to lat-lon (gets x_eve_lon_syn and z_eve_lat_syn)
call mesh_geo(nevent, x_eve_lon_syn, z_eve_lat_syn, x_eve_syn, z_eve_syn, UTM_PROJECTION_ZONE, IMESH2LONLAT)
! The distance error is due to the UTM conversion, but the lat-lon points are only
@@ -1066,6 +1153,22 @@
enddo
endif
+ ! write events for PRIOR MODEL to file
+ if( .not. READ_IN ) then
+ open(19,file=trim(out_dir2)//'events_prior_xz.dat',status='unknown')
+ do ievent = 1,nevent
+ write(19,'(3f20.10,1i12)') x_eve0_syn(ievent), z_eve0_syn(ievent), otime(ievent), ievent
+ enddo
+ close(19)
+ if(ISURFACE==1) then
+ open(18,file=trim(out_dir2)//'events_prior_lonlat.dat',status='unknown')
+ do ievent = 1,nevent
+ write(18,'(2f20.10,1i12)') x_eve_lon0_syn(ievent), z_eve_lat0_syn(ievent), ievent
+ enddo
+ close(18)
+ endif
+ endif
+
! write events for SYNTHETICS to file
if( .not. READ_IN ) then
open(19,file=trim(out_dir2)//'events_syn_xz.dat',status='unknown')
@@ -1100,28 +1203,41 @@
if( PERT_SOURCE_T == 1 .or. PERT_SOURCE_X == 1 ) then ! source perturbations
if(ISURFACE==1) then ! read in PERTURBED events from another file (for GJI-2007 simulations)
+ ! initialize to no perturbations from synthetics
+ !otime_dat(:) = otime_syn(:)
+ !x_eve_dat(:) = x_eve_syn(:)
+ !z_eve_dat(:) = z_eve_syn(:)
- ! perturbations generated from gji_src_inversion.m
- ! NOTE: OLD ordering scheme for source indexing
- open(19,file='INPUT/events_xyt_pert.dat',status='unknown')
- do ievent = 1,25
- read(19,*) temp1,temp2,temp3
- if( PERT_SOURCE_T == 0 ) temp3 = 0.0
- if( PERT_SOURCE_X == 0 ) then
- temp1 = 0.0 ; temp2 = 0.0
- endif
- x_eve_dat(ievent) = x_eve_syn(ievent) + temp1
- z_eve_dat(ievent) = z_eve_syn(ievent) + temp2
- otime_dat(ievent) = otime_syn(ievent) + temp3
- enddo
- close(19)
+ ! initialize to prior model
+ ! NOTE: these arrays have different lengths
+ otime_dat(1:nevent) = otime(1:nevent)
+ x_eve_dat(1:nevent) = x_eve0_syn(1:nevent)
+ z_eve_dat(1:nevent) = z_eve0_syn(1:nevent)
- ! read sigma values used to generate the source parameter perturbations
- ! NOTE: these over-write what is defined above
- open(20,file='INPUT/events_xyt_pert_sigmas.dat',status='unknown')
- read(20,*) temp1,src_pert_dist,src_pert_time
- close(20)
+ if (PERT_SOURCE_T == 1) then
+ otime_dat(:) = otime_dat(:) + otime_eve_dat_pert(:)
+ endif
+
+ if (PERT_SOURCE_X == 1) then
+ x_eve_dat(:) = x_eve_dat(:) + x_eve_dat_pert(:)
+ z_eve_dat(:) = z_eve_dat(:) + z_eve_dat_pert(:)
+ endif
+!!$ ! perturbations generated from gji_src_inversion.m
+!!$ ! NOTE: OLD ordering scheme for source indexing
+!!$ open(19,file='INPUT/events_xyt_pert.dat',status='unknown')
+!!$ do ievent = 1,25
+!!$ read(19,*) temp1,temp2,temp3
+!!$ if( PERT_SOURCE_T == 0 ) temp3 = 0.0
+!!$ if( PERT_SOURCE_X == 0 ) then
+!!$ temp1 = 0.0 ; temp2 = 0.0
+!!$ endif
+!!$ x_eve_dat(ievent) = x_eve_syn(ievent) + temp1
+!!$ z_eve_dat(ievent) = z_eve_syn(ievent) + temp2
+!!$ otime_dat(ievent) = otime_syn(ievent) + temp3
+!!$ enddo
+!!$ close(19)
+
!!$ open(19,file='OUTPUT/run_2550/events_xy.dat',status='unknown')
!!$ do ievent = 1,25
!!$ read(19,'(3f20.10,1i12)') x_eve_dat(ievent), z_eve_dat(ievent), otime_dat(ievent), itemp1
@@ -1604,11 +1720,12 @@
close(18)
! allocate model vector
- allocate(m0(nmod),mt(nmod),mdiff(nmod))
+ allocate(m0(nmod),mt(nmod),mdiff(nmod),mprior(nmod))
allocate(cov_imetric(nmod),icov_metric(nmod),cov_model(nmod),icov_model(nmod))
allocate(gradient(nmod),gradient_data(nmod),gradient_model(nmod))
- allocate(m0_vec(nmod),mt_vec(nmod),m0_vec_initial(nmod))
+ allocate(m0_vec(nmod),mt_vec(nmod))
allocate(g0(nmod),gt(nmod),gk(nmod),p0(nmod),pk(nmod))
+ !allocate(m0_vec_initial(nmod))
! gradient values for each parameter
allocate(gradient_data_all(nevent,NVAR))
@@ -1620,23 +1737,26 @@
mt(:) = 0.0
m0_vec(:) = 0.0
mt_vec(:) = 0.0
- !m0_prior(:) = 0.0
+ mprior(:) = 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),m_src_syn_vec_initial(nmod_src),m_src_syn_initial(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_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))
+ 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_syn_vec_initial(:) = 0.0
+ !m_src_syn_vec_initial(:) = 0.0
! fill source model vector for SYNTHETICS AND DATA
! NOTE: Only fill the source parameters that are being RUN;
@@ -1661,21 +1781,46 @@
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_prior(itemp1) = otime(ievent)
+ m_src_prior(itemp2) = x_eve0_syn(ievent)
+ m_src_prior(itemp3) = z_eve0_syn(ievent)
enddo
!------------------
- ! initial model vector (physical coordinates)
+ ! PRIOR and INITIAL model vectors
+ ! initial source parameters for SYNTHETICS
+ !m_src_syn_vec_initial(:) = m_src_syn_vec(:)
+ !m_src_syn_initial(:) = m_src_syn_vec(:) - m_src_syn_vec_initial(:) ! = 0.
+
+ ! initial source parameters for DATA and SYNTHETICS
+ ! 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)
+ !m0_vec_initial(:) = m0_vec(:)
- !m0_prior(:) = m0(:)
- 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)
+
+ ! write out initial model vector
if( .not. READ_IN ) then
open(88,file=trim(out_dir2)//'initial_model_vector.dat',status='unknown')
do i = 1,nmod
@@ -1684,35 +1829,21 @@
close(88)
endif
- !stop 'testing'
-
- ! apply uniform scaling
- !m_scale_src(:) = m_scale_src(:) / joint_str_src_grad_ratio
-
-!!$ do ievent = 1,nevent
-!!$ do i = 1,3
-!!$ itemp = (ievent-1)*3 + i
-!!$ m_scale_src_all(itemp) = m_scale_src(i)
-!!$ enddo
-!!$ enddo
-
- ! initial source parameters for SYNTHETICS
- m_src_syn_vec_initial(:) = m_src_syn_vec(:)
- m_src_syn_initial(:) = m_src_syn_vec(:) - m_src_syn_vec_initial(:) ! = 0.
-
- ! initial source parameters for DATA and SYNTHETICS
- ! 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(:)
-
! 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_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)
enddo
close(88)
+
+ 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
!------------------
@@ -1807,7 +1938,7 @@
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_total = (ugrad_ts/fac_ts) + (ugrad_xs/fac_xs) + (ugrad_ys/fac_ys)
elseif ( 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)
@@ -1816,21 +1947,20 @@
!ugrad_str = 0.2070999127d6 ; ugrad_ts = 0.3507740673d4 ; ugrad_xs = 0.1940368586d4 ; ugrad_ys = 0.2142118702d4 ! Nfac = 3, 25 events
!ugrad_str = 0.4247330856d6 ; ugrad_ts = 0.5896527525d4 ; ugrad_xs = 0.2643538184d4 ; ugrad_ys = 0.3342815391d4 ! Nfac = 3, 5 events
- ugrad_str = 0.1205146534d6 ; ugrad_ts = 0.2926409454d3 ; ugrad_xs = 0.9695606936d3 ; ugrad_ys = 0.7224889563d3 ! Nfac = 3, 1 event
+ !ugrad_str = 0.1205146534d6 ; ugrad_ts = 0.2926409454d3 ; ugrad_xs = 0.9695606936d3 ; ugrad_ys = 0.7224889563d3 ! Nfac = 3, 1 event
- !ugrad_str = ; ugrad_ts = ; ugrad_xs = ; ugrad_ys = ! Gaussians
- !ugrad_str = 0.8251756295d5 ; ugrad_ts = 0.4212723180d3 ; ugrad_xs = 0.5753336089d2 ; ugrad_ys = 0.4077345971d2 ! Gaussians, 1 event
- !ugrad_str = 0.3812017066d6 ; ugrad_ts = 0.2140202263d4 ; ugrad_xs = 0.2970951189d4 ; ugrad_ys = 0.4102594535d4 ! Gaussians, 5 events
- !ugrad_str = 0.1907942047d6 ; ugrad_ts = 0.1687351289d4 ; ugrad_xs = 0.2308837150d4 ; ugrad_ys = 0.2986802305d4 ! Gaussians, 25 events
-
+ ! 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.3812017066d6 ; ugrad_ts = 0.2140202263d4 ; ugrad_xs = 0.2970951189d4 ; ugrad_ys = 0.4102594535d4 ! Gaussians, 5 events
- !ugrad_str = 0.1907942047d6 ; ugrad_ts = 0.1687351289d4 ; ugrad_xs = 0.2308837150d4 ; ugrad_ys = 0.2986802305d4 ! Gaussians, 25 events
+ ! 3200/5000
+ 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/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.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)
endif
@@ -1857,10 +1987,7 @@
coverage_str = 1.0
coverage_src = 1.0
- ! structure part
cov_model(m_inds(1,1):m_inds(1,2)) = ( sigma_beta )**2 / da_local_vec(:) * AREA * coverage_str
-
- ! source part
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
@@ -1884,7 +2011,8 @@
call local2mvec(temp_local1, nmod_src, m_src_dat, nmod, mtarget)
! possible stopping criteria based on the target model
- ! NOTE: THIS IS NOT USED
+ ! NOTE 1: this depends on what you pick as your initial model (prior mean, or a sample)
+ ! NOTE 2: THIS IS NOT USED
!chi_model_stop = 0.5 * model_target_norm
chi_model_stop = 0.0
@@ -2106,7 +2234,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
+ m_src_syn_vec(:) = m_src_dat_vec(:) ! 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
@@ -2123,53 +2251,71 @@
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(:) - m_src_syn_vec_initial(:)
+ m0(nmod_str+1 : nmod) = m_src_syn_vec(:)
endif
+ !--------------------------
! compute model norm term in the misfit function
- ! NOTE 1: The prior model is always all ZEROS (B,Ts,Xs,Ys);
- ! hence, we replace (m0 - mprior) --> m0 and (mt - mprior) --> mt .
- ! NOTE 2: The initial model is NOT necessarily the prior model .
- ! NOTE 3: THE FACTOR OF 0.5 IS NOT INCLUDED HERE
+ ! NOTE 1: TWO VERSIONS using the balanced Cm (cov_imetric) and the original Cm (cov_model)
+ ! 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)
+
+ ! model vector, balanced Cm (cov_imetric)
+ 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)
-
+ 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_imetric, &
- model_norm, model_norm_struct, model_norm_source, model_norm_parts)
+ model_norm, model_norm_struct, model_norm_source, model_norm_parts, mprior)
endif
+ filename = trim(out_dir1)//'model_norm_all.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)
+
+ ! 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)
+
! compute the difference between the present model and the target model
- ! NOTE: We used the UN-BALANCED model covariance for the norms
mdiff(:) = 0.0
if(itest==0) then
- mdiff = m0 - mtarget
+ mdiff = mtarget - m0
elseif(itest==1) then
- mdiff = mt - mtarget
+ mdiff = mtarget - mt
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)
- ! write model norms to file
- filename1 = trim(out_dir1)//'model_norm_all.dat'
- filename2 = trim(out_dir1)//'model_diff_norm_all.dat'
- call write_norm_sq(filename1,model_norm_parts)
- call write_norm_sq(filename2,model_diff_norm_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)
- ! write unbalanced target model norm to file (cov_model)
- filename = trim(out_dir1)//'model_target_norm_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)
- call write_norm_sq(filename,model_target_norm_parts)
-
- ! write balanced target model norm to file (cov_imetric)
- filename = trim(out_dir1)//'model_target_norm_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)
- call write_norm_sq(filename,model_target_norm_parts)
-
!stop 'TESTING NORMS'
!--------------------------
@@ -2556,19 +2702,19 @@
write(91,*) 'istep, imod, ievent, irun : '
write(91,*) istep, imod, ievent, irun
write(91,*) 'SOURCE MODEL (xs, ys, t0) :'
- write(91,*) ' [m_src_syn_vec]:'
+ write(91,'(a12,3a18)') ' ','m_src_syn_vec','m_src_dat_vec','dat - syn'
write(91,'(a12,3f18.8)') ' t0, s : ', &
m_src_syn_vec(itemp1), &
m_src_dat_vec(itemp1), &
- m_src_syn_vec(itemp1) - m_src_dat_vec(itemp1)
+ m_src_dat_vec(itemp1) - m_src_syn_vec(itemp1)
write(91,'(a12,3f18.8)') ' xs, km : ', &
m_src_syn_vec(itemp2)/1000.0, &
m_src_dat_vec(itemp2)/1000.0, &
- (m_src_syn_vec(itemp2) - m_src_dat_vec(itemp2))/1000.0
+ (m_src_dat_vec(itemp2) - m_src_syn_vec(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_syn_vec(itemp3) - m_src_dat_vec(itemp3))/1000.0
+ (m_src_dat_vec(itemp3) - m_src_syn_vec(itemp3))/1000.0
endif
!--------------------------------------
@@ -2745,9 +2891,9 @@
! 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_syn_vec(itemp1) - m_src_dat_vec(itemp1), &
- m_src_syn_vec(itemp2) - m_src_dat_vec(itemp2), &
- m_src_syn_vec(itemp3) - m_src_dat_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)
enddo
close(19)
close(20)
@@ -3138,7 +3284,7 @@
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)') ' chi_model_stop(mtarget) : ', chi_model_stop
+ !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
endif
@@ -3180,12 +3326,17 @@
if (itest==0) then
- ! m_src_syn: w.r.t initial source values (m_src_syn_vec_initial)
m0(:) = 0.0
temp_local1(:,:,:) = log( beta_syn(:,:,:) / beta0 )
- m_src_syn(:) = m_src_syn_vec(:) - m0_vec_initial(nmod_str+1 : nmod)
+ 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)
+!!$ m0(:) = 0.0
+!!$ temp_local1(:,:,:) = log( beta_syn(:,:,:) / beta0 )
+!!$ m_src_syn(:) = m_src_syn_vec(:) - m0_vec_initial(nmod_str+1 : nmod)
+!!$ call local2mvec(temp_local1, nmod_src, m_src_syn, nmod, m0)
+
endif
!-----------------------------------------------------
@@ -3195,15 +3346,15 @@
if( itest==0 .or. POLY_ORDER==3 ) then
print *, ' Computing the gradient of the misfit function for a given model'
- gradient(:) = 0.0 ; gradient_data(:) = 0.0 ; gradient_model(:) = 0.0
+ !gradient(:) = 0.0 ; gradient_data(:) = 0.0 ; gradient_model(:) = 0.0 ! initialized above
+ ! source gradient
+ open(unit=19,file=trim(out_dir1)//'source_gradient.dat',status='unknown')
+ write(19,'(1e20.10)') (source_gradient(i), i = 1,nmod_src)
+ close(19)
+
if(INV_STRUCT_BETA == 1) then
- ! source gradient
- open(unit=19,file=trim(out_dir1)//'source_gradient.dat',status='unknown')
- write(19,'(1e20.10)') (source_gradient(i), i = 1,nmod_src)
- close(19)
-
! norm of parts of unbalanced gradient for all events
! NOTE: within each event directory there is also gradient_data_unbalanced.dat
open(unit=18,file=trim(out_dir1)//'gradient_data_all_unbalanced.dat',status='unknown')
@@ -3243,13 +3394,13 @@
! add the gradient term associated with the MODEL NORM in the misfit function
if( INCLUDE_MODEL_NORM ) then
if (itest == 0) then
- gradient_model(:) = m0(:) / cov_imetric(:)
+ gradient_model(:) = (m0(:) - mprior(:)) / cov_imetric(:)
else
- gradient_model(:) = mt(:) / cov_imetric(:)
+ gradient_model(:) = (mt(:) - mprior(:)) / cov_imetric(:)
endif
endif
- ! if NOT inverting for structure or source, then set those parts of the gradient to zero
+ ! KEY: if NOT inverting for structure or source, then set those parts of the gradient to zero
if(INV_STRUCT_BETA == 0) then
gradient_data(m_inds(1,1) : m_inds(1,2)) = 0.0
gradient_model(m_inds(1,1) : m_inds(1,2)) = 0.0
@@ -3259,7 +3410,7 @@
gradient_model(m_inds(2,1) : m_inds(2,2)) = 0.0
endif
if(INV_SOURCE_X == 0) then
- gradient_data( m_inds(3,1) : m_inds(4,2)) = 0.0
+ gradient_data(m_inds(3,1) : m_inds(4,2)) = 0.0
gradient_model(m_inds(3,1) : m_inds(4,2)) = 0.0
endif
@@ -3440,21 +3591,24 @@
!!$ endif
!!$ enddo
- ! STRUCTURE PARAMETERS (bulk, then beta)
- if(INV_STRUCT_BETA == 0) then
- mt_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
- else
- mt_vec(1 : NLOCAL) = beta0 * exp( mt(1 : NLOCAL) )
- !mt_vec(1 : NLOCAL) = bulk0 * exp( mt(1 : NLOCAL) )
- !mt_vec(NLOCAL+1 : nmod_str) = beta0 * exp( mt(NLOCAL+1 : nmod_str) )
- endif
+ mt_vec(1 : NLOCAL) = beta0 * exp( mt(1 : NLOCAL) )
+ mt_vec(nmod_str+1 : nmod) = mt(nmod_str+1 : nmod)
- ! SOURCE PARAMETERS
- if(INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0) then
- mt_vec(nmod_str+1 : nmod) = m0_vec_initial(nmod_str+1 : nmod) ! initial source
- else
- mt_vec(nmod_str+1 : nmod) = mt(nmod_str+1 : nmod) + m0_vec_initial(nmod_str+1 : nmod)
- endif
+!!$ ! STRUCTURE PARAMETERS (bulk, then beta)
+!!$ if(INV_STRUCT_BETA == 0) then
+!!$ mt_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
+!!$ else
+!!$ mt_vec(1 : NLOCAL) = beta0 * exp( mt(1 : NLOCAL) )
+!!$ !mt_vec(1 : NLOCAL) = bulk0 * exp( mt(1 : NLOCAL) )
+!!$ !mt_vec(NLOCAL+1 : nmod_str) = beta0 * exp( mt(NLOCAL+1 : nmod_str) )
+!!$ endif
+!!$
+!!$ ! SOURCE PARAMETERS
+!!$ if(INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0) then
+!!$ mt_vec(nmod_str+1 : nmod) = m0_vec_initial(nmod_str+1 : nmod) ! initial source
+!!$ else
+!!$ mt_vec(nmod_str+1 : nmod) = mt(nmod_str+1 : nmod) + m0_vec_initial(nmod_str+1 : nmod)
+!!$ endif
print *, 'lam_t_val = ', sngl(lam_t_val)
@@ -3556,21 +3710,24 @@
!!$ endif
!!$ enddo
- ! STRUCTURE PARAMETERS
- if(INV_STRUCT_BETA == 0) then
- m0_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
- else
- m0_vec(1 : NLOCAL) = beta0 * exp( m0(1 : NLOCAL) )
- !m0_vec(1 : NLOCAL) = bulk0 * exp( m0(1 : NLOCAL) )
- !m0_vec(NLOCAL+1 : nmod_str) = beta0 * exp( m0(NLOCAL+1 : nmod_str) )
- endif
+ m0_vec(1 : NLOCAL) = beta0 * exp( m0(1 : NLOCAL) )
+ m0_vec(nmod_str+1 : nmod) = m0(nmod_str+1 : nmod)
- ! SOURCE PARAMETERS
- if(INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0) then
- m0_vec(nmod_str+1 : nmod) = m0_vec_initial(nmod_str+1 : nmod) ! initial source
- else
- m0_vec(nmod_str+1 : nmod) = m0(nmod_str+1 : nmod) + m0_vec_initial(nmod_str+1 : nmod)
- endif
+!!$ ! STRUCTURE PARAMETERS
+!!$ if(INV_STRUCT_BETA == 0) then
+!!$ m0_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
+!!$ else
+!!$ m0_vec(1 : NLOCAL) = beta0 * exp( m0(1 : NLOCAL) )
+!!$ !m0_vec(1 : NLOCAL) = bulk0 * exp( m0(1 : NLOCAL) )
+!!$ !m0_vec(NLOCAL+1 : nmod_str) = beta0 * exp( m0(NLOCAL+1 : nmod_str) )
+!!$ endif
+!!$
+!!$ ! SOURCE PARAMETERS
+!!$ if(INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0) then
+!!$ m0_vec(nmod_str+1 : nmod) = m0_vec_initial(nmod_str+1 : nmod) ! initial source
+!!$ else
+!!$ m0_vec(nmod_str+1 : nmod) = m0(nmod_str+1 : nmod) + m0_vec_initial(nmod_str+1 : nmod)
+!!$ endif
print *, 'lam_0_val = ', sngl(lam_0_val)
@@ -3637,15 +3794,16 @@
deallocate(otime_syn,otime_dat)
! conjugate-gradient variables
- deallocate(m0,mt,mdiff,mtarget) ! m0_prior
+ deallocate(m0,mt,mdiff,mtarget,mprior)
deallocate(cov_imetric,icov_metric,cov_model,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,m0_vec_initial)
- deallocate(m_src_syn_vec,m_src_syn,m_src_syn_vec_initial,m_src_syn_initial)
+ deallocate(m0_vec,mt_vec)
+ deallocate(m_src_syn_vec,m_src_syn)
deallocate(m_src_dat_vec,m_src_dat)
+ !deallocate(m0_vec_initial,m_src_syn_vec_initial,m_src_syn_initial)
!deallocate(m_scale_src_all)
deallocate(cov_data,index_data,index_source)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90 2010-01-31 21:45:03 UTC (rev 16202)
@@ -99,6 +99,7 @@
! reference model and target model choice
integer, parameter :: IMODEL_SYN = 0
integer, parameter :: IMODEL_DAT = 2
+ logical, parameter :: M0ISMPRIOR = .true.
!-----------------------------------------------------------------------------------------
! 0 1 2 3
!-----------------------------------------------------------------------------------------
@@ -180,7 +181,7 @@
! VAR_RED_MIN : minimum variance reduction
! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
- integer, parameter :: NITERATION = 1
+ integer, parameter :: NITERATION = 2
double precision, parameter :: VAR_RED_MIN = 0.05d0
!double precision, parameter :: SIGMA_FAC = 2.0d0
!double precision, parameter :: CONV_STOP = 1.0d-4
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90 2010-01-31 21:45:03 UTC (rev 16202)
@@ -99,6 +99,7 @@
! reference model and target model choice
integer, parameter :: IMODEL_SYN = 0
integer, parameter :: IMODEL_DAT = 2
+ logical, parameter :: M0ISMPRIOR = .true.
!-----------------------------------------------------------------------------------------
! 0 1 2 3
!-----------------------------------------------------------------------------------------
@@ -180,7 +181,7 @@
! VAR_RED_MIN : minimum variance reduction
! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
- integer, parameter :: NITERATION = 1
+ integer, parameter :: NITERATION = 2
double precision, parameter :: VAR_RED_MIN = 0.05d0
!double precision, parameter :: SIGMA_FAC = 2.0d0
!double precision, parameter :: CONV_STOP = 1.0d-4
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90 2010-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90 2010-01-31 21:45:03 UTC (rev 16202)
@@ -99,6 +99,7 @@
! reference model and target model choice
integer, parameter :: IMODEL_SYN = 0
integer, parameter :: IMODEL_DAT = 2
+ logical, parameter :: M0ISMPRIOR = .true.
!-----------------------------------------------------------------------------------------
! 0 1 2 3
!-----------------------------------------------------------------------------------------
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-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 2010-01-31 21:45:03 UTC (rev 16202)
@@ -5,7 +5,7 @@
! Several lines need to be commented/uncommented to go from one to the other.
! index of the reference directory for the simulation output
- integer, parameter :: IRUNZ = 800
+ integer, parameter :: IRUNZ = 3200
!========================
! GRID, TIME-STEP, AND SOURCE PARAMETERS
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-01-31 03:16:47 UTC (rev 16201)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90 2010-01-31 21:45:03 UTC (rev 16202)
@@ -72,6 +72,8 @@
write(12,*) 'PRHO',PRHO
write(12,*) 'IMODEL_SYN',IMODEL_SYN
write(12,*) 'IMODEL_DAT',IMODEL_DAT
+ if(M0ISMPRIOR) write(12,*) 'M0ISMPRIOR 1'
+ if(.not. M0ISMPRIOR) write(12,*) 'M0ISMPRIOR 0'
write(12,*) 'ISMOOTH_EVENT_KERNEL',ISMOOTH_EVENT_KERNEL
write(12,*) 'ISMOOTH_MISFIT_KERNEL',ISMOOTH_MISFIT_KERNEL
write(12,*) 'ISMOOTH_INITIAL_MODEL',ISMOOTH_INITIAL_MODEL
@@ -1120,28 +1122,36 @@
subroutine compute_norm_sq(ievent_min, ievent_max, nevent, &
index_source, nmod, mvec, cov_model, &
- norm_total, norm_struct, norm_source, norm_parts)
+ norm_total, norm_struct, norm_source, norm_parts, &
+ mvec_prior)
! 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)
integer, intent(in) :: 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, 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
integer :: ievent, itemp1, itemp2, itemp3
!----------
norm_parts(:) = 0.0
+ mtemp = mvec
+ if (present(mvec_prior)) mtemp = mvec - mvec_prior
+
! structure part of the norm -- BETA only
- norm_parts(1) = sum( mvec(1 : NLOCAL)**2 / cov_model(1 : NLOCAL) )
+ norm_parts(1) = sum( mtemp(1 : NLOCAL)**2 / cov_model(1 : NLOCAL) )
! source part of the norm -- only events that you are inverting for
do ievent = ievent_min, ievent_max
@@ -1153,9 +1163,9 @@
itemp2 = NLOCAL + index_source(2,ievent)
itemp3 = NLOCAL + index_source(3,ievent)
- norm_parts(2) = norm_parts(2) + mvec(itemp1)**2 / cov_model(itemp1)
- norm_parts(3) = norm_parts(3) + mvec(itemp2)**2 / cov_model(itemp2)
- norm_parts(4) = norm_parts(4) + mvec(itemp3)**2 / cov_model(itemp3)
+ 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)
enddo
norm_struct = norm_parts(1)
More information about the CIG-COMMITS
mailing list