text:
function SE = functionPowerOptimization_maxmin(signal,interference,Pmax,prelogFactor)
%Compute DL power allocation that solves the max-min fairness problem,
%using the algorithm in Theorem 7.1.
%
%This function require additional software packages to be used, which
%need to be downloaded and installed separately. These packages are
%developed independently and are delivered with separate licenses.
%The implementation uses CVX (http://cvxr.com/cvx) and has been tested
%using CVX version 2.1. We recommend the use of the Mosek solver (we
%have tested using version 7.1.0.12).
%
%INPUT:
%signal       = K x L matrix where element (k,j) is a_jk in (7.2)
%interference = K x L x K x L matrix where (l,i,j,k) is b_lijk in (7.3)
%Pmax         = Maximum transmit power per BS
%prelogFactor = Prelog factor
%
%OUTPUT:
%SE = K x L matrix where element (k,j) is the downlink SE of UE k in cell j
%     using the max-min power allocation solution
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.01 (Last edited: 2019-03-16)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Extract number of UEs
K = size(signal,1);

%Extract number of cells
L = size(signal,2);

%Check which UEs that have non-zero channels, because these ones are
%excluded (they are considered inactive)
nonzero = reshape(signal,[K*L 1]);
nonzero = nonzero(nonzero>0);

%Initalize the gamma-variables in Algorithm 1
rateLower = 0;
rateUpper = log2(1+Pmax*min(nonzero));

%Set the accuracy of the bisection
delta = 0.01;

%Prepare to save the power solution
rhoBest = zeros(K,L);

%Solve the max-min problem by bisection - iterate until the difference
%between the lower and upper points in the interval is smaller than delta
while norm(rateUpper - rateLower) > delta

    %Compute the midpoint of the line. Note that we are performing the
    %bisection in the SE domain instead of the SINR domain as in Algorithm
    %1, since we can then specify delta as the SE difference
    rateCandidate = (rateLower+rateUpper)/2;

    %Transform the midpoints into SINR requirements
    gammaCandidate = 2.^(rateCandidate)-1;

    %Solve the feasibility problem using CVX
    [feasible,rhoSolution] = functionFeasibilityProblem_cvx(signal,interference,Pmax,K,L,gammaCandidate);

    %If the problem was feasible, then replace rateLower with
    %gammaCandidate and store rhoSolution as the new best solution.
    if feasible
        rateLower = rateCandidate;
        rhoBest = rhoSolution;
    else
        %If the problem was not feasible, then replace ratePoint with
        %gammaCandidate
        rateUpper = rateCandidate;
    end

end

%Compute the SEs using Theorem 4.6
SE = functionComputeSE_DL_poweralloc(rhoBest,signal,interference,prelogFactor);

function [feasible,rhoSolution] = functionFeasibilityProblem_cvx(signal,interference,Pmax,K,L,SINR)
%Solve the linear feasibility problem in Algorithm 1 using CVX, by adding
%an extra variable so that it becomes a minimization problem with better
%properties.

cvx_begin
cvx_quiet(true); % This suppresses screen output from the solver

variable rho(K,L);  %Variable for the K x L power allocation matrix
variable scaling    %Scaling parameter for power constraints

minimize scaling %Minimize the power indirectly by scaling the power constraints

subject to

for j = 1:L

    for k = 1:K

        if signal(k,j)>0

            %SINR constraints
            SINR*(sum(sum(rho.*interference(:,:,k,j))) + 1) - (rho(k,j)*signal(k,j)) <= 0

        end

        rho(k,j)>=0

    end

    sum(rho(:,j)) <= scaling*Pmax;

end

scaling >= 0; %Power constraints must be positive

cvx_end

%% Analyze the CVX output and prepare the output variables
if isempty(strfind(cvx_status,'Solved')) %Both the power minimization problem and the feasibility problem are infeasible
    feasible = false;
    rhoSolution = [];
elseif scaling>1 %Only the power minimization problem is feasible
    feasible = false;
    rhoSolution = rho;
else %Both the power minimization problem and feasibility problem are feasible
    feasible = true;
    rhoSolution = rho;
end

author: emilbjornson

repo: massivemimobook


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR:  Landmark Based Registration, affine transformation
%
% - load data (see setup2DhandData)
% - setup  viewer (viewImage2D), imgModelpolator (splineimgModel),
% - setup landmarks (LM)
% - run affine
%==============================================================================

clear, close all, help(mfilename)

%% setup hand data
setup2DhandData

if FAIRinput(mfilename,'set new landmarks ? ',0),
  [LM,fig] = getLandmarks(dataT,dataR,omega,m);
  close(fig);
end;

omegaT = omega(1,:);
omegaR = omega(end,:);
xT = getCellCenteredGrid(omegaT,m);
xR = getCellCenteredGrid(omegaR,m);
Tc = imgModel(dataT,omegaT,xT);
Rc = imgModel(dataR,omegaR,xR);

%% visualize data
FAIRfigure(1,'figname',mfilename); clf;
subplot(1,3,1); viewImage(Tc,omegaT,m); hold on;
ph = plotLM(LM(:,1:2),'numbering','on','color','r');
set(ph,'linewidth',2,'markersize',20);
title(sprintf('%s','T&LM'),'fontsize',20);

subplot(1,3,2); viewImage(Rc,omegaR,m); hold on;
ph = plotLM(LM(:,3:4),'numbering','on','color','g','marker','+');
set(ph,'linewidth',2,'markersize',20);
title(sprintf('%s','R&LM'),'fontsize',20);

%% compute landmark based registration
[yc,LM] = LMreg('linear',LM(:,1:4),xR);
TLM = imgModel(dataT,omegaT,yc);

subplot(1,3,3); cla; viewImage(TLM,omegaR,m); hold on;
ph = plotLM(LM(:,3:4),'numbering','off','color','g','marker','+');
qh = plotLM(LM(:,7:8),'numbering','off','color','m','marker','x');
rh = plot(LM(:,[3,7])',LM(:,[4,8])','m-','linewidth',3);
set([ph;qh;rh],'linewidth',2,'markersize',20);
title(sprintf('%s','T(y^{affine})&LM'),'fontsize',20);
%==============================================================================


text:
function C3 = C3f(epsi, C3x)
%C3F  Evaluate C_3
%
%   C3 = C3F(EPSI, C3X) evaluates C_{3,l} using Eq. (25) and the
%   coefficient vector C3X.  EPSI is a K x 1 array.  C3X is a 1 x 15 array.
%   C3 is a K x 5 array.

  nC3 = 6;
  nC3x = size(C3x, 2);
  j = nC3x;
  C3 = zeros(length(epsi), nC3 - 1);
  for k = nC3 - 1 : -1 : 1
    t = C3(:, k);
    for i = nC3 - k : -1 : 1
      t = epsi .* t + C3x(j);
      j = j - 1;
    end
    C3(:, k) = t;
  end
  mult = ones(length(epsi), 1);
  for k = 1 : nC3 - 1
    mult = mult .* epsi;
    C3(:, k) = C3(:, k) .* mult;
  end
end

author: Sable

repo: mcbench-benchmarks


text:
function [L,a,b] = rgb2lab(R,G,B)
% function [L, a, b] = RGB2Lab(R, G, B)
% RGB2Lab takes matrices corresponding to Red, Green, and Blue, and
% transforms them into CIELab.  This transform is based on ITU-R
% Recommendation  BT.709 using the D65 white point reference.
% The error in transforming RGB -> Lab -> RGB is approximately
% 10^-5.  RGB values can be either between 0 and 1 or between 0 and 255.
% By Mark Ruzon from C code by Yossi Rubner, 23 September 1997.
% Updated for MATLAB 5 28 January 1998.

if (nargin == 1)
  B = double(R(:,:,3));
  G = double(R(:,:,2));
  R = double(R(:,:,1));
end

if ((max(max(R)) > 1.0) | (max(max(G)) > 1.0) | (max(max(B)) > 1.0))
  R = R/255;
  G = G/255;
  B = B/255;
end

[M, N] = size(R);
s = M*N;

% Set a threshold
T = 0.008856;

RGB = [reshape(R,1,s); reshape(G,1,s); reshape(B,1,s)];

% RGB to XYZ
MAT = [0.412453 0.357580 0.180423;
       0.212671 0.715160 0.072169;
       0.019334 0.119193 0.950227];
XYZ = MAT * RGB;

X = XYZ(1,:) / 0.950456;
Y = XYZ(2,:);
Z = XYZ(3,:) / 1.088754;

XT = X > T;
YT = Y > T;
ZT = Z > T;

fX = XT .* X.^(1/3) + (~XT) .* (7.787 .* X + 16/116);

% Compute L
Y3 = Y.^(1/3);
fY = YT .* Y3 + (~YT) .* (7.787 .* Y + 16/116);
L  = YT .* (116 * Y3 - 16.0) + (~YT) .* (903.3 * Y);

fZ = ZT .* Z.^(1/3) + (~ZT) .* (7.787 .* Z + 16/116);

% Compute a and b
a = 500 * (fX - fY);
b = 200 * (fY - fZ);

L = reshape(L, M, N);
a = reshape(a, M, N);
b = reshape(b, M, N);

if ((nargout == 1) | (nargout == 0))
  L = cat(3,L,a,b);
end

author: Cloud-CV

repo: object-proposals


text:
function f = make(varargin)
%MAKE   Constructor shortcut for SINGFUN objects.
%   For SINGFUN methods implemented at the ONEFUN level, it is not possible to
%   call the class constructor file corresponding to a ONEFUN object directly.
%   F = MAKE(VARARGIN) allows us to get around this and construct a SINGFUN.
%
% See also SINGFUN.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Developer note: This is essentially a "factory method" in the sense of
% _Design Patterns_ by Gamma, Helm, Johnson, and Vlissides.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f = singfun(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

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

% Last Modified by GUIDE v2.5 13-Apr-2012 20:57:39

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

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

% --- Executes just before gui_reviewresults is made visible.
function gui_reviewresults_OpeningFcn(hObject, eventdata, handles, varargin)
% get data and options
data = varargin{1};
% set up default options for utl_summarize_structarray
opts = hlp_varargin2struct(varargin(2:end),'leftindent',[0 15],'rightindent',[0 15], ...
    'rewrite_fields',{'kld','KL divergence','nll','Negative log-likelihood', ...
                      'mcr','Error rate', 'mae','Mean abs error', 'mse','Mean square error',...
                      'max','Max absolut error','rms','Root mean square error','bias','Mean bias',...
                      'medse','Median square error','auc','Neg. area under ROC', ...
                      'cond_entropy','Conditional entropy','cross_entropy','Cross-Entropy', ...
                      'f_measure','Negative F-score','classes','Target classes','class_ratio','Class ratio', ...
                      'per_class','Per-class error','TP','True positive rate','TN','True negative rate',...
                      'FP','False positive rate','FN','False negative rate','measure','Measure',...
                      'targ','Targets','pred','Predictions','time','Runtime','target','Targets',...
                      'prediction','Predictions','expression','Reproducing expression','value','Value'});

% pick out the correct part of the array
if isscalar(data)
    sub_arrays = {};
    % check if some of the sub-fields are a struct arrays...
    for f=fieldnames(data)'
        field = data.(f{1});
        if isstruct(field) && ~isscalar(field)
            sub_arrays{end+1} = field; end
    end
    if length(sub_arrays) == 1
        % exactly one sub-array: take this as the data to summarize...
        data = sub_arrays{1};
    end
end

if isfield(data,'timestamp')
    data = rmfield(data,'timestamp'); end

% set the summary text
[text,handles.data] = summarize_structarray(data,opts);
set(handles.edit1,'String',[{''} text]);
set(handles.edit1,'Position',[0 0 1 1],'Units','normalized');

% create uitable and popuate with data
conts = handles.data.conts;
fullconts = handles.data.fullconts;
fnames = handles.data.fnames;
% no column is editable
columneditable = false(1,length(fnames));
% by default all columns are of type 'char'
columnformat = repmat({'char'},1,length(fnames));
for c=1:length(fullconts)
    col = fullconts{c};
    % check if all fields of a column are numeric or logical
    if all(cellfun(@isnumeric,conts{c}))
        columnformat{c} = 'numeric';
        for r=1:length(col)
            if isempty(col{r}) || ~isnumeric(col{r})
                col{r} = NaN; end
        end
    elseif all(cellfun(@islogical,conts{c}))
        columnformat{c} = 'logical';
    else
        % otherwise we convert to char what is not yet char
        for r=find(~cellfun('isclass',col,'char'))
            col{r} = hlp_tostring(col{r}); end
    end
    % ... and put it into the table data
    tabledata(:,c) = col;
end
handles.tabledata = tabledata;
handles.table1 = uitable('Data',tabledata,'ColumnName',fnames,'Parent',handles.uipanel2, ...
    'Units','normalized','Position',[0 0 1 1],'ColumnEditable',columneditable,'ColumnFormat',columnformat);
set(handles.table1,'Units','pixels')
wh = get(handles.table1,'Position');
colwidth = (wh(3)-35) / length(fnames);
set(handles.table1,'Units','normalized')
set(handles.table1,'ColumnWidth',repmat({round(colwidth)},1,length(fnames)));

% Choose default command line output for gui_reviewresults
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes gui_reviewresults wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = gui_reviewresults_OutputFcn(hObject, eventdata, handles)
try
    varargout{1} = handles.output;
    close(handles.figure1);
catch
end

% help button
function pushbutton5_Callback(hObject, eventdata, handles)
disp('coming soon...');

% ok button
function pushbutton1_Callback(hObject, eventdata, handles)
uiresume(handles.figure1);

function figure1_KeyPressFcn(hObject, eventdata, handles)
if strcmp(eventdata.Key,'return') || strcmp(eventdata.Key,'escape')
    pushbutton1_Callback(handles.pushbutton1, eventdata, handles); end

% Explore button
function pushbutton4_Callback(hObject, eventdata, handles)
function pushbutton2_Callback(hObject, eventdata, handles)
function pushbutton3_Callback(hObject, eventdata, handles)
function figure1_CreateFcn(hObject, eventdata, handles)
function uipanel1_CreateFcn(hObject, eventdata, handles)
function edit1_Callback(hObject, eventdata, handles)
function edit1_CreateFcn(hObject, eventdata, handles)
function pushbutton6_Callback(hObject, eventdata, handles)
io_mkdirs('home:/.bcilab/results/');
[FileName,PathName] = uiputfile('*.csv','Save table',env_translatepath('home:/.bcilab/results/untitled.csv'));
if FileName
    csvwrite([PathName,FileName],handles.tabledata);  end

function pushbutton6_ButtonDownFcn(hObject, eventdata, handles)

function [text,results] = summarize_structarray(data,varargin)
% Summarize a struct array's fields in text form.
% Text = summarize_structarray(Data,Options...)
%
% This function takes a struct array and calculates summary statistics for each of its numeric
% fields, and outputs a textual report of the results. This is for the review results GUI.
%
% In:
%   Data : A struct array for which a summary report shall be generated.
%
%   Options... : additional name-value pairs; possible names are:
%
%                 --- content control ---
%                 'retain_items' : cell array of checker functions (as in cellfun) that all retained
%                                  data items must satisfy (default: {@isvector,@isnumeric,@isfinite})
%
%                 'retain_fields' : cell array of checker functions (as in cellfun) that all
%                                   retained fields must satisfy (default: {@isvector})
%
%                 'rewrite_fields' : cell array of rewrite rules for field names
%                                   (e.g.{'mcr','Mis-classification rate','auc','Area under Curve'})
%
%                 --- formatting control ---
%
%                 'lhs' : format of the left-hand side of each line (default: '%s')
%
%                 'sep' : separator string between left-hand side and right-hand side (default: ' : ')
%
%                 'rhs' : expression for the right-hand side of each line (in terms of X, the data)
%                         (default: 'sprintf(''%.3f +/- $.3f (N=%d)'',mean(X),std(X),length(X))')
%
%                 'leftspacing' : number of characters to reserve on the left-hand-side;
%                                minimum marging to first character, and minimum spacing to separator
%                                (default: [5 30])
%
%                 'rightspacing': number of characters to reserve on the right-hand-side;
%                                minimum marging from last character, and minimum spacing from
%                                separator (default: [5 15])
%
% Out:
%   Text : a textual summary, akin to what disp(data(1)) produces, but with means and std. devs
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-08

% get the options
hlp_varargin2struct(varargin, 'retain_items', {@(x)isscalar(x),@(x)~isempty(x),@isnumeric,@isfinite}, 'retain_fields',{@(x)~isempty(x)}, ...
    'rewrite_fields',{}, 'lhs','%s', 'sep',' : ', 'rhs','sprintf(''%.3f +/- %.3f (N=%d)'',mean(X),std(X),length(X))', ...
    'leftspacing',[5 0], 'rightspacing', [5 0]);

% get the field names
fnames = fieldnames(data)';
% get the raw contents
raw = struct2cell(data(:));
% reformat the contents
for r=1:size(raw,1)
   conts{r} = {raw{r,:}}; end

% reduce the data items
for r=1:length(conts)
    retain{r} = 1:length(conts{r}); end
for checker = retain_items
    for r=1:length(conts)
        retain{r} = retain{r}(cellfun(checker{1},conts{r}(retain{r}))); end
end
% reduce conts
for r=1:length(conts)
    fullconts{r} = cell(1,length(conts{r}));
    conts{r} = conts{r}(retain{r});
    fullconts{r}(retain{r}) = conts{r};
end

% reduce the fields
mask = 1:length(fnames);
for checker = retain_fields
    mask = mask(cellfun(checker{1},conts(mask))); end
fnames = fnames(mask);
fullconts = fullconts(mask);
conts = conts(mask);

% rewrite field names
for f=1:length(fnames)
    fnames{f} = hlp_rewrite(fnames{f},rewrite_fields{:}); end %#ok<USENS>

% create the text
text = cell(1,length(conts));
for r=1:length(conts)
    X = [conts{r}{:}]; %#ok<NASGU>
    left{r} = sprintf(lhs,fnames{r});
    right{r} = eval(rhs);
end

leftspacing = max(leftspacing(2),leftspacing(1)+max(cellfun('length',left))); %#ok<NODEF>
rightspacing = max(rightspacing(2),rightspacing(1)+max(cellfun('length',right))); %#ok<NODEF>
for r=1:length(conts)
    left{r} = [repmat(' ',1,leftspacing - length(left{r})) left{r}];
    right{r} = [right{r} repmat(' ',1,rightspacing - length(right{r})) ];
    text{r} = [left{r} sep right{r}];
end

if nargout == 0
    fprintf('\n');
    % print the text to the console
    for r=1:length(text)
        fprintf([text{r} '\n']); end
    fprintf('\n');
else
    % create results structure
    results.text = text;
    results.fnames = fnames;
    results.conts = conts;
    results.fullconts = fullconts;
end

repo: REST


text:
function outputFilename = publish_mathjax(file, opts, varargin)
    %PUBLISH_MATHJAX  Publish to HTML, and use MathJax to render equations
    %
    % ## Input
    % * __file__ M-file to publish.
    % * __opts__ Structure of options. By default this function uses a custom
    %   XSL stylesheet to publish as HTML with MathJax rendered equations, i.e
    %   `struct('format','html', 'stylesheet','mxdom2mathjax.xsl')`
    %
    % ## Output
    % * __outputFilename__ path to the generated HTML document.
    %
    % ## Options
    % Accepts the same option as the publish function.
    %
    % ## Example
    %
    %     html = publish_mathjax('my_script.m');
    %     web(html, '-browser')
    %
    %     showdemo('my_script')
    %
    % ## References
    % Inspired by a MathWorks Support Team
    % [solution](https://www.mathworks.com/matlabcentral/answers/93851).
    %
    % See also: publish, grabcode
    %

    % path to this directory
    dname = fileparts(mfilename('fullpath'));

    % HTML publish options with custom XSL
    if nargin < 2, opts = struct(); end
    opts.format = 'html';
    opts.stylesheet = fullfile(dname, 'mxdom2mathjax.xsl');

    % publish
    outputFilename = publish(file, opts, varargin{:});

    % copy CSS file to output directory
    outputCSS = fullfile(fileparts(outputFilename), 'publish_custom.css');
    if true || exist(outputCSS, 'file') ~= 2
        copyfile(fullfile(dname, 'publish_custom.css'), outputCSS);
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function M = median_filter(im)
  % MEDIAN_FILTER simpler wrapper for calling medfilt2 on each channel
  %
  % M = median_filter(im)
  %
  % Input:
  %   im  w by h by c image
  % Output:
  %   M  median filtered image in each channel
  %

  % This could probably be a one-liner using num2cell and cellfun
  M = zeros(size(im));
  for c = 1 : size(im,3)
    M(:,:,c) = medfilt2(im(:,:,c));
  end
end

author: alecjacobson

repo: gptoolbox


text:
function submit(partId)
%SUBMIT Submit your code and output to the ml-class servers
%   SUBMIT() will connect to the ml-class server and submit your solution

  fprintf('==\n== [ml-class] Submitting Solutions | Programming Exercise %s\n==\n', ...
          homework_id());
  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid homework part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  [login password] = loginPrompt();
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to ml-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    % Get Challenge
    [login, ch, signature] = getChallenge(login);
    if isempty(login) || isempty(ch) || isempty(signature)
      % Some error occured, error string in first return element.
      fprintf('\n!! Error: %s\n\n', login);
      return
    end

    % Attempt Submission with Challenge
    ch_resp = challengeResponse(login, password, ch);
    [result, str] = submitSolution(login, ch_resp, partId, output(partId), ...
                                 source(partId), signature);

    fprintf('\n== [ml-class] Submitted Homework %s - Part %d - %s\n', ...
            homework_id(), partId, partNames{partId});
    fprintf('== %s\n', strtrim(str));
    if exist('OCTAVE_VERSION')
      fflush(stdout);
    end
  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '8';
end

function [partNames] = validParts()
  partNames = { 'Estimate Gaussian Parameters', ...
                'Select Threshold' ...
                'Collaborative Filtering Cost', ...
                'Collaborative Filtering Gradient', ...
                'Regularized Cost', ...
                'Regularized Gradient' ...
                };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'estimateGaussian.m' }, ...
           { 'selectThreshold.m' }, ...
           { 'cofiCostFunc.m' }, ...
           { 'cofiCostFunc.m' }, ...
           { 'cofiCostFunc.m' }, ...
           { 'cofiCostFunc.m' }, ...
           };
end

function out = output(partId)
  % Random Test Cases
  n_u = 3; n_m = 4; n = 5;
  X = reshape(sin(1:n_m*n), n_m, n);
  Theta = reshape(cos(1:n_u*n), n_u, n);
  Y = reshape(sin(1:2:2*n_m*n_u), n_m, n_u);
  R = Y > 0.5;
  pval = [abs(Y(:)) ; 0.001; 1];
  yval = [R(:) ; 1; 0];
  params = [X(:); Theta(:)];
  if partId == 1
    [mu sigma2] = estimateGaussian(X);
    out = sprintf('%0.5f ', [mu(:); sigma2(:)]);
  elseif partId == 2
    [bestEpsilon bestF1] = selectThreshold(yval, pval);
    out = sprintf('%0.5f ', [bestEpsilon(:); bestF1(:)]);
  elseif partId == 3
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 0);
    out = sprintf('%0.5f ', J(:));
  elseif partId == 4
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 0);
    out = sprintf('%0.5f ', grad(:));
  elseif partId == 5
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 1.5);
    out = sprintf('%0.5f ', J(:));
  elseif partId == 6
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 1.5);
    out = sprintf('%0.5f ', grad(:));
  end
