[cig-commits] r16258 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: PLOTTING matlab matlab/matlab_scripts src
carltape at geodynamics.org
carltape at geodynamics.org
Fri Feb 12 12:17:29 PST 2010
Author: carltape
Date: 2010-02-12 12:17:28 -0800 (Fri, 12 Feb 2010)
New Revision: 16258
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/gji_subspace_initial.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_diff_vec.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gnorm_sq.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/axes_expand.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/ridge_carl.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m
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_sub.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_variables.f90
Log:
Added option for Matlab-based conjugate gradient inversion using a full covariance matrix.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/gji_subspace_initial.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/gji_subspace_initial.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/gji_subspace_initial.pl 2010-02-12 20:17:28 UTC (rev 16258)
@@ -0,0 +1,308 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# gji_subspace_initial.pl
+# Carl Tape
+# 15-Jan-2010
+#
+# Copied from plot_surf_model_all.pl
+# This plots mprior, m00, mtar, and ln(mtar/m00)
+#
+# Example:
+# ../gji_subspace_initial.pl 7200 0/0 1/0/5 0.1/0.05 0
+#
+#==========================================================
+
+if (@ARGV < 5) {die("Usage: gji_subspace_initial.pl xxx\n");}
+($irun0,$imodels,$ievents,$pbar,$ifinite) = @ARGV;
+$comp = 1;
+
+# USER INPUT
+$ifig1 = 1;
+$ifig2 = 0;
+$itest = 0; # plot CG test models or not
+
+$icolor = 1;
+$ijpg = 0;
+$ipdf = 1;
+
+# base directory
+$pwd = $ENV{PWD};
+$dirplot = `dirname $pwd`; chomp($dirplot);
+$basedir = `dirname $dirplot`; chomp($basedir);
+#$basedir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work";
+
+($imodelmin,$imodelmax) = split("/",$imodels);
+($pmax,$ptick) = split("/",$pbar);
+($ievent_all,$ievent_one,$ievent) = split("/",$ievents);
+$rec_file1 = "sr.txt";
+$efile_syn = "events_syn_lonlat.dat";
+$mfile_syn = "structure_syn.dat";
+$mfile_dat = "structure_dat.dat";
+
+# directories
+$plotdir = "${basedir}/PLOTTING";
+$odir = "${basedir}/OUTPUT";
+$idir = "${basedir}/INPUT";
+$figdir = "${plotdir}/FIGURES";
+if (not -e $figdir) {die("Check if figdir $figdir exist or not\n");}
+
+$stirun0 = sprintf("%4.4i",$irun0);
+
+$dir0 = "$odir/run_${stirun0}";
+$dir2 = sprintf("$dir0/event_%3.3i",$ievent);
+$rec_file = "$dir2/${rec_file1}";
+$evefile = "$dir0/${efile_syn}"; # event file from the base directory
+if (not -f $evefile) { die("Check if $evefile exist or not\n") }
+
+# read in reference values: alpha0, beta0, per
+if (1==1) {
+ $fileref1 = "$dir0/reference_values.dat";
+ $fileref2 = "$dir0/reference_period.dat";
+ open(IN,$fileref1); @lines = <IN>; ($alpha0,$beta0) = split(" ",$lines[0]);
+ open(IN,$fileref2); $per = <IN>; chomp($per);
+} else {
+ $per = 20;
+ $beta0 = 3500;
+}
+$per_lab = sprintf("T = %3.3f s",$per);
+$beta0_lab = sprintf("beta0 = %3.3f km/s",$beta0/1000);
+#print "\n-- $per -- $per_lab -- $beta0_lab --\n"; die("testing");
+
+$plabel = "${plotdir}/gji_subspace_initial.pl";
+
+#die("\ntesting\n");
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "11";
+$fsize2 = "9";
+$fsize3 = "5";
+$fontno = "1"; # 1 or 4
+$tick = "0.1c";
+
+$tinc = 50; # tick increment (in seconds) for source time function
+
+# plot symbols for sources, receivers, and shelf
+if($ifinite==0) {$src0 = "-W0.75p -Sa0.20 -G255/0/0"} else {$src0 = "-Sc0.05"}
+$src = "-W0.75p -Sa0.20";
+$rec = "-W0.5p/0/0/0 -St0.10";
+$rec0 = "-Sc10p -W0.5p";
+$Wshelf = "-W1.0/0/0/0tap";
+
+# resolution of color plots
+$interp = "-I0.5m/0.5m -S4m"; # key information
+$interp = "-I2m/2m -S4m"; # key information
+$grdfile = "temp.grd";
+
+#-------------------------
+# color
+
+# KEY: scaling for color
+$scale_color = 21.0;
+$colorbar = "seis";
+
+#-------------------------
+
+$wid1 = 2.5;
+$wid2 = 2.0;
+
+$B1dat = "-B1:.\" \":WeSN";
+$B1syn = "-B1:.\" \":wESN";
+
+#===============================================
+$cshfile = "gji_subspace_initial.csh";
+print "\nWriting to CSH file ${cshfile}...\n";
+open(CSH,">$cshfile");
+print CSH "gmtset BASEMAP_TYPE plain PAPER_MEDIA letter TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 PLOT_DEGREE_FORMAT D HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1\n";
+#===============================================
+
+# make colorpoint file
+# $T1 = "-T3/4/0.1";
+# $pmax = 10;
+$dc = $pmax/10;
+$T1 = "-T-$pmax/$pmax/$dc";
+$cptfile1 = "color0.cpt";
+print CSH "makecpt -C$colorbar -D $T1 > $cptfile1\n";
+
+#$rec_file = "$edir/sr.txt";
+if (not -f $rec_file) {die("Check if rec_file $rec_file exist or not\n");}
+
+# sides to show the lat-lon tick marks
+$B0 = "-B1:.\" \":";
+ at Bopts1 = ("WesN","wEsN","Wesn","wEsn","Wesn","wEsn");
+ at Bopts2 = ("WesN","wesN","wEsN","Wesn","wesn","wEsn","Wesn","wesn","wEsn");
+
+# figure 1
+$dX = 1.1*$wid1; $dY = 0; $shift1 = "-X$dX -Y$dY";
+$dX = -1.1*$wid1; $dY = -1.1*$wid1; $shift2 = "-X$dX -Y$dY";
+ at shifts1 = ($origin,$shift1,$shift2,$shift1,$shift2,$shift1);
+
+# figure 2
+$dX = 1.1*$wid2; $dY = 0; $shift1 = "-X$dX -Y$dY";
+$dX = -2.2*$wid2; $dY = -1.1*$wid2; $shift2 = "-X$dX -Y$dY";
+ at shifts2 = ($origin,$shift1,$shift1,$shift2,$shift1,$shift1,$shift2,$shift1,$shift1);
+
+$fcol = 6;
+#@finds = (7,7,7,7,7,7,7,7,7);
+
+$J_title = "-JM1";
+$R_title = "-R0/1/0/1";
+
+#================
+
+# loop over models
+for ($m = $imodelmin; $m <= $imodelmax; $m = $m+1) {
+
+ $imodel = $m; # current model
+ $mtest = $m % 2; # =1 for a test model
+
+ $stirun = sprintf("%4.4i",$irun0+$imodel);
+ if ($mtest==1) {
+ $stmod = sprintf("m%2.2it",$imodel/2);
+ } else {
+ $stmod = sprintf("m%2.2i",$imodel/2);
+ }
+
+ # directories (current model and previos model)
+ $stirun = sprintf("%4.4i",$irun0+$imodel);
+ $dir1 = "$odir/run_${stirun}";
+
+ $name = "gji_subspace_initial";
+ $psfile = "$name.ps";
+ $jpgfile = "$name.jpg";
+
+# write plotting scripts
+$wid = $wid1;
+$J1 = "-JM${wid}i"; # in lat-lon
+$origin = "-X1.5 -Y7.25";
+$Dlen = 2.0; $Dx = $wid/2; $Dy = -0.35;
+$Dscale1 = "-D$Dx/$Dy/$Dlen/0.15h";
+#$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+#$Bscale2 = sprintf("-B%2.2e:\" Perturbation from m00 \": -E10p",$ptick);
+$Bscale = sprintf("-B%2.2e:\" \": -E10p",$ptick);
+$xlab = $wid/2;
+$ylab = -0.02*$wid;
+$olab = "-Xa$xlab -Ya$ylab";
+
+ @files = ("$dir0/${mfile_syn}","$dir0/${mfile_syn}",
+ "$dir1/${mfile_dat}","$dir0/${mfile_dat}");
+ @stlabs = ("(a) mprior","(b) $stmod","(c) mtar","(d) ln(mtar/$stmod)");
+
+ # only plot even models (m04 = iteration 8)
+ $kmaxt = 4;
+ $kmaxm = 4*(1 - $mtest); # =0 for test model
+ if ($itest == 1) {
+ $kmax = $kmaxt;
+ } else {
+ $kmax = $kmaxm;
+ }
+
+ # loop over subplots
+ for ($k = 0; $k < $kmax; $k = $k+1) {
+
+ $iB = $Bopts1[$k];
+ $shift = $shifts1[$k];
+ $B = "${B0}${iB}";
+ $find = $fcol;
+ $find2 = 2*$fcol;
+
+ $file0 = $files[0]; # initial model
+ $file1 = $files[$k];
+ if (not -f $file1) {
+ die("Check if $file1 exist or not\n");
+ }
+
+ # phase velocity map
+ if ($k==0 && $m==$imodelmin) {
+ # set bounds for the plotting
+ #$xmin = -121; $xmax = -115.0; $zmin = 31.75; $zmax = 36.5;
+ ($xmin,$xmax,$zmin,$zmax,$smin,$smax,$tmin,$tmax) = split(" ",`minmax -C $file1`);
+ $dinc = 0.25; # buffer around min/max
+ $xmin = $xmin-$dinc; $xmax = $xmax+$dinc;
+ $zmin = $zmin-$dinc; $zmax = $zmax+$dinc;
+ $R = "-R$xmin/$xmax/$zmin/$zmax";
+ print "\n$R\n";
+ }
+
+ if ($k==0) {
+ print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile\n"; # START
+ } else {
+ print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile\n";
+ }
+
+ if ($icolor == 1) {
+ #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ if ($k == 0) {
+ print CSH "awk '{print \$1,\$2,0}' $file0 | nearneighbor -G$grdfile $R $interp\n";
+ } elsif ($k == 1) {
+ print CSH "awk '{print \$1,\$2,\$${find}}' $file0 | nearneighbor -G$grdfile $R $interp\n";
+ } elsif ($k == 2) {
+ print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
+ } elsif ($k == 3) {
+ $file2 = "ftemp";
+ print CSH "paste $file1 $file0 > $file2\n";
+ print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
+ } else {
+ die("k ($k) must be 0,1,2,3");
+ }
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ }
+ print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+
+ # plot receivers
+ if (0==1) { # numbered large circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N $J1 $R -K -O -V $rec0 >> $psfile\n";
+ $rec_file2 = text_rec; $angle = 0; $just = "CM";
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3,$fsize3,$angle,$fontno,\"$just\",\$4}' $rec_file > $rec_file2\n";
+ print CSH "pstext $rec_file2 -N $J1 $R -K -O -V >> $psfile\n";
+
+ } else { # small circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile\n";
+
+ }
+
+ # plot sources
+ if ($ievent_one) {
+ print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile\n";
+ }
+ if ($ievent_all) {
+ print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile\n";
+ }
+
+ # plot label
+ $stlab = $stlabs[$k];
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 12 0 $fontno CM $stlab\nEOF\n";
+
+ # colorscale bars
+ #if ($k==2) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile\n";}
+ #if ($k==3) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile\n";}
+ if ($k==2 || $k==3) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale -K -O -V >> $psfile\n";}
+
+ # plot title and GMT header
+ if ($k==3) {
+ $plot_title = "Model m00, model $stmod, and target model (run_${stirun})";
+ $Utag = "-U/0/0.25/$plabel";
+ $shift = "-Xa-3.5 -Ya8.5";
+ print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
+
+ if ($ijpg == 1) {
+ print CSH "convert $psfile $jpgfile\n";
+ }
+ if ($ipdf == 1) {
+ print CSH "ps2pdf $psfile\n";
+ }
+
+ }
+ } # loop over subplots
+
+} # loop over models
+
+close (CSH);
+system("csh -f $cshfile");
+if($ifig1==1) {system("gv $psfile &")}
+if($ifig2==1) {system("gv $psfile2 &")}
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/gji_subspace_initial.pl
___________________________________________________________________
Name: svn:executable
+ *
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-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2010-02-12 20:17:28 UTC (rev 16258)
@@ -83,7 +83,7 @@
$icolor = 1; # ccc
$irun = $irun0_cg + $iter;
- at mods = ("0","0t","1","1t","2","2t","3","3t","4","4t","5","5t","6","6t","7","7t","8","8t","9","9t","10","10t","11","11t","12","12t","13","13t","14","14t","15","15t","16","16t");
+ at mods = ("0","0t","1","1t","2","2t","3","3t","4","4t","5","5t","6","6t","7","7t","8","8t","9","9t","10","10t","11","11t","12","12t","13","13t","14","14t","15","15t","16","16t","17","17t","18","18t","19","19t","20","20t","21","21t","22","22t","23","23t","24","24t","25","25t","26","26t","27","27t","28","28t","29","29t","30","30t","31","31t","32");
$mod = $mods[$iter];
$smod = "m\@+$mod\@+";
@@ -320,6 +320,7 @@
# scale for chi-vs-m plots (a1f2g1p)
$iter_tick = 2;
+ if ($a2x2 > 16) {$iter_tick = 4}
#$B3a = "-B${iter_tick}f1g2:\" k, model number \":/a1f3g3p:\" $misfitvar ( m ) ( $utype ) \":";
$B3a = "-B${iter_tick}f1g2:\" k, model number \":/a1f3g3p:\" $misfitvar ( m )\":";
$B3b = "-B${iter_tick}f1g2:\" k, model number \":/a1f3g3p:\" \":";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl 2010-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl 2010-02-12 20:17:28 UTC (rev 16258)
@@ -13,6 +13,7 @@
#
# Example:
# ../plot_surf_model_all.pl 800 0/32 0/1/5 0.1/0.05 0
+# ../plot_surf_model_all.pl 7200 64/64 1/0/5 0.1/0.05 0
#
#==========================================================
@@ -23,7 +24,7 @@
# USER INPUT
$ifig1 = 1;
$ifig2 = 1;
-$itest = 1; # plot CG test models or not
+$itest = 0; # plot CG test models or not
$icolor = 1;
$ijpg = 0;
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/axes_expand.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/axes_expand.m 2010-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/axes_expand.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -1,21 +1,19 @@
%
-% function ax1 = axes_expand(ax0,fac)
-% CARL TAPE, 16-Aug-2005
-% printed xxx
+% function ax1 = axes_expand(ax0,fac,iopt)
+% Carl Tape, 16-Aug-2005
%
-% This function inputs an axes box and outputs a new axes box,
-% expanded by the factor 'fac' in all 2,4,6 dimensions.
+% This function inputs an axes box and outputs a new expanded axes box.
%
% EXAMPLE:
-% ax0 = [-121 -114]; ax1 = axes_expand(ax0,2)
-% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,2)
-% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,0.30)
+% ax0 = [-121 -114]; ax1 = axes_expand(ax0,1.2,0)
+% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,1.2,0)
+% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,0.30,0)
%
% calls xxx
% called by xxx
%
-function ax1 = axes_expand(ax0,fac)
+function ax1 = axes_expand(ax0,fac,iopt)
% 1D, 2D, 3D
ndim = length(ax0)/2;
@@ -34,16 +32,18 @@
ymin = ax0(3);
ymax = ax0(4);
dy = ymax-ymin;
- ax1(3) = ymin - dy*(fac-1);
- ax1(4) = ymax + dy*(fac-1);
+ if iopt == 1, dinc = dy; else dinc = dx; end
+ ax1(3) = ymin - dinc*(fac-1);
+ ax1(4) = ymax + dinc*(fac-1);
if ax1(4) <= ax1(3), ax1 = ax0; end
end
if ndim == 3
zmin = ax0(3);
zmax = ax0(4);
dz = zmax-zmin;
- ax1(5) = zmin - dz*(fac-1);
- ax1(6) = zmax + dz*(fac-1);
+ if iopt == 1, dinc = dz; else dinc = dx; end
+ ax1(5) = zmin - dinc*(fac-1);
+ ax1(6) = zmax + dinc*(fac-1);
if ax1(6) <= ax1(5), ax1 = ax0; end
end
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/ridge_carl.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/ridge_carl.m 2010-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/ridge_carl.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -186,7 +186,7 @@
plot(log10(rssL),log10(mssL),'ko','markersize',8,'MarkerFaceColor','r');
plot(log10(rssG),log10(mssG),'kV','markersize',8,'MarkerFaceColor','g');
plot(log10(rssF),log10(mssF),'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1,1));
axy = axis; dx = axy(2)-axy(1);
for kk=1:nlab
ii = ilabs(kk);
@@ -200,7 +200,7 @@
plot(log10(lamL),kapL,'ko','markersize',8,'MarkerFaceColor','r');
plot(log10(lamG),kapG,'kV','markersize',8,'MarkerFaceColor','g');
plot(log10(lamF),kapF,'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1,1));
legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northwest');
xlabel(stx2); ylabel(sty2); grid on;
@@ -209,7 +209,7 @@
plot(log10(lamL),log10(FL),'ko','markersize',8,'MarkerFaceColor','r');
plot(log10(lamG),log10(FG),'kV','markersize',8,'MarkerFaceColor','g');
plot(log10(lamF),log10(FF),'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1,1));
legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northeast');
xlabel(stx3); ylabel(sty3); grid on;
@@ -218,7 +218,7 @@
plot(log10(lamL),log10(GL),'ko','markersize',8,'MarkerFaceColor','r');
plot(log10(lamG),log10(GG),'kV','markersize',8,'MarkerFaceColor','g');
plot(log10(lamF),log10(GF),'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1,1));
legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northwest');
xlabel(stx4); ylabel(sty4); grid on;
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m 2010-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -1,5 +1,5 @@
%
-% function
+% function [m,pk] = wave2d_cg(m0,gk,C,chi_k_val,mu_val,istep,imaketest,gfile,chitfile)
% Carl Tape, 21-Jan-2010
%
% This emulates the CG algorithm in wave2d.f90, but assumes a full
@@ -13,6 +13,8 @@
disp('------CG ALGORITHM---------');
+whos m0 gk C
+
nmod = length(m0);
mt = zeros(nmod,1);
mk = zeros(nmod,1);
@@ -29,8 +31,20 @@
if isinf(beta_val), error('beta_val is infinity'); end
end
pk = -C * gk + beta_val * p0;
-lam_t_val = 2.0*(mu_val - chi_k_val) / sum( gk .* pk ); % gk is hat, pk is non-hat
+lam_t_val_bot = sum( gk .* pk ); % gk is hat, pk is non-hat
+lam_t_val = 2.0*(mu_val - chi_k_val) / lam_t_val_bot;
+% check
+if 0==1
+ disp('checking values in the CG algorithm');
+ vals = [mu_val lam_t_val_bot lam_t_val beta_val chi_k_val];
+ %write(19,'(5e16.8)') mu_val, lam_t_val_bot, lam_t_val, beta_val, chi_k_val
+ dirbase = '/home/carltape/ADJOINT_TOMO/iterate_adj/';
+ vals0 = load([dirbase 'SEM2D_iterate_OUTPUT/run_9150/cg_test_vals.dat']);
+ for ii=1:length(vals), vals0(ii), vals(ii), end
+ error('testing here');
+end
+
if imaketest==1
% test model
mt = m0 + lam_t_val*pk;
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_diff_vec.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_diff_vec.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_diff_vec.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -0,0 +1,31 @@
+%
+% function
+% Carl Tape, 25-Jan-2010
+%
+%
+%
+% calls xxx
+% called by xxx
+%
+
+function wave2d_diff_vec(m1,m2,m_inds,m_labs,ifig)
+
+disp('wave2d_diff_vec: comparing two vectors:');
+
+for ii=1:4
+ inds = m_inds(ii,1):m_inds(ii,2);
+ disp(sprintf('%20s%10.3e%10.3e%10.3e%10.3e',m_labs{ii},...
+ norm(m1(inds)), norm(m2(inds)),...
+ norm(m1(inds)-m2(inds)),...
+ norm(m1(inds)-m2(inds)) / norm(m1(inds)) ))
+
+ if ifig==1
+ figure; nr=3; nc=1;
+ subplot(nr,nc,1); plot( m1(inds), '.'); title(stpars{ii});
+ subplot(nr,nc,2); plot( m2(inds), '.')
+ subplot(nr,nc,3); plot( m1(inds), m2(inds), '.')
+ end
+end
+
+
+%=========================================================
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gnorm_sq.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gnorm_sq.m 2010-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gnorm_sq.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -1,40 +0,0 @@
-%
-% function
-% Carl Tape, 25-Jan-2010
-%
-%
-%
-% calls xxx
-% called by xxx
-%
-
-function norm_parts = wave2d_gnorm_sq(g,C,m_inds)
-
-g = g(:);
-[a,b] = size(C);
-if a==1, error('check dimension of C'); end
-if b==1, icdiag = 1; end
-
-npar = length(m_inds);
-
-%stfmt = repmat('%16.4e',1,npar);
-
-disp('norm of parts:');
-norm_parts = zeros(npar,1);
-for ii=1:npar
- inds = m_inds(ii,1):m_inds(ii,2);
- if icdiag == 1
- norm_parts(ii) = sum( g(inds).^2 .* C(inds) );
- else
- norm_parts(ii) = g(inds)' * C(inds,inds) * g(inds);
- end
- disp(sprintf('%4i%16.4e',ii,norm_parts(ii)));
-end
-
-if 1==1
- disp(sprintf(' 1%16.4e',sum(norm_parts(1))));
- disp(sprintf(' 2-4%16.4e',sum(norm_parts(2:4))));
- disp(sprintf(' 1-4%16.4e',sum(norm_parts(1:4))));
-end
-
-%=========================================================
\ No newline at end of file
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m (from rev 16191, seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gnorm_sq.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -0,0 +1,127 @@
+%
+% function norm_parts = wave2d_norm_sq(m,mprior,C,m_inds,weights,imnorm)
+% Carl Tape, 25-Jan-2010
+%
+% Compute the norms of different parts of a model vector or gradient.
+%
+% calls xxx
+% called by xxx
+%
+
+function norm_parts = wave2d_norm_sq(m,mprior,C,m_inds,weights,imnorm)
+
+m = m(:);
+[a,b] = size(C);
+if a==1, error('check dimension of C'); end
+if b==1, icdiag = 1; else icdiag = 0; end
+
+npar = length(m_inds);
+
+%stfmt = repmat('%16.4e',1,npar);
+
+Cmat = zeros(a,b);
+mvec = zeros(length(m),1);
+npw = zeros(npar,1);
+if imnorm==0
+ Cmat = C;
+ mvec = m;
+ npw = 1 ./ weights;
+
+elseif imnorm==1
+ if icdiag==1, Cmat = 1./C; else Cmat = inv(C); end
+ mvec = m - mprior;
+ npw = weights;
+end
+
+disp('norm of parts:');
+norm_parts = zeros(npar,1);
+for ii=1:npar
+ inds = m_inds(ii,1):m_inds(ii,2);
+ if icdiag == 1
+ norm_parts(ii) = sum( mvec(inds).^2 .* Cmat(inds) );
+ else
+ norm_parts(ii) = mvec(inds)' * Cmat(inds,inds) * mvec(inds);
+ end
+ disp(sprintf('%4i%16.4e%16.4e',ii,norm_parts(ii),sqrt(norm_parts(ii)) ));
+end
+
+if 1==1
+ disp(sprintf(' 1%16.4e',sum(norm_parts(1))));
+ disp(sprintf(' 2-4%16.4e',sum(norm_parts(2:4))));
+ disp(sprintf(' 1-4%16.4e',sum(norm_parts(1:4))));
+end
+
+%=========================================================
+
+% subroutine compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+% mvec, mvec_prior, cov_model, norm_sq_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.
+% ! 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) :: imnorm, ievent_min, ievent_max, nevent, nmod
+% integer, dimension(NVAR_SOURCE, nevent), intent(in) ::index_source
+% 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, dimension(NVAR), intent(out) :: norm_sq_parts
+%
+% double precision, dimension(nmod) :: mtemp, ctemp
+% double precision, dimension(NVAR) :: npw
+% integer :: ievent, itemp1, itemp2, itemp3
+%
+% !----------
+%
+% norm_sq_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
+% ! NOTE 3: norm_parts_weight is related to the norm-squared weights (not norm)
+% if (imnorm == 0) then
+% ctemp(:) = 1.0 / cov_model(:)
+% mtemp(:) = mvec(:)
+% if (present(norm_parts_weight)) npw(:) = 1.0 / norm_parts_weight(:)
+% !if (present(norm_parts_weight)) npw(:) = 1.0 / norm_parts_weight(:)**2
+%
+% elseif (imnorm == 1) then
+% ctemp(:) = cov_model(:)
+% mtemp(:) = mvec(:) - mvec_prior(:)
+% if (present(norm_parts_weight)) npw(:) = norm_parts_weight(:)
+% !if (present(norm_parts_weight)) npw(:) = norm_parts_weight(:)**2
+%
+% else
+% stop 'imnorm must = 0 or 1'
+% endif
+%
+% ! structure part of the norm -- BETA only
+% ! NOTE: division corresponds to inversion of a diagonal covariance matrix
+% norm_sq_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 = NLOCAL + index_source(1,ievent)
+% itemp2 = NLOCAL + index_source(2,ievent)
+% itemp3 = NLOCAL + index_source(3,ievent)
+%
+% norm_sq_parts(2) = norm_sq_parts(2) + mtemp(itemp1)**2 / ctemp(itemp1)
+% norm_sq_parts(3) = norm_sq_parts(3) + mtemp(itemp2)**2 / ctemp(itemp2)
+% norm_sq_parts(4) = norm_sq_parts(4) + mtemp(itemp3)**2 / ctemp(itemp3)
+% enddo
+%
+% !norm_struct = norm_sq_parts(1)
+% !norm_source = sum( norm_sq_parts(2:4) )
+% !norm_total = sum( norm_sq_parts(1:4) )
+%
+% ! weight each part of the norm
+% norm_sq_parts(:) = norm_sq_parts(:) * npw(:)
+%
+% end subroutine compute_norm_sq
+%
+%
\ No newline at end of file
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m
___________________________________________________________________
Name: svn:mergeinfo
+
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-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -14,7 +14,7 @@
% axes limits
n1 = 0;
n2 = 8;
-n2 = 17;
+n2 = max([17 xlims(2)]);
%ax_chi = [n1-1 n2 10^0 axpoly(4)];
%ax_chi = [n1-1 n2 10.^[0 4]];
ax_chi = [n1-1 n2 10.^[-1 3] ];
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-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -60,7 +60,7 @@
% determine the number of iterations
k = 0;
-while k < 50
+while k < 100
irun = irun0 + k;
chifile = [dir0 odir 'run_' sprintf(stfm,irun) '/chi.dat'];
iexist = exist(chifile);
@@ -180,7 +180,8 @@
axi = xlim;
plot([lam0 lam0 axi(1)],[0 chi chi],'k');
plot(lam0,chi,'bo','markersize',msize,'MarkerFaceColor','b');
- if(lam0 > axi(2)), axis tight; end
+ %axis tight; ax1 = axis; axis([ax1(1:2) 0 ax1(4)]);
+ if lam0 > axi(2), axis tight; end
%end
end
@@ -364,8 +365,10 @@
% variance reduction
chi_before = chis(1:nchi-1);
chi_after = chis(2:nchi);
-var_red1 = 100 * ( 1 - ( (chi_after-chi_data_stop).^2 ./ (chi_before-chi_data_stop).^2 ) );
-var_red2 = 100 * ( 1 - ( (chi_after-chi_data_stop).^2 ./ (chi_before-chi_data_stop).^2 ) );
+%var_red1 = 100 * ( 1 - ( (chi_after-chi_data_stop).^2 ./ (chi_before-chi_data_stop).^2 ) );
+%var_red2 = 100 * ( 1 - ( (chi_after-chi_data_stop).^2 ./ (chi_before-chi_data_stop).^2 ) );
+var_red1 = log( chi_before ./ chi_after );
+var_red2 = log( chi_before ./ chi_after );
disp(' '); disp('VARIANCE REDUCTION');
disp([its(2:end) chi_before chi_after var_red1]);
@@ -373,7 +376,8 @@
% var red corresponding to best-fitting chi-vs-its curve
its2 = xticks;
chifit2 = interp1(its_smooth,chifit_smooth,its2);
-var_red1_fit = 100 * ( 1 - ( chifit2(2:end).^2 ./ chifit2(1:end-1).^2 ) );
+%var_red1_fit = 100 * ( 1 - ( chifit2(2:end).^2 ./ chifit2(1:end-1).^2 ) );
+var_red1_fit = log( chifit2(1:end-1) ./ chifit2(2:end) )
x_var = its(2:end);
x_var_fit = its2(2:end);
@@ -381,7 +385,7 @@
subplot(nr,nc,3); hold on;
plot(x_var_fit,var_red1_fit,'r.--','markersize',msize,'linewidth',lsize);
plot(x_var,var_red1,'b.-','markersize',msize,'linewidth',lsize);
-axis([xlims 0 100]); grid on; set(gca,'xtick',xticks);
+axis([xlims 0 1]); grid on; set(gca,'xtick',xticks);
xlabel(stx);
ylabel(' variance reduction');
title(' variance reduction between successive models');
@@ -471,7 +475,7 @@
figure; hold on;
specs = [2 14 18 12]; fac = 0.05;
- axpoly = axes_expand([x1 x2 0 max([y1 y2])],1.2);
+ axpoly = axes_expand([x1 x2 0 max([y1 y2])],1.2,1);
axpoly(3) = 0;
dy = axpoly(4) - axpoly(3);
dx = axpoly(2) - axpoly(1);
@@ -512,7 +516,7 @@
figure; hold on;
specs = [2 14 18 12]; fac = 0.05;
- axpoly = axes_expand([x1 x2 0 max([y1 y2])],1.2);
+ axpoly = axes_expand([x1 x2 0 max([y1 y2])],1.2,1);
axpoly(3) = 0;
dy = axpoly(4) - axpoly(3);
dx = axpoly(2) - axpoly(1);
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-02-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_run.m 2010-02-12 20:17:28 UTC (rev 16258)
@@ -20,7 +20,8 @@
colors;
stfm = '%4.4i';
-stpars = {'B = ln(beta/beta0)','Ts = ts = ts0','Xs = xs - xs0','Ys = ys - ys0'};
+%stpars = {'B = ln(beta/beta0)','Ts = ts = ts0','Xs = xs - xs0','Ys = ys - ys0'};
+stpars = {'B = ln(beta/beta0)','ts','xs','ys'};
npts = 100;
%---------------------------------------------------------
@@ -36,16 +37,16 @@
if ~exist(dirrand,'dir'), error(['dirrand ' dirrand ' does not exist']); end
% KEY
-parms = [2100 0];
+parms = [1000 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 = 0;
+iwrite = 1;
iplotmod = 1;
-iplotker = 0;
+iplotker = 1;
sticg = sprintf('%2.2i',icg);
%---------------------------------------------------------
@@ -116,23 +117,16 @@
% constants for model covariance
[stnames,stvals] = textread([idir1 'scaling_values_covm.dat'],'%s%f');
-for ii=1:length(stnames), eval([stnames{ii} ' = stvals(' num2str(ii) ')']); end
-% fac_str = stvals(1);
-% fac_ts = stvals(2);
-% fac_xs = stvals(3);
-% fac_ys = stvals(4);
-% fac_total = stvals(5);
-% ugrad_str = stvals(6);
-% ugrad_ts = stvals(7);
-% ugrad_xs = stvals(8);
-% ugrad_ys = stvals(9);
-% dnparm_src_run = stvals(10);
-% coverage_str = stvals(11);
-% coverage_src = stvals(12);
-% cov_imetric_fac_str = stvals(13);
-% cov_imetric_fac_ts = stvals(14);
-% cov_imetric_fac_xs = stvals(15);
-% cov_imetric_fac_ys = stvals(16);
+%for ii=1:length(stnames), eval([stnames{ii} ' = stvals(' num2str(ii) ')']); end
+covm_weight_parts = stvals(1:4);
+covg_weight_parts = stvals(6:9);
+ugsq_str = stvals(11);
+ugsq_ts = stvals(12);
+ugsq_xs = stvals(13);
+ugsq_ys = stvals(14);
+dnevent_run = stvals(15);
+coverage_str = stvals(16);
+coverage_src = stvals(17);
% constants for data covariance
[stnames,stvals] = textread([idir1 'scaling_values_covd.dat'],'%s%f');
@@ -193,26 +187,27 @@
%[m_str_lon,m_str_lat,m_str_kappa,m_str_mu,m_str_rho,m_str_B] = ...
% textread([idir1 'structure_syn.dat'],'%f%f%f%f%f%f%f');
-% INITIAL model vector (m00)
-%write(88,'(2e24.12)') m0_vec(i), m0(i)
-[m0_vec_initial, m0_initial] = textread([dir0 'initial_model_vector.dat'],'%f%f');
+% prior, initiral, and target models
+[mprior, m0_initial, mtarget] = textread([dir0 'prior_initial_target_models.dat'],'%f%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');
%[m0_B, m0_T, m0_X, m0_Y] = wave2d_splitm(m0,m_inds);
%[m0_vec_B, m0_vec_T, m0_vec_X, m0_vec_Y] = wave2d_splitm(m0_vec,m_inds);
+% initial
+
% load kernels
kernels_all = zeros(nlocal,nevent_run);
for ii=1:nevent_run
ievent = einds(ii)
dirK = [idir1 'event_' sprintf('%3.3i',ievent) '/'];
- %kernel = load([dirK 'kernel_basis']); Kbeta = kernel(:,7);
- kernel = load([dirK 'kernel_basis_smooth']); Kbeta = kernel(:,3);
+ if icg==3
+ kernel = load([dirK 'kernel_basis']); Kbeta = kernel(:,7);
+ else
+ kernel = load([dirK 'kernel_basis_smooth']); Kbeta = kernel(:,3);
+ end
kernels_all(:,ii) = Kbeta;
if iplotker==1
@@ -224,8 +219,8 @@
summed_ker = sum(kernels_all,2);
% source gradient
-%source_gradient = load([idir1 'source_gradient.dat']);
-source_gradient = zeros(nmod_src,1);
+source_gradient = load([idir1 'source_gradient.dat']);
+%source_gradient = zeros(nmod_src,1);
% data covariance
% nmeas_run = nevent_run * nrec * NCOMP
@@ -286,17 +281,18 @@
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;
- cov_model(indX) = sigma_xs^2 * dnparm_src_run * coverage_src;
- cov_model(indY) = sigma_zs^2 * dnparm_src_run * coverage_src;
+ cov_model(indB) = sigma_beta^2 * Atot ./ da_local_vec * coverage_str;
+ cov_model(indT) = sigma_ts^2 * dnevent_run * coverage_src;
+ cov_model(indX) = sigma_xs^2 * dnevent_run * coverage_src;
+ cov_model(indY) = sigma_zs^2 * dnevent_run * coverage_src;
cov_imetric = zeros(nmod,1);
- cov_imetric(indB) = cov_model(indB) * (fac_str / ugrad_str) * fac_total;
- cov_imetric(indT) = cov_model(indT) * (fac_ts / ugrad_ts) * fac_total;
- cov_imetric(indX) = cov_model(indX) * (fac_xs / ugrad_xs) * fac_total;
- cov_imetric(indY) = cov_model(indY) * (fac_ys / ugrad_ys) * fac_total;
+ cov_imetric(indB) = cov_model(indB) / covg_weight_parts(1);
+ cov_imetric(indT) = cov_model(indT) / covg_weight_parts(2);
+ cov_imetric(indX) = cov_model(indX) / covg_weight_parts(3);
+ cov_imetric(indY) = cov_model(indY) / covg_weight_parts(4);
else
+ % ONLY FOR TESTING THE COMPARISON
cov_model = zeros(nmod,1);
cov_model(indB) = sigma_beta^2;
cov_model(indT) = sigma_ts^2;
@@ -304,79 +300,86 @@
cov_model(indY) = sigma_zs^2;
cov_weight = zeros(nmod,1);
- %cov_weight(indB) = (fac_str / ugrad_str) * fac_total ./ da_local_vec * Atot * coverage_str;
- cov_weight(indB) = (fac_str / ugrad_str) * fac_total * coverage_str ./ Avec.^2;
- cov_weight(indT) = (fac_ts / ugrad_ts) * fac_total * dnparm_src_run * coverage_src;
- cov_weight(indX) = (fac_xs / ugrad_xs) * fac_total * dnparm_src_run * coverage_src;
- cov_weight(indY) = (fac_ys / ugrad_ys) * fac_total * dnparm_src_run * coverage_src;
+ %cov_weight(indB) = coverage_str / covg_weight_parts(1) * Atot ./ da_local_vec;
+ cov_weight(indB) = coverage_str / covg_weight_parts(1) ./ Avec.^2;
+ cov_weight(indT) = dnevent_run * coverage_src / covg_weight_parts(2);
+ cov_weight(indX) = dnevent_run * coverage_src / covg_weight_parts(3);
+ cov_weight(indY) = dnevent_run * coverage_src / covg_weight_parts(4);
cov_imetric = zeros(nmod,1);
cov_imetric(indB) = cov_model(indB) .* cov_weight(indB);
cov_imetric(indT) = cov_model(indT) .* cov_weight(indT);
cov_imetric(indX) = cov_model(indX) .* cov_weight(indX);
cov_imetric(indY) = cov_model(indY) .* cov_weight(indY);
+
+ wave2d_diff_vec(cov_imetric0,cov_imetric,m_inds,stpars,0);
+ break
end
% cov_model -- TO CHECK
- %write(19,'(4e20.10)') cov_imetric(i), icov_metric(i), cov_model(i), icov_model(i)
- cov_model_all = load([idir1 'cov_model_diagonal.dat']);
- cov_imetric0 = cov_model_all(:,1);
+ [cov_model0,cov_imetric0] = textread([idir1 'cov_model_diagonal.dat'],'%f%f');
% check
- for ii=1:4
- inds = m_inds(ii,1):m_inds(ii,2);
- disp('------');
- disp(sprintf('%.3e %.3e %.3e %.3e',...
- norm(cov_imetric0(inds)), norm(cov_imetric(inds)),...
- norm(cov_imetric0(inds)-cov_imetric(inds)),...
- norm(cov_imetric0(inds)-cov_imetric(inds)) / norm(cov_imetric0(inds))))
- end
+ disp('comparing cov_model and cov_imetric');
+ wave2d_diff_vec(cov_model0,cov_model,m_inds,stpars,0);
+ wave2d_diff_vec(cov_imetric0,cov_imetric,m_inds,stpars,0);
% gradient
%gradient_model(:) = m0(:) / cov_imetric(:)
- gradient_tot = zeros(nmod,1);
- gradient_model = m0 ./ cov_imetric;
- gradient_data = zeros(nmod,1);
- gradient_data(indB) = summed_ker .* da_local_vec;
- gradient_data(indTXY) = source_gradient;
+ gradient_tot = zeros(nmod,1);
+ gradient_model = (m0-mprior) ./ cov_imetric;
+ %gradient_model = (m0-mprior) ./ cov_model;
+ gradient_data = zeros(nmod,1);
+ gradient_data(indB) = summed_ker .* da_local_vec;
+ gradient_data(indTXY) = source_gradient;
+ % zero out parts
+ if INV_SOURCE_T==0
+ gradient_data(indT) = 0;
+ gradient_model(indT) = 0;
+ end
+ if INV_SOURCE_X==0
+ gradient_data(indX) = 0;
+ gradient_model(indX) = 0;
+ gradient_data(indY) = 0;
+ gradient_model(indY) = 0;
+ end
+ if INV_STRUCT_BETA==0
+ gradient_data(indB) = 0;
+ gradient_model(indB) = 0;
+ end
gradient_tot = gradient_data + gradient_model;
% check actual norms
% --> gradient_norm_all.dat, gradient_norm_data_all.dat, gradient_norm_model all.dat
- wave2d_gnorm_sq(gradient_tot,cov_imetric,m_inds);
- wave2d_gnorm_sq(gradient_data,cov_imetric,m_inds);
- wave2d_gnorm_sq(gradient_model,cov_imetric,m_inds);
+ %wave2d_gnorm_sq(gradient_tot,cov_imetric,m_inds);
+ %wave2d_gnorm_sq(gradient_data,cov_imetric,m_inds);
+ %wave2d_gnorm_sq(gradient_model,cov_imetric,m_inds);
+ wave2d_norm_sq(gradient_tot,mprior,cov_imetric,m_inds,ones(4,1),0);
+ wave2d_norm_sq(gradient_data,mprior,cov_imetric,m_inds,ones(4,1),0);
+ wave2d_norm_sq(gradient_model,mprior,cov_imetric,m_inds,ones(4,1),0);
% gradient -- TO CHECK
%write(19,'(3e20.10)') gradient(i), gradient_data(i), gradient_model(i)
[gradient_tot0,gradient_data0,gradient_model0] = textread([idir1 'gradient_vec.dat'],'%f%f%f');
% check
- for ii=1:4
- inds = m_inds(ii,1):m_inds(ii,2);
- disp(sprintf('%20s %.3e %.3e %.3e %.3e',stpars{ii},...
- norm(gradient_tot0(inds)), norm(gradient_tot(inds)),...
- norm(gradient_tot0(inds)-gradient_tot(inds)),...
- norm(gradient_tot0(inds)-gradient_tot(inds)) / norm(gradient_tot0(inds)) ))
- figure; nr=3; nc=1;
- subplot(nr,nc,1); plot( gradient_tot0(inds), '.'); title(stpars{ii});
- subplot(nr,nc,2); plot( gradient_tot(inds), '.')
- subplot(nr,nc,3); plot( gradient_tot0(inds), gradient_tot(inds), '.')
- end
-
+ disp('comparing gradient, gradient_model, and gradient_data');
+ wave2d_diff_vec(gradient_tot0,gradient_tot,m_inds,stpars,0);
+ wave2d_diff_vec(gradient_data0,gradient_data,m_inds,stpars,0);
+ wave2d_diff_vec(gradient_model0,gradient_model,m_inds,stpars,0);
+
% check chi model values
model_norm_parts = zeros(4,1);
- 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 );
+ model_norm_parts(1) = sum( (m0(indB)-mprior(indB)).^2 ./ cov_model(indB) * covm_weight_parts(1) );
+ model_norm_parts(2) = sum( (m0(indT)-mprior(indT)).^2 ./ cov_model(indT) * covm_weight_parts(2) );
+ model_norm_parts(3) = sum( (m0(indX)-mprior(indX)).^2 ./ cov_model(indX) * covm_weight_parts(3) );
+ model_norm_parts(4) = sum( (m0(indY)-mprior(indY)).^2 ./ cov_model(indY) * covm_weight_parts(4) );
+ model_norm = sum(model_norm_parts);
+ %model_norm = sum( (m0-mprior).^2 ./ cov_model );
% 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
@@ -403,34 +406,29 @@
pk = -cov_imetric .* gk + beta_val * p0;
% step for test model
mu_val = chi_data_stop;
- lam_t_val = 2.0*(mu_val - chi_k_val) / sum( gk .* pk );
+ lam_t_val_bot = sum( gk .* pk );
+ lam_t_val = 2.0*(mu_val - chi_k_val) / lam_t_val_bot;
- % % check
- % if 0==1
- % norm(p0), norm(g0), norm(pk)
- % vals = [mu_val lam_t_val beta_val chi_k_val];
- % %write(19,'(4e16.8)') mu_val, lam_t_val, beta_val, chi_k_val
- % vals0 = load([dirbase 'SEM2D_iterate_OUTPUT/run_1204/cg_test_vals.dat']);
- % for ii=1:length(vals), vals0(ii), vals(ii), end
- % end
+ % check
+ if 0==1
+ disp('checking values in the CG algorithm');
+ vals = [mu_val lam_t_val_bot lam_t_val beta_val chi_k_val];
+ %write(19,'(5e16.8)') mu_val, lam_t_val_bot, lam_t_val, beta_val, chi_k_val
+ vals0 = load([dirbase 'SEM2D_iterate_OUTPUT/run_9150/cg_test_vals.dat']);
+ for ii=1:length(vals), vals0(ii), vals(ii), end
+ break
+ end
if imaketest==1
% test model
mt = m0 + lam_t_val*pk;
mt_vec = zeros(nmod,1);
- if (INV_STRUCT_BETA == 0)
- mt_vec(1:nlocal) = m0_vec_initial(1:nlocal);
- else
- mt_vec(1:nlocal) = beta0 * exp( mt(1:nlocal) );
- end
+ %mt_vec(1 : NLOCAL) = beta0 * exp( mt(1 : NLOCAL) )
+ %mt_vec(nmod_str+1 : nmod) = mt(nmod_str+1 : nmod)
+ mt_vec(1:nlocal) = beta0 * exp( mt(1:nlocal) );
+ mt_vec(nlocal+1 : nmod) = mt(nlocal+1 : nmod);
- if and(INV_SOURCE_T == 0, INV_SOURCE_X == 0)
- mt_vec(nlocal+1 : nmod) = m0_vec_initial(nlocal+1 : nmod);
- else
- mt_vec(nlocal+1 : nmod) = mt(nlocal+1 : nmod) + m0_vec_initial(nlocal+1 : nmod);
- end
-
else
% load chi for test model
chi_t_val = load(chitfile);
@@ -463,25 +461,17 @@
lam_0_val = xmin;
mk = m0 + lam_0_val * pk;
- % STRUCTURE PARAMETERS
- if (INV_STRUCT_BETA == 0)
- mk_vec(1:nlocal) = m0_vec_initial(1:nlocal); % initial structure
- else
- mk_vec(1:nlocal) = beta0 * exp( mk(1:nlocal) );
- end
+ m0_vec = zeros(nmod,1);
+ %m0_vec(1 : NLOCAL) = beta0 * exp( m0(1 : NLOCAL) )
+ %m0_vec(nmod_str+1 : nmod) = m0(nmod_str+1 : nmod)
+ mk_vec(1:nlocal) = beta0 * exp( mk(1:nlocal) );
+ mk_vec(nlocal+1 : nmod) = mk(nlocal+1 : nmod);
- % SOURCE PARAMETERS
- if and(INV_SOURCE_T == 0, INV_SOURCE_X == 0)
- mk_vec(nlocal+1 : nmod) = m0_vec_initial(nlocal+1 : nmod); % initial source
- else
- mk_vec(nlocal+1 : nmod) = mk(nlocal+1 : nmod) + m0_vec_initial(nlocal+1 : nmod);
- end
-
end
%---------------------------------------------------------
% icg = 1 ==> text
-elseif icg == 2
+elseif or(icg == 2, icg==3)
% load Matlab parameterization of the fields (gaussian_2D.m)
load([dirrand 'matlab_vars']);
@@ -526,42 +516,76 @@
%--------------------------------------
% covariance matrix
-
+
% modified covariance matrix
%Cmod = diag(1./Avec) * C * diag(1./Avec);
- % full covariance matrix
nmodi = ncell + nmod_src;
Cfull = zeros(nmodi,nmodi);
- Cfull(indBi, indBi) = C;
- Cfull(indTi, indTi) = diag(sigma_ts^2 * ones(length(indTi),1) );
- Cfull(indXi, indXi) = diag(sigma_xs^2 * ones(length(indXi),1) );
- Cfull(indYi, indYi) = diag(sigma_zs^2 * ones(length(indYi),1) );
- %figure; spy(Cfull);
-
+ Cmod = zeros(nmodi,nmodi);
cov_weight_vec = zeros(nmodi,1);
- cov_weight_vec(indBi) = (fac_str / ugrad_str) * fac_total * coverage_str ./ Avec.^2;
- cov_weight_vec(indTi) = (fac_ts / ugrad_ts) * fac_total * dnparm_src_run * coverage_src;
- cov_weight_vec(indXi) = (fac_xs / ugrad_xs) * fac_total * dnparm_src_run * coverage_src;
- cov_weight_vec(indYi) = (fac_ys / ugrad_ys) * fac_total * dnparm_src_run * coverage_src;
- Cmod = diag(sqrt(cov_weight_vec)) * Cfull * diag(sqrt(cov_weight_vec));
- Cmoddiag = diag(diag(Cmod));
+ if 1==1
+ % full, unweighted covariance matrix
+ Cfull(indBi, indBi) = C;
+ Cfull(indTi, indTi) = diag(sigma_ts^2 * ones(length(indTi),1) );
+ Cfull(indXi, indXi) = diag(sigma_xs^2 * ones(length(indXi),1) );
+ Cfull(indYi, indYi) = diag(sigma_zs^2 * ones(length(indYi),1) );
+ %figure; spy(Cfull);
+
+ cov_weight_vec = zeros(nmodi,1);
+ cov_weight_vec(indBi) = coverage_str / covg_weight_parts(1) ./ Avec.^2;
+ cov_weight_vec(indTi) = dnevent_run * coverage_src / covg_weight_parts(2);
+ cov_weight_vec(indXi) = dnevent_run * coverage_src / covg_weight_parts(3);
+ cov_weight_vec(indYi) = dnevent_run * coverage_src / covg_weight_parts(4);
+
+ % KEY command for full covariance matrix
+ Cmod = diag(sqrt(cov_weight_vec)) * Cfull * diag(sqrt(cov_weight_vec));
+
+ else
+ % diagonal, unweighted covariance matrix
+ Cfull(indBi, indBi) = diag(diag(C) ./ Avec.^2 * coverage_str );
+ Cfull(indTi, indTi) = diag(sigma_ts^2 * dnevent_run * coverage_src * ones(length(indTi),1) );
+ Cfull(indXi, indXi) = diag(sigma_xs^2 * dnevent_run * coverage_src * ones(length(indXi),1) );
+ Cfull(indYi, indYi) = diag(sigma_zs^2 * dnevent_run * coverage_src * ones(length(indYi),1) );
+
+ % test with a diagonal matrix
+ Cmod(indBi, indBi) = Cfull(indBi, indBi) / covg_weight_parts(1);
+ Cmod(indTi, indTi) = Cfull(indTi, indTi) / covg_weight_parts(2);
+ Cmod(indXi, indXi) = Cfull(indXi, indXi) / covg_weight_parts(3);
+ Cmod(indYi, indYi) = Cfull(indYi, indYi) / covg_weight_parts(4);
+ end
% equivalent representations
%cdiagmod = diag( diag(sqrt(cov_weight_vec)) * diag(diag(Cfull)) * diag(sqrt(cov_weight_vec)) );
%cmoddiag = diag(Cmod);
- % approximate C-inverse by using diagonal only
- cinv = 1./diag(Cmod);
+ % C-inverse for the model norm computations
+ Cmodinv = diag(1./diag(Cmod));
+% if icg==2
+% %cinv = 1./diag(Cfull);
+% %cmodinv = 1./diag(Cmod);
+% %Cinv = diag(1./diag(Cfull));
+% Cmodinv = diag(1./diag(Cmod));
+% else
+% %Cinv = inv(Cfull);
+% Cmodinv = inv(Cmod);
+% end
+ if icg==2
+ Crun = diag(diag(Cmod));
+ else
+ Crun = Cmod;
+ end
+
%--------------------------------------
% model
m0i = wave2d_m_gll2cell(m0,xg,yg,x,y,nx,ny,nmod_src,1);
- m0i_vec = wave2d_m_gll2cell(m0_vec,xg,yg,x,y,nx,ny,nmod_src,1);
- m0i_initial = wave2d_m_gll2cell(m0_initial,xg,yg,x,y,nx,ny,nmod_src,1);
- m0i_vec_initial = wave2d_m_gll2cell(m0_vec_initial,xg,yg,x,y,nx,ny,nmod_src,1);
+ m0i_vec = wave2d_m_gll2cell(m0_vec,xg,yg,x,y,nx,ny,nmod_src,1);
+ mpriori = wave2d_m_gll2cell(mprior,xg,yg,x,y,nx,ny,nmod_src,1);
+ %m0i_initial = wave2d_m_gll2cell(m0_initial,xg,yg,x,y,nx,ny,nmod_src,1);
+ %m0i_vec_initial = wave2d_m_gll2cell(m0_vec_initial,xg,yg,x,y,nx,ny,nmod_src,1);
%m0i = zeros(nmodi,1);
%m0i(1:ncell) = m0(iGLL);
%m0i(indTXYi) = m0(indTXY);
@@ -571,14 +595,21 @@
% interpolate kernels onto cells
kernels_alli = zeros(ncell,nevent_run);
- for ii=1:nevent_run
- kernels_alli(:,ii) = kernels_all(iGLL,ii);
+ for ii=1:nevent_run
+ Kbeta = kernels_all(iGLL,ii); % crude interpolation
+ kernels_alli(:,ii) = Kbeta;
if iplotker==1
- Zp = reshape( kernels_alli(:,ii), ny, nx);
+ % gradient
+ Zp = reshape( Kbeta, ny, nx);
%Zp = griddata(xg,yg,,Xp,Yp,'nearest');
figure; imagesc(Zp); set(gca,'ydir','normal');
caxis([-1 1]*max(abs(Kbeta))*0.8);
+
+ % steepest ascent vector
+ Zp = reshape( C * Kbeta, ny, nx);
+ figure; imagesc(Zp); set(gca,'ydir','normal');
+ caxis([-1 1]*max(abs(C * Kbeta))*0.8);
end
end
summed_keri = sum(kernels_alli,2);
@@ -588,26 +619,46 @@
gradienti_tot = zeros(nmodi,1);
gradienti_data(indBi) = summed_keri .* dAvec;
gradienti_data(indTXYi) = source_gradient;
- gradienti_model = cinv .* m0i;
+ gradienti_model = Cmodinv * (m0i - mpriori);
+ %gradienti_model = Cinv * (m0i - mpriori);
+ % zero out parts
+ if INV_SOURCE_T==0
+ gradienti_data(indTi) = 0;
+ gradienti_model(indTi) = 0;
+ end
+ if INV_SOURCE_X==0
+ gradienti_data(indXi) = 0;
+ gradienti_model(indXi) = 0;
+ gradienti_data(indYi) = 0;
+ gradienti_model(indYi) = 0;
+ end
+ if INV_STRUCT_BETA==0
+ gradienti_data(indBi) = 0;
+ gradienti_model(indBi) = 0;
+ end
gradienti_tot = gradienti_data + gradienti_model;
% gradient -- TO CHECK
% NOTE: due to the new discretization of the kernels, we expect the
% norm of the structure gradient to be slightly different from before
- wave2d_gnorm_sq(gradienti_tot,diag(Cmod),m_indsi);
- wave2d_gnorm_sq(gradienti_data,diag(Cmod),m_indsi);
- wave2d_gnorm_sq(gradienti_model,diag(Cmod),m_indsi);
-
+ %wave2d_gnorm_sq(gradienti_tot,Crun,m_indsi);
+ %wave2d_gnorm_sq(gradienti_data,Crun,m_indsi);
+ %wave2d_gnorm_sq(gradienti_model,Crun,m_indsi);
+ wave2d_norm_sq(gradienti_tot,mpriori,Crun,m_indsi,ones(4,1),0);
+ wave2d_norm_sq(gradienti_data,mpriori,Crun,m_indsi,ones(4,1),0);
+ wave2d_norm_sq(gradienti_model,mpriori,Crun,m_indsi,ones(4,1),0);
+
% check parts of model norm
- model_norm_partsi = wave2d_gnorm_sq(m0i,1./diag(Cmod),m_indsi);
+ model_norm_partsi = wave2d_norm_sq(m0i,mpriori,Crun,m_indsi,covm_weight_parts,1);
model_normi = sum(model_norm_partsi);
% check
- model_norm0, model_normi
+ model_norm0, model_normi % sum( m0i(1:ncell).^2 .* cinv(1:ncell) )
% check chi model values
chi0 = load([idir1 'chi.dat'])
chi_k_val = 0.5*(data_norm + model_normi)
-
+ chi_k_val = chi0 % just use the real chi value
+
%---------------------------------------------------------
% emulate CG algorithm in wave2d.f90 -- THIS ASSUMES A FULL COVARIANCE MATRIX
@@ -615,45 +666,23 @@
gfile = [idir2 gfilemp];
mu_val = chi_data_stop;
gk = gradienti_tot;
- [M,pk] = wave2d_cg(m0i,gk,Cmoddiag,chi_k_val,mu_val,istep,imaketest,gfile,chitfile);
+ [M,pk] = wave2d_cg(m0i,gk,Crun,chi_k_val,mu_val,istep,imaketest,gfile,chitfile);
% convert model vector entries to physical values -- also ensure that
% entries are unchanged if you are not inverting for them
if imaketest==1
mti = M;
mti_vec = zeros(nmodi,1);
- % STRUCTURE PARAMETERS
- if (INV_STRUCT_BETA == 0)
- mti_vec(1:ncell) = m0_vec_initial(1:ncell);
- else
- mti_vec(1:ncell) = beta0 * exp( mti(1:ncell) );
- end
- % SOURCE PARAMETERS
- if and(INV_SOURCE_T == 0, INV_SOURCE_X == 0)
- mti_vec(ncell+1 : nmodi) = m0i_vec_initial(ncell+1 : nmodi);
- else
- mti_vec(ncell+1 : nmodi) = mti(ncell+1 : nmodi) + m0i_vec_initial(ncell+1 : nmodi);
- end
+ mti_vec(1:ncell) = beta0 * exp( mti(1:ncell) );
+ mti_vec(ncell+1 : nmodi) = mti(ncell+1 : nmodi);
else
mki = M;
mki_vec = zeros(nmodi,1);
- % STRUCTURE PARAMETERS
- if (INV_STRUCT_BETA == 0)
- mki_vec(1:ncell) = m0_vec_initial(1:ncell); % initial structure
- else
- mki_vec(1:ncell) = beta0 * exp( mki(1:ncell) );
- end
- % SOURCE PARAMETERS
- if and(INV_SOURCE_T == 0, INV_SOURCE_X == 0)
- mki_vec(ncell+1 : nmodi) = m0i_vec_initial(ncell+1 : nmodi); % initial source
- else
- mki_vec(ncell+1 : nmodi) = mki(ncell+1 : nmodi) + m0i_vec_initial(ncell+1 : nmodi);
- end
-
+ mki_vec(1:ncell) = beta0 * exp( mki(1:ncell) );
+ mki_vec(ncell+1 : nmodi) = mki(ncell+1 : nmodi);
end
-
else
error('invalid icg value');
@@ -700,8 +729,8 @@
% write structure model -- ONLY MU CHANGES (only m_B appears)
m_str_beta_new = beta0 * exp( m_B );
- m_str_rho_new = rho0*ones(nlocal,1);
- m_str_kappa_new = kappa0*ones(nlocal,1);
+ m_str_rho_new = rho0 * ones(nlocal,1);
+ m_str_kappa_new = kappa0 * ones(nlocal,1);
m_str_mu_new = m_str_rho_new .* m_str_beta_new.^2;
wave2d_write_str([odir 'structure_syn_m' stimo '.dat'],xg,yg,...
m_str_kappa_new,m_str_mu_new,m_str_rho_new,m_B);
@@ -713,10 +742,10 @@
% read sources for data
[m_src_lon,m_src_lat,m_src_ts_dat,m_src_xs_dat,m_src_ys_dat,m_src_ts_d,m_src_xs_d,m_src_ys_d] ...
= textread([idir1 'src_dat.dat'],'%f%f%f%f%f%f%f%f');
-
- m_src_ts_d_new = m_vec_ts - m_src_ts_dat;
- m_src_xs_d_new = m_vec_xs - m_src_xs_dat;
- m_src_ys_d_new = m_vec_ys - m_src_ys_dat;
+ % mdat - msrc : how far to still go
+ m_src_ts_d_new = m_src_ts_dat - m_vec_ts;
+ m_src_xs_d_new = m_src_xs_dat - m_vec_xs;
+ m_src_ys_d_new = m_src_ys_dat - m_vec_ys;
wave2d_write_src([odir 'src_syn_m' stimo '.dat'], m_src_lon, m_src_lat,...
m_vec_ts, m_vec_xs, m_vec_ys, m_src_ts_d_new, m_src_xs_d_new, m_src_ys_d_new);
@@ -734,8 +763,6 @@
end
-break
-
%---------------------------------------------------------
% compute CG update using a diagonal covariance matrix with interpolated cells for model parameters
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-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-02-12 20:17:28 UTC (rev 16258)
@@ -140,7 +140,7 @@
! weighting and smoothing
double precision :: gamma ! sigma
- double precision, dimension(NGLLX,NGLLZ,NSPEC) :: kbulk_smooth, kbeta_smooth ! scale_array
+ double precision, dimension(NGLLX,NGLLZ,NSPEC) :: kbeta_smooth ! kbulk_smooth
! double precision, dimension(NGLLX,NGLLZ,NSPEC) :: bulk_pert_syn, beta_pert_syn
double precision :: scale_bulk, scale_beta, scale_source, scale_struct, scale_struct_gji
double precision :: sigma_bulk0, sigma_beta0, sigma_bulk, sigma_beta, sigma_checker_scale
@@ -148,13 +148,14 @@
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_m, fac_total_g
- double precision :: ugrad_str, ugrad_ts, ugrad_xs, ugrad_ys, ugrad_src
+ double precision :: fac_str, fac_ts, fac_xs, fac_ys, fac_total_m ! fac_total_g
+ double precision :: ugsq_str, ugsq_ts, ugsq_xs, ugsq_ys, ugsq_src
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, 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 :: beta_val, beta_val_top, beta_val_bot, lam_0_val, lam_t_val, lam_t_val_bot
+ double precision :: chi_t_val, chi_k_val, mu_val, gnorm
double precision :: Pa,Pb,Pc,Pd,qfac,xx1,xx2,yy1,yy2,g1,g2,xmin
! arrays describing the structure and source parameters
@@ -206,19 +207,19 @@
!--------------------------------------
! stop program if there are certain unallowed paramter combinations
- if(READ_IN) then
+ if (READ_IN) then
print *, 'For READ_IN runs...'
- if(NITERATION /= 0) stop 'NITERATION = 0'
- if(ISRC_SPACE /= 6) stop 'ISRC_SPACE = 6'
- hmod = 1 ! iteration number for read-in models
+ if (NITERATION /= 0) stop 'NITERATION = 0'
+ if (ISRC_SPACE /= 6) stop 'ISRC_SPACE = 6'
+ hmod = 6 ! iteration number for read-in models
endif
- if( IKER/=0 .and. IKER/=1 .and. IKER/=2 ) stop 'IKER must = 0,1,2'
+ if ( IKER/=0 .and. IKER/=1 .and. IKER/=2 ) stop 'IKER must = 0,1,2'
- if( ADD_DATA_ERRORS .and. IKER/=1 ) stop 'for now, IKER = 1 for adding errors to data'
+ if ( ADD_DATA_ERRORS .and. IKER/=1 ) stop 'for now, IKER = 1 for adding errors to data'
! if you do not want kernels, then you cannot iterate
- if( (.not. COMPUTE_KERNELS) .and. NITERATION /= 0 ) stop 'NITERATION = 0'
+ if ( (.not. COMPUTE_KERNELS) .and. NITERATION /= 0 ) stop 'NITERATION = 0'
!--------------------------------------
@@ -235,8 +236,8 @@
!socal_dir = '/net/denali/scratch1/carltape/socal_modes/socal_modes/local_modes/'
! whether a point source or finite source is used
- if(ISRC_SPACE==1) FINITE_SOURCE = 0
- if(ISRC_SPACE/=1) FINITE_SOURCE = 1
+ 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
@@ -252,7 +253,7 @@
! load or specify events, depending on whether you are computing
! membrane surface waves (ISURFACE==1) or body waves
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
!----------
! read in target events (socal_quakes.m)
@@ -262,7 +263,7 @@
quake_file = 'socal_quakes_N025.dat'
open(55,file=trim(in_dir)//trim(quake_file),status='unknown')
read(55,*) nevent
- if(nevent > MAX_EVENT) stop 'nevent > MAX_EVENT (so increase MAX_EVENT)'
+ if (nevent > MAX_EVENT) stop 'nevent > MAX_EVENT (so increase MAX_EVENT)'
print *, nevent, ' target events (for synthetics)'
read(55,'(i14,3f14.7)') (socal_events_lab(i), &
socal_events_lon(i), socal_events_lat(i), socal_events_mag(i), i = 1,nevent)
@@ -280,7 +281,7 @@
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
+ if (0==1) then ! NEW perturbations
! perturbations generated from wave2d_sigmas.m
open(19,file='INPUT/events_txy_pert_initial.dat',status='old')
@@ -386,13 +387,13 @@
! all vectors should be this length
nevent = MAX_EVENT
- if(0==1) then ! read in target events (socal1D.m)
+ if (0==1) then ! read in target events (socal1D.m)
quake_file = 'socal_quakes_1D_rand10.dat' ! 200 km width
!quake_file = 'socal_quakes_1D_rand15.dat' ! 400 km width
open(55,file=trim(in_dir)//trim(quake_file),status='unknown')
read(55,*) nevent
- if(nevent > MAX_EVENT) stop 'nevent > MAX_EVENT (so increase MAX_EVENT)'
+ if (nevent > MAX_EVENT) stop 'nevent > MAX_EVENT (so increase MAX_EVENT)'
print *, nevent, ' target events (for synthetics)'
read(55,'(i10,3f24.12)') (junk1,x_eve0_dat(i),z_eve0_dat(i),junk2,i = 1,nevent)
close(55)
@@ -463,7 +464,7 @@
! if you are not computing kernels, loop over 25 TSVD models to read in
qmax = 25
- if(COMPUTE_KERNELS .or. (READ_IN .and. READ_SINGLE) ) qmax = 1
+ if (COMPUTE_KERNELS .or. (READ_IN .and. READ_SINGLE) ) qmax = 1
!============================================
! LOOP 1: different tomographic runs
@@ -473,17 +474,17 @@
! KEY COMMAND: scalelength of checker for velocity models (1,2,3)
Nfac = 3 ! use Nfac=3 for one-source examples
- if(READ_IN) then
+ if (READ_IN) then
irun0 = IRUNZ
else
irun0 = IRUNZ + 20*(iq-1) ! increment is 20
endif
! name the reference output directory for the optimization run
- if(READ_IN) then
+ if (READ_IN) then
! inversion steps done in Matlab
- if(READ_SINGLE) then
+ if (READ_SINGLE) then
write(out_dir2,'(a,i4.4,a,i4.4,a)') trim(out_dir3)//'run_',irun0,'/READ_IN/model_m',hmod,'/'
else
write(out_dir2,'(a,i4.4,a,i4.4,a,i4.4,a)') trim(out_dir3)//'run_',irun0,'/READ_IN/model_m',hmod,'/run_p',iq,'/'
@@ -536,7 +537,7 @@
sda_global_mean = sqrt(da_global_mean)
sda_local_mean = sqrt(da_local_mean)
- if(1==1) then
+ if (1==1) then
! global mesh
open(unit=15,file=trim(out_dir2)//'global_mesh.dat',status='unknown')
do iglob = 1,NGLOB
@@ -560,7 +561,7 @@
close(15)
endif
- if(0==1) then
+ if (0==1) then
! corner points for each element, and centerpoint (in km)
open(unit=15,file='elements.dat',status='unknown')
do ispec = 1,nspec
@@ -626,7 +627,7 @@
!--------------------------------------
! determine the UTM coordinates of your origin and corners
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
utm_xmin = 0.0 ; utm_zmin = 0.0
xtemp = LON_MIN ; ztemp = LAT_MIN
@@ -687,11 +688,11 @@
!--------------------------------------
! assign values for plotting grid
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
x_plot(:) = x_lon(:)
z_plot(:) = z_lat(:)
- elseif(ISURFACE==0) then
+ elseif (ISURFACE==0) then
x_plot(:) = x(:)
z_plot(:) = z(:)
@@ -704,11 +705,11 @@
!!$ ! homogeneous or heterogeneous phase velocity map
!!$
-!!$ if(IMODEL_DAT == 1 .or. IMODEL_SYN == 1) then ! read or compute a heterogeneous map
+!!$ if (IMODEL_DAT == 1 .or. IMODEL_SYN == 1) then ! read or compute a heterogeneous map
!!$ ! This calls model_files/get_socal_map.csh, which executes wave2d_socal.f90
!!$ ! from /net/denali/scratch1/carltape/socal_modes/local_modes
!!$
-!!$ if(0==1) then ! compute phase velocity map
+!!$ if (0==1) then ! compute phase velocity map
!!$
!!$ ! write lat-lon gridpoints to file
!!$ filename = trim(socal_dir) // 'socal_input.dat'
@@ -760,17 +761,17 @@
!!$
!!$ endif
!!$
-!!$ if(IMODEL_SYN == 1) c_glob_syn(:) = c_glob(:)
-!!$ if(IMODEL_DAT == 1) c_glob_dat(:) = c_glob(:)
+!!$ if (IMODEL_SYN == 1) c_glob_syn(:) = c_glob(:)
+!!$ if (IMODEL_DAT == 1) c_glob_dat(:) = c_glob(:)
! reference phase velocity (m/s) (should be obtained based on hdur)
! (This is only relevant for homogeneous models.)
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
alpha0 = 5800.0
beta0 = 3500.0 ! (3500 for GJI)
rho0 = DENSITY
- elseif(ISURFACE==0) then
+ elseif (ISURFACE==0) then
!alpha0 = a_layers(2) ! layer 2 of SoCal model (1D)
!beta0 = b_layers(2) ! layer 2 of SoCal model (1D)
alpha0 = (a_layers(2)+a_layers(1))/2
@@ -778,7 +779,7 @@
endif
! parameters for checkerboard pattern
- if(IMODEL_DAT==2 .or. IMODEL_SYN==2) then
+ if (IMODEL_DAT==2 .or. IMODEL_SYN==2) then
!Nfac = 3 ! scalelength of map = N * (wavelength of surface waves for hdur)
afac = 10.0 ! max PERCENT variation of synthetic model (10.0)
@@ -792,7 +793,8 @@
! NOTE: THESE ARE ONLY FOR A SPECIFIC MESH
if ((IMODEL_SYN == 3).and.(IMODEL_DAT == 3)) then
- idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_001/'
+ !idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_001/'
+ idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_002/'
! read in structure model for synthetics
write(filename2,'(a)') trim(idir)//'structure_syn_in.dat'
@@ -824,6 +826,14 @@
enddo
close(20)
+ ! read in sigma value describing the distribution
+ write(filename,'(a)') trim(idir)//'structure_sigma.dat'
+ open(unit=21,file=filename,status='old')
+ read(21,*) sigma_beta0
+ close(21)
+ sigma_bulk0 = sigma_beta0
+
+ ! compute beta for initial and target models
beta_syn = sqrt( mu_syn / rho_syn )
beta_dat = sqrt( mu_dat / rho_dat )
@@ -833,10 +843,6 @@
read(16,'(6e20.10)') alpha0, beta0, rho0, bulk0, kappa0, mu0
close(16)
- ! sigma values describing the distribution
- sigma_beta0 = 0.05
- sigma_bulk0 = 0.05
-
endif
!=========================================
@@ -853,7 +859,7 @@
! smooth the reference model
! (This was motivated for the imaging principle kernel tests, Feb 2007.)
- if(ISMOOTH_INITIAL_MODEL == 1) then
+ if (ISMOOTH_INITIAL_MODEL == 1) then
call smooth_function(kappa_syn, GAMMA_SMOOTH_MODEL, temp_local1) ; kappa_syn = temp_local1;
call smooth_function( mu_syn, GAMMA_SMOOTH_MODEL, temp_local1) ; mu_syn = temp_local1;
call smooth_function( rho_syn, GAMMA_SMOOTH_MODEL, temp_local1) ; rho_syn = temp_local1;
@@ -896,7 +902,7 @@
! if READ_IN option, then READ in the structure file (local level)
! NOTE: Assume that the mesh and scaling values are IDENTICAL
! to what was used in the base directory for the CG algorithm.
- if (READ_IN .and. INV_STRUCT_BETA == 1) then
+ if ( (READ_IN) .and. (INV_STRUCT_BETA==1) ) then
kappa_syn = 0.0 ; mu_syn = 0.0 ; rho_syn = 0.0
alpha_syn = 0.0 ; beta_syn = 0.0 ; bulk_syn = 0.0
@@ -928,7 +934,7 @@
close(19)
! smooth the input model
- if(ISMOOTH_MODEL_UPDATE == 1) then
+ if (ISMOOTH_MODEL_UPDATE == 1) then
call smooth_function(kappa_syn, GAMMA_SMOOTH_MODEL, temp_local1) ; kappa_syn = temp_local1;
call smooth_function( mu_syn, GAMMA_SMOOTH_MODEL, temp_local1) ; mu_syn = temp_local1;
call smooth_function( rho_syn, GAMMA_SMOOTH_MODEL, temp_local1) ; rho_syn = temp_local1;
@@ -940,7 +946,7 @@
endif
! unperturbed structure (set target map to be the map used to compute synthetics)
- if(PERT_STRUCT_BETA == 0) then
+ if (PERT_STRUCT_BETA == 0) then
mu_dat = mu_syn
beta_dat = beta_syn
kappa_dat = kappa_syn
@@ -1042,10 +1048,10 @@
! estimate the time step
dh = HEIGHT/dble((NGLLZ-1)*NEZ)
- if(dh > LENGTH/dble((NGLLX-1)*NEX)) dh = LENGTH/dble((NGLLX-1)*NEX)
+ if (dh > LENGTH/dble((NGLLX-1)*NEX)) dh = LENGTH/dble((NGLLX-1)*NEX)
print *
print *,' space step (km) :', sngl(dh/1000.0)
- if(ISURFACE==0) then
+ if (ISURFACE==0) then
print *,'time step est from courant = 0.2, Pmax : ',sngl(0.2*dh/alpha_max),' seconds'
print *,'time step est from courant = 0.2, Pmin : ',sngl(0.2*dh/alpha_min),' seconds'
print *,'time step est from courant = 0.2, Smax : ',sngl(0.2*dh/beta_max),' seconds'
@@ -1059,7 +1065,7 @@
!--------------------------------------
! events for synthetics (eglob_syn)
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! fill a portion of the [1:MAX_EVENT] vector with the SoCal events
x_eve_lon0_syn(1:nevent) = socal_events_lon(1:nevent)
z_eve_lat0_syn(1:nevent) = socal_events_lat(1:nevent)
@@ -1078,10 +1084,10 @@
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'
+ if (nevent < 1) stop 'Must have at least one event'
! allocate variables
- if(ISURFACE==1) allocate(x_eve_lon_syn(nevent),z_eve_lat_syn(nevent))
+ if (ISURFACE==1) allocate(x_eve_lon_syn(nevent),z_eve_lat_syn(nevent))
allocate(x_eve_syn(nevent),z_eve_syn(nevent))
allocate(eglob_syn(nevent))
allocate(ispec_eve_syn(nevent),xi_eve_syn(nevent),gamma_eve_syn(nevent))
@@ -1121,7 +1127,7 @@
enddo
! display target events and final events -- and the distance between (in meters)
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! 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)
@@ -1138,7 +1144,7 @@
write(*,'(i8, 5f17.10)') i, temp1, temp2, temp3, temp4, temp5*6371.0*1000.0
enddo
- elseif(ISURFACE==0) then
+ elseif (ISURFACE==0) then
print *
print *, 'events, km [x_eve0_syn, z_eve0_syn, x_eve_syn, x_eve_syn, dist]:'
do i = 1,nevent
@@ -1152,13 +1158,13 @@
endif
! write events for PRIOR MODEL to file
- if( .not. READ_IN ) then
+ 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
+ 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
@@ -1168,13 +1174,13 @@
endif
! write events for SYNTHETICS to file
- if( .not. READ_IN ) then
+ if ( .not. READ_IN ) then
open(19,file=trim(out_dir2)//'events_syn_xz.dat',status='unknown')
do ievent = 1,nevent
write(19,'(3f20.10,1i12)') x_eve_syn(ievent), z_eve_syn(ievent), otime_syn(ievent), ievent
enddo
close(19)
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
open(18,file=trim(out_dir2)//'events_syn_lonlat.dat',status='unknown')
do ievent = 1,nevent
!write(18,*) sngl(x_lon(eglob_syn(ievent))), sngl(z_lat(eglob_syn(ievent))), ievent
@@ -1192,15 +1198,15 @@
! events for data (eglob_dat)
! allocate variables
- if(ISURFACE==1) allocate(x_eve_lon_dat(nevent),z_eve_lat_dat(nevent))
+ if (ISURFACE==1) allocate(x_eve_lon_dat(nevent),z_eve_lat_dat(nevent))
allocate(x_eve_dat(nevent),z_eve_dat(nevent))
allocate(eglob_dat(nevent))
allocate(ispec_eve_dat(nevent),xi_eve_dat(nevent),gamma_eve_dat(nevent))
allocate(otime_dat(nevent))
- if( PERT_SOURCE_T == 1 .or. PERT_SOURCE_X == 1 ) then ! source perturbations
+ 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)
+ 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(:)
@@ -1226,8 +1232,8 @@
!!$ 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
+!!$ 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
@@ -1243,7 +1249,7 @@
!!$ close(19)
!!$ otime_dat(5) = otime_syn(5) + 0.8 ! testing for event 5 (0.8)
- elseif(ISURFACE==0) then
+ elseif (ISURFACE==0) then
! perturb the events to get target events for the synthetics
! perturbation is described as (r,theta) polar coordinates
@@ -1294,7 +1300,7 @@
eglob_dat(i), ispec_eve_dat(i), xi_eve_dat(i), gamma_eve_dat(i)
enddo
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! convert from mesh to lat-lon
call mesh_geo(nevent, x_eve_lon_dat, z_eve_lat_dat, x_eve_dat, z_eve_dat, UTM_PROJECTION_ZONE, IMESH2LONLAT)
@@ -1313,7 +1319,7 @@
else ! no source perturbation for target sources
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
x_eve_lon_dat(:) = x_eve_lon_syn(:)
z_eve_lat_dat(:) = z_eve_lat_syn(:)
endif
@@ -1328,13 +1334,13 @@
endif ! PERT_SOURCE
! write events for data to file
- if( .not. READ_IN ) then
+ if ( .not. READ_IN ) then
open(19,file=trim(out_dir2)//'events_dat_xz.dat',status='unknown')
do ievent = 1,nevent
write(19,'(3f20.10,1i12)') x_eve_dat(ievent), z_eve_dat(ievent), otime_dat(ievent), ievent
enddo
close(19)
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
open(19,file=trim(out_dir2)//'events_dat_lonlat.dat',status='unknown')
do ievent = 1,nevent
write(19,'(2f20.10,1i12)') x_eve_lon_dat(ievent), z_eve_lat_dat(ievent), ievent
@@ -1357,9 +1363,9 @@
! get the lat-lon of the TARGET RECEIVERS
- if(ISURFACE==1) then ! enter target lat-lon
+ if (ISURFACE==1) then ! enter target lat-lon
- if(IREC_SPACE==1) then ! individual receivers
+ if (IREC_SPACE==1) then ! individual receivers
! target receiver
!x_rec0(1) = 3 * LENGTH/4 ; z_rec0(1) = HEIGHT/2
@@ -1378,7 +1384,7 @@
nrec = 3
- elseif(IREC_SPACE==2) then ! actual station locations
+ elseif (IREC_SPACE==2) then ! actual station locations
! read in (target) receivers
recfile = trim(in_dir)//'STATION_149_full'
@@ -1390,7 +1396,7 @@
read(88,*) (x_rec_lon0(i),z_rec_lat0(i),i = 1,nrec)
close(88)
- elseif(IREC_SPACE==4) then ! 'regular' mesh of receivers
+ elseif (IREC_SPACE==4) then ! 'regular' mesh of receivers
! calculate mesh spacing
dx = LENGTH/NMESH_REC
@@ -1413,12 +1419,12 @@
endif ! IREC_SPACE
- elseif(ISURFACE==0) then ! enter target x-z
+ elseif (ISURFACE==0) then ! enter target x-z
- if(IREC_SPACE==1) then ! individual receivers
+ if (IREC_SPACE==1) then ! individual receivers
stop ' no individual receiver listed'
- elseif(IREC_SPACE==2) then ! actual station locations
+ elseif (IREC_SPACE==2) then ! actual station locations
! read in (target) receivers
recfile = trim(in_dir)//'STATIONS_socal1D_rand15' ! 200 km width
@@ -1428,7 +1434,7 @@
read(88,*) (x_rec0(i),z_rec0(i),i = 1,nrec)
close(88)
- elseif(IREC_SPACE==4) then ! 'regular' line of receivers (BODY waves)
+ elseif (IREC_SPACE==4) then ! 'regular' line of receivers (BODY waves)
! calculate mesh spacing
dx = LENGTH/NMESH_REC
@@ -1454,7 +1460,7 @@
endif
! make sure that there are fewer target points than the max allowed
- if(nrec > MAX_SR) then
+ if (nrec > MAX_SR) then
print *
print *, ' IREC_SPACE = ', IREC_SPACE
print *, ' nrec = ', nrec
@@ -1468,7 +1474,7 @@
! all vectors should be length MAX_SR
nrec = MAX_SR
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! convert from lat-lon to mesh coordinates (in meters)
call mesh_geo(nrec,x_rec_lon0,z_rec_lat0,x_rec0,z_rec0,UTM_PROJECTION_ZONE,ILONLAT2MESH)
endif
@@ -1478,10 +1484,10 @@
!call station_filter(nrec, x_rec, z_rec, dmin_trsh_r, ncoast, coast_x, coast_z)
!call station_filter_2(nrec, x_rec, z_rec, -1) ! -1 for left, +1 for right
- if(nrec < 1) stop 'Must have at least one receiver'
+ if (nrec < 1) stop 'Must have at least one receiver'
! allocate vectors
- if(ISURFACE==1) allocate(x_rec_lon(nrec),z_rec_lat(nrec))
+ if (ISURFACE==1) allocate(x_rec_lon(nrec),z_rec_lat(nrec))
allocate(x_rec(nrec),z_rec(nrec))
allocate(rglob(nrec))
allocate(ispec_rec(nrec),xi_rec(nrec),gamma_rec(nrec))
@@ -1502,7 +1508,7 @@
write(*,'(i8,2e18.8,2i8,2e18.8)') i, x_rec(i), z_rec(i), rglob(i), ispec_rec(i), xi_rec(i), gamma_rec(i)
enddo
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! convert from mesh to lat-lon
call mesh_geo(nrec, x_rec_lon, z_rec_lat, x_rec, z_rec, UTM_PROJECTION_ZONE, IMESH2LONLAT)
@@ -1534,7 +1540,7 @@
!deallocate(x_rec,z_rec)
! write receivers to file
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
open(18,file=trim(out_dir2)//'recs_lonlat.dat',status='unknown')
do irec = 1,nrec
write(18,'(2f20.10,1i12)') x_rec_lon(irec), z_rec_lat(irec), irec
@@ -1601,7 +1607,7 @@
! THIS PORTION COULD BE REMOVED -- IT WAS FOR COMPUTING THE SPECTRAL MAP
! FOR THE OCEAN MICROSEISM SIMULATIONS.
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! mesh of target points
x_recf0(:) = 0.0
@@ -1612,7 +1618,7 @@
do xmesh = 0.0,LENGTH,dx
do zmesh = 0.0,HEIGHT,dz
i = i+1
- if(i > MAX_SR_FAKE) stop 'i > MAX_SR_FAKE so change dx, dz, or MAX_SR_FAKE'
+ if (i > MAX_SR_FAKE) stop 'i > MAX_SR_FAKE so change dx, dz, or MAX_SR_FAKE'
x_recf0(i) = xmesh
z_recf0(i) = zmesh
enddo
@@ -1628,7 +1634,7 @@
! filter target points (utm-mesh) -- update nrecf
call station_filter(nrecf, x_recf0, z_recf0, ifilter_recf, STATION_GRID_BUFFER)
- if(nrecf < 1) stop 'Must have at least one fake (adjoint) receiver'
+ if (nrecf < 1) stop 'Must have at least one fake (adjoint) receiver'
! allocate vectors
allocate(x_recf(nrecf),z_recf(nrecf),fglob(nrecf))
@@ -1647,7 +1653,7 @@
! print *, sngl(x(fglob(i))/1000.0),sngl(z(fglob(i))/1000.0),sngl(x_recf(i)/1000.0),sngl(z_recf(i)/1000.0)
!enddo
- if(0==1) then
+ if (0==1) then
! display target receivers and final fake receivers -- distances in km
print *
print *, ' fake receivers [x_rec0, z_rec0, x_rec, x_rec, dist (km)]:'
@@ -1822,23 +1828,23 @@
temp_local1(:,:,:) = log( beta_dat(:,:,:) / beta0 )
call local2mvec(temp_local1, nmod_src, m_src_dat, nmod, mtarget)
- ! write out initial source vector
- !if( .not. READ_IN ) then
- open(89,file=trim(out_dir2)//'initial_structure.dat',status='unknown')
+ ! write out prior, initial, and target models
+ if ( .not. READ_IN ) then
+ open(89,file=trim(out_dir2)//'prior_initial_target_models.dat',status='unknown')
do i = 1,nmod
write(89,'(3e24.12)') mprior(i), m0(i), mtarget(i)
enddo
close(89)
- 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,'(3e24.12)') m_src_prior(i), m_src_syn(i), m_src_dat(i)
- enddo
- close(88)
- !endif
+!!$ 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,'(3e24.12)') m_src_prior(i), m_src_syn(i), m_src_dat(i)
+!!$ enddo
+!!$ close(88)
+ endif
!------------------
! MODEL COVARIANCE MATRIX = inverse metric tensor
@@ -1860,7 +1866,7 @@
endif
! scaling vector for SOURCE parameters
- if( GJI_PAPER == 1 ) then
+ if ( GJI_PAPER == 1 ) then
!joint_str_src_grad_ratio = (70000.0 / 5000.0) ! to match the new version
!joint_str_src_grad_ratio = 1.0
m_scale_src(1) = 2*hdur ! dts is scaled by period (20s)
@@ -1884,8 +1890,8 @@
!!$ ! A simple way is to incorporate a scaling factor in the sigma-structure terms.
!!$ ! This will control the relative weight of each set of parameters.
!!$ ! NOTE : We approximate our checkerboard-generated structure values by a Gaussian distribution.
-!!$ if( INV_SOURCE == 1 .and. INV_STRUCT_BETA == 1) then
-!!$ if( READ_IN ) then
+!!$ if ( INV_SOURCE == 1 .and. INV_STRUCT_BETA == 1) then
+!!$ if ( READ_IN ) then
!!$ joint_str_src_grad_ratio = 1.0
!!$ else
!!$ ! see also scale_struct_gji (F^2)
@@ -1924,6 +1930,7 @@
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
+ !fac_str = 0.85 ; fac_ts = 0.05 ; fac_xs = 0.05 ; fac_ys = 0.05
else
stop 'you must invert for something'
endif
@@ -1936,45 +1943,79 @@
covm_weight_parts(:) = covm_weight_parts(:) / sum(covm_weight_parts)
!covm_weight_parts(:) = 1.0 ! OLD VERSION
+ ! simple way to balance GRADIENT according to norm (NOT norm-squared)
+ covm_weight_parts = covm_weight_parts * covm_weight_parts
+
! unbalanced initial gradient values (gradient_norm_all_unbalanced.dat)
+ ! ugsq --> unbalanced gradient norm-squared
covg_weight_parts(:) = 1.0
- fac_total_g = 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)
+
+ ! NOTE: we obtain these values from either:
+ ! gradient_norm_data_all_unbalanced.dat
+ ! gradient_norm_all_unbalanced.dat
+
+ ugsq_str = 1.0 ; ugsq_ts = 1.0 ; ugsq_xs = 1.0 ; ugsq_ys = 1.0
+
! 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
+ !ugsq_str = 0.1296496361d6 ; ugsq_ts = 0.2706999550d4; ugsq_xs = 0.1334372694d4 ; ugsq_ys = 0.1743549898d4 ! Nfac = 2, 25 events
- !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
+ !ugsq_str = 0.2070999127d6 ; ugsq_ts = 0.3507740673d4 ; ugsq_xs = 0.1940368586d4 ; ugsq_ys = 0.2142118702d4 ! Nfac = 3, 25 events
+ !ugsq_str = 0.4247330856d6 ; ugsq_ts = 0.5896527525d4 ; ugsq_xs = 0.2643538184d4 ; ugsq_ys = 0.3342815391d4 ! Nfac = 3, 5 events
+ !ugsq_str = 0.1205146534d6 ; ugsq_ts = 0.2926409454d3 ; ugsq_xs = 0.9695606936d3 ; ugsq_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
- ! 3200/5000
- !ugrad_str = 0.3026988450d6 ; ugrad_ts = 0.8260031407d3 ; ugrad_xs = 0.1187469038d4 ; ugrad_ys = 0.1381078706d4 ! Gaussians, 25 events
+ ! 9200
+ ugsq_str = 0.8418200277d6 ; ugsq_ts = 0.4639935107d3 ; ugsq_xs = 0.1635906965d4 ; ugsq_ys = 0.1147402541d4 ! Gaussians, 1 event
+
+ ! 1200 - with new source
+ !ugsq_str = 0.1176265806d6 ; ugsq_ts = 0.1823448342d3 ; ugsq_xs = 0.6421163226d2 ; ugsq_ys = 0.2776564145d1 ! Gaussians, 1 event
+ ! 1200 - balance full gradient
+ !ugsq_str = 0.1176904168d6 ; ugsq_ts = 0.1443574694d3 ; ugsq_xs = 0.8206282043d2 ; ugsq_ys = 0.1632463053d2 ! Gaussians, 1 event
+
+ ! 1300 - with 0.7,0.1,0.1,0.1, old source, MOISMPRIOR
+ !ugsq_str = 0.1191920326d6 ; ugsq_ts = 0.4533353832d3 ; ugsq_xs = 0.6191223030d2 ; ugsq_ys = 0.4387672586d2 ! Gaussians, 1 event
+
+ ! 3200,3600
+ !ugsq_str = 0.3026988450d6 ; ugsq_ts = 0.8260031407d3 ; ugsq_xs = 0.1187469038d4 ; ugsq_ys = 0.1381078706d4 ! Gaussians, 25 events
+ ! 3400,3800
+ !ugsq_str = 0.3028636452d6 ; ugsq_ts = 0.8327738622d3 ; ugsq_xs =0.1195006162d4 ; ugsq_ys = 0.1381645351d4 ! Gaussians, 25 events
! 3300/1700
- !ugrad_str = 0.3352965529d6 ; ugrad_ts = 0.8068082292d3 ; ugrad_xs = 0.1209790150d4 ; ugrad_ys = 0.1391331978d4 ! Gaussians, 25 events
+ !ugsq_str = 0.3352965529d6 ; ugsq_ts = 0.8068082292d3 ; ugsq_xs = 0.1209790150d4 ; ugsq_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
+!!$ ! 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
+!!$
+!!$ covg_weight_parts(1) = ugsq_str/fac_str
+!!$ covg_weight_parts(2) = ugsq_ts/fac_ts
+!!$ covg_weight_parts(3) = ugsq_xs/fac_xs
+!!$ covg_weight_parts(4) = ugsq_ys/fac_ys
+!!$ covg_weight_parts(:) = covg_weight_parts(:) / sum(covg_weight_parts)
- 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(1) = ugsq_str / covm_weight_parts(1)
+ covg_weight_parts(2) = ugsq_ts / covm_weight_parts(2)
+ covg_weight_parts(3) = ugsq_xs / covm_weight_parts(3)
+ covg_weight_parts(4) = ugsq_ys / covm_weight_parts(4)
+
+ ! COMMENT OUT if you want the gradient-norm parts to sum to ONE
covg_weight_parts(:) = covg_weight_parts(:) / sum(covg_weight_parts)
+
endif
+ ! TESTING
+ !covm_weight_parts(:) = 1.0
+ !covg_weight_parts(:) = 1.0
+
! 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
- !if(fac_ts < EPS) fac_ts = 1.0
- !if(fac_xs < EPS) fac_xs = 1.0
- !if(fac_ys < EPS) fac_ys = 1.0
+ !if (fac_str < EPS) fac_str = 1.0
+ !if (fac_ts < EPS) fac_ts = 1.0
+ !if (fac_xs < EPS) fac_xs = 1.0
+ !if (fac_ys < EPS) fac_ys = 1.0
! Because we do not have perfect coverage, we need to adjust the normalization such that
! a perfectly recovered (source or structure) model gives a norm of about 1.0.
@@ -1995,7 +2036,7 @@
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
+ ! STRUCTURE: (fac_str / ugsq_str) * ugsq_str --> fac_str
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)
@@ -2013,7 +2054,7 @@
! possible stopping criteria based on fitting the data
chi_data_stop = 0.0
- if(ADD_DATA_ERRORS) chi_data_stop = 0.5
+ if (ADD_DATA_ERRORS) chi_data_stop = 0.5
! write model covariance values to file
open(unit=19,file=trim(out_dir2)//'sigma_values.dat',status='unknown')
@@ -2028,28 +2069,28 @@
do i=1,NVAR
write(19,*) 'covm_weight_parts', covm_weight_parts(i)
enddo
- write(19,*) 'sum(covm_weight_parts)', sum(covm_weight_parts)
+ 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,*) '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,*) 'ugsq_str', ugsq_str
+ write(19,*) 'ugsq_ts', ugsq_ts
+ write(19,*) 'ugsq_xs', ugsq_xs
+ write(19,*) 'ugsq_ys', ugsq_ys
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 / ugsq_str) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_ts', (fac_ts / ugsq_ts) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_xs', (fac_xs / ugsq_xs) * fac_total
+!!$ write(19,*) 'cov_imetric_fac_ys', (fac_ys / ugsq_ys) * fac_total
close(19)
open(unit=19,file=trim(out_dir2)//'scaling_values.dat',status='unknown')
@@ -2066,11 +2107,11 @@
!!$ 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,'(4e14.6)') ugsq_str, ugsq_ts, ugsq_xs, ugsq_ys
+!!$ !write(19,'(1e14.6)') ugsq_src
+!!$ write(19,'(4e14.6)') fac_str/ugsq_str, fac_ts/ugsq_ts, fac_xs/ugsq_xs, fac_ys/ugsq_ys
+!!$ write(19,'(4e14.6)') (fac_str/ugsq_str)*fac_total, (fac_ts/ugsq_ts)*fac_total, &
+!!$ (fac_xs/ugsq_xs)*fac_total, (fac_ys/ugsq_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
@@ -2078,7 +2119,7 @@
open(unit=19,file=trim(out_dir2)//'cov_model_diagonal.dat',status='unknown')
do i = 1,nmod
- write(19,'(4e20.10)') cov_model(i), cov_imetric(i)
+ write(19,'(2e20.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)
@@ -2129,15 +2170,15 @@
allocate(cov_data(nmeas))
cov_data(:) = 0.0
- if(IKER==0) then
+ if (IKER==0) then
cov_data(:) = SIGMA_WAVEFORM * SIGMA_WAVEFORM * nmeas_run
- elseif(IKER==1) then
+ elseif (IKER==1) then
cov_data(:) = SIGMA_DT * SIGMA_DT * nmeas_run
- elseif(IKER==2) then
+ elseif (IKER==2) then
cov_data(:) = SIGMA_DLNA * SIGMA_DLNA * nmeas_run
endif
- if(IKER==1) then
+ if (IKER==1) then
open(unit=19,file=trim(out_dir2)//'scaling_values_covd.dat',status='unknown')
write(19,*) 'ievent_min', ievent_min
write(19,*) 'ievent_max', ievent_max
@@ -2155,9 +2196,9 @@
close(19)
! load measurement perturbations
- if( ADD_DATA_ERRORS ) then
+ if ( ADD_DATA_ERRORS ) then
- if(nmeas .ge. 10000) stop 'measurement error file only contains 10000 lines'
+ if (nmeas .ge. 10000) stop 'measurement error file only contains 10000 lines'
allocate(measure_pert_vec(nmeas))
measure_pert_vec(:) = 0.0
@@ -2184,7 +2225,7 @@
!nrec = 1
- if(ISOURCE_LOG) open(91,file=trim(out_dir2)//'source_vector.log',status='unknown')
+ if (ISOURCE_LOG) open(91,file=trim(out_dir2)//'source_vector.log',status='unknown')
var_red_val = 100.0 ! initialize
chi_k_val = 1.0d10 ! initialize to a very high value
@@ -2199,9 +2240,9 @@
irun = irun0 + istep
print *,'=============================================================='
print *,' istep, imod, irun : ', istep, imod, irun
- if(INV_STRUCT_BETA==1) print *, ' inverting for structure parameters'
- if(INV_SOURCE_T==1) print *, ' inverting for source origin times'
- if(INV_SOURCE_X==1) print *, ' inverting for source locations'
+ if (INV_STRUCT_BETA==1) print *, ' inverting for structure parameters'
+ if (INV_SOURCE_T==1) print *, ' inverting for source origin times'
+ if (INV_SOURCE_X==1) print *, ' inverting for source locations'
print *,'=============================================================='
if (READ_IN) then
@@ -2219,15 +2260,15 @@
! KEY: obtain the structure and source arrays from the model vector:
! 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)
+ if (itest==0) then ! reference model (current model)
call mvec2local(nmod, nmod_src, m0_vec, beta_syn, m_src_syn)
- elseif(itest==1) then ! test model
+ elseif (itest==1) then ! test model
call mvec2local(nmod, nmod_src, mt_vec, beta_syn, m_src_syn)
endif
! for CG algorithm, update kappa_syn and mu_syn
- if( .not. READ_IN) then
+ if ( .not. READ_IN) then
!kappa_syn = rho_syn * bulk_syn * bulk_syn
mu_syn = rho_syn * beta_syn * beta_syn
endif
@@ -2235,7 +2276,7 @@
! read in the sources from another file
! NOTE: FOR NOW, WE DO NOT LOAD THE DATA SOURCES -- THEY SHOULD BE IDENTICAL.
! (In the future, we might want to modify this to read in ANY data sources.)
- if( READ_IN .and. (INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1)) then
+ if ( READ_IN .and. (INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1)) then
write(filename,'(a,i4.4,a)') trim(out_dir2)//'src_syn_m',hmod,'.dat'
open(unit=18,file=filename,status='unknown')
@@ -2263,96 +2304,68 @@
!--------------------------
! compute model norm term in the misfit function
- ! NOTE 1: TWO VERSIONS using the balanced Cm (cov_imetric) and the original Cm (cov_model)
+ ! NOTE 1: THREE possible outputs: cov_model, cov_model with covm weights, cov_model with covg weights
! NOTE 2: THE FACTOR OF 0.5 IS NOT INCLUDED HERE
- ! NOTE 3: This is norm-squared (no square root is taken)
+ ! NOTE 3: These variables are norm-squared (no square root is taken)
imnorm = 1 ! =1 for m-mprior ; =0 for g
! model vector, cov_model
- if(itest==0) then ! reference model (current model)
- call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ filename = trim(out_dir1)//'model_norm_all.dat'
+ if (itest==0) then ! reference model (current model)
+ call compute_norm_sq(filename, 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(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ elseif (itest==1) then ! test model
+ call compute_norm_sq(filename, 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_cov_model.dat'
- call write_norm_sq(filename,model_norm_parts)
+ !filename = trim(out_dir1)//'model_norm_all.dat'
+ !call write_norm_sq(filename,model_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)
+!!$ ! this output should exactly match gradient_norm_model_all.dat
+!!$ call compute_norm_sq(filename, imnorm, &
+!!$ ievent_min, ievent_max, nevent, index_source, nmod, &
+!!$ m0, mprior, cov_imetric, model_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, &
+ filename = trim(out_dir1)//'model_norm_target_all.dat'
+ call compute_norm_sq(filename, 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)
+ !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
+ if (itest==0) then
mdiff = mtarget - m0 + mprior
- elseif(itest==1) then
+ elseif (itest==1) then
mdiff = mtarget - mt + mprior
endif
! diff vector, cov_model
- call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ filename = trim(out_dir1)//'model_norm_diff_all.dat'
+ call compute_norm_sq(filename, 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)
+ !call write_norm_sq(filename,model_norm_diff_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'
+ !stop 'TESTING MODEL NORMS'
!--------------------------
@@ -2401,9 +2414,9 @@
! get the lat-lon of the TARGET SOURCES
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
- if(ISRC_SPACE <= 5) then ! if you do NOT want a point source from the event list
+ if (ISRC_SPACE <= 5) then ! if you do NOT want a point source from the event list
x_src_lon0(:) = 0.0 ; z_src_lat0(:) = 0.0
x_src0(:) = 0.0 ; z_src0(:) = 0.0
@@ -2487,7 +2500,7 @@
do xmesh = xcen-s_radius, xcen+s_radius, dx
do zmesh = zcen-s_radius, zcen+s_radius, dx
d = sqrt((xmesh - xcen)**2+(zmesh - zcen)**2)
- if(d < s_radius) then
+ if (d < s_radius) then
i = i+1
x_src0(i) = xmesh
z_src0(i) = zmesh
@@ -2507,7 +2520,7 @@
! make sure that there are fewer target points than the max allowed
! (this is not ideal, since you might have a file of points that extend far outside the grid)
- if(nsrc > MAX_SR) then
+ if (nsrc > MAX_SR) then
print *
print *, ' ISRC_SPACE = ', ISRC_SPACE
print *, ' nsrc = ', nsrc
@@ -2527,7 +2540,7 @@
! filter target points (utm-mesh) -- update nsrc
call station_filter(nsrc, x_src0, z_src0, ifilter_src, SOURCE_GRID_BUFFER)
- if(nsrc < 1) stop 'Must have at least one source'
+ if (nsrc < 1) stop 'Must have at least one source'
! allocate vectors
allocate(x_src(nsrc),z_src(nsrc),x_src_lon(nsrc),z_src_lat(nsrc))
@@ -2586,13 +2599,13 @@
endif
else
- if(ISRC_SPACE /= 6) stop ' for body waves, we only select a point source from the specified events'
+ if (ISRC_SPACE /= 6) stop ' for body waves, we only select a point source from the specified events'
nsrc = 1
endif ! ISURFACE==1
! allocate vectors (sources for data)
- if(ISURFACE==1) allocate(x_src_lon_dat(nsrc),z_src_lat_dat(nsrc))
+ if (ISURFACE==1) allocate(x_src_lon_dat(nsrc),z_src_lat_dat(nsrc))
allocate(x_src_dat(nsrc),z_src_dat(nsrc))
allocate(sglob_dat(nsrc))
allocate(ispec_src_dat(nsrc),xi_src_dat(nsrc),gamma_src_dat(nsrc))
@@ -2600,8 +2613,8 @@
! allocate 1-D Lagrange interpolators and derivatives
allocate(hxisd_store(nsrc,NGLLX),hgammasd_store(nsrc,NGLLZ))
- if(ISRC_SPACE <= 5) then ! if you do NOT want a point source from the event list
- if(ISURFACE==0) stop ' for body waves, we only select a point source from the specified events'
+ if (ISRC_SPACE <= 5) then ! if you do NOT want a point source from the event list
+ if (ISURFACE==0) stop ' for body waves, we only select a point source from the specified events'
! use same source for data and synthetics
sglob_dat(1:nsrc) = sglob(1:nsrc) ! closest gridpoint
@@ -2626,7 +2639,7 @@
x_src_dat(1) = x_eve_dat(ievent)
z_src_dat(1) = z_eve_dat(ievent)
origin_time_dat = otime_dat(ievent)
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
x_src_lon_dat(1) = x_eve_lon_dat(ievent)
z_src_lat_dat(1) = z_eve_lat_dat(ievent)
endif
@@ -2649,7 +2662,7 @@
itemp3 = index_source(3,ievent)
! allocate vectors (sources for synthetics)
- if(ISURFACE==1) allocate(x_src_lon(nsrc),z_src_lat(nsrc))
+ if (ISURFACE==1) allocate(x_src_lon(nsrc),z_src_lat(nsrc))
allocate(x_src(nsrc),z_src(nsrc))
allocate(sglob(nsrc))
allocate(ispec_src(nsrc),xi_src(nsrc),gamma_src(nsrc))
@@ -2662,7 +2675,7 @@
! filter target points (utm-mesh) -- update nsrc
call station_filter(nsrc, x_src, z_src, ifilter_src, SOURCE_GRID_BUFFER)
- if(nsrc /= 1) stop 'Must be a single point source'
+ if (nsrc /= 1) stop 'Must be a single point source'
! determine the (eta, xi) corresponding to the target points
! this UPDATES x_src, z_src; sglob is the index of the closest gridpoint
@@ -2691,7 +2704,7 @@
!-------------------
! source for data and source for synthetic
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! convert from mesh to lat-lon
call mesh_geo(nsrc, x_src_lon, z_src_lat, x_src, z_src, UTM_PROJECTION_ZONE, IMESH2LONLAT)
@@ -2706,7 +2719,7 @@
write(*,'(i8, 5f17.10)') i, temp1, temp2, temp3, temp4, temp5*6371.0*1000.0
enddo
- elseif(ISURFACE==0) then
+ elseif (ISURFACE==0) then
print *
print *, 'sources [x_src_dat, z_src_dat, x_src, z_src, dist (m)]:'
do i = 1,nsrc
@@ -2722,7 +2735,7 @@
endif ! ISRC_SPACE
! source log file
- if(ISOURCE_LOG .and. ISRC_SPACE==6) then
+ if (ISOURCE_LOG .and. ISRC_SPACE==6) then
! indices into source vector
itemp1 = index_source(1,ievent)
@@ -2757,12 +2770,12 @@
! source time function FOR DATA AND SYNTHETICS
! source magnitude (same for data and synthetics)
- if(ISURFACE==0) then
+ if (ISURFACE==0) then
! DEBUG ARRAY SIZE
f0(1) = FNORM * FOR_X
f0(2) = FNORM * FOR_Y
f0(3) = FNORM * FOR_Z
- elseif(ISURFACE==1) then
+ elseif (ISURFACE==1) then
f0(1) = FNORM
else
stop 'NCOMP must be 1 or 3'
@@ -2844,7 +2857,7 @@
deallocate(x_src,z_src)
deallocate(x_src_dat,z_src_dat)
- if(ISURFACE==1) deallocate(x_src_lon,z_src_lat,x_src_lon_dat,z_src_lat_dat)
+ if (ISURFACE==1) deallocate(x_src_lon,z_src_lat,x_src_lon_dat,z_src_lat_dat)
!stop 'testing'
@@ -2939,18 +2952,18 @@
endif ! ifirst_event
! source time functions for data and synthetics
- if(WRITE_STF_F) call write_seismogram(samp_dat, nsrc, trim(out_dir)//'stffor_dat')
- if(WRITE_STF_F) call write_seismogram(samp, nsrc, trim(out_dir)//'stffor_syn')
+ if (WRITE_STF_F) call write_seismogram(samp_dat, nsrc, trim(out_dir)//'stffor_dat')
+ if (WRITE_STF_F) call write_seismogram(samp, nsrc, trim(out_dir)//'stffor_syn')
! STOPPING CRITERIA : exit program if structure parameters are unrealistic
! for either the test model or present model
-!!$ if( maxval( m0*m0 / (SIGMA_FAC**2 * cov_imetric) ) > 1. ) then
+!!$ if ( maxval( m0*m0 / (SIGMA_FAC**2 * cov_imetric) ) > 1. ) then
!!$ print *, ' maxval( m0*m0 / (SIGMA_FAC**2 * cov_imetric) ) > 1. ) : ', &
!!$ maxval( m0*m0 / (SIGMA_FAC**2 * cov_imetric) )
!!$ print *, ' index : ', maxloc( m0*m0 / (SIGMA_FAC**2 * cov_imetric) )
!!$ stop ' STOP: model values are too extreme'
!!$ endif
-!!$ if( maxval( mt*mt / (SIGMA_FAC**2 * cov_imetric) ) > 1. ) then
+!!$ if ( maxval( mt*mt / (SIGMA_FAC**2 * cov_imetric) ) > 1. ) then
!!$ print *, ' maxval( mt*mt / (SIGMA_FAC**2 * cov_imetric) ) > 1. ) : ', &
!!$ maxval( mt*mt / (SIGMA_FAC**2 * cov_imetric) )
!!$ print *, ' index : ', maxloc( mt*mt / (SIGMA_FAC**2 * cov_imetric) )
@@ -2959,13 +2972,13 @@
!!$ vel_min = beta_min / 2.
!!$ vel_max = alpha_max * 2.
-!!$ if(itest==1) then
-!!$ if( minval(mt_vec(1:nmod_str)) < vel_min .or. maxval(mt_vec(1:nmod_str)) > vel_max ) then
+!!$ if (itest==1) then
+!!$ if ( minval(mt_vec(1:nmod_str)) < vel_min .or. maxval(mt_vec(1:nmod_str)) > vel_max ) then
!!$ print *, minval(mt_vec(1:nmod_str)), maxval(mt_vec(1:nmod_str))
!!$ stop ' STOP: test model for structure is too extreme'
!!$ endif
!!$ else
-!!$ if( minval(m0_vec(1:nmod_str)) < vel_min .or. maxval(m0_vec(1:nmod_str)) > vel_max ) then
+!!$ if ( minval(m0_vec(1:nmod_str)) < vel_min .or. maxval(m0_vec(1:nmod_str)) > vel_max ) then
!!$ print *, minval(m0_vec(1:nmod_str)), maxval(m0_vec(1:nmod_str))
!!$ stop ' STOP: current model for structure is too extreme'
!!$ endif
@@ -3001,7 +3014,7 @@
! write out seismograms at the receivers
data_tag = 'dat'
- if(WRITE_SEISMO_F) call write_seismogram(data, nrec, trim(out_dir)//data_tag)
+ if (WRITE_SEISMO_F) call write_seismogram(data, nrec, trim(out_dir)//data_tag)
!stop 'testing simulation for data'
@@ -3032,9 +3045,9 @@
! write out seismograms at the receivers
syn_tag = 'syn'
!syn_tag = 'forward'
- if(WRITE_SEISMO_F) call write_seismogram(syn, nrec, trim(out_dir)//syn_tag)
+ if (WRITE_SEISMO_F) call write_seismogram(syn, nrec, trim(out_dir)//syn_tag)
-!!$ if(WRITE_SPECTRAL_MAP_F) then
+!!$ if (WRITE_SPECTRAL_MAP_F) then
!!$ print *, 'compute and write out forward spectral map '
!!$ call write_spectral_map(syn, nrec, rglob, trim(out_dir)//'spectral_forward',WRITE_SPECTRA_F)
!!$ endif
@@ -3051,7 +3064,7 @@
tstart(:) = 0.0
tend(:) = NSTEP*DT
- if(ISURFACE==1) then
+ if (ISURFACE==1) then
! compute time windows for measurements based on a HOMOGENEOUS reference model
! Window with should be wide enough to capture synthetic and data waveforms,
! but narrow enough to exclude suprious boundary reflections.
@@ -3075,7 +3088,7 @@
!stop 'testing'
- if(WRITE_SEISMO_RECONSTRUCT) then
+ if (WRITE_SEISMO_RECONSTRUCT) then
allocate(data_recon(NSTEP,NCOMP,nrec))
data_recon = 0.0
endif
@@ -3084,31 +3097,31 @@
!call make_adjoint_source(nrec, syn, tstart, tend, adj_syn, data)
call mtm_adj(ievent, nrec, syn, tstart, tend, adj_syn, data, data_recon)
- if(WRITE_SEISMO_RECONSTRUCT) call write_seismogram(data_recon, nrec, trim(out_dir)//'dat_recon')
+ if (WRITE_SEISMO_RECONSTRUCT) call write_seismogram(data_recon, nrec, trim(out_dir)//'dat_recon')
! specify which components you want to send back
- if(ISURFACE==0) then
- if(REV_X == 0) adj_syn(:,1,:) = 0.0
- if(REV_Y == 0) adj_syn(:,2,:) = 0.0
- if(REV_Z == 0) adj_syn(:,3,:) = 0.0
+ if (ISURFACE==0) then
+ if (REV_X == 0) adj_syn(:,1,:) = 0.0
+ if (REV_Y == 0) adj_syn(:,2,:) = 0.0
+ if (REV_Z == 0) adj_syn(:,3,:) = 0.0
endif
! write out adjoint source time function at the receivers
stfadj_tag = 'stfadj'
- if(WRITE_STF_A) call write_seismogram(adj_syn, nrec, trim(out_dir)//stfadj_tag)
+ if (WRITE_STF_A) call write_seismogram(adj_syn, nrec, trim(out_dir)//stfadj_tag)
! OUTPUT ASCII FILES --> SAC FILES (make_sac_files.pl)
! (1) data, (2) synthetics, (3) adjoint source time function
-!!$ if(WRITE_SEISMO_F) then
+!!$ if (WRITE_SEISMO_F) then
!!$ filename1 = trim(script_dir)//'make_sac_files.csh'
!!$ filename2 = 'make_sac_files.pl'
!!$ open(19,file=filename1,status='unknown')
-!!$ if(IKER <= 4) write(19,'(7a,f12.6)') trim(script_dir)//trim(filename2),' ', &
+!!$ if (IKER <= 4) write(19,'(7a,f12.6)') trim(script_dir)//trim(filename2),' ', &
!!$ trim(out_dir),' ', trim(data_tag) ,' ','1', tshift
!!$ write(19,'(7a,f12.6)') trim(script_dir)//trim(filename2),' ', &
!!$ trim(out_dir),' ', trim(syn_tag) ,' ','1', tshift
-!!$ if(WRITE_STF_A) write(19,'(7a,f12.6)') trim(script_dir)//trim(filename2),' ', &
+!!$ if (WRITE_STF_A) write(19,'(7a,f12.6)') trim(script_dir)//trim(filename2),' ', &
!!$ trim(out_dir),' ', trim(stfadj_tag),' ','1', tshift
!!$ close(19)
!!$ call system('chmod 755 scripts/make_sac_files.csh ; scripts/make_sac_files.csh')
@@ -3118,7 +3131,7 @@
print *,'---------------------------'
print *,' misfit for IKER = ',IKER
print *,' istep, imod, irun, ievent : ', istep, imod, irun, ievent
- print *,' chi-data = ', sum(chi_data(ievent,:,:,1))
+ print *,' data_norm2(ievent) = ', sum(chi_data(ievent,:,:,1))
print *,'---------------------------'
!!$ ! write ALL adjoint sources to file
@@ -3147,7 +3160,7 @@
! we always evaluate the kernels for the present model
! we only evaluate the kernel for the test model if itest==1 or POLY_ORDER==3
- if( (itest==0 .or. POLY_ORDER==3) .and. COMPUTE_KERNELS) then
+ if ( (itest==0 .or. POLY_ORDER==3) .and. COMPUTE_KERNELS) then
print *
print *, 'compute the kernel via adjoint wavefield interaction'
print *
@@ -3197,7 +3210,7 @@
!----------------------
! smooth EACH event kernel for the data subspace method
- if(ISMOOTH_EVENT_KERNEL == 1) then
+ if (ISMOOTH_EVENT_KERNEL == 1) then
call smooth_function(btype_kernel, GAMMA_SMOOTH_KERNEL, kbeta_smooth)
!call smooth_function(atype_kernel, GAMMA_SMOOTH_KERNEL, kbulk_smooth)
@@ -3207,7 +3220,7 @@
endif
!!$ ! smooth EACH event kernel for the data subspace method
-!!$ if(ISMOOTH_EVENT_KERNEL == 1) then
+!!$ if (ISMOOTH_EVENT_KERNEL == 1) then
! write smoothed beta kernel to file
open(unit = 13, file = trim(out_dir)//'kernel_basis_smooth', status = 'unknown',iostat=ios)
@@ -3227,7 +3240,7 @@
enddo
close(13)
- ! store (unbalanced) norm of gradient
+ ! store UNBALANCED norm of gradient -- cov_model
itemp1 = index_source(1,ievent)
itemp2 = index_source(2,ievent)
itemp3 = index_source(3,ievent)
@@ -3250,13 +3263,13 @@
deallocate(atype_kernel,btype_kernel,rtype_kernel)
! write out adjoint seismograms at the original sources
- if(WRITE_SEISMO_A) call write_seismogram(samp, nsrc, trim(out_dir)//'synadj')
+ if (WRITE_SEISMO_A) call write_seismogram(samp, nsrc, trim(out_dir)//'synadj')
endif
! deallocate variables associated with one event
- if(WRITE_SEISMO_RECONSTRUCT) deallocate(data_recon)
+ if (WRITE_SEISMO_RECONSTRUCT) deallocate(data_recon)
deallocate(data, adj_syn, absorb_field)
deallocate(syn, tstart, tend)
@@ -3293,7 +3306,7 @@
! VARIANCE REDUCTION
! NOTE 1: if the new chi (chi_val) is larger than the previous chi (chi_k_val), then VR < 0.
! NOTE 2: consider whether chi_data_stop should be removed
- if(istep >= 2) then
+ if (istep >= 2) then
!var_red_val = 100.0 * ( 1.0 - ( (chi_val - chi_data_stop)**2 / (chi_k_val - chi_data_stop)**2 ) )
!var_red_val = -log( abs(chi_val - chi_data_stop) / abs(chi_k_val - chi_data_stop) )
var_red_val = log( chi_k_val / chi_val )
@@ -3314,14 +3327,14 @@
close(18)
deallocate(measure_vec)
- if(ISOURCE_LOG) then
+ if (ISOURCE_LOG) then
write(91,*) 'SUMMED MISFIT FOR ALL EVENTS'
write(91,'(a30,1i16)') ' N-data : ', nmeas_run
- 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) : ', &
+ write(91,'(a30,1e16.8)') ' data_norm2(m) : ', data_norm
+ write(91,'(a30,3e16.8)') ' model_norm2(m) : ', model_norm, model_norm_struct, model_norm_source
+ write(91,'(a30,3e16.8)') ' model_norm2_target(m) : ', &
model_norm_target, model_norm_target_struct, model_norm_target_source
- write(91,'(a30,3e16.8)') ' model_norm_diff(m) : ', &
+ write(91,'(a30,3e16.8)') ' model_norm2_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
@@ -3331,13 +3344,13 @@
! STOPPING CRITERIA for CG algorithm
! Program will stop after files are written (search stop_cg)
stop_cg = .false.
- if( data_norm <= 2.0 * chi_data_stop) then
+ if ( data_norm <= 2.0 * chi_data_stop) then
stop_cg = .true.
print *, ' data_norm <= 2(chi_data_stop) : ', data_norm, 2.0*chi_data_stop
endif
- !if( itest==0 .and. (chi_val .ge. chi_k_val) ) then
- if( itest==0 .and. var_red_val < VAR_RED_MIN ) then
+ !if ( itest==0 .and. (chi_val .ge. chi_k_val) ) then
+ if ( itest==0 .and. var_red_val < VAR_RED_MIN ) then
stop_cg = .true.
print *, ' chi_k_val : ', chi_k_val
print *, ' chi_val : ', chi_val
@@ -3346,8 +3359,8 @@
print *, ' var_red_val < VAR_RED_MIN'
endif
-!!$ if(istep == 0) chi_val_0 = chi_val
-!!$ if( chi_val <= chi_val_0 * CONV_STOP ) then
+!!$ if (istep == 0) chi_val_0 = chi_val
+!!$ if ( chi_val <= chi_val_0 * CONV_STOP ) then
!!$ print *, ' chi_val : ', chi_val
!!$ print *, ' chi_val_0 * CONV_STOP : ', chi_val_0 * CONV_STOP
!!$ print *, ' chi_val <= chi_val_0 * CONV_STOP'
@@ -3356,7 +3369,7 @@
! enter this section if you want to compute the gradient and/or iterate
- if( COMPUTE_KERNELS ) then
+ if ( COMPUTE_KERNELS ) then
print *, ' Now we compute the gradient of the misfit function (using the misfit kernel)'
@@ -3381,7 +3394,7 @@
! COMPUTE THE GRADIENT OF THE MISFIT FUNCTION, if the present model is not
! a test model or if the CG polynomial is a cubic function
- if( itest==0 .or. POLY_ORDER==3 ) then
+ 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 ! initialized above
@@ -3391,7 +3404,7 @@
write(19,'(1e20.10)') (source_gradient(i), i = 1,nmod_src)
close(19)
- if(INV_STRUCT_BETA == 1) then
+ if (INV_STRUCT_BETA == 1) then
! norm of parts of unbalanced gradient for all events
! NOTE: within each event directory there is also gradient_data_unbalanced.dat
@@ -3413,42 +3426,45 @@
! Smooth the summed event kernel (= misfit kernel) to remove spurious src-rec effects;
! this is done via Gaussian convolution.
- if( ISMOOTH_MISFIT_KERNEL == 1) then
+ if ( ISMOOTH_MISFIT_KERNEL == 1) then
call smooth_function(btype_kernel_sum, GAMMA_SMOOTH_KERNEL, kbeta_smooth)
else
kbeta_smooth = btype_kernel_sum
endif
- kbulk_smooth = 0.0 ! testing for S-wave only
+ !kbulk_smooth = 0.0 ! testing for S-wave only
endif ! INV_STRUCT_BETA
!------------------------
! assemble the gradient vector
- ! NOTE: structure gradient g_i = K_i * dA_i
+ ! NOTE 1: gradient of data term of misfit function
+ ! NOTE 2: structure gradient g_i = K_i * dA_i
temp_local1 = kbeta_smooth * da_local
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
+ ! KEY QUESTION: cov_model or cov_imetric? cov_model appears to perform slightly better
+ if ( INCLUDE_MODEL_NORM ) then
if (itest == 0) then
- gradient_model(:) = (m0(:) - mprior(:)) / cov_model(:)
+ gradient_model(:) = (m0(:) - mprior(:)) / cov_imetric(:)
+ !gradient_model(:) = (m0(:) - mprior(:)) / cov_model(:)
else
- gradient_model(:) = (mt(:) - mprior(:)) / cov_model(:)
+ gradient_model(:) = (mt(:) - mprior(:)) / cov_imetric(:)
+ !gradient_model(:) = (mt(:) - mprior(:)) / cov_model(:)
endif
endif
! KEY: if NOT inverting for structure or source, then set those parts of the gradient to zero
- if(INV_STRUCT_BETA == 0) then
+ 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
endif
- if(INV_SOURCE_T == 0) then
+ if (INV_SOURCE_T == 0) then
gradient_data(m_inds(2,1) : m_inds(2,2)) = 0.0
gradient_model(m_inds(2,1) : m_inds(2,2)) = 0.0
endif
- if(INV_SOURCE_X == 0) then
+ if (INV_SOURCE_X == 0) then ! both xs and ys
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
@@ -3456,6 +3472,19 @@
! overall gradient
gradient = gradient_data + gradient_model
+ ! NORMS OF THE GRADIENT (C^-1 is used, not C ; mprior is not used)
+ imnorm = 0
+!!$
+!!$ ! normalize -- no cov weights
+!!$ call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+!!$ gradient, mprior, cov_model, gradient_norm_parts)
+!!$ gnorm = sqrt( sum( gradient_norm_parts ) )
+!!$ print *, 'normalizing gradient by ', gnorm
+!!$ gradient(:) = gradient(:) / gnorm
+!!$ filename1 = trim(out_dir1)//'gradient_norm_all_unit.dat'
+!!$ call write_norm_sq(filename1,gradient_norm_parts)
+!!$ !stop 'testing gradient normalization'
+
! write gradient vector to file (source and structure parameters)
!write(19,'(1e20.10)') (gradient(i), i = 1,nmod)
open(unit=19,file=trim(out_dir1)//'gradient_vec.dat',status='unknown')
@@ -3464,30 +3493,36 @@
enddo
close(19)
- ! NORMS OF THE GRADIENT (C^-1 is used, not C ; mprior is not used)
- imnorm = 0
+ 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'
! compute the balanced gradient norm
- if (0==1) then
- call compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
+ ! NOTE: these two options should be equivalent (check to be sure)
+ if (1==1) then ! cov_imetric
+ call compute_norm_sq(filename1, 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, &
+ call compute_norm_sq(filename2, 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, &
+ call compute_norm_sq(filename3, 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, &
+ ! NOTE: gnorm( Cinv(m-mprior) ) should equal mnorm( m-mprior )
+
+ else ! cov_model plus weights
+ call compute_norm_sq(filename1, 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, &
+ call compute_norm_sq(filename2, 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, &
+ call compute_norm_sq(filename3, 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 )
@@ -3495,36 +3530,74 @@
!!$ 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_norm_data_parts)
- call write_norm_sq(filename3,gradient_norm_model_parts)
+!!$ ! 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_norm_data_parts)
+!!$ call write_norm_sq(filename3,gradient_norm_model_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, &
+ ! SAME AS ABOVE, using cov_model WITHOUT weights
+ 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 compute_norm_sq(filename1, 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, &
+ call compute_norm_sq(filename2, 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, &
+ call compute_norm_sq(filename3, 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_norm_data_parts)
+!!$ call write_norm_sq(filename3,gradient_norm_model_parts)
- ! 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_norm_data_parts)
- call write_norm_sq(filename3,gradient_norm_model_parts)
+ !stop 'TESTING GRADIENT NORMS'
+
+ ! update balance using unweighted gradient norm -- see cov_imetric above
+ ! NOTE 1: this will only use gradient for non-test models (itest=0)
+ ! NOTE 2: do you want this for the first step only (istep=0) or each step
+ if ( (istep==0) .and. (0==1) ) then
+ ! NOTE: choose to balance using the full unbalanced gradient or the data term only
+ gbalance_parts(:) = 0.0
+ gbalance_parts = gradient_norm_data_parts
+ !gbalance_parts = gradient_norm_parts
+
+ open(unit=19,file=trim(out_dir1)//'cov_imetric_weights.dat',status='unknown')
+ do i = 1,NVAR
+ write(19,'(3e16.6)') gbalance_parts(i), covm_weight_parts(i), &
+ gbalance_parts(i) / covm_weight_parts(i)
+ enddo
+ close(19)
+
+ cov_imetric(m_inds(1,1):m_inds(1,2)) = cov_model(m_inds(1,1):m_inds(1,2)) &
+ / (gbalance_parts(1) / covm_weight_parts(1))
+ cov_imetric(m_inds(2,1):m_inds(2,2)) = cov_model(m_inds(2,1):m_inds(2,2)) &
+ / (gbalance_parts(2) / covm_weight_parts(2))
+ cov_imetric(m_inds(3,1):m_inds(3,2)) = cov_model(m_inds(3,1):m_inds(3,2)) &
+ / (gbalance_parts(3) / covm_weight_parts(3))
+ cov_imetric(m_inds(4,1):m_inds(4,2)) = cov_model(m_inds(4,1):m_inds(4,2)) &
+ / (gbalance_parts(4) / covm_weight_parts(4))
+
+ ! write the new cov_imetric to file (over-write previous file)
+ open(unit=19,file=trim(out_dir2)//'cov_model_diagonal.dat',status='unknown')
+ do i = 1,nmod
+ write(19,'(2e20.10)') cov_model(i), cov_imetric(i)
+ enddo
+ close(19)
+ endif
+
!!$ ! scaling for joint inversions
-!!$ if( istep==0 .and. 0==1 ) then
+!!$ if ( istep==0 .and. 0==1 ) then
!!$
!!$ ! scaling for joint source-structure inversions
!!$
@@ -3584,7 +3657,7 @@
print *, ' Entering CG algorithm to compute new model or test model'
- if(itest==0) then ! if the present kernel is for a REAL model
+ if (itest==0) then ! if the present kernel is for a REAL model
chi_k_val = chi_val
gk(:) = gradient(:)
@@ -3592,7 +3665,7 @@
!-----------------------------------------
! update the current search direction p0 (p0 and pk are NON-hat)
- if(istep == 0) then
+ if (istep == 0) then
beta_val = 0.0
p0(:) = 0.0 ! initialize
else
@@ -3603,7 +3676,7 @@
pk(:) = -cov_imetric(:) * gk(:) + beta_val * p0(:)
! write gradient vectors (before re-assigning g0 and gk)
- if(1==1) then
+ if (1==1) then
open(unit=19,file=trim(out_dir1)//'cg_grad_vectors.dat',status='unknown')
do i = 1,nmod
write(19,'(4e16.6)') g0(i), gk(i), p0(i), pk(i)
@@ -3630,8 +3703,8 @@
!lam_t_val = -2.0*chi_k_val / dot_product(gk, pk) ! GJI paper
!istep_switch = 6
- !if(istep < istep_switch) lam_t_val = -2.0*chi_k_val / dot_product(gk, pk) ! quadratic extrapolation
- !if(istep >= istep_switch) lam_t_val = -chi_k_val / dot_product(gk, pk) ! linear extrapolation
+ !if (istep < istep_switch) lam_t_val = -2.0*chi_k_val / dot_product(gk, pk) ! quadratic extrapolation
+ !if (istep >= istep_switch) lam_t_val = -chi_k_val / dot_product(gk, pk) ! linear extrapolation
! compute test model
mt(:) = m0(:) + lam_t_val * pk(:)
@@ -3644,7 +3717,7 @@
!!$ itemp1 = nmod_str + (i-1)*3 + 1
!!$ itemp2 = nmod_str + (i-1)*3 + 2
!!$ itemp3 = nmod_str + (i-1)*3 + 3
-!!$ if( i < ievent_min .or. i > ievent_max ) then
+!!$ if ( i < ievent_min .or. i > ievent_max ) then
!!$ mt(itemp1) = 0.0
!!$ mt(itemp2) = 0.0
!!$ mt(itemp3) = 0.0
@@ -3655,7 +3728,7 @@
mt_vec(nmod_str+1 : nmod) = mt(nmod_str+1 : nmod)
!!$ ! STRUCTURE PARAMETERS (bulk, then beta)
-!!$ if(INV_STRUCT_BETA == 0) then
+!!$ 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) )
@@ -3664,7 +3737,7 @@
!!$ endif
!!$
!!$ ! SOURCE PARAMETERS
-!!$ if(INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0) then
+!!$ 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)
@@ -3678,10 +3751,10 @@
! write scalars to file
open(unit=19,file=trim(out_dir1)//'cg_test_vals.dat',status='unknown')
- write(19,'(4e16.8)') mu_val, lam_t_val, beta_val, chi_k_val
+ write(19,'(5e16.8)') mu_val, lam_t_val_bot, lam_t_val, beta_val, chi_k_val
close(19)
- elseif(itest==1) then ! if present kernel is for a test model
+ elseif (itest==1) then ! if present kernel is for a test model
chi_t_val = chi_val
@@ -3693,7 +3766,7 @@
g1 = sum(g0(:) * pk(:))
!g1 = dot_product(g0, pk) ! GJI paper
- if(POLY_ORDER == 3) then
+ if (POLY_ORDER == 3) then
! use cubic polynomial: six values gives an analytical minimum
! see Matlab scripts cubic_min_4.m and cubic_min.m
@@ -3711,9 +3784,9 @@
! get the analytical minimum
qfac = Pb**2.0 - 3.0*Pa*Pc
- if(Pa /= 0.0 .and. qfac >= 0.0) then
+ if (Pa /= 0.0 .and. qfac >= 0.0) then
xmin = (-Pb + sqrt(qfac)) / (3.0*Pa)
- elseif(Pa == 0.0 .and. Pb /= 0.0) then
+ elseif (Pa == 0.0 .and. Pb /= 0.0) then
xmin = -Pc/(2.0*Pb)
else
print *, 'Pa, Pb, Pc, Pd, qfac: ',Pa,Pb,Pc,Pd,qfac
@@ -3736,7 +3809,7 @@
Pc = yy1 - Pa*xx1**2 - Pb*xx1
! get the analytical minimum (the vertex)
- if(Pa /= 0.0) then
+ if (Pa /= 0.0) then
xmin = -Pb / (2.0*Pa)
else
print *, 'Pa, Pb, Pc: ',Pa,Pb,Pc
@@ -3763,7 +3836,7 @@
!!$ itemp1 = nmod_str + (i-1)*3 + 1
!!$ itemp2 = nmod_str + (i-1)*3 + 2
!!$ itemp3 = nmod_str + (i-1)*3 + 3
-!!$ if( i < ievent_min .or. i > ievent_max ) then
+!!$ if ( i < ievent_min .or. i > ievent_max ) then
!!$ m0(itemp1) = 0.0
!!$ m0(itemp2) = 0.0
!!$ m0(itemp3) = 0.0
@@ -3774,7 +3847,7 @@
m0_vec(nmod_str+1 : nmod) = m0(nmod_str+1 : nmod)
!!$ ! STRUCTURE PARAMETERS
-!!$ if(INV_STRUCT_BETA == 0) then
+!!$ 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) )
@@ -3783,7 +3856,7 @@
!!$ endif
!!$
!!$ ! SOURCE PARAMETERS
-!!$ if(INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0) then
+!!$ 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)
@@ -3797,7 +3870,7 @@
! write updated/test model to file (source + structure)
! write CG vectors to file
- if(1==1) then
+ if (1==1) then
open(unit=19,file=trim(out_dir1)//'cg_model_vectors.dat',status='unknown')
do i = 1,nmod
write(19,'(4e16.6)') m0(i), mt(i), m0_vec(i), mt_vec(i)
@@ -3812,7 +3885,7 @@
!====================
! stop CG algorithm (search stop_cg)
- if( stop_cg ) then
+ if ( stop_cg ) then
stop ' convergence criteria is met.'
endif
@@ -3828,28 +3901,28 @@
print *, ' Done with LOOP 2.'
- if(ISOURCE_LOG) close(91) ! source log file
+ if (ISOURCE_LOG) close(91) ! source log file
! deallocate event and receiver variables
- if(ISURFACE==1) deallocate(fglob)
+ if (ISURFACE==1) deallocate(fglob)
deallocate(rglob)
deallocate(ispec_rec,xi_rec,gamma_rec)
deallocate(x_rec,z_rec)
deallocate(hxir_store, hgammar_store)
- if(ISURFACE==1) deallocate(x_rec_lon,z_rec_lat)
+ if (ISURFACE==1) deallocate(x_rec_lon,z_rec_lat)
deallocate(eglob_syn)
deallocate(ispec_eve_syn,xi_eve_syn,gamma_eve_syn)
deallocate(x_eve_syn,z_eve_syn)
deallocate(hxie_store,hgammae_store)
- if(ISURFACE==1) deallocate(x_eve_lon_syn,z_eve_lat_syn)
+ if (ISURFACE==1) deallocate(x_eve_lon_syn,z_eve_lat_syn)
deallocate(eglob_dat)
deallocate(ispec_eve_dat,xi_eve_dat,gamma_eve_dat)
deallocate(x_eve_dat,z_eve_dat)
deallocate(hxied_store, hgammaed_store)
- if(ISURFACE==1) deallocate(x_eve_lon_dat,z_eve_lat_dat)
+ if (ISURFACE==1) deallocate(x_eve_lon_dat,z_eve_lat_dat)
deallocate(otime_syn,otime_dat)
@@ -3868,7 +3941,7 @@
!deallocate(m_scale_src_all)
deallocate(cov_data,index_data,index_source)
- if(ADD_DATA_ERRORS) deallocate(measure_pert_vec)
+ if (ADD_DATA_ERRORS) deallocate(measure_pert_vec)
! print output
print *
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-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90 2010-02-12 20:17:28 UTC (rev 16258)
@@ -1120,8 +1120,9 @@
!----------------------------------------------------
- subroutine compute_norm_sq(imnorm, ievent_min, ievent_max, nevent, index_source, nmod, &
- mvec, mvec_prior, cov_model, norm_parts, norm_parts_weight)
+ subroutine compute_norm_sq(filename, imnorm, &
+ ievent_min, ievent_max, nevent, index_source, nmod, &
+ mvec, mvec_prior, cov_model, norm_sq_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
@@ -1129,43 +1130,48 @@
! 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
+ character(len=200), intent(in) :: filename
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, 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, dimension(NVAR), intent(out) :: norm_parts
+ double precision, dimension(NVAR), intent(out) :: norm_sq_parts
double precision, dimension(nmod) :: mtemp, ctemp
- double precision, dimension(NVAR) :: npw
- integer :: ievent, itemp1, itemp2, itemp3
+ double precision, dimension(NVAR) :: npw, norm_sq_parts0
+ integer :: i, ievent, itemp1, itemp2, itemp3
!----------
- norm_parts(:) = 0.0
+ norm_sq_parts0(:) = 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
+ ! NOTE 3: norm_parts_weight is related to the norm-squared weights (not norm)
+ if (imnorm == 0) then ! norm of gradient
ctemp(:) = 1.0 / cov_model(:)
mtemp(:) = mvec(:)
if (present(norm_parts_weight)) npw(:) = 1.0 / norm_parts_weight(:)
+ !if (present(norm_parts_weight)) npw(:) = 1.0 / norm_parts_weight(:)**2
- elseif (imnorm == 1) then
+ elseif (imnorm == 1) then ! norm of model
ctemp(:) = cov_model(:)
mtemp(:) = mvec(:) - mvec_prior(:)
if (present(norm_parts_weight)) npw(:) = norm_parts_weight(:)
+ !if (present(norm_parts_weight)) npw(:) = norm_parts_weight(:)**2
+
else
stop 'imnorm must = 0 or 1'
endif
! structure part of the norm -- BETA only
! NOTE: division corresponds to inversion of a diagonal covariance matrix
- norm_parts(1) = sum( mtemp(1 : NLOCAL)**2 / ctemp(1 : NLOCAL) )
+ norm_sq_parts0(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
@@ -1173,37 +1179,74 @@
itemp2 = NLOCAL + index_source(2,ievent)
itemp3 = NLOCAL + index_source(3,ievent)
- 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)
+ norm_sq_parts0(2) = norm_sq_parts0(2) + mtemp(itemp1)**2 / ctemp(itemp1)
+ norm_sq_parts0(3) = norm_sq_parts0(3) + mtemp(itemp2)**2 / ctemp(itemp2)
+ norm_sq_parts0(4) = norm_sq_parts0(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_sq_parts0(1)
+ !norm_source = sum( norm_sq_parts0(2:4) )
+ !norm_total = sum( norm_sq_parts0(1:4) )
! weight each part of the norm
- norm_parts(:) = norm_parts(:) * npw(:)
+ norm_sq_parts(:) = norm_sq_parts0(:) * npw(:)
- end subroutine compute_norm_sq
+ !-----------------------------------------------------
+ ! write parts to file
- !-----------------------------------------------------
-
- subroutine write_norm_sq(filename,norm_parts)
-
- double precision, dimension(NVAR), intent(in) :: norm_parts
- character(len=200), intent(in) :: filename
- integer :: i
-
- ! norm of each component
+ ! norm-squared and norm of each parameter class
open(unit=19,file=filename,status='unknown')
do i = 1,NVAR
- write(19,'(1i10,3e20.10)') i, norm_parts(i)
+ write(19,'(1i4,6e16.6)') i, &
+ norm_sq_parts0(i), sqrt(norm_sq_parts0(i)), &
+ norm_sq_parts(i), sqrt(norm_sq_parts(i)), &
+ npw(i), sqrt(npw(i))
enddo
- write(19,'(1i10,1e20.10)') 1, norm_parts(1) ! structure
- write(19,'(1i10,1e20.10)') 2, sum( norm_parts(2:4) ) ! source
- write(19,'(1i10,1e20.10)') 1, sum( norm_parts(1:4) ) ! total
- end subroutine write_norm_sq
+ ! NOTE: THE SUMS ARE MORE RELEVANT IF BALANCING THE NORM-SQUARED
+ ! structure
+ write(19,'(1i4,6e16.6)') 1, &
+ norm_sq_parts0(1), sqrt(norm_sq_parts0(1)), &
+ norm_sq_parts(1), sqrt(norm_sq_parts(1)), 0.0, 0.0
+ ! source (sum norm-squared parts only)
+ write(19,'(1i4,6e16.6)') 2, &
+ sum(norm_sq_parts0(2:4)), sqrt(sum(norm_sq_parts0(2:4))), &
+ sum(norm_sq_parts(2:4)), sqrt(sum(norm_sq_parts(2:4))), 0.0, 0.0
+ ! total (sum norm-squared parts only)
+ write(19,'(1i4,6e16.6)') 1, &
+ sum(norm_sq_parts0(1:4)), sqrt(sum(norm_sq_parts0(1:4))), &
+ sum(norm_sq_parts(1:4)), sqrt(sum(norm_sq_parts(1:4))), 0.0, 0.0
+ close(19)
+
+ end subroutine compute_norm_sq
+
+!!$ !-----------------------------------------------------
+!!$
+!!$ subroutine write_norm_sq(filename,norm_sq_parts)
+!!$
+!!$ double precision, dimension(NVAR), intent(in) :: norm_sq_parts
+!!$ character(len=200), intent(in) :: filename
+!!$ integer :: i
+!!$
+!!$ ! norm-squared and norm of each parameter class
+!!$ open(unit=19,file=filename,status='unknown')
+!!$ do i = 1,NVAR
+!!$ write(19,'(1i10,2e20.10)') i, norm_sq_parts(i), sqrt(norm_sq_parts(i))
+!!$ enddo
+!!$
+!!$ ! NOTE: THE SUMS ARE MORE RELEVANT IF BALANCING THE NORM-SQUARED
+!!$ ! structure
+!!$ write(19,'(1i10,2e20.10)') 1, norm_sq_parts(1), sqrt(norm_sq_parts(1))
+!!$ ! source (sum norm-squared parts only)
+!!$ write(19,'(1i10,2e20.10)') 2, sum(norm_sq_parts(2:4)), &
+!!$ sqrt(sum(norm_sq_parts(2:4)))
+!!$ ! total (sum norm-squared parts only)
+!!$ write(19,'(1i10,2e20.10)') 1, sum(norm_sq_parts(1:4)), &
+!!$ sqrt(sum(norm_sq_parts(1:4)))
+!!$
+!!$ close(19)
+!!$
+!!$ end subroutine write_norm_sq
+
end module wave2d_sub
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-12 18:11:16 UTC (rev 16257)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_variables.f90 2010-02-12 20:17:28 UTC (rev 16258)
@@ -132,7 +132,8 @@
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, dimension(NVAR) :: gradient_norm_parts, gradient_norm_model_parts, &
+ gradient_norm_data_parts, gbalance_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
More information about the CIG-COMMITS
mailing list