[cig-commits] r13107 - in seismo/3D/ADJOINT_TOMO/iterate_adj/matlab: . MISFIT_in_development
carltape at geodynamics.org
carltape at geodynamics.org
Mon Oct 20 19:46:27 PDT 2008
Author: carltape
Date: 2008-10-20 19:46:27 -0700 (Mon, 20 Oct 2008)
New Revision: 13107
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/README
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt_run.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/plot_misfit.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/seis_norm.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/setup_misfit_dir.pl
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_misfit_dir.pl
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_subspace_dir.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
Log:
Updated output file capabilities for tabulating misfit values and values of the Hessian for the subspace method. Moved development scripts for misfit comparisons into MISFIT_in_development.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/README 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,7 @@
+Carl Tape, 20-Oct-2008
+
+/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj_work/matlab/MISFIT_in_development/README
+
+These scripts are in development. The objective is to have a way to combine synthetics and measurements from many different models into a series of plots.
+
+---------
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/compare_misfit.m (from rev 12980, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/compare_misfit.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/compare_misfit.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,302 @@
+%
+% compare_misfit.m
+% CARL TAPE, 23-July-2008
+% printed xxx
+%
+% This script reads in a set of window_chi output measurement files from
+% mt_measure_adj.f90, and it tabulates and plots misfit function values.
+%
+% /cig/seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/
+%
+% Copied from compute_misfit.m on 23-July-2008
+%
+% calls read_window_chi.m, plot_histo.m, readCMT.m, read_station_SPECFEM.m
+% called by xxx
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
+
+Tvec = [6 2];
+ % LOWEST period used for each measurement band-pass (upper is 40 s)
+comps = {'BHZ','BHR','BHT'};
+ % regardless of the component label for the DATA, the measurement code
+ % defaults so that the first two letters are always BH
+ncomp = length(comps);
+nper = length(Tvec);
+
+%DT = 0.011; % really this is only needed because it was left out from
+ % integrating the waveforms in mt-measure_adj.f90
+
+%-------------------------
+% USER INPUT
+
+iwrite = input(' Enter 1 to write files (0 otherwise) : ');
+
+% files: event IDs, CMTSOLUTIONS, stations
+dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_09/';
+%eid_file = [dir_source 'EIDs_only_loc'];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_m05';
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v09'];
+stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
+
+%-------------------------
+% read in list of event IDs and sources
+
+% load the event IDs corresponding to the kernels
+% load these as strings, since event IDs could have letters
+eids = textread(eid_file,'%s');
+eids_num = textread(eid_file,'%f'); % numeric version
+%eids_num = str2num(char(eids));
+nevent = length(eids);
+
+% read in CMT solutions
+[date,tshift,hdur,slat,slon,dep,M,eid,elabel] = readCMT(cmt_file_all, 13, 1);
+
+%-----------------------
+% read in stations file (used for synthetics)
+
+[rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(stations_file);
+nrec = length(stnm);
+for ii = 1:nrec, strec{ii} = [stnm{ii} '.' netwk{ii}]; end
+strec = strec(:);
+
+% make list of networks
+stnet = unique(netwk);
+nnet = length(stnet);
+
+%-----------------------
+
+iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
+
+imodel_min = input(' Enter minimum model index (0-5) : ');
+imodel_max = input(' Enter maximum model index (1-5) : ');
+nmodel = imodel_max - imodel_min;
+
+wdiff_array = zeros(nevent,nrec,ncomp,nper,nmodel);
+
+if iread == 1
+
+ % range of events
+ imin = 1; imax = nevent; % default
+ %imin = 4; imax = 7;
+ %imin = 4; imax = imin;
+
+ for imod = imodel_min:imodel_max
+
+ stmod = sprintf('m%2.2i',imod);
+ disp('====================================');
+ disp(['MODEL ' stmod]);
+
+ k2 = 0;
+ meas_array = [];
+ %for tt = 1:1
+ for tt = 1:nper
+ Tper = Tvec(tt);
+ sTper = sprintf('%3.3i',Tper);
+ if imod <= 2, sTper2 = sprintf('%2.2i',Tper); else sTper2 = sTper; end % CHT
+ disp('-----------------------------------------');
+
+ for ii = imin:imax % loop over events
+
+ disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
+ dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
+ mfile = [eids{ii} '_T' sTper2 '_window_chi'];
+ filename = [dir1 mfile];
+
+ if ~exist(filename,'file')
+ disp('--> file does not exist (or nwin = 0)');
+ else
+ % read the window_chi file
+ [stnet0,strec0,comp,iwin,iker,t1,t2,...
+ chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
+ measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
+ sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,...
+ wind2,wins2,win_chi0,windur, recd2,recs2,rec_chi0,recdur,...
+ tr_chi,am_chi] = read_window_chi(filename);
+
+ % waveform differences normalized by data-squared
+ win_chi = win_chi0 ./ wind2;
+ rec_chi = rec_chi0 ./ recd2;
+
+ nwin = length(strec0);
+ disp(['--> nwin = ' num2str(nwin)]);
+
+ % assign the string network an integer index
+ index_net = NaN * ones(nwin,1);
+ for inet = 1:nnet
+ ind_net = find( strcmp(stnet(inet), stnet0) == 1 );
+ index_net(ind_net) = inet;
+ end
+
+ % assign the string receiver an integer index
+ index_rec = NaN * ones(nwin,1);
+ for irec = 1:nrec
+ ind_rec = find( strcmp(strec(irec), strec0) == 1 );
+ index_rec(ind_rec) = irec;
+ end
+
+ % assign the string component an integer index
+ index_comp = NaN * ones(nwin,1);
+ for icomp = 1:ncomp
+ ind_comp = find( strcmp(comps(icomp), comp) == 1 );
+ index_comp(ind_comp) = icomp;
+ end
+
+ % measurement indices
+ k1 = k2+1;
+ k2 = k1+nwin-1;
+ kinds = [k1:k2]';
+
+ meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin ...
+ iker measCC_dT sigmaCC_dT measCC_dA sigmaCC_dA measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
+ % 1 kinds
+ % 2 index_T
+ % 3 index_event
+ % 4 index_network
+ % 5 index_rec
+ % 6 index_comp
+ % 7 iwin
+ % 8 iker
+ % 9 measCC_dT
+ % 10 sigmaCC_dT
+ % 11 measCC_dA
+ % 12 sigmaCC_dA
+ % 13 measMT_dT
+ % 14 sigmaMT_dT
+ % 15 win_chi
+ % 16 rec_chi
+ % 17 tr_chi
+ end
+
+ end
+ end
+
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ imatch = find( and( and( meas_array(:,2)==tt, meas_array(:,3)==ii), ...
+ and( meas_array(:,5)==jj, meas_array(:,6)==kk) ));
+ if ~isempty(imatch)
+ % take the first only
+ wdiff_array(ii,jj,kk,tt,imod+1) = meas_array(imatch(1),16);
+ end
+ end
+ end
+ end
+ end
+
+ end
+
+ save('wdiff_array.mat','wdiff_array');
+
+else
+ load('wdiff_array.mat');
+end
+
+whos wdiff_array
+
+% total number of windows
+N = length(meas_array);
+
+%----------------------------------------------
+% CODE IN PROGRESS
+% for each record for each model (m0, m1) that has at least one window picked,
+% we save the integrated waveform difference as the purest measure of misfit
+% between the data and synthetics
+
+[m1,m2,m3,m4,m5] = size(wdiff_array)
+
+if 0==1
+
+ ratio_array = []; ll = 0;
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ % check that there is a measurement for all models
+ if prod(wdiff_array(ii,jj,kk,tt,:)) > 0
+ ll = ll+1;
+ rat = wdiff_array(ii,jj,kk,tt,imodel_max+1) / wdiff_array(ii,jj,kk,tt,imodel_min+1);
+ ratio_array(ll,:) = [tt ii jj kk rat];
+ end
+ end
+ end
+ end
+ end
+
+ %ratio_sort = sortrows(ratio_array,5);
+ ratio_sort = sortrows(ratio_array,[-1 5]);
+
+ hmax = length(ratio_sort);
+ hmax = 100;
+ for hh = 1:hmax
+ tt = ratio_sort(hh,1);
+ ii = ratio_sort(hh,2);
+ jj = ratio_sort(hh,3);
+ kk = ratio_sort(hh,4);
+ disp(sprintf('%10s --> %6.1f %7s %5s %6.3f',char(eids{ii}),...
+ Tvec(tt),char(strec{jj}),char(comps{kk}),ratio_sort(hh,5) ));
+ end
+
+ %------------------------
+
+ rec_array = []; ll = 0;
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+
+ % check that there is a measurement for all components for all models
+ if prod( wdiff_array(ii,jj,:,tt,:) ) > 0
+ ll = ll+1;
+ rat = prod( wdiff_array(ii,jj,:,tt,imodel_max+1) ./ wdiff_array(ii,jj,:,tt,imodel_min+1) );
+ rec_array(ll,:) = [tt ii jj rat];
+ end
+ end
+
+ end
+ end
+ rec_sort = sortrows(rec_array,[-1 4]);
+
+ for hh = 1:length(rec_sort)
+ tt = rec_sort(hh,1);
+ ii = rec_sort(hh,2);
+ jj = rec_sort(hh,3);
+ disp(sprintf('%10s --> %6.1f %7s %5s %10.6f',char(eids{ii}),...
+ Tvec(tt),char(strec{jj}),rec_sort(hh,4) ));
+ end
+
+ %------------------------
+
+% % loop over all stations
+% tarray = []; ll = 0;
+% for jj = 1:nrec
+% imatch = find( ratio_array(:,3)==jj);
+% if length(imatch)==6
+% ll = ll+1;
+% tarray(ll,:) = [jj prod(ratio_array(imatch,5)) ];
+% end
+% end
+% tarray_sort = sortrows(tarray,2);
+% for hh = 1:20
+% jj = tarray_sort(hh,1);
+% disp(sprintf('%7s %6.3f',char(strec{jj}),tarray_sort(hh,2) ));
+% end
+
+ figure; N = plot_histo(ratio_array(:,end),[0:0.2:5]);
+ figure; N = plot_histo(log(ratio_array(:,end)),[-4:0.2:4]);
+
+ break
+end
+
+%=======================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/compare_misfit.m
___________________________________________________________________
Name: svn:mergeinfo
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,186 @@
+%
+% misfit_gmt.m
+% CARL TAPE, 03-Sept-2008
+% printed xxx
+%
+%
+% calls xxx
+% called by xxx
+%
+
+%-------------------------
+
+function norms = misfit_gmt(odir,Tmin,Tmax,smodel1,smodel2,eid,stnm,netwk)
+
+disp('running misfit_gmt.m...');
+
+
+%-------------------------
+
+ifigure = 0;
+
+Ttag = sprintf('T%3.3i_T%3.3i',Tmin,Tmax);
+ssuffix1 = ['semd.sac.d.' Ttag '.' smodel1];
+ssuffix2 = ['semd.sac.d.' Ttag '.' smodel2];
+dsuffix = ['sac.d.' Ttag];
+
+comps = {'Z','R','T'};
+%chans = {'HL','BL','HN','BN','HH','BH'};
+chans = {'HH','BH'};
+ncomp = length(comps);
+nchan = length(chans);
+
+%dir0 = pwd;
+%odir = [dir0 '/OUTPUT_MISFIT/EVENTS/' eid '/' stnm '.' netwk '/'];
+
+%-------------------------
+
+fhits = zeros(ncomp,3);
+ffiles = repmat(cellstr(''),ncomp,3);
+for icomp = 1:ncomp
+ for ichan = 1:nchan
+
+ comp = [chans{ichan} comps{icomp}];
+ dfile = [odir eid '.' netwk '.' stnm '.' comp '.' dsuffix];
+ sfile1 = [odir stnm '.' netwk '.' comp '.' ssuffix1];
+ sfile2 = [odir stnm '.' netwk '.' comp '.' ssuffix2];
+
+ if exist(dfile,'file')
+ fhits(icomp,1) = fhits(icomp,1)+1;
+ ffiles{icomp,1} = dfile;
+ end
+ if exist(sfile1,'file')
+ fhits(icomp,2) = fhits(icomp,2)+1;
+ ffiles{icomp,2} = sfile1;
+ end
+ if exist(sfile2,'file')
+ fhits(icomp,3) = fhits(icomp,3)+1;
+ ffiles{icomp,3} = sfile2;
+ end
+
+ end
+end
+fhits_sum = sum(fhits,2);
+
+if ifigure==1
+ figure; nr=5; nc=3;
+ st1 = ['syn(' smodel1 ')'];
+ st2 = ['syn(' smodel2 ')'];
+end
+
+norms = zeros(9,ncomp);
+
+for icomp = 1:ncomp
+%for icomp = 3
+
+ if fhits_sum(icomp) == 3
+ dfile = ffiles{icomp,1};
+ sfile1 = ffiles{icomp,2};
+ sfile2 = ffiles{icomp,3};
+
+ [seisD,HdrData,tnu,pobj,timeD] = readsac(dfile,0,'l');
+ [seisS1,HdrData,tnu,pobj,timeS1] = readsac(sfile1,0,'l');
+ [seisS2,HdrData,tnu,pobj,timeS2] = readsac(sfile2,0,'l');
+
+ ti = timeD(:);
+ n = length(ti);
+
+ % interpolate synthetics exactly onto the data time vector
+ % (In theory, they should already be within two time-steps at this point.)
+ s1 = interp1(timeS1,seisS1,ti,'linear','extrap');
+ s2 = interp1(timeS2,seisS2,ti,'linear','extrap');
+
+ % residuals for the two models; also perturbation to seismogram
+ res1 = s1 - seisD;
+ res2 = s2 - seisD;
+ spert = s2 - s1;
+
+ ymax = 1.05*max(abs([seisD' s1' s2' res1' res2']));
+ tmin = min(ti); tmax = max(ti);
+ ax0 = [tmin tmax -ymax ymax];
+
+ % compute nine different norms -- IS THIS EXACTLY WHAT WE WANT?
+ norms(1,icomp) = norm(seisD);
+ norms(2,icomp) = norm(s1);
+ norms(3,icomp) = norm(s2);
+ norms(4,icomp) = norm(res1);
+ norms(5,icomp) = norm(res2);
+ norms(6,icomp) = norm(s2-s1);
+ norms(7,icomp) = norm(res1)^2 / norm(seisD)^2; % chi^2
+ norms(8,icomp) = norm(res2)^2 / norm(seisD)^2; % chi^2
+ norms(9,icomp) = 100*(1 - (norm(res2)^2 / norm(res1)^2)); % VR
+
+ if ifigure==1
+ subplot(nr,nc,icomp);
+ plot(ti,s1,'r',ti,seisD,'b'); legend(st1,'data'); axis(ax0);
+ subplot(nr,nc,nc*1+icomp);
+ plot(ti,res1,'k'); legend([st1 ' - data']); axis(ax0);
+ subplot(nr,nc,nc*2+icomp);
+ plot(ti,s2,'r',ti,seisD,'b'); legend(st2,'data'); axis(ax0);
+ subplot(nr,nc,nc*3+icomp);
+ plot(ti,res2,'k'); legend([st2 ' - data']); axis(ax0);
+ subplot(nr,nc,nc*4+icomp);
+ plot(ti,spert,'k'); legend([st2 ' - ' st1]); axis(ax0);
+ end
+
+ % WRITE FILES FOR GMT PLOTTING
+ filename = [odir 'GMT_time_series_' Ttag '_' comps{icomp} '.dat'];
+ fid = fopen(filename,'w');
+ for ii = 1:n
+ fprintf(fid,[repmat('%12.4e',1,7) '\n'],...
+ ti(ii),seisD(ii),s1(ii),s2(ii),res1(ii),res2(ii),spert(ii) );
+ end
+ fclose(fid);
+
+ filename = [odir 'GMT_axes_' Ttag '_' comps{icomp} '.dat'];
+ fid = fopen(filename,'w');
+ fprintf(fid,[repmat('%12.4e',1,4) '\n'],...
+ ax0(1),ax0(2),ax0(3),ax0(4) );
+ fclose(fid);
+
+ filename = [odir 'GMT_norms_' Ttag '_' comps{icomp} '.dat'];
+ fid = fopen(filename,'w');
+ for kk = 1:9, fprintf(fid,'%12.4e\n',norms(kk,icomp) ); end
+ fclose(fid);
+
+ end
+
+end
+
+if ifigure == 1
+ orient tall, wysiwyg, fontsize(6)
+end
+
+return
+
+%-----------------------------------
+% if all three components exist, then compute additional time series for plots
+
+if sum(fhits_sum) == 9
+
+ % compute the time vector to use for all components
+ mint = zeros(ncomp,3);
+ maxt = zeros(ncomp,3);
+ for icomp = 1:ncomp
+ dfile = ffiles{icomp,1};
+ sfile1 = ffiles{icomp,2};
+ sfile2 = ffiles{icomp,3};
+ [seisD,HdrData,tnu,pobj,timeD] = readsac(dfile,0,'l');
+ [seisS1,HdrData,tnu,pobj,timeS1] = readsac(sfile1,0,'l');
+ [seisS2,HdrData,tnu,pobj,timeS2] = readsac(sfile2,0,'l');
+
+ mint(icomp,1) = min(timeD);
+ mint(icomp,2) = min(timeS1);
+ mint(icomp,3) = min(timeS2);
+ maxt(icomp,1) = max(timeD);
+ maxt(icomp,2) = max(timeS1);
+ maxt(icomp,3) = max(timeS2);
+ end
+
+ tmin = max(mint(:));
+ tmax = min(maxt(:));
+ ti = [tmin:0.05:tmax]';
+
+end
+
+%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt_run.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt_run.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/misfit_gmt_run.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,72 @@
+%
+% misfit_gmt_run.m
+% CARL TAPE, 03-Sept-2008
+% printed xxx
+%
+%
+% calls xxx
+% called by xxx
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+Tmin = 6;
+Tmax = 40;
+smodel1 = 'm00';
+smodel2 = 'm07';
+%stnm = 'USC'; netwk = 'CI';
+
+ncomp = 3;
+nper = 1;
+nmodel = 2;
+
+dir0 = pwd;
+dir1 = [dir0 '/OUTPUT_MISFIT/EVENTS/'];
+
+stations_file = [dir1 'STATIONS'];
+
+%-------------------------
+% read in list of event IDs and sources
+
+nevent = 1;
+eid = '14383980';
+dir2 = [dir1 eid '/'];
+
+%-----------------------
+% read in stations file (used for synthetics)
+
+[rlon0,rlat0,relev0,rburial0,stnm0,netwk0] = read_station_SPECFEM(stations_file);
+nrec = length(stnm0);
+%for ii = 1:nrec, strec{ii} = [stnm{ii} '.' netwk{ii}]; end
+%strec = strec(:);
+
+for irec = 1:nrec
+ disp(sprintf('%i out of %i -- %s.%s',irec,nrec,stnm0{irec},netwk0{irec}));
+end
+
+%-------------------------
+
+nevent = 1; ievent = 1;
+VR_array = zeros(nevent,nrec,ncomp,nper);
+
+for irec = 1:nrec
+
+ stnm = stnm0{irec};
+ netwk = netwk0{irec};
+ disp(sprintf('%i out of %i -- %s.%s',irec,nrec,stnm,netwk));
+
+ odir = [dir2 stnm '.' netwk '/'];
+
+ norms = misfit_gmt(odir,Tmin,Tmax,smodel1,smodel2,eid,stnm,netwk);
+
+ VR_array(ievent,irec,1,1) = norms(end,1); % Z
+ VR_array(ievent,irec,2,1) = norms(end,2); % R
+ VR_array(ievent,irec,3,1) = norms(end,3); % T
+
+end
+
+%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/plot_misfit.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/plot_misfit.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/plot_misfit.pl 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,349 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 02-Sept-2008
+# plot_misfit.pl
+#
+#
+# EXAMPLE:
+# plot_misfit.pl 9818433 USC/CI 6/40
+#
+#-----------------------------------
+
+#use Getopt::Std;
+#use POSIX;
+
+if (@ARGV < 3) {die("Usage: plot_misfit.pl eid sta/net Tmin/Tmax\n");}
+($eid,$stanet,$Ts) = @ARGV;
+
+#($sta,$net) = split(/\//,$stanet);
+($sta,$net) = split("/",$stanet);
+
+# bandpass range
+($Tmin,$Tmax) = split("/",$Ts);
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
+$Ttag = "${sTmin}_${sTmax}";
+
+#---------------------------
+# USER INPUT
+
+# models to compare
+$minm = 0;
+$maxm = 7;
+$smodel_min = sprintf("m%2.2i",$minm);
+$smodel_max = sprintf("m%2.2i",$maxm);
+
+# base directory (contains STATIONS and eid_list)
+#$dir0 = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj_work/matlab/OUTPUT_MISFIT";
+
+# CMTSOLUTION file
+$cmtfile = "../CMTSOLUTION";
+ at cmtvals = &get_cmt($cmtfile); # see subroutine below
+
+# STATIONS file
+$stafile = "../STATIONS_ADJOINT";
+
+# plot faults (southern California)
+$ifault = 1;
+
+# file to use to extract the SAC headers
+$Zsyn = "${sta}.${net}.BHZ.semd.sac.d.${Ttag}.${smodel_max}";
+
+$saclst = "/opt/seismo-util/bin/saclst";
+
+# axes positions
+$x0_title = 0.5; $y0_title = 10.5;
+$dY_under_title = 0.4;
+$dY_under_subtitle = 0.4;
+$dX_between_col = 0.2;
+$dY_under_seis = 0.3;
+
+# dimensions of seismogram subplots
+$xwid = 2.4;
+$ywid = 1.0;
+$J = "-JX${xwid}/${ywid}";
+
+# subtitles
+ at subs = ("(a) Vertical component, Z","(b) Radial component, R","(c) Transverse component, T","(d) Variance reduction for each component","(e) Source-receiver geometry");
+
+ at slabs = ("d and s($smodel_min)","d and s($smodel_max)","s($smodel_min) - d","s($smodel_max) - d","s($smodel_max) - s($smodel_min)");
+$nrows = 5; # <= 5
+
+$x0_map = 5.0; $y0_map = 2.3 + (5-$nrows)*$ywid;
+$x0_misfit = 1.5; $y0_misfit = 2.9 + (5-$nrows)*$ywid;
+
+#---------------------------
+
+# create psmeca file from CMTSOLUTION file
+ at M[0..5] = @cmtvals[12..17];
+$year = $cmtvals[0]; $month = $cmtvals[1]; $day = $cmtvals[2];
+$elat = $cmtvals[9]; $elon = $cmtvals[10]; $edep = $cmtvals[11];
+$eid = $cmtvals[18];
+$cmtpsmeca = "temp.txt";
+open(MIJ,">$cmtpsmeca");
+#printf MIJ "$elon $elat $edep $M[0] $M[1] $M[2] $M[3] $M[4] $M[5] 1 X Y $eid";
+printf MIJ "$elon $elat $edep $M[0] $M[1] $M[2] $M[3] $M[4] $M[5] 1";
+close(MIJ);
+
+# compute moment using Dahlen and Tromp, Eq. (5.91)
+# compute moment magnitude using Kanamori (1977) -- see Carl's Latex notes cmt.pdf
+$M0 = 1/sqrt(2) * sqrt($M[0]**2 + $M[1]**2 + $M[2]**2 + 2*($M[3]**2 + $M[4]**2 + $M[5]**2 ) );
+$k = -11.2 + (2/3)*log(5)/log(10);
+$Mw = (2/3)*log($M0)/log(10) + $k;
+$stM0 = sprintf("M0 = %.3e dyne-cm",$M0);
+$stMw = sprintf("Mw = %.2f",$Mw);
+print "\n moment : $stM0 \n moment magnitude : $stMw\n";
+
+# get the receiver location, azimuth, distance (from the $Zsyn file)
+# We want the azimuth for the file name, so that all the PDFs will be sorted by azimuth.
+(undef,$tmin0,$tmax0,$slon,$slat,$dist,$az,$Tchan) = split(" ",`$saclst b e stlo stla dist az kcmpnm f $Zsyn`);
+$edate = sprintf("%4.4i . %2.2i . %2.2i",$year,$month,$day);
+$stedep = sprintf("depth = %.1f km",$edep);
+$strec = "$sta . $net";
+$stdist = sprintf("dist = %.1f km",$dist);
+$staz = sprintf("az = %.1f deg",$az);
+
+#---------------------------
+
+$fontno = 1; # font number to use for text
+
+# set GMT defaults
+`gmtset BASEMAP_TYPE plain PAPER_MEDIA letter ANNOT_FONT_SIZE = 10p LABEL_FONT_SIZE = 10p PAGE_ORIENTATION = portrait PLOT_DEGREE_FORMAT D TICK_LENGTH 4p`;
+
+# name of the output files
+$azlabel = sprintf("%3.3i",$az);
+$name = "${azlabel}_${sta}_${net}_${Ttag}";
+$psfile = "$name.ps";
+$jpg_file = "$name.jpg";
+
+# overall title
+$title = "Event $eid, station ${sta}.${net}, bandpass ${Tmin}-${Tmax} s, models ${smodel_min} and ${smodel_max}";
+`pstext -R0/1/0/1 -JX1 -N -Xa${x0_title} -Ya${y0_title} -K -P > $psfile <<EOF\n 0 0 14 0 1 LM $title \nEOF\n`; # START (no -O)
+
+# subtitles
+$x0_sub = $x0_title;
+$y0_sub = $y0_title - $dY_under_title;
+
+for ($k = 1; $k <= 3; $k = $k+1) {
+ $x_sub = $x0_sub + ($k-1)*($xwid+$dX_between_col);
+ $origin_sub = "-Xa${x_sub} -Ya${y0_sub}";
+ `pstext -R0/1/0/1 -JX1 -N $origin_sub -K -O >> $psfile <<EOF\n 0 0 10 0 1 LM $subs[$k-1] \nEOF\n`;
+}
+
+# subtitle for misfit plot
+$x0_sub2 = $x0_misfit - 0.8;
+$y0_sub2 = $y0_title - $dY_under_title - $dY_under_subtitle - $nrows*$ywid - $dY_under_seis;
+$origin_sub = "-Xa${x0_sub2} -Ya${y0_sub2}";
+`pstext -R0/1/0/1 -JX1 -N $origin_sub -K -O >> $psfile <<EOF\n 0 0 10 0 1 LM $subs[3] \nEOF\n`;
+
+# subtitle for map
+$x0_sub3 = $x0_map - 0.8;
+$y0_sub3 = $y0_sub2;
+$origin_sub = "-Xa${x0_sub3} -Ya${y0_sub3}";
+`pstext -R0/1/0/1 -JX1 -N $origin_sub -K -O >> $psfile <<EOF\n 0 0 10 0 1 LM $subs[4] \nEOF\n`;
+
+#---------------------------------------------------------------------------
+# PLOT SEISMOGRAMS AND RESIDUALS
+
+$x0_seis = $x0_title;
+$y0_seis = $y0_title - $dY_under_title - $dY_under_subtitle - $ywid;
+
+$B0 = "-Ba50f10g50:\" \":/a1f1::wesN";
+$B = "-Ba50f10g50:\" \":/a1f1::wesn";
+
+#$slabinfo = "-C2p -W0/255/255o,0/255/255 -G0/0/0 -N";
+$slabinfo = "-C3p -W255/255/255o,0/0/0 -G0/0/0 -N";
+
+ at comps = ("Z","R","T");
+for ($k = 1; $k <= 3; $k = $k+1) {
+ $comp = $comps[$k-1];
+
+ $Pfile = "GMT_time_series_${Ttag}_${comp}.dat";
+ $bounds = "GMT_axes_${Ttag}_${comp}.dat";
+ $Nfile = "GMT_norms_${Ttag}_${comp}.dat";
+
+ # if the seismograms exist, then plot them
+ if ( (-f $Pfile) && (-f $bounds) ) {
+
+ open(IN,"$bounds"); @lines = <IN>;
+ ($tmin,$tmax,$ymin,$ymax) = split(" ",$lines[0]);
+ $tmin = -10;
+ if($Tmin==2) {$tmax = $tmax/2;} # USER: reduce the records for 2s
+ $R = "-R$tmin/$tmax/$ymin/$ymax";
+
+ # plot seismograms
+ $x_seis = $x0_seis + ($k-1)*($xwid+$dX_between_col);
+ $x_seis_lab = $x_seis + $xwid - 0.05;
+
+ for ($j = 1; $j <= $nrows; $j = $j+1) {
+ $y_seis = $y0_seis - ($j-1)*$ywid;
+ $originP = "-Xa${x_seis} -Ya${y_seis}";
+
+ $y_seis_lab = $y_seis + $ywid - 0.05;
+ $originPlab = "-Xa${x_seis_lab} -Ya${y_seis_lab}";
+
+ if($j==1) {`psbasemap $J $R $B0 -O -K $originP >> $psfile`}
+ else {`psbasemap $J $R $B -O -K $originP >> $psfile`}
+ if ($j==1) {
+ `awk \'{print \$1,\$2}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `awk \'{print \$1,\$3}\' $Pfile | psxy $J $R -W0.5p,255/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[0]\nEOF\n`;
+ } elsif ($j==2) {
+ `awk \'{print \$1,\$2}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `awk \'{print \$1,\$4}\' $Pfile | psxy $J $R -W0.5p,255/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[1]\nEOF\n`;
+ } elsif ($j==3) {
+ `awk \'{print \$1,\$5}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[2]\nEOF\n`;
+ } elsif ($j==4) {
+ `awk \'{print \$1,\$6}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[3]\nEOF\n`;
+ } elsif ($j==5) {
+ `awk \'{print \$1,\$7}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[4]\nEOF\n`;
+ }
+ }
+
+ } else {
+ print "seismogram files $Pfile and $bounds do not both exist\n";
+ }
+
+}
+
+#---------------------------------------------------------------------------
+# PLOT THE MISFIT
+
+$originN = "-Xa${x0_misfit} -Ya${y0_misfit}";
+
+$xwid = 1.5;
+$ywid = 1.3;
+$J = "-JX${xwid}/${ywid}";
+$B = "-Ba1f1:\"Model index\":/a50f10:\"Variance reduction\":WeSn";
+$kmax = $maxm + 3;
+$R = "-R-0.5/$kmax/-50/100";
+
+$Zinfo = "-Sc12p -W0.5p,0/0/0 -G255 -N";
+
+print "J option --> $J\n";
+print "R option --> $R\n";
+print "B option --> $B\n";
+
+ at comps = ("Z","R","T");
+for ($k = 1; $k <= 3; $k = $k+1) {
+ $comp = $comps[$k-1];
+
+ $Nfile = "GMT_norms_${Ttag}_${comp}.dat";
+
+ # if the norms file exists
+ if ( -f $Nfile ) {
+ # reduction in the misfit
+ open(IN,"$Nfile"); @Nlines = <IN>;
+ #$yred = 100*(1.0 - $Nlines[4]/$Nlines[3]);
+ $yred = $Nlines[8]; chomp($yred);
+ print "yred --> $yred\n";
+
+ # plot the misfit
+ if($k==1) {`psbasemap $J $R $B -O -K $originN >> $psfile`}
+ `psxy $R $J -W1p -K -O $originN >> $psfile <<EOF\n$minm 0\n$maxm $yred\nEOF\n`;
+ `psxy $R $J $Zinfo -K -O $originN >> $psfile <<EOF\n$minm 0\nEOF\n`;
+ `psxy $R $J $Zinfo -K -O $originN >> $psfile <<EOF\n$maxm $yred\nEOF\n`;
+ `pstext $R $J -N -K -O $originN >> $psfile <<EOF\n$maxm $yred 8 0 $fontno CM $comp\nEOF\n`;
+
+ } else {
+ print "norm file $Nfile does not exist\n";
+ }
+
+}
+
+#---------------------------------------------------------------------------
+# PLOT BASE MAP with stations, faults, CMT
+
+$originM = "-Xa${x0_map} -Ya${y0_map}";
+$proj = "-JM2.5";
+$textinfo = "-G0 -S1p,255/255/255";
+
+$xmin = -122; $xmax = -114;
+$ymin = 32; $ymax = 37;
+$bounds="-R$xmin/$xmax/$ymin/$ymax ";
+$xtick = 1; $ytick = 1;
+$tick="-B${xtick}/${ytick}::wesn";
+$cinfo="-W0.5 -Dh -A100";
+
+# plot the base map
+`psbasemap $bounds $proj $tick $originM -K -O >> $psfile`;
+`pscoast $bounds $proj $cinfo -S150/200/255 -W2 -G200/255/150 -N1/1p -N2/0.5p -O -K $originM >> $psfile`;
+
+# plot southern California faults
+if($ifault==1) {
+ $dir1 = "/home/carltape/gmt";
+ $fault_file = "$dir1/faults/jennings.xy";
+ $kcf_file = "$dir1/faults/kcf.xy";
+ $breck_file = "$dir1/faults/breck.xy";
+ $fault_infoK = "-M -W0.5p,0/0/0";
+ `psxy $fault_file $bounds $proj $fault_infoK -O -K $originM >> $psfile`;
+ `psxy $kcf_file $bounds $proj $fault_infoK -O -K $originM >> $psfile`;
+ `psxy $breck_file $bounds $proj $fault_infoK -O -K $originM >> $psfile`;
+}
+
+# plot the stations
+`awk \'{print \$4,\$3,0,4,B,C}\' $stafile | psxy $proj $bounds -M -St0.10 -W0.25p -G200/200/200 -N -O -K $originM >> $psfile`;
+
+#------------
+
+# plot labels next to the map
+ at mapstrings = ($edate,$eid,$stMw,$stedep,$strec,$stdist,$staz);
+$nstrings = @mapstrings;
+$x1 = $xmin - 4*$xtick;
+for ($i=0; $i < $nstrings; $i++) {
+ $y1 = $ymax - 0.5 - $i*$ytick*0.7;
+ `pstext $bounds $proj $tick -N -K -O $originM >> $psfile <<EOF\n $x1 $y1 10 0 $fontno LM $mapstrings[$i] \nEOF\n`;
+}
+
+# plot ray path, station, and CMT source
+`psxy $bounds $proj -W1p -K -O $originM >> $psfile <<EOF\n$elon $elat\n$slon $slat\nEOF\n`; # ray path
+#$cmtfont = 10; $cmt_scale = 0.3; $cmtinfo = "-Sm${cmt_scale}/$cmtfont -L0.5p/0/0/0 -G255/0/0";
+#`psmeca $cmtpsmeca $bounds $proj $cmtinfo -O -K $originM >> $psfile`; # CMT source
+$cmt_scale = 0.3; $cmtinfo = "-Sm${cmt_scale} -L0.5p/0/0/0 -G255/0/0";
+`psmeca $cmtpsmeca $bounds $proj $cmtinfo -O -K $originM >> $psfile`; # CMT source
+$elat_shift = $elat + 0.6;
+`pstext $bounds $proj $textinfo -K -O $originM >> $psfile <<EOF\n $elon $elat_shift 9 0 $fontno CM $eid\nEOF\n`; # station label
+`rm $cmtpsmeca`;
+
+`psxy $bounds $proj -St0.15 -W2 -G255/0/0 -K -O $originM >> $psfile <<EOF\n$slon $slat\nEOF\n`; # red station
+$slat_shift = $slat - 0.3;
+`pstext $bounds $proj $textinfo -K -O $originM >> $psfile <<EOF\n $slon $slat_shift 8 0 $fontno CM $sta \nEOF\n`; # station label
+
+#---------------------------------------------------------------------------
+
+`pstext -R0/1/0/1 -JX1 -N -Xa${x0_title} -Ya${y0_title} -O >> $psfile <<EOF\n 10 10 14 0 1 LM $title\nEOF\n`; # FINISH (no -K)
+
+#`convert $psfile $jpg_file`; `xv $jpg_file &`;
+`ps2pdf $psfile`;
+#`rm $psfile`; # remove PS file
+#system("xv $jpg_file &");
+
+#================================================
+sub get_cmt {
+ my ($cmtfile)=@_;
+ open(CMT,"$cmtfile") or die("Error opening cmtfile $cmtfile\n");
+ @cmt = <CMT>;
+ my($pde,$oyear,$omonth,$oday,$ohr,$omin,$osec1)=split(" ",$cmt[0]);
+ my($osec,$omsec)=split(/\./,$osec1); $omsec=$omsec*10;
+ my(undef,undef,$eid)=split(" ",$cmt[1]);
+ my(undef,undef,$tshift)=split(" ",$cmt[2]);
+ my(undef,undef,$hdur)=split(" ",$cmt[3]);
+ my(undef,$elat)=split(" ",$cmt[4]);
+ my(undef,$elon)=split(" ",$cmt[5]);
+ my(undef,$edep)=split(" ",$cmt[6]);
+ my(undef,$Mrr)=split(" ",$cmt[7]);
+ my(undef,$Mtt)=split(" ",$cmt[8]);
+ my(undef,$Mpp)=split(" ",$cmt[9]);
+ my(undef,$Mrt)=split(" ",$cmt[10]);
+ my(undef,$Mrp)=split(" ",$cmt[11]);
+ my(undef,$Mtp)=split(" ",$cmt[12]);
+ close(CMT);
+ return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp,$eid);
+}
+
+#================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/plot_misfit.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/seis_norm.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/seis_norm.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/seis_norm.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,24 @@
+%
+% function x = seis_norm(t,y)
+% Carl Tape, 03-Sept-2008
+%
+% This function computes the L2 norm of a seismogram.
+%
+% t2 2 N 2 N 2
+% INT s(t) dt dt SUM (s_i) SUM (s_i)
+% 2 t1 i=1 i=1
+% x = ------------------- = ---------------- = ---------
+% t2 - t1 dt*N N
+%
+% Note that the number will not depend on the length of the time series.
+%
+% calls xxx
+% called by xxx
+%
+
+function x = seis_norm(t,y)
+
+x = sqrt( sum(y.^2) / length(y) );
+
+%============================================================
+
\ No newline at end of file
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/setup_misfit_dir.pl (from rev 12980, seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_misfit_dir.pl)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/setup_misfit_dir.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/setup_misfit_dir.pl 2008-10-21 02:46:27 UTC (rev 13107)
@@ -0,0 +1,304 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 03-Sept-2008
+# setup_misfit_dir.pl
+#
+# This script xxx
+#
+# EXAMPLE: setup_misfit_dir.pl m00 m07 6/40
+#
+#-----------------------------------
+
+if (@ARGV < 3) {die("Usage: setup_misfit_dir.pl smodel_min smodel_max Tmin/Tmax\n")}
+($minm,$maxm,$Ts) = @ARGV;
+
+$pwd = $ENV{PWD};
+print "\n$pwd\n";
+
+# labels for bandpass-filtered data and synthetics
+($Tmin,$Tmax) = split("/",$Ts);
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
+$Ttag = "${sTmin}_${sTmax}";
+
+#-------------------------------------------
+# USER INPUT
+
+# directory containing all CMTSOLUTION files
+$dir_src = "/net/sierra/raid1/carltape/results/SOURCES/socal_09";
+$dir_cmt = "${dir_src}/CMT_files_pre_inverted";
+if (not -e $dir_src) {die("check if dir_src $dir_src exist or not\n")}
+if (not -e $dir_cmt) {die("check if dir_cmt $dir_cmt exist or not\n")}
+
+# list of event IDs to use
+$file_eid0 = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_m07";
+if (not -f ${file_eid0}) {die("check if file_eid ${file_eid0} exist or not\n")}
+
+# FULL stations file
+$file_stations0 = "/net/denali/home2/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
+if (not -f ${file_stations0}) {die("check if file_stations ${file_stations0} exist or not\n")}
+
+# data and synthetics directories
+$dir_data0 = "/net/sierra/raid1/carltape/socal/socal_3D/DATA/FINAL";
+$dir_syn01 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${minm}";
+$dir_syn02 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${maxm}";
+if (not -e ${dir_data0}) {die("check if dir_data ${dir_data0} exist or not\n")}
+if (not -e ${dir_syn01}) {die("check if dir_syn ${dir_syn01} exist or not\n")}
+if (not -e ${dir_syn02}) {die("check if dir_syn ${dir_syn02} exist or not\n")}
+
+# directory containing the output for windows, measurements, adjoint sources, etc
+$dir_run = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
+if (not -e ${dir_run}) {die("check if dir_run ${dir_run} exist or not\n")}
+
+#-------------------------------------------
+
+# make the base output directory if it does not exist
+$dir0 = "OUTPUT_MISFIT";
+`mkdir -p $dir0`;
+$dir1 = "$dir0/EVENTS";
+`mkdir -p $dir1`;
+
+# link various files
+$file_eid = "$dir0/eid_list";
+`ln -s ${file_eid0} ${file_eid}`;
+$file_stations = "$dir0/STATIONS";
+`ln -s ${file_stations0} ${file_stations}`;
+
+# open EID list
+open(IN,"${file_eid}"); @eids = <IN>; close(IN);
+$nevent = @eids;
+print "\n $nevent events in the list";
+
+# open STATIONS file
+open(IN,"${file_stations}"); @slines = <IN>; close(IN);
+$nrec = @slines - 1;
+print "\n $nrec stations in the list";
+
+#-------------------------------------------
+
+if (0==1) {
+ for ($ievent = 1; $ievent <= $nevent; $ievent++) {
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ }
+ die("testing");
+}
+
+#=============================================
+# MAKE DIRECTORIES
+
+#$imin = 1; $imax = $nevent;
+#$imin = 1; $imax = 100;
+$imin = 142; $imax = $imin;
+
+$isetup1 = 0;
+if ($isetup1 == 1) {
+
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+
+ # make base directory for each event
+ $dir2 = "$dir1/$eid";
+ `mkdir -p $dir2`;
+
+ # link CMTSOLUTION file
+ $cmtfile0 = "${dir_cmt}/CMTSOLUTION_${eid}";
+ $cmtfile = "$dir2/CMTSOLUTION";
+ `ln -s $cmtfile0 $cmtfile`;
+
+ # link STATIONS file from the windowing/measurement codes (different for each period range)
+ # --> use the max model to show the stations used
+ $stafile0 = "${dir_run}/$eid/$maxm/MEASURE_${sTmin}/ADJOINT_${sTmin}/STATIONS_ADJOINT";
+ if (not -f $stafile0) {die("check if stafile0 $stafile0 exist or not\n")}
+ $stafile = "$dir2/STATIONS_ADJOINT";
+ `ln -s $stafile0 $stafile`;
+
+ # make a directory for each station in the STATIONS list
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $dir_rec = "$dir2/$sta.$net";
+ print "\n $irec out of $nrec -- $dir_rec";
+ `mkdir -p $dir_rec`;
+ }
+ }
+
+}
+
+#=============================================
+# COPY IN DATA AND SYNTHETICS
+# NOTE: For now we are copying in all available records,
+# but really we only want BHZ,BHR,BHT (or HHZ,HHR,HHT).
+
+$isetup2 = 0;
+
+if ($isetup2 == 1) {
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+
+ # data and synthetics directories for each event
+ $dir_data = "${dir_data0}/$eid/PROCESSED/PROCESSED_${Ttag}";
+ $dir_syn1 = "${dir_syn01}/$eid/PROCESSED/PROCESSED_${Ttag}";
+ $dir_syn2 = "${dir_syn02}/$eid/PROCESSED/PROCESSED_${Ttag}";
+ if (not -e ${dir_data}) {die("check if dir_data ${dir_data} exist or not\n");}
+ if (not -e ${dir_syn1}) {die("check if dir_syn ${dir_syn1} exist or not\n");}
+ if (not -e ${dir_syn2}) {die("check if dir_syn ${dir_syn2} exist or not\n");}
+
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ #for ($irec = 128; $irec <= 128; $irec++) {
+
+ cd_directory($pwd);
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ # link plotting script
+ $plotfile0 = "${pwd}/plot_misfit.pl";
+ $plotfile = "${dir_rec}/plot_misfit.pl";
+ `rm $plotfile`;
+ `ln -s $plotfile0 $plotfile`;
+
+ $dtags = "${dir_data}/*$net.$sta*"; @dfiles = glob($dtags); `cp @dfiles ${dir_rec}`;
+
+ $stags1 = "${dir_syn1}/$sta.$net*"; @sfiles1 = glob($stags1); #`cp @sfiles1 ${dir_rec}`;
+ cd_directory($dir_rec);
+ `cp @sfiles1 .`;
+ `mv_files.pl -x "*semd.sac.d.${Ttag}" "*semd.sac.d.${Ttag}.$minm"`;
+
+ $stags2 = "${dir_syn2}/$sta.$net*"; @sfiles2 = glob($stags2); #`cp @sfiles2 ${dir_rec}`;
+ `cp @sfiles2 .`;
+ `mv_files.pl -x "*semd.sac.d.${Ttag}" "*semd.sac.d.${Ttag}.$maxm"`;
+ }
+ }
+
+}
+
+#=============================================
+# NOW RUN MATLAB CODE misfit_gmt_run.m
+
+#=============================================
+# GENERATE GMT FIGURES
+
+$isetup3 = 0;
+
+if ($isetup3 == 1) {
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ #for ($irec = 104; $irec <= 104; $irec++) {
+
+ cd_directory($pwd);
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ cd_directory($dir_rec);
+ @datfiles = glob("*.dat");
+ $ndat = @datfiles;
+
+ if($ndat > 0) {
+ `plot_misfit.pl $eid $sta/$net $Ts`;
+ }
+ }
+ }
+
+}
+
+#=============================================
+# CONCATENATE PDF FIGURES
+
+$isetup4 = 0;
+
+if ($isetup4 == 1) {
+
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ @pdffiles = glob("${dir_rec}/*${Ttag}.pdf");
+ $nfile = @pdffiles;
+ #$pdffile = "${dir_rec}/${sta}_${net}_${Ttag}.pdf";
+ if($nfile == 1) {
+ $pdffile = $pdffiles[0]; chomp($pdffile);
+ `cp $pdffile $dir2`;
+ }
+ }
+
+ # concatenate into a single pdf
+ $outfile = "$dir1/${eid}_${Ttag}_waveform_az.pdf";
+ `rm $outfile`;
+ @pdfall = glob("${dir2}/*${Ttag}.pdf");
+ print "\n /home/carltape/bin/pdcat -r @pdfall $outfile\n";
+ `/home/carltape/bin/pdcat -r @pdfall $outfile`;
+ `rm @pdfall`;
+ }
+}
+
+#=============================================
+# DELETE FILES
+
+$iclean = 0;
+
+if ($iclean == 1) {
+
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ #for ($irec = 128; $irec <= 128; $irec++) {
+
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ `rm ${dir_rec}/*pdf`; `rm ${dir_rec}/*ps`;
+ }
+ }
+}
+
+#================================================
+print "\n done with setup_misfit_dir.pl\n\n";
+#================================================
+
+sub cd_directory {
+ my($new) = @_;
+ my $old = `pwd`;
+ chomp($old);
+ check_directory($new);
+ #print "$prog: cd $new\n";
+ chdir $new;
+ return($old);
+}
+
+sub check_directory {
+ if(! -e $_[0] ) {
+ print "Directory not found: $_[0]\n";
+ exit(-1);
+ }
+}
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/setup_misfit_dir.pl
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2008-10-21 01:21:15 UTC (rev 13106)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -1,302 +0,0 @@
-%
-% compare_misfit.m
-% CARL TAPE, 23-July-2008
-% printed xxx
-%
-% This script reads in a set of window_chi output measurement files from
-% mt_measure_adj.f90, and it tabulates and plots misfit function values.
-%
-% /cig/seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/
-%
-% Copied from compute_misfit.m on 23-July-2008
-%
-% calls read_window_chi.m, plot_histo.m, readCMT.m, read_station_SPECFEM.m
-% called by xxx
-%
-
-%-------------------------
-
-clear
-close all
-format compact
-
-dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
-
-Tvec = [6 2];
- % LOWEST period used for each measurement band-pass (upper is 40 s)
-comps = {'BHZ','BHR','BHT'};
- % regardless of the component label for the DATA, the measurement code
- % defaults so that the first two letters are always BH
-ncomp = length(comps);
-nper = length(Tvec);
-
-%DT = 0.011; % really this is only needed because it was left out from
- % integrating the waveforms in mt-measure_adj.f90
-
-%-------------------------
-% USER INPUT
-
-iwrite = input(' Enter 1 to write files (0 otherwise) : ');
-
-% files: event IDs, CMTSOLUTIONS, ststions
-dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_8/';
-%eid_file = [dir_source 'EIDs_only_loc'];
-%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
-eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_m05';
-cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v8'];
-stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
-
-%-------------------------
-% read in list of event IDs and sources
-
-% load the event IDs corresponding to the kernels
-% load these as strings, since event IDs could have letters
-eids = textread(eid_file,'%s');
-eids_num = textread(eid_file,'%f'); % numeric version
-%eids_num = str2num(char(eids));
-nevent = length(eids);
-
-% read in CMT solutions
-[date,tshift,hdur,slat,slon,dep,M,eid,elabel] = readCMT(cmt_file_all, 13, 1);
-
-%-----------------------
-% read in stations file (used for synthetics)
-
-[rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(stations_file);
-nrec = length(stnm);
-for ii = 1:nrec, strec{ii} = [stnm{ii} '.' netwk{ii}]; end
-strec = strec(:);
-
-% make list of networks
-stnet = unique(netwk);
-nnet = length(stnet);
-
-%-----------------------
-
-iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
-
-imodel_min = input(' Enter minimum model index (0-5) : ');
-imodel_max = input(' Enter maximum model index (1-5) : ');
-nmodel = imodel_max - imodel_min;
-
-wdiff_array = zeros(nevent,nrec,ncomp,nper,nmodel);
-
-if iread == 1
-
- % range of events
- imin = 1; imax = nevent; % default
- %imin = 4; imax = 7;
- %imin = 4; imax = imin;
-
- for imod = imodel_min:imodel_max
-
- stmod = sprintf('m%2.2i',imod);
- disp('====================================');
- disp(['MODEL ' stmod]);
-
- k2 = 0;
- meas_array = [];
- %for tt = 1:1
- for tt = 1:nper
- Tper = Tvec(tt);
- sTper = sprintf('%3.3i',Tper);
- if imod <= 2, sTper2 = sprintf('%2.2i',Tper); else sTper2 = sTper; end % CHT
- disp('-----------------------------------------');
-
- for ii = imin:imax % loop over events
-
- disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
- dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
- mfile = [eids{ii} '_T' sTper2 '_window_chi'];
- filename = [dir1 mfile];
-
- if ~exist(filename,'file')
- disp('--> file does not exist (or nwin = 0)');
- else
- % read the window_chi file
- [stnet0,strec0,comp,iwin,iker,t1,t2,...
- chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
- measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
- sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,...
- wind2,wins2,win_chi0,windur, recd2,recs2,rec_chi0,recdur,...
- tr_chi,am_chi] = read_window_chi(filename);
-
- % waveform differences normalized by data-squared
- win_chi = win_chi0 ./ wind2;
- rec_chi = rec_chi0 ./ recd2;
-
- nwin = length(strec0);
- disp(['--> nwin = ' num2str(nwin)]);
-
- % assign the string network an integer index
- index_net = NaN * ones(nwin,1);
- for inet = 1:nnet
- ind_net = find( strcmp(stnet(inet), stnet0) == 1 );
- index_net(ind_net) = inet;
- end
-
- % assign the string receiver an integer index
- index_rec = NaN * ones(nwin,1);
- for irec = 1:nrec
- ind_rec = find( strcmp(strec(irec), strec0) == 1 );
- index_rec(ind_rec) = irec;
- end
-
- % assign the string component an integer index
- index_comp = NaN * ones(nwin,1);
- for icomp = 1:ncomp
- ind_comp = find( strcmp(comps(icomp), comp) == 1 );
- index_comp(ind_comp) = icomp;
- end
-
- % measurement indices
- k1 = k2+1;
- k2 = k1+nwin-1;
- kinds = [k1:k2]';
-
- meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin ...
- iker measCC_dT sigmaCC_dT measCC_dA sigmaCC_dA measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
- % 1 kinds
- % 2 index_T
- % 3 index_event
- % 4 index_network
- % 5 index_rec
- % 6 index_comp
- % 7 iwin
- % 8 iker
- % 9 measCC_dT
- % 10 sigmaCC_dT
- % 11 measCC_dA
- % 12 sigmaCC_dA
- % 13 measMT_dT
- % 14 sigmaMT_dT
- % 15 win_chi
- % 16 rec_chi
- % 17 tr_chi
- end
-
- end
- end
-
- for tt = 1:nper
- for ii = imin:imax
- disp(sprintf('Event %i out of %i',ii,nevent));
- for jj = 1:nrec
- for kk = 1:ncomp
- imatch = find( and( and( meas_array(:,2)==tt, meas_array(:,3)==ii), ...
- and( meas_array(:,5)==jj, meas_array(:,6)==kk) ));
- if ~isempty(imatch)
- % take the first only
- wdiff_array(ii,jj,kk,tt,imod+1) = meas_array(imatch(1),16);
- end
- end
- end
- end
- end
-
- end
-
- save('wdiff_array.mat','wdiff_array');
-
-else
- load('wdiff_array.mat');
-end
-
-whos wdiff_array
-
-% total number of windows
-N = length(meas_array);
-
-%----------------------------------------------
-% CODE IN PROGRESS
-% for each record for each model (m0, m1) that has at least one window picked,
-% we save the integrated waveform difference as the purest measure of misfit
-% between the data and synthetics
-
-[m1,m2,m3,m4,m5] = size(wdiff_array)
-
-if 0==1
-
- ratio_array = []; ll = 0;
- for tt = 1:nper
- for ii = imin:imax
- disp(sprintf('Event %i out of %i',ii,nevent));
- for jj = 1:nrec
- for kk = 1:ncomp
- % check that there is a measurement for all models
- if prod(wdiff_array(ii,jj,kk,tt,:)) > 0
- ll = ll+1;
- rat = wdiff_array(ii,jj,kk,tt,imodel_max+1) / wdiff_array(ii,jj,kk,tt,imodel_min+1);
- ratio_array(ll,:) = [tt ii jj kk rat];
- end
- end
- end
- end
- end
-
- %ratio_sort = sortrows(ratio_array,5);
- ratio_sort = sortrows(ratio_array,[-1 5]);
-
- hmax = length(ratio_sort);
- hmax = 100;
- for hh = 1:hmax
- tt = ratio_sort(hh,1);
- ii = ratio_sort(hh,2);
- jj = ratio_sort(hh,3);
- kk = ratio_sort(hh,4);
- disp(sprintf('%10s --> %6.1f %7s %5s %6.3f',char(eids{ii}),...
- Tvec(tt),char(strec{jj}),char(comps{kk}),ratio_sort(hh,5) ));
- end
-
- %------------------------
-
- rec_array = []; ll = 0;
- for tt = 1:nper
- for ii = imin:imax
- disp(sprintf('Event %i out of %i',ii,nevent));
- for jj = 1:nrec
-
- % check that there is a measurement for all components for all models
- if prod( wdiff_array(ii,jj,:,tt,:) ) > 0
- ll = ll+1;
- rat = prod( wdiff_array(ii,jj,:,tt,imodel_max+1) ./ wdiff_array(ii,jj,:,tt,imodel_min+1) );
- rec_array(ll,:) = [tt ii jj rat];
- end
- end
-
- end
- end
- rec_sort = sortrows(rec_array,[-1 4]);
-
- for hh = 1:length(rec_sort)
- tt = rec_sort(hh,1);
- ii = rec_sort(hh,2);
- jj = rec_sort(hh,3);
- disp(sprintf('%10s --> %6.1f %7s %5s %10.6f',char(eids{ii}),...
- Tvec(tt),char(strec{jj}),rec_sort(hh,4) ));
- end
-
- %------------------------
-
-% % loop over all stations
-% tarray = []; ll = 0;
-% for jj = 1:nrec
-% imatch = find( ratio_array(:,3)==jj);
-% if length(imatch)==6
-% ll = ll+1;
-% tarray(ll,:) = [jj prod(ratio_array(imatch,5)) ];
-% end
-% end
-% tarray_sort = sortrows(tarray,2);
-% for hh = 1:20
-% jj = tarray_sort(hh,1);
-% disp(sprintf('%7s %6.3f',char(strec{jj}),tarray_sort(hh,2) ));
-% end
-
- figure; N = plot_histo(ratio_array(:,end),[0:0.2:5]);
- figure; N = plot_histo(log(ratio_array(:,end)),[-4:0.2:4]);
-
- break
-end
-
-%=======================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2008-10-21 01:21:15 UTC (rev 13106)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -25,17 +25,22 @@
dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
-Tvec = [6 2];
- % LOWEST period used for each measurement band-pass (upper is 40 s)
+% min and max periods for the different bandpassed datasets
+Tminvec = [2 6]; Tmaxvec = [30 30];
+%Tminvec = [2 3 6]; Tmaxvec = [30 30 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(Tvec);
+nper = length(Tminvec);
-%DT = 0.011; % really this is only needed because it was left out from
- % integrating the waveforms in mt-measure_adj.f90
-
+% strings for labeling
+sTbp = repmat(cellstr(' '),1,nper);
+for tt = 1:nper
+ sTbp{tt} = sprintf('[%is,%is]',Tminvec(tt),Tmaxvec(tt));
+end
+
%-------------------------
% USER INPUT
@@ -53,16 +58,16 @@
stdcov = stdcovs{idatacov};
ftag = [stmod '_' stdcov];
-odir = ['OUTPUT/' stmod '/' stdcov '/'];
+odir = ['OUTPUT_SUBSPACE/' stmod '/' stdcov '/'];
iwrite = input(' Enter 1 to write files (0 otherwise) : ');
-% files: event IDs, CMTSOLUTIONS, ststions
-dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_9/';
+% files: event IDs, CMTSOLUTIONS, stations
+dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_10/';
eid_file = [dir_source 'EIDs_only_loc'];
-eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_' stmod];
-%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
-cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v9'];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_' stmod];
+eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v10'];
stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
%-------------------------
@@ -93,7 +98,7 @@
% if 0==1
% fid = fopen([odir ftag '_window_picks_by_' stks{kk}],'w');
% fprintf(fid,'%12s%5s%10s%10s%10s\n',stks{kk},' ',...
-% ['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
+% ['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);
@@ -162,16 +167,17 @@
meas_array = [];
%for tt = 1:1
for tt = 1:nper
- Tper = Tvec(tt);
- sTper = sprintf('%3.3i',Tper);
- if imod <= 2, sTper2 = sprintf('%2.2i',Tper); else sTper2 = sTper; end % CHT
+ Tmin = Tminvec(tt);
+ Tmax = Tmaxvec(tt);
+ Ttag = sprintf('T%3.3i_T%3.3i',Tmin,Tmax);
+ %Ttag = sprintf('T%3.3i',Tmin); % old version
disp('-----------------------------------------');
for ii = imin:imax % loop over events
disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
- dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
- mfile = [eids{ii} '_T' sTper2 '_window_chi'];
+ dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_' Ttag '/'];
+ mfile = [eids{ii} '_' Ttag '_window_chi'];
filename = [dir1 mfile];
if ~exist(filename,'file')
@@ -260,7 +266,7 @@
disp('====>');
itest = itestvec(ii);
meas_array(itest,:)
- display_meas(meas_array(itest,:),Tvec,eids,strec,comps);
+ display_meas(meas_array(itest,:),Tminvec,eids,strec,comps);
meas_array(itest,14)
0.5 * meas_array(itest,9)^2 / meas_array(itest,10)^2
0.5 * meas_array(itest,13)^2 / meas_array(itest,10)^2
@@ -273,28 +279,28 @@
% (2) network
% (3) station
% (4) component
-nwin_all_event = zeros(nevent,2);
+nwin_all_event = zeros(nevent,nper);
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);
end
end
-nwin_all_net = zeros(nnet,2);
+nwin_all_net = zeros(nnet,nper);
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);
+nwin_all_rec = zeros(nrec,nper);
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);
+nwin_all_comp = zeros(ncomp,nper);
for tt = 1:nper
for ii = 1:ncomp
imatch = find( and( meas_array(:,6)==ii, meas_array(:,2)==tt) );
@@ -303,18 +309,24 @@
end
%----------------------------------------------
-% tally stations by event (and period)
-nrec_all_event = zeros(nevent,3);
+% tally stations for each event (and period)
+nrec_all_event = zeros(nevent,nper+1);
for ii = 1:nevent
- imatch1 = find( and( meas_array(:,3)==ii, meas_array(:,2)==1) );
- imatch2 = find( and( meas_array(:,3)==ii, meas_array(:,2)==2) );
+ % 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,1) = length(unique(meas_array(imatch1,5)));
- nrec_all_event(ii,2) = length(unique(meas_array(imatch2,5)));
- nrec_all_event(ii,3) = length(unique(meas_array(imatch,5)));
+ nrec_all_event(ii,nper+1) = length(unique(meas_array(imatch,5)));
end
if iwrite == 1
+ sTfmti = repmat('%10i',1,nper);
+ sTfmts = repmat('%10s',1,nper);
+ sTdash = repmat(cellstr('--'),1,nper);
+
stks = {'event','network','receiver','component'};
for kk = 1:length(stks)
switch kk
@@ -326,29 +338,27 @@
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);
+ [junk, isort] = sortrows([[1:nw]' nwin_out ntot_2],-[2+nper]);
% write out the file
fid = fopen([odir ftag '_window_picks_by_' stks{kk}],'w');
- fprintf(fid,'%12s%5s%10s%10s%10s\n',stks{kk},' ',...
- ['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
- fprintf(fid,'%12s%5s%10i%10i%10i\n','TOTAL -->',' ',ntot_1,sum(ntot_1));
+ fprintf(fid,['%12s%5s' sTfmts '%10s\n'],stks{kk},' ',sTbp{:},'TOTAL');
+ fprintf(fid,['%12s%5s' sTfmti '%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));
+ fprintf(fid,['%12s%5i' sTfmti '%10i\n'],labs{jj},jj,nwin_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],-4);
- fid = fopen([odir ftag '_receiver_picks_by_event'],'w');
- fprintf(fid,'%12s%5s%10s%10s%10s\n','event',' ',...
- ['Tmin=' num2str(Tvec(1))],['Tmin=' num2str(Tvec(2))],'TOTAL');
+ [junk, isort] = sortrows([[1:nevent]' nrec_all_event],-[2+nper]);
+ fid = fopen([odir ftag '_receivers_used_by_event'],'w');
+ fprintf(fid,['%12s%5s' sTfmts '%10s\n'],'event',' ',sTbp{:},'TOTAL');
for ii = 1:nevent
jj = isort(ii);
- fprintf(fid,'%12s%5i%10i%10i%10i\n',eids{jj},jj,nrec_all_event(jj,:));
+ fprintf(fid,['%12s%5i' sTfmti '%10i\n'],eids{jj},jj,nrec_all_event(jj,:));
end
fclose(fid);
end
@@ -417,16 +427,17 @@
end
klabs = {'dnorm','nrec'};
- kind = [-6 -5];
- for kk=1:2
+ kind = -([4 3]+nper); % columns to sort
+ for kk = 1:2
dsort = sortrows(data_norms,kind(kk));
fid = fopen([odir ftag '_data_norms_sort_by_' klabs{kk} '.txt'],'w');
- fprintf(fid,'%12s%8s%6s%6s%6s%12s%12s%12s%12s\n','eid','Nwin',...
- ['Nrec' num2str(Tvec(1))],['Nrec' num2str(Tvec(2))],'Nrec','dnorm2*E','dnorm2','dnorm','weight');
- fprintf(fid,'%12s%8i%6s%6s%6s%12.2f%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq*nevent),sum(dnorm_sq),'--','--');
+ fprintf(fid,['%12s%8s' sTfmts '%10s%10s%10s%10s%10s\n'],'eid','Nwin',...
+ sTbp{:},'Nrec','dnorm2*E','dnorm2','dnorm','weight');
+ fprintf(fid,['%12s%8i' sTfmts '%10s%10.2f%10.4f%10s%10s\n'],...
+ 'TOTAL',sum(Ns),sTdash{:},'--',sum(dnorm_sq*nevent),sum(dnorm_sq),'--','--');
for ii = 1:nevent
if any(dsort(ii,1)==eids_save)
- fprintf(fid,'%12i%8i%6i%6i%6i%12.4f%12.4f%12.4f%12.4f\n',dsort(ii,:));
+ fprintf(fid,['%12i%8i' sTfmti '%10i%10.4f%10.4f%10.4f%10.4f\n'],dsort(ii,:));
end
end
%fprintf(fid,'%12s%8i%6s%6s%6s%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq),'--','--');
@@ -488,7 +499,7 @@
icheck = find( abs(meas_array(:,9)) >= DT_MIN )
[junk, isort] = sortrows( meas_array(icheck,:), -9 );
meas_disp = meas_array(icheck(isort),:);
-display_meas(meas_disp,Tvec,eids,strec,comps);
+display_meas(meas_disp,Tminvec,eids,strec,comps);
DA_MIN = -0.8;
DA_MAX = 3.0;
@@ -496,14 +507,14 @@
%[junk, isort] = sortrows( meas_array(icheck,:), 11 );
[junk, isort] = sortrows( meas_array(icheck,:), 3 );
meas_disp = meas_array(icheck(isort),:);
-display_meas(meas_disp,Tvec,eids,strec,comps);
+display_meas(meas_disp,Tminvec,eids,strec,comps);
CHI_MIN = 100;
DT_SIGMA_MIN = 0.19;
icheck = find( and( meas_array(:,17) >= CHI_MIN, meas_array(:,10) < DT_SIGMA_MIN) );
[junk, isort] = sortrows( meas_array(icheck,:), -15 );
meas_disp = meas_array(icheck(isort),:);
-display_meas(meas_disp,Tvec,eids,strec,comps);
+display_meas(meas_disp,Tminvec,eids,strec,comps);
ebads = unique( eids_num(meas_array(icheck,3)) )
%=======================================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_misfit_dir.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_misfit_dir.pl 2008-10-21 01:21:15 UTC (rev 13106)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_misfit_dir.pl 2008-10-21 02:46:27 UTC (rev 13107)
@@ -1,21 +0,0 @@
-#!/usr/bin/perl -w
-
-#-----------------------------------
-# Carl Tape, 03-Sept-2008
-# setup_misfit_dir.pl
-#
-# This script xxx
-#
-# EXAMPLE: setup_misfit_dir.pl m00 m07 2/40
-#
-#-----------------------------------
-
-if (@ARGV < 3) {die("Usage: setup_misfit_dir.pl smodel_min smodel_max Tmin/Tmax\n")}
-($minm,$maxm,$Tpers) = @ARGV;
-
-$dir0 = "OUTPUT_MISFIT";
-`mkdir -p $dir0`;
-
-# TO DO --->
-
-#================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_subspace_dir.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_subspace_dir.pl 2008-10-21 01:21:15 UTC (rev 13106)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/setup_subspace_dir.pl 2008-10-21 02:46:27 UTC (rev 13107)
@@ -14,12 +14,14 @@
#
# This should be run prior to running subspace_specfem.m.
#
-# EXAMPLE: setup_subspace_dir.pl 10
+# EXAMPLE:
+# setup_subspace_dir.pl 0 10
+# setup_subspace_dir.pl 11 11
#
#-----------------------------------
-if (@ARGV < 1) {die("Usage: setup_subspace_dir.pl smodel\n")}
-($maxm) = @ARGV;
+if (@ARGV < 2) {die("Usage: setup_subspace_dir.pl smodel_min smodel_max\n")}
+($minm,$maxm) = @ARGV;
# possible labels for choice of data covariance
@dcov_tags = ("event","window","none");
@@ -37,7 +39,7 @@
#$dir1 = "$dir0/$smodel";
#`mkdir -p $dir1`;
-for ($h = 0; $h <= $maxm; $h = $h+1) {
+for ($h = $minm; $h <= $maxm; $h = $h+1) {
$smodel = sprintf("m%2.2i",$h);
$dir1 = "${dir0}/${smodel}";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2008-10-21 01:21:15 UTC (rev 13106)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2008-10-21 02:46:27 UTC (rev 13107)
@@ -1,6 +1,6 @@
%
% subspace_specfem.m
-% CARL TAPE, 16-June-2008
+% CARL TAPE, 16-Oct-2008
% printed xxx
%
% This program assumes you have FIRST run compute_misfit.m, and also have
@@ -38,7 +38,7 @@
stdcovs = {'event','window','none'};
stdcov = stdcovs{idatacov};
ftag = [stmod '_' stdcov];
-odir0 = ['OUTPUT/' stmod '/' stdcov '/'];
+odir0 = ['OUTPUT_SUBSPACE/' stmod '/' stdcov '/'];
load([odir0 ftag '_data_norms']);
nsrc = length(eids);
@@ -130,7 +130,7 @@
% 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.
+ % 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
@@ -152,6 +152,13 @@
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');
@@ -170,6 +177,28 @@
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) :');
@@ -196,12 +225,13 @@
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([0 1e-10])
- map = colormap('gray'); colormap(flipud(map));
- colorbar; axis equal; axis tight;
+ 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
More information about the CIG-COMMITS
mailing list