end

function url = challenge_url()
  url = 'http://www.ml-class.org/course/homework/challenge';
end

function url = submit_url()
  url = 'http://www.ml-class.org/course/homework/submit';
end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          while ~feof(fid)
            line = fgets(fid);
            src = [src line];
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature] = getChallenge(email)
  str = urlread(challenge_url(), 'post', {'email_address', email});

  str = strtrim(str);
  [email, str] = strtok (str, '|');
  [ch, str] = strtok (str, '|');
  [signature, str] = strtok (str, '|');
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'homework', homework_id(), ...
            'part', num2str(part), ...
            'email', email, ...
            'output', output, ...
            'source', source, ...
            'challenge_response', ch_resp, ...
            'signature', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Password: ', 's');
end

function [str] = challengeResponse(email, passwd, challenge)
  salt = ')~/|]QMB3[!W`?OVt7qC"@+}';
  str = sha1([challenge sha1([salt email passwd])]);
  sel = randperm(numel(str));
  sel = sort(sel(1:16));
  str = str(sel);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function sig = signature(b,x0,y0)
%SIGNATURE Computes the signature of a boundary.
%   SIG = SIGNATURE(B,X0,Y0) computes the signature of boundary. B. A
%   signature is defined as the distance (DIST) from (X0,Y0) to the
%   boundary, as a function of angle. B is an np-by-2 matrix, with np >=
%   360, whose rows contain the (x,y) = (row,col) coordinates of the
%   boundary ordered in a clockwise or counterclockwise direction. If
%   (X0,Y0) is not included in the input argument, the centroid of the
%   boundary is used by default. SIG is a vector of size size 360-by-1,
%   indicating a signature resolution of one degree. The input must be a
%   one-pixel-thick boundary obtained, for example, by using function
%   bwboundaries.
%
%   If (X0,Y0) or the default centroid is outside the boundary, the
%   signature is not defined and an error is issued.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Check dimensions of b.
[np,nc] = size(b);
if (np < 360 || nc ~= 2)
   error('b must be of size (>=360)-by-2.');
