#include "TH1.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"

//TF1* get_mem_func(Int_t a, Int_t b, Int_t (*c)[15]), Double_t *d, Double_t *e);

void FLBM()
{
  gROOT->Reset();
  gROOT->SetStyle("Plain");
 
  const Int_t num_stages = 7;
 
  Double_t barrier_mean[] = {0.038, 0.044, 0.05, 0.07, 0.088, 0.111, 0.14, 0.176, 0.222, 0.279, 0.352, 0.443, 0.559, 0.704, 0.887};
  //I-, I, I+, L-, L, L+, M-, M, M+, H-, H, H+, VH-, VH, VH+
 
  //Double_t barrier_weight[] = {0.1832, 0.1832, 0.0611, 0.0916, 0.0366 , 0.1111, 0.0278, 0.2222, 0.0278, 0.0278, 0.0278};
  Double_t barWeight_cost[] = {130, 55, 5, 0.15, 0.14, 13.3, 1, 55, 1, 3.2, 1}; //Cost-based weights!
  Double_t barWeight_cost2[] = {0.5804, 0.1319, 0.0246, 0.0005, 0.0005, 0.0468, 0.0035, 0.1935, 0.0035, 0.0113, 0.0035};
  Double_t barWeight_ex1[] = {0.4317, 0.0457, 0.2137, 0.0944, 0.0188, 0.0393, 0.0077, 0.0427, 0.0274, 0.0910, 0.0056}; //Expert #1
  Double_t barWeight_ex3[] = {0.1832, 0.1832, 0.0611, 0.0916, 0.0366, 0.1111, 0.0278, 0.2222, 0.0278, 0.0278, 0.0278}; //Expert #3
  Double_t stage_weight[] = {2.21E-3, 1.177474403, 0.675927252, 0.440789474, 2.644396552, 8.814655172, 44.07327586};  //Stage weights, PWR once-through

  char* FN_LUT[] = {"I-", "I", "I+", "L-", "L", "L+", "M-", "M", "M+", "H-", "H", "H+", "VH-", "VH", "VH+"};

  Double_t PWR_cycle_params[7][15] = { 
    {337487.29, 0, 12.60117, 8.87731E-06, 0.103864232, 200, 0.001, 3.27E-09, 3, 100, 365, 0.562243396, 0.01, 10, 2},  
    {337487.29, 0.72, 12.60117, 8.87731E-06, 0.103864232, 200, 0.001, 1.74E-06, 3, 100, 365, 0.562243396, 0.01, 10, 2},
    {18864.78162, 3.5, 12.2581875, 1.02145E-05, 0.44817585, 200, 0.001, 3.58E-05, 3, 0.1, 180, 10.00806708, 0.01, 0.1, 1},
    {18864.78162, 3.5, 12.2581875, 1.02145E-05, 0.44817585, 200, 0.001, 4.67E-05, 3, 100, 365, 1.314088893, 0.1, 10, 0.5},
    {10, 47.24578648, 3326.837727, 0.122650464, 212.2926844, 3000, 100, 8.81E-02, 5, 1, 10, 24.54, 0, 2, 3},
    {10, 48.74887164, 3333.385524, 0.119703203, 205.3272656, 3000, 100, 8.81E-02, 5, 100, 180, 24.54, 0, 10, 10},
    {10, 50.91386492, 3230.72987, 9.19E-02, 153.6471609, 3000, 100, 8.81E-02, 5, 100, 100, 24.54, 0, 10, 50}
  };


  TF1* levelFunctions[7];
  Double_t stage_FN = 0;  
  Int_t PWR_barrier_level[7][15];
  Int_t num_barriers = 11;
  Double_t iso_level = 0;
 
  Double_t totEx1Weight = getWeight(barWeight_ex1, num_barriers);
  Double_t totEx3Weight = getWeight(barWeight_ex3, num_barriers);
  Double_t totCostWeight = getWeight(barWeight_cost, num_barriers);
  Double_t totCostWeight2 = getWeight(barWeight_cost2, num_barriers);

  Double_t total_stage_weight =  getWeight(stage_weight, num_stages);
  Double_t* barrier_weight = barWeight_cost; //Use the cost-based barriers for stage FN's 
  Double_t total_barrier_weight = totCostWeight; //Use cost-based for stage FN's
  

  for(Int_t i=0; i < num_stages; i++) {
    total_stage_weight += stage_weight[i];

    stage_FN = 0;
    iso_level = 0;

    //Idea for expansion - fold these function calls into an array - loop through an array of function calls!
    PWR_barrier_level[i][0] = bl_iso_cm(PWR_cycle_params[i][0]);
    PWR_barrier_level[i][1] = bl_iso_enrich(PWR_cycle_params[i][1]);
    PWR_barrier_level[i][2] = bl_iso_sfn(PWR_cycle_params[i][2]);
    PWR_barrier_level[i][3] = bl_iso_heat(PWR_cycle_params[i][3]);
    PWR_barrier_level[i][4] = bl_iso_gamma(PWR_cycle_params[i][4]);
    PWR_barrier_level[i][5] = bl_chem_sep(PWR_cycle_params[i][5]);
    PWR_barrier_level[i][6] = bl_rad_dose(PWR_cycle_params[i][6]);
    PWR_barrier_level[i][7] = bl_mass_bulk(PWR_cycle_params[i][7]);
    PWR_barrier_level[i][8] = bl_detectability(PWR_cycle_params[i][8]);
    PWR_barrier_level[i][9] = bl_fac_unattr(PWR_cycle_params[i][9]);
    PWR_barrier_level[i][10] = bl_fac_access(PWR_cycle_params[i][10]);
    PWR_barrier_level[i][11] = bl_avail_CM(PWR_cycle_params[i][11]);
    PWR_barrier_level[i][12] = bl_uncertainty(PWR_cycle_params[i][12]);
    PWR_barrier_level[i][13] = bl_skills_exp_knowl(PWR_cycle_params[i][13]);
    PWR_barrier_level[i][14] = bl_res_time(PWR_cycle_params[i][14]);
    
    //Need to calculate isotopic value first, then fold into below
    
    for (Int_t j=0;j<5;j++) {
      //Assign the max level of the isotopic barriers
      if(barrier_mean[PWR_barrier_level[i][j]] > iso_level) { 
	iso_level = barrier_mean[PWR_barrier_level[i][j]];
      }
    }
   
    stage_FN = barrier_weight[0]*iso_level;
  
    for (Int_t j = 1; j <  num_barriers; j++) {      
      stage_FN += barrier_weight[j]*barrier_mean[PWR_barrier_level[i][j+4]];
    }
    //cout << "Tot_bar_weight = " << total_bar_weight << endl;
    stage_FN = stage_FN/total_barrier_weight;
    cout << "Stage " << i+1 << ": Stage_FN = " << stage_FN << endl;      
  }

 
  //Assign gaussian functions to barrier levels

  Int_t* barrierPointers[7];
  for(Int_t i = 0; i < num_stages; i++) {
    barrierPointers[i] = PWR_barrier_level[i];
  }
  
  c1 = new TCanvas("c1","System Fuzzy Number Membership Function",200,10,800,600);
  pad1 = new TPad("pad1","LWR-OT Pad",0,0,1,1,-1,0,0);
  pad1->Draw();
  pad1->cd();    

  TF1* systemMember_cost = get_mem_func(num_barriers, num_stages, totCostWeight, total_stage_weight, barrierPointers, barWeight_cost, stage_weight, barrier_mean);
  //systemMember_cost->SetTitle("System Membership: Cost Weighting (Access)");
  systemMember_cost->SetTitle("LWR-OT System Membership: Cost (Access, Paper), Expert (#1 and #3) Weighting");  
  systemMember_cost->SetLineColor(kOrange);
  systemMember_cost->DrawCopy();

  TF1* systemMember_cost2 = get_mem_func(num_barriers, num_stages, totCostWeight2, total_stage_weight, barrierPointers, barWeight_cost2, stage_weight, barrier_mean);
  systemMember_cost2->SetTitle("System Membership: Cost Weighting (FLB Paper)"); 
  systemMember_cost2->SetLineColor(kRed);
  systemMember_cost2->DrawCopy("SAME");

  TF1* systemMember_ex1 = get_mem_func(num_barriers, num_stages, totEx1Weight, total_stage_weight, barrierPointers, barWeight_ex1, stage_weight, barrier_mean);
  systemMember_ex1->SetTitle("System Membership: Expert #1 Weighting");
  systemMember_ex1->SetLineColor(kBlue);
  systemMember_ex1->DrawCopy("SAME");
  
  TF1* systemMember_ex3 = get_mem_func(num_barriers, num_stages, totEx3Weight, total_stage_weight, barrierPointers, barWeight_ex3, stage_weight, barrier_mean);
  systemMember_ex3->SetTitle("System Membership: Expert #3 Weighting");
  systemMember_ex3->SetLineColor(kGreen);
  systemMember_ex3->DrawCopy("SAME");

  TLegend *legend = new TLegend(0.75,0.8,1,0.95,"");
  legend->AddEntry( systemMember_cost, "Cost-based (ACCESS)","l");
  legend->AddEntry( systemMember_cost2, "Cost-based (FLB Paper)","l");
  legend->AddEntry( systemMember_ex1, "Expert #1","l");
  legend->AddEntry( systemMember_ex3, "Expert #3","l"); 
  
  legend->Draw();

  //pad1->Print("C:\\Documents and Settings\\Steve Skutnik\\My Documents\\Research\\imgs\\pwr-ot.png");
  //pad1->Print("C:\\Documents and Settings\\Steve Skutnik\\My Documents\\Research\\imgs\\pwr-ot.pdf");
  
 //cout << "Eval(0.5) = " << myStageMember->Eval(0.5) << endl;
 


} //NEXT: Find centroid mean of system-level function for defuzzification

