text:
function [C, status] = parfileCondition(parfiles, nFrames);
%
% [C, status] = parfileCondition(parfiles, <nFrames=set by parfiles>);
%
% Estimate the condition number of an experimental paradigm
% specified by parfiles. This condition number determines
% whether or not time courses can be successfully deconvolved
% for data with this paradigm.
%
% parfiles: cell of parfile paths.
%
% nFrames: optional argument specifying # of frames to take from
% each parfile. If omitted, just takes max # of frames set by the
% parfiles.
%
% C: condition number, equal to sum(X' * X), where X is the
% design matrix constructed from the parfiles.
%
% status: 1 if time courses can be deconvolved, 0 if not.
% Basically, status just tests whether C is less than 10 million.
%
%
%
% ras 02/2006.
status = 0;

params = er_defaultParams;
params.glmHRF = 0;

% create Toeplitz matrix for deconvolution from the parfiles
S = delta_function_from_parfile(parfiles);
X = glm_deconvolution_matrix(S, -4:17);

% This part I've empirically determined from er_selxavg, but
% hope to understand soon:
% iteratively sum the matrix X' * W * X
% (where W is a whitening matrix which we basically always
% set to the identitiy matrix, so I'm only concerned with that
% case).
nFrames = size(X, 1);
nPred = size(X, 2);
nRuns = size(X, 3);
XtX = zeros(nPred, nPred);
for i = 1:nRuns
    XtX = XtX + [X(:,:,i)' * X(:,:,i)];
end

C = cond(XtX);
status = (C < 10^7);

if status==0
    disp('Paradigm is ill-conditioned.')
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [trans, points] = registerPoints3dAffine(points, target, varargin)
% Fit 3D affine transform using iterative algorithm.
%
%   TRANS = registerPoints3dAffine(POINTS, TARGET)
%   Computes the affine transform that maps the shape defines by POINTS
%   onto the shape defined by the points TARGET. Both POINTS and TARGET are
%   N-by-3 array of point coordinates, not necessarily the same size.
%   The result TRANS is a 4-by-4 affine transform.
%
%   TRANS = registerPoints3dAffine(POINTS, TARGET, NITER)
%   Specifies the number of iterations for the algorithm.
%
%   [TRANS, POINTS2] = registerPoints3dAffine(...)
%   Also returns the set of transformed points.
%
%   Example
%     registerPoints3dAffine
%
%   See also
%     transforms3d, fitAffineTransform3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2015-02-24,    using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015 INRA - Cepia Software Platform.

nIters = 10;
if ~isempty(varargin)
    nIters = varargin{1};
end

% keep original points to transform them at each
trans = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];

for i = 1:nIters
    % identify target points for each source point
    inds = findClosestPoint(points, target);
    corrPoints = target(inds, :);

    % compute transform for current iteration
    trans_i = fitAffineTransform3d(points, corrPoints);

    % apply transform, and update cumulated transform
    points = transformPoint3d(points, trans_i);
    trans = trans_i * trans;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function overlay_background(src_folder, dst_folder, bkgFilelist, bkgFolder, clutteredBkgRatio, single_thread)

if nargin < 6
    single_thread = 0;
end
if single_thread
    num_workers = 0;
else
    num_workers = 24;
end

image_files = rdir(fullfile(src_folder,'*/*.png'));
image_num = length(image_files);
fprintf('%d images in total.\n', image_num);
if image_num == 0
    return;
end
sunImageList = importdata(bkgFilelist);
rng('shuffle');

fprintf('Start overlaying images at time %s, it takes for a while...\n', datestr(now, 'HH:MM:SS'));
report_num = 80;
fprintf([repmat('.',1,report_num) '\n\n']);
report_step = floor((image_num+report_num-1)/report_num);
t_begin = clock;
%for i = 1:image_num
parfor(i = 1:image_num, num_workers)
    src_image_file = image_files(i).name;
    try
        [I, ~, alpha] = imread(src_image_file);
    catch
        fprintf('Failed to read %s\n', src_image_file);
    end

    s = size(I);
    fh = s(1); fw = s(2);
    mask = double(alpha) / 255;
    mask = repmat(mask,1,1,3);

    if rand() > clutteredBkgRatio
        I = uint8(double(I) .* mask + double(rand()*255) * (1 - mask));
    else
        while true
            id = randi(length(sunImageList));
            bg = imread(fullfile(bkgFolder, sunImageList{id}));
            s = size(bg);
            bh = s(1); bw = s(2);
            if bh < fh || bw < fw
                %fprintf(1, '.');
                continue;
            end
            if length(s) < 3
                continue;
            end
            break;
        end
        by = randi(bh - fh + 1);
        bx = randi(bw - fw + 1);
        bgcrop = bg(by:(by+fh-1), bx:(bx+fw-1), :);

        I = uint8(double(I) .* mask + double(bgcrop) .* (1 - mask));
    end

    if numel(I) == 0
        fprintf('Failed to overlay %s (empty image after crop)\n', src_image_file);
    else
        dst_image_file = strrep(src_image_file, src_folder, dst_folder);
        dst_image_file = strrep(dst_image_file, '.png', '.jpg');
        [dst_image_file_folder, ~, ~] = fileparts(dst_image_file);
        if ~exist(dst_image_file_folder, 'dir')
            mkdir(dst_image_file_folder);
        end
        %imwrite(I, dst_image_file, 'png', 'Alpha', alpha);
        imwrite(I, dst_image_file, 'jpg');
    end

    if mod(i, report_step) == 0
        fprintf('\b|\n');
    end
end
t_end = clock;
fprintf('%f seconds spent on background overlay!\n', etime(t_end, t_begin));

author: ShapeNet

repo: RenderForCNN


text:
function  omega = maxAngle(oR,varargin)
% get the maximum angle of the fundamental region
%
% Syntax
%   omega = maxAngle(oR)   % maximum angle in the orientation region
%   omega = maxAngle(oR,v) % maximum angle about axis v
%
% Description
% It is important that v is within the fundamental sector
%
% Input
%  oR - @orientationRegion
%  v  - @vector3d
%

if nargin>1 && isa(varargin{1},'vector3d')

  % ignore restrictions on the rotational axis
  %N = oR.N(oR.N.angle < pi<1e-3);
  N = oR.N;

  % project v to the fundamental sector
  dcs = properGroup(disjoint(oR.CS1,oR.CS2));
  if oR.antipodal, dcs = dcs.Laue; end
  h = varargin{1};
  h = reshape(project2FundamentalRegion(h,dcs),size(h));

  if isempty(N)
    omega = repmat(pi,size(h));
  else

    d = dot_outer(N.axis, normalize(h));
    dtan = repmat(-tan(N(:).angle./2),1,size(d,2));
    d(abs(d)<1e-2)=0;
    d = dtan .* d;

    %d = dot_outer(-tan(N.angle/2) .* N.axis, normalize(varargin{1}));

    d = 2*acot(d);
    d(d<0) = pi;

    omega = min(d,[],1);
    omega = reshape(omega,size(h));
    omega(omega<1e-4) = 0;
  end

elseif isempty(oR.V) || check_option(varargin,'complete')

  omega = pi;

else

  omega = max(angle(oR.V,'noSymmetry'));

end

author: mtex-toolbox

repo: mtex

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


text:
function [train_sample1, train_sample2, label1, label2] = gen_train_sample_xqda(label, cam, train_feat)

uni_label = unique(label);
ind1 = [];
ind2 = [];
label1 = [];
label2 = [];
for n = 1:length(uni_label)
    pos = find(label == uni_label(n));
    perm = randperm(length(pos));
    tmp1 = pos(perm(1:floor(length(pos)/2)));
    tmp2 = pos(perm(floor(length(pos)/2)+1:floor(length(pos)/2)*2));
    cam1 = cam(tmp1);
    cam2 = cam(tmp2);
    pos2 = find(cam1~=cam2);
    tmp1 = tmp1(pos2);
    tmp2 = tmp2(pos2);
    ind1 = [ind1; tmp1];
    ind2 = [ind2; tmp2];
    label1 = [label1; repmat(uni_label(n), [length(tmp1), 1])];
    label2 = [label2; repmat(uni_label(n), [length(tmp2), 1])];
end
train_sample1 = train_feat(:, ind1)';
train_sample2 = train_feat(:, ind2)';

repo: MARS-evaluation


text:
function B = dtimatrix(bvalues,bvectors)
% B = dtimatrix(bvalues,bvectors)
%
% Constructs a DTI design matrix from the bvalues and bvectors.
% bvalues is a vector of length N
% bvectors is a matrix either Nx3 or 3xN
%
% B will be N by 7
% The 7th is the mean (all ones)
% The tensor will be constructed using the following regressors
%     1 2 3
%     2 4 5
%     3 5 6
%
% $Id: dtimatrix.m,v 1.2 2011/03/02 00:04:12 nicks Exp $

%
% dtimatrix.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.2 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

if(nargin ~= 2)
  fprintf('B = dtimatrix(bvalues,bvectors)\n');
  return;
end

if(size(bvectors,2) ~= 3) bvectors = bvectors'; end
if(size(bvectors,2) ~= 3)
  fprintf('ERROR: bvectors must be Nx3 or 3xN\n');
  return;
end
Nb = size(bvectors,1);

if(Nb ~= length(bvalues))
  fprintf('ERROR: dimension mismatch between bvectors and bvalues\n');
  return;
end

bvalues = bvalues(:);

B = zeros(Nb,7);
B(:,1) =     bvalues .* bvectors(:,1).*bvectors(:,1);
B(:,2) = 2 * bvalues .* bvectors(:,1).*bvectors(:,2);
B(:,3) = 2 * bvalues .* bvectors(:,1).*bvectors(:,3);
B(:,4) =     bvalues .* bvectors(:,2).*bvectors(:,2);
B(:,5) = 2 * bvalues .* bvectors(:,2).*bvectors(:,3);
B(:,6) =     bvalues .* bvectors(:,3).*bvectors(:,3);
B(:,7) = 1;

return;


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

% Edit the above text to modify the response to help resultManager2
% Last Modified by GUIDE v2.5 11-Jul-2013 01:16:29
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

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

    handles.returns_daily = varargin{1} - 1;
    handles.returns = cumprod(varargin{1});
    handles.portfolio = varargin{2};
    handles.selectedAlgorithms = varargin{3};
    handles.frequency = varargin{4};
    [handles.r handles.c] = size(handles.portfolio{1});
    rd_cumulativeReturns_Callback(hObject, eventdata, handles);

     % Compute the statsand display the important numbers in table using the
    % library functions
    addpath('./lib');

    cumReturns = handles.returns;
    dailyReturns = handles.returns_daily;

    % Get final Values
    [r c] = size(cumReturns);
    handles.finalValues = cumReturns(r,:);

    % Get the mean returnfor every day - This is a simple average
    handles.meanReturns = mean(dailyReturns);

    % Get annualised returns
    denominator = 252/handles.frequency;
    Y = r/denominator;
    handles.annualisedReturns = handles.finalValues.^(1/Y)-1;

    % Get standard deviation- a measureof risk
    handles.standardDeviation = std(dailyReturns);

    % Get annualised standard deviation
    handles.annualisedStandardDeviation = handles.standardDeviation * sqrt(denominator);

    % Get sharpe ratios
    handles.sharpeRatios = sharpe(dailyReturns, handles.finalValues, handles.frequency);

    % Get Sortino ratios
    handles.sortinoRatios = sortino(dailyReturns, 0);

    % Get Value risks at level 5%
    handles.valueAtRisks = var5(dailyReturns);

    % Get Maximum draw down
    handles.mdds = maxDD_general(cumReturns);

    % Get Calmar ratios
    handles.calmars = calmar(cumReturns, handles.mdds, handles.frequency);

    % Fill up the tables
    set(handles.table,'ColumnName',handles.selectedAlgorithms);
    tableData = [handles.finalValues; handles.meanReturns; handles.annualisedReturns; handles.standardDeviation; handles.annualisedStandardDeviation; handles.sharpeRatios; handles.calmars; handles.sortinoRatios; handles.valueAtRisks; handles.mdds];
    set(handles.table,'Data',tableData);

    rmpath('./lib');

    guidata(hObject, handles);

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes resultManager2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = resultManager2_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

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

function rd_sharpe_Callback(hObject, eventdata, handles)
% hObject    handle to rd_sharpe (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rd_sharpe
    clearRadio(hObject, eventdata, handles);
    set(handles.rd_sharpe,'Value', 1);

    window = str2double(get(handles.windowSharpe,'string'));
    if isnan(window) || window < 2
        errorMessage('Window size must be a number greater than or equal to 2');
        return;
    end
    window = round(window);
    addpath('./lib');

    [r c] = size(handles.returns_daily);
    sharpeRatios = zeros(r,c);
    for i = 50:1:r
        start = max(i-window+1, 1);
        finish = i;
        sharpeRatios(i,:) = sharpe(handles.returns_daily(start:finish,:), handles.returns(finish,:)./handles.returns(start,:), handles.frequency);
    end

    rmpath('./lib');

    plot(handles.mainGraph, sharpeRatios);
    xlim(handles.mainGraph, [0 handles.r+1]);
    legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');

function windowSharpe_Callback(hObject, eventdata, handles)
% hObject    handle to windowSharpe (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of windowSharpe as text
%        str2double(get(hObject,'String')) returns contents of windowSharpe as a double

% --- Executes during object creation, after setting all properties.
function windowSharpe_CreateFcn(hObject, eventdata, handles)
% hObject    handle to windowSharpe (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in rd_calmar.
function rd_calmar_Callback(hObject, eventdata, handles)
% hObject    handle to rd_calmar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rd_calmar
    clearRadio(hObject, eventdata, handles);
    set(handles.rd_calmar,'Value', 1);

    window = str2double(get(handles.windowCalmar,'string'));
    if isnan(window) || window < 2
        errorMessage('Window size must be a number greater than or equal to 2');
        return;
    end
    window = round(window);
    addpath('./lib');

    [r c] = size(handles.returns_daily);
    calmarRatios = zeros(r,c);
    for i = 50:1:r
        start = max(i-window+1, 1);
        finish = i;
        mdd = maxDD_general(handles.returns(start:finish,:));
        calmarRatios(i,:) = calmar(handles.returns(start:finish,:), mdd,  handles.frequency);
    end

    rmpath('./lib');

    plot(handles.mainGraph, calmarRatios);
    xlim(handles.mainGraph, [0 handles.r+1]);
    legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');

% --- Executes on button press in rd_var.
function rd_var_Callback(hObject, eventdata, handles)
% hObject    handle to rd_var (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rd_var
    clearRadio(hObject, eventdata, handles);
    set(handles.rd_var,'Value', 1);

    window = str2double(get(handles.windowVar,'string'));
    if isnan(window) || window < 2
        errorMessage('Window size must be a number greater than or equal to 2');
        return;
    end
    window = round(window);
    addpath('./lib');

    [r c] = size(handles.returns_daily);
    vars = zeros(r,c);
    for i = 50:1:r
        start = max(i-window+1, 1);
        finish = i;
        vars(i,:) = var5(handles.returns_daily(start:finish,:));
    end

    rmpath('./lib');
    plot(handles.mainGraph, vars);
    xlim(handles.mainGraph, [0 handles.r+1]);
    legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');

function windowCalmar_Callback(hObject, eventdata, handles)
% hObject    handle to windowCalmar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of windowCalmar as text
%        str2double(get(hObject,'String')) returns contents of windowCalmar as a double

% --- Executes during object creation, after setting all properties.
function windowCalmar_CreateFcn(hObject, eventdata, handles)
% hObject    handle to windowCalmar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in rd_mdd.
function rd_mdd_Callback(hObject, eventdata, handles)
% hObject    handle to rd_mdd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rd_mdd
    clearRadio(hObject, eventdata, handles);
    set(handles.rd_mdd,'Value', 1);

    window = str2double(get(handles.windowMDD,'string'));
    if isnan(window) || window < 2
        errorMessage('Window size must be a number greater than or equal to 2');
        return;
    end
    window = round(window);
    addpath('./lib');

    [r c] = size(handles.returns_daily);
    mdds = zeros(r,c);
    for i = 50:1:r
        start = max(i-window+1, 1);
        finish = i;
        mdds(i,:) = maxDD_general(handles.returns(start:finish,:));
    end

    rmpath('./lib');

    plot(handles.mainGraph, mdds);
    xlim(handles.mainGraph, [0 handles.r+1]);
    legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');

% --- Executes on button press in rd_sortino.
function rd_sortino_Callback(hObject, eventdata, handles)
% hObject    handle to rd_sortino (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rd_sortino
    clearRadio(hObject, eventdata, handles);
    set(handles.rd_sortino,'Value', 1);

    window = str2double(get(handles.windowSortino,'string'));
    if isnan(window) || window < 2
        errorMessage('Window size must be a number greater than or equal to 2');
        return;
    end
    window = round(window);
    addpath('./lib');

    [r c] = size(handles.returns_daily);
    sortinos = zeros(r,c);
    for i = 50:1:r
        start = max(i-window+1, 1);
        finish = i;
        sortinos(i,:) = sortino(handles.returns_daily(start:finish,:),0);
    end

    rmpath('./lib');

    plot(handles.mainGraph, sortinos);
    xlim(handles.mainGraph, [0 handles.r+1]);
    legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');

function windowSortino_Callback(hObject, eventdata, handles)
% hObject    handle to windowSortino (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of windowSortino as text
%        str2double(get(hObject,'String')) returns contents of windowSortino as a double

% --- Executes during object creation, after setting all properties.
function windowSortino_CreateFcn(hObject, eventdata, handles)
% hObject    handle to windowSortino (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function windowVar_Callback(hObject, eventdata, handles)
% hObject    handle to windowVar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of windowVar as text
%        str2double(get(hObject,'String')) returns contents of windowVar as a double

% --- Executes during object creation, after setting all properties.
function windowVar_CreateFcn(hObject, eventdata, handles)
% hObject    handle to windowVar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function windowMDD_Callback(hObject, eventdata, handles)
% hObject    handle to windowMDD (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of windowMDD as text
%        str2double(get(hObject,'String')) returns contents of windowMDD as a double

% --- Executes during object creation, after setting all properties.
function windowMDD_CreateFcn(hObject, eventdata, handles)
% hObject    handle to windowMDD (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in rd_cumulativeReturns.
function rd_cumulativeReturns_Callback(hObject, eventdata, handles)
% hObject    handle to rd_cumulativeReturns (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rd_cumulativeReturns
    clearRadio(hObject, eventdata, handles);
    set(handles.rd_cumulativeReturns,'Value', 1);
    if get(handles.logCheck,'Value') ~=1
        plot(handles.mainGraph, handles.returns);
        xlim(handles.mainGraph, [0 handles.r+1]);
        legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');
        grid on;
    else
        logCheck_Callback(hObject, eventdata, handles);
    end

% --- Executes on button press in logCheck.
function logCheck_Callback(hObject, eventdata, handles)
% hObject    handle to logCheck (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of logCheck

    % Check if user wants to see cumulative returns
    selected  = get(handles.rd_cumulativeReturns,'Value');

    % only if user wants to see the cumulative returns
    if selected == 1
        logPlot = get(handles.logCheck,'Value');

        if logPlot == 1
            semilogy(handles.returns);
        else
            plot(handles.returns);
        end
        xlim(handles.mainGraph, [0 handles.r+1]);
        legend(handles.mainGraph , handles.selectedAlgorithms, 'Location', 'Best');
        grid on;
    end

function clearRadio(hObject, eventdata, handles)
% The function is a helper function to uncheck all radio buttons

    set(handles.rd_cumulativeReturns,'Value', 0);
    set(handles.rd_sharpe,'Value', 0);
    set(handles.rd_calmar,'Value', 0);
    set(handles.rd_sortino,'Value', 0);
    set(handles.rd_var,'Value', 0);
    set(handles.rd_mdd,'Value', 0);

% --- Executes on button press in dailyReturnsButton.
function dailyReturnsButton_Callback(hObject, eventdata, handles)
% hObject    handle to dailyReturnsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

    first = str2double(get(handles.dailyStart,'string'));
    last = str2double(get(handles.dailyEnd,'string'));

    if isnan(first) || isnan(last) || first < 1 || last - first < 1
        errorMessage('Error in entering start and end days. End day must be greater than start day. Start day must be positive integer.');
        return;
    end
    first = round(first);
    last = round(last);

    [r c] = size(handles.returns_daily);
    last = min(last,r);
    bar(handles.mainGraph, handles.returns_daily(first:last,:));
    xlim([0 last-first+2]);
    legend(handles.mainGraph,handles.selectedAlgorithms, 'Location', 'Best');

% --- Executes on button press in portfolioAllocationButton.
function portfolioAllocationButton_Callback(hObject, eventdata, handles)
% hObject    handle to portfolioAllocationButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    first = str2double(get(handles.portfolioStart,'string'));
    last = str2double(get(handles.portfolioEnd,'string'));

    if isnan(first) || isnan(last) || first < 1 || last - first < 1
        errorMessage('Error in entering start and end days. End day must be greater than start day. Start day must be positive integer.');
        return;
    end
    first = round(first);
    last = round(last);

    [r c] = size(handles.portfolio{1});
    last = min(last,r);

    % Number of algorithms
    [tmp n] = size(handles.portfolio);
    for i = 1:1:n
        portfolio = handles.portfolio{i};
        portfolio = portfolio(first:last,:);
        expected(:,i) = mean(portfolio)';
        deviation(:,i) = std(portfolio)';
    end
    errorbar(handles.mainGraph, expected, deviation,'x');
    xlim([0 c+1]);
    minMean = min(min(expected));
    maxMean = max(max(expected));
    maxDev = max(max(deviation));
    legend(handles.mainGraph, handles.selectedAlgorithms, 'Location', 'Best');

function portfolioStart_Callback(hObject, eventdata, handles)
% hObject    handle to portfolioStart (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of portfolioStart as text
%        str2double(get(hObject,'String')) returns contents of portfolioStart as a double

% --- Executes during object creation, after setting all properties.
function portfolioStart_CreateFcn(hObject, eventdata, handles)
% hObject    handle to portfolioStart (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function portfolioEnd_Callback(hObject, eventdata, handles)
% hObject    handle to portfolioEnd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of portfolioEnd as text
%        str2double(get(hObject,'String')) returns contents of portfolioEnd as a double

% --- Executes during object creation, after setting all properties.
function portfolioEnd_CreateFcn(hObject, eventdata, handles)
% hObject    handle to portfolioEnd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function dailyStart_Callback(hObject, eventdata, handles)
% hObject    handle to dailyStart (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of dailyStart as text
%        str2double(get(hObject,'String')) returns contents of dailyStart as a double

% --- Executes during object creation, after setting all properties.
function dailyStart_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dailyStart (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function dailyEnd_Callback(hObject, eventdata, handles)
% hObject    handle to dailyEnd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of dailyEnd as text
%        str2double(get(hObject,'String')) returns contents of dailyEnd as a double

% --- Executes during object creation, after setting all properties.
function dailyEnd_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dailyEnd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
% Determine sample time from position data
%
function sampleTime = sampleTimeFromData(positions)
    i = 3;
    while isnan(positions(i, 1)) || isnan(positions(i-1, 1))
        i = i + 1;
    end
    sampleTime = positions(i, 1) - positions(i-1, 1);
end

author: kavli-ntnu


text:
% mlp_classify
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [c, posterior] = mlp_classify(M, x0, Q0, raw)

if nargin < 3
    Q0 = [];
end

if nargin < 4
    raw = 0;
end

layers = M.structure.layers;
n_layers = length(layers);

posterior = x0;

if isfield(M, 'dbm') && M.dbm.use
    for l = 2:n_layers
        if M.dropout.use && l > 2
            posterior = posterior * bsxfun(@times, (M.W{l-1}), 1 - M.dropout.probs{l-1});
        else
            posterior = posterior * M.W{l-1};
        end

        if l < n_layers-1
            posterior = posterior + Q0{l+1} * (M.dbm.W{l})';
        end
        posterior = bsxfun(@plus, posterior, M.biases{l}');

        if l < n_layers
            posterior = sigmoid(posterior, M.hidden.use_tanh);
        end

        if l == n_layers && M.output.binary
            posterior = softmax(posterior);
        end
    end
else
    for l = 2:n_layers
        if M.dropout.use && l > 2
            posterior = bsxfun(@plus, posterior * (M.W{l-1}/2), M.biases{l}');
        else
            posterior = bsxfun(@plus, posterior * M.W{l-1}, M.biases{l}');
        end

        if l < n_layers
            posterior = sigmoid(posterior, M.hidden.use_tanh);
        end

        if l == n_layers && M.output.binary
            posterior = softmax(posterior);
        end
    end
end

if raw
    c = posterior;
else
    [maxp, c] = max(posterior, [], 2);
end

author: kyunghyuncho

repo: deepmat


text:
function frac = i4vec_frac ( n, a, k )

%*****************************************************************************80
%
%% I4VEC_FRAC searches for the K-th smallest entry in an I4VEC.
%
%  Discussion:
%
%    Hoare's algorithm is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Input, integer A(N), the array to search.
%
%    Input, integer K, the fractile to be sought.  If K = 1, the minimum
%    entry is sought.  If K = N, the maximum is sought.  Other values
%    of K search for the entry which is K-th in size.  K must be at
%    least 1, and no greater than N.
%
%    Output, integer FRAC, the value of the K-th fractile of A.
%
  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4VEC_FRAC - Fatal error!\n' );
    fprintf ( 1, '  Illegal nonpositive value of N = %d\n', n );
    error ( 'I4VEC_FRAC - Fatal error!' );
  end

  if ( k <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4VEC_FRAC - Fatal error!\n' );
    fprintf ( 1, '  Illegal nonpositive value of K = %d\n', k );
    error ( 'I4VEC_FRAC - Fatal error!' );
  end

  if ( n < k )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4VEC_FRAC - Fatal error!\n' );
    fprintf ( 1, '  Illegal N < K, K = %d\n', k );
    error ( 'I4VEC_FRAC - Fatal error!' );
  end

  left = 1;
  iryt = n;

  while ( 1 )

    if ( iryt <= left )
      frac = a(k);
      break;
    end

    x = a(k);
    i = left;
    j = iryt;

    while ( 1 )

      if ( j < i )
        if ( j < k )
          left = i;
        end
        if ( k < i )
          iryt = j;
        end
        break;
      end
%
%  Find I so that X <= A(I)
%
      while ( a(i) < x )
        i = i + 1;
      end
%
%  Find J so that A(J) <= X
%
      while ( x < a(j) )
        j = j - 1;
      end

      if ( i <= j )

        temp = a(i);
        a(i) = a(j);
        a(j) = temp;

        i = i + 1;
        j = j - 1;

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gX = rbfardjitKernDiagGradX(kern, X)

% RBFARDJITKERNDIAGGRADX Gradient of RBFARDJIT kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the automatic relevance determination radial basis function kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : rbfard2KernParamInit, kernDiagGradX, rbfard2kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Michalis K. Titsias, 2009

% KERN

gX = zeros(size(X));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ a, ifact ] = i4vec_red ( n, a, incx )

%*****************************************************************************80
%
%% I4VEC_RED divides out common factors in an I4VEC.
%
%  Discussion:
%
%    If A is a simple vector, then it has dimension N.
%
%    If A is a row of a matrix, then INCX will not be 1, and
%    the actual dimension of A is at least 1+(N-1)*INCX.
%
%    On output, the entries of A have no common factor
%    greater than 1.
%
%    If A is a simple vector, then INCX is 1, and we simply
%    check the first N entries of A.
%
%    If A is a row of a matrix, then INCX will be the number
%    of rows declared in the matrix, in order to allow us to
%    "skip" along the row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer A(*), the vector to be reduced.
%
%    Input, integer INCX, the distance between successive
%    entries of A that are to be checked.
%
%    Output, integer A(*), the reduced vector.
%
%    Output, integer IFACT, the common factor that was divided out.
%

%
%  Find the smallest nonzero value.
%
  ifact = 0;
  indx = 1;

  for i = 1 : n

    if ( a(indx) ~= 0 )

      if ( ifact == 0 )
        ifact = abs ( a(indx) );
      else
        ifact = min ( ifact, abs ( a(indx) ) );
      end

    end

    indx = indx + incx;

  end

  if ( ifact == 0 )
    return
  end
%
%  Find the greatest common factor of the entire vector.
%
  indx = 1;
  for i = 1 : n
    ifact = i4_gcd ( a(indx), ifact );
    indx = indx + incx;
  end

  if ( ifact == 1 )
    return
  end
%
%  Divide out the common factor.
%
  indx = 1;
  for i = 1 : n
    a(indx) = a(indx) / ifact;
    indx = indx + incx;
  end

  return
end


text:
function [ state ] = nonBlocking_movePTPHomeJointSpace( t , relVel)
%% This non-blocking function is used for performing point to point motion
%  in joint space to the home position of the robot, for the KUKA iiwa 7 R 800.

%% Syntax:
% [ state ] = movePTPHomeJointSpace( t , relVel)

%% About:
% This function is used to move the robot to home position, the motion is
% performed in joint space.

%% Arreguments:
% t: is the TCP/IP connection
% relVel: is a double, from zero to one, specifies the over-ride relative
% velocity.

% Copy right, Mohammad SAFEEA, 9th of May 2017

% check also the function: nonBlocking_isGoalReached()

    theCommand=['jRelVel_',num2str(relVel),'_']; % set over ride.
    fprintf(t, theCommand);
    message=fgets(t);

    jPos=[];
    for i=1:7
        jPos{i}=0;
    end

    global paramName;
    global paramVal;
    paramName='jpos';
    paramVal=jPos;

    sendJointsPositions( t ,jPos); % send destination joint positions.

    theCommand='doPTPinJS';
    fprintf(t, theCommand); % start the point to point motion.
    message=fgets(t);
    state=checkAcknowledgment(message);

end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [D,rho,dD,drho,d2Phi] = MI(Rc,Tc,Omega,m,varargin)
h   = Omega./m;
hd  = prod(h);
PARA = varargin{1};
doDerivative = (nargout > 3);

% D   = phi(rho(y))
% dD  = dPhi(res(y)) * dRes(y)
% d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider

% example MI:
%   phi   = res' * log(res + tol) + ...
%   dPhi  = log(res + tol) + res./(res + tol) + ...
%   d2Phi = (res + 2*tol)./(res + tol)^2 + ...
%   res   = rho(T,R)
%   dRes  = drho, see pdfestimate

tol         = PARA.entropyTol;
[rho,drho]  = pdfestimate(Rc,Tc,PARA,doDerivative);
[n1,n2]     = size(rho);

rhoR = sum(rho,2);
rhoT = sum(rho,1)';
rho  = rho(:);

D    = rhoR'*log(rhoR+tol)+rhoT'*log(rhoT+tol) - rho'*log(rho+tol);

if ~doDerivative, return; end;

SR    = sparse(kron(ones(1,n2),speye(n1,n1)));
ST    = sparse(kron(speye(n2,n2),ones(1,n1)));

dPhi  = ...
   (log(rhoR+tol)+rhoR./(rhoR+tol))'*SR ...
  +(log(rhoT+tol)+rhoT./(rhoT+tol))'*ST ...
  -(log(rho +tol)+rho ./(rho +tol))';

dD    = dPhi * drho;

d2Phi = ...
  SR'*sdiag((rhoR + 2*tol)./(rhoR+tol).^2)*SR ...
  +ST'*sdiag((rhoT + 2*tol)./(rhoT+tol).^2)*ST ...
  -sdiag((rho + 2*tol)./(rho+tol).^2);

a     = 1/sqrt(PARA.ngvR*PARA.ngvT);
a     = 1e0;
d2Phi = - a * d2Phi;

return;


text:
function [vals, timeaxis] = GetYamlVals(yamldata)
% this function converts data formatted in yaml style (cells containing timestamps and values)
% into matlab user friendly matrices.

% obtain number of samples
n = max(size(yamldata));

if n
    if not(iscell(yamldata{1}))
        timeaxis = double(yamldata{1});
        vals   = cell2mat(yamldata(2:end));
    else

        % create output matrices
        timeaxis =  NaN*ones(n,1);
        if n % only if there are some elements of timeaxis
            vals = NaN*ones(n,numel(yamldata{1})-1);
        end
        for i=1:n
            timeaxis(i) = double(yamldata{i}{1});
            vals(i,:)   = cell2mat(yamldata{i}(2:end));
        end
    end

end

end % end of function

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function r8poly3_root_test ( )

%*****************************************************************************80
%
%% R8POLY3_ROOT_TEST tests R8POLY3_ROOT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 4;

  a_test = [ 1.0, 9.0, 1.0, 1.0 ];
  b_test = [ -6.0, -36.0, -5.0, -8.0  ];
  c_test = [ 11.0, 54.0, 8.0, 25.0  ];
  d_test = [ -6.0, -27.0, -4.0, -26.0  ];
%
%  1: Three distinct real roots, 1, 2, 3.
%  2: One repeated real root, 1.5, 1.5, 1.5.
%  3: Two real roots, one repeated, 1, 2, 2.
%  4: One real root, a complex conjugate pair, 2, 3+2I, 3-2I.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8POLY3_ROOT_TEST\n' );
  fprintf ( 1, '  R8POLY3_ROOT finds roots of cubic equations.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    a = a_test(test);
    b = b_test(test);
    c = c_test(test);
    d = d_test(test);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Polynomial coefficients:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  A = %f, B = %f, C = %f, D = %f\n', a, b, c, d );

    [ r1, r2, r3 ] = r8poly3_root ( a, b, c, d );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Roots:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  %f\n', r1 );
    fprintf ( 1, '  %f\n', r2 );
    fprintf ( 1, '  %f\n', r3 );

  end

  return
end


text:
function ff = ifft(data)
% ff = ifft(data) does the inverse fourier transform of the data (power of 2 only)

% This m-file implementation of the fft is meant solely for use with adiff
% objects; it's far too inefficient for general use.

ff = fft(conj(data))/length(data);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [str] = printHyperOrbitsMFMSToTextbox(hHyperbolicOrbitsText, waypoints, hOrbitDepart, orbitsIn, orbitsOut, vInfArrive, OUnitVector, vInfMag, form, form2, paddLen)
%printHyperOrbitsMFMSToTextbox Summary of this function goes here
%   Detailed explanation goes here

    hRule = getHRule();

    [vInfRA,vInfDec,~] = cart2sph(OUnitVector(1),OUnitVector(2),OUnitVector(3));
    vInfRA = rad2deg(vInfRA);
    vInfDec = rad2deg(vInfDec);

    str = {};
    str{end+1} = ['Hyperbolic Departure Orbit from ', cap1stLetter(waypoints{1}.name)];
    str{end+1} = hRule;
    str{end+1} = [paddStr('Semi-major Axis = ',paddLen), num2str(hOrbitDepart(1), form), ' km'];
    str{end+1} = [paddStr('Eccentricity = ', paddLen), num2str(hOrbitDepart(2), form2)];
    str{end+1} = [paddStr('Inclination = ',paddLen), num2str(rad2deg(AngleZero2Pi(hOrbitDepart(3))), form), ' deg'];
    str{end+1} = [paddStr('Right Ascension of AN = ',paddLen), num2str(rad2deg(AngleZero2Pi(hOrbitDepart(4))), form), ' deg'];
    str{end+1} = [paddStr('Argument of Periapse = ',paddLen), num2str(rad2deg(AngleZero2Pi(hOrbitDepart(5))), form), ' deg'];
    str{end+1} = '---------------------';
    str{end+1} = [paddStr('Out. Hyp. Vel. Vect Rt. Asc. = ',paddLen), num2str(vInfRA, form), ' deg'];
    str{end+1} = [paddStr('Out. Hyp. Vel. Vect Declin. = ',paddLen), num2str(vInfDec, form), ' deg'];
    str{end+1} = [paddStr('Out. Hyp. Vel. Magnitude = ',paddLen), num2str(vInfMag, form2), ' km/s'];

    for(i=2:length(waypoints)) %#ok<*NO4LP>
        if(i < length(waypoints))
            [~, flyByRp] = computeApogeePerigee(orbitsIn(i-1,1), orbitsIn(i-1,2));

            str{end+1} = hRule;
            str{end+1} = ['Inbound Hyperbolic Flyby Orbit to ', cap1stLetter(waypoints{i}.name)];
            str{end+1} = hRule;
            str{end+1} = [paddStr('Semi-major Axis = ',paddLen), num2str(orbitsIn(i-1,1), form), ' km'];
            str{end+1} = [paddStr('Eccentricity = ', paddLen), num2str(orbitsIn(i-1,2), form2)];
            str{end+1} = [paddStr('Inclination = ',paddLen), num2str(rad2deg(AngleZero2Pi(orbitsIn(i-1,3))), form), ' deg'];
            str{end+1} = [paddStr('Right Ascension of AN = ',paddLen), num2str(rad2deg(AngleZero2Pi(orbitsIn(i-1,4))), form), ' deg'];
            str{end+1} = [paddStr('Argument of Periapse = ',paddLen), num2str(rad2deg(AngleZero2Pi(orbitsIn(i-1,5))), form), ' deg'];
            str{end+1} = [paddStr('Periapse Radius = ',paddLen), num2str(flyByRp, form), ' km'];

            hOrbit = orbitsOut(i-1,:);
            [~, OUnitVector, vInfMag] = computeHyperSVectOVect(hOrbit(1), hOrbit(2), hOrbit(3), hOrbit(4), hOrbit(5), 0.0, waypoints{i}.gm);
            OUnitVector = normVector(OUnitVector);
            [vInfRA,vInfDec,~] = cart2sph(OUnitVector(1),OUnitVector(2),OUnitVector(3));
            vInfRA = rad2deg(vInfRA);
            vInfDec = rad2deg(vInfDec);

            str{end+1} = hRule;
            str{end+1} = ['Outbound Hyperbolic Flyby Orbit from ', cap1stLetter(waypoints{i}.name)];
            str{end+1} = hRule;
            str{end+1} = [paddStr('Semi-major Axis = ',paddLen), num2str(orbitsOut(i-1,1), form), ' km'];
            str{end+1} = [paddStr('Eccentricity = ', paddLen), num2str(orbitsOut(i-1,2))];
            str{end+1} = [paddStr('Inclination = ',paddLen), num2str(rad2deg(AngleZero2Pi(orbitsOut(i-1,3))), form), ' deg'];
            str{end+1} = [paddStr('Right Ascension of AN = ',paddLen), num2str(rad2deg(AngleZero2Pi(orbitsOut(i-1,4))), form), ' deg'];
            str{end+1} = [paddStr('Argument of Periapse = ',paddLen), num2str(rad2deg(AngleZero2Pi(orbitsOut(i-1,5))), form), ' deg'];
            str{end+1} = [paddStr('Periapse Radius = ',paddLen), num2str(flyByRp, form), ' km'];
            str{end+1} = '---------------------';
            str{end+1} = [paddStr('Out. Hyp. Vel. Vect Rt. Asc. = ',paddLen), num2str(vInfRA, form), ' deg'];
            str{end+1} = [paddStr('Out. Hyp. Vel. Vect Declin. = ',paddLen), num2str(vInfDec, form), ' deg'];
            str{end+1} = [paddStr('Out. Hyp. Vel. Magnitude = ',paddLen), num2str(vInfMag, form2), ' km/s'];
        else
            [vInfRA,vInfDec,~] = cart2sph(vInfArrive(1),vInfArrive(2),vInfArrive(3));
            vInfRA = rad2deg(vInfRA);
            vInfDec = rad2deg(vInfDec);

            str{end+1} = hRule;
            str{end+1} = ['Inbound Hyperbolic Orbit to ', cap1stLetter(waypoints{i}.name)];
            str{end+1} = hRule;
            str{end+1} = [paddStr('Inb. Hyp. Vel. Vect Rt. Asc. = ',paddLen), num2str(vInfRA, form), ' deg'];
            str{end+1} = [paddStr('Inb. Hyp. Vel. Vect Declin. = ',paddLen), num2str(vInfDec, form), ' deg'];
            str{end+1} = [paddStr('Inb. Hyp. Vel. Magnitude = ',paddLen), num2str(norm(vInfArrive), form), ' km/s'];
%             str{end+1} = [paddStr('Hyperbolic Excess Vel. = ',paddLen), num2str(norm(vInfArrive), form), ' km/s'];
        end
    end

    set(hHyperbolicOrbitsText,'String',str);
end


text:
classdef PTKDicomImage < PTKImage
    % PTKDicomImage. A class for holding a 3D medical image volume.
    %
    %     PTKDicomImage is inherited from the fundamental PTKImage image class
    %     used by the TD MIM Toolkit. It provides additional routines and
    %     metadata associated with medical imaging data, most commonly imported
    %     from the DICOM standard.
    %
    %     In general, you do not need to create PTKDicomImages using the
    %     constructor. You should import data using function such as
    %     MimLoadImageFromDicomFiles() which correctly set the metadata.
    %     Thereafter, use .Copy() to make copies of the image, and .BlankCopy()
    %     followed by .ChangeRawImage() to create derived images. See PTKImage.m
    %     for an example of creating functions which modify images contained in
    %     PTKImage classes.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        Modality
        RescaleSlope
        RescaleIntercept
        RescaleUnits = ''
        IsCT = false
        IsMR = false
        MetaHeader
        StudyUid
    end

    methods (Static)
        function new_dicom_image = CreateDicomImageFromMetadata(original_image, metadata, slice_thickness, global_origin_mm, padding_value, reporting)

            % Our voxel size should be equal to SpacingBetweenSlices. However,
            % some scanners incorrectly set SpacingBetweenSlices as the
            % gap between slices, instead of the centre-to-centre distance.
            % To ensure correct reconstruction interval, should really use
            % patient image Position tag (0020,0032)
            % http://www.itk.org/pipermail/insight-users/2005-September/014
            % 711.html
            voxelsize_z = slice_thickness;

            if ndims(original_image.RawImage) > 3
                reporting.ShowWarning('PTKDicomImage:ColourDicomImageNotSupported', 'Colour Dicom images are not supported. Converting to greyscale', []);
                original_image.RawImage = mean(original_image.RawImage, 4);
            end

            if isfield(metadata, 'PixelSpacing')
                voxel_size = [metadata.PixelSpacing' voxelsize_z];
            else
                reporting.ShowWarning('PTKDicomImage:NoPixelSpacing', 'Could not determine the pixel size from the image. I am assuming 1x1mm');
                voxel_size = [1 1 voxelsize_z];
            end

            if isfield(metadata, 'ImageOrientationPatient')
                [new_dimension_order, flip] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation(metadata.ImageOrientationPatient, reporting);
                if isa(original_image, 'CoreWrapper')
                    original_image.RawImage = permute(original_image.RawImage, new_dimension_order);
                    for dimension_index = 1 : 3
                        if flip(dimension_index)
                            original_image.RawImage = flipdim(original_image.RawImage, dimension_index);
                        end
                    end
                else
                    original_image = permute(original_image, new_dimension_order);
                    for dimension_index = 1 : 3
                        if flip(dimension_index)
                            original_image = flipdim(original_image, dimension_index);
                        end
                    end
                end
                voxel_size = voxel_size(new_dimension_order);
            end
            if ~isempty(global_origin_mm)
                global_origin_mm = global_origin_mm([2 1 3]);
            end

            if isfield(metadata, 'RescaleSlope')
                rescale_slope = metadata.RescaleSlope;
            else
                rescale_slope = [];
            end

            if isfield(metadata, 'RescaleIntercept')
                rescale_intercept = metadata.RescaleIntercept;
            else
                rescale_intercept = [];
            end

            if isa(original_image, 'CoreWrapper')
                new_dicom_image = PTKDicomImage( ...
                    original_image.RawImage, rescale_slope, rescale_intercept, voxel_size, metadata.Modality, metadata.StudyInstanceUID, metadata ...
                    );
            else
                new_dicom_image = PTKDicomImage( ...
                    original_image, rescale_slope, rescale_intercept, voxel_size, metadata.Modality, metadata.StudyInstanceUID, metadata ...
                    );
            end
            patient_name = '';
            study_description = '';
            series_description = '';
            if isfield(metadata, 'PatientName')
                if isfield(metadata.PatientName, 'FamilyName')
                    patient_name = ['', metadata.PatientName.FamilyName, ' '];
                end
            end
            if isfield(metadata, 'StudyDescription')
                study_description = ['/ ', metadata.StudyDescription, ' '];
            end
            if isfield(metadata, 'SeriesDescription')
                series_description = ['/ ', metadata.SeriesDescription, ' '];
            end

            new_dicom_image.Title = [patient_name, study_description, series_description];
            new_dicom_image.GlobalOrigin = global_origin_mm;
            new_dicom_image.PaddingValue = padding_value;
        end
    end

    methods
        % The meta_data argument is optional, and is only used for saving DICOM
        % files. Other arguments are compulsory
        function obj = PTKDicomImage(original_image, rescale_slope, rescale_intercept, voxel_size, modality, study_uid, meta_data)

            obj = obj@PTKImage(original_image, PTKImageType.Grayscale, voxel_size);

            obj.StudyUid = study_uid;
            obj.Modality = modality;

            if exist('meta_data', 'var')
                obj.MetaHeader = meta_data;
            end

            if strcmp(modality, 'CT')
                obj.IsCT = true;
                obj.RescaleSlope = rescale_slope;
                obj.RescaleIntercept = rescale_intercept;
                obj.RescaleUnits = 'HU';
            elseif strcmp(modality, 'MR')
                obj.IsMR = true;
                obj.RescaleSlope = [];
                obj.RescaleIntercept = [];
            end
        end

        function [value, units] = GetRescaledValue(obj, global_coords)
            if obj.IsCT && CoreCompareUtilities.CompareEnumName(obj.ImageType, PTKImageType.Grayscale)
                value = obj.GreyscaleToHounsfield(obj.GetVoxel(global_coords));
                units = obj.RescaleUnits;
            else
                [value, units] = GetRescaledValue@PTKImage(obj, global_coords);
            end
        end

        function units_rescaled = GrayscaleToRescaled(obj, units_greyscale)
            if obj.IsCT
                units_rescaled = obj.GreyscaleToHounsfield(units_greyscale);
            else
                units_rescaled = units_greyscale;
            end
        end

        function units_greyscale = RescaledToGrayscale(obj, units_rescaled)
            if obj.IsCT
                units_greyscale = obj.HounsfieldToGreyscale(units_rescaled);
            else
                units_greyscale = units_rescaled;
            end
        end

        function units_greyscale = HounsfieldToGreyscale(obj, units_hu)
            % Conversion from Hounsfield units to image raw intensity values:
            if obj.IsCT
                if isempty(obj.RescaleIntercept) || isempty(obj.RescaleSlope)
                    units_greyscale = units_hu;
                else
                    units_greyscale = (units_hu - cast(obj.RescaleIntercept, 'like', units_hu)/cast(obj.RescaleSlope, 'like', units_hu));
                end
            else
               error('The HounsfieldToGreyscale() method was called, but this is not a CT image');
            end
        end

        function units_hu = GreyscaleToHounsfield(obj, units_greyscale)
            if isempty(obj.RescaleSlope) || isempty(obj.RescaleIntercept)
                units_hu = units_greyscale;
                return;
            end
            % Conversion from image raw intensity values to Hounsfield units:
            if obj.IsCT
                units_hu = int16(units_greyscale)*obj.RescaleSlope + obj.RescaleIntercept;
            else
               error('The HounsfieldToGreyscale() method was called, but this is not a CT image');
            end
        end

        function units_rescaled = RescaledToGreyscale(obj, units_rescaled)
            if obj.IsCT || obj.IsMR
                units_rescaled = obj.HounsfieldToGreyscale(units_rescaled);
            else
                error('The RescaledToGreyscale() method was called, but this is not a CT or MR image');
            end

        end

        function copy = Copy(obj)
            copy = PTKDicomImage(obj.RawImage, obj.RescaleSlope, obj.RescaleIntercept, obj.VoxelSize, obj.Modality, obj.StudyUid, obj.MetaHeader);
            metaclass = ?PTKDicomImage;
            property_list = metaclass.Properties;
            for i = 1 : length(property_list)
                property = property_list{i};
                if (~property.Dependent) && (~property.Constant)
                    copy.(property.Name) = obj.(property.Name);
                end
            end
        end

        function copy = BlankCopy(obj)
            copy = PTKDicomImage([], obj.RescaleSlope, obj.RescaleIntercept, obj.VoxelSize, obj.Modality, obj.StudyUid, obj.MetaHeader);
            metaclass = ?PTKDicomImage;
            property_list = metaclass.Properties;
            for i = 1 : length(property_list)
                property = property_list{i};
                if (~property.Dependent) && (~property.Constant) && (~strcmp(property.Name, 'RawImage'))
                    copy.(property.Name) = obj.(property.Name);
                end
            end
        end

        function is_equal = eq(obj, other)
            metaclass = ?PTKImage;
            property_list = metaclass.Properties;
            for i = 1 : length(property_list)
                property = property_list{i};
                if (~property.Dependent) && (~ismember(property.Name, obj.PropertiesToIgnoreOnComparison))
                    if ~isequal(other.(property.Name), obj.(property.Name))
                        is_equal = false;
                        return;
                    end
                end
            end
            is_equal = true;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function modX = convertToRelative(x,relativeParams)

bins = relativeParams.relativeBins;
while( any((bins(1:end-1)-bins(2:end))>0))
    ndx = find( (bins(1:end-1)-bins(2:end))>0);
    bins(ndx+1) = bins(ndx);
end

[~,modX] = histc(x,bins);
modX(modX>numel(bins))=numel(bins);

validX = ~isnan(modX) & (modX<length(bins)) & (modX>0);
extra = x(validX)-bins(modX(validX));
relExtra = extra./(bins(modX(validX)+1)-bins(modX(validX)));
modX(validX) = modX(validX) + relExtra;
modX = modX-1; % To make it go from 0 instead of 1

author: kristinbranson

repo: JAABA


text:
function determ = daub8_determinant ( n )

%*****************************************************************************80
%
%% DAUB8_DETERMINANT returns the determinant of the DAUB8 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = - 1.0;

  return
end


text:
function y = linearize( x )

% LINEARIZE    Linearize.
%    For real affine X, Y = X.
%    For convex X, Y is linear, and satisfies Y >= X.
%    For concave X, Y is linear, and satisfies Y <= X.
% This is used primarily within CVX functions to efficiently implement
% certain monotonic functions.

persistent P
if isempty( P ),
    P.map = cvx_remap( { 'r_affine' }, ...
        { 'convex' }, { 'concave' }, [2,3,4] );
    P.funcs = { [], @lin_affn, @lin_cnvx, @lin_cncv };
end
y = cvx_unary_op( P, x );

function y = lin_affn( x )
y = x;

function y = lin_cnvx( x ) %#ok
cvx_begin set
    variable y(size(x))
    x <= y; %#ok
cvx_end

function y = lin_cncv( x ) %#ok
cvx_begin set
    variable y(size(x))
    x >= y; %#ok
cvx_end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function legendre_poly_coef_test ( )

%*****************************************************************************80
%
%% LEGENDRE_POLY_COEF_TEST tests LEGENDRE_POLY_COEF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_POLY_COEF_TEST\n' );
  fprintf ( 1, '  LEGENDRE_POLY_COEF determines the Legendre \n' );
  fprintf ( 1, '  P(N) polynomial coefficients.\n' );

  c = legendre_poly_coef ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  P(%d)', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( j == 0 )
        fprintf ( 1, '  %f\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %f * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '  %f * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  return
end


text:
function [inv_map,bi_inv_map,logdet,invA] = invchol2(A)
% INVCHOL2
% Does a Cholesky decomposition on A and returns logdet, inverse and
% two function handles that respectively map X to A\X and A\X/A.
%

if nargin==0
    test_this();
    return;
end

if isreal(A)

    R = chol(A); %R'*R = A
    inv_map = @(X) R\(R'\X);

    % inv(A)*X*inv(A)
    bi_inv_map = @(X) (inv_map(X)/R)/(R');

    if nargout>2
        logdet = 2*sum(log(diag(R)));
        if nargout>3
            invA = inv_map(eye(size(A)));
        end
    end
else
    inv_map = @(X) A\X;

    % inv(A)*X*inv(A)
    bi_inv_map = @(X) (A\X)/A;

    if nargout>2
        logdet = log(det(A));
        if nargout>3
            invA = inv_map(eye(size(A)));
        end
    end
end

function test_this()
dim = 3;
r = randn(dim,2*dim);
A = r*r';
[inv_map,bi_inv_map,logdet,iA] = invchol2(A);
[logdet,log(det(A))]
X = randn(dim,3);
Y1 = A\X,
Y2 = inv_map(X)

Z1 = (A\X)/A
Z2 = bi_inv_map(X)

iA,
inv(A)

author: nesl


text:
function [] = LogScale( whichaxis,logbase,varargin)
%LogScale(whichaxis,logbase) renames the tick labels on your axes for a
%logarithmically scaled variable.
%
%INPUT
%   whichaxis   which axis did you log scale (can be 'x' 'y' or 'xy')
%   logbase     base of the logarithm you used for your points (ex 10)
%
%EXAMPLE
%   plot(log10(Xvals),Yvals)
%   LogScale('x',10)
%
%DLevenstein 2015
%%
p = inputParser;
addParameter(p,'exp',false)
addParameter(p,'nohalf',false)
parse(p,varargin{:})
expo = p.Results.exp;
nohalf = p.Results.nohalf;

%%
if strcmp(whichaxis,'y') || strcmp(whichaxis,'xy')
    range = get(gca,'YLim');
    switch nohalf
        case true
            range(1) = ceil(range(1)); range(2) = floor(range(2));
        case false
            range(1) = ceil(range(1).*2)./2; range(2) = floor(range(2).*2)./2;
    end
    ticks = [range(1):range(2)];
    if length(ticks)<=3 & ~nohalf
        ticks = [range(1):0.5:range(2)];
    end

    if length(ticks)>5
        ticks = ticks(1:2:end);
    end

    set(gca,'YTick',ticks)
    set(gca,'YTickLabels',round(logbase.^ticks,2,'significant'))

    if expo && logbase==10
        tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false);
        tickstrings = cellfun(@(X) replace(X,'','^'),tickstrings,'uniformoutput',false);
        tickstrings = cellfun(@(X) [num2str(logbase),X(1:end-1)],tickstrings,'uniformoutput',false);
        set(gca,'YTickLabels',tickstrings)
    end

end

if strcmp(whichaxis,'x') || strcmp(whichaxis,'xy')
    range = get(gca,'XLim');
    switch nohalf
        case true
            range(1) = ceil(range(1)); range(2) = floor(range(2));
        case false
            range(1) = ceil(range(1).*2)./2; range(2) = floor(range(2).*2)./2;
    end
    ticks = [range(1):range(2)];
    if length(ticks)<=3 & ~nohalf
        ticks = [range(1):0.5:range(2)];
    end

    if length(ticks)>5
        ticks = ticks(1:2:end);
    end

    set(gca,'XTick',ticks)
    set(gca,'XTickLabels',round(logbase.^ticks,2,'significant'))

    if expo && logbase==10
        tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false);
        tickstrings = cellfun(@(X) replace(X,'','^'),tickstrings,'uniformoutput',false);
        tickstrings = cellfun(@(X) [num2str(logbase),X(1:end-1)],tickstrings,'uniformoutput',false);
        set(gca,'XTickLabels',tickstrings)
    end
end

if strcmp(whichaxis,'z')
    range = get(gca,'ZLim');
    switch nohalf
        case true
            range(1) = ceil(range(1)); range(2) = floor(range(2));
        case false
            range(1) = ceil(range(1).*2)./2; range(2) = floor(range(2).*2)./2;
    end
    ticks = [range(1):range(2)];
    if length(ticks)<=3 & ~nohalf
        ticks = [range(1):0.5:range(2)];
    end

    if length(ticks)>5
        ticks = ticks(1:2:end);
    end

    set(gca,'ZTick',ticks)
    set(gca,'ZTickLabels',round(logbase.^ticks,2,'significant'))

    if expo
        tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false);
        tickstrings = cellfun(@(X) replace(X,'','^'),tickstrings,'uniformoutput',false);
        tickstrings = cellfun(@(X) [num2str(logbase),X(1:end-1)],tickstrings,'uniformoutput',false);
        set(gca,'ZTickLabels',tickstrings)
    end
end

if strcmp(whichaxis,'c')
    range = get(gca,'CLim');
    switch nohalf
        case true
            range(1) = ceil(range(1)); range(2) = floor(range(2));
        case false
            range(1) = ceil(range(1).*2)./2; range(2) = floor(range(2).*2)./2;
    end
    ticks = [range(1):range(2)];
    if length(ticks)<=3 & ~nohalf
        ticks = [range(1):0.5:range(2)];
    end

    if length(ticks)>=5
        ticks = ticks(1:2:end);
    end

    cb = get(gca,'colorbar');

    set(cb,'Ticks',ticks)
    set(cb,'TickLabels',round(logbase.^ticks,1,'significant'))

    if expo
        tickstrings = cellfun(@num2str,(num2cell(ticks)),'uniformoutput',false);
        tickstrings = cellfun(@(X) replace(X,'','^'),tickstrings,'uniformoutput',false);
        tickstrings = cellfun(@(X) [num2str(logbase),X(1:end-1)],tickstrings,'uniformoutput',false);
        set(cb,'TickLabels',tickstrings)
    end

end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [A_zV, a_zV, B_zV, b_zV] = GetPBConsts(DepthsV, EnergyC, PBDataS, flag);

%       - constantsAaBb(1) - depth
%       - constantsAaBb(2) - A
%       - constantsAaBb(3) - a
%       - constantsAaBb(4) - B
%       - constantsAaBb(5) - b
%
%     Example
%     using Ahnesjoe's notation
%     0.075  2.0958E-02  22.777100  1.8582E-05  0.058308
%     0.225  2.4948E-02  12.957500  2.1138E-05  0.067823
%
%LM:  JOD, 3 Nov 03, allow nearest neighbor indexing.
%JOD, 17 Nov 03, added scale factor, default at 100.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

scale = 100;
%For convenience, scale up A and B coefficients, so dose values are closer to unity.

if nargin == 4 & strcmp(flag,'nearest')

  Aahn6=PBDataS.aahn6b;
  Aahn18=flipud(PBDataS.aahn18b);

  MaxDepth=50;

  if EnergyC==6
   ParamMat=Aahn6;

  else
   if EnergyC==18
    ParamMat=Aahn18;
   end
  end

  DepthTableV=ParamMat(:,1);

  %our own linear interpolation (faster than Matlab).
  %compute two interpolation indices:
  if any([DepthsV>MaxDepth])
       disp('Warning!  Some depths exceed 50 cm, which is dosimetry limit!')
  end

  DepthsV=DepthsV.*[DepthsV<=MaxDepth]+MaxDepth.*[DepthsV>MaxDepth];

  IndexLow=floor((DepthsV+0.075)/0.15)+1;
  IndexHigh=ceil((DepthsV+0.075)/0.15)+1;

  DepthLow=DepthTableV(IndexLow);
  DepthHigh=DepthTableV(IndexHigh);

  indexV = round((DepthsV+0.075)/0.15)+1;

  A_zV= ParamMat(indexV,2);
  a_zV= ParamMat(indexV,3);
  B_zV= ParamMat(indexV,4);
  b_zV= ParamMat(indexV,5);

  A_zV = A_zV * scale;
  B_zV = B_zV * scale;

  return

end

%-----------code below not currently used (terrible!)---------------------%

Aahn6=PBDataS.aahn6b;
Aahn18=flipud(PBDataS.aahn18b);
Aahn8MCF=flipud(PBDataS.aahn8MMM);

MaxDepth=50;

if EnergyC==6
 ParamMat=Aahn6;
else
 if EnergyC==18
  ParamMat=Aahn18;
 end
end

if EnergyC==8
   ParamMat=Aahn8MCF;
   MaxDepth=20;
end

DepthTableV=ParamMat(:,1);

%our own linear interpolation (faster than Matlab).
%compute two interpolation indices:
if any([DepthsV>MaxDepth])
     warning('Some depths exceed 50 cm, which is dosimetry limit!')
 end

DepthsV=DepthsV.*[DepthsV<=MaxDepth]+MaxDepth.*[DepthsV>MaxDepth];

IndexLow=floor((DepthsV+0.075)/0.15)+1;
IndexHigh=ceil((DepthsV+0.075)/0.15)+1;

if EnergyC==8
    IndexLow=floor((DepthsV)/0.253) + 1;
    IndexHigh=ceil((DepthsV)/0.253) + 1;
    IndexLow(IndexLow > length(ParamMat(:,1))) = length(ParamMat(:,1));
    IndexHigh(IndexHigh > length(ParamMat(:,1))) = length(ParamMat(:,1));
end

DepthLow=DepthTableV(IndexLow);
DepthHigh=DepthTableV(IndexHigh);

A_zV=((ParamMat(IndexLow,2).*(DepthHigh-DepthsV)+...
    ParamMat(IndexHigh,2).*(DepthsV-DepthLow))./(DepthHigh-DepthLow+10*eps)).*...
      [DepthHigh~=DepthLow]+ParamMat(IndexLow,2).*[DepthHigh==DepthLow];
a_zV=((ParamMat(IndexLow,3).*(DepthHigh-DepthsV)+...
    ParamMat(IndexHigh,3).*(DepthsV-DepthLow))./(DepthHigh-DepthLow+10*eps)).*...
      [DepthHigh~=DepthLow]+ParamMat(IndexLow,3).*[DepthHigh==DepthLow];
B_zV=((ParamMat(IndexLow,4).*(DepthHigh-DepthsV)+...
   ParamMat(IndexHigh,4).*(DepthsV-DepthLow))./(DepthHigh-DepthLow+10*eps)).*...
      [DepthHigh~=DepthLow]+ParamMat(IndexLow,4).*[DepthHigh==DepthLow];
b_zV=((ParamMat(IndexLow,5).*(DepthHigh-DepthsV)+...
     ParamMat(IndexHigh,5).*(DepthsV-DepthLow))./(DepthHigh-DepthLow+10*eps)).*...
       [DepthHigh~=DepthLow]+ParamMat(IndexLow,5).*[DepthHigh==DepthLow];

%clear IndexLow IndexHigh DepthLow DepthHigh

%Scale final answers for convenience:

A_zV = A_zV * scale;
B_zV = B_zV * scale;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [Tshed,nShedsT] = BSLwatershed(regionPT,maskPT)
%"BSLwatershed"
%   BSL subfunction -- Returns watersheds for BSL estimate and the total
%   number of non-zero watersheds
%
% CRS 05/20/13
%
%Usage:
%   [Tshed,nShedsT] = BSLwatershed(regionPT,maskPT)
%       regionPT = subregion from PET images in SUV
%       maskPT   = VOI mask for restricting the watersheds
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.%
%
%% Smooth region to limit WS over-segmentation
smoothPT = fspecial('gaussian', 9, 3);
bPT = imfilter(regionPT,smoothPT,'replicate','same','conv');
bPT(bPT < 0) = 0;

%% Build 2D Watershed Transform of Region
Tshed = zeros(size(regionPT));
for i = 1:size(regionPT,3)

    % get slices
    PTtemp  = regionPT(:,:,i);
    bPTtemp = bPT(:,:,i);

    % get sheds
    % performs a 2D WS on an upsampled gradient image to remove boundary
    % elements
    hy = fspecial('sobel'); hx = hy';
    DelY = imfilter(double(PTtemp), hy, 'replicate');
    DelX = imfilter(double(PTtemp), hx, 'replicate');
    gradmagS = imresize(sqrt(DelX.^2 + DelY.^2),32*size(bPTtemp),'bicubic');
    Ls = double(watershed(gradmagS,8));
    % Down Samples the WS image to reduce the number of boundary elements
    L = imresize(Ls,size(bPTtemp),'nearest');

    %% Detemines the number of sheds for memory allocation
    Shed = []; ShedVal = []; meanShed = [];
    k = 0;
    for j = 1:max(L(:));
        if ( numel(find(maskPT(:,:,i).*L == j)) ~= 0 )
            k = k + 1;
        end
    end
    nSheds(i) = k;
    Shed = zeros([ size(PTtemp) nSheds(i)]);

    % reorders sheds from low to highest mean SUV and builds masks
    k = 0;
    for j = 1:max(L(:));
        if ( numel(find(maskPT(:,:,i).*L == j)) ~= 0 )
            k = k + 1;
            mShed = double(L);
            mShed(mShed ~= j) = 0;
            mShed(mShed > 0) = 1;
            Shed(:,:,k) = mShed.*PTtemp;
            ShedVal(k) = mean(nonzeros(Shed(:,:,k)));
            meanShed(:,:,k) = mShed*ShedVal(k);
        end
    end
    [~,Ival] = sort(ShedVal);
    Shed = Shed(:,:,Ival);
    meanShed = meanShed(:,:,Ival);

    % bulds a total shed from the means of each WS
    TmpShed = zeros(size(L));
    for j = 1:nSheds(i), TmpShed = meanShed(:,:,j) + TmpShed; end
    TmpShed1 = TmpShed;

    %% Unclassified voxel assignment
    % Finds unclassified voxels
    Cracks = double(L);
    Cracks(Cracks ~= 0) = 0.5; Cracks(Cracks == 0) = 1; Cracks(Cracks < 1)  = 0;
    iCrack = find(Cracks); Cracks = Cracks.*PTtemp;
    % Fill in unclassified voxels using the nearest mean shed value from
    % NEWS (k < 3) avoids infinite loop
    k = 0;
    while (sum(Cracks(:)) > 0  && k < 3)
        k = k + 1;
        [mSize, nSize] = size(L);
        for j = 1:numel(iCrack)
            q = floor( ( iCrack(j)-1 ) / mSize ) + 1;
            p = mod(iCrack(j),mSize);
            if (p == 0), p = mSize; end
            jT = p-1 + (q-1)*mSize; jB = p+1 + (q-1)*mSize;
            jL = p   + (q-2)*mSize; jR = p   + (q)*mSize;

            iVal = []; jCross = [];
            if (p > 1 && p < mSize)
                if (q > 1 && q < nSize), jCross = [ jT jB jL jR ];
                else
                    if (q > 1), jCross = [ jT jB jL ];
                    else jCross = [ jT jB jR ];
                    end
                end
            else
                if (p > 1)
                    if (q > 1 && q < nSize), jCross = [ jT jL jR ];
                    else
                        if (q > 1), jCross = [ jT jL ];
                        else jCross = [ jT jR ];
                        end
                    end
                else
                    if (q > 1), jCross = [ jB jL ];
                    else jCross = [ jB jR ];
                    end
                end
            end
            [~, iVal] = min(abs( TmpShed(jCross) - Cracks(iCrack(j)) ));
            if (TmpShed(jCross(iVal)) == 0)
                TmpShed1(iCrack(j)) = Cracks(iCrack(j));
            else
                TmpShed1(iCrack(j)) = TmpShed(jCross(iVal));
            end
        end
        TmpShed = TmpShed1;
        Cracks = TmpShed;
        Cracks(Cracks ~= 0) = 0.5; Cracks(Cracks == 0) = 1; Cracks(Cracks < 1)  = 0;
        iCrack = find(Cracks); Cracks = Cracks.*PTtemp;
    end
    Tshed(:,:,i) = TmpShed1;

end
nShedsT = sum(nSheds);


text:
function readBinary(obj,filename)
%READBINARY Load sequence from binary file.
%   READBINARY(seqObj, filename) Read the given filename and load sequence
%   data stored in the binary version of the Pulseq open file format. The
%   binary format is described in the specficiation available at
%   http://pulseq.github.io
%
%   Examples:
%   Load the sequence defined in gre.bin in sequences directory
%
%       readBinary(seqObj,'sequences/gre.bin')
%
% See also  writeBinary

binaryCodes = obj.getBinaryCodes();
fid=fopen(filename);
magicNum = fread(fid,8,'uchar');
assert(all(magicNum==binaryCodes.fileHeader(:)),'Not a binary file');
version_major = fread(fid,1,'int64');
version_minor = fread(fid,1,'int64');
version_revision = fread(fid,1,'int64');
assert(version_major==binaryCodes.version_major,'Unsupported version_major %d', version_major)
assert(version_minor==binaryCodes.version_minor,'Unsupported version_minor %d', version_minor)
assert(version_revision==binaryCodes.version_revision,'Unsupported version_revision %d', version_revision)

% set version
obj.version_major = version_major;
obj.version_minor = version_minor;
obj.version_revision = version_revision;

% Clear sequence data
obj.blockEvents={};
obj.definitions=containers.Map();
obj.gradLibrary=mr.EventLibrary();
obj.shapeLibrary=mr.EventLibrary();
obj.rfLibrary=mr.EventLibrary();
obj.adcLibrary=mr.EventLibrary();
obj.delayLibrary=mr.EventLibrary();

% Load data from file
while true
    section = int64(fread(fid,1,'int64'));
    if isempty(section)
        break
    end

    switch section
        case binaryCodes.section.definitions
            obj.definitions = readDefinitions(fid);
        case binaryCodes.section.blocks
            obj.blockEvents = readBlocks(fid);
        case binaryCodes.section.rf
            format = {'float64','int32','int32','int32','float64','float64'};
            obj.rfLibrary = readEvents(fid,format,[1 1 1 1e-6 1 1]);
        case binaryCodes.section.gradients
            format = {'float64','int32','int32'};
            obj.gradLibrary = readEvents(fid,format,[1 1 1e-6],'g',obj.gradLibrary);
        case binaryCodes.section.trapezoids
            format = {'float64','int64','int64','int64','int64'};
            obj.gradLibrary = readEvents(fid,format,[1 1e-6 1e-6 1e-6 1e-6],'t',obj.gradLibrary);
        case binaryCodes.section.adc
            format = {'int64','int64','int64','float64','float64'};
            obj.adcLibrary = readEvents(fid,format,[1 1e-9 1e-6 1 1]);
        case binaryCodes.section.delays
            format = {'int64'};
            obj.delayLibrary = readEvents(fid,format,1e-6);
        case binaryCodes.section.shapes
            obj.shapeLibrary = readShapes(fid);
        otherwise
            error('Unknown section code: %s',dec2hex(section));
    end
end
fclose(fid);

warning('readBinary() is not fully functional yet... as it needs to apply similar post-processing as read()');

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Helper functions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function def = readDefinitions(fid)
        %readDefinitions Read the [DEFINITIONS] section of a sequence file.
        %   defs=readDefinitions(fid) Read user definitions from file
        %   identifier of an open MR sequence file and return a map of
        %   key/value entries.

        def = containers.Map();
        numDefs = fread(fid,1,'int64');
        for i=1:numDefs
            c = fread(fid,1,'char');
            key=[];
            while c~=0 && length(key)<255
                key=[key c];
                c = fread(fid,1,'char');
            end
            key = char(key);
            count = fread(fid,1,'int8');
            type = fread(fid,1,'char');
            switch type
                case 'f'
                    values = double(fread(fid,count,'float64'));
                case 'i'
                    values = int(fread(fid,count,'int32'));
                case 'c'
                    values = char(fread(fid,count,'char')');
                    if values(end)==0
                        values=values(1:(end-1));
                    end
                otherwise
                    error('Unknown definition type');
            end
            def(key) = values;
        end
    end

    function eventTable = readBlocks(fid)
        %readBlocks Read the [BLOCKS] section of a sequence file.
        %   library=readBlocks(fid) Read blocks from file identifier of an
        %   open MR sequence file and return the event table.

        numBlocks = double(fread(fid,1,'int64'));
        eventIds = double(fread(fid,7*numBlocks,'int32'));
        eventTable_tmp = reshape(eventIds,7,numBlocks)';
        eventTable = cell(1, numBlocks);
        for ii = 1:numBlocks
            eventTable{ii} = eventTable_tmp(ii,:);
        end
    end

    function eventLibrary = readEvents(fid,format,scale,type,eventLibrary)
        %readEvents Read an event section of a sequence file.
        %   library=readEvents(fid,fmt,scale) Read event data of given
        %   format and scale elements according to column vector scale.
        %
        %   library=readEvents(fid,fmt,scale,type) Attach the type string
        %   to elements of the library.
        %
        %   library=readEvents(...,library) Append new events to the given
        %   library.

        if nargin<3
            scale=1;
        end
        if nargin<5
            eventLibrary=mr.EventLibrary();
        end
        numEvents = fread(fid,1,'int64');
        for i=1:numEvents
            id = double(fread(fid,1,'int32'));
            data = zeros(1,length(format));
            for j=1:length(format)
                data(j) = double(fread(fid,1,format{j}));
            end
            data=scale.*data;
            if nargin<4
                eventLibrary.insert(id,data);
            else
                eventLibrary.insert(id,data,type);
            end
        end
    end

    function shapeLibrary = readShapes(fid)
        %readShapes Read the [SHAPES] section of a sequence file.
        %   library=readShapes(fid) Read shapes from file identifier of an
        %   open MR sequence file and return a library of shapes.

        shapeLibrary=mr.EventLibrary();
        numShapes = fread(fid,1,'int64');
        for i=1:numShapes
            id = double(fread(fid,1,'int32'));
            numUncompressed = double(fread(fid,1,'int64'));
            numCompressed = double(fread(fid,1,'int64'));
            data = double(fread(fid,numCompressed,'float64'))';
            shapeData = [numUncompressed data];
            shapeLibrary.insert(id,shapeData);
        end
    end

end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function average = pointSetsAverage(pointSets, varargin)
%POINTSETSAVERAGE Compute the average of several point sets
%
%   AVERAGESET = pointSetsAverage(POINTSETS)
%   POINTSETS is a cell array containing several liste of points with the
%   same number of points. The function compute the average coordinate of
%   each vertex, and return the resulting average point set.
%
%   Example
%   pointSetsAverage
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-04-01,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% check input
if ~iscell(pointSets)
    error('First argument must be a cell array');
end

% number of sets
nSets   = length(pointSets);

% get reference size of coordinates array
set1    = pointSets{1};
refSize = size(set1);

% allocate memory for result
average = zeros(refSize);

% iterate on point sets
for i = 1:nSets
    % get current point set, and check its size
    set = pointSets{i};
    if sum(size(set) ~= refSize) > 0
        error('All point sets must have the same size');
    end

    % cumulative sum of coordinates
    average = average + set;
end

% normalize by the number of sets
average = average / nSets;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function SO3F = inv(SO3F)
% Define the inverse function $g$ of an SO3Fun $f$ by $g(R^{-1}) = f(R)$
% for all rotations $R\in SO(3)$.
%
% Syntax
%   SO3F = inv(F)
%
% Input
%  F - @SO3Fun
%
% Output
%  SO3F - @SO3FunHarmonic
%

SO3F = SO3FunHarmonic(SO3F);
SO3F = inv(SO3F);

% or alternative use:
% SO3F = SO3FunHandle(@(r) SO3F.eval(inv(r)),SO3F.SS,SO3F.CS);

end


text:
function y = FM(x,beta)

[r, c] = size(x);
if r*c == 0
    y = [];
    return;
end;
if (r == 1)
    x = x(:);
    len = c;
else
    len = r;
end;

Df = 2*pi*beta*100;

x = interp(x,20);
x = x*Df;

Fs = 20*8000;
Fc = 40000;
t = (0:1/Fs:((size(x,1)-1)/Fs))';
t = t(:, ones(1, size(x, 2)));

    x = 2 / Fs * pi * x;
    x = [zeros(1, size(x, 2)); cumsum(x(1:size(x,1)-1, :))];
    y = cos(2 * pi * Fc * t + x );
   plot(y)

author: Sable

repo: mcbench-benchmarks


text:
function h = mobile_text(varargin)
% MOBILE_TEXT('str1', 'str2', ...) places each string in a random position
% on the current axes.  The strings can be dragged around with the mouse.
% Returns handles to the text objects, for setting colors, etc.

strs = varargin;
n = length(strs);
ax = axis;
% random placement
x = rand(1,n)*(ax(2)-ax(1)) + ax(1);
y = rand(1,n)*(ax(4)-ax(3)) + ax(3);
h = [];
for i = 1:n
  h = [h text(x(i),y(i),strs{i})];
end
set(h,'ButtonDownFcn','move_obj(1)');
if nargout == 0
  clear h
end

author: andrewssobral

repo: lrslibrary


text:
%Author : Athi Narayanan S
%M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://sites.google.com/site/athisnarayanan/

function q5=mat_dec(data1,blkx)
[m,n]=size(data1);
r3=m/blkx;c3=n/blkx;q4=0;q1=0;
for i=1:r3
   for j=1:c3
      for s=1:blkx
         for k=1:blkx
            p3=s+q4;
            q2=k+q1;
            q5(s,k,i,j)=data1(p3,q2);
         end
      end
      q1=q1+blkx;
   end
   q4=q4+blkx;q1=0;
end


text:
function f = vl_ertr(f) ;
% VL_ERTR  Transpose exremal regions frames
%  F = VL_ERTR(F) transposes the frames F as returned by VL_MSER(). This
%  conversion is required as the VL_MSER algorithm considers the column
%  index I as the frist image index, while according standard image
%  convention the first coordinate is the abscissa X.
%
%  See also: VL_MSER(), VL_HELP().

if size(f,1) ~= 5
  error('F is not in the right format') ;
end

% adjust convention
f = f([2 1 5 4 3],:) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function lagrange2D
%% STARTVARIABLES
[X,Y]=meshgrid(-1:1:5,-1:1:5);
Z = sin(Y)*sin(X);

step = 0.2;
%% INTERPOLATION
% erst kurven in x richtung ausrechnen
xCurves={};
for i=1:size(X,1)
    % x vector must be a column vector for the lagrange function
    x = X(i,:)';
    % z vector must be a column vector for the lagrange function
    z = Z(i,:)';
    p=[];
    for j = x(1):step:x(end)
        % interpolate for every parameter value j and add it to p
        p = [p,lagrange(x,z,j)];
    end
    % save curves in x direction
    xCurves{i} = p;
end

y = Y(:,1);
% matrix for the graphical outpu
A=[];
% interpolate in y-direction
for i=1:length(xCurves{1})
    p=[];
    z=[];
    for l=1:length(y)
        z = [z;xCurves{l}(i)];
    end
    for j = y(1):step:y(end)
         % interpolate for every parameter value j and add it to p
        p = [p;lagrange(y,z,j)];
    end
    A = [A,p];
end

%% GRAPHICAL OUTPUT
% plot surface
surf(x(1):(x(end)-x(1))/(size(A,1)-1):x(end),y(1):(y(end)-y(1))/(size(A,1)-1):y(end),A);
hold on;
% plot points
for i=1:size(X,1)
    for j=1:size(Y,1)
        p = plot3(X(i,j),Y(i,j),Z(i,j));
        set(p,'Marker','.');
        set(p,'MarkerSize',30);
    end
end

end
%% nested functions
function res = lagrange(x,z,t)
% lagrange interpolation in 1d
% input: x...vektor, welcher die x (bzw. y) Werte beinhaltet

res = 0;
for i=1:length(x)
    % delete the component not needed
    temp = [x(1:i-1);x(i+1:end)];

    denominator = x(i)-temp;

    numerator = t-temp;
    res = res + z(i)*(prod(numerator))/(prod(denominator));
end
end


text:
function f = constructor(f, op, varargin)
%CONSTRUCTOR   CHEBFUN3 constructor.
%   Given a function OP of three variables, this code represents it as a
%   CHEBFUN3 object. A CHEBFUN3 object is a low-rank representation
%   expressing a function as a trilinear product of a discrete core tensor
%   and three quasimatrices consisting of univariate functions.
%
%
%   Since March 2023, chebfun3f is used by default to construct CHEBFUN3
%   objects. The legacy constructor chebfun3classic can be accessed by
%   providing the 'classic' flag as additional input argument.
%
% See also CHEBFUN2, CHEBFUN3T and CHEBFUN3V.

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

% Parse the inputs:
[op, dom, pref, fiberDim, vectorize, isEqui, fixedRank] = ...
    parseInputs(op, varargin{:});

% The 'equi' flag can be used only with numeric data:
if ( isEqui && ~isa(op, 'double') )
    error('CHEBFUN:CHEBFUN3:constructor:equi', ...
        'The EQUI flag is valid only when constructing from numeric data');
end

if ( isa(op, 'chebfun3') )     % CHEBFUN3( CHEBFUN3 )
    f = op;
    return
elseif ( isa(op, 'double') )   % CHEBFUN3( DOUBLE )
    f = chebfun3.chebfun3double(f, op, dom, pref, isEqui);
elseif ( strcmpi(pref.cheb3Prefs.constructor, 'classic') )
    f = chebfun3.chebfun3classic(f, op, pref, dom, vectorize, fiberDim);
else
    f = chebfun3.chebfun3f(f, op, pref, dom, vectorize);
end

if ( fixedRank )
    % Simplify the rank if requested.
    f = fixTheRank(f , fixedRank);
end
end
%% End of constructor

%%
function [op, dom, pref, fiberDim, vectorize, isEqui, ...
    fixedRank] = parseInputs(op, varargin)
vectorize = 0;
isEqui = 0;
isCoeffs = 0;
fixedRank = 0;
pref = chebfunpref();
fiberDim = [];
dom = [-1 1 -1 1 -1 1];

% Preferences structure given?
isPref = find(cellfun(@(p) isa(p, 'chebfunpref'), varargin));
if ( any(isPref) )
    pref = varargin{isPref};
    varargin(isPref) = [];
end

if ( isa(op, 'char') )     % CHEBFUN3( CHAR )
    op = str2op(op);
end

for k = 1:length(varargin)
    if any(strcmpi(varargin{k}, {'trig', 'periodic'}))
        pref.tech = @trigtech;
    elseif strcmpi(varargin{k}, 'eps')
        pref.cheb3Prefs.chebfun3eps = varargin{k+1};
    elseif strcmpi(varargin{k}, 'classic')
        pref.cheb3Prefs.constructor = 'classic';
    elseif strcmpi(varargin{k}, 'chebfun3f')
        pref.cheb3Prefs.constructor = 'chebfun3f';
    elseif strcmpi(varargin{k}, 'rank') % rank is specified.
        fixedRank = varargin{k+1};
    elseif ( isnumeric(varargin{k}) )
        if ( numel(varargin{k}) == 6 ) % domain is specified.
            dom = varargin{k};
        elseif ( numel(varargin{k}) == 3 ) % length is specified.
            if ( k > 1 && strcmpi(varargin{k-1}, 'rank') )
                % Rank is specified, not length. Don't confuse them.
                continue
            else
                % Interpret this as the user wants a fixed degree chebfun3
                % on the domain DOM.
                len = varargin{k};
                [xx, yy, zz] = chebpts3(len(1), len(2), len(3), dom);
                op = op(xx, yy, zz);
            end
        end
    elseif strcmpi(varargin{k}, {'fiberDim'})
        fiberDim = varargin{k+1};
    elseif any(strcmpi(varargin{k}, {'vectorize', 'vectorise'}))
        vectorize = true;
    elseif strcmpi(varargin{k}, 'coeffs')
        isCoeffs = 1;
    elseif strcmpi(varargin{k}, 'equi')
        isEqui = 1;
    end
end

if ( isCoeffs )
    op = chebfun3.coeffs2vals(op);
end

% If the vectorize flag is off, do we need to give user a warning?
if ( ~vectorize && ~isnumeric(op) ) % another check
    [vectorize, op] = vectorCheck(op, dom);
end

end

%%
function [vectorize, op] = vectorCheck(op, dom)
% Check for cases like op = @(x,y,z) x*y^2*z

vectorize = false;
[xx, yy, zz] = ndgrid(dom(1:2), dom(3:4), dom(5:6));
try
    A = feval(op, xx, yy, zz);
catch
    throwVectorWarning();
    vectorize = true;
    return
end

A = feval(op, xx, yy, zz);
if ( any(isinf(A(:) ) ) )
    error('CHEBFUN:CHEBFUN3:constructor:inf', ...
        'Function returned INF when evaluated');
elseif ( any(isnan(A(:)) ) )
    error('CHEBFUN:CHEBFUN3:constructor:nan', ...
        'Function returned NaN when evaluated');
end

if ( isscalar(A) )
    op = @(x,y,z) op(x,y,z) + 0*x + 0*y + 0*z;
end

end

%%
function throwVectorWarning()
warning('CHEBFUN:CHEBFUN3:constructor:vectorize',...
    ['Function did not correctly evaluate on an array.\n', ...
    'Turning on the ''vectorize'' flag. Did you intend this?\n', ...
    'Use the ''vectorize'' flag in the CHEBFUN3 constructor\n', ...
    'call to avoid this warning message.']);
end
%%
function op = str2op(op)
% OP = STR2OP(OP), finds independent variables in a string and returns an
% op handle than can be evaluated.

vars = symvar(op);        % Independent variables
numVars = numel(vars);
if ( numVars == 0 )
    op = @(x,y,z) eval (op);

elseif ( numVars == 1 )
    op = eval(['@(' vars{1} ', myVarBeta, myVarGamma)' op]);

elseif ( numVars == 2 )
    op = eval(['@(' vars{1} ',' vars{2} ', myVarGamma)' op]);

elseif ( numVars == 3 )
    op = eval(['@(' vars{1} ',' vars{2} ',' vars{3} ')' op]);

else
    error('CHEBFUN:CHEBFUN3:constructor:str2op:depvars', ...
        'Too many independent variables in string input.');
end

end

%%
function f = fixTheRank(f , fixedRank)
% Fix the rank of a CHEBFUN3. Used for calls to constructor with specified
% rank.

if ( any(fixedRank) < 0 )
    error('CHEBFUN:CHEBFUN3:constructor:fixTheRank:negative', ...
        'Ranks should all be nonnegative.')
elseif ( all(fixedRank) )
    [r1, r2, r3] = rank(f);
    t1 = fixedRank(1);
    t2 = fixedRank(2);
    t3 = fixedRank(3);

    % What to do with cols?
    if ( r1 > t1 )
        % Truncate cols:
        f.cols = f.cols(:, 1:t1);
        f.core = f.core(1:t1, :, :);
        r1 = t1; % New size of core
    elseif ( r1 < t1 )
        % Pad cols with approprate number of zero cols:
        zCols = chebfun(0, f.cols.domain);
        for jj = r1 : t1 - 1
            f.cols = [f.cols zCols];
        end
        % Pad mode 1 of the core tensor with zeros:
        tempCore = zeros(t1, r2, r3);
        tempCore(1:r1, :, :) = f.core;
        f.core = tempCore;
        r1 = t1; % New size of core
    end

    % What to do with rows?
    if ( r2 > t2 )
        % Truncate rows:
        f.rows = f.rows(:,1:t2);
        f.core = f.core(:, 1:t2, :);
        r2 = t2; % New size of core
    elseif ( r2 < t2 )
        % Pad rows with approprate number of zero rows:
        zRows = chebfun(0, f.rows.domain);
        for jj = r2 : t2 - 1
            f.rows = [f.rows zRows];
        end
        % Pad mode 2 of the core tensor with zeros:
        tempCore = zeros(r1, t2, r3);
        tempCore(:, 1:r2, :) = f.core;
        f.core = tempCore;
        r2 = t2; % New size of core
    end

    % What to do with tubes?
    if ( r3 > t3 )
        % Truncate tubes:
        f.tubes = f.tubes(:, 1:t3);
        f.core = f.core(:, :, 1:t3);
    elseif ( r3 < t3 )
        % Pad tubes with approprate number of zero tubes:
        zTubes = chebfun(0, f.tubes.domain);
        for jj = r3 : t3 - 1
            f.tubes = [f.tubes zTubes];
        end
        % Pad mode 3 of the core tensor with zeros:
        tempCore = zeros(r1, r2, t3);
        tempCore(:, :, 1:r3) = f.core;
        f.core = tempCore;
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function areBadVals = AreBadValues(A)
% Only returns true for bad array elements.
% Bad array elements include Inf (infinity), NA (missing), and NaN (not a
% number).
%
%    usage: areBadVals = AreBadValues(A)

    areBadVals = ~isfinite(A);
    % Inf, NA, and NaN are non-finite values.

end