//  Maurice.Clerc@WriteMe.com
//  Last update: 2006-08-14
// The only purpose of this algorithm is to be a
// _non_ _parametric_ optimiser for unimodal functions.
// (i.e. you have no parameter to tune)
//
//  It looks for the minimum by dichotomy (relaxation method
// one dimension at at time).
// In passing, check whether the function is unimodal or not
// Here it stops in that case, but you may try to continue
//
// WARNING: theoretically, when the algorithm says a function 
//					is not unimodal, it is not completely sure.
//					 It may stop by mistake.
// 					If you have an example, please send it to me
//
//         Conversely, it not always see "low" multimodality.
//				(for example Rosenbrock on [0, 10]^2)
//
//         So, sometimes it keeps searching even on a multimodal function.
//         In such a case, of course, it may find just a local optimum
//		(for Rosenbrock 2D it does find the global one though .... very very slowly)
//
// NOTE: It may not work at all for too flat functions, even unimodal ones, 
// if "accuracy" is too big, and/or "evalMax" too small.
//
//  For SciLab/MatLab
//
clear
global objective
global offset
//-------------------------------------- FUNCTION DEFINITIONS
function f=fit(X, func)
  select func
    case 1 then // Shifted Parabola/Sphere
      f=0;
      for d=1:Dim
        f=f+(X(d)-offset)*(X(d)-offset);
      end

    case 2 then // Partly shifted Parabola/Sphere
      f=0;
      noOffset=1;
      for d=1:Dim
        if noOffset==1 then
          f=f+X(d)*X(d);
        else
          f=f+(X(d)-offset/d)*(X(d)-offset/d);
        end
		noOffset=-noOffset;
      end

    case 3 then // Square root
      f=0;
    	for d=1:Dim
     	 f=f+sqrt(abs(X(d)-offset));
    	end

	case 4 then // Rosenbrock
		f=0;
		for d=1:Dim-1
			xd=X(d)-offset;
			xd1=X(d+1)-offset - xd*xd;
			f=f+100*xd1*xd1+(xd-1)*(xd-1);
		end	

	case 5 then // Rastrigin
		f=0;
		for d=1:Dim
			xd=X(d)-offset;
			f=f+xd*xd -10*cos(2*%pi*xd);
		end
		f=f+10*Dim;

	case 6 then // Schwefel
		f=0;
		for d=1:Dim
			s=0;
			for j=1:d
				s=s+X(j)-offset;
			end
			f=f+s*s;
		end

	case 7 then // 1/x + x -2.  Only with Dim=1 (and x>0)
		f=1/X(1)+X(1)-2;

	case 8 then // CEC 2005 benchmark Sphere (2D at max.)), on [-100,100]^D
		xd=X(1)+3.9311900e+001;
		f=xd*xd;
		if Dim>1 then
			xd=X(2)-5.8899900e+001;
			f=f+xd*xd;
		end
    	f=f-450;

	case 9 then // Benchmark Elliptic (2D at max.), on [-100, 100]^2
			    // Not exactly CEC 2005 one (no rotation)
		xd=X(1)+3.2201300e+001;
		f=xd*xd;
		if Dim>1 then
			xd=X(2)-6.4977600e+001;
			f=f+10^(6/Dim)*xd*xd;
		end
    	f=f-450;

	case 10 then // Alpine
		f=0;
      	for d=1:Dim
           f = f + abs(X(d) * sin(X(d)) + 0.1 * X(d));
      	end


  // Here put some other cases
  //

  end
  f=abs(f-objective);
endfunction

//-------------------------------------------  DICHOTOMY
Dim=2; // Dimension of the search space
for d=1:Dim // Bounds
  xMin(d)=-100;
  xMax(d)=100;
end

func=9 ;// Function
// List (* means "unimodal)
//* 1 => Shifted (on diagonal) Parabola/Sphere
//* 2 => Partly shifted Parabola
//* 3 => Square root For x>=0
//  4 => Rosenbrock.  
//  5 => Rastrigin
//* 6 => Schwefel 
//* 7 => 1/x + x -2  . Only with Dim=1, 
//      and don't use a symmetrical interval (in order to avoid x=0)
//* 8 => CEC 2005 Sphere (2D at max.). Don't forget to set objective=-450
//* 9 => Elliptic (2D at max.) . Objective=-450
// 10 => Alpine

