text:
%

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

function [ac, ... % the aerodynamic center of the aircraft in meters measured from the Tornado frame of reference origin.
    staticMarginMAC, ... % the static margin of the aircraft with the given center of gravity as a factor of the mean chord.
    staticMarginInMeters]...
    = computeStaticMargin(...
nameTornadoAircraftGeometry, ...
estimatedVelocity, ...
airDensity, ...
tornadoDirectory, ...
centerOfGravity)

% Define an initial Tornado state with given estimated velocity.
state.AS=      estimatedVelocity; %Airspeed m/s
state.alpha=   0.0;       %Angle of attack, radians
state.betha=   0;         %Angle of sideslip, radians
state.P=       0;         %Rollrate, rad/s
state.Q=       0;         %pitchrate, rad/s
state.R=       0;         %yawrate, rad/s
state.adot=    0;         %Alpha time derivative rad/s
state.bdot=    0;         %Betha time derivative rad/s
state.ALT=     0;         %Altitude, m
state.rho=     airDensity;%Desity, kg/m^3
state.pgcorr=  0;         %Apply prandtl glauert compressibility correction

% Get some properties from the aircraft geometry
cd(tornadoDirectory);
settings=config('startup');
% Load the aircraft
cd(settings.acdir);
load(nameTornadoAircraftGeometry); % aircraft geometry is now in variable 'geo'.
% Overwrite center of gravity with given value (!)
geo.CG = [centerOfGravity 0 0];
% Overwrite reference point for moment computation with center of gravity.
geo.ref_point = [centerOfGravity 0 0];

% Find static margin
% Since aerodynamic center for entire aircraft should correspond with
% neutral point, the static margin should add up (cg + margin = neutral
% point = aerodynamic center entire plane).
cd(tornadoDirectory);
out=fFindstaticmargin(geo,state);
staticMarginMAC = out.h(1);
mean_chord = geo.c(1);
staticMarginInMeters = out.h(1)*mean_chord;
ac = out.ac(1);

cd(fileparts(which(mfilename))); % Back to the directory of this script.

end

author: jrgenerative

repo: fixed-wing-sim


text:
function sendmail = spm_cfg_sendmail
% SPM Configuration file for sendmail
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_sendmail.m 6929 2016-11-14 13:07:31Z guillaume $

% ---------------------------------------------------------------------
% Recipient
% ---------------------------------------------------------------------
recipient         = cfg_entry;
recipient.tag     = 'recipient';
recipient.name    = 'Recipient';
recipient.help    = {'User to receive mail.'};
recipient.strtype = 's';
recipient.num     = [1 Inf];
% ---------------------------------------------------------------------
% Subject
% ---------------------------------------------------------------------
subject         = cfg_entry;
subject.tag     = 'subject';
subject.name    = 'Subject';
subject.val     = {['[SPM] [%DATE%] On behalf of ' spm('Ver')]};
subject.help    = {'The subject line of the message. %DATE% will be replaced by a string containing the time and date when the email is sent.'};
subject.strtype = 's';
subject.num     = [1 Inf];
% ---------------------------------------------------------------------
% Message
% ---------------------------------------------------------------------
message         = cfg_entry;
message.tag     = 'message';
message.name    = 'Message';
message.val     = {'Hello from SPM!'};
message.help    = {'A string containing the message to send.'};
message.strtype = 's';
message.num     = [1 Inf];
% ---------------------------------------------------------------------
% Attachments
% ---------------------------------------------------------------------
attachments         = cfg_files;
attachments.tag     = 'attachments';
attachments.name    = 'Attachments';
attachments.val     = {{}};
attachments.help    = {'List of files to attach and send along with the message.'};
attachments.filter  = '.*';
attachments.ufilter = '.*';
attachments.num     = [0 Inf];
% ---------------------------------------------------------------------
% SMTP Server
% ---------------------------------------------------------------------
smtp         = cfg_entry;
smtp.tag     = 'smtp';
smtp.name    = 'SMTP Server';
smtp.help    = {'Your SMTP server. If not specified, look for sendmail help.'};
smtp.strtype = 's';
try
    smtp.val = {getpref('Internet','SMTP_Server')};
end
smtp.num     = [1 Inf];
% ---------------------------------------------------------------------
% E-mail
% ---------------------------------------------------------------------
email         = cfg_entry;
email.tag     = 'email';
email.name    = 'E-mail';
email.help    = {'Your e-mail address. Look in sendmail help how to store it.'};
email.strtype = 's';
try
    email.val = {getpref('Internet','E_mail')};
end
email.num     = [1 Inf];
% ---------------------------------------------------------------------
% Zip attachments
% ---------------------------------------------------------------------
zip         = cfg_menu;
zip.tag     = 'zip';
zip.name    = 'Zip attachments';
zip.val     = {'No'};
zip.help    = {'Zip attachments before being sent along with the message.'};
zip.labels  = {'Yes' 'No'}';
zip.values  = {'Yes' 'No'}';
% ---------------------------------------------------------------------
% Parameters
% ---------------------------------------------------------------------
params         = cfg_branch;
params.tag     = 'params';
params.name    = 'Parameters';
params.val     = { smtp email zip};
params.help    = {'Preferences for your e-mail server (Internet SMTP server) and your e-mail address. MATLAB tries to read the SMTP mail server from your system registry. This should work flawlessly. If you encounter any error, identify the outgoing mail server for your electronic mail application, which is usually listed in the application''s preferences, or, consult your e-mail system administrator, and update the parameters. Note that this function does not support e-mail servers that require authentication.'};
% ---------------------------------------------------------------------
% Sendmail
% ---------------------------------------------------------------------
sendmail       = cfg_exbranch;
sendmail.tag   = 'sendmail';
sendmail.name  = 'Sendmail';
sendmail.val   = { recipient subject message attachments params};
sendmail.help  = {'Send a mail message (attachments optionals) to an address.'};
sendmail.prog  = @spm_sendmail;
%_______________________________________________________________________

%_______________________________________________________________________
function spm_sendmail(job)

try
    setpref('Internet','SMTP_Server',job.params.smtp);
    setpref('Internet','E_mail',job.params.email);

    subj = strrep(job.subject,'%DATE%',datestr(now));
    mesg = strrep(job.message,'%DATE%',datestr(now));
    mesg = [mesg 10 10 '-- ' 10 10 'Statistical Parametric Mapping'];

    if ~isempty(job.attachments)
        if strcmpi(job.params.zip,'Yes')
            zipfile = fullfile(tempdir,'spm_sendmail.zip');
            zip(zipfile,job.attachments);
            job.attachments = {zipfile};
        end
        sendmail(job.recipient,subj,mesg,job.attachments);
    else
        sendmail(job.recipient,subj,mesg);
    end
catch
    %- not an error to prevent an analysis to crash because of just that...
    fprintf('Sendmail failed...\n');
end

author: spm


text:
%Script to compile mex files stored in "private" folder
%Written by: Prateek Jain (pjain@cs.utexas.edu) and Raghu Meka (raghu@cs.utexas.edu)
%Last updated: October, 2009

cd private;
%mex -O -lmwlapack -lmwblas bdsqr_mex.c dbdqr.f -UDBDQR_IN_C -output bdsqr
%mex -O -lmwlapack -lmwblas reorth_mex.c reorth.f -output reorth
%mex -O -lmwlapack -lmwblas tqlb_mex.c tqlb.f -output tqlb
mex compute_X_Omega.c
mex compOmegaYx.c
mex compOmegaYtx.c
cd ..

author: andrewssobral

repo: lrslibrary


text:
function f2 = p02_f2 ( x )

%*****************************************************************************80
%
%% P02_F2 evaluates the second derivative for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    LE Scales,
%    Introduction to Non-Linear Optimization,
%    Springer, 1985.
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  f2 = 2.0 + exp ( -x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pdf = cosine_pdf ( x, a, b )

%*****************************************************************************80
%
%% COSINE_PDF evaluates the Cosine PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = ( 1 / ( 2 * PI * B ) ) * COS ( ( X - A ) / B )
%    for A - PI * B <= X <= A + PI * B
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a - pi * b )

    pdf = 0.0;

  elseif ( x <= a + pi * b )

    y = ( x - a ) / b;

    pdf = 1.0 / ( 2.0 * pi * b ) * cos ( y );

  elseif ( a + pi * b < x )

    pdf = 0.0;

  end

  return
end


text:
function Dij = TS_PairwiseDist(tsOrOps,whatData,distanceMetric,doSave)
% TS_PairwiseDist   Pairwise distances between either TimeSeries or Operations
%
% (If doSave), stores pairwise distances back the data back in the local .mat
% file. Warning: this can be huge for large datasets.
%
%---INPUTS:
%
% tsOrOps: Compute pairwise distances between all pairs of time series
%               ('ts', default), or operations ('ops')
%
% whatData: what hctsa data source to use (cf. TS_LoadData)
%
% distanceMetric: what distance metric to use, e.g., 'euclidean' (default for
%                 time series) or 'corr_fast' (default for operations)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check Inputs, set defaults:
% ------------------------------------------------------------------------------

if nargin < 1 || isempty(tsOrOps)
    tsOrOps = 'ts';
    fprintf(1,'Time series by default.\n');
end
if ~ismember(tsOrOps,{'ts','ops'})
    error('Invalid specifier ''%s'' (should be ''ts'' or ''ops'')',tsOrOps);
end

if nargin < 2
    whatData = 'norm';
end

if nargin < 3
    switch tsOrOps
    case 'ts'
        distanceMetric = 'euclidean';
    case 'ops'
        distanceMetric = 'corr_fast';
    end
end

% doSave
if nargin < 4
    % By default, save back to the file specified as whatData
    doSave = 1;
end

% ------------------------------------------------------------------------------
% Load data
% ------------------------------------------------------------------------------

if ischar(whatData)
    [TS_DataMat,~,~,theFile] = TS_LoadData(whatData);
else
    % If provided a matrix, assume it's from HCTSA_N.mat by default:
    TS_DataMat = whatData;
    theFile = 'HCTSA_N.mat';
end

% ------------------------------------------------------------------------------
% Compute pairwise distances between all objects
% ------------------------------------------------------------------------------
if strcmp(tsOrOps,'ops')
    % Transpose if we want columns ('ops')
    TS_DataMat = TS_DataMat';
end

fprintf(1,'Computing %s distances between %u objects...',distanceMetric,size(TS_DataMat,1));
Dij = BF_pdist(TS_DataMat,distanceMetric,1,[],1);
fprintf(1,' Done.\n');

% ------------------------------------------------------------------------------
% Save back to the local hctsa .mat file
% ------------------------------------------------------------------------------
if doSave
    switch tsOrOps
    case 'ts'
        theWhat = 'time series';
    case 'ops'
        theWhat = 'operations';
    end
    fprintf(1,'Saving pairwise distance information for %s back to %s...',...
                        theWhat,theFile);

    switch tsOrOps
    case 'ts'
        ts_clust = struct('Dij',Dij,'distanceMetric',distanceMetric);
        save(theFile,'ts_clust','-append')
    case 'ops'
        op_clust = struct('Dij',Dij,'distanceMetric',distanceMetric);
        save(theFile,'op_clust','-append')
    end

    fprintf(1,' Done.\n');
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

    Update = 0;
    for i=1:m,
        Update = Update + alpha/m * (theta' * X(i,:)' - y(i)) *  X(i,:)';
    end

    theta = theta - Update;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

repo: machine-learning-coursera


text:
function [AccountVal, PayoutVal, FeeVal] = calcValuePayoutAndFees_orig( ...
    Holdings_0, assetPrice, annualWithdrawalRate, annualFee)
