Re: [Caml-list] OCaml troll on Slashdot

[YANG Shouxun <yang.shx@fltrp.com>]

On Wednesday 16 March 2005 04:02, Yoann Padioleau wrote:
> YANG Shouxun <yang.shx@fltrp.com> writes:
> > No. My experiments show that the OCaml implementation performs far better
> > than the C++ implementation when the column and row get larger:
>
> Perhaps because you are a liar.

I'm 100% sure I'm not a liar. I don't know what's the benefit of being a liar, 
nor the benefits of calling others a liar.

I checked and get the same results, though the figures are of course not 
exactly the same.

> > C++ is compiled with -O3, not sure what is the better optimization level,
> > while OCaml version (actually I used Sys.argv and references to the two
> > parameters) is compiled with ocamlopt
> >
> > 4x4 c++
> > real    0m0.003s
> > user    0m0.002s
> > sys     0m0.002s
> >
> > 4x4 ocaml
> > real    0m0.177s
> > user    0m0.137s
> > sys     0m0.001s
> >
> > 8x8 c++
> > real    0m8.703s
> > user    0m7.000s
> > sys     0m0.018s
> >
> > 8x8 ocaml
> > real    0m0.747s
> > user    0m0.485s
> > sys     0m0.002s
>
> I dont know where you get those numbers ?

As I said, the C++ version was compiled with "-O3" and the OCaml version with 
ocamlopt. The compiled programs were timed using "time garden_XXX x y" where 
XXX is either cpp for C++ version and ml for OCaml version.

I'm using Debian GNU/Linux, with g++ version 3.3.5 and ocamlopt version 
3.08.2. The computer is an Intel P4 2.40GHz with 256MB memory and 512KB cache 
and 506MB swap.

Below is what I got from the site of C++ version:
---8<---

#include <math.h>
#include <stdio.h>


enum cell{ Empty=0, Planted=1 }; 


int columnSize = 5;

int numStates;
cell **allStates;

int *allCosts;

// array indicating when a center column of cells is covered by
// the columns on either side of it.
// indexed as isCovered[c1][c2][c3]   true if c1 and c3 cover c2,
// or false otherwise
char ***isCoveredArray;


char isCovered( int inC1, int inC2, int inC3 ) {
    cell *c1 = allStates[ inC1 ];
    cell *c2 = allStates[ inC2 ];
    cell *c3 = allStates[ inC3 ];

    for( int i=0; i<columnSize; i++ ) {
        // test for a cell that is not covered and return false
        if( c2[i] == Planted &&
            c1[i] == Planted &&
            c3[i] == Planted ) {

            if( i <=0 || c2[i-1] == Planted ) {
                if( i >=columnSize-1 || c2[i+1] == Planted ) {

                    return false;
                    }
                }
            }
        }

    // else all covered
    return true;
    }


void printState( int inStateIndex ) {
    cell *state = allStates[ inStateIndex ];

    for( int i=0; i<columnSize; i++ ) {
        if( state[i] == Empty ) {
            printf( "O" );
            }
        else {
            printf( "X" );
            }
        }
    printf( "\n" );
    }


struct costStatePair {
        int cost;
        // the index in the allStates array
        int stateIndex;
    };


// table for memoizing optLayout
struct costStatePair ***optLayoutTable;



// place where optLayout results are returned
struct costStatePair optLayoutResult;

/**
 * Computes the optimum layout for a given number of columns given
 * states for the first two columns.
 *
 * @param inNumColumns the number of columns.  Must be at least 3.
 * @param inFirstColumnStateIndex the index of the state (in allStates)
 *   of the first column.
 * @param inSecondColumnStateIndex the index of the state (in allStates)
 *   of the second column.
 *
 * @return result is set in the optLayoutResult structure and returned.
 *   Thus, this call is NOT threadsafe.   
 */
