Cross hole GPR tomopgraphy using Neural Networks
================================================

The example below requires the `SIPPI
toolbox <https://sippi.readthedocs.io/>`__.

::

    % grl_nn 
    clear all;close all

    createTrainingSet=1;
    createReferenceModel=1;
    useTargetDist=0;
    Ntrain=40000;

    Nr_modeling=6000; %size of sample for modeling error

    TrainSizes=[1000 5000 10000 20000 40000];

    splitData=3;
    epochs=30000; 
    hiddenLayerSize=80;

    %
    %Ntrain=2500;
    %splitData=0;
    %epochs=10000; 
    %hiddenLayerSize=20;
    %TrainSizes=[1000 2500];


    % mcmc 
    nite=2000000;;

    %% LOAD DATA AND CONFIGURATION
    D=load('AM13_data.mat');

    % REVERSE S and R for fd forward
    D2=D;
    D.S(352:end,:)=D2.R(352:end,:);
    D.R(352:end,:)=D2.S(352:end,:);
    D.S(:,1)=D.S(:,1)+1;
    D.R(:,1)=D.R(:,1)+1;
    clear D2
    i_use=1:1:size(D.S,1);
    id=1;
    data{id}.d_obs=D.d_obs(i_use);
    data{id}.d_std=D.d_std(i_use).*0+0.1;


    %% B: Define FD
    %forward_fd;
    forward.forward_function='sippi_forward_traveltime';
    forward.sources=D.S(i_use,:);
    forward.receivers=D.R(i_use,:);
    forward.type='fd';
    %forward.type='fat';forward.linear=1;forward.freq=0.1;
    %forward.type='eikonal';
    %forward.m0=prior{1}.m0;

    forward.fd.t=1e-7;
    forward.fd.addpar.Tg=100*10^6;
    forward.fd.dx_fwd=0.1;

    %% A: Define prior model
    im=1;
    dx=0.2;
    prior{im}.type='FFTMA';
    prior{im}.name='Velocity (m/ns)';
    prior{im}.m0=0.1431;
    prior{im}.Va='.000215 Sph(6)';
    prior{im}.x=[-.5:dx:6.5];
    prior{im}.x=[0:dx:7];
    prior{im}.y=[0:dx:13];
    prior{im}.cax=[-1 1].*.03+prior{im}.m0;

    if useTargetDist==1;
        d_target=[randn(1,100)*.003-0.01 randn(1,100)*.003+0.01]+prior{im}.m0;
        prior{im}.d_target=d_target;
        prior{im}.m0=0; %% MAKE SURE sippi_forward_traveltime tests for a non-zero velocity
    end


    %% C: Make Reference model
    if createReferenceModel==1
        rng(1);
        % Reference mo
        [m_ref,prior]=sippi_prior(prior);
        NM=prod(size(m_ref{1}));
        % reference data
        [d_ref,forward]=sippi_forward(m_ref,forward,prior);

        % data    
        data{1}.d_ref=d_ref{1};
        data{1}.d_noise=randn(size(d_ref{1})).*data{1}.d_std;
        data{1}.d_obs=data{1}.d_ref+data{1}.d_noise;



        % compute reference data in m0
        m_ref0=m_ref;
        m_ref0{1}=m_ref0{1}.*0+prior{1}.m0;
        disp('computing reference forward')
        [d_ref0,forward0]=sippi_forward(m_ref0,forward,prior);
        d0=d_ref0{1};

        save grl_ReferenceModel 
    else
        load grl_ReferenceModel
    end

    %% D: Create M-D training data set for forward model
    if createTrainingSet==1
        ATTS=zeros(length(m_ref{1}(:)),Ntrain);
        DATA=zeros(length(d_ref{1}(:)),Ntrain);

        iplot=1;
        t0=now;
        for i=1:Ntrain;

            if ((i/iplot)==round(i/iplot)&&(i>1));progress_txt(i,Ntrain,time_loop_end(t0,i-1,Ntrain));end
            m=sippi_prior(prior);
            try
                d=sippi_forward(m,forward,prior,data);

                ATTS(:,i)=m{1}(:);
                DATA(:,i)=d{1}(:)-d0;
            catch
                disp(sprintf('Something went wrong.'));
                i=i-1;
            end
        end
        save(sprintf('grl_%s_NM%d_NT%d',forward.type,NM,Ntrain));
    else
        load grl_eikonal_NM2376_NT300
        %load grl_eikonal_NM2376_NT1000.mat
    end


    %% E: SETUP FORWARD MODELS

    if ~exist('splitData'); 
        splitData=0; % SPLIT DATA FOR NN
    end
    if ~exist('epochs'); 
        epochs=100000; %   
    end
    if ~exist('hiddenLayerSize'); 
        hiddenLayerSize=80;
    end

    forward_nn.forward_function='sippi_forward_mynn';
    forward_nn.sources=forward.sources;
    forward_nn.receivers=forward.receivers;
    forward_nn.ATTS=ATTS;
    forward_nn.DATA=DATA;
    forward_nn.d0=d0;
    clear DATA ATTS
    forward_nn.splitData=splitData;
    forward_nn.epochs=epochs;;  
    forward_nn.hiddenLayerSize=hiddenLayerSize;
    forward_nn.max_nm=1e+10;
    txt=sprintf('grl_NM%d_DX%d_%s_NT%d_SD%d_NH%d',NM,dx*100,forward.type,Ntrain,forward_nn.splitData,forward_nn.hiddenLayerSize);
    forward_nn.mfunc_string=txt;
    disp(sprintf('%s: setting name ''%s''',mfilename,txt));

    % setup all forward models
    i_forward=0;
    if ~exist('TrainSizes'); 
        TrainSizes=[100 200];
    end
    for n_use=TrainSizes;
        if n_use<= size(forward_nn.ATTS,2);
            i_forward=i_forward+1;
            f_mul{i_forward}=forward_nn;
            f_mul{i_forward}.ATTS=forward_nn.ATTS(:,1:n_use);
            f_mul{i_forward}.DATA=forward_nn.DATA(:,1:n_use);       
            txt_use=sprintf('grl_NM%d_DX%d_%s_NT%d_SD%d_NH%d',NM,dx*100,forward.type,n_use,forward_nn.splitData,forward_nn.hiddenLayerSize);
            f_mul{i_forward}.mfunc_string=txt_use;
        end
    end

    % eikonal
    i_forward=i_forward+1;
    f_mul{i_forward}=forward;
    f_mul{i_forward}.type='eikonal';

    % ray_2d
    i_forward=i_forward+1;
    f_mul{i_forward}=forward;
    f_mul{i_forward}.type='ray_2d';
    f_mul{i_forward}.linear=1;
    f_mul{i_forward}.freq=0.1;
    f_mul{i_forward}.r=2;
    f_mul{i_forward}.normalize_vertical=0;


    %% F: EVALUATE forward models once to setup NN and Linear operators
    for i=1:length(f_mul);
        t1=now;
        [d_mul{i},f_mul{i}]=sippi_forward(m_ref,f_mul{i},prior);
        t2=now;
        time_mul{i}=(t2-t1)*3600*24;

    end

    save(sprintf('%s_forward',txt))

    %% G: Estimate modeling errors
    if ~exist('Nr_modeling'); 
        Nr_modeling=6000;
    end
    [Ct,dt,dd,d_full,d_app]=sippi_compute_modelization_forward_error(forward,f_mul,prior,Nr_modeling);

    % H: Setup one data structure per forward model, with the correct modeling error
    for i=1:length(f_mul);    

        %s=sum(abs(dd{5}));
        %ii=find(s<180);

        data_mul{i}=data;
        data_mul{i}{1}.dt=dt{i};
        data_mul{i}{1}.Ct=Ct{i};
    end

    save(sprintf('%s_modelerr',txt))


    %% I: Perform probabilistic inversion using extended Metropolis sampling
    if ~exist('nite'); 
        nite=1000000; % 
    end
    options.mcmc.m_ref=m_ref;
    options.mcmc.nite=nite;   % [1] : Number if iterations
    options.mcmc.i_sample=ceil(options.mcmc.nite/1000); % : Number of iterations between saving model to disk
    options.mcmc.i_plot=100000;  % [1]: Number of iterations between updating plots
    options.mcmc.i_save_workspace=100000;  % [1]: Number of iterations between

    i_burnin=options.mcmc.nite/30;
    prior{1}.seq_gibbs.i_update_step_max=i_burnin;

    options.mcmc.anneal.i_begin=1; % default, iteration number when annealing begins
    options.mcmc.anneal.i_end=ceil(i_burnin/2); %  iteration number when annealing stops
    options.mcmc.anneal.T_begin=5; % Start temperature for annealing
    options.mcmc.anneal.T_end=1; % End temperature for annealing

    %options.mcmc.n_chains=2; % set number of chains (def=1)
    %options.mcmc.T=[1 1.1 1.2];      % set temperature of chains [1:n_chains]

    % RUN MCMC

    rseed=1;
    for i=1:(length(f_mul));
        rng(rseed);
        if isfield(f_mul{i},'mfunc');
            options.txt=f_mul{i}.mfunc_string;        
        else
            try
                %options.txt=sprintf('grl_NM%d_DX%d_%s_SD%d_NH%d',NM,dx*100,f_mul{i}.type,forward_nn.splitData,forward_nn.hiddenLayerSize);
                options.txt=sprintf('grl_NM%d_DX%d_%s_SD%d_NH%d',NM,dx*100,forward_nn.type,forward_nn.splitData,forward_nn.hiddenLayerSize);
            catch
                options.txt=txt;
            end
        end
        t_start=now;
        [o]=sippi_metropolis(data_mul{i},prior,f_mul{i},options);
        options_out{i}.txt=o.txt;
        %sippi_plot_posterior_sample(options_out{i}.txt);
        sim_minutes(i)=(now-t_start)*60*24;
    end

    save(sprintf('%s_inverted',txt))