text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

[max_value, p] = max(sigmoid(X * all_theta'), [], 2);

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

end

author: gopaczewski

repo: coursera-ml


text:
% Extract the coefficient of a monomial t in the polynomial p
%
% Syntax:   >> c = PolynomialCoefficient(p,t)
% Syntax:   >> c = pcoef(p,t)
%
%   Input:    p --- (string) polynomial
%             t --- (string) monomial
%
%  Output:    c --- (numeric or string) the coefficient of t in p
%
% Example:   >> p = '-2*x + 9*x*y^2 - 9*x^2*y^2 + 8*x^3*y^2';
%            >> c = Coeff(p,'y^2')
%
%            c =
%
%            9*x - 9*x^2 + 8*x^3
%
%            >> d = Coeff(p,'x*y^2')
%
%            d =
%
%            9

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function x = PQgetData (WAV, i, N)
% Get data from internal buffer or file
% i - file position
% N - number of samples
% x - output data (scaled to the range -32768 to +32767)

% Only two files can be "active" at a time.
% N = 0 resets the buffer

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:34:10 $

persistent Buff

iB = WAV.iB + 1;
if (N == 0)
    Buff(iB).N = 20 * 1024;     % Fixed size
    Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N);
    Buff(iB).i = i;
end

if (N > Buff(iB).N)
    error ('>>> PQgetData: Request exceeds buffer size');
end

% Check if requested data is not already in the buffer
is = i - Buff(iB).i;
if (is < 0 | is + N - 1 > Buff(iB).N - 1)
    Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N);
    Buff(iB).i = i;
end

% Copy the data
Nchan = WAV.Nchan;
is = i - Buff(iB).i;
x = Buff(iB).x(1:Nchan,is+1:is+N-1+1);

%------
function x = PQ_ReadWAV (WAV, i, N)
% This function considers the data to extended with zeros before and
% after the data in the file. If the starting offset i is negative,
% zeros are filled in before the data starts at offset 0. If the request
% extends beyond the end of data in the file, zeros are appended.

Amax = 32768;
Nchan = WAV.Nchan;

x = zeros (Nchan, N);

Nz = 0;
if (i < 0)
    Nz = min (-i, N);
    i = i + Nz;
end

Ns = min (N - Nz, WAV.Nframe - i);
if (i >= 0 & Ns > 0)
    x(1:Nchan,Nz+1:Nz+Ns-1+1) = Amax * (wavread (WAV.Fname, [i+1 i+Ns-1+1]))';
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [ pn, dist ] = quad_point_near_2d ( q, p )

%*****************************************************************************80
%
%% QUAD_POINT_NEAR_2D computes the nearest point on a quadrilateral in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(2,4), the quadrilateral vertices.
%
%    Input, real P(2), the point whose nearest quadrilateral point
%    is to be determined.
%
%    Output, real PN(2), the nearest point to P.
%
%    Output, real DIST, the distance from the point to the
%    quadrilateral.
%
  dim_num = 2;
  nside = 4;
%
%  Find the distance to each of the line segments that make up the edges
%  of the quadrilateral.
%
  dist = Inf;
  pn(1:dim_num) = 0.0;

  for j = 1 : nside

    jp1 = i4_wrap ( j+1, 1, nside );

    [ pn2, dist2, tval ] = segment_point_near_2d ( q(1:dim_num,j)', ...
      q(1:dim_num,jp1)', p );

    if ( dist2 < dist )
      dist = dist2;
      pn(1:dim_num) = pn2(1:dim_num);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [index,numBeforeDeg]=monomial2DegNegLexOrderIdx(monomial,numBeforeDeg)
%%MONOMIAL2DEGNEGLEXORDERIDX Given the nonegative integer exponents of a
%          monomial term (e.g. x^5*y^3*z^6 is [5;3;6]), determine the index
%          of that monomial term in a list of all monomial terms in degree
%          negative lexicographic ordering. The degree of a monomial is the
%          sum of the exponents of the terms.
%
%INPUTS: monomial A numDimXnumMonomials set of numMonomials monomial terms
%                 consisting of nonnegative integer vectors of exponents of
%                 the elements Alternatively, if one just wants the
%                 numBeforeDegVec vector, then this should be an empty
%                 matrix and numBeforeDegVec should contain the
%                 dimensionality and the maximum degree desired.
%    numBeforeDeg An optional parameter. Providing this parameter will
%                 speed up repeated calls to this function by avoiding the
%                 recalculation of certain values. This vector is such that
%                 numBeforeDeg(i+1) is the number of monomial terms
%                 before degree i starting with degree 0, where
%                 numBeforeDeg(0+1)=0. This parameter is returned by the
%                 function. Note that numBeforeDeg depends on numDim, so
%                 vectors cannot be reused between cals where numDim
%                 changes. If one just wants this parameter (for example
%                 in advance to speed up many repeated calls), then set
%                 monomial=[] and numBeforeDeg=[numDim;maxDeg], where
%                 maxDeg the maximum degree desired, >=0.
%
%OUTPUTS: index A numMonomialsX1 vector of the position of the monomials in
%              the ordering of monomials in degree negative lexicographic
%              order, starting from 1. If monomial=[] and numBeforeDeg
%              is set to the dimensionality and maximum degree, then this
%              is an empty matrix.
% numBeforeDeg The same vector as the optional input. If the input is
%              not provided, then numBeforeDeg will be filled with
%              values up to the degree of the monomial provided having the
%              highest degree. If monomials=[] and numBeforeDeg on the
%              input is set to the maximum degree desired, then
%              numBeforeDeg will be filled all the way up to that
%              degree.
%
%Various monomial orderings are used in algorithms that manipulate
%polynomials, such as methods for solving or minimizing multivariate
%polynomials as in [1] and [2].
%
%The degree negative lexicographic ordering, for a fixed degree, coincides
%with the ordering of the terms provided by unrankTComposition with
%firstElMostSig=true. However, that rank is only an offset for a given
%degree. The number of terms before that degree is the sum of the number of
%compositions of lower degrees. Thus, one sees the need for the
%numBeforeDeg vector. The number of compositions of a given degree into
%n parts, (each part starting from 0, not 1) is
%binomial(degree+numDim-1,numDim-1).
%
%EXAMPLE:
%Consider monomials like x^1*y^0*z^0, x^12*y^4*z^8;, x*y*z, x^0*y^0*z^0(=1)
% monomial=[1,12,1,0;
%           0, 4,1,0;
%           0, 8,1,0];
% [index,numBeforeDeg]=monomial2DegNegLexOrderIdx(monomial)
% %One will get the indices index=[2;2687;15;1].
% %The value numBeforeDeg can be reused with monomials having the same
% %number of variables, such as
% monomial=[0;1;1];
% index=monomial2DegNegLexOrderIdx(monomial,numBeforeDeg)
% %On the other hand, if one only wants the numBeforeDeg vector for up to
% %degree 24 for three variables, then one
% [~,numBeforeDeg]=monomial2DegNegLexOrderIdx([],[3,24])
%
%REFERENCES:
%[1] P. Dreesen, "Back to the roots: Polynomial system solving using linear
%    algebra," Ph.D. dissertation, Katholieke Universiteit Leuven, Leuven,
%    Flanders, Belgium, Sep. 2013.
%[2] A. Papachristodoulou and S. Prajna, "A tutorial on sum of squares
%    techniques for systems analysis," in Proceedings of the American
%    Control Conference, Portland, OR, 8-10 Jun. 2005, pp. 2686-2700.
%
%February 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(monomial))
    numDim=numBeforeDeg(1);
    maxDeg=numBeforeDeg(2);
