Performance and stability improvements to wcs.all_world2pix - 2nd iterat...

This PR is a second iteration of PR #2373 with the main difference that now the changes are on a branch (instead of master), hopefully does not have (yet) merge conflicts, and addresses the issue of failing doctests from PR #2373 by using data file name relative to the package as done in PR #2807 and with updated output of the examples.

This (as did #2373) provides significant improvements to the performance (typically 1000-2000 times faster) and stability of `wcs.all_world2pix` function. The used algorithm is a completely vectorized implementation of the method of consecutive approximations (fixed point iterations). While this method is a low-order method, it is more than adequate for geometric distortions that are typical in astrophysical images (smooth, slowly varying polynomials).

Compared to the original implementation of `all_world2pix`, this implementation does not require users to have SciPy installed in order to run.

While there is a check to see if WCS contains higher-order distortions (and return the exact solution if it does not), the method itself does not *require* this check and will converge to the exact solution after one iteration if WCS does not contain higher order distortions. This check is mostly for performance purposes to avoid unnecessary iterative solution for images that do not contain distortions (and pix->sky requires only WCS-based, i.e., CD-only transformations).

The proposed replacement for the current implementation of `all_world2pix` avoids unnecessary use of `wcs_world2pix` and `all_pix2world` such as `wcs_world2pix(all_pix2world(pix))` by replacing this with `pix2foc`. This avoids unnecessary CD-matrix based transformation CD(CD^{-1})=I which in itself results in about 50% performance improvement over `wcs_world2pix(all_pix2world(pix))` and likely 25% over `all_pix2world` alone.

The function has additional code to check for *divergence* condition (as opposite to *slow convergence*) when consecutive approximations *diverge*. This can happen when the pixel coordinates are *far* away from the reference pixel (distortions are "severely" extrapolated), or for poorly-defined geometric distortions (at least for this numerical method).

A new and improved `NoConvergence` exception has attributes that allows users to retrieve best solution reached by the function, and indices of *slowly converging* solutions (points) that *did not converge* in the specified maximum number of iterations as well as indices of *diverging* solutions (likely an indication of poorly described distortions or very large distortions).

Finally, the new test function `test_all_world2pix_direct` (proposed as a replacement to the old test) provides a much more improved (versatile) and *thorough* test of the `all_world2pix` function.

These improvements address the issues raised in #1977, https://github.com/scipy/scipy/issues/3227 that plagued `all_world2pix`, and #2294 (in addition to performance improvements).

To keep things together, in the next comment I will outline the proposed method and code. For a mode detailed discussion and comments - see #2373.
>  Hi Leo,
>
>  Could you, please, take a look at doctests for
> https://github.com/astropy/astropy/pull/2816 - they are failing due to EOF
> after line 1978. I have already fixed a previous error but now this...
> 
>  Thank you!
>
>  Mihai

In doctests, I think that compound statements are continued to the next line with an ellipsis. Try changing this:

    >>> try:
    >>>   xy = w.all_world2pix(divradec, 1, maxiter=20,
    ...                        tolerance=1.0e-4, adaptive=False,
    ...                        detect_divergence=True,
    ...                        quiet=False)
    >>> except wcs.wcs.NoConvergence as e:
    >>>   print("Indices of diverging points: {{}}"
    ...         .format(e.divergent))
    >>>   print("Indices of poorly converging points: {{}}"
    ...         .format(e.slow_conv))
    >>>   print("Best solution:\\n{{}}".format(e.best_solution))
    >>>   print("Achieved accuracy:\\n{{}}".format(e.accuracy))
    >>>   raise e
    Indices of diverging points: [1]

to this:

    >>> try:
    ...   xy = w.all_world2pix(divradec, 1, maxiter=20,
    ...                        tolerance=1.0e-4, adaptive=False,
    ...                        detect_divergence=True,
    ...                        quiet=False)
    ... except wcs.wcs.NoConvergence as e:
    ...   print("Indices of diverging points: {{}}"
    ...         .format(e.divergent))
    ...   print("Indices of poorly converging points: {{}}"
    ...         .format(e.slow_conv))
    ...   print("Best solution:\\n{{}}".format(e.best_solution))
    ...   print("Achieved accuracy:\\n{{}}".format(e.accuracy))
    ...   raise e
        Indices of diverging points: [1]