if (expDomainIcon == null)
{
throw new RuntimeException("The icon for the experiment domain "+expdnames[0]+" is null.");
}
SLIconifiedAnnotation annot = new SLIconifiedAnnotation(null,null);
annot.addAnnotation(buffDomainIcon);
annot.addAnnotation("."+buffdnames[1]);
annot.addAnnotation("*");
annot.addAnnotation(expDomainIcon);
annot.addAnnotation("."+expdnames[1]);
annot.addAnnotation("*");
annot.addAnnotation(ssid);
annot.setText(buffername+"*"+exp+"*"+ssid);
legendStrings.add(annot);
//legendStrings.add(buffername+"*"+exp+"*"+ssid);
//stop = System.currentTimeMillis();
//System.out.println("Time it took3="+(stop-start));
//start = stop;
OBSNUM = this.transientStorage.getTitrationObservationNumbers(exp,ssid);
//Need the initial volume for this titration in order to
DblMatrix initialVolume = new DblMatrix(new Double(this.transientStorage.getInitialVolume(exp,ssid)));
//param = this.mkParameterName(new String[]{"initialVolExp",buffername,exp});
//initialpH_mean = initialpH.plus(initialpHExpdev).plus(initialpHExpSSdev);
DblMatrix currentVolume = initialVolume.plus(0.0);
DblMatrix DblCbThatShouldBeZero = new DblMatrix(OBSNUM.size());
DblMatrix DblpH4CbThatShouldBeZero = new DblMatrix(OBSNUM.size());
for (int k=0;k<OBSNUM.size();k++)
{
//long stop1 = System.currentTimeMillis();
//long start1 = stop1;
obsnum = (String)OBSNUM.get(k);
//Parameters not used should be set to reasonable default values (Usually 0.0) by now.
if (this.adjustForLatentpHLink())
{
String val = this.transientStorage.getpH(exp,ssid,obsnum);
latentpHi = new DblMatrix(new Double(val));
val = this.referenceStorage.getpH(exp,ssid,obsnum);
DblMatrix obspH = new DblMatrix(new Double(val));
//Implement the latent link equation...
DblMatrix B0 = interceptpH.plus(interceptpHExp).plus(interceptpHExpSS);
DblMatrix B1 = linearpH.plus(linearpHExp).plus(linearpHExpSS);
DblMatrix B2 = quadraticpH.plus(quadraticpHExp).plus(quadraticpHExpSS);
DblMatrix B3 = cubicpH.plus(cubicpHExp).plus(cubicpHExpSS);
//DblMatrix centeredDiff = latentpHi.minus(B0);
DblMatrix centeredDiff = latentpHi;
DblMatrix centeredDiffSqrd = centeredDiff.times(centeredDiff);
DblMatrix centeredDiffCubed = centeredDiffSqrd.times(centeredDiff);
DblMatrix obspH_hat = B0.plus(B1.times(centeredDiff).plus(B2.times(centeredDiffSqrd)).plus(B3.times(centeredDiffCubed)));
DblMatrix resid = obspH.minus(obspH_hat);
pHresids.setDblAt(resid,like_index);
this.normal.setMean(DblMatrix.ZERO);
this.normal.setStd(sigmaCalib);
DblMatrix likeval = this.normal.pdf(resid);
if (DblMatrix.test(DblMatrix.isNaN(likeval)))
{
resid.show("residual");
latentpHi.transpose().show("Latent pH");
obspH.transpose().show("Observed pH");
sigmaCalib.show("sigmaCalib");
LIKE.replicate(DblMatrix.INF,LIKE.Size);
//System.out.println("NaN in pH likelihood calc.");
return(LIKE);
}
if (DblMatrix.test(likeval.eq(0.0)))
{
resid.show("residual");
latentpHi.transpose().show("Latent pH");
obspH.transpose().show("Observed pH");
obspH_hat.transpose().show("Predicted pH");
sigmaCalib.show("sigmaCalib");
likeval.show("likeval","0.00E00");
LIKE.replicate(DblMatrix.INF,LIKE.Size);
//System.out.println("Zero in pH likelihood calc.");
return(LIKE);
}
LIKE.setDblAt(likeval,like_index);
like_index++;
}
//stop1 = System.currentTimeMillis();
//System.out.println("Time it took1a="+(stop1-start1));
//start1 = stop1;
//Bring the RhsSol up to date with the latest addition...
String titrant = this.transientStorage.getTitrant(exp,ssid,obsnum);
DblMatrix volAdded = new DblMatrix(new Double(this.transientStorage.getVolumeAdded(exp,ssid,obsnum)));
DblMatrix conc = null;
if (((Boolean)this.adjustForMap.get("adjustForTitrantConc")).booleanValue())
{
//Consider this the dilution of a unit concentration of titrant where the
//unit concentration is given by the latent "C" parameter for this titrant.
DblMatrix unitConcentration = new DblMatrix(1.0);
DblMatrix unitConcentrationExpdev = new DblMatrix(0.0);
DblMatrix unitConcentrationExpSSdev = new DblMatrix(0.0);
param = this.mkParameterName(new String[]{"T",titrant,buffername});
if (this.hasParameter(param))
{
unitConcentration = this.getParam(param);
}
param = this.mkParameterName(new String[]{"TExpdev",titrant,buffername,exp});
if (this.hasParameter(param))
{
unitConcentrationExpdev = this.getParam(param);
}
param = this.mkParameterName(new String[]{"TExpSSdev",titrant,buffername,exp,ssid});
if (this.hasParameter(param))
{
unitConcentrationExpSSdev = this.getParam(param);
}
conc = unitConcentration.plus(unitConcentrationExpdev).plus(unitConcentrationExpSSdev);
DblMatrix unitFactor = new DblMatrix(new Double(this.transientStorage.getTitrantConcentration(exp,ssid,obsnum)));
conc = conc.times(unitFactor);
}
else //Just take the recorded concentration as given...
{
conc = new DblMatrix(new Double(this.transientStorage.getTitrantConcentration(exp,ssid,obsnum)));
}
if (!this.titrantMap.containsKey(titrant))
{
throw new RuntimeException("Titrant map does not contain a mapping for the titrant "+titrant+".");
}
BufferSolution titr = (BufferSolution)this.titrantMap.get(titrant);
if (titr == null)
{
throw new RuntimeException("Titrant map contains a null mapping for the titrant "+titrant+".");
}
titr.setConcentration(titrant,conc);
//titr.show("titrant----------------");
//Update the RhsSol for the latest addition of titrant.
//System.out.println("Here is the repository of RhsSol:"+RhsSol.getRepository());
//stop1 = System.currentTimeMillis();
//System.out.println("Time it took2a="+(stop1-start1));
//start1 = stop1;
RhsSol = BufferSolution.mixSol(RhsSol,currentVolume,titr,volAdded);
//stop1 = System.currentTimeMillis();
//System.out.println("Time it took3a="+(stop1-start1));
//start1 = stop1;
//RhsSol.show("RhsSol after mix");
//Remember to update the total volume after an addition of titrant.
currentVolume = currentVolume.plus(volAdded);
if (this.isSamplingObject("RandomBufferSolution"))
{
if (((Boolean)this.adjustForMap.get("adjustForpKa")).booleanValue())
{
String[] components = MarginalSol.getComponents();
ComplexBufferPredictor cmplxPred = MarginalSol.getSummaryBuffer();
DblMatrix[] pKa = cmplxPred.getEstimatedpKa(buffername);
DblMatrix[] pKb = cmplxPred.getEstimatedpKb(buffername);
String ionstr = null;
String groupstr = null;
for (int ion=0;ion<pKa.length;ion++)
{
Integer ION = new Integer(ion);
ionstr = ION.toString();
if (!pKa[ion].isEmpty())
{
for (int group=0;group<pKa[ion].getN();group++)
{
if (group == 0)
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix meanpKa = pKa[ion].getDblAt(group);
DblMatrix pKaExpdev = new DblMatrix(0.0);
DblMatrix pKaExpSSdev = new DblMatrix(0.0);
if (!((Boolean)this.marginalizeMap.get("adjustForpKa")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKa",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKa = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKaExpdev",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKaExpdev = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKaExpSSdev",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKaExpSSdev = this.getParam(param);
}
}
pKa[ion].setDblAt(meanpKa.plus(pKaExpdev).plus(pKaExpSSdev),group);
}
else
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix pKaToLeft = pKa[ion].getDblAt(group-1);
DblMatrix meanpKadiff = pKa[ion].getDblAt(group).minus(pKaToLeft);
DblMatrix pKadiffExpdev = new DblMatrix(0.0);
DblMatrix pKadiffExpSSdev = new DblMatrix(0.0);
if (!((Boolean)this.marginalizeMap.get("adjustForpKa")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKadiff",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKadiff = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKadiffExpdev",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKadiffExpdev = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKadiffExpSSdev",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKadiffExpSSdev = this.getParam(param);
}
}
DblMatrix pKaDiff = meanpKadiff.plus(pKadiffExpdev).plus(pKadiffExpSSdev);
pKa[ion].setDblAt(pKaToLeft.plus(pKaDiff),group);
}
}
}
if (!pKb[ion].isEmpty())
{
for (int group=0;group<pKb[ion].getN();group++)
{
if (group == 0)
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix meanpKb = pKb[ion].getDblAt(group);
DblMatrix pKbExpdev = new DblMatrix(0.0);
DblMatrix pKbExpSSdev = new DblMatrix(0.0);
if (!((Boolean)this.marginalizeMap.get("adjustForpKa")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKb",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKb = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKb",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbExpdev = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKb",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbExpSSdev = this.getParam(param);
}
}
pKb[ion].setDblAt(meanpKb.plus(pKbExpdev).plus(pKbExpSSdev),group);
}
else
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix pKbToLeft = pKb[ion].getDblAt(group-1);
DblMatrix meanpKbdiff = pKb[ion].getDblAt(group).minus(pKbToLeft);
DblMatrix pKbdiffExpdev = new DblMatrix(0.0);
DblMatrix pKbdiffExpSSdev = new DblMatrix(0.0);
if (!((Boolean)this.marginalizeMap.get("adjustForpKa")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKbdiff",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKbdiff = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKbdiffExpdev",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbdiffExpdev = this.getParam(param);
}
}
if (!((Boolean)this.marginalizeMap.get("adjustForpKaExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"pKbdiffExpSSdev",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbdiffExpSSdev = this.getParam(param);
}
}
DblMatrix pKbDiff = meanpKbdiff.plus(pKbdiffExpdev).plus(pKbdiffExpSSdev);
pKb[ion].setDblAt(pKbToLeft.plus(pKbDiff),group);
}
}
}
}
//Notice that we will always get the last one -- we only can generate a single realization
//as a random sample...
//Update the pKa for this component....
if (((Boolean)this.adjustForMap.get("adjustForpKa")).booleanValue())
{
if (pKa.length>0)
{
// for (int w=0;w<pKa.length;w++)
// {
// pKa[w].show("Latest estimated pKa for MarginalSol "+w,"0.0000E00");
// }
cmplxPred.setEstimatedpKa(buffername,pKa);
cmplxPred.setEstimatedpKb(buffername,pKb);
}
}
MarginalSol.setSummaryBuffer(cmplxPred);
}
}//Finished updating pKa/pKb for marginal mean buffer random sample...
//Update the pKa/pKb for all components...
//We take care to maintain order-restricted inference on the pKa/pKb
//even in the presence of experiment and subsample effects.
//
// String[] components = RhsSol.getComponents();
// ComplexBufferPredictor cmplxPred = RhsSol.getSummaryBuffer();
// DblMatrix[] pKa = cmplxPred.getEstimatedpKa(buffername);
// DblMatrix[] pKb = cmplxPred.getEstimatedpKb(buffername);
//
String[] components = RhsSol.getComponents();
ComplexBufferPredictor cmplxPred = RhsSol.getSummaryBuffer();
DblMatrix[] pKa = cmplxPred.getEstimatedpKa(buffername);
DblMatrix[] pKb = cmplxPred.getEstimatedpKb(buffername);
if (pKa != null)
{
if (((Boolean)this.adjustForMap.get("adjustForpKa")).booleanValue())
{
String ionstr = null;
String groupstr = null;
for (int ion=0;ion<pKa.length;ion++)
{
Integer ION = new Integer(ion);
ionstr = ION.toString();
if (!pKa[ion].isEmpty())
{
for (int group=0;group<pKa[ion].getN();group++)
{
if (group == 0)
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix meanpKa = pKa[ion].getDblAt(group);
DblMatrix pKaExpdev = new DblMatrix(0.0);
DblMatrix pKaExpSSdev = new DblMatrix(0.0);
param = this.mkParameterName(new String[]{"pKa",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKa = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKaExpdev",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKaExpdev = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKaExpSSdev",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKaExpSSdev = this.getParam(param);
}
pKa[ion].setDblAt(meanpKa.plus(pKaExpdev).plus(pKaExpSSdev),group);
}
else
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix pKaToLeft = pKa[ion].getDblAt(group-1);
DblMatrix meanpKadiff = pKa[ion].getDblAt(group).minus(pKaToLeft);
DblMatrix pKadiffExpdev = new DblMatrix(0.0);
DblMatrix pKadiffExpSSdev = new DblMatrix(0.0);
param = this.mkParameterName(new String[]{"pKadiff",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKadiff = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKadiffExpdev",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKadiffExpdev = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKadiffExpSSdev",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKadiffExpSSdev = this.getParam(param);
}
DblMatrix pKaDiff = meanpKadiff.plus(pKadiffExpdev).plus(pKadiffExpSSdev);
pKa[ion].setDblAt(pKaToLeft.plus(pKaDiff),group);
}
}
}
if (!pKb[ion].isEmpty())
{
for (int group=0;group<pKb[ion].getN();group++)
{
if (group == 0)
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix meanpKb = pKb[ion].getDblAt(group);
DblMatrix pKbExpdev = new DblMatrix(0.0);
DblMatrix pKbExpSSdev = new DblMatrix(0.0);
param = this.mkParameterName(new String[]{"pKb",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKb = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKb",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbExpdev = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKb",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbExpSSdev = this.getParam(param);
}
pKb[ion].setDblAt(meanpKb.plus(pKbExpdev).plus(pKbExpSSdev),group);
}
else
{
Integer GROUP = new Integer(group);
groupstr = GROUP.toString();
DblMatrix pKbToLeft = pKb[ion].getDblAt(group-1);
DblMatrix meanpKbdiff = pKb[ion].getDblAt(group).minus(pKbToLeft);
DblMatrix pKbdiffExpdev = new DblMatrix(0.0);
DblMatrix pKbdiffExpSSdev = new DblMatrix(0.0);
param = this.mkParameterName(new String[]{"pKbdiff",buffername,ionstr,groupstr});
if (this.hasParameter(param))
{
meanpKbdiff = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKbdiffExpdev",buffername,exp,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbdiffExpdev = this.getParam(param);
}
param = this.mkParameterName(new String[]{"pKbdiffExpSSdev",buffername,exp,ssid,ionstr,groupstr});
if (this.hasParameter(param))
{
pKbdiffExpSSdev = this.getParam(param);
}
DblMatrix pKbDiff = meanpKbdiff.plus(pKbdiffExpdev).plus(pKbdiffExpSSdev);
pKb[ion].setDblAt(pKbToLeft.plus(pKbDiff),group);
}
}
}
}
}
}
//stop1 = System.currentTimeMillis();
//System.out.println("Time it took4a="+(stop1-start1));
//start1 = stop1;
//Update the pKa for this component....
if (pKa != null)
{
if (pKa.length>0)
{
//pKa[3].show("Latest estimated pKa...");
cmplxPred.setEstimatedpKa(buffername,pKa);
cmplxPred.setEstimatedpKb(buffername,pKb);
}
}
//RhsSol.show("RhsSol??");
//Now we hope that the Cb of the right hand side solution when
//evaluated at the latent pH is exactly equal to zero!
HashMap Cb = null;
if (this.adjustForLatentpHLink())
{
Cb = RhsSol.Cb(latentpHi,true);
DblpH4CbThatShouldBeZero.setDblAt(latentpHi,k);
if (this.adjustForSmoothingComponent())
{
SmoothingpH.setDblAt(latentpHi,smoothingDataIndex);
//if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
//{
// SmoothingpH.setDblAt(latentpHi,smoothingDataIndex);
//}
}
}
else
{
String val = this.transientStorage.getpH(exp,ssid,obsnum);
DblMatrix observedpHi = new DblMatrix(new Double(val));
Cb = RhsSol.Cb(observedpHi,true);
DblpH4CbThatShouldBeZero.setDblAt(observedpHi,k);
//Setting the XData for the smoother...
if (this.adjustForSmoothingComponent())
{
SmoothingpH.setDblAt(observedpHi,smoothingDataIndex);
//if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
//{
// SmoothingpH.setDblAt(observedpHi,smoothingDataIndex);
//}
}
}
//stop1 = System.currentTimeMillis();
//System.out.println("Time it took5a="+(stop1-start1));
//start1 = stop1;
// Set keys = Cb.keySet();
// Iterator iter = keys.iterator();
// while (iter.hasNext())
// {
// System.out.println(iter.next());
// }
// System.out.println("buffername="+buffername);
// ///
//DblMatrix thisUntitledCb = (DblMatrix)Cb.get(buffername);
//thisUntitledCb.show("thisUntitledCb","0.00E0");
///
//DblMatrix thisTitrantCb = (DblMatrix)Cb.get("com.mockturtlesolutions.HCl");
//thisTitrantCb.show("titrantCb","0.00E0");
//DblMatrix thiswaterCb = (DblMatrix)Cb.get("com.mockturtlesolutions.water");
//thiswaterCb.show("waterCb","0.00E0");
DblMatrix thisPointCb = (DblMatrix)Cb.get("total");
DblMatrix P0 = new DblMatrix(0.0);
DblMatrix P1 = new DblMatrix(0.0);
DblMatrix P2 = new DblMatrix(0.0);
DblMatrix P3 = new DblMatrix(0.0);
//Set the YData for the next go-round.
if (this.adjustForSmoothingComponent())
{
if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
{
SmoothingCbResidual.setDblAt(thisPointCb,smoothingDataIndex);
SmoothingSmoother.useSufficient(false);
DblMatrix complexCbContribution = SmoothingSmoother.predict(new DblMatrix[]{SmoothingpH.getDblAt(smoothingDataIndex)});
thisPointCb = thisPointCb.minus(complexCbContribution);
//System.out.println("Done with one!");
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.POLYNOMIAL_SMOOTHING)
{
DblMatrix pHcenter = SmoothingpH.getDblAt(smoothingDataIndex).minus(7.0);//Center about neutral pH as convention.
DblMatrix pHcenterSq = pHcenter.times(pHcenter);
DblMatrix pHcenterCd = pHcenter.times(pHcenterSq);
if (((Integer)this.adjustForMap.get("adjustForSmoothingModel")).intValue() == DefaultComplexBufferPredictorGenerator.LINEAR_SMOOTHER)
{
param = this.mkParameterName(new String[]{"smoothingIntercept",buffername});
P0 = this.getParam(param);
param = this.mkParameterName(new String[]{"smoothingSlope",buffername});
P1 = this.getParam(param);
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingModel")).intValue() == DefaultComplexBufferPredictorGenerator.QUADRATIC_SMOOTHER)
{
param = this.mkParameterName(new String[]{"smoothingIntercept",buffername});
P0 = this.getParam(param);
param = this.mkParameterName(new String[]{"smoothingSlope",buffername});
P1 = this.getParam(param);
param = this.mkParameterName(new String[]{"smoothingQuadratic",buffername});
P2 = this.getParam(param);
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingModel")).intValue() == DefaultComplexBufferPredictorGenerator.CUBIC_SMOOTHER)
{
param = this.mkParameterName(new String[]{"smoothingIntercept",buffername});
P0 = this.getParam(param);
param = this.mkParameterName(new String[]{"smoothingSlope",buffername});
P1 = this.getParam(param);
param = this.mkParameterName(new String[]{"smoothingQuadratic",buffername});
P2 = this.getParam(param);
param = this.mkParameterName(new String[]{"smoothingCubic",buffername});
P3 = this.getParam(param);
}
// P0.show("P0");
// P1.show("P1");
// P2.show("P2");
// P3.show("P3");
DblMatrix complexCbContribution = P0.plus(P1.times(pHcenter)).plus(P2.times(pHcenterSq)).plus(P3.times(pHcenterCd));
//complexCbContribution.show("Complex contribution");
thisPointCb = thisPointCb.plus(complexCbContribution);
}
}
DblCbThatShouldBeZero.setDblAt(thisPointCb,k);
this.normal.setMean(DblMatrix.ZERO);
this.normal.setStd(sigmaCb);
DblMatrix cblike = this.normal.pdf(thisPointCb);
if (DblMatrix.test(DblMatrix.isNaN(cblike)))
{
throw new RuntimeException("NaN for Cb likelihood.");
}
LIKE.setDblAt(cblike,like_index);
like_index++;
if (this.adjustForSmoothingComponent())
{
smoothingDataIndex++;
}
// //If we are being asked to produce a random draw we add it as a
// //report event.
// if (this.isSamplingObject("RandomBufferSolution"))
// {
// ObjectReporter rep = this.objectReporter;
// ObjectReportInstance ev = new ObjectReportInstance();
// ev.setObject(RhsSol);
// rep.addToReport(ev);
// }
//stop1 = System.currentTimeMillis();
//System.out.println("Time it took6a="+(stop1-start1));
//start1 = stop1;
}
//stop = System.currentTimeMillis();
//System.out.println("Time it took4="+(stop-start));
//start = stop;
CbThatShouldBeZero.add(DblCbThatShouldBeZero);
//If we are adjusting for calibration then the pH set here will be
//the latent pH otherwise they are the observed pH.
pH4CbThatShouldBeZero.add(DblpH4CbThatShouldBeZero);
}
}
//DblMatrix.std(pHresids).show("Here is the empirical STD of pH residuals","0.000E0");
//sigmaCalib.show("The current value of sigmaCalib:","0.000E0");
//DblMatrix.Mean(pHresids).show("Here is the empirical MEAN of pH residuals");
//DblMatrix EmpCbStd = DblMatrix.std(CbThatShouldBeZero);
//EmpCbStd.show("Here is the empirical STD of Cb","0.000E0");
//sigmaCb.show("The current value of sigmaCb:","0.000E0");
//
// DblMatrix EmpCbStd = DblMatrix.std(CbThatShouldBeZero);
// if (DblMatrix.test(EmpCbStd.gt(sigmaCb.times(10.0))))
// {
// System.out.println("sigmaCb is probably too small.");
// System.out.println("The empirical std. of Cb was found to be:"+EmpCbStd.getDoubleAt(0));
// System.out.println("sigmaCb is currently:"+sigmaCb.getDoubleAt(0));
// System.out.println("##");
//
// }
//
// if (DblMatrix.test(sigmaCb.gt(EmpCbStd.times(10.0))))
// {
// System.out.println("sigmaCb is probably too large.");
// System.out.println("The empirical std. of Cb was found to be:"+EmpCbStd.getDoubleAt(0));
// System.out.println("sigmaCb is currently:"+sigmaCb.getDoubleAt(0));
// System.out.println("##");
//
// }
//
// DblMatrix someCb = new DblMatrix(10);
// for (int i=0;i<someCb.getN();i++)
// {
// someCb.setDoubleAt(CbThatShouldBeZero.getDoubleAt(i),i);
// }
//someCb.show("CbThatShouldBeZero","0.00E00");
///////////////////////////////////////////////////////////////////////////////
//Report on the residuals!!
if ((this.likelihoodCalls % 1) == 0)
{
FigureReporter rep = this.getFigureReporter("Residuals");
SLFigure fig = rep.getNewFigure();
SLAxes ax = new SLAxes(fig);
ax.setXLabel("pH");
ax.setYLabel("Cb Residual");
SLLegend legend = new SLLegend(ax,ax,SLLegend.ALPHA_ORDER);
legend.setLegendLocale(SLLegend.LEGEND_SOUTH);
pHtoolsPrefs prefs = new pHtoolsPrefs();
prefs.initialize();
DblMatrix Xref = new DblMatrix(2);
//DblMatrix Yref;
for (int r=0;r<pH4CbThatShouldBeZero.size();r++)
{
DblMatrix pHdata = (DblMatrix)pH4CbThatShouldBeZero.get(r);
DblMatrix Cbdata = (DblMatrix)CbThatShouldBeZero.get(r);
Xref.setDblAt(pHdata.min(),0);
Xref.setDblAt(pHdata.max(),1);
//Yref.setDblAt(DblMatrix.min(Cbdata),0);
//Yref.setDblAt(DblMatrix.max(Cbdata),1);
SLLine line = new SLLine(ax);
line.setXData(pHdata);
line.setYData(Cbdata);
line.setMarkerEdgeColor(r);
line.setLineColor(r);
line.setMarkerType(r);
line.updatePreferences("default",prefs);
SLAnnotation annot = (SLAnnotation)legendStrings.get(r);
legend.addLegend(annot.getText(),annot,line);
}
//System.out.println("updatePreferences for legend of Residuals");
legend.updatePreferences("default",prefs);