text:
function [ component_num, c ] = i4block_components ( l, m, n, a )

%*****************************************************************************80
%
%% I4BLOCK_COMPONENTS assigns contiguous nonzero pixels to a common component.
%
%  Discussion:
%
%    On input, the A array contains values of 0 or 1.
%
%    The 0 pixels are to be ignored.  The 1 pixels are to be grouped
%    into connected components.
%
%    The pixel A(I,J,K) is "connected" to the pixels:
%
%      A(I-1,J,  K  ),  A(I+1,J,  K  ),
%      A(I,  J-1,K  ),  A(I,  J+1,K  ),
%      A(I,  J,  K-1),  A(I,  J,  K+1),
%
%    so most pixels have 6 neighbors.
%
%    On output, COMPONENT_NUM reports the number of components of nonzero
%    data, and the array C contains the component assignment for
%    each nonzero pixel, and is 0 for zero pixels.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer L, M, N, the order of the array.
%
%    Input, integer A(L,M,N), the pixel array.
%
%    Output, integer COMPONENT_NUM, the number of components
%    of nonzero data.
%
%    Output, integer C(L,M,N), the component array.
%

%
%  Initialization.
%
  c = zeros ( l, m, n );
  component_num = 0;
%
%  P is simply used to store the component labels.  The dimension used
%  here is, of course, usually an absurd overestimate.
%
  p = 1 : l * m * n;
%
%  "Read" the array one pixel at a time.  If a (nonzero) pixel's north or
%  west neighbor already has a label, the current pixel inherits it.
%  In case the labels disagree, we need to adjust the P array so we can
%  later deal with the fact that the two labels need to be merged.
%
  for i = 1 : l

    for j = 1 : m

      for k = 1 : n

        if ( i == 1 )
          north = 0;
        else
          north = c(i-1,j,k);
        end

        if ( j == 1 )
          west = 0;
        else
          west = c(i,j-1,k);
        end

        if ( k == 1 )
          up = 0;
        else
          up = c(i,j,k-1);
        end

        if ( a(i,j,k) ~= 0 )
%
%  New component?
%
          if ( north == 0 && west == 0 && up == 0 )
            component_num = component_num + 1;
            c(i,j,k) = component_num;
%
%  One predecessor is labeled.
%
          elseif ( north ~= 0 && west == 0 && up == 0 )
            c(i,j,k) = north;
          elseif ( north == 0 && west ~= 0 && up == 0 )
            c(i,j,k) = west;
          elseif ( north == 0 && west == 0 && up ~= 0 )
            c(i,j,k) = up;
%
%  Two predecessors are labeled.
%
          elseif ( north == 0 && west ~= 0 && up ~= 0 )
            c(i,j,k) = min ( west, up );
            c1 = min ( p(west), p(up) );
            p(west) = c1;
            p(up) = c1;
          elseif ( north ~= 0 && west == 0 && up ~= 0 )
            c(i,j,k) = min ( north, up );
            c1 = min ( p(north), p(up) );
            p(north) = c1;
            p(up) = c1;
          elseif ( north ~= 0 && west ~= 0 && up == 0 )
            c(i,j,k) = min ( north, west );
            c1 = min ( p(north), p(west) );
            p(north) = c1;
            p(west) = c1;
%
%  Three predecessors are labeled.
%
          elseif ( north ~= 0 && west ~= 0 && up ~= 0 )
            c(i,j,k) = min ( north, min ( west, up ) );
            c1 = min ( p(north), min ( p(west), p(up) ) );
            p(north) = c1;
            p(west) = c1;
            p(up) = c1;
          end

        end

      end

    end

  end
%
%  When a component has multiple labels, have the higher labels
%  point to the lowest one.
%
  for component = component_num : -1 : 1
    b = component;
    while ( p(b) ~= b )
      b = p(b);
    end
    p(component) = b;
  end
%
%  Locate the minimum label for each component.
%  Assign these mininum labels new consecutive indices.
%
  q = zeros ( 1, component_num );
  i = 0;
  for component = 1 : component_num
    if ( p(component) == component )
      i = i + 1;
      q(component) = i;
    end
  end

  component_num = i;
%
%  Replace the labels by consecutive labels.
%
  i = find ( c ~= 0 );
  c(i) = q ( p ( c(i) ) );

  return
end


text:
function test_bug1262

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_header ft_read_data ft_read_event

% Quick and dirty sanity check for reading of fcdc_buffer_offline data.

% Example data generated with sine example, saved with record.exe:
dirname = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1262/0001');

hdr = ft_read_header(dirname);
assert(hdr.nChans == 16);
assert(hdr.nSamples == 1216);

dat = ft_read_data(dirname);
assert(all(size(dat) == [hdr.nChans, hdr.nSamples]));

evt = ft_read_event(dirname);
assert(numel(evt) == 0);  % there are no events

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef (Abstract) mm_shared_pfcpf_acpi < mp.mm_shared_pfcpf_acp & mp.mm_shared_pfcpf_ac_i

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function ad = build_aux_data(obj, nm, dm, mpopt)
            %% call parent
            ad = build_aux_data@mp.mm_shared_pfcpf_acp(obj, nm, dm, mpopt);

            %% add data needed for current formulations
            ad = obj.build_aux_data_i(nm, ad);
        end

        function obj = add_system_vars_pf(obj, nm, dm, mpopt)
            %% get model variables
            vvars = nm.model_vvars();

            %% voltage angles
            obj.add_system_varset_pf(nm, vvars{1}, 'pv');
            obj.add_system_varset_pf(nm, vvars{1}, 'pq');

            %% reactive injections
            ad = obj.aux_data;
            v_ = ad.vm .* exp(1j * ad.va);
            z_ = ad.zr + 1j * ad.zi;
            Qpv = nm.C(ad.pv, :) * imag( nm.port_inj_power([v_; z_], 1) );
            Qg_pv = Qpv + ad.zi(ad.zi_idx);
            mmx_i1 = obj.var.N + 1;
            obj.add_var('Qg_pv', ad.npv, Qg_pv);
            mmx_iN = obj.var.N;
            if ad.npv
                obj.aux_data.var_map{end+1} = ...
                    {'zi', [], [], ad.zi_idx, mmx_i1, mmx_iN, []};
            end

            %% voltage magnitudes
            obj.add_system_varset_pf(nm, vvars{2}, 'pq');
        end

        function [f, J] = node_balance_equations(obj, x, nm)
            %% index vector
            ad = obj.aux_data;
            pvq = [ad.pv; ad.pq];

            %% update network model state ([v_; z_]) from math model state (x)
            [v_, z_] = obj.convert_x_m2n(x, nm, 1);

            %% incidence matrix
            C = nm.C;

            %% Jacobian
            if nargout > 1
                %% get port current injections with derivatives
                [I, dI.va, dI.vm, dI.zr, dI.zi] = nm.port_inj_current([v_; z_], 1);
                dI.va = C * dI.va;
                dI.vm = C * dI.vm;
                dI.zr = C * dI.zr;
                dI.zi = C * dI.zi;
                JJ = cell(2, length(ad.var_map));

                for k = 1:length(ad.var_map)
                    m = ad.var_map{k};
                    name = m{1};
                    if ~isempty(m{2})       %% i1:iN
                        i1 = m{2};
                        iN = m{3};
                        JJ{1, k} = real(dI.(name)(pvq, i1:iN));
                        JJ{2, k} = imag(dI.(name)(pvq, i1:iN));
                    elseif isempty(m{4})    %% :
                        JJ{1, k} = real(dI.(name)(pvq, :));
                        JJ{2, k} = imag(dI.(name)(pvq, :));
                    else                    %% idx
                        idx = m{4};
                        JJ{1, k} = real(dI.(name)(pvq, idx));
                        JJ{2, k} = imag(dI.(name)(pvq, idx));
                    end
                end
                J = vertcat( horzcat(JJ{1, :}), ...
                             horzcat(JJ{2, :})  );
            else
                %% get port current injections (w/o derivatives)
                I = nm.port_inj_current([v_; z_], 1);
            end

            %% nodal power balance
            II = C * I;
            f = [real(II(pvq)); imag(II(pvq))];
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [ x, out, opts ] = tfocs_TS( smoothF, affineF, projectorF, x0, opts )
% TFOCS_TS Tseng's modification of Nesterov's 2007 method.
% [ x, out, opts ] = tfocs_TS( smoothF, affineF, projectorF, x0, opts )
%   Implements Tseng's modification of the Nesterov 2007 method.
%   A variety of calling sequences are supported; type
%      help tfocs_help
%   for a full explanation.

% Nov 17 2016, hack for Matlab R2016b not allowing nargin/nargout in
% tfocs_initialize:
narginn = nargin; nargoutt = nargout;

% Initialization
alg = 'TS';
algorithm = 'Tseng''s single-projection modification of Nesterov''s 2007 method';
alpha = 0; beta = 0; mu = 0; L = 0; % Necessary due to MATLAB quirk
tfocs_initialize
if nargin == 0, return; end

while true,

      x_old =   x;
    A_x_old = A_x;
      z_old =   z;
    A_z_old = A_z;

    % The backtracking loop
    L_old      = L;
   	L          = L * alpha;
    theta_old  = theta;
    while true,

    	% Acceleration
        theta = 2 ./ ( 1 + sqrt( 1 + 4 * L / L_old / theta_old^2 ) );

        % Next iterate
        if theta < 1,
              y = ( 1 - theta ) *   x_old + theta *   z_old;
            A_y = ( 1 - theta ) * A_x_old + theta * A_z_old;
            f_y = Inf; g_Ay = []; g_y = []; C_y = Inf;
        end

        % Compute function values
        if isempty( g_y ),
            if isempty( g_Ay ), [ f_y, g_Ay ] = apply_smooth( A_y ); end
            g_y = apply_linear( g_Ay, 2 );
        end

        % Accumulated gradient
        if theta == 1,
            x_cent = x_old;
            g_a = g_y;
        else
            g_a = ( 1 - theta ) * g_a_old + theta * g_y;
        end
        step = 1 / ( theta^2 * L );
        [ C_z, z ] = apply_projector( x_cent - step * g_a, step );
        A_z = apply_linear( z, 1 );

        % New iterate
        if theta == 1,
            x   = z;
            A_x = A_z;
            C_x = C_z;
        else
            x   = ( 1 - theta ) *   x_old + theta *   z;
            A_x = ( 1 - theta ) * A_x_old + theta * A_z;
            C_x = Inf;
        end
        f_x = Inf; g_Ax = []; g_x = [];

        % Bactracking test
        tfocs_backtrack
        if do_break, break; end % new, for R2015b compatibility

    end

    % Collect data, evaluate stopping criteria, and print status
    tfocs_iterate
    if do_break, break; end % new, for R2015b compatibility

    g_a_old = g_a;

end

% Final processing
tfocs_cleanup

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
%The spherical K-means algorithm
%(C) 2007-2011 Nguyen Xuan Vinh
%Contact: vinh.nguyenx at gmail.com   or vinh.nguyen at monash.edu
%Reference:
%   [1] Xuan Vinh Nguyen: Gene Clustering on the Unit Hypersphere with the
%       Spherical K-Means Algorithm: Coping with Extremely Large Number of Local Optima. BIOCOMP 2008: 226-233
%Usage: Normalize the data set to have unit norm

function b=normalize_norm(a)
[n dim]=size(a);
for i=1:n
    a(i,:)=a(i,:)/norm(a(i,:));
end
b=a;

author: Sable

repo: mcbench-benchmarks


text:
function f = eightPoint(points1homo, points2homo)
% Normalize the points
num = size(points1homo, 2);
[points1homo, t1] = vision.internal.normalizePoints(points1homo, 2, 'double');
[points2homo, t2] = vision.internal.normalizePoints(points2homo, 2, 'double');
% unravel
m = coder.nullcopy(zeros(num, 9, 'double'));
m(:,1)=(points1homo(1,:).*points2homo(1,:))';
m(:,2)=(points1homo(2,:).*points2homo(1,:))';
m(:,3)=points2homo(1,:)';
m(:,4)=(points1homo(1,:).*points2homo(2,:))';
m(:,5)=(points1homo(2,:).*points2homo(2,:))';
m(:,6)=points2homo(2,:)';
m(:,7)=points1homo(1,:)';
m(:,8)=points1homo(2,:)';
m(:,9)=1;
% last eigen vector
[~, ~, vm] = svd(m, 0);
f = reshape(vm(:, end), 3, 3)';
[u, s, v] = svd(f);
s(end) = 0;
f = u * s * v';
% denormalize
f = t2' * f * t1;
f = f / norm(f);
if f(end) < 0
  f = -f;
