input boxes;
input mp-tool;
input marctool;

numeric csptextscale,           % scaling factor (relative to font size).
        cspvaluepenscale,       % penwidth for circles around values in domain.
        cspdomainpenscale,      % penwidth for ellipses around domains.
        cspdomaindashscale,     % scale for ellipses (dashed lines) around domains.
        cspconstraintpenscale,  % penwidth for lines between values in constraint.
        cspconstraintdashscale; % scale for (dashed) lines between values in constrain.

def scaleinitialisation =
  csptextscale           := 1.0;
  cspvaluepenscale       := 1pt;
  cspdomainpenscale      := 1.5pt;
  cspdomaindashscale     := 2.5;
  cspconstraintpenscale  := 1.1pt;
  cspconstraintdashscale := 2;
enddef;

boolean rotatetext;     % true if text has to be rotated with domain.
numeric nrvars;         % Records the number of variables.
numeric rotn[];         % rotation of domain of constraint.
numeric card[];         % card[ i ] is cardinality of i-th variable.
boolean uuniv[];        % uuniv[a] iff
                        % unary constrain on a-th variable is universal.
boolean buniv[][];      % buniv[a][b] iff
                        % constraint between a-th and b-th variable is universal.
boolean uctr[][];       % bctr[a][v] if 
                        %  v in unary constraint of a-th variable.
boolean bctr[][][][];   % bctr[a][b][v][w] if
                        %  (v,w) in constraint between a-th and b-th variable.
boolean alive[][];      % alive[a][v] iff
                        % v is still in domain of a-th variable.
picture names[];        % records labels of variable names.
picture values[][];     % values[a][v] records
                        % the picture of value v of a-th variable.
path    vcircles[][];   % vcircles[a][v] stores
                        % circle around value v of a-th variable.
color   fcircles[][];   % fcircles[a][v] stores
                        % the fillcolour of the circle of value v of variable a.
color   dcircles[][];   % dcircles[a][v] stores
                        % the drawcolour of the circle of value v of variable a.
color   tcircles[][];   % tcircles[a][v] stores
                        % the drawcolour of value v of variable a.
color   cline[][][][];  % cline[a][b][v][w] is the drawcolour of line
                        % (v,w) in constraint between a and b.
color   namecolour[];   % colour of the name of the domain.
color   bcolour[];      % colours of boundaries of domains.
path    boundaries[];   % boundaries[a] is boundary of domain of a-th variable.
pair    vcentres[][];   % vcentres[a][v] stores 
                        % centre of circle around value v of a-th variable.
pair    vshift;         % difference between to adjacent values in domain.
boolean hidden[][];     % hidden[a][b] is true if and only if the constraint

scaleinitialisation;
rotatetext := false;
nrvars     := 0;
vshift     := (3mm,0);

def initcsp( expr rotationmode ) text options =
  scaleinitialisation;
  rotatetext := rotationmode;
  nrvars     := 0;
  vshift     := (3mm,0);
  for option = options:
    if option = "lecture" :
      penrescaling( 1.25 );
      dashrescaling( 0.8 );
    fi;
  endfor;
enddef;

def penrescaling( expr s ) =
  cspvaluepenscale       := s * cspvaluepenscale;
  cspdomainpenscale      := s * cspdomainpenscale;
  cspconstraintpenscale  := s * cspconstraintpenscale;
enddef;

def dashrescaling( expr s ) =
  cspdomaindashscale     := s * cspdomaindashscale;
  cspconstraintdashscale := s * cspconstraintdashscale;
enddef;

vardef inconstraint( suffix $, $$ )( expr c ) =
  (bctr[ $.number ][ $$.number ][ xpart( c ) ][ ypart( c ) ])
enddef;

vardef boundaryofdomain( suffix $ )( expr h, rot ) =
  (numericboundaryofdomain( $.number, h, rot ))
enddef;

vardef numericboundaryofdomain( expr nr, h, rot ) =
  save c, p;
  numeric c;
  pair p[];
  c  = card[ nr ];
  p1 = vcentres[ nr ][1];
  p2 = vcentres[ nr ][c];
  p3 = (0,0.5h+0.7xpart( vshift )) rotated rot;
  (p1 + (p3 rotated 180)..
   p1 + (p3 rotated 90)..
   p1 + p3--
   p1 + p3--
   p2 + p3--
   p2 + p3..
   p2 + (p3 rotated 270)..
   p2 + (p3 rotated 180)--
   p2 + (p3 rotated 180)--
   cycle)
enddef;

boolean labelleft, labelright;
labelleft  = true;
labelright = false;

vardef domain@#( expr n, pos, rot, labelpos )( text vals_ ) =
  _n_ := str @#;
  generic_declare(numeric) _n.number;
  @#number := nrvars + 1;
  numericdomain( nrvars + 1, n, pos, rot, labelpos, vals_ );
enddef;

def numericdomain( expr nr, n, pos, rot, labelpos )( text vals_ ) =
  nrvars := nr;
  namecolour[ nr ] := black;
  bcolour[ nr ] := black;
  rotn[ nr ] := rot;
  begingroup
    save np_, p_, h_, c_, vv_, m_;
    picture vv_;
    pair np_, p_;
    numeric h_, c_, m_;
    p_ := (0,0);
    h_ := 0;
    c_ := 0;
    % First we need to know the distance required to display the values.
    for v_ = vals_:
      c_ := c_ + 1;
      m_ := 1.5 * csptextscale * max( heightofbox( v_ ), widthofbox( v_ ) );
      p_ := p_ + (m_, 0) + vshift;
      h_ := max( h_, m_ );
    endfor;
    card[ nr ] := c_;
    for i = 1 upto c_:
      uctr[ nr ][ i ] := true;
    endfor;
    c_  := 0;
    np_ := pos - (0.5(p_-vshift) rotated rot);
    p_  := np_;
    for v_ = vals_:
      vv_                  := v_ scaled csptextscale;
      c_                   := c_ + 1;
      fcircles[ nr ][ c_ ] := (-1,-1,-1);
      dcircles[ nr ][ c_ ] := black;
      tcircles[ nr ][ c_ ] := black;
      values[ nr ][ c_ ]   := vv_ rotated rot
                                  shifted
                                  (p_ + ((0.5*(h_-widthofbox( vv_ )),0)
                                         rotated rot
                                        ));
      vcentres[ nr ][ c_ ] := p_ + ((0.5h_,0.3h_) rotated rot);
      if rotatetext = false:
        values[ nr ][ c_ ] := values[ nr ][ c_ ]
                              rotatedabout( vcentres[ nr ][ c_ ], 360 - rot );
      fi;
      vcircles[ nr ][ c_ ] := thescaledcircle( vcentres[ nr ][ c_ ]
                                             , h_ + 0.4*xpart(vshift)
                                             );
      p_                   := p_ + (((h_, 0) + vshift) rotated rot);
      alive[ nr ][ c_ ]    := true;
    endfor;
    boundaries[ nr ] := numericboundaryofdomain( nr, h_, rot );
    names[ nr ] := nullpicture;
    if labelpos = labelleft:
      addto names[ nr ]
            also ((n scaled (1.5*csptextscale) rotated rot )
                  shifted (np_ + (((-1.5*csptextscale*widthofbox( n ), 0)
                                   - 1.5vshift
                                 ) rotated rot)));
    else:
      addto names[ nr ]
            also ((n scaled (1.5*csptextscale) rotated rot )
                  shifted (p_ + (((1.0*csptextscale*widthofbox( n ), 0)
                                  - 0.5vshift
                                 )
                                rotated rot)));
    fi;
    if rotatetext = false:
      names[ nr ] := names[ nr ]
                     rotatedabout( centreofbox( names[ nr ] ), 360 - rot );
    fi;
  endgroup;
  for i = 1 upto (nr-1):
    buniv[ nr ][ i  ]:= true;
    buniv[ i  ][ nr ]:= true;
  endfor;
  uuniv[ nr ] := true;
