function SchaetzeNowcasting (PDF)
% ------------------------------------------------------------------------------
%  SYNTAX
%SchaetzeNowcasting (PDF);
% ------------------------------------------------------------------------------
%  DESCRIPTION
%    Schätzt die für unsere Corona-Auswertung benötigten Nowcasting-Parameter
% ------------------------------------------------------------------------------
%  INPUT
%    PDF (optional): PDF = 0: keine PDF-Dateien erzeugen
%                    PDF > 0: PDF-Dateien erzeugen
% ------------------------------------------------------------------------------
%  OUTPUT
%    Datei mit den Nowcasting-Parametern für unsere Corona-Auswertung
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%  CHANGE TABLE
%  Date       Name       Modification
%  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%  17.09.2020 JMH        Initiale Version
%  11.12.2020 JMH        Programmkopf ergänzt
%  21.01.2021 JMH        Bildbeschriftungen aufgehübscht
%  20.02.2021 JMH        Nowcasting-Parameter median_Tage >= 0
%  28.02.2021 JMH        Nowcasting-Parameter median_Gesamt >= 0
%  23.04.2021 CN         Matlab-kompatibel gemacht
%  11.05.2021 JMH        Plot 1 aufgehübscht
%  09.10.2021 JMH        Plotten nur bei PDF > 0
%  28.10.2021 JMH        Toleranz gegen Datenausfälle eingebaut
% ------------------------------------------------------------------------------

if nargin < 1
  PDF = 0;
end

%
% Anzahl Tage, für die eine Nowcasting-Korrektur durchgeführt werden soll
%
AnzahlTage = 21;

%
% MAT-Datei mit akkumulierten Ergebnissen von corona3.m einlesen
%
if exist ('/home/jmh', 'dir')
  Pfad = '/bckdsk/tmp/RKI-Download/mat/';
else
  Pfad = '';
end

Dateiname = [Pfad, 'CoronaDatenRKI.mat'];

if exist (Dateiname, 'file') == 2
  load (Dateiname);
else
  SchreibeNowcastingParameter;
  return;
end

%
% Aktuelle (= letzte) RKI-Corona-Datei als Referenz verwenden
%
Referenz.Zeit_d = fix (CoronaDatenRKI (end).Zeit_d     (1:(end - 1)));
Referenz.Kranke =      CoronaDatenRKI (end).Kranke_oNC (1:(end - 1));

%
% Alte RKI-Corona-Dateien wochenweise verarbeiten, Alter: 2 bis 6 Wochen
%
LetzterDatensatz = size (CoronaDatenRKI, 2) - 14;
ErsterDatensatz  = LetzterDatensatz - 27;
Bereich          = ErsterDatensatz:LetzterDatensatz;

%
% Falls zu wenig Daten vorliegen, Standard-Nowcasting-Parameter verwenden
%
if ErsterDatensatz < 1
  SchreibeNowcastingParameter;
  return;
end

%
% Für alte RKI-Corona-Dateien die fehlenden Infizierten bestimmen
%
Gesamt = [];

for Index = Bereich

  Suche = Index;

  while Suche > 0
    try

      Dummy = CoronaDatenRKI (Suche).Zeit_d (end - 1);

      if Suche ~= Index
        CoronaDatenRKI (Index).Zeit_d     = CoronaDatenRKI (Suche).Zeit_d;
        CoronaDatenRKI (Index).Kranke_oNC = CoronaDatenRKI (Suche).Kranke_oNC;
      end

      Suche = 0;

    catch

      Suche = Suche - 7;

    end
  end

  Zeit_d    = ...
    fix (CoronaDatenRKI (Index).Zeit_d     ((end - AnzahlTage - 1):(end - 1)));

  Kranke    = ...
         CoronaDatenRKI (Index).Kranke_oNC ((end - AnzahlTage - 1):(end - 1));

  Ende      = find (Referenz.Zeit_d == Zeit_d (end));
  Soll      = diff (Referenz.Kranke ((Ende - AnzahlTage):Ende));
  Ist       = diff (Kranke);
  Vergleich = (Soll - Ist) ./ Ist;

  Gesamt    = [Gesamt; Vergleich];

