#include <stdio.h>
#include <string.h>
#include <math.h>
#include "Netica.h"

#define CHKERR  {if (GetError_ns (env, ERROR_ERR, NULL))  goto error;}

environ_ns* env;

/************************************************************************** 
Function: NextStates.
This is a utility function to return loop through rows of CPT
**************************************************************************/

int NextStates (state_bn* states, const nodelist_bn* nodes){
  int nn;
  for (nn = LengthNodeList_bn (nodes) - 1;  nn >= 0;  --nn) {
	if (++states[nn] < GetNodeNumberStates_bn (NthNode_bn (nodes, nn))) return 0;
	states[nn] = 0;
	}
  return 1;
}

/************************************************************************** 
Function: BhattDist.
This function Computes the Bhattarcharyya Distance between two probability
distributions
**************************************************************************/

double BhattDist(const prob_bn *probDist1, const prob_bn *probDist2, int n) {
  int i;
  double temp = 0.0;

  for(i=0; i<n; i++) {
	temp += (double) sqrt(probDist1[i]*probDist2[i]);
  }
  temp = -log(temp);
  printf("%.4lf\n", temp);
  return temp;
}

/************************************************************************** 
Function: doNets.
This function calls the Bhattarcharyya function for all probability 
distributions in the nodes of interest, which are hard coded in at the 
moment.  It identify probility distributions that have a significant difference
along with the actual distance 
**************************************************************************/

void doNets(net_bn *net1, net_bn *net2) {
  int i, j, n = 12;
  char interestNodes[80][80] = {
	{"Temp"},
	{"WQmain"},
	{"LowFlowSp"},
	{"CommCh"},
	{"Floodplain"},
	{"BiolPoten"},
	{"StrucHab"},
	{"Connectivity"},
	{"CurrentDiversity"},
	{"OverallFlow"},
	{"FutureAbundance"},
	{"FutureDiversity"}};
  node_bn *node1, *node2;
  state_bn* parent_states;
  const prob_bn *probDist1, *probDist2;

  for(i=0; i<n; i++) {
	node1 = NodeNamed_bn(interestNodes[i], net1);
	node2 = NodeNamed_bn(interestNodes[i], net2);
	printf("\nNode Name: %s\n", GetNodeName_bn(node1));

	parent_states = (state_bn*) calloc(LengthNodeList_bn(GetNodeParents_bn(node1)), sizeof(state_bn));
	for(j=0; j<LengthNodeList_bn(GetNodeParents_bn(node1)); j++) parent_states[j] = 0;
	printf("Parents: ");
	for(j=0; j<LengthNodeList_bn(GetNodeParents_bn(node1)); j++) {
		printf("%s\t", GetNodeName_bn(NthNode_bn(GetNodeParents_bn(node1), j)));
  	}
	printf("\nStates: ");
  	for(j=0; j<LengthNodeList_bn(GetNodeParents_bn(node1)); j++) {
		printf("%s\t", GetNodeStateName_bn(NthNode_bn(GetNodeParents_bn(node1), j), parent_states[j]));
  	}
	printf(":\t");
  	probDist1 = GetNodeProbs_bn(node1, parent_states);
  	probDist2 = GetNodeProbs_bn(node2, parent_states);
	BhattDist(probDist1, probDist2, GetNodeNumberStates_bn(node1));
	while(!NextStates(parent_states, GetNodeParents_bn(node1))) {
		printf("States: ");
  		for(j=0; j<LengthNodeList_bn(GetNodeParents_bn(node1)); j++) {
			printf("%s\t", GetNodeStateName_bn(NthNode_bn(GetNodeParents_bn(node1), j), parent_states[j]));
  		}
		printf(":\t");
  		probDist1 = GetNodeProbs_bn(node1, parent_states);
  		probDist2 = GetNodeProbs_bn(node2, parent_states);
		BhattDist(probDist1, probDist2, GetNodeNumberStates_bn(node1));
	}
	free(parent_states);
  }
}

/************************************************************************** 
main.
This function drives the overall program
**************************************************************************/

int main(void) {
  net_bn *net1, *net2;
  char mesg[MESG_LEN_ns];
  int res;
  char netName1[80];
  char netName2[80];
	
  env = NewNeticaEnviron_ns (NULL, NULL, NULL);
  res = InitNetica_bn (&env, mesg);
  printf ("%s\n", mesg);
  if (res < 0)  exit (-1);
	
  printf("Enter First Network Name: ");
  scanf("%s", netName1);
  strcat(netName1, ".dne");

  net1 = ReadNet_bn (NewStreamFile_ns (netName1, env, NULL), NO_WINDOW);
  CHKERR

  printf("Enter Second Network Name: ");
  scanf("%s", netName2);
  strcat(netName2, ".dne");

  net2 = ReadNet_bn (NewStreamFile_ns (netName2, env, NULL), NO_WINDOW);
  CHKERR

  printf("%s->%s\n", netName1, netName2);

  doNets(net1, net2);

end:
	FreeNet_bn (net1);
	FreeNet_bn (net2);
	res = CloseNetica_bn (env, mesg);
	printf ("%s\n", mesg);
	return (res < 0 ? -1 : 0);

error:
	fprintf (stderr, "Sensitivity Analysis: Error in %s\n", 
	         ErrorMessage_ns (GetError_ns (env, ERROR_ERR, NULL)));
	goto end;
}