Progress update

[adl1995]

@vissarion, @awulkiew - I will use this issue for providing updates on my progress. I believe this will make it convenient for you to provide feedback / review code.

During the community bonding period I

• Made a blog post describing my GSoC project and demonstrated the inaccuracy in current implementation of Boost Geometry for nearly antipodal points;
• Replicated the GenInverse method from GeographicLib to compute distance between two points. The code is present in adl1995/boost-geometry-proposal repository. This allowed me to get a bit familiar with the steps involved in the inverse problem;
• Started the implementation of direct geodesic problem (geometry/util/karney_direct.hpp) from Karney (2011) under the feature/geodesic_direct branch in this repository. I initially started implementing the inverse problem, but found that it was more complex and required additional steps, whereas, the direct problem is simpler and non-iterative. Since the inverse problem builds upon the direct problem, this implementation will allow me to better understand the mathematical details involved (such as series expansion). This is also how I proposed to do so in my project timeline;
• Got familiar with Maxima software and generated the series expansion for an integral (evaluate_series_A1) by making minor modifications in this script.

As I mentioned in #421, there is some preprocessing involved in the approach proposed by Karney. For now, I'm including these functions in geometry/util/math.hpp, as these can be reused in the inverse problem.

Currently, the code is available in karney_direct.hpp file, once this is complete, I will try and bundle this with Vincenty.

[vissarion]

@adl1995 Thanks for the update!

Note that there are formulas about series expansions / maxima scripts / Horner and Clenshaw summation in other places of the library, most notably here: area_formulas.hpp That could be useful to your project. Also feel free to propose a place for all of these utilities to avoid duplicates, geometry/util/math.hpp could be such a place or better geometry/util/series_expansions.hpp or similar.

[adl1995]

Update 19/05/2018:

• [utils] Add function for normalizing and rounding off an angle
• [formulas] Compute sin cos series using Clenshaw summation
• [util] Move series expansion functions to util/series_expansion.hpp

@vissarion Thanks for the suggestion. I looked through area_formulas.hpp, however, the series expansions I've currently used are different from the one present there. But, I will reuse them if I find a duplication.

Also, I have moved the series expansions from formulas/karney_direct.hpp to util/series_expansion.hpp.

[vissarion]

Indeed series for area are different but it is good to keep the same or similar format and eventually move everything related to series in the same folder. I just remembered that I have a prototype implementation for direct formula here. Feel free to use/modify/test it if you like.

[adl1995]

Update 25/05/2018:

Commits:

• [util] Evaluate coefficients for C1p using series expansion
• [formulas] Compute the latitude for second point following Karney's method
• [util] Add functions to evaluate coefficients for A3 and C3 using series expansion
• [formulas] Compute the longitude for second point using the longitudinal difference
• [util] Evaluate series and coefficients for A2 and C2 using series expansion
• [formulas] Compute the reduced length and geodesic scale using Karney's direct method
• [formulas] Make variable declarations constant

The implementation for the direct problem is complete. I'm now working on adding the test cases. The dataset I'm using is GeodTest associated with GeographicLib. It contains 500,000 entries for nearly antipodal points. However, I'm not sure where the current dataset in direct_cases.hpp is collected from. It contains three separate results for Thomas, Vincenty, and Karney's method. In this case, Thomas and Vincenty's method won't provide correct results, so should I only test this against Karney's method? @vissarion, @awulkiew - can you please comment on this?

Also, does Boost Geometry follow a convention for passing / returning data in degrees or in radians? I noticed in the test cases that returned data is converted into degrees before the actual comparison is performed. I'm not sure if user has the ability to specify this as part of the function call.

[vissarion]

The implementation for the direct problem is complete. I'm now working on adding the test cases. The dataset I'm using is GeodTest associated with GeographicLib. It contains 500,000 entries for nearly antipodal points. However, I'm not sure where the current dataset in direct_cases.hpp is collected from. It contains three separate results for Thomas, Vincenty, and Karney's method. In this case, Thomas and Vincenty's method won't provide correct results, so should I only test this against Karney's method? @vissarion, @awulkiew - can you please comment on this?

direct_cases.hpp contains data created by Adam. For the cases of nearly antipodal points there is no need to test against Thomas and Vincenty.

Also, does Boost Geometry follow a convention for passing / returning data in degrees or in radians? I noticed in the test cases that returned data is converted into degrees before the actual comparison is performed. I'm not sure if user has the ability to specify this as part of the function call.

