diff --git a/ATLAS/Learning_ini_chart.m b/ATLAS/Learning_ini_chart.m index 00c01ca..a38b207 100644 --- a/ATLAS/Learning_ini_chart.m +++ b/ATLAS/Learning_ini_chart.m @@ -1,46 +1,46 @@ -% Find the chart for randomly sampled initial points on the manifold. -%These initial points are evenly sampled from a long trajectory - -chart = cell(1,K_int); +% Find the chart for randomly sampled initial points on the manifold. +% These initial points are evenly sampled from a long trajectory. + +chart = cell(1, params.learning.K_int); Sim_parameter = struct( ... - 'T_max', T_one, ... - 'dt', dt , ... - 'gap', 1, ... - 'X_int', X_int ... + 'T_max', params.learning.T_one, ... + 'dt', params.RHS.dt , ... + 'gap', 1, ... + 'X_int', params.learning.X_int ... ); - - -disp(['Initial point: ', num2str(X_int), ' T_max: ', num2str(Sim_parameter.T_max)]) -% Simulate a long trajctory starting from this initial condition. -tic - -[output] = simulator_one(Sim_parameter); -t1=toc; -disp(['one long trajctory is simulated. The time spent is ', num2str(t1/60), ' mins.' ]) - - -% evenly sample N initial points from the trajctory +disp(['Initial point: ', num2str(params.learning.X_int), ' T_max: ', num2str(Sim_parameter.T_max)]) + +% Simulate a long trajectory starting from this initial condition. +tic + +[output] = params.simulator.one(Sim_parameter); +t1=toc; +disp(['one long trajectory is simulated. The time spent is ', num2str(t1/60), ' mins.' ]) + + + +% evenly sample N initial points from the trajectory output = output(10000:end,:); -[L, ~] = size(output); -X_int_sample = output(1: round(L/K_int):end, :); +[L, ~] = size(output); +X_int_sample = output(1: round(L/params.learning.K_int):end, :); tic -disp(['Parallelly run chart simulator starting from ', num2str(K_int) ,' initial points. ']) -for i = 1 : K_int +disp(['Parallelly run chart simulator starting from ', num2str(params.learning.K_int) ,' initial points. ']) +for i = 1 : params.learning.K_int disp(['Landmark',num2str(i), ' is learned.']) Sim_parameter = struct( ... - 'T_max', T_max, ... - 'dt', dt , ... - 'N', N, ... + 'T_max', params.RHS.T_max, ... + 'dt', params.RHS.dt , ... + 'N', params.RHS.N, ... 'X_int', X_int_sample(i,:) ... ); - [data, Cov_store, Mean_store] = simulator_par(Sim_parameter); - [chart{i}] = Learning_Slow_Manifold( data, Cov_store, Mean_store, D, d,modify ); - - + [data, Cov_store, Mean_store] = params.simulator.parallel(Sim_parameter); + [chart{i}] = Learning_Slow_Manifold( data, Cov_store, Mean_store, params.RHS.D, params.RHS.d, params.RHS.modify ); + + end t2 =toc; disp(['Initial learning stage is completed. The time spent is ', num2str(t2/60), ' mins.' ]) diff --git a/ATLAS/MSM_learning.m b/ATLAS/MSM_learning.m index 4af7649..a2ed28a 100644 --- a/ATLAS/MSM_learning.m +++ b/ATLAS/MSM_learning.m @@ -1,9 +1,9 @@ disp('Starting MSM learning stage') - [TranM] = MSM(chart, connectivity, MSM_parameter, weighted_dd2, d); - step = MSM_parameter.step; - N_state = MSM_parameter.N_state; - dt_s = MSM_parameter.dt_s; - save( TranM_fileName ,'TranM','step','N_state','dt_s') + [TranM] = MSM(chart, connectivity, params.MSM, weighted_dd2, params.RHS.d); + step = params.MSM.step; + N_state = params.MSM.N_state; + dt_s = params.MSM.dt_s; + save( params.paths.TranM ,'TranM','step','N_state','dt_s') disp('MSM learning stage is completed') \ No newline at end of file diff --git a/ATLAS/exploration_learning.m b/ATLAS/exploration_learning.m index 0ce87c6..fdcfb35 100644 --- a/ATLAS/exploration_learning.m +++ b/ATLAS/exploration_learning.m @@ -1,33 +1,32 @@ - disp('Starting exploration learning stage') % set to exploration mode mode = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh,nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, params.learning.option,mode ); - %% Use the exploration mode to further Learn + %% Use the exploration mode to further Learn explore_round = 1; - explore_round_max = 100; + explore_round_max = 100; K = length(chart); - random_start; + random_start; while bin_N~=1 && explore_round0)/N_IC * (-Q1(1,1,k)); - Q1(1,3,k) = sum(exit_curr1_ori<0)/N_IC * (-Q1(1,1,k)); - Q1(2,1,k) = sum(exit_curr2_ori>0)/N_IC * (-Q1(2,2,k)); - Q1(2,3,k) = sum(exit_curr2_ori<0)/N_IC * (-Q1(2,2,k)); - Q1(3,1,k) = sum(exit_curr3_ori<0)/N_IC * (-Q1(3,3,k)); - Q1(3,2,k) = sum(exit_curr3_ori>0)/N_IC * (-Q1(3,3,k)); + Q1(1,2,k) = sum(exit_curr1_ori>0)/params.MFPT.N_IC * (-Q1(1,1,k)); + Q1(1,3,k) = sum(exit_curr1_ori<0)/params.MFPT.N_IC * (-Q1(1,1,k)); + Q1(2,1,k) = sum(exit_curr2_ori>0)/params.MFPT.N_IC * (-Q1(2,2,k)); + Q1(2,3,k) = sum(exit_curr2_ori<0)/params.MFPT.N_IC * (-Q1(2,2,k)); + Q1(3,1,k) = sum(exit_curr3_ori<0)/params.MFPT.N_IC * (-Q1(3,3,k)); + Q1(3,2,k) = sum(exit_curr3_ori>0)/params.MFPT.N_IC * (-Q1(3,3,k)); Q2(1,1,k) = -1/mean(FPT_sim1); Q2(2,2,k) = -1/mean(FPT_sim2); Q2(3,3,k) = -1/mean(FPT_sim3); - Q2(1,2,k) = sum(exit_curr1_sim>0)/N_IC * (-Q2(1,1,k)); - Q2(1,3,k) = sum(exit_curr1_sim<0)/N_IC * (-Q2(1,1,k)); - Q2(2,1,k) = sum(exit_curr2_sim>0)/N_IC * (-Q2(2,2,k)); - Q2(2,3,k) = sum(exit_curr2_sim<0)/N_IC * (-Q2(2,2,k)); - Q2(3,1,k) = sum(exit_curr3_sim<0)/N_IC * (-Q2(3,3,k)); - Q2(3,2,k) = sum(exit_curr3_sim>0)/N_IC * (-Q2(3,3,k)); + Q2(1,2,k) = sum(exit_curr1_sim>0)/params.MFPT.N_IC * (-Q2(1,1,k)); + Q2(1,3,k) = sum(exit_curr1_sim<0)/params.MFPT.N_IC * (-Q2(1,1,k)); + Q2(2,1,k) = sum(exit_curr2_sim>0)/params.MFPT.N_IC * (-Q2(2,2,k)); + Q2(2,3,k) = sum(exit_curr2_sim<0)/params.MFPT.N_IC * (-Q2(2,2,k)); + Q2(3,1,k) = sum(exit_curr3_sim<0)/params.MFPT.N_IC * (-Q2(3,3,k)); + Q2(3,2,k) = sum(exit_curr3_sim>0)/params.MFPT.N_IC * (-Q2(3,3,k)); if k ==1 - disp(['The original simulator at trans state has ', num2str(mean(FPT_ori1),'%5.4f'),'+-',num2str(std(FPT_ori1)*1.96/sqrt(N_IC),'%5.4f')]) - disp(['The original simulator at cis-top state has ', num2str(mean(FPT_ori2),'%5.4f'),'+-',num2str(std(FPT_ori2)*1.96/sqrt(N_IC),'%5.4f')]) - disp(['The original simulator at cis-bottom state has ', num2str(mean(FPT_ori3),'%5.4f'),'+-',num2str(std(FPT_ori3)*1.96/sqrt(N_IC),'%5.4f')]) - disp(['The ATLAS simulator at trans state has ', num2str(mean(FPT_sim1),'%5.4f'),'+-',num2str(std(FPT_sim1)*1.96/sqrt(N_IC),'%5.4f')]) - disp(['The ATLAS simulator at cis-top state has ', num2str(mean(FPT_sim2),'%5.4f'),'+-',num2str(std(FPT_sim2)*1.96/sqrt(N_IC),'%5.4f')]) - disp(['The ATLAS simulator at cis-bottom state has ', num2str(mean(FPT_sim3),'%5.4f'),'+-',num2str(std(FPT_sim3)*1.96/sqrt(N_IC),'%5.4f')]) + disp(['The original simulator at trans state has ', num2str(mean(FPT_ori1),'%5.4f'),'+-',num2str(std(FPT_ori1)*1.96/sqrt(params.MFPT.N_IC),'%5.4f')]) + disp(['The original simulator at cis-top state has ', num2str(mean(FPT_ori2),'%5.4f'),'+-',num2str(std(FPT_ori2)*1.96/sqrt(params.MFPT.N_IC),'%5.4f')]) + disp(['The original simulator at cis-bottom state has ', num2str(mean(FPT_ori3),'%5.4f'),'+-',num2str(std(FPT_ori3)*1.96/sqrt(params.MFPT.N_IC),'%5.4f')]) + disp(['The ATLAS simulator at trans state has ', num2str(mean(FPT_sim1),'%5.4f'),'+-',num2str(std(FPT_sim1)*1.96/sqrt(params.MFPT.N_IC),'%5.4f')]) + disp(['The ATLAS simulator at cis-top state has ', num2str(mean(FPT_sim2),'%5.4f'),'+-',num2str(std(FPT_sim2)*1.96/sqrt(params.MFPT.N_IC),'%5.4f')]) + disp(['The ATLAS simulator at cis-bottom state has ', num2str(mean(FPT_sim3),'%5.4f'),'+-',num2str(std(FPT_sim3)*1.96/sqrt(params.MFPT.N_IC),'%5.4f')]) end diff --git a/model/Butane/butane_ori_simulation.m b/model/Butane/butane_ori_simulation.m index 9bbbb48..d7c6963 100644 --- a/model/Butane/butane_ori_simulation.m +++ b/model/Butane/butane_ori_simulation.m @@ -3,12 +3,12 @@ N_IC = 30000; Sim_parameter = struct( ... 'T_max', T_one*5, ... - 'dt', dt , ... + 'dt', params.RHS.dt , ... 'gap', 1, ... - 'X_int', X_int ... + 'X_int', params.learning.X_int ... ); - [output] = simulator_one(Sim_parameter); + [output] = params.simulator.one(Sim_parameter); output = output(1000:end,:); phi_ori = atan2(output(:,6),output(:,4)); disp(['one single traj is simulated with the original simulator']) - \ No newline at end of file + diff --git a/model/Butane/butane_relearning.m b/model/Butane/butane_relearning.m index f0b0847..de78931 100644 --- a/model/Butane/butane_relearning.m +++ b/model/Butane/butane_relearning.m @@ -1,27 +1,27 @@ -% previous relearning is to make sure it relaxes onto the manifold. -% Now we set N 10 times bigger. -% The upperbound and lowerbound are set to [10dt, 15dt]. This is the -% relaxation timescale. +% previous relearning is to make sure it relaxes onto the manifold. +% Now we set N 10 times bigger. +% The upper bound and lower bound are set to [10dt, 15dt]. This is the +% relaxation timescale. - RHS_parameter_relearn.T_max = 15*dt; - RHS_parameter_relearn.UpperBound = 15*dt; - RHS_parameter_relearn.LowerBound = 10*dt; - RHS_parameter_relearn.modify = 1; + params.relearn.RHS.T_max = 15*params.RHS.dt; + params.relearn.RHS.UpperBound = 15*params.RHS.dt; + params.relearn.RHS.LowerBound = 10*params.RHS.dt; + params.relearn.RHS.modify = 1; relearn_parameter = struct( ... - 'N_relearn', 10*N, ... + 'N_relearn', 10*params.RHS.N, ... 'iter', 2, ... - 'RHS_parameter', RHS_parameter_relearn, ... + 'RHS_parameter', params.relearn.RHS, ... 'relative_threshold', [0.005, 0.02], ... - 'simulator_par', @(Sim_parameter) simEuler_par(RHS_parameter_relearn,Sim_parameter) ... - ); + 'simulator_par', @(Sim_parameter) simEuler_par(params.relearn.RHS,Sim_parameter) ... + ); disp('Starting Relearning stage') - K = length(chart); + K = length(chart); index_to_learn = 1:K; [ chart,~ ] = relearn_chart(chart, relearn_parameter, index_to_learn); - [chart, connectivity, P] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - save( chart_fileName ,'chart','connectivity','P') - + [chart, connectivity, P] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + save( params.paths.chart ,'chart','connectivity','P') + disp('Relearning stage is completed.') - \ No newline at end of file + diff --git a/model/Butane/butane_simulation_learning.m b/model/Butane/butane_simulation_learning.m index 33d178e..43c3d41 100644 --- a/model/Butane/butane_simulation_learning.m +++ b/model/Butane/butane_simulation_learning.m @@ -3,22 +3,22 @@ for i = 1:500 disp(['traj No.',num2str(i) ]) K = length(chart); - chart_sim_parameter.nearest = max(round(rand * K),1); - chart_sim_parameter.connectivity = connectivity; - chart_sim_parameter.X_int = chart{chart_sim_parameter.nearest}.X_int; - chart_sim_parameter.Nstep = 2*10^3; - [~, ~, chart] = ATLAS_simulator2(weighted_dd2, chart_sim_parameter,RHS_parameter, simulator_par, chart); - [chart, connectivity, ~] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - if relearn_option == 1 - K_next = length(chart); + params.chart_sim_parameter.nearest = max(round(rand * K),1); + params.chart_sim_parameter.connectivity = connectivity; + params.chart_sim_parameter.X_int = chart{params.chart_sim_parameter.nearest}.X_int; + params.chart_sim_parameter.Nstep = 2*10^3; + [~, ~, chart] = ATLAS_simulator2(weighted_dd2, params.chart_sim_parameter,params.RHS, params.simulator.parallel, chart); + [chart, connectivity, ~] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + if params.relearn.option == 1 + K_next = length(chart); if K_next>K index_to_learn = K+1: K_next; - [ chart ] = relearn_chart(chart, relearn_parameter, index_to_learn); - [chart, connectivity, P, bin_N ] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - save( chart_fileName ,'chart','connectivity','P') + [ chart ] = relearn_chart(chart, params.relearn.settings, index_to_learn); + [chart, connectivity, P, bin_N ] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + save( params.paths.chart ,'chart','connectivity','P') end end - - save( chart_fileName ,'chart','connectivity','P') - end + + save( params.paths.chart ,'chart','connectivity','P') + end disp('Simulation stage is completed') diff --git a/model/Butane/set_parameter.m b/model/Butane/set_parameter.m index aa7a88b..e25d2e9 100644 --- a/model/Butane/set_parameter.m +++ b/model/Butane/set_parameter.m @@ -1,33 +1,40 @@ +function params = set_parameter() +%SET_PARAMETER Configure simulation, learning, and analysis settings for Butane. +% +% params = SET_PARAMETER() returns a struct with the configuration that was +% previously populated in the base workspace. The struct is organised into +% topical fields (RHS, simulator, learning, chart_sim_parameter, relearn, +% paths, etc.) so downstream workflows can reference the data explicitly. + + datapath = evalin('base','datapath'); + D = 6; % The dimension of the full system - d = 1; % The dimension of the slow manifold. Later we can estimate this dimenstion - chi_p = 1.96; % for 95% Confidence interval for 1D Gaussian distribution + d = 1; % The dimension of the slow manifold. Later we can estimate this dimension + chi_p = 1.96; % for 95% Confidence interval for 1D Gaussian distribution K_int = 2; % Number of initial points on the manifold. - N = 1600000; % Number of short traj starting from the same initial point - dt = 1*10^-6; % Simulation time-step - T_one = 10; % Time for single long trajectory to learn ini K charts. + N = 1600000; % Number of short trajectories starting from the same initial point + dt = 1*10^-6; % Simulation time-step + T_one = 10; % Time for single long trajectory to learn initial K charts. relearn_option = 1; % option to relearn - - t0 = 10*dt; % The time as the benchmark. Drift and diffusion varies small within t0 - T_max = 50*dt; %50*dt; % Time for short trajctories. + + t0 = 10*dt; % The time as the benchmark. Drift and diffusion varies small within t0 + T_max = 50*dt; % Time for short trajectories. X_int = [ -1.4461 0.5578 1.5361 1.4268 -2.0593 0.1615] + 0.1*randn(1,6); - threshold = [1.5, 0.001]; % The first one is R_max determined by the diameter of the slow manifold when setting the value of R and the second one threshold for occational prj - connectivity_threshold = 4; % two landmarks are connected if they are with in r*sqrt(t0) + threshold = [1.5, 0.001]; % First value R_max, second threshold for occasional projection + connectivity_threshold = 4; % Two landmarks are connected if within r*sqrt(t0) explore_threshold = 0.95; - option = 1; % 1 for fast mode prj, 2 for orthogonal prj - modify = 0; % 0 use the intercept(t=0) 1 use the tau_min point as the landmark - % The training window is larger than the relaxtion here. It is to ensure - % the landmark are sufficiently relaxed onto the manifold - LowerBound = 40*dt; %40*dt; % Lowerbound of the training window - UpperBound = 50*dt; %50*dt; % Upperbound of the training window - - - %parameter for simuation - Nstep = 2*10^2; - gap = 10; - - - - %parameter for RHS + option = 1; % 1 for fast mode projection, 2 for orthogonal projection + modify = 0; % 0 use the intercept(t=0) 1 use the tau_min point as the landmark + % The training window is larger than the relaxation here. It is to ensure + % the landmark is sufficiently relaxed onto the manifold + LowerBound = 40*dt; % Lower bound of the training window + UpperBound = 50*dt; % Upper bound of the training window + + % Parameters for simulation + Nstep = 2*10^2; + gap = 10; + + % Parameters for RHS l = 1.53; k2 = 319225; % 1.17*10^6; k3 = 62500; @@ -37,7 +44,7 @@ c3 = -3227.70; beta = 4*10^-3; sigma = sqrt(2*beta^(-1)); - + parameter = struct( ... 'l', l, ... 'k2', k2, ... @@ -48,10 +55,10 @@ 'c3', c3 , ... 'sigma', sigma ... ); - -drift = @(X) get_drift(X, parameter); -diffusion = @(X) get_diffusion(X, parameter); -RHS_parameter = struct( ... + + drift = @(X) get_drift(X, parameter); + diffusion = @(X) get_diffusion(X, parameter); + RHS_parameter = struct( ... 'dt', dt , ... 'T_max', T_max , ... 'N', N, ... @@ -69,56 +76,90 @@ 'LowerBound', LowerBound ... ); +% Setup simulator + simulator_one = @(Sim_parameter) simEuler_one_traj(RHS_parameter,Sim_parameter); + simulator_no_par = @(Sim_parameter) simEuler(RHS_parameter,Sim_parameter); + simulator_par = @(Sim_parameter) simEuler_par(RHS_parameter,Sim_parameter); -% setup simulator -simulator_one = @(Sim_parameter) simEuler_one_traj(RHS_parameter,Sim_parameter); -simulator_no_par = @(Sim_parameter) simEuler(RHS_parameter,Sim_parameter); -simulator_par = @(Sim_parameter) simEuler_par(RHS_parameter,Sim_parameter); - - -MSM_parameter = struct( ... + MSM_parameter = struct( ... 'step', 1, ... 'N_state', 1*10^6, ... 'dt_s', 10*dt ... ); - -chart_sim_parameter = struct ( ... + chart_sim_parameter = struct ( ... 'X_int', [], ... 'nearest', 1, ... 't0', t0, ... - 'dt_s', 10*dt, ... + 'dt_s', 10*dt, ... 'D', RHS_parameter.D, ... 'd', RHS_parameter.d, ... - 'Nstep', Nstep, ... + 'Nstep', Nstep, ... 'gap', gap, ... 'connectivity', [], ... 'explore_threshold', explore_threshold, ... 'N', N ... - ); - - - % parameter for relearn chart - - RHS_parameter_relearn = RHS_parameter; % one can change some RHS if you believe the invariant manifold is not changed. - relearn_parameter = struct( ... + ); + + % Parameters for relearn chart + + RHS_parameter_relearn = RHS_parameter; % one can change some RHS if you believe the invariant manifold is not changed. + relearn_parameter = struct( ... 'N_relearn', 6*N, ... - 'iter', 15, ... + 'iter', 15, ... 'RHS_parameter', RHS_parameter_relearn, ... 'relative_threshold', [0.01, 0.05], ... 'simulator_par', @(Sim_parameter) simEuler_par(RHS_parameter_relearn,Sim_parameter) ... ); - - % parameter for MFPTlength - N_IC = 30000; - %file path - butane_paths = build_butane_paths(datapath); - chart_fileName = butane_paths.chart; - chart_part_fileName = butane_paths.chart_part; - chart_plot_fileName = butane_paths.chart_plot; - TranM_fileName = butane_paths.TranM; - TranM_plot_fileName = butane_paths.TranM_plot; - chart_relearn_fileName = butane_paths.chart_relearn; - well_fileName = butane_paths.well; - FPT_fileName = butane_paths.FPT; + % Parameters for MFPT length + N_IC = 30000; + + % File paths + chart_fileName = [datapath,'chart.mat']; + chart_part_fileName = [datapath,'chart_part.mat']; + chart_plot_fileName = [datapath,'chart_plot.mat']; + TranM_fileName = [datapath,'TranM.mat']; + TranM_plot_fileName = [datapath,'TranM_plot.mat']; + chart_relearn_fileName = [datapath,'chart_relearn.mat']; + well_fileName = [datapath,'well.mat']; + FPT_fileName = [datapath,'Butane_FPT.mat']; + + params = struct(); + params.RHS = RHS_parameter; + params.simulator = struct( ... + 'one', simulator_one, ... + 'serial', simulator_no_par, ... + 'parallel', simulator_par ... + ); + params.MSM = MSM_parameter; + params.chart_sim_parameter = chart_sim_parameter; + params.relearn = struct( ... + 'option', relearn_option, ... + 'RHS', RHS_parameter_relearn, ... + 'settings', relearn_parameter ... + ); + params.learning = struct( ... + 'K_int', K_int, ... + 'T_one', T_one, ... + 'X_int', X_int, ... + 'explore_threshold', explore_threshold, ... + 'option', option ... + ); + params.simulation = struct( ... + 'Nstep', Nstep, ... + 'gap', gap ... + ); + params.paths = struct( ... + 'datapath', datapath, ... + 'chart', chart_fileName, ... + 'chart_part', chart_part_fileName, ... + 'chart_plot', chart_plot_fileName, ... + 'TranM', TranM_fileName, ... + 'TranM_plot', TranM_plot_fileName, ... + 'chart_relearn', chart_relearn_fileName, ... + 'well', well_fileName, ... + 'FPT', FPT_fileName ... + ); + params.MFPT = struct('N_IC', N_IC); +end diff --git a/model/Halfmoon/halfmoon_MFPT.m b/model/Halfmoon/halfmoon_MFPT.m index 8d3b9d6..454297d 100644 --- a/model/Halfmoon/halfmoon_MFPT.m +++ b/model/Halfmoon/halfmoon_MFPT.m @@ -1,59 +1,58 @@ - Mean_FPT = zeros(2*3,10); + Mean_FPT = zeros(6,10); relative_error_FPT = zeros(4,10); - t_final_ori = zeros(2,1); t_final = zeros(4,10); - - theta_s = atan(b2^2/(2*b1)); + + theta_s = atan(params.RHS.parameter.b2^2/(2*params.RHS.parameter.b1)); well_threshold = [-2*atan(4-sqrt(15)), -2*atan(4+sqrt(15))+2*pi]; - [FPT_ori1, t_final_ori(1)] = MFPT_halfmoon_ori(well_threshold, output, phi_ori,N_IC,RHS_parameter,chart_sim_parameter, nonlin_trans_inv); + [FPT_ori1, t_final_ori(1)] = MFPT_halfmoon_ori(well_threshold, output, phi_ori,params.MFPT.N_IC,params.RHS,params.chart_sim_parameter, params.RHS.nonlin_trans_inv); - well_threshold = [-2*atan(4+sqrt(15)), -2*atan(4-sqrt(15))]; - [FPT_ori2, t_final_ori(2)] = MFPT_halfmoon_ori(well_threshold, output, phi_ori,N_IC,RHS_parameter,chart_sim_parameter, nonlin_trans_inv); + well_threshold = [-2*atan(4+sqrt(15)), -2*atan(4-sqrt(15))]; + [FPT_ori2, t_final_ori(2)] = MFPT_halfmoon_ori(well_threshold, output, phi_ori,params.MFPT.N_IC,params.RHS,params.chart_sim_parameter, params.RHS.nonlin_trans_inv); %% ATLAS simulator for k = 1:10 disp([num2str(k),'th round']) - chart_fileName = [datapath,'chart',num2str(k),'.mat']; - FPT_fileName = [datapath,'Halfmoon_FPT',num2str(k),'.mat']; + chart_fileName = [params.paths.datapath,'chart',num2str(k),'.mat']; + FPT_fileName = [params.paths.datapath,'Halfmoon_FPT',num2str(k),'.mat']; load( chart_fileName ) K = length(chart); phi_store = zeros(1,K); for i = 1:K X_int = chart{i}.X_int; - phi_store(i) = atan2(X_int(2), X_int(1)); + phi_store(i) = atan2(X_int(2), X_int(1)); end - [chart, connectivity, P] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - chart_sim_parameter.connectivity = connectivity; - chart_sim_parameter.X_int = chart{chart_sim_parameter.nearest}.X_int; + [chart, connectivity, P] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + params.chart_sim_parameter.connectivity = connectivity; + params.chart_sim_parameter.X_int = chart{params.chart_sim_parameter.nearest}.X_int; mode = 2; well_threshold = [-2*atan(4-sqrt(15)), -2*atan(4+sqrt(15))+2*pi]; disp('Fast mode projection') option = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_sim1, t_final(1,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, N_IC, phi_store, chart,chart_sim_parameter); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_sim1, t_final(1,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, params.MFPT.N_IC, phi_store, chart,params.chart_sim_parameter); disp('Orthogonal projection') option = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_ort1, t_final(2,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, N_IC, phi_store, chart,chart_sim_parameter); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_ort1, t_final(2,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, params.MFPT.N_IC, phi_store, chart,params.chart_sim_parameter); well_threshold =[-2*atan(4+sqrt(15)), -2*atan(4-sqrt(15))]; disp('Fast mode projection') option = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_sim2, t_final(3,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, N_IC, phi_store, chart,chart_sim_parameter); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_sim2, t_final(3,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, params.MFPT.N_IC, phi_store, chart,params.chart_sim_parameter); disp('Orthogonal projection') option = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_ort2, t_final(4,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, N_IC, phi_store, chart, chart_sim_parameter); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_ort2, t_final(4,k)] = MFPT_halfmoon_sim(weighted_dd2, well_threshold, output, phi_ori, params.MFPT.N_IC, phi_store, chart, params.chart_sim_parameter); save(FPT_fileName,'FPT_ori1','FPT_ori2','FPT_sim1','FPT_sim2','FPT_ort1','FPT_ort2') @@ -66,14 +65,14 @@ if k ==1 disp(['Original simulator']) - disp(['The original simulator at right well has ',num2str(mean(FPT_ori1),'%5.0f'),'+-',num2str(std(FPT_ori1)*1.96/sqrt(N_IC),'%5.0f')]) - disp(['The original simulator at left well has ', num2str(mean(FPT_ori2),'%5.0f'),'+-',num2str(std(FPT_ori2)*1.96/sqrt(N_IC),'%5.0f')]) + disp(['The original simulator at right well has ',num2str(mean(FPT_ori1),'%5.0f'),'+-',num2str(std(FPT_ori1)*1.96/sqrt(params.MFPT.N_IC),'%5.0f')]) + disp(['The original simulator at left well has ', num2str(mean(FPT_ori2),'%5.0f'),'+-',num2str(std(FPT_ori2)*1.96/sqrt(params.MFPT.N_IC),'%5.0f')]) disp(['Oblique projection']) - disp(['The ATLAS simulator at right well has ',num2str(mean(FPT_sim1),'%5.0f'),'+-',num2str(std(FPT_sim1)*1.96/sqrt(N_IC),'%5.0f')]) - disp(['The ATLAS simulator at left well has ', num2str(mean(FPT_sim2),'%5.0f'),'+-',num2str(std(FPT_sim2)*1.96/sqrt(N_IC),'%5.0f')]) + disp(['The ATLAS simulator at right well has ',num2str(mean(FPT_sim1),'%5.0f'),'+-',num2str(std(FPT_sim1)*1.96/sqrt(params.MFPT.N_IC),'%5.0f')]) + disp(['The ATLAS simulator at left well has ', num2str(mean(FPT_sim2),'%5.0f'),'+-',num2str(std(FPT_sim2)*1.96/sqrt(params.MFPT.N_IC),'%5.0f')]) disp(['Orthogonal projection']) - disp(['The ATLAS simulator at right well has ',num2str(mean(FPT_ort1),'%5.0f'),'+-',num2str(std(FPT_ort1)*1.96/sqrt(N_IC),'%5.0f')]) - disp(['The ATLAS simulator at left well has ', num2str(mean(FPT_ort2),'%5.0f'),'+-',num2str(std(FPT_ort2)*1.96/sqrt(N_IC),'%5.0f')]) + disp(['The ATLAS simulator at right well has ',num2str(mean(FPT_ort1),'%5.0f'),'+-',num2str(std(FPT_ort1)*1.96/sqrt(params.MFPT.N_IC),'%5.0f')]) + disp(['The ATLAS simulator at left well has ', num2str(mean(FPT_ort2),'%5.0f'),'+-',num2str(std(FPT_ort2)*1.96/sqrt(params.MFPT.N_IC),'%5.0f')]) end @@ -93,4 +92,5 @@ disp(['Average relative error at Orthogonal ATLAS simulator at left state is ', num2str(mean(relative_error_FPT(4,:)),'%5.4f'),'+-',num2str(std(relative_error_FPT(4,:))*1.96/sqrt(10),'%5.4f')]) - save(Halfmoon_analysis_fileName,'Mean_FPT','relative_error_FPT','t_final_ori','t_final') \ No newline at end of file + Halfmoon_analysis_fileName = [params.paths.datapath,'Halfmoon_analysis.mat']; + save(Halfmoon_analysis_fileName,'Mean_FPT','relative_error_FPT','t_final_ori','t_final') diff --git a/model/Halfmoon/halfmoon_ori_simulation.m b/model/Halfmoon/halfmoon_ori_simulation.m index 5177796..755da53 100644 --- a/model/Halfmoon/halfmoon_ori_simulation.m +++ b/model/Halfmoon/halfmoon_ori_simulation.m @@ -1,15 +1,13 @@ - %% Original simulator - disp('Use original simulation to simulate one single traj.') Sim_parameter = struct( ... - 'T_max', T_one*3, ... - 'dt', dt , ... + 'T_max', params.learning.T_one*3, ... + 'dt', params.RHS.dt , ... 'gap', 1, ... - 'X_int', X_int ... + 'X_int', params.learning.X_int ... ); - [output] = simulator_one(Sim_parameter); + [output] = params.simulator.one(Sim_parameter); [L_out,~] = size(output); - temp = nonlin_trans_inv(output'); + temp = params.RHS.nonlin_trans_inv(output'); phi_ori = temp(1,:); disp(['one single traj is simulated with the original simulator']) diff --git a/model/Halfmoon/halfmoon_relearning.m b/model/Halfmoon/halfmoon_relearning.m index 66267ac..80b7df4 100644 --- a/model/Halfmoon/halfmoon_relearning.m +++ b/model/Halfmoon/halfmoon_relearning.m @@ -1,18 +1,18 @@ - RHS_parameter_relearn.modify = 1; + params.relearn.RHS.modify = 1; relearn_parameter = struct( ... - 'N_relearn', 2*N, ... + 'N_relearn', 2*params.RHS.N, ... 'iter', 1, ... - 'RHS_parameter', RHS_parameter_relearn, ... + 'RHS_parameter', params.relearn.RHS, ... 'relative_threshold', [0.005, 0.02], ... - 'simulator_par', @(Sim_parameter) simEulernonlintrans_par(RHS_parameter_relearn,Sim_parameter) ... - ); + 'simulator_par', @(Sim_parameter) simEulernonlintrans_par(params.relearn.RHS,Sim_parameter) ... + ); disp('Starting Relearning Chart') - K = length(chart); + K = length(chart); index_to_learn = 1:K; [ chart,~ ] = relearn_chart(chart, relearn_parameter, index_to_learn); - [chart, connectivity, P] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - save( chart_fileName ,'chart','connectivity','P') - + [chart, connectivity, P] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + save( params.paths.chart ,'chart','connectivity','P') + disp('Relearning stage is completed.') diff --git a/model/Halfmoon/halfmoon_simulation_learning.m b/model/Halfmoon/halfmoon_simulation_learning.m index 13a17e2..76cc575 100644 --- a/model/Halfmoon/halfmoon_simulation_learning.m +++ b/model/Halfmoon/halfmoon_simulation_learning.m @@ -1,24 +1,22 @@ - disp('Use ATLAS simulator to run long trajs') - for i = 1:200 disp(['traj No.',num2str(i) ]) K = length(chart); - chart_sim_parameter.nearest = max(round(rand * K),1); - chart_sim_parameter.connectivity = connectivity; - chart_sim_parameter.X_int = chart{chart_sim_parameter.nearest}.X_int; - chart_sim_parameter.Nstep = 1*10^5; - [~, ~, chart] = ATLAS_simulator2(weighted_dd2, chart_sim_parameter,RHS_parameter, simulator_par, chart); - [chart, connectivity, P] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - if relearn_option == 1 - K_next = length(chart); + params.chart_sim_parameter.nearest = max(round(rand * K),1); + params.chart_sim_parameter.connectivity = connectivity; + params.chart_sim_parameter.X_int = chart{params.chart_sim_parameter.nearest}.X_int; + params.chart_sim_parameter.Nstep = 1*10^5; + [~, ~, chart] = ATLAS_simulator2(weighted_dd2, params.chart_sim_parameter,params.RHS, params.simulator.parallel, chart); + [chart, connectivity, P] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + if params.relearn.option == 1 + K_next = length(chart); if K_next>K index_to_learn = K+1: K_next; - [ chart ] = relearn_chart(chart, relearn_parameter, index_to_learn); - [chart, connectivity, P, bin_N ] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - save( chart_fileName ,'chart','connectivity','P') + [ chart ] = relearn_chart(chart, params.relearn.settings, index_to_learn); + [chart, connectivity, P, bin_N ] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + save( params.paths.chart ,'chart','connectivity','P') end end - - save( chart_fileName ,'chart','connectivity','P') - end - disp('Simulation stage is completed') \ No newline at end of file + + save( params.paths.chart ,'chart','connectivity','P') + end + disp('Simulation stage is completed') diff --git a/model/Halfmoon/set_parameter.m b/model/Halfmoon/set_parameter.m index 07838c3..d2bea99 100644 --- a/model/Halfmoon/set_parameter.m +++ b/model/Halfmoon/set_parameter.m @@ -1,118 +1,118 @@ +function params = set_parameter() +%SET_PARAMETER Configure the Halfmoon benchmark parameters and metadata. +% +% params = SET_PARAMETER() returns a struct that exposes the configuration +% formerly initialised in the base workspace. Downstream scripts should +% consume the returned struct instead of depending on implicitly defined +% variables. + + datapath = evalin('base','datapath'); + D = 20; % The dimension of the full system - d = 1; % The dimension of the slow manifold. Later we can estimate this dimenstion - chi_p = 1.96; % for 95% Confidence interval for 1D Gaussian distribution - K_int = 10; % Number of initial points on the manifold. - N = 8000000; %6000000; % Number of short traj starting from the same initial point - dt = 0.05; % Simulation time-step - T_one = 800000; % Time for single long trajectory to learn ini K charts. - - t0 = 20*dt; %10*dt; % The time as the benchmark. Drift and diffusion varies small within t0 - T_max = 25*dt; % Time for short trajctories. + d = 1; % The dimension of the slow manifold. Later we can estimate this dimension + chi_p = 1.96; % for 95% Confidence interval for 1D Gaussian distribution + K_int = 10; % Number of initial points on the manifold. + N = 8000000; % Number of short trajectories starting from the same initial point + dt = 0.05; % Simulation time-step + T_one = 800000; % Time for single long trajectory to learn initial K charts. + + t0 = 20*dt; % The time as the benchmark. Drift and diffusion varies small within t0 + T_max = 25*dt; % Time for short trajectories. X_int = ones(1,D); X_int(2) = 0; - threshold = [0.5, 0.01]; % The first one is R_max determined by the diameter of the slow manifold when setting the value of R and the second one threshold for occational prj - connectivity_threshold = 4; % two landmarks are connected if they are with in r*sqrt(t0) + threshold = [0.5, 0.01]; % Thresholds for projection and occasional projection + connectivity_threshold = 4; % Two landmarks are connected if they are within r*sqrt(t0) explore_threshold = 0.95; - option = 1; % 1 for fast mode prj, 2 for orthogonal prj - modify = 0; % 0 use the intercept(t=0) 1 use the tau_min point as the landmark + option = 1; % 1 for fast mode projection, 2 for orthogonal projection + modify = 0; % 0 use the intercept(t=0) 1 use the tau_min point as the landmark relearn_option = 1; - - LowerBound = 20*dt; % Lowerbound of the training window - UpperBound = 25*dt; % Upperbound of the training window - -Nstep = 1*10^5; -gap=1; -a1 = 0; -a2 = 5*10^-3; -a3 = 2.5*10^-3; -a4 = 0.06; -eps = 0.01; -b1 = 0.04/eps; -b2 = 0.035/sqrt(eps); -b3 = 0.05/eps; -b4 = 0.02/sqrt(eps); -t = -1; - -Tran = diag(ones(1,D-2),1); -Tran(:,1)=ones(D-1,1); -Tran=Tran(1:D-2,:); - -Tran_inv = diag(ones(1,D-2),1); -Tran_inv(:,1)=ones(D-1,1)*(-1); -Tran_inv=Tran_inv(1:D-2,:); - - -parameter = struct( ... + LowerBound = 20*dt; % Lower bound of the training window + UpperBound = 25*dt; % Upper bound of the training window + + Nstep = 1*10^5; + gap = 1; + a1 = 0; + a2 = 5*10^-3; + a3 = 2.5*10^-3; + a4 = 0.06; + eps = 0.01; + b1 = 0.04/eps; + b2 = 0.035/sqrt(eps); + b3 = 0.05/eps; + b4 = 0.02/sqrt(eps); + t = -1; + + Tran = diag(ones(1,D-2),1); + Tran(:,1)=ones(D-1,1); + Tran=Tran(1:D-2,:); + + Tran_inv = diag(ones(1,D-2),1); + Tran_inv(:,1)=ones(D-1,1)*(-1); + Tran_inv=Tran_inv(1:D-2,:); + + parameter = struct( ... 'a1', a1, ... - 'a2', a2, ... + 'a2', a2, ... 'a3', a3, ... 'a4', a4, ... - 'b1', b1, ... - 'b2', b2, ... - 'b3', b3, ... - 'b4', b4, ... - 'Tran', Tran, ... - 'Tran_inv', Tran_inv, ... - 't', t ... + 'b1', b1, ... + 'b2', b2, ... + 'b3', b3, ... + 'b4', b4, ... + 'Tran', Tran, ... + 'Tran_inv', Tran_inv, ... + 't', t ... ); - -drift = @(X) get_drift(X, parameter); -diffusion = @(X) get_diffusion(X, parameter); -nonlin_trans = @(X) get_nonlin_trans(X, parameter); -nonlin_trans_inv = @(X) get_nonlin_trans_inv(X, parameter); -%weighted_dd = @(X_curr, chart,neigh) weighted_drift_diffusion(X_curr, chart, t0, chi_p, neigh, D,d, threshold, option); - -RHS_parameter = struct( ... - 'T_max', T_max,... - 'dt', dt, ... - 'N', N, ... - 't0', t0, ... - 'chi_p', chi_p, ... - 'threshold', threshold, ... + + drift = @(X) get_drift(X, parameter); + diffusion = @(X) get_diffusion(X, parameter); + nonlin_trans = @(X) get_nonlin_trans(X, parameter); + nonlin_trans_inv = @(X) get_nonlin_trans_inv(X, parameter); + + RHS_parameter = struct( ... + 'T_max', T_max, ... + 'dt', dt, ... + 'N', N, ... + 't0', t0, ... + 'chi_p', chi_p, ... + 'threshold', threshold, ... 'connectivity_threshold',connectivity_threshold, ... 'D', D, ... 'd', d, ... 'modify', modify, ... + 'parameter', parameter, ... 'drift', drift, ... 'diffusion', diffusion, ... 'nonlin_trans', nonlin_trans, ... 'nonlin_trans_inv', nonlin_trans_inv,... 'UpperBound', UpperBound,... 'LowerBound', LowerBound ... -); - - - + ); -simulator_one = @(Sim_parameter) simEulernonlintrans_one_traj(RHS_parameter,Sim_parameter); -simulator_no_par = @(Sim_parameter) simEulernonlintrans(RHS_parameter,Sim_parameter); -simulator_par = @(Sim_parameter) simEulernonlintrans_par(RHS_parameter,Sim_parameter); + simulator_one = @(Sim_parameter) simEulernonlintrans_one_traj(RHS_parameter,Sim_parameter); + simulator_no_par = @(Sim_parameter) simEulernonlintrans(RHS_parameter,Sim_parameter); + simulator_par = @(Sim_parameter) simEulernonlintrans_par(RHS_parameter,Sim_parameter); - - MSM_parameter = struct( ... + MSM_parameter = struct( ... 'step', 1, ... 'N_state', 1*10^6, ... 'dt_s', t0 ... ); - - - RHS_parameter_relearn = RHS_parameter; % one can change some RHS if you believe the invariant manifold is not changed. - relearn_parameter = struct( ... + RHS_parameter_relearn = RHS_parameter; % one can change some RHS if you believe the invariant manifold is not changed. + relearn_parameter = struct( ... 'N_relearn', N, ... - 'iter', 1, ... - 'RHS_parameter', RHS_parameter_relearn, ... - 'relative_threshold', [0.05, 0.1], ... + 'iter', 1, ... + 'RHS_parameter', RHS_parameter_relearn, ... + 'relative_threshold', [0.05, 0.1], ... 'simulator_par', @(Sim_parameter) simEulernonlintrans_par(RHS_parameter_relearn,Sim_parameter) ... ); - - -chart_sim_parameter = struct ( ... + chart_sim_parameter = struct ( ... 'X_int', [], ... 'nearest', 4, ... - 't', t, ... + 't', t, ... 't0', t0, ... 'dt_s', t0, ... 'D', RHS_parameter.D, ... @@ -122,15 +122,48 @@ 'connectivity', [], ... 'explore_threshold', explore_threshold, ... 'N', N ... - ); - - % parameter for MFPTlength - - N_IC = 30000; - - -chart_fileName = [datapath,'chart.mat']; -chart_part_fileName = [datapath,'chart_part.mat']; -TranM_fileName = [datapath,'TranM.mat']; -well_fileName = [datapath,'well.mat']; -FPT_fileName = [datapath,'Halfmoon_FPT.mat']; \ No newline at end of file + ); + + N_IC = 30000; % Parameter for MFPT length + + chart_fileName = [datapath,'chart.mat']; + chart_part_fileName = [datapath,'chart_part.mat']; + TranM_fileName = [datapath,'TranM.mat']; + well_fileName = [datapath,'well.mat']; + FPT_fileName = [datapath,'Halfmoon_FPT.mat']; + + params = struct(); + params.RHS = RHS_parameter; + params.simulator = struct( ... + 'one', simulator_one, ... + 'serial', simulator_no_par, ... + 'parallel', simulator_par ... + ); + params.MSM = MSM_parameter; + params.chart_sim_parameter = chart_sim_parameter; + params.relearn = struct( ... + 'option', relearn_option, ... + 'RHS', RHS_parameter_relearn, ... + 'settings', relearn_parameter ... + ); + params.learning = struct( ... + 'K_int', K_int, ... + 'T_one', T_one, ... + 'X_int', X_int, ... + 'explore_threshold', explore_threshold, ... + 'option', option ... + ); + params.simulation = struct( ... + 'Nstep', Nstep, ... + 'gap', gap ... + ); + params.paths = struct( ... + 'datapath', datapath, ... + 'chart', chart_fileName, ... + 'chart_part', chart_part_fileName, ... + 'TranM', TranM_fileName, ... + 'well', well_fileName, ... + 'FPT', FPT_fileName ... + ); + params.MFPT = struct('N_IC', N_IC); +end diff --git a/model/Peanut/peanut_MFPT.m b/model/Peanut/peanut_MFPT.m index 8e6b0ee..dc581f9 100644 --- a/model/Peanut/peanut_MFPT.m +++ b/model/Peanut/peanut_MFPT.m @@ -3,85 +3,82 @@ relative_error_FPT = zeros(8,10); for k = 1:10 - disp([num2str(k),'th round']) - chart_fileName = [datapath,'chart',num2str(k),'.mat']; - TranM_fileName = [datapath,'TranM',num2str(k),'.mat']; - FPT_fileName = [datapath,'Peanut_FPT',num2str(k),'.mat']; - + chart_fileName = [params.paths.datapath,'chart',num2str(k),'.mat']; + TranM_fileName = [params.paths.datapath,'TranM',num2str(k),'.mat']; + FPT_fileName = [params.paths.datapath,'Peanut_FPT',num2str(k),'.mat']; load(chart_fileName); load(TranM_fileName); mode = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - chart_sim_parameter.connectivity = connectivity; - chart_sim_parameter.X_int = chart{chart_sim_parameter.nearest}.X_int; - chart_sim_parameter.explore_threshold = 8; - chart_sim_parameter.gap = 1; - %chart_sim_parameter.Nstep = 2*10^7; - [X, nearest_store, chart] = ATLAS_simulator2(weighted_dd2, chart_sim_parameter,RHS_parameter, simulator_par, chart); disp(['One single trajectory is simulated of time ', num2str(chart_sim_parameter.Nstep*chart_sim_parameter.dt_s)]) + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0,chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, params.learning.option,mode ); + params.chart_sim_parameter.connectivity = connectivity; + params.chart_sim_parameter.X_int = chart{params.chart_sim_parameter.nearest}.X_int; + params.chart_sim_parameter.explore_threshold = 8; + params.chart_sim_parameter.gap = 1; + [X, nearest_store, chart] = ATLAS_simulator2(weighted_dd2, params.chart_sim_parameter, params.RHS, params.simulator.parallel, chart); disp(['One single trajectory is simulated of time ', num2str(params.chart_sim_parameter.Nstep*params.chart_sim_parameter.dt_s)]) %% Calculate MFPT - disp('Starting MFPT part') - N_IC = 15000; + disp('Starting MFPT part') + N_IC = params.MFPT.N_IC; set_well; %% disp('First Well: 0.05 level exit 0.02 level') - [FPT_ori1, t_final(1,k)] = MFPT_peanut_ori_general(well_threshold1, X_int_store1, RHS_parameter,chart_angle); + [FPT_ori1, t_final(1,k)] = MFPT_peanut_ori_general(well_threshold1, X_int_store1, params.RHS,chart_angle); disp('Fast mode projection') option = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_sim1, t_final(2,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold1, X_int_store1,nearest_store1, chart,chart_sim_parameter,chart_angle); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_sim1, t_final(2,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold1, X_int_store1,nearest_store1, chart,params.chart_sim_parameter,chart_angle); disp('Orthogonal projection') - option = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_ort1, t_final(3,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold1, X_int_store1,nearest_store1, chart,chart_sim_parameter,chart_angle); + option = 2; + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_ort1, t_final(3,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold1, X_int_store1,nearest_store1, chart,params.chart_sim_parameter,chart_angle); %% - disp('Second Well: -0.05 level exit -0.02 level') - [FPT_ori2, t_final(4,k)] = MFPT_peanut_ori_general(well_threshold2, X_int_store2,RHS_parameter,chart_angle); + disp('Second Well: -0.05 level exit -0.02 level') + [FPT_ori2, t_final(4,k)] = MFPT_peanut_ori_general(well_threshold2, X_int_store2,params.RHS,chart_angle); disp('Fast mode projection') option = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_sim2, t_final(5,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold2, X_int_store2,nearest_store2, chart,chart_sim_parameter,chart_angle); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_sim2, t_final(5,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold2, X_int_store2,nearest_store2, chart,params.chart_sim_parameter,chart_angle); disp('Orthogonal projection') - option = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_ort2, t_final(6,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold2, X_int_store2,nearest_store2, chart,chart_sim_parameter,chart_angle); + option = 2; + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_ort2, t_final(6,k)] = MFPT_peanut_sim_general(weighted_dd2, well_threshold2, X_int_store2,nearest_store2, chart,params.chart_sim_parameter,chart_angle); %% disp('Starting First Well: 0.05 level exit zero line') - [FPT_ori3, t_final(7,k)] = MFPT_peanut_ori_zeroline(X_int_store1,RHS_parameter,phi_zero, theta_zero); + [FPT_ori3, t_final(7,k)] = MFPT_peanut_ori_zeroline(X_int_store1,params.RHS,phi_zero, theta_zero); disp('Fast mode projection') option = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_sim3, t_final(8,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store1, nearest_store1, chart,chart_sim_parameter,phi_zero, theta_zero); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_sim3, t_final(8,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store1, nearest_store1, chart,params.chart_sim_parameter,phi_zero, theta_zero); option = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_ort3, t_final(9,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store1, nearest_store1, chart,chart_sim_parameter,phi_zero, theta_zero); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_ort3, t_final(9,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store1, nearest_store1, chart,params.chart_sim_parameter,phi_zero, theta_zero); %% disp('Starting Second Well: -0.05 level exit zero line') - [FPT_ori4, t_final(10,k)] = MFPT_peanut_ori_zeroline(X_int_store2,RHS_parameter,phi_zero, theta_zero); + [FPT_ori4, t_final(10,k)] = MFPT_peanut_ori_zeroline(X_int_store2,params.RHS,phi_zero, theta_zero); disp('Fast mode projection') option = 1; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_sim4, t_final(11,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store2, nearest_store2, chart,chart_sim_parameter,phi_zero, theta_zero); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_sim4, t_final(11,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store2, nearest_store2, chart,params.chart_sim_parameter,phi_zero, theta_zero); option = 2; - weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, t0, chi_p, D,d, threshold, option,mode ); - [FPT_ort4, t_final(12,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store2, nearest_store2, chart,chart_sim_parameter,phi_zero, theta_zero); + weighted_dd2 = @(X0, chart,neigh,nearest, connectivity_indices) weighted_drift_diffusion2( X0, chart, neigh, nearest, connectivity_indices, params.RHS.t0, params.RHS.chi_p, params.RHS.D,params.RHS.d, params.RHS.threshold, option,mode ); + [FPT_ort4, t_final(12,k)] = MFPT_peanut_sim_zeroline(weighted_dd2, X_int_store2, nearest_store2, chart,params.chart_sim_parameter,phi_zero, theta_zero); save(FPT_fileName,'FPT_ori1','FPT_ori2','FPT_ori3','FPT_ori4','FPT_sim1','FPT_sim2','FPT_sim3',... 'FPT_sim4','FPT_ort1','FPT_ort2','FPT_ort3','FPT_ort4') @@ -105,22 +102,22 @@ if k == 1 disp(['Original simulator']) - disp(['The original simulator at cyan state in case 1 has ', num2str(mean(FPT_ori1),'%5.1f'),'+-',num2str(std(FPT_ori1)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The original simulator at Red state in case 1 has ', num2str(mean(FPT_ori2),'%5.1f'),'+-',num2str(std(FPT_ori2)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The original simulator at cyan state in case 2 has ', num2str(mean(FPT_ori3),'%5.1f'),'+-',num2str(std(FPT_ori3)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The original simulator at Red state in case 2 has ', num2str(mean(FPT_ori4),'%5.1f'),'+-',num2str(std(FPT_ori4)*1.96/sqrt(N_IC),'%5.1f')]) + disp(['The original simulator at cyan state in case 1 has ', num2str(mean(FPT_ori1),'%5.1f'),'+-',num2str(std(FPT_ori1)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The original simulator at Red state in case 1 has ', num2str(mean(FPT_ori2),'%5.1f'),'+-',num2str(std(FPT_ori2)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The original simulator at cyan state in case 2 has ', num2str(mean(FPT_ori3),'%5.1f'),'+-',num2str(std(FPT_ori3)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The original simulator at Red state in case 2 has ', num2str(mean(FPT_ori4),'%5.1f'),'+-',num2str(std(FPT_ori4)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) disp(['Oblique projection']) - disp(['The ATLAS simulator at cyan state in case 1 has ', num2str(mean(FPT_sim1),'%5.1f'),'+-',num2str(std(FPT_sim1)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The ATLAS simulator at Red state in case 1 has ', num2str(mean(FPT_sim2),'%5.1f'),'+-',num2str(std(FPT_sim2)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The ATLAS simulator at cyan state in case 2 has ', num2str(mean(FPT_sim3),'%5.1f'),'+-',num2str(std(FPT_sim3)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The ATLAS simulator at Red state in case 2 has ', num2str(mean(FPT_sim4),'%5.1f'),'+-',num2str(std(FPT_sim4)*1.96/sqrt(N_IC),'%5.1f')]) + disp(['The ATLAS simulator at cyan state in case 1 has ', num2str(mean(FPT_sim1),'%5.1f'),'+-',num2str(std(FPT_sim1)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The ATLAS simulator at Red state in case 1 has ', num2str(mean(FPT_sim2),'%5.1f'),'+-',num2str(std(FPT_sim2)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The ATLAS simulator at cyan state in case 2 has ', num2str(mean(FPT_sim3),'%5.1f'),'+-',num2str(std(FPT_sim3)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The ATLAS simulator at Red state in case 2 has ', num2str(mean(FPT_sim4),'%5.1f'),'+-',num2str(std(FPT_sim4)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) disp(['Orthogonal Projection']) - disp(['The ATLAS simulator at cyan state in case 1 has ', num2str(mean(FPT_ort1),'%5.1f'),'+-',num2str(std(FPT_ort1)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The ATLAS simulator at Red state in case 1 has ', num2str(mean(FPT_ort2),'%5.1f'),'+-',num2str(std(FPT_ort2)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The ATLAS simulator at cyan state in case 2 has ', num2str(mean(FPT_ort3),'%5.1f'),'+-',num2str(std(FPT_ort3)*1.96/sqrt(N_IC),'%5.1f')]) - disp(['The ATLAS simulator at Red state in case 2 has ', num2str(mean(FPT_ort4),'%5.1f'),'+-',num2str(std(FPT_ort4)*1.96/sqrt(N_IC),'%5.1f')]) + disp(['The ATLAS simulator at cyan state in case 1 has ', num2str(mean(FPT_ort1),'%5.1f'),'+-',num2str(std(FPT_ort1)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The ATLAS simulator at Red state in case 1 has ', num2str(mean(FPT_ort2),'%5.1f'),'+-',num2str(std(FPT_ort2)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The ATLAS simulator at cyan state in case 2 has ', num2str(mean(FPT_ort3),'%5.1f'),'+-',num2str(std(FPT_ort3)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) + disp(['The ATLAS simulator at Red state in case 2 has ', num2str(mean(FPT_ort4),'%5.1f'),'+-',num2str(std(FPT_ort4)*1.96/sqrt(params.MFPT.N_IC),'%5.1f')]) end @@ -132,7 +129,7 @@ relative_error_FPT(5, k) = (Mean_FPT(9,k)-Mean_FPT(1,k))/Mean_FPT(1,k); relative_error_FPT(6, k) = (Mean_FPT(10,k)-Mean_FPT(2,k))/Mean_FPT(2,k); relative_error_FPT(7, k) = (Mean_FPT(11,k)-Mean_FPT(3,k))/Mean_FPT(3,k); - relative_error_FPT(8, k) = (Mean_FPT(12,k)-Mean_FPT(4,k))/Mean_FPT(4,k); + relative_error_FPT(8, k) = (Mean_FPT(12,k)-Mean_FPT(4,k))/Mean_FPT(4,k); end @@ -152,4 +149,5 @@ disp(['Average relative error at Orthogonal ATLAS simulator at Red state in case 2 is ', num2str(mean(relative_error_FPT(8,:)),'%5.4f'),'+-',num2str(std(relative_error_FPT(8,:))*1.96/sqrt(10),'%5.4f')]) - save(Peanut_analysis_fileName,'Mean_FPT','relative_error_FPT','t_final') \ No newline at end of file + Peanut_analysis_fileName = [params.paths.datapath,'Peanut_analysis.mat']; + save(Peanut_analysis_fileName,'Mean_FPT','relative_error_FPT','t_final') diff --git a/model/Peanut/peanut_simulation_learning.m b/model/Peanut/peanut_simulation_learning.m index 5e08c01..2b5fcac 100644 --- a/model/Peanut/peanut_simulation_learning.m +++ b/model/Peanut/peanut_simulation_learning.m @@ -1,6 +1,6 @@ disp('Use ATLAS simulator to simulate a single traj.') - chart_sim_parameter.Nstep = 2*10^7; - chart_sim_parameter.gap = 1000; - [~, ~, chart] = ATLAS_simulator2(weighted_dd2, chart_sim_parameter, RHS_parameter, simulator_par, chart); - [chart, connectivity, P] = landmark(chart, t0, chi_p,threshold(1), connectivity_threshold); - save( chart_fileName ,'chart','connectivity','P') \ No newline at end of file + params.chart_sim_parameter.Nstep = 2*10^7; + params.chart_sim_parameter.gap = 1000; + [~, ~, chart] = ATLAS_simulator2(weighted_dd2, params.chart_sim_parameter, params.RHS, params.simulator.parallel, chart); + [chart, connectivity, P] = landmark(chart, params.RHS.t0, params.RHS.chi_p,params.RHS.threshold(1), params.RHS.connectivity_threshold); + save( params.paths.chart ,'chart','connectivity','P') diff --git a/model/Peanut/set_parameter.m b/model/Peanut/set_parameter.m index 976aa3b..92e09e7 100644 --- a/model/Peanut/set_parameter.m +++ b/model/Peanut/set_parameter.m @@ -1,110 +1,102 @@ +function params = set_parameter() +%SET_PARAMETER Configure the Peanut benchmark parameters and metadata. +% +% params = SET_PARAMETER() returns a struct exposing the configuration that +% used to be scattered in the base workspace. Scripts should reference the +% returned struct to access settings, simulators, and file paths. - D = 3; % The dimension of the full system - d = 2; % The dimension of the slow manifold. Later we can estimate this dimenstion - chi_p = 5.991; % 95% Confidence interval for 2D Gaussian distribution - K_int = 100; % Number of initial points on the manifold. - N = 400000; % Number of short traj starting from the same initial point - epsilon = 0.005; % the small parameter - dt = 5*10^-4; % Simulation time-step - T_max = 0.1; % Time for short trajctories. - t0 = 0.1; %0.075; % The relaxtion time. - T_one = 1000; % Time for single long trajectory to learn ini K charts. - threshold = [1, 0.001]; % This is R_max determined by the diameter of the slow manifold when setting the value of R or the maximum possible layor - connectivity_threshold = 3; % two landmarks are connected if they are with in r*sqrt(t0) - explore_threshold = 1.05; - option = 1; % 1 for fast mode prj, 2 for orthogonal prj - modify = 0; % 0 use the intercept(t=0) 1 use the tau_min point as the landmark - LowerBound = 0.05; % Lowerbound of the training window - UpperBound = 0.10; % Upperbound of the training window - relearn_option = 0; - - - - %parameter for simuation - Nstep = 5*10^6; - gap = 10; - - % Parameters of RHS equations - a1 = 4; - a2 = 8; - c1 = 2; - c2 = 0.5; - c3 = 0.05; - c4 = 0.4; - c5 = 0.05; - c6 = 0.4; %0.5 - - % Specify the manifold and gradient - R = @(theta, phi) sqrt(a1+a2*(cos(theta)).^2); - Rprime = @(theta, phi) -a2*cos(theta).*sin(theta)./sqrt(a1+a2*cos(theta).^2); - - - % set_up an initial point to start - theta_int = rand * pi; - phi_int = rand * 2*pi; - R_int = R(theta_int, phi_int); - X_int = [R_int*sin(theta_int)*cos(phi_int), R_int*sin(theta_int)*sin(phi_int), R_int*cos(theta_int) ]; - - - % Specify the drift and diffusion in sph coordinate - r_RHS = @(r, theta, phi) -(c1./epsilon)./r.*( r - R(theta,phi) ); - theta_RHS = @(r, theta, phi) c3* cos(3*theta)./(r.*sin(theta)) ; % - c3* sin(2*theta + pi/4)./r; - phi_RHS = @(r, theta, phi) c5* sin(phi + theta)./r ; %c5 * sin(phi)./r ; - sigma_r = @(r, theta, phi) c2./(sqrt(epsilon).*r); - sigma_theta = @(r, theta, phi) c4.*sin(theta)./r; - sigma_phi = @(r, theta, phi) c6./r; - - b_sph = @(r, theta, phi) [ r_RHS(r,theta,phi); - theta_RHS(r,theta,phi); - phi_RHS(r,theta,phi) ]; - - Sigma_sph = @(r, theta, phi) [ sigma_r(r,theta,phi); - sigma_theta(r,theta,phi); - sigma_phi(r,theta,phi) ]; - - - - % Pre-load parameters to save some time. - drift = @(X) get_drift( X,epsilon,a1,a2,c1,c2,c3,c4,c5,c6 ); - diffusion = @(X) get_diffusion( X,epsilon,c2,c4,c6 ); - - - - - % Specify the drift and diffusion term on the reduced systems. The drift - % term has two parts b_stra and b_ito. The latter is the ito term. - - b_stra = @(theta, phi) [( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*cos(phi); - ( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*sin(phi); - ( Rprime(theta, phi).*cos(theta)-R(theta, phi).*sin(theta) ) ] .* theta_RHS( R(theta,phi), theta, phi ) + ... - [ -R(theta, phi).*sin(theta).*sin(phi); - R(theta, phi).*sin(theta).*cos(phi); - 0 ] .* phi_RHS( R(theta,phi), theta, phi ); - - H_true = @(theta, phi) [ ( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*cos(phi).*sigma_theta( R(theta,phi), theta, phi ), ... - -R(theta, phi).*sin(theta).*sin(phi).*sigma_phi( R(theta,phi), theta, phi ); - ( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*sin(phi).*sigma_theta( R(theta,phi), theta, phi ) , ... - R(theta, phi).*sin(theta).*cos(phi).*sigma_phi( R(theta,phi), theta, phi ); - ( Rprime(theta, phi).*cos(theta)-R(theta, phi).*sin(theta) ).*sigma_theta( R(theta,phi), theta, phi ), ... - 0 ] ; - - % this part is calculated by mathematica - b_ito = @(theta, phi) [ (1/2).*(a1+a2.*cos(theta).^2).^(-5/2).*cos(phi).*sin(theta).*((-1).*c6.^2.*( ... - a1+a2.*cos(theta).^2).^2+(-1).*c4.^2.*(a1+a2.*cos(theta).^2).*(a1+4.*a2.* ... - cos(theta).^2).*sin(theta).^2+a1.*a2.*c4.^2.*sin(theta).^4); - (-1/2).*(a1+a2.*cos(theta).^2).^(-5/2).*sin(theta).*(c6.^2.*(a1+a2.*cos(theta) ... - .^2).^2+c4.^2.*(a1+a2.*cos(theta).^2).*(a1+4.*a2.*cos(theta).^2).*sin(theta) ... - .^2+(-1).*a1.*a2.*c4.^2.*sin(theta).^4).*sin(phi); - (-1/4).*c4.^2.*cos(theta).*(a1+a2.*cos(theta).^2).^(-5/2).*(2.*a1.^2+ ... - a2.^2+2.*a2.*(3.*a1+a2).*cos(2.*theta)+a2.^2.*cos(4.*theta)).*sin(theta).^2; ]; - - - b_true = @(X) btrue(X,b_ito, b_stra); - Lambda_true = @(X) Lambdatrue(X,H_true); - R_true = @(X) Rtrue(X, R); - H_true_X = @(X) H_trueX(X, H_true); - - RHS_parameter = struct( ... + datapath = evalin('base','datapath'); + + D = 3; % The dimension of the full system + d = 2; % The dimension of the slow manifold + chi_p = 5.991; % 95% Confidence interval for 2D Gaussian distribution + K_int = 100; % Number of initial points on the manifold. + N = 400000; % Number of short trajectories starting from the same initial point + epsilon = 0.005; % Small parameter + dt = 5*10^-4; % Simulation time-step + T_max = 0.1; % Time for short trajectories. + t0 = 0.1; % The relaxation time. + T_one = 1000; % Time for single long trajectory to learn initial K charts. + threshold = [1, 0.001]; % Thresholds for projection and occasional projection + connectivity_threshold = 3; % Two landmarks are connected if within r*sqrt(t0) + explore_threshold = 1.05; + option = 1; % 1 for fast mode projection, 2 for orthogonal projection + modify = 0; % 0 use the intercept(t=0) 1 use the tau_min point as the landmark + LowerBound = 0.05; % Lower bound of the training window + UpperBound = 0.10; % Upper bound of the training window + relearn_option = 0; + + % Parameters for simulation + Nstep = 5*10^6; + gap = 10; + + % Parameters of RHS equations + a1 = 4; + a2 = 8; + c1 = 2; + c2 = 0.5; + c3 = 0.05; + c4 = 0.4; + c5 = 0.05; + c6 = 0.4; %0.5 + + % Specify the manifold and gradient + R = @(theta, phi) sqrt(a1+a2*(cos(theta)).^2); + Rprime = @(theta, phi) -a2*cos(theta).*sin(theta)./sqrt(a1+a2*cos(theta).^2); + + % Set up an initial point to start + theta_int = rand * pi; + phi_int = rand * 2*pi; + R_int = R(theta_int, phi_int); + X_int = [R_int*sin(theta_int)*cos(phi_int), R_int*sin(theta_int)*sin(phi_int), R_int*cos(theta_int) ]; + + % Specify the drift and diffusion in spherical coordinate + r_RHS = @(r, theta, phi) -(c1./epsilon)./r.*( r - R(theta,phi) ); + theta_RHS = @(r, theta, phi) c3* cos(3*theta)./(r.*sin(theta)) ; + phi_RHS = @(r, theta, phi) c5* sin(phi + theta)./r ; + sigma_r = @(r, theta, phi) c2./(sqrt(epsilon).*r); + sigma_theta = @(r, theta, phi) c4.*sin(theta)./r; + sigma_phi = @(r, theta, phi) c6./r; + + b_sph = @(r, theta, phi) [ r_RHS(r,theta,phi); + theta_RHS(r,theta,phi); + phi_RHS(r,theta,phi) ]; + + Sigma_sph = @(r, theta, phi) [ sigma_r(r,theta,phi); + sigma_theta(r,theta,phi); + sigma_phi(r,theta,phi) ]; + + % Pre-load parameters to save some time. + drift = @(X) get_drift( X,epsilon,a1,a2,c1,c2,c3,c4,c5,c6 ); + diffusion = @(X) get_diffusion( X,epsilon,c2,c4,c6 ); + + % Specify the drift and diffusion term on the reduced systems. The drift + % term has two parts b_stra and b_ito. The latter is the Ito term. + b_stra = @(theta, phi) [( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*cos(phi); + ( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*sin(phi); + ( Rprime(theta, phi).*cos(theta)-R(theta, phi).*sin(theta) ) ] .* theta_RHS( R(theta,phi), theta, phi ) + ... + [ -R(theta, phi).*sin(theta).*sin(phi); + R(theta, phi).*sin(theta).*cos(phi); + 0 ] .* phi_RHS( R(theta,phi), theta, phi ); + + H_true = @(theta, phi) [ ( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*cos(phi).*sigma_theta( R(theta,phi), theta, phi ), ... + -R(theta, phi).*sin(theta).*sin(phi).*sigma_phi( R(theta,phi), theta, phi ); + ( Rprime(theta, phi).*sin(theta)+R(theta, phi).*cos(theta) ).*sin(phi).*sigma_theta( R(theta,phi), theta, phi ) , ... + R(theta, phi).*sin(theta).*cos(phi).*sigma_phi( R(theta,phi), theta, phi ); + ( Rprime(theta, phi).*cos(theta)-R(theta, phi).*sin(theta) ).*sigma_theta( R(theta,phi), theta, phi ), ... + 0 ] ; + + % This part is calculated by Mathematica + b_ito = @(theta, phi) [ (1/2).*(a1+a2.*cos(theta).^2).^(-5/2).*cos(phi).*sin(theta).*((-1).*c6.^2.*(a1+a2.*cos(theta).^2).^2+(-1).*c4.^2.*(a1+a2.*cos(theta).^2).*(a1+4.*a2.*cos(theta).^2).*sin(theta).^2+a1.*a2.*c4.^2.*sin(theta).^4); + (-1/2).*(a1+a2.*cos(theta).^2).^(-5/2).*sin(theta).*(c6.^2.*(a1+a2.*cos(theta).^2).^2+c4.^2.*(a1+a2.*cos(theta).^2).*(a1+4.*a2.*cos(theta).^2).*sin(theta).^2+(-1).*a1.*a2.*c4.^2.*sin(theta).^4).*sin(phi); + (-1/4).*c4.^2.*cos(theta).*(a1+a2.*cos(theta).^2).^(-5/2).*(2.*a1.^2+a2.^2+2.*a2.*(3.*a1+a2).*cos(2.*theta)+a2.^2.*cos(4.*theta)).*sin(theta).^2; ]; + + b_true = @(X) btrue(X,b_ito, b_stra); + Lambda_true = @(X) Lambdatrue(X,H_true); + R_true = @(X) Rtrue(X, R); + H_true_X = @(X) H_trueX(X, H_true); + + RHS_parameter = struct( ... 'dt', dt , ... 'T_max', T_max , ... 'N', N, ... @@ -121,13 +113,12 @@ 'UpperBound', UpperBound, ... 'LowerBound', LowerBound ... ); - - - simulator_one = @(Sim_parameter) simEuler_one_traj(RHS_parameter,Sim_parameter); - simulator_no_par = @(Sim_parameter) simEuler(RHS_parameter,Sim_parameter); - simulator_par = @(Sim_parameter) simEuler_par(RHS_parameter,Sim_parameter); - - Exact_parameter = struct( ... + + simulator_one = @(Sim_parameter) simEuler_one_traj(RHS_parameter,Sim_parameter); + simulator_no_par = @(Sim_parameter) simEuler(RHS_parameter,Sim_parameter); + simulator_par = @(Sim_parameter) simEuler_par(RHS_parameter,Sim_parameter); + + Exact_parameter = struct( ... 'H_true_X', H_true_X, ... 'd', d, ... 'D', D, ... @@ -136,51 +127,81 @@ 'H_true', H_true, ... 'Lambda_true', Lambda_true ... ); - - - chart_sim_parameter = struct ( ... + + chart_sim_parameter = struct ( ... 'X_int', [], ... - 'nearest', 4, ... - 't0', t0, ... - 'dt_s', t0, ... + 'nearest', 4, ... + 't0', t0, ... + 'dt_s', t0, ... 'D', RHS_parameter.D, ... 'd', RHS_parameter.d, ... - 'Nstep', Nstep, ... + 'Nstep', Nstep, ... 'gap', gap, ... 'connectivity', [], ... 'explore_threshold', explore_threshold, ... 'N', N ... - ); - - MSM_parameter = struct( ... + ); + + MSM_parameter = struct( ... 'step', 1, ... 'N_state', 1*10^6, ... 'dt_s', t0 ... ); - - %parameter for relearning - - RHS_parameter_relearn = RHS_parameter; % one can change some RHS if you believe the invariant manifold is not changed. - relearn_parameter = struct( ... + + RHS_parameter_relearn = RHS_parameter; % one can change some RHS if you believe the invariant manifold is not changed. + relearn_parameter = struct( ... 'N_relearn', 10*N, ... 'iter', 2, ... 'RHS_parameter', RHS_parameter_relearn, ... 'relative_threshold', [0.005, 0.02], ... 'simulator_par', @(Sim_parameter) simHeun_par(RHS_parameter_relearn,Sim_parameter) ... ); - - % parameter for MFPT - N_IC = 15000; - - - -% file path -chart_fileName = [datapath,'chart.mat']; -chart_part_fileName = [datapath,'chart_part.mat']; -TranM_fileName = [datapath,'TranM.mat']; -FPT_fileName = [datapath,'Peanut_FPT.mat']; - - + N_IC = 15000; % Parameter for MFPT + chart_fileName = [datapath,'chart.mat']; + chart_part_fileName = [datapath,'chart_part.mat']; + TranM_fileName = [datapath,'TranM.mat']; + FPT_fileName = [datapath,'Peanut_FPT.mat']; + params = struct(); + params.RHS = RHS_parameter; + params.simulator = struct( ... + 'one', simulator_one, ... + 'serial', simulator_no_par, ... + 'parallel', simulator_par ... + ); + params.MSM = MSM_parameter; + params.chart_sim_parameter = chart_sim_parameter; + params.relearn = struct( ... + 'option', relearn_option, ... + 'RHS', RHS_parameter_relearn, ... + 'settings', relearn_parameter ... + ); + params.learning = struct( ... + 'K_int', K_int, ... + 'T_one', T_one, ... + 'X_int', X_int, ... + 'explore_threshold', explore_threshold, ... + 'option', option ... + ); + params.simulation = struct( ... + 'Nstep', Nstep, ... + 'gap', gap ... + ); + params.paths = struct( ... + 'datapath', datapath, ... + 'chart', chart_fileName, ... + 'chart_part', chart_part_fileName, ... + 'TranM', TranM_fileName, ... + 'FPT', FPT_fileName ... + ); + params.MFPT = struct('N_IC', N_IC); + params.exact = Exact_parameter; + params.manifold = struct( ... + 'R', R, ... + 'Rprime', Rprime, ... + 'b_sph', b_sph, ... + 'Sigma_sph', Sigma_sph ... + ); +end diff --git a/model/Peanut/set_well.m b/model/Peanut/set_well.m index 83f5e94..fa276a0 100644 --- a/model/Peanut/set_well.m +++ b/model/Peanut/set_well.m @@ -1,7 +1,7 @@ % store chart_angle K = length(chart); -chart_angle = zeros(d, K); +chart_angle = zeros(params.RHS.d, K); for i = 1:K X_int = chart{i}.X_int'; theta_X_int = atan2( sqrt( X_int(1)^2 + X_int(2)^2 ), X_int(3) ); @@ -70,10 +70,10 @@ index1 = find(ismember(nearest_store, well1)); L_IC = length(index1); - index1 = index1(round(linspace(1, L_IC, N_IC))); + index1 = index1(round(linspace(1, L_IC, params.MFPT.N_IC))); index2 = find(ismember(nearest_store, well2)); L_IC = length(index2); - index2 = index2(round(linspace(1, L_IC, N_IC))); + index2 = index2(round(linspace(1, L_IC, params.MFPT.N_IC))); X_int_store1 = X(index1,:); nearest_store1 = nearest_store(index1); X_int_store2 = X(index2,:); @@ -84,7 +84,7 @@ well_threshold2 =find(V(:,2)< -0.020); -chart_angle = zeros(d, K); +chart_angle = zeros(params.RHS.d, K); for i = 1:K X_int = chart{i}.X_int'; theta_X_int = atan2( sqrt( X_int(1)^2 + X_int(2)^2 ), X_int(3) );