struct costStatePair optLayout( int inNumColumns, int inFirstColumnStateIndex,
                                int inSecondColumnStateIndex ) {

    int c = inNumColumns-1;
    if( optLayoutTable[c]
                      [inFirstColumnStateIndex]
                      [inSecondColumnStateIndex].cost != -1 ) {
        // already have this result
        return optLayoutTable[c]
                             [inFirstColumnStateIndex]
                             [inSecondColumnStateIndex];
        }
    
    if( inNumColumns < 3 ) {
        optLayoutResult.cost =
            allCosts[ inFirstColumnStateIndex ] +
            allCosts[ inSecondColumnStateIndex ];
        optLayoutResult.stateIndex = -1;


        // save in table
        optLayoutTable[c]
            [inFirstColumnStateIndex]
            [inSecondColumnStateIndex].cost = optLayoutResult.cost;
        optLayoutTable[c]
            [inFirstColumnStateIndex]
            [inSecondColumnStateIndex].stateIndex = 
optLayoutResult.stateIndex;
        
        
        return optLayoutResult;
        }
            
    
    
    int bestCost = inNumColumns * columnSize + 1;
    int bestStateIndex = -1;

    // examin all possible states for the third column
    for( int i=0; i<numStates; i++ ) {

        // if i does its part to cover our second column
        if( isCovered( inFirstColumnStateIndex,
                       inSecondColumnStateIndex, i ) ) {

            char iCovered = false;
            
            int costOfRest = 0;
            if( inNumColumns > 3 ) {
                struct costStatePair result =
                    optLayout( inNumColumns - 1,
                               inSecondColumnStateIndex,
                               i );
                costOfRest = result.cost -
                    allCosts[ inSecondColumnStateIndex ];

                // optLayout will only pick a 3rd column that helps to
                // cover i
                iCovered = true;
                }
            else {
                // just the cost of i by itself
                costOfRest = allCosts[i];

                // make sure that i is covered by inSecondColumnStateIndex
                // put the all-planted column (0) on the other side of i
                iCovered = isCovered( inSecondColumnStateIndex, i, 0 );
                }

            if( iCovered && costOfRest < bestCost ) {
                bestCost = costOfRest;
                bestStateIndex = i;
                }
            }
        }
    
    optLayoutResult.cost = bestCost +
        allCosts[ inFirstColumnStateIndex ] +
        allCosts[ inSecondColumnStateIndex ];
    
    optLayoutResult.stateIndex = bestStateIndex;


    // save in table
    optLayoutTable[c]
        [inFirstColumnStateIndex]
        [inSecondColumnStateIndex].cost = optLayoutResult.cost;
    optLayoutTable[c]
        [inFirstColumnStateIndex]
        [inSecondColumnStateIndex].stateIndex = optLayoutResult.stateIndex;

    
    return optLayoutResult;
    }




void printOptLayout( int inNumColumns ) {
    // find best starting state for 10 columns by
    // running optLayout on 11 columns
    int bestFirstColumn = -1;
    int bestSecondColumn = -1;
    int bestCost = (inNumColumns+1) * columnSize + 1;

    
    
    for( int i=0; i<numStates; i++ ) {
        // state 0 is the all-planted state
        struct costStatePair result = optLayout( inNumColumns+1, 0, i );

        if( result.cost < bestCost ) {
            bestFirstColumn = i;
            bestSecondColumn = result.stateIndex;
            bestCost = result.cost;            
            }
        //printf( "best cost = %d\n", bestCost );
        }

    printf( "%d empty cells:\n", bestCost );
    printState( bestFirstColumn );
    printState( bestSecondColumn );
    
    for( int c=inNumColumns; c>=3; c-- ) {

        struct costStatePair restult =
            optLayout( c, bestFirstColumn, bestSecondColumn );
        
        printState( restult.stateIndex );

        bestFirstColumn = bestSecondColumn;
        bestSecondColumn = restult.stateIndex;
        }
    }