Double_t getWeight(Double_t *weightArray, Int_t length) 
{
  Double_t total_weight = 0;
  for(Int_t i=0; i < length; i++) {
    total_weight += weightArray[i];
  }
  return total_weight;
}

Int_t bl_iso_cm(Double_t CM) 
{
  Int_t barrier_level = 0;

  if(CM >= 100) {
    if(CM >= 150) {
      if(CM >= 200) {
	if(CM >= 300) {
	  if(CM >= 500) {
	    if(CM >= 800) {
	      barrier_level = 13;
	    }
	    else barrier_level = 10;
	  }
	  else barrier_level = 7;
	}
	else barrier_level = 4;
      }
      else barrier_level = 2;
    }
    else barrier_level = 1;
  } 
  else barrier_level = 0;

  return barrier_level;
}

Int_t bl_iso_sfn(Double_t SFN) 
{
  Int_t barrier_level = 0;

  if(SFN >= 1E+2) {
    if(SFN >= 1E+3) {
      if(SFN >= 1E+4) {
	if(SFN >= 1E+5) {
	  if (SFN >= 5E+5) {
	    if (SFN >= 1E+6) {
	      barrier_level = 10;
	    }
	    else barrier_level = 7;
	  }
	  else barrier_level = 4;
	}
	else barrier_level = 3;
      }
      else barrier_level = 2;
    }
    else barrier_level = 1;
  }
  else barrier_level = 0;
  
  return barrier_level;
}

