Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Adding a unified approach for sequentiual Gibbs resimulation - sippi_…
…prior_mps and sippi_prior_cholesky OK
  • Loading branch information
Thomas Mejer Hansen committed Feb 10, 2017
1 parent d2b7e23 commit ac8b209
Show file tree
Hide file tree
Showing 4 changed files with 47 additions and 66 deletions.
44 changes: 6 additions & 38 deletions sippi_get_resim_data.m
Expand Up @@ -4,48 +4,16 @@
%
% c_cond [n_cond,4]: col1: x, col2: y, col4: z, col4: d
%
% See also sippi_prior
% See also sippi_prior, sippi_sequential_gibbs_resim
%
function d_cond=sippi_get_resim_data(m_current,prior,ip);
if nargin<3
ip=1;
end

step=prior{ip}.seq_gibbs.step;
if prior{ip}.seq_gibbs.type==1;
if ~isfield(prior{ip}.seq_gibbs,'n_ite');
prior{ip}.seq_gibbs.n_ite=1;
end
n_ite=prior{ip}.seq_gibbs.n_ite;

if length(step)<2, step(2:3)=step(1);end
if length(step)<3, step(3)=step(2);end

ih=[];
used=prior{ip}.xx.*0+1;
for j=1:n_ite
%
pos(1)=min(prior{ip}.x)+rand(1)*(max(prior{ip}.x)-min(prior{ip}.x));
pos(2)=min(prior{ip}.y)+rand(1)*(max(prior{ip}.y)-min(prior{ip}.y));
pos(3)=min(prior{ip}.z)+rand(1)*(max(prior{ip}.z)-min(prior{ip}.z));
% get index of data to resimulate
i_resim=sippi_sequential_gibbs_resim(prior,ip);
% get index of hard conditional data
ih=setxor(1:(prod(prior{ip}.dim)),i_resim);

used(find( (abs(prior{ip}.xx-pos(1))<step(1)) & (abs(prior{ip}.yy-pos(2))<step(2)) ))=0;
end
ih=find(used);
%if j>1; ih=unique(ih);end
d_cond=[prior{ip}.xx(ih(:)) prior{ip}.yy(ih(:)) prior{ip}.zz(ih(:)) m_current{ip}(ih(:))];


elseif prior{ip}.seq_gibbs.type==2
% RANDOM SELECTION OF MODEL PARAMETERS FOR RESIM
N=prod(size(m_current{ip}));
n_resim=prior{1}.seq_gibbs.step(1);
if n_resim<=1;
% if n_resim is less than one
% n_resim defines the fraction of N to use
n_resim=ceil(n_resim.*N);
end
n_cond=N-n_resim;
ih=randomsample(N,n_cond);
d_cond=[prior{ip}.xx(ih(:)) prior{ip}.yy(ih(:)) prior{ip}.zz(ih(:)) m_current{ip}(ih(:))];
end
d_cond=[prior{ip}.xx(ih(:)) prior{ip}.yy(ih(:)) prior{ip}.zz(ih(:)) m_current{ip}(ih(:))];
36 changes: 10 additions & 26 deletions sippi_prior_cholesky.m
Expand Up @@ -28,7 +28,7 @@
% nd X nd size covariance matrix.
% (it is computed the first the sippi_prior_cholesky is called)
%
% See also: gaussian_simulation_cholesky
% See also: gaussian_simulation_cholesky, sippi_sequential_gibbs_resim
%
function [m_propose,prior]=sippi_prior_cholesky(prior,m_current,ip);

Expand All @@ -55,37 +55,21 @@

%% Sequential Gibbs?
if nargin>1

z_cur=prior{ip}.z_rand;
z_init=z_cur;
z_new=randn(size(z_cur));

if prior{ip}.seq_gibbs.type==1
disp(sprintf('%s : Box type resimulation not implemented for CHOL type prior',mfilename));
elseif prior{ip}.seq_gibbs.type==2
n_all=prod(size(z_cur));
if prior{ip}.seq_gibbs.step<=1
% use n_resim as a proportion of all random deviates
n_resim=prior{ip}.seq_gibbs.step.*n_all;
else
n_resim=prior{ip}.seq_gibbs.step;
end
n_resim=ceil(n_resim);
n_resim = min([n_resim n_all]);

ii=randomsample(n_all,n_resim);

ii = sippi_sequential_gibbs_resim(prior,ip);

z_cur(ii)=randn(size(z_cur(ii)));

end

i_resim = sippi_sequential_gibbs_resim(prior,ip);

z_cur(i_resim)=randn(size(z_cur(i_resim)));

%% linear combinartion of the perturbed paramaters
if (isfield(prior{ip},'gradual'))
if prior{ip}.gradual<1
if (isfield(prior{ip}.seq_gibbs,'gradual'))
if prior{ip}.seq_gibbs.gradual<1
if exist('gaussian_linear_combine','file')
i_perturbed=ii;%find((z_init-z_cur)~=0);
z_cur(i_perturbed) = gaussian_linear_combine(z_init(i_perturbed),z_cur(i_perturbed),prior{ip}.gradual,0);
i_perturbed=i_resim;%find((z_init-z_cur)~=0);
z_cur(i_perturbed) = gaussian_linear_combine(z_init(i_perturbed),z_cur(i_perturbed),prior{ip}.seq_gibbs.gradual,0);
end
end
end
Expand Down
2 changes: 1 addition & 1 deletion sippi_prior_mps.m
Expand Up @@ -26,7 +26,7 @@
% All properties for each algorithm are availale in the prior{ip}.MPS
% field
%
% See also: sippi_prior, ti, mps_cpp
% See also: sippi_prior, ti, mps_cpp, sippi_sequential_gibbs_resim

%% Sequential Gibbs sampling type 1 (box selection of pixels)
% prior{ip}.seq_gibbs.type=1;%
Expand Down
31 changes: 30 additions & 1 deletion sippi_sequential_gibbs_resim.m
Expand Up @@ -13,8 +13,37 @@

if prior{ip}.seq_gibbs.type==1;
% BOX TYPE


i_resim=[];
sippi_verbose(sprintf('%s: BOX type resim not implemented',mfilename))

step = prior{ip}.seq_gibbs.step;

if ~isfield(prior{ip}.seq_gibbs,'n_ite');
prior{ip}.seq_gibbs.n_ite=1;
end
n_ite=prior{ip}.seq_gibbs.n_ite;

if length(step)<2, step(2:3)=step(1);end
if length(step)<3, step(3)=step(2);end

ih=[];
used=prior{ip}.xx.*0+1;
for j=1:n_ite
%
pos(1)=min(prior{ip}.x)+rand(1)*(max(prior{ip}.x)-min(prior{ip}.x));
pos(2)=min(prior{ip}.y)+rand(1)*(max(prior{ip}.y)-min(prior{ip}.y));
pos(3)=min(prior{ip}.z)+rand(1)*(max(prior{ip}.z)-min(prior{ip}.z));

used(find( (abs(prior{ip}.xx-pos(1))<step(1)) & (abs(prior{ip}.yy-pos(2))<step(2)) ))=0;
end
i_resim=find(1-used);





%sippi_verbose(sprintf('%s: BOX type resim not implemented',mfilename))

elseif prior{ip}.seq_gibbs.type==2;
% RANDOM MODEL PARAMETERS
Expand Down

0 comments on commit ac8b209

Please sign in to comment.