end

median_Gesamt = max (median (Gesamt), 0);

median_Tage   = [];
Differenz     = [];

Wochentage    = shift (1:7, (mod ((LetzterDatensatz + 3), 7) + 1));

for Index = Wochentage

  Gesamt_WT    = Gesamt (Index:7:end, :);
  median_WT    = median (Gesamt_WT);
  Differenz_WT = Gesamt_WT - median_WT (ones (size (Gesamt_WT, 1), 1), :);

  median_Tage  = [median_Tage; median_WT];
  Differenz    = [Differenz;   Differenz_WT];

end

median_Tage = round (max (median_Tage, 0) * 1000) / 1000;

SchreibeNowcastingParameter (median_Tage);

%
% Median-Werte der in den RKI-Daten fehlenden Infizierten plotten
%
if PDF > 0

  figure (1);
  clf;

  if exist ('/home/jmh', 'dir')
    set (1, 'position', [7 664 560 420]);
  else
    set (1, 'position', [5  60 640 480]);
  end

  hold off;
  plot ([0, 11], [62, 0],                  'c-',   'LineWidth', 3);
  hold on;
  plot (fliplr (100 * median_Gesamt),      'k-o',  'LineWidth', 2);
  plot (fliplr (100 * median_Tage (1, :)), 'g--x', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Tage (2, :)), 'm--x', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Tage (3, :)), 'b--x', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Tage (4, :)), 'r--+', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Tage (5, :)), 'b--+', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Tage (6, :)), 'g--*', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Tage (7, :)), 'r--*', 'LineWidth', 1.5);
  plot (fliplr (100 * median_Gesamt),      'k-o',  'LineWidth', 2);
  hold off;
  legend ...
    ({'Nowcasting V1: Gerade';                                           ...
      'Nowcasting V2: Median der in den Daten fehlenden Infizierten';    ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Montag';     ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Dienstag';   ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Mittwoch';   ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Donnerstag'; ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Freitag';    ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Samstag';    ...
      'Nowcasting V3: Median der fehlenden Infizierten fuer Sonntag'}, ...
     'Location', 'NorthWest');
  grid ('on');
  set  (gca, 'xTick', 1:AnzahlTage);
  set  (gca, 'xDir',  'reverse');
  set  (gca, 'yTick', 0:10:max (101 * median_Gesamt));
  axis ([1, AnzahlTage, 0, max(101 * median_Gesamt)]);
  Titel  = 'Median-Werte der in den RKI-Daten fehlenden Infizierten';
  Titel  = sprintf ('%s fuer die letzten %d Tage', Titel, AnzahlTage);
  xLabel = ...
    sprintf ('letzte %d Tage (1 = gestern, 2 = vorgestern, ...)', AnzahlTage);
  title  (Titel);
  xlabel (xLabel);
  ylabel ...
    ('In den Daten fehlende Infizierte [%] (100% = unvollstaendige Daten)');

  %
  % Abweichungen von den Median-Werten für alte RKI-Corona-Datensätze plotten
  %
  figure (2);
  clf;

  if exist ('/home/jmh', 'dir')
    set (2, 'position', [521 664 560 420]);
  else
    set (2, 'position', [350 205 640 480]);
  end

  Schwelle  = 0.1;
  Differenz = min (Differenz,  Schwelle);
  Differenz = max (Differenz, -Schwelle);

  Inkrement = fix (size (Differenz, 1) / 7);

  colormap (jet (256));
  imagesc (fliplr (100 * Differenz));
  set  (gca, 'xTick', 1:AnzahlTage);
  set  (gca, 'xDir',  'reverse');
  set  (gca, 'yTick', Inkrement:Inkrement:size (Differenz, 1));
  axis ('xy');
  colorbar;
  Titel  = ...
    'Abweichungen von den Median-Werten fuer alte RKI-Corona-Datensaetze';
  Titel  = sprintf ('%s [%%]', Titel);
  yLabel = ...
    'Index RKI-Datensatz (sortiert nach Wochentagen: 1-4 = Mo, 5-8 = Di, ...)';
  title  (Titel);
  xlabel (xLabel);
  ylabel (yLabel);

  %
  % Mittelwert der in den RKI-Daten fehlenden Infizierten nach Nowcasting
  % plotten
  %
  Gesamt = [];

  for Index = Bereich

    Zeit_d    = ...
      fix (CoronaDatenRKI (Index).Zeit_d     ((end - AnzahlTage - 1):(end-1)));

    Kranke    = ...
           CoronaDatenRKI (Index).Kranke_oNC ((end - AnzahlTage - 1):(end-1));

    Ende      = find (Referenz.Zeit_d == Zeit_d (end));
    Soll      = diff (Referenz.Kranke ((Ende - AnzahlTage):Ende));
    Ist       = diff (Kranke);
    Vergleich = (Soll - Ist) ./ Soll;

    Gesamt    = [Gesamt; Vergleich];

  end

  mean_Gesamt_oNC = mean (Gesamt);

  Gesamt = [];

  figure (3);
  clf;

  if exist ('/home/jmh', 'dir')
    set (3, 'position', [1035 664 560 420]);
  else
    set (3, 'position',  [725  60 640 480]);
  end

  hold off;
  plot (fliplr (100 * mean_Gesamt_oNC), 'g-*', 'LineWidth', 2);
  hold on;
  plot ([NaN, NaN],                     'r-*', 'LineWidth', 2);
  plot ([NaN, NaN],                     'b-*', 'LineWidth', 2);

  for Index = Bereich

    Zeit_d    = ...
      fix (CoronaDatenRKI (Index).Zeit_d     ((end - AnzahlTage - 1):(end-1)));

    Kranke    = ...
           CoronaDatenRKI (Index).Kranke_oNC ((end - AnzahlTage - 1):(end-1));

    Zuwachs   = diff (Kranke);
    Fehlende  = ...
      round (cumsum (median_Tage ((mod ((Index + 3), 7) + 1), :) .* Zuwachs));

    Kranke (2:end) = Kranke (2:end) + Fehlende;

    Ende      = find (Referenz.Zeit_d == Zeit_d (end));
    Soll      = diff (Referenz.Kranke ((Ende - AnzahlTage):Ende));
    Ist       = diff (Kranke);
    Vergleich = (Soll - Ist) ./ Soll;

    Gesamt    = [Gesamt; Vergleich];

    plot (fliplr (100 * Vergleich));

  end

  mean_Gesamt  = mean (Gesamt);
  oben_Gesamt  = mean_Gesamt + std (Gesamt);
  unten_Gesamt = mean_Gesamt - std (Gesamt);

  plot (fliplr (100 * mean_Gesamt),     'r-*',  'LineWidth', 2);
  plot (fliplr (100 * oben_Gesamt),     'b--o', 'LineWidth', 2);
  plot (fliplr (100 * unten_Gesamt),    'b--o', 'LineWidth', 2);
  plot (fliplr (100 * mean_Gesamt_oNC), 'g-*',  'LineWidth', 2);
  hold off;
  legend ...
    ({'Mittelwert der in den Daten fehlenden Infizierten vor Nowcasting'; ...
      'Mittelwert der fehlenden Infizierten nach Nowcasting';             ...
      'Mittelwert +/- Standardabweichung'}, ...
     'Location', 'SouthWest');
  grid ('on');
  set  (gca, 'xTick', 1:AnzahlTage);
  set  (gca, 'xDir',  'reverse');
  set  (gca, 'yTick', -20:10:max (101 * mean_Gesamt_oNC));
  axis ([1, AnzahlTage, -20, max(101 * mean_Gesamt_oNC)]);
  Titel = ...
    'Mittelwert der in den RKI-Daten fehlenden Infizierten nach Nowcasting';
  Titel = sprintf ('%s V3 fuer die letzten %d Tage', Titel, AnzahlTage);
  title  (Titel);
  xlabel (xLabel);
  ylabel ('In den Daten fehlende Infizierte [%] (100% = vollstaendige Daten)');

  %
  % PDF-Dateien erzeugen
  %
  figure (1);
  if exist ('/home/jmh', 'dir')
    PDFFarbplot ('/home/jmh/schrott/tmp/Nowcasting_1', 1);
  end

  figure (2);
  if exist ('/home/jmh', 'dir')
    PDFFarbplot ('/home/jmh/schrott/tmp/Nowcasting_2', 1);
  end

  figure (3);
  if exist ('/home/jmh', 'dir')
    PDFFarbplot ('/home/jmh/schrott/tmp/Nowcasting_3', 1);
  end

