[cig-commits] r12573 - seismo/3D/ADJOINT_TOMO/iterate_adj/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Wed Aug 6 22:26:54 PDT 2008
Author: carltape
Date: 2008-08-06 22:26:54 -0700 (Wed, 06 Aug 2008)
New Revision: 12573
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m
Log:
Added a Matlab script to compare the misfit values for different models. Updated two other Matlab scripts for computing misfit values.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2008-08-07 05:26:54 UTC (rev 12573)
@@ -0,0 +1,302 @@
+%
+% compare_misfit.m
+% CARL TAPE, 23-July-2008
+% printed xxx
+%
+% This script reads in a set of window_chi output measurement files from
+% mt_measure_adj.f90, and it tabulates and plots misfit function values.
+%
+% /cig/seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/
+%
+% Copied from compute_misfit.m on 23-July-2008
+%
+% calls read_window_chi.m, plot_histo.m, readCMT.m, read_station_SPECFEM.m
+% called by xxx
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
+
+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 the waveforms in mt-measure_adj.f90
+
+%-------------------------
+% USER INPUT
+
+iwrite = input(' Enter 1 to write files (0 otherwise) : ');
+
+% files: event IDs, CMTSOLUTIONS, ststions
+dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_8/';
+%eid_file = [dir_source 'EIDs_only_loc'];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_m05';
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v8'];
+stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
+
+%-------------------------
+% 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
+eids = textread(eid_file,'%s');
+eids_num = textread(eid_file,'%f'); % numeric version
+%eids_num = str2num(char(eids));
+nevent = length(eids);
+
+% read in CMT solutions
+[date,tshift,hdur,slat,slon,dep,M,eid,elabel] = readCMT(cmt_file_all, 13, 1);
+
+%-----------------------
+% read in stations file (used for synthetics)
+
+[rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(stations_file);
+nrec = length(stnm);
+for ii = 1:nrec, strec{ii} = [stnm{ii} '.' netwk{ii}]; end
+strec = strec(:);
+
+% make list of networks
+stnet = unique(netwk);
+nnet = length(stnet);
+
+%-----------------------
+
+iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
+
+imodel_min = input(' Enter minimum model index (0-5) : ');
+imodel_max = input(' Enter maximum model index (1-5) : ');
+nmodel = imodel_max - imodel_min;
+
+wdiff_array = zeros(nevent,nrec,ncomp,nper,nmodel);
+
+if iread == 1
+
+ % range of events
+ imin = 1; imax = nevent; % default
+ %imin = 4; imax = 7;
+ %imin = 4; imax = imin;
+
+ for imod = imodel_min:imodel_max
+
+ stmod = sprintf('m%2.2i',imod);
+ disp('====================================');
+ disp(['MODEL ' stmod]);
+
+ k2 = 0;
+ meas_array = [];
+ %for tt = 1:1
+ for tt = 1:nper
+ Tper = Tvec(tt);
+ sTper = sprintf('%3.3i',Tper);
+ if imod <= 2, sTper2 = sprintf('%2.2i',Tper); else sTper2 = sTper; end % CHT
+ disp('-----------------------------------------');
+
+ for ii = imin:imax % loop over events
+
+ disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
+ dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
+ mfile = [eids{ii} '_T' sTper2 '_window_chi'];
+ filename = [dir1 mfile];
+
+ if ~exist(filename,'file')
+ disp('--> file does not exist (or nwin = 0)');
+ else
+ % read the window_chi file
+ [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)]);
+
+ % assign the string network an integer index
+ index_net = NaN * ones(nwin,1);
+ for inet = 1:nnet
+ ind_net = find( strcmp(stnet(inet), stnet0) == 1 );
+ index_net(ind_net) = inet;
+ end
+
+ % 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_comp = find( strcmp(comps(icomp), comp) == 1 );
+ index_comp(ind_comp) = 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_net index_rec index_comp iwin ...
+ iker measCC_dT sigmaCC_dT measCC_dA sigmaCC_dA measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
+ % 1 kinds
+ % 2 index_T
+ % 3 index_event
+ % 4 index_network
+ % 5 index_rec
+ % 6 index_comp
+ % 7 iwin
+ % 8 iker
+ % 9 measCC_dT
+ % 10 sigmaCC_dT
+ % 11 measCC_dA
+ % 12 sigmaCC_dA
+ % 13 measMT_dT
+ % 14 sigmaMT_dT
+ % 15 win_chi
+ % 16 rec_chi
+ % 17 tr_chi
+ end
+
+ end
+ end
+
+ for tt = 1:nper
+ for ii = imin:imax
+ 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,imod+1) = meas_array(imatch(1),16);
+ end
+ end
+ end
+ end
+ end
+
+ end
+
+ save('wdiff_array.mat','wdiff_array');
+
+else
+ load('wdiff_array.mat');
+end
+
+whos wdiff_array
+
+% total number of windows
+N = length(meas_array);
+
+%----------------------------------------------
+% CODE IN PROGRESS
+% 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
+
+[m1,m2,m3,m4,m5] = size(wdiff_array)
+
+if 0==1
+
+ ratio_array = []; ll = 0;
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ % check that there is a measurement for all models
+ if prod(wdiff_array(ii,jj,kk,tt,:)) > 0
+ ll = ll+1;
+ rat = wdiff_array(ii,jj,kk,tt,imodel_max+1) / wdiff_array(ii,jj,kk,tt,imodel_min+1);
+ ratio_array(ll,:) = [tt ii jj kk rat];
+ end
+ end
+ end
+ end
+ end
+
+ %ratio_sort = sortrows(ratio_array,5);
+ ratio_sort = sortrows(ratio_array,[-1 5]);
+
+ hmax = length(ratio_sort);
+ hmax = 100;
+ for hh = 1:hmax
+ 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
+
+ %------------------------
+
+ rec_array = []; ll = 0;
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+
+ % check that there is a measurement for all components for all models
+ if prod( wdiff_array(ii,jj,:,tt,:) ) > 0
+ ll = ll+1;
+ rat = prod( wdiff_array(ii,jj,:,tt,imodel_max+1) ./ wdiff_array(ii,jj,:,tt,imodel_min+1) );
+ rec_array(ll,:) = [tt ii jj rat];
+ end
+ end
+
+ end
+ end
+ rec_sort = sortrows(rec_array,[-1 4]);
+
+ for hh = 1:length(rec_sort)
+ tt = rec_sort(hh,1);
+ ii = rec_sort(hh,2);
+ jj = rec_sort(hh,3);
+ disp(sprintf('%10s --> %6.1f %7s %5s %10.6f',char(eids{ii}),...
+ Tvec(tt),char(strec{jj}),rec_sort(hh,4) ));
+ end
+
+ %------------------------
+
+% % loop over all stations
+% tarray = []; ll = 0;
+% for jj = 1:nrec
+% imatch = find( ratio_array(:,3)==jj);
+% if length(imatch)==6
+% 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
+
+ figure; N = plot_histo(ratio_array(:,end),[0:0.2:5]);
+ figure; N = plot_histo(log(ratio_array(:,end)),[-4:0.2:4]);
+
+ break
+end
+
+%=======================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2008-08-07 00:53:44 UTC (rev 12572)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2008-08-07 05:26:54 UTC (rev 12573)
@@ -1,6 +1,6 @@
%
% compute_misfit.m
-% CARL TAPE, 19-June-2008
+% CARL TAPE, 06-Aug-2008
% printed xxx
%
% This script reads in a set of window_chi output measurement files from
@@ -11,6 +11,8 @@
% The output files are stored within the directories generated by
% make_dirs.pl, which must be run first.
%
+% See also compare_misfit.m for comparing misfit values for multiple models.
+%
% calls read_window_chi.m, plot_histo.m, readCMT.m, read_station_SPECFEM.m
% called by xxx
%
@@ -22,7 +24,6 @@
format compact
dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
-ftag = '_window_chi';
Tvec = [6 2];
% LOWEST period used for each measurement band-pass (upper is 40 s)
@@ -33,7 +34,7 @@
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
+ % integrating the waveforms in mt-measure_adj.f90
%-------------------------
% USER INPUT
@@ -57,8 +58,11 @@
iwrite = input(' Enter 1 to write files (0 otherwise) : ');
% files: event IDs, CMTSOLUTIONS, ststions
-eid_file = ['/net/sierra/raid1/carltape/results/KERNELS/kernel_' stmod '/kernels_' stmod];
-cmt_file_all = '/net/sierra/raid1/carltape/results/SOURCES/socal_7/SOCAL_FINAL_CMT_v7';
+dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_9/';
+eid_file = [dir_source 'EIDs_only_loc'];
+eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_' stmod];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v9'];
stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
%-------------------------
@@ -66,13 +70,50 @@
% load the event IDs corresponding to the kernels
% load these as strings, since event IDs could have letters
-eids = textread(eid_file,'%s');
+eids = textread(eid_file,'%s');
+eids_num = textread(eid_file,'%f'); % numeric version
+%eids_num = str2num(char(eids));
nevent = length(eids);
-eids_num = str2num(char(eids)); % numeric version
% read in CMT solutions
-[date,tshift,hdur,slat,slon,dep,M,eid,elabel] = readCMT(cmt_file_all, 13, 1);
+[date,tshift,hdur,slat,slon,dep,M,eid_cmt,elabel] = readCMT(cmt_file_all, 13, 1);
+% get the indices of these events into the full list
+ievent_full = zeros(nevent,1);
+for ii=1:nevent
+ itemp = strmatch(eids(ii),eid_cmt);
+ if ~isempty(itemp)
+ ievent_full(ii) = itemp;
+ else
+ error(['event ' eids{ii} ' is not on the master list']);
+ end
+end
+disp('all input EIDs matched to the master list');
+
+% if 0==1
+% fid = fopen([odir ftag '_window_picks_by_' stks{kk}],'w');
+% fprintf(fid,'%12s%5s%10s%10s%10s\n',stks{kk},' ',...
+% ['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
+% fprintf(fid,'%12s%5s%10i%10i%10i\n','TOTAL -->',' ',ntot_1,sum(ntot_1));
+% for ii = 1:nw
+% jj = isort(ii);
+% fprintf(fid,'%12s%5i%10i%10i%10i\n',labs{jj},jj,nwin_out(jj,:),ntot_2(jj));
+% end
+% fclose(fid);
+%
+% fid = fopen('eids_reduced','w');
+% for ii=1:nevent
+% imatch = strmatch(eids(ii),eid_cmt);
+% slon_match = slon(imatch);
+% slat_match = slat(imatch);
+% if slon_match <= -117.2
+% fprintf(fid,'%s\n',eids{ii});
+% end
+% end
+% fclose(fid);
+% break
+% end
+
%-----------------------
% read in stations file (used for synthetics)
@@ -123,14 +164,14 @@
for tt = 1:nper
Tper = Tvec(tt);
sTper = sprintf('%3.3i',Tper);
- %if imod <= 2, sTper = sprintf('%2.2i',Tper); end % CHT
+ if imod <= 2, sTper2 = sprintf('%2.2i',Tper); else sTper2 = sTper; end % CHT
disp('-----------------------------------------');
for ii = imin:imax % loop over events
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'];
+ mfile = [eids{ii} '_T' sTper2 '_window_chi'];
filename = [dir1 mfile];
if ~exist(filename,'file')
@@ -178,7 +219,7 @@
kinds = [k1:k2]';
meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin ...
- iker measCC_dT sigmaCC_dT measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
+ iker measCC_dT sigmaCC_dT measCC_dA sigmaCC_dA measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
% 1 kinds
% 2 index_T
% 3 index_event
@@ -189,16 +230,17 @@
% 8 iker
% 9 measCC_dT
% 10 sigmaCC_dT
- % 11 measMT_dT
- % 12 sigmaMT_dT
- % 13 win_chi
- % 14 rec_chi
- % 15 tr_chi
+ % 11 measCC_dA
+ % 12 sigmaCC_dA
+ % 13 measMT_dT
+ % 14 sigmaMT_dT
+ % 15 win_chi
+ % 16 rec_chi
+ % 17 tr_chi
end
end
end
-
save('meas_array.mat','meas_array');
else
@@ -211,80 +253,6 @@
N = length(meas_array);
%----------------------------------------------
-% CODE IN PROGRESS
-% 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
-
-if 0==1
- 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
-
- 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
-end
-
-%----------------------------------------------
% check some of the output window rows
if 0==1
itestvec = [1 2];
@@ -293,9 +261,9 @@
itest = itestvec(ii);
meas_array(itest,:)
display_meas(meas_array(itest,:),Tvec,eids,strec,comps);
- meas_array(itest,12)
+ meas_array(itest,14)
0.5 * meas_array(itest,9)^2 / meas_array(itest,10)^2
- 0.5 * meas_array(itest,11)^2 / meas_array(itest,10)^2
+ 0.5 * meas_array(itest,13)^2 / meas_array(itest,10)^2
end
end
@@ -334,6 +302,18 @@
end
end
+%----------------------------------------------
+% tally stations by event (and period)
+nrec_all_event = zeros(nevent,3);
+for ii = 1:nevent
+ imatch1 = find( and( meas_array(:,3)==ii, meas_array(:,2)==1) );
+ imatch2 = find( and( meas_array(:,3)==ii, meas_array(:,2)==2) );
+ imatch = find(meas_array(:,3)==ii);
+ nrec_all_event(ii,1) = length(unique(meas_array(imatch1,5)));
+ nrec_all_event(ii,2) = length(unique(meas_array(imatch2,5)));
+ nrec_all_event(ii,3) = length(unique(meas_array(imatch,5)));
+end
+
if iwrite == 1
stks = {'event','network','receiver','component'};
for kk = 1:length(stks)
@@ -359,6 +339,18 @@
end
fclose(fid);
end
+
+ % write out the number of stations for each event
+ % THIS IS ALSO INCLUDED IN THE data_norms_sort TEXT FILE BELOW
+ [junk, isort] = sortrows([[1:nevent]' nrec_all_event],-4);
+ fid = fopen([odir ftag '_receiver_picks_by_event'],'w');
+ fprintf(fid,'%12s%5s%10s%10s%10s\n','event',' ',...
+ ['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
+ for ii = 1:nevent
+ jj = isort(ii);
+ fprintf(fid,'%12s%5i%10i%10i%10i\n',eids{jj},jj,nrec_all_event(jj,:));
+ end
+ fclose(fid);
end
%----------------------------------------------
@@ -389,7 +381,7 @@
dnorm_sq = zeros(nevent,1);
for ii = 1:nevent
imatch = find( meas_array(:,3)==ii );
- Sval = meas_array(imatch,15);
+ Sval = meas_array(imatch,17);
if sum(isnan(Sval)) > 0
disp(sprintf('i = %i, eid %s',ii,eids{ii}));
error('Sval has at least one NaN data entry');
@@ -410,27 +402,86 @@
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([odir ftag '_data_norms_sort.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),'--','--');
+ data_norms = [eids_num Ns nrec_all_event dnorm_sq*nevent dnorm_sq dnorm ws];
+
+ % list only records that fit criteria
+ if 0==1
+ eid_reject = load('/net/sierra/raid1/carltape/results/KERNELS/kernels_exclude');
+ X1 = 20;
+ X2 = 40;
+ D1 = 3;
+ isave = find( or(nrec_all_event(:,3) > X2, and(nrec_all_event(:,3) > X1, dnorm_sq*nevent > D1)) );
+ [eids_save,ijunk] = setdiff( eids_num(isave), eid_reject );
+ else
+ eids_save = eids_num;
+ end
+
+ klabs = {'dnorm','nrec'};
+ kind = [-6 -5];
+ for kk=1:2
+ dsort = sortrows(data_norms,kind(kk));
+ fid = fopen([odir ftag '_data_norms_sort_by_' klabs{kk} '.txt'],'w');
+ fprintf(fid,'%12s%8s%6s%6s%6s%12s%12s%12s%12s\n','eid','Nwin',...
+ ['Nrec' num2str(Tvec(1))],['Nrec' num2str(Tvec(2))],'Nrec','dnorm2*E','dnorm2','dnorm','weight');
+ fprintf(fid,'%12s%8i%6s%6s%6s%12.2f%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq*nevent),sum(dnorm_sq),'--','--');
+ for ii = 1:nevent
+ if any(dsort(ii,1)==eids_save)
+ fprintf(fid,'%12i%8i%6i%6i%6i%12.4f%12.4f%12.4f%12.4f\n',dsort(ii,:));
+ end
+ end
+ %fprintf(fid,'%12s%8i%6s%6s%6s%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq),'--','--');
+ fclose(fid);
+ end
+
+ fid1 = fopen([odir ftag '_eid_save.txt'],'w');
+ fid2 = fopen([odir ftag '_eid_exclude.txt'],'w');
for ii = 1:nevent
- fprintf(fid,'%12i%8i%12.4f%12.4f%12.4f\n',dsort(ii,:));
+ if any(eids_num(ii)==eids_save)
+ fprintf(fid1,'%i\n',eids_num(ii));
+ else
+ fprintf(fid2,'%i\n',eids_num(ii));
+ end
end
- fprintf(fid,'%12s%8i%12.4f%12s%12s\n','TOTAL',sum(Ns),sum(dnorm_sq),'--','--');
- fclose(fid);
+ fclose(fid2);
+
+ %km04 = load('/net/sierra/raid1/carltape/results/KERNELS/kernel_m04/kernels_m04');
+ %km05 = load('OUTPUT/m05/window/m05_window_eid_save.txt');
+ %eid_5_not_in_4 = setdiff(km05,km04)
+ %eid_4_not_in_5 = setdiff(km04,km05)
% save variables (for wave2d_subspace_3D.m)
save([odir ftag '_data_norms'],'idatacov','dcov_fac_e',...
- 'eids','eids_num','Ns','dnorm','ws','dmisfit','N');
+ 'eids','eids_num','Ns','dnorm','ws','dmisfit','N','nrec_all_event');
+
+ % write out all stations that have at least one measurement (for GMT)
+ iuse = find( sum(nwin_all_rec,2) > 0);
+ write_station_SPECFEM([odir 'STATIONS_' stmod],rlon(iuse),rlat(iuse),...
+ relev(iuse),rburial(iuse),stnm(iuse),netwk(iuse));
+
+ % write out a psmeca file for all events used (for GMT)
+ iuse = find( sum(nwin_all_event,2) > 0);
+ inds = ievent_full(iuse);
+ writeCMT_psmeca([odir 'socal_tomo_' stmod],date(inds),slat(inds),slon(inds),...
+ dep(inds),M(inds,:),eid_cmt(inds));
end
break
+% list the events that use at least 15 stations
+fid = fopen([odir ftag '_eid_good.txt'],'w');
+for ii = 1:nevent
+ if nrec_all_event(ii,3) > 15
+ fprintf(fid,'%i\n',eids_num(ii));
+ disp([eids_num(ii) nrec_all_event(ii,:)])
+ end
+end
+fclose(fid);
+
+break
+
%-------------------------------
+% find records that meet particular criteria
-% 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) )
@@ -439,31 +490,20 @@
meas_disp = meas_array(icheck(isort),:);
display_meas(meas_disp,Tvec,eids,strec,comps);
-% find records that meet particular criteria
+DA_MIN = -0.8;
+DA_MAX = 3.0;
+icheck = find( or( meas_array(:,11) <= DA_MIN, meas_array(:,11) >= DA_MAX) );
+%[junk, isort] = sortrows( meas_array(icheck,:), 11 );
+[junk, isort] = sortrows( meas_array(icheck,:), 3 );
+meas_disp = meas_array(icheck(isort),:);
+display_meas(meas_disp,Tvec,eids,strec,comps);
+
CHI_MIN = 100;
DT_SIGMA_MIN = 0.19;
-icheck = find( and( meas_array(:,15) >= CHI_MIN, meas_array(:,10) < DT_SIGMA_MIN) );
-[junk, isort] = sortrows( meas_array(icheck,:), -13 );
+icheck = find( and( meas_array(:,17) >= CHI_MIN, meas_array(:,10) < DT_SIGMA_MIN) );
+[junk, isort] = sortrows( meas_array(icheck,:), -15 );
meas_disp = meas_array(icheck(isort),:);
display_meas(meas_disp,Tvec,eids,strec,comps);
ebads = unique( eids_num(meas_array(icheck,3)) )
-break
-
-% tabulate results by different categories
-% (1) shortest period in bandpass
-% (2) event
-% (3) receiver
-% (4) component
-% (5) type of measurement (MT or CC)
-
-disp(' ');
-disp(['Tabulating results for ' num2str(nevent) ' selected events']);
-
-disp(' Binning by bandpass group:');
-for ii = 1:nper
- inds = find( meas_array(:,2) == ii );
- disp(sprintf(' %i out of %i -- Tshort = %.2f s -- %i windows',ii,nper,Tvec(ii),length(inds)));
-end
-
%=======================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m 2008-08-07 00:53:44 UTC (rev 12572)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m 2008-08-07 05:26:54 UTC (rev 12573)
@@ -13,7 +13,7 @@
X = meas_array;
[m,n] = size(X);
-if n ~= 15, error(' should be 15 columns in the input matrix'); end
+if n ~= 17, error(' should be 17 columns in the input matrix'); end
% 1 kinds
% 2 index_T
@@ -25,11 +25,13 @@
% 8 iker
% 9 measCC_dT
% 10 sigmaCC_dT
-% 11 measMT_dT
-% 12 sigmaMT_dT
-% 13 win_chi
-% 14 rec_chi
-% 15 tr_chi
+% 11 measCC_dA
+% 12 sigmaCC_dA
+% 13 measMT_dT
+% 14 sigmaMT_dT
+% 15 win_chi
+% 16 rec_chi
+% 17 tr_chi
index_win = X(:,1);
index_per = X(:,2);
@@ -41,18 +43,21 @@
iker = X(:,8);
measCC_dT = X(:,9);
sigmaCC_dT = X(:,10);
-measMT_dT = X(:,11);
-sigmaMT_dT = X(:,12);
-win_chi = X(:,13);
-rec_chi = X(:,14);
-tr_chi = X(:,15);
+measCC_dA = X(:,11);
+sigmaCC_dA = X(:,12);
+measMT_dT = X(:,13);
+sigmaMT_dT = X(:,14);
+win_chi = X(:,15);
+rec_chi = X(:,16);
+tr_chi = X(:,17);
disp('----------------');
disp('INDEX (window, period, event, receiver, component)');
for ii = 1:m
- 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),tr_chi(ii),index_win(ii) ));
+ disp(sprintf('%10s --> %6.1f %7s %5s %3i DT = %6.2f +- %6.2f DA = %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),measCC_dA(ii),tr_chi(ii),index_win(ii) ));
end
disp('INDEX (window, period, event, receiver, component)');
disp('====================');
More information about the cig-commits
mailing list