profilesMatlabProcessing/ecco_v4/insitu_cost.m

function []=insitu_cost(doComp,varargin);
%object:        computes or displays cost function statistics
%input:         doComp states whether to compute & save (1) or load & display (0)
%optional must take the following form {'name',param1,param2,...} Active ones are
%               {'dirData',dirData} is the data subdirectory name ('profiles/output/' by default)
%               {'listData',data1,data2,...} where data1 e.g. is 'argo_indian.nc'
%               {'listVar',var1,var2} where var1 and var2 are 'T' and 'S'
%               {'years',year0,year1} to limit computation to [year0 year1] interval
%               {'dirMat',dirMat} where dirMat is the output subdirectory ('mat/' by default)
%               {'suffMat',suffMat} is the mat file suffix ('all' by default)
%               {'dirTex',dirTex} where dirTex is the tex directory ('tex/' by default)
%               {'nameTex',nameTex} is the tex file name ('myPlots' by default)
%               {'addToTex',addToTex} states whether (1) or not (0 default) to 
%                        augment an ongoing tex file (see write2tex.m, dirTex)
%
%example:         insitu_cost('./',1,{'listData','argo_in*'});

gcmfaces_global; 
global myparms; 

listVar={'T','S'}; 
listBas={'atlExt','pacExt','indExt','arct'};
if sum([90 1170]~=mygrid.ioSize)>0; listBas={}; end;
listBasTxt=''; for bb=1:length(listBas); listBasTxt=[listBasTxt ' ' listBas{bb} '*']; end;
%shorter test case: listData={'argo_in*'}; listVar={'T'};
dirData='/profiles/output/';
dirMat='mat/'; suffMat='_all';
dirTex='tex/'; addToTex=0; nameTex='myPlots';
%set more optional paramaters to user defined values
for ii=1:nargin-1;
   if strcmp(varargin{ii}{1},'listData'); listData={varargin{ii}{2:end}}; 
   elseif strcmp(varargin{ii}{1},'listVar'); listVar={varargin{ii}{2:end}};
   elseif strcmp(varargin{ii}{1},'years'); year0=varargin{ii}{2}; year1=varargin{ii}{3};
   elseif strcmp(varargin{ii}{1},'dirMat'); dirMat=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];
   elseif strcmp(varargin{ii}{1},'dirTex'); dirTex=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'nameTex'); nameTex=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'dirData'); dirData=varargin{ii}{2};
   else;
       warning('inputCheck:insitu_cost',...
           ['unknown option ''' varargin{ii}{1} ''' was ignored']);
   end;
end;

if isempty(myparms)|isempty(mygrid);
  load([dirMat 'diags_grid_parms.mat']);
end;

if isempty(whos('year0'));
  year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);
end;

if isempty(whos('listData'));
  listData=dir([dirMat 'profiles/output/*.nc'])
  listData={listData(:).name};
  for ff=1:length(listData); listData{ff}=[listData{ff}(1:end-3) '*']; end;
end;

if doComp;

%time limits:
date0=datenum(year0,1,1); date1=datenum(year1,12,31);

for vv=1:length(listVar);
    varCur=listVar{vv};
   
    %get the cost square root:
    [MITprof]=MITprof_stats_load([dirData],listData,varCur);
    if ~isfield(MITprof,'prof_point');
      loc_tile=gcmfaces_loc_tile(90,90,MITprof.prof_lon,MITprof.prof_lat);
      MITprof.prof_point=loc_tile.point;
    end;
    %mask out values that are not in year range:
    ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);
    MITprof.prof(ii,:)=NaN;
    
    %mean and median:
    costCur=[nanmean(MITprof.prof(:).^2) nanmedian(MITprof.prof(:).^2)];
    fprintf('mean cost for %s: %0.3g (median: %0.3g) \n',varCur,costCur);
    eval(['costAve' varCur '.mean=costCur(1);']);
    eval(['costAve' varCur '.median=costCur(2);']);
    
    %depth/time mean cost:
    tmp_date=(MITprof.prof_date-date0)/365;
%    [x,y,z,n]=MITprof_stats(MITprof.prof_depth,MITprof.prof_depth,...
    [x,y,z,n]=MITprof_stats(MITprof.prof_depth,[0:200:1000 1500:500:6000],...
        tmp_date,[0:1/4:1+(year1-year0)],'mean',MITprof.prof.^2);
    z(n==0)=NaN; n(n==0)=NaN;
    eval(['depthTimeCost' varCur '.x=x;']); eval(['depthTimeCost' varCur '.y=y;']);
    eval(['depthTimeCost' varCur '.z=z;']); eval(['depthTimeCost' varCur '.n=n;']);
    
    %latitudinal distribution:
    [x,y,z,n]=MITprof_stats(MITprof.prof_lat,[-90:5:90],MITprof.prof,[-5:0.25:5]);
    z(n==0)=NaN; n(n==0)=NaN;
    eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);
    eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);

    %upper ocean temporal distribution
    for bb=1:length(listBas); 
    for ll=1:3;
      bbb=listBas{bb};
      if ll==1;     Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
      elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
      elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
      end;
      msk=v4_basin(bbb);
      msk(mygrid.YC<Lmin|mygrid.YC>Lmax)=0;
      %
      msk=convert2array(msk);
      MITprof.prof_msk=msk(MITprof.prof_point);     
      MITprofSub=MITprof_subset(MITprof,'msk',1,'depth',[0 700]);      
      %
      tmp_date=(MITprofSub.prof_date-date0)/365;
      [x,y,z,n]=MITprof_stats(tmp_date,[0:1/4:1+(year1-year0)],MITprofSub.prof,[-5:0.25:5]);
      z(n==0)=NaN; n(n==0)=NaN;
      eval([txt '.x=x;']); eval([txt '.y=y;']);
      eval([txt '.z=z;']); eval([txt '.n=n;']);
    end;
    end;

end;

if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat 'cost/insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* ' ...
      'listData listVar listBas ' listBasTxt ';']);

else;%display result

if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;

eval(['load ' dirMat 'insitu_cost' suffMat '.mat;']);

figureL;

for vv=1:length(listVar);
    varCur=listVar{vv};

    eval(['costCur=[costAve' varCur '.mean costAve' varCur '.median];']);

    eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);
    eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);
    z(n<1e2)=NaN; y=y+year0-1;
    subplot(2,length(listVar),vv); depthStretchPlot('pcolor',{y,x,z},[0:200:1000 1500:500:6000],[0 1000 4000]);
    caxis([0 8]); shading flat; colorbar; ylabel('depth (in m)');
    title(sprintf('mean cost for %s: %0.3g (median: %0.3g)',varCur,costCur));

    %latitudinal distribution:
    eval(['x=misfitDistrib' varCur '.x;']); eval(['y=misfitDistrib' varCur '.y;']);
    eval(['z=misfitDistrib' varCur '.z;']); eval(['n=misfitDistrib' varCur '.n;']);
    tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
    subplot(2,length(listVar),vv+length(listVar)); contourf(x,y,z,[0:0.05:0.5]); grid on;
    caxis([0 0.5]); colorbar; ylabel('normalized misfit'); xlabel('latitude'); title('pdf');

end;

myCaption={'Cost function (top) for in situ profiles, as a function of depth and time. ',...
        'Distribution of normalized misfits (bottom) as a function of latitude. For T (left) and S (right).'};
if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;

if ~isempty(listBas);

ii=0;
ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=3;
ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=3;
ii=ii+1; listPanels(ii).bas='arct'; listPanels(ii).ll=3;
ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=2;
ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=2;
ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=2;
ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=1;
ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=1;
ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=1;

for vv=1:length(listVar);
  varCur=listVar{vv};
  figureL;
  for ii=1:9;
      bbb=listPanels(ii).bas; ll=listPanels(ii).ll;
      if ll==1;     Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
      elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
      elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
      end;

    eval(['x=' txt '.x;']); eval(['y=' txt '.y;']);
    eval(['z=' txt '.z;']); eval(['n=' txt '.n;']);
    x=x+year0-1;
    tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
    subplot(3,3,ii); contourf(x,y,z,[0:0.05:0.5]); grid on;
    caxis([0 0.5]); colorbar; title(txt,'Interpreter','none');
    %title([txt ', normalized misfit pdf'],'Interpreter','none');
    if mod(ii,3)==1; ylabel('normalized misfit'); end;
    %if ii<=6; set(gca,'XTick',[]); end;
  end;

myCaption={'Distribution of normalized misfits per basin (panel) as ',...
        'a function of latitude, for ',varCur};
if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;

end;

end;

end;


1	gforget	1.5	function []=insitu_cost(doComp,varargin);
2	gforget	1.4	%object: computes or displays cost function statistics
3	gforget	1.5	%input: doComp states whether to compute & save (1) or load & display (0)
4	gforget	1.1	%optional must take the following form {'name',param1,param2,...} Active ones are
5	gforget	1.5	% {'dirData',dirData} is the data subdirectory name ('profiles/output/' by default)
6	gforget	1.1	% {'listData',data1,data2,...} where data1 e.g. is 'argo_indian.nc'
7			% {'listVar',var1,var2} where var1 and var2 are 'T' and 'S'
8			% {'years',year0,year1} to limit computation to [year0 year1] interval
9	gforget	1.5	% {'dirMat',dirMat} where dirMat is the output subdirectory ('mat/' by default)
10	gforget	1.4	% {'suffMat',suffMat} is the mat file suffix ('all' by default)
11	gforget	1.5	% {'dirTex',dirTex} where dirTex is the tex directory ('tex/' by default)
12			% {'nameTex',nameTex} is the tex file name ('myPlots' by default)
13	gforget	1.1	% {'addToTex',addToTex} states whether (1) or not (0 default) to
14	gforget	1.5	% augment an ongoing tex file (see write2tex.m, dirTex)
15	gforget	1.1	%
16	gforget	1.4	%example: insitu_cost('./',1,{'listData','argo_in*'});
17	gforget	1.1
18	gforget	1.12	gcmfaces_global;
19			global myparms;
20
21	gforget	1.2	listVar={'T','S'};
22	gforget	1.7	listBas={'atlExt','pacExt','indExt','arct'};
23	gforget	1.9	if sum([90 1170]~=mygrid.ioSize)>0; listBas={}; end;
24	gforget	1.7	listBasTxt=''; for bb=1:length(listBas); listBasTxt=[listBasTxt ' ' listBas{bb} '*']; end;
25	gforget	1.2	%shorter test case: listData={'argo_in*'}; listVar={'T'};
26	gforget	1.5	dirData='/profiles/output/';
27			dirMat='mat/'; suffMat='_all';
28			dirTex='tex/'; addToTex=0; nameTex='myPlots';
29	gforget	1.1	%set more optional paramaters to user defined values
30	gforget	1.5	for ii=1:nargin-1;
31	gforget	1.1	if strcmp(varargin{ii}{1},'listData'); listData={varargin{ii}{2:end}};
32			elseif strcmp(varargin{ii}{1},'listVar'); listVar={varargin{ii}{2:end}};
33			elseif strcmp(varargin{ii}{1},'years'); year0=varargin{ii}{2}; year1=varargin{ii}{3};
34	gforget	1.5	elseif strcmp(varargin{ii}{1},'dirMat'); dirMat=varargin{ii}{2};
35	gforget	1.1	elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];
36	gforget	1.5	elseif strcmp(varargin{ii}{1},'dirTex'); dirTex=varargin{ii}{2};
37	gforget	1.1	elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};
38	gforget	1.5	elseif strcmp(varargin{ii}{1},'nameTex'); nameTex=varargin{ii}{2};
39			elseif strcmp(varargin{ii}{1},'dirData'); dirData=varargin{ii}{2};
40	gforget	1.1	else;
41			warning('inputCheck:insitu_cost',...
42			['unknown option ''' varargin{ii}{1} ''' was ignored']);
43			end;
44			end;
45
46	gforget	1.12	if isempty(myparms)\|isempty(mygrid);
47			load([dirMat 'diags_grid_parms.mat']);
48			end;
49
50			if isempty(whos('year0'));
51			year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);
52			end;
53
54			if isempty(whos('listData'));
55			listData=dir([dirMat 'profiles/output/*.nc'])
56			listData={listData(:).name};
57			for ff=1:length(listData); listData{ff}=[listData{ff}(1:end-3) '*']; end;
58			end;
59
60	gforget	1.1	if doComp;
61
62			%time limits:
63			date0=datenum(year0,1,1); date1=datenum(year1,12,31);
64
65			for vv=1:length(listVar);
66			varCur=listVar{vv};
67	gforget	1.13
68	gforget	1.1	%get the cost square root:
69	gforget	1.5	[MITprof]=MITprof_stats_load([dirData],listData,varCur);
70	gforget	1.13	if ~isfield(MITprof,'prof_point');
71			loc_tile=gcmfaces_loc_tile(90,90,MITprof.prof_lon,MITprof.prof_lat);
72			MITprof.prof_point=loc_tile.point;
73			end;
74	gforget	1.1	%mask out values that are not in year range:
75			ii=find(MITprof.prof_date<date0\|MITprof.prof_date>date1);
76			MITprof.prof(ii,:)=NaN;
77
78			%mean and median:
79			costCur=[nanmean(MITprof.prof(:).^2) nanmedian(MITprof.prof(:).^2)];
80			fprintf('mean cost for %s: %0.3g (median: %0.3g) \n',varCur,costCur);
81			eval(['costAve' varCur '.mean=costCur(1);']);
82			eval(['costAve' varCur '.median=costCur(2);']);
83
84			%depth/time mean cost:
85			tmp_date=(MITprof.prof_date-date0)/365;
86			% [x,y,z,n]=MITprof_stats(MITprof.prof_depth,MITprof.prof_depth,...
87			[x,y,z,n]=MITprof_stats(MITprof.prof_depth,[0:200:1000 1500:500:6000],...
88			tmp_date,[0:1/4:1+(year1-year0)],'mean',MITprof.prof.^2);
89			z(n==0)=NaN; n(n==0)=NaN;
90			eval(['depthTimeCost' varCur '.x=x;']); eval(['depthTimeCost' varCur '.y=y;']);
91			eval(['depthTimeCost' varCur '.z=z;']); eval(['depthTimeCost' varCur '.n=n;']);
92
93			%latitudinal distribution:
94			[x,y,z,n]=MITprof_stats(MITprof.prof_lat,[-90:5:90],MITprof.prof,[-5:0.25:5]);
95			z(n==0)=NaN; n(n==0)=NaN;
96			eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);
97			eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);
98
99	gforget	1.7	%upper ocean temporal distribution
100			for bb=1:length(listBas);
101			for ll=1:3;
102			bbb=listBas{bb};
103			if ll==1; Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
104			elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
105			elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
106			end;
107			msk=v4_basin(bbb);
108			msk(mygrid.YC<Lmin\|mygrid.YC>Lmax)=0;
109			%
110			msk=convert2array(msk);
111			MITprof.prof_msk=msk(MITprof.prof_point);
112			MITprofSub=MITprof_subset(MITprof,'msk',1,'depth',[0 700]);
113			%
114			tmp_date=(MITprofSub.prof_date-date0)/365;
115			[x,y,z,n]=MITprof_stats(tmp_date,[0:1/4:1+(year1-year0)],MITprofSub.prof,[-5:0.25:5]);
116			z(n==0)=NaN; n(n==0)=NaN;
117			eval([txt '.x=x;']); eval([txt '.y=y;']);
118			eval([txt '.z=z;']); eval([txt '.n=n;']);
119			end;
120			end;
121
122	gforget	1.1	end;
123
124	gforget	1.10	if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
125			eval(['save ' dirMat 'cost/insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* ' ...
126	gforget	1.7	'listData listVar listBas ' listBasTxt ';']);
127	gforget	1.1
128			else;%display result
129
130	gforget	1.9	if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
131
132	gforget	1.5	eval(['load ' dirMat 'insitu_cost' suffMat '.mat;']);
133	gforget	1.1
134			figureL;
135
136			for vv=1:length(listVar);
137			varCur=listVar{vv};
138
139			eval(['costCur=[costAve' varCur '.mean costAve' varCur '.median];']);
140
141			eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);
142			eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);
143	gforget	1.7	z(n<1e2)=NaN; y=y+year0-1;
144	gforget	1.1	subplot(2,length(listVar),vv); depthStretchPlot('pcolor',{y,x,z},[0:200:1000 1500:500:6000],[0 1000 4000]);
145			caxis([0 8]); shading flat; colorbar; ylabel('depth (in m)');
146			title(sprintf('mean cost for %s: %0.3g (median: %0.3g)',varCur,costCur));
147
148			%latitudinal distribution:
149			eval(['x=misfitDistrib' varCur '.x;']); eval(['y=misfitDistrib' varCur '.y;']);
150			eval(['z=misfitDistrib' varCur '.z;']); eval(['n=misfitDistrib' varCur '.n;']);
151			tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
152			subplot(2,length(listVar),vv+length(listVar)); contourf(x,y,z,[0:0.05:0.5]); grid on;
153			caxis([0 0.5]); colorbar; ylabel('normalized misfit'); xlabel('latitude'); title('pdf');
154
155			end;
156
157	gforget	1.3	myCaption={'Cost function (top) for in situ profiles, as a function of depth and time. ',...
158	gforget	1.4	'Distribution of normalized misfits (bottom) as a function of latitude. For T (left) and S (right).'};
159	gforget	1.5	if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
160	gforget	1.1
161	gforget	1.8	if ~isempty(listBas);
162
163	gforget	1.7	ii=0;
164			ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=3;
165			ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=3;
166			ii=ii+1; listPanels(ii).bas='arct'; listPanels(ii).ll=3;
167			ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=2;
168			ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=2;
169			ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=2;
170			ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=1;
171			ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=1;
172			ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=1;
173
174			for vv=1:length(listVar);
175			varCur=listVar{vv};
176			figureL;
177			for ii=1:9;
178			bbb=listPanels(ii).bas; ll=listPanels(ii).ll;
179			if ll==1; Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
180			elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
181			elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
182			end;
183
184			eval(['x=' txt '.x;']); eval(['y=' txt '.y;']);
185			eval(['z=' txt '.z;']); eval(['n=' txt '.n;']);
186			x=x+year0-1;
187			tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
188			subplot(3,3,ii); contourf(x,y,z,[0:0.05:0.5]); grid on;
189			caxis([0 0.5]); colorbar; title(txt,'Interpreter','none');
190			%title([txt ', normalized misfit pdf'],'Interpreter','none');
191			if mod(ii,3)==1; ylabel('normalized misfit'); end;
192			%if ii<=6; set(gca,'XTick',[]); end;
193			end;
194
195			myCaption={'Distribution of normalized misfits per basin (panel) as ',...
196			'a function of latitude, for ',varCur};
197			if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
198
199			end;
200
201	gforget	1.1	end;
202
203	gforget	1.8	end;
204
205	gforget	1.1