public class Position implements Serializable {
    private static final long  serialVersionUID        = -4364498345989515044L;



    public static final double WGS84_EQUATORIAL_RADIUS = 6378137.0;        // ellipsoid equatorial getRadius, in meters
    public static final double WGS84_POLAR_RADIUS      = 6356752.3;        // ellipsoid polar getRadius, in meters
    public static final double WGS84_ES                = 0.00669437999013; // eccentricity squared, semi-major axis



    private boolean            _baseFormIsPolar;
    private Vector             _polar;     // Longutude, Latitude [, Elevation]
    private Vector             _cartesian; // X, Y, Z
    private boolean            _elevationUsed;

    /**
     * Create an object that represents a position on the surface of the Earth
     * using standard latitude and longitude. The resultant position will ignore
     * elevation.
     * 
     * @param longitude
     *            The longitude of the position, in degrees
     * @param latitude
     *            The latitude of the position, in degrees
     */
    public Position (double longitude, double latitude) {
        _polar = new Vector(AngleUtilities.intoRangeDegrees(0.0, longitude),
                            latitude);
        _cartesian = null;
        _elevationUsed = false;
        _baseFormIsPolar = true;
    }

    /**
     * Create an object that represents a position on the surface of the Earth
     * using standard latitude, longitude, and elevation.
     * 
     * @param longitude
     *            The longitude of the position, in degrees
     * @param latitude
     *            The latitude of the position, in degrees
     * @param elevation
     *            The elevation of the position, in meters above sea level
     */
    public Position (double longitude, double latitude, double elevation) {
        _polar = new Vector(AngleUtilities.intoRangeDegrees(0.0, longitude),
                            latitude, elevation);
        _cartesian = null;
        _elevationUsed = true;
        _baseFormIsPolar = true;
    }

    /**
     * Create an object that represents a position on the surface of the Earth,
     * using cartesian coordinates with the Y axis being the axis of rotation,
     * the XZ plane being the equator, with the Z axis pointed towards the
     * Greenwich meridian.
     * 
     * @param x
     *            The coordinate in the direction of (0 deg N, 90 deg E)
     * @param y
     *            The coordinate in the direction of the north pole
     * @param z
     *            The coordinate in the direction of (0 deg N, 0 deg E)
     * @param ignoreElevation
     *            if true, sets calculations using this position to ignore any
     *            elevation inherent in the position; if false, elevation will
     *            be taken into account
     */
    public Position (double x, double y, double z, boolean ignoreElevation) {
        _cartesian = new Vector(x, y, z);
        _polar = null;
        _elevationUsed = !ignoreElevation;
        _baseFormIsPolar = false;
    }

    /**
     * Copy a position, possibly changing from polar to cartesian or vice versa
     * 
     * @param oldForm
     *            The position being copied
     * @param polar
     *            True if the new position should be in polar coordinates, false
     *            if it should be in cartesian coordinates.
     */
    public Position (Position oldForm, boolean polar) {
        if (null == oldForm._polar) {
            _polar = null;
        } else {
            _polar = new Vector(oldForm._polar);
        }

        if (null == oldForm._cartesian) {
            _cartesian = null;
        } else {
            _cartesian = new Vector(oldForm._cartesian);
        }

        _elevationUsed = oldForm._elevationUsed;

        if (_baseFormIsPolar) updatePolar();
        else updateCartesian();
    }


    // ////////////////////////////////////////////////////////////////////////
    // Section: Conversion between polar and cartesian forms
    //
    private void updatePolar () {
        if (null != _polar) return;
        if (null == _cartesian) return;

        // According to
        // H. Vermeille, 
        // "An analytical method to transform geocentric into geodetic coordinates"
        // http://www.springerlink.com/content/3t6837t27t351227/fulltext.pdf
        // Journal of Geodesy, accepted 10/2010, not yet published

        // Note coordinates for this are rotation axis = Z, not rotation axis = Y.
        double X = _cartesian.coord(2);
        double Y = _cartesian.coord(0);
        double Z = _cartesian.coord(1);
        double XXpYY = X*X+Y*Y;
        double sqrtXXpYY = Math.sqrt(XXpYY);

        double a = WGS84_EQUATORIAL_RADIUS;
        double ra2 = 1/(a*a);
        double e2 = WGS84_ES;
        double e4 = e2*e2;

        // Step 1
        double p = XXpYY*ra2;
        double q = Z*Z*(1-e2)*ra2;
        double r = (p+q-e4)/6;

        double h;
        double phi;

        double evoluteBorderTest = 8*r*r*r+e4*p*q;
        if (evoluteBorderTest > 0 || q != 0) {
            double u;

            if (evoluteBorderTest > 0) {
                // Step 2: general case
                double rad1 = Math.sqrt(evoluteBorderTest);
                double rad2 = Math.sqrt(e4*p*q);
    
                // 10*e2 is my arbitrary decision of what Vermeille means by "near... the cusps of the evolute".
                if (evoluteBorderTest > 10*e2) {
                    double rad3 = Math.cbrt((rad1+rad2)*(rad1+rad2));
                    u = r + 0.5*rad3 + 2*r*r/rad3;
                } else {
                    u = r + 0.5*Math.cbrt((rad1+rad2)*(rad1+rad2))+0.5*Math.cbrt((rad1-rad2)*(rad1-rad2));
                }
            } else {
                // Step 3: near evolute
                double rad1 = Math.sqrt(-evoluteBorderTest);
                double rad2 = Math.sqrt(-8*r*r*r);
                double rad3 = Math.sqrt(e4*p*q);
                double atan = 2*Math.atan2(rad3, rad1+rad2)/3;
    
                u = -4*r*Math.sin(atan)*Math.cos(Math.PI/6+atan);
            }
    
            double v = Math.sqrt(u*u+e4*q);
            double w = e2*(u+v-q)/(2*v);
            double k = (u+v)/(Math.sqrt(w*w+u+v)+w);
            double D = k*sqrtXXpYY/(k+e2);
            double sqrtDDpZZ = Math.sqrt(D*D+Z*Z);

            h = (k+e2-1)*sqrtDDpZZ/k;
            phi = 2*Math.atan2(Z, sqrtDDpZZ+D);
        } else {
            // Step 4: singular disk
            double rad1 = Math.sqrt(1-e2);
            double rad2 = Math.sqrt(e2-p);
            double e = Math.sqrt(e2);

            h = -a*rad1*rad2/e;
            phi = rad2/(e*rad2+rad1*Math.sqrt(p));
        }

        // Compute lambda
        double lambda;
        double s2 = Math.sqrt(2);
        if ((s2-1)*Y < sqrtXXpYY+X) {
            // case 1 - -135deg < lambda < 135deg
            lambda = 2*Math.atan2(Y, sqrtXXpYY+X);
        } else if (sqrtXXpYY+Y < (s2+1)*X) {
            // case 2 - -225deg < lambda < 45deg
            lambda = -Math.PI*0.5+2*Math.atan2(X, sqrtXXpYY-Y);
        } else {
            // if (sqrtXXpYY-Y<(s2=1)*X) {  // is the test, if needed, but it's not
            // case 3: - -45deg < lambda < 225deg
            lambda = Math.PI*0.5 - 2*Math.atan2(X, sqrtXXpYY+Y);
        }

        if (_elevationUsed) {
            _polar = new Vector(AngleUtilities.intoRangeDegrees(0.0, Math.toDegrees(lambda)), Math.toDegrees(phi), h);
        } else {
            _polar = new Vector(AngleUtilities.intoRangeDegrees(0.0, Math.toDegrees(lambda)), Math.toDegrees(phi));
        }

        _polar.setPrecision(_cartesian.getPrecision());
    }

