text:
% fig2d size histogram for all clusters within one fish

% load all clusters (Fish6_woA_Master0.7)
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);
i_fish = 6;
[cIX,gIX,M] = LoadSingleFishDefault(i_fish,hfig);

%% fig4c: Distribution of cluster sizes
U = unique(gIX);
numU = length(U);
C = zeros(numU,1);
for i=1:numU,
    IX = find(gIX == U(i));
    C(i) = length(IX);
end
% [N,edges] = histcounts(C,10:10:2100);
%
% figure('Position',[500,500,350,300]);hold on;
% for i = 1:length(N),
%     plot([edges(i),edges(i)],[0,N(i)],'color',[1,0.5,0.5],'linewidth',6)
% end
% set(gca,'XScale','log')

% try making log-scaled bins for the histogram
% bins = 10:10:2100; % enough for Autoclus0.7
bins = 10:100:10700;
logbins = log(bins);
linbins = linspace(logbins(1),logbins(end),10);
xbins = exp(linbins);

[N,edges,bin] = histcounts(C,xbins);
%
figure('Position',[500,500,200,200]);hold on; % [500,500,350,300]
for i = 1:length(N),
    plot([edges(i),edges(i)],[0,N(i)],'color',[1,0.5,0.5],'linewidth',8)
end
% bar(edges(1:length(N)),N)
set(gca,'XScale','log')
xlim([5,10^4])
set(gca,'XTick',[10,100,1000,10000])
% ylim([0,70])
xlabel('cluster size')
ylabel('count')

%% fig4d: Distribution of within-cluster correlations
U = unique(gIX);
numU = length(U);
B = zeros(numU,3);
for i=1:numU,
    IX = find(gIX == U(i));
    coeffs = corr(M(IX,:)');
    m = coeffs(:);
%     B(i,1) = min(m);
    B(i,2) = mean(m);
%     B(i,3) = median(m);
end

figure('Position',[500,500,250,150]);
hist(B(:,2),0:0.05:1)
h = findobj(gca,'Type','patch');
h.FaceColor = [0.5 0.5 0.5];
h.EdgeColor = 'w';
xlim([0,1])
set(gca,'XTick', 0:0.2:1);
xlabel('average corr. within cluster')
ylabel('count')

%% fig4a... kmeans
i_ClusGroup = 2;% 3;
i_Cluster = 8; % 1;
[cIX,gIX] = LoadCluster_Direct(i_fish,i_ClusGroup,i_Cluster);
M = UpdateIndices_Manual(hfig,cIX,gIX,numU);

%% left plot
figure('Position',[50,100,800,1000]);
% isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
setappdata(hfig,'isPlotBehavior',1);
setappdata(hfig,'isStimAvr',0);
UpdateTimeIndex(hfig);
DrawTimeSeries(hfig,cIX,gIX);

%%

author: xiuyechen

repo: FishExplorer


text:
function vscl = vscale(f)
%VSCALE   Vertical scale of a CHEBFUN3.
%   VSCL = VSCALE(F) returns the vertical scale of a CHEBFUN3 object F as
%   determined by evaluating F on a coarse tensor-product grid.

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

% [TODO]: Should this also be taking the maximum along the edges when we are
% evaluating at 1st kind grids.

% If f is an empty CHEBFUN3, VSCL = 0:
if ( isempty(f) )
    vscl = 0;
    return
end

techCol = get(f.cols.funs{1}, 'tech');

% Get the degree of the CHEBFUN3:
[m, n, p] = length(f);

% If F is of low degree, then oversample:
m = min(max(m, 9), 41);
n = min(max(n, 9), 41);
p = min(max(p, 9), 41); % cannot afford to go over 41x41x41.

% Calculate values on a tensor grid:
dom = f.domain;
x = mypoints(m, dom(1:2), techCol);
y = mypoints(n, dom(3:4), techCol);
z = mypoints(p, dom(5:6), techCol);
[xx, yy, zz] = ndgrid(x, y, z);
vals = feval(f, xx, yy, zz);

% Take the absolute maximum:
vscl = max(abs(vals(:)));

end

%%
function x = mypoints(n, dom, tech)
% Get the sample points that correspond to the right grid for a particular
% technology.

% What tech am I based on?:
if ( isa(tech(), 'chebtech2') )
    x = chebpts(n, dom, 2);
elseif ( isa(tech(), 'chebtech1') )
    x = chebpts(n, dom, 1);
elseif ( isa(tech(), 'trigtech') )
    x = trigpts(n, dom);
else
    error('CHEBFUN:CHEBFUN3:vscale:mypoints:techType', ...
        'Unrecognized technology');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef TestEstimateRigidTransform
    %TestEstimateRigidTransform

    methods (Static)
        function test_1
            im1 = 255*uint8([...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 1 1 1 1 0 0 0;...
                0 0 0 1 0 1 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 0 0 0 1 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
            ]);
            im2 = circshift(im1, [0 1]);
            M = cv.estimateRigidTransform(im1, im2, 'FullAffine',false);
        end

        function test_images_input
            im1 = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale1.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            im2 = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale2.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            M = cv.estimateRigidTransform(im1, im2, 'FullAffine',true);
            validateattributes(M, {'double'}, {'size',[2 3]});
        end

        function test_points_input
            N = 50;
            pts1 = rand(2,N)*10;                  % a set of 2D points
            aff = randn(2,3);                     % true affine transformation
            pts2 = aff * [pts1; ones(1,N)];       % transform pts1
            pts2 = pts2 + randn(size(pts2))*0.01; % add noise
            M = cv.estimateRigidTransform(...
                num2cell(pts1,1), num2cell(pts2,1), 'FullAffine',true);
            validateattributes(M, {'double'}, {'size',[2 3]});
            err = norm(M - aff);
        end

        function test_error_argnum
            try
                cv.estimateRigidTransform();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [beadsAt,bestslice]=ExtractCoordinates(myres,s)
%Error: Undefined function or method 'ExtractCoordiantes' for input arguments of type 'dip_image'

[mymax,maxpos]=max(myres,[],3);
mymax
fprintf('Please select bead coordinates by left clicking. Finish with right click\n');
beadsAt=dipgetcoords(100);
fprintf('%d beads selected\n',size(beadsAt,1)-1);
beadsAt(end,:)=[];

for b=1:size(beadsAt,1)
    mypos=beadsAt(b,:);
    mybead=extract(squeeze(mymax),[s s],mypos);

    mymask=mybead >= max(mybead);
    mybeadslice=extract(squeeze(maxpos),[s s],mypos);
    bestslice(b)= round(mean(mybeadslice,mymask));

%    Version 1
%     [mm,mpx]=max(mybead,[],1); mpx=squeeze(mpx);
%     [mm,mppx]=max(squeeze(mm),[],1);
%     bestslice(b,1) = double(mpx(mppx));
%     [mm,mpy]=max(mybead,[],2);mpy=squeeze(mpy);
%     [mm,mppy]=max(squeeze(mm),[],1);
%     bestslice(b,2) = double(mpy(mppy));

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% script that demonstrates use of 2d scattering

clear; close all;
x = uiuc_sample;
filt_opt.J = 7;
filt_opt.L = 8;
scat_opt.oversampling = 0;
[Wop, filters] = wavelet_factory_2d(size(x), filt_opt, scat_opt);

%%
profile on;
tic;
[Sx, Ux] = scat(x, Wop);
toc;
profile off;
profile viewer;

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [pvec, pstruct] = tapas_softmax_binary_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1)    = exp(ptrans(1));         % be
pstruct.be = pvec(1);

return;

author: translationalneuromodeling

repo: tapas


text:
% Copyright (C) 2007 Peter Carbonetto. All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% Author: Peter Carbonetto
%         Dept. of Computer Science
%         University of British Columbia
%         May 19, 2007

function qR = bopt (K, C, f, Rv, Rf, Sv, Sf, NS, verbose)
  maxiter   = 100;         % The maximum number of iterations.
  tolerance = 1e-6;        % The convergence criterion.
  nr        = length(Rv);  % The number of large regions.
  ns        = length(Sv);  % The number of small regions.

  % Convert the input K to a row vector.
  K = K(:)';

  % Get the degrees of the separators.
  d = degrees(NS);

  % Initialize the marginals on the large and small regions to uniform
  % probability tables.
  qR = initmarginals(K,Rv);
  qS = initmarginals(K,Sv);

  % Compute the total number consistency constaints. We have one
  % constraint for every pair (R,S), then again for every possible
  % configuration xS.
  nc = numconsistencyconstraints(qS,d);

  % Run the IPOPT solver.
  qR             = vectorize(qR);
  qS             = vectorize(qS);
  nqr            = length(qR);
  nqs            = length(qS);
  [status qR qS] = ipopt({qR qS},{ repmat(eps,nqr,1) repmat(eps,nqs,1) },...
                  { repmat(inf,nqr,1) repmat(inf,nqs,1) },...
                  [ ones(1,nr) ones(1,ns) zeros(1,nc) ],...
                  [ ones(1,nr) ones(1,ns) zeros(1,nc) ],...
                  @computeJGObjective,@computeJGGradient,...
                  @computeJGConstraints,@computeJGJacobian,...
                  @computeJGHessian,{ K C f Rv Rf Sv Sf NS d },'',...
                  [],'mu_strategy','adaptive','max_iter',maxiter,...
                  'tol',tolerance,'jac_c_constant','yes',...
		  'jac_d_constant','yes','print_level',verbose*5);

  % Reshape the solution.
  qR = reshapemarginals(qR,Rv,K);
  qS = reshapemarginals(qS,Sv,K);

% ------------------------------------------------------------------
function d = degrees (NS)
  ns = length(NS);  % The number of separators.
  d  = zeros(1,ns);

  % Repeat for each separator.
  for s = 1:ns
    d(s) = length(NS{s});
  end

% ----------------------------------------------------------------
function qR = initmarginals (K, Rv)
  nr = length(Rv);  % The number of regions.
  qR = cell(1,nr);  % The return value.

  % Repeat for each region.
  for r = 1:nr
    is    = Rv{r};  % The variables nodes in the large region.
    table = ones([K(is) 1]);
    qR{r} = table / sum(table(:));
  end

% ----------------------------------------------------------------
function nc = numconsistencyconstraints (qS, d)
  ns = length(qS);  % The number of small regions.
  nc = 0;           % The return value.

  % Repeat for each small region.
  for s = 1:ns
    tablesize = numel(qS{s});
    nc        = nc + d(s) * tablesize;
  end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [xMid,fMid,exitCode]=bisectionRootFind(f,xSpan,AbsTol,maxIter)
%%BISECTIONROOTFIND Given two points that are of different signs
%          surrounding the root of a scalar funding, find the root of the
%          function via a bisection search.
%
%INPUTS: f The handle to the function whose root is desired. The function
%          is called as f(x), where x is a scalar value.
%    xSpan A 2X1 or 1X2 vector such that
%          sign(f(xSpan(1)))~=sign(f(xSpan(2)))
%   AbsTol The absolute tolerance on the value of f for converence. The
%          default if omitted or an empty matrix is passed is eps().
%  maxIter The maximum number of bisections to perform. The default if
%          omitted or an empty matrix is passed is 100.
%
%OUTPUTS: xMid The point found that is the value of the approximate zero.
%         fMid The value of the function at xMid.
%     exitCode A value indicating how the function terminated. Possible
%              values are:
%              0 Either AbsTol was satisfied or the bounds on the value
%                xMid went below finite precision constraints.
%              1 The maximum number of iterations occurred.
%
%The bisection algorithm is straightforward: Assuming that only a single
%root is surrounded by xSpan, one evaluates the point at the middle of the
%interval. One then discards the point in xSpan having the same sign as the
%point in the middle of the interval. That gives a new interval and the
%process continues until convergence. The method is discussed in [1].
%
%EXAMPLE:
% f=@(x)(x.^2-4);
% xSpan=[0;5];
% [xMid,fMid,exitCode]=bisectionRootFind(f,xSpan)
%In this example convergence is to the exact root, xMid=2.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Bisection." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/Bisection.html
%
%July 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(maxIter))
    maxIter=100;
