Field를 만들고, 특정 Field에서 원하는 값 가져오기 (interp2)

매트랩에서 시뮬레이션 하다보면 주어진 영역을 GRID로 만들어서 처리해야 할 일이 많다. 

이 때 유용하게 사용될 수 있는 것이 mesh와 grid이다. 

다음의 예를 통해서 설명해보겠다. 

 주어진 영역 (10*10)의 공간을 내가 임의의 resolution(20)으로 나눠서 처리할 일이 있다고 하자. Gaussian Process를 통한 prediction을 하거나 SVM 같은 classifier의 decision range를 눈으로 확인할 때 이러한 처리가 많이 사용될 것이다. 

 이번 포스팅에서는 주어진 영역을 적절한 grid로 나눈 후에 GP realization을 하고, 랜덤으로 설정한 위치에서 GP realization의 값을 읽은 것을 해보겠다. GP realization에 대해서는 이전 포스팅을 참조하도록 하자. 

 결과는 다음과 같다. 

data @(0.421, 0.819) => 1.621 

data @(1.509, 0.058) => 2.700 

data @(0.477, 0.733) => 3.156 

data @(4.297, 4.871) => 2.854 

data @(4.984, 2.768) => 2.577 

data @(1.653, 2.150) => 2.459 

data @(0.355, 4.439) => -0.458 

data @(0.186, 1.181) => 0.056 

data @(-0.345, -0.596) => -0.009 

data @(0.223, -0.091) => -1.519 

>> 

1. demo_gridstuff.m

%%

clc;

clear all;

close all;

%%

% 1. stage의 영역을 설정한다. 

stage.xmin = 0;

stage.xmax = 5;

stage.ymin = 0;

stage.ymax = 5;

stage.xres = 20;

stage.yres = 20;

stage.xgrid = linspace(stage.xmin, stage.xmax, stage.xres);

stage.ygrid = linspace(stage.ymin, stage.ymax, stage.yres);

[mesh_X, mesh_Y] = meshgrid(stage.xgrid, stage.ygrid); 

% 전체 grid의 (x, y) point가 col방향으로 들어가있다. 

stage.points = [mesh_X(:) mesh_Y(:)]; 

stage.nr_points = length(stage.points);

% 2. 이렇게 만든 stage에 GP realization을 통해서 값을 넣자. 

option_gp.useSE_K    = 1;

option_gp.sigma2f    = 1;

option_gp.sigma2w    = 1E-8;

option_gp.sigma2x    = 5E-1;

option_gp.Grad = 0; 

field = GP_realization(stage, option_gp);

%% 

% 3. 현재 주어진 field에서 임의의 점(x, y)의 데이터(z)를 얻어오는 방법 

nr_estimate = 10; 

x_rand = stage.xmin + (stage.xmax-stage.xmin)*rand(nr_estimate, 1);

y_rand = stage.ymin + (stage.ymax-stage.ymin)*rand(nr_estimate, 1);

z_rand = zeros(nr_estimate, 1);

for zi = 1:nr_estimate

    z_rand(zi) = interp2(stage.xgrid, stage.ygrid, field.zdata, ...

            x_rand(zi), y_rand(zi) );

end

fprintf('data @(%.3f, %.3f) \t=> %.3f \n', x_rand, y_rand, z_rand);

%%

% 4. 그림을 그려보자.  

% 먼저 Field를 그리고, 

figure();pcolor(stage.xgrid, stage.ygrid, field.zdata);

set(gca,'YDir','normal');  % 중요! <= Matrix와 Graph는 y축 방향이 반대이다. 

colorbar; caxis([field.zmin field.zmax]); hold on; 

shading interp;lighting phong; 

% 그리고 포인트와, 해당 field의 값을 그린다. 

for ei = 1:nr_estimate

    hold on;

    plot(x_rand(ei), y_rand(ei), 'o', 'LineWidth',2,...

        'MarkerEdgeColor','k',...

        'MarkerSize',10);

    str = sprintf('   %.3f', z_rand(ei));

    text(x_rand(ei), y_rand(ei), str);

    hold off;

end

title('Gaussain Process Realization and Estimates');

2. GP_realization (Gaussian Process를 통해서 realizaton하는 부분)

