/*
Include file: markov.inc
Macro's for Markov Chains of first and nth order transition matrices.

Author : Ingo Janssen
 
Version: 0.1 Alpha, 2021-06-28
         initial
*/

#version 3.8;
#ifndef(MARKOV_INC_TEMP)
  #declare MARKOV_INC_TEMP = version;
  #ifdef(View_POV_Include_Stack)
    #debug "including markov.inc\n"
  #end


/*Enumerate elements of Markov data structure*/
#declare _mc_order = 0;
#declare _mc_state = 1;
#declare _mc_tm = 2;
#declare _mc_mult = 3;
#declare _mc_stream = 4;


#macro MCMarkov(TM, Order, InitialHistory, Stream)
  /*: Create a Markov Chain structure. It keeps its own state.
  
  parameter:
    TM (array[n^order]{array[n],..}): Transition matrix containing probabilities for each possible transition.
    Order (int, [1,...]): The order of the Markov Chain.
    InitialHistory (int | arr[order]{int, ..}: Initial state of the chain. When order is 1
      an int is required. For higher orders matrices an array of int's with lengt order is required.
  
  result: 
    A Markov Chain data structure.
  */
  #if (Order < 1)
    #error "Order has to be 1 or bigger integer"
  #elseif (Order = 1)
    #local MC = _MCMarkovFirst(TM, InitialHistory, Stream)
  #else
    #local MC = _MCMarkovNth(TM, Order, InitialHistory, Stream)
  #end
  MC
#end


#macro _MCMarkovFirst(TM, InitialHistory, Stream)
  /*: Called by MCMarkov. Markov Chain data structure for first order matrix*/
  #local MC = array mixed[5];
  #local MC[_mc_order] = 1;
  #local MC[_mc_state] = InitialHistory;
  #local MC[_mc_tm] = TM;
  #local MC[_mc_stream] = Stream;
  MC
#end


#macro _MCMarkovNth(TM, Order, InitialHistory, Stream)
  /*: Called by MCMarkov. Markov Chain data structure for higher than 1 order matrix*/
  #local N = dimension_size(TM[0], 1);
  #local MC = array mixed[5];
  #local MC[_mc_order] = Order;
  #local MC[_mc_state] = InitialHistory;
  #local MC[_mc_tm] = TM;
  #local MC[_mc_mult] = array[N];
  #local MC[_mc_mult][0] = pow(N, Order - 1);
  #for (I, 1, Order-1) 
    #local MC[_mc_mult][I] = MC[_mc_mult][I-1]/N;   
  #end
  #local MC[_mc_stream] = Stream;
  MC
#end


#macro MCnextState(MC)
  /*:Progress the Markov Chain to its next state
    
    parameter:
      MC: a Markov Chain data structure
      Stream: a random number stream
    
    result:
      a new state
  */
  #if(MC[_mc_order] = 1)
    #local Result = _MCstateFirst(MC);
  #else
    #local Result = _MCstateNth(MC);
  #end
  Result
#end
 
 
#macro _MCstateFirst(MC)
  /*: Progress the first order Markov Chain to its next state*/
  #local Selector = rand(MC[_mc_stream]);
  #local RowIndex = MC[_mc_state];
  #local Row = MC[_mc_tm][RowIndex];
  #local ProbSum = 0;
  #for (I, 0, dimension_size(Row, 1)-1)
    #local ProbSum = ProbSum + Row[I];
    #if (ProbSum > Selector)
      #local State = I;
      #break  
    #end
  #end
  #local MC[_mc_state] = State;
  State
#end


#macro _MCstateNth(MC)
  /*: Progress the nth order Markov Chain to its next state*/
  #local Selector = rand(MC[_mc_stream]);
  #local RowIndex = 0;
  #for(I, 0, dimension_size(MC[_mc_state],1)-1)
    #local RowIndex = RowIndex + (MC[_mc_state][I] * MC[_mc_mult][I]);
  #end
  #local Row = MC[_mc_tm][RowIndex];
  #local ProbSum = 0;
  #local N = dimension_size(Row, 1);
  #for (I, 0, N-1)
    #local ProbSum = ProbSum + Row[I];
    #if (ProbSum > Selector)
      #local State = I;
      #break  
    #end
  #end
  #for(I, 0, dimension_size(MC[_mc_state],1)-1)
    #if(I = dimension_size(MC[_mc_state],1)-1)
      #local MC[_mc_state][I] = State;
    #else
      #local MC[_mc_state][I] = MC[_mc_state][I+1];        
    #end 
  #end 
  State
#end


#macro MCgetState(MC)
  /*: Get the current state of the Markov Chain*/
  #if(MC[_mc_order] = 1)
    #local State = MC[_mc_state];
  #else
    #local State = MC[dimension_size(MC[_mc_state],1)-1)];
  #end
  State
#end


#macro MCgetHist(MC)
  /*: get the current history of an nth order chain.
      A first order chain has no history
  */
  #if(MC[_mc_order] = 1)
    #error "A first order chain has no history, only a current state"
  #end
  #local RowIndex = 0;
  #for(I, 0, dimension_size(MC[_mc_state],1)-1)
    #local RowIndex = RowIndex + (MC[_mc_state][I] * MC[_mc_mult][I]);
  #end
  RowIndex  
#end


#if(input_file_name="markov.inc")

/*test and demo scene:
  The "probability blocks" in the transition matrices are intentionally kept
  identical. The random stream is reset for each order. The result should be
  that in all three orders the resulting new state states are identical. 
*/

