function ExtrahiereDatenRKI (PDF)
% ------------------------------------------------------------------------------
%  SYNTAX
%ExtrahiereDatenRKI (PDF);
% ------------------------------------------------------------------------------
%  DESCRIPTION
%    Extrahiert die Corona-Daten für unsere Auswertungen aus der RKI-Datei
% ------------------------------------------------------------------------------
%  INPUT
%    PDF (optional): PDF = 0: keine PDF-Dateien erzeugen
%                    PDF > 0: PDF-Dateien erzeugen
% ------------------------------------------------------------------------------
%  OUTPUT
%    BIN-Datei mit den Corona-Daten für unsere Auswertungen
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%  CHANGE TABLE
%  Date       Name       Modification
%  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%  06.08.2020 JMH        Initiale Version
%  08.08.2020 JMH        Konsistenzprüfungen und Korrekturen eingebaut
%  09.08.2020 JMH        Gestorbene und Genesene für aktuellen Tag korrigiert
%  19.08.2020 JMH        Meldedatum statt Erkrankungsdatum (kein Nowcasting!)
%  20.08.2020 JMH        Aktuellen Tag korrigiert
%  21.08.2020 JMH        Neues Download-Verzeichnis für jmh
%  26.08.2020 CN         Matlab-kompatibel gemacht
%  26.08.2020 JMH        Aufgeräumt
%  28.08.2020 JMH        Optionale Parameter mit Auswertemodus und Corona-Datei
%  05.09.2020 JMH        Verbessertes Nowcasting (Version 2)
%  07.09.2020 JMH        Verbesserte Parameter für Nowcasting Version 2
%  07.09.2020 JMH        Verbesserte Extrapolation für Nowcasting Version 2
%  08.09.2020 JMH        Genesene ab 14 Tage nach ihrem Erkrankungsdatum gesund
%  09.09.2020 JMH        Korrigierte Parameter für Nowcasting Version 2
%  10.09.2020 JMH        Verbesserte Parameter für Nowcasting Version 2
%  11.09.2020 JMH        Noch "optimalere" Parameter für Nowcasting Version 2
%  11.09.2020 JMH        Korrigierte Extrapolation für Nowcasting Version 2
%  12.09.2020 JMH        Update der Parameter für Nowcasting Version 2
%  12.09.2020 JMH        Fehlermeldungen verbessert
%  12.09.2020 JMH        Zeitstempel geändert: 00:00 => 18:00
%  14.09.2020 JMH        Nowcasting Version 1 verbessert
%  16.09.2020 JMH        Gestorbene ab 5 Tage nach ihrem Erkrankungsdatum tot
%  17.09.2020 JMH        Verbessertes Nowcasting (Version 3)
%  19.09.2020 JMH        Verbesserte Parameter für Nowcasting Version 3
%  20.09.2020 JMH        Uhrzeit aus Zeitstempel entfernt
%  24.09.2020 JMH        Update der Parameter für Nowcasting Version 3
%  26.09.2020 JMH        Update der Parameter für Nowcasting Version 3
%  01.10.2020 JMH        Verbesserte Parameter für Nowcasting Version 3
%  08.10.2020 JMH        Gestorbene ab 14 Tage nach ihrem Erkrankungsdatum tot
%  08.10.2020 JMH        Update der Parameter für Nowcasting Version 3
%  17.10.2020 JMH        Update der Parameter für Nowcasting Version 3
%  22.10.2020 JMH        y-Achsenbeschriftung der Histogramme mit Umlaut ä
%  24.10.2020 JMH        Update der Parameter für Nowcasting Version 3
%  25.10.2020 JMH        Median- und Mittelwert in Histogramme eingezeichnet
%  28.10.2020 JMH        Aktuelles Datum der RKI-Datei aus DAT-Datei einlesen
%  31.10.2020 JMH        Update der Parameter für Nowcasting Version 3
%  03.11.2020 JMH        Unstetigkeit am Ende der Genesenen und Toten reduziert
%  05.11.2020 JMH        Update der Parameter für Nowcasting Version 3
%  21.11.2020 JMH        Update der Parameter für Nowcasting Version 3
%  21.11.2020 JMH        Kranke_MD und Kranke_oNC ergänzt
%  21.11.2020 JMH        Verbesserungen von CN übernommen
%  21.11.2020 JMH        RKI-Datensätze nicht mehr nach Zeitstempeln sortieren
%  21.11.2020 JMH        Fortschrittsbalken eingebaut
%  21.11.2020 JMH        Dauerkranke nicht mehr anzeigen
%  21.11.2020 JMH        inPfad geändert
%  25.11.2020 JMH        Modus-Fehler korrigiert
%  26.11.2020 JMH        Update der Parameter für Nowcasting Version 3
%  03.12.2020 JMH        Update der Parameter für Nowcasting Version 3
%  10.12.2020 JMH        Update der Parameter für Nowcasting Version 3
%  10.12.2020 JMH        Modus und CoronaDatei eliminiert
%  23.01.2021 JMH        Histogramme in PDF-Dateien ausgeben
%  25.02.2021 JMH        Corona-Daten sollen mit dem 01.01.2020 beginnen
%  26.02.2021 JMH        Diagnoseverzug-Histogramme korrigiert
%  23.04.2021 CN         Histogramme aufgehübscht
%  24.04.2021 JMH        Aufgeräumt
%  03.05.2021 JMH        Auswertung der Altersgruppen ergänzt
%  17.05.2021 JMH        Hauptschleife in MEX-Datei, um Rechenzeit zu verkürzen
%  09.10.2021 JMH        Plotten nur bei PDF > 0
%  24.11.2021 JMH        Einlesen der RKI-Datei beschleunigt
%  01.12.2021 JMH        Binärdatei einlesen anstelle der beiden Textdateien
% ------------------------------------------------------------------------------

disp (' ');
disp ('Auswertung mit Erkrankungsdatum und Nowcasting sowie Alters-Auswertung');

if nargin < 1
  PDF = 0;
end

%
% Pfade bestimmen
%
if exist ('/home/jmh', 'dir')
  inPfad  = '/bckdsk/tmp/RKI-Download/';
  outPfad = '/home/jmh/octave/corona/';
  agePfad = '/home/jmh/octave/corona/Altersgruppen/';
else
  inPfad  = '';
  outPfad = '';
  agePfad = '';
end

%
% Dateinamen bestimmen
%
DateinameRKI = 'RKI_COVID19.bin';
DateinamePAR = 'NowcastingParameter.txt';
DateinameTXT = 'CoronaDatenRKI.txt';
DateinameAGE = 'CoronaDatenRKI_Alter.txt';

inDateiname  = [inPfad,  DateinameRKI];
parDateiname = [outPfad, DateinamePAR];
outDateiname = [outPfad, DateinameTXT];
ageDateiname = [agePfad, DateinameAGE];

%
% BIN-Datei mit den Corona-Daten vom RKI einlesen
%
inDatei      =  fopen (inDateiname, 'rb');

AnzahlFelder =  fread (inDatei, 1, 'int32');
AnzahlZeilen =  fread (inDatei, 1, 'int32');
Datum        = (fread (inDatei, 3, 'int32')).';

Daten        = (fread (inDatei, [AnzahlFelder, AnzahlZeilen], 'int16')).';

fclose (inDatei);

%
% Die Zeitachse ist in Tagen skaliert: 01.01.2020 = Tag Nr. 1
%
Tag_0 = datenum (2019, 12, 31);

%
% Aktuellen Tag bestimmen
%
%Datenstand: Datum, wann der Datensatz zuletzt aktualisiert worden ist
aktuellerTag = berechneZeit (Datum (3), Datum (2), Datum (1), 0, 0, Tag_0);

%
% Zeitstempel der RKI-Datensätze berechnen
%

% Meldedatum: Datum, wann der Fall dem Gesundheitsamt bekannt geworden ist
% (Stunden, Minuten und Sekunden sind immer Null!)
Meldezeit_d = ...
  berechneZeit (Daten (:, 3), Daten (:, 4), Daten (:,  5), 0, 0, Tag_0);

% Referenzdatum:
% Erkrankungsdatum bzw. wenn das nicht bekannt ist, das Meldedatum
% (Stunden, Minuten und Sekunden sind immer Null!)
Zeit_d      = ...
  berechneZeit (Daten (:, 8), Daten (:, 9), Daten (:, 10), 0, 0, Tag_0);

% IstErkrankungsbeginn: 1, wenn das Refdatum der Erkrankungsbeginn ist, 0 sonst
% (Bei IstErkrankungsbeginn = 0 sind Referenzdatum und Meldedatum identisch!)
istMeldedatum = find (Daten (:, 13) == 0);
if any (Zeit_d (istMeldedatum) ~= Meldezeit_d (istMeldedatum))
  disp  (Daten (istMeldedatum, [3, 4, 5, 8, 9, 10]));
  error ('%s: any (Zeit_d (istMeldedatum) ~= Meldezeit_d (istMeldedatum))', ...
         mfilename);
end

%
% Histogramme der Differenz von Meldedatum und Erkrankungsdatum (Diagnoseverzug)
%
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
  istErkrankungsbeginn = find (Daten (:, 13) == 1);
  AnzahlFall           = (Daten (istErkrankungsbeginn, 1)).';
  AnzahlFall_max       = max (AnzahlFall);
  Meldedatum           = Meldezeit_d (istErkrankungsbeginn);
  Erkrankungsdatum     = Zeit_d (istErkrankungsbeginn);
  Differenz            = [];
  for Index = 1:AnzahlFall_max
    Differenz        = [Differenz, (Meldedatum - Erkrankungsdatum)];
    AnzahlFall       = AnzahlFall - 1;
    groesserNull     = find (AnzahlFall > 0);
    AnzahlFall       = AnzahlFall       (groesserNull);
    Meldedatum       = Meldedatum       (groesserNull);
    Erkrankungsdatum = Erkrankungsdatum (groesserNull);
  end
  Medianwert = median (Differenz);
  Mittelwert = mean   (Differenz);
  Minimum    = -4;
  Maximum    = 30;
  [NN, XX]   = hist (Differenz, Minimum:Maximum, 100);
  hold off;
  bar (XX, NN, 1, 'facecolor', 'r', 'edgecolor', 'b');
  hold on;
  yLim = [0, (max(NN) + 0.5)];
  set  (gca, 'yLim', yLim);
  plot ([Medianwert, Medianwert], yLim, '-b', 'LineWidth',2);
  plot ([Mittelwert, Mittelwert], yLim, '-g', 'LineWidth',2);
  hold off;
  strMedianwert = sprintf ('Medianwert = %.1f Tage', Medianwert);
  strMittelwert = sprintf ('Mittelwert = %.1f Tage', Mittelwert);
  legend ({'Histogramm'; strMedianwert; strMittelwert});
  set    (gca, 'xLim', [(Minimum - 0.5), (Maximum + 0.5)]);
  grid ('on');
  title  ('Histogramm der Differenz von Meldedatum und Erkrankungsdatum');
  xlabel ('Diagnoseverzug [Tage]');
  ylabel ('Relative Haeufigkeit [%]');
  drawnow;

  figure (2);
  clf;
  if exist ('/home/jmh', 'dir')
    set (2, 'position', [521 664 560 420]);
  else
    set (2, 'position', [350 205 640 480]);
  end
  NN = cumsum (NN);
  hold off;
  bar (XX, NN, 1, 'facecolor', 'r', 'edgecolor', 'b');
  hold on;
  plot ([Medianwert, Medianwert], [0, 100], '-b', 'LineWidth',2);
  plot ([Mittelwert, Mittelwert], [0, 100], '-g', 'LineWidth',2);
  hold off;
  legend ({'Histogramm'; strMedianwert; strMittelwert}, ...
          'Location', 'SouthEast');
  axis   ([(Minimum - 0.5), (Maximum + 0.5), 0, 100]);
  grid ('on');
  title  ('Summenhistogramm der Differenz von Meldedatum und Erkrankungsdatum');
  xlabel ('Diagnoseverzug [Tage]');
  ylabel ('Akkumulierte relative Haeufigkeit [%]');
  drawnow;

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

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

end

%
% Corona-Daten für unsere Auswertung berechnen
%
if min (Zeit_d) < 1
  disp  (Daten (1, 8:10));
  error ('%s: min (Zeit_d) < 1', mfilename);
end

if max (Zeit_d) >= aktuellerTag
  disp  (Daten (end, 8:10));
  error ('%s: max (Zeit_d) >= aktuellerTag', mfilename);
end

Tag                   = 1:aktuellerTag;
Tag_Alter             = 1:(aktuellerTag - 1);

DatumUhrzeit          = berechneDatum ((Tag       + 18 /  24),       Tag_0);
DatumUhrzeit (end, :) = berechneDatum ((Tag (end) -  1 / (24 * 60)), Tag_0);

[Krankdauer_d, Kranke, Kranke_MD, Gesund, Tote, Kranke_Alter, Tote_Alter] = ...
  sortiere_RKI_Datensaetze_ein (length (Tag), Daten, Zeit_d, Meldezeit_d);

%
% Unstetigkeit am Ende der Genesenen und Toten reduzieren
%
dGesund = cumsum (fliplr (diff (Gesund (length (Kranke):end))));
dTote   = cumsum (fliplr (diff (Tote   (length (Kranke):end))));

dGesund = [dGesund, dGesund(end)];
dTote   = [dTote,   dTote(end)];

Gesund  = Gesund (1:length (Kranke));
Tote    = Tote   (1:length (Kranke));

Gesund ((end - Krankdauer_d):end) = Gesund ((end - Krankdauer_d):end) + dGesund;
Tote   ((end - Krankdauer_d):end) = Tote   ((end - Krankdauer_d):end) + dTote;

%
% Unstetigkeit am Ende der Toten reduzieren
%
dTote      = flipud (Tote_Alter ((length (Kranke_Alter) + 1):end, :));

Tote_Alter = Tote_Alter (1:length (Kranke_Alter), :);

Tote_Alter ((end - Krankdauer_d + 1):end, :) = ...
  Tote_Alter ((end - Krankdauer_d + 1):end, :) + dTote;

%
% Auswertung mit Erkrankungsdatum ohne Nowcasting
%
Kranke_oNC = Kranke;

%
% Nowcasting-Parameter bestimmen und einlesen
%
SchaetzeNowcasting (PDF);

Faktor = (load (parDateiname)).';

%
% Auswertung mit Erkrankungsdatum und Nowcasting
%
Wochentag = mod ((aktuellerTag + 1), 7) + 1;
Anzahl    = size (Faktor, 2);

Zuwachs   = diff (Kranke ((end - Anzahl - 1):(end - 1)));
Fehlende  = round (cumsum (Faktor (Wochentag, :) .* Zuwachs));

Kranke ((end - Anzahl):(end - 1)) = ...
  Kranke ((end - Anzahl):(end - 1)) + Fehlende;

%
% Neuinfizierte für aktuellen Tag
%
Kranke (end) = Kranke (end - 1);

%
% Corona-Daten sollen mit dem ersten Kranken beginnen
%
%Tag = Tag (find (Kranke > 0));

%
% TXT-Datei mit den Corona-Daten für unsere Corona-Auswertung speichern
%
outDatei = fopen (outDateiname, 'wt');

for Index = Tag
  fprintf (outDatei, '%4d %02d %02d %02d %02d ', DatumUhrzeit (Index, :));
  fprintf (outDatei, '%8d ',                     Kranke       (Index));
  fprintf (outDatei, '%8d ',                     Gesund       (Index));
  fprintf (outDatei, '%8d ',                     Tote         (Index));
  fprintf (outDatei, '%8d ',                     Kranke_MD    (Index));
  fprintf (outDatei, '%8d\n',                    Kranke_oNC   (Index));
end

fclose (outDatei);

%
% TXT-Datei mit den Corona-Daten für die Auswertung der Altersgruppen speichern
%
ageDatei = fopen (ageDateiname, 'wt');

for Index = Tag_Alter
  fprintf (ageDatei, '%4d %02d %02d ',            DatumUhrzeit (Index, 1:3));
  fprintf (ageDatei, '%8d %8d %8d %8d %8d %8d',   Kranke_Alter (Index, :));
  fprintf (ageDatei, '%8d %8d %8d %8d %8d %8d\n', Tote_Alter   (Index, :));
end

fclose (ageDatei);

end

function Zeit_d = berechneZeit (Jahr, Monat, Tag, Stunde, Minute, Tag_0)

Zeit_d = (datenum (Jahr, Monat, Tag) + (Stunde + Minute / 60) / 24 - Tag_0).';

end

function DatumUhrzeit = berechneDatum (Zeit_d, Tag_0)

DatumUhrzeit = zeros (numel (Zeit_d), 5);

for Index = 1:numel (Zeit_d)

  % Runden auf ganze Minuten:
  Vektor = datevec (round (24 * 60 * (Zeit_d (Index) + Tag_0)) / (24 * 60));

  DatumUhrzeit (Index, :) = Vektor (1:5);

end

end