text:
function figure_handle = PTKGraphContextVsMetric(table, metric, context_list, patient_uids)
    % PTKGraphContextVsMetric. Plots a graph showing measurement values against lung regions for one or more subjects
    %
    % PTKGraphContextVsMetric creates a figure and plots a graph showing measurements for particular lung regions, for one or more patients.
    %
    % The data to plot must be held in a MimResultsTable. You specify which contexts (regions), which metric (measurement) and which patients to plot.
    % The contexts (regions) will be shown along the x-axis, and the values of the metric will be shown on the y-axis.
    %
    % Labels and ticks will be automatically generated from the data.
    %
    % Note the figure is optimised for exporting to a file, not for displaying on-screen. The on-screen figure is not intended for viewing and may have font and graphic sizes in the wrong proportions. To view the figure in its correct proportions, you should export the figure as a graphics file (e.g. png or jpg) and then view the file.
    %
    % Syntax:
    %     figure_handle = PTKGraphContextVsMetric(table, metric, context_list, patient_uids)
    %
    % Inputs:
    %     table - a MimResultsTable containing the data to plot. Only a subset of data will be plotted, determined by the other parameters.
    %     metric - a string containing the id of the metric to plot. This must correspond to the metic id in the table.
    %     context_list - Set to [] to plot all contexts, otherwise specify a set of one or more contexts or context sets to plot. Contexts will appear on the x-axis.
    %     patient_uids - Set to [] to plot all patients, otherwise specify a set of one of more patient UIDs to appear on the graph.
    %         Each patient will appear with different markers and will appear in the legend.
    %
    % Output:
    %     figure_handle - the handle of the generated figure. Use this handle to export the figure to an image file
    %
    %     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.
    %

    if nargin < 4
        patient_uids = table.NameMaps{1}.keys;
    end

    y_label = table.NameMaps{2}(metric);

    figure_title = y_label;

    if numel(patient_uids) == 1
        figure_title = [figure_title, ' : ', table.NameMaps{1}(patient_uids{1})];
    end

    if nargin < 3 || isempty(context_list)
        context_list = table.IndexMaps{3}.keys;
    end

    context_list = SortContexts(context_list);

    metrix_index = table.IndexMaps{2}(metric);

    number_of_subjects = numel(patient_uids);

    colours = {[0, 0.7, 0.7], [0, 0, 1], [1, 0, 0], [0.7, 0, 0.7], [0, 0.7, 0], [0, 0, 0]};
    while numel(colours) < number_of_subjects
        colours = [colours, colours];
    end

    symbols = {'^', 'd', 'p', 'h', 's', 'o', '^', '*', 'v', '+'};
    while numel(symbols) < number_of_subjects
        symbols = [symbols, symbols];
    end

    label_font_size = 10;
    x_tick_label_size = 10;
    y_tick_label_size = 8;
    x_lung_label_size = 8;
    if numel(context_list) > 10
        x_tick_label_size = 5;
    end

    right_lung_colour = [0.7, 0, 0];
    left_lung_colour = [0, 0, 0.7];

    legend_font_size = 5;
    widthheightratio = 4/3;
    page_width_cm = 13;
    font_name = PTKSoftwareInfo.GraphFont;
    marker_line_width = 1;
    context_line_width = 1;
    marker_size = 8;
    legend_marker_size = 6;

    legend_strings = {};

    max_y = 0;
    min_y = [];

    context_labels = {};

    for context_list_index = 1 : numel(context_list)
        if ~isempty(context_list{context_list_index})
            context_text_label = table.NameMaps{3}(char(context_list{context_list_index}));
            if length(context_text_label) > 2
                if strcmp(context_text_label(1:2), 'R_') || strcmp(context_text_label(1:2), 'L_')
                    context_text_label = context_text_label(3:end);
                end
            end
            context_labels{end + 1} = context_text_label;
        end
    end

    figure_handle = figure;
    set(figure_handle, 'Units', 'centimeters');
    graph_size = [page_width_cm, (page_width_cm/widthheightratio)];

    axes_handle = gca;
    set(figure_handle, 'Name', figure_title);
    set(figure_handle, 'PaperPositionMode', 'auto');
    set(figure_handle, 'position', [0,0, graph_size]);
    hold(axes_handle, 'on');

    x_ticks = [];
    max_x = 0;

    plots = [];
    right_lung_segments_label_position = [];
    left_lung_segments_label_position = [];

    x_ticks_color_map = containers.Map('KeyType', 'uint32', 'ValueType', 'any');
    x_ticks_label_map = containers.Map('KeyType', 'uint32', 'ValueType', 'any');

    for patient_iterator = 1 : number_of_subjects
        patient_uid = patient_uids{patient_iterator};
        patient_index = table.IndexMaps{1}(patient_uid);

        x_coords = [];
        y_coords = [];

        x_coord = 1;

        for context_list_index = 1 : numel(context_list)
            this_context = char(context_list{context_list_index});
            if ~isempty(this_context)
                context_colour = [0.1, 0.1, 0.1];
                if length(this_context) > 1
                    if strcmp(this_context(1:2), 'R_')
                        context_colour = right_lung_colour;
                    elseif strcmp(this_context(1:2), 'L_')
                        context_colour = left_lung_colour;
                    end
                end

                if strcmp(this_context, 'R_M')
                    right_lung_segments_label_position = x_coord + 0.5;
                elseif strcmp(this_context, 'L_IL')
                    left_lung_segments_label_position = x_coord;
                end
                context_table_index = table.IndexMaps{3}(char(this_context));
                result = table.ResultsTable{patient_index, metrix_index, context_table_index, :};
                if ~isempty(result)
                    x_coords(end + 1) = x_coord;
                    x_ticks = union(x_ticks, x_coord);
                    y_coords(end + 1) = result;
                    x_ticks_color_map(x_coord) = {context_colour};

                    context_text_label = table.NameMaps{3}(this_context);
                    if length(context_text_label) > 2
                        if strcmp(context_text_label(1:2), 'R_') || strcmp(context_text_label(1:2), 'L_')
                            context_text_label = context_text_label(3:end);
                        end
                    end
                    context_labels{end + 1} = context_text_label;

                    x_ticks_label_map(x_coord) = {context_text_label};
                end
            end
            x_coord = x_coord + 1;
        end

        colour = colours{patient_iterator};
        symbol = symbols{patient_iterator};
        marker_face_colour = min(1, colour + 0.5);

        plots_info = [];
        plots_info.X = x_coords;
        plots_info.Y = y_coords;
        plots_info.Symbol = symbol;
        plots_info.Colour = colour;
        plots_info.MarkerFaceColour = marker_face_colour;

        plots{patient_iterator} = plots_info;

        % Find the maxima and maxima
        if isempty(min_y)
            min_y = min(y_coords(:));
        else
            min_y = min(min_y, min(y_coords(:)));
        end

        max_y = max(max_y, max(y_coords(:)));

        max_x = max(max_x, x_coord);

        legend_strings{patient_iterator} = table.NameMaps{1}(patient_uid);
    end

    [y_tick_spacing, min_y, max_y] = MimGraphUtilities.GetOptimalTickSpacing(min_y, max_y);

    y_ticks = min_y : y_tick_spacing : max_y;

    % Draw lines for each context
    for x_tick_index = 1 : numel(x_ticks)
        x_tick = x_ticks(x_tick_index);
        context_line_colour = x_ticks_color_map(x_tick);

        h_line = line('Parent', axes_handle, 'XData', [x_tick, x_tick], 'YData', [min_y, max_y], 'Color', context_line_colour{1}, 'LineStyle', ':', 'LineWidth', context_line_width);
        set(get(get(h_line,'Annotation'),'LegendInformation'), 'IconDisplayStyle','off'); % Exclude line from legend

        y_offset_ticklabel = abs(max_y - min_y)/60;
        context_label = x_ticks_label_map(x_tick);
        text(x_tick, min_y - y_offset_ticklabel, context_label{1}, 'FontSize', x_tick_label_size, 'clipping', 'off', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'Color', context_line_colour{1});
    end

    % Set the axes
    ylabel(axes_handle, y_label, 'FontName', font_name, 'FontSize', label_font_size);
    set(gca, 'XTick', [])
    set(gca, 'XTickLabel', '');
    set(gca, 'YTick', y_ticks, 'FontSize', y_tick_label_size)
    axis([0 max_x min_y max_y]);

    % Now draw the plots
    for patient_iterator = 1 : number_of_subjects
        plots_info = plots{patient_iterator};
        plot(axes_handle, plots_info.X, plots_info.Y, plots_info.Symbol, 'LineWidth', marker_line_width, 'MarkerEdgeColor', plots_info.Colour, 'Color', plots_info.Colour, 'MarkerSize', marker_size, 'MarkerFaceColor', plots_info.MarkerFaceColour);
    end

    if ~isempty(right_lung_segments_label_position) && ~isempty(left_lung_segments_label_position)
        y_offset = abs(max_y - min_y)/15;
        text(right_lung_segments_label_position, min_y - y_offset, 'Right lung', 'FontSize', x_lung_label_size, 'clipping', 'off', 'HorizontalAlignment', 'center', 'Color', right_lung_colour);
        text(left_lung_segments_label_position, min_y - y_offset, 'Left lung', 'FontSize', x_lung_label_size, 'clipping', 'off', 'HorizontalAlignment', 'center', 'Color', left_lung_colour);
    end

    % Create the legend
    legend_handle = legend(legend_strings, 'FontName', font_name, 'FontSize', legend_font_size, 'Location', 'NorthEast');
    legend_children = get(legend_handle, 'Children');
    for child = legend_children'
        if strcmp(get(child, 'Type'), 'line')
            if ~strcmp(get(child, 'Marker'), 'none')
                set(child, 'MarkerSize', legend_marker_size);
            end
        end
    end
end

function sorted_contexts = SortContexts(contexts)
    contexts_as_chars = [];
    for context = contexts
        contexts_as_chars{end + 1} = char(context);
    end
    contexts = contexts_as_chars;
    sorted_contexts = [];
    last_context_was_null = true;
    for context_set = GetContextLabels
        context = context_set{1};

        % We use null as a separator between groups of contexts
        if isempty(context)
            if ~last_context_was_null
                sorted_contexts{end + 1} = [];
                sorted_contexts{end + 1} = [];
            end
            last_context_was_null = true;
        else
            char_context = char(context);
            if ismember(char_context, contexts)
                sorted_contexts{end + 1} = char_context;
                contexts = setdiff(contexts, char_context);
                last_context_was_null = false;
            end
        end
    end
    sorted_contexts = [sorted_contexts, contexts];
end

function context_labels = GetContextLabels
    context_labels = {
        PTKContext.Lungs, PTKContext.RightLung, PTKContext.LeftLung, [], ...
        PTKContext.RightUpperLobe, PTKContext.RightMiddleLobe, PTKContext.RightLowerLobe, [], ...
        PTKContext.LeftUpperLobe, PTKContext.LeftLowerLobe, [], ...
        PTKContext.R_AP, PTKContext.R_P, PTKContext.R_AN, ... % Right upper lobe
        PTKContext.R_L, PTKContext.R_M, ... % Right middle lobe
        PTKContext.R_S, PTKContext.R_AB, PTKContext.R_LB, PTKContext.R_MB, PTKContext.R_PB, [], ... % Right lower lobe
        PTKContext.L_APP, PTKContext.L_APP2, PTKContext.L_AN, ... % Upper left lobe
        PTKContext.L_SL, PTKContext.L_IL, ... % Lingular of upper left lobe
        PTKContext.L_S, PTKContext.L_AMB, PTKContext.L_PB, PTKContext.L_LB}; % lower left lobe
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [score] = processOneImage_CNN(net,img_path)

I = imread(img_path);
filter_size = 7;
filter_ = gpuArray(fspecial('gaussian',filter_size,1.0));
I = localNormalization(I,filter_);
patchsize = size(net.meta.normalization.averageImage,1);

if strcmp(net.meta.modelType,'DirectTraining')
    step = patchsize / 2;
else
    step = 64;
end
[x1,x2,~] = size(I);
X = 1:step:x1-patchsize;
Y = 1:step:x2-patchsize;
data = zeros(patchsize,patchsize,3,length(X)*length(Y),'single');
cnt = 1;
for p = 1:length(X)
    for q = 1:length(Y)
        data(:,:,:,cnt) = I(X(p):X(p)+patchsize-1,Y(q):Y(q)+patchsize-1,:);
        cnt = cnt+1;
    end
end

data = bsxfun(@minus, single(data), net.meta.normalization.averageImage);
data = gpuArray(data) ;

net.eval({'data', data}) ;
index = find(arrayfun(@(a) strcmp(a.name, 'myLoss'), net.layers)==1);
tlayerName = net.layers(index).inputs{1};
score = net.vars(net.getVarIndex(tlayerName)).value ;
if strcmp(net.layers(end).block.lossType,'CE')
    score = squeeze(gather(vl_nnsoftmax(score)));
else
    score = squeeze(gather(score));
end

if size(score,2) == 1
    score = mean(score);
else
%     score_map = net.meta.classes.weights;
%     score = score_map * score;
    [score, ~,~] = svmpredict([1:size(score,2)]', double(score'), net.meta.svm_model1,'-q');
    score = mean(score);
end

end

author: HuiZeng


text:
% CLIENT connect to a server and read a message
%
% Usage - message = client(host, port, number_of_retries)
function message = client(host, port, number_of_retries)

    import java.net.Socket
    import java.io.*

    if (nargin < 3)
        number_of_retries = 20; % set to -1 for infinite
    end

    retry        = 0;
    input_socket = [];
    message      = [];

    while true

        retry = retry + 1;
        if ((number_of_retries > 0) && (retry > number_of_retries))
            fprintf(1, 'Too many retries\n');
            break;
        end

        try
            fprintf(1, 'Retry %d connecting to %s:%d\n', ...
                    retry, host, port);

            % throws if unable to connect
            input_socket = Socket(host, port);

            % get a buffered data input stream from the socket
            input_stream   = input_socket.getInputStream;
            d_input_stream = DataInputStream(input_stream);

            fprintf(1, 'Connected to server\n');

            % read data from the socket - wait a short time first
            pause(0.5);
            bytes_available = input_stream.available;
            fprintf(1, 'Reading %d bytes\n', bytes_available);

            message = zeros(1, bytes_available, 'uint8');
            for i = 1:bytes_available
                message(i) = d_input_stream.readByte;
            end

            message = char(message);

            % cleanup
            input_socket.close;
            break;

        catch
            if ~isempty(input_socket)
                input_socket.close;
            end

            % pause before retrying
            pause(1);
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function keys=hash_keys(hs)
% keys=hash_keys(hs)
%   Returns cell array of keys of the hash hs
%
%   Input:
%       hs  - hash clas
%
%   Output:
%       keys- cell array

% Dimitar Atanasov, 2010
% datanasov@nbu.bg

K = hs.keys;

keys = {};

for k=1:size(K,2)
    for k1 = K{k}
        keys = [keys k1];
    end;
end;


text:
function [selection,value] = listdlg2(varargin)
%LISTDLG  List selection dialog box.
%   [SELECTION,OK] = LISTDLG2('ListString',S) creates a modal dialog box
%   which allows you to select a string or multiple strings from a list.
%   SELECTION is a vector of indices of the selected strings (length 1 in
%   the single selection mode).  This will be [] when OK is 0.  OK is 1 if
%   you push the OK button, or 0 if you push the Cancel button or close the
%   figure.
%
%   Double-clicking on an item or pressing <CR> when multiple items are
%   selected has the same effect as clicking the OK button.  Pressing <CR>
%   is the same as clicking the OK button. Pressing <ESC> is the same as
%   clicking the Cancel button.
%
%   Inputs are in parameter,value pairs:
%
%   Parameter       Description
%   'ListString'    cell array of strings for the list box.
%   'SelectionMode' string; can be 'single' or 'multiple'; defaults to
%                   'multiple'.
%   'ListSize'      [width height] of listbox in pixels; defaults
%                   to [160 300].
%   'InitialValue'  vector of indices of which items of the list box
%                   are initially selected; defaults to the first item.
%   'Name'          String for the figure's title; defaults to ''.
%   'PromptString'  string matrix or cell array of strings which appears
%                   as text above the list box; defaults to {}.
%   'OKString'      string for the OK button; defaults to 'OK'.
%   'CancelString'  string for the Cancel button; defaults to 'Cancel'.
%
%   A 'Select all' button is provided in the multiple selection case.
%
%   Example:
%     d = dir;
%     str = {d.name};
%     [s,v] = listdlg('PromptString','Select a file:',...
%                     'SelectionMode','single',...
%                     'ListString',str)
 %
%  See also DIALOG, ERRORDLG, HELPDLG, INPUTDLG,
%    MSGBOX, QUESTDLG, WARNDLG.

%   Copyright 1984-2014 The MathWorks, Inc.

%   'uh'            uicontrol button height, in pixels; default = 22.
%   'fus'           frame/uicontrol spacing, in pixels; default = 8.
%   'ffs'           frame/figure spacing, in pixels; default = 8.

% simple test:
%
% d = dir; [s,v] = listdlg2('PromptString','Select a file:','ListString',{d.name});
%

% Generate a warning in -nodisplay and -noFigureWindows mode.
% warnfiguredialog('listdlg');
%  usedpath=strcat(matlabroot,'/toolbox/matlab/uitools/private');
%  addpath(genpath(usedpath));
narginchk(1,inf)

figname = '';
smode = 2;   % (multiple)
promptstring = {};
liststring = [];
listsize = [160 300];
initialvalue = [];
okstring = getString(message('MATLAB:uistring:popupdialogs:OK'));
cancelstring = getString(message('MATLAB:uistring:popupdialogs:Cancel'));
fus = 8;
ffs = 8;
uh = 22;

% custom change
pos = [];

if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(message('MATLAB:listdlg:InvalidArgument'))
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
     case 'name'
      figname = varargin{i+1};
     case 'promptstring'
      promptstring = varargin{i+1};
     case 'selectionmode'
      switch lower(varargin{i+1})
       case 'single'
        smode = 1;
       case 'multiple'
        smode = 2;
      end
     case 'listsize'
      listsize = varargin{i+1};
     case 'liststring'
      liststring = varargin{i+1};
     case 'initialvalue'
      initialvalue = varargin{i+1};
     case 'uh'
      uh = varargin{i+1};
     case 'fus'
      fus = varargin{i+1};
     case 'ffs'
      ffs = varargin{i+1};
     case 'okstring'
      okstring = varargin{i+1};
     case 'cancelstring'
      cancelstring = varargin{i+1};

	 % custom change
	 case 'position'
      pos = varargin{i+1};

	% original code
    otherwise
      error(message('MATLAB:listdlg:UnknownParameter', varargin{ i }))
    end
end

if ischar(promptstring)
    promptstring = cellstr(promptstring);
end

if isempty(initialvalue)
    initialvalue = 1;
end

if isempty(liststring)
    error(message('MATLAB:listdlg:NeedParameter'))
end

ex = get(0,'DefaultUicontrolFontSize')*1.7;  % height extent per line of uicontrol text (approx)

fp = get(0,'DefaultFigurePosition');
w = 2*(fus+ffs)+listsize(1);
h = 2*ffs+6*fus+ex*length(promptstring)+listsize(2)+uh+(smode==2)*(fus+uh);
fp = [fp(1) fp(2)+fp(4)-h w h];  % keep upper left corner fixed

fig_props = { ...
    'windowstyle'            'modal' ...
    'name'                   figname ...
    'color'                  get(0,'DefaultUicontrolBackgroundColor') ...
    'resize'                 'off' ...
    'numbertitle'            'off' ...
    'menubar'                'none' ...
    'visible'                'off' ...
    'createfcn'              ''    ...
    'position'               fp   ...
    'closerequestfcn'        'delete(gcbf)' ...
            };

liststring=cellstr(liststring);

fig = figure(fig_props{:});

if ~isempty(promptstring)
    prompt_text = uicontrol('Style','text','String',promptstring,...
        'HorizontalAlignment','left',...
        'Position',[ffs+fus fp(4)-(ffs+fus+ex*length(promptstring)) ...
        listsize(1) ex*length(promptstring)]); %#ok
end

btn_wid = (fp(3)-2*(ffs+fus)-fus)/2;

listbox = uicontrol('Style','listbox',...
                    'Position',[ffs+fus ffs+uh+4*fus+(smode==2)*(fus+uh) listsize],...
                    'String',liststring,...
                    'BackgroundColor','w',...
                    'Max',smode,...
                    'Tag','listbox',...
                    'Value',initialvalue, ...
                    'Callback', {@doListboxClick});

ok_btn = uicontrol('Style','pushbutton',...
                   'String',okstring,...
                   'Position',[ffs+fus ffs+fus btn_wid uh],...
                   'Tag','ok_btn',...
                   'Callback',{@doOK,listbox});

cancel_btn = uicontrol('Style','pushbutton',...
                       'String',cancelstring,...
                       'Position',[ffs+2*fus+btn_wid ffs+fus btn_wid uh],...
                       'Tag','cancel_btn',...
                       'Callback',{@doCancel,listbox});

if smode == 2
    selectall_btn = uicontrol('Style','pushbutton',...
                              'String',getString(message('MATLAB:uistring:popupdialogs:SelectAll')),...
                              'Position',[ffs+fus 4*fus+ffs+uh listsize(1) uh],...
                              'Tag','selectall_btn',...
                              'Callback',{@doSelectAll, listbox});

    if length(initialvalue) == length(liststring)
        set(selectall_btn,'Enable','off')
    end
    set(listbox,'Callback',{@doListboxClick, selectall_btn})
end

set([fig, ok_btn, cancel_btn, listbox], 'KeyPressFcn', {@doKeypress, listbox});

% custom change
if isempty(pos)
   set(fig,'Position',getnicedialoglocation(fp, get(fig,'Units'))); % only this line existed in original code
else
   set(fig, 'Position', pos); % set the figure position as desired using custom code
end

% Make ok_btn the default button.
setdefaultbutton(fig, ok_btn);

% make sure we are on screen
movegui(fig)
set(fig, 'Visible','on'); drawnow;

try
    % Give default focus to the listbox *after* the figure is made visible
    uicontrol(listbox);
    c = matlab.ui.internal.dialog.DialogUtils.disableAllWindowsSafely();
    uiwait(fig);
    delete(c);
catch
    if ishghandle(fig)
        delete(fig)
    end
end

if isappdata(0,'ListDialogAppData__')
    ad = getappdata(0,'ListDialogAppData__');
    selection = ad.selection;
    value = ad.value;
    rmappdata(0,'ListDialogAppData__')
else
    % figure was deleted
    selection = [];
    value = 0;
end
drawnow; % Update the view to remove the closed figure (g1031998)

%% figure, OK and Cancel KeyPressFcn
function doKeypress(src, evd, listbox) %#ok
switch evd.Key
 case 'escape'
  doCancel([],[],listbox);
end

%% OK callback
function doOK(ok_btn, evd, listbox) %#ok
if (~isappdata(0, 'ListDialogAppData__'))
    ad.value = 1;
    ad.selection = get(listbox,'Value');
    setappdata(0,'ListDialogAppData__',ad);
    delete(gcbf);
end

%% Cancel callback
function doCancel(cancel_btn, evd, listbox) %#ok
ad.value = 0;
ad.selection = [];
setappdata(0,'ListDialogAppData__',ad)
delete(gcbf);

%% SelectAll callback
function doSelectAll(selectall_btn, evd, listbox) %#ok
set(selectall_btn,'Enable','off')
set(listbox,'Value',1:length(get(listbox,'String')));

%% Listbox callback
function doListboxClick(listbox, evd, selectall_btn) %#ok
% if this is a doubleclick, doOK
if strcmp(get(gcbf,'SelectionType'),'open')
    doOK([],[],listbox);
elseif nargin == 3
    if length(get(listbox,'String'))==length(get(listbox,'Value'))
        set(selectall_btn,'Enable','off')
    else
        set(selectall_btn,'Enable','on')
    end
end

author: LimingShi


text:
function [gradient matrix] = pseudoGradient(skel, channels, gradientDirection, ind,rotInd2, ind2Compute)

% PSEUDOGRADIENT ths position function of a node is composed with 2
% terms. The first one is the position of the parent node, the second is an
% additional term to get the real position. This function computes the
% derivative of the second term for the node 'ind2Compute' in respect with
% the channels 'rotInd2' given a set of channels.
%
%	Description:
%
%	[gradient matrix] = PSEUDOGRADIENT(skel, channels,gradientDirection, ind,rotInd2, ind2Compute)
%
%	 Returns:
%	  GRADIENT - the coponents of the gradient (x,y,z)
%     MATRIX - the rotation matrix
%	 Arguments:
%	  SKEL - a skeleton structure
%	  CHANNELS - the channels where the gradient has to be computes
%     GRADIENTDIRECTION - the direction of the gradient
%     IND - a recursive paremeter... start with 0 in general
%     ROTIND2 - the gradient has to be computed in respect with this
%     channel
%     IND2COMPUTE - the node ID where the gradient has to be computed