else
    numDim=size(monomial,1);
    degs=sum(monomial,1);
    maxDeg=max(degs);
end

%Fill in the numBeforeDegVec vector.
if(isempty(monomial)||nargin<2||isempty(numBeforeDeg))
    numBeforeDeg=zeros(maxDeg+1,1);

    if(maxDeg>0)
        numBeforeDeg(1+1)=1;
    end

    for curDeg=2:maxDeg
        numBeforeDeg(curDeg+1)=numBeforeDeg(curDeg-1+1)+binomial(curDeg-1+numDim-1,numDim-1);
    end

    %If just the list numBeforeDegVec is desired.
    if(isempty(monomial))
        index=[];
        return;
    end
end

numMonomials=size(monomial,2);
%Allocate space.
index=zeros(1,numMonomials);

for curMonomial=1:numMonomials
    offset=numBeforeDeg(degs(curMonomial)+1);
    index(curMonomial)=rankTComposition(monomial(:,curMonomial)+1,true)+offset+1;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [L,W] = get_l(n,d)
%GET_L Compute discrete derivative operators.
%
% [L,W] = get_l(n,d)
%
% Computes the discrete approximation L to the derivative operator
% of order d on a regular grid with n points, i.e. L is (n-d)-by-n.
%
% L is stored as a sparse matrix.
%
% Also computes W, an orthonormal basis for the null space of L.

% Per Christian Hansen, IMM, 02/05/98.

% Initialization.
if (d<0), error ('Order d must be nonnegative'), end

% Zero'th derivative.
if (d==0), L = speye(n); W = zeros(n,0); return, end

% Compute L.
c = [-1,1,zeros(1,d-1)];
nd = n-d;
for i=2:d, c = [0,c(1:d)] - [c(1:d),0]; end
L = sparse(nd,n);
for i=1:d+1
  L = L + sparse(1:nd,(1:nd)+i-1,c(i)*ones(1,nd),nd,n);
end

% If required, compute the null vectors W via modified Gram-Schmidt.
if (nargout==2)
  W = zeros(n,d);
  W(:,1) = ones(n,1);
  for i=2:d, W(:,i) = W(:,i-1).*(1:n)'; end
  for k=1:d
     W(:,k) = W(:,k)/norm(W(:,k));
     W(:,k+1:d) = W(:,k+1:d) - W(:,k)*(W(:,k)'*W(:,k+1:d));
  end
end

repo: REST


text:
function ierror = perm0_check ( n, p )

%*****************************************************************************80
%
%% PERM0_CHECK checks a permutation of ( 0, ..., N-1).
%
%  Discussion:
%
%    The routine verifies that each of the integers from 0 to
%    to N-1 occurs among the N entries of the permutation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer P(N), the array to check.
%
%    Output, integer IERROR:
%    0, P is a legal permutation of (0,...,N-1).
%    1, P is not a legal permuation of ( 0,...,N-1).
%
  ierror = 0;

  for value = 0 : n - 1

    ierror = 1;

    for location = 1 : n
      if ( p(location) == value )
        ierror = 0;
        break;
      end
    end

    if ( ierror ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PERM0_CHECK - Warning!\n' );
      fprintf ( 1, '  Permutation is missing the value %d.\n', value );
      break;
    end

  end

  return
end


text:
%UT_WEIGHTS - Generate unscented transformation weights
%
% Syntax:
%   [WM,WC,c] = ut_weights(n,alpha,beta,kappa)
%
% In:
%   n     - Dimensionality of random variable
%   alpha - Transformation parameter  (optional, default 0.5)
%   beta  - Transformation parameter  (optional, default 2)
%   kappa - Transformation parameter  (optional, default 3-n)
%
% Out:
%   WM - Weights for mean calculation
%   WC - Weights for covariance calculation
%    c - Scaling constant
%
% Description:
%   Computes unscented transformation weights.
%
% See also UT_MWEIGHTS UT_TRANSFORM UT_SIGMAS
%

% Copyright (C) 2006 Simo S�rkk�
%
% $Id$
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [WM,WC,c] = ut_weights(n,alpha,beta,kappa)

%
% Check which arguments are there
%
if nargin < 1
  error('At least dimensionality n required.');
end
if nargin < 2
  alpha = [];
end
if nargin < 3
  beta = [];
end
if nargin < 4
  kappa = [];
end

%
% Apply default values
%
if isempty(alpha)
  alpha = 1;
end
if isempty(beta)
  beta = 0;
end
if isempty(kappa)
  kappa = 3 - n;
end

%
% Compute the normal weights
%
lambda = alpha^2 * (n + kappa) - n;

WM = zeros(2*n+1,1);
WC = zeros(2*n+1,1);
for j=1:2*n+1
  if j==1
    wm = lambda / (n + lambda);
    wc = lambda / (n + lambda) + (1 - alpha^2 + beta);
  else
    wm = 1 / (2 * (n + lambda));
    wc = wm;
  end
  WM(j) = wm;
  WC(j) = wc;
end

c = n + lambda;

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
filename = 'CantileverSquareSym';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
% m1 = 0.0101;
% m2 = 0.0101;
m1 = 0.8;
m2 = 0.8;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
%vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 50;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 500;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function p_xt= Motion_Model_Velocity_Test(u, x_p, a, gridsize)
%% Test to figure out which 'a' parameters suit you best
    % u: control command
    % x_p: initial position with orientation
    % a: [a_1 ... a_6 ] of tested parameters

    %the test will calculate the probability of the robot moving to every
    %possible cell in the grid generating a banana-shaped density function
    p_xt=zeros(gridsize(1), gridsize(2));
    for i=1:gridsize(1)
        for j=1:gridsize(2)-1
            p_xt(i,j) = Motion_Model_Velocity([j i x_p(3)], u, x_p, a) ;
        end
    end

    % Visualize results
    figure(66)
    clf;
    hold on;
    plot(x_p(1), x_p(2), 'dr');
    surf(p_xt);
    view([0 -90]);
    hold off;

end

author: Allopart

repo: rbpf-gmapping


text:
% VL_TEST_AIBHIST
function vl_test_aibhist

D   = 4 ;
K   = 20 ;

randn('state',0) ;
rand('state',0) ;

X1 = randn(2,300) ; X1(1,:) = X1(1,:) + 2 ;
X2 = randn(2,300) ; X2(1,:) = X2(1,:) - 2 ;
X3 = randn(2,300) ; X3(2,:) = X3(2,:) + 2 ;

C   = 1:K*K ;
Pcx = zeros(3,K*K) ;

f1 = quantize(X1,D,K) ;
f2 = quantize(X2,D,K) ;
f3 = quantize(X3,D,K) ;

Pcx(1,:) = vl_binsum(Pcx(1,:), ones(size(f1)), f1) ;
Pcx(2,:) = vl_binsum(Pcx(2,:), ones(size(f2)), f2) ;
Pcx(3,:) = vl_binsum(Pcx(3,:), ones(size(f3)), f3) ;

Pcx = Pcx / sum(Pcx(:)) ;

[parents_, cost_] = vl_aib(Pcx) ;
[parents,  cost ] = vl_aib(Pcx,'clusternull') ;

% find a null node for testing purposes
anull = min(find(parents_==0)) ;
f1 = [f1 repmat(anull,1,10)] ;

figure(100);  clf ;
subplot(1,2,1) ; hold on ;
plot(parents_, 'r.-') ;
plot(parents,  'g') ;
legend('signal null', 'cluster null') ;
subplot(1,2,2) ; hold on ;
plot(cost_, 'r.-') ;
plot(cost,  'g') ;
legend('signal null', 'cluster null') ;

range = [1 10 K*K-10 K*K] ;
for c=1:length(range)

  cut_size = range(c) ;

  % compare two methods of getting the same cut histogram
  [cut_,map_]  = vl_aibcut(parents_, cut_size) ;
  hist_        = vl_aibcuthist(map_, f1, 'nulls', 'append') ;
  histtree_    = vl_aibhist(parents_, f1) ;
  thist_       = histtree_(cut_) ;

  [cut,map]    = vl_aibcut(parents, cut_size) ;
  hist         = vl_aibcuthist(map, f1, 'nulls', 'append') ;
  histtree     = vl_aibhist(parents, f1) ;
  thist        = histtree(cut) ;

  figure(100 + c) ; clf ;
  subplot(2,2,1) ; hold on ; plot(hist_,'g.-') ; plot(thist_,'r') ;
  legend('cut+cuthist', 'hist+cut') ;
  title('vl_aibcuthist vs aibhist');
  subplot(2,2,2) ; hold on ; plot(histtree_) ; title('aibtree') ;

  subplot(2,2,3) ; hold on ; plot(hist,'g.-') ; plot(thist,'r') ;
  legend('cut+cuthist', 'hist+cut') ;
  title('vl_aibcuthist vs vl_aibhist (clust null)');
  subplot(2,2,4) ; hold on ; plot(histtree) ; title('aibtree (clust null)') ;
end

% --------------------------------------------------------------------
function f = quantize(X,D,K)
% --------------------------------------------------------------------
d = 2*D / K ;
j = round((X(1,:) + D) / d) ;
i = round((X(2,:) + D) / d) ;
j = max(min(j,K),1) ;
i = max(min(i,K),1) ;
f = sub2ind([K K],i,j) ;

author: jianxiongxiao

repo: ProfXkit


text:
function Write_ContinuousDynamics(FileWritingSetup)
% Write_ContinuousDynamics(FileWritingSetup)
%
% FUNCTION:
%
%   This function takes the symbolically derived equations of motion for
%   the retractable double pendulum biped model and writes them to function
%   files.
%
%
% INPUTS: FileWritingSetup.
%
%   States = [2N x 2] cell array where the first column is the name of
%       the states, and the second column is comments about that states.
%       The index of a state in this array will match the index in the
%       dynamics function to be written. This function assumes a second
%       order system that is being written in first order form, so if
%       States{i,1} = z, then States{i+N,1} = dz; Furthermore, the output
%       (dStates) will have the form: dStates(i,:) = dz,
%       and dStates(i+N,:) = ddz.
%
%   Actuators = [P x 2] cell array giving the names of each of the
%       actuators and their comments.
%
%   Contacts = [M x 2] cell array giving all of the contact forces to
%       output. The first column is names, second column is comments
%
%   Parameters = [K x 2] cell array with all of the input parameters. The
%       first column is for names, the second for comments
%
%   FileData = [L x 3+?] cell array with the names of each field in Dyn in
%       the first column, and the file name to write those dynamics
%       equations to. The final column is used for comments that are
%       specific to those dynamics. Seperate multiple lines by commas:
%       'comment line one', 'comment line two', ...
%
%   Dyn = A struct, with a field for each dynamics equation to write. Every
%       field must contain a symbolic expression for each of the state
%       accelerations and contact forces.
%
%   Directory = the name of the directory to write the files in.
%
%
% OUTPUTS:
%   This function will write a file for the dynamics equations in each
%   field of Dyn.
%
%
% NOTE:
%
%   1)  It is assumed that the symbols in all of the inputs are self
%       consistent; that is - if Contacts{1,1} = H1, then there must be a field
%       in Dyn.(FileNames{i,1}) called H1.
%
%
% Written by Matthew Kelly
% December 6, 2013
% Cornell University
%
% See also DERIVE_EOM

States = FileWritingSetup.States;
Dyn = FileWritingSetup.Dyn;
Contacts = FileWritingSetup.Contacts;
Actuators = FileWritingSetup.Actuators;
Parameters = FileWritingSetup.Parameters;
FileData = FileWritingSetup.FileData;
Directory = FileWritingSetup.Directory;

CurrDir = cd;
if ~exist(Directory,'dir')
    mkdir(Directory);
end
cd(Directory);

for idxPhase = 1:size(FileData,1)
    fid = fopen([FileData{idxPhase,2} '.m'],'w');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                         Function Header                             %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    fprintf(fid,['function [dStates, contactForces] = ' ...
        FileData{idxPhase,2} '(States, Actuators, Parameters)\n']);
    fprintf(fid,['%% function [dStates, Actuators] = ' ...
        FileData{idxPhase,2} '(States, Actuators, Parameters)\n']);
    fprintf(fid,'%%\n%% Computer Generated File -- DO NOT EDIT \n');
    fprintf(fid,['%%\n%% This function was created by the function '...
        'Write_ContinuousDynamics()\n']);
    fprintf(fid,['%% ' datestr(now) '\n%%\n']);
    for idxComment = 3:length(FileData(idxPhase,:))
        fprintf(fid,['%% ' FileData{idxPhase,idxComment} '\n']);
    end %idxComment
    fprintf(fid,'%%\n');
    fprintf(fid,'%% Matthew Kelly \n%% Cornell University \n%% \n\n');
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                           Parse Inputs                              %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    for idxState=1:size(States,1)
        fprintf(fid,[States{idxState,1} ' = States(' num2str(idxState) ',:); %% ' ...
            States{idxState,2} '\n']);
    end %idxState
    fprintf(fid,'\n');

    for idxAct=1:size(Actuators,1)
        rhs = [' = Actuators(' num2str(idxAct) ',:); %% ' Actuators{idxAct,2} '\n'];
        %Need to make sure that we disconnect the ankle motors when their
        %respective feed are not in contact with the ground!
        if strcmp(FileData{idxPhase,1},'Flight') || strcmp(FileData{idxPhase,1},'SingleTwo')
            if strcmp(Actuators{idxAct,1},'T1')
                rhs = ' = 0; %%Foot One not in contact with ground!\n';
            end
        end
        if strcmp(FileData{idxPhase,1},'Flight') || strcmp(FileData{idxPhase,1},'SingleOne')
            if strcmp(Actuators{idxAct,1},'T2')
                rhs = ' = 0; %%Foot Two not in contact with ground!\n';
            end
        end
        fprintf(fid,[Actuators{idxAct,1} rhs]);
    end %idxAct
    fprintf(fid,'\n');

    for idxParam=1:size(Parameters,1)
        fprintf(fid,[Parameters{idxParam,1} ' = Parameters.' ...
            Parameters{idxParam,1} '; %% ' ...
            Parameters{idxParam,2} '\n']);
    end %idxState
    fprintf(fid,'\n');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                   Explicitly write constraints                      %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    Ndof = size(States,1)/2;
    fprintf(fid,'%% Constraints for this phase: \n');
    for i=1:Ndof
        writeIndex = i + Ndof;
        writeName = ['d' States{writeIndex,1}];
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        if (writeData==0)
            fprintf(fid,[writeName ' = 0; %%' States{i,2} '\n']);
        end
    end
    for i=1:size(Contacts,1)
        writeName = Contacts{i,1};
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        if (writeData==0)
            fprintf(fid,[writeName ' = 0; %%' Contacts{i,2} '\n']);
        end
    end
    fprintf(fid,'\n');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                       Write Accelerations                           %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

     N = num2str(Ndof);
    fprintf(fid,'dStates = zeros(size(States));\n');
    fprintf(fid,['dStates(1:' N ',:) = '...
        'States((1+' N '):(' N '+' N '),:);\n']);
    for idxAcc=1:Ndof
        writeIndex = idxAcc + Ndof;
        writeName = ['d' States{writeIndex,1}];
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        vectorizedExpression = vectorize(writeData);
        fprintf(fid,['dStates(' num2str(writeIndex) ',:) = ' ...
            vectorizedExpression ';\n']);
    end %idxAcc
    fprintf(fid,'\n');

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
    %                       Write Contact Forces                          %
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    for idxContact=1:size(Contacts,1)
        fprintf(fid,['%% contactForces(' num2str(idxContact) ',:) == '...
            Contacts{idxContact,1} ' == ' Contacts{idxContact,2} '\n']);
    end %idxContact
    fprintf(fid,['contactForces = zeros(' num2str(size(Contacts,1))...
        ',size(States,2));\n']);
    for idxContact=1:size(Contacts,1)
        writeIndex = num2str(idxContact);
        writeName = Contacts{idxContact,1};
        writeData = Dyn.(FileData{idxPhase,1}).(writeName);
        vectorizedExpression = vectorize(writeData);
        fprintf(fid,['contactForces(' writeIndex ',:) = ' ...
            vectorizedExpression ';\n']);
    end %idxContact
    fprintf(fid,'\n');

    fprintf(fid,'end\n');
    fclose(fid);
end %idxPhase
cd(CurrDir);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  ~~~~~~~~~~~~~~~~~~~ Sub Functions  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function characterExpression = vectorize(symbolicExpression)

characterExpression = char(symbolicExpression);
characterExpression = strrep(characterExpression,'*','.*');
characterExpression = strrep(characterExpression,'/','./');
characterExpression = strrep(characterExpression,'^','.^');
characterExpression = strrep(characterExpression,...
    '0','zeros(1,size(States,2))');

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function varargout = opensignals_txt(filename, hdr, begsample, endsample, chanindx)

% OPENSIGNALS_TXT reads time series data from a Bitalino OpenSignals txt file
%
% Use as
%   hdr = opensignals_txt(filename);
%   dat = opensignals_txt(filename, hdr, begsample, endsample, chanindx);
%   evt = opensignals_txt(filename, hdr);
%
% See also FT_FILETYPE, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, QUALISYS_TSV, MOTION_C3D

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

needhdr = (nargin==1);
needevt = (nargin==2);
needdat = (nargin==5);

% use the full filename including path to distinguish between similarly named files in different directories
[p, f, x] = fileparts(filename);
if isempty(p)
  % no path was specified
  fullname = which(filename);
elseif startsWith(p, ['.' filesep])
  % a relative path was specified
  fullname = fullfile(pwd, p(3:end), [f, x]);
else
  fullname = filename;
end

fid = fopen_or_error(fullname, 'r');
line1 = fgetl(fid);
line2 = fgetl(fid);
fclose(fid);
% convert the json string into a structure
orig = jsondecode(line2(2:end));
fn = fieldnames(orig);
if numel(fn)==1
  % it contains a single field with the date and time
  orig = orig.(fn{1});
end

dat = readtable(fullname, 'FileType', 'text');
keep = false(1,size(dat,2));
for i=1:size(dat,2)
  keep(i) = isnumeric(dat{:,i});
end
dat = table2array(dat(:,keep))';

if needhdr
  %% parse the header
  hdr.Fs = orig.samplingRate;
  hdr.nChans = length(orig.column);
  hdr.nTrials = 1;
  hdr.nSamplesPre = 0;
  hdr.nSamples = size(dat,2);
  hdr.label = orig.column;

  % remember the details
  hdr.orig = orig;

  % return the header details
  varargout = {hdr};

elseif needdat
  %% parse the data
  dat = dat(chanindx, begsample:endsample);

  % return the data
  varargout = {dat};

elseif needevt
  %% parse the events
  ft_warning('cannot read events from %s', filename);

  varargout = {[]};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Generates state evolution and compares to observed MSEs for D-AMP and D-IT algorithms

addpath(genpath('..'));

%Parameters
denoiser='BM3D';
SamplingRate=.2;
noise_sig=0;
imsize=128;
filename='barbara.png';
iters=30;
N=5;%Number of tests to run.  Must be at least 2

Im=double(imread(filename));
x_0=imresize(Im,imsize/size(Im,1));
[height, width]=size(x_0);
n=length(x_0(:));
m=round(n*SamplingRate);

%Generate State Evolution and compute intermediate MSEs of D-AMP and D-IT
True_DAMP_MSE_array=zeros(N,iters);
True_DIT_MSE_array=zeros(N,iters);
Predicted_MSE_array=zeros(N,iters);
True_DAMP_MSE_array(:,1)=mean(x_0(:).^2);
True_DIT_MSE_array(:,1)=mean(x_0(:).^2);
Predicted_MSE_array(:,1)=mean(x_0(:).^2);
for i=1:N
    i
    M=randn(m,n);
    for j = 1:n
        M(:,j) = M(:,j) ./ sqrt(sum(abs(M(:,j)).^2));
    end
    y=M*x_0(:)+noise_sig*randn(m,1);
    x_t=zeros(n,1);x_t2=zeros(n,1);
    z_t=y;
    z_t2=y;
    for iter=2:iters
            [x_tplus1,z_tplus1,NA] = DAMP_oneIter(y,z_t,x_t,width,height,denoiser,M);
            z_t=z_tplus1; x_t=x_tplus1;
            [x_tplus12,z_tplus12,NA] = DIT_oneIter(y,z_t2,x_t2,width,height,denoiser,M);
            z_t2=z_tplus12; x_t2=x_tplus12;
            True_DAMP_MSE_array(i,iter)=mean((x_0(:)-x_tplus1(:)).^2);
            True_DIT_MSE_array(i,iter)=mean((x_0(:)-x_tplus12(:)).^2);
            Predicted_MSE_array(i,iter)= DAMP_SE_Prediction(x_0(:), Predicted_MSE_array(i,iter-1), m,n,noise_sig,denoiser,width,height);
    end
end
True_DAMP_MSE=mean(True_DAMP_MSE_array)';
True_DAMP_MSE_std=std(True_DAMP_MSE_array)';
True_DIT_MSE=mean(True_DIT_MSE_array)';
True_DIT_MSE_std=std(True_DIT_MSE_array)';
Predicted_MSE=mean(Predicted_MSE_array)';
Predicted_MSE_std=std(Predicted_MSE_array)';

%Plot Results
h=figure;
hold
errorbar(0:29,Predicted_MSE,Predicted_MSE_std,'-.b');
errorbar(0:29,True_DAMP_MSE,True_DAMP_MSE_std,'--g');
errorbar(0:29,True_DIT_MSE,True_DIT_MSE_std,'-r');
title([denoiser,'-AMP and ', denoiser '-IT State Evolution']);
legend(['Predicted ',num2str(100*SamplingRate),'%'], [denoiser,'-AMP ', num2str(100*SamplingRate),'%'], [denoiser,'-IT ', num2str(100*SamplingRate),'%']);
xlabel('Iteration');
ylabel('MSE');

author: ricedsp


text:
classdef StressVoigt3DRotator < StressVoigtRotator

    methods (Access = public)

        function obj = StressVoigt3DRotator(angle,dir)
            obj.compute(angle,dir)
        end
    end

    methods (Access = protected)

        function generateRotator(obj)
            theta = obj.angle;
            d  = obj.dir.getValue();
            u1 = d(1);
            u2 = d(2);
            u3 = d(3);

            A(1,1) = ((1 - cos(theta))*u1^2 + cos(theta))^2;
            A(2,2) = ((1 - cos(theta))*u2^2 + cos(theta))^2;
            A(3,3) = ((1 - cos(theta))*u3^2 + cos(theta))^2;

            A(4,4) = ((1 - cos(theta))*u2^2 + cos(theta))*((1 - cos(theta))*u3^2 + cos(theta)) - (u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u1*sin(theta) - u2*u3*(cos(theta) - 1));
            A(5,5) = ((1 - cos(theta))*u1^2 + cos(theta))*((1 - cos(theta))*u3^2 + cos(theta)) - (u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1));
            A(6,6) = ((1 - cos(theta))*u1^2 + cos(theta))*((1 - cos(theta))*u2^2 + cos(theta)) - (u3*sin(theta) + u1*u2*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1));

            A(1,2) = (u3*sin(theta) - u1*u2*(cos(theta) - 1))^2;
            A(2,1) = (u3*sin(theta) + u1*u2*(cos(theta) - 1))^2;

            A(1,3) = (u2*sin(theta) + u1*u3*(cos(theta) - 1))^2;
            A(3,1) = (u2*sin(theta) - u1*u3*(cos(theta) - 1))^2;

            A(1,4) = -2*(u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1));
            A(4,1) = -(u3*sin(theta) + u1*u2*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1));

            A(1,5) = -2*(u2*sin(theta) + u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));
            A(5,1) = (u2*sin(theta) - u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));

            A(1,6) = 2*(u3*sin(theta) - u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));
            A(6,1) = -(u3*sin(theta) + u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));

            A(2,3) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))^2;
            A(3,2) = (u1*sin(theta) + u2*u3*(cos(theta) - 1))^2;

            A(2,4) = 2*(u1*sin(theta) - u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));
            A(4,2) = -(u1*sin(theta) + u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));

            A(2,5) = -2*(u3*sin(theta) + u1*u2*(cos(theta) - 1))*(u1*sin(theta) - u2*u3*(cos(theta) - 1));
            A(5,2) = -(u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1));

            A(2,6) = -2*(u3*sin(theta) + u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));...
                A(6,2) = (u3*sin(theta) - u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));

            A(3,4) = -2*(u1*sin(theta) + u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));
            A(4,3) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));

            A(3,5) = 2*(u2*sin(theta) - u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));
            A(5,3) = -(u2*sin(theta) + u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));

            A(3,6) = -2*(u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1));
            A(6,3) = -(u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u1*sin(theta) - u2*u3*(cos(theta) - 1));

            A(4,5) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))*(u2*sin(theta) - u1*u3*(cos(theta) - 1)) - (u3*sin(theta) + u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));
            A(5,4) = (u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u2*sin(theta) + u1*u3*(cos(theta) - 1)) + (u3*sin(theta) - u1*u2*(cos(theta) - 1))*((1 - cos(theta))*u3^2 + cos(theta));

            A(4,6) = (u1*sin(theta) + u2*u3*(cos(theta) - 1))*(u3*sin(theta) + u1*u2*(cos(theta) - 1)) + (u2*sin(theta) - u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));
            A(6,4) = (u1*sin(theta) - u2*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1)) - (u2*sin(theta) + u1*u3*(cos(theta) - 1))*((1 - cos(theta))*u2^2 + cos(theta));

            A(5,6) = (u2*sin(theta) - u1*u3*(cos(theta) - 1))*(u3*sin(theta) - u1*u2*(cos(theta) - 1)) - (u1*sin(theta) + u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));
            A(6,5) = (u2*sin(theta) + u1*u3*(cos(theta) - 1))*(u3*sin(theta) + u1*u2*(cos(theta) - 1)) + (u1*sin(theta) - u2*u3*(cos(theta) - 1))*((1 - cos(theta))*u1^2 + cos(theta));

            obj.rotationMatrix = A;

        end
    end

