New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Introduce formula for Karney's direct geodesic method #486

Merged
merged 52 commits into from Jul 25, 2018

Conversation

Projects
None yet
3 participants
@adl1995
Copy link
Contributor

adl1995 commented May 31, 2018

This pull request introduces the direct geodesic method as proposed by Karney (2011). The following are the highlights of the changes introduced:

  • Implementation of the direct geodesic method.
  • Series expansions of integrals given in the paper. All of these were calculated using Maxima with the geod.mac script (associated with GeographicLib), except for co-efficients of C3. I noticed there was a variation in the results produced by GeographicLib and the function generated through Maxima, which caused lon2 to be inaccurate. Replacing the function with Geodesic::C3coeff() from GeographicLib removed this inaccuracy, so I have used that instead.
  • Utility functions in math.hpp, such as normalization, Clenshaw's summation, and Horner's evaluation. Some of these overlap with formulas/area_formulas.hpp. I will merge them in a single file as part of a separate PR.
  • Test cases with existing dataset and antipodal points dataset, collected from GeodTest, associated with GeographicLib. These were converted into C++ array format using this Python script. Geodesic scale (M12) is missing from the GeodTest dataset, so I generated it manually using this C++ script.

However, the tests fail when comparing geodesic scale (M12) with Boost implementation of the direct method. Also, I noticed that using the C++ header GeographicLib/Geodesic.hpp and GeodSolve CLI tool both provided different results. As an example, using these points:

lat1: 31.863784754278001, lon1: 0, azi1: 29.572313410210999, s12: 19993324.8682601

with GeographicLib, GeodSolve CLI, and Boost gave the following results:

GeographicLib version:      -0.99598315084972932
GeodSolve tool:             -0.99588800900114727
Boost version:              -0.99625199787331664

I'm not sure why this happens, because GeodSolve uses the same C++ code at the back-end. The difference might not seem much, but since only a tolerance within 1e-7 is acceptable, these tests fail.

adl1995 added some commits May 12, 2018

[formulas] Add draft of direct geodesic problem from Karney (2011)
The paper can be found at: https://arxiv.org/pdf/1109.4448.pdf
This commit also introduces the evaluate_series_A1 function
for evaluating the series expantion, which was generated
using Maxima: http://maxima.sourceforge.net
[formulas] Add function for evaluating coefficients for C1
- Add SED script for converting x to CT(x)
- Improve code documentation
[util] Evaluate coefficients for C1p using series expansion
- Fix conversion from degree to radian in sin_cos_degrees function
[util] Modify function for evaluting C3x coefficient
- Add separate function for evaluating C3 from C3x coefficient
[test] Add nearly antipodal points dataset for direct geodesic problem
Dataset is collected from:
https://zenodo.org/record/32156

It is then parsed using a Python script.
[test] Add geodesic length to antipodal points dataset
The geodesic length is calculated manually using GeographicLib/Geodesic.hpp
in C++. However, this value differs when calculated using the
CLI tool GeodSolve.

@adl1995 adl1995 referenced this pull request Jun 1, 2018

Open

Progress update #1

*/
template <
typename CT,
size_t SeriesOrder = 8,

This comment has been minimized.

@vissarion

vissarion Jun 4, 2018

Contributor

maybe place it at the end of parameters' list to be more consistent with similar strategies

This comment has been minimized.

@adl1995

adl1995 Jun 4, 2018

Contributor

Thanks! Just updated this.

CT const lat1 = la1;
Azi const azi12 = math::normalize_angle<CT>(azimuth12);

if (math::equals(distance, Dist(0)) || distance < Dist(0))

This comment has been minimized.

@vissarion

vissarion Jun 4, 2018

Contributor

why not creating a const for 0?

This comment has been minimized.

@adl1995

adl1995 Jun 4, 2018

Contributor

Thanks! Just updated this.

math::normalize_angle(lon12));
}

if (BOOST_GEOMETRY_CONDITION(CalcQuantities))