% CALCVALUEPAYOUTANDFEES_ORIG estimates the account value, the payout
% amounts, and the fees collected for a Monte Carlo simulation of the GMWB.
%
% Inputs:
% Holdings_0: An nAssets-by-1 vector of initial portfolio holdings, in
%   shares
% assetPrice: An nYears-by-nAssets-by-nTrials matrix of simulated asset
%   prices
% annualWithdrawalRate: A scalar indicating the fraction of the initial
%   portfolio value that the customer is guaranteed to withdraw each year
% annualFee: A scalar indicating the fraction of the current portfolio
%   value that the insurer may assess each year
%
% Outputs: All have size nYears+1-by-nTrials.  The first row is for the
%   initial year (t = 0)
% AccountVal: The total portofolio value after withdrawals.
% PayoutVal: The amount paid by the insurer to the customer if there is a
%   shortfall
% FeeVal: The amount collected by the insurer if funds remain in the
%   account

[nYears, ~, nTrials] = size(assetPrice);

AccountVal  = nan(nYears, nTrials);
PayoutVal   = nan(nYears, nTrials);
FeeVal      = nan(nYears, nTrials);

for iTrial = 1:nTrials
    % Initialize values at time 0.
    currentHoldings      = Holdings_0;
    currentPrice         = assetPrice(1, :, iTrial);
    AccountVal(1,iTrial) = currentPrice * currentHoldings;
    PayoutVal(1, iTrial) = 0;
    FeeVal(   1, iTrial) = 0;
    guaranteedWithdrawal = currentPrice * currentHoldings * annualWithdrawalRate;

    for iYear = 2:nYears
        currentPrice    = assetPrice(iYear, :, iTrial);
        curPortVal      = currentPrice * currentHoldings;
        totalWithdrawal = curPortVal * annualFee + guaranteedWithdrawal;

        FeeVal(iYear, iTrial) = curPortVal * annualFee;

        [currentHoldings, shortfall] = updateHoldings(currentPrice, ...
            currentHoldings, totalWithdrawal);

        PayoutVal(iYear, iTrial) = shortfall;
        AccountVal(iYear, iTrial) = currentPrice * currentHoldings;
    end
end
end

function [holdings, shortfall] = updateHoldings(currentPrice, ...
    currentHoldings, totalWithdrawal)
% Adjusts the current holdings and calculates the shortfall given the
% pre-withdrawal holdings, the withdrawal, and the current asset prices.

holdings = currentHoldings;
currentPrice = currentPrice';
curPortVal = currentHoldings .* currentPrice;

if any(curPortVal > 0)
    % We still have some funds in the portfolio: apply withdrawals in the
    % best way possible by distributing withdrawals across the holdings
    % proportionately to their values:
    newPortVal = curPortVal - curPortVal.*totalWithdrawal/sum(curPortVal);

    if sum(newPortVal) > 0      % Still have fund left in the account
        shortfall = 0;
        holdings = newPortVal ./ currentPrice;
        holdings(~isfinite(holdings)) = 0;
    else
        % We've just gone bankrupt.  Set all holdings to 0 and pay out the
        % shortfall.
        shortfall = -sum(newPortVal);
        holdings(:) = 0;
    end

else
    % Every asset has zero value-- the calculation is straightforward.
    shortfall   = totalWithdrawal;
    holdings(:) = 0;
end

end

author: Sable

repo: mcbench-benchmarks


text:
function cmap = CorrPlot(coeffs,isPlotText,ylabels)
im = coeffs;
im(isnan(im)) = 0;

% red-white-blue colormap
cmap = zeros(64,3);
cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)];
cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)];
cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)];
minlim = -1; %min(min(im));
maxlim = 1; %max(max(im));

RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap

image(RGB); axis equal; axis tight;
set(gca,'XTick',1:length(im),'YTick',1:length(im));
set(gca,'TickLength',[0,0]);

if exist('isPlotText','var'),
    if isPlotText,
        for i = 1:size(im,1),
            for j = 1:size(im,2),
                if im(j,i)~=0
                    text(i-0.3, j, num2str(round(im(j,i)*100)/100));%, 'Units', 'data') % flipped i/j necessary
                end
            end
        end
    end
end
if exist('ylabels','var'),
    set(gca,'YTickLabel',ylabels);
end
end

function RGB = ImageToRGB(im,cmap,minlim,maxlim)
L = size(cmap,1);
ix = round(interp1(linspace(minlim,maxlim,L),1:L,im,'linear','extrap'));
RGB = reshape(cmap(ix,:),[size(ix) 3]); % Make RGB image from scaled.
end

author: xiuyechen

repo: FishExplorer


text:
function [ node_num, tet_num ] = tet_mesh_order4_example_size ( )

%*****************************************************************************80
%
%% TET_MESH_ORDER4_EXAMPLE_SIZE sizes an example linear tet mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer TET_NUM, the number of tetrahedrons.
%
  node_num = 63;
  tet_num = 144;

  return
end


text:
% [readframe,nframes,fid,headerinfo] = get_readframe_fcn(filename)
% get_readframe_fcn inputs the name of a video, opens this video, and creates a
% function that can be used for random access to frames in the video. it can
% input files of type fmf, sbfmf, and avi. if videoio is installed on your
% machine, then get_readframe_fcn will use videoio to read all types of avi
% supported by videoio. otherwise, it uses VideoReader to open any type of avi
% supported by VideoReader. Run "help VideoReader" to determine what types of avi
% files are readable through VideoReader on your OS.
%
% Input:
% filename is the name of the video to read. The type of video is determined
% based on the extension, so be sure the extension of your movie is correct.
%
% Outputs:
% readframe is the handle to a function which inputs a frame number, reads this
% frame from filename, and outputs the corresponding frame.
% nframes is the number of frames in the video. For avis, this may be
% approximate, as it is computed from the duration and framerate.
% fid is the file identifier for filename. get_readframe_fcn will open the file
% filename. you must close fid yourself when you are done using the returned
% readframe.
% headerinfo is a struct with fields depending on the type of movie

function [readframe,nframes,fid,headerinfo] = get_readframe_fcn(filename,varargin)

i = find(strcmpi(varargin(1:2:end),'preload'));
if ~isempty(i),
  preload = varargin{i+1};
  varargin(i:i+1) = [];
else
  preload = false;
end
didpreload = false;

% allow videoio library to be used if it is installed and on the path
%CTRAX_ISVIDEOIO = exist('videoReader','file');
CTRAX_ISVIDEOIO = false;

if iscell(filename),

  readframes = cell(size(filename));
  nframes = inf;
  fid = nan(size(filename));
  headerinfo = cell(size(filename));
  for i = 1:numel(filename),
    [readframes{i},nframescurr,fid(i),headerinfo{i}] = get_readframe_fcn(filename{i},varargin{:});
    nframes = min(nframes,nframescurr);
  end

  readframe = @(f) multi_read_frame(f,readframes);
  return;

end

[~,ext] = splitext(filename);

if ispc && ~exist(filename,'file') && ~strcmpi(ext,'.seq'),
  [actualfilename,didfind] = GetPCShortcutFileActualPath(filename);
  if didfind,
    filename = actualfilename;
  end
end

if strcmpi(ext,'.fmf'),
  [header_size,version,nr,nc,bytes_per_chunk,nframes,data_format] = fmf_read_header(filename);
  fid = fopen(filename);
  readframe = @(f) fmfreadframe(fid,header_size+(f-1)*bytes_per_chunk,nr,nc,bytes_per_chunk,data_format);
  headerinfo = struct('header_size',header_size,'version',version,'nr',nr,'nc',nc,...
    'bytes_per_chunk',bytes_per_chunk,'nframes',nframes,'data_format',data_format,'type','fmf');
elseif strcmpi(ext,'.sbfmf'),
  [nr,nc,nframes,bgcenter,bgstd,frame2file] = sbfmf_read_header(filename);
  fid = fopen(filename);
  readframe = @(f) sbfmfreadframe(f,fid,frame2file,bgcenter);
  headerinfo = struct('nr',nr,'nc',nc,'nframes',nframes,'bgcenter',bgcenter,...
    'bgstd',bgstd,'frame2file',frame2file,'type','sbfmf');
elseif strcmpi(ext,'.ufmf'),
  headerinfo = ufmf_read_header(filename);
  readframe = ufmf_get_readframe_fcn(headerinfo,varargin{:});
  nframes = headerinfo.nframes;
  fid = headerinfo.fid;
elseif strcmpi(ext,'.mmf'),
  headerinfo = mmf_read_header(filename);
  readframe = mmf_get_readframe_fcn(headerinfo,varargin{:});
  nframes = headerinfo.nframes;
  fid = headerinfo.fid;
elseif strcmpi(ext,'.tif') || strcmp(ext,'.tiff'),
  info = imfinfo(filename);
  isimseq = false;
  if numel(info) == 1,
    filespec = regexprep(filename,['_\d+\.',ext,'$'],['_*.',ext]);
    imfiles = mydir(filespec);
    if numel(imfiles) > 1,
      imfiles = sort(imfiles);
      im = imread(imfiles{1});
      headerinfo = struct('nr',size(im,1),'nc',size(im,2),'ncolors',size(im,3),'nframes',numel(imfiles),...
        'type','imseq','imfiles',{imfiles});
      readframe = @(f) imseq_read_frame(f,imfiles);
      nframes = headerinfo.nframes;
      fid = -1;
      isimseq = true;
    end
  end
  if ~isimseq,
    headerinfo = struct('nr',info(1).Height,'nc',info(1).Width,'nframes',numel(info),'type','tif',...
      'bitdepth',info(1).BitDepth);
    readframe = @(f) deal(imread(filename,f),f);
    nframes = headerinfo.nframes;
    fid = -1;
  end
elseif strcmpi(ext,'.png'),
  info = imfinfo(filename);
  isimseq = false;
  filespec = regexprep(filename,'_\d+\.png$','_*.png');
  imfiles = mydir(filespec);
  if numel(imfiles) > 1,
    imfiles = sort(imfiles);
    im = imread(imfiles{1});
    headerinfo = struct('nr',size(im,1),'nc',size(im,2),'ncolors',size(im,3),'nframes',numel(imfiles),...
      'type','imseq','imfiles',{imfiles});
    readframe = @(f) imseq_read_frame(f,imfiles);
    nframes = headerinfo.nframes;
    fid = -1;
    isimseq = true;
  end
elseif strcmpi(ext,'.klb'),
  tfKLBLib = exist('readKLBslice','file')>0;
  if ~tfKLBLib
    error('get_readframe_fcn:klb','Cannot find readKLBslice function. Make sure KLB matlab library is installed.');
  end

  dim = myparse(varargin,'dim',3);
  headerinfo = readKLBheader(filename);
  readframe = klb_get_readframe_fcn(filename,varargin{:});
  nframes = headerinfo.xyzct(dim); % AL 201507: KLB doc a little unclear here but this seems right
  fid = 0;
elseif strcmpi(ext,'.mat'),
  videofiletype = load(filename,'videofiletype');
  switch videofiletype,

    case 'SingleLarvaTracker',
      videodata = load(filename);
      readframe = @(f) ReadSingleLarvaTrackerFrame(f,videodata.firstframeim,videodata.imraw,videodata.finalbbox,videodata.fps,varargin{:});
      nframes = numel(imraw);
      fid = 0;
      [nr,nc,~] = size(firstframeim);
      headerinfo = struct('nr',nr,'nc',nc,'nframes',nframes,'bgcenter',firstframeim,...
        'type','SingleLarvaTracker');
    otherwise
      error('Do not know how to parse mat file of type %s',videofiletype);
  end

