function compare_collocated_iwp_products
%
% See also: CollocatedNNIWP, NNTrainedProduct, InstrumentVisualiser,
%           CollocationRetrievalDatabaseProduct, compare_ice_retrievals

% which input-combination to include
inps_solar = {'avhrr_12', 'avhrr_12_angles', 'avhrr_12_lat', 'avhrr_12_lat_angles'};
inps_tir = {'avhrr_45', 'avhrr_45_tsurf', 'avhrr_45_angles', 'avhrr_45_lat', 'avhrr_45_lat_tsurf', 'avhrr_45_lat_tsurf_angles'};
inps_mw = {'mhs_345_only', 'mhs_345_lat', 'mhs_345_angles', 'mhs_345_lat_angles'};
% inps = inps_mw;
% inps = {'avhrr_12_lat_angles', 'avhrr_45_lat_tsurf', 'mhs_345_lat_angles', ...
%      'avhrr_1245_lat_angles_tsurf', 'avhrr_1245_mhs_345_lat_angles_tsurf', ...
%      'avhrr_12345_mhs_345_lat_angles', 'avhrr_12345_mhs_345_lat_angles_elevn'};
% inps = {'avhrr_12_mhs_345_lat_angles', 'avhrr_45_mhs_345_lat_angles_tsurf'}
%inps = {'avhrr_12345_mhs_345_lat_angles', 'avhrr_12345_mhs_345_lat_angles_elevn', 'avhrr_12345_mhs_345_lat_angles_elevn_dz'};
%inps = { 'avhrr_12345_angles_elevn_tsurfcfsr', 'avhrr_345_angles_elevn_tsurfcfsr', 'avhrr_12345_mhs_345_lat_angles_elevn_dz', 'avhrr_12345_lat_angles_elevn_dz', 'avhrr_12345_lat_angles_elevn', 'avhrr_12345_lat_lza_sza_saa_elevn'};
%inps = {'avhrr_345_angles_elevn_tsurfcfsr', 'avhrr_12345_mhs_345_lat_angles_elevn_dz', 'avhrr_12345_lat_angles_elevn_dz', 'avhrr_12345_lat_angles_elevn', 'avhrr_12345_lat_lza_sza_saa_elevn'};
%inps = {'avhrr_12345_mhs_345_lat_angles_elevn_tsurfcfsr', ...
%        'avhrr_12345_mhs_345_angles_elevn_tsurfcfsr', ...
%        'avhrr_12345_mhs_345_lat_angles_elevn', ...
%        'avhrr_345_mhs_345_lat_angles_elevn_tsurfcfsr', ...
%        'avhrr_345_mhs_345_angles_elevn_tsurfcfsr', ...
%        'avhrr_12345_mhs_345_lat_angles_tsurfcfsr', ...
%        'avhrr_12345_mhs_345_angles_tsurfcfsr', ...
%        'avhrr_12345_mhs_345_lat_angles', ...
%        'avhrr_12345_mhs_345_angles', ...
%        'avhrr_345_mhs_345_lat_angles_tsurfcfsr', ...
%        'avhrr_345_mhs_345_angles_tsurfcfsr', ...
%        'avhrr_345_mhs_345_lat_angles'};
% inps = {'avhrr_12345_mhs_345_lat_angles_elevn_tsurfcfsr', ...
%         'avhrr_12345_mhs_345_angles_elevn_tsurfcfsr', ...
%         'avhrr_12345_mhs_345_lat_angles_elevn', ...
%         'avhrr_345_mhs_345_lat_angles_elevn_tsurfcfsr', ...
%         'avhrr_345_mhs_345_angles_elevn_tsurfcfsr', ...
%         'avhrr_345_mhs_345_lat_angles_elevn', ...
%         'avhrr_345_mhs_345_lat_angles_tsurfcfsr', ...
%         'avhrr_345_mhs_345_angles_tsurfcfsr', ...
%         'avhrr_345_mhs_345_lat_angles'};
inps = {'avhrr_345_mhs_345_angles_tsurfcfsr'}; 
% inps = {'avhrr_12345_mhs_345_lat_angles_elevn_tsurfcfsr', ...
%         'avhrr_12345_mhs_345_angles_elevn_tsurfcfsr', ...
%         'avhrr_12345_mhs_345_lat_angles_elevn', ...
%         'avhrr_12345_mhs_345_angles_elevn', ...
%         'avhrr_345_mhs_345_lat_angles_elevn_tsurfcfsr', ...
%         'avhrr_345_mhs_345_angles_elevn_tsurfcfsr', ...
%         'avhrr_345_mhs_345_lat_angles_elevn', ...
%         'avhrr_345_mhs_345_angles_elevn', ...
%         'avhrr_12345_mhs_345_lat_angles_tsurfcfsr', ...
%         'avhrr_12345_mhs_345_angles_tsurfcfsr', ...
%         'avhrr_12345_mhs_345_lat_angles', ...
%         'avhrr_12345_mhs_345_angles', ...
%         'avhrr_345_mhs_345_lat_angles_tsurfcfsr', ...
%         'avhrr_345_mhs_345_angles_tsurfcfsr', ...
%         'avhrr_345_mhs_345_lat_angles', ...
%         'avhrr_345_mhs_345_angles'};

%inps = {'avhrr_12345_mhs_345_lat_angles_elevn_dz'}; % = spareice
collocation_limits = struct('B_BT', [100 400], ...
    'MEAN_AVHRR_Y', [0 400], ...
    'LAT1', [-90, 90], ...
    'B_SZA', [0 85]);

limname = 'global';
process_month = true;
no_nets = 5;
homog = true;

ylims = [0, 300];
%ylims = [0, 800]; % avhrr 45 global peak
%ylims = [0, 1200]; % avhrr 12 global peak
%ylims = [0, 450]; % mhs 345 global peak (not really)
% 
%             self.members.AVHRR_FLAG_3AB.type = 'int';
%             self.members.AVHRR_FLAG_3AB.atts.long_name = 'AVHRR-3 flag: 3A/3B';
%             self.members.AVHRR_FLAG_3AB.atts.valid_range = [0, 1];
%             self.members.AVHRR_FLAG_3AB.atts.origin = 'Based on AVHRR L1 data from NOAA CLASS archive';
% 
             

inputs.avhrr_12 = struct(...
        'MEAN_AVHRR_Y', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [0, 100], ...
            'chans', [1, 2]));

inputs.avhrr_12.MEAN_AVHRR_Y.stored.type = 'float';
inputs.avhrr_12.MEAN_AVHRR_Y.stored.atts.valid_min = 0;
inputs.avhrr_12.MEAN_AVHRR_Y.stored.long_name = 'AVHRR measurements averaged over MHS footprint';
inputs.avhrr_12.MEAN_AVHRR_Y.stored.atts.units = 'reflectance or BT [K]';
inputs.avhrr_12.MEAN_AVHRR_Y.stored.atts.origin = 'Based on AVHRR l1a data from NOAA CLASS archive';

inputs.avhrr_123 = struct(...
        'MEAN_AVHRR_Y', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [0, 100], ...
            'chans', [1, 2, 3], ...
            'stored', inputs.avhrr_12.MEAN_AVHRR_Y.stored));
        
inputs.avhrr_1245 = struct(...
        'MEAN_AVHRR_Y', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [0, 400], ...
            'chans', [1, 2, 4, 5], ...
            'stored', inputs.avhrr_12.MEAN_AVHRR_Y.stored));

inputs.avhrr_12345 = struct(...
        'MEAN_AVHRR_Y', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [0, 400], ...
            'chans', [1, 2, 3, 4, 5], ...
            'stored', inputs.avhrr_12.MEAN_AVHRR_Y.stored));

angles = struct(...
        'B_LZA', struct(...
            'transform', @cosd, ...
            'invtransform', @acosd, ...
            'lims', [0, 180]), ...
        'B_LAA', struct(...
            'transform', @cosd, ...
            'invtransform', @acosd, ...
            'lims', [-180, 180]), ...
        'B_SZA', struct(...
            'transform', @cosd, ...
            'invtransform', @acosd, ...
            'lims', [0, 85]), ...
        'B_SAA', struct(...
            'transform', @cosd, ...
            'invtransform', @acosd, ...
            'lims', [-180, 180]));
        
angles.B_LZA.stored.type = 'float';
angles.B_LZA.stored.atts.long_name = 'MHS local zenith angle';
angles.B_LZA.stored.atts.valid_range = [0, 180];
angles.B_LZA.stored.atts.units = 'degrees';
angles.B_LZA.stored.atts.origin = 'Copied from MHS L1 data from NOAA CLASS archive';
% 
angles.B_LAA.stored.type = 'float';
angles.B_LAA.stored.atts.units = 'degrees';
angles.B_LAA.stored.atts.long_name = 'MHS local azimuth angle';
angles.B_LAA.stored.atts.valid_range = [0, 360];
angles.B_LAA.stored.atts.origin = angles.B_LZA.stored.atts.origin;
% 
angles.B_SZA.stored.type = 'float';
angles.B_SZA.stored.atts.units = 'degrees';
angles.B_SZA.stored.atts.long_name = 'Solar zenith angle';
angles.B_SZA.stored.atts.valid_range = [0, 180];
angles.B_SZA.stored.atts.origin = angles.B_LZA.stored.atts.origin;
% 
angles.B_SAA.stored.type = 'float';
angles.B_SAA.stored.atts.units = 'degrees';
angles.B_SAA.stored.atts.long_name = 'Solar azimuth angle';
angles.B_SAA.stored.atts.valid_range = [0, 360];
angles.B_SAA.stored.atts.origin = angles.B_LZA.stored.atts.origin;
%
lat = struct(...
        'LAT1', struct(...
            'transform', @cosd, ...
            'invtransform', @acosd, ...
            'lims', [-90, 90]));
lat.LAT1.stored.type = 'float';
lat.LAT1.stored.atts.units = 'degrees_north';
lat.LAT1.stored.atts.long_name = 'Latitude (used as input)';
lat.LAT1.stored.atts.valid_range = [-90, 90];
lat.LAT1.stored.realname = 'LAT';


tsurf = struct(...
        'MEAN_ECMWF_Skin_temperature', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [200, 400]));
        
tsurfcfsr = struct(...
        'CFSR_Skin_temperature', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [200, 400], ...
            'process', ...
                {{@get_cfsr_skin_temperature, ...
                    struct(...
                        'LAT1', ...
                            struct(...,
                                'transform', @(x)x, ...
                                'invtransform', @(x)x, ...
                                'lims', [-90, 90]), ...
                        'LON1', ...
                            struct(...
                                'transform', @(x)x, ...
                                'invtransform', @(x)x, ...
                                'lims', [-180, 180]), ...
                        'TIME1', ...
                            struct(...
                                'transform', @(x)x, ...
                                'invtransform', @(x)x, ...
                                'lims', [-180, 180]))}}, ...
            'stored', struct(...
                'type', 'float', ...
                'atts', struct(...
                    'units', 'K', ...
                    'long_name', 'NCEP CFSR skin temperature', ...
                    'valid_range', [0, 400], ...
                    'origin', 'NCEP CFSR ds093.1'))));

