Skip to content

Commit

Permalink
Issue 108: API extensions to support interpolation
Browse files Browse the repository at this point in the history
  • Loading branch information
gwlucastrig committed Mar 18, 2024
1 parent 32a1a8e commit d7fb48c
Show file tree
Hide file tree
Showing 5 changed files with 449 additions and 46 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
import org.tinfour.common.IIncrementalTin;
import org.tinfour.common.IIncrementalTinNavigator;
import org.tinfour.common.IQuadEdge;
import org.tinfour.common.PolygonConstraint;
import org.tinfour.common.Thresholds;
import org.tinfour.common.Vertex;
import org.tinfour.interpolation.IInterpolatorOverTin;
Expand Down Expand Up @@ -101,7 +102,7 @@ public enum SelectionMethod {
private PrintStream reportSamples;

private final SurfaceModel surfaceModel;
private final List<Vertex>samplesFromRegression = new ArrayList<>();
private final List<Vertex> samplesFromRegression = new ArrayList<>();
private final VertexValuatorDefault defaultValuator = new VertexValuatorDefault();

/**
Expand Down Expand Up @@ -199,11 +200,11 @@ public double interpolate(final double x, final double y, IVertexValuator valuat
return computeRegression(x, y, valuator, edgeQuery, false, computeExtendedStatistics);
}



private double computeRegression(final double x, final double y,
IVertexValuator valuator, IQuadEdge edgeQuery,
boolean crossValidate, boolean computeExtendedStatistics) {
final IVertexValuator valuator,
final IQuadEdge edgeQuery,
final boolean crossValidate,
final boolean computeExtendedStatistics) {
// in the logic below, we access the Vertex x and y coordinates directly
// but we use the getZ() method to get the z value. Some vertices
// may actually be VertexMergerGroup instances and so the Z value must
Expand Down Expand Up @@ -267,14 +268,14 @@ private double computeRegression(final double x, final double y,
innerList.add(e.getB());
}

workTin.clear();
workTin.add(innerList, null);
if (!workTin.isBootstrapped()) {
workTin.clear();
return Double.NaN;
}
workNni.resetForChangeToTin();
NaturalNeighborElements innerElements = workNni.getNaturalNeighborElements(x, y);
workTin.clear();
if (innerElements == null) {
return Double.NaN;
}
Expand Down Expand Up @@ -302,15 +303,15 @@ private double computeRegression(final double x, final double y,
samplesFromRegression.add(C);
}
for (int i = 0; i < vArray.length; i++) {
samplesFromRegression.add(vArray[i]);
samplesFromRegression.add(vArray[i]);
samples[k][0] = vArray[i].getX() - x;
samples[k][1] = vArray[i].getY() - y;
samples[k][2] = vq.value(vArray[i]);
k++;
}
vArray = outerElements.getNaturalNeighbors();
for (int i = 0; i < vArray.length; i++) {
samplesFromRegression.add(vArray[i]);
samplesFromRegression.add(vArray[i]);
samples[k][0] = vArray[i].getX() - x;
samples[k][1] = vArray[i].getY() - y;
samples[k][2] = vq.value(vArray[i]);
Expand Down Expand Up @@ -370,6 +371,35 @@ public void resetForChangeToTin() {
nni.resetForChangeToTin();
}

/**
* Compute the parameters that would be estimated at the position of the
* vertex if it were not incorporated into the Delaunay triangulation.
*
* @param v a valid reference to a vertex instance currently incorporated
* into the Delaunay triangulation.
* @param computeExtendedStatistics compute the extended statistics such
* as the "hat" matrix and R-Student values.
* @return a valid floating-point value if the cross validation is successful;
* otherwise, NaN.
*/
public double crossValidate(Vertex v, boolean computeExtendedStatistics) {
clearResults();
xQuery = v.getX();
yQuery = v.getY();
double x = v.getX();
double y = v.getY();
IQuadEdge edgeQuery = this.getRelatedEdge(x, y);
if (edgeQuery == null) {
return Double.NaN;
}
Vertex A = edgeQuery.getA();
if (!A.equals(v)) {
// did not match a vertex in the TIN with the input
return Double.NaN;
}
return computeRegression(x, y, null, edgeQuery, true, computeExtendedStatistics);
}

/**
* Compute the parameters that would be estimated at the position of the
* vertex if it were not incorporated into the Delaunay triangulation.
Expand All @@ -394,7 +424,7 @@ public double crossValidate(Vertex v) {
// did not match a vertex in the TIN with the input
return Double.NaN;
}
return this.computeRegression(x, y, null, edgeQuery, true, true);
return computeRegression(x, y, null, edgeQuery, true, false);
}

/**
Expand Down Expand Up @@ -502,22 +532,40 @@ private NaturalNeighborElements getAdjacentElements(List<IQuadEdge> envelope) {
}

NaturalNeighborElements nne = null;
workTin.clear(); // not needed once everything else is correctly implemented
workTin.clear();
workTin.add(list, null);
if (workTin.isBootstrapped()) {
workNni.resetForChangeToTin();
nne = workNni.getNaturalNeighborElements(xQuery, yQuery);
}
workTin.clear();

if (nne == null && list.size() >= 3) {
// The NNI class couldn't find an envelope, so we conclude that
// the query point must be outside the convex hull of the working TIN.
// It was inside the original inner envelope (and main TIN),
// but is not inside the outer envelope that we are trying to build here.
// As a fall-back, we attempt to use the points we do have, assigning
// them equal weights. There might be a better way to do this based
// on distances or some other criteria.
List<IQuadEdge> perimeter = workTin.getPerimeter();
PolygonConstraint pc = new PolygonConstraint();
for (IQuadEdge iEdge : perimeter) {
pc.add(iEdge.getA());
}
pc.complete();

double areaOfConvexHull = pc.getArea();
Vertex[] neighbors = list.toArray(new Vertex[0]);
return new NaturalNeighborElements(xQuery, yQuery, new double[list.size()], neighbors, 1.0);
double[] lambda = new double[neighbors.length];
Arrays.fill(lambda, 1.0 / lambda.length);
return new NaturalNeighborElements(
xQuery, yQuery, lambda, neighbors, areaOfConvexHull);
}

return nne;
}

/**
/**
* Gets a safe copy of the computed polynomial coefficients from the
* regression (the "beta" parameters). These coefficients can be used
* for interpolation or surface modeling purposes. Developers
Expand All @@ -531,7 +579,7 @@ private NaturalNeighborElements getAdjacentElements(List<IQuadEdge> envelope) {
* array.
*/
public double[] getBeta() {
if (beta == null) {
if (beta == null) {
return new double[0];
}
return Arrays.copyOf(beta, beta.length);
Expand Down Expand Up @@ -606,7 +654,7 @@ public double getConfidenceIntervalHalfSpan(double populationFraction) {
* @param populationFraction a value in the range 0 &lt; factor &lt; 1.
* @return a positive floating-point value
*/
public double getPredictonIntervalHalfSpan(double populationFraction) {
public double getPredictionIntervalHalfSpan(double populationFraction) {
if (!(0 < populationFraction && populationFraction < 1)) {
throw new IllegalArgumentException("Population fraction is not in the range (0,1): " + populationFraction);
}
Expand Down Expand Up @@ -656,12 +704,13 @@ public RealMatrix getHatMatrix() {

/**
* Get the array of R-Student values corresponding to the input
* sample data. The R-Student values is often a useful tool in detecting
* anomalous points in the input same set.
* sample data. The R-Student values set is often a useful tool in detecting
* anomalous points in the input samples.
* The computation of the R-Student values is optional.
* If the extended-statistics option was not enabled,
* then this method will return a null.
* @return if enabled, a valid array; otherwise, a nukll.
*
* @return if enabled, a valid array; otherwise, a null.
*/
public double[] getRStudentValues() {
return olsSurface.rStudentValues;
Expand All @@ -678,16 +727,56 @@ public int getDegreesOfFreedom() {
return olsSurface.nDegOfFreedom;
}

/**
* Gets a safe copy of the computed standard error for the polynomial
* coefficients from the regression (the error of the "beta" parameters).
* <p>
*
* @return a if the last regression computation was successful, a valid array
* of standard-error values for the parameters corresponding
* to the selected surface model; otherwise, an array of length zero.
*/
public double[] getStandardErrorOfParameters() {
return olsSurface.getStandardErrorOfParameters();
}

/**
* Gets a list containing the vertices used in the most recent regression
* operation. The index of the vertices within the list will correspond
* operation. The index of the vertices within the list will correspond
* to the corresponding positions in the hat matrix and the R-Student Values
* array (if they were computed).
*
* @return a valid, potentially empty list.
*/
public List<Vertex>getVertices(){
ArrayList<Vertex>list = new ArrayList();
public List<Vertex> getVertices() {
ArrayList<Vertex> list = new ArrayList();
list.addAll(samplesFromRegression);
return list;
}

/**
* Get the r-squared value, also known as the coefficient of multiple
* determination.
*
* @return a value in the range 0 to 1.
*/
public double getRSquared() {
return olsSurface.getRSquared();
}

/**
* Get the estimated value at the specified coordinates using the
* parameters from the most recently completed regression operation.
*
* @param xEstimate the x coordinate of interest
* @param yEstimate the y coordinate of interest
* @return if successful, a valid floating-point value; otherwise,
* Double&#46;NaN.
*/
public double computeEstimatedValue(double xEstimate, double yEstimate) {
if (Double.isFinite(this.zResult)) {
return olsSurface.computeEstimatedValue(xEstimate - xQuery, yEstimate - yQuery);
}
return Double.NaN;
}
}
49 changes: 28 additions & 21 deletions analysis/src/main/java/org/tinfour/regression/OlsSurface.java
Original file line number Diff line number Diff line change
Expand Up @@ -762,26 +762,26 @@ public double computeRegression(
// Temporary code to compare results of this class
// (quadratic model only) with the Apache Math implementation
//
//double[] data = new double[nSamples * (5)];
//int k = 0;
//for (int i = 0; i < nSamples; i++) {
// data[k++] = samples[i][2];
// double x = samples[i][0];
// double y = samples[i][1];
// data[k++] = x;
// data[k++] = y;
// data[k++] = x * x;
// data[k++] = y * y;
//}
//OLSMultipleLinearRegression mls = new OLSMultipleLinearRegression();
//mls.newSampleData(data, nSamples, 4);
//double[] q = mls.estimateRegressionParameters();
//double[][] v = mls.estimateRegressionParametersVariance();
//double []erp = mls.estimateRegressionParametersStandardErrors();
//double erv = mls.estimateRegressandVariance();
//double eev = mls.estimateErrorVariance();
//double rss = mls.calculateResidualSumOfSquares();
//double dummy = rss; // to provide a break point for debugging
// int nVar = model.nCoefficients-1;
// double[] data = new double[nSamples * (nVar+1)];
// double[][] qm = computeDesignMatrix(model, nSamples, samples);
// int k = 0;
// for (int i = 0; i < nSamples; i++) {
// data[k++] = samples[i][2];
// for(int j=0; j<nVar; j++){
// data[k++] = qm[i][j+1];
// }
// }
// OLSMultipleLinearRegression mls = new OLSMultipleLinearRegression();
// mls.newSampleData(data, nSamples, nVar);
// double[] q = mls.estimateRegressionParameters();
// double[][] v = mls.estimateRegressionParametersVariance();
// double []erp = mls.estimateRegressionParametersStandardErrors();
// double erv = mls.estimateRegressandVariance();
// double eev = mls.estimateErrorVariance();
// double rss = mls.calculateResidualSumOfSquares();
// double er2 = mls.calculateRSquared();
// double dummy = rss; // to provide a break point for debugging
} catch (SingularMatrixException npex) {
return Double.NaN;
}
Expand Down Expand Up @@ -818,7 +818,7 @@ public double[] getParameters() {
* of standard-error values for the parameters corresponding
* to the selected surface model; otherwise, an array of length zero.
*/
public double[] getStandardErrorOfCoefficients() {
public double[] getStandardErrorOfParameters() {
if (standardErrorOfParameters == null) {
return new double[0];
}
Expand Down Expand Up @@ -1303,4 +1303,11 @@ RealMatrix getHatMatrix() {
return r2Sum;
}

/**
* Get the r-squared value, also known as the coefficient of multiple determination.
* @return a value in the range 0 to 1.
*/
public double getRSquared(){
return r2;
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -262,7 +262,7 @@ public double getAreaOfEnvelope(){
a = b;
}
}
return areaSum/2.0;
return Math.abs(areaSum/2.0);
}

/**
Expand Down
Loading

0 comments on commit d7fb48c

Please sign in to comment.