end


text:
function [ a, isplit ] = r8vec_split ( n, a, split )

%*****************************************************************************80
%
%% R8VEC_SPLIT "splits" an unsorted R8VEC based on a splitting value.
%
%  Discussion:
%
%    If the vector is already sorted, it is simpler to do a binary search
%    on the data than to call this routine.
%
%    The vector is not assumed to be sorted before input, and is not
%    sorted during processing.  If sorting is not needed, then it is
%    more efficient to use this routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Input, real A(N), the array to split.
%
%    Input, real SPLIT, the value used to split the vector.
%    It is not necessary that any value of A actually equal SPLIT.
%
%    Output, real A(N), all the entries of A that are less than or
%    equal to SPLIT are in A(1:ISPLIT).
%
%    Output, integer ISPLIT, indicates the position of the last
%    entry of the split vector that is less than or equal to SPLIT.
%

%
%  Partition the vector into A1, A2, A3, where
%    A1 = A(I1:J1) holds values <= SPLIT,
%    A2 = A(I2:J2) holds untested values,
%    A3 = A(I3:J3) holds values > SPLIT.
%
  i1 = 1;
  j1 = 0;

  i2 = 1;
  j2 = n;

  i3 = n + 1;
  j3 = n;
%
%  Pick the next item from A2, and move it into A1 or A3.
%  Adjust indices appropriately.
%
  for i = 1 : n

    if ( a(i2) <= split )

      i2 = i2 + 1;
      j1 = j1 + 1;

    else

      temp = a(i2);
      a(i2) = a(i3-1);
      a(i3-1) = temp;

      i3 = i3 - 1;
      j2 = j2 - 1;

    end

  end

  isplit = j1;

  return
