#ifdef MPI
 
#include "SA_MIRSolver.h"
#include "AllSchedulers.h"
#include "SA_MIRScheduler.h"
#include "my_math.h"
//#include <strstream.h>
#include <strstream>
#include <unistd.h>


// \******************************************************
// den gewaehlten SA-Algo ausfuehren
// \******************************************************
float	SA_MIRSolver::RunAnnealing()
{
	if ( Algorithm == NULL )
	{
		if (myrank == 0 || myrank == -1)
			std::cout << "No algorithm specified !\n";
		return -1;
	}		
	else if ( strcmp(Algorithm, "MIR") == 0 )
	{
		return MIR();
	}
	else if ( strcmp(Algorithm, "SequentialMIR") == 0 )
	{
		return SequentialMIR();
	}
	else if ( strcmp(Algorithm, "TestMIRScheduler") == 0 )
	{
//		TestMIRScheduler();
		return 1;
	}
	else if ( strcmp(Algorithm, "MIRTest") == 0 )
	{
//		return MIRTest();
	}
	else
	{
		if (myrank == 0 || myrank == -1)
			std::cout << "Unknown algorithm "<<Algorithm<<" specified!\n";
		return -1;
	}
}


// \******************************************************
//
// \******************************************************
float SA_MIRSolver::SequentialMIR()
{
//InitStates();	// Initialloesung erzeugen
//	Output();
//	WriteSolution();
//return(1);

	time_t starttime, endtime;
	starttime = time(NULL);
	Cooler->StartClock();
	float SAcost = 0;
	
	if ( myrank == 0 )
	{
		ShowConfig();
	}

	PhaseTwo(Minimum_RunLength, Maximum_RunLength, Betta_Runtime, Samples);
	
	N1_MIR=Maximum_RunLength;
 
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Bcast(&N1_MIR,1,MPI_INT,0,MPI_COMM_WORLD);
	
	if ( myrank == 0 )
	{
		ShowConfig();
		std::cout << "We'll make " << number_of_processors << " runs of length " << N1_MIR << std::endl;
	}
	
	SARun(N1_MIR);
	SATime = time(NULL) - starttime;
	BcastBestSolution(MPI_COMM_WORLD);
	
	PostprocTime = time(NULL);	
	P->Postprocessing(*BestSolution);
	PostprocTime = time(NULL) - PostprocTime;

	SAcost = P->GetCost(*State);
	std::cout << "Cost after " << N1_MIR << " steps on " << myrank << " " << SAcost << std::endl;
	Cooler->SetNewE(SAcost);
	// Bewirkt, dass in Output Mittelwert aller Endkosten aller durchgefuehrten Runs
	// geschrieben werden

	std::ofstream output;
	output.open(OutputFileName, std::ios::app);
	OutputStatistics(output,MPI_COMM_WORLD);
	WriteSolutionParallel(MPI_COMM_WORLD);

	return P->GetCost(*BestSolution);
}



float SA_MIRSolver::MIR()
{
	time_t starttime, endtime;
	starttime = time(NULL);

	if ( myrank == 0 )
	{
		ShowConfig();
	}

	ParallelRuns(Minimum_RunLength, Maximum_RunLength, Betta_Runtime, Samples);

	endtime = time(NULL);
	SATime = endtime - starttime;
	
	BcastBestSolution(MPI_COMM_WORLD);

	PostprocTime = time(NULL);	
	P->Postprocessing(*BestSolution);
	PostprocTime = time(NULL) - PostprocTime;

	std::ofstream output;
	output.open(OutputFileName, std::ios::app);
	OutputStatistics(output,MPI_COMM_WORLD);
	WriteSolution();

	return P->GetCost(*BestSolution);
}


// \******************************************************
//
// \******************************************************
float SA_MIRSolver::ParallelRun(int runlength)
{
time_t starttime, endtime;
starttime = time(NULL);
	(*P.*GetInitialSolution_fp)(*State);
	Cooler->SetStartE(P->GetCost(*State));
	//	std::cout << "Begin " << runlength << " " << P->GetCost(*State) << std::endl;
	((SA_MIRScheduler*)Cooler)->SetDesiredRunLength(runlength);	// sollte vermieden werden
	do
	{
		do
		{
			P->GetNeighbor(*State,Cooler->GetRelativeT());
			if ( Cooler->Accept(P->GetCost(*State)) )
			{
				P->UpdateSolution(*State);
				if ( Cooler->BestFound() )
				{
					P->Copy(*State,*BestSolution);
//					((MIRScheduler*)Cooler)->SetTbest(Cooler->GetT());
				}	
			}
			else
			{
				P->ResetSolution(*State);
			}
		} while (!Cooler->Equilibrium());
		Cooler->UpdateTemperature();
	} while (!Cooler->Frozen());
endtime = time(NULL);
	float Result_local = P->GetCost(*State);
	//	std::cout << std::endl << "Proc " << myrank << " after RunLength " << runlength
	//	 << " Cost " << Result_local << " in " << endtime-starttime << " seconds" << std::endl;
	float Result_global;
        int flag;
        MPI_Initialized(&flag);
        if (flag)
        {
		MPI_Barrier(MPI_COMM_WORLD);
		if ( OptType == Opt::MIN )
		{
			MPI_Reduce(&Result_local,&Result_global,1,MPI_FLOAT,MPI_MIN,0,MPI_COMM_WORLD);
		}
		else
		{
			MPI_Reduce(&Result_local,&Result_global,1,MPI_FLOAT,MPI_MAX,0,MPI_COMM_WORLD);		
		}
	}
	else
		Result_global = Result_local;
if ( myrank == 0 )
{
  int nproc;
  MPI_Comm_size(MPI_COMM_WORLD,&nproc);
  std::cout << nproc << " Processors " << "RunLength " << runlength
	 << " Objective " << Result_global << " in " << endtime-starttime << " seconds" << std::endl;
}
	return Result_global;
}

// \******************************************************
// input:
//	minimum_runlength, maximum_runlength, betta_runtime, samples
// purpose:
//	execute parallel(!) SA-runs of length from minimum_runlength to maximum_runlength by betta_runtime, each samples-time
//? output:
//?	Run_number - number of executed runs
//?	RunLength[Run_number]	- run legths
//?	EndCost[Run_number]		- SA-end costs
// \******************************************************
float SA_MIRSolver::ParallelRuns(int minimum_runlength, int maximum_runlength, float betta_runtime, int samples)
{
	int		actual_run_length;
//	int		Samples_Schranke = (int)ceil(float(samples)/number_of_processors);
	int		Samples_Schranke = samples;

	InitStates();													// create initial solution
	Cooler->StartClock();
	Cooler->WarmingUp(*P,*State,*BestSolution,MPI_COMM_WORLD);

	std::cout << "Barrier 1 " << minimum_runlength << " " << maximum_runlength << " " << samples << std::endl;	
  	MPI_Barrier(MPI_COMM_WORLD);
	
	for(int sample_number=0; sample_number < Samples_Schranke; sample_number++)
	{
		for( actual_run_length = minimum_runlength, Run_number = 0;
			 actual_run_length < maximum_runlength;
			 actual_run_length = (int)ceil(actual_run_length*betta_runtime) )
		{
			float SARunCost = ParallelRun(actual_run_length);			// execute a parallel SA-run of this length
			EndCost[Run_number] += SARunCost/Samples_Schranke;

			if ( myrank == 0 )
			{
				RunLength[Run_number] = actual_run_length;
			}
			Run_number++ ;
		}
	}

	std::cout << "Barrier 2" << std::endl;	
	// in AverageCost die Durchschnitte ueber alle Prozessoren sameln
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Reduce(EndCost,AverageCost,Run_number,MPI_FLOAT,MPI_SUM,0,MPI_COMM_WORLD);

	if (myrank == 0)
	  {
		for (int i=0; i<Run_number; i++)
		{
		  // AverageCost[i]	/= number_of_processors;
		
			int nproc;
			MPI_Comm_size(MPI_COMM_WORLD,&nproc);

			std::cout << "Average: " << nproc << " Processors " << "RunLength " << RunLength[i]
			     << " Objective " << AverageCost[i]	<< std::endl;
		}
	  }

	// in E_best wird auf Prozessor 0 das Minimum ueber alle bis jetzt gesehene Kosten gespeichert
	Collect_Ebest(MPI_COMM_WORLD);

	return P->GetCost(*BestSolution);
}


 
// \******************************************************
//
// \******************************************************
float SA_MIRSolver::oldMIR()
{
	time_t starttime, endtime;
	starttime = time(NULL);
	if ( myrank == 0 )
	{
		ShowConfig();
	}
	PhaseTwo(Minimum_RunLength, Maximum_RunLength, Betta_Runtime, Samples);
	MPI_Bcast(&MIR_Daten,2,MPI_INT,0,MPI_COMM_WORLD);
	if ( MIR_Daten.RunLaenge != 0 && MIR_Daten.AnzahlRuns != 0 )
	{
		ShowConfig();
		float SARunCost = 0;
		for (int i=0; i < MIR_Daten.AnzahlRuns; i++)
		{
			SARunCost += SARun(MIR_Daten.RunLaenge);
		}
		SARunCost /= MIR_Daten.AnzahlRuns;
 
		Cooler->SetNewE(SARunCost);
		// Bewirkt, dass in Output Mittelwert aller Endkosten aller durchgefuehrten Runs
		// geschrieben werden
	}
	
	endtime = time(NULL);
	SATime = endtime - starttime;
	
	BcastBestSolution(MPI_COMM_WORLD);
 
	PostprocTime = time(NULL);	
	P->Postprocessing(*BestSolution);
	PostprocTime = time(NULL) - PostprocTime;
	std::ofstream output;
	output.open(OutputFileName, std::ios::app);
	OutputStatistics(output,MPI_COMM_WORLD);
	WriteSolutionParallel(MPI_COMM_WORLD);
	return P->GetCost(*BestSolution);
}
// \******************************************************
// input:
//	minimum_runlength, maximum_runlength, betta_runtime, samples
// purpose:
//	execute SA-run of length from minimum_runlength to maximum_runlength by betta_runtime, each samples-time
// output:
//	Run_number - number of executed runs
//	RunLength[Run_number]	- run legths
//	EndCost[Run_number]		- SA-end costs
// \******************************************************
float SA_MIRSolver::PhaseTwo(int minimum_runlength, int maximum_runlength, float betta_runtime, int samples)
{
	int		actual_run_length;
	int		Samples_Schranke = (int)ceil(float(samples)/number_of_processors);
	if ( myrank == 0 && VerboseMode )
	{
		std::cout << "Process " << myrank << ": Phase Two begins..." << std::endl;
	}
 
	InitStates();													// create initial solution
	Cooler->WarmingUp(*P,*State,*BestSolution,MPI_COMM_WORLD);
	
  	MPI_Barrier(MPI_COMM_WORLD);
	for(int sample_number=0; sample_number < Samples_Schranke; sample_number++)
	{
		for( actual_run_length = minimum_runlength, Run_number = 0;
			 actual_run_length < maximum_runlength;
			 actual_run_length = (int)ceil(actual_run_length*betta_runtime) )
		{
			float SARunCost = SARun(actual_run_length);				// execute a SA-run of this length
			EndCost[Run_number] += SARunCost/Samples_Schranke;
			if ( myrank == 0 )
			{
				RunLength[Run_number] = actual_run_length;
			}
			Run_number++ ;
		}
	}
 
	// in AverageCost die Durchschnitte ueber alle Prozessoren sameln
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Reduce(EndCost,AverageCost,Run_number,MPI_FLOAT,MPI_SUM,0,MPI_COMM_WORLD);
 
	// in E_best wird auf Prozessor 0 das Minimum ueber alle bis jetzt gesehene Kosten gespeichert
	Collect_Ebest(MPI_COMM_WORLD);
 
	if ( myrank == 0 && VerboseMode )
	{
		std::cout << "Process " << myrank << ": Phase Two completed" << std::endl;
	}
 
	Setup();
 
	return P->GetCost(*BestSolution);
}


// \******************************************************
//
// \******************************************************
void SA_MIRSolver::Setup(void)
{
	ComputeDeltaCost();
	ComputePowerLawFit();
 
	N1_MIR = ComputeN1_MIR(Alpha_praxis, K_praxis);
	
	MIR_Daten.AnzahlRuns = 1;
	MIR_Daten.RunLaenge = (int)ceil(N1_MIR/Fraction);
}


// \******************************************************
//
// \******************************************************
float SA_MIRSolver::SARun(int runlength)
{
time_t starttime, endtime;
starttime = time(NULL);

	(*P.*GetInitialSolution_fp)(*State);
	((SA_MIRScheduler*)Cooler)->SetDesiredRunLength(runlength);	// sollte vermieden werden

	do
	{
		do
		{
			P->GetNeighbor(*State,Cooler->GetRelativeT());
			if ( Cooler->Accept(P->GetCost(*State)) )
			{
				P->UpdateSolution(*State);
				if ( Cooler->BestFound() )
				{
					P->Copy(*State,*BestSolution);
//					((MIRScheduler*)Cooler)->SetTbest(Cooler->GetT());
				}	
			}
			else
			{
				P->ResetSolution(*State);
			}
		} while (!Cooler->Equilibrium());
//		if ( myrank == 0 && VerboseMode )
//		{
//			((MIRScheduler*)Cooler)->OutputSubchainStatistics(std::cout); std::cout <<std::endl;
//		}
		Cooler->UpdateTemperature();

	} while (!Cooler->Frozen());
	if (Cooler->TimeExceeded()) {
	  std::cout << "TIME LIMIT EXPIRED!!! \n Annealing process stoppt! \n";
	}
	P->Reoptimize(*State);
	P->Postprocessing(*State);

endtime = time(NULL);
std::cout << std::endl << "Nr " << myrank << " hat bei RunLength " << runlength
	 << " Cost " << P->GetCost(*State) << " in " << endtime-starttime << " seconds" << std::endl;

	return P->GetCost(*State);
}



// \******************************************************
// Input:	array EndCost[0..Run_number-1] with end-costs of runs
// Output:	array DeltaCost[0..Run_number-1]  with DeltaCost[i]=AverageCost[i]-Ebest
//			array LogDeltaCost[0..Run_number-1]
// \******************************************************
void SA_MIRSolver::ComputeDeltaCost(void)
{
	if ( myrank == 0 )
	{
PowerLaw_average.open("pics/PowerLawAverage");
		int i;

		// compute DeltaCost per definition
		for (i=0; i<Run_number; i++)
		{
			AverageCost[i]	/= number_of_processors;
if ( PowerLaw_average )
PowerLaw_average << RunLength[i] << " " << AverageCost[i] << std::endl;
			
			if ( Ebest_user == 0 )
			{
				DeltaCost[i] = FABS(AverageCost[i] - Ebest);
			}
			else
			{
				DeltaCost[i] = FABS(AverageCost[i] - Ebest_user);
			}
		}

		// we'll fit in log-scale
		for (i=0; i<Run_number; i++)
		{
			if ( DeltaCost[i] != 0 )
			{
				LogDeltaCost[i] = log(DeltaCost[i]);
			}
			else
			{
				LogDeltaCost[i] = 0;
			}
		}
	}
}


