[cig-commits] r14497 - seismo/3D/ADJOINT_TOMO/iterate_adj/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Fri Mar 27 14:01:00 PDT 2009
Author: carltape
Date: 2009-03-27 14:00:59 -0700 (Fri, 27 Mar 2009)
New Revision: 14497
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
Log:
New and updated matlab scripts for misfit analysis of measurement output from ADJOINT_TOMO/measure_adj. Tested using measurement files from socal crustal models m00 and m16.
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 2009-03-27 21:00:59 UTC (rev 14497)
@@ -0,0 +1,576 @@
+%
+% compare_misfit.m
+% CARL TAPE, 06-Aug-2008
+% printed xxx
+%
+%
+% 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/';
+
+iVRlog = 1; % variance reduction formula (1 for VRL; 0 for VR)
+iwdiff = 1; % demoninator for waveform difference (1 for new; 0 for old "standard" version)
+wtag = sprintf('VR%i_wdiff%i',iVRlog,iwdiff);
+
+iwrite = 1;
+isciencehist = 1; % 0 or 1
+odirsci = '/home/carltape/science_paper/latex/figures/scripts/hist_data/';
+
+DTMIN = 1.0; % see description below (1.0 or 0.0 for socal)
+
+idataset = 3; % 1 (tomo), 2 (extra), 3 (simulation)
+dtags = {'tomo','extra','simulation'};
+dtag = dtags{idataset};
+
+%-------------------------
+
+% min and max periods for the different bandpassed datasets
+%Tminvec = [2 3 6]; Tmaxvec = [30 30 30];
+Tminvec = [2]; Tmaxvec = [30];
+
+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(Tminvec);
+stBtag = 'BP';
+for ii=1:nper
+ stBtag = [stBtag num2str(Tminvec(ii))];
+end
+
+% strings for labeling
+sTbp = repmat(cellstr(' '),1,nper);
+for tt = 1:nper
+ %sTbp{tt} = sprintf('[%is,%is]',Tminvec(tt),Tmaxvec(tt));
+ sTbp{tt} = sprintf('%i-%is',Tminvec(tt),Tmaxvec(tt));
+end
+
+%-------------------------
+% USER INPUT
+
+%imod1 = input(' Enter the model 1 (0): ');
+%imod2 = input(' Enter the model 2 (16): ');
+imod1 = 0;
+imod2 = 16;
+stmod = sprintf('m%2.2i',imod2); % put the VR files in the dir for the final model
+
+idatacov = 2;
+%idatacov = input(' Enter idatacov (1--by event, 2--by window, 3--none) : ');
+ % 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
+stdcovs = {'event','window','none'};
+stdcov = stdcovs{idatacov};
+
+ftag = [stmod '_' stdcov];
+odir = ['OUTPUT_SUBSPACE/' stmod '/' stdcov '/'];
+
+%iwrite = input(' Enter 1 to write files (0 otherwise) : ');
+
+% files: CMTSOLUTIONS and stations
+stsrcvers = '16'; % index for the set of sources (NOT the model iteration index)
+dir_source = ['/net/sierra/raid1/carltape/results/SOURCES/socal_' stsrcvers '/'];
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v' stsrcvers];
+stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
+
+% list of event IDs
+%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_run_' stmod];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+
+eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/eids_' dtag];
+
+%-------------------------
+% 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,seid_cmt,elabel] = readCMT(cmt_file_all, 13, 1);
+eid_cmt = str2num(char(seid_cmt));
+
+% get the indices of these events into the full list
+ievent_full = zeros(nevent,1);
+for ii=1:nevent
+ itemp = strmatch(eids(ii),seid_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');
+
+%-----------------------
+% 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);
+
+%--------------------------------------
+% create indexing array
+
+% % create indexing vector and array
+% ninds = nevent*nrec*ncomp*nper;
+% index_array2vec = zeros(nevent,nrec,ncomp,nper);
+% index_vec2array = zeros(ninds,4);
+%
+% nw = 0;
+% for ievent = 1:nevent
+% for irec = 1:nrec
+% for icomp = 1:ncomp
+% for iper = 1:nper
+% nw = nw+1;
+% index_vec2array(nw,:) = [ievent irec icomp iper];
+% index_array2vec(ievent,irec,icomp,iper) = nw;
+% end
+% end
+% end
+% end
+% nwin_tot = nw; % = nevent*nrec*ncomp*nper
+
+%-----------------------
+
+iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
+%iread = 1;
+if iread == 1
+ % range of events
+ imin = 1; imax = nevent; % default
+ %imin = 1; imax = 10;
+ %imin = 139; imax = imin;
+ %imin = 141; imax = imin; % Chino Hills
+
+ % model 1
+ stmod = sprintf('m%2.2i',imod1);
+ meas_array1 = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps);
+ save('meas_array1','meas_array1');
+
+ % model 2
+ stmod = sprintf('m%2.2i',imod2);
+ meas_array2 = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps);
+ save('meas_array2','meas_array2');
+else
+ load(['meas_array1.mat']);
+ load(['meas_array2.mat']);
+end
+
+% check at least that the station indices are identical
+if sum( meas_array1(:,5) - meas_array2(:,5) ) ~= 0
+ error('mismatch of the two window files');
+end
+
+% total number of windows
+N = length(meas_array1);
+
+% combine the columns that you want into one datafile
+% GET RID OF THE 0.5 FACTOR IN tr_chi, SO THAT IS IS DT^2 / SIGMA^2
+meas_all = [meas_array1(:,[1:7 9 10 13 14 11]) ... % first 11 columns
+ meas_array1(:,12) meas_array2(:,12) meas_array1(:,8) meas_array2(:,8) ...
+ meas_array1(:,15) meas_array2(:,15) meas_array1(:,16) meas_array2(:,16) ...
+ meas_array1(:,17) meas_array2(:,17) meas_array1(:,18) meas_array2(:,18) ...
+ meas_array1(:,19) meas_array2(:,19) meas_array1(:,23) meas_array2(:,23) ...
+ meas_array1(:,25) meas_array2(:,25)];
+% 1 kinds
+% 2 index_T
+% 3 index_event
+% 4 index_network
+% 5 index_rec
+% 6 index_comp
+% 7 iwin
+% 8 seisdur
+% 9 windur
+% 10 seisd2
+% 11 wind2
+% 12 T_pmax_dat
+% 13 m00 - T_pmax_syn
+% 14 m16 - T_pmax_syn
+% 15 m00 - iker
+% 16 m16 - iker
+%-------------------
+% waveform difference information
+% 17 m00 - seiss2
+% 18 m16 - seiss2
+% 19 m00 - wins2
+% 20 m16 - wins2
+% 21, m00 - seis_diff2
+% 22, m16 - seis_diff2
+% 23, m00 - win_diff2
+% 24, m16 - win_diff2
+%-------------------
+% 25, m00 - DT-CC
+% 26, m16 - DT-CC
+% 27, m00 - DT-MT
+% 28, m16 - DT-MT
+% 29, m00 - tr_chi
+% 30, m16 - tr_chi
+%-------------------
+% 31, m00 - seis_chi
+% 32, m16 - seis_chi
+% 33, m00 - win_chi
+% 34, m16 - win_chi
+% 35 --> VRseis (with NaN)
+% 36 --> VRwin (with NaN)
+
+% find the unique seismograms (value from first window in each record)
+bseis1 = (meas_all(:,7) == 1);
+nseis = sum(bseis1);
+disp(sprintf('%i out of %i windows are the only window on the seismogram',nseis,N));
+
+for kk = 0:1 % loop over both models
+ % find the number of MT and CC measurements used, and replace MT-DT with the CC value
+ iMT = find(meas_all(:,15+kk) == 1);
+ iCC = find(meas_all(:,15+kk) == 2);
+ nMT = length(iMT);
+ disp(sprintf('%i out of %i windows use multitaper measurment',nMT,N));
+
+ % compare DT for MT vs CC
+ figure; hold on; ymx = 4;
+ plot(meas_all(:,25+kk),meas_all(:,27+kk),'.');
+ plot(meas_all(iCC,25+kk),meas_all(iCC,27+kk),'ro');
+ plot([-1 1]*ymx,[-1 1]*ymx,'r--'); grid on; axis equal, axis([-1 1 -1 1]*ymx);
+ xlabel('Cross-correlation DT'); ylabel('Multitaper DT (when possible)');
+
+ % KEY: replace MT DT zero values with the CC DT values
+ meas_all(iCC,27+kk) = meas_all(iCC,25+kk);
+ plot(meas_all(iCC,25+kk),meas_all(iCC,27+kk),'ko');
+end
+
+%itemp = find(and(meas_all(:,23)==-2.9,meas_all(:,25) > 0))
+%display_meas(meas_array2(itemp,:),Tminvec,eids,strec,comps);
+%break
+
+%----------------------------------------
+% WAVEFORM DIFFERENCE MISFIT FUNCTION
+
+if iwdiff == 1
+ % new version
+ seis_chi1 = meas_all(:,21) ./ sqrt( meas_all(:,10) .* meas_all(:,17) );
+ seis_chi2 = meas_all(:,22) ./ sqrt( meas_all(:,10) .* meas_all(:,18) );
+ win_chi1 = meas_all(:,23) ./ sqrt( meas_all(:,11) .* meas_all(:,17) );
+ win_chi2 = meas_all(:,24) ./ sqrt( meas_all(:,11) .* meas_all(:,18) );
+else
+ % old "standard" version
+ seis_chi1 = meas_all(:,21) ./ meas_all(:,10);
+ seis_chi2 = meas_all(:,22) ./ meas_all(:,10);
+ win_chi1 = meas_all(:,23) ./ meas_all(:,11);
+ win_chi2 = meas_all(:,24) ./ meas_all(:,11);
+end
+
+% add misfit to meas_all
+meas_all = [meas_all seis_chi1 seis_chi2 win_chi1 win_chi2];
+
+%----------------------------------------
+% VARIANCE REDUCTION
+
+if iVRlog == 1
+ VRseis = log( seis_chi1 ./ seis_chi2 );
+ VRwin = log( win_chi1 ./ win_chi2 );
+ edges_vr = [-4:0.5:4]; ylims = [0 0.4];
+ VRtag = 'VRL';
+else
+ VRseis = (seis_chi1 - seis_chi2) ./ seis_chi1;
+ VRwin = (win_chi1 - win_chi2) ./ win_chi1;
+ edges_vr = [-4:0.25:1]; ylims = [0 0.4];
+ VRtag = 'VR';
+end
+
+
+% Remove windows that are have time shifts less than DTMIN both in the
+% initial model and in the final model -- we are not interested in the VR
+% for these records. But we are interested if the time shift gets WORSE
+% than DTMIN in going from m00 to m16.
+%DTMIN = 0.0;
+DTMIN_tag = sprintf('DTMIN_%3.3ims',DTMIN*1000);
+bbad = and( abs(meas_all(:,25)) < DTMIN, abs(meas_all(:,26)) < DTMIN );
+bgood = ~bbad;
+ngood = sum(bgood);
+VRwin(bbad) = NaN;
+
+if DTMIN == 0
+ VRseis(~bseis1) = NaN; % set VR for non-unique entries to NaN
+else
+ disp('excluding seismograms that have windows fit before and after');
+ if 0==1
+ % exclude seismograms in which ANY window has a time shifts less than DTMIN
+ % both in the initial model and in the final model
+ VRseis(~bseis1) = NaN; % set VR for non-unique entries to NaN
+ VRseis(bbad) = NaN; % set VR to NaN if ANY window is NaN
+ else
+ % exclude seismograms in which ALL windows have time shifts less than DTMIN
+ % both in the initial model and in the final model -- THIS IS SLOW
+
+ %for ii = 1:N % loop over all windows
+ ind1=3; ind2=5; ind3=6; ind4=2;
+ for ii = 1:N
+ %for ii = 10:10
+ if mod(ii,1000) == 0, disp(sprintf('%i out of %i',ii,N)); end
+ sind = meas_all(ii,[ind1 ind2 ind3 ind4]); % index into record
+
+ % find the other windows on the seismogram and assign NaN if all windows
+ % on the seismogram have VRwin NaN
+ imatch = find( and( and(meas_all(:,ind1)==sind(1),meas_all(:,ind2)==sind(2)),...
+ and(meas_all(:,ind3)==sind(3),meas_all(:,ind4)==sind(4)) ) == 1);
+ inans = isnan(VRwin(imatch));
+ if all(inans) % all windows are bad
+ VRseis(imatch) = NaN;
+ else % at least one window is good
+ ikeeper = imatch(~inans);
+ VRseis(setdiff(imatch,ikeeper(1))) = NaN;
+ end
+ end
+ end
+end
+
+%temp = [VRseis VRwin]; disp(temp(1:30,:)), sum(~isnan(VRseis)), sum(~isnan(VRwin))
+
+% good seismograms, unrepeated
+bgoodseis = ~isnan(VRseis);
+ngoodseis = sum(bgoodseis);
+
+disp(sprintf('%i out of %i windows have DT > %.2f s for m%2.2i or m%2.2i',ngood,N,DTMIN,imod1,imod2));
+disp(sprintf('%i out of %i seismograms have DT > %.2f s for m%2.2i or m%2.2i',ngoodseis,nseis,DTMIN,imod1,imod2));
+
+% Remove windows that are worse within records that are better
+%bbad2 = and( VRseis > 0, VRwin < 0); length(bbad2)
+%display_meas(meas_array1(bbad2==1,:),Tminvec,eids,strec,comps);
+%break
+
+% fraction of records and windows that have a worse waveform difference
+fworse_seis = length(find(VRseis(bgoodseis) < 0)) / ngoodseis;
+fworse_win = length(find(VRwin(bgood) < 0)) / ngood;
+
+% add VR to meas_all
+meas_all = [meas_all VRseis VRwin];
+
+VRvals = [mean(VRwin(bgood)) median(VRwin(bgood)) mean(VRseis(bgoodseis)) median(VRseis(bgoodseis))];
+disp(sprintf('Variance reduction values for the %i windows:',ngood));
+disp(sprintf('%12s%12s%12s%12s','VRwin-mean','VRwin-med','VRseis-mean','VRseis-med'));
+disp(sprintf('%12.4f%12.4f%12.4f%12.4f',VRvals));
+
+figure; nr=3; nc=2;
+edges_chi = [0:0.25:4];
+
+subplot(nr,nc,1);
+x = win_chi1(bgood);
+%x = log(win_chi1(bgood));
+Nbin = plot_histo(x,edges_chi);
+xlabel('chi2 for each window'); ylim(ylims); grid on;
+title(sprintf('chi2 for m%2.2i WINDOWS: mean %.3f, median %.3f',imod1,mean(x),median(x)));
+
+subplot(nr,nc,2);
+x = seis_chi1(bgoodseis);
+%x = log(seis_chi1(bgoodseis));
+Nbin = plot_histo(x,edges_chi);
+xlabel('chi2 for each seismogram'); ylim(ylims); grid on;
+title(sprintf('chi2 for m%2.2i SEIS: mean %.3f median %.3f',imod1,mean(x),median(x)));
+
+subplot(nr,nc,3);
+x = win_chi2(bgood);
+%x = log(win_chi2(bgood));
+Nbin = plot_histo(x,edges_chi);
+xlabel('chi2 for each window'); ylim(ylims); grid on;
+title(sprintf('chi2 for m%2.2i WIN: mean %.3f median %.3f',imod2,mean(x),median(x)));
+
+subplot(nr,nc,4);
+x = seis_chi2(bgoodseis);
+%x = log(seis_chi2(bgoodseis));
+Nbin = plot_histo(x,edges_chi);
+xlabel('chi2 for each seismogram'); ylim(ylims); grid on;
+title(sprintf('chi2 for m%2.2i SEIS: mean %.3f median %.3f',imod2,mean(x),median(x)));
+
+subplot(nr,nc,5);
+Nbin = plot_histo(VRwin(bgood),edges_vr);
+xlabel([VRtag ' for each window']); ylim(ylims); grid on;
+title(sprintf('WIN: %.2f > 0, %smean = %.3f, %smed = %.3f',1-fworse_win,VRtag,VRvals(1),VRtag,VRvals(2)));
+
+subplot(nr,nc,6);
+Nbin = plot_histo(VRseis(bgoodseis),edges_vr);
+xlabel([VRtag ' for each seismogram']); ylim(ylims); grid on;
+title(sprintf('SEIS: %.2f > 0, %smean = %.3f, %smed = %.3f',1-fworse_seis,VRtag,VRvals(3),VRtag,VRvals(4)));
+
+fontsize(10), orient tall, wysiwyg
+
+%---------------------
+% MT-DT and CC-DT for final model
+
+DTvals = zeros(1,4);
+x1 = meas_all(bgoodseis,26);
+x2 = meas_all(bgoodseis,28);
+DTvals = [mean(x1) std(x1) mean(x2) std(x2)];
+
+if 1==1
+ figure; nr=2; nc=1; ylims = [0 0.3];
+ edges_DT = [-3:0.25:3];
+ tlabs = {'CC','MT'};
+
+ for kk=[0 2]
+ % windows
+ subplot(nr,nc,1+kk/2);
+ x = meas_all(bgood,26+kk);
+ Nbin = plot_histo(x,edges_DT);
+ xlabel('DT for each window'); ylim(ylims); grid on;
+ title(sprintf('%s-DT, m%2.2i, %i events: mean %.3f pm %.3f',tlabs{1+kk/2},imod2,nevent,mean(x),std(x)));
+
+% % seismograms
+% subplot(nr,nc,2+kk);
+% x = meas_all(bgoodseis,26+kk);
+% Nbin = plot_histo(x,edges_DT);
+% xlabel('DT for each window'); ylim(ylims); grid on;
+% title(sprintf('%s-DT, m%2.2i, %i events: mean %.3f pm %.3f',tlabs{1+kk/2},imod2,nevent,mean(x),std(x)));
+ end
+ fontsize(10), orient tall, wysiwyg
+end
+
+% %----------------------------------------------
+% % check some of the output window rows
+% if 0==1
+% itestvec = [41];
+% for ii = 1:length(itestvec)
+% disp('====>');
+% itest = itestvec(ii);
+% meas_array1(itest,14)
+% meas_array2(itest,14)
+% display_meas(meas_array1(itest,:),Tminvec,eids,strec,comps);
+% display_meas(meas_array2(itest,:),Tminvec,eids,strec,comps);
+% end
+% end
+
+% add a couple more measures
+% 32 - Pwin, power in window relative to record
+% 33 - Pwin*VRwin, quality of VR for window (will favor single-window records)
+%meas_all_mod = [meas_all meas_all(:,11)./meas_all(:,10) meas_all(:,11)./meas_all(:,10).*meas_all(:,31)];
+meas_all_mod = meas_all;
+
+% if isciencehist==1
+%
+% % 1 - ii
+% % 2 - eid
+% % 3 - ind
+% % 4 - sta
+% % 5 - net
+% % 6 - comp
+% % 7 - iwin
+% % 8 - m16 - DT-CC
+% % 9 - m16 - DT-MT
+% % 10 - m00 - seis_chi
+% % 11 - m16 - seis_chi
+% % 12 - m00 - win_chi
+% % 13 - m16 - win_chi
+%
+% %dsort = meas_all_mod(bgood,:); % extract a subset
+% dsort = meas_all_mod;
+% klabs = {'win','seis'};
+% for kk = 1:2
+% fpre = 'hist_all_wdiff';
+% %fpre = 'hist_T006_T030_DT';
+%
+% fname = [odirsci fpre '_' DTMIN_tag '_' dtag '_' wtag '_' klabs{kk} '.dat'];
+% fid = fopen(fname,'w');
+% pp = 0;
+% for ii = 1:length(dsort)
+% ilist = 0;
+% if any(kk==[2 4]) % seismogram list
+% %if and(~isnan(dsort(ii,[33 34])),dsort(ii,7)==1), ilist = 1; end
+% if ~isnan(dsort(ii,35)), ilist = 1; end
+% else % window list
+% if ~isnan(dsort(ii,36)), ilist = 1; end
+% end
+%
+% if ilist==1
+% pp = pp+1;
+% fprintf(fid,['%6i%12i%6i%6s%6s%6s%6i' repmat('%8.3f',1,6) '\n'],...
+% pp,eids_num(dsort(ii,3)),dsort(ii,1),...
+% char(stnm{dsort(ii,5)}),char(stnet{dsort(ii,4)}),char(comps{dsort(ii,6)}),...
+% dsort(ii,[7 26 28 29:32]));
+% end
+% end
+% fclose(fid);
+% end
+% end
+
+if iwrite == 1
+
+ if isciencehist==0
+ klabs = {'DT2','DT_chi2','seis_chi1','seis_chi2','win_chi1','win_chi2','VRseis','VRwin'};
+ kind = [-26 30 31 32 33 34 -35 -36]; % columns to sort, corresponding to klabs
+ else
+ odir = odirsci;
+ klabs = {'seis_chi1','win_chi1'};
+ kind = [31 33];
+ end
+ numk = length(klabs);
+
+ disp('writing out sorted files...');
+ for kk = 1:numk
+ %for kk = numk:numk
+ %for kk = 1:1
+ dsort = sortrows(meas_all_mod,kind(kk));
+ fid = fopen([odir ftag '_' stBtag '_' DTMIN_tag '_' wtag '_' dtag '_sort_' klabs{kk} '.txt'],'w');
+ fprintf(fid,['%6s%12s%8s%6s%6s%6s%10s%6s' repmat('%8s',1,12) '\n'],' ','eid','ind',...
+ 'stnm','net','comp','BP','iwin','DTCC1','DTCC2','DTMT1','DTMT2',...
+ 'DTchi1','DTchi2','Fseis1','Fseis2','Fwin1','Fwin2','VRseis','VRwin');
+ pp = 0;
+ for ii = 1:N
+ ilist = 0;
+ if any(abs(kind(kk))==[31 32 35]) % seismogram list
+ %if and(~isnan(dsort(ii,[33 34])),dsort(ii,7)==1), ilist = 1; end
+ if ~isnan(dsort(ii,35)), ilist = 1; end
+ else % window list
+ if ~isnan(dsort(ii,36)), ilist = 1; end
+ end
+ if ilist == 1
+ pp = pp+1;
+ fprintf(fid,['%6i%12i%8i%6s%6s%6s%10s%6i' repmat('%8.3f',1,12) '\n'],...
+ pp,eids_num(dsort(ii,3)),...
+ dsort(ii,1),...
+ char(stnm{dsort(ii,5)}),...
+ char(stnet{dsort(ii,4)}),...
+ char(comps{dsort(ii,6)}),...
+ char(sTbp{dsort(ii,2)}),...
+ dsort(ii,[7 25:36]));
+ end
+ end
+ %fprintf(fid,'%12s%8i%6s%6s%6s%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq),'--','--');
+ fclose(fid);
+ end
+
+ % values to go with the histograms
+ % 1 - m16 DTCC mean
+ % 2 - m16 DTCC std
+ % 3 - m16 DTMT mean
+ % 4 - m16 DTMT std
+ % 5 - VRwin mean
+ % 6 - VRwin std
+ % 7 - VRseis mean
+ % 8 - VRseis mean
+ vtemp = [DTvals VRvals];
+ disp('writing out statistic values...');
+ fid = fopen([odir ftag '_' stBtag '_' DTMIN_tag '_' wtag '_' dtag '_values.txt'],'w');
+ fprintf(fid,[repmat('%12.3e',1,length(vtemp)) '\n'],vtemp);
+ fclose(fid);
+end
+
+%=======================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2009-03-27 20:51:08 UTC (rev 14496)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2009-03-27 21:00:59 UTC (rev 14497)
@@ -28,6 +28,7 @@
% min and max periods for the different bandpassed datasets
%Tminvec = [2 6]; Tmaxvec = [30 30];
Tminvec = [2 3 6]; Tmaxvec = [30 30 30];
+%Tminvec = [6]; Tmaxvec = [30];
comps = {'BHZ','BHR','BHT'};
% regardless of the component label for the DATA, the measurement code
@@ -63,14 +64,22 @@
iwrite = input(' Enter 1 to write files (0 otherwise) : ');
-% files: event IDs, CMTSOLUTIONS, stations
-dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_11/';
-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_v11'];
+% files: CMTSOLUTIONS and stations
+stsrcvers = '16'; % index for the set of sources (NOT the model iteration index)
+dir_source = ['/net/sierra/raid1/carltape/results/SOURCES/socal_' stsrcvers '/'];
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v' stsrcvers];
stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
+% list of event IDs
+%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_run_' stmod];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+
+eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/eids_simulation';
+%eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/eids_tomo';
+%eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/eids_extra';
+
%-------------------------
% read in list of event IDs and sources
@@ -82,12 +91,13 @@
nevent = length(eids);
% read in CMT solutions
-[date,tshift,hdur,slat,slon,dep,M,eid_cmt,elabel] = readCMT(cmt_file_all, 13, 1);
+[date,tshift,hdur,slat,slon,dep,M,seid_cmt,elabel] = readCMT(cmt_file_all, 13, 1);
+eid_cmt = str2num(char(seid_cmt));
% 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);
+ itemp = strmatch(eids(ii),seid_cmt);
if ~isempty(itemp)
ievent_full(ii) = itemp;
else
@@ -109,7 +119,7 @@
%
% fid = fopen('eids_reduced','w');
% for ii=1:nevent
-% imatch = strmatch(eids(ii),eid_cmt);
+% imatch = strmatch(eids(ii),seid_cmt);
% slon_match = slon(imatch);
% slat_match = slat(imatch);
% if slon_match <= -117.2
@@ -162,93 +172,10 @@
% range of events
imin = 1; imax = nevent; % default
- %imin = 72; imax = imin;
+ %imin = 139; imax = imin;
- k2 = 0;
- meas_array = [];
- %for tt = 1:1
- for tt = 1:nper
- Tmin = Tminvec(tt);
- Tmax = Tmaxvec(tt);
- Ttag = sprintf('T%3.3i_T%3.3i',Tmin,Tmax);
- Ttag2 = [Ttag '_' stmod]
- %Ttag = sprintf('T%3.3i',Tmin); % old version
- 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_' Ttag '/'];
- mfile = [eids{ii} '_' Ttag2 '_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
+ meas_array = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps);
+
save('meas_array.mat','meas_array');
else
@@ -260,56 +187,115 @@
% total number of windows
N = length(meas_array);
+% a waveform measure of the power in a window relative to the full record
+seisd2 = meas_array(:,13);
+wind2 = meas_array(:,14);
+seiss2 = meas_array(:,15);
+wins2 = meas_array(:,16);
+win_reld2 = wind2 ./ seisd2;
+win_rels2 = wins2 ./ seiss2;
+
%----------------------------------------------
% check some of the output window rows
if 0==1
- itestvec = [1 2];
+ itestvec = [1:10];
for ii = 1:length(itestvec)
disp('====>');
itest = itestvec(ii);
meas_array(itest,:)
display_meas(meas_array(itest,:),Tminvec,eids,strec,comps);
- meas_array(itest,14)
- 0.5 * meas_array(itest,9)^2 / meas_array(itest,10)^2
- 0.5 * meas_array(itest,13)^2 / meas_array(itest,10)^2
+ meas_array(itest,25)
+ 0.5 * meas_array(itest,19)^2 / meas_array(itest,20)^2
+ 0.5 * meas_array(itest,23)^2 / meas_array(itest,20)^2
end
end
%----------------------------------------------
-% tally windows by:
-% (1) event
-% (2) network
-% (3) station
-% (4) component
+% tally the following
+% (1) number of window picks
+% (2) window durations
+% (3) number of seismograms used
+% by these indices
+% (1) event
+% (2) network
+% (3) station
+% (4) component
nwin_all_event = zeros(nevent,nper);
+durwin_all_event = nwin_all_event;
+nseis_all_event = nwin_all_event;
for tt = 1:nper
for ii = 1:nevent
imatch = find( and( meas_array(:,3)==ii, meas_array(:,2)==tt) );
nwin_all_event(ii,tt) = length(imatch);
+ durwin_all_event(ii,tt) = sum(meas_array(imatch,10));
+
+ % pick the first window index only to count the records
+ imatch = find( and( meas_array(:,7)==1, and( meas_array(:,3)==ii, meas_array(:,2)==tt) ));
+ nseis_all_event(ii,tt) = length(imatch);
end
end
+
nwin_all_net = zeros(nnet,nper);
+durwin_all_net = nwin_all_net;
+nseis_all_net = nwin_all_net;
for tt = 1:nper
for ii = 1:nnet
imatch = find( and( meas_array(:,4)==ii, meas_array(:,2)==tt) );
nwin_all_net(ii,tt) = length(imatch);
+ durwin_all_net(ii,tt) = sum(meas_array(imatch,10));
+
+ % pick the first window index only to count the records
+ imatch = find( and( meas_array(:,7)==1, and( meas_array(:,4)==ii, meas_array(:,2)==tt) ));
+ nseis_all_net(ii,tt) = length(imatch);
end
end
+
nwin_all_rec = zeros(nrec,nper);
+durwin_all_rec = nwin_all_rec;
+nseis_all_rec = nwin_all_rec;
for tt = 1:nper
for ii = 1:nrec
imatch = find( and( meas_array(:,5)==ii, meas_array(:,2)==tt) );
nwin_all_rec(ii,tt) = length(imatch);
+ durwin_all_rec(ii,tt) = sum(meas_array(imatch,10));
+
+ % pick the first window index only to count the records
+ imatch = find( and( meas_array(:,7)==1, and( meas_array(:,5)==ii, meas_array(:,2)==tt) ));
+ nseis_all_rec(ii,tt) = length(imatch);
end
end
+
nwin_all_comp = zeros(ncomp,nper);
+durwin_all_comp = nwin_all_comp;
+nseis_all_comp = nwin_all_comp;
for tt = 1:nper
for ii = 1:ncomp
imatch = find( and( meas_array(:,6)==ii, meas_array(:,2)==tt) );
nwin_all_comp(ii,tt) = length(imatch);
+ durwin_all_comp(ii,tt) = sum(meas_array(imatch,10));
+
+ % pick the first window index only to count the records
+ imatch = find( and( meas_array(:,7)==1, and( meas_array(:,6)==ii, meas_array(:,2)==tt) ));
+ nseis_all_comp(ii,tt) = length(imatch);
end
end
+% quick computation of unique records
+nunique_vec = zeros(9,1);
+[junk,iupath] = unique([meas_array(:,2)],'rows'); nunique_vec(1) = length(iupath); % bandpasses
+[junk,iupath] = unique([meas_array(:,6)],'rows'); nunique_vec(2) = length(iupath); % components
+[junk,iupath] = unique([meas_array(:,4)],'rows'); nunique_vec(3) = length(iupath); % networks
+[junk,iupath] = unique([meas_array(:,3)],'rows'); nunique_vec(4) = length(iupath); % events
+[junk,iupath] = unique([meas_array(:,[4 5])],'rows'); nunique_vec(5) = length(iupath); % receivers
+[junk,iupath] = unique([meas_array(:,[3:5])],'rows'); nunique_vec(6) = length(iupath); % paths
+[junk,iupath] = unique([meas_array(:,[3:6])],'rows'); nunique_vec(7) = length(iupath); % seismograms
+nunique_vec(8) = length(find(meas_array(:,7)==1)); % bandpassed seismograms
+nunique_vec(9) = length(meas_array); % total windows
+
+fid = fopen([odir ftag '_nunique'],'w');
+for ii = 1:length(nunique_vec), fprintf(fid,'%i\n',nunique_vec(ii)); end;
+fclose(fid);
+
%----------------------------------------------
% tally stations for each event (and period)
nrec_all_event = zeros(nevent,nper+1);
@@ -325,44 +311,91 @@
end
if iwrite == 1
- sTfmti = repmat('%8i',1,nper);
- sTfmts = repmat('%8s',1,nper);
+ sTfmti = repmat('%6i',1,nper);
+ sTfmts = repmat('%6s',1,nper);
+ sTfmtf = repmat('%12.1f',1,nper);
+ sTfmts2 = repmat('%12s',1,nper);
sTdash = repmat(cellstr('--'),1,nper);
stks = {'event','network','receiver','component'};
for kk = 1:length(stks)
switch kk
- case 1, nwin_out = nwin_all_event; labs = eids;
- case 2, nwin_out = nwin_all_net; labs = stnet;
- case 3, nwin_out = nwin_all_rec; labs = strec;
- case 4, nwin_out = nwin_all_comp; labs = comps;
+ case 1, nwin_out = nwin_all_event; durwin_out = durwin_all_event; nseis_out = nseis_all_event; labs = eids;
+ case 2, nwin_out = nwin_all_net; durwin_out = durwin_all_net; nseis_out = nseis_all_net; labs = stnet;
+ case 3, nwin_out = nwin_all_rec; durwin_out = durwin_all_rec; nseis_out = nseis_all_rec; labs = strec;
+ case 4, nwin_out = nwin_all_comp; durwin_out = durwin_all_comp; nseis_out = nseis_all_comp; labs = comps;
end
+
+ %----------------------------
nw = length(nwin_out);
ntot_1 = sum(nwin_out, 1);
ntot_2 = sum(nwin_out, 2);
[junk, isort] = sortrows([[1:nw]' nwin_out ntot_2],-[2+nper]);
-
- % write out the file
+
fid = fopen([odir ftag '_window_picks_by_' stks{kk}],'w');
- fprintf(fid,['%12s%5s' sTfmts '%10s\n'],stks{kk},' ',sTbp{:},'TOTAL');
- fprintf(fid,['%12s%5s' sTfmti '%10i\n'],'TOTAL -->',' ',ntot_1,sum(ntot_1));
+ fprintf(fid,['%5s%12s%5s' sTfmts '%10s\n'],' ',stks{kk},' ',sTbp{:},'TOTAL');
+ fprintf(fid,['%5s%12s%5s' sTfmti '%10i\n'],' ','TOTAL -->',' ',ntot_1,sum(ntot_1));
for ii = 1:nw
jj = isort(ii);
- fprintf(fid,['%12s%5i' sTfmti '%10i\n'],labs{jj},jj,nwin_out(jj,:),ntot_2(jj));
+ fprintf(fid,['%5i%12s%5i' sTfmti '%10i\n'],ii,labs{jj},jj,nwin_out(jj,:),ntot_2(jj));
end
fclose(fid);
+
+ %----------------------------
+ ntot_1 = sum(durwin_out, 1);
+ ntot_2 = sum(durwin_out, 2);
+ [junk, isort] = sortrows([[1:nw]' durwin_out ntot_2],-[2+nper]);
+
+ fid = fopen([odir ftag '_window_duration_by_' stks{kk}],'w');
+ fprintf(fid,['%5s%12s%5s' sTfmts2 '%14s\n'],' ',stks{kk},' ',sTbp{:},'TOTAL');
+ fprintf(fid,['%5s%12s%5s' sTfmtf '%14.1f\n'],' ','TOTAL -->',' ',ntot_1,sum(ntot_1));
+ for ii = 1:nw
+ jj = isort(ii);
+ fprintf(fid,['%5i%12s%5i' sTfmtf '%14.1f\n'],ii,labs{jj},jj,durwin_out(jj,:),ntot_2(jj));
+ end
+ fclose(fid);
+
+ %----------------------------
+ ntot_1 = sum(nseis_out, 1);
+ ntot_2 = sum(nseis_out, 2);
+ [junk, isort] = sortrows([[1:nw]' nseis_out ntot_2],-[2+nper]);
+
+ fid = fopen([odir ftag '_seismograms_used_by_' stks{kk}],'w');
+ fprintf(fid,['%5s%12s%5s' sTfmts '%10s\n'],' ',stks{kk},' ',sTbp{:},'TOTAL');
+ fprintf(fid,['%5s%12s%5s' sTfmti '%10i\n'],' ','TOTAL -->',' ',ntot_1,sum(ntot_1));
+ for ii = 1:nw
+ jj = isort(ii);
+ fprintf(fid,['%5i%12s%5i' sTfmti '%10i\n'],ii,labs{jj},jj,nseis_out(jj,:),ntot_2(jj));
+ 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],-[2+nper]);
fid = fopen([odir ftag '_receivers_used_by_event'],'w');
- fprintf(fid,['%12s%5s' sTfmts '%10s\n'],'event',' ',sTbp{:},'TOTAL');
+ fprintf(fid,['%5s%12s%5s' sTfmts '%10s\n'],' ','event','ind',sTbp{:},'TOTAL');
for ii = 1:nevent
jj = isort(ii);
- fprintf(fid,['%12s%5i' sTfmti '%10i\n'],eids{jj},jj,nrec_all_event(jj,:));
+ fprintf(fid,['%5i%12s%5i' sTfmti '%10i\n'],ii,eids{jj},jj,nrec_all_event(jj,:));
end
fclose(fid);
+
+ % list the number of window picks per event, unsorted
+ % (This was used to computed the weighted coverage kernel for m16.)
+ nw = length(nwin_all_event);
+ ntot_2 = sum(nwin_all_event, 2);
+ fid0 = fopen([odir stmod '_window_picks_eids'],'w');
+ fid1 = fopen([odir stmod '_window_picks_nwin'],'w');
+ fid2 = fopen([odir stmod '_window_picks'],'w');
+ for ii = 1:nw
+ jj = isort(ii);
+ fprintf(fid0,['%s\n'],eids{ii});
+ fprintf(fid1,['%i\n'],ntot_2(ii));
+ fprintf(fid2,['%s %i\n'],eids{ii},ntot_2(ii));
+ end
+ fclose(fid0); fclose(fid1); fclose(fid2);
+
end
%----------------------------------------------
@@ -389,15 +422,15 @@
if sum(Ns) ~= length(meas_array), 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 );
- Sval = meas_array(imatch,17);
+ Sval = meas_array(imatch,25);
if sum(isnan(Sval)) > 0
disp(sprintf('i = %i, eid %s',ii,eids{ii}));
error('Sval has at least one NaN data entry');
else
+ % factor of 2 cancels the 0.5 used in computing the misfit (mt_measure_adj.f90)
dnorm_sq(ii) = sum(2 * Sval ./ dcov_fac(imatch) );
end
end
@@ -418,28 +451,31 @@
% This is useful if you want to only list events matching a certain set of criteria.
if 0==1
- eid_reject = load('/net/sierra/raid1/carltape/results/KERNELS/kernels_exclude');
- X1 = 20;
- X2 = 40;
- D1 = 3;
+ eid_reject = load('/net/sierra/raid1/carltape/results/EID_LISTS/eids_reject');
+ X1 = 20; X2 = 40; D1 = 0.4;
isave = find( or(nrec_all_event(:,3) > X2, and(nrec_all_event(:,3) > X1, dnorm_sq*nevent > D1)) );
+ %isave = find( and(nrec_all_event(:,3) > X1, dnorm_sq*nevent > D1));
[eids_save,ijunk] = setdiff( eids_num(isave), eid_reject );
else
- eids_save = eids_num;
+ %eids_save = load('/net/sierra/raid1/carltape/results/EID_LISTS/kernels_run_m14');
+ %eids_save = eids_num(find( or( slat(ievent_full) > 34.6, slat(ievent_full) < 33.8) )); % dm14
+ eids_save = eids_num; % default
end
- klabs = {'dnorm','nrec'};
- kind = -([5 4]+nper); % columns to sort
- for kk = 1:2
+ klabs = {'order','dnorm','nrec'};
+ kind = [2 -(nper+[5 4])]; % columns to sort, corresponding to klabs
+ for kk = 1:length(klabs)
dsort = sortrows(data_norms,kind(kk));
fid = fopen([odir ftag '_data_norms_sort_by_' klabs{kk} '.txt'],'w');
- fprintf(fid,['%12s%4s%8s' sTfmts '%8s%10s%10s%10s%10s\n'],'eid','ind','Nwin',...
+ fprintf(fid,['%4s%10s%4s%8s' sTfmts '%6s%10s%10s%10s%10s\n'],' ','eid','ind','Nwin',...
sTbp{:},'Nrec','dnorm2*E','dnorm2','dnorm','weight');
- fprintf(fid,['%12s%4i%8i' sTfmts '%8s%10.2f%10.4f%10s%10s\n'],...
+ fprintf(fid,['%14s%4i%8i' sTfmts '%6s%10.2f%10.4f%10s%10s\n'],...
'TOTAL',' ',sum(Ns),sTdash{:},'--',sum(dnorm_sq*nevent),sum(dnorm_sq),'--','--');
+ pp = 0;
for ii = 1:nevent
if any(dsort(ii,1)==eids_save)
- fprintf(fid,['%12i%4i%8i' sTfmti '%8i%10.4f%10.4f%10.4f%10.4f\n'],dsort(ii,:));
+ pp = pp+1;
+ fprintf(fid,['%4i%10i%4i%8i' sTfmti '%6i%10.4f%10.4f%10.4f%10.4f\n'],pp,dsort(ii,:));
end
end
%fprintf(fid,'%12s%8i%6s%6s%6s%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq),'--','--');
@@ -475,7 +511,7 @@
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));
+ dep(inds),M(inds,:),seid_cmt(inds));
end
break
@@ -495,28 +531,61 @@
%-------------------------------
% find records that meet particular criteria
-DT_MIN = 5;
+clc
+DT_MIN = 7;
+DT_MAX = 9;
%DT_SIGMA_MIN = 1;
-%icheck = find( and( meas_array(:,9) >= DT_MIN, meas_array(:,10) <= DT_SIGMA_MIN) )
-icheck = find( abs(meas_array(:,9)) >= DT_MIN )
-[junk, isort] = sortrows( meas_array(icheck,:), -9 );
-meas_disp = meas_array(icheck(isort),:);
+icheck = find( and( meas_array(:,19) >= DT_MIN, meas_array(:,19) <= DT_MAX) );
+if 0==1
+ [junk, isort] = sortrows( meas_array(icheck,:), -19 );
+ meas_disp = meas_array(icheck(isort),:);
+else
+ meas_disp = meas_array(icheck,:);
+end
display_meas(meas_disp,Tminvec,eids,strec,comps);
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 );
+icheck = find( or( meas_array(:,21) <= DA_MIN, meas_array(:,21) >= DA_MAX) );
[junk, isort] = sortrows( meas_array(icheck,:), 3 );
meas_disp = meas_array(icheck(isort),:);
display_meas(meas_disp,Tminvec,eids,strec,comps);
CHI_MIN = 100;
DT_SIGMA_MIN = 0.19;
-icheck = find( and( meas_array(:,17) >= CHI_MIN, meas_array(:,10) < DT_SIGMA_MIN) );
-[junk, isort] = sortrows( meas_array(icheck,:), -15 );
+icheck = find( and( meas_array(:,25) >= CHI_MIN, meas_array(:,20) < DT_SIGMA_MIN) );
+[junk, isort] = sortrows( meas_array(icheck,:), -16 );
meas_disp = meas_array(icheck(isort),:);
display_meas(meas_disp,Tminvec,eids,strec,comps);
ebads = unique( eids_num(meas_array(icheck,3)) )
+% records whos max freq content of window exceed the window duration
+% windur = meas_array(:,10);
+% icheck = find( or( meas_array(:,11) >= windur, meas_array(:,12) >= windur ) );
+% [junk, isort] = sortrows( [meas_array(icheck,:) meas_array(icheck,12)-windur(icheck)], -26 );
+% meas_disp = meas_array(icheck(isort),:);
+% display_meas(meas_disp(1:100,:),Tminvec,eids,strec,comps);
+
+% records whos max freq content of window exceed the window duration
+NCYCLE = 1;
+TPMIN = 20; TPMAX = 30;
+iTper = 1;
+windur = meas_array(:,10);
+icheck = find( and( win_rels2 > 0.1, ...
+ and( and(windur >= NCYCLE*TPMIN, meas_array(:,2) == iTper), ...
+ and( and(meas_array(:,11) >= TPMIN, meas_array(:,11) <= TPMAX), ...
+ and(meas_array(:,12) >= TPMIN, meas_array(:,12) <= TPMAX)) )));
+[junk, isort] = sortrows( [meas_array(icheck,:)], 3 );
+meas_disp = meas_array(icheck(isort),:);
+display_meas(meas_disp,Tminvec,eids,strec,comps);
+
+% long windows with prominant phases for 2-30s
+WINMIN = 30;
+iTper = 1;
+windur = meas_array(:,10);
+icheck = find( and( win_rels2 > 0.1, and(windur >= 30, meas_array(:,2) == iTper)));
+[junk, isort] = sortrows( [meas_array(icheck,:)], 3 );
+meas_disp = meas_array(icheck(isort),:);
+display_meas(meas_disp,Tminvec,eids,strec,comps);
+
%=======================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m 2009-03-27 20:51:08 UTC (rev 14496)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m 2009-03-27 21:00:59 UTC (rev 14497)
@@ -13,26 +13,8 @@
X = meas_array;
[m,n] = size(X);
-if n ~= 17, error(' should be 17 columns in the input matrix'); end
+if n ~= 25, error(' should be 25 columns in the input matrix'); end
-% 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
-
index_win = X(:,1);
index_per = X(:,2);
index_event = X(:,3);
@@ -41,20 +23,31 @@
index_comp = X(:,6);
isub_win = X(:,7);
iker = X(:,8);
-measCC_dT = X(:,9);
-sigmaCC_dT = X(:,10);
-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);
+%---------------------------
+seisdur = X(:,9);
+windur = X(:,10);
+T_pmax_dat = X(:,11);
+T_pmax_syn = X(:,12);
+%---------------------------
+seis_d2 = X(:,13);
+win_d2 = X(:,14);
+seis_s2 = X(:,15);
+win_s2 = X(:,16);
+seis_diff2 = X(:,17);
+win_diff2 = X(:,18);
+%---------------------------
+measCC_dT = X(:,19);
+sigmaCC_dT = X(:,20);
+measCC_dA = X(:,21);
+sigmaCC_dA = X(:,22);
+measMT_dT = X(:,23);
+sigmaMT_dT = X(:,24);
+tr_chi = X(:,25);
disp('----------------');
disp('INDEX (window, period, event, receiver, component)');
for ii = 1:m
- disp(sprintf('%10s --> %6.1f %7s %5s %3i DT = %6.2f +- %6.2f DA = %6.2f chi = %6.2f -- %i',...
+ disp(sprintf('%10s --> %6.1f %7s %6s %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) ));
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m 2009-03-27 20:51:08 UTC (rev 14496)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m 2009-03-27 21:00:59 UTC (rev 14497)
@@ -13,18 +13,18 @@
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)
+ wind2,wins2,windiff2,windur,...
+ recd2,recs2,recdiff2,recdur,...
+ tr_chi,am_chi,T_pmax_dat,T_pmax_syn] = 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');
+ wind2,wins2,windiff2,windur,...
+ recd2,recs2,recdiff2,recdur,...
+ tr_chi,am_chi,T_pmax_dat,T_pmax_syn] = ...
+ 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%f%f');
for kk = 1:length(sta), strec{kk} = [sta{kk} '.' netwk{kk}]; end
strec = strec(:);
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m 2009-03-27 21:00:59 UTC (rev 14497)
@@ -0,0 +1,128 @@
+%
+% read_window_chi_all.m
+% CARL TAPE, 14-March-2008
+% printed xxx
+%
+%
+%
+% calls read_window_chi.m
+% called by compute_misfit.m, compare_misfit.m
+%
+
+function meas_array = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps)
+
+nper = length(Tminvec);
+nevent = length(eids);
+nnet = length(stnet);
+nrec = length(strec);
+ncomp = length(comps);
+
+if imax > nevent, error('imax is greater than nevent'); end
+
+k2 = 0;
+meas_array = [];
+%for tt = 1:1
+for tt = 1:nper
+ Tmin = Tminvec(tt);
+ Tmax = Tmaxvec(tt);
+ Ttag = sprintf('T%3.3i_T%3.3i',Tmin,Tmax);
+ Ttag2 = [Ttag '_' stmod];
+ 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_' Ttag '/'];
+ mfile = [eids{ii} '_' Ttag2 '_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,windiff2,windur, ...
+ seisd2,seiss2,seisdiff2,seisdur,...
+ tr_chi,am_chi,T_pmax_dat,T_pmax_syn] = read_window_chi(filename);
+
+ % waveform differences normalized by data-squared
+ %win_chi = windiff2 ./ wind2; % window only
+ %seis_chi = seisdiff2 ./ seisd2; % full record (repeated for seismograms with more than one window)
+
+ % a waveform measure of the power in a window relative to the full record
+ %win_reld2 = wind2 ./ seisd2;
+ %win_rels2 = wins2 ./ seiss2;
+
+ 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]';
+
+ % NOTE: remove the factors of 0.5 that is used in the wdiff definition
+ meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin iker ...
+ seisdur windur T_pmax_dat T_pmax_syn ...
+ 2*seisd2 2*wind2 2*seiss2 2*wins2 2*seisdiff2 2*windiff2 ...
+ measCC_dT sigmaCC_dT measCC_dA sigmaCC_dA ...
+ measMT_dT sigmaMT_dT tr_chi];
+ % 1 kinds
+ % 2 index_T
+ % 3 index_event
+ % 4 index_network
+ % 5 index_rec
+ % 6 index_comp
+ % 7 iwin
+ % 8 iker
+ % 9 seisdur
+ % 10 windur
+ % 11 T_pmax_dat
+ % 12 T_pmax_syn
+ %-------------------
+ % waveform difference information
+ % 13 seisd2
+ % 14 wind2
+ % 15 seiss2
+ % 16 wins2
+ % 17 seisdiff2
+ % 18 windiff2
+ %-------------------
+ % other misfit measures
+ % 19 measCC_dT
+ % 20 sigmaCC_dT
+ % 21 measCC_dA
+ % 22 sigmaCC_dA
+ % 23 measMT_dT
+ % 24 sigmaMT_dT
+ % 25 tr_chi
+ end
+
+ end
+end
+
+%======================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2009-03-27 20:51:08 UTC (rev 14496)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2009-03-27 21:00:59 UTC (rev 14497)
@@ -1,6 +1,6 @@
%
% subspace_specfem.m
-% CARL TAPE, 16-Oct-2008
+% CARL TAPE, 01-Feb-2009
% printed xxx
%
% This program assumes you have FIRST run compute_misfit.m, and also have
More information about the CIG-COMMITS
mailing list