end


text:
function sparse_grid_mixed_weight_tests ( )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_WEIGHT_TESTS calls SPARSE_GRID_MIXED_WEIGHT_TEST with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    Local, real TOL, a tolerance for point equality.
%    A value of sqrt ( eps ) is reasonable, and will allow the code to
%    consolidate points which are equal, or very nearly so.  A value of
%    -1.0, on the other hand, will force the code to use every point, regardless
%    of duplication.
%
  addpath ( '../sandia_rules' );

  tol = sqrt ( eps );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TESTS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Call SPARSE_GRID_MIXED_WEIGHT with various arguments.\n' );
  fprintf ( 1, '  All tests will use a point equality tolerance of %e\n', tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 4 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 7 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 8 ]';
  alpha = [ 0.0, 1.5 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 2, 9 ]';
  alpha = [ 0.0, 0.5 ]';
  beta = [ 0.0, 1.5 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 6, 4 ]';
  alpha = [ 2.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 3;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 2, 5 ]';
  alpha = [ 0.0, 0.0, 0.0 ]';
  beta = [ 0.0, 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 3.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 3, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 13.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 13, 13 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 16.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 16, 16 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Level 4, Rule 17.
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 17, 17 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_weight_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_WEIGHT_TESTS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
function [Hx,Hxtest] = computeHx()
    Hx = eye(3,6);
    Hxtest = eye(6);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function [net, data] = gtminit2(net, options, data, init, samp_type, varargin)
%GTMINIT Initialise the weights and latent sample in a GTM given reference vectors.
%
%	Description
%	NET = GTMINIT(NET, OPTIONS, DATA, INIT, SAMPTYPE) takes a GTM NET and
%	generates a sample of latent data points and sets the centres (and
%	widths if appropriate) of NET.RBFNET.
%
% INIT defines initial locations of reference vectors (Gaussian means).
%
%	If the SAMPTYPE is 'REGULAR', then regular grids of latent data
%	points and RBF centres are created.  The dimension of the latent data
%	space must be 1 or 2.  For one-dimensional latent space, the
%	LSAMPSIZE parameter gives the number of latent points and the
%	RBFSAMPSIZE parameter gives the number of RBF centres.  For a two-
%	dimensional latent space, these parameters must be vectors of length
%	2 with the number of points in each of the x and y directions to
%	create a rectangular grid.  The widths of the RBF basis functions are
%	set by a call to RBFSETFW passing OPTIONS(7) as the scaling
%	parameter.
%
%	If the SAMPTYPE is 'UNIFORM' or 'GAUSSIAN' then the latent data is
%	found by sampling from a uniform or Gaussian distribution
%	correspondingly.  The RBF basis function parameters are set by a call
%	to RBFSETBF with the DATA parameter as dataset and the OPTIONS
%	vector.
%
%	Finally, the output layer weights of the RBF are initialised by
%	mapping the mean of the latent variable to the mean of the target
%	variable, and the L-dimensional latent variale variance to the
%	variance of the targets along the first L principal components.
%
%	See also
%	GTMINIT, GTM, GTMEM, PCA, RBFSETBF, RBFSETFW
%

%	Copyright (c) Ian T Nabney (1996-2001)
% Modified by Tommi Vatanen (2012)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

% Check type of sample
stypes = {'regular', 'uniform', 'gaussian', 'special'};
if (strcmp(samp_type, stypes)) == 0
  error('Undefined sample type.')
end

if net.dim_latent > size(data, 2)
  error('Latent space dimension must not be greater than data dimension')
end
nlatent = net.gmmnet.ncentres;
nhidden = net.rbfnet.nhidden;

% Create latent data sample and set RBF centres

switch samp_type
case 'regular'
   if nargin ~= 7
      error('Regular type must specify latent and RBF shapes');
   end
   l_samp_size = varargin{1};
   rbf_samp_size = varargin{2};
   if round(l_samp_size) ~= l_samp_size
      error('Latent sample specification must contain integers')
   end
   % Check existence and size of rbf specification
   if any(size(rbf_samp_size) ~= [1 net.dim_latent]) || ...
         prod(rbf_samp_size) ~= nhidden
      error('Incorrect specification of RBF centres')
   end
   % Check dimension and type of latent data specification
   if any(size(l_samp_size) ~= [1 net.dim_latent]) || ...
         prod(l_samp_size) ~= nlatent
      error('Incorrect dimension of latent sample spec.')
   end
   if net.dim_latent == 1
      net.X = [-1:2/(l_samp_size-1):1]';
      net.rbfnet.c = [-1:2/(rbf_samp_size-1):1]';
      net.rbfnet = rbfsetfw(net.rbfnet, options(7));
   elseif net.dim_latent == 2
      net.X = gtm_rctg(l_samp_size);
      net.rbfnet.c = gtm_rctg(rbf_samp_size);
      net.rbfnet = rbfsetfw(net.rbfnet, options(7));
   else
      error('For regular sample, input dimension must be 1 or 2.')
   end

case {'uniform', 'gaussian'}
   if strcmp(samp_type, 'uniform')
      net.X = 2 * (rand(nlatent, net.dim_latent) - 0.5);
   else
      % Sample from N(0, 0.25) distribution to ensure most latent
      % data is inside square
      net.X = randn(nlatent, net.dim_latent)/2;
   end
   net.rbfnet = rbfsetbf(net.rbfnet, options, net.X);

%   case {'special'}
%     net.X = [-0.5 sqrt(3)/2;0.5 sqrt(3)/2;0 0];
%     net.rbfnet.c = [-0.5 sqrt(3)/2;0.5 sqrt(3)/2;0 0];
%     net.rbfnet = rbfsetfw(net.rbfnet, options(7));
otherwise
   % Shouldn't get here
   error('Invalid sample type');
end

% use given unit vectors as the initial reference vectors
[~, Phi] = rbffwd(net.rbfnet, net.X);
% initparam are the reference vectors
net.rbfnet.w2 = Phi \ init;
net.rbfnet.b2 = zeros(1,size(data,2));

net.gmmnet.centres = rbffwd(net.rbfnet, net.X);
d = dist2(net.gmmnet.centres, net.gmmnet.centres) + ...
  diag(ones(net.gmmnet.ncentres, 1)*realmax);
sigma = mean(min(d))/2;

net.gmmnet.covars = sigma*ones(1, net.gmmnet.ncentres);

% Sub-function to create the sample data in 2d
function sample = gtm_rctg(samp_size)

xDim = samp_size(1);
yDim = samp_size(2);
% Produce a grid with the right number of rows and columns
[X, Y] = meshgrid([0:1:(xDim-1)], [(yDim-1):-1:0]);

% Change grid representation
sample = [X(:), Y(:)];

% Shift grid to correct position and scale it
maxXY= max(sample);
sample(:,1) = 2*(sample(:,1) - maxXY(1)/2)./maxXY(1);
sample(:,2) = 2*(sample(:,2) - maxXY(2)/2)./maxXY(2);

return;

author: ilarinieminen

repo: SOM-Toolbox


text:
classdef EnsembleClassifier < Algorithm

    properties (Access = public)
        shouldTrain = false;
        nLearners = 30;
    end

    properties (Access = private)
        classifier;
    end

    methods (Access = public)

        function obj = EnsembleClassifier()
            obj.name = 'Ensemble';
            obj.inputPort = DataType.kTable;
            obj.outputPort = DataType.kTable;
        end

        function dataOut = compute(obj,data)
            if obj.shouldTrain
                obj.train(data);
                dataOut = [];
            else
                dataOut = obj.test(data);
            end
        end

        function train(obj,table)
            template = templateTree(...
                'MaxNumSplits', 1852);

            if ~isempty(template)
                obj.classifier = fitcensemble(...
                    table.features, ...
                    table.label, ...
                    'Method', 'Bag', ...
                    'NumLearningCycles', obj.nLearners, ...
                    'Learners', template);
            end
        end

        function metrics = computeMetrics(obj,table)
            flops = timeit(@()obj.test(table)) / Constants.kReferenceComputingTime;
            memory = Helper.ComputeObjectSize(obj.classifier);
            outputSize = table.height * Constants.kClassificationResultBytes;
            metrics = Metric(flops,memory,outputSize);
        end

        function labels = test(obj,table)
            labels = predict(obj.classifier,table.features);
        end

        function str = toString(obj)
            str = sprintf('%s_%d_%d',obj.name,obj.nLearners);
        end

        function editableProperties = getEditableProperties(obj)
            editableProperties = Property('nLearners',obj.nLearners,1,100,PropertyType.kNumber);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function q=dcm2q(R)
% DCM2Q(R) converts direction cosine matrices into quaternions.
%
%     The resultant quaternion(s) will perform the equivalent vector
%     transformation as the input DCM(s), i.e.:
%
%       qconj(Q)*V*Q = R*V
%
%     where R is the DCM, V is a vector, and Q is the quaternion.  Note that
%     for purposes of quaternion-vector multiplication, a vector is treated
%     as a quaterion with a scalar element of zero.
%
%     If the input is a 3x3xN array, the output will be a vector of
%     quaternions where input direction cosine matrix R(:,:,k) corresponds
%     to the output quaternion Q(k,:).
%
% See also Q2DCM.

% Release: $Name: quaternions-1_2_2 $
% $Revision: 1.9 $
% $Date: 2002/01/21 06:46:20 $

% Copyright (C) 2000-2002, Jay A. St. Pierre.  All rights reserved.

if nargin~=1
  error('dcm2q() requires one input argument');
else
  size_R=size(R);
  if ( size_R(1)~=3 | size_R(2)~=3 | length(size_R)>3 )
    error(['Invalid input: must be a 3x3xN array'])
  end
end

q(1,:)=0.5*sqrt(1 + R(1,1,:) - R(2,2,:) - R(3,3,:)).*sgn(R(2,3,:)-R(3,2,:));
q(2,:)=0.5*sqrt(1 - R(1,1,:) + R(2,2,:) - R(3,3,:)).*sgn(R(3,1,:)-R(1,3,:));
q(3,:)=0.5*sqrt(1 - R(1,1,:) - R(2,2,:) + R(3,3,:)).*sgn(R(1,2,:)-R(2,1,:));
q(4,:)=0.5*sqrt(1 + R(1,1,:) + R(2,2,:) + R(3,3,:));

% Make quaternion vector a column of quaternions
q=q.';

q=real(q);

% Signum function (for our purposes)
function s=sgn(x)
  s=sign(x);
  s=s+(s==0);

author: christianwengert


text:
function [A_hat, E_hat, Y, mu, iter] = inexact_alm_rpca(D, lambda, tol, maxIter)

% Oct 2009
% This matlab code implements the inexact augmented Lagrange multiplier
% method for Robust PCA.
%
% D - m x n matrix of observations/data (required input)
%
% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%     - DEFAULT 1e-7 if omitted or -1.
%
% maxIter - maximum number of iterations
%         - DEFAULT 1000, if omitted or -1.
%
% Initialize A,E,Y,u
% while ~converged
%   minimize (inexactly, update A and E only once)
%     L(A,E,Y,u) = |A|_* + lambda * |E|_1 + <Y,D-A-E> + mu/2 * |D-A-E|_F^2;
%   Y = Y + \mu * (D - A - E);
%   \mu = \rho * \mu;
% end
%
% Minming Chen, October 2009. Questions? v-minmch@microsoft.com ;
% Arvind Ganesh (abalasu2@illinois.edu)
%
% Copyright: Perception and Decision Laboratory, University of Illinois, Urbana-Champaign
%            Microsoft Research Asia, Beijing

[m,n] = size(D);

if(nargin < 2) lambda = 1 / sqrt(m); end
if(nargin < 3) tol = 1e-7; elseif(tol == -1) tol = 1e-7; end
if(nargin < 4) maxIter = 1000; elseif(maxIter == -1) maxIter = 1000; end

% initialize
Y = D;
norm_two = norm(Y, 2);
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

A_hat = zeros( m, n);
E_hat = zeros( m, n);
mu = 1.25/norm_two; % this one can be tuned
mu_bar = mu * 1e7;
rho = 1.5;          % this one can be tuned
d_norm = norm(D, 'fro');

iter = 0;
total_svd = 0;
converged = false;
stopCriterion = 1;
sv = 10;
while ~converged
    iter = iter + 1;

    temp_T = D - A_hat + (1/mu)*Y;
    E_hat = max(temp_T - lambda/mu, 0);
    E_hat = E_hat+min(temp_T + lambda/mu, 0);

    %[U,S,V] = svd(D - E_hat + (1/mu)*Y, 'econ');
    [U,S,V] = svdecon(D - E_hat + (1/mu)*Y); % fastest

    diagS = diag(S);
    svp = length(find(diagS > 1/mu));
    if svp < sv
        sv = min(svp + 1, n);
    else
        sv = min(svp + round(0.05*n), n);
    end

    A_hat = U(:, 1:svp) * diag(diagS(1:svp) - 1/mu) * V(:, 1:svp)';

    total_svd = total_svd + 1;

    Z = D - A_hat - E_hat;

    Y = Y + mu*Z;
    mu = min(mu*rho, mu_bar);

    %% stop Criterion
    stopCriterion = norm(Z, 'fro') / d_norm;
    if stopCriterion < tol
        converged = true;
    end

    if mod( total_svd, 10) == 0
        disp(['#svd ' num2str(total_svd) ' r(A) ' num2str(rank(A_hat))...
            ' |E|_0 ' num2str(length(find(abs(E_hat)>0)))...
            ' stopCriterion ' num2str(stopCriterion)]);
    end

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end
end

author: andrewssobral

repo: lrslibrary


text:
function state = pbsGetJobState(scheduler, job, state)
%pbsGetJobState Gets the state of a job on a cluster.
%
% Set your schedulers's GetJobStateFcn to this function using the following
% command (see README):
%     set(sched, 'GetJobStateFcn', @pbsGetJobState);

%  Copyright 2006-2008 The MathWorks, Inc.

mlock;
persistent jobsToMonitorNames;

if isempty(jobsToMonitorNames)
    jobsToMonitorNames = {};
end

if isempty(scheduler.UserData)
    if ~iscell(scheduler.SubmitFcn) || length(scheduler.SubmitFcn) < 3
        error('distcomp:genericscheduler:SubmitFcnError',...
              'SubmitFcn must include clusterHost and remoteDataLocation as extra arguments.');
    end
    scheduler.UserData={ scheduler.SubmitFcn{2} ; scheduler.SubmitFcn{3} };
end

jobName = job.pGetEntityLocation;
if strcmp(state, 'finished')
% emc 6 March 2010 replace setxor with setdiff - behavior with empty set
      jobsToMonitorNames = setdiff(jobsToMonitorNames, {jobName});
else
%   timerName as suggested by Ray Norris TMW 5 March 2010
%   timerName = [job.Name '_' num2str(job.ID)];
%   timerName emc 6 March 2010
    timerName = [job.Name '_' strrep(job.SubmitTime, ' ','_')];
    if strcmp(state, 'queued') || strcmp(state, 'running')
        if ~ismember(jobName, jobsToMonitorNames) || ...
                                    isempty(timerfind('Name',timerName))
            jobsToMonitorNames = union(jobsToMonitorNames, {jobName});
            jobStateTimer = ...
                      timer('Period', 10.0, 'ExecutionMode', 'fixedRate');
            set(jobStateTimer, 'TimerFcn', ...
                         { @copyJobFilesIfFinished, scheduler, job }, ...
                                                        'Name',timerName);
            start(jobStateTimer);
        end
    end
end


text:
% Load data "bias60.mat".
clc;
clear;
load bias60;

% Create a figure.
figure('DefaultAxesFontSize',16);
ax = gca;

% Create a heatmap.
mat = bias{2}.*ones(61,144)+bias{3}'.*ones(61,144);
imagesc(mat);
c = colorbar;
colormap hot;
ylabel(c,'Speed (km/h)');
ax.XTick = [1,13,25,37,49,61,73,85,97,109,121,133,144];
ax.XTickLabel = ['00:00';'02:00';'04:00';'06:00';'08:00';'10:00';'12:00';...
    '14:00';'16:00';'18:00';'20:00';'22:00';'23:50'];
xtickangle(45);
ax.YTick = [1,6,11,16,21,26,31,36,41,46,51,56,61];
ax.YTickLabel = ['08-01';'08-06';'08-11';'08-16';'08-21';'08-26';'08-31';...
    '09-05';'09-10';'09-15';'09-20';'09-25';'09-30'];

% Save the figure as "heat_map60.pdf".
set(gcf, 'PaperSize', [6 4.5]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 6 4.5]);
saveas(gcf,'heat_map60','pdf');

