text:
function info = gpuinfo()
%GPUINFO  get basic details about the GPU being used
%
%   GPUINFO() gets hold of some basic info about the GPU device being used.
%
%   See also: CPUINFO

%   Copyright 2011-2012 The MathWorks, Inc.

if parallel.gpu.GPUDevice.isAvailable()
    gpu = gpuDevice();
info = struct( ...
    'Name', gpu.Name, ...
    'Clock', sprintf( '%u MHz', gpu.ClockRateKHz/1e3 ), ...
    'NumProcessors', gpu.MultiprocessorCount, ...
    'ComputeCapability', gpu.ComputeCapability, ...
    'TotalMemory', sprintf( '%1.2f GB', gpu.TotalMemory/2^30 ), ...
    'CUDAVersion', gpu.DriverVersion, ...
    'DriverVersion', parallel.internal.gpu.CUDADriverVersion );
else
    % No GPU, so create an empty structure
info = struct( ...
    'Name', '<no GPU available>', ...
    'Clock', '', ...
    'NumProcessors', 0, ...
    'ComputeCapability', '', ...
    'TotalMemory', 0, ...
    'CUDAVersion', '', ...
    'DriverVersion', '' );
end

repo: Aero-matlab


text:
function [surface_handles, colors] = isosurface(atlas_obj, varargin)
% Create a series of surfaces in different colors, one for each region
% - Options for single color
%
% surface_handles = isosurface(atlas_obj, [optional arguments])
%
% optional arguments:
% - Any optional inputs to imageCluster, e.g., 'alpha'
% - 'colors', followed by single color in { } or cell array of multiple colors
% - 'nomatchleftright', do not match colors across hemispheres (left/right)
%               Note: The default matches, and may override your colors.
%
% Examples:
% atlasfile = which('Morel_thalamus_atlas_object.mat');
% load(atlasfile)
%
% surface_handles = isosurface(atlas_obj);
% surface_handles = isosurface(atlas_obj, 'alpha', .5);
% surface_handles = isosurface(atlas_obj, 'alpha', .5, 'nomatchleftright');
%
% view(135, 30);
% lightRestoreSingle;
% lightFollowView;
%
% load(which('CIT168_MNI_subcortical_atlas_object.mat'));
%
% surface_handles = isosurface(atlas_obj, 'alpha', .5, 'color', {[.3 .6 .4] [.5 .4 .2]});
% surface_handles = isosurface(atlas_obj, 'alpha', .5, 'color', {[.3 .6 .4] [.5 .4 .2]}, 'nomatchleftright');
% p = addbrain('hires right');
% lightFollowView;

r = atlas2region(atlas_obj);

[surface_handles, colors] = isosurface(r, varargin{:});

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function gX = expKernDiagGradX(kern, X)

% EXPKERNDIAGGRADX Gradient of EXP kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the exponentiated kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : expKernParamInit, kernDiagGradX, expkernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

gX = kernDiagGradX(kern.argument, X);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function plot_McBwtime(catalog, sPar)
    % plot Mc and b as a function of time using the bootstrap approach
    % Uses the result matrix from calc_McBwtime
    % updated: 14.02.2005
    % J. Woessner
    % turned into function by Celso G Reyes 2017

    ZG = ZmapGlobal.Data; % used by get_zmap_globals

    % Get input

    % Initial values
    nSampleSize = 500;
    nOverlap    = 4;
    methodChoice = 1;
    nBstSample  = 200;
    minEventCount = 50;
    fBinning    = 0.1;
    nWindowSize = 5;
    fMcCorr     = 0;
    default_method       = McMethods.MaxCurvature; % default method
    figure_w_normalized_uicontrolunits(...
        'Name'          , 'Mc Input Parameter',...
        'NumberTitle'   , 'off',...
        'NextPlot'      , 'new',...
        'units'         , 'points',...
        'Visible','off' , ...
        'Position', [ 200 200 500 200]);
    axis off

    % Input parameters

    zdlg = ZmapDialog();
    zdlg.AddMcMethodDropdown('methodChoice', default_method);
    zdlg.AddEdit('nSampleSize'   , 'Sample window size' , nSampleSize   , 'number of events in each window');
    zdlg.AddEdit('minEventCount' , 'Min. # of events'   , minEventCount , 'minimum number of events in order to calc. window');
    zdlg.AddEdit('nOverlap'      , 'Window overlap (%)'     , nOverlap      , 'Samplesize/nOverlap determines overlap of moving windows');
    zdlg.AddEdit('nBstSample'    , 'Bootstraps'         , nBstSample    , 'Number of samples in bootstrap(?)');
    zdlg.AddEdit('fMcCorr'       , 'Mc correction'      , fMcCorr       , 'correction for the magnitude of completeness');
    zdlg.AddEdit('fBinning'      , 'Magnitude Binning'  , fBinning      , 'size of magnitude bins');
    zdlg.AddEdit('windowSize'    , 'Smooth plot'        , nWindowSize   , 'smooth plot');
    [zans,okPressed] = zdlg.Create('Name', 'Mc Input Parameter');

    if ~okPressed
        return
    end

    %% calculate  time series
    windowSize = zans.windowSize;
    catalog.sort('Date');

    [mResult] = calc_McBwtime(catalog, zans.nSampleSize, zans.nOverlap, zans.methodChoice, zans.nBstSample, zans.minEventCount, zans.fBinning, zans.fMcCorr);

    errbarParams = {'LineStyle', '-.', 'Linewidth', 2, 'Color', [0.5 0.5 0.5]};

    % Plot Mc time series
    if sPar == "mc"
        fig = figure('tag','Mc time series', 'visible','on', 'Name','Mc time series');
        resultCol   = [mResult.mcMeanWithTime];
        stdevcol    = [mResult.mcStdDevBoot];
        DisplayName = 'min Mag. Completeness (Mc)';
        errNames    = {'Mc - std', 'Mc + std'};
        myYLabel    = 'Mc';
        myLegend    = {'Mc','\delta Mc'};
    else
        % plot B-value
        fig     = figure('tag','b-value time series', 'visible','on', 'Name','b-value time series');
        resultCol   = [mResult.bMeanWithTime];
        stdevcol    = [mResult.bStdDevBoot];
        DisplayName = 'b-value';
        errNames    = {'b - std', 'b + std'};
        myYLabel    = 'b-value';
        myLegend    = {'b-value','\delta b'};
    end

    y    = filter(ones(1,windowSize)/windowSize, 1, resultCol(:)); %was mMc or mB
    yStd = filter(ones(1,windowSize)/windowSize, 1, stdevcol(:)); % was mMcstd1 or mBstd1

    if length(resultCol) > windowSize
        y(1:windowSize, 1)    = resultCol(1:windowSize);
        yStd(1:windowSize, 1) = stdevcol(1:windowSize);
    end

    x = [mResult.meanSampleTime];

    %% plot
    ax = axes(fig);
    plot(ax, x, y, '-', 'Linewidth', 2, 'Color', [0.2 0.2 0.2], 'DisplayName', DisplayName);
    ax.NextPlot = 'add';
    plot(ax, x, y-yStd, errbarParams{:}, 'DisplayName', errNames{1});
    plot(ax, x, y+yStd, errbarParams{:}, 'DisplayName', errNames{2});
    ax.NextPlot = 'replace';

    xlim(ax, bounds2(x))
    %y_bounds = [floor(min(y-yStd)), ceil(max(y+yStd))];
    %if ~any(ismissing(y_bounds))
    %    ylim(ax, [floor(min(y-yStd)), ceil(max(y+yStd))]);
    %end
    l1 = legend(ax,myLegend{:});
    set(l1, 'Fontweight', 'bold')
    set(ax, 'Fontweight', 'bold', 'FontSize', 10, 'Linewidth', 2, 'Tickdir', 'out')
    xlabel(ax, 'Time / [dec. year]', 'Fontweight', 'bold', 'FontSize',12)
    ylabel(ax, myYLabel, 'Fontweight', 'bold', 'FontSize', 12);

end

author: CelsoReyes


text:
function title = p08_title ( )

