/ Combinatorial design
/ The design here is unbalanced, but offers complete pairwise coverage.
/ This is also calledn a non-orthogonal array design of strength two.
/ Given a set of attributes with their values, generate all
/ possible experiments such that every combination of every pair
/ of attribute values is present in at least one experiment.
/ That goes into the table jakeout.
/ Here is a typical input format. The last line concerning pivots is optional.
/ foo 1 2 3
/ bar 1 2 3 4
/ zoo 1 2
/ dodo a b
/ zin c
/ fump c f
/ pivot: zoo 
/ May 6, 2003: added a preference for as much orthogonality as possible.
/ Notice that it is possible to design a no pivot and pivot at the
/ same time. e.g. a no pivot on six inputs gives 3 lights and 3 darks.
/ To turn this into a pivot, just convert the 3 dark ones to be identical
/ but have lights instead and similarly for the darks.



/ BASIC ROUTINES




/ returns one if x is a subset of y
subset:{[x; y]
  i: y ?/: x
  : ~ (#y) _in i
}

/ returns one if x is a subset of y
subset:{[x;y] (#y) > |/ y ?/: x}        

differ:{[x;y]  x[& ~ x _lin y]}


/finds intersection of two lists
/ fastest of all
intersect: {[x;y]
  x,: ()
  y,: ()
  if[(#x) < (#y)
    i: x ?/: y
    j: & i < #x
    :x[?i[j]]
  ]
  i: y ?/: x
  j: & i < #y
  :y[?i[j]]
}

/finds intersection of two lists
/ fastest of all
hasintersect: {[x;y]
  x,: ()
  y,: ()
  i: x ?/: y
  : (&/i) < #x
}




/finds indexes in x that intersect with y
intersectleftindexes: {[x;y]
  i: x ?/: y / where each y hits
  j: & i < #x / those ys that hit
  :i[j]
}

/finds intersection of two lists
/ and returns index pairs of matches. Assumes no duplicates
/ in either list
intersectbothindexes: {[x;y]
  x,: ()
  y,: ()
  i: x ?/: y
  pairs: (i ,' (!#y))
  k: & pairs[;0] < #x
  :pairs[k]
}

/ finds intersection of two lists that may have duplicates
bagintersectbothindexes:{[x;y]
  alreadyused: out: ()
  i: 0
  while[i < #x
	my: x[i]
	jj: & y ~\: my
	jj: differ[jj; alreadyused]
	if[0 < #jj
		out,: ,(i;*jj)
		alreadyused,: *jj
	]
	i+: 1
  ]
  :out
}




/ END BASICS


/ DUMP TABLE dumptable


/ formstring takes a list and makes a string
formstring:{[list]
  list,: ()
  : (-1) _ ,/ ($list) ,\: (" ")
}

formstringvertbar:{[list]
  list,: ()
  : (-1) _ ,/ ($list) ,\: ("|")
}

formstringcomma:{[list]
  list,: ()
  : (-1) _ ,/ ($list) ,\: (",")
}


/ Output a table (a variable) to a text file outfile (string)
/ e.g.   output[`guide; guide; "foobar"]
dumptable:{[tablename; table; outfile]
  out: ,("# "), ($tablename), ("|"), formstringvertbar[!table]
  first: *!table
  numofelements: . ("#"), ($tablename), ("."), ($first)
  i: 0
  while[i < numofelements
	list: table[;i]
	x: formstring'list
	out,: , ($x[0]), (":	"), (-1) _ ,/(1 _ x) ,\: (" ")
	i+: 1
  ]
  outfile 0: out
}

dumpmatcsv:{[mat; outfile]
	if[`disregard _in mat[0]
	  mat: (-1) _' mat
	]
	othernums: 1 _ !#mat[;0]
	i: 1
	while[i < #mat[;0]
	  mat[i]: (` $ $ i),mat[i]
	  i+: 1
	]
	mat[0]: (`num), mat[0]
  out: , formstringcomma[mat[0]]
  / first: *!table
  numofelements: #mat[;0]
  i: 1
  while[i < numofelements
	list: mat[i]
	x: formstring'list
	out,: , (-1) _ ,/x ,\: (",")
	i+: 1
  ]
  outfile 0: out
}

/ APPLICATION SPECIFIC

/ given attnames, attvalues, genpairs. Fill the exper table.
genexper:{[attpairs]
  while[0 < #attpairs
	pair: finddisjoint[attpairs]
	attpairs: pair[1]
	addexper[pair[0]]
  ]
}


/ given a set of pairs, find a disjoint set
/ we will process each such disjoint set
finddisjoint:{[x]
  / x: x[(#x) _draw -#x] / may want this in any case
  out: ()
  remaining: ()
  out,: ,*x
  current: ,/out
  x: 1 _ x
  while[0 < #x
	cands: *x
	flag: hasintersect[cands; current]
	if[flag; remaining,: ,*x]
	if[~flag
		out,: ,*x
		current,: *x
	]
	x: 1 _ x
  ]
  :(out; remaining)
}

permute:{[list]
  if[0 = #list; :list]
  :list[(#list) _draw -#list]
}

/ Given the pairs, generate the rows for the experiment by seeing
/ what is left to be done and then doing it.
/ Keep track of what can be collapsed.
addexper:{[pairs]
 origpairs: pairs
 x: differ[!#allatts; ,/pairs]
 pairs,: (_ (#x) % 2; 2) # x
 valsleftall: stilltodo'pairs / gives list of lists of values left
 maxcount: |/ #:'valsleftall
 if[0 < maxcount
  j: 0
  while[j < #valsleftall
	pair: pairs[j]
	valsleft: valsleftall[j]
	if[maxcount > #valsleft
		rest: (maxcount - #valsleft) # ,(`X;`X)
		valsleft,: rest
	]
	if[0 < #valsleft
	 vals:: valsleft[;0]
	 outmat[pair[0]],: vals
	 vals:: valsleft[;1]
	 outmat[pair[1]],: vals
	]
	j+: 1
  ]
 ]
}


/ collect attributes being compared
collect:{[indexpairs]
  out: ()
  i: 0
  while[i < #indexpairs
	pair: indexpairs[i]
	out,: ("["),($attnames[pair[0]]), (", "), ($attnames[pair[1]]), ("]")
	i+: 1
  ]
  :` $ out
}
  

/ faster routine that fails to eliminate some Xs.
stilltodo:{[pair]
  valuestodo: ,/attvalues[pair[0]] ,/:\: attvalues[pair[1]]
  valsdone: outmat[pair[0]] ,' outmat[pair[1]]
  valsleft: differ[valuestodo; ?valsdone]
  if[0 < #valsleft
        / try to compress based on pairs of values that are free
        jj: & valsdone ~\: (`X;`X)
        / change routine
        minnum: (#jj) & (#valsleft)
        if[0 < minnum
         jj@: !minnum
         valsdone[jj;0]: valsleft[!minnum; 0]
         outmat[pair[0]]:: valsdone[;0]
         valsdone[jj;1]: valsleft[!minnum; 1]
         outmat[pair[1]]:: valsdone[;1]
         valsleft: minnum _ valsleft
        ]
   ]
  :valsleft
}








/ find cross product of a list of value lists
/ It may have one element or many
crosssets:{[listofvallists]
  listofvallists,: ()
  out: listofvallists[0]
  if[1 = #listofvallists; :,:' listofvallists[0]]
  i: 1
  while[i < #listofvallists
        out: ,/ out ,/:\: listofvallists[i]
        i+: 1
  ]
  :out
}

/ INPUT GENERIC DATA


/ parses a field based on spaces
getfieldssymb :{[symbline]
  line: $symbline
  i: line = " "
  j1: &i
  j2: &~i
  line @:j2
  size: #j1
  x:` $ (0,(j1 - !size)) _ line
  ii: & ~ x = `
  :x[ii]
}

/ get rid of blanks at either end of the string
delendblanks:{[string]
  if[0 = #string; :""]
  if[string ~ ,"" ; :""]
  i: & ~ string = " "
  if[(#string) = (#i); :string]
  if[0 = (#i); :""]
  string: (- ((#string)  - (1 + *|i))) _ string
  :(*i) _ string
}
  
spit:{[line] (-2) _ ,/ ($line) ,\: (", ")}


/ Handles one line of input at a time 
/ in format 
/ attname value(s)
/ or
/ pivot: attname(s)
processline:{[line]
 emptyorblank:{[line] (0 = #delendblanks[line])}
 emptyflag: emptyorblank[line]
 if[~emptyflag
  myline: getfieldssymb[delendblanks[line]]
  pivotflag: ($myline[0]) _sm "*pivot*"
  if[~pivotflag
	attnames,: myline[0]
	attvalues,: ,1 _ myline
  ]
  if[pivotflag
	pivot,: 1 _ myline
	if[~ subset[pivot; attnames]
		x: (" The pivots ("), (spit[pivot]),(") are not a subset of ")
		x,: (" the attributes ("), (spit[attnames]), (").")
		` 0: ,x
		. "\\\\"
	]
  ]
 ]
}

even:{[num] (num % 2) = _ (num % 2)}


/ origmat has format att and values
/ the pivot attributes have willbepivot and X throughout.
/ These must be changed to the values in the cross-product.
/ attindexes are those indexes of attributes to be given the pivotvals.
insertpivots:{[origmat; pivotvals; attindexes]
  firsts: *:' origmat
  out: ,:' firsts
  otherindexes: differ[(!#origmat); attindexes] / non-pivots
  startmat: 1 _' origmat / withoutlabels
  mycross: crosssets[pivotvals] 
  i: 0
  while[i < #mycross
	thiscross: mycross[i]
	x: startmat
	j: 0
	while[j < #attindexes
		k: attindexes[j]
		out[k],: ((#x[k]) ) # thiscross[j]
		j+: 1
	]
	j: 0
 	while[j < #otherindexes
		k: otherindexes[j]
		out[k],: x[k]
		j+: 1
	]
	i+: 1
  ]
  :out
}
  

/ EXECUTION

attnames: ()
attvalues: () / will be list of lists
pivot: ()
file: _i[0] / should be a file with attributes and possibly a pivot

x: processline' 0: file

/ Number of attnames must be even. We add a dummy called disregard if needed. 
if[~ even[#attnames]
  attnames,: `disregard
  attvalues,: ,,`disregardval
]

origattvalues: attvalues
if[0 < #pivot
  ii: intersectleftindexes[attnames; pivot]
  pivotvals: (attvalues[ii]) ,()
  j: 0
  while[j < #ii
  	attvalues[ii[j]]: ,`willbepivot
	j+: 1
  ]
]

allatts: attnames



/ EXECUTION
x: _gtime _t
y: ((x[1]) ! 1003) _draw 300


if[0 < #pivot
  ii: intersectleftindexes[allatts; pivot]
  j: 0
  while[j < #ii
	allatts[ii[j]]: ` $ ("*"), ($allatts[ii[j]]), ("*")
	j+: 1
  ]
]
outmat:  ,:'(allatts)  / then we will fill in line as needed
	/ at the end we will do dumpmatcsv to get the csv file
bestmat: outmat
bestcountexper: 999999999
vals: () / used as a global
indexes: () / used as a global

genpairs: ,/(!#attnames) ,/:\: (!#attnames)


genpairs@: & genpairs[;0] < genpairs[;1]

/ ratio of fewest in a line with most
hilow:{[line]
  part: = line
  counts: #:' part
  : (&/counts) % (|/counts)
}
firsttimeflag: 1
firsttimeflag: 0
genexper[genpairs]
if[0 < |/#:'stilltodo'genpairs; !-1] / some genpairs are missing
countexper: #outmat[0]
if[countexper < bestcountexper
  bestcountexper: countexper / the number to beat
  outmatbest: outmat
  besthilow: +/hilow' 1 _' outmat
]
firsttimeflag: 0 / no randomness above
numiter: 1
if[1 < #_i; numiter: 0 $ _i[1]]
do[numiter
  outmat:  ,:'(allatts)  / then we will fill in line as needed
  genexper[genpairs]
  countexper: #outmat[0]
  myhilow: +/hilow' 1 _' outmat
  if[(~ countexper > bestcountexper)
	/ as good by number of experiments criterion
	/ now check for distribution of values
	if[myhilow > besthilow
          bestcountexper: countexper
  	  outmatbest: outmat
	  besthilow: myhilow
	]
  ]
]

outmat: outmatbest
if[0 < |/#:'stilltodo'genpairs; !-1] / some genpairs are missing

dumpmatcsv[(+outmat); "jakeout"]



\\