%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
% Tutorial for FAIR: regulariztion
%
% given a certain force field f on a 2D domain, this tutorial
% computes and visualizes the elastic displacement u, such that
%
%  B'*B u = f
%
% where B is the discrete elastic operator
% see also getElasticMatrixStg
%==============================================================================

clear, close all, help(mfilename);

                    % 2D example
omega  = [0,1,0,1]; % physical domain
m      = [16,12];   % number of discretization points
mu     = 1;         % Lame constants, control elasticity properties
lambda = 0;         % Youngs modulus and Poisson ratio
n      = prod(m);   % number of cells, and staggered dimensions
ns     = [(m(1)+1)*m(2),m(1)*(m(2)+1)];

% generate cell centered, staggered and nodal grids
% note: computation is staggered, others used for visualization
xc = getCellCenteredGrid(omega,m);
xn = getNodalGrid(omega,m);

% build elasticity operator on a staggered grid and visualize
B = getElasticMatrixStg(omega,m,mu,lambda);
FAIRfigure(1); clf; subplot(2,2,1); spy(B); title('B elastic on staggered grid')

% create a force field o cell centered grid and visualize
fc  = [0*xc(1:n);-10*sin(pi*xc(1:n))];
% computations on staggered grid
fs = grid2grid(fc,m,'centered','staggered');