elev = struct(...
        'Surface_elevation', struct(...
                'process', ...
                {{@get_surface_elevation, ...
                 struct(...
                    'LAT1', ...
                        struct(...
                            'transform', @(x)x, ...
                            'invtransform', @(x)x, ...
                            'lims', [-90, 90]), ...
                    'LON1', ...
                        struct(...
                            'transform', @(x)x, ...
                            'invtransform', @(x)x, ...
                            'lims', [-180, 180]))}}, ...
        'transform', @(x)max(x, 0)/1000, ... % NB: see below
        'invtransform', @(x)x*1000));
    
elev.Surface_elevation.stored.type = 'float';
elev.Surface_elevation.stored.atts.units = 'km';
elev.Surface_elevation.stored.atts.long_name = 'Surface_elevation';
elev.Surface_elevation.stored.atts.valid_range = [0, 10];
elev.Surface_elevation.stored.atts.origin = 'NOAA ETOPO1 1 Arc-Minute Global Relief Model';

% NB: rather than max(x, 0) as above, more correct would be to set to 0
% when ocean, leave as it is otherwise.  However, we can likely safely
% consider sub-ocean land surfaces as 0-elevation land surfaces for the
% present purpose.
    
dz = struct(...
        'Surface_elevation_std', struct(...
                'process', ...
                {{@(lat, lon) get_surface_elevation(lat, lon, 1), ...
                struct(...
                    'LAT1', elev.Surface_elevation.process{2}.LAT1, ...
                    'LON1', elev.Surface_elevation.process{2}.LON1)}}, ...
        'transform', elev.Surface_elevation.transform, ...
        'invtransform', elev.Surface_elevation.invtransform));

% FIXME: calculate one or two relevant angles from this?
%
% according to http://en.wikipedia.org/wiki/Great-circle_distance#Formulas

inputs.avhrr_12_angles = catstruct(...
    inputs.avhrr_12, angles);

inputs.avhrr_12_lat = catstruct(...
    inputs.avhrr_12, lat);

inputs.avhrr_12_lat_angles = catstruct(...
    inputs.avhrr_12_lat, angles);

inputs.avhrr_123_angles = catstruct(...
    inputs.avhrr_123, angles);

inputs.avhrr_1245_angles = catstruct(...
    inputs.avhrr_1245, angles);

inputs.avhrr_1245_lat_angles = catstruct(...
    inputs.avhrr_1245, lat, angles);

inputs.avhrr_1245_lat_angles_tsurf = catstruct(...
    inputs.avhrr_1245, lat, angles, tsurfcfsr);

inputs.avhrr_12345_angles = catstruct(...
    inputs.avhrr_12345, angles);

inputs.avhrr_12_lat_angles = catstruct(...
    inputs.avhrr_12_angles, lat);

inputs.avhrr_12345_lat_angles = catstruct(...
    inputs.avhrr_12345_angles, lat);

inputs.avhrr_12345_lat_angles_tsurf = catstruct(...
    inputs.avhrr_12345, lat, angles, tsurfcfsr);

% inputs.vis_angles = catstruct(...
%     inputs.vis, ...
%     struct(...
%         'B_SZA', struct(...
%             'transform', @cosd, ...
%             'invtransform', @acosd, ...
%             'lims', [0, 180])));

inputs.avhrr_345 = struct(...
        'MEAN_AVHRR_Y', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [0, 400], ...
            'chans', [3, 4, 5]));


inputs.avhrr_45 = struct(...
        'MEAN_AVHRR_Y', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [100, 400], ...
            'chans', [4, 5]));

avhf = {'123', '1234', '12345', '345', '45'};
for i = 1:length(avhf)
    a = sprintf('avhrr_%s', avhf{i});
    inputs.(a).MEAN_AVHRR_Y.stored = inputs.avhrr_12.MEAN_AVHRR_Y.stored;
end

%inputs.avhrr_45_tsurf = catstruct(inputs.avhrr_45, tsurf);
inputs.avhrr_45_tsurf = catstruct(inputs.avhrr_45, tsurfcfsr);
        
inputs.avhrr_45_lat = catstruct(inputs.avhrr_45, lat);

inputs.avhrr_45_angles = catstruct(inputs.avhrr_45, angles);
        
%inputs.avhrr_45_lat_tsurf = catstruct(inputs.avhrr_45, lat, tsurf);
inputs.avhrr_45_lat_tsurf = catstruct(inputs.avhrr_45, lat, tsurfcfsr);

%inputs.avhrr_45_lat_tsurf_angles = catstruct(inputs.avhrr_45, lat, tsurf, angles);
inputs.avhrr_45_lat_tsurf_angles = catstruct(inputs.avhrr_45, lat, tsurfcfsr, angles);

inputs.dir = struct(...
        'dt', struct(...
            'process', ...
                {{@(bt, surf) bsxfun(@minus, bt, surf), ...
                 struct(...
                    'MEAN_AVHRR_Y', ...
                        struct(...
                            'transform', @(x)x, ...
                            'invtransform', @(x)x, ...
                            'lims', [100 400], ...
                            'chans', [4, 5]), ...
                    'MEAN_ECMWF_Skin_temperature', ...
                        struct(...
                            'transform', @(x)x, ...
                            'invtransform', @(x)x, ...
                            'lims', [100 400]))}}));
                        
