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