Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Emax decay function for use with vaccine interventions #370

Draft
wants to merge 18 commits into
base: main
Choose a base branch
from
Draft

Add Emax decay function for use with vaccine interventions #370

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
ae0d54b
add biphasic decay function
acavelan Mar 8, 2023
43871b2
rename DecayFuncHet into DecayFunctionHet
acavelan May 10, 2023
f9f44b1
rename tMult into timeFactorHet
acavelan May 10, 2023
a9aebd5
make ConstantDecayFunction inherit BaseHetDecayFunction
acavelan May 10, 2023
ba7d773
move eval() to DecayFunctionHet
acavelan May 22, 2023
6c9173b
cleanup
acavelan May 22, 2023
cecd665
move timeFactorHet to DecayFunction
acavelan May 22, 2023
51c1b12
move hetFactor to DecayFunction specializations
acavelan May 23, 2023
63238b3
change shared_ptr to unique_ptr
acavelan May 23, 2023
9ccbe59
replace DecayFunctionHet by DecayFunction
acavelan May 23, 2023
2784fc7
remove NullDecayFunction
acavelan May 23, 2023
f292cad
rename complement to increasing
acavelan May 23, 2023
540874c
add addition decay function
acavelan May 23, 2023
95bc6c1
add plus, minus, divides and multiplies decay functions
acavelan May 23, 2023
3d306a4
set minimum value t 0.0 in the new operator decay function
acavelan May 24, 2023
ba7f99b
fix example_scenario
acavelan May 24, 2023
a168552
add Emax decay function
acavelan Jun 19, 2023
c7f1877
rebase
acavelan Jan 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 64 additions & 2 deletions model/util/DecayFunction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,6 @@ class LinearDecayFunction : public DecayFunction {
{}

double compute(double effectiveAge) const{
// cout << "Linear " << effectiveAge << " " << invL << " " << hetFactor <<endl;
if( effectiveAge * invL * hetFactor < 1.0 ){
return 1.0 - effectiveAge * invL * hetFactor;
}else{
Expand Down Expand Up @@ -251,7 +250,7 @@ class OperatorDecayFunction : public DecayFunction {
DecayFunction(elt.getIncreasing(), elt.getInitialEfficacy(), elt.getCV()) {
const scnXml::DecayFunction::DecaySequence &decaySequence = elt.getDecay();
if(decaySequence.size() != 2)
throw util::xml_scenario_error("Operator decay function expects two decay functions, " + to_string(decaySequence.size()) +" were given.");
throw util::xml_scenario_error("The operator decay function expects two decay functions, but " + to_string(decaySequence.size()) +" were given.");

f1 = makeObject(decaySequence[0], "Operator::f1");
f2 = makeObject(decaySequence[1], "Operator::f2");
Expand Down Expand Up @@ -281,6 +280,67 @@ class OperatorDecayFunction : public DecayFunction {
T op;
};

class EmaxFunction : public DecayFunction {
public:
EmaxFunction( const scnXml::DecayFunction& elt ) :
DecayFunction(elt.getIncreasing(), elt.getInitialEfficacy(), elt.getCV()) {
const scnXml::DecayFunction::DecaySequence &decaySequence = elt.getDecay();
if(decaySequence.size() != 1)
throw util::xml_scenario_error("The Emax function expects exactly one user function, but " + to_string(decaySequence.size()) +" were given.");
if(!elt.getEmax().present())
throw util::xml_scenario_error("Emax function: the Emax parameter is not declared");
if(!elt.getIC50().present())
throw util::xml_scenario_error("Emax function: the IC50 parameter is not declared");
if(!elt.getSlope().present())
throw util::xml_scenario_error("Emax function: the Slope parameter is not declared");
if(!elt.getInitialConcentration().present())
throw util::xml_scenario_error("Emax function: the InitialConcentration parameter is not declared");

f = makeObject(decaySequence[0], "Emax::f");

Emax = elt.getEmax().get();
IC50 = elt.getIC50().get();
slope = elt.getSlope().get();
initialConcentration = elt.getInitialConcentration().get();

IC50 = pow(IC50, slope); // compute once
}

EmaxFunction(const EmaxFunction &copy, unique_ptr<DecayFunction> f) :
DecayFunction(copy),
f(move(f)),
Emax(copy.Emax),
IC50(copy.IC50),
slope(copy.slope),
initialConcentration(copy.initialConcentration) {}

double compute(double effectiveAge) const {
double fcomp = pow(initialConcentration * f->eval(effectiveAge), slope);
cout << "t: " << effectiveAge << endl;
cout << "AB(t): " << f->eval(effectiveAge) << endl;
cout << "rho: " << initialConcentration << endl;
cout << "slope: " << slope << endl;
cout << "IC50^slope: " << IC50 << endl;
cout << "rho * AB(t)^slope: " << fcomp << endl;
cout << "Emax(t) = " << Emax * fcomp / (fcomp + IC50) << endl;
return max(min(Emax * fcomp / (fcomp + IC50), 1.0), 0.0);
}

SimTime sampleAgeOfDecay (LocalRng& rng) const {
return sim::roundToTSFromDays( f->sampleAgeOfDecay(rng) );
}

unique_ptr<DecayFunction> hetSample(double hetFactor) const {
unique_ptr<DecayFunction> fhetSample = f->hetSample(hetFactor);
unique_ptr<EmaxFunction> copy = make_unique<EmaxFunction>(*this, move(fhetSample));
return move(copy);
}

private:
unique_ptr<DecayFunction> f;
double Emax, IC50, slope, initialConcentration;
};

// ----- interface / static functions -----

unique_ptr<DecayFunction> DecayFunction::makeObject(
Expand Down Expand Up @@ -310,6 +370,8 @@ unique_ptr<DecayFunction> DecayFunction::makeObject(
return unique_ptr<DecayFunction>(new OperatorDecayFunction<std::divides<double>>( elt ));
}else if( func == "multiplies" ){
return unique_ptr<DecayFunction>(new OperatorDecayFunction<std::multiplies<double>>( elt ));
}else if( func == "Emax" ){
return unique_ptr<DecayFunction>(new EmaxFunction( elt ));
}else{
throw util::xml_scenario_error("decay function type " + string(func) + " of " + string(eltName) + " unrecognized");
}
Expand Down
23 changes: 14 additions & 9 deletions schema/util.xsd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@ Licence: GNU General Public Licence version 2 or later (see COPYING) -->
<xs:enumeration value="minus"/>
<xs:enumeration value="divides"/>
<xs:enumeration value="multiplies"/>
<xs:enumeration value="Emax"/>
</xs:restriction>
</xs:simpleType>
</xs:attribute>
Expand Down Expand Up @@ -213,26 +214,30 @@ Licence: GNU General Public Licence version 2 or later (see COPYING) -->
</xs:documentation>
</xs:annotation>
</xs:attribute>
<xs:attribute name="rho" type="xs:double" use="optional">
<xs:attribute name="initialConcentration" type="xs:double" use="optional">
<xs:annotation>
<xs:documentation>
biphasic: Proportion between 0 and 1, proportion of the response that is short-lived.</xs:documentation>
</xs:annotation>
</xs:attribute>
<xs:attribute name="halflife_short" type="xs:string" use="optional">
<xs:attribute name="Emax" type="xs:double" use="optional">
<xs:annotation>
<xs:documentation>
biphasic: halflife of short lived component (default to years).

</xs:documentation>
maximum efficacy (between 0 and 1)</xs:documentation>
</xs:annotation>
</xs:attribute>
<xs:attribute name="halflife_long" type="xs:string" use="optional">
<xs:attribute name="slope" type="xs:double" use="optional">
<xs:annotation>
<xs:documentation>
biphasic: halflife of long lived component (default to years).

</xs:documentation>
steepness of the sigmoidal curve
</xs:documentation>
</xs:annotation>
</xs:attribute>
<xs:attribute name="IC50" type="xs:double" use="optional">
<xs:annotation>
<xs:documentation>
IC 50
</xs:documentation>
</xs:annotation>
</xs:attribute>
</xs:complexType>
Expand Down
Loading