#global_settings{ assumed_gamma 1.0 }
#default{ finish{ ambient 0.1 diffuse 0.9 }} 

#declare State = array[3]{"A", "B", "C"}; 
#declare Runs = 25;

// First order Markov Chain
#debug "\nFirst order Markov Chain:\n"
#declare TransitionMatrix_1 = array[3]{
  array[3]{0.4, 0.3, 0.3},  //AA
  array[3]{0.6, 0.1, 0.3},  //AB
  array[3]{0.1, 0.7, 0.2},  //AC
};
#declare Order = 1;
#declare InitialState = 0;
#declare Stream = seed(8);
#declare Mc_1 = MCMarkov(TransitionMatrix_1, Order, InitialState, Stream);
#for (I, 0, Runs-1) 
  #debug concat("current state : ", State[MCgetState(Mc_1)], "\n")
  #declare NewState = MCnextState(Mc_1);
  #debug concat("new State : ", State[NewState], "\n\n")
#end


// Second order Markov Chain
#debug "\nSecond order Markov Chain:\n"
#declare Hist_2 = array[9]{
  "AA", "AB", "AC", 
  "BA", "BB", "BC", 
  "CA", "CB", "CC"
};
#declare TransitionMatrix_2 = array[9]{
  //        A    B    C
  array[3]{0.4, 0.3, 0.3},  //AA
  array[3]{0.6, 0.1, 0.3},  //AB
  array[3]{0.1, 0.7, 0.2},  //AC
  
  array[3]{0.4, 0.3, 0.3},  //BA
  array[3]{0.6, 0.1, 0.3},  //BB
  array[3]{0.1, 0.7, 0.2},  //BC
  
  array[3]{0.4, 0.3, 0.3},  //CA
  array[3]{0.6, 0.1, 0.3},  //CB
  array[3]{0.1, 0.7, 0.2},  //CC
};
#declare Order = 2;
#declare History = array[2]{0,0};
#declare Stream = seed(8);
#declare Mc_2 = MCMarkov(TransitionMatrix_2, Order, History, Stream);
#for (I, 0, Runs-1) 
  #debug concat("current history : ", Hist_2[MCgetHist(Mc_2)], "\n")
  #declare NewState = MCnextState(Mc_2);
  #debug concat("new State : ", State[NewState], "\n")
  #debug concat("new history : ", Hist_2[MCgetHist(Mc_2)], "\n\n")
#end


// Second order Markov Chain
#debug "\nThird order Markov Chain:\n"
#declare Hist_3 = array[27]{
  "AAA", "AAB", "AAC",
  "ABA", "ABB", "ABC",
  "ACA", "ACB", "ACC",
  "BAA", "BAB", "BAC",
  "BBA", "BBB", "BBC",
  "BCA", "BCB", "BCC",
  "CAA", "CAB", "CAC",
  "CBA", "CBB", "CBC",
  "CCA", "CCB", "CCC"
};
#declare TransitionMatrix_3 = array[27]{
  //        A    B    C
  array[3]{0.4, 0.3, 0.3},  //AAA
  array[3]{0.6, 0.1, 0.3},  //AAB
  array[3]{0.1, 0.7, 0.2},  //AAC
                            
  array[3]{0.4, 0.3, 0.3},  //ABA 
  array[3]{0.6, 0.1, 0.3},  //ABB 
  array[3]{0.1, 0.7, 0.2},  //ABC 
                            
  array[3]{0.4, 0.3, 0.3},  //ACA   
  array[3]{0.6, 0.1, 0.3},  //ACB   
  array[3]{0.1, 0.7, 0.2},  //ACC
     
  array[3]{0.4, 0.3, 0.3},  //BAA   
  array[3]{0.6, 0.1, 0.3},  //BAB   
  array[3]{0.1, 0.7, 0.2},  //BAC
                            
  array[3]{0.4, 0.3, 0.3},  //BBA   
  array[3]{0.6, 0.1, 0.3},  //BBB
  array[3]{0.1, 0.7, 0.2},  //BBC

  array[3]{0.4, 0.3, 0.3},  //BCA
  array[3]{0.6, 0.1, 0.3},  //BCB                         
  array[3]{0.1, 0.7, 0.2},  //BCC
  
  array[3]{0.4, 0.3, 0.3},  //CAA
  array[3]{0.6, 0.1, 0.3},  //CAB
  array[3]{0.1, 0.7, 0.2},  //CAC
  
  array[3]{0.4, 0.3, 0.3},  //CBA
  array[3]{0.6, 0.1, 0.3},  //CBB
  array[3]{0.1, 0.7, 0.2},  //CBC

  array[3]{0.4, 0.3, 0.3},  //CCA
  array[3]{0.6, 0.1, 0.3},  //CCB
  array[3]{0.1, 0.7, 0.2},  //CCC
};
#declare Order = 3;
#declare History = array[3]{0,0,0};
#declare Stream = seed(8);
#declare Mc_3 = MCMarkov(TransitionMatrix_3, Order, History, Stream);
#for (I, 0, Runs-1) 
  #debug concat("current history : ", Hist_3[MCgetHist(Mc_3)], "\n")
  #declare NewState = MCnextState(Mc_3);
  #debug concat("new State : ", State[NewState], "\n")
  #debug concat("new history : ", Hist_3[MCgetHist(Mc_3)], "\n\n")
#end


#end //#if(input_file_name="curve.inc")
#version MARKOV_INC_TEMP;
#end //#ifndef(MARKOV_INC_TEMP)