Update code (fix global vaccination error)

matlab scripts/global_vacc_pred.m

@@ -25,9 +25,9 @@
 
 %% Prepare vaccine data
 xx = load('vacc_data.mat');
-nidx = isnan(latest_vac); latest_vac(nidx) = popu(nidx).*(mean(latest_vac(~nidx))/mean(popu(~nidx)));
+%nidx = isnan(latest_vac); latest_vac(nidx) = popu(nidx).*(mean(latest_vac(~nidx))/mean(popu(~nidx)));
 vacc_lag_rate = ones(length(popu), 14);
-vacc_ad_dist(:) = 1; % This round specifies administration, rather than distribution. No need to model distribution rates
+vacc_ad_dist(:) = 1; 
 vacc_all = xx.g_vacc(:, 1:end); vacc_full = xx.g_vacc_full(:, 1:end); 
 vacc_all = fillmissing(vacc_all, 'previous', 2); vacc_all(isnan(vacc_all)) = 0;
 vacc_full = fillmissing(vacc_full, 'previous', 2); vacc_full(isnan(vacc_full)) = 0;

matlab scripts/quant_gen_base.m

@@ -3,7 +3,7 @@
 
 num_ahead = 4;
 quant_deaths = [0.01, 0.025, (0.05:0.05:0.95), 0.975, 0.99];
-quant_cases = [-0.975, 0.100, 0.250, 0.500, 0.750, 0.900, 0.975];
+quant_cases = [0.025, 0.100, 0.250, 0.500, 0.750, 0.900, 0.975];
 
 load latest_us_data.mat
 load us_hyperparam_latest.mat
@@ -111,6 +111,7 @@
     thisdata(all(thisdata==0, 2), :) = [];
     thisdata = diff(thisdata(:, 1:7:num_ahead)')';
     dt = data_4; gt_lidx = size(dt, 2); extern_dat = diff(dt(cid, gt_lidx-7:7:gt_lidx))';
+    extern_dat = [];
     thisdata = [thisdata; repmat(extern_dat, [1 size(thisdata, 2)])];
     thisdata = repmat(thisdata, [5 1]);
     quant_preds_cases(cid, :, :) = movmean(quantile(thisdata, quant_cases)', 5, 1);
@@ -121,6 +122,7 @@
     thisdata(all(thisdata==0, 2), :) = [];
     thisdata = diff(thisdata(:, 1:7:num_ahead)')';
     dt = deaths; gt_lidx = size(dt, 2); extern_dat = diff(dt(cid, gt_lidx-7:7:gt_lidx))';
+    %extern_dat = [];
     thisdata = [thisdata; repmat(extern_dat, [1 size(thisdata, 2)])];
     thisdata = repmat(thisdata, [5 1]);
     quant_preds_deaths(cid, :, :) = movmean(quantile(thisdata, quant_deaths)', 5, 1);
@@ -133,7 +135,7 @@
 
 %% Plot
 cidx = 1:56;
-sel_idx = 27; %sel_idx = contains(countries, 'Washington');
+sel_idx = 3; %sel_idx = contains(countries, 'Rhode Island');
 dt = deaths(cidx, :);
 dts = deaths_s(cidx, :);
 thisquant = squeeze(nansum(quant_preds_deaths(sel_idx, :, [1 12 23]), 1));

matlab scripts/smooth_epidata.m

@@ -1,22 +1,29 @@
-function [data_4_s] = smooth_epidata(data_4, smooth_factor)
+function [data_4_s] = smooth_epidata(data_4, smooth_factor, week_correction)
 %SMOOTH_EPIDATA removes outliers and smoothes
+
+if nargin < 3
+    week_correction = 0;
+end
+
 deldata = diff(data_4');
 deldata(deldata < 0) = 0;
 data_4_s = data_4;
 maxt = size(data_4, 2);
 date_map = ceil(((1:maxt-1) - mod(maxt-1, 7))/7);
 if isnumeric(smooth_factor)
-     cleandel = deldata;
-%     for cid = 1:size(data_4, 1)
-%         week_dat = diff(data_4(cid, mod(maxt-1, 7)+1:7:maxt)');
-%         [clean_week, TF] = filloutliers(week_dat, 'linear', 'movmedian', 20, 'ThresholdFactor', 3);
-%         [~, peak_idx] = findpeaks([week_dat; 0]);
-%         tf_vals = intersect(find(TF), peak_idx);
-%         for jj=1:length(tf_vals)
-%             cleandel(date_map==tf_vals(jj), cid) = clean_week(tf_vals(jj))/7;
-%         end
-%     end
-    deldata = cleandel;
+    cleandel = deldata;
+    if week_correction == 1
+        for cid = 1:size(data_4, 1)
+            week_dat = diff(data_4(cid, mod(maxt-1, 7)+1:7:maxt)');
+            [clean_week, TF] = filloutliers(week_dat, 'linear', 'movmedian', 20, 'ThresholdFactor', 3);
+            [~, peak_idx] = findpeaks([week_dat; 0]);
+            tf_vals = intersect(find(TF), peak_idx);
+            for jj=1:length(tf_vals)
+                cleandel(date_map==tf_vals(jj), cid) = clean_week(tf_vals(jj))/7;
+            end
+        end
+        deldata = cleandel;
+    end
     cleandel = filloutliers(deldata, 'center', 'movmean', smooth_factor*2, 'ThresholdFactor', 3);
     data_4_s = [data_4(:, 1) cumsum(movmean(cleandel', smooth_factor, 2), 2)];
 else