text:
function output = apply_net_with_mask(input, mask)
    global config mem;
    fake_output_for_test = config.NEW_MEM(zeros(config.output_size(1), config.output_size(2), config.output_size(3), config.batch_size));
    input_size = size(input);
    output = config.NEW_MEM((zeros(size(input, 1), size(input, 2), config.chs)));
    p_size = config.MEM.p_size;
    size_differ = config.size_differ;
    count = zeros(size(input, 1), size(input, 2));
    for v = 1 : length(config.MEM.start_rows)
        for h = 1 : length(config.MEM.start_cols)
            v_start = config.MEM.start_rows(v);
            v_end = v_start + p_size - 1;
            h_start = config.MEM.start_cols(h);
            h_end = h_start + p_size - 1;
            if(v_end > input_size(1))
                v_end = input_size(1);
                v_start = v_end - p_size + 1;
            end
            if(h_end > input_size(2))
                h_end = input_size(2);
                h_start = h_end - p_size + 1;
            end
            input_piece = input(v_start:v_end, h_start:h_end,:);
            m = mask(v_start:v_end, h_start:h_end,:);

            % make the mask list
            mask_li = {};
            mask_li{1} = m;

            %op_test_pipe(input_piece, fake_output_for_test);
            op_test_pipe_with_mask(input_piece, mask_li, fake_output_for_test);
            output_piece = mem.output;

            output(v_start+(size_differ(1)/2):v_end-(size_differ(1)/2), h_start+(size_differ(2)/2):h_end-(size_differ(2)/2),:) = ...
                    output(v_start+(size_differ(1)/2):v_end-(size_differ(1)/2), h_start+(size_differ(2)/2):h_end-(size_differ(2)/2),:) + output_piece;
            count(v_start+(size_differ(1)/2):v_end-(size_differ(1)/2), h_start+(size_differ(2)/2):h_end-(size_differ(2)/2)) = ...
                    count(v_start+(size_differ(1)/2):v_end-(size_differ(1)/2), h_start+(size_differ(2)/2):h_end-(size_differ(2)/2)) + 1;
        end
    end
    count = max(count, 1);
    output = bsxfun(@rdivide, output, count);
end

author: jimmy-ren


text:
function varargout = plot_directed_edges(V,E,varargin)
  % PLOT_DIRECTED_EDGES Plot direct edges using quiver
  %
  % p = plot_directed_edges(V,E)
  % p = plot_directed_edges(V,E,LINESPEC,...)
  %
  % Inputs:
  %   V  #V by dim list of edge vertiecs
  %   E  #E by 2 list of edges
  % Outputs:
  %   p  plot handle
  %
  % See also: plot_edges
  %
  if isempty(E)
    return;
  end

  switch size(V,2)
  case 2
    q = quiver( ...
      V(E(:,1),1), ...
      V(E(:,1),2), ...
      V(E(:,2),1)-V(E(:,1),1), ...
      V(E(:,2),2)-V(E(:,1),2), ...
      0, ...
      varargin{:});
   otherwise
    error('Unsupported dimension: %d',size(V,2));
   end

  if nargout >= 1
    varargout{1} = p;
  end

end

author: alecjacobson

repo: gptoolbox


text:
function prob_test035 ( )

%*****************************************************************************80
%
%% TEST035 tests CHI_SQUARE_NONCENTRAL_MEAN, ***_SAMPLE, ***_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST035\n' );
  fprintf ( 1, '  For the noncentral chi square PDF:\n' );
  fprintf ( 1, '  CHI_SQUARE_NONCENTRAL_SAMPLE samples.\n' );

  a = 3.0;
  b = 2.0;

  check = chi_square_noncentral_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST035 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = chi_square_noncentral_mean ( a, b );
  variance = chi_square_noncentral_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );
  fprintf ( 1, '  PDF mean =        %14f\n', mean );
  fprintf ( 1, '  PDF variance =    %14f\n', variance );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial seed =    %12d\n', seed );

  for i = 1 : nsample
    [ x(i), seed ] = chi_square_noncentral_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '  Final seed =      %12d\n', seed );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function geometry_test035 ( )

%*****************************************************************************80
%
%% TEST035 tests LINE_IMP_POINT_DIST_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  ntest = 3;

  atest = [ 2.0, 2.0, 2.0 ];
  btest = [ 5.0, 5.0, 5.0 ];
  ctest = [ 3.0, 3.0, 3.0 ];
  ptest = [ ...
    0.0, 6.0; ...
    0.0, 5.0; ...
    0.0, 4.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST035\n' );
  fprintf ( 1, '  LINE_IMP_POINT_DIST_2D finds the distance from\n' );
  fprintf ( 1, '    a point X Y to a line A * X + B * Y + C = 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   X       Y       A       B       C       DIST\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    a = atest(i);
    b = btest(i);
    c = ctest(i);
    p(1:dim_num) = ptest(1:dim_num,i);

    dist = line_imp_point_dist_2d ( a, b, c, p );

    fprintf ( 1, '  %8f  %8f  %8f  %8f  %8f  %8f\n', ...
      p(1:dim_num), a, b, c, dist );

  end

  return
end


%% NISwGSP-denny
% %{
imfolder = 'images\NISwGSP-denny';
im_n = 15;
imfile = cell(im_n,1);
for ii = 1:im_n
    imfile{ii} = sprintf('%s\\denny%02d.jpg', imfolder, ii-1);
end

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = zeros(im_n,2);
ei = 0;
for ii = 1:im_n-1
    ei = ei + 1;
    edge_list(ei,:) = [ii,ii+1];
end
edge_list(im_n,:) = [im_n,1];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 0, 'equi', 0.02, imfolder );
%}


text:
% ONE DIMENSIONAL EFG PROGRAM
% SET UP NODAL COORDINATES ALONG BAR, DETERMINE NUMBER OF CELLS
x = [0.0:.1:1.0];
nnodes = length(x);
ncells = nnodes-1;
nosDeDirichlet = 1;

% SET PARAMETERS FOR WEIGHT FUNCTION, MATERIAL PROPERITES
dmax = 2.0;
E=1.0; area=1.0;

% DETERMINE DMI FOR EACH NODE
dm = dmax*(x(2)-x(1))*ones(1,nnodes);

%SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL
gg = zeros(1,ncells);
jac = (x(2)-x(1))/2;
weight = 2;
gg = -.05:.1:0.95; gg(1) = 0.0;

% INITIALIZE MATRICES
k = zeros(nnodes);
f = zeros(nnodes,1);
GG = zeros(nnodes,nosDeDirichlet);

% LOOP OVER GAUSS POINTS
for j = 1:length(gg)
   xg = gg(j);

% DETERMINE DISTANCE BETWEEN NODES AND GAUSS POINT
dif = xg*ones(1,nnodes)-x;

% SET UP WEIGHTS W AND DW FOR EACH NODE
clear w dw
for i=1:nnodes
drdx = sign(dif(i))/dm(i);
r = abs(dif(i))/dm(i);
if r<=0.5
  w(i) = (2/3) - 4*r*r + 4*r^3;
  dw(i) = (-8*r + 12*r^2)*drdx;
elseif r<=1.0
  w(i) = (4/3)-4*r+4*r*r -(4/3)*r^3;
  dw(i) = (-4 + 8*r-4*r^2)*drdx;
elseif r>1.0
w(i) = 0.0;
dw(i) = 0.0;
end
end

%SET UP SHAPE FUNCTIONS AND DERIVATIVES
won = ones(1,nnodes);
p = [won;x];
B = p.*[w;w];
pp = zeros(2);
A = zeros(2);
dA = zeros(2);
for i=1:nnodes
   pp = p(1:2,i)*p(1:2,i)';
   A = A+w(1,i)*pp;
   dA = dA+dw(1,i)*pp;
end
Ainv = inv(A);
pg = [1 xg];
phi(j,:) = pg*Ainv*B;
db = p.*[dw;dw];
da = -Ainv*(dA*Ainv);
dphi(j,:) = [0 1]*Ainv*B+pg*(da*B+Ainv*db);

%ASSEMBLE DISCRETE EQUATIONS
if j == 1
   GG(1:3,1) = -phi(1:3)';
end
%else
    k = k+(weight*E*area*jac)*(dphi(j,:)'*dphi(j,:));
    fbody = area*xg;
    f = f+(weight*fbody*jac)*phi(j,:)';
%end
end

% ENFORCE BOUNDARY CONDITIONS USING LAGRANGE MULTIPLIERS
q = [0];
m = [k GG;GG' zeros(1)];

% SOLVE FOR NODAL PARAMETERS
d = m\[f' q]';
u = d(1:nnodes);

result = phi*u;
for i=1:length(x)
    exact(i) =  0.5*(x(i)) - x(i)^3/6.0;
end
%%
plot(x,u);
title('valores u');
figure();
plot(x,result);
title('valores u*phi');
figure();
plot(x,exact);
title('valores exatos');

repo: Aero-matlab


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: rate_ci.m
% computes confidence interval estimate rates.
%   x = number of events (numerator of r_hat)
%   A = test area (denominator of p_hat) - the CI will be in the same units as the rate,
%       which is in the same units as 1/A.
%   alpha = probability that the true value will fall outside the confidence interval
%   method =
%       1 for 1-sided confidence interval,
%          note that you get both the lower 1-sided and the upper 1-sided when method = 1.
%       2 for 2-sided, minimum length interval (default)
%       3 for 2-sided, symmetrical about the estimate (except when either end of the interval would fall outside the acceptable range [0,1]
%       4 for 2-sided, equal tail masses
%       5 for 2-sided, MatLab implementation
%       6 for 2-sided, Normal approximation
%   ci = (r_hat, lower bound of confidence interval, upper bound of confidence interval)
%       note that you get both the lower 1-sided and the upper 1-sided when sides = 1.
%
% 010129 tdr created from interval_est_r and prob_ci to simply input parameters
% 010315 tdr added method options > 2 and verbose.
% 010326 tdr added x=0 fix for Method 5
% 010412 tdr added method to verbose output
% 030110 tdr added warnings for extremes x's and alpha's
% 031223 tdr added test for NaN in inputs

function metrics = rate_ci(x,A,alpha,method, verbose)

% ------------------------------------------------------------
% start checking inputs

if nargin == 4,
    verbose = 0;
end

if nargin == 3,
    verbose = 0;
    method = 2;
end

if (nargin ~= 3) & (nargin ~= 4) & (nargin ~= 5),
    error('Requires 3, 4, or 5 input arguments x, A, alpha, method(default=2), verbose(default=0)');
end

% inputs must all be scalars
if (max(max([size(x); size(A); size(alpha)])) ~= 1)
    error('Non-scalar input');
end;

if isnan(x) | isnan(A) | isnan(alpha)
    warning('NaN input')
    metrics = [NaN, NaN, NaN];
    return
end;

% x must be integer
if (round(x) ~= x)
    x = round(x);
    warning('Non-integer input x');
end;

% A must be > 0
if (A <= 0),
   metrics = [ NaN NaN NaN];
   warning('A <= 0');
   return;
end;

% x must be >= 0
if ( x < 0),
   metrics = [ NaN NaN NaN];
   warning('x < 0');
   return;
end;

% results may be inaccurate for x too large or alpha too small
if ( x > 1e6), warning('x > 1e6, results may not be accurate.'); end;
if ( alpha < 1e-5), warning('alpha < 1e-5, results may not be accurate.'); end;

% alpha must be > 0.0 and < 1.0
if ( alpha < 0 | alpha > 1),
   metrics = [ NaN NaN NaN];
   warning('alpha < 0 or alpha > 1');
   return;
end;

% end checking inputs
% ------------------------------------------------------------

% call interval estimator with best method
tic;
switch method
case 1,
   ci = interval_est_r(x,A,alpha,'cs1_1s');
case 2,
   ci = interval_est_r(x,A,alpha,'ml');
case 3,
   ci = interval_est_r(x,A,alpha,'cs1');
case 4,
   ci = interval_est_r(x,A,alpha/2,'cs1_1s');
case 5,
   [lambda ci_tmp] = poissfit(x,alpha);
   ci = [lambda ci_tmp(1) ci_tmp(2)]/A;
   if (x == 0) ci(2) = 0; end;
case 6,
   ci = interval_est_r(x,A,alpha,'na');
otherwise
   error('Not a valid method: 1 (1-sided), 2 (min length), 3 (symm.), 4 (equal tail), 5 (MatLab-poissfit), 6 (Normal Approx)');
end;
rt = toc;

if verbose
    max_lambda = max(1000, x*1000);
    length = ci(3)-ci(2);
    lower_tail = integrate_poisspdf(x,0,A*ci(2));
    upper_tail = integrate_poisspdf(x,A*ci(3),max_lambda);
    if (method == 1)
        actual_alpha = (lower_tail + upper_tail)/2;
    else
        actual_alpha = lower_tail + upper_tail;
    end;
    if (abs(actual_alpha - alpha) > 0.00005) % close_enough is 0.00001 throughout, but spec is 0.0001.
        warning('Interval not to spec: abs(desired_alpha - actual_alpha) < 0.00005)');
        warning_stats = [alpha actual_alpha alpha-actual_alpha]
    end;
    disp('r_hat, Lower CI Bound, Upper CI Bound, x, A, Desired alpha, Method, Length, Lower Tail, Upper Tail, Actual alpha, Delta alpha, Run Time')
    metrics = [ci x A alpha method length lower_tail upper_tail actual_alpha (alpha - actual_alpha) rt];
    disp(sprintf('%8.4g, %8.4g, %8.4g, %d, %8.4g, %8.6g, %d, %8.4g, %8.6g, %8.6g, %8.6g, %8.6g, %8.2g',metrics));
else
    metrics = ci;
end;

author: Sable

repo: mcbench-benchmarks


text:
function [fi,ps] = VF(r,T,p)
% VF Conversion of relative humidity to volume fraction of water vapor
%
% This program transfers relative humidity, absolute temperature and
% atmospheric pressure to volume fraction of water vapor and saturation
% vapour pressure using Goff-Gratch equation over plane surfaces of pure
% ice. Relation is based on integration of Clausius-Clapeyron equation and
% on experimental data. Stated range of equation validity is -100 to 0 °C.
%
% [fi,ps] = VF(r,T,p)
%
% Inputs:
%   r is the relative humidity (%). r may be a vector or a scalar
%     variable in the interval [0,100].
%
%   T means the absolute temperature (degrees K). T must have
%     the same number of elements as r. A warning is generated
%     if any elements fall outside the range [173,273].
%
%   p is the atmospheric pressure (kPa). p must have the same number
%     of elements as r. All elements must be non-negative.
%
% Results:
%   fi is the volume fraction of water vapor (%).
%
%   ps is the saturation vapour pressure (kPa).
%
% Example: At an rh of 35 %, Temp = 265 °K, atmospheric pressure = 98.9 kPa
%  [fi,ps] = VF(0.35,265,98.9)
%
%   fi =
%       0.00108
%
%   ps =
%       0.305
%
% For citation ensue following format: Malec L. (2007). VF: Conversion
% of relative humidity to volume fraction. A MATLAB file. [WWW document].
% URL http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
% objectId=14165
%
% List R. J.: Smithsonian Meteorological Tables. Smithsonian Institution,
% Washington 1984.
%
% Copyright. February 26, 2007.

if any(r(:)<0) || any(r(:)>100)
  error 'Relative humidity(r) must be percent'
end
if any(T(:)<=0)
  error 'Temperature(T) must be positive (°K)'
elseif any(T(:)<173) || any(T(:)>273)
  warning 'Temperature(T) fell outside the working range of this relation'
end
if any(p(:)<0)
  error 'Pressure(p) must be positive (kPa)'
end

a1 = -9.09718;
a2 = -3.56654;
a3 = 0.876793;

% Induction of ice point and saturation vapour pressure.
To = 273.16;
eo = 6.1071;

% Saturation vapour pressure and water fraction computation.
ps = 10.^(a1 * (To ./ T - 1) + a2 * log10(To ./ T) + a3 * (1 - T / To) ...
    + log10(eo));
ps = ps / 10;
fi = r .* ps ./ p;


text:
function risultati = fhme(net, nodes_info, data, n)
%HMEFWD	Forward propagation through an HME model
%
% Each row of the (n x class_num) matrix 'risultati' containes the estimated class posterior prob.
%
% ----------------------------------------------------------------------------------------------------
% -> pierpaolo_b@hotmail.com   or   -> pampo@interfree.it
% ----------------------------------------------------------------------------------------------------
%
ns=net.node_sizes;
if nargin==3
    ndata=n;
else
    ndata=size(data, 1);
end
altezza=size(ns,2);
coeff=cell(altezza-1,1);
for m=1:ndata
    %- i=2 --------------------------------------------------------------------------------------
    s=struct(net.CPD{2});
    if nodes_info(1,2)==0,
        mu=[]; W=[]; predict=[];
        mu=s.mean(:,:);
        W=s.weights(:,:,:);
        predict=mu(:,:)+W(:,:,:)*data(m,:)';
        coeff{1,1}=predict';
    elseif nodes_info(1,2)==1,
        coeff{1,1}=fglm(s.glim{1}, data(m,:));
    else,
        coeff{1,1}=fmlp(s.mlp{1}, data(m,:));
    end
    %----------------------------------------------------------------------------------------------
    if altezza>3,
        for i=3:altezza-1,
            s=[]; f=[]; dpsz=[];
            f=family(net.dag,i); f=f(2:end-1); dpsz=prod(ns(f));
            s=struct(net.CPD{i});
            for j=1:dpsz,
                if nodes_info(1,i)==1,
                    coeff{i-1,1}(j,:)=coeff{i-2,1}(1,j)*fglm(s.glim{j}, data(m,:));
                else
                    coeff{i-1,1}(j,:)=coeff{i-2,1}(1,j)*fmlp(s.mlp{j}, data(m,:));
                end
            end
            app=cat(2, coeff{i-1,1}(:)); coeff{i-1,1}=app'; clear app;
        end
    end
    %- i=altezza ----------------------------------------------------------------------------------
    if altezza>2,
        i=altezza;
        s=[]; f=[]; dpsz=[];
        f=family(net.dag,i); f=f(2:end-1); dpsz=prod(ns(f));
        s=struct(net.CPD{i});
        if nodes_info(1,i)==0,
            mu=[]; W=[];
            mu=s.mean(:,:);
            W=s.weights(:,:,:);
        end
        for j=1:dpsz,
            if nodes_info(1,i)==0,
                predict=[];
                predict=mu(:,j)+W(:,:,j)*data(m,:)';
                coeff{i-1,1}(j,:)=coeff{i-2,1}(1,j)*predict';
            elseif nodes_info(1,i)==1,
                coeff{i-1,1}(j,:)=coeff{i-2,1}(1,j)*fglm(s.glim{j}, data(m,:));
            else
                coeff{i-1,1}(j,:)=coeff{i-2,1}(1,j)*fmlp(s.mlp{j}, data(m,:));
            end
        end
    end
    %----------------------------------------------------------------------------------------------
    risultati(m,:)=sum(coeff{altezza-1,1},1);
    clear coeff; coeff=cell(altezza-1,1);
end
return

%-------------------------------------------------------------------

function [y, a] = fglm(net, x)
%GLMFWD	Forward propagation through 1-layer net->GLM statistical model

ndata = size(x, 1);

a = x*net.w1 + ones(ndata, 1)*net.b1;

nout = size(a,2);
% Ensure that sum(exp(a), 2) does not overflow
maxcut = log(realmax) - log(nout);
% Ensure that exp(a) > 0
mincut = log(realmin);
a = min(a, maxcut);
a = max(a, mincut);
temp = exp(a);
y = temp./(sum(temp, 2)*ones(1,nout));

%-------------------------------------------------------------------

function [y, z, a] = fmlp(net, x)
%MLPFWD	Forward propagation through 2-layer network.

ndata = size(x, 1);

z = tanh(x*net.w1 + ones(ndata, 1)*net.b1);
a = z*net.w2 + ones(ndata, 1)*net.b2;
temp = exp(a);
nout = size(a,2);
y = temp./(sum(temp,2)*ones(1,nout));

%-------------------------------------------------------------------

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef ImgHash < handle
    %IMGHASH  Base class for Image Hashing algorithms
    %
    % This module brings implementations of different image hashing algorithms.
    % It provides algorithms to extract the hash of images and fast way to
    % figure out most similar images in huge data set.
    %
    % ### Supported Algorithms
    %
    % - Average hash (also called Different hash)
    % - PHash (also called Perceptual hash)
    % - Marr Hildreth Hash
    % - Radial Variance Hash
    % - Block Mean Hash (modes 0 and 1)
    % - Color Moment Hash (this is the one and only hash algorithm resist to
    %   rotation attack (-90~90 degree)
    %
    % You can study more about image hashing from the paper and websites
    % listed in the references section.
    %
    % ### Performance under different attacks
    %
    % ![Performance chart](https://docs.opencv.org/3.3.0/attack_performance.JPG)
    %
    % ### Speed comparison with PHash library (100 images from ukbench)
    %
    % ![Hash Computation chart](https://docs.opencv.org/3.3.1/hash_computation_chart.JPG)
    % ![Hash comparison chart](https://docs.opencv.org/3.3.1/hash_comparison_chart.JPG)
    %
    % As you can see, hash computation speed of img_hash module outperform
    % [PHash library](http://www.phash.org/) a lot.
    %
    % PS: I do not list out the comparison of Average hash, PHash and Color
    % Moment hash, because I cannot find them in PHash.
    %
    % ### Motivation
    %
    % Collects useful image hash algorithms into OpenCV, so we do not need to
    % rewrite them by ourselves again and again or rely on another 3rd party
    % library (ex: PHash library). BOVW or correlation matching are good and
    % robust, but they are very slow compare with image hash, if you need to
    % deal with large scale CBIR (content based image retrieval) problem,
    % image hash is a more reasonable solution.
    %
    % ### More info
    %
    % You can learn more about img_hash modules from following links, these
    % links show you how to find similar image from ukbench dataset, provide
    % thorough benchmark of different attacks (contrast, blur, noise
    % (gaussion, pepper and salt), jpeg compression, watermark, resize).
    %
    % * [Introduction to image hash module of OpenCV](http://qtandopencv.blogspot.my/2016/06/introduction-to-image-hash-module-of.html)
    % * [Speed up image hashing of OpenCV (img_hash) and introduce color moment hash](http://qtandopencv.blogspot.my/2016/06/speed-up-image-hashing-of-opencvimghash.html)
    %
    % ### Contributors
    %
    % [Tham Ngap Wei](mailto:thamngapwei@gmail.com)
    %
    % ## References
    % [lookslikeit]:
    % > Neal Krawetz.
    % > [Looks Like It](http://www.hackerfactor.com/blog/?/archives/432-Looks-Like-It.html).
    %
    % [tang2012perceptual]:
    % > Zhenjun Tang, Yumin Dai, Xianquan Zhang. "Perceptual Hashing for Color
    % > Images Using Invariant Moments. Appl. Math, 6(2S):643S-650S, 2012.
    %
    % [zauner2010implementation]:
    % > Christoph Zauner. "Implementation and benchmarking of perceptual image
    % > hash functions". 2010.
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    %% Constructor/destructor
    methods
        function this = ImgHash(alg, varargin)
            %IMGHASH  Constructor
            %
            %     obj = cv.ImgHash(alg)
            %     obj = cv.ImgHash(alg, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __alg__ image hash algorithm, one of:
            %   * __AverageHash__ Computes average hash value of the input
            %     image. This is a fast image hashing algorithm, but only work
            %     on simple case. For more details, please refer to
            %     [lookslikeit].
            %   * __BlockMeanHash__ Image hash based on block mean. See
            %     [zauner2010implementation] for details.
            %   * __ColorMomentHash__ Image hash based on color moments. See
            %     [tang2012perceptual] for details.
            %   * __MarrHildrethHash__ Marr-Hildreth Operator Based Hash,
            %     slowest but more discriminative. See
            %     [zauner2010implementation] for details.
            %   * __PHash__ Slower than average_hash, but tolerant of minor
            %     modifications. This algorithm can combat more variation than
            %     than AverageHash, for more details please refer to
            %     [lookslikeit].
            %   * __RadialVarianceHash__ Image hash based on Radon transform.
            %     See [tang2012perceptual] for details.
            %
            % ## Options
            % The following are options for the various algorithms:
            %
            % ### `BlockMeanHash`
            % * __Mode__ block mean hash mode. default 'Mode0'
            %
            % ### `MarrHildrethHash`
            % * __Alpha__ scale factor for marr wavelet. default 2
            % * __Scale__ level of scale factor. default 1
            %
            % ### `RadialVarianceHash`
            % * __Sigma__ Gaussian kernel standard deviation. default 1
            % * __NumOfAngleLine__ Number of angles to consider. default 180
            %
            % See also: cv.ImgHash.compute, cv.ImgHash.compare
            %
            this.id = ImgHash_(0, 'new', alg, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.ImgHash
            %
            if isempty(this.id), return; end
            ImgHash_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = ImgHash_(this.id, 'typeid');
        end
    end

    %% ImgHashBase
    methods
        function hash = compute(this, img)
            %COMPUTE  Computes hash of the input image
            %
            %     hash = obj.compute(img)
            %
            % ## Input
            % * __img__ input image wanting to compute its hash value.
            %
            % ## Output
            % * __hash__ hash of the image.
            %
            % See also: cv.ImgHash.compare
            %
            hash = ImgHash_(this.id, 'compute', img);
        end

        function val = compare(this, hashOne, hashTwo)
            %COMPARE  Compare two hash values
            %
            %     val = obj.compare(hashOne, hashTwo)
            %
            % ## Input
            % * __hashOne__ Hash value one.
            % * __hashTwo__ Hash value two.
            %
            % ## Output
            % * __val__ indicate similarity between the two hashes, the
            %   meaning of the value vary from algorithm to algorithm.
            %
            % See also: cv.ImgHash.compute
            %
            val = ImgHash_(this.id, 'compare', hashOne, hashTwo);
        end
    end

    %% Util functions
    methods (Static)
        function hash = averageHash(img)
            %AVERAGEHASH  Calculates the average hash in one call
            %
            %     hash = cv.ImgHash.averageHash(img)
            %
            % ## Input
            % * __img__ input image want to compute hash value, type should be
            %   `uint8` with 1/3/4 channels.
            %
            % ## Output
            % * __hash__ Hash value of input, it will contain 16 hex decimal
            %   number, return type is `uint8`
            %
            % See also: cv.ImgHash.ImgHash, cv.ImgHash.compute
            %
            hash = ImgHash_(0, 'averageHash', img);
        end

        function hash = blockMeanHash(img, varargin)
            %BLOCKMEANHASH  Computes block mean hash of the input image
            %
            %     hash = cv.ImgHash.blockMeanHash(img)
            %     hash = cv.ImgHash.blockMeanHash(img, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ input image want to compute hash value, type should be
            %   `uint8` with 1/3/4 channels.
            %
            % ## Output
            % * __hash__ Hash value of input, it will contain 16 hex decimal
            %   number, return type is `uint8`.
            %
            % ## Options
            % * __Mode__ block mean hash mode, one of:
            %   * __Mode0__ (default) use fewer blocks and generates 16*16/8
            %     `uint8` hash values.
            %   * __Mode1__ use block blocks (step_sizes/2) and generates
            %     `fix(31*31/8)+1` `uint8` hash values.
            %
            % See also: cv.ImgHash.ImgHash, cv.ImgHash.compute
            %
            hash = ImgHash_(0, 'blockMeanHash', img, varargin{:});
        end

        function hash = colorMomentHash(img)
            %COLORMOMENTHASH   Computes color moment hash of the input
            %
            %     hash = cv.ImgHash.colorMomentHash(img)
            %
            % ## Input
            % * __img__ input image want to compute hash value, type should be
            %   `uint8` with 1/3/4 channels.
            %
            % ## Output
            % * __hash__ 42 hash values with type `double`.
            %
            % The algorithm comes from the paper [tang2012perceptual].
            %
            % See also: cv.ImgHash.ImgHash, cv.ImgHash.compute
            %
            hash = ImgHash_(0, 'colorMomentHash', img);
        end

        function hash = marrHildrethHash(img, varargin)
            %MARRHILDRETHHASH  Computes average hash value of the input image
            %
            %     hash = cv.ImgHash.marrHildrethHash(img)
            %     hash = cv.ImgHash.marrHildrethHash(img, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ input image want to compute hash value, type should be
            %   `uint8` with 1/3/4 channels.
            %
            % ## Output
            % * __hash__ Hash value of input, it will contain 16 hex decimal
            %   number, return type is `uint8`.
            %
            % ## Options
            % * __Alpha__ scale factor for marr wavelet. default 2
            % * __Scale__ level of scale factor. default 1
            %
            % See also: cv.ImgHash.ImgHash, cv.ImgHash.compute
            %
            hash = ImgHash_(0, 'marrHildrethHash', img, varargin{:});
        end

        function hash = pHash(img)
            %PHASH  Computes pHash value of the input image
            %
            %     hash = cv.ImgHash.pHash(img)
            %
            % ## Input
            % * __img__ input image want to compute hash value, type should be
            %   `uint8` with 1/3/4 channels.
            %
            % ## Output
            % * __hash__ Hash value of input, it will contain 8 `uint8` values.
            %
            % See also: cv.ImgHash.ImgHash, cv.ImgHash.compute
            %
            hash = ImgHash_(0, 'pHash', img);
        end

        function hash = radialVarianceHash(img, varargin)
            %RADIALVARIANCEHASH  Computes radial variance hash of the input image
            %
            %     hash = cv.ImgHash.radialVarianceHash(img)
            %     hash = cv.ImgHash.radialVarianceHash(img, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ input image want to compute hash value, type should be
            %   `uint8`, with 1/3/4 channels.
            %
            % ## Output
            % * __hash__ Hash value of input, contains 40 `uint8` values.
            %
            % ## Options
            % * __Sigma__ Gaussian kernel standard deviation. default 1
            % * __NumOfAngleLine__ The number of angles to consider.
            %   default 180
            %
            % See also: cv.ImgHash.ImgHash, cv.ImgHash.compute
            %
            hash = ImgHash_(0, 'radialVarianceHash', img, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function pass = test_imag( pref )
% Test IMAG
if ( nargin == 0 )
    pref = chebfunpref;
end

tol = 100*pref.cheb2Prefs.chebfun2eps;

f = chebfun2(@(x,y) cos(x.*y));
g = chebfun2(@(x,y) sin(x+y.^2));

% Simple consistency check:
x = linspace(-1,1,3);
[xx, yy] = meshgrid(x);
h = f + 1i *g;
pass(1) = norm( imag( feval(h,xx,yy) ) - feval(g,xx,yy) ) < 10*tol;
pass(2) = norm( imag( h ) - g ) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function value = r8_si ( x )

%*****************************************************************************80
%
%% R8_SI computes an approximation to the value of the sine integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real X, the argument of the sine integral.
%
%    Output, real VALUE, the value of the sine integral.
%
  nsics = 12;
  sics = [ ...
    -0.1315646598184841929, -0.2776578526973601892, ...
     0.0354414054866659180, -0.0025631631447933978, ...
     0.0001162365390497009, -0.0000035904327241606, ...
     0.0000000802342123706, -0.0000000013562997693, ...
     0.0000000000179440722, -0.0000000000001908387, ...
     0.0000000000000016670, -0.0000000000000000122 ];

  eta = 0.1 * eps;
  nsi = inits ( sics, nsics, 0.1 * eps );

  i = find ( abs ( x ) < sqrt ( eps ) );
  value(i) = x(i);

  j = find ( abs ( x ) <= 4.0 );
  value(j) = csevl ( ( x(j).^2 - 8.0 ) * 0.125, sics, nsi );
  value(j) = x(j) .* ( 0.75 + value(j) );

  k = find ( 4.0 < abs ( x ) );
  [ f, g ] = r89sifg ( abs ( x(k) ) );
  value(k) = 0.5 * pi - f .* cos ( abs ( x(k) ) ) - g .* sin ( x(k) );
  l = find ( x(k) < 0.0 );
  value(l) = - value(l);

  return
end


text:
classdef PWMoutputDataClass < DataBaseClass

    %% Public 'read-only' properties
    properties (SetAccess = private)
        FileNameAppendage = '_PWMoutput.csv';
        AX0 = [];
        AX2 = [];
        AX4 = [];
        AX6 = [];
    end

    %% Public methods
    methods (Access = public)
        function obj = PWMoutputDataClass(fileNamePrefix)
            data = obj.ImportCSVnumeric(fileNamePrefix);
            obj.AX0 = data(:,2);
            obj.AX2 = data(:,3);
            obj.AX4 = data(:,4);
            obj.AX6 = data(:,5);
        end
    end
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
function bk = r4_besks ( xnu, x, nin )

%*****************************************************************************80
%
%% R4_BESKS evaluates a sequence of K Bessel functions at X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real XNU, the order of the first function.
%    |XNU| < 1.
%
%    Input, real X, the argument.
%
%    Input, integer NIN, the absolute value of NIN indicates
%    the number of terms to compute.
%    If NIN < 0, successive values of NU count DOWN from XNU.
%    If NIN > 0, successive values of NU count UP from XNU.
%
%    Output, real BK(abs(NIN)), the K Bessel functions.
%
  persistent xmax

  if ( isempty ( xmax ) )
    xmax = - log ( r4_mach ( 1 ) );
    xmax = xmax + 0.5 * log ( 3.14 * 0.5 / xmax );
  end

  bk = r4_beskes ( xnu, x, nin );

  expxi = exp ( - x );
  n = abs ( nin );

  bk(1:n) = expxi * bk(1:n);

  return
end


text:
function [chi, labels] = imEuler3dEstimate(img, varargin)
% Estimate Euler number in a 3D image.
%
%   CHIest = imEuler3dEstimate(IMG)
%   CHIest = imEuler3dEstimate(IMG, CONN)
%   Estimate Euler number in a 3D image, without taking into account the
%   contribution of the voxels located on image border. The result of this
%   function is usually divided by the volume the sampling window to obtain
%   an estimate of Euler number density.
%
%   Example
%     imEuler3dEstimate
%
%   See also
%     imEuler3dDensity, imEuler3d

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRAE - Cepia Software Platform.

%% Process input arguments

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    labels = unique(img);
    labels(labels==0) = [];
    chi = zeros(length(labels), 1);
    for i = 1:length(labels)
        chi(i) = imEuler3dEstimate(img==labels(i), varargin{:});
    end
    return;
end

% extract connectivity
conn = 6;
if ~isempty(varargin)
    conn = varargin{1};
end

% determines connectivity to use on faces
conn2d = 4;
if conn == 26
    conn2d = 8;
end

% in case of binary image, compute only one label...
labels = 1;

%% Main processing

% Euler-Poincare Characteristic of the binary structure in image
chi     = imEuler3d(img, varargin{:});

% compute EPC on each of the 12 border edge of image, and keep the average
chix    = mean([ ...
    imEuler1d(img(:,   1,   1)) ...
    imEuler1d(img(:, end,   1)) ...
    imEuler1d(img(:,   1, end)) ...
    imEuler1d(img(:, end, end)) ...
    ]);
chiy    = mean([ ...
    imEuler1d(img(  1, :,   1)) ...
    imEuler1d(img(end, :,   1)) ...
    imEuler1d(img(  1, :, end)) ...
    imEuler1d(img(end, :, end)) ...
    ]);
chiz    = mean([ ...
    imEuler1d(img(  1,   1, :)) ...
    imEuler1d(img(end,   1, :)) ...
    imEuler1d(img(  1, end, :)) ...
    imEuler1d(img(end, end, :)) ...
    ]);

% compute EPC on each of the 6 border faces, and keep the average
chixy    = mean([ ...
    imEuler2d(squeeze(img(:, :,   1)), conn2d) ...
    imEuler2d(squeeze(img(:, :, end)), conn2d) ...
    ]);
chixz    = mean([ ...
    imEuler2d(squeeze(img(:,  1,  :)), conn2d) ...
    imEuler2d(squeeze(img(:, end, :)), conn2d) ...
    ]);
chiyz    = mean([ ...
    imEuler2d(squeeze(img(  1, :, :)), conn2d) ...
    imEuler2d(squeeze(img(end, :, :)), conn2d) ...
    ]);

% compute EPC on each of the 8 corners of image, and keep the average
chixyz   = mean([ ...
    img(  1,   1,   1), ...
    img(end,   1,   1), ...
    img(  1, end,   1), ...
    img(end, end,   1), ...
    img(  1,   1, end), ...
    img(end,   1, end), ...
    img(  1, end, end), ...
    img(end, end, end), ...
    ]);

% estimate EPC in image using mean edge correction
chi = chi - (chix + chiy + chiz) + (chixy + chixz + chiyz) - chixyz;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
%Program for SVD-based Image Quality Measure

%Program Description
% This program outputs a graphical & numerical image quality measure
% based on Singular Value Decomposition.
% For details on the implementation, please refer
% Aleksandr Shnayderman, Alexander Gusev, and Ahmet M. Eskicioglu,
% "An SVD-Based Grayscale Image Quality Measure for Local and Global Assessment",
% IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 15, NO. 2, FEBRUARY 2006.
%
%Author : Athi Narayanan S
%Student, M.E, EST,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%s_athi1983@yahoo.co.in
%http://sites.google.com/site/athisnarayanan/

%Sample Usage

clc;
close all;
clear all;

%Read Input Reference Gray Image
refImg=imread('lena256.bmp');
figure;subplot(121);
imshow(refImg);title('Orginal Image');

%Read Input Distorted Gray Image
distImg = imread('gaussian_noise_10.bmp');
subplot(122);
imshow(distImg);title('Distorted Image');

blkSize = 8;

[graMeasure, scaMeasure] = SVDQualityMeasure(refImg, distImg, blkSize);

figure;imshow(graMeasure);title('Graphical Measure');
disp('Numerical Measure');
disp(scaMeasure);


text:
function B = thresholdLocally(A, blksz, varargin)
% Performs LOCAL Otsu thresholding on an image; user can specify blocksize
%
% SYNTAX: B = thresholdLocally(A,PV_Pairs)
%
% THRESHOLDLOCALLY processes an image, calling graythresh on LOCAL
% blocks in an image. This facilitates easy thresholding of images with
% uneven background illumination, for which global thresholding is
% inadequate. Uses the Image Processing Toolbox function BLOCKPROC
% (R2009b).
%
% INPUTS:
%    A:   Any image (or path/name of an image) suitable for processing
%         with im2bw()
%
%    BLKSZ: Blocksize ([M,N]) with which to process the image.
%           DEFAULT: [32 32]
%
%    (OPTIONAL):
%    PV_Pairs: Any valid parameter-value pairs accepted by blockproc.
%           DEFAULTS:
%             'BorderSize':       [6 6]
%             'PadPartialBlocks': true
%             'PadMethod':        'replicate'
%             'TrimBorder':       true
%             'Destination':      [NOT SPECIFIED] (See BLOCKPROC for usage)
%    FudgeFactor: As an additional PV_Pair, one may enter:
%             'FudgeFactor',       value (DEFAULT = 1),
%              You may provide a scalar multiplier for the local value
%              returned by graythresh.
%
% OUTPUT:
%    B:     Output image (Unless 'DESTINATION' output is specified.)
%
% USAGE NOTE: To specify any PV_Pairs, BLKSZ must be provided as the second
% input. If the default value of blksz is desired, an empty bracket ([]) may be
% provided. (See Example 3.)
%
% EXAMPLES:
%
% % NOTE: All examples use image 'rice.png', which ships with the Image
% %       Processing Toolbox.
%
% img = imread('rice.png');
%
% % EXAMPLE 1) Default usage:
% thresholded = thresholdLocally(img);
% imshow(thresholded)
%
% % EXAMPLE 2) Specifying non-default blocksize:
% thresholded = thresholdLocally(img,[16 16]);
%
% % EXAMPLE 3) Specifying non-default padmethod (using default blocksize):
% thresholded = thresholdLocally(img,[],'PadMethod','symmetric');
%
% % EXAMPLE 4) Comparing and local thresholding, and thresholding after
% %            background normalization using tophat filtering:
% figure
% subplot(2,2,1);imshow(img);title('Original Image');
% tmp = im2bw(img,graythresh(img));
% subplot(2,2,2);imshow(tmp);title('Globally Thresholded');
% tmp = imtophat(img,strel('disk',15));
% tmp = im2bw(tmp,graythresh(tmp));
% subplot(2,2,3);imshow(tmp);title('Globally Thresholded after Tophat')
% tmp = thresholdLocally(img);
% subplot(2,2,4);imshow(tmp);title('Locally Thresholded');
%
% Written by Brett Shoelson, PhD.
% 12/17/2010
% Modifications:
% * 12/20/2010   Modified significantly to accept as optional inputs all
%   parameter-value pairs accepted by BLOCKPROC, as well as an additional
%   "fudge factor" parameter that allows one to scale the local graythresh
%   value by a scalar multiple.
% * 02/08/2011   Modified default blocksize to that calculated by bestblk
%
% Copyright 2010 The MathWorks
%
% See also: blockproc, graythresh, im2bw

if ~nargin
    error('THRESHOLDLOCALLY: Requires at least one input argument.')
end

if ischar(A)
    A = imread(A);
end

% DEFAULTS
% M                = 32;
% N                = 32;
[M, N] = bestblk(size(A));
opts.BorderSize  = [6 6];
opts.PadPartialBlocks = true;
opts.PadMethod   = 'replicate';
opts.TrimBorder  = true;
opts.Destination = [];
opts.FudgeFactor = 1;

if nargin > 1 && ~isempty(blksz)
    if numel(blksz) == 1
        M = blksz; N = M;
    elseif numel(blksz) == 2
        M = blksz(1);N = blksz(2);
    else
        error('THRESHOLDLOCALLY: Improper specification of blocksize parameter');
    end
end

if nargin > 2
    opts = parsePV_Pairs(opts,varargin);
end

fun = @(block_struct) im2bw(block_struct.data,...
    min(max(opts.FudgeFactor*graythresh(block_struct.data),0),1));

if isempty(opts.Destination)
    B = blockproc(A,[M N],fun,...
        'BorderSize',opts.BorderSize,...
        'PadPartialBlocks',opts.PadPartialBlocks,...
        'PadMethod',opts.PadMethod,...
        'TrimBorder',opts.TrimBorder);
else
    B = [];
    blockproc(A,[M N],fun,...
        'BorderSize',opts.BorderSize,...
        'PadPartialBlocks',opts.PadPartialBlocks,...
        'PadMethod',opts.PadMethod,...
        'TrimBorder',opts.TrimBorder,...
        'Destination',opts.Destination);
end

end

function opts = parsePV_Pairs(opts,UserInputs)
    ind = find(strcmpi(UserInputs,'BorderSize'));
    if ~isempty(ind)
        opts.BorderSize = UserInputs{ind + 1};
    end
    ind = find(strcmpi(UserInputs,'PadPartialBlocks'));
    if ~isempty(ind)
        opts.PadPartialBlocks = UserInputs{ind + 1};
    end
    ind = find(strcmpi(UserInputs,'PadMethod'));
    if ~isempty(ind)
        opts.PadMethod = UserInputs{ind + 1};
    end
    ind = find(strcmpi(UserInputs,'TrimBorder'));
    if ~isempty(ind)
        opts.TrimBorder = UserInputs{ind + 1};
    end
    ind = find(strcmpi(UserInputs,'Destination'));
    if ~isempty(ind)
        opts.Destination = UserInputs{ind + 1};
    end
    ind = find(strcmpi(UserInputs,'FudgeFactor'));
    if ~isempty(ind)
        opts.FudgeFactor = UserInputs{ind + 1};
    end
end


text:
function e = end(q,i,n)
% overloads end function

if n==1
  e = numel(q.a);
else
  e = size(q,i);
end

author: mtex-toolbox

repo: mtex

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


text:
%MAPPING_TASK Find task of actual mapping call
%
%   STATE = MAPPING_TASK(ARGIN,TASK)
%
% INPUT
%   ARGIN   - Arguments of a mapping call
%   TASK    - 'definition','fixed execution','training',
%             'trained execution'
%
% OUTPUT
%   STATE   - TRUE or FALSE
%
% DESCRIPTION
% This routine facilitates the parsing of calls to PRTools mapping
% routines.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, DEFINE_MAPPING, TRAINED_CLASSIFIER, TRAINED_MAPPING

% Copyright: Robert P.W. Duin, r.p.w.duin@37steps.com

function state = mapping_task(argin,task)

if nargin < 2
  error('No task specified')
end

state = false;
switch lower(task)
  case {'definition'}
    if isempty(argin) || (iscell(argin) && isempty(argin{1}))
      state = true;
    end
  case {'fixed execution','fixed_execution','fixedexecution'}
    if (nargin > 0) && (isdouble(argin{1}) || isa(argin{1},'prdataset')) && ...
        (numel(argin) == 1 || ~ismapping(argin{2}))
      state = true;
    end
  case ('training')
    if (nargin > 0) && (isdouble(argin{1}) || isa(argin{1},'prdataset')) && ...
        (numel(argin)==1 || ~ismapping(argin{2}) || isuntrained(argin{2}))
      state = true;
    end
  case {'execution','trained execution','trained_execution','trainedexecution'}
    if (nargin > 1) && (isdouble(argin{1}) || isa(argin{1},'prdataset')) && ...
        ismapping(argin{2})
      state = true;
    end
  case {'combiner'}
    if (nargin > 0) && (ismapping(argin{1}))
      state = true;
    end
  otherwise
    error('Wrongly specified mapping task')
end

return

author: marianux

repo: ecg-kit

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


text:
function [x_trans,IDX,ntrans] = ...
  ComputePerFrameTrans(x,trans_types)

if ~exist('trans_types','var'),
  %trans_types = 'all';
  trans_types = uint8(15);
end

%usealltrans = ischar(trans_types) && strcmpi(trans_types,'all');

x_trans = x;
IDX.orig = 1;
IDX.abs = 0;
IDX.flip = 0;
%if usealltrans || ismember('abs',trans_types),
if bitand(2,trans_types),
  x_trans(end+1,:) = abs(x);
  IDX.abs = size(x_trans,1);
end
%if usealltrans || ismember('flip',trans_types),
if bitand(4,trans_types),
  x_trans(end+1,:) = -x;
  IDX.flip = size(x_trans,1);
end
ntrans = size(x_trans,1);

author: kristinbranson

repo: JAABA


text:
function D=Data3
Coord=100*[1 1 0;0 1 0; 0 0 0];   % coordinates of nodes
Con=[1 2; 1 3 ];
Re=[0 0 1;1 1 1; 1 1 1];
Load=zeros(size(Coord)); Load(1,:)=[0 -1e5 0];
% or:   Load=[0 0 0;0 -1e5 0;0 0 0;0 -1e5 0;0 0 0;0 0 0];
E=ones(1,size(Con,1))*1e7;      %*1e7;  % Elasticity ( Youngs Modulous)
% or:   E=[1 1 1 1 1 1 1 1 1 1]*1e7;
A=[10 10 ];
D=struct('Coord',Coord','Con',Con','Re',Re','Load',Load','E',E','A',A');


text:
function [ep,out] = VBA_groupBMC_btwConds(L,options,factors)
% group-level between-condition Bayesian model comparison
% function Ltilde = VBA_groupBMC_btwConds(L,options,factors)
% IN:
%   - L: nmXnsxnc array of log-model evidences (nm models; ns subjects; nc
%   conditions).
%   - options: a structure containing the following fields:
%       .DisplayWin: flag for display window
%       .verbose: flag for summary statistics display
%       .families: a cell array of size nf, which contains the indices of
%       the models that belong to each of the nf families.
%   - factors: n1Xn2X...Xnf factorial condition attribution matrix. Each
%   entry contains the index of the corresponding condition. This is used in
%   case of a factorial design, to gather conditions into orthogonal
%   factors.
% OUT:
%   - ep: a nfX1 vector containing the exceedance probability of no
%   difference in models along each of the nf dimensions of the factorial
%   design space (if design is not factorial, then nf=1).
%   - out: structure containing the following fields:
%       .dt: the algorithm execution time (in sec)
%       .date: date vector, in matlab format (see clock.m)
%       .families: a cell array of size 2Xnf, which contains the indices of
%       the models that belong to each of the tuples families 'equal' and
%       'not equal', for each dimension of the factorial design.
%       .options: this is useful when using default options
%       .C: ncXnt tuple label matrix (where the number of tuples is
%       nt=nm^nc). Each column identifies a nc-tuple by the index of models
%       in each of the  possible conditions.
%       .VBA: a nfX1 structure containing the output of VBA_groupBMC.m,
%       which has been called for each dimension of the factorial design.

options.tStart = tic;
[nm,ns,nc] = size(L);
try;options;catch;options=[];end
if nc==1
    disp('Error: only one condition detected!')
    ep = [];
    out = [];
    return
end
if ~isfield(options,'verbose')
    options.verbose = 1;
end
if ~isfield(options,'DisplayWin')
    options.DisplayWin = 1;
end
if ~isfield(options,'families')
    options.families = [];
end
if ~isempty(options.families)
    nfam = length(options.families);
    Cfam = zeros(nm,1);
    tmp = [];
    for i=1:nfam
        indf = options.families{i};
        if isempty(indf)
            disp(['Error: family #',num2str(i),' is empty!'])
            ep = [];
            out = [];
            return
        end
        Cfam(indf) = i;
        if ~isempty(intersect(tmp,indf))
            disp('Error: families are not mutually exclusive!')
            ep = [];
            out = [];
            return
        end
        tmp = [tmp; VBA_vec(indf)];
    end
    if ~isequal(VBA_vec(unique(tmp)),VBA_vec(1:nm))
        if numel(unique(tmp)) < nm
            disp('Error: families do not cover the entire set of models!')
        else
            disp('Error: families contain models that do not exist!')
        end
        ep = [];
        out = [];
        return
    end
else
    Cfam = VBA_vec(1:nm);
end

try;factors;catch;factors=VBA_vec(1:nc);end
sf = size(factors);
sf(sf<=1) = [];
nf = size(sf,2); % number of factors/dimensions across conditions
indf = cell(nf,1);
for f=1:nf
    nlevels = size(factors,f); % #levels in this dimension
    nconds = numel(factors)/nlevels; % #conditions per level
    indf{f} = zeros(nconds,nlevels);
    strpar = repmat([':,'],1,nf);
    strpar(end) = [];
    strpar(2*(f-1)+1) = 'k';
    for k=1:nlevels
        eval(['indf{f}(:,k) = VBA_vec(factors(',strpar,'));'])
    end
end
% form nc-tuples
[Ccon] = VBA_getNtuples(nm,nc,0);
nt = size(Ccon,2);
fam = cell(2,nf);
Lt = zeros(nt,ns);
if options.verbose
    fprintf(1,['Forming tuples families... ']);
    fprintf(1,'%6.2f %%',0)
end
for i=1:nt
    Ci = Ccon(:,i);
    for j=1:nc
        Lt(i,:) = Lt(i,:) + L(Ci(j),:,j);
    end
    % pool conditions across dimension
    for f=1:nf
        nconds = size(indf{f},1);
        flag = 1;
        for k=1:nconds
            flag = flag && isequal(length(unique(Cfam(Ci(indf{f}(k,:))))),1);
        end
        if flag
            fam{1,f} = [fam{1,f};i];
        else
            fam{2,f} = [fam{2,f};i];
        end
    end
    if options.verbose
        fprintf(1,repmat('\b',1,8))
        fprintf(1,'%6.2f %%',100*i/nt)
    end
end
if options.verbose
    fprintf(1,repmat('\b',1,8))
    fprintf(' OK.')
    fprintf('\n')
end

% perform per-condition BMS
for i=1:nc
    VBA_disp({'---',['[per-condition BMS: condition #',num2str(i),']']},options)
    [out.VBA.cond(i).posterior,out.VBA.cond(i).out] = VBA_groupBMC(L(:,:,i),options);
end

% perform btw-condition BMS
opt.verbose = 0;
opt.DisplayWin = options.DisplayWin;
ep = zeros(nf,1);
out.pep = zeros(nf,1);
for f=1:nf
    if nf >1
        VBA_disp({'---',['[factorial design: dimension #',num2str(f),' (',num2str(size(factors,f)),' levels)]'],'---'},options)
    else
        VBA_disp({'---',['[btw-condition stability assessment]'],'---'},options)
    end
    opt.families = fam(:,f);
    [out.VBA.btw(f).posterior,out.VBA.btw(f).out] = VBA_groupBMC(Lt,opt);

    %-- display summary statistics
    ep(f) = out.VBA.btw(f).out.families.ep(1);
    out.pep(f) = ep(f).*(1-out.VBA.btw(f).out.bor) + 0.5*out.VBA.btw(f).out.bor;

    if options.verbose
        if floor(out.VBA.btw(f).out.dt./60) == 0
            timeString = [num2str(floor(out.VBA.btw(f).out.dt)),' sec'];
        else
            timeString = [num2str(floor(out.VBA.btw(f).out.dt./60)),' min'];
        end
        n1 = length(opt.families{1});
        n2 = length(opt.families{2});
        fprintf(['VB converged in ',num2str(length(out.VBA.btw(f).out.F)),' iterations (took ~',timeString,').','\n'])
        fprintf(['Dimensions:','\n'])
        fprintf(['     - subjects: n=',num2str(ns),'\n'])
        fprintf(['     - conditions: c=',num2str(nc),'\n'])
        fprintf(['     - models: K=',num2str(nm),'\n'])
        fprintf(['     - ',num2str(nc),'-tuples: t=',num2str(nt),' (',num2str(n1),'+',num2str(n2),')','\n'])
        fprintf(['Posterior probabilities:','\n'])
        fprintf(['     - RFX: p(H1|y)= ','%4.3f','\n'],1-out.VBA.btw(f).out.bor)
        fprintf(['     - null: p(H0|y)= ','%4.3f','\n'],out.VBA.btw(f).out.bor)
        fprintf(['Assessment of model stability across conditions:','\n'])
        fprintf(['     - exc. prob.: ','%4.3f','\n'],ep(f))
        fprintf(['     - protected exc. prob.: ','%4.3f','\n'],out.pep(f))
        fprintf('\n')
    end
end

% wrap-up
out.dt = toc(options.tStart);
out.date = clock;
out.options = options;
out.families = fam;
out.factors = factors;
out.Ccon = Ccon;
out.L = L;

if options.DisplayWin
    VBA_displayGroupBMCbtw(ep,out)
end

author: MBB-team

repo: VBA-toolbox


text:
% GP_COV_SE - Squared exponential covariance function for Gaussian
% processes.
%
% [K, DK_LOGTHETA, DK_X2] = GP_COV_SE(X1, X2, LOGTHETA)

% Last modified 2010-10-06
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function [K, dK_logtheta, dK_x2] = gp_cov_se(x1, x2, logtheta)

% 2-norm distances
if ~isempty(x2)
  % Distances for covariances
  n1 = cols(x1);
  n2 = cols(x2);
  D = zeros(n1,n2);
  for i=1:n1
    for j=1:n2
      D(i,j) = norm(x1(:,i) - x2(:,j));
    end
  end
else
  % Distances for variances
  n = rows(x1);
  D = zeros(n,1);
end

% Covariance matrix
D2 = D.^2;
K = exp(logtheta(1)) * exp(-0.5*exp(-2*logtheta(2))*D2);

% Gradient for hyperparameters
if nargout >= 2
  dK_logtheta = zeros([size(D), 2]);
  dK_logtheta(:,:,1) = K;
  dK_logtheta(:,:,2) = K .* (-0.5*D2*exp(-2*logtheta(2))) .* (-2);
end

% Gradients for inputs x2
if nargout >= 3
  if isempty(x2)
    error('Can''t calculate gradient: x2 not given');
  end
  d = rows(x2); % dimensionality of inputs
  m = cols(x1); % number of other inputs
  n = cols(x2); % number of inputs
  dK_x2 = zeros([d,m,n]);
  for i=1:m
    for j=1:n
      dK_x2(:,i,j) = K(i,j) * (-0.5*exp(-2*logtheta(2))) * 2*(x2(:,j)-x1(:,i));
    end
  end
end

function [K, dK] = gpK_sqexp(D, p1, p2)

% [K, dK] = gpK_sqexp(D, p1, p2)
%
% Squared exponential covariance function for GP.
%
% p1 is the scale
% p2 is the length
%
K = p1^2*exp(-0.5*D.^2/(p2^2));

if nargout >= 2
  dK = zeros([size(D), 2]);
  dK(:,:,1) = K .* 2 / p1;
  dK(:,:,2) = K .* (-0.5*D.^2) .* (-2*p2^(-3));
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [ pred,model,prob ] = classf_( methodName,X,Y,param,Xtest )
%CLASSF_ Classification, such as KNN, SVM, ...
%	[ pred,model,prob ] = classf_( methodName,X,Y,param,Xtest )
%   A wrapper to call classf_methodName_tr/te. Returns the predicted labels
%   PRED for XTEST and its confidence PROB, the MODEL of the classifier for later
%   test data. The parameters PARAM is optional. PROB may only be available
%   for some methods.
%	METHODNAME is case insensitive.

%   Example:
%		[pred,model,prob] =	classf_('svm',X,Y,struct('c',10),Xtest);

%	Copyright 2015 Ke YAN, Tsinghua Univ. http://yanke23.tk, xjed09@gmail.com

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

methodName = lower(methodName);
code = sprintf('model=classf_%s_tr(X,Y,param);', methodName);
eval(code);
code = sprintf('[pred,prob]=classf_%s_te(model,Xtest);', methodName);
eval(code);

end

author: viggin

repo: domain-adaptation-toolbox


text:
function [h,h2] = logzplot(varargin)
%Surface plot with log-scaled color and z-axis
%
% SYNTAX
%
%   LOGZPLOT
%   LOGZPLOT(Z)
%   LOGZPLOT(Z,C)
%   LOGZPLOT(X,Y,Z)
%   LOGZPLOT(X,Y,Z,C)
%   LOGZPLOT(...,'PLOTFUN')
%   LOGZPLOT(TriRep,'PLOTFUN')
%   LOGZPLOT(TRI,X,Y,Z,'PLOTFUN')
%   LOGZPLOT(TRI,X,Y,Z,C,'PLOTFUN')
%   LOGZPLOT(...,@PLOTFUN)
%   LOGZPLOT colorbar or LOGZPLOT(...,'colorbar')
%   LOGZPLOT(HANDLE,...)
%   H = LOGZPLOT(...)
%   [H,H2] = LOGZPLOT(...)
%
% DESCRIPTION
%
% LOGZPLOT creates a plot with logarithmic scaling of the z-axis and color
% data.  The plot type can be any of IMAGE, MESH, PCOLOR, SURF, TRISURF or
% TRIMESH, specified as a string or function handle.  The default plot type
% is SURF.  If called without data inputs, LOGZPLOT applies the log-scale
% transformation to an existing surface, patch or image object.  LOGZPLOT
% called with the optional argument 'colorbar' will also create a
% log-scaled colorbar.
%
% LOGZPLOT(Z)
% LOGZPLOT(Z,C)
% LOGZPLOT(X,Y,Z)
% LOGZPLOT(X,Y,Z,C)   Create a log-scaled surface plot using the data in Z
% by first calling SURF() and then transforming the resulting surface
% object.  The optional arguments X and Y specify values for the X and Y
% locations of the elements of Z. If given, X and Y must be specified
% together, and be of the appropriate size for the data in Z.  If the
% optional argument C is given, the surface's color will be based on the
% values in C instead of Z.
%
% LOGZPLOT(...,'PLOTFUN')
% LOGZPLOT(TriRep,'PLOTFUN')
% LOGZPLOT(TRI,X,Y,Z,'PLOTFUN')
% LOGZPLOT(TRI,X,Y,Z,C,'PLOTFUN')
% LOGZPLOT(...,@PLOTFUN)   Use the plotting function specified by the
% string or function handle in PLOTFUN to plot the data. By default, SURF
% is used to generate the plot. Supported plotting functions are SURF,
% MESH, PCOLOR, TRISURF and TRIMESH. The number and type of the data inputs
% depends on the selected plotting function.
%
% For IMAGE plots, the data inputs must be of the form (C) or (X,Y,C),
% where C is the data used to determine the pixel colors, and X and Y
% specify the scales of the x and y axes.  See the documentation for IMAGE
% for more information.
%
% For SURF and MESH plots, the data inputs must be of the form
% (Z), (X,Y,Z), or (X,Y,Z,C).  See the documentation for SURF and MESH for
% more information.
%
% For PCOLOR plots, the data inputs must be of the form (C) or (X,Y,C).
% See the documentation for PCOLOR for more information.
%
% For TRISURF and TRIMESH plots, the data must be input either as a TriRep
% object or as (TRI,X,Y,Z,C) where TRI is the triangulation of the data in
% X and Y, Z is the height data, and C is optional color data to be used
% instead of Z to color the plot.  See the documentation for TRISURF and
% TRIMESH for more information
%
% LOGZPLOT(...,'colorbar')   Additionally creates a log-scaled colorbar.
%
% LOGZPLOT(AX_HANDLE,...)   Uses the existing axes specified by AX_HANDLE
% to create the plot.
%
% H = LOGZPLOT(...)   Returns the handle of the surface, patch or image
% object that has been transformed.
%
% [H,H2] = LOGZPLOT(...)   In addition to the plot object handle, returns
% the handle of the colorbar.
%
% LOGZPLOT   Calling LOGZPLOT with no arguments will apply a log-scale
% transformation to a surface or patch object located in the current axes.
% The first such object found will be used by LOGZPLOT.  If the object has
% already been transformed by LOGZPLOT, calling LOGZPLOT again will have no
% effect.
%
% LOGZPLOT colorbar
% LOGZPLOT('colorbar')   Additionally create a log-scaled colorbar, or
% scale an existing colorbar attached to the axes of the log-transformed
% plot. The resulting colorbar scale should be accurate to the log-scaled
% data.  See the REMARKS section below for more information.
%
% LOGZPLOT(HANDLE)
% LOGZPLOT(HANDLE,'colorbar')   Transform the surface or patch object
% specified by HANDLE, or, if HANDLE refers to an axes, transform a surface
% or patch object located in the specified axes.
%
% REMARKS
%
% If LOGZPLOT is used to transform an existing surface or patch object, the
% color data ('CData' for surface objects, 'FaceVertexCData' for trisurf
% and trimesh patch objects) must be in indexed form.  LOGZPLOT can not
% transform objects with truecolor (RGB) color data, and will exit with a
% warning if called on an object with truecolor color data.  LOGZPLOT sets
% the object's 'CDataMapping' property to 'scaled'.
%
% LOGZPLOT replaces the original color data, CData, with a transformed
% version, log10(CData).  An additional linear scaling is performed on the
% transformed data so that the scale on a colorbar will have the same range
% as the original data.
%
% Changing the value of the axes 'CLim' (color limits) property using
% either the CAXIS command or set(ax_handle,'CLim') will alter the mapping
% of the data to the colormap as described in the documentation for the
% CAXIS command. This is a useful technique to highlight different data
% ranges in the plot.  See Example 2 below.
%
% However, as noted above, the log-transformed data does not result in an
% accurate colorbar scale without an additional (linear) transformation.
% To accurately map the color data to the colorbar's scale after a change
% to the axes 'CLim' property, LOGZPLOT attaches a set of listener
% functions to the axes that correct the color data scaling whenever the
% 'CLim' property is changed.  This allows the user to specify 'CLim'
% values in the same units and range as the original data.
%
% The listener functions were written for MATLAB R2010a, and may not
% function correctly on older releases due to changes in the MATLAB
% graphics system. If problems related to the listener functions occur, set
% the parameter 'listenerEnable' to false in the first section of the
% LOGZPLOT code, just below the help text.  This will disable the listener
% functionality and prevent rescaling of the color data after a CLim
% change.
%
% The scaling of the color data performed by LOGZPLOT introduces numerical
% error. The magnitude of the error depends on the range of the original
% data as well as the values of the axes' CLim property.  The error is
% cumulative, so making many changes to the CLim values can potentially
% result in inaccurate color data.
%
% If an object created or modified by LOGZPLOT is located in an axes with
% non-log-scaled indexed color objects (surface, patch, or image), the
% colorbar scale will not be accurate.  To ensure accurate colorbar scales,
% do not combine a LOGZPLOT-scaled object with non-scaled indexed-color
% objects in the same axes.
%
% LOGZPLOT attempts to avoid re-transforming a previously log-scaled
% plot (from a previous call to LOGZPLOT).  This allows multiple calls
% to LOGZPLOT using the same object, e.g. to add a colorbar after creating
% the plot.
%
% Because MATLAB's OpenGL renderer does not support logarithmic axes, the
% figure's 'Renderer' property must be set to another renderer for proper
% display of the plot.  MATLAB should change the renderer automatically. In
% some cases when using an existing figure for the plot output, the
% renderer will not change automatically, and the log-scale z-axis will not
% display properly.  If this occurs, set the renderer manually using:
%    set(fig_handle,'Renderer','ZBuffer')
% where fig_handle is the handle of the figure in question.  See the
% documentation for 'Figure Properties' for more information.
%
% Compared to other high-level plotting functions, TRIMESH (and to a lesser
% extent TRISURF) offer incomplete support for the specification of an axes
% handle as a target for the plot output in place of gca().  In particular,
% TRIMESH cannot accept an axes handle if the data is specified as a TriRep
% object.  If LOGZPLOT is called with a TriRep object and TRIMESH as the
% plotting function, the plot will be created in the current axes, ignoring
% any axes handle input.
%
% EXAMPLES
%
% % Example 1 - Compare linear and log-scaled surface plots
%
% % Generate some Gaussian data with a small sinusoidal component:
% x = linspace(-10,10,101);
% [X ,Y] = meshgrid(x);
% Z = 5*exp(-(0.82*(X+3.5).^2 + 0.46*(Y-3.5).^2)) + ...
%     0.8*exp(-(0.45*(X-1.5).^2 + 0.95*(Y-1.5).^2)) + ...
%     0.2*exp(-(0.75*(X-4).^2 + 0.85*(Y+1).^2)+(0.7*(X-3)+.3*(Y-1)).^2)+...
%     0.06*exp(-(0.2*(X+2.5).^2 + 0.3*(Y+3.5).^2)) + ...
%     -0.45*exp(-(0.5*(X+3.3).^2 + 1.5*(Y-2.5).^2)) + ...
%     0.0015*sin(2.6*X+1.1*Y-0.2*X.*Y) + ...
%     0.0009*sin(1.3*X+2.1*Y+0.12*X.*Y);
%
% % Scale the data so its range is 1 to 30000:
% minz = 1;
% maxz = 3e4;
% Z = (maxz-minz)*(Z-min(Z(:)))./(max(Z(:))-min(Z(:))) + minz;
%
% % Add a large Gaussian that will swamp the rest of the data:
% Z = Z + 1e6*exp(-2*(X.^2 + (Y-4).^2));
%
% % Linear z-axis with linear color scale:
% figure(1)
% set(1,'position',[254 335 642 471])
% colormap(jet(64))
% h = surf(X,Y,Z); colorbar
% title('Linear Z-scale and Coloring')
% % Almost all of the surface's detail is hidden because both the height
% % scale and color scale are dominated by the single prominent feature of
% % the data.
%
% % Log-scale z-axis with linear color scale:
% % Use SURF, then change the z-axis scaling using set(gca,'ZScale','log').
% figure(2)
% set(2,'position',[254 335 642 471])
% colormap(jet(64))
% h = surf(X,Y,Z); colorbar
% set(gca,'ZScale','log')
% title('Log Z-scale with Linear Coloring')
% % Note how the log-scale z-axis improves the visibility of the small
% % elevation features of the surface. However, the coloring of the surface
% % is not ideal - the variation in color is still concentrated near the
% % top of the surface due to the linear color scale.
%
% % Log-scale plot using LOGZPLOT to achieve both log-scale z-axis and
% % log-scale color:
% figure(3)
% set(3,'position',[254 335 642 471])
% colormap(jet(64))
% logzplot(X,Y,Z,'colorbar')
% title('Logarithmic Z-scale and Coloring')
% % The plot created by LOGZPLOT shows the advantage of using log-scaled
% % color in addition to the log-scaled z-axis.
%
% % Example 2 - Multiple calls to LOGZPLOT, changing the 'CLim' property
%
% % (Using the X,Y,Z data from Example 1)
% figure(4)
% colormap(jet(64))
%
% % Make a surface using pcolor:
% pcolor(X,Y,Z); shading flat
%
% % Call logzplot to transform to the pcolor plot to log-scale:
% logzplot
%
% % Call logzplot again to add a log-scale colorbar:
% logzplot colorbar
%
% % Use the CAXIS command to set the axes CLim parameter (color limits) to
% % highlight the lower range of the data.  The colorbar scale should
% % adjust to the new limits:
% caxis([1 70])
%
% % Call CAXIS again with different CLim values, to highlight the middle
% % range of the data:
% caxis([40 8000])
%
% % One more call to CAXIS to reset the color limits:
% caxis auto
%
%
% See also SURF, MESH, PCOLOR, TRISURF, TRIMESH, CAXIS, IMAGE
%

% $$FileInfo
% $Filename: logzplot.m
% $Path: $toolboxroot/
% $Product Name: logzplot
% $Product Release: 1.2
% $Revision: 1.2.8
% $Toolbox Name: Custom Plots Toolbox
% $$
%
% Copyright (c) 2010-2012 John Barber.
%
% Release History:
% v 1.0 : 2010-Nov-08
%       - Initial release
% v 1.1 : 2010-Nov-30
%       - Added support for TRISURF and TRIMESH plots
%       - Improved memory performance
%       - Improved speed of axes CLim change listener function
%       - Improved colorbar support
% v 1.2 : 2012-May-30
%       - Added support for plots using IMAGE
%       - Added colorbar handle output
%       - Fixed CLim equality test bug
%       - Code cleanup (variable/function names changed, etc.)

%% Defaults and initial values

% Listener enable flag:
% Disable if listener function causes errors
% -- false = disabled
% -- true = enabled (default)
listenerEnable = true;
% listenerEnable = false;  % Uncomment this line to disable listener

% Supported plotting functions:
plotFunList = {'surf','pcolor','mesh','trisurf','trimesh','image'};

% Default values:
% hSurf = [];
% hAx = [];
% pcolorFlag = false;
inputHandle = [];
colorbarFlag = false;
plotFun = @surf;
x = [];
y = [];
z = [];
CData = [];
tr = [];
triRep = [];
badInput = false;
plotFlag = false;

%% Parse inputs
if ~isempty(varargin)

    % Get char arrays and make indices of other classes
    chars = varargin(cellfun(@ischar,varargin));
    numIdx = cellfun(@isnumeric,varargin);
    fHandleIdx = cellfun(@(x)(isa(x,'function_handle')),varargin);
    scalarIdx = cellfun(@isscalar,varargin);
    emptyIdx = cellfun(@isempty,varargin);
    cellIdx = cellfun(@iscell,varargin);
    triRepIdx = cellfun(@(x)(isa(x,'TriRep')),varargin);

    % Make an index of non-empty non-scalar arrays
    numIdx = numIdx & ~scalarIdx & ~emptyIdx;
    numIdx = find(numIdx);

    % Make a cell array with just the function handles
    fHandles = varargin(fHandleIdx);

    % Detect an HG handle if it was given
    handleList = cell2mat(varargin(scalarIdx & ~fHandleIdx & ...
        ~emptyIdx & ~cellIdx));
    inputHandle = handleList(find(ishghandle(handleList),1));
    inputHandle((end-length(inputHandle)+2):end) = [];

    % Test for 'colorbar' and 'plot function name' inputs
    k = 1;
    while k <= length(chars)
        switch lower(chars{k})
            case 'colorbar'
                colorbarFlag = true;
                k = k + 1;
            case plotFunList
                plotFun = str2func(chars{k});
                k = k + 1;
            otherwise
                k = k + 1;
        end
    end

    % If a valid function handle was input, take it
    if ~isempty(fHandles) && ...
            any(strcmp(func2str(fHandles{end}),plotFunList))
        plotFun = fHandles{end};
    end

    % Get the data to be plotted, otherwise, exit with an error
    switch func2str(plotFun)
        case 'image'
            % For image plots, valid inputs are 1 (CData) or 3 (x,y,CData)
            switch length(numIdx)
                case 1
                    z = varargin{numIdx};
                    x = 1:size(z,2);
                    y = 1:size(z,1);
                    plotFlag = true;
                case 3
                    x = varargin{numIdx(1)};
                    y = varargin{numIdx(2)};
                    z = varargin{numIdx(3)};
                    plotFlag = true;
                otherwise
                    badInput = true;
            end

        case {'surf','mesh','pcolor'}
            % For surface plots, valid inputs are 1 (z), 2 (z,CData),
            % 3 (x,y,z), or 4 (x,y,z,CData)
            switch length(numIdx)
                case 0
                    %
                    plotFlag = false;
                case 1
                    z = varargin{numIdx};
                    x = 1:size(z,2);
                    y = 1:size(z,1);
                    plotFlag = true;
                case 2
                    z = varargin{numIdx(1)};
                    CData = varargin{numIdx(2)};
                    x = 1:size(z,2);
                    y = 1:size(z,1);
                    plotFlag = true;
                case 3
                    x = varargin{numIdx(1)};
                    y = varargin{numIdx(2)};
                    z = varargin{numIdx(3)};
                    plotFlag = true;
                case 4
                    x = varargin{numIdx(1)};
                    y = varargin{numIdx(2)};
                    z = varargin{numIdx(3)};
                    CData = varargin{numIdx(4)};
                    plotFlag = true;
                otherwise
                    badInput = true;
            end

        case {'trisurf','trimesh'}
            % For trisurf/trimesh, valid inputs are 1 (TriRep object),
            % 4 (tr,x,y,z) or 5 (tr,x,y,z,CData)
            if find(triRepIdx)
                triRep = varargin{triRepIdx};
                plotFlag = true;
                % Force LOGZPLOT to use gca for trimesh because it does not
                % properly utilize an axes handle input when called with
                % TriRep input
                if ~isempty(inputHandle) && ...
                        strcmp(func2str(plotFun),'trimesh')
                    msgid = [mfilename ':IgnoreInputHandle'];
                    msgtext = ['trimesh does not properly support axes' ...
                        ' handle input.  Using gca() for plot output.'];
                    warning(msgid,msgtext)
                    inputHandle = gca;
                end
            else
                if length(numIdx) == 4
                    tr = varargin{numIdx(1)};
                    x = varargin{numIdx(2)};
                    y = varargin{numIdx(3)};
                    z = varargin{numIdx(4)};
                    plotFlag = true;
                elseif length(numIdx) == 5
                    tr = varargin{numIdx(1)};
                    x = varargin{numIdx(2)};
                    y = varargin{numIdx(3)};
                    z = varargin{numIdx(4)};
                    CData = varargin{numIdx(5)};
                    plotFlag = true;
                else
                    badInput = true;
                end
            end
    end

    % Exit with an error if the inputs didn't work out
    if badInput
        msgid = [mfilename ':InvalidInputs'];
        msgtext = ['Invalid input.  Type ' ...
            'help ' mfilename ' for correct syntax'];
        error(msgid,msgtext)
    end

end % Done parsing input arguments

%% Make/modify plot

% Get handles:
[hFig,hAx,hSurf,priorFlag] = getHandles(inputHandle,plotFlag);

if priorFlag
    % Surface is already scaled, so just call the colorbar subfunction,
    % then exit:
    if colorbarFlag
        hCbar = LogZPlotColorbar(hFig,hAx);
    end

    % Set return values and exit
    if nargout > 0
        h = hSurf;
    end
    if nargout > 1
        h2 = hCbar;
    end

    return
end

if ~plotFlag
    % Determine the plot type and get its CData

    switch get(hSurf,'Type')
        case 'surface'
            CDataName = 'CData';
            % To determine if it is a pcolor plot, check if the ZData is
            % all zeros
            pcolorFlag = all(all(get(hSurf,'ZData')==0));
        case 'patch'
            CDataName = 'FaceVertexCData';
            pcolorFlag = false;
        case 'image'
            CDataName = 'CData';
            pcolorFlag = true; % avoid setting ZScale
    end

    % Now get the CData to work on
    CData = get(hSurf,CDataName);

    % If the CData is RGB, we don't know what to do, so throw a warning
    if size(CData,3) == 3;
        msgid = [mfilename ':TrueColorCData'];
        msgtext = ['The CData of the surface object is truecolor. ' ...
            'The surface will not be transformed'];
        warning(msgid,msgtext)

        % Set return values and exit
        if nargout > 0
            h = [];
        end
        if nargout > 1
            h2 = [];
        end

        return
    end

else % plotFlag == true
    % Make a plot with the user-suppled data

    % If we didn't get any CData as an input, use the ZData
    if isempty(CData)
        CData = z;
    end

    % Handle different plot types
    switch func2str(plotFun)
        case {'surf','mesh'}
            CDataName = 'CData';
            pcolorFlag = false;
            hSurf = plotFun(hAx,x,y,z,CData);
        case 'pcolor'
            CDataName = 'CData';
            pcolorFlag = true;
            hSurf = plotFun(hAx,x,y,CData);
        case {'trisurf','trimesh'}
            CDataName = 'FaceVertexCData';
            pcolorFlag = false;
            if isempty(triRep)
                hSurf = plotFun(tr,x,y,z,CData(:),'Parent',hAx);
            elseif strcmp(func2str(plotFun),'trimesh')
                % Workaround for broken axes handle support in trimesh
                hSurf = plotFun(triRep);
            else
                hSurf = plotFun(triRep,'Parent',hAx);
            end
            CData = get(hSurf,CDataName);
        case 'image'
            CDataName = 'CData';
            pcolorFlag = true;
            if isvector(x), x2 = x; else x2 = x(1,:); end
            if isvector(y), y2 = y; else y2 = y(:,1); end
            hSurf = plotFun(x2,y2,CData,'Parent',hAx);
            set(hAx, 'YDir','normal')
    end

    % Clear x,y,z now that we're done with them.  We are clearing z
    % _before_ we modify CData.  This way, when CData <= z from the
    % assignment statement a few lines above, MATLAB won't need to create
    % a second copy of the data in memory.
    clear x y z tr triRep

end

% Get rid of (CData <= 0) and (CData == Inf)
CData(CData<=0) = NaN;
CData(CData==Inf) = NaN;

% Find min/max
minC = min(CData(:));
maxC = max(CData(:));

% Convert data to log space
CData = log10(CData);

% Normalize the log scaled CData to have the same range as the ZData so
% that the colorbar scale will be correct
CData = minC + (maxC-minC)*(CData-log10(minC))/(log10(maxC/minC));

% Now set the CData of the surface to the normalized CData
set(hSurf,CDataName,CData)

% Make sure we are in 'scaled' CDataMapping mode
set(hSurf,'CDataMapping','scaled')

% Set a flag so we'll know the surface has already been transformed
setappdata(hSurf,'LogZPlotTransformedCData',true)

% Set the Z scale to 'log' if it isn't a pcolor plot
if ~pcolorFlag
    set(hAx,'ZScale','log')
end

% Call the colorbar subfunction
hCbar = [];
if colorbarFlag
    hCbar = LogZPlotColorbar(hFig,hAx);
end

% Store the CData range and set up the listeners
if listenerEnable

    setappdata(hSurf,'LogZPlotCLim',[minC maxC])
    setappdata(hSurf,'LogZPlotCLimOriginal',[minC maxC]);
    try
        fh = @(src,event)CLimChange(src,event,hAx,hSurf,true,CDataName);

        hCLimModeListener = handle.listener(handle(hAx),...
            findprop(handle(hAx),'CLimMode'), 'PropertyPreSet',...
            {@CLimModeChange,hAx,fh});

        hCLimListener = handle.listener(handle(hAx), ...
            findprop(handle(hAx),'CLim'), 'PropertyPostSet', ...
            {@CLimChange,hAx,hSurf,false,CDataName});

        listenerHandles = [hCLimModeListener; hCLimListener];
        setappdata(hSurf,'LogZPlotListeners',listenerHandles)

    catch %#ok<CTCH>
    end

end

% Set return values and exit
if nargout > 0
    h = hSurf;
end
if nargout > 1
    h2 = hCbar;
end

end % End of function logzplot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CLimChange(src,evnt,hAx,hSurf,autoFlag,CDataName) %#ok<INUSL>
% Listener function to update hSurf.CData with the proper scaling whenever
% hAx.CLim changes.  Executed post-set.

% Get old and new CLim values
priorCLim = getappdata(hSurf,'LogZPlotCLim');
newCLim = get(hAx,'CLim');

% autoFlag is set when this function is called by the CLimMode pre-set
% listener, signifying that we are going back to 'auto' mode and need to
% scale the data to the original limits.  This handles the case where
% CLimMode is reset after having been set to manual.
if autoFlag
    newCLim = getappdata(hSurf,'LogZPlotCLimOriginal');
end

% If newCLim and priorCLim are the same, no need to do anything
if abs(priorCLim(1) - newCLim(1)) < 10*eps(priorCLim(1)) ...
   && abs(priorCLim(2) - newCLim(2)) < 10*eps(priorCLim(2))
    return
end

% Warn and don't scale data if CLim is <= 0
if any(newCLim<=0)
    msgid = [mfilename ':NegativeCAxisLimit'];
    msgtext = [mfilename ' only supports positive CLim values.  ' ...
               'Colorbar scale will not be accurate.'];
    warning(msgid,msgtext)
    return
end

% Chop up priorCLim
pL = priorCLim(1);
pH = priorCLim(2);

% Chop up newCLim
nL = newCLim(1);
nH = newCLim(2);

% There are two parts to the transformation: 1) Undo the previous
% normalization (restore CData to its original values) 2) Normalize the
% CData to the new CLim.  They are shown here separately, but computed
% in-place in one step by computing the constants K1 and K2 below.

% 1) Undo the previous normalization to priorCLim:
% CData = (CData-pL)*(log10(pH/pL)/(pH-pL)) + log10(pL);

% 2) Normalize CData to new CLim:
% CData = nL+(nH-nL)*(CData-log10(nL))/log10(nH/nL);

% Compute renormalization constants K1,K2
K1 = nL+(nH-nL)/(log10(nH/nL))*(log10(pL)-log10(nL)-pL*log10(pH/pL)/(pH-pL));
K2 = (nH-nL)*log10(pH/pL)/(log10(nH/nL)*(pH-pL));

% Update the surface with the new CData
set(hSurf,CDataName,K1+K2*get(hSurf,CDataName));

% Write the new CLim to appdata:
setappdata(hSurf,'LogZPlotCLim',newCLim);

end % End of function logzplot/CLimChange
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CLimModeChange(src,evnt,hAx,fhCLimChange)  %#ok<INUSL>
% Listener function to reset the CData using fhCLimChange whenever
% the CLimMode is changed back to 'auto'.  This needs to happen _before_
% the CLimMode is actually changed so that the auto-mode axes refresh uses
% the correct (new) limits instead of basing them on the old CData.

if strcmp(get(hAx,'CLimMode'),'manual')
    fhCLimChange(0,0)
end

end % End of function logzplot/CLimModeChange
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hFig,hAx,hSurf,priorFlag] = ...
                getHandles(inputHandle,plotFlag)
% Given an input handle and plotFlag, get figure, axes, surface handles and
% priorFlag.

% Default return value unless we find a flag on the surface object
priorFlag = false;

% Return empty for hSurf if plotFlag == true
hSurf = [];

if ~isempty(inputHandle)
    % If we were given a handle, determine its type and then get the other
    % handles we need
    switch lower(get(inputHandle,'Type'))
        case 'axes'
            hAx = inputHandle;
            hFig = ancestor(hAx,'figure');
        case {'surface','patch','image'}
            hSurf = inputHandle;
            hAx = ancestor(hSurf,'axes');
            hFig = ancestor(hSurf,'figure');
        otherwise
            msgid = [mfilename ':InvalidInputHandle'];
            msgtext = ['Input handle must refer to a valid axes or '...
                'surface object'];
            error(msgid,msgtext);
    end
else
    % No input handle so use the current figure and axes
    hFig = gcf;
    hAx = gca;
end

if ~plotFlag
    % If necessary, try to find a surface object
    if isempty(hSurf)
        hSurf = findobj(hAx,'Type','surface','-or','Type','patch',...
            '-or','Type','image');
    end

    % If we couldn't find one, or found more than one, exit
    if isempty(hSurf)
        msgid = [mfilename ':NoSurfaceObject'];
        msgtext = 'Could not locate a surface object to transform';
        error(msgid,msgtext);
    elseif length(hSurf) > 1
        msgid = [mfilename ':MultipleSurfaceObjects'];
        msgtext = [mfilename ' only supports one surface object per axes'];
        error(msgid,msgtext)
    end

    % Check for prior logzplot transformation
    priorFlag = getappdata(hSurf,'LogZPlotTransformedCData');

end

end % End of function logzplot/getHandles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hCbar = LogZPlotColorbar(hFig,hAx,hCbar)

if nargin == 2
    % See if there are existing colorbars in the figure
    hCbar = findobj(hFig,'tag','Colorbar');
end

% We only want to modify a colorbar whose peer axes is hAx.  To check this,
% use the not fully documented handle() and h.axes property.  Best to wrap
% this in a try block and do nothing if it fails.
try
    for k = length(hCbar):-1:1
        hTestAx = handle(hCbar(k));
        if (double(hTestAx.axes) ~= hAx)
            hCbar(k) = [];
        end
    end
catch   %#ok<CTCH>
    hCbar = [];
end

% For multiple colorbars, call this function recursively for each one
if length(hCbar) > 1
    for k = 1:length(hCbar)
        LogZPlotColorbar(hFig,hAx,hCbar(k));
    end
    return
end

% If there isn't a colorbar, make it now
if isempty(hCbar)
    hCbar = colorbar('peer',hAx);
end

switch get(hCbar,'Location')
    case {'East','West','EastOutside','WestOutside'}
        scaleName = 'YScale';
    case {'North','South','NorthOutside','SouthOutside'}
        scaleName = 'XScale';
end

% Set the colorbar axes scale to 'log'
set(hCbar,scaleName,'log')

end % End of function logzplot/LogZPlotColorbar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function [Zpca, U, mu, eigVecs] = PCA(Z,r)
%
% Syntax:       Zpca = PCA(Z,r);
%               [Zpca, U, mu] = PCA(Z,r);
%               [Zpca, U, mu, eigVecs] = PCA(Z,r);
%
% Inputs:       Z is an d x n matrix containing n samples of d-dimensional
%               data
%
%               r is the number of principal components to compute
%
% Outputs:      Zpca is an r x n matrix containing the r principal
%               components - scaled to variance 1 - of the input samples
%
%               U is a d x r matrix of coefficients such that
%               Zr = U * Zpca + repmat(mu,1,n);
%               is the r-dimensional PCA approximation of Z
%
%               mu is the d x 1 sample mean of Z
%
%               eigVecs is a d x r matrix containing the scaled
%               eigenvectors of the sample covariance of Z
%
% Description:  Performs principal component analysis (PCA) on the input
%               data
%
% Author:       Brian Moore
%               brimoor@umich.edu
%
% Date:         April 26, 2015
%               November 7, 2016
%

% Center data
[Zc, mu] = centerRows(Z);

% Compute truncated SVD
%[U, S, V] = svds(Zc,r); % Equivalent, but usually slower than svd()
[U, S, V] = svd(Zc,'econ');
U = U(:,1:r);
S = S(1:r,1:r);
V = V(:,1:r);

% Compute principal components
Zpca = S * V';
%Zpca = U' * Zc; % Equivalent but slower

if nargout >= 4
    % Scaled eigenvectors
    eigVecs = bsxfun(@times,U,diag(S)' / sqrt(size(Z,2)));
end

author: xiuyechen

repo: FishExplorer


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment2B
close all;
global robot
global parameters
global hfigures

%STOMP PARAMETERS
%conversion from cost to Prob factor
parameters.lambda = .4;
parameters.lambda_obstacles = .2;
%height of the obstacle
parameters.yo = 2.5;
%cost function starts at this distance
%must be below 0.3 for the 4 DOF robot
parameters.epsilon = 0.2;
%multiply noise by this facto
%parameters.noise_k = 5;
%parameters.noise_sigma_null_space = 0.01;
parameters.alpha=0.02;
parameters.time_step=0.01;

%number of waypoints
parameters.N = 12;
%number of particles
parameters.K = 20;
parameters.n_repeat = 30;
parameters.experiment_name = 'experiment2B_K20_N30.mat';

parameters.obstacles = [];

parameters.animate = 0;
close all
hfigures.hpaths = figure;
hfigures.hcosts = figure;
hfigures.hee = figure;
hfigures.htheta = figure;
hfigures.hbest_costs = figure;
hfigures.hdtheta = figure;

parameters.obstacles = [];
%LINE 1
x1 = -1.5;
y1 = .5; %m
x2 = 0;
y2 = 1.5; %m
phi = 3*pi/4;
p0 = [x1 y1 0]';
pf = [x2 y2 0]';
T0 = build_T_4dof(p0, phi);
parameters.obstacles{1}.line = [p0 pf];
parameters.obstacles{1}.T0 = T0;

%LINE 2
x1 = 0;
y1 = 1.5; %m
x2 = 1;
y2 = 1.5; %m
phi = pi/2;
p0 = [x1 y1 0]';
pf = [x2 y2 0]';
T0 = build_T_4dof(p0, phi);
parameters.obstacles{2}.line = [p0 pf];
parameters.obstacles{2}.T0 = T0;

%LINE 3
x1 = 1;
y1 = 1.5; %m
x2 = 2;
y2 = 0.5; %m
phi = pi/4;
p0 = [x1 y1 0]';
pf = [x2 y2 0]';
T0 = build_T_4dof(p0, phi);
parameters.obstacles{3}.line = [p0 pf];
parameters.obstacles{3}.T0 = T0;

%repeat the experiment E times
random_manips=[];
Gout = [];
for i=1:parameters.n_repeat
    close all
    [pk, final_manip] = stomp_null_space(robot);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
end

function T = build_T_4dof(p, phi)
T = [cos(phi) -sin(phi) 0 p(1);
     sin(phi) cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

 function T = build_T_sawyer(p, phi)
T = [1  -sin(phi) 0 p(1);
     0  cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

repo: ARTE


text:
function [Ef, Varf, lpy, Ey, Vary] = gpmc_loopreds(gp, x, y, varargin)
%GPMC_LOOPREDS  Leave-one-out predictions with Gaussian Process MCMC approximation
%
%  Description
%    [EFT, VARFT, LPYT, EYT, VARYT] = GPMC_LOOPREDS(RECGP, X, Y)
%    takes a Gaussian processes record structure RECGP (returned by
%    gp_mc) together with a matrix X of training inputs and vector
%    Y of training targets, and evaluates the leave-one-out
%    predictive distribution at inputs X and returns means EFT and
%    variances VARFT of latent variables, the logarithm of the
%    predictive densities PYT, and the predictive means EYT and
%    variances VARYT of observations at input locations X.
%
%    OPTIONS is optional parameter-value pair
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%    Note:
%      - the hyperparameters, hp, are sampled from the full
%        posterior p(hp|x,y) by gp_mc
%      - in case of Gaussian likelihood or non-Gaussian likelihood
%        and latent method Laplace/EP, the conditonal LOO-CV
%        distributions p(f_i | x_\i, y_\i, z_\i, hp_s) are computed
%        for each hyperparameter sample hp_s
%      - use gp_loopred to evaluate p(f_i | x_\i, y_\i, z_\i)
%
%  See also
%   GP_LOOPRED
%
% Copyright (c) 2008-2010,2012 Jarno Vanhatalo
% Copyright (c) 2010,2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Nothing to parse, but check the arguments anyway
ip=inputParser;
ip.FunctionName = 'GPMC_LOOPREDS';
ip.addRequired('gp',@isstruct);
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.parse(gp, x, y, varargin{:});
z=ip.Results.z;

if isfield(gp,'meanf') & ~isempty(gp.meanf)
  error('GPMC_LOOPREDS: Mean functions not yet supported');
end

nmc=size(gp.jitterSigma2,1);

if isfield(gp, 'latentValues') && ~isempty(gp.latentValues)
  % Non-Gaussian likelihood. The latent variables should be used in
  % place of observations
  lv = gp.latentValues';
end

for i1=1:nmc
  Gp = take_nth(gp,i1);
  if isfield(Gp,'latent_method') && isequal(Gp.latent_method,'MCMC')
    Gp = rmfield(Gp,'latent_method');
  end

  if isfield(gp, 'latentValues') && ~isempty(gp.latentValues)
    % latent values have been sampled with MCMC
    [Ef(:,i1), Varf(:,i1)] = gp_loopred(Gp, x, lv(:,i1), varargin{:});
    if nargout >=4
      [lpy(:,i1), Ey(:,i1), Vary(:,i1)] = Gp.lik.fh.predy(Gp.lik, Ef(:,i1), Varf(:,i1), y, z);
    else
      lpy(:,i1) = Gp.lik.fh.predy(Gp.lik, Ef(:,i1), Varf(:,i1), y, z);
    end
  else
    % Gaussian likelihood or Laplace/EP for latent values
    if nargout <= 2
      [Ef(:,i1), Varf(:,i1)] = gp_loopred(Gp, x, y, varargin{:});
    elseif nargout <=3
      [Ef(:,i1), Varf(:,i1), lpy(:,1)] = gp_loopred(Gp, x, y, varargin{:});
    else
      [Ef(:,i1), Varf(:,i1), lpy(:,1), Ey(:,i1), Vary(:,i1)] = gp_loopred(Gp, x, y, varargin{:});
    end
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function determ = cheby_diff1_determinant ( n )

%*****************************************************************************80
%
%% CHEBY_DIFF1_DETERMINANT returns the determinant of the CHEBY_DIFF1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 0.0;

  return
end

repo: jburkardt-m


text:
function circle_imp_print_2d ( r, center, title )

%*****************************************************************************80
%
%% CIRCLE_IMP_PRINT_2D prints an implicit circle in 2D.
%
%  Discussion:
%
%    An implicit circle in 2D satisfies:
%
%      ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 = R^2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Radius = %f\n', r );
  fprintf ( 1, '  Center = %f  %f\n', center(1:2,1) );

  return
end