[cig-commits] r11891 - seismo/3D/mt_measure_adj/scripts_tomo/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Thu May 1 17:17:00 PDT 2008
Author: carltape
Date: 2008-05-01 17:17:00 -0700 (Thu, 01 May 2008)
New Revision: 11891
Added:
seismo/3D/mt_measure_adj/scripts_tomo/matlab/hfile2hess.m
Modified:
seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m
seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl
seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m
Log:
Updates to matlab codes for subspace method.
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m 2008-05-01 23:54:26 UTC (rev 11890)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/compute_misfit.m 2008-05-02 00:17:00 UTC (rev 11891)
@@ -41,7 +41,7 @@
% -- N is the total number of windows: N = sum_e N_e
stdcov = {'event','window','none'};
-iwrite = 1;
+iwrite = input(' Enter 1 to write files (0 otherwise) : ');
%-------------------------
% read in list of event IDs and sources
@@ -72,6 +72,10 @@
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
@@ -118,13 +122,23 @@
mfile = [eids{ii} '_T' sTper '_window_chi'];
filename = [dir1 mfile];
- if exist(filename,'file')
- [strec0,cmp,iwin,iker,...
+ if ~exist(filename,'file')
+ disp('--> file does not exist (or nwin = 0)');
+ else
+ [stnet0,strec0,comp,iwin,iker,...
chi1,chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
meas1,measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
sigma1,sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,tr_chi,am_chi] = read_window_chi(filename);
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
@@ -135,8 +149,8 @@
% assign the string component an integer index
index_comp = NaN * ones(nwin,1);
for icomp = 1:ncomp
- ind_cmp = find( strcmp(comps(icomp), cmp) == 1 );
- index_comp(ind_cmp) = icomp;
+ ind_comp = find( strcmp(comps(icomp), comp) == 1 );
+ index_comp(ind_comp) = icomp;
end
% measurement indices
@@ -144,21 +158,22 @@
k2 = k1+nwin-1;
kinds = [k1:k2]';
- meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_rec index_comp iwin ...
+ meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin ...
iker chi1 measCC_dT sigmaCC_dT measMT_dT sigmaMT_dT tr_chi];
% 1 kinds
% 2 index_T
% 3 index_event
- % 4 index_rec
- % 5 index_comp
- % 6 iwin
- % 7 iker
- % 8 chi1
- % 9 measCC_dT
- % 10 sigmaCC_dT
- % 11 measMT_dT
- % 12 sigmaMT_dT
- % 13 tr_chi
+ % 4 index_network
+ % 5 index_rec
+ % 6 index_comp
+ % 7 iwin
+ % 8 iker
+ % 9 chi1
+ % 10 measCC_dT
+ % 11 sigmaCC_dT
+ % 12 measMT_dT
+ % 13 sigmaMT_dT
+ % 14 tr_chi
end
end
@@ -184,25 +199,75 @@
meas_array(itest,:)
display_meas(meas_array(itest,:),Tvec,eids,strec,comps);
meas_array(itest,13)
- 0.5 * meas_array(itest,9)^2 / meas_array(itest,10)^2
- 0.5 * meas_array(itest,11)^2 / meas_array(itest,10)^2
+ 0.5 * meas_array(itest,10)^2 / meas_array(itest,11)^2
+ 0.5 * meas_array(itest,12)^2 / meas_array(itest,11)^2
end
end
-% tally windows by event
-nwin_all = zeros(nevent,2);
-for ii = 1:nevent
- for tt = 1:nper
+%----------------------------------------------
+% tally windows by:
+% (1) event
+% (2) network
+% (3) station
+% (4) component
+nwin_all_event = zeros(nevent,2);
+for tt = 1:nper
+ for ii = 1:nevent
imatch = find( and( meas_array(:,3)==ii, meas_array(:,2)==tt) );
- nwin_all(ii,tt) = length(imatch);
+ nwin_all_event(ii,tt) = length(imatch);
end
end
-nwin_all_event = sum(nwin_all, 2);
-nwin_all_per = sum(nwin_all, 1);
-nwin_tot = [nwin_all nwin_all_event ; nwin_all_per sum(sum(nwin_all)) ];
-dtemp = [ [eids_num ; 0] nwin_tot];
-%sortrows(dtemp,3)
+nwin_all_net = zeros(nnet,2);
+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);
+ end
+end
+nwin_all_rec = zeros(nrec,2);
+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);
+ end
+end
+nwin_all_comp = zeros(ncomp,2);
+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);
+ end
+end
+if iwrite == 1
+ stks = {'event','network','receiver','component'};
+ for kk = 1:length(stks)
+ switch kk
+ case 1, nwin_out = nwin_all_event; labs = eids;
+ case 2, nwin_out = nwin_all_net; labs = stnet;
+ case 3, nwin_out = nwin_all_rec; labs = strec;
+ case 4, nwin_out = nwin_all_comp; labs = comps;
+ 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],-4);
+
+ % write out the file
+ fid = fopen(['window_picks_by_' stks{kk}],'w');
+ fprintf(fid,'%12s%10s%10s%10s\n',stks{kk},...
+ ['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
+ fprintf(fid,'%12s%10i%10i%10i\n','TOTAL',ntot_1,sum(ntot_1));
+ for ii = 1:nw
+ jj = isort(ii);
+ fprintf(fid,'%12s%10i%10i%10i\n',labs{jj},nwin_out(jj,:),ntot_2(jj));
+ end
+ fclose(fid);
+ end
+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);
@@ -222,19 +287,16 @@
end
% check number of windows
-if sum(Ns) ~= sum(sum(nwin_all)), error('inconsistent values for total windows'); end
+if sum(Ns) ~= length(meas_array), error('inconsistent values for total windows'); end
% compute data vector for subspace method
% factor of 2 cancels the 0.5 used in computing the misfit (mt_measure_adj.f90)
dnorm_sq = zeros(nevent,1);
for ii = 1:nevent
imatch = find( meas_array(:,3)==ii );
- dnorm_sq(ii) = sum(2 * meas_array(imatch,13) ./ dcov_fac(imatch) );
+ dnorm_sq(ii) = sum(2 * meas_array(imatch,14) ./ dcov_fac(imatch) );
end
-% check number of windows
-if sum(Ns) ~= sum(sum(nwin_all)), error('inconsistent values for total windows'); end
-
% weight vector for the subspace method
dnorm = sqrt(dnorm_sq);
ws = 1 ./ dnorm;
@@ -270,7 +332,7 @@
% % find records that meet particular criteria
% DT_MIN = 5;
% DT_SIGMA_MIN = 1;
-% icheck = find( and( meas_array(:,9) >= DT_MIN, meas_array(:,10) <= DT_SIGMA_MIN) )
+% icheck = find( and( meas_array(:,10) >= DT_MIN, meas_array(:,11) <= DT_SIGMA_MIN) )
% [junk, isort] = sortrows( meas_array(icheck,:), -9 )
% meas_disp = meas_array(icheck(isort),:);
% display_meas(meas_disp,Tvec,eids,strec,comps);
@@ -278,7 +340,7 @@
% find records that meet particular criteria
CHI_MIN = 100;
DT_SIGMA_MIN = 0.19;
-icheck = find( and( meas_array(:,13) >= CHI_MIN, meas_array(:,10) < DT_SIGMA_MIN) );
+icheck = find( and( meas_array(:,14) >= CHI_MIN, meas_array(:,11) < DT_SIGMA_MIN) );
[junk, isort] = sortrows( meas_array(icheck,:), -13 );
meas_disp = meas_array(icheck(isort),:);
display_meas(meas_disp,Tvec,eids,strec,comps);
Added: seismo/3D/mt_measure_adj/scripts_tomo/matlab/hfile2hess.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/hfile2hess.m (rev 0)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/hfile2hess.m 2008-05-02 00:17:00 UTC (rev 11891)
@@ -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
+
+%======================================================
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl 2008-05-01 23:54:26 UTC (rev 11890)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/make_dirs.pl 2008-05-02 00:17:00 UTC (rev 11891)
@@ -15,7 +15,7 @@
@dcov_tags = ("event","window","none");
# possible labels for choice of kernels
- at kers = ("mu_kernel","kappa_kernel","mu_kernel_smooth","kappa_kernel_smooth");
+ at kers = ("mu_kernel","kappa_kernel","mu_kernel_smooth","kappa_kernel_smooth","both","both_smooth");
$ndcov = @dcov_tags;
$nker = @kers;
@@ -35,7 +35,7 @@
if(-e $dir) {
print "--> dir exists -- now deleting and remaking\n";
- `rm -rf $dir`;
+ #`rm -rf $dir`;
} else {
print "--> dir does not exist -- now making\n";
}
Modified: seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m 2008-05-01 23:54:26 UTC (rev 11890)
+++ seismo/3D/mt_measure_adj/scripts_tomo/matlab/subspace_specfem.m 2008-05-02 00:17:00 UTC (rev 11891)
@@ -8,7 +8,7 @@
%
% See also wave2d_subspace.m
%
-% calls xxx
+% calls hfile2hess.m
% called by xxx
%
@@ -27,7 +27,7 @@
stker0 = '/net/sierra/raid1/carltape/results/KERNELS/';
-iwrite = 1;
+iwrite = input(' Enter 1 to write files (0 otherwise) : ');
%----------------------------------------------
% check the slick operation used in subspace_update.f90
@@ -76,11 +76,11 @@
%eids = load([stker0 'kernels_done_mod']);
%neid = length(eids);
-kmin = 3;
-kmax = 3;
+kmin = 6;
+kmax = 6;
nker = kmax - kmin + 1;
-stkers = {'mu_kernel','mu_kernel_smooth','kappa_kernel_smooth'};
-stkerlabs = {'unsmoothed beta kernel','smoothed beta kernel','smoothed bulk kernel'};
+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'};
@@ -93,41 +93,56 @@
for iker = kmin:kmax
stker = stkers{iker};
stkerlab = stkerlabs{iker};
+ stlab = [stker '_' stdcov{idatacov}];
- hfile = sprintf([stker0 'hessian_index_' stker '_%3.3i'],nsrc);
- if ~exist(hfile,'file'), error([hfile ' does not exist']); end
- [kall,iall,jall,Hij] = textread(hfile,'%f%f%f%f');
- n = length(kall);
- nhess = (-1 + sqrt(1+8*n))/2; % quadratic formula
+ if ~any(iker == [3 6])
+ % get the Hessian from the subspace_hessian.f90 output
+ hfile = sprintf([stker0 'hessian_index_' stker '_%3.3i'],nsrc);
+ H0 = hfile2hess(hfile,nsrc);
+ H = H0;
- if nhess ~= nsrc, error('number of events must be consistent'); end
- pinds = [1:nsrc]';
+ % include the weights determined from the data
+ if iw == 1
+ H = H .* Ws;
+ end
- % construct the Hessian matrix
- H0 = zeros(nsrc,nsrc);
- for i = 1:nsrc
- for j = i:nsrc
- ih = find( and( i == iall, j == jall) );
- H0(i,j) = Hij(ih);
+ % 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([stker0 'hessian_index_' stkers{iker-2} '_%3.3i'],nsrc);
+ H0_beta = hfile2hess(hfile_beta,nsrc);
+ hfile_bulk = sprintf([stker0 'hessian_index_' 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 = H_beta + H_bulk;
+
+ dtemp = [[1:nsrc]' eids_num diag(H_beta) diag(H_bulk) diag(H) diag(H_bulk)./diag(H_beta) ];
+ dtemp_sort = sortrows(dtemp,-6);
+
+ if iwrite == 1
+ filename = ['hessian_' stlab];
+ fid = fopen(filename,'w');
+ fprintf(fid,'Hessian diagonal contributions from beta and bulk:\n');
+ fprintf(fid,'%6s%10s%10s%10s%10s%10s\n','index','eid','beta','bulk','total','bulk/beta');
+ for ix = 1:nsrc
+ fprintf(fid,'%6i%10i%10.2e%10.2e%10.2e%10.4f\n',dtemp_sort(ix,:));
+ 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
- H0 = H0 + H0' - diag(diag(H0)); % fill the lower triangle
- H = H0;
- % include the weights determined from the data
- if iw == 1
- H = H .* Ws;
- end
-
- % When we computed the kernels, we did not include the
- % normalization factor based on the number of windows picked
- % that is included within the data covariance (and used in computing Ws).
- % Thus, we compute the MATRIX Cd_fac above to account for this.
- H = H ./ Cd_fac;
-
- % put the matrix onto order-1 values for computation purposes
- %Hnorm = norm(H); H = H / norm(H);
-
%H = H + eye(nsrc);
% mutemp = inv(H + eye(nsrc)) * dnorm;
@@ -135,6 +150,10 @@
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(:)) );
@@ -166,8 +185,10 @@
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(diag(H),'descend')),'ro');
+ 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--');
@@ -179,7 +200,7 @@
fontsize(11), orient tall, wysiwyg
%----------------------------------------------
-
+
itsvd = 0;
if itsvd == 0
@@ -192,15 +213,16 @@
lamvec = 10.^lampwr * sum(dH); % based on trace of H
%lamvec = 10.^lampwr * sqrt(sum(diag(H'*H))); % based on trace of H'*H
else
- minlampwr = -23; maxlampwr = -3; % 1
- minlampwr = -12; maxlampwr = -5; % 3
+ minlampwr = -16; maxlampwr = -9; % 1
+ %minlampwr = -12; maxlampwr = -5; % 3
lampwr = linspace(minlampwr,maxlampwr,numlam);
+ lampwr = lampwr + 16 - 2;
lamvec = 10.^lampwr;
end
[f_h, rss, mss, Gvec, Fvec, dof, kap, iL, iGCV, iOCV] = ridge_carl(dnorm,H,lamvec);
- title(stkerlab);
-
+ title({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);
@@ -238,8 +260,12 @@
subplot(nr,nc,3); plot( dnorm, H*mtemp, '.');
xlabel('dnorm obs'); ylabel('dnorm pred'); grid on; axis equal;
- subplot(nr,nc,4); plot( 100*(dnorm - H*mtemp)./dnorm ,'.'); grid on;
- xlabel('event index'); ylabel('dnorm residual = 100(pred - obs)/obs');
+ 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;
fontsize(10), orient tall, wysiwyg
@@ -256,11 +282,12 @@
% write mu vectors to file
if iwrite == 1
- odir = ['mu_dirs/mu_' stkers{iker} '_' stdcov{idatacov} ];
+ otag = ['mu_' stlab];
+ odir = ['mu_dirs/' otag];
for ip = 1:numlam
lam = lamvec(ip);
mutemp = inv(H'*H + lam^2*eye(nsrc,nsrc))*H'*dnorm;
- filename = sprintf([odir '_p%3.3i'],ip);
+ filename = sprintf([otag '_p%3.3i'],ip);
fid = fopen([odir '/' filename],'w');
for ii = 1:nsrc
fprintf(fid,'%24.12e\n',mutemp(ii));
@@ -322,7 +349,7 @@
% write mu vectors to file
if iwrite == 1
- odir = ['mu_dirs/mu_' stkers{iker} '_' stdcov{idatacov} ];
+ odir = ['mu_dirs/mu_' stlab ];
for ip = 1:length(tpinds);
pmax = tpinds(ip);
filename = sprintf([odir '_p%3.3i'],pmax);
@@ -348,7 +375,7 @@
% write data-norm vector to file
fid = fopen('nwin_tot','w');
- fprintf(fid,'%i\n',N);
+ fprintf(fid,'%.0f\n',N);
fclose(fid);
% write data-norm vector to file
@@ -359,9 +386,10 @@
fclose(fid);
% write the factors for the data covariance (nevent x 1)
+ % WRITE THIS IN DOUBLE PRECISION, JUST FOR SIMPLICITY
fid = fopen(['dcov_fac_' stdcov{idatacov}],'w');
for isrc = 1:nsrc
- fprintf(fid,'%12i\n',dcov_fac_e(isrc));
+ fprintf(fid,'%24.12e\n',dcov_fac_e(isrc));
end
fclose(fid);
end
More information about the cig-commits
mailing list