%*****************************************************************************80
%
%% P08_TITLE returns a title for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Penalty Function #1.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
classdef Radio
    methods(Static)
        function [jp, bg]=Panel(varargin)
            [~, jsa, ~, defaultIdx]=getMsgTypeAndOptions(...
                javax.swing.JOptionPane.QUESTION_MESSAGE, varargin);
            [jp,bg]=Radio.Panel2(jsa, defaultIdx, 11);
        end

        function [jp, bg]=PanelAndCallback(callback, noDefault, varargin)
            [~, jsa, ~, defaultIdx]=getMsgTypeAndOptions(...
                javax.swing.JOptionPane.QUESTION_MESSAGE, varargin);
            if noDefault
                defaultIdx=0;
            end
            [jp,bg]=Radio.Panel2(jsa, defaultIdx, 11, false, callback);
        end

        function [outerPnl, bg, checkBoxes, innerPnl]=Panel2(javaArray, ...
                dfltIdxs, addScrollAfter, checkBox, callback, ...
                priorCmps, sortIdxs, visibleIdxs)
            if nargin<4 || ~checkBox
                checkBox=false;
                jObj='javax.swing.JRadioButton';
            else
                jObj='javax.swing.JCheckBox';
            end
            N=javaArray.length;
            isJPanels=N>0 && isa(javaArray(1),'javax.swing.JPanel');
            innerPnl=javaObjectEDT('javax.swing.JPanel');
            innerPnl.setLayout(javax.swing.BoxLayout(innerPnl, 1));
            if checkBox
                bg=javaObjectEDT('javax.swing.JPanel');
                if ~isJPanels
                    columns=1;
                    horizontal=true;
                    for i=1:N
                        if javaArray(i).length>2
                            horizontal=false;
                            break;
                        else
                        end
                    end
                else
                    columns=1;
                    if N>1
                        columns=javaArray(1).getComponentCount;
                    end
                    horizontal=false;
                end
                if ~horizontal
                    if ~isJPanels
                        innerPnl.setLayout(javaObjectEDT(...
                            'java.awt.GridLayout', N, 1));
                    else
                        innerPnl.setLayout(javaObjectEDT(...
                            'java.awt.GridBagLayout'));
                        gbc=javaObjectEDT('java.awt.GridBagConstraints');
                        gbc.fill=2;
                    end
                else
                    innerPnl.setLayout(javaObjectEDT(...
                        'java.awt.GridLayout', 1, N));
                end
            else
                bg=javax.swing.ButtonGroup;
            end
            checkBoxes=java.util.ArrayList;
            idx=1;
            scrollToHeight=0;
            first=true;
            if ~isJPanels
                %javaArray is char or strings
                for i=1:N
                    rb=javaObjectEDT(jObj, javaArray(i));
                    rObj=rb;
                    pb=handle(rb, 'CallbackProperties');
                    if nargin>=5
                        set(pb,'ActionPerformedCallback', callback);
                    end
                    checkBoxes.add(rb);
                    bg.add(rb);
                    innerPnl.add(rObj);
                    if ~isempty(find(dfltIdxs==idx, 1))
                        rb.setSelected(true);
                        if first
                            d=rb.getPreferredSize;
                            scrollToHeight=d.height*(i-1);
                            first=false;
                        end
                    end
                    idx=idx+1;
                end
            else
                try
                    if nargin>5
                        priorCnt=length(priorCmps);
                        columns=priorCnt/N;
                    end
                    if nargin>6
                        if nargin<8 || isempty(visibleIdxs)
                            for i=1:N
                                idx2=(2*sortIdxs(i))-1;
                                rb=priorCmps{idx2};
                                checkBoxes.add(rb);
                                gbc.gridy=i-1;
                                gbc.gridx=0;
                                gbc.weightx=.8;
                                innerPnl.add(rb, gbc);
                                for j=2:columns
                                    gbc.gridx=j-1;
                                    gbc.weightx=.2;
                                    jp=priorCmps{(idx2+j)-1};
                                    innerPnl.add(jp, gbc);
                                end
                                idx=idx+1;
                            end
                        else
                            for i=1:N
                                si=sortIdxs(i);
                                if isempty(find(visibleIdxs==si, 1))
                                    continue;
                                end
                                idx2=(2*si)-1;
                                rb=priorCmps{idx2};
                                checkBoxes.add(rb);
                                gbc.gridy=i-1;
                                gbc.gridx=0;
                                gbc.weightx=.8;
                                innerPnl.add(rb, gbc);
                                for j=2:columns
                                    gbc.gridx=j-1;
                                    gbc.weightx=.2;
                                    jp=priorCmps{(idx2+j)-1};
                                    innerPnl.add(jp, gbc);
                                end
                                idx=idx+1;
                            end
                        end
                    else
                    for i=1:N
                            rb=javaObjectEDT(jObj, javaArray(i).getComponent(0).getText);
                            pb=handle(rb, 'CallbackProperties');
                            if nargin>=5
                                set(pb,'ActionPerformedCallback', callback);
                            end
                            checkBoxes.add(rb);
                            gbc.gridy=i-1;
                            gbc.gridx=0;
                            gbc.weightx=.8;
                            innerPnl.add(rb, gbc);
                            for j=2:columns
                                gbc.gridx=j-1;
                                gbc.weightx=.2;
                                jp=javaObjectEDT('javax.swing.JPanel');
                                jp.add(javaArray(i).getComponent(j-1));
                                innerPnl.add(jp, gbc);
                            end
                            if ~isempty(find(dfltIdxs==idx, 1))
                                rb.setSelected(true);
                                if first
                                    d=rb.getPreferredSize;
                                    scrollToHeight=d.height*(i-1);
                                    first=false;
                                end
                            end
                            idx=idx+1;
                        end
                    end
                catch ex
                    disp(['If javaArray contains JPanels then '...
                        'Radio.Panel2' newline 'expects component 0'...
                        ' to be JLabel and any 2nd component!' newline]);
                    ex.getReport
                end
            end
            outerPnl=javaObjectEDT('javax.swing.JPanel');
            outerPnl.setBorder(...
                javax.swing.BorderFactory.createEmptyBorder(1,8,1,4));
            if addScrollAfter<N
                d=innerPnl.getPreferredSize;
                d.width=floor(d.width*1.1);
                d.height=floor( (addScrollAfter/N)*d.height);
                innerPnl.setBorder([]);
                scroll=javaObjectEDT('javax.swing.JScrollPane', innerPnl);
                if scrollToHeight>0
                    scroll.getViewport.setViewPosition(...
                        java.awt.Point(0, scrollToHeight));
                end
                scroll.setPreferredSize(d);
                outerPnl=scroll;
            else
                cb=javax.swing.BorderFactory.createCompoundBorder(...
                    javax.swing.BorderFactory.createBevelBorder(...
                    javax.swing.border.BevelBorder.RAISED), ...
                    javax.swing.BorderFactory.createEmptyBorder(8, 8, 8, 8));
                innerPnl.setBorder(cb);
                outerPnl.add(innerPnl);
            end
        end

        function [outerPnl, bg, btns, sortGui]=Panel3(strs, dfltIdxs, ...
                addScrollAfter, callback, checkBox, advice)
            isRadio=nargin<5||~checkBox;
            bg=[];
            sortGui=[];
            if isRadio
                jObj='javax.swing.JRadioButton';
            else
                jObj='javax.swing.JCheckBox';
            end
            N=length(strs);
            innerPnl=javaObjectEDT('javax.swing.JPanel');
            innerPnl.setLayout(javax.swing.BoxLayout(innerPnl, 1));
            if isRadio
                bg=javax.swing.ButtonGroup;
            end
            btns=java.util.ArrayList;
            scrollToHeight=0;
            for i=1:N
                rb=javaObjectEDT(jObj, strs{i});
                rObj=rb;
                pb=handle(rb, 'CallbackProperties');
                set(pb,'ActionPerformedCallback', @(h,e)innerCallback(h,e,i));
                btns.add(rb);
                if isRadio
                    bg.add(rb);
                end
                innerPnl.add(rObj);
            end
            first=true;
            if ~isempty(dfltIdxs)
                for i=1:N
                    if ~isempty(find(dfltIdxs==i, 1))
                        rb=btns.get(i-1);
                        rb.setSelected(true);
                        if first
                            d=rb.getPreferredSize;
                            scrollToHeight=d.height*(i-1);
                            first=false;
                        end
                    end
                end
            end
            outerPnl=javaObjectEDT('javax.swing.JPanel');
            outerPnl.setBorder(...
                javax.swing.BorderFactory.createEmptyBorder(1,8,1,4));
            if addScrollAfter<N
                d=innerPnl.getPreferredSize;
                d.width=floor(d.width*1.1);
                d.height=floor( (addScrollAfter/N)*d.height);
                innerPnl.setBorder([]);
                scroll=javaObjectEDT('javax.swing.JScrollPane', innerPnl);
                if scrollToHeight>0
                    scroll.getViewport.setViewPosition(...
                        java.awt.Point(0, scrollToHeight));
                end
                scroll.setPreferredSize(d);
                outerPnl=scroll;
            else
                cb=javax.swing.BorderFactory.createCompoundBorder(...
                    javax.swing.BorderFactory.createBevelBorder(...
                    javax.swing.border.BevelBorder.RAISED), ...
                    javax.swing.BorderFactory.createEmptyBorder(8, 8, 8, 8));
                innerPnl.setBorder(cb);
                outerPnl.add(innerPnl);
            end
            if ~isRadio
                bp=Gui.BorderPanel;
                north=Gui.BorderPanel;
                bp.add(north, 'North');
                bp.add(outerPnl, 'Center');
                if nargin==6
                    south=Gui.Panel;
                    bp.add(south, 'South');
                    south.add(Gui.Label(advice));
                end
                allChb=Gui.CheckBox('All', true, [],[], @(h,e)doAll(h,e));
                north.add(allChb, 'West');
                outerPnl=bp;
                sortGui=SortGui([], allChb, 'All', north, strs, btns, innerPnl);
                sortGui.setAllChbText;
            end

            function innerCallback(h,e,i)
                if ~isempty(sortGui)
                    sortGui.setAllChbText;
                end
                feval(callback, h, e, i);
            end

            function doAll(h, e)
                state=h.isSelected;
                if state
                    Gui.ShowBusy(sortGui.dlg, 'Showing all', ...
                        'demoicon.gif', 7);
                else
                    Gui.ShowBusy(sortGui.dlg, 'Hiding all', 'hide.png', 5);
                end
                it=btns.iterator;
                idx=1;
                while it.hasNext
                    b=it.next;
                    if ~isempty(Gui.WindowAncestor(b))
                        if state ~= b.isSelected
                            b.setSelected(state);
                            try
                                priorActionCommand=b.getActionCommand;
                                b.setActionCommand(Plots.DOING_ALL);
                                ok=feval(callback, b, e, idx);
                                b.setActionCommand(priorActionCommand);
                                if ~ok
                                    break;
                                end
                            catch
                            end
                        end
                    end
                    idx=idx+1;
                end
                drawnow;
                Gui.HideBusy(sortGui.dlg)
            end

        end

        function [pnl, checkBoxes]=VerticalCheckBoxes(title, ...
                labels, dfltIdxs)
            jObj='javax.swing.JCheckBox';
            N=length(labels);
            pnl=Gui.SetTitledBorder(title);
            pnl.setLayout(javaObjectEDT('java.awt.GridLayout', N, 1));
            checkBoxes={};
            for i=1:N
                cbx=javaObjectEDT(jObj, labels{i});
                if ~isempty(find(dfltIdxs==i, 1))
                    cbx.setSelected(true);
                end
                checkBoxes{end+1}=cbx;
                pnl.add(cbx);
            end
        end

        function [pnl, bg, radioButtons]=Vertical(title, ...
                labels, dfltIdx)
            jObj='javax.swing.JRadioButton';
            N=length(labels);
            pnl=Gui.SetTitledBorder(title);
            pnl.setLayout(javaObjectEDT('java.awt.GridLayout', N, 1));
            radioButtons={};
            bg=javax.swing.ButtonGroup;
            for i=1:N
                rb=javaObjectEDT(jObj, labels{i});
                if dfltIdx==i
                    rb.setSelected(true);
                end
                radioButtons{end+1}=rb;
                pnl.add(rb);
                bg.add(rb);
            end
        end

        function [pnl, bg, radioButtons]=Horizontal(title, labels,...
                dfltIdx)
            jObj='javax.swing.JRadioButton';
            N=length(labels);
            if ~isempty(title)
                pnl=Gui.SetTitledBorder(title);
            else
                pnl=Gui.Panel;
            end

            bg=javax.swing.ButtonGroup;
            radioButtons={};
            for i=1:N
                rb=javaObjectEDT(jObj, labels{i});
                radioButtons{end+1}=rb;
                bg.add(rb);
                pnl.add(rb);
                if dfltIdx==i
                    rb.setSelected(true);
                end
            end
        end

        function Clear(bg)
            e=bg.getElements;
            while e.hasMoreElements
                rb=e.nextElement;
                if rb.isSelected
                    rb.setSelected(false);
                end
            end
        end

        function idx=Choice(bg)
            idx=1;
            e=bg.getElements;
            while e.hasMoreElements
                rb=e.nextElement;
                if rb.isSelected
                    return;
                end
                idx=idx+1;
            end
            idx=0;
        end

        function Disable(bg, idx)
            idx_=1;
            e=bg.getElements;
            while e.hasMoreElements
                rb=e.nextElement;
                if idx==idx_
                    rb.setEnabled(false);
                    return;
                end
                idx_=idx_+1;
            end
            idx_=0;
        end

        function [outerPnl, bg, checkBoxes, innerPnl]=Panel4(...
                javaArray, dfltIdxs, addScrollAfter, checkBox, gridBag)
            if ~checkBox
                checkBox=false;
                jObj='javax.swing.JRadioButton';
            else
                jObj='javax.swing.JCheckBox';
            end
            usingBag=false;
            N=javaArray.length;
            isJPanels=N>0 && isa(javaArray(1),'javax.swing.JPanel');
            innerPnl=javaObjectEDT('javax.swing.JPanel');
            innerPnl.setLayout(javax.swing.BoxLayout(innerPnl, 1));
            columns=1;
            if checkBox
                bg=javaObjectEDT('javax.swing.JPanel');
                if ~isJPanels
                    horizontal=true;
                    for i=1:N
                        if javaArray(i).length>2
                            horizontal=false;
                            break;
                        else
                        end
                    end
                else
                    columns=javaArray(1).getComponentCount;
                    horizontal=false;
                end
                if ~horizontal
                    if ~gridBag && ~isJPanels
                        innerPnl.setLayout(javaObjectEDT(...
                            'java.awt.GridLayout', N, 1));
                    else
                        innerPnl.setLayout(javaObjectEDT(...
                            'java.awt.GridBagLayout'));
                        gbc=javaObjectEDT('java.awt.GridBagConstraints');
                        %gbc.weighty=1;
                        gbc.anchor=gbc.WEST;
                        gbc.fill=0;
                        usingBag=true;
                    end
                else
                    innerPnl.setLayout(javaObjectEDT(...
                        'java.awt.GridLayout', 1, N));
                end
            else
                bg=javax.swing.ButtonGroup;
            end
            checkBoxes=java.util.ArrayList;
            idx=1;
            scrollToHeight=0;
            first=true;
            if ~checkBox || ~gridBag
                %javaArray is char or strings
                for i=1:N
                    if isJPanels
                        rb=javaObjectEDT(jObj, javaArray(i).getComponent(0).getText);
                    else
                        rb=javaObjectEDT(jObj, javaArray(i));
                    end
                    rObj=rb;
                    checkBoxes.add(rb);
                    bg.add(rb);
                    if ~usingBag
                        innerPnl.add(rObj);
                    else
                        gbc.gridy=i-1;
                        gbc.gridx=0;
                        innerPnl.add(rObj, gbc);
                    end
                    for j=2:columns
                        gbc.gridx=j-1;
                        gbc.weightx=.2;
                        jp=javaObjectEDT('javax.swing.JPanel');
                        jp.add(javaArray(i).getComponent(j-1));
                        innerPnl.add(jp, gbc);
                    end
                    if ~isempty(find(dfltIdxs==idx, 1))
                        rb.setSelected(true);
                        if first
                            d=rb.getPreferredSize;
                            scrollToHeight=d.height*(i-1);
                            first=false;
                        end
                    end
                    idx=idx+1;
                end
            else
                try
                    for i=1:N
                        if isJPanels
                            rb=javaObjectEDT(jObj, javaArray(i).getComponent(0).getText);
                        else
                            rb=javaObjectEDT(jObj, javaArray(i));
                        end
                        checkBoxes.add(rb);
                        gbc.gridy=i-1;
                        gbc.gridx=0;
                        gbc.weightx=.8;
                        fprintf('y=%d x=%d wx=%d wy=%d fill=%d\n', ...
                            gbc.gridy, gbc.gridx, gbc.weighty, ...
                            gbc.weightx, gbc.fill);
                        innerPnl.add(rb, gbc);
                        for j=2:columns
                            gbc.gridx=j-1;
                            gbc.weightx=.2;
                            jp=javaObjectEDT('javax.swing.JPanel');
                            jp.add(javaArray(i).getComponent(j-1));
                            innerPnl.add(jp, gbc);
                        end
                        if ~isempty(find(dfltIdxs==idx, 1))
                            rb.setSelected(true);
                            if first
                                d=rb.getPreferredSize;
                                scrollToHeight=d.height*(i-1);
                                first=false;
                            end
                        end
                        idx=idx+1;
                    end
                catch ex
                    disp(['If javaArray contains JPanels then '...
                        'Radio.Panel2' newline 'expects component 0'...
                        ' to be JLabel and any 2nd component!' newline]);
                    ex.getReport
                end
            end
            outerPnl=javaObjectEDT('javax.swing.JPanel');
            outerPnl.setBorder(...
                javax.swing.BorderFactory.createEmptyBorder(1,8,1,4));
            if addScrollAfter<N
                d=innerPnl.getPreferredSize;
                d.width=floor(d.width*1.1);
                d.height=floor( (addScrollAfter/N)*d.height);
                innerPnl.setBorder([]);
                scroll=javaObjectEDT('javax.swing.JScrollPane', innerPnl);
                if scrollToHeight>0
                    scroll.getViewport.setViewPosition(...
                        java.awt.Point(0, scrollToHeight));
                end
                scroll.setPreferredSize(d);
                outerPnl=scroll;
            else
                cb=javax.swing.BorderFactory.createCompoundBorder(...
                    javax.swing.BorderFactory.createBevelBorder(...
                    javax.swing.border.BevelBorder.RAISED), ...
                    javax.swing.BorderFactory.createEmptyBorder(8, 8, 8, 8));
                innerPnl.setBorder(cb);
                outerPnl.add(innerPnl);
            end
        end

    end
