// bibliothèque de fonctions scilab pour l'examen optimisation
// Bruno Pinçon

function [y] = eval_func(func,x)
   // f est soit une fonction, soit une liste contenant 
   // une fonction (de x) suivie de ses paramètres
   //  list( f, p1, p2, p3 )
   if typeof(func) == "function" then
      y = func(x)
   elseif typeof(func) == "list" then
      y = func(1)(x,func(2:$))
   else
      error("=> mauvais argument dans eval_func")
   end
endfunction

function [x1,f1] = recherche_section_doree(x0, f0, d, f, pas, ...
					   tol, tolbis, verb)
//     x0 : point de départ de la recherche
//     f0 : valeur de f en x0
//      d : direction de recherche
//      f : fonction à minimiser dans la direction d (g(t) = f(x0 + t*d))
//    pas : pas initial pour la phase 1
//    tol : tolérance sur l'intervalle [a,c]
// tolbis : tolérance sur les variations de g
//   verb : booléen (vrai pour afficher des informations)
//
// Ce code n'a pas de garde-fou : vous pouvez rajouter un test
// limitant le nombre maximum d'itérations.
//
  tau = 0.5*(1 + sqrt(5))
  tau2 = tau^2
  
  // 1- phase de recherche du triplet (a,b,c)
  if verb then printf("\n recherche lin phase 1"), end
  a = 0; ga = f0; c = pas; gc = eval_func(f,x0 + c*d);
  if gc >= ga then
     while %t
	b = c/tau2; gb = eval_func(f,x0 + b*d);
	if gb < ga then, break, end
	c = b; gc = gb
     end
  else // gc < ga
     b = c; gb = gc; 
     while %t
	c = a + (b-a)*tau2; gc = eval_func(f,x0 + c*d);
	if gc > gb then, break, end
	a = b; ga = gb
	b = c; gb = gc
     end
  end
  
  // 2 - phase de réduction de l'intervalle par la section dorée
  if verb then printf("\n recherche lin phase 2"), end
  while  c-a > tol  |  max(ga,gc)-gb > tolbis*max(1,abs(gb))
     delta = c + a - b
     gdelta = eval_func(f,x0 + delta*d);
     if gdelta < gb then
	if delta < b then, c = b; gc = gb; else, a = b; ga = gb; end
	b = delta; gb = gdelta
     else
	if delta < b then, a = delta; ga = gdelta; else, c = delta; gc = gdelta; end
     end
  end
  
  // mise à jour sortie
  x1 = x0 + b*d; f1 = gb
endfunction


function [xopt, fopt, gopt, iter] = gradient_conjugue(x0, f, grdf, pas, ...
						      tol, itermax, verb)
   iter = 0
   xopt = x0; fopt = eval_func(f,x0); gopt = eval_func(grdf,x0); d = -gopt; 
   norm2_gopt = gopt'*gopt; 
   
   while sqrt(norm2_gopt) > tol  & iter < itermax then
      [xopt,fopt] = recherche_section_doree(xopt, fopt, d, f, pas, ...
					    1e-9, 1e-11, %f)
      gopt_new = eval_func(grdf,xopt)
      Beta = (gopt_new - gopt)'*gopt_new / norm2_gopt
      gopt = gopt_new; norm2_gopt = gopt'*gopt; 
      d = -gopt + Beta*d
      iter = iter + 1
      if verb then
	 printf("\n iteration = %d", iter) 
	 printf("\n  ||g(x)|| = %g",  sqrt(norm2_gopt))	 
	 printf("\n      f(x) = %g", fopt)
      end
   end
endfunction

function [y] = gauss(x,p)
   // mélange de 2 gaussiennes
   // x peut être un vecteur et alors y est un vecteur
   // de même format, la fct étant évaluée composante par
   // composante
   alpha = p(1)
   mu1 =  p(2)
   sigma1 = p(3)
   mu2 = p(4)
   sigma2 = p(5)
   coef = 1/sqrt(2*%pi);
   y = coef*(  (1-alpha)*exp(-0.5*((x-mu1)/sigma1).^2)/sigma1 ...
	        +  alpha*exp(-0.5*((x-mu2)/sigma2).^2)/sigma2 )
endfunction
 
function [dy] = dgaussdp(x,p)
   // gradient du mélange de 2 gaussiennes (sous forme colonne)
   // x peut être un vecteur COLONNE et alors y est une matrice
   //  5 x length(x), la colonne i correspondant au gradient
   //  (colonne) en x(i)
   [m,n] = size(x)
   if n ~= 1 then
      error(" dgaussdp : le 1 er argument doit être un vecteur colonne")
   end
   alpha = p(1)
   mu1 =  p(2)
   sigma1 = p(3)
   mu2 = p(4)
   sigma2 = p(5)
   coef = 1/sqrt(2*%pi);
   Q1 = exp(-0.5*((x-mu1)/sigma1).^2)
   Q2 = exp(-0.5*((x-mu2)/sigma2).^2) 
   dy = coef* [-Q1/sigma1 + Q2/sigma2, ...
	      (1-alpha)*Q1.*(x-mu1)/sigma1^3,...
	      (1-alpha)*Q1/sigma1^2 .*(((x-mu1)/sigma1).^2 -1),...
	          alpha*Q2.*(x-mu2)/sigma2^3,...
	          alpha*Q2/sigma2^2 .*(((x-mu2)/sigma2).^2 -1)]'
endfunction

function [x,y] = charger_donnees()
   datas =  fscanfMat("datas.exam")
   x = datas(:,1)
   y = datas(:,2)
endfunction

//function [popt, fopt, gopt, iter] = gauss_newton(p0, f, grdf, tol, itermax, ...
//						 x, y)
//   iter = 0
//   popt = p0; 
//   fopt = f(p0,x,y); 
//   gopt = grdf(p0,x,y); 
//   
//   while norm(gopt) > tol  & iter < itermax then
//      // calcul de l'approximation de la hessienne
//      ................
//      ................
//      // calcul de la direction de descente
//      ................
//
//   
//      // recherche unidimentionnelle (commencant par le
//      // pas de Newton cad 1, pas divise par 2 tant que
//      // la fonction de decroit pas)
//      tau = 1;
//      while %t
//	     pnew = popt + tau*d;
//	     fnew = f(pnew,x,y)  
//	     if fnew < fopt then
//	         break
//	     else
//	         tau = tau/2
//	     end
//      end
//      popt = pnew; fopt = fnew
//      gopt = grdf(popt,x,y)
//      iter = iter + 1
//   end
//endfunction