%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Converting harmonised meta-data for tu, ts, uu, us and fates into mass
%%% flows using WOT1.2 POM and unit weights
%%%
%%% Dmitry Yumashev, 12/12/2019
%%% Dmitry Yumashev, 21/12/2019
%%% Dmitry Yumashev, 22/12/2019
%%% Dmitry Yumashev, 23/12/2019
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all %#ok<*CLALL>

%%% general settings

step_fitting = 0.5; % required new timestep after re-binning

scaling_factor_samples = 1e+6; % this needs to be large to improve the accuracy of fitting

%%% options for the fmincon solver (nonlinear programming)

options = optimoptions(@fmincon,'Algorithm','interior-point','OptimalityTolerance',1e-10,'StepTolerance',1e-12,'MaxFunctionEvaluations',1e+4,'MaxIterations',1e+4,'PlotFcn','optimplotx','UseParallel',1);


%%%

current_year = 2019;

%%% new timespan for the time in use and time in storage distributions

new_pdf_timestep = 1;
new_pdf_timespan_start = 0;
new_pdf_timespan_end = 40;
total_new_pdf_timesteps = 1 + floor((new_pdf_timespan_end - new_pdf_timespan_start) / new_pdf_timestep);

new_pdf_timespan = transpose(linspace(new_pdf_timespan_start,new_pdf_timespan_end,total_new_pdf_timesteps)); % include t=0 in the new entries; ANNUAL time step!!!
% NOTE: transpose to match with the structure of the time distribution data


%%% new range for the time series results (historic POM and WG)

new_timeseries_timestep = 1;
new_timeseries_timerange_start = 1950;
new_timeseries_timerange_end = 2021;

total_new_timeseries_timesteps = 1 + floor((new_timeseries_timerange_end - new_timeseries_timerange_start) / new_timeseries_timestep);

new_timeseries_timerange = linspace(new_timeseries_timerange_start,new_timeseries_timerange_end,total_new_timeseries_timesteps); % NOTE: do not transpose 


%%% input and output files

filename_output_prelim = 'flows_and_stocks_meta_prelim.xlsx';
filename_output_advanced = 'flows_and_stocks_meta_advanced.xlsx';

%%% input files

% time
folder_time = 'time';
filename_time = [folder_time '/time_final_pdfs.csv'];
filename_time_availability = [folder_time '/UNU_cat_timedata.csv'];

% stock
folder_stock = 'stock';
filename_stock = [folder_stock '/stock_final_pdfs.csv'];
filename_stock_availability = [folder_stock '/UNU_cat_stockdata.csv'];

% fate
folder_fate = 'fate';
filename_fate = [folder_fate '/final_fate.csv'];

% WOT1.2 UK POM and WG
folder_wot = 'wot'; 
filename_wot_pom = [folder_wot '/Summary_POM_from_WOT12.xlsx'];
filename_wot_wg = [folder_wot '/Summary_WG_from_WOT12.xlsx'];

% ONS UK households data
folder_ons = 'ons';
filename_ons = [folder_ons '/ons_households_data.xlsx'];


%%% read the data sample

% time
[nums_time, ~, raw_time] = xlsread(filename_time); 
nums_time_avail = xlsread(filename_time_availability);

% stock
[nums_stock, ~, raw_stock] = xlsread(filename_stock);
nums_stock_avail = xlsread(filename_stock_availability);

% fate
[nums_fate, ~, raw_fate] = xlsread(filename_fate); 

% WOT

% POM units
sheet_pom_units = 'GBR WOT1.2 POM Units';
range_pom_units = 'A2:AS57';
[nums_wot_pom_units,~,raw_wot_pom_units] = xlsread(filename_wot_pom, sheet_pom_units, range_pom_units);

% POM weight per unit
sheet_pom_weightpu = 'GBR WOT1.2 Weight Per Unit';
range_pom_weightpu = 'A2:AS57';
nums_wot_pom_weightpu = xlsread(filename_wot_pom, sheet_pom_weightpu, range_pom_weightpu);

% WG tonnes
sheet_wg_tonnes = 'GBR WOT1.2 WG Tonnes';
range_wg_tonnes = 'A2:AS57';
nums_wot_wg_tonnes = xlsread(filename_wot_wg, sheet_wg_tonnes, range_wg_tonnes);

% ONS

sheet_ons = 'estimates';
nums_ons = xlsread(filename_ons, sheet_ons);


%%% extract unu keys for the available datasets

% time

unu_keys_time = nums_time_avail(:,2);
%
flags_time_use = nums_time_avail(:,3);
flags_time_storage = nums_time_avail(:,4);
%
indx_time_use = find(flags_time_use==1); % flag=1 means data is available
unu_keys_time_use = unu_keys_time(indx_time_use);
total_unu_keys_time_use = length(unu_keys_time_use);
%
indx_time_storage = find(flags_time_storage==1); % flag=1 means data is available
unu_keys_time_storage = unu_keys_time(indx_time_storage);
total_unu_keys_time_storage = length(unu_keys_time_storage);

% stock

unu_keys_stock = nums_stock_avail(:,2);
%
flags_stock_use = nums_stock_avail(:,3);
flags_stock_storage = nums_stock_avail(:,4);
%
indx_stock_use = find(flags_stock_use==1); % flag=1 means data is available
unu_keys_stock_use = unu_keys_stock(indx_stock_use);
total_unu_keys_stock_use = length(unu_keys_stock_use);
%
indx_stock_storage = find(flags_stock_storage==1); % flag=1 means data is available
unu_keys_stock_storage = unu_keys_stock(indx_stock_storage);
total_unu_keys_stock_storage = length(unu_keys_stock_storage);

% fate: use the single date file for fates

all_entries_unu_keys_fate = nums_fate(:,2); % NOTE the column number
%
unu_keys_fate = unique(all_entries_unu_keys_fate);
total_unu_keys_fate = length(unu_keys_fate);

