%% The Stochastic potential well with driving and red noise
% This script is a tutorial script for stochastic integration and
% visualization using red noise realized using the method of Bartosch.
% The governing ODEs of the system are:
% x'=v
% v'=-V'(x)-damping*v+f_deterministic+f_stochastic;
% where V(x) is the potential of interest.
% A tilted double well potential is implemented in this code.
%% This section is all the book keeping for initialization
% Initialize the random number generator and seed it.
rng('shuffle');
%%
% Specify the global parameters for the pendulum
L=5;
x0deep=-3;
x0shal=3;
wellslope=2;
dxnum=1e-8;
myV=@(xnow) ((xnow-x0deep).*(xnow-x0shal)).^2+wellslope*xnow;
myVp=@(xnow) (myV(xnow+dxnum)-myV(xnow-dxnum))./(2*dxnum);
damp=8e-1;
myamp=15; myper=2;
myforce=@(t) myamp*sin(2*pi*t/myper);
%% Bi-stable Potential Well
figure(1)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','bold');
wellxs=linspace(-5,5,1001);
plot(wellxs,myV(wellxs),'k-')
grid on
xlabel('x')
ylabel('Potential')

%%
% Specify the time step and stochastic parameters
numsteps=10; % Sets how many time steps between saves.
numouts=200; % Sets the number of outputs.
dt=1e-3; % sets the time step.  Stochastic methods are low order so small time steps are required.
numens=100; % Sets the number of ensemble members.
%% The parameters for the red noise as realized by Bartosch's algorithm
memt=0.2; % Sets the memory of the noise.  Need at least ten time steps in memt.
rho=exp(-dt/memt); rhoc=sqrt(1-rho*rho); % Parameters for Bartosch's method.
%mysig=40; % Set the standard deviation of the noise. This is strong noise
mysig=30; % Set the standard deviation of the noise. This is medium noise
prevrand=zeros(numens,1); % Start g at g0.
%% Initial conditions
x=x0shal*ones(numens,1); % Start near the shallow well location.
% x=1e-4*pi*ones(numens,1); % Start nearthe stable state.
v=zeros(numens,1); % Initially assume zero velocity.
t=0;
%% Set up arrays for storage and store initial state
xs=zeros(numouts+1,numens);
vs=zeros(numouts+1,numens);
ts=zeros(numouts+1,1);
xs(1,:)=x;
vs(1,:)=v;
ts(1)=t;
%% The main work horse loops.  The double loop can be used view to intermediate results
for ii=1:numouts
    for jj=1:numsteps
        nowrand=randn(numens,1)*mysig; % Generate a new bunch of random numbers.
        pertnow=(nowrand*rhoc+rho*prevrand);
        prevrand=pertnow; % Use Bartosch's memory to get the right "memory".
        %%
        % Define the rhs for the DE.
        rhth=v;
        rhv=-myVp(x)-damp*v+pertnow+myforce(t)*ones(numens,1);
        %%
        % Step the DE forward. If g had a white noise component this would
        % get trickier to account for the sqrt(dt) in Ito's calculus.
        x=x+dt*rhth;
        v=v+dt*rhv;
        t=t+dt;
    end
    xs(ii+1,:)=x;
    vs(ii+1,:)=v;
    ts(ii+1)=t;
end
%% Now that we have the data we can make a variety of pictures
figure(2)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(ts,xs(:,1))
xlabel('t')
ylabel('x')
title('one realization')
subplot(2,1,2)
plot(ts,vs(:,1))
xlabel('t')
ylabel('v')
%%
figure(3)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(xs(:,1),vs(:,1))
ylabel('v')
xlabel('x')
title('one phase space plot')
subplot(2,1,2)
plot(mean(xs,2),mean(vs,2))
ylabel('v')
xlabel('x')
title('mean phase space plot')
%%
figure(4)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(ts,std(xs,[],2))
ylabel('standard deviation of x')
xlabel('t')
subplot(2,1,2)
plot(ts,std(vs,[],2))
ylabel('standard deviation of v')
xlabel('t')
%%
figure(5)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','bold');
subplot(2,1,1)
mnth=min(xs(end,:));
mxth=max(xs(end,:));
bins=linspace(mnth,mxth,floor(numens/10)+1);
hist(xs(end,:),bins)
title('histogram of xs final time')
subplot(2,1,2)
mnv=min(vs(end,:));
mxv=max(vs(end,:));
bins=linspace(mnv,mxv,floor(numens/10)+1);
hist(vs(end,:),bins)
title('histogram of vs final time')
%% Waiting Time Computation
% We now compute waiting times.  These are a standard part of analysis for
% many stochastic systems.  In this case the waiting times can be used to
% characterize the phenomenon of stochastic resonance (basically by looking
% at the histogram of waiting times and seeing if multiples of the period
% of the deterministic driving are evident).  The waiting time is defined
% as the time it takes to move from one well to the other.  That means we
% can use the unstable equilibrium point between the two wells to calculate
% the waiting time.
waits=[];
%%
% Set how close to the unstable equilibrium is close enough.
myeps=0.1;
%%
% Now find the unstable equilibrium point.
xnow=linspace(-1,1,1001);
vv=myV(xnow);
[aa bb]=max(xnow);
xtop=xnow(bb);
for ii=1:numens
    xnow=xs(:,ii);
    crossings=find(abs(xnow-xtop)<myeps);
    %%
    % Find the crossings and make sure the set isn't empty.
    if isempty(crossings)==0
        waitnow=[];
        waitnow(1)=crossings(1);
        dummy=diff(crossings);
        tooshorts=find(dummy<2);
        %%
        % Only keep data if we're not fidgeting around the equilibrium.
        dummy(tooshorts)=[];
        waitnow=[waitnow;dummy];
        waits=[waits;waitnow];
    end
end
waits=waits*dt*numsteps; % Scale.
%%
figure(6)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','bold');
hist(waits,floor(numens/10)+1)
title('histogram of waiting times')