// \******************************************************
// input:	LogDeltaCost[0..Run_number-1] and RunLength[0..Run_number-1]
// output:	
//
// Our Situation is:							     DeltaCost = (K_praxis/RunLength)^Alpha_praxis
// least squares method estimation of the line	log(DeltaCost) = Alpha_praxis*LogK_praxis + MinusAplha_praxis*log(RunLength)
// RunLength is our Variable
// we search for LogK_praxis and MinusAlpha_praxis
// MinusAlpha = Sum(i=[0..Run_number])[(LogDeltaCost[i]-AverageDeltaCost)(log(RunLength[i])-AverageLogRunLength)] /
//				Sum(i=[0..Run_number])[log(RunLength[i])-AverageLogRunLength]^2
// LogK_praxis = AverageDeltaCost - MinusAlpha*AverageLogRunLength
// K_praxis   = exp(LogK_praxis)
// \******************************************************
void SA_MIRSolver::ComputePowerLawFit(void)
{
	if ( myrank == 0 )
	{
		float LogK_praxis=0;
		
		// Durchschnittsberechnung
		float AverageDeltaCost=0, AverageLogRunLength=0;
		for (int j=0; j<Run_number; j++)
		{
			AverageDeltaCost += LogDeltaCost[j];
			AverageLogRunLength += log(RunLength[j]);
		}
		AverageDeltaCost /= Run_number;
		AverageLogRunLength /= Run_number;
		
		float Summe1=0,Summe2=0;
		for (int k=0; k<Run_number; k++)
		{
			Summe1 += (LogDeltaCost[k]-AverageDeltaCost)*(log(RunLength[k])-AverageLogRunLength);
			Summe2 += pow(log(RunLength[k])-AverageLogRunLength,2);
		}
		
		if ( Summe2 != 0 )
		{
			Alpha_praxis = - Summe1 / Summe2;
			std::cout << "Alpha_praxis = " << Alpha_praxis << std::endl;
		}
		else
		{
			std::cerr << "SA_MIRSolver: K_praxis can't be calculated" << std::endl;
			Alpha_praxis = 0;
		}

		if ( Alpha_praxis != 0 )
		{
			LogK_praxis = (AverageDeltaCost + Alpha_praxis*AverageLogRunLength)/Alpha_praxis;
			K_praxis = exp(LogK_praxis);
			std::cout << "LogK_praxis = " << LogK_praxis << std::endl;
			std::cout << "K_praxis = " << K_praxis << std::endl;
		}
		else
		{
			std::cerr << "SA_MIRSolver: K_praxis can't be calculated" << std::endl;
		}

// Paare (log(DeltaCost),log(run_length)) bestimmen und ausgeben
for (int l=0; l<Run_number; l++)
{
	PowerLaw_outputfile << log(RunLength[l]) << " " << LogDeltaCost[l] << std::endl;
	PowerLaw_fit << log(RunLength[l]) << " " << Alpha_praxis*LogK_praxis - Alpha_praxis*log(RunLength[l]) << std::endl;
}
//Plot_Test(K,Alpha);
	}	
}


// \******************************************************
//
// \******************************************************
int	SA_MIRSolver::ComputeN1_MIR(float Alpha_pr, double K_pr)
{	
	int N1_mir=0;
	float ebest;
	
	if ( myrank == 0 )
	{
		float Solution_Quality;
		if ( Ebest_user == 0 )
		{
			ebest = Ebest;
		}
		else
		{
			ebest = Ebest_user;
		}
		
		Solution_Quality = Epsilon * ebest;
		
		if ( Alpha_pr != 0 && pow(Solution_Quality,1/Alpha_pr) != 0 )
		{
			N1_mir = (int)floor(K_pr * pow(Solution_Quality,-1/Alpha_pr));

			std::cout << "N1_MIR = " << N1_mir << std::endl;
			if ( N1_mir !=0 )
			{
std::cout << "Loesungsqualitaet nach " << N1_mir << " Iterationen " << pow(K_pr/N1_mir,Alpha_pr)/ebest << std::endl;
				EndSolutionQuality = pow(K_pr/N1_mir,Alpha_pr)/ebest;
			}
		}
		else
		{
			std::cerr << "SA_MIRSolver: N1_MIR can't be calculated" << std::endl;
			N1_mir = (int)P->GetLocalN();
		}
	}
	return N1_mir;
}


// \************************************************************************************************************



// \******************************************************
//
// \******************************************************
void SA_MIRSolver::Collect_Ebest(MPI_Comm comm)
{
	float Ebest_local = Cooler->GetBestE();

	if ( myrank >=0 )	// d.h. falls MPI initialisiert ist
	{
		MPI_Barrier(MPI_COMM_WORLD);

		if ( OptType == Opt::MIN )
		{
			MPI_Reduce(&Ebest_local,&Ebest,1,MPI_FLOAT,MPI_MIN,0,comm);
		}
		else
		{
			MPI_Reduce(&Ebest_local,&Ebest,1,MPI_FLOAT,MPI_MAX,0,comm);		
		}
	}
	else
		Ebest = Ebest_local;
	if (myrank ==0 && VerboseMode) {std::cout << "Ebest  " << Ebest << std::endl;}
}



// \******************************************************
//
// \******************************************************
SA_MIRSolver::SA_MIRSolver(int *argc, char ***argv, SA_Problem & prb):
	SA_ParSolver(prb),
	Alpha_MIR(0),
	MinusAlpha_MIR(0),
	K_MIR(0),
	B_MIR(0),
	LogB_MIR(0),
	N1_MIR(0),
	Epsilon(0.1),
	s_MIR(0),
	Ns_MIR(0),
	q_MIR(0),
	Nq_MIR(0),
	Minimum_RunLength(0),
	Maximum_RunLength(0),
	Run_number(0),
	Betta_Runtime(1.1),
	RunFactor(5),
	Samples(0),
	Steptime(0),
	Runtime(0),
	Estar(0),
	Ebest_user(0),
	Fraction(1),
	samirsolvoutput("SA_MIRSolver","SA_Output.rsc")
{
 
 
	int flag=0;
	MPI_Initialized(&flag);
	if ( !flag )
	{
		std::cout << "MPI nicht initialisiert! Initialisiere MPI..." << std::endl;
		std::cout << "argc " << *argc << std::endl;
		std::cout << "argv " << *argv[0] << std::endl;
		MPI_Init(argc, argv);
		std::cout << "MPI erfolgreich initialisiert!" << std::endl;
	}
	else
	{
		std::cout << "MPI bereits initialisiert!" << std::endl;		
	}
 
	MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
	MPI_Comm_size(MPI_COMM_WORLD,&number_of_processors);
	
	// Falls Multiprozessorbetrieb : jeder Proc. benutzt eigenen Zufallszahlstream
	if( myrank >=0 )
	{
		Select_Stream(myrank);
	}


	int dummy;
	MPI_Get_processor_name(processor_name,&dummy);
	std::cout << "Process " << myrank << " is alive on " << processor_name << std::endl; //" with pid " << getpid() << endl;
	
  	MPI_Barrier(MPI_COMM_WORLD);
  	
	MIR_Daten.AnzahlRuns = 1;
	MIR_Daten.RunLaenge = (int)ceil(P->GetLocalN());

	EndCost = new float[MAX_Run_number];
	for (int j=0; j<MAX_Run_number; j++)
	{
		EndCost[j] = 0;
	}

	if( myrank == 0 )
	{
		AverageCost	= new float[MAX_Run_number];
		DeltaCost	= new float[MAX_Run_number];
		LogDeltaCost	= new float[MAX_Run_number];
		RunLength	= new int[MAX_Run_number];

		for (int j=0; j<MAX_Run_number; j++)
		{
			AverageCost[j] = DeltaCost[j] = LogDeltaCost[j] = RunLength[j] = 0;
		}
	}
	
	PowerLaw_outputfile.open("pics/PowerLaw");
	PowerLaw_fit.open("pics/PowerLawFit");
	PowerLaw_fit1.open("pics/PowerLawFit1");
 
	samirsolvoutput.EnableLevelNumbers(false);
	samirsolvoutput.EnableIdentifier(true);
	samirsolvoutput.SetOutputLevel(10);
 
	samirsolvoutput.AddVariable(1,SA_INT, &number_of_processors);
	samirsolvoutput.AddVariable(2,SA_FLOAT, &Alpha_MIR);
	samirsolvoutput.AddVariable(3,SA_FLOAT, &MinusAlpha_MIR);
	samirsolvoutput.AddVariable(4,SA_FLOAT, &Alpha_MIR_old);
	samirsolvoutput.AddVariable(5,SA_FLOAT, &MinusAlpha_MIR_old);
	samirsolvoutput.AddVariable(6,SA_FLOAT, &Alpha_praxis);
	samirsolvoutput.AddVariable(7,SA_FLOAT, &K_MIR);
	samirsolvoutput.AddVariable(8,SA_FLOAT, &K_MIR_old);
	samirsolvoutput.AddVariable(9,SA_DOUBLE, &K_praxis);
	samirsolvoutput.AddVariable(10,SA_FLOAT, &Ebest);    
	samirsolvoutput.AddVariable(11,SA_FLOAT, &Ebest_user);    
	samirsolvoutput.AddVariable(12,SA_FLOAT, &Estar);    
	samirsolvoutput.AddVariable(13,SA_FLOAT, EndCost);    
	samirsolvoutput.AddVariable(14,SA_FLOAT, AverageCost);    
	samirsolvoutput.AddVariable(15,SA_FLOAT, DeltaCost);    
	samirsolvoutput.AddVariable(16,SA_FLOAT, LogDeltaCost);    
	samirsolvoutput.AddVariable(17,SA_FLOAT, &B_MIR);    
	samirsolvoutput.AddVariable(18,SA_FLOAT, &LogB_MIR);    
	samirsolvoutput.AddVariable(19,SA_FLOAT, &LogK_MIR);    
	samirsolvoutput.AddVariable(20,SA_FLOAT, &B_MIR_old);    
	samirsolvoutput.AddVariable(21,SA_FLOAT, &LogB_MIR_old);    
	
	samirsolvoutput.SetPreviousOutput(SA_ParSolver::GetSolverOutput());
}


// \******************************************************
//
// \******************************************************
SA_MIRSolver::~SA_MIRSolver()
{
	delete [] EndCost;

	if( myrank == 0 )
	{
		delete [] AverageCost;
		delete [] DeltaCost;
		delete [] LogDeltaCost;
		delete [] RunLength;
	}
	
	PowerLaw_outputfile.close();
	PowerLaw_fit.close();
	PowerLaw_fit1.close();

	MPI_Finalize();
}

SA_Output *SA_MIRSolver::GetSolverOutput() {
  return(&samirsolvoutput);
}
// \******************************************************
// Methode ReadConfig - Liest die Konfigurations-
// daten ( geklammert mit '{' '}'  ) aus dem istream
// \******************************************************
int SA_MIRSolver::ReadConfig(std::istream & config)
{
	char text[BUFLEN];
	config >> text;
	if ( strcmp(text,"{") != 0)
	{
		std::cout << "Config file section for SA_MIRSolver does not begin with {!\n";
		return FALSE;
	}
	
	while (config >> text && ( strcmp(text,"}") != 0 ))
	{
		if (strcmp(text, "SA_Solver") == 0)
		{
			if ( ! SA_Solver::ReadConfig(config) )
			{
				std::cout <<"ReadConfing for Solver failed!\n";
				return 0;
			}
		}
		else if (strcmp(text, "SA_Scheduler") == 0)
		{
			// Instantiierung des Schedulers
			config >> text;
			if (strcmp(text, "SA_MIRScheduler") == 0)
			{
				SetScheduler(new SA_MIRScheduler());
			}
			else if (strcmp(text, "SA_EasyScheduler") == 0)
			{
				config >> text;
				SetScheduler(new SA_EasyScheduler());
			}
			//			else if (strcmp(text, "FloodScheduler") == 0)
			//	{
			//		Cooler = new FloodScheduler();
			//	}
			// else if (strcmp(text, "OttenScheduler") == 0)
			//{
			//	config >> text;
			//	Cooler = new OttenScheduler();
			//}
			else
			{
				std::cout << "Wrong Scheduler " << text << "!\n";
				return 0; 
			}
			
			// Einlesen der Konfigurationsdaten des Schedulers
			if ( ! Cooler->ReadConfig(config) )
			{
				std::cout <<"ReadConfing for Scheduler failed!\n";
				return 0;
			}
		}
		else if (strcmp(text, "RunFactor") == 0)
		{
			config >> text;
			RunFactor = atof(text);
		}
		else if (strcmp(text, "Betta_Runtime") == 0)
		{
			config >> text;
			Betta_Runtime = atof(text);
		}
		else if (strcmp(text, "Maximum_RunLength") == 0)
		{
			config >> text;
			Maximum_RunLength = atoi(text);
		}
		else if (strcmp(text, "Minimum_RunLength") == 0)
		{
			config >> text;
			Minimum_RunLength = atoi(text);
		}
		else if (strcmp(text, "Samples") == 0)
		{
			config >> text;
			Samples = atoi(text);
		}
		else if (strcmp(text, "Epsilon") == 0)
		{
			config >> text;
			Epsilon = atof(text);
		}
		else if (strcmp(text, "Derivation") == 0)
		{
			config >> text;
			Derivation = atof(text);
		}
		else if (strcmp(text, "Fraction") == 0)
		{
			config >> text;
			Fraction = atoi(text);
		}
		else if (strcmp(text, "Ebest") == 0)
		{
			config >> text;
			Ebest_user = atof(text);
		}
		else if (strcmp(text, "algorithm") == 0)
		{
			config >> text;
			delete Algorithm;
			Algorithm = new char[strlen(text)+1];
			(strcpy(Algorithm, text) == 0);
		}
		else 
		{
			std::cout <<"In SA_MIRSolver::ReadConfig unrecognized keyword \""<<text<<"\""<<std::endl;
		}
	}
	
	if ( OptType == Opt::MIN )
	{
		Ebest = INFINITY;
	}
	else
	{
		Ebest = 0;
	}
	
	if ( Betta_Runtime == 0 )
	{
		Betta_Runtime = 1.1;
	};
	if ( RunFactor == 0 )
	{
		RunFactor = 5;
	};
	if ( Minimum_RunLength == 0 )
	{
		Minimum_RunLength = (int)ceil(P->GetLocalN());
	};
	if ( Maximum_RunLength == 0 )
	{
		Maximum_RunLength = (int)ceil(RunFactor*P->GetLocalN());
	};
	if ( Samples == 0 )
	{
		Samples = 10;
	};

	
	return (strcmp(text,"}") == 0 );
}


// \******************************************************
// Methode ShowConfig - Gibt die aktuelle Konfiguration
// aus (nur zur Kontrolle).
// \******************************************************
void SA_MIRSolver::ShowConfig()
{
	SA_Solver::ShowConfig();
	std::cout << "SA_MIRSolver Konfiguration:" << std::endl
		<< "\tEpsilon:\t\t" << Epsilon << std::endl
		<< "\tDerivation:\t\t" << Derivation << std::endl
		<< "\tMinimum_RunLength:\t" << Minimum_RunLength << std::endl
		<< "\tMaximum_RunLength:\t" << Maximum_RunLength << std::endl
		<< "\tBetta_Runtime:\t\t" << Betta_Runtime << std::endl
		<< "\tSamples:\t\t" << Samples << std::endl
		<< "\tAnzahlRuns:\t\t" << MIR_Daten.AnzahlRuns << std::endl
		<< "\tRunLaenge:\t\t" << MIR_Daten.RunLaenge << std::endl;
}