elseif strcmpi(ext,'.seq'),
  [indexfilename,seqtype] = myparse(varargin,'indexfilename',0,'seqtype',0);

  if ischar(indexfilename),
    seqtype = 'shay';
  end

  % get file names
  if ispc && ~exist(filename,'file'),

    [actualfilename,didfind] = GetPCShortcutFileActualPath(filename);
    if ~didfind,
      error('Could not find movie file %s',filename);
    end
    filename0 = filename;
    % use actualfilename instead
    filename = actualfilename;

    % try the index file using the original filename
    if ~ischar(indexfilename),
      indexfilename = regexprep(filename0,'seq$','mat');
      [indexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
      if ~didfind && ~strcmp(filename0,filename),
        % try to find the index file using the source filename
        indexfilename = regexprep(filename,'seq$','mat');
        [indexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
        % if didn't find index file, then maybe this is a piotr seq file
        if ~didfind,
          if ischar(seqtype) && strcmpi(seqtype,'shay'),
            error('Could not find index file for %s',filename);
          else
            seqtype = 'piotr';
          end
        end
      end
    end

  else

    % set seqtype to piotr if indexfilename doesn't exist
    if ~ischar(seqtype) && ~ischar(indexfilename),
      indexfilename = regexprep(filename,'seq$','mat');
      if ~exist(indexfilename,'file'),
        seqtype = 'piotr';
      end
    end

  end

  % get actual filename for shortcuts
  if strcmpi(seqtype,'piotr'),
    warning('Closing file currently not implemented for Piotr''s seq files...');
    headerinfo = seqIo(filename,'getinfo');
    headerinfo.nr = headerinfo.height;
    headerinfo.nc = headerinfo.width;
    headerinfo.nframes = headerinfo.numFrames;
    headerinfo.type = 'seq_piotr';
    fid = 0;
    nframes = headerinfo.numFrames;
    sr = seqIo(filename,'r');
    readframe = @(f) seq_read_frame_piotr(f,sr);
  else

    if ischar(indexfilename),
      headerinfo = r_readseqinfo(filename,indexfilename);
    else
      headerinfo = r_readseqinfo(filename);
    end
    if numel(headerinfo.m_aiSeekPos) < headerinfo.m_iNumFrames,
      headerinfo.m_iNumFrames = numel(headerinfo.m_aiSeekPos);
    end
    headerinfo.nr = headerinfo.m_iHeight;
    headerinfo.nc = headerinfo.m_iWidth;
    headerinfo.nframes = headerinfo.m_iNumFrames;
    headerinfo.type = 'seq';
    fid = 0;
    nframes = headerinfo.m_iNumFrames;
    readframe = @(f) read_seq_frame(headerinfo,f);
  end

elseif strcmpi(ext,'.h5'),

  [dims,datasetname,fixdims,fixidx] = myparse(varargin,'slicedims',3,'datasetname',[],'fixdims',[],'fixidx',[]);
  if isempty(datasetname),
    headerinfo = h5info(filename);
    datasetname = headerinfo.Name;
    headerinfo = headerinfo.Datasets;
  else
    headerinfo = h5info(filename,datasetname);
  end
  headerinfo.datasetname = datasetname;
  headerinfo.dims = dims;
  headerinfo.fixdims = fixdims;
  headerinfo.fixidx = fixidx;
  headerinfo.sz = headerinfo.Dataspace.Size;
  headerinfo.ndims = numel(headerinfo.sz);
  headerinfo.nframes = headerinfo.sz(headerinfo.dims);

  readframe = h5_get_readframe_fcn(filename,headerinfo);
  nframes = headerinfo.nframes;
  fid = 0;

else
  fid = 0;

  isindexedmjpg = false;
  if strcmpi(ext,'.mjpg'),
    [indexfilename] = myparse(varargin,'indexfilename',[]);
    if isempty(indexfilename),
      % see if this looks like an indexed mjpg file
      [p,n] = fileparts(filename);
      indexfilename = fullfile(p,[n,'.txt']);
      if ~exist(indexfilename,'file') && ispc,
        [actualindexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
        if didfind,
          indexfilename = actualindexfilename;
        end
      end
      if exist(indexfilename,'file'),
        try
          ReadIndexedMJPGHeader(filename,indexfilename);
          isindexedmjpg = true;
        catch
          indexfilename = 0;
        end
      else
        indexfilename = 0;
      end
    elseif ischar(indexfilename),
      isindexedmjpg = true;
    end
  end

  if isindexedmjpg,
    % get file names
    if ispc && ~exist(filename,'file'),
      [actualfilename,didfind] = GetPCShortcutFileActualPath(filename);
      if ~didfind,
        error('Could not find movie file %s',filename);
      end
      filename0 = filename;
      % use actualfilename instead
      filename = actualfilename;

      [actualindexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
      if ~didfind,
        error('Could not find index file %s',indexfilename);
      end
      indexfilename0 = indexfilename;
      % use actualfilename instead
      indexfilename = actualindexfilename;
    end

    headerinfo = ReadIndexedMJPGHeader(filename,indexfilename);
    headerinfo.fid = 0;
    nframes = headerinfo.nframes;
    readframe = @(f) read_mjpg_frame(headerinfo,f);

  else

    if CTRAX_ISVIDEOIO,
      readerobj = videoReader(filename,'preciseFrames',30,'frameTimeoutMS',5000);
      info = getinfo(readerobj);
      nframes = info.numFrames;
      seek(readerobj,0);
      seek(readerobj,1);
      readframe = @(f) videoioreadframe(readerobj,f);
      headerinfo = info;
      headerinfo.type = 'avi';
    else

      % optional argument to not use the read function of VideoReader
      [useRead,neednframes] = myparse(varargin,'useRead',true,'neednframes',true);

      try
        readerobj = VideoReader(filename);
        headerinfo = get(readerobj);

        headerinfo.type = 'avi';
        headerinfo.nr = headerinfo.Height;
        headerinfo.nc = headerinfo.Width;

        if preload,

          nframes = 0;
          preloaddata = struct;
          preloaddata.ims = {};
          preloaddata.ts = [];

          while true,

            if ~hasFrame(readerobj),
              break;
            end
            preloaddata.ts(end+1) = get(readerobj,'CurrentTime');
            preloaddata.ims{end+1} = readFrame(readerobj);
            nframes = nframes+1;

          end

          assert(numel(preloaddata.ims) == nframes);

        else
          if neednframes && useRead,
            nframes = get(readerobj,'NumberOfFrames');
          end
          if ~neednframes || ~useRead || isempty(nframes),
            % approximate nframes from duration
            nframes = round(get(readerobj,'Duration')*get(readerobj,'FrameRate'));
          end
          headerinfo.readerobj = readerobj;

        end

        headerinfo.nframes = nframes;
        frame_rate = headerinfo.FrameRate ;  % Hz
        dt = 1/frame_rate ;  % s
        headerinfo.timestamps = dt * (0:(nframes-1)) ;  % s

        if preload,

          readframe = @(f) preload_read_frame(preloaddata,f);
          didpreload = true;

        else

          if useRead,
            readframe = @(f) avi_read_frame(readerobj,headerinfo,f);
          else
            readframe = @(f) avi_read_frame_useReadFrame(readerobj,headerinfo,f);
          end

        end
      catch ME_videoreader
        error('Could not open file %s with VideoReader: %s',...
          filename,getReport(ME_videoreader));
      end
    end
  end
end

if preload && ~didpreload,

  preloaddata = struct;
  preloaddata.ims = cell(1,headerinfo.nframes);
  preloaddata.ts = nan(1,headerinfo.nframes);

  for f = 1:headerinfo.nframes,
    [preloaddata.ims{f},preloaddata.ts(f)] = readframe(f);
  end

  readframe = @(f) preload_read_frame(preloaddata,f);
  didpreload = true;

end

function varargout = multi_read_frame(f,readframes)

[im,timestamp] = readframes{1}(f);

for i = 2:numel(readframes),
  im = cat(2,im,readframes{i}(f));
end
varargout{1} = im;
if nargout >= 2,
  varargout{2} = timestamp;
end

function [im,timestamp] = seq_read_frame_piotr(f,sr)
im = [];
timestamp = [];

if ~sr.seek(f-1),
  warning('Could not seek to frame %d',f);
  return;
end

[im,timestamp] = sr.getframe();

function [im,timestamp] = avi_read_frame_useReadFrame(readerobj,headerinfo,f)

timestamp = (f-1)/headerinfo.FrameRate;
oldtimestamp = get(readerobj,'CurrentTime');
if oldtimestamp ~= timestamp,
  set(readerobj,'CurrentTime',timestamp);
end
timestamp = get(readerobj,'CurrentTime');
im = readFrame(readerobj);

function [im,timestamp] = preload_read_frame(preloaddata,f)

try
  im = preloaddata.ims{f};
  timestamp = preloaddata.ts(f);
catch ME,
  if ~isempty(preloaddata.ims),
    im = preloaddata.ims{end};
    timestamp = preloaddata.ts(end);
  end
  warning('Error reading preloaded frame %d: %s',f,getReport(ME));
end

function [im,timestamp] = avi_read_frame(readerobj,headerinfo,f)

try
  im = read(readerobj,f);
catch ME,
  warning('Error reading the first try: %s',getReport(ME));
  pause(.01);
  im = read(readerobj,f);
end
timestamp = (f-1)/headerinfo.FrameRate;

function [im,f] = imseq_read_frame(f,imfiles)

im = imread(imfiles{f});

author: kristinbranson

repo: JAABA


function [ ds ] = TrainValidateMarket(ds, learn_algs, X, X_test, X_query, idxa, idxb, flag)
%
% Input:
%
%  ds - data struct that stores the result
%  learn_algs - algorithms that are used for cross validation
%  X - input matrix, each column is an input vector [DxN*2]. N is the
%  number of pairs.
%  idxa - index of image A in X [1xN]
%  idxb - index of image B in X [1xN]
%
% Output:
%
%  ds - contains the result

c = 1;
clc;
% train
for aC=1:length(learn_algs)
    cHandle = learn_algs{aC};
    fprintf('    training %s ',upper(cHandle.type));
    if aC > 3
        s = learnPairwise(cHandle,X,idxa(1:10:end),idxb(1:10:end),logical(flag(1:10:end)));
    else
        s = learnPairwise(cHandle,X,idxa,idxb,logical(flag));
    end
    if ~isempty(fieldnames(s))
        fprintf('... done in %.4fs\n',s.t);
        ds(c).(cHandle.type) = s;
    else
        fprintf('... not available');
    end
end

% test
testID = importdata('data\testID.mat'); % identity label of all testing bboxes
testCAM = importdata('data\testCAM.mat'); % camera label of all testing bboxes
queryID = importdata('data\queryID.mat'); % identity label of all query bboxes
queryCAM = importdata('data\queryCAM.mat'); % camera label of all query bboxes
names = fieldnames(ds(c));
for nameCounter=1:length(names)
    fprintf('    evaluating %s ',upper(names{nameCounter}));
    [ds(c).(names{nameCounter}).mAP, ds(c).(names{nameCounter}).cmc] = ...
        calcmAP(ds(c).(names{nameCounter}).M, X_test, X_query, testID, testCAM, queryID, queryCAM);
    fprintf('... done \n');
end

end


text:
function [obj, p, g] = runICPInit(obj, p, g)

% Start ------------------------------------------------------------------------

% Start
msg('O', 'SetLogLevel', p.LogLevel);

% Log to file
diary off; fclose all;
p2report = fullfile(obj.OutputFolder, 'ICPLog.txt');
if exist(p2report)==2 % delete old report file if any
    try delete(p2report); end
end
if ~exist(obj.OutputFolder, 'dir')
    mkdir(obj.OutputFolder);
end
diary(p2report);

g.proc = {'GLOBALICP' 'RUNICP'};
msg('S', g.proc, 'LogLevel', 'basic');

% Initialize variables ---------------------------------------------------------

g.nPC = numel(obj.PC);
g.nItICP = 1;
obj.D.stats{g.nItICP}.normdx = Inf;

% Initialize homogeneous transformation matrices
obj.D.H = repmat({eye(4)}, 1, g.nPC);

% Estimate input parameters set to auto ----------------------------------------

if strcmpi(p.HullVoxelSize          , 'auto') || ...
   strcmpi(p.UniformSamplingDistance, 'auto') || ...
   strcmpi(p.MaxDistance            , 'auto') || ...
   strcmpi(p.MaxRoughness           , 'auto')

    msg('S', {g.proc{:} 'ESTIMATE INPUT PARAMETERS'}, 'LogLevel', 'basic');

    if strcmpi(p.HullVoxelSize, 'auto')
        firstTryHullVoxelSize = 5*p.PlaneSearchRadius;
        validVoxelSizes = [0.01 0.02 0.04 0.05 0.08 0.1 0.16 0.2 0.25 0.4 0.5 0.8 1 1.25 2 2.5 4 5 6.25 10 12.5 20 25 50 100]'; % mod(100, voxelSize) has to be zero, see pointCloud.getVoxelSize
        idx = knnsearch(validVoxelSizes, firstTryHullVoxelSize);
        p.HullVoxelSize = validVoxelSizes(idx);
        msg('V', p.HullVoxelSize, 'HullVoxelSize', 'Prec', 2, 'LogLevel', 'basic');
    end

    if strcmpi(p.UniformSamplingDistance, 'auto')
        p.UniformSamplingDistance = 2*p.PlaneSearchRadius;
        msg('V', p.UniformSamplingDistance, 'UniformSamplingDistance', 'Prec', 2, 'LogLevel', 'basic');
    end

    if strcmpi(p.MaxDistance, 'auto')
        p.MaxDistance = 3*p.PlaneSearchRadius;
        msg('V', p.MaxDistance, 'MaxDistance', 'Prec', 2, 'LogLevel', 'basic');
    end

    if strcmpi(p.MaxRoughness, 'auto')
        p.MaxRoughness = 0.2*p.PlaneSearchRadius;
        msg('V', p.MaxRoughness, 'MaxRoughness', 'Prec', 2, 'LogLevel', 'basic');
    end

    msg('E', {g.proc{:} 'ESTIMATE INPUT PARAMETERS'}, 'LogLevel', 'basic');

end

% Assign point cloud dependent parameters --------------------------------------

% Uniform sampling for each point cloud
if numel(p.UniformSamplingDistance) == 1
    p.UniformSamplingDistance = ones(1,g.nPC)*p.UniformSamplingDistance;
end

% Random subsampling for each point cloud
if numel(p.RandomSubsampling) == 1
    p.RandomSubsampling = repmat(p.RandomSubsampling, 1, g.nPC);
end

% Normal subsampling for each point cloud
if numel(p.NormalSubsampling) == 1
    p.NormalSubsampling = repmat(p.NormalSubsampling, 1, g.nPC);
end

% Max leverage subsampling for each point cloud
if numel(p.MaxLeverageSubsampling) == 1
    p.MaxLeverageSubsampling = repmat(p.MaxLeverageSubsampling, 1, g.nPC);
end

% Report input point clouds ----------------------------------------------------

msg('T', '-------------------------------------------------------------------------------', 'LogLevel', 'basic');

msg('T', 'INPUT POINT CLOUDS:', 'LogLevel', 'basic');

msg('T', sprintf('%12s %12s %8s %8s %8s %8s    file', 'point_cloud', 'fix/loose', 'USD', 'RS', 'NS', 'MLS'), 'LogLevel', 'basic');

for i = 1:g.nPC
    if ismember(i, p.IdxFixedPointClouds), fixOrLoose = 'fix  '; else fixOrLoose = 'loose'; end
    if p.RandomSubsampling(i)     , RS  = 'true'; else RS  = 'false'; end
    if p.NormalSubsampling(i)     , NS  = 'true'; else NS  = 'false'; end
    if p.MaxLeverageSubsampling(i), MLS = 'true'; else MLS = 'false'; end

    [~, file] = fileparts(obj.PC{i});

    msg('T', sprintf('%12s %12s %8.2f %8s %8s %8s    %s', ...
        ['[' num2str(i) ']'], ...
        fixOrLoose, ...
        p.UniformSamplingDistance(i), ...
        RS, ...
        NS, ...
        MLS, ...
        file), ...
        'LogLevel', 'basic');

end

msg('T', 'where USD = UniformSamplingDistance', 'LogLevel', 'basic');
msg('T', '      RS  = RandomSubsampling'      , 'LogLevel', 'basic');
msg('T', '      NS  = NormalSubsampling'      , 'LogLevel', 'basic');
msg('T', '      MLS = MaxLeverageSubsampling' , 'LogLevel', 'basic');

msg('T', '-------------------------------------------------------------------------------', 'LogLevel', 'basic');

end

author: pglira


text:
classdef PostProcessColumn < handle

    properties (Access = public)
        m
    end

    properties (Access = private)
        polygon
    end

    properties (Access = private)
        designVariable
        dim
        mesh
        scale
    end

    methods (Access = public)

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

        function plotColumn(obj)
            obj.createPolygon();
            obj.createMesh();
            obj.plotMesh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.designVariable = cParams.designVariable;
            obj.dim            = cParams.dim;
            obj.mesh           = cParams.mesh;
            obj.scale          = cParams.scale;
        end

        function createPolygon(obj)
            z = obj.designVariable.getColumnArea;
            % d     = obj.dim;
            scl   = obj.scale;
            coord = obj.mesh.coord;
            nnod     = obj.mesh.nelem+1;
            dimFig   = 2;
            vertElem = 4;
            vertex = zeros(vertElem*obj.mesh.nelem+1,dimFig);
            for iNod = 1:nnod-1
                vertex(2*iNod-1,:)   = [scl*z(iNod)/2 coord(iNod)];
                vertex(2*iNod,:) = [scl*z(iNod)/2 coord(iNod+1)];
            end
            vertex = obj.flip(vertex,vertElem,dimFig);
            vertex(end,:) = vertex(1,:);
            obj.polygon = polyshape(vertex);
        end

        function createMesh(obj)
            pgon = obj.polygon;
            tr = triangulation(pgon);
            model = createpde;
            tnodes = tr.Points';
            telements = tr.ConnectivityList';
            geometryFromMesh(model,tnodes,telements);
            generateMesh(model,'GeometricOrder','linear','Hmax',0.1)
            coord = model.Mesh.Nodes';
            connec = model.Mesh.Elements';
            s.coord = coord;
            s.connec = connec;
            obj.m = Mesh(s);
        end

        function vertex = flip(obj,vertex,vertElem,dimFig)
            nElem = obj.mesh.nelem;
            nodes = dimFig*nElem+1:vertElem*nElem;
            vertex(nodes,1) = - fliplr(vertex(1:dimFig*nElem,1)')';
            vertex(dimFig*nElem+1:vertElem*nElem,2) = fliplr(vertex(1:dimFig*nElem,2)')';
        end

        function plotMesh(obj)
            figure
            clf
            obj.m.plot();
            grid on
            grid minor
            title('Optimal clamped-clamped column','Interpreter', 'latex','FontSize',20, 'fontweight','b');
            xlabel('A(x)','Interpreter', 'latex','fontsize',14,'fontweight','b');
            ylabel('x','Interpreter', 'latex','fontsize',14,'fontweight','b');
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function pass = test_roots_slow( pref )
% Check that the marching squares and Bezoutian agree with each other.
%%

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e3 * pref.cheb2Prefs.chebfun2eps;
j = 1;

%% (slow one)
f = chebfun2(@(x,y)exp(x-2*x.^2-y.^2).*sin(10*(x+y+x.*y.^2)));
g = chebfun2(@(x,y)exp(-x+2*y.^2+x.*y.^2).*sin(10*(x-y-2*x.*y.^2)));
r1 = roots([f;g],'ms');
r2 = roots([f;g],'resultant');
pass(j) = ( norm(sort(r1(:,1))-sort(r2(:,1))) < tol ); j = j + 1;
pass(j) = ( norm(sort(r1(:,2))-sort(r2(:,2))) < tol ); j = j + 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pass = test_inherited( )
% Test inherited commands, abs, cos, cosh, conj, etc.

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% abs
f = spherefun( @(x,y,z) -1+0*x);
h = spherefun( @(x,y,z) 1+0*x);
g = abs( f );
pass(1) = norm( g - h ) < tol;

% cos
f = spherefun( @(x,y,z) cos(x));
h = spherefun( @(x,y,z) cos(cos(x)));
g = cos( f );
pass(2) = norm( g - h ) < tol;

% cosh
f = spherefun( @(x,y,z) cos(x));
h = spherefun( @(x,y,z) cosh(cos(x)));
g = cosh( f );
pass(3) = norm( g - h ) < tol;

% conj
f = spherefun( @(x,y,z) x);
g = conj( f );
pass(4) = norm( f - g ) < tol;

% ctranspose
f = spherefun( @(x,y,z) y);
g = ctranspose( f );
pass(5) = norm( f - g ) < tol;

% diag
f = spherefun( @(lam,theta) sin(theta).*cos(lam));
h = chebfun( @(x) sin(x).*cos(x), [0 pi]);
g = diag( f );
pass(6) = norm( g - h ) < tol;

% exp
f = spherefun( @(x,y,z) cos(x)+sin(y));
h = spherefun( @(x,y,z) exp(cos(x)+sin(y)));
g = exp( f );
pass(7) = norm( g - h ) < tol;

% flipdim, fliplr, flipud
f = spherefun( @(lam, theta) sin(theta).*cos(lam));
g1 = spherefun( @(lam, theta) sin(theta).*cos(-lam));
g2 = spherefun( @(lam, theta) sin(-theta).*cos(lam));
h1 = fliplr( f );
h2 = flipud( f );
h3 = flipdim( f, 1);
h4 = flipdim( f, 2);
pass(8) = ( norm( g1 - h1 ) < tol );
pass(9) = ( norm( g2 - h2 ) < tol );
pass(10) = ( norm( g2 - h3 ) + norm( g1 - h4 ) < tol );

% imag
f = spherefun( @(x,y,z) x );
g = imag( f );
pass(11) = norm( g  ) < tol;

% isequal
f = spherefun( @(x,y,z) cos(x));
g = spherefun( @(x,y,z) cos(-x) );
pass(12) = isequal( f, g );

f = spherefun( @(x,y,z) cos(x));
g = spherefun( @(x,y,z) cos(x) + 1 );
pass(13) = ~isequal( f, g );

% isreal
f = spherefun( @(x,y,z) x);
pass(14) = isreal(f);

f = spherefun( @(x,y,z) cos(x.*y.*z));
pass(15) = isreal(f);

% iszero
f = spherefun( @(x,y,z) 0*cos(x.*y.*z));
pass(16) = iszero( f );
f = spherefun( zeros(10) );
pass(17) = iszero( f );

% length
f = spherefun( @(x,y,z) cos(z));
[m, n] = length( f );
pass(18) = (m == 1);

% log
f = spherefun( @(x,y,z) exp(z));
h = spherefun( @(x,y,z) z);
g = log( f );
pass(19) = norm( g - h ) < tol;

% rank

% sin
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) sin(cos(x.*y)));
g = sin( f );
pass(20) = norm( g - h ) < tol;

% sinh
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) sinh(cos(x.*y)));
g = sinh( f );
pass(21) = norm( g - h ) < tol;

% size
pass(22) = (size( f,1 ) == inf && size( f,2 ) == inf);

% sqrt
f = spherefun( @(x,y,z) cos(x.*y).^2);
h = spherefun( @(x,y,z) cos(x.*y));
g = sqrt( f );
pass(23) = norm( g - h ) < tol;

% tan
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) tan(cos(x.*y)));
g = tan( f );
pass(24) = norm( g - h ) < 10*tol;

