diff --git a/TestSampling2PDF.m b/TestSampling2PDF.m new file mode 100644 index 0000000..e4aeda2 --- /dev/null +++ b/TestSampling2PDF.m @@ -0,0 +1,47 @@ + +fun = @(x,mu,sigma) 1/sqrt(2*pi*sigma^2) * exp(-(x-mu).^2/(2*sigma^2)); + +mu=[1 1]; +s=[1 1]; +w=[1 1]; + +fun2 = @(x) fun(x,mu(1),s(1)).^w(1) .* fun(x,mu(2),s(2)).^w(2); +fun3 = @(u) integral(fun2,-Inf,u) / integral(fun2,-Inf,Inf); + + +funA = @(u) integral( @(x) (1/sqrt(2*pi*s(1)^2) * exp(-(x-mu(1)).^2/(2*s(1)^2))).^w(1) .* (1/sqrt(2*pi*s(2)^2) * exp(-(x-mu(2)).^2/(2*s(2)^2))).^w(2),-Inf,u,'RelTol',0.1,'ArrayValued',true); + + +opt = optimset('Display','off','TolFun',.01); +n=(513*65)*1; +u=nan(1,n); +tic +for i=1:n + u(i) = fsolve(funA,rand(1),opt); +end +toc + +xx=-3:0.1:3; + +figure; hold on; +plot(xx,fun(xx,mu(1),s(1))) +plot(xx,fun(xx,mu(2),s(2))) +plot(xx,fun2(xx)) + + + +syms x + +for i=1:n +mu=[1 1]; +s=[1 1]; +w=[1 1]; + +fun = (1/sqrt(2*pi*s(1)^2) * exp(-(x-mu(1)).^2/(2*s(1)^2))).^w(1) .* (1/sqrt(2*pi*s(2)^2) * exp(-(x-mu(2)).^2/(2*s(2)^2))).^w(2); + +intfun = int(fun,x); + +% fplot(intfun); + +solve(intfun == 0); +end \ No newline at end of file diff --git a/script_BSS.m b/script_BSS.m index 7d3c6c0..ca83c36 100644 --- a/script_BSS.m +++ b/script_BSS.m @@ -101,8 +101,6 @@ kern.dens = reshape(mean(jpdf,3),numel(kern.axis_prim), numel(kern.axis_sec)); - - save(['result-BSGS/' file],'-append','kern'); %% BSGS diff --git a/script_BSS_old.m b/script_BSS_old.m new file mode 100644 index 0000000..8700ed2 --- /dev/null +++ b/script_BSS_old.m @@ -0,0 +1,287 @@ +clc; % clear the command line +addpath(genpath('./.')); % Add folder and sub-folder to path +dbstop if error % activate debug in error + +%% DATA CREATION +% This section gather all possible way to create the data. |gen| struct +% store the parameter and |data_generation.m| compute everything + +% Grid size +gen.xmax = 240; %total length in unit [m] +gen.ymax = 20; %total hight in unit [m] + +% Scale define the subdivision of the grid (multigrid). At each scale, the +% grid size is $(2^gen.sx+1) \times (2^gen.sy+1)$ +gen.sx = 10; +gen.sy = 7; + +% Generation Method: All method generate with FFTMA a gaussian field. +% 'Normal' with normal distribution \sim N(gen.mu,gen.std) +% 'LogNormal' +% 'fromRho': log transform it with the parameter defined below +% 'fromK': generate with FFTMA a field of Hyraulic conductivity and log transform it with the parameter defined below +gen.method = 'fromPhi'; + +% 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 = 4; % number of well or number of point +gen.covar(1).model = 'exponential'; +gen.covar(1).range0 = [27 2.7]; +gen.covar(1).azimuth = 0; +gen.covar(1).c0 = 1; +gen.covar = kriginginitiaite(gen.covar); +gen.mu = 0.27; % parameter of the first field. +gen.std = .05; +gen.Rho.method = 'R2'; % 'Paolo' (default for gen.method Paolo), 'noise', 'RESINV3D' + +% Electrical inversion +gen.Rho.grid.nx = 240; +gen.Rho.grid.ny = 20; % log-spaced grid. +gen.Rho.elec.spacing = 2; % in grid spacing unit. +gen.Rho.elec.config_max = 6000; % number of configuration of electrode maximal +gen.Rho.dmin.res_matrix = 2; % resolution matrix: 1-'sensitivity' matrix, 2-true resolution matrix or 0-none +gen.Rho.dmin.tolerance = 10; + +% Other parameter +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 = 'test'; +gen.seed = 'default'; + +% Run the function +data_generation(gen); +%[fieldname, grid_gen, K_true, phi_true, sigma_true, K, sigma, Sigma, gen] = data_generation(gen); + +%% Create joint-pdf +file='GEN-Run_1_2017-05-07_14-37'; +files={'GEN-Run_2_2017-05-07_21-56','GEN-Run_3_2017-05-06_18-21','GEN-Run_4_2017-05-07_17-17','GEN-Run_5_2017-05-07_22-56','GEN-Run_6_2017-05-07_21-08','GEN-Run_7_2017-05-07_21-04'}; + +%load(['result-BSS/' file],'K_true','Sigma'); + +% Correct +% grid_G=gen.Rho.grid; +% for i_files = 1: numel(files) +% data=dlmread(['Y:\BSGS\result-BSS\data_gen\Run_' num2str(i_files+1) '\IO-file\f001_res.dat']); +% output.res=flipud(reshape(data(:,3),grid_G.ny,grid_G.nx)); +% Rho.d_raw = flipud(output.res); +% f = griddedInterpolant({grid_G.y,grid_G.x},Rho.d_raw,'nearest','nearest'); +% Rho.d = f({grid_gen.y,grid_gen.x}); +% Sigma.d = 1000./Rho.d; +% Sigma.d_raw = 1000./Rho.d_raw; +% save(['result-BSS/' files{i_files}],'-append','Sigma'); +% end + +for i_files = 1: numel(files) + load(['result-BSGS/' files{i_files}],'K_true','Sigma'); + K_true_log_list(:,i_files) = log(K_true(:)); + Sigma_d_list(:,i_files) =Sigma.d(:); +end + + +% Scott's rules +% dS = 3.5*std(Sigma_d_list(:))*numel(Sigma_d_list)^(-1/3); +% dK = 3.5*std(K_true_log_list(:))*numel(K_true_log_list(:))^(-1/3); +% kern.axis_sec = min(Sigma_d_list(:)):dS:max(Sigma_d_list(:)); +% kern.axis_prim = min(K_true_log_list(:)):dK:max(K_true_log_list(:)); + +[kern.prior, kern.axis_prim] = ksdensity(K_true_log_list(:)); +[~ , kern.axis_sec] = ksdensity(Sigma_d_list(:)); + +[X,Y] = meshgrid(kern.axis_sec, kern.axis_prim); +kern.XY = [X(:),Y(:)]; + +for i_files = 1: numel(files) + jpdf(:,:,i_files) = ksdensity([Sigma_d_list(:,i_files) K_true_log_list(:,i_files)],kern.XY); +end + +for i_files = 1: numel(files) + figure; + imagesc( reshape(jpdf(:,:,i_files),numel(kern.axis_prim), numel(kern.axis_sec))); +end + +kern.dens = reshape(mean(jpdf,3),numel(kern.axis_prim), numel(kern.axis_sec)); + + + +save(['result-BSGS/' file],'-append','kern'); + +%% BSGS +% Generation of the high resolution electrical conductivity (sSigma) from +% scarse electrical data (sigma) and large scale inverted ERt (Sigma). +clear all; close all; addpath(genpath('./.')); +load('result-BSS/GEN-Run_1_2017-05-07_14-37'); + +Nscore = nscore(kern, struct('nscore', 1), 0); %Prim.d, kern.axis_prim, 'pchip', 'pchip', parm.plot.ns + +[~,s_id]=min(bsxfun(@minus,kern.axis_sec,Sigma.d(:)).^2,[],2); +sec_pdf = kern.dens(:,s_id); +sec.pdf = bsxfun(@times, sec_pdf, 1./sum(sec_pdf)); +sec.axis = Nscore.forward(kern.axis_prim); + +parm.k.covar = gen.covar; +parm.k.covar.range0 = fliplr(gen.covar.range0) ./ [grid_gen.dy grid_gen.dx]; + +parm.saveit = false; + +parm.seed_path = 'shuffle'; +parm.seed_search = 'shuffle'; +parm.seed_U = 'shuffle'; + +parm.k.wradius = 3; +parm.k.lookup = false; +parm.k.nb = 80; +parm.n_real = 11*2; +parm.mg=1; + +parm.aggr.sum = 1; +parm.aggr.T = .5; +parm.aggr.method = 'cst'; + +% use the log of hyd. cond. +hd = sampling_pt(struct('x',1:grid_gen.nx,'y',1:grid_gen.ny),K_true,2,0); +hd.d = Nscore.forward(hd.d); + +f0=kern.prior ./ sum(kern.prior); +nx = grid_gen.nx; +ny = grid_gen.ny; + + +[Res] = BSS(nx,ny,hd,f0,sec,parm); + +% image +figure(1); clf; +n=min([5 parm.n_real+1]); +subplot(n,1,1); imagesc(log(K_true));colorbar; +for i_real=1:n-1 + subplot(n,1,1+i_real); hold on; + imagesc(Res(:,:,i_real));colorbar; + scatter(hd.x,hd.y,[],hd.d,'filled') + axis tight; +end + + +% variogram +figure(2); clf; hold on; +id = grid_gen.xi_scale) hd.x(hd.scale>i_scale)]; + + for i_pt = start(i_scale)+(1:nb(i_scale)) + + % Compute distance to the hard data + hd_XY_d = bsxfun(@minus,hd_XY_s,XY_i(i_pt,:)); + hd_XY_d = hd_XY_d(hd_XY_d(:,1)0 && ijt(2)>0 && ijt(1)<=ny && ijt(2)<=nx + if Path(ijt(1),ijt(2)) < i_pt % check if it,jt exist + n=n+1; + neigh(n) = nn; + NEIGH_1(n) = ijt(1); + NEIGH_2(n) = ijt(2); + if n >= k_nb + break; + end + end + end + end + + + % Covariance computation + if n==0 + S(i_pt) = k_covar_c0; + else + NEIGH(i_pt,:) = NEIGH_1 + (NEIGH_2-1)* ny; + D = pdist([0 0; ss_hd_XY_s_s(neigh(1:n),:)]*k.covar.cx); + C = k.covar.g(D); + + if n==1 + a0_C = C; + ab_C = 1; + else + a0_C = C(1:n)'; + % Equivalent to : squareform(C(n+1:end)); + ab_C = diag(ones(n,1))*0.5; + ab_C(tril(true(n),-1)) = C(n+1:end); + ab_C = ab_C + ab_C'; + end + + % Weights and variance error + l = ab_C \ a0_C; + LAMBDA(i_pt,:) = [l; nan(k.nb-n,1)]; + S(i_pt) = k_covar_c0 - l'*a0_C; + end + end + % disp(['scale: ' num2str(i_scale) '/' num2str(sn)]) +end +t.weight = toc(tik.weight); +disp(['Weights Computed in ' num2str(t.weight*60)] ) + + +%% Prepare OF +id = grid_gen.xthr),l(l(:)>thr),'.k') +end + + +figure(3); clf; +for i_s=1:sn + subplot(2,5,i_s); hold on; xlabel(num2str(i_s)) + tis = reshape(t(:,i_s,:),nsamples*n_pool,1); + histogram(tis(acc(:)==1)) + axis tight; +end + +%% Metroplis 2 + + +parm.n_real=4; +var=1; + +%n_pool=48; +%parpool(n_pool); + +n_parm=10; + +T=cell(n_parm,1); +OF=cell(n_parm,1); +OFr=cell(n_parm,1); +A=cell(n_parm,1); + +n_last_iter=6; + + +for i1=1:n_parm % number of param to calibrate + + T{i1} = cell(i1,1); + OF{i1} = cell(i1,1); + OFr{i1} = cell(i1,1); + A{i1} = cell(i1,1); + + if i1==1 + T{i1}{1} = .5 ; + [OF{i1}{1}, OFr{i1}{1}] = fmin(T{i1}{1},parm,ny,nx,sn,start,nb,LAMBDA,NEIGH,S,sec,path,f0,id,kern,Gamma_t_id,Sigma_d,XY,dens,hd); + A{i1}{1} = 1 ; + disp('Initialized the first test case') + else + T{i1}{1} = T{i1-1}{end}(end); % upodate the size of T_best + disp(['Changed to new level with ' num2str(i1) ' params']) + end + + for i2 = 1:i1 % each of param to calibrate + + if i2~=1 + disp(['Changed to the ' num2str(i2) 'th parameters of the ' num2str(i1) ' from this level']) + T{i1}{i2} = T{i1}{i2-1}(end); + OF{i1}{i2} = OF{i1}{i2-1}(end); + OFr{i1}{i2} = OFr{i1}{i2-1}(end); + A{i1}{i2} = 1; + end + + condition=1; + + last=1; + cur=2; + i_var=0; + var = exp(i_var); + while condition + disp(['New test nb ' num2str(cur) ' | parameters: ' num2str(i2) '/' num2str(i1)]) + if (mod(cur,n_last_iter)==0) + dacc_ratio = mean(A{i1}{i2}(cur-n_last_iter+1:cur-1)); + if dacc_ratio > .5 % accept too much -> reduce variance + i_var = i_var-1; + elseif dacc_ratio < .1 % reject too much -> increase variance + i_var = i_var-1; % + end + var = exp(i_var); + disp(['dacc_ratio=' num2str(dacc_ratio) ' | var=' num2str(var)]) + if i_var<-5 + condition=0; + end + end + + % Perturbation of the ith param with a var var and block + % between 0 and 1 + T{i1}{i2}(cur,:) = mod(normrnd(T{i1}{i2}(last,:),var),1); + disp(['New T=' num2str(T{i1}{i2}(cur,:)) ' | Last T =' num2str(T{i1}{i2}(last,:))]) + + + [OF{i1}{i2}(cur), OFr{i1}{i2}(cur)] = fmin(T{i1}{i2}(cur,:),parm,ny,nx,sn,start,nb,LAMBDA,NEIGH,S,sec,path,f0,id,kern,Gamma_t_id,Sigma_d,XY,dens,hd); + + while abs(OF{i1}{i2}(cur) - OF{i1}{i2}(last)) < abs(OFr{i1}{i2}(cur))+abs(OFr{i1}{i2}(last)) + parm.n_real = parm.n_real*2; + disp(['Update the number of real to: ' num2str(parm.n_real)]) + if parm.n_real >500 + condition = 0; + disp('Too many real, stop the current parameter'); + continue + end + [OF{i1}{i2}(last),OFr{i1}{i2}(last)] = fmin(T{i1}{i2}(last,:),parm,ny,nx,sn,start,nb,LAMBDA,NEIGH,S,sec,path,f0,id,kern,Gamma_t_id,Sigma_d,XY,dens,hd); + [OF{i1}{i2}(cur),OFr{i1}{i2}(cur)] = fmin(T{i1}{i2}(cur,:),parm,ny,nx,sn,start,nb,LAMBDA,NEIGH,S,sec,path,f0,id,kern,Gamma_t_id,Sigma_d,XY,dens,hd); + end + + + if OF{i1}{i2}(cur) < OF{i1}{i2}(last) + last = cur; + A{i1}{i2}(cur)=1; + disp(['Accepted with OF=' num2str(OF{i1}{i2}(cur)) ' +/- ' num2str(OFr{i1}{i2}(cur))]); + else + A{i1}{i2}(cur)=0; + end + + +% if sum(A{i1}{i2}) >= n_last_iter +% y=OF{i1}{i2}(A{i1}{i2}); +% y=y(en-n_last_iter:end); +% x=(1:n_last_iter)'; +% slope = x'\y; +% ttest = ( (slope-0)*sqrt(n_last_iter-2) )/sqrt( sum((y-slope).^2) / sum( ( x-mean(x) ).^2 )); +% +% if (ttest< threasold) +% condition=0; +% end +% end + + cur=cur+1; + end + + end + +end + + + + + + + + + + + + + + + + + + + +