Home > atmlab > retrieval > spareice > NNTrainedProduct.m

NNTrainedProduct

PURPOSE ^

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

NNTrainedProduct.m

SOURCE CODE ^

0001 classdef NNTrainedProduct < RetrievalDatabaseProduct
0002     % base class for any ANN-trained product
0003     %
0004     % NNTrainedProduct properties:
0005     %
0006     %   targets - fields to be used as training targets
0007     %
0008     % NNTrainedProduct Methods:
0009     %
0010     %   make_and_train_network -
0011     %   loadnet -
0012     %   storenet -
0013     %   evaluate_test_data -
0014     %   evaluate_and_summarise_testdata -
0015     %
0016     % See also: RetrievalDatabaseProduct, IceMusicNNIWP, CollocatedNNIWP,
0017     % compare_ice_retrievals, compare_collocated_iwp_products,
0018     % InstrumentVisualiser
0019     
0020     % Changes in SPARE-ICE:
0021     %
0022     % version 0.1: fixed bug in storage P_CLOUD: should be float, not int
0023     % version 0.2: fixed bug in lacking filtering for various flagged values
0024     % version 0.3: added surface elevation
0025     % version 0.4: set surface elevation to 0 if z<0
0026     % version 0.5: include some measure of std.dev of z
0027     % version 0.5b: fixed re-introduced bug lacking filtering, new network
0028     % version 0.6: use less info for patternnet than for fitnet
0029     % version 0.7: add corrected meta-info
0030     % version 0.8:
0031    
0032     properties
0033         %selection
0034         
0035         net
0036         trdat
0037                 
0038         % options for plotting scatter-plot
0039         opts = struct('trans', @(x)log10(x), ...
0040            'invtrans', @(x)10.^(x), ...
0041            'axprops', struct(...
0042                 'xscale', 'log', ...
0043                 'yscale', 'log', ...
0044                 'xgrid', 'on', ...
0045                 'xminorgrid', 'off', ...
0046                 'ygrid', 'on', ...
0047                 'yminorgrid', 'off'), ...
0048            'bins', linspace(-1, 4, 50), ...
0049            'medprops', struct('LineWidth', 2, 'LineStyle', '-', 'Color', [0 0 0]), ...
0050            'diagonal', struct('LineWidth', 1, 'LineStyle', ':', 'Color', [0 0 0]))
0051 
0052 %        chanstr;
0053 %        freqstr;
0054 
0055         % hackish way; if an earlier version of the product had what's
0056         % needed, don't need to recollocate
0057         altfilename;
0058         altbasedir;
0059         
0060         % different versions with different network paths.
0061         % may be set to cell array of strings.
0062         networks;        
0063         
0064         % to homog or not to homog, that is the question
0065         homog;
0066     end
0067     
0068     properties (Abstract)
0069         % How many entries to use at most in the training?  If you have
0070         % more data, the rest will go into the testing data.  The data in
0071         % maxsize will be divided between training and validation.
0072         maxsize
0073         
0074         % Structure with one field for each target.
0075         % Must have:
0076         %   - transform     transformation to apply, i.e. log10
0077         %   - invtransform  inverse therefor, i.e. @(x)10.^x
0078         %   - lims          limits to apply before doing anything
0079         % May have:
0080         %   - classify        Also train a classification net with this as
0081         %                     the cloudy/cloudfree cutoff.
0082         %   - bins            Default bins over which to do statistics.
0083         %   - relevant        Minimum x-value to consider in bestnet.
0084         %   - optimise_field  What field (as returned by get_performance)
0085         %                     to optimise for in 'bestnet'.
0086         %
0087         % targets.IWP.transform = @log10
0088         % targets.IWP.invtransform = @(x)10.^(x)
0089         % targets.IWP.lims = [realmin realbax]
0090         % targets.IWP.classify = 10;
0091         targets
0092         
0093         % Structure with one field for each input.
0094         %
0095         %
0096         inputs
0097     end
0098         
0099     methods
0100         function self = NNTrainedProduct(varargin)
0101             self = self@RetrievalDatabaseProduct(varargin{:});
0102                        
0103             self.setmembers();
0104             self.setcols();
0105 %             self.freqstr = sprintf('%s ', self.freqs{self.chans});
0106 %             self.chanstr = sprintf('%d', self.chans);
0107 
0108             if isempty(self.homog)
0109                 self.homog = false;
0110             end
0111         end
0112          
0113         
0114         function make_and_train_network(self, varargin)
0115             % make_and_train_network Create and train NN with self.data
0116             %
0117             % FORMAT
0118             %
0119             %   nniwp.make_and_train_network([noise_level[, count]])
0120             %
0121             % IN
0122             %
0123             %   noise_level     double      Gaussian noise, σ.  If
0124             %                               given, constant for all inputs.
0125             %                               If not given or [], take from
0126             %                               self.inputs.(field).noise.
0127             %
0128             %   count           integer     Number of networks to train.
0129             %                               Defaults to 1.  Can be set
0130             %                               higher to choose e.g. 'best of
0131             %                               5'.
0132             
0133             [noise_level, count] = optargs(varargin, {[], 1});
0134             if isempty(self.data)
0135                 error(['atmlab:' mfilename ':nodata'], ...
0136                     'No data found. Did you run .getdata(...)?');
0137             end
0138             
0139             
0140             %% prepare data, including noise (perhaps)
0141             
0142             logtext(atmlab('OUT'), 'Preparing data\n');
0143             self.limit_data();
0144             dat = self.data; %(self.selection, :);
0145             if self.homog
0146                 dat = self.homog_inputs(dat, self.localcols);
0147             end
0148             [x_all, x_all_cols] = self.select_and_prepare_data(dat, self.localcols, self.inputs);
0149             [y, y_cols] = self.select_and_prepare_data(dat, self.localcols, self.targets);
0150             try
0151                 y_class = self.get_classification(dat);
0152                 classify = true;
0153             catch ME
0154                 switch ME.identifier
0155                     case {['atmlab:' mfilename ':noclassify']}
0156                         classify = false;
0157                     otherwise
0158                         ME.rethrow();
0159                 end
0160             end
0161             % get rid of data ALREADY nan before flagging anything
0162             valid_overall = ~any(isnan(x_all), 1);
0163             % limit further to finite data
0164             % test using imag~=0 because isreal will return a scalar when
0165             % given an array
0166             x_all(imag(x_all)~=0) = NaN;
0167             y(imag(y)~=0) = NaN;
0168             %fin = all(isfinite(x_all), 1) & all(isfinite(y), 1);
0169             fin = all(isfinite(x_all), 1) & ...
0170                   self.select_for_regression(y, y_cols);
0171             all(isfinite(y), 1);
0172             if ~all(fin)
0173                 logtext(atmlab('OUT'), 'Further reducing to %d for finiteness\n', sum(fin));
0174             end
0175             %x = x(:, fin);
0176             %y = y(:, fin);
0177             
0178             % Seed random number generator so I can later obtain the exact
0179             % some noise, should I want to
0180             seed = now();
0181             rng(seed);
0182             if isempty(noise_level) % apply per field
0183                 noise = bsxfun(@times, randn(size(x_all)), self.get_noise_vector(x_all_cols).');
0184             elseif isscalar(noise_level)
0185                 noise = randn(size(x_all)) * noise_level;
0186             else
0187                 noise = bsxfun(@times, randn(size(x_all)), noise_level.');
0188             end
0189             % make sure I keep track of noiseless data, store this along
0190             % with seed
0191             x_orig = x_all;
0192             x_all = x_all + noise;
0193             
0194             % might not use same inputs for both nets
0195             % any field in self.inputs that has self.inputs.(field).net =
0196             % 'fit' will only be used for fitting, and any field that has
0197             % 'pat' will only be used for pattern-matching
0198             [x_fit, x_fit_cols] = self.subselect_inputs(x_all, x_all_cols, self.inputs, 'fit');
0199             [x_orig_fit, ~] = self.subselect_inputs(x_orig, x_all_cols, self.inputs, 'fit');
0200 
0201             if classify
0202                 [x_pat, x_pat_cols] = self.subselect_inputs(x_all, x_all_cols, self.inputs, 'pat');
0203                 [x_orig_pat, ~] = self.subselect_inputs(x_orig, x_all_cols, self.inputs, 'pat');
0204             end
0205            
0206             logtext(atmlab('OUT'), 'Getting %d network(s)\n', count);
0207 %            loaded = false;
0208             for net_count = 0:count-1
0209                 %% Bail out if no need for any network
0210                 
0211                 f = self.calculate_network_file_name(noise_level, net_count);
0212                 if exist(f, 'file') && self.overwrite==0
0213 %                    if loaded
0214 %                        continue
0215 %                    else
0216                     logtext(atmlab('OUT'), 'Already exists, loading from file %s\n', f);
0217                     self.loadnets(f);
0218                     if isa(self.net, 'network') || (classify && ~isfield(self.net, 'pat'))
0219                         logtext(atmlab('OUT'), ...
0220                             'On second thought, I lack a pattern net, retraining!\n');
0221 %                         elseif length(self.net.pat.divideParam.trainInd) + length(self.net.pat.divideParam.valInd) + length(self.net.pat.divideParam.testInd) < length(self.data)
0222 %                             logtext(atmlab('OUT'), ...
0223 %                                 'On second thought, we have different data now, retraining"\n');
0224                     else
0225                         loaded = true;
0226                         continue;
0227                     end
0228 %                    end
0229                 end
0230                 
0231                 %% parallelise
0232                 
0233                 if matlabpool('size') == 0
0234                     matlabpool('open', 6);
0235                 end
0236                 
0237                 %% Set up neural networks
0238                 
0239                 % Create a Fitting and a Pattern Network
0240                 hiddenLayerSize = 10; %min([size(x, 1) 10]);
0241                 logtext(atmlab('OUT'), 'Creating networks to store at %s, hidden layer %d nodes\n', f, hiddenLayerSize);
0242                 fitter = fitnet(hiddenLayerSize);
0243 
0244                 % Choose Input and Output Pre/Post-Processing Functions
0245                 % For a list of all processing functions type: help nnprocess
0246                 %self.net.inputs{1}.processFcns = {'removeconstantrows','mapminmax'};
0247                 %self.net.outputs{2}.processFcns = {'removeconstantrows','mapminmax'};
0248                 fitter.inputs{1}.processFcns = {'removeconstantrows','mapstd'};
0249                 fitter.outputs{2}.processFcns = {'removeconstantrows','mapstd'};
0250                 
0251                 
0252                 % Setup Division of Data for Training, Validation, Testing
0253                 % For a list of all data division functions type: help nndivide
0254                 
0255                 fitter.divideFcn = 'divideind';
0256                 fitter.divideMode = 'sample';  % Divide up every sample
0257 
0258                 % make my own division
0259                 divall = zeros(1, size(y, 2));
0260                 rr = rand(1, size(y, 2));
0261                 trfrac = min([(self.maxsize * (2/3)) / size(y, 2), 0.5]);
0262                 divall(rr <= trfrac) = 1; % train
0263                 divall(rr > trfrac & rr <= trfrac*1.5) = 2; % val
0264                 divall(rr > trfrac * 1.5) = 3; % test (not used in training)
0265                 % but for the fitter, exclude nonfinite values
0266                 fitter.divideParam.trainInd = find(divall==1 & fin);
0267                 fitter.divideParam.valInd = find(divall==2 & fin);
0268                 fitter.divideParam.testInd = find(divall==3 & fin);
0269 %                 fitter.divideParam.trainRatio = 50/100;
0270 
0271                 if classify
0272                     patter = patternnet(hiddenLayerSize);
0273                     patter.divideFcn = 'divideind';
0274                     patter.divideMode = 'sample';
0275                     % for the patter, I should include cloud-free cases
0276                     patter.divideParam.trainInd = find(divall==1 & valid_overall);
0277                     patter.divideParam.valInd = find(divall==2 & valid_overall);
0278                     patter.divideParam.testInd = find(divall==3 & valid_overall);
0279                     patter.trainParam.show = 10;
0280                 end
0281 
0282                 % For help on training function 'trainlm' type: help trainlm
0283                 % For a list of all training functions type: help nntrain
0284                 % See http://www.mathworks.se/help/nnet/ug/speed-and-memory-comparison-for-training-multilayer-networks.html
0285                 fitter.trainFcn = 'trainlm';  % Levenberg-Marquardt
0286                 % choose trainscg, lower on memory.
0287                 % self.net.trainFcn = 'trainscg';
0288                 %            self.net.trainFcn = 'trainbr';  % Bayesian regularisation
0289                 %            self.net.trainFcn = 'trainbfg';  % BFGS quasi-Newton backpropagation
0290                 
0291                 fitter.trainParam.show = 10;
0292                 fitter.efficiency.memoryReduction = 50;
0293                 
0294                 % Choose a Performance Function
0295                 % For a list of all performance functions type: help nnperformance
0296                 fitter.performFcn = 'mse';  % Mean squared error
0297                 
0298                 fitter.userdata.inputs = self.inputs;
0299                 fitter.userdata.targets = self.targets;
0300                 
0301                 % Train the Networks
0302                 
0303                 logtext(atmlab('OUT'), 'Training fitting network (%d inputs, %d hidden, %d target; %d/%d/%d training elements)\n', ...
0304                     size(x_fit, 1), fitter.layers{1}.size, ...
0305                     size(y, 1), length(fitter.divideParam.trainInd), ...
0306                     length(fitter.divideParam.valInd), ...
0307                     length(fitter.divideParam.testInd));
0308                 [fitter, fit_trdat] = train(fitter, double(x_fit), double(y), ...
0309                     'useParallel', 'yes', ...
0310                     'showResources', 'yes', ...
0311                     'reduction', 100);
0312                 
0313                 if classify
0314 
0315                     logtext(atmlab('OUT'), 'Training cloud-filter pattern network (%d inputs, %d hidden, %d target; %d/%d/%d training elements)\n', ...
0316                         size(x_pat, 1), patter.layers{1}.size, ...
0317                         size(y_class, 1), length(patter.divideParam.trainInd), ...
0318                         length(patter.divideParam.valInd), ...
0319                         length(patter.divideParam.testInd));
0320                     
0321                     % for some reason, training a patternnet really doesn't
0322                     % like to see NaNs even in part of the data that aren't
0323                     % used, but set aside for testing or even not used at all.
0324                     % As a workaround, put zeroes here instead, they aren't
0325                     % used anyway...!
0326                     new_y_class = zeros(size(y_class), 'double');
0327                     new_y_class(:, patter.divideParam.trainInd) = y_class(:, patter.divideParam.trainInd);
0328                     new_y_class(:, patter.divideParam.valInd) = y_class(:, patter.divideParam.valInd);
0329                     new_x_pat = zeros(size(x_pat), 'double');
0330                     new_x_pat(:, patter.divideParam.trainInd) = x_pat(:, patter.divideParam.trainInd);
0331                     new_x_pat(:, patter.divideParam.valInd) = x_pat(:, patter.divideParam.valInd);
0332                     
0333                     [patter, pat_trdat] = train(patter, new_x_pat, new_y_class, ...
0334                         'useParallel', 'yes', ...
0335                         'showResources', 'yes', ...
0336                         'reduction', 60);
0337 
0338                 end
0339                 
0340                 %[self.net, self.trdat] = train(self.net, x, y, 'useParallel', 'no', 'showResources', 'yes');
0341                 
0342                 % storing full traindata, valdata, testdata with network.
0343                 % The trainInd, valInd, testInd are only useful if I'm 100%
0344                 % sure that I get the same data, and with making
0345                 % subselections or getting more collocations, this is not
0346                 % the case
0347                 fitter.userdata.alldata = [x_orig_fit; y];
0348 %                fitter.userdata.traindata = [x_orig_fit(:, fit_trdat.trainInd); y(:, fit_trdat.trainInd)];
0349 %                fitter.userdata.valdata = [x_orig_fit(:, fit_trdat.valInd); y(:, fit_trdat.valInd)];
0350 %                fitter.userdata.testdata = [x_orig_fit(:, fit_trdat.testInd); y(:, fit_trdat.testInd)];
0351                 fitter.userdata.x_cols = x_fit_cols;
0352                 fitter.userdata.y_cols = ...
0353                     structfun(@(x)x+size(x_orig_fit, 1), ...
0354                               y_cols, 'UniformOutput', false);
0355                 
0356                 if classify
0357                     patter.userdata.traindata = ...
0358                         [x_orig_pat(:, pat_trdat.trainInd); ...
0359                          y_class(:, pat_trdat.trainInd)];
0360                     patter.userdata.valdata = ...
0361                         [x_orig_pat(:, pat_trdat.valInd); ...
0362                          y_class(:, pat_trdat.valInd)];
0363                     patter.userdata.testdata = ...
0364                         [x_orig_pat(:, pat_trdat.testInd); ...
0365                          y_class(:, pat_trdat.testInd)];
0366                     patter.userdata.cols = x_pat_cols;
0367 
0368                     patter.userdata.randseed = seed;
0369                     patter.userdata.noiselevel = noise_level;
0370                 end
0371 
0372                 fitter.userdata.randseed = seed;
0373                 fitter.userdata.noiselevel = noise_level;
0374 
0375                 logtext(atmlab('OUT'), 'Storing network\n');
0376                 logtext(atmlab('OUT'), 'Output: %s\n', f);
0377                 fitter.userdata.stored = f;
0378                 fitter.userdata.no = net_count;
0379                 fitter.userdata.note = self.name;
0380 
0381                 self.net = struct();
0382                 self.net.fit = fitter;
0383                 self.trdat.fit = fit_trdat;
0384 
0385                 if classify
0386                     patter.userdata.stored = f;
0387                     patter.userdata.no = net_count;
0388                     patter.userdata.note = fitter.userdata.note;
0389                     self.net.pat = patter;
0390                     self.trdat.pat = pat_trdat;
0391                 end
0392                 
0393                 self.storenets(f);
0394             end
0395             
0396         end
0397         
0398         function [y_class, y_dat] = get_classification(self, dat)
0399             % from data subselection, get classification targets
0400             %
0401             % must have exactly one classification target
0402             %
0403             % two outputs:
0404             % y_class   [2 x N] matrix, cloud/cloudfree
0405             % y_data    [1 x N] vector, corresponding measurement
0406             self.localcols;
0407             tgs = fieldnames(self.targets);
0408             y_class = zeros(2, size(dat, 1));
0409             done = false;
0410             for i = 1:length(tgs)
0411                 tg = tgs{i};
0412                 if isfield(self.targets.(tg), 'classify') && self.targets.(tg).classify
0413                     if done
0414                         error(['atmlab:' mfilename ':twofields'], ...
0415                             'Two fields with classify?!');
0416                     end
0417                     y_class(1, :) = dat(:, self.localcols.(tg))<=self.targets.(tg).classify;
0418                     y_class(2, :) = dat(:, self.localcols.(tg))>self.targets.(tg).classify;
0419                     y_dat = dat(:, self.localcols.(tg));
0420                     % FIXME: why do I get any where both are false?
0421                     done = true;
0422                 end
0423             end
0424             if ~done
0425                 error(['atmlab:' mfilename ':noclassify'], ...
0426                     'No fields with classify?!');
0427             end
0428         end
0429 
0430         
0431         function storenets(self, path)
0432             % store artificial neural network and training data to 'path'
0433            nt = self.net; %#ok<NASGU>
0434            tr = self.trdat; %#ok<NASGU>
0435            save(path, 'nt', 'tr');
0436         end        
0437         
0438         function loadnets(self, varargin)
0439             % Loads previously-trained neural-net
0440             %
0441             % FORMAT
0442             %
0443             %   nn.loadnets(p)
0444             %
0445             % IN
0446             %
0447             %   p    either full path or string corresponding to one of the
0448             %        members of self.networks, which translates to a full
0449             %        path.  Can also be omitted completely, will then try
0450             %        to load the latest version based on self.version.
0451             if isempty(self.version)
0452                 narginchk(2, 2);
0453                 p = varargin{1};
0454             else
0455                 narginchk(1, 2);
0456                 p = optargs(varargin, {['version_', strrep(self.version, '.', '_')]});
0457             end
0458             if isfield(self.networks, p)
0459                 path = self.networks.(p);
0460             else
0461                 path = p;
0462             end
0463             logtext(atmlab('OUT'), 'Loading network and training info from %s\n', ...
0464                 path);
0465             self.net = loadvar(path, 'nt');
0466             self.trdat = loadvar(path, 'tr');
0467         end
0468         
0469         function [best, perfs] = bestnet(self, noise_level)
0470             % when multiple nets were trained, select the best one
0471             
0472             logtext(atmlab('OUT'), 'Selecting best-performing network\n');
0473             % consider all targets
0474             targs = fieldnames(self.targets);
0475 
0476             % infinite loop to count how many networks we have
0477             total = 0;
0478             while true
0479                 f = self.calculate_network_file_name(noise_level, total);
0480                 if ~exist(f, 'file')
0481                     break
0482                 end
0483                 total = total + 1;
0484             end
0485 
0486             perfs = cell(total, length(targs));
0487             %perf = zeros(total, length(targs));
0488             perf_fit = zeros(total, length(targs));
0489             classify = zeros(length(targs), 1);
0490             for i = 0:total-1
0491                 f = self.calculate_network_file_name(noise_level, i);
0492                 self.loadnets(f);
0493                 % What input to test the performance for?  How to
0494                 % weigh 'best' when there are different inputs?
0495                 for t = 1:length(targs)
0496                     targ = targs{t};
0497                     [xx, p] = self.get_performance(noise_level, targ, ...
0498                         self.targets.(targ).bins);
0499                     perfs{i+1, t} = {xx, p};
0500                     perf(i+1, t) = p;
0501                     classify(t) = self.targets.(targ).classify;
0502                     if classify(t)
0503                         [~, fp, fn] = self.test_patternnet_performance();
0504                         classify(t) = true;
0505                         perf_fit(i+1, t) = mean(fp+fn);
0506                     end
0507                 end
0508             end
0509             quals = cell(length(targs), 1);
0510             for t = 1:length(targs)
0511                 targ = targs{t};
0512                 relevant = (all([perf(:, t).count]>100, 2));
0513                 if isfield(self.targets.(targ), 'relevant')
0514                     relevant = relevant & (perfs{1, t}{1}>self.targets.(targ).relevant);
0515                 end
0516                 allperf = [perf(:, t).(safegetfield(self.targets.(targ), 'optimise_field', 'rel_medfracerr'))];
0517                 quals{t} = mean(allperf(relevant, :));
0518                 if classify(t)
0519                     quals{t} = quals{t} + vec2row((perf_fit(:, t)./median(perf_fit(:, t)) > 1.2) * 100); % this means the cloud filter went horribly wrong
0520                 end
0521             end
0522             [~, best] = min(mean(vertcat(quals{:}), 1));
0523             best = best - 1; % count from 0
0524             % load the best net
0525             self.loadnets(self.calculate_network_file_name(noise_level, best));
0526         end
0527         
0528         
0529         function [x, y_ref, y_retr, x_cols, y_cols, y_retr_nonoise] ...
0530                 = evaluate_test_data(self, varargin)
0531             % evaluates test-data for NN
0532             %
0533             % FORMAT
0534             %
0535             %   [x, y_ref, y_retr, x_cols, y_cols, y_retr_nonoise] = prod.evaluate_test_data(noise_level)
0536             %
0537             % IN
0538             %
0539             %   noise_level     double      noise, optional, defaults to
0540             %                               training noise
0541             %
0542             % OUT
0543             %
0544             %   x
0545             %   y_ref
0546             %   y_retr
0547             %   x_cols
0548             %   y_cols
0549             
0550             ud = self.net.fit.userdata;
0551             indy_dat = self.data(self.trdat.fit.testInd, :);
0552             [y_ref, y_cols] = get_columns(indy_dat, ...
0553                                           self.localcols, ...
0554                                           fieldnames(ud.targets));
0555             [x, x_cols] = self.select_and_prepare_data(indy_dat, self.localcols, self.inputs);
0556             x = x.';
0557             %[x, x_cols] = get_columns(indy_dat, ...
0558             %                          self.localcols, ...
0559             %                          structfun(@(X) safegetfield(X, 'chans', -1), ud.inputs, 'UniformOutput', false));
0560 
0561             if isfield(ud, 'alldata')
0562                 x_cols = ud.x_cols;
0563                 seed = ud.randseed;
0564                 [xx, xx_cols] = self.select_and_prepare_data(...
0565                     ud.alldata.', ud.x_cols, ud.inputs);
0566                 xx = xx(:, self.trdat.fit.testInd).';
0567                 [yy_ref, yy_cols] = self.select_and_prepare_data(...
0568                     ud.alldata.', ud.y_cols, ud.targets);
0569                 yy_ref = yy_ref(:, self.trdat.fit.testInd).';
0570                 assert(isequal(xx, x));
0571 %                targs = fieldnames(y_cols);
0572 %                for t=1:length(targs)
0573 %                    targ = targs{t};
0574 %                    y_ref(:, t) = self.targets.(targ).invtransform(y_ref(:, t));
0575 %                end
0576             end
0577             rng('shuffle');
0578             noise_level = optargs(varargin, {ud.noiselevel});
0579             % FIXME: vertical correlation in error...
0580             if isempty(x)
0581                 noise = x;
0582             elseif isscalar(noise_level)
0583                 noise = randn(size(x)) * noise_level;
0584             elseif isempty(noise_level)
0585                 noise = bsxfun(@times, randn(size(x)), self.get_noise_vector(x_cols));
0586             else
0587                 noise = bsxfun(@times, randn(size(x)), noise_level);
0588             end
0589             y_retr_raw = sim(self.net.fit, (x + noise).').';
0590             y_retr_raw_nonoise = sim(self.net.fit, x.').';
0591             % invert for each retrieved field
0592             y_retr = y_retr_raw;
0593             y_retr_nonoise = y_retr_raw_nonoise;
0594             fns = fieldnames(ud.targets);
0595             for i = 1:length(fns)
0596                 fn = fns{i};
0597                 y_retr(:, y_cols.(fn)) = ud.targets.(fn).invtransform(y_retr_raw(:, y_cols.(fn)));
0598                 y_retr_nonoise(:, y_cols.(fn)) = ud.targets.(fn).invtransform(y_retr_raw_nonoise(:, y_cols.(fn)));
0599             end
0600             
0601             %y_retr = invtrans(y_retr_raw);
0602             %y_retr = y_retr(:);
0603         end
0604            
0605         function dev_selftest = evaluate_and_summarise_testdata(self, varargin)
0606            % evaluate test-data and perform several statistical tests
0607            %
0608            %
0609            % - make a scatter density plot
0610            %
0611            % IN
0612            %
0613            %    noise_level
0614            %
0615            %        optional, defaults to same as training.
0616            %        Can be either a scalar, or a
0617            %        vector corresponding to the number of channels, or []
0618            %        which means same as training.
0619            %
0620            % See also: scatter_density_plot, plot_scatter_perf
0621            
0622            noise_simul = optargs(varargin, {self.net.userdata.noiselevel});
0623            
0624            %% self-test
0625            
0626            figure('visible', 'off');
0627            [x, y] = self.evaluate_test_data(noise_simul);      
0628            filt = all(x>0, 2) & ...
0629                   all(y>0, 2) & ...
0630                   all(isfinite(x), 2) & ...
0631                   all(isfinite(y), 2);
0632 
0633            [h_ax, h_cb] = self.plot_scatter_perf(x(filt), y(filt));
0634                      
0635            ylabel(h_cb, 'Number of retrievals');
0636            xlabel(h_ax, 'Independent IWP [g/m^2]');
0637            ylabel(h_ax, 'Retrieved IWP [g/m^2]');
0638  
0639            if isfield(self.net.userdata, 'realnoise') && self.net.userdata.realnoise
0640                train_noise_str = 'realnoise';
0641            elseif isscalar(self.net.userdata.noiselevel)
0642                train_noise_str = sprintf('%.2fK', self.net.userdata.noiselevel); 
0643            else
0644                train_noise_str = [sprintf('%.2f,', self.net.userdata.noiselevel) 'K'];
0645            end
0646            
0647            if isequal(noise_simul, self.net.userdata.noiselevel)
0648                simul_noise_str = 'trainnoise';
0649            elseif isscalar(noise_simul)
0650                simul_noise_str = sprintf('%.2fK', noise_simul);
0651            else
0652                simul_noise_str = [sprintf('%.2f,', noise_simul) 'K'];
0653            end
0654            
0655            title(sprintf('NN-IWP self-test, Evans data, noise %s/%s\nchannels %s', train_noise_str, simul_noise_str, self.freqstr));
0656            print(gcf(), '-depsc', fullfile(cscol('plot_base'), sprintf('nniwp_evansdb_selftest_noise%s_%s_chans%s.eps', ...
0657                train_noise_str, simul_noise_str, self.chanstr)));
0658 
0659            %% calculate median fractional error per bin, store results in plotdata
0660            
0661            xbins = self.opts.bins;
0662            dev_selftest = self.binned_median_fractional_error(x(filt), y(filt), xbins);
0663            M = [xbins.', dev_selftest.']; %#ok<NASGU>
0664            of = fullfile(cscol('plotdata_base'), sprintf('nniwp_evansdb_selftest_fracerr_noise%s_%s_chans%s.dat', ...
0665                train_noise_str, simul_noise_str, self.chanstr));
0666            save(of, 'M', '-ascii');
0667         end
0668        
0669         function [xx, S] = get_performance(self, noise_level, varargin)
0670             % get many performance statistics
0671             %
0672             % Default, old-style, is hard-coded logarithmically for IWP.
0673             % For new behaviour, pass appropriate variables.
0674             %
0675             % IN
0676             %
0677             %   noise_level
0678             %   variable        (must be in self.targets)
0679             %   bins
0680 
0681             [x, y_ref, y_retr, x_cols, y_cols] = self.evaluate_test_data(noise_level);
0682             filt = (y_ref>0) & (y_retr>0) & isfinite(y_ref) & isfinite(y_retr);
0683             
0684             [variable, bins] = optargs(varargin, {'IWP', []});
0685 
0686             y_ref = y_ref(all(filt, 2), y_cols.(variable));
0687             y_retr = y_retr(all(filt, 2), y_cols.(variable));
0688             
0689             if isempty(bins)
0690                 bins = floor(min(log10(y_ref(:)))):.2:ceil(max(log10(y_ref(:))));
0691             end
0692 
0693             bins = vec2col(bins);
0694 
0695             func = self.targets.(variable).transform;
0696             defunc = self.targets.(variable).invtransform;
0697             % calculate different ways to get performance
0698             S.abs_medad = cellfun(@(yy) mad(yy, 1), bin(func(y_ref), y_retr, bins));
0699             S.abs_meanad = cellfun(@(yy) mad(yy, 0), bin(func(y_ref), y_retr, bins));
0700             S.abs_rmeanse = cellfun(@(dy) rms(dy), bin(func(y_ref), y_retr-y_ref, bins));
0701             S.abs_rmedse = cellfun(@(dy) sqrt(median(dy.^2)), bin(func(y_ref), y_retr-y_ref, bins));
0702             S.abs_ir = cellfun(@iqr, bin(func(y_ref), y_retr, bins));
0703             %S.rel_medfracerr = cellfun(@median, bin(func(y_ref), abs(y_retr-y_ref)./y_ref, bins));
0704             S.rel_medfracerr = cellfun(@median, bin(func(y_ref), defunc(abs(func(y_retr) - func(y_ref)))-1, bins));
0705             S.abs_medfracerr = S.rel_medfracerr .* defunc(bins);
0706             S.rel_rmeanse_frac = cellfun(@(dy) rms(dy), bin(func(y_ref), abs(y_retr-y_ref)./y_ref, bins));
0707             S.abs_rmeanse_frac = S.rel_rmeanse_frac .* defunc(bins);
0708             S.count = cellfun(@(yy) size(yy, 1), bin(func(y_ref), y_retr, bins));
0709             xx = defunc(bins);
0710         end
0711 
0712         function [bn, fp, fn, fp_med_ref, fp_med_retr, fn_med_ref, fn_med_retr] = test_patternnet_performance(self)
0713             % Test how well the cloudy/cloudfree recognition works
0714             %
0715             %p = self.net.pat;
0716             %
0717             % Returns:
0718             %
0719             % bins, false positives, false negatives
0720             [y_class, y_class_data] = self.get_classification(self.data);
0721             ref_clear = y_class(1, self.net.pat.divideParam.testInd);
0722             ref_cloudy = y_class(2, self.net.pat.divideParam.testInd);
0723             ref_data = y_class_data(self.net.pat.divideParam.testInd);
0724             [x, x_cols] = self.select_and_prepare_data(self.data, self.localcols, self.inputs);
0725             x_pat = self.subselect_inputs(x, x_cols, self.inputs, 'pat');
0726             x_fit = self.subselect_inputs(x, x_cols, self.inputs, 'fit');
0727 
0728             yclrt = self.net.pat(x_pat(:, self.net.pat.divideParam.testInd));
0729             % use the pattern testInd for the fit network here, because I'm
0730             % interested in retrieved IWP in particular where the reference
0731             % is 0
0732             ydatrt = 10.^self.net.fit(x_fit(:, self.net.pat.divideParam.testInd));
0733             
0734             valid = (isfinite(ref_data)&ref_data>=0).';
0735             
0736             bn = linspace(0, 1, 70);
0737             fp = zeros(size(bn));
0738             fn = zeros(size(bn));
0739             fp_med_ref = zeros(size(bn));
0740             fp_med_retr = zeros(size(bn));
0741             fn_med_ref = zeros(size(bn));
0742             fn_med_retr = zeros(size(bn));
0743             for i = 1:length(bn)
0744                 frc = bn(i);
0745                 retr_cloudy = yclrt(2, :) > frc;
0746                 retr_clear = ~retr_cloudy;
0747                 fp(i) = sum(retr_cloudy & ref_clear & valid)./length(find(valid));
0748                 fn(i) = sum(retr_clear & ref_cloudy & valid)./length(find(valid));
0749                 fp_med_ref(i) = median(ref_data(valid & retr_cloudy & ref_clear));
0750                 fp_med_retr(i) = median(ydatrt(valid & retr_cloudy & ref_clear));
0751                 fn_med_ref(i) = median(ref_data(valid & ref_cloudy & retr_clear));
0752                 fn_med_retr(i) = median(ydatrt(valid & ref_cloudy & retr_clear));
0753             end
0754         end
0755         
0756         function [ptls, ptls_nonoise] = iwc_profile_performance(self, ptiles)
0757             % get profile performance, return percentiles for ratio iwc/ref
0758             %
0759             % profile_performace([1, 25, 50, 75, 99])
0760             %
0761             % See also: prctile
0762             [x, y_ref, y_retr, x_cols, y_cols, y_retr_nonoise] = self.evaluate_test_data();
0763             fy = y_ref(:, y_cols.IWC) ./ y_retr(:, y_cols.IWC);
0764             fy_nonoise = y_ref(:, y_cols.IWC) ./ y_retr_nonoise(:, y_cols.IWC);
0765             ptls = prctile(fy, ptiles);
0766             ptls_nonoise = prctile(fy_nonoise, ptiles);
0767         end
0768         
0769         function dy = akm(self)
0770             % get averaging kernel matrix for all testdata
0771             %
0772             % OUT
0773             %
0774             %   dy  matrix of size (N, p, q) with averaging kernel matrix
0775             %       per measurement, perturbed per noise-level.  Resulting
0776             %       fractional change in IWC.
0777             
0778             % construct a N x p x p matrix where each page is
0779             n_inputs = self.net.inputs{1}.size;
0780             n_outputs = self.net.outputs{end}.size;
0781             n_ret = length(self.trdat.testInd);
0782             [x_unpert, x_cols] = get_columns(self.data(self.trdat.testInd, :), ...
0783                                              self.localcols, ...
0784                                              fieldnames(self.inputs));
0785             dx = self.get_noise_vector(x_cols);
0786             bsxfun(@times, eye(n_inputs), dx);
0787             x_pert = bsxfun(@plus, ...
0788                             repmat(x_unpert, [1, 1, n_inputs]), ...
0789                             reshape(bsxfun(@times, eye(n_inputs), dx), ...
0790                                     [1, n_inputs, n_inputs]));
0791                             
0792 %             x_pert = bsxfun(@times, ...
0793 %                             repmat(x_unpert, [1, 1, n_inputs]), ...
0794 %                             reshape(1 + eye(n_inputs)*pert_frac, [1, n_inputs, n_inputs]));
0795             y_unpert = sim(self.net, x_unpert.').';
0796             y_pert = permute(reshape(self.net(reshape(permute(x_pert, [2, 3, 1]), ...
0797                                                      [n_inputs, n_inputs*n_ret, 1])), ...
0798                                      [n_outputs, n_inputs, n_ret]), ...
0799                              [3, 1, 2]);
0800 
0801             dy = 10.^bsxfun(@minus, y_pert, y_unpert)-1;
0802         end
0803         
0804         function plot_akm(self, dy)
0805             % plot akm
0806             
0807             % FIXME: move to other class
0808             n_inputs = size(dy, 2);
0809             n_outputs = size(dy, 1);
0810             %pcolor(1:n_inputs, linspace(4.5, 15, n_outputs), squeeze(dy(1, :, :)));
0811             sanepcolor(1:n_inputs, linspace(4.5, 15, n_outputs), dy);
0812             shading faceted;
0813             colormap(drywet(0.1));
0814             cb = colorbar();
0815             ylabel(cb, '%IWC change per noise level measurement');
0816             zerobright();
0817             hold on;
0818             %plot([12, 12], [4.5, 15], 'LineWidth', 2.5, 'Color', 'black');
0819             xlabel(sprintf('Channel (%s)', strjoin(vec2row(fieldnames(self.inputs)), ', ')));
0820             ylabel('Height [km]');
0821             title(['Sensitivity analysis ' strrep(self.name, '_', ' ')]);            
0822             save_figure_multi(gcf(), fullfile(cscol('plot_base'), [self.name '_mean_akm']), 'png', 'eps', 'fig');
0823         end
0824         
0825         function nm = calculate_network_file_name(self, noise_level, i)
0826             inps = fieldnames(self.inputs);
0827             tgts = fieldnames(self.targets);
0828 
0829             inpcell = cell(size(inps));
0830             for k = 1:length(inpcell)
0831                 inp = inps{k};
0832                 if isfield(self.inputs.(inp), 'chans')
0833                     inpcell{k} = sprintf('%s,%s', ...
0834                         inp, ...
0835                         strjoin(arrayfun(@num2str, ...
0836                                         self.inputs.(inp).chans, ...
0837                                         'UniformOutput', false), ...
0838                                 ','));
0839                 else
0840                     inpcell{k} = inp;
0841                 end
0842             end
0843 
0844             tgtcell = cell(size(tgts));
0845             for k = 1:length(tgtcell)
0846                 tgt = tgts{k};
0847                 if isfield(self.targets.(tgt), 'regression_range')
0848                     tgtcell{k} = sprintf('%s,%.4g-%.4g', ...
0849                         tgt, ...
0850                         self.targets.(tgt).regression_range(1), ...
0851                         self.targets.(tgt).regression_range(2));
0852                 else
0853                     tgtcell{k} = tgt;
0854                 end
0855             end
0856             
0857             if length(unique(noise_level))==1
0858                 noise_str = sprintf('%.2f', unique(noise_level));
0859             else
0860                 noise_str = sprintf('%2f,', noise_level);
0861             end
0862             nm = sprintf('/storage3/user_data/gerrit/neuralnets/%s_%s_to_%s_noise%s_v%s_%d.mat', ...
0863                 self.name, ...
0864                 sprintf('%s,', strjoin(vec2row(inpcell), ',')), ...
0865                 sprintf('%s,', strjoin(vec2row(tgtcell), ',')), ...
0866                 noise_str, ...
0867                 self.version, ...
0868                 i);
0869             
0870         end
0871         
0872         function res = retrieve(self, data, cols, varargin)
0873             % retrieve data
0874             %
0875             % iwp = ccniwp.retrieve(data, cols[, cutoff])
0876             %
0877             % This method:
0878             %
0879             %  - applies limitations (lat/lon)
0880             %  - applies transforms and invtransforms
0881             %  - applies both neural networks
0882             %  - decides cloud/cloudfree with 'cutoff'
0883             %  - returns retrieved IWP in g/m^2
0884             %
0885             % IN
0886             %
0887             %   data    matrix with data, must have at least lat/lon and
0888             %           all the inputs
0889             %   cols    description of columns. NB: must match 'inputs'
0890             %   cutoff  cloud/nocloud?  Defaults to 0.5
0891             %
0892             % OUT
0893             %
0894             %   res - structure with many fields, depends on the exact
0895             %   product
0896            
0897             cutoff = optargs(varargin, {0.5});
0898             
0899             if isequal(self.net, [])
0900                 logtext(atmlab('OUT'), 'No networks loaded yet, will load default one for this version\n');
0901                 self.loadnets();
0902             end
0903             
0904             ud = self.net.fit.userdata;
0905             
0906             [x, x_cols] = self.select_and_prepare_data(data, cols, ud.inputs);
0907             
0908             flds = fieldnames(x_cols);
0909             for i = 1:length(flds)
0910                 fld = flds{i};
0911                 if isfield(ud.inputs.(fld), 'lims')
0912                     lims.(fld) = ud.inputs.(fld).lims;
0913                 end
0914             end
0915 
0916             % FIXME still hardcoded; somtimes wrong in neural net
0917             lims.B_SZA = [0, 85]; 
0918             lims.B_LAA = [-180, 180];
0919             lims.B_SAA = [-180, 180];
0920             %limmer = collocation_restrain(inp, limstruct2limmat(lims, inp_cols));
0921             limmer = collocation_restrain(x.', limstruct2limmat(getfields(lims, intersect(fieldnames(lims), fieldnames(x_cols))), x_cols));
0922             %inp = inp(limmer, :);
0923             x = x(:, limmer);
0924             data = data(limmer, :); % for lat/lon/time still ekep track
0925             lat = data(:, cols.LAT1);
0926             lon = data(:, cols.LON1);
0927 
0928             %cloudpat = self.net.pat(inp.').';
0929             cloudpat = self.net.pat(self.subselect_inputs(x, x_cols, ud.inputs, 'pat')).';
0930             cloudy = cloudpat(:, 2)>cutoff;
0931             
0932             %iwp = self.net.fit(inp.').';
0933             iwp = self.net.fit(self.subselect_inputs(x, x_cols, ud.inputs, 'fit')).';
0934 
0935             % FIXME: not so hardcoded to MEAN_IWP_2C please?
0936             res.raw_iwp = ud.targets.MEAN_IWP_2C.invtransform(iwp);
0937             res.iwp = res.raw_iwp;
0938             res.iwp(~cloudy) = 0;
0939             res.cloud_probability = cloudpat(:, 2);
0940             res.lat = lat;
0941             res.lon = lon;
0942             res.time = data(:, cols.TIME1);
0943             
0944             % also put inputs here
0945             
0946             inps = fieldnames(self.inputs);
0947             for i = 1:length(inps)
0948                 inp = inps{i};
0949                 if isfield(cols, inp)
0950                     res.(inp) = data(:, cols.(inp));
0951                 else
0952                     res.(inp) = x(x_cols.(inp), :);
0953                 end
0954             end
0955             
0956 %             res.MEAN_AVHRR_Y = data(:, cols.MEAN_AVHRR_Y);
0957 %             res.AVHRR_FLAG_3AB = data(:, cols.AVHRR_FLAG_3AB);
0958 %             res.MHS = data(:, cols.B_BT(3:5));
0959 %             res.B_LZA = data(:, cols.B_LZA);
0960 %             res.B_LAA = data(:, cols.B_LAA);
0961 %             res.B_SZA = data(:, cols.B_SZA);
0962 %             res.B_SAA = data(:, cols.B_SAA);
0963 %             res.Surface_elevation = x(x_cols.Surface_elevation, :);
0964 %             if isfield(x_cols, 'Surface_elevation_std')
0965 %                 res.Surface_elevation_std = x(x_cols.Surface_elevation_std, :);
0966 %             end
0967            
0968         end
0969         
0970         function [data, info] = retrieve_granule(self, granule, spec)
0971             % res = retrieve_granule(granule, spec)
0972             %
0973             % See also: retrieve, retrieve_and_store_granule
0974             
0975             altfile = self.get_alt_file(granule, spec);
0976             if isempty(altfile)
0977                 gotdata = false;
0978             else
0979                 S = self.read_alt_granule(granule, spec);
0980                 if isempty(S.TIME)
0981                     logtext(atmlab('OUT'), 'Read no data, skipping\n')
0982                     data = [];
0983                     info = struct();
0984                     return
0985                 else
0986                     gotdata = true;
0987  %                   [M, cc] = self.struct2mat(getfields(S, intersect(fieldnames(S), fieldnames(self.cols))));
0988                     [M, cc] = self.struct2mat(S);
0989                 end
0990             end
0991             
0992             if ~gotdata
0993                 D = datasets();
0994                 
0995                 [result, additional_results, ~] = D.collocation_mhs_avhrr.collocate_granule(granule, spec, spec, {D.associated_mhs_avhrr, D.collapsed_mhs_avhrr}, false);
0996 %                 lockfile = get_lock(atmlab('WORK_AREA'), 'collocations.lock');
0997 %                 cl1 = onCleanup(@()delete(lockfile));
0998                 if isempty(result)
0999                     logtext(atmlab('ERR'), 'Found no data, skipping\n');
1000                     data = [];
1001                     info = struct();
1002                     return
1003                 end
1004                 [M, cc] = D.associated_mhs_avhrr.merge_matrix(result, D.collocation_mhs_avhrr.cols, additional_results{1}, D.associated_mhs_avhrr.cols);
1005                 clear result;
1006                 [M, cc] = D.collapsed_mhs_avhrr.merge_matrix(M, cc, additional_results{2}, D.collapsed_mhs_avhrr.cols);
1007                 clear additional_results;
1008 %                 logtext(atmlab('OUT'), 'Clearing lockfile\n');
1009 %                 delete(lockfile);
1010             end
1011             
1012             logtext(atmlab('OUT'), 'Retrieving %d values\n', size(M, 1));
1013             res = self.retrieve(M, cc, 0.5);
1014             logtext(atmlab('OUT'), 'Done\n');
1015             
1016             
1017             % FIXME: make this more flexible, for going back to older
1018             % versions etc.
1019             data = zeros(size(res.iwp, 1), max(structfun(@(X)max(X), self.cols)));
1020             data(:, self.cols.LAT) = res.lat;
1021             data(:, self.cols.LON) = res.lon;
1022             data(:, self.cols.TIME) = res.time;
1023             data(:, self.cols.P_CLOUD) = res.cloud_probability;
1024             data(:, self.cols.IWP_RAW) = res.raw_iwp;
1025             data(:, self.cols.IWP) = res.iwp;
1026             %data(:, self.cols.MEAN_AVHRR_Y) = res.MEAN_AVHRR_Y;
1027             
1028             % and the inputs
1029             inps = fieldnames(self.inputs);
1030             for i = 1:length(inps)
1031                 inp = inps{i};
1032                 if isfield(self.cols, inp)
1033                     if (min(size(res.(inp))) ~= size(self.cols.(inp), 2)) && max(size(res.(inp))) > 3
1034                         data(:, self.cols.(inp)) = res.(inp)(:, self.inputs.(inp).chans);
1035                     else
1036                         data(:, self.cols.(inp)) = res.(inp);
1037                     end
1038                 end
1039             end
1040             
1041 %             data(:, self.cols.AVHRR_FLAG_3AB) = res.AVHRR_FLAG_3AB;
1042 %             data(:, self.cols.MHS) = res.MHS;
1043 %             data(:, self.cols.B_LZA) = res.B_LZA;
1044 %             data(:, self.cols.B_LAA) = res.B_LAA;
1045 %             data(:, self.cols.B_SZA) = res.B_SZA;
1046 %             data(:, self.cols.B_SAA) = res.B_SAA;
1047 %             data(:, self.cols.Surface_elevation) = res.Surface_elevation;
1048 %             if isfield(res, 'Surface_elevation_std') && isfield(self.cols, 'Surface_elevation_std')
1049 %                 data(:, self.cols.Surface_elevation_std) = res.Surface_elevation_std;
1050 %             end
1051             
1052             % FIXME: less hard-coded please
1053             info.Conventions = 'CF-1.6';
1054             info.title = 'AVHRR+MHS IWP trained by collocations with 2C-ICE';
1055             info.version = self.version;
1056             %info.history = self.changelog;
1057             info.source = 'Retrieved from collocated CPR/MHS 2C-ICE-trained ANN';
1058             info.nnet = self.net.pat.userdata.stored;
1059             info.references = 'Holl, G., S. Eliasson, J. Mendrok, and S. A. Buehler (submitted 2013), SPARE-ICE: synergistic IWP from passive operational sensors, J. Geophys. Res.';
1060             
1061             %            info.stored_lims = struct2string_compact(storelims);
1062             % FIXME: more additional info
1063             
1064             
1065         end
1066         
1067         function retrieve_and_store_granule(self, granule, spec)
1068             % ccniwp.retrieve_and_store_granule(granule, spec)
1069             %
1070             % Apply neural network retrieval to 'granule' / 'spec' and
1071             % store result according to self.find_granule_by_datetime.
1072             % If a retrieval has been performed previously for this
1073             % granule, this method reads the input data from the earlier
1074             % version.  Otherwise, it will read all inputs needed, perform
1075             % collocations (this will take some minutes).  Then, it applies
1076             % the neural network (loaded via self.loadnets) and stores the
1077             % result.
1078             %
1079             % IN
1080             %
1081             %   granule     datetime-vector describing granule start time
1082             %   spec        string describing satellite (e.g. 'noaa18')
1083             %
1084             % See also: retrieve, loadnets
1085             %
1086             % Example: D.col_syn_iwp.retrieve_and_store_granule([2007, 8, 1, 1, 48], 'noaa18');
1087             outfile = self.find_granule_by_datetime(granule, spec);
1088             
1089             if exist(outfile, 'file') && ~self.overwrite
1090                 logtext(atmlab('OUT'), 'Already exists: %s, continue\n', outfile);
1091                 return
1092             end
1093             
1094 %             if exist(outfile, 'file')
1095 %                 switch self.overwrite
1096 %                     case 0
1097 %                     case 1 % continue normally
1098 %                         gotdata = false;
1099 %                     case 2
1100 %                         S = self.read_granule(granule, spec, fieldnames(self.cols));
1101 %                         if isempty(S.time) % no data
1102 %                             logtext(atmlab('OUT'), 'Read no data, skipping\n')
1103 %                             return
1104 %                         end
1105 %                         [M, cc] = self.struct2mat(S);
1106 %                         gotdata = true;
1107 %                 end
1108 %             else
1109 %                 altfile = self.get_alt_file(granule, spec);
1110 %                 if isempty(altfile)
1111 %                     gotdata = false;
1112 %                 else
1113 %                     S = self.read_alt_granule(granule, spec);
1114 %                     if isempty(S.TIME)
1115 %                         logtext(atmlab('OUT'), 'Read no data, skipping\n')
1116 %
1117 %                         return
1118 %                     else
1119 %                         gotdata = true;
1120 %                         [M, cc] = self.struct2mat(S);
1121 %                     end
1122 %                 end
1123 %             end
1124 %             if ~gotdata
1125 %                 D = datasets();
1126 %
1127 %                 [result, additional_results, also] = D.collocation_mhs_avhrr.collocate_granule(granule, spec, spec, {D.associated_mhs_avhrr, D.collapsed_mhs_avhrr}, false);
1128 %                 lockfile = get_lock(atmlab('WORK_AREA'), 'collocations.lock');
1129 %                 cl1 = onCleanup(@()delete(lockfile));
1130 %                 if isempty(result)
1131 %                     logtext(atmlab('ERR'), 'Found no data, skipping\n');
1132 %                     return
1133 %                 end
1134 %                 [M, cc] = D.associated_mhs_avhrr.merge_matrix(result, D.collocation_mhs_avhrr.cols, additional_results{1}, D.associated_mhs_avhrr.cols);
1135 %                 clear result;
1136 %                 [M, cc] = D.collapsed_mhs_avhrr.merge_matrix(M, cc, additional_results{2}, D.collapsed_mhs_avhrr.cols);
1137 %                 clear additional_results;
1138 %                 logtext(atmlab('OUT'), 'Clearing lockfile\n');
1139 %                 delete(lockfile);
1140 %             end
1141 %
1142 %             logtext(atmlab('OUT'), 'Retrieving %d values\n', size(M, 1));
1143 %             res = self.retrieve(M, cc, 0.5);
1144 %             logtext(atmlab('OUT'), 'Done\n');
1145 %
1146 %
1147 %             % FIXME: make this more flexible, for going back to older
1148 %             % versions etc.
1149 %             data = zeros(size(res.iwp, 1), 20);
1150 %             data(:, self.cols.LAT) = res.lat;
1151 %             data(:, self.cols.LON) = res.lon;
1152 %             data(:, self.cols.TIME) = res.time;
1153 %             data(:, self.cols.P_CLOUD) = res.cloud_probability;
1154 %             data(:, self.cols.IWP_RAW) = res.raw_iwp;
1155 %             data(:, self.cols.IWP) = res.iwp;
1156 %             data(:, self.cols.MEAN_AVHRR_Y) = res.MEAN_AVHRR_Y;
1157 %             data(:, self.cols.AVHRR_FLAG_3AB) = res.AVHRR_FLAG_3AB;
1158 %             data(:, self.cols.MHS) = res.MHS;
1159 %             data(:, self.cols.B_LZA) = res.B_LZA;
1160 %             data(:, self.cols.B_LAA) = res.B_LAA;
1161 %             data(:, self.cols.B_SZA) = res.B_SZA;
1162 %             data(:, self.cols.B_SAA) = res.B_SAA;
1163 %             data(:, self.cols.Surface_elevation) = res.Surface_elevation;
1164 %             data(:, self.cols.Surface_elevation_std) = res.Surface_elevation_std;
1165 %
1166 %             % FIXME: less hard-coded please
1167 %             info.title = 'AVHRR+MHS IWP trained by collocations with 2C-ICE';
1168 %             info.version = 1;
1169 %             %info.history = self.changelog;
1170 %             info.source = 'Retrieved from collocated CPR/MHS 2C-ICE-trained ANN';
1171 %             info.nnet = self.net.pat.userdata.stored;
1172 % %            info.stored_lims = struct2string_compact(storelims);
1173 %             % FIXME: more additional info
1174             
1175             [data, info] = self.retrieve_granule(granule, spec);
1176             if isempty(data) % no data
1177                 logtext(atmlab('OUT'), 'Got no data, therefore not storing\n');
1178                 return
1179             end
1180             logtext(atmlab('OUT'), 'Storing\n');
1181             % ugly hack; I don't want to extend, but rewrite this one,
1182             % therefore temporarily changing self.overwrite
1183             if self.overwrite == 2
1184                 self.overwrite = 1;
1185                 self.store(granule, spec, data, info);
1186                 self.overwrite = 2;
1187             else
1188                 self.store(granule, spec, data, info);
1189             end
1190         end
1191         
1192         function retrieve_and_store_period(self, date1, date2, spec)
1193             D = datasets();
1194             % weirdness, see http://stackoverflow.com/q/19406551/974555
1195             m = D.mhs;
1196             allgrans = cached_evaluation(@m.find_granules_for_period, date1, date2, spec, 'EXTRA', self.name, m.name);
1197             
1198             for i = 1:size(allgrans, 1)
1199                 logtext(atmlab('OUT'), 'Processing: %s\n', num2str(allgrans(i, :)));
1200                 
1201                 try
1202                     self.retrieve_and_store_granule(allgrans(i, :), spec);
1203                 catch ME
1204                     switch ME.identifier
1205                         case {'atmlab:find_granule_by_datetime', 'atmlab:collocate', 'atmlab:CollocatedDataset:noother', 'MATLAB:nomem'}
1206                             logtext(atmlab('ERR'), 'Cannot retrieve for %s %s: %s\n', ...
1207                                 spec, num2str(allgrans(i, :)), ME.message);
1208                             continue
1209                         otherwise
1210                             ME.rethrow();
1211                     end
1212                 end
1213             end
1214             
1215         end
1216         
1217         function setcols(self)
1218             fs = fieldnames(self.members);
1219             N = 1;
1220             for i = 1:length(fs)
1221                 f = fs{i};
1222                 if isfield(self.members.(f), 'dims')
1223                     L = self.members.(f).dims{2};
1224                 else
1225                     L = 1;
1226                 end
1227                 self.cols.(f) = N:(N+L-1);
1228                 N = N + L;
1229             end
1230         end
1231         
1232         function setmembers(self)
1233             % FIXME: this should be more flexible...
1234             
1235             self.members.LAT.type = 'float';
1236             self.members.LAT.atts.long_name = 'Latitude';
1237             self.members.LAT.atts.valid_range = [-90 90];
1238             self.members.LAT.atts.units = 'degrees_north';
1239             
1240             self.members.LON.type = 'float';
1241             self.members.LON.atts.long_name = 'Longitude';
1242             self.members.LON.atts.units = 'degrees_east';
1243             self.members.LON.atts.valid_range = [-180 180];
1244             
1245             self.members.TIME.type = 'int';
1246             self.members.TIME.atts.long_name = 'Measurement time';
1247             self.members.TIME.atts.units = 'seconds since 1970-01-01 00:00:00';
1248            
1249             self.members.P_CLOUD.type = 'float';
1250             self.members.P_CLOUD.atts.long_name = 'Cloud probability';
1251             self.members.P_CLOUD.atts.valid_range = [0, 1];
1252 
1253             self.members.IWP_RAW.type = 'float';
1254             self.members.IWP_RAW.atts.long_name = 'Raw retrieved Ice Water Path, not considering cloud probability';
1255             self.members.IWP_RAW.atts.valid_min = 0;
1256             self.members.IWP_RAW.atts.units = 'g/m^2';
1257             self.members.IWP_RAW.atts.standard_name = 'atmosphere_cloud_ice_content';
1258 
1259             self.members.IWP.type = 'float';
1260             self.members.IWP.atts.long_name = 'Ice Water Path';
1261             self.members.IWP.atts.valid_min = 0;
1262             self.members.IWP.atts.units = 'g/m^2';
1263             self.members.IWP.atts.standard_name = 'atmosphere_cloud_ice_content';
1264 
1265             inp_fields = setdiff(fieldnames(self.inputs), {'LAT1'});
1266             for i = 1:length(inp_fields)
1267                 f = inp_fields{i};
1268                 if isfield(self.inputs, f) && isfield(self.inputs.(f), 'stored')
1269                     self.members.(f) = self.inputs.(f).stored;
1270                     if isfield(self.inputs.(f), 'chans')
1271                         self.members.(f).dims = {[f, '_CHANS'], length(self.inputs.(f).chans)};
1272                     end
1273                 end
1274             end
1275             
1276 %             self.members.MEAN_AVHRR_Y.type = 'float';
1277 %             self.members.MEAN_AVHRR_Y.atts.long_name = 'AVHRR measurement averaged over MHS footprint';
1278 %             self.members.MEAN_AVHRR_Y.atts.valid_min = 0;
1279 %             self.members.MEAN_AVHRR_Y.atts.units = 'reflectance or BT [K]';
1280 %             self.members.MEAN_AVHRR_Y.atts.origin = 'Based on AVHRR l1a data from NOAA CLASS archive';
1281 %             self.members.MEAN_AVHRR_Y.dims = {'AVHRR_CHANS', 5};
1282 %
1283 %             self.members.AVHRR_FLAG_3AB.type = 'int';
1284 %             self.members.AVHRR_FLAG_3AB.atts.long_name = 'AVHRR-3 flag: 3A/3B';
1285 %             self.members.AVHRR_FLAG_3AB.atts.valid_range = [0, 1];
1286 %             self.members.AVHRR_FLAG_3AB.atts.origin = 'Based on AVHRR L1 data from NOAA CLASS archive';
1287 %
1288 %             self.members.MHS.type = 'float';
1289 %             self.members.MHS.atts.long_name = 'MHS radiance, channels 3-5';
1290 %             self.members.MHS.atts.valid_range = [100, 400];
1291 %             self.members.MHS.atts.units = 'K';
1292 %             self.members.MHS.atts.origin = 'Copied from MHS L1 data from NOAA CLASS archive';
1293 %             self.members.MHS.dims = {'MHS_HUM_CHANS', 3};
1294 %
1295 %             self.members.B_LZA.type = 'float';
1296 %             self.members.B_LZA.atts.long_name = 'MHS local zenith angle';
1297 %             self.members.B_LZA.atts.valid_range = [0, 180];
1298 %             self.members.B_LZA.atts.units = 'degrees';
1299 %             self.members.B_LZA.atts.origin = self.members.MHS.atts.origin;
1300 %
1301 %             self.members.B_LAA.type = 'float';
1302 %             self.members.B_LAA.atts.units = 'degrees';
1303 %             self.members.B_LAA.atts.long_name = 'MHS local azimuth angle';
1304 %             self.members.B_LAA.atts.valid_range = [0, 360];
1305 %             self.members.B_LAA.atts.origin = self.members.MHS.atts.origin;
1306 %
1307 %             self.members.B_SZA.type = 'float';
1308 %             self.members.B_SZA.atts.units = 'degrees';
1309 %             self.members.B_SZA.atts.long_name = 'Solar zenith angle';
1310 %             self.members.B_SZA.atts.valid_range = [0, 180];
1311 %             self.members.B_SZA.atts.origin = self.members.MHS.atts.origin;
1312 %
1313 %             self.members.B_SAA.type = 'float';
1314 %             self.members.B_SAA.atts.units = 'degrees';
1315 %             self.members.B_SAA.atts.long_name = 'Solar azimuth angle';
1316 %             self.members.B_SAA.atts.valid_range = [0, 360];
1317 %             self.members.B_SAA.atts.origin = self.members.MHS.atts.origin;
1318 %
1319 %             self.members.Surface_elevation.type = 'float';
1320 %             self.members.Surface_elevation.atts.units = 'km';
1321 %             self.members.Surface_elevation.atts.long_name = 'Surface_elevation';
1322 %             self.members.Surface_elevation.atts.valid_range = [0, 10];
1323 %             self.members.Surface_elevation.atts.origin = 'NOAA ETOPO1 1 Arc-Minute Global Relief Model';
1324             
1325 %             self.members.Surface_elevation_std.type = 'float';
1326 %             self.members.Surface_elevation_std.atts.units = 'km';
1327 %             self.members.Surface_elevation_std.atts.long_name = 'Surface_elevation_std';
1328 %             self.members.Surface_elevation_std.atts.description = 'Standard deviation of 9 surface elevations in gridbox + neighbours';
1329 %             self.members.Surface_elevation_std.atts.valid_min = 0;
1330 %             self.members.Surface_elevation.atts.origin = 'Based on NOAA ETOPO1 1 Arc-Minute Global Relief Model';
1331            
1332         end
1333         
1334         function cc = fix_cols_in_postprocessing(~, cc)
1335         end
1336         
1337         function [x_new, x_new_cols] = subselect_inputs(~, x_all, x_all_cols, inps, mod)
1338             flds = fieldnames(inps);
1339             x_new = zeros(size(x_all));
1340             k = 0;
1341             for i = 1:length(flds)
1342                 fld = flds{i};
1343                 if isfield(inps.(fld), 'net')
1344                     if ~isequal(inps.(fld).net, mod)
1345                         continue
1346                     end
1347                 end
1348                 % chans are already taken care of in
1349                 % select_and_prepare_data!
1350 %                if isfield(inps.(fld), 'chans')
1351 %                    chans = inps.(fld).chans;
1352 %                else
1353 %                    chans = 1:length(x_all_cols.(fld));
1354 %                end
1355                 %newslc = (k+1):(k+length(x_all_cols.(fld)));
1356                 newslc = (k+1):(k+length(x_all_cols.(fld)));
1357                 %oldslc = par(x_all_cols.(fld), chans);
1358                 x_new(newslc, :) = x_all(x_all_cols.(fld), :);
1359                 x_new_cols.(fld) = newslc;
1360                 k = k + length(x_all_cols.(fld));
1361             end
1362             x_new(k+1:end, :) = [];
1363         end
1364         
1365 
1366 %         function [M, cc] = struct2mat(self, S)
1367 %             % copy it all over to the matrix
1368 %             flds = fieldnames(self.cols);
1369 %             M = zeros(size(S.time, 1), max(cellfun(@(X)max(self.cols.(X)), fieldnames(self.cols))));
1370 %             for i = 1:length(flds)
1371 %                 fld = flds{i};
1372 %                 M(:, self.cols.(fld)) = S.(fld);
1373 %             end
1374 %             cc = self.cols;
1375 %             % column names are different, channels etc...
1376 %             % expand this
1377 %             cc = self.fix_cols_in_postprocessing(cc);
1378 %         end
1379         
1380         % attempt (so far) to make this more generic
1381         %
1382         function [M, cc] = struct2mat(self, S)
1383             % copy it all over to the matrix
1384             flds = fieldnames(S);
1385             L = max(structfun(@(X) size(X, 1), S));
1386             %M = zeros(size(S.(flds{1}), 1), sum(structfun(@(X) size(X, 2), S)));
1387             N = 1;
1388             for i = 1:length(flds)
1389                 fld = flds{i};
1390                 if ~(size(S.(fld), 1)==L)
1391                     continue;
1392                 end
1393                 nn = size(S.(fld), 2);
1394                 if isfield(self.inputs, fld) && isfield(self.inputs.(fld), 'chans') && length(self.inputs.(fld).chans) < size(S.(fld), 2)
1395                     nn = length(self.inputs.(fld).chans);
1396                 end
1397                 nc = N:(N+nn-1);
1398                 N = N + nn;
1399                 cc.(fld) = nc;
1400                 if isfield(self.inputs, fld) && isfield(self.inputs.(fld), 'chans') && length(self.inputs.(fld).chans) < size(S.(fld), 2)
1401                     M(:, nc) = S.(fld)(:, self.inputs.(fld).chans);
1402                 else
1403                     M(:, nc) = S.(fld);
1404                 end
1405             end
1406             %cc = self.cols;
1407             % column names are different, channels etc...
1408             % expand this
1409             cc = self.fix_cols_in_postprocessing(cc);
1410         end
1411         
1412         function [fn, s, altfile, altbase] = get_alt_file(self, granule, spec)
1413             % return path to earlier file, if available, and version to
1414             % which it belongs
1415             
1416             
1417             cr_fn = self.filename;
1418             cr_bd = self.basedir;
1419             co1 = onCleanup(@()setfield(self, 'filename', cr_fn)); %#ok
1420             co2 = onCleanup(@()setfield(self, 'basedir', cr_bd)); %#ok
1421             
1422             if ~isempty(self.altfilename)
1423                 self.filename = self.altfilename;
1424             end
1425             
1426             if ~isempty(self.altbasedir)
1427                 self.basedir = self.altbasedir;
1428             end
1429             
1430             % try old versions
1431            
1432             curver = self.version;
1433             co = onCleanup(@()setfield(self, 'version', curver)); %#ok
1434             for s = {'0.7', '0.6', '0.5b', '0.5', '0.4', '0.3', '0.2', '0.1'}
1435                 s = s{1};
1436                 self.version = s;
1437                 fn = self.find_granule_by_datetime(granule, spec);
1438                 if exist(fn, 'file') % success
1439                     break
1440                 else % reset
1441                     fn = [];
1442                 end
1443             end
1444             self.version = curver;
1445             altfile = self.filename;
1446             altbase = self.basedir;
1447             
1448 %             if isempty(fn) % no luck yet
1449 %
1450 %                 cr_fn = self.filename;
1451 %                 cr_bd = self.basedir;
1452 %                 co1 = onCleanup(@()setfield(self, 'filename', cr_fn));
1453 %                 co2 = onCleanup(@()setfield(self, 'basedir', cr_bd));
1454 %                 self.filename = self.altfilename;
1455 %                 self.basedir = self.altbasedir;
1456 %                 %             cand =
1457 %                 fn = self.find_granule_by_datetime(granule, spec);
1458 %                 if ~exist(fn, 'file')
1459 %                     fn = [];
1460 %                 end
1461             self.filename = cr_fn;
1462             self.basedir = cr_bd;
1463             self.pos2re(); % fix regexp that may have been messed up
1464 %                 s = 0;
1465 %                 altfile = self.altfilename;
1466 %                 altbase = self.altbasedir;
1467 %             end
1468 %             self.filename = cr_fn;
1469 %             if ~b % try with old versions
1470 %                 curver = self.version;
1471 %                 co2 = onCleanup(@()setfield(self, 'version', curver));
1472 %                 for s = {'0.5b', '0.5', '0.4', '0.3', '0.2', '0.1'}
1473 %                     self.version = s;
1474 %                 end
1475 %             end
1476         end
1477         
1478         function S = read_alt_granule(self, granule, spec)
1479             
1480 %             cr_fn = self.filename;
1481 %             self.filename = self.altfilename;
1482 %             co = onCleanup(@()setfield(self, 'filename', cr_fn)); %#ok<SFLD>
1483 %             ud = self.net.fit.userdata;
1484 %             fields = fieldnames(ud.inputs);
1485 %             needsproc = cellfun(@(X)isfield(ud.inputs.(X), 'process'), fieldnames(ud.inputs));
1486             [~, vers, altfile, altbase] = self.get_alt_file(granule, spec);
1487             crvers = self.version;
1488             crfile = self.filename;
1489             crbase = self.basedir;
1490             self.version = vers;
1491             self.filename = altfile;
1492             self.basedir = altbase;
1493             co1 = onCleanup(@()setfield(self, 'version', crvers)); %#ok
1494             co2 = onCleanup(@()setfield(self, 'filename', crfile)); %#ok
1495             co3 = onCleanup(@()setfield(self, 'basedir', crbase)); %#ok
1496             self.pos2re();
1497             try
1498                 S = self.read_granule(granule, spec, fieldnames(self.members));
1499             catch ME
1500                 switch ME.identifier
1501                     case 'atmlab:SatDataset:cannotread'
1502                         % perhaps it has additional inputs compared to
1503                         % 'back then'.  Try to read only the ones I really
1504                         % need, obtain the rest through post-processing
1505                         inpfields = fieldnames(self.inputs);
1506 %                        flds = [intersect(fieldnames(self.members), ...
1507                         flds = cellfun(@(X) ...
1508                                    safegetfield(self.inputs.(X).stored, 'realname', X), ...
1509                                    inpfields(structfun(@(X) ~isfield(X, 'process'), self.inputs)), ...
1510                                    'UniformOutput', false);
1511                         %inpfields(structfun(@(X) ~isfield(X, 'process'), self.inputs))); 'LAT'];
1512                         S = self.read_granule(granule, spec, flds);
1513                         S = self.reader_processor(self, S, flds);
1514                     otherwise
1515                         ME.rethrow();
1516                 end
1517             end
1518             %S = self.read_granule(granule, spec);
1519             self.version = crvers;
1520             self.filename = crfile;
1521             self.basedir = crbase;
1522             self.pos2re(); % fix regexp that may have been messed up
1523             %S = self.read_granule(granule, spec, fieldnames(self.members));
1524             %S.Surface_elevation = get_surface_elevation(S.lat, S.lon);
1525 %             self.filename = cr_fn;
1526         end
1527     end
1528     
1529     
1530     methods (Access = {?SatDataset})
1531         
1532         function [S, strattr] = read_homemade_granule(self, file, varargin)
1533             fields = optargs(varargin, {{}});
1534             fields = union(vec2row(fields), {'LAT', 'LON', 'TIME'});
1535             [S, strattr] = read_homemade_granule@RetrievalDatabaseProduct(self, file, fields);
1536             S.lat = S.LAT;
1537             S.lon = S.LON;
1538             S.time = S.TIME;
1539             S = rmfield(S, setdiff({'LAT', 'LON', 'TIME'}, fields));
1540         end
1541         
1542     end
1543     
1544     methods (Access = protected)
1545         function noisevec = get_noise_vector(self, cols)
1546             % get a noise-vector responding to the columns of input data
1547             %
1548             % takes a structure 'cols', each field should correspond to a
1549             % self.inputs.(field) which must have a 'noise' vector
1550             inp_fields = fieldnames(self.inputs);
1551             %noisevec = zeros(self.net.inputs{1}.size, 1);
1552             for i = 1:length(inp_fields)
1553                 fn = inp_fields{i};
1554                 if isfield(self.inputs.(fn), 'chans')
1555                     noisevec(cols.(fn)) = self.inputs.(fn).noise(self.inputs.(fn).chans);
1556                 else
1557                     noisevec(cols.(fn)) = self.inputs.(fn).noise;
1558                 end
1559             end
1560         end
1561 
1562         function mask = select_for_regression(self, y, y_cols)
1563             % get a mask of those y-targets we want for regression
1564             %
1565             % If self.targets.(field).regression_range is set, select those
1566             % inside the range.  Otherwise, select all that are finite.
1567             %
1568             % WARNING! Input data from y is expected to have passed
1569             % through self.select_and_prepare_data, meaning that
1570             % 'transform' has been applied!
1571             %
1572             % IN
1573             %
1574             %   y
1575             %   y_cols
1576             mask = all(isfinite(y), 1);
1577             fields = fieldnames(y_cols);
1578             for i = 1:length(fields)
1579                 field = fields{i};
1580                 if isfield(self.targets.(field), 'regression_range')
1581                     rng = self.targets.(field).regression_range;
1582                     % unapply transformation
1583                     dat = self.targets.(field).invtransform(y(y_cols.(field), :));
1584                     mask = mask & (dat > rng(1));
1585                     mask = mask & (dat < rng(2));
1586                 end
1587             end
1588         end
1589     end
1590     
1591     methods (Static)
1592         function dat = homog_inputs(dat, cols)
1593             % homogeonise in latitude (equal area-density)
1594             
1595             % Calculate how many entries to choose per
1596             % 1°-latitude band.
1597             bins = -89.5:1:89.5;
1598             vals = bin(dat(:, cols.LAT1), 1:size(dat, 1), bins);
1599             %[C, X] = hist(dat(:, cols.LAT1), floor(min(dat(:, cols.LAT2)))-0.5:1:ceil(max(dat(:, cols.LAT2)))+0.5);
1600             no_entries_per_lat = round(floor(par(cellfun(@length, vals), 45)/1000)*1000*cosd(bins));
1601             inds = arrayfun(@(i) vals{i}(1:min([length(vals{i}), no_entries_per_lat(i)])), 1:180, 'UniformOutput', false);
1602             dat = dat(sort(vertcat(inds{:})), :);
1603         end
1604     end
1605 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005