enddef;

def removefromdomain( suffix $ )( text vals_ ) =
  numericremovefromdomain( $.number, vals_ );
enddef;

def numericremovefromdomain( expr nr )( text vals_ ) =
  for v = vals_:
    if (0 < v) and (v <= card[ nr ]):
      alive[ nr ][ v ] := false;
    fi;
  endfor;
enddef;

def newconstraint( suffix $, $$ ) =
  numericnewconstraint( $.number, $$.number );
enddef;

def numericnewconstraint( expr na, nb ) =
  buniv[ na ][ nb ] := false;
  buniv[ nb ][ na ] := false;
  for v = 1 upto card[ na ]:
    for w = 1 upto card[ nb ]:
      bctr[ na ][ nb ][v][w] := true;
      bctr[ nb ][ na ][w][v] := true;
      cline[ na ][ nb ][v][w] := black;
      cline[ nb ][ na ][w][v] := black;
    endfor;
  endfor;
enddef;

def equality( suffix $, $$ ) =
  numericequality( $.number, $$.number );
enddef;

def numericequality( expr na, nb) =
  begingroup
    save e;
    boolean e;
    if buniv[ na ][ nb ] = true:
      numericnewconstraint( na, nb );
    fi;
    for v = 1 upto card[ na ]:
      for w = 1 upto card[ nb ]:
        e := bctr[ na ][ nb ][v][w] and (v = w);
        bctr[ na ][ nb ][v][w] := e;
        bctr[ nb ][ na ][w][v] := e;
        cline[ na ][ nb ][v][w] := black;
        cline[ nb ][ na ][w][v] := black;
      endfor;
    endfor;
  endgroup;
enddef;

def inequality( suffix $, $$ ) =
  numericinequality( $.number, $$.number );
enddef;

def numericinequality( expr na, nb ) =
  begingroup
    save e;
    boolean e;
    if buniv[ na ][ nb ] = true:
      numericnewconstraint( na, nb );
    fi;
    for v = 1 upto card[ na ]:
      for w = 1 upto card[ nb ]:
        e := bctr[ na  ][ nb ][v][w] and (v <> w);
        bctr[ na ][ nb ][v][w] := e;
        bctr[ nb ][ na ][w][v] := e;
        cline[ na ][ nb ][v][w] := black;
        cline[ nb ][ na ][w][v] := black;
      endfor;
    endfor;
  endgroup;
enddef;

def addunaryconstraint( suffix $ )( text values_ ) =
  numericaddunaryconstraint( $.number, values_ );
enddef;

def numericaddunaryconstraint( expr na )( text values_ ) =
  uuniv[ na ] := false;
  begingroup
    for v = 1 upto card[ na ]:
      uctr[ na ][ v ] := false;
    endfor;
    for vv = values_:
      if (vv < 1) or (card[ na ] < vv):
        error( "numericaddbinaryconstraint" );
      fi;
      for v = 1 upto card[ na ]:
        if vv = v: uctr[ na ][ v ] := true; fi;
      endfor;
    endfor;
  endgroup;
enddef;

def addconstraint( suffix $, $$ )( text pairs_ ) =
  numericaddconstraint( $.number, $$.number, pairs_ );
enddef;

def numericaddconstraint( expr na, nb )( text pairs_ ) =
  begingroup
    save v, w, f;
    boolean f;
    if buniv[ na ][ nb ] = true:
      numericnewconstraint( na, nb );
    fi;
    for vv = 1 upto card[ na ]:
      for ww = 1 upto card[ nb ]:
        if bctr[ na ][ nb ][ vv ][ ww ]:
          f := false;
          for p = pairs_:
            v := xpart( p );
            w := ypart( p );
            if (v < 1) or (card[ na ] < v) or (w < 1) or (card[ nb ] < w):
               error( "numericaddconstraint" );
            fi;
            if (v=vv) and (w=ww):
              f := true;
            fi;
          endfor;
          bctr[ na ][ nb ][vv][ww] := f;
          bctr[ nb ][ na ][ww][vv] := f;
          cline[ na ][ nb ][vv][ww] := black;
          cline[ nb ][ na ][ww][vv] := black;
        fi;
      endfor;
    endfor;
  endgroup;
