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 <= limit → no preference), Instructor.DiscouragedLimit (distance <= limit → discouraged), 018 * Instructor.ProhibitedLimit (distance <= limit → 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}