This comment has been minimized.

@vissarion

This comment has been minimized.

@adl1995

adl1995 Jun 4, 2018

Contributor

I actually tried using that code. However, almost all the test cases fail if I do. I find that odd, since it seems to be working well with Vincenty and Thomas. I will try and dig a little deeper into the code and see if I'm missing out on a step.


namespace boost { namespace geometry { namespace series_expansion {

/*

This comment has been minimized.

@vissarion

vissarion Jun 4, 2018

Contributor

there is a lot of duplication in the maxima scripts inside comments, could you try to decrease it?

This comment has been minimized.

@adl1995

adl1995 Jun 5, 2018

Contributor

I just removed the duplicated code in this commit. Please let me know if I should further reduce the comments.

This comment has been minimized.

@vissarion

vissarion Jun 5, 2018

Contributor

It could be more useful to have all maxima script in one place. Because for example if someone want to reproduce series for I2 then the definition of ataylor is missing.

Also the following script (taken from area formula, right?) is repeated 7 times!

To replace each number x by CT(x) the following script can be used: sed -e 's/[0-9]\+/CT(&)/g; s/\[CT/\[/g; s/)\]/\]/g; s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;'

Finally, is the script taken from http://geographiclib.sourceforge.net/html/geod.mac without modifications? If so you can just add the link in the comments.

This comment has been minimized.

@adl1995

adl1995 Jun 5, 2018

Contributor

Also the following script (taken from area formula, right?) is repeated 7 times!
...

Yes, it matches the one from area formula, except the last part s/epsCT(2)/eps2/g;. Before that the script was converting eps2 -> epsCT(2). So, I will just keep a single copy of this script and remove the duplicated ones.

Finally, is the script taken from http://geographiclib.sourceforge.net/html/geod.mac without modifications? If so you can just add the link in the comments.

The scripts are indeed taken from (http://geographiclib.sourceforge.net/html/geod.mac), however, I had to make minor adjustments e.g. change the function header and data types to keep it consistent with Boost Geometry guidelines. Is it okay to discard those changes? If so, I will only provide a link to geod.mac, otherwise, I will add them in a separate file.

This comment has been minimized.

@vissarion

vissarion Jun 5, 2018

Contributor

The scripts are indeed taken from (http://geographiclib.sourceforge.net/html/geod.mac), however, I had to make minor adjustments e.g. change the function header and data types to keep it consistent with Boost Geometry guidelines. Is it okay to discard those changes? If so, I will only provide a link to geod.mac, otherwise, I will add them in a separate file.

Ok, since they contain changes, I think it is useful to keep them in Boost.Geometry.

This comment has been minimized.

@adl1995

adl1995 Jun 5, 2018

Contributor

I have moved the Maxima scripts to: doc/other/maxima/geod.mac.

@vissarion

This comment has been minimized.

Copy link
Contributor

vissarion commented Jun 4, 2018

Thanks for this PR!

Series expansions of integrals given in the paper. All of these were calculated using Maxima with the geod.mac script (associated with GeographicLib), except for co-efficients of C3. I noticed there was a variation in the results produced by GeographicLib and the function generated through Maxima, which caused lon2 to be inaccurate. Replacing the function with Geodesic::C3coeff() from GeographicLib removed this inaccuracy, so I have used that instead.

lon2 is inaccurate with respect to the data set in https://zenodo.org/record/32156 ? How do you know that Geodesic::C3coeff() is more accurate than the function generated from Maxima?

Test cases with existing dataset and antipodal points dataset, collected from GeodTest, associated with GeographicLib. These were converted into C++ array format using this Python script. Geodesic scale (M12) is missing from the GeodTest dataset, so I generated it manually using this C++ script.

Could you please add this information in the code as a comment for future reference?

adl1995 added some commits Jun 4, 2018

adl1995 added some commits Jun 25, 2018

[formulas] Use assignment operator on the same line for consistency
Other changes include the update of series expansion function
calls, as the template arguments are reversed.
[util] Avoid passing array size using std::vector
Previously, the array size was passed in as a
separate parameter.
[util] Add BOOST_GEOMETRY_ASSERT in series expansion and normalizatio…
…n function

Modified functions are:
- evaluate_coeffs_C3x
- normalize_values

math::normalize_values<CT>(sin_beta1, cos_beta1);

cos_beta1 = std::max(math::sqrt(c0), cos_beta1);

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

std::max won't compile with msvc. Take it in brackets. Furthermore sqrt(0) is 0. So:

(std::max)(c0, cos_beta1);

or am I missing something?

CT const lat1 = la1;

Azi azi12 = azimuth12;
math::normalize_angle<degree, Azi>(azi12);

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

So if you called it normalize_azimuth() then it'd be clear that it has to be in range (-180, 180].

// Convert to radians.
remainder *= d2r<T>();

T s = std::sin(remainder), c = std::cos(remainder);

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

Before you used math functions without std::. If you specify std:: then the compiler will not use functions from the namespace where a user-defined numeric type is defined. So please remove std::. If you like you can bring std:: math functions into scope:

using std::floor;
using std::sin;
using std::cos;

but they are not used like that in the whole library so I'd say this is not needed.

// the argument to the range [-45, 45] before converting it to radians.
T remainder; int quotient;

remainder = std::fmod(x, T(360));

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

Use math::mod()

return 0;
}

T y = std::abs(x);

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

math::abs()

template <typename Coeffs, typename CT>
inline void evaluate_coeffs_C3x(Coeffs &c, size_t const& SeriesOrder, CT const& n) {
size_t const coeff_size = Coeffs::static_size;
BOOST_GEOMETRY_ASSERT(coeff_size == (SeriesOrder * (SeriesOrder - 1)) / 2);

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

If you took SeriesOrder as template parameter you could check this at compile time:

static size_t const coeff_size = Coeffs::static_size;
static size_t const expected_size = (SeriesOrder * (SeriesOrder - 1)) / 2;
BOOST_MPL_ASSERT((coeff_size == expected_size));

This comment has been minimized.

@adl1995

adl1995 Jun 27, 2018

Contributor

I'm having some trouble passing SeriesOrder as the template parameter.

I changed the evaluate_coeffs_C3x function header to:

template <typename Coeffs, typename CT, size_t SeriesOrder>
inline void evaluate_coeffs_C3x(Coeffs &c, CT const& n) {

and called it as:

evaluate_coeffs_C3x(*this, n);

but I get this error:

error: ‘evaluate_coeffs_C3x’ was not declared in this scope

I also tried changing the order of parameters as:

template <size_t SeriesOrder, typename CT, typename Coeffs>
inline void evaluate_coeffs_C3x(CT const& n, Coeffs &c) {

but calling the function as:

evaluate_coeffs_C3x(n, *this);

or:

evaluate_coeffs_C3x<SeriesOrder, CT>(n, *this);

doesn't work either.

It produces the error:

error: expected primary-expression before ‘>’ token

and:

couldn't deduce template parameter ‘SeriesOrder’

respectively.

This comment has been minimized.

@awulkiew

awulkiew Jun 27, 2018

Member

If you're calling these functions inside structs' constructors then they should be defined before structs.

Either of these should work:

template <size_t SeriesOrder, typename CT, typename Coeffs>
inline void evaluate_coeffs_C3x(CT const& n, Coeffs &c){}

// this is more strict, it takes only the struct you like
template <size_t SeriesOrder, typename CT>
inline void evaluate_coeffs_C3x(CT const& n, coeffs_C3x<SeriesOrder, CT> &c){}

and either of these should work:

evaluate_coeffs_C3x<SeriesOrder>(n, *this);
evaluate_coeffs_C3x<SeriesOrder, CT>(n, *this);
evaluate_coeffs_C3x<SeriesOrder, CT, coeffs_C3x>(n, *this); // assuming it's called in coeffs_C3x's constructor

This comment has been minimized.

@adl1995

adl1995 Jun 28, 2018

Contributor

Thanks, this worked.

Out of curiosity, why was this working before for the previous function calls? Does this have something to do with passing explicit template parameters? As I only got the error with evaluate_coeffs_C3x<SeriesOrder>.

This comment has been minimized.

@awulkiew

awulkiew Jun 28, 2018

Member

No idea. You'd have to know how exactly type deduction is done in your compiler. Are you using Visual Studio? Its compiler allows things it shouldn't sometimes.

This comment has been minimized.

@adl1995

adl1995 Jun 28, 2018

Contributor

No, I'm using g++ version 7.2.0 on Ubuntu.

int m = Coeffs1::static_size - l - 1;
mult *= eps;

std::vector<CT> coeffs2_slice(coeffs2.begin(), coeffs2.begin() + offset);

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

You don't have to create a container here, one each iteration btw. This will be horribly slow. Pass range into math::polyval().

\brief Evaluate the polynomial.
*/
template<typename CT>
inline CT polyval(std::vector<CT> coeff,

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member
template <typename It>
inline CT polyval(It first, It last, std::iterator_traits<It>::value_type const& eps)

or just

template <typename It, typename T>
inline CT polyval(It first, It last, T const& eps)
inline CT polyval(std::vector<CT> coeff,
const CT eps)
{
int N = boost::size(coeff) - 1;

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

std::distance(first, last) - 1

int N = boost::size(coeff) - 1;
int index = 0;

CT y = N < 0 ? 0 : coeff[index++];

This comment has been minimized.

@awulkiew

awulkiew Jun 26, 2018

Member

*(first + (index++))

adl1995 added some commits Jun 27, 2018

[util][formulas] Rename normalize_angle function to normalize_azimuth
For normalizing longitudes, the normalize_longitude function is
used instead.
[util] Pass SeriesOrder as template parameter in evaluate_coeffs_C3x(…
…) function

The coefficient container structs are moved to the
bottom of the file.
[util] Pass range into math::polyval() instead of std::vector
This is done to avoid creating a separate container in each
iteration.

@adl1995 adl1995 force-pushed the BoostGSoC18:feature/geodesic_direct branch from ab45a02 to 6219503 Jun 28, 2018

Merge branch 'develop' into feature/geodesic_direct
Conflicts:
	include/boost/geometry/util/math.hpp
	test/formulas/direct.cpp

The conflicting files have been updated.
@adl1995

This comment has been minimized.

Copy link
Contributor

adl1995 commented Jun 29, 2018

Thank you for the detailed reviews. I think all the requested changes have been addressed. I have also merged develop into this branch, and resolved all conflicting files.

Please let me know if there are any further changes required.

@awulkiew

This comment has been minimized.

Copy link
Member

awulkiew commented Jul 2, 2018

Thanks! I'm ok with merging. I'd prefer to wait until after the 1.68 release. Or do you need this PR merged in order to work on the next step of the project?

@adl1995

This comment has been minimized.

Copy link
Contributor

adl1995 commented Jul 2, 2018

That's fine. I can work on local branches for the next part i.e. the inverse method.

Btw, when is the next release scheduled for?

@awulkiew

This comment has been minimized.

Copy link
Member

awulkiew commented Jul 3, 2018

Here is the schedule: https://www.boost.org/development/index.html
The branches for final release are closed 1 Aug however after beta is released we only do bugfixing.

I think the whole month is too long to keep you wait so I'm ok with merging it now. If needed I'll create a branch for bugfixing based on previous develop state.

@adl1995

This comment has been minimized.

Copy link
Contributor

adl1995 commented Jul 3, 2018

Okay, as you deem fit. If it is merged before August 14th, I can mention this in my final evaluation report.

@awulkiew awulkiew added this to the 1.69 milestone Jul 25, 2018

@awulkiew awulkiew merged commit 79ef70f into boostorg:develop Jul 25, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment