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

double grad;
double devi;
int action;
int count;

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

environ_ns* env;

/************************************************************************** 
Function: getUserNode.
This is a utility function to prompt user for a node from a list
**************************************************************************/

node_bn* getUserNode(const nodelist_bn* nodeList, char* type){
  int nodeNum, i;

  printf("Nodes in network:\n");
  for(i=0; i < LengthNodeList_bn(nodeList); i++) {
	printf("%d:\t%s\n", i+1, GetNodeName_bn (NthNode_bn (nodeList, i)));
  }

  printf("\nEnter %s Node Number: ", type);
  scanf("%d", &nodeNum);
  while(NthNode_bn(nodeList, nodeNum-1) == NULL) {
	printf("Invalid Node, try again: ");
	scanf("%d", &nodeNum);
  }
  return NthNode_bn(nodeList, nodeNum-1);
}

/************************************************************************** 
Function: getUserObservations.
This is a utility function to prompt user for a subset of a node list
**************************************************************************/

nodelist_bn* getUserObservations(nodelist_bn* nodeList, nodelist_bn* obsList){
  int nodeNum, i;

  printf("\nNodes in network:\n");
  for(i=0; i < LengthNodeList_bn(nodeList); i++) {
	printf("%d:\t%s\n", i+1, GetNodeName_bn(NthNode_bn(nodeList, i)));
  }
  printf("%d:\tStop\n", i+1);

  printf("\nEnter Observation Node Number (or option): ");
  scanf("%d", &nodeNum);

  if(nodeNum == i+1) return obsList;

  while(NthNode_bn(nodeList, nodeNum-1) == NULL) {
	printf("Invalid Node, try again: ");
	scanf("%d", &nodeNum);
  }

  AddNodeToList_bn(RemoveNthNode_bn(nodeList, nodeNum-1), obsList, 0);
  return getUserObservations(nodeList, obsList);
}

/************************************************************************** 
Function: userDone.
This is a utility function to prompt user to continue or finish
**************************************************************************/

int userDone() {
  char userAction[80];

  if(!action) {
	action = 1;
	return 0;
  }
  printf("Continue? (y/n): ");
  scanf("%s", userAction);
  if(!strcmp(userAction,"n")) return 1;
  return 0;
}

/************************************************************************** 
Function: userEnd.
This is a utility function just waits for the user to finish
**************************************************************************/

void userEnd() {
  char userAction[80];

  printf("Finish? (y/n): ");
  scanf("%s", userAction);
}

/************************************************************************** 
Function: listEvidence.
This is a utility function to generate list of evidence set
**************************************************************************/

void listEvidence(char list[200], nodelist_bn *evidenceSet) {
  int i;
  int num_evidence = LengthNodeList_bn(evidenceSet);
  node_bn *tempNode;
  char temp[20];

  strcpy(list, "");
  if(num_evidence == 0) strcat(list, "Empty");
  for(i=0; i<num_evidence; i++) {
	tempNode = NthNode_bn(evidenceSet, i);
	sprintf(temp, "%s = %s", GetNodeName_bn(tempNode), GetNodeStateName_bn(tempNode, GetNodeFinding_bn(tempNode)));
	strcat(list, temp);
	if(i!=num_evidence-1) strcat(list, ", ");
  }
}

/************************************************************************** 
Function: listParentStates.
This is a utility function to generate list of parent state of probability distribution
**************************************************************************/

void listParentStates(char list[200], const nodelist_bn* parents, state_bn *parent_states) {
  int i;
  int num_parents = LengthNodeList_bn(parents);

  strcpy(list, "");
  for(i=0; i<num_parents; i++) {
	strcat(list, GetNodeName_bn(NthNode_bn(parents, i)));
	strcat(list, "=");
	strcat(list, GetNodeStateName_bn(NthNode_bn(parents, i), parent_states[i]));
	if(i!=(num_parents-1)) strcat(list, ", ");
  }
}

/************************************************************************** 
Function: copyProb.
This is a utility function to create copy of a probability distribution
**************************************************************************/

