/*============================================================================

 Function      vfStep

 Usage         vfStep(klo, khi, nksub, kslo, kshi, ksmid1, ksmid2, w, wmax,
                      windmax, w1, w2, w3, i, j, l, ydepK, howard, K, Z, P,
		      Exp, V0, V, G)

 Arguments     klo:     reference to integer which represents the index
                        corresponding to the lowest value of the capital grid
			over which to maximize.

               khi:     reference to integer which represents the index
                        corresponding to the highest value of the capital
			grid over which to maximize.

	       nksub:   reference to the integer length of the subgrid of
	                capital (beginning at klo) over which to maximize.

	       kslo:    reference to integer representing the index
	                corresponding to the lower bound of the region within
		        the subgrid of capital which the binary search method
		        isolates for maximization.

	       kshi:    reference to integer representing the index
	                corresponding to the upper bound of the region within
		        the subgrid of capital which the binary search method
		        isolates for maximization.

	       ksmid1:  reference to integer representing the index
	                corresponding to one of two midpoints of the region
		        within the subgrid of capital which the binary search
		        method isolates for maximization.

	       ksmid1:  reference to integer representing the index
	                corresponding to one of two midpoints of the region
		        within the subgrid of capital which the binary search
		        method isolates for maximization. ksmid2 = ksmid1 + 1.

	       w:       reference to double used to store value of utility
	                at each value of capital subgrid (grid search).

	       wmax:    reference to double that tracks the max of utility
	                over the capital grid (grid search).

	       windmax: reference to integer that tracks the index of the
                        capital grid that corresponds to the argmax (grid
			search).

	       w1:      reference to double used to store the value of utility
	                at a particular value of capital during binary search.

	       w2:      reference to double used to store the value of utility
	                at a particular value of capital during binary search.

	       w3:      reference to double used to store the value of utility
	                at a particular value of capital during binary search.
		        This variable is only used when the subgrid of capital
		        has three values.

	       i:       reference to integer value used for incrementing.

	       j:       reference to integer value used for incrementing.

	       l:       reference to integer value used for incrementing.

	       ydepK:   reference to double which stores the value of output
	                plus capital net of depreciation.

	       howard:  reference to constant boolean that indicates if the
	                current iteration of the value function will perform
			a maximization (false) or if it will simply compute
			the new value function using the concurrent policy
			function (true).

	       K:       pointer to array of constant doubles which stores the
	                grid of capital values.

	       Z:       pointer to array of constant doubles which stores the
	                TFP grid values.
	               
	       P:       pointer to array of constant doubles which stores the
	                TFP transition matrix values.
	              
	       Exp:     pointer to array of doubles which stores the
	                continuation value of utility at each value of capital
			subgrid (grid search).		     

	       V0:      pointer to the array of constant doubles which
	                represents the initial value function.

	       V:       pointer to the array of doubles which represents the 
	                value function.

	       G:       pointer to the array of doubles which represents the 
	                policy function.

 Description   This function performs one iteration of the value function
               iteration algorithm, using V0 as the current value function and
	       either maximizing the LHS of the Bellman if howard = false or
	       using the concurrent policy function as the argmax if
	       howard = true. Monotonicity of the policy function is exploited
	       to restrict the region of maximization. Maximization is
	       performed by either a grid search or binary search algorithm.

 Dependencies  Global variables: eta, beta, alpha, delta, nk, nz,
                                 maxtype, howard (globalvars.h).

               Functions:        pow (math.h);
				 cblas_dgemm, cblas_ddot (cblas.h);
	                         binary_val, grid_max, binary_max (auxfuncs.h).

 Return value  void.

 =============================================================================

 Author:       Eric M. Aldrich

 Contact:      ealdrich@gmail.com

 Date:         20 April 2010

 ============================================================================*/

#include "globalvars.h"
#include "auxfuncs.h"
#include <math.h>
#include <iostream>
extern "C" {
  #include "cblas.h"
}

using namespace std;

void vfStep(int& klo, int& khi, int& nksub, int& kslo, int& kshi, int& ksmid1,
	    int& ksmid2, double& w, double& wmax, int& windmax, double& w1,
	    double& w2, double& w3, int& i, int& j, int& l, double& ydepK,
	    const bool& howard, const double* K,const double* Z,
	    const double* P, double* Exp, const double* V0, double* V,
	    double* G)
{
  for(i = 0 ; i < nk ; ++i){
    for(j = 0 ; j < nz ; ++j){

      // output and depreciated capital
      ydepK = Z[j]*pow(K[i],alpha) + (1-delta)*K[i];

      // maximize on non-howard steps
      if(howard == false){

	// impose constraints on grid for future capital
	klo = 0;
	khi = binary_val(ydepK, nk, K); // consumption nonnegativity
	if(K[khi] > ydepK) khi -= 1;

	// further restrict capital grid via monotonicity (CPU only)
	if(i > 0){
	  if(G[(i-1)*nz+j] > klo & G[(i-1)*nz+j] < khi) klo = (int)G[(i-1)*nz+j];
	}
	nksub = khi-klo+1;

	// continuation value for subgrid
	Exp = (double*)realloc(Exp, nksub*sizeof(double));
	cblas_dgemv(CblasRowMajor, CblasNoTrans, nksub, nz, 1.0, (V0+klo*nz),
		    nz, (P+j*nz), 1, 0.0, Exp, 1);

	// maximization either via grid (g), of binary search (b)
	// if binary, turn off policy iteration (to preserve concavity)
	if(maxtype == 'g'){
	  grid_max(klo, nksub, l, w, wmax, windmax, ydepK, K, Exp, V+i*nz+j,
		   G+i*nz+j);
	} else if (maxtype == 'b'){
	  binary_max(klo, nksub, kslo, kshi, ksmid1, ksmid2, w1, w2, w3,
		     ydepK, K, Exp, V+i*nz+j, G+i*nz+j);
	}

      // iterate on the policy function on non-howard steps
      } else {
	Exp = (double*)realloc(Exp, sizeof(double));
	Exp[0] = cblas_ddot(nz, (V0+(int)G[i*nz+j]*nz), 1, (P+j*nz), 1);
	V[i*nz+j] = pow(ydepK-K[(int)G[i*nz+j]],1-eta)/(1-eta) + beta*Exp[0];
      }
    }
  }
}