// \******************************************************
// Methode OutputStatistics - Schreibt die Ergebnisse in die Logdatei.
// \******************************************************
void SA_MIRSolver::CreateLog(std::ostream & output, MPI_Comm comm)
{
	float in,out;
	in = Cooler->GetE();
	MPI_Reduce(&in,&out,1,MPI_FLOAT,MPI_SUM,0,comm);

	if ( myrank == 0 )
	{
		float EndE = out/number_of_processors;
		double SolQual;
		if ( Ebest_user == 0 )
		{
			SolQual = FABS(EndE-Ebest)/Ebest;
		}
		else
		{
			SolQual = FABS(EndE-Ebest_user)/Ebest_user;			
		}
		  
	// Ausgabe der Daten: Solv= Solver, NPROC=Anzahl der Prozessoren, F=Dateiname, Z=Lese-/Rechenzeit 
		output
			<< "Solv=MIR NPROC=" << number_of_processors
			<< " F=" << filename(DataFileName)
			<< " Z=" << FileReadTime << "/" << InitialSolutionTime << "/" << SATime << "/" << PostprocTime
			<< " StrtSl=" << ( GetInitialSolution_fp == & SA_Problem::GetRandomSolution ? "Rnd":"Init" )

//std::cout << "TargetSolQual=" << EndSolutionQuality << " " << "RealSolutionQuality=" << FABS(Cooler->GetE()-Ebest_user)/Ebest_user << std::endl;

			<< " Samples=" << Samples
			<< " Epsilon=" << Epsilon
			<< " Derivation=" << Derivation
			<< " TargetSolQual=" << EndSolutionQuality
			<< " SolQual=" << SolQual
			<< " MinRL=" << Minimum_RunLength
			<< " MaxRL=" << Maximum_RunLength
			<< " RunB=" << Betta_Runtime
			<< " Alpha_MIR=" << Alpha_praxis
			<< " K_MIR=" << K_praxis
			<< " N1_MIR=" << N1_MIR
			<< " Frac=" << Fraction
//			<< " Ns_MIR=" << Ns_MIR
//			<< " s_MIR=" << s_MIR
//			<< " Nq_MIR=" << Nq_MIR
//			<< " q_MIR=" << q_MIR
			<< " ";
			
	    Cooler->CreateLog(output,comm);
	    output << std::endl;
	}
	else
	{
		Cooler->CreateLog(std::cout,comm);
	}
}

