Resources

Balanced Incomplete Block Design (clpe13)

/***************************************************************/
/*                                                             */
/*          S A S   S A M P L E   L I B R A R Y                */
/*                                                             */
/*    NAME: clpe13                                             */
/*   TITLE: Balanced Incomplete Block Design (clpe13)          */
/* PRODUCT: OR                                                 */
/*  SYSTEM: ALL                                                */
/*    KEYS: OR                                                 */
/*   PROCS: CLP                                                */
/*    DATA:                                                    */
/*                                                             */
/* SUPPORT:                             UPDATE:                */
/*     REF:                                                    */
/*    MISC: Example 13 from the CLP Procedure chapter of the   */
/*          Constraint Programming book.                       */
/*                                                             */
/***************************************************************/


%macro bibd(v, b, r, k, lambda, out=bibdout);
   /* Arrange v objects into b blocks such that:
         (i) each object occurs in exactly r blocks,
         (ii) each block contains exactly k objects,
         (iii) every pair of objects occur together in exactly lambda blocks.

      Equivalently, create a binary matrix with v rows and b columns,
      with r 1s per row, k 1s per column,
      and scalar product lambda between any pair of distinct rows.
   */

   /* Check necessary conditions */
   %if (%eval(&r * &v) ne %eval(&b * &k)) or
      (%eval(&lambda * (&v - 1)) ne %eval(&r * (&k - 1))) or
      (&v > &b) %then %do;
      %put BIBD necessary conditions are not met.;
      %goto EXIT;
   %end;

   proc clp out=&out(keep=x:) domain=[0,1] varselect=FIFO;
      /* Decision variables: */
      /* Decision variable X_i_c = 1 iff object i occurs in block c. */
      var (
           %do i=1 %to &v;
              x&i._1-x&i._&b.
           %end;
          ) = [0,1];

      /* Mandatory constraints: */
      /* (i) Each object occurs in exactly r blocks. */
      %let q = %eval(&b.-&r.);  /* each row has &q 0s and &r 1s */
      %do i=1 %to &v;
         gcc( x&i._1-x&i._&b. ) = ((0,0,&q.) (1,0,&r.));
      %end;

      /* (ii) Each block contains exactly k objects. */
      %let h = %eval(&v.-&k.);  /* each column has &h 0s and &k 1s */
      %do c=1 %to &b;
         gcc(
             %do i=1 %to &v;
                x&i._&c.
             %end;
            ) = ((0,0,&h.) (1,0,&k.));
      %end;

      /* (iii) Every pair of objects occurs in exactly lambda blocks. */
      %let t = %eval(&b.-&lambda.);
      %do i=1 %to %eval(&v.-1);
         %do j=%eval(&i.+1) %to &v;
            /* auxiliary variable p_i_j_c =1 iff both i and j occur in c */
            var ( p&i._&j._1-p&i._&j._&b. ) = [0,1];
            %do c=1 %to &b;
               reify p&i._&j._&c.: (x&i._&c. + x&j._&c. = 2);
            %end;

            gcc(p&i._&j._1-p&i._&j._&b.) = ((0,0,&t.) (1,0,&lambda.));
         %end;
      %end;

      /* Symmetry breaking constraints: */
      /* Break row symmetry via lexicographic ordering constraints. */
      %do i = 2 %to &v.;
         %let i1 = %eval(&i.-1);
         lexico( (x&i._1-x&i._&b.) LEX_LT (x&i1._1-x&i1._&b.) );
      %end;

      /* Break column symmetry via lexicographic ordering constraints. */
      %do c = 2 %to &b.;
         %let c1 = %eval(&c.-1);
         lexico( ( %do i = 1 %to &v.;
                      x&i._&c.
                   %end; )
                 LEX_LE
                 ( %do i = 1 %to &v.;
                      x&i._&c1.
                   %end; ) );
      %end;
   run;
   %put &_orclp_;
%EXIT:
%mend bibd;


%bibd(15,15,7,7,3);


%macro bibd_out(v, b, r, k, lambda, out=bibdout, transpose=0);
   /* Create a binary matrix with v rows and b columns from the solution
      of the bibd macro.  If transpose = 1 then the matrix will be transposed
      for convenience of display.
   */
   data bibdmat;
      set &out;
      array Block{&b.};
      %do i = 1 %to &v.;
         %do j = 1 %to &b.;
            Block[&j.] = x&i._&j.;
         %end;

         output;
      %end;

      drop x:;
   run;

   %if &transpose %then %do;
      /* Transposes the rows and columns of the binary matrix for
         convenience of display.  */
      proc transpose data=bibdmat
                     out=bibdmat2(rename=(_NAME_=Block))
                     prefix=Object;
      run;

      /* Print the solution */
      proc print data=bibdmat2;
         title "Balanced Incomplete Block Design Problem";
         title2 "(&v, &b, &r, &k, &lambda)";
      run;
   %end;
   %else %do;
      /* Print the solution */
      proc print data=bibdmat;
         title "Balanced Incomplete Block Design Problem";
         title2 "(&v, &b, &r, &k, &lambda)";
      run;
%end;
%mend bibd_out;

%bibd_out(15,15,7,7,3);