text:
function[]=makefigs_morsearea
%MAKEFIGS_MORSEAREA  Makes a sample figure for MORSEAREA.

C1=logspace(0.1,3,100)';
be1=linspace(1,10,101)';
[be,C]=meshgrid(be1,C1);
ga=2;

figure
A = morsearea(C,ga(1),be);
contourf(C1,be1,log10(A'),20)
xlog
xlabel('Parameter C')
ylabel('Parameter \beta')
title('Log10 Area with \gamma =2')
colorbar

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
classdef task_pf < mp.task
%MP.TASK_PF  MATPOWER task for power flow (PF).
%   MP.TASK_PF provides implementation for power flow problem.
%
%   Properties
%       dc
%       iterations
%       ref
%       ref0
%       va_ref0
%       fixed_q_idx
%       fixed_q_qty
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2020-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        tag = 'PF';
        name = 'Power Flow';
        dc      %% true if DC network model (cached in run_pre(), from mpopt)
        iterations
        ref
        ref0
        va_ref0
        fixed_q_idx
        fixed_q_qty
    end

    methods
        %%-----  task methods  -----
        function [d, mpopt] = run_pre(obj, d, mpopt)
            [d, mpopt] = run_pre@mp.task(obj, d, mpopt);    %% call parent

            %% cache DC model flag
            obj.dc = strcmp(upper(mpopt.model), 'DC');
        end

        function dm = next_dm(obj, mm, nm, dm, mpopt, mpx)
            if ~obj.dc && mpopt.pf.enforce_q_lims
                %% adjust iteration count for previous runs
                obj.iterations = obj.iterations + mm.soln.output.iterations;
                mm.soln.output.iterations = obj.iterations;

                %% enforce Q limits
                [success, dm] = obj.enforce_q_lims(nm, dm, mpopt);
                if ~success                 %% entire task fails if Q lim
                    obj.success = success;  %% enforcement indicates failure
                end
            else        %% don't enforce generator Q limits, once is enough
                dm = [];
            end
        end

        function [success, dm] = enforce_q_lims(obj, nm, dm, mpopt);
            gen_dme = dm.elements.gen;
            [mn, mx, both] = gen_dme.violated_q_lims(dm, mpopt);

            if ~isempty(both)   %% we have some Q limit violations
                if isempty(mn) && isempty(mx)   %% infeasible
                    if mpopt.verbose
                        fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(both));
                    end
                    dm = [];
                    success = 0;
                else
                    if mpopt.verbose && ~isempty(mx)
                        fprintf('Gen %d at upper Q limit, converting to PQ bus\n', gen_dme.on(mx));
                    end
                    if mpopt.verbose && ~isempty(mn)
                        fprintf('Gen %d at lower Q limit, converting to PQ bus\n', gen_dme.on(mn));
                    end

                    %% save corresponding limit values
                    obj.fixed_q_qty(mx) = gen_dme.qg_ub(mx);
                    obj.fixed_q_qty(mn) = gen_dme.qg_lb(mn);
                    mx = [mx;mn];

                    %% set qg to binding limit
                    gen_dme.tab.qg(gen_dme.on(mx)) = ...
                        obj.fixed_q_qty(mx) * dm.base_mva;

                    %% convert to PQ bus
                    bus_dme = dm.elements.bus;
                    ref0 = find(bus_dme.type == mp.NODE_TYPE.REF);
                    bidx = bus_dme.i2on(gen_dme.bus(gen_dme.on(mx)));   %% bus of mx
                    if length(ref0) > 1 && any(bus_dme.type(bidx) == mp.NODE_TYPE.REF)
                        error('mp.task_pf/enforce_q_lims: Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.');
                    end
                    %% set bus type to PQ
                    bus_dme.set_bus_type_pq(dm, bidx);
                    %% potentially pick new slack bus
                    ntv = nm.node_types(nm, dm);        %% node type vector
                    [i1, iN] = nm.get_node_idx('bus');  %% bus node indices
                    btv = ntv(i1:iN);                   %% bus type vector

                    %% indicate if there's been a change in slack bus
                    ref = find(btv == mp.NODE_TYPE.REF);    %% new ref bus indices
                    if mpopt.verbose && ref ~= ref0
                        fprintf('Bus %d is new slack bus\n', ...
                            bus_dme.ID(bus_dme.on(ref)));
                    end

                    %% save indices to list of Q limited gens
                    obj.fixed_q_idx = [obj.fixed_q_idx; mx];

                    %% update dm for next step
                    dm.initialize();
                    dm.update_status();
                    dm.build_params();
                    success = 1;
                end
            else                %% no more Q violations
                dm = [];
                success = 1;
            end
        end

        %%-----  data model methods  -----

        %%-----  network model methods  -----
        function nm_class = network_model_class_default(obj, dm, mpopt)
            switch upper(mpopt.model)
                case 'AC'
                    if mpopt.pf.v_cartesian
                        nm_class = @mp.net_model_acc;
                    else
                        nm_class = @mp.net_model_acp;
                    end
                case 'DC'
                    nm_class = @mp.net_model_dc;
            end
        end

        function nm = network_model_build_post(obj, nm, dm, mpopt)
            %% initialize task data, if non-empty AC case with Q lim enforced
            if ~obj.dc && mpopt.pf.enforce_q_lims ~= 0 && nm.np ~= 0
                if obj.i_nm == 1
                    [ref, ~, ~] = nm.node_types(obj, dm);
                    gen_dme =  dm.elements.gen;
                    obj.iterations = 0;
                    obj.ref0 = ref;             %% initial ref node indices
                    obj.ref = ref;              %% current ref node indices
                    obj.va_ref0 = nm.get_va(ref);%% initial ref node voltage angles
                    obj.fixed_q_idx = [];       %% indices of fixed Q gens
                    obj.fixed_q_qty = zeros(gen_dme.n, 1);  %% Q output of fixed Q gens
                else        %% update index of ref bus
                    [obj.ref, ~, ~] = nm.node_types(obj, dm);
                end
            end
        end

        function nm = network_model_x_soln(obj, mm, nm)
            nm = network_model_x_soln@mp.task(obj, mm, nm);

            %% if ref node has been changed, adjust voltage angles
            %% to make angle at original ref node = specified value
            if ~obj.dc && obj.i_nm > 1 && obj.ref ~= obj.ref0
                vm = abs(nm.soln.v);
                va = angle(nm.soln.v);
                va = va - va(obj.ref0) + obj.va_ref0;
                nm.soln.v = vm .* exp(1j * va);
            end
        end

        %%-----  mathematical model methods  -----
        function mm_class = math_model_class_default(obj, nm, dm, mpopt)
            switch upper(mpopt.model)
                case 'AC'
                    if mpopt.pf.v_cartesian
                        if mpopt.pf.current_balance
                            mm_class = @mp.math_model_pf_acci;
                        else
                            mm_class = @mp.math_model_pf_accs;
                        end
                    else
                        if mpopt.pf.current_balance
                            mm_class = @mp.math_model_pf_acpi;
                        else
                            mm_class = @mp.math_model_pf_acps;
                        end
                    end
                case 'DC'
                    mm_class = @mp.math_model_pf_dc;
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function d = gamkl(x,logx,apost,bpost,aprior,bprior)
warning('This function is deprecated')

% KL divergence for Gamma distributions.
% A lower bound for a Gamma distributed X:
% p(X) = G(aprior,bprior)
% q(X) = G(apost,bpost)
% <X> = x
% <log X> = logx
%
% KL(q||p) = <log q(X)> - <log p(X)> = -H(q) - <log p(x)>
% where expectation is taken over q(X) and H() is entropy

% Cost from prior
d = -aprior*log(bprior) + gammaln(aprior) - (aprior-1)*logx + bprior*x;
% Cost from q-posterior
d = d - gament(apost,bpost);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% A = riemann_mean(B,epsilon,tol)
%
% Calcul du barycentre des matrice de covariances.
% A : baricentre des K matrices NxN
%
% B : Matrice NxNxK
% epsilon : Pas de la descente de gradient
% tol : arret de la descente si le crit�re < tol

function [A, critere, niter] = riemann_mean(B,args)

N_itermax = 100;
if (nargin<2)||(isempty(args))
    tol = 10^-5;
    A = mean(B,3);
else
    tol = args{1};
    A = args{2};
end

niter = 0;
fc = 0;

while (niter<N_itermax)
    niter = niter+1;
    % Tangent space mapping
    T = Tangent_space(B,A);
    % sum of the squared distance
    fcn = sum(sum(T.^2));
    % improvement
    conv = abs((fcn-fc)/fc);
    if conv<tol % break if the improvement is below the tolerance
       break;
    end
    % arithmetic mean in tangent space
    TA = mean(T,2);
    % back to the manifold
    A = UnTangent_space(TA,A);
    fc = fcn;
end

if niter==N_itermax
    disp('Warning : Nombre d''it�rations maximum atteint');
end

critere = fc;

author: alexandrebarachant

repo: covariancetoolbox


text:
function [equivalence_classes, counter, newA, newl] = palette_wl(A, labels)
%  Usage: palette_wl for enclosing subgraph labeling in WLNM
%  --Input--
%  -A: original adjacency matrix of the enclosing subgraph
%  -labels: initial colors
%  --Output--
%  -equivalence_classes: final labels
%
%  Original author: Roman Garnett
%  Original paper:
%    Kersting, K., Mladenov, M., Garnett, R., and Grohe, M. Power
%    Iterated Color Refinement. (2014). AAAI Conference on Artificial
%    Intelligence (AAAI 2014).
%
%  *author: Muhan Zhang, Washington University in St. Louis
%%
  original_labels = labels;
  % if no labels given, use initially equal labels
  if (nargin < 2)
    labels = ones(size(A, 1), 1);
  end

  equivalence_classes = zeros(size(labels));

  % iterate WL transformation until stability
  counter = 0;
  while (~isequal(labels, equivalence_classes))
    counter = counter + 1;
    equivalence_classes = labels;
    labels = wl_transformation(A, labels);
  end

  [~, I] = sort(equivalence_classes, 'descend');
  newA = A(I, I);
  newl = original_labels(I);

end

function new_labels = wl_transformation(A, labels)

  % the ith entry is equal to the 2^(i-1)'th prime
  primes_arguments_required = [2, 3, 7, 19, 53, 131, 311, 719, 1619, ...
                      3671, 8161, 17863, 38873, 84017, 180503, 386093, ...
                      821641, 1742537, 3681131, 7754077, 16290047];

  num_labels = max(labels);

  % generate enough primes to have one for each label in the graph
  primes_argument = primes_arguments_required(ceil(log2(num_labels)) + 1);
  p = primes(primes_argument);
  log_primes = log(p(1:num_labels))';

  Z = ceil(sum(log_primes(labels)));

  signatures = labels + A * log_primes(labels) / Z;

  % map signatures to integers counting from 1
  [~, ~, new_labels] = uniquetol(signatures);

end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function grains = cat(dim,varargin)
% implements [grain1, grain2]
%
% Syntax
%   g = [grains_1, grains_2, grains_n]
%   g = [grains('fe') grains('mg')]
%   g = [grains(1:100) grains(500:end)]
%
% See also
% grain2d/vertcat

grains = cat@dynProp(1,varargin{:});

for k = 2:numel(varargin)

  ng = varargin{k};

  if isempty(ng), continue; end
  grains.id = [grains.id; ng.id];
  grains.phaseId = [grains.phaseId; ng.phaseId];
  grains.grainSize = [grains.grainSize; ng.grainSize];
  grains.poly = [grains.poly; ng.poly];
  grains.inclusionId = [grains.inclusionId; ng.inclusionId];
  grains.boundary = [grains.boundary; ng.boundary];
  grains.innerBoundary = [grains.innerBoundary; ng.innerBoundary];

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function incube()
    % This subroutine assigns Parameter values for the grid
    % which is used to calculate a Max Z value map.
    % dx, dy is the grid spacing in degrees.
    % For each one of the grid points, Ni events are counted.
    %

    global step
    report_this_filefun();

    % initial values
    %
    dx = 1.00;
    dy = 1.00 ;
    ni = 100;

    %
    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Alarm Cube Input Parameter',...
        'NumberTitle','off', ...
        'NextPlot','new', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ ZG.welcome_pos + [200, -200], 450, 250]);
    axis off

    % creates a dialog box to input grid parameters
    %
    freq_field  =uicontrol('Style','edit',...
        'Position',[.60 .50 .22 .10],...
        'Units','normalized','String',num2str(step),...
        'callback',@callbackfun_001);

    freq_field2=uicontrol('Style','edit',...
        'Position',[.60 .40 .22 .10],...
        'Units','normalized','String',num2str(dx),...
        'callback',@callbackfun_002);

    uicontrol('Units','normal',...
        'Position',[.1 .70 .20 .12] ,'String','Load Alarm',...
        'callback',@callbackfun_003)

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.80 .05 .15 .12 ],...
        'Units','normalized', 'Callback', @callbackfun_004,'String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.10 .05 .35 .12 ],...
        'Units','normalized',...
        'callback',@callbackfun_005,...
        'String','Rubberband');

    go_button2=uicontrol('Style','Pushbutton',...
        'Position',[.50 .05 .15 .12 ],...
        'Units','normalized',...
        'callback',@callbackfun_006,...
        'String','LTA');

    txt3 = text(...
        'Position',[0.30 0.84 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'Color','r' ,...
        'FontWeight','bold',...
        'String',' Alarm Cube Parameters ');

    txt5 = text(...
        'Position',[0. 0.42 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Time length of winlen_days in years');

    txt6 = text(...
        'Position',[0. 0.20 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Select one of the statistics below ');

    txt1 = text(...
        'Position',[0. 0.53 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m,...
        'FontWeight','bold',...
        'String','Step width (Integer): ');

    set(gcf,'visible','on');
    watchoff

    function allhist(sta)
        % ZMAP script allhist.m.
        % calculates a histogram of all z-value in space and time
        % using the rubberband function with a window length of 1.5 years
        %
        % stefan wiemer 11/94
        %
        % sta : either rub (rubber band) or lta

        global iala

        % Input Rubberband
        %
        report_this_filefun();

        tre2 = 4.00;

        clear abo;
        abo=[];

        % initial parameter
        winlen_days = floor(ZG.compare_window_dur/ZG.bin_dur);
        [len, ncu] = size(cumuall); len = len-2;
        var1 = zeros(1,ncu);
        var2 = zeros(1,ncu);
        mean1 = zeros(1,ncu);
        mean2 = zeros(1,ncu);
        as = zeros(1,ncu);
        n2 = [];

        % loop over all grid points for percent
        %
        %

        % loop over all point for rubber band
        %
        wai = waitbar(0,' Please Wait ...  ');
        set(wai,'NumberTitle','off','Name','Allhist - Percent done');
        pause(1)

        drawnow

        n2 = zeros(1,length(-15:0.1:15));
        if sta == 'lta'

            for ti = winlen_days:step:len-winlen_days
                cu = [cumuall(1:ti-1,:) ; cumuall(ti+winlen_days+1:len,:)];
                mean1 = mean(cu);
                mean2 = mean(cumuall(ti:ti+winlen_days,:));
                for i = 1:ncu
                    var1(i) = cov(cu(:,i));
                    var2(i) = cov(cumuall(ti:ti+winlen_days,i));
                end     % for i
                as = (mean1 - mean2)./(sqrt(var1/(len-winlen_days)+var2/winlen_days));
                [m,n] = size(as);
                reall = reshape(as,1,m*n);

                % set values gretaer ZG.tresh_km = nan
                %
                %s = cumuall(len,:);
                %r = reshape(as,length(gy),length(gx));
                l = reall > tre2;
                s = [  loc(1,l) loc(2,l) loc(3,l)   reall(l) ];
                s = [reshape(s,length(s)/4,4) ones(length(s)/4,1)*ti];
                abo = [abo ;  s];
                [n,x] =histogram(reall,(-15:0.10:15));
                n2 = n2 + n;
                waitbar((ti-winlen_days)/(len-2*winlen_days))
            end   % for ti
        end % if lta

        if sta == 'rub'
            for ti = winlen_days:step:len-winlen_days
                for i = 1:ncu
                    mean1(i) = mean(cumuall(1:ti,i));
                    mean2(i) = mean(cumuall(ti+1:ti+winlen_days,i));
                    var1(i) = cov(cumuall(1:ti,i));
                    var2(i) = cov(cumuall(ti+1:ti+winlen_days,i));
                end %  for i ;
                as = (mean1 - mean2)./(sqrt(var1/ti+var2/winlen_days));

                [m,n] = size(as);
                reall = reshape(as,1,m*n);

                % set values gretaer ZG.tresh_km = nan
                %
                s = cumuall(len,:);
                %r = reshape(s,length(gy),length(gx));
                l = reall > tre2;
                s = [  loc(1,l) loc(2,l) loc(3,l)   reall(l) ];
                s = [reshape(s,length(s)/4,4) ones(length(s)/4,1)*ti];
                abo = [abo ;  s];

                [n,x] =histogram(reall,(-15:0.10:15));
                n2 = n2 + n;
                waitbar((ti-winlen_days)/(len-2*winlen_days))
            end   % for ti
        end   % if riub

        close(wai)
        figure
        bar(x,n2,'k');
        grid
        xlabel('z-value ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Number ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        watchoff

        set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.l,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on')

        set(gca,'Color',[1 1 0.7])

        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Print ',...
            'callback',@callbackfun_007)

        uicontrol('Units','normal',...
            'Position',[.0 .75 .08 .06],'String','Close ',...
            'callback',@callbackfun_008)

        abo2 = abo;
        iala = ZG.compare_window_dur;
        try
            msg.infodisp('  ','Save Alarm Cube?');
            [file1,path1] = uiputfile(fullfile(ZmapGlobal.Data.Directories.data, '*.mat'), 'Save Alarm Cube Datafile?',400,400);
            if length(file1) > 1
                save([path1 file1],'cumuall','abo','loc','abo2','iala','winlen_years' )
            end
        catch ME
            warning(ME)
        end

        % plot the cube
        plotala()
    end

    function loadala()
        % Load the alarm data set
        % TODO Define what is an alarm datafile
        report_this_filefun();

        cupa = cd;

        try
            delete(pd)
        catch ME
            error_handler(ME, ' ');
        end

        [file1,path1] = uigetfile(['*.mat'],'Alarm Data File?');

        if length(path1) > 1
            load([path1 file1])
            plotala()
        else
            return
        end

    end
end
function callbackfun_001(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    step=str2double(freq_field.String);
    freq_field.String=num2str(step);
end

function callbackfun_002(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    ZG.compare_window_dur=years(str2double(mysrc.String));
end

function callbackfun_003(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    close;
    loadala;
end

function callbackfun_004(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    close;

end

function callbackfun_005(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    close;
    allhist('rub');
end

function callbackfun_006(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    close;
    allhist('lta');
end

function callbackfun_007(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    printdlg;
end

function callbackfun_008(mysrc,myevt)

    callback_tracker(mysrc,myevt,mfilename('fullpath'));
    f1=gcf;
    f2=gpf;
    set(f1,'Visible','off');
end

author: CelsoReyes


text:
function determ = ijfact1_determinant ( n )

%*****************************************************************************80
%
%% IJFACT1_DETERMINANT returns the determinant of the IJFACT1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  for i = 1 : n - 1
    determ = determ * r8_factorial ( i + 1 ) * r8_factorial ( n - i );
  end

  determ = determ * r8_factorial ( n + 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function imo = cnn_imagenet_get_batch(images, varargin)
% CNN_IMAGENET_GET_BATCH  Load, preprocess, and pack images for CNN evaluation

opts.imageSize = [227, 227] ;
opts.border = [29, 29] ;
opts.keepAspect = true ;
opts.numAugments = 1 ;
opts.transformation = 'none';
opts.affine = false;
opts.averageImage = [] ;
opts.rgbVariance = zeros(0,3,'single') ;
opts.interpolation = 'bilinear' ;
opts.numThreads = 1 ;
opts.prefetch = false ;
opts = vl_argparse(opts, varargin);

% fetch is true if images is a list of filenames (instead of
% a cell array of images)
fetch = numel(images) >= 1 && ischar(images{1}) ;

% prefetch is used to load images in a separate thread
prefetch = fetch & opts.prefetch ;

if prefetch
  vl_imreadjpeg(images, 'numThreads', opts.numThreads, 'prefetch') ;
  imo = [] ;
  return ;
end
if fetch
  im = vl_imreadjpeg(images,'numThreads', opts.numThreads) ;
else
  im = images ;
end

tfs = [] ;
switch opts.transformation
  case 'none'
    tfs = [
      .5 ;
      .5 ;
       0 ] ;
  case 'f5'
    tfs = [...
      .5 0 0 1 1 .5 0 0 1 1 ;
      .5 0 1 0 1 .5 0 1 0 1 ;
       0 0 0 0 0  1 1 1 1 1] ;
  case 'f25'
    [tx,ty] = meshgrid(linspace(0,1,5)) ;
    tfs = [tx(:)' ; ty(:)' ; zeros(1,numel(tx))] ;
    tfs_ = tfs ;
    tfs_(3,:) = 1 ;
    tfs = [tfs,tfs_] ;
  case 'stretch'
  otherwise
    error('Uknown transformations %s', opts.transformation) ;
end
[~,transformations] = sort(rand(size(tfs,2), numel(images)), 1) ;

if ~isempty(opts.rgbVariance) && isempty(opts.averageImage)
  opts.averageImage = zeros(1,1,3) ;
end
if numel(opts.averageImage) == 3
  opts.averageImage = reshape(opts.averageImage, 1,1,3) ;
end

imo = zeros(opts.imageSize(1), opts.imageSize(2), 3, ...
            numel(images)*opts.numAugments, 'single') ;

si = 1 ;
for i=1:numel(images)

  % acquire image
  if isempty(im{i})
    imt = imread(images{i}) ;
    imt = single(imt) ; % faster than im2single (and multiplies by 255)
  else
    imt = im{i} ;
  end
  if size(imt,3) == 1
    imt = cat(3, imt, imt, imt) ;
  end

  % resize
  w = size(imt,2) ;
  h = size(imt,1) ;
  factor = [(opts.imageSize(1)+opts.border(1))/h ...
            (opts.imageSize(2)+opts.border(2))/w];

  if opts.keepAspect
    factor = max(factor) ;
  end
  if any(abs(factor - 1) > 0.0001)
    imt = imresize(imt, ...
                   'scale', factor, ...
                   'method', opts.interpolation) ;
  end

  if opts.affine && ~strcmp(opts.transformation, 'none'),
      sz0 = size(imt);
      % affine
      fac = rand * .3 + 1;
      shear = rand * 0.8 - 0.4;
      aspect = rand * .5 - .25 + 1;
      theta = randi(20)-10;

      if aspect > 1 ,
          A = [fac, 0, 0; fac*shear*aspect, fac*aspect, 0; 0, 0, 1];
      else
          A = [fac/aspect, 0, 0; fac*shear, fac, 0; 0, 0, 1];
      end
      R = [cosd(theta), -sind(theta), 0; sind(theta), cosd(theta), 0; 0,0,1];
      A = R*A;
      T = affine2d(A);
      imt = imwarp(imt, T);
      w0 = size(imt,2);
      h0 = size(imt,1);
      if h0 > sz0(1) ,
          imt = imt(ceil((h0-sz0(1))/2):sz0(1)+ceil((h0-sz0(1))/2)-1,...
              ceil((w0-sz0(2))/2):sz0(2)+ceil((w0-sz0(2))/2)-1,:);
      elseif h0 < sz0(1) ,
          imt = zeros(sz0(1), sz0(2), 3);
          imt(ceil((-h0+sz0(1))/2):h0+ceil((-h0+sz0(1))/2)-1,...
              ceil((-w0+sz0(2))/2):w0+ceil((-w0+sz0(2))/2)-1,:) = imt;
      end
  end

  % crop & flip
  w = size(imt,2) ;
  h = size(imt,1) ;
  for ai = 1:opts.numAugments
    switch opts.transformation
      case 'stretch'
        sz = round(min(opts.imageSize(1:2)' .* (1-0.1+0.2*rand(2,1)), [h;w])) ;
        dx = randi(w - sz(2) + 1, 1) ;
        dy = randi(h - sz(1) + 1, 1) ;
        flip = rand > 0.5 ;
      otherwise
        tf = tfs(:, transformations(mod(ai-1, numel(transformations)) + 1)) ;
        sz = opts.imageSize(1:2) ;
        dx = floor((w - sz(2)) * tf(2)) + 1 ;
        dy = floor((h - sz(1)) * tf(1)) + 1 ;
        flip = tf(3) ;
    end
    sx = round(linspace(dx, sz(2)+dx-1, opts.imageSize(2))) ;
    sy = round(linspace(dy, sz(1)+dy-1, opts.imageSize(1))) ;
    if flip, sx = fliplr(sx) ; end

    if ~isempty(opts.averageImage)
      offset = opts.averageImage ;
      if ~isempty(opts.rgbVariance)
        offset = bsxfun(@plus, offset, reshape(opts.rgbVariance * randn(3,1), 1,1,3)) ;
      end
      imo(:,:,:,si) = bsxfun(@minus, imt(sy,sx,:), offset) ;
    else
      imo(:,:,:,si) = imt(sy,sx,:) ;
    end
    si = si + 1 ;
  end
end

repo: ResNet-Matconvnet


text:
function [vp f] = vanish_from_minevidence(lines, imgsize)
%
% find three orthogonal vanishing points by
% 1. findinhg vertical vp.
% 2. selecting one additional line.
% 3. pick the best f for the current triplet of vanishing points.
% 4. go to 2 and pick the best additional line.

%%
THRES_THETA = 10;
THRES_THETA_GUESS = 20;

%%
for i = 1:length(lines)
    lines(i).length = norm(lines(i).point1 - lines(i).point2);
end
global maxlinelength;
maxlinelength = max([lines.length]);

%% find vertical vanishing point
guessvp = [0 -1e5]; % initial guess of the vertical vp
lc1_guess = line_belongto_vp(lines, guessvp, THRES_THETA_GUESS);

vp1 = find_single_vp(lines(lc1_guess), THRES_THETA);
lc1 = line_belongto_vp(lines, vp1, THRES_THETA);

% global img; disp_lines(img, lines);
% global img; disp_lines(img, lines(lc1_guess));
% global img; disp_lines(img, lines(lc1)); plot(vp1(1), vp1(2), 'x', 'LineWidth', 5, 'MarkerSize', 20);

%% find the other two vanishing points by
[vp2 vp3 f] = find_other_vp(lines(~lc1), vp1, THRES_THETA, imgsize);

%%
INFDIST = 1e7; % distance to a far away point...
if ~isreal(vp1), vp1 = INFDIST * vp1 / j; end
if ~isreal(vp2), vp2 = INFDIST * vp2 / j; end
if ~isreal(vp3), vp3 = INFDIST * vp3 / j; end

%% order vp
imgcenter = [imgsize(2) imgsize(1)]/2;
if norm(vp2-imgcenter) < norm(vp3-imgcenter)
    temp = vp2;
    vp2 = vp3;
    vp3 = temp;
end

vp = {vp1; vp2; vp3};

%%
function [vp2 vp3 f] = find_other_vp(lines_other, vp1, THRES_THETA, imgsize)

NUMITER = 200; % number of lines to sample
FSWEEP = 50:50:1000; % focal length to test

if length(lines_other) > NUMITER
    linesample = randsample(length(lines_other), NUMITER);
else
    linesample = 1:length(lines_other);
end

scorearray = zeros(length(linesample), length(FSWEEP));
for i = 1:length(linesample)
for j = 1:length(FSWEEP)
    line = lines_other(linesample(i));
    f = FSWEEP(j);

    [vp2 vp3] = vpset_from_vp1_line_f(vp1, line, f, imgsize);

    % line assignment...
    [vv2 vv2angle] = line_belongto_vp(lines_other, vp2, THRES_THETA);
    [vv3 vv3angle] = line_belongto_vp(lines_other, vp3, THRES_THETA);
    vvboth = vv2 & vv3;
    vvboth2 = vvboth & (vv2angle<vv3angle);
    vvboth3 = vvboth & (vv2angle>=vv3angle);
    vvvv2 = (vv2 & ~vv3) | vvboth2;
    vvvv3 = (vv3 & ~vv2) | vvboth3;

    % score
    score = getvote(vp2, lines_other(vvvv2), THRES_THETA) + ...
            getvote(vp3, lines_other(vvvv3), THRES_THETA);

    scorearray(i,j) = score;
end
end

[c idx] = max(scorearray(:));
[row col] = ind2sub(size(scorearray), idx);

line = lines_other(linesample(row));
f = FSWEEP(col);
[vp2 vp3] = vpset_from_vp1_line_f(vp1, line, f, imgsize);

%%
function [vp2 vp3] = vpset_from_vp1_line_f(vp1, line, f, imgsize)

[vl p1 p2] = vp2vl(vp1, f, imgsize);

% intersection of chosen line and vl
% [intpt degen] = line_intersect(p1, p2, line.point1, line.point2);
intpt = line_intersect_forvp(p1, p2, line.point1, line.point2);

% second vp is the intpt
% third vp is orthogonal to the 1st and 2nd
vp2 = intpt;
% Is3 = cross([vp1(:)-imgsize([2 1])'/2; f], [vp2(:)-imgsize([2 1])'/2; f]);
Is3 = cross( xy2is(vp1, f, imgsize), xy2is(vp2, f, imgsize) );
vp3 = Is3(1:2)' / Is3(3) * f + imgsize([2 1])/2;

%%
function [vl p1 p2] = vp2vl(vp, f, imgsize)
% returns the equation for the vanishing line of planes
% that are orthogonal to the lines which are
% associated with the vanishing point.
%
% vl [3x1]: vl(1)*x + vl(2)*y + vl(3) = 0
% p1, p2: two points on line

% principal point (assume at image center)
cx = imgsize(2)/2;
cy = imgsize(1)/2;

% if vp is at infinity...
if ~isreal(vp)
%     vl is a line passing through image center.
    a = vp(1) / j;
    b = vp(2) / j;
    c = - a*cx - b*cy;
    vl = [a b c];
    p1 = [cx cy];
    p2 = [cx+b cy-a];
    return;
end

% normalized coord of vp
nx = vp(1) - cx;
ny = vp(2) - cy;

% ray
% r = [nx ny f]';

% ortho plane
% p = r / norm(r);

% direction of intersecting line
% [nx ny f]' x [0 0 1]' = [ny -nx 0]

% one point on line
% on imaging plane: z0 = f
% on ortho plane:   nx*nx0 + ny*ny0 + f*z0 = 0
if abs(nx) > abs(ny)
    z0 = f;
    ny0 = 0;
    nx0 = - (ny * ny0 + f * f) / nx;
else
    z0 = f;
    nx0 = 0;
    ny0 = - (nx * nx0 + f * f) / ny;
end
x0 = nx0 + cx;
y0 = ny0 + cy;

% line eq in 3D
% [nx0 ny0 z0]' + lambda * [ny -nx 0]

% line eq in 2D
% [x0 y0]' + lambda * [ny -nx]
% <==>
% ax + by + c = 0;
a = nx;
b = ny;
c = - a * x0 - b * y0;

vl = [a b c]';

p1 = [x0 y0];
p2 = [x0 y0] + [ny -nx];

%%
function is = xy2is(xy, f, imgsize)

if ~isreal(xy)
    is = [xy(1)/j; xy(2)/j; 0];
    is = is / norm(is);
else
    is = [xy(1)-imgsize(2)/2; xy(2)-imgsize(1)/2; f];
    is = is / norm(is);
end

author: zouchuhang

repo: LayoutNet


text:
function u = fem2d_bvp_linear ( nx, ny, a, c, f, x, y )

%*****************************************************************************80
%
%% FEM2D_BVP_LINEAR solves boundary value problem on a rectangle.
%
%  Discussion:
%
%    The program uses the finite element method, with piecewise linear basis
%    functions to solve a 2D boundary value problem over a rectangle
%
%    The following differential equation is imposed inside the region:
%
%      - d/dx a(x,y) du/dx - d/dy a(x,y) du/dy + c(x,y) * u(x,y) = f(x,y)
%
%    where a(x,y), c(x,y), and f(x,y) are given functions.
%
%    On the boundary, the solution is constrained to have the value 0.
%
%    The finite element method will use a regular grid of NX nodes in X, and
%    NY nodes in Y.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of X and Y grid values.
%
%    Input, function A(X,Y), evaluates a(x,y);
%
%    Input, function C(X,Y), evaluates c(x,y);
%
%    Input, function F(X,Y), evaluates f(x,y);
%
%    Input, real X(NX), Y(NY), the mesh points.
%
%    Output, real U(NX,NY), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%

%
%  Quadrature definitions.
%
  quad_num = 3;
  abscissa(1) = -0.774596669241483377035853079956;
  abscissa(2) = 0.000000000000000000000000000000;
  abscissa(3) = 0.774596669241483377035853079956;
  weight(1) = 0.555555555555555555555555555556;
  weight(2) = 0.888888888888888888888888888889;
  weight(3) = 0.555555555555555555555555555556;
%
%  Make room for the matrix A and right hand side b.
%
  mn = nx * ny;
  A = zeros ( mn, mn );
  b = zeros ( mn, 1 );
%
%  Compute the matrix entries by integrating over each element.
%
  for ex = 1 : nx - 1

    w = ex;
    e = ex + 1;

    xw = x(w);
    xe = x(e);

    for ey = 1 : ny - 1

      s = ey;
      n = ey + 1;

      ys = y(s);
      yn = y(n);

      sw = ( ey - 1 ) * nx + ex;
      se = ( ey - 1 ) * nx + ex + 1;
      nw =   ey       * nx + ex;
      ne =   ey       * nx + ex + 1;

      for qx = 1 : quad_num

        xq = ( ( 1.0 - abscissa(qx) ) * xw   ...
             + ( 1.0 + abscissa(qx) ) * xe ) ...
             / 2.0;

        for qy = 1 : quad_num

          yq = ( ( 1.0 - abscissa(qy) ) * ys   ...
               + ( 1.0 + abscissa(qy) ) * yn ) ...
               / 2.0;

          wq = weight(qx) * ( xe - xw ) / 2.0 ...
             * weight(qy) * ( yn - ys ) / 2.0;

          vsw  = ( xe - xq ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
          vswx = (    -1.0 ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
          vswy = ( xe - xq ) / ( xe - xw ) * (    -1.0 ) / ( yn - ys );

          vse  = ( xq - xw ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
          vsex = ( 1.0     ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
          vsey = ( xq - xw ) / ( xe - xw ) * (    -1.0 ) / ( yn - ys );

          vnw  = ( xe - xq ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
          vnwx = (    -1.0 ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
          vnwy = ( xe - xq ) / ( xe - xw ) * ( 1.0     ) / ( yn - ys );

          vne  = ( xq - xw ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
          vnex = ( 1.0     ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
          vney = ( xq - xw ) / ( xe - xw ) * ( 1.0     ) / ( yn - ys );

          aq = a ( xq, yq );
          cq = c ( xq, yq );
          fq = f ( xq, yq );

          A(sw,sw) = A(sw,sw) + wq * ( vswx * aq * vswx ...
                                     + vswy * aq * vswy ...
                                     + vsw  * cq * vsw );
          A(sw,se) = A(sw,se) + wq * ( vswx * aq * vsex ...
                                     + vswy * aq * vsey ...
                                     + vsw  * cq * vse );
          A(sw,nw) = A(sw,nw) + wq * ( vswx * aq * vnwx ...
                                     + vswy * aq * vnwy ...
                                     + vsw  * cq * vnw );
          A(sw,ne) = A(sw,ne) + wq * ( vswx * aq * vnex ...
                                     + vswy * aq * vney ...
                                     + vsw  * cq * vne );
          b(sw) =    b(sw)    + wq * ( vsw * fq );

          A(se,sw) = A(se,sw) + wq * ( vsex * aq * vswx ...
                                     + vsey * aq * vswy ...
                                     + vse  * cq * vsw );
          A(se,se) = A(se,se) + wq * ( vsex * aq * vsex ...
                                     + vsey * aq * vsey ...
                                     + vse  * cq * vse );
          A(se,nw) = A(se,nw) + wq * ( vsex * aq * vnwx ...
                                     + vsey * aq * vnwy ...
                                     + vse  * cq * vnw );
          A(se,ne) = A(se,ne) + wq * ( vsex * aq * vnex ...
                                     + vsey * aq * vney ...
                                     + vse  * cq * vne );
          b(se) =    b(se)    + wq * ( vse * fq );

          A(nw,sw) = A(nw,sw) + wq * ( vnwx * aq * vswx ...
                                     + vnwy * aq * vswy ...
                                     + vnw  * cq * vsw );
          A(nw,se) = A(nw,se) + wq * ( vnwx * aq * vsex ...
                                     + vnwy * aq * vsey ...
                                     + vnw  * cq * vse );
          A(nw,nw) = A(nw,nw) + wq * ( vnwx * aq * vnwx ...
                                     + vnwy * aq * vnwy ...
                                     + vnw  * cq * vnw );
          A(nw,ne) = A(nw,ne) + wq * ( vnwx * aq * vnex ...
                                     + vnwy * aq * vney ...
                                     + vnw  * cq * vne );
          b(nw) =    b(nw)    + wq * ( vnw * fq );

          A(ne,sw) = A(ne,sw) + wq * ( vnex * aq * vswx ...
                                     + vney * aq * vswy ...
                                     + vne  * cq * vsw );
          A(ne,se) = A(ne,se) + wq * ( vnex * aq * vsex ...
                                     + vney * aq * vsey ...
                                     + vne  * cq * vse );
          A(ne,nw) = A(ne,nw) + wq * ( vnex * aq * vnwx ...
                                     + vney * aq * vnwy ...
                                     + vne  * cq * vnw );
          A(ne,ne) = A(ne,ne) + wq * ( vnex * aq * vnex ...
                                     + vney * aq * vney ...
                                     + vne  * cq * vne );
          b(ne) =    b(ne)    + wq * ( vne * fq );

        end
      end
    end
  end
%
%  Where a node is on the boundary,
%  replace the finite element equation by a boundary condition.
%
  k = 0;
  for y = 1 : ny
    for x = 1 : nx
      k = k + 1;
      if ( x == 1 | x == nx | y == 1 | y == ny )
        A(k,1:mn) = 0.0;
        A(1:mn,k) = 0.0;
        A(k,k) = 1.0;
        b(k) = 0.0;
      end
    end
  end

  if ( 0 )
    spy ( A );
    pause
  end
%
%  Solve the linear system.
%
  u = A \ b;
%
%  Make the vector U into a matrix.
%  Hope that of the various possible reorderings, you and MATLAB agree.
%
%    U(matrix)   U(vector)
%   ----------  ----------
%    U11 U12  =>  U11
%    U21 U22      U21
%                 U12
%                 U22
%
  u = reshape ( u, nx, ny );

  return
end


text:
function [ a, b, c ] = sphere_triangle_sides_to_angles ( r, as, bs, cs )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_SIDES_TO_ANGLES computes spherical triangle angles in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real AS, BS, CS, the (geodesic) length of the
%    sides of the triangle.
%
%    Output, real A, B, C, the spherical angles of the triangle.
%    Angle A is opposite the side of length AS, and so on.
%
  asu = as / r;
  bsu = bs / r;
  csu = cs / r;
  ssu = ( asu + bsu + csu ) / 2.0;

  tan_a2 = sqrt ( ( sin ( ssu - bsu ) * sin ( ssu - csu ) ) / ...
                  ( sin ( ssu ) * sin ( ssu - asu )     ) );

  a = 2.0 * atan ( tan_a2 );

  tan_b2 = sqrt ( ( sin ( ssu - asu ) * sin ( ssu - csu ) ) / ...
                  ( sin ( ssu ) * sin ( ssu - bsu )     ) );

  b = 2.0 * atan ( tan_b2 );

  tan_c2 = sqrt ( ( sin ( ssu - asu ) * sin ( ssu - bsu ) ) / ...
                  ( sin ( ssu ) * sin ( ssu - csu )     ) );

  c = 2.0 * atan ( tan_c2 );

  return
end


text:
function prediction = ResNet(varargin)
%ResNet Returns a ResNet-50/101/152/custom model for ImageNet
%   M = models.ResNet() returns the ResNet-50 model proposed in:
%
%     He et al., "Deep Residual Learning for Image Recognition", CVPR 2016.
%
%   models.ResNet(..., 'option', value, ...) accepts the following
%   options:
%
%   `variant`:: '50'
%     Model variant: 50, 101, 152 or 'custom'.
%
%   `pretrained`:: false
%     If true, returns a model pre-trained on ImageNet (using the
%     MatConvNet example code).
%
%   `input`:: default input
%     Specifies an input (images) layer for the network. If unspecified, a
%     new one is created.
%
%   `numClasses`:: 1000
%     Number of output classes.
%
%   `blocksPerSection`:: depends on variant
%     Vector with the number of blocks in each section. This is set
%     automatically depending on the variant, unless it is 'custom'. The
%     standard variants have only 4 sections (thus this is a 4-elements
%     vector), but different numbers of sections are supported.
%
%   `downsampleSection`:: [0, 1, 1, 1]
%     Whether to downsample spatially at each section (binary).
%
%   `channelsPerSection`:: [64, 128, 256, 512]
%     Base number of output channels in each section.
%
%   `channelsPerBlock`:: [1, 1, 4]
%     Modifying factor for number of channels per conv layer in a block.
%
%   `initialSize`:: [7, 7, 3, 64]
%     Filter size of the input convolution (4D tensor).
%
%   Any other options will be passed to models.ConvBlock(), and can be used
%   to change the activation function, weight initialization, etc.
%
%   Suggested SGD training options are also returned in the struct M.meta.

% Copyright (C) 2018 Joao F. Henriques, Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  % parse options. unknown arguments will be passed to ConvBlock (e.g.
  % activation).
  opts.pretrained = false ;  % whether to fetch a pre-trained model
  opts.variant = '50';  % model variant: 50, 101, 152 or 'custom'
  opts.input = Input('name', 'images', 'gpu', true) ;  % default input layer
  opts.numClasses = 1000 ;  % number of predicted classes
  opts.blocksPerSection = [] ;  % vector with #blocks per section (depends on variant)
  opts.downsampleSection = [0, 1, 1, 1] ;  % whether to downsample in each section
  opts.channelsPerSection = [64, 128, 256, 512] ;  % #channels per section
  opts.channelsPerBlock = [1, 1, 4] ;  % modifying factor for #channels per conv in a block
  opts.initialSize = [7, 7, 3, 64] ;  % filter size of the input convolution (4D tensor)
  [opts, convBlockArgs] = vl_argparse(opts, varargin, 'nonrecursive') ;

  % set the number of blocks per section for standard ResNet variants
  if isnumeric(opts.variant)  % accept variant 50 instead of '50'
    opts.variant = int2str(opts.variant) ;
  end
  switch opts.variant
  case '18'
    opts.blocksPerSection = [2, 2, 2, 2] ;
    opts.channelsPerBlock(3:end) = [];  % only two convolutions per block
  case '34'
    opts.blocksPerSection = [3, 4, 6, 3] ;
    opts.channelsPerBlock(3:end) = [];  % only two convolutions per block
  case '50'
    opts.blocksPerSection = [3, 4, 6, 3] ;
  case '101'
    opts.blocksPerSection = [3, 4, 23, 3] ;
  case '152'
    opts.blocksPerSection = [3, 8, 36, 3] ;
  case 'custom'
    assert(~isempty(opts.blocksPerSection), ...
      'To build a custom ResNet, blocksPerSection must be specified.') ;
  otherwise
    error('Unknown variant.') ;
  end

  assert(any(numel(opts.channelsPerBlock) == [2, 3]), ...
    'Can only specify 2 or 3 channels per block.') ;

  % default training options for this network (returned as output.meta)
  meta.batchSize = 64 ;  % 256 with subbatches
  meta.imageSize = [224, 224, 3] ;
  meta.augmentation.crop = 224 / 256;
  meta.augmentation.location = true ;
  meta.augmentation.flip = true ;
  meta.augmentation.brightness = 0.1 ;
  meta.augmentation.aspect = [3/4, 4/3] ;
  meta.augmentation.scale = [0.4, 1.1] ;
  meta.weightDecay = 0.0001 ;
  meta.momentum = 0.9 ;

  % the default learning rate schedule
  if ~opts.pretrained
    meta.learningRate = [0.1 * ones(1,30), 0.01 * ones(1,30), 0.001 * ones(1,30)] ;
    meta.numEpochs = numel(meta.learningRate) ;
  else  % fine-tuning has lower LR
    meta.learningRate = 1e-5 ;
    meta.numEpochs = 20 ;
  end

  % return a pre-trained model
  if opts.pretrained
    if opts.numClasses ~= 1000
      warning('Model options are ignored when loading a pre-trained model.') ;
    end
    prediction = models.pretrained(['imagenet-resnet-' opts.variant '-dag']) ;

    % return prediction layer (not softmax)
    assert(isequal(prediction{1}.func, @vl_nnsoftmax)) ;
    prediction = prediction{1}.inputs{1} ;

    % replace input layer with the given one
    input = prediction.find('Input', 1) ;
    prediction.replace(input, opts.input) ;

    % in pre-trained networks, the last pooling is fixed to have size 7x7.
    % to make it compatible with any input size, replace it with mean().
    pool = prediction.inputs{1} ;
    assert(isequal(pool.func, @vl_nnpool));  % sanity check

    poolIn = pool.inputs{1} ;  % input to pool
    prediction.inputs{1} = mean(mean(poolIn, 1), 2) ;  % replace pooling

    prediction.meta = meta ;
    return
  end

  % get conv block generator with the given options (can be overriden)
  conv = models.ConvBlock('hasBias', false, 'pad', 'same', 'batchNorm', true, ...
    'preActivationBatchNorm', true, 'batchNormEpsilon', 1e-5, convBlockArgs{:}) ;

  % build network
  images = opts.input ;

  % add input section
  x = conv(images, 'size', opts.initialSize, 'stride', 2) ;
  x = vl_nnpool(x, 3, 'stride', 2, 'pad', 1)  ;
  channelsOut = opts.initialSize(4) ;

  for s = 1 : numel(opts.blocksPerSection)  % iterate sections
    for l = 1 : opts.blocksPerSection(s)  % iterate segments
      % output channels of each conv in this segment
      ch = opts.channelsPerSection(s) * opts.channelsPerBlock ;

      % downsample if first segment in section from section 2 onwards
      stride = 1 ;
      if l == 1 && opts.downsampleSection(s)
        stride = 2 ;
      end

      % this will be connected to the sum later
      sumInput = x ;

      if numel(ch) == 3
        % 3 convs for deeper variants: 1x1, 3x3, 1x1
        x = conv(x, 'size', [1, 1, channelsOut, ch(1)]) ;
        x = conv(x, 'size', [3, 3, ch(1), ch(2)], 'stride', stride) ;
        x = conv(x, 'size', [1, 1, ch(2), ch(3)], 'activation', 'none') ;
      else
        % 2 convs for shallower variants: 3x3, 3x3
        x = conv(x, 'size', [3, 3, channelsOut, ch(1)], 'stride', stride) ;
        x = conv(x, 'size', [3, 3, ch(1), ch(2)], 'activation', 'none') ;
      end

      % optional adapter layer, parallel to the above 3 convs
      if l == 1 && channelsOut ~= ch(end)
        sumInput = conv(sumInput, 'size', [1, 1, channelsOut, ch(end)], ...
          'stride', stride, 'activation', 'none') ;
      end

      % sum layer
      x = sumInput + x ;
      x = vl_nnrelu(x) ;

      channelsOut = ch(end);
    end
  end

  % average pooling
  x = mean(mean(x, 1), 2) ;

  % prediction layer
  prediction = conv(x, 'size', [1, 1, channelsOut, opts.numClasses], ...
    'batchNorm', false, 'activation', 'none', 'hasBias', true) ;

  prediction.meta = meta ;

end

repo: MFT


text:
function results_image = PTKColourBranchesBelowLobe(start_branches, airway_tree, template)
    % PTKColourBranchesBelowLobe. Given a set of labelled branches, this creates an output image with all the subtrees coloured by lobe
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    left_upper_startindices = start_branches.LeftUpper;
    left_lower_startindices = start_branches.LeftLower;
    right_upper_startindices = start_branches.RightUpper;
    right_mid_startindices = start_branches.RightMid;
    right_lower_startindices = start_branches.RightLower;
    uncertain_segments = start_branches.LeftUncertain;

    left_upper_startindices = GetOptionalFieldArray(left_upper_startindices, 'Children');
    left_lower_startindices = GetOptionalFieldArray(left_lower_startindices, 'Children');
    right_upper_startindices = GetOptionalFieldArray(right_upper_startindices, 'Children');
    right_mid_startindices = GetOptionalFieldArray(right_mid_startindices, 'Children');
    right_lower_startindices = GetOptionalFieldArray(right_lower_startindices, 'Children');
    uncertain_segments = GetOptionalFieldArray(uncertain_segments, 'Children');

    start_segments = [left_upper_startindices, left_lower_startindices, right_upper_startindices, ...
        right_mid_startindices, right_lower_startindices, uncertain_segments];

    results_image = PTKDivideAirwayTreeByCentrelineBranches(start_segments, airway_tree, template);
end

function field_array = GetOptionalFieldArray(structure, field_name)
    if isfield(structure, field_name)
        field_array = [structure.(field_name)];
    else
        field_array = [];
    end
end

function results_image = PTKDivideAirwayTreeByCentrelineBranches(start_segments, airway_tree, template)
    number_startsegments = length(start_segments);

    start_segment_voxels = [];
    extended_start_segment_voxels = [];

    results_image = zeros(template.ImageSize, 'uint8');

    % Get centreline voxels for each starting segment
    for index = 1 : number_startsegments
        start_segment_voxels{index} = PTKTreeUtilities.GetCentrelineVoxelsForTheseBranches(start_segments(index), template);
        extended_start_segment_voxels{index} = GetCentrelineVoxels(start_segments(index), template);
    end

    segments_to_do = airway_tree;

    while ~isempty(segments_to_do)
        segment = segments_to_do(end);
        segments_to_do(end) = [];
        voxel_indices_in_segment = segment.GetAllAirwayPoints;
        voxel_indices_in_segment = template.GlobalToLocalIndices(voxel_indices_in_segment);

        segment_indices = [];

        for index = 1 : number_startsegments
            if any(ismember(extended_start_segment_voxels{index}, voxel_indices_in_segment))
                segment_indices = [segment_indices, index];
            end
        end

        if length(segment_indices) == 1
            results_image(voxel_indices_in_segment) = segment_indices;
        end

        segments_to_do = [segments_to_do segment.Children];
    end
end

function voxel_indices = GetCentrelineVoxels(start_indices, template)
    voxel_indices = PTKTreeUtilities.GetCentrelineVoxelsForTheseBranches(start_indices, template);

    % Add nearest neighbours to the list of voxels, otherwise it is possible for
    % a diagnoally-connected centreline segment to pass through a
    % diagnoally-connected airway segment
    voxel_indices = MimImageCoordinateUtilities.AddNearestNeighbours(voxel_indices, template);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [errorCode, ima] = pco_sdk_example_single(exposure_time,triggermode)
% set variables grab and display a single images
%
%   [errorCode] = pco_sdk_example_single(exposure_time,triggermode)
%
% * Input parameters :
%    exposure_time           camera exposure time (default=10ms)
%    triggermode             camera trigger mode (default=AUTO)
%
% * Output parameters :
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%
%grab images from a recording pco camera
%using script function pco_camera_stack
%display the grabbed images
%
%function workflow
%open camera
%set variables
%start camera
%grab image
%show image
%stop camera
%close camera
%

glvar=struct('do_libunload',0,'do_close',0,'camera_open',0,'out_ptr',[]);

if(~exist('exposure_time','var'))
 exposure_time = 10;
end

if(~exist('triggermode','var'))
 triggermode = 0;
end

%reduce_display_size=1: display only top-left corner 800x600Pixel
reduce_display_size=1;

pco_camera_load_defines();
subfunc=pco_camera_subfunction();
subfunc.fh_lasterr(0);

[errorCode,glvar]=pco_camera_open_close(glvar);
pco_errdisp('pco_camera_setup',errorCode);
disp(['camera_open should be 1 is ',int2str(glvar.camera_open)]);
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

out_ptr=glvar.out_ptr;

try

subfunc.fh_stop_camera(out_ptr);

cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,~,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
pco_errdisp('PCO_GetCameraDescription',errorCode);

if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
 subfunc.fh_enable_timestamp(out_ptr,TIMESTAMP_MODE_BINARYANDASCII);
end

%enable MetaData if available
if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_METADATA))
 subfunc.fh_set_metadata_mode(out_ptr,1);
end

subfunc.fh_set_exposure_times(out_ptr,exposure_time,2,0,2)
subfunc.fh_set_triggermode(out_ptr,triggermode);

%if PCO_ArmCamera does fail no images can be grabbed

hwio_sig=libstruct('PCO_Signal');
set(hwio_sig,'wSize',hwio_sig.structsize);
[errorCode,~,hwio_sig] = calllib('PCO_CAM_SDK', 'PCO_GetHWIOSignal', out_ptr,0,hwio_sig);
pco_errdisp('PCO_GetHWIOSignal',errorCode);
hwio_sig.wEnabled = 1;
[errorCode,~,~] = calllib('PCO_CAM_SDK', 'PCO_SetHWIOSignal', out_ptr,0,hwio_sig);
pco_errdisp('PCO_SetHWIOSignal',errorCode);

errorCode = calllib('PCO_CAM_SDK', 'PCO_ArmCamera', out_ptr);
if(errorCode~=PCO_NOERROR)
 pco_errdisp('PCO_ArmCamera',errorCode);
 ME = MException('PCO_ERROR:ArmCamera','Cannot continue script with not armed camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

%adjust transfer parameter if necessary
subfunc.fh_set_transferparameter(out_ptr);
subfunc.fh_get_triggermode(out_ptr);
subfunc.fh_show_frametime(out_ptr);

disp('get single image');
%subfunc.fh_start_camera(out_ptr);
[errorCode,ima,metadata,glvar]=pco_camera_stack(1,glvar);
if(errorCode==PCO_NOERROR)
 if(reduce_display_size~=0)
  [ys,xs]=size(ima);
  xmax=800;
  ymax=600;
  if((xs>xmax)&&(ys>ymax))
   imar=ima(1:ymax,1:xmax);
  elseif(xs>xmax)
   imar=ima(:,1:xmax);
  elseif(ys>ymax)
   imar=ima(1:ymax,:);
  end
 else
  imar=ima;
 end
 m=max(max(imar(10:end-10,10:end-10)));
% imshow(ima',[0,m+100]);
 draw_image(imar,[0 m+100]);
 disp(['found max ',int2str(m)]);
end

subfunc.fh_stop_camera(out_ptr);

if(~isempty(metadata))
 [metadatasize,~]=size(metadata);
 [errorCode_m,metastruct]=subfunc.fh_get_struct_metadata(metadata,metadatasize);
 if(errorCode_m)
  pco_errdisp('get_struct_metadata',errorCode_m);
 else
  txt=subfunc.fh_print_meta_timestamp(metastruct);
  disp(['Timestamp of meta (',num2str(1,'%04d'),'):  ',txt]);
  subfunc.fh_print_meta_struct(metastruct);
 end
end

disp('Press "Enter" to proceed')
pause();
close()

catch ME
 errorCode=subfunc.fh_lasterr();
 txt=blanks(101);
 txt=calllib('PCO_CAM_SDK','PCO_GetErrorTextSDK',pco_uint32err(errorCode),txt,100);

 clearvars -except ME glvar;

 if(glvar.camera_open==1)
  glvar.do_close=1;
  glvar.do_libunload=1;
  pco_camera_open_close(glvar);
 end

 if strfind(ME.identifier,'PCO_ERROR:')
  msg=[ME.identifier,' ',ME.message];
  disp(txt);
  warning('off','backtrace')
  warning(msg)
  for k=1:length(ME.stack)
    disp(['from file ',ME.stack(k).file,' at line ',num2str(ME.stack(k).line)]);
  end

  clearvars;
  commandwindow;
  return;
 else
  clearvars -except ME;
  rethrow(ME)
 end
end

clearvars -except glvar ima errorCode;

if(glvar.camera_open==1)
 glvar.do_close=1;
 glvar.do_libunload=1;
 pco_camera_open_close(glvar);
end

clearvars glvar;
commandwindow;
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function varargout = vals2coeffs( U, S, V )
%VALS2COEFFS    Convert matrix of values to Chebyshev coefficients.
%
% V = VALS2COEFFS( C ) converts a matrix C of values representing samples of a
% function from a tensor Chebyshev grid and converts them to a matrix V of
% bivariate Chebyshev coefficients for the corresponding interpolant.
%
% [U, S, V] = VALS2COEFFS( C, D, R ) the same as above but keeps everything in
% low rank form. Here, C*D*R.' is a sample on the Chebyshev-Chebyshev
% tensor grid.
%
% See also COEFFS2VALS

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

tech = chebfunpref().tech();

if ( nargin == 1 )
    U = tech.vals2coeffs( U );
    U = tech.vals2coeffs( U.' ).';
    varargout = {U};
elseif ( nargin == 3 )
    U = tech.vals2coeffs( U );
    V = tech.vals2coeffs( V );
    varargout = {U S V};
else
    error('CHEBFUN:CHEBFUN2:vals2coeffs:inputs', ...
        'The number of input arguments should be one or two.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Figures for first whole cell paper
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 12/19/2011
classdef PaperFigures
    methods (Static)
        function run()
            %import classes
            import edu.stanford.covert.cell.sim.analysis.Figures34;
            import edu.stanford.covert.cell.sim.analysis.SingleGeneDeletions;

            %base directory for output
            baseDir = 'documentation/paper/figures/';

            %% Main paper
            %figure 1: designed with Adobe Illustrator
            %- documentation/paper/figures/model overview.ai
            %- documentation/paper/figures/algorithm overview with process-state network.ai

            %figure 2, 3, 5: plotted entirely with MATLAB
            Figures34.run('documentation/paper/figures');

            %figure 4: combination of MATLAB, illustrator
            Figures34.run('documentation/paper/figures');

            %figure 6: plotted entirely with MATLAB
            SingleGeneDeletions.run();                %run analysis and cache results
            SingleGeneDeletions.plotOverview(...      %plot overview figure from cached results
                [baseDir filesep 'figure6.pdf'], ...
                [baseDir filesep 'figure6.xls']);

            %figure 7: combination of MATLAB, illustrator

            %% Supplement
            %state, process schematics: designed with Adobe Illustrator
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function s = sum(A)
	%Get first non-singleton dimension
	sz = size(A);
	d = find(size(A)~=1);
	if isempty(d)
		s = A;
		return
	end
	d=d(1);

	idx = cell(1,ndims(A));
	for k = 1:ndims(A)
		idx{k} = ':';
	end

	ssz = sz;
	ssz(d) = 1;
	s = zeros(ssz);
	for k = 1:size(A,d)
		idx{d} = k;
		s = s + A(idx{:});
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function [p,pq,xx] = lgpdens(x,varargin)
%LGPDENS Logistic-Gaussian Process density estimate for 1D and 2D data
%
%  Description
%    LGPDENS(X,OPTIONS) Compute and plot LGP density estimate. X is
%    1D or 2D point data. For 1D data plot the mean and 95% region.
%    For 2D data plot the density contours.
%
%    [P,PQ,XT] = LGPDENS(X,OPTIONS) Compute LGP density estimate
%    and return mean density P, 2.5% and 97.5% percentiles PQ, and
%    grid locations.
%
%    [P,PQ,XT] = LGPDENS(X,XT,OPTIONS) Compute LGP density estimate
%    in the given grid locations XT.
%
%    OPTIONS is optional parameter-value pair
%      gridn     - optional number of grid points used in each axis direction
%                  default is 400 for 1D, 20 for 2D.
%      range     - tells the estimation range, default is
%                  [min(min(x),mean(x)-3*std(x)), max(max(x),mean(x)+3*std(x))]
%                  for 1D [XMIN XMAX]
%                  for 2D [X1MIN X1MAX X2MIN X2MAX]
%      gpcf      - optional function handle of a GPstuff covariance function
%                      (default is @gpcf_sexp)
%      latent_method - optional 'Laplace' (default) or 'MCMC'
%      int_method    - optional 'mode' (default), 'CCD' or 'grid'
%                      if latent_method is 'MCMC' then int_method is 'MCMC'
%      display   - defines if messages are displayed.
%                  'off' (default) displays no output
%                  'on' gives some output
%                  'iter' displays output at each iteration
%      speedup   - defines if speed-up is used.
%                  'off' (default) no speed-up is used
%                  'on' With SEXP or EXP covariance function in 2D case
%                  uses Kronecker product structure and approximates the
%                  full posterior with a low-rank approximation. Otherwise
%                  with SEXP, EXP, MATERN32 and MATERN52 covariance
%                  functions in 1D and 2D cases uses FFT/FFT2 matrix-vector
%                  multiplication speed-up in the Newton's algorithm.
%      cond_dens - defines if conditional density estimate is computed.
%                  'off' (default) no conditional density
%                  'on' computes for 2D the conditional median density
%                  estimate p(x2|x1) when the matrix [x1 x2] is given as
%                  input.
%      basis_function - defines if basis functions are used.
%                       'on' (default) uses linear and quadratic basis
%                       functions
%                       'off' no basis functions

% Copyright (c) 2011-2012 Jaakko Riihimäki and Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LGPDENS';
  ip.addRequired('x', @(x) isnumeric(x) && size(x,2)==1 || size(x,2)==2);
  ip.addOptional('xt',NaN, @(x) isnumeric(x) && size(x,2)==1 || size(x,2)==2);
  ip.addParamValue('gridn',[], @(x) isnumeric(x));
  ip.addParamValue('range',[], @(x) isempty(x)||isreal(x)&&(length(x)==2||length(x)==4));
  ip.addParamValue('gpcf',@gpcf_sexp,@(x) ischar(x) || isa(x,'function_handle'));
  ip.addParamValue('latent_method','Laplace', @(x) ismember(x,{'EP' 'Laplace' 'MCMC'}))
  %ip.addParamValue('latent_method','Laplace', @(x) ismember(x,{'EP' 'Laplace'}))
  ip.addParamValue('int_method','mode', @(x) ismember(x,{'mode' 'CCD', 'grid'}))
  ip.addParamValue('normalize',false, @islogical);
  ip.addParamValue('display', 'off', @(x) islogical(x) || ...
                   ismember(x,{'on' 'off' 'iter'}))
  ip.addParamValue('speedup',[], @(x) ismember(x,{'on' 'off'}));
  ip.addParamValue('cond_dens',[], @(x) ismember(x,{'on' 'off'}));
  ip.addParamValue('basis_function',[], @(x) ismember(x,{'on' 'off'}));

  ip.parse(x,varargin{:});
  x=ip.Results.x;
  xt=ip.Results.xt;
  gridn=ip.Results.gridn;
  xrange=ip.Results.range;
  gpcf=ip.Results.gpcf;
  latent_method=ip.Results.latent_method;
  int_method=ip.Results.int_method;
  normalize=ip.Results.normalize;
  display=ip.Results.display;
  speedup=ip.Results.speedup;
  cond_dens=ip.Results.cond_dens;
  basis_function=ip.Results.basis_function;

  [n,m]=size(x);

  switch m
    case 1 % 1D
      if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
        error('LGPDENS: the input x must be 2D if cond_dens option is ''on''.')
      end
      % Parameters for a grid
      if isempty(gridn)
        % number of points
        gridn=400;
      end
      xmin=min(x);xmax=max(x);
      if ~isempty(xrange)
        % extend given range to include min(x) and max(x)
        xmin=min(xmin,xrange(1));
        xmax=max(xmax,xrange(2));
      elseif ~isnan(xt)
        % use xt to define range and
        % extend it to include min(x) and max(x)
        xmin=min(xmin,min(xt));
        xmax=max(xmax,max(xt));
      else
        xmin=min(xmin,mean(x)-3*std(x));
        xmax=max(xmax,mean(x)+3*std(x));
      end

      % Discretize the data
      if isnan(xt)
        xx=linspace(xmin,xmax,gridn)';
      else
        xx=xt;
        gridn=numel(xt);
      end

      xd=xx(2)-xx(1);
      yy=hist(x,xx)';

      % normalise, so that same prior is ok for different scales
      xxn=(xx-mean(xx))./std(xx);

      %[Ef,Covf]=gpsmooth(xxn,yy,[xxn; xtn],gpcf,latent_method,int_method);
      [Ef,Covf]=gpsmooth(xxn,yy,xxn,gpcf,latent_method,int_method,display,speedup,gridn,cond_dens,basis_function);

      if strcmpi(latent_method,'MCMC')
        PJR=zeros(size(Ef,1),size(Covf,3));
        for i1=1:size(Covf,3)
          qr=bsxfun(@plus,randn(1000,size(Ef,1))*chol(Covf(:,:,i1),'upper'),Ef(:,i1)');
          qjr=exp(qr)';
          pjr=bsxfun(@rdivide,qjr,sum(qjr));
          pjr=pjr./xd;
          PJR(:,i1)=mean(pjr,2);
        end
        pjr=PJR;
      else
        qr=bsxfun(@plus,randn(1000,size(Ef,1))*chol(Covf,'upper'),Ef');
        qjr=exp(qr)';
        pjr=bsxfun(@rdivide,qjr,sum(qjr(1:gridn,:)));
        pjr=pjr./xd;
      end
      pp=mean(pjr')';
      ppq=prctile(pjr',[2.5 97.5])';

      if nargout<1
        % no output, do the plot thing
        newplot
        hp=patch([xx; xx(end:-1:1)],[ppq(:,1); ppq(end:-1:1,2)],[.8 .8 .8]);
        set(hp,'edgecolor',[.8 .8 .8])
        xlim([xmin xmax])
        line(xx,pp,'linewidth',2);
      else
        p=pp;
        pq=ppq;
      end

    case 2 % 2D

      if ~isempty(cond_dens) && strcmpi(cond_dens,'on') && ~isempty(speedup) && strcmp(speedup, 'on')
        warning('No speed-up option available with the cond_dens option. Using full covariance instead.')
        speedup='off';
      end
      % Find unique points
      [xu,I,J]=unique(x,'rows');
      % and count number of repeated x's
      counts=crosstab(J);
      nu=length(xu);

      % Parameters for a grid
      if isempty(gridn)
        % number of points in each direction
        gridn=20;
      end
      if numel(gridn)==1
        gridn(2)=gridn(1);
      end
      x1min=min(x(:,1));x1max=max(x(:,1));
      x2min=min(x(:,2));x2max=max(x(:,2));
      if ~isempty(xrange)
        % extend given range to include min(x) and max(x)
        x1min=min(x1min,xrange(1));
        x1max=max(x1max,xrange(2));
        x2min=min(x2min,xrange(3));
        x2max=max(x2max,xrange(4));
      elseif ~isnan(xt)
        % use xt to define range and
        % extend it to include min(x) and max(x)
        x1min=min(x1min,min(xt(:,1)));
        x1max=max(x1max,max(xt(:,1)));
        x2min=min(x2min,min(xt(:,2)));
        x2max=max(x2max,max(xt(:,2)));
      else
        x1min=min(x1min,mean(x(:,1))-3*std(x(:,1)));
        x1max=max(x1max,mean(x(:,1))+3*std(x(:,1)));
        x2min=min(x2min,mean(x(:,2))-3*std(x(:,2)));
        x2max=max(x2max,mean(x(:,2))+3*std(x(:,2)));
      end

      % Discretize the data
      if isnan(xt)
        % Form regular grid to discretize the data
        zz1=linspace(x1min,x1max,gridn(1))';
        zz2=linspace(x2min,x2max,gridn(2))';
        [z1,z2]=meshgrid(zz1,zz2);
        z=[z1(:),z2(:)];
        nz=length(z);
        xx=z;
        if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
          % use ntx2 times more grid points for predictions
          if gridn(2)>10
            ntx2=3;
          else
            ntx2=10;
          end
          zzt1=linspace(x1min,x1max,gridn(1))';
          zzt2=linspace(x2min,x2max,gridn(2)*ntx2)';
          [zt1,zt2]=meshgrid(zzt1,zzt2);
          zt=[zt1(:),zt2(:)];
          %nzt=length(zt);
          xt=zt;
        end
      else
        xx=xt;
        gridn=[length(unique(xx(:,1))) length(unique(xx(:,2)))];
      end
      yy=zeros(nz,1);
      zi=interp2(z1,z2,reshape(1:nz,gridn(2),gridn(1)),xu(:,1),xu(:,2),'nearest');
      for i1=1:nu
        yy(zi(i1),1)=yy(zi(i1),1)+counts(i1);
      end
      %ye=ones(nz,1)./nz.*n;

      unx1=unique(xx(:,1));
      unx2=unique(xx(:,2));
      xd=(unx1(2)-unx1(1))*(unx2(2)-unx2(1));

      % normalise, so that same prior is ok for different scales
      xxn=bsxfun(@rdivide,bsxfun(@minus,xx,mean(xx,1)),std(xx,1));
      if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
        xxtn=bsxfun(@rdivide,bsxfun(@minus,xt,mean(xx,1)),std(xx,1));
      end

      % [Ef,Covf]=gpsmooth(xxn,yy,[xxn; xtn],gpcf,latent_method,int_method);
      if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
        [Ef,Covf]=gpsmooth(xxn,yy,xxtn,gpcf,latent_method,int_method,display,speedup,gridn,cond_dens,basis_function);
      else
        [Ef,Covf]=gpsmooth(xxn,yy,xxn,gpcf,latent_method,int_method,display,speedup,gridn,cond_dens,basis_function);
      end

      if strcmpi(latent_method,'MCMC')

        if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
          unx2=(unique(xt(:,2)));
          xd2=(unx2(2)-unx2(1));
          PJR=zeros(size(Ef,1),size(Covf,3));
          for i1=1:size(Covf,3)
            qr=bsxfun(@plus,randn(1000,size(Ef,1))*chol(Covf(:,:,i1),'upper'),Ef(:,i1)');
            qjr=exp(qr)';
            %pjr=bsxfun(@rdivide,qjr,sum(qjr));
            pjr=qjr;
            pjr2=reshape(pjr,[gridn(2)*ntx2 gridn(1) size(pjr,2)]);
            for j1=1:size(pjr2,3)
              pjr2(:,:,j1)=bsxfun(@rdivide,pjr2(:,:,j1),sum(pjr2(:,:,j1)))./xd2;
            end
            pjr=reshape(pjr2,[gridn(2)*ntx2*gridn(1) size(pjr,2)]);
            PJR(:,i1)=mean(pjr,2);
          end
          pjr=PJR;
          %qp=median(pjr2,3);
          %qp=bsxfun(@rdivide,qp,sum(qp,1));

        else

          PJR=zeros(size(Ef,1),size(Covf,3));
          for i1=1:size(Covf,3)
            qr=bsxfun(@plus,randn(1000,size(Ef,1))*chol(Covf(:,:,i1),'upper'),Ef(:,i1)');
            qjr=exp(qr)';
            pjr=bsxfun(@rdivide,qjr,sum(qjr));
            pjr=pjr./xd;
            PJR(:,i1)=mean(pjr,2);
          end
          pjr=PJR;
          %pjr=mean(PJR,2);
        end

      else
        if strcmpi(speedup,'on') && length(Covf)==2
          qr1=bsxfun(@plus,bsxfun(@times,randn(1000,size(Ef,1)),sqrt(Covf{1})'),Ef');
          qr2=randn(1000,size(Covf{2},1))*Covf{2};
          qr=qr1+qr2;
        else
          qr=bsxfun(@plus,randn(1000,size(Ef,1))*chol(Covf,'upper'),Ef');
        end
        qjr=exp(qr)';
        if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
          pjr=zeros(size(qjr));
          unx2=unique(xt(:,2));
          xd2=(unx2(2)-unx2(1));
          for k1=1:size(qjr,2)
            qjrtmp=reshape(qjr(:,k1),[gridn(2)*ntx2 gridn(1)]);
            qjrtmp=bsxfun(@rdivide,qjrtmp,sum(qjrtmp));
            qjrtmp=qjrtmp./xd2;
            pjr(:,k1)=qjrtmp(:);
          end
        else
          pjr=bsxfun(@rdivide,qjr,sum(qjr));
          pjr=pjr./xd;
        end
      end

      %if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
      %  pp=median(pjr')';
      %else
      pp=mean(pjr')';
      %end
      ppq=prctile(pjr',[2.5 97.5])';

      if nargout<1
        % no output, do the plot thing
        if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
          pjr2=reshape(pjr,[gridn(2)*ntx2 gridn(1) size(pjr,2)]);
          %qp=median(pjr2,3);
          qp=mean(pjr2,3);
          qp=bsxfun(@rdivide,qp,sum(qp,1));
          qpc=cumsum(qp,1);
          PL=[.05 .1 .2 .5 .8 .9 .95];
          for i1=1:gridn(1)
            pc=qpc(:,i1);
            for pli=1:numel(PL),
              qi(pli)=find(pc>PL(pli),1);
            end,
            ql(:,i1)=unx2(qi);
          end
          hold on
          h1=plot(zz1,ql(4,:)','-', 'color', [0 0 255]./255,'linewidth',2);
          h2=plot(zz1,ql([3 5],:)','--', 'color', [0 127 0]./255,'linewidth',1);
          h3=plot(zz1,ql([2 6],:)','-.', 'color', [255 0 0]./255,'linewidth',1);
          h4=plot(zz1,ql([1 7],:)',':', 'color', [0 0 0]./255,'linewidth',1);
          hold off
          legend([h1 h2(1) h3(1) h4(1)],'.5','.2/.8','.1/.9','.05/.95')
          %plot(zz1,ql','linewidth',1)
          %legend('.05','.1','.2','.5','.8','.9','.95')
          xlim([x1min x1max])
          ylim([x2min x2max])
        else
          G=zeros(size(z1));
          G(:)=prctile(pjr',50);
          %contour(z1,z2,G);
          pp=G(:);
          p1=pp./sum(pp);
          pu=sort(p1,'ascend');
          pc=cumsum(pu);
          PL=[.05 .1 .2 .5 .8 .9 .95];
          qi=[];
          for pli=1:numel(PL)
            qi(pli)=find(pc>PL(pli),1);
          end
          pl=pu(qi).*sum(pp);
          contour(z1,z2,G,pl);
          %hold on, plot(x(:,1),x(:,2),'kx')
          %colorbar
        end
      else
        p=pp;
        pq=ppq;
      end
    otherwise
      error('X has to be Nx1 or Nx2')
  end
end

function [Ef,Covf] = gpsmooth(xx,yy,xxt,gpcf,latent_method,int_method,display,speedup,gridn,cond_dens,basis_function)
% Make inference with log Gaussian process and EP or Laplace approximation

  % gp_mc and gp_ia still uses numeric display option
  if strcmp(display,'off')
    displ=0;
  else
    displ=1;
  end

  nin = size(xx,2);
  % init gp
  if strfind(func2str(gpcf),'ppcs')
    % ppcs still have nin parameter...
    gpcf1 = gpcf('nin',nin);
  else
    gpcf1 = gpcf();
  end

  % weakly informative prior
  pm = prior_logunif();
  pl = prior_t('s2', 10^2, 'nu', 4);
  pa = prior_t('s2', 10^2, 'nu', 4);
  %pm = prior_sqrtt('s2', 10^2, 'nu', 4);
  %pl = prior_t('s2', 1^2, 'nu', 4);
  %pa = prior_t('s2', 10^2, 'nu', 4);
  % different covariance functions have different parameters
  if isfield(gpcf1,'magnSigma2')
     gpcf1 = gpcf(gpcf1, 'magnSigma2', .5, 'magnSigma2_prior', pm);
  end
  if isfield(gpcf1,'lengthScale')
     gpcf1 = gpcf(gpcf1, 'lengthScale', .5, 'lengthScale_prior', pl);
  end
  if isfield(gpcf1,'alpha')
    gpcf1 = gpcf(gpcf1, 'alpha', 20, 'alpha_prior', pa);
  end
  if isfield(gpcf1,'biasSigma2')
    gpcf1 = gpcf(gpcf1, 'biasSigma2', 10, 'weightSigma2', 10,'biasSigma2_prior',prior_logunif(),'weightSigma2_prior',prior_logunif());
  end

  if ~isempty(cond_dens) && strcmp(cond_dens, 'on')
    lik=lik_lgpc;
    lik.gridn=gridn;
  else
    lik=lik_lgp;
  end

  % Create the GP structure
  if ~isempty(basis_function) && strcmp(basis_function, 'off')
    gp = gp_set('lik', lik, 'cf', {gpcf1}, 'jitterSigma2', 1e-4);
  else
    %gpmfco = gpmf_constant('prior_mean',0,'prior_cov',100);
    gpmflin = gpmf_linear('prior_mean',0,'prior_cov',100);
    gpmfsq = gpmf_squared('prior_mean',0,'prior_cov',100);
    gp = gp_set('lik', lik, 'cf', {gpcf1}, 'jitterSigma2', 1e-4, 'meanf', {gpmflin,gpmfsq});
  end
  % First optimise hyperparameters using Laplace approximation
  gp = gp_set(gp, 'latent_method', 'Laplace');
  opt=optimset('TolFun',1e-2,'TolX',1e-3,'Display',display);

  if ~isempty(speedup) && strcmp(speedup, 'on')
    gp.latent_opt.gridn=gridn;
    gp.latent_opt.pcg_tol=1e-12;
    if size(xx,2)==2 && (strcmp(gp.cf{1}.type,'gpcf_sexp') || strcmp(gp.cf{1}.type,'gpcf_exp'))
      % exclude eigenvalues smaller than 1e-6 or take 50%
      % eigenvalues at most
      gp.latent_opt.eig_tol=1e-6;
      gp.latent_opt.eig_prct=0.5;
      gp.latent_opt.kron=1;
      opt.LargeScale='off';
      if norm(xx-xxt)~=0
        warning('In the low-rank approximation the grid locations xx are used instead of xxt in predictions.')
        xxt=xx;
      end
    elseif strcmp(gp.cf{1}.type,'gpcf_sexp') || strcmp(gp.cf{1}.type,'gpcf_exp') || strcmp(gp.cf{1}.type,'gpcf_matern32') || strcmp(gp.cf{1}.type,'gpcf_matern52')
      gp.latent_opt.fft=1;
    end
  end

  if exist('fminunc')
    gp=gp_optim(gp,xx,yy,'opt',opt, 'optimf', @fminunc);
  else
    gp=gp_optim(gp,xx,yy,'opt',opt, 'optimf', @fminlbfgs);
  end
  %gradcheck(gp_pak(gp), @gpla_nd_e, @gpla_nd_g, gp, xx, yy);

  if strcmpi(latent_method,'MCMC')
    gp = gp_set(gp, 'latent_method', 'MCMC');

    %if ~isempty(cond_dens) && strcmpi(cond_dens,'on')
    if size(xx,2)==2
      % add more jitter for 2D cases with MCMC
      gp = gp_set(gp, 'jitterSigma2', 1e-2);
      %error('LGPDENS: MCMC is not implemented if cond_dens option is ''on''.')
    end

    % Here we use two stage sampling to get faster convergence
    hmc_opt=hmc2_opt;
    hmc_opt.steps=10;
    hmc_opt.stepadj=0.05;
    hmc_opt.nsamples=1;
    latent_opt.display=0;
    latent_opt.repeat = 20;
    latent_opt.sample_latent_scale = 0.5;
    hmc2('state', sum(100*clock))

    % The first stage sampling
    [r,g,opt]=gp_mc(gp, xx, yy, 'hmc_opt', hmc_opt, 'latent_opt', latent_opt, 'nsamples', 1, 'repeat', 15, 'display', displ);
    %[r,g,opt]=gp_mc(gp, xx, yy, 'latent_opt', latent_opt, 'nsamples', 1, 'repeat', 15);

    % re-set some of the sampling options
    hmc_opt.steps=4;
    hmc_opt.stepadj=0.05;
    %latent_opt.repeat = 5;
    hmc2('state', sum(100*clock));

    % The second stage sampling
    % Notice that previous record r is given as an argument
    [rgp,g,opt]=gp_mc(gp, xx, yy, 'hmc_opt', hmc_opt, 'nsamples', 500,'latent_opt', latent_opt, 'record', r,  'display', displ);
    % Remove burn-in
    rgp=thin(rgp,102,4);

    [Ef, Covf] = gpmc_jpreds(rgp, xx, yy, xxt);

  else
    if strcmpi(int_method,'mode')
      % Just make prediction for the test points
      [Ef,Covf] = gp_pred(gp, xx, yy, xxt);
    else
      % integrate over the hyperparameters
      %[~, ~, ~, Ef, Covf] = gp_ia(opt, gp, xx, yy, xt, param);
      gpia=gp_ia(gp, xx, yy, 'int_method', int_method, 'display', displ);
      [Ef, Covf]=gpia_jpred(gpia, xx, yy, xxt);
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip