#include utils.gau;

proc grids_fcts(h1,h2,startt,maxt,thmin1,thmax1,dth1,thmin2,thmax2,dth2);
local res, fct, a,b, adr,ngr1,ngr2,ptrue,probs,pp1,pp2, theta1,theta2;

probs=calprobs1(h1,startt,maxt);
pp1=probs[.,1]|probs[.,2];
pp1=pp1|(1-sumc(pp1));

probs=calprobs1(h2,startt,maxt);
pp2=probs[.,1]|probs[.,2];
pp2=pp2|(1-sumc(pp2));

ptrue=pp1|pp2;

ngr1=(thmax1-thmin1)/dth1+2;
ngr1=floor(ngr1);
theta1=seqa(thmin1,dth1,ngr1);

ngr2=(thmax2-thmin2)/dth2+2;
ngr2=floor(ngr2);
theta2=seqa(thmin2,dth2,ngr2);

res=1000+zeros(1,3);

for ia (1,rows(theta1),1);
a=theta1[ia];
for ib (1,rows(theta2),1);
b=theta2[ib];

adr=hsec;
fct= func_2mult(a,b,maxt,ptrue,2);
/* print /flush fct~a~b~(hsec-adr); */
res=res|(fct~a~b);

endfor;
endfor;

retp(res[2:rows(res),.]);
endp;


proc lp_max_value(aa,bb,nmet);

/*  This routibe returns the function value for the LP problem that 
    determines whether aa*x=bb has a solution with x(i) >= 0 for all i */

local fct,arows,acols,xm,xn,xmp,xnp,xm1,xm2,xm3,xica,eps,cc,startv,ii;


if (nmet.eq 1);

    arows=rows(aa); acols=cols(aa);
    aa=eye(arows)~aa;
    /* here we add the v's */
    
    cc=-ones(arows,1)|zeros(acols,1);
    startv=bb|zeros(acols,1)|zeros(arows,1);
    #ifos2win

           
        struct LP lp0;
        lp0=createLP;
        lp0.a=aa;
        lp0.b=bb;
        lp0.c=cc;
        lp0.l=0;
        lp0.u=1e200;
        
        struct LPcontrol c0;
        c0=createLPcontrol;
        c0.output=0;
        c0.constraintType=3;
        c0.minimize=1;
        c0.start=startv;
        
        struct LPout out0;
        out0=lpmt(lp0,c0);
        fct=out0.optimumvalue;
        print "return code   ";; out0.returncode; 
        print "function value";; fct; 
   
        
    #else
       print /flush "OS not recognized"; 
       stop;
    #endif
elseif (nmet.eq 2);
    aa=bb~(-aa);
    
    arows=rows(aa); acols=cols(aa);
    aa=aa~-eye(arows);
    xm=rows(aa); xn=cols(aa)-1;
    
    aa=(0~zeros(1,acols-1)~(-ones(1,arows)))|aa;
    aa=aa|zeros(1,cols(aa));
    
    xmp=rows(aa); xnp=cols(aa); xm1=0; xm2=0; xm3=rows(bb);
    xica=0; 
    eps=1.0e-8;
    aa=aa';
    
    #IFUNIX
       dlibrary /home/honore/adriana/simplex1.so;
       dllcall simplex_(aa,xm,xn,xmp,xnp,xm1,xm2,xm3,xica);  
    #else
    #ifos2win
       dlibrary c:\research\forgau95\sim_new\release\sim_new.dll;
       dllcall isimpl(aa,xm,xn,xmp,xnp,xm1,xm2,xm3,xica,eps); 
    #else
       print /flush "OS not recognized"; 
       stop;
    #endif
    #endif
    if (xica.le -100);
        print "icase=";;
        print xica;;
        print " in lp_max_val";
        stop;
    endif; 
    fct=aa[1,1];
else;
    print /flush "nmet not valid in lp_max_val";
endif;
retp(fct);
endp;