1/* wuenic_ver_3.pl version 3.
    2
    3Implements WHO & UNICEF rules for estimating national
    4infant immunization coverage. Includes explanations and
    5grade of confidence in estimate.
    6
    7Based on methods described in:
    8
    9  Burton A, Monash R, Lautenbach B, Gacic-Dobo M, Neill M, Karimov
   10  R, Wolfson L, Jones G, Birmingham M. WHO and UNICEF estimates of
   11  national infant immunization coverage: methods and processes.
   12  Bull World Health Organ 2009; 87:535-541.
   13  http://www.who.int/bulletin/volumes/87/7/08-053819.pdf
   14
   15  Burton A, Gacic-Dobo M, Karimov R, Kowalski R.
   16  A computational logic-based representation of the WHO and UNICEF
   17  estimates of national immunization coverage. DRAFT 23 January 2011.
   18
   19  Articles and code available at: http://sites.google.com/site/wuenic/
   20
   21Author: Tony BURTON
   22        System Analyst
   23        Strategic Information Group
   24        Expanded Pogramme on Immunization
   25        Department of Immunization, Vaccines, and Biologicals
   26        World Health Organization
   27        1211 Geneva 27
   28        Switzerland
   29        BurtonA@who.int
   30
   31depends upon: xsb 3.1 prolog compiler; library: lists, standard.
   32http://xsb.sourceforge.net/
   33
   34Created:       6 November 2011.
   35Last update:   6 May 2016. ipv1 & rcv1 added, GoC for rcv1 based on mcv1 or mcv2 GoC
   36
   37Naming conventions:
   38    predicates names: lower_case_underscore_separated
   39	function names:   lower_case_underscore_separated
   40    variable names:   InterCaps
   41    literals:         camelCase
   42
   43Note that code is constructed "top down" with higher level goals appearing first.
   44Auxaliarity predicates and utility routines appear at the end of the module.
   45
   46%TBD: verify meaning of "new vaccine" for modified right handed sawtooth criteria.
   47%TBD: add assign indirect value for country vaccine year (anchor? wuenic?, both?)
   48%TBD: Delete (or fix) bubble up GoC.
   49%TBD: Add green line
   50
   51working group decisions:	comment
   52							ignoreGov, acceptGov, *modifyGov,
   53							ignoreAdmin, acceptAdmin, *modifyAdmin,
   54							ignoreReported, acceptReported, *modifiyReported,
   55							ignoreSurvey, acceptSurvey, modifySurvey,
   56							interpolate, calibrate, *use_reported,
   57							assignAnchor, assignWUENIC, assignGoC.
   58
   59Input atomic sentences:
   60    admin(country,vaccine,year,percent coverage).
   61    gov(country,vaccine,year,percent coverge).
   62    legacy(country,vaccine,year,percent coverage).
   63    survey_results(country,vaccine,year,surveyid,description,percent coverage).
   64    estimate_required(country,vaccine,year,presentation,comment).
   65    wgd(country,vaccine,year1,year2,action,explanation1,coverageID1,coverageAssign1,coverageID2,coverageAsign2)
   66	births_UNPD(country,year,births)
   67	si_UNPD(country,year,births_surviving_to_age_one)
   68
   69WHO and UNICEF working group memebers as of 1 January 2010 - 30 April 2012:
   70        Dr David BROWN, UNICEF/New York (dbrown@unicef.org)
   71        Mr Tony BURTON, WHO/Geneva  (burtona@who.int)
   72        Ms Marta GACIC-DOBO, WHO/Geneva  (gacicdobom@who.int)
   73        Mr Rouslan KARIMOV, UNICEF/NEW YORK (rkarimov@unicef.org)
   74        Dr Robert KOWALSKI, Imperial College London (r.kowalski@imperial.ac.uk)
   75 */

   76
   77% use xsb tabling feature to increases inference time.
   78% ---------------------------------------------------
   79:- table anchor_point/6.
   80:- table survey/5.
   81:- table reported/5.
   82:- table wuenic_I/6.
   83
   84:- op(500,xfy,:).
   85
   86% sawtooth_threshold : difference in increase/decrease between Y+/-1 and Y in reported data
   87% survey_reported_threshold : differerence between survey results and reported data.
   88% reported_calibrated_threshold : if mixed anchor points swithc between using reported and calibarated.
   89% -----------------------------------------------------------------------------------------------------
   90sawtooth_threshold(10).
   91survey_reported_threshold(10).
   92confidence_survey_scope(2).
   93confidence_survey_threshold(10).
   94confidence_UNPD_threshold(10).
   95
   96% establish relationship between name of first and third dose of vaccine.
   97% used in estimating recall bias.
   98% ------------------------------
   99vaccine(dtp3,dtp1).
  100vaccine(pol3,pol1).
  101vaccine(hib3,hib1).
  102vaccine(hepb3,hepb1).
  103vaccine(pcv3,pcv1).
  104
  105% Imported from data
  106%
  107% admin(country,vaccine,year,coverage).
  108% gov(country,vaccine,year,coverage).
  109% vaccinated(country,vaccine,year,vaccinated).
  110% target(country,vaccine,year,vaccinated).
  111% legacy(country,vaccine,year,coverage).
  112% survey_results(country,vaccine,year,id,description,coverage).
  113% wgd(country,vaccine,year1,year2,action,explanation,covid1,covass1,covid2,covass2).
  114% births_UNPD(country,year,births).
  115% si_UNPD(country,year,surviving_infants).
  116
  117% Load country-specific predicates describing data, survey_results,
  118% working group decisions and whether an estimate is required.
  119% Call estimate/0 to create country-specific estimates.
  120% -----------------------------------------------------
  121% :- ['data.pl'].
  122% :- estimate.
  123
  124% ==================================================
  125% top level predicate. Creates and outputs estimates.
  126% ==================================================
  127estimate :-
  128
  129	% Unify country code, country name and date for output.
  130	% ----------------------------------------------------
  131	country(CountryCode,CountryName),
  132	date(Date),
  133
  134	% Collect the set of all final estimates in list Estimate.
  135      % -------------------------------------------------------
  136	setof([
  137		CountryName,
  138		Date,
  139		CountryCode,
  140		Vaccine,
  141		Year,
  142		Coverage,
  143		PrevRev,
  144		GC,
  145		Admin,
  146		Gov,
  147		Reported,
  148		Vaccinated,
  149		Target,
  150		UnpdBirths,
  151		UnpdSI,
  152		SeriesValue,
  153		Source,
  154		SurveyInfo,
  155		Rule,
  156		Explanation],
  157
  158	   wuenic(
  159		CountryCode,
  160		Vaccine,
  161		Year,
  162		Rule,
  163		Explanation,
  164		Coverage,
  165		PrevRev,
  166		GC,
  167		Admin,
  168		Gov,
  169		Reported,
  170		Vaccinated,
  171		Target,
  172		UnpdBirths,
  173		UnpdSI,
  174		Source,
  175		SeriesValue,
  176		SurveyInfo),
  177
  178		Estimates),
  179
  180	open_out_file(OutFile,'wuenic.out',	'Country\tProductionDate\tISOCountryCode\tVaccine\tYear\tWUENIC\tWUENICPreviousRevision\tGradeOfConfidence\tAdministrativeCoverage\tGovernmentEstimate\tReportedCoverage\tChildrenVaccinated\tChildrenInTarget\tBirthsUNPD\tSurvivingInfantsUNPD\tReportedTimeSeries\tReportedTimeSeriesSource\tSurveyInformation\tRule\tComment\t'),
  181	output_results(Estimates,OutFile), close(OutFile),
  182	crypto_file_hash('wuenic.out', Hash, []),
  183	writeln(Hash).
  184
  185% Final estimate where there are wuenic values
  186% --------------------------------------------
  187wuenic(C,V,Y,Rule,Explanation,Coverage,PrevRev,GC,Admin,Gov,Reported,Vaccinated,Target,UnpdBirths,UnpdSI,Source,SeriesValue,SurveyInfo) :-
  188	estimate_required(C,V,Y,_,_),
  189	wuenic_I(C,V,Y,Rule,Explain,Cov),
  190	bound_0_100(Cov,Coverage),
  191
  192	confidence(C,V,Y,Rule,Coverage,GoCExplanation,GC),
  193	%assign_GoC(C,V,Y,Rule,Coverage,GoCExplanation,GC),
  194
  195	collect_data(C,V,Y,PrevRev,Admin,Gov,Reported,Vaccinated,Target,UnpdBirths,UnpdSI,SeriesValue,Source,SurveyInfo),
  196	change_from_previous_revision(C,V,Y,Coverage,Change),
  197	collect_explanations(C,V,Y,Text),
  198	my_concat_atom([Explain,' ',Text,' ',Change,' ',GoCExplanation],Explanation).
  199
  200% ------------------------------------------------------------------------------
  201%  End of wuenic top level routine.
  202
  203confidence(C, rcv1, Y, Rule, Coverage, Explanation, GC),
  204    estimate_required(C, rcv1, Y, _, mcv2)
  205 => assign_GoC(C, mcv2, Y, Rule, Coverage, Explanation, GC).
  206
  207
  208confidence(C, rcv1, Y, Rule, Coverage, Explanation, GC)
  209 => assign_GoC(C, mcv1, Y, Rule, Coverage, Explanation, GC).
  210
  211confidence(C, V, Y, Rule, Coverage, Explanation, GC)
  212 => assign_GoC(C, V, Y, Rule, Coverage, Explanation, GC).
  213
  214% GoC = 1 is low confidence (1 star)
  215% GoC = 3 is high confidence (3 stars)
  216assign_GoC(C, V, Y, _Rule, _, Support, GoC),
  217    workingGroupDecision(C, V, Y, assignGoC, Explanation, _, G)
  218 => GoC = G,
  219    my_concat_atom(['GoC=Assigned by working group. ', Explanation], Support).
  220
  221% Three stars: supported by reported data, survey and denominator
  222assign_GoC(C, V, Y, _Rule, _Coverage, Support, GoC),
  223    goc_reported_condition(C, V, Y, 'R+'),
  224    goc_survey_condition(C, V, Y, 'S+'),
  225    goc_denominator_condition(C, V, Y, 'D+')
  226 => GoC = '3',
  227    Support = 'GoC=R+ S+ D+'.
  228
  229% Two stars
  230assign_GoC(C, V, Y, _Rule, _Coverage, Support, GoC),
  231    two_sources(C, V, Y, _, S),
  232    not(challenge(C, V, Y, _, _, _))
  233 => GoC = '2',
  234    Support = S.
  235
  236assign_GoC(C, V, Y, _Rule, _Coverage, Support, GoC),
  237    one_source(C, V, Y, _, S),
  238    not(challenge(C, V, Y, _, _, _))
  239 => GoC = '2',
  240    Support = S.
  241
  242assign_GoC(C, V, Y, Rule, Coverage, Support, GoC),
  243    challenge(C, V, Y, Rule, Coverage, _)
  244 => GoC = '1',
  245    setof(Evidence, challenge(C, V, Y, Rule, Coverage, Evidence), List),
  246    my_concat_atom(['Estimate challenged by: ', List], Support).
  247
  248assign_GoC(C, V, Y, Rule, Coverage, Support, GoC),
  249    no_data(C, V, Y, Rule, Coverage, S)
  250 => Support = S,
  251    GoC = '1'.
  252
  253goc_reported_condition(C, V, Y, Support) :-
  254    reported(C, V, Y, _, _),
  255    wuenic_I(C, V, Y, R, _, _),
  256    (	memberchk(R, ['R:', 'R: AP'])
  257    ->	Support = 'R+'
  258    ;	Support = 'R-'
  259    ).
  260
  261% Is there any challenging survey?
  262goc_survey_condition(C, V, Y, Support),
  263    challenging_survey_in_scope(C, V, Y)
  264 => Support = 'S-'.
  265
  266% If there isn't, is there any supporting survey?
  267goc_survey_condition(C, V, Y, Support),
  268    supporting_survey_in_scope(C, V, Y)
  269 => Support = 'S+'.
  270
  271% No survey within scope
  272goc_survey_condition(_C, _V, _Y, _Support)
  273 => fail.
  274
  275% Old comment: simplify previous rule to
  276% supporting_survey_in_scope(C,V,Y,Rule) :-
  277%     survey(C,V,Y,_,_),
  278%     wuenic_I(C,V,Y,'S: AP',_,_).
  279
  280% rewrite rule to look at relationship between estimate rule and
  281% surveys in scope rule rather than difference in COV and SurveyCoverage
  282supporting_survey_in_scope(C, V, Y) :-
  283    estimate_required(C, V, Year, _, _),
  284    survey(C, V, Year, _, Survey),
  285    confidence_survey_scope(Scope),
  286    abs(Y - Year) =< Scope,
  287    confidence_survey_threshold(Threshold),
  288    wuenic_I(C, V, Y, _, _Explanation, Coverage),
  289    abs(Coverage - Survey) =< Threshold.
  290
  291% rewrite to look for surveys that challenge surveys.
  292% Modified 2017-05-05 Burton
  293%
  294% todo: merge support and challenge
  295challenging_survey_in_scope(C, V, Y) :-
  296    estimate_required(C, V, Year, _, _),
  297    survey(C, V, Year, _, Survey),
  298    confidence_survey_scope(Scope),
  299    abs(Y - Year) =< Scope,
  300    confidence_survey_threshold(Threshold),
  301    wuenic_I(C, V, Y, _, _Explanation, Coverage),
  302    abs(Coverage - Survey) > Threshold.
  303
  304goc_denominator_condition(C, V, Y, Support) :-
  305    goc_unpd_recal(C, V, Y),
  306    recal_unpd(C, V, Y, Recalc),
  307    confidence_UNPD_threshold(Threshold),
  308    wuenic_I(C, V, Y, _Rule, _Explain, Coverage),
  309    (	abs(Coverage - Recalc) < Threshold
  310    ->	Support = 'D+'
  311    ;	Support = 'D-'
  312    ).
  313
  314% Ensure unpd data exist.
  315goc_unpd_recal(C, V, Y) :-
  316    vaccinated(C, V, Y, _),
  317    births_UNPD(C, Y, _),
  318    si_UNPD(C, Y, _).
  319
  320% Recalculate coverage using reported number of children vaccinated
  321% and births and surviving infants from UNPD estimates.
  322% Births used for bcg and hepb birth dose, surviving infants
  323% for remaining vaccines.
  324recal_unpd(C, V, Y, Coverage),
  325    member(V, [bcg, hepbb])
  326 => vaccinated(C, V, Y, Vaccinated),
  327    births_UNPD(C, Y, Births),
  328    Coverage is Vaccinated / Births * 100. % Todo: remove 100
  329
  330recal_unpd(C, V, Y, Coverage)
  331%    member(V,
  332%      [dtp1, dtp3, pol3, ipv1, mcv1, mcv2, rcv1, hepb3, hib3, rotac,
  333%       pcv3, yfv]),
  334 => vaccinated(C, V, Y, Vaccinated),
  335    si_UNPD(C, Y, SI),
  336    Coverage is Vaccinated / SI * 100.
  337
  338	% Supported by two data sources, no challenge
  339	% --------------------------------------------
  340	two_sources(C,V,Y,_Rule,'GoC=R+ S+') :-
  341		goc_reported_condition(C,V,Y,'R+'),
  342		goc_survey_condition(C,V,Y,'S+').
  343	two_sources(C,V,Y,_Rule,'GoC=S+ D+') :-
  344		goc_survey_condition(C,V,Y,'S+'),
  345		goc_denominator_condition(C,V,Y,'D+').
  346	two_sources(C,V,Y,_Rule,'GoC=R+ D+') :-
  347		goc_reported_condition(C,V,Y,'R+'),
  348		goc_denominator_condition(C,V,Y,'D+').
  349
  350	% Supported by single source, no challenge
  351	% -------------------------------------
  352	one_source(C,V,Y,_Rule,'GoC=R+') :-
  353		not(two_sources(C,V,Y,_,_)),
  354		goc_reported_condition(C,V,Y,'R+').
  355	one_source(C,V,Y,_Rule,'GoC=S+') :-
  356		not(two_sources(C,V,Y,_,_)),
  357		goc_survey_condition(C,V,Y,'S+').
  358	one_source(C,V,Y,_Rule,'GoC=D+') :-
  359		not(two_sources(C,V,Y,_,_)),
  360		goc_denominator_condition(C,V,Y,'D+').
  361
  362	% Empirical evidence challenges estimate.
  363	% --------------------------------------
  364	challenge(C,V,Y,_Rule,_Coverage,Condition) :-
  365		(goc_reported_condition(C,V,Y,'R-'),Condition = 'R-');
  366		(goc_survey_condition(C,V,Y,'S-'),Condition = 'S-');
  367		(goc_denominator_condition(C,V,Y,'D-'),Condition = 'D-').
  368
  369	% No empirical supporting evidence.
  370	% --------------------------------
  371	no_data(C,V,Y,_Rule,_Coverage,'GoC=No accepted empirical data') :-
  372		not(goc_reported_condition(C,V,Y,_)),
  373		not(goc_survey_condition(C,V,Y,_)),
  374		not(goc_denominator_condition(C,V,Y,_)).
  375
  376	change_from_previous_revision(C,V,Y,Coverage,'') :-
  377		legacy(C,V,Y,PreviousCoverage),
  378		PreviousCoverage = Coverage.
  379
  380	change_from_previous_revision(C,V,Y,Coverage,Change) :-
  381		legacy(C,V,Y,PreviousCoverage),
  382		not(PreviousCoverage = Coverage),
  383		my_concat_atom(['Estimate of ',Coverage,' percent changed from previous revision value of ',PreviousCoverage,' percent. '],Change).
  384
  385	change_from_previous_revision(C,V,Y,_,'') :-
  386		not(legacy(C,V,Y,_)).
  387
  388% Estimate for non-DTP1 & RCV1 vaccines.
  389% ---------------------------------------
  390wuenic_I(C,V,Y,Rule,Explanation,Coverage) :-
  391	wuenic_II(C,V,Y,Rule,Explanation,Coverage),
  392	not(workingGroupDecision(C,V,Y,assignWUENIC,_,_,_)),
  393	not(member(V,['dtp1','rcv1'])).
  394
  395% Estimate for DTP1 where DTP3 estimate <= DTP1 estimate.
  396% --------------------------------------------------------
  397wuenic_I(C,dtp1,Y,Rule,Explanation,Coverage) :-
  398	wuenic_II(C,dtp1,Y,Rule,Explanation,Coverage),
  399	not(workingGroupDecision(C,dtp1,Y,assignWUENIC,_,_,_)),
  400	wuenic_II(C,dtp3,Y,_,_,DTP3Coverage),
  401	DTP3Coverage =< Coverage.
  402
  403% Estimate for DTP1 where DTP3 > DTP1: RMF
  404% -----------------------------------------
  405wuenic_I(C,dtp1,Y,'RMF:',Explanation,RMFCoverage) :-
  406	wuenic_II(C,dtp1,Y,_,_,DTP1Coverage),
  407	not(workingGroupDecision(C,dtp1,Y,assignWUENIC,_,_,_)),
  408	wuenic_II(C,dtp3,Y,_,_,DTP3),
  409	DTP3 > DTP1Coverage,
  410	RMFCoverage is round(18.258 - 0.6088*DTP3 - 0.0058*DTP3*DTP3),
  411	my_concat_atom(['DTP1 coverage estimated based on DTP3 coverage of ',DTP3,'. '],Explanation).
  412
  413% Estimate for DTP1 where DTP1 not reported: RMF
  414% -----------------------------------------
  415wuenic_I(C,dtp1,Y,'RMF',Explanation,RMFCoverage) :-
  416	wuenic_II(C,dtp3,Y,_,_,DTP3),
  417	not(wuenic_II(C,dtp1,Y,_,_,_)),
  418	not(workingGroupDecision(C,dtp1,Y,assignWUENIC,_,_,_)),
  419	RMFCoverage is round(18.258 - 0.6088*DTP3 - 0.0058*DTP3*DTP3),
  420	my_concat_atom(['Estimate based on DTP3 coverage of ',DTP3,'. '],Explanation).
  421
  422% ==================
  423% Estimate for RCV1 where RCV1 given at MCV1.
  424% -------------------------------------------
  425wuenic_I(C,rcv1,Y,Rule,'Estimate based on estimated MCV1. ',Coverage) :-
  426	estimate_required(C,rcv1,Y,_,FirstRubellaDose),
  427	wuenic_II(C,mcv1,Y,Rule,_Explanation,Coverage),
  428	not(workingGroupDecision(C,rcv1,Y,assignWUENIC,_,_,_)),
  429	not(firstRubellaAtSecondMCV(C,rcv1,Y,FirstRubellaDose)).
  430
  431% Estimate for RCV1 where RCV1 given at MCV2.
  432% -------------------------------------------
  433wuenic_I(C,rcv1,Y,Rule,'First dose of rubella vaccine given with second dose of measles containing vaccine. Estimate based on MCV2 estimate',Coverage) :-
  434	estimate_required(C,rcv1,Y,_,FirstRubellaDose),
  435	wuenic_II(C,mcv2,Y,Rule,_Explanation,Coverage),
  436	not(workingGroupDecision(C,rcv1,Y,assignWUENIC,_,_,_)),
  437	firstRubellaAtSecondMCV(C,rcv1,Y,FirstRubellaDose).
  438
  439% ===============
  440% Estimate assigned by working group.
  441% -----------------------------------
  442wuenic_I(C,V,Y,'W:',Explanation,Coverage) :-
  443	workingGroupDecision(C,V,Y,assignWUENIC,Explanation,_,Coverage).
  444
  445% First rubella given with second measles dose
  446% ---------------------------------------------
  447firstRubellaAtSecondMCV(_C, rcv1, _Y, mcv2).
  448
  449%   Preliminary estimates.
  450% ==============================================
  451
  452% No anchor points for any year. Reported data only.
  453% ------------------------------------------------
  454wuenic_II(C,V,Y,'R:',Explain,Coverage) :-
  455	estimate_required(C,V,Y,_,_),
  456	reported_time_series(C,V,Y,Source,Coverage),
  457	not(anchor_point(C,V,_,_,_,_)),
  458	explain(ro,Source,Explain).
  459
  460% Estimate at anchor points.
  461% --------------------------
  462wuenic_II(C,V,Y,Rule,Explanation,Coverage) :-
  463	estimate_required(C,V,Y,_,_),
  464	anchor_point(C,V,Y,Rule,Explanation,Coverage).
  465
  466% Estimate between anchor points: between two reported anchors, gov data.
  467% ----------------------------------------------
  468% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on coverage reported by national government. ',Coverage) :-
  469wuenic_II(C,V,Y,'R:','Estimate informed by reported data. ',Coverage) :-
  470	reported_time_series(C,V,Y,Source,Coverage),
  471	member(Source,['gov']),
  472	between_anchor_points(C,V,Y,_,RuleBefore,_,_,RuleAfter,_),
  473	both_anchors_resolved_to_reported(RuleBefore,RuleAfter),
  474	not(workingGroupDecision(C,V,Y,interpolate,_,_,_)),
  475	not(workingGroupDecision(C,V,Y,calibrate,_,_,_)).
  476
  477% Estimate between anchor points: between two reported anchors, admin data.
  478% -------------------------------------------------------------------------
  479% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on reported administrative data. ',Coverage) :-
  480wuenic_II(C,V,Y,'R:','Estimate informed by reported administrative data. ',Coverage) :-
  481	reported_time_series(C,V,Y,Source,Coverage),
  482	member(Source,['admin']),
  483	between_anchor_points(C,V,Y,_,RuleBefore,_,_,RuleAfter,_),
  484	both_anchors_resolved_to_reported(RuleBefore,RuleAfter),
  485	not(workingGroupDecision(C,V,Y,interpolate,_,_,_)),
  486	not(workingGroupDecision(C,V,Y,calibrate,_,_,_)).
  487
  488% Estimate between anchor points: between two reported anchors, interpolation.
  489% -----------------------------------------------------------------------------
  490wuenic_II(C,V,Y,'R:',Explain,Coverage) :-
  491	reported_time_series(C,V,Y,Source,Coverage),
  492	member(Source,['interpolated']),
  493	between_anchor_points(C,V,Y,_YrBefore,RuleBefore,_,_YrAfter,RuleAfter,_),
  494	both_anchors_resolved_to_reported(RuleBefore,RuleAfter),
  495	not(workingGroupDecision(C,V,Y,interpolate,_,_,_)),
  496	not(workingGroupDecision(C,V,Y,calibrate,_,_,_)),
  497	% DWB 2023-APR concat_atom(['Estimate based on interpolation between coverage reported by national government. '],Explain).
  498	my_concat_atom(['Estimate informed by interpolation between reported data. '],Explain).
  499
  500% Estimate between anchor points: calibrated
  501% ------------------------------------------
  502wuenic_II(C,V,Y,'C:',Explanation,Coverage) :-
  503	reported_time_series(C,V,Y,_,_ReportedCoverage), % MG: _ReportedCoverage because the line seems to be a guard
  504	between_anchor_points(C,V,Y,YrBefore,RuleBefore,_,YrAfter,RuleAfter,_),
  505	not(both_anchors_resolved_to_reported(RuleBefore,RuleAfter)),
  506	not(workingGroupDecision(C,V,Y,interpolate,_,_,_)),
  507	calibrate(C,V,YrBefore,YrAfter,Y,Coverage),
  508	my_concat_atom(['Reported data calibrated to ',YrBefore,' and ',YrAfter,' levels. '],Explanation).
  509
  510% Estimate between anchor points: interpolation forced by working group.
  511% ---------------------------------------------------------------------
  512wuenic_II(C,V,Y,'W-I:',Explanation,Coverage) :-
  513	between_anchor_points(C,V,Y,YrBefore,_,CoverageBefore,
  514							   YrAfter,_,CoverageAfter),
  515	workingGroupDecision(C,V,Y,interpolate,WGD_E,_,_),
  516	interpolate(YrBefore,CoverageBefore,YrAfter,CoverageAfter,Y,Coverage),
  517	% DWB 2023-APR concat_atom(['Estimate based on interpolation between ',YrBefore,' and ',YrAfter,' levels. ',WGD_E],Explanation).
  518	my_concat_atom(['Estimate informed by interpolation between ',YrBefore,' and ',YrAfter,' levels. ',WGD_E],Explanation).
  519
  520% Estimate before earliest anchor: reported
  521% ------------------------------------------
  522% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on reported data. ',ReportedCoverage) :-
  523wuenic_II(C,V,Y,'R:','Estimate informed by reported data. ',ReportedCoverage) :-
  524	reported_time_series(C,V,Y,Source,ReportedCoverage),
  525	member(Source,['admin','gov']),
  526	not(anchor_point(C,V,Y,_,_,_)),
  527
  528	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  529	Y < AnchorYear,
  530	not(anchor_point_before(C,V,AnchorYear)),
  531	member(AnchorRule,['R: AP']).
  532
  533% Estimate before earliest anchor: reported-extrapolated / interpolated
  534% ------------------------------------------
  535% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on interpolation between data reported by national government. ',ReportedCoverage) :-
  536wuenic_II(C,V,Y,'R:','Estimate informed by interpolation between reported data. ',ReportedCoverage) :-
  537	reported_time_series(C,V,Y,Source,ReportedCoverage),
  538	member(Source,['interpolated']),
  539	not(anchor_point(C,V,Y,_,_,_)),
  540
  541	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  542	Y < AnchorYear,
  543	not(anchor_point_before(C,V,AnchorYear)),
  544	member(AnchorRule,['R: AP']).
  545
  546% Estimate before earliest anchor: reported-extrapolated / interpolated
  547% ------------------------------------------
  548wuenic_II(C,V,Y,'R:','Estimate based on extrapolation from data reported by national government. ',ReportedCoverage) :-
  549	reported_time_series(C,V,Y,Source,ReportedCoverage),
  550	member(Source,['extrapolated']),
  551	not(anchor_point(C,V,Y,_,_,_)),
  552
  553	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  554	Y < AnchorYear,
  555	not(anchor_point_before(C,V,AnchorYear)),
  556	member(AnchorRule,['R: AP']).
  557
  558% Estimate before earliest anchor: calibrated
  559% ---------------------------------------------
  560wuenic_II(C,V,Y,'C:',Explanation,Coverage) :-
  561	reported_time_series(C,V,Y,_,ReportedCoverage),
  562	not(anchor_point(C,V,Y,_,_,_)),
  563
  564	anchor_point(C,V,AnchorYear,AnchorRule,_,AnchorCoverage),
  565	Y < AnchorYear,
  566	not(anchor_point_before(C,V,AnchorYear)),
  567	not(member(AnchorRule,['R: AP'])),
  568
  569	reported_time_series(C,V,AnchorYear,_,ReportedCoverageAtAnchor),
  570	Adj is AnchorCoverage - ReportedCoverageAtAnchor,
  571	Coverage is round(ReportedCoverage + Adj),
  572	my_concat_atom(['Reported data calibrated to ',AnchorYear,' levels. '],Explanation).
  573
  574% Estimate after latest anchor: reported
  575% --------------------------------------
  576% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on coverage reported by national government.',ReportedCoverage) :-
  577wuenic_II(C,V,Y,'R:','Estimate informed by reported data.',ReportedCoverage) :-
  578	reported_time_series(C,V,Y,gov,ReportedCoverage),
  579	not(anchor_point(C,V,Y,_,_,_)),
  580
  581	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  582	Y > AnchorYear,
  583	not(anchor_point_after(C,V,AnchorYear)),
  584	member(AnchorRule,['R: AP']).
  585
  586% Estimate after latest anchor: reported
  587% --------------------------------------
  588% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on reported administrative data. ',ReportedCoverage) :-
  589wuenic_II(C,V,Y,'R:','Estimate informed by reported administrative data. ',ReportedCoverage) :-
  590	reported_time_series(C,V,Y,admin,ReportedCoverage),
  591	not(anchor_point(C,V,Y,_,_,_)),
  592
  593	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  594	Y > AnchorYear,
  595	not(anchor_point_after(C,V,AnchorYear)),
  596	member(AnchorRule,['R: AP']).
  597
  598% Estimate after latest anchor: reported-extrapolated / interpolated
  599% ------------------------------------------------------------------
  600% DWB 2023-APR wuenic_II(C,V,Y,'R:','Estimate based on interpolation between data reported by national government. ',ReportedCoverage) :-
  601wuenic_II(C,V,Y,'R:','Estimate informed by interpolation between reported data. ',ReportedCoverage) :-
  602	reported_time_series(C,V,Y,Source,ReportedCoverage),
  603	member(Source,['interpolated']),
  604	not(anchor_point(C,V,Y,_,_,_)),
  605
  606	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  607	Y > AnchorYear,
  608	not(anchor_point_after(C,V,AnchorYear)),
  609	member(AnchorRule,['R: AP']).
  610
  611% Estimate after latest anchor: reported-extrapolated / interpolated
  612% ------------------------------------------------------------------
  613wuenic_II(C,V,Y,'R:','Estimate based on extrapolation from data reported by national government. ',ReportedCoverage) :-
  614	reported_time_series(C,V,Y,Source,ReportedCoverage),
  615	member(Source,['extrapolated']),
  616	not(anchor_point(C,V,Y,_,_,_)),
  617
  618	anchor_point(C,V,AnchorYear,AnchorRule,_,_),
  619	Y > AnchorYear,
  620	not(anchor_point_after(C,V,AnchorYear)),
  621	member(AnchorRule,['R: AP']).
  622
  623% Estimate after latest anchor: calibrated
  624% ----------------------------------------
  625wuenic_II(C,V,Y,'C:',Explanation,Coverage) :-
  626	reported_time_series(C,V,Y,_,ReportedCoverage),
  627	not(anchor_point(C,V,Y,_,_,_)),
  628
  629	anchor_point(C,V,AnchorYear,AnchorRule,_,AnchorCoverage),
  630	Y > AnchorYear,
  631	not(anchor_point_after(C,V,AnchorYear)),
  632	not(member(AnchorRule,['R: AP'])),
  633
  634	reported_time_series(C,V,AnchorYear,_,ReportedCoverageAtAnchor),
  635	%not(reported_reason_to_exclude(C,V,Y,_,_)),
  636	Adj is AnchorCoverage - ReportedCoverageAtAnchor,
  637	Coverage is round(ReportedCoverage + Adj),
  638	my_concat_atom(['Reported data calibrated to ',AnchorYear,' levels.'],Explanation).
  639
  640both_anchors_resolved_to_reported(RuleBefore,RuleAfter) :-
  641  member(RuleBefore,['R: AP']),
  642  member(RuleAfter,['R: AP']).
  643
  644% DWB 2023-APR	explain(ro,gov,'Estimate based on coverage reported by national government. ').
  645explain(ro,gov,'Estimate informed by reported data. ').
  646% DWB 2023-APR	explain(ro,admin,'Estimate based on reported administrative estimate. ').
  647explain(ro,admin,'Estimate informed by reported administrative data. ').
  648% DWB 2023-APR	explain(ro,interpolated,'Estimate based on interpolation between reported values. ').
  649explain(ro,interpolated,'Estimate informed by interpolation between reported data. ').
  650% DWB 2023-APR	explain(ro,extrapolated,'Estimate based on extrapolation from data reported by national government. ').
  651explain(ro,extrapolated,'Estimate informed by extrapolation from reported data. ').
  652
  653% =====================
  654% Level two processing:
  655%   Determine coverage value at anchor points
  656%   defined as years where there are multiple
  657%   data points (reported | survey | wgd).
  658%
  659%	Anchor point assignments:
  660%		Reported (gov,admin,extrapolated / interpolated timer series) supported by survey.
  661%		Reported not supported by survey.
  662%		Reported value "anchored" by working group.
  663%		Working group assigns anchor point value.
  664% =================================================
  665
  666% Anchor point set to reported value by working group
  667anchor_point(C, V, Y, Type, Explanation, Coverage),
  668    reported_time_series(C, V, Y, _, Reported),
  669    workingGroupDecision(C, V, Y, assignAnchor, Expl, _, Assigned),
  670    Reported = Assigned
  671 => Type = 'R: AP',
  672    Coverage = Assigned,
  673    Explanation = Expl.
  674
  675% Anchor point value NOT set to reported value by working group
  676anchor_point(C, V, Y, Type, Explanation, Coverage),
  677    reported_time_series(C, V, Y, _, Reported),
  678    workingGroupDecision(C, V, Y, assignAnchor, Expl, _, Assigned),
  679    Reported \= Assigned
  680 => Type = 'W: AP',
  681    Coverage = Assigned,
  682    my_concat_atom(
  683	['Estimate of ', Assigned, ' percent assigned by working group. ',
  684	 Expl], Explanation).
  685
  686% Survey results support reported
  687anchor_point(C, V, Y, Type, Explanation, Coverage),
  688    reported_time_series(C, V, Y, Source, Reported),
  689%   not(workingGroupDecision(C, V, Y, assignAnchor, _, _, _)),
  690    survey(C, V, Y, Expl, Survey),
  691    survey_reported_threshold(Threshold),
  692    abs(Reported - Survey) =< Threshold
  693 => Type = 'R: AP',
  694    Coverage = Reported,
  695    explain(ap, Source, Expl, Explanation).
  696
  697% Survey results challenge reported
  698anchor_point(C, V, Y, Type, Explanation, Coverage),
  699%   reported_time_series(C, V, Y, _, Reported),
  700%   not(workingGroupDecision(C,V,Y,assignAnchor,_,_,_)),
  701    survey(C, V, Y, Expl, Survey)
  702%   not(survey_supports_reported(ReportedCoverage,SurveyCoverage))
  703 => Type = 'S: AP',
  704    Coverage = Survey,
  705    my_concat_atom(
  706	['Survey evidence does not support reported data. Estimate based on survey results. ',
  707	 Expl, ' '], Explanation).
  708
  709anchor_point(_C, _V, _Y, _Type, _Explanation, _Coverage)
  710 => fail.
  711
  712explain(ap, gov, Expl, Explanation)
  713 => my_concat_atom(
  714	['Estimate informed by reported data supported by survey. ',
  715	 Expl], Explanation).
  716
  717explain(ap, admin, Expl, Explanation)
  718 => my_concat_atom(
  719	['Estimate informed by reported administrative data supported by survey. ',
  720	 Expl], Explanation).
  721
  722explain(ap, interpolated, Expl, Explanation)
  723 => my_concat_atom(
  724	['Estimate informed by interpolation between reported data supported by survey. ',
  725	 Expl], Explanation).
  726
  727explain(ap, extrapolated, Expl, Explanation)
  728 => my_concat_atom(
  729	['Estimate based on extrapolation from data reported by national government supported by survey. ',
  730	 Expl], Explanation).
  731
  732% ==============================================
  733% Level one processing:
  734%
  735%   Review data reported by national authorities
  736%   and survey results. Reported and survey data
  737%   are accepted, ignored, or modified. Modified
  738%   data are accepted in the modified form.
  739% ==============================================
  740
  741% Final survey information. If multiple survey in the
  742% same year, accepted and modified results are averaged.
  743% ----------------------------------
  744survey(C,V,Y,Explanation,Coverage) :-
  745%	bagof(Cov,survey_accepted(C,V,Y,_,_,Cov),CoverageList),
  746	bagof(Cov,Dist^SID^survey_accepted(C,V,Y,SID,Dist,Cov),CoverageList),
  747
  748	length(CoverageList,N),
  749	sum_list(CoverageList,SumCov),
  750	Coverage is round(SumCov / N),
  751	my_concat_atom(['Survey evidence of ',Coverage,' percent based on ',N, ' survey(s). '],Explanation).
  752
  753% Unmodified survey results accpeted.
  754% -----------------------------------
  755survey_accepted(C,V,Y,SurveyID,_,Coverage) :-
  756	survey_results_for_analysis(C,V,Y,SurveyID,_,Coverage),
  757	not(survey_results_modified(C,V,Y,SurveyID,_,_)),
  758	not(survey_reason_to_exclude(C,V,Y,SurveyID,_)).
  759
  760% Modified survey results accepted.
  761% ---------------------------------
  762survey_accepted(C,V,Y,SurveyID,_,ModifiedCoverage) :-
  763	survey_results_for_analysis(C,V,Y,SurveyID,_,_Coverage),
  764	survey_results_modified(C,V,Y,SurveyID,_,ModifiedCoverage),
  765	not(survey_reason_to_exclude(C,V,Y,SurveyID,_)).
  766
  767% Survey results modified for recall bias.
  768% ---------------------------------------
  769survey_results_modified(C,V,Y,SurveyID,Explanation,ModifiedCoverage) :-
  770	recall_bias_modified(C,V,Y,SurveyID,Explanation,ModifiedCoverage).
  771
  772% Adjust third dose for recall bias. Apply dropout observed
  773% between 1st and 3 dose documenented by card to history doses.
  774% Recalculate "card or history" based on adjustment to history.
  775% -------------------------------------------------------------
  776recall_bias_modified(C,V,Y,SurveyID,Explanation,ModifiedCoverage) :-
  777	member(V,['dtp3','pol3','hib3','hepb3','pcv3']),
  778
  779	% Associate third dose with first dose.
  780	vaccine(V,FirstDose),
  781
  782	% Third dose, card only
  783	survey_results(C,V,Y,SurveyID,DescriptionCard3Dose,C3Cov),
  784	member(confirm:MethodCard3Dose,DescriptionCard3Dose), member(MethodCard3Dose,['card']),
  785	member(age:AgeCohortCard3Dose,DescriptionCard3Dose),
  786	member(AgeCohortCard3Dose,['12-23 m','18-29 m','15-26 m', '24-35 m']),
  787
  788	% First dose, card or history
  789	survey_results(C,FirstDose,Y,SurveyID,DescriptionCoH1Dose,CoH1Cov),
  790	member(confirm:MethodCoH1Dose,DescriptionCoH1Dose), member(MethodCoH1Dose,['card or history']),
  791	member(age:AgeCohortCoH1,DescriptionCoH1Dose),
  792	member(AgeCohortCoH1,['12-23 m','18-29 m','15-26 m','24-35 m']),
  793
  794	% First dose, card only
  795	survey_results(C,FirstDose,Y,SurveyID,DescriptionCard1Dose,C1Cov), C1Cov > 0,
  796	member(confirm:MethodCard1Dose,DescriptionCard1Dose), member(MethodCard1Dose,['card']),
  797	member(age:AgeCohortCard1Dose,DescriptionCard1Dose),
  798	member(AgeCohortCard1Dose,['12-23 m','18-29 m','15-26 m','24-35 m']),
  799
  800	Adj is C3Cov / C1Cov,
  801	ThirdHistoryAdj is ((CoH1Cov - C1Cov)*Adj),
  802	CovAdjustedRecall is C3Cov + ThirdHistoryAdj,
  803
  804	bound_0_100(CovAdjustedRecall,ModifiedCoverage),
  805
  806	survey_results_for_analysis(C,V,Y,SurveyID,_,Coverage),
  807
  808	Coverage \= ModifiedCoverage,
  809
  810	SurveyCoverage is round(Coverage),
  811	CH1 is round(CoH1Cov),
  812	C1 is round(C1Cov),
  813	C3 is round(C3Cov),
  814
  815	survey_results_for_analysis(C,V,Y,SurveyID,Description,_),
  816	member(title:Survey,Description),
  817
  818	my_concat_atom([Survey,' card or history results of ',SurveyCoverage,' percent modifed for recall bias to ',
  819			ModifiedCoverage,' percent based on 1st dose card or history coverage of ',
  820			CH1,' percent, 1st dose card only coverage of ',C1,' percent and 3rd dose card only coverage of ',
  821			C3,' percent. '],Explanation).
  822
  823% -------------------------------------------------
  824% Survey results accepted if no reason to exclude.
  825%
  826% Reasons to exclude a survey include:
  827%    Sample size < 300,
  828%    The working group decides to exclude the survey.
  829
  830% --------------------------------------------------
  831% Reason to exclude survey: sample size < 300.
  832% --------------------------------------------
  833survey_reason_to_exclude(C,V,Y,SurveyID,Explanation) :-
  834	survey_results_for_analysis(C,V,Y,SurveyID,Description,_),
  835	member(ss:SampleSize,Description),
  836	SampleSize < 300,
  837	not(workingGroupDecision(C,V,Y,acceptSurvey,Explanation,SurveyID,_)),
  838	my_concat_atom(['Survey results ignored. Sample size ',SampleSize,' less than 300. '],Explanation).
  839
  840% Reason to exclude survey: working group decision to exclude survey identified by survey id.
  841% -------------------------------------------------------------------------------------------
  842survey_reason_to_exclude(C,V,Y,SurveyID,Explain) :-
  843	survey_results_for_analysis(C,V,Y,SurveyID,Description,_),
  844	workingGroupDecision(C,V,Y,ignoreSurvey,Explanation,SurveyID,_),
  845	member(title:Survey,Description),
  846	my_concat_atom([Survey,' results ignored by working group. ',Explanation],Explain).
  847
  848% Reason to exclude survey: working group decision to delate all surveys identified by country, vaccine, year.
  849% ------------------------------------------------------------------------------------------------------------
  850survey_reason_to_exclude(C,V,Y,SurveyID,Explain) :-
  851	survey_results_for_analysis(C,V,Y,SurveyID,Description,_),
  852	workingGroupDecision(C,V,Y,ignoreSurvey,Explanation,na,_),
  853	member(title:Survey,Description),
  854	my_concat_atom([Survey,' results ignored by working group. ',Explanation],Explain).
  855
  856% Survey results passed for inclusion in the analysis include:
  857% card or history results for cohorts 12-23, 18-29, 15-26, 24-35 months of age
  858% --------------------------------------------------------
  859survey_results_for_analysis(C,V,Y,SurveyID,Description,Coverage) :-
  860	survey_results(C,V,Y,SurveyID,Description,Coverage),
  861	member(confirm:Method,Description),
  862	member(Method,['card or history']),
  863	member(age:AgeCohort,Description),
  864	member(AgeCohort,['12-23 m','18-29 m','15-26 m','24-35 m']).
  865
  866% =============================================
  867% Create complete time series of reported data.
  868%
  869% Timeseries consists of:
  870%    coverage reported (gov or admin) by government
  871%    interpolated values if no reported or reported excluded between 2 years of reported.
  872%    value of nearest year for years no reported or reported excluded
  873%        before/after earliest/latest reported value (extrapolated)
  874% ==============================================
  875% Reported data.
  876% --------------
  877reported_time_series(C,V,Y,Source,Coverage) :-
  878		estimate_required(C,V,Y,_,_),
  879		reported(C,V,Y,Source,Coverage),
  880		not(reported_reason_to_exclude(C,V,Y,_,_)).
  881
  882% Interpolation, reported data excluded between two years
  883% -------------------------------------------------------
  884reported_time_series(C,V,Y,interpolated,Coverage) :-
  885		estimate_required(C,V,Y,_,_),
  886		reported(C,V,Y,_,_),
  887		once(reported_reason_to_exclude(C,V,Y,_,_Reason)),
  888		year_between_reported(C,V,Y,YearBefore,CoverageBefore,YearAfter,CoverageAfter),
  889		interpolate(YearBefore,CoverageBefore,YearAfter,CoverageAfter,Y,Coverage).
  890
  891% Interpolation, no reported data between two years
  892% -------------------------------------------------
  893reported_time_series(C,V,Y,interpolated,Coverage) :-
  894		estimate_required(C,V,Y,_,_),
  895		not(reported(C,V,Y,_,_)),
  896		year_between_reported(C,V,Y,YearBefore,CoverageBefore,YearAfter,CoverageAfter),
  897		interpolate(YearBefore,CoverageBefore,YearAfter,CoverageAfter,Y,Coverage).
  898
  899% Extrapolation, reported data for year beyond earliest / latest required estimate excluded
  900% ---------------------------------------------------------------------------------------
  901reported_time_series(C,V,Y,extrapolated,CoverageNearest) :-
  902		estimate_required(C,V,Y,_,_),
  903		reported(C,V,Y,_,_),
  904		once(reported_reason_to_exclude(C,V,Y,_,_Reason)),
  905		not(year_between_reported(C,V,Y,_,_,_,_)),
  906		nearest_reported(C,V,Y,_YearNearest,CoverageNearest).
  907
  908% Extrapolation, no reported data for year beyond earliest / latest required estimate
  909% ---------------------------------------------------------------------------------
  910reported_time_series(C,V,Y,extrapolated,CoverageNearest) :-
  911	estimate_required(C,V,Y,_,_),
  912	not(reported(C,V,Y,_,_)),
  913	not(year_between_reported(C,V,Y,_,_,_,_)),
  914	nearest_reported(C,V,Y,_YearNearest,CoverageNearest).
  915
  916% =====================================
  917% Reasons to exclude reported data are:
  918%  1. Working group decision.
  919%  2. Coverage > 100%
  920%  2. Inconsistent temporal changes (sawtooth or sudden change most recent year)
  921% =====================================
  922% MG, todo: short version that just checks if excluded
  923%
  924% Working group decision
  925reported_reason_to_exclude(C, V, Y, wdg, Explanation) :-
  926    workingGroupDecision(C, V, Y, ignoreReported, Expl, _, _),
  927    my_concat_atom(
  928	['Reported data excluded. ', Expl], Explanation).
  929
  930reported_reason_to_exclude(C, V, Y, Reason, Explanation) :-
  931    not(workingGroupDecision(C, V, Y, acceptReported, _, _, _)),
  932    reported(C, V, Y, _, Coverage),
  933    reported_reason_to_exclude(C, V, Y, Coverage, Reason, Explanation).
  934
  935% Reported coverage > 100%
  936reported_reason_to_exclude(_C, _V, _Y, Coverage, 100, Explanation) :-
  937    Coverage > 100,
  938    my_concat_atom(
  939	['Reported data excluded because ', Coverage,
  940	 ' percent greater than 100 percent. '], Explanation).
  941
  942% Sawtooth - inconsistent temporal change
  943reported_reason_to_exclude(C, V, Y, Coverage, sawtooth, Explanation) :-
  944    Before is Y - 1,
  945    reported(C, V, Before, _, CovBefore),
  946    After is Y + 1,
  947    reported(C, V, After, _, CovAfter),
  948    sawtooth_threshold(Threshold),
  949    Coverage - CovBefore > Threshold,
  950    Coverage - CovAfter  > Threshold,
  951    my_concat_atom(
  952	['Reported data excluded due to an increase from ', CovBefore,
  953	 ' percent to ', Coverage, ' percent with decrease ', CovAfter,
  954	 ' percent. '],	Explanation).
  955
  956reported_reason_to_exclude(C, V, Y, Coverage, sawtooth, Explanation) :-
  957    Before is Y - 1,
  958    reported(C, V, Before, _, CovBefore),
  959    After  is Y + 1,
  960    reported(C, V, After, _, CovAfter),
  961    sawtooth_threshold(Threshold),
  962    CovBefore - Coverage  > Threshold,
  963    CovAfter  - Coverage  > Threshold,
  964    my_concat_atom(
  965	['Reported data excluded due to decline in reported coverage from ',
  966	 CovBefore,' percent to ', Coverage, ' percent with increase to ',
  967	 CovAfter,' percent. '], Explanation).
  968
  969% Sudden decline in most recently reported data for new vaccines
  970reported_reason_to_exclude(C, V, Y, Coverage, temporalChange, Explanation) :-
  971    member(V, [pcv3, rotac]),
  972    not(reported_later(C, V, Y)),
  973    Before is Y - 1,
  974    reported(C, V, Before, _, CovBefore),
  975    sawtooth_threshold(Threshold),
  976    CovBefore - Coverage > Threshold,
  977    my_concat_atom(
  978	['Reported data excluded due to decline in reported coverage from ',
  979	 CovBefore,' level to ', Coverage,' percent. '], Explanation).
  980
  981% Sudden change in most recently reported data for classic vaccines
  982reported_reason_to_exclude(C, V, Y, Coverage, temporalChange, Explanation) :-
  983    not(member(V, [pcv3, rotac])),
  984    not(reported_later(C, V, Y)),
  985    Before is Y - 1,
  986    reported(C, V, Before, _, CovBefore),
  987    sawtooth_threshold(Threshold),
  988    abs(CovBefore - Coverage) > Threshold,
  989    my_concat_atom(
  990	['Reported data excluded due to sudden change in coverage from ',
  991	 CovBefore,' level to ', Coverage,' percent. '], Explanation).
  992
  993reported_later(C, V, Year) :-
  994    reported(C, V, After, _, _),
  995    After > Year.
  996
  997% Level 0:
  998% Reported to WHO and UNICEF is government estimate. If government
  999% estimate missing, then reported is administrative data.
 1000reported(C, V, Y, Source, Coverage),
 1001    gov(C, V, Y, Cov),
 1002    not(workingGroupDecision(C, V, Y, ignoreGov, _, _, _))
 1003 => Source = gov,
 1004    Coverage = Cov.
 1005
 1006reported(C, V, Y, Source, Coverage),
 1007    admin(C, V, Y, Cov),
 1008    not(workingGroupDecision(C, V, Y, ignoreAdmin, _, _, _))
 1009 => Source = admin,
 1010    Coverage = Cov.
 1011
 1012reported(_C, _V, _Y, _Source, _Coverage)
 1013 => fail.
 1014
 1015% ==========================================
 1016% Utilities and general auxiliary predicates
 1017% -------------------------------------------
 1018% Working group decision applies for a given year
 1019workingGroupDecision(C, V, Y, Action, Explanation, Survey, Coverage) :-
 1020    wgd(C, V, Begin, End, Action, Explanation, Survey, Coverage, _, _),
 1021    Begin =< Y, Y =< End.
 1022
 1023% Anchor points to given year
 1024% MG, issue: return closest anchor point instead of generate & test
 1025anchor_point_before(C, V, Before) :-
 1026    anchor_point(C, V, Y, _, _, _),
 1027    Y < Before.
 1028
 1029anchor_point_after(C, V, After) :-
 1030    anchor_point(C, V, Y, _, _, _),
 1031    Y > After.
 1032
 1033anchor_point_between(C, V, Before, After) :-
 1034    anchor_point(C, V, Y, _, _, _),
 1035    Before < Y, Y < After.
 1036
 1037% routines for interpolated points
 1038% --------------------------------
 1039year_between_reported(C,V,Y,YearBefore,CoverageBefore,YearAfter,CoverageAfter) :-
 1040	reported(C,V,YearBefore,_,CoverageBefore),
 1041	YearBefore < Y,
 1042	not(reported_reason_to_exclude(C,V,YearBefore,_,_)),
 1043
 1044	reported(C,V,YearAfter,_,CoverageAfter),
 1045	YearAfter > Y,
 1046	not(reported_reason_to_exclude(C,V,YearAfter,_,_)),
 1047	not(reported_data_between(C,V,YearBefore,YearAfter)).
 1048
 1049reported_data_between(C,V,YearBefore,YearAfter) :-
 1050	reported(C,V,YearBetween,_,_),
 1051	not(reported_reason_to_exclude(C,V,YearBetween,_,_)),
 1052	YearBetween > YearBefore,
 1053	YearBetween < YearAfter.
 1054
 1055% Extrapolation using nearest neighbor
 1056% ----------------------------------
 1057nearest_reported(C,V,Y,YearNearest,CoverageNearest) :-
 1058	reported(C,V,YearNearest,_,CoverageNearest),
 1059	not(reported_reason_to_exclude(C,V,YearNearest,_,_)),
 1060	not(reported_closer(C,V,Y,YearNearest)).
 1061
 1062reported_closer(C,V,Y,YearNearest) :-
 1063	reported(C,V,YearTest,_,_),
 1064	not(reported_reason_to_exclude(C,V,YearTest,_,_)),
 1065	abs(Y - YearTest) < abs(Y - YearNearest).
 1066
 1067% Get values of nearest surrounding anchor points.
 1068% -------------------------------------------------
 1069between_anchor_points(C,V,Y,PreceedingAnchorYear,PreceedingRule,PreceedingCov,
 1070                          SucceedingAnchorYear,SucceedingRule,SucceedingCov) :-
 1071	not(anchor_point(C,V,Y,_,_,_)),
 1072	anchor_point(C,V,PreceedingAnchorYear,PreceedingRule,_,PreceedingCov),
 1073	anchor_point(C,V,SucceedingAnchorYear,SucceedingRule,_,SucceedingCov),
 1074	PreceedingAnchorYear < Y,
 1075	SucceedingAnchorYear > Y,
 1076	not(anchor_point_between(C,V,PreceedingAnchorYear,SucceedingAnchorYear)).
 1077
 1078% Interpolate between two years
 1079interpolate(Early, EarlyCov, Late, LateCov, Year, Coverage) :-
 1080    Slope is (LateCov - EarlyCov) / (Late - Early),
 1081    Coverage is round(EarlyCov + (Year - Early) * Slope).
 1082
 1083% Calibrate reported data to anchor points
 1084% MG, issue: why are values rounded?
 1085calibrate(C, V, Prec, Succ, Y, Coverage),
 1086    reported_time_series(C, V, Prec, _, PrecRep),
 1087    reported_time_series(C, V, Succ, _, SuccRep)
 1088 => anchor_point(C, V, Prec, _, _, PrecCov),
 1089    anchor_point(C, V, Succ, _, _, SuccCov),
 1090    reported_time_series(C, V, Y, _, Reported),
 1091    interpolate(Prec, PrecRep, Succ, SuccRep, Y, RepInterp),
 1092    interpolate(Prec, PrecCov, Succ, SuccCov, Y, AnchInterp),
 1093    Adj is AnchInterp - RepInterp,
 1094    Coverage is round(Reported + Adj).
 1095
 1096% Reported data only at preceeding anchor
 1097calibrate(C, V, Prec, _, Y, Coverage)
 1098 => reported_time_series(C, V, Prec, _, PrecRep),
 1099    anchor_point(C, V, Prec, _, _, PrecCov),
 1100    Adj is PrecCov - PrecRep,
 1101    reported_time_series(C, V, Y, _, Reported),
 1102    Coverage is round(Reported + Adj).
 1103
 1104% Ensure estimates are between 0 and 99.
 1105% MG, issue: why are values rounded?
 1106bound_0_100(X, Y) :-
 1107    Y is max(0, min(99, round(X))).
 1108
 1109% Add underlying data to each C, V, Y estimate
 1110% MG, issue: is this extra iteration really needed?
 1111collect_data(C, V, Y, PrevRev, Admin, Gov, Reported, Vaccinated,
 1112        Target, UnpdBirths, UnpdSI, SeriesValue, Source, SurveyInfo)
 1113 => legacy_estimate(C, V, Y, PrevRev),
 1114    admin_data(C, V, Y, Admin),
 1115    gov_data(C, V, Y, Gov),
 1116    reported_data(C, V, Y, Reported),
 1117    vaccinated_data(C, V, Y, Vaccinated),
 1118    target_data(C, V, Y, Target),
 1119    time_series_data(C, V, Y, Source, SeriesValue),
 1120    unpd_births_data(C, Y, UnpdBirths),
 1121    unpd_si_data(C, Y, UnpdSI),
 1122    survey_data(C, V, Y, SurveyInfo).
 1123
 1124legacy_estimate(C, V, Y, PrevRev),
 1125    legacy(C, V, Y, D)
 1126 => PrevRev = D.
 1127
 1128legacy_estimate(_C, _V, _Y, PrevRev)
 1129 => PrevRev = ''.
 1130
 1131admin_data(C, V, Y, Admin),
 1132    admin(C, V, Y, D)
 1133 => Admin = D.
 1134
 1135admin_data(_C, _V, _Y, Admin)
 1136 => Admin = ''.
 1137
 1138gov_data(C, V, Y, Gov),
 1139    gov(C, V, Y, D)
 1140 => Gov = D.
 1141
 1142gov_data(_C, _V, _Y, Gov)
 1143 => Gov = ''.
 1144
 1145reported_data(C, V, Y, Reported),
 1146    reported(C, V, Y, _, D)
 1147 => Reported = D.
 1148
 1149reported_data(_C, _V, _Y, Reported)
 1150 => Reported = ''.
 1151
 1152vaccinated_data(C, V, Y, Vaccinated),
 1153    vaccinated(C, V, Y, D)
 1154 => Vaccinated = D.
 1155
 1156vaccinated_data(_C, _V, _Y, Vaccinated)
 1157 => Vaccinated = ''.
 1158
 1159target_data(C, V, Y, Target),
 1160    target(C, V, Y, D)
 1161 => Target = D.
 1162
 1163target_data(_C, _V, _Y, Target)
 1164 => Target = ''.
 1165
 1166unpd_births_data(C, Y, Births),
 1167    births_UNPD(C, Y, D)
 1168 => Births = D.
 1169
 1170unpd_births_data(_C, _Y, Births)
 1171 => Births = ''.
 1172
 1173unpd_si_data(C, Y, UnpdSI),
 1174    si_UNPD(C, Y, D)
 1175 => UnpdSI = D.
 1176
 1177unpd_si_data(_C, _Y, UnpdSI)
 1178 => UnpdSI = ''.
 1179
 1180time_series_data(C, V, Y, Source, Value),
 1181    reported_time_series(C, V, Y, S, D)
 1182 => Source = S,
 1183    Value = D.
 1184
 1185time_series_data(_C, _V, _Y, Source, Value)
 1186 => Source = '',
 1187    Value = ''.
 1188
 1189% MG: Check if this can be removed
 1190survey_data(C, V, Y, Info),
 1191    survey(C, V, Y, _, D)
 1192 => Info = D.
 1193
 1194survey_data(_C, _V, _Y, Info)
 1195 => Info = ''.
 1196
 1197% Collect explanations in natural language terms
 1198%
 1199% MG, todo: store information in the first place
 1200collect_explanations(C, V, Y, Explanations) :-
 1201    findall(Text, explanation(C, V, Y, Text), Explanations).
 1202
 1203collect_explanations(C, V, Y, Explanations) :-
 1204    findall(T, explanation(C, V, Y, T), Explanations).
 1205
 1206explanation(C, V, Y, Text) :-
 1207    survey_reason_to_exclude(C, V, Y, _, Text).
 1208
 1209explanation(C, V, Y, Text) :-
 1210    survey_results_modified(C, V, Y, _, Text, _).
 1211
 1212explanation(C, V, Y, Text) :-
 1213    reported_reason_to_exclude(C, V, Y, _, Text).
 1214
 1215explanation(C, V, Y, Text) :-
 1216    workingGroupDecision(C, V, Y, comment, Text, _, _).
 1217
 1218explanation(C, V, Y, Text) :-
 1219    workingGroupDecision(C, V, Y, acceptSurvey, Text, _, _).
 1220
 1221explanation(C, V, Y, Text) :-
 1222    workingGroupDecision(C, V, Y, acceptReported, Text, _, _).
 1223
 1224explanation(C, V, Y, Text) :-
 1225    workingGroupDecision(C, V, Y, ignoreGov, Text, _, _).
 1226
 1227% MG, issue: Why has this been commented out?
 1228% explanation(C, V, Y, Text) :-
 1229%     workingGroupDecision(C, V, Y, interpolate, Text, _, _).
 1230
 1231open_out_file(Out,File,Header) :-
 1232	open(File,write,Out),
 1233	write(Out,Header),
 1234	nl(Out).
 1235
 1236output_results([],_).
 1237output_results([H|T],Out) :- output_fields(H,Out), output_results(T,Out).
 1238
 1239output_fields([],Out) :- nl(Out).
 1240output_fields([H|T],Out) :- write(Out,H),write(Out,'\t'),output_fields(T,Out).
 1241
 1242% MG, temporary: concatenate string representation also of non-atoms
 1243my_concat_atom(List, String) :-
 1244    maplist(term_string, List, Strings),
 1245    atomics_to_string(Strings, String)