model ModelName
  options noimplicit
  uses "mmxprs"

  parameters
    N=5
  end-parameters

  declarations
    POINTS=1..N
    SIDES=1..4 ! 1 = top, 2 = right, 3 = bottom, 4 = left
    x: array(POINTS) of mpvar
    y: array(POINTS) of mpvar
    side: array(POINTS,SIDES) of mpvar
    distances: array(POINTS,POINTS) of mpvar
    Objective:linctr
  end-declarations

  procedure linearization(x: linctr, y: linctr, b: mpvar, M: real)
    y <= M*b
    y <= x
    x + M*b - M <= y
  end-procedure

  procedure linearization2(x: linctr, y: linctr, b: mpvar, M: real)
    y <= M*b
    y >= -M*b
    x + M*b - M <= y
    x - M*b + M >= y
  end-procedure

  procedure absmodel(x: linctr, y: mpvar, M: real)

    declarations
      b: mpvar
      xTimesB: mpvar
    end-declarations

    b is_binary

    y = 2*xTimesB - x
    xTimesB >= 0
    x <= b*M

    linearization2(x,xTimesB,b,M)

  end-procedure

  forall(p in POINTS, s in SIDES) side(p,s) is_binary

  forall(p in POINTS) do
    sum(s in SIDES) side(p,s) = 1
  end-do

  forall(p in POINTS) do
    x(p) <= 1
    y(p) <= 1
  end-do

  declarations
    side1TimesY: array(POINTS) of mpvar
    side2TimesX: array(POINTS) of mpvar
    side3TimesY: array(POINTS) of mpvar
    side4TimesX: array(POINTS) of mpvar
  end-declarations

  forall(p in POINTS) do

    ! top:
    side1TimesY(p) = side(p,1)
    linearization(y(p),side1TimesY(p),side(p,1),1)

    ! right:
    side2TimesX(p) = side(p,2)
    linearization(x(p),side2TimesX(p),side(p,2),1)

    ! bottom:
    side3TimesY(p) = 0
    linearization(y(p),side3TimesY(p),side(p,3),1)

    ! left:
    side4TimesX(p) = 0
    linearization(x(p),side4TimesX(p),side(p,4),1)

  end-do

  declarations
    xDifference: array(POINTS,POINTS) of mpvar
    yDifference: array(POINTS,POINTS) of mpvar
  end-declarations

  forall(p1 in POINTS, p2 in POINTS) do
    distances(p1,p2) = xDifference(p1,p2)+yDifference(p1,p2)

    absmodel(x(p1)-x(p2), xDifference(p1,p2), 1)
    absmodel(y(p1)-y(p2), yDifference(p1,p2), 1)
  end-do

  Objective := sum(p1 in POINTS, p2 IN POINTS) distances(p1,p2)

  maximise(Objective)

  writeln("Total distances: ", getobjval)

  forall(p in POINTS) do
    writeln("Point ", p, ": (", getsol(x(p)), " ", getsol(y(p)), ")")
  end-do

end-model