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 Tomas Muller<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    
121    private final ReentrantReadWriteLock iLock = new ReentrantReadWriteLock();
122    
123    /** Default properties */
124    public DistanceMetric() {
125    }
126    
127    /** With provided ellipsoid 
128     * @param model ellipsoid model
129     **/
130    public DistanceMetric(Ellipsoid model) {
131        iModel = model;
132        if (iModel == Ellipsoid.LEGACY) {
133            iSpeed = 100.0 / 15;
134            iInstructorDiscouragedLimit = 5.0;
135            iInstructorProhibitedLimit = 20.0;
136        }
137    }
138
139    /** With provided ellipsoid and student speed
140     * @param model ellipsoid model
141     * @param speed student speed in meters per minute 
142     **/
143    public DistanceMetric(Ellipsoid model, double speed) {
144        iModel = model;
145        iSpeed = speed;
146    }
147    
148    /** Configured using properties 
149     * @param properties solver configuration
150     **/
151    public DistanceMetric(DataProperties properties) {
152        if (Ellipsoid.LEGACY.name().equals(properties.getProperty("Distances.Ellipsoid",Ellipsoid.LEGACY.name()))) {
153            //LEGACY MODE
154            iModel = Ellipsoid.LEGACY;
155            iSpeed = properties.getPropertyDouble("Student.DistanceLimit", 1000.0 / 15) / 10.0;
156            iInstructorNoPreferenceLimit = properties.getPropertyDouble("Instructor.NoPreferenceLimit", 0.0);
157            iInstructorDiscouragedLimit = properties.getPropertyDouble("Instructor.DiscouragedLimit", 5.0);
158            iInstructorProhibitedLimit = properties.getPropertyDouble("Instructor.ProhibitedLimit", 20.0);
159            iNullDistance = properties.getPropertyDouble("Distances.NullDistance", 1000.0);
160            iMaxTravelTime = properties.getPropertyInt("Distances.MaxTravelDistanceInMinutes", 60);
161        } else {
162            iModel = Ellipsoid.valueOf(properties.getProperty("Distances.Ellipsoid", Ellipsoid.WGS84.name()));
163            if (iModel == null) iModel = Ellipsoid.WGS84;
164            iSpeed = properties.getPropertyDouble("Distances.Speed", properties.getPropertyDouble("Student.DistanceLimit", 1000.0 / 15));
165            iInstructorNoPreferenceLimit = properties.getPropertyDouble("Instructor.NoPreferenceLimit", iInstructorNoPreferenceLimit);
166            iInstructorDiscouragedLimit = properties.getPropertyDouble("Instructor.DiscouragedLimit", iInstructorDiscouragedLimit);
167            iInstructorProhibitedLimit = properties.getPropertyDouble("Instructor.ProhibitedLimit", iInstructorProhibitedLimit);
168            iNullDistance = properties.getPropertyDouble("Distances.NullDistance", iNullDistance);
169            iMaxTravelTime = properties.getPropertyInt("Distances.MaxTravelDistanceInMinutes", 60);
170        }
171        iComputeDistanceConflictsBetweenNonBTBClasses = properties.getPropertyBoolean(
172                "Distances.ComputeDistanceConflictsBetweenNonBTBClasses", iComputeDistanceConflictsBetweenNonBTBClasses);
173        iInstructorLongTravelInMinutes = properties.getPropertyDouble("Instructor.InstructorLongTravelInMinutes", 30.0);
174    }
175
176    /** Degrees to radians 
177     * @param deg degrees
178     * @return radians
179     **/
180    protected double deg2rad(double deg) {
181        return deg * Math.PI / 180;
182    }
183    
184    /** Compute distance between the two given coordinates
185     * @param lat1 first coordinate's latitude
186     * @param lon1 first coordinate's longitude
187     * @param lat2 second coordinate's latitude
188     * @param lon2 second coordinate's longitude
189     * @return distance in meters
190     * @deprecated Use @{link {@link DistanceMetric#getDistanceInMeters(Long, Double, Double, Long, Double, Double)} instead (to include travel time matrix when available).
191     */
192    @Deprecated
193    public double getDistanceInMeters(Double lat1, Double lon1, Double lat2, Double lon2) {
194        if (lat1 == null || lat2 == null || lon1 == null || lon2 == null)
195            return iNullDistance;
196        
197        if (lat1.equals(lat2) && lon1.equals(lon2)) return 0.0;
198        
199        // legacy mode -- euclidian distance, 1 unit is 10 meters
200        if (iModel == Ellipsoid.LEGACY) {
201            if (lat1 < 0 || lat2 < 0 || lon1 < 0 || lon2 < 0) return iNullDistance;
202            double dx = lat1 - lat2;
203            double dy = lon1 - lon2;
204            return Math.sqrt(dx * dx + dy * dy);
205        }
206        
207        String id = null;
208        if (lat1 < lat2 || (lat1 == lat2 && lon1 <= lon2)) {
209            id =
210                Long.toHexString(Double.doubleToRawLongBits(lat1)) +
211                Long.toHexString(Double.doubleToRawLongBits(lon1)) +
212                Long.toHexString(Double.doubleToRawLongBits(lat2)) +
213                Long.toHexString(Double.doubleToRawLongBits(lon2));
214        } else {
215            id =
216                Long.toHexString(Double.doubleToRawLongBits(lat1)) +
217                Long.toHexString(Double.doubleToRawLongBits(lon1)) +
218                Long.toHexString(Double.doubleToRawLongBits(lat2)) +
219                Long.toHexString(Double.doubleToRawLongBits(lon2));
220        }
221        
222        iLock.readLock().lock();
223        try {
224            Double distance = iDistanceCache.get(id);
225            if (distance != null) return distance;
226        } finally {
227            iLock.readLock().unlock();
228        }
229        
230        iLock.writeLock().lock();
231        try {
232            Double distance = iDistanceCache.get(id);
233            if (distance != null) return distance;
234
235            double a = iModel.a(), b = iModel.b(),  f = iModel.f();  // ellipsoid params
236            double L = deg2rad(lon2-lon1);
237            double U1 = Math.atan((1-f) * Math.tan(deg2rad(lat1)));
238            double U2 = Math.atan((1-f) * Math.tan(deg2rad(lat2)));
239            double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
240            double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
241            
242            double lambda = L, lambdaP, iterLimit = 100;
243            double cosSqAlpha, cos2SigmaM, sinSigma, cosSigma, sigma, sinLambda, cosLambda;
244            do {
245              sinLambda = Math.sin(lambda);
246              cosLambda = Math.cos(lambda);
247              sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
248                (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
249              if (sinSigma==0) return 0;  // co-incident points
250              cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
251              sigma = Math.atan2(sinSigma, cosSigma);
252              double sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
253              cosSqAlpha = 1 - sinAlpha*sinAlpha;
254              cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
255              if (Double.isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (�6)
256              double C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
257              lambdaP = lambda;
258              lambda = L + (1-C) * f * sinAlpha *
259                (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
260            } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);
261            if (iterLimit==0) return Double.NaN; // formula failed to converge
262           
263            double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
264            double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
265            double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
266            double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
267              B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
268            
269            // initial & final bearings
270            // double fwdAz = Math.atan2(cosU2*sinLambda, cosU1*sinU2-sinU1*cosU2*cosLambda);
271            // double revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
272            
273            // s = s.toFixed(3); // round to 1mm precision
274
275            distance = b*A*(sigma-deltaSigma);
276            iDistanceCache.put(id, distance);
277            return distance;
278        } finally {
279            iLock.writeLock().unlock();
280        }
281    }
282    
283    /**
284     * Compute distance in minutes.
285     * 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).
286     * @param lat1 first coordinate's latitude
287     * @param lon1 first coordinate's longitude
288     * @param lat2 second coordinate's latitude
289     * @param lon2 second coordinate's longitude
290     * @return distance in minutes
291     * @deprecated Use @{link {@link DistanceMetric#getDistanceInMinutes(Long, Double, Double, Long, Double, Double)} instead (to include travel time matrix when available).
292     */
293    @Deprecated
294    public int getDistanceInMinutes(double lat1, double lon1, double lat2, double lon2) {
295        return (int) Math.round(getDistanceInMeters(lat1, lon1, lat2, lon2) / iSpeed);
296    }
297    
298    /**
299     * Converts minutes to meters.
300     * Property Distances.Speed (in meters per minute) is used, defaults to 1000 meters per 15 minutes.
301     * @param min minutes to travel
302     * @return meters to travel
303     */
304    public double minutes2meters(int min) {
305        return iSpeed * min;
306    }
307    
308
309    /** Back-to-back classes in rooms within this limit have neutral preference 
310     * @return limit in meters
311     **/
312    public double getInstructorNoPreferenceLimit() {
313        return iInstructorNoPreferenceLimit;
314    }
315
316    /** Back-to-back classes in rooms within this limit have discouraged preference 
317     * @return limit in meters
318     **/
319    public double getInstructorDiscouragedLimit() {
320        return iInstructorDiscouragedLimit;
321    }
322
323    /** Back-to-back classes in rooms within this limit have strongly discouraged preference, it is prohibited to exceed this limit.
324     * @return limit in meters 
325     **/
326    public double getInstructorProhibitedLimit() {
327        return iInstructorProhibitedLimit;
328    }
329    
330    /**
331     * When Distances.ComputeDistanceConflictsBetweenNonBTBClasses is enabled, distance limit (in minutes)
332     * for a long travel.
333     * @return travel time in minutes
334     */
335    public double getInstructorLongTravelInMinutes() {
336        return iInstructorLongTravelInMinutes;
337    }
338    
339    /** True if legacy mode is used (Euclidian distance where 1 unit is 10 meters) 
340     * @return true if the ellipsoid model is the old one
341     **/
342    public boolean isLegacy() {
343        return iModel == Ellipsoid.LEGACY;
344    }
345    
346    /** Maximal travel distance between rooms when no coordinates are given 
347     * @return travel time in minutes
348     **/
349    public int getMaxTravelDistanceInMinutes() {
350        return iMaxTravelTime;
351    }
352
353    /** Add travel time between two locations 
354     * @param roomId1 first room's id
355     * @param roomId2 second room's id
356     * @param travelTimeInMinutes travel time in minutes 
357     **/
358    public void addTravelTime(Long roomId1, Long roomId2, Integer travelTimeInMinutes) {
359        iLock.writeLock().lock();
360        try {
361            if (roomId1 == null || roomId2 == null) return;
362            if (roomId1 < roomId2) {
363                Map<Long, Integer> times = iTravelTimes.get(roomId1);
364                if (times == null) { times = new HashMap<Long, Integer>(); iTravelTimes.put(roomId1, times); }
365                if (travelTimeInMinutes == null)
366                    times.remove(roomId2);
367                else
368                    times.put(roomId2, travelTimeInMinutes);
369            } else {
370                Map<Long, Integer> times = iTravelTimes.get(roomId2);
371                if (times == null) { times = new HashMap<Long, Integer>(); iTravelTimes.put(roomId2, times); }
372                if (travelTimeInMinutes == null)
373                    times.remove(roomId1);
374                else
375                    times.put(roomId1, travelTimeInMinutes);
376            }            
377        } finally {
378            iLock.writeLock().unlock();
379        }
380    }
381    
382    /** Return travel time between two locations. 
383     * @param roomId1 first room's id
384     * @param roomId2 second room's id
385     * @return travel time in minutes
386     **/
387    public Integer getTravelTimeInMinutes(Long roomId1, Long roomId2) {
388        iLock.readLock().lock();
389        try {
390            if (roomId1 == null || roomId2 == null) return null;
391            if (roomId1 < roomId2) {
392                Map<Long, Integer> times = iTravelTimes.get(roomId1);
393                return (times == null ? null : times.get(roomId2));
394            } else {
395                Map<Long, Integer> times = iTravelTimes.get(roomId2);
396                return (times == null ? null : times.get(roomId1));
397            }
398        } finally {
399            iLock.readLock().unlock();
400        }
401    }
402    
403    /** Return travel time between two locations. Travel times are used when available, use coordinates otherwise. 
404     * @param roomId1 first room's id
405     * @param lat1 first room's latitude
406     * @param lon1 first room's longitude
407     * @param roomId2 second room's id
408     * @param lat2 second room's latitude
409     * @param lon2 second room's longitude
410     * @return distance in minutes
411     **/
412    public Integer getDistanceInMinutes(Long roomId1, Double lat1, Double lon1, Long roomId2, Double lat2, Double lon2) {
413        Integer distance = getTravelTimeInMinutes(roomId1, roomId2);
414        if (distance != null) return distance;
415        
416        if (lat1 == null || lat2 == null || lon1 == null || lon2 == null)
417            return getMaxTravelDistanceInMinutes();
418        else 
419            return (int) Math.min(getMaxTravelDistanceInMinutes(), Math.round(getDistanceInMeters(lat1, lon1, lat2, lon2) / iSpeed));
420    }
421    
422    /** Return travel distance between two locations.  Travel times are used when available, use coordinates otherwise
423     * @param roomId1 first room's id
424     * @param lat1 first room's latitude
425     * @param lon1 first room's longitude
426     * @param roomId2 second room's id
427     * @param lat2 second room's latitude
428     * @param lon2 second room's longitude
429     * @return distance in meters
430     **/
431    public double getDistanceInMeters(Long roomId1, Double lat1, Double lon1, Long roomId2, Double lat2, Double lon2) {
432        Integer distance = getTravelTimeInMinutes(roomId1, roomId2);
433        if (distance != null) return minutes2meters(distance);
434        
435        return getDistanceInMeters(lat1, lon1, lat2, lon2);
436    }
437    
438    /** Return travel times matrix
439     * @return travel times matrix
440     **/
441    public Map<Long, Map<Long, Integer>> getTravelTimes() { return iTravelTimes; }
442    
443    /**
444     * True if distances should be considered between classes that are NOT back-to-back. Distance in minutes is then 
445     * 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.
446     * @return true if distances should be considered between classes that are NOT back-to-back
447     **/
448    public boolean doComputeDistanceConflictsBetweenNonBTBClasses() {
449        return iComputeDistanceConflictsBetweenNonBTBClasses;
450    }
451    
452    /** Few tests 
453     * @param args program arguments
454     **/
455    public static void main(String[] args) {
456        System.out.println("Distance between Prague and Zlin: " + new DistanceMetric().getDistanceInMeters(50.087661, 14.420535, 49.226736, 17.668856) / 1000.0 + " km");
457        System.out.println("Distance between ENAD and PMU: " + new DistanceMetric().getDistanceInMeters(40.428323, -86.912785, 40.425078, -86.911474) + " m");
458        System.out.println("Distance between ENAD and ME: " + new DistanceMetric().getDistanceInMeters(40.428323, -86.912785, 40.429338, -86.91267) + " m");
459        System.out.println("Distance between Prague and Zlin: " + new DistanceMetric().getDistanceInMinutes(50.087661, 14.420535, 49.226736, 17.668856) / 60 + " hours");
460        System.out.println("Distance between ENAD and PMU: " + new DistanceMetric().getDistanceInMinutes(40.428323, -86.912785, 40.425078, -86.911474) + " minutes");
461        System.out.println("Distance between ENAD and ME: " + new DistanceMetric().getDistanceInMinutes(40.428323, -86.912785, 40.429338, -86.91267) + " minutes");
462    }
463
464}