%             'transform', @(x)x, ...
%             'invtransform', @(x)x, ...
%             'lims', [-100, 400], ...
%             'chans', [4, 5]));
        
% inputs.lmw = struct(...
%         'B_BT', struct(...
%             'transform', @(x)x, ...
%             'invtransform', @(x)x, ...
%             'lims', [100 400], ...
%             'chans', 1));

inputs.mhs_345 = struct(...
        'B_BT', struct(...
            'transform', @(x)x, ...
            'invtransform', @(x)x, ...
            'lims', [100 400], ...
            'chans', 3:5, ...
            'net', 'fit')); % only for the fitnet

inputs.mhs_345.B_BT.stored.realname = 'MHS';
inputs.mhs_345.B_BT.stored.type = 'float';
inputs.mhs_345.B_BT.stored.atts.long_name = 'MHS radiance, channels 3-5';
inputs.mhs_345.B_BT.stored.atts.valid_range = [100, 400];
inputs.mhs_345.B_BT.stored.atts.units = 'K';
inputs.mhs_345.B_BT.stored.atts.origin = 'Copied from MHS L1 data from NOAA CLASS archive';
inputs.mhs_345.B_BT.stored.dims = {'MHS_HUM_CHANS', 3};
%
inputs.mhs_345_only = inputs.mhs_345;
inputs.mhs_345_only.B_BT = rmfield(inputs.mhs_345_only.B_BT, 'net');
% 
inputs.mhs_345_lat = catstruct(inputs.mhs_345_only, lat);

inputs.mhs_345_angles = catstruct(inputs.mhs_345_only, angles);

inputs.mhs_345_lat_angles = catstruct(inputs.mhs_345_only, lat, angles);

%inputs.vis_dir = catstruct(inputs.vis, inputs.dir);

%inputs.vis_dir_ang = catstruct(inputs.vis_dir, angles);

inputs.avhrr_12_mhs_345 = catstruct(inputs.avhrr_12, inputs.mhs_345);

inputs.avhrr_12_mhs_345_angles = catstruct(inputs.avhrr_12_mhs_345, angles);

inputs.avhrr_12_mhs_345_lat_angles = catstruct(inputs.avhrr_12_mhs_345, angles, lat);

inputs.avhrr_12_mhs_345_lat_angles_tsurf = catstruct(inputs.avhrr_12_mhs_345, angles, lat, tsurf);

inputs.avhrr_45_mhs_345 = catstruct(inputs.avhrr_45, inputs.mhs_345);

inputs.avhrr_45_mhs_345_lat = catstruct(inputs.avhrr_45_mhs_345, lat);

inputs.avhrr_45_mhs_345_lat_angles = catstruct(inputs.avhrr_45_mhs_345, lat, angles);

inputs.avhrr_45_mhs_345_lat_angles_tsurf = catstruct(inputs.avhrr_45_mhs_345, lat, angles, tsurfcfsr);

%inputs.dir_mw = catstruct(inputs.dir, inputs.mw);

inputs.avhrr_1245_mhs_345 = catstruct(inputs.avhrr_1245, inputs.mhs_345);

inputs.avhrr_1245_mhs_345_angles = catstruct(inputs.avhrr_1245_mhs_345, angles);

inputs.avhrr_1245_mhs_345_lat_angles = catstruct(inputs.avhrr_1245_mhs_345, lat, angles);

inputs.avhrr_1245_mhs_345_lat_angles_tsurf = catstruct(inputs.avhrr_1245_mhs_345, lat, angles, tsurfcfsr);

inputs.avhrr_12345_mhs_345 = catstruct(inputs.avhrr_12345, inputs.mhs_345);

inputs.avhrr_12345_mhs_345_angles = catstruct(inputs.avhrr_12345_mhs_345, angles);

inputs.avhrr_12345_mhs_345_lat_angles = catstruct(inputs.avhrr_12345_mhs_345_angles, lat);

inputs.avhrr_12345_mhs_345_lat_angles_elevn = catstruct(inputs.avhrr_12345_mhs_345_angles, lat, elev);

%inputs.avhrr_12345_mhs_345_lat_angles_elevn_dz = catstruct(inputs.avhrr_12345_mhs_345_angles, lat, elev, dz);

inputs.avhrr_12345_mhs_345_lat_angles_tsurf = catstruct(inputs.avhrr_12345_mhs_345_angles, lat, tsurfcfsr);

inputs.avhrr_12345_lat_angles_elevn = catstruct(inputs.avhrr_12345, lat, angles, elev);

inputs.avhrr_12345_lat_lza_sza_saa_elevn = catstruct(inputs.avhrr_12345, lat, getfields(angles, 'B_LZA', 'B_SZA', 'B_SAA'), elev);

%inputs.avhrr_12345_lat_angles_elevn_dz = catstruct(inputs.avhrr_12345, lat, angles, elev, dz);

%inputs.avhrr_12345_angles_elevn_tsurfcfsr = catstruct(inputs.avhrr_12345, angles, elev, tsurfcfsr);

%inputs.avhrr_345_angles_elevn_tsurfcfsr = catstruct(inputs.avhrr_12345, angles, elev, tsurfcfsr);

