Added main_bimodal.m for sampling bimodal distribution

main_bimodal.m

@@ -0,0 +1,273 @@
+function main_bimodal(option,varargin)
+
+switch option
+    case 'plot_only'
+        plot_bimodal(varargin{1},[2]);
+        return
+    case 'analysis_only'
+        post_process_bimodal(varargin{1});
+        return
+    case 'sample_path'
+        gen_sample_path_bimodal;
+    otherwise
+        exit_time;
+end
+end
+
+
+function gen_sample_path_bimodal(~)
+
+%p.ntrials = 8*3; %for parallel processing,use multiples of 8
+p.target = 'bimodal';
+p.sigma = 0.32;
+p.dt = 1e-3; %1e-2 is no good for calculating Riesz derivative
+T = 1e4;
+
+%AIM: 1. Sample trajectory. 2. histogram
+
+a = [1.2 2];
+p.location = 2.5; %modal location
+solver = {'Hamiltonian2','Langevin'}; %'Hamiltonian','Underdamped'};
+%for accurate results don't need the other two... too slow
+
+X = zeros(T/p.dt,length(solver),length(a)); %save 1 example is ok
+
+disp('Generating raw data...')
+for i = 1:length(a) %group fractional/nonfractional together as this is less familiar to people
+    aa = a(i);
+    tic
+    for j = 1:length(solver)
+        p.methods.solver = solver{j};
+
+        disp('Simulating...')
+        %parfor k = 1:p.ntrials
+        [temp,t] = fHMC(T,aa,p);
+        %end
+
+        clc
+        fprintf('Solver: %u/%u\n',j,length(solver))
+        fprintf('alpha: %u/%u\n',i,length(a))
+        toc
+
+        X(:,j,i) = temp;
+    end
+
+end
+
+n = floor(T/p.dt); %number of samples
+t = (0:n-1)*p.dt;
+
+save main_bimodal_raw_data_sample_path.mat T a solver X t p
+end
+
+function exit_time()
+
+%x0 = linspace(0.5,2.5,21); %locations of bimodal distribution
+%a = [1.2 2];
+%p.T = 1e3;
+
+x0 = linspace(0.5,1.5,21);
+a = [2]; %do gaussian only
+p.T = 1e4; %do more time since time is longer
+
+solver = {'Hamiltonian2','Langevin'};
+p.target = 'bimodal';
+p.ntrials = 8*3;
+p.dt = 1e-3;
+p.sigma = 0.32;
+
+flag = true(size(a));
+tic
+
+
+
+binedge = linspace(-5,5,101);
+dx = (binedge(2)-binedge(1))*0.5;
+bin = binedge(2:end) - dx;
+
+H = zeros(length(bin), length(x0),length(solver),length(a)); %histogram
+T = zeros(length(x0),length(solver),length(a),p.ntrials);
+
+tic
+disp('Beginning simulation...')
+for m = 1:length(solver)
+    p.methods.solver = solver{m};
+
+    for i = 1:length(x0) %modal separation
+        p.location = x0(i); %2.5;
+
+        for k = 1:length(a)
+
+            if ~flag(k)
+                T(i,m,k,:) = NaN;
+                continue
+            end
+
+            aa = a(k);
+            TT = p.T;
+
+            h = zeros(length(bin), p.ntrials);
+            texit = zeros(1,p.ntrials);
+
+            parfor j = 1:p.ntrials
+                [X,t] = fHMC(TT,aa,p);
+                h(:,j) = histcounts(X,binedge,'normalization','pdf');
+                try
+                    texit(j) = mean([TrapTime(X>0);TrapTime(X<0)]); %mean exit time
+                catch ME
+                    texit(j) = NaN;
+                    disp('Unknown error!')
+                end
+            end
+            H(:,i,m,k) = mean(h,2);
+            T(i,m,k,:) = texit;
+
+            if all(isnan(texit)) %T(k,i,j)>1e5 %if mean trap time is too long, skip subsequent simulations
+                %nan => TrapTime output empty array
+                flag(k) = false;
+            end
+
+            clc
+            fprintf('Solver: %u/%u\n',m,length(solver))
+            fprintf('alpha: %u/%u\n',k,length(a))
+            fprintf('Mode location: %u/%u\n',i,length(x0))
+            toc
+
+        end
+
+    end
+end
+
+%save main_bimodal_mean_exit_time.mat H bin T p x0 a solver
+save main_bimodal_mean_exit_time_gaussian.mat H bin T p x0 a solver
+
+end
+
+
+function plot_bimodal(Q,flag)
+
+close all
+clc
+
+if any(flag==0)
+load main_bimodal_mean_exit_time_gaussian.mat T x0 p
+%load exit_time_HMC_vs_FHMC_corrected.mat T x0
+T(T==0)=nan;
+T = mean(T,4)*p.dt;
+x0 = x0*2;x0=x0(:);
+xx = linspace(x0(1),x0(end),101);
+
+f0 = myfigure;
+for i =2:-1:1
+    indx = ~isnan(T(:,i));
+    linestyle = {'--','-'};
+    markerstyle = {'x','o'};
+    try
+        f = fit(x0(indx),T(indx,i),'exp1');
+        yy = f(xx);        
+        plot(xx(yy<100),yy(yy<100),linestyle{i},'color',mycolor(1,i+2),'HandleVisibility','off','linewidth',1)
+        hold on
+    catch
+    end
+    plot(x0(indx),T(indx,i),markerstyle{i},'color',mycolor(1,i+2),'linewidth',1);
+end
+
+load main_bimodal_mean_exit_time.mat T x0 p
+T(T==0)=nan;
+T = mean(T,4)*p.dt;
+x0 = x0*2;x0=x0(:);
+xx = linspace(x0(1),x0(end),101);
+
+for i =2:-1:1
+    f = fit(x0(:),T(:,i,1),'poly1');
+    linestyle = {'--','-'};
+    markerstyle = {'x','o'};
+    plot(xx,f(xx),linestyle{i},'color',mycolor(1,i),'HandleVisibility','off','linewidth',1)
+    hold on
+    %indx = 1:2:length(x0);
+    plot(x0,T(:,i,1),markerstyle{i},'color',mycolor(1,i),'linewidth',1);    
+end
+
+
+subplotmod;
+
+xlim([1 5])
+ylim([0 100])
+xlabel('Modal Separation')
+ylabel('Mean Exit Time')
+
+legend('LD','HD','FLD','FHD','box','off','location','ne')
+offsetAxes;
+export_fig(f0,'figures/fig_main_bimodal_met.pdf','-pdf','-nocrop','-transparent','-painters');
+end
+
+if any(flag==1)
+    %load main_bimodal_raw_data_sample_path.mat
+    f1 = myfigure;
+    ttl = {'Fractional Hamiltonian Dynamics','Fractional Langevin Dynamics','Hamiltonian Dynamics','Langevin Dynamics'};
+
+    Tmax = 2e2;
+    tindx = (1:40:(Tmax/Q.p.dt))+2e4;
+
+    for k = 1:4
+        subplot(4,1,k)
+        [j,i]=ind2sub([2 2],k);
+        plot(Q.t(tindx)-Q.t(tindx(1)),Q.X(tindx,j,i),'.','color',mycolor(1,k),'markersize',1);
+        title(ttl(k))
+
+        if k<4
+            set(gca,'xtick',[],'XColor','none');
+        else
+            xlabel('t')
+        end
+        ylabel('x_t')
+        %subplotmod;
+        box off
+        set(gca,'TickDir','out')
+        %hold on
+        %plot(xx,yy,'k--','linewidth',1.5)        
+        ylim([-4 4])
+        xlim([0 Tmax])
+        set(gca,'linewidth',1);
+
+    end
+    pos =get(f1,'Position');
+    set(gcf,'position',[pos(1) pos(2) pos(3) 600])
+
+    export_fig(f1,'figures/fig_main_bimodal_sample_path.pdf','-pdf','-nocrop','-transparent','-painters');
+end
+
+if any(flag ==2) %Histogram
+    f2 = myfigure;
+    gs1 = makedist('normal','mu',2.5,'sigma',0.32);
+    gs2 = makedist('normal','mu',-2.5,'sigma',0.32);
+
+    xx = linspace(-4,4,101);
+    yy = 0.5*pdf(gs1,xx)+0.5*pdf(gs2,xx);
+
+    binedge = linspace(-4,4,101);
+    bin = (binedge(2)-binedge(1))*0.5 + binedge(1:end-1);
+
+    for k = 1:4
+        subplot(4,1,k)
+        [j,i]=ind2sub([2 2],k);
+        hc = histcounts(Q.X(:,j,i),binedge,'normalization','pdf');
+        %bar(Q.bin,Q.H(:,j,i),1,'FaceColor',mycolor(1,k),'EdgeColor',mycolor(1,k),'FaceAlpha',0.5);
+        bar(bin,hc,1,'FaceColor',mycolor(1,k),'EdgeColor',mycolor(1,k),'FaceAlpha',0.5);hold on
+        set(gca,'xtick',[])
+        set(gca,'ytick',[])
+        hold on
+        plot(xx,yy,'k','linewidth',1)        
+        if k<0;%k <2.5
+            ylim([0 1.5/2])
+        else
+            ylim([0 1.5])
+        end
+        xlim([-4 4])
+        set(gcf,'position',[100 100 200 600])
+    end
+    %export_fig(f2,'figures/fig_main_unimodal_postprocess_histogram.pdf','-pdf','-nocrop','-transparent','-painters');
+    print(gcf, '-dpdf', 'figures\fig_main_bimodal_histogram.pdf'); 
+end
+
+end

