function sampler(NO,NI,NT,AT,LS,PT,EC)

%clear all;  
path(path,'../COMMON/');
warning off;

%this script is generally pretty stable, and should be used
%for batch processes

numsteps     =  NO^2 ;
innersteps   =  NI   ;
num_trials   =  NT   ;
accept_type  =  AT   ;  % 1)  for (p3*t31) / (p1*t13)
                        % 2)  for (p3+p2)  / (p1+p4)               
landscape_type  = LS ;
pert_type       = PT ;
pe_flag         = EC ;

B = 45 ;
x_lower    = -B;          x_upper = B;
y_lower    = -B;          y_upper = B;
Volume     = (2*B)^2;
%
num_paths=1;
switch pert_type
    case 1
% bias balanced orthogonal  
      dx_12   = 0.5 ;   dx_max  = 0.0 ;
      dy_12   = 0.0;    dy_max  = 0.2; 
    inner_parms_dx = [x_lower x_upper y_lower y_upper...
    innersteps dx_max dy_max dx_12 dy_12 num_trials num_paths];
      dx_12   = 0.0 ;   dx_max  = 0.2 ;
      dy_12   = 0.5;    dy_max  = 0.0; 
    inner_parms_dy = [x_lower x_upper y_lower y_upper...
    innersteps dx_max dy_max dx_12 dy_12 num_trials num_paths];
%
    case 2
% bias balanced hybrid        
      dx_12   = 0.5 ;    dx_max  = 0.2 ;
      dy_12   = 0.0 ;    dy_max  = 0.2; 
    inner_parms_dx = [x_lower x_upper y_lower y_upper...
    innersteps dx_max dy_max dx_12 dy_12 num_trials num_paths]; 
      dx_12   = 0.0 ;   dx_max  = 0.2 ;
      dy_12   = 0.5;    dy_max  = 0.2; 
    inner_parms_dy = [x_lower x_upper y_lower y_upper...
    innersteps dx_max dy_max dx_12 dy_12 num_trials num_paths];
%
    case 3
%        
% diagonal (already balanced)
      dx_12   = 0.5 ;   dx_max  = 0.2 ;
      dy_12   = 0.5;    dy_max  = 0.2; 
      inner_parms_dx = [x_lower x_upper y_lower y_upper...
      innersteps dx_max dy_max dx_12 dy_12 num_trials num_paths];
      inner_parms_dy=inner_parms_dx;
%
end
%

fname=[ 'NO_',num2str(NO),...
       '-NI_',num2str(NI),...
       '-NT_',num2str(NT),...
       '-AT_',num2str(AT),...
       '-LS_',num2str(LS),...
       '-PT_',num2str(PT),...
       '-EC_',num2str(EC)];
%     
run landscapes; G_fac=1;
%
switch landscape_type
    case 1
        %single basin
        parms = landscape.SB;  G_fac=0.2;  
    case 2
        %orthogonal disjoint 
        parms = landscape.OB;
    case 3
        %orthogonal overlapping
        parms = landscape.OB; G_fac=0.2;
    case 4
        %long axes at acute angles, accessible
        parms  = landscape.AB; 
    case 5
        %broad, unimodal, asymmetric potential
        parms  = landscape.AB2; 
    case 6
        %asymmetric egg crate
        parms = landscape.MB;  G_fac=3;
end
%

dx=sqrt( Volume/numsteps );
dx=2*dx;
x3D=x_lower : dx : x_upper-dx;
y3D = x3D;
parms(:,6)=parms(:,6)*G_fac;
ratio=1;failcount=0;
%this first script will test the following:

%**********************************************************

%visualization of potential:
e_log=[]; ew=1;
 for i=1:length(x3D)
  for j=1:length(y3D)
   f(i,j) = e_V( y3D(i)  , x3D(j)  ,parms);
    f_ij   = e_V( y3D(i)  , x3D(j)  ,parms);
    e_log  = [ e_log -log( f_ij ) ];
   %e_log  = [e_log -log(f(i,j) )];
  end
 end
