[cig-commits] r18230 - seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts
carltape at geodynamics.org
carltape at geodynamics.org
Wed Apr 13 09:44:53 PDT 2011
Author: carltape
Date: 2011-04-13 09:44:53 -0700 (Wed, 13 Apr 2011)
New Revision: 18230
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_VR.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_compute_hessian.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_dsub.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m2mvec.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_plot_hessian.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_tsvd_hessian.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_hessian.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/diff_files.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/quad_min_4.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cell2gll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gll2cell.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m_gll2cell.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_splitm.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_src.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_str.m
Log:
added and updated matlab testing scripts for iterative 2D inversion
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/diff_files.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/diff_files.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/diff_files.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -34,10 +34,9 @@
if 0==1
bdir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_OUTPUT/';
- file1 = [bdir 'run_1500/READ_IN_CGF90_CDIAG/model_m0001/src_syn_m0001.dat'];
- file2 = [bdir 'run_1500/READ_IN/model_m0001/src_syn_m0001.dat'];
- diff_files(file1,file2)
-
+ file1 = [bdir 'run_1500/READ_IN_CGF90_CDIAG/model_m0001/src_syn_m0001.dat'];
+ file2 = [bdir 'run_1500/READ_IN/model_m0001/src_syn_m0001.dat'];
+ diff_files(file1,file2)
end
%=========================================================
\ No newline at end of file
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/quad_min_4.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/quad_min_4.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/quad_min_4.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -46,7 +46,7 @@
stlabs = [stlabs(1:2) fliplr(stlabs(3:5))]
end
- axpoly = axes_expand([x1plot x2plot 0 max([y1 y2])],1.2);
+ axpoly = axes_expand([x1plot x2plot 0 max([y1 y2])],1.2,1);
axpoly(3) = 0;
dy = axpoly(4) - axpoly(3);
dx = axpoly(2) - axpoly(1);
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_VR.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_VR.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_VR.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,60 @@
+%
+% function wave2d_VR
+% Carl Tape, 25-Jan-2010
+%
+%
+%
+% calls xxx
+% called by wave2d_cg_run.m
+%
+
+function [VR_i, VR_s, VR_tot] = wave2d_VR(DT1,DT2,sigmaDT,Cd,Dinds)
+
+% D_inds is nsrc x 2
+nsrc = length(Dinds);
+n = length(DT1);
+
+chi1 = DT1.^2 ./ Cd;
+chi2 = DT2.^2 ./ Cd;
+
+VR_i = log( chi1 ./ chi2 );
+
+% find indices of measurements for which both DT1 and DT2 are less than sigmaDT
+% --> VR = 0 for these
+ivr = find( and( abs(DT1) < sigmaDT, abs(DT2) < sigmaDT ) );
+disp(sprintf('%i / %i measurements with DT < %.2f',length(ivr),n,sigmaDT));
+%for ii=1:length(ivr)
+% disp(sprintf('%10i%10.4f%10.4f',ivr(ii),DT1(ivr(ii)),DT2(ivr(ii))))
+%end
+VR_i(ivr) = 0;
+
+VR_s = zeros(nsrc,1);
+for ii=1:nsrc
+ inds = [Dinds(ii,1) : Dinds(ii,2)];
+ VR_s(ii) = sum( VR_i(inds) );
+end
+
+VR_tot = log( sum(chi1) / sum(chi2) );
+
+if 1==1
+ figure; nr=3; nc=1;
+
+ subplot(nr,nc,1); hold on;
+ plot( DT1, DT2, '.');
+ xlabel('DT - model 1'); ylabel('DT - model 2');
+ axis equal, axis tight, grid on;
+ ax1 = axis; val = min(ax1([2 4]));
+ plot([-1 1]*val,[-1 1]*val,'r--','linewidth',3);
+
+ subplot(nr,nc,2); hold on;
+ plot( DT2 - DT1, '.'); grid on; xlim([0 n])
+ xlabel('index'); ylabel('DT2 - DT1')
+
+ subplot(nr,nc,3); hold on;
+ plot( VR_i, '.'); grid on; xlim([0 n])
+ xlabel('index'); ylabel('VR')
+
+ orient tall, wysiwyg
+end
+
+%=========================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cell2gll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cell2gll.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cell2gll.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -11,14 +11,14 @@
% called by xxx
%
-function fg = wave2d_cell2gll(xg,yg,xc,yc,fc,nxc,nyc)
-%function fg = wave2d_cell2gll(xg,yg,X,Y,F)
+%function fg = wave2d_cell2gll(xg,yg,xc,yc,fc,nxc,nyc)
+function fg = wave2d_cell2gll(xg,yg,X,Y,F)
% reshape xc, yc, fc into matrices
% NOTE: this assumes that they were ORIGINALLY made using meshgrid
-X = reshape(xc,nyc,nxc);
-Y = reshape(yc,nyc,nxc);
-F = reshape(fc,nyc,nxc);
+%X = reshape(xc,nyc,nxc);
+%Y = reshape(yc,nyc,nxc);
+%F = reshape(fc,nyc,nxc);
fg = interp2(X,Y,F,xg,yg,'cubic');
%fg = interp2(X,Y,F,xg,yg,'cubic');
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_cg.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -24,7 +24,10 @@
p0 = zeros(nmod,1);
else
% load p0 and g0, the PREVIOUS gradient vectors (pk and gk)
- if ~exist(gfile,'file'), error('gfile does not exist'); end
+ if ~exist(gfile,'file')
+ disp(gfile);
+ error('gfile does not exist');
+ end
[g0,p0] = textread(gfile,'%f%f');
beta_val = ( (gk - g0)' * C * gk ) / ( g0' * C * g0 );
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_compute_hessian.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_compute_hessian.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_compute_hessian.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,72 @@
+%
+% function
+% Carl Tape, 25-Jan-2010
+%
+% Compute the Hessian for the source subspace projection method.
+%
+% INPUT
+% Gstr gradient w.r.t. structure parameters
+% Gsrc gradient w.r.t. source parameters
+% Cstr covariance matrix or diagonal part (as a vector)
+% csrc diagonal covariance matrix (as a vector)
+%
+% calls xxx
+% called by xxx
+%
+
+function [Hstr,Hsrc] = wave2d_compute_hessian(Gstr,Gsrc,Cstr,csrc)
+
+[nsrc,nstr] = size(Gstr);
+[a,b] = size(Cstr);
+if or(a==1,b==1)
+ icdiag = 1;
+ Cstr = Cstr(:)'; % ensure that C is a row vector
+else
+ icdiag = 0;
+end
+
+% Hessian for structure parameters (from event kernels)
+Hstr = zeros(nsrc,nsrc);
+if icdiag == 1
+ disp('constructing the Hessian with diagonal Cm for structural part');
+ for ii = 1:nsrc
+ for jj = 1:nsrc
+ Hstr(ii,jj) = dot(Gstr(ii,:), Cstr.*Gstr(jj,:));
+ end
+ end
+
+else
+ disp('constructing the Hessian with full Cm for structural part');
+ Hstr = Gstr * Cstr * Gstr';
+end
+
+% Hessian for source parameters
+Hsrc = zeros(nsrc,nsrc);
+Hsrc = Gsrc * diag(csrc) * Gsrc'; % csrc is Msrc x 1
+
+disp('norm(Hstr), norm(Hsrc):');
+norm(Hstr), norm(Hsrc)
+
+%=========================================================
+
+if 0==1
+ S = 5; % number of sources
+ Mstr = 20;
+ Msrc = S*3;
+ inds_str = [1:Mstr]';
+ inds_src = [1:Msrc]';
+
+ % fill G
+ Gstr = rand(S,Mstr);
+ Gsrc = zeros(S,Msrc);
+ %Gsrc = repmat(diag(rand(S,1)),1,Msrc1+Msrc2);
+ Cstr = rand(Mstr,Mstr); % full matrix
+ csrc = rand(Msrc,1); % vector representing diagonal matrix
+
+ % testing
+ Hstr = wave2d_compute_hessian(Gstr,Gsrc,Cstr,csrc)
+ Hstr = wave2d_compute_hessian(Gstr,Gsrc,diag(diag(Cstr)),csrc)
+ Hstr = wave2d_compute_hessian(Gstr,Gsrc,diag(Cstr),csrc)
+end
+
+%=========================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_dsub.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_dsub.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_dsub.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,48 @@
+%
+% function [dsub, indmat] = wave2d_dsub(d,covd,nvec)
+% Carl Tape, 25-Feb-2010
+%
+% This function splits a wave2d.f90 model vector into constituent parts.
+%
+% INPUT
+% d N x 1 data vector
+% covd N x 1 diagonal terms of data covariance matrix
+% nvec S x 1 vector of number of measurments per source
+%
+% OUTPUT
+% dsub S x 1 provected data vector
+%
+% calls xxx
+% called by xxx
+%
+
+function [dsub, indmat] = wave2d_dsub(d,covd,nvec)
+
+d = d(:);
+covd = covd(:);
+n = sum(nvec);
+nsrc = length(nvec);
+
+if length(d) ~= length(covd), error('d and covd not same lengths'); end
+if length(d) ~= n, error('inconsistent number of measurements'); end
+
+% indexing for measurements
+cnvec = cumsum(nvec);
+indmat = [ [1 ; 1+cnvec(1:end-1)] cnvec ];
+disp([indmat nvec]);
+
+dsub = zeros(nsrc,1);
+for isrc = 1:nsrc
+ inds = [indmat(isrc,1) : indmat(isrc,2)];
+ dtemp = d(inds);
+ ctemp = covd(inds);
+
+ % L2-norm-squared
+ %dsub(isrc) = dtemp' * diag(1./ctemp) * dtemp;
+ dsub(isrc) = sum( dtemp.^2 ./ ctemp );
+end
+
+% check
+if any(dsub==0), error(' For at least one source, there is perfect fit.'); end
+
+%=========================================================
\ No newline at end of file
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gll2cell.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gll2cell.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_gll2cell.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -5,26 +5,29 @@
% Given an input grid (xg,yg) for which a function is defined, this
% function returns the index of the nearest (xg,yg) point to each input
% (x,y) point. This provides a quick vector for nearest neighbor
-% interpolation, which is of course crude.
+% interpolation, which is of course CRUDE.
%
% calls xxx
% called by xxx
%
-function [inds, dmins] = wave2d_gll2cell(xg,yg,xc,yc)
+%function [inds, dmins] = wave2d_gll2cell(xg,yg,xc,yc)
+function Fc = wave2d_gll2cell(xg,yg,fg,Xc,Yc)
-nc = length(xc); % number of input points (regular mesh)
-ng = length(xg); % number of GLL gridpoints (iregugular mesh)
+Fc = griddata(xg,yg,fg,Xc,Yc,'cubic');
-% loop over input points
-inds = zeros(nc,1);
-dmin = zeros(nc,1);
-for ii=1:nc
- dtemp = sqrt( (xc(ii)-xg).^2 + (yc(ii)-yg).^2 );
- [dmin,imin] = min( dtemp );
- inds(ii) = imin(1);
- dmins(ii) = dmin(1);
-end
+% nc = length(xc); % number of input points (regular mesh)
+% ng = length(xg); % number of GLL gridpoints (irregular mesh)
+%
+% % loop over input points
+% inds = zeros(nc,1);
+% dmin = zeros(nc,1);
+% for ii=1:nc
+% dtemp = sqrt( (xc(ii)-xg).^2 + (yc(ii)-yg).^2 );
+% [dmin,imin] = min( dtemp );
+% inds(ii) = imin(1);
+% dmins(ii) = dmin(1);
+% end
iplot = 0;
if iplot==1
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m2mvec.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m2mvec.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m2mvec.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,21 @@
+%
+% function mvec = wave2d_m2mvec(m,m_inds,beta0)
+% Carl Tape, 25-Jan-2010
+%
+% Convert a model vector from 'inversion parameters' to 'physical
+% parameters'. In our case, only the structural parameters differ.
+%
+% calls xxx
+% called by xxx
+%
+
+function mvec = wave2d_m2mvec(m,m_inds,beta0)
+
+nmod = length(m);
+nmod_str = m_inds(1,2);
+
+mvec = zeros(nmod,1);
+mvec(1:nmod_str) = beta0 * exp( m(1:nmod_str) );
+mvec(nmod_str+1 : nmod) = m(nmod_str+1 : nmod);
+
+%=========================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m_gll2cell.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m_gll2cell.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_m_gll2cell.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -4,30 +4,35 @@
%
%
%
-% calls xxx
+% calls wave2d_gll2cell.m, wave2d_cell2gll.m
% called by xxx
%
-function m_out = wave2d_m_gll2cell(m_in,xg,yg,xc,yc,nxc,nyc,nmod_src,iopt)
+function m_out = wave2d_m_gll2cell(m_in,xg,yg,Xc,Yc,nxc,nyc,nmod_src,iopt)
-nc = length(xc); % regular mesh of cells
-ng = length(xg); % irregular mesh GLL points
+nc = length(Xc(:)); % regular mesh of cells
+ng = length(xg); % irregular mesh GLL points
if iopt==1 % gll2cell
nmod_out = nc+nmod_src;
if length(m_in) ~= ng+nmod_src, error('check dimensions of input'); end
+ m_out = zeros(nmod_out,1);
- m_out = zeros(nmod_out,1);
- iGLL = wave2d_gll2cell(xg,yg,xc,yc);
- m_out(1:nc) = m_in(iGLL);
+ %iGLL = wave2d_gll2cell(xg,yg,xc,yc);
+ %m_out(1:nc) = m_in(iGLL);
+ m_out(1:nc) = wave2d_gll2cell(xg,yg,m_in(1:ng),Xc,Yc);
+
m_out(nc+1:nmod_out) = m_in(ng+1:ng+nmod_src);
else % cell2gll
nmod_out = ng+nmod_src;
if length(m_in) ~= nc+nmod_src, error('check dimensions of input'); end
+ m_out = zeros(nmod_out,1);
- m_out = zeros(nmod_out,1);
- m_out(1:ng) = wave2d_cell2gll(xg,yg,xc,yc,m_in(1:nc),nxc,nyc);
+ %m_out(1:ng) = wave2d_cell2gll(xg,yg,xc,yc,m_in(1:nc),nxc,nyc);
+ M_in = reshape(m_in(1:nc),nyc,nxc);
+ m_out(1:ng) = wave2d_cell2gll(xg,yg,Xc,Yc,M_in);
+
m_out(ng+1:nmod_out) = m_in(nc+1:nc+nmod_src);
end
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_norm_sq.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -28,7 +28,8 @@
npw = 1 ./ weights;
elseif imnorm==1
- if icdiag==1, Cmat = 1./C; else Cmat = inv(C); end
+ %if icdiag==1, Cmat = 1./C; else Cmat = inv(C); end
+ Cmat = diag(1./diag(C)); % APPROXIMATION
mvec = m - mprior;
npw = weights;
end
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_plot_hessian.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_plot_hessian.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_plot_hessian.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,64 @@
+%
+% function wave2d_plot_hessian(Hstr,Hsrc,H)
+% Carl Tape, 16-Feb-2010
+%
+% Plot the Hessian for the source subspace projection method.
+%
+% calls xxx
+% called by wave2d_subspace.m, wave2d_cg_run.m
+%
+
+function wave2d_plot_hessian(Hstr,Hsrc,H)
+
+disp(' properties of Hessian (min, median, mean(abs), max, std):');
+stH = sprintf('min %.2e, median %.2e, mean(abs) %.2e, max %.2e, std %.2e',...
+ min(H(:)), median(H(:)), mean(abs(H(:))), max(H(:)), std(H(:)));
+disp(stH);
+
+%if ~and(INV_STRUCT==1, INV_SOURCE==1)
+if ~and(norm(Hstr)==0, norm(Hsrc)==0)
+
+ figure;
+ %pcolor(H); shading flat;
+ imagesc(H);
+ xlabel('Source index'); ylabel('Source index');
+ title({'Hessian (symmetric matrix)',stH});
+ map = colormap('gray'); colormap(flipud(map));
+ colorbar; axis equal; axis tight;
+
+else
+ figure; nr=3; nc=1;
+ subplot(nr,nc,1);
+ %pcolor(H); shading flat;
+ imagesc(H);
+ xlabel('Source index'); ylabel('Source index');
+ title({'Hessian (symmetric matrix)',stH});
+ map = colormap('gray'); colormap(flipud(map));
+ colorbar; axis equal; axis tight;
+
+ subplot(nr,nc,2);
+ %pcolor(Hstr); shading flat;
+ imagesc(Hstr);
+ xlabel('Source index'); ylabel('Source index');
+ title('Hessian for structure (symmetric matrix)');
+ map = colormap('gray'); colormap(flipud(map));
+ colorbar; axis equal; axis tight;
+
+ subplot(nr,nc,3);
+ %pcolor(Hsrc); shading flat;
+ imagesc(Hsrc);
+ xlabel('Source index'); ylabel('Source index');
+ title('Hessian for source (diagonal matrix)');
+ map = colormap('gray'); colormap(flipud(map));
+ colorbar; axis equal; axis tight;
+ orient tall, wysiwyg
+end
+
+% spectrum of singular values
+[U,S,V] = svd(H);
+s = diag(S);
+figure; semilogy(s,'.-','markersize',18);
+xlabel('Index'); ylabel('Singular value');
+title(sprintf('singular values of H, cond(H) = %.2e',cond(H)));
+
+%=========================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_splitm.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_splitm.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_splitm.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -1,5 +1,5 @@
%
-% function
+% function [m_str, m_ts, m_xs, m_ys] = wave2d_splitm(m,m_inds)
% Carl Tape, 25-Jan-2010
%
% This function splits a wave2d.f90 model vector into constituent parts.
@@ -10,7 +10,11 @@
function [m_str, m_ts, m_xs, m_ys] = wave2d_splitm(m,m_inds)
-if length(m) ~= m_inds(4,2), error('incompatible indexing'); end
+if length(m) ~= m_inds(4,2)
+ whos m
+ m_inds
+ error('incompatible indexing');
+end
m_str = m(m_inds(1,1):m_inds(1,2));
m_ts = m(m_inds(2,1):m_inds(2,2));
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_tsvd_hessian.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_tsvd_hessian.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_tsvd_hessian.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,68 @@
+%
+% function wave2d_tsvd_hessian(H,dnorm)
+% Carl Tape, 25-Jan-2010
+%
+% This function examines possible regularization of H via truncated
+% singular value decomposition.
+%
+% calls tsvd.m
+% called by wave2d_cg_run.m
+%
+
+function [mu_all] = wave2d_tsvd_hessian(H,dnorm)
+
+nsrc = length(dnorm);
+
+% Matlab SVD
+[~,S,~] = svd(H);
+s = diag(S);
+
+pinds = [1:nsrc]';
+[mu_all,rss,f_r_ss] = tsvd(dnorm,H,pinds); % KEY: TSVD
+
+% norms of mu vectors
+mu_norm = zeros(nsrc,1);
+for ip = 1:nsrc
+ mu_norm(ip) = norm(mu_all(:,ip));
+end
+
+figure; nr=2; nc=2;
+xlab1 = 'p, singular value index';
+xlab2 = 'p, singular value truncation index';
+ylab1 = 'singular value';
+ylab2 = 'misfit : dot[ d - H*mu(p), d - H*mu(p) ]';
+
+subplot(nr,nc,1); plot(pinds,s,'.','markersize',20);
+grid on; xlabel(xlab1); ylabel(ylab1);
+subplot(nr,nc,2); semilogy(pinds,s,'.','markersize',20);
+grid on; xlabel(xlab1); ylabel(ylab1);
+subplot(nr,nc,3); plot(pinds,rss,'.','markersize',20);
+grid on; xlabel(xlab2); ylabel(ylab2);
+subplot(nr,nc,4); semilogy(pinds,rss,'.','markersize',20);
+grid on; xlabel(xlab2); ylabel(ylab2);
+orient tall, wysiwyg
+
+figure; nr=2; nc=2;
+ylab3 = 'norm of mu vector';
+
+subplot(nr,nc,1); semilogy(pinds,s,'.','markersize',20);
+grid on; xlabel(xlab1); ylabel(ylab1);
+subplot(nr,nc,2); semilogy(pinds,rss,'.','markersize',20);
+grid on; xlabel(xlab2); ylabel(ylab2);
+subplot(nr,nc,3); semilogy(pinds,mu_norm,'.','markersize',20);
+grid on; xlabel(xlab2); ylabel(ylab3);
+subplot(nr,nc,4); loglog(mu_norm,rss,'.-','markersize',20);
+grid on; xlabel(ylab3); ylabel(ylab2);
+orient tall, wysiwyg
+
+% mu vectors
+figure; nr=2; nc=1;
+subplot(nr,nc,1); plot(pinds,mu_all,'.-');
+grid on; xlabel('event index'); ylabel('elements of mu vectors');
+xlim([0 nsrc+1]);
+subplot(nr,nc,2); plot(mu_all','.-');
+grid on; xlabel(xlab2); ylabel('elements of mu vectors');
+xlim([0 nsrc+1]);
+orient tall, wysiwyg
+
+%=========================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_hessian.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_hessian.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_hessian.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -0,0 +1,61 @@
+%
+% wave2d_write_hessian.m
+% Carl Tape, 01-Feb-2009
+%
+%
+%
+% calls xxx
+% called by xxx
+%
+
+function wave2d_write_hessian(H,eids,odir,nprint)
+
+disp(['wave2d_write_hessian.m: writing files to ' odir]);
+
+[a,b] = size(H);
+if a ~= b, error('H must b square'); end
+nsrc = a;
+
+nprintmax = (nsrc-1)*nsrc/2;
+if nprint > nprintmax, error(sprintf('nprint (%i) must be less than nprintmax (%i)',nprint,nprintmax)); end
+
+% construct the matrix of data-covariance normalization
+i_ind = zeros(nsrc,nsrc);
+j_ind = zeros(nsrc,nsrc);
+for i = 1:nsrc
+ for j = 1:nsrc
+ i_ind(i,j) = i;
+ j_ind(i,j) = j;
+ end
+end
+
+% 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(:) ];
+bkeep = sum(Hcols,2)~=0; Hcols = Hcols(bkeep,:);
+[Hsort,iHsort_good] = sortrows(Hcols, -1);
+[Hsort,iHsort_bad] = sortrows(Hcols, 1);
+
+% make a list of the largest N positive and negative Hessian elements
+filename = [odir 'hessian_elements_good'];
+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 = [odir 'hessian_elements_bad'];
+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);
+
+%===================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_src.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_src.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_src.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -10,6 +10,8 @@
function wave2d_write_src(filename,xs0,ys0,ts,xs,ys,ts_res,xs_res,ys_res)
+disp(['writing to ' filename]);
+
% ! sources for synthetics
% write(20,'(8e20.10)') xtemp, ztemp, &
% m_src_syn_vec(itemp1), m_src_syn_vec(itemp2), m_src_syn_vec(itemp3), &
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_str.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_str.m 2011-04-12 00:05:16 UTC (rev 18229)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/matlab_scripts/wave2d_write_str.m 2011-04-13 16:44:53 UTC (rev 18230)
@@ -11,6 +11,8 @@
function wave2d_write_str(filename,x,y,kappa,mu,rho,B)
+disp(['writing to ' filename]);
+
% ! CURRENT MODEL (synthetics)
% write(19,'(6e20.10)') x_plot(iglob), z_plot(iglob), &
% kappa_syn(i,j,ispec), mu_syn(i,j,ispec), rho_syn(i,j,ispec), &
More information about the CIG-COMMITS
mailing list