This is an automated email from the ASF dual-hosted git repository. erans pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-geometry.git
commit 1c54b6ca49a920435f5aa742db3cc52b1e6f7d2b Author: Matt Juntunen <[email protected]> AuthorDate: Tue Jul 10 22:06:01 2018 -0400 GEOEMTRY-7: removing spherical coordinate code for working with jacobians and hessians since I'm not sure how it would be used; it can be added back later if needed --- .../threed/SphericalDerivativeConverter.java | 241 --------------------- .../threed/SphericalDerivativeConverterTest.java | 187 ---------------- 2 files changed, 428 deletions(-) diff --git a/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/SphericalDerivativeConverter.java b/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/SphericalDerivativeConverter.java deleted file mode 100644 index 45b8794..0000000 --- a/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/SphericalDerivativeConverter.java +++ /dev/null @@ -1,241 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.geometry.euclidean.threed; - -/** Class containing methods for converting gradients and Hessian - * matrices from spherical to Cartesian coordinates. - */ -public class SphericalDerivativeConverter { - - /** Spherical coordinates. */ - private final SphericalCoordinates spherical; - - /** Cartesian vector equivalent to spherical coordinates. */ - private final Vector3D vector; - - /** Jacobian of (r, θ Φ). */ - private double[][] jacobian; - - /** Hessian of radius. */ - private double[][] rHessian; - - /** Hessian of azimuthal angle in the x-y plane θ. */ - private double[][] thetaHessian; - - /** Hessian of polar (co-latitude) angle Φ. */ - private double[][] phiHessian; - - public SphericalDerivativeConverter(SphericalCoordinates spherical) { - this.spherical = spherical; - this.vector = spherical.toVector(); - - computeJacobian(); - } - - /** Return the {@link SphericalCoordinates} for this instance. - * @return spherical coordinates for this instance - */ - public SphericalCoordinates getSpherical() { - return spherical; - } - - /** Return the {@link Vector3D} for this instance. This vector is - * equivalent to the spherical coordinates. - * @return vector for this instance - */ - public Vector3D getVector() { - return vector; - } - - /** Convert a gradient with respect to spherical coordinates into a gradient - * with respect to Cartesian coordinates. - * @param sGradient gradient with respect to spherical coordinates - * {df/dr, df/dθ, df/dΦ} - * @return gradient with respect to Cartesian coordinates - * {df/dx, df/dy, df/dz} - */ - public double[] toCartesianGradient(final double[] sGradient) { - // compose derivatives as gradient^T . J - // the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0 - return new double[] { - sGradient[0] * jacobian[0][0] + sGradient[1] * jacobian[1][0] + sGradient[2] * jacobian[2][0], - sGradient[0] * jacobian[0][1] + sGradient[1] * jacobian[1][1] + sGradient[2] * jacobian[2][1], - sGradient[0] * jacobian[0][2] + sGradient[2] * jacobian[2][2] - }; - } - - /** Convert a Hessian with respect to spherical coordinates into a Hessian - * with respect to Cartesian coordinates. - * <p> - * As Hessian are always symmetric, we use only the lower left part of the provided - * spherical Hessian, so the upper part may not be initialized. However, we still - * do fill up the complete array we create, with guaranteed symmetry. - * </p> - * @param sHessian Hessian with respect to spherical coordinates - * {{d<sup>2</sup>f/dr<sup>2</sup>, d<sup>2</sup>f/drdθ, d<sup>2</sup>f/drdΦ}, - * {d<sup>2</sup>f/drdθ, d<sup>2</sup>f/dθ<sup>2</sup>, d<sup>2</sup>f/dθdΦ}, - * {d<sup>2</sup>f/drdΦ, d<sup>2</sup>f/dθdΦ, d<sup>2</sup>f/dΦ<sup>2</sup>} - * @param sGradient gradient with respect to spherical coordinates - * {df/dr, df/dθ, df/dΦ} - * @return Hessian with respect to Cartesian coordinates - * {{d<sup>2</sup>f/dx<sup>2</sup>, d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dxdz}, - * {d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dy<sup>2</sup>, d<sup>2</sup>f/dydz}, - * {d<sup>2</sup>f/dxdz, d<sup>2</sup>f/dydz, d<sup>2</sup>f/dz<sup>2</sup>}} - */ - public double[][] toCartesianHessian(final double[][] sHessian, final double[] sGradient) { - computeHessians(); - - // compose derivative as J^T . H_f . J + df/dr H_r + df/dtheta H_theta + df/dphi H_phi - // the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0 - // and H_theta is only a 2x2 matrix as it does not depend on z - final double[][] hj = new double[3][3]; - final double[][] cHessian = new double[3][3]; - - // compute H_f . J - // beware we use ONLY the lower-left part of sHessian - hj[0][0] = sHessian[0][0] * jacobian[0][0] + sHessian[1][0] * jacobian[1][0] + sHessian[2][0] * jacobian[2][0]; - hj[0][1] = sHessian[0][0] * jacobian[0][1] + sHessian[1][0] * jacobian[1][1] + sHessian[2][0] * jacobian[2][1]; - hj[0][2] = sHessian[0][0] * jacobian[0][2] + sHessian[2][0] * jacobian[2][2]; - hj[1][0] = sHessian[1][0] * jacobian[0][0] + sHessian[1][1] * jacobian[1][0] + sHessian[2][1] * jacobian[2][0]; - hj[1][1] = sHessian[1][0] * jacobian[0][1] + sHessian[1][1] * jacobian[1][1] + sHessian[2][1] * jacobian[2][1]; - // don't compute hj[1][2] as it is not used below - hj[2][0] = sHessian[2][0] * jacobian[0][0] + sHessian[2][1] * jacobian[1][0] + sHessian[2][2] * jacobian[2][0]; - hj[2][1] = sHessian[2][0] * jacobian[0][1] + sHessian[2][1] * jacobian[1][1] + sHessian[2][2] * jacobian[2][1]; - hj[2][2] = sHessian[2][0] * jacobian[0][2] + sHessian[2][2] * jacobian[2][2]; - - // compute lower-left part of J^T . H_f . J - cHessian[0][0] = jacobian[0][0] * hj[0][0] + jacobian[1][0] * hj[1][0] + jacobian[2][0] * hj[2][0]; - cHessian[1][0] = jacobian[0][1] * hj[0][0] + jacobian[1][1] * hj[1][0] + jacobian[2][1] * hj[2][0]; - cHessian[2][0] = jacobian[0][2] * hj[0][0] + jacobian[2][2] * hj[2][0]; - cHessian[1][1] = jacobian[0][1] * hj[0][1] + jacobian[1][1] * hj[1][1] + jacobian[2][1] * hj[2][1]; - cHessian[2][1] = jacobian[0][2] * hj[0][1] + jacobian[2][2] * hj[2][1]; - cHessian[2][2] = jacobian[0][2] * hj[0][2] + jacobian[2][2] * hj[2][2]; - - // add gradient contribution - cHessian[0][0] += sGradient[0] * rHessian[0][0] + sGradient[1] * thetaHessian[0][0] + sGradient[2] * phiHessian[0][0]; - cHessian[1][0] += sGradient[0] * rHessian[1][0] + sGradient[1] * thetaHessian[1][0] + sGradient[2] * phiHessian[1][0]; - cHessian[2][0] += sGradient[0] * rHessian[2][0] + sGradient[2] * phiHessian[2][0]; - cHessian[1][1] += sGradient[0] * rHessian[1][1] + sGradient[1] * thetaHessian[1][1] + sGradient[2] * phiHessian[1][1]; - cHessian[2][1] += sGradient[0] * rHessian[2][1] + sGradient[2] * phiHessian[2][1]; - cHessian[2][2] += sGradient[0] * rHessian[2][2] + sGradient[2] * phiHessian[2][2]; - - // ensure symmetry - cHessian[0][1] = cHessian[1][0]; - cHessian[0][2] = cHessian[2][0]; - cHessian[1][2] = cHessian[2][1]; - - return cHessian; - } - - /** Evaluates (r, θ, φ) Jacobian. */ - private void computeJacobian() { - - // intermediate variables - final double r = spherical.getRadius(); - final double x = vector.getX(); - final double y = vector.getY(); - final double z = vector.getZ(); - final double rho2 = x * x + y * y; - final double rho = Math.sqrt(rho2); - final double r2 = rho2 + z * z; - - jacobian = new double[3][3]; - - // row representing the gradient of r - jacobian[0][0] = x / r; - jacobian[0][1] = y / r; - jacobian[0][2] = z / r; - - // row representing the gradient of theta - jacobian[1][0] = -y / rho2; - jacobian[1][1] = x / rho2; - // jacobian[1][2] is already set to 0 at allocation time - - // row representing the gradient of phi - jacobian[2][0] = x * z / (rho * r2); - jacobian[2][1] = y * z / (rho * r2); - jacobian[2][2] = -rho / r2; - } - - /** Lazy evaluation of Hessians. */ - private void computeHessians() { - if (rHessian == null) { - - // intermediate variables - final double r = spherical.getRadius(); - final double x = vector.getX(); - final double y = vector.getY(); - final double z = vector.getZ(); - final double x2 = x * x; - final double y2 = y * y; - final double z2 = z * z; - final double rho2 = x2 + y2; - final double rho = Math.sqrt(rho2); - final double r2 = rho2 + z2; - final double xOr = x / r; - final double yOr = y / r; - final double zOr = z / r; - final double xOrho2 = x / rho2; - final double yOrho2 = y / rho2; - final double xOr3 = xOr / r2; - final double yOr3 = yOr / r2; - final double zOr3 = zOr / r2; - - // lower-left part of Hessian of r - rHessian = new double[3][3]; - rHessian[0][0] = y * yOr3 + z * zOr3; - rHessian[1][0] = -x * yOr3; - rHessian[2][0] = -z * xOr3; - rHessian[1][1] = x * xOr3 + z * zOr3; - rHessian[2][1] = -y * zOr3; - rHessian[2][2] = x * xOr3 + y * yOr3; - - // upper-right part is symmetric - rHessian[0][1] = rHessian[1][0]; - rHessian[0][2] = rHessian[2][0]; - rHessian[1][2] = rHessian[2][1]; - - // lower-left part of Hessian of azimuthal angle theta - thetaHessian = new double[2][2]; - thetaHessian[0][0] = 2 * xOrho2 * yOrho2; - thetaHessian[1][0] = yOrho2 * yOrho2 - xOrho2 * xOrho2; - thetaHessian[1][1] = -2 * xOrho2 * yOrho2; - - // upper-right part is symmetric - thetaHessian[0][1] = thetaHessian[1][0]; - - // lower-left part of Hessian of polar (co-latitude) angle phi - final double rhor2 = rho * r2; - final double rho2r2 = rho * rhor2; - final double rhor4 = rhor2 * r2; - final double rho3r4 = rhor4 * rho2; - final double r2P2rho2 = 3 * rho2 + z2; - phiHessian = new double[3][3]; - phiHessian[0][0] = z * (rho2r2 - x2 * r2P2rho2) / rho3r4; - phiHessian[1][0] = -x * y * z * r2P2rho2 / rho3r4; - phiHessian[2][0] = x * (rho2 - z2) / rhor4; - phiHessian[1][1] = z * (rho2r2 - y2 * r2P2rho2) / rho3r4; - phiHessian[2][1] = y * (rho2 - z2) / rhor4; - phiHessian[2][2] = 2 * rho * zOr3 / r; - - // upper-right part is symmetric - phiHessian[0][1] = phiHessian[1][0]; - phiHessian[0][2] = phiHessian[2][0]; - phiHessian[1][2] = phiHessian[2][1]; - } - } -} diff --git a/commons-geometry-euclidean/src/test/java/org/apache/commons/geometry/euclidean/threed/SphericalDerivativeConverterTest.java b/commons-geometry-euclidean/src/test/java/org/apache/commons/geometry/euclidean/threed/SphericalDerivativeConverterTest.java deleted file mode 100644 index 4b9bcbd..0000000 --- a/commons-geometry-euclidean/src/test/java/org/apache/commons/geometry/euclidean/threed/SphericalDerivativeConverterTest.java +++ /dev/null @@ -1,187 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.geometry.euclidean.threed; - -import org.apache.commons.geometry.core.Geometry; -import org.junit.Assert; -import org.junit.Test; - -public class SphericalDerivativeConverterTest { - - private static final double EPS = 1e-10; - - @Test - public void testConstructor() { - // arrange - SphericalCoordinates sc = SphericalCoordinates.of(2, Geometry.PI, Geometry.HALF_PI); - - // act - SphericalDerivativeConverter conv = new SphericalDerivativeConverter(sc); - - // assert - checkSpherical(conv.getSpherical(), 2, Geometry.PI, Geometry.HALF_PI); - checkVector(conv.getVector(), -2, 0, 0); - } - - @Test - public void testToCartesianGradient() { - // NOTE: The following set of test data is taken from the original test for this code in commons-math. - // The test in that project generated and checked the inputs on the fly using the commons-math differentiation - // classes. However, since we don't have the benefit of those here, we're using some selected data points - // from that test. - - // act/assert - checkToCartesianGradient( - SphericalCoordinates.of(0.2, 0.1, 0.1), - new double[] { 3.1274095413292105E-4, -1.724542757978006E-6, 1.5102449769881866E-4 }, - new double[] { 0.0007872851, -0.0000078127, 0.0002357921 }); - - checkToCartesianGradient( - SphericalCoordinates.of(0.2, 0.1, 1.6), - new double[] { -7.825830329191124E-8, 7.798528724837122E-10, -4.027286034178383E-7 }, - new double[] { -0.0000000197, 0.0000000019, 0.0000020151 }); - - checkToCartesianGradient( - SphericalCoordinates.of(0.2, 1.6, 0.1), - new double[] { -9.075903886546823E-6, -1.5573157416535893E-5, -4.352284221940998E-6 }, - new double[] { 0.0007802833, 0.0000002252, -0.000006858 }); - - checkToCartesianGradient( - SphericalCoordinates.of(0.2, 2.4, 2.4), - new double[] { 6.045188551967462E-4, 2.944844493772992E-5, 5.207279563401837E-5 }, - new double[] { -0.0003067696, -0.0000146129, -0.0006216347 }); - - checkToCartesianGradient( - SphericalCoordinates.of(9.2, 5.5, 2.4), - new double[] { 27.09285722408859, 327.829199283976, 422.53939642005736 }, - new double[] { 26.1884919572, 48.3685006936, -51.0009075025 }); - } - - private void checkToCartesianGradient(SphericalCoordinates spherical, double[] sGradient, double[] cGradient) { - SphericalDerivativeConverter conv = new SphericalDerivativeConverter(spherical); - - double[] result = conv.toCartesianGradient(sGradient); - - Assert.assertArrayEquals(cGradient, result, EPS); - } - - @Test - public void testToCartesianHessian() { - // NOTE: The following set of test data is taken from the original test for this code in commons-math. - // The test in that project generated and checked the inputs on the fly using the commons-math differentiation - // classes. However, since we don't have the benefit of those here, we're using some selected data points - // from that test. - // - // The NaN values in the input spherical Hessians are only present to ensure that the upper-right - // part of the matrix is not used in the calculation. - - // act/assert - checkToCartesianHessian( - SphericalCoordinates.of(0.2, 0.0, 0.1), - new double[] { 3.147028015595093E-4, -1.5708927954007288E-7, 1.5209020574753025E-4 }, - new double[][] { - { 0.004720542023392639, Double.NaN, Double.NaN }, - { -3.927231988501822E-6, -1.5732003526076452E-5, Double.NaN }, - { 0.0030418041149506037, -3.0840214797113795E-6, -1.56400962465978E-4 } - }, - new double[][] { - { 0.0, -3.940348984959686E-4, 0.011880399467047453 }, - { -3.940348984959686E-4, 7.867570038987733E-6, -1.1860608699245036E-4 }, - { 0.011880399467047453, -1.1860608699245036E-4, 0.002384031969540735 } - }); - - checkToCartesianHessian( - SphericalCoordinates.of(0.2, 0.2, 1.7), - new double[] { -6.492205616890373E-6, 9.721055406032577E-8, -7.490005649457144E-6 }, - new double[][] { - { -9.660140526063848E-5, Double.NaN, Double.NaN }, - { 2.087263937942704E-6, 3.0135301759512823E-7, Double.NaN }, - { -1.4908056742242714E-4, 2.228225255291761E-6, -1.1271700251178201E-4 } - }, - new double[][] { - { 0.0, 8.228328248729827E-7, 1.9536195257978514E-4 }, - { 8.228328248729827E-7, -1.568516517220037E-7, -1.862033454396115E-5 }, - { 1.9536195257978514E-4, -1.862033454396115E-5, -0.0029473017314775615 } - }); - - checkToCartesianHessian( - SphericalCoordinates.of(0.2, 1.6, 0.1), - new double[] { -9.075903886546686E-6, -1.5573157416535897E-5, -4.352284221940931E-6 }, - new double[][] { - { -1.3557892633841054E-4, Double.NaN, Double.NaN }, - { -3.106944464923055E-4, 4.4143436330613375E-7, Double.NaN }, - { -8.660889278565699E-5, -1.489922640116937E-4, 5.374400993902801E-6 } - }, - new double[][] { - { 0.0, -3.862868527078941E-4, 0.011763015339492582 }, - { -3.862868527078941E-4, -2.229868350965674E-7, 3.395142163599996E-6 }, - { 0.011763015339492582, 3.395142163599996E-6, -6.892478835391066E-5 } - }); - - checkToCartesianHessian( - SphericalCoordinates.of(0.2, 2.4, 2.5), - new double[] { 6.911538590806891E-4, 3.344602742543664E-5, 3.330643810411849E-5 }, - new double[][] { - { 0.010200457858547542, Double.NaN, Double.NaN }, - { 6.695363800209198E-4, -3.070347513695088E-5, Double.NaN }, - { 6.68380906286568E-4, 3.001744637007274E-5, -2.273032055462482E-4 } - }, - new double[][] { - { 0.0, 1.9000713243497378E-4, 0.007402721147059207 }, - { 1.9000713243497378E-4, 1.6118798431431763E-5, 3.139960286869248E-4 }, - { 0.007402721147059207, 3.139960286869248E-4, 0.008155571186075681 } - }); - - checkToCartesianHessian( - SphericalCoordinates.of(9.2, 5.6, 2.5), - new double[] { 41.42645719593436, 859.1407583470807, 939.7112322238082 }, - new double[][] { - { 11.642163255436742, Double.NaN, Double.NaN }, - { 54.8154280776715, 5286.1651942531325, Double.NaN }, - { 60.370567966140726, 4700.570567363823, 4929.996883244262 } - }, - new double[][] { - { 0.0, 36.772022140868714, -22.087375306566134 }, - { 36.772022140868714, 212.8111723550033, -63.91326828897971 }, - { -22.087375306566134, -63.91326828897971, 25.593304575600133 } - }); - } - - private void checkToCartesianHessian(SphericalCoordinates spherical, double[] sGradient, - double[][] sHessian, double[][] cHessian) { - SphericalDerivativeConverter conv = new SphericalDerivativeConverter(spherical); - - double[][] result = conv.toCartesianHessian(sHessian, sGradient); - - Assert.assertEquals(cHessian.length, result.length); - for (int i=0; i<cHessian.length; ++i) { - Assert.assertArrayEquals("Hessians differ at row " + i, cHessian[i], result[i], EPS); - } - } - - private void checkSpherical(SphericalCoordinates c, double radius, double azimuth, double polar) { - Assert.assertEquals(radius, c.getRadius(), EPS); - Assert.assertEquals(azimuth, c.getAzimuth(), EPS); - Assert.assertEquals(polar, c.getPolar(), EPS); - } - - private void checkVector(Vector3D v, double x, double y, double z) { - Assert.assertEquals(x, v.getX(), EPS); - Assert.assertEquals(y, v.getY(), EPS); - Assert.assertEquals(z, v.getZ(), EPS); - } -}