end

end

function SchreibeNowcastingParameter (Faktor)

if nargin < 1

  Faktor = [0.017, 0.009, 0.001, 0.005, 0.019, 0.026, 0.025; ...
            0.010, 0.002, 0.004, 0.010, 0.033, 0.027, 0.018; ...
            0.005, 0.006, 0.011, 0.024, 0.034, 0.020, 0.010; ...
            0.008, 0.014, 0.028, 0.026, 0.025, 0.011, 0.005; ...
            0.012, 0.032, 0.034, 0.020, 0.013, 0.005, 0.008; ...
            0.035, 0.044, 0.028, 0.025, 0.006, 0.008, 0.012; ...
            0.054, 0.038, 0.030, 0.010, 0.009, 0.014, 0.037; ...
            0.039, 0.033, 0.018, 0.018, 0.018, 0.040, 0.058; ...
            0.043, 0.027, 0.022, 0.028, 0.055, 0.063, 0.041; ...
            0.039, 0.029, 0.044, 0.057, 0.079, 0.046, 0.045; ...
            0.039, 0.057, 0.082, 0.098, 0.061, 0.050, 0.042; ...
            0.075, 0.111, 0.137, 0.082, 0.083, 0.052, 0.046; ...
            0.146, 0.178, 0.110, 0.113, 0.081, 0.059, 0.087; ...
            0.243, 0.148, 0.162, 0.111, 0.097, 0.100, 0.167; ...
            0.209, 0.211, 0.159, 0.145, 0.173, 0.201, 0.273; ...
            0.276, 0.232, 0.205, 0.225, 0.323, 0.339, 0.242; ...
            0.318, 0.275, 0.297, 0.405, 0.551, 0.302, 0.330; ...
            0.344, 0.351, 0.538, 0.776, 0.456, 0.402, 0.363; ...
            0.497, 0.597, 1.018, 0.644, 0.531, 0.412, 0.354; ...
            0.835, 1.447, 0.729, 0.682, 0.491, 0.401, 0.573; ...
            5.024, 1.725, 1.966, 1.375, 1.714, 1.595, 3.891].';

end

if exist ('/home/jmh', 'dir')
  Pfad = '/home/jmh/octave/corona/';
else
  Pfad = '';
end

Datei  = fopen ([Pfad, 'NowcastingParameter.txt'], 'wt');

Format = '%5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n';

for Index = 1:size (Faktor, 2)
  fprintf (Datei, Format, Faktor (:, Index));
end

fclose (Datei);

for xIndex = 1:size (Faktor, 2)

  if xIndex == 1
    fprintf ('\nFaktor = [');
  else
    fprintf ('          ');
  end

  for yIndex = 1:size (Faktor, 1)

    fprintf ('%5.3f', Faktor (yIndex, xIndex));

    if yIndex < size (Faktor, 1)
      fprintf (', ');
    end

  end

  if xIndex < size (Faktor, 2)
    fprintf ('; ...\n');
  else
    fprintf ('].'';\n\n');
  end

end

end