Home > atmlab > retrieval > spareice > InstrumentVisualiser.m

InstrumentVisualiser

PURPOSE ^

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

InstrumentVisualiser.m

SOURCE CODE ^

0001 classdef InstrumentVisualiser < handle
0002     
0003     % mixin with functionality to visualise performance etc.
0004     %
0005     % InstrumentVisualiser Methods:
0006     %
0007     %   plot_performance_with_requirements -
0008     %   plot_channels_bt -
0009     %
0010     % See also: IceMusicNNIWP, RetrievalDatabaseProducte, NNTrainedProduct
0011 
0012     
0013     properties (Abstract)
0014         opts
0015     end
0016     
0017 %     properties (Constant, Abstract)
0018 %         detect_limit
0019 %         accuracy
0020 %     end
0021     
0022     methods (Abstract)
0023         evaluate_test_data(self, noise_level)
0024         get_performance(self, noise_level)
0025     end    
0026     
0027     methods
0028         %         function self = InstrumentVisualiser(varargin)
0029         %             self = self@DBIWP(varargin{:});
0030         %         end
0031         
0032         
0033         function plot_performance_with_requirements(self, noise_level, varargin)
0034             
0035             % Plots performance with requirements.
0036             %
0037             % Makes two plots; absolute and relative scales.
0038             %
0039             % Single argument is basename for output files.
0040             %
0041             % if given, detect_limit and accuracy are structures describing
0042             % the detection limit and accuracy
0043             
0044             [detect_limit, accuracy] = optargs(varargin, {[], []});
0045             doreq = nargin > 2;
0046             
0047 %             [x, y] = self.evaluate_test_data(0.1);
0048 %             filt = (x>0) & (y>0) & isfinite(x) & isfinite(y);
0049 %
0050 %             x = x(filt);
0051 %             y = y(filt);
0052 %
0053 %             % calculate different ways to get performance
0054 %             abs_medad = cellfun(@(yy) mad(yy, 1), bin(log10(x), y, self.opts.bins));
0055 %             abs_meanad = cellfun(@(yy) mad(yy, 0), bin(log10(x), y, self.opts.bins));
0056 %             abs_rmeanse = cellfun(@(dy) rms(dy), bin(log10(x), y-x, self.opts.bins));
0057 %             abs_rmedse = cellfun(@(dy) sqrt(median(dy.^2)), bin(log10(x), y-x, self.opts.bins));
0058 %             abs_ir = cellfun(@iqr, bin(log10(x), y, self.opts.bins));
0059 %             rel_medfracerr = cellfun(@median, bin(log10(x), abs(y-x)./x, self.opts.bins));
0060 %             abs_medfracerr = rel_medfracerr .* (10.^self.opts.bins);
0061 %             rel_rmeanse_frac = cellfun(@(dy) rms(dy), bin(log10(x), abs(y-x)./x, self.opts.bins));
0062 %             abs_rmeanse_frac = rel_rmeanse_frac .* (10.^self.opts.bins);
0063 %
0064 %             xx = 10.^(self.opts.bins);
0065             %           self.binned_mad(x(filt), y(filt), self.opts.bins);
0066             
0067             %% make "absolute" plot
0068             
0069 %            [xx, S] = self.get_performance_with_requirements(0.1);
0070             [xx, S] = self.get_performance(noise_level);
0071           
0072             figure();
0073             % plot various ways to visualise error
0074             plot(xx, S.abs_medad, ...
0075                 xx, S.abs_meanad, ...
0076                 xx, S.abs_rmeanse, ...
0077                 xx, S.abs_rmedse, ...
0078                 xx, S.abs_ir, ...
0079                 xx, S.abs_medfracerr, ...
0080                 xx, S.abs_rmeanse_frac);
0081             legend('median absolute deviation', ...
0082                 'mean absolute deviation', ...
0083                 'root mean square error', ...
0084                 'root median square error', ...
0085                 'inter-quantile range', ...
0086                 'median fractional error', ...
0087                 'root mean square fractional error', ...
0088                 'Location', 'NorthWest');
0089             
0090             hold on;
0091             
0092             if doreq
0093                 crit_treshhold = detect_limit.treshhold / accuracy.treshhold;
0094                 crit_target = detect_limit.target / accuracy.target;
0095                 
0096                 % make sure target/treshhold lines are continuous
0097                 xxx = sort([xx; crit_target; crit_treshhold]);
0098                 
0099                 % calculate upper treshhold = 10 g/m^2
0100                 acc_tresh = xxx.*accuracy.treshhold;
0101                 acc_targ = xxx.*accuracy.target;
0102                 
0103                 % plot target and treshhold
0104                 plot(xxx(acc_tresh>=detect_limit.treshhold), acc_tresh(acc_tresh>=detect_limit.treshhold), 'k--')
0105                 plot(xxx(acc_targ>=detect_limit.target), acc_targ(acc_targ>=detect_limit.target), 'k--');
0106                 
0107                 % plot detection limits
0108                 plot([1e-1, crit_treshhold], [detect_limit.treshhold detect_limit.treshhold], 'k--');
0109                 plot([1e-1, crit_target], [detect_limit.target detect_limit.target], 'k--');
0110             end
0111             
0112             % xx*(1+self.accuracy.target)=10
0113             
0114             set(gca(), 'YScale', 'log', 'XScale', 'log', 'XLim', [5*1e-1 1e3], 'YLim', [5*1e-1 1e3])
0115             
0116             title([self.name ' performance, trop. nad. all chans']);
0117             xlabel('IWP [g/m^2]');
0118             ylabel('IWP error [g/m^2]');
0119             grid('on');
0120             
0121             [~, b] = fileparts(self.net.userdata.stored);
0122             save_figure_multi(gcf(), ...
0123                 fullfile(cscol('plot_base'), [b '_abs']), ...
0124                 'png', 'eps', 'fig');
0125             
0126             %% make "relative" plot
0127             
0128             figure();
0129             plot(xx, 100*S.rel_medfracerr);
0130             
0131             % plot target and treshhold
0132             %           plot(xxx(acc_tresh>=self.detect_limit.treshhold), self.detect_limit.treshhold, 'k--')
0133             %           plot(xxx(acc_targ>=self.detect_limit.target), acc_targ(acc_targ>=self.detect_limit.target), 'k--');
0134             
0135             legend('median fractional error');
0136             ylim([0 100]);
0137             xlim([0 3000]);
0138             xlabel('IWP [g/m^2]');
0139             ylabel('IWP error [%]');
0140             title([ self.name ' performance, trop. nad. all chans']);
0141             grid('on');
0142             save_figure_multi(gcf(), ...
0143                               fullfile(cscol('plot_base'), [b '_rel']), ...
0144                               'png', 'eps', 'fig');
0145             
0146         end
0147         
0148         function plot_channels_bt(self)
0149             % for all channels, plot BT vs. log10(IWP)
0150             
0151             ff = self.data(:, self.localcols.IWP)>0;
0152             op.transx = @(x)log10(x);
0153             op.invtransx = @(x)10.^(x);
0154             op.axprops.xscale = 'log';
0155             
0156             for i = self.chans
0157                 figure();
0158                 scatter_density_plot(self.data(ff, self.localcols.IWP), self.data(ff, self.localcols.BT(i)), op);
0159                 xlabel('IWP [g/m^2]');
0160                 ylabel('BT [K]');
0161                 title(sprintf('BT/IWP Ch. %d (%s GHz)', self.chans(i), self.freqs{i}));
0162                 print('-depsc', fullfile(cscol('plot_base'), sprintf('icemusic_btiwp_ch%d.eps', i)));
0163             end
0164         end
0165         
0166         function iwc_profile_errorhist(self, varargin)
0167             % visualise IWC error as a function of height, statistically
0168             
0169             % get data to go into histogram plot
0170             [x, y_ref, y_retr, x_cols, y_cols] = self.evaluate_test_data(varargin{:});
0171             dlogy = log10(y_ref(:, y_cols.IWC)) - log10(y_retr(:, y_cols.IWC));
0172             bins = linspace(-4, 4, 40);
0173             y_axis = linspace(4.5, 15, 21); % FIXME: read from data
0174             [N, X] = arrayfun(@(i) hist(dlogy(:, i), bins), 1:size(dlogy, 2), 'UniformOutput', false);
0175             
0176             % plot everything
0177             pcolor(10.^(X{1}), y_axis, vertcat(N{:}));
0178             shading('flat');
0179             hold('all');
0180             grid('on');
0181             plot(10.^median(dlogy, 1), y_axis, 'k-', 'LineWidth', 2);
0182             set(gca, 'XScale', 'log', ...
0183                      'layer', 'top', ...
0184                      'XLim', [1e-2 1e2]); 
0185             title(strrep([self.name ' IWC error'], '_', ' '));
0186             xlabel('IWC retrieved / reference');
0187             ylabel('Height [km]');
0188             cb = colorbar;
0189             ylabel(cb, 'No. profiles');
0190             [~, b] = fileparts(self.net.userdata.stored);
0191             save_figure_multi(gcf(), fullfile(cscol('plot_base'), [b '_iwchist']), 'eps', 'png', 'fig');
0192 
0193         end
0194     end
0195     
0196 end

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