[cig-commits] r15486 - in seismo/3D/ADJOINT_TOMO/iterate_adj/matlab: . matlab_scripts
carltape at geodynamics.org
carltape at geodynamics.org
Thu Jul 30 07:23:34 PDT 2009
Author: carltape
Date: 2009-07-30 07:23:31 -0700 (Thu, 30 Jul 2009)
New Revision: 15486
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/CMT2m0.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/axes_expand.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/curvature.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/defval.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/fontsize.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/gcvfctn.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/hfile2hess.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/m02mw.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/nounder.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ocv_carl.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/osdep.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/plot_histo.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readCMT.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_station_SPECFEM.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi_all.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readsac.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ridge_carl.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/suf.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/wysiwyg.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/year.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/axes_expand.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/curvature.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/defval.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/fontsize.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/gcvfctn.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/hfile2hess.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/nounder.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ocv_carl.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/osdep.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/plot_histo.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readCMT.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_station_SPECFEM.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readsac.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ridge_carl.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/suf.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/wysiwyg.m
Log:
Added matlab_scripts directory to better organize the scripts. Now the main Matlab scripts will load the path to matlab_scripts.
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/axes_expand.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/axes_expand.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/axes_expand.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,50 +0,0 @@
-%
-% function ax1 = axes_expand(ax0,fac)
-% CARL TAPE, 16-Aug-2005
-% printed xxx
-%
-% This function inputs an axes box and outputs a new axes box,
-% expanded by the factor 'fac' in all 2,4,6 dimensions.
-%
-% EXAMPLE:
-% ax0 = [-121 -114]; ax1 = axes_expand(ax0,2)
-% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,2)
-% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,0.30)
-%
-% calls xxx
-% called by xxx
-%
-
-function ax1 = axes_expand(ax0,fac)
-
-% 1D, 2D, 3D
-ndim = length(ax0)/2;
-ax1 = zeros(1,ndim*2);
-
-% return original axes if new axes are non-sensical
-
-xmin = ax0(1);
-xmax = ax0(2);
-dx = xmax-xmin;
-ax1(1) = xmin - dx*(fac-1);
-ax1(2) = xmax + dx*(fac-1);
-if ax1(2) <= ax1(1), ax1 = ax0; end
-
-if ndim >= 2
- ymin = ax0(3);
- ymax = ax0(4);
- dy = ymax-ymin;
- ax1(3) = ymin - dy*(fac-1);
- ax1(4) = ymax + dy*(fac-1);
- if ax1(4) <= ax1(3), ax1 = ax0; end
-end
-if ndim == 3
- zmin = ax0(3);
- zmax = ax0(4);
- dz = zmax-zmin;
- ax1(5) = zmin - dz*(fac-1);
- ax1(6) = zmax + dz*(fac-1);
- if ax1(6) <= ax1(5), ax1 = ax0; end
-end
-
-%===================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,576 +0,0 @@
-%
-% 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
-
-%=======================================================================
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-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,583 @@
+%
+% compare_misfit.m
+% CARL TAPE, 30-July-2009
+% printed xxx
+%
+% This program reads in two sets of window_chi files with identical
+% windows, computes the misfit for both, and then compute the variance
+% reduction within each window.
+%
+% calls read_window_chi.m, plot_histo.m, readCMT.m, read_station_SPECFEM.m
+% called by xxx
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+% add path to additional matlab scripts
+path(path,'./matlab_scripts')
+
+dir0 = '/home/carltape/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 = 0;
+isciencehist = 0; % 0 or 1
+odirsci = '/home/carltape/manuscripts/2009/tomo_gji/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 = ['/home/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 = ['/home/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 IT 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 wdiff (with NaN)
+% 36 --> VRwin wdiff (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 median
+ % 7 - VRseis mean
+ % 8 - VRseis median
+ 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
+
+%=======================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,591 +0,0 @@
-%
-% compute_misfit.m
-% CARL TAPE, 06-Aug-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/
-%
-% 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
-%
-
-%-------------------------
-
-clear
-close all
-format compact
-
-dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
-
-% 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
- % defaults so that the first two letters are always BH
-ncomp = length(comps);
-nper = length(Tminvec);
-
-% 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
-
-imod = input(' Enter the current model number (0, 1, ...): ');
-stmod = sprintf('m%2.2i',imod);
-
-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_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
-
-% 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');
-
-% 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(Tminvec(1))],['Tmin=' num2str(Tminvec(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),seid_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)
-
-[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) : ');
-
-if iread == 1
-
- % range of events
- imin = 1; imax = nevent; % default
- %imin = 139; imax = imin;
-
- meas_array = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps);
-
- save('meas_array.mat','meas_array');
-
-else
- load('meas_array.mat');
-end
-
-whos meas_array
-
-% 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: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,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 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);
-for ii = 1:nevent
- % matches, considering each bandpass set separately
- for tt = 1:nper
- imatch = find( and( meas_array(:,3)==ii, meas_array(:,2)==tt) );
- nrec_all_event(ii,tt) = length(unique(meas_array(imatch,5)));
- end
- % matches, considering all bandpass sets together
- imatch = find(meas_array(:,3)==ii);
- nrec_all_event(ii,nper+1) = length(unique(meas_array(imatch,5)));
-end
-
-if iwrite == 1
- 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; 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]);
-
- fid = fopen([odir ftag '_window_picks_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,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,['%5s%12s%5s' sTfmts '%10s\n'],' ','event','ind',sTbp{:},'TOTAL');
- for ii = 1:nevent
- jj = isort(ii);
- 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
-
-%----------------------------------------------
-
-% construct data covariance normalization terms -- the sigma estimates have
-% already been folded into the misfit function value 0.5 DT^2 / sigma^2
-Ns = zeros(nevent,1);
-dcov_fac = zeros(N,1);
-dcov_fac_e = zeros(nevent,1);
-for ii = 1:nevent
- imatch = find( meas_array(:,3)==ii );
- Ns(ii) = length(imatch);
- if idatacov == 1
- dcov_fac_e(ii) = N;
- elseif idatacov == 2
- dcov_fac_e(ii) = Ns(ii) * nevent;
- elseif idatacov == 3
- dcov_fac_e(ii) = 1;
- end
- dcov_fac(imatch) = dcov_fac_e(ii);
-end
-
-% check number of windows
-if sum(Ns) ~= length(meas_array), error('inconsistent values for total windows'); end
-
-% compute data vector for subspace method
-dnorm_sq = zeros(nevent,1);
-for ii = 1:nevent
- imatch = find( meas_array(:,3)==ii );
- 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
-
-% weight vector for the subspace method
-dnorm = sqrt(dnorm_sq);
-ws = 1 ./ dnorm;
-
-% total data misfit -- this is like a PER WINDOW measure of misfit
-dmisfit = sum( dnorm_sq );
-disp(' total data misfit = d^T Cd^-1 d and its square-root:');
-disp(dmisfit);
-disp(sqrt(dmisfit));
-
-if iwrite == 1
- % display sorted from greatest to least norm
- data_norms = [eids_num [1:nevent]' Ns nrec_all_event dnorm_sq*nevent dnorm_sq dnorm ws];
-
- % 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/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 = 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 = {'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,['%4s%10s%4s%8s' sTfmts '%6s%10s%10s%10s%10s\n'],' ','eid','ind','Nwin',...
- sTbp{:},'Nrec','dnorm2*E','dnorm2','dnorm','weight');
- 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)
- 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),'--','--');
- fclose(fid);
- end
-
- fid1 = fopen([odir ftag '_eid_save.txt'],'w');
- fid2 = fopen([odir ftag '_eid_exclude.txt'],'w');
- for ii = 1:nevent
- if any(eids_num(ii)==eids_save)
- fprintf(fid1,'%i\n',eids_num(ii));
- else
- fprintf(fid2,'%i\n',eids_num(ii));
- end
- end
- 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','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,:),seid_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
-
-clc
-DT_MIN = 7;
-DT_MAX = 9;
-%DT_SIGMA_MIN = 1;
-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(:,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(:,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);
-
-%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,632 @@
+%
+% compute_misfit.m
+% CARL TAPE, 06-Aug-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/
+%
+% 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
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+dir0 = '/home/carltape/RUNS/';
+
+% add path to additional matlab scripts
+path(path,'./matlab_scripts')
+
+% 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
+ % defaults so that the first two letters are always BH
+ncomp = length(comps);
+nper = length(Tminvec);
+
+% 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
+
+imod = input(' Enter the current model number (0, 1, ...): ');
+stmod = sprintf('m%2.2i',imod);
+
+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 = ['/home/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 = ['/home/carltape/results/EID_LISTS/syn_run_' stmod];
+%eid_file = ['/home/carltape/results/EID_LISTS/kernels_run_' stmod];
+%eid_file = ['/home/carltape/results/EID_LISTS/kernels_use_' stmod];
+
+eid_file = '/home/carltape/results/EID_LISTS/eids_simulation';
+%eid_file = '/home/carltape/results/EID_LISTS/eids_tomo';
+%eid_file = '/home/carltape/results/EID_LISTS/eids_extra';
+
+%-------------------------
+% 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');
+
+% 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(Tminvec(1))],['Tmin=' num2str(Tminvec(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),seid_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)
+
+[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) : ');
+
+if iread == 1
+ % range of events
+ imin = 1; imax = nevent; % default
+ %imin = 139; imax = imin;
+
+ meas_array = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps);
+ save('meas_array.mat','meas_array');
+else
+ load('meas_array.mat');
+end
+
+whos meas_array
+
+% comparing multitaper measurments with cross-correlation measurements
+if 0==1
+ iMT = find( meas_array(:,25) ~= 0 );
+
+ figure; plot( meas_array(iMT,21),meas_array(iMT,25),'.');
+ xlabel('DA - CC'); ylabel('DA - MT');
+ axis equal; grid on
+
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1);
+ x = meas_array(iMT,21);
+ edges = [-4:0.25:4]; ylims = [0 0.5];
+ Nbin = plot_histo(x,edges);
+ xlabel('DA - CC for each window'); ylim(ylims); grid on;
+
+ subplot(nr,nc,2);
+ x = meas_array(iMT,25);
+ Nbin = plot_histo(x,edges);
+ xlabel('DA - MT for each window'); ylim(ylims); grid on;
+
+ %----------------------------------------
+
+ figure; plot( meas_array(iMT,19),meas_array(iMT,23),'.');
+ xlabel('DT - CC'); ylabel('DT - MT');
+ axis equal; grid on
+
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1);
+ x = meas_array(iMT,19);
+ edges = [-4:0.5:4]; ylims = [0 0.5];
+ Nbin = plot_histo(x,edges);
+ xlabel('DT - CC for each window'); ylim(ylims); grid on;
+
+ subplot(nr,nc,2);
+ x = meas_array(iMT,23);
+ Nbin = plot_histo(x,edges);
+ xlabel('DT - MT for each window'); ylim(ylims); grid on;
+
+ break
+end
+
+% 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: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,27)
+ 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 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);
+for ii = 1:nevent
+ % matches, considering each bandpass set separately
+ for tt = 1:nper
+ imatch = find( and( meas_array(:,3)==ii, meas_array(:,2)==tt) );
+ nrec_all_event(ii,tt) = length(unique(meas_array(imatch,5)));
+ end
+ % matches, considering all bandpass sets together
+ imatch = find(meas_array(:,3)==ii);
+ nrec_all_event(ii,nper+1) = length(unique(meas_array(imatch,5)));
+end
+
+if iwrite == 1
+ 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; 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]);
+
+ fid = fopen([odir ftag '_window_picks_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,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,['%5s%12s%5s' sTfmts '%10s\n'],' ','event','ind',sTbp{:},'TOTAL');
+ for ii = 1:nevent
+ jj = isort(ii);
+ 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
+
+%----------------------------------------------
+
+% construct data covariance normalization terms -- the sigma estimates have
+% already been folded into the misfit function value 0.5 DT^2 / sigma^2
+Ns = zeros(nevent,1);
+dcov_fac = zeros(N,1);
+dcov_fac_e = zeros(nevent,1);
+for ii = 1:nevent
+ imatch = find( meas_array(:,3)==ii );
+ Ns(ii) = length(imatch);
+ if idatacov == 1
+ dcov_fac_e(ii) = N;
+ elseif idatacov == 2
+ dcov_fac_e(ii) = Ns(ii) * nevent;
+ elseif idatacov == 3
+ dcov_fac_e(ii) = 1;
+ end
+ dcov_fac(imatch) = dcov_fac_e(ii);
+end
+
+% check number of windows
+if sum(Ns) ~= length(meas_array), error('inconsistent values for total windows'); end
+
+% compute data vector for subspace method
+dnorm_sq = zeros(nevent,1);
+for ii = 1:nevent
+ imatch = find( meas_array(:,3)==ii );
+ Sval = meas_array(imatch,27);
+ 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
+
+% weight vector for the subspace method
+dnorm = sqrt(dnorm_sq);
+ws = 1 ./ dnorm;
+
+% total data misfit -- this is like a PER WINDOW measure of misfit
+dmisfit = sum( dnorm_sq );
+disp(' total data misfit = d^T Cd^-1 d and its square-root:');
+disp(dmisfit);
+disp(sqrt(dmisfit));
+
+if iwrite == 1
+ % display sorted from greatest to least norm
+ data_norms = [eids_num [1:nevent]' Ns nrec_all_event dnorm_sq*nevent dnorm_sq dnorm ws];
+
+ % This is useful if you want to only list events matching a certain set of criteria.
+ if 0==1
+ eid_reject = load('/home/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 = load('/home/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 = {'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,['%4s%10s%4s%8s' sTfmts '%6s%10s%10s%10s%10s\n'],' ','eid','ind','Nwin',...
+ sTbp{:},'Nrec','dnorm2*E','dnorm2','dnorm','weight');
+ 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)
+ 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),'--','--');
+ fclose(fid);
+ end
+
+ fid1 = fopen([odir ftag '_eid_save.txt'],'w');
+ fid2 = fopen([odir ftag '_eid_exclude.txt'],'w');
+ for ii = 1:nevent
+ if any(eids_num(ii)==eids_save)
+ fprintf(fid1,'%i\n',eids_num(ii));
+ else
+ fprintf(fid2,'%i\n',eids_num(ii));
+ end
+ end
+ fclose(fid2);
+
+ %km04 = load('/home/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','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,:),seid_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
+
+clc
+DT_MIN = 7;
+DT_MAX = 9;
+%DT_SIGMA_MIN = 1;
+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(:,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(:,27) >= 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);
+
+%=======================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/curvature.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/curvature.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/curvature.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,25 +0,0 @@
-%
-% function [i0,kap] = curvature(x,y)
-% CARL TAPE, 25-July-2005
-%
-% This function compute the curvature of a 1-D function
-% and returns the computed curvature, as well as the
-% point of maximum POSITIVE curvature.
-%
-% This was written for quantifying L-curves.
-%
-% calls xxx
-% called by spline_wang_D.m, test_del2.m
-%
-
-function [i0,kap] = curvature(x,y)
-
-f1 = gradient(y,x);
-f2 = gradient(f1,x);
-kap = f2 ./ (1 + f1.^2).^(3/2); % see Mathworld
-
-%[kap0,i0] = max(abs(kap));
-[kap0,i0] = max(kap);
-
-%============================================================
-
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/defval.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/defval.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/defval.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,24 +0,0 @@
-function defval(name,value)
-% DEFVAL(name,value)
-%
-% Checks if a variable 'name' exists in the caller's workspace,
-% if nonexistent or if empty 'name' is set to value 'value' in caller's workspace
-%
-% Won't work for an unassigned structure variable
-%
-% Last modified by fjsimons-at-alum.mit.edu, 09/12/2007
-
-if ~ischar(name),
- error('The first argument of defval has to be a string (variable name)');
-end
-
-si=1;
-if evalin('caller',[ 'exist(''' name ''')']);,
- si=evalin('caller',[ 'isempty(' name ')']);
-end
-if si,
- assignin('caller',name,value);
- na=dbstack;
- %disp(['Default value used in ' na(2).name ': ' name '=' num2str(value(1,:))])
-end
-
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,58 +0,0 @@
-%
-% display_meas.m
-% CARL TAPE, 18-Mar-2008
-% printed xxx
-%
-% /cig/seismo/3D/mt_measure_adj_work/scripts_tomo/matlab/
-%
-% calls xxx
-% called by xxx
-%
-
-function display_meas(meas_array,Tvec,eids,strec,comps)
-
-X = meas_array;
-[m,n] = size(X);
-if n ~= 25, error(' should be 25 columns in the input matrix'); end
-
-index_win = X(:,1);
-index_per = X(:,2);
-index_event = X(:,3);
-index_net = X(:,4);
-index_rec = X(:,5);
-index_comp = X(:,6);
-isub_win = X(:,7);
-iker = X(:,8);
-%---------------------------
-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 %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) ));
-end
-disp('INDEX (window, period, event, receiver, component)');
-disp('====================');
-
-%======================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/fontsize.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/fontsize.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/fontsize.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,24 +0,0 @@
-%
-% function function fontsize(fs,h)
-% Juliette Artru, 15-Sept-2004
-% printed xxxx
-%
-% This changes the font on an entire figure.
-%
-% fs: font size
-% h: handle graphics (optional)
-%
-
-function fontsize(fs,h);
-
-if(nargin<2);h=gcf;end
-hc=get(gcf,'child');
-hall=hc;
-for k=1:length(hc);
- if(strcmp(get(hc(k),'Type'),'axes'))
- hall=[hall; get(hc(k),'XLabel') ; get(hc(k),'YLabel') ; get(hc(k),'Title')];
- end
-end
-set(hall,'fontsize',fs);
-
-
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/gcvfctn.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/gcvfctn.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/gcvfctn.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,47 +0,0 @@
-%
-% function G = gcvfctn(h, s2, fc, rss0, dof0)
-% Carl Tape (Tapio Schneider, ACM 118)
-% 09-Nov-2004
-%
-% GCVFCTN Evaluate generalized cross-validation function.
-%
-% gcvfctn(h, s2, fc, rss0, dof0) is the value of the GCV function
-% for ridge regression with regularization parameter h.
-%
-% INPUT:
-% h regularization parameter
-% s2 squared singular value of the design matrix
-% fc coefficients fc = U(:, 1:q)'*g
-% rss0 the residual sum of squares of the (non-regularized)
-% least squares estimate
-% dofo0 the number of residual degrees of freedom of
-% the (non-regularized) least squares estimate
-%
-% U matrix of left singular vectors
-% g vector of response variables
-% q the smaller of the numbers of rows and columns of
-% the design matrix, q = min(n,p)
-%
-% Auxiliary function for GCV.
-%
-% Adapted from GCVFUN in Per-Christian Hansen's Regularization Toolbox.
-%
-% See Schneider (2001) for details.
-%
-% calls xxx
-% called by gcv.m
-%
-
-function G = gcvfctn(h, s2, fc, rss0, dof0)
-
-% SIMILAR TO filter factor for Tikhonov regularization:
-% 1 - f = F, where F are the filter factors
-f = h^2 ./ (s2 + h^2);
-
-RSS = norm(f.*fc)^2 + rss0;
-T2 = (dof0 + sum(f))^2;
-
-G = RSS / T2;
-
-%=============================================================
-
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/hfile2hess.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/hfile2hess.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/hfile2hess.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,31 +0,0 @@
-%
-% display_meas.m
-% CARL TAPE, 18-April-2008
-% printed xxx
-%
-% /cig/seismo/3D/mt_measure_adj_work/scripts_tomo/matlab/
-%
-% calls xxx
-% called by xxx
-%
-
-function H = hfile2hess(hfile,nsrc)
-
-if ~exist(hfile,'file'), error([hfile ' does not exist']); end
-[kall,iall,jall,Hij] = textread(hfile,'%f%f%f%f');
-
-n = length(kall);
-nhess = (-1 + sqrt(1+8*n))/2; % quadratic formula
-if nhess ~= nsrc, error('number of events must be consistent'); end
-
-% construct the Hessian matrix
-H = zeros(nsrc,nsrc);
-for i = 1:nsrc
- for j = i:nsrc
- ih = find( and( i == iall, j == jall) );
- H(i,j) = Hij(ih);
- end
-end
-H = H + H' - diag(diag(H)); % fill the lower triangle
-
-%======================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/CMT2m0.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/CMT2m0.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/CMT2m0.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,64 @@
+%
+% function M0 = CMT2m0(im0,M)
+% CARL TAPE, 02-Feb-2007
+% printed xxx
+%
+% This function converts from a (CMT) moment tensor to scalar seismic
+% moment. The units of M0 are the same as the elements of Mij, which
+% should be Newton-meter (N-m), although it does not affect the function
+% here.
+%
+% See Ekstrom email (11-Oct-2006) and corresponding Latex notes.
+%
+% moment tensor M = [Mrr Mtt Mpp Mrt Mrp Mtp]
+%
+% calls xxx
+% called by richter.m, denaliquake.m
+%
+
+function M0 = CMT2m0(im0,M)
+
+% make sure M is N by 6
+[a,b] = size(M); if b ~= 6, M = M'; end
+[a,b] = size(M); if b ~= 6, error('dimension of M must be N by 6'); end
+
+ncmt = a;
+Mrr = M(:,1); Mtt = M(:,2); Mpp = M(:,3);
+Mrt = M(:,4); Mrp = M(:,5); Mtp = M(:,6);
+
+M0 = zeros(ncmt,1);
+for ii=1:ncmt
+
+ % if you need to compute the eigenvalues
+ if or(im0==2,im0==3)
+ % convention: r (up), theta (south), phi (east)
+ Mcmt = zeros(3,3);
+ Mcmt = [ Mrr(ii) Mrt(ii) Mrp(ii) ;
+ Mrt(ii) Mtt(ii) Mtp(ii) ;
+ Mrp(ii) Mtp(ii) Mpp(ii) ];
+
+ [V, D] = eig(Mcmt);
+ lams = diag(D)';
+ isign = sign(lams);
+ %(inv(V)*Mcmt*V - D) / norm(Mcmt)
+
+ % adjust the order of eigenvectors
+ % [lamsort, isort] = sort(abs(lams),'descend');
+ % Vsort = V(:,isort)
+ % Dsort = diag(lamsort.*isign(isort))
+ % (inv(Vsort)*Mcmt*Vsort - Dsort) / norm(Mcmt) % check
+ end
+
+ % formula to convert M to M0
+ % (1) Dahlen and Tromp, eq. 5.91
+ % (2) Harvard CMT way -- double-couple moment
+ % (3) Caltech way
+ switch im0
+ case 1, M0(ii) = 1/sqrt(2) * sqrt( Mrr(ii)^2 + Mtt(ii)^2 + Mpp(ii)^2 ...
+ + 2*(Mrt(ii)*Mrt(ii) + Mrp(ii)*Mrp(ii) + Mtp(ii)*Mtp(ii) ) );
+ case 2, M0(ii) = (abs(min(lams)) + abs(max(lams)) ) / 2;
+ case 3, M0(ii) = max(abs(lams));
+ end
+end
+
+%=====================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/CMT2m0.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/axes_expand.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/axes_expand.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/axes_expand.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/axes_expand.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,50 @@
+%
+% function ax1 = axes_expand(ax0,fac)
+% CARL TAPE, 16-Aug-2005
+% printed xxx
+%
+% This function inputs an axes box and outputs a new axes box,
+% expanded by the factor 'fac' in all 2,4,6 dimensions.
+%
+% EXAMPLE:
+% ax0 = [-121 -114]; ax1 = axes_expand(ax0,2)
+% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,2)
+% ax0 = [-121 -114 31 37]; ax1 = axes_expand(ax0,0.30)
+%
+% calls xxx
+% called by xxx
+%
+
+function ax1 = axes_expand(ax0,fac)
+
+% 1D, 2D, 3D
+ndim = length(ax0)/2;
+ax1 = zeros(1,ndim*2);
+
+% return original axes if new axes are non-sensical
+
+xmin = ax0(1);
+xmax = ax0(2);
+dx = xmax-xmin;
+ax1(1) = xmin - dx*(fac-1);
+ax1(2) = xmax + dx*(fac-1);
+if ax1(2) <= ax1(1), ax1 = ax0; end
+
+if ndim >= 2
+ ymin = ax0(3);
+ ymax = ax0(4);
+ dy = ymax-ymin;
+ ax1(3) = ymin - dy*(fac-1);
+ ax1(4) = ymax + dy*(fac-1);
+ if ax1(4) <= ax1(3), ax1 = ax0; end
+end
+if ndim == 3
+ zmin = ax0(3);
+ zmax = ax0(4);
+ dz = zmax-zmin;
+ ax1(5) = zmin - dz*(fac-1);
+ ax1(6) = zmax + dz*(fac-1);
+ if ax1(6) <= ax1(5), ax1 = ax0; end
+end
+
+%===================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/axes_expand.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/curvature.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/curvature.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/curvature.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/curvature.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,25 @@
+%
+% function [i0,kap] = curvature(x,y)
+% CARL TAPE, 25-July-2005
+%
+% This function compute the curvature of a 1-D function
+% and returns the computed curvature, as well as the
+% point of maximum POSITIVE curvature.
+%
+% This was written for quantifying L-curves.
+%
+% calls xxx
+% called by spline_wang_D.m, test_del2.m
+%
+
+function [i0,kap] = curvature(x,y)
+
+f1 = gradient(y,x);
+f2 = gradient(f1,x);
+kap = f2 ./ (1 + f1.^2).^(3/2); % see Mathworld
+
+%[kap0,i0] = max(abs(kap));
+[kap0,i0] = max(kap);
+
+%============================================================
+
\ No newline at end of file
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/curvature.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/defval.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/defval.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/defval.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/defval.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,24 @@
+function defval(name,value)
+% DEFVAL(name,value)
+%
+% Checks if a variable 'name' exists in the caller's workspace,
+% if nonexistent or if empty 'name' is set to value 'value' in caller's workspace
+%
+% Won't work for an unassigned structure variable
+%
+% Last modified by fjsimons-at-alum.mit.edu, 09/12/2007
+
+if ~ischar(name),
+ error('The first argument of defval has to be a string (variable name)');
+end
+
+si=1;
+if evalin('caller',[ 'exist(''' name ''')']);,
+ si=evalin('caller',[ 'isempty(' name ')']);
+end
+if si,
+ assignin('caller',name,value);
+ na=dbstack;
+ %disp(['Default value used in ' na(2).name ': ' name '=' num2str(value(1,:))])
+end
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/defval.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/display_meas.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,58 @@
+%
+% display_meas.m
+% CARL TAPE, 18-Mar-2008
+% printed xxx
+%
+% /cig/seismo/3D/mt_measure_adj_work/scripts_tomo/matlab/
+%
+% calls xxx
+% called by xxx
+%
+
+function display_meas(meas_array,Tvec,eids,strec,comps)
+
+X = meas_array;
+[m,n] = size(X);
+if n ~= 25, error(' should be 25 columns in the input matrix'); end
+
+index_win = X(:,1);
+index_per = X(:,2);
+index_event = X(:,3);
+index_net = X(:,4);
+index_rec = X(:,5);
+index_comp = X(:,6);
+isub_win = X(:,7);
+iker = X(:,8);
+%---------------------------
+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 %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) ));
+end
+disp('INDEX (window, period, event, receiver, component)');
+disp('====================');
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/fontsize.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/fontsize.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/fontsize.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/fontsize.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,24 @@
+%
+% function function fontsize(fs,h)
+% Juliette Artru, 15-Sept-2004
+% printed xxxx
+%
+% This changes the font on an entire figure.
+%
+% fs: font size
+% h: handle graphics (optional)
+%
+
+function fontsize(fs,h);
+
+if(nargin<2);h=gcf;end
+hc=get(gcf,'child');
+hall=hc;
+for k=1:length(hc);
+ if(strcmp(get(hc(k),'Type'),'axes'))
+ hall=[hall; get(hc(k),'XLabel') ; get(hc(k),'YLabel') ; get(hc(k),'Title')];
+ end
+end
+set(hall,'fontsize',fs);
+
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/fontsize.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/gcvfctn.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/gcvfctn.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/gcvfctn.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/gcvfctn.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,47 @@
+%
+% function G = gcvfctn(h, s2, fc, rss0, dof0)
+% Carl Tape (Tapio Schneider, ACM 118)
+% 09-Nov-2004
+%
+% GCVFCTN Evaluate generalized cross-validation function.
+%
+% gcvfctn(h, s2, fc, rss0, dof0) is the value of the GCV function
+% for ridge regression with regularization parameter h.
+%
+% INPUT:
+% h regularization parameter
+% s2 squared singular value of the design matrix
+% fc coefficients fc = U(:, 1:q)'*g
+% rss0 the residual sum of squares of the (non-regularized)
+% least squares estimate
+% dofo0 the number of residual degrees of freedom of
+% the (non-regularized) least squares estimate
+%
+% U matrix of left singular vectors
+% g vector of response variables
+% q the smaller of the numbers of rows and columns of
+% the design matrix, q = min(n,p)
+%
+% Auxiliary function for GCV.
+%
+% Adapted from GCVFUN in Per-Christian Hansen's Regularization Toolbox.
+%
+% See Schneider (2001) for details.
+%
+% calls xxx
+% called by gcv.m
+%
+
+function G = gcvfctn(h, s2, fc, rss0, dof0)
+
+% SIMILAR TO filter factor for Tikhonov regularization:
+% 1 - f = F, where F are the filter factors
+f = h^2 ./ (s2 + h^2);
+
+RSS = norm(f.*fc)^2 + rss0;
+T2 = (dof0 + sum(f))^2;
+
+G = RSS / T2;
+
+%=============================================================
+
\ No newline at end of file
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/gcvfctn.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/hfile2hess.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/hfile2hess.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/hfile2hess.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/hfile2hess.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,31 @@
+%
+% display_meas.m
+% CARL TAPE, 18-April-2008
+% printed xxx
+%
+% /cig/seismo/3D/mt_measure_adj_work/scripts_tomo/matlab/
+%
+% calls xxx
+% called by xxx
+%
+
+function H = hfile2hess(hfile,nsrc)
+
+if ~exist(hfile,'file'), error([hfile ' does not exist']); end
+[kall,iall,jall,Hij] = textread(hfile,'%f%f%f%f');
+
+n = length(kall);
+nhess = (-1 + sqrt(1+8*n))/2; % quadratic formula
+if nhess ~= nsrc, error('number of events must be consistent'); end
+
+% construct the Hessian matrix
+H = zeros(nsrc,nsrc);
+for i = 1:nsrc
+ for j = i:nsrc
+ ih = find( and( i == iall, j == jall) );
+ H(i,j) = Hij(ih);
+ end
+end
+H = H + H' - diag(diag(H)); % fill the lower triangle
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/hfile2hess.m
___________________________________________________________________
Name: svn:mergeinfo
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/m02mw.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/m02mw.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/m02mw.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,50 @@
+%
+% function Mw = m02mw(imag,M0)
+% CARL TAPE, 31-Oct-2007
+% printed xxx
+%
+% This function inputs a vector of moments (M0) in units of N-m and outputs
+% a column vector of moment magnitudes (Mw).
+%
+% See Latex notes cmt.pdf.
+%
+% EXAMPLE (Denali earthquake):
+% M0 = 7.48*1e20; Mw = m02mw(1,M0)
+%
+% calls xxx
+% called by richter.m
+%
+
+function Mw = m02mw(imag,M0)
+
+% convert moment tensor from N-m to dyne-cm, since these formulas are
+% designed for dyne-cm
+M0 = 1e7 * M0(:);
+
+if imag==1
+ % Kanamori 1978
+ %k = -10.7;
+
+ % Kanamori 1977
+ %k = -(2/3)*(11.8 - log10(1/2e4)); % for M0 in units of dyne-cm
+ %k = -(2/3)*(11.8 - log10(5e-5)); % dyne-cm
+ %k = -(2/3)*(16.8 - log10(5)); % dyne-cm
+ %k = -11.2 + (2/3)*log10(5); % dyne-cm
+
+ % Kanamori 1977 or 1978
+ %Mw = (2/3) * log10(M0) + k;
+
+ % Kanamori 1977 (modified form, but exact)
+ A = 2/(3*log(10));
+ K = 0.2*10^16.8;
+ Mw = A*log(M0/K);
+
+elseif imag == 2
+ % Harvard CMT
+ Mw = (2/3) * (log10(M0) - 16.1);
+
+else
+ error('imag must be 1 or 2.');
+end
+
+%=====================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/m02mw.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/nounder.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/nounder.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/nounder.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/nounder.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,9 @@
+function sin=nounder(sin)
+% sin=NOUNDER(sin)
+%
+% Changes underscores to dashes in a string
+%
+% Last modified by fjsimons-at-alum.mit.edu, Feb 06th, 2004
+
+sin(find(abs(sin)==95))='-';
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/nounder.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ocv_carl.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ocv_carl.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ocv_carl.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ocv_carl.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,87 @@
+%
+% function rss_vec = ocv_carl(g, A, lamvec)
+% Carl Tape, 30-March-2006
+%
+% Copied from gcv_carl.m on 27-March-2006.
+%
+% Returns the ordinary cross-validation (OCV) function corresponding to a
+% set of input regularization (or damping) parameters.
+%
+% TWO ALGORITHMS ARE SHOWN:
+% (1) Brute force, which involves ndata inversions per lambda
+% (2) Elegant formalisum, which involves one inversion per lambda
+% --> See Latex notes
+% /home/carltape/classes/acm118/2006_handouts/hw3/hw3sol_2006_prob3.pdf
+%
+% Using some GPS data, I checked that these approaches are identical.
+%
+% calls xxx
+% called by ridge_carl.m
+%
+
+function rss_vec = ocv_carl(d, A, lamvec)
+
+% Size of inputs
+[ndata, nparm] = size(A);
+numlam = length(lamvec);
+
+if (min(lamvec) < 0)
+ error('Impossible regularization parameter lambda.')
+end
+
+rss_vec = zeros(numlam,1);
+
+% loop over regularization parameters
+for ii=1:numlam
+ lam = lamvec(ii);
+ disp([' ii = ' num2str(ii) '/' num2str(numlam) ', lam = ' num2str(lam)]);
+
+ if 1==1
+ H = A*inv(A'*A + lam^2*eye(nparm))*A';
+ dhat = H*d;
+
+ % OCV residual
+ res = (d - dhat) ./ (1 - diag(H));
+
+ % sum the residuals
+ rss_vec(ii) = sum(res.^2);
+
+ else
+ % loop over datapoints
+ for jj=1:ndata
+ %disp([' jj = ' num2str(jj) '/' num2str(ndata) ]);
+
+ % indices for which you compute the model parameters
+ switch jj
+ case 1, einds = [2:ndata];
+ case ndata, einds = [1:ndata-1];
+ otherwise, einds = [1:jj-1 jj+1:ndata];
+ end
+
+ % indices to estimate the RSS
+ oinds = jj;
+
+ % reduced matrices
+ X = A(einds,:);
+ g = d(einds);
+
+ % note: regularization matrix is identity matrix
+ f_h = inv(X'*X + lam^2*eye(nparm))*X'*g;
+
+ % estimate the model at the datapoints NOT used
+ res = A(oinds,:) * f_h - d(oinds);
+
+ % sum the residuals
+ rss_vec(ii) = rss_vec(ii) + sum(res.^2);
+ end
+ end
+
+end
+
+rss_vec = rss_vec / ndata; % normalization
+
+figure; loglog(lamvec, rss_vec, '.'); grid on;
+xlabel(' Regularization parameter, log (\lambda)');
+ylabel(' RSS at datapoints NOT used in computing model');
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ocv_carl.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/osdep.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/osdep.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/osdep.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/osdep.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,16 @@
+function of=osdep
+% of=OSDEP
+%
+% Returns the value of the read option to be used
+% on the LOCAL operating system for files created
+% on the LOCAL operating system.
+%
+% Last modified by fjsimons-at-alum.mit.edu, 23.11.2004
+
+
+if strcmp(getenv('OSTYPE'),'linux')
+ of= 'l';
+end
+if strcmp(getenv('OSTYPE'),'solaris')
+ of= 'b';
+end
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/osdep.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/plot_histo.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/plot_histo.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/plot_histo.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/plot_histo.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,22 @@
+%
+% function N = plot_histo(hdat,edges)
+% CARL TAPE, 30-May-2007
+% printed xxx
+%
+% This function plots a histogram with cyan bars and black boundaries.
+%
+% calls xxx
+% called by xxx
+%
+
+function N = plot_histo(hdat,edges)
+
+Ntotal = length(hdat);
+[N,bin] = histc(hdat,edges);
+bar(edges,N/Ntotal,'histc');
+xlim([min(edges) max(edges)]);
+ylabel([' Fraction (N = ' num2str(Ntotal) ')']);
+
+h = findobj(gca,'Type','patch'); set(h,'FaceColor',[0 1 1],'EdgeColor','k');
+
+%=======================================================================
\ No newline at end of file
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/plot_histo.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readCMT.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readCMT.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readCMT.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readCMT.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,85 @@
+%
+% readCMT.m
+% CARL TAPE, 26-June-2007
+% printed xxx
+%
+% This file reads a concatenated CMTSOLUTION file and outputs the data.
+% The CMT tshift parameter is added to the origin time, which is the
+% time used for the initial CMT solution.
+%
+% INPUT:
+% filename file containing all CMTSOLUTION files concatenated together
+% nlines_cmt number of lines per CMT solution (13)
+% nspace_between_entries number of spaces between CMTSOLUTION blocks (0 or 1)
+%
+% OUTPUT:
+% M moment tensor in N-m units
+%
+% calls xxx
+% called by slab_readCMT.m, socal_quakes_2007.m
+%
+
+function [date,tshift,hdur,lat,lon,dep,M,eid,elabel] = ...
+ readCMT(filename, nlines_cmt, nspace_between_entries)
+
+% read in concatenated CMTSOLUTION files
+lines = textread(filename,'%s','delimiter','\n');
+nlines = length(lines);
+
+% number of events
+enum = (nlines + nspace_between_entries) / (nlines_cmt + nspace_between_entries);
+if mod(enum,1) ~= 0, error(' enum should be an integer'); end
+disp([' File : ' filename ]);
+disp([' Number of events : ' num2str(enum) ]);
+
+% initialize vectors
+date = zeros(enum,1); tshift = zeros(enum,1); hdur = zeros(enum,1);
+lat = zeros(enum,1); lon = zeros(enum,1); dep = zeros(enum,1);
+Mrr = zeros(enum,1); Mtt = zeros(enum,1); Mpp = zeros(enum,1);
+Mrt = zeros(enum,1); Mrp = zeros(enum,1); Mtp = zeros(enum,1);
+
+for kk = 1:enum
+ in1 = (kk-1)*(nlines_cmt + nspace_between_entries);
+
+ % first full line of CMTSOLUTION, as a string
+ ltemp = lines{in1+1};
+
+ if length(ltemp) == 0, error('the CMTSOLUTION line is length zero'); end
+
+ % replace W with space for more recent CMT solutions
+ % SWEQ2006
+ % PDEW2006
+ if or(strcmp(ltemp(4),'W'),strcmp(ltemp(4),'Q')), ltemp(4) = ' '; end
+
+ [j1,yr,mo,dy,hr,min,sec,lat_pde(kk),lon_pde(kk),dep_pde(kk),j5,j6,name1,name2,name3,name4,name5,name6] = ...
+ strread(ltemp,'%s%f%f%f%f%f%f%f%f%f%f%f%s%s%s%s%s%s');
+ [j1,j2,eid(kk)] = strread(lines{in1+2},'%s%s%s');
+ [j1,j2,tshift(kk)] = strread(lines{in1+3},'%s%s%f');
+ [j1,j2,hdur(kk)] = strread(lines{in1+4},'%s%s%f');
+ [j1,lat(kk)] = strread(lines{in1+5},'%s%f');
+ [j1,lon(kk)] = strread(lines{in1+6},'%s%f');
+ [j1,dep(kk)] = strread(lines{in1+7},'%s%f');
+ [j1,Mrr(kk)] = strread(lines{in1+8},'%s%f');
+ [j1,Mtt(kk)] = strread(lines{in1+9},'%s%f');
+ [j1,Mpp(kk)] = strread(lines{in1+10},'%s%f');
+ [j1,Mrt(kk)] = strread(lines{in1+11},'%s%f');
+ [j1,Mrp(kk)] = strread(lines{in1+12},'%s%f');
+ [j1,Mtp(kk)] = strread(lines{in1+13},'%s%f');
+
+ % NOTE the inclusion of the time shift here for the CMT solution
+ date(kk) = datenum(yr,mo,dy,hr,min, sec + tshift(kk));
+
+ elabel{kk} = [char(name1) ' ' char(name2) ' ' char(name3) ...
+ ' ' char(name4) ' ' char(name5) ' ' char(name6)];
+end
+
+eid = eid(:);
+elabel = elabel(:);
+
+% moment tensor elements
+M = [Mrr Mtt Mpp Mrt Mrp Mtp];
+
+% convert moment tensor from dyne-cm to N-m
+M = 1e-7 * M;
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readCMT.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_station_SPECFEM.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_station_SPECFEM.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_station_SPECFEM.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_station_SPECFEM.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,23 @@
+%
+% function [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(filename)
+% CARL TAPE, 15-March-2007
+% printed xxx
+%
+% This function inputs a STATIONS file formatted for SPECFEM3D.f90.
+% See also write_station_SPECFEM.m
+%
+% In Perl, this can be done as follows:
+% open(IN,${stations_seis}); @temp = <IN>; $nrec = $temp[0]; chomp($nrec);
+% print CSH "tail -n $nrec ${stations_seis} > stemp0 \n";
+% print CSH "awk '{print \$4,\$3}' stemp0 > temp \n";
+%
+% calls xxx
+% called by xxx
+%
+
+function [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(filename)
+
+% read in the STATIONS file
+[stnm,netwk,rlat,rlon,relev,rburial] = textread(filename,'%s%s%f%f%f%f','headerlines',1);
+
+%=======================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_station_SPECFEM.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -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,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,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(:);
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi_all.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi_all.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi_all.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,131 @@
+%
+% 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 measMT_dA sigmaMT_dA tr_chi am_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 measMT_dA
+ % 26 sigmaMT_dA
+ % 27 tr_chi
+ % 28 am_chi
+ end
+
+ end
+end
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/read_window_chi_all.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readsac.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readsac.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readsac.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readsac.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,110 @@
+function varargout=readsac(filename,plotornot,osd)
+% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename)
+% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename,1)
+% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename,1,'l')
+% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename,1,'b')
+%
+% READSAC(...)
+%
+% Reads SAC-formatted data (and makes a plot)
+%
+% INPUT:
+%
+% filename The filename, full path included
+% plotornot 1 Makes a plot
+% 0 Does not make a plot [default]
+% osd 'b' for data saved on Solaris read into Linux
+% 'l' for data saved on Linux read into Linux
+%
+% OUPUT:
+%
+% SeisData The number vector
+% HdrData The header structure array
+% tnu The handle to the plot title; or appropriate string
+% pobj The handle to the plot line and the xlabel
+% tims The times on the x-axis
+%
+% See SETVAR.PL out of SACLIB.PM
+%
+% Last modified by fjsimons-at-alum.mit.edu, 07/07/2008
+
+defval('plotornot',0)
+defval('osd',osd)
+
+fid=fopen(filename,'r',osd);
+if fid==-1
+ error([ 'File ',filename,' does not exist in current path ',pwd])
+end
+HdrF=fread(fid,70,'float32');
+HdrN=fread(fid,15,'int32');
+HdrI=fread(fid,20,'int32');
+HdrL=fread(fid,5,'int32');
+HdrK=str2mat(fread(fid,[8 24],'char'))';
+SeisData=fread(fid,HdrN(10),'float32');
+fclose(fid);
+
+% Meaning of LCALDA from http://www.llnl.gov/sac/SAC_Manuals/FileFormatPt2.html
+% LCALDA is TRUE if DIST, AZ, BAZ, and GCARC are to be calculated from
+% station and event coordinates. But I don't know how to set this
+% variable to a logical yet using setvar.pl. What is the TRUE/FALSE
+% format in SAC? 18.1.2006
+
+% If you change any of this, change it in WRITESAC as well!
+HdrData=struct(...
+ 'AZ',HdrF(52),...
+ 'B',HdrF(6),...
+ 'BAZ',HdrF(53),...
+ 'DIST',HdrF(51),...
+ 'DELTA',HdrF(1),...
+ 'E',HdrF(7),...
+ 'EVDP',HdrF(39),...
+ 'EVEL',HdrF(38),...
+ 'EVLA',HdrF(36),...
+ 'EVLO',HdrF(37),...
+ 'GCARC',HdrF(54),...
+ 'IFTYPE',HdrI(1),...
+ 'KCMPNM',HdrK(21,HdrK(21,:)>64),...
+ 'KINST',HdrK(24,HdrK(24,:)>64),...
+ 'KSTNM',HdrK(1,HdrK(1,:)>64),...
+ 'KUSER0',HdrK(18,HdrK(18,:)>64),...
+ 'LCALDA',HdrL(4),...
+ 'MAG',HdrF(40),...
+ 'NPTS',HdrN(10),...
+ 'NVHDR',HdrN(7),...
+ 'NZHOUR',HdrN(3),...
+ 'NZJDAY',HdrN(2),...
+ 'NZMIN',HdrN(4),...
+ 'NZMSEC',HdrN(6),...
+ 'NZSEC',HdrN(5),...
+ 'NZYEAR',HdrN(1),...
+ 'SCALE',HdrF(4),...
+ 'STDP',HdrF(35),...
+ 'STEL',HdrF(34),...
+ 'STLA',HdrF(32),...
+ 'STLO',HdrF(33),...
+ 'T0',HdrF(11),...
+ 'T1',HdrF(12),...
+ 'T2',HdrF(13),...
+ 'T3',HdrF(14));
+
+tims=linspace(HdrData.B,HdrData.E,HdrData.NPTS);
+
+if plotornot==1
+ pobj=plot(tims,SeisData);
+ filename(find(abs(filename)==95))='-';
+ tito=nounder(suf(filename,'/'));
+ if isempty(tito)
+ tito=nounder(filename);
+ end
+ tnu=title(tito);
+ axis tight
+ pobj(2)=xlabel([ 'Time (s)']);
+else
+ tnu=nounder(suf(filename,'/'));
+ pobj=0;
+end
+
+nam={'SeisData' 'HdrData' 'tnu' 'pobj','tims'};
+for index=1:nargout
+ varargout{index}=eval(nam{index});
+end
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/readsac.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ridge_carl.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ridge_carl.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ridge_carl.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ridge_carl.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,227 @@
+%
+% function
+% Carl Tape (Tapio Schneider, ACM 118)
+% 06-Nov-2006
+%
+% This function inputs a design matrix, a data vector, and a vector of
+% regularization parameters, and it returns three different curves that may
+% be used to select the best parameter:
+% (1) L-curve and curvature
+% (2) generalized cross-validation curve (GCV)
+% (3) ordinary cross-validation (OCV), also known as 'leave-one-out' CV
+%
+% It is best to input a large number of regularization parameters, so that
+% the min and max of the respective functions can be easily obtained.
+%
+% This program is copied in part from ridge_tapio.m
+%
+% NOTE THE PLOTTING OPTIONS.
+%
+%-------------------------------------------------
+% RIDGE Ridge regression estimates.
+%
+% Given a vector g, a design matrix X, and
+% a regularization parameter h,
+%
+% [m, rss, mss, dof] = ridge_tapio(g, X, h)
+%
+% returns the ridge regression estimate of the vector f in the
+% linear regression model
+%
+% g = X*f + noise.
+%
+% Also returned are the residual sum of squares rss, the sum of
+% squares mss of the elements of the ridge regression estimate
+% m (the squared norm of m), and the effective number of
+% residual degrees of freedom dof.
+%
+% If h is a vector of regularization parameters, the i-th column
+% m(:,i) is the ridge regression estimate for the regularization
+% parameter h(i); the i-th elements of rss and mss are the
+% associated residual sum of squares and estimate sum of squares.
+%
+% If no regularization parameter h is given, generalized
+% cross-validation is used to determine the regularization
+% parameter. The chosen regularization parameter h and the value of
+% the GCV function are then returned as the fifth and sixth
+% output arguments
+%
+% [m, rss, mss, dof, h, G] = ridge_tapio(g, X);
+%
+% Adapted from various routines in Per Christian Hansen's
+% Regularization Toolbox.
+%
+% calls gcvfctn.m, curvature.m
+% called by xxx
+%
+
+function [m, rss, mss, Gvec, Fvec, dof, kap, iL, iGCV, iOCV] = ridge_carl(dvec, X, hvec)
+
+% Size of inputs
+[n, p] = size(X);
+q = min(n, p);
+nh = length(hvec);
+if (min(hvec) < 0)
+ error('Impossible regularization parameter h.')
+end
+
+% Initialize outputs
+m = zeros(p, nh);
+rss = zeros(nh, 1);
+mss = zeros(nh, 1);
+dof = zeros(nh, 1);
+
+% Compute SVD of X
+[U, S, V] = svd(X, 0);
+s = diag(S); % vector of singular values
+s2 = s.^2;
+
+% Coefficients in expansion of solution in terms of right singular vectors
+fc = U(:, 1:q)'*dvec;
+zeta = s .* fc;
+
+% Treat each regularization parameter separately
+for j = 1:nh
+ m(:, j) = V(:, 1:q) * (zeta ./ (s2 + hvec(j)^2));
+ mss(j) = sum(m(:, j).^2);
+ rss(j) = hvec(j)^4 * sum(fc.^2 ./ (s2 + hvec(j)^2).^2);
+ dof(j) = n - sum(s2./(s2 + hvec(j)^2));
+end
+
+% In overdetermined case, add rss of least-squares problem
+if (n > p)
+ rss = rss + sum((dvec - U(:, 1:q)*fc).^2);
+end
+
+%-----------------------
+% determine the Lcurve pick (max curvature)
+
+x1 = log10(rss);
+y1 = log10(mss);
+
+% % smooth curvature interpolation to get h_L
+% num = 1000;
+% xsmooth = linspace(x1(1),x1(end),1000);
+% ysmooth = interp1(x1,y1,xsmooth,'cubic');
+% [i0,kap_smooth] = curvature(xsmooth,ysmooth);
+% rss_L = 10^xsmooth(i0);
+% mss_L = 10^ysmooth(i0);
+% h_L = 10^interp1(x1,log10(hvec),xsmooth(i0),'cubic');
+
+% curvature, based on input h values alone
+[iL,kap] = curvature(x1,y1);
+%h_L = hvec(iL);
+%rss_L = 10^x1(iL);
+%mss_L = 10^y1(iL);
+
+%-----------------------
+% obtain GCV `best' solution and GCV curve
+
+% GCV minimum -- 'exact' in the sense a minimization method is used
+% [hmin, Gmin] = gcv_tapio(U, s, dvec, 'ridge');
+
+% GCV minimum -- 'crude' in the sense that we coarsely sample the function
+dof0 = n-q;
+rss0 = sum((dvec - U(:, 1:q)*fc).^2);
+Gvec = zeros(nh,1);
+for j = 1:nh
+ Gvec(j) = gcvfctn(hvec(j), s2, fc, rss0, dof0);
+end
+[Gmin,iGCV] = min(Gvec);
+%hmin = hvec(iGCV);
+
+% GCV best model and L-curve point for h_GCV (= hmin)
+%mod_min = inv(X'*X + hmin^2*eye(p))*X'*dvec;
+%res = X*mod_min - dvec;
+%rss_min = sum(res.^2);
+%mss_min = sum(mod_min.^2);
+
+% compute G for the Lcurve pick
+%G_L = gcvfctn(h_L, s2, fc, rss0, dof0);
+
+%-----------------------
+% ordinary (leave-one-out) cross-validation
+
+Fvec = ocv_carl(dvec, X, hvec);
+[Fmin,iOCV] = min(Fvec);
+
+%======================================================
+% PLOTTING
+
+lamL = hvec(iL); GL = Gvec(iL); rssL = rss(iL); mssL = mss(iL); kapL = kap(iL); FL = Fvec(iL);
+lamF = hvec(iOCV); GF = Gvec(iOCV); rssF = rss(iOCV); mssF = mss(iOCV); kapF = kap(iOCV); FF = Fvec(iOCV);
+lamG = hvec(iGCV); GG = Gvec(iGCV); rssG = rss(iGCV); mssG = mss(iGCV); kapG = kap(iGCV); FG = Fvec(iGCV);
+
+x1 = log10(rss);
+y1 = log10(mss);
+x2 = log10(hvec);
+y2 = kap;
+x3 = log10(hvec);
+y3 = log10(Fvec);
+x4 = log10(hvec);
+y4 = log10(Gvec);
+
+stx1 = ' Misfit norm, log10 RSS';
+sty1 = ' Model norm, log10 MSS';
+stx2 = ' Regularization parameter, log10 \lambda';
+sty2 = ' Curvature of L-curve, \kappa(\lambda)';
+stx3 = stx2;
+sty3 = ' OCV function, log10 F(\lambda)';
+stx4 = stx2;
+sty4 = ' GCV function, log10 G(\lambda)';
+
+%stfm = '%.4f';
+stfm = '%.2e';
+stlam_L = [' \lambda-L = ' num2str(sprintf(stfm, lamL))];
+stlam_ocv = [' \lambda-OCV = ' num2str(sprintf(stfm, lamF))];
+stlam_gcv = [' \lambda-GCV = ' num2str(sprintf(stfm, lamG))];
+
+%------------------------
+figure; nr=2; nc=2;
+msize = 8;
+nlab = 10; ilabs = round(linspace(1,nh,nlab));
+
+subplot(nr,nc,1); hold on;
+plot(x1,y1,'.');
+plot(log10(rssL),log10(mssL),'ko','markersize',8,'MarkerFaceColor','r');
+plot(log10(rssG),log10(mssG),'kV','markersize',8,'MarkerFaceColor','g');
+plot(log10(rssF),log10(mssF),'k^','markersize',8,'MarkerFaceColor','c');
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+axy = axis; dx = axy(2)-axy(1);
+for kk=1:nlab
+ ii = ilabs(kk);
+ text(x1(ii)+dx*0.02,y1(ii),[num2str(sprintf(stfm, hvec(ii)))],'fontsize',8,'color','b');
+end
+legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv);
+xlabel(stx1); ylabel(sty1); grid on;
+
+subplot(nr,nc,2); hold on;
+plot(x2,y2,'.');
+plot(log10(lamL),kapL,'ko','markersize',8,'MarkerFaceColor','r');
+plot(log10(lamG),kapG,'kV','markersize',8,'MarkerFaceColor','g');
+plot(log10(lamF),kapF,'k^','markersize',8,'MarkerFaceColor','c');
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northwest');
+xlabel(stx2); ylabel(sty2); grid on;
+
+subplot(nr,nc,3); hold on;
+plot(x3,y3,'.');
+plot(log10(lamL),log10(FL),'ko','markersize',8,'MarkerFaceColor','r');
+plot(log10(lamG),log10(FG),'kV','markersize',8,'MarkerFaceColor','g');
+plot(log10(lamF),log10(FF),'k^','markersize',8,'MarkerFaceColor','c');
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northeast');
+xlabel(stx3); ylabel(sty3); grid on;
+
+subplot(nr,nc,4); hold on;
+plot(x4,y4,'.');
+plot(log10(lamL),log10(GL),'ko','markersize',8,'MarkerFaceColor','r');
+plot(log10(lamG),log10(GG),'kV','markersize',8,'MarkerFaceColor','g');
+plot(log10(lamF),log10(GF),'k^','markersize',8,'MarkerFaceColor','c');
+axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
+legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northwest');
+xlabel(stx4); ylabel(sty4); grid on;
+
+orient tall, wysiwyg, fontsize(9)
+
+%======================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/ridge_carl.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/suf.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/suf.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/suf.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/suf.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,39 @@
+function suffix=suf(string,delim,iter)
+% suffix=SUF(string,delim,iter)
+%
+% Finds the suffix or extension in a string, say a filename.
+%
+% INPUT:
+%
+% string The string to be parsed
+% delim The string that separates prefix from suffix
+% iter 0 Don't iterate
+% 1 Iterate on the suffix to get the last bit
+%
+% See also PREF
+%
+% Last modified by fjsimons-at-alum.mit.edu, 05/06/2008
+
+defval('delim','.');
+defval('iter',1)
+
+ldel=length(delim);
+[prefix,suffix]=strtok(string,delim);
+
+% This added 17.11.2004
+if iscell(suffix)
+ suffix=cell2mat(suffix);
+end
+
+suffix=suffix(ldel+1:end);
+
+if iter==1
+ % This suffix might not be the last one
+ while findstr(suffix,delim)
+ suffix=suf(suffix,delim);
+ end
+end
+
+
+
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/suf.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/wysiwyg.m (from rev 15485, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/wysiwyg.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/wysiwyg.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/wysiwyg.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,17 @@
+function wysiwyg
+%WYSIWYG -- this function is called with no args and merely
+% changes the size of the figure on the screen to equal
+% the size of the figure that would be printed,
+% according to the papersize attribute. Use this function
+% to give a more accurate picture of what will be
+% printed.
+% Dan(K) Braithwaite, Dept. of Hydrology U.of.A 11/93
+
+unis = get(gcf,'units');
+ppos = get(gcf,'paperposition');
+set(gcf,'units',get(gcf,'paperunits'));
+pos = get(gcf,'position');
+pos(3:4) = ppos(3:4);
+set(gcf,'position',pos);
+set(gcf,'units',unis);
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/wysiwyg.m
___________________________________________________________________
Name: svn:mergeinfo
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/year.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/year.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/year.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,24 @@
+function y = year(d)
+%YEAR Year of date.
+% Y = YEAR(D) returns the year of a serial date number or a date string, D.
+%
+% For example, y = year(728647) or y = year('19-Dec-1994') returns y = 1994.
+%
+% See also DATEVEC, DAY, MONTH.
+
+% Copyright 1995-2006 The MathWorks, Inc.
+% $Revision: 1.6.2.2 $ $Date: 2006/04/20 17:47:27 $
+
+if nargin < 1
+ error('finance:year:missingInputs','Please enter D.')
+end
+if ischar(d)
+ d = datenum(d);
+end
+
+c = datevec(d(:)); % Generate date vectors from dates
+y = c(:,1); % Extract years
+if ~ischar(d)
+ y = reshape(y,size(d));
+end
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/year.m
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/nounder.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/nounder.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/nounder.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,9 +0,0 @@
-function sin=nounder(sin)
-% sin=NOUNDER(sin)
-%
-% Changes underscores to dashes in a string
-%
-% Last modified by fjsimons-at-alum.mit.edu, Feb 06th, 2004
-
-sin(find(abs(sin)==95))='-';
-
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ocv_carl.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ocv_carl.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ocv_carl.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,87 +0,0 @@
-%
-% function rss_vec = ocv_carl(g, A, lamvec)
-% Carl Tape, 30-March-2006
-%
-% Copied from gcv_carl.m on 27-March-2006.
-%
-% Returns the ordinary cross-validation (OCV) function corresponding to a
-% set of input regularization (or damping) parameters.
-%
-% TWO ALGORITHMS ARE SHOWN:
-% (1) Brute force, which involves ndata inversions per lambda
-% (2) Elegant formalisum, which involves one inversion per lambda
-% --> See Latex notes
-% /home/carltape/classes/acm118/2006_handouts/hw3/hw3sol_2006_prob3.pdf
-%
-% Using some GPS data, I checked that these approaches are identical.
-%
-% calls xxx
-% called by ridge_carl.m
-%
-
-function rss_vec = ocv_carl(d, A, lamvec)
-
-% Size of inputs
-[ndata, nparm] = size(A);
-numlam = length(lamvec);
-
-if (min(lamvec) < 0)
- error('Impossible regularization parameter lambda.')
-end
-
-rss_vec = zeros(numlam,1);
-
-% loop over regularization parameters
-for ii=1:numlam
- lam = lamvec(ii);
- disp([' ii = ' num2str(ii) '/' num2str(numlam) ', lam = ' num2str(lam)]);
-
- if 1==1
- H = A*inv(A'*A + lam^2*eye(nparm))*A';
- dhat = H*d;
-
- % OCV residual
- res = (d - dhat) ./ (1 - diag(H));
-
- % sum the residuals
- rss_vec(ii) = sum(res.^2);
-
- else
- % loop over datapoints
- for jj=1:ndata
- %disp([' jj = ' num2str(jj) '/' num2str(ndata) ]);
-
- % indices for which you compute the model parameters
- switch jj
- case 1, einds = [2:ndata];
- case ndata, einds = [1:ndata-1];
- otherwise, einds = [1:jj-1 jj+1:ndata];
- end
-
- % indices to estimate the RSS
- oinds = jj;
-
- % reduced matrices
- X = A(einds,:);
- g = d(einds);
-
- % note: regularization matrix is identity matrix
- f_h = inv(X'*X + lam^2*eye(nparm))*X'*g;
-
- % estimate the model at the datapoints NOT used
- res = A(oinds,:) * f_h - d(oinds);
-
- % sum the residuals
- rss_vec(ii) = rss_vec(ii) + sum(res.^2);
- end
- end
-
-end
-
-rss_vec = rss_vec / ndata; % normalization
-
-figure; loglog(lamvec, rss_vec, '.'); grid on;
-xlabel(' Regularization parameter, log (\lambda)');
-ylabel(' RSS at datapoints NOT used in computing model');
-
-%======================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/osdep.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/osdep.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/osdep.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,16 +0,0 @@
-function of=osdep
-% of=OSDEP
-%
-% Returns the value of the read option to be used
-% on the LOCAL operating system for files created
-% on the LOCAL operating system.
-%
-% Last modified by fjsimons-at-alum.mit.edu, 23.11.2004
-
-
-if strcmp(getenv('OSTYPE'),'linux')
- of= 'l';
-end
-if strcmp(getenv('OSTYPE'),'solaris')
- of= 'b';
-end
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/plot_histo.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/plot_histo.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/plot_histo.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,22 +0,0 @@
-%
-% function N = plot_histo(hdat,edges)
-% CARL TAPE, 30-May-2007
-% printed xxx
-%
-% This function plots a histogram with cyan bars and black boundaries.
-%
-% calls xxx
-% called by xxx
-%
-
-function N = plot_histo(hdat,edges)
-
-Ntotal = length(hdat);
-[N,bin] = histc(hdat,edges);
-bar(edges,N/Ntotal,'histc');
-xlim([min(edges) max(edges)]);
-ylabel([' Fraction (N = ' num2str(Ntotal) ')']);
-
-h = findobj(gca,'Type','patch'); set(h,'FaceColor',[0 1 1],'EdgeColor','k');
-
-%=======================================================================
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readCMT.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readCMT.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readCMT.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,85 +0,0 @@
-%
-% readCMT.m
-% CARL TAPE, 26-June-2007
-% printed xxx
-%
-% This file reads a concatenated CMTSOLUTION file and outputs the data.
-% The CMT tshift parameter is added to the origin time, which is the
-% time used for the initial CMT solution.
-%
-% INPUT:
-% filename file containing all CMTSOLUTION files concatenated together
-% nlines_cmt number of lines per CMT solution (13)
-% nspace_between_entries number of spaces between CMTSOLUTION blocks (0 or 1)
-%
-% OUTPUT:
-% M moment tensor in N-m units
-%
-% calls xxx
-% called by slab_readCMT.m, socal_quakes_2007.m
-%
-
-function [date,tshift,hdur,lat,lon,dep,M,eid,elabel] = ...
- readCMT(filename, nlines_cmt, nspace_between_entries)
-
-% read in concatenated CMTSOLUTION files
-lines = textread(filename,'%s','delimiter','\n');
-nlines = length(lines);
-
-% number of events
-enum = (nlines + nspace_between_entries) / (nlines_cmt + nspace_between_entries);
-if mod(enum,1) ~= 0, error(' enum should be an integer'); end
-disp([' File : ' filename ]);
-disp([' Number of events : ' num2str(enum) ]);
-
-% initialize vectors
-date = zeros(enum,1); tshift = zeros(enum,1); hdur = zeros(enum,1);
-lat = zeros(enum,1); lon = zeros(enum,1); dep = zeros(enum,1);
-Mrr = zeros(enum,1); Mtt = zeros(enum,1); Mpp = zeros(enum,1);
-Mrt = zeros(enum,1); Mrp = zeros(enum,1); Mtp = zeros(enum,1);
-
-for kk = 1:enum
- in1 = (kk-1)*(nlines_cmt + nspace_between_entries);
-
- % first full line of CMTSOLUTION, as a string
- ltemp = lines{in1+1};
-
- if length(ltemp) == 0, error('the CMTSOLUTION line is length zero'); end
-
- % replace W with space for more recent CMT solutions
- % SWEQ2006
- % PDEW2006
- if or(strcmp(ltemp(4),'W'),strcmp(ltemp(4),'Q')), ltemp(4) = ' '; end
-
- [j1,yr,mo,dy,hr,min,sec,lat_pde(kk),lon_pde(kk),dep_pde(kk),j5,j6,name1,name2,name3,name4,name5,name6] = ...
- strread(ltemp,'%s%f%f%f%f%f%f%f%f%f%f%f%s%s%s%s%s%s');
- [j1,j2,eid(kk)] = strread(lines{in1+2},'%s%s%s');
- [j1,j2,tshift(kk)] = strread(lines{in1+3},'%s%s%f');
- [j1,j2,hdur(kk)] = strread(lines{in1+4},'%s%s%f');
- [j1,lat(kk)] = strread(lines{in1+5},'%s%f');
- [j1,lon(kk)] = strread(lines{in1+6},'%s%f');
- [j1,dep(kk)] = strread(lines{in1+7},'%s%f');
- [j1,Mrr(kk)] = strread(lines{in1+8},'%s%f');
- [j1,Mtt(kk)] = strread(lines{in1+9},'%s%f');
- [j1,Mpp(kk)] = strread(lines{in1+10},'%s%f');
- [j1,Mrt(kk)] = strread(lines{in1+11},'%s%f');
- [j1,Mrp(kk)] = strread(lines{in1+12},'%s%f');
- [j1,Mtp(kk)] = strread(lines{in1+13},'%s%f');
-
- % NOTE the inclusion of the time shift here for the CMT solution
- date(kk) = datenum(yr,mo,dy,hr,min, sec + tshift(kk));
-
- elabel{kk} = [char(name1) ' ' char(name2) ' ' char(name3) ...
- ' ' char(name4) ' ' char(name5) ' ' char(name6)];
-end
-
-eid = eid(:);
-elabel = elabel(:);
-
-% moment tensor elements
-M = [Mrr Mtt Mpp Mrt Mrp Mtp];
-
-% convert moment tensor from dyne-cm to N-m
-M = 1e-7 * M;
-
-%======================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_station_SPECFEM.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_station_SPECFEM.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_station_SPECFEM.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,23 +0,0 @@
-%
-% function [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(filename)
-% CARL TAPE, 15-March-2007
-% printed xxx
-%
-% This function inputs a STATIONS file formatted for SPECFEM3D.f90.
-% See also write_station_SPECFEM.m
-%
-% In Perl, this can be done as follows:
-% open(IN,${stations_seis}); @temp = <IN>; $nrec = $temp[0]; chomp($nrec);
-% print CSH "tail -n $nrec ${stations_seis} > stemp0 \n";
-% print CSH "awk '{print \$4,\$3}' stemp0 > temp \n";
-%
-% calls xxx
-% called by xxx
-%
-
-function [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(filename)
-
-% read in the STATIONS file
-[stnm,netwk,rlat,rlon,relev,rburial] = textread(filename,'%s%s%f%f%f%f','headerlines',1);
-
-%=======================================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,32 +0,0 @@
-%
-% 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,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,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(:);
-
-%======================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/read_window_chi_all.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,128 +0,0 @@
-%
-% 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
-
-%======================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readsac.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readsac.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/readsac.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,110 +0,0 @@
-function varargout=readsac(filename,plotornot,osd)
-% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename)
-% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename,1)
-% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename,1,'l')
-% [SeisData,HdrData,tnu,pobj,tims]=READSAC(filename,1,'b')
-%
-% READSAC(...)
-%
-% Reads SAC-formatted data (and makes a plot)
-%
-% INPUT:
-%
-% filename The filename, full path included
-% plotornot 1 Makes a plot
-% 0 Does not make a plot [default]
-% osd 'b' for data saved on Solaris read into Linux
-% 'l' for data saved on Linux read into Linux
-%
-% OUPUT:
-%
-% SeisData The number vector
-% HdrData The header structure array
-% tnu The handle to the plot title; or appropriate string
-% pobj The handle to the plot line and the xlabel
-% tims The times on the x-axis
-%
-% See SETVAR.PL out of SACLIB.PM
-%
-% Last modified by fjsimons-at-alum.mit.edu, 07/07/2008
-
-defval('plotornot',0)
-defval('osd',osd)
-
-fid=fopen(filename,'r',osd);
-if fid==-1
- error([ 'File ',filename,' does not exist in current path ',pwd])
-end
-HdrF=fread(fid,70,'float32');
-HdrN=fread(fid,15,'int32');
-HdrI=fread(fid,20,'int32');
-HdrL=fread(fid,5,'int32');
-HdrK=str2mat(fread(fid,[8 24],'char'))';
-SeisData=fread(fid,HdrN(10),'float32');
-fclose(fid);
-
-% Meaning of LCALDA from http://www.llnl.gov/sac/SAC_Manuals/FileFormatPt2.html
-% LCALDA is TRUE if DIST, AZ, BAZ, and GCARC are to be calculated from
-% station and event coordinates. But I don't know how to set this
-% variable to a logical yet using setvar.pl. What is the TRUE/FALSE
-% format in SAC? 18.1.2006
-
-% If you change any of this, change it in WRITESAC as well!
-HdrData=struct(...
- 'AZ',HdrF(52),...
- 'B',HdrF(6),...
- 'BAZ',HdrF(53),...
- 'DIST',HdrF(51),...
- 'DELTA',HdrF(1),...
- 'E',HdrF(7),...
- 'EVDP',HdrF(39),...
- 'EVEL',HdrF(38),...
- 'EVLA',HdrF(36),...
- 'EVLO',HdrF(37),...
- 'GCARC',HdrF(54),...
- 'IFTYPE',HdrI(1),...
- 'KCMPNM',HdrK(21,HdrK(21,:)>64),...
- 'KINST',HdrK(24,HdrK(24,:)>64),...
- 'KSTNM',HdrK(1,HdrK(1,:)>64),...
- 'KUSER0',HdrK(18,HdrK(18,:)>64),...
- 'LCALDA',HdrL(4),...
- 'MAG',HdrF(40),...
- 'NPTS',HdrN(10),...
- 'NVHDR',HdrN(7),...
- 'NZHOUR',HdrN(3),...
- 'NZJDAY',HdrN(2),...
- 'NZMIN',HdrN(4),...
- 'NZMSEC',HdrN(6),...
- 'NZSEC',HdrN(5),...
- 'NZYEAR',HdrN(1),...
- 'SCALE',HdrF(4),...
- 'STDP',HdrF(35),...
- 'STEL',HdrF(34),...
- 'STLA',HdrF(32),...
- 'STLO',HdrF(33),...
- 'T0',HdrF(11),...
- 'T1',HdrF(12),...
- 'T2',HdrF(13),...
- 'T3',HdrF(14));
-
-tims=linspace(HdrData.B,HdrData.E,HdrData.NPTS);
-
-if plotornot==1
- pobj=plot(tims,SeisData);
- filename(find(abs(filename)==95))='-';
- tito=nounder(suf(filename,'/'));
- if isempty(tito)
- tito=nounder(filename);
- end
- tnu=title(tito);
- axis tight
- pobj(2)=xlabel([ 'Time (s)']);
-else
- tnu=nounder(suf(filename,'/'));
- pobj=0;
-end
-
-nam={'SeisData' 'HdrData' 'tnu' 'pobj','tims'};
-for index=1:nargout
- varargout{index}=eval(nam{index});
-end
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ridge_carl.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ridge_carl.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/ridge_carl.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,227 +0,0 @@
-%
-% function
-% Carl Tape (Tapio Schneider, ACM 118)
-% 06-Nov-2006
-%
-% This function inputs a design matrix, a data vector, and a vector of
-% regularization parameters, and it returns three different curves that may
-% be used to select the best parameter:
-% (1) L-curve and curvature
-% (2) generalized cross-validation curve (GCV)
-% (3) ordinary cross-validation (OCV), also known as 'leave-one-out' CV
-%
-% It is best to input a large number of regularization parameters, so that
-% the min and max of the respective functions can be easily obtained.
-%
-% This program is copied in part from ridge_tapio.m
-%
-% NOTE THE PLOTTING OPTIONS.
-%
-%-------------------------------------------------
-% RIDGE Ridge regression estimates.
-%
-% Given a vector g, a design matrix X, and
-% a regularization parameter h,
-%
-% [m, rss, mss, dof] = ridge_tapio(g, X, h)
-%
-% returns the ridge regression estimate of the vector f in the
-% linear regression model
-%
-% g = X*f + noise.
-%
-% Also returned are the residual sum of squares rss, the sum of
-% squares mss of the elements of the ridge regression estimate
-% m (the squared norm of m), and the effective number of
-% residual degrees of freedom dof.
-%
-% If h is a vector of regularization parameters, the i-th column
-% m(:,i) is the ridge regression estimate for the regularization
-% parameter h(i); the i-th elements of rss and mss are the
-% associated residual sum of squares and estimate sum of squares.
-%
-% If no regularization parameter h is given, generalized
-% cross-validation is used to determine the regularization
-% parameter. The chosen regularization parameter h and the value of
-% the GCV function are then returned as the fifth and sixth
-% output arguments
-%
-% [m, rss, mss, dof, h, G] = ridge_tapio(g, X);
-%
-% Adapted from various routines in Per Christian Hansen's
-% Regularization Toolbox.
-%
-% calls gcvfctn.m, curvature.m
-% called by xxx
-%
-
-function [m, rss, mss, Gvec, Fvec, dof, kap, iL, iGCV, iOCV] = ridge_carl(dvec, X, hvec)
-
-% Size of inputs
-[n, p] = size(X);
-q = min(n, p);
-nh = length(hvec);
-if (min(hvec) < 0)
- error('Impossible regularization parameter h.')
-end
-
-% Initialize outputs
-m = zeros(p, nh);
-rss = zeros(nh, 1);
-mss = zeros(nh, 1);
-dof = zeros(nh, 1);
-
-% Compute SVD of X
-[U, S, V] = svd(X, 0);
-s = diag(S); % vector of singular values
-s2 = s.^2;
-
-% Coefficients in expansion of solution in terms of right singular vectors
-fc = U(:, 1:q)'*dvec;
-zeta = s .* fc;
-
-% Treat each regularization parameter separately
-for j = 1:nh
- m(:, j) = V(:, 1:q) * (zeta ./ (s2 + hvec(j)^2));
- mss(j) = sum(m(:, j).^2);
- rss(j) = hvec(j)^4 * sum(fc.^2 ./ (s2 + hvec(j)^2).^2);
- dof(j) = n - sum(s2./(s2 + hvec(j)^2));
-end
-
-% In overdetermined case, add rss of least-squares problem
-if (n > p)
- rss = rss + sum((dvec - U(:, 1:q)*fc).^2);
-end
-
-%-----------------------
-% determine the Lcurve pick (max curvature)
-
-x1 = log10(rss);
-y1 = log10(mss);
-
-% % smooth curvature interpolation to get h_L
-% num = 1000;
-% xsmooth = linspace(x1(1),x1(end),1000);
-% ysmooth = interp1(x1,y1,xsmooth,'cubic');
-% [i0,kap_smooth] = curvature(xsmooth,ysmooth);
-% rss_L = 10^xsmooth(i0);
-% mss_L = 10^ysmooth(i0);
-% h_L = 10^interp1(x1,log10(hvec),xsmooth(i0),'cubic');
-
-% curvature, based on input h values alone
-[iL,kap] = curvature(x1,y1);
-%h_L = hvec(iL);
-%rss_L = 10^x1(iL);
-%mss_L = 10^y1(iL);
-
-%-----------------------
-% obtain GCV `best' solution and GCV curve
-
-% GCV minimum -- 'exact' in the sense a minimization method is used
-% [hmin, Gmin] = gcv_tapio(U, s, dvec, 'ridge');
-
-% GCV minimum -- 'crude' in the sense that we coarsely sample the function
-dof0 = n-q;
-rss0 = sum((dvec - U(:, 1:q)*fc).^2);
-Gvec = zeros(nh,1);
-for j = 1:nh
- Gvec(j) = gcvfctn(hvec(j), s2, fc, rss0, dof0);
-end
-[Gmin,iGCV] = min(Gvec);
-%hmin = hvec(iGCV);
-
-% GCV best model and L-curve point for h_GCV (= hmin)
-%mod_min = inv(X'*X + hmin^2*eye(p))*X'*dvec;
-%res = X*mod_min - dvec;
-%rss_min = sum(res.^2);
-%mss_min = sum(mod_min.^2);
-
-% compute G for the Lcurve pick
-%G_L = gcvfctn(h_L, s2, fc, rss0, dof0);
-
-%-----------------------
-% ordinary (leave-one-out) cross-validation
-
-Fvec = ocv_carl(dvec, X, hvec);
-[Fmin,iOCV] = min(Fvec);
-
-%======================================================
-% PLOTTING
-
-lamL = hvec(iL); GL = Gvec(iL); rssL = rss(iL); mssL = mss(iL); kapL = kap(iL); FL = Fvec(iL);
-lamF = hvec(iOCV); GF = Gvec(iOCV); rssF = rss(iOCV); mssF = mss(iOCV); kapF = kap(iOCV); FF = Fvec(iOCV);
-lamG = hvec(iGCV); GG = Gvec(iGCV); rssG = rss(iGCV); mssG = mss(iGCV); kapG = kap(iGCV); FG = Fvec(iGCV);
-
-x1 = log10(rss);
-y1 = log10(mss);
-x2 = log10(hvec);
-y2 = kap;
-x3 = log10(hvec);
-y3 = log10(Fvec);
-x4 = log10(hvec);
-y4 = log10(Gvec);
-
-stx1 = ' Misfit norm, log10 RSS';
-sty1 = ' Model norm, log10 MSS';
-stx2 = ' Regularization parameter, log10 \lambda';
-sty2 = ' Curvature of L-curve, \kappa(\lambda)';
-stx3 = stx2;
-sty3 = ' OCV function, log10 F(\lambda)';
-stx4 = stx2;
-sty4 = ' GCV function, log10 G(\lambda)';
-
-%stfm = '%.4f';
-stfm = '%.2e';
-stlam_L = [' \lambda-L = ' num2str(sprintf(stfm, lamL))];
-stlam_ocv = [' \lambda-OCV = ' num2str(sprintf(stfm, lamF))];
-stlam_gcv = [' \lambda-GCV = ' num2str(sprintf(stfm, lamG))];
-
-%------------------------
-figure; nr=2; nc=2;
-msize = 8;
-nlab = 10; ilabs = round(linspace(1,nh,nlab));
-
-subplot(nr,nc,1); hold on;
-plot(x1,y1,'.');
-plot(log10(rssL),log10(mssL),'ko','markersize',8,'MarkerFaceColor','r');
-plot(log10(rssG),log10(mssG),'kV','markersize',8,'MarkerFaceColor','g');
-plot(log10(rssF),log10(mssF),'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
-axy = axis; dx = axy(2)-axy(1);
-for kk=1:nlab
- ii = ilabs(kk);
- text(x1(ii)+dx*0.02,y1(ii),[num2str(sprintf(stfm, hvec(ii)))],'fontsize',8,'color','b');
-end
-legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv);
-xlabel(stx1); ylabel(sty1); grid on;
-
-subplot(nr,nc,2); hold on;
-plot(x2,y2,'.');
-plot(log10(lamL),kapL,'ko','markersize',8,'MarkerFaceColor','r');
-plot(log10(lamG),kapG,'kV','markersize',8,'MarkerFaceColor','g');
-plot(log10(lamF),kapF,'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
-legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northwest');
-xlabel(stx2); ylabel(sty2); grid on;
-
-subplot(nr,nc,3); hold on;
-plot(x3,y3,'.');
-plot(log10(lamL),log10(FL),'ko','markersize',8,'MarkerFaceColor','r');
-plot(log10(lamG),log10(FG),'kV','markersize',8,'MarkerFaceColor','g');
-plot(log10(lamF),log10(FF),'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
-legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northeast');
-xlabel(stx3); ylabel(sty3); grid on;
-
-subplot(nr,nc,4); hold on;
-plot(x4,y4,'.');
-plot(log10(lamL),log10(GL),'ko','markersize',8,'MarkerFaceColor','r');
-plot(log10(lamG),log10(GG),'kV','markersize',8,'MarkerFaceColor','g');
-plot(log10(lamF),log10(GF),'k^','markersize',8,'MarkerFaceColor','c');
-axis tight; ax1 = axis; axis(axes_expand(ax1,1.1));
-legend(' \lambda',stlam_L,stlam_gcv,stlam_ocv,'location','northwest');
-xlabel(stx4); ylabel(sty4); grid on;
-
-orient tall, wysiwyg, fontsize(9)
-
-%======================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,462 +0,0 @@
-%
-% subspace_specfem.m
-% CARL TAPE, 01-Feb-2009
-% printed xxx
-%
-% This program assumes you have FIRST run compute_misfit.m, and also have
-% generated the corresponding Hessian matrix to load here.
-%
-% See also wave2d_subspace.m
-%
-% calls hfile2hess.m
-% called by xxx
-%
-
-clear
-close all
-format short
-format compact
-
-ax1 = [-121 -114 31 37];
-stfm = '%4.4i';
-
-xlab1 = 'p, singular value index';
-xlab2 = 'p, singular value truncation index';
-ylab1 = 'log10( singular value )';
-ylab2 = 'log10( misfit : dot[ d - G*dm(p), d - G*dm(p) ] )';
-
-imod = input(' Enter the current model number (0, 1, ...): ');
-iwrite = input(' Enter 1 to write files (0 otherwise) : ');
-
-stmod = sprintf('m%2.2i',imod);
-
-% directory containing the Hessian file that was output from the cluster
-sthess = ['/net/sierra/raid1/carltape/results/HESSIANS/'];
-
-% load the data and event weights (compute_misfit.m)
-idatacov = input(' Enter idatacov (1--by event, 2--by window, 3--none) : ');
-stdcovs = {'event','window','none'};
-stdcov = stdcovs{idatacov};
-ftag = [stmod '_' stdcov];
-odir0 = ['OUTPUT_SUBSPACE/' stmod '/' stdcov '/'];
-load([odir0 ftag '_data_norms']);
-nsrc = length(eids);
-
-%----------------------------------------------
-% check the slick operation used in subspace_update.f90
-
-if 0==1
- npar = 10;
- nsrc = 4;
- G = rand(npar,nsrc);
- d = rand(nsrc,1);
-
- x = zeros(npar,1);
- for isrc=1:nsrc % loop over kernels (hundreds)
- for ipar=1:npar % loop over gridpoints (millions)
- x(ipar) = x(ipar) + G(ipar,isrc) * d(isrc);
- end
- end
- G*d, x
-end
-
-%----------------------------------------------
-
-% construct the matrix of data-covariance normalization
-Cd_fac = zeros(nsrc,nsrc);
-for i = 1:nsrc
- for j = 1:nsrc
- Cd_fac(i,j) = dcov_fac_e(i) * dcov_fac_e(j);
- i_ind(i,j) = i;
- j_ind(i,j) = j;
- end
-end
-
-% % OBSOLETE: construct the matrix of weight factors and data-covariance normalization
-% if 0==1
-% Ws = zeros(nsrc,nsrc);
-% Cd_fac = zeros(nsrc,nsrc);
-% for i = 1:nsrc
-% for j = 1:nsrc
-% Ws(i,j) = ws(i) * ws(j);
-% Cd_fac(i,j) = dcov_fac_e(i) * dcov_fac_e(j);
-% i_ind(i,j) = i;
-% j_ind(i,j) = j;
-% end
-% end
-% figure; pcolor(i_ind,j_ind,Ws); shading flat;
-% xlabel('Row index'); ylabel('Column index');
-% title('Symmetric matrix of data weight terms');
-% %caxis([0 200]);
-% colorbar; axis equal; axis tight;
-% end
-
-% load the event IDs corresponding to the kernels
-%eids = load([sthess 'kernels_done_mod']);
-%neid = length(eids);
-
-kmin = 6;
-kmax = 6;
-nker = kmax - kmin + 1;
-stkers = {'mu_kernel','kappa_kernel','both','mu_kernel_smooth','kappa_kernel_smooth','both_smooth'};
-stkerlabs = {'unsmoothed beta kernel','unsmoothed beta kernel','beta + bulk','smoothed beta kernel','smoothed bulk kernel','smoothed beta + bulk'};
-
-%stwts = {'with data weights'};
-%stwts = {'with data weights','with NO data weights'};
-nker = length(stkers);
-%stcs = {'b.','r.'};
-%iw = 1;
-
-%for iw = 1:length(stwts)
-
-for iker = kmin:kmax
- stker = stkers{iker};
- stkerlab = stkerlabs{iker};
- stlab = [ftag '_' stker];
-
- dir1 = [odir0 stker '/'];
-
- if ~any(iker == [3 6])
- % get the Hessian from the subspace_hessian.f90 output
- hfile = sprintf([sthess 'hessian_index_' stmod '_' stker '_%3.3i'],nsrc);
- H0 = hfile2hess(hfile,nsrc);
- H = H0;
-
- % include the weights determined from the data
- %if iw == 1
- % H = H .* Ws;
- %end
-
- % When we computed the kernels, we did not include the
- % normalization factor based on the number of windows picked
- % that is included within the data covariance (and used in computing Ws).
- % Thus we compute the MATRIX Cd_fac above to account for this.
- H = H ./ Cd_fac;
-
- % put the matrix onto order-1 values for computation purposes
- %Hnorm = norm(H); H = H / norm(H);
-
- else
- hfile_beta = sprintf([sthess 'hessian_index_' stmod '_' stkers{iker-2} '_%3.3i'],nsrc);
- H0_beta = hfile2hess(hfile_beta,nsrc);
- hfile_bulk = sprintf([sthess 'hessian_index_' stmod '_' stkers{iker-1} '_%3.3i'],nsrc);
- H0_bulk = hfile2hess(hfile_bulk,nsrc);
-
- %H_beta = H0_beta .* Ws ./ Cd_fac;
- %H_bulk = H0_bulk .* Ws ./ Cd_fac;
- H_beta = H0_beta ./ Cd_fac;
- H_bulk = H0_bulk ./ Cd_fac;
- H = H_beta + H_bulk;
-
- dtemp = [[1:nsrc]' eids_num diag(H_beta) diag(H_bulk) diag(H) diag(H_bulk)./diag(H_beta) Ns];
- dtemp_sort1 = sortrows(dtemp,-6); % sort by relative importance of bulk-to-shear
- dtemp_sort2 = sortrows(dtemp,-5); % sort by diagonal
-
- % sort all entries of the Hessian -- note that some are negative
- % we take the upper traingular elements, excluding the diagonal
- Htri = triu(H,1); itri = triu(i_ind,1); jtri = triu(j_ind,1);
- Hcols = [Htri(:) itri(:) jtri(:) ];
- [Hsort,iHsort_good] = sortrows(Hcols, -1);
- [Hsort,iHsort_bad] = sortrows(Hcols, 1);
-
- if iwrite == 1
- filename = [dir1 'hessian_bulk-shear_' stlab];
- fid = fopen(filename,'w');
- fprintf(fid,'Hessian diagonal contributions from beta and bulk:\n');
- fprintf(fid,'%6s%10s%10s%10s%10s%10s%6s\n','index','eid','beta','bulk','total','bulk/beta','Nwin');
- for ix = 1:nsrc
- fprintf(fid,'%6i%10i%10.2e%10.2e%10.2e%10.4f%6i\n',dtemp_sort1(ix,:));
- end
- fclose(fid);
-
- filename = [dir1 'hessian_overall_' stlab];
- fid = fopen(filename,'w');
- fprintf(fid,'Hessian diagonal contributions from beta and bulk:\n');
- fprintf(fid,'%6s%10s%10s%10s%10s%10s%6s\n','index','eid','beta','bulk','total','bulk/beta','Nwin');
- for ix = 1:nsrc
- fprintf(fid,'%6i%10i%10.2e%10.2e%10.2e%10.4f%6i\n',dtemp_sort2(ix,:));
- end
- fclose(fid);
-
- % make a list of the largest N positive and negative Hessian elements
- nprint = 100; % must be less than nsrc*nsrc
- filename = [dir1 'hessian_elements_good_' stlab];
- fid = fopen(filename,'w');
- fprintf(fid,'Largest positive non-diagonal elements of Hessian:\n');
- fprintf(fid,'%10s%6s%6s%12s%12s\n','Hij','i','j','eid-i','eid-j');
- for ix = 1:nprint
- iH = iHsort_good(ix);
- fprintf(fid,'%10.2e%6i%6i%12s%12s\n',Hcols(iH,:),eids{Hcols(iH,2)},eids{Hcols(iH,3)});
- end
- fclose(fid);
-
- filename = [dir1 'hessian_elements_bad_' stlab];
- fid = fopen(filename,'w');
- fprintf(fid,'Largest negative non-diagonal elements of Hessian:\n');
- fprintf(fid,'%10s%6s%6s%12s%12s\n','Hij','i','j','eid-i','eid-j');
- for ix = 1:nprint
- iH = iHsort_bad(ix);
- fprintf(fid,'%10.2e%6i%6i%12s%12s\n',Hcols(iH,:),eids{Hcols(iH,2)},eids{Hcols(iH,3)});
- end
- fclose(fid);
- end
-
- disp(' gradient balance for the Hessian (subspace) inversion (mean of last column) :');
- disp(mean( diag(H_bulk)./diag(H_beta) ))
- end
-
- %H = H + eye(nsrc);
- % mutemp = inv(H + eye(nsrc)) * dnorm;
-
- dH = diag(H);
- trH = sum(dH);
- trHH = sum(diag(H'*H));
-
- %------------------------
-
- pinds = [1:nsrc]';
-
- disp(' properties of Hessian (min, median, mean(abs), max, std):');
- stH1 = sprintf('min %.1e, median %.1e, mean(abs) %.1e, max %.1e, std %.1e, sum %.1e',...
- min(H(:)), median(H(:)), mean(abs(H(:))), max(H(:)), std(H(:)), sum(H(:)) );
- stH2 = sprintf('DIAGONAL : min %.1e, median %.1e, mean(abs) %.1e, max %.1e, std %.1e, sum %.1e',...
- min(dH), median(dH), mean(abs(dH)), max(dH), std(dH), sum(dH) );
- disp(stH1);
- disp(stH2);
-
- % plot H
- Hfac = 0.5; % adjust to saturate the color scale as you wish
- figure; pcolor(i_ind,j_ind,H); shading flat;
- xlabel('Row index'); ylabel('Column index');
- title({['Hessian (symmetric matrix) for ' stkerlab],stH1,stH2});
- caxis([-1 1]*max(H(:))*Hfac), colormap('jet'), colorbar;
- %map = colormap('seis'); colormap(flipud(map));
- axis equal; axis tight;
- fontsize(10), orient tall, wysiwyg
-
- % sort the diagonal values of H
- Hdiag = diag(H);
- %Hdiag = diag(H0);
- [jk,isort] = sortrows([pinds eids_num Hdiag],-3);
- npick = min([20 nsrc]);
- npick = nsrc;
- disp(' '); disp([' first ' num2str(npick) ' sorted diagonal entries of H: ' stkerlab]);
- for ii = 1:npick
- jj = isort(ii);
- disp(sprintf('%3i %10i (%3i), Ne = %5i, Hkk = %10.4e',...
- ii,eids_num(jj),pinds(jj),Ns(jj),Hdiag(jj)));
- end
-
- [jtemp,isort] = sort(Hdiag,'descend');
-
- figure(20+iker); hold on;
- plot(pinds,log10(sort(Hdiag,'descend')),'ro');
- %plot(pinds,log10(sort(diag(H'*H),'descend')),'b.');
- plot([1 nsrc],log10(trH)*[1 1],'r--');
- %plot([1 nsrc],log10(trHH)*[1 1],'b--');
- %legend('diagonal element of H','diagonal element of HtH','trace (H)','trace (HtH)');
- legend('diagonal element of H','trace (H)');
- xlabel('sorted event index'); ylabel(' log10 H-kk');
- title(sprintf(' trace of H = %6.2e --- trace of HtH = %6.2e',trH,trHH));
- xlim([1 nsrc]); grid on;
- fontsize(11), orient tall, wysiwyg
-
- %----------------------------------------------
-
- itsvd = 0;
-
- if itsvd == 0
-
- % regularization choices
- numlam = 100;
- if 0==1
- minlampwr = -2.25; maxlampwr = 1;
- lampwr = linspace(minlampwr,maxlampwr,numlam);
- lamvec = 10.^lampwr * sum(dH); % based on trace of H
- %lamvec = 10.^lampwr * sqrt(sum(diag(H'*H))); % based on trace of H'*H
- else
- %minlampwr = -2; maxlampwr = 5; % m03
- minlampwr = -3; maxlampwr = 3; % m04 - m11
- minlampwr = -5; maxlampwr = 0; % m12
- lampwr = linspace(minlampwr,maxlampwr,numlam);
- lamvec = 10.^lampwr;
- end
-
- [f_h, rss, mss, Gvec, Fvec, dof, kap, iL, iGCV, iOCV] = ridge_carl(dnorm,H,lamvec);
- title({['model ' stmod],stkerlab,sprintf('iL = %i, iOCV = %i, iGCV = %i',iL,iOCV,iGCV)});
-
- mu_all = zeros(nsrc,numlam);
- mu_norm = zeros(numlam,1);
- mu_res = zeros(numlam,1);
- for ip = 1:numlam
- lam = lamvec(ip);
- %if ip==iGCV, disp('iGCV -- GCV lambda pick'); end
- %if ip==iOCV, disp('iOCV -- OCV lambda pick'); end
- %if ip==iL, disp('iL -- L-curve lambda pick'); end
- disp(sprintf('%i out of %i, lam = %.2e',ip,numlam,lam));
-
- mu_all(:,ip) = inv(H + lam*eye(nsrc))*dnorm;
- %mu_all(:,ip) = inv(H'*H + lam^2*eye(nsrc))*H'*dnorm;
-
- % f = inv(G'*diag(Wvec)*G + lam0^2*eye(ngrid*ndim))*G'*diag(Wvec)*d;
- mu_norm(ip) = norm(mu_all(:,ip));
- mu_res(ip) = norm( dnorm - H*mu_all(:,ip) );
- end
-
- ipick = iOCV; % KEY: select a damping parameter (iL, iOCV, iGCV)
- ltemp = lamvec(ipick)
- mtemp = mu_all(:,ipick);
- mtemp_norm = mu_norm(ipick)
-
- figure; nr=2; nc=2;
-
- subplot(nr,nc,1); hold on; grid on;
- plot( log10(lamvec), log10(mu_norm),'b.','markersize',10);
- plot( log10(lamvec(ipick)), log10(mu_norm(ipick)),'ro','markersize',10);
- xlabel('log10( lambda )'); ylabel('log10( norm of mu )');
-
- subplot(nr,nc,2); hold on; grid on;
- plot( log10( lamvec ), mu_res, '.');
- plot( log10( lamvec(ipick) ), mu_res(ipick),'ro','markersize',10);
- xlabel('log10( lambda )'); ylabel('dnorm residual');
-
- subplot(nr,nc,3); plot( dnorm, H*mtemp, '.');
- xlabel('dnorm obs'); ylabel('dnorm pred'); grid on; axis equal;
-
- subplot(nr,nc,4); stdres = std(dnorm - H*mtemp);
- plot_histo(dnorm - H*mtemp,[-4:0.5:4]*stdres);
- %plot( 100*(dnorm - H*mtemp)./dnorm ,'.');
- %ylabel('dnorm residual = 100(pred - obs)/obs');
- title(' residuals : d - H*mu');
- xlabel('event index'); grid on;
-
- orient tall, wysiwyg, fontsize(10)
-
-% iCV = input(' Enter 1 for OCV or 0 for GCV : ');
-% if iCV == 1
-% lam = lamvec(iOCV);
-% elseif iCV == 0
-% lam = lamvec(iGCV);
-% else
-% error(' iCV must be 1 or 0');
-% end
-% mu = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm; % KEY vector
-% stp = sprintf('lambda = %.3f',lam);
-
- % write mu vectors to file
- if iwrite == 1
- otag = ['mu_' stlab];
- odir = [dir1 'mu_all/'];
- for ip = 1:numlam
- lam = lamvec(ip);
- mutemp = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm;
- filename = sprintf([otag '_p%3.3i'],ip);
- fid = fopen([odir filename],'w');
- for ii = 1:nsrc
- fprintf(fid,'%24.12e\n',mutemp(ii));
- end
- fclose(fid);
- end
-
- fid = fopen([odir '/lambda_' stdcov],'w');
- for ip = 1:numlam
- lam = lamvec(ip);
- fprintf(fid,'%24.12e%24.12e%6i\n',lamvec(ip),1/lamvec(ip),ip);
- end
- fclose(fid);
- end
-
- else
-
- % analyze the singular values of H
- [U,S,V] = svd(H);
- s = diag(S);
-
- pmin = min([13 nsrc]);
- ifit = [(pmin+1):60];
- %ifit = pinds;
- [xf,yf,mf,bf,rms] = linefit(pinds(ifit), log10(s(ifit)));
-
- %subplot(nr,nc,iker);
- figure; hold on;
- plot(xf,yf,'g','linewidth',6);
- plot(pinds,log10(s),'b.','markersize',10);
- plot(pmin*[1 1],[min(log10(s)) max(log10(s))],'k','linewidth',2)
- grid on; xlabel(xlab1); ylabel(ylab1);
- title(['singular value spectra for Hessian : ' stker]);
- fontsize(10), orient tall, wysiwyg
-
- %-------------------------
-
- %tpinds = [1:10:100];
- tpinds = pinds;
- [mu_all,rss,f_r_ss] = tsvd(dnorm,H,tpinds);
-
- % norms of mu vectors
- mu_norm = zeros(nsrc,1);
- for ip = 1:nsrc
- mu_norm(ip) = norm(mu_all(:,ip));
- end
-
- figure;
- subplot(2,2,1); hold on; grid on;
- plot(pinds,log10(s),'b.','markersize',10);
- xlabel(xlab1); ylabel(ylab1);
- subplot(2,2,2); hold on; grid on;
- plot( tpinds, log10(rss),'b.','markersize',10)
- xlabel(xlab2); ylabel(ylab2);
- subplot(2,2,3); hold on; grid on;
- plot( tpinds, log10(mu_norm),'b.','markersize',10)
- xlabel(xlab2); ylabel('log10( norm of mu )');
- fontsize(10), orient tall, wysiwyg
-
- % write mu vectors to file
- if iwrite == 1
- odir = [dir1 'mu_all/'];
- for ip = 1:length(tpinds);
- pmax = tpinds(ip);
- filename = sprintf([odir '_p%3.3i'],pmax);
- fid = fopen([odir '/' filename],'w');
- for ii = 1:nsrc
- fprintf(fid,'%24.12e\n',mu_all(ii,ip));
- end
- fclose(fid);
- end
- end
- end
-
-
-end
-
-%end
-
-% write misfit values to file -- these do not depend on the choice of model
-% variable used in constructing the Hessian
-if iwrite == 1
- % write data-norm vector to file
- fid = fopen([odir0 ftag '_dmisfit'],'w');
- fprintf(fid,'%24.12e\n',dmisfit);
- fclose(fid);
-
- % write data-norm vector to file
- fid = fopen([odir0 ftag '_nwin_tot'],'w');
- fprintf(fid,'%.0f\n',N);
- fclose(fid);
-
- % write data-norm vector to file
- fid = fopen([odir0 ftag '_data_norm'],'w');
- for isrc = 1:nsrc
- fprintf(fid,'%24.12e\n',dnorm(isrc));
- end
- fclose(fid);
-
- % write the factors for the data covariance (nevent x 1)
- % WRITE THIS IN DOUBLE PRECISION, JUST FOR SIMPLICITY
- fid = fopen([odir0 ftag '_dcov_fac'],'w');
- for isrc = 1:nsrc
- fprintf(fid,'%24.12e\n',dcov_fac_e(isrc));
- end
- fclose(fid);
-end
-
-%======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -0,0 +1,465 @@
+%
+% subspace_specfem.m
+% CARL TAPE, 01-Feb-2009
+% printed xxx
+%
+% This program assumes you have FIRST run compute_misfit.m, and also have
+% generated the corresponding Hessian matrix to load here.
+%
+% See also wave2d_subspace.m
+%
+% calls hfile2hess.m
+% called by xxx
+%
+
+clear
+close all
+format short
+format compact
+
+% add path to additional matlab scripts
+path(path,'./matlab_scripts')
+
+ax1 = [-121 -114 31 37];
+stfm = '%4.4i';
+
+xlab1 = 'p, singular value index';
+xlab2 = 'p, singular value truncation index';
+ylab1 = 'log10( singular value )';
+ylab2 = 'log10( misfit : dot[ d - G*dm(p), d - G*dm(p) ] )';
+
+imod = input(' Enter the current model number (0, 1, ...): ');
+iwrite = input(' Enter 1 to write files (0 otherwise) : ');
+
+stmod = sprintf('m%2.2i',imod);
+
+% directory containing the Hessian file that was output from the cluster
+sthess = ['/home/carltape/results/HESSIANS/'];
+
+% load the data and event weights (compute_misfit.m)
+idatacov = input(' Enter idatacov (1--by event, 2--by window, 3--none) : ');
+stdcovs = {'event','window','none'};
+stdcov = stdcovs{idatacov};
+ftag = [stmod '_' stdcov];
+odir0 = ['OUTPUT_SUBSPACE/' stmod '/' stdcov '/'];
+load([odir0 ftag '_data_norms']);
+nsrc = length(eids);
+
+%----------------------------------------------
+% check the slick operation used in subspace_update.f90
+
+if 0==1
+ npar = 10;
+ nsrc = 4;
+ G = rand(npar,nsrc);
+ d = rand(nsrc,1);
+
+ x = zeros(npar,1);
+ for isrc=1:nsrc % loop over kernels (hundreds)
+ for ipar=1:npar % loop over gridpoints (millions)
+ x(ipar) = x(ipar) + G(ipar,isrc) * d(isrc);
+ end
+ end
+ G*d, x
+end
+
+%----------------------------------------------
+
+% construct the matrix of data-covariance normalization
+Cd_fac = zeros(nsrc,nsrc);
+for i = 1:nsrc
+ for j = 1:nsrc
+ Cd_fac(i,j) = dcov_fac_e(i) * dcov_fac_e(j);
+ i_ind(i,j) = i;
+ j_ind(i,j) = j;
+ end
+end
+
+% % OBSOLETE: construct the matrix of weight factors and data-covariance normalization
+% if 0==1
+% Ws = zeros(nsrc,nsrc);
+% Cd_fac = zeros(nsrc,nsrc);
+% for i = 1:nsrc
+% for j = 1:nsrc
+% Ws(i,j) = ws(i) * ws(j);
+% Cd_fac(i,j) = dcov_fac_e(i) * dcov_fac_e(j);
+% i_ind(i,j) = i;
+% j_ind(i,j) = j;
+% end
+% end
+% figure; pcolor(i_ind,j_ind,Ws); shading flat;
+% xlabel('Row index'); ylabel('Column index');
+% title('Symmetric matrix of data weight terms');
+% %caxis([0 200]);
+% colorbar; axis equal; axis tight;
+% end
+
+% load the event IDs corresponding to the kernels
+%eids = load([sthess 'kernels_done_mod']);
+%neid = length(eids);
+
+kmin = 6;
+kmax = 6;
+nker = kmax - kmin + 1;
+stkers = {'mu_kernel','kappa_kernel','both','mu_kernel_smooth','kappa_kernel_smooth','both_smooth'};
+stkerlabs = {'unsmoothed beta kernel','unsmoothed beta kernel','beta + bulk','smoothed beta kernel','smoothed bulk kernel','smoothed beta + bulk'};
+
+%stwts = {'with data weights'};
+%stwts = {'with data weights','with NO data weights'};
+nker = length(stkers);
+%stcs = {'b.','r.'};
+%iw = 1;
+
+%for iw = 1:length(stwts)
+
+for iker = kmin:kmax
+ stker = stkers{iker};
+ stkerlab = stkerlabs{iker};
+ stlab = [ftag '_' stker];
+
+ dir1 = [odir0 stker '/'];
+
+ if ~any(iker == [3 6])
+ % get the Hessian from the subspace_hessian.f90 output
+ hfile = sprintf([sthess 'hessian_index_' stmod '_' stker '_%3.3i'],nsrc);
+ H0 = hfile2hess(hfile,nsrc);
+ H = H0;
+
+ % include the weights determined from the data
+ %if iw == 1
+ % H = H .* Ws;
+ %end
+
+ % When we computed the kernels, we did not include the
+ % normalization factor based on the number of windows picked
+ % that is included within the data covariance (and used in computing Ws).
+ % Thus we compute the MATRIX Cd_fac above to account for this.
+ H = H ./ Cd_fac;
+
+ % put the matrix onto order-1 values for computation purposes
+ %Hnorm = norm(H); H = H / norm(H);
+
+ else
+ hfile_beta = sprintf([sthess 'hessian_index_' stmod '_' stkers{iker-2} '_%3.3i'],nsrc);
+ H0_beta = hfile2hess(hfile_beta,nsrc);
+ hfile_bulk = sprintf([sthess 'hessian_index_' stmod '_' stkers{iker-1} '_%3.3i'],nsrc);
+ H0_bulk = hfile2hess(hfile_bulk,nsrc);
+
+ %H_beta = H0_beta .* Ws ./ Cd_fac;
+ %H_bulk = H0_bulk .* Ws ./ Cd_fac;
+ H_beta = H0_beta ./ Cd_fac;
+ H_bulk = H0_bulk ./ Cd_fac;
+ H = H_beta + H_bulk;
+
+ dtemp = [[1:nsrc]' eids_num diag(H_beta) diag(H_bulk) diag(H) diag(H_bulk)./diag(H_beta) Ns];
+ dtemp_sort1 = sortrows(dtemp,-6); % sort by relative importance of bulk-to-shear
+ dtemp_sort2 = sortrows(dtemp,-5); % sort by diagonal
+
+ % sort all entries of the Hessian -- note that some are negative
+ % we take the upper traingular elements, excluding the diagonal
+ Htri = triu(H,1); itri = triu(i_ind,1); jtri = triu(j_ind,1);
+ Hcols = [Htri(:) itri(:) jtri(:) ];
+ [Hsort,iHsort_good] = sortrows(Hcols, -1);
+ [Hsort,iHsort_bad] = sortrows(Hcols, 1);
+
+ if iwrite == 1
+ filename = [dir1 'hessian_bulk-shear_' stlab];
+ fid = fopen(filename,'w');
+ fprintf(fid,'Hessian diagonal contributions from beta and bulk:\n');
+ fprintf(fid,'%6s%10s%10s%10s%10s%10s%6s\n','index','eid','beta','bulk','total','bulk/beta','Nwin');
+ for ix = 1:nsrc
+ fprintf(fid,'%6i%10i%10.2e%10.2e%10.2e%10.4f%6i\n',dtemp_sort1(ix,:));
+ end
+ fclose(fid);
+
+ filename = [dir1 'hessian_overall_' stlab];
+ fid = fopen(filename,'w');
+ fprintf(fid,'Hessian diagonal contributions from beta and bulk:\n');
+ fprintf(fid,'%6s%10s%10s%10s%10s%10s%6s\n','index','eid','beta','bulk','total','bulk/beta','Nwin');
+ for ix = 1:nsrc
+ fprintf(fid,'%6i%10i%10.2e%10.2e%10.2e%10.4f%6i\n',dtemp_sort2(ix,:));
+ end
+ fclose(fid);
+
+ % make a list of the largest N positive and negative Hessian elements
+ nprint = 100; % must be less than nsrc*nsrc
+ filename = [dir1 'hessian_elements_good_' stlab];
+ fid = fopen(filename,'w');
+ fprintf(fid,'Largest positive non-diagonal elements of Hessian:\n');
+ fprintf(fid,'%10s%6s%6s%12s%12s\n','Hij','i','j','eid-i','eid-j');
+ for ix = 1:nprint
+ iH = iHsort_good(ix);
+ fprintf(fid,'%10.2e%6i%6i%12s%12s\n',Hcols(iH,:),eids{Hcols(iH,2)},eids{Hcols(iH,3)});
+ end
+ fclose(fid);
+
+ filename = [dir1 'hessian_elements_bad_' stlab];
+ fid = fopen(filename,'w');
+ fprintf(fid,'Largest negative non-diagonal elements of Hessian:\n');
+ fprintf(fid,'%10s%6s%6s%12s%12s\n','Hij','i','j','eid-i','eid-j');
+ for ix = 1:nprint
+ iH = iHsort_bad(ix);
+ fprintf(fid,'%10.2e%6i%6i%12s%12s\n',Hcols(iH,:),eids{Hcols(iH,2)},eids{Hcols(iH,3)});
+ end
+ fclose(fid);
+ end
+
+ disp(' gradient balance for the Hessian (subspace) inversion (mean of last column) :');
+ disp(mean( diag(H_bulk)./diag(H_beta) ))
+ end
+
+ %H = H + eye(nsrc);
+ % mutemp = inv(H + eye(nsrc)) * dnorm;
+
+ dH = diag(H);
+ trH = sum(dH);
+ trHH = sum(diag(H'*H));
+
+ %------------------------
+
+ pinds = [1:nsrc]';
+
+ disp(' properties of Hessian (min, median, mean(abs), max, std):');
+ stH1 = sprintf('min %.1e, median %.1e, mean(abs) %.1e, max %.1e, std %.1e, sum %.1e',...
+ min(H(:)), median(H(:)), mean(abs(H(:))), max(H(:)), std(H(:)), sum(H(:)) );
+ stH2 = sprintf('DIAGONAL : min %.1e, median %.1e, mean(abs) %.1e, max %.1e, std %.1e, sum %.1e',...
+ min(dH), median(dH), mean(abs(dH)), max(dH), std(dH), sum(dH) );
+ disp(stH1);
+ disp(stH2);
+
+ % plot H
+ Hfac = 0.5; % adjust to saturate the color scale as you wish
+ figure; pcolor(i_ind,j_ind,H); shading flat;
+ xlabel('Row index'); ylabel('Column index');
+ title({['Hessian (symmetric matrix) for ' stkerlab],stH1,stH2});
+ caxis([-1 1]*max(H(:))*Hfac), colormap('jet'), colorbar;
+ %map = colormap('seis'); colormap(flipud(map));
+ axis equal; axis tight;
+ fontsize(10), orient tall, wysiwyg
+
+ % sort the diagonal values of H
+ Hdiag = diag(H);
+ %Hdiag = diag(H0);
+ [jk,isort] = sortrows([pinds eids_num Hdiag],-3);
+ npick = min([20 nsrc]);
+ npick = nsrc;
+ disp(' '); disp([' first ' num2str(npick) ' sorted diagonal entries of H: ' stkerlab]);
+ for ii = 1:npick
+ jj = isort(ii);
+ disp(sprintf('%3i %10i (%3i), Ne = %5i, Hkk = %10.4e',...
+ ii,eids_num(jj),pinds(jj),Ns(jj),Hdiag(jj)));
+ end
+
+ [jtemp,isort] = sort(Hdiag,'descend');
+
+ figure(20+iker); hold on;
+ plot(pinds,log10(sort(Hdiag,'descend')),'ro');
+ %plot(pinds,log10(sort(diag(H'*H),'descend')),'b.');
+ plot([1 nsrc],log10(trH)*[1 1],'r--');
+ %plot([1 nsrc],log10(trHH)*[1 1],'b--');
+ %legend('diagonal element of H','diagonal element of HtH','trace (H)','trace (HtH)');
+ legend('diagonal element of H','trace (H)');
+ xlabel('sorted event index'); ylabel(' log10 H-kk');
+ title(sprintf(' trace of H = %6.2e --- trace of HtH = %6.2e',trH,trHH));
+ xlim([1 nsrc]); grid on;
+ fontsize(11), orient tall, wysiwyg
+
+ %----------------------------------------------
+
+ itsvd = 0;
+
+ if itsvd == 0
+
+ % regularization choices
+ numlam = 100;
+ if 0==1
+ minlampwr = -2.25; maxlampwr = 1;
+ lampwr = linspace(minlampwr,maxlampwr,numlam);
+ lamvec = 10.^lampwr * sum(dH); % based on trace of H
+ %lamvec = 10.^lampwr * sqrt(sum(diag(H'*H))); % based on trace of H'*H
+ else
+ %minlampwr = -2; maxlampwr = 5; % m03
+ minlampwr = -3; maxlampwr = 3; % m04 - m11
+ minlampwr = -5; maxlampwr = 0; % m12
+ lampwr = linspace(minlampwr,maxlampwr,numlam);
+ lamvec = 10.^lampwr;
+ end
+
+ [f_h, rss, mss, Gvec, Fvec, dof, kap, iL, iGCV, iOCV] = ridge_carl(dnorm,H,lamvec);
+ title({['model ' stmod],stkerlab,sprintf('iL = %i, iOCV = %i, iGCV = %i',iL,iOCV,iGCV)});
+
+ mu_all = zeros(nsrc,numlam);
+ mu_norm = zeros(numlam,1);
+ mu_res = zeros(numlam,1);
+ for ip = 1:numlam
+ lam = lamvec(ip);
+ %if ip==iGCV, disp('iGCV -- GCV lambda pick'); end
+ %if ip==iOCV, disp('iOCV -- OCV lambda pick'); end
+ %if ip==iL, disp('iL -- L-curve lambda pick'); end
+ disp(sprintf('%i out of %i, lam = %.2e',ip,numlam,lam));
+
+ mu_all(:,ip) = inv(H + lam*eye(nsrc))*dnorm;
+ %mu_all(:,ip) = inv(H'*H + lam^2*eye(nsrc))*H'*dnorm;
+
+ % f = inv(G'*diag(Wvec)*G + lam0^2*eye(ngrid*ndim))*G'*diag(Wvec)*d;
+ mu_norm(ip) = norm(mu_all(:,ip));
+ mu_res(ip) = norm( dnorm - H*mu_all(:,ip) );
+ end
+
+ ipick = iOCV; % KEY: select a damping parameter (iL, iOCV, iGCV)
+ ltemp = lamvec(ipick)
+ mtemp = mu_all(:,ipick);
+ mtemp_norm = mu_norm(ipick)
+
+ figure; nr=2; nc=2;
+
+ subplot(nr,nc,1); hold on; grid on;
+ plot( log10(lamvec), log10(mu_norm),'b.','markersize',10);
+ plot( log10(lamvec(ipick)), log10(mu_norm(ipick)),'ro','markersize',10);
+ xlabel('log10( lambda )'); ylabel('log10( norm of mu )');
+
+ subplot(nr,nc,2); hold on; grid on;
+ plot( log10( lamvec ), mu_res, '.');
+ plot( log10( lamvec(ipick) ), mu_res(ipick),'ro','markersize',10);
+ xlabel('log10( lambda )'); ylabel('dnorm residual');
+
+ subplot(nr,nc,3); plot( dnorm, H*mtemp, '.');
+ xlabel('dnorm obs'); ylabel('dnorm pred'); grid on; axis equal;
+
+ subplot(nr,nc,4); stdres = std(dnorm - H*mtemp);
+ plot_histo(dnorm - H*mtemp,[-4:0.5:4]*stdres);
+ %plot( 100*(dnorm - H*mtemp)./dnorm ,'.');
+ %ylabel('dnorm residual = 100(pred - obs)/obs');
+ title(' residuals : d - H*mu');
+ xlabel('event index'); grid on;
+
+ orient tall, wysiwyg, fontsize(10)
+
+% iCV = input(' Enter 1 for OCV or 0 for GCV : ');
+% if iCV == 1
+% lam = lamvec(iOCV);
+% elseif iCV == 0
+% lam = lamvec(iGCV);
+% else
+% error(' iCV must be 1 or 0');
+% end
+% mu = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm; % KEY vector
+% stp = sprintf('lambda = %.3f',lam);
+
+ % write mu vectors to file
+ if iwrite == 1
+ otag = ['mu_' stlab];
+ odir = [dir1 'mu_all/'];
+ for ip = 1:numlam
+ lam = lamvec(ip);
+ mutemp = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm;
+ filename = sprintf([otag '_p%3.3i'],ip);
+ fid = fopen([odir filename],'w');
+ for ii = 1:nsrc
+ fprintf(fid,'%24.12e\n',mutemp(ii));
+ end
+ fclose(fid);
+ end
+
+ fid = fopen([odir '/lambda_' stdcov],'w');
+ for ip = 1:numlam
+ lam = lamvec(ip);
+ fprintf(fid,'%24.12e%24.12e%6i\n',lamvec(ip),1/lamvec(ip),ip);
+ end
+ fclose(fid);
+ end
+
+ else
+
+ % analyze the singular values of H
+ [U,S,V] = svd(H);
+ s = diag(S);
+
+ pmin = min([13 nsrc]);
+ ifit = [(pmin+1):60];
+ %ifit = pinds;
+ [xf,yf,mf,bf,rms] = linefit(pinds(ifit), log10(s(ifit)));
+
+ %subplot(nr,nc,iker);
+ figure; hold on;
+ plot(xf,yf,'g','linewidth',6);
+ plot(pinds,log10(s),'b.','markersize',10);
+ plot(pmin*[1 1],[min(log10(s)) max(log10(s))],'k','linewidth',2)
+ grid on; xlabel(xlab1); ylabel(ylab1);
+ title(['singular value spectra for Hessian : ' stker]);
+ fontsize(10), orient tall, wysiwyg
+
+ %-------------------------
+
+ %tpinds = [1:10:100];
+ tpinds = pinds;
+ [mu_all,rss,f_r_ss] = tsvd(dnorm,H,tpinds);
+
+ % norms of mu vectors
+ mu_norm = zeros(nsrc,1);
+ for ip = 1:nsrc
+ mu_norm(ip) = norm(mu_all(:,ip));
+ end
+
+ figure;
+ subplot(2,2,1); hold on; grid on;
+ plot(pinds,log10(s),'b.','markersize',10);
+ xlabel(xlab1); ylabel(ylab1);
+ subplot(2,2,2); hold on; grid on;
+ plot( tpinds, log10(rss),'b.','markersize',10)
+ xlabel(xlab2); ylabel(ylab2);
+ subplot(2,2,3); hold on; grid on;
+ plot( tpinds, log10(mu_norm),'b.','markersize',10)
+ xlabel(xlab2); ylabel('log10( norm of mu )');
+ fontsize(10), orient tall, wysiwyg
+
+ % write mu vectors to file
+ if iwrite == 1
+ odir = [dir1 'mu_all/'];
+ for ip = 1:length(tpinds);
+ pmax = tpinds(ip);
+ filename = sprintf([odir '_p%3.3i'],pmax);
+ fid = fopen([odir '/' filename],'w');
+ for ii = 1:nsrc
+ fprintf(fid,'%24.12e\n',mu_all(ii,ip));
+ end
+ fclose(fid);
+ end
+ end
+ end
+
+
+end
+
+%end
+
+% write misfit values to file -- these do not depend on the choice of model
+% variable used in constructing the Hessian
+if iwrite == 1
+ % write data-norm vector to file
+ fid = fopen([odir0 ftag '_dmisfit'],'w');
+ fprintf(fid,'%24.12e\n',dmisfit);
+ fclose(fid);
+
+ % write data-norm vector to file
+ fid = fopen([odir0 ftag '_nwin_tot'],'w');
+ fprintf(fid,'%.0f\n',N);
+ fclose(fid);
+
+ % write data-norm vector to file
+ fid = fopen([odir0 ftag '_data_norm'],'w');
+ for isrc = 1:nsrc
+ fprintf(fid,'%24.12e\n',dnorm(isrc));
+ end
+ fclose(fid);
+
+ % write the factors for the data covariance (nevent x 1)
+ % WRITE THIS IN DOUBLE PRECISION, JUST FOR SIMPLICITY
+ fid = fopen([odir0 ftag '_dcov_fac'],'w');
+ for isrc = 1:nsrc
+ fprintf(fid,'%24.12e\n',dcov_fac_e(isrc));
+ end
+ fclose(fid);
+end
+
+%======================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/suf.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/suf.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/suf.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,39 +0,0 @@
-function suffix=suf(string,delim,iter)
-% suffix=SUF(string,delim,iter)
-%
-% Finds the suffix or extension in a string, say a filename.
-%
-% INPUT:
-%
-% string The string to be parsed
-% delim The string that separates prefix from suffix
-% iter 0 Don't iterate
-% 1 Iterate on the suffix to get the last bit
-%
-% See also PREF
-%
-% Last modified by fjsimons-at-alum.mit.edu, 05/06/2008
-
-defval('delim','.');
-defval('iter',1)
-
-ldel=length(delim);
-[prefix,suffix]=strtok(string,delim);
-
-% This added 17.11.2004
-if iscell(suffix)
- suffix=cell2mat(suffix);
-end
-
-suffix=suffix(ldel+1:end);
-
-if iter==1
- % This suffix might not be the last one
- while findstr(suffix,delim)
- suffix=suf(suffix,delim);
- end
-end
-
-
-
-
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/wysiwyg.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/wysiwyg.m 2009-07-29 20:16:08 UTC (rev 15485)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/wysiwyg.m 2009-07-30 14:23:31 UTC (rev 15486)
@@ -1,17 +0,0 @@
-function wysiwyg
-%WYSIWYG -- this function is called with no args and merely
-% changes the size of the figure on the screen to equal
-% the size of the figure that would be printed,
-% according to the papersize attribute. Use this function
-% to give a more accurate picture of what will be
-% printed.
-% Dan(K) Braithwaite, Dept. of Hydrology U.of.A 11/93
-
-unis = get(gcf,'units');
-ppos = get(gcf,'paperposition');
-set(gcf,'units',get(gcf,'paperunits'));
-pos = get(gcf,'position');
-pos(3:4) = ppos(3:4);
-set(gcf,'position',pos);
-set(gcf,'units',unis);
-
More information about the CIG-COMMITS
mailing list