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 HashMap<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 Double distance = iDistanceCache.get(id);
192
193 if (distance == null) {
194 double a = iModel.a(), b = iModel.b(), f = iModel.f(); // ellipsoid params
195 double L = deg2rad(lon2-lon1);
196 double U1 = Math.atan((1-f) * Math.tan(deg2rad(lat1)));
197 double U2 = Math.atan((1-f) * Math.tan(deg2rad(lat2)));
198 double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
199 double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
200
201 double lambda = L, lambdaP, iterLimit = 100;
202 double cosSqAlpha, cos2SigmaM, sinSigma, cosSigma, sigma, sinLambda, cosLambda;
203 do {
204 sinLambda = Math.sin(lambda);
205 cosLambda = Math.cos(lambda);
206 sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
207 (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
208 if (sinSigma==0) return 0; // co-incident points
209 cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
210 sigma = Math.atan2(sinSigma, cosSigma);
211 double sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
212 cosSqAlpha = 1 - sinAlpha*sinAlpha;
213 cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
214 if (Double.isNaN(cos2SigmaM)) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (�6)
215 double C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
216 lambdaP = lambda;
217 lambda = L + (1-C) * f * sinAlpha *
218 (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
219 } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);
220 if (iterLimit==0) return Double.NaN; // formula failed to converge
221
222 double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
223 double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
224 double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
225 double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
226 B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
227
228 // initial & final bearings
229 // double fwdAz = Math.atan2(cosU2*sinLambda, cosU1*sinU2-sinU1*cosU2*cosLambda);
230 // double revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
231
232 // s = s.toFixed(3); // round to 1mm precision
233
234 distance = b*A*(sigma-deltaSigma);
235 iDistanceCache.put(id, distance);
236 }
237
238 return distance;
239 }
240
241 /**
242 * Compute distance in minutes.
243 * 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).
244 * @deprecated Use @{link {@link DistanceMetric#getDistanceInMinutes(Long, Double, Double, Long, Double, Double)} instead (to include travel time matrix when available).
245 */
246 @Deprecated
247 public int getDistanceInMinutes(double lat1, double lon1, double lat2, double lon2) {
248 return (int) Math.round(getDistanceInMeters(lat1, lon1, lat2, lon2) / iSpeed);
249 }
250
251 /**
252 * Converts minutes to meters.
253 * Property Distances.Speed (in meters per minute) is used, defaults to 1000 meters per 15 minutes.
254 */
255 public double minutes2meters(int min) {
256 return iSpeed * min;
257 }
258
259
260 /** Back-to-back classes in rooms within this limit have neutral preference */
261 public double getInstructorNoPreferenceLimit() {
262 return iInstructorNoPreferenceLimit;
263 }
264
265 /** Back-to-back classes in rooms within this limit have discouraged preference */
266 public double getInstructorDiscouragedLimit() {
267 return iInstructorDiscouragedLimit;
268 }
269
270 /** Back-to-back classes in rooms within this limit have strongly discouraged preference, it is prohibited to exceed this limit. */
271 public double getInstructorProhibitedLimit() {
272 return iInstructorProhibitedLimit;
273 }
274
275 /**
276 * When Distances.ComputeDistanceConflictsBetweenNonBTBClasses is enabled, distance limit (in minutes)
277 * for a long travel.
278 */
279 public double getInstructorLongTravelInMinutes() {
280 return iInstructorLongTravelInMinutes;
281 }
282
283 /** True if legacy mode is used (Euclidian distance where 1 unit is 10 meters) */
284 public boolean isLegacy() {
285 return iModel == Ellipsoid.LEGACY;
286 }
287
288 /** Maximal travel distance between rooms when no coordinates are given */
289 public int getMaxTravelDistanceInMinutes() {
290 return iMaxTravelTime;
291 }
292
293 /** Add travel time between two locations */
294 public void addTravelTime(Long roomId1, Long roomId2, Integer travelTimeInMinutes) {
295 if (roomId1 == null || roomId2 == null) return;
296 if (roomId1 < roomId2) {
297 Map<Long, Integer> times = iTravelTimes.get(roomId1);
298 if (times == null) { times = new HashMap<Long, Integer>(); iTravelTimes.put(roomId1, times); }
299 if (travelTimeInMinutes == null)
300 times.remove(roomId2);
301 else
302 times.put(roomId2, travelTimeInMinutes);
303 } else {
304 Map<Long, Integer> times = iTravelTimes.get(roomId2);
305 if (times == null) { times = new HashMap<Long, Integer>(); iTravelTimes.put(roomId2, times); }
306 if (travelTimeInMinutes == null)
307 times.remove(roomId1);
308 else
309 times.put(roomId1, travelTimeInMinutes);
310 }
311 }
312
313 /** Return travel time between two locations. */
314 public Integer getTravelTimeInMinutes(Long roomId1, Long roomId2) {
315 if (roomId1 == null || roomId2 == null) return null;
316 if (roomId1 < roomId2) {
317 Map<Long, Integer> times = iTravelTimes.get(roomId1);
318 return (times == null ? null : times.get(roomId2));
319 } else {
320 Map<Long, Integer> times = iTravelTimes.get(roomId2);
321 return (times == null ? null : times.get(roomId1));
322 }
323 }
324
325 /** Return travel time between two locations. Travel times are used when available, use coordinates otherwise. */
326 public Integer getDistanceInMinutes(Long roomId1, Double lat1, Double lon1, Long roomId2, Double lat2, Double lon2) {
327 Integer distance = getTravelTimeInMinutes(roomId1, roomId2);
328 if (distance != null) return distance;
329
330 if (lat1 == null || lat2 == null || lon1 == null || lon2 == null)
331 return getMaxTravelDistanceInMinutes();
332 else
333 return (int) Math.min(getMaxTravelDistanceInMinutes(), Math.round(getDistanceInMeters(lat1, lon1, lat2, lon2) / iSpeed));
334 }
335
336 /** Return travel distance between two locations. Travel times are used when available, use coordinates otherwise. */
337 public double getDistanceInMeters(Long roomId1, Double lat1, Double lon1, Long roomId2, Double lat2, Double lon2) {
338 Integer distance = getTravelTimeInMinutes(roomId1, roomId2);
339 if (distance != null) return minutes2meters(distance);
340
341 return getDistanceInMeters(lat1, lon1, lat2, lon2);
342 }
343
344 /** Return travel times matrix */
345 public Map<Long, Map<Long, Integer>> getTravelTimes() { return iTravelTimes; }
346
347 /**
348 * True if distances should be considered between classes that are NOT back-to-back. Distance in minutes is then
349 * 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.
350 **/
351 public boolean doComputeDistanceConflictsBetweenNonBTBClasses() {
352 return iComputeDistanceConflictsBetweenNonBTBClasses;
353 }
354
355 /** Few tests */
356 public static void main(String[] args) {
357 System.out.println("Distance between Prague and Zlin: " + new DistanceMetric().getDistanceInMeters(50.087661, 14.420535, 49.226736, 17.668856) / 1000.0 + " km");
358 System.out.println("Distance between ENAD and PMU: " + new DistanceMetric().getDistanceInMeters(40.428323, -86.912785, 40.425078, -86.911474) + " m");
359 System.out.println("Distance between ENAD and ME: " + new DistanceMetric().getDistanceInMeters(40.428323, -86.912785, 40.429338, -86.91267) + " m");
360 System.out.println("Distance between Prague and Zlin: " + new DistanceMetric().getDistanceInMinutes(50.087661, 14.420535, 49.226736, 17.668856) / 60 + " hours");
361 System.out.println("Distance between ENAD and PMU: " + new DistanceMetric().getDistanceInMinutes(40.428323, -86.912785, 40.425078, -86.911474) + " minutes");
362 System.out.println("Distance between ENAD and ME: " + new DistanceMetric().getDistanceInMinutes(40.428323, -86.912785, 40.429338, -86.91267) + " minutes");
363 }
364
365 }