[cig-commits] r11845 - in seismo/3D/mt_measure_adj: . PLOTS scripts_tomo scripts_tomo/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Tue Apr 22 11:56:23 PDT 2008
Author: carltape
Date: 2008-04-22 11:56:23 -0700 (Tue, 22 Apr 2008)
New Revision: 11845
Added:
seismo/3D/mt_measure_adj/PLOTS/RECON/
seismo/3D/mt_measure_adj/PLOTS/RECON_TEST/
seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl
seismo/3D/mt_measure_adj/scripts_tomo/matlab/mu_dirs/
seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m
Modified:
seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl
seismo/3D/mt_measure_adj/PLOTS/plot_win_adj_all.pl
seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl
seismo/3D/mt_measure_adj/mt_measure_adj.f90
seismo/3D/mt_measure_adj/mt_sub.f90
seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m
seismo/3D/mt_measure_adj/scripts_tomo/matlab/display_meas.m
seismo/3D/mt_measure_adj/scripts_tomo/matlab/kernels_done_mod
seismo/3D/mt_measure_adj/scripts_tomo/run_mt_cc_plot.pl
seismo/3D/mt_measure_adj/scripts_tomo/run_mt_measure_adj.pl
seismo/3D/mt_measure_adj/scripts_tomo/run_tomo.pl
Log:
Updated code to output cross-correlation-shifted synthetics for plotting purposes. Now, the plot scripts shows data, synthetics, and CC-shifted (traveltime and amplitude) synthetics, in addition to the adjoint sources and windows. Also added several Matlab scripts in scripts_tomo/matlab for computing the misfit function, and for computing the Hessian for the subspace method.
Modified: seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -11,12 +11,12 @@
# (Search for "USER".)
#
# EXAMPLE:
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT.WINDOWS
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -l 0/200 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT.WINDOWS
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -r -l 0/200 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT.WINDOWS
+# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -r RECON -w MEASUREMENT.WINDOWS
+# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -l 0/200 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -r RECON -w MEASUREMENT.WINDOWS
+# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -r -l 0/200 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -r RECON -w MEASUREMENT.WINDOWS
#
-# plot_win_adj.pl -m ../CMTSOLUTION_9818433 -k 1 -l 0/150 -n CLC/CI -a STATIONS_ADJOINT -d DATA_TEST -s SYN_TEST -w MEASUREMENT.WINDOWS
-# plot_win_adj.pl -m ../CMTSOLUTION_14186612 -k 2 -r -l 0/200 -n SCI2/CI -a STATIONS_ADJOINT -d DATA_TEST -s SYN_TEST -w MEASUREMENT.WINDOWS
+# plot_win_adj.pl -m ../CMTSOLUTION_9818433 -k 1 -l 0/150 -n CLC/CI -a STATIONS_ADJOINT -d DATA_TEST -s SYN -r RECON_TEST -w MEASUREMENT.WINDOWS
+# plot_win_adj.pl -m ../CMTSOLUTION_14186612 -k 2 -r -l 0/200 -n SCI2/CI -a STATIONS_ADJOINT -d DATA_TEST -s SYN -r RECON_TEST -w MEASUREMENT.WINDOWS
#
# From the "EXAMPLE" directory:
# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -l 0/200 -n CLC/CI -a STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT.WINDOWS
@@ -29,7 +29,7 @@
sub Usage{
print STDERR <<END;
- plot_win_adj.pl -m CMTFILE -n station_name/network_name -b iboth -r -a STATION_FILE -d data_dir -s syn_dir -w winfile
+ plot_win_adj.pl -m CMTFILE -n station_name/network_name -b iboth -r -a STATION_FILE -d data_dir -s syn_dir -c recon_dir -w winfile
with
-m -- use CMTSOLUTION file to plot event focal mechanism (beach ball) on the map
@@ -41,13 +41,14 @@
-a -- use the station file to plot all the stations on the map
-d -- data directory
-s -- synthetics directory
+ -c -- reconstructed directory
-w -- the file with the windows
END
exit(1);
}
if (@ARGV == 0) { Usage(); }
-if (!getopts('k:m:l:n:b:ra:d:s:w:')) {die('Check input arguments\n');}
+if (!getopts('k:m:l:n:b:ra:d:s:c:w:')) {die('Check input arguments\n');}
print "STATIONS file is $opt_a\n";
if ($opt_m) {@cmtvals = &get_cmt($opt_m);}
#if ($opt_k) {$iker = $opt_k;} # does not work for iker = 0
@@ -58,6 +59,7 @@
if(!$opt_b) {$iboth = 0} else {$iboth = $opt_b}
if (!$opt_d) {$opt_d = "DATA"}
if (!$opt_s) {$opt_s = "SYN"}
+if (!$opt_c) {$opt_c = "RECON"}
# name of file containing measurements (see mt_sub.f90)
$meas_file = "window_chi";
@@ -119,7 +121,7 @@
# plot the base map
$origin = "-X3.0";
-`psbasemap $bounds $proj $tick $origin -K > $ps_file`; # START
+`psbasemap $bounds $proj $tick $origin -K -P > $ps_file`; # START
`pscoast $bounds $proj $cinfo -S150/200/255 -W2 -G200/255/150 -N1/1p -N2/0.5p -O -K >> $ps_file`;
# plot southern California faults
@@ -138,22 +140,22 @@
`awk \'{print \$4,\$3,0,4,B,C}\' $opt_a | psxy $proj $bounds -M -St0.15 -W2 -G200/200/200 -N -O -K >>$ps_file`;
#------------
-# data and synthetics
+# data and synthetics and reconstructed synthetics
$suffix = ".d"; # USER MUST CHANGE
# DATA
$d_tag = "${opt_d}/$eid.$net.$sta";
-$Tdat=`ls ${d_tag}.*T.sac${suffix}`; chomp($Tdat);
-$Rdat=`ls ${d_tag}.*R.sac${suffix}`; chomp($Rdat);
-$Zdat=`ls ${d_tag}.*Z.sac${suffix}`; chomp($Zdat);
-#$Tdat="$opt_d/$eid.$net.$sta.BHT.sac${suffix}";
-#$Rdat="$opt_d/$eid.$net.$sta.BHR.sac${suffix}";
-#$Zdat="$opt_d/$eid.$net.$sta.BHZ.sac${suffix}";
+$Tdat = `ls ${d_tag}.*T.sac${suffix}`; chomp($Tdat);
+$Rdat = `ls ${d_tag}.*R.sac${suffix}`; chomp($Rdat);
+$Zdat = `ls ${d_tag}.*Z.sac${suffix}`; chomp($Zdat);
+#$Tdat = "$opt_d/$eid.$net.$sta.BHT.sac${suffix}";
+#$Rdat = "$opt_d/$eid.$net.$sta.BHR.sac${suffix}";
+#$Zdat = "$opt_d/$eid.$net.$sta.BHZ.sac${suffix}";
if (0==1) {
- ($Tdat,$Tchan)=split(" ",`$saclst kcmpnm f ${d_tag}*T.sac${suffix}`);
- ($Rdat,$Rchan)=split(" ",`$saclst kcmpnm f ${d_tag}*R.sac${suffix}`);
- ($Zdat,$Zchan)=split(" ",`$saclst kcmpnm f ${d_tag}*Z.sac${suffix}`);
+ ($Tdat,$Tchan) = split(" ",`$saclst kcmpnm f ${d_tag}*T.sac${suffix}`);
+ ($Rdat,$Rchan) = split(" ",`$saclst kcmpnm f ${d_tag}*R.sac${suffix}`);
+ ($Zdat,$Zchan) = split(" ",`$saclst kcmpnm f ${d_tag}*Z.sac${suffix}`);
print "\n $Tdat $Tchan \n $Rdat $Rchan \n $Zdat $Zchan \n";
}
if(-f $Tdat) {$Tflag = 1} else {$Tflag = 0}
@@ -167,15 +169,22 @@
$Tlab = "$sta.$net.BHT";
$Rlab = "$sta.$net.BHR";
$Zlab = "$sta.$net.BHZ";
-$Tsyn="$opt_s/${Tlab}.semd.sac${suffix}";
-$Rsyn="$opt_s/${Rlab}.semd.sac${suffix}";
-$Zsyn="$opt_s/${Zlab}.semd.sac${suffix}";
+$Tsyn = "${opt_s}/${Tlab}.semd.sac${suffix}";
+$Rsyn = "${opt_s}/${Rlab}.semd.sac${suffix}";
+$Zsyn = "${opt_s}/${Zlab}.semd.sac${suffix}";
+
+# SYNTHETICS
+$suffix = ".recon.cc"; # USER MUST CHANGE
+$Trecon = "${opt_c}/${Tlab}${suffix}";
+$Rrecon = "${opt_c}/${Rlab}${suffix}";
+$Zrecon = "${opt_c}/${Zlab}${suffix}";
+
#------------
-#print "\n $Tdat $Rdat $Zdat \n $Tsyn $Rsyn $Zsyn \n"; die("testing");
+#print "\n $Tdat $Rdat $Zdat \n $Tsyn $Rsyn $Zsyn \n $Trecon $Rrecon $Zrecon \n"; die("testing");
# get the receiver location, azimuth, distance (from the $Zsyn file)
-(undef,$tmin0,$tmax0,$slon,$slat,$dist,$az,$Tchan)=split(" ",`$saclst b e stlo stla dist az kcmpnm f $Zsyn`);
+(undef,$tmin0,$tmax0,$slon,$slat,$dist,$az,$Tchan) = split(" ",`$saclst b e stlo stla dist az kcmpnm f $Zsyn`);
$edate = sprintf("%4.4i . %2.2i . %2.2i",$year,$month,$day);
$stedep = sprintf("depth = %.2f km",$edep);
$strec = "$sta . $net";
@@ -222,76 +231,25 @@
$xtextpos4 = 0.95*$tmax; # final window measurement label
$ytextpos5 = 1.35*$Srange; # titles
-$ytextpos4 = 1.00*$Srange; # dT label
+$ytextpos4 = 1.03*$Srange; # dT label
$ytextpos3 = 0.90*$Srange; # dlnA label
$ytextpos2 = 0.50*$Srange; # Z, R, T label
$ytextpos1 = 0.15*$Srange; # ymax label
# pstext info for measurement labels
-$textinfo = "-C1p -W255/255/255o,255/255/255";
+$textinfo = "-C2p -W0/255/255o,0/255/255 -G0/0/0";
$proj = "-JX3/1";
$dX = 4;
$dY = 1.5;
-#-----------------------------------
-
-# pssac2 scaling
-# -M vertical scaling in sacfile_unit/MEASURE_UNIT = size<required>
-# size: each trace will normalized to size (in MEASURE_UNIT)
-# scale = (1/size) * [data(max) - data(min)]
-
-# limits for data records
-$maxT = 0; $maxR = 0; $maxZ = 0;
-if($Tflag==1) {(undef,$minT,$maxT)=split(" ",`$saclst depmin depmax f $Tdat`); $maxT=max(abs($minT),abs($maxT));}
-if($Rflag==1) {(undef,$minR,$maxR)=split(" ",`$saclst depmin depmax f $Rdat`); $maxR=max(abs($minR),abs($maxR));}
-if($Zflag==1) {(undef,$minZ,$maxZ)=split(" ",`$saclst depmin depmax f $Zdat`); $maxZ=max(abs($minZ),abs($maxZ));}
-
-# limits for synthetic records
-(undef,$minT2,$maxT2) = split(" ",`$saclst depmin depmax f $Tsyn`);
-$maxT2 = max(abs($minT2),abs($maxT2));
-(undef,$minR2,$maxR2) = split(" ",`$saclst depmin depmax f $Rsyn`);
-$maxR2 = max(abs($minR2),abs($maxR2));
-(undef,$minZ2,$maxZ2) = split(" ",`$saclst depmin depmax f $Zsyn`);
-$maxZ2 = max(abs($minZ2),abs($maxZ2));
-
-# overall max values
-$maxT3 = max($maxT,$maxT2);
-$maxR3 = max($maxR,$maxR2);
-$maxZ3 = max($maxZ,$maxZ2);
-$maxD = max($maxT,max($maxR,$maxZ));
-$maxS = max($maxT2,max($maxR2,$maxZ2));
-$max = max($maxD,$maxS);
-#print "\n Overall max values: $maxD (data), $maxS (syn) $max (overall)\n";
-
-$bounds = "-R$tmin/$tmax/$Smin/$Smax";
-$scale = 1.5; # scaling for plotting seismograms
-$size = 1/$scale;
-
-if ($opt_r) {
- # normalize by the same value
- $sizeT = $size*$maxT/$max;
- $sizeT2 = $size*$maxT2/$max;
- $sizeR = $size*$maxR/$max;
- $sizeR2 = $size*$maxR2/$max;
- $sizeZ = $size*$maxZ/$max;
- $sizeZ2 = $size*$maxZ2/$max;
-
-} else {
- # normalize by different values
- $sizeT = $size*$maxT/$maxT3;
- $sizeT2 = $size*$maxT2/$maxT3;
- $sizeR = $size*$maxR/$maxR3;
- $sizeR2 = $size*$maxR2/$maxR3;
- $sizeZ = $size*$maxZ/$maxZ3;
- $sizeZ2 = $size*$maxZ2/$maxZ3;
-}
-
$tick = "-Ba50f10:\"Time (s)\":/a1f1S";
$tick2 = "-Ba50f10/a1f1S";
-$red_pen = "-W1/250/0/0";
-$black_pen = "-W1/0/0/0";
+$pen = "0.7p";
+$red_pen = "-W${pen},250/0/0";
+$black_pen = "-W${pen},0/0/0";
+$recon_pen = "-W${pen},0/255/255";
my(@Twinb); my(@Twine);
my(@Rwinb); my(@Rwine);
@@ -371,14 +329,97 @@
}
}
+#==========================================
+
+# pssac2 scaling
+# -M vertical scaling in sacfile_unit/MEASURE_UNIT = size<required>
+# size: each trace will normalized to size (in MEASURE_UNIT)
+# scale = (1/size) * [data(max) - data(min)]
+
+# limits for data records
+$maxT = 0; $maxR = 0; $maxZ = 0;
+if($Tflag==1) {(undef,$minT,$maxT)=split(" ",`$saclst depmin depmax f $Tdat`); $maxT=max(abs($minT),abs($maxT));}
+if($Rflag==1) {(undef,$minR,$maxR)=split(" ",`$saclst depmin depmax f $Rdat`); $maxR=max(abs($minR),abs($maxR));}
+if($Zflag==1) {(undef,$minZ,$maxZ)=split(" ",`$saclst depmin depmax f $Zdat`); $maxZ=max(abs($minZ),abs($maxZ));}
+
+# limits for synthetic records
+(undef,$minT2,$maxT2) = split(" ",`$saclst depmin depmax f $Tsyn`);
+$maxT2 = max(abs($minT2),abs($maxT2));
+(undef,$minR2,$maxR2) = split(" ",`$saclst depmin depmax f $Rsyn`);
+$maxR2 = max(abs($minR2),abs($maxR2));
+(undef,$minZ2,$maxZ2) = split(" ",`$saclst depmin depmax f $Zsyn`);
+$maxZ2 = max(abs($minZ2),abs($maxZ2));
+
+# overall max values
+$maxTDS = max($maxT,$maxT2);
+$maxRDS = max($maxR,$maxR2);
+$maxZDS = max($maxZ,$maxZ2);
+$maxD = max($maxT,max($maxR,$maxZ));
+$maxS = max($maxT2,max($maxR2,$maxZ2));
+$max = max($maxD,$maxS);
+#print "\n Overall max values: $maxD (data), $maxS (syn) $max (overall)\n";
+
+$bounds = "-R$tmin/$tmax/$Smin/$Smax";
+$scale = 1.5; # scaling for plotting seismograms
+$size = 1/$scale;
+
+if ($opt_r) {
+ # normalize by the same value (max)
+ $sizeT = $size*$maxT/$max;
+ $sizeT2 = $size*$maxT2/$max;
+ $sizeR = $size*$maxR/$max;
+ $sizeR2 = $size*$maxR2/$max;
+ $sizeZ = $size*$maxZ/$max;
+ $sizeZ2 = $size*$maxZ2/$max;
+
+} else {
+ # normalize by different values
+ $sizeT = $size*$maxT/$maxTDS;
+ $sizeT2 = $size*$maxT2/$maxTDS;
+ $sizeR = $size*$maxR/$maxRDS;
+ $sizeR2 = $size*$maxR2/$maxRDS;
+ $sizeZ = $size*$maxZ/$maxZDS;
+ $sizeZ2 = $size*$maxZ2/$maxZDS;
+}
+
+# limits for reconstructed records
+if ($ntline) {
+ (undef,$minT3,$maxT3) = split(" ",`$saclst depmin depmax f $Trecon`);
+ $maxT3 = max(abs($minT3),abs($maxT3));
+ if ($opt_r) {
+ $sizeT3 = $size*$maxT3/$max;
+ } else {
+ $sizeT3 = $size*$maxT3/$maxTDS;
+ }
+}
+if ($nrline) {
+ (undef,$minR3,$maxR3) = split(" ",`$saclst depmin depmax f $Rrecon`);
+ $maxR3 = max(abs($minR3),abs($maxR3));
+ if ($opt_r) {
+ $sizeR3 = $size*$maxR3/$max;
+ } else {
+ $sizeR3 = $size*$maxR3/$maxRDS;
+ }
+}
+if ($nzline) {
+ (undef,$minZ3,$maxZ3) = split(" ",`$saclst depmin depmax f $Zrecon`);
+ $maxZ3 = max(abs($minZ3),abs($maxZ3));
+ if ($opt_r) {
+ $sizeZ3 = $size*$maxZ3/$max;
+ } else {
+ $sizeZ3 = $size*$maxZ3/$maxZDS;
+ }
+}
+
# strings for plotting -- show the yaxis limits for each trace (based on synthetics)
$strT = sprintf('%.2e',$maxT2/$sizeT2);
$strR = sprintf('%.2e',$maxR2/$sizeR2);
$strZ = sprintf('%.2e',$maxZ2/$sizeZ2);
+#==========================================
+
# TRANSVERSE component: data, synthetics, and windows
-#`pssac2 -X-2 -Y4.5 $proj $Tdat $bounds -Ent-3 -M${sizeT} -P $black_pen -N -K -O $tick >> $ps_file`;
-`pssac2 -X-2 -Y4.5 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} -P $red_pen -N -K -O $tick >> $ps_file`;
+`pssac2 -X-2 -Y4.5 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick >> $ps_file`; # synthetics
if ($ntline) {
for ($i=0; $i<$nT; $i++) {
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
@@ -397,14 +438,14 @@
`echo \"$xtext $ytextpos3 7 0 1 $just $stlabT_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
}
}
-if($Tflag==1) {`pssac2 $proj $Tdat $bounds -Ent-3 -M${sizeT} -P $black_pen -N -K -O >> $ps_file`;}
-`pssac2 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} -P $red_pen -N -K -O >> $ps_file`;
+if($Tflag==1) {`pssac2 $proj $Tdat $bounds -Ent-3 -M${sizeT} $black_pen -N -K -O >> $ps_file`;} # data
+`pssac2 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`; # synthetics
+if($ntline) {`pssac2 $proj $Trecon $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strT\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# RADIAL component: data, synthetics, and windows
-#`pssac2 -Y$dY $proj $Rdat $bounds -Ent-3 -M${sizeR} -P $black_pen -N -K -O $tick2 >> $ps_file`;
-`pssac2 -Y$dY $proj $Rsyn $bounds -Ent-3 -M${sizeR2} -P $red_pen -N -K -O $tick2 >> $ps_file`;
+`pssac2 -Y$dY $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nrline) {
for ($i=0; $i<$nR; $i++) {
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
@@ -419,14 +460,14 @@
`echo \"$xtext $ytextpos3 7 0 1 $just $stlabR_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
}
}
-if($Rflag==1) {`pssac2 $proj $Rdat $bounds -Ent-3 -M${sizeR} -P $black_pen -N -K -O >> $ps_file`;}
-`pssac2 $proj $Rsyn $bounds -Ent-3 -M${sizeR2} -P $red_pen -N -K -O >> $ps_file`;
+if($Rflag==1) {`pssac2 $proj $Rdat $bounds -Ent-3 -M${sizeR} $black_pen -N -K -O >> $ps_file`;}
+`pssac2 $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O >> $ps_file`;
+if($nrline) {`pssac2 $proj $Rrecon $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strR\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# VERTICAL component: data, synthetics, and windows
-#`pssac2 -Y$dY $proj $Zdat $bounds -Ent-3 -M${sizeZ} -P $black_pen -N -K -O $tick2 >> $ps_file`;
-`pssac2 -Y$dY $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} -P $red_pen -N -K -O $tick2 >> $ps_file`;
+`pssac2 -Y$dY $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nzline) {
for ($i=0; $i<$nZ; $i++) {
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
@@ -441,8 +482,9 @@
`echo \"$xtext $ytextpos3 7 0 1 $just $stlabZ_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
}
}
-if($Zflag==1) {`pssac2 $proj $Zdat $bounds -Ent-3 -M${sizeZ} -P $black_pen -N -K -O $tick2 >> $ps_file`;}
-`pssac2 $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} -P $red_pen -N -O -K >> $ps_file`;
+if($Zflag==1) {`pssac2 $proj $Zdat $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $tick2 >> $ps_file`;}
+`pssac2 $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $ps_file`;
+if($nzline) {`pssac2 $proj $Zrecon $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strZ\" | pstext $proj $bounds -N -O -K >> $ps_file`;
@@ -455,8 +497,8 @@
$tick="-Ba50f10:\"Time (s)\":/a1f1S";
$tick2="-Ba50f10/a1f1S";
$proj="-JX3/1";
-$red_pen = "-W1/250/0/0";
-$black_pen = "-W1/0/0/0";
+#$red_pen = "-W1p,250/0/0";
+#$black_pen = "-W1p,0/0/0";
# adjoint sources (always have the BH_ channel, even if the data do not)
$Zadj="ADJOINT_SOURCES/$sta.$net.BHZ.${klab}.adj";
@@ -506,7 +548,7 @@
# VERTICAL component: adjoint source and windows
$bounds="-R$tmin/$tmax/-$maxZ/$maxZ";
-`psbasemap $proj $bounds $tick2 -P -K -O >> $ps_file`;
+`psbasemap $proj $bounds $tick2 -K -O >> $ps_file`;
if ($nzline) {
for($i =0; $i<$nZ; $i++){
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
@@ -518,15 +560,15 @@
}}
if (-f $Zadj) {
$ytextpos2 = 0; $ytextpos1 = -$maxZ + 0.15*(2*$maxZ);
- `psxy $proj $Zadj $bounds -P $black_pen -N -K -O >> $ps_file`;
- if($iboth==1) {`psxy $proj $Zadj2 $bounds -P $red_pen -N -K -O >> $ps_file`;}
+ `psxy $proj $Zadj $bounds $black_pen -N -K -O >> $ps_file`;
+ if($iboth==1) {`psxy $proj $Zadj2 $bounds $red_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strZ\" | pstext $proj $bounds -N -K -O >> $ps_file`;
}
# RADIAL component: adjoint source and windows
$bounds="-R$tmin/$tmax/-$maxR/$maxR";
-`psbasemap -Y-$dY $proj $bounds $tick2 -P -K -O >> $ps_file`;
+`psbasemap -Y-$dY $proj $bounds $tick2 -K -O >> $ps_file`;
if ($nrline) {
for($i =0; $i<$nR; $i++){
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
@@ -538,15 +580,15 @@
}}
if (-f $Radj) {
$ytextpos2 = 0; $ytextpos1 = -$maxR + 0.15*(2*$maxR);
- `psxy $proj $Radj $bounds -P $black_pen -N -K -O >> $ps_file`;
- if($iboth==1) {`psxy $proj $Radj2 $bounds -P $red_pen -N -K -O >> $ps_file`;}
+ `psxy $proj $Radj $bounds $black_pen -N -K -O >> $ps_file`;
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ if($iboth==1) {`psxy $proj $Radj2 $bounds $red_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strR\" | pstext $proj $bounds -N -K -O >> $ps_file`;
}
# TRANSVESE component: adjoint source and windows
$bounds="-R$tmin/$tmax/-$maxT/$maxT";
-`psbasemap -Y-$dY $proj $bounds $tick -P -K -O >> $ps_file`;
+`psbasemap -Y-$dY $proj $bounds $tick -K -O >> $ps_file`;
if ($ntline) {
for($i =0; $i<$nT; $i++){
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
@@ -558,8 +600,8 @@
}}
if (-f $Tadj) {
$ytextpos2 = 0; $ytextpos1 = -$maxT + 0.15*(2*$maxT);
- `psxy $proj $Tadj $bounds -P $black_pen -N -K -O >> $ps_file`;
- if($iboth==1) {`psxy $proj $Tadj2 $bounds -P $red_pen -N -K -O >> $ps_file`;}
+ `psxy $proj $Tadj $bounds $black_pen -N -K -O >> $ps_file`;
+ if($iboth==1) {`psxy $proj $Tadj2 $bounds $red_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -K -O >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strT\" | pstext $proj $bounds -N -K -O >> $ps_file`;
}
Modified: seismo/3D/mt_measure_adj/PLOTS/plot_win_adj_all.pl
===================================================================
--- seismo/3D/mt_measure_adj/PLOTS/plot_win_adj_all.pl 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/PLOTS/plot_win_adj_all.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -14,14 +14,14 @@
sub Usage{
print STDERR <<END;
- plot_all_win_adj.pl -m CMTFILE -a STATION_FILE -b iboth -l tmin/tmax -k iker -d data_dir -s syn_dir -w winfile
+ plot_all_win_adj.pl -m CMTFILE -a STATION_FILE -b iboth -l tmin/tmax -k iker -d data_dir -s syn_dir -c recon_dir -w winfile
END
exit(1);
}
if (@ARGV == 0) { Usage(); }
-if (!getopts('m:a:b:d:k:l:s:w:')) {die('Check input arguments\n');}
+if (!getopts('m:a:b:d:k:l:s:c:w:')) {die('Check input arguments\n');}
if($opt_m) {$cmtfile=$opt_m;}
if(!$opt_b) {$iboth = 0} else {$iboth = $opt_b}
#if($opt_k) {$iker = $opt_k;} # does not work for iker = 1
@@ -30,6 +30,7 @@
if($opt_w) {$winfile=$opt_w;}
if($opt_d) {$data_dir=$opt_d;}
if($opt_s) {$syn_dir=$opt_s;}
+if($opt_c) {$recon_dir=$opt_c;}
$ns = `sed -n 1p $opt_a`;
$nline = 1;
@@ -41,7 +42,7 @@
print "$sta $net\n";
# include the -r command here if you want the traces scaled by the max Z-R-T value
- #print "\n plot_win_adj.pl -m $cmtfile -n $sta/$net -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -w $winfile\n";
- `plot_win_adj.pl -m $cmtfile -n $sta/$net -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -w $winfile`;
+ print "\n plot_win_adj.pl -m $cmtfile -n $sta/$net -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -c $recon_dir -w $winfile\n";
+ `plot_win_adj.pl -m $cmtfile -n $sta/$net -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -c $recon_dir -w $winfile`;
}
Modified: seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl
===================================================================
--- seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -1,7 +1,7 @@
#!/usr/bin/perl -w
#-----------------------------------
-# Carl Tape, 15-Nov-2007
+# Carl Tape, 06-April-2008
# plot_win_stats_all.pl
#
# This script calls plot_win_stats.pl to plot histograms of measurements for each event.
@@ -12,6 +12,8 @@
# EXAMPLE:
# plot_win_stats_all.pl 6/40 m0 0
# plot_win_stats_all.pl 6/40 m0 1 (plot histogram)
+# plot_win_stats_all.pl 2/40 m0 0
+# plot_win_stats_all.pl 2/40 m0 1 (plot histogram)
#
#-----------------------------------
Modified: seismo/3D/mt_measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/mt_measure_adj/mt_measure_adj.f90 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/mt_measure_adj.f90 2008-04-22 18:56:23 UTC (rev 11845)
@@ -13,7 +13,7 @@
character(len=150) :: out_dir,datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
integer :: num_files, num_meas, i, j, k, iker, iker0, ios, is_mtm, is_mtm0, npt1, npt2, npts, nn
- double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src
+ double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, tseis_recon_cc
double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, chi, tt, dtt, df
double precision, dimension(3*N_MEASUREMENT) :: window_chi
double precision :: fend0, fstart0, fend, fstart
@@ -124,6 +124,7 @@
syn(:) = 0.0
adj_syn_all(:) = 0.0
+ recon_cc_all(:) = 0.0
read(11,'(a)',iostat=ios) datafile
if (ios /= 0) stop 'Error reading datafile'
@@ -224,8 +225,9 @@
!call mt_measure_select_0(nlen,dt,i_left,i_right,fstart,fend,use_trace)
! make measurements
+ ! CHT: also compute reconstructed synthetics for CC (and MT, if specified) measurements
call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,dzr20,dzr0,nlen,&
- tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,&
+ tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc,&
i_pmax,i_right0,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
i_right = i_right0
i_left = 1
@@ -249,15 +251,15 @@
endif
print *, ' Status: ', trim(measure_file_prefix), use_trace
- ! CHT: if window is rejected by mt_measure_select, then use cross-correlation
- ! --> What about quality checks for the cross-correlation measurements?
+ ! CHT: If MT measurement window is rejected by mt_measure_select, then use a CC measurement.
+ ! --> But what about quality checks for the cross-correlation measurements?
if(.not. use_trace) then
!stop 'why was this MT measurement rejected?'
print *, 'Reverting from multitaper measurement to cross-correlation measurement'
iker = 2
is_mtm = 3
call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,dzr20,dzr0,nlen,&
- tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,&
+ tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc,&
i_pmax,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA)
use_trace = .true.
endif
@@ -274,9 +276,7 @@
close(71)
endif
- !stop 'testing'
-
- ! compute adjoint sources and misfit function values
+ ! compute adjoint sources and misfit function values and also the CC-reconstructed records
if (use_trace) then
tr_chi = 0. ; am_chi = 0. ! must be initialized
@@ -298,9 +298,12 @@
! combine adjoint sources from different measurement windows
adj_syn_all(:) = adj_syn_all(:) + tr_adj_src(:)
+ ! combine CC-reconstructed records
+ recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + tseis_recon_cc(1:nlen)
+
endif
- ! CHT: (re-)set to multitaper parameters
+ ! CHT: (re-)set to multitaper parameters, if originally specified
if (iker0 == 1) then
iker = iker0
is_mtm = is_mtm0
@@ -308,6 +311,9 @@
enddo ! nmeas
+ !----------------------------
+ ! write out the adjoint source
+
! filter the adjoint source (higher frequencies could enter from the tapering operations)
! (See also time_window in mt_adj.f90, which tapers the windows.)
call xapiir(adj_syn_all,npts,'BU',TRBDNDW,APARM,IORD,'BP',fstart0,fend0,dt,PASSES)
@@ -323,8 +329,18 @@
! we can have choices of outputing the adjoint source as ASCII or SAC format
call dwrite_ascfile_f(trim(adj_file_prefix)//'.adj',tt,dtt,nn,adj_syn_all)
- ! call dwrite_sacfile_f(trim(file_prefix)//'.adj',t0,dt,npts,adj_syn_all)
+ !call dwrite_sacfile_f(trim(datafile),trim(adj_file_prefix)//'.adj',tt,nn,adj_syn_all)
+ !----------------------------
+ ! write out the CC-reconstructed data from synthetics
+
+ ! cut and interpolate, then write ASCII file
+ !call interpolate_syn(recon_cc_all,t0,dt,npts,tt,dtt,nn)
+ !call dwrite_ascfile_f(trim(file_prefix2)//'.recon.cc',tt,dtt,nn,recon_cc_all)
+
+ ! write SAC file
+ call dwrite_sacfile_f(trim(datafile),trim(file_prefix2)//'.recon.cc',t0,npts,recon_cc_all)
+
enddo ! npairs
close(11) ! read: MEASUREMENT.WINDOWS
Modified: seismo/3D/mt_measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/mt_measure_adj/mt_sub.f90 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/mt_sub.f90 2008-04-22 18:56:23 UTC (rev 11845)
@@ -40,7 +40,7 @@
! =================================================================================================
subroutine mt_measure(is_mtm,iker,datafile,file_prefix,data,syn,t0,dt,npts,tstart,tend, &
- istart,dzr20,dzr0,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,&
+ istart,dzr20,dzr0,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc, &
i_pmax,i_right,trans_w,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
integer, intent(in) :: is_mtm, iker, npts
@@ -140,8 +140,8 @@
dat_displ(1:nlen) = dzr20(1:nlen)
do i = 2, nlen-1
- syn_veloc(i) = (syn_displ(i+1) - syn_displ(i-1)) / (2.*dt)
- dat_veloc(i) = (dat_displ(i+1) - dat_displ(i-1)) / (2.*dt)
+ syn_veloc(i) = (syn_displ(i+1) - syn_displ(i-1)) / (2.0*dt)
+ dat_veloc(i) = (dat_displ(i+1) - dat_displ(i-1)) / (2.0*dt)
enddo
syn_veloc(1) = (syn_displ(2) - syn_displ(1)) / dt
syn_veloc(nlen) = (syn_displ(nlen) - syn_displ(nlen-1)) /dt
@@ -149,10 +149,10 @@
dat_veloc(nlen) = (dat_displ(nlen) - dat_displ(nlen-1)) /dt
do i = 2, nlen-1
- syn_accel(i) = (syn_veloc(i+1) - syn_veloc(i-1)) / (2.*dt)
+ syn_accel(i) = (syn_veloc(i+1) - syn_veloc(i-1)) / (2.0*dt)
enddo
syn_accel(1) = (syn_veloc(2) - syn_veloc(1)) / dt
- syn_accel(nlen) = (syn_veloc(nlen) - syn_veloc(nlen-1)) /dt
+ syn_accel(nlen) = (syn_veloc(nlen) - syn_veloc(nlen-1)) / dt
! note that the records are UN-SHIFTED
dlnA = sqrt( (dt * sum( dat_displ(1:nlen) * dat_displ(1:nlen) )) &
@@ -597,7 +597,7 @@
do i = i_left, i_right ! CHT: 1 --> i_left
! type of filter in the freq domain
!w_taper(i) = 1. ! boxcar
- !w_taper(i) = 1. - (2./nw)**2 * ((i-1) - nw/2.)**2 ! welch
+ !w_taper(i) = 1. - (2.0/nw)**2 * ((i-1) - nw/2.0)**2 ! welch
w_taper(i) = 1. - cos(PI*(i-i_left)/(i_right-i_left))**ipwr_w ! cosine
enddo
@@ -650,7 +650,7 @@
ipwr_t = 10
do i = 1,nlen
fac = 1. ! boxcar window
- !fac = 1 - sfac2*((i-1) - dble(nlen1)/2.)**2 ! welch window
+ !fac = 1 - sfac2*((i-1) - dble(nlen1)/2.0)**2 ! welch window
!fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t ! cosine window
time_window(i) = fac
enddo
@@ -661,7 +661,7 @@
! compute synthetic velocity
do i = 2, nlen-1
- vzr0(i) = (dzr0(i+1) - dzr0(i-1)) / (2.*dt)
+ vzr0(i) = (dzr0(i+1) - dzr0(i-1)) / (2.0*dt)
enddo
vzr0(1) = (dzr0(2) - dzr0(1)) / dt
vzr0(nlen) = (dzr0(nlen) - dzr0(nlen-1)) / dt
@@ -701,7 +701,7 @@
! note: cannot directly apply tapers on vzr0, which is not the same as vzr
do i = 2, nlen-1
- vzr(i) = (dzr(i+1) - dzr(i-1)) / (2.*dt)
+ vzr(i) = (dzr(i+1) - dzr(i-1)) / (2.0*dt)
enddo
vzr(1) = (dzr(2) - dzr(1)) / dt
vzr(nlen) = (dzr(nlen) - dzr(nlen-1)) /dt
@@ -807,7 +807,8 @@
endif
enddo
- ! CHT: ompute misfit function value and measurement value
+ ! CHT: compute misfit function value and measurement value
+ ! Note: The taper functions for MT may include error estimates.
! 1: waveform difference
! 2: multitaper, TT
! 3: multitaper, dlnA
@@ -892,12 +893,12 @@
endif
!write(*,'(a8,4f12.6)') 'fstart :', fstart, NCYCLE_IN_WINDOW/(Wlen), NCYCLE_IN_WINDOW, Wlen
- !write(*,'(a8,4f12.6)') 'fend :', fend, 1./(2.*dt), dt
+ !write(*,'(a8,4f12.6)') 'fend :', fend, 1./(2.0*dt), dt
! DECREASE the frequency range of the measurement (and adjoint source)
! note NCYCLE_IN_WINDOW and window length
fstart = max(fstart, NCYCLE_IN_WINDOW/Wlen)
- fend = min(fend, 1./(2.*dt))
+ fend = min(fend, 1./(2.0*dt))
if (fstart >= fend) then
print *, 'rejecting trace for frequency range (NCYCLE_IN_WINDOW):'
print *, ' fstart >= fend : ', fstart, fend
@@ -1069,7 +1070,7 @@
!!$
!!$ ! compute synthetic velocity (unshifted)
!!$ do i = 2, nlen-1
-!!$ vzr(i) = (dzr(i+1) - dzr(i-1)) / (2.*dt)
+!!$ vzr(i) = (dzr(i+1) - dzr(i-1)) / (2.0*dt)
!!$ enddo
!!$ vzr(1) = (dzr(2) - dzr(1)) / dt
!!$ vzr(nlen) = (dzr(nlen) - dzr(nlen-1)) / dt
@@ -1121,7 +1122,7 @@
! compute synthetic velocity (unshifted)
do i = 2, nlen-1
- syn_v_dt_dlnA(i) = (syn_d_dt_dlnA(i+1) - syn_d_dt_dlnA(i-1)) / (2.*dt)
+ syn_v_dt_dlnA(i) = (syn_d_dt_dlnA(i+1) - syn_d_dt_dlnA(i-1)) / (2.0*dt)
enddo
syn_v_dt_dlnA(1) = (syn_d_dt_dlnA(2) - syn_d_dt_dlnA(1)) / dt
syn_v_dt_dlnA(nlen) = (syn_d_dt_dlnA(nlen) - syn_d_dt_dlnA(nlen-1)) / dt
@@ -1250,9 +1251,9 @@
if (i > 1 .and. i < i_right) then
! check the smoothness (2nd-order derivative) by 2*pi changes
- smth = phi_wt(i+1) + phi_wt(i-1) - 2. * phi_wt(i)
- smth1 = (phi_wt(i+1) + TWOPI) + phi_wt(i-1) - 2. * phi_wt(i)
- smth2 = (phi_wt(i+1) - TWOPI) + phi_wt(i-1) - 2. * phi_wt(i)
+ smth = phi_wt(i+1) + phi_wt(i-1) - 2.0 * phi_wt(i)
+ smth1 = (phi_wt(i+1) + TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
+ smth2 = (phi_wt(i+1) - TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
if(abs(smth1).lt.abs(smth).and.abs(smth1).lt.abs(smth2).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
if (DISPLAY_DETAILS) print *, 'phase correction : 2 pi', sngl(fr(i)), sngl(phi_wt(i) - phi_wt(i+1))
do j = i+1, i_right
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m 2008-04-22 18:56:23 UTC (rev 11845)
@@ -8,7 +8,7 @@
%
% /cig/seismo/3D/mt_measure_adj_work/scripts_tomo/matlab/
%
-% calls readCMT.m, readSTATIONS.m
+% calls xxx
% called by xxx
%
@@ -16,7 +16,6 @@
clear
close all
-clc
format compact
dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
@@ -29,17 +28,31 @@
comps = {'BHZ','BHR','BHT'};
% regardless of the component label for the DATA, the measurement code
% defaults so that the first two letters are always BH
-
+
ncomp = length(comps);
nper = length(Tvec);
+idatacov = input(' Enter idatacov (1, 2, 3) : ');
+ % 1: weight by N
+ % 2: weight by E N_e in blocks that are N_e each in length
+ % 3: weight by 1
+ % -- N_e is the number of windows for event e
+ % -- E is the number of events
+ % -- N is the total number of windows: N = sum_e N_e
+stdcov = {'event','window','none'};
+
+iwrite = 1;
+
%-------------------------
% read in list of event IDs and sources
+% load the event IDs corresponding to the kernels
% load these as strings, since event IDs could have letters
-kerlist = 'kernels_done_mod';
-eids = textread(kerlist,'%s');
+stker0 = '/net/sierra/raid1/carltape/results/KERNELS/';
+%eids = textread([stker0 'kernels_test'],'%s');
+eids = textread([stker0 'kernels_done_mod'],'%s');
nevent = length(eids);
+eids_num = str2num(char(eids)); % numeric version
dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_5/';
%eid_file = [dir_source 'SOCAL_FINAL_CMT_v5_eid'];
@@ -85,62 +98,191 @@
%windows_array = zeros(nevent,nrec,ncomp,nper); % number of windows
-k2 = 0;
-meas_array = [];
-for tt = 1:1
-%for tt = 1:nper
- Tper = Tvec(tt);
- sTper = sprintf('%2.2i',Tper);
-
- %for ii = 1:20
- for ii = 1:nevent
-
- disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
- dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
- mfile = [eids{ii} '_T' sTper '_window_chi'];
- filename = [dir1 mfile];
-
- if exist(filename,'file')
- [strec0,cmp,iwin,iker,...
- chi1,chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
- meas1,measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
- sigma1,sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,tr_chi,am_chi] = read_window_chi(filename);
- nwin = length(strec0)
+iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
- % assign the string receiver an integer index
- index_rec = NaN * ones(nwin,1);
- for irec = 1:nrec
- ind_rec = find( strcmp(strec(irec), strec0) == 1 );
- index_rec(ind_rec) = irec;
- end
+if iread == 1
- % assign the string component an integer index
- index_comp = NaN * ones(nwin,1);
- for icomp = 1:ncomp
- ind_cmp = find( strcmp(comps(icomp), cmp) == 1 );
- index_comp(ind_cmp) = icomp;
+ k2 = 0;
+ meas_array = [];
+ %for tt = 1:1
+ for tt = 1:nper
+ Tper = Tvec(tt);
+ sTper = sprintf('%2.2i',Tper);
+ disp('-----------------------------------------');
+
+ %for ii = 1:20
+ for ii = 1:nevent
+
+ disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
+ dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
+ mfile = [eids{ii} '_T' sTper '_window_chi'];
+ filename = [dir1 mfile];
+
+ if exist(filename,'file')
+ [strec0,cmp,iwin,iker,...
+ chi1,chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
+ meas1,measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
+ sigma1,sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,tr_chi,am_chi] = read_window_chi(filename);
+ nwin = length(strec0);
+
+ % assign the string receiver an integer index
+ index_rec = NaN * ones(nwin,1);
+ for irec = 1:nrec
+ ind_rec = find( strcmp(strec(irec), strec0) == 1 );
+ index_rec(ind_rec) = irec;
+ end
+
+ % assign the string component an integer index
+ index_comp = NaN * ones(nwin,1);
+ for icomp = 1:ncomp
+ ind_cmp = find( strcmp(comps(icomp), cmp) == 1 );
+ index_comp(ind_cmp) = icomp;
+ end
+
+ % measurement indices
+ k1 = k2+1;
+ k2 = k1+nwin-1;
+ kinds = [k1:k2]';
+
+ meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_rec index_comp iwin ...
+ iker chi1 measCC_dT sigmaCC_dT measMT_dT sigmaMT_dT tr_chi];
+ % 1 kinds
+ % 2 index_T
+ % 3 index_event
+ % 4 index_rec
+ % 5 index_comp
+ % 6 iwin
+ % 7 iker
+ % 8 chi1
+ % 9 measCC_dT
+ % 10 sigmaCC_dT
+ % 11 measMT_dT
+ % 12 sigmaMT_dT
+ % 13 tr_chi
end
- % measurement indices
- k1 = k2+1;
- k2 = k1+nwin-1;
- kinds = [k1:k2]';
-
- meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_rec index_comp iwin ...
- iker chi1 measCC_dT sigmaCC_dT measMT_dT sigmaMT_dT tr_chi];
end
+ end
+
+ save('meas_array.mat','meas_array');
+else
+ load('meas_array.mat');
+end
+
+whos meas_array
+
+% total number of windows
+N = length(meas_array);
+
+% check some of the output window rows
+if 0==1
+ itestvec = [1 2];
+ for ii = 1:length(itestvec)
+ disp('====>');
+ itest = itestvec(ii);
+ meas_array(itest,:)
+ display_meas(meas_array(itest,:),Tvec,eids,strec,comps);
+ meas_array(itest,13)
+ 0.5 * meas_array(itest,9)^2 / meas_array(itest,10)^2
+ 0.5 * meas_array(itest,11)^2 / meas_array(itest,10)^2
end
end
+% tally windows by event
+nwin_all = zeros(nevent,2);
+for ii = 1:nevent
+ for tt = 1:nper
+ imatch = find( and( meas_array(:,3)==ii, meas_array(:,2)==tt) );
+ nwin_all(ii,tt) = length(imatch);
+ end
+end
+nwin_all_event = sum(nwin_all, 2);
+nwin_all_per = sum(nwin_all, 1);
+nwin_tot = [nwin_all nwin_all_event ; nwin_all_per sum(sum(nwin_all)) ];
+dtemp = [ [eids_num ; 0] nwin_tot];
+%sortrows(dtemp,3)
+
+% construct data covariance normalization terms -- the sigma estimates have
+% already been folded into the misfit function value 0.5 DT^2 / sigma^2
+Ns = zeros(nevent,1);
+dcov_fac = zeros(N,1);
+dcov_fac_e = zeros(nevent,1);
+for ii = 1:nevent
+ imatch = find( meas_array(:,3)==ii );
+ Ns(ii) = length(imatch);
+ if idatacov == 1
+ dcov_fac_e(ii) = N;
+ elseif idatacov == 2
+ dcov_fac_e(ii) = Ns(ii) * nevent;
+ elseif idatacov == 3
+ dcov_fac_e(ii) = 1;
+ end
+ dcov_fac(imatch) = dcov_fac_e(ii);
+end
+
+% check number of windows
+if sum(Ns) ~= sum(sum(nwin_all)), error('inconsistent values for total windows'); end
+
+% compute data vector for subspace method
+% factor of 2 cancels the 0.5 used in computing the misfit (mt_measure_adj.f90)
+dnorm_sq = zeros(nevent,1);
+for ii = 1:nevent
+ imatch = find( meas_array(:,3)==ii );
+ dnorm_sq(ii) = sum(2 * meas_array(imatch,13) ./ dcov_fac(imatch) );
+end
+
+% check number of windows
+if sum(Ns) ~= sum(sum(nwin_all)), error('inconsistent values for total windows'); end
+
+% weight vector for the subspace method
+dnorm = sqrt(dnorm_sq);
+ws = 1 ./ dnorm;
+
+% total data misfit -- this is like a PER WINDOW measure of misfit
+dmisfit = sum( dnorm_sq );
+disp(' total data misfit = d^T Cd^-1 d and its square-root:');
+disp(dmisfit);
+disp(sqrt(dmisfit));
+
+if iwrite == 1
+ % display sorted from greatest to least norm
+ data_norms = [eids_num Ns dnorm_sq dnorm ws];
+ dsort = sortrows(data_norms,5);
+ fid = fopen(['data_norms_sort_dcov_' stdcov{idatacov} '.txt'],'w');
+ fprintf(fid,'%12s%8s%12s%12s%12s\n','eid','Nwin','dnorm2','dnorm','weight');
+ fprintf(fid,'%12s%8i%12.4f%12s%12s\n','TOTAL',sum(Ns),sum(dnorm_sq),'--','--');
+ for ii = 1:nevent
+ fprintf(fid,'%12i%8i%12.4f%12.4f%12.4f\n',dsort(ii,:));
+ end
+ fprintf(fid,'%12s%8i%12.4f%12s%12s\n','TOTAL',sum(Ns),sum(dnorm_sq),'--','--');
+ fclose(fid);
+
+ % save variables (for wave2d_subspace_3D.m)
+ save(['data_norms_' stdcov{idatacov}],'idatacov','dcov_fac_e',...
+ 'eids','eids_num','Ns','dnorm','ws','dmisfit','N');
+end
+
+break
+
+%-------------------------------
+
+% % find records that meet particular criteria
+% DT_MIN = 5;
+% DT_SIGMA_MIN = 1;
+% icheck = find( and( meas_array(:,9) >= DT_MIN, meas_array(:,10) <= DT_SIGMA_MIN) )
+% [junk, isort] = sortrows( meas_array(icheck,:), -9 )
+% meas_disp = meas_array(icheck(isort),:);
+% display_meas(meas_disp,Tvec,eids,strec,comps);
+
% find records that meet particular criteria
-DT_MIN = 5;
-DT_SIGMA_MIN = 1;
-icheck = find( and( meas_array(:,9) >= DT_MIN, meas_array(:,10) <= DT_SIGMA_MIN) )
-[junk, isort] = sortrows( meas_array(icheck,:), -9 )
+CHI_MIN = 100;
+DT_SIGMA_MIN = 0.19;
+icheck = find( and( meas_array(:,13) >= CHI_MIN, meas_array(:,10) < DT_SIGMA_MIN) );
+[junk, isort] = sortrows( meas_array(icheck,:), -13 );
meas_disp = meas_array(icheck(isort),:);
-
display_meas(meas_disp,Tvec,eids,strec,comps);
+ebads = unique( eids_num(meas_array(icheck,3)) )
break
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/display_meas.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/display_meas.m 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/display_meas.m 2008-04-22 18:56:23 UTC (rev 11845)
@@ -47,9 +47,9 @@
disp('----------------');
disp('INDEX (window, period, event, receiver, component)');
for ii = 1:m
- disp(sprintf('%10s --> %6.1f %7s %5s %3i -- %6.2f +- %6.2f -- %i',char(eids{index_event(ii)}),...
+ disp(sprintf('%10s --> %6.1f %7s %5s %3i -- %6.2f +- %6.2f -- chi %6.2f -- %i',char(eids{index_event(ii)}),...
Tvec(index_per(ii)),char(strec{index_rec(ii)}),char(comps{index_comp(ii)}),isub_win(ii),...
- measCC_dT(ii),sigmaCC_dT(ii),index_win(ii)));
+ measCC_dT(ii),sigmaCC_dT(ii),tr_chi(ii),index_win(ii) ));
end
disp('INDEX (window, period, event, receiver, component)');
disp('====================');
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/kernels_done_mod
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/kernels_done_mod 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/kernels_done_mod 2008-04-22 18:56:23 UTC (rev 11845)
@@ -6,7 +6,6 @@
10100053
10063349
14095540
-14094992
9644345
9674213
9686565
@@ -38,7 +37,6 @@
9171679
14186612
14138080
-3321595
9122706
9114812
9113909
@@ -56,10 +54,10 @@
14077668
9753489
9941081
-9755013
9753485
9753497
9753949
+9173374
13938812
13936432
13936812
@@ -69,8 +67,11 @@
9775765
9117942
3320736
+9111353
13945908
+9114763
9703873
+9173365
14116972
13692644
14116920
@@ -88,10 +89,14 @@
9818433
9735129
9644101
+14239184
10148421
10148369
9655209
+10187953
14118096
+9627557
+12659440
14151344
9718013
10223765
@@ -102,26 +107,28 @@
14255632
14133048
9722633
-9722529
9774569
9826789
14183744
-10230869
+10226877
14178184
14179736
14178248
14178188
14178236
-14179292
10207681
+14263544
14263768
-13966672
+14263712
+13966396
9154092
14137160
14181056
10215753
10186185
-14072464
+12456160
+9146641
3321590
10964587
10992159
+14178212
Added: seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl (rev 0)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -0,0 +1,47 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 22-April-2008
+# make_dirs.pl
+#
+# This script empties the files within the mu directories.
+# These files are primarily the mu vectors used in the subspace method.
+#
+# This should be run prior to running subspace_specfem.m.
+#
+#-----------------------------------
+
+# possible labels for choice of data covariance
+ at dcov_tags = ("event","window","none");
+
+# possible labels for choice of kernels
+ at kers = ("mu_kernel","kappa_kernel","mu_kernel_smooth","kappa_kernel_smooth");
+
+$ndcov = @dcov_tags;
+$nker = @kers;
+
+$dir0 = "mu_dirs";
+if(not -e $dir0) {die("dir dir0 $dir0 does not exist\n");}
+
+for ($i = 1; $i <= $nker; $i = $i+1) {
+
+ for ($j = 1; $j <= $ndcov; $j = $j+1) {
+
+ $ker = $kers[$i-1];
+ $dcov = $dcov_tags[$j-1];
+ $dir = "$dir0/mu_${ker}_${dcov}";
+
+ print "$dir \n";
+
+ if(-e $dir) {
+ print "--> dir exists -- now deleting and remaking\n";
+ `rm -rf $dir`;
+ } else {
+ print "--> dir does not exist -- now making\n";
+ }
+ `mkdir -p $dir`;
+
+ }
+}
+
+#================================================
Property changes on: seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m (rev 0)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m 2008-04-22 18:56:23 UTC (rev 11845)
@@ -0,0 +1,369 @@
+%
+% subspace_specfem.m
+% CARL TAPE, 16-April-2008
+% printed xxx
+%
+% This program assumes you have FIRST run compute_misfit.m, and also have
+% generated the corresponding Hessian matrix to load here.
+%
+% See also wave2d_subspace.m
+%
+% calls xxx
+% called by xxx
+%
+
+clear
+close all
+format short
+format compact
+
+ax1 = [-121 -114 31 37];
+stfm = '%4.4i';
+
+xlab1 = 'p, singular value index';
+xlab2 = 'p, singular value truncation index';
+ylab1 = 'log10( singular value )';
+ylab2 = 'log10( misfit : dot[ d - G*dm(p), d - G*dm(p) ] )';
+
+stker0 = '/net/sierra/raid1/carltape/results/KERNELS/';
+
+iwrite = 1;
+
+%----------------------------------------------
+% check the slick operation used in subspace_update.f90
+
+if 0==1
+ npar = 10;
+ nsrc = 4;
+ G = rand(npar,nsrc);
+ d = rand(nsrc,1);
+ G*d
+
+ x = zeros(npar,1);
+ for isrc=1:nsrc
+ for ipar=1:npar
+ x(ipar) = x(ipar) + G(ipar,isrc) * d(isrc);
+ end
+ end
+end
+
+%----------------------------------------------
+
+% load the data and event weights (compute_misfit.m)
+idatacov = input(' Enter idatacov (1, 2, 3) : ');
+stdcov = {'event','window','none'};
+load(['data_norms_' stdcov{idatacov}]);
+nsrc = length(eids);
+
+% construct the matrix of weight factors and data-covariance normalization
+Ws = zeros(nsrc,nsrc);
+Cd_fac = zeros(nsrc,nsrc);
+for i = 1:nsrc
+ for j = 1:nsrc
+ Ws(i,j) = ws(i) * ws(j);
+ Cd_fac(i,j) = dcov_fac_e(i) * dcov_fac_e(j);
+ i_ind(i,j) = i;
+ j_ind(i,j) = j;
+ end
+end
+figure; pcolor(i_ind,j_ind,Ws); shading flat;
+xlabel('Row index'); ylabel('Column index');
+title('Symmetric matrix of data weight terms');
+%caxis([0 200]);
+colorbar; axis equal; axis tight;
+
+% load the event IDs corresponding to the kernels
+%eids = load([stker0 'kernels_done_mod']);
+%neid = length(eids);
+
+kmin = 3;
+kmax = 3;
+nker = kmax - kmin + 1;
+stkers = {'mu_kernel','mu_kernel_smooth','kappa_kernel_smooth'};
+stkerlabs = {'unsmoothed beta kernel','smoothed beta kernel','smoothed bulk kernel'};
+
+stwts = {'with data weights'};
+%stwts = {'with data weights','with NO data weights'};
+nker = length(stkers);
+%stcs = {'b.','r.'};
+iw = 1;
+
+%for iw = 1:length(stwts)
+
+for iker = kmin:kmax
+ stker = stkers{iker};
+ stkerlab = stkerlabs{iker};
+
+ hfile = sprintf([stker0 'hessian_index_' stker '_%3.3i'],nsrc);
+ if ~exist(hfile,'file'), error([hfile ' does not exist']); end
+ [kall,iall,jall,Hij] = textread(hfile,'%f%f%f%f');
+ n = length(kall);
+ nhess = (-1 + sqrt(1+8*n))/2; % quadratic formula
+
+ if nhess ~= nsrc, error('number of events must be consistent'); end
+ pinds = [1:nsrc]';
+
+ % construct the Hessian matrix
+ H0 = zeros(nsrc,nsrc);
+ for i = 1:nsrc
+ for j = i:nsrc
+ ih = find( and( i == iall, j == jall) );
+ H0(i,j) = Hij(ih);
+ end
+ end
+ H0 = H0 + H0' - diag(diag(H0)); % fill the lower triangle
+ H = H0;
+
+ % include the weights determined from the data
+ if iw == 1
+ H = H .* Ws;
+ end
+
+ % When we computed the kernels, we did not include the
+ % normalization factor based on the number of windows picked
+ % that is included within the data covariance (and used in computing Ws).
+ % Thus, we compute the MATRIX Cd_fac above to account for this.
+ H = H ./ Cd_fac;
+
+ % put the matrix onto order-1 values for computation purposes
+ %Hnorm = norm(H); H = H / norm(H);
+
+ %H = H + eye(nsrc);
+ % mutemp = inv(H + eye(nsrc)) * dnorm;
+
+ dH = diag(H);
+ trH = sum(dH);
+ trHH = sum(diag(H'*H));
+
+ disp(' properties of Hessian (min, median, mean(abs), max, std):');
+ stH1 = sprintf('min %.1e, median %.1e, mean(abs) %.1e, max %.1e, std %.1e, sum %.1e',...
+ min(H(:)), median(H(:)), mean(abs(H(:))), max(H(:)), std(H(:)), sum(H(:)) );
+ stH2 = sprintf('DIAGONAL : min %.1e, median %.1e, mean(abs) %.1e, max %.1e, std %.1e, sum %.1e',...
+ min(dH), median(dH), mean(abs(dH)), max(dH), std(dH), sum(dH) );
+ disp(stH1);
+ disp(stH2);
+
+ % plot H
+ figure; pcolor(i_ind,j_ind,H); shading flat;
+ xlabel('Row index'); ylabel('Column index');
+ title({['Hessian (symmetric matrix) for ' stkerlab],stH1,stH2});
+ %caxis([0 1e-10])
+ map = colormap('gray'); colormap(flipud(map));
+ colorbar; axis equal; axis tight;
+ fontsize(10), orient tall, wysiwyg
+
+ % sort the diagonal values of H
+ Hdiag = diag(H);
+ %Hdiag = diag(H0);
+ [jk,isort] = sortrows([pinds eids_num Hdiag],-3);
+ npick = min([20 nsrc]);
+ npick = nsrc;
+ disp(' '); disp([' first ' num2str(npick) ' sorted diagonal entries of H: ' ...
+ stkerlab ' ' stwts{iw}]);
+ for ii = 1:npick
+ jj = isort(ii);
+ disp(sprintf('%3i %10i (%3i), Ne = %5i, Hkk = %10.4e',...
+ ii,eids_num(jj),pinds(jj),Ns(jj),Hdiag(jj)));
+ end
+
+ figure(20+iker); hold on;
+ plot(pinds,log10(sort(diag(H),'descend')),'ro');
+ %plot(pinds,log10(sort(diag(H'*H),'descend')),'b.');
+ plot([1 nsrc],log10(trH)*[1 1],'r--');
+ %plot([1 nsrc],log10(trHH)*[1 1],'b--');
+ %legend('diagonal element of H','diagonal element of HtH','trace (H)','trace (HtH)');
+ legend('diagonal element of H','trace (H)');
+ xlabel('sorted event index'); ylabel(' log10 H-kk');
+ title(sprintf(' trace of H = %6.2e --- trace of HtH = %6.2e',trH,trHH));
+ xlim([1 nsrc]); grid on;
+ fontsize(11), orient tall, wysiwyg
+
+ %----------------------------------------------
+
+ itsvd = 0;
+
+ if itsvd == 0
+
+ % regularization choices
+ numlam = 100;
+ if 0==1
+ minlampwr = -2.25; maxlampwr = 1;
+ lampwr = linspace(minlampwr,maxlampwr,numlam);
+ lamvec = 10.^lampwr * sum(dH); % based on trace of H
+ %lamvec = 10.^lampwr * sqrt(sum(diag(H'*H))); % based on trace of H'*H
+ else
+ minlampwr = -23; maxlampwr = -3; % 1
+ minlampwr = -12; maxlampwr = -5; % 3
+ lampwr = linspace(minlampwr,maxlampwr,numlam);
+ lamvec = 10.^lampwr;
+ end
+
+ [f_h, rss, mss, Gvec, Fvec, dof, kap, iL, iGCV, iOCV] = ridge_carl(dnorm,H,lamvec);
+ title(stkerlab);
+
+ mu_all = zeros(nsrc,numlam);
+ mu_norm = zeros(numlam,1);
+ mu_res = zeros(numlam,1);
+ for ip = 1:numlam
+ lam = lamvec(ip);
+ %if ip==iGCV, disp('iGCV -- GCV lambda pick'); end
+ %if ip==iOCV, disp('iOCV -- OCV lambda pick'); end
+ %if ip==iL, disp('iL -- L-curve lambda pick'); end
+ disp(sprintf('%i out of %i, lam = %.2e',ip,numlam,lam));
+
+ mu_all(:,ip) = inv(H + lam*eye(nsrc))*dnorm;
+ %mu_all(:,ip) = inv(H'*H + lam^2*eye(nsrc))*H'*dnorm;
+
+ % f = inv(G'*diag(Wvec)*G + lam0^2*eye(ngrid*ndim))*G'*diag(Wvec)*d;
+ mu_norm(ip) = norm(mu_all(:,ip));
+ mu_res(ip) = norm( dnorm - H*mu_all(:,ip) );
+ end
+
+ ltemp = lamvec(iGCV)
+ mtemp = mu_all(:,iGCV);
+ mtemp_norm = mu_norm(iGCV)
+
+ figure; nr=2; nc=2;
+
+ subplot(nr,nc,1); hold on; grid on;
+ plot( log10(lamvec), log10(mu_norm),'b.','markersize',10);
+ plot( log10(lamvec(iGCV)), log10(mu_norm(iGCV)),'ro','markersize',10);
+ xlabel('log10( lambda )'); ylabel('log10( norm of mu )');
+
+ subplot(nr,nc,2); hold on; grid on;
+ plot( log10( lamvec ), mu_res, '.');
+ plot( log10( lamvec(iGCV) ), mu_res(iGCV),'ro','markersize',10);
+ xlabel('log10( lambda )'); ylabel('dnorm residual');
+
+ subplot(nr,nc,3); plot( dnorm, H*mtemp, '.');
+ xlabel('dnorm obs'); ylabel('dnorm pred'); grid on; axis equal;
+
+ subplot(nr,nc,4); plot( 100*(dnorm - H*mtemp)./dnorm ,'.'); grid on;
+ xlabel('event index'); ylabel('dnorm residual = 100(pred - obs)/obs');
+
+ fontsize(10), orient tall, wysiwyg
+
+% iCV = input(' Enter 1 for OCV or 0 for GCV : ');
+% if iCV == 1
+% lam = lamvec(iOCV);
+% elseif iCV == 0
+% lam = lamvec(iGCV);
+% else
+% error(' iCV must be 1 or 0');
+% end
+% mu = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm; % KEY vector
+% stp = sprintf('lambda = %.3f',lam);
+
+ % write mu vectors to file
+ if iwrite == 1
+ odir = ['mu_dirs/mu_' stkers{iker} '_' stdcov{idatacov} ];
+ for ip = 1:numlam
+ lam = lamvec(ip);
+ mutemp = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm;
+ filename = sprintf([odir '_p%3.3i'],ip);
+ fid = fopen([odir '/' filename],'w');
+ for ii = 1:nsrc
+ fprintf(fid,'%24.12e\n',mutemp(ii));
+ end
+ fclose(fid);
+ end
+
+ fid = fopen([odir '/lambda_' stdcov{idatacov}],'w');
+ for ip = 1:numlam
+ lam = lamvec(ip);
+ fprintf(fid,'%24.12e%24.12e%6i\n',lamvec(ip),1/lamvec(ip),ip);
+ end
+ fclose(fid);
+ end
+
+ else
+
+ % analyze the singular values of H
+ [U,S,V] = svd(H);
+ s = diag(S);
+
+ pmin = min([13 nsrc]);
+ ifit = [(pmin+1):60];
+ %ifit = pinds;
+ [xf,yf,mf,bf,rms] = linefit(pinds(ifit), log10(s(ifit)));
+
+ %subplot(nr,nc,iker);
+ figure; hold on;
+ plot(xf,yf,'g','linewidth',6);
+ plot(pinds,log10(s),'b.','markersize',10);
+ plot(pmin*[1 1],[min(log10(s)) max(log10(s))],'k','linewidth',2)
+ grid on; xlabel(xlab1); ylabel(ylab1);
+ title(['singular value spectra for Hessian : ' stker]);
+ fontsize(10), orient tall, wysiwyg
+
+ %-------------------------
+
+ %tpinds = [1:10:100];
+ tpinds = pinds;
+ [mu_all,rss,f_r_ss] = tsvd(dnorm,H,tpinds);
+
+ % norms of mu vectors
+ mu_norm = zeros(nsrc,1);
+ for ip = 1:nsrc
+ mu_norm(ip) = norm(mu_all(:,ip));
+ end
+
+ figure;
+ subplot(2,2,1); hold on; grid on;
+ plot(pinds,log10(s),'b.','markersize',10);
+ xlabel(xlab1); ylabel(ylab1);
+ subplot(2,2,2); hold on; grid on;
+ plot( tpinds, log10(rss),'b.','markersize',10)
+ xlabel(xlab2); ylabel(ylab2);
+ subplot(2,2,3); hold on; grid on;
+ plot( tpinds, log10(mu_norm),'b.','markersize',10)
+ xlabel(xlab2); ylabel('log10( norm of mu )');
+ fontsize(10), orient tall, wysiwyg
+
+ % write mu vectors to file
+ if iwrite == 1
+ odir = ['mu_dirs/mu_' stkers{iker} '_' stdcov{idatacov} ];
+ for ip = 1:length(tpinds);
+ pmax = tpinds(ip);
+ filename = sprintf([odir '_p%3.3i'],pmax);
+ fid = fopen([odir '/' filename],'w');
+ for ii = 1:nsrc
+ fprintf(fid,'%24.12e\n',mu_all(ii,ip));
+ end
+ fclose(fid);
+ end
+ end
+ end
+
+
+end
+
+%end
+
+if iwrite == 1
+ % write data-norm vector to file
+ fid = fopen('dmisfit','w');
+ fprintf(fid,'%24.12e\n',dmisfit);
+ fclose(fid);
+
+ % write data-norm vector to file
+ fid = fopen('nwin_tot','w');
+ fprintf(fid,'%i\n',N);
+ fclose(fid);
+
+ % write data-norm vector to file
+ fid = fopen(['data_norm_' stdcov{idatacov}],'w');
+ for isrc = 1:nsrc
+ fprintf(fid,'%24.12e\n',dnorm(isrc));
+ end
+ fclose(fid);
+
+ % write the factors for the data covariance (nevent x 1)
+ fid = fopen(['dcov_fac_' stdcov{idatacov}],'w');
+ for isrc = 1:nsrc
+ fprintf(fid,'%12i\n',dcov_fac_e(isrc));
+ end
+ fclose(fid);
+end
+
+%======================================================================
Modified: seismo/3D/mt_measure_adj/scripts_tomo/run_mt_cc_plot.pl
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/run_mt_cc_plot.pl 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/scripts_tomo/run_mt_cc_plot.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -26,6 +26,13 @@
if($Tmin==2) {$Lend = 120;}
$lcut = "$Lstart/$Lend";
+if ($itest == 0) {
+ $dir_syn = "SYN"; $dir_data = "DATA"; $dir_meas = "MEASURE"; $dir_recon = "RECON";
+} else {
+ $dir_syn = "SYN_TEST"; $dir_data = "DATA_TEST"; $dir_meas = "MEASURE_TEST"; $dir_recon = "RECON_TEST";
+}
+$plot_recon_dir = "PLOTS/${dir_recon}";
+
# delete all figures in PLOTS
`rm PLOTS/*ps PLOTS/*jpg PLOTS/*pdf`;
@@ -70,11 +77,8 @@
print CSH "\\rm -rf ${plot_adj_dir}\n mkdir ${plot_adj_dir}\n";
print CSH "cp OUTPUT_FILES/*adj ${plot_adj_dir}\n";
-if ($itest == 0) {
- $dir_syn = "SYN"; $dir_data = "DATA"; $dir_meas = "MEASURE";
-} else {
- $dir_syn = "SYN_TEST"; $dir_data = "DATA_TEST"; $dir_meas = "MEASURE_TEST";
-}
+# copy reconstructed records into PLOTS
+ print CSH "cp OUTPUT_FILES/*recon.cc ${plot_recon_dir}\n";
#////////////////////////////////////////////////////////
@@ -104,7 +108,7 @@
print CSH "\\mv STATIONS_ADJOINT PLOTS\n";
print CSH "cd PLOTS\n";
-print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b 1 -k 1 -a $stafile2 -d $dir_data -s $dir_syn -w MEASUREMENT.WINDOWS\n";
+print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b 1 -k 1 -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS\n";
print CSH "\\rm test*.pdf\n";
print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n";
print CSH "cd -\n";
Modified: seismo/3D/mt_measure_adj/scripts_tomo/run_mt_measure_adj.pl
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/run_mt_measure_adj.pl 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/scripts_tomo/run_mt_measure_adj.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -57,6 +57,13 @@
print "\n-- $cmtfile --\n";
+if ($itest == 0) {
+ $dir_syn = "SYN"; $dir_data = "DATA"; $dir_meas = "MEASURE"; $dir_recon = "RECON";
+} else {
+ $dir_syn = "SYN_TEST"; $dir_data = "DATA_TEST"; $dir_meas = "MEASURE_TEST"; $dir_recon = "RECON_TEST";
+}
+$plot_recon_dir = "PLOTS/${dir_recon}";
+
#=============================================
# write the C-shell script to file
$cshfile = "run_mt_measure_adj.csh";
@@ -82,8 +89,11 @@
if($iboth==0) {print CSH "\\rm -rf ${plot_adj_dir}\n mkdir ${plot_adj_dir}\n";}
print CSH "cp OUTPUT_FILES/*adj ${plot_adj_dir}\n";
+ # copy reconstructed records into PLOTS
+ print CSH "cp OUTPUT_FILES/*recon.cc ${plot_recon_dir}\n";
+
# copy chi values into PLOTS
- print CSH "cp window_chi PLOTS\n";
+ print CSH "cp window_chi PLOTS\n";
# convert to SAC files (creates *.sac files)
#print CSH "ascii2sac.csh ${plot_adj_dir}*.adj\n";
@@ -97,12 +107,6 @@
print CSH "cp STATIONS_ADJOINT ${adj_dir}\n";
print CSH "\\mv STATIONS_ADJOINT PLOTS\n";
-if ($itest == 0) {
- $dir_syn = "SYN"; $dir_data = "DATA"; $dir_meas = "MEASURE";
-} else {
- $dir_syn = "SYN_TEST"; $dir_data = "DATA_TEST"; $dir_meas = "MEASURE_TEST";
-}
-
# make plots of (filtered) data, synthetics, windows, and adjoint sources
if ($iplot == 1) {
# remove any plot files
@@ -115,7 +119,7 @@
$lcut = "$Lstart/$Lend"; # time interval to cut for plotting
print CSH "cd PLOTS\n";
- print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b $iboth -k $iker -a $stafile2 -d $dir_data -s $dir_syn -w MEASUREMENT.WINDOWS\n";
+ print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b $iboth -k $iker -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS\n";
#print CSH "\\rm mt_cc_all*.pdf\n";
print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n";
print CSH "cd -\n";
Modified: seismo/3D/mt_measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/run_tomo.pl 2008-04-22 18:23:03 UTC (rev 11844)
+++ seismo/3D/mt_measure_adj/scripts_tomo/run_tomo.pl 2008-04-22 18:56:23 UTC (rev 11845)
@@ -102,18 +102,23 @@
print CSH "run_mt_cc_plot.pl $Ts $tvec 0 $iparbools $par1 $par2 $par3\n";
}
- $infile1 = "window_chi"; $outfile1 = "$dir_meas/${ftag}_${infile1}";
- $infile2 = "run_file"; $outfile2 = "$dir_meas/${ftag}_${infile2}";
- $infile3 = "mt_cc_all.pdf"; $outfile3 = "$dir_meas/${ftag}_${infile3}";
- #$infile4 = $outfile1; $outfile4 = "$dir_meas/${ftag}_time_shifts";
+ $infile1 = "window_chi"; $outfile1 = "$dir_meas/${ftag}_${infile1}";
+ $infile2 = "run_file"; $outfile2 = "$dir_meas/${ftag}_${infile2}";
+ $infile3 = "mt_cc_all.pdf"; $outfile3 = "$dir_meas/${ftag}_${infile3}";
+ $infile4 = "MEASUREMENT.PAR"; $outfile4 = "$dir_meas/${ftag}_${infile4}";
+ $infile5 = "MEASUREMENT.WINDOWS"; $outfile5 = "$dir_meas/${ftag}_${infile5}";
+ #$infile6 = $outfile1; $outfile6 = "$dir_meas/${ftag}_time_shifts";
print CSH "rm -rf $dir_meas; mkdir $dir_meas\n";
print CSH "cp $infile1 $outfile1\n";
print CSH "cp $infile2 $outfile2\n";
print CSH "cp PLOTS/$infile3 $outfile3\n";
+ print CSH "cp $infile4 $outfile4\n";
+ print CSH "cp $infile5 $outfile5\n";
+
print CSH "rm -rf $dir_adj\n";
print CSH "cp -r ADJOINT_SOURCES $dir_adj\n";
- #print CSH "awk '{print \$15}' $infile4 > $outfile4\n";
+ #print CSH "awk '{print \$15}' $infile6 > $outfile6\n";
print CSH "echo done with Event $eid -- $ievent out of $imax/$nevent\n";
# copy another copy to a single directory
More information about the cig-commits
mailing list