$ontext


   Erwin Kalvelagen, March 2010

   Reference: Bosch, R.A.: Painting by Numbers. Optima 65, 16–17 (2001)

   Note: This one is wrong:

                 1   1      wrong
                 1   1    / /
             2 1 1 1 1 1 2 3
             3 2 1 2 1 2 3 4 8 9
       3 6   X X X   X X X X X X
       1 4   X           X X X X
     1 1 3       X   X     X X X
         2                   X X
       3 3       X X X     + X X     <<< wrong
       1 4   X           + + X X     <<< wrong
       2 5   X X       X X X X X
       2 5   X X       X X X X X
       1 1         X           X
         3       X X X


    I added the + cells to be painted black

$offtext

$if not set paintinput $set paintinput paintinput.gdx
$set havesol 0
$if set solution $set havesol 1

sets
  i 'rows'
  j 'columns'
  t 'clusters'
;

parameter ccluster(t,j);
parameter rcluster(i,t);


$gdxin %paintinput%
$load i j t
display i,j,t;
$loaddc ccluster rcluster
option ccluster:0,rcluster:0;
display ccluster,rcluster;


alias(t,tt),(i,ii),(j,jj);
set rd(t,tt) 'rev diag';
rd(t,tt)$(ord(t)+ord(tt)=card(t)+1) = yes;
display rd;


* quick check
abort$(sum((t,j),ccluster(t,j))<>sum((i,t),rcluster(i,t))) "clusters don't add up";

parameters
   rclo(i,t)  'first cell cluster (i,t) can be placed'
   rcup(i,t)  'last cell cluster (i,t) can be placed'
   cclo(j,t)  'first cell cluster (j,t) can be placed'
   ccup(j,t)  'last cell cluster (j,t) can be placed'
;

rclo(i,t) = 0; rcup(i,t) = 0;
cclo(j,t) = 0; ccup(j,t) = 0;
loop(rd(t,tt),
   rclo(i,t)$rcluster(i,t) = rclo(i,t-1)+rcluster(i,t-1)+1;
   rcup(i,tt)$rcluster(i,tt) = rcup(i,tt+1)-rcluster(i,tt)-1
                               +(card(j)+2)$(rcup(i,tt+1)=0);

   cclo(j,t)$ccluster(t,j) = cclo(j,t-1)+ccluster(t-1,j)+1;
   ccup(j,tt)$ccluster(tt,j) = ccup(j,tt+1)-ccluster(tt,j)-1
                               +(card(i)+2)$(ccup(j,tt+1)=0);
);

option rclo:0,rcup:0,cclo:0,ccup:0; display rclo,rcup,cclo,ccup;

sets
  rok(i,t,j) 'allowed values for rc(i,t,j)'
  cok(j,t,i) 'allowed values for cc(j,t,i)'
;
rok(i,t,j)$(ord(j)>=rclo(i,t) and ord(j)<=rcup(i,t)) = yes;
cok(j,t,i)$(ord(i)>=cclo(j,t) and ord(i)<=ccup(j,t)) = yes;


sets
   rcover(i,t,j,jj) 'rc(i,t,j) covers x(i,jj)'
   ccover(j,t,i,ii) 'cc(j,t,i) covers x(ii,j)'
;

rcover(rok(i,t,j),jj)$(ord(jj)>=ord(j) and ord(jj)<=ord(j)+rcluster(i,t)-1) = yes;
ccover(cok(j,t,i),ii)$(ord(ii)>=ord(i) and ord(ii)<=ord(i)+ccluster(t,j)-1) = yes;

sets
   rgap(i,t,j,jj) 'given rc(i,t,j)=1 this is where rc(i,t+1,jj) can start'
   cgap(j,t,i,ii) 'given cc(j,t,i)=1 this is where cc(j,t+1,i) can start'
   rgapok(i,t,j)  'there is at least one rgap(i,t,j,jj)'
   cgapok(j,t,i)  'there is at least one cgap(j,t,i,ii)'
;

rgap(i,t,j,jj)$(rok(i,t,j) and rok(i,t+1,jj) and ord(jj)>ord(j)+rcluster(i,t)) = yes;
cgap(j,t,i,ii)$(cok(j,t,i) and cok(j,t+1,ii) and ord(ii)>ord(i)+ccluster(t,j)) = yes;
rgapok(i,t,j)$sum(rgap(i,t,j,jj),1) = yes;
cgapok(j,t,i)$sum(cgap(j,t,i,ii),1) = yes;

variables
   x(i,j) 'cell: 0-white, 1-black'
   rc(i,t,j) 'start of row cluster (i,t)'
   cc(j,t,i) 'start of column cluster (j,t)'
   dummy
;

binary variables x,rc,cc;

equations
   rowcluster1(i,t)
   rowcluster2(i,t,j)
   colcluster(j,t)
   colcluster2(j,t,i)
   cover1(i,j) 'if not covered by row cluster => white'
   cover2(i,j) 'if not covered by col cluster => white'
   cover3(i,t,j,jj)  'if covered by row cluster => black'
   cover4(j,t,i,ii)  'if covered by col cluster => black'
   extra
   edummy
;

rowcluster1(i,t)$rcluster(i,t)..
     sum(rok(i,t,j),rc(rok)) =e= 1;

rowcluster2(i,t,j)$rgapok(i,t,j)..
     rc(i,t,j) =L= sum(rgap(i,t,j,jj),rc(i,t+1,jj));

colcluster(j,t)$ccluster(t,j)..
     sum(cok(j,t,i),cc(cok)) =e= 1;

colcluster2(j,t,i)$cgapok(j,t,i)..
     cc(j,t,i) =L= sum(cgap(j,t,i,ii),cc(j,t+1,ii));

cover1(i,jj)..
     x(i,jj) =L= sum(rcover(i,t,j,jj),rc(i,t,j));

cover2(ii,j)..
     x(ii,j) =L= sum(ccover(j,t,i,ii),cc(j,t,i));

cover3(rcover(i,t,j,jj))..
     x(i,jj) =G= rc(i,t,j);

cover4(ccover(j,t,i,ii))..
     x(ii,j) =G= cc(j,t,i);

extra..
     sum((i,j),x(i,j)) =e= sum((i,t),rcluster(i,t));

edummy..
     dummy =e= 0;


$ifthen %havesol%==1
parameter sol(i,j);
$gdxin %solution%
$load sol
option sol:0;
display sol;
x.fx(i,j) = sol(i,j);
$endif


model m /all/;
solve m minimizing dummy using mip;

option rc:0,x:0;
display x.l,rc.l;

parameter xsol(i,j);
xsol(i,j) = round(x.l(i,j));
execute_unload "paintsolution",xsol;