end

if(nargin<3||isempty(AbsTol))
    AbsTol=eps();
end

xLeft=xSpan(1);
xRight=xSpan(2);

fLeft=f(xLeft);
fRight=f(xRight);

if(sign(fLeft)==sign(fRight))
    error('bisectionRootFind:NoChange','There is no sign change in the initial set. ');
end

xMid=xRight;
fMid=fRight;

exitCode=1;
for curIter=1:maxIter
    xMid=(xLeft+xRight)/2;
    fMid=f(xMid);

    if(abs(fMid)<AbsTol||xMid==xRight||xMid==xLeft)
        exitCode=0;
        break;
    end

    if(sign(fLeft)~=sign(fMid))
        xRight=xMid;
    else
        fLeft=fMid;
        xLeft=xMid;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function digest = md5(message)
    % digest = md5(message)
    %  Compute the MD5 digest of the message, as a hexadecimal digest.
    %  based on the source code from https://rosettacode.org/wiki/MD5

    % Follow the MD5 algorithm from RFC 1321 [1] and Wikipedia [2].
    %  [1] http://tools.ietf.org/html/rfc1321
    %  [2] http://en.wikipedia.org/wiki/MD5

    % m is the modulus for 32-bit unsigned arithmetic.
    m = 2 ^ 32;

    % s is the shift table for circshift(). Each shift is negative
    % because it is a left shift.
    s = [-7, -12, -17, -22
         -5,  -9, -14, -20
         -4, -11, -16, -23
         -6, -10, -15, -21];

    % t is the sine table. Each sine is a 32-bit integer, unsigned.
    t = floor(abs(sin(1:64)) .* m);

    % Initialize the hash, as a row vector of 32-bit integers.
    digest = [hex2dec('67452301') ...
              hex2dec('EFCDAB89') ...
              hex2dec('98BADCFE') ...
              hex2dec('10325476')];

    % If message contains characters, convert them to ASCII values.
    message = double(message);
    bytelen = numel(message);

    % Pad the message by appending a 1, then appending enough 0s to make
    % the bit length congruent to 448 mod 512. Because we have bytes, we
    % append 128 '10000000', then append enough 0s '00000000's to make
    % the byte length congruent to 56 mod 64.
    message = [message, 128, zeros(1, mod(55 - bytelen, 64))];

    % Convert the message to 32-bit integers, little endian.
    % For little endian, first byte is least significant byte.
    message = reshape(message, 4, numel(message) / 4);
    message = message(1,:) + ...            % least significant byte
              message(2,:) * 256 + ...
              message(3,:) * 65536 + ...
              message(4,:) * 16777216;      % most significant byte

    % Append the bit length as a 64-bit integer, little endian.
    bitlen = bytelen * 8;
    message = [message, mod(bitlen, m), mod(bitlen / m, m)];

    % Process each 512-bit block. Because we have 32-bit integers, each
    % block has 16 elements, message(k + (0:15)).
    for k = 1:16:numel(message)
        % Copy hash.
        a = digest(1); b = digest(2); c = digest(3); d = digest(4);

        % Do 64 operations.
        for i = (1:64)
            % Convert b, c, d to row vectors of bits (0s and 1s).
            bv = dec2bin(b, 32) - '0';
            cv = dec2bin(c, 32) - '0';
            dv = dec2bin(d, 32) - '0';

            % Find f  = mix of b, c, d.
            %      ki = index in 0:15, to message(k + ki).
            %      sr = row in 1:4, to s(sr, :).
            if i <= 16          % Round 1
                f = (bv & cv) | (~bv & dv);
                ki = i - 1;
                sr = 1;
            elseif i <= 32      % Round 2
                f = (bv & dv) | (cv & ~dv);
                ki = mod(5 * i - 4, 16);
                sr = 2;
            elseif i <= 48      % Round 3
                f = xor(bv, xor(cv, dv));
                ki = mod(3 * i + 2, 16);
                sr = 3;
            else                % Round 4
                f = xor(cv, bv | ~dv);
                ki = mod(7 * i - 7, 16);
                sr = 4;
            end

            % Convert f, from row vector of bits, to 32-bit integer.
            f = bin2dec(char(f + '0'));

            % Do circular shift of sum.
            sc = mod(i - 1, 4) + 1;
            sum = mod(a + f + message(k + ki) + t(i), m);
            sum = dec2bin(sum, 32);
            sum = circshift(sum, [0, s(sr, sc)]);
            sum = bin2dec(sum);

            % Update a, b, c, d.
            temp = d;
            d = c;
            c = b;
            b = mod(b + sum, m);
            a = temp;
        end %for i

        % Add hash of this block to hash of previous blocks.
        digest = mod(digest + [a, b, c, d], m);
    end %for k

    % Convert hash from 32-bit integers, little endian, to bytes.
    digest = [digest                % least significant byte
              digest / 256
              digest / 65536
              digest / 16777216];   % most significant byte
    digest = reshape(mod(floor(digest), 256), 1, numel(digest));

    % Convert hash to hexadecimal.
    digest = dec2hex(digest);
    digest = lower(reshape(transpose(digest), 1, numel(digest)));
end %md5

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [Com] = Data_Input(Time,ComandoC,ComandoC0,startCC,ComandoR,ComandoR0,startCR,ComandoS,FreqS,ComandoS0,FreqS0,startCS)

if(Time < startCC)
    ComC=ComandoC0;
else
    ComC=ComandoC;
end
if(Time < startCR)
    ComR=ComandoR0;
    startR=0;
    K=0;
else
    ComR=ComandoR;
    startR=startCR;
    K=ComandoR0*startCR;
end
if(Time < startCS)
    ComS=ComandoS0;
    Freq=FreqS0;
    T=0;
else
    ComS=ComandoS;
    Freq=FreqS;
    T=startCS;
end
Com=ComC+K+ComR*(Time-startR)+ComS*sin(Freq*(Time-T));
end

author: Sable

repo: mcbench-benchmarks


text:
function epsilon = epsm(units)

    %EPSM  Calculate the accuracy of the map computations
    %
    %  e = EPSM returns the accuracy of computations performed in
    %  the Mapping Toolbox.  The accuracy returned is in degrees.
    %
    %  e = EPSM('units') returns the accuracy in the units specified
    %  by the string 'units'.  If omitted, 'degrees' are assumed.
    %
    %  See also EPS

    %  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
    %  Written by:  E. Byrns, E. Brown
    %   $Revision: 1399 $    $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    %  Define the limit in degrees

    %report_this_filefun(mfilename('fullpath'));

    degepsilon = 1.0E-6;

    if nargin == 0
        epsilon = degepsilon;   return

        %  Speed up function with special unit string tests

    elseif strcmp(units,'degrees')
        epsilon = degepsilon;   return

    elseif strcmp(units,'radians')
        epsilon = degepsilon*pi/180;   return

    else
        [units,msg] = unitstr(units,'angles');
        if ~isempty(msg);   error(msg);   end
        epsilon = angledim(degepsilon,'degrees',units);
    end
end

author: CelsoReyes


text:
%% Export Pole Figure Data
%
%%
%

author: mtex-toolbox

repo: mtex

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


text:
function a = r8cc_inc ( m, n, nz_num, colptr, rowind, a, i, j, aij )

%*****************************************************************************80
%
%% R8CC_INC increments a value of a R8CC matrix.
%
%  Discussion:
%
%    The R8CC format is the double precision sparse compressed column
%    format.  Associated with this format, we have an M by N matrix
%    with NZ_NUM nonzero entries.  We construct the column pointer
%    vector COL of length N+1, such that entries of column J will be
%    stored in positions COL(J) through COL(J+1)-1.  This indexing
%    refers to both the ROW and A vectors, which store the row indices
%    and the values of the nonzero entries.  The entries of the
%    ROW vector corresponding to each column are assumed to be
%    ascending sorted.
%
%    The R8CC format is equivalent to the MATLAB "sparse" format,
%    and the Harwell Boeing "real unsymmetric assembled" (RUA) format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%
%    Input, integer N, the number of columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero entries.
%
%    Input, integer COLPTR(N+1), indicate where each column's data begins.
%
%    Input, integer ROWIND(NZ_NUM), the row indices.
%
%    Input, real A(NZ_NUM), the nonzero entries.
%
%    Input, integer I, J, the indices of the value to retrieve.
%
%    Input, real AIJ, the value to be added to A(I,J).
%
%    Output, real A(NZ_NUM), entry (I,J) has been incremented.
%

%
%  Seek sparse index K corresponding to full index (I,J).
%
  k = r8cc_ijk ( m, n, nz_num, colptr, rowind, i, j );