function [field] ...

    = GP_realization(stage, option_gp)

x_grid_resol    = stage.xres;

y_grid_resol    = stage.yres;        

x_min           = stage.xmin;

x_max           = stage.xmax;

y_min           = stage.ymin;

y_max           = stage.ymax;

%%

x       = linspace(x_min, x_max, x_grid_resol);

y       = linspace(y_min, y_max, y_grid_resol);

[X,Y]   = meshgrid(x, y); 

mesh    = [X(:) Y(:)]; % 2-D mesh

n       = size(mesh, 1);

%%

K       = GP_kernel(mesh, mesh, option_gp) + sqrt(option_gp.sigma2w)*eye(n, n);

z       = chol(K)'*randn(n,1);

Z       = reshape(z, y_grid_resol, x_grid_resol);

field       = [];

field.xgrid = X;

field.ygrid = Y;

field.zdata = Z - mean(mean(Z));

field.xres = x_grid_resol;

field.yres = y_grid_resol;

field.xmin = x_min;

field.xmax = x_max;

field.ymin = y_min;

field.ymax = y_max; 

field.zmin = min(min(field.zdata ));

field.zmax = max(max(field.zdata ));

end

%%

% pred.pos  = [reshape(field.xgrid, [], 1), ...

%     reshape(flipud(field.ygrid), [], 1)];

function K = GP_kernel(X1, X2, option)

%

% X1: 데이터는 col방향으로 들어있다.

% X2: 데이터는 col방향으로 들어있다.

%

% 여기에는 sigma2w는 포함되어있지 않다!! 

%

% <SE kernel 예제 1>

% option_gp.useSE_K    = 1;

% option_gp.sigma2f    = 1;

% option_gp.sigma2w    = 1E-3;

% option_gp.sigma2x    = 5E-1;

% option_gp.Grad = 0; <== 미분을 안한다. 할 경우는 일단 scalar

% 

% <D kernel 예제 2>

% option_gp.useD_K     = 1;

% option_gp.sigma2f    = 1;

% option_gp.sigma2w    = 1E-3;

% option_gp.r          = 0.5;

% option_gp.Grad = 0;

%

% 

%

if ~isfield(option, 'Grad')

    option.Grad = 0;

end

if ~isfield(option, 'useSE_K')

    option.useSE_K = 0;

end

if ~isfield(option, 'useD_K')

    option.useD_K = 0;

end

len_X1  = size(X1, 1);

len_X2  = size(X2, 1);

K       = zeros(len_X1, len_X2);

for i = 1:len_X1

    for j  = 1:len_X2

        

        if option.Grad

            % X1에서 미분한 gradient를 구한다. (이 때는 앞에서 정의한 K를 사용하지 않을 것이다.)

            % 미분을 일반적으로 하나의 X1과 X2를 가정한다. (예: 두 위치, len_X1 = len_X2 = 1)

            if option.useSE_K

                K = se_kernel(X1(i,:), X2(j,:), option)*-(X1(i,:)-X2(i,:))/option.sigma2x;

            elseif option.useD_K

                diff  = X1(i,:) - X2(j,:);

                dist  = norm(diff);

                sig2f = option.sigma2f;

                r     = option.r;

                if dist == 0

                    K   = zeros(size(diff));

                elseif dist > r

                    K   = zeros(size(diff));

                else

                    K   = sig2f*pi*(dist/r-1)*sin(pi*dist/r)*diff/dist/r;                    

                end

            end

        else

            % kernel의 값을 구한다. 

            if option.useSE_K

                K_data = se_kernel(X1(i,:), X2(j,:), option);

            elseif option.useD_K

                K_data = distribute_kernel(X1(i,:), X2(j,:), option);

            end

            K(i, j) = K_data; 

        end

    end

end

end

%%

function output = se_kernel(x1, x2, option)

    output = option.sigma2f*exp(-norm(x1-x2)^2/option.sigma2x/2);

end

%%

function output = distribute_kernel(x1, x2, option)

    output = option.sigma2f*Lambda(norm(x1-x2)/option.r);

end

function output_lambda = Lambda(h)

    if h <= 1 

        output_lambda = (1-h)*cos(pi*h) + 1/pi*sin(pi*h);

    else 

        output_lambda = 0;

    end

end