offset=0;
objective=-450;
accuracy=0.000001;
evalMax=1000;
dfMin=accuracy; //Threshold as stop criterion on each dimension
								// and also to check unimodality

// ------------ Dichotomy
evalCurrent=0;
unimodalGlob=Dim; // Hypothesis "unimodal=TRUE"


// First interval [X1, X2]
dFree=1;
for d=1:Dim X1(d)=xMin(dFree); end
X2=X1; X2(dFree)=xMax(dFree);
f1=fit(X1,func); evalCurrent=evalCurrent+1;
f2=fit(X2,func); evalCurrent=evalCurrent+1;

// First middle
X3=(X1+X2)/2;
f3=fit(X3,func); evalCurrent=evalCurrent+1;

// First best (on current dimension)
Xbest=X3;
fBest=f3;

// First best best
XbestBest=Xbest;
fBestBest=fBest;

//---------------------------------------------------------- OPTIMISATION
while evalCurrent<=evalMax & fBestBest > accuracy & unimodalGlob>0 

//dFree

	fBestPrev=fBest;
	df=dfMin+1; // Arbitrary value to begin
	unimodal=1;

	while evalCurrent<=evalMax & df > dfMin & unimodal==1 & fBest>accuracy  // Optimise along dFree
		obvious=0;
		// Sometimes the next interval is obvious
		if f1>f3 & f3>f2 then
			// same X2
			X1=X3;
			f1=f3;
			obvious=1;
		else
			if f1<f3 & f3<f2 then
				// same X1
				X2=X3;
				f2=f3;
				obvious=1;
			end
		end

		if obvious==0 then
	 		// New positions
			X4=(X1+X3)/2;
			X5=(X3+X2)/2;
			f4=fit(X4,func); evalCurrent=evalCurrent+1;
			f5=fit(X5,func); evalCurrent=evalCurrent+1;

			// Check if improvement
			if f4<fBest then Xbest=X4; fBest=f4; end;
			if f5<fBest then Xbest=X5; fBest=f5; end;

	//-- Display
	//	 X1 
  //   X2 
	//		"f1, f4, f3, f5, f2"
	//	 [f1, f4, f3, f5, f2] 
	// --

			// Check unimodality
			if f1<f4 & f3<f4 then unimodal=0; end
			if f4<f3 & f5<f3 then unimodal=0; end
			if f3<f5 & f2<f5 then unimodal=0; end
			if unimodal==0 then unimodalGlob=unimodalGlob-1; end
		
			if unimodalGlob<=0 then 
 				"*** WARNING. Probably multimodal function. I stop here"
			else
				// Prepare the next interval
				if f4<f3 then 
					X2=X3;
					f2=f3;
					// Same X1
				else
					if f5<f3 then
						X1=X3;
						f1=f3;
						// Same X2
					else
						X1=X4;
						X2=X5;
						f1=f4;
						f2=f5;
					end
				end
			end // End of "if unimodalGlob<=0"
		end// End of "if obvious==0"

		// New middle
		X3=(X1+X2)/2;
		f3=fit(X3,func); evalCurrent=evalCurrent+1;
//[f1,f2,f3]

		// Check if improvement
		if f3<fBest then Xbest=X3; fBest=f3; end

		// Update the best best
		if fBest<fBestBest then
			fBestBest=fBest;
			XbestBest=Xbest;
		end

		df=max(abs(f1-f3),abs(f2-f3)); // For stop criterion

//df, fBest, fBestPrev
		fBestPrev=fBest;	
	end // End of Optimisation on one dimension


//Xbest
//fBest
//[f1, f4, f3, f5, f2]

	// Next dimension (cyclicly)
	if dFree==Dim then dFree=1;
	else dFree=dFree+1; 		
	end

	// New interval
	X1=Xbest; X1(dFree)=xMin(dFree);
	X2=Xbest; X2(dFree)=xMax(dFree);
	f1=fit(X1,func); evalCurrent=evalCurrent+1;
	f2=fit(X2,func); evalCurrent=evalCurrent+1;
end // End of "while non stop
evalCurrent, fBestBest, XbestBest