% tand
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) tand(cos(x.*y)));
g = tand( f );
pass(25) = norm( g - h ) < tol;

% tanh
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) tanh(cos(x.*y)));
g = tanh( f );
pass(26) = norm( g - h ) < tol;

% uminus
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) -(cos(x.*y)));
g = uminus( f );
pass(27) = norm( g - h ) < tol;

% uplus
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) +(cos(x.*y)));
g = uplus( f );
pass(28) = norm( g - h ) < tol;

end


text:
function meta=recover_meta(in, dirac)

%meta.order
%meta.scale
%meta.orientation
%meta.dirac_norm
%meta.ave

J=in{1}.meta.j;
L=max(in{2}.meta.theta);

meta.order(1)=1;
meta.scale(1)=-1;
meta.orientation(1)=0;
meta.dirac_norm(1)=norm(dirac{1}.signal{1}(:));
meta.ave(1)=mean(in{1}.signal{1}(:));
r=2;
for m=2:size(in,2)
for l=1:size(in{m}.signal,2)
meta.order(r)=m;
meta.scale(r)=simplestack(in{m}.meta.j(1:end-1,l),J,m-1);
meta.orientation(r)=simplestack(in{m}.meta.theta(:,l)-1,L,m-1);
meta.dirac_norm(r)=norm(dirac{m}.signal{l}(:));
meta.ave(r)=mean(in{m}.signal{l}(:));
r=r+1;
end
end

end

function out=simplestack(code, base, len)
out=0;
for l=1:len
out=base*out+code(l);
end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [E,J]=SynthMeasIsoGPD(d, protocol)
% Computes signals and their derivatives for isotropic free diffusion with
% diffusivity d for the protocol.
%
% [E,J]=SynthMeasIsoGPD(d, protocol)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameter d.
%
% d is the diffusivity.
%
% protocol is the object containing the acquisition protocol.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

if(strcmp(protocol.pulseseq, 'PGSE') || strcmp(protocol.pulseseq, 'STEAM'))

    GAMMA = 2.675987E8;
    modQ = GAMMA*protocol.smalldel'.*protocol.G';
    modQ_Sq = modQ.^2;
    difftime = (protocol.delta'-protocol.smalldel'/3);

    E = exp(-difftime.*modQ_Sq*d);
    if(nargout>1)
        J = -difftime.*modQ_Sq.*E;
    end
elseif(strcmp(protocol.pulseseq, 'FullSTEAM'))

    GAMMA = 2.675987E8;

    FullG = protocol.grad_dirs.*repmat(protocol.G', [1 3]);
    GdGd = sum(FullG.*FullG,2)';
    GcGc = sum(protocol.cG.*protocol.cG,2)';
    GrGr = sum(protocol.rG.*protocol.rG,2)';
    GdGc = sum(FullG.*protocol.cG,2)';
    GdGr = sum(FullG.*protocol.rG,2)';
    GcGr = sum(protocol.cG.*protocol.rG,2)';

    tdd = protocol.gap1 + protocol.gap2 + protocol.TM + 2*protocol.sdelc + 2*protocol.smalldel/3 + 2*protocol.sdelr;
    tcc = protocol.TM + 2*protocol.sdelc/3 + 2*protocol.sdelr;
    trr = protocol.TM + 2*protocol.sdelr/3;
    tdc = protocol.TM + protocol.sdelc + 2*protocol.sdelr;
    tdr = protocol.TM + protocol.sdelr;
    tcr = protocol.TM + protocol.sdelr;

    gqt = GAMMA^2*(GdGd.*protocol.smalldel.^2.*tdd + ...
    GcGc.*protocol.sdelc.^2.*tcc + ...
    GrGr.*protocol.sdelr.^2.*trr + ...
    2*GdGc.*protocol.smalldel.*protocol.sdelc.*tdc + ...
    2*GdGr.*protocol.smalldel.*protocol.sdelr.*tdr + ...
    2*GcGr.*protocol.sdelc.*protocol.sdelr.*tcr);

    E=exp(-gqt'*d);

    if(nargout>1)
        J = -gqt.*E;
    end
elseif(strcmp(protocol.pulseseq, 'DSE'))
    bValue = GetB_ValuesDSE(protocol.G', protocol.delta1', protocol.delta2', protocol.delta3', protocol.t1', protocol.t2', protocol.t3');
    E = exp(-bValue*d);
    if(nargout>1)
        J = -bValue.*E;
    end
elseif(strcmp(protocol.pulseseq, 'OGSE'))
    bValue = GetB_Values(protocol)';
    E = exp(-bValue*d);
    if(nargout>1)
        J = -bValue.*E;
    end
else
    error('Not implemented for pulse sequence: %s', protocol.pulseseq);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function varargout=randiscr(varargin)
%
% For calling details please see v_randiscr.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_randiscr.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: randiscr.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_randiscr(varargin{:});
else
    v_randiscr(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function mri = ft_defacevolume(cfg, mri)

% FT_DEFACEVOLUME allows you to de-identify an anatomical MRI by erasing specific
% regions, such as the face and ears. The graphical user interface allows you to
% position a box over the anatomical data inside which all anatomical voxel values will
% be replaced by zero. You might have to call this function multiple times when both
% face and ears need to be removed. Following defacing, you should check the result
% with FT_SOURCEPLOT.
%
% Use as
%   mri = ft_defacevolume(cfg, mri)
%
% The configuration can contain the following options
%   cfg.method     = 'box', 'spm' (default = 'box')
%
% If you specify the box method, the following options apply
%   cfg.translate  = initial position of the center of the box (default = [0 0 0])
%   cfg.scale      = initial size of the box along each dimension (default is automatic)
%   cfg.rotate     = initial rotation of the box (default = [0 0 0])
%   cfg.selection  = which voxels to keep, can be 'inside' or 'outside' (default = 'outside')
%   cfg.smooth     = 'no' or the FWHM of the gaussian kernel in voxels (default = 'no')
%   cfg.keepbrain  = 'no' or 'yes', segment and retain the brain (default = 'no')
%   cfg.feedback   = 'no' or 'yes', whether to provide graphical feedback (default = 'no')
%
% If you specify no smoothing, the selected area will be zero-masked. If you
% specify a certain amount of smoothing (in voxels FWHM), the selected area will
% be replaced by a smoothed version of the data.
%
% The spm method does not have any options, it uses SPM_DEFACE from the
% SPM12 toolbox.
%
% See also FT_ANONYMIZEDATA, FT_DEFACEMESH, FT_ANALYSISPIPELINE, FT_SOURCEPLOT

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    mri
ft_preamble provenance mri
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% set the defaults
cfg.method    = ft_getopt(cfg, 'method', 'box');
cfg.rotate    = ft_getopt(cfg, 'rotate', [0 0 0]);
cfg.scale     = ft_getopt(cfg, 'scale'); % the automatic default is determined further down
cfg.translate = ft_getopt(cfg, 'translate', [0 0 0]);
cfg.selection = ft_getopt(cfg, 'selection', 'outside');
cfg.smooth    = ft_getopt(cfg, 'smooth', 'no');
cfg.keepbrain = ft_getopt(cfg, 'keepbrain', 'no');
cfg.feedback  = ft_getopt(cfg, 'feedback', 'no');

ismri    = ft_datatype(mri, 'volume') && isfield(mri, 'anatomy');
ismesh   = isfield(mri, 'pos'); % triangles are optional

if ismri
  % check if the input data is valid for this function
  mri = ft_checkdata(mri, 'datatype', 'volume', 'feedback', 'yes');
end

switch cfg.method
  case 'spm'
    % this requires SPM12 on the path
    ft_hastoolbox('spm12', 1);

    % defacing relies on coregistration, which relies on the MRI being reasonably aligned for SPM
    mri = ft_checkdata(mri, 'hascoordsys', 'yes');

    % remember the original transformation matrix and coordinate system
    original = [];
    original.transform = mri.transform;
    original.coordsys  = mri.coordsys;
    mri = ft_convert_coordsys(mri, 'acpc');

    filename1 = {[tempname '.nii']};
    ft_write_mri(filename1{1}, mri, 'dataformat', 'nifti');

    % % apply a least squares pre-alignment step in order to make spm_deface more robust
    % % this could be done conditional on the modality/contrast, which is part of the BIDS filename
    % template = spm_vol(fullfile(spm('Dir'),'canonical','avg152PD.nii'));
    % template = spm_vol(fullfile(spm('Dir'),'canonical','avg152T1.nii'));
    % template = spm_vol(fullfile(spm('Dir'),'canonical','avg152T2.nii'));
    % filevol = spm_vol(filename1{1});
    % M = spm_affreg(template, filevol);
    % spm_get_space(filename1{1}, M * filevol.mat);

    filename2 = spm_deface(filename1);
    mri = ft_read_mri(filename2{1});

    % put the original transformation matrix and coordinate system back
    mri.transform = original.transform;
    mri.coordsys = original.coordsys;

    % clean up the temporary files
    delete(filename1{1});
    delete(filename2{1});

  case 'box'
    % determine the size of the "unit" sphere in the origin and the length of the axes
    switch mri.unit
      case 'mm'
        axmax = 150;
        rbol  = 5;
      case 'cm'
        axmax = 15;
        rbol  = 0.5;
      case 'm'
        axmax = 0.15;
        rbol  = 0.005;
      otherwise
        ft_error('unknown units "%s"', unit);
    end

    figHandle = figure;
    set(figHandle, 'CloseRequestFcn', @cb_close);

    % clear persistent variables to ensure fresh figure
    clear ft_plot_slice

    if ismri
      % the volumetric data needs to be interpolated onto three orthogonal planes
      % determine a resolution that is close to, or identical to the original resolution
      [corner_vox, corner_head] = cornerpoints(mri.dim, mri.transform);
      diagonal_head = norm(range(corner_head));
      diagonal_vox  = norm(range(corner_vox));
      resolution    = diagonal_head/diagonal_vox; % this is in units of "mri.unit"

      % create a contrast enhanced version of the anatomy
      mri.anatomy = double(mri.anatomy);
      dum = unique(mri.anatomy(:));
      clim(1) = dum(round(0.05*numel(dum)));
      clim(2) = dum(round(0.95*numel(dum)));
      anatomy = (mri.anatomy-clim(1))/(clim(2)-clim(1));

      ft_plot_ortho(anatomy, 'transform', mri.transform, 'unit', mri.unit, 'resolution', resolution, 'style', 'intersect');
    elseif ismesh
      if isfield(mri, 'hex')
        ft_plot_mesh(mri, 'surfaceonly', 'yes');
      else
        ft_plot_mesh(mri);
      end
    end

    axis vis3d
    view([110 36]);

    % shift the axes to the left
    ax = get(gca, 'position');
    ax(1) = 0;
    set(gca, 'position', ax);

    % get the xyz-axes
    xdat  = [-axmax 0 0; axmax 0 0];
    ydat  = [0 -axmax 0; 0 axmax 0];
    zdat  = [0 0 -axmax; 0 0 axmax];

    % get the xyz-axes dotted
    xdatdot = (-axmax:(axmax/15):axmax);
    xdatdot = xdatdot(1:floor(numel(xdatdot)/2)*2);
    xdatdot = reshape(xdatdot, [2 numel(xdatdot)/2]);
    n       = size(xdatdot,2);
    ydatdot = [zeros(2,n) xdatdot zeros(2,n)];
    zdatdot = [zeros(2,2*n) xdatdot];
    xdatdot = [xdatdot zeros(2,2*n)];

    % plot axes
    hl = line(xdat, ydat, zdat);
    set(hl(1), 'linewidth', 1, 'color', 'r');
    set(hl(2), 'linewidth', 1, 'color', 'g');
    set(hl(3), 'linewidth', 1, 'color', 'b');
    hld = line(xdatdot, ydatdot, zdatdot);
    for k = 1:n
      set(hld(k    ), 'linewidth', 3, 'color', 'r');
      set(hld(k+n*1), 'linewidth', 3, 'color', 'g');
      set(hld(k+n*2), 'linewidth', 3, 'color', 'b');
    end

    if isempty(cfg.scale)
      cfg.scale = [axmax axmax axmax]/2;
    end

    guidata(figHandle, cfg);

    % add the GUI elements
    cb_creategui(gca);
    cb_redraw(gca);

    uiwait(figHandle);
    cfg = guidata(figHandle);
    delete(figHandle);
    drawnow
    fprintf('keeping all voxels from MRI that are %s the box\n', cfg.selection)

    % the order of application is scale, rotate, translate
    S = cfg.S;
    R = cfg.R;
    T = cfg.T;

    if ismri
      % it is possible to convert the box to headcoordinates, but it is more efficient the other way around
      [X, Y, Z] = ndgrid(1:mri.dim(1), 1:mri.dim(2), 1:mri.dim(3));
      voxpos = ft_warp_apply(mri.transform, [X(:) Y(:) Z(:)]);  % voxel positions in head coordinates
      voxpos = ft_warp_apply(inv(T*R*S), voxpos);               % voxel positions in box coordinates

      remove = ...
        voxpos(:,1) > -0.5 & ...
        voxpos(:,1) < +0.5 & ...
        voxpos(:,2) > -0.5 & ...
        voxpos(:,2) < +0.5 & ...
        voxpos(:,3) > -0.5 & ...
        voxpos(:,3) < +0.5;

    elseif ismesh || issource
      meshpos = ft_warp_apply(inv(T*R*S), mri.pos);               % mesh vertex positions in box coordinates

      remove = ...
        meshpos(:,1) > -0.5 & ...
        meshpos(:,1) < +0.5 & ...
        meshpos(:,2) > -0.5 & ...
        meshpos(:,2) < +0.5 & ...
        meshpos(:,3) > -0.5 & ...
        meshpos(:,3) < +0.5;
    end

    if strcmp(cfg.selection, 'inside')
      % invert the selection, i.e. keep the voxels inside the box
      remove = ~remove;
    end

    if ismri
      if istrue(cfg.keepbrain)
        tmpcfg = [];
        tmpcfg.output = {'brain'};
        seg = ft_volumesegment(tmpcfg, mri);
        fprintf('keeping voxels in brain segmentation\n');
        % keep the tissue of the brain
        remove(seg.brain) = 0;
        clear seg
      end

      if istrue(cfg.feedback)
        tmpmri = keepfields(mri, {'anatomy', 'transform', 'coordsys', 'units', 'dim'});
        tmpmri.remove = remove;
        tmpcfg = [];
        tmpcfg.funparameter = 'remove';
        ft_sourceplot(tmpcfg, tmpmri);
      end

      if isequal(cfg.smooth, 'no')
        fprintf('zero-filling %.0f%% of the volume\n', 100*mean(remove));
        mri.anatomy(remove) = 0;
      else
        tmp = mri.anatomy;
        tmp = (1 + 0.5.*randn(size(tmp))).*tmp; % add 50% noise to each voxel
        tmp = volumesmooth(tmp, cfg.smooth, 'anatomy');
        fprintf('smoothing %.0f%% of the volume\n', 100*mean(remove));
        mri.anatomy(remove) = tmp(remove);
      end

    elseif ismesh
      % determine all fields that might need to be defaced
      fn = setdiff(fieldnames(mri), ignorefields('deface'));
      dimord = cell(size(fn));
      for i=1:numel(fn)
        dimord{i} = getdimord(mri, fn{i});
      end
      % this applies to headshapes and meshes in general
      fprintf('keeping %d and removing %d vertices in the mesh\n', sum(remove==0), sum(remove==1));
      if isfield(mri, 'tri')
        [mri.pos, mri.tri] = remove_vertices(mri.pos, mri.tri, remove);
      elseif isfield(mri, 'tet')
        [mri.pos, mri.tet] = remove_vertices(mri.pos, mri.tet, remove);
      elseif isfield(mri, 'hex')
        [mri.pos, mri.hex] = remove_vertices(mri.pos, mri.hex, remove);
      else
        mri.pos = mri.pos(~remove,1:3);
      end
      for i=1:numel(fn)
        dimtok = tokenize(dimord{i}, '_');
        % do some sanity checks
        if any(strcmp(dimtok, '{pos}'))
          ft_error('not supported');
        end
        if numel(dimtok)>5
          ft_error('too many dimensions');
        end
        % remove the same positions from each matching dimension
        if numel(dimtok)>0 && strcmp(dimtok{1}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(~remove,:,:,:,:);
        end
        if numel(dimtok)>1 && strcmp(dimtok{2}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,~remove,:,:,:);
        end
        if numel(dimtok)>2 && strcmp(dimtok{3}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,:,~remove,:,:);
        end
        if numel(dimtok)>3 && strcmp(dimtok{4}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,:,:,~remove,:);
        end
        if numel(dimtok)>4 && strcmp(dimtok{5}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,:,:,:,~remove);
        end
      end % for fn
      mri = removefields(mri, {'dim', 'transform'}); % these fields don't apply any more
    end % ismesh

    % remove the temporary fields from the configuration, keep the rest for provenance
    cfg = removefields(cfg, {'R', 'S', 'T'});
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous mri
ft_postamble provenance mri
ft_postamble history mri
ft_postamble savevar mri

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_redraw(figHandle, varargin)
persistent p
% define the position of each GUI element
figHandle = get(figHandle, 'parent');
cfg = guidata(figHandle);

rx = str2double(get(findobj(figHandle, 'tag', 'rx'), 'string'));
ry = str2double(get(findobj(figHandle, 'tag', 'ry'), 'string'));
rz = str2double(get(findobj(figHandle, 'tag', 'rz'), 'string'));
tx = str2double(get(findobj(figHandle, 'tag', 'tx'), 'string'));
ty = str2double(get(findobj(figHandle, 'tag', 'ty'), 'string'));
tz = str2double(get(findobj(figHandle, 'tag', 'tz'), 'string'));
sx = str2double(get(findobj(figHandle, 'tag', 'sx'), 'string'));
sy = str2double(get(findobj(figHandle, 'tag', 'sy'), 'string'));
sz = str2double(get(findobj(figHandle, 'tag', 'sz'), 'string'));

% remember the user specified transformation
cfg.rotate    = [rx ry rz];
cfg.translate = [tx ty tz];
cfg.scale     = [sx sy sz];

R = rotate   (cfg.rotate);
T = translate(cfg.translate);
S = scale    (cfg.scale);

% remember the transformation matrices
cfg.R = R;
cfg.T = T;
cfg.S = S;

% start with a cube of unit dimensions
x1 = -0.5;
y1 = -0.5;
z1 = -0.5;
x2 = +0.5;
y2 = +0.5;
z2 = +0.5;

plane1 = [
  x1 y1 z1
  x2 y1 z1
  x2 y2 z1
  x1 y2 z1];

plane2 = [
  x1 y1 z2
  x2 y1 z2
  x2 y2 z2
  x1 y2 z2];

plane3 = [
  x1 y1 z1
  x1 y2 z1
  x1 y2 z2
  x1 y1 z2];

plane4 = [
  x2 y1 z1
  x2 y2 z1
  x2 y2 z2
  x2 y1 z2];

plane5 = [
  x1 y1 z1
  x2 y1 z1
  x2 y1 z2
  x1 y1 z2];

plane6 = [
  x1 y2 z1
  x2 y2 z1
  x2 y2 z2
  x1 y2 z2];

plane1 = ft_warp_apply(T*R*S, plane1);
plane2 = ft_warp_apply(T*R*S, plane2);
plane3 = ft_warp_apply(T*R*S, plane3);
plane4 = ft_warp_apply(T*R*S, plane4);
plane5 = ft_warp_apply(T*R*S, plane5);
plane6 = ft_warp_apply(T*R*S, plane6);

if all(ishandle(p))
  delete(p);
end

p(1) = patch(plane1(:,1), plane1(:,2), plane1(:,3), 'y');
p(2) = patch(plane2(:,1), plane2(:,2), plane2(:,3), 'y');
p(3) = patch(plane3(:,1), plane3(:,2), plane3(:,3), 'y');
p(4) = patch(plane4(:,1), plane4(:,2), plane4(:,3), 'y');
p(5) = patch(plane5(:,1), plane5(:,2), plane5(:,3), 'y');
p(6) = patch(plane6(:,1), plane6(:,2), plane6(:,3), 'y');
set(p, 'FaceAlpha', 0.3);

guidata(figHandle, cfg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_creategui(figHandle, varargin)
% define the position of each GUI element
figHandle = get(figHandle, 'parent');
cfg = guidata(figHandle);

% constants
CONTROL_WIDTH   = 0.05;
CONTROL_HEIGHT  = 0.08;
CONTROL_HOFFSET = 0.68;
CONTROL_VOFFSET = 0.20;

% rotateui
uicontrol('tag', 'rotateui', 'parent', figHandle, 'units', 'normalized', 'style', 'text', 'string', 'rotate', 'callback', [])
uicontrol('tag', 'rx', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.rotate(1)), 'callback', @cb_redraw)
uicontrol('tag', 'ry', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.rotate(2)), 'callback', @cb_redraw)
uicontrol('tag', 'rz', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.rotate(3)), 'callback', @cb_redraw)
ft_uilayout(figHandle, 'tag', 'rotateui', 'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET,                 'vpos', CONTROL_VOFFSET);
ft_uilayout(figHandle, 'tag', 'rx',       'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET);
ft_uilayout(figHandle, 'tag', 'ry',       'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET);
ft_uilayout(figHandle, 'tag', 'rz',       'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET);

% scaleui
uicontrol('tag', 'scaleui', 'parent', figHandle, 'units', 'normalized', 'style', 'text', 'string', 'scale', 'callback', [])
uicontrol('tag', 'sx', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.scale(1)), 'callback', @cb_redraw)
uicontrol('tag', 'sy', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.scale(2)), 'callback', @cb_redraw)
uicontrol('tag', 'sz', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.scale(3)), 'callback', @cb_redraw)
ft_uilayout(figHandle, 'tag', 'scaleui', 'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET,                 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'sx',      'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'sy',      'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'sz',      'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);

% translateui
uicontrol('tag', 'translateui', 'parent', figHandle, 'units', 'normalized', 'style', 'text', 'string', 'translate', 'callback', [])
uicontrol('tag', 'tx', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.translate(1)), 'callback', @cb_redraw)
uicontrol('tag', 'ty', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.translate(2)), 'callback', @cb_redraw)
uicontrol('tag', 'tz', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.translate(3)), 'callback', @cb_redraw)
ft_uilayout(figHandle, 'tag', 'translateui', 'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET,                 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'tx',          'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'ty',          'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'tz',          'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);

% somehow the toolbar gets lost in 2012b
set(figHandle, 'toolbar', 'figure');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_close(figHandle, varargin)
% the figure will be closed in the main function after collecting the guidata
uiresume;


text:
function [metabolite_structure] =createNewMetaboliteStructure(input,source,metabolite_structure_rBioNet,metab,rxnDB)
% This function creates a metabolite structure using the provided input. If
% no VMHId is provided in the header of the input, then VMHId are
% generated.
%
% INPUT
% Input                          Cell array containing the metabolites
%                                   The information provided must be as follows:
%                                   metList={
%                                   'VMH ID' 'metabolite_name' 'HMDB' 'inchistring' 'neutral_formula' 'charged_formula' 'charge'
%                                   'cot' 'Cotinine' 'HMDB0001046'  '' '' '' ''
%                                   'coto' 'Cotinine n-oxide' 'HMDB0001411'  '' '' '' ''
% source                            source of the information contained in metArray
%                                   (e.g., 'Manually assembled by IT')
% metabolite_structure_rBioNet      (optional and only necessary if the
%                                   input has no VMHId)
%
% OUTPUT
% metabolite_structure              metabolite structure
%
% Ines Thiele, 09/2021

%load rbionet data
rBioNetPath =  fileparts(which('tutorial_MetaboAnnotator'));
if exist([rBioNetPath filesep 'cache' filesep 'metab.mat'],'file')
    load([rBioNetPath filesep 'cache' filesep 'metab.mat']);
elseif exist([rBioNetPath filesep 'data' filesep 'metab.mat'],'file')
    load([rBioNetPath filesep 'data' filesep 'metab.mat']);
else
    %TODO
end

if exist([rBioNetPath filesep 'cache' filesep 'rxn.mat'],'file')
    load([rBioNetPath filesep 'cache' filesep 'rxn.mat']);
elseif  exist([rBioNetPath filesep 'data' filesep 'rxn.mat'],'file')
    load([rBioNetPath filesep 'data' filesep 'rxn.mat']);
else
    %TODO
end

RAW = input;

% load the designated field names associated with each metabolites
metaboliteStructureFieldNames;
% find VMH ID in the input file
vmh_col = find(contains(lower(RAW(1,:)),'vmh'));
% clean up some known issues:
name_col = find(contains(lower(RAW(1,:)),'name'));

% test if metabolite name column exists
if isempty(name_col)
    RAW{1} = 'metNames';
    RAW = cat(1, RAW, metList);
else
    RAW(1,name_col) = regexprep(RAW(1,name_col),RAW(1,name_col),'metNames');
end

if isempty(vmh_col) || vmh_col ==0
    if ~exist('metabolite_structure_rBioNet','var')
        load met_strc_rBioNet;
    end
end

% populate metabolite structure
for i = 2 : size(RAW,1)
    name = RAW(i,find(contains(lower(RAW(1,:)),'name')));
    if ~isempty(vmh_col) && vmh_col ~=0 % VMH ID exists
        Ori = RAW{i,vmh_col};

        RAW{i,vmh_col} = regexprep(RAW{i,vmh_col},'-','_minus_');
        RAW{i,vmh_col} = regexprep(RAW{i,vmh_col},'(','_parentO_');
        RAW{i,vmh_col} = regexprep(RAW{i,vmh_col},')','_parentC_');
        VMHId = RAW{i,vmh_col};
    else
        [VMHId] = generateVMHMetAbbr(RAW{i,name_col},metabolite_structure_rBioNet,metab,rxnDB);
        Ori = VMHId;
        vmh_col = 0;
    end
    % adjust the abbreviations to the requirements of Matlab field
    % names
    % the original designated VMH ids listed in the field VMHId

    %make a temporary VMHId name in case the existing one is an invalid field name
    VMHId_suggestion = VMHId;
    if  isvarname(strcat('VMHId_', VMHId))
        fieldname = strcat('VMH_',VMHId);
    else

        VMHId_suggestion = VMHId;
        %VMHId = datestr(now,'yyyymmddTHHMMSS');

        fieldname = strcat('VMH_',VMHId);
    end

    metabolite_structure.(fieldname) = struct();
    metabolite_structure.(fieldname).VMHId_suggestion = VMHId_suggestion;

    for j = 1:size(RAW,2)
        if j~=vmh_col
            % check whether a columnn header is part of the field2Add
            % definition
            colName1 =  regexprep(RAW(1,j),' ','_');
            colName2 =  regexprep(RAW(1,j),' ','');
            if ~isempty(find(contains(lower(field2Add),lower(colName1)))) || ~isempty(find(contains(lower(field2Add),lower(colName2))))
                clear F
                % assign correct field name
                if ~isempty(find(contains(lower(field2Add),lower(colName1))))
                    F=  field2Add(find(contains(lower(field2Add),lower(colName1))));
                elseif ~isempty(find(contains(lower(field2Add),lower(colName2))))
                    F=  field2Add(find(contains(lower(field2Add),lower(colName2))));
                end
                % make sure that the novel entry is not NaN or empty
                if isempty(find(isnan(RAW{i,j}))) && ~isempty(RAW{i,j}) && ~isempty(F)
                    metabolite_structure.(strcat('VMH_',VMHId)).(F{1}) = RAW{i,j};
                    metabolite_structure.(strcat('VMH_',VMHId)).([F{1},'_source']) = [source ': ' datestr(now)];
                else
                    metabolite_structure.(strcat('VMH_',VMHId)).(F{1}) = NaN;
                    metabolite_structure.(strcat('VMH_',VMHId)).([F{1},'_source']) = NaN;
                end
                metabolite_structure.(strcat('VMH_',VMHId)).VMHId = Ori;
            end
        end
    end
end
% add the remaining fields to the metabolite structure
metabolite_structure= addField2MetStructure(metabolite_structure);
% clean up known, potential issues with the input data
[metabolite_structure] = cleanUpMetabolite_structure(metabolite_structure);
% check for known, potentially remaining issues with the metabolite IDs
% any errors will be removed.
removeErrors = 1;
[metabolite_structure,errorFlag] = sanityCheckMetIds(metabolite_structure,removeErrors);

author: opencobra

repo: cobratoolbox


text:
%The field of definition of the quadratic function
%have been known,
%solve the range of values.
%rov(a,b,c,x1,x2),(x1<x2)
function rov(a,b,c,x1,x2)
a=sym(a);b=sym(b);c=sym(c);x1=sym(x1);x2=sym(x2);
dd=simple((-b)/(2*a));jj=simple((4*a*c-b^2)/(4*a));
fx1=simple(a*(x1^2)+b*x1+c);fx2=simple(a*(x2^2)+b*x2+c);
if double(a)>0
   if double(x2)<=double(dd)
      fmax=fx1
      fmin=fx2
   elseif double(x1)>=double(dd)
      fmax=fx2
      fmin=fx1
   elseif double(x1)<double(dd)&double(x2)>double(dd)
      fmax=big(fx1,fx2)
      fmin=jj
   end
else
   if double(x2)<=double(dd)
      fmax=fx2
      fmin=fx1
   elseif double(x1)>=double(dd)
      fmax=fx1
      fmin=fx2
   elseif double(x1)<double(dd)&double(x2)>double(dd)
      fmax=jj
      fmin=small(fx1,fx2)
   end
end


text:
function r8vec_midspace_test ( )

%*****************************************************************************80
%
%% R8VEC_MIDSPACE_TEST tests R8VEC_MIDSPACE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 June 2012
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_MIDSPACE_TEST\n' );
  fprintf ( 1, '  For a R8VEC:\n' );
  fprintf ( 1, '  R8VEC_MIDSPACE: evenly spaced midpoints between A and B\n' );

  a = 10.0;
  b = 20.0;

  x = r8vec_midspace ( n, a, b );
  r8vec_print ( n, x, '  r8vec_midspace ( 5, 10, 20 )' );

  return
end


text:
function output = DesignBlitEqualizer(raw_BLIT_response)
% DesignBlitEqualizer -- Equalize raw BLIT (Band Limitted Impulse Train)
% to have flat gain using FIR equalizer.
% This equalizer assumes Nuttall-based BLIT with zero at 0.55 times
% sampling frequency.
% This is a supporting function for  fractional pitch allocation of
% excitation pulse
%
% output = DesignBlitEqualizer(raw_BLIT_response)
%
% Input argument
%   raw_BLIT_response: impulse response of antialiasing BLIT filter
%
% Output
%   output : FIR equalizer to flatten frequency response

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

fftl = 8192; % Finally the central 121 points are used.
% The truncation length 121 is decided for the response to decrease
% relatively small enough.
amp = abs(fft(raw_BLIT_response / sum(raw_BLIT_response), fftl));
eqAmp = 1.0 ./ (amp + 10.0 ^ (-45 / 20));
output = real(fftshift(ifft(eqAmp)));
output = output(fftl / 2 + 1 + (-60:60));
end

author: google


text:
function blas1_d_test12 ( )

%*****************************************************************************80
%
%% TEST12 tests IDAMAX, DAXPY, DDOT, DSCAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST12\n' );
  fprintf ( 1, '  Use IDAMAX, DAXPY, DDOT and DSCAL\n' );
  fprintf ( 1, '  in a Gauss elimination routine.\n' );
  fprintf ( 1, '\n' );
%
%  Set the matrix.
%
  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = 2.0;
      elseif ( i == j + 1 )
        a(i,j) = - 1.0;
      elseif ( i == j - 1 )
        a(i,j) = - 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end
%
%  Set the right hand side.
%
  b(1:n-1) = 0.0;
  b(n) = n + 1.0;

  info = 0;

  for k = 1 : n - 1

    l = idamax ( n-k+1, a(k:n,k), 1 ) + k - 1;
    ipvt(k) = l;

    if ( a(l,k) == 0.0 )

      info = k;

    else

      if ( l ~= k )
        t = a(l,k);
        a(l,k) = a(k,k);
        a(k,k) = t;
      end

      t = -1.0 / a(k,k);
      a(k+1:n,k) = dscal ( n-k, t, a(k+1:n,k), 1 );

      for j = k+1 : n
        t = a(l,j);
        if ( l ~= k )
          a(l,j) = a(k,j);
          a(k,j) = t;
        end
        a(k+1:n,j) = daxpy ( n-k, t, a(k+1:n,k), 1, a(k+1:n,j), 1 );
      end

    end

  end

  ipvt(n) = n;
  if ( a(n,n) == 0.0 )
    info = n;
  end

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The matrix is singular.\n' );
    return
  end

  for k = 1 : n-1
    l = ipvt(k);
    t = b(l);
    if ( l ~= k )
      b(l) = b(k);
      b(k) = t;
    end
    b(k+1:n) = daxpy ( n-k, t, a(k+1:n,k)', 1, b(k+1:n), 1 );
  end

  for k = n : -1 : 1
    b(k) = b(k) / a(k,k);
    t = - b(k);
    b(1:k-1) = daxpy ( k-1, t, a(1:k-1,k)', 1, b(1:k-1), 1 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution (should be (1,2,3,...,n)):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12f\n', b(i) );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Trdata,Vldata] = getTrVlData(TrsRGB_dir, TrXYZ_dir, ...
    VlsRGB_dir, VlXYZ_dir, aug, Tr_ImageNum, Vl_ImageNum, imgSize)

%get training imgs
Tr_in_images = imageDatastore(TrsRGB_dir, 'ReadFcn', @inRead);

%get training imgs GT
Tr_gt_images = imageDatastore(TrsRGB_dir);

Tr_gt_images.ReadFcn = @(filename)gtRead(filename, ...
    TrsRGB_dir, TrXYZ_dir);

%get random images instead of using the entire dataset
if Tr_ImageNum ~=0  && Tr_ImageNum < length(Tr_in_images.Files)
    inds = randperm(Tr_ImageNum);
    Tr_in_images.Files = Tr_in_images.Files(inds(1:Tr_ImageNum));
    Tr_gt_images.Files = Tr_gt_images.Files(inds(1:Tr_ImageNum));
end

%get validation imgs
Vl_in_images = imageDatastore(VlsRGB_dir, 'ReadFcn', @inRead);

%get validation imgs GT
Vl_gt_images = imageDatastore(VlsRGB_dir);

Vl_gt_images.ReadFcn = @(filename)gtRead(filename, ...
    VlsRGB_dir, VlXYZ_dir);

%get random images instead of using the entire dataset
if Vl_ImageNum ~=0  && Vl_ImageNum < length(Vl_in_images.Files)
    inds = randperm(Vl_ImageNum);
    Vl_in_images.Files = Vl_in_images.Files(inds(1:Vl_ImageNum));
    Vl_gt_images.Files = Vl_gt_images.Files(inds(1:Vl_ImageNum));
end

if aug == 1
    augmenter = imageDataAugmenter( ...
        'RandScale',[1 1.2], ... %scaling
        'RandXReflection',true, ... %X reflection
        'RandYReflection',true); %Y reflection

    Trdata = randomPatchExtractionDatastore(Tr_in_images,Tr_gt_images, ...
        imgSize,'DataAugmentation',augmenter,'PatchesPerImage',...
        1);
elseif aug == 0
    Trdata = randomPatchExtractionDatastore(Tr_in_images,Tr_gt_images, ...
        imgSize,'PatchesPerImage', 1);
end

Vldata = randomPatchExtractionDatastore(Vl_in_images,Vl_gt_images, ...
        imgSize,'PatchesPerImage', 2);

author: mahmoudnafifi


text:
function x = VBA_spm_invBcdf(F,v,w,tol)
% Inverse Cumulative Distribution Function (CDF) of Beta distribution
% FORMAT x = spm_invBcdf(F,v,w,tol)
%
% F   - CDF (lower tail p-value)
% v   - Shape parameter (v>0)
% w   - Shape parameter (w>0)
% x   - Beta ordinates at which CDF F(x)=F
% tol - tolerance [default 10^-6]
%__________________________________________________________________________
%
% spm_invBcdf implements the inverse Cumulative Distribution Function
% for Beta distributions.
%
% Returns the Beta-variate x such that Pr{X<x} = F for X a Beta
% random variable with shape parameters v and w.
%
% Definition:
%--------------------------------------------------------------------------
% The Beta distribution has two shape parameters, v and w, and is
% defined for v>0 & w>0 and for x in [0,1] (See Evans et al., Ch5).
% The Cumulative Distribution Function (CDF) F(x) is the probability
% that a realisation of a Beta random variable X has value less than
% x. F(x)=Pr{X<x}: This function is usually known as the incomplete Beta
% function. See Abramowitz & Stegun, 26.5; Press et al., Sec6.4 for
% definitions of the incomplete beta function.
%
% Variate relationships:
%--------------------------------------------------------------------------
% Many: See Evans et al., Ch5
%
% Algorithm:
%--------------------------------------------------------------------------
% Interval bisection to find the inverse CDF to within tol.
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%__________________________________________________________________________
% Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_invBcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Parameters
%--------------------------------------------------------------------------
Dtol  = 10^-8;
maxIt = 10000;

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<4, tol = Dtol; end
if nargin<3, error('Insufficient arguments'), end

ad = [ndims(F);ndims(v);ndims(w)];
rd = max(ad);
as = [[size(F),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))];...
      [size(w),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size');
end

%-Computation - Undefined and special cases
%--------------------------------------------------------------------------
%-Initialise result to zeros
x = zeros(rs);

%-Only defined for F in [0,1] & strictly positive v & w.
% Return NaN if undefined.
md = ( F>=0  &  F<=1  &  v>0  &  w>0 );
if any(~md(:))
    x(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Special cases: x=0 when F=0, x=1 when F=1
x(md & F==1) = 1;

%-Compute where defined & not special case
%--------------------------------------------------------------------------
Q  = find( md  &  F>0  &  F<1 );
if isempty(Q), return, end
if xa(1), FQ=F(Q); FQ=FQ(:); else FQ=F*ones(length(Q),1); end
if xa(2), vQ=v(Q); vQ=vQ(:); else vQ=v*ones(length(Q),1); end
if xa(3), wQ=w(Q); wQ=wQ(:); else wQ=w*ones(length(Q),1); end

%-Interval bisection
%--------------------------------------------------------------------------
a  = zeros(length(Q),1); fa = a-FQ;
b  =  ones(length(Q),1); fb = b-FQ;
i  = 0;
xQ = a+1/2;
QQ = 1:length(Q);

while ~isempty(QQ) &&  i<maxIt
    fxQQ        = betainc(xQ(QQ),vQ(QQ),wQ(QQ))-FQ(QQ);
    mQQ         = fa(QQ).*fxQQ > 0;
    a(QQ(mQQ))  = xQ(QQ(mQQ));   fa(QQ(mQQ))  = fxQQ(mQQ);
    b(QQ(~mQQ)) = xQ(QQ(~mQQ));  fb(QQ(~mQQ)) = fxQQ(~mQQ);
    xQ(QQ)      = a(QQ) + (b(QQ)-a(QQ))/2;
    QQ          = QQ( (b(QQ)-a(QQ))>tol );
    i           = i+1;
end

if i==maxIt
    warning('convergence criteria not reached - maxIt reached');
end

x(Q) = xQ;

author: MBB-team

repo: VBA-toolbox


text:
% DEMSWISSROLLLLE4 Demonstrate LLE on the oil data.

% MLTOOLS

[Y, lbls] = lvmLoadData('swissRoll');

options = lleOptions(32);
model = lleCreate(2, size(Y, 2), Y, options);
model = lleOptimise(model, 2);

lvmScatterPlotColor(model, model.Y(:, 2));

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, model.Y(:, 2), 'SwissRoll', 4, true);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function plot_comparison(P1,P2,fig,ms1,ms2)

% Plots two point clouds "P1" and "P2" so that those points of "P2" which are
% not in "P1" are plotted in red whereas the common points are plotted in
% blue. "fig" and "ms1" and "ms2" are the figure number and marker sizes.

if nargin == 3
    ms1 = 3;
    ms2 = 3;
elseif nargin == 4
    ms2 = 3;
end

if ms1 == 0
    ms1 = 3;
end
if ms2 == 0
    ms2 = 3;
end

P2 = setdiff(P2,P1,'rows');

figure(fig)
if size(P1,2) == 3
    plot3(P1(:,1),P1(:,2),P1(:,3),'.b','Markersize',ms1)
    hold on
    plot3(P2(:,1),P2(:,2),P2(:,3),'.r','Markersize',ms2)
elseif size(P1,2) == 2
    plot(P1(:,1),P1(:,2),'.b','Markersize',ms1)
    hold on
    plot(P2(:,1),P2(:,2),'.r','Markersize',ms2)
end
hold off
axis equal

author: InverseTampere

repo: TreeQSM


text:
function b = diag(a,k)
% function B = diag(A,K)
%
% DESCRIPTION
%   Diagonal polynomial matrices and diagonals of polynomial matrices.
%
% INPUTS
%   A: polynomial
%   K: integer
%
% OUTPUTS
%   B: polynomial
%
% SYNTAX
%   B=diag(A,K)
%     If A is a polynomial vector of length N then B is an M-by-M
%     square matrix with M=N+ABS(K) and the elements of A on the K^th
%     diagonal.  K=0 is the main diagonal, K>0 is above the main diagonal
%     and K<0 is below the main diagonal
%   B=diag(A,K)
%     If A is a polynomial matrix then B is a column vector whose elements
%     are the entries of the K^th diagonal of A.
%   B=diag(A)
%     This is the same as B=diag(A,0).
%
% See also triu, tril, blkdiag

% 12/3/2009: PJS  Initial Coding

if nargin==1
    k=0;
end

sza = size(a);
if any(sza==1)
    % a is a vector --> b is a matrix

    % Use double/diag to find right dimensions and indices
    la = length(a);
    bmat = diag(1:la,k);
    idx = find(bmat);

    % Assign values
    b = polynomial( zeros(size(bmat)) );
    L.type = '()';
    L.subs = {idx};
    b = subsasgn(b,L,a);
else
    % a is a matrix --> b is a vector
    if k>=sza(2) || k <= -sza(1)
        b = polynomial(zeros(0,1));
        return
    end

    % Use double/diag to find right indices
    amat = reshape( 1:(sza(1)*sza(2)), sza );
    idx = diag(amat,k);

    % Assign values
    L.type = '()';
    L.subs = {idx};
    b = subsref(a,L);
end

author: yu-jiang

repo: radpbook

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


text:
function g = lfmKernGradient(kern, t, varargin)

% LFMKERNGRADIENT Gradient of LFM kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the latent force
% model kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in t.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t1 : the input locations associated with the rows of the
% kernel matrix.
% ARG t2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as t1 and the same number of columns
% as t2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO lfmKernParamInit, kernGradient, lfmKernDiagGradient, kernGradX,
% lfmXlfmKernGradient
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

if length(varargin)<2
  t2 = t;
else
  t2 = varargin{1};
end

[g1, g2] = lfmXlfmKernGradient(kern, kern, t, t2, varargin{end});

g = real(g1 + g2);

%g = real(g1);


text:
classdef TestAgastFeatureDetector
    %TestAgastFeatureDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
    end

    methods (Static)
        function test_detect_img
            obj = cv.AgastFeatureDetector('Threshold',10, 'NonmaxSuppression',true);
            assert(obj.Threshold == 10);
            assert(obj.NonmaxSuppression == true);
            typename = obj.typeid();

            img = imread(TestAgastFeatureDetector.im);
            kpts = obj.detect(img);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestAgastFeatureDetector.kfields, fieldnames(kpts))));
        end

        function test_detect_imgset
            img = imread(TestAgastFeatureDetector.im);
            imgs = {img, img};
            obj = cv.AgastFeatureDetector();

            kpts = obj.detect(imgs);
            validateattributes(kpts, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(kpt) validateattributes(kpt, {'struct'}, {'vector'}), kpts);
            cellfun(@(kpt) assert(all(ismember(TestAgastFeatureDetector.kfields, fieldnames(kpt)))), kpts);
        end

        function test_detect_mask
            img = imread(TestAgastFeatureDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image
            obj = cv.AgastFeatureDetector();
            kpts = obj.detect(img, 'Mask',mask);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestAgastFeatureDetector.kfields, fieldnames(kpts))));
            xy = cat(1, kpts.pt);
            assert(all(xy(:,1) <= ceil(size(img,2)/2)));
        end

        function test_detect_types
            types = {'AGAST_5_8', 'AGAST_7_12d', 'AGAST_7_12s', 'OAST_9_16'};
            img = imread(TestAgastFeatureDetector.im);
            obj = cv.AgastFeatureDetector();
            for i=1:numel(types)
                obj.Type = types{i};
                kpts = obj.detect(img);
                validateattributes(kpts, {'struct'}, {'vector'});
            end
        end

        function test_error_1
            try
                cv.AgastFeatureDetector('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv