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}