%
%  If no K was found, we fail.
%
  if ( k == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8CC_INC - Fatal error!\n' );
    fprintf ( 1, '  R8CC_IJK could not find the entry.\n' );
    fprintf ( 1, '  Row I = %d\n', i );
    fprintf ( 1, '  Col J = %d\n', j );
    error ( 'R8CC_INC - Fatal error!' );
  end

  a(k) = a(k) + aij;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function disp_lines(rectimg, lines, label, showimage)

if ~exist('showimage','var') || showimage==1
	figure; imshow(rectimg); hold on;
else
	hold on;
end

if nargin<3
    for k = 1:length(lines)
        xy = [lines(k).point1; lines(k).point2];
%         col = [0 0 1];
        col = rand(1,3);
        plot(xy(:,1),xy(:,2), 'LineWidth',2, 'Color', col);
        plot(xy(1,1), xy(1,2), 'x', 'LineWidth',2, 'Color', col);
        text(mean(xy(:,1)),mean(xy(:,2)), num2str(k), 'Color', col);
    end
else
    for k = 1:length(lines)
        xy = [lines(k).point1; lines(k).point2];
text(xy(1,1), xy(1,2), sprintf('%d',k));
        if label(k) == 1
            plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','red');
        elseif label(k) == -1
            plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','blue');
        elseif label(k) == 0
            plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
        end
    end
end

author: zouchuhang

repo: LayoutNet


text:
%RNA
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 4/22/2011
classdef RNA
    methods (Static = true)
        function run(sim, nSample, nIter, fileName)
            import edu.stanford.covert.cell.sim.analysis.RNA;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.PrintUtil;

            %simulate
            [expRnaSyn, expRnaExp] = RNA.calcExpectedSynthesisExpression(sim);
            [simRnaSyn, simRnaExp] = RNA.sampleRNASynthesis(sim, nSample, nIter);

            %print
            [content, colLabels, indentation] = RNA.tabulateRNASynthesisExpression(sim, expRnaSyn, expRnaExp, simRnaSyn, simRnaExp);
            if nargin == 3
                PrintUtil.printToStdIO(content, colLabels, struct('indentation', indentation));
            else
                PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'Synthesis-Expression', struct('indentation', indentation));
            end

            %plot
            if nargin == 3
                RNA.plotRnaSynthesis(PlotUtil.newAxesHandle(), sim, expRnaSyn, simRnaSyn);
                RNA.plotRnaExpression(PlotUtil.newAxesHandle(), sim, expRnaExp, simRnaExp);
            else
                [axesHandle, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();

                RNA.plotRnaSynthesis(axesHandle, sim, expRnaSyn, simRnaSyn);
                saveas(figHandle, [fileName '-Synthesis.pdf']);

                RNA.plotRnaExpression(axesHandle, sim, expRnaExp, simRnaExp);
                saveas(figHandle, [fileName '-Expression.pdf']);

                close(figHandle);
            end
        end
    end

    methods (Static = true)
        function [content, colLabels, indentation] = tabulateRNASynthesisExpression(sim, expRnaSyn, expRnaExp, simRnaSyn, simRnaExp)
            rna = sim.state('Rna');
            matIdxs = rna.matureIndexs;

            colLabels = [
                'ID', 'Name', 'Transcription Unit', ...
                'Expected Synthesis (molecules / cell cycle)', 'Mean Simulated Synthesis (molecules / cell cycle)', 'Std Simulated Synthesis (molecules / cell cycle)', ...
                cellfun(@(x) ['Synthesis - Simulation #' num2str(x)], num2cell(1:size(simRnaSyn, 2)), 'UniformOutput', false), ...
                'Expected Expression (molecules * s / cell cycle)', 'Mean Simulated Expression (molecules * s / cell cycle)', 'Std Simulated Expression (molecules * s / cell cycle)', ...
                cellfun(@(x) ['Expression - Simulation #' num2str(x)], num2cell(1:size(simRnaExp, 2)), 'UniformOutput', false), ...
                ];

            [~, tuIdxs] = find(rna.nascentRNAMatureRNAComposition);

            content = [
                rna.wholeCellModelIDs(matIdxs) rna.names(matIdxs) rna.wholeCellModelIDs(tuIdxs) ...
                num2cell([
                expRnaSyn mean(simRnaSyn, 2) std(simRnaSyn, 0, 2) simRnaSyn ...
                expRnaExp mean(simRnaExp, 2) std(simRnaExp, 0, 2) simRnaExp ...
                ])];

            indentation = zeros(size(content, 1), 1);
        end

        function plotRnaSynthesis(axesHandle, sim, expRnaSyn, simRnaSyn)
            rna = sim.state('Rna');

            hold(axesHandle, 'on');

            h = [
                plot(axesHandle, expRnaSyn(rna.matureMRNAIndexs), mean(simRnaSyn(rna.matureMRNAIndexs, :), 2), 'r.')
                plot(axesHandle, expRnaSyn(rna.matureRRNAIndexs), mean(simRnaSyn(rna.matureRRNAIndexs, :), 2), 'g.')
                plot(axesHandle, expRnaSyn(rna.matureSRNAIndexs), mean(simRnaSyn(rna.matureSRNAIndexs, :), 2), 'b.')
                plot(axesHandle, expRnaSyn(rna.matureTRNAIndexs), mean(simRnaSyn(rna.matureTRNAIndexs, :), 2), 'c.')
                ];
            line([0 max(max(expRnaSyn), max(mean(simRnaSyn, 2)))], [0 max(max(expRnaSyn), max(mean(simRnaSyn, 2)))], 'Parent', axesHandle);
            xlim(axesHandle, [0 max(max(expRnaSyn), max(mean(simRnaSyn, 2)))]);
            ylim(axesHandle, [0 max(max(expRnaSyn), max(mean(simRnaSyn, 2)))]);
            axis(axesHandle, 'square');
            box(axesHandle, 'on');
            xlabel(axesHandle, 'Expected Synthesis', 'FontSize', 12);
            ylabel(axesHandle, 'Mean Simulated Synthesis', 'FontSize', 12);
            legend(h, {'mRNA', 'rRNA', 'sRNA', 'tRNA'}, 'Location', 'NorthEastOutside');
        end

        function plotRnaExpression(axesHandle, sim, expRnaExp, simRnaExp)
            rna = sim.state('Rna');

            hold(axesHandle, 'on');

            h = [
                plot(axesHandle, expRnaExp(rna.matureMRNAIndexs), mean(simRnaExp(rna.matureMRNAIndexs, :), 2), 'r.')
                plot(axesHandle, expRnaExp(rna.matureRRNAIndexs), mean(simRnaExp(rna.matureRRNAIndexs, :), 2), 'g.')
                plot(axesHandle, expRnaExp(rna.matureSRNAIndexs), mean(simRnaExp(rna.matureSRNAIndexs, :), 2), 'b.')
                plot(axesHandle, expRnaExp(rna.matureTRNAIndexs), mean(simRnaExp(rna.matureTRNAIndexs, :), 2), 'c.')
                ];
            line([0 max(max(expRnaExp), max(mean(simRnaExp, 2)))], [0 max(max(expRnaExp), max(mean(simRnaExp, 2)))], 'Parent', axesHandle);
            xlim(axesHandle, [0 max(max(expRnaExp), max(mean(simRnaExp, 2)))]);
            ylim(axesHandle, [0 max(max(expRnaExp), max(mean(simRnaExp, 2)))]);
            axis(axesHandle, 'square');
            box(axesHandle, 'on');
            xlabel(axesHandle, 'Expected Expression', 'FontSize', 12);
            ylabel(axesHandle, 'Mean Simulated Expression', 'FontSize', 12);
            legend(h, {'mRNA', 'rRNA', 'sRNA', 'tRNA'}, 'Location', 'NorthEastOutside');
        end

        function [rnaSyn, rnaExp] = calcExpectedSynthesisExpression(sim)
            import edu.stanford.covert.util.ConstantUtil;

            time = sim.state('Time');
            mass = sim.state('Mass');
            rna = sim.state('Rna');

            exp = rna.expression(rna.matureIndexs);
            mws = rna.molecularWeights(rna.matureIndexs) / ConstantUtil.nAvogadro;
            drs = rna.decayRates(rna.matureIndexs);

            initCnts = exp * mass.cellInitialDryWeight * mass.dryWeightFractionRNA / (exp' * mws);

            rnaExp = initCnts * time.cellCycleLength / log(2);
            rnaSyn = initCnts + initCnts .* drs * time.cellCycleLength / log(2);
        end

        function [rnaSyn, rnaExp] = sampleRNASynthesis(sim, nSample, nIter)
            import edu.stanford.covert.cell.sim.analysis.RNA;

            rna = sim.state('Rna');
            nRna = size(rna.matureIndexs, 1);

            rnaSyn = zeros(nRna, nSample);
            rnaExp = zeros(nRna, nSample);

            for i = 1:nSample
                [rnaSyn(:, i), rnaExp(:, i)] = RNA.simulateRNASynthesis(nIter, i);
            end
        end

        function [rnaSyn, rnaExp] = simulateRNASynthesis(nIter, seed)
            %% references
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], {
                'Metabolism'
                'Transcription'
                'RNAProcessing'
                'RNAModification'
                'tRNAAminoacylation'
                'RNADecay'});

            g = sim.gene;
            rna = sim.state('Rna');
            time = sim.state('Time');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');

            %% seed
            sim.applyOptions('seed', seed);
            sim.seedRandStream();

            for i = 1:numel(sim.states)
                o = sim.states{i};
                o.seed = seed;
                o.seedRandStream();
            end

            for i = 1:numel(sim.processes)
                o = sim.processes{i};
                o.seed = seed;
                o.seedRandStream();
            end

            %% turn off decay
            rna.decayRates(setdiff(1:end, rna.intergenicIndexs)) = 0;

            %% keep track of initial state
            initRNAs = rna.counts;
            initMonomers = pm.counts;
            initComplexs = pc.counts;

            %% simulate
            allRnaSyn = zeros(size(rna.counts));
            allRnaExp = zeros(size(rna.counts));
            for i = 1:nIter
                %mock protein synthesis
                pm.counts = ceil(initMonomers * exp(i * log(2) / time.cellCycleLength));
                pc.counts = ceil(initComplexs * exp(i * log(2) / time.cellCycleLength));

                %remember old RNA counts
                oldRnaCounts = rna.counts;

                %simulate
                sim.evolveState();

                %accounting of new RNAs
                allRnaSyn = allRnaSyn + max(0, rna.counts - oldRnaCounts);

                %accounting of expressed RNAs
                allRnaExp = allRnaExp + rna.counts;
                allRnaExp(rna.matureIndexs(setdiff(1:end, rna.matureMRNAIndexs)), sim.compartment.cytosolIndexs) = ...
                    allRnaExp(rna.matureIndexs(setdiff(1:end, rna.matureMRNAIndexs)), sim.compartment.cytosolIndexs) + ...
                    sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * sum(initComplexs(pc.matureIndexs, :) + initComplexs(pc.boundIndexs, :), 2);
            end

            rnaSyn = allRnaSyn(rna.matureIndexs, sim.compartment.cytosolIndexs);
            aminoacylatedTfs = allRnaSyn(rna.aminoacylatedIndexs, sim.compartment.cytosolIndexs) ~= 0;
            rnaSyn(aminoacylatedTfs) = allRnaSyn(rna.aminoacylatedIndexs(aminoacylatedTfs), sim.compartment.cytosolIndexs) - ...
                initRNAs(rna.matureIndexs(aminoacylatedTfs), sim.compartment.cytosolIndexs);

            rnaExp = ...
                allRnaExp(rna.processedIndexs, sim.compartment.cytosolIndexs) + ...
                allRnaExp(rna.matureIndexs, sim.compartment.cytosolIndexs) + ...
                allRnaExp(rna.aminoacylatedIndexs, sim.compartment.cytosolIndexs) + ...
                allRnaExp(rna.boundIndexs, sim.compartment.cytosolIndexs) + ...
                allRnaExp(rna.misfoldedIndexs, sim.compartment.cytosolIndexs) + ...
                allRnaExp(rna.damagedIndexs, sim.compartment.cytosolIndexs);

            %% scale to cell cycle
            rnaSyn = rnaSyn * 1 / (exp(nIter * log(2) / time.cellCycleLength) - 1);
            rnaExp = rnaExp * 1 / (exp(nIter * log(2) / time.cellCycleLength) - 1);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% Semantic Soft Segmentation
% This function implements the soft segmentation approach described in
% Yagiz Aksoy, Tae-Hyun Oh, Sylvain Paris, Marc Pollefeys, Wojciech Matusik
% "Semantic Soft Segmentation", ACM TOG (Proc. SIGGRAPH) 2018

function [softSegments, initSoftSegments, Laplacian, affinities, features, superpixels, eigenvectors, eigenvalues] = SemanticSoftSegmentation(image, features)

    disp('Semantic Soft Segmentation')
    % Prepare the inputs and superpixels
    image = im2double(image);
    if size(features, 3) > 3 % If the features are raw, hyperdimensional, preprocess them
        features = preprocessFeatures(features, image);
    else
        features = im2double(features);
    end
    superpixels = Superpixels(image);
    [h, w, ~] = size(image);

    disp('     Computing affinities')
    % Compute the affinities and the Laplacian
    affinities{1} = mattingAffinity(image);
    affinities{2} = superpixels.neighborAffinities(features); % semantic affinity
    affinities{3} = superpixels.nearbyAffinities(image); % non-local color affinity
    Laplacian = affinityMatrixToLaplacian(affinities{1} + 0.01 * affinities{2} + 0.01 * affinities{3}); % Equation 6

    disp('     Computing eigenvectors')
    % Compute the eigendecomposition
    eigCnt = 100; % We use 100 eigenvectors in the optimization
    [eigenvectors, eigenvalues] = eigs(Laplacian, eigCnt, 'SM');

    disp('     Initial optimization')
    % Compute initial soft segments
    initialSegmCnt = 40;
    sparsityParam = 0.8;
    iterCnt = 40;
    % feeding features to the function below triggers semantic intialization
    initSoftSegments = softSegmentsFromEigs(eigenvectors, eigenvalues, Laplacian, ...
                                            h, w, features, initialSegmCnt, iterCnt, sparsityParam, [], []);

    % Group segments w.r.t. their mean semantic feature vectors
    groupedSegments = groupSegments(initSoftSegments, features);

    disp('     Final optimization')
    % Do the final sparsification
    softSegments = sparsifySegments(groupedSegments, Laplacian, imageGradient(image, false, 6));

    disp('     Done.')
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ll,lls]=stdtloglik(x,mu,sigma2,nu)
% Log likelihood of the Standardized T distribution
%
% USAGE:
%   [LL,LLS]=stdtloglik(X,MU,SIGMA2,NU)
%
% INPUTS:
%   X      - Standardized T random variables, either scalar or column vector
%   MU     - Mean of X, either scalar or size(x)
%   SIGMA2 - Variance of X, either scalar or size(x)
%   V      - Degree of freedom parameters, either scalar or size(x)
%
% OUTPUTS:
%   LL    - Log-likelihood evaluated at X
%   LLS   - Vector of log-likelihoods corresponding to X
%
% COMMENTS:
%   V>2
%
% REFERENCES:
%   [1] Cassella and Berger (1990) 'Statistical Inference'
%
% See also STDTCDF, STDTINV, STDTRND, STDTPDF

% Copyright:
% Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 9/1/2004

[T,K]=size(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if K~=1
    error('X must be a column vector');
end

if nargin==4
    if length(mu)~=1 && ~all(size(mu)==[T K])
        error('mu must be either a scalar or the same size as X');
    end
    if any(sigma2<=0)
        error('sigma2 must contain only positive elements')
    end
    if length(sigma2)==1
        sigma2=sigma2*ones(T,K);
    elseif size(sigma2,1)~=T || size(sigma2,2)~=1
        error('sigma2 must be a scalar or a vector with the same dimensions as X');
    end
    if length(nu)>1 || nu<=2
        error('nu must be a scalar greater than 2');
    end
    x=x-mu;
else
    error('Only 4 inputs supported');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Compute the log likelihood
ll = T*gammaln(0.5*(nu+1)) - T*gammaln(nu/2) - T/2*log(pi*(nu-2)) ...
    - 0.5*sum(log(sigma2)) - ((nu+1)/2)*sum(log(1 + (x.^2)./(sigma2*(nu-2))));

%Compute the individual log likelihoods if needed
if nargout>1
    lls = gammaln(0.5*(nu+1)) - gammaln(nu/2) - 1/2*log(pi*(nu-2))...
        - 0.5*(log(sigma2)) - ((nu+1)/2)*(log(1 + (x.^2)./(sigma2*(nu-2))));
end

author: bashtage

repo: mfe-toolbox


text:
function gp_init_kron(covfun1, covfun2)

% low-level approach..?

cov1 = gp_cov_pp();
cov2 = gp_cov_noise_iso();
cov = gp_cov_sum(cov1,cov2);

gp_moments()

function [mu, Cov] = gp_moments(y,Ky,Kxy,Kx)

% mu = K_xy * (K_y \ y);
% Cov = K_x - K_xy * inv(K_y) * K_yx

if isnumeric(Kxy)
  Kxy = @(x) Kxy*x;
end
if isnumeric(Ky)
  Ky = @(x) linsolve_cov(Ky,x);
end

mu = Kxy( Ky( y ) );
if nargout >= 2
  if isvector(Kx)
    % Evaluate variances
    Cov = Kx -
  else
    % Evaluate covariance
    Kxy = Kxy(speye(size(Ky1)));
    Cov = Kx - Kxy
  end
end

end

% regular gp:

% set up covariance function
%
% hyperparameter stuff should be removed.. GP is always conditional on
% the given hyperparams?
cov1 = gp_cov_pp(pdf_gamma(2,3), pdf_gamma(0.5,0.5), 'jitter', 1e-6); % struct
gp1 = gp_init(cov1);
cov2 = gp_cov_noise_iso(pdf_gamma(1,1), 'jitter', 0); % struct
gp2 = gp_init(cov2);
% cov = gp_cov_sum_sparse(cov1, cov2);

gp_y = gp_sum(gp1,gp2); % ?

gp_y.set_data_inputs(x);

gp = gp_init(cov)

gp.set_data_inputs(x)

gp.set_hyperparameters(logtheta, ..)
sampler = gp.get_hyperparameter_sampler(..)

gp.get_posterior_moments(logtheta,..)
sampler = gp.get_posterior_sampler(logtheta, ..)
f = sampler(y, ..)
gp.get_hyperparameter_optimizer(..)

% or:
%gp = gp_init(prior_covariance, likelihood_covariance);

%%%%

gp.set_data_inputs(..)

gp.optimize_hyperparameters(..) % ?
gp.draw_hyperparameters(..) %

gp.evaluate_prior_covariance(hyperparameters)
gp.evaluate_likelihood_covariance(noise_covariance)
gp.draw_sample(y); % without data draws from prior

%
% CONSTRUCT THE GP
%

% Create covariance structs
cov1 = gp_init_cov_pp(); % struct
cov1.params.prior = pdf_product(pdf_gamma(2,3), pdf_gamma(0.5,0.5)); % struct
cov1.jitter = 1e-6;
% functions:
% cov1.data_inputs(x) - evaluates, e.g., distance matrix
% cov1.evaluate
% cov1.cholesky()
% cov1.linsolve()
% cov1.multiply()
cov2 = gp_init_cov_pp(); % struct
cov2.params.prior = pdf_product(pdf_gamma(2,3), pdf_gamma(0.5,0.5)); % struct
cov2.jitter = 1e-6;

% Create Kronecker product covariance struct
cov = gp_cov_kron(cov1, cov2);

% Also, there could be, e.g., mean function struct?

gp = gp_init(cov); % ????

%
% ASSIGN DATA
%

gp = gp_set_data_inputs(gp, {x1,x2}); % makes it possible to pre-evaluate,
                                      % e.g., distance matrix

gp_set_hyperparameters() % ??

gp_evaluate_prior_covariance(); % evaluate K and possibly L_K ?

gp_observations(Y,s^2);

% when to evaluate preconditioner?

%
% DO INFERENCE
%

gp_sample_cg(); %% ??

f = gp_sample();

% you could also have functions such as gp_mean, gp_variance for VB
% inference..?i

N1 = [];
N2 = [];

K1 = [];
K2 = [];

LD1 = [];
LD2 = [];

L1 = [];
L2 = [];

N1 = 500;
x1 = 1:N1;
K1 = gp_cov_pp(log(10), x1, x1);

N2 = 500;
x2 = 1:N2;
K2 = gp_cov_pp(log(10), x2, x2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function set_parameters(logt1,logt2,x1,x2)
  logtheta1 = logt1;
  X1 = x1;
  N1 = size(x1,2);

  logtheta2 = logt2;
  X2 = x2;
  N2 = size(x2,2);

  evaluate_covariance();
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function rndfunc = get_posterior_sampling_function(Y,noise_std)

  s = noise_std;

  Afun = @(x) (kronprod(K1,K2,x,'vector')+s^2*x);
  Mfun = get_preconditioner();

  rndfunc = @() draw_sample(Y,Afun,Mfun);

  end

  function F = draw_sample()

  Z0 = randn(N1,N2);
  Z = randn(N1,N2);
  F0 = kronprod(L2,L1,Z0);
  G = F0 + s*Z - Y;
  X = reshape( pcg(Afun, G(:), 1e-4, 20, Mfun), N1, N2 );
  F = F0 - kronprod(K1,K2,X);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function precondfun = get_preconditioner()
  LD_C1 = ldlchol(K1+s*speye(N1));
  LD_C2 = ldlchol(K1+s*speye(N1));
  precondfun = @(x) kronsolve_for_cg2(LD_C1, LD_C2, x);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function evaluate_covariance()
  K1 = covfun1(logtheta1, X1, X1);
  LD1 = ldlchol(K1);
  [L,D] = ldlsplit(LD1);
  L1 = L*sqrt(D);

  K2 = covfun2(logtheta2, X2, X2);
  LD2 = ldlchol(K2);
  [L,D] = ldlsplit(LD2);
  L2 = L*sqrt(D);
  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function meshnd_example
%MESHND_EXAMPLE example usage of meshnd and meshsparse.
%
% Example:
%   meshnd_example
%
% See also meshnd.

% Copyright 2007, Timothy A. Davis, Univ. of Florida

help meshnd

% 2D mesh, compare with Cleve Moler's demos

m = 7 ;
n = 7 ;

[G p pinv Gnew] = meshnd (m,n) ;
fprintf ('Original mesh:\n') ;
disp (G) ;
fprintf ('Permuted node numbers using meshnd.m (nested dissection):\n') ;
disp (Gnew) ;

Moler = nested (n+2) ;
Moler = Moler (2:n+1,2:n+1) ;
fprintf ('Cleve Moler''s nested dissection ordering, using nested.m\n') ;
disp (Moler) ;
fprintf ('Difference between nested.m and meshnd.m:\n') ;
disp (Gnew-Moler) ;

% 2D and 3D meshes

stencils = [5 9 7 27] ;
mm = [7 7 7 7] ;
nn = [7 7 7 7] ;
kk = [1 1 7 7] ;

for s = 1:4

    m = mm (s) ;
    n = nn (s) ;
    k = kk (s) ;
    [G p] = meshnd (mm (s), nn (s), kk (s)) ;
    A = meshsparse (G, stencils (s)) ;
    C = A (p,p) ;
    parent = etree (C) ;
    try
        L = chol (C, 'lower') ;
    catch
        % old version of MATLAB
        L = chol (C)' ;
    end
    subplot (4,5,(s-1)*5 + 1) ;
    do_spy (A) ;
    if (k > 1)
	title (sprintf ('%d-by-%d-by-%d mesh, %d-point stencil', ...
	    m, n, k, stencils (s))) ;
    else
	title (sprintf ('%d-by-%d mesh, %d-point stencil', ...
	    m, n, stencils (s))) ;
    end
    subplot (4,5,(s-1)*5 + 2) ;
    do_spy (C) ;
    title ('nested dissection') ;
    subplot (4,5,(s-1)*5 + 3) ;
    treeplot (parent) ;
    title ('etree') ;
    xlabel ('') ;
    subplot (4,5,(s-1)*5 + 4) ;
    do_spy (L) ;
    title (sprintf ('Cholesky with nd, nnz %d', nnz (L))) ;
    try
        % use the built-in AMD
        p = amd (A) ;
    catch
        try
            % use AMD from SuiteSparse
            p = amd2 (A) ;
        catch
            % use the older built-in SYMAMD
            p = symamd (A) ;
        end
    end
    try
        L = chol (A (p,p), 'lower') ;
    catch
        % old version of MATLAB
        L = chol (A (p,p))' ;
    end
    subplot (4,5,(s-1)*5 + 5) ;
    do_spy (L) ;
    title (sprintf ('Cholesky with amd, nnz %d', nnz (L))) ;

end

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

function do_spy (A)
%DO_SPY use cspy(A) to plot a matrix, or spy(A) if cspy not installed.
try
    % This function is in CSparse.  It generates better looking plots than spy.
    cspy (A) ;
catch
    spy (A) ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%   SCRIPT TO COMPUTE THE TORQUES AT EACH JOINT FOR DIFFERENT MOTION STATES OF
%   a 3 DOF robotic arm.
%
%   In order to select an actuator (i.e. an electric motor, brushless for
%   example), we may need.
%   Nominal torque and speed: torque and speed for the 80% of the use of the motor.
%   Peak torque and speed: torque and speed for short periods of time. That
%   is: a higher torque that can be exerted at higher speeds during a
%   maximum of a 20% of the time.
%
%   Of course, the 80-20% are just bare numbers and should be given by the
%   robot manufacturer. Actually the torque in any motor depends directly
%   on the quantity of current that can be driven into the motor.
%   Typically, the peak torque is associated with a peak current. If the
%   peak current is maintained for a long time the motor will not be able
%   to dissipate the heat inside, thus generating high temperatures that
%   could melt the coils, conductors... etc.
%
%   The script uses the inverse dynamic model of the robot to simulate
%   different motion states and compute the torques for each situation.
%   The torques at each joint, as well as the torques at each motor are computed
%   (reduced by the gear ratio). A trapezoidal speed profile can be
%   selected at the parameters section of this file. This trapezoidal
%   profile should meet the desired features of the robot, such as maximum
%   joint speed, maximum acceleration/deceleration.
%   Please note that the movement of the robot is not simulated. The reader
%   should imagine that the manipulator is fixed at a given initial
%   position and the inverse dynamic function returns the torques at each
%   joint necessary to bring the arm to that motion state (defined by position q,
%   speed qd and acceleration qdd).
%
%   As a reference, below you can find the call to the inverse dynamics
%   function.
%
%   TAU = inversedynamic(robot, Q, QD, QDD, GRAV, FEXT)
%
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

function motor_selection_3dofplanar

close all
global robot
robot.dynamics.friction = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS SECTION
%   Feel free to change the values of q, maximum_speeds and
%   maximum_accels.
%
%   This script tries to allow the student to test any mechanism
%   at the worst case. In this sense, q should be adjusted
%   as the pose where each joint would (statically) be needing a
%   higher torque.
%   The maximum_speeds and maximum_acceleration define a trapezoidal
%   speed profile. This trapezoidal speed is used by most machines
%   to command changes in speed in any of their joints.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% robot pose: experiment by changing the pose while observing the different
%             torques at each joint
q=[0 0 0]; %rad
%maximum speeds at for each joint
maximum_speeds=[pi pi pi];%rad/second
%maximum acceleration/deceleration for each joint
maximum_accels=[pi/4 pi/4 pi/4]; %rad/second^2

% time of the trapezoidal profile that the joint moves at maximum speed
time_at_constant_speed=2; %seconds

%load robot parameters. Just uncomment this line
robot=load_robot('example', '3dofplanar');
drawrobot3d(robot, q)

% CUIDADO: los resultados se calculan con estas ratios.
robot.motors.G = [1 1 1]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  FIRST, COMPUTE TRAPEZOIDAL PROFILES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%compute acceleration plus deceleration times for every joint
time_acc = 2*maximum_speeds./maximum_accels+time_at_constant_speed;

%compute the total time for the slowest joint
total_time=max(time_acc);

% Trapezoidal speed profiles for each joint
[input_speeds, input_accels, time]=build_trapezoidal_speed_profile(maximum_speeds, maximum_accels, total_time);

% EJERCICIO:
% REALIZA UN PLOT DE input_speeds e input_accels.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FINALLY, COMPUTE TORQUES FOR EACH MOTION STATE.
% Please note that we consider that the robot is placed at a fixed position and consider
% different motion situations when we change the acceleration and speed at
% each joint. For each motion state, the inverse dynamic model returns the
% torques at each joint that would bring the robot to that motion.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
compute_inverse_dynamics(q, input_speeds, input_accels, time);

%Computes a trapezoidal speed profile for every joint given maximum
%permitted accelerations and maximum joint speeds
function [input_speeds, input_accelerations, time]=build_trapezoidal_speed_profile(maximum_speeds, maximum_accels, total_time)

delta_time=0.01;

%build time vector: twice acceleration time plus time at constant speed
time = 0:delta_time:total_time;

input_speeds=[];
input_accelerations=[];

for j=1:length(maximum_speeds),
    vel_row=[];
    acc_row=[];
    for i=1:length(time),
        [vel acc] = compute_values(time(i), maximum_speeds(j), maximum_accels(j), total_time);
        vel_row = [vel_row vel];
        acc_row = [acc_row acc];
    end
    input_speeds = [input_speeds; vel_row];
    input_accelerations = [input_accelerations; acc_row];
end

%returns the values of velocity and speed corresponding to a given time
function [vel acc]=compute_values(time_i, vel_max, acc_max, total_time)

tacc = vel_max/acc_max;
tdec = total_time-tacc;

if time_i < tacc
    vel = time_i.*acc_max;
    acc = acc_max;
    return;
elseif (time_i >= tacc) & (time_i < tdec)
    vel = vel_max;
    acc = 0;
    return;
else % time_i> tdec
    vel = vel_max-(time_i-tdec)*acc_max;
    acc = -acc_max;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   COMPUTE THE INVERSE DYNAMICS FOR EACH MOTION STATE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function compute_inverse_dynamics(q, input_speeds, input_accels, time)
global robot

%adjust_view(robot)
torques=[];
for j=1:length(time)
    fprintf('\nComputing time %d out of %d', j, length(time));
    % compute the torque to bring the robot instantaneously to this motion
    % state. change M=1  to add the effects of a 1kg mass load at the end effector
    M=0.1;
    %please note that the force due to the load acts on the z axis of
    tau=inversedynamic(robot, q, input_speeds(:,j), input_accels(:,j), [0  -9.81 0]', [0 -M*9.81 0 0 0 0]');
    torques=[torques tau];
end

%plot trapezoidal profiles
figure, hold, xlabel('time (s)'), ylabel('Input reference speeds (rad/s)')
plot(time, input_speeds(1,:), time, input_speeds(2,:), time, input_speeds(3,:) );
legend('Speed for joint 1 (qd1)','Speed for joint 2 (qd2)', 'Speed for joint 3 (qd3)')
%plot trapezoidal profiles, acceleration
figure, hold, xlabel('time (s)'), ylabel('Input reference acceleration (rad/s)')
plot(time, input_accels(1,:), time, input_accels(2,:), time, input_accels(3,:));
legend('Acceleration for joint 1 (qd1)','Acceleration for joint 2 (qd2)', 'Acceleration for joint 3 (qd3)')

% plot results. First, torques at each joint
figure, hold, xlabel('time (s)'), ylabel('Join Torques (N m)')
plot(time, torques(1,:), time, torques(2,:));
legend('Torque for joint 1 ','Torque  for joint 2 ')

% OTRA INFORMACION QUE ES ESENCIAL QUE PLOTEES:
% - Los pares en cada motor (utiliza robot.motors.G(1), .G(2)...

% - La potencia instantánea que realiza cada motor (la potencia es tau*w en
% cada articulación). La potencia pico del motor nos dará una idea de qué
% motor buscar en el catálogo.

% - Plotea la velocidad de cada motor en rpm: en efecto, los fabricantes
% tienen la costumbre de expresar la velocidad del motor en rpm y no en
% rad/s. Los fabricantes de robots tienen la costumbre de expresar las
% velocidades articulares en grados/s. Se te sugiere que plotees las
% velocidades de cada motor en rpm.

% Finalmente, en base al resultado, de cada motor, expresa:
% - El par pico (valor máximo en valor absoluto del par realizado por el motor).
% - El par nominal. Es el par en el punto intermedio del resultado del vector torques.
% - Velocidad nominal del motor en rpm.
% - La velocidad máxima en rpm del motor.

repo: ARTE


text:
function env_release_cluster
% Release any currently acquired cluster resources.
%
% This stops the heartbeat signal; it is the workers' responsibility to clean themselves up when
% they are no longer needed.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-04-12

disp('Releasing cluster resources...');
global tracking;
% this will lead to the deletion of any involved heartbeat timer.
tracking.cluster_requested = [];
par_globalsetting('engine','local');
if tracking.autokill_workers
    par_clearworkers([],true); end
disp('Set default compute scheduling to local.');

repo: REST


text:
function c = roots_to_i4poly ( n, x )

%*****************************************************************************80
%
%% ROOTS_TO_I4POLY converts polynomial roots to polynomial coefficients.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of roots specified.
%
%    Input, integer X(N), the roots.
%
%    Output, integer C(1:N+1), the coefficients of the polynomial.
%

%
%  Initialize C to (0, 0, ..., 0, 1).
%  Essentially, we are setting up a divided difference table.
%
  c(1:n) = 0;
  c(n+1) = 1;
%
%  Convert to standard polynomial form by shifting the abscissas
%  of the divided difference table to 0.
%
  for j = 1 : n
    for i = 1 : n+1-j
      c(n-i) = c(n-i) - x(n+1-i-j+1) * c(n-i+1);
    end
  end

  return
end


text:
% IS_VALID_STRUCT2 - Determine if the input structure is "valid"
%
%     IS_VALID_STRUCT2 determines whether or not the input
%     structure contains additional field data beyond the
%     minimum valid fields checked for by  IS_VALID_STRUCT.
%     Recall that a valid FEM_GRID_STRUCT contains (atleast)
%     the following (NOT EMPTY) fields:
%
%       .name,.e,.x,.y,.z,.bnd
%
%     IS_VALID_STRUCT2 verifies the existence of .A, .B, .A0,
%     and .T which are filled by BELINT, and of .ar filled by
%     EL_AREAS.  These additional FEM fields are needed by
%     OPNML/MATLAB element-finding, basis  and interpotion
%     routines.
%
% CALL: errflag=is_valid_struct2(fem_grid_struct)
%
% Written by : Brian O. Blanton
% Summer 1998
%
function errflag=is_valid_struct2(fem_grid_struct)

errflag=0;

% Make sure input argument is actually a structure
%
if ~isstruct(fem_grid_struct)
   disp('    Argument to IS_VALID_STRUCT2 must be a structure.');return
end

% now, make sure the structure contains the additional fields,
% as above
if ~isfield(fem_grid_struct,'A')
   disp('    "A" field not part of fem_grid_struct; run BELINT');return
elseif ~isfield(fem_grid_struct,'B')
   disp('    "B" field not part of fem_grid_struct; run BELINT');return
elseif ~isfield(fem_grid_struct,'A0')
   disp('    "A0" field not part of fem_grid_struct; run BELINT');return
elseif ~isfield(fem_grid_struct,'T')
   disp('    "T" field not part of fem_grid_struct; run BELINT');return
elseif ~isfield(fem_grid_struct,'ar')
   disp('    "ar" field not part of fem_grid_struct; run ');return
end

% now, make sure these additional fields are NOT EMPTY
%
if isempty(fem_grid_struct.A)
   disp('    "A" field in fem_grid_struct is EMPTY');return
elseif isempty(fem_grid_struct.B)
   disp('    "B" field in fem_grid_struct is EMPTY');return
elseif isempty(fem_grid_struct.A0)
   disp('    "A0" field in fem_grid_struct is EMPTY');return
elseif isempty(fem_grid_struct.T)
   disp('    "T" field in fem_grid_struct is EMPTY');return
elseif isempty(fem_grid_struct.ar)
   disp('    "ar" field in fem_grid_struct is EMPTY');return
end

errflag=1;

%
%        Brian O. Blanton
%        Curr. in Marine Science
%        15-1A Venable Hall
%        CB# 3300
%        Uni. of North Carolina
%        Chapel Hill, NC
%                 27599-3300
%
%        919-962-4466
%        blanton@marine.unc.edu
%
%        Summer 1998

author: CHLNDDEV


text:
function outstruct=memmapstream(bytes, format, varargin)
%
%    outstruct=memmapstream(bytes, format)
%
%    Map a byte-array (in char array or uint8/int8 array) into a structure
%    using a dictionary (format is compatible with memmapfile in MATLAB)
%
%    This function is compatible with both MATLAB and GNU Octave.
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        bytes: a char, int8 or uint8 vector or array
%        format: a 3-column cell array in the format compatible with the
%              'Format' parameter of memmapfile in MATLAB. It has the
%              following structure
%
%             column 1: data type string, it can be one of the following
%                'int8','int16','int32','int64',
%                'uint8','uint16','uint32','uint64',
%                'single','double','logical'
%             column 2: an integer vector denoting the size of the data
%             column 3: a string denoting the fieldname in the output struct
%
%             For example format={'int8',[1,8],'key'; 'float',[1,1],'value'}
%             reads the first 8 bytes from 'bytes' as the first subfield
%             'key' and the following 4 bytes as the floating point 'value'
%             subfield.
%
%    output:
%        outstruct: a structure containing the required field
%
%    example:
%        bytestream=['Andy' 5 'JT'];
%        format={'uint8', [1,4], 'name',
%              'uint8', [1,1], 'age',
%              'uint8', [1,2], 'school'};
%        data=memmapstream(bytestream,format);
%
%    this file is part of JNIfTI specification: https://github.com/fangq/jnifti
%
%    License: Apache 2.0, see https://github.com/fangq/jnifti for details
%

if(nargin<2)
   error('must provide bytes and format as inputs');
end

if(~ischar(bytes) && ~isa(bytes,'int8') && ~isa(bytes,'uint8') || isempty(bytes))
   error('first input, bytes, must be a char-array or uint8/int8 vector');
end

if(~iscell(format) || size(format,2)<3 || size(format,1)==0 || ~ischar(format{1,1}))
   error('second input, format, must be a 3-column cell array, in a format described by the memmapfile Format field.');
end

bytes=bytes(:)';

datatype=struct('int8',1,'int16',2,'int32',4,'int64',8,'uint8',1,'uint16',2,'uint32',4,'uint64',8,'single',4,'double',8);

opt=varargin2struct(varargin{:});
opt.usemap=jsonopt('usemap',0,opt) && exist('containers.Map');

if(opt.usemap)
    outstruct=containers.Map();
else
    outstruct=struct();
end
len=1;
for i=1:size(format,1)
    bytelen=datatype.(format{i,1})*prod(format{i,2});
    if(opt.usemap)
        outstruct(format{i,3})=reshape(typecast(uint8(bytes(len:bytelen+len-1)),format{i,1}),format{i,2});
    else
        outstruct.(format{i,3})=reshape(typecast(uint8(bytes(len:bytelen+len-1)),format{i,1}),format{i,2});
    end
    len=len+bytelen;
    if(len>length(bytes))
        break;
    end
end


text:
function seg = segBuildMesh(seg, savePath, varargin);
%
% Build a new mesh structure for a mrVista2 segmentation object.
%
% seg = segBuildMesh(seg, [savePath=dialog], [options]);
%
% INPUTS:
%
% seg: segmentation structure (see segCreate).
%
% savePath: optional path to save the mesh; if entered as 'dialog', will
% pop up a user dialog to choose the save location.
%
% Optional input arguments:
%   'visualize', [val]: set visualize flag to 1 or 0 [default 1]. If 1,
%                       will open up the new mesh in mrMesh.
%   'numGrayLayers', [val]: set the # of gray layers to grow. [Default 0]
%   'meshName', [name]: set the name of the meah. [Default: same name as
%                       segmentation, plus a number if there are already
%                       other meshes attached]
%
%
%
% OUTPUTS:
%
% seg: modified segmentation with a new mesh structure attached to the
% seg.mesh{} field.
%
% ras, 10/2006.
if notDefined('savePath'), savePath = ''; end

%% params / defaults
visualize = 1;
numGrayLayers = 0;
meshName = seg.name;
if ~isempty(seg.mesh), meshName = [meshName num2str(length(seg.mesh)+1)]; end

%% parse options
for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'visualize', visualize = varargin{i+1};
        case {'numgraylayers' 'nlayers'}, numGrayLayers = varargin{i+1};
        case {'name' 'meshname'}, meshName = varargin{i+1};
    end
end

%% get the WM classification
C = segGet(seg, 'classification');

% Create a volume which designates the white matter locations
wm = uint8( zeros(C.header.xsize, C.header.ysize, C.header.zsize) );
xx = C.header.voi(1):C.header.voi(2);
yy = C.header.voi(3):C.header.voi(4);
zz = C.header.voi(5):C.header.voi(6);
wm(xx,yy,zz) = wm(xx,yy,zz) | (C.data == C.type.white);

% add gray layers if specified
if numGrayLayers > 0
    [nodes edges C] = mrgGrowGray(C, numGrayLayers);
    wm = uint8( (C.data(xx,yy,zz) == C.type.white) | ...
                (C.data(xx,yy,zz) == C.type.gray) );
end

%% build the new mesh
[newMesh tenseMesh] = mrmBuild(wm, seg.voxelSize, 1, visualize);

%% attach the 'tense' (unsmoothed) mesh to the new mesh as the
%% 'initVertices' field (very useful to hold on to)
initVertices = meshGet(tenseMesh, 'vertices');
newMesh = meshSet(newMesh, 'initialvertices', initVertices);
if isempty(seg.nodes) | isempty(seg.edges)
    [seg.nodes seg.edges] = segGet(seg, 'gray');
end
v2g = mrmMapVerticesToGray(initVertices, seg.nodes, seg.voxelSize, seg.edges);

newMesh = meshSet(newMesh, 'vertexGrayMap', v2g);
newMesh = meshSet(newMesh, 'name', meshName);
newMesh = meshSet(newMesh, 'nGrayLayers', numGrayLayers);

%% attach new mesh to segmentation
seg.mesh{end+1} = newMesh;

% set the new mesh as the selected one
seg.settings.mesh = length(seg.mesh);

%% save if selected
if ~isempty(savePath)
    if isequal(lower(savePath), 'dialog')
        savePath = mrSelectDataFile('stayput', 'w', '*.mat',...
                                    'Save Mesh File As...');
        if isempty(savePath)
            disp('Not Saving Mesh')
            return
        end
    end

    try
        mrmWriteMeshFile(newMesh, savePath);
    catch
        msg = sprintf(['Couldn''t save file %s. Not erroring, but you ' ...
                       'should try saving again.'], savePath);
        myWarnDlg(msg);
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function rot = byMatrix(M,varargin)
% define rotations by matrices
%
% Syntax
%   M = eye(3)
%   rot = rotation.byMatrix(M)
%
% Input
%  M - 3x3 rotation matrix
%
% Output
%  rot - @rotation
%
% See also
% rotation/rotentation rotation/byEuler rotation/byAxisAngle

% negative determinant indicates improper rotation
isInv = false(size(M,3),1);
for i = 1:size(M,3)
  isInv(i) = det(M(:,:,i))<0;
end

% ensure M is proper
M(:,:,isInv) = -M(:,:,isInv);

rot = rotation(mat2quat(M));
rot.i = isInv;


text:
function img = imresizecrop(img, M, METHOD)
%
% img = imresizecrop(img, M, METHOD);
%
% Output an image of size M(1) x M(2).
%
% URL: http://labelme.csail.mit.edu/Release3.0/browserTools/php/matlab_toolbox.php

if nargin < 3
    METHOD = 'bilinear';
end

if length(M) == 1
    M = [M(1) M(1)];
end

scaling = max([M(1)/size(img,1) M(2)/size(img,2)]);

%scaling = M/min([size(img,1) size(img,2)]);

newsize = round([size(img,1) size(img,2)]*scaling);
img = imresize(img, newsize, METHOD);

[nr nc cc] = size(img);

sr = floor((nr-M(1))/2);
sc = floor((nc-M(2))/2);

img = img(sr+1:sr+M(1), sc+1:sc+M(2),:);

author: adikhosla

repo: feature-extraction


text:
function [nm] = A2nm(A)
% Convert length from angstroms to nanometers.
% Chad A. Greene 2012
nm = A*0.1;


text:
function ret=writesdpa(fname,A,b,c,K,pars)
%WRITESDPA Take a problem in SeDuMi format and writes it in SDPA format
%
%  Usage:
%
%  ret=writesdpa(fname,A,b,c,K,pars)
%
%      fname           Name of SDPpack file, in quotes
%      A,b,c,K         Problem in SeDuMi form
%      pars            Optional parameters.
%                      pars.printlevel=0           No printed output
%                      pars.printlevel=1 (default) Some printed output.
%
%      ret             ret=0 on success, ret=1 on failure.
%
%  Notes:
%
%     Problems with complex data are not allowed.
%
%     Quadratic cone and rotated cone constraints are not supported.
%
%     Nonsymmetric A.s and C.s matrices are symmetrized with A=(A+A')/2
%     a warning is given when this happens.
%
%     Free variables are not supported.
%
%     Floating point numbers are written out with 18 decimal digits for
%     accuracy.
%
%  Please contact the author (Brian Borchers, borchers@nmt.edu) with any
%  questions or bug reports.
%
%  Third Version: 3/3/06.  Corrected a bug in the handling of
%                          nonsymmetric constraint matrices.
%
%  Second Version: 7/14/04.  Modified to vastly speed up operations on sparse
%                 matrices.  On some problems, this version is 100 times
%                 faster!
%
%  First Version: 7/14/03.  Modified from old writesdp.m which wrote problems
%                 in SDPPack format.
%
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

%
% First, check to see whether or not we should be quiet.
%
if (nargin > 5)
  if (isfield(pars,'printlevel'))
    if (pars.printlevel == 0)
      quiet=1;
    else
      quiet=0;
    end
  else
    quiet=0;
  end
else
  pars.printlevel=1;
  quiet=0;
end
%
%  First, check for complex numbers in A, b, or c.
%
if (isreal(A) ~= 1)
  if (quiet == 0)
    fprintf('A is not real!\n');
  end
  ret=1;
  return
end
if (isreal(b) ~= 1)
  if (quiet == 0)
    fprintf('b is not real!\n');
  end
  ret=1;
  return
end
if (isreal(c) ~= 1)
  if (quiet == 0)
    fprintf('c is not real!\n');
  end
  ret=1;
  return
end
%
%  Check for any quadratic cone constraints.
%
if (isfield(K,'q'))
  if ((~isempty(K.q)) & (K.q ~= 0))
    if (quiet == 0)
      fprintf('quadratic cone constraints are not supported.\n');
    end
    ret=1;
    return
  end
end
%
%  Check for any rotated cone constraints.
%
if (isfield(K,'r'))
  if ((~isempty(K.r)) & (K.r ~= 0))
    if (quiet == 0)
      fprintf('rotated cone constraints are not supported.\n');
    end
    ret=1;
    return
  end
end
%
% Check for any free variables.
%
if (isfield(K,'f'))
  if ((~isempty(K.f)) & (K.f ~= 0))
    if (quiet == 0)
      fprintf('Free variables are not supported.\n');
    end
    ret=1;
    return
  end
end
%
%  Find the number of constraints.
%
m=length(b);
%
%  Deal with the following special case.  If A is transposed, transpose
%  it again so that it is of the right size.
%
[Am,An]=size(A);
if (Am ~= m)
  if (An == m)
    if (quiet==0)
      fprintf('Transposing A to match b \n');
    end
    AT=A;
    A=A';
%
% Also swap Am and An so that they're correct.
%
    temp=Am;
    Am=An;
    An=temp;
  else
%
% In this case, A is just plain the wrong size.
%
    if (quiet==0)
      fprintf('A is not of the correct size to match b \n');
    end
    ret=1;
    return
  end
else
%
% No need to transpose A, but we'll need AT.
%
  AT=A';
end
%
%  Deal with the following special case:  if c==0, then c should really
%  be a zero vector of the appropriate size.
%
if (c == 0)
  if (quiet==0)
    fprintf('Expanding c to the appropriate size\n');
  end
  c=sparse(An,1);
end
%
% If c is empty, then act as if it was zero.
%
if (isempty(c))
  if (quiet==0)
    fprintf('Expanding empty c to zeros of the appropriate size\n');
  end
  c=sparse(An,1);
end
%
% If c is a row vector, make it a column vector.
%
[cm,cn]=size(c);
if (cn > cm)
  c=c';
end
%
%  Get the size data.
%
%
% First, the size of the LP block.
%
if (isfield(K,'l'))
  nlin=K.l;
  sizelin=nlin;
  if (isempty(sizelin))
    sizelin=0;
    nlin=0;
  end
  if (K.l == 0)
    nlin=0;
    sizelin=0;
  end
else
  nlin=0;
  sizelin=0;
end
%
% Get the sizes of the SDP blocks.
%
if (isfield(K,'s'))
  nsdpblocks=length(K.s);
  sizesdp=sum((K.s).^2);
  if (isempty(sizesdp))
    sizesdp=0;
    nsdpblocks=0;
  end
  if (K.s == 0)
    nsdpblocks=0;
    sizesdp=0;
  end
else
  sizesdp=0;
  nsdpblocks=0;
end
%
% Figure out the number of blocks.
%
nblocks=nsdpblocks;
if (nlin>0)
  nblocks=nblocks+1;
end
%
%  print out some size information
%
if (quiet==0)
  fprintf('Number of constraints: %d \n',m);
  fprintf('Number of SDP blocks: %d \n',nsdpblocks);
  fprintf('Number of LP vars: %d \n',nlin);
end
%
%  Open up the file for writing.
%
fid=fopen(fname,'w');
if (fid==-1)
  if (quiet==0)
    fprintf('Could not open file for output!');
  end
  ret=1;
  return
end
%
%  Print out m, the number of constraints.
%
fprintf(fid,'%d \n',m);
%
% Next, the number of blocks.
%
fprintf(fid,'%d \n',nblocks);
%
% Print out the block structure.
%
if (K.s > 0)
  fprintf(fid,'%d ',K.s);
end
if (nlin > 0)
  fprintf(fid,'%d ',-nlin);
end
fprintf(fid,'\n');
%
%  Next, b, with all on one line.
%
fprintf(fid,'%.18e  ',full(b));
fprintf(fid,'\n');
%
% First, the C matrix.
%
%
%  First, calculate where in c things start.
%
base=sizelin+1;
%
%  Next, work through the SDP blocks.
%
for i=1:nsdpblocks
%
% Get the current block of C.
%
  I=find(c);
  II=find(I>=base);
  I=I(II);
  II=find(I<=base+K.s(i)^2-1);
  I=I(II);
  II=I-(base-1)*ones(size(I));
  work=sparse(II,ones(size(II)),c(I),K.s(i)^2,1);
  work=reshape(work,K.s(i),K.s(i));

%
% Check this block for symmetry.
%
if (norm(work-work','fro') ~= 0)
  if (quiet==0)
    fprintf('Non symmetric C.s matrix!\n');
  end
  work=(work+work')/2;
end
%
% Write out the C.s matrix.
%
  work=triu(work);
  [II,JJ,V]=find(work);
  cnt=length(II);
  if (cnt ~= 0)
    fprintf(fid,'%d %d %d %d %.18e\n',[zeros(size(II)) i*ones(size(II)) II JJ -V]');
  end

%
%  Next, update to the next base.
%
  base=base+K.s(i)^2;
end
%
% Print out the coefficients for the linear block of C.
%
for i=1:nlin
  if (c(i) ~= 0.0)
    fprintf(fid,'%d %d %d %d %.18e\n',[0 nsdpblocks+1 i i -full(c(i))]);
  end
end

%
%  Now, loop through the constraints, one at a time.
%
for cn=1:m
%
%  Print out the SDP part of constraint cn.
%
  base=sizelin+1;
  rowcn=AT(:,cn);
  for i=1:nsdpblocks
    I=find(rowcn);
    II=find(I>=base);
    I=I(II);
    II=find(I<=(base+K.s(i)^2-1));
    I=I(II);
    II=I-(base-1)*ones(size(I));
    work=sparse(II,ones(size(II)),rowcn(I),K.s(i)^2,1);

    work=reshape(work,K.s(i),K.s(i));

    if (norm(work-work','fro') ~= 0)
      if (quiet==0)
        fprintf('Non symmetric A.s matrix! \n');
      end
      work=(work+work')/2;
    end
%
% Ignore the lower left triangle.
%
    work=triu(work);
%
% Get the nonzero entries.
%
    [II,JJ,V]=find(work);
    cnt=length(II);
    if (cnt ~= 0)
      fprintf(fid,'%d %d %d %d %.18e\n',[cn*ones(size(II)) i*ones(size(II)) II JJ V]');
    end
%
%  Next, update to the next base.
%
    base=base+K.s(i)^2;
  end
%
% Finally, the linear part.
%
  I=find(rowcn);
  II=find(I<=nlin);
  I=I(II);
  workrow=sparse(I,ones(size(I)),rowcn(I),nlin,1);
  [II,JJ,V]=find(workrow);
  if (length(II) > 0)
    fprintf(fid,'%d %d %d %d %.18e\n',[cn*ones(length(II),1) (nsdpblocks+1)*ones(length(II),1) II II V]');
  end
end
%
% Close the file.
%
fclose(fid);
%
% Return success
%
ret=0;
return

repo: cluster


text:
function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

h = sigmoid(X*theta);
% J = (1/m)*sum(-y .* log(h) - (1 - y) .* log(1-h));
shift_theta = theta(2:size(theta));
theta_reg = [0;shift_theta];

J = (1/m)*(-y'* log(h) - (1 - y)'*log(1-h))+(lambda/(2*m))*theta_reg'*theta_reg;

% grad_zero = (1/m)*X(:,1)'*(h-y);
% grad_rest = (1/m)*(shift_x'*(h - y)+lambda*shift_theta);
% grad      = cat(1, grad_zero, grad_rest);

grad = (1/m)*(X'*(h-y)+lambda*theta_reg);

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

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

H = sigmoid(X*theta);
T = y.*log(H) + (1 - y).*log(1 - H);
J = -1/m*sum(T) + lambda/(2*m)*sum(theta(2:end).^2);

ta = [0; theta(2:end)];
grad = X'*(H - y)/m + lambda/m*ta;

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
clear all; close all; clc;
%% Auxiliary Classifier Generative Adversarial Network
%% Load Data
load('mnistAll.mat')
trainX = preprocess(mnist.train_images);
trainY = mnist.train_labels;
testX = preprocess(mnist.test_images);
testY = mnist.test_labels;
%% Settings
settings.latentDim = 100; settings.num_labels = 10;
settings.batch_size = 32; settings.image_size = [28,28,1];
settings.lrD = 0.0002; settings.lrG = 0.0002; settings.beta1 = 0.5;
settings.beta2 = 0.999; settings.maxepochs = 50;

%% Initialization
%% Generator
paramsGen.EMW1 = dlarray(...
    initializeGaussian([settings.latentDim,...
    settings.num_labels]));
paramsGen.FCW1 = dlarray(initializeGaussian([128*7*7,...
    settings.latentDim]));
paramsGen.FCb1 = dlarray(zeros(128*7*7,1,'single'));
paramsGen.BNo1 = dlarray(zeros(128,1,'single'));
paramsGen.BNs1 = dlarray(ones(128,1,'single'));
paramsGen.TCW1 = dlarray(initializeGaussian([3,3,128,128]));
paramsGen.TCb1 = dlarray(zeros(128,1,'single'));
paramsGen.BNo2 = dlarray(zeros(128,1,'single'));
paramsGen.BNs2 = dlarray(ones(128,1,'single'));
paramsGen.TCW2 = dlarray(initializeGaussian([3,3,64,128]));
paramsGen.TCb2 = dlarray(zeros(64,1,'single'));
paramsGen.BNo3 = dlarray(zeros(64,1,'single'));
paramsGen.BNs3 = dlarray(ones(64,1,'single'));
paramsGen.CNW1 = dlarray(initializeGaussian([3,3,64,1]));
paramsGen.CNb1 = dlarray(zeros(1,1,'single'));
stGen.BN1 = []; stGen.BN2 = []; stGen.BN3 = [];

%% Discriminator
paramsDis.CNW1 = dlarray(initializeGaussian([3,3,1,16]));
paramsDis.CNb1 = dlarray(zeros(16,1,'single'));
paramsDis.CNW2 = dlarray(initializeGaussian([3,3,16,32]));
paramsDis.CNb2 = dlarray(zeros(32,1,'single'));
paramsDis.BNo1 = dlarray(zeros(32,1,'single'));
paramsDis.BNs1 = dlarray(ones(32,1,'single'));

paramsDis.CNW3 = dlarray(initializeGaussian([3,3,32,64]));
paramsDis.CNb3 = dlarray(zeros(64,1,'single'));
paramsDis.BNo2 = dlarray(zeros(64,1,'single'));
paramsDis.BNs2 = dlarray(ones(64,1,'single'));

paramsDis.CNW4 = dlarray(initializeGaussian([3,3,64,128]));
paramsDis.CNb4 = dlarray(zeros(128,1,'single'));

paramsDis.FCW1 = dlarray(initializeGaussian([1+settings.num_labels,128*4*4]));
paramsDis.FCb1 = dlarray(zeros(1+settings.num_labels,1,'single'));
stDis.BN1 = []; stDis.BN2 = [];

% average Gradient and average Gradient squared holders
avgG.Dis = []; avgGS.Dis = []; avgG.Gen = []; avgGS.Gen = [];
%% Train
numIterations = floor(size(trainX,4)/settings.batch_size);
out = false; epoch = 0; global_iter = 0;
while ~out
    tic;
    shuffleid = randperm(size(trainX,4));
    trainXshuffle = trainX(:,:,:,shuffleid);
    trainYshuffle = trainY(shuffleid);
    fprintf('Epoch %d\n',epoch)
    for i=1:numIterations
        global_iter = global_iter+1;
        noise = gpdl(randn([settings.latentDim,...
            settings.batch_size]),'CB');
        idx = (i-1)*settings.batch_size+1:i*settings.batch_size;
        XBatch=gpdl(single(trainXshuffle(:,:,:,idx)),'SSCB');
        YBatch=gpdl(single(trainYshuffle(idx)),'B');

        [GradGen,GradDis,stGen,stDis] = ...
                dlfeval(@modelGradients,XBatch,YBatch,noise,...
                paramsGen,paramsDis,stGen,stDis,settings);

        % Update Discriminator network parameters
        [paramsDis,avgG.Dis,avgGS.Dis] = ...
            adamupdate(paramsDis, GradDis, ...
            avgG.Dis, avgGS.Dis, global_iter, ...
            settings.lrD, settings.beta1, settings.beta2);

        % Update Generator network parameters
        [paramsGen,avgG.Gen,avgGS.Gen] = ...
            adamupdate(paramsGen, GradGen, ...
            avgG.Gen, avgGS.Gen, global_iter, ...
            settings.lrG, settings.beta1, settings.beta2);

        if i==1 || rem(i,20)==0
            progressplot(paramsGen,stGen,settings);
%             if i==1 || (epoch>=0 && i==1)
%                 h = gcf;
%                 % Capture the plot as an image
%                 frame = getframe(h);
%                 im = frame2im(frame);
%                 [imind,cm] = rgb2ind(im,256);
%                 % Write to the GIF File
%                 if epoch == 0
%                   imwrite(imind,cm,'DCGANmnist.gif','gif', 'Loopcount',inf);
%                 else
%                   imwrite(imind,cm,'DCGANmnist.gif','gif','WriteMode','append');
%                 end
%             end
        end

    end

    elapsedTime = toc;
    disp("Epoch "+epoch+". Time taken for epoch = "+elapsedTime + "s")
    epoch = epoch+1;
    if epoch == settings.maxepochs
        out = true;
    end
end
%% Helper Functions
%% preprocess
function x = preprocess(x)
x = double(x)/255;
x = (x-.5)/.5;
x = reshape(x,28,28,1,[]);
end
%% extract data
function x = gatext(x)
x = gather(extractdata(x));
end
%% gpu dl array wrapper
function dlx = gpdl(x,labels)
dlx = gpuArray(dlarray(x,labels));
end
%% Weight initialization
function parameter = initializeGaussian(parameterSize,sigma)
if nargin < 2
    sigma = 0.05;
end
parameter = randn(parameterSize, 'single') .* sigma;
end
%% Generator
function [dly,st] = Generator(dlx,labels,params,st)
dly = embedding(dlx,labels,params);

% fully connected
dly = fullyconnect(dly,params.FCW1,params.FCb1);
dly = relu(dly);
dly = gpdl(reshape(dly,7,7,128,[]),'SSCB');
if isempty(st.BN1)
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,...
        params.BNo1,params.BNs1,'MeanDecay',0.8);
else
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,...
        params.BNs1,st.BN1.mu,st.BN1.sig,...
        'MeanDecay',.8);
end
% transposed convolution
dly = dltranspconv(dly,params.TCW1,params.TCb1,...
    'Stride',2,'Cropping','same');
dly = relu(dly);
if isempty(st.BN2)
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,...
        params.BNo2,params.BNs2,'MeanDecay',0.8);
else
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,...
        params.BNs2,st.BN2.mu,st.BN2.sig,...
        'MeanDecay',.8);
end
dly = dltranspconv(dly,params.TCW2,params.TCb2,...
    'Stride',2,'Cropping','same');
dly = relu(dly);
if isempty(st.BN3)
    [dly,st.BN3.mu,st.BN3.sig] = batchnorm(dly,...
        params.BNo3,params.BNs3,'MeanDecay',0.8);
else
    [dly,st.BN3.mu,st.BN3.sig] = batchnorm(dly,params.BNo3,...
        params.BNs3,st.BN3.mu,st.BN3.sig,...
        'MeanDecay',.8);
end

dly = dlconv(dly,params.CNW1,params.CNb1,...
            'Padding','same');
% tanh
dly = tanh(dly);
end
%% Discriminator
function [dly,st] = Discriminator(dlx,params,st)
% convolution
%1
dly = dlconv(dlx,params.CNW1,params.CNb1,...
            'Stride',2,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);
%2
dly = dlconv(dly,params.CNW2,params.CNb2,...
            'Stride',2,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);
if isempty(st.BN1)
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,...
        params.BNo1,params.BNs1,'MeanDecay',0.8);
else
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,...
        params.BNs1,st.BN1.mu,st.BN1.sig,...
        'MeanDecay',0.8);
end

%3
dly = dlconv(dly,params.CNW3,params.CNb3,...
            'Stride',2,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);
if isempty(st.BN2)
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,...
        params.BNo2,params.BNs2,'MeanDecay',0.8);
else
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,...
        params.BNs2,st.BN2.mu,st.BN2.sig,...
        'MeanDecay',0.8);
end

%4
dly = dlconv(dly,params.CNW4,params.CNb4,...
            'Stride',1,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);

% Fully connected
dly = gpdl(reshape(dly,4*4*128,[]),'CB');
dly = fullyconnect(dly,params.FCW1,params.FCb1);
% sigmoid
dly(1,:) = sigmoid(dly(1,:));
% softmax
dly(2:end,:) = softmax(dly(2:end,:));
end
%% modelGradients
function [GradGen,GradDis,stGen,stDis]=modelGradients(x,y,z,paramsGen,...
    paramsDis,stGen,stDis,settings)
y0 = randi([0,9],[settings.batch_size,1]);
ohey0 = onehotencoding(y0,settings.num_labels);
ohey = onehotencoding(y,settings.num_labels);

[fake_images,stGen] = Generator(z,y0,paramsGen,stGen);
d_output_real = Discriminator(x,paramsDis,stDis);
[d_output_fake,stDis] = Discriminator(fake_images,paramsDis,stDis);

% Loss due to true or not
% d_loss = -mean(.9*log(d_output_real+eps)+log(1-d_output_fake+eps));
d_loss_real = -.5*mean(log(d_output_real(1,:)+eps))-...
    .5*mean(log(sum(d_output_real(2:end,:).*ohey,1))+eps);
d_loss_fake = -.5*mean(log(1-d_output_fake(1,:)+eps))-...
    .5*mean(log(sum(d_output_fake(2:end,:).*ohey0,1))+eps);
d_loss = (d_loss_real+d_loss_fake)/2;
g_loss = -.5*mean(log(d_output_fake(1,:)+eps))-...
    .5*mean(log(sum(d_output_fake(2:end,:).*ohey0,1))+eps);

% For each network, calculate the gradients with respect to the loss.
GradGen = dlgradient(g_loss,paramsGen,'RetainData',true);
GradDis = dlgradient(d_loss,paramsDis);
end
%% progressplot
function progressplot(paramsGen,stGen,settings)
r = 2; c = 5;
labels = gpdl(single([0:9]'),'B');
noise = gpdl(randn([settings.latentDim,r*c]),'CB');
gen_imgs = Generator(noise,labels,paramsGen,stGen);
gen_imgs = reshape(gen_imgs,28,28,[]);

fig = gcf;
if ~isempty(fig.Children)
    delete(fig.Children)
end

I = imtile(gatext(gen_imgs));
I = rescale(I);
imagesc(I)
title("Generated Images")
colormap gray

drawnow;
end
%% dropout
function dly = dropout(dlx,p)
if nargin < 2
    p = .3;
end
[n,d] = rat(p);
mask = randi([1,d],size(dlx));
mask(mask<=n)=0;
mask(mask>n)=1;
dly = dlx.*mask;
end
%% embedding
function dly = embedding(dlx,labels,params)
% params EM W (latentDim,num_labels)
%               / (img_elements,num_labels)
maskW = params.EMW1(:,labels+1);
dly = dlx.*maskW;
end
%% one hot encoding
function ohe = onehotencoding(labels,numLabels)
numBatch = length(labels);
ohe = zeros(numLabels,numBatch);

for i = 1:numBatch
    ohe(labels(i)+1,i)=1;
end

end

author: zcemycl

repo: Matlab-GAN

save_path: github-repos/MATLAB/zcemycl-Matlab-GAN


text:
% EEGFILT -  (high|low|band)-pass filter data using two-way least-squares
%              FIR filtering. Optionally uses the window method instead of
%              least-squares. Multiple data channels and epochs supported.
%              Requires the MATLAB Signal Processing Toolbox.
% Usage:
%  >> [smoothdata] = eegfilt(data,srate,locutoff,hicutoff);
%  >> [smoothdata,filtwts] = eegfilt(data,srate,locutoff,hicutoff, ...
%                                    epochframes,filtorder,revfilt,firtype,causal);
% Inputs:
%   data        = (channels,frames*epochs) data to filter
%   srate       = data sampling rate (Hz)
%   locutoff    = low-edge frequency in pass band (Hz)  {0 -> lowpass}
%   hicutoff    = high-edge frequency in pass band (Hz) {0 -> highpass}
%   epochframes = frames per epoch (filter each epoch separately {def/0: data is 1 epoch}
%   filtorder   = length of the filter in points {default 3*fix(srate/locutoff)}
%   revfilt     = [0|1] reverse filter (i.e. bandpass filter to notch filter). {default 0}
%   firtype     = 'firls'|'fir1' {'firls'}
%   causal      = [0|1] use causal filter if set to 1 (default 0)
%
% Outputs:
%    smoothdata = smoothed data
%    filtwts    = filter coefficients [smoothdata <- filtfilt(filtwts,1,data)]
%
% See also: FIRLS, FILTFILT

% Author: Scott Makeig, Arnaud Delorme, Clemens Brunner SCCN/INC/UCSD, La Jolla, 1997

% Copyright (C) 4-22-97 from bandpass.m Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.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.

% 05-08-97 fixed frequency bound computation -sm
% 10-22-97 added MINFREQ tests -sm
% 12-05-00 added ERROR calls -sm
% 01-25-02 reformated help & license, added links -ad
% 03-20-12 added firtype option -cb

function [smoothdata,filtwts] = eegfilt(data,srate,locutoff,hicutoff,epochframes,filtorder,revfilt,firtype,causal)

if nargin<4
    fprintf('');
    help eegfilt
    return
end

%if ~exist('firls')
%   error('*** EEGFILT requires the signal processing toolbox. ***');
%end

[chans frames] = size(data);
if chans > 1 && frames == 1,
    help eegfilt
    error('input data should be a row vector.');
end
nyq            = srate*0.5;  % Nyquist frequency
%MINFREQ = 0.1/nyq;
MINFREQ = 0;

minfac         = 3;    % this many (lo)cutoff-freq cycles in filter
min_filtorder  = 15;   % minimum filter length
trans          = 0.15; % fractional width of transition zones

if locutoff>0 && hicutoff > 0 && locutoff > hicutoff,
    error('locutoff > hicutoff ???\n');
end
if locutoff < 0 || hicutoff < 0,
    error('locutoff | hicutoff < 0 ???\n');
end

if locutoff>nyq,
    error('Low cutoff frequency cannot be > srate/2');
end

if hicutoff>nyq
    error('High cutoff frequency cannot be > srate/2');
end

if nargin<6
    filtorder = 0;
end
if nargin<7
    revfilt = 0;
end
if nargin<8
    firtype = 'firls';
end
if nargin<9
    causal = 0;
end

if strcmp(firtype, 'firls')
    warning('Using firls to estimate filter coefficients. We recommend that you use fir1 instead, which yields larger attenuation. In future, fir1 will be used by default!');
end

if isempty(filtorder) || filtorder==0,
    if locutoff>0,
        filtorder = minfac*fix(srate/locutoff);
    elseif hicutoff>0,
        filtorder = minfac*fix(srate/hicutoff);
    end

    if filtorder < min_filtorder
        filtorder = min_filtorder;
    end
end

if nargin<5
    epochframes = 0;
end
if epochframes ==0,
    epochframes = frames;    % default
end
epochs = fix(frames/epochframes);
if epochs*epochframes ~= frames,
    error('epochframes does not divide frames.\n');
end

if filtorder*3 > epochframes,   % Matlab filtfilt() restriction
    fprintf('eegfilt(): filter order is %d. ',filtorder);
    error('epochframes must be at least 3 times the filtorder.');
end
if (1+trans)*hicutoff/nyq > 1
    error('high cutoff frequency too close to Nyquist frequency');
end

if locutoff > 0 && hicutoff > 0,    % bandpass filter
    if revfilt
        fprintf('eegfilt() - performing %d-point notch filtering.\n',filtorder);
    else
        fprintf('eegfilt() - performing %d-point bandpass filtering.\n',filtorder);
    end
    fprintf('            If a message, ''Matrix is close to singular or badly scaled,'' appears,\n');
    fprintf('            then Matlab has failed to design a good filter. As a workaround, \n');
    fprintf('            for band-pass filtering, first highpass the data, then lowpass it.\n');
    if strcmp(firtype, 'firls')
        f=[MINFREQ (1-trans)*locutoff/nyq locutoff/nyq hicutoff/nyq (1+trans)*hicutoff/nyq 1];
        fprintf('eegfilt() - low transition band width is %1.1g Hz; high trans. band width, %1.1g Hz.\n',(f(3)-f(2))*srate/2, (f(5)-f(4))*srate/2);
        m=[0       0                      1            1            0                      0];
    elseif strcmp(firtype, 'fir1')
        filtwts = fir1(filtorder, [locutoff, hicutoff]./(srate/2));
    end
elseif locutoff > 0                % highpass filter
    if locutoff/nyq < MINFREQ
        error(sprintf('eegfilt() - highpass cutoff freq must be > %g Hz\n\n',MINFREQ*nyq));
    end
    fprintf('eegfilt() - performing %d-point highpass filtering.\n',filtorder);
    if strcmp(firtype, 'firls')
        f=[MINFREQ locutoff*(1-trans)/nyq locutoff/nyq 1];
        fprintf('eegfilt() - highpass transition band width is %1.1g Hz.\n',(f(3)-f(2))*srate/2);
        m=[   0             0                   1      1];
    elseif strcmp(firtype, 'fir1')
        filtwts = fir1(filtorder, locutoff./(srate/2), 'high');
    end
elseif hicutoff > 0                %  lowpass filter
    if hicutoff/nyq < MINFREQ
        error(sprintf('eegfilt() - lowpass cutoff freq must be > %g Hz',MINFREQ*nyq));
    end
    fprintf('eegfilt() - performing %d-point lowpass filtering.\n',filtorder);
    if strcmp(firtype, 'firls')
        f=[MINFREQ hicutoff/nyq hicutoff*(1+trans)/nyq 1];
        fprintf('eegfilt() - lowpass transition band width is %1.1g Hz.\n',(f(3)-f(2))*srate/2);
        m=[     1           1              0                 0];
    elseif strcmp(firtype, 'fir1')
        filtwts = fir1(filtorder, hicutoff./(srate/2));
    end
else
    error('You must provide a non-0 low or high cut-off frequency');
end
if revfilt
    if strcmp(firtype, 'fir1')
        error('Cannot reverse filter using ''fir1'' option');
    else
        m = double(~m);
    end
end

if strcmp(firtype, 'firls')
    filtwts = firls(filtorder,f,m); % get FIR filter coefficients
end

smoothdata = zeros(chans,frames);
for e = 1:epochs                % filter each epoch, channel
    for c=1:chans
        try
            if causal
                 smoothdata(c,(e-1)*epochframes+1:e*epochframes) = filter(  filtwts,1,data(c,(e-1)*epochframes+1:e*epochframes));
            else smoothdata(c,(e-1)*epochframes+1:e*epochframes) = filtfilt(filtwts,1,data(c,(e-1)*epochframes+1:e*epochframes));
            end
        catch,
            if causal
                 smoothdata(c,(e-1)*epochframes+1:e*epochframes) = filter(  filtwts,1,double(data(c,(e-1)*epochframes+1:e*epochframes)));
            else smoothdata(c,(e-1)*epochframes+1:e*epochframes) = filtfilt(filtwts,1,double(data(c,(e-1)*epochframes+1:e*epochframes)));
            end
        end
        if epochs == 1
            if rem(c,20) ~= 0, fprintf('.');
            else               fprintf('%d',c);
            end
        end
    end
end
fprintf('\n');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab