function out = COVID19Modelingv2(country,daysforecast,clearfig,gotopast)
%COVID19MODELINGV2 fits the last COVID-19 cases (updated daily) for one country, a selection of entities or all entities (169) and proposes forecasts
% FORECAST results are based on the assumption of logistic distribution, they are proposed "AS IS"
% The logistic distribution looks as the normal distribution in shape but has heavier tails (higher kurtosis), see https://en.wikipedia.org/wiki/Logistic_regression
% (when the dataset is too small or when the numbers of cases did not enter in the % exponential phase, the statistics is very poor)
%
% Syntax: COVID19Modeling(country,daysforecast)
%   Default country = "selection", where selection is a list of 24 entities
%   (selection = ["France" "Italy" "Spain" "Germany" "UK" "Belgium" "Portugal" "Croatia" "Denmark" "US" "Canada" "Australia" "China" "India" "Japan" "Korea, South" "Iran" "Russia" "Lebanon" "Netherlands" "Austria" "Denmark" "Turkey" "Israel"])
%   Default number days of forecast = 14
%
%   LIST ALL entitIES/PROVINCES/REGIONS/DEPARTMENTS AVAILABLE (the list depends on database updates)
%   CURRENT DATABASES CONTAIN:   **180** coutries and up to **428** sub-regions.
%          COVID19Modelingv2('list')
%
%   Examples:
%       COVID19Modeling("France") or COVID19Modeling('France')
%       COVID19Modeling(["France" "Italy"]) or COVID19Modeling({'France' 'Italy'})
%       COVID19Modeling() or COVID19Modeling("selection") --> a selection of 24 entities
%            COVID19Modeling("all") or COVID19Modeling('all') --> do the analysis for all entities in the list (169)
%
% ==== GLOBAL DATA SOURCE
%   2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository by Johns Hopkins CSSE
%   https://github.com/CSSEGISandData/COVID-19
% Overview
%   Novel Corona Virus (COVID-19) epidemiological data since 22 January 2020.
%   Data are compiled by the Johns Hopkins University Center for Systems Science and Engineering (JHU CCSE)
%   from various sources including the World Health Organization (WHO), DXY.cn. Pneumonia. 2020, BNO News, National 
%   Health Commission of the People’s Republic of China (NHC), China CDC (CCDC), Hong Kong Department of Health, Macau Government, 
%   Taiwan CDC, US CDC, Government of Canada, Australia Government Department of Health, European Centre for Disease Prevention and Control (ECDC), 
%   Ministry of Health Singapore (MOH). JSU CCSE maintains the data on the 2019 Novel Coronavirus COVID-19 (2019-nCoV)
%   Data Repository on github: https://github.com/CSSEGISandData/COVID-19
%
% ==== FRENCH DATA SOURCE
%   Organisation informelle issue de la société civile dont l'objet est de consolider des données et de proposer des outils de visualisation concernant l'épidémie de COVID19 en France.
%   https://www.data.gouv.fr/fr/organizations/opencovid19-fr/
%   https://github.com/opencovid19-fr
%                       |
%                       v
%   +-----------------------------------------+
%   | The section below keeps tracks of the   |
%   | many updates occurring with time. The   |
%   | outbreak is very active, so the changes |
%   | are only briefly reported.              |
%   +-----------------------------------------+
%
% ==== PROJECT HISTORY ===== (the current code is a fork of a fork)
%   Coronavirus Tracker - Country Tracker - Joshua McGee
%   Created to track the simulate the spread of Coronavirus (COVID-19)
%   Data is stored online and is provided via JHU CSSE from various sources including:
%   "the World Health Organization (WHO), DXY.cn. Pneumonia. 2020, BNO News,
%   National Health Commission of the People? Republic of China (NHC),
%   China CDC (CCDC), Hong Kong Department of Health, Macau Government, Taiwan CDC, US CDC,
%   Government of Canada, Australia Government Department of Health,
%   European Centre for Disease Prevention and Control (ECDC) and Ministry of
%   Health Singapore (MOH)"
%
% ==== DATA HISTORY ====
%   Three new time series tables will be added for the US. The first two will be the confirmed cases and deaths,
%   reported at the county level. The third, number of tests conducted, will be reported at the state level.
%   These new tables will be named time_series_covid19_confirmed_US.csv, time_series_covid19_deaths_US.csv, 
%   time_series_covid19_testing_US.csv, respectively.
% 
%   Changes to the current time series include the removal of the US state and county-level entries, which will be
%   replaced with a new single country level entry for the US. The tables will be renamed time_series_covid19_confirmed_global.csv
%   and time_series_covid19_deaths_global.csv, and time_series_covid19_testing_global.csv, respectively.
%   The ISO code will be added in the global time series tables.
%   The FIPS code will be added in the new US time series tables.
%
% 
% Fields available in the data include Province/State, Country/Region, Last Update, Confirmed, Suspected, Recovered, Deaths.
% 
% On 23/03/2020, a new data structure was released. The current resources are:
% 
% time_series_covid19_confirmed_global.csv
% time_series_covid19_deaths_global.csv for the latest time series data
% ---DEPRECATION WARNING---
% The resources below ceased being updated on 22/03/2020:
% 
% time_series_19-covid-Confirmed.csv
% time_series_19-covid-Deaths.csv
% time_series_19-covid-Recovered.csv
%
%
% ============= LIST OF REVISIONS BY Olivier Vitrac (INRAe)
%   2020-03-20 fork from https://fr.mathworks.com/matlabcentral/fileexchange/74632-covid19-modeling-fitvirus-update (version 2.1)
%   2020-03-20 syntax COVID19Modeling(country,daysforecast)
%   2020-03-20 print PDF, PNG, FIG France Data by defaut (remove DOM-TOM)
%   2020-03-21 PREFETCH file for acceleration
%   2020-03-21 Print the predictions in a text file
%   2020-03-21 fix DOM-TOM, run for a selection a entities by default and store the results
%   2020-03-21 embbed internal functions for printing (low level ones)
%   2020-03-22 remove empty folder, fix warnings and other bugs, expand help
%   2020-03-23 add forecast to text report
%   2020-03-23 add death toll projections to the figure (two estimates)
%   2020-03-23 add analysis for each state/province (China, US, Canada, Australia, UK, France)
%   2020-03-24 fix the calculation of death rate based on maximum, fix country,state name (e.g. France, France)
%   2020-03-25 Major update to follow the new datasets from Johns Hopkins University Center for Systems Science and Engineering
%   2020-03-25 add maps (with good resolution)
%   2020-03-26 change map layout (force maximized) for detailed view
%   2020-03-26 fix how death toll are represented on the map (max lowered from 50000 down to 20000, less scary)
%   2020-03-28 [MAJOR UPDATE] full implementation of US data (grouped by state) // Data start from March 10, 2020 (before cases were counted very differently)
%   2020-03-29 [MAJOR UPDATE] add France details data France/Departements (96 entries) and France/Regions (13 entries), please DOM-TOM are stored in the main folder simulation/ 
%   2020-03-30 fix delay in updating France results, causing the last data (zeros) to be taken and finally the loss of the map
%   2020-03-31 add specific display when World, US, French databases are updated (it was required to identify when data are updated)
%   2020-03-31 check names, elapsed and remaining time
%   2020-04-01 add confidence intervals to cases, cases/day, prepare animations
%   2020-04-02 animations, more robust confidence interval calculations
%   2020-04-02 convert AVI into light MP4 videos (h264), fix animations (scales and colors)
%   2020-04-03 new global sources (unification of sources), pepare MAJOR UPDATE
%   2020-04-04 prerelease, a lot of nomenclature to fix
%   2020-04-05 [MAJOR UPDATE] v 2.60, full consistency check based on daily data
%   2020-04-05 h264 videos are replaced by h265 (MP4) ones, along VP9 (WEBM) ones to maximze compatibility with browsers
%   2020-04-06 fix US maps, cities and ships were without Lat/Lon data (such data have been removed, kinetics and maps are useless now)
%   2020-04-11 interrupts COVID19Modelingv2() when 1/b(2) is too large, unrealistic time constant (threshold: 25 years)
%   2020-04-15 several consolidations for French data (EHPAD data are now correctly managed)
%              NOTE: EHPAD data are counted only at national level, no geographical information
%   2020-04-18 upper bound of death rate in maps increased up to 40000
%   2020-04-22 discard new entries US/Federal Bureau of Prisons, US/US Military, US/Veteran Hospitals added on April, 21
%   2020-04-25 use raw global date with using the proxy data.humdate.org
%              [MAJOR UPDATE] v2.70, compatibility with Matlab R2020a, automatic color scheme for death toll in maps
%   2020-04-28 add structure output (array if several countries), change color, add ParisSubur, fix 'List'
%   2020-05-15 add Germany/Unknown Lat,Lon with NaN, NaN
%   2020-05-29 fix Port Quarantine in Japan
%   2020-11-05 remove missing entries in tmp
%   2020-11-13 fix date in France national results (v 2.76)
%   2020-11-28 fix Repatriated Travellers,Canada - starting from 2020-11-24 (v. 2.77)
%   2021-02-25 fix unknown for China (v. 2.771)