inputs.avhrr_12345_mhs_345_lat_angles_elevn_tsurfcfsr = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, lat, angles, elev, tsurfcfsr);
inputs.avhrr_12345_mhs_345_angles_elevn_tsurfcfsr = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, angles, elev, tsurfcfsr);
inputs.avhrr_12345_mhs_345_lat_angles_elevn = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, lat, angles, elev);
inputs.avhrr_12345_mhs_345_angles_elevn = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, angles, elev);
inputs.avhrr_345_mhs_345_lat_angles_elevn_tsurfcfsr = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, lat, angles, elev, tsurfcfsr);
inputs.avhrr_345_mhs_345_angles_elevn_tsurfcfsr = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, angles, elev, tsurfcfsr);
inputs.avhrr_345_mhs_345_lat_angles_elevn = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, lat, angles, elev);
inputs.avhrr_345_mhs_345_angles_elevn = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, angles, elev);
inputs.avhrr_12345_mhs_345_lat_angles_tsurfcfsr = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, lat, angles, tsurfcfsr);
inputs.avhrr_12345_mhs_345_angles_tsurfcfsr = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, angles, tsurfcfsr);
inputs.avhrr_12345_mhs_345_lat_angles = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, lat, angles);
inputs.avhrr_12345_mhs_345_angles = ...
    catstruct(inputs.avhrr_12345, inputs.mhs_345, angles);
inputs.avhrr_345_mhs_345_lat_angles_tsurfcfsr = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, lat, angles, tsurfcfsr);
inputs.avhrr_345_mhs_345_angles_tsurfcfsr = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, angles, tsurfcfsr);
inputs.avhrr_345_mhs_345_lat_angles = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, lat, angles);
inputs.avhrr_345_mhs_345_angles = ...
    catstruct(inputs.avhrr_345, inputs.mhs_345, angles);

%inputs.vis_dir_mw = catstruct(inputs.vis_dir, inputs.mw);

%inputs.vis_dir_mw_ang = catstruct(inputs.vis_dir_mw, angles);


collocation_filters = {...
    {@(avhrr, surf) avhrr(:, 5) - surf < -20, ...
    {'MEAN_AVHRR_Y', 'MEAN_ECMWF_Skin_temperature'}}, ...
    {@(avhrr) avhrr(:, 5) - avhrr(:, 4) > -2, ...
    {'MEAN_AVHRR_Y'}}, ...
    };

%            'MEAN_ECMWF_Skin_temperature', struct(...
%                'transform', @(x)x, ...
%                'invtransform', @(x)x, ...
%                'lims', [100 400]), ...
%            'B_SZA', struct(...
%                'transform', @sind, ...
%                'invtransform', @asind, ...
%                'lims', [0 80]));
%        inputs = struct(...
%            'B_BT', struct(...
%                 'transform', @(x)x, ...
%                 'invtransform', @(x)x, ...
%                 'lims', [100 400], ...
%                 'chans', [3, 4, 5]));
            
% figure(1);
% clf();
% hold all;
% grid on;

set(0,'DefaultAxesLineStyleOrder',{'-','--', '-.', ':'});
%set(0,'DefaultFigureWindowStyle','docked');

filters.all = {};
filters.dsurf = collocation_filters(1);
filters.dchan = collocation_filters(2);
filters.dsurf_dchan = collocation_filters;
% using https://en.wikipedia.org/wiki/Extreme_points_of_Greenland#Extreme_points
filters.NoGreenland = {{@(lat, lon) ~((lat > 59.75) & (lat < 83.5) & (lon > -73) & (lon < -11)), {'LAT1', 'LON1'}}};

%inps = fieldnames(inputs);
%inps = {'vis', 'avhrr_123', 'vis_angles', 'avhrr_123_angles'};
%inps = {'vis', 'avhrr_123', 'vis_angles', 'avhrr_123_angles', 'ir', 'vis_ir', 'vis_ir_ang', 'avhrr_12345', 'avhrr_12345_angles'};
%inps = {'vis', 'vis_angles', 'vis_ir', 'vis_ir_ang', 'vis_mw', 'vis_mw_ang', 'vis_ir_mw', 'vis_ir_mw_ang'};
%inps = {'vis_angles', 'ir', 'dir', 'mw', 'vis_ir_ang', 'ir_mw', 'vis_mw_ang', 'vis_ir_mw_ang'};
%inps = {'avhrr_12', 'avhrr_123', 'avhrr_45', 'avhrr_1245', 'avhrr_12345'};
%inps = {'avhrr_12_angles', 'avhrr_45', 'avhrr_45_tsurf', 'dir', ...
%    'avhrr_12345_angles', 'avhrr_12_mhs_345_angles', ...
%    'avhrr_45_mhs_345', 'avhrr_45_mhs_345_lat', 'avhrr_12345_mhs_345_angles', ...
%    'mhs_345', 'mhs_345_lat', 'avhrr_12_lat_angles', ...
%    'avhrr_12345_lat_angles', 'avhrr_12345_mhs_345_lat_angles'};
%inps = {'avhrr_12', 'avhrr_12_angles', 'avhrr_12_lat', 'avhrr_12_lat_angles'};
%filts = fieldnames(filters);
%filts = {'all', 'NoGreenland'};
filts = {'all'};%, 'NoGreenland'};

cloud_filters = figure('visible', 'off');
hold all;
grid on;
cloud_filters_leg = {};

figs_inps_retr = struct();
leg_inps_retr = struct();
figs_filts_retr = struct();
leg_filts_retr = struct();