author: xinychen

repo: academic-drawing


text:
init_fig;

% c
axes('position', [.05, .57, .95, .37]);
hold on;
plot(y, 'color', col{8}/255);
alpha(.7);
plot(true_c, 'color', col{3}/255, 'linewidth', 1.5);
plot(c_oasis, '-.', 'color', col{1}/255);
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', []);
set(gca, 'ytick', 0:2);
ylabel('Fluor.');
box off;
legend('Data', 'Truth', 'OASIS', 'location', 'northeast', 'orientation', 'horizental');
% s
axes('position', [.05, .18, .95, .37]);
hold on;
plot(true_s, 'color', col{3}/255, 'linewidth', 1.5);
plot(s_oasis, '-.', 'color', col{1}/255);
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', get(gca, 'xtick')/30);
set(gca, 'ytick', [0,1]);
xlabel('Time [s]');
ylabel('Activity.');

author: zhoupc


text:
% This file load eq catalogs in hypoellipes format. Since the formated reading
% is fast, this is the best way to load data for
% large catalogs.
% Earthquakes within a rectengular box (latmin latmax lonmin lonmax)
% a time window and a magnitude threshold can be selected.
% Note: lat/Lon are assumed to be in minutes in the original catalog
%
% last modified: March 96   stefan wiemer

