<commit_before>import numpy as np
import scipy.integrate
import scipy.optimize

DEG2RAD = np.pi/180

def _toArray(x):
    """Convert x to array if needed

    Returns:
        array(x), boolean if x was an array before
    """
    if hasattr(x, '__iter__'):
        return np.array(x), True
    return np.array([x], dtype=np.double), False

def ellipticity(a, b):
    """Returns 1-abs(b/a)"""
    return 1-np.abs(b/a)
def meanDistortion(a, b):
    """Returns average `ellipticity` over all `a`,`b`"""
    return np.mean(ellipticity(a,b))
def maxDistortion(a,b):
    """Returns max `ellipticity` over all `a`,`b`"""
    return np.max(ellipticity(a,b))
def stdDistortion(a,b):
    """Returns `std(b/a)`"""
    return (b/a).std() # include the sign
def stdScale(a,b):
    """Returns `std(a*b)`

    This is useful for conformal projections.
    """
    return (a*b).std()
def stdDistortionScale(a,b):
    """Retruns sum of `stdScale` and `stdDistortion`.

    This is useful for a compromise between equal-area and conformal projections.
    """
    return stdScale(a,b) + stdDistortion(a,b)

def _optimize_objective(x, proj_cls, lon_type, lon, lat, crit):
    """Construct projections from parameters `x` and compute `crit` for `lon, lat`"""
    proj = proj_cls(*x, lon_type=lon_type)
    a, b = proj.distortion(lon, lat)
    return crit(a,b)

def _optimize(proj_cls, x0, lon_type, lon, lat, crit, bounds=None):
    """Determine parameters for `proj_cls` that minimize `crit` over `lon, lat`.

    Args:
        proj_cls: projection class
        x0: arguments for projection class `__init__`
        lon_type: type of longitude, "lon" or "ra" (see `BaseProjection`)
        lon: list of rectascensions
        lat: list of declinations
        crit: optimization criterion
            needs to be function of semi-major and semi-minor axes of the Tissot indicatix
        bounds: list of upper and lower bounds on each parameter in `x0`

    Returns:
        optimized projection of class `proj_cls`
    """
    print ("optimizing parameters of %s to minimize %s" % (proj_cls.__name__, crit.__name__))
    x, fmin, d = scipy.optimize.fmin_l_bfgs_b(_optimize_objective, x0, args=(proj_cls, lon_type, lon, lat, crit), bounds=bounds, approx_grad=True)
    res = proj_cls(*x, lon_type=lon_type)
    print ("best objective %.6f at %r" % (fmin, res))
    return res

def _dist(radec, proj, xy):
    return np.sum((xy - np.array(proj(radec[0], radec[1])))**2)


