String exp;
String ssid;
String obsnum;
String param;
String buffername = this.transientStorage.getNickname();
DblMatrix val1,val2;
DblMatrix interceptpH = null;
DblMatrix interceptpHExp = null;
DblMatrix interceptpHExpSigma = null;
DblMatrix interceptpHExpSS = null;
DblMatrix interceptpHExpSSSigma = null;
DblMatrix linearpH = null;
DblMatrix linearpHExp = null;
DblMatrix linearpHExpSigma = null;
DblMatrix linearpHExpSS = null;
DblMatrix linearpHExpSSSigma = null;
DblMatrix quadraticpH = null;
DblMatrix quadraticpHExp = null;
DblMatrix quadraticpHExpSigma = null;
DblMatrix quadraticpHExpSS = null;
DblMatrix quadraticpHExpSSSigma = null;
DblMatrix cubicpH = null;
DblMatrix cubicpHExp = null;
DblMatrix cubicpHExpSigma = null;
DblMatrix cubicpHExpSS = null;
DblMatrix cubicpHExpSSSigma = null;
DblMatrix sigmaCalib = null;
DblMatrix sigmaCb = null;
//DblMatrix initialpH_mean = null;
DblMatrix latentpHi = null;
if (this.numParams()==0)
{
throw new RuntimeException("No parameters have been established.");
}
BufferSolution MarginalSol = null;
param = this.mkParameterName(new String[]{"sigmaCb",buffername});
sigmaCb = this.getParam(param);
sigmaCb.show("sigmaCb","0.00000E0");
//System.out.println("adjustForLatentpHLinkIntercept is "+(((Boolean)this.adjustForMap.get("adjustForLatentpHLinkIntercept")).booleanValue()));
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkIntercept")).booleanValue())
{
param = this.mkParameterName(new String[]{"interceptpH",buffername});
if (this.hasParameter(param))
{
//System.out.println("Has intercept!");
interceptpH = this.getParam(param);
}
else
{
interceptpH = new DblMatrix(0.0);
}
}
else
{
interceptpH = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkLinear")).booleanValue())
{
param = this.mkParameterName(new String[]{"linearpH",buffername});
if (this.hasParameter(param))
{
//System.out.println("Has linear!");
linearpH = this.getParam(param);
}
else
{
linearpH = new DblMatrix(0.0);
}
}
else
{
linearpH = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkQuadratic")).booleanValue())
{
param = this.mkParameterName(new String[]{"quadraticpH",buffername});
if (this.hasParameter(param))
{
//System.out.println("Has quadratic!");
quadraticpH = this.getParam(param);
}
else
{
quadraticpH = new DblMatrix(0.0);
}
}
else
{
quadraticpH = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkCubic")).booleanValue())
{
param = this.mkParameterName(new String[]{"cubicpH",buffername});
if (this.hasParameter(param))
{
//System.out.println("Has cubic!");
cubicpH = this.getParam(param);
}
else
{
cubicpH = new DblMatrix(0.0);
}
}
else
{
cubicpH = new DblMatrix(0.0);
}
this.normal.setMean(new DblMatrix(0.0));
DblMatrix LIKE = DblMatrix.J(2*this.totalNumberOfTitrationPoints);
DblMatrix pHresids = new DblMatrix(2*this.totalNumberOfTitrationPoints);
Vector CbThatShouldBeZero = new Vector();
Vector pH4CbThatShouldBeZero = new Vector();
Vector legendStrings = new Vector();
BufferSolution RhsSol = null;
////
int smoothingDataIndex = 0;
int like_index = 0;
//int cb_index = 0;
if (((Boolean)this.adjustForMap.get("adjustForSmoothing")).booleanValue())
{
smoothingDataIndex = 0;
}
DblMatrix SmoothingpH = null;
DblMatrix SmoothingCbResidual = null;
Lpreg SmoothingSmoother = null;
for (int i=0;i<Exps.size();i++)
{
exp = (String)Exps.get(i);
//RhsSol.add(this.transientStorage,thisConc);
SSIDS = this.transientStorage.getTitrationSubsampleIDs(exp);
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkInterceptExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"interceptpHExp",buffername,exp});
if (this.hasParameter(param))
{
interceptpHExp = this.getParam(param);
}
else
{
interceptpHExp = new DblMatrix(0.0);
}
}
else
{
interceptpHExp = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkLinearExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"linearpHExp",buffername,exp});
if (this.hasParameter(param))
{
linearpHExp = this.getParam(param);
}
else
{
linearpHExp = new DblMatrix(0.0);
}
}
else
{
linearpHExp = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkQuadraticExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"quadraticpHExp",buffername,exp});
if (this.hasParameter(param))
{
quadraticpHExp = this.getParam(param);
}
else
{
quadraticpHExp = new DblMatrix(0.0);
}
}
else
{
quadraticpHExp = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkCubicExp")).booleanValue())
{
param = this.mkParameterName(new String[]{"cubicpHExp",buffername,exp});
if (this.hasParameter(param))
{
cubicpHExp = this.getParam(param);
}
else
{
cubicpHExp = new DblMatrix(0.0);
}
}
else
{
cubicpHExp = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForSmoothingExp")).booleanValue())
{
smoothingDataIndex = 0;
}
for (int j=0;j<SSIDS.size();j++)
{
ssid = (String)SSIDS.get(j);
if (this.adjustForLatentpHLink())//Regardless -- now we know we have at least one titration.
{
param = this.mkParameterName(new String[]{"sigmaCalib",buffername});
sigmaCalib = this.getParam(param);
//sigmaCalib.show(param);
this.normal.setStd(sigmaCalib);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkInterceptExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"interceptpHExpSS",buffername,exp,ssid});
if (this.hasParameter(param))
{
interceptpHExpSS = this.getParam(param);
}
else
{
interceptpHExpSS = new DblMatrix(0.0);
}
}
else
{
interceptpHExpSS = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkLinearExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"linearpHExpSS",buffername,exp,ssid});
if (this.hasParameter(param))
{
linearpHExpSS = this.getParam(param);
}
else
{
linearpHExpSS= new DblMatrix(0.0);
}
}
else
{
linearpHExpSS = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkQuadraticExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"quadraticpH",buffername,exp,ssid});
if (this.hasParameter(param))
{
quadraticpHExpSS = this.getParam(param);
}
else
{
quadraticpHExpSS = new DblMatrix(0.0);
}
}
else
{
quadraticpHExpSS = new DblMatrix(0.0);
}
if (((Boolean)this.adjustForMap.get("adjustForLatentpHLinkCubicExpSS")).booleanValue())
{
param = this.mkParameterName(new String[]{"cubicpH",buffername,exp,ssid});
if (this.hasParameter(param))
{
cubicpHExpSS = this.getParam(param);
}
else
{
cubicpHExpSS = new DblMatrix(0.0);
}
}
else
{
cubicpHExpSS = new DblMatrix(0.0);
}
//This is probably a big time drag...
//stop = System.currentTimeMillis();
//System.out.println("Time it took1="+(stop-start));
//start = stop;
RhsSol = new BufferSolution();
RhsSol.setConnection(this.connection);
RhsSol.setRepository(this.repository);
RhsSol.initialize();
if (this.isSamplingObject("RandomBufferSolution"))
{
MarginalSol = new BufferSolution();
MarginalSol.setConnection(this.connection);
MarginalSol.setRepository(this.repository);
MarginalSol.initialize();
}
//stop = System.currentTimeMillis();
//System.out.println("Time it took2="+(stop-start));
//start = stop;
DblMatrix thisConc = new DblMatrix(1.0);
DblMatrix marginalConc = thisConc.copy();
if (((Boolean)this.adjustForMap.get("adjustForBufferConc")).booleanValue())
{
param = this.mkParameterName(new String[]{"C",buffername});
thisConc = this.getParam(param);
if (this.isSamplingObject("RandomBufferSolution"))
{
if (!((Boolean)this.marginalizeMap.get("adjustForBufferConc")).booleanValue())
{
marginalConc = thisConc.copy();
}
}
}
//thisConc.show("This conc");
//Even though thisConc would not be updated unless the parameters were defined we
//check the boolean flags for self-documentation purposes and clarity.
if (((Boolean)this.adjustForMap.get("adjustForBufferExpConc")).booleanValue())
{
param = this.mkParameterName(new String[]{"CExpdev",buffername,exp});
if (this.hasParameter(param))
{
thisConc = thisConc.plus(this.getParam(param));
if (this.isSamplingObject("RandomBufferSolution"))
{
if (!((Boolean)this.marginalizeMap.get("adjustForBufferExpConc")).booleanValue())
{
marginalConc = thisConc.copy();
}
}
}
}
if (((Boolean)this.adjustForMap.get("adjustForBufferExpSSConc")).booleanValue())
{
param = this.mkParameterName(new String[]{"CExpSSdev",buffername,exp,ssid});
if (this.hasParameter(param))
{
thisConc = thisConc.plus(this.getParam(param));
if (this.isSamplingObject("RandomBufferSolution"))
{
if (!((Boolean)this.marginalizeMap.get("adjustForBufferExpSSConc")).booleanValue())
{
marginalConc = thisConc.copy();
}
}
}
}
//thisConc.show("This concentration");
RhsSol.add(this.transientStorage,thisConc);
if (this.isSamplingObject("RandomBufferSolution"))
{
MarginalSol.add(this.transientStorage,marginalConc);
}
DblMatrix smoothingParameter = null;
if (this.adjustForSmoothingComponent())
{
String mapkey1 = null;
String mapkey2 = null;
String mapkey3 = null;
//These initializations of the maps should only need to be run once under any new parameter configuration.
if (((Boolean)this.adjustForMap.get("adjustForSmoothing")).booleanValue())
{
if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
{
mapkey1 = this.mkParameterName(new String[]{"pH"});
mapkey2 = this.mkParameterName(new String[]{"CbResidual"});
mapkey3 = this.mkParameterName(new String[]{"Lpreg"});
if (!this.smoothingpHMap.containsKey(mapkey1))
{
this.smoothingpHMap.put(mapkey1,DblMatrix.span(0.0,14.0,totalNumberOfTitrationPoints));
this.smoothingCbResidualMap.put(mapkey2,new DblMatrix(totalNumberOfTitrationPoints));
this.smoothingSmootherMap.put(mapkey3,null);
}
param = this.mkParameterName(new String[]{"smoothingParameter",buffername});
smoothingParameter = this.getParam(param);
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.POLYNOMIAL_SMOOTHING)
{
mapkey1 = this.mkParameterName(new String[]{"pH"});
if (!this.smoothingpHMap.containsKey(mapkey1))
{
this.smoothingpHMap.put(mapkey1,DblMatrix.span(0.0,14.0,totalNumberOfTitrationPoints));
}
}
}
else if (((Boolean)this.adjustForMap.get("adjustForSmoothingExp")).booleanValue())
{
if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
{
mapkey1 = this.mkParameterName(new String[]{"pH",exp});
mapkey2 = this.mkParameterName(new String[]{"CbResidual",exp});
mapkey3 = this.mkParameterName(new String[]{"Lpreg",exp});
if (!this.smoothingpHMap.containsKey(mapkey1))
{
int numOfPoints = 0;
for (int ii=0;ii<SSIDS.size();ii++)
{
String ss = (String)SSIDS.get(ii);
OBSNUM = this.transientStorage.getTitrationObservationNumbers(exp,ss);
numOfPoints = numOfPoints + OBSNUM.size();
}
this.smoothingpHMap.put(mapkey1,DblMatrix.span(0.0,14.0,numOfPoints));
this.smoothingCbResidualMap.put(mapkey2,new DblMatrix(numOfPoints));
this.smoothingSmootherMap.put(mapkey3,null);
}
param = this.mkParameterName(new String[]{"smoothingParameterExp",buffername,exp});
smoothingParameter = this.getParam(param);
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.POLYNOMIAL_SMOOTHING)
{
mapkey1 = this.mkParameterName(new String[]{"pH"});
if (!this.smoothingpHMap.containsKey(mapkey1))
{
this.smoothingpHMap.put(mapkey1,DblMatrix.span(0.0,14.0,totalNumberOfTitrationPoints));
}
}
}
else if (((Boolean)this.adjustForMap.get("adjustForSmoothingExpSS")).booleanValue())
{
if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
{
mapkey1 = this.mkParameterName(new String[]{"pH",exp,ssid});
mapkey2 = this.mkParameterName(new String[]{"CbResidual",exp,ssid});
mapkey3 = this.mkParameterName(new String[]{"Lpreg",exp,ssid});
if (!this.smoothingpHMap.containsKey(mapkey1))
{
int numOfPoints = 0;
OBSNUM = this.transientStorage.getTitrationObservationNumbers(exp,ssid);
numOfPoints = numOfPoints + OBSNUM.size();
this.smoothingpHMap.put(mapkey1,DblMatrix.span(0.0,14.0,numOfPoints));
this.smoothingCbResidualMap.put(mapkey2,new DblMatrix(numOfPoints));
this.smoothingSmootherMap.put(mapkey3,null);
}
param = this.mkParameterName(new String[]{"smoothingParameterExpSS",buffername,exp,ssid});
smoothingParameter = this.getParam(param);
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.POLYNOMIAL_SMOOTHING)
{
mapkey1 = this.mkParameterName(new String[]{"pH"});
if (!this.smoothingpHMap.containsKey(mapkey1))
{
this.smoothingpHMap.put(mapkey1,DblMatrix.span(0.0,14.0,totalNumberOfTitrationPoints));
}
}
}
if (smoothingDataIndex == 0) //We are starting anew...
{
if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.LOWESS_SMOOTHING)
{
SmoothingpH = (DblMatrix)this.smoothingpHMap.get(mapkey1);
SmoothingCbResidual = (DblMatrix)this.smoothingCbResidualMap.get(mapkey2);
SmoothingSmoother = (Lpreg)this.smoothingSmootherMap.get(mapkey3);
if (SmoothingSmoother == null)
{
SmoothingSmoother = new Lpreg();
this.smoothingSmootherMap.put(mapkey3,SmoothingSmoother);
}
DblMatrix[] XDATA = new DblMatrix[]{SmoothingpH};
DblMatrix YDATA = SmoothingCbResidual;
SmoothingSmoother.setData(XDATA,YDATA);
NNBandwidth BW = (NNBandwidth)SmoothingSmoother.getBandwidthMethod();
FixedSmoothMethod SM = (FixedSmoothMethod)BW.getSmoothMethod();
DblMatrix SP = smoothingParameter.divideBy(smoothingParameter.plus(1.0));
SM.setSmoothParameter(SP);
BW.setSmoothMethod(SM);
}
else if (((Integer)this.adjustForMap.get("adjustForSmoothingMethod")).intValue() == DefaultComplexBufferPredictorGenerator.POLYNOMIAL_SMOOTHING)
{
SmoothingpH = (DblMatrix)this.smoothingpHMap.get(mapkey1);
}
}
// ComplexBufferPredictor cmplx = RhsSol.getComplexBufferPredictor();
// Smoother smoother = cmplx.getSmoother();
// DblMatrix[] XDATA = new DblMatrix[]{SmoothingpH};
// DblMatrix YDATA = SmoothingCbResidual;
//
// smoother.setData(XDATA,YDATA);
}
if (((Boolean)this.adjustForMap.get("adjustForSmoothingExpSS")).booleanValue())
{
smoothingDataIndex = 0;
}
String[] buffdnames = IconServer.splitDomainName(buffername);
String[] expdnames = IconServer.splitDomainName(exp);
ImageIcon buffDomainIcon = this.iconServer.getIconForDomain(buffdnames[0]);
ImageIcon expDomainIcon = this.iconServer.getIconForDomain(expdnames[0]);
if (buffDomainIcon == null)
{
throw new RuntimeException("The icon for the buffer domain "+buffdnames[0]+" is null.");
}
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);
SLLine refline = new SLLine(ax);
refline.setXData(Xref);
refline.setYData(new DblMatrix(2));
refline.setMarkerType(-1);//No Marker
refline.setLineColor(Color.RED);
FigureReportInstance ev = new FigureReportInstance();
ev.setFigure(fig);
////
ax.updatePreferences("default",prefs);
fig.updatePreferences("default",prefs);
refline.updatePreferences("default",prefs);
////
rep.addToReport(ev);
//////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
//Report on the predicted Cb!!
rep = this.getFigureReporter("Predicted");
if (rep != null)
{
fig = rep.getNewFigure();
ax = new SLAxes(fig);
ax.setXLabel("pH");
ax.setYLabel("Cb");
legend = new SLLegend(ax,ax,SLLegend.ALPHA_ORDER);
legend.setLegendLocale(SLLegend.LEGEND_SOUTH);
prefs = new pHtoolsPrefs();
prefs.initialize();
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);
param = this.mkParameterName(new String[]{"C",buffername});
DblMatrix thisConc = this.getParam(param);
BufferSolution sol = new BufferSolution();
sol.setConnection(this.connection);
sol.setRepository(this.repository);
sol.initialize();
sol.add(this.transientStorage,thisConc);
DblMatrix Cbdata = (DblMatrix)((HashMap)sol.Cb(pHdata)).get("total");//We take a peak at the Cb that would be given by the current parameterization.
Cbdata = DblMatrix.abs(Cbdata);
//We take absolute value here since generally we can't distinguish between +/- acid/base.
Xref.setDblAt(pHdata.min(),0);
Xref.setDblAt(pHdata.max(),1);