function [HistM, food]=acoli_hist(initpos, foodfn, rttp, ttrp, nruns, runlength, xbins, ybins)
% [HistM, food]=acoli_hist(initpos, foodfn, rttp, ttrp, nruns, runlength, xbins, ybins)
%
% Computes the number of times a bacterium enters a region. The regions are
% defined as rectangular bins which cover the space in a regular grid. The
% experiment consists of several trials (nruns) during each of which the
% bacterium is allowed to move for several steps.
%
% Behavior: the bacterium either performs a run -- advancing two units in
% the current direction -- or a tumble -- rotating the current direction by
% a random angle between 20 and 90 degrees and advancing one unit.
% Switching between the two modes is controlled by the two probabilities
% rttp and ttrp. In addition, if the bacterium senses that it has moved in
% the direction of increasing food concentration, then the random angle
% during a tumble is between 10 and 45 degrees, i.e. the bacterium
% does not turn much if it goes toward food.
%
%
% Input
%
% initpos = initial position of the bacterium, e.g., [0 0]
% foodfn = selects food function, must be either 1 or 2
% rttp = run-to-tumble probability, 
% ttrp = tumble-to-run probability,
% nruns = number of trials
% runlength = number of steps taken for each trial
% xbins/ybins = number of bins in the x and y axis, resp.
%
% Output
%
% HistM = matrix of bins containing the number of times a bacterium has
%         passed through the corresponding region, i.e. histogram matrix.
% food = matrix containing the concentration of food, i.e. the food function
 

% run
% t1=clock; [HistM,bloc,Food]=ecoli_hist(20000/2, 100, 30,30);
% elapsed=etime(clock,t1)

showfigs=0;
record_hist = 1;
sp=0.1;
xmin=-2; xmax = 2;
ymin=-2; ymax = 2;
[X,Y] = meshgrid(xmin:sp:xmax, ymin:sp:ymax);

xstp = (xmax-xmin)/xbins;
ystp = (ymax-ymin)/ybins;


xlocs = xmin:xstp:xmax; xstp = (xmax-xmin) / length(xlocs);
ylocs = ymin:ystp:ymax; ystp = (ymax-ymin) / length(ylocs);
HistM = zeros(length(xlocs), length(ylocs));

%%%%%
foodx=0; foody=0;

if foodfn==2
    Food.func = @foodfunc2;
elseif foodfn == 1
    Food.func = @foodfunc1;
else
   error('foodfn must be either 1 or 2.') 
end
    
Food.pts = [0.75, 0.75; -0.5 -0.5];
[XX,YY] = meshgrid(xmin:0.01:xmax, ymin:0.01:ymax);
food = Food.func(XX,YY);

dt = sp;

% bacterium frames = 300
% gradient ascent frames = 40

for runi = 1:nruns
    bloc.pos = initpos; %[-0.5,-0.5]; % this is the bacterium [x,y; vx,vy]
    bloc.oldpos = bloc.pos;
    bloc.vec =  [rand()-0.5, rand()-0.5]; %[1,-1];
    bloc.fgrad=0;
    bloc.mode = 0; % 0 - run; 1 -- tumble
    bloc.hist = bloc.pos; % history of position

    for i=1:runlength
        bloc = update_bacterium(bloc,Food,dt); 
         %   bloc = update_rand_gradascent(bloc,Food,dt);
        px = bloc.pos(1);
        py = bloc.pos(2);
        if px >= xmin & px <= xmax & py >= ymin & py <= ymax
            ptj = floor((px - xmin)/xstp)+1;
            pti = floor((py - ymin)/ystp)+1;
            HistM(pti,ptj) = HistM(pti,ptj) + 1;
        end
        
%         if showfigs~=0
%             display_status(X,Y,Food, bloc);
%             getframe;
%         end
    end
end

% axis xy

if showfigs~=0
    figure(2)
    [C,h]  = contour(X,Y,Food.func(X,Y));
    set(h,'ShowText','on','LabelSpacing',144*5);
    colormap(0.7*[1 1 1])
    hold on
    display_trajectory(bloc);
    plot(Food.pts(1,2), Food.pts(1,2),'*')
    hold off
end
%%
    function c=foodfunc1(x,y)
        s=.5;
        c = exp(-s*((x-Food.pts(1,1)).^2 + (y-Food.pts(1,2)).^2));
        
    end
%%
    function c=foodfunc2(x,y)
        s=5;
        c = exp(-s*((x-Food.pts(1,1)).^2 + (y-Food.pts(1,2)).^2)) +...
            exp(-s*((x-Food.pts(2,1)).^2 + (y-Food.pts(2,2)).^2));
        
    end

%%
    function nv=rotatevec(v,a)
        % rotate vector v by angle a
       rotM = [
           cos(a) sin(a)
           -sin(a) cos(a)
           ];
       
       nv = (rotM * v')';
    end

%% update_bacterium(bloc,Food, dt)
    function newbloc=update_bacterium(bloc,Food, dt)
        runtotumble_prob = 1-rttp; 0.1;%*2;
        tumbletorun_prob = 1-ttrp; 0.3;
        sec1 = 1; % # of steps corresponding to 1 sec
        sec3 = 1; sec1*3; % # of steps corresponding to 3 secs
        
        fdist = Food.func; % food distribution function
        
        % compute concentration gradient
        fgrad = fdist(bloc.pos(1), bloc.pos(2)) ...
            - fdist(bloc.oldpos(1), bloc.oldpos(2));

        newvec = bloc.vec;
        newpos = bloc.pos;
        newbloc = bloc;
        
        actualruntotumble_prob = runtotumble_prob;
        rotupgrad = 1; % power for tumble rotation
        
        if length(bloc.hist) >= sec3+sec1 | 1;
           %cvals1=Food.func(bloc.hist(end-sec1+1:end,1), bloc.hist(end-sec1+1:end,1)); 
           %cvals_prev3=Food.func(bloc.hist(end-sec1-sec3+1:end-sec1,1), bloc.hist(end-sec1-sec3+1:end-sec1,1));

           %temporal_diff = mean(cvals1) - mean(cvals_prev3);
           
           tdiff = Food.func(bloc.pos(1),bloc.pos(2)) - Food.func(bloc.oldpos(1),bloc.oldpos(2));
           %if temporal_diff > 0
           if tdiff > 0
               actualruntotumble_prob = runtotumble_prob; %/(1e1);  %*temporal_diff^2);
               rotupgrad =  2; 1.2;
           else
               actualruntotumble_prob = runtotumble_prob;
           end
        end
        
        % ignoring grad measurement
        %actualruntotumble_prob = runtotumble_prob; rotupgrad = 1;

        if bloc.mode == 0 % run
            if actualruntotumble_prob < rand()
                newbloc.mode = 1; % change to tumble mode
            end
                       
            vec = bloc.vec / norm(bloc.vec);
%            newvec = 2*rotatevec(vec, rand()*8 / (2*pi)); % two units
            newvec = 2*vec;%2*rotatevec(vec, 0*rand()*8 / (2*pi)); % two units
        
        end

        if bloc.mode == 1 % tumble
           if tumbletorun_prob < rand()
               newbloc.mode = 0; % change to run mode
           end
           
           %%%%
           %newbloc.mode = 0; % change to run mode
           %%%%

           vec = bloc.vec / norm(bloc.vec);
           rnum = rand();
           sgn = 1;sign(rand() - 0.5);
           newvec = 1*rotatevec(vec, sgn*(20*(1-rnum) + rnum*90)/rotupgrad / 360*(2*pi)); % one unit
        end
        
%         if fgrad < 0.0 || fgrad*(0.1*rand()) > 0.1
%             angle = pi/4 * (rand());% * (1-abs(fgrad));
%             newvec=rotatevec(bloc.vec,angle);                        
%         end
        
        newpos = bloc.pos + dt*newvec;
                
        newbloc.oldpos = bloc.pos;
        newbloc.pos = newpos;
        newbloc.vec = newvec;
        newbloc.fgrad = fgrad;
        
        if record_hist ~=0
            newbloc.hist = [bloc.hist; newpos];
        end
    end

%%
    function display_trajectory(bloc)
        % draw trajectory
        plot(bloc.hist(:,1), bloc.hist(:,2), ':');

        % draw end location
        plot(bloc.pos(1),bloc.pos(2), 'ko','MarkerFaceColor','r')

        % draw initial location
        plot(bloc.hist(1,1),bloc.hist(1,2), 'ks','MarkerFaceColor','k')

        axis equal tight
        axis([xmin xmax ymin ymax])

    end
end