    private void updateCartesian () {
        if (null != _cartesian) return;
        if (null == _polar) return;

        // from WorldWind: ElipsoidalGlobe.geodeticToCartesian
        double latitude = getLatitudeRadians();
        double cosLat = Math.cos(latitude);
        double sinLat = Math.sin(latitude);
        double longitude = getLongitudeRadians();
        double cosLon = Math.cos(longitude);
        double sinLon = Math.sin(longitude);
        double elevation = (_elevationUsed ? getElevation() : 0.0);

        // Get the radius (in meters) of the vertical in prime meridian
        double rpm = WGS84_EQUATORIAL_RADIUS / Math.sqrt(1.0 - WGS84_ES * sinLat * sinLat);

        // Y axis is the axis of rotation, just like in WorldWind.
        double x = (rpm + elevation) * cosLat * sinLon;
        double z = (rpm + elevation) * cosLat * cosLon;
        double y = (rpm * (1.0 - WGS84_ES) + elevation) * sinLat;

        _cartesian = new Vector(x, y, z);
        _cartesian.setPrecision(_polar.getPrecision());
    }



    // ////////////////////////////////////////////////////////////////////////
    // Section: accessors
    //
    public void setPrecision (double precision) {
        if (null != _polar) _polar.setPrecision(precision);
        if (null != _cartesian) _cartesian.setPrecision(precision);
    }

    public double getPrecision () {
        if (_baseFormIsPolar) return _polar.getPrecision();
        else return _cartesian.getPrecision();
    }

    public double getLongitude () {
        updatePolar();
        return _polar.coord(0);
    }

    public double getLongitudeRadians () {
        return Math.toRadians(getLongitude());
    }

    public double getLatitude () {
        updatePolar();
        return _polar.coord(1);
    }

    public double getLatitudeRadians () {
        return Math.toRadians(getLatitude());
    }

    /**
     * Gets the elevation of this position.
     * 
     * @return The elevation of this position in meters
     * @throws UnsupportedOperationException
     *             if this position doesn't specify an elevation.
     */
    public double getElevation () {
        if (!_elevationUsed)
            throw new UnsupportedOperationException("Attempt to get elevation of elevationless position");
        updatePolar();
        return _polar.coord(2);
    }

    public boolean hasElevation () {
        return _elevationUsed;
    }

    public Vector getAsCartesian () {
        updateCartesian();
        return _cartesian;
    }



    // ////////////////////////////////////////////////////////////////////////
    // Section: Calculations
    //
    public double getDistanceMeters (Position p) {
        double dR = Math.toRadians(getAngularDistance(p));
        // estimate radius based on average latitude
        double radius = getEarthRadius((getLatitude()+p.getLatitude())/2.0);
        return dR*radius;
    }
    /**
     * Gets the angular distance between this point and the other, in degrees
     * 
     * @param p
     *            The other point
     * @return The angular distance, in degrees, between this and p
     */
    public double getAngularDistance (Position p) {
        double lon1 = getLongitudeRadians();
        double lat1 = getLatitudeRadians();
        double lon2 = AngleUtilities.intoRangeDegrees(lon1, p.getLongitudeRadians());
        double lat2 = p.getLatitudeRadians();

        // Taken from WorldWind's LatLon class:
        double epsilon = getPrecision();
        if (Math.abs(lon1 - lon2) < epsilon && Math.abs(lat1 - lat2) < epsilon)
            return 0.0;

        // Taken from "Map Projections - A Working Manual", page 30, equation
        // 5-3a.
        // The traditional d=2*asin(a) form has been replaced with
        // d=2*atan2(sqrt(a), sqrt(1-a))
        // to reduce rounding errors with large distances.
        double a = Math.sin((lat2 - lat1) / 2.0)
                   * Math.sin((lat2 - lat1) / 2.0) + Math.cos(lat1)
                   * Math.cos(lat2) * Math.sin((lon2 - lon1) / 2.0)
                   * Math.sin((lon2 - lon1) / 2.0);
        double distanceRadians = 2.0 * Math.atan2(Math.sqrt(a),
                                                  Math.sqrt(1 - a));

        if (Double.isNaN(distanceRadians))
            return 0.0;
        else
            return Math.toDegrees(distanceRadians);
    }

    public double getAngularDistanceSquared (Position other) {
        double d = getAngularDistance(other);
        return d*d;
    }

    public double getCartesianDistance (Position p) {
        return getAsCartesian().getDistance(p.getAsCartesian());
    }

    public double getCartesianDistanceSquared (Position p) {
        return getAsCartesian().getDistanceSquared(p.getAsCartesian());
    }

    /**
     * Gets the angle from this point to the other.
     * 
     * @param p
     *            The other point.
     * @return The angle from this to p, in degree, with north=0, east=90
     */
    public double getAzimuth (Position p) {
        double lon1 = getLongitudeRadians();
        double lat1 = getLatitudeRadians();
        double lon2 = AngleUtilities.intoRangeRadians(lon1, p.getLongitudeRadians());
        double lat2 = p.getLatitudeRadians();

        // Taken from WorldWind's LatLon class:
        double epsilon = getPrecision();
        if (Math.abs(lon1-lon2) < epsilon && Math.abs(lat1-lat2) < epsilon)
            return 0.0;

        if (Math.abs(lon1-lon2) < epsilon) {
            if (lat1 > lat2) return 180.0;
            else return 0.0;
        }

        // Taken from "Map Projections - A Working Manual", page 30, equation 5-4b.
        // The atan2() function is used in place of the traditional atan(y/x) to simplify the case when x==0.
        double y = Math.cos(lat2) * Math.sin(lon2 - lon1);
        double x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1)
                   * Math.cos(lat2) * Math.cos(lon2 - lon1);
        double azimuthRadians = Math.atan2(y, x);