% ===================== Modifications OV
% current version
versn = 2.771;  % numbering rule (2 because v2)
% Original data source (see https://data.humdata.org/dataset/5dff64bc-a671-48da-aa87-2ca40d7abf02) - 
% Starting from April, 3rd, original data started to be corrupted (e.g. mising data for an entire country such as France)
% As it was initiated for US starting from April, 23, raw data are used instead
% global source: https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data/csse_covid_19_daily_reports
% URL below are still load for back compatibility between formats but the data are originating from raw data
% Before April 25th, 2020 (using official proxy)
% resultURL = 'https://data.humdata.org/hxlproxy/api/data-preview.csv?url=https%3A%2F%2Fraw.githubusercontent.com%2FCSSEGISandData%2FCOVID-19%2Fmaster%2Fcsse_covid_19_data%2Fcsse_covid_19_time_series%2Ftime_series_covid19_confirmed_global.csv&filename=time_series_covid19_confirmed_global.csv';
% deathURL = 'https://data.humdata.org/hxlproxy/api/data-preview.csv?url=https%3A%2F%2Fraw.githubusercontent.com%2FCSSEGISandData%2FCOVID-19%2Fmaster%2Fcsse_covid_19_data%2Fcsse_covid_19_time_series%2Ftime_series_covid19_deaths_global.csv&filename=time_series_covid19_deaths_global.csv';
% After April 25th, 2020 (John Hoppkins updated lately its data for Fri 24th on Sat 25th)
resultURL = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv';
deathURL = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv';

% RAW data: CSSE COVID-19 Dataset
% https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data
% File naming convention: MM-DD-YYYY.csv in UTC.
% Field description
% Province/State: China - province name; US/Canada/Australia/ - city name, state/province name; Others - name of the event (e.g., "Diamond Princess" cruise ship); other entities - blank.
% Country/Region: country/region name conforming to WHO (will be updated).
% Last Update: MM/DD/YYYY HH:mm (24 hour format, in UTC).
% Confirmed: the number of confirmed cases. For Hubei Province: from Feb 13 (GMT +8), we report both clinically diagnosed and lab-confirmed cases. For lab-confirmed cases only (Before Feb 17), please refer to who_covid_19_situation_reports. For Italy, diagnosis standard might be changed since Feb 27 to "slow the growth of new case numbers." (Source)
% Deaths: the number of deaths.
% Recovered: the number of recovered cases
URLrawresultGLOBAL = @(d) sprintf('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_daily_reports/%s.csv',datestr(d,'mm-dd-yyyy'));


% FLAGS for debugging and production
% for debugging (one country), use noglobalprefetch = true and noprefetch = true;
% for production, use noglobalprefetch = flase and noprefetch = true;
noglobalprefetch = false; % set it to true to avoid prefetch
noprefetch = true; % set it to true to avoid country prefetch (usually fast enven if true)
nodailyreportprefetch = false; % force the daily report to be read again (mainly for debugging or for implementing new rules)
GIFon = true; % generate map animations (GIF and MP4)
offline = false; % offline debugging (github went off-line, see https://www.githubstatus.com/incidents/80d0cs6kpsps)
listofmaps = {'World' 'Asia' 'America' 'Europe' 'France'};
mapon = true;

% arg check
if nargin<1, country = ""; end
if nargin<2, daysforecast = []; end
if nargin<3, clearfig = false; end
if nargin<4 || isempty(gotopast), gotopast = 0; end

% pseudo recursion to go several days in the past
% to be used to renegerate missing data on the server
% set today = today() in base, if you launch calculations crossing midnight
% EXAMPLE: today=today(); COVID19Modelingv2('all',[],true,1:4)
if any(gotopast<0), error('gotopast must be positive or zero'), end
if length(gotopast)>1
    for i=1:length(gotopast)
        COVID19Modelingv2(country,daysforecast,clearfig,gotopast(i))
    end
    return
end

% Date initialization (based on base), added on April 5h, 2020
today = evalin('base','today') - gotopast; % today in 'base' can be either a variable or the function today()

% DEFAULT values
cT0 = clock; % clock
selection = ["France" "Italy" "Spain" "Germany" "UK" "Belgium" "Portugal" "Croatia" "Denmark" "US" "Canada" "Australia" "China" "India" "Japan" "Korea, South" "Iran" "Russia" "Lebanon" "Netherlands" "Austria" "Turkey" "Israel"];
parissuburb = [ "France" "France/Regions/Région Île-de-France" "France/Departements/DPT 91,Essonne" "France/Departements/DPT 92,Hauts-de-Seine" "France/Departements/DPT 93,Seine-Saint-Denis" "France/Departements/DPT 94,Val-de-Marne" "France/Departements/DPT 77,Seine-et-Marne" "France/Departements/DPT 75,Paris"];
country_default = "France";
dp_default = 14; % forecasting
rootfolder = pwd;
datafolder = 'data';
outputfolder = 'simulation';
mapfolder = 'map';
globaloutputpath = fullfile(rootfolder,outputfolder); % default in first version, then the name of the country is appended
outputpath = globaloutputpath;
mappath = fullfile(rootfolder,mapfolder);
datetxt = datestr(today,'yyyy-mm-dd');
yesterday = today - 1;
deathratemax = 0.25; % to prevent unrealistic death rates (added 2020-03-30)
pauseduration = 0.25; % pause duration during update

% VIDEO converters from GIF (compatible iOS, Chrome, Firefox, Opera)
% Source: https://medium.com/abraia/video-transcoding-and-optimization-for-web-with-ffmpeg-made-easy-511635214df0
% Example (to refresh all video rapidly)
%gifs = explore('*.gif','map',0,'fullabbreviate'); mp4s = regexprep(gifs,'.gif$','.mp4'); webms = regexprep(gifs,'.gif$','.webm'); delete(mp4s{:}); cellfun(gif2mp4,gifs,mp4s); delete(webms{:}); cellfun(gif2webm,gifs,webms);
gif2mp4 =  @(gif,mp4)  system(sprintf('ffmpeg -f gif -i ''%s'' -c:v libx265 -crf 23 -tag:v hvc1 -pix_fmt yuv420p -color_primaries 1 -color_trc 1 -colorspace 1 -movflags +faststart -an ''%s''',gif,mp4));
gif2webm = @(gif,webm) system(sprintf('ffmpeg -f gif -i ''%s'' -c:v libvpx-vp9 -crf 30 -speed 3 -pix_fmt yuv420p -color_primaries 1 -color_trc 1 -colorspace 1 -movflags +faststart -an ''%s''',gif,webm));


% GLOBAL prefetch
listvars = {'infectedtable','countrytable','countrytable1','Countrytotalinfected','Countrytotaldead','infected','deathrate'};
for d={datafolder outputfolder mapfolder}
    if ~exist(fullfile(rootfolder,d{1}),'dir'), mkdir(rootfolder,d{1}); end
end
globalprefetchfile = fullfile(rootfolder,datafolder,sprintf('PREFETCH_COVID_%s.mat',datetxt));
if exist(globalprefetchfile,'file'), load(globalprefetchfile,'listentities'); else, listentities={}; end

% CHECK countries
if ischar(country) || iscellstr(country), country=string(country); end
if  isempty(country) || (length(country)==1 && ((country == "") || ismissing(country))), country = country_default; end
% fix country names // added on 2020-03-31 (check conversion with country=listentities; country(country~=listentities)
% this update enables the user to enter a country without kowing the case (files are case senstive)
country = regexprep(country,...
    {'(^.)(.*$)'               ,'([\s\,\/].)' ,'(/DPT|^US|^Uk|\(.)','-(.)','\,(\s?.)','''(.)'},...
    {'${upper($1)}${lower($2)}','${upper($1)}','${upper($1)}'      ,'-${upper($1)}',',${upper($1)}','''${upper($1)}'},...
    'ignorecase');
prepositions = {' The ',' And ',' Of ',' Le ',' La ',' Et ','-Et-',' Du ','-Du-',' De ','-De-','D'''};
country = regexprep(country,prepositions,lower(prepositions),'ignorecase');
exceptions = {'UKraine' 'Ukraine';'DPT 2a','DPT 2A';'DPT 2B','DPT 2B';'Ms Zaandam','MS Zaandam';'Dom-Tom','DOM-TOM'};
country = regexprep(country,exceptions(:,1),exceptions(:,2),'ignorecase');
% end fix entities (all 429 rules)
if isempty(daysforecast), daysforecast = dp_default; end
if lower(country)=="selection", country = selection; end
if lower(country)=="parissuburb", country = parissuburb; end
if lower(country)=="all"
    if ~exist(globalprefetchfile,'file') || noglobalprefetch
        fprintf('\n\n *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* \n\n')
        fprintf('\n\n Initialization of the database with de default country ("%s") \n\n',country_default)
        fprintf('\n\n *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* \n\n')
        COVID19Modelingv2(country_default,daysforecast,clearfig,gotopast)
        fprintf('\n\n *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* \n\n')
        load(globalprefetchfile,'listentities')
    end
    country = listentities; 
end
if isempty(country), country = country_default; end
if length(country)>1 % pseudo-recursion (if several entities are requested)
    nentities = length(country);
    if nentities>10, close all; clearfig = true; end
    for icountry = 1:nentities
        fprintf('\n\n[%d/%d] %s (v. %0.3g) analyzing %s\n',icountry,nentities,mfilename,versn,country(icountry))
        try
            outtmp = COVID19Modelingv2(country(icountry),daysforecast,clearfig,gotopast);
        catch
            fprintf('Last error:\n%s\n',lasterr); %#ok<LERR>
            fprintf('no convergence for %s\n',country(icountry))
            countrypath(country(icountry))
            cleancountrypath()
        end
        done = 100*icountry/nentities; elapsed = etime(clock,cT0); remaining = (100-done)*elapsed/done;
        if icountry<nentities
            fprintf('\n\n[done=%0.3g%%, elapased time=%0.4g s, remaining time %0.4g s] %s analysed (%d of %d)\n',done,elapsed,remaining,country(icountry),icountry,nentities);
        end
        if nargout
            if icountry==1, out=repmat(outtmp,nentities,1); end
            out(icountry) = outtmp;
        end
    end
    fprintf('\n\n[done=%0.3g%%, elapased time=%0.4g s] %d entities or regions have been analysed\n',done,elapsed,nentities);
    return
end

% list entities if needed
if any(lower(country)=="list"), showentitiesintab(); return, end

% output and prefetch files
countrypath()
datafile = @(f,e) fullfile(rootfolder,datafolder,sprintf('COVID_%s_%s.%s',datetxt,f,e));
outputfile = sprintf('COVID_%s_%s',regexprep(country,'^.*/',''),datetxt);
prefetchfile = fullfile(rootfolder,datafolder,sprintf('PREFETCH_COVID_%s_%s.mat',regexprep(country,'/','-'),datetxt));

% GLOBAL RESULTS
GLOBALrawresultfile = @(d) fullfile(rootfolder,datafolder,sprintf('COVID_GLOBAL_rawresult_%s.csv',datestr(d,'yyyy-mm-dd')));
GLOBALprefetchfile = fullfile(rootfolder,datafolder,'PREFETCH_COVID_GLOBAL_rawdata.mat'); % specific prefetch for global data 
GLOBALfirstdate = datenum('01-22-2020','mm-dd-yyyy'); % first date in COVID-19 daily reports

% detailed FRANCE data (ongoing work, not fully satisfactory, gaps)
% Source: https://github.com/opencovid19-fr/data/blob/master/dist/chiffres-cles.csv
URLrawresultFR = 'https://raw.githubusercontent.com/opencovid19-fr/data/master/dist/chiffres-cles.csv';
FRrawresultfile = fullfile(rootfolder,datafolder,'COVID_FR_rawresult.csv');
FRprefetchfile = fullfile(rootfolder,datafolder,sprintf('PREFETCH_COVID_FR_rawdata_%s.mat',datetxt)); % specific prefetch for France

% workaround for USA data (details by state and county removed from the main stream)
%URLpreviewresultUS = sprintf('https://github.com/CSSEGISandData/COVID-19/blob/master/csse_covid_19_data/csse_covid_19_daily_reports/%s.csv',datestr(yesterday,'mm-dd-yyyy'));
URLrawresultUS = @(d) sprintf('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_daily_reports/%s.csv',datestr(d,'mm-dd-yyyy'));
USrawresultfile = @(d) fullfile(rootfolder,datafolder,sprintf('COVID_USA_rawresult_%s.csv',datestr(d,'yyyy-mm-dd')));
USprefetchfile = fullfile(rootfolder,datafolder,'PREFETCH_COVID_USA_rawdata.mat'); % specific prefetch for US
USfirstdate = datenum('03-10-2020','mm-dd-yyyy'); % before this date, data are per county, not state

%% Obtaining and formating data - courtesy of Toshi Takeuchi - https://www.mathworks.com/matlabcentral/profile/authors/951521
if ~exist(prefetchfile,'file') || noprefetch
    if ~exist(globalprefetchfile,'file') || noglobalprefetch
        %% global prefetch
        resultfile = datafile('result','txt');
        deathresultfile = datafile('deathresult','txt');
        % depreciated: see https://data.humdata.org/dataset/5dff64bc-a671-48da-aa87-2ca40d7abf02
        %result=webread('https://proxy.hxlstandard.org/api/data-preview.csv?url=https%3A%2F%2Fraw.githubusercontent.com%2FCSSEGISandData%2FCOVID-19%2Fmaster%2Fcsse_covid_19_data%2Fcsse_covid_19_time_series%2Ftime_series_19-covid-Confirmed.csv&filename=time_series_2019-ncov-Confirmed.csv','options','table');
        %deathresult = webread('https://proxy.hxlstandard.org/api/data-preview.csv?url=https%3A%2F%2Fraw.githubusercontent.com%2FCSSEGISandData%2FCOVID-19%2Fmaster%2Fcsse_covid_19_data%2Fcsse_covid_19_time_series%2Ftime_series_19-covid-Deaths.csv&filename=time_series_2019-ncov-Deaths.csv','options','table');
        % updated result
        fprintf('\n\n%1$080d\n%2$46s\n%3$55s\n%1$080d\n',0,'UPDATE WORLD DATA',sprintf('Pause during %0.3g s (CTRL+C to cancel)',pauseduration)), pause(pauseduration)
        if ~offline && (~exist(resultfile,'file') || ~exist(deathresultfile,'file'))
            result  = webread(resultURL);
            deathresult = webread(deathURL);
            if istable(result)
                writetable(result,resultfile,'WriteVariableNames',false);           %  global cases
            else
                fid=fopen(resultfile,'w'); fprintf(fid,'%s',result); fclose(fid);
            end
            if istable(deathresult)
                writetable(deathresult,deathresultfile,'WriteVariableNames',false); %  global death toll
            else
                fid=fopen(deathresultfile,'w'); fprintf(fid,'%s',deathresult); fclose(fid);
            end
        end
        opts = detectImportOptions(resultfile, "TextType","string");
        opts1 = detectImportOptions(deathresultfile, "TextType","string");
        opts.VariableNamesLine = 1; opts.DataLines = [2,inf]; opts.PreserveVariableNames = true;
        opts1.VariableNamesLine = 1; opts1.DataLines = [2,inf]; opts1.PreserveVariableNames = true;
        times_conf = readtable(resultfile,opts);
        times_conf1 = readtable(deathresultfile,opts1);
        
        %
        % +---------------------------------------------------------------------------------------+
        % | FIX Country names for better identification                                           |
        % |                                                                                       |
        % |     ==> unchanged rules since 2020-03-20                                              |
        % |                                                                                       |
        % +---------------------------------------------------------------------------------------+
        %     times_conf.("Country/Region")(times_conf.("Country/Region") == "China") = "Mainland China";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Czechia") = "Czech Republic";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Iran (Islamic Republic of)") = "Iran";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Republic of Korea") = "Korea, South";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Republic of Moldova") = "Moldova";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Russian Federation") = "Russia";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Taipei and environs") = "Taiwan";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Taiwan*") = "Taiwan";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "United Kingdom") = "UK";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "Viet Nam") = "Vietnam";
        times_conf.("Country/Region")(times_conf.("Country/Region") == "occupied Palestinian territory") = "Palestine";
        times_conf.("Country/Region")(times_conf.("Province/State") == "St Martin") = "France";
        times_conf.("Country/Region")(times_conf.("Province/State") == "Saint Barthelemy") = "France";
        % ---------------- remove DOM/TOM (OV), fix 2020-03-21, remove Repatriated Travellers,Canada (fix 2020-11-28)
        times_conf.("Country/Region")(...
            (times_conf.("Country/Region") == "France") & ~ismissing(times_conf.("Province/State")) & (times_conf.("Province/State") ~= "France") ...
            )="France DOM-TOM";
        times_conf(times_conf.("Province/State")=="Repatriated Travellers",:) = [];
        % ---------------- end OV
        %     times_conf1.("Country/Region")(times_conf1.("Country/Region") == "China") = "Mainland China";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Czechia") = "Czech Republic";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Iran (Islamic Republic of)") = "Iran";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Republic of Korea") = "Korea, South";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Republic of Moldova") = "Moldova";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Russian Federation") = "Russia";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Taipei and environs") = "Taiwan";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Taiwan*") = "Taiwan";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "United Kingdom") = "UK";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "Viet Nam") = "Vietnam";
        times_conf1.("Country/Region")(times_conf1.("Country/Region") == "occupied Palestinian territory") = "Palestine";
        times_conf1.("Country/Region")(times_conf1.("Province/State") == "St Martin") = "France";
        times_conf1.("Country/Region")(times_conf1.("Province/State") == "Saint Barthelemy") = "France";
        % ---------------- remove DOM/TOM (OV), , aggregate China
        times_conf1.("Country/Region")(...
            (times_conf1.("Country/Region") == "France") & ~ismissing(times_conf1.("Province/State")) & (times_conf1.("Province/State") ~= "France") ...
              )="France DOM-TOM";
        times_conf1(times_conf1.("Province/State")=="Repatriated Travellers",:) = [];
        % ---------------- end OV
        
        % fix for Matlab R2020a
        times_conf.("Province/State")(ismissing(times_conf.("Province/State"))) = "";
        times_conf1.("Province/State")(ismissing(times_conf1.("Province/State"))) = "";
        % end fix
        
        % +-----------------------------------------------------------------------------------------+
        % | GLOBAL DATA all countries and regions except US states (specific aggregation procedure) |
        % | and detailed French data (regions and departments) which are coming from a different    |
        % | project.                                                                                |
        % |                                                                                         |
        % |     ==> The many changes in the management reflect the expansion of COVID itself        |
        % |                                                                                         |
        % +-----------------------------------------------------------------------------------------+
        % load all data and update with the last available day
        lastGLOBALrawresult = loadGLOBALresults(yesterday);
        % default object
        defautGLOBAL  =  times_conf(find(times_conf.("Country/Region")=="China",1,'first'),:); 
        defautGLOBAL1 = times_conf1(find(times_conf.("Country/Region")=="China",1,'first'),:); 
        defautGLOBAL{1,1} =""; defautGLOBAL{1,2} ="to be set"; defautGLOBAL{1,3:end} =0;
        defautGLOBAL1{1,1}=""; defautGLOBAL1{1,2}="to be set"; defautGLOBAL1{1,3:end}=0;
        % list of Province/State     Country/Region
        tmp = arrayfun(@(r) r.data(:,["Province/State" "Country/Region" "lat" "lon"]),lastGLOBALrawresult,'UniformOutput',false);
        tmp = cat(1,tmp{:});
        tmp(ismissing(tmp(:,2)),:)=[]; % removing missing entries
        % entities
        E= fliplr(unique(tmp(:,[2 1])));
        nE = height(E);
        % times_conf
        hT = height(times_conf);
        TC = times_conf.("Country/Region");
        TS = times_conf.("Province/State"); TS(TS=="")="none";
        Tc = lower(TC);
        Ts = lower(TS);
        % times_conf1
        hT1 = height(times_conf1);
        TC1 = times_conf1.("Country/Region");
        TS1 = times_conf.("Province/State"); TS1(TS1=="")="none";
        Tc1 = lower(TC1);
        Ts1 = lower(TS1);
        % missing coordinates
        missingE = table(... 
            ["North Ireland" ; "Montserrat" ; "none"     ; "Western Australia"; "External territories" ;"Jervis Bay Territory" ; "none"       ;"none"  ;"Falkland Islands (Islas Malvinas)"],...
            ["UK"            ;  "UK"        ; "Australia"; "Australia"        ; "Australia"            ;"Australia"            ; "Madagascar" ;"Malawi";"UK"                               ],...
            [ 53.1424        ; 16.742498    ; -26.4391   ; -25.760321         ; -10.4912               ; -35.1413              ; -18.7793     ;-13.2512;-51.563412                         ],... 
            [ 7.6921         ; -62.187366   ; 133.2813   ; 122.805176         ; 105.6230               ; 150.6913              ; 46.8345      ;34.3015 ;-59.820557                         ],...
            'VariableNames',{'Province/State' 'Country/Region' 'lat' 'lon'});
        % Building correspondence list:: listentities --> times_conf, times_conf1
        list2conf = NaN(nE,4); % [ row # times_conf ][ row # times_conf1 ][ lat ][ lon ]
        for ie = 1:nE
            je = find(ismember(tmp(:,[1,2]),E(ie,:)) & ~isnan(tmp{:,3}) & (tmp{:,3}~=0),1,'last');
            if isempty(je)
                fprintf('\t-->no Lat,Lon for %s/%s\n',C,S);
            else
                list2conf(ie,[3 4]) = tmp{je,[3 4]};
            end
            C = E.("Country/Region")(ie); c = lower(C);
            S = E.("Province/State")(ie); s = lower(S);
            ic  = find( (Tc==c) & (Ts==s) );
            ic1 = find( (Tc1==c) & (Ts1==s) );
            nic=max(length(ic),length(ic1));
            nocoord = isnan(list2conf(ie,3));
            if nic==0
                ke = find(ismember(missingE(:,[1 2]),E(ie,:)));
                % for the special case of an Uknown region in a Country, the Lat,Long values of the country apply
                % the table missingE is updated dynamically (no more need to add the entity by hand) - added 2020-05-16, OV
                % 'quarantine' added on May 29, 2020 dur Port Quarantine in Japan
                if nocoord && isempty(ke) && ( (S=="Unknown") || ~isempty(regexpi(S,'quarantine','once')) )
                    fprintf('\t-->add dynamically the missing entry (if it is not a success, read the following messages) %s/%s\n',C,S);
                    ic  = find( (Tc==c) & ( (Ts=="none") | (Ts=="beijing") ) ); % for China, beijing is used
                    if ~isempty(ic)
                        missingE(end+1,:) = table(S,C,times_conf.Lat(ic(1)),times_conf.Long(ic(1)),'VariableNames',{'Province/State' 'Country/Region' 'lat' 'lon'}); %#ok<AGROW>
                        ke = find(ismember(missingE(:,[1 2]),E(ie,:)));
                    end
                end               
                if nocoord && isempty(ke)
                    error('\t-->no Lat/Lon data for the missing entry %s/%s\n',C,S);
                else
                    fprintf('\t-->create an entry for the missing entry %s/%s\n',C,S);
                    hT = hT + 1;
                    hT1 = hT1 + 1;
                    times_conf(hT,:) = defautGLOBAL;
                    times_conf1(hT1,:) = defautGLOBAL1;
                    times_conf{hT,2} =C;  times_conf{hT,1} =regexprep(S,"^none",""); 
                    times_conf1{hT1,2}=C; times_conf1{hT1,1}=regexprep(S,"^none","");
                    list2conf(ie,1:2) = [hT hT1];
                    if nocoord, list2conf(ie,[3 4])  = missingE{ke,[3 4]}; end
                    times_conf{hT,3:4}   = list2conf(ie,[3 4]);
                    times_conf1{hT1,3:4} = list2conf(ie,[3 4]);
                end
            elseif nic==1
                if nocoord
                    ke = find(ismember(missingE(:,[1 2]),E(ie,:)));
                    if isempty(ke), error('\t-->no Lat/Lon data for the existing entry %s/%s, please add them\n',C,S); end
                    fprintf('\t-->add Lat,Lon to the existing entry %s/%s\n',C,S);
                    list2conf(ie,[3 4])  = missingE{ke,[3 4]};
                    tobefixed = ismember(times_conf(:,1:2),E(ie,:));
                    times_conf{tobefixed,3} = list2conf(ie,3);
                    times_conf{tobefixed,4} = list2conf(ie,4);
                    tobefixed = ismember(times_conf1(:,1:2),E(ie,:));
                    times_conf1{tobefixed,3} = list2conf(ie,3);
                    times_conf1{tobefixed,4} = list2conf(ie,4);
                end
                list2conf(ie,1:2) = [ic ic1];
            else
                for i=1:nic
                    fprintf('\t-->%d entries for %s/%s :\n%s',nic,C,S,fprintf('\t\t%s/%s\n',TC(ic(i)),TS(ic(i))));
                end
                error('Please check the results above, they are not consistent')
            end
        end % next ientity
        
        
        % listofdates
        [~,ind] = sort(cat(1,lastGLOBALrawresult.date),'descend');
        lastGLOBALrawresult = lastGLOBALrawresult(ind); nlast = length(lastGLOBALrawresult);
        tolerance = 0.2;
        for ilast = 1:nlast
            d = lastGLOBALrawresult(ilast).datestr;
            if ~ismember(d,times_conf.Properties.VariableNames), fprintf('WARNING:: %s is newer than the one in the database\n',d); end
            tmp = lastGLOBALrawresult(ilast).data(:,1:2);
            tmp.("Province/State")(tmp.("Province/State")=="none")="";
            tmp(ismissing(tmp(:,2)),:)=[]; % remove missing data on 2020-11-05
            [E,ie,je] = intersect(tmp,times_conf(:,1:2),'stable');
            nE = height(E);
            if nE<height(tmp), error('global database does not fit updated data, please fix it'), end
            for i = 1:nE
                if isempty(E.("Province/State"){i})
                    Etxt = E.("Country/Region"){i};
                else
                    Etxt = sprintf('%s/%s',E.("Country/Region"){i},E.("Province/State"){i});
                end
                % times_conf
                currentvalue = times_conf.(d)(je(i));
                updatedvalue = lastGLOBALrawresult(ilast).data.cases(ie(i));
                if updatedvalue>currentvalue, warn='INCR CASES'; else, warn='DECR CASES'; end %#ok<NASGU>
                if (isnan(currentvalue) || currentvalue==0) && (updatedvalue>0)
                    fprintf('[%s:%s] assign %0.6g CASES to (no cases before)\n',d,Etxt,updatedvalue);
                    updt = true;
                elseif (abs(updatedvalue-currentvalue)/currentvalue) > 15*tolerance
                    if updatedvalue>currentvalue, warn='INCR Confirmed CASES'; else, warn='DECR Confirmed CASES'; end
                    warning('[%s:%s] %s:: variation of %d CASES (%d->%d) (ERROR, far above tolerance %0.3g %%)',...
                        d,Etxt,warn,updatedvalue-currentvalue,currentvalue,updatedvalue,100*tolerance);
                    updt = true;
                elseif (abs(updatedvalue-currentvalue)/currentvalue) > tolerance
                    if updatedvalue>currentvalue, warn='INCR Death toll'; else, warn='DECR Death toll'; end
                    fprintf('[%s:%s] %s:: variation of %d CASES (%d->%d) (above tolerance %0.3g %%)\n',...
                        d,Etxt,warn,updatedvalue-currentvalue,currentvalue,updatedvalue,100*tolerance);
                    updt = true;
                else
                    updt = false;
                end
                if updt, times_conf.(d)(je(i)) = max(currentvalue,updatedvalue); end
                % times_conf1
                currentvalue = times_conf1.(d)(je(i));
                updatedvalue = lastGLOBALrawresult(ilast).data.deaths(ie(i));
                if updatedvalue>currentvalue, warn='INCR Death toll'; else, warn='DECR Death toll'; end %#ok<NASGU>
                if (isnan(currentvalue) || currentvalue==0) && (updatedvalue>0)
                    fprintf('[%s:%s] assign %0.6g DEATHS to (no cases before)\n',d,Etxt,updatedvalue);
                    updt = true;
                elseif (abs(updatedvalue-currentvalue)/currentvalue) > 15*tolerance
                    if updatedvalue>currentvalue, warn='INCR Death toll'; else, warn='DECR Death toll'; end
                    warning('[%s:%s] %s:: variation of %d DEATHS (%d->%d) (ERROR, far above tolerance %0.3g %%)',...
                        d,Etxt,warn,updatedvalue-currentvalue,currentvalue,updatedvalue,100*tolerance);
                    updt = true;
                elseif (abs(updatedvalue-currentvalue)/currentvalue) > tolerance
                    if updatedvalue>currentvalue, warn='INCR Death toll'; else, warn='DECR Death toll'; end
                    fprintf('[%s:%s] %s:: variation of %d DEATHS (%d->%d) (above tolerance %0.3g %%)\n',...
                        d,Etxt,warn,updatedvalue-currentvalue,currentvalue,updatedvalue,100*tolerance);
                    updt = true;
                else
                    updt = false;
                end
                if updt, times_conf1.(d)(je(i)) = max(currentvalue,updatedvalue); end
            end % next i
        end % next ilast
        
        
        % +--------------------------------------------------------------------------------------+
        % | DETAILED DATA FOR METROPOLITAN FRANCE (DOM-TOM managed through WHO, more consistent) |
        % |                                                                                      |
        % |     ==> two levels of granularity are considered: departments (96) and regions (13)  |
        % |                                                                                      |
        % +--------------------------------------------------------------------------------------+
        FRrawresults = loadFRANCEresults();
        % consolidation of WHOS data (a few mistakes) by national ones
        iFR = find(times_conf.("Country/Region")=="France");
        iFR1 = find(times_conf1.("Country/Region")=="France");        
        FRthresh = 0.3;
        [current_who,current_fr] = deal(struct('deaths',0,'cases',0));
        if any(diff(times_conf{iFR,5:end})), warning('Inconsistent WHOS Cases data for FFRANCE (%d consecutive days with lowering numbers)',length(find(diff(times_conf{iFR,5:end})))); end
        times_conf{iFR1,5:end} = cummax(times_conf{iFR1,5:end});
        
        if any(diff(times_conf1{iFR1,5:end})), warning('Inconsistent WHOS Deaths data for FFRANCE (%d consecutive days with lowering numbers)',length(find(diff(times_conf{iFR,5:end})))); end
        times_conf1{iFR1,5:end} = cummax(times_conf1{iFR1,5:end});

        for i = 1:length(FRrawresults.national)
            d = FRrawresults.national(i).datestr;
            if ismember(d,times_conf.Properties.VariableNames)
                current_who = struct('deaths',max(current_who.deaths,times_conf1{iFR1,d}),...
                    'cases' ,max(current_who.cases, times_conf{iFR,d}));
                current_fr  = struct('deaths',max(current_fr.deaths,FRrawresults.national(i).data.deaths),...
                    'cases',max(current_fr.cases,  FRrawresults.national(i).data.cases));
                if current_fr.deaths > current_who.deaths
                    if (current_fr.deaths>10) && (current_fr.deaths/current_who.deaths - 1)>FRthresh
                        warning('[FRANCE %s - DEATH TOLL check]\t threshold (%d %%) exceeded:: WHOS data (%d) vs FRANCE (%d)',d,100*FRthresh,current_who.deaths,current_fr.deaths);
                    else
                        fprintf(1,'[FRANCE %s - DEATH TOLL check]\t WHOS data (%d) underestimate deaths (%d)\n',d,current_who.deaths,current_fr.deaths);
                    end
                    times_conf1{iFR1,FRrawresults.national(i).datestr} = current_fr.deaths;
                    current_who.deaths = current_fr.deaths;
                end
                if current_fr.cases>current_who.cases
                    if (current_fr.cases>100) && (current_fr.cases/current_who.cases - 1)>FRthresh
                        warning('[FRANCE %s - CASES TOLL check]\t threshold (%d %%) exceeded:: WHOS data (%d) vs FRANCE (%d)',d,100*FRthresh,current_who.cases,current_fr.cases);
                    else
                        fprintf(1,'[FRANCE %s - CASES TOLL check]\t WHOS data (%d) underestimate cases (%d)\n',d,current_who.cases,current_fr.cases);
                    end
                    times_conf{iFR,FRrawresults.national(i).datestr} = current_fr.cases;
                    current_who.cases = current_fr.cases;
                end
            end % date exist (difference between France and Who submissions)
        end % next i
            
        % default content
        defautFR = times_conf(iFR,:); defautFR{1,5:end}=0;
        defautFR1 = times_conf1(iFR1,:); defautFR1{1,5:end}=0;
        defautFR.("Country/Region")='France/Regions/';
        defautFR1.("Country/Region")='France/Regions/';
        % French regions
        lastFRrawresult = FRrawresults.region;
        listofstates = lastFRrawresult(end).data.state; % same construction of US
        listofdates = {lastFRrawresult.datestr};
        [~,ilistofdates,jlistofdates] = intersect(defautFR.Properties.VariableNames,listofdates);
        ietat = 0;
        for etat = listofstates'
            ietat = ietat+1;
            times_conf(end+1,:) = defautFR; %#ok<AGROW>
            times_conf.("Province/State")(end) = "Région "+string(etat);
            times_conf.("Lat")(end) = lastFRrawresult(end).data.lat(ietat);
            times_conf.("Long")(end) = lastFRrawresult(end).data.lon(ietat);
            times_conf1(end+1,:) = defautFR1; %#ok<AGROW>
            times_conf1.("Province/State")(end) = "Région "+string(etat);
            times_conf1.("Lat")(end) = lastFRrawresult(end).data.lat(ietat);
            times_conf1.("Long")(end) = lastFRrawresult(end).data.lon(ietat);
            for idate = 1:length(ilistofdates) %#ok<FXUP>
                jetat = ismember(lastFRrawresult(jlistofdates(idate)).data.state,etat);
                if ismember(etat,lastFRrawresult(jlistofdates(idate)).data.state)
                    times_conf{end,ilistofdates(idate)} = lastFRrawresult(jlistofdates(idate)).data.cases(jetat);
                    times_conf1{end,ilistofdates(idate)} = lastFRrawresult(jlistofdates(idate)).data.deaths(jetat);
                end
            end
        end
        % French departments
        defautFR.("Country/Region")='France/Departements/';
        defautFR1.("Country/Region")='France/Departements/';
        lastFRrawresult = FRrawresults.departement;
        listofstates = lastFRrawresult(end).data.state; % same construction of US
        listofdates = {lastFRrawresult.datestr};
        [~,ilistofdates,jlistofdates] = intersect(defautFR.Properties.VariableNames,listofdates);
        ietat = 0;
        for etat = listofstates'
            ietat = ietat+1;
            times_conf(end+1,:) = defautFR; %#ok<AGROW>
            times_conf.("Province/State")(end) = "DPT "+string(etat);
            times_conf.("Lat")(end) = lastFRrawresult(end).data.lat(ietat);
            times_conf.("Long")(end) = lastFRrawresult(end).data.lon(ietat);
            times_conf1(end+1,:) = defautFR1; %#ok<AGROW>
            times_conf1.("Province/State")(end) = "DPT "+string(etat);
            times_conf1.("Lat")(end) = lastFRrawresult(end).data.lat(ietat);
            times_conf1.("Long")(end) = lastFRrawresult(end).data.lon(ietat);
            for idate = 1:length(ilistofdates) %#ok<FXUP>
                jetat = ismember(lastFRrawresult(jlistofdates(idate)).data.state,etat);
                if ismember(etat,lastFRrawresult(jlistofdates(idate)).data.state)
                    times_conf{end,ilistofdates(idate)} = lastFRrawresult(jlistofdates(idate)).data.cases(jetat);
                    times_conf1{end,ilistofdates(idate)} = lastFRrawresult(jlistofdates(idate)).data.deaths(jetat);
                end
            end
        end
        %
        % +---------------------------------------------------------------------------------------+
        % | USA data, aggregated by state (workaround, the procedure has been changed 2020-03-22) |
        % |                                                                                       |
        % |     ==> data are aggregated by state (same source as GLOBAL data)                     |
        % |                                                                                       |
        % +---------------------------------------------------------------------------------------+
        lastUSrawresult = loadUSresults(yesterday);
        defautUS = times_conf(times_conf.("Country/Region")=="US",:); defautUS{1,5:end}=0;
        defautUS1 = times_conf1(times_conf.("Country/Region")=="US",:); defautUS1{1,5:end}=0;
        listofstates = lastUSrawresult(end).data.state;
        listofdates = {lastUSrawresult.datestr};
        [~,ilistofdates,jlistofdates] = intersect(defautUS.Properties.VariableNames,listofdates);
        ietat = 0;
        for etat = listofstates'
            ietat = ietat+1;
            times_conf(end+1,:) = defautUS; %#ok<AGROW>
            times_conf.("Province/State")(end) = string(etat);
            times_conf.("Lat")(end) = lastUSrawresult(end).data.lat(ietat);
            times_conf.("Long")(end) = lastUSrawresult(end).data.lon(ietat);
            times_conf1(end+1,:) = defautUS1; %#ok<AGROW>
            times_conf1.("Province/State")(end) = string(etat);
            times_conf1.("Lat")(end) = lastUSrawresult(end).data.lat(ietat);
            times_conf1.("Long")(end) = lastUSrawresult(end).data.lon(ietat);
            for idate = 1:length(listofdates) %#ok<FXUP>
                jetat = ismember(lastUSrawresult(jlistofdates(idate)).data.state,etat);
                if ismember(etat,lastUSrawresult(jlistofdates(idate)).data.state)
                    times_conf{end,ilistofdates(idate)} = lastUSrawresult(jlistofdates(idate)).data.cases(jetat);
                    times_conf1{end,ilistofdates(idate)} = lastUSrawresult(jlistofdates(idate)).data.deaths(jetat);
                end
            end
        end

        %
        % +----------------------------------------------------------------------------------------------------------+
        % | Create "pseudo" entities for states and provinces, e.g. US,Connecticut (added 2020-03-23)                |
        % | Prepare map data (added 2020-03-25), Note that mapresult2 contains time resolved data (added 2020-04-02) |                                                                                      |
        % |                                                                                                          |
        % |     ==> to be updated according to the evolution of the outbreak                                         |                                          |
        % +----------------------------------------------------------------------------------------------------------+
        %before April 25th: death scale needs to be updated manually
        %deathbounds = [0 10 50 100 250 500 1000 5000 10000 20000 30000 40000 50000 70000];
        %deathcategories = arrayfun(@(x) sprintf('less than %d',x),deathbounds(2:end),'UniformOutput',false);
        ndeathcategories_ = 15;
        deathbounds_ = [0 10:10:40  50:50:250 300:100:800 1000:500:3000 4000:1000:10000 15000:5000:30000 40000:1e4:200000 25e4:5e4:1e6];
        deathdecade_ = @(dmax) 10^floor(log10(dmax));
        deathnext_   = @(dmax) deathdecade_(dmax)*ceil(dmax/deathdecade_(dmax));
        removefirst_ = @(x) x(2:end);
        deathbounds = @(dmax) deathbounds_(interp1(deathbounds_,1:length(deathbounds_),linspace(0,sqrt(deathnext_(dmax)),ndeathcategories_).^2,'nearest'));
        deathcategories = @(dmax) arrayfun(@(x) sprintf('less than %d',x),removefirst_(deathbounds(dmax)),'UniformOutput',false);
        % nentries x ndates (note that the title of the map of today shows now the date of the results)
        listofdates = cellstr(datestr(datenum(times_conf.Properties.VariableNames(5:end),'mm/dd/yy'),'yyyy-mm-dd'));
        ndates = length(listofdates);
        nentries = height(times_conf);
        mapresult = cell2table(repmat({"" NaN NaN NaN NaN NaN},nentries,1),...
            'VariableNames',{'country','lat','lon','cases','deathtoll','deathcategory'});
        mapresult2 = cell2table(repmat({"" NaN NaN zeros(1,ndates) zeros(1,ndates) NaN(1,ndates)},nentries,1),...
            'VariableNames',{'country','lat','lon','cases','deathtoll','deathcategory'});
        % for each entry
        firstdate = zeros(nentries,1)+ndates;
        for i=1:nentries
            pays = times_conf.('Country/Region')(i);
            etat = times_conf.('Province/State')(i);
            if ~ismissing(etat) && etat ~="" && etat~="The"
                if any(pays{1}=='/')
                    pseudo = pays+etat;
                else
                    pseudo = regexprep(sprintf('%s,%s',pays,etat),{'\,\s+' ''''},{',',''});
                end
                fprintf('--> Create %s\n',pseudo)
                times_conf(end+1,:) = times_conf(i,:); %#ok<AGROW>
                times_conf1(end+1,:) = times_conf1(i,:); %#ok<AGROW>
                times_conf.('Country/Region')(end)=pseudo;
                times_conf1.('Country/Region')(end)=pseudo;
                mapresult.country(i) = pseudo;
                mapresult2.country(i) = pseudo;                
            else
                mapresult.country(i) = pays;
                mapresult2.country(i) = pays;
            end
            mapresult{i,2:3} = [times_conf.Lat(i) times_conf.Long(i)];
            mapresult{i,4} = sum(max(cummax(times_conf{i,5:end},2),[],2),1); %sum(times_conf{i,end});   % fix 2020/03/30
            mapresult{i,5} = sum(max(cummax(times_conf1{i,5:end},2),[],2),1); %sum(times_conf1{i,end}); % fix 2020/03/30
            mapresult2{i,2:3} = [times_conf.Lat(i) times_conf.Long(i)];            
            mapresult2{i,4} = max(0,sum(cummax(times_conf{i,5:end},2),1));
            mapresult2{i,5} = max(0,sum(cummax(times_conf1{i,5:end},2),1));
            itmp = find(mapresult2{i,4}>0,1,'first');
            if ~isempty(itmp), firstdate(i) = itmp; end
        end
        if any(isnan(mapresult.lat) | isnan(mapresult.lon)), error('some geographical data are missing, please check'), end
        mapresult.cases = max(0,mapresult.cases);
        dmax_ = max(mapresult.deathtoll);
        mapresult.deathcategory = discretize(max(0,mapresult.deathtoll),deathbounds(dmax_),'categorical',deathcategories(dmax_));
        %mapresult = mapresult(mapresult.cases>0,:); % 0 data are removed (added 2020-03-26)
        %mapresult((mapresult.lat==0) & (mapresult.lon==0),:) = [];
        badmaps = (mapresult.cases==0) | ((mapresult.lat==0) & (mapresult.lon==0));
        mapresult(badmaps,:) = [];
        mapresult2(badmaps,:) = [];
        nentries = height(mapresult2);

        % ##############################################################
        % ##############################################################
        % ##                 R E N D E R    M A P S                   ##
        % ##############################################################
        % ##############################################################

        % debug
        %idebug = find(~cellfun(@isempty,regexp(mapresult.country,'^France,DPT')));
        %mapresult = mapresult(idebug,:);
        %% Define maps (added 2020-03-25)
        maptype = struct(...
            'World',struct('Lat',[-23.8426   76.1121],'Lon',[-166.5742  192.4487],...
            'PaperPosition',[-0.7442    7.1413   22.4725   15.4399],...
            'PaperOrientation','Portrait','dpi',240,'title','World Map','MapLayout','normal',...
            'excluded','France[/\,]','BubbleWidthRange',[3 30]),...
            'Europe',struct('Lat',[32.6780   64.6604],'Lon',[-14.4820   49.5634],'title','Europe',...
            'PaperPosition',[0.6345    0.6505   28.4084   19.7150],...
            'PaperOrientation','Landscape','dpi',200,'MapLayout','maximized',...
            'excluded','France[/\,]','BubbleWidthRange',[7 70]));
        maptype.America = maptype.Europe;
        maptype.America.Lat=[-43.8426   76.1121];
        maptype.America.Lon=[-166.5742  -30];
        maptype.America.title = 'America';
        maptype.Asia = maptype.Europe;
        maptype.Asia.Lat=[-43.8426   76.1121];
        maptype.Asia.Lon=[30.5742  166];
        maptype.Asia.title = 'Asia';
        maptype.Asia.BubbleWidthRange = [7 40];
        maptype.Asia.excluded='^US|^France[/\,]';
        maptype.France = maptype.Europe;
        maptype.France.Lat=[41.9877   51.4434];
        maptype.France.Lon=[-5.1139   13.3431];
        maptype.France.title = 'France';
        maptype.France.excluded='^(?!France).*$'; % does not start with France (see: https://superuser.com/questions/655715/regex-does-not-begin-by-pattern)
        if ~isempty(setdiff(listofmaps,fieldnames(maptype))), error('several requested maps are not defined, please define them first'), end
            
        % PLOT MAPS
        if mapon && ~isempty(listofmaps)            
            maptitle = @(where,when) sprintf('COVID-19 - %s on %s',where,when);
            hfigmap = figure;
            mapresult = flipud(mapresult); % entities are plotted last
            mapresult2 = flipud(mapresult2); % entities are plotted last
            cases =  mapresult.cases;
            
            for mtp=listofmaps(:)'
                % figure defintion
                figname = sprintf('COVID_MAP_%s_%s',mtp{1},datetxt);
                GIFfile = fullfile(mappath,sprintf('%s.gif',figname));
                MP4file = fullfile(mappath,sprintf('%s.mp4',figname));   %    H265 video
                WEBMfile = fullfile(mappath,sprintf('%s.webm',figname)); %    VP9  video
                if GIFon && exist(GIFfile,'file'), delete(GIFfile), end % to be cleaned before adding new frames
                % figure layout
                included = cellfun(@isempty,regexp(mapresult.country,maptype.(mtp{1}).excluded));
                if ismember(mtp{1},{'France' 'Europe'})
                    included = included & ...
                        (mapresult.lat>=maptype.(mtp{1}).Lat(1)) & (mapresult.lat<=maptype.(mtp{1}).Lat(2)) & ...
                        (mapresult.lon>=maptype.(mtp{1}).Lon(1)) & (mapresult.lon<=maptype.(mtp{1}).Lon(2));
                end
                if strcmpi(mtp{1},'France')
                    included = included & ~cellfun(@isempty,regexp(mapresult.country,'^France/'));
                end
                dmax_ = max(mapresult.deathtoll(included));
                col = repmat([0.8 0.8 0.8],length(deathbounds(dmax_))-1,1);
                ncolvalid = find(deathbounds(dmax_)<max(mapresult.deathtoll(included)),1,'last'); % old code kept even if ncolvalid = ndeathcategories_
                col(1:ncolvalid,:)=flipud(tooclear(cbrewer('div','RdYlGn',ncolvalid)));
                % map for all dates (GIF), last state (static PNG) --> static PNG to be used for video Poster
                if GIFon, date2map=min(firstdate(included)):ndates; else, date2map=ndates; end %#ok<UNRCH>
                maxcases = max(cases(included));
                for idate=date2map %#ok<FXUP>
                    clf(hfigmap)
                    % populate mapresult (static) from mapresult2 (with time resolution)
                    mapresult{:,4:5} = zeros(nentries,2);
                    for ientry=1:nentries
                        mapresult{ientry,4:5} = [mapresult2{ientry,4}(idate) mapresult2{ientry,5}(idate)];
                    end % next entry
                    current_included = included & (mapresult{:,4}>0);
                    mapresult.deathcategory = discretize(max(0,mapresult.deathtoll),deathbounds(dmax_),'categorical',deathcategories(dmax_));
                    % do the plot
                    hg=geobubble(mapresult(current_included,:),'lat','lon',...
                        'SizeVariable','cases',...
                        'SizeLimits',[1 maxcases],...
                        'ColorVariable','deathcategory',...
                        'MapLayout',maptype.(mtp{1}).MapLayout,...
                        'BaseMap','grayterrain',...
                        'Title',maptitle(maptype.(mtp{1}).title,listofdates{idate}),...
                        'BubbleColorList', col,...
                        'BubbleWidthRange',maptype.(mtp{1}).BubbleWidthRange,...
                        'ColorLegendTitle','Death Toll',...
                        'SizeLegendTitle','Confirmed Cases'); %#ok<NASGU>
                    % rename figure and apply attributes
                    formatfig(hfigmap,...
                        'figname',figname,...
                        'paperposition',maptype.(mtp{1}).PaperPosition,...
                        'paperorientation',maptype.(mtp{1}).PaperOrientation);
                    % change view
                    geolimits(maptype.(mtp{1}).Lat,maptype.(mtp{1}).Lon)
                    % add title if maximized
                    if strcmpi(maptype.(mtp{1}).MapLayout,'maximized')
                        drawnow
                        hax = axes('position',[ 0.1300    0.1100    0.7750    0.8150 ],'visible','off');
                        ht= titles(hax,maptitle(maptype.(mtp{1}).title,listofdates{idate}),'x',.2,'y',-.06,'fontsize',10,'color',rgb('Beige'));
                        delete(hax), set(ht.text,'visible','on')
                    end
                    drawnow
                    % grab the frame if GIFon
                    if GIFon
                        fprintf('... grab %s on %s\n',maptype.(mtp{1}).title,listofdates{idate});
                        gif_add_frame(hfigmap,GIFfile,5);
                    end
                end % next idate
                % print maps
                print_png(maptype.(mtp{1}).dpi,figname,mappath,'',false,0,0,false)
                copyfile(fullfile(mappath,sprintf('%s.png',figname)),fullfile(mappath,sprintf('LAST_update_%s.png',mtp{1})))
                % Videos and naming rules for
                if GIFon % generate h265, VP9 movies
                    % intermediate files include date; they are renamed at the end of the process to enable static links and save space
                    animationfile = @(ext) fullfile(mappath,sprintf('LAST_animation_%s.%s',mtp{1},ext));
                    % MP4 (h265) video (old procedure: system(sprintf('ffmpeg -f gif -i ''%s'' ''%s''',GIFfile,MP4file));)
                    if exist(MP4file,'file'), delete(MP4file), end, gif2mp4(GIFfile,MP4file);     % -> for IOS
                    % webm (VP9) video (old procedure: system(sprintf('ffmpeg -f gif -i ''%s'' -vcodec libvpx ''%s''',GIFfile,WEBMfile)); )
                    if exist(WEBMfile,'file'), delete(WEBMfile), end, gif2webm(GIFfile,WEBMfile); % -> for Chrome
                    % move all intermediate files
                    if exist(animationfile('mp4'),'file'), delete(animationfile('mp4')), end, movefile(MP4file,animationfile('mp4'),'f');
                    fprintf('\n---> MP4 (h265) animation generated\n\t for ''%s''\n',maptitle(maptype.(mtp{1}).title,listofdates{idate}));
                    fileinfo(animationfile('mp4'))
                    if exist(animationfile('webm'),'file'), delete(animationfile('webm')), end, movefile(WEBMfile,animationfile('webm'),'f');
                    fprintf('\n---> WEBM (VP9) animation generated\n\t for ''%s''\n',maptitle(maptype.(mtp{1}).title,listofdates{idate}));
                    fileinfo(animationfile('webm'))
                    if exist(animationfile('gif'),'file'), delete(animationfile('gif')), end, movefile(GIFfile,animationfile('gif'))
                    fprintf('\n---> GIF animation generated\n\t for ''%s''\n',maptitle(maptype.(mtp{1}).title,listofdates{idate}));
                    fileinfo(animationfile('gif'))
                end
            end % next map
        else
            maptype = [];
        end 
        
        % end pseudo country
        listentities = unique(times_conf.("Country/Region"));
        save(globalprefetchfile,'times_conf','times_conf1','listentities','mapresult','maptype'); 
                % history: 'BubbleColorList'
        %         [
        %                      0    0.4078    0.2157
        %                 0.1020    0.5961    0.3137
        %                 0.4000    0.7412    0.3882
        %                 0.5333    0.8000    0.4000
        %                 0.5010    0.7010    0.2657
        %                 0.7010    0.7873    0.3951
        %                 0.8500    0.8500    0.5990
        %                 0.8461    0.7284    0.3951
        %                 0.8422    0.5324    0.2304
        %                 0.8304    0.4069    0.1676
        %                 0.9569    0.4275    0.2627
        %                 0.8431    0.1882    0.1529
        %                 0.6471         0    0.1490
        %                 ]
    else %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        load(globalprefetchfile,'times_conf','times_conf1','listentities')
        
    end % end globalprefetch ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        % ##############################################################
        % ##############################################################
        % ##############################################################
        % ####            F  O  R  E  C  A  S  T  I  N  G           ####
        % ##############################################################
        % ##############################################################
        % ##############################################################

    
    % check country (added 2020-03-21)
    checkcountry()
    % parse data
    vars = times_conf.Properties.VariableNames;
    vars1 = times_conf1.Properties.VariableNames;
    times_conf_country = groupsummary(times_conf,"Country/Region",{'sum','mean'},vars(3:end));
    times_conf_country1 = groupsummary(times_conf1,"Country/Region",{'sum','mean'},vars1(3:end));
    vars = times_conf_country.Properties.VariableNames;
    vars = regexprep(vars,"^(sum_)(?=L(a|o))","remove_");
    vars = regexprep(vars,"^(mean_)(?=[0-9])","remove_");
    vars = erase(vars,{'sum_','mean_'});
    times_conf_country.Properties.VariableNames = vars;
    vars1 = times_conf_country1.Properties.VariableNames;
    vars1 = regexprep(vars1,"^(sum_)(?=L(a|o))","remove_");
    vars1 = regexprep(vars1,"^(mean_)(?=[0-9])","remove_");
    vars1 = erase(vars1,{'sum_','mean_'});
    times_conf_country1.Properties.VariableNames = vars;
    infectedtable = removevars(times_conf_country,[{'GroupCount'},vars(contains(vars,"remove_"))]);
    countrytable = infectedtable(strcmp(infectedtable.("Country/Region"),country), :);
    if isempty(countrytable), error('unknown country ''%s''',country); end
    deathtable = removevars(times_conf_country1,[{'GroupCount'},vars1(contains(vars1,"remove_"))]);
    countrytable1 = deathtable(strcmp(deathtable.("Country/Region"),country), :);
    countrytable = countrytable(:,4:end);
    countrytable1 = countrytable1(:,4:end);
    cols1 = size(countrytable);
    cols2 = size(countrytable1);
    Countrytotaldead = zeros(1,cols2(2));
    Countrytotalinfected = zeros(1,cols1(2));
    for i = 1:cols1(2)
        Countrytotalinfected(i) = table2array(countrytable(1,i));
        Countrytotaldead(i) = table2array(countrytable1(1,i));
    end
    % remove trailing 0 dur to not updated results (added 2020/03/30), French data are not updated early in the morning
    lastnonzero = max(find(Countrytotalinfected>0,1,'last'),find(Countrytotaldead>0,1,'last'));
    if ~isempty(lastnonzero) && lastnonzero<cols1(2)
        fprintf('\t [%s] remove %d trailing zero(s) (data not fully up-to-date)\n',country,cols1(2)-lastnonzero);
        Countrytotalinfected = Countrytotalinfected(1:lastnonzero);
        Countrytotaldead = Countrytotaldead(1:lastnonzero);
        lastdaywithdata = 'LAST DAY WITH DATA';
    else
        lastdaywithdata = 'TODAY';
    end
    % end remove
    infected = Countrytotalinfected;
    deathrate = min(deathratemax,Countrytotaldead(end)/Countrytotalinfected(end));
    warning('off')
    save(prefetchfile,listvars{:})
