diagnose_zonal_feedback_lags.m 11.6 KB
Newer Older
Claas Faber's avatar
Claas Faber committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                    %
%   DIAGNOSE FEEDBACK LAGS FOR EACH LONGITUDE        %
%   MAIN: CALCULATION FOR A GIVEN FEEDBACK ELEMENT   %
%   Plotting is done in plot_lagged_feedbacks.m      %
%                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Equatorial indices at lag 0 can be calculated with
% zonal_feedback_strengths_seasonal_regressions.m

diagnose_zonal_feedback_lags_settings;

%  if load_data
%      clear all;
%      diagnose_zonal_feedback_lags_settings;
%  end

disp(sprintf('\n%s %s-%s lag relationship (%s), %d-%d\n',...
        basin_name,upper(reference_name),upper(lag_name),...
        analysis_title,start_year,stop_year));

if load_data
% ------------------ %
% LOAD AND CROP DATA %
% ------------------ %

    disp('loading data');
    data = struct();

    for var_name = {all_vars{:}}

        var_name = char(var_name);

        % file definition
        in_file = strrep(file_dummy,'@DATASET_NAME@',all_data.(var_name));
        in_file = strrep(in_file,'@VAR_NAME@',var_name);

        % actual loading
        data.(var_name) = load(in_file);

    % load all data
    end

    % check that both date vectors are equal
    equal_years      = ( data.(reference_name).date_vector(:,1) == data.(lag_name).date_vector(:,1) );
    equal_months     = ( data.(reference_name).date_vector(:,2) == data.(lag_name).date_vector(:,2) );
    n_reference      = size(data.(reference_name).date_vector,1);
    n_lag            = size(data.(lag_name).date_vector,1);

    equal_lengths = ( n_reference == n_lag );
    equal_entries = ( sum( equal_years & equal_months ) == n_reference );

    if ( ~ ( equal_lengths & equal_entries ) )
        error('Reference and lag data do not have equal time axis.');
    end

% if: load data
end

if calculate_feedback_strengths
% ----------------------------- %
% CALCULATE LAGGED CORRELATIONS %
% ----------------------------- %

    % switch off iteration limit warning for robust regression
    if ( analysis_id == 3 )
        warning('off',warn_id);
    end

    % both variables use the same date vector
    date_vector = data.(reference_name).date_vector;

    % basin zonal extend
    n_lon = length(data.(reference_name).central_longitudes);

    % how many lags in total?
    all_lags    = [-maximum_lag:1:maximum_lag];
    n_lags      = length(all_lags);

    feedback_strengths = zeros(n_lon,12,n_lags);

    if save_lagged_indices

        % saving the lagged time series
        all_reference_indices = zeros(n_lon,12,n_lags,1);
        all_lag_indices       = zeros(n_lon,12,n_lags,1);

        reference_index = zeros(n_lon,12,1);
        lag_index       = zeros(n_lon,12,1);

    end

    % loop over all longitudes, lags, and months
    textprogressbar('Calculating lags: ')

    for c_lon = 1:n_lon
    textprogressbar(c_lon/n_lon*100);

        for ref_month = 1:12

            for c_lag = 1:n_lags

                lag = all_lags(c_lag);

                % ADAPT START AND STOP YEARS?
                % +++++++++++++++++++++++++++

                    % initial state
                    lag_month     = ref_month + lag;
                    negative_wrap = ( lag_month <  1);
                    positive_wrap = ( lag_month > 12 );
                    n_wraps       = 0;

                    % start and stop years for both time series
                    start_ref = start_year;
                    stop_ref  = stop_year;
                    start_lag = start_year;
                    stop_lag  = stop_year;

                    % hard-wire whether to shift up or down
                    % the lag time
                    lag_up   = positive_wrap;
                    lag_down = negative_wrap;

                    % wrap around the lag month
                    while ( negative_wrap | positive_wrap );

                        if negative_wrap
                            lag_month = lag_month + 12;
                        else
                            lag_month = lag_month - 12;
                        end

                        % keep track of how many times you wrapped around
                        % the true month
                        n_wraps = n_wraps + 1;

                        % wrap finished?
                        negative_wrap = ( lag_month <  1);
                        positive_wrap = ( lag_month > 12 );

                    end

                    % match the range where both time series exist
                    % lag_month > 12: lag time series needs to start later
                    % lag_month < 1 : ref time series needs to start later
                    if lag_up
                        stop_ref  = stop_ref  - n_wraps;
                        start_lag = start_lag + n_wraps;
                    elseif lag_down
                        start_ref = start_ref + n_wraps;
                        stop_lag  = stop_lag  - n_wraps;
                    end

                % SELECT CORRECT TIME
                % +++++++++++++++++++

                    % reference
                    is_ref_year  = ( ( date_vector(:,1) >= start_ref ) & ...
                                     ( date_vector(:,1) <= stop_ref ) );
                    is_ref_month = ( date_vector(:,2) == ref_month );
                    is_ref       = ( is_ref_year & is_ref_month );

                    % lag
                    is_lag_year  = ( ( date_vector(:,1) >= start_lag ) & ...
                                     ( date_vector(:,1) <= stop_lag ) );
                    is_lag_month = ( date_vector(:,2) == lag_month );
                    is_lag       = ( is_lag_year & is_lag_month );

                % CALCULATE FEEDBACK STRENGTHS
                % ++++++++++++++++++++++++++++

                    ref_series = squeeze(data.(reference_name).indices(c_lon,is_ref));
                    lag_series = squeeze(data.(lag_name).indices(c_lon,is_lag));

                    % calculate feedback strengths
                    switch analysis_id

                        % correlation
                        case 1

                            % acc
                            current_feedback = corrcoef(ref_series,lag_series);
                            current_feedback = current_feedback(1,2);

                        % normal linear regression
                        case 2

                            regression_parameters = polyfit(ref_series,lag_series,1);
                            current_feedback = regression_parameters(1);

                        % robust regression
                        case 3

                            rregression_parameters = robustfit(ref_series,lag_series);
                            current_feedback = rregression_parameters(2);  % yes, the order is different from polyfit ...

                    % switch: calculate feedback strength
                    end

                    % save feedback strength data
                    feedback_strengths(c_lon,ref_month,c_lag) = current_feedback;

                    % save indices
                    if save_lagged_indices

                        n_series = length(ref_series);

                        all_reference_indices(c_lon,ref_month,c_lag,1:n_series) = ref_series;
                        all_lag_indices(c_lon,ref_month,c_lag,1:n_series)       = lag_series;

                    end

            % for: calculate feedback_strengths for all lags
            end

        % for: calculate feedback_strengths for all longitudes
        end

    % for: calculate feedback_strengths for all months
    end
    textprogressbar('');

    % replace zeros by NaNs
    all_reference_indices(all_reference_indices == 0.) = NaN;
    all_lag_indices(all_lag_indices == 0.)             = NaN;

    % switch rreg warning back on
    if ( analysis_id == 3 )
        warning('on',warn_id);
    end

% --------------------------- %
% FIND LAGS AND SAVE THE DATA %
% --------------------------- %

    % mark peak feedback_strengths
    % find the x-indices of the peak feedback_strengths for each month (i.e. in the y-dimension)
    % no_peak: value, is NaN when all feedback_strengths are insignificant
    % peak_index: index value
    % use the sign-sensitive maximum, not the maximum of the absolute feedback strengths,
    % since a strongly negative (damped) feedback is not what you are interested in here.
    [peak_strengths,peak_index] = max(feedback_strengths,[],3);

    % actual lag
    lag_values = all_lags(peak_index);

    % transpose everything for correct later use in plotting
    lag_values     = lag_values';
    peak_strengths = peak_strengths';

% if: calculate feedback_strengths
end

if save_lagged_indices
% --------- %
% SAVE DATA %
% --------- %

    disp('saving data');

    n_series           = size(all_reference_indices,4);
    central_longitudes = data.(reference_name).central_longitudes;

    % for each month and longitude, select the correct pair of reference
    % and lagged variable
    for c_lon = 1:n_lon

        for ref_month = 1:12

            use_lag = peak_index(c_lon,ref_month);

            reference_index(c_lon,ref_month,1:n_series) = ...
                    squeeze(all_reference_indices(c_lon,ref_month,use_lag,:));
            lag_index(c_lon,ref_month,1:n_series) = ...
                    squeeze(all_lag_indices(c_lon,ref_month,use_lag,:));

        % for: select lagged time series for all longitudes
        end

    % for: select lagged time series for all longitudes
    end

    % actual saving
    % create a readme
    README = sprintf(['VARIABLE LIST\n\n',...
                      'reference_name:     reference variable\n',...
                      'lag_name:           variable that is lagged relative to the reference\n',...
                      'all_data:           used datasets\n',...
                      'start_year\n',...
                      'stop_year\n',...
                      'central_longitudes: longitudes along the equator\n',...
                      'date_vector\n',...
                      'all_lags:           list of all lags that have been diagnosed\n',...
                      'maximum_lag\n\n',...
                      'feedback_strengths: for lon x month x lag\n',...
                      'peak_strengths:     maximum correlation along lag-dimension, month x lon\n',...
                      'peak_index:         index of the maximum correlation, month x lon\n',...
                      'lag_values:         actual lag, month x lon\n\n',...
                      'indices:            contains for each longitude and each month\n',...
                      '                    the pair of reference and lagged variable that corresponds\n',...
                      '                    to the peak correlation. lon x month x time\n',...
                      'all_indices:        all pairs of reference/lag time series, for all lags.\n',...
                      '                    lon x month x lag x time']);

    % put data into the same order as in the other index-file
    indices     = struct(reference_name,reference_index,...
                         lag_name,lag_index);
    all_indices = struct(reference_name,all_reference_indices,...
                         lag_name,all_lag_indices);

    % actual saving
    save(out_file,'reference_name','lag_name','all_data','start_year','stop_year',...
                  'central_longitudes','date_vector','all_lags','maximum_lag',...
                  'feedback_strengths','peak_strengths','peak_index','lag_values',...
                  'indices','all_indices',...
                  'README');

% if: save data
end

% load previous data
if ~perform_calculations
    load(out_file); 
end

% ------------- %
% PRODUCE PLOTS %
% ------------- %

    diagnose_zonal_feedback_lags_plots;