text:
function yp = p33_fun ( neqn, t, y )

%*****************************************************************************80
%
%% P33_FUN evaluates the function for problem P33.
%
%  Discussion:
%
%    2 equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wayne Enright, John Pryce,
%    Algorithm 648,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 1, pages 28-34.
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the derivative
%    function.
%
%    Output, real YP(NEQN), the value of the derivative function.
%
  yp = zeros ( neqn, 1 );

  delta = p33_param ( 'GET', 'DELTA' );
  yp(1) = y(2);
  yp(2) = delta * ( 1.0 - y(1).^2 ) * y(2) - y(1);

  return
end


text:
function meta = metaAddData(fieldName, value, viewType);
%
%  meta = metaAddData(fieldName, value, [viewType='Inplane']);
%
% Add metadata to a session.
%
% ras, 11/02/06: wrote it.
if nargin<2, error('Not enough Input args.'); end
if notDefined('viewType'), viewType = 'Inplane'; end

switch lower(viewType)
    case 'inplane'

    case 'volume'

    case 'gray'

    case 'mr'

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function hole_point = p14_hole_point ( hole_index, m )

%*****************************************************************************80
%
%% P14_HOLE_POINT returns a point inside a given hole in problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer HOLE_INDEX, the index of the hole.
%
%    Input, integer M, the spatial dimension.
%
%    Output, real HOLE_POINT(M), a point in the hole
%
  hole_point(1:2) = [ 257.068, 278.225 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function disp(obj)
% Display test results
%
% disp(obj)
%
%   obj             MOxUnitTestReport instance.
%
% Notes:
%   - this function writes results to the file descriptor provided when
%     obj was instantiated, and with the provided verbosity level
%
% NNO 2015

    output_format='text';
    fprintf(obj.stream,'%s\n',getSummaryStr(obj,output_format));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef TP6 < PROBLEM
% <multi> <real> <large/none> <robust>
% Test problem for robust multi-objective optimization
% delta --- 0.05 --- Maximum disturbance degree
% H     ---   50 --- Number of disturbances

%------------------------------- Reference --------------------------------
% A. Gaspar-Cunha, J. Ferreira, and G. Recio, Evolutionary robustness
% analysis for multi-objective optimization: benchmark problems, Structural
% and Multidisciplinary Optimization, 2014, 49: 771-793.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        delta;      % Maximum disturbance degree
        H;          % Number of disturbances
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.delta,obj.H] = obj.ParameterSet(0.05,50);
            obj.M = 2;
            if isempty(obj.D); obj.D = 5; end
            obj.lower    = [0,-ones(1,obj.D-1)];
            obj.upper    = [1, ones(1,obj.D-1)];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(~,PopDec)
            PopObj(:,1) = PopDec(:,1);
            h = 1 - PopDec(:,1).^2;
            g = sum(10-10*cos(4*pi*(PopDec(:,2:end)))+PopDec(:,2:end).^2,2);
            S = 1./(0.2+PopDec(:,1)) + PopDec(:,1).^2;
            PopObj(:,2) = h + g.*S;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^2 ;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case {'Mean_IGD','Mean_HV','Worst_IGD','Worst_HV'}
                    score = feval(metName,Population,obj);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Perturb solutions multiple times
        function PopX = Perturb(obj,PopDec,N)
            if nargin < 3; N = obj.H; end
            Delta = repmat(obj.delta.*(obj.upper-obj.lower),N*size(PopDec,1),1);
            w     = UniformPoint(N,obj.D,'Latin');
            Dec   = 2*Delta.*w(reshape(repmat(1:end,size(PopDec,1),1),1,[]),:) + repmat(PopDec,N,1) - Delta;
            Dec   = obj.CalDec(Dec);
            PopX  = SOLUTION(Dec,obj.CalObj(Dec),obj.CalCon(Dec));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function rule_num = wandzura_rule_num ( )

%*****************************************************************************80
%
%% WANDZURA_RULE_NUM returns the number of Wandzura rules available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wandzura, Hong Xiao,
%    Symmetric Quadrature Rules on a Triangle,
%    Computers and Mathematics with Applications,
%    Volume 45, Number 12, June 2003, pages 1829-1840.
%
%  Parameters:
%
%    Output, integer RULE_NUM, the number of rules available.
%
  rule_num = 6;

  return
end


text:
function radians = degrees_to_radians ( degrees )

%*****************************************************************************80
%
%% DEGREES_TO_RADIANS converts an angle measure from degrees to radians.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real DEGREES, the angle measure in degrees.
%
%    Output, real RADIANS, the angle measure in radians.
%
  radians = ( degrees / 180.0 ) * 3.141592653589793;

  return
end


text:
% Fig. 6.42   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

alpha=logspace(-2,0,100);
N=size(alpha);
ialpha=ones(N)./alpha;
dum=(ialpha-ones(N))./(ialpha+ones(N));
phimax=asin(dum);
ff=180/pi;
phimax=ff*phimax;
semilogx(ialpha,phimax);
xlabel('1/\alpha');
ylabel('Maximum phase lead (deg)');
grid;
title('Fig. 6.53 Maximum phase increase for lead compensation.');

author: Sable

repo: mcbench-benchmarks


text:
 This Software computes the Robust Principal Component Analysis or Robust Singular Value Decomposition
 explained in the references:

 REFERENCES:
 De la Torre, F. and Black, M. J., Robust principal component analysis for computer vision. Int. Conf. on Computer Vision 2001, Vancouver.
 De la Torre, F. and Black, M. J., A Framework for Robust Subspace Learning. In preparation.

This is EXPERIMENTAL software, so with high probability is possible to have bugs. If you find any bug, please
send an e-mail to  ftorre@salleURL.edu.
Use at your own risk. No warranty is implied by this distribution.

There are three Matlab(it has been tested in version 5.3) programs:

rpca: An example program of how to run RPCA
rob_pca: Computes Robust PCA
weighted_pca: Computes Weighted PCA

data_2_days is a file containing the data gathered during 2 days with a static camera.
Beware! you need at least a 800M RAM memory to run the algorithm in Matlab with this data!
It takes several hours to converge!

If you do not have this memory space, just compute the algorithm with half of the data
(e.g. add this line after loading the data in rpca.m Data=Data(:,1:end/2);)

Link to download original source code + data: https://gofile.io/?c=QIsyrz

author: andrewssobral

repo: lrslibrary


