001    package net.sf.cpsolver.ifs.util;
002    
003    import java.util.HashMap;
004    import 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 Tomas Muller<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     */
038    public 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    }