%%% subsets of unus belonging to multiple datasets

unu_keys_time_use_and_storage_simult = intersect(unu_keys_time_use, unu_keys_time_storage);
total_unu_keys_time_use_and_storage_simult = length(unu_keys_time_use_and_storage_simult);

unu_keys_stock_use_and_storage_simult = intersect(unu_keys_stock_use, unu_keys_stock_storage); 
total_unu_keys_stock_use_and_storage_simult = length(unu_keys_stock_use_and_storage_simult);

unu_keys_all_variables_simult = intersect(unu_keys_fate, intersect(unu_keys_time_use_and_storage_simult, unu_keys_stock_use_and_storage_simult));
total_unu_keys_all_variables_simult = length(unu_keys_all_variables_simult);

%%% extended pool of unus featuring at least in one dataset

unu_keys_entire_pool = union(unu_keys_fate, union(union(unu_keys_time_use, unu_keys_time_storage), union(unu_keys_stock_use, unu_keys_stock_storage)));
total_unu_keys_entire_pool = length(unu_keys_entire_pool);

unu_centuries_entire_pool = 100 * unique(floor(unu_keys_entire_pool/100));
total_unu_centuries_entire_pool = length(unu_centuries_entire_pool);

%%% extract cdfs for times, frequencies for stocks and frequencies for fates, wherever the data is available 

flags_all_variables = zeros(total_unu_keys_entire_pool, 8); % binary flags for data availability; columns: #, unu key, data available for time in use, time in storage, stock in use, stock in storage, fate, total variables available

% normalised pdfs for times in use and storage
all_entries_unu_keys_time = nums_time(:,2); % NOTE the column number
all_entries_timespan = nums_time(:,3); % NOTE the column number
all_entries_cdf_time_use = nums_time(:,5); % NOTE the column number
all_entries_cdf_time_storage = nums_time(:,7); % NOTE the column number

% normalised freqs for units in use and storage
all_entries_unu_keys_stock = nums_stock(:,2); % NOTE the column number
all_entries_stockspan = nums_stock(:,3); % NOTE the column number
all_entries_freq_stock_use = nums_stock(:,4); % NOTE the column number
all_entries_freq_stock_storage = nums_stock(:,6); % NOTE the column number

% normalised freqs for fates
all_entries_unu_keys_fate = nums_fate(:,2); % NOTE the column number
all_entries_cats_primary_fates = raw_fate(2:end,3); % NOTE the column and the row numbers! (skip first row when using raw cell array data)
all_entries_freq_fates = nums_fate(:,4); % NOTE the column number

cats_primary_fate = unique(all_entries_cats_primary_fates);
total_cats_primary_fate = numel(cats_primary_fate); % assuming the categorie are the same for each unu!


%%% extract unus and the relevant time series for POM and WG from WOT

pom_timeseries_years = nums_wot_pom_units(1,4:end); % NOTE the row and column numbers
%
all_entries_unu_keys_pom_units = nums_wot_pom_units(3:end,2); % NOTE the row and column numbers
all_entries_timeseries_pom_units = nums_wot_pom_units(3:end,4:end); % NOTE the row and column numbers
%
all_entries_unu_keys_pom_weightpu = nums_wot_pom_weightpu(3:end,2); % NOTE the row and column numbers
all_entries_timeseries_pom_weightpu = nums_wot_pom_weightpu(3:end,4:end); % NOTE the row and column numbers; UNITS: kg!!!
%
all_entries_unu_keys_wg_tonnes = nums_wot_wg_tonnes(3:end,2); % NOTE the row and column numbers
all_entries_timeseries_wg_tonnes = nums_wot_wg_tonnes(3:end,4:end); % NOTE the row and column numbers