prob_bn* copyProb(const prob_bn* orig, int size) {
  int i;
  prob_bn* new = (prob_bn*) malloc (size * sizeof (prob_bn));

  for(i=0; i<size; i++) {
	new[i] = orig[i];
  }
  return new;
}

/************************************************************************** 
Function: greater.
This is a utility function to return greater of two doubles
**************************************************************************/

double greater(double a, double b) {
  if(a>=b) return a;
  else return b;
}

/************************************************************************** 
Function: normalise.
This is a function to normalise the probability distribution after a parameter has been altered
**************************************************************************/

void normalise (prob_bn* probDist, int distSize, int fixed, double new) {
  int i;
  double old = probDist[fixed];

  probDist[fixed] = new;
  for(i=0; i<distSize; i++) {
	if(i!=fixed) probDist[i] = probDist[i]*(1-new)/(1-old);
  }
}

/************************************************************************** 
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: SolveLinear.
This is a function to solve set of linear equations
**************************************************************************/

int SolveLinear(double *A, double *B, int n) {
  int i, j, k;
  double temp;

  for(j=0; j < n-1; j++) {
	temp = A[j*n+j];
	k = j;
	for(i=j+1; i < n; i++) {
		if(A[i*n+j] > temp) {
			temp = A[i*n+j];
			k = i;
		}
	}
	if(k != j) {
		for(i=j; i < n; i++) {
			temp = A[j*n+i];
			A[j*n+i] = A[k*n+i];
			A[k*n+i] = temp;
		}
		temp = B[j];
		B[j] = B[k];
		B[k] = temp;
	}
	if(A[j*n+j] == 0) return 0;
	for(i=j+1; i < n; i++) {
		temp = A[i*n+j]/A[j*n+j];
		for(k=j+1; k < n; k++) {
			A[i*n+k] = A[i*n+k] - temp*A[j*n+k]; 
		}
		B[i] = B[i] - temp*B[j];
	}
  }
  if(A[(n-1)*n+(n-1)] == 0) return 0;
  B[n-1] = B[n-1]/A[(n-1)*n+(n-1)];
  for(i=n-2; i >= 0; i--) {
	temp = 0;
	for(j=i+1; j < n; j++) temp += A[i*n+j]*B[j];
	B[i] = (B[i] - temp)/A[i*n+i];
  }
  return 1;
}

/************************************************************************** 
Function: descendEvideCheck.
This is a function check if a node has an observed descendent
**************************************************************************/

int descendEvideCheck(node_bn *current){
  int i;
  const nodelist_bn *children = GetNodeChildren_bn(current);
  int num_children = LengthNodeList_bn(children);

  if(GetNodeFinding_bn(current) != NO_FINDING) return 1;

  for(i=0; i<num_children; i++) {
	if(descendEvideCheck(NthNode_bn(children, i))) return 1;
  }
  return 0;
}

/************************************************************************** 
Function: parentDSep.
This function is part of simple d-separation test algorithm it searchs though ancestor nodes
**************************************************************************/

int parentDSep(node_bn *current, node_bn *find, nodelist_bn* path) {
  int i, j;
  int flag = 0;
  const nodelist_bn *parents, *children;
  int num_parents, num_children, num_path;

  AddNodeToList_bn(current, path, 0);
  num_path = LengthNodeList_bn(path);

  if(!strcmp(GetNodeName_bn(current), GetNodeName_bn(find))) return 0;
  if(GetNodeFinding_bn(current) != NO_FINDING) return 1;

  parents = GetNodeParents_bn(current);
  num_parents = LengthNodeList_bn(parents);
  children = GetNodeChildren_bn(current);
  num_children = LengthNodeList_bn(children);

  for(i=0; i<num_parents; i++) {
	for(j=0; j<num_path; j++) {
		if(!strcmp(GetNodeName_bn(NthNode_bn(parents, i)), GetNodeName_bn(NthNode_bn(path, j)))) {
			flag = 1;
			break;
		}
	}
	if(flag) flag = 0;
	else if(!parentDSep(NthNode_bn(parents, i), find, path)) return 0;
  }

  for(i=0; i<num_children; i++) {
	for(j=0; j<num_path; j++) {
		if(!strcmp(GetNodeName_bn(NthNode_bn(children, i)), GetNodeName_bn(NthNode_bn(path, j)))) {
			flag = 1;
			break;
		}
	}
	if(flag) flag = 0;
	else if(!childDSep(NthNode_bn(children, i), find, path)) return 0;
  }
  RemoveNthNode_bn(path, 0);
  return 1;
}

