// The Lowess fitter requires a strictly increasing sequence
// on the domain (i.e. JD values), i.e. no duplicates.
Map<Double, Double> jdToMagMap = new TreeMap<Double, Double>();
for (int i = 0; i < obs.size(); i++) {
ValidObservation ob = obs.get(i);
// This means that the last magnitude for a JD wins!
jdToMagMap.put(ob.getJD(), ob.getMag());
}
double[] xvals = new double[jdToMagMap.size()];
double[] yvals = new double[jdToMagMap.size()];
int index = 0;
for (Double jd : jdToMagMap.keySet()) {
xvals[index] = jd;
yvals[index++] = jdToMagMap.get(jd);
}
try {
final LoessInterpolator interpolator = new LoessInterpolator();
function = interpolator.interpolate(xvals, yvals);
fit = new ArrayList<ValidObservation>();
residuals = new ArrayList<ValidObservation>();
double sumSqResiduals = 0;
String comment = "From Lowess fit";
// Create fit and residual observations and
// compute the sum of squares of residuals for
// Akaike and Bayesean Information Criteria.
for (int i = 0; i < xvals.length && !interrupted; i++) {
double jd = xvals[i];
double mag = yvals[i];
double y = function.value(jd);
ValidObservation fitOb = new ValidObservation();
fitOb.setDateInfo(new DateInfo(jd));
fitOb.setMagnitude(new Magnitude(y, 0));
fitOb.setBand(SeriesType.Model);
fitOb.setComments(comment);
fit.add(fitOb);
ValidObservation resOb = new ValidObservation();
resOb.setDateInfo(new DateInfo(jd));
double residual = mag - y;
resOb.setMagnitude(new Magnitude(residual, 0));
resOb.setBand(SeriesType.Residuals);
resOb.setComments(comment);
residuals.add(resOb);
sumSqResiduals += (residual * residual);
}