report_this_filefun(mfilename('fullpath'));

%  close the file in case it is still open
safe_fclose(fid);

% selection criteria
tmin = 60.0001;
tmax = 97.00900;
lonmin = -180.0;
lonmax =  180.0;
latmin = -90.0;
latmax = 90.0 ;
Mmin = 3.1
deperr = 1000.5;

a = []; b = [];
n2 = 0;

% open the file and read 2000 lines at a time
[file1,path1] = uigetfile([ '*.mat'],' Earthquake Datafile');
fid = fopen([path1 file1],'r') ;;

while  ferror(fid) == ''
    l =fscanf(fid,'%2d%2d%2d%2d%2d%4d%2d%*1c%4d%3d%*1c%4d%5f%2d%*31c%2d%*15c%4d',[14 2000]) ;
    n2 = n2 +  length(l(1,:));

    b = [ -l(9,:)-l(10,:)/6000 ; l(7,:)+l(8,:)/6000 ;...
        l(1,:);l(2,:);l(3,:);l(13,:)/10;l(11,:)/100;l(4,:);...
        l(5,:); l(14,:)/100;l(12,:)/10];
    b = b';
    l =  b(:,6) >= Mmin & b(:,1) >= lonmin & b(:,1) <= lonmax & ...
        b(:,2) >= latmin & b(:,2) <= latmax & b(:,3) <= tmax  & ...
        b(:,3) >= tmin & b(:,10) <= deperr;
    a = [a ; b(l,:)];
    disp([ num2str(n2) ' earthquakes scanned, ' num2str(length(a)) ' EQ found'])
    if max(b(:,3)) >  tmax ; break; end
