/* Richard A. DeVenezia
 * www.devenezia.com
 *
 * Temperature difference matrix
 *
 * On July 3, 2003 Ken Keung posed an interesting question to SAS-L

Subject: A computation question


Hello,

I have a question.
Here is a dataset with 6 variables.
Let me explain the variables first.

OBS : Observation number
C1: Component 1 (dummy variable: 1="In the pot", and 0="NOT in the pot")
C2: Component 2 (dummy variable: 1="In the pot", and 0="NOT in the pot")
C3: Component 3 (dummy variable: 1="In the pot", and 0="NOT in the pot")
HOW_MANY: How many components in the pot.
TEMP : Temperature

OBS C1 C2 C3 HOW_MANY  TEMP
1   1  0  0  1         4
2   0  1  0  1         1
3   0  0  1  1         3
4   1  1  0  2         16
5   1  0  1  2         13
6   0  1  1  2         9
7   1  1  1  3         22

As you can see, all the combinations of 3 components are listed.

If we put the Component 1 into a empty pot,
then the temperature goes to 4 from 0.

When Component 2 is added to the pot (so it has Component 1 and 2),
the temperature increases to 16.

The increment of temperature is 12 (=16 - 4), which is due to the addition
of Component 2 to the existing pot.

One thing to note is that the temperature always goes up when a component
is added.

I want to produce an output table like this.

OBS C1 C2 C3 HOW_MANY
1   4  .  .  1
2   .  1  .  1
3   .  .  3  1
4   15 12 .  2
5   10 .  9  2
6   .  6  8  2
7   13 9  6  3

For clarification, consider C1 in OBS 7, its value is 13.
13 is temperature increase (from 9 to 22), when Component 1
is added to the pot when it already contained Component 2 and 3.

-----

My challenge is that I have 20 components, NOT 3 as shown above.
That means there are 1,048,575 observations
(2 to the power of 20 minus 1) in the dataset.

After approximately 20 hours (YES! hours) waiting, my computer (1.7GHZ and
512MB memory) couldn't produce the output.

I know that I should not blame the computer.
Anyway, I was able to obtain the sas output with 15 Components case in an
hour.

Would you kindly help me to code it more efficient way ?

Thank you very much for your help in advance.
--------------------------------------------------------------------------------
 */

/*
 * Right away, the two state values, 0 and 1 of the variables C1-C20
 * smelled of binary.  Any attempt to deal with C1-C20 in a non binary sense
 * would be foolhardy.
 *
 * Ken states he has temperature data for all cases of c1-c20
 * Let x = binary value represented by in/out state values found in c1-c20
 * Thus T[x] is known for 0<x<2^20
 *
 * Kens problem is essentially finding temperature
 * differences between two elements located an array of 2^20 temperatures.
 *
 * For each T[x], 0<x<=2*20,
 * compute D[x,1]...D[x,20]
 * where D[x,i] = T[x] - T[y]
 * where y = x with the ith bit zeroed
 *
 * The following program solves the 20 component problem in about 60 seconds.
 * It takes another 8 seconds to generate realistic fake data.
 * The testbed was P4 3.06GHz, 1g RAM, 80g HD w 8m buffer, running Win2K Pro.
 */


%macro fakedata (N=);

  /*
   * apriori knowledge is that temperature always increases by a little per bit added
   * lets us create some fake data similar to what T[x] might really be
   */

  %local i j;

  data Temperatures;

    attrib c1-c&N length=3 format=1.;

    retain seed 12345;

    array c c1-c&N;
    retain c1-c&N 0;
    array bit bit1-bit&N;

    bit0 = 0;
    temperature = 0;

    oldmax = 0;
    newmax = 0;

    %do i = 1 %to &N;

      %do j = 1 %to &i;

        do bit&j = bit%eval(&j-1)+1 to &N - %eval(&i - &j);

          c[bit&j] = 1;

      %end;

      temperature = oldmax + 1 + int(15*ranuni(seed));

      newmax = max (newmax, temperature);

      output;

      %do j = 1 %to &i;

          c[bit%eval(&i-&j+1)] = 0;

        end;

      %end;

      oldmax = newmax;
    %end;

    keep c1-c&N temperature; * oldmax newmax;
  run;

%mend;

options mprint;

*/*;
%let N = 20;
%*let N = 15;
%fakedata (N=&N);
*/;


 /*
%let N = 3;
data Temperatures;
input c1 c2 c3 temperature;
cards4;
1 0 0 4
0 1 0 1
0 0 1 3
1 1 0 16
1 0 1 13
0 1 1 9
1 1 1 22
;;;;
*/;


/*
 * Compute D[x,1] ... D[x,&N]
 */

%let nStates = %sysevalf (2**&N-1);

data ContributionsToStateTemperature;

  array T [ 0:&nStates, -1:&N ] _temporary_;

  length state 8 stateR $&N. temperature 8 nbits 3;
  array c c1-c&N;

  format state binary&N..;
  format temperature c1-c&N 4. nbits 2.;

  keep state stateR temperature nbits c1-c&N;

  * store original temperatures in T[*,0];

  _n_ = 0;

  do while (not endOfData);
    set Temperatures end=endOfData;
    _n_ + 1;
    state = 0;
    do i = 0 to &N-1;
      state = bor (state, c[i+1]*2**i);
    end;
    T [ state, 0 ] = temperature;
    T [ _n_, -1 ] = state;  * inverse map for outputting in original order;
  end;

  T [ 0, 0 ] = 0;

  * loop through each state;

  do state = 0 to &nStates;

    * at each state see what state is arrived at when
    * each not in the pot component is in turn placed in the pot;

    do i = 0 to &N-1;

      * determine state that would occur if
      * some of component i added to pot;

      newstate = bor (state, 2**i);

      if newstate ne state then do;

        * added something that was not there, so
        * now we can compute D (how much the temperature will go up
        * by adding the component to the pot);

        T [ newstate, i+1 ] = T [ newstate, 0 ] - T [ state, 0 ];
      end;
    end;
  end;

  * output the delta matrix;

  do _n_ = 1 to &nStates;
    state = T [ _n_, -1 ] ;
    if state = . then continue;
    nbits = 0;
    do i = 1 to &N;
      c[i] = T [ state, i ];
      nbits + ( c[i] ne . );
    end;
    temperature = T [ state, 0 ];
    stateR = reverse(put(state,binary&N..));
    output;
  end;

  stop;
run;