D = datasets();
%pcd = PersistentCachedData('/scratch/uni/u237/users/gholl/cache');
pcd = PersistentCachedData('/scratch/uni/u237/users/gholl/cache');
for i = 1:length(inps)
    inp = inps{i};
    figs_inps_retr.(inp) = figure('visible', 'off');
    hold all;
    grid on;
    
    for j = 1:length(filts)
        filt = filts{j};
    
        if homog
            nm = [inp '_' filt '_' limname '_h'];
        else
            nm = [inp '_' filt '_' limname];
        end
        
        if isfield(D, nm)
            datasets('delete', D.(nm));
        end
        
        ccn = CollocatedNNIWP('name', nm, ...
            'inputs', inputs.(inp), ...
            'collocation_filters', filters.(filt), ...
            'collocation_limits', collocation_limits, ...
            'version', '0.7', ...
            'homog', homog);
        
        % FIXME: disable this one later
        %ccn.overwrite = true;
        
        ccn.getdata();
%        if isfield(ccn.inputs, 'CFSR_Skin_temperature')
%            ccn.overwrite = true;
%        end
        ccn.make_and_train_network(0, no_nets);
        % plot ALL performances for the various networks --> some kind of
        % error
        [best, allperf] = ccn.bestnet(0);
        N = length(allperf);
        f_nets = figure('visible', 'off');
        hold all; grid on;
        for k = 1:N
            plot(allperf{k}{1}, allperf{k}{2}.rel_medfracerr*100);
        end
        title(['Net var ' strrep(nm, '_', ' ')]);
        xlabel('IWP [g/m^2]');
        ylabel('med frac err');
        set(gca(), 'XScale', 'log');
        xlim([1e-1, 1e4]);
        ylim(ylims);
        save_figure_multi(f_nets, ...
            fullfile(cscol('plot_base'), ['allnet_h' num2str(homog) '_' nm '_log']), ...
            'png', 'eps', 'fig');
        % and in a table
        pp = arrayfun(@(i) allperf{i}{2}.rel_medfracerr, 1:5, 'UniformOutput', false);
        M = [allperf{1}{1}, horzcat(pp{:})];
        save(fullfile(cscol('plotdata_base'), ['allnet_h' num2str(homog) '_' nm '_log.dat']), '-ascii', 'M');

        % and consider performance for best net
        [xx, perf] = ccn.get_performance(0);
        
        % plot false positive / false negative for colud filter
        [bn, fp, fn, fp_med_ref, fp_med_retr, fn_med_ref, fn_med_retr] = ccn.test_patternnet_performance();
        figure('visible', 'off');
        plot(bn, fp, bn, fn, bn, fp+fn);
        legend({'false positive', 'false negative', 'false total'}, ...
            'Location', 'NorthWest');
        title(['Cloud filter performance ' strrep(nm, '_', ' ')]);
        ylabel('Error fraction');
        xlabel('Probability cutoff');
        save_figure_multi(gcf(), fullfile(cscol('plot_base'), ['cloud_filter_' nm]), 'png', 'eps', 'fig');
        
        % store to file
        M = [bn; fp; fn; fp+fn; fp_med_ref; fp_med_retr; fn_med_ref; fn_med_retr].'; %#ok<NASGU>
        save(fullfile(cscol('plotdata_base'), ['cloud_filter_' num2str(homog) '_' ccn.name '.dat']), '-ascii', 'M');
        
        % plot to filters plot for technique
        figure(figs_inps_retr.(inp)); set(gcf(), 'Visible', 'off');
        plot(xx, 100*perf.rel_medfracerr);
        leg_inps_retr.(inp){j} = strrep(filt, '_', ' ');
        
        % plot cloud filter performance overall
        figure(cloud_filters); set(gcf(), 'Visible', 'off');
        plot(bn, fp+fn);
        cloud_filters_leg = [cloud_filters_leg, strrep(nm, '_', ' ')]; %#ok<AGROW>
        
        % plot to techniques plot for this filter
        if ~isfield(figs_filts_retr, filt)
            figs_filts_retr.(filt) = figure('visible', 'off');
            hold all;
            grid on;
        end
        figure(figs_filts_retr.(filt)); set(gcf(), 'Visible', 'off');
        plot(xx, 100*perf.rel_medfracerr);
        leg_filts_retr.(filt){i} = strrep(inp, '_', ' ');
        
        % store performance to file
        M = [xx, 100*perf.rel_medfracerr]; %#ok<NASGU>
        save(fullfile(cscol('plotdata_base'), [ccn.name '.dat']), '-ascii', 'M');
        
        % plot scatter density plot
        f = figure('visible', 'off');
        [x, y_ref, y_retr] = ccn.evaluate_test_data(0);
        [h_ax, h_cb] = scatter_density_plot(...
            y_ref, ...
            y_retr, ...
            struct(...
                'trans', @(x)log10(x), ...
                'invtrans', @(x)10.^x, ...
                'axprops', struct(...
                    'xscale', 'log', ...
                    'yscale', 'log', ...
                    'xlim', [1e-1, 1e4], ...
                    'ylim', [1e-1, 1e4]), ...
                'medprops', struct(...
                    'LineWidth', 3, ...
                    'LineStyle', '-', ...
                    'Color', [.5, .5, .5]), ...
                'diagonal', struct(...
                    'LineWidth', 1, ...
                    'LineStyle', '-', ...
                    'Color', [0, 0, 0])));
        xlabel(h_ax, '2C-ICE testing IWP [g/m^2]');
        ylabel(h_ax, 'Passive retrieved IWP [g/m^2]');
        title(['Retrieved vs. reference, ' strrep(nm, '_', ' ')]);
        bn = ['sdp_h' num2str(homog) '_' nm '_2cice'];
        save_figure_multi(f, fullfile(cscol('plot_base'), bn), 'png', 'eps', 'fig');
        
        % second sdp version without axis etc., for use inside pgfplots
        axis(gca(), 'off');
        set(h_cb, 'Visible', 'off');
        title('');
        save_figure_multi(f, fullfile(cscol('plot_base'), [bn '_bare']), 'png', 'eps', 'fig');
        f2 = figure('Visible', 'off');
        set(h_cb, 'YTick', []);
        set(h_cb, 'Visible', 'on');
        copyobj(h_cb, f2);
        set(f2, 'Colormap', get(f, 'Colormap'));
        save_figure_multi(f2, fullfile(cscol('plot_base'), [bn '_bare_legend']), 'png', 'eps', 'fig');
        
        if process_month
            ccn.basedir = fullfile('/storage3/user_data/gerrit/products/experimental', ccn.name);
            if ~exist(ccn.basedir, 'file')
                mkdir(ccn.basedir);
            end
            ccn.subdir = '$YEAR4/$MONTH/$DAY';
            ccn.filename = 'granule_$VERSION_$SAT_$HOUR_$MINUTE.nc.gz';
            ccn.altbasedir = '/storage3/user_data/gerrit/products/passive_syn_iwp';
