From b87f3f0fee9401cc33eb3c544fc6c00d4ea67084 Mon Sep 17 00:00:00 2001 From: qiyangku Date: Mon, 30 Nov 2020 15:46:38 +0800 Subject: [PATCH] Added main_bimodal.m for sampling bimodal distribution --- main_bimodal.m | 273 +++++++++++++++++++++++++++++++++++++++ mains/subroutines/fHMC.m | 24 +++- 2 files changed, 295 insertions(+), 2 deletions(-) create mode 100644 main_bimodal.m diff --git a/main_bimodal.m b/main_bimodal.m new file mode 100644 index 0000000..ff8844f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/mains/subroutines/fHMC.m b/mains/subroutines/fHMC.m index 7f8444e..b4c6a0b 100644 --- a/mains/subroutines/fHMC.m +++ b/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)