end

author: yihui-he


text:
function r = hammersley_sequence ( n )

%*****************************************************************************80
%
%% HAMMERSLEY_SEQUENCE computes N elements of a leaped Hammersley subsequence.
%
%  Discussion:
%
%    The DIM_NUM-dimensional Hammersley sequence is really DIM_NUM separate
%    sequences, each generated by a particular base.  If the base is
%    greater than 1, a standard 1-dimensional
%    van der Corput sequence is generated.  But if the base is
%    negative, this is a signal that the much simpler sequence J/(-BASE)
%    is to be generated.  For the standard Hammersley sequence, the
%    first spatial coordinate uses a base of (-N), and subsequent
%    coordinates use bases of successive primes (2, 3, 5, 7, 11, ...).
%    This program allows the user to specify any combination of bases,
%    included nonprimes and repeated values.
%
%    This routine selects elements of a "leaped" subsequence of the
%    Hammersley sequence.  The subsequence elements are indexed by a
%    quantity called STEP, which starts at 0.  The STEP-th subsequence
%    element is simply element
%
%      SEED(1:DIM_NUM) + STEP * LEAP(1:DIM_NUM)
%
%    of the original Hammersley sequence.
%
%
%    This routine "hides" a number of input arguments.  To specify these
%    arguments explicitly, use I4_TO_HAMMERSLEY instead.
%
%    All the arguments have default values.  However, if you want to
%    examine or change them, you may call the appropriate routine first.
%
%    * DIM_NUM, the spatial dimension,
%      Default: DIM_NUM = 1;
%      Required: 1 <= DIM_NUM is required.
%
%    * STEP, the subsequence index.
%      Default: STEP = 0.
%      Required: 0 <= STEP.
%
%    * SEED(1:DIM_NUM), the Hammersley sequence element corresponding to STEP = 0.
%      Default SEED = (0, 0, ... 0).
%      Required: 0 <= SEED(1:DIM_NUM).
%
%    * LEAP(1:DIM_NUM), the succesive jumps in the Hammersley sequence.
%      Default: LEAP = (1, 1, ..., 1).
%      Required: 1 <= LEAP(1:DIM_NUM).
%
%    * BASE(1:DIM_NUM), the Hammersley bases.
%      Default: BASE = (2, 3, 5, 7, 11, ... ) or (-N,2,3,5,7,11,...) if N is known.
%      Required: 1 < BASE(I) for any van der Corput dimension, or BASE(I) < 0
%      to generate the fractional sequence J/|BASE(I)|.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    J M Hammersley,
%    Monte Carlo methods for solving multivariable problems,
%    Proceedings of the New York Academy of Science,
%    Volume 86, 1960, pages 844-874.
%
%    Ladislav Kocis and William Whiten,
%    Computational Investigations of Low-Discrepancy Sequences,
%    ACM Transactions on Mathematical Software,
%    Volume 23, Number 2, 1997, pages 266-294.
%
%  Parameters:
%
%    Input, integer N, the number of elements to compute.
%
%    Output, real R(DIM_NUM,N), the SEED-th through SEED+N-1th elements of
%    the Hammersley sequence.
%
  dim_num = hammersley_dim_num_get ( );
  step = hammersley_step_get ( );
  seed = hammersley_seed_get ( );
  leap = hammersley_leap_get ( );
  base = hammersley_base_get ( );

  r = i4_to_hammersley_sequence ( dim_num, n, step, seed, leap, base );

  step = step + n;

  hammersley_step_set ( step );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Partial fraction expansion of a rational function
%
%

%        s^2+3s+1
% X(s)=----------------
%       s^3+5s^2+2s-8

A=[1 5 2 -8];
rt=roots(A)
syms s
X=(s^2+3*s+1)/(s^3+5*s^2+2*s-8);

c1=limit((s-rt(1))*X,s,rt(1))
c2=limit( (s-rt(2))*X,s,rt(2))
c3=limit( (s-rt(3))*X,s,rt(3))

X1=simplify( (s-rt(1))*X );
X2=simplify( (s-rt(2))*X );
X3=simplify( (s-rt(3))*X );
c1=subs(X1,s,rt(1))
c2=subs(X2,s,rt(2))
c3=subs(X3,s,rt(3))

%alternative conversion
num=[ 1 3 1];
den=[ 1 5 2 -8];
[R,P,K]=residue(num,den)
syms s
X=R(1)/(s-P(1))+R(2)/(s-P(2))+R(3)/(s-P(3))
pretty(X)

%
%        s^2+3s+1
% X(s)=------------
%       s^3-3s+2

A=[1 0 -3 2];
rt=roots(A)
syms s
X=(s^2+3*s+1)/(s^3-3*s+2);

c1=limit((s-rt(1))*X,s,rt(1))

r=2;
i=1;
f= ((s-1)^r )*X;
d=diff(f,s,r-i);
fact=1/factorial(r-i);
c2=limit(fact*d,s,1)

i=2;
d=diff(f,s,r-i);
fact=1/factorial(r-i);
c3=limit(fact*d,s,1)

%alternative conversion
num=[ 1 3 1];
den=[1 0 -3 2];
[R,P,K]=residue(num,den)

%
%        s^2-12s+11
% X(s)=------------
%        s^2+4s+3

num=[1 -12 11];
den=[ 1 4 3];
[k,g]=deconv(num,den)
rt=roots(den)

syms s
GA=(-16*s+8)/(s^2+4*s+3);
c1=limit((s-rt(1))*GA,s,rt(1))
c2=limit((s-rt(2))*GA,s,rt(2))

%alternative conversion
num=[ 1 -12 11];
den=[ 1 4 3];
[R,P,K]=residue(num,den)


text:
function h_ax = plotSphFunctionGrid(F, aziRes, polarRes, realComplex, h_ax)
%PLOTSPHFUNCTIONGRID Plots a spherical function defined on a grid
%
%   F:  matrix of function values on the grid points
%   aziRes: grid resolution at azimuth (degrees)
%   polarRes: grid resolution in elevation (degrees)
%   realComplex: {'real','complex'} if the function is real then it is
%                plotted with one surface for its positive part and one for
%                its negative part. If it is complex, the magnitude
%                function is plotted, with its phase mapped on the colormap
%   h_ax: optional argument to define an axis handle for the plot,
%         otherwise the new axis handle is returned
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Archontis Politis, 20/02/2015
%   archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==4
    figure
    h_ax = axes;
elseif nargin==3
    figure
    h_ax = axes;
    realComplex = 'complex';
else
    axes(h_ax);
end

% construct grid
azi = 0:aziRes:360;
elev = 0:polarRes:180;
azi_rad = azi*pi/180;
elev_rad = elev*pi/180;
[Az, El] = meshgrid(azi_rad, elev_rad);

% construct real positive and negative parts if real function
D_x = cos(Az).*sin(El).*squeeze(abs(F));
D_y = sin(Az).*sin(El).*squeeze(abs(F));
D_z = cos(El).*squeeze(abs(F));

if isequal(realComplex, 'real')
    Dp_x = D_x.*(F>=0);
    Dp_y = D_y.*(F>=0);
    Dp_z = D_z.*(F>=0);
    Dn_x = D_x.*(F<0);
    Dn_y = D_y.*(F<0);
    Dn_z = D_z.*(F<0);
elseif isequal(realComplex, 'complex')
    Dm_x = D_x.* abs(F);
    Dm_y = D_y.* abs(F);
    Dm_z = D_z.* abs(F);
end

% plot 3d axes
maxF = max(max(abs(F)));
line([0 1.5*maxF],[0 0],[0 0],'color',[1 0 0])
line([0 0],[0 1.5*maxF],[0 0],'color',[0 1 0])
line([0 0],[0 0],[0 1.5*maxF],'color',[0 0 1])

% plot function
hold on
if isequal(realComplex, 'real')
    Hp = surf(Dp_x, Dp_y, Dp_z);
    Hn = surf(Dn_x, Dn_y, Dn_z);
    set(Hp, 'FaceColor', 'b')
    set(Hp, 'EdgeAlpha', 1)
    set(Hn, 'FaceColor', 'r')
    set(Hn, 'EdgeAlpha', 1)
elseif isequal(realComplex, 'complex')
    Hm = surf(Dm_x, Dm_y, Dm_z, angle(F));
    set(Hm, 'EdgeAlpha', 1)
end
xlabel('x')
ylabel('y')
zlabel('z')
light('Position',[0 0 1],'Style','infinite');
light('Position',[-1 -1 -1],'Style','infinite');
light('Position',[0 0 -1],'Style','infinite');
material shiny
axis equal
grid

end

author: polarch

repo: Spherical-Harmonic-Transform


text:
function F = besselk(nu, F, scale, pref)
%BESSELK   Modified Bessel function of second kind of a CHEBFUN.
%   K = BESSELK(NU, F) computes the modified Bessel function of second kind
%   K_NU(F) of the nonzero CHEBFUN F. If F passes through the origin in its
%   domain, then an error is returned. The order NU need not be an integer but
%   must be real. The argument F can be complex.
%
%   K = BESSELK(NU, F, SCALE) returns a scaled K_NU(F) specified by SCALE:
%         0 - (default) is the same as BESSELK(NU, F),
%         1 - scales K_NU(F) by exp(F)).
%
%   K = BESSELK(NU, F, SCALE, PREF) uses the CHEBFUNPREF object PREF when
%   building the CHEBFUN K.
%
% See also AIRY, BESSELH, BESSELI, BESSELJ, BESSELY.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( nargin < 4 )
    pref = chebfunpref();
end
if ( nargin < 3 )
    scale = 0;
end

% Loop over the columns:
for k = 1:numel(F)
    F(k) = columnBesselk(nu, F(k), scale, pref);
end

end

function g = columnBesselk(nu, f, scale, pref)

% Check for roots:
r = roots(f, 'nojump', 'nozerofun');
if ( numel(r) > 0 )
    error('CHEBFUN:CHEBFUN:besselk:zero', 'F has roots in its domain.');
end

% Compose:
g = compose(f, @(x) besselk(nu, x), pref);

% Scale (as described in help documentation):
if ( scale == 1 )
    scl = exp(f);
    g = scl.*g;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function c = acos(a)
% ACOS for adiff objects.

c = adiff(acos(a.x),rowmult(-1./sqrt(1-a.x.^2), a.dx),a.root);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Globally Optimal Solutions,
%           Solvability/Observability Analysis and Uncertainty Description.
%           IEEE Transactions on Robotics.
%           https://doi.org/10.1109/TRO.2022.3155880
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com
%
%
% test_stewart.m: The QPEP illustration of forwart kinematics of hexapod
%                 Stewart platform

clear all
close all
clc

if(verLessThan('matlab', '8.0.0'))
   error('The MATLAB version is too old to be supported.');
end

format long g

addpath('func_files');
addpath('solvers');
addpath('utils');
addpath('homotopy');

R0 = angle2dcm(-8 * pi / 180, 12 * pi / 180, -15 * pi / 180, 'XYZ');
q0 = dcm2quat(R0).';
if(q0(1) < 0)
    q0 = - q0;
end
t0 = 1e-2 * randn(3, 1);
X0 = inv([R0, t0;
     zeros(1, 3), 1]);

base = [
     0.1448888739433600, 1,  0.0388228567653781;
    -0.0388228567653781, 1,  0.1448888739433600;
    -0.1060660171779820, 1,  0.1060660171779820;
    -0.1060660171779820, 1, -0.1060660171779820;
    -0.0388228567653781, 1, -0.1448888739433600;
     0.1448888739433600, 1, -0.0388228567653781;
    ].';
plat = [
     0.0707106781186548, 1, 0.07071067811865480;
     0.0258819045102521, 1, 0.09659258262890680;
    -0.0965925826289068, 1, 0.02588190451025210;
    -0.0965925826289068, 1, -0.0258819045102521;
     0.0258819045102521, 1, -0.0965925826289068;
     0.0707106781186548, 1, -0.0707106781186548;
    ].';

conv = [
    1, 0, 0;
    0, 0, 1;
    0, -1, 0;
    ];
height = 0.15;
base = conv * base;
plat = conv * plat;
plat(3, :) = plat(3, :) + height;
plat00 = plat;
base00 = base;