%            ccn.altfilename = '$SAT_$HOUR_$MINUTE.nc.gz';
            ccn.altfilename = 'ccniwp_$VERSION_$SAT_$HOUR_$MINUTE.nc.gz';
            ccn.reader_processor = @ccniwp_proc;
            logtext(atmlab('OUT'), 'Processing data\n');
            
            start = [2007, 1, 1];
            finish = [2007, 12, 31];
            sat = 'noaa18';
            if ~isempty(ccn.find_granules_by_date(start, sat)) && ~isempty(ccn.find_granules_by_date(finish, sat))
                logtext(atmlab('OUT'), 'It seems I''m already done.\n');
            else
                ccn.retrieve_and_store_period(start, finish, sat);
            end
            
            % and plot it
            %gridmean = pcd.evaluate(1, @ccn.level3, start, finish, sat, {'IWP'}, 1, @(X)(X(X<5e5)), 'EXTRA', ccn.name);
            S = gmt_get_struct_gridded_map_iwp;
            S.header = [ccn.name datestr(datenum(start), 'YYYY-mm-dd'), '--', datestr(datenum(finish), 'YYYY-mm-dd'), ' v' ccn.version];
            S.header_fontsize = '18p';
            S.filename = ['map_gridded_IWP_' ccn.name '_' datestr(datenum(start), 'YYYYmmdd'), '-', datestr(datenum(finish), 'YYYYmmdd'), 'v' strrep(ccn.version, '.', '_') '.pdf'];
            S.outdir = cscol('plot_base');
            if exist(fullfile(S.outdir, S.filename), 'file')
                logtext(atmlab('OUT'), 'Already exists: %s.  Skipping.\n', fullfile(S.outdir, S.filename));
            else
                S.data = pcd.evaluate(1, @ccn.level3, start, finish, sat, {'IWP'}, 1, @(X)(X(X<5e5)), 'EXTRA', ccn.name);            
                gmt_plot(S);
            end

        end
        
        datasets('delete', ccn);
    end
end

% finalise all figures

allfigs_retr = vec2row(cell2mat(struct2cell(catstruct(figs_filts_retr, figs_inps_retr))));
for fig = allfigs_retr
    figure(fig);set(fig, 'Visible', 'off');
    xlabel('IWP [g/m^2]');
    ylabel('IWP difference to reference [%]');
    ylim(ylims);
    xlim([0, 3000]);
end

base_plotfile = fullfile(cscol('plot_base'), 'comparison_collocated_iwp_products');
for i = 1:length(inps)
    inp = inps{i};    
    figure(figs_inps_retr.(inp)); set(gcf(), 'Visible', 'off');
    legend(leg_inps_retr.(inp){:}, 'Location', 'SouthWest');
    title([strrep(inp, '_', ' '), ' ', limname]);
    nm = [base_plotfile, '_h', num2str(homog), '_', limname, '_', inp, '_', strjoin(filts, ','), '_filters'];
    save_figure_multi(figs_inps_retr.(inp), ...
        [nm '_lin'], ...
        'png', 'eps', 'fig');
    set(gca(), 'XScale', 'log');
    xlim([1e-1, 1e4]);
    save_figure_multi(figs_inps_retr.(inp), ...
        [nm '_log'], ...
        'png', 'eps', 'fig');
end

for i = 1:length(filts)
    filt = filts{i};
    figure(figs_filts_retr.(filt)); set(gcf(), 'Visible', 'off');
    legend(leg_filts_retr.(filt){:}, 'Location', 'SouthWest');
    title([strrep(filt, '_', ' '), ' ', limname]);
    nm = [base_plotfile, '_h', num2str(homog), '_', limname, '_', filt, '_', strjoin(vec2row(inps), ','), '_techniques'];
    % check if name too long?
    [~, bnm] = fileparts(nm);
    if length(bnm) >= 220 % 252 because '.png' or '.eps' will be added still
        logtext(atmlab('OUT'), 'Name too long, shortening %s...\n', nm);