class BaseProjection(object):
    """Projection base class

    Every projection needs to implement three methods:
    * `transform(self, lon, lat)`: mapping from lon/lat to map x/y
    * `invert(self, x, y)`: the inverse mapping from x/y to lon/lat

    All methods accept either single number or arrays and return accordingly.
    """
    def __init__(self, lon_0=0, lon_type="ra"):
        """Initialize projection

        Args:
            lon_0 (int, float):  reference longitude
            lon_type (string): type of longitude
                "lon" for a standard coordinate system (right-handed, -180..180 deg)
                "ra" for n equatorial coordinate system (left-handed, 0..360 deg)
        """
        assert lon_type in ['ra', 'lon']
        self.lon_0 = lon_0
        self.lon_type = lon_type
        if self.lon_type == "ra" and self.lon_0 < 0:
            self.lon_0 += 360
        elif self.lon_type == "lon" and self.lon_0 > 180:
            self.lon_0 -= 360

    def __call__(self, lon, lat):
        """Shorthand for `transform`. Works also with single coordinates

        Args:
            lon (float, array): longitude
            lat (float, array): latitude
        """
        lon_, isArray = _toArray(lon)
        lat_, isArray = _toArray(lat)
        assert len(lon_) == len(lat_)

        x, y = self.transform(lon_, lat_)
        if isArray:
            return x, y
        else:
            return x[0], y[0]

    def transform(self, lon, lat):
        """Convert longitude/latitude into map coordinates

        Note: Unlike `__call__`, lon/lat need to be arrays!

        Args:
            lon (array): longitudes
            lat (array): latitudes

        Returns:
            x,y with the same format as lon,lat
        """
        pass

    def inv(self, x, y):
        """Shorthand for `invert`. Works also with single coordinates

        Args:
            x (float, array): horizontal map coordinate
            y (float, array): vertical map coordinate
        """
        x_, isArray = _toArray(x)
        y_, isArray = _toArray(y)
        lon, lat = self.invert(x_, y_)
        if isArray:
            return lon, lat
        else:
            return lon[0], lat[0]

    def invert(self, x, y):
        """Convert map coordinates into longitude/latitude

        Args:
            x (array): horizontal map coordinates
            y (array): vertical map coordinates

        Returns:
            lon,lat with the same format as x,y
        """
        # default implementation for non-analytic inverses
        assert len(x) == len(y)

        bounds = ((None,None), (-90, 90)) # lon/lat limits
        start = (self.lon_0,0) # lon/lat of initial guess: should be close to map center
        lon, lat = np.empty(len(x)), np.empty(len(y))
        i = 0
        for x_,y_ in zip(x, y):
            xy = np.array([x_,y_])
            radec, fmin, d = scipy.optimize.fmin_l_bfgs_b(_dist, start, args=(self, xy), bounds=bounds, approx_grad=True)
            if fmin < 1e-6: # smaller than default tolerance of fmin
                lon[i], lat[i] = radec
            else:
                lon[i], lat[i] = -1000, -1000
            i += 1

        return lon, lat

    @property
    def poleIsPoint(self):
        """Whether the pole is mapped onto a point"""
        try:
            return self._poleIsPoint
        except AttributeError:
            self._poleIsPoint = {}
            N = 10
            # run along the poles from the left to right outer meridian
            rnd_meridian = -180 + 360*np.random.rand(N) + self.lon_0
            for deg in [-90, 90]:
                line = self.transform(rnd_meridian, deg*np.ones(N))
                if np.unique(line[0]).size > 1 and np.unique(line[1]).size > 1:<commit_msg>fixed bug that prevented meridian labels at +-90<commit_after>import numpy as np
import scipy.integrate
import scipy.optimize

DEG2RAD = np.pi/180

def _toArray(x):
    """Convert x to array if needed

    Returns:
        array(x), boolean if x was an array before
    """
    if hasattr(x, '__iter__'):
        return np.array(x), True
    return np.array([x], dtype=np.double), False

def ellipticity(a, b):
    """Returns 1-abs(b/a)"""
    return 1-np.abs(b/a)
def meanDistortion(a, b):
    """Returns average `ellipticity` over all `a`,`b`"""
    return np.mean(ellipticity(a,b))
def maxDistortion(a,b):
    """Returns max `ellipticity` over all `a`,`b`"""
    return np.max(ellipticity(a,b))
def stdDistortion(a,b):
    """Returns `std(b/a)`"""
    return (b/a).std() # include the sign
def stdScale(a,b):
    """Returns `std(a*b)`

    This is useful for conformal projections.
    """
    return (a*b).std()
def stdDistortionScale(a,b):
    """Retruns sum of `stdScale` and `stdDistortion`.

    This is useful for a compromise between equal-area and conformal projections.
    """
    return stdScale(a,b) + stdDistortion(a,b)

def _optimize_objective(x, proj_cls, lon_type, lon, lat, crit):
    """Construct projections from parameters `x` and compute `crit` for `lon, lat`"""
    proj = proj_cls(*x, lon_type=lon_type)
    a, b = proj.distortion(lon, lat)
    return crit(a,b)

def _optimize(proj_cls, x0, lon_type, lon, lat, crit, bounds=None):
    """Determine parameters for `proj_cls` that minimize `crit` over `lon, lat`.

    Args:
        proj_cls: projection class
        x0: arguments for projection class `__init__`
        lon_type: type of longitude, "lon" or "ra" (see `BaseProjection`)
        lon: list of rectascensions
        lat: list of declinations
        crit: optimization criterion
            needs to be function of semi-major and semi-minor axes of the Tissot indicatix
        bounds: list of upper and lower bounds on each parameter in `x0`

    Returns:
        optimized projection of class `proj_cls`
    """
    print ("optimizing parameters of %s to minimize %s" % (proj_cls.__name__, crit.__name__))
    x, fmin, d = scipy.optimize.fmin_l_bfgs_b(_optimize_objective, x0, args=(proj_cls, lon_type, lon, lat, crit), bounds=bounds, approx_grad=True)
    res = proj_cls(*x, lon_type=lon_type)
    print ("best objective %.6f at %r" % (fmin, res))
    return res