else % use the PREFETCH (OV)
    load(prefetchfile,listvars{:})
end


%algorithm to determine how many slow-growth phase points to include
startidx = zeros(length(infected));
for i = 1:length(infected)
    if i == length(infected)
        break
    end
    if infected(i+1) > 1.5*infected(i) && infected(i) > 0
        startidx(i) = i;
    end
end
startidx = find(startidx~=0, 1, 'first');
first_day=datenum('2020/01/22'); % spread start date - China (do not change)
first_day = first_day+startidx;
fprintf('Virus Spread Prediction for: %s\n',country)
sampledTime = 0:1:length(infected(startidx:end))-1;
sampledDate = first_day + sampledTime;
[b0] = iniGuess(infected(startidx:end));
if isempty(b0) || b0(2) == Inf || b0(3) == Inf
    fprintf('***Warning: Fail to calculate initial quess. Use default.\n');
    b0 = [max(infected) 0.5 max(infected)]';
end
K0 = b0(1);
r0 = b0(2);
A0 = b0(3);
fprintf('  Initial guess K = %g  r = %g  A = %g\n',K0,r0,A0);
infected = infected(startidx:end);
infectedidx = length(infected);

K      = NaN(length(infected),1);
r      = NaN(length(infected),1);
A      = NaN(length(infected),1);
C0     = NaN(length(infected),1);
tpeak  = NaN(length(infected),1);
dpeak  = NaN(length(infected),1);
dend   = NaN(length(infected),1);
dCpeak = NaN(length(infected),1);
tau    = NaN(length(infected),1);
err    = NaN(length(infected),1);
R2     = NaN(length(infected),1);