end
ferror(fid)
fclose(fid);

dep1 = 0.3*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep2 = 0.6*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep3 = max(a.Depth);

stri1 = [file1];
tim1 = minti;
tim2 = maxti;
minma2 = minma;
maxma2 = maxma;
minde = min(a.Depth);
maxde = max(a.Depth);
rad = 50.;
ic = 0;
ya0 = 0.;
xa0 = 0.;
iwl3 = 1.;
step = 3;
t1p(1) = 80.;
t2p(1) = 85.;
t3p(1) = 90.;
t4p(1) = 93.;
tresh = 10;

% save the catalog
save newcat.mat a

author: CelsoReyes


text:
% VL_DEMO_QUICKSHIFT  Demo: Quick shift: basic functionality

pfx = fullfile(vl_root,'figures','demo') ;
randn('state',0) ;
rand('state',0) ;
figure(1) ; clf ;

% --------------------------------------------------------------------
%                                                        Load a figure
% --------------------------------------------------------------------
I = imread(fullfile(vl_root,'data','roofs1.jpg')) ;

image(I) ;
axis equal off tight ;
vl_demo_print('quickshift_image') ;

% --------------------------------------------------------------------
%                                     Create a quickshift segmentation
% --------------------------------------------------------------------

ratio = 0.5;
kernelsize = 2;
maxdist = 10;

Iseg = vl_quickseg(I, ratio, kernelsize, maxdist);
image(Iseg);
axis equal off tight;
vl_demo_print('quickshift_qseg_1') ;

kernelsize = 2;
maxdist = 20;
Iseg = vl_quickseg(I, ratio, kernelsize, maxdist);
image(Iseg);
axis equal off tight;
vl_demo_print('quickshift_qseg_2') ;

maxdist = 50;
ndists = 10;
Iedge = vl_quickvis(I, ratio, kernelsize, maxdist, ndists);
imagesc(Iedge);
axis equal off tight;
colormap gray;
vl_demo_print('quickshift_qvis') ;

author: Cloud-CV

repo: object-proposals