/************************************************************************** 
Function: childDSep.
This function is part of simple d-separation algorithm it searchs though descendant nodes
**************************************************************************/

int childDSep(node_bn *current, node_bn *find, nodelist_bn* path) {
  int i, j;
  int flag = 0;
  int action = 1;
  int desendCheck;
  const nodelist_bn *parents, *children;
  int num_parents, num_children, num_path;

  AddNodeToList_bn(current, path, 0);
  num_path = LengthNodeList_bn(path);

  if(!strcmp(GetNodeName_bn(current), GetNodeName_bn(find))) return 0;
  if(GetNodeFinding_bn(current) != NO_FINDING) action = 0;

  parents = GetNodeParents_bn(current);
  num_parents = LengthNodeList_bn(parents);
  children = GetNodeChildren_bn(current);
  num_children = LengthNodeList_bn(children);

  if(action) {
  	for(i=0; i<num_children; i++) {
		for(j=0; j<num_path; j++) {
			if(!strcmp(GetNodeName_bn(NthNode_bn(children, i)), GetNodeName_bn(NthNode_bn(path, j)))) {
				flag = 1;
				break;
			}
		}
		if(flag) flag = 0;
		else {
			desendCheck = childDSep(NthNode_bn(children, i), find, path);
			if(!desendCheck) return 0;
			if(desendCheck == -1) action = 0;
		}
  	}
  }

  if(!action) {
  	for(i=0; i<num_parents; i++) {
		for(j=0; j<num_path; j++) {
			if(!strcmp(GetNodeName_bn(NthNode_bn(parents, i)), GetNodeName_bn(NthNode_bn(path, j)))) {
				flag = 1;
				break;
			}
		}
		if(flag) flag = 0;
		else if(!parentDSep(NthNode_bn(parents, i), find, path)) return 0;
  	}
  }

  RemoveNthNode_bn(path, 0);
  if(!action) return -1;
  else return 1;
}

/************************************************************************** 
Function: dSep.
This function is part of simple d-separation algorithm
it searchs though ancestor and descendents nodes trying to
find a path to the end node
**************************************************************************/

int dSep(node_bn *start, node_bn *end) {
  int i;
  nodelist_bn *path = NewNodeList_bn(0, env);
  const nodelist_bn *parents, *children;
  int num_parents, num_children;

  AddNodeToList_bn(start, path, LAST_ENTRY);

  if(!strcmp(GetNodeName_bn(start), GetNodeName_bn(end))) return 0;
  if(GetNodeFinding_bn(start) != NO_FINDING) return 1;

  parents = GetNodeParents_bn(start);
  num_parents = LengthNodeList_bn(parents);
  children = GetNodeChildren_bn(start);
  num_children = LengthNodeList_bn(children);

  for(i=0; i<num_parents; i++) {
	if(!parentDSep(NthNode_bn(parents, i), end, path)) return 0;
  }
  for(i=0; i<num_children; i++) {
	if(!childDSep(NthNode_bn(children, i), end, path)) return 0;
  }
  return 1;
}

/************************************************************************** 
Function: testIfSensitParam.
This is a function to determine if a node is in the sensitivity set
it uses parts of the dSep algorithm to this 
**************************************************************************/

int testIfSensitParam(node_bn *start, node_bn *end) {
  int i;
  nodelist_bn *path = NewNodeList_bn(0, env);

  if(!childDSep(start, end, path)) return 1;
  return 0;
}