%        nm = [bnm(1:100), '........', bnm(end-100:end)];
%        nm = strrep(nm, 'avhrr_', 'a');
        [dirname, bnm] = fileparts(nm);
        bnm = [bnm(1:100) '...' num2str(length(inps)) '...' bnm(end-100:end)];
        nm = fullfile(dirname, bnm);
%        if length(bnm) >= 250
%            nm = strrep(nm, 'mhs_', 'm');
%            nm = strrep(nm, 'angles', 'g');
%            [~, bnm] = fileparts(nm);
%            if length(bnm) >= 250
%                [core, bnm] = fileparts(nm);
%                bnm = regexprep(bnm, '\d*', '');
%                bnm = strrep(bnm, 'comparison_collocated_iwp_products_', 'ccip_');
%                nm = fullfile(core, bnm);
%                %nm = regexprep(nm, '\d*', '');
%                if length(bnm) >= 250
%                    [core, bnm] = fileparts(nm);
%                    bnm = [bnm(1:100) '.........' bnm(end-100:end)];
%                    nm = fullfile(core, bnm);
%                end
%            end
%        end
    end
    save_figure_multi(figs_filts_retr.(filt), ...
        [nm '_lin'], ...
        'png', 'eps', 'fig');
    set(gca(), 'XScale', 'log');
    xlim([1e-1, 1e4]);
    save_figure_multi(figs_filts_retr.(filt), ...
        [nm '_log'], ...
        'png', 'eps', 'fig');    
end

figure(cloud_filters); set(gcf(), 'Visible', 'off');
xlabel('NN prob. cutoff');
ylabel('Total error rate (fp + fn)');
title(['Cloud filter quality, ', strrep(limname, '_', ' ')]);
legend(cloud_filters_leg{:}, 'Location', 'NorthWest');
save_figure_multi(cloud_filters, fullfile(cscol('plot_base'), ['cloud_filters_h', num2str(homog), '_' limname]), 'png', 'eps', 'fig');
% 
%     allnames = {[inp '_all'], [inp '_filt_dsurf'], [inp '_filt_dchan'], [inp '_filt_dsurf_dchan']};
%     nm = [inp '_all'];
%     ccn.(nm) = CollocatedNNIWP('name', nm, ...
%         'inputs', inputs.(inp), ...
%         'collocation_filters', {}, ...
%         'collocation_limits', collocation_limits);
%     nm = [inp '_filt_dsurf'];
%     ccn.(nm) = CollocatedNNIWP('name', nm, ...
%         'inputs', inputs.(inp), ...
%         'collocation_filters', collocation_filters(1), ...
%         'collocation_limits', collocation_limits);
%     nm = [inp '_filt_dchan'];
%     ccn.(nm) = CollocatedNNIWP('name', nm, ...
%         'inputs', inputs.(inp), ...
%         'collocation_filters', collocation_filters(2), ...
%         'collocation_limits', collocation_limits);
%     nm = [inp '_filt_dsurf_dchan'];
%     ccn.(nm) = CollocatedNNIWP('name', nm, ...
%         'inputs', inputs.(inp), ...
%         'collocation_filters', collocation_filters, ...
%         'collocation_limits', collocation_limits);
%     
%     for k = 1:length(allnames)
%         ff = allnames{k};
%         if isfield(figs_filts.(ff))
%             figure(figs_filts.(ff));
%         else
%             figs_filts.(ff) = figure();
%             hold all;
%             grid on;
%         end
%         % for each combination technique / filter, plot both in
%         % all-techniques plot, and in all-filters plot
%         ccn.(ff).getdata();
%         ccn.(ff).make_and_train_network(1, 5);
%         ccn.(ff).plot_performance_with_requirements(1);
%         
%         [xx, perf.(ff)] = ccn.(ff{1}).get_performance(1);
%         
%         % plot to all-techniques plot                
%         figure(figs_inps.(inp));
%         plot(xx, 100*perf.(ff).rel_medfracerr);
%         leg_inps.(inp){i} = strrep(ff{1}, '_', ' ');
%     end
% end
% figure(1);
% xlabel('IWP [g/m^2]');
% ylabel('IWP difference to reference [%]');
% ylim([0, 100]);
% xlim([0, 3000]);
% %nms = cellfun(@(s) strrep(s, '_', ' '), allf, 'UniformOutput', false);
% %legend(nms{:});
% legend(leg{:}, 'Location', 'SouthWest');
% save_figure_multi(1, fullfile(cscol('plot_base'), 'comparison_collocated_iwp_products_lin'), 'png', 'eps', 'fig');
% set(gca(), 'XScale', 'log');
% xlim([5e0, 1e4]);
% save_figure_multi(1, fullfile(cscol('plot_base'), 'comparison_collocated_iwp_products_log'), 'png', 'eps', 'fig');
end

function T = get_cfsr_skin_temperature(lat, lon, time)
D = datasets();
[ye, mo, da, ho, mi, se] = unixsecs2date(time);
dv = [ye, mo, da, ho, mi, se];
[dv_sorted, I] = sortrows(dv);
cf = D.CFSR;
pcd = PersistentCachedData('/scratch/uni/u237/users/gholl/cache');
S = pcd.evaluate(1, @cf.read_from_grid, lat(I), lon(I), dv(I, :), {'TMP_L1'});
T(I) = S.TMP_L1;
end