mains/subroutines/fHMC.m

@@ -1,4 +1,4 @@
-function [X,t] = fHMC(T,a,p)
+function [X,t,V] = fHMC(T,a,p)
 %Levy Monte Carlo
 %INPUT: T = time span (s), a = Levy characteristic exponent
 % p = other parameters
@@ -22,7 +22,7 @@
 x = p.x0; %initial condition
 r = p.r0;
 X = zeros(n,1); %position
-%R = zeros(n,1); %momentum
+V = zeros(n,1); %momentum
 
 switch p.target
     case 'unimodal'
@@ -96,6 +96,26 @@
                 end
             end
             X(i) = x;
+            V(i) = r;
+        end
+    case 'Experimental'
+        for i = 1:n
+            while true
+                %diffusion term
+                gam = 1;
+                beta = 1;
+                g = stblrnd(a,0,1,0); %Levy r.v. wt expo a and scale 1                
+                xnew = x + gam*ba(x)*dt + beta*r*dt + (gam^(1/a))*g*dta;
+                rnew = r + beta*ba(x)*dt;
+                %accept criterion
+                if abs(xnew)<p.xmax %5
+                    x = xnew;
+                    r = rnew;
+                    break %accept
+                end
+            end
+            X(i) = x;
+            V(i) = r;
         end
     case 'Underdamped'
         sas = makedist('Stable','alpha',1/a,'beta',0,'gam',1,'delta',0);