Skip to content

Commit

Permalink
ORC v1.1.1: try catch in ucld scale operator
Browse files Browse the repository at this point in the history
  • Loading branch information
jordandouglas committed Oct 20, 2022
1 parent e6c3198 commit 1187ce0
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 79 deletions.
11 changes: 3 additions & 8 deletions build.xml
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
</description>

<!-- set global properties for this build -->
<property name="versionORC" value="1.1.0" />
<property name="versionORC" value="1.1.1" />
<property name="srcORC" location="src" />
<property name="docORC" location="doc" />
<property name="buildORC" location="build" />
Expand Down Expand Up @@ -189,7 +189,6 @@
<mkdir dir="${Add_on_dir}/lib" />
<mkdir dir="${Add_on_dir}/doc" />
<mkdir dir="${Add_on_dir}/examples" />
<mkdir dir="${Add_on_dir}/templates" />
<mkdir dir="${Add_on_dir}/fxtemplates" />

<copy todir="${Add_on_dir}">
Expand All @@ -211,10 +210,6 @@
</copy>


<copy todir="${Add_on_dir}/templates">
<fileset file="templates/*.xml" />
</copy>


<copy todir="${Add_on_dir}/fxtemplates">
<fileset file="fxtemplates/*.xml" />
Expand Down Expand Up @@ -279,8 +274,8 @@
productVersion="${version_number}.0"
txtProductVersion="${versionORC}" />

<copy todir="${Windows_package_dirORC}/templates/">
<fileset dir="templates/" />
<copy todir="${Windows_package_dirORC}/fxtemplates/">
<fileset dir="fxtemplates/" />
</copy>

<zip destfile="${Windows_dirORC}/${CladeAge_name} v${versionORC}.zip">
Expand Down
153 changes: 83 additions & 70 deletions src/orc/consoperators/UcldScalerOperator.java
Original file line number Diff line number Diff line change
Expand Up @@ -80,29 +80,35 @@ public double proposal() {
LN = (LogNormalDistributionModel) distribution;
}

// Step1: calculate rates/quantiles for real rates given the current ucldStdev
if (this.usingRates) {

// Rates -> quantiles
for (int idx = 0; idx < this.nrates; idx++) {
double r = rates.getValue(idx);
real_rates[idx] = r;
double q = getRateQuantiles(r);
quantiles[idx] = q;

try {

// Step1: calculate rates/quantiles for real rates given the current ucldStdev
if (this.usingRates) {

// Rates -> quantiles
for (int idx = 0; idx < this.nrates; idx++) {
double r = rates.getValue(idx);
real_rates[idx] = r;
double q = getRateQuantiles(r);
quantiles[idx] = q;

}

}else {

// Quantiles -> rates
for (int idx = 0; idx < this.nrates; idx++) {
double q = this.quantilesParam.getValue(idx);
quantiles[idx] = q;
double r = getRealRate(q);
real_rates[idx] = r;

}

}

}else {

// Quantiles -> rates
for (int idx = 0; idx < this.nrates; idx++) {
double q = this.quantilesParam.getValue(idx);
quantiles[idx] = q;
double r = getRealRate(q);
real_rates[idx] = r;

}

}catch(Exception e) {
return Double.NEGATIVE_INFINITY;
}

// Step2: use scale operation to propose a new_ucldStdev
Expand All @@ -122,15 +128,11 @@ public double proposal() {
new_stdev = stdev * scale;
if (new_stdev <= 0) return Double.NEGATIVE_INFINITY;

// set the new value
ucldStdev.setValue(new_stdev);


// return the log hastings ratio for scale operation





// use mean (M=1) and variance (stdev square) to calculate miu and sigma square of lognormal in log space
double variance = FastMath.exp(stdev * stdev) -1;
double new_variance = FastMath.exp(new_stdev * new_stdev) -1; // new sigma square of lognormal
Expand All @@ -139,54 +141,65 @@ public double proposal() {

// Step3: calculate the new rates or quantiles under the proposed new_ucldStdev

// Rates
if (this.usingRates) {
for (int idx = 0; idx < this.nrates; idx++) {
double r_ = getRealRate(quantiles[idx]);

// to avoid numerical issues
if (r_== 0.0 || r_ == Double.POSITIVE_INFINITY) {
return Double.NEGATIVE_INFINITY;
}

// partial derivative of icdf_s'(cdf_s(r)) / r
hastingsRatio = hastingsRatio + getDicdf(real_rates[idx],quantiles[idx],miu,new_miu,stdev,new_stdev);

try {

// Rates
if (this.usingRates) {
for (int idx = 0; idx < this.nrates; idx++) {
double r_ = getRealRate(quantiles[idx]);

// to avoid numerical issues
if (r_== 0.0 || r_ == Double.POSITIVE_INFINITY) {
return Double.NEGATIVE_INFINITY;
}

// partial derivative of icdf_s'(cdf_s(r)) / r
hastingsRatio = hastingsRatio + getDicdf(real_rates[idx],quantiles[idx],miu,new_miu,stdev,new_stdev);

rates.setValue(idx, r_);
}

// Quantiles
}else {

// Step3: calculate the new quantiles such that the rates remain constant

for (int idx = 0; idx < this.nrates; idx++) {

// Original quantile
double q = this.quantilesParam.getArrayValue(idx);

// Original rate (want the new rate to be the same)
double r = real_rates[idx];

// Propose q_ such that the rate is constant
double q_ = this.getRateQuantiles(r);
if (q_ <= 0 || q_ >= 1) return Double.NEGATIVE_INFINITY;
quantilesParam.setValue(idx, q_);


// logHR
double logHRBranch = this.getQuantHR(q, miu, new_miu, stdev, new_stdev);


//System.out.println(stdev + " -> " + new_stdev + "; " + q + " -> " + q_ + "; HR: " + logHRBranch);
hastingsRatio += logHRBranch;

}


rates.setValue(idx, r_);
}

// Quantiles
}else {

// Step3: calculate the new quantiles such that the rates remain constant

for (int idx = 0; idx < this.nrates; idx++) {

// Original quantile
double q = this.quantilesParam.getArrayValue(idx);

// Original rate (want the new rate to be the same)
double r = real_rates[idx];

// Propose q_ such that the rate is constant
double q_ = this.getRateQuantiles(r);
if (q_ <= 0 || q_ >= 1) return Double.NEGATIVE_INFINITY;
quantilesParam.setValue(idx, q_);


// logHR
double logHRBranch = this.getQuantHR(q, miu, new_miu, stdev, new_stdev);


//System.out.println(stdev + " -> " + new_stdev + "; " + q + " -> " + q_ + "; HR: " + logHRBranch);
hastingsRatio += logHRBranch;

}



}

}catch(Exception e) {
return Double.NEGATIVE_INFINITY;
}


// set the new value
ucldStdev.setValue(new_stdev);

//System.out.println("HR " + hastingsRatio);

return hastingsRatio;
Expand Down
2 changes: 1 addition & 1 deletion version.xml
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
<package name='ORC' version='1.1.0'>
<package name='ORC' version='1.1.1'>
<depends on='BEAST.base' atleast='2.7.0'/>
<depends on='BEAST.app' atleast='2.7.0'/>
<depends on='BEASTLabs' atleast='2.0.0'/>
Expand Down

0 comments on commit 1187ce0

Please sign in to comment.