text:
function [T1, M0] = Compute_M0_T1_OnSPGR(data, flipAngles, TR, b1Map, roi, verbose)
%Perform a simple linear-least squares data fit on variable flip angle SPGR data
%
% function [M0, T1] = Compute_M0_T1_OnSPGR(data, flipAngles, TR [, b1Map, roi, verbose])
% -----------------------------------------------------------
% INPUTS:
%   data: width x length x slices x flipAngles matrix
%   flipAngles: vector of flip angles (in degrees) corresponding to last dimension of 'data'
%   TR: Repetition time in seconds
%   b1Map: width x length x slices matrix containing relative flip angle
%         (i.e. if nominal alpha is 60 and measured alpha is 61, then b1Map = 61/60
%   roi: width x length x slices binary mask
%   verbose: logical - for debugging

if ndims(data)<3, data = permute(data(:),[2 3 4 1]); end
[nX, nY, nZ, nFlip] = size(data);
nVox = nX*nY*nZ;

if ~exist('b1Map', 'var') || isempty(b1Map)
    b1Map = ones(nX, nY, nZ);
end

if ~exist('roi', 'var') || isempty(roi)
    roi = true(nX, nY, nZ);
end

if ~exist('verbose', 'var')
    verbose = 0;
end

if length(b1Map) ~= length(data(:,:,:,1)), error('B1 size is different from data size'); end
if ~islogical(roi), roi = logical(roi); end

% Reshape data into 2D array so that future steps are simpler
data = reshape(data,[nVox nFlip])'; % Transpose because MATLAB is column-major
data = data(:, roi(:));
nVox = sum(roi(:));

% Large matrix with redundant values will make computations faster,
% but will also consume more memory
alpha = deg2rad(flipAngles);
if isrow(alpha), alpha = alpha'; end
alpha = repmat(alpha,[1 nVox]) .* repmat(b1Map(roi(:))', [nFlip 1]);

% Do the linear least squares fit and estimate M0 & T1
y = data ./ sin(alpha);
x = data ./ tan(alpha);
[fittedSlope, fittedIntercept] = LinLeastSquares(x,y);
estM0 = fittedIntercept ./ (1-fittedSlope);
estT1 = -TR./log(fittedSlope);

% Assign arbitrary M0 & T1 value if fitted value is unphysical.
% Might be better to set this to NaN so users can identify voxels where fit fails
failedFit = isnan(fittedSlope)      | fittedSlope<0     | ...
            isnan(fittedIntercept)  | fittedIntercept<0 ;
failedFitValue = NaN;
estM0(failedFit) = failedFitValue;
estT1(failedFit) = failedFitValue;

% Assign estimated M0 and T1 values to correct voxel
[T1, M0] = deal(zeros(nX,nY,nZ));
M0(roi(:)) = estM0;
T1(roi(:)) = estT1;
end % END OF Compute_M0_T1_OnSPGR

function [fittedSlope, fittedIntercept] = LinLeastSquares(x,y)
% Simple linear least squares fit
% Inputs:
%   - x and y, arrays of equal sizes
% The first dimension should contain the measurements
% The second dimension could be different samples (e.g. voxels)

if size(x)~=size(y), error('X and Y must have same size for linear fitting'); end
% Use the fact that slope = cov(x,y) / cov(x,x)
lengthX = size(x,1);
numerator = sum(x.*y) - sum(x).*sum(y) / lengthX;
denominator = sum(x.^2) - sum(x).^2 / lengthX;
fittedSlope = numerator ./ denominator;
% Line of best fit has to pass through point (meanX, meanY)
% Use this fact to get intecept
fittedIntercept = mean(y) - fittedSlope .* mean(x);
end % END OF LinLeastSquares


text:
%FIGURE_sliding_distance

%%%% EXPERIMENT %%%%
%
% For a slendar rod toppling from rest, plot how the slip distance varies
% as a function of coefficient of friction.
%
%

clc; clear;

setup.Tspan = [0,10];  %only used for timeout of the simulation
setup.tol = 1e-12;   %Accuracy of the intergation method
setup.dataFreq = 750;   %How much data to return?
setup.solver = @ode45;   %Tell the simulation to use a still solver.

%%%% HACK %%%%
% If the perturbation is too small, then numerical errors start to creep
% into the results. If perturbation is larger than the critical angle of
% the smallest non-zero coefficient of friction, then it will cause a
% direct error by violating the initial phase assumption.
setup.perturbation = 1e-2;
%%%% DONE %%%%

MoI = logspace(-3,0,250);
Mu = [0,0.05,0.1,0.5,1,5,10,inf];
P.m = 1;
P.g = 1;
P.L = 1;

%%%% Run big experiment %%%%
Data = runSlipExperiment(setup,MoI,Mu);

%%%% Plotting %%%%
N_mu = length(Data);
xBnd = [min(Data(1).moi), max(Data(1).moi)];
Mu = zeros(N_mu,1);
style = {'k--','b--','r--','m--','g--',...
    'k-','b-','r-','m-','g-'...
    'k:','b:','r:','m:','g:'};
LINEWIDTH = 3;
FontSize.Title = 16;
FontSize.label = 12;
names = cell(1,N_mu);
for i=1:N_mu
    Mu(i) = Data(i).mu;
    names{i} = num2str(['u = ' num2str(Mu(i))]);
end

%%%% Distance %%%%
H_dist = figure(100); clf; hold on;
set(H_dist,'Name','SlipDist','NumberTitle','off')
IDX = false(2,length(Mu));
for i=1:length(Data)
    sty = style{mod(i-1,length(style))+1};
    subplot(2,1,1); hold on;
    idx = Data(i).pos~=0;  IDX(1,i) = sum(idx)~=0;
    semilogx(Data(i).moi(idx),Data(i).pos(idx),sty,'LineWidth',LINEWIDTH);
    title('Backwards Slip (falling forward)','FontSize',FontSize.Title)
    xlabel('Moment of Inertia','FontSize',FontSize.label)
    ylabel('Slip Distance','FontSize',FontSize.label)
    set(gca,'Xscale','log')
    subplot(2,1,2); hold on;
    idx = Data(i).neg~=0;  IDX(2,i) = sum(idx)~=0;
    semilogx(Data(i).moi(idx),Data(i).neg(idx),sty,'LineWidth',LINEWIDTH);
    title('Forwards Slip (falling forward)','FontSize',FontSize.Title)
    xlabel('Moment of Inertia','FontSize',FontSize.label)
    ylabel('Slip Distance','FontSize',FontSize.label)
    set(gca,'Xscale','log')
end
%subplot(2,1,1); legend(names(IDX(1,:)),'Location','NorthEast');
% extents = [xBnd,-1,0];  axis(extents);
subplot(2,1,2); legend(names(IDX(2,:)),'Location','SouthWest');
% extents = [xBnd,0,0.14];

%%%% Critical Angle %%%%
H_angle = figure(101); clf; hold on;
set(H_angle,'Name','CriticalAngle','NumberTitle','off')
IDX = false(2,length(Mu));
for i=1:length(Data)
    sty = style{mod(i-1,length(style))+1};

    idx1 = Data(i).pos~=0;  IDX(1,i) = sum(idx1)~=0;
    subplot(2,1,1); hold on
    semilogx(Data(i).moi(idx1),-Data(i).angle(idx1)*180/pi,sty,'LineWidth',LINEWIDTH);
    title('Backwards Slip (falling forward)','FontSize',FontSize.Title)
    xlabel('Moment of Inertia','FontSize',FontSize.label)
    ylabel('Critical Angle (deg)','FontSize',FontSize.label)
    set(gca,'Xscale','log')

    idx2 = Data(i).neg~=0; IDX(2,i) = sum(idx2)~=0;
    idx2a = idx2 & idx1;   %In these special cases, it slips backwards and then forwards. Only plot the FIRST critical angle.
    idx2b = idx2 & ~idx1;
    subplot(2,1,2); hold on;
%     semilogx(Data(i).moi(idx2a),-Data(i).angle(idx2a)*180/pi,sty,'LineWidth',LINEWIDTH);
    semilogx(Data(i).moi(idx2b),-Data(i).angle(idx2b)*180/pi,sty,'LineWidth',LINEWIDTH);
    title('Forwards Slip (falling forward)','FontSize',FontSize.Title)
    xlabel('Moment of Inertia','FontSize',FontSize.label)
    ylabel('Critical Angle (deg)','FontSize',FontSize.label)
    set(gca,'Xscale','log')
end

subplot(2,1,1); legend(names(IDX(1,:)),'Location','NorthEast');
% extents = [xBnd,0,45];  axis(extents);
subplot(2,1,2); legend(names(IDX(2,:)),'Location','NorthWest');
% extents = [xBnd,30,60];

%%%% SAVE %%%%
save('DATA_Slip_vs_Mu.mat','Data');
save2pdf('../WriteUp/Figures/Slip_vs_Mu__Distance.pdf',H_dist,600);
save2pdf('../WriteUp/Figures/Slip_vs_Mu__Angle.pdf',H_angle,600);

author: MatthewPeterKelly

repo: dscTutorials


text:
function value = r8_gamma ( x )

%*****************************************************************************80
%
%% R8_GAMMA evaluates Gamma(X) for a real argument.
%
%  Discussion:
%
%    This routine calculates the gamma function for a real argument X.
%
%    Computation is based on an algorithm outlined in reference 1.
%    The program uses rational functions that approximate the gamma
%    function to at least 20 significant decimal digits.  Coefficients
%    for the approximation over the interval (1,2) are unpublished.
%    Those for the approximation for 12 <= X are from reference 2.
%
%    MATLAB provides a GAMMA function, which is likely to be faster, more
%    accurate, and which vectorizes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by William Cody, Laura Stoltz.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    William Cody,
%    An Overview of Software Development for Special Functions,
%    in Numerical Analysis Dundee, 1975,
%    edited by GA Watson,
%    Lecture Notes in Mathematics 506,
%    Springer, 1976.
%
%    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
%    Charles Mesztenyi, John Rice, Henry Thatcher,
%    Christoph Witzgall,
%    Computer Approximations,
%    Wiley, 1968,
%    LC: QA297.C64.
%
%  Parameters:
%
%    Input, real X, the argument of the function.
%
%    Output, real VALUE, the value of the function.
%

%
%  Coefficients for minimax approximation over (12, INF).
%
  c = [ ...
   -1.910444077728E-03, ...
    8.4171387781295E-04, ...
   -5.952379913043012E-04, ...
    7.93650793500350248E-04, ...
   -2.777777777777681622553E-03, ...
    8.333333333333333331554247E-02, ...
    5.7083835261E-03 ];
%
%  Mathematical constants
%
  sqrtpi = 0.9189385332046727417803297;
%
%  Machine dependent parameters
%
  xbig = 171.624E+00;
  xminin = 2.23E-308;
  eps = 2.22E-16;
  xinf = 1.79E+308;
%
%  Numerator and denominator coefficients for rational minimax
%  approximation over (1,2).
%
  p = [ ...
   -1.71618513886549492533811E+00, ...
    2.47656508055759199108314E+01, ...
   -3.79804256470945635097577E+02, ...
    6.29331155312818442661052E+02, ...
    8.66966202790413211295064E+02, ...
   -3.14512729688483675254357E+04, ...
   -3.61444134186911729807069E+04, ...
    6.64561438202405440627855E+04 ];

  q = [ ...
   -3.08402300119738975254353E+01, ...
    3.15350626979604161529144E+02, ...
   -1.01515636749021914166146E+03, ...
   -3.10777167157231109440444E+03, ...
    2.25381184209801510330112E+04, ...
    4.75584627752788110767815E+03, ...
   -1.34659959864969306392456E+05, ...
   -1.15132259675553483497211E+05 ];

  parity = 0;
  fact = 1.0;
  n = 0;
  y = x;
%
%  Argument is negative.
%
  if ( y <= 0.0 )

    y = - x;
    y1 = floor ( y );
    res = y - y1;

    if ( res ~= 0.0 )

      if ( y1 ~= floor ( y1 * 0.5 ) * 2.0 )
        parity = 1;
      end

      fact = - pi / sin ( pi * res );
      y = y + 1.0;

    else

      res = xinf;
      value = res;
      return

    end

  end
%
%  Argument is positive.
%
  if ( y < eps )
%
%  Argument < EPS.
%
    if ( xminin <= y )
      res = 1.0 / y;
    else
      res = xinf;
      value = res;
      return
    end

  elseif ( y < 12.0 )

    y1 = y;
%
%  0.0 < argument < 1.0.
%
    if ( y < 1.0 )

      z = y;
      y = y + 1.0;
%
%  1.0 < argument < 12.0.
%  Reduce argument if necessary.
%
    else

      n = floor ( y ) - 1;
      y = y - n;
      z = y - 1.0;

    end
%
%  Evaluate approximation for 1.0 < argument < 2.0.
%
    xnum = 0.0;
    xden = 1.0;
    for i = 1 : 8
      xnum = ( xnum + p(i) ) * z;
      xden = xden * z + q(i);
    end

    res = xnum / xden + 1.0;
%
%  Adjust result for case  0.0 < argument < 1.0.
%
    if ( y1 < y )

      res = res / y1;
%
%  Adjust result for case 2.0 < argument < 12.0.
%
    elseif ( y < y1 )

      for i = 1 : n
        res = res * y;
        y = y + 1.0;
      end

    end

  else
%
%  Evaluate for 12.0 <= argument.
%
    if ( y <= xbig )

      ysq = y * y;
      sum = c(7);
      for i = 1 : 6
        sum = sum / ysq + c(i);
      end
      sum = sum / y - y + sqrtpi;
      sum = sum + ( y - 0.5 ) * log ( y );
      res = exp ( sum );

    else

      res = xinf;
      value = res;
      return

    end

  end
%
%  Final adjustments and return.
%
  if ( parity )
    res = - res;
  end

  if ( fact ~= 1.0 )
    res = fact / res;
  end

  value = res;

  return
end


text:
function o = boxoverlap(a, b)
% Compute the symmetric intersection over union overlap between a set of
% bounding boxes in a and a single bounding box in b.
%
% a  a matrix where each row specifies a bounding box
% b  a single bounding box

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

x1 = max(a(:,1), b(1));
y1 = max(a(:,2), b(2));
x2 = min(a(:,3), b(3));
y2 = min(a(:,4), b(4));

w = x2-x1+1;
h = y2-y1+1;
inter = w.*h;
aarea = (a(:,3)-a(:,1)+1) .* (a(:,4)-a(:,2)+1);
barea = (b(3)-b(1)+1) * (b(4)-b(2)+1);
% intersection over union overlap
o = inter ./ (aarea+barea-inter);
% set invalid entries to 0 overlap
o(w <= 0) = 0;
o(h <= 0) = 0;

repo: PRW-baseline


text:
function calculateQualityMeasures(inputDir, datasetConfig, resultDir, binningFactors, numberOfFrames, numberOfFrames_val, sliding_val, sr_method_val, binning_val, scenes_val, compressions_val, measure_val, border)

    % Get configuration for this dataset.
    load([inputDir,'/', datasetConfig, '.mat']);

    % Iterate over compression settings.
    for compressIdx = compressions_val

        if isnan(compressIdx)
            compress_folder = 'Uncoded';
        else
            compress_affix = ['QP', num2str(compressIdx)];
            compress_folder = ['H265', compress_affix];
        end

        % Iterate over all scenes.
        for datasetIdx = scenes_val
            % Iterate over the different binning factors.
            for binningFactorIdx = binning_val
                % Iterate over different SR methods.
                for sr_method = sr_method_val

                    if sr_method == 0
                        continue;
                    end

                    % Iterate over different sequence lengths.
                    for numberOfFramesIdx = numberOfFrames_val
                        % Iterate over different sliding windows.
                        for slidingIdx = sliding_val

                            gtFilename = [resultDir, filesep, evalData.scenes{datasetIdx}, filesep, 'Uncoded', filesep, 'mat', filesep, evalData.motionTypes{datasetIdx}, '_bin', num2str(binningFactors(binningFactorIdx)), '_sr0', '_f', num2str(numberOfFrames(numberOfFramesIdx,binningFactorIdx)), '_win', num2str(slidingIdx,'%02d'), '.mat'];
                            srFilename = [resultDir ,filesep, evalData.scenes{datasetIdx}, filesep, compress_folder,filesep, 'mat', filesep, evalData.motionTypes{datasetIdx}, '_bin', num2str(binningFactors(binningFactorIdx)), '_sr', num2str(sr_method), '_f', num2str(numberOfFrames(numberOfFramesIdx,binningFactorIdx)), '_win', num2str(slidingIdx,'%02d'), '.mat'];
                            if ~exist(gtFilename, 'file') || ~exist(srFilename, 'file')
                                fprintf('Ground truth and result not available (%s, %s)\n', gtFilename, srFilename);
                            else
                                % Calculate the different quality measures.
                                groundTruth = [];
                                SR = [];
                                SR_aligned = [];
                                qualityMeasuresVec = qualityMeasures;
                                for measureIdx = measure_val
                                    % Prepare result directory for current quality metric
                                    if ~exist([resultDir, filesep,evalData.scenes{datasetIdx}, filesep,compress_folder, filesep,'quality_qm', num2str(measureIdx)], 'dir')
                                        mkdir([resultDir, filesep,evalData.scenes{datasetIdx}, filesep,compress_folder], ['quality_qm', num2str(measureIdx)]);
                                    end

                                    if exist([resultDir, filesep, evalData.scenes{datasetIdx}, filesep, compress_folder, filesep, 'quality_qm', num2str(measureIdx), filesep,evalData.motionTypes{datasetIdx}, '_bin', num2str(binningFactors(binningFactorIdx)), '_sr',num2str(sr_method), '_f', num2str(numberOfFrames(numberOfFramesIdx,binningFactorIdx)), '_win', num2str(slidingIdx,'%02d'), '_qm', num2str(measureIdx), '.mat'],'file')
                                        disp('Quality metric already calculated!');
                                    else
                                        if isempty(groundTruth)
                                            % Load the ground truth.
                                            load(gtFilename);
                                            groundTruth = srImages.groundTruth;
                                        end

                                        if isempty(SR)
                                            % Load the SR image.
                                            load(srFilename);
                                            srMethodName = fieldnames(srImages);
                                            srMethodName = srMethodName{1};
                                            SR = getfield(srImages, srMethodName);
                                        end

                                        % Align SR image to the ground truth if
                                        % desired (depends on the SR method).
                                        if isempty(SR_aligned)
                                            SR_aligned = alignToGroundTruth(SR, sr_method, binningFactors(binningFactorIdx));
                                        end

                                        % Calculate the quality measure.
                                        qualityMeasure = calculateQualityMeasureForImage(SR_aligned, groundTruth, qualityMeasuresVec{measureIdx}, border);
                                        save([resultDir, filesep, evalData.scenes{datasetIdx}, filesep, compress_folder, filesep, 'quality_qm', num2str(measureIdx), filesep, evalData.motionTypes{datasetIdx}, '_bin', num2str(binningFactors(binningFactorIdx)), '_sr',num2str(sr_method), '_f', num2str(numberOfFrames(numberOfFramesIdx,binningFactorIdx)), '_win', num2str(slidingIdx,'%02d'), '_qm', num2str(measureIdx), '.mat'], 'qualityMeasure');
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

end

function sr = alignToGroundTruth(sr, srMethodIdx, binningFactor)

    srMethods = SRMethods;
    if ~ismatrix(sr)
        sr = sr(:,:,1);
    end
    sr = double(sr);

    if strcmp(srMethods(srMethodIdx).name, 'vsrnet')
        % Scale intensities to [0, 1].
        sr = sr / 255;
    end

    if strcmp(srMethods(srMethodIdx).pixelSampling, 'topLeft')
        % Geometric alignment of SR image to ground truth for algorithms
        % that use the top left corner for pixel sampling.
        [x_mesh, y_mesh] = meshgrid(1:size(sr, 2), 1:size(sr,1));
        x_mesh_shift = x_mesh - (binningFactor - 1) / 2;
        y_mesh_shift = y_mesh - (binningFactor - 1) / 2;
        sr = griddata(x_mesh, y_mesh, sr, x_mesh_shift, y_mesh_shift, 'cubic');
    end

end

function qualityMeasure = calculateQualityMeasureForImage(SR, groundTruth, qualityMeasureName, border)

    % Crop boundary of ground truth and super-resolved image.
    groundTruth = groundTruth( (border+1):(end-border), (border+1):(end-border) );
    SR = SR( (border+1):(end-border), (border+1):(end-border) );

    % Calculate the desired measure.
    qualityMeasure.(qualityMeasureName) = callQualityMeasure(SR, groundTruth, qualityMeasureName);

end

author: thomas-koehler

repo: SupER


text:
function [U,S,output] = lmlra_minf(T,U0,S0,options)
%LMLRA_MINF LMLRA by unconstrained nonlinear optimization.
%   [U,S,output] = lmlra_minf(T,U0,S0) computes the factor matrices U{1},
%   ..., U{N} and core tensor S belonging to a low multilinear rank
%   approximation of the N-th order tensor T by minimizing
%   0.5*frob(T-lmlragen(U,S))^2. Each term U{r} is a cell array of N factor
%   matrices U{r}{n}, followed by a core tensor U{r}{N+1}. The algorithm is
%   initialized with the factor matrices U0{n} and core tensor S0. The
%   structure output returns additional information:
%
%      output.Name  - The name of the selected algorithm.
%      output.<...> - The output of the selected algorithm.
%
%   lmlra_minf(T,U0,S0,options) may be used to set the following options:
%
%      options.Algorithm =     - The desired optimization method.
%      [{@minf_lbfgsdl}|...
%       @minf_lbfgs|@minf_ncg]
%      options.<...>           - Parameters passed to the selected method,
%                                e.g., options.TolFun, options.TolX and
%                                options.LineSearchOptions. See also help
%                                [options.Algorithm].
%
%   See also lmlra_nls.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.
%   [2] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables," SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.

% Forward problem to btd_minf as a one-term BTD.
if nargin < 4, options = struct; end
[U,output] = btd_minf(T,{[U0(:).',S0]},options);
S = U{1}{end};
U = U{1}(1:end-1);


text:
function M = stiefelfactory(n, p, k, gpuflag)
% Returns a manifold structure to optimize over orthonormal matrices.
%
% function M = stiefelfactory(n, p)
% function M = stiefelfactory(n, p, k)
% function M = stiefelfactory(n, p, k, gpuflag)
%
% The Stiefel manifold is the set of orthonormal nxp matrices. If k
% is larger than 1, this is the Cartesian product of the Stiefel manifold
% taken k times. The metric is such that the manifold is a Riemannian
% submanifold of R^nxp equipped with the usual trace inner product, that
% is, it is the usual metric.
%
% Points are represented as matrices X of size n x p x k (or n x p if k=1,
% which is the default) such that each n x p matrix is orthonormal,
% i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
% i = 1 : k if k > 1. Tangent vectors are represented as matrices the same
% size as points.
%
% The default retraction is QR-based: it is only a first-order retraction.
% To use the polar retraction (which is second order), run
%    M.retr = M.retr_polar;
% after creating M with this factory. This can be reverted with
%    M.retr = M.retr_qr;
% If used, you may also want to update M.invretr similarly.
%
% Set gpuflag = true to have points, tangent vectors and ambient vectors
% stored on the GPU. If so, computations can be done on the GPU directly.
%
% By default, k = 1 and gpuflag = false.
%
% See also: grassmannfactory rotationsfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%  July  5, 2013 (NB) : Added ehess2rhess.
%  Jan. 27, 2014 (BM) : Bug in ehess2rhess corrected.
%  June 24, 2014 (NB) : Added true exponential map and changed the randvec
%                       function so that it now returns a globally
%                       normalized vector, not a vector where each
%                       component is normalized (this only matters if k>1).
%  July 17, 2018 (NB) : Now both QR (default) and polar retractions are
%                       directly accessible, and their inverses are also
%                       implemented.
%  Aug.  2, 2018 (NB) : Added GPU support: just set gpuflag = true.
%  June 18, 2019 (NB) : Using qr_unique for retr and rand.
%  July  9, 2019 (NB) : Added a comment about QR retraction being first
%                       order only.
%  Jan.  8, 2021 (NB) : Added tangent2ambient+tangent2ambient_is_identity.

    assert(n >= p, 'The dimension n must be larger than the dimension p.');

    if ~exist('k', 'var') || isempty(k)
        k = 1;
    end
    if ~exist('gpuflag', 'var') || isempty(gpuflag)
        gpuflag = false;
    end

    % If gpuflag is active, new arrays (e.g., via rand, randn, zeros, ones)
    % are created directly on the GPU; otherwise, they are created in the
    % usual way (in double precision).
    if gpuflag
        array_type = 'gpuArray';
    else
        array_type = 'double';
    end

    if k == 1
        M.name = @() sprintf('Stiefel manifold St(%d, %d)', n, p);
    elseif k > 1
        M.name = @() sprintf('Product Stiefel manifold St(%d, %d)^%d', n, p, k);
    else
        error('k must be an integer no less than 1.');
    end

    M.dim = @() k*(n*p - .5*p*(p+1));

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d(:));

    M.dist = @(x, y) error('stiefel.dist not implemented yet.');

    M.typicaldist = @() sqrt(p*k);

    M.proj = @projection;
    function Up = projection(X, U)

        XtU = multiprod(multitransp(X), U);
        symXtU = multisym(XtU);
        Up = U - multiprod(X, symXtU);

% The code above is equivalent to, but faster than, the code below.
%
%     Up = zeros(size(U));
%     function A = sym(A), A = .5*(A+A'); end
%     for i = 1 : k
%         Xi = X(:, :, i);
%         Ui = U(:, :, i);
%         Up(:, :, i) = Ui - Xi*sym(Xi'*Ui);
%     end

    end

    M.tangent = M.proj;

    M.tangent2ambient_is_identity = true;
    M.tangent2ambient = @(X, U) U;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)
        XtG = multiprod(multitransp(X), egrad);
        symXtG = multisym(XtG);
        HsymXtG = multiprod(H, symXtG);
        rhess = projection(X, ehess - HsymXtG);
    end

    M.retr_qr = @retraction_qr;
    function Y = retraction_qr(X, U, t)
        % It is necessary to call qr_unique rather than simply qr to ensure
        % this is a retraction, to avoid spurious column sign flips.
        if nargin < 3
            Y = qr_unique(X + U);
        else
            Y = qr_unique(X + t*U);
        end
    end

    M.invretr_qr = @invretr_qr;
    function U = invretr_qr(X, Y)
        XtY = multiprod(multitransp(X), Y);
        R = zeros(p, p, k, array_type);
        H = 2*eye(p, array_type);
        for kk = 1 : k
            % For each slice, assuming the inverse retraction is well
            % defined for the given inputs, we have:
            %   X + U = YR
            % Left multiply with X' to get
            %   I + X'U = X'Y R
            % Since X'U is skew symmetric for a tangent vector U at X, add
            % up this equation with its transpose to get:
            %   2I = (X'Y) R + R' (X'Y)'
            % Contrary to the polar factorization, here R is not symmetric
            % but it is upper triangular. As a result, this is not a
            % Sylvester equation and we must solve it differently.
            R(:, :, kk) = solve_for_triu(XtY(:, :, kk), H);
            % Then,
            %   U = YR - X
            % which is what we compute below.
        end
        U = multiprod(Y, R) - X;
    end

    M.retr_polar = @retraction_polar;
    function Y = retraction_polar(X, U, t)
        if nargin < 3
            Y = X + U;
        else
            Y = X + t*U;
        end
        for kk = 1 : k
            [u, s, v] = svd(Y(:, :, kk), 'econ'); %#ok
            Y(:, :, kk) = u*v';
        end
    end

    M.invretr_polar = @invretr_polar;
    function U = invretr_polar(X, Y)
        XtY = multiprod(multitransp(X), Y);
        MM = zeros(p, p, k, array_type);
        H = 2*eye(p, array_type);
        for kk = 1 : k
            % For each slice, assuming the inverse retraction is well
            % defined for the given inputs, we have:
            %   X + U = YM
            % Left multiply with X' to get
            %   I + X'U = X'Y M
            % Since X'U is skew symmetric for a tangent vector U at X, add
            % up this equation with its transpose to get:
            %   2I = (X'Y) M + M' (X'Y)'
            %      = (X'Y) M + M (X'Y)'   since M is symmetric.
            % Solve for M symmetric with a call to sylvester:
            MM(:, :, kk) = sylvester_nochecks(XtY(:, :, kk), XtY(:, :, kk)', H);
            % Note that the above is really a Lyapunov equation: it could
            % be solved faster by exploiting the fact the same matrix
            % appears twice on the left, with one the transpose of the
            % other. Then,
            %   U = YM - X
            % which is what we compute below.
        end
        U = multiprod(Y, MM) - X;
    end

    % By default, we use the QR retraction
    M.retr = M.retr_qr;
    M.invretr = M.invretr_qr;

    M.exp = @exponential;
    function Y = exponential(X, U, t)
        if nargin == 2
            tU = U;
        else
            tU = t*U;
        end
        Y = zeros(size(X), array_type);
        I = eye(p, array_type);
        Z = zeros(p, array_type);
        for kk = 1 : k
            % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
            Xkk = X(:, :, kk);
            Ukk = tU(:, :, kk);
            Y(:, :, kk) = [Xkk Ukk] * ...
                         expm([Xkk'*Ukk , -Ukk'*Ukk ; I , Xkk'*Ukk]) * ...
                         [ expm(-Xkk'*Ukk) ; Z ];
        end

    end

    M.hash = @(X) ['z' hashmd5(X(:))];

    M.rand = @() qr_unique(randn(n, p, k, array_type));

    M.randvec = @randomvec;
    function U = randomvec(X)
        U = projection(X, randn(n, p, k, array_type));
        U = U / norm(U(:));
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(n, p, k, array_type);

    M.transp = @(x1, x2, d) projection(x2, d);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, p, k]);
    M.vecmatareisometries = @() true;

    % Automatically convert a number of tools to support GPU.
    if gpuflag
        M = factorygpuhelper(M);
    end

end

author: NicolasBoumal

repo: manopt


text:
function HTotal=spherAngUvCrossHessian(uv,systemType,Ms,Muv)
%%SPHERANGUVCROSSHESSIAN Determine second partial derivative matrices of 3D
%               spherical angular components with respect to u-v direction
%               cosines.
%
%INPUTS: uv A 2XN or 3XN (if the third components of the unit vector) set
%          of direction [u;v;w] cosines values in 3D. If the third
%          component of the unit vector is omitted, it is assumed to be
%          positive.
% systemType An optional parameter specifying the axes from which the
%          angles for the spherical coordinate system are measured in
%          radians. Possible vaues are
%          0 (The default if omitted) Azimuth is measured counterclockwise
%            from the x-axis in the x-y plane. Elevation is measured up
%            from the x-y plane (towards the z-axis). This is consistent
%            with common spherical coordinate systems for specifying
%            longitude (azimuth) and geocentric latitude (elevation).
%          1 Azimuth is measured counterclockwise from the z-axis in the
%            z-x plane. Elevation is measured up from the z-x plane
%            (towards the y-axis). This is consistent with some spherical
%            coordinate systems that use the z axis as the boresight
%            direction of the radar.
%          2 This is the same as 0 except instead of being given
%            elevation, one desires the angle away from the z-axis, which
%            is (pi/2-elevation).
%   Ms,Muv If either the spherical coordinate system or the u-v coordinate
%          system is rotated compared to the global Cartesian coordinate
%          system, these optional 3X3 matrices provide the rotations. Ms
%          is a 3X3 matrix to go from the alignment of a global
%          Cartesian coordinate system to that in which the spherical
%          coordinates are computed. Similarly, Muv is a rotation matrix
%          to go from the alignment of a global Cartesian cordinate system
%          to that in which the u-v(-w) coordinates are computed. If
%          either of these in omitted or an empty matrix is passed, then
%          the missing one is replaced with the identity matrix.
%
%OUTPUTS: HTotal A 2X2X2XN matrix of second derivatives where
%                HTotal(:,:,i,j) is the Hessian matrix of the ith
%                component of [azimuth;elevation] evaluated at the jth
%                point in uv. The ordering of the second derivatives in the
%                i,jth Hessian matrix is [d/du^2, d/(dudv);
%                                         d/(dudv), d/dv^2]
%
%EXAMPLE:
%Here, we verify that the Hessian matrix computed by this function is close
%to the computed using forward differencing of the gradient.
% systemType=0;
% uv=[0.1;-0.2];
%
% H=spherAngUvCrossHessian(uv,systemType);
% J=spherAngUvCrossGrad(uv,systemType);
% epsVal=1e-8;
% uv1=uv+[epsVal;0];
% J1=spherAngUvCrossGrad(uv1,systemType);
% dAz=(J1-J)/epsVal;
%
% uv1=uv+[0;epsVal];
% J1=spherAngUvCrossGrad(uv1,systemType);
% dEl=(J1-J)/epsVal;
%
% HNumDiff=zeros(2,2,2,1);
%
% %Derivatives of azimuth
% HNumDiff(1,1,1)=dAz(1,1);%dAz/du^2
% HNumDiff(1,2,1)=dAz(1,2);%dAz/(dudv)
% HNumDiff(2,1,1)=HNumDiff(1,2,1);
% HNumDiff(2,2,1)=dEl(1,2);%dAz/dv^2
%
% %Derivatives of elevation
% HNumDiff(1,1,2)=dAz(2,1);%dEl/du^2
% HNumDiff(1,2,2)=dAz(2,2);%dEl/(dudv)
% HNumDiff(2,1,2)=HNumDiff(1,2,2);
% HNumDiff(2,2,2)=dEl(2,2);%dEl/dv^2
%
% max(abs(HNumDiff(:)-H(:)))
%One will see that the difference between the true Hessian and the Hessian
%from numeric differentiation is on the order of 8e-7, which indicates good
%agreement.
%
%June 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(Muv))
    Muv=eye(3,3);
end

if(nargin<3||isempty(Ms))
    Ms=eye(3,3);
end

if(nargin<2||isempty(systemType))
   systemType=0;
end

hasW=size(uv,1)>2;
N=size(uv,2);

M=Ms*Muv';
m11=M(1,1);
m12=M(1,2);
m13=M(1,3);
m21=M(2,1);
m22=M(2,2);
m23=M(2,3);
m31=M(3,1);
m32=M(3,2);
m33=M(3,3);

HTotal=zeros(2,2,2,N);
for curPoint=1:N
    if(hasW==false)
        uvwCur=[uv(:,curPoint);sqrt(1-uv(1,curPoint)^2-uv(2,curPoint)^2)];
    else
        uvwCur=uv(:,curPoint);
    end
    u=uvwCur(1);
    v=uvwCur(2);
    w=uvwCur(3);

    uvwCur=M*uvwCur;
    u1=uvwCur(1);
    v1=uvwCur(2);
    w1=uvwCur(3);

    %First derivatives
    du1du=m11-m13*u/w;
    du1dv=m12-m13*v/w;
    dv1du=m21-m23*u/w;
    dv1dv=m22-m23*v/w;
    dw1du=m31-m33*u/w;
    dw1dv=m32-m33*v/w;

    w3=w^3;
    %Second derivatives
    d2u1dudu=m13*(v^2-1)/w3;
    d2u1dudv=-m13*u*v/w3;
    d2u1dvdv=m13*(u^2-1)/w3;
    d2v1dudu=m23*(v^2-1)/w3;
    d2v1dudv=-m23*u*v/w3;
    d2v1dvdv=m23*(u^2-1)/w3;
    d2w1dudu=m33*(v^2-1)/w3;
    d2w1dudv=-m33*u*v/w3;
    d2w1dvdv=m33*(u^2-1)/w3;

    switch(systemType)
        case 0
            denom1=(u1^2+v1^2)^2;
            denom2=sqrt(1-w1^2)^3;

            dazdu2=(-2*u1*(2*du1du*dv1du*u1-du1du^2*v1+dv1du^2*v1)+(2*du1du*dv1du+d2v1dudu*u1-d2u1dudu*v1)*(u1^2+v1^2))/denom1;
            dazdudv=(d2v1dudv*u1^3+v1^2*(du1dv*dv1du+du1du*dv1dv-d2u1dudv*v1)-u1^2*(du1dv*dv1du+du1du*dv1dv+d2u1dudv*v1)+u1*v1*(2*du1du*du1dv-2*dv1du*dv1dv+d2v1dudv*v1))/denom1;
            dazdv2=(-2*u1*(2*du1dv*dv1dv*u1-du1dv^2*v1+dv1dv^2*v1)+(2*du1dv*dv1dv+d2v1dvdv*u1-d2u1dvdv*v1)*(u1^2+v1^2))/denom1;
            deldu2=(d2w1dudu+dw1du^2*w1-d2w1dudu*w1^2)/denom2;
            deldudv=(d2w1dudv+dw1du*dw1dv*w1-d2w1dudv*w1^2)/denom2;
            deldv2=(d2w1dvdv+dw1dv^2*w1-d2w1dvdv*w1^2)/denom2;

%If there were no rotations, it would be:
%             u2v2=u^2+v^2;
%             dazdu2=(2*u*v)/u2v2^2;
%             dazdudv=(-u^2+v^2)/u2v2^2;
%             dazdv2=-((2*u*v)/u2v2^2);
%             deldu2=-((u^4+v^2-v^4)/(w*sqrt(u2v2))^3);
%             deldudv=-((u*v*(-1+2*u^2+2*v^2))/(w*sqrt(u2v2))^3);
%             deldv2=-((u^2-u^4+v^4)/(w*sqrt(u2v2))^3);
        case 1
            denom1=(u1^2+w1^2)^2;
            denom2=sqrt(1-v1^2)^3;

            dazdu2=(2*u1*(2*du1du*dw1du*u1-du1du^2*w1+dw1du^2*w1)+(-2*du1du*dw1du-d2w1dudu*u1+d2u1dudu*w1)*(u1^2+w1^2))/denom1;
            dazdudv=(u1^2*(du1dv*dw1du+du1du*dw1dv-d2w1dudv*u1)+u1*(-2*du1du*du1dv+2*dw1du*dw1dv+d2u1dudv*u1)*w1-(du1dv*dw1du+du1du*dw1dv+d2w1dudv*u1)*w1^2+d2u1dudv*w1^3)/denom1;
            dazdv2=(2*u1*(2*du1dv*dw1dv*u1-du1dv^2*w1+dw1dv^2*w1)+(-2*du1dv*dw1dv-d2w1dvdv*u1+d2u1dvdv*w1)*(u1^2+w1^2))/denom1;
            deldu2=(d2v1dudu+dv1du^2*v1-d2v1dudu*v1^2)/denom2;
            deldudv=(d2v1dudv+dv1du*dv1dv*v1-d2v1dudv*v1^2)/denom2;
            deldv2=(d2v1dvdv+dv1dv^2*v1-d2v1dvdv*v1^2)/denom2;

%If there were no rotations, it would be:
%             dazdu2=u/w^3;
%             dazdudv=v/w^3;
%             dazdv2=-((u*(-1+u^2+(-1+u^2)*v^2+2*v^4))/(w^3*(-1+v^2)^2));
%             deldu2=0;
%             deldudv=0;
%             deldv2=v/(1-v^2)^(3/2);
        case 2
            denom1=(u1^2+v1^2)^2;
            denom2=sqrt(1-w1^2)^3;

            dazdu2=(-2*u1*(2*du1du*dv1du*u1-du1du^2*v1+dv1du^2*v1)+(2*du1du*dv1du+d2v1dudu*u1-d2u1dudu*v1)*(u1^2+v1^2))/denom1;
            dazdudv=(d2v1dudv*u1^3+v1^2*(du1dv*dv1du+du1du*dv1dv-d2u1dudv*v1)-u1^2*(du1dv*dv1du+du1du*dv1dv+d2u1dudv*v1)+u1*v1*(2*du1du*du1dv-2*dv1du*dv1dv+d2v1dudv*v1))/denom1;
            dazdv2=(-2*u1*(2*du1dv*dv1dv*u1-du1dv^2*v1+dv1dv^2*v1)+(2*du1dv*dv1dv+d2v1dvdv*u1-d2u1dvdv*v1)*(u1^2+v1^2))/denom1;
            deldu2=(-dw1du^2*w1+d2w1dudu*(-1+w1^2))/denom2;
            deldudv=(-dw1du*dw1dv*w1+d2w1dudv*(-1+w1^2))/denom2;
            deldv2=(-dw1dv^2*w1+d2w1dvdv*(-1+w1^2))/denom2;

%If there were no rotations, it would be:
%             u2v2=u^2+v^2;
%             dazdu2=(2*u*v)/u2v2^2;
%             dazdudv=(-u^2+v^2)/u2v2^2;
%             dazdv2=-((2*u*v)/u2v2^2);
%             deldu2=(u^4+v^2-v^4)/(w*sqrt(u2v2))^3;
%             deldudv=(u*v*(-1+2*u^2+2*v^2))/(w*sqrt(u2v2))^3;
%             deldv2=(u^2-u^4+v^4)/(w*sqrt(u2v2))^3;
        case 3
            denom1=(u1^2+v1^2)^2;
            denom2=sqrt(1-w1^2)^3;

            dazdu2=(2*u1*(2*du1du*dv1du*u1-du1du^2*v1+dv1du^2*v1)+(-2*du1du*dv1du-d2v1dudu*u1+d2u1dudu*v1)*(u1^2+v1^2))/denom1;
            dazdudv=(u1^2*(du1dv*dv1du+du1du*dv1dv-d2v1dudv*u1)+u1*(-2*du1du*du1dv+2*dv1du*dv1dv+d2u1dudv*u1)*v1-(du1dv*dv1du+du1du*dv1dv+d2v1dudv*u1)*v1^2+d2u1dudv*v1^3)/denom1;
            dazdv2=(2*u1*(2*du1dv*dv1dv*u1-du1dv^2*v1+dv1dv^2*v1)+(-2*du1dv*dv1dv-d2v1dvdv*u1+d2u1dvdv*v1)*(u1^2+v1^2))/denom1;
            deldu2=(d2w1dudu+dw1du^2*w1-d2w1dudu*w1^2)/denom2;
            deldudv=(d2w1dudv+dw1du*dw1dv*w1-d2w1dudv*w1^2)/denom2;
            deldv2=(d2w1dvdv+dw1dv^2*w1-d2w1dvdv*w1^2)/denom2;

%If there were no rotations, it would be:
%             u2v2=u^2+v^2;
%             dazdu2=-(2*u*v)/u2v2^2;
%             dazdudv=(u-v)*(u+v)/u2v2^2;
%             dazdv2=(2*u*v)/u2v2^2;
%             deldu2=-((u^4+v^2-v^4)/(w*sqrt(u2v2))^3);
%             deldudv=-((u*v*(-1+2*u^2+2*v^2))/(w*sqrt(u2v2))^3);
%             deldv2=-((u^2-u^4+v^4)/(w*sqrt(u2v2))^3);
        otherwise
            error('Invalid system type specified.')
    end

    H=zeros(2,2,2);
    %Derivatives of azimuth
    H(1,1,1)=dazdu2;
    H(1,2,1)=dazdudv;
    H(2,1,1)=H(1,2,1);
    H(2,2,1)=dazdv2;
    %Derivatives of elevation
    H(1,1,2)=deldu2;
    H(1,2,2)=deldudv;
    H(2,1,2)=H(1,2,2);
    H(2,2,2)=deldv2;

    HTotal(:,:,:,curPoint)=H;
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [R] = spm_resels(FWHM,L,SPACE)
% Returns the RESEL counts of a search volume
% FORMAT [R] = spm_resels(FWHM,L,SPACE)
% FWHM       - smoothness of the component fields {FWHM - voxels}
% L          - space definition            {in voxels}
%                L = radius                {Sphere}
%                L = [height width length] {Box}
%                L = XYZ pointlist         {Discrete voxels}
%                L = Mapped image volume   {Image}
% SPACE      - Search space
%               'S' - Sphere
%               'B' - Box
%               'V' - Discrete voxels
%               'I' - Image VOI
%
% R          - RESEL counts {adimensional}
%
%__________________________________________________________________________
% For one or two dimensional spaces the appropriate manifold is
% used (e.g. sphere -> disc -> line).
%
% Reference : Worsley KJ et al 1996, Hum Brain Mapp. 4:58-73
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston & Matthew Brett
% $Id: spm_resels.m 3899 2010-05-25 15:36:40Z guillaume $

% Dimensionality
%--------------------------------------------------------------------------
switch SPACE

case 'S'                                                           % Sphere
    %----------------------------------------------------------------------
    s     = L(:)./FWHM(:);
    s     = s(s > 0);
    if length(s) == 2,  SPACE = 'D';  end
    if length(s) == 1,  SPACE = 'L';  end

case 'B'                                                              % Box
    %----------------------------------------------------------------------
    s     = L(:)./FWHM(:);
    s     = s(s > 0);
    if length(s) == 2,  SPACE = 'R';  end
    if length(s) == 1,  SPACE = 'L';  end
end

% Default {sphere - assuming L = volume i.e. number of voxels)
%--------------------------------------------------------------------------
if nargin < 3
    SPACE = 'S';
    L     = (L*(3/4)/pi)^(1/3);
end

% RESEL Counts (R)
%==========================================================================

switch SPACE

case 'S'                                                           % Sphere
    %----------------------------------------------------------------------
    s     = prod(s).^(1/3);
    R     = [1 4*s 2*pi*s^2 (4/3)*pi*s^3];

case 'D'                                                             % Disc
    %----------------------------------------------------------------------
    s     = prod(s).^(1/2);
    R     = [1 pi*s pi*s^2 0];

case 'B'                                                              % Box
    %----------------------------------------------------------------------
    R     = [1 sum(s) (s(1)*s(2) + s(2)*s(3) + s(1)*s(3)) prod(s)];

case 'R'                                                        % Rectangle
    %----------------------------------------------------------------------
    R     = [1 sum(s) prod(s) 0];

case 'L'                                                             % Line
    %----------------------------------------------------------------------
    R     = [1 s 0 0];

case 'V'                                                           % Voxels
    %----------------------------------------------------------------------
    %R     = spm_Pec_resels(L,FWHM);
    V     = zeros(max(L,[],2)');
    LL    = mat2cell(L,ones(1,size(L,1)),size(L,2));
    V(sub2ind(size(V),LL{:})) = 1;
    R     = spm_resels_vol(V,FWHM)';

case 'I'                                                            % Image
    %----------------------------------------------------------------------
    R     = spm_resels_vol(L,FWHM)';

end

author: spm


text:
function x = ndmax(x,dim)
% NDMAX    Multi-dimensional maximization.
% NDMAX(X,DIM) takes the maximum element along the dimensions in DIM,
% and squeezes the result.

% Written by Thomas P Minka
% (c) Microsoft Corporation. All rights reserved.

sz = size(x);
for i=1:length(dim)
  x = max(x, [], dim(i));
end
x = mysqueeze(x);
addflops(2*(prod(sz)-prod(size(x))));


text:
function [Y,x,DCM] = spm_dcm_generate(syn_model,SNR,show_graphics)
% Generate synthetic data from a DCM specification
% FORMAT [Y,x,DCM] = spm_dcm_generate(syn_model,SNR)
%
% syn_model     - Name of synthetic DCM file
% SNR           - Signal to noise ratio [default: 1]
% show_graphics - Whether to plot each timeseries [default: true]
%
% This routine will update the DCM.Y field as follows:
%           Y.y    - synthetic BOLD data
%           Y.secs - overall number of seconds
%           Y.Q    - components of error precision
%
% and will enter neuronal activity (first hidden var in each region) into
% DCM.x
%
% Y           - Simulated (Noisy) BOLD data
% x           - Simulated neuronal activity (first hidden variable in each region)
% DCM         - Full generative model
%
%__________________________________________________________________________
% Copyright (C) 2002-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny & Klaas Enno Stephan
% $Id: spm_dcm_generate.m 6716 2016-02-08 18:21:37Z peter $

% Check parameters and load specified DCM
%--------------------------------------------------------------------------
if isstruct(syn_model)
    DCM       = syn_model;
    syn_model = ['DCM-' date '.mat'];
else
    load(syn_model)
end
if nargin < 2 || isempty(SNR)
    SNR = 1;
end
if nargin < 3 || isempty(show_graphics)
    show_graphics = true;
end

% Unpack
%--------------------------------------------------------------------------
U     = DCM.U;        % inputs
v     = DCM.v;        % number of scans
n     = DCM.n;        % number of regions
m     = size(U.u,2);  % number of inputs

% Check whether the model is stable by examining the eigenvalue
% spectrum for the intrinsic connectivity matrix
%--------------------------------------------------------------------------
[is_stable, eigval] = spm_dcm_check_stability(DCM);
if ~is_stable
    fprintf('Modelled system is potentially unstable:\n');
    fprintf('Lyapunov exponent of combined connectivity matrix is %f\n',max(eigval));
    fprintf('Check the output to ensure that values are in a normal range.\n')
end

% check whether this is a nonlinear DCM
%--------------------------------------------------------------------------
if ~isfield(DCM,'d') || isempty(DCM.d)
    DCM.d = zeros(n,n,0);
    M.IS  = 'spm_int';
else
    M.IS  = 'spm_int_D';
end

% priors
%--------------------------------------------------------------------------
[pE,pC] = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d);

% complete model specification
%--------------------------------------------------------------------------
M.f     = 'spm_fx_fmri';
M.g     = 'spm_gx_state_fmri';
M.x     = sparse(n,5);
M.pE    = pE;
M.pC    = pC;
M.m     = size(U.u,2);
M.n     = size(M.x(:),1);
M.l     = size(M.x,1)*2; % twice as many "outputs" (as hemo and neuronal)
M.N     = 32;
M.dt    = 16/M.N;
M.ns    = v;

% fMRI slice time sampling
%--------------------------------------------------------------------------
try, M.delays = DCM.delays; end
try, M.TE     = DCM.TE;     end

% twice as many "outputs" (hemo and neuronal) - need this to coerce spm_int
%--------------------------------------------------------------------------
M.delays = [M.delays; M.delays];

% Integrate and compute hemodynamic response at v sample points
%--------------------------------------------------------------------------
y      = feval(M.IS,DCM.Ep,M,U);

% Compute required r: standard deviation of additive noise, for all areas
%--------------------------------------------------------------------------
if isinf(SNR)
    r = zeros(n,n);
else
    r = diag(std(y(:,1:n))/SNR);
end

% Add noise
%--------------------------------------------------------------------------
p      = 1;
a      = 1/16;
a      = [1 -a];
K      = inv(spdiags(ones(v,1)*a,-[0:p],v,v));
K      = K*sqrt(v/trace(K*K'));
z      = randn(v,n);
e      = K*z;
Y      = DCM.Y;
Y.Q    = spm_Ce(v*ones(1,n));

Y.y    = y(:,1:n) + e*r;
x      = y(:,n+1:2*n);
Y.secs = Y.dt*v;

% Save synthetic DCM
%--------------------------------------------------------------------------
DCM.Y  = Y;                                    % simulated data
DCM.y  = y(:,1:n);                             % simulated BOLD
DCM.x  = x;                                    % simulated neuronal
DCM.M  = M;                                    % model

save(syn_model, 'DCM', spm_get_defaults('mat.format'));

if spm('CmdLine') || ~show_graphics, return; end

% Display the time series generated
%--------------------------------------------------------------------------
spm_figure('GetWin','Simulated BOLD time series');
t     = Y.dt*[1:1:v];
for i = 1:n,
    subplot(n,1,i);
    plot(t,Y.y(:,i));
    title(sprintf('Region %s', Y.name{i}));
    if i<n set(gca,'XTickLabel',[]); end
end
xlabel('secs');

spm_figure('GetWin','Simulated Neuronal Activity');

t     = Y.dt*[1:1:v];
for i = 1:n,
    subplot(n,1,i);
    plot(t,x(:,i));
    title(sprintf('Region %s', Y.name{i}));
    if i<n set(gca,'XTickLabel',[]); end
end
xlabel('secs');


text:
function [y, gene, times, scale, rawExp] = gpsimLoadEcoliData

% GPSIMLOADECOLIDATA Load in E. coli data for the represion case.
% FORMAT

% RETURN y : the normalised expression levels.
% RETURN yvar : the variance of the normalised expression levels.
% RETURN gene : the gene names and Affymetrix array tags.
% RETURN times : the times of the expression measurements.
% RETURN scale : the scaling factor applied to normalise.
% RETURN rawExp : the raw gene expresion level.
%
% SEEALSO : demEcoliMap1
%
% COPYRIGHT : Pei Gao, Neil D. Lawrence and Magnus Rattary, 2008

% SHEFFIELDML

if exist('./data/ecoliData.mat') == 2
  load('./data/ecoliData.mat');
else
  expData = importdata('./data/ecoliNormalisedData.txt');
  rawExp.data = expData';
  rawExp.genes = {'dinF', 'dinI', 'lexA', 'recA', 'recN', 'ruvA', 'ruvB', ...
                  'sbmC', 'sulA', 'umuC', 'umuD', 'uvrB', 'yebG', 'yjiW' };

  targetInd = [2 3 6 7 9 11 12 14];

% Perform log-normal transformation.
  yFull = exp(rawExp.data(:, targetInd));  % Logs are normally distributed
                                   % ... recover mean in exp space.

  % Rescale so that average standard deviation of curves is 1.
  scale = sqrt(var(yFull));
  scaleMat = ones(size(yFull, 1), 1)*scale;
  y{1} = yFull./scaleMat;
  times = [0 5 10 20 40 60]';
  gene = rawExp.genes(targetInd);
%  save('./data/ecoliData.mat', 'y', 'gene', 'times', 'scale', 'rawExp');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function tf=isNonFailure(obj)
% return whether this test was a non-failure
%
% tf=isSuccess(obj)
%
% Input:
%   obj                     MoxUnitPassedTestOutcome object
%
% Output:
%   tf                      set to true, because a passed test is not a
%                           failure
%
% See also: @MoxUnitTestOutcome

    tf=true;


text:
function varargout = crossentropy(varargin)
% CROSSENTROPY
%
% y = CROSSENTROPY(x,y)
%
% Computes/declares cross entropy -sum(x.*log(y))
%
% See also ENTROPY, KULLBACKLEIBLER

switch class(varargin{1})

    case {'sdpvar','ndsdpvar'}

        varargin{1} = reshape(varargin{1},[],1);
        varargin{2} = reshape(varargin{2},[],1);

        if length(varargin{1})~=length(varargin{2})
            if length(varargin{1})==1
                varargin{1} = repmat(varargin{1},length(varargin{2}),1);
            elseif  length(varargin{2})==1
                varargin{2} = repmat(varargin{2},length(varargin{1}),1);
            else
                error('Dimension mismatch in crossentropy')
            end
        end

        varargout{1} = yalmip('define','crossentropy_internal',[varargin{1};varargin{2}]);

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function L = levin(ntscIm)
% LEVIN Compute an image Laplacian according to the *code* accompanying
% "Colorization Using Optimization" [Levin et al. 2004]. Note that this is not
% the same as what is described in the *text* of the paper. This matrix is not
% symmetric and should be treated as a system matrix rather than a quadratic
% coefficient matrix of an energy.
%
% L = levin(ntscIm)
%
% Inputs:
%   ntscIm  h by w  Luminance image
% Outputs:
%   L  w*h by w*h spares Laplacian
%
% Note: Anat Levin's code, Alec's interpretation

n=size(ntscIm,1); m=size(ntscIm,2);
imgSize=n*m;

%nI(:,:,1)=ntscIm(:,:,1);

indsM=reshape([1:imgSize],n,m);

wd=1;

len=0;
consts_len=0;
col_inds=zeros(imgSize*(2*wd+1)^2,1);
row_inds=zeros(imgSize*(2*wd+1)^2,1);
vals=zeros(imgSize*(2*wd+1)^2,1);
gvals=zeros(1,(2*wd+1)^2);

for j=1:m
   for i=1:n
      consts_len=consts_len+1;

      %if (~colorIm(i,j))
        tlen=0;
        for ii=max(1,i-wd):min(i+wd,n)
           for jj=max(1,j-wd):min(j+wd,m)

              if (ii~=i)|(jj~=j)
                 len=len+1; tlen=tlen+1;
                 row_inds(len)= consts_len;
                 col_inds(len)=indsM(ii,jj);
                 gvals(tlen)=ntscIm(ii,jj,1);
              end
           end
        end
        t_val=ntscIm(i,j,1);
        gvals(tlen+1)=t_val;
        c_var=mean((gvals(1:tlen+1)-mean(gvals(1:tlen+1))).^2);
        csig=c_var*0.6;
        mgv=min((gvals(1:tlen)-t_val).^2);
        if (csig<(-mgv/log(0.01)))
	   csig=-mgv/log(0.01);
	end
	if (csig<0.000002)
	   csig=0.000002;
        end

        gvals(1:tlen)=exp(-(gvals(1:tlen)-t_val).^2/csig);
        gvals(1:tlen)=gvals(1:tlen)/sum(gvals(1:tlen));
        vals(len-tlen+1:len)=-gvals(1:tlen);
      %end

      len=len+1;
      row_inds(len)= consts_len;
      col_inds(len)=indsM(i,j);
      vals(len)=1;

   end
end

vals=vals(1:len);
col_inds=col_inds(1:len);
row_inds=row_inds(1:len);

L=sparse(row_inds,col_inds,vals,consts_len,imgSize);

author: alecjacobson

repo: gptoolbox


text:
function U = initializeSubspace(inputFolder, imgDir, param)
    % initial lize the  subspace by svd

    trainSize = min(length(imgDir), param.trainSize);
    sampleSize = min(param.sampleSize, trainSize);

    randIndex =   randperm(trainSize);

    imageIndex =  randIndex(1:sampleSize);

    %imageIndex =  randperm(trainSize, sampleSize);

    U0=[];
    for i=1:sampleSize
        index = imageIndex(i);
        imgPath = [inputFolder, imgDir(index).name];
        im = imread(imgPath);
        im = im(:);
        U0 = [U0 im];
    end

    [U, D] = svds(double(U0), param.rank);
end


text:
function x = randdirichlet(a)
% RANDDIRICHLET   Sample from Dirichlet distribution
%
% X = RANDDIRICHLET(A) returns a matrix, the same size as A
%       where X(:,j) is sampled from a Dirichlet(A(:,j)) distribution.
% Note: This means the *columns* of X sum to one

x = randgamma(a);
x = bsxfun( @rdivide, x, sum(x,1) );
%Z = sum(x,1);
%x = x./Z(ones(size(a,1),1),:);

author: michaelchughes

repo: NPBayesHMM


text:
function fk=dirft2d3(nj,xj,yj,cj,iflag,nk,sk,tk)
%DIRFT2D3: Direct (slow) computation of nonuniform FFT in R^2 - Type 3.
%
%  FK = DIRFT2D3(NJ,XJ,YJ,CJ,IFLAG,NK,SK,TK);
%
%                 1  nj
%     fk(k)    = -- SUM cj(j) exp(+/-i s(k) xj(j)) exp(+/-i t(k) yj(j))
%                nj j=1
%
%     If (iflag .ge.0) the + sign is used in the exponential.
%     If (iflag .lt.0) the - sign is used in the exponential.
%
%  Input parameters:
%
%     nj     number of sources   (integer)
%     xj,yj  location of sources (real *8)
%
%            on interval [-pi,pi].
%
%     cj     strengths of sources (complex *16)
%     iflag  determines sign of FFT (see above)
%     nk     number of (noninteger) Fourier modes computed
%     sk,tk  k-values (locations) of desired Fourier modes
%
%  Output parameters:
%
%     fk     Fourier transform values (complex *16)
%

fk=zeros(nk,1)+1i*zeros(nk,1);

mex_id_ = 'dirft2d3(i int[x], i double[], i double[], i dcomplex[], i int[x], i int[x], i double[], i double[], io dcomplex[])';
[fk] = nufft2d(mex_id_, nj, xj, yj, cj, iflag, nk, sk, tk, fk, 1, 1, 1);

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function labels=cookLabels(labelsRaw)
% Make sure the labels only contains the fields we want from the raw labels.

% ALTODO: remove me or put in Labels

nLabels=length(labelsRaw);
labels = ...
  struct('t0s',cell(1,nLabels), ...
         't1s',cell(1,nLabels), ...
         'names',cell(1,nLabels), ...
         'flies',cell(1,nLabels), ...
         'off',cell(1,nLabels), ...
         'timestamp',cell(1,nLabels), ...
         'imp_t0s',cell(1,nLabels), ...
         'imp_t1s',cell(1,nLabels));
fieldNames=fieldnames(labels);
for i=1:nLabels
  for j=1:length(fieldNames)
    fieldName=fieldNames{j};
    labels(i).(fieldName)=labelsRaw(i).(fieldName);
  end
end

end

author: kristinbranson

repo: JAABA


text:
function d = distance(u, v)
% Euclidean distance between u and v
%
% :Usage:
% ::
%
%     d = distance(u, v)
%
% :Inputs:
%
%   **u** and **v:**
%        should be column vectors or matrices in k-dim space, where k is
%        the number of columns of u and v.
%
%        if u is a single row, replicates u to size(v,1)

if size(u,1) < size(v,1)
    u = repmat(u,size(v,1),1);
    if size(u,1) < size(v,1), error('Inputs must have same number of rows, or one row for 1st input!'),end
end

delt = u - v;
d = (sum(delt .^ 2,2)) .^.5;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef testPrintingDescriptor < handle

    properties (Access = protected, Abstract)
        fileOutputName
        iter
    end

    methods (Access = protected)

        function hasPassed = hasPassed(obj)
            hasMshChanged = obj.compareFile('.msh');
            hasResChanged = obj.compareFile('.res');
            filesHaveChanged = hasMshChanged || hasResChanged;
            hasPassed = ~filesHaveChanged;
        end

        function selectComputedVar(obj)
        end

    end

    methods (Access = private)

       function sF = obtainSavedPrintedFile(obj,ext)
           out   = obj.fileOutputName;
           fullOutSavedName = [out,'.flavia',ext];
           sF = fullfile('tests','PrintingTests','PrintedFiles',fullOutSavedName);
        end

        function oF = obtainOutPutFile(obj,ext)
           out   = obj.fileOutputName;
           fullOutName = [out,num2str(obj.iter),'.flavia',ext];
           oF = fullfile('Output',out,fullOutName);
        end

        function hasChanged = compareFile(obj,ext)
            sF = obj.obtainSavedPrintedFile(ext);
            oF = obj.obtainOutPutFile(ext);
            hasChanged = FileComperator().areFilesDifferent(sF,oF);
        end

    end

end


text:
function [IDVarNorm,IDVar,IDVarMax]=identVarGauss(w,mu,Sigma,xDim)
%%IDENTVARGAUSS Given a Gaussian mixture representing uncertain target
%               states, determine the normalized and regular identity (ID)
%               variances, as defined in [1]. The identity variance is a
%               measure of uncertainty in the identities of which target is
%               which given the fact that the exact positions of the
%               targets are uncertain (Gaussian). The identity
%               variance provides a metric of how confident one is in
%               which target is which. This function provides an explicit
%               solution under the assumption that the joint probability
%               distribution function (PDF) of the targets is a Gaussian
%               mixture.
%
%INPUTS: w A numHypX1 or 1XnumHyp vector of weights of the hypotheses. Note
%          that w(i)>=0 for all i and sum(w)=1.
%       mu This is an (xDim*numTar)XnumHyp set of stacked state vectors for
%          all of the targets for each hypothesis. Alternatively, if the
%          target states are all uncorrelated, this can be an
%          xDimXnumTarXnumHyp hypermatrix.
%    Sigma This is an (xDim*numTar)XxDim*numTar)XnumHyp set of covariance
%          matrices for all of the stacked mean values in mu for each
%          hypothesis. If all of the target states are uncorrelated, an
%          xDimXxDimXnumTarXnumHyp set of hypermatrices for each target and
%          hypothesis individually can be passed.
%     xDim If the targets are correlated (mu is (xDim*numTar)XnumHyp in
%          size), then the state dimensions size xDim must be explicitly
%          provided.
%
%OUTPUTS: IDVarNorm The normalized ID variance. This is a value between 0
%               and 1. Zero means that the identities of the targets are
%               completely unknown, and one means that they are completely
%               certain.
%         IDVar The non-normalized ID variance.
%      IDVarMax The normalizing constant for the ID variance.
%
%This implements the algorithm of [1].
%
%EXAMPLE:
%Two targets, two hypotheses, three dimensional states.
% mu=zeros(3,2,2);
% mu(:,1,1)=[20;-30;0;];
% mu(:,1,2)=[-15;20;1];
% mu(:,2,1)=[-15;20;3];
% mu(:,2,2)=[-15;20;3];
% Sigma(:,:,1,1)=4*eye(3);
% Sigma(:,:,1,2)=Sigma(:,:,1,1);
% Sigma(:,:,2,1)=[1,  0.5, -0.5;
%                 0.5,  2,  0.5;
%                -0.5,0.5,  3];
% Sigma(:,:,2,2)=2*Sigma(:,:,2,1);
% w=[0.5;0.5];
% IDVarNorm0=identVarGauss(w,mu,Sigma)
% %One will get an ID variance of about 0.8161. However, if all first target
% %hypothese are moved far away from the second target, then the ambiguity is
% %reduced.
% mu(:,1,1)=mu(:,1,1)+500;
% mu(:,1,2)=mu(:,1,2)+500;
% IDVarNorm1=identVarGauss(w,mu,Sigma)
% %Here, one gets an ID variance of essentially 1, meaning that the
% %identities are very clear. on the other hand, if the targets are made to
% %coincide, then the identity variance becomes zero.
% mu(:,2,1)=mu(:,1,1);
% mu(:,2,2)=mu(:,1,2);
% Sigma(:,:,2,1)=Sigma(:,:,1,1);
% Sigma(:,:,2,2)=Sigma(:,:,1,2);
% IDVarNorm2=identVarGauss(w,mu,Sigma)
% %Now, the identity variance is zero.
%
%REFERENCES:
%[1] D. F. Crouse and P. Willett, "Identity variance for multi-object
%    estimation," in Proceedings of SPIE: Signal and Data Processing of
%    Small Targets, vol. 8137, San Diego, CA, 21 Aug. 2011.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numHyp=length(w);

if(nargin<4||isempty(xDim))
    xDim=size(mu,1);
end

if(ndims(mu)==3)
    numTar=size(mu,2);
    totalDim=xDim*numTar;
    mu=reshape(mu,totalDim,numHyp);

    SigmaNew=zeros(xDim*numTar,xDim*numTar,numHyp);
    for curHyp=1:numHyp
        span=1:xDim;
        for curTar=1:numTar
            SigmaNew(span,span,curHyp)=Sigma(:,:,curTar,curHyp);
            span=span+xDim;
        end
    end
    Sigma=SigmaNew;
else
    totalDim=size(mu,1);
    numTar=totalDim/xDim;
end

%The sizes of the adjusted inputs:
%w is numHypX1 or 1XnumHyp
%mu is xDimXnumTarXnumHyp
%Sigma is (xDim*numTar)X(xDim*numTar)XnumHyp

w=w(:);

%The total number of permutations.
numTarPerm=factorial(numTar);

%Here, we implement Equation 22. Note that the matrix H in Equation 31c is
%the same if i and j are swapped.
val1=0;%For the value of the first sum in Equation 22.
val2=0;%For the value of the second sum in 22 and the value of Equation 23.
for curI=1:numTarPerm
    curTerm=w'*calcHTilde(w,mu,Sigma,numTar,xDim,curI,curI)*w;
    val1=val1+curTerm;
    val2=val2+curTerm;

    for curJ=curI+1:numTarPerm
        curTerm=w'*calcHTilde(w,mu,Sigma,numTar,xDim,curI,curJ)*w;
        val2=val2+2*curTerm;%The 2 is for the ordering i,j as well as j,i
    end
end

IDVar=val1/numTarPerm-val2/numTarPerm^2;
IDVarMax=val2*(numTarPerm-1)/(numTarPerm^2);
IDVarNorm=IDVar/IDVarMax;
end

function Ht=calcHTilde(w,mu,Sigma,numTar,xDim,i,j)
%This function implements Equation 31c.

    numHyp=length(w);
    Ht=zeros(numHyp,numHyp);

    idxI=getPermIndices(i-1,numTar,xDim);
    idxJ=getPermIndices(j-1,numTar,xDim);

    for m=1:numHyp
        mumi=mu(idxI,m);
        Sigmami=Sigma(idxI,idxI,m);

        for n=m:numHyp
            munj=mu(idxJ,n);
            Sigmanj=Sigma(idxJ,idxJ,n);

            SigmamiInv=inv(Sigmami);
            SigmanjInv=inv(Sigmanj);

            Sigmamnij=inv(SigmamiInv+SigmanjInv);
            mumnij=Sigmamnij*(SigmamiInv*mumi+SigmanjInv*munj);

            val=0;
            for k=1:numHyp
                val=val+w(k)*GaussianD.PDF(mumnij,mu(:,k),Sigmamnij+Sigma(:,:,k));
            end
            val=val*w(m)*w(n)*GaussianD.PDF(mumi,munj,Sigmami+Sigmanj);

            Ht(m,n)=val;
            Ht(n,m)=val;
        end
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   KUKA KR30 L16 2
%
%   Author: Arturo Gil. Universidad Miguel Hern�ndez de Elche.
%   email: arturo.gil@umh.es date:   08/01/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

robot.DH.theta= '[-q(1) q(2)-pi/2 q(3)   q(4) q(5) q(6)]';
robot.DH.d='[0.815          0     0    1.545 0 0.158]';
robot.DH.a='[0.350         1.2    0.145    0 0 0]';
robot.DH.alpha= '[-pi/2 0 -pi/2 pi/2 -pi/2 0]';

robot.J=[];
robot.name= 'KUKA_KR30_L16_2';

robot.inversekinematic_fn = 'inversekinematic_kuka_kr_30_L16_2(robot, T)';

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-185) deg2rad(185); %Axis 1, minimum, maximum
                deg2rad(-135) deg2rad(35); %Axis 2, minimum, maximum
                deg2rad(-120) deg2rad(158); %Axis 3
                deg2rad(-350) deg2rad(350); %Axis 4
                deg2rad(-130) deg2rad(130); %Axis 5
                deg2rad(-350) deg2rad(350)]; %Axis 6

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(100); %Axis 1, rad/s
                deg2rad(80); %Axis 2, rad/s
                deg2rad(80); %Axis 3, rad/s
                deg2rad(230); %Axis 4, rad/s
                deg2rad(165); %Axis 5, rad/s
                deg2rad(249)];%Axis 6, rad/s
% end effectors maximum velocity
robot.linear_velmax = 2.0; %m/s
robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [200 20 40]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-2 2 -2 2 0 2.5];
%read graphics files
robot = read_graphics(robot);

%DYNAMICS
robot.has_dynamics=0;

repo: ARTE


text:
function [out] = evap_18(p1,p2,p3,S,Ep)
%evap_18

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Exponentially declining evaporation from deficit store
% Constraints:  -
% @(Inputs):    p1   - linear scaling parameter [-]
%               p2   - linear scaling parameter [-]
%               p3   - storage scaling parameter [mm]
%               S    - current storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]

out = p1.*exp(-1.*p2.*S./p3).*Ep;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function test0076 ( )

%*****************************************************************************80
%
%% TEST0076 tests JED_TO_RD and RD_TO_JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0076\n' );
  fprintf ( 1, '  For the RD:\n' );
  fprintf ( 1, '  JED_TO_RD: JED -> RD.\n' );
  fprintf ( 1, '  RD_TO_JED: RD -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    RD                JED (out)\n' );
  fprintf ( 1, '\n' );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    rd2 = jed_to_rd ( jed1 );

    jed3 = rd_to_jed ( rd2 );

    fprintf ( 1, '  %11.2f     %11.2f     %11.2f\n', jed1, rd2, jed3 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

function y = lininterp(t1, x, t2)

if (t2(1) < t1(1))
    t1 = [t2(1); t1];
    x = [x(1,:); x];
end
if (t2(end) > t1(end))
    t1 = [t1; t2(end)];
    x = [x; x(end,:)];
end

y = interp1(t1, x, t2);

author: alex-golts

repo: Video-Stabilization


text:
% function generating toy dataset of 3 manifolds in the 2D space
% X = banana3(n);
% n: number of points per manifold
% X: [3N x 2] matrix of data points
%
% Authors: A. Iscen, G. Tolias, Y. Avrithis, T. Furon, O. Chum. 2017.
% Use rng(155); n = 100; to identically fit our Figure 1 in the CVPR 2017 paper.
function X = banana3(n)

	if nargin < 1, n = 100; end

	a = .25;
	s = 1:n;
	m = n + 1;
	R = [1 0; 0 -1];
	A = banana(m) + repmat([1-a, 0], [m 1]);
	B = banana(m) + repmat([3+a, 0], [m 1]);
	C = banana(n) * R + repmat([2, a], [n 1]);
	X = [A(s+1,:); B(s,:); C];
end

function X = banana(n, a, b)

	if nargin < 2, a = .12; end
	if nargin < 3, b = .1; end

	t = linspace(0, pi, n)';
	u = (1 - 2 * b) * t + b * pi;
	r = a * randn(n, 1) .* sin(u) + 1;
	X = [r .* cos(t) r .* sin(t)];

end

author: ahmetius

repo: diffusion-retrieval


text:
function [gKern, gVarmeans, gVarcovars] = linard2VardistPsi0Gradient(linard2Kern, vardist, covGrad)

% LINARD2VARDISTPSI0GRADIENT description.

% VARGPLVM

A = linard2Kern.inputScales;
gKern = covGrad*sum((vardist.means.*vardist.means) + vardist.covars,1);

gVarmeans = 2*(vardist.means*sparse(diag(A)));
%gVarmeans1 = 2*(repmat(A,size(vardist.means,1),1).*vardist.means);

gVarcovars = ones(size(vardist.means,1),1)*A;

gVarmeans = covGrad*gVarmeans(:)';
gVarcovars = covGrad*gVarcovars(:)';

repo: vargplvm


text:
%% Machine Learning Online Class
%  Exercise 8 | Anomaly Detection and Collaborative Filtering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     estimateGaussian.m
%     selectThreshold.m
%     cofiCostFunc.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% =============== Part 1: Loading movie ratings dataset ================
%  You will start by loading the movie ratings dataset to understand the
%  structure of the data.
%
fprintf('Loading movie ratings dataset.\n\n');

%  Load data
load ('ex8_movies.mat');

%  Y is a 1682x943 matrix, containing ratings (1-5) of 1682 movies on
%  943 users
%
%  R is a 1682x943 matrix, where R(i,j) = 1 if and only if user j gave a
%  rating to movie i

%  From the matrix, we can compute statistics like average rating.
fprintf('Average rating for movie 1 (Toy Story): %f / 5\n\n', ...
        mean(Y(1, R(1, :))));

%  We can "visualize" the ratings matrix by plotting it with imagesc
imagesc(Y);
ylabel('Movies');
xlabel('Users');

%fprintf('\nProgram paused. Press enter to continue.\n');
%pause;

%% ============ Part 2: Collaborative Filtering Cost Function ===========
%  You will now implement the cost function for collaborative filtering.
%  To help you debug your cost function, we have included set of weights
%  that we trained on that. Specifically, you should complete the code in
%  cofiCostFunc.m to return J.

%  Load pre-trained weights (X, Theta, num_users, num_movies, num_features)
load ('ex8_movieParams.mat');

%  Reduce the data set size so that this runs faster
num_users = 4; num_movies = 5; num_features = 3;
X = X(1:num_movies, 1:num_features);
Theta = Theta(1:num_users, 1:num_features);
Y = Y(1:num_movies, 1:num_users);
R = R(1:num_movies, 1:num_users);

%  Evaluate cost function
J = cofiCostFunc([X(:) ; Theta(:)], Y, R, num_users, num_movies, ...
               num_features, 0);

fprintf(['Cost at loaded parameters: %f '...
         '\n(this value should be about 22.22)\n'], J);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============== Part 3: Collaborative Filtering Gradient ==============
%  Once your cost function matches up with ours, you should now implement
%  the collaborative filtering gradient function. Specifically, you should
%  complete the code in cofiCostFunc.m to return the grad argument.
%
fprintf('\nChecking Gradients (without regularization) ... \n');

%  Check gradients by running checkNNGradients
checkCostFunction;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ========= Part 4: Collaborative Filtering Cost Regularization ========
%  Now, you should implement regularization for the cost function for
%  collaborative filtering. You can implement it by adding the cost of
%  regularization to the original cost computation.
%

%  Evaluate cost function
J = cofiCostFunc([X(:) ; Theta(:)], Y, R, num_users, num_movies, ...
               num_features, 1.5);

fprintf(['Cost at loaded parameters (lambda = 1.5): %f '...
         '\n(this value should be about 31.34)\n'], J);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ======= Part 5: Collaborative Filtering Gradient Regularization ======
%  Once your cost matches up with ours, you should proceed to implement
%  regularization for the gradient.
%

%
fprintf('\nChecking Gradients (with regularization) ... \n');

%  Check gradients by running checkNNGradients
checkCostFunction(1.5);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============== Part 6: Entering ratings for a new user ===============
%  Before we will train the collaborative filtering model, we will first
%  add ratings that correspond to a new user that we just observed. This
%  part of the code will also allow you to put in your own ratings for the
%  movies in our dataset!
%
movieList = loadMovieList();

%  Initialize my ratings
my_ratings = zeros(1682, 1);

% Check the file movie_idx.txt for id of each movie in our dataset
% For example, Toy Story (1995) has ID 1, so to rate it "4", you can set
my_ratings(1) = 4;

% Or suppose did not enjoy Silence of the Lambs (1991), you can set
my_ratings(98) = 2;

% We have selected a few movies we liked / did not like and the ratings we
% gave are as follows:
my_ratings(7) = 3;
my_ratings(12)= 5;
my_ratings(54) = 4;
my_ratings(64)= 5;
my_ratings(66)= 3;
my_ratings(69) = 5;
my_ratings(183) = 4;
my_ratings(226) = 5;
my_ratings(355)= 5;

fprintf('\n\nNew user ratings:\n');
for i = 1:length(my_ratings)
    if my_ratings(i) > 0
        fprintf('Rated %d for %s\n', my_ratings(i), ...
                 movieList{i});
    end
end

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ================== Part 7: Learning Movie Ratings ====================
%  Now, you will train the collaborative filtering model on a movie rating
%  dataset of 1682 movies and 943 users
%

fprintf('\nTraining collaborative filtering...\n');

%  Load data
load('ex8_movies.mat');

%  Y is a 1682x943 matrix, containing ratings (1-5) of 1682 movies by
%  943 users
%
%  R is a 1682x943 matrix, where R(i,j) = 1 if and only if user j gave a
%  rating to movie i

%  Add our own ratings to the data matrix
Y = [my_ratings Y];
R = [(my_ratings ~= 0) R];

%  Normalize Ratings
[Ynorm, Ymean] = normalizeRatings(Y, R);

%  Useful Values
num_users = size(Y, 2);
num_movies = size(Y, 1);
num_features = 10;

% Set Initial Parameters (Theta, X)
X = randn(num_movies, num_features);
Theta = randn(num_users, num_features);

initial_parameters = [X(:); Theta(:)];

% Set options for fmincg
options = optimset('GradObj', 'on', 'MaxIter', 100);

% Set Regularization
lambda = 10;
theta = fmincg (@(t)(cofiCostFunc(t, Y, R, num_users, num_movies, ...
                                num_features, lambda)), ...
                initial_parameters, options);

% Unfold the returned theta back into U and W
X = reshape(theta(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(theta(num_movies*num_features+1:end), ...
                num_users, num_features);

fprintf('Recommender system learning completed.\n');

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ================== Part 8: Recommendation for you ====================
%  After training the model, you can now make recommendations by computing
%  the predictions matrix.
%

p = X * Theta';
my_predictions = p(:,1) + Ymean;

movieList = loadMovieList();

[r, ix] = sort(my_predictions, 'descend');
fprintf('\nTop recommendations for you:\n');
for i=1:10
    j = ix(i);
    fprintf('Predicting rating %.1f for movie %s\n', my_predictions(j), ...
            movieList{j});
end

fprintf('\n\nOriginal ratings provided:\n');
for i = 1:length(my_ratings)
    if my_ratings(i) > 0
        fprintf('Rated %d for %s\n', my_ratings(i), ...
                 movieList{i});
    end
end

author: worldveil

repo: coursera-ml


text:
classdef LevelSetSmoothRectangleInclusion < ...
        LevelSetAbstractRectangleInclusion & ...
        PnormDescriptor

    methods (Access = public)

        function obj = LevelSetSmoothRectangleInclusion(cParams)
            obj.m1 = cParams.widthH;
            obj.m2 = cParams.widthV;
            obj.pnorm  = cParams.pnorm;
            obj.compute(cParams);
        end

    end

end


text:
function root_num = p12_root_num ( )

%*****************************************************************************80
%
%% P12_ROOT_NUM returns the number of known roots for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer ROOT_NUM, the number of known roots.
%
  root_num = 1;

  return
end


text:
function value = legendre_3d_monomial_integral ( a, b, p )

%*****************************************************************************80
%
%% LEGENDRE_3D_MONOMIAL_INTEGRAL the Legendre integral of a monomial.
%
%  Discussion:
%
%    The Legendre integral to be evaluated has the form
%
%      I(f) = integral ( z1 <= z <= z2 )
%             integral ( y1 <= y <= y2 )
%             integral ( x1 <= x <= x2 ) x^i y^j z^k dx dy dz
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(3), the lower limits of integration.
%
%    Input, real B(3), the upper limits of integration.
%
%    Input, integer P(3), the exponents of X and Y.
%
%    Output, real VALUE, the value of the exact integral.
%
  value = ( b(1) ^ ( p(1) + 1 ) - a(1) ^ ( p(1) + 1 ) ) / ( p(1) + 1 ) ...
        * ( b(2) ^ ( p(2) + 1 ) - a(2) ^ ( p(2) + 1 ) ) / ( p(2) + 1 ) ...
        * ( b(3) ^ ( p(3) + 1 ) - a(3) ^ ( p(3) + 1 ) ) / ( p(3) + 1 );

  return
end


text:
%-fanDTasia ToolBox------------------------------------------------------------------
% This Matlab script is part of the fanDTasia ToolBox: a Matlab library for Diffusion
% Weighted MRI (DW-MRI) Processing, Diffusion Tensor (DTI) Estimation, High-order
% Diffusion Tensor Analysis, Tensor ODF estimation, Visualization and more.
%
% A Matlab Tutorial on DW-MRI can be found in:
% http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/tutorial.php
%
%----------------------------------------------
%openFDT
%This function opens an image volume in fanDTasia(FDT) file format.
%fanDTasia is an on-line tool for DTI DW-MRI processing.
%
%example 1: A=openFDT('DTI.fdt')
%A is a 4-dimensional matrix that contains a set of volume images.
%
%example 2: [A,B]=openFDT('DTI.fdt')
%A is a 4-dimensional matrix that contains a set of volume images.
%B is a Nx4 matrix that contains the gradient orientations and bvalues.
%B is read from a .txt file with the same name as the fdt file. (In example2: DTI.txt)
%
%----------------------------------------------
%
%Downloaded from Angelos Barmpoutis' web-page.
%
%This program is freely distributable without
%licensing fees and is provided without guarantee
%or warrantee expressed or implied.
%
%Copyright (c) 2007 Angelos Barmpoutis.
%
%VERSION : 20110611
%
%-----------------------------------------------

function [data,bval]=openFDT(name)

fid=fopen(name,'r','b');

for i=1:4
  sz(i)=fread(fid,1,'int');
end

	data=zeros(sz(1),sz(2),sz(3),sz(4));
	for c=1:sz(4)
       for z=1:sz(3)
          for y=1:sz(2)
              data(:,y,z,c)=fread(fid,sz(1),'float');
          end
       end
	end

fclose(fid);

if nargout==2
    bname=name;
    bname(length(name))='t';
    bname(length(name)-1)='x';
    bname(length(name)-2)='t';
    bval=open_TXT(bname);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bval=open_TXT(fnm)
%this simple function opens a txt file that contain 4 columns of numbers
%An example of such a text file is here:
%0.1 100.23 24.3 0.334
%10.345 0.2 1.23 123.4
%... as many lines you want.

fid=fopen(fnm,'r');

if fid==-1
    fprintf(1,'ERROR: FILE CANNOT BE OPENED\n');
    return;
end

bval=[];
b=zeros(1,4);
while feof(fid)==0
    b=fscanf(fid,'%f',4)';
    bval=[bval;b];
end

fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
function this = subsasgn(this, selector, value)
%SUBSASGN   Modify parameters of BA optical flow algorithm
%   THIS = SUBSASGN(THIS, SEL, VAL) sets parameter of a HS optical
%   flow algorithm THIS using given selector SEL to value VAL.  The
%   modified object is returned.
%   Available selectors:
%     - ".optical_flow":   Flow base object
%     - ".pyramid_levels": Levels of the image pyramid
%     - ".warping_mode":   'forward' for forward warping and 'backward'
%                          for backward warping
%     - ".solver":         Selector for the linear equation solver
%                          ('backslash' for the builtin solver and 'sor'
%                          for a MEX-based SOR solver)
%     - ".max_iters":      Maximum number of iterations at each pyramid
%                          level
%     - ".limit_update":   Flag whether to restrict the incremental flow
%                          at each level to [-1, 1]
%
%   This is a member function of the class 'ba_optical_flow'.

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%   $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Allow full access to all fields of the structure
  this = builtin('subsasgn', this, selector, value);