end

% Some boundary tracing programs, such as bwboundaries.m, result in a
% sequence in which the coordinates of the first and last points are the
% same. If this is the case, in b, eliminate the last point.
if isequal(b(1,:),b(np,:))
   b = b(1:np - 1,:);
   np = np - 1;
end

% If (x0,y0) is not specified, set the origin from which the signature
% is computed as the centroid.
if nargin == 1
   % Coordinates of the centroid.
   x0 = sum(b(:,1))/np;
   y0 = sum(b(:,2))/np;
end

% Check to see that (x0,y0) is inside the boundary.
IN = inpolygon(x0,y0,b(:,1),b(:,2));
if ~IN
   error('(x0,y0) or centroid is not inside the boundary.')
end

% Shift origin of coordinate system to (x0,y0).
b(:,1) = b(:,1) - x0;
b(:,2) = b(:,2) - y0;

% Convert the coordinates to polar. But first have to convert the given
% image coordinates, (x,y), to the coordinate system used by MATLAB for
% conversion between Cartesian and polar cordinates. Designate these
% coordinates by (xcart,ycart). The two coordinate systems are related
% as follows:  xcart = y and ycart = -x, where (x,y) = (row,col).
xcart = b(:,2);
ycart = -b(:,1);
[theta,rho] = cart2pol(xcart,ycart);

% Convert angles to degrees.
theta = theta.*(180/pi);

% Convert to all nonnegative angles.
j = theta == min(theta(:));
k = theta == max(theta(:));
theta = theta.*(0.5*abs(1 + sign(theta)))...
   - 0.5*(-1 + sign(theta)).*(360 + theta);

% Set the smallest angle to 0 and the largest angle to 359.
theta(j) = 0;
theta(k) = 359;

% Sort angles in ascending order and pair them with the corresponding
% rho
[theta,idx] = sort(theta);
rho(idx) = rho;
n = size(theta,1);

% Subsample the signature into 360 angle increments.
n = size(rho,1);
idx = round(linspace(1,n,360));
rho = rho(idx);
sig = rho;

author: dipum

repo: dipum-toolbox

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


text:
function handles = plotData3D(hObject, handles)

plotDimensions = get(handles.lstPlotDimensions, 'Value');

FigureTitle = 'Data';
[handles, handles.figData] = openPlotFigure(hObject, handles, ...
    'Data Plot (3D)', FigureTitle);

view(3);

scatter3(handles.Data(plotDimensions(1), :), ...
         handles.Data(plotDimensions(2), :), ...
         handles.Data(plotDimensions(3), :), ...
         getPlotMarkerSize(), ...
         ['k' getPlotMarkerStyle()]);

hold off;

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [w,run] = train_is(x,y,w,lambda)
% Iterative scaling

% Written by Thomas P Minka

if any(x(:) < 0)
  error('Iterative scaling must have x > 0')
  %x = x - min(x(:));
end

if nargin < 4
  lambda = 0;
end
[d,n] = size(x);
flops(0);
step = 1/max(sum(x,1));
i1 = (y > 0)';
% sum of all data
xs = row_sum(x);
% sum of data in class 1
x1 = x*i1;
% sum of data in class 2
x2 = x*(1-i1);
% ratio
r12 = x1./x2;
flops(flops + n + flops_row_sum(x));
if nargout > 1
  run.w = [];
  run.flops = [];
  run.e = [];
