// Pierre L. Douillet
// www.douillet.info
// version 21
// module planx

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

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

// 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(['lire','traiter','aleas','better'],..
        ['PLANX.20: menu principal'])
    select item
    case 0 then return
    case 1 then readdata(); 
    case 2 then menu_traiter();
    case 3 then sig=x_mdialog('recherche aléatoire',['nb d''essais'],['30']);
      if sig ~= [] then goodplan(eval(sig(1))); end//if
    case 4 then sig=x_mdialog('recherche dirigée',  ['nb de blocs'],['15']);
      if sig ~= [] then better(eval(sig(1))); end//if
    end//select
  end//while
endfunction

function menu_traiter
global fun maa
  while %t do
  item=x_choose(['mambms','draw','map','plages', 'xy+yz', 'xx+xy', 'lenteur', 'plot3d'],..
      ['traiter: choisir un item'])
  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
    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');
  case 5 then
    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
    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)
  case 7 then
    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)
  case 8 then sig=x_mdialog('response surface',['grid','s'],['30','1']);
    drawfig(eval(sig(1)), eval(sig(2)) )
  end//select
  end//while
endfunction

//*********************************************************
//* lecture des données
//*********************************************************

function rep=readdata()
global shortnam datas 

item=x_choose(['xper','usxper','caout','clair','eddy','nine','eval07b','baie01','baie02', 'baie03'],..
    ['choisir le fichier à lire'])

select item
case 0 then return
case 1 then 
  nam = "experimentique.txt"; 	shortnam="xper";
case 2 then 
  nam = "us-experimentique.txt"; 	shortnam="usxper";
case 3 then 
  nam="caoutchouc.txt"; 	shortnam="caout"
case 4 then 
  nam="claire.txt"; 		shortnam="clair"
case 5 then 
  nam="nist-eddy.was.txt"; 	shortnam="eddy"
case 6 then 
  nam="nine3.txt"; 	shortnam="nine"
case 7 then 
  nam="eval07b.txt"; 	shortnam="eval07b"
case 8 then 
  nam="baie01.txt"; 	shortnam="baie"
case 9 then 
  nam="baie02.txt"; 	shortnam="baie"
case 10 then 
  nam="baie03.csv"; 	shortnam="baie"
end//select

// 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 shortnam
case "xper" then 
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, '°C', '') ;
case "usxper" then 
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, '°C', '') ;
case "clair" then
  rep = strsubst( rep, ' ', '_') ; rep = strsubst( rep, ',', '.') ; rep = strsubst( rep, ';', ' ') ;
case "eval07b" 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');
end//select

// format et test de longueur du format
toks=tokenpos(rep(1)) ; toksx=size(toks,'r') ; // [" ", ascii(9)]
frmt='%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s' ; 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


select shortnam
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]) 
end//select

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

readmoda

endfunction

function readmoda
global ix 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 !!!

// liste des modalités
moda=list() ; siza=1:ix ; zona=1:ix+1 ;
for i =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)
  itm=[1:ix];
  for i =1:ix do
  ici=vectorfind(moda(i), txt(i), 'c'); 
  if ici == [] then error (sprintf("==> : undefined code %s for_modality %s", txt(i), datas(1,i) )); end//if
  itm(i)= ici-1 ; end//for
endfunction

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

function cod=itm2cod(itm)
  cod=[1] ; for i=1:ix do
  ici=itm(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
endfunction

function num = itm2num(itm)
  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
  itm=[] ; 
  for i =1:ix-1 do [num, tmp]= euclid(num, 1+siza(i)); itm=[itm, tmp]; end//for
  itm=[itm,num]
endfunction

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

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

function mat=data2mat(data)
  mat=[]
  for j =1:jx do mat=[mat;itm2cod(txt2itm(data(j+1,:)))]; end//for  //; special
endfunction

// calculs de base
function frv=mambms ()
global ma mb ms mx mdelta
  ma=data2mat(datas) ; mb= eval(datas(2:$, $)) ;
  ms = ma' * ma ; mx= (1/ms) * ma' * mb ; mdelta= ma*mx-mb ;
  redvar= variance(mdelta); rawvar= variance(mb);
  if lx==jx 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 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 [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
  command='[x,lagr,f]=quapro(-qqq,-ppp,[],[],-1.2*ones(ix,1),+1.2*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 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

function rr=recmap(tox,toy,DAT)
// DAT est soit une matrice of strings (=datas)
//         soit une matrice d'items
  txy=[tox,toy]
  function num = itm2numxy(itm)//inside
    num=0
    for i=txy do num = num*(1+siza(i))+itm(i) ; end//for
  endfunction//inside
  function nums=plan2numxy(plan, txt)//inside
    nums=ones(1,jx)
    if txt then
      for j =1:jx do nums(j)= itm2numxy(txt2itm(plan(j+1,:))); end//for
    else
      for j =1:jx do nums(j)= itm2numxy(plan(j,:)); end//for
    end//if
  endfunction//inside

  rr=zeros(prod(siza(toy)+1), prod(siza(tox)+1)) ;
  ptn=plan2numxy(DAT, type(DAT(1))==10); 
  for j=ptn do rr(j+1)=rr(j+1)+1; end//for
endfunction

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

function [Datas,mA,mS]=randplan()
  Datas=datas(1,1:ix);
  if jx<sqrt(Omega) then 
    truc=grand(1,jx, 'uin', 0, Omega-1); //ligne pour boucle for
    for j=truc do Datas=[Datas;itm2txt(num2itm(j))]; end//for //; special
    mA=data2mat(Datas); mS=mA'*mA;
  else
    lesj=0:Omega-1; for k=Omega:-1:Omega-jx+1 do
      loc=grand(1,1,'uin',1,k); j=lesj(loc); lesj=[lesj(1:loc-1),lesj(loc+1:$)];
      Datas=[Datas;itm2txt(num2itm(j))]; 
    end//for
    mA=data2mat(Datas); mS=mA'*mA;
  end//if
if det(mS)<0.00001 then [Datas,mA,mS]=randplan(); end//if
endfunction


function goodplan(tx)
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("Datas, mA, mS generated in %d essais, laps= %6.4f sec., min(spec(mS))= %6.4f\n", tx, toc(), memscore)
endfunction

//------------------------------------le package better-----------------------------
// 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=[];
  for j =1:Omega do mT=[mT;itm2cod(dT(j,:))]; end//for
endfunction

function [err,qui]=maxerr(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()
//scoped mT
  ims=(1/mS)
  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; [score, nouv]= randerr(); score0=score; mAug=[mA;mT(nouv,:)];
  mtry=mAug(2:$,:); scotry=maxerr(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=maxerr(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 qo= better(tx)
global mA mS Datas
//scoping mT
  tic(); clear('mT'); mT=tous(); qo=%inf;
  for q=1:tx do
    [Datas,mA,mS]=randplan(); score=stepbetter();
    if score<qo then 
      printf("%3d----------------------------- %f\n",q,score); qo=score; qa=mA; qs=mS;
    else
      printf("%3d-- %f\n",q,score);
    end//if
  end//for
  mA=qa; mS=qs; Datas=[datas(1,1:$-1);itm2txt(cod2itm(mA))] 
  printf(">>>-------------------maxerr(mS)= %f\n",qo);
  printf(">>>-------------------maxerr(ms)= %f\n",maxerr(ms));
  printf("\nDatas, 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 drawfig(qq, mys)
top=1.3
 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);
 clf(); curfig=gcf(); curfig.color_map=jetcolormap(32); 
 surf(yy,xx,zz,'facecol','interp'); 
 curax=gca(); curax.rotation_angles=[0,-90];  
endfunction
//-----------------------------------------------------------------


menu

function rapport
global fun maa
  readdata()
  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