From 6d4c10a90e18103334d98bf41ab992f03a3fc353 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Wed, 18 Jul 2018 11:04:06 +0200 Subject: [PATCH] update --- .gitignore | 2 + ERT/R2/Matlat2R2.m | 15 +-- ERT/R2/readOutput.m | 4 +- ERT/data_generation.m | 2 +- ERT/script_ERT.m | 253 ++++++++++++++++++++++++-------------- HT/covar2W.m | 25 ++++ HT/data_generation.m | 10 +- HT/script_HT.m | 280 ++++++++++++++++++++++++++++++------------ 8 files changed, 409 insertions(+), 182 deletions(-) create mode 100644 HT/covar2W.m diff --git a/.gitignore b/.gitignore index 63f6e73..4790b9d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,8 @@ ERT/data_gen/* HT/data_gen/* +HT/GW/* *.asv +*.eps *result/ *.mat diff --git a/ERT/R2/Matlat2R2.m b/ERT/R2/Matlat2R2.m index 38d6db7..d9a4af1 100644 --- a/ERT/R2/Matlat2R2.m +++ b/ERT/R2/Matlat2R2.m @@ -62,13 +62,14 @@ d.value = d.rho(:)' ; elseif d.job_type == 1 % for inversion, only one average value is given... - d.elem_1 = [1 ];% 3461 3509 3557 3605 3653 3701 3749 3797]; - d.elem_2 = [(d.numnp_y-1)*(d.numnp_x-1)];% 3472 3520 3568 3616 3664 3712 3760 3808]; - d.value = [d.rho_avg ];% 10 10 10 10 10 10 10 10] ; +% d.elem_1 = [1 ];% 3461 3509 3557 3605 3653 3701 3749 3797]; +% d.elem_2 = [(d.numnp_y-1)*(d.numnp_x-1)];% 3472 3520 3568 3616 3664 3712 3760 3808]; +% d.value = [d.rho_avg ];% 10 10 10 10 10 10 10 10] ; -% d.elem_1 = 1:((d.numnp_y-1)*(d.numnp_x-1)); -% d.elem_2 = 1:((d.numnp_y-1)*(d.numnp_x-1)); -% d.value = d.rho(:)' ; +warning('check this') + d.elem_1 = 1:((d.numnp_y-1)*(d.numnp_x-1)); + d.elem_2 = 1:((d.numnp_y-1)*(d.numnp_x-1)); + d.value = d.rho(:)' ; end d.num_regions = numel(d.value); % number of resistivity regions @@ -78,7 +79,7 @@ %% INVERSE SOLUTION if d.job_type==1 % inverse solution - % d.inverse_type = 1; % Inverse type: 0-pseudo-Marquardt, 1-regularised solution with linear filter (usual mode), + d.inverse_type = 1; % Inverse type: 0-pseudo-Marquardt, 1-regularised solution with linear filter (usual mode), d.target_decrease = 0; % 2-regularised type with quadratic filter, 3-qualitative solution or 4-blocked linear regularised type if d.mesh_type==4 || d.mesh_type==5 % quadrilateral mesh diff --git a/ERT/R2/readOutput.m b/ERT/R2/readOutput.m index 78c234b..a46ea9b 100644 --- a/ERT/R2/readOutput.m +++ b/ERT/R2/readOutput.m @@ -73,8 +73,8 @@ end end end - - % Data Weight + + % Data Weight data=importdata([d.filepath dataset '_err.dat']); diagWd=zeros(n_obs,1); diagWd(sum(abs(output.J))~=0) = data.data(:,5); diff --git a/ERT/data_generation.m b/ERT/data_generation.m index c0aaabf..dfdfb6e 100644 --- a/ERT/data_generation.m +++ b/ERT/data_generation.m @@ -256,7 +256,7 @@ A2=[A{:}]; fprintf(fid,'%d\n',A2(1,1)); for u=2:size(A2,1) - fprintf(fid,'%d %d %d %d %d %f %f\n',A2(u,:)); + fprintf(fid,'%d %d %d %d %d %.10e %.5f\n',A2(u,:)); end fclose(fid); diff --git a/ERT/script_ERT.m b/ERT/script_ERT.m index 74eba55..d1f8375 100644 --- a/ERT/script_ERT.m +++ b/ERT/script_ERT.m @@ -122,6 +122,21 @@ Sec.dnorm(i) = Res(~gen.i.grid.inside)'*Sec.d(~gen.i.grid.inside); end +% U = zeros(numel(Sec.d), numel(Prim.d)); +% for i=1:numel(Sec.d) +% Res = zeros(numel(Sec.y),numel(Sec.x)); +% Res(i)=1; +% F = griddedInterpolant({Sec.y,Sec.x},Res,'linear'); +% res_t = F({Prim.y,Prim.x}); +% U(i,:) = res_t(:) ./sum(res_t(:)); +% end +% U(isnan(U))=0; +% +% +% G2 =Sigma.res*U; +% imagesc(reshape(G * Prim.d(:), numel(Sec.y), numel(Sec.x))) +% figure; imagesc(reshape(G2 * Prim.d(:), numel(Sec.y), numel(Sec.x))) +% figure; imagesc(reshape((G2*Prim.d(:)-G * Prim.d(:))./(G * Prim.d(:)), numel(Sec.y), numel(Sec.x))) % View Resolution matrix % OLD large grid: i=[1838 4525 8502]; th=.05; x_lim=[16 25; 35 65; 90 99]; y_lim=[-.12 7; -.12 19.57; -.12 8]; ccaxis=[-.02 .02; -.002 .002; -.01 .01]; @@ -195,15 +210,15 @@ W(ij,:) = CCainv * CCb(:,ij); S(ij) = Cz0 - W(ij,:) * CCb(:,ij); end -save(['result/' fieldname '_cond'],'W','S','Prim_pt','G','Nscore','Sec','Prim') +% save(['result/' fieldname '_cond'],'W','S','Prim_pt','G','Nscore','Sec','Prim') % load(['result/' fieldname '_cond'],'W','S','Prim_pt','G','Nscore','Sec','Prim') % save(['result/' fieldname '_cond_noHD'],'W','S','Prim_pt','G','Nscore','Sec','Prim') % load(['result/' fieldname '_cond_noHD'],'W','S','Prim_pt','G','Nscore','Sec','Prim') % Compute the Kriging map -zh = reshape( W * [Sec.d(:)-Sec.dnorm(:) ; Prim_pt.d], numel(Prim.y), numel(Prim.x));numel(Prim.y), numel(Prim.x); -zhtest = reshape( W * [Test_Sec_d(:) ; Prim_pt.d], numel(Prim.y), numel(Prim.x)); +zh = reshape( W * [Sec.d(:)-Sec.dnorm(:) ; Prim_pt.d], numel(Prim.y), numel(Prim.x)); +zhtest = reshape( W * [Test_Sec_d(:)-Sec.dnorm(:) ; Prim_pt.d], numel(Prim.y), numel(Prim.x)); figure(5); clf; colormap(viridis()) subplot(2,1,1);surf(Prim.x,Prim.y,zh,'EdgeColor','none','facecolor','flat'); caxis([-3 3]) @@ -212,6 +227,14 @@ view(2); axis equal tight; set(gca,'Ydir','reverse'); xlabel('x');ylabel('y'); colorbar('southoutside'); title('Kriging Estimate Error Variance') %export_fig -eps 'Krig' +figure(55); clf; colormap(viridis()) +subplot(3,1,1);surf(Prim.x,Prim.y,zh,'EdgeColor','none','facecolor','flat'); caxis([-3 3]) +view(2); axis equal tight; set(gca,'Ydir','reverse'); xlabel('x');ylabel('y'); title('Kriging Estimate') +subplot(3,1,2);surf(Prim.x,Prim.y,zhtest,'EdgeColor','none','facecolor','flat'); +view(2); axis equal tight; set(gca,'Ydir','reverse'); xlabel('x');ylabel('y'); title('Kriging Estimate ') +subplot(3,1,3);surf(Prim.x,Prim.y,zh-zhtest,'EdgeColor','none','facecolor','flat'); +view(2); axis equal tight; set(gca,'Ydir','reverse'); xlabel('x');ylabel('y'); colorbar('southoutside'); title('Error Kriging Estimate') + figure(51); clf; S(S<=eps)=eps; histogram( (Prim.d(:)-zh(:))./ sqrt(S(:)) ) @@ -226,14 +249,13 @@ %% Simulation of the Area-to-point Kriging rng('shuffle'); -parm.n_real = 500; +parm.n_real = 30; covar = kriginginitiaite(gen.covar); zcs=nan(numel(Prim.y), numel(Prim.x),parm.n_real); z=nan(numel(Sec.y), numel(Sec.x),parm.n_real); for i_real=1:parm.n_real - zs(:,:,i_real) = fftma_perso(covar, struct('x',Prim.x,'y',Prim.y)); -end - %zs = fftma_perso(gen.covar, grid_gen); + % zs(:,:,i_real) = fftma_perso(covar, struct('x',Prim.x,'y',Prim.y)); + zs = fftma_perso(gen.covar, grid_gen); zhs = W * [G * zs(:) ; zs(Prim_pt.id)./1.3]; r = zh(:) + (zs(:) - zhs(:)); zcs(:,:,i_real) = reshape( r, numel(Prim.y), numel(Prim.x)); @@ -292,17 +314,17 @@ figure(9);clf; subplot(2,1,1); hold on; -h1=plot(Prim.x(1:2:end),vario_x(:,1:2:end)','b','color',[.5 .5 .5]); -h2=plot(Prim.x,vario_prim_x,'-r','linewidth',2); -h3=plot(Prim.x,1-covar.g(Prim.x*covar.cx(1)),'--k','linewidth',2); -xlim([0 60]); xlabel('Lag-distance h_x ');ylabel('Variogram \gamma(h_x)') +h1=plot(Prim.x(1:2:end)-Prim.x(1),vario_x(:,1:2:end)','b','color',[.5 .5 .5]); +h2=plot(Prim.x-Prim.x(1),vario_prim_x,'-r','linewidth',2); +h3=plot(Prim.x-Prim.x(1),1-covar.g((Prim.x-Prim.x(1))*covar.cx(1)),'--k','linewidth',2); +xlim([0 60]); xlabel('Lag-distance h_x ');ylabel('Variogram \gamma(h_x)'); ylim([0 1.5]) legend([h1(1) h2 h3],'500 realizations','True field','Theorical Model') subplot(2,1,2); hold on; -h1=plot(Prim.y,vario_y','b','color',[.5 .5 .5]); -h2=plot(Prim.y,vario_prim_y,'-r','linewidth',2); -h3=plot(Prim.y,1-covar.g(Prim.y*covar.cx(4)),'--k','linewidth',2); +h1=plot(Prim.y(1:2:end)-Prim.y(1),vario_y(:,1:2:end)','b','color',[.5 .5 .5]); +h2=plot(Prim.y-Prim.y(1),vario_prim_y,'-r','linewidth',2); +h3=plot(Prim.y-Prim.y(1),1-covar.g((Prim.y-Prim.y(1))*covar.cx(4)),'--k','linewidth',2); xlim([0 6]); xlabel('Lag-distance h_y ');ylabel('Variogram \gamma(h_y)') -legend([h1(1) h2 h3],'500 realizations','True field','Theorical Model') +legend([h1(1) h2 h3],'500 realizations','True field','Theorical Model'); ylim([0 1.2]) % export_fig -eps 'Vario' @@ -330,7 +352,7 @@ % parm.n_real=12; fsim_pseudo=nan(numel(gen.f.output.pseudo),parm.n_real); fsim_resistance=nan(numel(gen.f.output.resistance),parm.n_real); -rho = 1000./Nscore.inverse(zs); +rho = 1000./Nscore.inverse(zcs); % fieldname= 'AsSumbitted_2018-06-27_10-07'; @@ -346,7 +368,7 @@ for i_real=1:parm.n_real f.filepath = ['data_gen/IO-file-' num2str(i_real) '/']; mkdir(f.filepath) - copyfile('R2/R2.exe',[f.filepath 'R2.exe']) + %copyfile('R2/R2.exe',[f.filepath 'R2.exe']) copyfile(['result/' fieldname '_IO-file/forward/electrodes.dat'],[f.filepath 'electrodes.dat']) copyfile(['result\' fieldname '_IO-file\forward\protocol.dat'],[f.filepath 'protocol.dat']) end @@ -374,22 +396,16 @@ figure(11); hold on; -histogram(WRMSE); histogram(WRMSE_uncond); plot([1 1],[0 100]) +histogram(WRMSE); histogram(WRMSE_uncond); plot([1 1],[0 100]); set(gca, 'XScale', 'log') xlabel('WRMSE'); ylabel('Histogram'); legend('Conditional realizations', 'Unconditional realizations', 'True initial field') % export_fig -eps 'misfit-hist' figure(12);clf; colormap(viridis());c_axis=[min(gen.f.output.pseudo(:)) max(gen.f.output.pseudo(:))]; clf; -subplot(3,1,1); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],gen.f.output.pseudo,'filled');set(gca,'Ydir','reverse');caxis(c_axis); xlim([0 100]); ylim([0 16]); colorbar('southoutside'); -subplot(3,1,2); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],mean(fsim_pseudo,2),'filled');set(gca,'Ydir','reverse');caxis(c_axis); colorbar('southoutside');xlim([0 100]); ylim([0 16]) -subplot(3,1,3); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],std(fsim_pseudo,[],2)./mean(fsim_pseudo,2),'filled');set(gca,'Ydir','reverse'); colorbar('southoutside');xlim([0 100]); ylim([0 16]) -% export_fig -eps 'pseudo-sec' - -figure(13);clf;colormap(viridis()); c_axis=[min(gen.f.output.pseudo(:)) max(gen.f.output.pseudo(:))]; clf; -subplot(2,1,1); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],(mean(fsim_pseudo,2)-gen.f.output.pseudo)./gen.f.output.pseudo,'filled');set(gca,'Ydir','reverse'); colorbar('southoutside');xlim([0 100]); ylim([0 16]) -caxis([-.05 .1]) -subplot(2,1,2); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],(gen.i.output.pseudoCalc-gen.f.output.pseudo)./gen.f.output.pseudo,'filled');set(gca,'Ydir','reverse'); colorbar('southoutside');xlim([0 100]); ylim([0 16]) -caxis([-.05 .1]) +subplot(4,1,1); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],gen.f.output.pseudo,'filled');set(gca,'Ydir','reverse');caxis(c_axis); xlim([0 100]); ylim([0 16]); colorbar('southoutside'); +subplot(4,1,2); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],mean(fsim_pseudo,2),'filled');set(gca,'Ydir','reverse');caxis(c_axis); colorbar('southoutside');xlim([0 100]); ylim([0 16]) +subplot(4,1,3); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],std(fsim_pseudo,[],2)./mean(fsim_pseudo,2),'filled');set(gca,'Ydir','reverse'); colorbar('southoutside');xlim([0 100]); ylim([0 16]) +subplot(4,1,4); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],(mean(fsim_pseudo,2)-gen.f.output.pseudo)./gen.f.output.pseudo,'filled');set(gca,'Ydir','reverse'); colorbar('southoutside');xlim([0 100]); ylim([0 16]) % export_fig -eps 'pseudo-sec-err' figure(23); clf; hold on; axis equal tight; @@ -400,8 +416,8 @@ scatter(mean(fsim_pseudo,2),gen.f.output.pseudo,'.g'); x=[floor(min(fsim_pseudo(:))) ceil(max(fsim_pseudo(:)))]; plot(x,x,'-r'); -plot(x,x-x*3*gen.i.b_wgt,'--r'); -plot(x,x+x*3*gen.i.b_wgt,'--r'); +plot(x,x-x*gen.i.b_wgt,'--r'); +plot(x,x+x*gen.i.b_wgt,'--r'); xlabel('Apparent resistivity measured from simulated fields'); ylabel('Apparent resistivity measured from true fields'); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log') @@ -409,85 +425,80 @@ -%% Forward of upscaled +%% Invertion on fine scale grid -filepath='.\data_gen\IO-file-inversion-scaleoftrue\'; +i = gen.i; +i.grid = gen.f.grid; +i.grid.nx = numel(i.grid.x); +i.grid.ny = numel(i.grid.y); +i.elec_spacing = gen.f.elec_spacing; +i.elec_id = gen.f.elec_id; +i_fine_scale = Matlat2R2(i,gen.elec); -f={}; -f.res_matrix = gen.f.res_matrix; -f.grid = gen.i.grid; -f.header = 'Forward'; % title of up to 80 characters -f.job_type = 0; -f.filepath = ['data_gen/IO-file-inversion-CoarseScale/']; -f.readonly = 0; -f.alpha_aniso = gen.f.alpha_aniso; -f.elec_spacing = gen.f.elec_spacing; -f.elec_id = gen.i.elec_id; -f.rho = i.output.res; -f.num_regions = 1+numel(f.rho); -f.rho_min = gen.f.rho_min; -f.rho_avg = gen.f.rho_avg; -f.rho_max = gen.f.rho_max; +% save(['result/' fieldname '_i_fine_scale'],'i_fine_scale','-v7.3') +% load(['result/' fieldname '_i_fine_scale']) -mkdir(f.filepath) -f = Matlat2R2(f,gen.elec); % write file and run forward modeling -misfit = sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) -misfit = sqrt( mean (((f.output.pseudo - i.output.pseudoObs) ./ (i.output.pseudoObs)).^2 )) +%% Inversion base on true model +i = gen.i; +i.grid = gen.f.grid; +i.grid.nx = numel(i.grid.x); +i.grid.ny = numel(i.grid.y); +i.elec_spacing = gen.f.elec_spacing; +i.elec_id = gen.f.elec_id; +fid = fopen([gen.f.filepath 'R2_forward.dat'],'r'); +A = textscan(fid,'%f %f %f %f %f %f %f');fclose(fid); +A{end-1}(2:end) = (1+.02*randn(numel(gen.f.output.resistance),1)).*gen.f.output.resistancewitherror; +fid=fopen([gen.f.filepath 'R2_forward.dat'],'w'); +A2=[A{:}]; +fprintf(fid,'%d\n',A2(1,1)); +for u=2:size(A2,1) + fprintf(fid,'%d %d %d %d %d %.9e %.5f\n',A2(u,:)); +end +fclose(fid); -%% Inverting true model 1 iter vs inverting on fine scale grid +i.max_iterations = 1; +i.tolerance = 1.36; -i={}; -i.grid = gen.f.grid; % -i.grid.nx = numel(i.grid.x); -i.grid.ny = numel(i.grid.y); -i.n_plus = 10; -i.res_matrix = 3; +i.rho = mean(1000./sigma_true(:))*ones( numel(i.grid.y), numel(i.grid.x)); +i.rho(i.grid.inside)= 1000./sigma_true; % f({grid_Rho.y,grid_Rho.x}); -i.elec_spacing = gen.f.elec_spacing; -i.elec_id = gen.f.elec_id; +i_true_mod = Matlat2R2(i,gen.elec); +% save(['result/' fieldname '_i_true_mod'],'i_true_mod','-v7.3') +% load(['result/' fieldname '_i_true_mod']) -i.a_wgt = 0;%0.01; -i.b_wgt = 0.02;%0.02; -% fid = fopen([gen.f.filepath 'R2_forward.dat'],'r'); -% A = textscan(fid,'%f %f %f %f %f %f %f');fclose(fid); -% A{end-1}(2:end) = gen.f.output.resistancewitherror; -% fid=fopen([gen.f.filepath 'R2_forward.dat'],'w'); -% A2=[A{:}]; -% fprintf(fid,'%d\n',A2(1,1)); -% for u=2:size(A2,1) -% fprintf(fid,'%d %d %d %d %d %f %f\n',A2(u,:)); -% end -% fclose(fid); - -i.header = 'Inverse'; % title of up to 80 characters -i.job_type = 1; -i.filepath = 'data_gen/IO-file/'; -i.readonly = 0; -i.alpha_aniso = gen.f.alpha_aniso; -i.tolerance = 1; -i.rho_avg = gen.f.rho_avg; -i.inverse_type = 1; - -% on fine scale -i.max_iterations = 10; -i_fine_scale = Matlat2R2(i,gen.elec); + +%% forward response based on the fine-scale inversion -% on true model -i.max_iterations = 1; -i.tolerance = 1.4; -i.inverse_type = 0; -i.rho = mean(1000./sigma_true(:))*ones( numel(i.grid.y), numel(i.grid.x)); -i.rho(i.grid.inside)= 1000./sigma_true; % f({grid_Rho.y,grid_Rho.x}); +i = i_fine_scale; +i = i_true_mod; -i_true_mod = Matlat2R2(i,gen.elec); +f=gen.f; +f.rho = 1000./flipud(i.output.res); +f = Matlat2R2(f,gen.elec); % write file and run forward modeling +sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) + +Nscore.forward = @(x) (log(x./43)/1.4-gen.mu)./gen.std; +Nscore.inverse = @(x) 43.*exp(1.4*(x.*gen.std+gen.mu)); +% Sec_d = Nscore.forward(1000./flipud(i.output.res)); +Prim.d_full = Nscore.forward(1000./gen.f.rho); +Prim.d = reshape(Prim.d_full(gen.f.grid.inside),grid_gen.ny,grid_gen.nx); +% Prim.x = grid_gen.x; Prim.y = grid_gen.y; +% [Prim.X, Prim.Y] = meshgrid(Prim.x, Prim.y); + + +ztrue = reshape(i.output.Res * Prim.d_full(:),i.grid.ny,i.grid.nx); +f.rho = 1000./Nscore.inverse(ztrue); +f = Matlat2R2(f,gen.elec); % write file and run forward modeling +sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) -Jrfs = reshape(sum(abs(i_fine_scale.output.J),2),size(Prim.d_full)); + +Jrfs = reshape(sum(abs(i.output.J),2),size(Prim.d_full)); Jrtm = reshape(sum(abs(i_true_mod.output.J),2),size(Prim.d_full)); figure; @@ -499,6 +510,62 @@ subplot(2,1,2); scatter(gen.f.pseudo_x,gen.f.pseudo_y,[],mean(i_true_mod.output.J,1),'filled'); set(gca,'ydir','reverse') +%% Forward of upscaled + +filepath='.\data_gen\IO-file-inversion-scaleoftrue\'; + +f=gen.f; +f.grid = gen.i.grid; +f.elec_id = gen.i.elec_id; + +% f(Z^{est}) +f.rho = 1000./Sigma.d; +f.filepath = ['data_gen/IO-file-forward-tomogram/']; +mkdir(f.filepath) +f = Matlat2R2(f,gen.elec); % write file and run forward modeling +sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) +sqrt(mean(((f.output.pseudo - gen.i.output.pseudoObs)./(gen.i.output.pseudoObs) ).^2)) + + +% f( RUz^{true} ) +f.rho = 1000./Nscore.inverse(Test_Sec_d); +f = Matlat2R2(f,gen.elec); % write file and run forward modeling +sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) +sqrt(mean(((f.output.pseudo - gen.i.output.pseudoObs)./(gen.i.output.pseudoObs) ).^2)) + + +% f(Uz^{true}) +U = zeros(numel(Sec.d), numel(Prim.d)); +for i=1:numel(Sec.d) + Res = zeros(numel(Sec.y),numel(Sec.x)); + Res(i)=1; + F = griddedInterpolant({Sec.y,Sec.x},Res,'linear'); + res_t = F({Prim.y,Prim.x}); + U(i,:) = res_t(:) ./sum(res_t(:)); +end +U(isnan(U))=0; +Uztrue = reshape(U*Prim.d(:),numel(Sec.y),numel(Sec.x)); +f.rho = 1000./Sigma.d; +f.rho(f.grid.inside) = 1000./Nscore.inverse(Uztrue(f.grid.inside)); +f = Matlat2R2(f,gen.elec); % write file and run forward modeling +sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) +sqrt(mean(((f.output.pseudo - gen.i.output.pseudoObs)./(gen.i.output.pseudoObs) ).^2)) + + + + +%% Forward true field + +f=gen.f; +F = griddedInterpolant({Sigma.y,Sigma.x},1000./Sigma.d,'linear'); +f.rho = F({f.grid.y, f.grid.x}); +f.rho(f.grid.inside)= 1000./sigma_true; +f = Matlat2R2(f,gen.elec); % write file and run forward modeling +sqrt( mean (((f.output.resistance - gen.f.output.resistancewitherror) ./ (gen.f.output.resistancewitherror)).^2 )) +sqrt(mean(((f.output.pseudo - gen.i.output.pseudoObs)./(gen.i.output.pseudoObs) ).^2)) + + + %% Figure for Synthetic schema figure(199);clf; n=5; diff --git a/HT/covar2W.m b/HT/covar2W.m new file mode 100644 index 0000000..37dfb78 --- /dev/null +++ b/HT/covar2W.m @@ -0,0 +1,25 @@ +function W = covar2W(covar,Prim,Sec,G,Prim_pt,Cmt) + +% Compute the covariance of the spatial model +Cz = covar.g(squareform(pdist([Prim.X(:) Prim.Y(:)]*covar.cx))); +Cz=sparse(Cz); +Czd = Cz * G'; +Cd = G * Czd; + + +% Combine both covariance an built the Kriging System +Cd2 = Cd+Cmt; +Czh = Cz(Prim_pt.id,Prim_pt.id); +Czzh = Cz(Prim_pt.id,:); +Czhd = Czd( Prim_pt.id ,:); +CCa = [ Cd2 Czhd' ; Czhd Czh ]; +CCb = [ Czd' ; Czzh' ]; + +% Solve the kriging system +W=zeros(numel(Prim.x)*numel(Prim.y),numel(Sec.d(:))+Prim_pt.n); +% parpool(48) +parfor ij=1:numel(Prim.x)*numel(Prim.y) + W(ij,:) = CCa \ CCb(:,ij); +end + +end \ No newline at end of file diff --git a/HT/data_generation.m b/HT/data_generation.m index 5c0eef5..2bb1646 100644 --- a/HT/data_generation.m +++ b/HT/data_generation.m @@ -42,12 +42,13 @@ assert(isfield(gen, 'ny')) if ~isfield(gen, 'method'); gen.method = 'Random'; end if ~isfield(gen, 'samp'); gen.samp = 2; end -if ~isfield(gen, 'samp_n'); gen.samp_n = 1/100 * (2^gen.sx+1)*(2^gen.sy+1); end +if ~isfield(gen, 'samp_n'); gen.samp_n = round(1/100 * gen.nx*gen.ny); end if ~isfield(gen, 'mu'); gen.mu = 0; end if ~isfield(gen, 'std'); gen.std = 1; end % Other if ~isfield(gen, 'plotit'); gen.plotit = 0; end if ~isfield(gen, 'saveit'); gen.saveit = 1; end +if ~isfield(gen, 'forwardonly'); gen.saveit = 0; end if ~isfield(gen, 'name'); gen.name = ''; end if ~isfield(gen, 'seed'); gen.seed = 'default'; end if ~isfield(gen, 'plot'); gen.plot = true; end @@ -174,16 +175,19 @@ f = Matlat2R2(f,elec); % write file and run forward modeling -if 0 && gen.saveit +if gen.forwardonly gen.Rho.f = f; gen.Rho.elec = elec; filename = ['data_gen/data/FOR-', gen.name ,'_', datestr(now,'yyyy-mm-dd_HH-MM'), '.mat']; + if gen.saveit save(filename, 'K_true', 'grid_gen', 'gen') %'sigma', + end + keyboard end % Add some error to the observation i.a_wgt = 0;%0.01; -i.b_wgt = 0.01; +i.b_wgt = 0.05; % var(R) = (a_wgt*a_wgt) + (b_wgt*b_wgt) * (R*R) f.output.resistancewitherror = i.a_wgt.*randn(numel(f.output.resistance),1) + (1+i.b_wgt*randn(numel(f.output.resistance),1)).*f.output.resistance; %f.output.resistancewitherror(f.output.resistancewitherror>0) = -f.output.resistancewitherror(f.output.resistancewitherror>0); diff --git a/HT/script_HT.m b/HT/script_HT.m index df56f3f..d508b0d 100644 --- a/HT/script_HT.m +++ b/HT/script_HT.m @@ -16,15 +16,14 @@ % Generation parameter gen.samp = 1; % Method of sampling of K and g | 1: borehole, 2:random. For fromK or from Rho only gen.samp_n = 0; % number of well or number of point -gen.covar(1).model = 'k-bessel'; +gen.covar(1).model = 'matern'; gen.covar(1).alpha = 0.5; -gen.covar(1).range0 = [5 20]; +gen.covar(1).range0 = [4 16]; gen.covar(1).azimuth = 0; gen.covar(1).c0 = 1; gen.covar = kriginginitiaite(gen.covar); gen.mu = 0.25;%0.27; % parameter of the first field. gen.std = 0.05;%.05; -gen.Rho.method = 'R2'; % 'Paolo' (default for gen.method Paolo), 'noise', 'RESINV3D' % Electrical inversion gen.Rho.f ={}; @@ -33,6 +32,8 @@ gen.Rho.elec.bufzone_y = 2; % number of electrod to skip gen.Rho.elec.x_t = 25; % in unit [m] | adapted to fit the fine grid gen.Rho.elec.x_r = [10 40]; % in unit [m] | adapted to fit the fine grid +% gen.Rho.elec.x_t = 10; % in unit [m] | adapted to fit the fine grid +% gen.Rho.elec.x_r = [11 12 15 20 25 30 40]; % in unit [m] | adapted to fit the fine grid gen.Rho.elec.config_max = 6000; % number of configuration of electrode maximal gen.Rho.i.grid.nx = (gen.nx-1)/2+1; % | adapted to fit the fine grid gen.Rho.i.grid.ny = (gen.ny-1)/2+1; % log-spaced grid | adapted to fit the fine grid @@ -42,7 +43,7 @@ gen.plotit = true; % display graphic or not (you can still display later with |script_plot.m|) gen.saveit = true; % save the generated file or not, this will be turn off if mehod Paolo or filename are selected gen.name = '60x20'; -gen.seed = 8; +gen.seed = 11; % Run the function fieldname = data_generation(gen); @@ -52,15 +53,51 @@ %% Reproducing Theim equation +gen.xmax = 100; %total length in unit [m] +gen.ymax = 100; %total hight in unit [m] +gen.nx = gen.xmax*2+1; +gen.ny = gen.ymax*2+1; + +% Generation parameter +gen.covar(1).model = 'matern'; +gen.covar(1).alpha = 0.5; +gen.covar(1).range0 = [4 16]; +gen.covar(1).azimuth = 0; +gen.covar(1).c0 = 1; +gen.covar = kriginginitiaite(gen.covar); +gen.mu = (log10(10^-5)+4.97)/6.66; +gen.std = 0; + +% Electrical inversion +gen.Rho.f ={}; +gen.Rho.f.res_matrix = 0; +gen.Rho.elec.spacing_y = 1; % in unit [m] | adapted to fit the fine grid +gen.Rho.elec.bufzone_y = 30; % number of electrod to skip +gen.Rho.elec.x_t = 50; % in unit [m] | adapted to fit the fine grid +gen.Rho.elec.x_r = [51 52 55 60 75 90 100]; % in unit [m] | adapted to fit the fine grid +gen.Rho.elec.config_max = 6000; % number of configuration of electrode maximal +gen.Rho.i.grid.nx = (gen.nx-1)/2+1; % | adapted to fit the fine grid +gen.Rho.i.grid.ny = (gen.ny-1)/2+1; % log-spaced grid | adapted to fit the fine grid +gen.Rho.i.res_matrix = 3; + +% Other parameter +gen.forwardonly = 1; +gen.plotit = true; % display graphic or not (you can still display later with |script_plot.m|) +gen.saveit = true; % save the generated file or not, this will be turn off if mehod Paolo or filename are selected +gen.name = '100x100'; +gen.seed = 11; + +% Run the function +fieldname = data_generation(gen); + % Uniform K, 1 pumping well and 3 measuring well separated by 10 m each. -load('data_gen/data/FOR-60x20_2018-02-27_09-34.mat') +load('result/FOR-100x100_2018-04-04_14-23.mat') % Thiem equation % % $$h_{2}-h_{1}={\frac {Q}{2\pi b K}}\ln \left({\frac {r_{1}}{r_{2}}}\right)$$ % - uex = unique(gen.Rho.elec.X(gen.Rho.elec.data(:,1))); id = bsxfun(@eq,gen.Rho.elec.X(gen.Rho.elec.data(:,1))',uex ); h = id* gen.Rho.f.output.resistance./sum(id,2); @@ -83,56 +120,65 @@ % $$h_{2}-h_{1}={\frac {Q}{4\pi K}}\ln \left({\frac {1}{r_{2}} - \frac {1}{r_{1}}}\right)$$ -[~,id] = min((gen.Rho.elec.Y(gen.Rho.elec.data(:,3))-30).^2); -id = gen.Rho.elec.Y(gen.Rho.elec.data(:,3))==gen.Rho.elec.Y(gen.Rho.elec.data(id,3)); +[~,id0] = min((gen.Rho.elec.Y(gen.Rho.elec.data(:,3))-50).^2); +id = gen.Rho.elec.Y(gen.Rho.elec.data(:,3))==gen.Rho.elec.Y(gen.Rho.elec.data(id0,3)); h = gen.Rho.f.output.resistance(id); x=gen.Rho.elec.X(gen.Rho.elec.data(id,1)); y=gen.Rho.elec.Y(gen.Rho.elec.data(id,1)); -x0=gen.Rho.elec.X(gen.Rho.elec.data(id,3)); -y0=gen.Rho.elec.Y(gen.Rho.elec.data(id,3)); +x0=gen.Rho.elec.X(gen.Rho.elec.data(id0,3)); +y0=gen.Rho.elec.Y(gen.Rho.elec.data(id0,3)); figure(2); hold on; scatter(x, y,[],log(h),'filled') plot(x0,y0,'xk') +rectangle('Position',[0 0 gen.xmax gen.ymax]) + Q=1; K=mean(K_true(:)); r=sqrt( (x-x0).^2 + (y-y0).^2 ); -R= 5; H=mean(h(r==R)); -ht = -Q/(4*pi*K).*(-1./r+1/R) + H; + +theim = @(Q,K,hmax,r) -Q/(4*pi*K).*(1./Inf-1./r) + hmax; + +rmse = @(x) sqrt( sum( ( theim(Q,K,x,r) - h ).^2 ) ); +h0 = fminsearch(rmse,min(h)); + figure; hold on; scatter3(x,y,h,'.') -scatter3(x,y,ht,'.') +scatter3(x,y,theim(Q,K,h0,r),'.') +rectangle('Position',[0 0 gen.xmax gen.ymax]); view(3) + + +R = sqrt( (grid_gen.X-x0).^2 + (grid_gen.Y-y0).^2 ); + +figure; hold on; +s=surf(grid_gen.X, grid_gen.Y, theim(Q,K,h0,R)); +scatter3(x,y,h,'.k') +s.EdgeColor='none'; +set(gca,'zscale','log') +xlabel('X [m]'); ylabel('Y [m]'); zlabel('Head [log-m]'); +legend({'Analytical solution','Observations'}) + -%% Inversion -fieldname='GEN-60x20_2018-03-12_15-49'; -load(['data_gen/data/' fieldname '.mat']) + + +%% Testing other geostat param + +fieldname='GEN-60x20_2018-03-13_14-26'; +load(['result/' fieldname '.mat']) addpath('../functions','R2'); % Normal Score based on known distribution of Prim and Sec Nscore.forward = @(x) ( (log10(x)+4.97)./6.66 - gen.mu)./gen.std; -Nscore.inverse = @(x) 6.66*10^(x.*gen.std+gen.mu)-4.97; +Nscore.inverse = @(x) 10.^(6.66*(x.*gen.std+gen.mu)-4.97); Sec=K; Sec.d = Nscore.forward(Sec.d); Prim.d = Nscore.forward(K_true); Prim.x = grid_gen.x; Prim.y = grid_gen.y; Prim.X = grid_gen.X; Prim.Y = grid_gen.Y; - -% Figure of intial data -figure(1); clf;c_axis=[ min(Prim.d(:)) max(Prim.d(:)) ]; -subplot(2,1,1);imagesc(grid_gen.x, grid_gen.y, Prim.d); caxis(c_axis); title('True field'); axis tight equal; -subplot(2,1,2); hold on; -imagesc(Sec.x, Sec.y, Sec.d); title('Inverted field'); box on;colormap(viridis()); axis tight equal;set(gca,'Ydir','reverse') -plot(gen.Rho.elec.x_t,gen.Rho.elec.y,'or') -plot(gen.Rho.elec.x_r(1),gen.Rho.elec.y,'xb') -plot(gen.Rho.elec.x_r(2),gen.Rho.elec.y,'xb') -% export_fig -eps 'Prim_and_sec' - - - % Built the matrix G which link the true variable Prim.d to the measured coarse scale d G = zeros(numel(Sec.d), numel(Prim.d)); for i=1:numel(Sec.d) @@ -142,66 +188,140 @@ G(i,:) = res_t(:) ./sum(res_t(:)); end -Test_Sec_d = reshape(G * Prim.d(:), numel(Sec.y), numel(Sec.x)); -figure(4); clf; colormap(viridis()) -subplot(3,1,1); -surface(Sec.X,Sec.Y,Sec.d,'EdgeColor','none','facecolor','flat'); -view(2); set(gca,'Ydir','reverse'); caxis([-1.5 1.5]); axis equal tight; box on; xlabel('x');ylabel('y'); colorbar; title('Inverted field Z^{est}') -subplot(3,1,2);surface(Sec.X,Sec.Y,Test_Sec_d,'EdgeColor','none','facecolor','flat'); -view(2); set(gca,'Ydir','reverse'); caxis([-1.5 1.5]);axis equal tight; box on; xlabel('x');ylabel('y'); colorbar; title('G-transform of the true field Gz^{true}') -subplot(3,1,3);surface(Sec.X,Sec.Y,Test_Sec_d-Sec.d,'EdgeColor','none','facecolor','flat'); -view(2); set(gca,'Ydir','reverse'); axis equal tight; box on; xlabel('x');ylabel('y');colorbar; title('Error Gz^{true}-Z^{est}') -% export_fig -eps 'SecvsTestSecD' +covar = gen.covar; +covar.range0 = covar.range0/2; +covar.alpha = covar.alpha/2; +covar = kriginginitiaite(covar); -% Generate a sampling -Prim_pt = sampling_pt(Prim,Prim.d,2,0); +plot(Prim.x,1-covar.g(Prim.x*covar.cx(1))) +xlim([0 50]) +% Compute the weightings matrix +W = covar2W(covar,Prim,Sec,G,Prim_pt,Cmt); -% Compute the covariance of the data error -Cm = inv(sqrtm(full(gen.Rho.i.output.R(gen.Rho.i.output.inside,gen.Rho.i.output.inside)))); -Cmt = (eye(size(Sec.res))-Sec.res)*Cm; +% save(['result/' fieldname '_cond_covar_rx2'],'W','Prim_pt','G','Nscore','Sec','Prim','covar') +save(['result/' fieldname '_cond'],'W','Prim_pt','G','Nscore','Sec','Prim','covar') -% Compute the covariance of the spatial model -covar = kriginginitiaite(gen.covar); -Cz = covar.g(squareform(pdist([Prim.X(:) Prim.Y(:)]*covar.cx))); -Cz=sparse(Cz); -Czd = Cz * G'; -Cd = G * Czd; +% Load all files +listing = dir(['result/' fieldname '_cond*']); +parm.n_real = 30; +err=nan(numel(listing),parm.n_real); -% Combine both covariance an built the Kriging System -Cd2 = Cd+Cmt; -Czh = Cz(Prim_pt.id,Prim_pt.id); -Czzh = Cz(Prim_pt.id,:); -Czhd = Czd( Prim_pt.id ,:); -CCa = [ Cd2 Czhd' ; Czhd Czh ]; -CCb = [ Czd' ; Czzh' ]; - +for i_l = 1:numel(listing) + + % Load + load([listing(i_l).folder '\' listing(i_l).name]); + Nscore.inverse = @(x) 10.^(6.66*(x.*gen.std+gen.mu)-4.97); + + % Compute the Kriging map + zh = reshape( W * [Sec.d(:) ; Prim_pt.d], numel(Prim.y), numel(Prim.x)); + + % Compute the Realizations + zcs=nan(numel(Prim.y), numel(Prim.x),parm.n_real); + z=nan(numel(Sec.y), numel(Sec.x),parm.n_real); + for i_real=1:parm.n_real + zs = fftma_perso(covar, struct('x',Prim.x,'y',Prim.y)); + zhs = W * [G * zs(:) ; zs(Prim_pt.id)./1.3]; + r = zh(:) + (zs(:) - zhs(:)); + zcs(:,:,i_real) = reshape( r, numel(Prim.y), numel(Prim.x)); + z(:,:,i_real)=reshape(G*r(:), numel(Sec.y), numel(Sec.x) ); + end + + % Compute the Vario & Histo + vario_x=nan(parm.n_real,numel(Prim.x)); + vario_y=nan(parm.n_real,numel(Prim.y)); + xi = -3:.1:3; + f=nan(parm.n_real,numel(xi)); + for i_real=1:parm.n_real + r = zcs(:,:,i_real); + %r = (r(:)-mean(r(:)))./std(r(:)); + [vario_x(i_real,:),vario_y(i_real,:)]=variogram_gridded_perso(reshape( r(:), numel(Prim.y), numel(Prim.x))); + + r=zcs(:,:,i_real); + f(i_real,:) = ksdensity(r(:),xi); + end + [vario_prim_x,vario_prim_y]=variogram_gridded_perso(Prim.d); + + + % Compute the forward response + fsim_pseudo=nan(numel(gen.Rho.f.output.pseudo),parm.n_real); + fsim_resistance=nan(numel(gen.Rho.f.output.resistance),parm.n_real); + rho = 1./Nscore.inverse(zcs); + parfor i_real=1:parm.n_real + f={}; + f.res_matrix = gen.Rho.f.res_matrix; + f.grid = gen.Rho.f.grid; + f.header = 'Forward'; % title of up to 80 characters + f.job_type = 0; + f.filepath = ['data_gen/IO-file-' num2str(i_real) '/']; + f.readonly = 0; + f.alpha_aniso = gen.Rho.f.alpha_aniso; + f.elec_X_id = gen.Rho.f.elec_X_id; + f.elec_Y_id = gen.Rho.f.elec_Y_id; + f.rho = rho(:,:,i_real); + f.num_regions = gen.Rho.f.num_regions; + f.rho_min = gen.Rho.f.rho_min; + f.rho_avg = gen.Rho.f.rho_avg; + f.rho_max = gen.Rho.f.rho_max; + + mkdir(f.filepath) + f = Matlat2R2(f,gen.Rho.elec); % write file and run forward modeling + fsim_pseudo(:,i_real) = f.output.pseudo; + fsim_resistance(:,i_real) = f.output.resistance; + end + + % Compute the misfit + fsim_misfit=nan(numel(gen.Rho.f.output.resistancewitherror),parm.n_real); + for i_real=1:parm.n_real + fsim_misfit(:,i_real) = (fsim_resistance(:,i_real) - gen.Rho.f.output.resistancewitherror) ./ (gen.Rho.i.a_wgt + gen.Rho.i.b_wgt*gen.Rho.f.output.resistancewitherror); + err(i_l, i_real) = sqrt(mean(fsim_misfit(:,i_real).^2)); + end -% Solve the kriging system -W=zeros(numel(Prim.x)*numel(Prim.y),numel(Sec.d(:))+Prim_pt.n); -parfor ij=1:numel(Prim.x)*numel(Prim.y) - W(ij,:) = CCa \ CCb(:,ij); end -% save(['ERT/result/' fieldname '_cond'],'W','Prim_pt','G','Nscore','Sec','Prim') -% Compute the Kriging map -zh = reshape( W * [Sec.d(:) ; Prim_pt.d], numel(Prim.y), numel(Prim.x)); -%zhtest = reshape( W * [Test_Sec_d(:) ; Prim_pt.d], ny, nx); -figure(5); clf; colormap(viridis()) -surf(Prim.x,Prim.y,zh,'EdgeColor','none','facecolor','flat'); caxis([-3 3]) -view(2); axis equal tight; set(gca,'Ydir','reverse'); xlabel('x');ylabel('y'); colorbar('southoutside'); title('Kriging Estimate') -% export_fig -eps 'Krig' +figure; boxplot(err','Labels', {listing.name},'LabelOrientation','inline') + + + + + +%% GW + +kn=numel(K_true); + +%hydraulic conductivity +k=[(1:kn)' ones(kn,1)*0.1 K_true(:) zeros(kn,1) K_true(:) zeros(kn,2) K_true(:) ones(kn,1)*0.1e-5]; +dlmwrite('k.dat',k,'delimiter',' ','precision','%6.12f'); + +%mass properties +mp=[(1:kn)' ones(kn,1)*-4 ones(kn,1) ones(kn,1)*0.1 zeros(kn,1) ones(kn,1)*2e-9 zeros(kn,2) zeros(kn,1)]; +dlmwrite('mp.dat',mp,'delimiter',' ','precision','%6.12f'); + +% run gw +!gw3.1.6_Win32.exe + + + + + + + +%% ------------------------------------------------------------------------ +% Previous Script +%-------------------------------------------------------------------------- + + %% Simulation of the Area-to-point Kriging rng('shuffle'); parm.n_real = 30; -covar = kriginginitiaite(gen.covar); +% covar = kriginginitiaite(gen.covar); zcs=nan(numel(Prim.y), numel(Prim.x),parm.n_real); z=nan(numel(Sec.y), numel(Sec.x),parm.n_real); for i_real=1:parm.n_real @@ -223,9 +343,9 @@ subplot(4,1,4);surf(Prim.x, Prim.y, std(zcs,[],3),'EdgeColor','none','facecolor','flat'); view(2); axis tight equal; set(gca,'Ydir','reverse'); colorbar; %export_fig -eps 'PrimOverview' -figure(65);clf; colormap(viridis()) +figure;clf; colormap(viridis()) c_axis=[ -3 3]; -subplot(4,1,1);surf(Prim.x, Prim.y, Prim.d,'EdgeColor','none','facecolor','flat'); caxis(c_axis);view(2); axis tight equal; set(gca,'Ydir','reverse'); colorbar; +subplot(4,1,1);surf(Prim.x, Prim.y, Prim.d,'EdgeColor','none','facecolor','flat'); caxis(c_axis);view(2); axis tight equal; set(gca,'Ydir','reverse'); % hold on; scatter(Prim_pt.x,Prim_pt.y,'filled','r') subplot(4,1,2);surf(Prim.x, Prim.y, zcs(:,:,1),'EdgeColor','none','facecolor','flat'); caxis(c_axis);view(2); axis tight equal; set(gca,'Ydir','reverse'); subplot(4,1,3);surf(Prim.x, Prim.y, zcs(:,:,2),'EdgeColor','none','facecolor','flat'); caxis(c_axis);view(2); axis tight equal; set(gca,'Ydir','reverse'); @@ -240,6 +360,14 @@ subplot(4,1,4);surf(Sec.x, Sec.y, std(z,[],3),'EdgeColor','none','facecolor','flat'); title('d Average of True field');view(2); axis tight equal; set(gca,'Ydir','reverse'); colorbar; %export_fig -eps 'SecOverview' +figure;clf; colormap(viridis()) +c_axis=[ -1 1]; +subplot(4,1,1);surf(Sec.x, Sec.y, Sec.d,'EdgeColor','none','facecolor','flat'); caxis(c_axis); view(2); axis tight equal; set(gca,'Ydir','reverse'); +subplot(4,1,2);surf(Sec.x, Sec.y, z(:,:,1),'EdgeColor','none','facecolor','flat'); caxis(c_axis); view(2); axis tight equal; set(gca,'Ydir','reverse'); +subplot(4,1,3);surf(Sec.x, Sec.y, z(:,:,2),'EdgeColor','none','facecolor','flat'); caxis(c_axis); view(2); axis tight equal; set(gca,'Ydir','reverse'); +subplot(4,1,4);surf(Sec.x, Sec.y, z(:,:,3),'EdgeColor','none','facecolor','flat'); caxis(c_axis); view(2); axis tight equal; set(gca,'Ydir','reverse'); +%export_fig -eps 'SecOverview' + figure(8);clf;colormap(viridis()) subplot(2,1,1);surface(Sec.X,Sec.Y,Test_Sec_d-Sec.d,'EdgeColor','none','facecolor','flat'); view(2); set(gca,'Ydir','reverse'); axis equal tight; box on; xlabel('x');ylabel('y'); caxis([-.6 .6]);colorbar('southoutside'); subplot(2,1,2);surf(Sec.x, Sec.y, mean(z,3)-Sec.d,'EdgeColor','none','facecolor','flat'); view(2); axis tight equal; set(gca,'Ydir','reverse'); colorbar('southoutside');caxis([-.6 .6]) @@ -261,7 +389,7 @@ h1=plot(Prim.x(1:2:end),vario_x(:,1:2:end)','b','color',[.5 .5 .5]); h2=plot(Prim.x,vario_prim_x,'-r','linewidth',2); h3=plot(Prim.x,1-covar.g(Prim.x*covar.cx(1)),'--k','linewidth',2); -xlim([0 60]); xlabel('Lag-distance h_x ');ylabel('Variogram \gamma(h_x)') +xlim([0 30]); xlabel('Lag-distance h_x ');ylabel('Variogram \gamma(h_x)') legend([h1(1) h2 h3],'500 realizations','True field','Theorical Model') subplot(2,1,2); hold on; h1=plot(Prim.y,vario_y','b','color',[.5 .5 .5]); @@ -269,7 +397,7 @@ h3=plot(Prim.y,1-covar.g(Prim.y*covar.cx(4)),'--k','linewidth',2); xlim([0 6]); xlabel('Lag-distance h_y ');ylabel('Variogram \gamma(h_y)') legend([h1(1) h2 h3],'500 realizations','True field','Theorical Model') -export_fig -eps 'Vario' +% export_fig -eps 'Vario' figure(10);clf; hold on; colormap(viridis())