end
for iter = 1:10000
  old_w = w;
  s1 = 1./(1+exp(-(w'*x)))';
  xs1 = x*s1;
  %delta1 = x1./xs1;
  %w = w + step*log(delta1);
  %delta2 = (xs-xs1)./x2;
  %w = w + step*log(delta2);
  delta = xs./xs1 - 1;
  r12 = (x1 - lambda*w)./(x2 - lambda*w);
  w = w + step*log(r12.*delta);
  if iter == 1
    fl = flops_mul(w',x)+n*(flops_exp+2) + ...
	flops_mul(x,s1)+2*d + ...
	d*(1+flops_exp+2);
  end
  flops(flops + fl);

  if nargout > 1 & rem(iter,100) == 1
    run.w(:,end+1) = w;
    run.flops(end+1) = flops;
    run.e(end+1) = logProb(scale_cols(x,y),w) -0.5*lambda*w'*w;
  end
  if rem(iter,1000) == 0
    fprintf('IS iter %d\n', iter)
  end

  if max(abs(w - old_w)) < 1e-5
    break
  end
end
figure(2)
plot(run.e)
if iter == 10000
  warning('not enough iters')
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%kPlot3D '3D Plot Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Plot3D.pane file
%
% Parameters:
% MultiChoice: plot3DXOrientation 'X Orientation', default: 1: 'x orientation'
%    Choices are:
%   1: 'Width'
%   2: 'Height'
%   3: 'Depth'
%   4: 'Time'
%   5: 'Elements'
% MultiChoice: plot3DYOrientation 'Y Orientation', default: 2: 'y orientation'
%    Choices are:
%   1: 'Width'
%   2: 'Height'
%   3: 'Depth'
%   4: 'Time'
%   5: 'Elements'
% Integer: plot3DWidthOffset 'Width Offset', default: 0: 'width offset'
% Integer: plot3DHeightOffset 'Height Offset', default: 0: 'height offset'
% Integer: plot3DDepthOffset 'Depth Offset', default: 0: 'depth offset'
% Integer: plot3DTimeOffset 'Time Offset', default: 0: 'time offset'
% Integer: plot3DElementsOffset 'Elements Offset', default: 0: 'elements offset'
%
% Example: kPlot3D( {'plot3DXOrientation',1;'plot3DYOrientation',2;'plot3DWidthOffset',0;'plot3DHeightOffset',0;'plot3DDepthOffset',0;'plot3DTimeOffset',0;'plot3DElementsOffset',0})
%
% Khoros helpfile follows below:
%
% The 3D plot object supports the display of a 3D plot.
% .begin tagged
% .item "Plot Type"
% This list selection lets you set the plot type that is used by the
% 3D plot object.  Choices include: "Line Plot," "Wireframe," "Mesh,"
% "Horizon," "Scatter," "Impulse," "Contour 3D," "Contour 2D,"  or
% "Constant Shade.
% .item "Line Type"
% There are seven line types available: "Solid," "Dotted," "Dot Dash,"
% "Short Dash," "Long Dash," "Odd Dash," and "Grid Dotted."  The line type
% is not used with scatter plots or shaded plots.
% .item "Marker Type"
% This attribute only applies to plots that use a marker, ie, scatter plots.
% There is a wide variety of marker types that may be used with scatter plots.
% .item "Surface Shade Type"
% This attribute only applies to the constant shade plot type.
% It specifies what part of the data is to dictate the shading.
% "Imagery" specifies that shading is to be done on imagery (or color)
% ie, using the "d" value of each data point as defined by its C
% oord structure.  "Elevation" specifies that shading is to be done on
% elevation, ie, using the "z" value of each data point as defined
% by its Coord structure.  When "Normals" is used, the normal to each
% data point is computed, and shading is done using the normals.
% .item "Color Origination"
% When this attribute is set to "Use Foreground Color," the plot is colored in
% the foreground color.   When set to "Use Data Values," the plot is colored
% according to the plot data (see manual for more details).
% .item "Foreground Color"
% This stringlist selection lets you set the foreground color of the 3D
% plot, using the color name.  It is only used when "Color Origination"
% is set to "Use Foreground Color."
% .item "Fill Color"
% This stringlist selection lets you set the color that is used with
% filled plots.  The only filled plots are the mesh and the horizon plots.
% .item "X and Y Orientation"
% These toggles only take effect when the plot data input is provided via a
% data object (as opposed to an array of Coords).  Furthermore, they only apply
% when the plot data is stored in the value segment of the data object (as
% opposed to the location segment of the data object).  They dictate how the
% 3D plot data is to be extracted from the value segment of the data object.
% The 3D plot data is always extracted as a surface, but that surface may be
% defined over any two of the dimensions supported by the value segment. For
% example, the default surface orientation is (width x height), but you might
% specify it to be (height x depth),  (width x time), or any other combination
% of the 5 supported dimensions.
% .item "Width, Height, Depth, Time, And Elements Offsets"
% These integer selections may be set in conjunction with the "X and Y
% Orientation" toggles.  When neither the X or Y orientation" set to "Width,"
% the "Width Offset" may be used to specify the width offset at which the
% surface of plot data is to be extracted.  Similarly, "Height Offset,"
% "Depth Offset," "Time Offset" and "Elements Offset" may be used to specify
% offsets into height, depth, time, and elements whenever that dimension is
% "not" being plotted.  Note that an offset may only be used when that
% dimension of the value data is greater than 1.
% .item "Rotation About X"
% This parameter allows you to rotate the camera (the eye) about the X axis.
% Values are given in degrees, from 0 to 360.
% .item "Rotation About Y"
% This parameter allows you to rotate the camera (the eye) about the Y axis.
% Values are given in degrees, from 0 to 360.
% .item "Rotation About Z"
% This parameter allows you to rotate the camera (the eye) about the Z axis.
% Values are given in degrees, from 0 to 360.
% .item "Eye Distance"
% This parameter specifies the distance between the plot and the eye.
% The default value of 0 is right up against the bounding box of the
% plot.
% .end tagged

function varargout = kPlot3D(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kPlot3D(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'plot3DXOrientation', 1;'plot3DYOrientation', 2;'plot3DWidthOffset', 0;'plot3DHeightOffset', 0;'plot3DDepthOffset', 0;'plot3DTimeOffset', 0;'plot3DElementsOffset', 0};
maxval={0,0,1,1,1,1,1};
minval={0,0,1,1,1,1,1};
istoggle=[0,0,0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'MultiChoice','MultiChoice','Integer','Integer','Integer','Integer','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'Plot3D" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function pass = test_laplacian(pref)
% Test LAPLACIAN

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 50*pref.cheb3Prefs.chebfun3eps;

% Check definition:
F = chebfun3(@(x,y,z) cos(x.*z)+sin(y.*z));
lapF = laplacian(F);
lapF1 = diff(F, 2, 1) + diff(F, 2, 2) + diff(F, 2, 3);
pass(1) = norm(lapF1 - lapF) < tol;

% Check definition:
F = chebfun3(@(x,y,z) cos(x) + y.*z + z.^2);
lapF = laplacian(F);
lapF1 = diff(F, 2, 1) + diff(F, 2, 2) + diff(F, 2, 3);
pass(2) = norm(lapF1 - lapF) < tol;

end


text:
% Computing most probable explanation.

% If you don't break ties consistently, loopy can give wrong mpe
% even though the graph has no cycles, and even though the max-marginals are the same.
% This example was contributed by Wentau Yih <wtyih@yahoo.com> 29 Jan 02.

% define loop-free graph structure (all edges point down)
%
% Xe1   Xe2
%  |    |
%  E1   E2
%    \ /
%     R1
%     |
%    Xr1

N = 6;
dag = zeros(N,N);
Xe1 = 1; Xe2 = 2; E1 = 3; E2 = 4; R1 = 5; Xr1 = 6;
dag(Xe1, E1) = 1;
dag(Xe2, E2) = 1;
dag([E1 E2], R1) = 1;
dag(R1, Xr1) = 1;

node_sizes = [ 1 1 2 2 2 1 ];

% create BN

bnet = mk_bnet(dag, node_sizes, 'observed', [Xe1 Xe2 Xr1]);

% fill in CPT

bnet.CPD{Xe1} = tabular_CPD(bnet, Xe1, [1]);
bnet.CPD{Xe2} = tabular_CPD(bnet, Xe2, [1]);
bnet.CPD{E1} = tabular_CPD(bnet, E1, [0.2 0.8]);
bnet.CPD{E2} = tabular_CPD(bnet, E2, [0.3 0.7]);
bnet.CPD{R1} = tabular_CPD(bnet, R1, [1 1 1 0.8 0 0 0 0.2]);
bnet.CPD{Xr1} = tabular_CPD(bnet, Xr1, [0.15 0.85]);

clear engine;
engine{1} = belprop_inf_engine(bnet);
engine{2} = jtree_inf_engine(bnet);
engine{3} = global_joint_inf_engine(bnet);
engine{4} = var_elim_inf_engine(bnet);

evidence = cell(1,N);
evidence{Xe1} = 1;  evidence{Xe2} = 1;  evidence{Xr1} = 1;

mpe = find_mpe(engine{1}, evidence, 'break_ties', 0) % gives wrong results
mpe = find_mpe(engine{1}, evidence)
for i=2:4
  mpe = find_mpe(engine{i}, evidence)
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function newdata=jdatadecode(data,varargin)
%
% newdata=jdatadecode(data,opt,...)
%
% Convert all JData object (in the form of a struct array) into an array
% (accepts JData objects loaded from either loadjson/loadubjson or
% jsondecode for MATLAB R2016b or later)
%
% This function implements the JData Specification Draft 3 (Jun. 2020)
% see http://github.com/NeuroJSON/jdata for details
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      data: a struct array. If data contains JData keywords in the first
%            level children, these fields are parsed and regrouped into a
%            data object (arrays, trees, graphs etc) based on JData
%            specification. The JData keywords are
%               "_ArrayType_", "_ArraySize_", "_ArrayData_"
%               "_ArrayIsSparse_", "_ArrayIsComplex_",
%               "_ArrayZipType_", "_ArrayZipSize", "_ArrayZipData_"
%      opt: (optional) a list of 'Param',value pairs for additional options
%           The supported options include
%               Recursive: [1|0] if set to 1, will apply the conversion to
%                            every child; 0 to disable
%               Base64: [0|1] if set to 1, _ArrayZipData_ is assumed to
%                         be encoded with base64 format and need to be
%                         decoded first. This is needed for JSON but not
%                         UBJSON data
%               Prefix: ['x0x5F'|'x'] for JData files loaded via loadjson/loadubjson, the
%                         default JData keyword prefix is 'x0x5F'; if the
%                         json file is loaded using matlab2018's
%                         jsondecode(), the prefix is 'x'; this function
%                         attempts to automatically determine the prefix;
%                         for octave, the default value is an empty string ''.
%               FullArrayShape: [0|1] if set to 1, converting _ArrayShape_
%                         objects to full matrices, otherwise, stay sparse
%               MaxLinkLevel: [0|int] When expanding _DataLink_ pointers,
%                         this sets the maximum level of recursion
%               FormatVersion: [2|float]: set the JSONLab output version;
%                         since v2.0, JSONLab uses JData specification Draft 1
%                         for output format, it is incompatible with all
%                         previous releases; if old output is desired,
%                         please set FormatVersion to 1
%
% output:
%      newdata: the covnerted data if the input data does contain a JData
%               structure; otherwise, the same as the input.
%
% examples:
%      obj={[],{'test'},true,struct('sparse',sparse(2,3),'magic',uint8(magic(5)))}
%      jdata=jdatadecode(jdataencode(obj))
%      isequaln(obj,jdata)
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

    newdata=data;
    opt=struct;
    if(nargin==2)
        opt=varargin{1};
    elseif(nargin>2)
        opt=varargin2struct(varargin{:});
    end
    opt.fullarrayshape=jsonopt('FullArrayShape',0,opt);
    opt.maxlinklevel=jsonopt('MaxLinkLevel',0,opt);

    %% process non-structure inputs
    if(~isstruct(data))
        if(iscell(data))
            newdata=cellfun(@(x) jdatadecode(x,opt),data,'UniformOutput',false);
        elseif(isa(data,'containers.Map'))
            newdata=containers.Map('KeyType',data.KeyType,'ValueType','any');
            names=data.keys;
            for i=1:length(names)
                newdata(names{i})=jdatadecode(data(names{i}),opt);
            end
        end
        return;
    end

    %% assume the input is a struct below
    fn=fieldnames(data);
    len=length(data);
    needbase64=jsonopt('Base64',0,opt);
    format=jsonopt('FormatVersion',2,opt);
    if(isoctavemesh)
        prefix=jsonopt('Prefix','',opt);
    else
        prefix=jsonopt('Prefix','x0x5F',opt);
    end
    if(~isfield(data,N_('_ArrayType_')) && isfield(data,'x_ArrayType_'))
        prefix='x';
        opt.prefix='x';
    end

    %% recursively process subfields
    if(jsonopt('Recursive',1,opt)==1)
      for i=1:length(fn) % depth-first
        for j=1:len
            if(isstruct(data(j).(fn{i})) || isa(data(j).(fn{i}),'containers.Map'))
                newdata(j).(fn{i})=jdatadecode(data(j).(fn{i}),opt);
            elseif(iscell(data(j).(fn{i})))
                newdata(j).(fn{i})=cellfun(@(x) jdatadecode(x,opt),newdata(j).(fn{i}),'UniformOutput',false);
            end
        end
      end
    end

    %% handle array data
    if(isfield(data,N_('_ArrayType_')) && (isfield(data,N_('_ArrayData_')) || isfield(data,N_('_ArrayZipData_'))))
      newdata=cell(len,1);
      for j=1:len
        if(isfield(data,N_('_ArrayZipSize_')) && isfield(data,N_('_ArrayZipData_')))
            zipmethod='zip';
            if(isstruct(data(j).(N_('_ArrayZipSize_'))))
                data(j).(N_('_ArrayZipSize_'))=jdatadecode(data(j).(N_('_ArrayZipSize_')),opt);
            end
	    dims=data(j).(N_('_ArrayZipSize_'))(:)';
            if(length(dims)==1)
                 dims=[1 dims];
            end
            if(isfield(data,N_('_ArrayZipType_')))
                zipmethod=data(j).(N_('_ArrayZipType_'));
            end
            if(ismember(zipmethod,{'zlib','gzip','lzma','lzip','lz4','lz4hc','base64'}))
                decompfun=str2func([zipmethod 'decode']);
                arraytype=data(j).(N_('_ArrayType_'));
                chartype=0;
                if(strcmp(arraytype,'char') || strcmp(arraytype,'logical'))
                    chartype=1;
                    arraytype='uint8';
                end
                if(needbase64 && strcmp(zipmethod,'base64')==0)
                    ndata=reshape(typecast(decompfun(base64decode(data(j).(N_('_ArrayZipData_')))),arraytype),dims);
                else
                    ndata=reshape(typecast(decompfun(data(j).(N_('_ArrayZipData_'))),arraytype),dims);
                end
                if(chartype)
                    ndata=char(ndata);
                end
            else
                error('compression method is not supported');
            end
        else
            if(isstruct(data(j).(N_('_ArrayData_'))))
                data(j).(N_('_ArrayData_'))=jdatadecode(data(j).(N_('_ArrayData_')),opt);
            end
            if(isstruct(data(j).(N_('_ArrayData_'))) && isfield(data(j).(N_('_ArrayData_')),N_('_ArrayType_')))
                data(j).(N_('_ArrayData_'))=jdatadecode(data(j).(N_('_ArrayData_')),varargin{:});
            end
            if(iscell(data(j).(N_('_ArrayData_'))))
                data(j).(N_('_ArrayData_'))=cell2mat(cellfun(@(x) double(x(:)),data(j).(N_('_ArrayData_')),'uniformoutput',0)).';
            end
            ndata=cast(data(j).(N_('_ArrayData_')),char(data(j).(N_('_ArrayType_'))));
        end
        if(isfield(data,N_('_ArrayZipSize_')))
            if(isstruct(data(j).(N_('_ArrayZipSize_'))))
                data(j).(N_('_ArrayZipSize_'))=jdatadecode(data(j).(N_('_ArrayZipSize_')),opt);
            end
	    dims=data(j).(N_('_ArrayZipSize_'))(:)';
            if(iscell(dims))
                dims=cell2mat(dims);
            end
            if(length(dims)==1)
                 dims=[1 dims];
            end
            ndata=reshape(ndata(:),fliplr(dims));
            ndata=permute(ndata,ndims(ndata):-1:1);
        end
        iscpx=0;
        if(isfield(data,N_('_ArrayIsComplex_')) && isstruct(data(j).(N_('_ArrayIsComplex_'))) )
                data(j).(N_('_ArrayIsComplex_'))=jdatadecode(data(j).(N_('_ArrayIsComplex_')),opt);
        end
        if(isfield(data,N_('_ArrayIsComplex_')) && data(j).(N_('_ArrayIsComplex_')) )
                iscpx=1;
        end
        iscol=0;
        if(isfield(data,N_('_ArrayOrder_')))
            arrayorder=data(j).(N_('_ArrayOrder_'));
            if(~isempty(arrayorder) && (arrayorder(1)=='c' || arrayorder(1)=='C'))
                iscol=1;
            end
        end
        if(isfield(data,N_('_ArrayIsSparse_')) && isstruct(data(j).(N_('_ArrayIsSparse_'))) )
                data(j).(N_('_ArrayIsSparse_'))=jdatadecode(data(j).(N_('_ArrayIsSparse_')),opt);
        end
        if(isfield(data,N_('_ArrayIsSparse_')) && data(j).(N_('_ArrayIsSparse_')))
                if(isfield(data,N_('_ArraySize_')))
                    if(isstruct(data(j).(N_('_ArraySize_'))))
                        data(j).(N_('_ArraySize_'))=jdatadecode(data(j).(N_('_ArraySize_')),opt);
                    end
                    dim=data(j).(N_('_ArraySize_'))(:)';
                    if(iscell(dim))
                        dim=cell2mat(dim);
                    end
                    dim=double(dim);
		    if(length(dim)==1)
			dim=[1 dim];
		    end
                    if(iscpx)
                        ndata(end-1,:)=complex(ndata(end-1,:),ndata(end,:));
                    end
                    if isempty(ndata)
                        % All-zeros sparse
                        ndata=sparse(dim(1),prod(dim(2:end)));
                    elseif dim(1)==1
                        % Sparse row vector
                        ndata=sparse(1,ndata(1,:),ndata(2,:),dim(1),prod(dim(2:end)));
                    elseif dim(2)==1
                        % Sparse column vector
                        ndata=sparse(ndata(1,:),1,ndata(2,:),dim(1),prod(dim(2:end)));
                    else
                        % Generic sparse array.
                        ndata=sparse(ndata(1,:),ndata(2,:),ndata(3,:),dim(1),prod(dim(2:end)));
                    end
                else
                    if(iscpx && size(ndata,2)==4)
                        ndata(3,:)=complex(ndata(3,:),ndata(4,:));
                    end
                    ndata=sparse(ndata(1,:),ndata(2,:),ndata(3,:));
                end
        elseif(isfield(data,N_('_ArrayShape_')))
                if(isstruct(data(j).(N_('_ArrayShape_'))))
                    data(j).(N_('_ArrayShape_'))=jdatadecode(data(j).(N_('_ArrayShape_')),opt);
                end
                if(iscpx)
                    if(size(ndata,1)==2)
                        dim=size(ndata);
                        dim(end+1)=1;
                        arraydata=reshape(complex(ndata(1,:),ndata(2,:)),dim(2:end));
                    else
                        error('The first dimension must be 2 for complex-valued arrays');
                    end
                else
                    arraydata=data.(N_('_ArrayData_'));
                end
                shapeid=data.(N_('_ArrayShape_'));
                if(isfield(data,N_('_ArrayZipSize_')))
                        datasize=data.(N_('_ArrayZipSize_'));
                        if(iscell(datasize))
                            datasize=cell2mat(datasize);
                        end
                        datasize=double(datasize);
                        if(iscpx)
                            datasize=datasize(2:end);
                        end
                else
                        datasize=size(arraydata);
                end
                if(isstruct(data(j).(N_('_ArraySize_'))))
                    data(j).(N_('_ArraySize_'))=jdatadecode(data(j).(N_('_ArraySize_')),opt);
                end
                arraysize=data.(N_('_ArraySize_'));

                if(iscell(arraysize))
                    arraysize=cell2mat(arraysize);
                end
                arraysize=double(arraysize);
                if(ischar(shapeid))
                        shapeid={shapeid};
                end
                arraydata=double(arraydata).';
                if(strcmpi(shapeid{1},'diag'))
                        ndata=spdiags(arraydata(:),0,arraysize(1),arraysize(2));
                elseif(strcmpi(shapeid{1},'upper') || strcmpi(shapeid{1},'uppersymm'))
                        ndata=zeros(arraysize);
                        ndata(triu(true(size(ndata)))')=arraydata(:);
                        if(strcmpi(shapeid{1},'uppersymm'))
                            ndata(triu(true(size(ndata))))=arraydata(:);
                        end
                        ndata=ndata.';
                elseif(strcmpi(shapeid{1},'lower') || strcmpi(shapeid{1},'lowersymm'))
                        ndata=zeros(arraysize);
                        ndata(tril(true(size(ndata)))')=arraydata(:);
                        if(strcmpi(shapeid{1},'lowersymm'))
                            ndata(tril(true(size(ndata))))=arraydata(:);
                        end
                        ndata=ndata.';
                elseif(strcmpi(shapeid{1},'upperband') || strcmpi(shapeid{1},'uppersymmband'))
                        if(length(shapeid)>1 && isvector(arraydata))
                            datasize=double([shapeid{2}+1, prod(datasize)/(shapeid{2}+1)]);
                        end
                        ndata=spdiags(reshape(arraydata,min(arraysize),datasize(1)),-datasize(1)+1:0,arraysize(2),arraysize(1)).';
                        if(strcmpi(shapeid{1},'uppersymmband'))
                            diagonal=diag(ndata);
                            ndata=ndata+ndata.';
                            ndata(1:arraysize(1)+1:end)=diagonal;
                        end
                elseif(strcmpi(shapeid{1},'lowerband') || strcmpi(shapeid{1},'lowersymmband'))
                        if(length(shapeid)>1 && isvector(arraydata))
                            datasize=double([shapeid{2}+1, prod(datasize)/(shapeid{2}+1)]);
                        end
                        ndata=spdiags(reshape(arraydata,min(arraysize),datasize(1)),0:datasize(1)-1,arraysize(2),arraysize(1)).';
                        if(strcmpi(shapeid{1},'lowersymmband'))
                            diagonal=diag(ndata);
                            ndata=ndata+ndata.';
                            ndata(1:arraysize(1)+1:end)=diagonal;
                        end
                elseif(strcmpi(shapeid{1},'band'))
                        if(length(shapeid)>1 && isvector(arraydata))
                            datasize=double([shapeid{2}+shapeid{3}+1, prod(datasize)/(shapeid{2}+shapeid{3}+1)]);
                        end
                        ndata=spdiags(reshape(arraydata,min(arraysize),datasize(1)),double(shapeid{2}):-1:-double(shapeid{3}),arraysize(1),arraysize(2));
                elseif(strcmpi(shapeid{1},'toeplitz'))
                        arraydata=reshape(arraydata,flipud(datasize(:))');
                        ndata=toeplitz(arraydata(1:arraysize(1),2),arraydata(1:arraysize(2),1));
                end
                if(opt.fullarrayshape && issparse(ndata))
                        ndata=cast(full(ndata),data(j).(N_('_ArrayType_')));
                end
        elseif(isfield(data,N_('_ArraySize_')))
            if(isstruct(data(j).(N_('_ArraySize_'))))
                data(j).(N_('_ArraySize_'))=jdatadecode(data(j).(N_('_ArraySize_')),opt);
            end
            if(iscpx)
                ndata=complex(ndata(1,:),ndata(2,:));
            end
            if(format>1.9 && iscol==0)
                data(j).(N_('_ArraySize_'))=data(j).(N_('_ArraySize_'))(end:-1:1);
            end
            dims=data(j).(N_('_ArraySize_'))(:)';
            if(iscell(dims))
                dims=cell2mat(dims);
            end
	    if(length(dims)==1)
		dims=[1 dims];
	    end
            ndata=reshape(ndata(:),dims(:)');
            if(format>1.9 && iscol==0)
                ndata=permute(ndata,ndims(ndata):-1:1);
            end
        end
        newdata{j}=ndata;
      end
      if(len==1)
          newdata=newdata{1};
      end
    end

    %% handle table data
    if(isfield(data,N_('_TableRecords_')))
        newdata=cell(len,1);
        for j=1:len
            ndata=data(j).(N_('_TableRecords_'));
            if(iscell(ndata))
                if(iscell(ndata{1}))
                    rownum=length(ndata);
                    colnum=length(ndata{1});
                    nd=cell(rownum, colnum);
                    for i1=1:rownum
                        for i2=1:colnum
                            nd{i1,i2}=ndata{i1}{i2};
                        end
                    end
                    newdata{j}=cell2table(nd);
                else
                    newdata{j}=cell2table(ndata);
                end
            else
                newdata{j}=array2table(ndata);
            end
            if(isfield(data(j),N_('_TableRows_'))&& ~isempty(data(j).(N_('_TableRows_'))))
                newdata{j}.Properties.RowNames=data(j).(N_('_TableRows_'))(:);
            end
            if(isfield(data(j),N_('_TableCols_')) && ~isempty(data(j).(N_('_TableCols_'))))
                newdata{j}.Properties.VariableNames=data(j).(N_('_TableCols_'));
            end
        end
        if(len==1)
            newdata=newdata{1};
        end
    end

    %% handle map data
    if(isfield(data,N_('_MapData_')))
        newdata=cell(len,1);
        for j=1:len
            key=cell(1,length(data(j).(N_('_MapData_'))));
            val=cell(size(key));
            for k=1:length(data(j).(N_('_MapData_')))
                key{k}=data(j).(N_('_MapData_')){k}{1};
                val{k}=jdatadecode(data(j).(N_('_MapData_')){k}{2},opt);
            end
            ndata=containers.Map(key,val);
            newdata{j}=ndata;
        end
        if(len==1)
            newdata=newdata{1};
        end
    end

    %% handle graph data
    if(isfield(data,N_('_GraphNodes_')) && exist('graph','file') && exist('digraph','file'))
        newdata=cell(len,1);
        isdirected=1;
        for j=1:len
            nodedata=data(j).(N_('_GraphNodes_'));
            if(isstruct(nodedata))
                nodetable=struct2table(nodedata);
            elseif(isa(nodedata,'containers.Map'))
                nodetable=[keys(nodedata);values(nodedata)];
                if(strcmp(nodedata.KeyType,'char'))
                    nodetable=table(nodetable(1,:)',nodetable(2,:)','VariableNames',{'Name','Data'});
                else
                    nodetable=table(nodetable(2,:)','VariableNames',{'Data'});
                end
            else
                nodetable=table;
            end

            if(isfield(data,N_('_GraphEdges_')))
                edgedata=data(j).(N_('_GraphEdges_'));
            elseif(isfield(data,N_('_GraphEdges0_')))
                edgedata=data(j).(N_('_GraphEdges0_'));
                isdirected=0;
            elseif(isfield(data,N_('_GraphMatrix_')))
                edgedata=jdatadecode(data(j).(N_('_GraphMatrix_')),varargin{:});
            end

            if(exist('edgedata','var'))
                if(iscell(edgedata))
                    endnodes=edgedata(:,1:2);
                    endnodes=reshape([endnodes{:}],size(edgedata,1),2);
                    weight=cell2mat(edgedata(:,3:end));
                    edgetable=table(endnodes,[weight.Weight]','VariableNames',{'EndNodes','Weight'});

                    if(isdirected)
                        newdata{j}=digraph(edgetable,nodetable);
                    else
                        newdata{j}=graph(edgetable,nodetable);
                    end
                elseif(ndims(edgedata)==2 && isstruct(nodetable))
                    newdata{j}=digraph(edgedata,fieldnames(nodetable));
                end
            end
        end
        if(len==1)
            newdata=newdata{1};
        end
    end

    %% handle bytestream and arbitrary matlab objects
    if(isfield(data,N_('_ByteStream_')) && isfield(data,N_('_DataInfo_')))
        newdata=cell(len,1);
        for j=1:len
            if(isfield(data(j).(N_('_DataInfo_')),'MATLABObjectClass'))
                if(needbase64)
                    newdata{j}=getArrayFromByteStream(base64decode(data(j).(N_('_ByteStream_'))));
                else
                    newdata{j}=getArrayFromByteStream(data(j).(N_('_ByteStream_')));
                end
            end
        end
        if(len==1)
            newdata=newdata{1};
        end
    end

    %% handle data link
    if(opt.maxlinklevel>0 && isfield(data,N_('_DataLink_')))
        if(ischar(data.(N_('_DataLink_'))))
            datalink=data.(N_('_DataLink_'));
            if(regexp(datalink,'\:\$'))
                ref=regexp(datalink, '^(?<proto>[a-zA-Z]+://)*(?<path>.+)(?<delim>\:)()*(?<jsonpath>(?<=:)\$\d*\.*.*)*', 'names');
            else
                ref=regexp(datalink, '^(?<proto>[a-zA-Z]+://)*(?<path>.+)(?<delim>\:)*(?<jsonpath>(?<=:)\$\d*\..*)*', 'names');
            end
            if(~isempty(ref.path))
                uripath=[ref.proto ref.path];
                [fpath, fname, fext]=fileparts(uripath);
                opt.maxlinklevel=opt.maxlinklevel-1;
                switch(lower(fext))
                    case {'.json','.jnii','.jdt','.jdat','.jmsh','.jnirs'}
                        newdata=loadjson(uripath, opt);
                    case {'.bjd' ,'.bnii','.jdb','.jbat','.bmsh','.bnirs', '.jamm'}
                        newdata=loadbj(uripath, opt);
                    case {'.ubj'}
                        newdata=loadubjson(uripath, opt);
                    case {'.msgpack'}
                        newdata=loadmsgpack(uripath, opt);
                    case {'.h5','.hdf5','.snirf'}  % this requires EasyH5 toolbox
                        newdata=loadh5(uripath, opt);
                    otherwise
                        % _DataLink_ url does not specify type, assuming JSON format
                        if(regexpi(datalink,'^\s*(http|https|ftp|file)://'))
                            newdata=loadjson(uripath, opt);
                        else
                            warning('_DataLink_ url is not supported')
                        end
                end
                if(~isempty(ref.jsonpath))
                    newdata=getfromjsonpath(newdata,ref.jsonpath);
                end
            end
        end
    end

    %% subfunctions
    function escaped=N_(str)
      escaped=[prefix str];
    end
end

author: fangq

repo: jsonlab

save_path: github-repos/MATLAB/fangq-jsonlab


text:
function dirText = dimFlip(dirText);
% flip direction text 'a <--> b' to read 'b <--> a'.
%
% Usage:
% dirText = dimFlip(dirText);
%
% ras, 11/05
I = findstr(dirText, ' <--> ');
if isempty(I), return; end
lhs = dirText(1:I);
rhs = dirText(I+5:end);
dirText = [rhs '<-->' lhs];
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
  function [xs, precon] = eml_osem(x, Gb, yi, ci, ri, varargin)
%|function xs = eml_osem(x, Gb, yi, ci, ri, [options])
%| E-ML-OSEM algorithm for image reconstruction from Poisson emission data
%| (ordered subsets expectation maximization)
%| model: Y_i ~ Poisson(c_i [G x]_i + r_i)
%|
%| in
%|	x	[np 1]		initial estimate
%|	Gb	[nd np]		Gblock object (see eml_osem_test.m)
%|	yi,ci,ri [nb na]	see em_fbp.m (for model too)
%|
%| option
%|	'niter', int		# iterations (default: 1+1)
%|	'isave', []		which iterations to archive (default: 'all')
%|	'pixmax', float		upper constraint for pixel values
%|	'precon [np {1|nblock}]	preconditioners; recommend: 'classic'
%|	'relax0' [1] or [2]	relax0 or (relax0, relax_rate)
%|	userfun @		user-defined function handle
%| out
%|	xs [np length(isave)]	updated image vectors each iteration
%|
%| Copyright 2002-3-14, Jeff Fessler, University of Michigan

if nargin < 3, ir_usage, end

% backward compatability with old usage
% eml_osem(x, Gb, yi, ci, ri, niter, pixmax, precon, relax0)
if length(varargin) && isnumeric(varargin{1})
	[xs, precon] = eml_osem_orig(x, Gb, yi, ci, ri, varargin{:});
return
end

% defaults
arg.niter = 1;
arg.isave = 'all';
arg.pixmax = inf;
arg.chat = false;
arg.relax0 = 1;
arg.relax_rate = 0;
arg.precon = 'classic';
arg.userfun = [];
arg = vararg_pair(arg, varargin);

arg.isave = iter_saver(arg.isave, arg.niter);

Gb = block_ob(Gb, 'ensure'); % make it a block object (if not already)
nblock = block_ob(Gb, 'n');
starts = subset_start(nblock);

if ~isvar('ci') || isempty(ci)
	ci = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end

if length(arg.relax0) == 1
	arg.relax_rate = 0;
elseif length(arg.relax0) == 2
	arg.relax_rate = arg.relax0(2);
	arg.relax0 = arg.relax0(1);
else
	error relax
end

eml_check(yi, ci, ri, 'os', nblock);
[nb na] = size(yi);

%
% precompute the preconditioners for classic OSEM
%
if streq(arg.precon, 'classic')
	precon = zeros(size(Gb,2), nblock);
	for iset=1:nblock
		ticker([mfilename ' : precon'], iset, nblock)
		istart = starts(iset);
		ia = istart:nblock:na;
		Asum = Gb{istart}' * col(ci(:,ia));
		Asum(Asum == 0) = Inf; % avoid divide by 0
		precon(:, iset) = 1 ./ Asum;
	end, clear Asum

%
% This 'fast' approach uses the same preconditioner for each subset
% which helps convergence (if diminishing relaxation is used).
% Of course, convergence is not really essential in the unregularized case,
% but this saves a bit of memory by using just one preconditioner, which
% ought to work when the subsets are reasonably balanced.
% But, if diminishing relaxation is not used, then this can cause problems
% because pixels can get stuck at zero. So I no longer recommend this.
%
elseif streq(arg.precon, 'fast')
	warning 'Using fast preconditioner rather than classic OSEM'
	precon = Gb' * ci(:); % complete backprojection
	precon(precon == 0) = Inf; % avoid divide by 0
	precon = nblock ./ precon;

else % user-supplied precon
	precon = arg.precon;
	if ncol(precon) ~= nblock && ncol(precon) ~= 1
		error 'precon columns must be 1 or nblock'
	end
end

%
% loop over iterations
%
xs = zeros(numel(x), length(arg.isave));
if any(x <= 0), error 'need x > 0', end
x = min(x, arg.pixmax);
if any(arg.isave == 0)
	xs(:,find(arg.isave == 0)) = x;
end

for iter = 1:arg.niter
	relax = arg.relax0 / (1 + arg.relax_rate * (iter-2));

	%
	% loop over subsets
	%
	for iset=1:nblock
		ticker(mfilename, [iter iset], [arg.niter nblock])

		istart = starts(iset);
		ia = istart:nblock:na;

		% predicted measurements
		yp = reshape(Gb{istart} * x, nb, length(ia));
		yp = ci(:,ia) .* yp + ri(:, ia);
		yp(yp == 0) = inf;	% avoids /0 error

		pre = precon(:,min(iset,ncol(precon)));	% 1 or iset

		dhi = ci(:,ia) .* (yi(:,ia) ./ yp - 1);
		grad = Gb{istart}' * dhi(:);

		x = x + relax * (x .* pre) .* grad;

		% fix: implement ahn's shift trick? no, just use classic.
		nneg = sum(x < 0);
		if nneg, printm('oh no! %d negatives', nneg), minmax(x), end

		x = max(x,0);	% caution: if x becomes zero, it is always zero!
		x = min(x,arg.pixmax);
	end

	if arg.chat, printm('Range %g %g', min(x), max(x)), end
	if any(arg.isave == iter)
		xs(:,find(arg.isave == iter)) = x;
	end

	if ~isempty(arg.userfun)
		feval(arg.userfun);
	end
end

%
% eml_osem_orig()
% arguments: pixmax, precon, relax0
%
function [xs, precon] = eml_osem_orig(x, Gb, yi, ci, ri, niter, varargin)
arg = {'niter', niter};
opt = {'pixmax', 'precon', 'relax0'};
for ii=1:length(varargin)
	arg = {arg{:}, opt{ii}, varargin{ii}};
end
[xs, precon] = eml_osem(x, Gb, yi, ci, ri, arg{:});

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [rVectECI, vVectECI] = getECIVectFromCR3BPVect(ut, rVectCR3BP, vVectCR3BP, secBodyInfo, celBodyData)

    priBodyInfo = secBodyInfo.getParBodyInfo(celBodyData);
    gmu = priBodyInfo.gm;

    rotMatEci2Cr3bp = getCR3BPRotMat(ut, secBodyInfo, celBodyData);
    zVectCR3BP = rotMatEci2Cr3bp(:,3);

    rotMatCr3bp2Eci = rotMatEci2Cr3bp;
    rVectECI = rotMatCr3bp2Eci * rVectCR3BP;

    secPeriod = computePeriod(secBodyInfo.sma, gmu);
    rotRateRadSec = 2*pi/secPeriod;
    omegaRI = zVectCR3BP * rotRateRadSec;

    vVectECI = rotMatCr3bp2Eci*(vVectCR3BP + crossARH(omegaRI, rVectCR3BP));
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function quadrule_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests GAUSS_LEGENDRE_RULE_COMPUTE and GAUSS_LEGENDRE_RULE_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  GAUSS_LEGENDRE_RULE_COMPUTE computes\n' );
  fprintf ( 1, '  a Clenshaw-Curtis quadrature rule.\n' );
  fprintf ( 1, '  GAUSS_LEGENDRE_RULE_SET sets\n' );
  fprintf ( 1, '  a Clenshaw-Curtis quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compare:\n' );
  fprintf ( 1, '    (X1,W1) from GAUSS_LEGENDRE_RULE_SET\n' );
  fprintf ( 1, '    (X2,W2) from GAUSS_LEGENDRE_RULE_COMPUTE\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '     Order        W1            W2            X1            X2\n' );
  fprintf ( 1, '\n' );

  for order = 1 : 3 : 10

    [ x1, w1 ] = gauss_legendre_rule_set ( order );
    [ x2, w2 ] = gauss_legendre_rule_compute ( order );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %8d\n', order );

    for i = 1 : order
      fprintf ( 1, '            %12f  %12f  %12f  %12f\n', ...
        w1(i), w2(i), x1(i), x2(i) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [perimeterEdges,eulerCondition] = findLegalPerimeters2(mesh,perimDist)
%
%  [perimeterEdges,eulerCondition] = findLegalPerimeters2(mesh,perimDist)
%
% Given a mesh structure (nodes, edges, distances from start point) and a
% threshold distance this routine will return a list of edges that
% constitute separate perimeters at a distance of perimDist from the start
% point.
%
% Note that because of intrinsic curvature, there can be more than one
% perimeter at the required distance. This routine makes sure that it
% returns separate perimeters (ones with no common nodes).

% Find perims with simple thold
insideNodes=find(mesh.dist<=perimDist);
insideNodes=insideNodes(:);

insideNodes = removeHangingNodes(mesh,insideNodes); % Cleans up the mesh
[perimeterEdges,eulerCondition]=findGroupPerimeter(mesh,insideNodes); % Find perimeter(s) - there may be more than 1
badPerimNodes = findBadPerimNodes(mesh,perimeterEdges); % See if some perimeters are joined up.

numBadNodes=9999999; % Hope we don't get more than this...

while (numBadNodes>0)
    [perimeterEdges,eulerCondition]=findGroupPerimeter(mesh,insideNodes);
	length(perimeterEdges);
	length(unique(perimeterEdges,'rows'));
    fprintf('Euler number=%d\n',eulerCondition);

    badPerimNodes=findBadPerimNodes(mesh,perimeterEdges);
    numBadNodes=length(badPerimNodes);
    fprintf('There are %d bad perim nodes.\n',numBadNodes);

    if(numBadNodes)
		[insideNodes]=correctBadNodes(mesh,insideNodes,badPerimNodes); % Splits up joined perimeters
		insideNodes=removeHangingNodes(mesh,insideNodes); % Cleans up mesh again
    end
end

return;


text:
function write_neuralynx_nse(filename, nse)

% WRITE_NEURALYNX_NSE writes spike timestamps and waveforms to a NSE file
% The input data should be scaled in uV.
%
% Use as
%   write_neuralynx_nse(filename, nse)
%
% See also READ_NEURALYNX_NSE

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

if ~isa(nse.TimeStamp, 'uint64')
  ft_error('timestamps should be uint64');
end

% convert the data from uV into V
nse.dat = nse.dat * 1e-6;
% scale the data and convert to 16 bits,
% this has to be done prior to writing the header
ADMaxValue = double(intmax('int16'));
ADMaxVolts = max(abs(nse.dat(:)));
if ADMaxVolts>0
  ADBitVolts = ADMaxVolts / ADMaxValue;
else
  ADBitVolts = 1;
end
nse.dat = int16(nse.dat / ADBitVolts);
% update the header with the calibration values
nse.hdr.ADBitVolts = ADBitVolts;
nse.hdr.ADMaxValue = ADMaxValue;

% construct the header
buf  = [];
buf  = [buf sprintf('######## Neuralynx Data File Header\r\n')];
buf  = [buf sprintf('## File Name: %s\r\n', filename)];
buf  = [buf sprintf('## Time Opened: (m/d/y): %s  At Time: %s\r\n', datestr(clock, 'mm/dd/yy'), datestr(clock, 'HH:MM:SS'))];
f = fieldnames(nse.hdr);
for i=1:length(f)
  v = getfield(nse.hdr, f{i});
  switch class(v)
    case 'char'
      buf = [buf sprintf('-%s\t%s\r\n', f{i}, v)];
    case 'double'
      buf = [buf sprintf('-%s\t%g\r\n', f{i}, v)];
    otherwise
      ft_error('unknown class in writing header');
  end
end

% pad the rest of the header with zeros
buf((end+1):16384) = 0;

% write the header to file
fid  = fopen(filename, 'wb', 'ieee-le');
fwrite(fid, buf);

% The format of a single electrode record is
%   int64 TimeStamp
%   int32 ScNumber
%   int32 CellNumber
%   int32 Param[0] ... Param[7]
%   int16 Samp[0] ... int16 Samp[31]

for i=1:size(nse.dat,2)
  % write a single continuous data record
  fwrite(fid, nse.TimeStamp(i) , 'int64');
  fwrite(fid, nse.ScNumber(i)  , 'int32');
  fwrite(fid, nse.CellNumber(i), 'int32');
  fwrite(fid, nse.Param(:,i)   , 'int32');
  fwrite(fid, nse.dat(:,i)     , 'int16');
end

% close the file
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [] = doc_examples
% File:      doc_examples.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.14 - 2012.06.16
% Language:  MATLAB R2012a
% Purpose:   documentation examples of fig2u3d
% Copyright: Ioannis Filippidis, 2012-

%fname = two_variable_func;
%fname = three_variable_level_set;
%fname = molecule_structure;
%fname = robotic_manipulator;
fname = fluid_flow;
fig2u3d(gca, fname, '-pdf')

function [fname] = two_variable_func
%f(x,y) surface, gradient, contours
ax = newax;
hold(ax, 'on')

dom = 10 *[-1, 1, -1, 1];
res = [30, 40];

q = domain2vec(dom, res);
x = q(1, :);
y = q(2, :);
f = sin(sqrt(0.5 .*x.^2 +y.^2) ) +2;
v = bsxfun(@times, cos(sqrt(0.5 .*x.^2 +y.^2) ) .*(0.5 .*x.^2 +y.^2).^(-0.5), [0.5 .*x; y] );

vsurf(ax, q, f, res)
vcontour(ax, q, f, res)
quivermd(ax, q, v)
axis tight
view(ax, 3)

fname = 'two_var_func';

function [fname] = three_variable_level_set
% Elliptic supercylide level set surface, gradient and 2D field section
example_supercyclide
fname = 'example_supercyclide';

function [fname] = molecule_structure
draw_crystal_lattice
fname = 'diamond';

function [fname] = robotic_manipulator
% requires the Robotics Toolbox by Peter I. Corke
mdl_puma560
p560.plot(zeros(1,6) )
fname = 'puma560';

function [fname] = fluid_flow
Wind
fname = 'wind';

author: Sable

repo: mcbench-benchmarks


text:
function [Y, X] = minandmax2( f )
%MINANDMAX2   Find global minimum and maximum of a SPHEREFUN.
%   M = minandmax2(F) returns the minimum and maximum value of a SPHEREFUN
%   over its domain. M is a vector of length 2 such that
%   M(1) = min(f(lambda,theta)) and M(2) = max(f(lambda,theta)).
%
%   [M, LOC] = minandmax2(F) also returns the position of the minimum and
%   maximum. For example,
%
%       F(LOC(1,1),LOC(1,2)) = M(1)  and  F(LOC(2,1),LOC(2,2)) = M(2)
%
% See also SPHEREFUN/MAX2, SPHEREFUN/MIN2, SPHEREFUN/NORM.

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

% check for empty SPHEREFUN.
if ( isempty( f ) )
    Y = [];
    X = [];
    return
end

[C, D, R] = cdr(f);
f = C*D*R.';
% Maximum possible sample matrix size:
maxsize = 4e3;

% Is the function the zero function?
if ( iszero( f )  )
    dom = f.domain;
    X = [ (dom(2) + dom(1))/2 (dom(4) + dom(3))/2 ];
    X = [ X ; X ];
    Y = [0 ; 0];
    return;
end

% Extract low rank representation:
frows = f.rows;
fcols = f.cols;
piv = f.pivotValues;
dom = f.domain;

% Share out scaling:
sgn = sign( piv ).';
sq = 1 ./ sqrt( abs( piv ) );
frows = frows * diag( sq.'.*sgn );
fcols = fcols * diag( sq );

if ( length(f) == 1 ) % Rank-1 is easy:
    % We can find it from taking maximum and minimum in x and y direction.

    % Find minandmax of rows and columns:
    [yr, xr] = minandmax( frows );
    [yc, xc] = minandmax( fcols );
    % All possible combinations:
    vv = [yr(1)*yc(1), yr(1)*yc(2), yr(2)*yc(1), yr(2)*yc(2)];

    [Y(2), indmx] = max( vv );
    [Y(1), indmn] = min( vv );

    % Work out the location of the maximum.
    X = zeros(2);
    X(1,1) = xr(2);
    X(1,2) = xc(2);

    if ( indmn <= 2 )
        X(1,1) = xr(1);
    end
    if ( mod(indmn,2) == 1 )
        X(1,2) = xc(1);
    end
    X(2,1) = xr(2);
    X(2,2) = xc(2);

    if ( indmx <= 2 )
        X(2,1) = xr(1);
    end

    if ( mod(indmx,2) == 1 ),
        X(2,2) = xc(1);
    end

elseif ( length(f) <= maxsize )

    % We seek a fast initial guess. So we first truncate the SPHEREFUN.
    ypts = [trigpts(length(fcols), fcols.domain);pi];
    xpts = [trigpts(length(frows), frows.domain);pi];
    cvals = feval(fcols, ypts);
    rvals = feval(frows, xpts);

    A = cvals*rvals.';
    % Maximum entry in discretisation.
    [ignored, ind] = min( A(:) ); %#ok<ASGLU>
    [row, col] = ind2sub(size(A), ind);
    X(1,1) = xpts( col );
    X(1,2) = ypts( row );
    Y(1) = feval( f, X(1,1), X(1,2) );
    % Minimum entry in discretisation.
    [ignored, ind] = max(A(:)); %#ok<ASGLU>
    [row, col] = ind2sub(size(A), ind);
    X(2,1) = xpts(col);
    X(2,2) = ypts(row);
    Y(2) = feval(f, X(2,1), X(2,2));

    % Get more digits with optimisation algorithms.
    lb = [ dom(1) ; dom(3) ];
    ub = [ dom(2) ; dom(4) ];

    try

        % If the optimization toolbox is available then use it to get a better maximum.

        warnstate = warning;
        warning('off'); %#ok<WNOFF> % Disable verbose warnings from fmincon.
        options = optimset('Display', 'none', 'TolFun', eps, 'TolX', eps, ...
            'algorithm', 'active-set');
        [mn, Y(1)] = fmincon(@(x,y) feval(f, x(1), x(2)), X(1, :), ...
            [], [], [], [], lb, ub, [], options);
        [mx, Y(2)] = fmincon(@(x) -feval(f, x(1), x(2)), X(2,:), ...
            [], [], [], [], lb, ub, [], options);
        Y(2) = -Y(2);
        X(1,:) = mn;
        X(2,:) = mx;
        warning(warnstate);

    catch

        try
            % Try converting to an unconstrained problem and using built-in solver.

            % Maps from [-1, 1] to [dom(1:2)] and [dom(3:4)], respectively.
            map1 = bndfun.createMap(dom(1:2));
            map2 = bndfun.createMap(dom(3:4));
            % Unconstrained initial guesses:
            Z(:,1) = asin(map1.Inv(X(:,1)));
            Z(:,2) = asin(map2.Inv(X(:,2)));
            % Maps from R to [dom(1), dom(2)] and [dom(3), dom(4)], respectively.
            map1 = @(x) map1.For(sin(x));
            map2 = @(x) map2.For(sin(x));
            % Set options:
            options = optimset('Display', 'off', 'TolFun', eps, 'TolX', eps);
            warnstate = warning;
            warning('off'); %#ok<WNOFF> % Disable verbose warnings from fminsearch.
            f_mapped = @(x) feval(f, map1(x(1)), map2(x(2)));
            [mn, Y(1)] = fminsearch(@(x) f_mapped(x), Z(1, :), options);
            [mx, Y(2)] = fminsearch(@(x) -f_mapped(x), Z(2, :), options);
            Y(2) = -Y(2);
            X(1:2,1) = map1([mn(1) ; mx(1)]);
            X(1:2,2) = map2([mn(2) ; mx(2)]);
            warning(warnstate);

        catch

            % Nothing is going to work so initial guesses will have to do.

        end

    end

elseif ( length(f) >= maxsize )

    error('CHEBFUN:SPHEREFUN:minandmax2:length', 'Rank is too large.');

end

% If the location of the max/min is outside of [-pi,pi]^2, which can happen
% because Spherefun is based on the double Fourier sphere method, then
% translate back to [0,pi]x[-pi,pi] using BMC symmetry:

% check minimum:
if ( X(1,1) < 0 )
    X(1,1) = X(1,1)+pi;
    X(1,2) = -X(1,2);
end

% check maximum:
if ( X(2,1) < 0 )
    X(2,1) = X(2,1)+pi;
    X(2,2) = -X(2,2);
end

end

%%%
% Use the approach below when bivariate rootfinding is fully implemented.
%%%
% % Use bivariate rootfinding to find all the local extrema:
% F = gradient( f );
% r = roots( F );
% if ( ~isempty( r ) )
%     [inMax, idInMax] = max( feval(f, r(:,1), r(:,2) ) );
%     [inMin, idInMin] = min( feval(f, r(:,2), r(:,2) ) );
% else
%     inMax = inf;   % max and min must occur on boundary.
%     inMin = inf;
% end
%
% % Search along boundary:
% dom = f.domain;
% left = feval(f, dom(1), ':');
% right = feval(f, dom(2), ':');
% down = feval(f, ':', dom(3));
% up = feval(f, ':', dom(4));
% [Yleft, Xleft] = minandmax( left );
% [Yright, Xright] = minandmax( right );
% [Yup, Xup] = minandmax( up );
% [Ydown, Xdown] = minandmax( down );
%
% % Store Min/Max location for later:
% BcMinLocations = [ Xleft(1,:) ; Xright(1,:) ; Xup(1,:) ; Xdown(1,:) ];
% BcMaxLocations = [ Xleft(2,:) ; Xright(2,:) ; Xup(2,:) ; Xdown(2,:) ];
%
% [ BcMax, idBcMax ] = max( [ Yleft(2), Yright(2), Yup(2), Ydown(2) ].' );
% [ BcMin, idBcMin ] = min( [ Yleft(1), Yright(1), Yup(1), Ydown(1) ].' );
%
% % What is the global min and max?:
% [Ymax, inOrOutMax] = max( [inMax, BcMax].' );
% [Ymin, inOrOutMin] = min( [inMin, BcMin].' );
% Y = [Ymin, Ymax];
% X = zeros(2);
%
% % Unravel to find locations:
% if ( inOrOutMin == 1 )
%     X(1,:) = [ r(idInMin,1), r(idInMin,2) ];
% else
%     X(1,:) = BcMinLocations(idBcMin, :);
% end
% if ( inOrOutMax == 1 )
%     X(2,:) = [ r(idInMax,1), r(idInMax,2) ];
% else
%     X(2,:) = BcMaxLocations(idBcMax, :);
% end
%
% end

repo: chebfun


text:
% TIME = EPOCH2DATENUM(EPOCH) translates Unix epoch date format into Matlab
% numeric date format.

% Author: Michael West, Modified by Glenn Thompson
% $Date$
% $Revision$

function time = epoch2datenum(epoch)
    if admin.antelope_exists()
        % check ANTELOPE installation if this line segfaults
        time = datenum(epoch2str(epoch,'%m %d %Y %H %M %S.%s'),'mm dd yyyy HH MM SS.FFF');
    else % does not handle leap seconds
        time = datenum(1970,1,1,0,0,epoch);
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function accuracy = rksr_dsk_classifier(TrainSet, TestSet, options)
% Riemannian kernelized sparse representation classification (R-KSRC) with discriminative Stein kernel algorithm
%
% Inputs:
%       TrainSet            train sets of size dxn, where d is dimension and n is number of sets
%       TestSet             test sets of size dxn, where d is dimension and n is number of sets
%       options             options
% Output:
%       accuracy            classification accurary
%
% References:
%       M. Harandi, R. Hartley, B. Lovell and C. Sanderson,
%       "Sparse coding on symmetric positive definite manifolds using bregman divergences,"
%       IEEE Transactions on Neural Networks and Learning Systems, vol.27, no.6, pp.1294-1306, 2016.
%
%       M. Harandi, C. Sanderson, R. Hartley and B. Lovell,
%       "Sparse coding and dictionary learning for symmetric positive definite matrices: a kernel approach,"
%       European Conference on Computer Vision (ECCV), 2012.
%
%       J. Zhang, L. Wang. L. Zhou, and W. Li,
%       "Learning discriminative Stein kernel for SPD matrices and its applications,"
%       IEEE Transactions on Neural Networks and Learning Systems, vol.27, no.5, pp.1020-1033, 2015.
%
% Originally created by Mehrtash Harandi (mehrtash.harandi at gmail dot com).
% Originally created by J. Zhang (jz163@uowmail.edu.au).
% Modified by H. Kasai on July 10, 2017.

    % retrieve dimension of the SPD matrices
    dim = size(TrainSet.X_cov, 1);
    test_num = length(TestSet.y);
    class_num = length(unique(TestSet.y));

    % calculate eigen decomposition
    train_decomp = Decomposite_eig_new(TrainSet);
    test_decomp = Decomposite_eig_new(TestSet);

    % learn alpha of kernel
    if ~options.original_alpha
        initial_alpha = 1*ones(1,dim); % the initial alpha corresponding to the original Stein kernel
        LB = 0.01*initial_alpha;

        fmincon_opt = optimset('Algorithm', 'interior-point'); % run interior-point algorithm

        if options.verbose
            fmincon_opt.Display = 'iter';
        else
            fmincon_opt.Display = 'off';
        end
        fmincon_opt.MaxIter = 100;
        fmincon_opt.TolFun = 1e-5;
        %tic
        optimal_alpha = fmincon(@(alpha) objfun_ff_new(alpha,TrainSet.y,train_decomp, options.lambda,initial_alpha,options.obj_method,options.theta),initial_alpha,[],[],[],[],LB,[],[],fmincon_opt);
        %toc
    else
        optimal_alpha = ones(1,dim);
    end

    % compute the Stein divergence with the obtained adjustment parameter optimal_alpha
    S_test          = EigComp2SD_power_new(train_decomp, test_decomp, optimal_alpha);
    S_train         = EigComp2SD_power_new(train_decomp, train_decomp, optimal_alpha);
    train_kernel    = exp(-1 * options.theta * S_train);
    test_kernel     = exp(-1 * options.theta * S_test);

    % normalize dictionary
    [KD, ~] = data_normalization(train_kernel, [], 'std');
    [KX, ~] = data_normalization(test_kernel, [], 'std');

    [KD_U, KD_D, ~] = svd(KD);
    A = diag(sqrt(diag(KD_D))) * KD_U';
    D_Inv = KD_U * diag(1./sqrt(diag(KD_D)));
    KX = D_Inv' * KX;

    % perform lasso
    param.lambda = options.lambda;
    param.lambda2 =  0;
    param.mode = 2;
    scX = full(mexLasso(KX,A,param));

    % prepare class array
    classes = unique(TrainSet.y);
    % prepare predicted label array
    identity = zeros(1, test_num);

    for i = 1 : test_num
        % prepare residual array
        residuals = zeros(1, class_num);

        % calculate residual for each class
        for j = 1 : class_num
            idx = find(TrainSet.y == classes(j));
            %residuals(j) = norm(KX(:, i) - A(:,idx)*scX(idx, i))/sum(scX(idx, i) .* scX(idx, i));
            residuals(j) = norm(KX(:, i) - A(:,idx)*scX(idx, i));

            if strcmp(options.mode, 'src')
                residuals(j) = norm(KX(:, i) - A(:,idx)*scX(idx, i));
            elseif strcmp(options.mode, 'ip_linear')
                scX_i = scX(idx, i);
                %residuals(j) = sum(abs(scX_i));
                residuals(j) = sum(scX_i);
            elseif strcmp(options.mode, 'ip_max')
                scX_i = scX(idx, i);
                residuals(j) = max(abs(scX_i));
            end
        end

        % calculate the predicted label
        [~, label] = min(residuals);
        if strcmp(options.mode, 'src')
            [~, label] = min(residuals);
        elseif strcmp(options.mode, 'ip_linear') || strcmp(options.mode, 'ip_max')
            [~, label] = max(residuals);
        end
        identity(i) = label;

        if options.verbose
            correct = (label == TestSet.y(1, i));
            fprintf('# RSR-DSK-%s (with %s metric): test:%03d, predict class: %03d --> ground truth :%03d (%d)\n', options.obj_method, options.mode, i, label, TestSet.y(1, i), correct);
        end
    end

    % calculate accuracy
    correct_num = sum(identity == TestSet.y);
    accuracy = correct_num/test_num;
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% Profiler extension for Kernel Least-Mean-Square algorithm
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef klms_profiler < klms

    properties (GetAccess = 'public', SetAccess = 'private')
        elapsed = 0; % elapsed time
        prev_dict_size = 0; % previous dictionary size for growth check
    end

    methods

        function kaf = klms_profiler(parameters) % constructor
            if nargin<1, parameters = struct(); end
            kaf = kaf@klms(parameters);
        end

        function flops = lastflops(kaf) % flops for last iteration
            m = size(kaf.dict,1);
            if kaf.prev_dict_size < m % growing
                m1 = m - 1;
                floptions = struct(...
                    'sum', m1 - 1, ...
                    'mult', m1 + 1, ...
                    sprintf('%s_kernel',kaf.kerneltype), [m1, 1, size(kaf.dict,2)]);

                flops = kflops(floptions);
            else
                flops = 0;
            end
        end

        %% flops breakdown

        % k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
        % kernel: m1

        % y_est = k'*kaf.alpha;
        % sum: m1 - 1
        % mult: m1

        % err = y - y_est;
        % sum: 1

        % kaf.alpha = [kaf.alpha; kaf.mu*err];
        % mult: 1

        %%

        function train_profiled(kaf,x,y)
            kaf.prev_dict_size = size(kaf.dict,1);
            t1 = tic;
            kaf.train(x,y);
            t2 = toc(t1);
            kaf.elapsed = kaf.elapsed + t2;
        end

        function bytes = lastbytes(kaf) % bytes used in last iteration
            m = size(kaf.dict,1);
            bytes = 8*(m + m*size(kaf.dict,2)); % 8 bytes for double precision
            % alpha, dict
        end

    end
end

repo: kafbox