n0 = ceil(0.25*length(infected));

opts = optimoptions('lsqcurvefit','Display','off','SpecifyObjectiveGradient',true);
for n = n0:length(infected) %#ok<FXUP>
    [b,resnorm,residual,~,~,~,Jacobian] = lsqcurvefit(@fun,b0,sampledTime(1:n),infected(1:n),[0 0 0],[],opts);    
    err(n) = resnorm;
    R2(n) = calcR2(sampledTime(1:n),infected(1:n),b);
    K(n)   = fix(b(1));
    r(n)   = b(2);
    A(n)   = b(3);
    C0(n)  = fix(K(n)/(A(n) + 1));
    tpeak(n)  = fix(log(A(n))/r(n));
    dpeak(n)  = tpeak(n) + first_day;
    dend(n)   = 2*tpeak(n) + first_day;
    dCpeak(n) = fix(r(n)*K(n)/4);
    tau(n)    = fix(4/r(n));
end
if issparse(Jacobian), Jacobian = full(Jacobian); end

fprintf('Regression parameters for the full dataset\n')
mdl = fitnlm(sampledTime(1:n),infected(1:n),@fun,b,...
    'CoefficientNames',{'K','r','A'});
coef = mdl.Coefficients.Estimate;
if abs(coef(1)/K(infectedidx)) > 2
    fprintf('***Warning: results of lsqcurvefit and fitnlm differs significantly. \n');
    fprintf('   Knlm/Klsq = %g\n',coef(1)/K(infectedidx));
    fprintf('   rnlm/rlsq = %g\n',coef(2)/r(infectedidx));
    fprintf('   Anlm/Alsq = %g\n',coef(3)/A(infectedidx));
end

fprintf('\nEvaluation of model parameters for %s\n',country)
fprintf('%4s %10s %8s %8s %7s %5s %5s %10s %5s %5s %10s %5s\n',...
    'day','date','C','K','r','C0','Tau','end','dCdt','tpeak','peak','R2')
fprintf('%4s %10s %8s %8s %7s %5s %5s %10s %5s %5s %10s\n',...
    ' ',' ','(cases)','(cases)','(1/day)','(cases)','(day)',' ','(c/day)','(day)',' ')
for n = n0:length(infected) %#ok<FXUP>
    fprintf('%4d %10s %8d %8d %7.3f %5d %5d %10s %5d %5d %10s %5.3f\n',...
        n,datestr(sampledDate(n)),infected(n),K(n),r(n),C0(n),...
        tau(n),datestr(dend(n)),dCpeak(n),tpeak(n),datestr(dpeak(n)),R2(n));
end

time = 0:1:max(ceil(2.75*tpeak(n)),ceil(now-first_day)); % extended period: 
date = first_day + time;

%% additional check after Palestine data made Matlab entering within an infinite loop on April, 11, 2020
lowest_b = 1/(25*365); % time scale, 25 years
if b(2) < lowest_b
    fprintf(1,'No convergence for %s (too large time constant %0.6g years)\n',country,1/(365*b(2)));
    return
end

%% updated death rate (added OV)
spdead = spaps(Countrytotalinfected,Countrytotaldead,0.01*var(Countrytotaldead)*length(Countrytotaldead));
deathrate2 = fnval(fnder(spdead),Countrytotalinfected);
if deathrate2(end)>0
    deathrate2 = deathrate2(end);
elseif any(deathrate2>0)
    deathrate2 = max(deathrate2);
else
    deathrate2 = deathrate;
end
deathrate2 = min(deathratemax,max(deathrate,min(deathrate2,7*deathrate))); % avoid unrealistic value (in particular on entities with low death toll)


% ##############################################################
% ##############################################################
% ##               R E N D E R    P L O T S                   ##
% ##############################################################
% ##############################################################
%% figure
if clearfig
    clf; currentfigure = gcf;
else
    currentfigure = figure;