def _dist(radec, proj, xy):
    return np.sum((xy - np.array(proj(radec[0], radec[1])))**2)


class BaseProjection(object):
    """Projection base class

    Every projection needs to implement three methods:
    * `transform(self, lon, lat)`: mapping from lon/lat to map x/y
    * `invert(self, x, y)`: the inverse mapping from x/y to lon/lat

    All methods accept either single number or arrays and return accordingly.
    """
    def __init__(self, lon_0=0, lon_type="ra"):
        """Initialize projection

        Args:
            lon_0 (int, float):  reference longitude
            lon_type (string): type of longitude
                "lon" for a standard coordinate system (right-handed, -180..180 deg)
                "ra" for n equatorial coordinate system (left-handed, 0..360 deg)
        """
        assert lon_type in ['ra', 'lon']
        self.lon_0 = lon_0
        self.lon_type = lon_type
        if self.lon_type == "ra" and self.lon_0 < 0:
            self.lon_0 += 360
        elif self.lon_type == "lon" and self.lon_0 > 180:
            self.lon_0 -= 360

    def __call__(self, lon, lat):
        """Shorthand for `transform`. Works also with single coordinates

        Args:
            lon (float, array): longitude
            lat (float, array): latitude
        """
        lon_, isArray = _toArray(lon)
        lat_, isArray = _toArray(lat)
        assert len(lon_) == len(lat_)

        x, y = self.transform(lon_, lat_)
        if isArray:
            return x, y
        else:
            return x[0], y[0]

    def transform(self, lon, lat):
        """Convert longitude/latitude into map coordinates

        Note: Unlike `__call__`, lon/lat need to be arrays!

        Args:
            lon (array): longitudes
            lat (array): latitudes

        Returns:
            x,y with the same format as lon,lat
        """
        pass

    def inv(self, x, y):
        """Shorthand for `invert`. Works also with single coordinates

        Args:
            x (float, array): horizontal map coordinate
            y (float, array): vertical map coordinate
        """
        x_, isArray = _toArray(x)
        y_, isArray = _toArray(y)
        lon, lat = self.invert(x_, y_)
        if isArray:
            return lon, lat
        else:
            return lon[0], lat[0]

    def invert(self, x, y):
        """Convert map coordinates into longitude/latitude

        Args:
            x (array): horizontal map coordinates
            y (array): vertical map coordinates

        Returns:
            lon,lat with the same format as x,y
        """
        # default implementation for non-analytic inverses
        assert len(x) == len(y)

        bounds = ((None,None), (-90, 90)) # lon/lat limits
        start = (self.lon_0,0) # lon/lat of initial guess: should be close to map center
        lon, lat = np.empty(len(x)), np.empty(len(y))
        i = 0
        for x_,y_ in zip(x, y):
            xy = np.array([x_,y_])
            radec, fmin, d = scipy.optimize.fmin_l_bfgs_b(_dist, start, args=(self, xy), bounds=bounds, approx_grad=True)
            if fmin < 1e-6: # smaller than default tolerance of fmin
                lon[i], lat[i] = radec
            else:
                lon[i], lat[i] = -1000, -1000
            i += 1

        return lon, lat

    @property
    def poleIsPoint(self):
        """Whether the pole is mapped onto a point"""
        try:
            return self._poleIsPoint
        except AttributeError:
            self._poleIsPoint = {}
            N = 10
            # run along the poles from the left to right outer meridian
            rnd_meridian = -180 + 360*np.random.rand(N) + self.lon_0
            for deg in [-90, 90]:
                line = self.transform(rnd_meridian, deg*np.ones(N))
                if np.unique(line[0]).size > 1 or np.unique(line[1]).size > 1:
{'fraction_non_alphanumeric': 0.11767531219980788, 'fraction_numerical': 0.013288504642971502, 'mean_word_length': 2.2764227642276422, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1038117', 'n_tokens_mistral': 4089, 'n_tokens_neox': 3858, 'n_words': 1256}