leg0 = zeros(6, 1);
plat0 = zeros(3, 6);
for i = 1 : 6
    plat0(:, i) = R0 * plat(:, i)  + t0;
    res = base(:, i) - plat0(:, i);
    leg0(i) = sqrt(res.' * res);
end
colors = linspecer(8);

figure(1);
subplot(1, 2, 1);
plot3(base(1, :), base(2, :), base(3, :), 'LineStyle', 'None', 'Marker', '.', 'MarkerSize', 10); hold on
plot3(base(1, 1 : 6), base(2, 1 : 6), base(3, 1 : 6), 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
plot3([base(1, 1), base(1, 6)], [base(2, 1), base(2, 6)], [base(3, 1), base(3, 6)], 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
plot3(plat0(1, :), plat0(2, :), plat0(3, :), 'LineStyle', 'None', 'Marker', '.', 'MarkerSize', 10); hold on
plot3(plat0(1, 1 : 6), plat0(2, 1 : 6), plat0(3, 1 : 6), 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
plot3([plat0(1, 1), plat0(1, 6)], [plat0(2, 1), plat0(2, 6)], [plat0(3, 1), plat0(3, 6)], 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
for i = 1 : 6
    plot3([base(1, i), plat0(1, i)], [base(2, i), plat0(2, i)], [base(3, i), plat0(3, i)], 'LineStyle', '-', 'LineWidth', 4, 'Marker', 'None'); hold on
end
hold on
fill3(base(1, :), base(2, :), base(3, :), colors(8, :)); hold on
fill3(plat0(1, :), plat0(2, :), plat0(3, :), colors(5, :)); hold off
grid on
grid minor
title('Reference Result', 'Interpreter', 'LaTeX', 'FontSize', 14);

base_ = base.';
plat_ = plat.';
counter = 1;
base = [];
plat = [];
for i = 1 : 6
    base = [base; base_(i, :)];
    plat = [plat; plat_(i, :)];
    if(i < 6)
        tmp = (base_(i, :) + base_(i + 1, :)) / 2;
        base = [base; tmp];

        tmp = (plat_(i, :) + plat_(i + 1, :)) / 2;
        plat = [plat; tmp];
    end

end
base = base.';
plat = plat.';
len = size(base, 2);

leg0 = zeros(len, 1);
plat0 = zeros(3, len);
for i = 1 : len
    plat0(:, i) = R0 * plat(:, i)  + t0;
    res = base(:, i) - plat0(:, i);
    leg0(i) = sqrt(res.' * res);
end

syms q0 q1 q2 q3
q = [q0; q1; q2; q3];
syms t1 t2 t3;
t = [t1; t2; t3];
R = q2R(q);
syms r1 r2 r3 r4
r = [r1; r2; r3; r4];
rr = r(1 : 3);
eqs = sym(zeros(len + 3, 1));
for i = 1 : len
    eqs(i) = base(:, i).' * base(:, i) - 2 * plat(:, i).' * R.' * base(:, i) + ...
             plat(:, i).' * plat(:, i) - 2 * base(:, i).' * t + 2 * plat(:, i).' * rr + r4 - leg0(i)^2;
end
eqs(len + 1) = q.' * q - 1;
eqs(len + 2) = rr.' * rr - r4;
eqs(len + 3) = t.' * t - r4;
eqs = expand(eqs);
x = [q; t; r];
H = expand(jacobian(eqs, x).' * eqs);
assumeAlso(q.' * q == 1);
assumeAlso(rr.' * rr - r4 == 0);
assumeAlso(t.' * t - r4 == 0);
eq = vpa(expand(simplify(H)), 32);

eq_ = eq(1 : 4);
G = jacobian(eq(5 : 11), [t; r]);
ss = - pinv(G) * (eq(5 : 11) - G * [t; r]);
ss = neglect_tiny_terms(ss, 32);
ss = ss.';
eq_ = eq(1 : 4);
eq_ = subs(eq_, t1, ss(1));
eq_ = subs(eq_, t2, ss(2));
eq_ = subs(eq_, t3, ss(3));
eq_ = subs(eq_, r1, ss(4));
eq_ = subs(eq_, r2, ss(5));
eq_ = subs(eq_, r3, ss(6));
eq_ = subs(eq_, r4, ss(7));
t_func = matlabFunction([ss(1); ss(2); ss(3)], 'Vars', {q});
r_func = matlabFunction([ss(4); ss(5); ss(6)], 'Vars', {q});
r4_func = matlabFunction(ss(7), 'Vars', {q});
eq_ = vpa(expand(eval(eq_)), 32);
syms lambda
eq_ = [
    neglect_tiny_terms(eq_, 32).';
    q.' * q - 1;
    ]

eqs = [
    eq_(1 : 4) + lambda * q;
    eq_(5);
    ];
str = '';
for i = 1 : length(eqs)
    str = strcat(str, sprintf(' PP{%d} = char(vpa(%%s, 32));', i));
end

str_ = sprintf(str, char(eqs(1)), ...
                    char(eqs(2)), ...
                    char(eqs(3)), ...
                    char(eqs(4)), ...
                    char(eqs(5)));
eval(str_);
[S, vars] = psolve(PP);
S = S.';
SS = S;
for i = 1 : length(vars)
    if(strcmp(vars{i}, 'q0'))
        SS(:, 1) = S(:, i);
    elseif(strcmp(vars{i}, 'q1'))
        SS(:, 2) = S(:, i);
    elseif(strcmp(vars{i}, 'q2'))
        SS(:, 3) = S(:, i);
    elseif(strcmp(vars{i}, 'q3'))
        SS(:, 4) = S(:, i);
    elseif(strcmp(vars{i}, 'lambda'))
        SS(:, 5) = S(:, i);
    end
end
S = real(SS);
xs_ = S;
sols = SS.';

num = size(sols, 2);
sol = zeros(4, num);
ts = zeros(3, num);
Ls = 1e50 * ones(num, 1);
q_true = dcm2quat(R0).';
if(q_true(1) < 0)
    q_true = - q_true;
end
for i = 1 : num
    sol(:, i) = real(sols(1 : 4, i));
    sol(:, i) = sol(:, i) ./ norm(sol(:, i));
    if(sol(1, i) < 0)
        sol(:, i) = - sol(:, i);
    end
    C = q2R(sol(:, i));
    t = t_func(sol(:, i));
    r4 = r4_func(sol(:, i));
    ts(:, i) = t;
    res = abs(q_true - sol(:, i));
    loss = res.' * res;
    Ls(i) = loss;
end
[~, idx] = sort(Ls);

q_ = sol(:, idx(1)).'
q_true_ = q_true.'

R_ = quat2dcm(sol(:, idx(1)).');
t_ = t0;
plat0 = zeros(3, 6);
for i = 1 : 6
    plat0(:, i) = R_ * plat00(:, i)  + t_;
end
base = base00;

subplot(1, 2, 2);
plot3(base(1, :), base(2, :), base(3, :), 'LineStyle', 'None', 'Marker', '.', 'MarkerSize', 10); hold on
plot3(base(1, 1 : 6), base(2, 1 : 6), base(3, 1 : 6), 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
plot3([base(1, 1), base(1, 6)], [base(2, 1), base(2, 6)], [base(3, 1), base(3, 6)], 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
plot3(plat0(1, :), plat0(2, :), plat0(3, :), 'LineStyle', 'None', 'Marker', '.', 'MarkerSize', 10); hold on
plot3(plat0(1, 1 : 6), plat0(2, 1 : 6), plat0(3, 1 : 6), 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
plot3([plat0(1, 1), plat0(1, 6)], [plat0(2, 1), plat0(2, 6)], [plat0(3, 1), plat0(3, 6)], 'LineStyle', '-', 'LineWidth', 2, 'Marker', 'None'); hold on
for i = 1 : 6
    plot3([base(1, i), plat0(1, i)], [base(2, i), plat0(2, i)], [base(3, i), plat0(3, i)], 'LineStyle', '-', 'LineWidth', 4, 'Marker', 'None'); hold on
end
hold on
fill3(base(1, :), base(2, :), base(3, :), colors(8, :)); hold on
fill3(plat0(1, :), plat0(2, :), plat0(3, :), colors(5, :)); hold off
grid on
grid minor
title('QPEP Result', 'Interpreter', 'LaTeX', 'FontSize', 14);

if(~ispc())
    set(gcf, 'Position', [634 780 1159 320])
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function diginto(thestruct, level)

% DIGINTO  Pretty print entire structure trees.
%   DIGINTO prints the fields tree of all structures in the workspace.
%
%   DIGINTO('all') is the same as DIGINTO.
%
%   DIGINTO(STR) pretty prints the fields' tree of structure STR.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargin == 0 || strcmp(thestruct,'all')
    disp('Printing all structures in workspace ''base''.')
    w = evalin('base','whos');
    for i = 1:numel(w)
        if strcmp(w(i).class,'struct')
            disp(w(i).name)
            diginto(evalin('base',w(i).name))
        end
    end
    return
end

if nargin < 2
    level = 0;
end
if level == 0
    disp(inputname(1))
end

fn = fieldnames(thestruct);
for n = 1:length(fn)
    tabs = '';
    for m = 0:level
        tabs = [tabs '|   '];
    end
    disp([tabs '.' fn{n}])
    fn2 = thestruct.(fn{n});
    if isstruct(fn2)
        diginto(fn2, level+1);
    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
%  generate a random mxn matrix with entries in [-1,1]
%
%  syntax  >> A = Srand(m,n)


text:
function [X, r] = odePinv(A,varargin)
%
% the standard pinv has been modified to also return the rank of the matrix
%
%PINV   Pseudoinverse.
%   X = PINV(A) produces a matrix X of the same dimensions
%   as A' so that A*X*A = A, X*A*X = X and A*X and X*A
%   are Hermitian. The computation is based on SVD(A) and any
%   singular values less than a tolerance are treated as zero.
%   The default tolerance is MAX(SIZE(A)) * NORM(A) * EPS(class(A)).
%
%   PINV(A,TOL) uses the tolerance TOL instead of the default.
%
%   Class support for input A:
%      float: double, single
%
%   See also RANK.

%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 5.12.4.2 $  $Date: 2004/12/06 16:35:27 $

if isempty(A)     % quick return
  X = zeros(size(A'),class(A));
  return
end

[m,n] = size(A);

if n > m
   [X, r] = odePinv(A',varargin{:})';
else
   [U,S,V] = svd(A,0);
   if m > 1, s = diag(S);
      elseif m == 1, s = S(1);
      else s = 0;
   end
   if nargin == 2
      tol = varargin{1};
   else
      tol = max(m,n) * eps(max(s));
   end
   r = sum(s > tol);
   if (r == 0)
      X = zeros(size(A'),class(A));
   else
      s = diag(ones(r,1)./s(1:r));
      X = V(:,1:r)*s*U(:,1:r)';
   end
end


text:
function condition_test01 ( )

%*****************************************************************************80
%
%% CONDITION_TEST01 tests CONDITION_LINPACK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   12 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONDITION_TEST01\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  CONDITION_LINPACK estimates the L1 condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix               Order   Condition         Linpack\n' );
  fprintf ( 1, '\n' );
%
%  Combinatorial matrix.
%
  name = 'Combinatorial';
  n = 4;
  alpha = 2.0;
  beta = 3.0;
  a = combin ( alpha, beta, n );
  a_inverse = combin_inverse ( alpha, beta, n );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1 = a_norm_l1 * a_inverse_norm_l1;
  cond = condition_linpack ( n, a );
  fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
%
%  CONEX1
%
  name = 'CONEX1';
  n = 4;
  alpha = 100.0;
  a = conex1 ( alpha );
  a_inverse = conex1_inverse ( alpha );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1 = a_norm_l1 * a_inverse_norm_l1;
  cond = condition_linpack ( n, a );
  fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
%
%  CONEX2
%
  name = 'CONEX2';
  n = 3;
  alpha = 100.0;
  a = conex2 ( alpha );
  a_inverse = conex2_inverse ( alpha );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1 = a_norm_l1 * a_inverse_norm_l1;
  cond = condition_linpack ( n, a );
  fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
%
%  CONEX3
%
  name = 'CONEX3';
  n = 5;
  a = conex3 ( n );
  a_inverse = conex3_inverse ( n );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1 = a_norm_l1 * a_inverse_norm_l1;
  cond = condition_linpack( n, a );
  fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
%
%  CONEX4
%
  name = 'CONEX4';
  n = 4;
  a = conex4 ( );
  a_inverse = conex4_inverse ( );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1 = a_norm_l1 * a_inverse_norm_l1;
  cond = condition_linpack ( n, a );
  fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
%
%  KAHAN
%
  name = 'KAHAN';
  n = 4;
  alpha = 0.25;
  a = kahan ( alpha, n, n );
  a_inverse = kahan_inverse ( alpha, n );
  a_norm_l1 = norm ( a, 1 );
  a_inverse_norm_l1 = norm ( a_inverse, 1 );
  cond_l1 = a_norm_l1 * a_inverse_norm_l1;
  cond = condition_linpack ( n, a );
  fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
%
%  Random
%
  seed = 123456789;

  for i = 1 : 5
    name = 'RANDOM';
    n = 4;
    [ a, seed ] = r8mat_uniform_01 ( n, n, seed );
    a_inverse = inv ( a );
    a_norm_l1 = norm ( a, 1 );
    a_inverse_norm_l1 = norm ( a_inverse, 1 );
    cond_l1 = a_norm_l1 * a_inverse_norm_l1;
    cond = condition_linpack ( n, a );
    fprintf ( 1, '  %20s  %4d  %14.6g  %14.6g\n', name, n, cond_l1, cond );
  end

  return
end


text:
function v = zUpDown

v = 2*strcmpi(getMTEXpref('zAxisDirection'),'outOfPlane')-1;

end

author: mtex-toolbox

repo: mtex

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


text:
function V = LP3cvx(J, model)

n = size(model.S,2);

cvx_begin quiet

    variable v(n);

    maximize (ones(1,numel(J)) * v(J) );

    model.S*v==0;
    v>=model.lb;
    v<=model.ub;

cvx_end

V = v;

author: opencobra

repo: cobratoolbox


text:
function [cost,costs] = rcrossvalidate(model, L, wfun, estfct,combinefct)

% Estimate the model performance of a model with robust [$ l$] -fold crossvalidation

% CAUTION!! Use this function only to obtain the value of the rcrossvalidation score
% function given the tuning parameters. Do not use this function together with
% 'tunelssvm', but use 'rcrossvalidatelssvm' instead. The latter is a faster
% implementation which uses previously computed results.
%
%
% >> cost = rcrossvalidate({Xtrain,Ytrain,type,gam,sig2})
% >> cost = rcrossvalidate( model)
%
% Robustness in the $l$-fold crossvalidation score function is obtained by
% iteratively reweighting schemes.
%
% This routine is computational intensive.
%
%
% Some commonly used criteria are:
%
% >> cost = rcrossvalidate(model, 10, 'whuber', 'mae')
% >> cost = rcrossvalidate(model, 10, 'whampel', 'mae')
% >> cost = rcrossvalidate(model, 10, 'wlogistic', 'mae')
% >> cost = rcrossvalidate(model, 10, 'wmyriad', 'mae')
%
% Full syntax
%
%     1. Using LS-SVMlab with the functional interface:
%
% >> [cost, costs] = rcrossvalidate({X,Y,type,gam,sig2,kernel,preprocess}, L, wfun, estfct, combinefct)
%
%       Outputs
%         cost          : Cost estimation of the L-fold cross validation
%         costs(*)      : L x 1 vector with costs estimated on the L different folds
%       Inputs
%         X             : Training input data used for defining the LS-SVM and the preprocessing
%         Y             : Training output data used for defining the LS-SVM and the preprocessing
%         type          : 'function estimation' ('f') or 'classifier' ('c')
%         gam           : Regularization parameter
%         sig2          : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         kernel(*)     : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess'(*) or 'original'
%         L(*)          : Number of folds (by default 10)
%         wfun(*)       : weighting scheme (by default: whuber)
%         estfct(*)     : Function estimating the cost based on the residuals (by default mse)
%         combinefct(*) : Function combining the estimated costs on the different folds (by default mean)
%
%
%     2. Using the object oriented interface:
%
% >> [cost, costs] = crossvalidate(model, L, wfun, estfct, combinefct)
%
%       Outputs
%         cost          : Cost estimation of the L-fold cross validation
%         costs(*)      : L x 1 vector with costs estimated on the L different folds
%         ec(*)         : N x 1 vector with residuals of all data
%       Inputs
%         model         : Object oriented representation of the LS-SVM model
%         Xval          : Nt x d matrix with the inputs of the validation points used in the procedure
%         Yval          : Nt x m matrix with the outputs of the validation points used in the procedure
%         L(*)          : Number of folds (by default 10)
%         wfun(*)       : weighting scheme (by default: whuber)
%         estfct(*)     : Function estimating the cost based on the residuals (by default mse)
%         combinefct(*) : Function combining the estimated costs on the different folds (by default mean)
%
% See also:
% mae,whuber,wlogistic,whampel,wmyriad, crossvalidate, trainlssvm, robustlssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% LS-SVMlab
eval('model = initlssvm(model{:});',' ');
eval('L;','L=min(ceil(sqrt(model.nb_data)),10);');
eval('estfct;','estfct=''mae'';');
eval('combinefct;','combinefct=''mean'';');
eval('wfun;','wfun=''whuber'';');
%
% initialisation and defaults
%
nb_data = size(model.ytrain,1);

if L==nb_data, p = 1:nb_data; else p = randperm(nb_data); end
px = model.xtrain(p,:);
py = model.ytrain(p,:);

[~,Y] = postlssvm(model,[],py);

%initialize: no incremental  memory allocation
costs = zeros(L,length(model.gam));
block_size = floor(nb_data/L);

S = ones(nb_data,1);
Atot = kernel_matrix(px,model.kernel_type,model.kernel_pars)+eye(nb_data)./model.gam;
%
%
% start loop over l validations
%
for l = 1:L,

    % divide in data and validation set, trainings data set is a copy
    % of permutated_data, validation set is just a logical index
    if l==L,
        train = 1:block_size*(l-1);
        validation = block_size*(l-1)+1:nb_data;
    else
        train = [1:block_size*(l-1) block_size*l+1:nb_data];
        validation = block_size*(l-1)+1:block_size*l;
    end

    A = [0 S(train)';S(train) Atot(train,train)];
    b = [0;py(train)];

    % Solve linear system
    sol = linsolve(A,b,struct('SYM',true));

    % Determine residuals ek
    ek = sol(2:end)./model.gam;
    g = model.gam;
    %for i=2:size(A,1), A(i,i) = A(i,i) - 1/g; end
    A = A-eye(size(train,2)+1)./g; A(1,1)=0;
    Ah = A;
    %
    % robust estimation of the variance
    %
    for k = 1:20
        vare = 1.483*median(abs((ek)-median(ek)));
        alphaold = sol(2:end);
        %
        % robust re-estimation of the alpha's and the b
        %
        cases = reshape((ek./vare),1,size(ek,1));
        W = g*weightingscheme(cases,wfun);

        for t=1:size(train,2), A(t+1,t+1) = A(t+1,t+1)+1./W(t); end

        sol = linsolve(A,b,struct('SYM',true));

        ek = sol(2:end)./W';
        A = Ah;
        if norm(abs(alphaold-sol(2:end)),'fro')<=1e-4,
            %fprintf('\n Converged after %.0f iteration(s)', k);
            k = inf;
        end
        model.status = 'changed';
    end

    % regression
    % Simulate system on validation data
    yh = Atot(train,validation)'*sol(2:end) + ones(numel(validation),1)*sol(1);
    [~,yh] = postlssvm(model,[],yh);
    costs(l,1) = feval(estfct,yh - Y(validation,:));
end
cost = feval(combinefct, costs);

author: peterhcharlton

repo: RRest


text:
function j = jacobi_symbol ( q, p )

%*****************************************************************************80
%
%% JACOBI_SYMBOL evaluates the Jacobi symbol (Q/P).
%
%  Definition:
%
%    If P is prime, then
%
%      Jacobi Symbol (Q/P) = Legendre Symbol (Q/P)
%
%    Else
%
%      let P have the prime factorization
%
%        P = Product ( 1 <= I <= N ) P(I)^E(I)
%
%      Jacobi Symbol (Q/P) =
%
%        Product ( 1 <= I <= N ) Legendre Symbol (Q/P(I))^E(I)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996, pages 86-87.
%
%  Parameters:
%
%    Input, integer Q, an integer whose Jacobi symbol with
%    respect to P is desired.
%
%    Input, integer P, the number with respect to which the Jacobi
%    symbol of Q is desired.  P should be 2 or greater.
%
%    Output, integer J, the Jacobi symbol (Q/P).
%    Ordinarily, J will be -1, 0 or 1.
%    -2, not enough factorization space.
%    -3, an error during Legendre symbol calculation.
%

%
%  P must be greater than 1.
%
  if ( p <= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_SYMBOL - Fatal error!\n' );
    fprintf ( 1, '  P must be greater than 1.\n' );
    j = -2;
    return
  end
%
%  Decompose P into factors of prime powers.
%
  [ nfactor, factor, power, nleft ] = i4_factor ( p );

  if ( nleft ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_SYMBOL - Fatal error!\n' );
    fprintf ( 1, '  Not enough factorization space.\n' );
    j = -2;
    error ( 'JACOBI_SYMBOL - Fatal error!' );
  end
%
%  Force Q to be nonnegative.
%
  qq = q;

  while ( qq < 0 )
    qq = qq + p;
  end
%
%  For each prime factor, compute the Legendre symbol, and
%  multiply the Jacobi symbol by the appropriate factor.
%
  j = 1;
  for i = 1 : nfactor
    pp = factor(i);
    l = legendre_symbol ( qq, pp );
    if ( l < -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'JACOBI_SYMBOL - Fatal error!\n' );
      fprintf ( 1, '  Error during Legendre symbol calculation.\n' );
      j = -3;
      error ( 'JACOBI_SYMBOL - Fatal error!' );
    end
    j = j * l^power(i);
  end

  return
end


text:
%  Figure 7.23      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 7.23
% fig7_23.m
% Nyquist for LQR
clf;
f=[0 1;0 0];
g=[0;1];
h=[1 0];
j=0;
r=1;
rho=1;
q=rho*h'*h;
[k,s]=lqr(f,g,q,r);
sys=ss(f,g,k,0);
% circle of radius one
plot(-1+sin(0:.1:2*pi),cos(0:.1:2*pi),'r--');
axis equal;
text(-1,0,'x');
hold on;
w=logspace(0,2);
nyquist(sys,w);
title('Fig. 7.23: Nyquist diagram');


text:
function pnamefull = fullpath(pname, style)

pnamefull = '';
currdir = pwd;

if ~exist('pname','var')
    return;
end
if ~exist('style','var')
    style = 'linux';
end

if strcmp(style, 'linux')
    sep = '/';
else
    sep = filesep;
end

% If path to file wasn't specified at all, that is, if only the filename was
% provided without an absolute or relative path, the add './' prefix to file name.
if isempty(fileparts(pname))
    pname = ['./', pname];
end

% Fix performance issue: its expensive to call cd so avoid it if we can. A valid path is either
% relative or absolute. If pname is already absolute then no need to change folder or anything else
% we are done and simply exist function. We check full path by prepending './'. If a valid path is
% absolute then prepending a './' will make it invalid. JD, Mar 21, 2022
if ispathvalid(pname) && ~ispathvalid(['./', pname])
    pnamefull = pname;
    pnamefull(pnamefull=='/' | pnamefull=='\') = sep;
    return;
end

% If path is file, extract pathname
try

    [p0,f0,e0] = fileparts(pname);
    if ispathvalid(pname, 'file') || ~isempty(strfind(f0,'*')) || ~isempty(strfind(e0,'*')) %#ok<*STREMP>
        % We have a valid file path
        cd(p0);
        f = f0;
        e = e0;
    elseif ispathvalid(pname, 'dir')
        cd(pname);
        f = '';
        e = '';
    else
        return
    end

catch ME

    cd(currdir)
    % rethrow(ME)
    return

end
p = pwd;

pnamefull = [p,sep,f,e];
if ~ispathvalid(pnamefull)
    pnamefull = '';
    cd(currdir);
    return;
end

pnamefull(pnamefull=='/' | pnamefull=='\') = sep;

cd(currdir);

author: BUNPC


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%
for i=1:p
	X_poly(:,i)=X .^ i;
end

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

end

author: loserChen

repo: Coursera-MachineLearning


text:
function E = r8mat_expm2 ( A )

%*****************************************************************************80
%
%% R8MAT_EXPM2 uses the Taylor series for the matrix exponential.
%
%  Discussion:
%
%    Formally,
%
%      exp ( A ) = I + A + 1/2 A^2 + 1/3! A^3 + ...
%
%    This function sums the series until a tolerance is satisfied.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2011
%
%  Author:
%
%    Cleve Moler, Charles Van Loan
%
%  Reference:
%
%    Cleve Moler, Charles VanLoan,
%    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
%    Twenty-Five Years Later,
%    SIAM Review,
%    Volume 45, Number 1, March 2003, pages 3-49.
%
%  Parameters:
%
%    Input, real A(N,N), the matrix.
%
%    Output, real E(N,N), the estimate for exp(A).
%
  E = zeros ( size ( A ) );
  F = eye ( size ( A ) );
  k = 1;

  while ( 0 < norm ( E + F - E, 1 ) )
    E = E + F;
    F = A * F / k;
    k = k + 1;
  end

  return
end


text:
function saveplots(fmri_model, varargin)

if isempty(varargin) % Output dir
    savedir = pwd;
else
savedir = varargin{1};
end

if ~exist(savedir, 'dir'), mkdir(savedir); end

% Register of possible plot names associated with this object type

plotnames = {'basis sets' ...
    'design matrix image' ...
    'design matrix' ...
    };

for j = 1:length(plotnames)

    han = findobj('Name', plotnames{j});
    if ~isempty(han) && ishandle(han) && strcmp(get(han, 'Type'), 'figure')

        fprintf('Saving figure: %s\n', plotnames{j});
        figure(han);
        scn_export_papersetup(500);

        name = get(han, 'Name');
        name(name == ' ') = '_';

        name = fullfile(savedir, name);

        saveas(han, name, 'png');

    end

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% POP_RMDAT - Keep or remove continuous data around specific events
%
% Usage:
%   >> OUTEEG = pop_rmdat( EEG); % pop-up a data entry window
%   >> OUTEEG = pop_rmdat( EEG, typerange, timelimits, invertselection);
%
% Graphic interface:
%   "Time-locking event type(s)" - [edit box] Select 'Edit > Event values'
%                to see a list of event.type values; else use the push button.
%                To use event types containing spaces, enter in single-quotes.
%   "..."      - [push button] scroll event types.
%   "Time limits" - [edit box] epoch latency range [start, end] in seconds relative
%                to the event type latency.
%
% Inputs:
%   EEG        - Input dataset. Data may already be epoched; in this case,
%                extract (shorter) subepochs time locked to epoch events.
%   typerange  - Cell array of event types to time lock to.
%                (Note: An event field called 'type' must be defined in
%                the 'EEG.event' structure).
%   timelimits - Epoch latency limits [start end] in seconds relative to the time-locking event
%                {default: [-1 2]}%
%   invertselection - [0|1] Invert selection. 0 is to keep data around events
%                and 1 is to remove them {default:1 to remove them}
%
% Outputs:
%   OUTEEG     - output dataset
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, 2009-

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG, com] = pop_rmdat( EEG, events, timelims, invertsel )

if nargin < 1
   help pop_rmdat;
	return;
end
com = '';
if nargin < 4
    invertsel = 1;
    warning('This function behavior has been fixed and data is removed by default')
end
if isempty(EEG(1).event)
    error( [ 'No event. This function removes data' 10 'based on event latencies' ]);
end
if isempty(EEG(1).trials)
    error( [ 'This function only works with continuous data' ]);
end
if ~isfield(EEG(1).event, 'latency')
    error( 'Absent latency field in event array/structure: must name one of the fields ''latency''');
end

if nargin < 3
	% popup window parameters
	% -----------------------
   promptstr    = { strvcat('Time-locking event type(s) ([]=all):', ...
                    'Select ''Edit > Event values'' to see type values.'), ...
                    'Epoch limits [start, end] in seconds:', ...
                    'Name for the new dataset:', ...
					'Out-of-bounds EEG rejection limits ([min max], []=none):'  };

   cbevent = ['tmpEEG = get(gcbf, ''userdata'');' ...
           'if ~isfield(tmpEEG.event, ''type'')' ...
				   '   errordlg2(''No type field'');' ...
				   'else' ...
                   '   tmpevent = tmpEEG.event;' ...
                   '   if isnumeric(tmpEEG.event(1).type),' ...
				   '        [tmps,tmpstr] = pop_chansel(unique([ tmpevent.type ]));' ...
				   '   else,' ...
                   '        [tmps,tmpstr] = pop_chansel(unique({ tmpevent.type }));' ...
                   '   end;' ...
				   '   if ~isempty(tmps)' ...
				   '       set(findobj(''parent'', gcbf, ''tag'', ''events''), ''string'', tmpstr);' ...
				   '   end;' ...
				   'end;' ...
				   'clear tmpEEG tmps tmpv tmpstr tmpevent tmpfieldnames;' ];

   geometry = { [2 1 1.2] [2 1 1.2] };
   uilist = { { 'style' 'text'       'string' 'Event type(s) ([]=all)' } ...
              { 'style' 'edit'       'string' '' 'tag' 'events' } ...
              { 'style' 'pushbutton' 'string' '...' 'callback' cbevent } ...
              { 'style' 'text'       'string' 'Time limits [start, end] in sec.' } ...
              { 'style' 'edit'       'string' '-1 1' } ...
              { 'style' 'popupmenu'  'string' 'Keep selected|Remove selected' } };

   result = inputgui( 'geometry', geometry, 'uilist', uilist, 'helpcom', 'pophelp(''pop_rmdat'')', 'title', 'Remove data portions around events - pop_rmdat()', 'userdata', EEG);
   if isempty(result) return; end

   if strcmpi(result{1}, '[]'), result{1} = ''; end
   if ~isempty(result{1})
       if strcmpi(result{1}(1),'''')   % If event type appears to be in single-quotes, use comma
                                       % and single-quote as delimiter between event types. toby 2.24.2006
                                       % fixed Arnaud May 2006
            events = eval( [ '{' result{1} '}' ] );
       else events = parsetxt( result{1});
       end
   else events = {};
   end
   timelims = eval( [ '[' result{2} ']' ] );
   invertsel = result{3}-1;

end

% process multiple datasets
% -------------------------
if length(EEG) > 1
    if nargin < 3
        [ EEG, com ] = eeg_eval( 'pop_rmdat', EEG, 'warning', 'on', 'params', { events, timelims, invertsel } );
    else
        [ EEG, com ] = eeg_eval( 'pop_rmdat', EEG, 'params', { events, timelims, invertsel } );
    end
    return;
end

tmpevent = EEG.event;

% Checking for numeric values in fieldnames and changing them to string
 for i = 1: length(tmpevent)
     checknum(i) = isnumeric(tmpevent(i).type);
 end
  if sum(checknum)~=0
       tmpindx = find(checknum == 1);
       for i = 1: length(tmpindx)
           tmpevent(tmpindx(i)).type = num2str(tmpevent(tmpindx(i)).type);
       end
  end

% compute event indices
% ---------------------
allinds = [];
for index = 1:length(events)
    inds = strmatch(events{index},{ tmpevent.type }, 'exact');
    allinds = [allinds(:); inds(:) ]';
end
allinds = sort(allinds);
if isempty(allinds)
    disp('No event found');
    return;
end

% compute time limits
% -------------------
array = [];
bnd = strmatch('boundary', lower({tmpevent.type }));
bndlat = [ tmpevent(bnd).latency ];
for bind = 1:length(allinds)
    evtlat = EEG.event(allinds(bind)).latency;
    evtbeg = evtlat+EEG.srate*timelims(1);
    evtend = evtlat+EEG.srate*timelims(2);
    if any(bndlat > evtbeg & bndlat < evtend)
        % find the closer upper and lower boundaries
        bndlattmp = bndlat(bndlat > evtbeg & bndlat < evtend);
        diffbound = bndlattmp-evtlat;
        allneginds = find(diffbound < 0);
        allposinds = find(diffbound > 0);
        if ~isempty(allneginds), evtbeg = bndlattmp(allneginds(1)); end
        if ~isempty(allposinds), evtend = bndlattmp(allposinds(1)); end;
        fprintf('Boundary found: time limits for event %d reduced from %3.2f to %3.2f\n', allinds(bind), ...
            (evtbeg-evtlat)/EEG.srate, (evtend-evtlat)/EEG.srate);
    end
    if ~isempty(array) && evtbeg < array(end)
        array(end) = evtend;
    else
        array = [ array; evtbeg  evtend];
    end
end

if ~isempty(array) && array(1) < 1, array(1) = 1; end
if ~isempty(array) && array(end) > EEG.pnts, array(end) = EEG.pnts; end
if isempty(array)
    disp('No event found');
    return;
end

if invertsel
    EEG = pop_select(EEG, 'notime', (array-1)/EEG.srate);
else
    EEG = pop_select(EEG, 'time', (array-1)/EEG.srate);
end

% generate output command
% -----------------------
com = sprintf('EEG = pop_rmdat( EEG, %s);', vararg2str( { events timelims invertsel } ));

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function out = lsb_dec( wavin, pass )
%LSB_DEC LSB Coding technique for WAV files using one bit to hide
%
%   INPUT VARIABLES
%       wavin : Path of stego signal (wav input)
%       pass  : Password for decryption
%
%   OUTPUT VARIABLES
%       out   : Retrieved message
%
%   Kadir Tekeli (kadir.tekeli@outlook.com)

if (nargin<2)
    pass = 'password123';
end

%Header = 1:40, Length = 41:43, Data = 44:end
fid = fopen(wavin,'r');
header = fread(fid,40,'uint8=>char');
dsize  = fread(fid,1,'uint32');
stego  = fread(fid,inf,'uint16');
fclose(fid);
%Control variable is read (1:8)
control = bitget(stego(1:8),1)';

if b2d(control) == mod(sum(double(pass)),256)
    %Length of message is read (9:48)
    m = bitget(stego(9:48),1);
    len = b2d(m')*8;
    %Hidden message is read and decrypted (49:len)
    dat = xor(bitget(stego(49:48+len),1),prng(pass,len));
    bin = reshape(dat,len/8,8);
    out = char(b2d(bin))';
else
    warning('Password is wrong or message is corrupted!');
    out = [];
end
end

function d = b2d(b)
%B2D Minimal implentation of bi2de function
  if isempty(b)
    d = [];
  else
    d = b * (2 .^ (0:length(b(1,:)) - 1)');
  end
end

function out = prng( key, L )
%PRNG Pseudorandom number generator for encryption/decryption
pass = sum(double(key).*(1:length(key)));
rand('seed', pass);
out = (rand(L, 1)>0.5);
end

author: ktekeli

repo: audio-steganography-algorithms


text:
function model = fgplvmExpandParam(model, params)

% FGPLVMEXPANDPARAM Expand a parameter vector into a GP-LVM model.
% FORMAT
% DESC takes an FGPLVM structure and a vector of parameters, and
% fills the structure with the given parameters. Also performs any
% necessary precomputation for likelihood and gradient
% computations, so can be computationally intensive to call.
% ARG model : the FGPLVM structure to put the parameters in.
% ARG params : parameter vector containing the parameters to put in
% the FGPLVM structure.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2009
%
% MODIFICATION : Carl Henrik Ek, 2009
%
% SEEALSO : fgplvmCreate, fgplvmExtractParam, modelExpandParam

% FGPLVM

startVal = 1;
if isfield(model, 'back') & ~isempty(model.back)
  % update modelParameters
  endVal = model.back.numParams;
  model.back = modelExpandParam(model.back, params(startVal:endVal));

  % update latent locations
  tmp = modelOut(model.back,model.y);
  tmp_dim = 1;
  for(i = 1:1:model.q)
    if(length(find(model.back.indexOut==i))~=0)
      model.X(:,model.back.indexOut(tmp_dim)) = tmp(:,tmp_dim);
      tmp_dim = tmp_dim + 1;
    else
      startVal = endVal + 1;
      endVal = endVal + model.N;
      model.X(:,i) = reshape(params(startVal:endVal),model.N,1);
    end
  end
  clear tmp tmp_dim;
else
  endVal = model.N*model.q;
  model.X = reshape(params(startVal:endVal), model.N, model.q);
end
startVal = endVal+1;
endVal = endVal + model.kern.nParams;

switch model.approx
 case 'ftc'
  endVal = endVal;
 case {'dtc', 'dtcvar', 'fitc', 'pitc'}
  if model.fixInducing
    % account for beta attached to the end.
    endVal = endVal + 1;
    % X_u values are taken from X values.
    model.X_u = model.X(model.inducingIndices, :);
  else
    % Parameters include inducing variables and beta.
    endVal = endVal + model.q*model.k + 1;
  end

 otherwise
  error('Unknown approximation type.')
end
if model.learnScales
  endVal = endVal + model.d;
end
model = gpExpandParam(model, params(startVal:endVal));

% Give parameters to dynamics if they are there.
if isfield(model, 'dynamics') & ~isempty(model.dynamics)
  startVal = endVal + 1;
  endVal = length(params);

  % Fill the dynamics model with current latent values.
  model.dynamics = modelSetLatentValues(model.dynamics, model.X);

  % Update the dynamics model with parameters (thereby forcing recompute).
  model.dynamics = modelExpandParam(model.dynamics, params(startVal:endVal));
end

% Constraints
if(isfield(model,'constraints')&&~isempty(model.constraints))
  for(i = 1:1:model.constraints.numConstraints)
    model.constraints.comp{i} = constraintExpandParam(model.constraints.comp{i},model.X);
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function X=MergeMargCop(W,F,U)

[J,K]=size(W);
X=0*U;
for k=1:K
    dd = interp1(F(:,k),W(:,k),U(:,k),'linear','extrap');
    X(:,k) = dd;
end

author: Sable


text:
function [xV, yV, segNum] = redoSection(hCERR)
%function [xV, yV, segNum] = redoSection(hCERR)
%'Snip' & Redo a section of the active contour.
%Directions:  click (Button 1) once where the redone section should begin,
%             then click again where the section should end.  After the section
%             is deleted, fill in from the first point to the second, with
%             Button 1 creating points and button 2 closing the contour.
%
%Called by contourControls.
%LM:  2 April 02, JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

figure(hCERR)

%Get the start point:

hPlot = getActiveSegment;
segS = get(hPlot,'userdata');
xV = segS.GUI_xV;
yV = segS.GUI_yV;
segNum = segS.segNum;

hold on
[x, y, button]  = ginput(1);

if isempty(button)
  xV = [];
  yV = [];
  return
end

switch button  %Maybe more options in the future?

  case 1

    rV = ((x - xV).^2 + (y - yV).^2 ).^0.5;
    indV = find(rV == min(rV));
    ind_a = indV(1);

end

%Now get the ending point:
[x, y, button]  = ginput(1);

if isempty(button)
  xV = [];
  yV = [];
  return
end

switch button

  case 1 %get nearest point

    rV = ((x - xV).^2 + (y - yV).^2 ).^0.5;
    indV = find(rV == min(rV));
    ind_b = indV(1);

end

%reorganize xV and yV:
%ind_b becomes the first point, and ind_a becomes the last point.
%Also: snip the shortest path in terms of number of points.
len = length(xV);
if ind_b > ind_a
  %which path is smallest?
  path1 = ind_b - ind_a;
  path2 = ind_a + (len - ind_b);
  if path1 <= path2
    xV = [xV(ind_b:end),xV(1:ind_a)];
    yV = [yV(ind_b:end),yV(1:ind_a)];
  else
    xV = [xV(ind_a:ind_b)];
    yV = [yV(ind_a:ind_b)];
    xV = fliplr(xV);
    yV = fliplr(yV);
  end
else  %ind_a > ind_b
  %which path is smallest?
  path1 = ind_a - ind_b;
  path2 = ind_b + (len - ind_a);
  if path1 <= path2
    xV = [xV(ind_a:end),xV(1:ind_b)];
    yV = [yV(ind_a:end),yV(1:ind_b)];
    xV = fliplr(xV);
    yV = fliplr(yV);
  else
    xV = [xV(ind_b:ind_a)];
    yV = [yV(ind_b:ind_a)];
  end
end

%Plot remaining segment
delete(hPlot)
hPlot2 = plot(xV,yV,'rs-');  %plot closed contour
set(hPlot2,'markersize',4)

%Go into new segment mode:
xLast = xV(end);
yLast = yV(end);

finished = 0;
hLineV = [];
while ~finished
  [x, y, button] = ginput(1);
  if button == 3
    finished = 1;
  else
    xV = [xV, x];
    yV = [yV, y];
    h = plot([xLast, x],[yLast, y],'rs-');
    set(h,'markersize',4)
    hLineV = [hLineV, h];
    xLast = x;
    yLast = y;
  end

end

delete(hLineV)
delete(hPlot2)

x2V = [xV, xV(1)];
y2V = [yV, yV(1)];

hPlot = plot(x2V,y2V,'rs-');
set(hPlot,'markersize',4)

hold off

%segmentsS contains all the information about the plotted segments - it gets written into
%the userdata portion of any segment.

segS.active = 1;
segS.hPlot  = hPlot;
segS.GUI_xV = xV;
segS.GUI_yV = yV;

set(hPlot,'tag','newSegment')
set(hPlot,'userdata',segS)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function refresh()

%REFRESH refresh all paths in GISMO
% REFRESH refresh all paths in GISMO. This is useful when adding new
% functions and packages. These are not always recognized on the fly by
% Matlab.
%
% See also admin.which admin.remove admin.getpath

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

import admin.*

which(getpath);

author: geoscience-community-codes

repo: GISMO


text:
function decomp = SFC_3D(path, obs, boundary)
    % Initialize MultipleSegments
	local_bbox = boundary; offset_x = 0;
    decomp = MultipleSegments(local_bbox);
    decomp.set_obs(obs);
    decomp.dilate(path, offset_x);
end

author: LenaShengzhen

repo: AerialRobotics


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

function source = tldInitFirstFrame(tld,source,min_win)

% load the first frame into memory
source.im0  = img_get(source,source.idx(1));

% set the initial bounding box:
% - from file
if source.camera == 0 && exist([source.input '/init.txt'],'file')
    bb = dlmread([source.input '/init.txt']);
    source.bb = bb(:);

    % check
    if isempty(source.bb) || min(bb_size(source.bb)) < min_win
        exit('Error: bounding box is incorrectly defined or too small');
    end
% - by mouse
else
    source.bb = bb_click(tld,source.im0.input);

    % check
    if isempty(source.bb) || min(bb_size(source.bb)) < min_win
       source = [];
    end
end

repo: OpenTLD


text:
function i4_uniform_ab_test ( )

%*****************************************************************************80
%
%% I4_UNIFORM_AB_TEST tests I4_UNIFORM_AB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2014
%
%  Author:
%
%    John Burkardt
%
  a = -100;
  b = 200;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_UNIFORM_AB_TEST\n' );
  fprintf ( 1, '  I4_UNIFORM_AB computes pseudorandom values\n' );
  fprintf ( 1, '  in an interval [A,B].\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The lower endpoint A = %f\n', a );
  fprintf ( 1, '  The upper endpoint B = %f\n', b );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 20

    [ j, seed ] = i4_uniform_ab ( a, b, seed );

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

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This M file creates the "LossSpec_IGBT_Library.mat " file
% containing loss and thermal IGBT models provided in the
% "LossModelib.mdl' libray.
% IGBT specifications are saved in the "LossSpec_IGBT" structure.
%
% Data provided below for three different commercial
% IGBT/Diode half-bridge packs are extracted from manufacturer data sheets
% (see pdf files provided in the same directory).
%
% How to to add your own IGBT specifications
% ------------------------------------------
%   1)  Open the "LossModelib.mdl' library
%   2)  Duplicate a cell containing specifications of one existing IGBT
%       and insert it at the end of this file
%       (see "Add your own specifications here" cell)
%   3)  Modify this new set of data according to your specifications
%   4)  Run this script file in order to update the "LossSpec_IGBT_Library.mat " file
%   5)  In the "LossModelib.mdl' library
%            - select the IGBT or half-bridge model to update
%            - open the Mask Editor
%            - select "IGBT_Type" parameter
%            - In the "Popups" window, enter a new description
%              after existing IGBT descriptions
%   6) Save the "LossModelib.mdl' library.
% Note: Future release of SPS will provide an easy-to-use GUI to add or
% edit device characteristics.

% Pierre Giroux, Hydro-Quebec (IREQ), March 2012

% Ic=  Collector current
% Vce= Collector-emitter voltage
% Eon= Turn-on switching energy (mJ)
% Eoff= Turn-off switching energy (mJ)
% Tj= Junction temperature (degrees Celsius)
% Vcc= Supply voltage (V)
% Note1: IGBT switching energies are function of Vcc, Ic and Tj
% Note2: IGBT on-state losses = Vce (function of Ic and Tj) * Ic
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
clear all

%% Fuji Electric IGBT Module, 2 in one-package, 600V, 150A
Type='IGBT';
Manufacturer= 'Fuji Electric';
PartNo= '2MBI150U2A 060';
Description= 'IGBT Module, 2 in one-package, 600V, 150A';
%
% Typical turn-on (Eon)switching energies vs Ic
Vcc_Eon= 300;
Tj_Eon=  [  25 125 ];
Ic_Eon=  [  0  25    50    75    100   125   150   175   200   225  ];
Eon=     [  0  0.85  1.62  2.40  3.18  4.02  4.87  5.83  6.81  7.84
            0  1.12  2.15  3.20  4.27  5.49  6.67  8.18  9.85  11.6 ];
% Typical turn-off (Eoff)switching energies vs Ic
Vcc_Eoff= 300;
Tj_Eoff=  [ 25 125 ];
Ic_Eoff=  [ 0  25    50    75    100   125   150   175   200   225  ];
Eoff=     [ 0  0.66  1.33  2.13  3.18  4.37  5.46  6.90  8.34  9.92
            0  0.87  1.97  2.99  4.27  5.49  6.90  8.43  10.1  11.8 ];
% Typical on-state characteristics
Tj_OnState=  [ 25 125 ];
Ic_OnState = [ 0   1      5     10    15    33    50    75    100   150   200   300 ];
Vce_OnState= [ 0   0.64   0.74  0.80  0.87  1.08  1.19  1.37  1.51  1.77  2.03  2.48
               0   0.48   0.66  0.75  0.87  1.08  1.26  1.50  1.68  2.02  2.35  3.01 ];
%-------------------------------------------------------------------------
%
% Thermal Impedance:
%
Rth_jc=0.25;    % Junction-to-Case thermal resistance(K/W)
% Based on transient thermal impedance vs. time curve
Cth_j=0.27;   % Junction thermal capacitance (Joule/Kelvin)
%
%-------------------------------------------------------------------------
LossSpec_IGBT(1)= struct('Type',Type, 'Manufacturer',Manufacturer, 'PartNo',PartNo, ...
    'Description',Description, 'Vcc_Eon',Vcc_Eon, ...
    'Tj_Eon',Tj_Eon, 'Ic_Eon',Ic_Eon, 'Eon',Eon, 'Vcc_Eoff',Vcc_Eoff, ...
    'Tj_Eoff',Tj_Eoff, 'Ic_Eoff',Ic_Eoff, 'Eoff',Eoff, ...
    'Tj_OnState',Tj_OnState, 'Ic_OnState',Ic_OnState, 'Vce_OnState',Vce_OnState, ...
    'Rth_jc',Rth_jc, 'Cth_j',Cth_j);
clear Vcc_Eon Tj_Eon Ic_Eon Eon Vcc_Eoff Tj_Eoff Ic_Eoff Eoff Tj_OnState ...
      Ic_OnState Vce_OnState Manufacturer PartNo Description Type ...
      Rth_jc Cth_j
%=========================end of specifications===========================

%% ABB IGBT Module, ABB HiPak, 1700V, 800A

Type='IGBT';
Manufacturer= 'ABB';
PartNo= '5SNE 0800M170100';
Description= 'IGBT Module, ABB HiPak, 1700V, 800A';

% For Tj=25 deg. C, manufacturer specifications provide Eon and Eoff at only one current (800 A).
% The Eon and and Eoff values @ Tj=25 and 125 deg. C, 800 A specified below
% are used to deduce the Eon(Ic) and Eoff(Ic) curves @ Tj=25 deg. C,
% from the Eon(Ic) and Eoff(Ic) curves @ Tj=125 deg. C
% assuming that Eon and Eoff stay proportionnal to their values @ Tj= 125 deg. C

Tj_800A=   [ 25   125 ];  % Junction temperature (degrees Celsius)
Eon_800A=  [ 160  250 ];  % Turn-on energy at Tj_800A (mJ) for Ic=800A
Eoff_800A= [ 220  300 ];  % Turn-off energy at Tj_800A (mJ) for Ic=800A
%
% Typical turn-on (Eon)switching energies vs Ic
Vcc_Eon= 900;
Tj_Eon=  [ 25 125 ];
Ic_Eon=    [ 0 198  399 501 686 800 900 999 1100 1200 1300 1440 1550 1620 ];
% Eon @ 125 C (per data sheet curves)
Eon(2,:)=  [ 0 73.4 108 136 202 250 300 350 413  480  551  652  750  810  ];
% Eon @ 25 C (calculated using data sheet spec at 25 C)
Eon(1,:)=Eon_800A(1)/Eon_800A(2) * Eon(2,:);
%
% Typical turn-off (Eoff)switching energies vs Ic
Vcc_Eoff= 900;
Tj_Eoff=  [ 25 125 ];
Ic_Eoff=    [ 0 198  399 501 686 800 900 999 1100 1200 1300 1440 1550 1620 ];
% Eon @ 125 C (per data sheet curves)
Eoff(2,:)=  [ 0 99.7 161 192 255 300 335 372 413  451  497  558  608  640  ];
% Eon @ 25 C (calculated using data sheet spec at 25 C)
Eoff(1,:)=Eoff_800A(1)/Eoff_800A(2) * Eoff(2,:);
%
% IGBT - Typical on-state characteristics
Tj_OnState=  [ 25 125 ];      % (C)
Ic_OnState = [ 0  0.01  46   102  194  299  399  498  601  700  797  999  1200 1500 ];
Vce_OnState= [ 0  0.70  0.99 1.19 1.38 1.57 1.72 1.86 2.00 2.13 2.25 2.50 2.74 3.11
               0  0.70  0.87 1.12 1.38 1.65 1.85 2.04 2.23 2.41 2.58 2.93 3.28 3.82 ];
%-------------------------------------------------------------------------
% Thermal properties:
%
Rth_jc=0.021;   % Junction-to-Case thermal resistance(K/W)
%
% Based on transient thermal impedance vs. time curve
Cth_j=6.5095;     % Junction thermal capacitance (Joule/Kelvin)
%
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
LossSpec_IGBT(2)= struct('Type',Type, 'Manufacturer',Manufacturer, 'PartNo',PartNo, ...
    'Description',Description, 'Vcc_Eon',Vcc_Eon, ...
    'Tj_Eon',Tj_Eon, 'Ic_Eon',Ic_Eon, 'Eon',Eon, 'Vcc_Eoff',Vcc_Eoff, ...
    'Tj_Eoff',Tj_Eoff, 'Ic_Eoff',Ic_Eoff, 'Eoff',Eoff, ...
    'Tj_OnState',Tj_OnState, 'Ic_OnState',Ic_OnState, 'Vce_OnState',Vce_OnState, ...
    'Rth_jc',Rth_jc, 'Cth_j',Cth_j);
%
clear Vcc_Eon Tj_Eon Ic_Eon Eon Vcc_Eoff Tj_Eoff Ic_Eoff Eoff Tj_OnState ...
      Ic_OnState Vce_OnState Manufacturer PartNo Description Type ...
      Rth_jc Cth_j Tj_800A Eon_800A Eoff_800A
%=========================end of specifications============================

%% ABB  Half-bridge IGBT Module, 3300V, 250A
Type='IGBT';
Manufacturer= 'ABB';
PartNo= '5SNG 0250P330300';
Description= 'Half-bridge IGBT Module, 3300V, 250A';

% For Tj=25 deg. C, manufacturer specifications provide Eon and Eoff at only one current (250 A).
% The Eon and and Eoff values @ Tj=25 and 125 deg. C, 250 A specified below
% are used to deduce the Eon(Ic) and Eoff(Ic) curves @ Tj=25 deg. C,
% from the Eon(Ic) and Eoff(Ic) curves @ Tj=125 deg. C
% assuming that Eon and Eoff stay proportionnal to their values @ Tj= 125 deg. C

Tj_250A=   [ 25    125 ];  % Junction temperature (degrees Celsius)
Eon_250A=  [ 330   425 ];  % Turn-on energy at Tj_250A (mJ) for Ic=250A
Eoff_250A= [ 330   450 ];  % Turn-off energy at Tj_250A (mJ) for Ic=250A
%
% Typical turn-on (Eon)switching energies vs Ic
Vcc_Eon= 1800;
Tj_Eon=  [ 25 125 ];
Ic_Eon=   [ 0 37.8  99   199  250  300  349  400  451  499  ];
% Eon @ 125 C (per data sheet curves)
Eon(2,:)= [ 0 83.6  160  318  425  522  643  775  919  1070 ];
% Eon @ 25 C (calculated using data sheet spec at 25 C)
Eon(1,:)=Eon_250A(1)/Eon_250A(2) * Eon(2,:);
%
% Typical turn-off (Eoff)switching energies vs Ic
Vcc_Eoff= 1800;
Tj_Eoff=  [ 25 125 ];
Ic_Eoff=    [ 0 37.8  99   199  250  300  349  400  451  499  ];
% Eon @ 125 C (per data sheet curves)
Eoff(2,:)=  [ 0 118   218  381  450  543  622  699  782  861  ];
% Eon @ 25 C (calculated using data sheet spec at 25 C)
Eoff(1,:)=Eoff_250A(1)/Eoff_250A(2) * Eoff(2,:);
%
% IGBT - Typical on-state characteristics
Tj_OnState=  [ 25 125 ];      % (C)
Ic_OnState=  [ 0   1      16    77    131   174   276   375   500  ];
Vce_OnState= [ 0   0.51   1.00  1.50  1.81  2.04  2.49  2.88  3.37
               0   0.50   1.00  1.71  2.18  2.50  3.11  3.77  4.52 ];
%-------------------------------------------------------------------------
% Thermal properties:
%
Rth_jc=0.051;   % Junction-to-Case thermal resistance(K/W)
%
% Based on transient thermal impedance vs. time curve
Cth_j=2.6804;     % Junction thermal capacitance (Joule/Kelvin)
%
%-------------------------------------------------------------------------
LossSpec_IGBT(3)= struct('Type',Type, 'Manufacturer',Manufacturer, 'PartNo',PartNo, ...
    'Description',Description, 'Vcc_Eon',Vcc_Eon, ...
    'Tj_Eon',Tj_Eon, 'Ic_Eon',Ic_Eon, 'Eon',Eon, 'Vcc_Eoff',Vcc_Eoff, ...
    'Tj_Eoff',Tj_Eoff, 'Ic_Eoff',Ic_Eoff, 'Eoff',Eoff, ...
    'Tj_OnState',Tj_OnState, 'Ic_OnState',Ic_OnState, 'Vce_OnState',Vce_OnState, ...
    'Rth_jc',Rth_jc, 'Cth_j',Cth_j);
%
clear Vcc_Eon Tj_Eon Ic_Eon Eon Vcc_Eoff Tj_Eoff Ic_Eoff Eoff Tj_OnState ...
      Ic_OnState Vce_OnState Manufacturer PartNo Description Type ...
      Rth_jc Cth_j Tj_250A Eon_250A Eoff_250A
%=========================end of specifications===========================

%% Add your own specifications here
% Type='IGBT';
% Manufacturer= 'My Manufacturer';
% PartNo= 'XXX';
% Description= 'My IGBT Module, YYY V, ZZZ A;
% ....

%% Save IBGT thermal data structure "LossSpec_IGBT" in "LossSpec_IGBT_Library.mat" file
save LossSpec_IGBT_Library LossSpec_IGBT
%----------------------------END OF SCRIPT FILE--------------------------------;

repo: mcbench-benchmarks


text:
function c=samplebase
% samplebase constructor

c=class(struct('dummy',[]),'samplebase');


text:
%% Analyze Sentiment Using FinBERT
% This example shows how to analyze sentiment of financial data using a
% pretrained FinBERT model.
%
% FinBERT is a sentiment analysis model trained on financial text data and
% fine-tuned for sentiment analysis. The model is based on the BERT-Base
% architecture.
%
% This example shows how to classify the sentiment of financial news
% reports using a pretrained FinBERT model.

%% Load Pretrained FinBERT Model
% Load a pretrained FinBERT model using the |finbert| function. The model
% consists of a tokenizer that encodes text as sequences of integers, and
% a structure of parameters.
mdl = finbert

%%
% View the FinBERT model tokenizer. The tokenizer encodes text as sequences
% of integers and holds the details of padding, start, separator and mask
% tokens.
tokenizer = mdl.Tokenizer

%% Analyze Sentiment in Text
% Create an string array of text data.
str = [
    "In an unprecendented move the stock market has hit new records today following news of a new vaccine."
    "Businesses in this sector suffer dramatic losses on the back of the pandemic."
    "The ship unloader is totally enclosed along the entire conveying line to the storage facilities."
    "Experts estimate the value of its remaining stake in the company at $ 27 million."
    "The company said that sales in the three months to the end of March slid to EUR86 .4 m US$ 113.4 m from EUR91 .2 m last year."
    "Finance experts calculate that it has lost EUR 4mn in the failed project."
    "They signed a large deal with an international industrial group which will deliver automation solutions and connectivity services."
    "Operating profit rose to EUR 5mn from EUR 2.8 mn in the fourth quarter of 2008"];

%%
% Encode the text data as a sequence of tokens using the FinBERT model
% tokenizer.
tokens = encode(tokenizer,str);

%%
% Pad the sequences to have the same length using the |padsequences|
% function. Specify the padding value to match the FinBERT model tokenizer.
X = padsequences(tokens,2,"PaddingValue",mdl.Tokenizer.PaddingCode);

%%
% Evaluate the sentiment and sentiment scores using the
% |finbert.sentimentModel| function.
[sentiment,scores] = finbert.sentimentModel(X,mdl.Parameters)

author: matlab-deep-learning

repo: transformer-models


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function sPb_thin = spectralPb_fast(ws_wt2, nvec, ic_gamma, dthresh)
% function sPb_thin = spectralPb_fast(ws_wt2, nvec, ic_gamma, dthresh)
%
% description:
%   fast spectral gradient contours
%
% Jon Barron and Pablo Arbelaez
% <arbelaez@berkeley.edu>
% Jan 2014

if nargin<4, dthresh = 2; end
if nargin<3, ic_gamma = 0.12; end
if nargin<2, nvec = 6; end

[tx2, ty2] = size(ws_wt2);
tx=(tx2-1)/2; ty=(ty2-1)/2;

l{1} = ws_wt2(1:2:end,2:2:end);
l{2}= ws_wt2(2:2:end,1:2:end);

% build the pairwise affinity matrix
[val,I,J] = buildW(l{1},l{2}, dthresh, ic_gamma);
W = sparse(val,I,J);

[EigVect, EVal] =  ncuts_downsample3(W, nvec, 2, 2, [ty, tx]);

clear D W opts;

EigVal = diag(EVal);
clear Eval;

EigVal(1:end) = EigVal(end:-1:1);
EigVect(:, 1:end) = EigVect(:, end:-1:1);

vect = zeros(tx, ty, nvec);
for v = 2 : nvec,
    vect(:, :, v) = reshape(EigVect(:, v), [ty tx])';
end
clear EigVect;

%% spectral Pb
for v=2:nvec,
    vect(:,:,v)=(vect(:,:,v)-min(min(vect(:,:,v))))/(max(max(vect(:,:,v)))-min(min(vect(:,:,v))));
end

sPb_thin = zeros(2*tx+1, 2*ty+1);
for v = 1 : nvec
    if EigVal(v) > 0,
        vec = vect(:,:,v)/sqrt(EigVal(v));
        sPb_thin = sPb_thin + seg2bdry_wt(vec, 'doubleSize');
     end
end
sPb_thin = sPb_thin.^(1/sqrt(2));

%%

function [EV, EVal] = ncuts_downsample3(A, NVEC, N_DOWNSAMPLE, DECIMATE, SZ)
% A = affinity matrix
% NEVC = number of eigenvectors (set to 16?)
% N_DOWNSAMPLE = number of downsampling operations (2 seems okay)
% DECIMATE = amount of decimation for each downsampling operation (set to 2)
% SZ = size of the image corresponding to A

A_down = A;
SZ_down = SZ;

Bs = cell(N_DOWNSAMPLE,1);
for di = 1:N_DOWNSAMPLE

    % Create a binary array of the pixels that will remain after decimating
    % every other row and column
    [i,j] = ind2sub(SZ_down, 1:size(A_down,1));
    do_keep = (mod(i, DECIMATE) == 0) & (mod(j, DECIMATE) == 0);

    % Downsample the affinity matrix
    A_sub = A_down(:,do_keep)';

    % Normalize the downsampled affinity matrix
    d = (sum(A_sub,1) + eps);
    B = bsxfun(@rdivide, A_sub, d)';

    % "Square" the affinity matrix, while downsampling
    A_down = A_sub*B;
    SZ_down = floor(SZ_down / 2);

    % Hold onto the normalized affinity matrix for bookkeeping
    Bs{di} = B;
end

% Get the eigenvectors of the Laplacian
%EV = ncuts(A_down, NVEC);
[EV, EVal] = ncuts2(A_down, NVEC);

% diag(EVal)

% Upsample the eigenvectors
for di = N_DOWNSAMPLE:-1:1
    EV = Bs{di} * EV;
end

% whiten the eigenvectors, as they can get scaled weirdly during upsampling
EV = whiten(EV, 1, 0);

function [EV, EVal] = ncuts2(A, n_ev)
[wx, wy] = size(A);
x = 1 : wx;
S = full(sum(A, 1));
D = sparse(x, x, S, wx, wy);
clear S x;

opts.issym=1;
opts.isreal = 1;
opts.disp=0;
[EV, EVal] = eigs((D - A) + (10^-10) * speye(size(D)), D, n_ev, 'sm', opts);
clear D A opts;

v = diag(EVal);
[sorted, sortidx] = sort(v, 'descend');
EV = EV(:,sortidx);
EVal = diag(sorted);

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function [DEM] = spm_ADEM_update(DEM,COV)
% Updates ADEM structure using conditional expectations
% FORMAT [DEM] = spm_ADEM_update(DEM,COV)
%
% DEM - DEM structure
% COV - Covariance of parameter (P) fluctuations (E): P(i + 1) = P(i) + E
%     - where cov(E) = COV*pC
%
% This routine updates posterior expectations about states and parameters
% by replacing prior expectations with posterior expectations (and
% similarly updating hidden states and causes to the final iteration). If
% called with an extra argument, the posterior variances of the
% parameters are also updated.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_ADEM_update.m 6506 2015-07-24 10:26:51Z karl $

% preliminaries
%--------------------------------------------------------------------------
if nargin < 2, COV = 0; end

% update states and parameters (model)
%--------------------------------------------------------------------------
n     = length(DEM.M);
C     = DEM.qP.C;
for i = 1:(n - 1)

    % states
    %----------------------------------------------------------------------
    qE          = DEM.qU.x{i};
    if ~isempty(qE)
        DEM.M(i).x = spm_unvec(qE(:,end),DEM.M(i).x);
    end

    % parameters
    %----------------------------------------------------------------------
    qE          = spm_vec(DEM.qP.P{i});
    DEM.M(i).pE = spm_unvec(qE,DEM.M(i).pE);

    if nargin > 1

        % parameter covariance
        %------------------------------------------------------------------
        pC          = DEM.M(i).pC;
        np          = length(pC);
        qC          = C(1:np,1:np);
        C           = C(np + 1:end,np + 1:end);
        DEM.M(i).pC = qC + COV*pC;

    end

end

for i = 1:n
    if ~isempty(DEM.M(i).v)
        DEM.M(i).v  = spm_unvec(DEM.qU.v{i}(:,end),DEM.M(i).v);
    end
end

% update states and action (process)
%--------------------------------------------------------------------------
n     = length(DEM.G);
for i = 1:(n - 1)
    DEM.G(i).x = spm_unvec(DEM.pU.x{i}(:,end),DEM.G(i).x);
end
for i = 1:n
    if ~isempty(DEM.G(i).v)
        DEM.G(i).v  = spm_unvec(DEM.pU.v{i}(:,end),DEM.G(i).v);
    end
end
if isfield(DEM.G,'a')
    DEM.G(n).a = spm_unvec(DEM.qU.a{n}(:,end),DEM.G(n).a);
end

author: spm


text:
function phiFaceAverage = upwindMean1D(phi, u)
% This function gets the value of the field variable phi defined
% over the MeshStructure and calculates the upwind average on
% the cell faces, based on the direction of the velocity vector for a uniform mesh.
%
% SYNOPSIS:
%   phiFaceAverage = upwindMean1D(phi, u)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;

% check the size of the variable and the mesh dimension
Nx = phi.domain.dims(1);

% assign to a temp variable for boundary corrections
phi_tmp = phi.value;

% correct the value of phi at the boundary (calculation trick)
% assign the value of the left boundary to the left ghost cell
phi_tmp(1) = (phi.value(1)+phi.value(2))/2;
% assign the value of the right boundary to the right ghost cell
phi_tmp(end) = (phi.value(end)+phi.value(end-1))/2;

% calculate the average value
xvalue = (ux>0).*phi_tmp(1:Nx+1)+ ...
                        (ux<0).*phi_tmp(2:Nx+2)+ ...
                        0.5*(ux==0).*(phi.value(1:Nx+1)+phi.value(2:Nx+2));
phiFaceAverage=FaceVariable(phi.domain, xvalue, [], []);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%% Training face landmark detector
%
% This demo helps to train your own face landmark detector. You can train your
% own face landmark detection by just providing the paths for directory
% containing the images and files containing their corresponding face
% landmarks. As this landmark detector was originally trained on
% <http://www.ifp.illinois.edu/~vuongle2/helen/ HELEN dataset>, the training
% follows the format of data provided in HELEN dataset.
%
% The dataset consists of |.txt| files whose first line contains the image
% name which then follows the annotations. The format of a file containing
% annotations should be the following:
%
%  /data/helen/100032540_1.jpg
%  565.86 , 758.98
%  564.27 , 781.14
%  ...
%
% The above format is similar to HELEN dataset which is used for training the
% model.
%
% For a description of training parameters used in |configFile|, see the demo
% |facemark_kazemi_train_config_demo.m|.
%
% <<https://docs.opencv.org/3.4.0/2.jpg>>
%
% You can also download a pre-trained model |face_landmark_model.dat|,
% see the demo |facemark_kazemi_detect_img_demo|.
% (that way you can skip training and simply load the model).
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.4.0/modules/face/samples/sample_train_landmark_detector.cpp>
% * <https://github.com/opencv/opencv_contrib/blob/3.4.0/modules/face/tutorials/face_landmark/face_landmark_trainer.markdown>
%

%% Options

% [INPUT] path to the directory containing all text and image files
dname = fullfile(mexopencv.root(),'test','facemark','helen');
assert(isdir(dname), 'missing data directory');

% [INPUT] path to configuration xml file containing parameters for training
% https://github.com/opencv/opencv_contrib/raw/3.4.0/modules/face/samples/sample_config_file.xml
configFile = fullfile(mexopencv.root(),'test','facemark','config.xml');
assert(exist(configFile, 'file') == 2, 'missing train config file');

% [OUTPUT] path for saving the trained model
modelFile = fullfile(tempdir(), 'model_kazemi.dat');

% [INPUT] path to the cascade xml file for the face detector
xmlFace = fullfile(mexopencv.root(),'test','lbpcascade_frontalface.xml');
download_classifier_xml(xmlFace);

% name of user-defined face detector function
faceDetectFcn = 'myFaceDetector';
assert(exist([faceDetectFcn '.m'], 'file') == 2, 'missing face detect function');

% width/height which you want all images to get to scale the annotations.
% larger images are slower to process
scale = [460 460];

%% Data

% get names of files in which annotations and image names are found
filenames = cv.glob(fullfile(dname, '*.txt'));

% load image names and their corresponding landmarks
disp('Loading data...')
[imgFiles, pts] = cv.Facemark.loadTrainingData3(filenames);

% load images
imgs = cell(size(imgFiles));
for i=1:numel(imgFiles)
    if true
        % HELEN dataset annotations only store image basename
        fname = fullfile(dname, [imgFiles{i} '.jpg']);
    else
        fname = imgFiles{i};
    end
    imgs{i} = cv.imread(fname);
end

%% Init
% create instance of the face landmark detection class,
% and set the face detector function
obj = cv.FacemarkKazemi('ConfigFile',configFile);
obj.setFaceDetector(faceDetectFcn);

%% Train
% perform training
disp('Training...')
tic
success = obj.training(imgs, pts, configFile, scale, 'ModelFilename',modelFile);
toc
if success
    disp('Training successful')
else
    disp('Training failed')
end

%%
% In the above call, |scale| is passed to scale all images and their
% corresponding landmarks, as it takes greater time to process large images.
% After scaling data it calculates mean shape of the data which is used as
% initial shape while training. It trains the model and stores the trained
% model file with the specified filename. As the training starts, you will see
% something like this:
%
% <<https://docs.opencv.org/3.4.0/train1.png>>
%
% The error rate on trained images depends on the number of images used for
% training:
%
% <<https://docs.opencv.org/3.4.0/train.png>>
%
% The error rate on test images depends on the number of images used for
% training:
%
% <<https://docs.opencv.org/3.4.0/test.png>>
%

%% Helper functions

function download_classifier_xml(fname)
    if exist(fname, 'file') ~= 2
        % attempt to download trained Haar/LBP/HOG classifier from Github
        url = 'https://cdn.rawgit.com/opencv/opencv/3.4.0/data/';
        [~, f, ext] = fileparts(fname);
        if strncmpi(f, 'haarcascade_', length('haarcascade_'))
            url = [url, 'haarcascades/'];
        elseif strncmpi(f, 'lbpcascade_', length('lbpcascade_'))
            url = [url, 'lbpcascades/'];
        elseif strncmpi(f, 'hogcascade_', length('hogcascade_'))
            url = [url, 'hogcascades/'];
        else
            error('File not found');
        end
        urlwrite([url f ext], fname);
    end
end

% The facemark API provides the functionality to the user to use their own
% face detector. The code below implements a sample face detector. This
% function must be saved in its own M-function to be used by the facemark API.
function faces = myFaceDetector(img)
    persistent obj
    if isempty(obj)
        obj = cv.CascadeClassifier();
        obj.load(xmlFace);
    end

    if size(img,3) > 1
        gray = cv.cvtColor(img, 'RGB2GRAY');
    else
        gray = img;
    end
    gray = cv.equalizeHist(gray);
    faces = obj.detect(gray, 'ScaleFactor',1.4, 'MinNeighbors',2, ...
        'ScaleImage',true, 'MinSize',[30 30]);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function varargout = or(varargin)
%OR (overloaded)

% Models OR using a nonlinear operator definition
if nargin == 1
    varargout{1} = varargin{1};
else
    varargout{1} = (yalmip('define','lmior',varargin{:}) == 1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% The COBRAToolbox: testAddKeyToKnownHosts.m
%
% Purpose:
%     - test if a key can be added to known hosts
%
% Authors:
%     - Laurent Heirendt, February 2018
%

global CBTDIR

% save the current path
currentDir = pwd;

% test for a failure of an unknown host
assert(verifyCobraFunctionError('addKeyToKnownHostsCT','inputs',{'github123.co'}));

% find the keyscan
[status_keyscan, result_keyscan] = system('ssh-keyscan');

if status_keyscan == 1 && contains(result_keyscan, 'usage:')
    % get the user directory
    if ispc
        homeDir = getenv('userprofile');
    else
        homeDir = getenv('HOME');
    end

    % define the name of the file
    khFile = [homeDir, filesep, '.ssh', filesep, 'known_hosts'];
    khFile_bk = [khFile, '_bk'];

    % backup the existing known hosts file if it exists
    [status, ~, ~] = copyfile(khFile, khFile_bk);

    % remove the key first from the known hosts file
    [status, result] = system('ssh-keygen -R github.com');

    % run the function to add the github.com key
    statusAddKey = addKeyToKnownHostsCT();

    % check if the key has been added succesfully
    assert(statusAddKey);

    % verify that the key has been added to the known hosts file
    [~, result_grep] = system(['grep "^github.com " ', khFile]);

    % test whether the resulting string is not empty (site has been found)
    assert(~strcmp(result_grep, ''))

    % test if the host is already known
    assert(addKeyToKnownHostsCT('github.com'));

    % remove the old host file generated by running ssh-keygen -R
    delete([khFile '.old']);

    % remove the changed file and restore the backup file
    delete(khFile);
    [status, ~, ~] = movefile(khFile_bk, khFile);
    if status
        fprintf(' > Original known hosts file restored.\n');
    end

    % print a success message
    fprintf(' > Test to add a key passed.\n');
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
disp('Hello world!')

name = input('whats your name? ','s');

disp(['hello ',name])

repo: Aero-matlab