end
formatfig(currentfigure,'figname',outputfile,'paperposition',[0.6345    0.6505   28.4084   19.7150],'paperorientation','landscape')
hold on
hp = zeros(6,1);
% confidence interval
CI = 0.95;
% cases
%hp(1) = plot(date,fun(b,time)/1000,'k-','LineWidth',2);
[ypred,delta] = nlpredci(@fun,time,b,residual,'Jacobian',Jacobian,'alpha',1-CI);
delta = max(0,delta); % remove NaN
ypredCI = cummax([ypred-delta;ypred+delta]',1);
ypredCImin = min(min(ypredCI));
hp(1) = plot(date,ypred/1000,'k-','LineWidth',2);
if (ypredCImin>=0) && max(delta)<0.9*max(ypred)
    plot(date,max(0,ypredCI/1000),'-','color',rgb('SlateGray'),'linewidth',1)
else
    ypredCI = ypredCI+NaN; %make NaN to for output
    dispf('\t--->no confidence interval for CASES in ''%s''',country)
end
ylm = get(gca,'Ylim');
xlm = get(gca,'Xlim');
www = xlm(2);
hhh = ylm(2);
tp2 = tpeak(infectedidx) - fix(2/r(infectedidx)) + first_day;
tp3 = tpeak(infectedidx) + fix(2/r(infectedidx)) + first_day;
tp4 = 2*tpeak(infectedidx) + first_day;
tp = tpeak(infectedidx) + first_day;
h = plot([tp,tp],[0,hhh],'r','LineWidth',1);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';

h = fill([tp2,tp3,tp3,tp2],[0 0 hhh hhh],rgb('Crimson'),'FaceAlpha',0.15,'EdgeColor','none');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';

h = fill([tp3,tp4,tp4,tp3],[0 0 hhh hhh],rgb('Orange'),'FaceAlpha',0.15,'EdgeColor','none');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';

h = fill([tp4,www,www,tp4],[0 0 hhh hhh],rgb('ForestGreen'),'FaceAlpha',0.15,'EdgeColor','none');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';

hp(2) = scatter(sampledDate,infected/1000,50,'k','filled');
h = scatter(sampledDate,infected/1000,30,'w','filled');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
datetick('x',2,'keepticks')
ylabel('Infected or Death toll \fontsize{14}\bf(\times1000)','fontsize',12)
xlabel('Date','fontsize',12)

% Infection rate
yyaxis right
pdfpred = dfun(b,time);
bint = nlparci(b,residual,'Jacobian',Jacobian,'alpha',1-CI);
pdfCI = [dfun(bint(:,1),time);dfun(bint(:,2),time)]';
pdfCImin = min(min(pdfCI));
hp(3) = plot(date,pdfpred,'b','LineWidth',1);
if (pdfCImin>=0) && max(max(pdfCI,[],2)-min(pdfCI,[],2))<2.5*max(pdfpred)
    plot(date,pdfCI,'-','LineWidth',0.3,'color',rgb('DodgerBlue'));
else
    dispf('\t--->no confidence interval for INFECTION RATE for ''%s''',country)
    pdfCI = pdfCI + NaN;
end
ylabel('Infected/day','fontsize',12)

% death toll
caserate = [0,diff(fun(b,time))];
caserate(date<=now) = deathrate * caserate(date<=now);
caserate(date>now) = deathrate2 * caserate(date>now);
caserate = cumsum(caserate);
yyaxis left
hp(4:6) = [
    plot(date,deathrate*fun(b,time)/1000,'k:','LineWidth',2)
    plot(date,caserate/1000,'k--','LineWidth',2)
    plot(sampledDate,Countrytotaldead(startidx:end)/1000,'s','markeredgecolor','k','markerfacecolor','none','markersize',6)
    ];

ax = gca;
ax.YAxis(1).Color = 'k';
ax.YAxis(2).Color = 'b';
legend(hp,{
    'Predicted Cases'
    'Actual Cases'
    'Infection Rate (cases/day)'
    'Death Estimates (average)'
    'Death Estimates (updated)'
    'Death Cases (actual)'},'Location','best')
title(sprintf('Coronavirus epidemic in %s',country))
grid on
title(sprintf('COVID-19 Epidemic Simulation for %s',country))
str1 = sprintf('Total Projected: \\bf%0.0f\\rm | Total Projected Death Toll: \\bf%0.0f\\rm | R^2 = %0.3f \n Date: \\bf%s\\rm | Total Cases: %0.0f | Total Deaths: %0.0f \n \\fontsize{9}\\itConfidence intervals at %0.3g %% ',...
    K(end),deathrate*K(end),mdl.Rsquared.Adjusted,datestr(sampledDate(infectedidx)),Countrytotalinfected(end),Countrytotaldead(end),CI*100);
T = text(min(get(gca,'xlim')), max(get(gca,'ylim')), str1);
set(T, 'fontsize', 10, 'verticalalignment', 'top', 'horizontalalignment', 'left');
hold off

%% prepare output: R
R = struct(...
             'country',country,...
             'datestr',datestr(sampledDate(infectedidx)),...
             'date',date,'ypred',ypred,'ypredCI',ypredCI,...
             'tp',[tp tp2 tp3 tp4],'www',www,'hhh',hhh,...
             'sampledDate',sampledDate,...
             'infected',infected,...
             'death',Countrytotaldead(startidx:end),...
             'pdfpred',pdfpred,...
             'pdfCI',pdfCI,...
             'deathrate',deathrate,...
             'deathest',deathrate*fun(b,time),...
             'deathest2',caserate,...
             'projectedinfected',K(end),...
             'projecteddeaths',deathrate*K(end),...
             'currentinfected',Countrytotalinfected(end),...
             'currentdeaths',Countrytotaldead(end) ...
             );
         
outtmp = struct(...
    'entity',regexprep(country,'^.*/',''),...
    'country',country,...
    'datestr',datestr(sampledDate(infectedidx)),...
    'data',R,...
    'fill',@() [
        fill([R.tp(2),R.tp(3),R.tp(3),R.tp(2)],[0 0 R.hhh R.hhh],rgb('Crimson'),'FaceAlpha',0.15,'EdgeColor','none')
        fill([R.tp(3),R.tp(4),R.tp(4),R.tp(3)],[0 0 R.hhh R.hhh],rgb('Orange'),'FaceAlpha',0.15,'EdgeColor','none')
        fill([R.tp(4),R.www,R.www,R.tp(4)],[0 0 R.hhh R.hhh],rgb('ForestGreen'),'FaceAlpha',0.15,'EdgeColor','none')
    ],...
    'scatter',@() [
         plot(R.sampledDate,R.infected/1000,'o','markersize',7.5,'markeredgecolor','k','markerfacecolor','w')
         plot(R.sampledDate,R.death/1000,'s','markeredgecolor','k','markerfacecolor','none','markersize',6)
            ],...
    'plot',@() [
        plot(R.date,R.ypred/1000,'k-','LineWidth',2)
        plot(R.date,max(0,R.ypredCI/1000),'-','color',rgb('SlateGray'),'linewidth',1)
        plot([1,1]*R.tp(1),[0,R.hhh],'r','LineWidth',1)
        plot(R.date,R.deathest/1000,'k:','LineWidth',2)
        plot(R.date,R.deathest2/1000,'k--','LineWidth',2)
                ],...
    'plot2',@() [
        plot(date,pdfpred/1000,'b','LineWidth',1)
        plot(R.date,R.pdfCI/1000,'-','LineWidth',0.3,'color',rgb('DodgerBlue'))
        ] ...
        );
        
         
%% print (aded OV)
formatax(gca,'fontsize',10), drawnow
print_pdf(600,get(gcf,'filename'),outputpath,'nocheck')
print_png(300,get(gcf,'filename'),outputpath,'',false,0,0,false)
saveas(gcf,fullfile(outputpath,[get(gcf,'filename') '.fig']),'fig')

%% Summary and predictions
fidout = fopen(fullfile(outputpath,sprintf('%s.txt',outputfile)),'w');
warning('on')
for fid=[1 fidout] %#ok<FXUP>
    fprintf(fid,'Summary on %s\n',datetxt);
    fprintf(fid,'---------------------\n');
    fprintf(fid,' country: %s\n',country);     % add country
    fprintf(fid,' date: %10s  day: %3d\n',datestr(sampledDate(infectedidx)),infectedidx);
    fprintf(fid,' start date: %s \n',datestr(first_day));
    fprintf(fid,' number of cases: %d\n',infected(infectedidx));
    fprintf(fid,' number of deaths: %d\n',max(Countrytotaldead));
    fprintf(fid,' estimated epidemic size (cases): %d\n',K(infectedidx));
    fprintf(fid,' estimated epidemic rate (1/day): %0.3g\n',r(infectedidx));
    fprintf(fid,' estimated initial state (cases): %d\n',C0(infectedidx));
    fprintf(fid,' estimated initial doubling time (day): %3.1f\n',round(log(2)/r(infectedidx),1));
    fprintf(fid,' estimated duration of fast growth phase (day): %d\n',tau(infectedidx));
    fprintf(fid,' estimated peak date: %s  day: %d \n',datestr(dpeak(infectedidx)),tpeak(infectedidx));
    fprintf(fid,' estimated peak rate (cases/day): %d\n',dCpeak(infectedidx));
    fprintf(fid,' estimated end of transition phase: %s  day: %d\n',datestr(dend(infectedidx)),2*tpeak(infectedidx));    
    
    tp2 = tpeak(infectedidx) - fix(2/r(infectedidx));
    tp3 = tpeak(infectedidx) + ceil(2/r(infectedidx));
    tp4 = 2*tpeak(infectedidx);  % end date
    tp = tpeak(infectedidx) ;
    if infectedidx < tp2
        phase = 1;
        txt = 'slow growth';
    elseif infectedidx < tp
        phase = 2;
        txt = 'fast growth acceleration phase';
    elseif infectedidx < tp3
        phase = 3;
        txt = 'fast growth deceleration phase';
    elseif infectedidx < tp4
        phase = 4;
        txt = 'slow growth transition phase';
    else
        phase = 5;
        txt = 'ending phase';
    end
    fprintf(fid,' epidemic phase: %g/5 %s\n\n',phase, txt);
    fprintf(fid,' estimated death rate[*]: %0.3g%% (averaged), %0.3g%% (projection)\n\n\n',deathrate*100,deathrate2*100);
    fprintf(fid,' [*] The value is highly indicative, based on reported deaths / reported cases\n');
    fprintf(fid,'     Since cases are counted and reported differently between entities,\n');
    fprintf(fid,'     death rates are expected to vary significantly between entities.\n\n');
    predict(fid) % forecasting
    fprintf(fid,'\n\n\n\n\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
    fprintf(fid,' DISCLAIMER: the presented results are derived from an automatized procedure\n');
    fprintf(fid,' without any human supervision. The results are provided "AS IS".\n');
    fprintf(fid,' Insufficient data and entities not entered in the exponential phase \n');
    fprintf(fid,' are known to lead to unreliable predictions.\n');
    fprintf(fid,' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n');
    fprintf(fid,'%s - %s@%s - %s version %0.4g\n',datestr(now),username,localname,mfilename,versn);
end
fclose(fidout);
if any(country{1}==',') || any(country{1}=='/')
    typ = 'province/region';
else
    typ = 'country';
end
fprintf('>> %s (v%0.4g):: %s ''%s'' analyzed in %0.4g s\n',mfilename,versn,typ,country,etime(clock,cT0));


%%
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    %       N E S T E D   F U N C T I O N S
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function [y,J] = fun(b,t) %-------------------------- logistic model
        % Logistic model
        y = b(1)./(1 + b(3)*exp(-b(2)*t));
        if nargout > 1
            nt = length(t);
            J = zeros(3,nt);
            J(1,1:nt) = 1./(1 + b(3)*exp(-b(2)*t));
            J(3,1:nt) = -b(1)*exp(-b(2)*t).*J(1,1:nt).^2;
            J(2,1:nt) = -b(3)*t.*J(3,1:nt);                                                                                                                                             
            J = J';
        end
    end

    function y = dfun(b,t)
        y = b(1)*b(2)*b(3)*exp(-b(2)*t)./(1 + b(3)*exp(-b(2)*t)).^2;
    end

    function predict(fid)  %-------------------------- logistic model
        if nargin<1, fid = 1; end
        K0 = K(infectedidx);
        r0 = r(infectedidx);
        A0 = A(infectedidx);
        
        fprintf(fid,'\nShort-term forecasting for %s\n',country);
        fprintf(fid,'%4s %12s %10s %10s %10s %10s %10s  %10s\n',...
            'day','date','actual','predict','error %','c./day act.','c./day pred.','error %');

        for n = infectedidx-5:infectedidx+daysforecast %#ok<FXUP>
            if n <= infectedidx
                if n==infectedidx, flag=sprintf('  <---- %s',lastdaywithdata); else, flag=''; end
                c0 = infected(n) - infected(n-1);
                c1 = round(funC(n-1),0)-round(funC(n-2),0);
                fprintf(fid,'%4d %12s %10d %10d %10.2f %10d %10d  %10.2f %s\n',n,datestr(n-1+first_day),infected(n),round(funC(n-1),0),...
                    abs(round(funC(n-1),0)/infected(n)-1)*100,...
                    c0,c1,abs(c1/c0-1)*100,flag) ;
            else
                fprintf(fid,'%4d %12s %10s %10d %10s %10s %10d\n',n,datestr(n-1+first_day),'-',...
                    round(funC(n-1),0),'-','-',round(funC(n-1),0)-round(funC(n-2),0));
            end
        end
        
        function C = funC(t)
            C = K0/(1+A0*exp(-r0*t));
        end
    end % ========== end predict

    function R2 = calcR2(t,C,b)
        K00 = b(1);
        r00 = b(2);
        A00 = b(3);
        zbar = sum(C)/length(C);
        SStot = sum((C - zbar).^2);
        SSres = sum((C - funC(t)).^2);
        R2 = 1 - SSres/SStot;
        function Ca = funC(t)
            Ca = K00./(1+A00*exp(-r00*t));
        end
    end

    function [b0] = iniGuess(C)
        b0 = [];
        n = length(C);
        if mod(n,2) == 0
            k1 = 1;
            k3  = n-1;
        else
            k1 = 1;
            k3 = n;
        end
        k2 = (k1 + k3)/2;
        m = k2 - k1 -1;
        q = C(k2)^2 - C(k3)*C(k1);
        if q <= 0
            return
        end
        p = C(k1)*C(k2) - 2*C(k1)*C(k3) + C(k2)*C(k3);
        if p <= 0
            return
        end
        K = C(k2)*p/q;
        r = log(C(k3)*(C(k2) - C(k1))/C(k1)/(C(k3) - C(k2)));
        if r < 0
            return
        end
        A = (C(k3) - C(k2))*(C(k2) - C(k1))/q*...
            (C(k3)*(C(k2) - C(k1))/C(k1)/(C(k3) - C(k2)))^((k3-m)/m);
        if A <= 0
            return
        end
        b0 = [K r A]';
    end

    function countrypath(currentcountry) % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-added OV 2020-03-22
        if nargin<1
            currentcountry = country;
            if isempty(listentities)
                country = regexprep(currentcountry,{'(^.)(.*$)' 'Us|Uk','Korea\,\ssouth'},{'${upper($1)}${lower($2)}' '${upper($0)}' 'Korea, South'});
            else
                currentcountry = country;
                country = listentities(ismember(lower(listentities),lower(currentcountry)));
                missingentities = setdiff(currentcountry,country);
                for imissing=1:length(missingentities)
                    fprintf('WARNING: the country ''%s'' does not exist\n',missingentities(imissing));
                end
                if isempty(country), error('no valid entities, nothing to do'), end
            end
        end
        outputpath = fullfile(globaloutputpath,char(currentcountry));
        if ~exist(outputpath,'dir'), mkdir(outputpath); end
    end % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
    
    function cleancountrypath % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-added OV 2020-03-22
        if exist(outputpath,'dir') && length(dir(outputpath))<3 % remove empty folder
            rmdir(outputpath)
        end
    end % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

    function checkcountry() % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- added OV 2020-03-21
        if ~ismember(country,listentities)
            fprintf('\nList of entities (%d)\n',length(listentities))
            for s = 'A':'Z'
                c = listentities(~cellfun(@isempty,regexp(listentities,['^' s],'once')));
                if ~isempty(c)
                    fprintf('[%s]: | %s\n',s,sprintf('%s | ',c{:}))
                end
            end
            cleancountrypath()
            error('the country (%s) does not exist, or no data is available. See list above',country)
        end
    end %+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

    function showentitiesintab() % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- added OV 2020-03-30
        mfile = mfilename;
        C = regexprep(listentities,'[,/].*$',''); %  entities without region/province/department
        L = regexprep(listentities,'(.).*$','$1'); % first country letter
        Lu = unique(L);
        [ntotalentities,ntotalentries]=deal(0);
        for iL=1:length(Lu)
            okC = C(L==Lu(iL));
            Cu = unique(okC);
            nCu = length(Cu);
            fprintf('%1$s %2$-54s------ %1$s ------%3$60s %1$s\n',sprintf('[%s]',Lu(iL)),sprintf('--- %d entities ---',nCu),'---');
            for iC=1:nCu
                ntotalentities = ntotalentities+1;
                ic = find(C==Cu(iC));
                nic = length(ic);
                if nic==1
                    fprintf('\t%-64s%s("%s")\n',listentities(ic),mfile,listentities(ic));
                    ntotalentries = ntotalentries+1;
                else
                    fprintf('\t** %s (%d entries) **\n',Cu(iC),nic);
                    for jc=1:length(ic)
                        fprintf('\t\t>  %-53s%s("%s")\n',listentities(ic(jc)),mfile,listentities(ic(jc)));
                        ntotalentries = ntotalentries+1;
                    end
                end
            end
        end
        fprintf('\n>>>> %s() can analyze %d entities (incl. ships), and %d regions/provinces/departments on %s\n',mfile,ntotalentities,ntotalentries,datetxt);
    end

    function GLOBALresults=loadGLOBALresults(when,initialize) % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- added OV 2020-03-27
        if nargin<1, when = yesterday; end
        if nargin<2, initialize = true; end
        % initialization
        if initialize
            if exist(GLOBALprefetchfile,'file') && ~nodailyreportprefetch
                load(GLOBALprefetchfile,'GLOBALresults') % load GLOBALresults
                lastday = GLOBALresults(end).date;
            else
                fprintf('\n\n%1$080d\n%2$45s\n%3$55s\n%1$080d\n',0,'UPDATE GLOBAL DATA',sprintf('Pause during %0.3g s (CTRL+C to cancel)',pauseduration)), pause(pauseduration)
                GLOBALresults = [];
                lastday = GLOBALfirstdate-1;
            end
            when = max(when);
            if when>lastday % if something to load
                res = loadGLOBALresults(lastday+1:when,false);
                if isempty(GLOBALresults)
                    GLOBALresults = res;
                else
                    GLOBALresults = [res;GLOBALresults];
                    [~,ind] = sort([GLOBALresults.date],'ascend');
                    GLOBALresults = GLOBALresults(ind);
                end
                save(GLOBALprefetchfile,'GLOBALresults')
            end
            return % end initialization
        end
        % pseudo-recursion
        nwhen = length(when);
        if nwhen>1
            GLOBALresults = repmat(loadGLOBALresults(when(1),false),nwhen,1);
            for ires=2:nwhen
                GLOBALresults(ires) = loadGLOBALresults(when(ires),false);
            end
            return
        end
        % result corresponding to a single date
        fres = GLOBALrawresultfile(when);
        if ~offline || ~exist(fres,'file')
            tmp = webread(URLrawresultGLOBAL(when));
            fid=fopen(fres,'w'); fprintf(fid,'%s\n',tmp); fclose(fid); % raw GLOBAL results
        end
        opts2 = detectImportOptions(fres, "TextType","string");
        opts2.VariableNamesLine = 1;
        opts2.DataLines = [2,inf];
        opts2.PreserveVariableNames = true;
        tmp = readtable(fres,opts2);
        tmp.Properties.VariableNames = opts2.VariableNames; % added 2020-04-25 / fix R2020a
        tokeep = true(height(tmp),1);
        isUS = tmp.('Country_Region')~="US";
        tokeep(isUS) = tmp.('Province_State')(isUS)~="US";
        tmp = tmp(tokeep,:); % original data
        tmp2=tmp;            % filtered data (not the same size)
        % fix names (variable nomenclature)
        tmp2.Province_State(ismissing(tmp2.Province_State) | tmp2.Province_State=="") = "none";
        tmp2.Province_State(tmp2.Province_State ==tmp2.Country_Region) = "none";
        % special BUG in data
        if (when==737873) %'23-Mar-2020'
            if tmp2.Confirmed(tmp2.Province_State=="French Polynesia")>1000 
                warning('Fix French Polynesia confused with Metropolitan France')
                tmp2.Province_State(tmp2.Province_State=="French Polynesia") = "";
            end
        end
        % main fixes
        tmp2(tmp2.Province_State == "Repatriated Travellers",:) = [];        
        tmp2(tmp2.Province_State == "Recovered",:) = [];
        tmp2(tmp2.Province_State == "Bavaria",:) = [];
        tmp2(tmp2.Province_State == "From Diamond Princess",:) = [];
        tmp2(tmp2.Country_Region == "MS Zaandam",:) = [];
        tmp2(tmp2.Country_Region == "Diamond Princess",:) = [];        
        tmp2(tmp2.Province_State == "Diamond Princess",:) = [];
        tmp2(tmp2.Province_State == "Diamond Princess cruise ship",:) = [];
        tmp2(tmp2.Province_State == "Cruise ship",:) = [];
        tmp2(tmp2.Province_State == "Cruise Ship",:) = [];
        tmp2(tmp2.Province_State == "London, ON",:) = [];
        tmp2(tmp2.Province_State == "Toronto, ON",:) = [];
        tmp2(tmp2.Province_State == "Montreal, QC",:) = [];
        tmp2(tmp2.Province_State == " Montreal, QC",:) = [];
        tmp2(tmp2.Province_State == "Edmonton, Alberta",:) = [];
        tmp2(tmp2.Province_State == "Calgary, Alberta",:) = [];
        tmp2.Country_Region(tmp2.Country_Region == "Cape Verde") = "Cabo Verde";        
        tmp2.Country_Region(tmp2.Country_Region == "Gambia, The") = "Gambia";
        tmp2.Country_Region(tmp2.Country_Region == "The Gambia") = "Gambia";
        tmp2.Country_Region(tmp2.Country_Region == "Bahamas, The") = "Bahamas";
        tmp2.Country_Region(tmp2.Country_Region == "The Bahamas") = "Bahamas";
        tmp2.Country_Region(tmp2.Country_Region == "Vatican City") = "Holy See";
        tmp2.Country_Region(tmp2.Country_Region == "Czechia") = "Czech Republic";
        tmp2.Country_Region(tmp2.Country_Region == "Iran (Islamic Republic of)") = "Iran";
        tmp2.Country_Region(tmp2.Country_Region == "Republic of Korea") = "Korea, South";
        tmp2.Country_Region(tmp2.Country_Region == "South Korea") = "Korea, South";
        tmp2.Country_Region(tmp2.Country_Region == "Republic of Ireland") = "Ireland";
        tmp2.Country_Region(tmp2.Country_Region == "Republic of Moldova") = "Moldova";
        tmp2.Country_Region(tmp2.Country_Region == "Russian Federation") = "Russia";
        tmp2.Country_Region(tmp2.Country_Region == "Ivory Coast") = "Cote d'Ivoire";
        tmp2.Country_Region(tmp2.Country_Region == "Republic of the Congo") = "Congo (Brazzaville)";
        tmp2.Country_Region(tmp2.Country_Region == "Viet Nam") = "Vietnam";
        % Palestine
        tmp2.Country_Region(tmp2.Country_Region == "occupied Palestinian territory") = "Palestine";

        % China
        tmp2.Country_Region(tmp2.Country_Region == "Taipei and environs") = "Taiwan";
        tmp2.Country_Region(tmp2.Country_Region == "Taiwan*") = "Taiwan";
        tmp2.Country_Region(tmp2.Country_Region == "Mainland China") = "China";
        tmp2.Province_State(tmp2.Country_Region == "Hong Kong") = "Hong Kong";
        tmp2.Country_Region(tmp2.Province_State == "Hong Kong") = "China";
        tmp2.Province_State(tmp2.Country_Region == "Hong Kong SAR") = "Hong Kong";
        tmp2.Country_Region(tmp2.Province_State == "Hong Kong SAR") = "China";
        tmp2.Province_State(tmp2.Country_Region == "Macau") = "Macau";
        tmp2.Country_Region(tmp2.Province_State == "Macau") = "China";
        tmp2.Province_State(tmp2.Country_Region == "Macao") = "Macau";
        tmp2.Country_Region(tmp2.Province_State == "Macao") = "China";
        % UK fixes
        tmp2.Country_Region(tmp2.Country_Region == "United Kingdom") = "UK";
        tmp2.Province_State(tmp2.Country_Region == "Gibraltar") = "Gibraltar";
        tmp2.Country_Region(tmp2.Province_State == "Gibraltar") = "UK";
        tmp2.Province_State(tmp2.Country_Region == "Jersey") = "Jersey";
        tmp2.Country_Region(tmp2.Province_State == "Jersey") = "UK";
        tmp2.Province_State(tmp2.Country_Region == "Guernsey") = "Guernsey";
        tmp2.Country_Region(tmp2.Province_State == "Guernsey") = "UK";
        tmp2.Province_State(tmp2.Country_Region == "Channel Islands") = "Channel Islands";
        tmp2.Country_Region(tmp2.Province_State == "Channel Islands") = "UK";        
        tmp2.Province_State(tmp2.Country_Region == "North Ireland") = "North Ireland";
        tmp2.Country_Region(tmp2.Province_State == "North Ireland") = "UK";
        tmp2.Province_State(tmp2.Country_Region == "Cayman Islands") = "Cayman Islands";
        tmp2.Country_Region(tmp2.Province_State == "Cayman Islands") = "UK";
        tmp2.Province_State(tmp2.Country_Region == "Falkland Islands (Islas Malvinas)") = "Falkland Islands (Islas Malvinas)";
        tmp2.Country_Region(tmp2.Province_State == "Falkland Islands (Islas Malvinas)") = "UK";
        % Netherlands
        tmp2.Province_State(tmp2.Country_Region == "Aruba") = "Aruba";
        tmp2.Country_Region(tmp2.Province_State == "Aruba") = "Netherlands";
        tmp2.Province_State(tmp2.Country_Region == "Curacao") = "Curacao";
        tmp2.Country_Region(tmp2.Province_State == "Curacao") = "Netherlands";
        % Denmark
        tmp2.Province_State(tmp2.Country_Region == "Greenland") = "Greenland";
        tmp2.Country_Region(tmp2.Province_State == "Greenland") = "Denmark";
        tmp2.Province_State(tmp2.Country_Region == "Faroe Islands") = "Faroe Islands";
        tmp2.Country_Region(tmp2.Province_State == "Faroe Islands") = "Denmark";
        % US
        tmp2.Province_State(tmp2.Country_Region == "Guam") = "Guam";
        tmp2.Country_Region(tmp2.Province_State == "Guam") = "US";
        tmp2.Province_State(tmp2.Country_Region == "Puerto Rico") = "Puerto Rico";
        tmp2.Country_Region(tmp2.Province_State == "Puerto Rico") = "US";
        % general fix
        tmp2.Province_State(tmp2.Province_State ==tmp2.Country_Region) = "none";
        % France
        tmp2.Province_State(tmp2.Province_State == "Saint Martin") = "St Martin";
        tmp2.Country_Region(tmp2.Country_Region == "Saint Martin") = "St Martin";
        tmp2.Country_Region(tmp2.Country_Region == "St. Martin") = "St Martin";
        tmp2.Province_State(tmp2.Country_Region == "St Martin") = "St Martin";
        tmp2.Province_State(tmp2.Country_Region == "Saint Barthelemy") = "Saint Barthelemy";
        tmp2.Province_State(tmp2.Country_Region == "Guadeloupe") = "Guadeloupe";
        tmp2.Province_State(tmp2.Country_Region == "Martinique") = "Martinique";
        tmp2.Province_State(tmp2.Country_Region == "Reunion") = "Reunion";
        tmp2.Province_State(tmp2.Country_Region == "Mayotte") = "Mayotte";
        tmp2.Province_State(tmp2.Country_Region == "French Guiana") = "French Guiana";
        tmp2.Province_State(tmp2.Country_Region == "Fench Guiana") = "French Guyana";
        tmp2.Country_Region(tmp2.Province_State == "St Martin") = "France";
        tmp2.Country_Region(tmp2.Province_State == "Saint Barthelemy") = "France";
        tmp2.Country_Region(tmp2.Province_State == "Guadeloupe") = "France";        
        tmp2.Country_Region(tmp2.Province_State == "Martinique") = "France";
        tmp2.Country_Region(tmp2.Province_State == "Reunion") = "France";
        tmp2.Country_Region(tmp2.Province_State == "Mayotte") = "France";
        tmp2.Country_Region(tmp2.Province_State == "French Guiana") = "France";
        tmp2.Country_Region(tmp2.Province_State == "Fench Guiana") = "France";        
        tmp2.Country_Region(  (tmp2.Province_State ~= "none") & ...
                               (tmp2.Country_Region == "France") ...
                            )="France DOM-TOM";
        tmp2.Province_State(  (tmp2.Province_State ~= "none") & ...
                               (tmp2.Country_Region == "US") ...
                            )="none";  % prepare the merging of US data (they are not cumulated in the original data)
        tmp2.Province_State = regexprep(tmp2.Province_State,{'^\s*','\s*$','None','Fench'},{'','','none','French'});
        tmp2.Country_Region = regexprep(tmp2.Country_Region,{'^\s*','\s*$'},{'',''});
        isothers = (tmp2.Country_Region == "others"); % %Diamond Princess, Grand Princess
        tmp2.Country_Region(isothers) = tmp2.Province_State(isothers);
        tmp2.Province_State(isothers) = "none";
        iserr = (tmp2.Country_Region == "French Guiana"); % mistake (inversion)
        tmp2.Province_State(iserr) = tmp2.Country_Region(iserr);
        tmp2.Country_Region(iserr) = "France";
        % unique names
        tmpid = tmp2(:,{'Province_State','Country_Region'});
        listGLOBALentities = unique(tmpid,'stable');
        nentities = height(listGLOBALentities);
        res = [listGLOBALentities repmat({NaN NaN NaN NaN},nentities,1)]; % 'Lat' 'Long_' 'Confirmed' 'Deaths'
        res.Properties.VariableNames={'Province/State','Country/Region','lat','lon','cases','deaths'};
        % fieldnames are variables according to the date, a regular expression secure naming
        itmpfields = cellfun(@(key) find(~cellfun(@isempty,regexp(tmp2.Properties.VariableNames,key))) ,...
            {'^[Ll]at' '^[Ll]on' '^[Cc]onf' '^[Dd]ea'},'UniformOutput',false);
        for ientity=1:nentities %#ok<FXUP>
            j=ismember(tmpid,listGLOBALentities(ientity,:));
            % Lat and Lon if supplied
            if ~isempty(itmpfields{1}), res{ientity,3} = mean(tmp2{j,itmpfields{1}}); end
            if ~isempty(itmpfields{2}), res{ientity,4} = mean(tmp2{j,itmpfields{2}}); end
            % Confirmed cases and deaths
            res(ientity,5:6) = {sum(tmp2{j,itmpfields{3}}) sum(tmp2{j,itmpfields{4}})};
        end
        
        % main output
        GLOBALresults = struct(...
            'date',when,...
            'datestr',regexprep(datestr(when,'mm/dd/yy'),'(^)0|(/)0','$1'),...
            'data',res ...
            );

    end %+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

    function USresults=loadUSresults(d,initialize) % +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- added OV 2020-03-27
        if nargin<1, d = yesterday; end
        if nargin<2, initialize = true; end
        % initialization
        if initialize
            if exist(USprefetchfile,'file') && ~nodailyreportprefetch
                load(USprefetchfile,'USresults') % load USresults
                lastday = USresults(end).date;
            else
                fprintf('\n\n%1$080d\n%2$45s\n%3$55s\n%1$080d\n',0,'UPDATE US DATA',sprintf('Pause during %0.3g s (CTRL+C to cancel)',pauseduration)), pause(pauseduration)
                USresults = [];
                lastday = USfirstdate-1;
            end
            d = max(d);
            if d>lastday % if something to load
                res = loadUSresults(lastday+1:d,false);
                if isempty(USresults)
                    USresults = res;
                else
                    USresults = [res;USresults];
                    [~,ind] = sort([USresults.date],'ascend');
                    USresults = USresults(ind);
                end
                save(USprefetchfile,'USresults')
            end
            return % end initialization
        end
        % pseudo-recursion
        nd = length(d);
        if nd>1
            USresults = repmat(loadUSresults(d(1),false),nd,1);
            for ires=2:length(d)
                USresults(ires) = loadUSresults(d(ires),false);
            end
            return
        end
        % result corresponding to a single date
        fres = USrawresultfile(d);
        if ~offline && ~exist(fres,'file')
            tmp = webread(URLrawresultUS(d));
            fid=fopen(fres,'w'); fprintf(fid,'%s\n',tmp); fclose(fid); % raw US results
        end
        opts2 =  detectImportOptions(fres, "TextType","char");
        opts2.VariableNamesLine = 1;
        opts2.DataLines = [2,inf];
        opts2.PreserveVariableNames = true;
        tmp = readtable(fres,opts2);
        tmp.Properties.VariableNames = opts2.VariableNames; % added 2020-04-25 / fix R2020a
        tmp = tmp(tmp.('Country_Region')=="US",:);
        % remove data with no geographical attibutes (time evolution has non sense, no posssibility to generate 
        listUSstates = setdiff(unique(tmp.Province_State),{'Recovered' 'Diamond Princess' 'Grand Princess' 'Federal Bureau of Prisons' 'US Military','Veteran Hospitals','Repatriated Travellers'});
        nentities = length(listUSstates);
        res = [listUSstates repmat({NaN NaN NaN NaN},nentities,1)]; % 'Lat' 'Long_' 'Confirmed' 'Deaths'
        % fieldnames are variables according to the date, a regular expression secure naming
        itmpfields = cellfun(@(key) find(~cellfun(@isempty,regexp(tmp.Properties.VariableNames,key))) , {'^[Ll]at' '^[Ll]on' '^[Cc]onf' '^[Dd]ea'});
        for ientity=1:nentities
            j=ismember(tmp.Province_State,listUSstates{ientity});
            res(ientity,2:5) = {nanmean(tmp{j,itmpfields(1)}) nanmean(tmp{j,itmpfields(2)}) sum(tmp{j,itmpfields(3)}) sum(tmp{j,itmpfields(4)})};
            if isnan(res{ientity,2}) || isnan(res{ientity,3}), error('[%d/%d] no Lat/Lon data for %s/%s, please check',ientity,nentities,tmp.Country_Region{find(j,1)},tmp.Province_State{find(j,1)}); end
        end
        USresults = struct(...
            'date',d,...
            'datestr',regexprep(datestr(d,'mm/dd/yy'),'(^)0|(/)0','$1'),...
            'data',cell2table(res,'VariableNames',{'state','lat','lon','cases','deaths'}) ...
            );

    end %+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

    function FRresults=loadFRANCEresults() % ------------------- LOADFRANCERESULTS, added OV 2020-03-29
        if exist(FRprefetchfile,'file') && ~noglobalprefetch
            load(FRprefetchfile,'FRresults')
        else % regenerate data (few seconds required)
            fprintf('\n\n%1$080d\n%2$45s\n%3$53s\n%1$080d\n',0,'UPDATE FRANCE DATA',sprintf('Pause during %0.3g s (CTRL+C to cancel)',pauseduration)), pause(pauseduration)
            % France centers of gravity of departements
            % see: http://www.ign.fr/publications-de-l-ign/Institut/Publications/IGN_Magazine/82/IGN_MAG_82.pdf#15
            % source: https://fr.wikipedia.org/wiki/Liste_des_centres_de_gravit%C3%A9_des_d%C3%A9partements_m%C3%A9tropolitains_fran%C3%A7ais
            GPSFRdpt = cell2table({
                1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	'2A'	'2B'	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
                46.0994	49.5594	46.3936	44.1061	44.6636	43.9375	44.7517	49.6156	42.9208	48.3044	43.1033	44.2803	43.5433	49.0997	45.0511	45.7181	45.7808	47.0647	45.3569	41.8636	42.3942	47.4247	48.4411	46.0903	45.1042	47.1653	44.6842	49.1136	48.3875	48.2611	43.9933	43.3586	43.6928	44.8253	43.5797	48.1544	46.7778	47.2581	45.2633	46.7283	43.9656	47.6167	45.7269	45.1281	47.3614	47.9119	44.6242	44.3675	44.5172	47.3908	49.0794	48.9492	48.1094	48.1467	48.7869	48.9894	47.8464	49.0372	47.1153	50.4472	49.4103	48.6236	50.4936	45.7258	43.2567	43.0531	42.6	48.6708	47.8586	45.8703	47.6411	46.6447	47.9944	45.4775	46.0344	48.8567	49.655	48.6267	48.815	46.5556	49.9581	43.7853	44.0858	43.4606	43.9939	46.6747	46.5639	45.8917	48.1967	47.8397	47.6317	48.5222	48.8472	48.9175	48.7775	49.0828
                5.3489	3.5583	3.1883	6.2439	6.2631	7.1164	4.4247	4.6408	1.5039	4.1617	2.4142	2.6797	5.0864	-0.3636	2.6686	0.2017	-0.6744	2.4911	1.8769	8.9881	9.2064	4.7722	-2.8642	2.0189	0.7414	6.3617	5.1681	0.9961	1.3703	-4.0589	4.1803	1.1728	0.4533	-0.5753	3.3672	-1.6386	1.5758	0.6914	5.5761	5.6978	-0.7839	1.4294	4.1658	3.8064	-1.6822	2.3442	1.6047	0.4603	3.5003	-0.5642	1.3275	4.2386	5.2264	-0.6581	6.165	5.3817	-2.81	6.6633	3.5047	3.2206	2.4253	0.1289	2.2886	3.1408	-0.7614	0.1639	2.5222	7.5514	7.2742	4.6414	6.0861	4.5422	0.2222	6.4436	6.4281	2.3422	1.0264	2.9333	1.8417	-0.3172	2.2778	2.1661	1.2819	6.2181	5.1861	-1.2978	0.4603	1.2353	6.3806	3.5644	6.9286	2.2431	2.2458	2.4781	2.4689	2.1311
                }','VariableNames',{'code','lat','lon'}); % Code Lat Lon
            GPSFRdpt.code(cellfun(@isnumeric,GPSFRdpt.code))  = cellfun(@(c) sprintf('%02d',c),GPSFRdpt.code(cellfun(@isnumeric,GPSFRdpt.code)),'UniformOutput',false);
            
            % France centers of gravity of regions
            %Source: http://www.ign.fr/institut/actus/lign-calcule-position-centre-geographique-13-nouvelles-regions-metropolitaines
            GPSFRreg = cell2table({
                'Grand Est' 'Nouvelle-Aquitaine' 'Auvergne-Rhône-Alpes' 'Bourgogne-Franche-Comté' 'Bretagne' 'Centre-Val de Loire' 'Corse' 'Île-de-France' 'Occitanie' 'Hauts-de-France' 'Normandie' 'Pays de la Loire' 'Provence-Alpes-Côte d''Azur'
                1 2 3 4 5 6 7 8 9 10 11 12 13
                48.6892000000000 45.1922000000000 45.5158000000000 47.2353000000000 48.1797000000000 47.4806000000000 42.1497000000000 48.7092000000000 43.7022000000000 49.9661000000000 49.1211000000000 47.4747000000000 43.9550000000000
                5.61940000000000 0.197800000000000 4.53810000000000 4.80920000000000 -2.83860000000000 1.68530000000000 9.10530000000000 2.50470000000000 2.13720000000000 2.77530000000000 0.106700000000000 -0.823900000000000 6.05330000000000
                }','VariableNames',{'name','code','lat','lon'}); % Name  Dpt Lat Lon
            
            % load data
            if ~offline && ~exist(URLrawresultFR,'file')
                tmp = webread(URLrawresultFR);
                fid=fopen(FRrawresultfile,'w'); fprintf(fid,'%s\n',tmp); fclose(fid); % raw US results
            end
            opts2 =  detectImportOptions(FRrawresultfile, "TextType","char");
            opts2.VariableNamesLine = 1;
            opts2.DataLines = [2,inf];
            opts2.PreserveVariableNames =   true;
            tmp = readtable(FRrawresultfile,opts2);
            tmp.Properties.VariableNames = opts2.VariableNames; % added 2020-04-25 / fix R2020a
            % cleaning new fields added on April 1st, 2020
            for checkprop = {'cas_confirmes' 'cas_ehpad' 'cas_confirmes_ehpad' 'cas_possibles_ehpad' 'deces' 'deces_ehpad' 'reanimation' 'hospitalises' 'gueris' 'depistes'}
                if iscellstr(tmp.(checkprop{1})) || isstring(tmp.(checkprop{1}))
                    tmp.(checkprop{1}) = str2double(tmp.(checkprop{1}));
                end
                tmp.(checkprop{1}) = max(0,tmp.(checkprop{1}));
            end
%             tmp.depistes = max(0,str2double(tmp.depistes));
%             tmp.gueris = max(0,str2double(tmp.gueris));
%             tmp.hospitalises = max(0,str2double(tmp.hospitalises));
%             tmp.reanimation = max(0,str2double(tmp.reanimation));
%             tmp.deces = max(0,str2double(tmp.deces));
% 
%             tmp.cas_ehpad = max(0,str2double(tmp.cas_ehpad));
%             tmp.cas_confirmes_ehpad = max(0,str2double(tmp.cas_confirmes_ehpad));
%             tmp.cas_possibles_ehpad = max(0,str2double(tmp.cas_possibles_ehpad));
%             tmp.deces_ehpad = max(0,str2double(tmp.deces_ehpad));
            
            % parsing
            tmp1 = tmp(ismember(tmp.granularite,'departement'),:);
            tmp1 = tmp1(cellfun(@length,tmp1{:,3})<=6,:); % remove DOM-TOM considered in global results from WHO
            [rawcodedepartement,ind,idpt] = unique(tmp1.maille_code); %#ok<ASGLU>
            codedepartement = regexprep(rawcodedepartement,'^DEP-','');%<-------- departmental data
            listdepartements = tmp1.maille_nom(ind);
            listdepartementswithcode = cellfun(@(c,d) sprintf('%s,%s',c,d),codedepartement,listdepartements,'UniformOutput',false);
            ndepartements = length(listdepartements);
            tmp2 = tmp(ismember(tmp.granularite,'region'),:); %<----------------- regional data
            tmp2 = tmp2(ismember(tmp2.maille_nom,GPSFRreg.name),:);
            listregions = unique(tmp2.maille_nom);
            nregions=length(listregions);
            tmp3 = tmp(ismember(tmp.granularite,'pays'),:);  %<------------------ national data incl. EHPAD data
            % prefetch
            tmp1.datenum = datenum(tmp1.date);
            tmp2.datenum = datenum(tmp2.date);
            tmp3.datenum = datenum(tmp3.date);            
            d1 = unique(tmp1.datenum); nd1 = length(d1);
            d2 = unique(tmp2.datenum); nd2 = length(d2);
            d3 = unique(tmp3.datenum); nd3 = length(d3);
            [~,~,jdpt]=intersect(codedepartement,GPSFRdpt.code,'stable');
            [~,~,jreg]=intersect(listregions,GPSFRreg.name,'stable');
            resdptmt = cell2table(...
                [listdepartementswithcode num2cell(GPSFRdpt.lat(jdpt),2) num2cell(GPSFRdpt.lon(jdpt),2) repmat({0 0},ndepartements,1)],...
                'VariableNames',{'state','lat','lon','cases','deaths'});
            resreg = cell2table(...
                [listregions num2cell(GPSFRreg.lat(jreg),2) num2cell(GPSFRreg.lon(jreg),2) repmat({0 0},nregions,1)],...
                'VariableNames',{'state','lat','lon','cases','deaths'});
            resnat = cell2table({'France' NaN NaN 0 0},'VariableNames',{'state','lat','lon','cases','deaths'});
            FRresults = struct(...
                'departement',repmat(struct('date',[],'datestr','','data',resdptmt),nd1,1),...
                'region',repmat(struct('date',[],'datestr','','data',resreg),nd2,1),...
                'national',repmat(struct('date',[],'datestr','','data',resnat),nd3,1) ...
                );
            % populate departements (note that some data are missing)
            % EHPAD data not including in this set
            [previouscases, previousdeaths] = deal(zeros(ndepartements,1));
            for idate=1:nd1 %#ok<FXUP>
                okdate = (tmp1.datenum==d1(idate));
                for idpt=1:ndepartements
                    okdpt = find(okdate & ismember(tmp1.maille_code,rawcodedepartement(idpt)));
                    FRresults.departement(idate).date = d1(idate);
                    FRresults.departement(idate).datestr = regexprep(datestr(d1(idate),'mm/dd/yy'),'(^)0|(/)0','$1');
                    if any(okdpt)
                        currentcases = max(previouscases(idpt),max(tmp1(okdpt,:).cas_confirmes));
                        currentdeaths = max(previousdeaths(idpt),max(tmp1(okdpt,:).deces));
                        estimatedcases = max(tmp1(okdpt,:).deces) + max(tmp1(okdpt,:).hospitalises) + max(tmp1(okdpt,:).gueris) +...
                            max(tmp1(okdpt,:).depistes); % new field added on Apr 1st, 2020
                    else
                        currentcases = previouscases(idpt);
                        currentdeaths = previousdeaths(idpt);
                        estimatedcases = 0;
                    end
                    % record
                    FRresults.departement(idate).data.cases(idpt) = max(estimatedcases,currentcases);
                    FRresults.departement(idate).data.deaths(idpt) = currentdeaths;
                    % update
                    previouscases(idpt) =  FRresults.departement(idate).data.cases(idpt);
                    previousdeaths(idpt) =  FRresults.departement(idate).data.deaths(idpt);
                end
            end
            % populate regions (note that some data are missing)
            % EHPAD data not including in this set
            [previouscases, previousdeaths] = deal(zeros(nregions,1));
            for idate=1:nd2 %#ok<FXUP>
                okdate = (tmp2.datenum==d2(idate));
                for ireg=1:nregions
                    okreg = find(okdate & ismember(tmp2.maille_nom,listregions(ireg)));
                    FRresults.region(idate).date = d2(idate);
                    FRresults.region(idate).datestr = regexprep(datestr(d2(idate),'mm/dd/yy'),'(^)0|(/)0','$1');
                    if any(okreg)
                        currentcases = max(previouscases(ireg),max(tmp2(okreg,:).cas_confirmes));
                        currentdeaths = max(previousdeaths(ireg),max(tmp2(okreg,:).deces));
                        estimatedcases = max(tmp2(okreg,:).deces) + max(tmp2(okreg,:).hospitalises) + max(tmp2(okreg,:).gueris)+...
                            max(tmp2(okreg,:).depistes); % added on Apr 1st, 2020
                    else
                        currentcases = previouscases(ireg);
                        currentdeaths = previousdeaths(ireg);
                        estimatedcases = 0;
                    end
                    % record
                    FRresults.region(idate).data.cases(ireg) = max(estimatedcases,currentcases);
                    FRresults.region(idate).data.deaths(ireg) = currentdeaths;
                    % update
                    previouscases(ireg) = FRresults.region(idate).data.cases(ireg);
                    previousdeaths(ireg) =  FRresults.region(idate).data.deaths(ireg);
                end
            end
            % populate national results
             [previouscases, previousdeaths] = deal(0);
             for idate=1:nd3 %#ok<FXUP>
                 okdate = (tmp3.datenum==d3(idate));
                 if any(okdate)
                     FRresults.national(idate).date = d3(idate);
                     FRresults.national(idate).datestr = regexprep(datestr(d3(idate),'mm/dd/yy'),'(^)0|(/)0','$1');
                     if any(okdate)
                         currentcases = max(previouscases, max(tmp3(okdate,:).cas_confirmes) + max(tmp3(okdate,:).cas_confirmes_ehpad));
                         currentdeaths = max(previousdeaths,max(tmp3(okdate,:).deces) + max(tmp3(okdate,:).deces_ehpad));
                         estimatedcases = max(tmp3(okdate,:).deces) +...
                             max(tmp3(okdate,:).hospitalises) +...
                             max(tmp3(okdate,:).gueris) + ...
                             max(tmp3(okdate,:).depistes); % added on Apr 1st, 2020
                     else
                         currentcases = previouscases;
                         currentdeaths = previousdeaths;
                         estimatedcases = 0;
                     end
                     % record
                     FRresults.national(idate).data.cases = max(previouscases,max(estimatedcases,currentcases));
                     FRresults.national(idate).data.deaths = max(previousdeaths,currentdeaths);
                     % update
                     previouscases = FRresults.national(idate).data.cases;
                     previousdeaths =  FRresults.national(idate).data.deaths;
                 else
                     warning('[%d/%d] France results:: bad date',idate,nd3)
                 end
             end
       
            % save
            save(FRprefetchfile,'FRresults')
        end

    end % ------------------- LOADFRANCERESULTS

% output
if nargout, out = outtmp; end

end

%% PRIVATE FUNCTIONS (from INRA/MS, Olivier Vitrac)

function print_pdf(resolution,fichier,chemin,options,varargin)
%PRINT_PDF print interactively the current figure as PDF document
%   print_PDF([resolution,filename,fullpath,options])
%   print_PDF generates a standard pdf based on figure properties
%   print_PDF(resolution,fichier,options)
%   print_PDF(resolution,fichier,chemin,options)
%
%   default options = '-loose'
%   options 'nocheck' forces printing without any dialog
%   options 'append' is equivalent to nocheck but with an append feature (via a first PostScript 2.1 printing)
%   for 3D without PLOTOBJ, options = '-opengl' is recommened (see PRINT for details)
%
% Advance use of 'append' (PDF conversion from PostScript is performed for efficiency only when the last page is printed)
%       print_pdf(resolution,fichier,chemin,options,'pagenum',currentpage,'pagestart',1,'pagestop',finalpage)
%           currentpage = page index to print (when currentpage = pagestart value, the existing postscript file is removed)
%             finalpage = page to initiate PS-to-PDF printing
%   Example of use (print all figures in a single PDF)
%       listoffigs = sort(findobj('type','figure'))'; numfig = length(listoffigs);
%       for i=1:length(listoffigs)
%           figure(listoffigs(i)), print_pdf(600,'myPDFfile.pdf',pwd,'append','pagenum',i,'pagestart',1,'pagestop',numfig)
%       end

% Woodox 1.0/MS 2.1 - 13/08/07 - Olivier Vitrac - rev. 16/02/18

% revision
% 25/07/07 add options and fix chemin option
% 13/08/07 fix path ambiguity when both fichier and chemin contain a path
% 11/12/07 add fileinfo, preview and dialog
% 24/06/09 add nocheck
% 12/12/09 add evince &  for unix systems (in replacement of winopen)
% 20/07/11 fix chemin with 'nocheck'
% 21/07/11 new fix for 'nocheck'
% 22/07/11 fix extension with 'nocheck'
% 03/03/12 add append
% 12/03/12 help improvements, add example using append
% 04/06/17 fix figure number for R2013b and later
% 16/02/18 fix example


% definitions
resolution_default = 600; % PDF
ext_default = '.pdf';
currentfig = gcf;
if verLessThan('matlab','9.0'), currentfignum = currentfig; else, currentfignum = currentfig.Number; end
options_default = '-loose';
pmode = get(currentfig,'paperpositionmode');
psdefault = struct('pagenum',NaN,'pagestart',1,'pagestop',Inf);

% arg check
if nargin<1, resolution = []; end
if nargin<2, fichier = ''; end
if nargin<3, chemin = ''; end
if nargin<4, options = ''; end
psoptions = argcheck(varargin,psdefault);
if isempty(resolution), resolution = resolution_default; end
if isempty(fichier), [~,fichier] = fileparts(get(currentfig,'filename')); end
if isempty(fichier), [~,fichier] = fileparts(get(currentfig,'name')); end
if isempty(fichier), fichier = sprintf('Figure_%d',currentfignum); end
if isempty(options), options = options_default; end
if any(chemin) 
    if chemin(1)=='-', options = chemin; chemin = ''; end
    [pathstr,name,ext] = fileparts(fichier);
    if ~isempty(pathstr) && exist(fullfile(chemin,pathstr),'dir')
        chemin = fullfile(chemin,pathstr);
    end
else
    [chemin,name,ext] = fileparts(fichier);
end
if isempty(chemin), chemin = pwd; end
if ~exist(chemin,'dir'), error('the path ''%s'' does not exist',chemin), end
if ~strcmp(ext,'.pdf'), ext = ext_default; end
fichier = fullfile(chemin,[name ext]);

% printing with no check
if strcmpi(options,'nocheck')
    print(['-r' int2str(resolution)],'-dpdf','',fichier), return
elseif strcmpi(options,'append')
    tmppsfile = regexprep(fichier,['\' ext_default '$'],'.ps','ignorecase');
    okps2pdf = true;
    if ~isnan(psoptions.pagenum) && (psoptions.pagenum>=psoptions.pagestart)
        if isinf(psoptions.pagestop) || isnan(psoptions.pagestop)
             psmsg = sprintf('page %d (first page=page %d)',psoptions.pagenum,psoptions.pagestart);
        else, psmsg = sprintf('page %d/%d (first page=page %d)',psoptions.pagenum,psoptions.pagestop,psoptions.pagestart);
             okps2pdf = psoptions.pagenum>=psoptions.pagestop;
        end
    else,     psmsg = '(specify ''pagenum'', ''pagestart'', ''pagestop'' for optimization)';
    end
    if exist(tmppsfile,'file') && ~isnan(psoptions.pagenum) && (psoptions.pagenum<=psoptions.pagestart), delete(tmppsfile), end
    % PS print
    dispf('PRINTPDF %s with options=''append'': Postscript level 2 printing...',psmsg)
    print('-append',['-r' int2str(resolution)],'-dpsc2','',tmppsfile)
    % PS to PDF conversion (if required)
    if okps2pdf
        dispf('PRINTPDF %s with options=''append'': PS->PDF...',psmsg)
        if exist(fichier,'file'), delete(fichier), end
        ps2pdf('psfile', tmppsfile, 'pdffile', fichier, 'gspapersize', 'a4');
        if ~isnan(psoptions.pagestop) && ~isinf(psoptions.pagestop)  && (psoptions.pagenum>=psoptions.pagestop)
            delete(tmppsfile)
            fileinfo(fichier)
        end
    else
        dispf('PRINTPDF %s with options=''append'': conversion PS->PDF postponed...',psmsg)
        dispf('\t ps2pdf(''psfile'',''%s'',''pdffile'',''%s'',''gspapersize'',''a4'')\n',tmppsfile,fichier);
    end
    return
end

% printing with controls
printon = true;
ok = ~exist(fichier,'file');
while ~ok
    answer = questdlg({sprintf('the file ''%s'' already exist',[name ext]),sprintf('in ''%s''',chemin)},'Overwrite an existing PDF file','overwrite','new file','cancel','overwrite');
    if strcmp(answer,'new file')
        [fichier,chemin] = uiputfile('*.pdf','Choose a new filename for your pdf');
        [chemin,name,ext] = fileparts(fullfile(chemin,fichier));
        if ~strcmp(ext,'.pdf'), ext = ext_default; end
        fichier = fullfile(chemin,[name ext]);
        ok = ~exist(fichier,'file');
    else
        printon = ~strcmp(answer,'cancel');
        ok = true;
    end
end
if printon && ~strcmp(pmode,'auto')
    answer=questdlg({ sprintf('Current ''paperpositionmode'' is ''%s''',pmode)
            ' '
            'If you set it as ''auto'', you can change the paper layout'
            'by resizing the figure and you can check the result by refreshing'
            'the preview in the printpreview panel.'
            'With ''auto'' printing starts only after you close the the preview panel'
            'You do not need to change the options in the preview panel.'
            ' '
            'Select an option' },'PAPERPOSITIONMODE: auto or not ?','auto','no change','cancel','no change');
    if strcmp(answer,'auto')
        set(currentfig,'paperpositionmode','auto');
    elseif strcmp(answer,'cancel')
        printon = false;
    else
        dispf('no change, paperpositionmode: %s',pmode)
    end
end
 
if printon
    if get(currentfig,'paperpositionmode')
        disp('Printing starts only after you close the the preview panel...')
        uiwait(printpreview)
    end
    start = clock;
    dispf('printing ''%s''....',fichier)
    print(['-r' int2str(resolution)],'-dpdf',options, fichier)
    dispf('... end in %0.3g s',etime(clock,start))
    fileinfo(fichier)
    if isunix
       system(['evince ' fichier ' &']);
    else
        winopen(fichier)
    end
else
    disp('PRINT_PDF canceled')
end
end

function print_png(resolution,filename,filepath,options,flipon,margin,negativeon,horizontalon)
%PRINT_PNG  print active window as a PNG image
%   print_png(resolution,filename[,filepath,options])
%   print_png(resolution,filename,options)
%   print_png(resolution,filename,options [,flipon,margin,negativeon])
%       options must be followed by '-' (ex. '-opengl')
%       imfile,flipon,margin,negativeon are parameters of pngtruncateim

% INRA\MS 2.0 - 28/01/01 - Olivier Vitrac - 10/05/14

% revision
% 25/07/07 add options and fix filepath option
% 13/08/07 fix path ambiguity when both filename and filepath contain a path
% 09/09/12 add pngtruncateim
% 10/05/14 force flipon when paperorientation is set to landscape

% Default
flipon_default = false;
margin_default = 0;
negativeon_default = false;
horizontalon_default = false;

% arg check
if nargin<2, error('Two arguments are required'); end
if nargin<3, filepath = ''; end
if nargin<4, options = ''; end
cropon = (nargin>=5);
if nargin<5, flipon = []; end
if nargin<6, margin = []; end
if nargin<7, negativeon= []; end
if nargin<8, horizontalon = []; end
if isempty(flipon), flipon = flipon_default; end
if ~flipon && strcmpi(get(gcf,'PaperOrientation'),'landscape') && verLessThan('matlab','9.5'), flipon = true; end
if isempty(margin), margin = margin_default; end
if isempty(negativeon), negativeon = negativeon_default; end
if isempty(horizontalon), horizontalon = horizontalon_default; end

% do the print
if any(filepath)
    if filepath(1)=='-', options = filepath; filepath = ''; end
    [~,name,ext] = fileparts(filename);
else
    [filepath,name,ext] = fileparts(filename);
end
if isempty(ext), ext = '.png'; end
if isempty(filepath), filepath = cd; end
if ~exist(filepath,'dir'), error('the path ''%s'' does not exist',filepath), end
filename = [name ext];
print (['-r' int2str(resolution)],'-dpng',options, fullfile(filepath,filename))
if cropon
    pngtruncateim(fullfile(filepath,filename),flipon,margin,negativeon,horizontalon) % varargin{:})
else
    dispf('to truncate the generated PNG image, use the following code:\npngtruncateim(''%s'',0,0,0)',fullfile(filepath,filename))
end
end

function out = pngtruncateim(imfile,flipon,margin,negativeon,horizontalon)
%PNGTRUNCATEIM  crop and rotate PNG already saved images
%   Syntax: pngtruncateim(imfile [,flip,margin,negative,horizontal])
%       default values: flip = false
%                       margin = 100
%                       negative = false
%                       horizontal = false

% Migration 2.0 - 24/05/11 - INRA\Olivier Vitrac - rev. 15/10/15

% Revision history
% 27/05/11 guess extension
% 07/11/11 add negative
% 29/01/12 accept image data
% 29/03/12 fix negativeon
% 02/10/12 add horizontalon
% 16/03/15 iterate when imfile is a cell array
% 15/10/15 enable imfile as directory

% default values
flipon_default = false;
margin_default = 100;
negativeon_default = false;
horizontalon_default = false;

% arg check
if nargin<2, flipon = []; end
if nargin<3, margin = []; end
if nargin<4, negativeon= []; end
if nargin<5, horizontalon = []; end
if isempty(flipon), flipon = flipon_default; end
if isempty(margin), margin = margin_default; end
if isempty(negativeon), negativeon = negativeon_default; end
if isempty(horizontalon), horizontalon = horizontalon_default; end

% main section
if ischar(imfile)
    if exist(imfile,'file')~=2
        if exist(imfile,'dir')
            pngtruncateim(explore('*.png',imfile,0,'fullabbreviate'),0,0,0,0)
            return
        else
            error('the file ''%s'' does not exist',imfile)
        end
    end
    im = imread(imfile);
    nofile = false;
elseif isnumeric(imfile)
    im = imfile;
    nofile = true;
elseif iscellstr(imfile)
    out = cellfun(@(f) pngtruncateim(f,flipon,margin,negativeon,horizontalon),imfile);
    return
else
    error('the argument must be a valid image filename or image data')
end

siz = size(im);
if negativeon, im = 255-im; end
imb = min(im,[],3);
lim = zeros(2,2);
dimlist = 1:2;
for dim = dimlist
    lim(dim,1) = find(min(imb,[],dimlist(mod(dim,2)+1))<255,1,'first');
    lim(dim,2) = find(min(imb,[],dimlist(mod(dim,2)+1))<255,1,'last');
end
lim(:,1) = max(lim(:,1)-margin,1);
lim(:,2) = min(lim(:,2)+margin,siz(1:2)');
im = im(lim(1,1):lim(1,2),lim(2,1):lim(2,2),:);
if flipon, im = flipdim(permute(im,[2 1 3]),2); end
if negativeon, im = 255-im; end
if horizontalon && size(im,1)>size(im,2), im = flipdim(permute(im,[2 1 3]),2); end

% save final result
if nofile
    out = im;
else
    ext = uncell(regexp(imfile,'\.([^\.]+)$','tokens'));
    if isempty(ext)
        imwrite(im,imfile,'png');
    else
        imwrite(im,imfile,ext{1});
    end
    if nargout, out = true; end
end
end

function [param,remaining]=argcheck(list,propdefault,keywordlist,varargin)
%ARGCHECK check arguments passed as keywords or as pairs property/value, possible conflicts (override) are managed via precedence rules
%    Syntax: param = argcheck(list,propertylist [,keywordlist,'case','property'])
%   Options: [param,remaining] = argcheck(...)
%    INPUTS
%           list: list of arguments coded in a cell array (Note that structure arguments are exapanded in list)
%                 Since 28/02/12, structures are not anymore expanded at the end of the list and the original order is kept.
%                 Note that this behavior can have effects on some existing functions using argcheck (please check).
%                 Use the keyword 'nosort' to returns to the previous behavior.
%                 e.g. argcheck({struct('a',false,'b',true) 'a' true},struct('a',false,'b',true)) returns a=FALSE
%                      argcheck({struct('a',false,'b',true) 'a' true},struct('a',false,'b',true),'','nosort') returns a=TRUE
%    propdefault: list of pair property/defaultvalue coded in a structure (propdefault.property=value)
%          propdefault.keyword = defaultvalue
%                 Note that in case of multiple definitions, the latter are ignored (this behavior is modified with 'property'):
%                 e.g. argcheck({'a',true,'a',false},struct('a',false)) returns a=TRUE
%    keywordlist: list (cell array of strings) of accepted keywords (use [] to omit it)
%         'case': force case sensitive
%     'property': force property precedence (instead of keyword)
%                 Note that in this case the last value overrides any previous definition:
%                 e.g. argcheck({'a',false,'a',true},struct('a',false),'a','property') returns a=TRUE
%        'order': order fields
%         'keep': keep undefined properties setup in list (e.g. argcheck({'d','D','e','E','f','F'},struct('a','A','b','B','c','C'),'','keep'))
%       'nosort': returns to the behavior before 28/02/12: structures are expanded at the end of the list
%'nostructexpand': prevent structure to be expanded as a list of parameters
%'NaNequalmissing': NaN are interpreted as missing values
%    OUTPUTS
%          param: recognized keywords and pair property/value (empty values are replaced by their default values)
%            param.keyword = true (if keyword was definded) and false otherwise
%            param.property=value (if property exist), use propertylist.defaultvalue if needed (e.g. when value is set to [])
%      remaining: cell array containing un-used/not recognized parameters
%
% NB: keywords have a higher precedence on properties (a missing keyword is always set to false whatever the value assigned in propdefault)
%     e.g. argcheck({'a',false},struct('a',false),'a') returns TRUE
%
% TIP: Do not define the same parameter as both a parameter and a keyword with 'property' active, it will create a bug (it is a normal behavior, perhaps to be fixed in the future)

% MS 2.1 - 25/12/09 - INRA/Olivier Vitrac - rev. 12/08/18

% Revision history
% 27/12/09: case insensitive, several fixes
% 04/01/10: fix ambiguities between keywords and properties
% 01/02/10: fix behavior when no keywords were defined
% 05/02/10: fix syntax typing error
% 09/02/10: fix list as a structure and not a cell
% 09/02/10: fix argcheck([],[],'keyword')
% 09/02/10: add casesensitive
% 12/02/10 fix [a,b]=argcheck({'remaining input'},[],{'cmd','trun'})
% 13/02/10 add 'order'
% 14/02/10 add 'inherit','keep'
% 17/02/10 add 'nostructexpand'
% 24/02/10 fix case for propertylist and when option 'keep' is used
% 08/04/11 add 'NaNequalmissing'
% 28/02/12 add three notes on conflict resolution when multiple definitions are found and when 'property' is used
% 28/02/12 sort mixed structures list to keep precedence rules when user override (with multiple definitions) are used
% 18/02/16 add isnan2() for cell inputs combined with 'NaNequalmissing'
% 12/08/18 cosmetic for better compatibility with R2018 and a new TIP added

% definitions
options_list = {'case' 'property' 'order' 'keep' 'inherit' 'nostructexpand' 'NaNequalmissing' 'nosort'}; % inherit is a private property

% argcheck
if nargin<1,  error('SYNTAX: [param,remaining]=argcheck(list,propdefault [,keywordlist,''case'',''property''])'), end
if nargin<2, propdefault = []; end
if nargin<3, keywordlist = []; end
if isempty(varargin)
    options = struct('case',false,'property',false,'order',false,'keep',false,'inherit',false,'nostructexpand',false,'NaNequalmissing',false,'nosort',false);
else
    options = argcheck(varargin,[],options_list);
end
list = list(:); nlist = length(list);
if isempty(propdefault)
    %provisional implementation (OV 12/02/10)
    %dispf('\tto be used only in DEBUG mode')
    %if (mod(nlist,2)>0), error('list of property/value is unpaired'), end
    if ~isempty(list) && (mod(nlist,2)==0) && iscellstr(list(1:2:end)) 
        if ~iscellstr(list(1:2:end-1)), error('improper key values in list'), end
        propdefault = cell2struct(repmat({[]},nlist/2,1),list(1:2:end-1),1);
    end
end
if options.inherit, param = list; else param = struct([]); end
if ~iscell(list), list = {list}; end

% expand all structures
if ~options.nostructexpand
    isastruct = cellfun(@isstruct,list);
    nistruct = ~isastruct; %find(~istruct);
    istruct = find(isastruct);
    nstruct = length(istruct);
    if nstruct>0
        [keys,values,pos] = deal(cell(nstruct,1));
        for i=1:nstruct
            keys{i} = fieldnames(list{istruct(i)});
            pos{i} = ones(2*length(keys{i}),1)*istruct(i); % orginal position is set for the keyword and value
            values{i} = struct2cell(list{istruct(i)});
        end
        tab = [ cat(1,keys{:}) cat(1,values{:}) ]';
        list = cat(1,list(nistruct),tab(:)); % unsorted (before 28/02/12)
        if ~options.nosort
            [~,order] = sort(cat(1,find(~isastruct),cat(1,pos{:})),'ascend'); % all positions are assembled as list is and then sorted
            list = list(order); % sorted (after 28/02/12)
        end
        nlist = length(list);
    end
end

% character inputs
ichar   = find(cellfun(@(x)ischar(x),list));
if ~options.case, list(ichar) = lower(list(ichar)); end
iremaining = true(1,nlist);
    
% search keywords
if ~isempty(keywordlist)
    if ~iscell(keywordlist), keywordlist={keywordlist}; end
    if ~iscellstr(keywordlist), error('list of keywords must be a cell array of char'), end
    for keyword = keywordlist(:)'
        if options.case
            ifound = ismember(list(ichar),keyword{1});
        else
            ifound = ismember(list(ichar),lower(keyword{1}));
        end
        if any(ifound)
            param(1).(keyword{1})=true;
            iremaining(ichar(ifound)) = false;
        else
            param(1).(keyword{1})=false;
        end
    end
end

% search pairs property/value
if ~isempty(propdefault)
    propertylist = fieldnames(propdefault)';
    if options.case, proptosearch = propertylist; else, proptosearch = lower(propertylist); end
    if ~iscellstr(propertylist), error('list of properties must be a cell array of char'), end %#ok<ISCLSTR>
    iprop = 0;
    for prop = proptosearch
        iprop = iprop + 1;
        ifound = find(ismember(list(ichar),prop));
        if any(ifound)
            for j=1:length(ifound)
                if ~isfield(param,propertylist{iprop}) || options.property
                    if ichar(ifound(j))<nlist && ~(options.NaNequalmissing && any(isnan2(list{ichar(ifound(j))+1})))
                        param(1).(propertylist{iprop}) = list{ichar(ifound(j))+1};
                        iremaining(ichar(ifound(j))+[0 1]) = false;
                    elseif j<2
                        param(1).(propertylist{iprop}) =[];
                        iremaining(ichar(ifound(j))) = false;
                    end
                end
            end
        else
            param(1).(propertylist{iprop}) = [];
        end
        if isempty(param(1).(propertylist{iprop})), param(1).(propertylist{iprop}) = propdefault.(propertylist{iprop}); end
    end
end

% keep properties if required
if options.keep
    if options.case % create a structure from remaining parameters
        paramroot = argcheck(list(iremaining),[],'','nostructexpand','case'); % case is used
    else
        paramroot = argcheck(list(iremaining),[],'','nostructexpand'); % no case
    end
    param = argcheck(param,paramroot,'','case','inherit');
end
    
% ouput
if options.order, param = orderfields(param); end
if nargout>1, remaining = list(iremaining); end

end % end function

%% %%%% PRIVATE FUNCTIONS %%%%
function b=isnan2(x)
    if isnumeric(x), b=isnan(x);
    elseif ischar(x), b=false;
    elseif iscell(x), b = cellfun(@isnan2,x);
    else b = false(size(x));
    end
end

function u=username()
%USERNAME return the username using the system command whoami
% MS 2.1 - 02/02/08 - INRA\Olivier  - rev. 08/06/09
    if isunix
        [r,u]=system('whoami');
        u = u((u>='A') & (u<='z'));
    else
        u = getenv('USERNAME');
    end
end

function name = localname()
% LOCALNAME return the netbios name of the local machine
% MS 1.0 - 26/10/04 - INRA\Olivier Vitrac - rev. 03/02/08
    if isunix
        [r,name] = system('hostname');
    else
        name = getenv('COMPUTERNAME');
    end
    name = name(name>32);
end

function [rout,strout] = fileinfo(filename,ppath,dispon,stringonly,noerror)
    %FILEINFO returns a structure that contains the main information of a given file
    %  Syntax: info = fileinfo(filename,[path])
    %  Options: [info,str] = fileinfo(filename,[path,dispon,stringonly,noerror])
    %
    % INPUTS
    %   filename can include a full path or not
    %   path is used if filename does not include it (the current path is used otherwise)
    %   dispon (flag, default=true) prints on screen if true
    %   stringonly (flag, default=false) returns a formated single string instead of a structure if true 
    %   noerror (flag, default=false) no error message if true, corresponding files appear with sizes 0
    %   
    % OUTPUTS
    %   info = structure matching dir/ls properties (string if stringonly is true)
    %   str = formated string
    %
    %   See also: explore, isformat, isformat_multiple, ls, dir

    % MS-MATLAB 1.0 - 20/04/04 - INRA\Olivier Vitrac - rev. 30/12/12

    % revision history
    % 19/08/04: display if no ouput, added str
    % 11/09/07 fix filename as cell array
    % 12/01/11 fix versn in FILEPARTS for Matlab later than 7.11 (not supported any more)
    % 30/12/12 cosmetic modifications
    % 09/02/16 improved help, add noerror

    % arg check
    if nargin<2, ppath = ''; end
    if nargin<3, dispon = []; end
    if nargin<4, stringonly = []; end
    if nargin<5, noerror = []; end
    if ~nargout || isempty(dispon), dispon = true; end
    if isempty(stringonly), stringonly = false; end
    if isempty(noerror), noerror = false; end
    oldmatlab = verLessThan('matlab','7.11');
    issomethingmissing = false;

    if iscell(filename)
        m = length(filename);
        str = cell(m,1);
        [r,str{1}] = fileinfo(filename{1},ppath,dispon);
        for i=2:m, [r(i),str{i}] = fileinfo(filename{i},ppath,dispon); end
        dispon = false;
    else
        % other checks
        if oldmatlab
            [pathstr,name,ext,versn] = fileparts(filename);
        else
            [pathstr,name,ext] = fileparts(filename);
            versn = [];
        end
        if isempty(pathstr)
            if isempty(ppath)
                pathstr = cd;
            else
                pathstr = ppath;
            end
        end
        if ~exist(pathstr,'dir')
            if noerror
                issomethingmissing = true;
            else
                error('the directory ''%s'' does not exist',pathstr)
            end
        end
        fullfilename = fullfile(pathstr,[name ext versn]);
        if ~exist(fullfilename,'file')
            if noerror
                issomethingmissing = true;
            else
                error('the file ''%s'' does not exist in ''%s''',[name ext versn],pathstr)
            end
        end

        if issomethingmissing % faked info
            f = struct('date',datestr(now),'bytes',0);
        else % true info
            f = dir(fullfilename);
        end
        r =struct(   'filename', [name ext versn], ...
            'name',     name, ...
            'ext',      ext, ...
            'ver',      versn, ...
            'date',     f.date, ...
            'bytes',    f.bytes, ...
            'path',     pathstr...
            );

        if f.bytes>1024*1024
            str = sprintf('\t%s%s\t\t%s\t\t%0.1f MBytes\t\t%s',name,ext,f.date,f.bytes/(1024*1024),pathstr);
        elseif f.bytes>1024
            str = sprintf('\t%s%s\t\t%s\t\t%0.1f kBytes\t\t%s',name,ext,f.date,f.bytes/1024,pathstr);
        else
            str = sprintf('\t%s%s\t\t%s\t\t%d Bytes\t\t%s',name,ext,f.date,f.bytes,pathstr);
        end
    end

    % output
    if nargout>0, rout = r; end
    if stringonly, rout = str; return, end 
    if nargout>1, strout = str; end
    if dispon, fprintf('%s\n',str), end
end