f=f/sum(sum(f))/dx^2;
 
nbins=200;  
 de       = ( 8 -min(e_log) )/ (nbins-1);
 e_plt    = min(e_log):abs(de):8;
 h_e.true = histc(e_log,e_plt);
 e_log=[];
% %h_e.true = h_e.true.*exp(-(e_plt-max(e_log)) );
% %h_e.true = h_e.true/sum(h_e.true)/de;
% V3D=-log(f);  
% trip=innersteps

%multiscale settings()()()()()()()()()()()()()()()()()()
y_test=(y_upper-y_lower)*rand + y_lower;
x_test=(x_upper-x_lower)*rand + x_lower;
p1=e_V(x_test, y_test, parms);
q1=[x_test y_test]';
num_paths=5;
simple_parms = [x_lower x_upper y_lower y_upper 0.5];
inner_parms= [x_lower x_upper y_lower y_upper...
    innersteps dx_max dy_max dx_12 dy_12 num_trials num_paths];
xy_trials=[]; xy_accept=[]; N_acc=0;
h_e.sim=e_plt*0; e_log=[];
AR=0;
P1=1; 
x_plt=x_lower:dx:x_upper-dx;
y_plt=y_lower:dx:y_upper-dx;
x_pmf=x_plt*0; y_pmf=y_plt*0;
nbins_work=50;S=1;L=[0 0 0]';L_plt=[]; pe_arr=1;
ew_r=1; pe=1; e_Wd_r=[];
PP0=[];PP1=[];PP2=[];L1=0;L2=0;
accept=1;
pe_plt1=[]; pe_plt2=[]; pe1=1; pe2=2;
pp0=1; pp1=1; pp2=1;
xy_pmf=x_pmf*y_pmf';
h_x=x_pmf; h_y=y_pmf;
h_xy=xy_pmf;


pmf.y = Landau(y3D,parms,'x');
pmf.y = pmf.y/sum(pmf.y)/(y3D(2)-y3D(1));
pmf.x = Landau(x3D,parms,'y');
pmf.x = pmf.x/sum(pmf.x)/(x3D(2)-x3D(1));
%=============================================



for i=1:numsteps
     if (i==1)
        PP0(1)=1; PP1(1)=1; PP2(1)=1;
        pe_plt1(1)=1;pe_plt2(1)=1;
        L_plt(:,1)=L;
     else
       if (accept)
           %only the accepted 
           %forward trajectories are in equilibrium
           PP0=[PP0 pp0]; PP1=[PP1 pp1];  PP2=[PP2 pp2];
           L_plt=[L_plt L];
           pe_plt1=[pe_plt1 pe1];
           pe_plt2=[pe_plt2 pe2];
%        else
%           ls=length(PP0);
%            PP0=[PP0 PP0(ls)]; PP1=[PP1 PP1(ls)];  PP2=[PP2 PP2(ls)];
%            L_plt=[L_plt L_plt(:,ls) ];
%            pe_plt1=[pe_plt1 pe_plt1(ls) ]; 
%            pe_plt2=[pe_plt2 pe_plt2(ls) ];
       end
       maxwin=5000;
       if (length(PP0)==maxwin+1)
         tmp = PP0(2:maxwin+1);  PP0=tmp; 
         tmp = PP1(2:maxwin+1); PP1=tmp;
         tmp = PP2(2:maxwin+1); PP2=tmp;
         tmp = pe_plt1(2:maxwin+1); pe_plt1 = tmp;
         tmp = pe_plt2(2:maxwin+1); pe_plt2 = tmp;
         tmp = L_plt(:,2:maxwin+1) ; L_plt = tmp;
       end
     end

     %if (mod(i,100)==0)
     %    disp([num2str(i), ' pe: ',num2str(pe) ] );
     %end
     if (mod(i,1000)==0)
