    001package net.sf.cpsolver.ifs.util;
    002
    003import java.util.HashMap;
    004import java.util.Map;
    005
    006/**
    007 * Common class for computing distances and back-to-back instructor / student conflicts.
    008 *
    009 * When property Distances.Ellipsoid is set, the distances are computed using the given (e.g., WGS84, see {@link Ellipsoid}).
    010 * In the legacy mode (when ellipsoid is not set), distances are computed using Euclidian distance and 1 unit is considered 10 meters.
    011 * <br><br>
    012 * For student back-to-back conflicts, Distances.Speed (in meters per minute) is considered and compared with the break time
    013 * of the earlier class.
    014 * <br><br>
    015 * For instructors, the preference is computed using the distance in meters and the three constants
    016 * Instructor.NoPreferenceLimit (distance <= limit -> no preference), Instructor.DiscouragedLimit (distance <= limit -> discouraged),
    017 * Instructor.ProhibitedLimit (distance <= limit -> strongly discouraged), the back-to-back placement is prohibited when the distance is over the last limit.
    018 *
    019 * @version IFS 1.2 (Iterative Forward Search)<br>
    020 *          Copyright (C) 2006 - 2010 Tomáš Müller<br>
    021 *          <a href="mailto:muller@unitime.org">muller@unitime.org</a><br>
    022 *          <a href="http://muller.unitime.org">http://muller.unitime.org</a><br>
    023 * <br>
    024 *          This library is free software; you can redistribute it and/or modify
    025 *          it under the terms of the GNU Lesser General Public License as
    026 *          published by the Free Software Foundation; either version 3 of the
    027 *          License, or (at your option) any later version. <br>
    028 * <br>
    029 *          This library is distributed in the hope that it will be useful, but
    030 *          WITHOUT ANY WARRANTY; without even the implied warranty of
    031 *          MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
    032 *          Lesser General Public License for more details. <br>
    033 * <br>
    034 *          You should have received a copy of the GNU Lesser General Public
    035 *          License along with this library; if not see
    036 *          <a href='http://www.gnu.org/licenses/'>http://www.gnu.org/licenses/</a>.
    037 */
    038public class DistanceMetric {
    039    public static enum Ellipsoid {
    040        LEGACY ("Euclidean metric (1 unit equals to 10 meters)", "X-Coordinate", "Y-Coordinate", 0, 0, 0),
    041        WGS84 ("WGS-84 (GPS)", 6378137, 6356752.3142, 1.0 / 298.257223563),
    042        GRS80 ("GRS-80", 6378137, 6356752.3141, 1.0 / 298.257222101),
    043        Airy1830 ("Airy (1830)", 6377563.396, 6356256.909, 1.0 / 299.3249646),
    044        Intl1924 ("Int'l 1924", 6378388, 6356911.946, 1.0 / 297),
    045        Clarke1880 ("Clarke (1880)", 6378249.145, 6356514.86955, 1.0 / 293.465),
    046        GRS67 ("GRS-67", 6378160, 6356774.719, 1.0 / 298.25);
    047
    048        private double iA, iB, iF;
    049        private String iName, iFirstCoord, iSecondCoord;
    050
    051        Ellipsoid(String name, double a, double b) {
    052            this(name, "Latitude", "Longitude", a, b, (a - b) / a);
    053        }
    054        Ellipsoid(String name, double a, double b, double f) {
    055            this(name, "Latitude", "Longitude", a, b, f);
    056        }
    057        Ellipsoid(String name, String xCoord, String yCoord, double a, double b, double f) {
    058            iName = name;
    059            iFirstCoord = xCoord; iSecondCoord = yCoord;
    060            iA = a; iB = b; iF = f;
    061        }
    062
    063        /** Major semiaxe A */
    064        public double a() { return iA; }
    065        /** Minor semiaxe B */
    066        public double b() { return iB; }
    067        /** Flattening (A-B) / A */
    068        public double f() { return iF; }
    069
    070        /** Name of this coordinate system */
    071        public String getEclipsoindName() { return iName; }
    072        /** Name of the fist coordinate (e.g., Latitude) */
    073        public String getFirstCoordinateName() { return iFirstCoord; }
    074        /** Name of the second coordinate (e.g., Longitude) */
    075        public String getSecondCoordinateName() { return iSecondCoord; }
    076    }
    077
    078    /** Elliposid parameters, default to WGS-84 */
    079    private Ellipsoid iModel = Ellipsoid.WGS84;
    080    /** Student speed in meters per minute (defaults to 1000 meters in 15 minutes) */
    081    private double iSpeed = 1000.0 / 15;
    082    /** Back-to-back classes: maximal distance for no preference */
    083    private double iInstructorNoPreferenceLimit = 0.0;
    084    /** Back-to-back classes: maximal distance for discouraged preference */
    085    private double iInstructorDiscouragedLimit = 50.0;
    086    /**
    087     * Back-to-back classes: maximal distance for strongly discouraged preference
    088     * (everything above is prohibited)
    089     */
    090    private double iInstructorProhibitedLimit = 200.0;
    091    /**
    092     * When Distances.ComputeDistanceConflictsBetweenNonBTBClasses is enabled, distance limit (in minutes)
    093     * for a long travel.
    094     */
    095    private double iInstructorLongTravelInMinutes = 30.0;
    096
    097    /** Default distance when given coordinates are null. */
    098    private double iNullDistance = 10000.0;
    099    /** Maximal travel time in minutes when no coordinates are given. */
    100    private int iMaxTravelTime = 60;
    101    /** Travel times overriding the distances computed from coordintaes */
    102    private Map<Long, Map<Long, Integer>> iTravelTimes = new HashMap<Long, Map<Long,Integer>>();
    103    /** Distance cache  */
    104    private Map<String, Double> iDistanceCache = new HashMap<String, Double>();
    105    /** True if distances should be considered between classes that are NOT back-to-back */
    106    private boolean iComputeDistanceConflictsBetweenNonBTBClasses = false;
    107
    108    /** Default properties */
    109    public DistanceMetric() {
    110    }
    111
    112    /** With provided ellipsoid */
    113    public DistanceMetric(Ellipsoid model) {
    114        iModel = model;
    115        if (iModel == Ellipsoid.LEGACY) {
    116            iSpeed = 100.0 / 15;
    117            iInstructorDiscouragedLimit = 5.0;
    118            iInstructorProhibitedLimit = 20.0;
    119        }
    120    }
    121
    122    /** With provided ellipsoid and student speed */
    123    public DistanceMetric(Ellipsoid model, double speed) {
    124        iModel = model;
    125        iSpeed = speed;
    126    }
    127
    128    /** Configured using properties */
    129    public DistanceMetric(DataProperties properties) {
    130        if (Ellipsoid.LEGACY.name().equals(properties.getProperty("Distances.Ellipsoid",Ellipsoid.LEGACY.name()))) {
    131            //LEGACY MODE
    132            iModel = Ellipsoid.LEGACY;
    133            iSpeed = properties.getPropertyDouble("Student.DistanceLimit", 1000.0 / 15) / 10.0;
    134            iInstructorNoPreferenceLimit = properties.getPropertyDouble("Instructor.NoPreferenceLimit", 0.0);
    135            iInstructorDiscouragedLimit = properties.getPropertyDouble("Instructor.DiscouragedLimit", 5.0);
    136            iInstructorProhibitedLimit = properties.getPropertyDouble("Instructor.ProhibitedLimit", 20.0);
    137            iNullDistance = properties.getPropertyDouble("Distances.NullDistance", 1000.0);
    138            iMaxTravelTime = properties.getPropertyInt("Distances.MaxTravelDistanceInMinutes", 60);
    139        } else {
    140            iModel = Ellipsoid.valueOf(properties.getProperty("Distances.Ellipsoid", Ellipsoid.WGS84.name()));
    141            if (iModel == null) iModel = Ellipsoid.WGS84;
    142            iSpeed = properties.getPropertyDouble("Distances.Speed", properties.getPropertyDouble("Student.DistanceLimit", 1000.0 / 15));
    143            iInstructorNoPreferenceLimit = properties.getPropertyDouble("Instructor.NoPreferenceLimit", iInstructorNoPreferenceLimit);
    144            iInstructorDiscouragedLimit = properties.getPropertyDouble("Instructor.DiscouragedLimit", iInstructorDiscouragedLimit);
    145            iInstructorProhibitedLimit = properties.getPropertyDouble("Instructor.ProhibitedLimit", iInstructorProhibitedLimit);
    146            iNullDistance = properties.getPropertyDouble("Distances.NullDistance", iNullDistance);
    147            iMaxTravelTime = properties.getPropertyInt("Distances.MaxTravelDistanceInMinutes", 60);
    148        }
    149        iComputeDistanceConflictsBetweenNonBTBClasses = properties.getPropertyBoolean(
    150                "Distances.ComputeDistanceConflictsBetweenNonBTBClasses", iComputeDistanceConflictsBetweenNonBTBClasses);
    151        iInstructorLongTravelInMinutes = properties.getPropertyDouble("Instructor.InstructorLongTravelInMinutes", 30.0);
    152    }
    153
    154    /** Degrees to radians */
    155    protected double deg2rad(double deg) {
    156        return deg * Math.PI / 180;
    157    }
    158
    159    /** Compute distance between the two given coordinates
    160     * @deprecated Use @{link {@link DistanceMetric#getDistanceInMeters(Long, Double, Double, Long, Double, Double)} instead (to include travel time matrix when available).
    161     */
    162    @Deprecated
    163    public double getDistanceInMeters(Double lat1, Double lon1, Double lat2, Double lon2) {
    164        if (lat1 == null || lat2 == null || lon1 == null || lon2 == null)
    165            return iNullDistance;
    166
    167        if (lat1.equals(lat2) && lon1.equals(lon2)) return 0.0;
    168
    169        // legacy mode -- euclidian distance, 1 unit is 10 meters
    170        if (iModel == Ellipsoid.LEGACY) {
    171            if (lat1 < 0 || lat2 < 0 || lon1 < 0 || lon2 < 0) return iNullDistance;
    172            double dx = lat1 - lat2;
    173            double dy = lon1 - lon2;
    174            return Math.sqrt(dx * dx + dy * dy);
    175        }
    176
    177        String id = null;
    178        if (lat1 < lat2 || (lat1 == lat2 && lon1 <= lon2)) {
    179            id =
    180                Long.toHexString(Double.doubleToRawLongBits(lat1)) +
    181                Long.toHexString(Double.doubleToRawLongBits(lon1)) +
    182                Long.toHexString(Double.doubleToRawLongBits(lat2)) +
    183                Long.toHexString(Double.doubleToRawLongBits(lon2));
    184        } else {
    185            id =
    186                Long.toHexString(Double.doubleToRawLongBits(lat1)) +
    187                Long.toHexString(Double.doubleToRawLongBits(lon1)) +
    188                Long.toHexString(Double.doubleToRawLongBits(lat2)) +
    189                Long.toHexString(Double.doubleToRawLongBits(lon2));
    190        }
    191        synchronized (iDistanceCache) {
    192            Double distance = iDistanceCache.get(id);
    193
    194            if (distance == null) {
    195                double a = iModel.a(), b = iModel.b(),  f = iModel.f();  // ellipsoid params
    196                double L = deg2rad(lon2-lon1);
    197                double U1 = Math.atan((1-f) * Math.tan(deg2rad(lat1)));
    198                double U2 = Math.atan((1-f) * Math.tan(deg2rad(lat2)));
    199                double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
    200                double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
    201
    202                double lambda = L, lambdaP, iterLimit = 100;
    203                double cosSqAlpha, cos2SigmaM, sinSigma, cosSigma, sigma, sinLambda, cosLambda;
    204                do {
    205                  sinLambda = Math.sin(lambda);
    206                  cosLambda = Math.cos(lambda);
    207                  sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
    208                    (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
    209                  if (sinSigma==0) return 0;  // co-incident points
    210                  cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
    211                  sigma = Math.atan2(sinSigma, cosSigma);
    212                  double sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    213                  cosSqAlpha = 1 - sinAlpha*sinAlpha;
    214                  cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    215                  if (Double.isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (???6)
    216                  double C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
    217                  lambdaP = lambda;
    218                  lambda = L + (1-C) * f * sinAlpha *
    219                    (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
    220                } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);
    221                if (iterLimit==0) return Double.NaN; // formula failed to converge
    222
    223                double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
    224                double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
    225                double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
    226                double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
    227                  B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
    228
    229                // initial & final bearings
    230                // double fwdAz = Math.atan2(cosU2*sinLambda, cosU1*sinU2-sinU1*cosU2*cosLambda);
    231                // double revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
    232
    233                // s = s.toFixed(3); // round to 1mm precision
    234
    235                distance = b*A*(sigma-deltaSigma);
    236                iDistanceCache.put(id, distance);
    237            }
    238
    239            return distance;
    240        }
    241    }
    242
    243    /**
    244     * Compute distance in minutes.
    245     * Property Distances.Speed (in meters per minute) is used to convert meters to minutes, defaults to 1000 meters per 15 minutes (that means 66.67 meters per minute).
    246     * @deprecated Use @{link {@link DistanceMetric#getDistanceInMinutes(Long, Double, Double, Long, Double, Double)} instead (to include travel time matrix when available).
    247     */
    248    @Deprecated
    249    public int getDistanceInMinutes(double lat1, double lon1, double lat2, double lon2) {
    250        return (int) Math.round(getDistanceInMeters(lat1, lon1, lat2, lon2) / iSpeed);
    251    }
    252
    253    /**
    254     * Converts minutes to meters.
    255     * Property Distances.Speed (in meters per minute) is used, defaults to 1000 meters per 15 minutes.
    256     */
    257    public double minutes2meters(int min) {
    258        return iSpeed * min;
    259    }
    260
    261
    262    /** Back-to-back classes in rooms within this limit have neutral preference */
    263    public double getInstructorNoPreferenceLimit() {
    264        return iInstructorNoPreferenceLimit;
    265    }
    266
    267    /** Back-to-back classes in rooms within this limit have discouraged preference */
    268    public double getInstructorDiscouragedLimit() {
    269        return iInstructorDiscouragedLimit;
    270    }
    271
    272    /** Back-to-back classes in rooms within this limit have strongly discouraged preference, it is prohibited to exceed this limit. */
    273    public double getInstructorProhibitedLimit() {
    274        return iInstructorProhibitedLimit;
    275    }
    276
    277    /**
    278     * When Distances.ComputeDistanceConflictsBetweenNonBTBClasses is enabled, distance limit (in minutes)
    279     * for a long travel.
    280     */
    281    public double getInstructorLongTravelInMinutes() {
    282        return iInstructorLongTravelInMinutes;
    283    }
    284
    285    /** True if legacy mode is used (Euclidian distance where 1 unit is 10 meters) */
    286    public boolean isLegacy() {
    287        return iModel == Ellipsoid.LEGACY;
    288    }
    289
    290    /** Maximal travel distance between rooms when no coordinates are given */
    291    public int getMaxTravelDistanceInMinutes() {
    292        return iMaxTravelTime;
    293    }
    294
    295    /** Add travel time between two locations */
    296    public void addTravelTime(Long roomId1, Long roomId2, Integer travelTimeInMinutes) {
    297        synchronized (iTravelTimes) {
    298            if (roomId1 == null || roomId2 == null) return;
    299            if (roomId1 < roomId2) {
    300                Map<Long, Integer> times = iTravelTimes.get(roomId1);
    301                if (times == null) { times = new HashMap<Long, Integer>(); iTravelTimes.put(roomId1, times); }
    302                if (travelTimeInMinutes == null)
    303                    times.remove(roomId2);
    304                else
    305                    times.put(roomId2, travelTimeInMinutes);
    306            } else {
    307                Map<Long, Integer> times = iTravelTimes.get(roomId2);
    308                if (times == null) { times = new HashMap<Long, Integer>(); iTravelTimes.put(roomId2, times); }
    309                if (travelTimeInMinutes == null)
    310                    times.remove(roomId1);
    311                else
    312                    times.put(roomId1, travelTimeInMinutes);
    313            }
    314        }
    315    }
    316
    317    /** Return travel time between two locations. */
    318    public Integer getTravelTimeInMinutes(Long roomId1, Long roomId2) {
    319        synchronized (iTravelTimes) {
    320            if (roomId1 == null || roomId2 == null) return null;
    321            if (roomId1 < roomId2) {
    322                Map<Long, Integer> times = iTravelTimes.get(roomId1);
    323                return (times == null ? null : times.get(roomId2));
    324            } else {
    325                Map<Long, Integer> times = iTravelTimes.get(roomId2);
    326                return (times == null ? null : times.get(roomId1));
    327            }
    328        }
    329    }
    330
    331    /** Return travel time between two locations. Travel times are used when available, use coordinates otherwise. */
    332    public Integer getDistanceInMinutes(Long roomId1, Double lat1, Double lon1, Long roomId2, Double lat2, Double lon2) {
    333        Integer distance = getTravelTimeInMinutes(roomId1, roomId2);
    334        if (distance != null) return distance;
    335
    336        if (lat1 == null || lat2 == null || lon1 == null || lon2 == null)
    337            return getMaxTravelDistanceInMinutes();
    338        else
    339            return (int) Math.min(getMaxTravelDistanceInMinutes(), Math.round(getDistanceInMeters(lat1, lon1, lat2, lon2) / iSpeed));
    340    }
    341
    342    /** Return travel distance between two locations.  Travel times are used when available, use coordinates otherwise. */
    343    public double getDistanceInMeters(Long roomId1, Double lat1, Double lon1, Long roomId2, Double lat2, Double lon2) {
    344        Integer distance = getTravelTimeInMinutes(roomId1, roomId2);
    345        if (distance != null) return minutes2meters(distance);
    346
    347        return getDistanceInMeters(lat1, lon1, lat2, lon2);
    348    }
    349
    350    /** Return travel times matrix */
    351    public Map<Long, Map<Long, Integer>> getTravelTimes() { return iTravelTimes; }
    352
    353    /**
    354     * True if distances should be considered between classes that are NOT back-to-back. Distance in minutes is then
    355     * to be compared with the difference between end of the last class and start of the second class plus break time of the first class.
    356     **/
    357    public boolean doComputeDistanceConflictsBetweenNonBTBClasses() {
    358        return iComputeDistanceConflictsBetweenNonBTBClasses;
    359    }
    360
    361    /** Few tests */
    362    public static void main(String[] args) {
    363        System.out.println("Distance between Prague and Zlin: " + new DistanceMetric().getDistanceInMeters(50.087661, 14.420535, 49.226736, 17.668856) / 1000.0 + " km");
    364        System.out.println("Distance between ENAD and PMU: " + new DistanceMetric().getDistanceInMeters(40.428323, -86.912785, 40.425078, -86.911474) + " m");
    365        System.out.println("Distance between ENAD and ME: " + new DistanceMetric().getDistanceInMeters(40.428323, -86.912785, 40.429338, -86.91267) + " m");
    366        System.out.println("Distance between Prague and Zlin: " + new DistanceMetric().getDistanceInMinutes(50.087661, 14.420535, 49.226736, 17.668856) / 60 + " hours");
    367        System.out.println("Distance between ENAD and PMU: " + new DistanceMetric().getDistanceInMinutes(40.428323, -86.912785, 40.425078, -86.911474) + " minutes");
    368        System.out.println("Distance between ENAD and ME: " + new DistanceMetric().getDistanceInMinutes(40.428323, -86.912785, 40.429338, -86.91267) + " minutes");
    369    }
    370
    371}