//// \******************************************************
////
//// \******************************************************
//float SA_MIRSolver::MIRTest()
//{
//	time_t starttime, endtime;
//	starttime = time(NULL);
//	
//std::cout << "Process " << myrank << ": Phase Two begins..." << std::endl;
//	PhaseTwo();
//	Setup();
//std::cout << "Process " << myrank << ": Phase Two completed" << std::endl;
//
//	if ( myrank == 0 )
//	{
//		if ( N1_MIR == 0 ) //|| N1_MIR == NaN )
//		{
//			N1_MIR = (int)ceil(P->GetLocalN());
//		}
//	}
//	
//	MPI_Barrier(MPI_COMM_WORLD);
//	MPI_Bcast(&N1_MIR,1,MPI_INT,0,MPI_COMM_WORLD);
//	
//	if ( myrank == 0 )
//	{
//		Cooler->ShowConfig(std::cout);
//		ShowConfig();
//		std::cout << "We'll make " << number_of_processors << " runs of length " << N1_MIR << std::endl;
//	}
//
//	int Real_N1_MIR = N1_MIR;
//// For the use with CC added:
//	const int MaxProcN = 128;
//	if (number_of_processors > MaxProcN)
//	{
//		std::cout << "Too many procs, please specify greater MaxProcN or do dynamic memory ! \n";
//		exit(1);
//	}
//	float	SA_Costs[100][MaxProcN];
//	for (int ii=0; ii<100; ii++)
//	{
//		for( int jj=0; jj<number_of_processors; jj++)
//		{
//			SA_Costs[ii][jj] = 0;
//		}
//	}
//	
//	int real_size=0;
//	
//	for (int i=2;Real_N1_MIR>=N1_MIR/Fraction;Real_N1_MIR=N1_MIR/i,i++)
//	{
//		SARun(Real_N1_MIR);
//		P->Postprocessing(*State);
//		SA_Costs[i-2][myrank] = P->GetCost(*State);
//		std::cout << "Cost after " << Real_N1_MIR << " steps on " << myrank << " " << SA_Costs[i-2][myrank] << std::endl;
//		real_size++;
//	}
//	MPI_Barrier(MPI_COMM_WORLD);
//
//	float	SA_Costs_global[100][MaxProcN];
//	for (int i1=0; i1<100; i1++)
//	{
//		for( int j1=0; j1<number_of_processors; j1++)
//		{
//			SA_Costs_global[i1][j1] = 0;
//		}
//	}
//
//
//	MPI_Reduce(SA_Costs,SA_Costs_global,100*number_of_processors,MPI_FLOAT,MPI_SUM,0,MPI_COMM_WORLD);
//
//	if (myrank ==0)
//	{
//		for (int l=0; l<real_size; l++)
//		{
////			std::cout << "SA_Costs_global " << " by RunLength " << N1_MIR/(l+1) << std::endl;
////			for (int k=0; k < number_of_processors; k++)
////			{
////				std::cout << SA_Costs_global[l][k] << " ";
////			}
//			std::cout << std::endl;
//			std::cout << "Best by RunLength " << N1_MIR/(l+1) << " " << BestOf(SA_Costs_global[l],number_of_processors) << std::endl;
//			std::cout << "Average by RunLength " << N1_MIR/(l+1) << " " << AverageOf(SA_Costs_global[l],number_of_processors) << std::endl;
//		}
//	}
//	return 0;
//}
//
//// \******************************************************
////
//// \******************************************************
//void SA_MIRSolver::Setup(void)
//{
////	ComputeDeltaCost();
////	float Alpha_mir,B_mir,K_mir;
////	ComputePowerLawFit(Alpha_mir,B_mir,K_mir);
////
////	
//////	K_MIR = ComputeK_MIR(Alpha_MIR);
////
////	// Der Ansatz von Lee(1995)
////	int N1_mir = ComputeN1_MIR(Alpha_mir,B_mir,K_mir);
////
////	int anzahlruns_Lee,runlaenge_Lee;
////	float runtime_Lee;
////	Ns_MIR = ComputeNs_MIR(Alpha_mir,K_mir,N1_mir,s_MIR,anzahlruns_Lee,runlaenge_Lee,runtime_Lee);
////	if( myrank == 0 )
////	{
////		std::cout << "Lee: We'll make " << s_MIR << " runs of length " << runlaenge_Lee
////			 << " ~ " << K_mir << " ; each processor performs " << anzahlruns_Lee
////			 << " runs" << std::endl;
////	}
////
////	// neuer Ansatz
////	int anzahlruns_PB,runlaenge_PB,q_mir;
////	float runtime_PB;
////	Nq_MIR = ComputeNq_MIR(Alpha_mir,K_mir,N1_mir,q_mir,anzahlruns_PB,runlaenge_PB,runtime_PB);
////	if( myrank == 0 )
////	{
////		std::cout << "Paderborn: We'll make " << q_mir << " runs of length " << runlaenge_PB
////			 << " ~ " << exp(1)*K_mir << " ; each processor performs " << anzahlruns_PB
////			 << " runs" << std::endl;
////	}
////	q_MIR = q_mir;
////	
////	// Derivation-Ansatz
////	N1_mir = ComputeN1_MIR_Derivation(Alpha_mir,B_mir);
////	N1_MIR = N1_mir;
////	
////	if ( N1_MIR < P->GetLocalN() )
////	{
////		N1_MIR = (int)ceil(P->GetLocalN());
////	}
////	
//////	MIR_Daten.AnzahlRuns = s_MIR/number_of_processors;
//////	MIR_Daten.RunLaenge = (int)(exp(1)*K_mir);
////	MIR_Daten.AnzahlRuns = number_of_processors;
////	MIR_Daten.RunLaenge = N1_MIR;
////	
////	Alpha_MIR = Alpha_mir;
////	B_MIR = B_mir;
////	K_MIR = K_mir;
////
//////	Plot(B_MIR,Alpha_MIR);	
////
//////	if ( myrank == 0 )
//////	{
//////		Prepare_fitsess(Alpha_mir, B_mir, Ebest);
//////		fitsess(Alpha_mir, B_mir, Run_number);
//////	}
//}
//
//// \******************************************************
////
//// \******************************************************
////void SA_MIRSolver::ComputeK(float n, float kappa_n, double emin, float K, float alpha, float deviation)
////{
////	float myK1=0, myK2=0;
////
////	Collect_Ebest(MPI_COMM_WORLD);
////	
////	// first estimation
////	std::cout << "Ebest  " << Ebest << std::endl;
////	myK1 = K/pow((kappa_n-Ebest)/Ebest,1/alpha);
////	std::cout << "myK1 " << myK1 << std::endl;
////	
////	// second estimation
////	myK2 = n*pow(1-pow(deviation,2)/(2*pow(kappa_n-Ebest,2)),1/alpha);
////	std::cout << "myK2 " << myK2 << std::endl;
////}
//
//
//// \******************************************************
////
//// \******************************************************
////void SA_MIRSolver::ComputeK_first(float K_praxis, float Alpha_praxis)
////{
////	float Ebest_my = 0;
////	if ( Ebest_user == 0 )
////	{
////		Collect_Ebest(MPI_COMM_WORLD);
////		Ebest_my = Ebest;
////	}
////	else
////	{
////		Ebest_my = Ebest_user;
////	}
////
////}
//
//// \******************************************************
//// see page 83, formel (24)
//// Our Situation is:							     DeltaCost = (Estar-Ebest)(K_MIR/RunLength)^Alpha_MIR
////														B_MIR := (Estar-Ebest)K_MIR^Alpha_MIR
////													 DeltaCost = B_MIR * RunLength^(-Alpha_MIR)
//// least squares method estimation of the line	log(DeltaCost) = LogB_MIR + MinusAplha_MIR*log(RunLength)
//// RunLength is our Variable
//// we search for LogB_MIR and MinusAlpha_MIR
//// MinusAlpha_MIR = Sum(i=[0..Run_number])[(LogDeltaCost[i]-AverageDeltaCost)(log(RunLength[i])-AverageLogRunLength)] /
////				Sum(i=[0..Run_number])[log(RunLength[i])-AverageLogRunLength]^2
//// LogB_MIR = AverageDeltaCost - MinusAlpha_MIR*AverageLogRunLength
//// K_MIR = (B_MIR/(Estar-Ebest))^(1/Alpha_MIR)
//// \******************************************************
//void SA_MIRSolver::ComputePowerLawFit(float& Alpha_mir, float& B_mir, float& K_mir)
//{
//	Alpha_mir = B_mir = K_mir = 0;
//	if ( myrank == 0 )
//	{
//		// Durchschnittsberechnung
//		float AverageDeltaCost=0, AverageLogRunLength=0;
//		for (int j=0; j<Run_number; j++)
//		{
//			AverageDeltaCost += LogDeltaCost[j];
//			AverageLogRunLength += log(RunLength[j]);
//		}
//		AverageDeltaCost /= Run_number;
//		AverageLogRunLength /= Run_number;
//		
//		float Summe1=0,Summe2=0;
//		for (int k=0; k<Run_number; k++)
//		{
//			Summe1 += (LogDeltaCost[k]-AverageDeltaCost)*(log(RunLength[k])-AverageLogRunLength);
//			Summe2 += pow(log(RunLength[k])-AverageLogRunLength,2);
//		}
//		MinusAlpha_MIR = Summe1 / Summe2;
//		Alpha_mir = -MinusAlpha_MIR;
//		std::cout << "Alpha_MIR = " << Alpha_mir << std::endl;
//
//		if ( Alpha_mir != 0 )
//		{
//			LogB_MIR = AverageDeltaCost - MinusAlpha_MIR*AverageLogRunLength;
//			B_mir = exp(LogB_MIR);
//			std::cout << "B_MIR = " << B_mir << std::endl;
//		}
//		else
//		{
//			std::cerr << "SA_MIRSolver: B_MIR can't be calculated" << std::endl;
//		}
//		
//		Estar = AverageCost[Run_number-1];
//				
//		if ( Ebest_user == 0 )
//		{
//			if ( Alpha_mir != 0 && B_mir != 0 && Estar != Ebest )
//			{
//				K_mir = pow(B_mir/FABS(Estar - Ebest), 1/Alpha_mir);
//				std::cout << "K_MIR (first)= " << K_mir << std::endl;
//			}
//			else
//			{
//				std::cerr << "SA_MIRSolver: K_MIR (first) can't be calculated" << std::endl;
//			}
//		}
//		else 
//		{
//			if ( Alpha_mir != 0 && B_mir != 0 && Estar != Ebest_user )
//			{
//				K_mir = pow(B_mir/FABS(Estar - Ebest_user), 1/Alpha_mir);
//				std::cout << "K_MIR (first)= " << K_mir << std::endl;		
//			}
//			else
//			{
//				std::cerr << "SA_MIRSolver: K_MIR (first) can't be calculated" << std::endl;
//			}
//		}
//			
//float DiffSum = 0;
//// Paare (log(DeltaCost),log(run_length)) bestimmen und ausgeben
//for (int l=0; l<Run_number; l++)
//{
//DiffSum += pow( (B_mir*pow(RunLength[l],-Alpha_mir)+Ebest) - AverageCost[l], 2);
//PowerLaw_outputfile << log(RunLength[l]) << " " << LogDeltaCost[l] << std::endl;
//PowerLaw_fit << log(RunLength[l]) << " " << LogB_MIR - Alpha_mir*log(RunLength[l]) << std::endl;
//PowerLaw_fit1 << log(RunLength[l]) << " " << log(FABS(Estar-Ebest) * pow(K_mir/RunLength[l],Alpha_mir)) << std::endl;
////PowerLaw_fit1 << log(RunLength[l]) << " " << Alpha_mir*log(K_mir/RunLength[l]) << std::endl;
//}
//std::cout << "With linear fit: DiffSum = " << DiffSum << std::endl;
//	}
//}
//
//
//
//
//
//
//
//
//
//
//
//
//
//
//// \******************************************************
////
//// \******************************************************
////int	SA_MIRSolver::ComputeN1_MIR(float Alpha_mir, float B_mir, float K_mir)
////{	
////	int N1_mir=0;
////	if ( myrank == 0 )
////	{
////		float Solution_Quality;
////		if ( Ebest_user == 0 )
////		{
////			Solution_Quality = Epsilon * Ebest;
////		}
////		else
////		{
////			Solution_Quality = Epsilon * Ebest_user;		
////		}
////		
////		if ( VerboseMode )
////		std::cout << "SolutionQuality=" << Epsilon << std::endl;
////		
////		if ( Alpha_mir != 0 && pow(Solution_Quality,1/Alpha_mir) != 0 )
////		{
////			double buf1 = B_mir/Solution_Quality;
////			double buf2 = 1/Alpha_mir;
////			double buf3 = pow(buf1,buf2);
////			std::cout << buf1 << " hoch " << buf2 << " ist " << buf3 << std::endl;
////			N1_mir = (int)floor(buf3);
////			
//////			N1_mir = (int)floor(pow(B_mir/Solution_Quality,1/Alpha_mir));
////			std::cout << "N1_MIR = " << N1_mir << std::endl;
////			std::cout << "Wahrscheinlichkeit[XN1_MIR != Ebest] = " << pow(K_mir/N1_mir, Alpha_mir) << std::endl;
////
//////			double EndSolutionValue = B_mir*pow(N1_mir,-Alpha_mir)+Ebest_user;
//////			std::cout << "EndSolutionValue1=" << EndSolutionValue << std::endl;
//////			double EndSolutionQuality = FABS(EndSolutionValue-Ebest_user)/Ebest_user;
//////			std::cout << "EndSolutionQuality1=" << EndSolutionQuality << std::endl;		
////
//////			PowerLaw_average.open("pics/PowerLawAverage",ios::app);
//////			PowerLaw_average << N1_mir << " " << Ebest << std::endl;
//////			PowerLaw_average.close();
//////			PowerLaw_outputfile << log(N1_mir) << " " << log(Solution_Quality) << std::endl;
////		}
////		else
////		{
////			std::cerr << "SA_MIRSolver: N1_MIR can't be calculated" << std::endl;
////			N1_mir = (int)P->GetLocalN();
////		}
////	}
////	return N1_mir;
////}
//
//// \******************************************************
//// Eingabe: K_MIR, N1_MIR, Alpha_MIR
//// Ausgabe: Ns_MIR, s_MIR, MIR_Daten.AnzahlRuns, MIR._Daten.RunLaenge, Runtime
//// \******************************************************
//int	SA_MIRSolver::ComputeNs_MIR(float Alpha_mir,
//							 float K_mir,
//							 float N1_mir,
//							 int& s_mir,
//							 int& anzahlruns,
//							 int& runlaenge,
//							 float& runtime)
//{	
//	int Ns_mir = 0;
//	anzahlruns = runlaenge = s_mir = 0;
//	runtime = 0;
//	
//	if ( myrank == 0 )
//	{
//		int Ns_mir_two,Ns_mir_three;
//		
//		if ( exp(1)*K_mir != 0 )
//		{
//			s_mir = (int)ceil(N1_mir/(exp(1)*K_mir));
//		}
//		else
//		{
//			std::cerr << "SA_MIRSolver: s_MIR can't be calculated" << std::endl;
//			s_mir = 1;
//		}
//		
//		if ( s_mir != 0 && exp(1)*pow(K_mir,1/s_mir) != 0 )
//		{
//			Ns_mir = (int)( pow(N1_mir,1+(1/s_mir)) / (exp(1) * pow(K_mir,1/s_mir)) );
//			Ns_mir_two = (int)( pow(N1_mir,1/s_mir) * s_mir * pow(K_mir,1-1/s_mir) );
//			float Optimum_Prob = pow(K_mir/N1_mir, Alpha_mir);
//			Ns_mir_three = (int)( s_mir * K_mir * pow(Optimum_Prob,-1/(s_mir*Alpha_mir)) );
//		}
//		else
//		{
//			std::cerr << "SA_MIRSolver: Ns_MIR can't be calculated" << std::endl;
//			Ns_mir = (int)P->GetLocalN();
//		}
//
//		// check
//		std::cout << "it must be zero : " <<
//			pow( s_mir*K_mir/Ns_mir, Alpha_mir*s_mir ) -
//			pow( K_mir/N1_mir, Alpha_mir ) << std::endl;
//
//		if ( s_mir != 0 && Ns_mir !=0 )
//		{
//			anzahlruns = 1+(s_mir-1)/number_of_processors;
//			runlaenge = 1+(Ns_mir-1)/s_mir;
//		}
//		else
//		{
//			anzahlruns = 1;
//			runlaenge = (int)P->GetLocalN();
//		}
//		runtime = anzahlruns * runlaenge* Steptime;
//
//		std::cout << "Ns_MIR = " << Ns_mir << " Ns_mir_2 " << Ns_mir_two
//			 << " Ns_mir_3 " << Ns_mir_three << " s_MIR = " << s_mir << std::endl;
//	}
//	return Ns_mir;
//}
//
//
//// \******************************************************
//// Eingabe: K_MIR, Alpha_MIR
//// Ausgabe: Nq_MIR, q_MIR, MIR_Daten.AnzahlRuns, MIR._Daten.RunLaenge, Runtime
//// \******************************************************
//int	SA_MIRSolver::ComputeNq_MIR(float Alpha_mir,
//							 float K_mir,
//							 int N1_mir,
//							 int& q_mir,
//							 int& anzahlruns,
//							 int& runlaenge,
//							 float& runtime)
//{
//	int Nq_mir = 0;
//	anzahlruns = runlaenge = q_mir = 0;
//	runtime = 0;
//
//	float Optimum_Prob = pow(K_mir/N1_mir, Alpha_mir);
//
//	if ( myrank == 0 )
//	{
//		if ( Alpha_mir != 0 )
//		{
//			Nq_mir = (int)(FABS(log(Optimum_Prob))*exp(1)*K_mir / Alpha_mir);
////			Nq_mir = (int)(FABS(log(pow(Optimum_Prob,1/Alpha_mir))*exp(1)*K_mir));
////			Nq_mir = (int)(FABS(log(Epsilon*Epsilon))*exp(1)*K_mir / Alpha_mir);
//		}
//		else
//		{
//			std::cerr << "SA_MIRSolver: Nq_MIR can't be calculated" << std::endl;
//			Nq_mir = (int)P->GetLocalN();
//		}
//
//		
//		if ( exp(1)*K_mir != 0 )
//		{
//			q_mir = (int)ceil(Nq_mir/(exp(1)*K_mir));
//		}
//		else
//		{
//			std::cerr << "SA_MIRSolver: q_MIR can't be calculated" << std::endl;
//			q_mir = 1;
//		}
//		
//		if ( q_mir != 0 && Nq_mir !=0 )
//		{
//			anzahlruns = 1+(q_mir-1)/number_of_processors;
//			runlaenge = 1+(Nq_mir-1)/q_mir;
//		}
//		else
//		{
//			anzahlruns = 1;
//			runlaenge = (int)P->GetLocalN();
//		}
//		runtime = anzahlruns * runlaenge* Steptime;
//
//		std::cout << "Nq_MIR = " << Nq_mir << " q_MIR = " << q_mir << std::endl;
//
//		// check
//		std::cout << "it must be zero (two): " <<
//			pow( q_mir*K_mir/Nq_mir, Alpha_mir*q_mir ) -
//			pow( K_mir/N1_mir, Alpha_mir ) << std::endl;
//	}
//	return Nq_mir;
//}
//
//// \******************************************************
////
//// \******************************************************
//void	SA_MIRSolver::TestMIRScheduler(void)
//{
//  	MPI_Barrier(MPI_COMM_WORLD);
//	// Initialloesung erzeugen
//	InitStates();
//
//	std::cout << "SA_MIRSolver: WarmingUp ..." << std::endl;
//	double time_accept=0,time_reset=0;
//	//Setzen der Cooling-Parameter mittels einer Funktion von MIRScheduler
//	((MIRScheduler*)Cooler)->WarmingUp(*P,*State,*BestSolution,MPI_COMM_WORLD);
////	Steptime = (time_accept+time_reset)/2;
//	std::cout << "SA_MIRSolver: WarmingUp completed" << std::endl;
//
//	for( int actual_run_length = Minimum_RunLength;
//		 actual_run_length < Maximum_RunLength;
//		 actual_run_length = (int)ceil(actual_run_length*Betta_Runtime) )
//	{
//		((MIRScheduler*)Cooler)->SetDesiredRunLength(actual_run_length);
//	}
//}
//
//
//// \******************************************************
////
//// \******************************************************
//void SA_MIRSolver::Plot(float B_mir, float Alpha_mir)
//{
//	int i,j;
//	if ( myrank == 0 )
//	{
//		std::ofstream	MIRPlot("pics/MIRPlot"),
//					MIRDeltaCost("pics/MIRDeltaCost");
//
//		for( i = Minimum_RunLength; i < 30000; i=i+2 )
//		{
////			MIRPlot << i << " " << FABS(Estar-Ebest)*pow(K_mir/i,Alpha_mir)+Ebest << std::endl;
//			MIRPlot << i << " " << Ebest_user-B_mir*pow(i,-Alpha_mir) << std::endl;
//		}
//
//		for( j = 0; j < Run_number; j++ )
//		{
//			MIRDeltaCost << RunLength[j] << " " << FABS(AverageCost[j] - Ebest) << std::endl;
//		}
//
////		MIRPlot.close();
////		MIRDeltaCost.close();
//	}
//}
//
//// \******************************************************
////
//// \******************************************************
//void SA_MIRSolver::Plot_Test(double K, float Alpha)
//{
//	int i,j;
//	if ( myrank == 0 )
//	{
//		std::ofstream	MIRPlot("pics/MIRPlot");
//
//		for( i = Minimum_RunLength; i < 30000; i=i+2 )
//		{
//			MIRPlot << i << " " << Ebest_user-pow(K/i,Alpha) << std::endl;
//		}
//	}
//}
//
//
//
//
//
//
//float	SA_MIRSolver::BestOf(float* array,int size)
//{
//	if ( array != 0 && size != 0 )
//	{
//		float best = array[0];
//		for (int i=0; i<size; i++)
//		{
//			if ( Cooler->Better(array[i],best) )
//			{
//				best = array[i];
//			}
//		}
//		return best;
//	}
//	else
//	return 0;
//}
//
//
//float	SA_MIRSolver::AverageOf(float* array,int size)
//{
//	if ( array != 0 && size != 0 )
//	{
//		float Sum = 0;
//		for (int i=0; i<size; i++)
//		{
//			Sum += array[i];
//		}
//		return Sum/size;
//	}
//	else
//	return 0;
//}
//
//
//// \******************************************************
////
//// \******************************************************
//float	SA_MIRSolver::ComputeK_MIR(float Alpha_mir)
//{
//	float K_mir = 0;
//	float k = 0;
//	Cooler->CollectSubchainStatistics(MPI_COMM_WORLD);
//	std::cout << "Cooler->GetRunSigma()= " << Cooler->GetRunSigma() << std::endl;
//	
//	if ( myrank == 0 )
//	{
//		if ( Cooler->GetRunSigma() != 0 )
//		{
//			k = FABS(Estar-Ebest)/Cooler->GetRunSigma();
//			if ( k*k <= 0.5 )
//			{
//				std::cout << "k nicht OK: " << k << std::endl;
//				K_mir = RunLength[Run_number-1];
//			}
//			else
//			{
//				std::cout << "k OK: " << k << std::endl;
//				K_mir = RunLength[Run_number-1]*pow(1-1/(2*k*k),1/Alpha_mir);
//			}
//		}
//		else
//		{
//			std::cerr << "SA_MIRSolver: K_MIR (second) can't be calculated" << std::endl;
//			K_mir = 0;
//		}
//		std::cout << "k= " << k << std::endl;
//		std::cout << "K_MIR (second)= " << K_mir << std::endl;
//	}
//	return K_mir;
//}
//
//
//
//// \******************************************************
////
//// \******************************************************
//int	SA_MIRSolver::ComputeN1_MIR_Derivation(float Alpha_mir, float B_mir)
//{
//	// f(x) = B_mir*x^(-Alpha_mir)+Ebest
//	// f'(x) = -Alpha_mir*B_mir / x^(Alpha_mir+1)
//	// |f'(x)| < Derivation <=> x > (Derivation/(Alpha_mir*B_mir))^(Alpha_mir+1)
//
//	double N1_mir=0;
//	if ( myrank == 0 )
//	{
////		N1_mir = (int)ceil(pow(Derivation/(Alpha_mir*B_mir),Alpha_mir+1));
////		N1_mir = (int)ceil(pow(P->GetLocalN()*Alpha_mir*B_mir/log(Ebest),1/(Alpha_mir+1)));
////		N1_mir = (int)ceil(Derivation*pow(P->GetLocalN()*Alpha_mir*B_mir,1/(Alpha_mir+1)));
//		N1_mir = (int)ceil(pow(P->GetLocalN()*Alpha_mir*B_mir/Derivation,1/(Alpha_mir+1)));
////		N1_mir = (int)ceil(log(P->GetLocalN())*pow(P->GetLocalN()*Alpha_mir*B_mir,1/(Alpha_mir+1)));
//		std::cout << P->GetLocalN()*Alpha_mir*B_mir << "^" << 1/(Alpha_mir+1) << std::endl;
//		std::cout << "N1_mir with Derivation-method is " << N1_mir << std::endl;
//		
//		double EndSolutionValue = B_mir*pow(N1_mir,-Alpha_mir)+Ebest_user;
//		std::cout << "TargetSolValue=" << EndSolutionValue << std::endl;
//		EndSolutionQuality = FABS(EndSolutionValue-Ebest_user)/Ebest_user;
//	}
//	MPI_Bcast(&EndSolutionQuality,1,MPI_FLOAT,0,MPI_COMM_WORLD);
//	return (int)N1_mir;
//}
//
//
//
//// \******************************************************
////
//// \******************************************************
//float SA_MIRSolver::MIR()
//{
//	time_t starttime, endtime;
//	starttime = time(NULL);
//
//std::cout << "Process " << myrank << ": Phase Two begins..." << std::endl;
//	PhaseTwo();
//	Setup();
//std::cout << "Process " << myrank << ": Phase Two completed" << std::endl;
//	
//	if ( myrank == 0 )
//	{
//		Cooler->ShowConfig(std::cout);
//		ShowConfig();	
//	}
//
//	MPI_Bcast(&MIR_Daten,2,MPI_INT,0,MPI_COMM_WORLD);
//	MPI_Bcast(&s_MIR,1,MPI_INT,0,MPI_COMM_WORLD);
//	MPI_Bcast(&K_MIR,1,MPI_INT,0,MPI_COMM_WORLD);
//
//	if ( MIR_Daten.RunLaenge != 0 && MIR_Daten.AnzahlRuns != 0 )
//	{
//		std::cout << "We'll make " << MIR_Daten.AnzahlRuns*number_of_processors
//			 << " runs of length " << MIR_Daten.RunLaenge
//			 << " ; each processor performs " << MIR_Daten.AnzahlRuns
//			 << " runs" << std::endl;
//
//		if ( myrank == 0 )
//		{
//			std::cout << "Estimated Runtime: " << Runtime << " seconds" << std::endl;
//			std::cout << "this are " << Runtime/60 << " minutes" << std::endl;
//		}
//
//		std::cout << "Process " << myrank << ": Phase Three  begins..." << std::endl;
//		float SARunCost = 0;
//		for (int i=0; i < MIR_Daten.AnzahlRuns; i++)
//		{
//			SARunCost += SARun(MIR_Daten.RunLaenge);
//		}
//		SARunCost /= MIR_Daten.AnzahlRuns;
//
//		Cooler->SetNewE(SARunCost);
//		// Bewirkt, dass in Output Mittelwert aller Endkosten aller durchgefuehrten Runs
//		// geschrieben werden
//	}
//	
//	endtime = time(NULL);
//	SATime = endtime - starttime;
//	std::cout << "Process " << myrank << ": Phase Three completed" << std::endl;
//	
//	PostprocTime = time(NULL);	
//	P->Postprocessing(*BestSolution);
//	PostprocTime = time(NULL) - PostprocTime;
//
//	Output();
//	WriteSolution();
//
//	return P->GetCost(*BestSolution);
//}
//
#endif