// Pierre L. Douillet
// www.douillet.info
// version 26

// ligne 526 : test mauvais conditionnement
// utiliser better2 (50, 100)

ident='planx.26';
if MSDOS then 
  cd "c:/MapleV6/Users/Douillet/planx" ; baserep=pwd(); 
  mydir=baserep+"/"; myfig=baserep+"/"; 
else
  cd "~/docs/Cherche/doe" ; baserep=pwd(); 
  mydir=baserep+"/datas/";
  myfig=baserep+"/figures/";
end//if

global shortnam datas ix ifx jx lx moda siza zona Omega 
global ma maa mb mc ms mx mdelta Datas mA mS mX mDelta fun idents

// Une fonction peut lire les variables de la fonction appelante (scoped)
// Au moment où une fonction va écrire, la variable devient locale
// Pour qu'une fonction puisse écrire dans une variable de la fonction
// appelante, il faut que les DEUX aient déclaré la variable comme globale
// OU BIEN utiliser le mécanisme clear/return (scoping)

//*********************************************************
//* menus
//*********************************************************

function menu
  while %t do
    item=x_choose(['read','deal','pick','optimize'],..
        [ident+ ' main menu'])
    select item
    case 0 then return
    case 1 then menu_read(); 
    case 2 then menu_deal();
    case 3 then sig=x_mdialog('random search',['nb of trials'],['30']);
      if sig ~= [] then iterplan(eval(sig(1))); end//if
    case 4 then sig=x_mdialog('directed search',  ['nb of blocs'],['15']);
      if sig ~= [] then better(eval(sig(1))); end//if
    end//select
  end//while
endfunction

function menu_deal
global fun maa
  while %t do
  item=x_choose(['mambms','draw','map','plages', 'xy+yz', 'xx+xy', 'lenteur', 'plot3d'],..
      ['deal menu'])
  select item
  case 0 then return
  case 1 then mambms();  
  case 2 then draw_residus(); draw_influences(); 
  case 3 then sig=x_mdialog('factors to map',['x','y'],['1','2']);
    disp(recmap([eval(sig(1))], [eval(sig(2))], datas)); printf("");
  case 4 then plages();
  case 5 then
    printf('\nxy+yz\n')
    maa=wcenter(eval(datas(2:$,1:$-1)), 'r');
    for u=1:ix do 
      for v=u+1:ix do aug(u,v); end//for
    end//for
    mAmBmS ([ones(jx,1), maa]);
  case 6 then
    printf('\nxx+xy\n')
    traiter67(eval(datas(2:$,1:$-1)));
  case 7 then
    traiter67([eval(datas(2:$,1)) .^(-1), eval(datas(2:$,2:3)) ]);
  case 8 then sig=x_mdialog('response surface',['grid','s'],['30','1']);
    drawfig(eval(sig(1)), eval(sig(2)) )
  end//select
  end//while
endfunction

//*********************************************************
//* reading data
//*********************************************************

function rep=menu_read(quiquoi)
global describe shortnam datas ifx idents mc

// nam                    shortnam       french                        
describe=[
  "experimentique-fr.txt" "xperfr"     "xperfr"
  "experimentique-us.txt" "xperus"     "xperfr"
  "caoutchouc.txt"        "caout"      "none"
  "claire.txt"            "clair"      "clair"
  "nist-eddy.was.txt"     "eddy"       "eddy"
  "nine1.txt"             "nine1"       "none"
  "nine3.txt"             "nine3"       "none"
  "eval07a.txt"           "eval07a"    "eval07"
  "baie03.txt"            "baie3"      "baie"
  "tess.csv"              "tess"       "none"
  "TestColoris06.csv"     "color"      "color"
  "fillitude.txt"         "fille"      "fille"
  "pharma.txt"            "pharma"     "pharma"
];
describe_old=[
  "eval07b.txt"           "eval07b"    "eval07"  
  "eval07c.txt"           "eval07c"    "eval07"
  "baie01.txt"            "baie1"      "baie"
  "baie02.txt"            "baie2"      "baie"
  "waeghe.csv"            "waeghe"     "none"
]

describe=gsort(describe, 'lr', 'i');
if isdef('quiquoi') then 
  item=vectorfind(describe(:,1), quiquoi,'r'); 
  if item == [] then item=0 ; end//if
