//
// Take a representative set for domain parameter mu
// Want to find choice of mu that gives largest error estimator

real errorestimator0=0;

real bestk0;

real ds=1;// particion fina
//
cout<<"Start Search... "<< n <<endl;
for(int j=1;j<=(k0max-k0min);j++)
{ 

//
// TO MAKE SURE DO NOT REPITE PARAMETER AND INDEX
// I AM SURE IT CAN BE MADE BETTER!
//
//cout<<" Parameter index j = "<< j <<" Best parameter index jbest "<< jbest[cont]<<endl;
//cout<< " Check  (j!=jbest[cont]) "<<(j!=jbest[cont])<<endl;

// 
//IDEA
//     WHEN sum=n THE INDEX j HAS NOT BEEN USED BEFORE
//
	real sum=0;
	   //for(int s=0;s<=cont;s++)
	for(int s=0;s<n;s++)
	 {//cout<<" j=  "<<j <<" jbest[s]=  "<<jbest[s]<< " bool "<<(j!=jbest[s])<<endl;
	 sum=sum+(j!=jbest[s]); //sum=cont ---> j!=jbest[s] cierto (igual a 1) para todos
	 // sum<cont---> j!=jbest[s] falso (igual a 0) para alguno
	 //cout<<" Valor de sum "<<sum <<endl;
	      } 
	  if (sum==n) // NOW WE USE INDEX j
	  {
	  k0=k0min+ds*j;  // PARAMETER ACCORDING TO j

//
// FIX ALL THE PARAMETERS ACCORDING TO THIS ONE
//
       include "parameters.idp";  
//
//SOLVE THE PROBLEM ACCORDING TO THIS PARAMETER WITH REDUCED BASIS
//
include "ReduceBasisSolver.idp"
//
// SEE SOLUTION
//
// isovalues
//
for (int i=0;i<viso.n;i++)
viso[i]=RBu[].min+i*(RBu[].max-RBu[].min)/mm;

plot(cmm= " Search "+n+ ": RB Solution for k0 = "+k0,RBu,viso=viso(0:viso.n-1),hsv=colorhsv,value=1);

//
// COMPUTE ERROR ESTIMATE:
// RIESZ REPRESENTATION OF L(*)-a(RBu,*) en Vh. 
//
// WE HAVE ALREADY  Riesz Representation of L(*) 
//
// AND Riesz Representation of a(RBu,*) is given by 
//  
//    (Riesz_a(RBu,*), v)_V=a(RBu,v), for all v en Vh
//
// COMPUTATION OF  a(RBu,*)
//
varf aa0(unused,w)=int3d(Th3,0,2,4,6,8,10)(
                             dx(RBu)*dx(w)+dy(RBu)*dy(w)+dz(RBu)*dz(w)    
                            );
real[int] aa0vec=aa0(0,Vh);

varf aa1(unused,w)=int3d(Th3,1,5,9)(
                             dx(RBu)*dx(w)+dy(RBu)*dy(w)+dz(RBu)*dz(w)    
                            );
real[int] aa1vec=aa1(0,Vh);

varf aa2(unused,w)=int3d(Th3,3,7,11)(
                             dx(RBu)*dx(w)+dy(RBu)*dy(w)+dz(RBu)*dz(w)    
                            );
real[int] aa2vec=aa2(0,Vh);

//
varf aa3(unused,w)=int2d(Th3,10,12,14,16,18,110,11,15,19,13,17,111,30,32
                   ,34,36,38,310,31,35,39,33,37,311,40,44,48,41,45,49
                   ,20,24,28,23,27,211,199,201,299,301,399,401,499,501
                   ,599,601,699,700,701)(RBu*w);

real[int] aa3vec=aa3(0,Vh);
//
varf aa4(unused,w)=int3d(Th3)(RBu*w);
real[int] aa4vec=aa4(0,Vh);
//
// VECTOR ASSOCIATED TO a(RBu,*) ACTING ON FEM BASIS
//
Vh rhsforRieszaRBu;
rhsforRieszaRBu[]=k0*aa0vec;
rhsforRieszaRBu[]=rhsforRieszaRBu[]+kleft*aa1vec;
rhsforRieszaRBu[]=rhsforRieszaRBu[]+kright*aa2vec;
rhsforRieszaRBu[]=rhsforRieszaRBu[]+aa*aa4vec;
rhsforRieszaRBu[]=rhsforRieszaRBu[]+Bi*aa3vec;
// 
// RIESZ REPRESENTATION COMES FROM INVERTING THE MATRIX
// OF THE SCALAR PRODUCT 
//
Vh RieszaRBu;
RieszaRBu[]=Ascalarproduct^-1*rhsforRieszaRBu[];
//
// VECTOR THAT HAS RIESZ REPRESENTATION OF RESIDUAL
//
Vh diff=RieszaRBu-Rieszf;
//
// ERROR ESTIMATOR IS THE NORM OF DIFF 
//
real part1= int3d(Th3)( dx(diff)^2+dy(diff)^2+dz(diff)^2); 
real part2=int3d(Th3)(diff^2);                   
real normdiff=sqrt(part1+part2*invvol); 
real errorestimator= normdiff/k0min;                 

cout <<"Search "<<n<<" and iter j = "<<j<<". For k0 = "<<k0 <<": error estimator "<<errorestimator<<endl;
//
// SEARCH FOR THE LARGEST ERROR ESTIMATOR
//
if(errorestimator>errorestimator0)
    {   
	errorestimator0=errorestimator;
	bestk0=k0;
	jbest[n]=j;
    }
} 
}
cout <<"Search "<<n<<". Best k0 = "<<bestk0 <<": error estimator "<<errorestimator0<<endl;