[cig-commits] r11922 - in seismo/3D/mt_measure_adj: . PLOTS scripts_tomo scripts_tomo/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Wed May 7 06:27:12 PDT 2008
Author: carltape
Date: 2008-05-07 06:27:12 -0700 (Wed, 07 May 2008)
New Revision: 11922
Added:
seismo/3D/mt_measure_adj/PLOTS/make_pdf_all.pl
seismo/3D/mt_measure_adj/scripts_tomo/matlab/read_window_chi.m
Modified:
seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl
seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl
seismo/3D/mt_measure_adj/mt_constants.f90
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/prepare_mt_measure_adj.pl
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:
Changed output for window_chi files to include a waveform difference measure for the entire record containing each measurement window. Also modified the output plotting order to allow for sorting due to station name, distance, or azimuth. Also updated Matlab scripts for the subspace method.
Added: seismo/3D/mt_measure_adj/PLOTS/make_pdf_all.pl
===================================================================
--- seismo/3D/mt_measure_adj/PLOTS/make_pdf_all.pl (rev 0)
+++ seismo/3D/mt_measure_adj/PLOTS/make_pdf_all.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -0,0 +1,51 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# make_pdf_all.pl
+# Carl Tape
+# 06-May-2008
+#
+# This sorts the output PDF files in PLOTS into order by
+# 1 station
+# 2 distance
+# 3 azimuth
+#
+# EXAMPLE (from PLOTS): make_pdf_all.pl 3 mt_cc_all
+#
+#==========================================================
+
+if (@ARGV < 2) {die("Usage: make_pdf_all.pl xxx\n")}
+($isort,$otag) = @ARGV;
+
+ at tags = ("sta","dist","az");
+$tag = $tags[$isort-1];
+
+$ofile = "${otag}_${tag}.pdf";
+
+$station_sort = "STATIONS_${tag}";
+if (not -f "${station_sort}") {die("check if station_sort ${station_sort} exist or not\n");}
+ at pdffiles = glob("*pdf");
+$npdf = @pdffiles;
+
+print "\n $npdf files\n";
+
+for ($ik = 1; $ik <= $npdf; $ik = $ik+1) {
+
+ $pdffile0 = $pdffiles[$ik-1];
+ ($pdffile) = split(" ",`basename $pdffile0`);
+ ($sta,$net,undef,undef,undef) = split("_",$pdffile);
+ $stanet = "$sta.$net";
+ $iline=`grep -n \"$stanet\" ${station_sort} | awk -F: '{print \$1}'`; chomp($iline);
+ print "$ik $stanet $iline\n";
+
+ $sti = sprintf("%3.3i",$iline);
+ $pdffile_out = "${sti}_${pdffile}";
+ `\\mv $pdffile0 ${pdffile_out}`;
+ #print "move $iline $sti $pdffile0 ${pdffile_out}\n";
+}
+
+# concatenate into one file
+`/home/carltape/bin/pdcat -r *.pdf $ofile`;
+
+#=================================================================
Property changes on: seismo/3D/mt_measure_adj/PLOTS/make_pdf_all.pl
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/PLOTS/plot_win_adj.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -271,7 +271,7 @@
push @Twinb, "$winb";
push @Twine, "$wine";
- (undef,undef,undef,undef,undef,$kplotT,$chi1T,$chi2T,$chi3T,$chi4T,$chi5T,$meas1T,$meas2T,$meas3T,$meas4T,$meas5T,$sigma1T,$sigma2T,$sigma3T,$sigma4T,$sigma5T) = split(" ",$mlines[$i]);
+ (undef,undef,undef,undef,undef,$kplotT,undef,undef,$chi2T,$chi3T,$chi4T,$chi5T,$meas2T,$meas3T,$meas4T,$meas5T,$sigma2T,$sigma3T,$sigma4T,$sigma5T) = split(" ",$mlines[$i]);
$stlabT_dT[$i] = sprintf("dT = %.2f \\261 %.2f",$meas4T,$sigma4T);
$stlabT_dA[$i] = sprintf("dA = %.2f \\261 %.2f",$meas5T,$sigma5T);
}
@@ -297,7 +297,7 @@
push @Rwinb, "$winb";
push @Rwine, "$wine";
- (undef,undef,undef,undef,undef,$kplotR,$chi1R,$chi2R,$chi3R,$chi4R,$chi5R,$meas1R,$meas2R,$meas3R,$meas4R,$meas5R,$sigma1R,$sigma2R,$sigma3R,$sigma4R,$sigma5R) = split(" ",$mlines[$i]);
+ (undef,undef,undef,undef,undef,$kplotR,undef,undef,$chi2R,$chi3R,$chi4R,$chi5R,$meas2R,$meas3R,$meas4R,$meas5R,$sigma2R,$sigma3R,$sigma4R,$sigma5R) = split(" ",$mlines[$i]);
$stlabR_dT[$i] = sprintf("dT = %.2f \\261 %.2f",$meas4R,$sigma4R);
$stlabR_dA[$i] = sprintf("dA = %.2f \\261 %.2f",$meas5R,$sigma5R);
}
@@ -323,7 +323,7 @@
push @Zwinb, "$winb";
push @Zwine, "$wine";
- (undef,undef,undef,undef,undef,$kplotZ,$chi1Z,$chi2Z,$chi3Z,$chi4Z,$chi5Z,$meas1Z,$meas2Z,$meas3Z,$meas4Z,$meas5Z,$sigma1Z,$sigma2Z,$sigma3Z,$sigma4Z,$sigma5Z) = split(" ",$mlines[$i]);
+ (undef,undef,undef,undef,undef,$kplotZ,undef,undef,$chi2Z,$chi3Z,$chi4Z,$chi5Z,$meas2Z,$meas3Z,$meas4Z,$meas5Z,$sigma2Z,$sigma3Z,$sigma4Z,$sigma5Z) = split(" ",$mlines[$i]);
$stlabZ_dT[$i] = sprintf("dT = %.2f \\261 %.2f",$meas4Z,$sigma4Z);
$stlabZ_dA[$i] = sprintf("dA = %.2f \\261 %.2f",$meas5Z,$sigma5Z);
}
Modified: seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl
===================================================================
--- seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/PLOTS/plot_win_stats_all.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -33,7 +33,7 @@
# directories
$odir = "/net/sierra/raid1/carltape/results/MEASUREMENTS";
-$dir_all = "$odir/ALL_${sTmin}";
+$dir_all = "$odir/$smodel/ALL_${sTmin}";
$dir_run = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
if (not -e $dir_all) {die("check if $dir_all exist or not\n")}
if (not -e $dir_run) {die("check if $dir_run exist or not\n")}
Modified: seismo/3D/mt_measure_adj/mt_constants.f90
===================================================================
--- seismo/3D/mt_measure_adj/mt_constants.f90 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/mt_constants.f90 2008-05-07 13:27:12 UTC (rev 11922)
@@ -38,6 +38,9 @@
! total number of possible measuement types (waveform, MT-tt, MT-amp, XC-tt, XC-amp)
integer, parameter :: N_MEASUREMENT = 5
+! number of elements in window_chi
+ integer, parameter :: NCHI = 3*(N_MEASUREMENT-1) + 8
+
integer, parameter :: IKER_WV = 0
integer, parameter :: IKER_MT = 1
integer, parameter :: IKER_CC = 2
Modified: seismo/3D/mt_measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/mt_measure_adj/mt_measure_adj.f90 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/mt_measure_adj.f90 2008-05-07 13:27:12 UTC (rev 11922)
@@ -15,7 +15,7 @@
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, 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, dimension(NCHI) :: window_chi
double precision :: fend0, fstart0, fend, fstart
!double precision :: TSHORT, TLONG ! mtm_variables.f90
@@ -220,6 +220,18 @@
print *
print *, ' Measurements No.', j, ' ... '
+ ! initialize the measurements
+ window_chi(:) = 0.
+
+ ! compute integrated waveform difference, normalized by duration of the record
+ ! NOTE: (1) for comparison with waveform_chi, we include the 0.5 factor
+ ! (2) we might want to include dt as an integration factor (also for waveform_chi),
+ ! but the ratio (d-s)^2 / d^2 avoids the need for dt, nstep, or length of record
+ window_chi(17) = 0.5 * sum( data**2 )
+ window_chi(18) = 0.5 * sum( syn**2 )
+ window_chi(19) = 0.5 * sum( (data-syn)**2 )
+ window_chi(20) = npts*dt
+
! get the starting frequency to avoid measuring long periods not present
!fstart = fstart0 ; fend = fend0
!call mt_measure_select_0(nlen,dt,i_left,i_right,fstart,fend,use_trace)
@@ -289,10 +301,13 @@
endif
! KEY: write misfit function values to file (two for each window)
- ! 1: waveform chi, 2: MT-TT chi, 3: MT-dlnA chi, 4: XC-TT chi, 5: XC-dlnA chi
- ! 6: waveform chi, 7: MT-TT meas, 8: MT-dlnA meas, 9: XC-TT meas, 10: XC-dlnA meas
- ! 11: waveform chi, 12: MT-TT error, 13: MT-dlnA error, 14: XC-TT error, 15: XC-dlnA error
- write(13,'(a14,a8,a3,a5,i4,i4,15e14.6,2e14.6)') file_prefix0,sta,net,chan,j,iker,window_chi(:),tr_chi,am_chi
+ ! 1: MT-TT chi, 2: MT-dlnA chi, 3: XC-TT chi, 4: XC-dlnA chi
+ ! 5: MT-TT meas, 6: MT-dlnA meas, 7: XC-TT meas, 8: XC-dlnA meas
+ ! 9: MT-TT error, 10: MT-dlnA error, 11: XC-TT error, 12: XC-dlnA error
+ ! WINDOW : 13: data power, 14: syn power, 15: (data-syn) power, 16: window duration
+ ! FULL RECORD: 17: data power, 18: syn power, 19: (data-syn) power, 20: record duration
+ write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6)') file_prefix0,sta,net,chan,j,iker,&
+ tstart,tend,window_chi(:),tr_chi,am_chi
print *, ' tr_chi = ', sngl(tr_chi), ' am_chi = ', sngl(am_chi)
! combine adjoint sources from different measurement windows
Modified: seismo/3D/mt_measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/mt_measure_adj/mt_sub.f90 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/mt_sub.f90 2008-05-07 13:27:12 UTC (rev 11922)
@@ -74,17 +74,19 @@
filename = trim(file_prefix)
- !--------------------------------------------------------------------------
- ! window and interpolate data and synthetics
- !--------------------------------------------------------------------------
-
if (DISPLAY_DETAILS) then
call dwrite_sacfile_f(datafile,trim(file_prefix)//'_data.sac',t0,npts,data)
call dwrite_sacfile_f(datafile,trim(file_prefix)//'_syn.sac',t0,npts,syn)
+ !call dwrite_sacfile_f(datafile,trim(file_prefix)//'_diff_dms.sac',t0,npts,data-syn)
!call dwrite_ascfile_f(trim(file_prefix)//'_data.txt',t0,dt,npts,data)
!call dwrite_ascfile_f(trim(file_prefix)//'_syn.txt',t0,dt,npts,syn)
endif
+
+ !--------------------------------------------------------------------------
+ ! window and interpolate data and synthetics
+ !--------------------------------------------------------------------------
+
! interpolate data and synthetics, and also extract time-windowed records
call interpolate_data_and_syn(data, syn, tstart, tend, t0, dt, NPT, dzr2, dzr, nlen, istart)
@@ -135,7 +137,7 @@
! call dwrite_sacfile_f(datafile,'windowed_shifted_data.sac',tstart,nlen,dzr20)
!endif
- ! amplitude MEASUREMENT (definition of Dahlen and Baig (2002), Eq. 3,17,18 : dlnA = Aobs/Asyn - 1)
+ ! compute velocity of synthetics and data
syn_displ(1:nlen) = dzr0(1:nlen)
dat_displ(1:nlen) = dzr20(1:nlen)
@@ -154,6 +156,7 @@
syn_accel(1) = (syn_veloc(2) - syn_veloc(1)) / dt
syn_accel(nlen) = (syn_veloc(nlen) - syn_veloc(nlen-1)) / dt
+ ! amplitude MEASUREMENT (definition of Dahlen and Baig (2002), Eq. 3,17,18 : dlnA = Aobs/Asyn - 1)
! note that the records are UN-SHIFTED
dlnA = sqrt( (dt * sum( dat_displ(1:nlen) * dat_displ(1:nlen) )) &
/ (dt * sum( syn_displ(1:nlen) * syn_displ(1:nlen) )) ) - 1.
@@ -565,7 +568,7 @@
double precision, dimension(:), intent(out) :: tr_adj_src
double precision, intent(out) :: tr_chi
- double precision, dimension(3*N_MEASUREMENT), intent(out) :: window_chi
+ double precision, dimension(NCHI), intent(inout) :: window_chi
double precision, dimension(:), intent(out), optional :: am_adj_src
double precision, intent(out), optional :: am_chi
@@ -579,7 +582,7 @@
complex*16, dimension(NPT) :: pwc_adj,qwc_adj
double precision, dimension(NPT) :: dtau_pj_t, dlnA_qj_t
double precision :: dtau_wtr, dlnA_wtr, err_t, err_A
- double precision :: waveform_chi, waveform_temp
+ double precision :: waveform_chi, waveform_d2, waveform_s2, waveform_temp1, waveform_temp2, waveform_temp3
! waveform adjoint source is passed by tr_adj_src and tr_chi
if (iker == 0 .and. (present(am_adj_src) .or. present(am_chi))) stop 'am_adj_src and am_chi are not needed for iker = 0 (waveform adjoint source case)'
@@ -762,9 +765,10 @@
ishift = int(tshift_xc/dt)
- ! CHT: compute the integrated waveform difference squared,
- ! no matter what measurement you make
- waveform_temp = 0.
+ ! CHT: Compute the integrated waveform difference squared,
+ ! no matter what measurement you make.
+ ! To do this, we must UNSHIFT the data.
+ waveform_temp1 = 0. ; waveform_temp2 = 0. ; waveform_temp3 = 0.
do i = 1,nlen
i1 = istart + i
if (i-ishift < 1 .or. i-ishift > nlen) then
@@ -772,9 +776,14 @@
else
i2 = i-ishift
endif
- waveform_temp = waveform_temp + (( dzr0(i) - dzr20(i2) ) * time_window(i))**2
+ waveform_temp1 = waveform_temp1 + ( dzr20(i2) * time_window(i) )**2
+ waveform_temp2 = waveform_temp2 + ( dzr0(i) * time_window(i) )**2
+ waveform_temp3 = waveform_temp3 + (( dzr20(i2) - dzr0(i) ) * time_window(i) )**2
enddo
- waveform_chi = 0.5 * dt * waveform_temp
+ ! NOTE: does not include DT factor or normalization by duration of window
+ waveform_d2 = 0.5 * waveform_temp1
+ waveform_s2 = 0.5 * waveform_temp2
+ waveform_chi = 0.5 * waveform_temp3
! compute traveltime and amplitude adjoint sources
do i = 1,nlen
@@ -809,34 +818,36 @@
! 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
- ! 4: cross-correlation, TT
- ! 5: cross-correlation, dlnA
- window_chi(:) = 0.
+ ! 1: multitaper, TT
+ ! 2: multitaper, dlnA
+ ! 3: cross-correlation, TT
+ ! 4: cross-correlation, dlnA
+ !window_chi(:) = 0.
! misfit function value
- window_chi(1) = waveform_chi
- if(iker==1) window_chi(2) = 0.5 * 2.0 * df * sum( (dtau_w(1:i_right))**2 * wp_taper(1:i_right) )
- if(iker==1) window_chi(3) = 0.5 * 2.0 * df * sum( (dlnA_w(1:i_right))**2 * wq_taper(1:i_right) )
- window_chi(4) = 0.5 * (tshift_xc/sigma_dt_cc)**2
- window_chi(5) = 0.5 * (dlnA/sigma_dlnA_cc)**2
+ if(iker==1) window_chi(1) = 0.5 * 2.0 * df * sum( (dtau_w(1:i_right))**2 * wp_taper(1:i_right) )
+ if(iker==1) window_chi(2) = 0.5 * 2.0 * df * sum( (dlnA_w(1:i_right))**2 * wq_taper(1:i_right) )
+ window_chi(3) = 0.5 * (tshift_xc/sigma_dt_cc)**2
+ window_chi(4) = 0.5 * (dlnA/sigma_dlnA_cc)**2
! measurement (no uncertainty estimates)
- window_chi(6) = window_chi(1)
- if(iker==1) window_chi(7) = sum( dtau_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
- if(iker==1) window_chi(8) = sum( dlnA_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
- window_chi(9) = tshift_xc
- window_chi(10) = dlnA
+ if(iker==1) window_chi(5) = sum( dtau_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
+ if(iker==1) window_chi(6) = sum( dlnA_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
+ window_chi(7) = tshift_xc
+ window_chi(8) = dlnA
! estimated mesurement uncertainties
- window_chi(11) = window_chi(1)
- if(iker==1) window_chi(12) = sigma_dt
- if(iker==1) window_chi(13) = sigma_dlnA
- window_chi(14) = sigma_dt_cc
- window_chi(15) = sigma_dlnA_cc
+ if(iker==1) window_chi(9) = sigma_dt
+ if(iker==1) window_chi(10) = sigma_dlnA
+ window_chi(11) = sigma_dt_cc
+ window_chi(12) = sigma_dlnA_cc
+ ! for normalization, divide by duration of window
+ window_chi(13) = waveform_d2
+ window_chi(14) = waveform_s2
+ window_chi(15) = waveform_chi
+ window_chi(16) = nlen*dt
+
!!$ open(88,file='testing.dat')
!!$ do i = 1,i_right
!!$ write(88,'(5e18.6)') df*i, dtau_w(i), dlnA_w(i), wp_taper(i), wq_taper(i)
@@ -844,15 +855,15 @@
!!$ close(88)
if(iker==0) then ! waveform
- tr_chi = window_chi(1)
+ tr_chi = waveform_chi
elseif(iker==1) then ! multitaper
- tr_chi = window_chi(2)
- am_chi = window_chi(3)
+ tr_chi = window_chi(1)
+ am_chi = window_chi(2)
elseif(iker==2) then ! cross_correlation
- tr_chi = window_chi(4)
- am_chi = window_chi(5)
+ tr_chi = window_chi(3)
+ am_chi = window_chi(4)
endif
@@ -1321,7 +1332,7 @@
if (ishift > 0) tseis_recon_cc_dt(1:ishift+1) = tseis_recon_cc_dt(ishift+2)
if (ishift < 0) tseis_recon_cc_dt(nlen+ishift:nlen) = tseis_recon_cc_dt(nlen+ishift-1)
- ! apply transfer function to synthetics (phase and amplitude)
+ ! apply cross-correlation amplitude measurement
tseis_recon_cc(:) = 0.
tseis_recon_cc(:) = tseis_recon_cc_dt * (1. + dlnA)
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m 2008-05-07 13:27:12 UTC (rev 11922)
@@ -14,24 +14,30 @@
%-------------------------
-clear
+%clear
close all
format compact
dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
ftag = '_window_chi';
-stmod = 'm0';
- % iteration label (m0, m1, m2, ...)
Tvec = [6 2];
% LOWEST period used for each measurement band-pass (upper is 40 s)
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);
+%DT = 0.011; % really this is only needed because it was left out from
+ % integrating hte waveforms in mt-Measure_adj.f90
+
+%-------------------------
+% user input
+
+imod = input(' Enter model number (0, 1, ...): ');
+stmod = sprintf('m%i',imod);
+
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
@@ -48,7 +54,7 @@
% load the event IDs corresponding to the kernels
% load these as strings, since event IDs could have letters
-stker0 = '/net/sierra/raid1/carltape/results/KERNELS/';
+stker0 = '/net/sierra/raid1/carltape/results/KERNELS/kernel_m0/';
%eids = textread([stker0 'kernels_test'],'%s');
eids = textread([stker0 'kernels_done_mod'],'%s');
nevent = length(eids);
@@ -100,22 +106,20 @@
%-----------------------
-%windows_array = zeros(nevent,nrec,ncomp,nper); % number of windows
-
iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
if iread == 1
k2 = 0;
meas_array = [];
- %for tt = 1:1
- for tt = 1:nper
+ 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
+ for ii = 3:3
+ %for ii = 1:nevent
disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
@@ -125,10 +129,17 @@
if ~exist(filename,'file')
disp('--> file does not exist (or nwin = 0)');
else
- [stnet0,strec0,comp,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);
+ [stnet0,strec0,comp,iwin,iker,t1,t2,...
+ chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
+ measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
+ sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,...
+ wind2,wins2,win_chi0,windur, recd2,recs2,rec_chi0,recdur,...
+ tr_chi,am_chi] = read_window_chi(filename);
+
+ % waveform differences normalized by data-squared
+ win_chi = win_chi0 ./ wind2;
+ rec_chi = rec_chi0 ./ recd2;
+
nwin = length(strec0);
disp(['--> nwin = ' num2str(nwin)]);
@@ -159,7 +170,7 @@
kinds = [k1:k2]';
meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin ...
- iker chi1 measCC_dT sigmaCC_dT measMT_dT sigmaMT_dT tr_chi];
+ iker measCC_dT sigmaCC_dT measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
% 1 kinds
% 2 index_T
% 3 index_event
@@ -168,12 +179,13 @@
% 6 index_comp
% 7 iwin
% 8 iker
- % 9 chi1
- % 10 measCC_dT
- % 11 sigmaCC_dT
- % 12 measMT_dT
- % 13 sigmaMT_dT
- % 14 tr_chi
+ % 9 measCC_dT
+ % 10 sigmaCC_dT
+ % 11 measMT_dT
+ % 12 sigmaMT_dT
+ % 13 win_chi
+ % 14 rec_chi
+ % 15 tr_chi
end
end
@@ -190,6 +202,80 @@
% total number of windows
N = length(meas_array);
+%----------------------------------------------
+
+% for each record for each model (m0, m1) that has at least one window picked,
+% we save the integrated waveform difference as the purest measure of misfit
+% between the data and synthetics
+nmodel = 2;
+%wdiff_array = zeros(nevent,nrec,ncomp,nper,nmodel);
+mm = imod+1;
+
+for tt = 1:nper
+ for ii = 1:nevent
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ imatch = find( and( and( meas_array(:,2)==tt, meas_array(:,3)==ii), ...
+ and( meas_array(:,5)==jj, meas_array(:,6)==kk) ));
+ if ~isempty(imatch)
+ % take the first only
+ wdiff_array(ii,jj,kk,tt,mm) = meas_array(imatch(1),14);
+ end
+ end
+ end
+ end
+end
+
+break
+
+ratio_array = []; ll = 0;
+for tt = 1:nper
+ for ii = 1:nevent
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ if prod(wdiff_array(ii,jj,kk,tt,:)) > 0
+ ll = ll+1;
+ rat = wdiff_array(ii,jj,kk,tt,2) / wdiff_array(ii,jj,kk,tt,1);
+ ratio_array(ll,:) = [tt ii jj kk rat];
+ end
+ end
+ end
+ end
+end
+
+ratio_sort = sortrows(ratio_array,5);
+
+tarray = []; ll = 0;
+for jj = 1:nrec
+ imatch = find( ratio_array(:,3)==jj);
+ if length(imatch)==3
+ ll = ll+1;
+ tarray(ll,:) = [jj prod(ratio_array(imatch,5)) ];
+ end
+end
+tarray_sort = sortrows(tarray,2);
+for hh = 1:20
+ jj = tarray_sort(hh,1);
+ disp(sprintf('%7s %6.3f',char(strec{jj}),tarray_sort(hh,2) ));
+end
+
+for hh = 1:20
+ tt = ratio_sort(hh,1);
+ ii = ratio_sort(hh,2);
+ jj = ratio_sort(hh,3);
+ kk = ratio_sort(hh,4);
+ disp(sprintf('%10s --> %6.1f %7s %5s %6.3f',char(eids{ii}),...
+ Tvec(tt),char(strec{jj}),char(comps{kk}),ratio_sort(hh,5) ));
+end
+
+figure; N = plot_histo(ratio_array(:,end),[0:0.2:5]);
+figure; N = plot_histo(log(ratio_array(:,end)),[-4:0.2:4]);
+
+break
+
+%----------------------------------------------
% check some of the output window rows
if 0==1
itestvec = [1 2];
@@ -198,9 +284,9 @@
itest = itestvec(ii);
meas_array(itest,:)
display_meas(meas_array(itest,:),Tvec,eids,strec,comps);
- meas_array(itest,13)
- 0.5 * meas_array(itest,10)^2 / meas_array(itest,11)^2
- 0.5 * meas_array(itest,12)^2 / meas_array(itest,11)^2
+ meas_array(itest,12)
+ 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
@@ -255,12 +341,12 @@
% write out the file
fid = fopen(['window_picks_by_' stks{kk}],'w');
- fprintf(fid,'%12s%10s%10s%10s\n',stks{kk},...
+ fprintf(fid,'%12s%5s%10s%10s%10s\n',stks{kk},' ',...
['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
- fprintf(fid,'%12s%10i%10i%10i\n','TOTAL',ntot_1,sum(ntot_1));
+ fprintf(fid,'%12s%5s%10i%10i%10i\n','TOTAL -->',' ',ntot_1,sum(ntot_1));
for ii = 1:nw
jj = isort(ii);
- fprintf(fid,'%12s%10i%10i%10i\n',labs{jj},nwin_out(jj,:),ntot_2(jj));
+ fprintf(fid,'%12s%5i%10i%10i%10i\n',labs{jj},jj,nwin_out(jj,:),ntot_2(jj));
end
fclose(fid);
end
@@ -294,7 +380,7 @@
dnorm_sq = zeros(nevent,1);
for ii = 1:nevent
imatch = find( meas_array(:,3)==ii );
- dnorm_sq(ii) = sum(2 * meas_array(imatch,14) ./ dcov_fac(imatch) );
+ dnorm_sq(ii) = sum(2 * meas_array(imatch,15) ./ dcov_fac(imatch) );
end
% weight vector for the subspace method
@@ -332,7 +418,7 @@
% % find records that meet particular criteria
% DT_MIN = 5;
% DT_SIGMA_MIN = 1;
-% icheck = find( and( meas_array(:,10) >= DT_MIN, meas_array(:,11) <= DT_SIGMA_MIN) )
+% 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);
@@ -340,7 +426,7 @@
% find records that meet particular criteria
CHI_MIN = 100;
DT_SIGMA_MIN = 0.19;
-icheck = find( and( meas_array(:,14) >= CHI_MIN, meas_array(:,11) < DT_SIGMA_MIN) );
+icheck = find( and( meas_array(:,15) >= 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);
Added: seismo/3D/mt_measure_adj/scripts_tomo/matlab/read_window_chi.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/read_window_chi.m (rev 0)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/read_window_chi.m 2008-05-07 13:27:12 UTC (rev 11922)
@@ -0,0 +1,32 @@
+%
+% read_window_chi.m
+% CARL TAPE, 14-March-2008
+% printed xxx
+%
+% This file reads in the output file from mt_measure_adj.f90.
+%
+% calls xxx
+% called by compute_misfit.m
+%
+
+function [netwk,strec,cmp,iwin,iker,t1,t2,...
+ chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
+ measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
+ sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,...
+ wind2,wins2,winchi,windur,...
+ recd2,recs2,recchi,recdur,...
+ tr_chi,am_chi] = read_window_chi(filename)
+
+[flab,sta,netwk,cmp,iwin,iker,t1,t2,...
+ chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
+ measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
+ sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,...
+ wind2,wins2,winchi,windur,...
+ recd2,recs2,recchi,recdur,...
+ tr_chi,am_chi] = ...
+ textread(filename,'%s%s%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f');
+
+for kk = 1:length(sta), strec{kk} = [sta{kk} '.' netwk{kk}]; end
+strec = strec(:);
+
+%======================================================
Modified: seismo/3D/mt_measure_adj/scripts_tomo/prepare_mt_measure_adj.pl
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2008-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -108,23 +108,33 @@
# NOTE: NO EXTENSION IS PUT ON THE BAND-PASSED VERSIONS
if ($iplot == 1) {
# synthetics
- $syn_dir_plot = "PLOTS/${dir_syn}";
- print CSH "\\rm -rf $syn_dir_plot\n";
- print CSH "mkdir $syn_dir_plot\n";
+ ${syn_dir_plot} = "PLOTS/${dir_syn}";
+ print CSH "\\rm -rf ${syn_dir_plot}\n";
+ print CSH "mkdir ${syn_dir_plot}\n";
print CSH "cd ${dir_meas_syn}\n";
print CSH "process_trinet_syn_new.pl -S -t $Ts -d ../${syn_dir_plot} *sac.d\n";
print CSH "cd -\n";
# data
$data_dir_plot = "PLOTS/${dir_data}";
- print CSH "\\rm -rf $data_dir_plot\n";
- print CSH "mkdir $data_dir_plot\n";
+ print CSH "\\rm -rf ${data_dir_plot}\n";
+ print CSH "mkdir ${data_dir_plot}\n";
print CSH "cd ${dir_meas_data}\n";
print CSH "process_cal_data.pl -t $Ts -d ../${data_dir_plot} *sac.d\n";
print CSH "cd -\n";
- # remove any plots in PLOTS
+ # remove any old plots in PLOTS directory
print CSH "\\rm PLOTS/*pdf PLOTS/*jpg PLOTS/*ps\n";
+
+ # make a set of SYNTHETIC stations lists sorted by: station, distance, azimuth
+ @tags = ("sta","dist","az");
+ $fall = "${syn_dir_plot}/*sac.d";
+ for ($ij = 1; $ij <= 3; $ij = $ij+1) {
+ $tag = $tags[$ij-1];
+ $ofile = "PLOTS/STATIONS_${tag}";
+ print CSH "\\rm -rf $ofile\n";
+ print CSH "saclst kstnm knetwk dist az f $fall | awk '{print \$2\".\"\$3,\$4,\$5}' | uniq | sort -g -k $ij > $ofile\n";
+ }
}
## this version allows for flexibility for the exact name of the window file
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-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/scripts_tomo/run_mt_cc_plot.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -18,6 +18,8 @@
if (@ARGV < 7) {die("Usage: run_mt_cc_plot.pl Tmin/Tmax tstart/dt/npt itest iparbools par1 par2 par3\n")}
($Ts,$tvec,$itest,$iparbools,$par1,$par2,$par3) = @ARGV;
+$pwd = $ENV{PWD};
+
($Tmin,$Tmax) = split("/",$Ts);
# time interval to cut for plotting
@@ -110,7 +112,13 @@
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 -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";
+
+# make a single PDF file
+$isort = 3;
+$otag = "mt_cc_all";
+print CSH "make_pdf_all.pl $isort $otag\n";
+
+#print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n"; # replace with make_pdf_all.pl
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-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/scripts_tomo/run_mt_measure_adj.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -120,8 +120,14 @@
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 -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 "\\rm test*.pdf\n";
+
+ # make a single PDF file
+ $isort = 3;
+ $otag = "mt_cc_all";
+ print CSH "make_pdf_all.pl $isort $otag\n";
+
+ #print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n"; # replace with make_pdf_all.pl
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-05-07 00:35:54 UTC (rev 11921)
+++ seismo/3D/mt_measure_adj/scripts_tomo/run_tomo.pl 2008-05-07 13:27:12 UTC (rev 11922)
@@ -22,7 +22,7 @@
$sTmin = sprintf("T%2.2i",$Tmin);
$dir0 = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
-$dir_all = "/net/sierra/raid1/carltape/results/MEASUREMENTS/ALL_${sTmin}";
+$dir_all = "/net/sierra/raid1/carltape/results/MEASUREMENTS/${smodel}/ALL_${sTmin}";
if (not -e $dir0) {die("check if $dir0 exist or not\n")}
if (not -e $dir_all) {die("check if $dir_all exist or not\n")}
@@ -45,8 +45,8 @@
#die("testing");
$imin = 1; $imax = $nevent;
-#$imin = 1; $imax = 172;
-#$imin = 133; $imax = $imin;
+#$imin = 1; $imax = 20;
+#$imin = 43; $imax = $imin;
# EXAMPLE lines from EIDs_simulation_duration
# 5 10222753 300.0 27300 0.26 332.0 4.010
@@ -82,8 +82,10 @@
$window_file = "${dir_win}/MEASUREMENT_WINDOWS_${ftag}";
- if (-e $dir_adj) {
- print "--> DO NOT RUN: adjoint source directory already exists ($dir_adj)\n";
+ #if (-e $dir_adj) {
+ if (-e $dir_meas) {
+ #print "--> DO NOT RUN: adjoint source directory already exists ($dir_adj)\n";
+ print "--> DO NOT RUN: measurement directory already exists ($dir_meas)\n";
} else {
@@ -104,7 +106,9 @@
$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}";
+ #$tfile = `ls PLOTS/*mt_cc*`; $infile3 = `basename $tfile`;
+ $infile3 = "mt_cc_all_az.pdf"; # NEED to make this flexible (make make_pdf_all.pl)
+ $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";
@@ -119,10 +123,11 @@
print CSH "rm -rf $dir_adj\n";
print CSH "cp -r ADJOINT_SOURCES $dir_adj\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
print CSH "cp -r $outfile3 $dir_all\n";
+
+ print CSH "echo done with Event $eid -- $ievent out of $imax/$nevent\n";
}
}
More information about the cig-commits
mailing list