Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Initial support for 2D mixed point simulation (mixsim_2D)
  • Loading branch information
cultpenguin committed May 16, 2017
1 parent 30bf1ed commit eae33a4
Show file tree
Hide file tree
Showing 5 changed files with 172 additions and 2 deletions.
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions examples/prior_tests/prior_reals_fftma.m
Expand Up @@ -4,8 +4,8 @@
%% Define an FFTMA type a priori model
im=1;
prior{im}.type='FFTMA';
prior{im}.x=[0:.1:10]; % X array
prior{im}.y=[0:.1:20]; % Y array
prior{im}.x=[0:.1:10];
prior{im}.y=[0:.1:20];
prior{im}.m0=10;
prior{im}.Va='1 Sph(10,90,.25)';
prior{im}.fftma_options.constant_C=0;
Expand Down
96 changes: 96 additions & 0 deletions examples/prior_tests/prior_reals_mixsim.m
@@ -0,0 +1,96 @@
% prior_reals_mps: Sampling MIXSIM type prior models
rng('default');rng(1);
clear all,close all;


ip=1;
prior{ip}.type='mixsim'; % MPS type
%prior{ip}.x=[0:.5:10]; % X array
%prior{ip}.y=[0:.5:20]; % Y array
prior{ip}.x=[1:1:20]; % X array
prior{ip}.y=[1:1:40]; % Y array

prior{ip}.x=[1:1:60]; % X array % X array (CURRENTLY ONLY FOR INTEGEGER X STARTING FROM 1, WITH DX=1)
prior{ip}.y=[1:1:20]; % Y array % Y array (CURRENTLY ONLY FOR INTEGEGER Y STARTING FROM 1, WITH DY=1)



%% set TI
% Strebelles channels
TI=channels;
TI=TI(3:3:100,1:1:10)+1;
% 1D logs
%clear TI;load Logs1D.mat;for mm=1:size(Logs1D,2);TI{mm}=Logs1D(1:end,mm);end

prior{ip}.TI=TI;





%%
%% UNCONDITIONAL
figure(10);clf
for i=1:5;
[m,prior]=sippi_prior(prior);
subplot(1,5,i);
imagesc(prior{1}.x,prior{1}.y,m{1});
axis image
%xlabel('X')
%ylabel('Y')
title(sprintf('t=%3.2g',prior{1}.time))
drawnow;
end
caxis([1 prior{1}.options.sV])
colormap(sippi_colormap(1));
colorbar_shift;
print_mul(sprintf('prior_reals_%ss',prior{1}.type));
s=suptitle(sprintf('Independent realizations from a %s type prior',prior{1}.type))
set(s,'interpreter','none')


%% SEQUENTIAL GIBBS

[m2,prior]=sippi_prior(prior,m);

%% SEQUENTIAL GIBBS
figure(11);clf
prior{1}.seq_gibbs.type=[2];
prior{1}.seq_gibbs.step=[.9];

%prior{1}.seq_gibbs.type=[1];
%prior{1}.seq_gibbs.step=[10];

for i=1:5;
[m,prior]=sippi_prior(prior,m);
figure(11);
subplot(1,5,i);
imagesc(prior{1}.x,prior{1}.y,m{1});
axis image
colormap(sippi_colormap(1));
caxis([1 prior{1}.options.sV])
axis tight
set(gca,'nextplot','replacechildren');
drawnow;
end
print_mul(sprintf('prior_reals_%s_%d',prior{1}.type,prior{1}.seq_gibbs.type));
s=suptitle(sprintf('Independent realizations from a %s type prior (%d)',prior{1}.type,prior{1}.seq_gibbs.type))
set(s,'interpreter','none')

return
%% PRIOR MOVIE
fclose all;
fname=[mfilename,'_target','.mp4'];
vidObj = VideoWriter(fname,'MPEG-4');
open(vidObj);
N=200;
for i=1:N;
[m,prior]=sippi_prior(prior,m);
sippi_plot_prior(prior,m);
colormap(sippi_colormap(1));
caxis([1 prior{1}.options.sV]);
title(sprintf('Realizations %03d/%03d',i,N));
currFrame = getframe;
writeVideo(vidObj,currFrame);
end
close(vidObj);
74 changes: 74 additions & 0 deletions sippi_prior_mixsim.m
@@ -0,0 +1,74 @@
% sippi_prior_mixsim : MIXsim prior for SIPPI (only 2D)
%
% Example:
% ip=1;
% prior{ip}.type='mixsim'; % MIXSIM type
% prior{ip}.x=[1:1:20]; % X array
% prior{ip}.y=[1:1:40]; % Y array
% TI=channels;
% TI=TI(3:3:100,1:1:10)+1;
% prior{1}.TI=TI; % must be integer values starting with 1! (no zero values)
% m=sippi_prior(prior);
%
% sippi_plot_prior_sample(prior);
%
% See also sippi_prior, mixsim_2D
%

%
% Please refer to the following paper when using this algorithm:
% Cordua, K.S., T.M. Hansen, M.L. Gulbrandsen, C. Barnes, and
% K. Mosegaard, 2016. Mixed-point geostatistical simulation: A
% combination of two- and multiple-point geostatistics. Geophysical
% Research Letters.
%
% K.S. Cordua and T.M. Hansen (2016)
%
% See also: sippi_sequential_gibbs_resim
%
function [m_propose,prior]=sippi_prior_mixsim(prior,m_current,ip);

if nargin<3;
ip=1;
end

if ~isfield(prior{ip},'init')
prior=sippi_prior_init(prior);
end

if ~isfield(prior{ip},'TI');
sippi_verbose(sprintf('%s: No TI field (training image) set!',mfilename))
m_propose=[];
return
end


prior{ip}.options.do_cond=0; % no conditional simulation
%% Sequential gibbs resampling
if nargin>1
prior{ip}.options.do_cond=1;
d_cond=sippi_get_resim_data(m_current,prior,ip);
% set hard data
if ~isempty(d_cond);
prior{ip}.options.data=[d_cond(:,4), d_cond(:,2), d_cond(:,1)];
end
end


%%
prior{ip}.options.null=[];
[m_propose{ip},prior{ip}.options]=mixsim_2D(prior{ip}.dim(2),prior{ip}.dim(1),1,prior{ip}.TI,prior{ip}.options);

debug=0;
if (debug==1)
if nargin>1
figure(1);clf;
subplot(1,3,1);imagesc(m_current{ip});axis image;
try
subplot(1,3,2);scatter(prior{ip}.options.data(:,3),prior{ip}.options.data(:,2),20,prior{ip}.options.data(:,1),'filled');axis image
set(gca,'ydir','revers')
end
subplot(1,3,3);imagesc(m_propose{ip});axis image;
drawnow;
end
end

0 comments on commit eae33a4

Please sign in to comment.