[cig-commits] r18231 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: PLOTTING matlab src
carltape at geodynamics.org
carltape at geodynamics.org
Fri Apr 15 12:59:49 PDT 2011
Author: carltape
Date: 2011-04-15 12:59:48 -0700 (Fri, 15 Apr 2011)
New Revision: 18231
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_all.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_mean.pl
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.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
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/matlab/wave2d_subspace.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
Log:
utility scripts updates for iterative inversion for 2D code
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 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2011-04-15 19:59:48 UTC (rev 18231)
@@ -26,7 +26,7 @@
#
#==========================================================
-if (@ARGV < 6) {die("Usage: plot_kernels.pl colors iker irun0 ipoly qmax HESSIAN\n");}
+if (@ARGV < 6) {die("Usage: plot_kernels.pl colors iker irun0 ipoly qmax READ_IN\n");}
($colors,$iker,$irun0,$ipoly,$qmax,$ihessian) = @ARGV;
$iter = 0;
@@ -58,6 +58,7 @@
$ddir = "${pdir}/DATA_FILES";
$idir = "${basedir}/INPUT";
$odir = "${basedir}/OUTPUT/run_";
+#$odir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_OUTPUT_OLD_2/run_"; # CHT
if (not -e $ddir) {die("Check if ddir $ddir exist or not\n");}
# misfit function variable
@@ -413,7 +414,7 @@
# replace source and structure files with those from the Hessian run
if($ihessian == 1) {
for ($h = 1; $h <= $qmax; $h = $h+1) {
- $dir = sprintf("$dir0/HESSIAN/model_m%2.2i",$h);
+ $dir = sprintf("$dir0/READ_IN/model_m%4.4i",$h);
$file1syn_str = "$dir/$mfile_syn_str";
$file1syn_src = "$dir/$mfile_syn_src";
if (not -f $file1syn_str) {die("Check if $file1syn_str exist or not\n");}
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_all.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_all.pl 2011-04-15 19:59:48 UTC (rev 18231)
@@ -0,0 +1,256 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_gauss_all.pl
+# Carl Tape
+# 15-Jan-2010
+#
+#
+# Set of 1000 models
+# ../plot_gauss_all.pl
+#
+#==========================================================
+
+#if (@ARGV < 5) {die("Usage: plot_gauss_all.pl xxx\n");}
+#($im,$imodel,$pbar) = @ARGV;
+
+$imodel = 0;
+$im = 1;
+$irun0 = $im;
+
+# 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";
+
+$pbar = "0.1/0.05";
+($pmax,$ptick) = split("/",$pbar);
+$mfile_syn = "structure_syn.dat";
+
+# directories
+$plotdir = "${basedir}/PLOTTING";
+$odir = "${basedir}/OUTPUT"; # might need to be modified
+$idir = "${basedir}/INPUT";
+$figdir = "${plotdir}/FIGURES";
+if (not -e $figdir) {die("Check if figdir $figdir exist or not\n");}
+
+$stimodel = sprintf("m%.1f",$imodel/2);
+$stirun0 = sprintf("%4.4i",$irun0);
+$stirun = sprintf("%4.4i",$irun0+$imodel);
+
+$dir0 = "$odir/run_${stirun0}";
+
+if(1==1) {
+ $dir1 = "$odir/run_${stirun}";
+} else {
+ $dir1 = sprintf("$dir0/READ_IN_CG/model_m%4.4i",$imodel);
+}
+
+$file1syn = "$dir1/${mfile_syn}";
+if (not -f $file1syn) { die("Check if $file1syn exist or not\n") }
+#if (not -f $evefile) { die("Check if $evefile exist or not\n") }
+
+# source and receivers
+$rfile = "$dir0/recs_lonlat.dat";
+$efile = "$dir0/events_syn_lonlat.dat";
+if (not -f $efile) { die("Check if efile $efile exist or not\n") }
+if (not -f $rfile) { die("Check if rfile $rfile 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");
+
+#$edir = sprintf("$dir/event_%3.3i",$event);
+#print "\n $dir,$edir,$mfile_dat,$mfile_syn,$beta0,$per,$ifinite,$iopt \n";
+
+$plabel = "${plotdir}/plot_gauss_all.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
+$src = "-W0.75p -Sa0.20";
+$rec = "-W0.5p,0/0/0 -Sc3p";
+$Wshelf = "-W1.0/0/0/0tap";
+
+# resolution of color plots
+$interp = "-I0.5m/0.5m -S4m"; # key information
+$grdfile = "temp.grd";
+
+#-------------------------
+# color
+
+# KEY: scaling for color
+$scale_color = 21.0;
+$colorbar = "seis";
+
+#-------------------------
+
+# write plotting scripts
+$wid = 3.5;
+$J1 = "-JM${wid}i"; # in lat-lon
+$origin = "-X1.5 -Y2.0";
+
+$B1dat = "-B1:.\" \":WeSN";
+$B1syn = "-B1:.\" \":wESN";
+
+$Dlen = 2.0;
+$Dx = $wid/2;
+$Dy = -0.35;
+
+$Dscale1 = "-D$Dx/$Dy/$Dlen/0.10h";
+
+#$Bscale1 = sprintf("-B%2.2e:\" Phase Velocity ( km s\@+-1\@+ )\": -E10p",$ptick);
+$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+
+#-------------------------
+
+$title = "Membrane Wave Speed";
+#$title = "Rayleigh Wave Phase Velocity";
+
+# 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 $file1syn`);
+#($tt) = sort{$b <=> $a} ($tt,abs($tmin),abs($tmax));
+$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 $smin $smax \n";
+
+# plot title
+$xtx = $xmin+0.5*($xmax-$xmin);
+$ztx = $zmin+1.10*($zmax-$zmin);
+
+ #===============================================
+ $cshfile = "plot_gauss_all.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";
+ #print CSH "makecpt -C$colorbar -D $T1 > temp1\n";
+ #print CSH "sed 's/^B.*/B 170 0 0 /' temp1 > temp2\n";
+ #print CSH "sed 's/^F.*/F 0 0 255 /' temp2 > $cptfile1\n";
+
+#=============================================
+
+ $origin = "-X1.5 -Y5.0";
+ $B1syn = "-B1:.\" \":WESN";
+
+#=============================================
+
+$pmin = 1; $pmax = 1000;
+#$pmin = 95; $pmax = 97;
+
+$ibest = 1; $ik = 0;
+
+$idir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_set_0001";
+
+for ($p = $pmin; $p <= $pmax; $p ++ ) {
+
+ $im = $p;
+ $stim = sprintf("%4.4i",$im);
+ $file1syn = "$idir/structure_model_${stim}.dat";
+ if (not -f $file1syn) {
+ die("Check if file1syn $file1syn exist or not\n");
+ }
+
+ $stirun = sprintf("%4.4i",$irun0+$im-1);
+ $chifile = "$odir/run_${stirun}/chi.dat";
+ if (not -f $chifile) {die("Check if chifile $chifile exist or not\n");}
+ open(IN,"$chifile");
+ $chi = <IN>;
+ if ($im==933) {$chi = 0} # target model
+ #$schi = sprintf("F(m) = %.1f",$chi);
+ $schi = sprintf("%.1f",$chi);
+ close(IN);
+
+ if ( $ibest==0 || ($ibest==1 && $chi <= 100) ) {
+ #print "$p, $schi\n";
+ #print CSH "echo $p, $schi\n";
+
+ if($ibest==1) {
+ $ik = $ik + 1;
+ $name = "gaussian_cov_model_ibest_${stim}";
+ #$name = "gaussian_cov_model_ibest_$ik";
+ } else {
+ $name = "gaussian_cov_model_${stim}";
+ }
+ $psfile = "$name.ps";
+ $jpgfile = "$name.jpg";
+
+ # phase velocity map
+ print CSH "psbasemap $B1syn $R $J1 -P -K -V $origin > $psfile\n"; # START
+ #print CSH "awk '{print \$1,\$2,\$6}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+ #print CSH "awk '{print \$2,\$1}' $idir/BOUNDARIES/oms_shelf |psxy $J1 $R $Wshelf -K -O -V >> $psfile\n";
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";
+
+ if (1==1) {
+ # faults
+ $faultfile = "/home/carltape/gmt/faults/jennings_more.xy";
+ if (not -f $faultfile) {
+ die("Check if faultfile $faultfile exist or not\n");
+ }
+ print CSH "psxy $faultfile $J1 $R -W0.75p -M -K -V -O >> $psfile\n";
+
+ # sources and receivers, plus chi label
+ if ($ibest==1) {
+ print CSH "pstext -N $J1 $R -K -O -V -Xa-0 -Ya-0.9 >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM $schi\nEOF\n";
+ print CSH "awk '{print \$1,\$2}' $efile | psxy $R $J1 $src -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2}' $rfile | psxy $R $J1 $rec -K -O -V >> $psfile\n";
+ }
+ }
+
+ # plot title and GMT header
+ $Utag = "-U/0/0.25/$plabel";
+ $shift = "-Xa-0.5 -Ya-0.5";
+
+ $stim = sprintf("%4.4i",$im);
+ print CSH "pstext -N $J1 $R -K -O -V -Xa-0 -Ya-0.4 >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM Model\nEOF\n";
+ print CSH "pstext -N $J1 $R -O -V -Xa-0 -Ya-0.65 >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM $stim\nEOF\n"; # FINISH
+
+ } # ibest = 1
+}
+
+#-----------------------------
+# print CSH "convert $psfile $jpgfile\n";
+
+ close (CSH);
+ system("csh -f $cshfile");
+ system("gv $psfile &")
+ #if($iopt <= 2) {system("xv $jpgfile &")}
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_all.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_mean.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_mean.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_mean.pl 2011-04-15 19:59:48 UTC (rev 18231)
@@ -0,0 +1,274 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_gauss_mean.pl
+# Carl Tape
+# 15-June-2010
+#
+# Figure to present at Albert Tarantola symposium
+#
+# ../plot_gauss_mean.pl
+#
+#==========================================================
+
+$imodel = 0;
+$im = 1;
+$irun0 = $im;
+
+# 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";
+
+$mfile_syn = "structure_syn.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");}
+
+$stimodel = sprintf("m%.1f",$imodel/2);
+$stirun0 = sprintf("%4.4i",$irun0);
+$stirun = sprintf("%4.4i",$irun0+$imodel);
+
+$dir0 = "$odir/run_${stirun0}";
+
+if(1==1) {
+ $dir1 = "$odir/run_${stirun}";
+} else {
+ $dir1 = sprintf("$dir0/READ_IN_CG/model_m%4.4i",$imodel);
+}
+
+$file1syn = "$dir1/${mfile_syn}";
+if (not -f $file1syn) { die("Check if $file1syn exist or not\n") }
+#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");
+
+#$edir = sprintf("$dir/event_%3.3i",$event);
+#print "\n $dir,$edir,$mfile_dat,$mfile_syn,$beta0,$per,$ifinite,$iopt \n";
+
+$plabel = "${plotdir}/plot_gauss_mean.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
+$src0 = "-W0.75p -Sa0.20 -G255/0/0";
+$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
+$grdfile = "temp.grd";
+
+#-------------------------
+# color
+
+# KEY: scaling for color
+$scale_color = 21.0;
+$colorbar = "seis";
+
+#-------------------------
+
+# write plotting scripts
+$wid = 2.5;
+$J1 = "-JM${wid}i"; # in lat-lon
+$origin = "-X1.5 -Y2.0";
+
+$B1dat = "-B1:.\" \":WeSN";
+$B1syn = "-B1:.\" \":wESN";
+
+$Dlen = 2.0;
+$Dx = $wid/2;
+$Dy = -0.35;
+
+$Dscale1 = "-D$Dx/$Dy/$Dlen/0.10h";
+
+$ptick = 0.05;
+#$Bscale1 = sprintf("-B%2.2e:\" Phase Velocity ( km s\@+-1\@+ )\": -E10p",$ptick);
+$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+
+#-------------------------
+
+$title = "Membrane Wave Speed";
+#$title = "Rayleigh Wave Phase Velocity";
+
+# 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 $file1syn`);
+#($tt) = sort{$b <=> $a} ($tt,abs($tmin),abs($tmax));
+$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 $smin $smax \n";
+
+# plot title
+$xtx = $xmin+0.5*($xmax-$xmin);
+$ztx = $zmin+1.10*($zmax-$zmin);
+
+ #===============================================
+ $cshfile = "plot_gauss_mean.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
+ $pmax = 0.1;
+ $dc = $pmax/10;
+ $T1 = "-T-$pmax/$pmax/$dc";
+ $cptfile1 = "color1.cpt";
+ print CSH "makecpt -C$colorbar -D $T1 > $cptfile1\n";
+
+ $pmax = 0.1/10;
+ $dc = $pmax/10;
+ $T1 = "-T-$pmax/$pmax/$dc";
+ $cptfile2 = "color2.cpt";
+ print CSH "makecpt -C$colorbar -D $T1 > $cptfile2\n";
+
+ $pmax = 0.0025;
+ $dc = $pmax/20;
+ $T1 = "-T0/$pmax/$dc";
+ $cptfile3 = "color3.cpt";
+ print CSH "makecpt -Chot -D -I $T1 > $cptfile3\n";
+
+ $Cfac = 10;
+ $pmax = 0.0025/$Cfac;
+ $dc = $pmax/20;
+ $T1 = "-T0/$pmax/$dc";
+ $cptfile4 = "color4.cpt";
+ print CSH "makecpt -Chot -D -I $T1 > $cptfile4\n";
+
+#=============================================
+
+ $origin = "-X1 -Y7.5";
+ $B1syn = "-B1:.\" \":WESN";
+
+#=============================================
+
+$idir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_set_0001";
+
+$name = "gaussian_mean";
+$psfile = "$name.ps";
+$jpgfile = "$name.jpg";
+
+ $B1syn = "-B1:.\" \":WesN";
+ $file1syn = "$idir/mean.dat";
+ if (not -f $file1syn) {die("Check if file1syn $file1syn exist or not\n");}
+ $olab = "-Xa1.3 -Ya-0.15";
+
+ # mean model
+ print CSH "psbasemap $B1syn $R $J1 -P -K -V $origin > $psfile\n"; # START
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+ #print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";
+ $faultfile = "/home/carltape/gmt/faults/jennings_more.xy";
+ if (not -f $faultfile) {die("Check if faultfile $faultfile exist or not\n");}
+ print CSH "psxy $faultfile $J1 $R -W0.75p -M -K -V -O >> $psfile\n";
+ print CSH "pstext -N $J1 $R -K -O -V $olab >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM Mean Model\nEOF\n";
+
+ $B1syn = "-B1:.\" \":Wesn";
+ $file1syn = "$idir/mean_est.dat";
+ if (not -f $file1syn) {die("Check if file1syn $file1syn exist or not\n");}
+ $dY0 = 3.2;
+ $shift = "-Y-$dY0";
+
+ # estimated mean model
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+ #print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";
+ print CSH "psxy $faultfile $J1 $R -W0.75p -M -K -V -O >> $psfile\n";
+ print CSH "pstext -N $J1 $R -K -O -V $olab >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM Estimated Mean Model\nEOF\n";
+
+ $B1syn = "-B1:.\" \":Wesn";
+
+ # estimated mean model -- zoomed
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile2 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+ #print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";
+ print CSH "psxy $faultfile $J1 $R -W0.75p -M -K -V -O >> $psfile\n";
+ print CSH "pstext -N $J1 $R -K -O -V $olab >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM Estimated Mean Model (x10)\nEOF\n";
+
+#------------------------------------------------------------------------
+# covariance matrices -- NOTE plotting j,i,Cij
+
+ $file1syn = "$idir/C_Cest.dat";
+ if (not -f $file1syn) {die("Check if file1syn $file1syn exist or not\n");}
+ $B1syn = "-Ba256f32:.\" \":WesN";
+ $J1 = "-JX${wid}i/-${wid}i";
+ $dY2 = 2*$dY0;
+ $shift = "-X3.5i -Y${dY2}i";
+ $xmin = 0; $xmax = 1025; $zmin = 0; $zmax = 1025;
+ $R = "-R$xmin/$xmax/$zmin/$zmax";
+ $interp = "-I2/2 -S2"; # key information
+
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
+ #print CSH "awk '{print \$1,\$2,\$3}' $file1syn | psxy -C$cptfile3 -Sc1p $R $J1 -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$2,\$1,\$3}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile3 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pstext -N $J1 $R -K -O -V $olab >>$psfile<<EOF\n $xmin $zmax $fsize0 0 $fontno CM Covariance Matrix\nEOF\n";
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V >> $psfile\n";
+
+ $shift = "-Y-$dY0";
+
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
+ #print CSH "awk '{print \$1,\$2,\$3}' $file1syn | psxy -C$cptfile3 -Sc1p $R $J1 -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$2,\$1,\$4}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile3 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pstext -N $J1 $R -K -O -V $olab >>$psfile<<EOF\n $xmin $zmax $fsize0 0 $fontno CM Maximum Likelihood Matrix\nEOF\n";
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V >> $psfile\n";
+
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
+ #print CSH "awk '{print \$1,\$2,\$3}' $file1syn | psxy -C$cptfile4 -Sc1p $R $J1 -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$2,\$1,\$4-\$3}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile4 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "pstext -N $J1 $R -K -O -V $olab >>$psfile<<EOF\n $xmin $zmax $fsize0 0 $fontno CM Residual Matrix (x$Cfac)\nEOF\n";
+ print CSH "psbasemap $B1syn $R $J1 -K -O -V >> $psfile\n";
+
+
+#-----------------------------
+# print CSH "convert $psfile $jpgfile\n";
+
+ print CSH "pstext -N $J1 $R -O -V -Xa-0 -Ya-0.65 >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno CM \nEOF\n"; # FINISH
+
+ close (CSH);
+ system("csh -f $cshfile");
+ system("gv $psfile &")
+ #if($iopt <= 2) {system("xv $jpgfile &")}
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_gauss_mean.pl
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2011-04-15 19:59:48 UTC (rev 18231)
@@ -25,6 +25,12 @@
# ../plot_surf_model.pl 4 700 0 1/0/5 0.1/0.05 stf_00001_1 0
# ../plot_surf_model.pl 4 800 0 0/1/5 0.1/0.05 stf_00001_1 0
#
+# ../plot_surf_model.pl 4 5000 0 1/0/5 0.1/0.05 stf_00001_1 0
+# ../plot_surf_model.pl 4 5000 2 1/0/5 0.1/0.05 stf_00001_1 0
+#
+# Set of 1000 models
+# ../plot_surf_model.pl 1 0001 0 1/0/1 0.1/0.05 stf_00001_1 0
+#
#==========================================================
if (@ARGV < 7) {die("Usage: plot_surf_model.pl xxx\n");}
@@ -51,11 +57,18 @@
$figdir = "${plotdir}/FIGURES";
if (not -e $figdir) {die("Check if figdir $figdir exist or not\n");}
+$stimodel = sprintf("m%.1f",$imodel/2);
$stirun0 = sprintf("%4.4i",$irun0);
$stirun = sprintf("%4.4i",$irun0+$imodel);
$dir0 = "$odir/run_${stirun0}";
-$dir1 = "$odir/run_${stirun}";
+
+if(1==1) {
+ $dir1 = "$odir/run_${stirun}";
+} else {
+ $dir1 = sprintf("$dir0/READ_IN_CG/model_m%4.4i",$imodel);
+}
+
$dir2 = sprintf("$dir1/event_%3.3i",$ievent);
$stf_file = "$dir2/${stf_file1}";
$rec_file = "$dir2/${rec_file1}";
@@ -182,6 +195,7 @@
$plot_title = "Model for synthetics";
$origin = "-X1.5 -Y5.0";
+ $B1syn = "-B1:.\" \":WESN";
# phase velocity map
print CSH "psbasemap $B1syn $R $J1 -P -K -V $origin > $psfile\n"; # START
@@ -192,9 +206,15 @@
#print CSH "awk '{print \$2,\$1}' $idir/BOUNDARIES/oms_shelf |psxy $J1 $R $Wshelf -K -O -V >> $psfile\n";
print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";
+ if (1==1) {
+ $faultfile = "/home/carltape/gmt/faults/jennings_more.xy";
+ if (not -f $faultfile) {die("Check if faultfile $faultfile exist or not\n");}
+ print CSH "psxy $faultfile $J1 $R -W0.75p -M -K -V -O >> $psfile\n";
+ }
+
# plot title and GMT header
$Utag = "-U/0/0.25/$plabel";
- $shift = "-Xa-$dX -Ya4";
+ $shift = "-Xa0 -Ya4";
print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
#=============================================
@@ -259,7 +279,7 @@
if($iopt == 4) {
# model for data (optional)
- $plot_title = "Target model, initial model, and source time function (run_${stirun})";
+ $plot_title = "Target model, model $stimodel, and source time function (run_${stirun})";
$shift = "-Y4.0 -X-0.9";
# phase velocity map
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 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl 2011-04-15 19:59:48 UTC (rev 18231)
@@ -22,7 +22,8 @@
$comp = 1;
# USER INPUT
-$ifig1 = 1;
+$ifig0 = 0;
+$ifig1 = 0;
$ifig2 = 1;
$itest = 0; # plot CG test models or not
@@ -168,34 +169,177 @@
$stmod = sprintf("m%2.2i",$imodel/2);
}
+ # target model
+ $filedat = "$dir0/${mfile_dat}";
+
# directories (current model and previos model)
$stirun = sprintf("%4.4i",$irun0+$imodel);
$dir1 = "$odir/run_${stirun}";
- $name = "model_all_${stirun}";
- $psfile = "$name.ps";
- $jpgfile = "$name.jpg";
+ #-----------------------------
+ # FIGURE 0
+ if ($ifig0 == 1) {
+
+ $name = "model_all_${stirun}_f0";
+ $psfile0 = "$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 mtar\": -E10p",$ptick);
+ $xlab = $wid/2;
+ $ylab = -0.02*$wid;
+ $olab = "-Xa$xlab -Ya$ylab";
+
+ print "\n----FIGURE 1--model $stmod----\n";
+
+ @files = ("$dir0/${mfile_syn}","$dir0/${mfile_syn}",
+ "$dir1/${mfile_syn}","$dir1/${mfile_syn}",
+ "$dir0/${mfile_dat}","$dir0/${mfile_dat}");
+ @stlabs = ("m00","mtar - m00",$stmod,"mtar - $stmod","mtar","mtar - mtar");
+
+ # only plot even models (m04 = iteration 8)
+ $kmaxt = 6;
+ $kmaxm = 6*(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;
+
+ $file1 = $files[$k];
+ if (not -f $file1) {
+ die("Check if $file1 exist or not\n");
+ }
+
+ # column-concatenate target model and current model
+ if ($k % 2 == 1) {
+ $file2 = "ftemp";
+ print CSH "paste $filedat $file1 > $file2\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 > $psfile0\n"; # START
+ } else {
+ print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile0\n";
+ }
+
+ if ($icolor == 1) {
+ #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile0\n";
+ if ($k % 2 == 0) {
+ print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
+ } else {
+ print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
+ }
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile0\n";
+ }
+ print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile0\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 >> $psfile0\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 >> $psfile0\n";
+
+ } else { # small circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile0\n";
+
+ }
+
+ # plot sources
+ if ($ievent_one) {
+ print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile0\n";
+ }
+ if ($ievent_all) {
+ print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile0\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 >>$psfile0<<EOF\n0 0 12 0 $fontno CM $stlab\nEOF\n";
+
+ # colorscale bars
+ if ($k==4) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile0\n";
+ }
+ if ($k==5) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile0\n";
+ }
+
+ # plot title and GMT header
+ if ($k==5) {
+ $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 >>$psfile0<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
+
+ if ($ijpg == 1) {
+ print CSH "convert $psfile0 $jpgfile\n";
+ }
+ if ($ipdf == 1) {
+ print CSH "ps2pdf $psfile0\n";
+ }
+
+ }
+ } # loop over subplots
+ }
+
+ #-----------------------------
+ # FIGURE 1
+
if ($ifig1 == 1) {
-# 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);
-$xlab = $wid/2;
-$ylab = -0.02*$wid;
-$olab = "-Xa$xlab -Ya$ylab";
+ $name = "model_all_${stirun}_f1";
+ $psfile1 = "$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);
+ $xlab = $wid/2;
+ $ylab = -0.02*$wid;
+ $olab = "-Xa$xlab -Ya$ylab";
+
print "\n----FIGURE 1--model $stmod----\n";
@files = ("$dir0/${mfile_syn}","$dir0/${mfile_syn}",
- "$dir1/${mfile_syn}","$dir1/${mfile_syn}",
- "$dir0/${mfile_dat}","$dir0/${mfile_dat}");
- @stlabs = ("m00","ln(m00/m00)",$stmod,"ln($stmod/m00)","mtar","ln(mtar/m00)");
+ "$dir1/${mfile_syn}","$dir1/${mfile_syn}",
+ "$dir0/${mfile_dat}","$dir0/${mfile_dat}");
+ @stlabs = ("m00","m00 - m00",$stmod,"$stmod - m00","mtar","mtar - m00");
# only plot even models (m04 = iteration 8)
$kmaxt = 6;
@@ -240,53 +384,53 @@
}
if ($k==0) {
- print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile\n"; # START
+ print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile1\n"; # START
} else {
- print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile\n";
+ print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile1\n";
}
if ($icolor == 1) {
- #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile1\n";
if ($k % 2 == 0) {
print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
} else {
print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
}
- print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile1\n";
}
- print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
+ print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile1\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";
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N $J1 $R -K -O -V $rec0 >> $psfile1\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";
+ print CSH "pstext $rec_file2 -N $J1 $R -K -O -V >> $psfile1\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";
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile1\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";
+ print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile1\n";
}
if ($ievent_all) {
- print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile1\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";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile1<<EOF\n0 0 12 0 $fontno CM $stlab\nEOF\n";
# colorscale bars
if ($k==4) {
- print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile\n";
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile1\n";
}
if ($k==5) {
- print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile\n";
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile1\n";
}
# plot title and GMT header
@@ -294,13 +438,13 @@
$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
+ print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile1<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
if ($ijpg == 1) {
- print CSH "convert $psfile $jpgfile\n";
+ print CSH "convert $psfile1 $jpgfile\n";
}
if ($ipdf == 1) {
- print CSH "ps2pdf $psfile\n";
+ print CSH "ps2pdf $psfile1\n";
}
}
@@ -327,7 +471,7 @@
print "\n--------FIGURE 2---models $stmodp, $stmodt, $stmod---\n";
- $name2 = "model_pert_${stirun}";
+ $name2 = "model_all_${stirun}_f2";
$psfile2 = "$name2.ps";
$jpgfile2 = "$name2.jpg";
@@ -351,9 +495,9 @@
@files1 = ("$dir1p/${mfile_syn}","$dir1p/${mfile_syn}","$dir0/${mfile_dat}",
"$dir1t/${mfile_syn}","$dir1t/${mfile_syn}","$dir0/${mfile_dat}",
"$dir1/${mfile_syn}","$dir1/${mfile_syn}","$dir0/${mfile_dat}");
- @stlabs = ($stmodp,"ln($stmodp/$stmodp)","ln(mtar/$stmodp)",
- $stmodt,"ln($stmodt/$stmodp)","ln(mtar/$stmodt)",
- $stmod,"ln($stmod/$stmodp)","ln(mtar/$stmod)");
+ @stlabs = ($stmodp,"$stmodp - $stmodp","mtar - $stmodp",
+ $stmodt,"$stmodt - $stmodp","mtar - $stmodt",
+ $stmod,"$stmod - $stmodp","mtar - $stmod");
$kmax = 9;
@@ -476,7 +620,8 @@
close (CSH);
system("csh -f $cshfile");
-if($ifig1==1) {system("gv $psfile &")}
+if($ifig0==1) {system("gv $psfile0 &")}
+if($ifig1==1) {system("gv $psfile1 &")}
if($ifig2==1) {system("gv $psfile2 &")}
#=================================================================
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 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2011-04-15 19:59:48 UTC (rev 18231)
@@ -23,7 +23,7 @@
clear
% add path to additional matlab scripts
-path(path,[pwd '/matlab_scripts']);
+path([pwd '/matlab_scripts'],path);
colors;
npts = 100;
@@ -43,6 +43,8 @@
%parms = [500 0 3]; % joint inversion: 25 sources, Nfac=3
%parms = [600 0 3]; % structure inversion: 25 sources, Nfac=2
+parms = [6000 0 3];
+
%---------------------------------------------------------
irun0 = parms(1); % irun for m0
@@ -55,14 +57,16 @@
%---------------------------------------------------------
-odir = 'OUTPUT/';
+odir0 = 'OUTPUT/';
+odir = [dir0 odir0];
+%odir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_OUTPUT_OLD_2/'; % CHT
strun0 = ['irun0 = ' num2str(irun0) ];
% determine the number of iterations
k = 0;
while k < 100
irun = irun0 + k;
- chifile = [dir0 odir 'run_' sprintf(stfm,irun) '/chi.dat'];
+ chifile = [odir 'run_' sprintf(stfm,irun) '/chi.dat'];
iexist = exist(chifile);
if iexist == 2
chi = load(chifile);
@@ -95,7 +99,7 @@
% chits = zeros(niter,1);
%
% % load initial chi value
-% stf = [dir0 odir 'run_' sprintf(stfm,irun0) '/'];
+% stf = [odir 'run_' sprintf(stfm,irun0) '/'];
% ww = 'summed_chi_all';
% load([stf ww '.dat']); chi = eval(ww); chis(1) = chi;
%
@@ -103,8 +107,8 @@
% for ii=1:niter
% irun = irun_vec(ii);
%
-% stft = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
-% stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
+% stft = [odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+% stf = [odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
% chit = load([stft 'summed_chi_all.dat']);
% chi = load([stf 'summed_chi_all.dat']);
% chis(ii+1) = chi;
@@ -135,7 +139,7 @@
%if or(ifinal_test==1, ii < niter)
if icubic == 1
ww = 'cubic_poly';
- stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
@@ -143,7 +147,7 @@
g1 = temp(5); g2 = temp(6);
else
ww = 'quad_poly';
- stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
@@ -387,8 +391,8 @@
plot(x_var,var_red1,'b.-','markersize',msize,'linewidth',lsize);
axis([xlims 0 1]); grid on; set(gca,'xtick',xticks);
xlabel(stx);
-ylabel(' variance reduction');
-title(' variance reduction between successive models');
+ylabel(' Logarithmic Variance Reduction');
+title(' Variance reduction between successive models');
% subplot(nr,nc,4); plot(its(2:end),var_red2,'.-');
% xlim(xlims); grid on; set(gca,'xtick',its);
@@ -405,7 +409,7 @@
plot(its_smooth,converge_fit,'r--',converge_its,converge_pts,'.','markersize',msize,'linewidth',lsize);
axis([xlims 0 0.7]); grid on; set(gca,'xtick',xticks);
xlabel(stx);
-ylabel(' order of convergence rate, | \Delta log10(S) / \Delta k |');
+ylabel(' Order of convergence rate, | \Delta log10(S) / \Delta k |');
fontsize(9), orient tall, wysiwyg
@@ -416,7 +420,7 @@
% specify number of iterations to plot; load initial chi
niter = 1;
stirun = [' irun0 = ' num2str(irun0) '; niter = ' num2str(niter) ];
-stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun0)) '/'];
+stf = [odir 'run_' num2str(sprintf(stfm,irun0)) '/'];
ww = 'chi';
load([stf ww '.dat']); chi = eval(ww); chis(1) = chi;
@@ -428,14 +432,14 @@
if icubic == 1
ww = 'cubic_poly';
- stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
y1 = temp(3); y2 = temp(4);
g1 = temp(5); g2 = temp(6);
else
ww = 'quad_poly';
- stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
y1 = temp(3); y2 = temp(4);
@@ -554,7 +558,7 @@
% load actual chi-value computed from the next run
% wave2d: xx1,xx2,yy1,yy2,g1,g2,Pa,Pb,Pc,Pd,xmin
ww = 'chi';
- stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
+ stf = [odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
load([stf ww '.dat']);
chi = eval(ww);
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 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_run.m 2011-04-15 19:59:48 UTC (rev 18231)
@@ -16,7 +16,7 @@
clear
% add path to additional matlab scripts
-path(path,[pwd '/matlab_scripts']);
+path([pwd '/matlab_scripts'],path);
colors;
stfm = '%4.4i';
@@ -33,24 +33,33 @@
if ~exist(dirbase,'dir'), error(['dirbase ' dirbase ' does not exist']); end
% directory with Matlab parameterization for full covariance matrix
-dirrand = [dirbase '/SEM2D_iterate_INPUT/random_fields/model_001/'];
+im = 1; % KEY: index for model (im=1 for 50km, im=2 for 25km)
+stim = sprintf('%3.3i',im);
+dirrand = [dirbase '/SEM2D_iterate_INPUT/random_fields/model_' stim '/'];
if ~exist(dirrand,'dir'), error(['dirrand ' dirrand ' does not exist']); end
% KEY
-parms = [1000 0];
-%parms = [1500 4];
+%parms = [4000 31];
+%parms = [6000 31];
+parms = [6000 8];
-% 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;
+% icg = 1: emulate wave2d.f90 algorithm (diagonal Cm)
+% icg = 2: emulate wave2d.f90 algorithm, but use 32 x 32 cells (diagonal Cm)
+% icg = 3: CG algorithm with full Cm and 32 x 32 cells
+% icg = 4: source subspace algorithm with full Cm and 32 x 32 cells
+icg = 4;
iwrite = 1;
iplotmod = 1;
-iplotker = 1;
+iplotker = 0;
sticg = sprintf('%2.2i',icg);
+nsrc = 25;
+
%---------------------------------------------------------
+% whether to use smoothed event kernels
+if any(icg == [3 4]), ismooth = 0; else ismooth = 1; end
+
irun0 = parms(1);
istep = parms(2);
imaketest = ~mod(istep,2);
@@ -87,16 +96,14 @@
if ~exist(dirm0,'dir'), error(['dirm0 ' dirm0 ' does not exist']); end
% assign the input and output directories
+idir1 = dirm0;
+idir2 = dirmp;
if imaketest==1
imo = imt;
- idir1 = dirm0;
- idir2 = dirmp;
- odir = dirmt;
-
+ odir = dirmt;
+ if icg==4, imo = imk; odir = dirmk; end
else
imo = imk;
- idir1 = dirm0;
- idir2 = dirmp; % directory for previous gradient
odir = dirmk;
if ~exist(dirmt,'dir'), error(['dirmt ' dirmt ' does not exist']); end
end
@@ -187,7 +194,7 @@
%[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');
-% prior, initiral, and target models
+% prior, initial, and target models
[mprior, m0_initial, mtarget] = textread([dir0 'prior_initial_target_models.dat'],'%f%f%f');
% CURRENT model vector
@@ -196,32 +203,6 @@
%[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) '/'];
- 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
- Zp = griddata(xg,yg,Kbeta,Xp,Yp,'nearest');
- figure; imagesc(Zp); set(gca,'ydir','normal');
- caxis([-1 1]*max(abs(Kbeta))*0.8);
- end
-end
-summed_ker = sum(kernels_all,2);
-
-% source gradient
-source_gradient = load([idir1 'source_gradient.dat']);
-%source_gradient = zeros(nmod_src,1);
-
% data covariance
% nmeas_run = nevent_run * nrec * NCOMP
% double precision, parameter :: SIGMA_DT = 0.10
@@ -251,8 +232,41 @@
data_norm0 = load([idir1 'data_norm.dat']);
model_norm0 = load([idir1 'model_norm.dat']);
chi0 = load([idir1 'chi.dat'])
-data_norm0, data_norm, sum(chi_data_vec)
+data_norm0, data_norm, sum(chi_data_vec), DTvec' * diag(1./covd_diag) * DTvec
+% vector of number of measurement windows per event
+Ns_vec = nrec * ones(nsrc,1); % KEY: same measurements per event
+Nsi_vec = zeros(nmeas_run,1);
+cNs = cumsum(Ns_vec);
+Ns_inds = [ [1 ; cNs(1:end-1)+1] cNs ];
+
+% KEY: compute projected data vector
+[dnorm2, indmat] = wave2d_dsub(DTvec,covd_diag,Ns_vec);
+whos chi_data_vec dnorm2
+data_norm0, data_norm, sum(chi_data_vec), sum(dnorm2)
+
+% compute variance reduction for EVERY SINGLE WINDOW
+if and(imaketest==1, istep >= 2)
+ % load measurements for PREVIOUS model
+ meas_allp = load([idir2 'measure_vec.dat']);
+ DTvecp = meas_allp(:,1);
+
+ [VR_i, VR_s, VR_tot] = wave2d_VR(DTvecp,DTvec,sigma_DT,covd_diag,Ns_inds);
+end
+
+% KEY: stopping criterion
+if and(imaketest==1, istep >= 2)
+ % previous chi value
+ chip = load([dirmp 'chi.dat']);
+ disp('----------------');
+ chip, chi0, log( chip / chi0 ), VAR_RED_MIN
+ if log( chip / chi0 ) < VAR_RED_MIN
+ error('VR stopping criteria has been met');
+ end
+end
+
+%----------------------------------
+
% sigma values
% write(19,'(2e20.10)') sigma_beta, m_scale_str(1)
% write(19,'(2e20.10)') sigma_ts, m_scale_src(1)
@@ -268,6 +282,33 @@
[alpha0, beta0, rho0, bulk0, kappa0, mu0] = ...
textread([idir1 'reference_values.dat'],'%f%f%f%f%f%f');
+% 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_smooth']);
+ if ismooth == 0
+ Kbeta = kernel(:,4); cfac = 0.2;
+ else
+ Kbeta = kernel(:,3); cfac = 0.8;
+ end
+ kernels_all(:,ii) = Kbeta;
+
+ if and(iplotker==1,icg==1)
+ Zp = griddata(xg,yg,Kbeta,Xp,Yp,'nearest');
+ figure; imagesc(Zp); set(gca,'ydir','normal');
+ axis equal, axis tight
+ colormap(seis); caxis([-1 1]*max(abs(Kbeta))*cfac);
+ title(sprintf('Event kernel %i/%i',ii,nevent_run));
+ end
+end
+summed_ker = sum(kernels_all,2);
+
+% source gradient
+source_gradient = load([idir1 'source_gradient.dat']);
+%source_gradient = zeros(nmod_src,1);
+
%========================================================================
% ONE STEP OF A CONJUGATE GRADIENT ALGORITHM
@@ -291,6 +332,7 @@
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);
@@ -423,12 +465,6 @@
% test model
mt = m0 + lam_t_val*pk;
- mt_vec = zeros(nmod,1);
- %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);
-
else
% load chi for test model
chi_t_val = load(chitfile);
@@ -438,7 +474,7 @@
xx2 = lam_t_val
yy1 = chi_k_val
yy2 = chi_t_val
- g1 = sum(gk .* pk)
+ g1 = sum(gk .* pk);
%g1 = sum(g0 .* pk);
% coefficients of the quadratic polynomial (ax^2 + bx + c)
@@ -461,17 +497,12 @@
lam_0_val = xmin;
mk = m0 + lam_0_val * pk;
- 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);
-
end
%---------------------------------------------------------
-% icg = 1 ==> text
-elseif or(icg == 2, icg==3)
+
+elseif any(icg == [2 3 4])
+
% load Matlab parameterization of the fields (gaussian_2D.m)
load([dirrand 'matlab_vars']);
@@ -484,7 +515,7 @@
% compute the GLL gridpoint (local) closest to each nx x ny cell;
% this will give a quick mapping from kernels to cells
- [iGLL, dGLL] = wave2d_gll2cell(xg,yg,x,y);
+ %[iGLL, dGLL] = wave2d_gll2cell(xg,yg,x,y);
% indexing for cell-based model
m_indsi = m_inds - nlocal + ncell; m_indsi(1,1) = 1;
@@ -501,7 +532,7 @@
Avec = sqrt( dAvec/Atot );
n*dA - Atot, sum(Avec.^2) % check
- if 1==1
+ if 0==1
% plot full covariance matrix
figure; imagesc(C);
title('prior model covariance matrix');
@@ -539,7 +570,7 @@
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
+ % KEY: weighting for full covariance matrix
Cmod = diag(sqrt(cov_weight_vec)) * Cfull * diag(sqrt(cov_weight_vec));
else
@@ -574,42 +605,50 @@
if icg==2
Crun = diag(diag(Cmod));
- else
+ else % icg = 3,4
Crun = Cmod;
end
%--------------------------------------
- % model
+ % models
- 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);
- 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);
+ m0i = wave2d_m_gll2cell(m0,xg,yg,Xp,Yp,nx,ny,nmod_src,1);
+ mpriori = wave2d_m_gll2cell(mprior,xg,yg,Xp,Yp,nx,ny,nmod_src,1);
+ mtargeti = wave2d_m_gll2cell(mtarget,xg,yg,Xp,Yp,nx,ny,nmod_src,1);
%--------------------------------------
% gradient
% interpolate kernels onto cells
kernels_alli = zeros(ncell,nevent_run);
- for ii=1:nevent_run
- Kbeta = kernels_all(iGLL,ii); % crude interpolation
- kernels_alli(:,ii) = Kbeta;
+ for ii=1:nevent_run
+ % kernels_all is computed above
+ Kbeta0 = kernels_all(:,ii);
+ %Kbeta = Kbeta0(iGLL); % crude interpolation
+ Kbeta = wave2d_gll2cell(xg,yg,Kbeta0,Xp,Yp);
+ kernels_alli(:,ii) = Kbeta(:);
if iplotker==1
- % 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);
+ % plot ghat and g = Cm*ghat
+ figure; nr=2; nc=1;
+ Zp = Kbeta;
+ %Zp = reshape( Kbeta, ny, nx); % gradient
+ subplot(nr,nc,1); imagesc(Zp); set(gca,'ydir','normal'); axis equal, axis tight;
+ caxis([-1 1]*max(abs(Kbeta(:)))*0.8); colormap(seis); colorbar;
+ title(sprintf('kernel %i out of %i',ii,nevent_run));
- % 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);
+ Zp = reshape( C * Kbeta(:), ny, nx); % steepest ascent vector
+ subplot(nr,nc,2); imagesc(Zp); set(gca,'ydir','normal'); axis equal, axis tight;
+ caxis([-1 1]*max(abs(C * Kbeta(:)))*0.8); colormap(seis); colorbar;
+ title(sprintf('kernel %i out of %i',ii,nevent_run));
+ orient tall, wysiwyg;
+
+ if 0==1
+ %print(gcf,'-depsc',sprintf('kernel_%2.2i',ii) );
+ %print(gcf,'-dpsc',sprintf('kernel_%2.2i',ii) );
+ print(gcf,'-dpdf',sprintf('kernel_%2.2i',ii) );
+ pause(1);
+ end
end
end
summed_keri = sum(kernels_alli,2);
@@ -655,34 +694,217 @@
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
- % compute new model vector (mt or mk)
- gfile = [idir2 gfilemp];
- mu_val = chi_data_stop;
- gk = gradienti_tot;
- [M,pk] = wave2d_cg(m0i,gk,Crun,chi_k_val,mu_val,istep,imaketest,gfile,chitfile);
+ if icg ~=4
+ % emulate CG algorithm in wave2d.f90 -- THIS ASSUMES A FULL COVARIANCE MATRIX
+
+ % compute new model vector (mt or mk)
+ gfile = [idir2 gfilemp];
+ mu_val = chi_data_stop;
+ gk = gradienti_tot;
+ [M,pk] = wave2d_cg(m0i,gk,Crun,chi_k_val,mu_val,istep,imaketest,gfile,chitfile);
+
+ if imaketest==1
+ mti = M;
+ else
+ mki = M;
+ end
- % 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);
- mti_vec(1:ncell) = beta0 * exp( mti(1:ncell) );
- mti_vec(ncell+1 : nmodi) = mti(ncell+1 : nmodi);
+ else
+ % blocks from Cm
+ Cstr = Cmod(indBi,indBi);
+ csrc = diag(Cmod(indTXYi,indTXYi)); % diagonal
- else
- mki = M;
- mki_vec = zeros(nmodi,1);
- mki_vec(1:ncell) = beta0 * exp( mki(1:ncell) );
- mki_vec(ncell+1 : nmodi) = mki(ncell+1 : nmodi);
+ % gradient for structure
+ Gstr = zeros(nevent_run,ncell);
+ wvec = sqrt( dnorm2 );
+ Gstr = -kernels_alli' .* repmat(dAvec',nevent_run,1) .* repmat(1./wvec, 1, ncell);
+
+ % gradient for source
+ Gsrc = zeros(nevent_run,nmod_src);
+ % QUESTION: SHOULD THIS BE THE FULL GRADIENT OR JUST THE DATA TERM?
+ %Gsrc = [diag(gradienti_tot(indTi)) diag(gradienti_tot(indXi)) diag(gradienti_tot(indYi))];
+
+ % partial derivatives matrix (nsrc x nmodi)
+ G = [Gstr Gsrc];
+
+ % modified data vector
+ dmod = sqrt(dnorm2) + G*(m0i-mpriori);
+
+ % compute subspace Hessian
+ [Hstr,Hsrc] = wave2d_compute_hessian(Gstr,Gsrc,Cstr,csrc);
+ H = Hstr + Hsrc + eye(nsrc,nsrc);
+ wave2d_plot_hessian(Hstr,Hsrc,H);
+
+ % write out files for Hessian
+ %if iwrite==1
+ % nprint = 20;
+ % wave2d_write_hessian(H,eids,dir2,nprint);
+ %end
+
+ % examining H\dnorm2 operation with TSVD
+ %mu_all = wave2d_tsvd_hessian(H,dnorm2);
+ %mu_all = wave2d_tsvd_hessian(H,sqrt(dnorm2));
+ mu_all = wave2d_tsvd_hessian(H,dmod);
+ if 0==1
+ lambda = 1;
+ ndiff = zeros(nsrc,1);
+ for ii=1:nsrc
+ %for ii=1:5
+ mu = mu_all(:,ii);
+ mplot = lambda * Cstr * Gstr' * mu;
+ pmax = max(abs(mplot));
+ pmax = 0.1;
+
+ ndiff(ii) = norm( mplot - (mtargeti(indBi) - m0i(indBi)) );
+
+ if 1==1
+ Zp = griddata(x,y,mplot,Xp,Yp,'nearest');
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1); imagesc(Zp); set(gca,'ydir','normal');
+ title(sprintf('dm -- TSVD truncated at index %i / %i, norm = %.2e',ii,nsrc,ndiff(ii)));
+ caxis([-1 1]*pmax); colormap(seis); colorbar;
+ axis equal, axis tight
+
+ subplot(nr,nc,2); hold on;
+ if ii > 1, plot(mu_all(:,ii-1),'ro--','markersize',8,'linewidth',1); end
+ if ii < nsrc, plot(mu_all(:,ii+1),'ko--','markersize',8,'linewidth',1); end
+ plot(mu,'b.-','markersize',20,'linewidth',2);
+ axis([0 nsrc+1 max(abs(mu))*[-1 1]]);
+ grid on; xlabel('event index'); ylabel('elements of mu vector');
+ title('\mu_{-1} (RED) --> \mu_0 (BLUE) --> \mu_{+1} (BLACK)');
+ orient tall, wysiwyg
+
+ print(gcf,'-dpdf',sprintf('%s_m%s_dm_tsvd_%2.2i',stirun0,stim0,ii) );
+ end
+ end
+
+ % difference between target model and new model
+ figure; plot(ndiff,'.-','markersize',20,'linewidth',2); grid on;
+ xlabel('TSVD truncation index'); ylabel('norm( dm0 - (mtar-m0) )');
+ print(gcf,'-dpdf',sprintf('%s_m%s_dm_ndiff',stirun0,stim0) );
+ end
+
+ % initialize to no update
+ dm_str = zeros(nmod_str,1);
+ dm_src = zeros(nmod_src,1);
+ dm_str1 = zeros(nmod_str,1); dm_str2 = zeros(nmod_str,1);
+ dm_src1 = zeros(nmod_src,1); dm_src2 = zeros(nmod_src,1);
+
+ lambda = 1;
+ itest = 2;
+ switch itest
+ case 1
+ % KEY: compute model update
+ mu = H \ sqrt(dnorm2); % why not dmod?
+ %mu = mu_all(:,15);
+
+ % compute the STRUCTURE model update
+ if INV_STRUCT_BETA == 1
+ dm_str = Cstr * Gstr' * mu;
+ end
+ % compute the SOURCE model update
+ if and(INV_SOURCE_T==1, INV_SOURCE_X==1)
+ dm_src = csrc .* (Gsrc' * mu);
+ end
+ dm = [dm_str ; dm_src];
+ mki = m0i + dm;
+
+ case 2 % actual iterative algorithm
+ mu = H \ dmod;
+ if INV_STRUCT_BETA == 1
+ dm_str = Cstr * Gstr' * mu;
+ end
+ if and(INV_SOURCE_T==1, INV_SOURCE_X==1)
+ dm_src = csrc .* (Gsrc' * mu);
+ end
+ dm = [dm_str ; dm_src];
+ mki = mpriori + dm;
+
+ case 3 % actual iterative algorithm
+ if INV_STRUCT_BETA == 1
+ mu1 = H \ sqrt(dnorm2);
+ mu2 = H \ (Gstr*(m0i(indBi)-mpriori(indBi)));
+ dm_str1 = Cstr * Gstr' * mu1;
+ dm_str2 = Cstr * Gstr' * mu2;
+ end
+ dm1 = [dm_str1 ; dm_src1];
+ dm2 = [dm_str2 ; dm_src2];
+ dm = dm1 + dm2;
+ mki = mpriori + dm; % mpriori or m0i
+ figure; plot( dm1, dm2 ,'.')
+
+ case 4
+ if INV_STRUCT_BETA == 1
+ dm_str1 = (Gstr'*Gstr + Cmodinv(indBi,indBi) ) \ (Gstr'*sqrt(dnorm2));
+ dm_str2 = -(diag(diag(Cstr))*Gstr'*Gstr + eye(ncell,ncell) ) \ (m0i(indBi)-mpriori(indBi));
+ end
+ dm1 = [dm_str1 ; dm_src1];
+ dm2 = [dm_str2 ; dm_src2];
+ dm = dm1 + dm2;
+ mki = m0i + dm;
+ figure; plot( dm1, dm2 ,'.')
+ end
+
+ % linear combination of kernels
+ figure; plot(mu,'.-'); xlabel('event index');
+ title('mu vector for source subspace approach');
+
+ figure; nr=3; nc=3; cmax = 0.1;
+ for ii=1:6
+ switch ii
+ case 1, mplot = mpriori(1:ncell); smt = 'mprior';
+ case 2, mplot = mtargeti(1:ncell); smt = 'mtarget';
+ case 3, mplot = m0i(1:ncell); smt = 'm0';
+ case 4, mplot = mki(1:ncell); smt = 'mnew';
+ if itest==3, smt = 'mnew = mpriori + dm'; else smt = 'mnew = m0 + dm'; end
+ case 5, mplot = mki(1:ncell) - m0i(1:ncell); smt = 'mnew - m0';
+ case 6, mplot = mtargeti(1:ncell) - m0i(1:ncell); smt = 'mtarget - m0';
+ end
+ Zp = griddata(x,y,mplot,Xp,Yp,'nearest');
+ subplot(nr,nc,ii);
+ imagesc(Zp); set(gca,'ydir','normal'); title(smt)
+ caxis([-1 1]*cmax); colormap(seis); axis equal, axis tight
+ end
+
+ if or(itest==3, itest==4)
+ for ii=7:9
+ switch ii
+ case 7, mplot = dm(1:ncell); smt = 'dm';
+ case 8, mplot = dm1(1:ncell); smt = 'dm(term 1)';
+ case 9, mplot = dm2(1:ncell); smt = 'dm(term 2)';
+ end
+ Zp = griddata(x,y,mplot,Xp,Yp,'nearest');
+ subplot(nr,nc,ii);
+ imagesc(Zp); set(gca,'ydir','normal'); title(smt)
+ caxis([-1 1]*cmax); colormap(seis); axis equal, axis tight
+ end
+ orient tall, wysiwyg
+ else
+ % uniform sum to make 'misfit kernel' -- this is the pattern for the CG algorithm
+ mplot = Cstr * Gstr' * ones(nsrc,1);
+ pmax = max(abs(mplot));
+ Zp = griddata(x,y,mplot,Xp,Yp,'nearest');
+ subplot(nr,nc,7); imagesc(Zp); set(gca,'ydir','normal'); title('dm (CG)');
+ caxis([-1 1]*pmax); colormap(seis); axis equal, axis tight
+ orient tall, wysiwyg
+ end
+
+% disp(' norm, min, max of structure update:');
+% norm_dm_str = dm_str' * diag(1./diag(Cstr)) * dm_str
+% min_dm_str = min( dm_str )
+% max_dm_str = max( dm_str )
+% disp(' norm, min, max of source update:');
+% norm_dm_src = sum( dm_src .* dm_src ./ csrc )
+% min_dm_src = min( dm_src )
+% max_dm_src = max( dm_src )
end
+
else
error('invalid icg value');
@@ -691,40 +913,85 @@
%========================================================================
% WRITE NEW MODEL TO FILE
-% write current gradient and search direction to file
-if iwrite==1
- wave2d_write_grad([idir1 gfilem0],gk,pk); % note idir1
-end
+if icg ~= 4 % CG method
+ % write current gradient and search direction to file
+ if iwrite==1
+ wave2d_write_grad([idir1 gfilem0],gk,pk); % note idir1
+ end
-if icg==1
- if imaketest==1, M_vec = mt_vec; M = mt; else M_vec = mk_vec; M = mk; end
- M0 = m0; M0_vec = m0_vec;
- x0 = xg; y0 = yg;
-else
- if imaketest==1, M_vec = mti_vec; M = mti; else M_vec = mki_vec; M = mki; end
- M0 = m0i; M0_vec = m0i_vec;
+ if icg==1
+ if imaketest==1, M = mt; else M = mk; end
+ M0 = m0; Mtar = mtarget;
+ x0 = xg; y0 = yg;
+ else
+ if imaketest==1, M = mti; else M = mki; end
+ M0 = m0i; Mtar = mtargeti;
+ x0 = x; y0 = y;
+ m_inds = m_indsi;
+ end
+
+else % source subspace method
+ M = mki;
+ M0 = m0i; Mtar = mtargeti;
+ x0 = x; y0 = y;
m_inds = m_indsi;
- x0 = x; y0 = y;
end
-% initial model
+% convert to physical parameters (only differes for structure parameters)
+% --> current, new, and target models
+M0_vec = wave2d_m2mvec(M0,m_inds,beta0);
+M_vec = wave2d_m2mvec(M,m_inds,beta0);
+Mtar_vec = wave2d_m2mvec(Mtar,m_inds,beta0);
+
+% split into components
[m0_B, m0_ts, m0_xs, m0_ys] = wave2d_splitm(M0,m_inds);
+[m_B, m_ts, m_xs, m_ys] = wave2d_splitm(M,m_inds);
+[mtar_B, mtar_ts, mtar_xs, mtar_ys] = wave2d_splitm(Mtar,m_inds);
+
[m0_vec_B, m0_vec_ts, m0_vec_xs, m0_vec_ys] = wave2d_splitm(M0_vec,m_inds);
-% new model (test or update)
-[m_B, m_ts, m_xs, m_ys] = wave2d_splitm(M,m_inds);
[m_vec_B, m_vec_ts, m_vec_xs, m_vec_ys] = wave2d_splitm(M_vec,m_inds);
+[mtar_vec_B, mtar_vec_ts, mtar_vec_xs, mtar_vec_ys] = wave2d_splitm(Mtar_vec,m_inds);
if iplotmod==1
% change w.r.t. previous model
- Zp = griddata(x0,y0,log(m_vec_B ./ m0_vec_B),Xp,Yp,'nearest');
- figure; imagesc(Zp); set(gca,'ydir','normal');
- colormap(seis); caxis([-1 1]*0.1); colorbar;
+ %Zp = griddata(x0,y0,log(m_vec_B ./ m0_vec_B),Xp,Yp,'nearest');
+ %figure; imagesc(Zp); set(gca,'ydir','normal');
+ %colormap(seis); caxis([-1 1]*0.1); colorbar;
+
+ figure; nr=3; nc=2; cmax = 0.1;
+ if imaketest==1
+ tlabs = {'ln(m0/beta0)','mtar - m0','ln(mt/beta0)','mtar - mt','ln(mtar/beta0)','mt - m0'};
+ else
+ tlabs = {'ln(m0/beta0)','mtar - m0','ln(mk/beta0)','mtar - mk','ln(mtar/beta0)','mk - m0'};
+ end
+ for kk = 1:6
+ switch kk
+ case 1, mplot = log(m0_vec_B/beta0);
+ case 2, mplot = log(mtar_vec_B./m0_vec_B);
+ case 3, mplot = log(m_vec_B/beta0);
+ case 4, mplot = log(mtar_vec_B./m_vec_B);
+ case 5, mplot = log(mtar_vec_B/beta0);
+ case 6, mplot = log(m_vec_B./m0_vec_B);
+ end
+
+ Zp = griddata(x0,y0,mplot,Xp,Yp,'nearest');
+ %mnorm = sum( mplot.^2 ./ cov_beta' );
+
+ subplot(nr,nc,kk); hold on;
+ imagesc(Zp); set(gca,'ydir','normal');
+ caxis([-1 1]*cmax); colormap(seis); axis equal, axis tight
+ title(tlabs{kk});
+ %plot(rlon,rlat,'k.','markersize',16);
+ %plot(slon,slat,'p','markersize',18,'markeredgecolor','k','linewidth',2);
+ %title([tlabs{kk} sprintf(', norm = %.2f',mnorm)]);
+ end
+ orient tall, wysiwyg
end
if iwrite==1
% interpolate onto SEM mesh
if icg ~= 1
- m_B = wave2d_cell2gll(xg,yg,x,y,m_B,nx,ny);
+ m_B = wave2d_cell2gll(xg,yg,Xp,Yp,reshape(m_B,ny,nx));
end
% write structure model -- ONLY MU CHANGES (only m_B appears)
@@ -763,11 +1030,4 @@
end
-%---------------------------------------------------------
-% compute CG update using a diagonal covariance matrix with interpolated cells for model parameters
-
-
-
-%/home/carltape/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_001/
-
%=============================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m 2011-04-15 19:59:48 UTC (rev 18231)
@@ -1,20 +1,24 @@
%
% wave2d_subspace.m
-% CARL TAPE, 26-Jan-2010
+% Carl Tape, 26-Jan-2010
%
% This program implements the subspace method inversion notes developed by
-% Jeroen Tromp, Malcolm Sambridge, and Carl Tape.
+% Malcolm Sambridge, Jeroen Tromp, and Carl Tape.
+%
+% This assumes a diagonal model covariance matrix and requires
+% regularization -- several choices are explored here.
%
% calls xxx
% called by xxx
%
-clc, clear, close all
+clc, clear
+%close all
format short, format compact
%warning off
% add path to additional matlab scripts
-path(path,[pwd '/matlab_scripts']);
+path([pwd '/matlab_scripts'],path);
colors;
@@ -26,122 +30,14 @@
dir_run = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_OUTPUT/';
-TSVD = 0; % truncated singular value decomposition
-
+TSVD = 1; % truncated singular value decomposition
inew = 0; % subspace AND parameter class inversions
%----------------------------------------------
-if 0==1
- temp = load('/net/denali/scratch1/carltape/OUTPUT/run_7000/structure_dat.dat');
- %temp = load('/net/denali/scratch1/carltape/OUTPUT/run_7000/READ_IN/structure_dat.dat');
- lon = temp(:,1); lat = temp(:,2); beta = temp(:,7);
-
- [X,Y,Z] = griddataXB(lon,lat,beta,100,'nearest');
- figure; cmax = 0.1; hold on;
- pcolor(X,Y,Z); shading interp;
- caxis([-1 1]*cmax); colormap(seis);
-end
-
-% check the relative dimensions of all the matrices we need
-if 0==1
- S = 5; % number of sources
- Cstr = 1; % number of structure parameter classes
- Csrc = 2; % number of source parameter classes
- C = Cstr + Csrc;
- Mstr1 = 20; % number of structure parameters
- Msrc1 = 1; % number of source parameters in class 1 (origin time)
- Msrc2 = 2; % number of source parameters in class 2 (xy, ys)
- Mstr = Mstr1;
- Msrc = S*(Msrc1 + Msrc2);
- M = Mstr + Msrc; % number of rows in P
- N = C*S; % number of columns in P
-
- P = zeros(M,N);
- P(1:Mstr1,1:S) = rand(Mstr1,S);
- for s = 1:S
- P(Mstr1+s,S+s) = rand;
- P(Mstr1+S+s,2*S+s) = rand;
- P(Mstr1+2*S+s,2*S+s) = rand;
- end
-
- % fill G
- G = zeros(S,M);
- Gstr = rand(S,Mstr);
- Gsrc = zeros(S,Msrc);
- if 1==1
- Gsrc = repmat(diag(rand(S,1)),1,Msrc1+Msrc2);
- else
- itemp = (Msrc1 + Msrc2)*[1:S]';
- indmat2 = [ [1 ; 1+itemp(1:end-1)] itemp ];
- for isrc = 1:S
- inds = indmat2(isrc,1) : indmat2(isrc,2);
- %Gsrc(isrc,inds) = -ws(isrc) * grad_src(inds);
- Gsrc(isrc,inds) = rand(1,length(inds));
- end
- end
- G = [Gstr Gsrc];
-
- Cm = diag(rand(M,1));
- Cmi = P'*Cm*P;
- GG = G*Cm*P;
-
- d = zeros(S,1);
- Hnew = GG'*GG + Cmi;
- dnew = GG'*d;
- mu = Hnew*dnew;
- whos d Hnew dnew mu
-
- figure; nr=3; nc=4; msize=3;
- subplot(nr,nc,2); spy(G,msize); title(sprintf('G is %i by %i',size(G)));
- subplot(nr,nc,3); spy(G',msize); title(sprintf('G^T is %i by %i',size(G')));
- subplot(nr,nc,4); spy(G*G',msize); title(sprintf('G G^T is %i by %i',size(G*G')));
-
- subplot(nr,nc,5); spy(P',msize); title(sprintf('P^T is %i by %i',size(P')));
- subplot(nr,nc,6); spy(Cm,msize); title(sprintf('Cm is %i by %i',size(Cm)));
- subplot(nr,nc,7); spy(P,msize); title(sprintf('P is %i by %i',size(P)));
- subplot(nr,nc,8); spy(Cmi,msize); title(sprintf('P^T Cm P is %i by %i',size(Cmi)));
-
- subplot(nr,nc,9); spy(G,msize); title(sprintf('G is %i by %i',size(G)));
- subplot(nr,nc,10); spy(Cm,msize); title(sprintf('Cm is %i by %i',size(Cm)));
- subplot(nr,nc,11); spy(P,msize); title(sprintf('P is %i by %i',size(P)));
- subplot(nr,nc,12); spy(GG,msize); title(sprintf('G Cm P is %i by %i',size(GG)));
- orient tall, wysiwyg
-
- break
-end
-
-%---------------------------------------------
-
-if 0==1
- nk = 100;
- kvec = 10.^linspace(-4,4,nk);
-
- figure; nr=2; nc=1;
- for jj=1:2
- switch jj
- case 1, ndata = 5; nparm = 10;
- case 2, ndata = 10; nparm = 5;
- end
- Cm = diag(rand(nparm,1));
- G = rand(ndata,nparm);
- d = rand(ndata,1);
-
- norm_dm = zeros(nk,1);
- for kk = 1:nk
- C = Cm * kvec(kk);
- norm_dm(kk) = norm( C * G' * inv(G*C*G') * d );
- end
- subplot(nr,nc,jj);
- semilogx(kvec, norm_dm, '.');
- end
-end
-
-%----------------------------------------------
-
irun0 = input(' Enter irun0 : ');
iread = input(' Enter 1 to read in new event kernels (0 otherwise) : ');
-hmod = input(' Enter next model number (hmod) : ');
+hmod = input(' Enter next model number (hmod = 2,4,6,8,...) : ');
INV_STRUCT = input(' Enter 1 to invert for STRUCTURE : ');
INV_SOURCE = input(' Enter 1 to invert for SOURCE : ');
iwrite = input(' Enter 1 to write out files : ');
@@ -175,9 +71,12 @@
%sval_cut = SIGMA_DT * sqrt(nmeas);
sval_cut = 5;
+% event IDs
+for ii=1:nsrc, eids(ii) = cellstr(sprintf(stfm,ii)); end
+
%----------------
-stm = sprintf(stfmt,hmod);
+stm = sprintf(stfm,hmod);
% directories
dir0 = [dir_run 'run_' sprintf(stfm,irun0) '/'];
@@ -188,15 +87,13 @@
if ~exist(dir2), mkdir(dir2); end
dir2lab = ['run-' sprintf(stfm,irun0) ', READ-IN -- toward model m' stm ];
-if hmod == 1
+if hmod == 2
dir1 = dir0;
else
- dir1 = [dir0 'READ_IN/model_m' sprintf(stfmt,hmod-1) '/'];
+ dir1 = [dir0 'READ_IN/model_m' sprintf(stfm,hmod-2) '/'];
end
-% indexing for measurements
-itemp = [nrec:nrec:nmeas]';
-indmat1 = [ [1 ; 1+itemp(1:end-1)] itemp ];
+% indexing for source parameters
itemp = NPARM_SOURCE*[1:25]';
indmat2 = [ [1 ; 1+itemp(1:end-1)] itemp ];
@@ -227,19 +124,10 @@
rlon = rec_lonlat(:,1);
rlat = rec_lonlat(:,2);
-% load the data covariance matrix
-cov_data = load([dir1 'cov_data_diagonal.dat']);
-
-% load all the measurements and partition into matrix
-meas_all = load([dir1 'measure_vec.dat']);
-dT_all = meas_all(:,1); % WITH ERRORS ADDED
-dT_all_norm = dT_all.^2 ./ cov_data;
-
% load the covariance matrix
% --> cov_imetric(NLOCAL+1 : nmod_str) = ( sigma_beta )**2 / da_local_vec(:) * AREA
%cov_model = load('/net/denali/scratch1/carltape/OUTPUT_2/run_9100/cov_imetric_diagonal.dat');
-cov_model_all = load([dir1 'cov_model_diagonal.dat']);
-cov_model = cov_model_all(:,1);
+[cov_model, cov_imetric] = textread([dir1 'cov_model_diagonal.dat'],'%f%f');
cov_beta = cov_model(inds_B)';
cov_src = cov_model(inds_src)';
clear cov_model_all
@@ -259,6 +147,9 @@
m_all = mtemp(:,1);
m_src = mtemp(inds_src)';
+% prior, initial, and target models
+[mprior, m0_initial, mtarget] = textread([dir0 'prior_initial_target_models.dat'],'%f%f%f');
+
% load the structure files
mtemp = load([dir1 'structure_syn.dat']);
m_str_lon = mtemp(:,1);
@@ -306,30 +197,21 @@
data_norm = load([dir1 'data_norm.dat']);
model_norm = load([dir1 'model_norm.dat']);
-disp(' '); disp(' CHECKING VARIOUS NORMS:');
-disp(' model norm:');
-model_norm, sum( m_all.^2 ./ cov_model )
-disp(' data norm:');
-data_norm, sum(dT_all_norm)
-disp(' misfit function value:');
-chi, 0.5*( sum( dT_all_norm ) + sum( m_all.^2 ./ cov_model ))
-
%====================================================
-% partition DATA vector of traveltime measurements into a matrix,
-% since in this case we know that there are the same number of picks for each event
-dT_mat = zeros(nsrc,nrec);
-dT_norm_mat = zeros(nsrc,nrec);
-for isrc = 1:nsrc
- inds = [indmat1(isrc,1) : indmat1(isrc,2)];
- dT_mat(isrc,:) = dT_all(inds);
- dT_norm_mat(isrc,:) = dT_all_norm(inds);
-end
+% load the data covariance matrix
+cov_data = load([dir1 'cov_data_diagonal.dat']);
-% compute the new data vector
-dnorm = sum( dT_norm_mat, 2);
-if any(dnorm==0), error(' For at least one source, there is perfect fit.'); end
+% load all the measurements and partition into matrix
+meas_all = load([dir1 'measure_vec.dat']);
+dT_all = meas_all(:,1); % WITH ERRORS ADDED
+% vector of number of measurement windows per event
+Ns_vec = nrec * ones(nsrc,1); % KEY: same measurements per event
+
+% KEY: compute projected data vector
+[dnorm, indmat] = wave2d_dsub(dT_all,cov_data,Ns_vec);
+
% compute the weights (SIGN OR NOT?)
%ws = zeros(nsrc,1);
%ws = 1 ./ sqrt( sum( dT_norm_mat, 2) );
@@ -337,9 +219,21 @@
% compute the new data vector
%dnorm = 1 ./ ws;
-% check
-sum(sum(dT_norm_mat)), sum(dT_all_norm), sum( dnorm )
+disp(' '); disp(' CHECKING VARIOUS NORMS:');
+disp(' model norm:');
+model_norm, sum( (m_all-mprior).^2 ./ cov_model )
+disp(' data norm:');
+data_norm, sum(dnorm)
+disp(' misfit function value:');
+chi, 0.5*( sum( dnorm ) + sum( (m_all-mprior).^2 ./ cov_model ))
+%--------------------------------------------------------------------------
+% ASSEMBLING THE GRADIENTS
+
+% initialize
+Gstr = zeros(nsrc,NLOCAL);
+Gsrc = zeros(nsrc,nmod_src);
+
if INV_STRUCT == 1
% load the jacobian for constructing the "event gradient" from the event kernel
@@ -350,20 +244,19 @@
%iread = 1;
efile = 'wave2d_kernel';
if iread==1
- % THE SPECIFICATION FOR SMOOTHED KERNELS IS DONE IN WAVE2D.F90
- % EITHER WAY, YOU LOAD THE FILES kernel_basis_smooth, WHICH MAY OR
- % MAY NOT BE SMOOTHED.
ismooth = input(' Enter 1 to read smoothed event kernels (0 otherwise) : ');
disp('reading in the event kernels...');
Kall = zeros(nsrc,NLOCAL);
for isrc = 1:nsrc
isrc
- dirK = [dir1 'event_' sprintf(stfmt,isrc) '/'];
- if ismooth == 1
- kernel = load([dirK 'kernel_basis_smooth']); Kbeta = kernel(:,3)';
- else
- kernel = load([dirK 'kernel_basis']); Kbeta = kernel(:,7)';
- end
+ dirK = [dir1 'event_' sprintf('%3.3i',isrc) '/'];
+ kernel = load([dirK 'kernel_basis_smooth']); Kbeta = kernel(:,3)';
+ if ismooth==0, Kbeta = kernel(:,4); else Kbeta = kernel(:,3); end
+ %if ismooth == 1
+ % kernel = load([dirK 'kernel_basis_smooth']); Kbeta = kernel(:,3)';
+ %else
+ % kernel = load([dirK 'kernel_basis']); Kbeta = kernel(:,7)';
+ %end
lon = kernel(:,1);
lat = kernel(:,2);
Kall(isrc,:) = Kbeta;
@@ -376,7 +269,6 @@
% construct G (plot event kernels too)
%nsrc = 5;
- Gstr = zeros(nsrc,NLOCAL);
for isrc = 1:nsrc
Kbeta = Kall(isrc,:);
%Gstr(isrc,:) = -ws(isrc) * Kbeta .* Ai; % SIGN OR NOT?
@@ -396,43 +288,55 @@
end
end
- % Hessian for structure parameters (from event kernels)
- disp('constructing the Hessian...');
- Hstr = zeros(nsrc,nsrc);
- for ii = 1:nsrc
- for jj = 1:nsrc
- Hstr(ii,jj) = dot(Gstr(ii,:), cov_beta.*Gstr(jj,:));
- end
- end
+% % Hessian for structure parameters (from event kernels)
+% disp('constructing the Hessian...');
+% Hstr = zeros(nsrc,nsrc);
+% for ii = 1:nsrc
+% for jj = 1:nsrc
+% Hstr(ii,jj) = dot(Gstr(ii,:), cov_beta.*Gstr(jj,:));
+% end
+% end
-else
- Hstr = zeros(nsrc,nsrc);
end
-% Hessian and projected gradient for source parameters
-Hsrc = zeros(nsrc,nsrc);
+% projected gradient for source parameters
+%Hsrc = zeros(nsrc,nsrc);
%Hsrc_vec = zeros(nsrc,1);
if INV_SOURCE == 1
- Gsrc = zeros(nsrc,nmod_src);
+ % QUESTION: SHOULD THIS BE THE FULL GRADIENT OR JUST THE DATA TERM?
Gsrc = [diag(gradient(inds_ts)) diag(gradient(inds_xs)) diag(gradient(inds_ys))];
- %for isrc = 1:nsrc
- % inds = indmat2(isrc,1) : indmat2(isrc,2);
- % %Gsrc(isrc,inds) = -ws(isrc) * grad_src(inds);
- % Gsrc(isrc,inds) = grad_src(inds);
- %end
- Hsrc = Gsrc * diag(cov_src) * transpose(Gsrc);
- %for isrc = 1:nsrc
- % inds = indmat2(isrc,1) : indmat2(isrc,2);
- % Hsrc_vec(isrc) = ws(isrc)^2 * dot(grad_src(inds), cov_src(inds).*grad_src(inds) );
- %end
- %Hsrc = diag(Hsrc_vec);
+% %for isrc = 1:nsrc
+% % inds = indmat2(isrc,1) : indmat2(isrc,2);
+% % %Gsrc(isrc,inds) = -ws(isrc) * grad_src(inds);
+% % Gsrc(isrc,inds) = grad_src(inds);
+% %end
+% Hsrc = Gsrc * diag(cov_src) * transpose(Gsrc);
+%
+% %for isrc = 1:nsrc
+% % inds = indmat2(isrc,1) : indmat2(isrc,2);
+% % Hsrc_vec(isrc) = ws(isrc)^2 * dot(grad_src(inds), cov_src(inds).*grad_src(inds) );
+% %end
+% %Hsrc = diag(Hsrc_vec);
+
end
+%--------------------------------------------------------------------------
+% COMPUTE THE HESSIAN IN THE SUBSPACE OF SOURCES
+
+[Hstr,Hsrc] = wave2d_compute_hessian(Gstr,Gsrc,cov_beta,cov_src);
+
% overall Hessian
+% NOTE: identity matrix term
H = Hstr + Hsrc + eye(nsrc,nsrc);
+% write out files for Hessian
+if iwrite==1
+ nprint = 20;
+ wave2d_write_hessian(H,eids,dir2,nprint);
+end
+
% construct projection matrix (and assign gradient)
if INV_STRUCT == 1
G = Gstr;
@@ -478,43 +382,9 @@
disp(' gradient balance for the Hessian (subspace) inversion (mean of last column) :');
disp(mean( diag(Hsrc)./diag(Hstr) ))
-disp(' properties of Hessian (min, median, mean(abs), max, std):');
-stH = sprintf('min %.2e, median %.2e, mean(abs) %.2e, max %.2e, std %.2e',...
- min(H(:)), median(H(:)), mean(abs(H(:))), max(H(:)), std(H(:)));
-disp(stH);
+% plot Hessian matrices
+wave2d_plot_hessian(Hstr,Hsrc,H);
-if ~and(INV_STRUCT==1, INV_SOURCE==1)
- figure;
- %pcolor(H); shading flat; xlabel('Row index'); ylabel('Column index');
- imagesc(H); ylabel('Row index'); xlabel('Column index');
- title({'Hessian (symmetric matrix)',stH});
- map = colormap('gray'); colormap(flipud(map));
- colorbar; axis equal; axis tight;
-else
- figure; nr=3; nc=1;
- subplot(nr,nc,1);
- %pcolor(H); shading flat; xlabel('Row index'); ylabel('Column index');
- imagesc(H); ylabel('Row index'); xlabel('Column index');
- title({'Hessian (symmetric matrix)',stH});
- map = colormap('gray'); colormap(flipud(map));
- colorbar; axis equal; axis tight;
-
- subplot(nr,nc,2);
- %pcolor(Hstr); shading flat; xlabel('Row index'); ylabel('Column index');
- imagesc(Hstr); ylabel('Row index'); xlabel('Column index');
- title('Hessian for structure (symmetric matrix)');
- map = colormap('gray'); colormap(flipud(map));
- colorbar; axis equal; axis tight;
-
- subplot(nr,nc,3);
- %pcolor(Hsrc); shading flat; xlabel('Row index'); ylabel('Column index');
- imagesc(Hsrc); ylabel('Row index'); xlabel('Column index');
- title('Hessian for source (diagonal matrix)');
- map = colormap('gray'); colormap(flipud(map));
- colorbar; axis equal; axis tight;
- orient tall, wysiwyg
-end
-
if INV_SOURCE == 1
figure; nr=2; nc=1;
subplot(nr,nc,1); spy(Gsrc); title('Gsrc');
@@ -602,6 +472,8 @@
%H = H / 0.5; % 9550
+ error('replace the following code with wave2d_tsvd_hessian.m');
+
pinds = [1:nsrc]';
[mu_all,rss,f_r_ss] = tsvd(dnorm,H,pinds); % KEY: TSVD
@@ -672,7 +544,7 @@
% KEY LOOP OVER DIFFERENT MODELS (EACH OBTAINED FROM THE SAME EVENT KERNELS)
for ip = 1:nump
- stip = sprintf(stfmt,ip);
+ stip = sprintf(stfm,ip);
% OPTIONS FOR INVERTING HESSIAN:
% (1) truncated singular value decomposition
@@ -755,8 +627,10 @@
dm_str = cov_beta' .* (transpose(Gstr) * mu);
end
end
- %disp(' norm of structure update:');
- %sum( dm_str .* dm_str ./ cov_beta' )
+ disp(' norm, min, max of structure update:');
+ norm_dm_str = sum( dm_str .* dm_str ./ cov_beta' )
+ min_dm_str = min( dm_str )
+ max_dm_str = max( dm_str )
% KEY: solve for the SOURCE model update (ts1, ts2, ..., xs1, xs2, ..., ys1, ys2, ...)
dm_src = zeros(nmod_src,1);
@@ -771,8 +645,10 @@
dm_src = cov_src' .* (transpose(Gsrc) * mu);
end
end
- %disp(' norm of source update:');
- %sum( dm_src .* dm_src ./ cov_src' )
+ disp(' norm, min, max of source update:');
+ norm_dm_src = sum( dm_src .* dm_src ./ cov_src' )
+ min_dm_src = min( dm_src )
+ max_dm_src = max( dm_src )
%-------------------
@@ -848,7 +724,7 @@
%disp(' norm of new structure model:');
%sum( m_str_B_new .* m_str_B_new ./ cov_beta' )
- % plot structure model update and new model
+ % PLOT structure model update and new model
if and(INV_STRUCT == 1, ifigure == 1)
tlabs = {{dir2lab,'data'},'m00',['dm (' stplab ')'],['m01 (' stplab ')']};
cmax = 0.1;
@@ -861,13 +737,14 @@
case 4, mplot = m_str_B_new;
end
[X,Y,Z] = griddataXB(lon,lat,mplot,100,'nearest');
+ mnorm = sum( mplot.^2 ./ cov_beta' );
subplot(nr,nc,kk); hold on;
pcolor(X,Y,Z); shading interp;
caxis([-1 1]*cmax); colormap(seis); axis equal, axis tight
plot(rlon,rlat,'k.','markersize',16);
plot(slon,slat,'p','markersize',18,'markeredgecolor','k','linewidth',2);
- title(tlabs{kk});
+ title([tlabs{kk} sprintf(', norm = %.2f',mnorm)]);
end
% TSVD plots
@@ -973,6 +850,9 @@
end
end
+%======================================================================
+% EXTRA CODE FOR TESTING
+
%------------------------------------------------
% extra code for plotting misfit values
@@ -983,14 +863,14 @@
for h = 1:4
strun0 = num2str(run0);
- stm = sprintf(stfmt,h);
+ stm = sprintf(stfm,h);
dir_run = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_OUTPUT/';
idir = [dir_run 'run_' strun0 '/READ_IN/model_m' stm '/'];
load([idir 'wave2d_subspace_matlab_all']);
clear data_norm model_norm chi_total pmax_vec
for ii = 1:25
- stp = sprintf(stfmt,ii);
+ stp = sprintf(stfm,ii);
dir1 = [ idir 'run_p' stp '/'];
chi_file = [dir1 'data_norm.dat'];
if exist(chi_file)
@@ -1031,4 +911,173 @@
end
+%----------------------------------------------
+
+if 0==1
+ temp = load('/net/denali/scratch1/carltape/OUTPUT/run_7000/structure_dat.dat');
+ %temp = load('/net/denali/scratch1/carltape/OUTPUT/run_7000/READ_IN/structure_dat.dat');
+ lon = temp(:,1); lat = temp(:,2); beta = temp(:,7);
+
+ [X,Y,Z] = griddataXB(lon,lat,beta,100,'nearest');
+ figure; cmax = 0.1; hold on;
+ pcolor(X,Y,Z); shading interp;
+ caxis([-1 1]*cmax); colormap(seis);
+end
+
+%----------------------------------------------
+
+% check the relative dimensions of matrices and plot a schematic figure
+if 0==1
+ S = 5; % number of sources
+ Cstr = 1; % number of structure parameter classes
+ Csrc = 2; % number of source parameter classes
+ C = Cstr + Csrc;
+ Mstr1 = 20; % number of structure parameters
+ Msrc1 = 1; % number of source parameters in class 1 (origin time)
+ Msrc2 = 2; % number of source parameters in class 2 (xy, ys)
+ Mstr = Mstr1;
+ Msrc = S*(Msrc1 + Msrc2);
+ M = Mstr + Msrc; % number of rows in P
+ N = C*S; % number of columns in P
+
+ P = zeros(M,N);
+ P(1:Mstr1,1:S) = rand(Mstr1,S);
+ for s = 1:S
+ P(Mstr1+s,S+s) = rand;
+ P(Mstr1+S+s,2*S+s) = rand;
+ P(Mstr1+2*S+s,2*S+s) = rand;
+ end
+
+ % fill G
+ G = zeros(S,M);
+ Gstr = rand(S,Mstr);
+ Gsrc = zeros(S,Msrc);
+ if 1==1
+ Gsrc = repmat(diag(rand(S,1)),1,Msrc1+Msrc2);
+ else
+ itemp = (Msrc1 + Msrc2)*[1:S]';
+ indmat2 = [ [1 ; 1+itemp(1:end-1)] itemp ];
+ for isrc = 1:S
+ inds = indmat2(isrc,1) : indmat2(isrc,2);
+ %Gsrc(isrc,inds) = -ws(isrc) * grad_src(inds);
+ Gsrc(isrc,inds) = rand(1,length(inds));
+ end
+ end
+ G = [Gstr Gsrc];
+
+ Cm = diag(rand(M,1));
+ Cmi = P'*Cm*P;
+ GG = G*Cm*P;
+
+ d = zeros(S,1);
+ Hnew = GG'*GG + Cmi;
+ dnew = GG'*d;
+ mu = Hnew*dnew;
+ whos d Hnew dnew mu
+
+ figure; nr=3; nc=4; msize=3;
+ %subplot(nr,nc,2); spy(G,msize); title(sprintf('G is %i by %i',size(G)));
+ %subplot(nr,nc,3); spy(G',msize); title(sprintf('G^T is %i by %i',size(G')));
+ %subplot(nr,nc,4); spy(G*G',msize); title(sprintf('G G^T is %i by %i',size(G*G')));
+ subplot(nr,nc,1); spy(G,msize); title(sprintf('G is %i by %i',size(G)));
+ subplot(nr,nc,2); spy(Cm,msize); title(sprintf('Cm is %i by %i',size(Cm)));
+ subplot(nr,nc,3); spy(G',msize); title(sprintf('G^T is %i by %i',size(G')));
+ subplot(nr,nc,4); spy(G*Cm*G',msize); title(sprintf('G Cm G^T is %i by %i',size(G*Cm*G')));
+
+ subplot(nr,nc,5); spy(P',msize); title(sprintf('P^T is %i by %i',size(P')));
+ subplot(nr,nc,6); spy(Cm,msize); title(sprintf('Cm is %i by %i',size(Cm)));
+ subplot(nr,nc,7); spy(P,msize); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,8); spy(Cmi,msize); title(sprintf('P^T Cm P is %i by %i',size(Cmi)));
+
+ subplot(nr,nc,9); spy(G,msize); title(sprintf('G is %i by %i',size(G)));
+ subplot(nr,nc,10); spy(Cm,msize); title(sprintf('Cm is %i by %i',size(Cm)));
+ subplot(nr,nc,11); spy(P,msize); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,12); spy(GG,msize); title(sprintf('G Cm P is %i by %i',size(GG)));
+ orient tall, wysiwyg
+
+ % projection for data
+ S = 5;
+ N = 15;
+ Ns = N/S;
+
+ P = zeros(S,N);
+ for s = 1:S
+ P(s, Ns*(s-1) + [1:Ns]) = rand;
+ end
+ d = randn(N,1);
+ C = rand(N,N);
+ Cblock = P'*P;
+ Cdiag = diag(rand(N,1));
+
+ figure; nr=4; nc=4;
+ subplot(nr,nc,1); spy(P); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,2); spy(P'); title(sprintf('P^T is %i by %i',size(P')));
+ subplot(nr,nc,3); spy(P*P'); title(sprintf('P*P^T is %i by %i',size(P*P')));
+
+ subplot(nr,nc,5); spy(P'); title(sprintf('P^T is %i by %i',size(P')));
+ subplot(nr,nc,6); spy(P); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,7); spy(P'*P); title(sprintf('P^T*P is %i by %i',size(P'*P)));
+
+ subplot(nr,nc,9); spy(C'); title(sprintf('C^T is %i by %i',size(C')));
+ subplot(nr,nc,10); spy(P'*P); title(sprintf('P^T*P is %i by %i',size(P'*P)));
+ subplot(nr,nc,11); spy(C); title(sprintf('C is %i by %i',size(C)));
+ subplot(nr,nc,12); spy(C'*P'*P*C); title(sprintf('C^T*P^T*P*C is %i by %i',size(C'*P'*P*C)));
+
+ subplot(nr,nc,13); spy(Cdiag'); title(sprintf('C^T is %i by %i',size(Cdiag')));
+ subplot(nr,nc,14); spy(P'*P); title(sprintf('P^T*P is %i by %i',size(P'*P)));
+ subplot(nr,nc,15); spy(Cdiag); title(sprintf('C is %i by %i',size(Cdiag)));
+ subplot(nr,nc,16); spy(Cdiag'*P'*P*Cdiag); title(sprintf('C^T*P^T*P*C is %i by %i',size(Cdiag'*P'*P*Cdiag)));
+ orient tall, wysiwyg
+ %print(gcf,'-dpsc', 'test')
+
+ figure; nr=4; nc=4;
+ subplot(nr,nc,1); spy(P); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,2); spy(C); title(sprintf('C is %i by %i',size(C)));
+ subplot(nr,nc,3); spy(d); title(sprintf('d is %i by %i',size(d)));
+ subplot(nr,nc,4); spy(P*C*d); title(sprintf('P*C*d is %i by %i',size(P*C*d)));
+
+ subplot(nr,nc,5); spy(P); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,6); spy(Cblock); title(sprintf('C is %i by %i',size(Cblock)));
+ subplot(nr,nc,7); spy(d); title(sprintf('d is %i by %i',size(d)));
+ subplot(nr,nc,8); spy(P*Cblock*d); title(sprintf('P*C*d is %i by %i',size(P*Cblock*d)));
+
+ subplot(nr,nc,9); spy(P); title(sprintf('P is %i by %i',size(P)));
+ subplot(nr,nc,10); spy(Cdiag); title(sprintf('C is %i by %i',size(Cdiag)));
+ subplot(nr,nc,11); spy(d); title(sprintf('d is %i by %i',size(d)));
+ subplot(nr,nc,12); spy(P*Cdiag*d); title(sprintf('P*C*d is %i by %i',size(P*Cdiag*d)));
+
+ subplot(nr,nc,14); spy(P*Cdiag); title(sprintf('P*C is %i by %i',size(P*Cdiag)));
+ subplot(nr,nc,15); spy(Cdiag*d); title(sprintf('C*d is %i by %i',size(Cdiag*d)));
+ orient tall, wysiwyg
+ %print(gcf,'-dpsc', 'test')
+
+ break
+end
+
+%---------------------------------------------
+
+if 0==1
+ nk = 100;
+ kvec = 10.^linspace(-4,4,nk);
+
+ figure; nr=2; nc=1;
+ for jj=1:2
+ switch jj
+ case 1, ndata = 5; nparm = 10;
+ case 2, ndata = 10; nparm = 5;
+ end
+ Cm = diag(rand(nparm,1));
+ G = rand(ndata,nparm);
+ d = rand(ndata,1);
+
+ norm_dm = zeros(nk,1);
+ for kk = 1:nk
+ C = Cm * kvec(kk);
+ norm_dm(kk) = norm( C * G' * inv(G*C*G') * d );
+ end
+ subplot(nr,nc,jj);
+ semilogx(kvec, norm_dm, '.');
+ end
+end
+
%======================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2011-04-13 16:44:53 UTC (rev 18230)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2011-04-15 19:59:48 UTC (rev 18231)
@@ -211,7 +211,7 @@
print *, 'For READ_IN runs...'
if (NITERATION /= 0) stop 'NITERATION = 0'
if (ISRC_SPACE /= 6) stop 'ISRC_SPACE = 6'
- hmod = 6 ! iteration number for read-in models
+ hmod = 1 ! iteration number for read-in models
endif
if ( IKER/=0 .and. IKER/=1 .and. IKER/=2 ) stop 'IKER must = 0,1,2'
@@ -469,6 +469,7 @@
!============================================
! LOOP 1: different tomographic runs
do iq = 1,qmax
+ !do iq = 1,250 ! GAUSS
!============================================
! KEY COMMAND: scalelength of checker for velocity models (1,2,3)
@@ -480,6 +481,8 @@
irun0 = IRUNZ + 20*(iq-1) ! increment is 20
endif
+ !irun0 = IRUNZ + iq ! GAUSS
+
! name the reference output directory for the optimization run
if (READ_IN) then
@@ -793,11 +796,11 @@
! 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_002/'
+ idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/'
! read in structure model for synthetics
- write(filename2,'(a)') trim(idir)//'structure_syn_in.dat'
+ write(filename2,'(a)') trim(idir)//'model_0001/structure_syn_in.dat'
+ !write(filename2,'(a,i4.4,a)') trim(idir)//'model_set_0001/structure_model_',iq,'.dat'
open(unit=19,file=filename2,status='old')
do ispec = 1,NSPEC
do j = 1,NGLLZ
@@ -812,7 +815,7 @@
close(19)
! read in structure model for data
- write(filename2,'(a)') trim(idir)//'structure_dat_in.dat'
+ write(filename2,'(a)') trim(idir)//'model_0001/structure_dat_in.dat'
open(unit=20,file=filename2,status='old')
do ispec = 1,NSPEC
do j = 1,NGLLZ
@@ -827,7 +830,7 @@
close(20)
! read in sigma value describing the distribution
- write(filename,'(a)') trim(idir)//'structure_sigma.dat'
+ write(filename,'(a)') trim(idir)//'model_0001/structure_sigma.dat'
open(unit=21,file=filename,status='old')
read(21,*) sigma_beta0
close(21)
@@ -838,7 +841,7 @@
beta_dat = sqrt( mu_dat / rho_dat )
! read reference values from file -- these were used to make the perturbed field
- idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/meshes/mesh_001/'
+ idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/meshes/mesh_0001/'
open(unit=16,file=trim(idir)//'reference_values.dat', status='unknown')
read(16,'(6e20.10)') alpha0, beta0, rho0, bulk0, kappa0, mu0
close(16)
@@ -3231,7 +3234,8 @@
do j = 1,NGLLZ
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
- write(13,'(3e16.6)') x_plot(iglob), z_plot(iglob), kbeta_smooth(i,j,ispec)
+ write(13,'(4e16.6)') x_plot(iglob), z_plot(iglob), &
+ kbeta_smooth(i,j,ispec), btype_kernel(i,j,ispec)
ktemp(k) = kbeta_smooth(i,j,ispec) ! used to compute norm of gradient
k = k+1
@@ -3935,7 +3939,7 @@
deallocate(g0,gt,gk,p0,pk)
deallocate(source_gradient)
deallocate(m0_vec,mt_vec)
- deallocate(m_src_syn,m_src_dat)
+ deallocate(m_src_syn,m_src_dat,m_src_prior)
!deallocate(m_src_syn_vec,m_src_dat_vec)
!deallocate(m0_vec_initial,m_src_syn_vec_initial,m_src_syn_initial)
!deallocate(m_scale_src_all)
More information about the CIG-COMMITS
mailing list