enddef;

def drawdomain( suffix $ ) =
  numericdrawdomain( $.number );
enddef;

def numericdrawdomain( expr nr ) =
  draw thenumericdomain( nr );
enddef;

vardef thenumericdomain( expr nr ) =
  save pic;
  picture pic;
  pic := thenumericvalues( nr );
  addto pic also thenumericboundary( nr );
  addto pic also thenumericname( nr );
  pic
enddef;

vardef thevaluecircle( suffix $ )( expr val ) =
  thenumericvaluecircle( $.number )( val )
enddef;

vardef thenumericvaluecircle( expr num )( expr val ) =
  (vcircles[ num ][ val ])
enddef;

vardef thenumericvalues( expr nr ) =
  save pic;
  picture pic;
  pic := nullpicture;
  for val = 1 upto card[ nr ]:
    if alive[ nr ][ val ]:
      if fcircles[ nr ][ val ] <> (-1,-1,-1):
        addto pic contour
              vcircles[ nr ][ val ]
              withcolor fcircles[ nr ][ val ];
      fi;
      addto pic doublepath
            vcircles[ nr ][ val ]
            withcolor dcircles[ nr ][ val ]
            withpen pencircle scaled cspvaluepenscale;
      addto pic also values[ nr ][ val ]
            withcolor tcircles[ nr ][ val ];
    fi;
  endfor;
  pic
enddef;

vardef thenumericboundary( expr nr ) =
  save pic;
  picture pic;
  pic := nullpicture;
  addto pic doublepath
        boundaries[ nr ]
        withpen pencircle scaled cspdomainpenscale
        withcolor bcolour[ nr ]
        dashed evenly scaled cspdomaindashscale;
  pic
enddef;

vardef thenumericname( expr nr ) =
  save pic;
  picture pic;
  pic := nullpicture;
  addto pic also names[ nr ]
        withcolor namecolour[ nr ];
  pic
enddef;

def drawbinaryelation( suffix $, $$ ) =
  numericdrawbinaryrelation( $.number, $$.number );
enddef;

def numericdrawbinaryrelation( expr na, nb ) =
  draw thenumericbinaryrelation( na, nb );
enddef;

vardef thenumericbinaryrelation( expr na, nb ) =
  save pic, v, w;
  numeric v, w;
  picture pic;
  pic := nullpicture;
  if not buniv[ na ][ nb ]:
    for v = 1 upto card[ na ]:
      if alive[ na ][v]:
        for w = 1 upto card[ nb ]:
          if alive[ nb ][w] and bctr[ na ][ nb ][v][w]:
            addto pic doublepath
                  ( (vcentres[ na ][v]--vcentres[ nb ][w])
                    cutbefore vcircles[ na ][v]
                    cutafter  vcircles[ nb ][w] )
                    withpen pencircle scaled cspconstraintpenscale
                    withcolor cline[na][nb][v][w]
                    dashed evenly scaled cspconstraintdashscale;
          fi;
        endfor;
      fi;
    endfor;
  fi;
  pic
enddef;

def drawunaryconstraint( suffix $ ) =
  numericdrawunaryconstraint( $.number );
enddef;

def numericdrawunaryconstraint( expr number ) =
  draw thenumericunaryconstraint( number );
enddef;