Int_t bl_iso_enrich(Double_t enrich) 
{
  Int_t barrier_level = 0;

  if(enrich <= 80) {
    if(enrich <= 50) {
      if(enrich <= 45) {
	if(enrich <= 40) {
	  if(enrich <= 30) {
	    if(enrich <= 20) {
	      if(enrich <= 10) {
		if(enrich <= 5) {
		  if(enrich <= 1) {
		    barrier_level = 13; //VH
		  }
		  else barrier_level = 10; //H
		}
		else barrier_level = 9; //H-
	      }
	      else barrier_level = 7; //M
	    }
	    else barrier_level = 6; //M-
	  }
	  else barrier_level = 5; //L+
	}	 
	else barrier_level = 4; //L
      }
      else barrier_level = 3; //L-
    }
    else barrier_level = 2; //I+
  }
  else barrier_level = 0; //I- 

  return barrier_level;
}

Int_t bl_iso_heat(Double_t hr) 
{
  Int_t barrier_level = 0;

  if(hr > 1) {
    if(hr > 10) {
      if(hr > 100) {
	if(hr > 500) {
	  barrier_level = 10; //H
	} 
	else barrier_level = 7; //M
      }
      else barrier_level = 4; //L
    }
    else barrier_level = 3; //L- 
  }
  else barrier_level = 0; //I-

  return barrier_level;
}

Int_t bl_iso_gamma(Double_t gamma) 
{
  Int_t barrier_level = 0;

  if(gamma> 1E+3) {
    if(gamma >= 1E+4) {
      if(gamma >= 1E+5) {
	if(gamma >= 1E+6) {
	  barrier_level = 10; //H
	}
	else barrier_level = 7; //M
      } 
      else barrier_level = 4; //L
    }
    else barrier_level = 1;
  }
  else barrier_level = 0;

  return barrier_level;
}