/************************************************************************** 
Function: getHypVars.
This is a function compute the coeffiecents of a hyperbolic function 
**************************************************************************/

int getHypVars(double *vars, node_bn *qNode, node_bn *tNode, int index, state_bn *parent_states, net_bn *net) {
  int i, j;
  int t_num_states = GetNodeNumberStates_bn(tNode);
  int q_num_states = GetNodeNumberStates_bn(qNode);
  const prob_bn *bel;
  const prob_bn *origProbDist;
  prob_bn* tempProbDist;
  double *A = (double*) calloc(3*3, sizeof(double));
  double *B = (double*) calloc(3, sizeof(double));
  double temp;

  bel = GetNodeBeliefs_bn (qNode);
  origProbDist = copyProb(GetNodeProbs_bn (tNode, parent_states), t_num_states);

  for(j=0; j<q_num_states; j++) {
	A[0*3+0] = origProbDist[index];
	A[0*3+1] = 1.0;
	A[0*3+2] = -bel[j];
	B[0] = bel[j]*origProbDist[index];

	for(i=1; i<3; i++) {
		temp = origProbDist[index] + i*0.33;
		if(temp > 1) temp = temp-1.0;
  		tempProbDist = copyProb(origProbDist, t_num_states);
  		normalise(tempProbDist, t_num_states, index, temp);
  		SetNodeProbs_bn (tNode, parent_states, tempProbDist);
  		CompileNet_bn (net);
		A[i*3+0] = temp;
		A[i*3+1] = 1.0;
		A[i*3+2] = -bel[j];
		B[i] = bel[j]*temp;
  		free(tempProbDist);
	}
	SetNodeProbs_bn(tNode, parent_states, origProbDist);
	if(!SolveLinear(A, B, 3)) return 0;
	for(i=0; i<3; i++) {
		vars[i*3+j] = B[i];
		if(!finite(vars[i*3+j])) return 0;
	}
  }
  return 1;
}

/************************************************************************** 
Function: needHypPlot.
This is a function test if hyperbolic function will be of interest to the user 
**************************************************************************/

int needHypPlot(double *vars, int num_plots, int index, node_bn *tNode, state_bn *parent_states) {
  int i;
  int t_num_states = GetNodeNumberStates_bn(tNode);
  const prob_bn *origProbDist = copyProb(GetNodeProbs_bn (tNode, parent_states), t_num_states);
  double temp;

  for(i=0; i<num_plots; i++) {
	temp = sqrt(fabs(-vars[1*3+i]+vars[0*3+i]*vars[2*3+i]));
	if(fabs(fabs(greater(-vars[2*3+i]+temp, -vars[2*3+i]-temp))-origProbDist[index])<=devi) return 1;
	else if(fabs((vars[0*3+i]*vars[2*3+i]-vars[1*3+i])/(pow(origProbDist[index]+vars[2*3+i],2)))>=grad) return 1;
  }
  return 0;
}

/************************************************************************** 
Function: getHypPlot.
This is a function uses the preceeding functions to compute a hyperbolic
function, determine if it is of interest and display it using gnuplot
if it is 
**************************************************************************/

