/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/ecco_v4/insitu_cost.m
ViewVC logotype

Contents of /MITgcm_contrib/gael/profilesMatlabProcessing/ecco_v4/insitu_cost.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.5 - (show annotations) (download)
Sun Dec 30 00:35:53 2012 UTC (12 years, 6 months ago) by gforget
Branch: MAIN
Changes since 1.4: +22 -21 lines
- revise handling of directories.
- update listData.
- add nameTex optional input parameter.

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

  ViewVC Help
Powered by ViewVC 1.1.22