vardef thenumericunaryconstraint( expr nr ) =
  save pic;
  picture pic;
  pic := nullpicture;
  if not uuniv[ nr ]:
    for v = 1 upto card[ nr ]:
      if alive[ nr ][ v ] and uctr[ nr ][ v ]:
        addto pic doublepath
              (((vcentres[ nr ][ v ]){dir (rotn[ nr ] + 50)}..
                  {dir (rotn[ nr ] + 180)}(vcentres[ nr ][ v ] +
                                           ((0,csptextscale * cm)
                                            rotated rotn[ nr ]
                                           ))
                  {dir (rotn[ nr ] + 180)}..
                  {dir (rotn[ nr ] + 310)}(vcentres[ nr ][ v ]))
                  cutbefore vcircles[ nr ][ v ]
                  cutafter  vcircles[ nr ][ v ])
                  dashed evenly
                  withpen pencircle scaled cspconstraintpenscale;
      fi;
    endfor;
  fi;
  pic
enddef;

vardef thecsp =
  save first, last, pic;
  picture pic;
  numeric first, last;
  pic := nullpicture;
  first := 1;
  last  := nrvars;
  for a = first upto last:
    addto pic also thenumericboundary( a );
  endfor;
  for a = first upto last:
    addto pic also thenumericunaryconstraint( a );
    for b = a + 1 upto nrvars:
      addto pic also thenumericbinaryrelation( a, b );
    endfor;
  endfor;
  for a = first upto last:
    addto pic also thenumericvalues( a );
    addto pic also thenumericname( a );
  endfor;
  pic
enddef;

vardef drawcsp =
  draw thecsp;
enddef;

def colourvalues( suffix $ )( expr fc, dc, tc )( text vs_ ) =
  numericcolourvalues( $.number )( fc, dc, tc )( vs_ );
enddef;

def numericcolourvalues( expr a, fc, dc, tc)( text vs_ ) =
  for v = vs_:
    fcircles[ a ][ v ] := fc;
    dcircles[ a ][ v ] := dc;
    tcircles[ a ][ v ] := tc;
  endfor;
enddef;

def uncolourvalues( suffix $ )( text vs_ ) =
  numericuncolourvalues( $.number )( vs_ );
enddef;

def numericuncolourvalues( expr a )( text vs_ ) =
  for v = vs_:
    fcircles[ a ][ v ] := (-1,-1,-1);
    dcircles[ a ][ v ] := black;
    tcircles[ a ][ v ] := black;
  endfor;
enddef;

def colourlines( suffix $, $$ )( expr c )( text ps_ ) =
  numericcolourlines( $.number, $$.number, c )( ps_ );
enddef;

def numericcolourlines( expr a, b, c )( text ps_ ) =
  for p = ps_:
    cline[ a ][ b ][ xpart( p ) ][ ypart( p ) ] := c;
    cline[ b ][ a ][ ypart( p ) ][ xpart( p ) ] := c;
  endfor;
enddef;

def colourname( suffix $ )( expr c ) =
  numericcolourname( $.number, c );
enddef;

def numericcolourname( expr a, c ) =
  namecolour[ a ] := c;
enddef;

def colourboundary( suffix $ )( expr c ) =
  numericcolourboundary( $.number, c );
enddef;

def numericcolourboundary( expr a, c ) =
  bcolour[ a ] := c;
enddef;

%
% Crude way---believe me, I know---to inforce arc-consistency.
%
def arcconsistency =
  begingroup
    save changed, supported;
    boolean changed, supported;
    changed := true;
    forever: exitunless changed;
      changed := false;
      for a = 1 upto nrvars:
        for b = 1 upto nrvars:
          if (a <> b):
            if (not (buniv[a][b])):
              for v = 1 upto card[a]:
                if alive[a][v]:
                  supported := false;
                  for w = 1 upto card[b]:
                    if alive[b][w] and bctr[a][b][v][w]:
                      supported := true;
                    fi;
                  endfor;
                  if not supported:
                    changed     := true;
                    alive[a][v] := false;
                  fi;
                fi;
              endfor;
            fi;
          fi;
        endfor;
      endfor;
    endfor;
  endgroup
enddef;