% compute the displacement from (B'*B)*us = fs
% note B is derivative operator and has null space (constants)

% remove constant parts: E'*(f-E*w) = 0
E = [ones(ns(1),1)*[1,0];ones(ns(2),1)*[0,1]];
fs = fs - E*((E'*E)\(E'*fs));

% solve for displacements
warning off % matrix is singular, but MATLAB doesn't care
us = (B'*B)\fs;
warning on

% visualization on cell centered and nodal grids
uc = grid2grid(us,m,'staggered','centered');
un = grid2grid(us,m,'staggered','nodal');

% shortcuts for visualization,
vecNorm   = @(v) sqrt(sum(reshape(v,[],2).^2,2));
viewField = @(v) viewImage2Dsc(vecNorm(v),omega,m);
PG        = @(x) plotGrid(x,omega,m);

% visualize force field
FAIRfigure(1); subplot(2,2,3);
viewField(fc); hold on; colormap(gray); colorbar; PG(xn);
qh = quiver(xc(1:n),xc(n+1:end),fc(1:n),fc(n+1:end),1);
set(qh,'color','r','linewidth',1)
title('a force field')

% visualize displacement field
FAIRfigure(1); subplot(2,2,4);
viewField(uc); hold on; colormap(gray); colorbar; PG(xn);
qh = quiver(xc(1:n),xc(n+1:end),uc(1:n),uc(n+1:end),1);
set(qh,'color','g','linewidth',1)
title('the displacement field')

% visuaize the displaced grid
FAIRfigure(1); subplot(2,2,2);
xh = PG(xn); hold on; axis equal image; yh = PG(xn+un);
set(xh,'linewidth',2); set(xh,'linewidth',2,'color','g');
title('initial and displaced grids')

%==============================================================================


text:
function maximise( varargin )

%MAXIMISE Specifiies a concave (or affine) objective to be maximized.

cvx_problem = evalin( 'caller', 'cvx_problem', '[]' );
if isempty( cvx_problem ) || ~isa( cvx_problem, 'cvxprob' ),
    error( 'A cvx problem does not exist in this scope.' );
end
if ~isempty( cvx_problem.direction ),
    if isequal( cvx_problem.direction, 'find' ),
        error( 'Objective functions cannot be added to sets.' );
    else
        error( 'An objective function has already been supplied.' );
    end
end

if nargin < 1,
    error( 'Objective expression missing.' );
elseif iscellstr( varargin ),
    arg = evalin( 'caller', sprintf( '%s ', varargin{:} ) );
elseif nargin > 1,
    error( 'Too many input arguments.' );
else
    arg = varargin{1};
end

if ~isa( arg, 'cvx' ) && ~isa( arg, 'double' ) && ~isa( arg, 'sparse' ),
    error( 'Cannot accept an objective of type ''%s''.', class( arg ) );
end
persistent remap
if isempty( remap ),
    remap = cvx_remap( 'concave', 'log-concave' );
end
vx = remap( cvx_classify( arg ) );
if ~all( vx ),
    error( 'Disciplined convex programming error:\n   Cannot maximise a(n) %s expression.', cvx_class(arg(vx==0),false,true) );
end

newobj( cvx_problem, 'maximize', arg );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function pth = voxDataDir(view);
%
% pth = voxDataDir(view);
%
% Return the directory in which voxel-wise
% data is kept for the current view type/
% data type. Makes it if necessary.
%
% See: er_voxelData.m
%
% ras, 04/05.
pth = fullfile(dataDir(view),'VoxelData');

if ~exist(pth)
    % make it
    [p f] = fileparts(pth);
    try
        [success, msg] = mkdir(p,f);
    catch
        errmsg = sprintf('Couldn''t create voxel data dir: %s',msg)
        error(errmsg);
    end
    fprintf('Made directory %s for voxel data.\n',pth);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef PTKGetMaskForContext < PTKPlugin
    % PTKGetMaskForContext. Plugin for returning a mask corresponding to the
    %     current context
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Mask for context'
        ToolTip = 'Shows a mask corresponding to the current context'
        Category = 'Context'

        Context = PTKContextSet.Any
        AllowResultsToBeCached = false
        AlwaysRunPlugin = true
        PluginType = 'DoNothing'
        HidePluginInDisplay = true
        FlattenPreviewImage = false
        PTKVersion = '2'
        ButtonWidth = 6
        ButtonHeight = 1
        GeneratePreview = false
    end

    methods (Static)
        function results = RunPlugin(dataset, context, reporting)
            % Get a mask for the current region to analyse
            results = dataset.GetTemplateImage(context);
            if (context == PTKContext.LungROI) || (context == PTKContext.OriginalImage)
                results = dataset.GetTemplateImage(PTKContext.Lungs);
            end

            results.CropToFit;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout = fdsn_param_dialog(varargin)
    % FDSN_PARAM_DIALOG MATLAB code for fdsn_param_dialog.fig
    %      FDSN_PARAM_DIALOG, by itself, creates a new FDSN_PARAM_DIALOG or raises the existing
    %      singleton*.
    %
    %      H = FDSN_PARAM_DIALOG returns the handle to a new FDSN_PARAM_DIALOG or the handle to
    %      the existing singleton*.
    %
    %      FDSN_PARAM_DIALOG('CALLBACK',hObject,eventData,handles,...) calls the local
    %      function named CALLBACK in FDSN_PARAM_DIALOG.M with the given input arguments.
    %
    %      FDSN_PARAM_DIALOG('Property','Value',...) creates a new FDSN_PARAM_DIALOG or raises the
    %      existing singleton*.  Starting from the left, property value pairs are
    %      applied to the GUI before fdsn_param_dialog_OpeningFcn gets called.  An
    %      unrecognized property name or invalid value makes property application
    %      stop.  All inputs are passed to fdsn_param_dialog_OpeningFcn via varargin.
    %
    %      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
    %      instance to run (singleton)".
    %
    % See also: GUIDE, GUIDATA, GUIHANDLES

    % Edit the above text to modify the response to help fdsn_param_dialog

    % Last Modified by GUIDE v2.5 29-Jun-2018 16:25:23

    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
        'gui_Singleton',  gui_Singleton, ...
        'gui_OpeningFcn', @fdsn_param_dialog_OpeningFcn, ...
        'gui_OutputFcn',  @fdsn_param_dialog_OutputFcn, ...
        'gui_LayoutFcn',  [] , ...
        'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT

    % --- Executes just before fdsn_param_dialog is made visible.
function fdsn_param_dialog_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to fdsn_param_dialog (see VARARGIN)

    % Choose default command line output for fdsn_param_dialog
    handles.output = hObject;

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes fdsn_param_dialog wait for user response (see UIRESUME)
    % uiwait(handles.fdsn_import_dialog);

    % --- Outputs from this function are returned to the command line.
function varargout = fdsn_param_dialog_OutputFcn(hObject, eventdata, handles)

    % Get default command line output from handles structure
    varargout{1} = handles.output;

    % --- Executes on selection change in data_provider.
function data_provider_Callback(hObject, eventdata, handles)
    % Hints: contents = cellstr(get(hObject,'String')) returns data_provider contents as cell array
    %        contents{get(hObject,'Value')} returns selected item from data_provider

    currprovider = hObject.UserData(hObject.Value);
    set(handles.provider_details,'String',...
        sprintf('%s\n%s\n%s',...
        currprovider.description,...
        currprovider.location,...
        currprovider.serviceURLs.eventService));
    if hObject.Value==1
        hObject.BackgroundColor = [1.0 0.95 0.95];
        msg.infodisp('Importing FDSN data - First choose a data provider...','FDSN Fetch');

    else
        hObject.BackgroundColor = [0.95 1.0 0.95];
        msg.infodisp('Importing FDSN data - Choose the desired catalog constraints (time, magnitude, etc..)','FDSN Fetch');

    end
     if ~ handles.catalog_name.UserData.touched
         handles.catalog_name.String = [currprovider.name,'_fdsn']
     end

    % --- Executes during object creation, after setting all properties.
function data_provider_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    datacenter_details = webread('http://service.iris.edu/irisws/fedcatalog/1/datacenters');

    %dump datacenters with no event catalog access
    i=1;
    while i <= numel(datacenter_details)
        fldnm = fieldnames(datacenter_details(i).serviceURLs);
        if ~ismember('eventService',fldnm)
            datacenter_details(i)=[];
        else
            i=i+1;
        end
    end
    datacenter_details(2:end+1) = datacenter_details;
    datacenter_details(1).name = '<none>';
    datacenter_details(1).description = 'choose a datacenter';
    datacenter_details(1).website = ' ';
    datacenter_details(1).location = ' ';
    datacenter_details(1).lastUpdate = ' ';
    datacenter_details(1).serviceURLs.eventService = '';

    if exist('fdsnservices.json','file')
        try
            %TODO add a way to modify data centers
            % this is the datacenter_details structure, saved as a json file in the resrc directory
            jj=jsondecode(fileread('fdsnservices.json')); % get additional services
            for i=1:numel(jj)
                % only include datacenters that are not already retrieved by the querying fedcatalog
                if ~ismember(jj(i).name,{datacenter_details.name})
                    datacenter_details(end+1)=jj(i);
                end
            end
        catch ME
            disp(['unable to access additional datacenter information: ', ME.message]);
        end
    end
    hObject.UserData = datacenter_details;
    currprovider = datacenter_details(hObject.Value);
    set(hObject,'string',{datacenter_details.name});

    % --- Executes on button press in rect_geoselect.
function rect_geoselect_Callback(hObject, eventdata, handles)
    % Hint: get(hObject,'Value') returns toggle state of rect_geoselect
    set([handles.minlatitude, handles.maxlatitude, ...
        handles.minlongitude, handles.maxlongitude],'enable','on');

    set([handles.latitude,handles.longitude,...
        handles.minradius,handles.maxradius],'enable','off');

    enforce_gt_edit_relationship(handles.minlatitude, handles.maxlatitude);

    % --- Executes on button press in rad_geoselect.
function rad_geoselect_Callback(hObject, eventdata, handles)
    % Hint: get(hObject,'Value') returns toggle state of rad_geoselect
    set([handles.minlatitude, handles.maxlatitude, ...
        handles.minlongitude, handles.maxlongitude],'enable','off');
    set([handles.latitude,handles.longitude,...
        handles.minradius,handles.maxradius],'enable','on');
    enforce_gt_edit_relationship(handles.minradius, handles.maxradius);

function maxlatitude_Callback(hObject, eventdata, handles)

    update_editfield_value(hObject);
    enforce_gt_edit_relationship(handles.minlatitude, hObject);

    % --- Executes during object creation, after setting all properties.
function maxlatitude_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function minlatitude_Callback(hObject, eventdata, handles)
    update_editfield_value(hObject);
    enforce_gt_edit_relationship(hObject, handles.maxlatitude);

    % --- Executes during object creation, after setting all properties.
function minlatitude_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function minlongitude_Callback(hObject, eventdata, handles)
    update_editfield_value(hObject);
    if isnan(hObject.Value) || isempty(hObject.Value)
        set(hObject,'BackgroundColor',[1.0 1.0 1.0]);
    else
        set(hObject, 'BackgroundColor',[0.95 1.0 0.95]);
    end

    % --- Executes during object creation, after setting all properties.
function minlongitude_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function maxlongitude_Callback(hObject, eventdata, handles)

    update_editfield_value(hObject);
    if isnan(hObject.Value) || isempty(hObject.Value)
        set(hObject,'BackgroundColor',[1.0 1.0 1.0]);
    else
        set(hObject, 'BackgroundColor',[0.95 1.0 0.95]);
    end

    % --- Executes during object creation, after setting all properties.
function maxlongitude_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value=nan;

function starttime_Callback(hObject, eventdata, handles)

    try
        d = datetime(get(hObject,'String'));
        set(hObject,'backgroundcolor',[0.95 1.0 0.95]);
        set(hObject,'Value',datenum(d));
    catch
        d = str2double(get(hObject,'String'));
        if d >= 1800 && d < 3000
            %treat as year
            set(hObject,'String',[get(hObject,'String'), '-01-01']);
            starttime_Callback(hObject, eventdata, handles);
            return
        end

        set(hObject,'backgroundcolor',[1.0 0.95 0.95]);
        set(hObject,'Value',nan);
    end

    % --- Executes during object creation, after setting all properties.
function starttime_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function endtime_Callback(hObject, eventdata, handles)
    try
        d = datetime(get(hObject,'String'));
        set(hObject,'backgroundcolor',[0.95 1.0 0.95]);
        set(hObject,'Value',datenum(d));
    catch
        d = str2double(get(hObject,'String'));
        if d >= 1800 && d < 3000
            %treat as year
            set(hObject,'String',[get(hObject,'String'), '-01-01']);
            endtime_Callback(hObject, eventdata, handles);
            return
        end
        set(hObject,'backgroundcolor',[1.0 0.95 0.95]);
        set(hObject,'Value',nan);
    end

    % --- Executes during object creation, after setting all properties.
function endtime_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function mindepth_Callback(hObject, eventdata, handles)

    update_editfield_value(hObject);
    enforce_gt_edit_relationship(hObject, handles.maxdepth);

    % --- Executes during object creation, after setting all properties.
function mindepth_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function maxdepth_Callback(hObject, ~, handles)

    update_editfield_value(hObject);
    enforce_gt_edit_relationship(handles.mindepth, hObject);

    % --- Executes during object creation, after setting all properties.
function maxdepth_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function minmagnitude_Callback(hObject, ~, handles)

    update_editfield_value(hObject);
    enforce_gt_edit_relationship(hObject, handles.maxmagnitude)

    % --- Executes during object creation, after setting all properties.
function minmagnitude_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function maxmagnitude_Callback(hObject, ~, handles)
    update_editfield_value(hObject);
    enforce_gt_edit_relationship(handles.minmagnitude, hObject)

    % --- Executes during object creation, after setting all properties.
function maxmagnitude_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    hObject.Value = nan;

function magnitudetype_Callback(~, ~, ~)

    % --- Executes during object creation, after setting all properties.
function magnitudetype_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

function latitude_Callback(hObject, ~, ~)

    update_editfield_value(hObject);
    % --- Executes during object creation, after setting all properties.
function latitude_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

function longitude_Callback(hObject, ~, ~)

    update_editfield_value(hObject);
    % --- Executes during object creation, after setting all properties.
function longitude_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

function minradius_Callback(hObject, ~, handles)
    update_editfield_value(hObject);
    enforce_gt_edit_relationship(hObject, handles.maxradius);

    % --- Executes during object creation, after setting all properties.
function minradius_CreateFcn(hObject, ~, ~)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

function maxradius_Callback(hObject, ~, handles)

    update_editfield_value(hObject);
    enforce_gt_edit_relationship(handles.minradius, hObject);

    % --- Executes during object creation, after setting all properties.
function maxradius_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

    % --- Executes on button press in ignore_geo.
function ignore_geo_Callback(~, ~, handles)
    % hObject    handle to ignore_geo (see GCBO)

    % Hint: get(hObject,'Value') returns toggle state of ignore_geo
    set([handles.minlatitude, handles.maxlatitude, ...
        handles.minlongitude, handles.maxlongitude],'enable','off');
    set([handles.latitude,handles.longitude,...
        handles.minradius,handles.maxradius],'enable','off');

    % --- Executes on button press in Fetch.
function Fetch_Callback(~, ~, handles)
    % hObject    handle to Fetch (see GCBO)
    % assemble the actual query
    % TODO: do this
    ZG=ZmapGlobal.Data;
    if (handles.data_provider.Value == 1)
        % no datacenter has been chosen
        beep;
        h=errordlg('Incomplete Request: You must first choose a data provider.', ...
            ['Error:', get(get(gco,'Parent'),'Name')],'modal');
        waitfor(h);
        errorflash(handles.data_provider);
        return
    end
    queryset=handles.data_provider.String{handles.data_provider.Value};

    queryset = add_datetime(handles, queryset, 'starttime');
    queryset = add_datetime(handles, queryset, 'endtime');
    if get(handles.rect_geoselect,'Value')
        queryset = add_numeric(handles, queryset, 'minlatitude');
        queryset = add_numeric(handles, queryset, 'maxlatitude');
        queryset = add_numeric(handles, queryset, 'minlongitude');
        queryset = add_numeric(handles, queryset, 'maxlongitude');

    elseif get(handles.rad_geoselect,'Value')
        queryset = add_numeric(handles, queryset, 'latitude');
        queryset = add_numeric(handles, queryset, 'longitude');
        queryset = add_numeric(handles, queryset, 'minradius');
        queryset = add_numeric(handles, queryset, 'maxradius');
    end

    queryset = add_numeric(handles, queryset, 'minmagnitude');
    queryset = add_numeric(handles, queryset, 'maxmagnitude');
    queryset = add_numeric(handles, queryset, 'mindepth');
    queryset = add_numeric(handles, queryset, 'maxdepth');
    queryset = add_string(handles, queryset, 'magnitudetype');

    msg.infodisp('Importing FDSN data from the web. This might take a minute','FDSN Fetch');

    set(handles.fdsn_import_dialog,'Visible','off');
    drawnow;
    % FETCH THE DATA
    if ~iscell(queryset)
        h=errordlg('Incomplete Request: You must first choose some filter constraints', ...
            ['Error:', get(get(gco,'Parent'),'Name')],'modal');
        waitfor(h)
        set(handles.fdsn_import_dialog,'Visible','on');
        return
    end
    watchon;
    whos queryset

    m=msgbox_nobutton('Please wait while requested data is downloaded','Downloading');
    m.ButtonString='wait';
    %% do the import
    tmp=import_fdsn_event(1, queryset{:});

    %%
    % make it OK to close dialog box.
    if isvalid(m)
        m.Name='Done';
        m.String=sprintf('Done Downloading. Found %d events',tmp.Count);
        m.delay_for_close(1);
        delete(m);
    end
    watchoff;

    % CONVERT
    msg.infodisp('Converting to a ZmapCatalog','FDSN Fetch');
    if ~isa(tmp,'ZmapCatalog')
        ZG.primeCatalog=ZmapCatalog.from(tmp);
    else
        ZG.primeCatalog=tmp;
    end
    ZG.primeCatalog.sort('Date')
    assert(ZG.primeCatalog.Date(1)<=ZG.primeCatalog.Date(end));
    % provider_details=handles.data_provider.UserData(handles.data_provider.Value);
    if isempty(ZG.primeCatalog.Name) %TODO move this functionality into import_fdsn_event
        ZG.primeCatalog.Name = handles.catalog_name.String;
    end
    pause(.1)

    %name catalog
    %sdlg.prompt='Provide a catalog name (used in plots, files)';
    assert(ZG.primeCatalog.Date(1)<=ZG.primeCatalog.Date(end));

    if isvalid(m)
        close(m)
    end
    clear tmp

    uimemorize_catalog();

    % close(hObject.Parent); % or set Visibility to 'off' ?

function [val] = getvalue(handles, label)
    val = get(handles.(label),'Value');

function [optlist] = add_numeric(handles, optlist, label)
    % optlist = add_numeric(handles, optlist, 'something')
    % adds only if 'something is not nan or empty
    val = get(handles.(label),'Value');
    if ~isempty(val) && ~isnan(val)
        optlist = [optlist, {label, val}];
    end

function [optlist] = add_string(handles, optlist, label)
    % optlist = add_numeric(handles, optlist, 'something')
    % adds only if 'something is not nan or empty
    val = get(handles.(label),'String');
    if ~isempty(val) && ~any(val==' ')
        optlist = [optlist, {label, val}];
    end

function [optlist] = add_datetime(handles, optlist, label)
    % optlist = add_numeric(handles, optlist, 'something')
    % adds only if 'something is not nan or empty
    val = get(handles.(label),'Value');
    if ~isempty(val) && ~isnan(val)
        val=datestr(val,'yyyy-mm-dd hh:MM:SS');
        val(val==' ') = 'T';
        optlist = [optlist, {label, val}];
    end

function enforce_gt_edit_relationship(hSmaller, hBigger)
    nan_bgcolor = [1.0 1.0 1.0];
    good_bgcolor = [0.95 1.0 0.95];
    bad_bgcolor = [1.0 0.95 0.95];
    hasSmaller = ~isnan(hSmaller.Value) && ~isempty(hSmaller.Value);
    hasBigger = ~isnan(hBigger.Value) && ~isempty(hBigger.Value);
    if hasBigger && ~hasSmaller
        hSmaller.BackgroundColor = nan_bgcolor;
        hBigger.BackgroundColor = good_bgcolor;
        %
    elseif hasSmaller && ~hasBigger
        hSmaller.BackgroundColor = good_bgcolor;
        hBigger.BackgroundColor = nan_bgcolor;
    elseif ~hasBigger && ~hasSmaller
        hSmaller.BackgroundColor = nan_bgcolor;
        hBigger.BackgroundColor = nan_bgcolor;
    else
        if hBigger.Value >= hSmaller.Value
            hSmaller.BackgroundColor = good_bgcolor;
            hBigger.BackgroundColor = good_bgcolor;
        else
            hSmaller.BackgroundColor = bad_bgcolor;
            hBigger.BackgroundColor = bad_bgcolor;
        end
    end

function catalog_name_Callback(hObject, ~, ~)
% Hints: get(hObject,'String') returns contents of catalog_name as text
%        str2double(get(hObject,'String')) returns contents of catalog_name as a double
hObject.UserData.touched = true;

% --- Executes during object creation, after setting all properties.
function catalog_name_CreateFcn(hObject, ~, ~)
hObject.UserData.touched = false;

% --- Executes on button press in cancelbutton.
function cancelbutton_Callback(hObject, eventdata, handles)
delete(handles.fdsn_import_dialog);

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over cancelbutton.
function cancelbutton_ButtonDownFcn(hObject, eventdata, handles)

author: CelsoReyes


text:
function [ data ] = vis_wcf( label_list, obb_list )

linewidth = 2;
corlist = {[1,0.5,0.5,0.3],[0.5,0.7843,0.5,0.3],[0.5,0.5,1,0.3],[1,0.5,0,0.3],[0.1961,0.7843,1,0.3]};

front = [0,0,1];
up = [0,1,0];
axes = [1,0,0];

obblist = [obb_list(1:6,:);zeros(3,size(obb_list,2));obb_list(7:9,:)];
obblist(8,:) = 1;

for i = 1:length(label_list)
    if(i<3)
        % skip the visualization of floor and ceiling
        continue
    end
    p = obblist(:,i);

    center = p(1:3);
    lengths = p(10:12);

    dir_1 = p(4:6);
    dir_2 = p(7:9);
    dir_1 = dir_1/norm(dir_1);
    dir_2 = dir_2/norm(dir_2);
    dir_3 = cross(dir_1,dir_2);
    dir_3 = dir_3/norm(dir_3);
    d1 = 0.5*lengths(1)*dir_1;
    d2 = 0.5*lengths(2)*dir_2;
    d3 = 0.5*lengths(3)*dir_3;

    cornerpoints(1,:) = center-d1-d2-d3;
    cornerpoints(2,:) = center+d1-d2-d3;
    cornerpoints(3,:) = center-d1-d2+d3;
    cornerpoints(4,:) = center+d1-d2+d3;
    cornerpoints(:,1) = 1 - cornerpoints(:,1);

    cor = corlist{mod(i-1,length(corlist))+1};
    plot([cornerpoints(1,1),cornerpoints(2,1)],[cornerpoints(1,3),cornerpoints(2,3)],'Color',cor(1:3),'linewidth',linewidth);
    hold on
    plot([cornerpoints(2,1),cornerpoints(4,1)],[cornerpoints(2,3),cornerpoints(4,3)],'Color',cor(1:3),'linewidth',linewidth);
    plot([cornerpoints(4,1),cornerpoints(3,1)],[cornerpoints(4,3),cornerpoints(3,3)],'Color',cor(1:3),'linewidth',linewidth);
    plot([cornerpoints(3,1),cornerpoints(1,1)],[cornerpoints(3,3),cornerpoints(1,3)],'Color',cor(1:3),'linewidth',linewidth);

	label = label_list{i};
	min_x = min(cornerpoints(:,1));
	max_z = max(cornerpoints(:,3));
	t = text(min_x+0.01,max_z-0.02,[label]);
	t.BackgroundColor = cor;
	t.FontSize = 8;
end

axis equal
axis off
fig = gcf;
fig.PaperUnits = 'points';
fig.PaperPosition = [0 0 224 224];%168 %108

f=getframe(gcf);
data=f.cdata;

end

repo: GRAINS


text:
function turn_imagesc()
%TURN_IMAGES corrects the Y axis for imagesc plots

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Main ============================================================
set(gca,'YDir','normal');

author: sfstoolbox

repo: sfs-matlab


text:
function dz = dynAcc(z, u, p)
% dz = dynAcc(z, u, p)
%
% This function computes the dynamics of a simple planar quad-rotor
% helicopter, including chain integrator for acceleration cost function.
%
% INPUTS:
%   z = [9, n] = [X; V1; V2] = state matrix
%   u = [5, n] = [U1; U2] = control matrix
%   p = parameter struct:
%       .g = acceleration due to gravity
%       .d = half distance between rotors
%       .m = mass of each rotor (half-mass of the quad-rotor)
%
% OUTPUTS:
%   dz = derivative of state matrix
%

% Unpack the inputs
X = z(1:3,:);   % configuration
V1 = z(4:6,:);  % rates (dynamics)
% V2 = z(7:9,:);  % rates (integrator)    %unused
U1 = u(1:2,:);  % actuators
U2 = u(3:5,:);  % kinematics

% Call to the actual dynamics function
dV1 = dynQuadRotor(X, U1, p);

% Chain integrators
dX = V1;
dV2 = U2;

% Pack up the outputs
dz = [dX; dV1; dV2];

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [indpeak twfPer indtrough] = peakfinder(twf,fmin,fmax)
% [indpeak twfPer indtrough] = peakfinder(twf,<fmin>,<fmax>)
%
% Finds local peaks in a near-periodic waveform. This fails if the
% peaks are closer than 1/2 period. twfPer is the period of twf in
% samples. If the 1st or last peak are less than 90% of the mean of
% the rest of the peaks, they are exluded.
%
% fmin,fmax are constraints on the fundamental frequency and are
% given in units of items per time point (NOT IN Hz!).
%
% $Id: peakfinder.m,v 1.9 2010/05/18 19:31:44 greve Exp $

indpeak = [];

if(nargin < 1 | nargin > 3)
  fprintf('[indpeak twfPer indtrough] = peakfinder(twf,<fmin>,<fmax>)\n');
  return;
end

Ntp = length(twf);
nn = 1:Ntp;

% detrend - necessary?
X = fast_polytrendmtx(1,Ntp,1,3);
%twf = twf - X*(inv(X'*X)*(X'*twf));

% Get major period of waveform
[fftaxis, deltafreq] = fast_fftaxis(Ntp,1);
nfft = length(fftaxis);
nnfft = 1:nfft;
twffft = abs(fft(twf-mean(twf)));
if(exist('fmin','var'))
  indok = find(fftaxis >= fmin & fftaxis <= fmax);
else
  indok = [1:length(fftaxis)];
end
[tmp k] = max(twffft(indok));
twfFreq = fftaxis(indok(k));
twfPer = 1/twfFreq;
twfPerSamp      = round(twfPer);
twfHalfPerSamp  = round(twfPerSamp/2);
twfQuartPerSamp = round(twfPerSamp/4);

%plot(fftaxis,twffft(nnfft))
%keyboard

% Assume global peak is a local peak
[tmp k0] = max(twf);
%fprintf('global peak at %d %f\n',k0,k0/25);

% Look ahead, starting at global peak
indpeak = k0;
kprev = k0;
while(1)
  % Find next max by searhing over time starting at 1/2
  % period beyond the previous max and ending one period
  % later. This fails if the peaks are closer than 1/2
  % period.
  kstart = kprev + twfHalfPerSamp;
  k = kstart + [0:twfPerSamp-1];
  indok = find(k < Ntp);
  if(length(indok) < twfQuartPerSamp) break; end
  k = k(indok);
  [tmp mmax] = max(twf(k));
  kmax = kstart + mmax - 1;
  indpeak = [indpeak kmax];
  kprev = kmax;
end

% Look behind (reverse and look ahead)
twfrev = flipud(twf(:));
k0rev = Ntp - k0 + 1; % DONT = max(twfrev);
indpeakrev = []; % dont include k0rev here
kprev = k0rev;
while(1)
  kstart = kprev + twfHalfPerSamp;
  k = kstart + [0:twfPerSamp-1];
  indok = find(k < Ntp);
  if(length(indok) < twfQuartPerSamp) break; end
  k = k(indok);
  [tmp mmax] = max(twfrev(k));
  kmax = kstart + mmax - 1;
  indpeakrev = [indpeakrev kmax];
  kprev = kmax;
end

% Convert reversed indices to forard indices
indpeakrevfor = Ntp - indpeakrev + 1;
indpeak = sort([indpeak indpeakrevfor]);
npeaks = length(indpeak);

% Decide whether to eliminate the first peak
% Compute mean of closest 3 peaks
peakfirst = twf(indpeak(1));
indpm = [2:min(4,npeaks)];
peakmean = mean(twf(indpeak(indpm)));
% Must be greater than 0.7 times this mean
if(peakfirst < .7*peakmean)  indpeak = indpeak(2:end); end

% Decide whether to eliminate the last peak
% Compute mean of closest 3 peaks
peaklast  = twf(indpeak(end));
indpm = [max(npeaks-3,1),max(npeaks-1,1)];
peakmean = mean(twf(indpeak(indpm)));
% Must be greater than 0.7 times this mean
if(peaklast  < .7*peakmean)  indpeak = indpeak(1:end-1); end

% Make sure they are unique (why not always?)
indpeak = unique(indpeak);

npeaks = length(indpeak);
indtrough = zeros(size(indpeak));
for nthpeak = 1:npeaks
  i1 = indpeak(nthpeak);
  if(nthpeak < npeaks)
    i2 = indpeak(nthpeak+1);
  else
    i2 = Ntp;
  end
  if(i2 > i1 + 1) i1 = i1 + 1; end
  [mmin imin] = min(twf(i1:i2));
  indtrough(nthpeak) = imin+i1-1;
end

if(0)
nn = 1:Ntp;
plot(nn,twf,nn(indpeak),twf(indpeak),'*',nn(indtrough),twf(indtrough),'o');
fprintf('Period %f\n',twfPerSamp);
keyboard
end

return;


text:
function Hd = bandpass_filter_design(Fs)

%
% MATLAB Code
% Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
%
% Generated on: 23-Apr-2013 18:34:23
%

% Fstop1 = 0.001;  % First Stopband Frequency
% Fpass1 = 5;      % First Passband Frequency
% Fpass2 = 35;     % Second Passband Frequency
% Fstop2 = 50;     % Second Stopband Frequency
% Astop1 = 20;     % First Stopband Attenuation (dB)
% Apass  = 1;      % Passband Ripple (dB)
% Astop2 = 20;     % Second Stopband Attenuation (dB)
% % Fs     = 1000;   % Sampling Frequency
%
% h = fdesign.bandpass('fst1,fp1,fp2,fst2,ast1,ap,ast2', Fstop1, Fpass1, ...
%     Fpass2, Fstop2, Astop1, Apass, Astop2, Fs);
%
% Hd = design(h, 'butter', ...
%     'MatchExactly', 'stopband');

N     = 100;  % Order
Fpass = 35;   % Passband Frequency
Fstop = 50;   % Stopband Frequency
Wpass = 1;    % Passband Weight
Wstop = 100;  % Stopband Weight

% Calculate the coefficients using the FIRLS function.
b  = firls(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop]);
Hd1 = dfilt.dffir(b);

Fstop = 0.01;        % Stopband Frequency
Fpass = 1;           % Passband Frequency
Astop = 120;         % Stopband Attenuation (dB)
Apass = 1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
Hd2 = design(h, 'butter', 'MatchExactly', match);

Hd = dfilt.cascade(Hd2,Hd2);

% [EOF]

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [AuuCondensed,AuxCondensed,AxxCondensed] = OCP_Auu_Aux_Axx_Condensed(lambda,mu,u,x,z,p) %#codegen
    [AuuCondensed,AuxCondensed,AxxCondensed] = OCP_GEN_Auu_Aux_Axx_Condensed(lambda,mu,u,x,z,p);
end

author: deng-haoyang

repo: ParNMPC


text:
%rf_combineRF.m
%Peter Truong, Sunnybrook Research Institute 2021.
%
% USAGE:
% [rf,AMPINT]=rf_combineRF(rf_struct1,rf_struct2);
%
% DESCRIPTION:
% Combines two similar rf waveforms (e.g. same pulse, but one frequency
% shifted), in FID-A rf pulse structure format.
%
% INPUTS:
% rf_struct1 = RF pulse definition structure.
% rf_struct2 = RF pulse definition structure, similar to rf_struct1

%
% OUTPUTS:
% rf         = Output rf waveform of the combined rf pulse, in FID-A rf
%              pulse structure format.
% AMPINT     = Calculated amplitude integral (for use in Siemens .pta files).
function [rf,AMPINT]=rf_combineRF(rf_struct1,rf_struct2)

%Convert waveforms to complex form
rf1_waveform=rf_struct1.waveform(:,2).*exp(1i.*rf_struct1.waveform(:,1)*pi/180);
rf2_waveform=rf_struct2.waveform(:,2).*exp(1i.*rf_struct2.waveform(:,1)*pi/180);

%Here we compute the Amplitude integral(AMPINT), which is used by magnetom
%to calculate the transmitter power that is required in order to achieve
%the desired flip angle
rf1_waveform_scaled=rf1_waveform./max(abs(rf1_waveform));
AI=sum(rf1_waveform_scaled);

combined_waveform=rf1_waveform + rf2_waveform;
combined_waveform_scaled=combined_waveform./max(abs(combined_waveform));
AMPINT=AI./max(abs(combined_waveform_scaled));

rf=rf_struct1;
rf.waveform(:,1)=phase(combined_waveform_scaled).*180/pi;
rf.waveform(:,2)=abs(combined_waveform_scaled);

% updating time-bandwidth product & time-w1max product
rf.tbw=rf_struct1.tbw+rf_struct2.tbw;
rf.tw1=rf_struct1.tw1+rf_struct2.tw1;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%RECALL   query the current discretized problem in the workspace
%IFISS scriptfile: DJS; 27 April 2012.
%Copyright (c) 2009 D.J. Silvester, H.C. Elman, A. Ramage
%
%For more help, type 'helpme'.
if exist('pde','var')==0,
error('No discrete problem found in the workspace!'),
end
%
if pde>10 & pde <20, fprintf(' Unsteady solution\n'), pdex=pde-10;
else, pdex=pde; end
%
if pdex<3,
fprintf(' scalar problem       : domain       : approximation')
     if pdex==1, fprintf('\n Poisson Equation     :'),
 elseif pdex==2, fprintf('\n Convection-Diffusion :'),
    end
    if domain==1, fprintf(' square       :'),
elseif domain==2, fprintf(' L-shaped     :'),
            else, fprintf(' non-standard :'),
    end
    if qmethod==1, fprintf(' Q1\n'),
elseif qmethod==2, fprintf(' Q2\n'),
			else,  fprintf(' non-standard\n'),
    end
else
fprintf(' flow problem  :        geometry          : approximation')
    if pdex==3, fprintf('\n Stokes        :'),
elseif pdex==4, fprintf('\n Navier-Stokes :'),
elseif pdex==5, fprintf('\n Boussinesq    :'),
    end
	if domain==1,  fprintf(' enclosed square          :'),
elseif domain==3,  fprintf(' step with inflow/outflow :'),
elseif domain==4,  fprintf(' channel with obstacle    :'),
elseif domain==7,  fprintf(' enclosed cavity          :'),
elseif domain==10, fprintf(' simple channel           :'),
	         else, fprintf(' non-standard'),
	end
	if qmethod==0, fprintf(' Q1-Q1\n'),
elseif qmethod==1, fprintf(' Q1-P0\n'),
elseif qmethod==2, fprintf(' Q2-Q1\n'),
elseif qmethod==3, fprintf(' Q2-P1\n'),
elseif qmethod==12, fprintf(' Q2-Q1-Q2\n'),
else, fprintf('non-standard\n'),
    end
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = PlotComet_3D(x_data, y_data, z_data, varargin)

%function creates a moving comet plot in 3 dimensions using the specificed
%x, y, and z line data
%
%PlotComet_3D(x_data, y_data, z_data)
%
%Optional Input Arguements:
%cFigure, handle to the current figure or axes to create the plot in (will
%   add the line data to the plot regardless of whether hold is on
%
%Frequency, playback frequency, if not provided assumed to be 10 Hz
%
%blockSize, number of points in the comet tail, assumed to be 1/20 of
%   length(x_data) if not provided
%
%PlotComet_3D(x_data, y_data, z_data, 'cFigure',FigureHandle,...
%     'Frequency',freq, 'blockSize',num_points)
%
%tailFormat, a structured variable containing the fomrating requirements of
%   the tail (if different than the default red solid line)
%   ex: Tail = struct('LineWidth',2,'Color','g','LineStyle','*');
%   PlotComet_3D(x_data,y_data,z_data,'tailFormat',Tail);
%
%headFormat, a structued vairable containing the formating requirements of
%   the head (if different than the default blue circle)
%   ex: Head = struct('LineStyle','none','MarkerSize',8,'Marker','^','Color','k');
%   PlotComet_3D(x_data,y_data,z_data,'headFormat',head);
%   Note: 'LineSytle','none' is required since the head is plotted as a
%   single point!
%
%Optional Output Arguements
%hFigure = figure handle to the comet plot
%
%Example: Create a surface plot and construct a moving line plot through it
%with connected cyan stars and dashed line on a 50 points long tail
%and large green square for a head.
%
%%create some data for the surface
%[X, Y] = meshgrid(linspace(-1, 1, 25), linspace(-1, 1, 25));
%Z = X .* exp(-X.^2 - Y.^2);
%%open a new figure and plot the surface
%surf(X, Y, Z);
%fig = gcf;
%%create the line data
%t = -pi:pi/500:pi;
%x = sin(5*t);
%y = cos(3*t);
%z = t;
%%define the tail formated structure
%Tail = struct('LineStyle','--','Marker','*','Color','c',...
%   'LineWidth',2,'MarkerSize',4);
%%define the head formated structure
%Head = struct('LineStyle','none','Marker','s','Color','g','MarkerSize',10);
%
%%Invoke the comet plot with 50 Hz playback, and a tail length of 50 pts
%PlotComet_3D(x,y,z,'cFigure',fig,'blockSize',50,'Frequency',50,...
%   'headFormat',Head,'tailFormat',Tail);
%
%Nick Hansford
%09/26/2013

%initialize the inputs
freq = 10;
blockSize = floor(1/20*length(x_data));
tailFormat = struct('LineWidth',1,'Color','r','LineStyle','-');
headFormat = struct('LineStyle','none','Marker','o','MarkerSize',6,...
    'Color','b');

%parse out the inputs
argCount = nargin - 3;

for i = 1:2:argCount
    %         caseVar = varargin{i}
    switch varargin{i}
        case 'cFigure'
            cFigure = varargin{i+1};
            %get the original size:
            cAxes = get(cFigure,'CurrentAxes');
            xLim = get(cAxes,'XLim');
            yLim = get(cAxes,'YLim');
            zLim = get(cAxes,'ZLim');

            %resize if the new plot will exceed them
            if xLim(1) > min(x_data)
                xLim(1) = min(x_data);
            end

            if xLim(2) < max(x_data)
                xLim(2) = max(x_data);
            end

            if yLim(1) > min(y_data)
                yLim(1) = min(y_data);
            end

            if yLim(2) < max(y_data)
                yLim(2) = max(y_data);
            end

            if zLim(1) > min(z_data)
                zLim(1) = min(z_data);
            end

            if zLim(2) < max(z_data)
                zLim(2) = max(z_data);
            end

            axis([xLim, yLim, zLim]);

        case 'blockSize'
            blockSize = varargin{i+1};
        case 'Frequency'
            freq = varargin{i+1};

        case 'tailFormat'
            tailFormat = varargin{i+1};

        case 'headFormat'
            headFormat = varargin{i+1};
    end
end

%make sure the figure exists, if not, create it
if ~exist('cFigure')
    cFigure = figure;
    view(3);
    axis([min(x_data), max(x_data),...
        min(y_data), max(y_data),...
        min(z_data), max(z_data)]);
end

%user can pass in current axes, if so, get the parent for the figure window
%and use that instead...should make this compatible with GUIs?
if strcmp(get(cFigure,'Type'),'axes')
    cFigure = get(cFigure,'Parent');
end

%activate the figure window
figure(cFigure);
cAxes = get(cFigure,'CurrentAxes');

oldNextPlot = get(cAxes,'NextPlot');
set(cAxes,'NextPlot','add');

pauseTime = 1./freq;

n_start = 1;
n_stop = 1;
%put on the starting point
plot3(x_data(n_start:n_stop),...
    y_data(n_start:n_stop),...
    z_data(n_start:n_stop),tailFormat);
plot3(x_data(n_stop), y_data(n_stop), z_data(n_stop),headFormat);

%playback!
for n = 1:1:length(x_data)
    a = tic;
    %delete the previous plot
    hChild = get(cAxes,'Children');
    delete(hChild(1:2));

    if n <= blockSize
        n_start = 1;
        n_stop = n;
    else
        n_start = n_start + 1;
        n_stop = n_stop + 1;
    end

    %new plot
    plot3(x_data(n_start:n_stop),...
        y_data(n_start:n_stop),...
        z_data(n_start:n_stop),tailFormat);
    plot3(x_data(n_stop), y_data(n_stop), z_data(n_stop),headFormat);
    drawnow;

    %update playback refresh rate
    b = toc(a);
    pause(pauseTime-b);
end

set(cAxes,'NextPlot',oldNextPlot);
% fprintf('Finished\n');

if nargout == 1
    varargout{1} = cFigure;
end

author: Sable

repo: mcbench-benchmarks


text:
function mfig = montagefigurefiles(x,y)
% by Nick Sinclair (nws5@pitt.edu)  Jul 08, 2010

% figurearray is row-vector of handles to figures that you'd like
% included in the montage

% x is horizontal number of subplots
% y is vertical number of subplots

% props array specifies properties that you'd like copied over to the new
% subplots. I did this because copying all of the properties gave me errors
% because I didn't know how to exclude the read-only properties. So for now
% you have to specify each one. Check 'Axes Properties' in help for more
% properties to copy over.

% The space between subplots can be adjusted with 'spacerx' and 'spacery'

% The axis labels and Ticks are only displayed on the first column and
% last row. If it's not appropriate, it's easy to change. For example, it's
% particularly ugly if you have an odd number of plots or varying xlim/ylim
% Just set LabelOnlyPerimeter = 0;

% returns handle to new figure mfig

% IMPORTANT ADJUSTABLE PARAMETERS
props ={'XLim','YLim','CLim','Xscale','Yscale', ...
    'YTickLabelMode','XTickLabelMode','View','YDir', ...
    'XDir','DataAspectRatio','DataAspectRatioMode'};
legprops = {'String','FontAngle','FontName','FontSize','FontUnits','FontWeight'};
    %properties to copy from legends

left = 0.03;  % margins (add it in if you like, in the same way as left)
top = .98;    %
bottom = 0.05;%
spacerx = 1/10; % size of space between subplots compared to subplot width
spacery = 1/5;  %
newfigposition = [.1 .1 .75 .75]; %set size and position of newfig
LabelOnlyPerimeter = 1;
% END IMPORTANT ADJUSTABLE PARAMETERS

[FileName,PathName] = uigetfile('*.*','Select Figure Files to Montage', 'MultiSelect','on');

for jj = 1:length(FileName),
    figurearray(jj) = openfig(strcat(PathName,char(FileName(jj))),'new','invisible');
end

mfig = figure('units', 'normalized', 'outerposition', newfigposition);

j = 1;

count = 0;

for ii = 1:length(figurearray),
   fig(ii).axhandles = findobj(get(figurearray(ii),'Children'),'Type','Axes','-not','Tag','Colorbar','-not','Tag','legend');
   count = count+length(fig(ii).axhandles); %count up axes excluding colorbars, incl subplots
end

% disp(['found: ',num2str(count),' plots']);

wxfactor = 1/spacerx;                   % convenience
wyfactor = 1/spacery;
dx = (1-left)/(((wxfactor+1)*x)+1);     %calc so that it all adds up to 1
dy = (top-bottom)/(((wyfactor+1)*y)+1); %
sizex = wxfactor*dx;                    %
sizey = wyfactor*dy;                    %

for i = 1:length(figurearray),
    for subaxis = 1:length(fig(i).axhandles)
        row = floor((j-1)/x) + 1;
        col = mod((j-1),x)+1;
        set(0,'CurrentFigure',mfig);
        pos(1) = left + col*(dx) + (col-1)*sizex; %set subplot position
        pos(3) = sizex;                           %
        pos(2) = top - row*dy - (row)*sizey;      %
        pos(4) = sizey;                           %
        s(j) = subplot('Position',pos); %handle to axis j in subplot (target)
%         set(0,'CurrentFigure',i);
        origaxis = fig(i).axhandles(subaxis);

        %This label complication is because the xlabel property of the axis is
        %only a handle to a child string of the axis. hence
        %get(axishandle,'xlabel') gives you the handle to the label which you
        %can then set properties for.

        copyobj(get(origaxis,'Children'),s(j)) %send all children to subplot j
        set(get(s(j),'Title'),'String',get(get(origaxis,'Title'),'String'));
        set(get(s(j),'Title'),'FontSize',(0.7*get(get(origaxis,'Title'),'FontSize')));
        set(s(j),props,get(origaxis,props));    %set the specified properties of axis j
        if LabelOnlyPerimeter,
            if (col==1),            % if first column set labels/ticks
                set(get(s(j),'ylabel'),'String',get(get(origaxis,'ylabel'),'String'));
                set(get(s(j),'ylabel'),'FontSize',(0.7*get(get(origaxis,'ylabel'),'FontSize')));
            else
                set(s(j),'YTickLabelMode','manual');
                set(s(j),'YTick',[]);
            end

            if (row == ceil(length(figurearray)/x)), %if last row, labels/ticks
                set(get(s(j),'xlabel'),'String',get(get(origaxis,'xlabel'),'String'));
                set(get(s(j),'xlabel'),'FontSize',(0.7*get(get(origaxis,'xlabel'),'FontSize')));
            else
                set(s(j),'XTickLabelMode','manual');
                set(s(j),'XTick',[]);
            end
        else
                set(get(s(j),'ylabel'),'String',get(get(origaxis,'ylabel'),'String'));
                set(get(s(j),'ylabel'),'FontSize',(0.7*get(get(origaxis,'ylabel'),'FontSize')));
                set(get(s(j),'xlabel'),'String',get(get(origaxis,'xlabel'),'String'));
                set(get(s(j),'xlabel'),'FontSize',(0.7*get(get(origaxis,'xlabel'),'FontSize')));
        end

        if ~isempty(legend(origaxis)),
            legend(s(j),'toggle');
            set(legend(s(j)),legprops,get(legend(origaxis),legprops));
        end

        j = j+1;
    end
end

for jj = 1:length(FileName),
    close(figurearray(jj),'new','invisible');
end


text:
classdef Kernel < handle
    %{
        Kernel function

            INPUT
        X         data (n*d)
        Y         data (m*d)

        OUTPUT
        K         kernel matrix (n*m)

    ----------------------------------------------------------------------

        type   -

        linear      :  k(x,y) = x'*y
        polynomial  :  k(x,y) = (γ*x'*y+c)^d
        gaussian    :  k(x,y) = exp(-γ*||x-y||^2)
        sigmoid     :  k(x,y) = tanh(γ*x'*y+c)
        laplacian   :  k(x,y) = exp(-γ*||x-y||)

        degree -  d
        offset -  c
        gamma  -  γ

    ----------------------------------------------------------------------

        Version 1.1, 11-MAY-2021
        Email: iqiukp@outlook.com
    ----------------------------------------------------------------------
    %}

    properties
        type = 'gaussian'
        offset = 0
        gamma = 0.1
        degree = 2
    end

    methods
        function obj = Kernel(varargin)
            inputValue = varargin;
            nParameter = size(inputValue, 2)/2;
            supportedKernelFunc = {'linear', 'gaussian', 'polynomial', 'sigmoid', 'laplacian'};
            for n = 1:nParameter
                parameter = inputValue{(n-1)*2+1};
                value = inputValue{(n-1)*2+2};
                if strcmp(parameter, 'type')
                    if ~any(strcmp(value, supportedKernelFunc))
                    errorText = sprintf([
                        'Unsupported kernel function.\n',...
                        'Use one of these kernel functions:\n', ...
                        'linear, gaussian, polynomial, sigmoid, laplacian.']);
                    error(errorText)
                    end
                end
                obj.(parameter) = value;
            end
        end

        function K = computeMatrix(obj, x, y)
            K = zeros(size(x, 1), size(y, 1));
            % compute the kernel matrix
            switch obj.type
                case 'linear' % linear kernel function
                    K = x*y';

                case 'gaussian' % gaussian kernel function
                    try
                        K = exp(-obj.gamma*pdist2(x, y, 'squaredeuclidean'));
                    catch
                        sx = sum(x.^2, 2);
                        sy = sum(y.^2, 2);
                        xy = 2*x*y';
                        K = exp((bsxfun(@minus, bsxfun(@minus, xy, sx), sy'))*obj.gamma);
                    end

                case 'polynomial' % polynomial kernel function
                    K = (obj.gamma*x*y'+obj.offset).^double(obj.degree);

                case 'sigmoid' % sigmoid kernel function
                    K = tanh(obj.gamma*x*y'+obj.offset);

                case 'laplacian' % laplacian kernel function
                    K = exp(-obj.gamma*pdist2(x, y, 'cityblock'));
            end
        end
    end
end

author: iqiukp

repo: KPCA-MATLAB

save_path: github-repos/MATLAB/iqiukp-KPCA-MATLAB


text:
%--------------------
function [EbN, Es] = PQ_excitCB (X2)

persistent W2 EIN

NF = 2048;
Version = 'Basic';
if (isempty (W2))
    Fs = 48000;
    f = linspace (0, Fs/2, NF/2+1);
    W2 = PQWOME (f);
    [Nc, fc] = PQCB (Version);
    EIN = PQIntNoise (fc);
end

% Allocate storage
XwN2 = zeros (1, NF/2+1);

% Outer and middle ear filtering
Xw2(1,:) = W2 .* X2(1,1:NF/2+1);
Xw2(2,:) = W2 .* X2(2,1:NF/2+1);

% Form the difference magnitude signal
for (k = 0:NF/2)
    XwN2(k+1) = (Xw2(1,k+1) - 2 * sqrt (Xw2(1,k+1) * Xw2(2,k+1)) ...
               + Xw2(2,k+1));
end

% Group into partial critical bands
Eb(1,:) = PQgroupCB (Xw2(1,:), Version);
Eb(2,:) = PQgroupCB (Xw2(2,:), Version);
EbN     = PQgroupCB (XwN2, Version);

% Add the internal noise term => "Pitch patterns"
E(1,:) = Eb(1,:) + EIN;
E(2,:) = Eb(2,:) + EIN;

% Critical band spreading => "Unsmeared excitation patterns"
Es(1,:) = PQspreadCB (E(1,:), Version);
Es(2,:) = PQspreadCB (E(2,:), Version);

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
classdef EchoServer < WebSocketServer
    %ECHOSERVER Summary of this class goes here
    %   Detailed explanation goes here

    properties
    end

    methods
        function obj = EchoServer(varargin)
            %Constructor
            obj@WebSocketServer(varargin{:});
        end
    end

    methods (Access = protected)
        function onOpen(obj,conn,message)
            fprintf('%s\n',message)
        end

        function onTextMessage(obj,conn,message)
            % This function sends an echo back to the client
            conn.send(message); % Echo
        end

        function onBinaryMessage(obj,conn,bytearray)
            % This function sends an echo back to the client
            conn.send(bytearray); % Echo
        end

        function onError(obj,conn,message)
            fprintf('%s\n',message)
        end

        function onClose(obj,conn,message)
            fprintf('%s\n',message)
        end
    end
end


text:
function nav_o = gnss_outage(nav, times_out, ACTION)
% gnss_outage: removes or gets navigation data from navigation data structure
% in GNSS outages periods.
%
% INPUT
%   nav, navigation data structure.
%   times_out, Nx1 vector with outage times where times_out(1) is the start
%     time of 1st GNSS OUTAGE, times_out(2) is the end of 1st GNSS OUTAGE,
%     times_out(3) is the start of the 2nd GNSS OUTAGE, and so on. N should
%     be even.
%   ACTION, 'REMOVE' (default) or 'GET' (string).
%
% OUTPUT
%   nav, same navigation data structure.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Version: 001
% Date:    2021/12/14
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

if (nargin < 3)

    ACTION = 'REMOVE';
end

OUTAGE_N = max(size(times_out)) / 2;

%% INIT DATA

if (strcmp(ACTION,'REMOVE'))

    fprintf('gnss_outage: %d GNSS outages are forced... \n', OUTAGE_N)

    nav_o = nav;

elseif(strcmp(ACTION,'GET'))

    nav_o.t = [];
    nav_o.lat = [];
    nav_o.lon = [];
    nav_o.h = [];
    nav_o.vel = [];

    if (isfield(nav,'roll'))
        nav_o.roll = [];
    end

    if (isfield(nav,'pitch'))
        nav_o.pitch = [];
    end

    if (isfield(nav,'yaw'))
        nav_o.yaw = [];
    end
else

    error('gnss_outage: action is not defined')
end

%% REMOVE OR GET

for k = 1:OUTAGE_N

    if (strcmp(ACTION,'REMOVE'))

        idx  = find(nav_o.t >= times_out(2*k-1), 1, 'first' );
        fdx  = find(nav_o.t <= times_out(2*k), 1, 'last' );

        if(isempty(idx) || isempty(fdx))
            error('gnss_outage: outage empty index')
        end

        if(idx >= fdx)
            error('gnss_outage: indexes inconsistency')
        end

        nav_o.t  (idx:fdx)    = [];
        nav_o.lat(idx:fdx)    = [];
        nav_o.lon(idx:fdx)    = [];
        nav_o.h  (idx:fdx)    = [];
        nav_o.vel(idx:fdx, :) = [];

        if (isfield(nav,'roll'))
            nav_o.roll(idx:fdx)  = [];
        end

        if (isfield(nav,'pitch'))
            nav_o.pitch(idx:fdx) = [];
        end

        if (isfield(nav,'yaw'))
            nav_o.yaw(idx:fdx)   = [];
        end

    elseif(strcmp(ACTION,'GET'))

        idx  = find(nav.t >= times_out(2*k-1), 1, 'first' );
        fdx  = find(nav.t <= times_out(2*k), 1, 'last' );

        if(isempty(idx) || isempty(fdx))
            error('gnss_outage: outage empty index')
        end

        if(idx >= fdx)
            error('gnss_outage: indexes inconsistency')
        end

        nav_o.t     = [nav_o.t   ; nav.t(idx:fdx)] ;
        nav_o.lat   = [nav_o.lat ; nav.lat(idx:fdx)];
        nav_o.lon   = [nav_o.lon ; nav.lon(idx:fdx)];
        nav_o.h     = [nav_o.h   ; nav.h(idx:fdx) ];
        nav_o.vel   = [nav_o.vel ; nav.vel(idx:fdx, :)];

        if (isfield(nav,'roll'))
            nav_o.roll = [nav_o.roll ; nav.roll(idx:fdx)];
        end

        if (isfield(nav,'pitch'))
            nav_o.pitch = [nav_o.pitch ; nav.pitch(idx:fdx)];
        end

        if (isfield(nav,'yaw'))
            nav_o.yaw = [nav_o.yaw ; nav.yaw(idx:fdx)];
        end

    else

        error('gnss_outage: action is not defined')
    end

end

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function optval = PW_Vega_CallPut(S,K,C,r,sigma,T)
    % S = NSim x 1 matrix of simulated prices
    % K = Strike price
    % C = 1 -> Call; C = 0 -> Put

    if(C==1)
        Indicator = (S(:,end)>K);

    else
        Indicator = (S(:,end) <= K);
    end
    Sval = S(:,end) .* ((log(S(:,end)./S(1,1)) - (r+0.5*sigma^2)*T)/sigma);
    optval =exp(-r*T) * mean(Indicator .*Sval );
end


text:
%% Anaylsis_SSVEP
clear all; clc; close all;
%% initialization
DATADIR = 'WHERE\IS\DATA';
%% SSVEP
SSVEPDATA = 'EEG_SSVEP.mat';
STRUCTINFO = {'EEG_SSVEP_train', 'EEG_SSVEP_test'};
SESSIONS = {'session1', 'session2'};
TOTAL_SUBJECTS = 54;

%% INITIALIZATION
FS = 100;

%init
params = {'time', 4;...
'freq' , 60./[5, 7, 9, 11];...
'fs' , FS;...
'band' ,[0.5 40];...
'channel_index', [23:32]; ...
'time_interval' ,[0 4000]; ...
'marker',  {'1','up';'2', 'left';'3', 'right';'4', 'down'}; ...
};

%% validation
for sessNum = 1:length(SESSIONS)
    session = SESSIONS{sessNum};
    fprintf('\n%s validation\n',session);
    for subNum = 1:TOTAL_SUBJECTS
        subject = sprintf('s%d',subNum);
        fprintf('LOAD %s ...\n',subject);

        data = importdata(fullfile(DATADIR,session,subject,SSVEPDATA));

        CNT{1} = prep_resample(data.(STRUCTINFO{1}), FS,{'Nr', 0});
        CNT{2} = prep_resample(data.(STRUCTINFO{2}), FS,{'Nr', 0});
        ACC.SSVEP(subNum,sessNum) = SSVEP_performance(CNT, params);
        fprintf('%d = %f\n',subNum, ACC.SSVEP(subNum,sessNum));
        clear CNT
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function [ y, j, ierror ] = yj_check_gregorian ( y, j )

%*****************************************************************************80
%
%% YJ_CHECK_GREGORIAN checks a Gregorian YJ date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, the YJ date.
%
%    Output, integer IERROR, is 0 if no serious error was found in
%    the date, and 1 otherwise.
%

%
%  Check the year.
%
  [ y, ierror ] = y_check_gregorian ( y );

  if ( ierror ~= 0 )
    return
  end
%
%  Make sure J is not too small or too big.
%
  [ y, j ] = j_borrow_gregorian ( y, j );

  [ y, j ] = j_carry_gregorian ( y, j );

  return
end


text:
function cluster_center = cluster_initialize_1 ( dim_num, point_num, ...
  cluster_num, point )

%*****************************************************************************80
%
%% CLUSTER_INITIALIZE_1 initializes the clusters to data points.
%
%  Discussion:
%
%    The cluster centers are simply chosen to be the first data points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates
%    of the points.
%
%    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the coordinates of the cluster centers.
%
  cluster_center(1:dim_num,1:cluster_num) = point(1:dim_num,1:cluster_num)

  return
end


text:
function so = streamCdr(si)
% STREAMCDR  Returns the tail of a stream
so = forceEval(si{2});


text:
%% triSurfSlice
% Below is a demonstration of the features of the |triSurfSlice| function

%%
clear; close all; clc;

%% Syntax
% |[F,V,C,logicSide,Eb]=triSurfSlice(F1,V1,C1,P,n,snapTolerance);|

%% Description
% This function slices the input mesh (defined by faces F, vertices V and
% optional color data C) using a plane specified by a point P and a normal
% vector to the plane n. Edges crossing the plane are cut at the plane, and
% new points are introduced at the intersection, and new triangles are
% created. Points within snapTolerance from the plane are snapped to the
% plane.

%% Examples

%%
% Plot settings
fontSize=15;

%% Example 1: Slicing a surface
% Create example surface data consisting of faces F and vertices V

[F,V]=stanford_bunny; %Get faces and vertices
C=[]; %No color data used in this example so left empty
snapTolerance=mean(patchEdgeLengths(F,V))/100;
n=vecnormalize([0 0 1]); %Normal direction to plane
P=mean(V,1); %Point on plane

%Slicing surface (note 3rd color data output is supressed)
[Fc,Vc,~,logicSide,Eb]=triSurfSlice(F,V,C,P,n,snapTolerance);

%%
% Visualize slices mesh

cFigure;
subplot(1,2,1); hold on;
hp1=gpatch(Fc(~logicSide,:),Vc,'bw','k',1);
hp2=gpatch(Fc(logicSide,:),Vc,'rw','k',1);
legend([hp1 hp2],{'Surface above plane','Surface below plane'})
axisGeom; axis manual; camlight headligth;
colormap gjet;
set(gca,'FontSize',fontSize);

subplot(1,2,2); hold on;
gpatch(Fc(logicSide,:),Vc,'w','k',1);
gpatch(Fc(~logicSide,:),Vc,'w','none',0.25);
hp1=gpatch(Eb,Vc,'none','b',1,3);
hp2=quiverVec(P,n,50,'k');
legend([hp1 hp2],{'Intersection curve','Plane normal vector'})
axisGeom; axis manual; camlight headligth;
set(gca,'FontSize',fontSize);
gdrawnow;

%% Example 2: Slicing while keeping track of surface color data
% Create example surface data consisting of faces F, vertices V, and
% (optionally) color data C.

[F,V]=stanford_bunny; %Get faces and vertices

%Create example color data
FV=patchCentre(F,V)/5; %Face centre coordinates
C=sin(FV(:,1)); %Optional surface color data to "keep hold off"

snapTolerance=mean(patchEdgeLengths(F,V))/100;
n=vecnormalize([0 0 1]); %Normal direction to plane
P=mean(V,1); %Point on plane

%Slicing surface
[Fc,Vc,Cc,logicSide,Eb]=triSurfSlice(F,V,C,P,n,snapTolerance);

%%
% Visualize color data handling

cFigure;
subplot(1,2,1); hold on;
title('original color data');
gpatch(F,V,C,'k',1);
axisGeom; axis manual; camlight headligth;
set(gca,'FontSize',fontSize);
colormap gjet;

subplot(1,2,2); hold on;
title('color data on "cut" surface');
gpatch(Fc,Vc,Cc,'k',1);
gpatch(Eb,Vc,'none','k',1,3);
axisGeom; axis manual; camlight headligth;
set(gca,'FontSize',fontSize);
gdrawnow;

%% Example 3: Slicing a surface and keeping one side
% Create example surface data consisting of faces F, vertices V

[F,V]=stanford_bunny; %Get faces and vertices
snapTolerance=mean(patchEdgeLengths(F,V))/100;
n=vecnormalize([0 0 1]); %Normal direction to plane
P=mean(V,1); %Point on plane

%Slicing surface
[Fc,Vc,Cc,logicSide,Eb]=triSurfSlice(F,V,[],P,n,snapTolerance);

%Keep only part under plane and clean up
[Fcc,Vcc,indFix]=patchCleanUnused(Fc(logicSide,:),Vc);
Eb_cc=indFix(Eb);

%%
% Visualize slices mesh

cFigure; hold on;
hp1=gpatch(F,V,'kw','none',0.25);
hp2=gpatch(Fcc,Vcc,'rw','k',1);
gpatch(Eb_cc,Vcc,'none','b',1,3);
legend([hp1 hp2],{'Original','Cut surface'})
axisGeom; axis manual; camlight headligth;
colormap gjet;
set(gca,'FontSize',fontSize);
gdrawnow;

%% Example 4: Demo animation for changing orientation
% Visualizing slicing operation for varying angles

hf=cFigure;
subplot(1,2,1); hold on;
hp4=gpatch(Fc,Vc,logicSide,'k',1);
axisGeom; axis manual; camlight headligth;
colormap gjet;

subplot(1,2,2); hold on;
hp1=gpatch(Fc(logicSide,:),Vc,'w','k',1);
hp2=gpatch(Fc(~logicSide,:),Vc,'w','none',0.25);
hp3=gpatch(Eb,Vc,'none','b',1,3);
hp5=quiverVec(P,n,25,'k');
axisGeom; axis manual; camlight headligth;
gdrawnow;

vn=hp5.Vertices;

nSteps=50; %Number of animation steps

%Create the time vector
animStruct.Time=linspace(0,1,nSteps);

%The vector lengths
a=linspace(0,2*pi,nSteps);
b=linspace(0,2*pi,nSteps);
for q=1:1:nSteps
    R=euler2DCM([a(q) b(q) 0]);
    nn=n*R;
    [Fc,Vc,~,logicSide,Eb]=triSurfSlice(F,V,C,P,nn,snapTolerance);

    vnn=(vn-P)*R+P;
    %Set entries in animation structure
    animStruct.Handles{q}=[hp1 hp1 hp2 hp2 hp3 hp3 hp4 hp4 hp4 hp5]; %Handles of objects to animate
    animStruct.Props{q}={'Vertices','Faces',...
                         'Vertices','Faces',...
                         'Vertices','Faces',...
                         'Vertices','Faces','CData',...
                         'Vertices'...
                         }; %Properties of objects to animate
    animStruct.Set{q}={Vc,Fc(logicSide,:),Vc,Fc(~logicSide,:),Vc,Eb,Vc,Fc,double(logicSide),vnn}; %Property values for to set in order to animate
end
anim8(hf,animStruct);

%% Example 5: Demo animation illustrating "sharp triangle" problem

[F,V]=quadBox([1 1 1],[5 5 5]);
[F,V]=quad2tri(F,V,'x');
FV=patchCentre(F,V)/5;
C=sin(FV(:,1))+sin(FV(:,2))+sin(FV(:,3));
snapTolerance=mean(patchEdgeLengths(F,V))/50; %Snapping tolerance
n=vecnormalize([0 0 1]); %Plane normal vector

P=mean(V,1); %Point on plane

[Fc,Vc,Cc,logicSide]=triSurfSlice(F,V,C,P,n,snapTolerance);
Eb=patchBoundary(Fc(logicSide,:));

%%
% Visualizing slicing operation and sharp triangles. Animate this figure to
% view effect.

hf=cFigure;
subplot(1,2,1); hold on;
hp4=gpatch(Fc,Vc,logicSide,'k',1,3);
axisGeom; axis manual; camlight headligth;
colormap gjet; axis off;
view(-90,0); zoom(1.25);

subplot(1,2,2); hold on;
hp1=gpatch(Fc(logicSide,:),Vc,'w','k',1,3);
hp2=gpatch(Fc(~logicSide,:),Vc,'w','none',0.25);
hp3=gpatch(Eb,Vc,'none','b',1,4);

axisGeom; axis manual; camlight headligth;
view(-90,0); zoom(1.25);axis off;
gdrawnow;

nSteps=75; %Number of animation steps
animStruct.Time=linspace(0,1,nSteps); %Create the time vector
z=linspace(min(V(:,3)),max(V(:,3)),nSteps);
for q=1:1:nSteps
    P=[0 0 z(q)];
    [Fc,Vc,~,logicSide,Eb]=triSurfSlice(F,V,C,P,n,snapTolerance);

    %Set entries in animation structure
    animStruct.Handles{q}=[hp1 hp1 hp2 hp2 hp3 hp3 hp4 hp4 hp4]; %Handles of objects to animate
    animStruct.Props{q}={'Vertices','Faces',...
                         'Vertices','Faces',...
                         'Vertices','Faces',...
                         'Vertices','Faces','CData'...
                         }; %Properties of objects to animate
    animStruct.Set{q}={Vc,Fc(logicSide,:),Vc,Fc(~logicSide,:),Vc,Eb,Vc,Fc,double(logicSide)}; %Property values for to set in order to animate
end
anim8(hf,animStruct);

%% Example 6: Cutting non-closed surfaces

[F,V]=graphicsModels(10); %Get faces and vertices

snapTolerance=mean(patchEdgeLengths(F,V))/100;
n=vecnormalize([0 0 1]); %Normal direction to plane
P=mean(V,1); %Point on plane

%Slicing surface
[Fc,Vc,Cc,logicSide,Eb]=triSurfSlice(F,V,C,P,n,snapTolerance);

%%
% Visualize slices mesh

cFigure; hold on;
gpatch(Fc,Vc,logicSide,'k',1);
gpatch(Eb,Vc,'none','g',1,3);
axisGeom; axis manual; camlight headligth;
colormap gjet;
set(gca,'FontSize',fontSize);
gdrawnow;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = ihashsum(varargin)
% VL_IHASHSUM
[varargout{1:nargout}] = vl_ihashsum(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
function test_failed=test_dft
Lr=[1, 19, 20];

test_failed=0;

disp(' ===============  TEST_DFT ==============');

for jj=1:length(Lr)
  L=Lr(jj);
    for n = 1:2

    if (n==1)
       type = 'complex';
       f=tester_crand(L,1);
    elseif (n==2)
       type = 'real';
       f=tester_rand(L,1);
    end

    c1=dft(f);
    c2=ref_dft(f);

    res=norm(c1-c2);
    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf('DFT %6s  L:%3i %0.5g %s',type,L,res,fail);
    disp(s);
    end
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% This file plot a summary of the seismicity on one volcano
% on one wheet of paper.

% Stefan Wiemer   08/96
report_this_filefun(mfilename('fullpath'));

% lets start with a map

[existFlag,figNumber]=figure_exists('Seismicity Summary',1);
newMapWindowFlag=~existFlag;

% Set up the Seismicity Map window Enviroment
%
if newMapWindowFlag
    mapI = figure_w_normalized_uicontrolunits( ...
        'Name','Seismicity Summary',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'backingstore','on',...
        'NextPlot','add', ...
        'Visible','off', ...
        'Position',[ (fipo(3:4)-[600 600]), (ZmapGlobal.Data.map_len - [200 80])]);

    stri1 = [file1];

    matdraw

else
end

% show the figure
%
figure_w_normalized_uicontrolunits(mapI)
reset(gca)
delete(gca);delete(gca);delete(gca);
delete(gca);delete(gca);delete(gca);
delete(gca);delete(gca);delete(gca);
watchon;
set(gca,'visible','off','SortMethod','childorder')
hold off

%set(set_ni3,'String',num2str(ni));
% find min and Maximum axes points
s1_east = max(a.Longitude);
s2_west = min(a.Longitude);
s3_north = max(a.Latitude);
s4_south = min(a.Latitude);
%ni = 100;
orient landscape
set(gcf,'PaperPosition',[ 1.0 1.0 8 6])
rect = [0.55,  0.50, 0.25, 0.30];
axis('equal')
axes('position',rect)
%
% find start and end time of catalogue "a"
%
t0b = min(a.Date);
n = a.Count;
teb = a(n,3) ;
tdiff =round(teb - t0b)*365/par1;

n = length(a);

% plot earthquakes (differnt colors for varous depth layers) as
% defined in "startzmap"
%
hold on

dep1 = -1;
dep2 = 1;
dep3 = 4;
deplo1 =plot(a(a.Magnitude<=dep1,1),a(a.Magnitude<=dep1,2),'.b');
set(deplo1,'MarkerSize',ms6,'Marker',ty1,'era','normal')
deplo2 =plot(a(a.Magnitude<=dep2&a.Magnitude>dep1,1),a(a.Magnitude<=dep2&a.Magnitude>dep1,2),'.g');
set(deplo2,'MarkerSize',ms6,'Marker',ty2,'era','normal');
deplo3 =plot(a(a.Magnitude<=dep3&a.Magnitude>dep2,1),a(a.Magnitude<=dep3&a.Magnitude>dep2,2),'.r');
set(deplo3,'MarkerSize',ms6,'Marker',ty3,'era','normal')
ls1 = sprintf('M < %3.1f ',dep1);
ls2 = sprintf('M < %3.1f ',dep2);
ls3 = sprintf('M < %3.1f ',dep3);

%le =legend('+b',ls1,'og',ls2,'xr',ls3);
%set(le,'position',[ 0.75 0.50 0.12 0.07],'FontSize',6)
axis([ s2_west s1_east s4_south s3_north])
%xlabel('Longitude [deg]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
%ylabel('Latitude [deg]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
%strib = [  ' Map of   '  name '; '  num2str(t0b) ' to ' num2str(teb) ];
%title2(strib,'FontWeight','normal',...
% 'FontSize',ZmapGlobal.Data.fontsz.s,'Color','k')

h1 = gca;
 set(gca,'Color',color_bg);

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'normal','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

% next we plot a magnitude stem-plot
rect = [0.15, 0.10, 0.30, 0.20];
axes('position',rect)
h2 = gca;

hs = stem(a.Date,a.Magnitude);
set(hs,'MarkerSize',4)
hold on

set(gca,'XLIM',[min(a.Date) max(a.Date)+0.01])
xl = get(gca,'Xlim');

xlabel('Time in Years ]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.m)
ylabel('Magnitude','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.m)

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'normal','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

 set(gca,'Color',color_bg);
%grid
hold off

% how about a histogram next

rect = [0.15, 0.30, 0.30, 0.15];
axes('position',rect)
h3 = gca;

[n,x] =histogram(a.Date,(min(a.Date):par1:max(a.Date)));
fillbar(x,n,'k')
ylabel('# per day','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.m)

set(gca,'XTickLabels',[])
set(gca,'Xlim',xl)
 set(gca,'Color',color_bg);

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'normal','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
%grid

% now a time-depth plot
rect = [0.15, 0.60, 0.30, 0.20];
axes('position',rect)
h4 = gca;

deplo1 =plot(a(a.Depth<=dep1,3),-a(a.Depth<=dep1,7),'.b');
set(deplo1,'MarkerSize',ms6,'Marker',ty1,'era','normal')
hold on
deplo2 =plot(a(a.Depth<=dep2&a.Depth>dep1,3),-a(a.Depth<=dep2&a.Depth>dep1,7),'.g') ;
set(deplo2,'MarkerSize',ms6,'Marker',ty2,'era','normal');
deplo3 =plot(a(a.Depth<=dep3&a.Depth>dep2,3),-a(a.Depth<=dep3&a.Depth>dep2,7),'.r') ;
set(deplo3,'MarkerSize',ms6,'Marker',ty3,'era','normal')

hold on

ylabel('Depth in [km] ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.m)
%grid
set(gca,'XTickLabels',[])
set(gca,'XLim',xl)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'normal','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

 set(gca,'Color',color_bg);

% next a cumulative number plot

rect = [0.15, 0.45, 0.30, 0.15];
axes('position',rect)
hold on
h3 = gca;

[n,x] =histogram(a.Date,(min(a.Date):par1:max(a.Date)));
fillbar(x,cumsum(n),'k')
ylabel('Cumulative # ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.m)
set(gca,'Xlim',xl)

set(gca,'XTickLabels',[])
set(gca,'Xlim',xl)
 set(gca,'Color',color_bg);

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'normal','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
%grid

% now plot a b-value distribution

rect = [0.55, 0.15, 0.25, 0.25];
axes('position',rect)

newcat =a;

avob

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'normal','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

 set(gca,'Color',color_bg);

watchoff
set(gcf,'PaperPosition',[ 0.1 0.1 11 9])
rect=[0 0 1 1];
h2 =axes('position',rect);
set(h2,'visible','off');
l = a.Count;
txt1=text(.15, .85,[' ' num2str(floor(a(1,3))) '/'  num2str(a(1,4))  '/' ,...
    num2str(a(1,5)) ' - '  num2str(floor(a(l,3))) '/'  num2str(a(l,4))  '/'  num2str(a(l,5)) ]);
set(txt1,'FontWeight','bold','FontSize',14);
c = fix(clock);
txt1=text(.15, .82,['created: ' date '  ' num2str(c(4)) '.' num2str(c(5))  ] )   ;
set(txt1,'FontWeight','normal','FontSize',12)


text:
function [dip] = ft_datatype_dip(dip, varargin)

% FT_DATATYPE_DIP descripts the FieldTrip MATLAB structure for dip data
%
% The dip structure reprents a dipole model that has been fitted to
% ERP or ERF data using a non-linear optimization approach. It is
% usually generated by the FT_DIPOLEFITTING function.
%
% FIXME more information should be added here
%
% See also FT_DATATYPE, FT_DATATYPE_SOURCE, FT_DATATYPE_VOLUME

% Copyright (C) 2011, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isempty(dip)
  return;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function hyperSaveFigure(h, filename, fmt)
% HYPERSAVEFIGURE Writes a figure to disk as an image.
%
% Usage
%   hyperSaveFigure(gcf, 'filename.png');
% Inputs
%   h - Handle to figure
%   filename - Filename for output file.  Extension determines image type.
%   fmt - Format of output image.  'wysiwyg' or wysiwyp' for 'what you see
%         is what you get' and what 'you see is what you print' respectively.
% Outputs
%   none

if (nargin == 2)
    fmt = 'wysiwyg';
end

fmt = lower(fmt);

if strcmp(fmt, 'wysiwyp')
    saveas(h, filename);
elseif strcmp(fmt, 'wysiwyg')
    set(h, 'Color',[1 1 1]);
    frame = getframe(h);
    [X,map] = frame2im(frame);
    imwrite(X ,filename);
else
    error('Bad format string specified.');
end
%print(h, '-dpng', filename);

author: davidkun

repo: HyperSpectralToolbox


text:
%DEMO_MULTICLASS  Classification problem demonstration for 3 classes
%                 using Gaussian process prior
%
%  Description
%    The data used in the demonstration program is the same used by
%    Radford M. Neal in his three-way classification example in
%    Software for Flexible Bayesian Modeling
%    (http://www.cs.toronto.edu/~radford/fbm.software.html) The
%    data consists of 1000 4-D vectors which are classified into
%    three classes. The data is generated by drawing the components
%    of vector, x1, x2, x3 and x4, uniformly form (0,1). The class
%    of each vector is selected according to the first two
%    components of the vector, x_1 and x_2. After this a Gaussian
%    noise with standard deviation of 0.1 has been added to every
%    component of the vector. Because there are two irrelevant
%    components in the input vector a prior with ARD should be of
%    help.
%
%    The data is divided into two parts, trainig set of 400 units
%    and test set of 600 units.
%
%    The latent values for N training points and C classes are
%    f=(f1_1,f2_1,...,fN_1,f1_2,f2_2,...,fN_2,...,f1_C,f2_C,...,fN_C)^T,
%    and are given a zero mean Gaussian process prior
%
%      f ~ N(0, K),
%
%    where K is a block diagonal covariance matrix with blocks
%    K_1,...,K_C whose elements are given by K_ij = k(x_i, x_j |
%    th). The function k(x_i, x_j | th) is covariance function and
%    th its parameters.
%
%    In this demo we approximate the posterior distribution with
%    Laplace approximation.
%

% Copyright (c) 2010 Jaakko Riihim�ki, Jarno Vanhatalo, 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.

% Load the data
S = which('demo_multiclass');
L = strrep(S,'demo_multiclass.m','demodata/cdata.txt');
x=load(L);
y=repmat(0,size(x,1),3);
y(x(:,5)==0,1) = 1;
y(x(:,5)==1,2) = 1;
y(x(:,5)==2,3) = 1;
x(:,end)=[];

% Divide the data into training and test parts.
xt = x(401:end,:);
x=x(1:400,:);
yt=y(401:end,:);
y=y(1:400,:);

[n, nin] = size(x);

% Create covariance functions
gpcf1 = gpcf_sexp('lengthScale', ones(1,nin), 'magnSigma2', 1);
% Set the prior for the parameters of covariance functions
pl = prior_t('s2',10,'nu',10);
pm = prior_sqrtt('s2',10,'nu',10);
gpcf1 = gpcf_sexp(gpcf1, 'lengthScale_prior', pl,'magnSigma2_prior', pm);

% Create the GP structure
gp = gp_set('lik', lik_softmax, 'cf', gpcf1, 'jitterSigma2', 1e-2);

% ------- Laplace approximation --------
fprintf(['Softmax model with Laplace integration over the latent\n' ...
         'values and MAP estimate for the parameters\n'])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'Laplace');

% We could also use own covariance function for each output the following
% way
% gp2 = gp_set('lik', lik_softmax, 'cf', {gpcf1 gpcf1 gpcf1}, 'jitterSigma2', 1e-2);
% gp2 = gp_set(gp2, 'latent_method', 'Laplace');
% gp2.comp_cf = {1 2 3};
% [Eft2, Varft2, ~, ~, pyt2] = gp_pred(gp2, x, y, xt, 'yt', ones(size(yt)));

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3,'Display','iter');
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% make the prediction for test points
[Eft, Varft, lpyt] = gp_pred(gp, x, y, xt, 'yt', ones(size(yt)));
Eft = reshape(Eft, size(xt,1), size(yt,2));
lpyt = reshape(lpyt, size(xt,1), size(yt,2));

% calculate the percentage of misclassified points
tt = exp(lpyt)==repmat(max(exp(lpyt),[],2),1,size(exp(lpyt),2));
missed = (sum(sum(abs(tt-yt)))/2)/size(yt,1)

% grid for making prediction
xtg1 = meshgrid(linspace(min(x(:,1))-.1, max(x(:,1))+.1, 30));
xtg2 = meshgrid(linspace(min(x(:,2))-.1, max(x(:,2))+.1, 30))';
xtg=[xtg1(:) xtg2(:) repmat(mean(x(:,3:4)), size(xtg1(:),1),1)];

[Eft, Covft, pg] = gp_pred(gp, x, y, xtg, 'yt', ones(size(xtg,1),3));
Eft = reshape(Eft, size(xtg,1), 3);
pg=reshape(pg,size(xtg,1),3);

% plot the train data o=0, x=1
figure, set(gcf, 'color', 'w'), hold on
plot(x(y(:,1)==1,1),x(y(:,1)==1,2),'ro', 'linewidth', 2);
plot(x(y(:,2)==1,1),x(y(:,2)==1,2),'x', 'linewidth', 2);
plot(x(y(:,3)==1,1),x(y(:,3)==1,2),'kd', 'linewidth', 2);
axis([-0.4 1.4 -0.4 1.4])
contour(xtg1, xtg2, reshape(exp(pg(:,1)),30,30),'r', 'linewidth', 2)
contour(xtg1, xtg2, reshape(exp(pg(:,2)),30,30),'b', 'linewidth', 2)
contour(xtg1, xtg2, reshape(exp(pg(:,3)),30,30),'k', 'linewidth', 2)

% MCMC approach

% Set the approximate inference method
% Note that MCMC for latent values requires often more jitter
lat = gp_pred(gp, x, y, x);
gp = gp_set(gp, 'latent_method', 'MCMC', 'jitterSigma2', 1e-4);
gp = gp_set(gp, 'latent_opt', struct('method',@scaled_mh));
gp.latentValues = lat(:);

gp_e(gp_pak(gp), gp, x,y)
gp_g(gp_pak(gp), gp, x,y)
gradcheck(randn(size(gp_pak(gp))), @gp_e, @gp_g, gp, x, y);

% Set the parameters for MCMC...
hmc_opt.steps=10;
hmc_opt.stepadj=0.001;
hmc_opt.nsamples=1;
latent_opt.display=0;
latent_opt.repeat = 20;
latent_opt.sample_latent_scale = 0.05;
hmc2('state', sum(100*clock))

% Sample
[r,g,opt]=gp_mc(gp, x, y, 'hmc_opt', hmc_opt, 'latent_opt', latent_opt, 'nsamples', 1, 'repeat', 15);

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

% Sample
[rgp,g,opt]=gp_mc(gp, x, y, 'nsamples', 400, 'hmc_opt', hmc_opt, 'latent_opt', latent_opt, 'record', r);
% Remove burn-in
rgp=thin(rgp,102);

% Make predictions
[Efs_mc, Varfs_mc, pgs_mc] = gpmc_preds(rgp, x, y, xtg, 'yt', ones(size(xtg,1),3));

Ef_mc = reshape(mean(Efs_mc,2),900,3);
pg_mc = reshape(mean(exp(pgs_mc),2),900,3);

figure, set(gcf, 'color', 'w'), hold on
plot(x(y(:,1)==1,1),x(y(:,1)==1,2),'ro', 'linewidth', 2);
plot(x(y(:,2)==1,1),x(y(:,2)==1,2),'x', 'linewidth', 2);
plot(x(y(:,3)==1,1),x(y(:,3)==1,2),'kd', 'linewidth', 2);
axis([-0.4 1.4 -0.4 1.4])
contour(xtg1, xtg2, reshape(pg_mc(:,1),30,30),'r', 'linewidth', 2)
contour(xtg1, xtg2, reshape(pg_mc(:,2),30,30),'b', 'linewidth', 2)
contour(xtg1, xtg2, reshape(pg_mc(:,3),30,30),'k', 'linewidth', 2)

% With scaled HMCS

gp2 = gp_set(gp, 'latent_opt', struct('method',@scaled_hmc));
gp2.latentValues = lat(:);

% Set the parameters for MCMC...
hmc_opt.steps=10;
hmc_opt.stepadj=0.001;
hmc_opt.nsamples=1;

% latent opt
latent_opt.nsamples=1;
latent_opt.nomit=0;
latent_opt.persistence=0;
latent_opt.repeat=20;
latent_opt.steps=20;
latent_opt.stepadj=0.15;
latent_opt.window=5;
hmc2('state', sum(100*clock))

% Here we make an initialization with
% slow sampling parameters

[rgp2,gp2,opt]=gp_mc(gp2, x, y, 'hmc_opt', hmc_opt, 'latent_opt', latent_opt, 'nsamples', 1, 'repeat', 15);

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

% Sample
[rgp2,g2,opt2]=gp_mc(gp2, x, y, 'nsamples', 200, 'hmc_opt', hmc_opt, 'latent_opt', latent_opt, 'record', rgp2);
% Remove burn-in
rgp2=thin(rgp2,102);

% Make predictions
[Efs_mc, Varfs_mc, pgs_mc] = gpmc_preds(rgp2, x, y, xtg, 'yt', ones(size(xtg,1),3));

Ef_mc2 = reshape(mean(Efs_mc,2),900,3);
pg_mc2 = reshape(mean(exp(pgs_mc),2),900,3);

figure, set(gcf, 'color', 'w'), hold on
plot(x(y(:,1)==1,1),x(y(:,1)==1,2),'ro', 'linewidth', 2);
plot(x(y(:,2)==1,1),x(y(:,2)==1,2),'x', 'linewidth', 2);
plot(x(y(:,3)==1,1),x(y(:,3)==1,2),'kd', 'linewidth', 2);
axis([-0.4 1.4 -0.4 1.4])
contour(xtg1, xtg2, reshape(pg_mc2(:,1),30,30),'r', 'linewidth', 2)
contour(xtg1, xtg2, reshape(pg_mc2(:,2),30,30),'b', 'linewidth', 2)
contour(xtg1, xtg2, reshape(pg_mc2(:,3),30,30),'k', 'linewidth', 2)