As far as I understand, there is no issue here since formulas are internal i.e. users have access to algorithms. For the inverse case there is distance algorithm that takes a strategy as input and the strategy calls either Andoyer or Thomas or Vincenty formula. For the direct case a similar algorithm is needed.

[adl1995]

Update 01/06/2018:

Commits:

• [formulas][util] Add missing import and function return type
• [test] Add direct Karney's method to direct.cpp test cases
• [util] Modify function for evaluting C3x coefficient
• [util] Add functions for normalizing and evaluating polynomial
• [formulas] Fix direct geodesic method by performing normalization
• [test] Add nearly antipodal points dataset for direct geodesic problem
• [test] Add geodesic length to antipodal points dataset
• [test] Test Karney's method on antipodal points dataset
• [util] Update series expansion for C3x

Other work:

• Made a blog post on using variadic templates with lambda expressions.
• Opened up a pull request introducing Karney's direct method.

As far as I understand, there is no issue here since formulas are internal i.e. users have access to algorithms. For the inverse case there is distance algorithm that takes a strategy as input and the strategy calls either Andoyer or Thomas or Vincenty formula. For the direct case a similar algorithm is needed.

@vissarion Thanks for pointing that out! So, for each direct case i.e. Thomas, Vincenty, and Karney, I will introduce a strategy and add a coordinates algorithm(?) which will call the respective strategy. Which strategy should I select as default?

Also, is there a way to run CircleCI builds for pull requests in the main Boost Geometry repository? I wasn't able to resolve the issue related with running builds on this fork.

[vissarion]

@vissarion Thanks for pointing that out! So, for each direct case i.e. Thomas, Vincenty, and Karney, I will introduce a strategy and add a coordinates algorithm(?) which will call the respective strategy. Which strategy should I select as default?

To be consistent with the distance case (i.e. inverse formula) the default should be Andoyer. However, there is no Andoyer equivalent formula for the direct case. To closest choice could be Thomas with 1st order series (instead of 2nd order) see https://github.com/boostorg/geometry/pull/431

For now I would suggest to focus on inverse problem for nearly antipodal points and when is done return back on this. I will try to remove some formulas that appear to be redundant in https://github.com/boostorg/geometry/pull/431 and submit a new PL.

Also, is there a way to run CircleCI builds for pull requests in the main Boost Geometry repository? I wasn't able to resolve the issue related with running builds on this fork.

I think you cannot run on main BG repo. One workaround (a bit inconvenient) is to fork on your account and run CircleCI tests on that fork.

@awulkiew what do you think?

[awulkiew]

Are you talking about your PR: https://github.com/boostorg/geometry/pull/486? It seems that the test was performed even though you used feature/... branch name (tests should be disabled for branches other than develop and master, see below). I don't know why it is the case, it should simply not exist. But anyway, it seems that the test is performed but the result is failure due to the lack of g++-4.9 package. I don't know why is that since in the main repository everything is fine. If you want to play with it to figure out why it is the case you can add some lines to circle.yml script printing some diagnostics like system version or some packages info etc.

