module FileFormat.TSPLIB(
loadTSPFile
)where
import CombinatorialOptimisation.TSP
import Data.Maybe
import qualified Data.Map as M
import qualified Data.Array as A
import Data.List
import CombinatorialOptimisation.TSP.FixedPoint
euclidianDistance :: (Double,Double)->(Double,Double)->Double
euclidianDistance (a,b) (c,d) = sqrt ((ac)*(ac)+(bd)*(bd))
euclidianDistanceCeil :: (Double,Double)->(Double,Double)->Double
euclidianDistanceCeil a b = fromIntegral . ceiling $ euclidianDistance a b
pseudoEuclideanDistance :: (Double,Double)->(Double,Double)->Double
pseudoEuclideanDistance (a,b) (c,d) = let e = sqrt (((ac)*(ac)+(bd)*(bd)) /10 )
f = fromIntegral (floor e)
in if f<e then f+1 else f
geoDistance :: (Double,Double)->(Double,Double)->Double
geoDistance (x1,y1) (x2,y2) = encodeFloat (floor dij) 0
where
q1 = cos (lon1 lon2)
q2 = cos (lat1 lat2)
q3 = cos (lat1 + lat2)
lon1 = degConvert y1
lon2 = degConvert y2
lat1 = degConvert x1
lat2 = degConvert x2
dij = 6378.388 * (acos( 0.5*((1.0+q1)*q2 ((1.0q1)*q3) )) ) + 1.0
degConvert m = let deg = encodeFloat (floor m) 0
miN = m deg
in 3.141592 * (deg + (5.0 * miN/3.0))/180.0
data Specification = IGNORE String | USEFUL String String | ENDSPEC String | FAIL String deriving Show
isUsefulSpec (USEFUL _ _) = True
isUsefulSpec _ = False
readSpecificationLine :: String->Specification
readSpecificationLine s
| likeString "NAME" s = IGNORE s
| likeString "TYPE" s = USEFUL "Type" (trim s)
| likeString "NODE_COORD_SECTION" s = ENDSPEC "NODE COORD"
| likeString "EDGE_WEIGHT_SECTION" s = ENDSPEC "EDGE WEIGHT"
| likeString "COMMENT" s = IGNORE s
| likeString "DIMENSION" s = USEFUL "Dimension" (trim s)
| likeString "DISPLAY_DATA_TYPE" s = IGNORE s
| likeString "EDGE_WEIGHT_TYPE" s = USEFUL "EdgeWeightType" (trim s)
| likeString "EDGE_WEIGHT_FORMAT" s = USEFUL "EdgeWeightFormat" (trim s)
| otherwise = FAIL $ "unrecognised field in specification : "++s
where
likeString q s = take (length q) s == q
trim s = let s' = (dropWhile (==' ')) . (drop 1) . (dropWhile (/=':')) $ s
in reverse . (dropWhile (==' ')) . reverse $ s'
readSpecification :: [String]->([Specification],[String])
readSpecification [] = ([FAIL "seem to have run out of data, without ending the specification phase"],[])
readSpecification (s:ss) = let p = readSpecificationLine s
(rs,es) = readSpecification ss
in case p of
ENDSPEC k -> ([USEFUL "DATA PART TYPE" k],ss)
IGNORE _ -> (p:rs,es)
FAIL _ -> (p:rs,es)
USEFUL _ _ -> (p:rs,es)
loadTSPFile :: InternalStorage->FilePath->IO TSPProblem
loadTSPFile storageType fName
= do rawContents<-readFile fName
let (spec,remainder) = readSpecification $ lines rawContents
let specList = map (\(USEFUL a b)->(a,b)) . filter isUsefulSpec $ spec
let dataPart = fromJust $ lookup "DATA PART TYPE" specList
let numNodes = read $ fromJust $ lookup "Dimension" specList
case dataPart of
"NODE COORD" -> return $ loadFromNodePositions storageType numNodes specList (takeWhile (/="EOF") remainder)
"EDGE WEIGHT" -> return $ loadFromMatrix storageType numNodes specList (takeWhile (/="EOF") remainder)
_ -> error "Unsupported data section"
loadFromNodePositions :: InternalStorage->Int->[(String,String)]->[String]->TSPProblem
loadFromNodePositions storageType numCities spec d
= let d' = map ((\[a,b]->(a,b)) . (map readF) . (drop 1) . words) d
posArr = A.listArray (0 , numCities1) d'
cities = [0 ..(numCities1)]
explicit = A.listArray (0,numCities*numCities1) [realToFrac $ distFunc (posArr A.! a) (posArr A.! b) | a<-cities,b<-cities]
triangular = A.listArray (0,sum [0..numCities]) [realToFrac $ distFunc (posArr A.! a) (posArr A.! b) | a<-cities,b<-[0..a]]
p = case storageType of
ExplicitMatrix -> \x y->explicit A.! (x * numCities + y)
TriangularMatrix -> (\x y->let x' = min x y; y' = max x y in triangular A.! (div (y'*y'+y') 2 + x'))
Recomputation -> \a b->realToFrac $ distFunc (posArr A.! a) (posArr A.! b)
in TSPProblem 0 M.empty p numCities M.empty M.empty
where
readF = read :: (String->Double)
distanceFunctionSpec = fromJust $ lookup "EdgeWeightType" spec
distFunc = case distanceFunctionSpec of
"GEO" -> geoDistance
"EUC_2D" -> euclidianDistance
"ATT" -> pseudoEuclideanDistance
"CEIL_2D"-> euclidianDistanceCeil
_ -> error $ "unsupported distance function : "++distanceFunctionSpec
loadFromMatrix :: InternalStorage->Int->[(String,String)]->[String]->TSPProblem
loadFromMatrix Recomputation _ _ _ = error "Cannot load matrix and store in recomputation form"
loadFromMatrix storageType numCities spec d
| distanceFunctionSpec /= "EXPLICIT" = error "loading from non explicit matrix? (does not make sense)"
| storageType == TriangularMatrix = TSPProblem 0 M.empty (\x y->let x' = min x y; y' = max x y in triangularArray A.! (div (y'*y'+y') 2 + x')) numCities M.empty M.empty
| storageType == ExplicitMatrix = TSPProblem 0 M.empty (\x y->explicitArray A.! (y * numCities + x)) numCities M.empty M.empty
where
readF = realToFrac . (read :: (String->Double))
distanceFunctionSpec = fromJust $ lookup "EdgeWeightType" spec
inputNumberSequence = concatMap ((map readF) . words) d
cities = [0 ..(numCities1)]
blankMatrix = foldl' (\m d->M.insert (d,d) 0 m) M.empty cities
fillPair m ((a,b),c) = M.insert (b,a) c $ (M.insert (a,b) c m)
makeFromTri = foldl' fillPair blankMatrix
loadedMap = case fromJust $ lookup "EdgeWeightFormat" spec of
"FULL_MATRIX" -> foldl' (\m (d,c)->M.insert d c m) M.empty (zip [(a,b) | a<-cities,b<-cities] inputNumberSequence)
"UPPER_ROW" -> makeFromTri $ zip [(a,b) | a<-cities,b<-[a..(numCities 1)],a/=b] inputNumberSequence
"LOWER_ROW" -> makeFromTri $ zip [(a,b) | a<-cities,b<-[0..a],a/=b] inputNumberSequence
"UPPER_DIAG_ROW" -> makeFromTri $ zip [(a,b) | a<-cities,b<-[a..(numCities 1)]] inputNumberSequence
"LOWER_DIAG_ROW" -> makeFromTri $ zip [(a,b) | a<-cities,b<-[0..a]] inputNumberSequence
"UPPER_COL" -> makeFromTri $ zip [(b,a) | a<-cities,b<-[a..(numCities 1)],a/=b] inputNumberSequence
"LOWER_COL" -> makeFromTri $ zip [(b,a) | a<-cities,b<-[0..a],a/=b] inputNumberSequence
"UPPER_DIAG_COL" -> makeFromTri $ zip [(b,a) | a<-cities,b<-[a..(numCities 1)]] inputNumberSequence
"LOWER_DIAG_COL" -> makeFromTri $ zip [(b,a) | a<-cities,b<-[0..a]] inputNumberSequence
explicitArray = A.listArray (0,numCities*numCities1) [loadedMap M.! (a,b) | a<-cities,b<-cities]
triangularArray = A.listArray (0,sum [0..numCities]) [loadedMap M.! (a,b) | a<-cities,b<-[0..a]]