end


text:
classdef GradientVariationWithBoundaryComputer < handle

    properties (Access = public)
       error
       mean
       desv
    end

    properties (Access = private)
        levelSet
        curvature
        rPerimeter
        backgroundMesh
        boundaryMesh
        domainLength
        circleCase
        filePlotName
        circleMesh
        unfittedMesh
        integrator

        gradientCircunf
        nEpsilon
    end

    methods (Access = public)

        function obj = GradientVariationWithBoundaryComputer(cParams)
            obj.init(cParams)
        end

        function compute(obj)
            obj.createIntegrator();
            for iEpsilon = 1:obj.nEpsilon
                obj.computeGradientInCircunference(iEpsilon);
                obj.computeVariables(iEpsilon);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.filePlotName   = cParams.filePlotName;
            obj.backgroundMesh = cParams.backgroundMesh;
            obj.circleMesh     = cParams.circleMesh;
            obj.levelSet       = cParams.levelSet;
            obj.curvature      = cParams.curvature;
            obj.rPerimeter     = cParams.regularizedPerimeter;
            obj.domainLength   = cParams.domainLength;
            obj.circleCase     = cParams.circleCase;
            obj.nEpsilon       = size(obj.rPerimeter.epsilons,2);
        end

        function computeVariables(obj,iepsilon)
           int = obj.integrator;
           M = int.computeLHS();
           k = obj.curvature;
           g = obj.gradientCircunf;
           obj.error(iepsilon) = (g-k)'*M*(g-k)/(k*sum(M(:))*k);
           obj.mean(iepsilon) = sum(int.integrate(g))/sum(int.integrate(ones(size(g))));
           obj.desv(iepsilon) = (g-obj.mean(iepsilon))'*M*(g-obj.mean(iepsilon))/(sum(M(:)));
        end

        function createIntegrator(obj)
           m = obj.circleMesh.mesh;
           m = m.computeCanonicalMesh();
           s.type = 'SIMPLE';
           s.mesh = m;
           s.npnod = m.nnodes;
           s.globalConnec = m.connec;
           obj.integrator = Integrator.create(s);
        end

        function computeGradientInCircunference(obj,iepsilon)
            nCell = obj.circleMesh.mesh.nelem;
            gB = zeros(nCell,1);
            gD = obj.rPerimeter.perimetersGradient(:,iepsilon);
            switch obj.circleCase
                case 'interior'
                    for icell = 1:nCell
                        gNodalInCell  = obj.computeFnodalInCell(icell,gD);
                        xPos          = obj.computeXpos(icell);
                        gInXpos       = obj.interpolateValue(xPos,gNodalInCell);
                        gB(icell) = gInXpos;
                    end
                case 'exterior'
                    m = obj.circleMesh;
                    gB = gD(m.nodesInBoxFaces,1);
            end
            obj.gradientCircunf = gB;
        end

        function fNodalInCell = computeFnodalInCell(obj,icell,fNodal)
            m = obj.circleMesh;
            cellGlobal   = m.cellContainingSubcell(icell);
            nodeTriangle = obj.backgroundMesh.connec(cellGlobal,:);
            fNodalInCell = fNodal(nodeTriangle(:),:);
        end

        function x1pos = computeXpos(obj,icell)
            m = obj.circleMesh;
            x1pos = m.xCoordsIso(:,1,icell);
            x1pos = squeeze(x1pos);
        end

        function fInterp = interpolateValue(obj,xpos,fNodal)
            interp = Interpolation.create(obj.backgroundMesh,'LINEAR');
            interp.computeShapeDeriv(xpos);
            shape = interp.shape;
            fInterp = 0;
            for inode = 1:length(shape)
                fInterp = fInterp + fNodal(inode)*shape(inode);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function guv = evaluate_log_posterior_grad(this, uv)
%EVALUATE_LOG_POSTERIOR computes the gradient of the log-posterior
%   (negative energy) wrt the flow fields UV
%   Actually only proportional to the log posterior since the variance of neither the
%   spatial nor the data terms is considered
%
%   This is a member function of the class 'ba_optical_flow'.
%
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

% Spatial term
S       = this.spatial_filters;
gu1     = zeros(size(uv,1), size(uv,2));
gv1     = gu1;