int getHypPlot(node_bn *qNode, node_bn *tNode, int index, state_bn *parent_states, nodelist_bn *evidenceSet, net_bn *net, FILE* plots) {
  int i;
  int q_num_states = GetNodeNumberStates_bn(qNode);
  int t_num_states = GetNodeNumberStates_bn(tNode);
  const prob_bn *origProbDist = copyProb(GetNodeProbs_bn (tNode, parent_states), t_num_states);
  double *vars = (double*) calloc(3*q_num_states, sizeof(double));
  char label[200];
  char list[200];
  double x[2];
  double y[2];
  gnuplot_ctrl * h1;

  if(!getHypVars(vars, qNode, tNode, index, parent_states, net)) return 0;
  if(needHypPlot(vars, q_num_states, index, tNode, parent_states)) {
	if(userDone()) exit(1);
  	h1 = gnuplot_init();
  	gnuplot_resetplot(h1);

	listEvidence(list, evidenceSet);
	sprintf(label,"Pr(%s | %s)", GetNodeName_bn(qNode), list);
	gnuplot_set_ylabel(h1, label);

	listParentStates(list, GetNodeParents_bn(tNode), parent_states);
	sprintf(label,"Pr(%s=%s | %s)", GetNodeName_bn(tNode), GetNodeStateName_bn(tNode, index), list);
	gnuplot_set_xlabel(h1, label);

  	gnuplot_setstyle(h1, "points");
  	x[0] = 0; x[1] = 1; y[0] = 1; y[1] = 0;
  	gnuplot_plot_xy(h1, x, y, 2, "");

  	gnuplot_setstyle(h1, "lines");
  	x[0] = x[1] = origProbDist[index]; y[0] = 0; y[1] = 1;
  	gnuplot_plot_xy(h1, x, y, 2, "Initial");

	for(i=0; i<q_num_states; i++) {
  		sprintf(list,"(%lf*x+%lf)/(x+%lf)", vars[0*3+i], vars[1*3+i], vars[2*3+i]);
		strcpy(label, GetNodeStateName_bn(qNode, i));
  		gnuplot_plot_equation(h1, list, label);
	}
  }
  return 1;
}

/************************************************************************** 
Function: HypFuncs.
This is a function runs the hyperbolic function generator for all 
parameters in the CPT
**************************************************************************/

void HypFuncs (node_bn *iNode, node_bn *tNode, nodelist_bn *evidenceSet, net_bn *net, FILE* plots) {
  int i, j;
  int t_num_states = GetNodeNumberStates_bn(tNode);
  const nodelist_bn* parents = GetNodeParents_bn (tNode);
  int num_parents = LengthNodeList_bn (parents);
  state_bn* parent_states = (state_bn*) calloc (num_parents, sizeof (state_bn));

  for(j=0; j<num_parents; j++) {
	parent_states[j] = 0;
  }
  for(i=0; i<t_num_states; i++) {
	if(!getHypPlot(iNode, tNode, i, parent_states, evidenceSet, net, plots)) getLinPlot(iNode, tNode, i, parent_states, evidenceSet, net, plots);
  }
  while(!NextStates (parent_states, parents)) {
  	for(i=0; i<t_num_states; i++) {
		if(!getHypPlot(iNode, tNode, i, parent_states, evidenceSet, net, plots)) getLinPlot(iNode, tNode, i, parent_states, evidenceSet, net, plots);
  	}
  }
}

/************************************************************************** 
Function: getLinVars.
This is a function compute the coeffiecents of a linear function 
**************************************************************************/

int getLinVars(double *vars, node_bn *qNode, node_bn *tNode, int index, state_bn *parent_states, net_bn *net) {
  int i, j;
  int t_num_states = GetNodeNumberStates_bn(tNode);
  int q_num_states = GetNodeNumberStates_bn(qNode);
  const prob_bn *bel;
  const prob_bn *origProbDist;
  prob_bn* tempProbDist;
  double *A = (double*) calloc(2*2, sizeof(double));
  double *B = (double*) calloc(2, sizeof(double));
  double temp;

  bel = GetNodeBeliefs_bn (qNode);
  origProbDist = copyProb(GetNodeProbs_bn (tNode, parent_states), t_num_states);
  for(j=0; j<q_num_states; j++) {
	A[0*2+0] = origProbDist[index];
	A[0*2+1] = 1.0;
	B[0] = bel[j];

	temp = origProbDist[index]+0.5;
	if(temp > 1) temp = temp-1.0;
  	tempProbDist = copyProb (origProbDist, t_num_states);
  	normalise(tempProbDist, t_num_states, index, temp);
  	SetNodeProbs_bn (tNode, parent_states, tempProbDist);
  	CompileNet_bn (net);
	A[1*2+0] = temp;
	A[1*2+1] = 1.0;
	B[1] = bel[j];
  	free(tempProbDist);
	SetNodeProbs_bn (tNode, parent_states, origProbDist);

	if(!SolveLinear(A, B, 2)) return 0;
	for(i=0; i<2; i++) {
		vars[i*2+j] = B[i];
		if(!finite(vars[i*2+j])) return 0;
	}
  }
  return 1;
}