int main( int inNumArgs, char **inArgs ) {

    int numColumns = 10;

    if( inNumArgs > 1 ) {
        sscanf( inArgs[1], "%d", &columnSize );
        }
    if( inNumArgs > 2 ) {
        sscanf( inArgs[2], "%d", &numColumns );
        }
    
    
    numStates = (int)( pow( 2, columnSize ) );
    allStates = new cell*[ numStates ];
    allCosts = new int[ numStates ];


    
    
    
    int i, j, k;
    for( i=0; i<numStates; i++ ) {

        allCosts[i] = 0;

        allStates[i] = new cell[ columnSize ];
        for( int b=0; b<columnSize; b++ ) {

            if( ( ( i >> b ) & 0x1 ) == 1 ) {
                allStates[i][b] = Empty;
                allCosts[i] ++;
                }
            else {
                allStates[i][b] = Planted;
                }
            }
        }

    int numEntries = 0;
    optLayoutTable = new struct costStatePair **[numColumns +1];
    int c;
    for( c=0; c<numColumns+1; c++ ) {
        optLayoutTable[c] = new struct costStatePair *[numStates];
        for( int i=0; i<numStates; i++ ) {
            optLayoutTable[c][i] = new struct costStatePair[ numStates ];

            for( int j=0; j<numStates; j++ ) {
                optLayoutTable[c][i][j].cost = -1;
                optLayoutTable[c][i][j].stateIndex = -1;

                numEntries ++;
                }
            }            
        }
    printf( "Done constructing memoization table with %d entries\n",
            numEntries );
    
    /*
    isCoveredArray = new char**[ numStates ];
    for( i=0; i<numStates; i++ ) {
        isCoveredArray[i] = new char*[numStates];
        for( j=0; j<numStates; j++ ) {
            isCoveredArray[i][j] = new char[numStates];
            for( k=0; k<numStates; k++ ) {
                isCoveredArray[i][j][k] = isCovered( i, j, k );
                }
            }
        }
    */

    /*
    for( i=0; i<numStates; i++ ) {
        for( j=0; j<numStates; j++ ) {
            delete [] isCoveredArray[i][j];
            }
        delete [] isCoveredArray[i];
        }
    delete [] isCoveredArray;
    */



    for( c=2; c<=numColumns; c++ ) {
        printf( "\n%d by %d\n", columnSize, c );

        printOptLayout( c );
        }

    
    


    for( c=0; c<numColumns+1; c++ ) {
        
        for( int i=0; i<numStates; i++ ) {
            delete [] optLayoutTable[c][i];
            }
        delete [] optLayoutTable[c];    
        }
    delete [] optLayoutTable;
    
    for( i=0; i<numStates; i++ ) {
        delete [] allStates[i];
        }
    delete [] allStates;
    delete [] allCosts;
    
    return 0;
    }

---8<---

The original OCaml version:

---8<---
(* Computes optimal garden layouts using dynamic programming*)

(* the width of a state in our dynamic programming table *)
(* the algorithm will run in O( 2^columnSize ) time *)
let columnSize = 4;;

(* How many rows to compute solutions for *)
let numRows = 10;;


(* generate all possible states for columns of height h *)
type cell = Empty 
       | Planted;;
    

let printCell c =
  match c with
 | Empty   -> "O"
 | Planted -> "X";;


let rec printCells cs = 
  match cs with
 | []        -> ""
 | (c::rest) -> (printCell c) ^ (printCells rest);;


let rec printStates states =
  match states with
 | []        -> ""
 | (s::[])   -> printCells s
 | (s::rest) -> (printCells s) ^ " " ^ (printStates rest);;



(* apply f to every element of list *)
let rec map f list =
  match list with
 | []      -> []
 | x::rest -> (f x) :: (map f rest);;



let rec generateStates numStates =
  match numStates with 
 | 0 -> [[]]
 | h ->
  let 
   shorterStates = generateStates (h-1) and
   prepend c s = c::s
  in
    ( map ( prepend Empty ) shorterStates ) @
    ( map ( prepend Planted ) shorterStates );;



let allStates = generateStates columnSize;;

let rec genIndexList currentIndex maxIndex =
  if( currentIndex <= maxIndex ) 
  then currentIndex :: ( genIndexList (currentIndex+1) maxIndex )
  else [];;


let rec length list =
  match list with
 | []      -> 0
 | x::rest -> 1 + length rest;;


let stateIndices = genIndexList 0 ((length allStates) - 1);;


let rec findIndexRec s states startIndex =
 match states with
   | []      -> -1
   | (x::xs) ->
    if( s = x ) then startIndex
    else findIndexRec s xs (startIndex+1);;