%
        disp( [ num2str(i) ,...
            ' <e^Wd>P: ' ,num2str( mean(PP0) ),...
            ' <e^Wd>F: ' ,num2str( mean(PP1) ),...
            ' <e^Wd>R: ' ,num2str( mean(PP2) ) ] );
        disp( ['corr',...
            ' <e^Wd>P: ', num2str( mean(PP0.*(pe_plt1)) ),...
            ' <e^Wd>F: ', num2str( mean(PP1.*pe_plt1) ),...
            ' <e^Wd>R: ', num2str( mean(PP2.*pe_plt2) ) ] ) ;
        
            disp([    ' AR: ', num2str(AR),...
            ' L(1): ',num2str(L(1)),...
            ' L(2): ',num2str(L(2)),...
            ' L(3): ',num2str(L(3)),...
            ' S:', num2str(S) ]);
        disp(['    ratio: ',num2str(ratio),...
               ' pe: ',num2str(pe), ' ratio*pe: ', num2str(ratio*pe)])
%       disp([num2str(ratio1),'  ',num2str(ratio)  ]);
%       disp(['ewT:  ', num2str(ew_T)]);
        h_e.sim=h_e.sim+histc(e_log,e_plt);
        e_log=[];
        h_x = h_x + hist(xy_accept(:,1),x_plt);
        h_y = h_y + hist(xy_accept(:,2),y_plt);
        h_xy =h_xy +d2hist(xy_accept,x_plt,y_plt); 
        xy_trials=[];
        spot_plot;
        spotplot_2D;
        figure(10)
        clf;
        title(['\lambda(i):  i= ', num2str(i) ]);
        plot(L_plt(1:2,:)'); legend('1','2');
        drawnow;      
    end
    
%     [ q3, q2, t_23, t_41,W_23,W_41]=...
%     inner_t41_cb_1(q1,inner_parms, parms);
%     q4=q1+(q3-q2);

%keyboard;

    if (rand>=0.5)
        inner_parms=inner_parms_dx;
    else
        inner_parms=inner_parms_dy;
    end
%

%    inner_parms
    [ q3, q2, t_23, t_41,W_23,W_41]=...
     inner_t41_cb_true_rvs(q1,inner_parms, parms);
     q4=q1+(q3-q2);  

    P1 = e_V(q1(1),  q1(2),  parms)   ; 
    P2 = e_V(q2(1),  q2(2),  parms)   ;  
    P3 = e_V(q3(1),  q3(2),  parms)   ; 
    P4 = e_V(q4(1),  q4(2),  parms)   ;  
    e_log=[e_log -log(P1)]            ;     

% acceptance probs ========================

    pp0 = P3/P1*prod(t_41./t_23); pp0=1/pp0;
    pp1 = P2/P1;   % pp1=1/pp1;
    pp2 = P4/P3;   % pp2=1/pp2;
%   pp2 = pp0;
                
    PPT_0    = [PP0 pp0]; 
    PPT_1    = [PP1 pp1]; 
    PPT_2    = [PP2 pp2];
    
    if (i>500)
    %L=[0 0 0]';
    %[pe,L,S] = find_pe_3parm(PPT_0,PPT_1,PPT_2,L);
    
     if (pe_flag==1)
        %[pe,L,S] = find_pe_3parm(PPT_0,PPT_1,PPT_2,L);
        
        L=[L(1) L(2)]';
        [pe,L,S,fail]=find_pe_2parm(PPT_1,PPT_2,L);
        if(fail==1)
            disp([ 'fail at: ' num2str(i), 'L= ', num2str(L(1)),...
                ' ',num2str(L(2))]);
            failcount=failcount+1;
            if (failcount==5)
             disp(['trying new lambda:',...  
             num2str(L(1)),' ',num2str(L(2))]);
%         
              failcount=0;
              L=L+(rand(2,1)-0.5)*1e-2;
              disp(['    new lambda:',...  
             num2str(L(1)),' ',num2str(L(2))]);           
            end
        else
            failcount=0;
        end
        L=[L(1) L(2) 0]';  
        pe1=pe; pe2=pe;      
        
%         [pe1,L1,S]=find_pe(PPT_1,L1);    
%         [pe2,L2,S]=find_pe(PPT_2,L2);
%         %L=[L 0 0]';
%         L=[L1 L2 0]';
%         pe=pe2/pe1;
     end

    else
        pe=1;
    end
%   ratio=P3/P1*prod(t_41)/prod(t_23);
%   pe=1;L=1;
    if(pe_flag==0)
       pe=1;
    end
    if(pe==0)
       pe=1e-10;
    end

    switch accept_type
        case 1
         ratio  = P3/P1 * prod(t_41)/prod(t_23);
         pe_tmp = exp( L(2)/P3 - L(1)/P1 );
         %ratio  = ratio * pe_tmp;
    %     ratio=ratio*pe_tmp;
        case 2
         ratio=(P3+P2)/(P1+P4);
         %pe_tmp=exp(L(2)/P3-L(1)/P1);
         %ratio=ratio/pe;
        case 3
         ratio=P3/P1;
    end
    %ratio=(P3+P2)/(P1+P4);
    %ratio=ratio*pe_tmp;
    
    %ratio=ratio/pe;
%=========================================

    accept=0.0;
    if (ratio >=1)
        accept=1;
    elseif ( rand < ratio )
         accept=1; 
    end
%
    if (accept==1)
         q1_0=q1;          q1=q3;
         N_acc = N_acc+1; AR=N_acc/i;
    end   
%   xy_trials=[xy_trials; q3(1) q3(2)];
    xy_accept=[xy_accept; q1(1) q1(2)];
end

%disp('hit enter for plots');
%keyboard;

%% statistics
C   = cov(x_pmf,x_plt);
S.x = [C(1,2) C(1,1) C(2,2)];
C   = cov(y_pmf,y_plt);
S.y = [C(1,2) C(1,1) C(2,2)];

x_pmf = h_x/sum(h_x)/dx;
y_pmf = h_y/sum(h_y)/dx;
xy_pmf=h_xy/sum(sum(h_xy))/dx^2;
h_e.true=h_e.true/sum(h_e.true)/de;
h_e.sim=h_e.sim/sum(h_e.sim)/de;
%
HE=[e_plt'; h_e.true'; h_e.sim'];
HQ=[x_plt' y_plt' pmf.x' pmf.y' x_pmf' y_pmf'];
%keyboard;
H2D=[f xy_pmf];
%H2D=[f; xy_pmf];
%; ]; 

Lx=length(x_plt); Ly=length(y_plt);
chi_sq(1)  = mean( ( pmf.x-x_pmf ).^2 );
chi_sq(2)  = mean( ( pmf.y-y_pmf ).^2 );
chi_sq(3) = 1/Lx/Ly * sum(sum((xy_pmf-f').^2));
chi_sq(4) = mean((h_e.true-h_e.sim).^2);



%% plotting and saving

spot_plot;
spotplot_2D;
figure(200)
plot(e_plt, h_e.true, e_plt,h_e.sim);

%keyboard;
mkdir(fname)
print('-f500','-dpng','-r300',[fname,'\',fname,'-pmfs.png']);
print('-f600','-dpng','-r300',[fname,'\',fname,'-main.png']);
print('-f602','-dpng','-r300',[fname,'\',fname,'-2d_t.png']);
print('-f602','-dpng','-r300',[fname,'\',fname,'-2d_s.png']);
print('-f603','-dpng','-r300',[fname,'\',fname,'-2d_d.png']);

%keyboard;
disp(['saving: ', fname])
%
file=[fname,'/',fname ,'.chi_sq'];
save(file,'-ascii', 'chi_sq')
file=[fname,'/',fname ,'.he'];
save(file,'-ascii', 'HE')
file=[fname,'/',fname ,'.hq'];
save(file,'-ascii', 'HQ')
file=[fname,'/',fname ,'.h2d'];
save(file,'-ascii', 'H2D')
%
save( [fname,'/',fname] );