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

ident='planx.30';
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 mT
global ma maa mb mc ms ims mx mdelta Datas mA mS imS mX mDelta fun idents describe

global qx ; qx= 50;  //qx used by XchgOUT1
verb=%T;

// A Scilab function can read all of the variables of the calling function (scoped)
// When the function tries to write into a variable, this variable becomes 'local'
// To enable a function to write into a variable of the caller,
// either BOTH must have declared this variable as global
// or the clear/return (aka scoping) trick must be used

//*********************************************************
//* 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 iterrandplan(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"
  "deux_deux.txt"         "22"         "none"
  "deux_deux_deux.txt"    "222"        "none"
  "trois_trois.txt"       "33"         "none"
  "trois_trois_trois.txt" "333"        "none"
  "trois_quatre.txt"      "34"         "none"
  "trois_cinq.txt"        "35"         "none"
  "fractional.txt"        "frac"       "none" 
  ];
describe_old=[
  "eval07b.txt"           "eval07b"    "eval07"  
  "eval07c.txt"           "eval07c"    "eval07"
  "baie01.txt"            "baie1"      "baie"
  "baie02.txt"            "baie2"      "baie"
  "waeghe.csv"            "waeghe"     "none"
]

describe=describe(:,[2,1,3]);
describe=gsort(describe, 'lr', 'i');describe=describe(:,[2,1,3]);
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);

// reading each line of the file (as a single string)
fich = mydir + nam ; printf('source: %s\n',fich)
hndl = mopen(fich, "r") ; rep = mgetl(hndl) ; mclose(hndl) ;

// preformatting the data
// ... cela change selon les fichiers !!!!
// caveat : this changes when file changes !!!!
select french
case "xperfr" then 
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, code2str(276)+"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


// testing the length of the 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;

// cutting each line into words (giving a rectangular matrix... or an error)
datas = msscanf(-1, rep, frmt) ;
if size(datas,'r')==1 then disp(rep(1:3)); error 'something happens' ; end//if

ifx=0; // ifx is the number of continuous parameters (left columns)
select french
case "eddy" then
// traitement special : les resultats sont en colonne 1, le # en derniere
  tmp=['0';datas(2:$,$)]; [qqq,trad]=gsort(eval(tmp),'r','i') //row, increasing
  datas=datas(trad,[2,3,4,1]) 
case "baie" then
// traitement special : 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()
tous();

endfunction

function readmoda
global ix ifx jx lx moda siza zona Omega
  // first row contains titles
  // last column contains results
  [jx,ix]= size(datas); ix=ix-1; jx=jx-1 // decremente ix et jx !!!
  // ifx set in menu_read, according to french
  
  // listing the modalities of each factor
  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)))' ;
   if size(moda(i), '*') == 1 then printf("\nparameter %d takes only onve value \n")
			error "stopped"; 
		end//if
   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)// size of the product space

endfunction

//-------------- gestion des modalites, codage/decodage -----------------------
// txt = ligne de texte
// itm = ligne de numeros, de 0 à  siza(i)
// num = numero unique, de 0 à  Omega-1
// cod = le codage en variables indepeNdantes 

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 is a column vector
  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 is either a string matrix (=datas)
//        or     an item  matrix
  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 ims 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 ; ims=inv(ms); mx= ims * ma' * mb ; mdelta= ma*mx-mb ;
  redvar= variance(mdelta); rawvar= variance(mb);
  if (lx==jx) | (abs(redvar)<%eps) then 
    frv= %nan; fqual=%nan;
  else 
    frv= rawvar/redvar*(jx-lx)/(jx-1); fqual= (rawvar-redvar)/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 imS mX mDelta
  mA=maa ; mb= eval(datas(2:$, $)) ;
  mS = mA' * mA ; imS=inv(mS); mX= imS * 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-etoupe');
  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))';
  command='[x,lagr,f]=qld(-qqq,-ppp,[],[],-top*ones(ix,1),+top*ones(ix,1),0)';
  printf('%s\n\n', command);
//  [x,lagr,info]=qld(Q,p,C,b,ci,cs,me [,tol])
  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

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

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


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

// global mA mS mT   (deja declares)   

function tous()
global mT
  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,ran,qui]=randerr(mSS)