if (ind == 0)
    ind = ind2Compute;
end
matrix =eye(3);
parent = skel.tree(ind).parent;

rotVal = zeros(1, 3);
derive =0 ;
for j = 1:length(skel.tree(ind).rotInd)
    rind = skel.tree(ind).rotInd(j);
    if (rind == rotInd2)
        derive = 1;
    end
    if rind
        rotVal(j) = channels(rind);
    else
        rotVal(j) = 0;
    end
end

if derive
    tdof = rotationMatrixGradient(deg2rad(rotVal(1)), ...
        deg2rad(rotVal(2)), ...
        deg2rad(rotVal(3)), ...
        skel.tree(ind).order,gradientDirection);

else
    tdof = rotationMatrix(deg2rad(rotVal(1)), ...
        deg2rad(rotVal(2)), ...
        deg2rad(rotVal(3)), ...
        skel.tree(ind).order);

end

torient = rotationMatrix(deg2rad(skel.tree(ind).axis(1)), ...
    deg2rad(skel.tree(ind).axis(2)), ...
    deg2rad(skel.tree(ind).axis(3)), ...
    skel.tree(ind).axisOrder);
torientInv = rotationMatrix(deg2rad(-skel.tree(ind).axis(1)), ...
    deg2rad(-skel.tree(ind).axis(2)), ...
    deg2rad(-skel.tree(ind).axis(3)), ...
    skel.tree(ind).axisOrder(end:-1:1));

if (parent~=0)

   [useless matrixTmp] = pseudoGradient(skel, channels, gradientDirection, parent,rotInd2, ind2Compute);
    matrix = torientInv*tdof*torient*matrixTmp;
else

    derive = 0;
    rotVal = skel.tree(1).orientation;
    for i = 1:length(skel.tree(1).rotInd)
        rind = skel.tree(1).rotInd(i);

        if (rind == rotInd2)
            derive = 1;

        end
        if rind
            rotVal(i) = rotVal(i) + channels(rind);
        end
    end
    if derive

        matrix2 = rotationMatrixGradient(deg2rad(rotVal(1)), ...
            deg2rad(rotVal(2)), ...
            deg2rad(rotVal(3)), ...
            skel.tree(1).order,gradientDirection);

    else

        matrix2 = rotationMatrix(deg2rad(rotVal(1)), ...
            deg2rad(rotVal(2)), ...
            deg2rad(rotVal(3)), ...
            skel.tree(1).axisOrder);

    end

matrix = matrix * matrix2;
end
gradient = skel.tree(ind2Compute).offset*matrix;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
 function [xs, ni] = tml_sps(x, Gt, yi, bi, ri, niter, pixmax, curv)
%function [xs, ni] = tml_sps(x, Gt, yi, bi, ri, niter, pixmax, curv)
%	One iteration of the ML-SPS algorithm for transmission Poisson problem
%	(separable paraboloidal surrogates)
%	model: Y_i ~ Poisson(b_i exp(-[G x]_i) + r_i)
%	Input
%		x	[np,1]	initial guess
%		Gt		transpose of system matrix
%		yi		transmission sinogram
%		bi		blank scan factors
%		ri		background (randoms, scatter, crosstalk, etc)
%		bi,ri:		optional (can use empty matrices)
%		yi,bi,ri	must have identical dimensions
%		pixmax		upper constraint for pixel values
%			can be scalar (e.g. 'inf') or an array the size of x
%		curv		'oc' for erdogan's optimal curvatures
%				'pc' for erdogan's fast precomputed curvatures,
%					which usually gives faster convergence,
%					but can be nonmonotone
%				'nr' newton curvatures, can be nonmonotone
%	Output
%		x [np,niter]	updated image vectors each iteration
%
%	fix: the really slick way to do this would be to use
%	the fast precomputed denominator and just backtrack to the monotone
%	version on those rare occasions when it goes downhill
%
%	Copyright 2000-3-01	Jeff Fessler	The University of Michigan

if nargin < 3, ir_usage, end

[nb, na] = size(yi);

if (nargin < 4 || isempty(bi))
	bi = ones(size(yi));
end
if (nargin < 5 || isempty(ri))
	ri = zeros(size(yi));
end
if (nargin < 6 || isempty(niter))
	niter = 2;
end
if (nargin < 7 || isempty(pixmax))
	pixmax = inf;
end
if (nargin < 8 || isempty(curv))
	curv = 'oc';
end

trl_check(yi, bi, ri);

	gi = sum(Gt)';	% g_i = sum_j g_ij

	if strcmp(curv, 'pc')
		ni = trl_curvature_pre(yi, bi, ri);	% precomputed
		denom = Gt * (gi .* ni(:));
	elseif strcmp(curv, 'nr')
		warning 'newton curvatures can be non-monotone'
	elseif ~strcmp(curv, 'oc')
		error 'curv not implemented'
	end

xs = zeros(length(x), niter);
x = max(x,0);
x = min(x,pixmax);
xs(:,1) = x;

