//Algorithme du simplexe
//www.douillet.info
// version 06

myfuncprot=funcprot(); funcprot(0);
cd "~/docs/Ensait/oprea/scilab";


//next comming are gone 
//global lesms  //la suite des matrices ms

global myseed cramerflag
global vue fvu ms fff bor oldpiv newpiv


function simplex_one()
  global cramerflag
  global vue fvu ms fff bor oldpiv newpiv
  cramerflag=[]
  myieee=ieee(); ieee(2); 
  mymprintf("seed= %d\n", myseed);
  fvu="vue= "; for j=1:size(vue,'c') do; fvu=fvu+"%d "; end//do
  fvu=fvu+"\n";
  mymprintf(fvu, vue);
  ttt([1:size(ms,'c')-1])="%2.0f"; ttt(vue)="%+10f"; fff=msprintf("%s ",ttt)+"%+13f\n";
  mymprintf(fff, ms)
  id=eye(ms);

  while %t do
    myprintf("------------------------------\n");
    qq=find(ms(1:$-1,$)<=sqrt(%eps)); cramerflag=sort([qq,vue]);
    dd=find(cramerflag(1:$-1)==cramerflag(2:$)); cramerflag([dd,dd+1])=[];
    if cramerflag<>[] then
      mymprintf("\n--- not a Cramer system -- check eqn : %d\n\n", cramerflag');
      if ~toscr then 
        mprintf("\n--- not a Cramer system -- check eqn : %d\n\n", cramerflag');
      end//if
      break
    end//if
    
    coe=ms($,1:$-1); piv=find(coe>0);
    if size(piv,'*')==0 then 
      break; 
    end//fi
    
    oldpiv=piv(grand(1,'uin',1,size(piv,'*')));

    ccd=ms(1:$-1,oldpiv); ccd(vue)=ones(vue');
    bor1=-ms(1:$-1,$) ./ ccd; bor=bor1;
    tmp=find(bor<=0); 
    bor(tmp)=ones(tmp')*%inf;
    if %f then
      tmp=[ms(1:$-1,oldpiv),ccd,bor1,bor];
      tmr=msprintf("%+12f %+12f %+12f %+012f\n", tmp);
      tmr=strsubst(tmr,"Inf","+Inf--------");
      mymprintf("%s\n",tmr)
    end//if print décision
    [qqq,newpiv]=min(bor);

    pp=zeros(ms); pp(vue,:)=id(vue,:);
    ligne=ms(newpiv,:);
    fac=ligne(oldpiv); ligne(oldpiv)=0; ligne(newpiv)=-1;
    pp(oldpiv,:)=-ligne/fac;
    pp($,:)=id($,:);

    tmq=find(vue==oldpiv); vue(tmq)=newpiv; vue=gsort(vue,'g','i');
    mymprintf(fvu, vue);
    ms=ms*pp;
    ttt([1:size(ms,'c')-1])="%2.0f"; ttt(vue)="%+10f"; fff=msprintf("%s ",ttt)+"%+13f\n";
    mymprintf(fff, ms)
    
  end//while %t
  ieee(myieee)
endfunction//simplex_one 



function ms=simplexe(ma,mb,mza,mzb,seed, toscr,tofil,filna)
  global myseed cramerflag
  global vue fvu ms fff bor oldpiv newpiv

  printf("\n");

  if isdef('seed') then 
    myseed=seed
  elseif myseed <> [] then //une variable globale est toujours créée. Elle vaut alors [].
    myseed=myseed+1
  else
    myseed=12345
  end//if
  grand('setsd', myseed);
    
  if ~isdef('toscr') then toscr=%f; end//if

  if ~isdef('tofil') then 
    tofil=%f;
  elseif tofil==%t
    if ~isdef('filna') then filna="msg"; end//if
  end//if

  if isdef('filna') then 
    tofil=%t;
    if length(filna)<4 | part(filna,[length(filna)-3:length(filna)]) <> ".txt" then
      filna=filna+".txt";
    end//if
  else
    filna="none";
  end//if
  
//printf("6 seed=%s, toscr=%s, tofil=%s, filna=%s\n", string(isdef('seed')), string(isdef('toscr')), string(isdef('tofil')), string(isdef('filna')) )
//printf("myseed=%d, toscr=%s, tofil=%s, filna=%s\n", myseed, string(toscr), string(tofil), filna)
  

  if tofil then hd=mopen(filna,'w'); end//if
  tmpa=[""]; tmpb=[""];
  if toscr then tmpa=[tmpa;"mprintf(aa)"];     tmpb=[tmpb;"mprintf(aa,bb)"]; end//if
  if tofil then tmpa=[tmpa;"mfprintf(hd,aa)"]; tmpb=[tmpb;"mfprintf(hd,aa,bb)"]; end//if
  deff("myprintf(aa)", tmpa); deff("mymprintf(aa,bb)", tmpb);

  // on commence à travailler

  [neq,nxx]=size(ma); vue=1:nxx;

  ms=[eye(nxx,nxx),zeros(nxx,neq+1)];
  ms=[ms;[-ma;mza],zeros(neq+1,neq), [mb;mzb]];

  spc=find(ms(1:$-1,$)<0); nsp=size(spc,'*'); 

  if nsp >0 then
    spid=eye(neq+nxx+nsp+1, neq+nxx+1);
    finms=ms($,:); ms=ms(1:$-1,:);
    for vv=spc do
      ms=[ms;-ms(vv,:)]; ms(vv,:)=spid(vv,:);
      ms($,vv)=1; vue=[vue, vv];
    end//for
    ms=[ms(:,1:$-1),zeros(neq+nxx+nsp,nsp),ms(:,$)];
    ms=[ms; -sum(ms(neq+nxx+1:$,:),'r')];
  end//if

  simplex_one()

  if nsp>0 & cramerflag == [] then
    if vue(nxx+1:$)<> [nxx+neq+1:nxx+neq+nsp] then
      myprintf("\n--- overdeterminated, not solvable system ---\n\n");
      if ~toscr then 
        mprintf("\n--- overdeterminated, not solvable system ---\n\n");
      end//if
    else
      ms($-nsp:$,:)=[]; ms(:,nxx+neq+1:$-1)=[];
      ms=[ms; [mza,zeros(1,neq)]*ms]; ms($,$)=ms($,$)+mzb;
      vue=vue(1:nxx);
      simplex_one()
    end//if vue
  end//if deux phases

  myprintf("--done------------------------\n");
  tmp=sprintf("myseed=%d, toscr=%s, tofil=%s, filna=%s\n", myseed, string(toscr), string(tofil), filna);
  //mymprintf("%s", tmp)
  if tofil then hd=mclose(hd); end//if  
    
endfunction//simplexe


//select 4
//case 1
  ma=[1 0; 0 1; 3 6];
  mb=[1000 500 4200]';
  mza=[4 12] ;
  mzb=[0];
simplexe(ma,mb,mza,mzb, seed=18, toscr=%t, filna="msg1")
//case 2
  ma=[1 0 0; 0 1 0; 0 0 1; 3 6 2];
  mb=[1000 5000 1500 6750]';
  mza=[4 12 3]; 
  mzb=[0];
simplexe(ma,mb,mza,mzb, toscr=%t, filna="msg2")
//case 3
  ma=[1 0; 0 1; 3 6; -1 -1; -1 -2];
  mb=[1000 500 4200 -1000 -1200]';
  mza=[4 12] ;
  mzb=[0];
simplexe(ma,mb,mza,mzb, seed=24, toscr=%t, filna="msg3")
//case 4
  ma=[1 0; 0 1; 3 6; -1 -1; -1 -2];
  mb=[1000 500 4200 -1000+0.00001 -1200]';
  mza=[4 12] ;
  mzb=[0];
simplexe(ma,mb,mza,mzb, toscr=%t, filna="msg4")
//case 5
  ma=[1 0; 3 -6; -1 1; -1 2];
  mb=[1000 1200  -500  -200]';
  mza=[4 -12] ;
  mzb=[6000];
simplexe(ma,mb,mza,mzb, toscr=%t, filna="msg5")
//case 6
  ma=[1 0; 0 1; 3 6; -1 -1];
  mb=[1000 500 4200 -1500]';
  mza=[4 12] ;
  mzb=[0];
simplexe(ma,mb,mza,mzb, seed=18, toscr=%t, filna="m6")

//end//select