/************************************************************************** 
Function: needLinPlot.
This is a function test if linear function will be of interest to the user 
**************************************************************************/

int needLinPlot(double *vars, int num_plots, node_bn *tNode, state_bn *parent_states) {
  int i;
  int t_num_states = GetNodeNumberStates_bn(tNode);
  const prob_bn *origProbDist = copyProb(GetNodeProbs_bn (tNode, parent_states), t_num_states);

  for(i=0; i<num_plots; i++) {
	if(fabs(vars[0*2+i])>=grad) return 1;
  }
  return 0;
}

/************************************************************************** 
Function: getLinPlot.
This is a function uses the preceeding functions to compute a linear
function, determine if it is of interest and display it using gnuplot
if it is 
**************************************************************************/

int getLinPlot(node_bn *qNode, node_bn *tNode, int index, state_bn *parent_states, nodelist_bn *evidenceSet, net_bn *net, FILE *plots) {
  int i;
  int q_num_states = GetNodeNumberStates_bn(qNode);
  int t_num_states = GetNodeNumberStates_bn(tNode);
  const prob_bn *origProbDist = copyProb(GetNodeProbs_bn (tNode, parent_states), t_num_states);
  double *vars = (double*) calloc(2*q_num_states, sizeof(double));
  char label[200];
  char list[200];
  double x[2];
  double y[2];
  gnuplot_ctrl * h1;

  if(!getLinVars(vars, qNode, tNode, index, parent_states, net)) return 0;
  if(needLinPlot(vars, q_num_states, tNode, parent_states)) {
	if(userDone()) exit(1);
  	h1 = gnuplot_init();
  	gnuplot_resetplot(h1);

	listEvidence(list, evidenceSet);
	sprintf(label,"Pr(%s|%s)", GetNodeName_bn(qNode), list);
	gnuplot_set_ylabel(h1, label);

	listParentStates(list, GetNodeParents_bn(tNode), parent_states);
	sprintf(label,"Pr(%s=%s | %s)", GetNodeName_bn(tNode), GetNodeStateName_bn(tNode, index), list);
	gnuplot_set_xlabel(h1, label);

  	gnuplot_setstyle(h1, "points");
  	x[0] = 0; x[1] = 1; y[0] = 1; y[1] = 0;
  	gnuplot_plot_xy(h1, x, y, 2, "");

  	gnuplot_setstyle(h1, "lines");
  	x[0] = x[1] = origProbDist[index]; y[0] = 0; y[1] = 1;
  	gnuplot_plot_xy(h1, x, y, 2, "Initial");

	for(i=0; i<q_num_states; i++) {
  		sprintf(list,"%lf*x+%lf", vars[0*2+i], vars[1*2+i]);
		strcpy(label, GetNodeStateName_bn (qNode, i));
  		gnuplot_plot_equation(h1, list, label);
	}
  }
  return 1;
}

/************************************************************************** 
Function: LinFuncs.
This is a function runs the Linear function generator for all 
parameters in the CPT
**************************************************************************/

void LinFuncs (node_bn *iNode, node_bn *tNode, nodelist_bn *evidenceSet, net_bn *net, FILE* plots) {
  int i, j;
  int t_num_states = GetNodeNumberStates_bn(tNode);
  const nodelist_bn* parents = GetNodeParents_bn (tNode);
  int num_parents = LengthNodeList_bn (parents);
  state_bn* parent_states = (state_bn*) calloc (num_parents, sizeof (state_bn));

  for(j=0; j<num_parents; j++) {
	parent_states[j] = 0;
  }
  for(i=0; i<t_num_states; i++) {
	getLinPlot(iNode, tNode, i, parent_states, evidenceSet, net, plots);
  }
  while(!NextStates (parent_states, parents)) {
  	for(i=0; i<t_num_states; i++) {
		getLinPlot(iNode, tNode, i, parent_states, evidenceSet, net, plots);
  	}
  }
}