for i = 1:length(S)

    u_ = conv2(uv(:,:,1), S{i}, 'valid');
    v_ = conv2(uv(:,:,2), S{i}, 'valid');

    Si = reshape(S{i}(end:-1:1), size(S{i}));

    if isa(this.rho_spatial_u{i}, 'robust_function')
        u_ = -reshape(deriv(this.rho_spatial_u{i}, u_(:)), size(u_));
        v_ = -reshape(deriv(this.rho_spatial_v{i}, v_(:)), size(v_));
    elseif isa(this.rho_spatial_u{i}, 'gsm_density')
        u_ = reshape(evaluate_log_grad(this.rho_spatial_u{i}, u_(:)'), size(u_));
        v_ = reshape(evaluate_log_grad(this.rho_spatial_v{i}, v_(:)'), size(v_));
    else
        error('evaluate_log_posterior: unknown rho function!');
    end;

    gu1 = gu1+conv2(u_, Si, 'full');
    gv1 = gv1+conv2(v_, Si, 'full');
end;

gu2     = zeros(size(uv,1), size(uv,2));
gv2     = gu2;

% Data term
[It Ix Iy] = partial_deriv(this.images, uv, this.interpolation_method, this.deriv_filter, this.blend);

if isa(this.rho_data, 'robust_function')
    temp   = -reshape(deriv(this.rho_data, It(:)), size(It));

elseif isa(this.rho_data, 'gsm_density')
    temp   = reshape(evaluate_log_grad(this.rho_data, It(:)'), size(It));
else
    error('evaluate_log_posterior: unknown rho function!');
end;

gu2     = temp.*Ix;
gv2     = temp.*Iy;

guv = cat(3, gu2+this.lambda*gu1, gv2+this.lambda*gv1);

author: kristinbranson

repo: JAABA


text:
function [Eft, Covft, ljpyt, Eyt, Covyt] = gp_jpred(gp, x, y, varargin)
%GP_JPRED  Joint predictions with Gaussian process
%
%  Description
%    [EFT, COVFT] = GP_JPRED(GP, X, Y, XT, OPTIONS)
%    takes a GP structure together with matrix X of training
%    inputs and vector Y of training targets, and evaluates the
%    joint predictive distribution at test inputs XT. Returns a posterior
%    mean EFT and covariance COVFT of latent variables.
%
%        Eft =  E[f | xt,x,y,th]  = K_fy*(Kyy+s^2I)^(-1)*y
%      Covft = Var[f | xt,x,y,th] = K_fy - K_fy*(Kyy+s^2I)^(-1)*K_yf.
%
%    Each row of X corresponds to one input vector and each row of
%    Y corresponds to one output vector.
%
%    [EFT, COVFT, LJPYT] = GP_JPRED(GP, X, Y, XT, 'yt', YT, ...)
%    returns also logarithm of the predictive joint density LJPYT of
%    the observations YT at test input locations XT. This can be
%    used for example in the cross-validation. Here Y has to be
%    vector.
%
%    [EFT, COVFT, LJPYT, EYT, COVYT] = GP_JPRED(GP, X, Y, XT, 'yt', YT, ...)
%    returns also the posterior predictive mean and covariance.
%
%    [EF, COVF, LJPY, EY, VARY] = GP_JPRED(GP, X, Y, OPTIONS)
%    evaluates the predictive distribution at training inputs X
%    and logarithm of the predictive density PY of the training
%    observations Y.
%
%    OPTIONS is optional parameter-value pair
%      predcf - an index vector telling which covariance functions are
%                 used for prediction. Default is all (1:gpcfn).
%                 See additional information below.
%      tstind - a vector/cell array defining, which rows of X belong
%               to which training block in *IC type sparse models.
%               Default is []. In case of PIC, a cell array
%               containing index vectors specifying the blocking
%               structure for test data. IN FIC and CS+FIC a
%               vector of length n that points out the test inputs
%               that are also in the training set (if none, set
%               TSTIND = [])
%      yt     - optional observed yt in test points (see below)
%
%    NOTE! In case of FIC and PIC sparse approximation the
%    prediction for only some PREDCF covariance functions is just
%    an approximation since the covariance functions are coupled in
%    the approximation and are not strictly speaking additive
%    anymore.
%
%    For example, if you use covariance such as K = K1 + K2 your
%    predictions Ef1 = gp_pred(GP, X, Y, X, 'predcf', 1) and Ef2 =
%    gp_pred(gp, x, y, x, 'predcf', 2) should sum up to Ef =
%    gp_pred(gp, x, y, x). That is Ef = Ef1 + Ef2. With FULL model
%    this is true but with FIC and PIC this is true only
%    approximately. That is Ef \approx Ef1 + Ef2.
%
%    With CS+FIC the predictions are exact if the PREDCF covariance
%    functions are all in the FIC part or if they are CS
%    covariances.
%
%    NOTE! When making predictions with a subset of covariance
%    functions with FIC approximation the predictive variance can
%    in some cases be ill-behaved i.e. negative or
%    unrealistically small. This may happen because of the
%    approximative nature of the prediction.
%
%  See also
%    GP_SET, GP_OPTIM, DEMO_REGRESSION*
%

% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2008 Jouni Hartikainen
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2010,2012 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.

if iscell(gp) || numel(gp.jitterSigma2)>1 || isfield(gp,'latent_method')
  % use an inference specific method
  if iscell(gp)
    fh_pred=@gpia_jpred;
  elseif numel(gp.jitterSigma2)>1
    fh_pred=@gpmc_jpred;
  elseif isfield(gp,'latent_method')
    fh_pred=gp.fh.jpred;
  else
    error('Logical error by coder of this function!')
  end
  switch nargout
    case 1
      [Eft] = fh_pred(gp, x, y, varargin{:});
    case 2
      [Eft, Covft] = fh_pred(gp, x, y, varargin{:});
    case 3
      [Eft, Covft, ljpyt] = fh_pred(gp, x, y, varargin{:});
    case 4
      [Eft, Covft, ljpyt, Eyt] = fh_pred(gp, x, y, varargin{:});
    case 5
      [Eft, Covft, ljpyt, Eyt, Covyt] = fh_pred(gp, x, y, varargin{:});
  end
  return
end

ip=inputParser;
ip.FunctionName = 'GP_JPRED';
ip.addRequired('gp',@isstruct);
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                 isvector(x) && isreal(x) && all(isfinite(x)&x>0))
ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                 (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
ip.parse(gp, x, y, varargin{:});
xt=ip.Results.xt;
yt=ip.Results.yt;
predcf=ip.Results.predcf;
tstind=ip.Results.tstind;
if isempty(xt)
  xt=x;
  if isempty(tstind)
    if iscell(gp)
      gptype=gp{1}.type;
    else
      gptype=gp.type;
    end
    switch gptype
      case {'FULL' 'VAR' 'DTC' 'SOR'}
        tstind = [];
      case {'FIC' 'CS+FIC'}
        tstind = 1:size(x,1);
      case 'PIC'
        if iscell(gp)
          tstind = gp{1}.tr_index;
        else
          tstind = gp.tr_index;
        end
    end
  end
  if isempty(yt)
    yt=y;
  end
end

tn = size(x,1);
if isfield(gp.lik, 'nondiagW') && ~ismember(gp.lik.type, {'LGP' 'LGPC'})
  switch gp.lik.type
    case {'Softmax', 'Multinom'}
      nout=size(y(:),1)/tn;
    otherwise
      if ~isfield(gp.lik,'xtime') && size(y,1)~=size(x,1)
        y=reshape(y,size(x,1),size(y,1)/size(x,1));
      end
      nout=length(gp.comp_cf);

      % Indices for looping over latent processes
      if ~isfield(gp.lik, 'xtime')
        nl=[0 repmat(n,1,nout)];
        nl=cumsum(nl);
      else
        xtime=gp.lik.xtime;
        ntime = size(xtime,1);
        n=tn-ntime;
        nl=[0 ntime n];
        nl=cumsum(nl);
      end
  end
  y=reshape(y,tn,nout);

  if isfield(gp, 'comp_cf')  % own covariance for each ouput component
    multicf = true;
    if length(gp.comp_cf) ~= nout
      error('GP2_PRED: the number of component vectors in gp.comp_cf must be the same as number of outputs or latent processes.')
    end
    if ~isempty(predcf)
      if ~iscell(predcf) || length(predcf)~=nout
        error(['GP2_PRED: if own covariance for each output or latent process component is used,'...
          'predcf has to be cell array and contain nout (vector) elements.   '])
      end
    else
      predcf = gp.comp_cf;
    end
  else
    multicf = false;
    for i1=1:nout
      predcf2{i1} = predcf;
    end
    predcf=predcf2;
  end
end

if nargout > 2 && isempty(yt)
  ljpyt=[];
end

% Evaluate this if sparse model is used
switch gp.type
  case 'FULL'

    if ~isfield(gp.lik, 'nondiagW') && ismember(gp.lik.type, {'LGP', 'LGPC'})
        %evaluate a = C\y;
        % -------------------
        [c, C]=gp_trcov(gp,x);

        if issparse(C)
          LD = ldlchol(C);
          a = ldlsolve(LD,y);
        elseif isempty(C)
          C=0;
          L=[];
          a = zeros(length(y),1);
        else
          L = chol(C)';
          a = L'\(L\y);
        end

        % evaluate K*a
        % -------------------
        K=gp_cov(gp,x,xt,predcf);
        Eft = K'*a;

        if  isfield(gp,'meanf')
          if issparse(C)
            [RB RAR] = mean_jpredf(gp,x,xt,K,LD,a,'gaussian',[]);    % terms with non-zero mean -prior
          else
            [RB RAR] = mean_jpredf(gp,x,xt,K,L,a,'gaussian',[]);    % terms with non-zero mean -prior
          end
          Eft = Eft + RB;
        end

        % Evaluate variance
        % Vector of diagonal elements of covariance matrix
        if nargout > 1

          V = gp_trcov(gp,xt,predcf);
          if issparse(C)
            Covft = (V - (K'*ldlsolve(LD,K)));
          else
            v = L\K;
            Covft = (V - (v'*v));
          end

          % If there are specified mean functions
          if  isfield(gp,'meanf')
            Covft = Covft + RAR;
          end
        end

        if nargout > 2
          % Scale mixture model in lik_smt is a special case
          % handle it separately
          if ~strcmp(gp.lik.type, 'lik_smt')
            % normal case
            [V, Cv] = gp_trvar(gp,xt,predcf);
            Eyt = Eft;
            apu = Cv - V;
            Covyt = Covft + diag(apu); % Utilize the Covft calculated above (faster!?) dimensions did not match here earlier!
            % Log joint predictive density
            if ~isempty(yt)
              ljpyt = mnorm_lpdf(yt', Eyt', Covyt);
            end
          else
            % scale mixture case
            nu = gp.lik.nu;             % Not working at the moment. Probably.
            sigma2 = gp.lik.tau2.*gp.lik.alpha.^2;
            sigma = sqrt(sigma2);

            Eyt = Eft;
            Covyt = (nu./(nu-2).*sigma2);

            for i2 = 1:length(Eft)
              mean_app = Eft(i2);
              sigm_app = sqrt(Covft(i2));

              pd = @(f) t_pdf(yt(i2), nu, f, sigma).*norm_pdf(f,Eft(i2),sqrt(Covft(i2)));
              pyt(i2) = quadgk(pd, mean_app - 12*sigm_app, mean_app + 12*sigm_app);
            end
          end
        end
    else
      % Likelihoods with non-diagonal Hessian
      L = zeros(tn,tn,nout);
      ntest=size(xt,1);
      K_nf = zeros(ntest,tn,nout);
      if multicf
        for i1=1:nout
          [tmp,C] = gp_trcov(gp, x, gp.comp_cf{i1});
          L(:,:,i1) = chol(C)';
          K_nf(:,:,i1) = gp_cov(gp,xt,x,predcf{i1});
        end
      else
        for i1=1:nout
          [tmp,C] = gp_trcov(gp, x);
          L(:,:,i1) = chol(C)';
          K_nf(:,:,i1) = gp_cov(gp,xt,x,predcf{i1});
        end
      end

      Eft = zeros(ntest,nout);
      for i1=1:nout
        Eft(:,i1) = K_nf(:,:,i1)*(L(:,:,i1)'\(L(:,:,i1)\y(:,i1)));
      end

      Covft = zeros(ntest*nout,ntest*nout);
      if nargout > 1
        if multicf
          for i1=1:nout
            v = L(:,:,i1)\K_nf(:,:,i1)';
            V = gp_trcov(gp,xt,gp.comp_cf{i1});
            Covft((i1-1)*ntest+1:i1*ntest,(i1-1)*ntest+1:i1*ntest) = V - v'*v;
          end
        else
          for i1=1:nout
            v = L(:,:,i1)\K_nf(:,:,i1)';
            V = gp_trcov(gp,xt,predcf{i1});
            Covft((i1-1)*ntest+1:i1*ntest,(i1-1)*ntest+1:i1*ntest) = V - v'*v;
          end
        end
      end
      if nargout > 2
        % normal case
        Eyt = Eft;
        Covyt=Covft;
        ljpyt=0;
        for i1=1:nout
          [V, Cv] = gp_trvar(gp,xt,predcf{i1});
          % Diagonal indices
          i2=(i1-1)*(nout*ntest^2+ntest)+1:nout*ntest+1:i1*(nout*ntest^2+ntest);
          Covyt(i2) = Covyt(i2) + Cv' - V';
          if ~isempty(yt)
            ljpyt = ljpyt + mnorm_lpdf(yt(:,i1)', Eyt(:,i1)', Covyt((i1-1)*ntest+1:i1*ntest,(i1-1)*ntest+1:i1*ntest));
          end
        end
        Eyt=Eyt(:);
      end
      Eft=Eft(:);
    end

  case 'FIC'
    % Check the tstind vector
    if nargin > 5
        if ~isempty(tstind) && length(tstind) ~= size(x,1)
            error('tstind (if provided) has to be of same lenght as x.')
        end
    else
        tstind = [];
    end

    u = gp.X_u;
    m = size(u,1);
    % Turn the inducing vector on right direction
    if size(u,2) ~= size(x,2)
        u=u';
    end
    % Calculate some help matrices
    [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
    K_fu = gp_cov(gp, x, u);   % f x u
    K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
    K_nu = gp_cov(gp,xt,u);       % n x u
    Luu = chol(K_uu,'lower');

    % Evaluate the Lambda (La) for specific model
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    B=Luu\(K_fu');
    Qv_ff=sum(B.^2)';
    Lav = Cv_ff-Qv_ff;   % 1 x f, Vector of diagonal elements
                         % iLaKfu = diag(inv(Lav))*K_fu = inv(La)*K_fu
    iLaKfu = zeros(size(K_fu));  % f x u,
    n=size(x,1)
    for i=1:n
        iLaKfu(i,:) = K_fu(i,:)./Lav(i);  % f x u
    end
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;

    L = iLaKfu/chol(A);
    p = y./Lav - L*(L'*y);

    % Prediction matrices formed with only subset of cf's.
    if ~isempty(predcf)
        K_fu = gp_cov(gp, x, u, predcf);   % f x u
        K_uu = gp_trcov(gp, u, predcf);     % u x u, noiseles covariance K_uu
        K_nu = gp_cov(gp,xt,u,predcf);       % n x u
    end
    Eft = K_nu*(K_uu\(K_fu'*p));

    % if the prediction is made for training set, evaluate Lav also for prediction points
    if ~isempty(tstind)
        [Kv_ff, Cv_ff] = gp_trvar(gp, xt(tstind,:), predcf);
        Luu = chol(K_uu,'lower');
        B=Luu\(K_fu');
        Qv_ff=sum(B.^2)';
        Lav2 = zeros(size(Eft));
        Lav2(tstind) = Kv_ff-Qv_ff;
        Eft(tstind) = Eft(tstind) + Lav2(tstind).*p;
    end

    if nargout > 1
        [Knn_v, Cnn_v] = gp_trvar(gp,xt,predcf);
        Luu = chol(K_uu, 'lower');
        B=Luu\(K_fu');
        B2=Luu\(K_nu');
        Lav_n = Knn_v - sum(B2.^2)';
        BL = B*L;

        if isempty(tstind)
            error('tstind not provided.')
        end

        % if the prediction is made for training set, evaluate Lav also for prediction points
        if ~isempty(tstind)
          K = B'*B2 + diag(Lav_n);
          K2 = B2'*B2 + diag(Lav_n);
          C = B'*B + diag(Lav);

          L = chol(C,'lower');
          a = L'\(L\y);
          v = L\K;

          Covft = K2-v'*v;
        end
    end

    if nargout > 2
        Eyt = Eft;
        Covyt = Covft + diag(Cnn_v) - diag(Knn_v);
        if ~isempty(yt)
          ljpyt = mnorm_lpdf(yt', Eyt', Covyt);
        end
    end

  case {'PIC' 'PIC_BLOCK'}
    u = gp.X_u;
    ind = gp.tr_index;
    if size(u,2) ~= size(x,2)
        u=u';
    end

    % Calculate some help matrices
    [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
    K_fu = gp_cov(gp, x, u);         % f x u
    K_nu = gp_cov(gp, xt, u);         % n x u
    K_uu = gp_trcov(gp, u);    % u x u, noiseles covariance K_uu
    Luu = chol(K_uu)';

    % Evaluate the Lambda (La) for specific model
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    B=Luu\K_fu';
    iLaKfu = zeros(size(K_fu));  % f x u
    for i=1:length(ind)
        Qbl_ff = B(:,ind{i})'*B(:,ind{i});
        [Kbl_ff, Cbl_ff] = gp_trcov(gp, x(ind{i},:));
        La{i} = Cbl_ff - Qbl_ff;
        iLaKfu(ind{i},:) = La{i}\K_fu(ind{i},:);
    end
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;            % Ensure symmetry
    L = iLaKfu/chol(A);

    tyy = y;
    % From this on evaluate the prediction
    % See Snelson and Ghahramani (2007) for details
    p=iLaKfu*(A\(iLaKfu'*tyy));
    for i=1:length(ind)
        p2(ind{i},:) = La{i}\tyy(ind{i},:);
    end
    p= p2-p;

    % Prediction matrices formed with only subsetof cf's.
    if ~isempty(predcf)
        K_fu = gp_cov(gp, x, u, predcf);        % f x u
        K_nu = gp_cov(gp, xt, u, predcf);         % n x u
        K_uu = gp_trcov(gp, u, predcf);          % u x u, noiseles covariance K_uu
    end

    iKuuKuf = K_uu\K_fu';
    w_bu=zeros(length(xt),length(u));
    w_n=zeros(length(xt),1);
    B2=Luu\(K_nu');
    for i=1:length(ind)
        w_bu(tstind{i},:) = repmat((iKuuKuf(:,ind{i})*p(ind{i},:))', length(tstind{i}),1);
        K_nf = gp_cov(gp, xt(tstind{i},:), x(ind{i},:),predcf);              % n x u
        w_n(tstind{i},:) = K_nf*p(ind{i},:);
        Qbl_ff = B2(:,tstind{i})'*B2(:,tstind{i});
        [Kbl_ff, Cbl_ff] = gp_trcov(gp, xt(tstind{i},:));
        La2{i} = Kbl_ff - Qbl_ff;
    end

    Eft = K_nu*(iKuuKuf*p) - sum(K_nu.*w_bu,2) + w_n;

    if nargout > 1
        % Form iLaKfu again if a subset of cf's is used for making predictions
        if ~isempty(predcf)
            iLaKfu = zeros(size(K_fu));  % f x u
            for i=1:length(ind)
                iLaKfu(ind{i},:) = La{i}\K_fu(ind{i},:);
            end
        end

%         kstarstar = gp_trvar(gp, xt, predcf);
%         KnuiKuu = K_nu/K_uu;
%         KufiLaKfu = K_fu'*iLaKfu;
%         QnfL = KnuiKuu*(K_fu'*L);
%         Covft1 = zeros(size(xt,1),size(xt,1));
%         Covft2 = zeros(size(xt,1),size(xt,1));
%         Covft3 = zeros(size(xt,1),size(xt,1));
%         C = B'*B;
%         KKnn = gp_trcov(gp,xt,predcf);
%         Knn = B2'*B2;
%         Knf = B2'*B;
%         KKnf = gp_cov(gp,xt,x,predcf);
%         for i=1:length(ind)
%             KubiLaKbu = K_fu(ind{i},:)'/La{i}*K_fu(ind{i},:);
%             nonblock = KufiLaKfu - KubiLaKbu;
%             Covft1(tstind{i}) = diag(KnuiKuu(tstind{i},:)*nonblock*KnuiKuu(tstind{i},:)');
%
%             Knb = gp_cov(gp, xt(tstind{i},:), x(ind{i},:), predcf);
%             Covft2(tstind{i}) = diag(Knb/La{i}*Knb');
%
%             KnbL = Knb*L(ind{i},:);
%             QnbL = KnuiKuu(tstind{i},:)*(K_fu(ind{i},:)'*L(ind{i},:));
%             %Covft3(tstind{i}) = sum(QnfL(tstind{i},:) - QnbL + KnbL,2);
%             Covft3(tstind{i}) = diag((QnfL(tstind{i},:) - QnbL + KnbL)*(QnfL(tstind{i},:) - QnbL + KnbL)');
%             C(ind{i},ind{i}) = C(ind{i},ind{i}) + La{i};
%             Knn(ind{i},ind{i}) = Knn(ind{i},ind{i}) + KKnn(tstind{i},tstind{i}) - B2(:,tstind{i})'*B2(:,tstind{i});
%             Knf(tstind{i},ind{i}) = Knf(tstind{i},ind{i}) + KKnf(tstind{i},ind{i}) - B2(:,tstind{i})'*B(:,ind{i});
%         end
%         L = chol(C,'lower');
%         v = L\Knf;
%         Covft = Knn-v'*v;
% %         Covft = kstarstar - (Covft1 + Covft2 - Covft3);       %MUUTETTU
%     end

        B2=Luu\(K_nu');
        C = B'*B;
        KKnn = gp_trcov(gp,xt,predcf);
        Knn = B2'*B2;
        Knf = B2'*B;
        KKnf = gp_cov(gp,xt,x,predcf);
        for i=1:length(ind)
            if (size(ind{i},1) ~= size(tstind{i},1))
                error('Testing and training block cell array vectors differ in size.');     % gp.tr_index and tstind
            end
            C(ind{i},ind{i}) = C(ind{i},ind{i}) + La{i};
            Knn(ind{i},ind{i}) = Knn(ind{i},ind{i}) + KKnn(tstind{i},tstind{i}) - B2(:,tstind{i})'*B2(:,tstind{i});
            Knf(tstind{i},ind{i}) = Knf(tstind{i},ind{i}) + KKnf(tstind{i},ind{i}) - B2(:,tstind{i})'*B(:,ind{i});
        end

        L = chol(C)';

        v = L\Knf';
        Covft = (Knn) - (v'*v);
    end

    if nargout > 2
        Eyt = Eft;
        [Knn_v, Cnn_v] = gp_trvar(gp,xt,predcf);
        Covyt = Covft + diag(Cnn_v) - diag(Knn_v);
        if ~isempty(yt)
          ljpyt = mnorm_lpdf(yt', Eyt', Covyt);
        end
    end
  case 'CS+FIC'
    % Here tstind = 1 if the prediction is made for the training set
    if nargin > 5
        if ~isempty(tstind) && length(tstind) ~= size(x,1)
            error('tstind (if provided) has to be of same lenght as x.')
        end
    else
        tstind = [];
    end

    n = size(x,1);
    n2 = size(xt,1);

    u = gp.X_u;
    m = size(u,1);
    ncf = length(gp.cf);

    % Indexes to all non-compact support and compact support covariances.
    cf1 = [];
    cf2 = [];
    % Indexes to non-CS and CS covariances, which are used for predictions
    predcf1 = [];
    predcf2 = [];

    % Loop through all covariance functions
    for i = 1:ncf
        % Non-CS covariances
        if ~isfield(gp.cf{i},'cs')
            cf1 = [cf1 i];
            % If used for prediction
            if ~isempty(find(predcf==i))
                predcf1 = [predcf1 i];
            end
        % CS-covariances
        else
            cf2 = [cf2 i];
            % If used for prediction
            if ~isempty(find(predcf==i))
                predcf2 = [predcf2 i];
            end
        end
    end
    if isempty(predcf1) && isempty(predcf2)
        predcf1 = cf1;
        predcf2 = cf2;
    end

    % Determine the types of the covariance functions used
    % in making the prediction.
    if ~isempty(predcf1) && isempty(predcf2)       % Only non-CS covariances
        ptype = 1;
        predcf2 = cf2;
    elseif isempty(predcf1) && ~isempty(predcf2)   % Only CS covariances
        ptype = 2;
        predcf1 = cf1;
    else                                           % Both non-CS and CS covariances
        ptype = 3;
    end

    % First evaluate needed covariance matrices
    % v defines that parameter is a vector
    [Kv_ff, Cv_ff] = gp_trvar(gp, x, cf1);  % f x 1  vector
    K_fu = gp_cov(gp, x, u, cf1);         % f x u
    K_uu = gp_trcov(gp, u, cf1);    % u x u, noiseles covariance K_uu
    K_uu = (K_uu+K_uu')./2;     % ensure the symmetry of K_uu

    Luu  = chol(K_uu)';
    K_nu = gp_cov(gp, xt, u, cf1);         % n x u

    % Evaluate the Lambda (La)
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    B=Luu\(K_fu');       % u x f
    Qv_ff=sum(B.^2)';
    Lav = Cv_ff-Qv_ff;   % f x 1, Vector of diagonal elements

    K_cs = gp_trcov(gp,x,cf2);
    Kcs_nf = gp_cov(gp, xt, x, predcf2);
    Kcs_nn = gp_trcov(gp, xt, predcf2);
    La = sparse(1:tn,1:tn,Lav,tn,tn) + K_cs;

    iLaKfu = La\K_fu;
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;     % Ensure symmetry
    L = iLaKfu/chol(A);

    p = La\y - L*(L'*y);

    %p2 = y./Lav - iLaKfu*(A\(iLaKfu'*y));
    %    Knf = K_nu*(K_uu\K_fu');

    K_fu = gp_cov(gp, x, u, predcf1);       % f x u
    K_uu = gp_trcov(gp, u, predcf1);         % u x u, noiseles covariance K_uu
    K_uu = (K_uu+K_uu')./2;                  % ensure the symmetry of K_uu
    K_nu = gp_cov(gp, xt, u, predcf1);        % n x u

    % Calculate the predictive mean according to the type of
    % covariance functions used for making the prediction
    if ptype == 1
        Eft = K_nu*(K_uu\(K_fu'*p));
    elseif ptype == 2
        Eft = Kcs_nf*p;
    else
        Eft = K_nu*(K_uu\(K_fu'*p)) + Kcs_nf*p;
    end

    % evaluate also Lav2 if the prediction is made for training set
    if ~isempty(tstind)
        [Kv_ff, Cv_ff] = gp_trvar(gp, xt(tstind,:), predcf1);
        Luu = chol(K_uu)';
        B=Luu\(K_fu');
        Qv_ff=sum(B.^2)';
        Lav2 = zeros(size(Eft));
        Lav2(tstind) = Kv_ff-Qv_ff;
    end

    % Add also Lav2 if the prediction is made for training set
    % and non-CS covariance function is used for prediction
    if ~isempty(tstind) && (ptype == 1 || ptype == 3)
        Eft(tstind) = Eft(tstind) + Lav2(tstind).*p;
    end

    if nargout > 1
        [Knn_v, Cnn_v] = gp_trvar(gp,xt,predcf);
        Luu = chol(K_uu, 'lower');
        B=Luu\(K_fu');
        B2=Luu\(K_nu');
        iLaKfu = La\K_fu;

        % Calculate the predictive variance according to the type
        % covariance functions used for making the prediction
        if ptype == 1 || ptype == 3
            % FIC part of the covariance
%             Covft = Knn_v - sum(B2'.*(B*(La\B')*B2)',2) +
%             sum((K_nu*(K_uu\(K_fu'*L))).^2, 2);             % MUUTETTU
            tmpmatr = (K_nu*(K_uu\(K_fu'*L)));
            Covft = B2'*B2 + Kcs_nn + sparse(1:n2,1:n2,Knn_v - sum(B2.^2)',n2,n2) - B2'*(B*(La\B')*B2) + tmpmatr*tmpmatr';
%             Covft = B2'*B2 + Kcs_nn - B2'*(B*(La\B')*B2) + tmpmatr*tmpmatr';

            % Add Lav2 if the prediction is made for the training set
            if  ~isempty(tstind)
                % Non-CS covariance
                if ptype == 1
                    Kcs_nf = sparse(tstind,1:n,Lav2(tstind),n2,n);
                % Non-CS and CS covariances
                else
                    Kcs_nf = Kcs_nf + sparse(tstind,1:n,Lav2(tstind),n2,n);
                end
                 % Add Lav2 and possibly Kcs_nf
                 tmpmatr = Kcs_nf/chol(La);
                 Covft = Covft - tmpmatr*tmpmatr';
                 tmpmatr = Kcs_nf*L;
                 Covft = Covft + tmpmatr*tmpmatr';
                 Covft = Covft - 2.*(Kcs_nf*iLaKfu)*(K_uu\K_nu') + 2.*(Kcs_nf*L)*(L'*K_fu*(K_uu\K_nu'));
            % In case of both non-CS and CS prediction covariances add
            % only Kcs_nf if the prediction is not done for the training set
            elseif ptype == 3
                tmpmatr = Kcs_nf/chol(La);
                Covft = Covft - tmpmatr*tmpmatr';
                tmpmatr = Kcs_nf*L;
                Covft = Covft + tmpmatr*tmpmatr';
                Covft = Covft - 2*(Kcs_nf*iLaKfu)*(K_uu\K_nu') + 2.*(Kcs_nf*L)*(L'*K_fu*(K_uu\K_nu'));
            end
        % Prediction with only CS covariance
        elseif ptype == 2
            Covft = Knn_v - sum((Kcs_nf/chol(La)).^2,2) + sum((Kcs_nf*L).^2, 2) ;

        end
    end

% $$$     Lav_pr = Kv_ff-Qv_ff;
% $$$     K2 = B'*B + K_cs + diag(Lav_pr);
% $$$     C = B'*B + K_cs + diag(Lav);
% $$$     K = B'*B + Kcs_nf + diag(Lav_pr);
% $$$
% $$$     L = chol(C)';
% $$$     %    y=K'*(C\y);
% $$$     a = L'\(L\y);
% $$$     Eft = K'*a;
% $$$
% $$$     v = L\K;
% $$$     V = gp_trvar(gp,xt,predcf);
% $$$     Covft = V - diag(v'*v);

    if nargout > 2
        Eyt = Eft;
        Covyt = Covft + diag(Cnn_v) - diag(Knn_v);
        if ~isempty(yt)
          ljpyt = mnorm_lpdf(yt', Eyt', Covyt);
        end
    end

  case {'VAR' 'DTC' 'SOR'}
    % Check the tstind vector
    if nargin > 5
        if ~isempty(tstind) && length(tstind) ~= size(tx,1)
            error('tstind (if provided) has to be of same lenght as tx.')
        end
    else
        tstind = [];
    end

    u = gp.X_u;
    m = size(u,1);
    % Turn the inducing vector on right direction
    if size(u,2) ~= size(x,2)
        u=u';
    end
    % Calculate some help matrices
    [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
    K_fu = gp_cov(gp, x, u);   % f x u
    K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
    K_nu = gp_cov(gp,xt,u);       % n x u
    Luu = chol(K_uu, 'lower');

    % Evaluate the Lambda (La) for specific model
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    B=Luu\(K_fu');
    Qv_ff=sum(B.^2)';
    Lav2 = diag(Cv_ff-Kv_ff);
    Lav = Cv_ff-Kv_ff;   % 1 x f, Vector of diagonal elements
                         % iLaKfu = diag(inv(Lav))*K_fu = inv(La)*K_fu
    iLaKfu = zeros(size(K_fu));  % f x u,
    n=size(x,1)
    for i=1:n
        iLaKfu(i,:) = K_fu(i,:)./Lav(i);  % f x u
    end
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;

    L = iLaKfu/chol(A);
    p = y./Lav - L*(L'*y);

    % Prediction matrices formed with only subset of cf's.
    if ~isempty(predcf)
        K_fu = gp_cov(gp, x, u, predcf);   % f x u
        K_uu = gp_trcov(gp, u, predcf);     % u x u, noiseles covariance K_uu
        K_nu = gp_cov(gp,xt,u,predcf);       % n x u
    end
    Eft = K_nu*(K_uu\(K_fu'*p));

    if nargout > 1
%         [Knn_v, Cnn_v] = gp_trvar(gp,xt,predcf);
        [Knn_v, Cnn_v] = gp_trcov(gp,xt,predcf);
        Luu = chol(K_uu,'lower');
        B=Luu\(K_fu');
        B2=Luu\(K_nu');

        Covftr = sum(B2'.*(B*bsxfun(@ldivide,Lav2,B')*B2)',2) - sum((K_nu*(K_uu\(K_fu'*L))).^2, 2);
%         Covftr =
        switch gp.type
          case {'VAR' 'DTC'}
            Covft = Knn_v - Covftr;         % Knn_v = diag(K_*,*)
          case  'SOR'
            Covft = sum(B2.^2,1)' - Covftr;     % sum(B2.^2,1)' = diag(Q_*,*)
        end

    end
    if nargout > 2
        Eyt = Eft;
        switch gp.type
          case {'VAR' 'DTC'}
            Covyt = Covft + Cnn_v - Knn_v;
          case 'SOR'
            Covyt = Covft + Cnn_v - sum(B2.^2,1)';
        end
    end
    if nargout > 4
        pyt = norm_pdf(y, Eyt, sqrt(Covyt));
    end

  case 'SSGP'
    if nargin > 4
        error(['Prediction with a subset of original ' ...
               'covariance functions not currently implemented with SSGP']);
    end

    [Phi_f, S] = gp_trcov(gp, x);
    Phi_a = gp_trcov(gp, xt);
    m = size(Phi_f,2);
    ns = eye(m,m)*S(1,1);

    L = chol(Phi_f'*Phi_f + ns,'lower');
    Eft = Phi_a*(L'\(L\(Phi_f'*y)));

    if nargout > 1
        Covft = sum(Phi_a/L',2)*S(1,1);
    end
    if nargout > 2
        error('gp_pred with three output arguments is not implemented for SSGP!')
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function gX = whiteKernGradX(kern, X, X2)

% WHITEKERNGRADX Gradient of WHITE kernel with respect to input locations.
% FORMAT
% DESC computes the gradident of the white noise
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO whiteKernParamInit, kernGradX, whiteKernDiagGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

if nargin<3
    X2 = X;
end

gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Anti-smearing function                     %
%       code by Fabrizio Conso, university of pavia, student       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                  %
 function fout = antismear(fin,fs,nfft)
 % fin  = normalized signal frequency
 % fs   = normalized sampling frequency
 % nfft = number of fft points
 format long;
 ts=1/fs;
 c=primes(nfft*ts*fin);
 maxprimes=max(c);
 fout=maxprimes/(nfft*ts);
%                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable

repo: mcbench-benchmarks


text:
% [LEV,IND] = wpyrLev(PYR,INDICES,LEVEL)
%
% Access a level from a separable QMF/wavelet pyramid.
% Return as an SxB matrix, B = number of bands, S = total size of a band.
% Also returns an Bx2 matrix containing dimensions of the subbands.

% Eero Simoncelli, 6/96.

function [lev,ind] =  wpyrLev(pyr,pind,level)

if ((pind(1,1) == 1) | (pind(1,2) ==1))
  nbands = 1;
else
  nbands = 3;
end

if ((level > wpyrHt(pind)) | (level < 1))
  error(sprintf('Level number must be in the range [1, %d].', wpyrHt(pind)));
end

firstband = 1 + nbands*(level-1)
firstind = 1;
for l=1:firstband-1
  firstind = firstind + prod(pind(l,:));
end

ind = pind(firstband:firstband+nbands-1,:);
lev  = pyr(firstind:firstind+sum(prod(ind'))-1);

repo: SelfExSR


text:
function Archive = UpdateArchive(Population,NA)
% Update the archive in MOEA/IGD-NS

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Detect the non-dominated solutions
    Population = Population(NDSort(Population.objs,1)==1);

    %% Select the extreme solutions
    Choose          = false(1,length(Population));
    [~,extreme]     = max(Population.objs,[],1);
    Choose(extreme) = true;

    %% Select other solutions by truncation
    if sum(Choose) > NA
        selected = find(Choose);
        Choose   = selected(randperm(length(selected),NA));
    else
        Cosine = 1 - pdist2(Population.objs,Population.objs,'cosine');
        Cosine(logical(eye(length(Cosine)))) = 0;
        while sum(Choose) < NA && ~all(Choose)
            unSelected = find(~Choose);
            [~,x]      = min(max(Cosine(~Choose,Choose),[],2));
            Choose(unSelected(x)) = true;
        end
    end
    Archive = Population(Choose);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


%% REW_city
% %{
imfolder = 'images\REW_city';
im_n = 5;
imfile = cell(im_n,1);
for ii = 1:im_n
    imfile{ii} = sprintf('%s\\%01d.jpg', imfolder, ii);
end

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = zeros(im_n-1,2);
for ei = 1:im_n-1
    edge_list(ei,:) = [ei,ei+1];
end

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 0, 'equi', 0, imfolder );
%}


text:
function sparse_grid_mixed_write_tests ( )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_WRITE_TESTS calls SPARSE_GRID_MIXED_WRITE with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    Local, real TOL, a tolerance for point equality.
%    A value of sqrt ( eps ) is reasonable, and will allow the code to
%    consolidate points which are equal, or very nearly so.  A value of
%    -1.0, on the other hand, will force the code to use every point, regardless
%    of duplication.
%
  addpath ( '../sandia_rules' );

  tol = sqrt ( eps );

  fprintf ( 1, '\n' );
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WRITE_TESTS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Call SPARSE_GRID_MIXED_WRITE_TEST with various arguments.\n' );
  fprintf ( 1, '  All tests will use a point equality tolerance of %e\n', tol );

  dim_num = 2;
  level_max = 2;
  rule = [ 1, 1 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l2_ccxcc';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 2;
  level_max = 3;
  rule = [ 1, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l3_ccxgp';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 2;
  level_max = 2;
  rule = [ 1, 4 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l2_ccxgl';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 2;
  level_max = 2;
  rule = [ 1, 7 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l2_ccxlg';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 2;
  level_max = 2;
  rule = [ 1, 8 ]';
  alpha = [ 0.0, 1.5 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l2_ccxglg';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 2;
  level_max = 2;
  rule = [ 2, 9 ]';
  alpha = [ 0.0, 0.5 ]';
  beta = [ 0.0, 1.5 ]';
  file_name = 'sparse_grid_mixed_d2_l2_f2xgj';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 2;
  level_max = 2;
  rule = [ 6, 4 ]';
  alpha = [ 2.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l2_gghxgl';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );

  dim_num = 3;
  level_max = 2;
  rule = [ 1, 2, 5 ]';
  alpha = [ 0.0, 0.0, 0.0 ]';
  beta = [ 0.0, 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d3_l2_ccxf2xgh';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );
%
%  Dimension 2, Level 4, Rule 3
%
  dim_num = 2;
  level_max = 4;
  rule = [ 3, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l4_gpxgp';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );
%
%  Dimension 2, Level 4, Rule 13
%
  dim_num = 2;
  level_max = 4;
  rule = [ 13, 13 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l4_gpsexgpse';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );
%
%  Dimension 2, Level 4, Rule 16
%
  dim_num = 2;
  level_max = 4;
  rule = [ 16, 16 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l4_gpmexgpme';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );
%
%  Dimension 2, Level 4, Rule 17
%
  dim_num = 2;
  level_max = 4;
  rule = [ 17, 17 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  file_name = 'sparse_grid_mixed_d2_l4_ccnxccn';
  sparse_grid_mixed_write_test ( dim_num, level_max, rule, alpha, beta, tol, ...
    file_name );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WRITE_TESTS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
% another simple 1D convection example with periodic BC
% diffusion coefficients are defined but not used
% It compares the results of upwind with TVD and shows how diffusive the upwind
% See how diffusive the upwind scheme can be althou it is so bad
% everywhere. Play with flux limiters, time steps, initial condition and time steps
% and see the difference between schemes in action
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
% define a 1D domain and mesh
W = 1;
Nx = 500;
mesh1 = createMesh1D(Nx, W);
x = mesh1.cellcenters.x;
% define the boundaries
BC = createBC(mesh1); % all Neumann
BC.left.periodic=1;
BC.right.periodic =1;
% Initial values
phi_old = createCellVariable(mesh1, 0.0, BC);
phi_old.value(20:120) = 1;
phi_old.value(180:400)= sin(x(180:400)*10*pi());
% initial guess for phi
phi = phi_old;
% initial values for upwind scheme
phiuw = phi;
phiuw_old=phi;
% keep the initial values for visualization
phiinit=phi_old;
% velocity field
u = 0.3;
uf = createFaceVariable(mesh1, u);
% diffusion field
D = 1e-2;
Df = createFaceVariable(mesh1, D);
% transient term coefficient
alfa = createCellVariable(mesh1,1.0);
% upwind convection term
Mconvuw = convectionUpwindTerm1D(uf);
% define the BC term
[Mbc, RHSbc] = boundaryCondition(BC);
% choose a flux limiter
FL = fluxLimiter('Superbee');
% solver
dt = 0.0005; % time step
final_t = W/u;
t = 0;
while t<final_t
    t = t+dt;
    % inner loop for TVD scheme
    [Mt, RHSt] = transientTerm(phi_old, dt, alfa);
    for j = 1:5
        [Mconv, RHSconv] = convectionTvdTerm1D(uf, phi, FL);
        M = Mconv+Mt+Mbc;
        RHS = RHSt+RHSbc+RHSconv;
        phi = solvePDE(mesh1, M, RHS);
    end
    [Mtuw, RHStuw] = transientTerm(phiuw_old, dt, alfa);
    Muw = Mconvuw+Mtuw+Mbc;
    RHSuw = RHStuw+RHSbc;
    phiuw = solvePDE(mesh1, Muw, RHSuw);
    phiuw_old = phiuw;
    phi_old = phi;
    figure(1);plot(x, phiinit.value(2:Nx+1), x, phi.value(2:Nx+1), '-o', x, ...
        phiuw.value(2:Nx+1)); drawnow;

end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function R2 = rsquare(y,yhat)
% PURPOSE:  calculate r square using data y and estimates yhat
% -------------------------------------------------------------------
% USAGE: R2 = rsquare(y,yhat)
% where:
%        y are the original values as vector or 2D matrix and
%        yhat are the estimates calculated from y using a regression, given in
%        the same form (vector or raster) as y
% -------------------------------------------------------------------------
% OUTPUTS:
%        R2 is the r square value calculated using 1-SS_E/SS_T
% -------------------------------------------------------------------
% Note: NaNs in either y or yhat are deleted from both sets.
%
% Felix Hebeler, Geography Dept., University Zurich, Feb 2007

if nargin ~= 2
    error('This function needs some exactly 2 input arguments!');
end

% reshape if 2d matrix
yhat=reshape(yhat,1,size(yhat,1)*size(yhat,2));
y=reshape(y,1,size(y,1)*size(y,2));

% delete NaNs
while sum(isnan(y))~=0 || sum(isnan(yhat))~=0
    if sum(isnan(y)) >= sum(isnan(yhat))
        yhat(isnan(y))=[];
        y(isnan(y))=[];
    else
        y(isnan(yhat))=[];
        yhat(isnan(yhat))=[];
    end
end

% 1 - SSe/SSt
R2 = 1 - ( sum( (y-yhat).^2 ) / sum( (y-mean(y)).^2 ) );

% SSr/SSt
% R2 = sum((yhat-mean(y)).^2) / sum( (y-mean(y)).^2 ) ;

if R2<0 || R2>1
    error(['R^2 of ',num2str(R2),' : yhat does not appear to be the estimate of y from a regression.'])
end


text:
function [xUpdate,PUpdate,innov,Pzz,W]=cubKalUpdate(xPred,PPred,z,R,h,xi,w,innovTrans,measAvgFun,stateDiffTrans,stateTrans)
%CUBKALUPDATE Perform the measurement update step in the cubature Kalman
%             filter with additive measurement noise. Unlike the square
%             root version, one can use cubature points with negative
%             weights.
%
%INPUTS: xPred The xDim X 1 predicted target state.
%        PPred The xDim X xDim predicted state covariance matrix.
%            z The zDim X 1 vector measurement.
%            R The zDim X zDim measurement covariance matrix in the native
%              coordinate system of the measurement.
%            h A function handle for the measurement function that takes
%              the state as its argument.
%           xi An xDim X numCubPoints matrix of cubature points. If this
%              and the next parameter are omitted or empty matrices are
%              passed, then fifthOrderCubPoints(xDim) is used. It is
%              suggested that xi and w be provided to avoid needless
%              recomputation of the cubature points.
%            w A numCubPoints X 1 vector of the weights associated with the
%              cubature points.
%   innovTrans An optional function handle that computes and optionally
%              transforms the value of the difference between the
%              observation and any predicted points. This is called as
%              innovTrans(a,b) and the default if omitted or an empty
%              matrix is passed is @(a,b)bsxfun(@minus,a,b). This must be
%              able to handle sets of values. For a zDimX1 measurement,
%              either of the inputs could be zDimXN in size while one of
%              the inputs could be zDimX1 in size.  This only needs to be
%              supplied when a measurement difference must be restricted
%              to a certain range. For example, the innovation between two
%              angles will be 2*pi if one angle is zero and the other
%              2*pi, even though they are the same direction. In such an
%              instance, a function handle to the
%              wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%              appropriate parameters should be passed for innovTrans.
%   measAvgFun An optional function handle that, when given N measurement
%              values with weights, produces the weighted average. This
%              function only has to be provided if the domain of the
%              measurement is not linear. For example, when averaging
%              angular values, then the function meanAng should be used.
% stateDiffTrans An optional function handle that takes an xDimXN matrix of
%              N differences between states and transforms them however
%              might be necessary. If not transformation is necessary, this
%              parameter can be omitted or an empty matrix passed.
%   stateTrans An optional function that takes a state estimate and
%              transforms it. This is useful if one wishes the elements of
%              the state to be bound to a certain domain. For example, if
%              an element of the state is an angle, one might generally
%              want to bind it to the region +/-pi.
%
%OUTPUTS: xUpdate The xDim X 1 updated state vector.
%         PUpdate The updated xDim X xDim state covariance matrix.
%      innov, Pzz The zDimX1 innovation and the zDimXzDim innovation
%                 covariance matrix are returned in case one wishes to
%                 analyze the consistency of the estimator or use those
%                 values in gating or likelihood evaluation.
%               W The xDimXzDim gain used in the update. This can be
%                 useful when gating and using the function
%                 calcMissedGateCov.
%
%If the function h needs additional parameters beyond the state, then the
%parameters can be passed by using an anonymous function as the function
%handle. For example, suppose that the measurement function is measFunc and
%it needs the additional parameters param1 and param2. In this instance,
%rather than using
%h=@measFunc
%one should use
%h=@(x)measFunc(x,param1,param2)
%This way, every time cubKalUpdate calls measFunc (via h) with a
%different x, those two parameters are always passed.
%
%The mathematics behind the function cubKalUpdate are described in
%detail in Section IX of [1] and in [2]. Note that this is essentially the
%"unscented Kalman filter" with additive noise. One simply has to provide
%the filter with the appropriate cubature points and weights.
%
%The optional parameters innovTrans and measAvgFun are not described in
%references [1] and [2], but allow for possible modifications to the filter
%as described in [3]. The parameters have been added to allow the filter to
%be used with angular quantities. For example, if the measurement consisted
%of range and angle, z=[r;theta], then
%innovTrans=@(a,b)[bsxfun(@minus,a(1,:),b(1,:));
%                  wrapRange(bsxfun(@minus,a(2,:),b(2,:)),-pi,pi)];
%measAvgFun=@(z,w)[calcMixtureMoments(z(1,:),w);
%                  meanAng(z(2,:),w')];
%should be used to approximately deal with the circular nature of the
%measurements.
%
%REFERENCES:
%[1] D. F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems Magazine,
%    vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[2] I. Arasaratnam and S. Haykin, "Cubature Kalman filters," IEEE
%    Transactions on Automatic Control, vol. 54, no. 6, pp. 1254-1269,
%    Jun. 2009.
%[3] D. F. Crouse, "Cubature/ unscented/ sigma point Kalman filtering with
%    angular measurement models," in Proceedings of the 18th International
%    Conference on Information Fusion, Washington, D.C., 6-9 Jul. 2015.
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    xDim=size(xPred,1);

    if(nargin<6||isempty(xi))
        [xi,w]=fifthOrderCubPoints(xDim);
    end

    if(nargin<8||isempty(innovTrans))
        %The function just returns the input difference.
        innovTrans=@(a,b)bsxfun(@minus,a,b);
    end

    if(nargin<9||isempty(measAvgFun))
        measAvgFun=@(zPoints,w)calcMixtureMoments(zPoints,w);
    end

    if(nargin<10||isempty(stateDiffTrans))
        stateDiffTrans=@(x)x;
    end

    if(nargin<11||isempty(stateTrans))
        stateTrans=@(x)x;
    end

    zDim=size(z,1);

    numCubPoints=size(xi,2);

    %cholSemiDef is used instead of chol in case a positive semi-definite
    %covariance matrix is passed.
    SPred=cholSemiDef(PPred,'lower');
    %Predicted cubature state points
    xPredPoints=stateTrans(transformCubPoints(xi,xPred,SPred));

    %Predicted cubature measurement points
    zPredPoints=zeros(zDim,numCubPoints);
    for curP=1:numCubPoints
        zPredPoints(:,curP)=h(xPredPoints(:,curP));
    end

    %Measurement prediction.
    zPred=measAvgFun(zPredPoints,w);

    %The innovation, transformed as necessary to keep values in a desired
    %range.
    innov=innovTrans(z,zPred);

    xPredCenPoints=stateDiffTrans(bsxfun(@minus,xPredPoints,xPred));
    %Centered, predicted cubature measurement points, transformed as
    %necessary to keep the values within a desired range.
    zPredCenPoints=innovTrans(zPredPoints,zPred);
    Pzz=R;
    Pxz=zeros(xDim,zDim);
    for curP=1:numCubPoints
        diff=zPredCenPoints(:,curP);
        Pzz=Pzz+w(curP)*(diff*diff');
        Pxz=Pxz+w(curP)*xPredCenPoints(:,curP)*diff';
    end

    %The filter gain
    W=Pxz/Pzz;

    %Updated state estimate
    xUpdate=stateTrans(xPred+W*innov);

    %Updated state covariance matrix
    %We could just do a simple one-line solution as in [1] and [2].
    %However, that does not guarantee that PUpdate will always be position
    %(semi)definite. Thus, we use an equivalent but more complicated update
    %formula based on Equation in Appendix C of [1].
    PUpdate=W*R*W';
    for curP=1:numCubPoints
        diff=stateDiffTrans(xPredCenPoints(:,curP)-W*zPredCenPoints(:,curP));
        PUpdate=PUpdate+w(curP)*(diff*diff');
    end
    %Ensure symmetry
    PUpdate=(PUpdate+PUpdate')/2;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef TestTransientAreasSegmentationModule
    %TestTransientAreasSegmentationModule

    properties (Constant)
        filename = fullfile(mexopencv.root(),'test','balloon.jpg');
    end

    methods (Static)
        function test_run
            img = cv.imread(TestTransientAreasSegmentationModule.filename, 'Color',true);
            sz = [size(img,2) size(img,1)];

            retina = cv.Retina(sz);
            retina.run(img);
            magno = retina.getMagnoRAW();

            seg = cv.TransientAreasSegmentationModule(sz);
            seg.clearAllBuffers();
            seg.run(magno);

            sz = seg.getSize();
            validateattributes(sz, {'numeric'}, {'vector', 'integer', 'numel',2});

            transientAreas = seg.getSegmentationPicture();
            validateattributes(transientAreas, {'uint8'}, {'size',[sz(2) sz(1)]});
            transientAreas = logical(transientAreas);
        end

        function test_params
            img = cv.imread(TestTransientAreasSegmentationModule.filename, 'Color',true);
            seg = cv.TransientAreasSegmentationModule([size(img,2) size(img,1)]);

            fname = [tempname() '.xml'];
            cObj = onCleanup(@() TestTransientAreasSegmentationModule.deleteFile(fname));
            seg.write(fname);
            seg.setup(fname);

            seg.setupParameters('LocalEnergyTemporalConstant',0.5);

            str = seg.printSetup();
            validateattributes(str, {'char'}, {'row', 'nonempty'});

            params = seg.getParameters();
            validateattributes(params, {'struct'}, {'scalar'});
        end
    end

    %% helper functions
    methods (Static)
        function deleteFile(fname)
            if exist(fname, 'file') == 2
                delete(fname);
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% Load example Data
clear all;
load('/home/guillaume/MatlabWorkSpace/Gaussian/MergeGMMs/Data/X1.mat');
load('/home/guillaume/MatlabWorkSpace/Gaussian/MergeGMMs/Data/X2.mat');
%% Plot data

close all;
figure; grid on; box on; hold on;
plot(X1(:,1),X1(:,2),'-b');
plot(X2(:,1),X2(:,2),'-r');
legend('X1','X2');
axis equal;
title('training data for GMMs');

%% GMM model for X1
options = statset('Display','final','MaxIter',1000);

GMModel = fitgmdist(X1,80,'Options',options,'CovarianceType','full','RegularizationValue',1e-03,'Replicates',5,'Start','plus');

gmm_x1.Priors = GMModel.ComponentProportion;
gmm_x1.Mu     = GMModel.mu';
gmm_x1.Sigma  = GMModel.Sigma;

GMModel       = fitgmdist(X2,80,'Options',options,'CovarianceType','full','RegularizationValue',1e-03,'Replicates',5,'Start','plus');
gmm_x2.Priors = GMModel.ComponentProportion;
gmm_x2.Mu     = GMModel.mu';
gmm_x2.Sigma  = GMModel.Sigma;

%% Plot Two orignial GMM

close all;
figure; grid on; box on; hold on;
plot_gmm_contour(gca,gmm_x1.Priors,gmm_x1.Mu,gmm_x1.Sigma,[0 0 1],3);
plot_gmm_contour(gca,gmm_x2.Priors,gmm_x2.Mu,gmm_x2.Sigma,[1 0 0],3);
title('GMM (1) & (2)');
axis equal;

%% Merge GMM (1) & (2)
options     = statset('Display','off','MaxIter',1000);
[gmm_x3] = merge_gmms(gmm_x1,X1,gmm_x2,X2,0.1,options);

%% Plot merged GMM

figure; grid on; box on; hold on;
plot_gmm_contour(gca,gmm_x3.Priors,gmm_x3.Mu,gmm_x3.Sigma,[0 0.8 0],3);
title('GMM (3)');
axis equal;

%% Learn GMM

options = statset('Display','final','MaxIter',1000);
K       = 100;
BIC     = zeros(K,1);
AIC     = zeros(K,1);

tic
parfor k = 1:K
    disp(['k(' num2str(k) ')']);
    GMModel = fitgmdist(X2,k,'Options',options,'CovarianceType','full','RegularizationValue',1e-03,'Replicates',1,'Start','plus');
    BIC(k)  = GMModel.BIC;
    AIC(k)  = GMModel.AIC;
end
toc

%%

figure; hold on;
plot(BIC,'-r');
plot(AIC,'-b');
legend('BIC','AIC');
xlabel('k');
title(' AIC & BIC vs K');

author: epfl-lasa