% unique unu keys for pom and wg (just in case there are repetitions
unu_keys_pom_units = unique(all_entries_unu_keys_pom_units);
unu_keys_pom_weightpu = unique(all_entries_unu_keys_pom_weightpu);
unu_keys_wg_tonnes = unique(all_entries_unu_keys_wg_tonnes);

if sum(abs(unu_keys_pom_weightpu - unu_keys_pom_units)) > 0  || sum(abs(unu_keys_wg_tonnes - unu_keys_pom_units)) > 0
    error('Error: WOT1.2 UNU entries for POM units, POM weight per unit and WG tonnes do not match...')   
end
    
%
unu_labels_pom_units = raw_wot_pom_units(3:end,3); % NOTE the row and column numbers


%%% extract ONS number of households time series

households_timeseries_years = nums_ons(end-3,:); % data timepoints for uk households
uk_households = nums_ons(end,:); % time_series data for the households

% extrapolate the number of households beyond the last available year (2017
% in the currently available ONS dataset; could be updated later on)
uk_households_extrap = interp1(households_timeseries_years,uk_households,new_timeseries_timerange,'linear','extrap');

current_uk_households = uk_households_extrap(new_timeseries_timerange==current_year);


%%% execute subroutine for extracting unu-level data from raw entries

extract_unu_data


%%% execute subroutine for regridding time distributions (-> annual step)

regrid_time_distrib


%%% execute subroutine for extrapolating time in use, time in storage and
%%% fate to unus that do not have the relevant data; the extrapolations are
%%% scored according to the accuracy (how far the available data is from
%%% the current unu)

time_fate_unu_extrap


%%% finally, execute subroutine for estimating secondary fates (by removing storage and re-distributing its frequency among the other non-trivial fates) 

secondary_fate


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% now we can start working out the flows

allunudata(total_unu_keys_entire_pool).results.current_house_mean_stock_use_units = 0;
allunudata(total_unu_keys_entire_pool).results.current_uk_mean_stock_use_tonnes = 0;

allunudata(total_unu_keys_entire_pool).results.current_house_mean_stock_storage_units = 0;
allunudata(total_unu_keys_entire_pool).results.current_uk_mean_stock_storage_tonnes = 0;

allunudata(total_unu_keys_entire_pool).results.ts_pom_units(new_timeseries_timerange) = zeros;

allunudata(total_unu_keys_entire_pool).results.ts_pom_tonnes(new_timeseries_timerange) = zeros;
allunudata(total_unu_keys_entire_pool).results.ts_discarded_after_use_tonnes(new_timeseries_timerange) = zeros;
allunudata(total_unu_keys_entire_pool).results.ts_discarded_after_storage_tonnes(new_timeseries_timerange) = zeros;


for k = total_unu_keys_entire_pool:-1:1 % backwards to allocate structure arrays  
    
    % analysing the metadata, one variable at a time
    
    % stock in use => work out the mean
    
    if flags_all_variables(k,5) == 1 % stock in use 
        
        sss = allunudata(k).distributions.stockspan;
        fff = allunudata(k).distributions.freq_stock_use;
        
        fff(isnan(fff)) = 0; % failsafe for pdf
                      
        % use matrix multiplication to work out the sum representing the
        % mean
                
        allunudata(k).results.current_house_mean_stock_use_units = transpose(sss) * fff;
        
        % convert mean stock per household into total uk weight
        current_weightpu = (10^-3) * allunudata(k).timeseries.pom_weightpu(pom_timeseries_years==current_year); % kg->tonnes; current year only
        
        allunudata(k).results.current_uk_mean_stock_use_tonnes = current_uk_households * current_weightpu *...
            allunudata(k).results.current_house_mean_stock_use_units;
                        
    end
    
    % stock in storage => work out the mean
    
    if flags_all_variables(k,6) == 1 % stock in storage 
       
        sss = allunudata(k).distributions.stockspan;
        fff = allunudata(k).distributions.freq_stock_storage;
        
        fff(isnan(fff)) = 0; % failsafe for pdf
                      
        % use matrix multiplication to work out the sum representing the
        % mean per household
        allunudata(k).results.current_house_mean_stock_storage_units = transpose(sss) * fff;
        
        % convert mean stock per household into total uk weight
        current_weightpu = (10^-3) * allunudata(k).timeseries.pom_weightpu(pom_timeseries_years==current_year); % kg->tonnes; current year only
        
        allunudata(k).results.current_uk_mean_stock_storage_tonnes = current_uk_households * current_weightpu *...
            allunudata(k).results.current_house_mean_stock_storage_units;
    
    end   
   
    % time in use => work out time series for discarded units
    
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
        
    pdf = allunudata(k).distributions.new_pdf_time_use;
           
    %%% proceed with working out the discarded flow for all unus
    
    timespan_pdf = new_pdf_timespan;
    %
    flow_units = allunudata(k).timeseries.pom_units;
    %
    pom_weightpu = (10^-3) * allunudata(k).timeseries.pom_weightpu; % kg into tonnes
    %
    flow_tonnes = flow_units .* pom_weightpu; % element by element multiplication (not matrix multiplication)
    %
    timespan_flow = pom_timeseries_years;
     
    timespan_out = new_timeseries_timerange;
        
    % work out the historic probability integral, both for the flows in units
    % and tonnes
    [integral_units, flow_units_extrap] = probability_integral(flow_units,timespan_flow,pdf,timespan_pdf,timespan_out);
    [integral_tonnes, flow_tonnes_extrap] = probability_integral(flow_tonnes,timespan_flow,pdf,timespan_pdf,timespan_out);

    allunudata(k).results.ts_pom_units = flow_units_extrap;
    allunudata(k).results.ts_pom_tonnes = flow_tonnes_extrap;
    
    allunudata(k).results.ts_discarded_after_use_tonnes = integral_tonnes;
    
    % now apply the fates data to the discarded flow
    
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
       
    freq = allunudata(k).distributions.freq_fate;
               
    % loop over fate splits to estimate primary fate flows 
    
    ts_primary_fate = zeros(total_cats_primary_fate,total_new_timeseries_timesteps);
    
    for j = 1:total_cats_primary_fate
    
        for i = 1:total_new_timeseries_timesteps
                    
            ts_primary_fate(j,i) = freq(j) .* allunudata(k).results.ts_discarded_after_use_tonnes(i);
    
        end
        
%         if unu_keys_entire_pool(k) == 303 && j==6
%            
%             fprintf('discarded after use %f, into hoarding %f\n', allunudata(k).results.ts_discarded_after_use_tonnes(end-2),ts_primary_fate(j,end-2))
%                 
%             disp('debugging...')
%             
%         end
        
        
    end
        
    allunudata(k).results.ts_primary_fate_tonnes = ts_primary_fate;
    
    % now use the primary flow into storage calculated above to estimate the discarded flow after storage 
       
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
    
    pdf = allunudata(k).distributions.new_pdf_time_storage;
            
    %%% proceed with working out the discarded flow for all unus
    
    timespan_pdf = new_pdf_timespan;
    
    cats_primary_fate_current_unu = allunudata(k).distributions.cats_primary_fate;
    [~,j_storage] = ismember('Storage',cats_primary_fate_current_unu);
        
    flow_tonnes = ts_primary_fate(j_storage,:);
    
    timespan_flow = new_timeseries_timerange; % NOTE: the discarded flow already has the correct time range
     
    timespan_out = new_timeseries_timerange;
        
    % work out the historic probability integral
    [integral_tonnes, ~] = probability_integral(flow_tonnes,timespan_flow,pdf,timespan_pdf,timespan_out);
     
    allunudata(k).results.ts_discarded_after_storage_tonnes = integral_tonnes;  
      
    %%% finally, apply the approximate secondary fate frequencies to the flow discarded after storage to obtain secondary fate flows 
    
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
               
    freq = allunudata(k).distributions.freq_fate_secondary;
         
    % loop over fate splits to estimate primary fate flows 
    
    ts_secondary_fate = zeros(total_cats_primary_fate,total_new_timeseries_timesteps); % NOTE: we kept the number of categories the same and assigned zero freq to storage in the secondary fates (re-distributing the primary storage freq equally among all other categories)
    
    for j = 1:total_cats_primary_fate
    
        for i = 1:total_new_timeseries_timesteps
                                
            ts_secondary_fate(j,i) = freq(j) .* allunudata(k).results.ts_discarded_after_storage_tonnes(i); % NOTE the variable name: discarded after storage!!!!
    
        end
        
    end
        
    allunudata(k).results.ts_secondary_fate_tonnes = ts_secondary_fate;
    
    % combined fate flows
    
    ts_combined_fate = ts_primary_fate + ts_secondary_fate;
    
    allunudata(k).results.ts_combined_fate_tonnes = ts_combined_fate;
    
    
%     figure
%     hold on
%     plot(new_yearrange,ts_combined_fate(1,:))
%     plot(new_yearrange,ts_combined_fate(2,:))
%     plot(new_yearrange,ts_combined_fate(3,:))
%     plot(new_yearrange,ts_combined_fate(4,:))
%     plot(new_yearrange,ts_combined_fate(5,:))
%     plot(new_yearrange,ts_combined_fate(6,:))
%     plot(new_yearrange,ts_combined_fate(7,:))
%     plot(new_yearrange,ts_combined_fate(8,:))
%     hold off
%     title(['discarded after storage, UNU=' num2str(flags_all_variables(k,2))])

end


%%% write the results into aux arrays into order to save in excel

aux_array_results_pom = NaN * ones(total_unu_keys_entire_pool,3);
aux_array_results_stock = NaN * ones(total_unu_keys_entire_pool,4);
aux_array_results_weee_totals = NaN * ones(total_unu_keys_entire_pool,3);
aux_array_results_weee_fates = NaN * ones(total_unu_keys_entire_pool,total_cats_primary_fate);

aux_unu_labels = cell(total_unu_keys_entire_pool,1);

ind_current_year_original = find(pom_timeseries_years==current_year);
ind_current_year_extended = find(new_timeseries_timerange==current_year);

for k = 1:total_unu_keys_entire_pool
    
    % this is available for all UNUs
    aux_array_results_pom(k,1) = allunudata(k).timeseries.pom_weightpu(ind_current_year_original); %#ok<*FNDSB> % NOTE: this time series uses a different timescale (original from WOT rather than extended for the subsequent analysis)
    aux_array_results_pom(k,2) = allunudata(k).results.ts_pom_units(ind_current_year_extended);
    aux_array_results_pom(k,3) = allunudata(k).results.ts_pom_tonnes(ind_current_year_extended);
    
    %
    
    % NOTE: we do not extrapolate the stock data to other UNUs at this
    % stage -- further analysis is required
    if flags_all_variables(k,5) == 1 % stock in use
    
        aux_array_results_stock(k,1) = allunudata(k).results.current_house_mean_stock_use_units;
        aux_array_results_stock(k,2) = allunudata(k).results.current_uk_mean_stock_use_tonnes;

    end
    
    if flags_all_variables(k,6) == 1 % stock in storage
        
        aux_array_results_stock(k,3) = allunudata(k).results.current_house_mean_stock_storage_units; 
        aux_array_results_stock(k,4) = allunudata(k).results.current_uk_mean_stock_storage_tonnes; 

    end
    
    %
    
    aux_array_results_weee_totals(k,1) = allunudata(k).results.ts_discarded_after_use_tonnes(ind_current_year_extended); 
    aux_array_results_weee_totals(k,2) = allunudata(k).results.ts_discarded_after_storage_tonnes(ind_current_year_extended); 

    aux_array_results_weee_fates(k,1:total_cats_primary_fate) = allunudata(k).results.ts_combined_fate_tonnes(1:total_cats_primary_fate,ind_current_year_extended);

    aux_array_results_weee_totals(k,3) = sum(aux_array_results_weee_fates(k,:));
    
    ind_unu = find(unu_keys_pom_units == unu_keys_entire_pool(k)); % locate current unu key in the extended set available for pom
    
    aux_unu_labels{k} = unu_labels_pom_units{ind_unu};
    
end

%%% add unu number and labels as the first two colums and write into excel

% flows and flags

aux_array_flags = flags_all_variables(:,3:8); % flags indicating data availability for time in use, time in storage, stock in use, stock in storage, fate, total variables available

aux_cell = [num2cell(unu_keys_entire_pool) aux_unu_labels... 
    num2cell(aux_array_results_pom)...
    num2cell(aux_array_results_weee_totals) num2cell(aux_array_results_weee_fates)... 
    num2cell(aux_array_flags)...
    num2cell(scores_time_and_fate_unu_extrap)];

columns_labels_unu = {'UNU', 'Description'};

columns_labels_results_pom = {'Average weight per unit (kg)',  'POM (units/yr)', 'POM (tonnes/yr)'};
columns_labels_results_weee_totals = {
    'Flow discarded after use (tonnes/yr)',...
    'Flow discarded after hoarding (tonnes/yr)',...
    'Total WEEE (tonnes/yr)'};
columns_labels_results_weee_fates = {
    'Flow "Donation or re-use" (tonnes/yr)'...
    'Flow "General bin" (tonnes/yr)'...
    'Flow "Other" (tonnes/yr)'...
    'Flow "Recycling" (tonnes/yr)'...
    'Flow "Sold" (tonnes/yr)'...
    'Flow "Hoarding" (tonnes/yr)'...
    'Flow "Take-back scheme" (tonnes/yr)'...
    'Flow "Unknown" (tonnes/yr)'};
columns_labels_flags = {'avail time in use (0/1)', 'avail time hoarded (0/1)',...
    'avail stock in use (0/1)', 'avail stock hoarded (0/1)',...
    'avail fates (0/1)', 'total avail (0/5)'};
columns_labels_scores_unu_extrap = {...
    'score extrap time use (of 5)'...
    'score extrap time storage (of 5)'...
    'score extrap fate (of 5)'};

columns_labels = [columns_labels_unu...
    columns_labels_results_pom...
    columns_labels_results_weee_totals...
    columns_labels_results_weee_fates...
    columns_labels_flags...
    columns_labels_scores_unu_extrap];

output_matrix = [columns_labels; aux_cell];

sheet = ['flows_' num2str(current_year)];

xlswrite(filename_output_prelim,output_matrix,sheet)

% stocks

aux_cell = [num2cell(unu_keys_entire_pool) aux_unu_labels... 
    num2cell(aux_array_results_stock)];

columns_labels_results_stock = {'Stock in use per household (units)', 'Stock in use UK (tonnes)'...
    'Stock hoarded per household (units)', 'Stock hoarded UK (tonnes)'};

columns_labels = [columns_labels_unu columns_labels_results_stock];

output_matrix = [columns_labels; aux_cell];

sheet = ['stocks_' num2str(current_year)];

xlswrite(filename_output_prelim,output_matrix,sheet)


% return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Advanced analysis: combining times and stocks to reconstruct their 
%%% evolution based on the conservation of units 
%%%

timespan_out = new_timeseries_timerange;

timesteps_pdf_scale_factor = 1980:1:2015; % cover all past decades but last one 
total_timeperiods_pdf_scale_factor = numel(timesteps_pdf_scale_factor) - 1; % periods WITHIN the specified steps => one less

pdf_scale_factor_base = ones(1,total_timeperiods_pdf_scale_factor); 

% allunudata(total_unu_keys_entire_pool).results_advanced.current_house_mean_stock_use_units = 0;
% allunudata(total_unu_keys_entire_pool).results_advanced.current_uk_mean_stock_use_tonnes = 0;
% 
% allunudata(total_unu_keys_entire_pool).results_advanced.current_house_mean_stock_storage_units = 0;
% allunudata(total_unu_keys_entire_pool).results_advanced.current_uk_mean_stock_storage_tonnes = 0;
% 
% allunudata(total_unu_keys_entire_pool).results_advanced.ts_discarded_after_use_tonnes(new_timeseries_timerange) = zeros;
% allunudata(total_unu_keys_entire_pool).results_advanced.ts_discarded_after_storage_tonnes(new_timeseries_timerange) = zeros;
% 
% allunudata(total_unu_keys_entire_pool).results_advanced.ts_uk_mean_stock_use_units(new_timeseries_timerange) = zeros;
% allunudata(total_unu_keys_entire_pool).results_advanced.ts_uk_mean_stock_storage_units(new_timeseries_timerange) = zeros;

timespan_pdf = new_pdf_timespan;
timespan_flow = pom_timeseries_years;

year_target_stock = current_year;

aux_stock_match_labels = cell(total_unu_keys_entire_pool,1);

for k=1:total_unu_keys_entire_pool
   
    %%% for all unus: extract pdf for time in use and pom in units
    
    pdf = allunudata(k).distributions.new_pdf_time_use;
        
    flow_units = allunudata(k).timeseries.pom_units;
            
    %%% perform advanced analysis for the unus with stock in use data
    
    if flags_all_variables(k,5) == 1 % stock in use
        
        %

        target_stock = current_uk_households * allunudata(k).results.current_house_mean_stock_use_units;
        fprintf('target_stock = %e',target_stock)

        %

        %%% now use nonlinear programming optimisation to minimise the misfit between
        %%% between the reconstructed and observed stock

        % x <=> pdf_scale_factor

        fun = @(x) stock_objective_function(flow_units,timespan_flow,...
        pdf,timespan_pdf,x,timesteps_pdf_scale_factor,timespan_out,...
        target_stock,year_target_stock);

        % use function handle to define the function of x to be minimised

        % define initial condition for the optimisation: unit scaling factors
        % for all the scaling periods (i.e. use present-day distribution)

        x0 = pdf_scale_factor_base;

        % set the upper and lower boundaries for x
        % NOTE: these are rows
        lb = 0.5 * x0; % no less than *** of the present-day pdf timespan 
        ub = 2.0 * x0; % no more than *** times the present-day pdf timespan 


        fprintf('\n')
        text = ['running fmincon for UNU=' num2str(unu_keys_entire_pool(k)) '...'];
        fprintf(text) 

        [x_opt, fmin] = fmincon(fun,x0,[],[],[],[],lb,ub,[],options); % NOTE the excluded fields (must be replaced with [] for the solver to work)

        fprintf('done!\n')

    %     disp(x_opt)

        fprintf('UNU %i: min misfit relative to target stock = %e\n',unu_keys_entire_pool(k),fmin/target_stock)

        pdf_scale_factor_opt = x_opt; % optimal scaling for the pdf time in use
        % pdf_scale_factor_base is set outside the loop
        
        [discarded_base, ~] = probability_integral_advanced(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_base,timesteps_pdf_scale_factor,timespan_out);

        [discarded_opt, flow_extrap] = probability_integral_advanced(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_opt,timesteps_pdf_scale_factor,timespan_out);

        [stock_base] = stock_ode(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_base,timesteps_pdf_scale_factor,timespan_out);

        [stock_opt] = stock_ode(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_opt,timesteps_pdf_scale_factor,timespan_out);

    
        %%% plotting the reconstructed discarded units and stocks

        fig = figure('Position', get(0, 'Screensize'));
        set(gca,'FontSize',14)
        set(gcf,'color','w')
        set(fig, 'Visible', 'off');

        xmin = 1980; 
        xmax = 2020;
        ymin = 0;

        text = ['UNU ' num2str(unu_keys_entire_pool(k))]; 

        subplot(1,2,1)
        hold on
        l1=plot(new_timeseries_timerange,flow_extrap,'b');
        l2=plot(new_timeseries_timerange,discarded_base,'k', 'LineWidth',1.5);
        l3=plot(new_timeseries_timerange,discarded_opt,'r', 'LineWidth',1.5);
        hold off

        legend([l1 l2 l3],{'POM WOT1.2','discarded basic est.','discarded advanced est.'},'Location','northwest')

        set(gca,'xlim',[xmin xmax])
        set(gca,'ylim',[ymin inf])

        xlabel('year')
        ylabel('units')

        title(['Flows POM & discarded after use, ' text]) 

        subplot(1,2,2)
        hold on
        l1=plot(year_target_stock,target_stock,'m*', 'LineWidth',1.5);
        l2=plot(new_timeseries_timerange,stock_base,'k', 'LineWidth',1.5);
        l3=plot(new_timeseries_timerange,stock_opt,'r', 'LineWidth',1.5);
        hold off

        legend([l1 l2 l3],{'stock data','stock basic est.','stock advanced est.'},'Location','northwest')

        set(gca,'xlim',[xmin xmax])
        set(gca,'ylim',[ymin inf])

        xlabel('year')
        ylabel('units')

        title(['Stock in use, ' text]) 

        filename = ['plots/' text '.png']; 
        saveas(fig, filename, 'png')

        %%% based on the plots, we assign the following labels:
        
        if ismember(unu_keys_entire_pool(k),[301 901])
            
            aux_stock_match_labels{k} = 'mismatch';
            
        else
            
            aux_stock_match_labels{k} = 'match';
            
        end
        
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% put forward composite product for discarded flows for those unus
        %%% that achieved good fit; for others, discard te stock data and use
        %%% the estimated stock instead
               
        % advanced analysis
        allunudata(k).results_advanced.ts_uk_mean_stock_use_units(:) = stock_opt(:);
        allunudata(k).results_advanced.ts_discarded_after_use_units(:) = discarded_opt(:);
        
        % basic analysis
        allunudata(k).results.ts_uk_mean_stock_use_units(:) = stock_base(:);
        allunudata(k).results.ts_discarded_after_use_units(:) = discarded_base(:);

        
        % composite results: average of the two above
        allunudata(k).results_composite.ts_uk_mean_stock_use_units = ...
            0.5*(allunudata(k).results_advanced.ts_uk_mean_stock_use_units + allunudata(k).results.ts_uk_mean_stock_use_units);
        %
        allunudata(k).results_composite.ts_discarded_after_use_units =...
            0.5*(allunudata(k).results_advanced.ts_discarded_after_use_units + allunudata(k).results.ts_discarded_after_use_units);
        
        
    else    
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% for all other unus, execute stock calculation routine with unit
        %%% scaling and put this forward as "composite"

        [discarded_base, ~] = probability_integral_advanced(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_base,timesteps_pdf_scale_factor,timespan_out);

        [stock_base] = stock_ode(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_base,timesteps_pdf_scale_factor,timespan_out);

    
        %%% recording the basic analysis results for thes unus into the composite structure
        
        allunudata(k).results_composite.ts_uk_mean_stock_use_units(:) = stock_base(:);
        allunudata(k).results_composite.ts_discarded_after_use_units(:) = discarded_base(:);
    
        %%%
        
    end
    
    %%% now we have the data for uk-wide stock in use and discarded after
    %%% use as number of units; convert them to number of units per
    %%% household and total uk tonnage
   
    %%% results in the current year (could add full time series later)
    
    allunudata(k).results_composite.current_house_mean_stock_use_units =...
        allunudata(k).results_composite.ts_uk_mean_stock_use_units(new_timeseries_timerange==current_year) ./ current_uk_households;        
       
    current_weightpu = (10^-3) * allunudata(k).timeseries.pom_weightpu(pom_timeseries_years==current_year); % kg -> tonnes
     
    allunudata(k).results_composite.current_uk_mean_stock_use_tonnes =... 
        allunudata(k).results_composite.ts_uk_mean_stock_use_units(new_timeseries_timerange==current_year) * current_weightpu;
    
end % end of the loop over unus



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% now we can proceed with the remaining steps of the analysis: work out
%%% primary storage flow, storage stock and tonnage, and finally other
%%% fates, all based on the composite primary discarded after use flow 
%%% calculated above; use the template from the preliminary results, 
%%% including the spreadsheet summary

%%% NOTE, unlike in the preliminary runs, work everything out in units and
%%% only convert to tonnes in the end (for the current year)

for k = 1:total_unu_keys_entire_pool  
              
    aux_stock_match_labels{k} = 'NA'; % no pairs of time in use and stock in use data for the rest of the UNUs

    % now apply the fates data to the discarded flow after use
    
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
       
    freq = allunudata(k).distributions.freq_fate;
               
    % loop over fate splits to estimate primary fate flows 
    
    ts_primary_fate_units = zeros(total_cats_primary_fate,total_new_timeseries_timesteps);
    
    for j = 1:total_cats_primary_fate
    
        for i = 1:total_new_timeseries_timesteps
                    
            % NOTE: use the composite primary discared flow calculated above!!!
            ts_primary_fate_units(j,i) = freq(j) .* allunudata(k).results_composite.ts_discarded_after_use_units(i);
    
        end
        
%         if unu_keys_entire_pool(k) == 303 && j==6
%            
%             fprintf('discarded after use %f, into hoarding %f\n', allunudata(k).results_composite.ts_discarded_after_use_tonnes(end-2),ts_primary_fate(j,end-2))
%                 
%             disp('debugging...')
%             
%         end
        
        
    end
        
    allunudata(k).results_composite.ts_primary_fate_units = ts_primary_fate_units;
    
    % now use the primary flow into storage calculated above to estimate the discarded flow after storage 
       
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
      
    cats_primary_fate_current_unu = allunudata(k).distributions.cats_primary_fate;
    [~,j_storage] = ismember('Storage',cats_primary_fate_current_unu);
        
    flow_units = ts_primary_fate_units(j_storage,:);
    timespan_flow = new_timeseries_timerange; % NOTE: the discarded flow already has the correct time range
       
    %
    
    pdf = allunudata(k).distributions.new_pdf_time_storage;      
    timespan_pdf = new_pdf_timespan;
         
    timespan_out = new_timeseries_timerange;
    
    [discarded_base, ~] = probability_integral_advanced(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_base,timesteps_pdf_scale_factor,timespan_out);

    [stock_base] = stock_ode(flow_units,timespan_flow,...
        pdf,timespan_pdf,pdf_scale_factor_base,timesteps_pdf_scale_factor,timespan_out);
    
    %%%
    
    allunudata(k).results_composite.ts_discarded_after_storage_units(:) = discarded_base(:);  
    allunudata(k).results_composite.ts_uk_mean_stock_storage_units(:) = stock_base(:);
    
    %%% now we have the data for uk-wide stock in use and discarded after
    %%% use as number of units; convert them to number of units per
    %%% household and total uk tonnage
   
    %%% results in the current year (could add full time series later)
    
    allunudata(k).results_composite.current_house_mean_stock_storage_units =...
        allunudata(k).results_composite.ts_uk_mean_stock_storage_units(new_timeseries_timerange==current_year) ./ current_uk_households;        
       
    current_weightpu = (10^-3) * allunudata(k).timeseries.pom_weightpu(pom_timeseries_years==current_year); % kg -> tonnes
     
    allunudata(k).results_composite.current_uk_mean_stock_storage_tonnes =... 
        allunudata(k).results_composite.ts_uk_mean_stock_storage_units(new_timeseries_timerange==current_year) * current_weightpu;
    
    
    
    %%% finally, apply the approximate secondary fate frequencies to the flow discarded after storage to obtain secondary fate flows 
    
    % NOTE: the required pdf for the given unu is either available from the data or has been extrapolated
    % from other unus
               
    freq = allunudata(k).distributions.freq_fate_secondary;
         
    % loop over fate splits to estimate primary fate flows 
    
    ts_secondary_fate_units = zeros(total_cats_primary_fate,total_new_timeseries_timesteps); % NOTE: we kept the number of categories the same and assigned zero freq to storage in the secondary fates (re-distributing the primary storage freq equally among all other categories)
    
    for j = 1:total_cats_primary_fate
    
        for i = 1:total_new_timeseries_timesteps
                                
            ts_secondary_fate_units(j,i) = freq(j) .* allunudata(k).results_composite.ts_discarded_after_storage_units(i); % NOTE the variable name: discarded after storage!!!!
    
        end
        
    end
        
    allunudata(k).results_composite.ts_secondary_fate_units = ts_secondary_fate_units;
    
    % combined fate flows
    
    ts_combined_fate_units = ts_primary_fate_units + ts_secondary_fate_units;
    
    allunudata(k).results_composite.ts_combined_fate_units = ts_combined_fate_units;
    
    %%%%
    
    
    
    
%     figure
%     hold on
%     plot(new_yearrange,ts_combined_fate(1,:))
%     plot(new_yearrange,ts_combined_fate(2,:))
%     plot(new_yearrange,ts_combined_fate(3,:))
%     plot(new_yearrange,ts_combined_fate(4,:))
%     plot(new_yearrange,ts_combined_fate(5,:))
%     plot(new_yearrange,ts_combined_fate(6,:))
%     plot(new_yearrange,ts_combined_fate(7,:))
%     plot(new_yearrange,ts_combined_fate(8,:))
%     hold off
%     title(['discarded after storage, UNU=' num2str(flags_all_variables(k,2))])

end


%%% write the results into aux arrays into order to save in excel

aux_array_results_pom = NaN * ones(total_unu_keys_entire_pool,3);
aux_array_results_stock = NaN * ones(total_unu_keys_entire_pool,4);
% aux_array_results_stock_change = NaN * ones(total_unu_keys_entire_pool,2);
aux_array_results_weee_totals = NaN * ones(total_unu_keys_entire_pool,3);
aux_array_results_weee_fates = NaN * ones(total_unu_keys_entire_pool,total_cats_primary_fate);

aux_unu_labels = cell(total_unu_keys_entire_pool,1);

ind_current_year_original = find(pom_timeseries_years==current_year);
ind_current_year_extended = find(new_timeseries_timerange==current_year);

for k = 1:total_unu_keys_entire_pool
    
    % this is available for all UNUs
    aux_array_results_pom(k,1) = allunudata(k).timeseries.pom_weightpu(ind_current_year_original); % NOTE: this time series uses a different timescale (original from WOT rather than extended for the subsequent analysis)
    aux_array_results_pom(k,2) = allunudata(k).results.ts_pom_units(ind_current_year_extended);
    aux_array_results_pom(k,3) = allunudata(k).results.ts_pom_tonnes(ind_current_year_extended);
    
    % NOTE: results for current stock in use and stock in storage are now
    % recorded for all unus
    
    aux_array_results_stock(k,1) = allunudata(k).results_composite.current_house_mean_stock_use_units;
    aux_array_results_stock(k,2) = allunudata(k).results_composite.current_uk_mean_stock_use_tonnes;

    aux_array_results_stock(k,3) = allunudata(k).results_composite.current_house_mean_stock_storage_units; 
    aux_array_results_stock(k,4) = allunudata(k).results_composite.current_uk_mean_stock_storage_tonnes; 

   
    % NOTE: we now need to apply units -> weights conversion (unlike in the
    % preliminary analysis)
    
    current_weightpu = (10^-3) * allunudata(k).timeseries.pom_weightpu(pom_timeseries_years==current_year); % kg -> tonnes
      
    aux_array_results_weee_totals(k,1) = current_weightpu * allunudata(k).results_composite.ts_discarded_after_use_units(ind_current_year_extended); 
    aux_array_results_weee_totals(k,2) = current_weightpu * allunudata(k).results_composite.ts_discarded_after_storage_units(ind_current_year_extended); 

    aux_array_results_weee_fates(k,1:total_cats_primary_fate) = current_weightpu * allunudata(k).results_composite.ts_combined_fate_units(1:total_cats_primary_fate,ind_current_year_extended);

    aux_array_results_weee_totals(k,3) = sum(aux_array_results_weee_fates(k,:)); % total discarded
    
    ind_unu = find(unu_keys_pom_units == unu_keys_entire_pool(k)); % locate current unu key in the extended set available for pom
    
    aux_unu_labels{k} = unu_labels_pom_units{ind_unu};
    
end

%%% add unu number and labels as the first two colums and write into excel

% flags

aux_array_flags = flags_all_variables(:,3:8); % flags indicating data availability for time in use, time in storage, stock in use, stock in storage, fate, total variables available

% new output array -- re-arranged order (triplet "POM, Delta in use,
% Discarded after use"; triplet "Flow into hoarding, Delta hoarded, Discarded after hoarding"; Total discarded; 7 Fates (excluding hoarding))

% also, change from ton to kton in the mass flow entries; NaNs stay as
% they are, giving empty columns after executing num2cell

aux_array_rearranged = NaN * ones(total_unu_keys_entire_pool,20); % NaNs will be converted into empty columns after executing num2cell

% cols 1, 2: POM weight per unit, POM units

aux_array_rearranged(:,1:2) = aux_array_results_pom(:,1:2); % no scaling for kg/unit and number of units!

% col 3: empty
% cols 4, 5, 6: POM, Delta in use, Discarded after use

aux_array_rearranged(:,4) = (1e-3) * aux_array_results_pom(:,3); % POM tonnes
% NOTE the columns order!!!
aux_array_rearranged(:,6) = (1e-3) * aux_array_results_weee_totals(:,1); % Discarded after use
%
aux_array_rearranged(:,5) = aux_array_rearranged(:,4) - aux_array_rearranged(:,6); % Delta in use

% col 7: empty
% cols 8, 9, 10: Flow into hoarding, Delta hoarded, Discarded after hoarding

aux_array_rearranged(:,8) = (1e-3) * aux_array_results_weee_fates(:,6); % Flow into hoarding (fate #6)
% NOTE the columns order!!!
aux_array_rearranged(:,10) = (1e-3) * aux_array_results_weee_totals(:,2); % Discarded after hoarding
%
aux_array_rearranged(:,9) = aux_array_rearranged(:,8) - aux_array_rearranged(:,10); % Delta hoarded

% col 11: empty
% col 12: Total discarded 
aux_array_rearranged(:,12) = (1e-3) * aux_array_results_weee_totals(:,3); % Total discarded

% col 13: empty
% cols 14-20: 7 Fates (excluding hoarding, which is fate #6) => use fates
% 1-5 follwed by 7-8

aux_array_rearranged(:,14:18) = (1e-3) * aux_array_results_weee_fates(:,1:5); % first 5 fates 
aux_array_rearranged(:,19:20) = (1e-3) * aux_array_results_weee_fates(:,7:8); % remaining 5 fates 

% empty column before flags and scores

aux_array_empty_column = NaN * ones(total_unu_keys_entire_pool,1);

% combine all entries together in a single cell; 

aux_cell = [num2cell(unu_keys_entire_pool) aux_unu_labels... 
    num2cell(aux_array_rearranged)...
    num2cell(aux_array_empty_column)...
    num2cell(aux_array_flags)...
    num2cell(scores_time_and_fate_unu_extrap)];

columns_labels_unu = {'UNU', 'Description'};
%
columns_labels_results_rearranged = {...
    'Average weight per unit (kg)', 'POM (units/yr)'...
    'Flows in use ->'...
    'POM (kton/yr)', 'Rate of change for stock in use (kton/yr)', 'Flow discarded after use (kton/yr)'...
    'Flows in hoarding ->'...
    'Flow into "Hoarding" (kton/yr)' 'Rate of change for stock in hoarding (kton/yr)' 'Flow discarded after hoarding (kton/yr)'...
    'Totals ->'...
    'Total WEEE (kton/yr)'...
    'EEE & WEEE Fates ->'...
    'Flow "Donation or re-use" (kton/yr)'...
    'Flow "General bin" (kton/yr)'...
    'Flow "Other" (kton/yr)'...
    'Flow "Recycling" (kton/yr)'...
    'Flow "Sold" (kton/yr)'...
    'Flow "Take-back scheme" (kton/yr)'...
    'Flow "Unknown" (kton/yr)'};   
%
columns_labels_empty = {'Data availability and confidence ->'};
%     
columns_labels_flags = {'avail time in use (0/1)', 'avail time hoarded (0/1)',...
    'avail stock in use (0/1)', 'avail stock hoarded (0/1)',...
    'avail fates (0/1)', 'total avail (0/5)'};
% 
columns_labels_scores_unu_extrap = {...
    'score extrap time use (of 5)'...
    'score extrap time storage (of 5)'...
    'score extrap fate (of 5)'};

columns_labels = [columns_labels_unu...
    columns_labels_results_rearranged...
    columns_labels_empty...
    columns_labels_flags...
    columns_labels_scores_unu_extrap];

output_matrix = [columns_labels; aux_cell];

sheet = ['flows_' num2str(current_year)];

xlswrite(filename_output_advanced,output_matrix,sheet)

% stocks

aux_cell = [num2cell(unu_keys_entire_pool) aux_unu_labels... 
    num2cell(aux_array_results_stock)...
    aux_stock_match_labels];

columns_labels_results_stock = {'Stock in use per household (units)', 'Stock in use UK (tonnes)'...
    'Stock hoarded per household (units)', 'Stock hoarded UK (tonnes)'};

columns_labels_stock_match = {'stock in use match (reconstr. vs metadata)'};

columns_labels = [columns_labels_unu...
    columns_labels_results_stock...
    columns_labels_stock_match];

output_matrix = [columns_labels; aux_cell];

sheet = ['stocks_' num2str(current_year)];

xlswrite(filename_output_advanced,output_matrix,sheet)