let findIndex s = findIndexRec s allStates 0;;
  
 


let rec generateSolidState length cellValue =
  match length with 
 | 0 -> []
 | x -> cellValue::(generateSolidState (x-1) cellValue );;

let allPlantedState = generateSolidState columnSize Planted;;
let allEmptyState = generateSolidState columnSize Empty;;


(* computes the min cost pair in a list of cost-state pairs *)
let rec minCost pairs =
  match pairs with
 | []          -> ( 100000, allEmptyState )  (* error *)
 | (c,s)::[]   -> (c,s)
 | (c,s)::rest ->
  let
   ( minCostRest, minStateRest ) = minCost rest
  in
    if( c < minCostRest ) then (c,s)
    else ( minCostRest, minStateRest );;

(* computes the min cost pair in a list of cost-state-state triples *)
let rec minCost3 triples =
  match triples with
 | []          -> ( 100000, allEmptyState, allEmptyState )  (* error *)
 | (c,s1,s2)::[]   -> (c,s1,s2)
 | (c,s1,s2)::rest ->
  let
   ( minCostRest, minState1Rest, minState2Rest ) = minCost3 rest
  in
    if( c < minCostRest ) then (c,s1,s2)
    else ( minCostRest, minState1Rest, minState2Rest );;


(* adds c to the cost of the minimum-cost pair in a list *)
let addToMin c pairList = 
  let (minC, minS) = minCost pairList in
 (c + minC, minS );;
 


(* computes the cost of a state*)
let rec cost state =
  match state with
 | []  -> 0
 | Empty::rest -> 1 + cost rest
 | Planted::rest -> cost rest;;
    
let costList = map cost allStates;;


(* Set up an associative list for memoization *)
let lookup key table = List.assoc key !table;;
let insert key value table = table := (key, value) :: !table;;


(* memoize any 3-parameter function *)
let memoize3 table f x y z =
    try
        lookup (x,y,z) table
    with
    | Not_found ->
        let result = f x y z in
        insert (x, y, z) result table;
        result;;

(* table for memoizing optLayout *)
let isCovered_table = ref [];;