Int_t bl_chem_sep(Double_t sc) 
{

  Int_t barrier_level = 0;
  
  if(sc >= 10) {
    if(sc >= 50) {
      if(sc >= 150) {
	if(sc >= 250) {
	  if(sc >= 400) {
	    if(sc >= 600) {
	      if(sc >= 800) {
		if(sc >= 1000) {
		  if(sc >= 1500) {
		    if(sc >= 2000) {
		      if(sc >= 2500) {
			if(sc >= 3000) {
			  if(sc >= 5000) {
			    barrier_level = 14; //VH+
			  }
			  else barrier_level = 13; //VH
			}
			else barrier_level = 12; //VH-
		      }
		      else barrier_level = 11; //H+
		    }
		    else barrier_level = 10; //H
		  }
		  else barrier_level = 9; //H-
		}
		else barrier_level = 8; //M+
	      }
	      else barrier_level = 7; //M
	    }
	    else barrier_level = 6; //M-
	  }
	  else barrier_level = 5; //L+
	}	 
	else barrier_level = 4; //L
      }
      else barrier_level = 3; //L-
    }
    else barrier_level = 2; //I+
  }
  else barrier_level = 1; //I

  return barrier_level;
}

Int_t bl_rad_dose(Double_t dose_rate) 
{
  Int_t barrier_level = 0;

  if(dose_rate >= 1E-5) {
    if(dose_rate >= 5E-5) {
      if(dose_rate >= 1) {
	if(dose_rate >= 10) {
	  barrier_level = 13; //VH
	}
	else barrier_level = 10; //H
      }
      else barrier_level = 7; //M
    }
    else barrier_level = 4; //L
  }
  else barrier_level = 1; //I

  return barrier_level;
}

Int_t bl_mass_bulk(Double_t conc) 
{
  Int_t barrier_level = 0;
  
  if(conc < 0.1) {
    if(conc < 0.02) {
      if(conc < 0.01) {
	if(conc < 0.005) {
	  barrier_level = 13; //VH
	}
	else barrier_level = 10; //H 
      }
      else barrier_level = 7; //M
    }
    else barrier_level = 4; //L
  }
  else barrier_level = 1; //I

  return barrier_level;
}

Int_t bl_fac_unattr(Double_t mod_time) 
{
  Int_t barrier_level = 0;

  if (mod_time >= 0.1) {
    if(mod_time >= 1) {
      if(mod_time >= 5) {
	if(mod_time >= 30) {
	  if(mod_time >= 100) {
	    barrier_level = 14; //VH+
	  }
	  else barrier_level = 13; //VH
	}
	else barrier_level = 10; //H
      }
      else barrier_level = 7; //M
    }
    else barrier_level = 4; //L 
  }
  else barrier_level = 1; //I
  
  return barrier_level;
}


Int_t bl_detectability(Int_t detect) {
  Int_t barrier_level = 0;
  switch(detect) {
  case 1: 
    barrier_level = 1; //I
    break;
  case 2:
    barrier_level = 4; //L
    break;
  case 3:
    barrier_level = 7; //M
    break;
  case 4: 
    barrier_level = 10; //H
    break;
  case 5: 
    barrier_level = 13; //VH
    break;
  default: 
    barrier_level = 1;
    break;
  }
  return barrier_level;
}
  
Int_t bl_fac_access(Double_t acc_time) 
{
  Int_t barrier_level = 0;

  if(acc_time < 365) {
    if(acc_time < 243) {
      if (acc_time < 180) {
	  if(acc_time < 100) {
	    if(acc_time < 10) {
	      barrier_level = 14; //VH+
	    }
	    else barrier_level = 13; //VH
	  }
	  else barrier_level = 10; //H
	}
	else barrier_level = 7; //M
      }
      else barrier_level = 4; //L
  }
  else barrier_level = 1; //I

  return barrier_level;
}

Int_t bl_avail_CM(Double_t CM)
{

  Int_t barrier_level = 0;
  
  if(CM < 10) {
    if(CM < 1) {
      if(CM < 0.1) {
	if(CM < 0.01) {
	  barrier_level = 13; //VH
	}
	else barrier_level = 10; //H
      }
      else barrier_level = 7; //M
    }
    else barrier_level = 4; //L
  }
  else barrier_level = 1; //I

  return barrier_level;
}

Int_t bl_uncertainty(Double_t CM_uncert) 
{
  Int_t barrier_level = 0;

  if(CM_uncert < 0.96) {
    if(CM_uncert < 0.82) {
      if(CM_uncert < 0.65) {
	if(CM_uncert < 0.29) {
	  barrier_level = 14; //VH+
	}
	else barrier_level = 13; //VH
      }
      else barrier_level = 10; //H
    }
    else barrier_level = 7; //M
  }
  else barrier_level = 1; //I

  return barrier_level;
}

Int_t bl_skills_exp_knowl(Double_t time)
{
  Int_t barrier_level = 0;
  
  if(time >= 0.1) {
    if(time >= 0.5) {
      if(time >= 2) {
	if(time >= 10) {
	  barrier_level = 13; //VH
	}
	else barrier_level = 10; //H
      }
      else barrier_level = 7; //M
    }
    else barrier_level = 4; //L
  }
  else barrier_level = 1;

  return barrier_level;
}
  
Int_t bl_res_time(Double_t time) {
  Int_t barrier_level = 0;
  
  if(time < 10) {
    if(time < 5) {
      if(time < 2) {
	if(time < 0.5) {
	  barrier_level = 13; //VH
	}
	else barrier_level = 10; //H
      }
      else barrier_level = 7; //M
    }
    else barrier_level = 4; //L
  }
  else barrier_level = 1; //I
 
  return barrier_level;
}

TF1* get_mem_func(Int_t num_barriers, Int_t num_stages, Double_t total_bar_weight, Double_t total_stage_weight, Int_t **barriers, Double_t *barrier_weight, Double_t *stage_weight, Double_t *mean_value)
{
  Int_t num_params = 4+(num_barriers)*(num_stages+1) + num_stages;
  //cout << "Num_params = " << num_params << endl;
  TF1* smember = new TF1(membership_func,0,1,num_params);
  smember->SetParameter(0,num_barriers);
  smember->SetParameter(1,num_stages);
  smember->SetParameter(2,total_bar_weight);
  smember->SetParameter(3,total_stage_weight);

  Int_t par_iter = 4;
  
  //Set barrier weights
  for(Int_t i=0; i < num_barriers; i++) {
    smember->SetParameter(par_iter+i, barrier_weight[i]);
  }
  
  for(Int_t k = 0; k < num_stages; k++) {      
    par_iter = 4 + num_barriers;
    
    //Set stage weights
    smember->SetParameter(par_iter+k, stage_weight[k]);
    
    par_iter += num_stages; //Skip to next block
    
    //Assign the max level of the isotopic barriers
    Double_t iso_level = 0;      
    for (Int_t j=0;j<5;j++) {    
      if( mean_value[barriers[k][j]] > iso_level) { 
	iso_level = mean_value[barriers[k][j]];
      }
    }
    
    par_iter += k*num_barriers;
    // cout << "Iso par_iter = " << par_iter << endl;
    smember->SetParameter(par_iter, iso_level);

    for(Int_t i = 1; i < num_barriers; i++) {
      smember->SetParameter(par_iter+i, mean_value[barriers[k][i+4]]);     
    }
  }
  return smember;
}


//Need to figure out why this explodes if # of stages > 1...
Double_t membership_func(Double_t *x, Double_t *par) 
{
  
  Double_t val = 0;
  Double_t tval = 0; //Temp val

  Double_t xx  = x[0];

  //Structure of params:
  //params[0] = num_barriers 
  //params[1] = num_stages
  //params[2] = total_bar_weight
  //params[3] = total_stage_weight
  //params[4+num_barriers] = barrier weights
  //paramrs[4+num_barriers+num_stages] = stage_weights
  //params[4+num_barriers+num_barriers+num_stages*num_barriers] = mean_values[]

  //par_iter+num_stages+k*num_barriers+i
  Int_t numBar = (Int_t) par[0];
  Int_t numStages = (Int_t) par[1];
  
  Double_t par_iter = 4+numBar+numStages;

  for(Int_t j = 0; j < numStages; j++) {
    tval = 0;
    for (Int_t i = 0; i < numBar; i++) {    
      //Fix the iterator here...
      par_iter = 4 + numBar + numStages + j*numBar;
      //cout << "par_iter = " << par_iter+i << endl;
      
      tval += par[i+4]*TMath::Gaus(xx, par[par_iter+i], (par[par_iter+i]/10),kFALSE);
    }    
    val += par[4+numBar+j]*tval/par[2]; //i.e., level weight * level value / (total barrier weight)   
    //cout << "Total Barrier Weight: " << par[1] << endl;
  }
  return (val/par[3]); //Normalize by total stage weight
}