else
  item=x_choose(describe(:,2), ['read menu'])
end//if

if item == 0 then rep=[]; return; end//if
nam=describe(item,1);
shortnam=describe(item,2);
french=describe(item,3);

// lire les lignes du fichier (en tant que chaînes de caractères)
fich = mydir + nam ; printf('source: %s\n',fich)
hndl = mopen(fich, "r") ; rep = mgetl(hndl) ; mclose(hndl) ;

// mise en forme préalable des données
// ... cela change selon les fichiers !!!!
select french
case "xperfr" then 
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, '°C', '') ;
  // #x0bc !!!  et pas unicode
case "clair" then
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, ',', '.') ; rep = strsubst( rep, ';', ' ') ;
case "eval07" then
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, '_;_', ' ') ; 
case "baie" then
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, ',', '.') ;  rep = strsubst( rep, ';', ' ') ;
  rep = strsubst(rep, 'With_yellow_fabric_in_API', 'yellow');
case "color" then
  rep = strsubst( rep, ' ', '_') ;
  rep = strsubst(rep, ascii(9)+ascii(9), ascii(9)+'999'+ascii(9));
  rep = strsubst(rep, ascii(9)+ascii(9), ascii(9)+'999'+ascii(9));
end//select


// format et test de longueur du format
toks=tokenpos(rep(1)) ; toksx=size(toks,'r') ; // [" ", ascii(9)]
frmt=''; for j=1:68 do frmt=frmt+'%s'; end//do 
frmtx=size(strindex(frmt,'%'),'c') ;
if frmtx < toksx then error(sprintf("too short format : %d < %d", frmtx, toksx)); end//if;

// interprétation des lignes de texte (donne une matrice)
datas = msscanf(-1, rep, frmt) ;
if size(datas,'r')==1 then disp(rep(1:3)); error 'something happens' ; end//if

ifx=0; // le nombre de paramètres à valeurs continues (placés d'abord)
select french
case "eddy" then
// traitement spécial : les résultats sont en colonne 1, le # en dernière
  tmp=['0';datas(2:$,$)]; [qqq,trad]=gsort(eval(tmp),'r','i') //row, increasing
  datas=datas(trad,[2,3,4,1]) 
case "baie" then
// traitement spécial : le # en col3
  tmp=['0';datas(2:$,3)]; [qqq,trad]=gsort(eval(tmp),'r','i') //row, increasing
  datas=datas(trad,[5,6,9,12]) 
case "color" then
  datas=datas(:,[2, 7,8,9,10, 3,4,5,6, 11]);
  datas=datas(find(and(datas<>'999','c')), :);
  ifx=5
case "fille" then
  idents=datas(2:$,1)
//  datas=datas(:,[2,3,4,6]); ifx=3;
  datas=datas(:,[2,3,4,5,6]); ifx=4;
case "pharma" then
  //objectif double
  mc=datas(2:$,$); datas=datas(:,1:$-1);
end//select

printf('%s  ', datas(1,:)'); printf('\n\n') 

readmoda

endfunction

function readmoda
global ix ifx jx lx moda siza zona Omega
  // la première ligne est une ligne de titres
  // la dernière colonne est une colonne de résultats
  [jx,ix]= size(datas); ix=ix-1; jx=jx-1 // décrémente ix et jx !!!
  // ifx est positionné par french dans menu_read
  
  // liste des modalités
  moda=list() ; siza=1:ix ; zona=1:ix+1 ;
  for i =1:ifx do
   moda(i)= size(unique(datas(2:$,i)),'*') ;
   siza(i)= 1;
   zona(i+1)=zona(i)+siza(i);  
  end//for
  for i =ifx+1:ix do
   moda(i)= sort(unique(datas(2:$,i)))' ;
   siza(i)= size(moda(i), 'c')-1;
   zona(i+1)=zona(i)+siza(i);
  end//for
  zona=[zona(1:$-1)+1;zona(2:$)]
  lx=sum(siza)+1 ; Omega=prod(siza+1)// les différentes expériences possibles

endfunction

//-------------- gestion des modalités, codage/decodage -----------------------
// txt = ligne de texte
// itm = ligne de numéros, de 0 à  siza(i)
// num = numéro unique, de 0 à  Omega-1
// cod = le codage en variables indépeNdantes 

function itm=txt2itm(txt)
  jj=size(txt,'r');  itm=ones(jj,ix)*17;
  for j=1:jj do
    for i =1:ifx do
      itm(j,i)= eval(txt(j,i)) ;
    end//for
    for i =ifx+1:ix do
      ici=vectorfind(moda(i), txt(j,i), 'c');
      if ici == [] then error (sprintf("==> : undefined modality %s for_factor %s", txt(jj,i), datas(1,i) )); end//if
      itm(j,i)= ici-1 ;
    end//for
  end//for;
endfunction

function txt=itm2txt(itm)
  jj=size(itm,'r')
  txt= string(itm);
  for j=1:jj do 
    for i =1:ifx    do txt(j,i)=string(itm(j,i))    ; end//for
    for i =ifx+1:ix do txt(j,i)=moda(i)(itm(j,i)+1) ; end//for
  end//for
endfunction

function num = itm2num(itm)
  if ifx > 0 then error (sprintf('attempt to itm2num while ifx= %d', ifx)) ; end//if
  jj=size(itm,'r'); num=zeros(jj,1)
  for j=1:jj do 
    for i=ix:-1:1 do num(j) = num(j)*(1+siza(i))+itm(j,i); end//for
  end//for
endfunction

function [q,r] = euclid(a,b)
  q= floor(a/b); r=a-b*q
endfunction

function itm = num2itm(num) //num est un vecteur colonne
  if ifx > 0 then error (sprintf('attempt to num2itm while ifx= %d', ifx)) ; end//if
  itm=[] ; 
  for i =1:ix-1 do [num, tmp]= euclid(num, 1+siza(i)); itm=[itm, tmp]; end//for
  itm=[itm,num]
endfunction

function codd=itm2cod(itm)
  jj=size(itm,'r'); codd=[];
  for j=1:jj do
    cod=[1,itm(j,1:ifx)] ; 
    for i=ifx+1:ix do
    ici=itm(j,i); if ici==0 then co=-ones(1,siza(i))
      else co=[zeros(1,ici-1),1,zeros(1,siza(i)-ici)]; end//if
    cod=[cod, co]; end//for
  codd=[codd;cod]; end//for
endfunction

function itm = cod2itm(cod)
  jj=size(cod,'r')
  itm=zeros(jj, ix);
  for j=1:jj do 
    for i=1:ifx do
      itm(j,i)=cod(j,i)
    end//for
    for i=ifx+1:ix do
      itm(j,i)=max(0,cod(j,zona(1,i):zona(2,i))*(1:siza(i))')
    end//for
  end//for
endfunction

function rr=recmap(tox,toy,DAT)
// DAT est soit une matrice of strings (=datas)
//         soit une matrice d'items
  txy=[tox,toy]
  fac=zeros(1,ix); tmp=prod(siza(txy)+1);
  for i=txy do tmp = tmp/(1+siza(i)); fac(i)=tmp; end//for

  if type(DAT(1))==10 then
    ptn= fac*txt2itm(DAT(2:$,:))';
  else
    ptn= fac*DAT';
  end//if
  rr=zeros(prod(siza(toy)+1), prod(siza(tox)+1)) ;
  for j=ptn do rr(j+1)=rr(j+1)+1; end//for
endfunction


//-------------- mise en oeuvre-------------------------------------------

// calculs de base
function frv=mambms ()
global ma mb ms mx mdelta
  ma=itm2cod(txt2itm(datas(2:$,:))) ; mb= eval(datas(2:$, $)) ;
  if mb == zeros(jx,1) then acp() ; frv=1; return ; end//if
  ms = ma' * ma ; mx= (1/ms) * ma' * mb ; mdelta= ma*mx-mb ;
  redvar= variance(mdelta); rawvar= variance(mb);
  if (lx==jx) | (abs(redvar)<%eps) then frv= %nan
   else frv= rawvar/redvar*(jx-lx)/(jx-1);
  end//if
  printf('lx= %2d, rawvar= %f, redvar= %f, est(FRV)= %f\n', lx, rawvar, redvar, frv)
endfunction

function FRV=mAmBmS (maa)
global mA mb mS mX mDelta
  mA=maa ; mb= eval(datas(2:$, $)) ;
  mS = mA' * mA ; mX= (1/mS) * mA' * mb ; mDelta= mA*mX-mb ;
  Redvar= variance(mDelta); rawvar= variance(mb);
  Lx=size(mA,'c'); if (Lx==jx) | (abs(Redvar)<%eps) then FRV= %nan
   else FRV= rawvar/Redvar*(jx-Lx)/(jx-1);
  end//if
  printf('Lx= %2d, rawvar= %f, Redvar= %f, est(FRV)= %f\n', Lx, rawvar, Redvar, FRV)
endfunction

//-------------- acp --------------------------------------------------

function acp ()
global ma mA mS 
ma=ma(:,2:$); //on vire la colonne de 1
mA=wcenter(ma,'r'); mS=mA'*mA/(jx-1); 
[pp,dd]=spec(mS); vec1=pp(:,$); vec2=pp(:,$-1); vp1=dd($,$); vp2=dd($-1,$-1);

xxx=mA*vec1; yyy=mA*vec2

curfig= scf(2); clf();
curfig.figure_size=[500, 500] ;
plot2d([min(xxx),max(xxx)],[0,0])
plot2d([0,0],[min(yyy),max(yyy)])
fac=0.1;
for j=1:jx do
 xstring(xxx(j)+fac*(2*rand()-1),yyy(j)+fac*(2*rand()-1), idents(j))
end//do

qq=pp*dd^0.5; 
fac=max(xxx)/0.9; 

for i=1:ix do
 xstring(qq(i,ix)*fac, qq(i,ix-1)*fac, datas(1,i))
 curax=gca();
 curax.children(1).font_foreground=5;
end//do
zangle=-%pi:%pi/20:+%pi;
plot2d(cos(zangle)*fac, sin(zangle)*fac)
curax=gca();
curax.children(1).children.foreground=5;

curax=gca();
curax.sub_ticks=[0,0] ;
curax.isoview='on';

endfunction

//-------------- les dessins --------------------------------------------------

function draw_residus
  curfig= scf(1); clf();
  curfig.figure_size=[1000, 500] ;
  
  plot(mb, mdelta, 'xk')
  
  curax=get('current_axes') ;
  curax.children.children.thickness=3 ;
  curax.x_location='middle' ;
  curax.box='off' ;
  curax.sub_ticks=[0,0] ;
  curax.font_size=4 ;
  curax.title.font_size=4 ;
  curax.data_bounds(3)=curax.data_bounds(3)*1.02 ; // pour faire apparaitre le dernier point
  xtitle('discrepancies versus values')
  try
    curax.margins=[0.05,0.01,0.15,0.05] ;
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(1, myfig+shortnam+"_residus")
  end//try
endfunction


// influences
function draw_influences
  curfig= scf(2) ; clf();
  curfig.figure_size=[1000, 500] ;

  base=2 ; xti=ones(1:ix)' ; influences=[] ;
  for i=1:ix do
  msg=mx(base:base+siza(i)-1);
  msg=[-sum(msg);msg] ; //; special
  influences(i)=max(msg)-min(msg) ;
  plot([base+2*i:base+2*i+siza(i)],msg) ; plot([base+2*i:base+2*i+siza(i)],msg, 'o')
  xti(i)= base+2*i+siza(i)/2 ;
  base=base+siza(i) ;
  end//for

  curax=get('current_axes') ;
  curax.box='off' ;
  curax.sub_ticks=[0,0] ;
  xtidat=strsubst(datas(1,1:$-1)', 'Presse-étoupe', 'P-étoupe');
  curax.x_ticks=tlist(curax.x_ticks(1), xti, xtidat) ; 
//  curax.data_bounds(1,1)=3 ; //marge à gauche
  curax.data_bounds(1,1)=curax.data_bounds(1,1)-1 ;
  curax.data_bounds(2,1)=curax.data_bounds(2,1)+1 ;
  curax.data_bounds(1,2)=curax.data_bounds(1,2)*1.05 ;
  curax.data_bounds(2,2)=curax.data_bounds(2,2)*1.05 ;
  curax.tight_limits="on" ;

  curax.x_location='bottom' ;   curax.y_location='left' ;
  xtitle('influences : '+sprintf('%4.2f  ', influences)) ;
  curax.font_size=4 ;   curax.title.font_size=4 ;
  try
    curax.margins=[0.07,0.03,0.13,0.1] ;
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(2, myfig+shortnam+"_influ")
  end//try
endfunction


//-------------- le cas continu -----------------------------------------------

function aug(u,v)
  global maa 
  maa=[maa, maa(:,u) .* maa(:,v)]
endfunction

function plages()
    if mA==[] then printf('plages : only in the continuous case\n'); continue ; end//if
    printf('%12s', datas(1,1:$-1)');printf('\n');
    printf('%12f', max(mA(:,[2:ix+1]),'r')')   ;printf('\n');
    printf('%12f', min(mA(:,[2:ix+1]),'r')')   ;printf('\n\n');
endfunction

function traiter67(maaX);
  global maa fun
  maa=wcenter(maaX, 'r');
  for u=1:ix do 
    for v=u:ix do aug(u,v); end//for
  end//for
  mAmBmS ([ones(jx,1), maa]);
  modele='y=%f+%f*s+%f*t+%f*a+ (s.*(%f*s+%f*t+%f*a))+(t.*(%f*t+%f*a))+%f*a.*a';
  azor=msprintf( strsubst(modele, '%f', '%16.12f'), mX')
  deff('y=fun(s,t,a)', azor)
  optima(maaX)
endfunction

function [x,lagr,f]=optima(maaX)
  clear('ppp'); clear('qqq');
  ppp=mX(2:ix+1)
  qqq=zeros(ix,ix); offset=ix+1;
  for j=1:ix do
    offset=offset-j+1;
    for k=1:ix do
      if j>k then qqq(j,k)=qqq(k,j)
      else qqq(j,k)=mX(offset+ix*(j-1)+k)
      end//if
    end//for
    qqq(j,j)=2*qqq(j,j)
  end//for
  top=1.4;
  command='[x,lagr,f]=quapro(-qqq,-ppp,[],[],-top*ones(ix,1),+top*ones(ix,1))';
  //printf('%s\n\n', command);
  execstr(command); f=mX(1)-f;
  xx=mean(maaX,'r')+stdev(maaX,'r') .* x';
  printf('%12.5f ',  x ); printf('%12.5f\n',f);
  printf('%12s ', datas(1,:)'); printf('\n');
  printf('%12.5f ', xx'); printf('%12.5f\n\n',f);
  [ppp,qqq]=return(ppp,qqq);
endfunction

function drawfig(qq, mys)
  curfig= scf(3) ; clf();
  curfig.figure_size=[800, 800] ;
  curfig.color_map=jetcolormap(32); 

  top=1.4
  xx=ones(qq,1)*linspace(-top,+top,qq);
  yy=linspace(-top,+top,qq)'*ones(1,qq);
  ss=ones(qq,qq)*mys;
  zz=fun(ss,xx,yy);
  surf(xx,yy,zz,'facecol','interp');
  curax=gca(); curax.rotation_angles=[0,-90];
  xtitle(sprintf('file= %s, R%s= %f, x= R%s, y= R%s', shortnam, datas(1,1), mys, datas(1,2), datas(1,3)))
endfunction

//-------------- selection d'un design ----------------------------------------

function [Datas,mA,mS]=randplan()
  nums=zeros(jx,1);
  while %t do
    lesj=0:Omega-1; 
    for j=1:jx do
      loc=grand(1,1,'uin',1,Omega+1-j); nums(j)=lesj(loc); lesj=[lesj(1:loc-1),lesj(loc+1:$)];
    end//for
    itms=num2itm(nums); mA=itm2cod(itms); mS=mA'*mA;
    if det(mS)>0.5 then Datas= [datas(1,1:ix); itm2txt(itms)]; break; end//if
  end//while
endfunction

function iterplan(tx, varargin)
global Datas mA mS
  memscore=-1; tic()
  for t=1:tx do [dat,mAA,mSS]=randplan(); score=min(spec(mSS));
    if score>memscore then memscore=score, Datas=dat, mA=mAA, mS=mSS; end//if
  end//for
  printf("\n Datas, mA, mS generated in %d essais, laps= %6.4f sec., min(spec(mS))= %6.4f\n", tx, toc(), memscore)
  if size(varargin)>0 then mAmBmS(mA); end//if
endfunction

// scoped mT
// global mA mS   (déjà  déclarés)   

function mT=tous()
  dT=(0:siza(1))'; 
  for i=2:ix do
    toux0=dT; jj=size(dT,'r'); dT=[]; 
    for k=0:siza(i) do dT=[dT;[toux0,k*ones(jj,1)]] ; end//for
  end//for
  mT=itm2cod(dT);
endfunction

function [err,qui]=maxerr_old(mS)
//scoped mT
  if det(mS)<0.00001 then [err,qui]=(%inf,%nan); return ; end//if
  ims=(1/mS)
  gr=[1:Omega]; for k=1:Omega do
    lig=mT(k,:); gr(k)=lig*ims*lig'
  end//for
  [err,qui]=maxi(gr)
endfunction

function [err,ran,qui]=randerr(mSS)
// err, un lieu de err, tous les err
//scoped mT
  if det(mSS)<0.5 then [err,ran,qui]=(%inf,%nan,%nan); return ; end//if
  ims=(1/mSS)
  gr=[1:Omega]; for k=1:Omega do
    lig=mT(k,:); gr(k)=lig*ims*lig'
  end//for
  err=maxi(gr)
  qui=vectorfind(gr,err,'c')
  quix=grand(1,1,'uin',1,size(qui,'*')); ran=qui(quix)
endfunction

function score=stepbetter()
global mA mS
//scoped mT
  mAA=mA; 
  while %t do [score, nouv]= randerr(mS);
    if ~isnan(nouv) then break; end//if
  end//while
  score0=score; mAug=[mA;mT(nouv,:)];
  mtry=mAug(2:$,:); scotry=randerr(mtry'*mtry);
  if scotry<score then score=scotry; mAA=mtry; end//if
  for j=2:jx do
  mtry=[mAug(1:j-1,:);mAug(j+1:$,:)]; scotry=randerr(mtry'*mtry);
   if scotry<score then score=scotry; mAA=mtry; end//if
  end//for
  if score <> score0 then
    mA=mAA; mS=mA'*mA; score= stepbetter() 
  end//if
endfunction

function score=stepbetter0()
global mA mS
//scoped mT
  mAA=mA; 
  while %t do [score, nouv]= randerr(mS);
    if ~isnan(nouv) then break; end//if
  end//while
  score0=score; mAug=[mA;mT(nouv,:)];
  mtry=mAug(2:$,:); scotry=randerr(mtry'*mtry);
  if scotry<score then score=scotry; mAA=mtry; end//if
  for j=2:jx do
  mtry=[mAug(1:j-1,:);mAug(j+1:$,:)]; scotry=randerr(mtry'*mtry);
   if scotry<score then score=scotry; mAA=mtry; end//if
  end//for
  if score <> score0 then
    mA=mAA; mS=mA'*mA; 
  end//if
endfunction


function [qo,loc]= better(tx)
global mA mS Datas
//scoping mT
  tic(); clear('mT'); mT=tous(); qo=%inf; msmo=-%inf; loc=0;
  for q=1:tx do
    [Datas,mA,mS]=randplan(); score=stepbetter(); msm=min(spec(mS)); 
    if (score<qo) | ((score==qo)&(msm>msmo)) then 
      printf("%3d----------------------------- %f -- %f\n",q,score,msm); qo=score; qa=mA; qs=mS; msmo=msm; loc=q;
    // else  printf("%3d-- %f\n",q,score);
    end//if
  end//for
  mA=qa; mS=qs; Datas=[datas(1,1:$-1);itm2txt(cod2itm(mA))] 
  printf(">>>-------------------randerr(mS)= %f\n",qo);
  printf(">>>-------------------randerr(ms)= %f\n",randerr(ms));
  printf("\n Datas, mA, mS generated in %d essais, laps= %6.4f sec., min(spec(mS))= %6.4f\n", tx, toc(), min(spec(mS)) )
mT=return(mT)
endfunction

function [qscore, loc]=better2(tx, qx)
global mA mS Datas 
//scoping mT
  tic(); clear('mT'); mT=tous();
  qscore=%inf;  qmsm=-%inf;
  for q=1:qx do
    [Datas,mA,mS]=randplan(); score=randerr(mS); // mS est inversible
    for t=1:tx do
      kill=grand(1,'uin',1,jx);
      mRed=[mA(1:kill-1,:);mA(kill+1:$,:)];
      
      [qqq, nouv]=randerr(mRed'*mRed);
      if isnan(nouv) then qdet=0; 
        for it=1:Omega do tmp=[mRed;mT(it,:)]; idet=det(tmp'*tmp);
         if idet>qdet then qdet=idet; nouv=it; end//if
        end//do
      printf("%d %d ======= nan encountered ; nouv=%d\n", q,t,nouv)  
      end//if      

      mAA=[mRed;mT(nouv,:)]; scotry=randerr(mAA'*mAA);
      if scotry<score then
        mA=mAA; mS=mA'*mA; score=scotry;
        //printf("%3d-- %f\n",t,score);
        // vieux printf("%3d----------------------------- %f\n",t,score);
        // vieux else printf("%3d-- %f\n",t,scotry);
      end//if
    end//do
    msm=min(spec(mS));
    if score < qscore  | ((score==qscore)&(msm>qmsm))  then 
      qA=mA; qS=mS; qscore=score; qmsm=msm; loc=q;
      printf("%3d----------------------------- %f -- %f\n",q,qscore,qmsm);
      // else printf("%3d-------- %f\n",q,score);
    end//if
  end//do
  mA=qA; mS=qS;
  Datas=[datas(1,1:$-1);itm2txt(cod2itm(mA))]; 
  printf(">>>-------------------randerr(mS)= %f\n",randerr(mS));
  printf(">>>-------------------randerr(ms)= %f\n",randerr(ms));
  printf("\n Datas, mA, mS generated in %d essais, laps= %6.4f sec., min(spec(mS))= %6.4f\n", tx*qx, toc(), min(spec(mS)) )
mT=return(mT)
endfunction


function loop_better
global jx
 resu=mopen('loop_resu_'+shortnam+'.txt','w');

 for jx=lx:13 do [qo,loc]=better(300)
 fprintf(resu, "%2d\t%6.4f\t%6.4f\t%d\t%d\n", jx, qo, min(spec(mS)), trace(mS),loc );
 end//for

 for jx=14:15 do [qo,loc]=better2(100,300)
 fprintf(resu, "%2d\t%6.4f\t%6.4f\t%d\t%d\n", jx, qo, min(spec(mS)), trace(mS),loc );
 end//for

 for jx=16:16 do [qo,loc]=better2(150,300)
 fprintf(resu, "%2d\t%6.4f\t%6.4f\t%d\t%d\n", jx, qo, min(spec(mS)), trace(mS),loc );
 end//for

 for jx=17:Omega do [qo,loc]=better2(150,200)
 fprintf(resu, "%2d\t%6.4f\t%6.4f\t%d\t%d\n", jx, qo, min(spec(mS)), trace(mS),loc );
 end//for
 mclose(resu);
endfunction //loop_better

function loop_better_claire
global jx
 resu=mopen('loop_resu_'+shortnam+'.txt','w');

 for jx=lx:lx+1 do [qo,loc]=better(300)
 fprintf(resu, "%2d\t%6.4f\t%6.4f\t%d\t%d\n", jx, qo, min(spec(mS)), trace(mS),loc );
 end//for

 for jx=lx+2:Omega/2 do [qo,loc]=better2(100,300)
 fprintf(resu, "%2d\t%6.4f\t%6.4f\t%d\t%d\n", jx, qo, min(spec(mS)), trace(mS),loc );
 end//for

 mclose(resu);
endfunction //loop_better

function draw_better
global bet
  clf();
  resu=mopen('loop_resu_'+shortnam+'.txt','r'); 
  req = mgetl(resu) ; mclose(resu) ;
  bet = eval(msscanf(-1, req, "%s%s%s%s")) ;
  xxx=bet(:,1);
  yyy=bet(:,2).*xxx;
  zzz=(bet(:,3).^(-1)).*xxx;
  fac=min(yyy)/min(zzz); zzz=zzz*fac; 
  plot2d(xxx,yyy);
  plot2d(bet(:,1),zzz,5);
  plot2d(xxx, min(yyy)*ones(xxx),2)
endfunction


//-----------------------------------------------------------------

function rapport
global fun maa
  menu_read()
  mambms()

  maa=wcenter(eval(datas(2:$,1:$-1)), 'r');
  for u=1:ix do 
    for v=u+1:ix do aug(u,v); end//for
  end//for
  mAmBmS ([ones(jx,1), maa]);

  maaX=eval(datas(2:$,1:$-1)); maa=wcenter(maaX, 'r'); 
  for u=1:ix do 
    for v=u:ix do aug(u,v); end//for
  end//for
  mAmBmS ([ones(jx,1), maa]);
  modele='y=%f+%f*s+%f*t+%f*a+ (s.*(%f*s+%f*t+%f*a))+(t.*(%f*t+%f*a))+%f*a.*a'; 
  deff('y=fun(s,t,a)', msprintf( strsubst(modele, '%f', '%16.12f'),mX'))
  optima(maaX)

  maaX=[eval(datas(2:$,1)) .^(-1), eval(datas(2:$,2:3)) ];maa=wcenter(maaX, 'r'); 
  for u=1:ix do 
    for v=u:ix do aug(u,v); end//for
  end//for
  mAmBmS ([ones(jx,1), maa]);
  modele='y=%f+%f*s+%f*t+%f*a+ (s.*(%f*s+%f*t+%f*a))+(t.*(%f*t+%f*a))+%f*a.*a'; 
  deff('y=fun(s,t,a)', msprintf( strsubst(modele, '%f', '%16.12f'),mX'))
  optima(maaX)
endfunction

function res=mat2tex(ma,varargin)
tmp=texprint(ma)
tmp=strsubst(tmp,'{\pmatrix{','')
tmp=strsubst(tmp,'}}','')
deb='$$\left(\begin{array}{';
for j=1:size(ma,'c') do deb=deb+'c'; end//for
deb=deb+'}'
res=[deb;tmp;'\end{array}\right)$$']
if size(varargin)>0 then
  fich=myfig+varargin(1)+'.tex'
  disp(fich)
  hndl=mopen(fich, 'w')
//  mputstr(res(1),hndl) ; mputstr(res(2),hndl) ; mputstr(res(3),hndl)
  fprintf(hndl, "%s\n%s\n%s", res(1),res(2),res(3))
  mclose(hndl)
end//if

endfunction

function xchgplan(tx)
global mA mS Datas 
//scoping mT
  tic(); clear('mT'); mT=tous();
  bluk=zeros(jx,Omega);
  tscore=%inf; tmsm=-%inf;
  for t=1:tx do
    [Datas,mA,mS]=randplan(); 
    score=randerr(mS); qscore=%inf; essais=4;
    printf("\n-----     starting with score= %8.5f\n", score);
    while (score<qscore) | (essais<3) do
      if score==qscore then essais=essais+1; end//if
      qscore=score
      for j=1:jx do
        for tt=1:Omega do
        me=[mA(1:j-1,:);mT(tt,:);mA(j+1:$,:)]; me=me'*me;
        if det(me)<0.5 then bluk(j,tt)=%inf; continue; end//if
        bluk(j,tt)=max(mT*inv(me)*mT')
        end//do
      end//do
      lieu=find(matrix(bluk, 1,-1)<min(bluk)*1.0001)
      lieu=lieu(grand(1,'uin',1,size(lieu,'*')))
      xxx=pmodulo(lieu-1,jx)+1
      yyy=(lieu-xxx)/jx+1; score=bluk(xxx,yyy);
      printf("Xchange %4d with %4d, score= %8.5f\n",xxx,yyy,score);
      mA=[mA(1:xxx-1,:);mT(yyy,:);mA(xxx+1:$,:)]; mS=mA'*mA;
    end//while
  
    msm=min(spec(mS));
    if score < tscore  | ((score==tscore)&(msm>tmsm))  then 
      tA=mA; tS=mS; tscore=score; tmsm=msm;
      printf("%3d---------------------------- %f -- %f\n",t,tscore,tmsm);
      else printf("%3d-------- %f\n",t,score);
    end//if

    
  end//do
  mA=tA; mS=tS;
  Datas=[datas(1,1:$-1);itm2txt(cod2itm(mA))]; 
  printf(">>>-------------------randerr(mS)= %f\n",randerr(mS));
  printf(">>>-------------------randerr(ms)= %f\n",randerr(ms));
  printf("\n Datas, mA, mS generated in %d essais, laps= %6.4f sec., min(spec(mS))= %6.4f\n", tx, toc(), min(spec(mS)) )
mT=return(mT)
endfunction



menu()