Side note. Running tests on branches different than develop and master is disabled (see: https://github.com/boostorg/geometry/blob/develop/circle.yml#L10). It's because running all of the Boost.Geometry tests takes long time and we share CircleCI containers with all other Boost libraries so we could interfere with others. But also because the sript is sending info about code coverage manually to a single Coveralls account using $COVERALLS_REPO_TOKEN defined in CircleCI by me. This is why the tests for PRs are not performed. This could be changed, i.e. we should execute a lighter test for other branches and not send info about coverage.

[awulkiew]

I see on the PR page that you pushed more commits yet the tests were not performed. Why the test were performed before I don't know. Maybe it was some CircleCI anomaly.

[adl1995]

@awulkiew I actually ran that build manually from CircleCI dashboard. I've tried running tests on BoostGSoC18/geometry and adl1995/geometry before, but they both fail due to the lack of g++-4.9 package.

I will try tweaking the circle.yml script and see if I can find out what's causing the issue.

[adl1995]

@vissarion For the inverse problem, should I handle the case for meridional and equatorial points? Or, is this functionality already implemented in formulas/elliptic_arc_length.hpp?

[vissarion]

It would be better to implement a version of meridional and equatorial that is consistent with the new general method, if possible. Then we can test against the current implementation.

[adl1995]

Update 09/06/2018:

Commits (direct problem):

• [test] Update geodesic scale (M12) in GeographicLib dataset
• [test] Use series order 2 for Karney's direct method
• [util] Use Maxima generated function for computing C3x coefficients
• [util] Move sin_cos_series function to series_expansion.hpp
• [doc] Move Maxima scripts for geodesics to doc/other/maxima/geod.mac
• [test] Add comment providing dataset source and how it is parsed
• [util] Remove duplicated Maxima code from series_expansion.hpp
• [formulas] Use constant type variables for comparison
• [formulas] Move SeriesOrder to the end of template parameter list

Commits (inverse problem):

• [formulas] Add function for computing the length at the meridians
• [formulas] Arrange points in canonical form for inverse geodesic problem
• [util] Add function for returning NaN (not a number)
• [util] Add functions to normalize / sum two given values (angles)

[vissarion]

@adl1995 I have created the PR mentioned above with some direct formulas picked from an old PR. You can have a look here: https://github.com/boostorg/geometry/pull/489

[adl1995]

Update 15/06/2018:

Commits (direct problem):

• [doc][util][formulas][test] Add copyright information in updated files
• [formulas][util] Improve code formatting to conform with guidelines
• [formulas] Use updated functions for normalization
• [util] Use existing normalize_spheroidal_coordinates class for normalizing an angle
• [formulas] Use template argument CT instead of double
• [formulas] Fix calculation of t for finding the geodesic scale (M12)
• [test] Calculate geodesic scale (M12) using high precision arithmetic

Commits (inverse problem):

• [formulas] Add function lambda12 to regulate bracketing range in Karney inverse
• [formulas] Perform normalization on starting guess if it passes the sanity check
• [formulas] Solve the astroid equation for inverse problem
• [formulas] Flip sign of cos_lam12
• [formulas] Find starting point for inverse problem (short lines)
• [formulas] Handle case for equatorial points in inverse problem

[adl1995]

Update 22/06/2018:

Commits (direct problem):

• [formula][util] Pass boost::array to series expansion functions

Commits (inverse problem):

• [formulas] Ensure reduced length and geodesic scale are computed in canonical form
• [formulas] Use midpoint of bracket when value lies outside of range
• [formulas] Update bracketing values in Newton's method
• [util] Move difference_angle function to normalize_spheroidal_coordinates.hpp
• [formulas][util] Improve code formatting to conform with guidelines
• [formulas] Use updated functions for normalization
• [util] Use existing normalize_spheroidal_coordinates class for normalizing an angle

[adl1995]

Update 29/06/2018:

Commits (direct problem):

• [util] Pass range into math::polyval() instead of std::vector
• [util] Pass SeriesOrder as template parameter in evaluate_coeffs_C3x() function
• [util] Rename math::normalize_values to math::normalize_unit_vector
• [util] Use functions from math namespace instead of std
• [util][formulas] Rename normalize_angle function to normalize_azimuth
• [util] Change static inline to inline in series_expansion.hpp file
• [util] Add BOOST_GEOMETRY_ASSERT in series expansion and normalization function
• [util] Avoid passing array size using std::vector
• [formulas] Use assignment operator on the same line for consistency
• [util] Reverse template argument order for series expansion functions
• [formulas] Use namespace alias se for series_expansion
• [formaulas][util] Define coefficient containers for computing series expansions

[adl1995]

Update 06/07/2018:

Commits (inverse problem):

• [formulas] Store values from Karney's inverse method in result_inverse structure
• [formulas][test] Add Karney's inverse method in inverse test cases
• [test] Add nearly antipodal points dataset for inverse geodesic problem
• [formulas] Resolve inaccuracy in starting point for Newton's method

The implementation of the inverse method is nearly complete. I'm currently trying to resolve inaccuracies in the calculation, which leads some test cases to failure.

[adl1995]

Update 13/07/2018:

Commits (inverse problem):

• [formulas][util] Reformat code in karney_inverse to use coefficient containers
• [formulas] Use namespace alias se for series_expansion in karney_inverse
• [formulas] Fix incorrect argument to meridian_length() function
• [formulas] Remove unused variable a12 in karney_inverse

This completes the implementation of Karney's inverse method. I will now add a strategy for this method, similar to currently available methods in Boost Geometry i.e. Vincenty, Thomas, and Andoyer.

[adl1995]

Update 20/07/2018:

Commits (inverse problem):

• [formulas][test] Update copyright information
• [strategies] Add distance strategy for Karney's inverse formula

[adl1995]

@vissarion, @awulkiew - Since the Karney's direct and inverse method are now implemented (Cf. #486, #500), I'm trying to find other useful features to add. Listed below are the features I think might benefit the Boost Geometry library:

• Benchmark the performance of Karney's method (direct, inverse) against other available formulas
• Add examples for distance strategies (this could be a general example)
• Add strategies for direct formulas (currently only inverse formulas have distance strategies)
• Identify if the points passed to the distance function are antipodal, if so, call the appropriate formula i.e. Karney's method (similar to what is currently done for meridional points here)

Do you think either of these would be a useful addition to the library? If you have other feature requests that might be useful, please let me know.

[vissarion]

What you propose make sense as potential useful features. Regarding strategies for direct formulas I am not sure if it is needed since strategies are used by algorithms, so there should be an algorithm that need such a strategy.

Maybe the most useful feature related to your work is to make andoyer and/or thomas formulas work with antipodal points. I do not know whether this is a hard task but I think it worth a try.

[adl1995]

Maybe the most useful feature related to your work is to make andoyer and/or thomas formulas work with antipodal points. I do not know whether this is a hard task but I think it worth a try.

Since Karney's method is tightly coupled with rest of the code, figuring out which portion to replace in Andoyer, Thomas, or Vincenty's method might be challenging, but I will check the paper and see if it provides a solution to modify existing formulas.

Instead, would it be beneficial if the general distance function (for inverse methods) checks if the points passed to it are antipodal (we can find a general threshold through trial-and-error)? And if so, it can call Karney's method to produce the correct result.

[vissarion]

Instead, would it be beneficial if the general distance function (for inverse methods) checks if the points passed to it are antipodal (we can find a general threshold through trial-and-error)? And if so, it can call Karney's method to produce the correct result.

That would cause inconsistencies. Assume that you use distance with andoyer formula and you use the method you propose for antipodal. Then the distances for two points nearly antipodal and two points close to the threshold of being nearly antipodal but out of the "nearly antipodal region" will have a large difference in accuracy.

[adl1995]

Update 27/07/2018:

Commits:

• [example] Add example on distance formula

@vissarion, @awulkiew - Regarding the issue for updating Andoyer, Thomas, and Vincenty's method to work well with nearly antipodal points. From what I understand, the implementation of Karney's method starts off differently from other methods i.e. it involves computing series expansion of integrals and relies on modified equations from numerous geodesic algorithms. Given this, I don't think there is a plug-in solution for existing methods in Boost Geometry to produce accurate results for nearly antipodal points.

I think the documentation in the library could be improved by mentioning the pros and cons for using either of the geodesic algorithms, and notifying the users to only use Karney's method for nearly antipodal points. Would this be useful? Please let me know if this is already mentioned somewhere.

As only two weeks are remaining for this GSoC, I want to utilize them by doing something productive. During the next week, I will do benchmarking between the available geodesic algorithms and report their efficiency and accuracy. If you have any other suggestions, please let me know.

[vissarion]

Sure, you can work on documentation (see also the relevant discussion here: https://github.com/boostorg/geometry/pull/490). Benchmarks sounds very interesting to me. Looking forward to see the results. You can find some old benchmarks here: https://github.com/boostorg/geometry/pull/431

[adl1995]

@vissarion - Thanks for the links.

I looked through boostorg#431. I am first benchmarking the inverse methods against Karney. Following your benchmarks, for the error you report in the table, is the difference between the expected distance and the resulting distance using the algorithm (maximized over all runs)?

Is there any other useful performance metric that I could report except for accuracy and time?

Also, do you think (https://github.com/awulkiew/benchmark-geometry) would fit this use-case? I could export the results to JSON format and create graphs similar to this: https://awulkiew.github.io/benchmark-geometry.

[adl1995]

Update 03/08/2018:

This week I performed performance benchmarks of geodesic methods (Cf. #3). For now, I have mainly compared the execution times of direct and inverse methods and showed them in tabular format.

During the next week, I will try to visualize the benchmarks using graphs, and also prepare my GSoC final evaluation report. I plan on making a blog post, similar to what I did for the previous GSoC.