(* checks if each cell in center colum is covered by an empty cell *)
let rec isCovered c1 c2 c3 =
  let 
   memoized_isCovered = memoize3 isCovered_table isCovered 
  in 
  match (c1, c2, c3) with
   (* if center runs out of cells first, we are covered *)
 | (_, [], _) -> true
   (* if others run out first, we are not covered *)
 | ([], _, _) -> false
 | (_, _, []) -> false
   (* Empty followed by Planted in center colum
      skip this cell and next cell *)
 | ( (_::(_::rest1)), (Empty::(Planted::rest2)), (_::(_::rest3)) ) -> 
  true && ( isCovered rest1 rest2 rest3 )
 
 | ( (Empty::rest1), (_::rest2), (_::rest3) ) -> 
  true && ( isCovered rest1 rest2 rest3 )

 | ( (_::rest1), (_::rest2), (Empty::rest3) ) -> 
    true && ( isCovered rest1 rest2 rest3 )

   (* Empty followed by Empty in center
      This cell is covered, but don't skip next cell *)
 | ( (_::rest1), (Empty::rest2), (_::rest3) ) -> 
  true && ( isCovered rest1 rest2 rest3 )
   (* this cell is covered by next cell *)
 | ( (_::rest1), (Planted::(Empty::rest2)), (_::rest3) ) -> 
  true && ( isCovered rest1 (Empty::rest2) rest3 )

   (* default:  we reach this if our current cell is not covered *)
 | (_, _, _) -> false;;


let memoized_isCovered = memoize3 isCovered_table isCovered;;

(*

-- this is a lazy array of arrays of arrays.
-- first index is number of columns (-1... so 1 colum result is
--   at index 0
-- Second index is the state index (an index into allStates).  This
--    is the state of the first column
-- Third index is the state index (an index into allStates).  This
--    is the state of the column before the first column
-- So, first dimension is infinite, while second and third dimensions are
-- finite.
optSolutionArray = [  list2 | 
       f <- (map optLayout [1..]), 
       list <- [ map f stateIndices ], 
                      list2 <- [
        [ map f2 stateIndices | f2 <- list ] 
          ]
       ]

-- lazy array memoizing isCovered
isCoveredArray = [  list2 | 
     f <- (map isCovered allStates), 
     list <- [ map f allStates ], 
                    list2 <- [
         [ map f2 allStates | f2 <- list ] 
        ]
     ]
*)

(* applies a list of functions to an argument *)
let rec applyFunctionList functions argument =
  match functions with
 | [] -> []
 | f::rest -> (f argument)::(applyFunctionList rest argument);; 


exception Length_mismatch__makeTriples;;

let rec makeTriples listOfPairs listOfSingles =
  match (listOfPairs, listOfSingles) with
 | ([],[]) -> []
 | ((p1,p2)::pRest, s::sRest) -> (p1,p2,s)::(makeTriples pRest sRest)
 | _ -> raise Length_mismatch__makeTriples;;  


(* discard middle element of each triple *)
let rec discardMiddle listOfTriples =
  match listOfTriples with
 | [] -> []
 | (t1,_,t3)::tRest -> (t1,t3)::(discardMiddle tRest);;


let rec filterTriples f list =
  match list with
 | [] -> []
 | (x,y,z)::rest -> 
  if( f x y z ) then (x,y,z)::( filterTriples f rest )
  else filterTriples f rest;;

let rec filter f list =
  match list with
 | [] -> []
 | x::rest -> 
  if( f x ) then x::( filter f rest )
  else filter f rest;;



(* Simple memo fib *)
(*
let rec fib = function
    | 0 -> 0
    | 1 -> 1
    | n -> memo_fib (n-1) + memo_fib (n-2)
and memo_fib n = memoize fib n;;
*)

(* table for memoizing optLayout *)
let optLayout_table = ref [];;


(*
-- computes the optimum layout of i columns given that the first column
-- has state s and the column before the first column has state s1 
-- Ensures that first colum (s) is covered by s1 and the optimal solution
--   to the remaining colums.  Does not ensure that s1 is covered (since
--   s is fixed, this function has no control over the coverage of s1). *)
let rec optLayout numColumns s s1 =
  (*print_string "Working on optLayout ";
  print_int numColumns;
  print_string (" s=" ^ (printCells s));
  print_string (" s1=" ^ (printCells s1));
  print_newline ();
  *)
  match numColumns with
 | 1 -> ( cost s, allPlantedState )
 | i ->
  let 
   memoized_optLayout = memoize3 optLayout_table optLayout
  in
  (* so much cleaner with list comprehensions
    addToMin ( cost s )
     [ (c,s2) |
     s2 <- allStates, 
     (c,s3) <- [ optLayout (i-1) s2 s ],
     isCovered s s2 s3,
     isCovered s1 s s2 ]
  *)
  let 
   removeUncoveredS2 s2 = (isCovered s1 s s2)
  in
  let
   allGoodStates = List.filter removeUncoveredS2 allStates 
  in
  let
   (* (opt i-1) applied to all states, a list of
      functions waiting to be applied to s *)
   optAppliedToAllStates = 
           List.map (memoized_optLayout (i-1)) allGoodStates
  in
  let
   (* these are all optimal pairs for all choices of s2 *)
   optAppliedToAllS2 = applyFunctionList optAppliedToAllStates s
  in
  let
   (* triples of (c,s3,s2) *)
   optC_S3_S2 = makeTriples optAppliedToAllS2 allGoodStates
  in
  let
   (* a filter function to remove uncovered triples from our list *)
   removeUncovered (c, s3, s2) = (isCovered s s2 s3)  (*&& 
                                (isCovered s1 s s2)*)
  in
  let
   covered_optC_S3_S2 = List.filter removeUncovered optC_S3_S2
  in

  let
   (minC, minS3, minS2) = minCost3 covered_optC_S3_S2
  in
    ( (cost s) + minC, minS2 )
  
(*
-- look up result in optSolutionArray instead of computing
-- This memoization dramatically speeds up computation
-- And hey... it was very elegant to express this in Haskell
memoizedOptLayout 1 s s1 = optLayout 1 s s1
memoizedOptLayout i s s1 = optSolutionArray !! (i-1) !! s !! s1 
*)

let memoized_optLayout = memoize3 optLayout_table optLayout;;


(* computes the first column of the optimum x-column layout *)
let computeOpt x =
    (* force an extra planted column before the first column
       force an extra empty colum before that *)
  optLayout (x + 1) allPlantedState allEmptyState;;


let rec printOptGivenState numColumns s s1 = 
  match numColumns with
 | 1 -> printCells s
 | x ->
  let 
   (c, optNextState) = optLayout x s s1
  in
    ( printCells s ) ^ "\n" ^ 
   printOptGivenState (x - 1) optNextState s;;
  

(* hey... it actually works!
 this is the main function *)
let printOpt x =
  let
   (c, optStartingState) = computeOpt x
  in
 print_int c;
 print_string 
   (" empty cells:\n" ^ 
   printOptGivenState x optStartingState allPlantedState );;

let main () = 
  (* first, fill up the memo table *)
(*
  for i=1 to 10 do
 print_string( "Filling table row " );
 print_int( i );
 print_newline();
 List.map 
   (applyFunctionList ( List.map (memoized_optLayout i) allStates )) 
   allStates;
    print_string( "   done with table row " );
 print_int( i );
 print_newline();
  done;
*)
  for i = 2 to numRows do
 print_int( columnSize );
 print_string( " by " );
 print_int( i );
 print_newline();
 printOpt i;
 print_newline();
 print_newline();
  done;;

main ();;

---8<---

Since I don't want to change the source code and recompile each time I want to 
pass a different argument list, I made the following changes, without 
handling some possible exceptions:

---8<---
--- Garden.ml   2005/03/15 09:38:16     1.1
+++ Garden.ml   2005/03/16 01:22:08
@@ -2,10 +2,10 @@

 (* the width of a state in our dynamic programming table *)
 (* the algorithm will run in O( 2^columnSize ) time *)
-let columnSize = 4;;
+let columnSize = ref 4;;

 (* How many rows to compute solutions for *)
-let numRows = 10;;
+let numRows = ref 10;;


 (* generate all possible states for columns of height h *)
@@ -54,7 +54,7 @@



-let allStates = generateStates columnSize;;
+let allStates = generateStates !columnSize;;

 let rec genIndexList currentIndex maxIndex =
   if( currentIndex <= maxIndex )
@@ -88,8 +88,8 @@
        | 0 -> []
        | x -> cellValue::(generateSolidState (x-1) cellValue );;

-let allPlantedState = generateSolidState columnSize Planted;;
-let allEmptyState = generateSolidState columnSize Empty;;
+let allPlantedState = generateSolidState !columnSize Planted;;
+let allEmptyState = generateSolidState !columnSize Empty;;


 (* computes the min cost pair in a list of cost-state pairs *)
@@ -372,7 +372,7 @@

 let main () =
   (* first, fill up the memo table *)
-(*
+  (*
   for i=1 to 10 do
        print_string( "Filling table row " );
        print_int( i );
@@ -385,8 +385,11 @@
        print_newline();
   done;
 *)
-  for i = 2 to numRows do
-       print_int( columnSize );
+  let len = Array.length Sys.argv in
+    if len > 1 then columnSize := int_of_string (Sys.argv.(1));
+    if len > 2 then numRows := int_of_string (Sys.argv.(2));
+  for i = 2 to !numRows do
+       print_int( !columnSize );
        print_string( " by " );
        print_int( i );
        print_newline();
---8<---

> I tried the code of the "troll" and the ocaml version performs far _worse_
> than the c++ implementation when the column and row get larger.
>
> > 12x12 c++
> > the process was killed by the OS
> >
> > 12x12 ocaml
> > real    0m1.210s
> > user    0m0.886s
> > sys     0m0.001s

I admit I didn't look into the C++ or OCaml code yet, but that's not what I'm 
trying to say here.

Everybody is welcome to point out where I might be mistaken in producing the 
incredible results (at least to Yoann Padioleau) and confirm whether the 
experiments can be repeated or not.