profilesMatlabProcessing/ecco_v4/insitu_cost.m

function []=insitu_cost(dirData,doComp,varargin);
%[]=insitu_cost(dirData,listData,listVar,varargin);
%object:        computes or displays cost function statistics
%input:         dirData is the data main directory name
%               doComp states whether to compute & save (1) or load & display (0)
%optional must take the following form {'name',param1,param2,...} Active ones are
%               {'subdirData',subdirData} 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
%               {'subdirMat',subdirMat} where subdirMat is the output subdirectory ('mat/' by default)
%               {'suffMat',suffMat} is the mat file suffix ('all' by default)
%               {'subdirTex',subdirTex} where subdirTex is the tex directory ('tex/' by default)
%               {'addToTex',addToTex} states whether (1) or not (0 default) to 
%                        augment an ongoing tex file (see write2tex.m, subdirTex)
%
%example:         insitu_cost('./',1,{'listData','argo_in*'});

%load paths and grid:
%--------------------
global mygrid; if isempty(mygrid); grid_load('GRID/',5,'compact'); end;

%set more optional paramaters to default values
global myparms; year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);
%year0=1992; year1=2008; 
listData={'WOD09_XBT*','WOD09_CTD*','argo_in*','argo_at*','argo_pa*','seals*'}; 
listVar={'T','S'}; 
%shorter test case: listData={'argo_in*'}; listVar={'T'};
subdirData='/profiles/output/';
subdirMat='mat/'; suffMat='_all';
subdirTex='tex/'; addToTex=0;
%set more optional paramaters to user defined values
for ii=1:nargin-2;
   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},'subdirMat'); subdirMat=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];
   elseif strcmp(varargin{ii}{1},'subdirTex'); subdirTex=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};
   elseif strcmp(varargin{ii}{1},'subdirData'); subdirData=varargin{ii}{2};
   else;
       warning('inputCheck:insitu_cost',...
           ['unknown option ''' varargin{ii}{1} ''' was ignored']);
   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 '/' subdirData],listData,varCur);
    
    %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;']);

end;

eval(['save ' dirData '/' subdirMat 'insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* listData listVar;']);

else;%display result

eval(['load ' dirData '/' subdirMat '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;
    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)');
    xlabel(['date (in years since ' num2str(year0) ')']);
    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([dirData '/' subdirTex '/myPlots.tex'],2,myCaption,gcf); end;

end;


1	function []=insitu_cost(dirData,doComp,varargin);
2	%[]=insitu_cost(dirData,listData,listVar,varargin);
3	%object: computes or displays cost function statistics
4	%input: dirData is the data main directory name
5	% doComp states whether to compute & save (1) or load & display (0)
6	%optional must take the following form {'name',param1,param2,...} Active ones are
7	% {'subdirData',subdirData} is the data subdirectory name ('profiles/output/' by default)
8	% {'listData',data1,data2,...} where data1 e.g. is 'argo_indian.nc'
9	% {'listVar',var1,var2} where var1 and var2 are 'T' and 'S'
10	% {'years',year0,year1} to limit computation to [year0 year1] interval
11	% {'subdirMat',subdirMat} where subdirMat is the output subdirectory ('mat/' by default)
12	% {'suffMat',suffMat} is the mat file suffix ('all' by default)
13	% {'subdirTex',subdirTex} where subdirTex is the tex directory ('tex/' by default)
14	% {'addToTex',addToTex} states whether (1) or not (0 default) to
15	% augment an ongoing tex file (see write2tex.m, subdirTex)
16	%
17	%example: insitu_cost('./',1,{'listData','argo_in*'});
18
19	%load paths and grid:
20	%--------------------
21	global mygrid; if isempty(mygrid); grid_load('GRID/',5,'compact'); end;
22
23	%set more optional paramaters to default values
24	global myparms; year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);
25	%year0=1992; year1=2008;
26	listData={'WOD09_XBT','WOD09_CTD','argo_in','argo_at','argo_pa','seals'};
27	listVar={'T','S'};
28	%shorter test case: listData={'argo_in*'}; listVar={'T'};
29	subdirData='/profiles/output/';
30	subdirMat='mat/'; suffMat='_all';
31	subdirTex='tex/'; addToTex=0;
32	%set more optional paramaters to user defined values
33	for ii=1:nargin-2;
34	if strcmp(varargin{ii}{1},'listData'); listData={varargin{ii}{2:end}};
35	elseif strcmp(varargin{ii}{1},'listVar'); listVar={varargin{ii}{2:end}};
36	elseif strcmp(varargin{ii}{1},'years'); year0=varargin{ii}{2}; year1=varargin{ii}{3};
37	elseif strcmp(varargin{ii}{1},'subdirMat'); subdirMat=varargin{ii}{2};
38	elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];
39	elseif strcmp(varargin{ii}{1},'subdirTex'); subdirTex=varargin{ii}{2};
40	elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};
41	elseif strcmp(varargin{ii}{1},'subdirData'); subdirData=varargin{ii}{2};
42	else;
43	warning('inputCheck:insitu_cost',...
44	['unknown option ''' varargin{ii}{1} ''' was ignored']);
45	end;
46	end;
47
48	if doComp;
49
50	%time limits:
51	date0=datenum(year0,1,1); date1=datenum(year1,12,31);
52
53	for vv=1:length(listVar);
54	varCur=listVar{vv};
55
56	%get the cost square root:
57	[MITprof]=MITprof_stats_load([dirData '/' subdirData],listData,varCur);
58
59	%mask out values that are not in year range:
60	ii=find(MITprof.prof_date<date0\|MITprof.prof_date>date1);
61	MITprof.prof(ii,:)=NaN;
62
63	%mean and median:
64	costCur=[nanmean(MITprof.prof(:).^2) nanmedian(MITprof.prof(:).^2)];
65	fprintf('mean cost for %s: %0.3g (median: %0.3g) \n',varCur,costCur);
66	eval(['costAve' varCur '.mean=costCur(1);']);
67	eval(['costAve' varCur '.median=costCur(2);']);
68
69	%depth/time mean cost:
70	tmp_date=(MITprof.prof_date-date0)/365;
71	% [x,y,z,n]=MITprof_stats(MITprof.prof_depth,MITprof.prof_depth,...
72	[x,y,z,n]=MITprof_stats(MITprof.prof_depth,[0:200:1000 1500:500:6000],...
73	tmp_date,[0:1/4:1+(year1-year0)],'mean',MITprof.prof.^2);
74	z(n==0)=NaN; n(n==0)=NaN;
75	eval(['depthTimeCost' varCur '.x=x;']); eval(['depthTimeCost' varCur '.y=y;']);
76	eval(['depthTimeCost' varCur '.z=z;']); eval(['depthTimeCost' varCur '.n=n;']);
77
78	%latitudinal distribution:
79	[x,y,z,n]=MITprof_stats(MITprof.prof_lat,[-90:5:90],MITprof.prof,[-5:0.25:5]);
80	z(n==0)=NaN; n(n==0)=NaN;
81	eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);
82	eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);
83
84	end;
85
86	eval(['save ' dirData '/' subdirMat 'insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* listData listVar;']);
87
88	else;%display result
89
90	eval(['load ' dirData '/' subdirMat 'insitu_cost' suffMat '.mat;']);
91
92	figureL;
93
94	for vv=1:length(listVar);
95	varCur=listVar{vv};
96
97	eval(['costCur=[costAve' varCur '.mean costAve' varCur '.median];']);
98
99	eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);
100	eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);
101	z(n<1e2)=NaN;
102	subplot(2,length(listVar),vv); depthStretchPlot('pcolor',{y,x,z},[0:200:1000 1500:500:6000],[0 1000 4000]);
103	caxis([0 8]); shading flat; colorbar; ylabel('depth (in m)');
104	xlabel(['date (in years since ' num2str(year0) ')']);
105	title(sprintf('mean cost for %s: %0.3g (median: %0.3g)',varCur,costCur));
106
107	%latitudinal distribution:
108	eval(['x=misfitDistrib' varCur '.x;']); eval(['y=misfitDistrib' varCur '.y;']);
109	eval(['z=misfitDistrib' varCur '.z;']); eval(['n=misfitDistrib' varCur '.n;']);
110	tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
111	subplot(2,length(listVar),vv+length(listVar)); contourf(x,y,z,[0:0.05:0.5]); grid on;
112	caxis([0 0.5]); colorbar; ylabel('normalized misfit'); xlabel('latitude'); title('pdf');
113
114	end;
115
116	myCaption={'Cost function (top) for in situ profiles, as a function of depth and time. ',...
117	'Distribution of normalized misfits (bottom) as a function of latitude. For T (left) and S (right).'};
118	if addToTex; write2tex([dirData '/' subdirTex '/myPlots.tex'],2,myCaption,gcf); end;
119
120	end;
121
122