(!******************************************************
   Mosel Example Problems
   ======================

   file j5tax.mos
   ``````````````
   Choice of locations for income tax offices
 
   A graph representing the cities and major roads connecting
   them is provided. The population of each city and the 
   distance between two cities directly connected are also 
   known from the graph. Where should tax offices be located 
   to minimize the average distance per inhabitant to the 
   closest tax office?

   Given the graph, the distance matrix is calculated using a
   specialized algorithm (Floyd-Warshall), implemented by the
   procedure 'calculatedist'. These distances are used 
   in the MIP model where two families of decision variables 
   are defined.

   (c) 2008-2022 Fair Isaac Corporation
       author: S. Heipcke, Mar. 2002, rev. Mar. 2022
*******************************************************!)

model "J-5 Tax office location"
 uses "mmxprs"

 forward procedure calculatedist

 declarations
  CITIES = 1..12                        ! Set of cities

  DIST: array(CITIES,CITIES) of integer ! Distance matrix
  POP: array(CITIES) of integer         ! Population of cities
  LEN: dynamic array(CITIES,CITIES) of integer ! Road lengths
  NUMLOC: integer                       ! Desired number of tax offices
  
  build: array(CITIES) of mpvar         ! 1 if office in city, 0 otherwise
  depend: array(CITIES,CITIES) of mpvar ! (c,d) 1 if city c depends on office
                                        ! in city d, 0 otherwise
 end-declarations

 initializations from 'j5tax.dat'
  LEN POP NUMLOC
 end-initializations

! Calculate the distance matrix
 calculatedist    

! Objective: weighted total distance
 TotDist:= sum(c,d in CITIES) POP(c)*DIST(c,d)*depend(c,d)

! Assign cities to offices
 forall(c in CITIES) sum(d in CITIES) depend(c,d) = 1
 
! Limit total number of offices
 sum(c in CITIES) build(c) <= NUMLOC
 
! Relations between dependencies and offices built
 forall(c,d in CITIES) depend(c,d) <= build(d)  
 
 forall(c in CITIES) build(c) is_binary

! Solve the problem
 minimize(TotDist)
 
! Solution printing
 writeln("Total weighted distance: ", getobjval, 
        " (average per inhabitant: ", getobjval/sum(c in CITIES) POP(c), ")")
 forall(c in CITIES | getsol(build(c))>0) do
  write("Office in ", c, ": ")
  forall(d in CITIES | getsol(depend(d,c))>0) write(" ", d)
  writeln
 end-do

!-----------------------------------------------------------------

! Calculate the distance matrix using Floyd-Warshall algorithm
 procedure calculatedist
 ! Initialize all distance labels with a sufficiently large value
  BIGM:=sum(c,d in CITIES | exists(LEN(c,d))) LEN(c,d)
  forall(c,d in CITIES) DIST(c,d):=BIGM

  forall(c in CITIES) DIST(c,c):=0    ! Set values on the diagonal to 0
 
! Length of existing road connections
  forall(c,d in CITIES | exists(LEN(c,d))) do
   DIST(c,d):=LEN(c,d)
   DIST(d,c):=LEN(c,d)
  end-do
 
! Update shortest distance for every node triple 
  forall(b,c,d in CITIES | c<d )
   if DIST(c,d) > DIST(c,b)+DIST(b,d) then
    DIST(c,d):= DIST(c,b)+DIST(b,d)
    DIST(d,c):= DIST(c,b)+DIST(b,d)
   end-if 
 end-procedure 
 
end-model