        if (Double.isNaN(azimuthRadians))
            return 0.0;
        else
            return Math.toDegrees(azimuthRadians);
    }

    /**
     * Get the point offset from this point by the given distance at the given
     * azimuth.
     * 
     * @param azimuth
     *            The angle, in degree, with N=0, E=90, in which the desired
     *            offset point lies
     * @param distance
     *            The angular distance, in degrees, by which to offset this
     *            point.
     * @return The point offset from this one as described.
     */
    public Position offset (double azimuthDegrees, double distanceDegrees) {
        double lon = getLongitudeRadians();
        double lat = getLatitudeRadians();
        double azimuth = Math.toRadians(azimuthDegrees);
        double distance = Math.toRadians(distanceDegrees);

        // Taken from WorldWind's LatLon class:
        if (distance == 0)
            return this;

        // Taken from "Map Projections - A Working Manual", page 31, equation 5-5 and 5-6.
        double endLatRadians = Math.asin(Math.sin(lat) * Math.cos(distance)
                                         + Math.cos(lat) * Math.sin(distance)
                                         * Math.cos(azimuth));
        double endLonRadians = lon
                               + Math.atan2(Math.sin(distance)
                                                    * Math.sin(azimuth),
                                            Math.cos(lat) * Math.cos(distance)
                                                    - Math.sin(lat)
                                                    * Math.sin(distance)
                                                    * Math.cos(azimuth));

        if (Double.isNaN(endLatRadians) || Double.isNaN(endLonRadians))
            return this;

        Position p = new Position(Math.toDegrees(endLonRadians),
                                  Math.toDegrees(endLatRadians));
        p.setPrecision(getPrecision());
        return p;
    }

    /**
     * Get the point on the great circle from a to b at the given latitude
     * 
     * @param a
     *            The segment start point
     * @param b
     *            The segment end point
     * @param latitude
     *            The latitude of the desired result (in degrees)
     * @return A point on the great circle between A and B at the given
     *         latitude, or null if the shortest-path segment of the great
     *         circle from a to b doesn't intersect that latitude. If the given
     *         latitude is intersected multiple times, the first is returned.
     */
    public static Position greatCircleAtLatitude (Position p1, Position p2, double latitude) {
        // from http://williams.best.vwh.net/avform.htm#Par
        double lon1 = p1.getLongitudeRadians();
        double lat1 = p1.getLatitudeRadians();
        double lon2 = AngleUtilities.intoRangeRadians(lon1, p2.getLongitudeRadians());
        double lat2 = p2.getLatitudeRadians();

        double lat3 = Math.toRadians(latitude);

        double sinlat1 = Math.sin(lat1);
        double coslat1 = Math.cos(lat1);
        double sinlat2 = Math.sin(lat2);
        double coslat2 = Math.cos(lat2);
        double sinlat3 = Math.sin(lat3);
        double coslat3 = Math.cos(lat3);

        double dl = lon1-lon2;
        double sindl = Math.sin(dl);
        double cosdl = Math.cos(dl);

        double A = sinlat1 * coslat2 * coslat3 * sindl;
        double B = sinlat1 * coslat2 * coslat3 * cosdl - coslat1 * sinlat2 * coslat3;
        double C = coslat1 * coslat2 * sinlat3 * sindl;

        double baseLon = Math.atan2(B,  A);

        if (C*C > (A*A + B*B)) {
            // No crossing at that latitude
            return null;
        } else {
            double dlon = Math.acos(C/Math.sqrt(A*A+B*B));
            double centerLon = (lon1+lon2)/2.0;
            double targetLon1 = AngleUtilities.intoRangeRadians(centerLon, (lon1+dlon+baseLon));
            Position r1 = new Position(Math.toDegrees(targetLon1), latitude);
            double targetLon2 = AngleUtilities.intoRangeRadians(centerLon, (lon1-dlon+baseLon));
            Position r2 = new Position(Math.toDegrees(targetLon2), latitude);

            boolean good1, good2;
            if (lon1 < lon2) {
                good1 = (lon1 <= targetLon1 && targetLon1 <= lon2);
                good2 = (lon1 <= targetLon2 && targetLon2 <= lon2);
            } else {
                good1 = (lon2 <= targetLon1 && targetLon1 <= lon1);
                good2 = (lon2 <= targetLon2 && targetLon2 <= lon1);
            }

            if (good1 && good2) {
                double d1 = p1.getAngularDistance(r1);
                double d2 = p1.getAngularDistance(r2);
                if (d1 < d2) {
                    return r1;
                } else {
                    return r2;
                }
            } else if (good1) {
                return r1;
            } else if (good2) {
                return r2;
            } else {
                return null;
            }
        }
    }

    public static Position greatCircleAtLongitude (Position a, Position b,
                                                  double longitude) {
        double abAz = a.getAzimuth(b);
        double azTan = Math.tan(Math.toRadians(abAz));
        double latA = a.getLatitudeRadians();
        double dLon = Math.toRadians(longitude)-a.getLongitudeRadians();
        double y = Math.sin(dLon) + azTan*Math.sin(latA)*Math.cos(dLon);
        double x = azTan*Math.cos(latA);
        double latitudeRad = Math.atan2(y, x);
        double latitude = Math.toDegrees(latitudeRad);
        return new Position(longitude, latitude);
    }

    public static Rectangle2D getGreatCircleLonLatBounds (Position p1, Position p2) {
        double epsilon = p1.getPrecision();

        // from http://williams.best.vwh.net/avform.htm#Par
        double lon1 = p1.getLongitude();
        double lat1 = p1.getLatitude();
        double lon2 = AngleUtilities.intoRangeDegrees(lon1, p2.getLongitude());
        double lat2 = p2.getLatitude();

        // Check for the case where the two points define the equator
        if (Math.abs(lat1) < epsilon && Math.abs(lat2) < epsilon) {
            double minLon = Math.min(lon1, lon2);
            double maxLon = Math.max(lon1, lon2);
            return new Rectangle2D.Double(minLon, 0.0, maxLon-minLon, 0.0);
        }

        // Not the equator; therefore it crosses the equator. Figure out where -
        // max and min will be half way in between
        //
        // The rest of this code is adapted from  greatCircleAtLatitude 
        // (i.e., http://williams.best.vwh.net/avform.htm#Par)
        lon1 = Math.toRadians(lon1);
        lat1 = Math.toRadians(lat1);
        lon2 = Math.toRadians(lon2);
        lat2 = Math.toRadians(lat2);
        double lat3 = 0.0;

        double sinlat1 = Math.sin(lat1);
        double coslat1 = Math.cos(lat1);
        double sinlat2 = Math.sin(lat2);
        double coslat2 = Math.cos(lat2);
        double sinlat3 = Math.sin(lat3);
        double coslat3 = Math.cos(lat3);

        double dl = lon1-lon2;
        double sindl = Math.sin(dl);
        double cosdl = Math.cos(dl);

        double A = sinlat1 * coslat2 * coslat3 * sindl;
        double B = sinlat1 * coslat2 * coslat3 * cosdl - coslat1 * sinlat2 * coslat3;
        double C = coslat1 * coslat2 * sinlat3 * sindl;
        assert(0.0 == C);

        // C should be 0, therefore dlon should be 0. lon1+baseLon should
        // therefore be an equatorial crossing point. The max is therefore PI/2
        // away in the direction of either point.
        double baseLon = Math.atan2(B,  A);
        Position equatorialCrossing = new Position(Math.toDegrees(lon1 + baseLon - Math.PI/2.0), 0.0);
        double direction = equatorialCrossing.getAzimuth(p1);
        Position maxPt, minPt;
        if (direction > 0) {
            maxPt = equatorialCrossing.offset(-direction, 90);
            minPt = equatorialCrossing.offset(-direction, -90);
        } else {
            maxPt = equatorialCrossing.offset(direction, 90);
            minPt = equatorialCrossing.offset(direction, -90);
        }

        double minLon = Math.min(p1.getLongitude(), p2.getLongitude());
        double maxLon = Math.max(AngleUtilities.intoRangeDegrees(minLon, p1.getLongitude()),
                                 AngleUtilities.intoRangeDegrees(minLon, p2.getLongitude()));
        double minLat = Math.min(p1.getLatitude(), p2.getLatitude());
        double maxLat = Math.max(p1.getLatitude(), p2.getLatitude());

        double centerLon = (minLon+maxLon)/2.0;
        double minPtLon = AngleUtilities.intoRangeDegrees(centerLon, minPt.getLongitude());
        if (minLon < minPtLon && minPtLon < maxLon) {
            minLat = minPt.getLatitude();
        }
        double maxPtLon = AngleUtilities.intoRangeDegrees(centerLon, maxPt.getLongitude());
        if (minLon < maxPtLon && maxPtLon < maxLon) {
            maxLat = maxPt.getLatitude();
        }
        return new Rectangle2D.Double(minLon, minLat, maxLon-minLon, maxLat-minLat);
    }



    @Override
    public boolean equals (Object obj) {
        if (this == obj) return true;
        if (null == obj) return false;
        if (!(obj instanceof Position)) return false;
        Position p = (Position) obj;

        if (_baseFormIsPolar) {
            double epsilon = getPrecision();
            if (Math.abs(getLongitude() - p.getLongitude()) >= epsilon)
                return false;
            if (Math.abs(getLatitude() - p.getLatitude()) >= epsilon)
                return false;
            if (_elevationUsed) {
                if (!p._elevationUsed)
                    return false;
                if (Math.abs(getElevation() - p.getElevation()) >= epsilon)
                    return false;
            }
        } else {
            if (!_cartesian.equals(p.getAsCartesian())) return false;
        }
        return true;
    }

    @Override
    public String toString () {
        if (_baseFormIsPolar) {
            double lon = getLongitude();
            String lonMarker = "E";
            if (lon < 0) {
                lon = -lon;
                lonMarker = "W";
            }
            int lonDegrees = (int) Math.floor(lon);
            int lonMinutes = (int) Math.floor((lon-lonDegrees)*60.0);
            int milliLonSeconds = (int) Math.round(((lon-lonDegrees)*60.0-lonMinutes)*60000.0);
            if (60000 == milliLonSeconds) {
                milliLonSeconds = 0;
                lonMinutes++;
                if (60 <= lonMinutes) {
                    lonMinutes -= 60;
                    ++lonDegrees;
                }
            }
    
            double lat = getLatitude();
            String latMarker = "N";
            if (lat < 0) {
                lat = -lat;
                latMarker = "S";
            }
            int latDegrees = (int) Math.floor(lat);
            int latMinutes = (int) Math.floor((lat-latDegrees)*60.0);
            int milliLatSeconds = (int) Math.round(((lat-latDegrees)*60.0-latMinutes)*60000.0);
            if (60000 == milliLatSeconds) {
                milliLatSeconds = 0;
                latMinutes++;
                if (60 <= latMinutes) {
                    latMinutes -= 60;
                    ++latDegrees;
                }
            }
    
            return String.format("%d%d'%.3f\"%s, %d%d'%.3f\"%s",
                                 lonDegrees, lonMinutes, milliLonSeconds/1000.0, lonMarker,
                                 latDegrees, latMinutes, milliLatSeconds/1000.0, latMarker);
        } else {
            return _cartesian.toString();
        }
    }

    /**
     * Get the radius of the earth at this point.
     */
    public double getEarthRadius () {
        return getEarthRadius(getLatitude());
    }

    public static double getEarthRadius (double latitude) {
        // From WorldWind, EllipsoidalGlobe.getRadiusAt(...)
        return new Position(0.0, latitude).getAsCartesian().vectorLength();
    }
}
 public class PositionCalculationParameters {
    private PositionCalculationType _type;
    private double                  _allowedError;
    private double                  _epsilon;
    private boolean                 _ignoreDirection;

    /**
     * Create a set of parameters that determine how geodetic calculations are
     * performed
     * 
     * @param type
     *            The way in which the calculations are to be performed (i.e.,
     *            using what system and approcimation)
     * @param allowedError
     *            The proportion of error allowed in approximations
     * @param precision
     *            The amount of error allowed in equality tests
     */
    public PositionCalculationParameters (PositionCalculationType type,
                                          double allowedError, double precision,
                                          boolean ignoreDirection) {
        _type = type;
        _allowedError = allowedError;
        _epsilon = precision;
        _ignoreDirection = ignoreDirection;
    }

    /**
     * Get the way in which the calculations are to be performed (i.e., using
     * what system and approcimation)
     * 
     * @see PositionCalculationType
     */
    public PositionCalculationType getCalculationType () {
        return _type;
    }

    /**
     * Get the proportion of error allowed in approximations
     */
    public double getAllowedError () {
        return _allowedError;
    }

    /**
     * Get the amount of error allowed in equality tests
     */
    public double getPrecision () {
        return _epsilon;
    }

    /**
     * If true, the direction in which a track goes should be ignored when
     * calculating distance between tracks.
     */
    public boolean ignoreDirection () {
        return _ignoreDirection;
    }
} abstract public class Track {
    private PositionCalculationParameters _parameters;
    private List<Position>                _points;
    private List<Double>                  _parameterization;
    private double                        _length;
    private Track                         _reverse;
    private Map<String, Double>           _statistics;

    protected Track (PositionCalculationParameters parameters, Position... points) {
        this(parameters, Arrays.asList(points));
    }

    protected Track (PositionCalculationParameters parameters,
                     List<Position> points) {
        _parameters = parameters;
        _points = new ArrayList<Position>(points);
        _reverse = null;
        _statistics = null;

        calculateLengthParameterization();

        reduce();
    }

    protected Track (PositionCalculationParameters parameters,
                     List<Position> points,
                     List<Double> parameterization) {
        _parameters = parameters;
        _points = new ArrayList<Position>(points);
        _parameterization = new ArrayList<Double>(parameterization);
        _reverse = null;
        _statistics = null;

        calculateLength();

        reduce();
    }

    protected Track (Track oldTrack, PositionCalculationParameters newParameters) {
        _parameters = newParameters;
        _points = new ArrayList<Position>();
        _parameterization = new ArrayList<Double>();
        _reverse = null;
        _statistics = null;

        for (int i=0; i<oldTrack._points.size(); ++i) {
            Position p = new Position(oldTrack._points.get(i),
                                      !_parameters.getCalculationType().equals(PositionCalculationType.Cartesian3D));
            p.setPrecision(_parameters.getPrecision());
            _points.add(p);
            _parameterization.add(oldTrack._parameterization.get(i));
        }

        calculateLength();

        reduce();
    }

    public Track reverse () {
        if (null == _reverse) {
            List<Position> reversePoints = new ArrayList<Position>(_points);
            Collections.reverse(reversePoints);
            List<Double> reverseParameterization = new ArrayList<Double>(_parameterization);
            Collections.reverse(reverseParameterization);
            for (int i=0; i<reverseParameterization.size(); ++i)
                reverseParameterization.set(i, 1.0-reverseParameterization.get(i));
            _reverse = createTrack(reversePoints, reverseParameterization);
            _reverse._reverse = this;
        }
        return _reverse;
    }

    public PositionCalculationParameters getParameters () {
        return _parameters;
    }

    public List<Position> getPoints () {
        return _points;
    }

    public double getLength () {
        return _length;
    }

    protected List<Double> getParameterization () {
        return _parameterization;
    }

    /*
     * This is currently used only for filling in intermediate points for
     * cartesian tracks. As such, the length always stays the same. If it is
     * ever used to actually change the underlying geometry, it would have to
     * take a new length too.
     */
    protected void updatePoints (List<Position> points, List<Double> parameterization) {
        _points = points;
        _parameterization = parameterization;
    }

    /*
     * Universal precalculation step - figure out the length parameterization of the trajectory
     */
    private void calculateLengthParameterization () {
        _parameterization = new ArrayList<Double>(_points.size());

        calculateLength();
        double cumulativeLength = 0.0;

        Position lastPoint = null;
        for (Position point: _points) {
            if (0 == _length) {
                _parameterization.add(0.0);
            } else {
                if (null != lastPoint) {
                    cumulativeLength += getSegmentDistance(lastPoint, point);
                }
                _parameterization.add(cumulativeLength/_length);
                lastPoint = point;
            }
        }
    }

    /** Get the distance from the start to the end position */
    abstract protected double getSegmentDistance (Position start, Position end);
    /** Interpolate 100*t percent of the way from start to end */
    abstract protected Position interpolate (Position start, Position end, double t);
    /**
     * Get the amount of error that would be introduced by removing point b,
     * lieing between points a and c from this track
     */
    abstract protected double getRelativeError (Position a, Position b, Position c);
    
    /** Create a track of the current type from the listed points */
    abstract protected Track createTrack (List<Position> points);
    /** Create a track of the current type from the listed points, whose parameterization is already calculated. */
    abstract protected Track createTrack (List<Position> path, List<Double> parameterization);

    /*
     * Get the relative importance of point b relative to points a and c.
     * Basically, this is a measure of how far b is from each; it will be 1 when
     * b is equidistant from each, and will go up from there. We do this so as
     * to preserve corners.
     * 
     * The relative importance should be 0 if b is a duplicate of a or c.
     */
    private double getRelativeImportance (Position a, Position b, Position c) {
        double dab = getSegmentDistance(a, b);
        double dbc = getSegmentDistance(b, c);

        if (dab < getParameters().getPrecision()) return 0; // duplicate point, doesn't matter at all
        if (dbc < getParameters().getPrecision()) return 0; // duplicate point, doesn't matter at all

        return Math.max(dab/dbc, dbc/dab);
    }



    private boolean closerToReverse (Track them) {
        // First try matching endpoints
        Position ourStart = _points.get(0);
        Position ourEnd = _points.get(_points.size()-1);

        List<Position> theirPoints = them.getPoints();
        Position theirStart = theirPoints.get(0);
        Position theirEnd = theirPoints.get(theirPoints.size()-1);

        double dss = getSegmentDistance(ourStart, theirStart);
        double dse = getSegmentDistance(ourStart, theirEnd);
        double des = getSegmentDistance(ourEnd, theirStart);
        double dee = getSegmentDistance(ourEnd, theirEnd);

        double requiredConfidence = 0.5;
        return (dse/dee < requiredConfidence && des / dss < requiredConfidence);
    }

    public double getDistance (Track them) {
        if (_parameters.ignoreDirection()) {
            if (closerToReverse(them))
                return getDistanceWithDirection(them.reverse());
            else
                return getDistanceWithDirection(them);
        } else {
            return getDistanceWithDirection(them);
        }
    }

    private double getDistanceWithDirection (Track them) {
        List<Double> ourParameterization = getParameterization();
        List<Double> theirParameterization = them.getParameterization();
        List<Double> joinedParameterization = ListUtilities.joinLists(ourParameterization,
                                                                      theirParameterization,
                                                                      _parameters.getPrecision());

        Position pALast = null;
        Position pBLast = null;
        double dLast = 0;
        double totalDistance = 0.0;
        for (double d : joinedParameterization) {
            Position pA = getLengthParamterizedPoint(d);
            Position pB = them.getLengthParamterizedPoint(d);
            if (null != pALast) {
                double startDistance = getSegmentDistance(pALast, pBLast);
                double endDistance = getSegmentDistance(pA, pB);
                totalDistance += (startDistance + endDistance) / 2
                                 * (d - dLast);
            }

            dLast = d;
            pALast = pA;
            pBLast = pB;
        }

        return totalDistance / ((_length + them._length) / 2.0);
    }



    public Track weightedAverage (Track them, double ourWeight, double theirWeight) {
        if (_parameters.ignoreDirection()) {
            if (closerToReverse(them))
                return weightedAverageWithDirection(them.reverse(), ourWeight, theirWeight);
            else
                return weightedAverageWithDirection(them, ourWeight, theirWeight);
        } else {
            return weightedAverageWithDirection(them, ourWeight, theirWeight);
        }
    }
    protected Track weightedAverageWithDirection (Track them, double ourWeight, double theirWeight) {
        double theirRelWeight = theirWeight/(ourWeight+theirWeight);

        // Get all parameterization points
        List<Double> ourParameterization = getParameterization();
        List<Double> theirParameterization = them.getParameterization();
        List<Double> joinedParameterization = ListUtilities.joinLists(ourParameterization,
                                                                      theirParameterization,
                                                                      _parameters.getPrecision());

        // Average the tracks along each parameterization point
        List<Position> meanPath = new ArrayList<Position>();
        for (double d: joinedParameterization) {
            Position pUs = getLengthParamterizedPoint(d);
            Position pThem = them.getLengthParamterizedPoint(d);
            Position weightedMean = interpolate(pUs, pThem, theirRelWeight);
            meanPath.add(weightedMean);
        }

        return createTrack(meanPath);
    }

    protected String getLabel () {
        return "Trajectory";
    }
    @Override
    public String toString () {
        StringBuffer result = new StringBuffer();
        result.append(getLabel());
        result.append("<");
        result.append(getParameters().getCalculationType());
        result.append(">[");
        for (int i=0; i<_points.size(); ++ i) {
            if (0 < i) result.append(", ");
            result.append(_points.get(i));
        }
        result.append("]");
        return result.toString();
    }

    @Override
    public boolean equals (Object that) {
        if (this == that) return true;
        if (null == that) return false;
        if (!(that instanceof Track)) return false;

        Track track = (Track) that;
        return getDistance(track) < getParameters().getPrecision();
    }



    /*
     * Remove points that don't add anything significant to this path.
     * Insignificance is defined by _parameter.getAllowedError(), and the
     * definition of {@link #getRelativeError}.
     */
    private void reduce () {
        double minImportance = Math.sqrt(1.0/_parameters.getAllowedError());
        // Remove points that don't contribute much
        int i=0;
        while (i<_points.size()-2) {
            Position a = _points.get(i);
            Position b = _points.get(i+1);
            Position c = _points.get(i+2);
            boolean remove = (a.equals(b) || c.equals(b));
            if (!remove) {
                double relativeError = getRelativeError(a, b, c);
                double relativeImportance = getRelativeImportance(a, b, c);
                remove = (relativeError < _parameters.getAllowedError());
                if (remove && minImportance < relativeImportance)
                    remove = false;
                if (remove) {
                    getRelativeError(a, b, c);
                    getRelativeImportance(a, b, c);
                }
            }
            if (remove) {
                _points.remove(i+1);
                _parameterization.remove(i+1);
            } else {
                ++i;
            }
        }

        // finally, check for duplication in the last two points
        int n = _points.size();
        if (n>1) {
            Position a = _points.get(n-2);
            Position b = _points.get(n-1);
            if (a.equals(b)) {
                _points.remove(n-2);
                _parameterization.remove(n-2);
            }
        }
    }



    // Should only be called for initialization of _length; after initial
    // setting, that should be trusted
    private void calculateLength () {
        _length = 0.0;
        if (null == _points || _points.size() < 2)
            return;

        Position lastPoint = null;
        for (Position point : _points) {
            if (null != lastPoint)
                _length += getSegmentDistance(point, lastPoint);
            lastPoint = point;
        }
    }

    protected Position getLengthParamterizedPoint (double parameter) {
        if (parameter < 0.0 || 1.0 < parameter)
            throw new IllegalArgumentException("Length paramterization parameter must be between 0 and 1");

        // Only one point; return it.
        int N = _parameterization.size();
        if (0 == N) return null;
        if (1 == N) return _points.get(0);

        int n;
        for (n = 0; n < N - 1; ++n) {
            if (_parameterization.get(n + 1) > parameter)
                break;
        }

        double startD = _parameterization.get(n);
        Position start = _points.get(n);
        if (Math.abs(parameter - startD) < _parameters.getPrecision())
            return start;

        double endD = _parameterization.get(n + 1);
        Position end = _points.get(n+1);
        if (Math.abs(parameter-endD) < _parameters.getPrecision())
            return end;

        double pSeg = (parameter - startD) / (endD - startD);
        return interpolate(start, end, pSeg);
    }


    public void addStatistic (String statName, double stat) {
        if (null == _statistics)
            _statistics = new HashMap<String, Double>();

        _statistics.put(statName, stat);
    }

    public Map<String, Double> getStatistics () {
        return _statistics;
    }
}
 public class GeodeticTrack extends Track {
    public GeodeticTrack (PositionCalculationParameters parameters,
                          Position... points) {
        super(parameters, points);
    }

    public GeodeticTrack (PositionCalculationParameters parameters,
                          List<Position> points) {
        super(parameters, points);
    }

    public GeodeticTrack (PositionCalculationParameters parameters,
                          List<Position> points,
                          List<Double> parameterization) {
        super(parameters, points, parameterization);
    }

    public GeodeticTrack (Track oldTrack) {
        super(oldTrack,
              new PositionCalculationParameters(PositionCalculationType.Geodetic,
                                                oldTrack.getParameters().getAllowedError(),
                                                oldTrack.getParameters().getPrecision(),
                                                oldTrack.getParameters().ignoreDirection()));
    }

    public GeodeticTrack (Track oldTrack, double allowedError, double precision, boolean ignoreDirection) {
        super(oldTrack,
              new PositionCalculationParameters(PositionCalculationType.Geodetic, allowedError, precision, ignoreDirection));
    }

    @Override
    protected double getSegmentDistance (Position start, Position end) {
        return start.getAngularDistance(end);
    }

    @Override
    protected Position interpolate (Position start, Position end, double t) {
        if (null == start || null == end)
            return null;

        double lenSeg = start.getAngularDistance(end);
        double azSeg = start.getAzimuth(end);
        return start.offset(azSeg, lenSeg * t);
//        Vector vs = start.getAsCartesian();
//        Vector ve = end.getAsCartesian();
//        Vector vavg = vs.add(ve.subtract(vs).scale(t));
//        Position p = new Position(vavg.coord(0), vavg.coord(1), vavg.coord(2), true);
//        if (start.hasElevation()) {
//            double es = start.getElevation();
//            double ee = end.getElevation();
//            double e = es + (ee-es)*t;
//            p = new Position(p.getLongitude(), p.getLatitude(), e);
//        } else {
//            p = new Position(p.getLongitude(), p.getLatitude());
//        }
//        return p;
    }

    @Override
    protected double getRelativeError (Position a, Position b, Position c) {
        // Get the area of the triangle subtended by the three, divided by
        // the square of the length of AC
        double dab = a.getAzimuth(b);
        double dac = AngleUtilities.intoRangeDegrees(dab, a.getAzimuth(c));
        double da = Math.abs(dab-dac);
        double dbc = b.getAzimuth(c);
        double dba = AngleUtilities.intoRangeDegrees(dbc, b.getAzimuth(a));
        double db = Math.abs(dbc-dba);
        double dca = c.getAzimuth(a);
        double dcb = AngleUtilities.intoRangeDegrees(dca, c.getAzimuth(b));
        double dc = Math.abs(dca-dcb);

        double triangleArea = Math.toRadians(da+db+dc-180); // area of triangle/r62
        double longSide = Math.toRadians(a.getAngularDistance(c)); // length of long side
        return triangleArea/(longSide*longSide);
    }

    @Override
    protected Track createTrack (List<Position> points) {
        return new GeodeticTrack(getParameters(), points);
    }

    @Override
    protected Track createTrack (List<Position> points,
                                 List<Double> parameterization) {
        return new GeodeticTrack(getParameters(), points,
                                 parameterization);
    }
}
 public class DataSet implements Serializable, Iterable<Instance> {
	private static final long serialVersionUID = -544003860939601075L;
	
	private final Map<String, Instance> map = new LinkedHashMap<String, Instance>();

	private class DataSetIterator implements Iterator<Instance> {
		Iterator<Entry<String, Instance>> iterator;
		
		public DataSetIterator(Map<String, Instance> map) {
			this.iterator = map.entrySet().iterator();
		}
		
		@Override
		public boolean hasNext() {
			return iterator.hasNext();
		}
		@Override
		public Instance next() {
			return iterator.next().getValue();
		}
		@Override
		public void remove() {
			iterator.remove();
		}
	};
	
	public Set<String> getKeys() {
		return map.keySet();
	}
	
	/***
	 * Add an Instance to the Dataset
	 * @param inst the Instance to add
	 * @return whether an instance was replace with the same id
	 */
	public boolean add(Instance inst) {
		return (map.put(inst.getId(), inst) != null);
	}
	
	/***
	 * Remove an instance from the Dataset
	 * @param inst the Instance to remove
	 * @return the removed Instance
	 */
	public Instance remove(Instance inst) {
		return map.remove(inst.getId());
	}

	/***
	 * Get an instance by Id
	 * @param id the id of the Instance to return
	 * @return the Instance with id or null if no Instance exists
	 */
	public Instance get(String id) {
		return map.get(id);
	}
	
	/***
	 * Return the number of Instances in the DataSet
	 * @return
	 */
	public int size() {
		return map.size();
	}
	
	/***
	 * Return a DataSet with random fraction number of Instances as this DataSet
	 * @param fraction fraction of Instances to return in the range of 0 and 1
	 * @return resulting sample DataSet
	 */
	public DataSet sample(double fraction) {
		DataSet sample = new DataSet();
		
		// max fraction is 1
		if (fraction > 1) fraction = 1;
		
		// no samples are requested - return an empty data set
		if (fraction <= 0) return sample;
		
		long numToSample = Math.round( map.size() * fraction );
		
		ArrayList<String> keys = new ArrayList<String>(map.keySet());
		
		// randomly pick k instances as the initial k means
		ArrayList<String> indexes = new ArrayList<String>(keys.size()); 
		for (int i = 0; i < indexes.size(); i++) {
			indexes.add( keys.get(i) );
		}
		Collections.shuffle(indexes);  // permutate the indexes
		for (int i = 0; i < numToSample; i++) {
			sample.add( get(indexes.get(i)) );
		}
		return sample;
	}
	
	private void swap(Instance a, Instance b) {
		Instance tmp = a;
		a = b;
		b = tmp;
	}
	
	private void shuffle(Instance[] array) {
		Random rnd = new Random();
		int n = array.length;

		while (n > 1) {
			int k = rnd.nextInt(n);
			n--;
			swap(array[n], array[k]);
		}
	}
	
	/***
	 * Randomly split this DataSet into n similarly sized DataSets.
	 * @param n the number of folds to split the DataSet - value must be greater than 1 and not greater than the number of instances in DataSet
	 * @return a list of n DataSets
	 */
	public List<DataSet> nFolds(int n) {
		// Make sure n is valid: each fold must have at least one instance!
		if (n > size() || n < 1) return null;
		
		List<DataSet> folds = new LinkedList<DataSet>();
		
		Instance[] instances = new Instance[size()];
		instances = map.values().toArray(instances);
		shuffle(instances);
	
		int sliceSize = size() / n;
		int extra = size() % n;
		
		// create n folds
		for (int i=0; i < n; i++) {
			DataSet fold = new DataSet();
			folds.add(fold);
			
			int start = i*sliceSize;
			int end = start + sliceSize;
			for (int j=start; j < end; j++) {
				fold.add(instances[j]);
			}
		}
		
		// evenly distribute any extra instances
		for (int i=0; i < extra; i++) {
			int offset = instances.length - extra + i;
			folds.get(i).add(instances[offset]);
		}
		
		return folds;
	}
	
	/***
	 * Normalize the specified Feature for all Instances in this DataSet.
	 * 
	 * Currently only NumericVectorFeature types are supported.
	 * 
	 * @param featureName the name of the feature to normalize
	 */
	public void normalizeInstanceFeature(String featureName) {
		List<Feature> allFeatures = new ArrayList<Feature>();
		
		// gather up all matching features
		for (Instance inst : this) {
			if (inst.containsFeature(featureName)) {
				allFeatures.add( inst.getFeature(featureName) );
			}
		}
		
		if (allFeatures.isEmpty()) return;
		
		double N = allFeatures.size();
		
		// currently only support normalizing numeric vector features
		if ( (allFeatures.get(0) instanceof NumericVectorFeature) == false ) return;
		
		double[] meanVector = ((NumericVectorFeature)allFeatures.get(0)).getValue().clone();
		
		// compute mean of feature
		for (int i=1; i < allFeatures.size(); i++) {
			NumericVectorFeature v = (NumericVectorFeature)allFeatures.get(i);
			double[] vals = v.getValue();
			
			for (int j=0; j < meanVector.length; j++) {
				meanVector[j] += vals[j];
			}
		}
		for (int i=0; i < meanVector.length; i++) {
			meanVector[i] /= N;
		}
		
		double[] stdevVector = new double[meanVector.length];
		
		// compute stdev of feature
		for (int i=0; i < allFeatures.size(); i++) {
			NumericVectorFeature v = (NumericVectorFeature)allFeatures.get(i);
			double[] vals = v.getValue();
			
			for (int j=0; j < meanVector.length; j++) {
				stdevVector[j] += (vals[j] - meanVector[j])*(vals[j] - meanVector[j]);
			}
		}
		for (int i=0; i < stdevVector.length; i++) {
			stdevVector[i] = Math.sqrt( stdevVector[i] / (N-1) );
		}
		
		// normalize each feature vector
		for (int i=0; i < allFeatures.size(); i++) {
			NumericVectorFeature v = (NumericVectorFeature)allFeatures.get(i);
			double[] vals = v.getValue();
			
			for (int j=0; j < vals.length; j++) {
				vals[j] = (vals[j] - meanVector[j]) / stdevVector[j];
			}
		}
	}
	 
	@Override
	public Iterator<Instance> iterator() {
		return new DataSetIterator(map);
	}

	/***
	 * Add a collection of Instances to the DataSet
	 * @param c the Instances to add
	 * @return true if an existing Instance with a matching id in the DataSet was replaced 
	 */
	public boolean addAll(Collection<Instance> c) {
		boolean altered = false;
		
		for (Instance i : c) {
			if ( add(i) ) altered = true;
		}
		return altered;
	}

	/***
	 * Remove all Instances from this DataSet
	 */
	public void clear() {
		map.clear();
	}

	/***
	 * Return true if the DataSet contains the specified Instance
	 * @param inst the Instance to test
	 * @return true if the Instance is a member of the DataSet
	 */
	public boolean contains(Instance inst) {
		return map.containsKey(inst);
	}

	/***
	 * Return true if the DataSet contains all the specified Instances
	 * @param c the Instances to test
	 * @return true if the Instances are all members of the DataSet
	 */
	public boolean containsAll(Collection<Instance> c) {
		for (Instance i : c) {
			if (map.containsKey(i.getId()) == false) return false;
		}
		return true;
	}

	/***
	 * Returns true if this DataSet contains no Instances
	 * @return true if this DataSet contains no Instances
	 */
	public boolean isEmpty() {
		return map.isEmpty();
	}
}
 public class Instance implements Serializable {
	private static final long serialVersionUID = -8781788906032267606L;
	
	protected String 	id;
	protected String	classLabel;
	protected Map<String, Feature> features = new LinkedHashMap<String, Feature>();
	
	public Instance() { 
		this(UUID.randomUUID().toString());
	}
	
	/***
	 * Constructor to specify your own id's
	 * 
	 * NOTE: id's MUST be unique!  It's is up to the caller to ensure this.
	 * 
	 * @param id
	 */
	public Instance(String id) {
		this.id = id;
	}
	
	/***
	 * Return whether this Instance has the specified label
	 * 
	 * @param label label to test
	 * @return true if the Instance has a matching class label
	 */
	public boolean hasClassLabel(String label) {
		return classLabel.equalsIgnoreCase(label);
	}
	
	/***
	 * Set the class label for the Instance
	 * @param label label to apply to the Instance
	 */
	public void setClassLabel(String label) {
		classLabel = label;
	}
	
	/***
	 * Return the class label given to this Instance or null if none 
	 * @return class label or null if none
	 */
	public String getClassLabel() {
		return classLabel;
	}
	
	/***
	 * The unique id for this Instance
	 * @return the Instance id
	 */
	public String getId() {
		return id;
	}
	
	/***
	 * Set the unique Instance id
	 * @param id id of Instance
	 */
	public void setId(String id) {
		this.id = id;
	}
	
	@Override
	public int hashCode() {
		return id.hashCode();
	}
	
	@Override
	public String toString() {
		StringBuilder str = new StringBuilder();
		str.append("\"id:" + id + "\",");
		int i = 1;
		for (Feature feature : features.values()) {
			str.append("\"" + feature.toString() + "\"");
			if (i < features.size()) str.append(",");
			i++;
		}
		return str.toString();
	}
	
	/***
	 * Add a Feature to this Instance.  Features are typed data attributes that describe the Instance.
	 * 
	 * Each Feature must have a unique id to distinguish them.
	 * 
	 * @param feature feature to add to the Instance
	 */
	public void addFeature(Feature feature) {
		features.put(feature.getId(), feature);
	}

	/***
	 * Add a collection of Features to this Instance.  Features are typed data attributes that describe the Instance.
	 * 
	 * Each Feature must have a unique id to distinguish them.
	 * 
	 * @param features features to add to the Instance
	 */
	public void addFeatures(Collection<? extends Feature> features) {
		for (Feature f : features) {
			addFeature(f);
		}
	}

	/***
	 * Return true if this Instance has a Feature with the specified name
	 * @param featureName featureName to check
	 * @return true if this Instance contain a Feature with name
	 */
	public boolean containsFeature(String featureName) {
		return features.containsKey(featureName);
	}
	
	/***
	 * Return the Feature with the specified name or null if the Instance doen't have a matching Feature
	 * @param featureName name of the Feature
	 * @return the Feature with the specified name or null if the Instance don't have a matching Feature
	 */
	public Feature getFeature(String featureName) {
		return features.get(featureName);
	}

	/***
	 * Return true if the Instance contains no Features
	 * @return true if the Instance contains no Features
	 */
	@JsonIgnore
	public boolean isEmpty() {
		return features.isEmpty();
	}

	/***
	 * Return all the Features associated with this Instance
	 * @return collection of Features
	 */
	@JsonIgnore
	public Collection<Feature> getAllFeatures() {
		return features.values();
	}

	/***
	 * Return a map of the Features associated with this Instance keyed by Feature name
	 * @return a Map of Features
	 */
	public Map<String, Feature> getFeatures() {
		return features;
	}

	/***
	 * Add all Features in a Map to the Instance whether the Map is keyed by Feature name
	 * @param feature Map to add
	 */
	public void setFeatures(Map<String, Feature> features) {
		features.putAll(features);
	}

	/***
	 * Remove the Feature with the specified name from the Instance
	 * @param name of the Feature to remove
	 */
	public void removeFeature(String featureName) {	
		features.remove(featureName);
	}

	/***
	 * Remove all the Features associated with this Instance
	 */
	public void removeAllFeatures() {
		features.clear();
	}

	/***
	 * Return the number of Features associated with this Instance
	 * @return the number of Features
	 */
	public int numFeatures() {
		return features.size();
	}	
}
 public class TrackFeature extends Feature {
    private static final long serialVersionUID = -7923265691042687996L;



    private Track track;
    
    @Override
    public String toString() {
        return (this.getName() + " : " + track.toString());
    }
    
    public TrackFeature() {
        super();
    }
    
    public TrackFeature(String name) {
        super(name);
    }

    public void setValue(Track track) {
        this.track = track;
    }
    
    public Track getValue() {
        return track;
    }
}
 public class TrackCentroid implements Centroid<TrackFeature>{
    private static final long serialVersionUID = 1L;

    private String name;
    private double _weight;
    private Track  _centroid;
    private int    _changes;

    public TrackCentroid () {
        _weight = 0;
        _centroid = null;
    }

    public int getChanges () {
        return _changes;
    }

    @Override
    public void add(TrackFeature feature) {
        ++_changes;
        Track newTrack = feature.getValue();
        // add the track to the list of members
        if (null == _centroid) {
            _centroid = newTrack;
            _weight = feature.getWeight();
        } else {
            // Weight according to how many are already there
            _centroid = _centroid.weightedAverage(newTrack, _weight, feature.getWeight());
            _weight = _weight + feature.getWeight();
        }
    }

    @Override
    public void remove (TrackFeature feature) {
        ++_changes;
        Track toRemove = feature.getValue();

        //remove the track from the list of members
        double newWeight = _weight - feature.getWeight();
        if (newWeight <= 0) {
            _centroid = null;
            _weight = 0.0;
        } else {
            _centroid = _centroid.weightedAverage(toRemove, _weight, -feature.getWeight());
            _weight = newWeight;
        }
    }

    @Override
    public void setName(String name) {
        this.name = name;
    }

    @Override
    public String getName() {
        return this.name;
    }

    // Erasure makes this definitionally ok; type checking at compile time
    // should still happen, and nothing would be different at run time if we got
    // something that the compiler though really was of the right type, properly
    // generified.
    @SuppressWarnings({ "unchecked", "rawtypes" })
    @Override
    public Class<TrackFeature> getType() {
        return (Class) TrackFeature.class;
    }

    @Override
    public Collection<TrackFeature> getAggregatableCentroid () {
        return Collections.singleton(getCentroid());
    }

    @Override
    public TrackFeature getCentroid() {
        TrackFeature centroid = new TrackFeature(name);
        centroid.setValue(_centroid);
        centroid.setWeight(_weight);
        return centroid;
    }

	@Override
	public void reset() {
		_weight = 0;
		_centroid = null;
		_changes = 0;
	}
}
 public class TrackDistance extends DistanceFunction<TrackFeature>{
    private static final long serialVersionUID = -6853731905216910630L;



    public TrackDistance(double weight) {
        super(weight);
    }
    
    @Override
    public double distance(TrackFeature x, TrackFeature y) {
        return x.getValue().getDistance(y.getValue());
    }
}
 public interface ClusterResult extends Iterable<Cluster>, Serializable {

	public boolean isEmpty();
	
	public int size();
} public class BCubed {
	private double precision = 0, recall = 0, f = 0;
	
	/***
	 * BCubed recall is the proportion of instances with class label are in this cluster.
	 * @param label is the class label we are evaluating
	 * @param cluster that we are evaluating
	 * @param clusters is a collection of all the clusters
	 * @return the proportion of instances with class label that are members of cluster
	 */
	public double recall(String label, Cluster cluster, Collection<Cluster> clusters) {
		double count = 0, total = 0, recall = 0;
		
		// first determine how many instances have class label in this cluster
		for (Instance i : cluster.getMembers()) {
			if (i.hasClassLabel(label)) count++;
		}
		
		// find out how many instances in total have this class label
		for (Cluster c : clusters) {
			for (Instance i : c.getMembers()) {
				if (i.hasClassLabel(label)) total++;
			}
		}
		
		if (total > 0) recall = count / total;
		
		return recall;
	}
	
	/***
	 * BCubed precision is the proportion of instances in a cluster that share the same label.
	 * 
	 * @param label is the class label we are evaluating
	 * @param cluster that we are evaluating
	 * @return the proportion of instances in cluster that have the specified label
	 */
	public double precision(String label, Cluster cluster) {
		double count = 0, precision = 0;
		
		for (Instance i : cluster.getMembers()) {
			if (i.hasClassLabel(label)) count++; 
		}
		
		if (count > 0) precision = count / cluster.getMembers().size(); 
		
		return precision;
	}
	
	public double validate(Collection<Cluster> clusters) {
		int num = 0;
		double p = 0, r = 0;
		
		for (Cluster c : clusters) {
			for (Instance inst : c.getMembers()) {
				p += precision(inst.getClassLabel(), c) ;
				r += recall(inst.getClassLabel(), c, clusters);
				num++;
			}
		}
		
		// calculate average precision and recall
		precision = p / num;
		recall = r / num;
		
		// calculate and return the f score
		f = 2 * precision * recall / (precision + recall);
		
		return f;
	}
	
	public double getPrecision() {
		return precision;
	}
	
	public double getRecall() {
		return recall;
	}
	
	public double getBCubed() {
		return f;
	}
}
import java.awt.BasicStroke;
import java.awt.Color;
import java.awt.Dimension;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.Point;
import java.awt.geom.Rectangle2D;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
import java.util.logging.Logger;
import junit.framework.Assert;
import org.junit.Ignore;
import org.junit.Test;
import com.oculusinfo.geometry.cartesian.CubicBSpline;
import com.oculusinfo.geometry.geodesic.Position;
import com.oculusinfo.geometry.geodesic.PositionCalculationParameters;
import com.oculusinfo.geometry.geodesic.PositionCalculationType;
import com.oculusinfo.geometry.geodesic.Track;
import com.oculusinfo.geometry.geodesic.tracks.Cartesian3DTrack;
import com.oculusinfo.geometry.geodesic.tracks.GeodeticTrack;
import com.oculusinfo.math.linearalgebra.Vector;
import com.oculusinfo.ml.DataSet;
import com.oculusinfo.ml.Instance;
import com.oculusinfo.ml.feature.spatial.TrackFeature;
import com.oculusinfo.ml.feature.spatial.centroid.TrackCentroid;
import com.oculusinfo.ml.feature.spatial.distance.TrackDistance;
import com.oculusinfo.ml.unsupervised.cluster.Cluster;
import com.oculusinfo.ml.unsupervised.cluster.ClusterResult;
import com.oculusinfo.ml.unsupervised.cluster.kmeans.KMeans;
import com.oculusinfo.ml.validation.unsupervised.external.BCubed;

/**
 * Copyright (c) 2013 Oculus Info Inc.
 * http://www.oculusinfo.com/
 *
 * Released under the MIT License.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy of
 * this software and associated documentation files (the "Software"), to deal in
 * the Software without restriction, including without limitation the rights to
 * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
 * of the Software, and to permit persons to whom the Software is furnished to do
 * so, subject to the following conditions:

 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.

 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */
package com.oculusinfo.ml.tracks;





public class TrackClusteringTests {
    private static final Logger LOGGER = Logger.getLogger(TrackClusteringTests.class.getName());

    private static final PositionCalculationParameters GEODETIC_PARAMETERS =