/************************************************************************** 
Function: sensit.
This is a function drives the sensitivity analysis for an evidence instance
**************************************************************************/

void sensit(node_bn *iNode, node_bn *tNode, nodelist_bn *evidenceSet, net_bn* net, FILE* plots) {
  char list[200];

listEvidence(list, evidenceSet);
printf("");
  if(testIfSensitParam(tNode, iNode)) {
	if(!descendEvideCheck(tNode)) LinFuncs(iNode, tNode, evidenceSet, net, plots);
	else HypFuncs(iNode, tNode, evidenceSet, net, plots);
  }
}

/************************************************************************** 
Function: evideSet.
This is a function to loop though the evidence set
**************************************************************************/

void evideSet(node_bn *iNode, node_bn *tNode, nodelist_bn *remNodes, nodelist_bn *evidenceSet, net_bn* net, FILE* plots) {
  int i, j;
  node_bn *tempNode;
  int num_nodes = LengthNodeList_bn(remNodes);
  int num_states;

  for(i=0; i<num_nodes; i++) {
	if(!dSep(NthNode_bn(remNodes, i), iNode)) {
		tempNode = RemoveNthNode_bn(remNodes, i);
		evideSet(iNode, tNode, remNodes, evidenceSet, net, plots);
		num_states = GetNodeNumberStates_bn (tempNode);
		for(j=0; j<num_states; j++) {
			EnterFinding_bn(tempNode, j);
			AddNodeToList_bn(tempNode, evidenceSet, LAST_ENTRY);
			sensit(iNode, tNode, evidenceSet, net, plots);
			evideSet(iNode, tNode, remNodes, evidenceSet, net, plots);
			RetractNodeFindings_bn(tempNode);
			RemoveNthNode_bn(evidenceSet, LAST_ENTRY);
		}
		AddNodeToList_bn(tempNode, remNodes, i);
		break;
	}
  }
}

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

int main (void){
  net_bn* net;
  char mesg[MESG_LEN_ns];
  int res;
  char netName[80];
  double temp;
  node_bn *iNode, *tNode;
  action = 0;
  count = 0;
  nodelist_bn *allNodes;
  nodelist_bn *evidenceSet;
  nodelist_bn *obsSet;
  FILE* plots;
	
  env = NewNeticaEnviron_ns (NULL, NULL, NULL);
  res = InitNetica_bn (&env, mesg);
  printf ("%s\n", mesg);
  if (res < 0)  exit (-1);

  printf("Enter Network Name: ");
  scanf("%s", netName);
  strcat(netName, ".dne");
  net = ReadNet_bn (NewStreamFile_ns (netName, env, NULL), NO_WINDOW);
  CHKERR

  iNode = getUserNode(GetNetNodes_bn(net), "Interest");
  printf("Enter Minimum Variation of Interest Node Probabilies (0.0 - 100.0): ");
  scanf("%lf", &temp);
  devi = (double) temp/100.0;

  tNode = getUserNode(GetNetNodes_bn(net), "Test");
  printf("Enter Maximum Variation of Parameters in Test Node (0.0 - 100.0): ");
  scanf("%lf", &temp);
  grad = (double) devi/(temp/100);

  printf("Gradient of interest: %lf\n", grad);
  printf("Deviation of interest: %lf\n", devi); 

  CompileNet_bn (net);

  allNodes = DupNodeList_bn(GetNetNodes_bn(net));
  evidenceSet = NewNodeList_bn(0, env);
  obsSet = NewNodeList_bn(0, env);
  RemoveNthNode_bn(allNodes, IndexOfNodeInList_bn(iNode, allNodes, 0));
  getUserObservations(allNodes, obsSet);

  sensit(iNode, tNode, evidenceSet, net, plots);
  evideSet(iNode, tNode, obsSet, evidenceSet, net, plots);
  userEnd();

end:
	FreeNet_bn (net);
	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;
}