function sampler(NO,NI,NT,AT,LS,NV,LD,CF)

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

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

numsteps        =  NO ;
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       = 3 ;
lambda          = LD  ; % 1) is true or typical reverse
                        % 2) is concurrent reverse

pe_flag=0;
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 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=inner_parms_dx;
%
end
%

fname=[ 'NO-',num2str(NO),...
       '-NI-',num2str(NI),...
       '-NT-',num2str(NT),...
       '-AT-',num2str(AT),...
       '-LS-',num2str(LS),...
       '-NV-',num2str(NV),...
       '-LD-',num2str(LD),...
       '-CF-',num2str(CF)];
%     
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
%
parms_R2=landscape.OB;

dx=sqrt( Volume/numsteps );
dx=2*dx;  dx=0.5;
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:



M = 0.5*[ 1 0; 0 1 ]; 
dt = 0.1  %what is the scale?
gamma_fric = 0;  bbk_parms{1} = M;
bbk_parms{2} = dt;  bbk_parms{3} = gamma_fric;

hmc_nsteps=NV;
hmc_parms{1}=parms;
hmc_parms{2}=bbk_parms;
hmc_parms{3}=hmc_nsteps;

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

%visualization of potential:
e_log=[]; ew=1;
%keyboard;
 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 CF];
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));
%=============================================
PP0(1)=1;
for i=1:numsteps
     if (mod(i,500)==0)
        disp([num2str(i),' ','<T23>= ' , num2str(mean(pp0)),' AR: ',num2str(AR) ]);
       % h_e.sim=h_e.sim+histc(e_log,e_plt);
        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;
        drawnow;    
        disp(['pcorr:  ',num2str(pcorr)] )
     end
%if (reverse_path==2)     
     [ q3, q2, pcorr]=...
              inner_hmc_lincorr(q1,inner_parms, parms,hmc_parms);
 %    inner_hmc_lambda(q1,inner_parms, parms,hmc_parms,lambda);
     
%end
    P1 = e_V(q1(1),  q1(2),  parms)   ; 
    P2 = e_V(q2(1),  q2(2),  parms)   ;  
    P3 = e_V(q3(1),  q3(2),  parms)   ; 
% acceptance probs ========================

%    pp0 = prod(t_23); %pp0=1/pp0;
    PPT_0    = [PP0 pp0];
 %   w_23=prod(t_41);
    switch accept_type
        case 1          
          ratio=P3/P1;
        case 2 
            ratio=P3/P1*pcorr;
    end
    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] );