From 76aa174562f39be876c1099c34245ca0a839c5c5 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Tue, 19 Jan 2021 14:30:55 +0100
Subject: [PATCH] FRAP KL div works except for scaling factor. Normalization?
---
prob_laplace.m | 46 ++++++++++++++++++++++++++++------------------
1 file changed, 28 insertions(+), 18 deletions(-)
diff --git a/prob_laplace.m b/prob_laplace.m
index 55325cf..bb79134 100644
--- a/prob_laplace.m
+++ b/prob_laplace.m
@@ -348,7 +348,7 @@ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = linspace(0, 5, 300);
T_mov = Ternary_model(0, 'FRAP', {-5, b(7/3, 10^-6), 0.5, e(7/3),...
- 0, 1, 10, 7, 0, 'Constituent'}, t, 0.2);
+ 0, 1, 10, 7, 0, 'Constituent', 0}, t, 0.2);
T_mov.solve_tern_frap();
norm_fac = 1/sum(diff(T_mov.x)...
.*(T_mov.sol(1, 1:end-1)+T_mov.sol(1, 2:end))/2);
@@ -400,31 +400,37 @@ end
%% %%%%%%%%%%%%%%%%%%%%%% FRAP JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params = {-5, b(7/3, 10^-6), 0.5, e(7/3), 0, 1, 10, 7, 0, 'Constituent', 0};
-t = linspace(0, 5, 51);
+t = linspace(0, 2, 101);
bp = 5;
direc = -2; % 1: jump from left to right, -1: jump from right to left.
%%
-[ls, p_i_i] = frap_distro(t, params, 2, bp);
-[~, p_o_o] = frap_distro(t, params, -2, bp);
-[~, p_i_o] = frap_distro(t, params, 1, bp);
-[~, p_o_i] = frap_distro(t, params, -1, bp);
+[ls, p_i_i] = frap_distro(t, params, 2, bp, 1);
+[~, p_o_o] = frap_distro(t, params, -2, bp, 1);
+[~, p_i_o] = frap_distro(t, params, 1, bp, 1);
+[~, p_o_i] = frap_distro(t, params, -1, bp, 1);
+%%
+p_i_i(isnan(p_i_i)) = 0;
+p_o_o(isnan(p_o_o)) = 0;
+p_i_o(isnan(p_i_o)) = 0;
+p_o_i(isnan(p_o_i)) = 0;
%%
figure; hold on;
plot(ls, sum(p_i_i, 2));
plot(ls, sum(p_o_o, 2));
plot(ls, sum(p_i_o, 2));
plot(ls, sum(p_o_i, 2));
+save('/data/biophys/hubatsch/MatlabWorkspaces/prob_laplace_190121.mat');
%% Save time points for comparison with python
csvwrite('lb_in_out.csv', ls)
csvwrite('prob_in_out.csv', p);
%%
-KL = zeros(1, 20);
-for i = 1:20
+KL = zeros(1, 100);
+for i = 1:68
i_t = 1+i;
j_t = 2+i;
integrand = @(p, q) abs(p.*log(p./q));
-integrate = @(x, y) sum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
+integrate = @(x, y) nansum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
n_i = @(p1, p2, p3, p4) abs(integrate(ls, p1) + integrate(ls, p2)+...
integrate(ls, p3) + integrate(ls, p4));
N_i = n_i(p_i_i(:, i_t), p_o_o(:, i_t), p_i_o(:, i_t), p_o_i(:, i_t));
@@ -435,18 +441,22 @@ KL(i) = integrate(ls, integrand(p_i_i(:, i_t)/N_i, p_i_i(:, j_t)/N_j))+...
integrate(ls, integrand(p_i_o(:, i_t)/N_i, p_o_i(:, j_t)/N_j))+...
integrate(ls, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
end
-plot(t(2:21), abs(KL));
-ylim([0, 1])
+plot(t(2:101), abs(KL)/0.02);
+% ylim([0, 1])
%%
figure; hold on;
-plot(ls, p_i_i(:, i_t)/N_i);
-plot(ls, p_o_o(:, i_t)/N_i);
-plot(ls, p_i_o(:, i_t)/N_i);
-plot(ls, p_o_i(:, i_t)/N_i);
+plot(ls, p_i_i(:, j_t)/N_j, 'b');
+plot(ls, p_o_o(:, j_t)/N_j, 'b');
+plot(ls, p_i_o(:, j_t)/N_j, 'b');
+plot(ls, p_o_i(:, j_t)/N_j, 'b');
+%%
+plot(T_mov.t, s_dot)
+hold on; xlim([0, 2]);
+plot(t(2:101), 0.1908/0.1885*0.61*abs(KL)/0.02);
%% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [ls, p] = frap_distro(t, params, direc, bp)
+function [ls, p] = frap_distro(t, params, direc, bp, ind_delta)
x0 = sort(5-abs(direc)/direc*(0:0.002:4.01));
% Run simulations for 'delta' IC across outside
@@ -463,12 +473,12 @@ end
T_mov = Ternary_model(0, 'FRAP', params, t, 0.2);
T_mov.solve_tern_frap();
ls = -abs(direc)/direc*(0.001:0.04:4);
-n_T=45;
+n_T=70;
p = nan(length(ls), n_T);
for j = 1:n_T
tic
parfor i = 1:length(ls)
- p(i, j) = int_prob(ls(i), F, x0, direc, j, bp, 5, T_mov);
+ p(i, j) = int_prob(ls(i), F, x0, direc, j, bp, ind_delta, T_mov);
end
toc
end
--
GitLab