// parse_func.c
/* This file contains the data readin routines for the simann.c program. */

#include "lsq.h"

extern char errstr[200];  // error string that all routines have access to

// file lines that will be searched for during the simulated annealing parameters read-in process
#define number_of_iterations_line "Niterations"
#define number_of_iterations_per_batch_line "Nperbatch"
#define prescription_filename_line "prescFile"
#define beamlet_header_line "beamletHeaderFileName"
#define initial_beam_weights_filename_line "initBeamWeightFile"
#define dose_batch_base_name_line "doseBatchBaseName"
#define dose_batch_extension_line "doseBatchExtension"
#define weight_batch_base_name_line "weightBatchBaseName"
#define weight_batch_extension_line "weightBatchExtension"
#define obj_func_name_line "objFuncFileName"
#define mod_factor_line "modFactor"

int read_in_opt_parameters(OPT_PARAM *op, char *opt_header)
// Loads the optimization parameters that define how the optimization will take place.
{
	int num_iterations_flag = 0;
	int num_per_batch_flag = 0;
	int prescription_filename_flag = 0;
	int beamlet_header_flag = 0;
	int init_beam_weights_flag = 0;
	int dose_batch_base_name_flag = 0;
	int dose_batch_extension_flag = 0;
	int weight_batch_base_name_flag = 0;
	int weight_batch_extension_flag = 0;
	int obj_func_name_flag = 0;
	int mod_factor_flag = 0;
		
	char str[200];

	FILE *fid;

	// open the optimization header
	if ( (fid = fopen(opt_header,"r")) == NULL)
	{
		printf("%s\n",opt_header);
	  sprintf(errstr,"Could not open the optimization header file.\n");
	  return(FAILURE);
	}

	// read in the optimization parameters
	while (fscanf(fid,"%s\n",str) != EOF)  // read until the end of the file
	{
	  if (!strcmp(str,number_of_iterations_line))
	    if (fscanf(fid,"%d\n",&op->Niterations) != 1)
		{
			sprintf(errstr,"Failed to read-in the number of iterations from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      num_iterations_flag = 1;
	  else if (!strcmp(str,number_of_iterations_per_batch_line))
	    if (fscanf(fid,"%d\n",&op->Nperbatch) != 1)
		{
			sprintf(errstr,"Failed to read-in the number of iterations per batch from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      num_per_batch_flag = 1;
	  else if (!strcmp(str,prescription_filename_line))
	    if (fscanf(fid,"%s\n",op->prescription_filename) != 1)
		{
			sprintf(errstr,"Failed to read-in the prescription filename line from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      prescription_filename_flag = 1;
	  else if (!strcmp(str,beamlet_header_line))
	    if (fscanf(fid,"%s\n",op->beamlet_header_file) != 1)
		{
			sprintf(errstr,"Failed to read-in the beamlet header file line from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      beamlet_header_flag = 1;
	  else if (!strcmp(str,initial_beam_weights_filename_line))
	    if (fscanf(fid,"%s\n",op->initial_beam_weights_filename) != 1)
		{
			sprintf(errstr,"Failed to read-in the initial beam weights file from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      init_beam_weights_flag = 1;
	  else if (!strcmp(str,dose_batch_base_name_line))
	    if (fscanf(fid,"%s\n",op->dose_batch_base_name) != 1)
		{
			sprintf(errstr,"Failed to read-in the dose batch base name from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      dose_batch_base_name_flag = 1;
	  else if (!strcmp(str,dose_batch_extension_line))
	    if (fscanf(fid,"%s\n",op->dose_batch_extension) != 1)
		{
			sprintf(errstr,"Failed to read-in the dose batch base extension from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      dose_batch_extension_flag = 1;
	  else if (!strcmp(str,weight_batch_base_name_line))
	    if (fscanf(fid,"%s\n",op->weight_batch_base_name) != 1)
		{
			sprintf(errstr,"Failed to read-in the weight batch base name from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      weight_batch_base_name_flag = 1;
	  else if (!strcmp(str,weight_batch_extension_line))
	    if (fscanf(fid,"%s\n",op->weight_batch_extension) != 1)
		{
			sprintf(errstr,"Failed to read-in the weight batch base extension from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      weight_batch_extension_flag = 1;
	  else if (!strcmp(str,obj_func_name_line))
	    if (fscanf(fid,"%s\n",op->obj_func_name) != 1)
		{
			sprintf(errstr,"Failed to read-in the objective function filename from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      obj_func_name_flag = 1;
	  else if (!strcmp(str,mod_factor_line))
	    if (fscanf(fid,"%f\n",&op->modFactor) != 1)
		{
			sprintf(errstr,"Failed to read-in the modulation factor from %s file.",opt_header);
			return(FAILURE);
		}
	    else
	      mod_factor_flag = 1;
	  else
	  {   
			sprintf(errstr,"Unrecognized string in the %s file:\n%s\n",opt_header,str);
			return(FAILURE);   
	  }
	}
	fclose(fid);

    // confirm that all of the required lines have been found
	if (num_iterations_flag == 0)
	{
	  sprintf(errstr,"Unable to find the number of iterations line in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	if (num_per_batch_flag == 0)
	{
	  sprintf(errstr,"Unable to find the number of iterations per batch line in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (prescription_filename_flag == 0)
	{
	  sprintf(errstr,"Unable to find the prescription filename file in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (beamlet_header_flag == 0)
	{
	  sprintf(errstr,"Unable to find the beamlet header line in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (init_beam_weights_flag == 0)
	{
	  sprintf(errstr,"Unable to find the initial beam weights filename in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (dose_batch_base_name_flag == 0)
	{
	  sprintf(errstr,"Unable to find the dose batch base name in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (dose_batch_extension_flag == 0)
	{
	  sprintf(errstr,"Unable to find the dose batch extension in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (weight_batch_base_name_flag == 0)
	{
	  sprintf(errstr,"Unable to find the weight batch base name in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (weight_batch_extension_flag == 0)
	{
	  sprintf(errstr,"Unable to find the weight batch extension in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (obj_func_name_flag == 0)
	{
	  sprintf(errstr,"Unable to find the objective function file name in the file %s.\n",opt_header);
	  return(FAILURE);
	}
	else if (mod_factor_flag == 0)
	{
	  sprintf(errstr,"Unable to find the modulation factor in the file %s.\n",opt_header);
	  return(FAILURE);
	}

	return(SUCCESS);
}

// file lines that will be searched-for during the geometry read-in process
#define Ntissue_line			"Ntissue"
#define tissueNum_line			"tissueNum"
#define name_line				"name"
#define alpha_line				"alpha"
#define betaVPlus_line			"betaVPlus"
#define dVPlus_line				"dVPlus"
#define vPlus_line				"vPlus"
#define betaVMinus_line			"betaVMinus"
#define dVMinus_line			"dMinus"
#define vMinus_line				"vMinus"
#define betaPlus_line			"betaPlus"
#define dosePlusFilename_line	"dosePlusFilename"
#define betaMinus_line			"betaMinus"
#define doseMinusFilename_line  "doseMinusFilename"

int read_in_prescription(PRESC *presc, char *presc_filename)
// Prescription structure is read in
{
	char str[200];

	int i,j,dims[3],Nind;

	FILE *fid, *doseMinusFile, *dosePlusFile;

	// read in the prescription file
	if ( (fid = fopen(presc_filename,"r")) == NULL)  // open the prescription file
	{
	  sprintf(errstr,"Could not open the prescription file, %s",presc_filename);
	  return(FAILURE);
	}

	// read in the number of tissue types
	if (fscanf(fid,"%s\n",str) != 1)
	{
		sprintf(errstr,"Failed to read in first line of %s",presc_filename);
		return(FAILURE);
	}
	else if (strcmp(str,Ntissue_line))
	{
		sprintf(errstr,"First line of %s must be %s",presc_filename,Ntissue_line);
		return(FAILURE);
	}
	else if (fscanf(fid,"%d\n",&presc->Ntissue) != 1)
	// found Ntissue_line, so read-in next line
	{
		sprintf(errstr,"Failed to read in the number of tissues from %s",presc_filename);
		return(FAILURE);
	}

	// allocate memory for the tissue type structure
	if ((presc->array = (TISSUE *)malloc(sizeof(TISSUE)*presc->Ntissue)) == NULL)
	{
		sprintf(errstr,"Failed to allocate space for prescription tissues",presc_filename);
		return(FAILURE);
	}
	printf("allocated space for %d tissues\n",presc->Ntissue);

	// read-in the prescription parameters for each tissue type
	for (i=0;i<presc->Ntissue;i++)
	{
		// pop off the tissue_num line
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,tissueNum_line))
		{
			sprintf(errstr,"Missing tissueNum line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue number
		if (fscanf(fid,"%d\n",&j) != 1)
		{
			sprintf(errstr,"Could not read in tissue number for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}
		else if (i != j)
		{
			sprintf(errstr,"Tissue indexed by %d in file %s has an inconsistent tissue number of %d.\n",i,presc_filename,j);
			return(FAILURE);
		}

		printf("tissue number is %d\n",j);

		// pop off the name line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,name_line))
		{
			sprintf(errstr,"Missing name line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}

		// get the tissue name
		if (fgets(presc->array[i].name,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read in tissue name for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue name is %s\n",presc->array[i].name);

		// pop off the alpha line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,alpha_line))
		{
			sprintf(errstr,"Missing alpha line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue alpha (importance)
		if (fscanf(fid,"%f\n",&presc->array[i].alpha) != 1)
		{
			sprintf(errstr,"Could not read in tissue alpha for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue alpha is %f\n",presc->array[i].alpha);

		// pop off the betaVPlus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,betaVPlus_line))
		{
			sprintf(errstr,"Missing betaVPlus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue betaVPlus
		if (fscanf(fid,"%f\n",&presc->array[i].betaVPlus) != 1)
		{
			sprintf(errstr,"Could not read in tissue betaVPlus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue betaVPlus is %f\n",presc->array[i].betaVPlus);

		// pop off the dVPlus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,dVPlus_line))
		{
			sprintf(errstr,"Missing dVPlus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue dVPlus
		if (fscanf(fid,"%f\n",&presc->array[i].dVPlus) != 1)
		{
			sprintf(errstr,"Could not read in tissue dVPlus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue dVPlus is %f\n",presc->array[i].dVPlus);

		// pop off the vPlus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,vPlus_line))
		{
			sprintf(errstr,"Missing vPlus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue vPlus
		if (fscanf(fid,"%f\n",&presc->array[i].vPlus) != 1)
		{
			sprintf(errstr,"Could not read in tissue vPlus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue vPlus is %f\n",presc->array[i].vPlus);

		// pop off the betaVMinus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,betaVMinus_line))
		{
			sprintf(errstr,"Missing betaVMinus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}

		// get the tissue beta_minus
		if (fscanf(fid,"%f\n",&presc->array[i].betaVMinus) != 1)
		{
			sprintf(errstr,"Could not read in tissue betaVMinus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue betaVMinus is %f\n",presc->array[i].betaVMinus);

		// pop off the dVMinus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,dVMinus_line))
		{
			sprintf(errstr,"Missing dVMinus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue dVMinus
		if (fscanf(fid,"%f\n",&presc->array[i].dVMinus) != 1)
		{
			sprintf(errstr,"Could not read in tissue dVMinus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue dVMinus is %f\n",presc->array[i].dVMinus);

		// pop off the vMinus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,vMinus_line))
		{
			sprintf(errstr,"Missing vMinus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue vMinus
		if (fscanf(fid,"%f\n",&presc->array[i].vMinus) != 1)
		{
			sprintf(errstr,"Could not read in tissue vMinus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue vMinus is %f\n",presc->array[i].vMinus);

		// pop off the betaPlus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,betaPlus_line))
		{
			sprintf(errstr,"Missing betaPlus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue betaPlus
		if (fscanf(fid,"%f\n",&presc->array[i].betaPlus) != 1)
		{
			sprintf(errstr,"Could not read in tissue betaPlus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue betaPlus is %f\n",presc->array[i].betaPlus);

		// pop off the dosePlusFileName
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,dosePlusFilename_line))
		{
			sprintf(errstr,"Missing dosePlusfilename line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the dose plus prescription filename
		if (fscanf(fid,"%s\n",presc->array[i].dosePlusFilename) != 1)
		{
			sprintf(errstr,"Could not read in dosePlusFilename for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		// read in the dose prescription file, which is stored in the same format as an individual beamlet:
		// Xcount Ycount Zcount Nind non_zero_indices non_zero_values
		//  int    int    int    int    int array       single array
		
		// open the overdose penalty file
		if ((dosePlusFile = fopen(presc->array[i].dosePlusFilename,"rb")) == NULL)
		{
			sprintf(errstr,"Could not read in the file %s.\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}
		
		// read-in the first three indices, which are the int values (Xcount, Ycount, Zcount)
		if (fread(dims,sizeof(int),3,dosePlusFile) != 3)
		{
			sprintf(errstr,"Could not read in Xcount, Ycount, Zcount data from %s\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}
		
		if (i == 0)
		{
			// record tissue mask dimensions
			presc->x_count = dims[0];
			presc->y_count = dims[1];
			presc->z_count = dims[2];
			presc->Nvox = presc->x_count*presc->y_count*presc->z_count;
		}  // ensure that all tissue masks have the same dimensions
		else
			if ((int)dims[0] != presc->x_count || dims[1] != presc->y_count || dims[2] != presc->z_count)
			{
				sprintf(errstr,"Dose prescription file %d of %d has dimensions of (%d,%d,%d), when they should be (%d,%d,%d).\n",
					i+1,presc->Ntissue,dims[0],dims[1],dims[2],
					presc->x_count,presc->y_count,presc->z_count);
				return(FAILURE);
			}

		// read in the number of non-zero indices
		if (fread(&Nind,sizeof(int),1,dosePlusFile) != 1)
		{
			sprintf(errstr,"Could not read in Nind data from %s\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}
		else
			presc->array[i].Nind = Nind;
		
		// allocate space for the indices
		if ((presc->array[i].ind = (int *)malloc(sizeof(int)*Nind)) == NULL)
		{
			sprintf(errstr,"Failed to allocate space for mask indices from %s\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}

		// allocate space for the dose prescription data
		if ((presc->array[i].dPlus = (float *)malloc(sizeof(int)*Nind)) == NULL)
		{
			sprintf(errstr,"Failed to allocate space for data indices from %s\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}

		// read in the non-zero indices from the mask file
		if (fread(presc->array[i].ind,sizeof(int),Nind,dosePlusFile) != Nind)
		{
			sprintf(errstr,"Failed to read in all mask indices from %s\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}

		// read in the non-zero data from the mask file
		if (fread(presc->array[i].dPlus,sizeof(float),Nind,dosePlusFile) != Nind)
		{
			sprintf(errstr,"Failed to read in all non-zero data from %s\n",presc->array[i].dosePlusFilename);
			return(FAILURE);
		}

		fclose(dosePlusFile);

		// pop off the betaMinus line
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,betaMinus_line))
		{
			sprintf(errstr,"Missing betaMinus line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the tissue betaMinus
		if (fscanf(fid,"%f\n",&presc->array[i].betaMinus) != 1)
		{
			sprintf(errstr,"Could not read in tissue betaMinus for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		printf("Tissue betaMinus is %f\n",presc->array[i].betaMinus);

		// pop off the doseMinusFileName
		if (fgets(str,100,fid) == NULL)
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}
		else if (!strcmp(str,doseMinusFilename_line))
		{
			sprintf(errstr,"Missing doseMinusfilename line from tissue %d in file %s\n",i,presc_filename);
			return(FAILURE);
		}
		
		// get the dose plus prescription filename
		if (fscanf(fid,"%s\n",presc->array[i].doseMinusFilename) != 1)
		{
			sprintf(errstr,"Could not read in doseMinusFilename for tissue %d in file %s.\n",i,presc_filename);
			return(FAILURE);
		}

		// read in the dose prescription file, which is stored in the same format as an individual beamlet:
		// Xcount Ycount Zcount Nind non_zero_indices non_zero_values
		//  int    int    int    int    int array       single array
		
		// open the underdose penalty file
		if ((doseMinusFile = fopen(presc->array[i].doseMinusFilename,"rb")) == NULL)
		{
			sprintf(errstr,"Could not read in the file %s.\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}
		
		// read-in the first three indices, which are the int values (Xcount, Ycount, Zcount)
		if (fread(dims,sizeof(int),3,doseMinusFile) != 3)
		{
			sprintf(errstr,"Could not read in Xcount, Ycount, Zcount data from %s\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}
		
		if (i == 0)
		{
			// record tissue mask dimensions
			presc->x_count = dims[0];
			presc->y_count = dims[1];
			presc->z_count = dims[2];
		}  // ensure that all tissue masks have the same dimensions
		else
			if ((int)dims[0] != presc->x_count || dims[1] != presc->y_count || dims[2] != presc->z_count)
			{
				sprintf(errstr,"Dose prescription file %d of %d has dimensions of (%d,%d,%d), when they should be (%d,%d,%d).\n",
					i+1,presc->Ntissue,dims[0],dims[1],dims[2],
					presc->x_count,presc->y_count,presc->z_count);
				return(FAILURE);
			}

		// read in the number of non-zero indices
		if (fread(&Nind,sizeof(int),1,doseMinusFile) != 1)
		{
			sprintf(errstr,"Could not read in Nind data from %s\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}
		else
			presc->array[i].Nind = Nind;
		
		// allocate space for the indices
		if ((presc->array[i].ind = (int *)malloc(sizeof(int)*Nind)) == NULL)
		{
			sprintf(errstr,"Failed to allocate space for mask indices from %s\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}

		// allocate space for the dose prescription data
		if ((presc->array[i].dMinus = (float *)malloc(sizeof(int)*Nind)) == NULL)
		{
			sprintf(errstr,"Failed to allocate space for data indices from %s\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}

		// read in the non-zero indices from the mask file
		if (fread(presc->array[i].ind,sizeof(int),Nind,doseMinusFile) != Nind)
		{
			sprintf(errstr,"Failed to read in all mask indices from %s\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}

		// read in the non-zero data from the mask file
		if (fread(presc->array[i].dMinus,sizeof(float),Nind,doseMinusFile) != Nind)
		{
			sprintf(errstr,"Failed to read in all non-zero data from %s\n",presc->array[i].doseMinusFilename);
			return(FAILURE);
		}

		fclose(doseMinusFile);

		// diagnostic lines for ensuring that indices were properly read in

		/* printf("successfully read-in indices\n");
		sprintf(str2,"tissue%d.bin",i);
		printf("str2 = %s\n",str2);
		test_file = fopen(str2,"wb"); 
		fwrite(presc->array[i].ind,sizeof(int),presc->array[i].Nind,test_file);
		fclose(test_file); */
	}

	fclose(fid);

 return SUCCESS;
}

int read_in_beamlets(GRID_ARRAY *beamlets, char *beamlet_header, char *init_beam_weights_file, 
					 int *usedVoxels, OPT_PARAM *op, GRID *doseInit)
// This function is lazier than read_in_geometry since it does not check to ensure that
// all of the required lines in the beamlet_header file are present.  This is something
// that can be added later, as it is not necessary for the function of the program.
{
 int k,j,m,n,Nbeamlets,num,Nind,beamletNum;
 int dims[3];
 int *beamlet_tally;  // tally to ensure that all beamlets are read-in
 int Nuv;
 int *ind;
 float *data;

 char bixbatchfile[200];
 char str[200];
 char bix_dir[200];  // directory containing beamlet batches
 char bix_batch_base_filename[200];
 char bix_batch_extension[200];

 FILE *fid;

	// read in beamlet header file
	if( (fid = fopen(beamlet_header,"rb")) == NULL)
	{
		sprintf(errstr,"Could not open beamlet header file.");
		return(FAILURE);
	}
	else
	{
		// pop off the next line of the header file
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}

		// read in the dimensions of the beamlets
		if (fscanf(fid,"%d %d %d\n",&beamlets->x_count,&beamlets->y_count,&beamlets->z_count) != 3)
		{
			sprintf(errstr,"Could not read-in the beamlet dimensions.\n");
			return(FAILURE);
		}
		else
			beamlets->Nvox = beamlets->x_count*beamlets->y_count*beamlets->z_count;

		// pop off the next line of the header file
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}

		// read in the number of beamlets
		if (fscanf(fid,"%d\n",&beamlets->num_beamlets) != 1)
		{
			sprintf(errstr,"Could not read-in the number of beamlets.\n");
			return(FAILURE);
		}

		// pop off the next line of the header file
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}

		// read in the number of beamlet batches
		if (fscanf(fid,"%d\n",&beamlets->num_beamlet_batches) != 1)
		{	
			sprintf(errstr,"Could not read-in the number of beamlet batches.\n");
			return(FAILURE);
		} 

		// pop off the next line of the header file
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}

		// read in the name of the beamlet directory
		if (fscanf(fid,"%s\n",bix_dir) != 1)
		{	
			sprintf(errstr,"Could not read-in the name of the beamlet directory.\n");
			return(FAILURE);
		} 

		strcpy(op->bix_dir,bix_dir); // save bixel directory

		// pop off the next line of the header file
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}

		// read in the name of the beamlet directory
		if (fscanf(fid,"%s\n",bix_batch_base_filename) != 1)
		{	
			sprintf(errstr,"Could not read-in the base filename of the beamlet batches.\n");
			return(FAILURE);
		} 

		strcpy(op->bix_batch_base_filename,bix_batch_base_filename); // save bixel base filename

		// pop off the next line of the header file
		if (fgets(str,100,fid) == NULL) 
		{
			sprintf(errstr,"Could not read from from header file.");
			return(FAILURE);
		}

		// read in the beamlet batch extension name
		if (fscanf(fid,"%s\n",bix_batch_extension) != 1)
		{	
			sprintf(errstr,"Could not read-in the beamlet batch extension name.\n");
			return(FAILURE);
		}

		strcpy(op->bix_batch_extension,bix_batch_extension); // save bixel extension
	}

 // read-in beamlet data
 if((beamlet_tally = (int *)malloc(beamlets->num_beamlets*sizeof(int))) == NULL)
 {
	sprintf(errstr,"Unable to allocate memory for beamlet_tally.");
	return(FAILURE);
 }

 // allocate space for sparse beamlets
 if ((beamlets->array = (SPARSE_ARRAY *)malloc(beamlets->num_beamlets*sizeof(SPARSE_ARRAY))) == NULL)
 {
	sprintf(errstr,"Unable to allocate memory for sparse beamlet matrix.\n");
	return(FAILURE);
 } 
 
  // allocate space for the other arrays in beamlets
  if ((beamlets->beam_weight = (float *)malloc(beamlets->num_beamlets*sizeof(float))) == NULL)
  {
	sprintf(errstr,"Unable to allocate memory for beam_weight array in beamlet structure.\n");
	return(FAILURE);
  } 

  if ((beamlets->initial_beam_weight = (float *)malloc(beamlets->num_beamlets*sizeof(float))) == NULL)
  {
	sprintf(errstr,"Unable to allocate memory for initial_beam_weight array in beamlet structure.\n");
	return(FAILURE);
  } 

 // load the initial beamlet weights
 if((fid = fopen(init_beam_weights_file,"rb")) == NULL)
 {
	 sprintf(errstr,"Failed to open initial beam intensity file.\n");
	 return(FAILURE);
 }

 if((int)fread(beamlets->initial_beam_weight,sizeof(float),beamlets->num_beamlets,fid) != beamlets->num_beamlets)
 {
	 sprintf(errstr,"Failed to read in initial beamlet intensities.\n");
	 return(FAILURE);
 }
 else
	 printf("Successfully read in initial beamlet weights.\n");

 // copy initial beam weights to the variable beam weights
 for (k=0;k<beamlets->num_beamlets;k++) 
	 beamlets->beam_weight[k] = beamlets->initial_beam_weight[k];

 // truncate beamlet weights that are less than zero
 for (k=0;k<beamlets->num_beamlets;k++) 
	 if (beamlets->beam_weight[k] < 0)
		 beamlets->beam_weight[k] = 0.0;

 // truncate the beamlet weights that violate the modulation requirements
 truncateBeamWeights(beamlets->beam_weight,beamlets->num_beamlets, op->modFactor);

 // start tally at all zeros
 for (k=0;k<beamlets->num_beamlets;k++) beamlet_tally[k] = 0;

 // start initial dose grid with all zeros
 for (k=0;k<beamlets->Nvox;k++)
	 doseInit->matrix[k] = 0.0;
 
 beamletNum = 0;

 for (k=0;k<beamlets->num_beamlet_batches;k++)
 {
	// open the beamlet batch file
	sprintf(bixbatchfile,"%s/%s%d.%s",bix_dir,bix_batch_base_filename,k,bix_batch_extension);
	if ((fid = fopen(bixbatchfile,"rb")) == NULL)
	{
		sprintf(errstr,"Unable to open file for beamlet batch %d.\n",k);
		return(FAILURE);
	}

	// read-in the number of beamlets that are in this file
    if((int)fread(&Nbeamlets,sizeof(int),1,fid) != 1)
	{
		sprintf(errstr,"Unable to read number of beamlets in file:\n%s",bixbatchfile);
		return(FAILURE);
	}

	// read-in all of the beamlets in the current batch
	for (j=0;j<Nbeamlets;j++)
	{
		if((int)fread(&num,sizeof(int),1,fid) != 1)
		{
			sprintf(errstr,"Unable to read beamlet number for beamlet %d in file:\n%s",j,bixbatchfile);
			return(FAILURE);
		}

		// read the beamlet dimensions
		if((int)fread(dims,sizeof(int),3,fid) != 3)
		{
			sprintf(errstr,"Unable to read beamlet dimensions for beamlet %d in file:\n%s",j,bixbatchfile);
			return(FAILURE);
		}

		// ensure that the beamlet has the same grid size as the prescription
		if((dims[0] != beamlets->x_count) || (dims[1] != beamlets->y_count) || (dims[2] != beamlets->z_count))
		{
			sprintf(errstr,"Mismatch in grid dimensions for beamlet %d in file:\n%s",j,bixbatchfile);
			return(FAILURE);
		}

		// read the number of non-zero values in the current beamlet
		if((int)fread(&Nind,sizeof(int),1,fid) != 1)
		{
			sprintf(errstr,"Unable to read number of non-zero values for beamlet %d in file:\n%s",j,bixbatchfile);
			return(FAILURE);
		}

		// allocate memory and read-in indices for the current beamlet
		if((ind = (int *)malloc(sizeof(int)*Nind)) == NULL)
		{
			sprintf(errstr,"Unable to allocate memory for indices of beamlet %d\n",num);
			return(FAILURE);
		}
		if((int)fread(ind,sizeof(int),Nind,fid) != Nind)
		{
			sprintf(errstr,"Unable to read-in indices for beamlet %d\n",num);
			return(FAILURE);
		}

		// allocate memory and read-in indices for the current beamlet
		if((data = (float *)malloc(sizeof(float)*Nind)) == NULL)
		{
			sprintf(errstr,"Unable to allocate memory for data of beamlet %d\n",num);
			return(FAILURE);
		}	
		if((int)fread(data,sizeof(float),Nind,fid) != Nind)
		{
			sprintf(errstr,"Unable to read-in dose data for beamlet %d\n",num);
			return(FAILURE);
		}

		if (beamletNum < 0 || beamletNum >= beamlets->num_beamlets)
		{
			sprintf(errstr,"%d is an invalid beamlet number for beamlet %d in batch %d.",num,j,k);
			return(FAILURE);
		}
		else
			beamlet_tally[beamletNum] = 1; // tally that beamlet "beamletNum" was just read in

		// add the beamlet to the initial dose distribution
		for (m=0;m<Nind;m++)
			doseInit->matrix[ind[m]] += beamlets->beam_weight[beamletNum]*data[m];

		// Count number of voxels in the current beamlet that will actually be used
		// in the optimization.
		Nuv = 0;
		for (m=0;m<Nind;m++)
			if (usedVoxels[ind[m]] == 1)
				Nuv++;

		beamlets->array[beamletNum].size = Nuv;

		// allocate space for the cleaned-beamlets
		if((beamlets->array[beamletNum].indices = (int *)malloc(sizeof(int)*Nuv)) == NULL)
		{
			sprintf(errstr,"Unable to allocate memory for indices of beamlet %d\n",num);
			return(FAILURE);
		}

		// allocate memory and read-in indices for the current beamlet
		if((beamlets->array[beamletNum].data = (float *)malloc(sizeof(float)*Nuv)) == NULL)
		{
			sprintf(errstr,"Unable to allocate memory for data of beamlet %d\n",num);
			return(FAILURE);
		}

		// pull out only the voxels that are actually used
		n = 0;
		for (m=0;m<Nind;m++)
			if (usedVoxels[ind[m]] == 1)
			{
				beamlets->array[beamletNum].indices[n] = ind[m];
				beamlets->array[beamletNum].data[n] = data[m];
				n++;
			}

		// free indices and data for the full beamlet
		free(ind);
		free(data);

		beamletNum++;  // increment the beamlet number
	}
	printf("Successfully loaded beamlet batch %d of %d.\n",k+1,beamlets->num_beamlet_batches);

	fclose(fid);

}

  
  // check beamlet tally to ensure that all beamlets were read in
  for (j=0;j<beamlets->num_beamlet_batches;j++)
	  if (beamlet_tally[j] == 0)
	  {
		 sprintf(errstr,"Missing beamlet %d.",j);
		 return(FAILURE);
	  }

  return(SUCCESS);
}