// err value, RAND err loc, all err loc
  if det(mSS)<0.5 then [err,ran,qui]=(%inf,%nan,%nan); return ; end//if
  ims=inv(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 [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 jscore=XchgIN1()
global mA mS
  qscore=%inf; jscore=randerr(mS) ; essais=1;
  if verb then printf("\n-----     starting with score= %8.5f\n", jscore); end//if 
  while (jscore<qscore) | (0<essais) do
    if jscore==qscore then essais=essais-1; end//if
    qscore=jscore
    [jscore, k_in]= randerr(mS); jloc=0;
    score0=jscore;
    for j=1:jx do
      meA=[mA(1:j-1,:);mT(k_in,:);mA(j+1:$,:)]; meS=meA'*meA;
      score=randerr(meS);
      if score<jscore then jscore=score; jmA=meA; jloc=j; end//if
    end//for
    if jscore<score0 then mA=jmA; mS=mA'*mA; 
      if verb then printf("Xchange %4d with %4d, score= %8.5f\n", jloc, k_in, jscore); end//if
    end//if  
  end//while  
endfunction

qx=50;

function [qscore, qloc]=XchgOUT1()
global mA mS Datas
  qscore=randerr(mS); qloc=0; // mS est inversible
  if verb then printf("\n-----     starting with score= %8.5f\n", qscore); end//if
  for q=1:qx 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
    if verb then printf("%d %d ======= nan encountered ; nouv=%d\n", t,nouv); end//if   
    end//if      

    qmA=[mA(1:kill-1,:);mT(nouv,:);mA(kill+1:$,:)]; qmS=qmA'*qmA; score=randerr(qmS);
    if score<qscore then
      mA=qmA; mS=qmS; qscore=score;
      if verb then printf("Xchange %4d with %4d, score= %8.5f\n", kill, nouv, score); end//if
      //printf("%3d-- %f\n",t,score);
      // vieux else printf("%3d-- %f\n",t,scotry);
    end//if
  end//do
endfunction


function score=XchgBOTH1();
  global mA mS
// part d'un randplan donne
// amelioration par double echange
// continue tant que cela ameliore
  bluk=zeros(jx,Omega);
  score=randerr(mS); qscore=%inf; essais=1;
  if verb then printf("\n-----     starting with score= %8.5f\n", score); end//if
  while (score<qscore) | (0<essais) 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);
    if verb then printf("Xchange %4d with %4d, score= %8.5f\n",xxx,yyy,score); end//if
    mA=[mA(1:xxx-1,:);mT(yyy,:);mA(xxx+1:$,:)]; mS=mA'*mA;
  end//while
endfunction


deff('score=fibre_one()','score=randerr(mS)');

function [qscore, loc]=old_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 [tscore, tloc]=monte_carlo(fibre, tx)
global mA mS Datas 
  tic();  tscore=%inf; tmsm=-%inf;
  for t=1:tx do
    [Datas,mA,mS]=randplan();
    score=fibre(); msm=min(spec(mS));
    if score < tscore  | ((score==tscore)&(msm>tmsm))  then 
      printf("%3d---------------------------- %8.5f -- %8.6f\n",t,score,msm);
      tscore=score; tmA=mA; tmsm=msm; tloc=t;
      else printf("%3d-------- %8.5f\n",t,score);
    end//if
  end//for
  mA=tmA; mS=mA'*mA; Datas=[datas(1,1:$-1);itm2txt(cod2itm(mA))]; 
  printf(">>>-------------------randerr(mS)= %f, min(spec(mS))= %f\n", tscore,      tmsm);
  printf(">>>-------------------randerr(ms)= %f, min(spec(ms))= %f\n", randerr(ms), min(spec(ms)));
  printf("\n Datas, mA, mS generated in %d rand_blocs, laps= %6.4f sec.\n", tx, toc())
endfunction

function iterrandplan(tx)
  monte_carlo(fibre_one, tx)
endfunction

// monte_carlo(fibre_one, 10)
// monte_carlo(XchgBOTH1, 10)
// monte_carlo(XchgOUT1, 10)
// monte_carlo(XchgIN1, 10)

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

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//do

 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//do

 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






menu()

// menu_read('nine3.txt'); mambms(); jx=20; [rDatas,rmA,rmS]=randplan();

if %f then
  for ii=1:10 do
  [rDatas,rmA,rmS]=randplan();
  mA=rmA; mS=rmS; err1=XchgBOTH1();
  mA=rmA; mS=rmS; err2=XchgIN1();
  printf("start=%9.6f IN=%9.6f BOTH=%9.6f\n",randerr(rmS),err2,err1);
  if err2<err1 then break; end//if
  end//do
end//if