function sampler(NO,NI,NW,CT,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   ;
num_trials=2;
%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
%        
      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),...
       '-NW-',num2str(NW),...
       '-CT-',num2str(CT),...
       '-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 CT];
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(['eW:  ',num2str(eWtot)] )
       % keyboard;
        xy_accept=[];
     end

%----Generate large trial move:
     DXtot = (rand-0.5)*dx_12*(x_upper-x_lower);
     DYtot = (rand-0.5)*dy_12*(y_upper-y_lower);
     dx_i = DXtot/NW;  dy_i = DYtot/NW;
     %xflag=0; if (rand>0.5), xflag=1;end
     %q3=q1;
     eWtot=1;
     qi=q1;
     for ni=1:NW
      qis=qi+[dx_i dy_i]';    
        %if (xflag==1), q2=q3+[dx_i 0]';,end
        %if (xflag==0), q2=q3+[0 dx_i]';, end
      [ qi, eW]=...
              inner_stdmc(qi,qis,inner_parms, parms);
       eWtot=eWtot*eW;
      %  keyboard;
     end
     q3=qi;
% acceptance probs ========================
    ratio=eWtot;

    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;var
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')
file=[fname,'/',fname ,'.ar'];
save(file,'-ascii', 'AR')

%
save( [fname,'/',fname] );