%
%	loop over iterations
%
for ii=2:niter
	li = reshape(Gt' * x, size(yi));	% l=G*x "line integrals"
	bel = bi .* exp(-li);
	yb = bel + ri;			% predicted measurement means

	dothi = (1 - yi ./ yb) .* bel;

	if strcmp(curv, 'oc')
		%	optimal curvatures (for ensured monotone increase)
		ni = trl_curvature(yi, bi, ri, li, 'oc');
		denom = Gt * (gi .* ni(:));
	elseif strcmp(curv, 'nr')
		ni = (1 - ri.*yi./yb.^2) .* bel;
		denom = Gt * (gi .* ni(:));
	end

	x = x + (Gt * dothi(:)) ./ denom;	% there's the update!
	x = max(x,0);				% enforce nonnegativity
	x = min(x,pixmax);			% enforce upper bound constraint

	xs(:,ii) = x;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%--- help for checksiz ---
%
% CHECKSIZ Test size consistency among input arguments.
%
%    [ECODE, ERRMSG] = checksiz(SIZES, FUN)
%
%    Inputs:
%     SIZES - An Mx2 matrix containing the size of each input argument.
%
%       FUN - The name of the function from which CHECKSIZ is invoked.
%
%    Outputs:
%     ECODE - Specifies if any size consistencies are violated.
%             1 - Yes
%             0 - No
%
%    ERRMSG - Displays an error message as an output (instead of to the command
%                  window) if boundary conditions are not met.
%
%    For example:
%       To check for input argument size consistency in the function FOO which
%       takes the inputs A and B.
%
%       A = rand(5); B = rand(7);
%       [ecode, errMsg] = checksiz([size(A); size(B)], 'FOO')
%
%       ecode =
%            1
%
%       errMsg =
%       Dimensions of inputs are inconsistent.
%
%    See also CHECKRNG, CHECKTYP.
%

author: jmaih


text:
function rec = VOCreadrecxml(path)

x=VOCreadxml(path);
x=x.annotation;

rec.folder=x.folder;
rec.filename=x.filename;
rec.source.database=x.source.database;
rec.source.annotation=x.source.annotation;
rec.source.image=x.source.image;

rec.size.width=str2double(x.size.width);
rec.size.height=str2double(x.size.height);
rec.size.depth=str2double(x.size.depth);

rec.segmented=strcmp(x.segmented,'1');

rec.imgname=[x.folder '/JPEGImages/' x.filename];
rec.imgsize=str2double({x.size.width x.size.height x.size.depth});
rec.database=rec.source.database;

for i=1:length(x.object)
    rec.objects(i)=xmlobjtopas(x.object(i));
end

function p = xmlobjtopas(o)

p.class=o.name;

if isfield(o,'pose')
    if strcmp(o.pose,'Unspecified')
        p.view='';
    else
        p.view=o.pose;
    end
else
    p.view='';
end

if isfield(o,'truncated')
    p.truncated=strcmp(o.truncated,'1');
else
    p.truncated=false;
end

if isfield(o,'occluded')
    p.occluded=strcmp(o.occluded,'1');
else
    p.occluded=false;
end

if isfield(o,'difficult')
    p.difficult=strcmp(o.difficult,'1');
else
    p.difficult=false;
end

p.label=['PAS' p.class p.view];
if p.truncated
    p.label=[p.label 'Trunc'];
end
if p.occluded
    p.label=[p.label 'Occ'];
end
if p.difficult
    p.label=[p.label 'Diff'];
end

p.orglabel=p.label;

p.bbox=str2double({o.bndbox.xmin o.bndbox.ymin o.bndbox.xmax o.bndbox.ymax});

p.bndbox.xmin=str2double(o.bndbox.xmin);
p.bndbox.ymin=str2double(o.bndbox.ymin);
p.bndbox.xmax=str2double(o.bndbox.xmax);
p.bndbox.ymax=str2double(o.bndbox.ymax);

if isfield(o,'polygon')
    %warning('polygon unimplemented');
    p.polygon=[];
else
    p.polygon=[];
end

if isfield(o,'mask')
    warning('mask unimplemented');
    p.mask=[];
else
    p.mask=[];
end

if isfield(o,'part')&&~isempty(o.part)
    p.hasparts=true;
    for i=1:length(o.part)
        p.part(i)=xmlobjtopas(o.part(i));
    end
else
    p.hasparts=false;
    p.part=[];
end

if isfield(o,'actions')
    p.hasactions=true;
    fn=fieldnames(o.actions);
    for i=1:numel(fn)
        p.actions.(fn{i})=strcmp(o.actions.(fn{i}),'1');
    end
else
    p.hasactions=false;
    p.actions=[];
end

author: quantombone

repo: exemplarsvm


text:
function [x, options, flog, pointlog, scalelog] = scg(f, x, options, gradf, varargin)
%SCG	Scaled conjugate gradient optimization.
%
%	Description
%	[X, OPTIONS] = SCG(F, X, OPTIONS, GRADF) uses a scaled conjugate
%	gradients algorithm to find a local minimum of the function F(X)
%	whose gradient is given by GRADF(X).  Here X is a row vector and F
%	returns a scalar value. The point at which F has a local minimum is
%	returned as X.  The function value at that point is returned in
%	OPTIONS(8).
%
%	[X, OPTIONS, FLOG, POINTLOG, SCALELOG] = SCG(F, X, OPTIONS, GRADF)
%	also returns (optionally) a log of the function values after each
%	cycle in FLOG, a log of the points visited in POINTLOG, and a log of
%	the scale values in the algorithm in SCALELOG.
%
%	SCG(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional arguments to
%	be passed to F() and GRADF().     The optional parameters have the
%	following interpretations.
%
%	OPTIONS(1) is set to 1 to display error values; also logs error
%	values in the return argument ERRLOG, and the points visited in the
%	return argument POINTSLOG.  If OPTIONS(1) is set to 0, then only
%	warning messages are displayed.  If OPTIONS(1) is -1, then nothing is
%	displayed.
%
%	OPTIONS(2) is a measure of the absolute precision required for the
%	value of X at the solution.  If the absolute difference between the
%	values of X between two successive steps is less than OPTIONS(2),
%	then this condition is satisfied.
%
%	OPTIONS(3) is a measure of the precision required of the objective
%	function at the solution.  If the absolute difference between the
%	objective function values between two successive steps is less than
%	OPTIONS(3), then this condition is satisfied. Both this and the
%	previous condition must be satisfied for termination.
%
%	OPTIONS(9) is set to 1 to check the user defined gradient function.
%
%	OPTIONS(10) returns the total number of function evaluations
%	(including those in any line searches).
%
%	OPTIONS(11) returns the total number of gradient evaluations.
%
%	OPTIONS(14) is the maximum number of iterations; default 100.
%
%	See also
%	CONJGRAD, QUASINEW
%

%	Copyright (c) Ian T Nabney (1996-2001)

%  Set up the options.
if length(options) < 18
  error('Options vector too short')
end

if(options(14))
  niters = options(14);
else
  niters = 100;
end

display = options(1);
gradcheck = options(9);

% Set up strings for evaluating function and gradient
f = fcnchk(f, length(varargin));
gradf = fcnchk(gradf, length(varargin));

nparams = length(x);

%  Check gradients
if (gradcheck)
  feval('gradchek', x, f, gradf, varargin{:});
end

sigma0 = 1.0e-4;
fold = feval(f, x, varargin{:});	% Initial function value.
fnow = fold;
options(10) = options(10) + 1;		% Increment function evaluation counter.
gradnew = feval(gradf, x, varargin{:});	% Initial gradient.
gradold = gradnew;
options(11) = options(11) + 1;		% Increment gradient evaluation counter.
d = -gradnew;				% Initial search direction.
success = 1;				% Force calculation of directional derivs.
nsuccess = 0;				% nsuccess counts number of successes.
beta = 1.0;				% Initial scale parameter.
betamin = 1.0e-15; 			% Lower bound on scale.
betamax = 1.0e100;			% Upper bound on scale.
j = 1;					% j counts number of iterations.
if nargout >= 3
  flog(j, :) = fold;
  if nargout == 4
    pointlog(j, :) = x;
  end
end

% Main optimization loop.
while (j <= niters)

  % Calculate first and second directional derivatives.
  if (success == 1)
    mu = d*gradnew';
    if (mu >= 0)
      d = - gradnew;
      mu = d*gradnew';
    end
    kappa = d*d';
    if kappa < eps
      options(8) = fnow;
      return
    end
    sigma = sigma0/sqrt(kappa);
    xplus = x + sigma*d;
    gplus = feval(gradf, xplus, varargin{:});
    options(11) = options(11) + 1;
    theta = (d*(gplus' - gradnew'))/sigma;
  end

  % Increase effective curvature and evaluate step size alpha.
  delta = theta + beta*kappa;
  if (delta <= 0)
    delta = beta*kappa;
    beta = beta - theta/kappa;
  end
  alpha = - mu/delta;

  % Calculate the comparison ratio.
  xnew = x + alpha*d;
  fnew = feval(f, xnew, varargin{:});
  options(10) = options(10) + 1;
  Delta = 2*(fnew - fold)/(alpha*mu);
  if (Delta  >= 0)
    success = 1;
    nsuccess = nsuccess + 1;
    x = xnew;
    fnow = fnew;
  else
    success = 0;
    fnow = fold;
  end

  if nargout >= 3
    % Store relevant variables
    flog(j) = fnow;		% Current function value
    if nargout >= 4
      pointlog(j,:) = x;	% Current position
      if nargout >= 5
	scalelog(j) = beta;	% Current scale parameter
      end
    end
  end
  if display > 0
    fprintf(1, 'Cycle %4d  Error %11.6f  Scale %e\n', j, fnow, beta);
  end

  if (success == 1)
    % Test for termination

    if (max(abs(alpha*d)) < options(2) & max(abs(fnew-fold)) < options(3))
      options(8) = fnew;
      return;

    else
      % Update variables for new position
      fold = fnew;
      gradold = gradnew;
      gradnew = feval(gradf, x, varargin{:});
      options(11) = options(11) + 1;
      % If the gradient is zero then we are done.
      if (gradnew*gradnew' == 0)
	options(8) = fnew;
	return;
      end
    end
  end

  % Adjust beta according to comparison ratio.
  if (Delta < 0.25)
    beta = min(4.0*beta, betamax);
  end
  if (Delta > 0.75)
    beta = max(0.5*beta, betamin);
  end

  % Update search direction using Polak-Ribiere formula, or re-start
  % in direction of negative gradient after nparams steps.
  if (nsuccess == nparams)
    d = -gradnew;
    nsuccess = 0;
  else
    if (success == 1)
      gamma = (gradold - gradnew)*gradnew'/(mu);
      d = gamma*d - gradnew;
    end
  end
  j = j + 1;
end

% If we get here, then we haven't terminated in the given number of
% iterations.

options(8) = fold;
if (options(1) >= 0)
  disp(maxitmess);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function test_ft_freqdescriptives

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqdescriptives
% DATATYPE freq

freq.label = {'1';'2'};
freq.freq  = 1:10;
freq.time  = 1:5;
freq.dimord = 'rpt_chan_freq_time';
freq.powspctrm = randn(100,2,10,5);

cfg = [];
cfg.variance = 'yes';
freqout = ft_freqdescriptives(cfg, freq);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% imStats(IM1,IM2)
%
% Report image (matrix) statistics.
% When called on a single image IM1, report min, max, mean, stdev, skew,
% and kurtosis (4th moment about the mean, divided by squared variance)
%
% When called on two images (IM1 and IM2), report min, max, mean,
% stdev of the difference, and also SNR (relative to IM1).

% Eero Simoncelli, 6/96.

function [] = imStats(im1,im2)

if (~isreal(im1))
  error('Args must be real-valued matrices');
end

if (exist('im2') == 1)
  difference = im1 - im2;
  [mn,mx] = range2(difference);
  mean = mean2(difference);
  v = var2(difference,mean);
  if (v < realmin)
    snr = Inf;
  else
    snr = 10 * log10(var2(im1)/v);
  end
  fprintf(1, 'Difference statistics:\n');
  fprintf(1, '  Range: [%c, %c]\n',mn,mx);
  fprintf(1, '  Mean: %f,  Stdev (rmse): %f,  SNR (dB): %f\n',...
      mean,sqrt(v),snr);
else
  [mn,mx] = range2(im1);
  mean = mean2(im1);
  var = var2(im1,mean);
  stdev = sqrt(real(var))+sqrt(imag(var));
  sk = skew2(im1, mean, stdev^2);
  kurt = kurt2(im1, mean, stdev^2);
  fprintf(1, 'Image statistics:\n');
  fprintf(1, '  Range: [%f, %f]\n',mn,mx);
  fprintf(1, '  Mean: %f,  Stdev: %f,  Skew: %f,  Kurt: %f\n',mean,stdev,sk,kurt);
end

repo: SelfExSR


text:
%% Classic Monte Carlo simualtion
% Vincent Leclercq, The MathWorks, 2007, vincent.leclercq@mathworks.fr
%

clear all;
close all;

NbTrials = 10000;

%RunMode = 'LogNomrality';

RunMode = 'OptionPricing';
%% Load Data (retrieved originally from Thomson Datastream)

load Equities.mat
PastDate = today()-1 * 365;

%% Retrieve the dates in a numeric format

Dates   = cellfun(@(x)(datenum(x,'yyyy-mm-ddTHH:MM:SS')),Equities.DATE);
AssetPrices = cellfun(@(x) (str2double(x)),Equities.P );

%% Plot the Series

plot(Dates,AssetPrices);set(gcf,'WindowStyle','Docked');
legend(Equities.DISPNAME);
datetick('x','mmmyy');
xlim([PastDate, today()]);
grid on;

%% Compute the returns
% We can compute the returns for one or many stocks at the same time using
% matrix computation and matlab easy syntax

Suez_Returns = tick2ret(AssetPrices);
DailyVol = std(Suez_Returns);
AnnualVol = DailyVol * sqrt(252);

SpotPrice = AssetPrices(end);

InterestRate = 0.0375;

%% Portfolio simulation (Monte Carlo) using the Financial toolbox function
% For help, one can use the doc portsim function
% We call the portfolio simulation using Financial toolbox to
% simulate 10000 scenarios. Of course, correlation are preserved
% We assume an horizon of 6 * 22 trading days, ie 6 month maturity

%% Using Annual statistics

SimulatedRetsAnnual = portsim(InterestRate,AnnualVol^2, 12, 1/12, NbTrials,'Expected'); % NbStep * TimeStep = 1 (in years !!!)

%% Using Daily statistics
NumberOfSimulationSteps = 12;
SimulatedRetsDaily = portsim(InterestRate./252, DailyVol^2, 12, 252./12, NbTrials,'Expected');% NbStep * TimeStep = 252 (in days!!!)

%% Generate the Prices and plot them

SimulatedPricesAnnual = ret2tick(squeeze(SimulatedRetsAnnual) ,SpotPrice);
SimulatedPricesDaily  = ret2tick(squeeze(SimulatedRetsDaily)  ,SpotPrice );

figure;hist(SimulatedPricesAnnual(end,:),40);title('Prices, annual timestep used');set(gcf,'WindowStyle','Docked');
figure;hist(SimulatedPricesDaily(end,:),40);title('Prices, daily timestep used');set(gcf,'WindowStyle','Docked');

%% Check For LogNormality of the Price series

ExpectedVariance = (SpotPrice^2)  * (exp(AnnualVol^2) - 1)* exp(2*InterestRate);

disp(['Mean Price (Annual Parameters) -> ', num2str(mean(SimulatedPricesAnnual(end,:))) ' , Theoric value (Hull) :' num2str(SpotPrice * exp(InterestRate))]);
disp(['Mean Price (Daily Parameters) -> ', num2str(mean(SimulatedPricesDaily(end,:))) ' , Theoric value (Hull) :' num2str(SpotPrice * exp(InterestRate))]);

disp(['Expected Variance (Annual Parameters) -> ', num2str(var(SimulatedPricesAnnual(end,:))) ' , Theoric value (Hull) :' num2str(ExpectedVariance)]);
disp(['Expected Variance (Daily Parameters) -> ', num2str(var(SimulatedPricesDaily(end,:))) ' , Theoric value (Hull) :' num2str(ExpectedVariance)]);

%% Parameter sweep
% Now that we have done this, we can ccompute the same thing for different
% Exercise prices
if strcmp(RunMode,'OptionPricing')
    k = 1;
    NumberOfSteps = 400;
    ExercisePrices= linspace(0.8 * SpotPrice,1.2 * SpotPrice,NumberOfSteps);

    VanillaPriceAnnual    = zeros(NumberOfSteps,1);
    VanillaPriceDaily    = zeros(NumberOfSteps,1);

    ProbabilityITMAnnual  = zeros(NumberOfSteps,1);
    ProbabilityITMDaily  = zeros(NumberOfSteps,1);
    CIAnnual    =     zeros(NumberOfSteps,2);
    CIDaily     =     zeros(NumberOfSteps,2);
    BLSPrices       = zeros(NumberOfSteps,1);

    %%
    TimeInYear = 1;
    for i = 1 : NumberOfSteps
        [BLSPrices(i),dummy]   = blsprice(SpotPrice, ExercisePrices(i), InterestRate, TimeInYear, AnnualVol, 0);
        [VanillaPriceAnnual(i), ProbabilityITMAnnual(i) ,CIAnnual(i,:)] = GetOptionPrice(SimulatedPricesAnnual,ExercisePrices(i),TimeInYear,InterestRate,'Vanilla');
        [VanillaPriceDaily(i), ProbabilityITMDaily(i) ,  CIDaily(i,:)] = GetOptionPrice(SimulatedPricesDaily,ExercisePrices(i),TimeInYear,InterestRate,'Vanilla');

    end;

%%

    h =    figure;
    [AX,H1,H2] = plotyy(ExercisePrices,[VanillaPriceAnnual BLSPrices CIAnnual] , ExercisePrices,ProbabilityITMAnnual);

    xlabel('Exercise Price');
    title('Option prices for a Vanilla option using a 1 year - Annual volatility');
    Axes_YLabels = get(AX,'Ylabel');
    set(Axes_YLabels{1},'String','Option Price') ;

    set(H1(1),'LineStyle','-');
    set(H1(1),'Color','r');
    set(H1(1),'LineWidth',2);

    set(H1(2),'LineStyle','-');
    set(H1(2),'Color','b');
    set(H1(2),'LineWidth',2);

    set(H1(3),'LineStyle','--');
    set(H1(3),'Color','r');
    set(H1(3),'LineWidth',1);

    set(H1(4),'LineStyle',':');
    set(H1(4),'Color','r');
    set(H1(4),'LineWidth',1);

    set(H2,'LineStyle','-');
    set(H2,'Color','g');
    set(H2,'LineWidth',2);

    set(get(AX(2),'Ylabel'),'String','Probability of being In the Money');
    legend(H1,{['Option Price (Monte Carlo)'], ['Option Price (Black Scholes)'], ['99% Confidence interval (Lower)'],['99% Confidence interval (Upper)']},'Location','NorthEast');
    legend(H1(2),{'Option Price (Black Shcoles)'},'Location','NorthEast');
    legend(H2,{'Probability'},'Location','SouthWest');
    grid on;
set(h,'WindowStyle','Docked');

%%
h=     figure;
    [AX,H1,H2] = plotyy(ExercisePrices,[VanillaPriceDaily BLSPrices CIDaily] , ExercisePrices,ProbabilityITMDaily);

    xlabel('Exercise Price');
    title('Option prices for a Vanilla option using a 1 year - Daily volatility');
    Axes_YLabels = get(AX,'Ylabel');
    set(Axes_YLabels{1},'String','Option Price') ;

    set(H1(1),'LineStyle','-');
    set(H1(1),'Color','r');
    set(H1(1),'LineWidth',2);

    set(H1(2),'LineStyle','-');
    set(H1(2),'Color','b');
    set(H1(2),'LineWidth',2);

    set(H1(3),'LineStyle','--');
    set(H1(3),'Color','r');
    set(H1(3),'LineWidth',1);

    set(H1(4),'LineStyle',':');
    set(H1(4),'Color','r');
    set(H1(4),'LineWidth',1);

    set(H2,'LineStyle','-');
    set(H2,'Color','g');
    set(H2,'LineWidth',2);

    set(get(AX(2),'Ylabel'),'String','Probability of being In the Money');
    legend(H1,{['Option Price (Monte Carlo)'], ['Option Price (Black Scholes)'], ['99% Confidence interval (Lower)'],['99% Confidence interval (Upper)']},'Location','NorthEast');
    legend(H2,{'Probability'},'Location','SouthWest');
    grid on;
    set(h,'WindowStyle','Docked');
 end;


text:
classdef MOCell < ALGORITHM
% <multi> <real/integer/label/binary/permutation> <constrained/none>
% Cellular genetic algorithm

%------------------------------- Reference --------------------------------
% A. J. Nebro, J. J. Durillo, F. Luna, B. Dorronsoro, and E. Alba, MOCell:
% A cellular genetic algorithm for multiobjective optimization,
% International Journal of Intelligent Systems, 2009, 24(7): 726-746.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Neighborhood of each cell
            cSize    = floor(sqrt(Problem.N));
            n1       = reshape(1:cSize^2,cSize,cSize);
            n2       = n1([end,1:end-1],[end,1:end-1]);
            n3       = n1([end,1:end-1],:);
            n4       = n1([end,1:end-1],[2:end,1]);
            n5       = n1(:,[end,1:end-1]);
            n6       = n1(:,[2:end,1]);
            n7       = n1([2:end,1],[end,1:end-1]);
            n8       = n1([2:end,1],:);
            n9       = n1([2:end,1],[2:end,1]);
            Neighbor = [n1(:),n2(:),n3(:),n4(:),n5(:),n6(:),n7(:),n8(:),n9(:)];

            %% Generate random population
            Population = Problem.Initialization(size(Neighbor,1));
            Archive    = [];

            %% Optimization
            while Algorithm.NotTerminated(Archive)
                % Update the fitness of solutions in the population
                FrontNo  = NDSort(Population.objs,Population.cons,inf);
                CrowdDis = CrowdingDistance(Population.objs,FrontNo);

                % Generate one offspring for each cell
                Offspring(1:length(Population)) = SOLUTION();
                for i = 1 : length(Population)
                    parents      = Neighbor(i,TournamentSelection(2,2,FrontNo(Neighbor(i,:)),-CrowdDis(Neighbor(i,:))));
                    Offspring(i) = OperatorGAhalf(Problem,Population(parents));
                end

                % Replace the old population
                ViolationP = sum(max(0,Population.cons),2);
                ViolationO = sum(max(0,Offspring.cons),2);
                Replace    = ViolationO < ViolationP | ViolationO == ViolationP & any(Population.objs>Offspring.objs,2);
                Population(Replace) = Offspring(Replace);

                % Update the archive
                Archive = [Archive,Offspring];
                Archive = Archive(NDSort(Archive.objs,Archive.cons,1)==1);
                if length(Archive) > Problem.N
                    [~,rank] = sort(CrowdingDistance(Archive.objs,ones(1,length(Archive))),'descend');
                    Archive  = Archive(rank(1:Problem.N));
                end

                % Feedback
                nRep = min(min(20,length(Population)),length(Archive));
                Population(randperm(length(Population),nRep)) = Archive(randperm(length(Archive),nRep));
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function stroud_test46 ( )

%*****************************************************************************80
%
%% TEST46 tests TORUS_5S2, TORUS_6S2, TORUS_14S.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_3D_INDEX;

  r1 = 0.5;
  r2 = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST46\n' );
  fprintf ( 1, '  For the interior of a torus,\n' );
  fprintf ( 1, '  TORUS_5S2,\n' );
  fprintf ( 1, '  TORUS_6S2, and\n' );
  fprintf ( 1, '  TORUS_5S2 approximate integrals.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Inner radius = %f\n', r1 );
  fprintf ( 1, '  Outer radius = %f\n', r2 );
  fprintf ( 1, '  Volume = %f\n', torus_volume_3d ( r1, r2 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Rule:        #5S2          #6S2          #14S\n' );
  fprintf ( 1, '    F(X)\n' );
  fprintf ( 1, '\n' );

  num = function_3d_num ( );

  for i = 1 : num

    FUNC_3D_INDEX = i;

    result1 = torus_5s2 ( 'function_3d', r1, r2 );
    result2 = torus_6s2 ( 'function_3d', r1, r2 );
    result3 = torus_14s ( 'function_3d', r1, r2 );

    fname = function_3d_name ( i );

    fprintf ( 1, '  %s  %12f  %12f  %12f\n', fname, result1, result2, result3 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = perform_wavortho_transf(f,Jmin,dir,options)

% perform_wavortho_transf - compute orthogonal wavelet transform
%
%   fw = perform_wavortho_transf(f,Jmin,dir,options);
%
%   You can give the filter in options.h.
%
%   Works in arbitrary dimension.
%
%   Copyright (c) 2009 Gabriel Peyre

options.null = 0;
h = getoptions(options,'h', compute_wavelet_filter('Daubechies',4) );
g = [0 h(length(h):-1:2)] .* (-1).^(1:length(h));

n = size(f,1);
Jmax = log2(n)-1;

if dir==1
    %%% FORWARD %%%
    for j=Jmax:-1:Jmin
        sel = 1:2^(j+1);
        a = subselect(f,sel);
        for d=1:nb_dims(f)
            a = cat(d, subsampling(cconv(a,h,d),d), subsampling(cconv(a,g,d),d) );
        end
        f = subassign(f,sel,a);
    end
else
    %%% FORWARD %%%
    for j=Jmin:Jmax
        sel = 1:2^(j+1);
        a = subselect(f,sel);
        for d=1:nb_dims(f)
            w = subselectdim(a,2^j+1:2^(j+1),d);
            a = subselectdim(a,1:2^j,d);
            a = cconv(upsampling(a,d),reverse(h),d) + cconv(upsampling(w,d),reverse(g),d);
        end
        f = subassign(f,sel,a);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = subselect(f,sel)
switch nb_dims(f)
    case 1
        f = f(sel);
    case 2
        f = f(sel,sel);
    case 3
        f = f(sel,sel,sel);
    case 4
        f = f(sel,sel,sel,sel);
    case 5
        f = f(sel,sel,sel,sel,sel);
    case 6
        f = f(sel,sel,sel,sel,sel,sel);
    case 7
        f = f(sel,sel,sel,sel,sel,sel,sel);
    case 8
        f = f(sel,sel,sel,sel,sel,sel,sel,sel);
    otherwise
        error('Not implemented');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = subselectdim(f,sel,d)
switch d
    case 1
        f = f(sel,:,:,:,:,:,:,:);
    case 2
        f = f(:,sel,:,:,:,:,:,:);
    case 3
        f = f(:,:,sel,:,:,:,:,:);
    case 4
        f = f(:,:,:,sel,:,:,:,:);
    case 5
        f = f(:,:,:,:,sel,:,:,:);
    case 6
        f = f(:,:,:,:,:,sel,:,:);
    case 7
        f = f(:,:,:,:,:,:,sel,:);
    case 8
        f = f(:,:,:,:,:,:,:,sel);
    otherwise
        error('Not implemented');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = subassign(f,sel,g)
switch nb_dims(f)
    case 1
        f(sel) = g;
    case 2
        f(sel,sel) = g;
    case 3
        f(sel,sel,sel) = g;
    case 4
        f(sel,sel,sel,sel) = g;
    case 5
        f(sel,sel,sel,sel,sel) = g;
    case 6
        f(sel,sel,sel,sel,sel,sel) = g;
    case 7
        f(sel,sel,sel,sel,sel,sel,sel) = g;
    case 8
        f(sel,sel,sel,sel,sel,sel,sel,sel) = g;
    otherwise
        error('Not implemented');
end

author: gpeyre

repo: matlab-toolboxes


text:
function [b0, dt6, mmPerVox] = dtiComputeTensorFromRaw(firstDataFile, gradsFile)
% Computes diffusion tensor from raw diffusion-weighted images.
%
% [b0, dt6, mmPerVox] = dtiComputeTensorFromRaw(firstDataFile, gradsFile)
%
% Loads a whole directory of raw image files, combines them according to
% the gradient direction information in gradsFile, then computes the
% diffusion tensor. Also returns the xform that will rotate the images into
% our cannonical orientation (as close to axial as possible- see
% computeCannonicalXformFromIfile).
%
% TODO:
%
%
% HISTORY:
% 2005.05.25: RFD (bob@sirl.stanford.edu) wrote it.
%

if(~exist('gradsFile','var')), gradsFile = []; end
if(~exist('firstDataFile','var')), firstDataFile = []; end

% Load raw data, averaging repeats
%
% *** TODO: we should NOT average repeated b>0 images. Rather, we should
% fit the tensor to all the individual data points.
[img, mmPerVox, gradDirs, bVals, xformToCan] = dtiAverageRawData(firstDataFile, gradsFile);

% Compute the tensor
%
% Maybe try code in AFNI's 3dDWItoDT?
% http://afni.nimh.nih.gov/sscc/presentations/Diffusion%20Tensor.pdf

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function fid = logfileinit(logfile)
% fid = logfileinit(logfile)
fid = fopen(logfile,'wt');
fprintf(fid,'Computation: %s\n\n', datestr(now));

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: October 8th, 2018
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs / non-invariant beams *)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: helps input previous .vtk information for restarting a
%           simulation that has ended because of power failure, etc.
%
%      NOTE: for restart protocol, need to have .vtk data for:
%                       1. lagPts (Lagrangian positions)
%                       2. u (velocity field)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [current_time,cter,ctsave,U,V,xLag,yLag,xLag_P,yLag_P,path_to_data] = help_Me_Restart(dt)

%
% NEEDS TO BE HARDCODED PER SIMULATION BASIS
%
ctsave = 3;               % Last time-step of data saved (# at end of .vtk file);
print_dump = 40;           % Print_dump interval as given in input2d

% Path to Simulation Data (e.g., including viz_IB2d folder)
path_to_data = '/Users/battistn/Desktop/IB2d/matIB2d/Examples/Example_Test_Restart_Protocol/Standard_Rubberband_Restart/viz_IB2d/';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ------------- DO NOT CHANGE BELOW --------------- %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[current_time,cter,ctsave,U,V,xLag,yLag,xLag_P,yLag_P] = pass_Back_Data_For_Restart(dt,ctsave,print_dump,path_to_data);

author: nickabattista


text:
function K = besselj_correlation ( s, t )

%*****************************************************************************80
%
%% BESSELJ_CORRELATION evaluates the Bessel J correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S(*), T(*), pairs of argument values.
%
%    Output, real K(*), the correlation function values
%
  K = besselj ( 0, abs ( s - t ) );

  return
end


text:
function intVal=monomialIntSpherSurf(alpha)
%%MONOMIALINTSPHERSURF Evaluate the integral of prod_{i=1}^Nx(i)^alpha(i)
%          over the surface defined by sum(x.^2)=1. That is, this function
%          gives the value of the desired monomial integral taken over the
%          surface of an n-dimensional hypersphere. The region is
%          designated as U_n in [1].
%
%INPUTS: alpha An NX1 or 1XN vector of the integer exponents of the
%              monomial term. All elements must be >=0.
%
%OUTPUTS: intVal The value of the specified integral.
%
%The formula is taken from Chapter 7.6 of [1]. These types of explicit
%moment formulae are useful for testing cubature integration points.
%
%EXAMPLE:
%Here we verify that the moment value produced here equals that obtained
%using cubature points of an appropriate order.
% [xi,w]=seventhOrderSpherSurfCubPoints(3);
% alpha=[0;4;2];
% theMoment=findMomentFromSamp(alpha,xi,w)
% intVal=monomialIntSpherSurf(alpha)
%One will see that theMoment and intVal are both about 0.3590.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(all(mod(alpha,2)==0))
    n=length(alpha);

    intVal=2*exp(sum(gammaln((alpha+1)/2))-gammaln((n+sum(alpha))/2));
else
    intVal=0;
end
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:
function [hFig, iDS, iFig] = view_channels_3d(FileNames, Modality, SurfaceType, is3DElectrodes, isDetails)
% VIEW_CHANNELS_3D: Display channel files on top of subject anatomy.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2019

global GlobalData;
% Parse inputs
if (nargin < 5) || isempty(isDetails)
    isDetails = 0;
end
if (nargin < 4) || isempty(is3DElectrodes)
    is3DElectrodes = 0;
end
if (nargin < 3) || isempty(SurfaceType)
    SurfaceType = 'scalp';
end
if ischar(FileNames)
    FileNames = {FileNames};
end
hFig = [];
iDS = [];
iFig = [];
% Coils or channel markers?
isShowCoils = ismember(Modality, {'Vectorview306', 'CTF', '4D', 'KIT', 'KRISS', 'BabyMEG', 'NIRS-BRS', 'RICOH'});

% === DISPLAY SURFACE ===
% Get study
[sStudy, iStudy] = bst_get('ChannelFile', FileNames{1});
if isempty(sStudy)
    return
end
% Get subject
sSubject = bst_get('Subject', sStudy.BrainStormSubject);
% View surface if available
if ~isempty(sSubject)
    % If displaying MEG coils: remove completely the transparency for nicer display with Matlab >= 2014b
    if isShowCoils
        opaqueAlpha = 0;
    else
        opaqueAlpha = .1;
    end
    % If passing a filename
    if ~isempty(strfind(SurfaceType, '.mat'))
        SurfaceFile = SurfaceType;
        SurfaceType = file_gettype(SurfaceFile);
    else
        SurfaceFile = [];
    end
    % Display surface
    switch lower(SurfaceType)
        case 'cortex'
            if ~isempty(sSubject.iCortex) && (sSubject.iCortex <= length(sSubject.Surface))
                if isempty(SurfaceFile)
                    SurfaceFile = sSubject.Surface(sSubject.iCortex).FileName;
                end
                switch (Modality)
                    case 'SEEG',  SurfAlpha = .8;
                    case 'ECOG',  SurfAlpha = .2;
                    otherwise,    SurfAlpha = opaqueAlpha;
                end
                hFig = view_surface(SurfaceFile, SurfAlpha, [], 'NewFigure');
            end
        case 'innerskull'
            if ~isempty(sSubject.iInnerSkull) && (sSubject.iInnerSkull <= length(sSubject.Surface))
                if isempty(SurfaceFile)
                    SurfaceFile = sSubject.Surface(sSubject.iInnerSkull).FileName;
                end
                switch (Modality)
                    case 'SEEG',  SurfAlpha = .5;
                    case 'ECOG',  SurfAlpha = .2;
                    otherwise,    SurfAlpha = opaqueAlpha;
                end
                hFig = view_surface(SurfaceFile, SurfAlpha, [], 'NewFigure');
            end
        case 'scalp'
            if ~isempty(sSubject.iScalp) && (sSubject.iScalp <= length(sSubject.Surface))
                if isempty(SurfaceFile)
                    SurfaceFile = sSubject.Surface(sSubject.iScalp).FileName;
                end
                switch (Modality)
                    case 'SEEG',  SurfAlpha = .8;
                    case 'ECOG',  SurfAlpha = .8;
                    case 'NIRS-BRS'
                        if isDetails
                            SurfAlpha = .8;
                        else
                            SurfAlpha = opaqueAlpha;
                        end
                    otherwise
                        SurfAlpha = opaqueAlpha;
                end
                hFig = view_surface(SurfaceFile, SurfAlpha, [], 'NewFigure');
            end
        case {'anatomy', 'subjectimage'}
            if ~isempty(sSubject.iAnatomy) && (sSubject.iAnatomy <= length(sSubject.Anatomy))
                if isempty(SurfaceFile)
                    SurfaceFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
                end
                SurfAlpha = .1;
                hFig = view_mri_3d(SurfaceFile, [], SurfAlpha, 'NewFigure');
            end
    end
end
% Warning if no surface was found
if isempty(hFig)
    disp('BST> Warning: The anatomy of this subject was not imported properly.');
end

% === DISPLAY CHANNEL FILES ===
% Only one: Markers and labels
if (length(FileNames) == 1)
    if is3DElectrodes
        isLabels = 0;
    else
        isLabels = 1;
    end
    isMarkers = ~isShowCoils || isDetails;
    [hFig, iDS, iFig] = view_channels(FileNames{1}, Modality, isMarkers, isLabels, hFig, is3DElectrodes);
    % SEEG and ECOG: Open tab "iEEG"
    if ismember(Modality, {'SEEG', 'ECOG', 'ECOG+SEEG'})
        gui_brainstorm('ShowToolTab', 'iEEG');
    end
% Multiple: Markers only
else
    ColorTable = [1,0,0; 0,1,0; 0,0,1];
    for i = 1:length(FileNames)
        color = ColorTable(mod(i-1,size(ColorTable,1))+1,:);
        % View channel file
        [hFig, iDS, iFig] = view_channels(FileNames{i}, Modality, 1, 0, hFig);
        % Rename objects so that they are not deleted by the following call
        hPatch = findobj(hFig, 'tag', 'SensorsPatch');
        set(hPatch, 'FaceColor', color, 'FaceAlpha', .3, 'SpecularStrength', 0, ...
                    'EdgeColor', color, 'EdgeAlpha', .2, 'LineWidth', 1, ...
                    'Marker', 'none', 'Tag', 'MultipleSensorsPatches');
        if (i ~= length(FileNames))
            hPatch = copyobj(hPatch, get(hPatch,'Parent'));
            % Delete loaded channel information to force it to be reloaded by the view_channels function
            GlobalData.DataSet(iDS).Channel = [];
        end
    end
end

author: brainstorm-tools


text:
function X = dissect(Y)

X = dissect_internal(Y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains2(tpl)
% create_restrictions_and_markov_chains2 -- creates restrictions and
% markov chains for the SVAR model in which Coefficients and variances have
% different chains, different regimes, and different durations
%
% ::
%
%
%   [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains2(tpl)
%
% Args:
%
%    - **tpl** [struct]: template created for SVAR objects
%
% Returns:
%    :
%
%    - **lin_restr** [cell]: two column-cell (see below). The first column
%    contains COEF objects or linear combinations of COEF objects, which are
%    themselves COEF objects.
%
%    - **nonlin_restr** [cell]: one column-cell containing inequality or
%    nonlinear restrictions
%
%    - **tpl** [struct]: modified template
%
% Note:
%
%    - The syntax to construct an advanced COEF object is
%    a=coef(eqtn,vname,lag,chain_name,state)
%      - **eqtn** [integer|char]: integer or variable name
%      - **vname** [integer|char]: integer or variable name
%      - **lag** [integer]: integer or variable name
%      - **chain_name** [char]: name of the markov chain
%      - **state** [integer]: state number
%
%    - RISE sorts the endogenous variables alphabetically and use this order
%    to tag each equation in SVAR and RFVAR models.
%
%    - The lag coefficients are labelled a0, a1, a2,...,ak, for a model with k
%    lags. Obviously, a0 denotes the contemporaneous coefficients.
%
%    - The constant terms labelled c_1_1, c_2_2,...,c_n_n, for a model with n
%    endogenous variables.
%
%    - The standard deviations labelled sig_1_1, sig_2_2,...,sig_n_n, for a
%    model with n endogenous variables.
%
% Example:
%    coef('pi','ygap',0,'policy',1)
%
%    coef(2,3,0,'policy',1)
%
%    coef(2,'ygap',0,'policy',1)
%
%    coef('pi',3,0,'policy',1)
%
%    See also:

% We borrow both the restrictions and the markov chains from the model in
% which all coefficients across all equations switch in lockstep.
%--------------------------------------------------------------------------
[lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains1(tpl);

% Then we add another chain controling all variances across all equations
%-------------------------------------------------------------------------
% We just make sure we do not change the restrictions
last=numel(tpl.markov_chains);
tpl.markov_chains(last+1)=struct('name','syncvol',...
    'states_expected_duration',[2+1i,2+1i,2+1i],...
    'controlled_parameters',{{'sig'}});

end


text:
function [latout,lonout] = my_interpm(lat,lon,maxdiff)
% [latout,lonout] = my_interpm(lat,lon,maxdiff)
% This function  fills in any gaps in latitude (lat) or longitude (lon) data vectors
% that are greater than a defined tolerance maxdiff apart in either dimension.
% lat and lon should be row vectors.
%
% latout and lonout are the new latitude and longitude data vectors, in which any gaps
% larger than maxdiff in the original vectors have been filled with additional points.
%
% by Jindong Wang on March 1, 2018
%
if isrow(lat)
    lat=lat';
end
if isrow(lon)
    lon=lon';
end
ny=size(lat,1);
nx=size(lon,1);
if nx~=ny || nx<2 || ny<2 || maxdiff<10^-8
    error('Error: Wrong input!');
end

dlat=abs(lat(2:end)-lat(1:end-1));
dlon=abs(lon(2:end)-lon(1:end-1));
nin=ceil(max([dlat,dlon],[],2)/maxdiff)-1;
sumnin=sum(nin,'omitnan');
if sumnin==0
    disp('No incertion needed.');
    latout=lat;
    lonout=lon;
    return
end
nout=sumnin+nx;
latout=nan(nout,1);
lonout=nan(nout,1);

n=1;
for i=1:nx-1;
    ni=nin(i);
    if ni==0 || isnan(ni)
        latout(n)=lat(i);
        lonout(n)=lon(i);
        nstep=1;
    else
        ilat=linspace(lat(i),lat(i+1),ni+2);
        ilon=linspace(lon(i),lon(i+1),ni+2);
        latout(n:n+ni)=ilat(1:ni+1);
        lonout(n:n+ni)=ilon(1:ni+1);
        nstep=ni+1;
    end
    n=n+nstep;
end
latout(end)=lat(end);
lonout(end)=lon(end);

author: CHLNDDEV


text:
function [L] = ilogit(t)
% Calculate the inverse logit function corresponding to the value t
%
% :Usage:
% ::
%
%     function [L] = ilogit(t)
%
% :Output:
%
%   **L:**
%        exp(t)./(1+exp(t));
%
% ..
%    By Martin Lindquist and Tor Wager
%    Edited 12/12/06
% ..

L = exp(t)./(1+exp(t));

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% fig_taper_cos3
% compare 3 different interpolation kernels:
%	min-max kernel with uniform scaling factors
%	min-max kernel with optimized scaling factors
%	cos^3-tapered dirichlet kernel of magnusson

J = 6;
N = 1000;
K = 2 * N;
k = linspace(-J/2,J/2,101)';
fm = diric(2*pi*k/J, J) .* cos((2*pi*k/J)/2).^3;	% tapered
fn = nufft1_kernel(k, N, J, K, 1, 0);			% min-max unif
fo = nufft1_kernel(k, N, J, K, 'best');
fo = fo / max(fo);

%fm = 0*fm;
%fn = nufft1_kernel(k, N, J, 9*N, 1, 0);
%fo = nufft1_kernel(k, N, J, 2*N, 1, 0);

plot(k, fn, '-', k, fm, '--', k, fo, '-.')
axisy(-0.3,1.1)
legend('min-max', 'Magnusson', 'optimized')
title(sprintf('J=%d K/N=%g', J, K/N))


text:
function trans = homothecy(point, ratio)
%HOMOTHECY create a homothecy as an affine transform.
%
%   TRANS = homothecy(POINT, K);
%   POINT is the center of the homothecy, K is its factor.
%
%   See also:
%   transforms2d, transformPoint, createTranslation

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2005-01-20
% Copyright 2005 INRA - TPV URPOI - BIA IMASTE

% deprecation warning
warning('geom2d:deprecated', ...
    '''homothecy'' is deprecated, use ''createHomothecy'' instead');

% call current implementation
trans = createHomothecy(point, ratio);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function model = createToyModelForConstraints()
% createToyModelForConstraints creates a toy model for
% including coupling constraints for functional testing.
% The model created looks as follows:
%
%
%   <-> A -----> B ---> C --> E <->
%        \       ^      ^
%         \     /      /
%          -> D  ---> F
%
%      H -> I -> J
%       ^        /
%        \      /
%           ---
%
model = createModel();
%Reactions in {Rxn Name, MetaboliteNames, Stoichiometric coefficient} format
Reactions = {'R1',{'A','B'},[-1 1],-1000,1000;...
             'R2',{'A', 'D'},[-1 2],0,1000;...
             'R3',{'D','B'},[-1 1 ],0,30;...
             'R4',{'B','C'},[-1 1],0,20;...
             'R5',{'C','E'},[-1 1],0,1000;...
             'R6',{'D','F'},[-1 1],0,20;...
             'R7',{'F','C'},[-1 1],0,30;...
             'R8',{'H','I'},[-1 1],0,1000;...
             'R9',{'I','J'},[-1 1],0,1000;...
             'R10',{'J','H'},[-1 1],0,1000;...
             };
ExchangedMets = {'A','E'};
%Add Reactions
for i = 1:size(Reactions,1)
    %All reactions are irreversible
    model = addReaction(model,Reactions{i,1},'metaboliteList',Reactions{i,2},'stoichCoeffList',Reactions{i,3},'lowerBound',Reactions{i,4},'upperBound',Reactions{i,5});
end
model = changeGeneAssociation(model,'R3','G1 or G2');
model = changeGeneAssociation(model,'R4','G3 or G5');
model = changeGeneAssociation(model,'R7','G4 and (G1 or G2)');
model = changeGeneAssociation(model,'R4','(G3 or G5) and G6');

%Add Exchangers
model = addExchangeRxn(model,ExchangedMets,-1000*ones(numel(ExchangedMets),1),1000*ones(numel(ExchangedMets),1));
model = changeObjective(model,'EX_E',1);
model = addCOBRAConstraints(model,{'R3','R7'},5, 'c',[1,1],'dsense','G','ConstraintID','R3_and_R7_above_5');
model = addCOBRAConstraints(model,{'R4','R6'},30, 'c',[1,1],'dsense','L','ConstraintID','R4_and_R6_lowerthan_30');

author: opencobra

repo: cobratoolbox


text:
%% Estimate the quality of QRS complex detections
% Calculate the quality of QRS detections based on the the likelihood of
% measurements from RR interval series with respect to a trained model. See
% for further details. See 'A Pattern-Recognition Approach for
% Lead-Selection in Heartbeat Detection' CINC 2014
%
% Example
%
%   payload_out = calculateSeriesQuality(payload_out, ECG_header, start_end_this_segment)
%
% Arguments:
%   + payload_out: structure with the QRS detections series to estimate the
%       quality. One fieldname contains the QRS detections series in a
%       'time' field. For example
%       payload_out.('QRS_detector1_name_lead_I').time contains the
%       detections of QRS_detector1 in lead1. See ECGtask_QRS_detection for
%       examples.%
%   + ECG_header is a structure (ECG_header prop. in ECGwrapper object)
%       describing the signal.
%   + start_end_this_segment: array samples where starts and finish this
%       signal segment [start_sample end_sample].
%
% Output:
%   + payload_out: The same input structure with extra field
%       'series_quality' describing the quaility of the QRS detection
%       series in three fields: ratios, a measure of quality of the QRS
%       detections, estimated_labs, for internal use, and AnnNames which
%       matches the previous fields to the fieldname found in payload_out.
%
% See also ECGtask_QRS_detection
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Birthdate: 01/01/2012
% Last update: 18/10/2014
% Copyright 2008-2015
function payload_out = calculateSeriesQuality(payload_out, ECG_header, start_end_this_segment)

    if( isempty(payload_out) )
        return
    end

    % estimate quality of QRS detections performed
    [AnnNames, all_annotations] = getAnnNames(payload_out);

    [ ratios, estimated_labs ] = CalcRRserieRatio(all_annotations, ECG_header, start_end_this_segment);
    [ratios, best_detections_idx] = sort(ratios, 'descend');
    AnnNames = AnnNames(best_detections_idx,:);

    payload_out.series_quality.ratios = ratios;
    payload_out.series_quality.estimated_labs = estimated_labs(best_detections_idx);
    payload_out.series_quality.AnnNames = AnnNames;
    payload_out.series_quality.sampfreq = ECG_header.freq;

author: marianux

repo: ecg-kit

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


text:
function C = VBA_spm_mesh_curvature(M)
% Compute a crude approximation of the curvature of a surface mesh
% FORMAT C = spm_mesh_curvature(M)
% M        - a patch structure
%
% C        - curvature vector
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_curvature.m 3135 2009-05-19 14:49:42Z guillaume $

A = VBA_spm_mesh_adjacency(M);
A = sparse(1:size(M.vertices,1),1:size(M.vertices,1),1./sum(A,2)) * A;

C = (A-speye(size(A))) * double(M.vertices);
N = VBA_spm_mesh_normals(M);
C = sign(sum(N.*C,2)) .* sqrt(sum(C.*C,2));

author: MBB-team

repo: VBA-toolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This test compares the result of the inverse dynamic function using both:
%   a) the Newton-Euler algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_forwarddynamics()
close all
global robot g q0 qd0 tau total_simulation_time

g=[0  -9.81 0]'; %y0 axis
%initial position and joint speed
q0 = [0 0]';
qd0 = [0 0]';
tau = [0 0]';%no torques applied
total_simulation_time = 5;
%select friction or not
robot.dynamics.friction = 1;

fprintf('\nTHE SIMULATION PRESENTS THE ROBOT AT AN INITIAL POSITION WHEN NO TORQUES ARE APPLIED\n')

%load robot parameters
robot=load_robot('example', '2dofplanar');

% drawrobot3d(robot, q0);
% adjust_view(robot);

[t q qd] = forward_test_A();

figure, plot(t, q, 'r'), hold

[t q qd] = forward_test_B();

plot(t, q, 'b')

speed = 10
animate(robot,[q(1,1:speed:length(q)); q(2,1:speed:length(q))])

function [t q qd] = forward_test_A()
global robot g q0 qd0 tau total_simulation_time

fprintf('\nCOMPUTING FORWARD DYNAMICS (this may take a while)')

%this may take a while, since it requires integration
%of the acceleration at each time step
[t q qd] = forwarddynamic(robot, total_simulation_time, q0, qd0, tau, g, []);

function [t q qd] = forward_test_B()
global total_simulation_time

[t, y] = runge_kutta(@forward_dynamic_robot2, [0 0 0 0]', [0 total_simulation_time], 0.01);

%Assign joint variables and joint speeds.
q = y(1:2,:);
qd = y(3:4,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to simulate the movement of a 2 DOF robot
% Called from function exerciseD()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = forward_dynamic_robot2(t, y)
global tau g robot
%we must return the solution of
% [dx1/dt; dx2/dt], in this case [qd; qdd]
t
qdd=forwarddynamics_2dofplanar(robot, y(1:2,1), y(3:4,1), tau, -sum(g), [0 0 0 0 0 0]);
xd = [y(3:4,1); qdd];

repo: ARTE


text:
function [I,p] = LEMS1im_read()
% LEMS1im_read: Read one image for LEMS correction
%
% [I,p] = LEMS1im_read()
% I: selected image
% p: parameters
%         info: 'N/A'   % 'N/A' if not dicom
%         Imin: 2       % for contrast
%         Imax: 255
%         rect: [10.2942 13.0805 113.4776 118.8813] % for croping I0
%     filename: 'arlie.pd.im256gl4.tif' % I0 location
%     pathname: [1x64 char]
%
% Olivier Salvado, Case, 15-nov-06
% v0: from stepA_1 in LEMS 2D pipeline

disp(' ')
disp([' --- Read one image for LEMS1im...'])
% colormap(gray(256)),clf
% set(gcf,'BackingStore','on')
% debug = 0;
% warning off MATLAB:divideByZero

%% ask for the file
wiam = cd;
% cd('C:\Documents and Settings\oxs13\My Documents\0 Main backup\Matlab Codes\Data for experiment')
cd('C:\')
[filename, filepath] = uigetfile('*.*', 'Select image to correct');
cd(wiam);
disp(['     In folder:' filepath])

%% read the file
try
    % --- try to read a Dicom
    info = dicominfo([filepath filename]);
    I0 = dicomread([filepath filename]);
catch
    % --- if there is an error, it must be an image.
    info = 'N/A';
    I0 = imread([filepath filename]);
end

%% check if it is a rgb image
if size(I0,3)>1,
    disp('This image has multiple channels, try to convert to gray')
    I0 = rgb2gray(I0);
end

%% convert to [0...255]
I0 = single(I0);
Imin = min(I0(:));
Imax = max(I0(:));
I0 = I0-Imin;
I = I0/(Imax-Imin)*255;

%% update the outputs
p.info = info;
p.Imin = Imin;
p.Imax = Imax;
p.filename = filename;
p.filepath = filepath;

disp('     ... Data loaded')

author: Sable

repo: mcbench-benchmarks


text:
% Test file for trigtech/diffmat.m

function pass = test_diffmat(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();
valsDiscretizationclass = trigcolloc();

%% 1st-order derivative:

%% Force n to be odd:
op = @(x) exp(cos(pi*x));
f = testclass.make(op, [], pref);
n = length(f);
if ( ~rem(n,2) )
    n = n+1;
    x = trigpts(n);
    v = op(x);
else
    x = trigpts(n);
    v = get(f, 'values');
end
D = valsDiscretizationclass.diffmat(n, 1);
df = D*v;
df_exact = @(x) -pi*sin(pi*x).*exp(cos(pi*x));
err = df_exact(x) - df;
pass(1) = ( norm(err, inf) < 1e3*vscale(f)*eps );

a = 10; b = 20;
f = testclass.make(@(x) cos(a*pi*sin(b*pi*x)), [], pref);
n = length(f);
if ( ~rem(n,2) )
    n = n+1;
    x = trigpts(n);
    v = op(x);
else
    x = trigpts(n);
    v = get(f, 'values');
end
D = valsDiscretizationclass.diffmat(n, 1);
df = D*v;
df_exact = @(x) -pi^2*a*b*cos(b*pi*x).*sin(a*pi*sin(b*pi*x));
err = df_exact(x) - df;
pass(2) = ( norm(err, inf) < 1e7*vscale(f)*eps );

%% Force n to be even:
op = @(x) exp(-50*x.^2);
f = testclass.make(op, [], pref);
n = length(f);
if ( rem(n,2) )
    n = n+1;
    x = trigpts(n);
    v = op(x);
else
    x = trigpts(n);
    v = get(f, 'values');
end
D = valsDiscretizationclass.diffmat(n, 1);
df = D*v;
df_exact = @(x) - 100*x.*exp(-50*x.^2);
err = df_exact(x) - df;
pass(3) = ( norm(err, inf) < 1e3*vscale(f)*eps );

a1 = 4; b1 = 3; a2 = 6; b2 = 4;
op = @(x) cos(a1*pi*sin(b1*pi*x)) + 1i*cos(a2*pi*sin(b2*pi*x));
f = testclass.make(op, [], pref);
n = length(f);
if ( rem(n,2) )
    n = n+1;
    x = trigpts(n);
    v = op(x);
else
    x = trigpts(n);
    v = get(f, 'values');
end
D = valsDiscretizationclass.diffmat(n, 1);
df = D*v;
df_exact = @(x) -pi^2*a1*b1*cos(b1*pi*x).*sin(a1*pi*sin(b1*pi*x)) ...
    - 1i*pi^2*a2*b2*cos(b2*pi*x).*sin(a2*pi*sin(b2*pi*x));
err = df_exact(x) - df;
pass(4) = ( norm(err, inf) < 1e5*vscale(f)*eps );

%% 2nd-order derivative (odd n & even n):
op = @(x) exp(cos(4*pi*x))-1;
f = testclass.make(op, [], pref);
n = length(f);
x = trigpts(n);
v = get(f, 'values');
D2 = valsDiscretizationclass.diffmat(n, 2);
df2 = D2*v;
df2_exact = @(x) -16*pi^2*exp(cos(4*pi*x)).*(cos(4*pi*x) + cos(4*pi*x).^2 - 1);
err = df2_exact(x) - df2;
pass(5) = ( norm(err, inf) < 5e5*vscale(f)*eps );

n = n+1;
x = trigpts(n);
v = op(x);
D2 = valsDiscretizationclass.diffmat(n, 2);
df2 = D2*v;
err = df2_exact(x) - df2;
pass(6) = ( norm(err, inf) < 5e5*vscale(f)*eps );

%% Check higher-order derivatives:

% Odd p:
op = @(x) sin(pi*x);
f = testclass.make(op, [], pref);
n = length(f);
x = trigpts(n);
v = get(f, 'values');
D5 = valsDiscretizationclass.diffmat(n, 5);
df5 = D5*v;
df5_exact = @(x) pi^5*cos(pi*x);
err = df5_exact(x) - df5;
pass(7) = ( norm(err, inf) < 1e3*vscale(f)*eps );

n = n+1;
x = trigpts(n);
v = op(x);
D5 = valsDiscretizationclass.diffmat(n, 5);
df5 = D5*v;
err = df5_exact(x) - df5;
pass(8) = ( norm(err, inf) < 1e3*vscale(f)*eps );

% Even p:
op = @(x) sin(pi*x);
f = testclass.make(op, [], pref);
n = length(f);
x = trigpts(n);
v = get(f, 'values');
D6 = valsDiscretizationclass.diffmat(n, 6);
df6 = D6*v;
df6_exact = @(x) -pi^6*sin(pi*x);
err = df6_exact(x) - df6;
pass(9) = ( norm(err, inf) < 1e4*vscale(f)*eps );

n = n+1;
x = trigpts(n);
v = op(x);
D6 = valsDiscretizationclass.diffmat(n, 6);
df6 = D6*v;
err = df6_exact(x) - df6;
pass(10) = ( norm(err, inf) < 1e5*vscale(f)*eps );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pot = recursive_combine_pots(pot1, pot2)
% RECURSIVE_COMBINE_POTS recursive combine two potentials
% pot = recursive_combine_pots(pot1, pot2)

pot1 = reduce_pot(pot1);
pot2 = reduce_pot(pot2);
% Recursion is stopped, if recusive-combination is defined by direct combination,
% i.e. if the domain of one potential is disjoint from the head of the other.
if (isempty(myintersect(pot1.domain,pot2.cheaddom))|...
    isempty(myintersect(pot1.cheaddom,pot2.domain)))
    pot = direct_combine_pots(pot1,pot2);
else
    % Test wether one of the set-differences is not empty
    % as defined in Lauritzen99 "Stable Local Computation with Conditional Gaussian Distributions"
    % on page 9
    D12 = mysetdiff(pot1.cheaddom, pot2.domain);
    D21 = mysetdiff(pot2.cheaddom, pot1.domain);
    if (isempty(D12) & isempty(D21))
       assert(0,'Recursive combination is not defined');
    end

    if ~isempty(D12)
        % Calculate the complementary potential for the set
        % D1\D12 as defined in Lauritzen 99, page 9
    keep = mysetdiff(pot1.domain,D12);
        [margpot, comppot] = complement_pot(pot1,keep);
        margpot = reduce_pot(margpot);
        comppot = reduce_pot(comppot);
        pot = direct_combine_pots( recursive_combine_pots(margpot, pot2), comppot);
    elseif ~isempty(D21)
        keep = mysetdiff(pot2.domain,D21);
        [margpot, comppot] = complement_pot(pot2,D21);
        margpot = reduce_pot(margpot);
        comppot = reduce_pot(comppot);
        pot = direct_combine_pots( recursive_combine_pots(pot1, margpot), comppot);
    end
end

author: bayesnet

repo: bnt


text:
function [stepsize newx storedb lsmem lsstats] = ...
                linesearch(problem, x, d, f0, df0, options, storedb, lsmem)
% Standard line search algorithm (step size selection) for descent methods.
%
% function [stepsize newx storedb lsmem lsstats] =
%               linesearch(problem, x, d, f0, df0, options, storedb, lsmem)
%
% Base linesearch algorithm for descent methods, based on a simple
% backtracking method. The search direction provided has to be a descent
% direction, as indicated by a negative df0 = directional derivative of f
% at x along d.
%
% The algorithm is invariant under positive scaling of the cost function
% and under offsetting, that is: if the cost function f is replaced by
% 8*f+3 for example, the returned step size will be the same. Furthermore,
% the returned step size is independent of the norm of the search direction
% vector d: only the direction of d is important.
%
% Below, the step will be constructed as alpha*d, and the step size is the
% norm of that vector, thus: stepsize = alpha*norm_d. The step is executed
% by retracting the vector alpha*d from the current point x, giving newx.
%
% Inputs
%
%  problem : structure holding the description of the optimization problem
%  x : current point on the manifold problem.M
%  d : tangent vector at x (descent direction) -- its norm is irrelevant
%  f0 : cost value at x
%  df0 : directional derivative at x along d
%  options : options structure (see in code for usage)
%  storedb : store database structure for caching purposes
%  lsmem : a special memory holder to give the linesearch the opportunity
%          to "remember" what happened in the previous calls. See below.
%
% Outputs
%
%  stepsize : norm of the vector retracted to reach newx from x.
%  newx : next iterate suggested by the line search algorithm, such that
%         the retraction at x of the vector alpha*d reaches newx.
%  storedb : the (possibly updated) store database structure.
%  lsmem : the (possibly updated) lsmem memory holder.
%  lsstats : statistics about the line search procedure (stepsize, number
%            of trials etc).
%
% About lsmem : It can be anything, and will typically be a matrix or a
%               structure. When first calling the line search, it should be
%               passed as the empty matrix []. For subsequent calls
%               (pertaining to the same solver call), the previously
%               returned lsmem should be passed as input. This memory
%               holder gives the line search a chance to exploit knowledge
%               of previous decisions to make new decisions.
%
% See also: steepestdescent conjugategradients linesearch_adaptive

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%     Sept. 13, 2013 (NB) : User control over the parameters of the
%                           linesearch via the options
%                           ls_contraction_factor, ls_optimism,
%                           ls_suff_decr and ls_max_steps. See in code for
%                           the effect of those.
%     Sept. 13, 2013 (NB) : The automatic direction reversal feature was
%                           removed (it triggered when df0 > 0). Direction
%                           reversal is a decision that needs to be made by
%                           the solver, so it can know about it.
%     Sept. 13, 2013 (NB) : The linesearch is now invariant under rescaling
%                           of the cost function f. That is, if f is
%                           replaced by 8*f (and hence the directional
%                           derivatives of f are scaled accordingly), the
%                           stepsizes computed will not change.
%     Nov. 7, 2013 (NB) :   The linesearch is now invariant under rescaling
%                           of the search direction d. The meaning of
%                           stepsize is also more clear in the comments.
%                           Added a parameter ls_initial_stepsize to give
%                           users control over the first step size trial.

    % Backtracking default parameters. These can be overwritten in the
    % options structure which is passed to the solver.
    default_options.ls_contraction_factor = .5;
    default_options.ls_optimism = 1/.5;
    default_options.ls_suff_decr = 1e-4;
    default_options.ls_max_steps = 25;
    default_options.ls_initial_stepsize = 1;

    options = mergeOptions(default_options, options);

    contraction_factor = options.ls_contraction_factor;
    optimism = options.ls_optimism;
    suff_decr = options.ls_suff_decr;
    max_ls_steps = options.ls_max_steps;
    initial_stepsize = options.ls_initial_stepsize;

    % Compute the norm of the search direction.
    % This is useful to make the linesearch algorithm invariant under the
    % scaling of d. The rationale is that the important information is the
    % search direction, not the size of that vector. The question of how
    % far we should go is precisely what the linesearch algorithm is
    % supposed to answer: the calling algorithm should not need to care.
    norm_d = problem.M.norm(x, d);

    % At first, we have no idea of what the step size should be.
    alpha = NaN;

    % If we know about what happened at the previous step, we can leverage
    % that to compute an initial guess for the step size, as inspired from
    % Nocedal&Wright, p59.
    if isstruct(lsmem) && isfield(lsmem, 'f0')
        % Pick initial step size based on where we were last time,
        alpha = 2*(f0-lsmem.f0)/df0;
        % and go look a little further (or less far), just in case.
        alpha = optimism*alpha;
    end

    % If we have no information about the previous iteration (maybe this is
    % the first one?) or if the above formula gave a too small step size
    % (perhaps it is even negative), then fall back to a user supplied
    % suggestion for the first step size (the "a priori").
    % At any rate, the choice should be invariant under rescaling of the
    % cost function f and of the search direction d, and it should be
    % bounded away from zero for convergence guarantees. We must allow it
    % to be close to zero though, for fine convergence.
    if isnan(alpha) || alpha*norm_d <= eps
        alpha = initial_stepsize/norm_d;
    end

    % Make the chosen step and compute the cost there.
    newx = problem.M.retr(x, d, alpha);
    [newf storedb] = getCost(problem, newx, storedb);
    cost_evaluations = 1;

    % Backtrack while the Armijo criterion is not satisfied
    while newf > f0 + suff_decr*alpha*df0

        % Reduce the step size,
        alpha = contraction_factor * alpha;

        % and look closer down the line
        newx = problem.M.retr(x, d, alpha);
        [newf storedb] = getCost(problem, newx, storedb);
        cost_evaluations = cost_evaluations + 1;

        % Make sure we don't run out of budget
        if cost_evaluations >= max_ls_steps
            break;
        end

    end

    % If we got here without obtaining a decrease, we reject the step.
    if newf > f0
        alpha = 0;
        newx = x;
        newf = f0; %#ok<NASGU>
    end

    % As seen outside this function, stepsize is the size of the vector we
    % retract to make the step from x to newx. Since the step is alpha*d:
    stepsize = alpha * norm_d;

    % Save the situtation faced now so that at the next iteration, we will
    % know something about the previous decision.
    lsmem.f0 = f0;
    lsmem.df0 = df0;
    lsmem.stepsize = stepsize;

    % Save some statistics also, for possible analysis.
    lsstats.costevals = cost_evaluations;
    lsstats.stepsize = stepsize;
    lsstats.alpha = alpha;

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function pde = biharmonicdata
%% BIHARMONICDATA trigonometric data for biharmonic equation
%
%     f = 64*pi^4*sin(2*pi*x)*cos(2*pi*y);
%     u = sin(2*pi*x)*cos(2*pi*y);
%     Du = (2*pi*cos(2*pi*x)*cos(2*pi*y), -2*pi*sin(2*pi*x)*sin(2*pi*y));
%     w = Delta u = -8*pi^2*sin(2*pi*x)*cos(2*pi*y);
%     du/dn = g_N on [0,1]^2.
%
% Added by Jie Zhou.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'exactu',@exactu,'Du',@Du,'exactw',@exactw,'Dw',@Dw,...
             'g_D',@g_D,'g_N',@g_N);

% load data (right hand side function)
function rhs =  f(p)
    x = p(:,1); y = p(:,2);
    rhs =  64*pi^4*sin(2*pi*x).*cos(2*pi*y);
end

% exact solution
function u =  exactu(p)
    x = p(:,1); y = p(:,2);
    u =  sin(2*pi*x).*cos(2*pi*y);
end

% derivative of the exact solution
function uprime =  Du(p)
    x = p(:,1); y = p(:,2);
    uprime(:,1) = 2*pi*cos(2*pi*x).*cos(2*pi*y);
    uprime(:,2) = -2*pi*sin(2*pi*x).*sin(2*pi*y);
end

% exact solution of w = -Delta u
function w =  exactw(p)
    x = p(:,1); y = p(:,2);
    w = -8*pi^2*sin(2*pi*x).*cos(2*pi*y);
end

% derivative of the exact solution w
function wprime =  Dw(p)
    x = p(:,1); y = p(:,2);
    wprime(:,1) = -16*pi^3*cos(2*pi*x).*cos(2*pi*y);
    wprime(:,2) = 16*pi^3*sin(2*pi*x).*sin(2*pi*y);
end

% Dirichlet boundary condition
function u = g_D(p)
    u = exactu(p);
end

% Neumann boundary condition
function f = g_N(p)
    f = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2);
    uprime = [2*pi*cos(2*pi*p(:,1)).*cos(2*pi*p(:,2)) -2*pi*sin(2*pi*p(:,1)).*sin(2*pi*p(:,2))];
    leftbd = (abs(x)<eps);  % n = (-1,0);
    f(leftbd) = - uprime(leftbd,1);
    rightbd = (abs(x-1)<eps); % n = (1,0);
    f(rightbd) = uprime(rightbd,1);
    topbd = (abs(y-1)<eps);   % n = (0,1)
    f(topbd) = uprime(topbd,2);
    bottombd = (abs(y)<eps);% n = (0,-1)
    f(bottombd) = - uprime(bottombd,2);
end
end

repo: ifem


text:
function FN = flipped_normals(V,F)
%FLIPPED_NORMALS Compute the flipped normals of a triangle mesh
%
% FN = flipped_normals(V,F);
%
% Inputs:
%  V,F  triangle mesh
% Outputs:
%  FN  flipped normals of the mesh V,N

%Compute per-face normals.
N = normals(V,F);
N = N ./ normrow(V,F);

%Flip the per-face normals.
FN = -N;

end

author: odedstein

repo: sgi-introduction-course