text:
function mosaic = sift_mosaic(im1, im2)
% SIFT_MOSAIC Demonstrates matching two images using SIFT and RANSAC
%
%   SIFT_MOSAIC demonstrates matching two images based on SIFT
%   features and RANSAC and computing their mosaic.
%
%   SIFT_MOSAIC by itself runs the algorithm on two standard test
%   images. Use SIFT_MOSAIC(IM1,IM2) to compute the mosaic of two
%   custom images IM1 and IM2.

% AUTORIGHTS

if nargin == 0
  im1 = imread(fullfile(vl_root, 'data', 'river1.jpg')) ;
  im2 = imread(fullfile(vl_root, 'data', 'river2.jpg')) ;
end

% make single
im1 = im2single(im1) ;
im2 = im2single(im2) ;

% make grayscale
if size(im1,3) > 1, im1g = rgb2gray(im1) ; else im1g = im1 ; end
if size(im2,3) > 1, im2g = rgb2gray(im2) ; else im2g = im2 ; end

% --------------------------------------------------------------------
%                                                         SIFT matches
% --------------------------------------------------------------------

[f1,d1] = vl_sift(im1g) ;
[f2,d2] = vl_sift(im2g) ;

[matches, scores] = vl_ubcmatch(d1,d2) ;

numMatches = size(matches,2) ;

X1 = f1(1:2,matches(1,:)) ; X1(3,:) = 1 ;
X2 = f2(1:2,matches(2,:)) ; X2(3,:) = 1 ;

% --------------------------------------------------------------------
%                                         RANSAC with homography model
% --------------------------------------------------------------------

clear H score ok ;
for t = 1:100
  % estimate homograpyh
  subset = vl_colsubset(1:numMatches, 4) ;
  A = [] ;
  for i = subset
    A = cat(1, A, kron(X1(:,i)', vl_hat(X2(:,i)))) ;
  end
  [U,S,V] = svd(A) ;
  H{t} = reshape(V(:,9),3,3) ;

  % score homography
  X2_ = H{t} * X1 ;
  du = X2_(1,:)./X2_(3,:) - X2(1,:)./X2(3,:) ;
  dv = X2_(2,:)./X2_(3,:) - X2(2,:)./X2(3,:) ;
  ok{t} = (du.*du + dv.*dv) < 6*6 ;
  score(t) = sum(ok{t}) ;
end

[score, best] = max(score) ;
H = H{best} ;
ok = ok{best} ;

% --------------------------------------------------------------------
%                                                  Optional refinement
% --------------------------------------------------------------------

function err = residual(H)
 u = H(1) * X1(1,ok) + H(4) * X1(2,ok) + H(7) ;
 v = H(2) * X1(1,ok) + H(5) * X1(2,ok) + H(8) ;
 d = H(3) * X1(1,ok) + H(6) * X1(2,ok) + 1 ;
 du = X2(1,ok) - u ./ d ;
 dv = X2(2,ok) - v ./ d ;
 err = sum(du.*du + dv.*dv) ;
end

if exist('fminsearch') == 2
  H = H / H(3,3) ;
  opts = optimset('Display', 'none', 'TolFun', 1e-8, 'TolX', 1e-8) ;
  H(1:8) = fminsearch(@residual, H(1:8)', opts) ;
else
  warning('Refinement disabled as fminsearch was not found.') ;
end

% --------------------------------------------------------------------
%                                                         Show matches
% --------------------------------------------------------------------

dh1 = max(size(im2,1)-size(im1,1),0) ;
dh2 = max(size(im1,1)-size(im2,1),0) ;

figure(1) ; clf ;
subplot(2,1,1) ;
imagesc([padarray(im1,dh1,'post') padarray(im2,dh2,'post')]) ;
o = size(im1,2) ;
line([f1(1,matches(1,:));f2(1,matches(2,:))+o], ...
     [f1(2,matches(1,:));f2(2,matches(2,:))]) ;
title(sprintf('%d tentative matches', numMatches)) ;
axis image off ;

subplot(2,1,2) ;
imagesc([padarray(im1,dh1,'post') padarray(im2,dh2,'post')]) ;
o = size(im1,2) ;
line([f1(1,matches(1,ok));f2(1,matches(2,ok))+o], ...
     [f1(2,matches(1,ok));f2(2,matches(2,ok))]) ;
title(sprintf('%d (%.2f%%) inliner matches out of %d', ...
              sum(ok), ...
              100*sum(ok)/numMatches, ...
              numMatches)) ;
axis image off ;

drawnow ;

% --------------------------------------------------------------------
%                                                               Mosaic
% --------------------------------------------------------------------

box2 = [1  size(im2,2) size(im2,2)  1 ;
        1  1           size(im2,1)  size(im2,1) ;
        1  1           1            1 ] ;
box2_ = inv(H) * box2 ;
box2_(1,:) = box2_(1,:) ./ box2_(3,:) ;
box2_(2,:) = box2_(2,:) ./ box2_(3,:) ;
ur = min([1 box2_(1,:)]):max([size(im1,2) box2_(1,:)]) ;
vr = min([1 box2_(2,:)]):max([size(im1,1) box2_(2,:)]) ;

[u,v] = meshgrid(ur,vr) ;
im1_ = vl_imwbackward(im2double(im1),u,v) ;

z_ = H(3,1) * u + H(3,2) * v + H(3,3) ;
u_ = (H(1,1) * u + H(1,2) * v + H(1,3)) ./ z_ ;
v_ = (H(2,1) * u + H(2,2) * v + H(2,3)) ./ z_ ;
im2_ = vl_imwbackward(im2double(im2),u_,v_) ;

mass = ~isnan(im1_) + ~isnan(im2_) ;
im1_(isnan(im1_)) = 0 ;
im2_(isnan(im2_)) = 0 ;
mosaic = (im1_ + im2_) ./ mass ;

figure(2) ; clf ;
imagesc(mosaic) ; axis image off ;
title('Mosaic') ;

if nargout == 0, clear mosaic ; end

end

author: Cloud-CV

repo: object-proposals


text:
function v = ompver(history)
%OMPVER OMP toolbox version information.
%  OMPVER displays the current OMP toolbox version information.
%
%  OMPVER('history') also displays history information about the previous
%  versions of the OMP toolbox and their change logs.
%
%  V = OMPVER returns the version number of the current OMP toolbox, and
%  does not display any information.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  October 2009

ver = 10;

if (nargout>0)
  if (nargin>0)
    error('Invalid number of parameters.');
  end
  v = ver;

else

  if (nargin==0 || (nargin==1 && strcmpi(history,'history')))

    disp(' ');
    disp('---------------------------------');
    printf('    OMP Toolbox version %d        ',ver);
    disp('---------------------------------');
    disp(' ');

  else
    error('Unknown parameters.');
  end

  if (nargin>0)

    %% changes only displayed since first public release %%

%     disp(' ');
%     disp('OMP Toolbox version 1, 28.10.07');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Initial release.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('OMP Toolbox version 2, 30.10.07');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. OMP and OMP2 can now operate without specifying G, computing it on the fly.');
%     disp('  2. Improved documentation.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. OMP2 uses much less memory, by allocating space for at most N coefficients');
%     disp('     for each signal rather than M in v1.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('OMP Toolbox version 3, 16.12.07');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. Support for sparse output, and an optional parameter for selecting');
%     disp('     between sparse and full output.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. When selected atoms were dependent, OMP would cause a NaN result.');
%     disp('     Added a stopping criterion to handle this case.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. Dynamic memory allocations now done using mxMalloc/mxCalloc.');
%     disp('  2. Integer types converted to mwSize/mwIndex to allow compilation');
%     disp('     under 64-bit systems using the -largeArrayDims parameter.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('OMP Toolbox version 4, 03.01.08');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. Sparse Gamma''s are now represented correctly, with ascending row numbers.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. Added MAKE function for automated compilation of all OMP functions.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('OMP Toolbox version 5, 18.01.08');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. Significantly improved performance by using BLAS/LAPACK functions.');
%     disp('  2. More meaningful profiling reports.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. Changes to code structure: extracted core omp code to separate');
%     disp('     callable functions for code re-use.');
%     disp('  2. Improved MAKE file identifies compilation platform automatically.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('OMP Toolbox version 6, 31.01.08');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. Implemented efficient incremental update of the residual norm in OMP2,');
%     disp('     eliminating the need for explicit computation of the residual itself.');
%     disp('  2. Added ''maxatoms'' parameter to OMP2.');
%     disp('  3. Updated call syntax to support the changes in (1) and (2).');
%     disp('  4. Slightly improved documentation.');
%     disp('  5. Added the OMPVER function.');
%     disp('  6. Added the Contents.m file.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. Better memory management in OMP2. Much less memory is consumed by using');
%     disp('     dynamic reallocation of temporary helper matrices which were originaly');
%     disp('     preallocated to a size proportional to G, even when G was not specified.');
%     disp('  2. OMP2: atom selection loop limited to N atoms (loop could continue beyond');
%     disp('     this limit if the error target was not achieved, causing a runtime error.');
%     disp('     This could occur when the error limit is very close to 0).');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('OMP Toolbox version 7, 15.04.08');
%     disp('--------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. OMP2: Improved incremental computation of the residual norm requires less');
%     disp('     memory and is more efficient.');
%     disp('  2. Eliminated separate profiling OMP functions, now unified into the OMP');
%     disp('     and OMP2 functions.');
%     disp('  3. New ''messages'' option for displaying coding progress and remaining time.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. New code structure: unified all omp and omp2 variants into ompmex.c and')
%     disp('     omp2mex.c. Parameter parsing is now done in the Matlab omp.m and omp2.m');
%     disp('     functions, allowing simpler and more advanced parsing.');
%     disp('  2. Much simpler MAKE function, and new create_blasfuncs.m function for');
%     disp('     managing the BLAS function names.');
%     disp(' ');
%     disp(' ');
%     disp(' ');

    disp(' ');
    disp(' ');
    disp('OMP Toolbox version 7, 15.04.08');
    disp('--------------------------------');
    disp(' ');
    disp('First public release.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('OMP Toolbox version 8, 29.05.08');
    disp('--------------------------------');
    disp(' ');
    disp('Bug fixes:');
    disp(' ');
    disp('  1. Added OMP stopping criterion: when maximal inner product between the');
    disp('     residual and the dictionary is too small.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('OMP Toolbox version 9, 21.05.09');
    disp('--------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. Source code now available for compilation on any Matlab platform.');
    disp('  2. Dictionary normalization is now verified before performing OMP.');
    disp('  3. Added OMPDEMO and OMPSPEEDTEST.');
    disp('  4. Improved profiling information.');
    disp(' ');
    disp('Internal changes:');
    disp(' ');
    disp('  1. Completely restructured code.');
    disp('  2. When G is not specified, OMP-Cholesky is used rather than computing G');
    disp('     and using Batch-OMP (faster).');
    disp('  3. Estimated remaining time now computed with sub-second accuracy.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('OMP Toolbox version 10, 18.10.09');
    disp('---------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. When the sparse-coding threshold is such that no atoms need to be');
    disp('     selected, OMP now terminates without computing D''*x, leading to');
    disp('     a significant speed increase.');
    disp('  2. New ''checkdict'' parameter in OMP and OMP2 allows disabling dictionary');
    disp('     normalization verification for accelerated performance.');
    disp(' ');
    disp('Internal changes:');
    disp(' ');
    disp('  1. Restructured and updated C libraries for compatibility with the new');
    disp('     Sparse OMP Toolbox.');
    disp('  2. Slightly modified parameters in OMPSPEEDTEST for consistency with');
    disp('     OMPSPEEDCOMPARE in the Sparse OMP Toolbox.');
    disp(' ');
    disp(' ');

  end
end

author: thfylsty


text:
function [Gs, op, nodes] = mk_nbrs_of_digraph2(G0)
% MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal
% [Gs, op, nodes] = mk_nbrs_of_digraph(G0)
% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor.
% nodes(i,1:2) are the head and tail of the operated-on arc.

[I,J] = find(G0);
G0bar = setdiag(~G0, 0); % exclude self loops in graph complement
[Ibar,Jbar] = find(G0bar);
nnbrs = 2*length(I) + length(Ibar);
Gs = cell(1, nnbrs);
op = cell(1, nnbrs);
nodes = zeros(nnbrs, 2);

nbr = 1;
% all single edge deletions
for e=1:length(I)
  i = I(e); j = J(e);
  G = G0;
  G(i,j) = 0;
  Gs{nbr} = G;
  op{nbr} = 'del';
  nodes(nbr, :) = [i j];
  nbr = nbr + 1;
end

% all single edge reversals
for e=1:length(I)
  i = I(e); j = J(e);
  G = G0;
  G(i,j) = 0;
  G(j,i) = 1;
  Gs{nbr} = G;
  op{nbr} = 'rev';
  nodes(nbr, :) = [i j];
  nbr = nbr + 1;
end

[I,J] = find(~G0);
% all single edge additions
for e=1:length(I)
  i = I(e); j = J(e);
  G = G0;
  if i ~= j % don't add self loops
    G(i,j) = 1;
    Gs{nbr} = G;
    op{nbr} = 'add';
    nodes(nbr, :) = [i j];
    nbr = nbr + 1;
  end
end

assert(nnbrs == nbr-1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef MMMOP3 < PROBLEM
% <multi/many> <real> <multimodal>
% Multi-modal multi-objective optimization problem
% kA --- 1 --- Number of decision variables in XA
% c  --- 2 --- Parameter c
% d  --- 2 --- Parameter d

%------------------------------- Reference --------------------------------
% Y. Liu, G. G. Yen, and D. Gong, A multi-modal multi-objective
% evolutionary algorithm using two-archive and recombination strategies,
% IEEE Transactions on Evolutionary Computation, 2019, 23(4): 660-674.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(Access = private)
        kA;     % Number of decision variables in XA
        c;      % Parameter c
        d;      % Parameter d
        POS;    % Pareto optimal set for IGDX calculation
    end
	methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 2; end
            if isempty(obj.D); obj.D = obj.M + 1; end
            [obj.kA,obj.c,obj.d] = obj.ParameterSet(1,2,2);
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            y      = PopDec(:,1:obj.M-1).*obj.d-floor(PopDec(:,1:obj.M-1).*obj.d);
            g      = obj.kA+sum(cos(2.*pi.*obj.c.*PopDec(:,obj.M:obj.M-1+obj.kA)),2)+sum((PopDec(:,obj.M+obj.kA:end)-0.5).^2,2);
            PopObj = repmat(1+g,1,obj.M).*fliplr(cumprod([ones(size(g,1),1),cos(y*pi/2)],2)).*[ones(size(g,1),1),sin(y(:,obj.M-1:-1:1)*pi/2)];
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            XA      = Grid(0.5/obj.c:1/obj.c:(1-0.5/obj.c),obj.kA);
            X       = UniformPoint(N/size(XA,1),obj.M-1,'grid');
            obj.POS = [repmat(X,size(XA,1),1),XA(repmat(1:end,size(X,1),1),:),zeros(size(X,1)*size(XA,1),obj.D-size(X,2)-size(XA,2))+0.5];
            % Generate points on Pareto front
            R = UniformPoint(N,obj.M);
            R = R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                a = linspace(0,pi/2,100)';
                R = [sin(a),cos(a)];
            elseif obj.M == 3
                a = linspace(0,pi/2,10)';
                R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
            else
                R = [];
            end
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            XA = Grid(0.5/obj.c:1/obj.c:(1-0.5/obj.c),obj.kA);
            X  = Grid(0.5/obj.d:1/obj.d:(1-0.5/obj.d),obj.M-1);
            R  = [repmat(X,size(XA,1),1),XA(repmat(1:end,size(X,1),1),:)];
            [~,Label]  = min(pdist2(PopDec(:,1:obj.M+obj.kA-1),R),[],2);
            tempStream = RandStream('mlfg6331_64','Seed',2);
            if obj.M == 2
                for i = 1 : size(R,1)
                    color = rand(tempStream,1,3);
                    Draw(Population(Label==i).objs+(i-1)*0.05,'o','MarkerSize',6,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2',[]});
                    Draw(obj.PF+(i-1)*0.05,'-','LineWidth',1,'Color',color);
                end
            elseif obj.M == 3
                for i = 1 : size(R,1)
                    color = rand(tempStream,1,3);
                    ax = Draw(Population(Label==i).objs+(i-1)*0.05,'o','MarkerSize',8,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2','\it f\rm_3'});
                    surf(ax,obj.PF{1}+(i-1)*0.05,obj.PF{2}+(i-1)*0.05,obj.PF{3}+(i-1)*0.05,'EdgeColor',color,'FaceColor','none');
                end
            else
                for i = 1 : size(R,1)
                    Draw(Population(Label==i).objs,'-','Color',rand(tempStream,1,3),'LineWidth',2);
                end
            end
        end
    end
end

function W = Grid(gap,M)
    if M < 1
        W = zeros(1,0);
    else
        eval(sprintf('[%s]=ndgrid(gap);',sprintf('c%d,',1:M)))
        eval(sprintf('W=[%s];',sprintf('c%d(:),',1:M)))
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [tt_chg, waveform_chg] = restoreAdditionalShapeSamples(tt,waveform,first,last,gradRasterTime,iBlock)
% restore shape: if we had a
% trapezoid converted to shape we have to find
% the "corners" and we can eliminate internal
% samples on the straight segments
% but first we have to restore samples on the
% edges of the gradient raster intervals
% for that we need the first sample
    max_abs=max(abs(waveform));
    odd_step1=[first 2*waveform'];
    odd_step2=odd_step1.*(mod(1:length(odd_step1),2)*2-1);
    waveform_odd_rest=(cumsum(odd_step2).*(mod(1:length(odd_step2),2)*2-1))';
    waveform_odd_interp=[first; 0.5*(waveform(1:end-1)+waveform(2:end)); last];
    if abs(waveform_odd_rest(end)-last)>2e-5*max_abs % what's the reasonable threshold?
        blInfo='';
        if exist('iBlock')
            blInfo=['[block ' num2str(iBlock) '] '];
        end
        warning('mr:restoreShape',[blInfo 'Last restored point ' ...
                 'differs too much from the recorded last, skipping the shape restoration step; ' ...
                 'deviation: ' num2str(abs(waveform_odd_rest(end)-last)) 'Hz/m (' num2str(abs(waveform_odd_rest(end)-last)/max_abs*100) '%%); ' ...
                 'No reason to panic, this is typical for spirals. ' ...
                 'To disable this warning (assuming you know what you are doing) add \n' ...
                 '`warning(''OFF'', ''mr:restoreShape'')´ to your sequence-generating script.']);
        tt_chg=[0 tt' tt(end)+gradRasterTime/2];
        waveform_chg = [first waveform' last];
        return;
    end
    %figure; plot([0,10e-6+grad.t'],waveform_odd_rest-waveform_odd_interp);
    waveform_odd_mask=abs(waveform_odd_rest-waveform_odd_interp)<=eps+2e-5*max_abs; % threshold ???
    waveform_odd=waveform_odd_interp.*waveform_odd_mask+waveform_odd_rest.*(1-waveform_odd_mask);

    % combine odd & even
    comb=[ 0 waveform' ; waveform_odd' ];
    waveform_os=comb(2:end)';

    tt_odd=(0:(length(waveform_odd_rest)-1))*gradRasterTime;
    tt_os=(0:(length(waveform_os)-1))*gradRasterTime*0.5;

    waveform_even_reint=0.5*(waveform_odd_rest(1:end-1)+waveform_odd_rest(2:end));

    maskChanges = abs([1; diff(waveform_os,2); 1])>1e-8;   % TRUE if values change
    waveform_chg = waveform_os(maskChanges)';                     % Elements without repetitions
    tt_chg=tt_os(maskChanges);
    %figure;plot(grad.tt,grad.waveform);hold on; plot(tt_chg,waveform_chg); plot(tt_chg,waveform_chg,'o');
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function X = QLiftRec2MinMin(C, S)
%-----------------------------------------------------------------------------
% QLiftRec2MinMin
% Multilevel 2-D reconstruction by inverting the lifting scheme and using
% quincunx grids.
%
% The MinMin scheme has been proposed by Heijmans and Goutsias, see e.g.
%    H.J.A.M. Heijmans, J. Goutsias,
%    Multiresolution signal decomposition schemes.
%    Part 2: morphological wavelets.
%    CWI Report PNA-R9905, Amsterdam, 1999.
%    http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04625D.pdf
%
% Calls for:
%            retrieveQ1001, retrieveR,
%            getcolor01, getcolor10, getcolor00, getcolor11,
%            putcolor01, putcolor10, putcolor00, putcolor11.
% See also: QLiftDec2MaxMin, QLiftRec2
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: February 3, 2003.
% (c) 1999-2003 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
% Note: argument list of this function might be extended with an argument that
%       points to another level than 1, nargin could be checked for this.
% Firstly, check input data
%
if isempty(C)
  error(' QLiftRec2MinMin - empty decomposition ');
else
  if isempty(S)
    error(' QLiftRec2MinMin - empty bookkeeping ');
  end
end
%
N = numoflevs(S);
if mod(N, 2) == 1
  error(' QLiftRec2MinMin - only an even number of levels is accepted ');
end
%
%Secondly, start reconstruction
%
for lev=N:-2:1
%
%  The Inverse Scheme proceeds from rectangular grid to quincunx grid.
   if lev >= N
     APPROX00 = retrieveR(lev, 'a', C, S);     % "even slots"
   end
   DETAIL11 = retrieveR(lev, 'd', C, S);       % "odd slots"
   minO = min( min(min(APPROX00)), min(min(DETAIL11)) );
   maxO = max( max(max(APPROX00)), max(max(DETAIL11)) );
%  cmin = minO-(maxO-minO);
   cmax = maxO+(maxO-minO);
%
%  Stage: undo update
   sizeA00 = size(APPROX00);
   Q0011A00 = APPROX00 - ...
              min(zeros(sizeA00), synA00Qmin(DETAIL11, sizeA00, cmax));
   clear APPROX00 sizeA00;
%
%  Stage: undo predict
   Q0011A11 = DETAIL11 + synA11Qmin(Q0011A00, size(DETAIL11), cmax);
   clear DETAIL11;
%
%  Merge
%  The union Q0011 of Q0011A00 & Q0011A11 now contains the approximation on
%  the next scale (with index lev-1).
%
%  The Inverse Scheme proceeds from quincunx grid to rectangular grid.
%
%  The "even slots" are in the colours 00 and 11,
%  the "odd slots"  are in the colours 10 and 01.
%
%  Stage: undo update
   [DETAIL10, DETAIL01] = retrieveQ1001(lev-1, 'd', C, S);
   A00 = Q0011A00 - ...
         min(zeros(size(Q0011A00)), synA00min(DETAIL10, DETAIL01, cmax));
   clear Q0011A00;
   A11 = Q0011A11 - ...
         min(zeros(size(Q0011A11)), synA11min(DETAIL10, DETAIL01, cmax));
   clear Q0011A11;
%
%  Stage: undo predict
   A10 = DETAIL10 + synA10min(A11, A00, cmax);
   A01 = DETAIL01 + synA01min(A11, A00, cmax);
   sizeR = size(DETAIL10) + size(DETAIL01);
   clear DETAIL10 DETAIL01;
%
%  Merge
   APPROX00 = putcolor00(A00, sizeR) + putcolor11(A11, sizeR) + ...
              putcolor10(A10, sizeR) + putcolor01(A01, sizeR);
%  APPROX00 is used in the next iteration of the "for"-loop!
   clear A00 A11 A10 A01;
%
%  Warning: not yet verified whether sizeQ and sizeR will be consistent for all
%           possible griddimensions of original image.
%
   if (lev-1) == 1
     X = APPROX00; clear APPROX00;
   end
end
%-----------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function f1 = p11_f1 ( x )

%*****************************************************************************80
%
%% P11_F1 evaluates the first derivative for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  if ( 3.0 * x - 1.0 < 0.0 )
    f1 = - 3.0;
  else
    f1 = + 3.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [EEG, params, computationTimes] = prepPipeline(EEG, params)
% Run the PREP pipeline for standardized EEG data preparation
%
% Input parameters:
%  EEG                       An EEGLAB structure with the data and chanlocs
%  params                    A structure with at least the following:
%
%     name                   A string with a name identifying dataset
%                            [default: EEG.setname]
%     referenceChannels      A vector of channels to be used for
%                            rereferencing [default: all channels]
%     evaluationChannels     A vector of channels to be used for
%                            EEG interpolation [default: all channels]
%     rereferencedChannels   A vector of channels to be
%                            line-noise removed, and referenced
%                            [default: all channels]
%     lineFrequencies        A list of line frequencies to be removed
%                            [default: 60, 120, 180, 240]
%
% Prep allows many other parameters to be overridden, but is meant to
% be used in a fully automated fashion. For a full listing of the
% defaults for each step in the pipeline, execute:
%     showPrepDefaults(EEG).
%
% Output parameters:
%   EEG                      An EEGLAB structure with the data processed
%                            and status written in EEG.etc.noiseDetection
%   computationTimes         Time in seconds for each stage
%
% Additional setup:
%    EEGLAB should be in the path.
%    The EEG-Clean-Tools/PrepPipeline directory and its subdirectories
%    should be in the path.
%
% Author:  Kay Robbins, UTSA, March 2015
%
% Full documentation is available in a user manual distributed with this
% source.

%% Setup the output structures and set the input parameters
    computationTimes= struct( 'detrend', 0, 'lineNoise', 0, 'reference', 0);
    errorMessages = struct('status', 'good', 'boundary', 0, ...
                   'detrend', 0, 'lineNoise', 0, 'reference', 0, ...
                   'postProcess', 0);
    [backupOptionsFile, currentOptionsFile, warningsState] = setupForEEGLAB();
    finishup = onCleanup(@() cleanup(backupOptionsFile, currentOptionsFile, ...
        warningsState));
    if isfield(EEG.etc, 'noiseDetection')
        warning('EEG.etc.noiseDetection already exists and will be cleared\n')
    end
    if ~exist('params', 'var')
        params = struct();
    end
    if ~isfield(params, 'name')
        params.name = ['EEG' EEG.filename];
    end
    EEG.etc.noiseDetection = ...
        struct('name', params.name, 'version', getPrepVersion, ...
               'originalChannelLabels', [], ...
               'errors', [], 'boundary', [], 'detrend', [], ...
               'lineNoise', [], 'reference', [], 'postProcess', [], ...
               'interpolatedChannelNumbers', [], 'removedChannelNumbers', [], ...
               'stillNoisyChannelNumbers', [], 'fullReferenceInfo', false);
    EEG.data = double(EEG.data);   % Don't monkey around -- get into double
    EEG.etc.noiseDetection.originalChannelLabels = {EEG.chanlocs.labels};

%% Check for the general defaults
    try
        defaults = getPrepDefaults(EEG, 'general');
        [params, errors] = checkPrepDefaults(params, params, defaults);
        if ~isempty(errors)
            error('prepPipeline:GeneralDefaultError', ['|' sprintf('%s|', errors{:})]);
        end
    catch mex
        warning('[%s] %s: Prep could not begin processing the EEG', ...
           mex.identifier, getReport(mex, 'extended', 'hyperlinks', 'on'));
        return;
    end

%% Check for boundary events
    fprintf('Checking for boundary events\n');
    try
        defaults = getPrepDefaults(EEG, 'boundary');
        [boundaryOut, errors] = checkPrepDefaults(params, struct(), defaults);
        if ~isempty(errors)
            error('boundary:BadParameters', ['|' sprintf('%s|', errors{:})]);
        end
        EEG.etc.noiseDetection.boundary = boundaryOut;
        if ~boundaryOut.ignoreBoundaryEvents && ...
                isfield(EEG, 'event') && ~isempty(EEG.event)
            eTypes = find(strcmpi({EEG.event.type}, 'boundary'));
            if ~isempty(eTypes)
                error('boundary:UnremovedBoundary', ['Dataset ' params.name  ...
                    ' has boundary events: [' getListString(eTypes) ...
                    '] which are treated as discontinuities unless ' ...
                    'set to ignore. Prep cannot continue']);
            end
        end
    catch mex
        errorMessages.boundary = ['prepPipeline bad boundary events: ' ...
             getReport(mex, 'basic', 'hyperlinks', 'off')];
        errorMessages.status = 'unprocessed';
        EEG.etc.noiseDetection.errors = errorMessages;
        if strcmpi(params.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
               getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
    end

%% Part II:  HP the signal for detecting bad channels
    fprintf('Preliminary detrend to compute reference\n');
    try
        tic
        [EEGNew, detrend] = removeTrend(EEG, params);
        EEG.etc.noiseDetection.detrend = detrend;
        % Make sure detrend defaults are available for referencing
        defaults = getPrepDefaults(EEG, 'detrend');
        params = checkPrepDefaults(detrend, params, defaults);
        computationTimes.detrend = toc;
    catch mex
        errorMessages.detrend = ['prepPipeline failed removeTrend: ' ...
             getReport(mex, 'basic', 'hyperlinks', 'off')];
        errorMessages.status = 'unprocessed';
        EEG.etc.noiseDetection.errors = errorMessages;
        if strcmpi(params.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
                getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
    end

%% Part III: Remove line noise
    fprintf('Line noise removal\n');
    try
        tic
        [EEGClean, lineNoise] = removeLineNoise(EEGNew, params);
        EEG.etc.noiseDetection.lineNoise = lineNoise;
        lineChannels = lineNoise.lineNoiseChannels;
        EEG.data(lineChannels, :) = EEG.data(lineChannels, :) ...
             - EEGNew.data(lineChannels, :) + EEGClean.data(lineChannels, :);
        clear EEGNew;
        computationTimes.lineNoise = toc;
    catch mex
        errorMessages.lineNoise = ['prepPipeline failed removeLineNoise: ' ...
            getReport(mex, 'basic', 'hyperlinks', 'off')];
        errorMessages.status = 'unprocessed';
        EEG.etc.noiseDetection.errors = errorMessages;
        if strcmpi(params.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
                getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
    end

%% Part IV: Find reference
    fprintf('Find reference\n');
    try
        tic
        [EEG, referenceOut] = performReference(EEG, params);
        EEG.etc.noiseDetection.reference = referenceOut;
        EEG.etc.noiseDetection.fullReferenceInfo = true;
        EEG.etc.noiseDetection.interpolatedChannelNumbers = ...
            referenceOut.interpolatedChannels.all;
        EEG.etc.noiseDetection.stillNoisyChannelNumbers = ...
            referenceOut.noisyStatistics.noisyChannels.all;
        computationTimes.reference = toc;
    catch mex
        errorMessages.reference = ['prepPipeline failed to perform reference: ' ...
            getReport(mex, 'basic', 'hyperlinks', 'off')];
        errorMessages.status = 'unprocessed';
        EEG.etc.noiseDetection.errors = errorMessages;
        if strcmpi(params.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
                getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
    end

    %% Report that there were no errors in the prep running
    EEG.etc.noiseDetection.errors = errorMessages;
end

%% Cleanup callback
function cleanup(backupFile, currentFile, warningsState)
% Restore EEGLAB options file and warning settings
   restoreEEGOptions(backupFile, currentFile);
   warning(warningsState);
end % cleanup

author: VisLab

repo: EEG-Clean-Tools


text:
function [rsmap] = computeTSResStd(view,scanNum,sliceNum,freq)
%
% function [rsmap] = computeTSResStd(view,scanNum,sliceNum,ncycles)
%
% Computes residual tseries std map removing freq and its harmonics
%
%  typical freq=6 cycles/scan
%
% 05/05/99 rmk
% 9/27/99 dbr  Removed hard-coded "72" length.
% djh, 11/00  Calls new version of percentTSeries that has the option of dividing
%             by spatialGradient (estimate of intensity inhomogeneity) instead
%             of dividing by mean at each pixel. Changed input arguments
%             correspondingly.

% Call percentTSeries to load it, and remove dc and trend
view = percentTSeries(view, scanNum, sliceNum);
ptSeries = view.tSeries;
tl=size(ptSeries,1);

% remove only 1f, 2f, 3f, & 5f
freqs=[freq,2*freq,3*freq,5*freq];

% remove all the harmonics
%freqs=[freq:freq:floor(tl/2)-1];

% remove everything but freq+1 and freq-1
%freqs=[1:freq-2,freq,freq+2:floor(tl/2)-1];

nyquistFreq = floor(tl/2)+1;
plusHarmonics=freqs+1;
minusHarmonics=tl-freqs+1;
allFreqs=[plusHarmonics,minusHarmonics];

% Compute fourier transform, zero out freqs, and inverse transform
ft=fft(ptSeries);
ft(allFreqs,:)=0;

% Compute std
rsmap = std(real(ifft(ft)));

% Compensate for the number of zero'd freqs so that this is an estimate of
% the std of the noise in the time series as if there were no signal modulation.
rsmap = rsmap * sqrt((tl-1)/(tl-length(allFreqs)-1));

% This last step converts it so that it is an estimate of the reliability across
% repeated scans. The reliability of the amplitude estimate depends of course
% on the number of time samples in each scan. For a sinusoid plus noise:
%     std(amps across repeated scans) = std(noise in time series) * sqrt(2/N)
% See simulation in debugging code below.
rsmap = rsmap*sqrt(2/tl);

return

% Debug/test

% Generate simulated data set (1 cyc/scan plus noise)
tl=20; freq = 1;
%tl=30; freq = 2;
t = linspace(0,2*pi,tl+1);
t = t(1:tl)';
sine = sin(freq*t);
ptSeries = repmat(sine,[1,1e4]);
ptSeries = ptSeries + randn(size(ptSeries));

% Compute resStd (code copied from above);
freqs = [freq,2*freq,3*freq,5*freq];
nyquistFreq = floor(tl/2)+1;
plusHarmonics = freqs+1;
minusHarmonics = tl-freqs+1;
allFreqs = [plusHarmonics,minusHarmonics];
ft = fft(ptSeries);
ft(allFreqs,:) = 0;
rsmap = std(real(ifft(ft)));
rsmap = rsmap * sqrt((tl-1)/(tl-length(allFreqs)-1));
rsmap = rsmap * sqrt(2/tl);

% Compute std of amp estimates
ft = fft(ptSeries);
amps = 2*abs(ft(freq+1,:))/tl;

% These two numbers should be the same
mean(rsmap)
std(amps)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [X, meta] = partial_FDCT(x, is_real, finest, nbscales, nbangles_coarse)
%Apply a forward fast discrete curvelet transformation to x.
%The input image is of the size MxN
%The resulting pyramidial structure of the curvelet coefficients is sorted by size
%and truncated to the size MN. This way, the transformed image X and the original image x
%have same size. The original pyramidial structure will be reconstructed while the iFDCT is performed, filling
%neglected coefficients with 0.

[M, N] 	= size(x);
C 	= fdct_wrapping(x, is_real, finest, nbscales, nbangles_coarse);

stacked = [];
for s=1:length(C)
  for w=1:length(C{s})
    stacked = [stacked; C{s}{w}(:)];
  end
end

[sorted, indicees] 	= sort(abs(stacked), 'descend');
wholeLength 		= M*N;
X                   = reshape(stacked(indicees(1:wholeLength)),M,N);
stackedPos          = reshape(indicees(1:wholeLength),M,N);

meta 	= struct('wholeTransform', {C}, 'wholeLength', wholeLength, 'stackedPos', stackedPos); % meta data needed while computing the iFDCT
end

author: thomaskuestner


text:
function varargout = hikmeanshist(varargin)
% VL_HIKMEANSHIST  Compute histogram of quantized data
%  H = VL_HIKMEANSHIST(TREE,PATH) computes the histogram of the HIKM tree
%  nodes activated by the root-to-leaf paths PATH. PATH is usually
%  obtained by quantizing data by means of VL_HIKMEANSPUSH().
%
%  The histogram H has one bin for each node of the HIKM tree TREE.
%  The tree has K = TREE.K nodes and depth D = TREE.DEPTH.  Therefore
%  there are M = (K^(D+1) - 1) / (K - 1) nodes in the tree (not
%  counting the root which carries no information). Nodes are stacked
%  into a vector of bins in breadth first order.
%
%  Example::
%    The following relations illustrate the structure of PATH:
%      H(1)   = # of paths such that PATH(1,:) = 1
%      H(K)   = # of paths such that PATH(1,:) = K
%      H(K+1) = # of paths such that PATH(1:2,:) = [1 ; 1]
%      H(K+K) = # of paths such that PATH(1:2,:) = [1 ; K]
%
%  See also: VL_HIKMEANS(), VL_HIKMEANSPUSH(), VL_HELP().
[varargout{1:nargout}] = vl_hikmeanshist(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% 	Problem 8- Discrete time system

% a)delay of two units
num=[1 0 0];
den=[1 .2 .01];
Ts=0.1;
H1=tf(num,den,Ts);
H2=tf(num,den,Ts,'inputdelay',2)

% b)impulse response of H1(s)
n=0:8;
h=dimpulse(num,den,n);
stairs(n,h)
legend('h[n]')

% c)impulse response of H2(s)
figure
num=1;
h2=dimpulse(num,den,n);
stairs(n,h2)
legend('h[n-2]');


text:
function snew = convertHDF5StrToMatlabStr(s, type)
snew = '';

if ~exist('type', 'var') || isempty(type)
    type = 'char';
end

% Convert muti-row char array to cell string array
if size(s,1)==1 && strcmp(type,'char')
    snew = strtrim_improve(s);
else
    snew = cell(size(s,1),1);
    for ii = 1:size(s,1)
        snew{ii} = strtrim_improve(s(ii,:));
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function embedding = optimize_layout5(head_embedding, tail_embedding, ...
    head, tail, n_epochs, n_vertices, epochs_per_sample, a, b,...
    gamma, initial_alpha, negative_sample_rate, verbose, min_dist)
%OPTIMIZE_LAYOUT2 Another version of optimize_layout.m that slightly
% optimizes the speed for MATLAB; in particular, several quantities are
% vectorized. This function is only called if the UMAP method is 'MATLAB
% vectorized'.
%
% embedding = OPTIMIZE_LAYOUT2(head_embedding, tail_embedding, head, tail,
% n_epochs, n_vertices, epochs_per_sample, a, b)
%
% Parameters
% ----------
% head_embedding: array of size (n_samples, n_components)
%     The initial embedding to be improved by SGD.
%
% tail_embedding: array of size (source_samples, n_components)
%     The reference embedding of embedded points. If not embedding new
%     previously unseen points with respect to an existing embedding this
%     is simply the head_embedding (again); otherwise it provides the
%     existing embedding to embed with respect to.
%
% head: array of size (n_1_simplices, 1)
%     The indices of the heads of 1-simplices with non-zero membership.
%
% tail: array of size (n_1_simplices, 1)
%     The indices of the tails of 1-simplices with non-zero membership.
%
% n_epochs: double
%     The number of training epochs to use in optimization.
%
% n_vertices: double
%     The number of vertices (0-simplices) in the dataset.
%
% epochs_per_samples: array of size (n_1_simplices, 1)
%     A double value of the number of epochs per 1-simplex. 1-simplices with
%     weaker membership strength will have more epochs between being sampled.
%
% a: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% b: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% gamma: double (optional, default 1)
%     Weight to apply to negative samples.
%
% initial_alpha: double (optional, default 1)
%     Initial learning rate for the SGD.
%
% negative_sample_rate: double (optional, default 5)
%     Number of negative samples to use per positive sample.
%
% verbose: boolean (optional, default false)
%     Whether to report information on the current progress of the algorithm.
%
% Returns
% -------
% embedding: array of size (n_samples, n_components)
%     The optimized embedding.
%
% See also: OPTIMIZE_LAYOUT
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
    if nargin < 14
        min_dist = 0.3;
        if nargin < 13
            verbose = false;
            if nargin < 12
                negative_sample_rate = 5;
                if nargin < 11
                    initial_alpha = 1;
                    if nargin < 10
                        gamma = 1;
                    end
                end
            end
        end
    end
    [n_samples, dim] = size(head_embedding);
    ONES=ones(1, dim);
    BG2S=2*gamma* b;
    FOURS=4*ONES;
    ABNEG2=-2.0*a*b;
    BNEG1=b-1;

    same_embedding = isequal(head_embedding, tail_embedding);
    alpha = initial_alpha;

    n_edges = size(epochs_per_sample, 1);
    weights = ones(n_edges,1)./epochs_per_sample;

%     full_weights = full(sparse(double(head), double(tail), double(weights), n_samples, size(tail_embedding, 1)));
%     if same_embedding
%         full_weights = full_weights + eye(n_samples);
%     end
%     repulsion_weights = 1-full_weights;
%     [rep_head, rep_tail, rep_weights] = find(repulsion_weights);
    rep_head = head;
    rep_tail = tail;
    rep_weights = 1-weights;
    rep_epochs = ones(n_edges,1)./(rep_weights);

    %repulsion_epochs = ones(size(rep_weights,1),1)./(rep_weights);
    epochs_per_sample = vertcat(epochs_per_sample, rep_epochs);
    aughead = vertcat(head, rep_head);
    augtail = vertcat(tail, rep_tail);

    epochs_per_negative_sample = epochs_per_sample(1:n_edges) / single(negative_sample_rate);
    epoch_of_next_negative_sample = epochs_per_negative_sample;
    epoch_of_next_sample = epochs_per_sample;

    %tc=tic;

            dists = sqrt(sum((head_embedding(head,:) - tail_embedding(tail,:)).^2, 2));
            %spread = 1; %Should pass in spread from UMAP object.

            result = cross_entropy(double(head_embedding), double(tail_embedding), double(head), double(tail), double(weights), a, b, same_embedding);

            disp(['The current cross entropy is ' num2str(result)]);

    for n = 1:n_epochs
        l=epoch_of_next_sample <= n;
        nl = epoch_of_next_sample(1:n_edges) <= n;
        idxs= find(l);
        N=length(idxs);
        nNegSamples = floor((single(n) - epoch_of_next_negative_sample(nl)) ./ epochs_per_negative_sample(nl));
        mxNeg=max(nNegSamples);
        randy=int32(randi(n_vertices, N, mxNeg));
        next=epoch_of_next_sample(l) + epochs_per_sample(l);
        next_neg=epoch_of_next_negative_sample(nl)+(nNegSamples .* epochs_per_negative_sample(nl));
        %Can NOT vectorize dist_squared OUTSIDE next for loop since head_embedding and tail_embedding
        %change WITHIN next for loop since head and tail contain repeating values
        %   dist_squareds=sqrt(sum((head_embedding(head(l),:)-tail_embedding(tail(l),:)).^2, 2)).^2;

        for lIdx = 1:N
            i=idxs(lIdx);

            j = aughead(i);
            k = augtail(i);

            current = head_embedding(j,:);
            other = tail_embedding(k,:);

            dist_squared = norm(current - other).^2;
            %Can NOT vectorize dist_squared OUTSIDE this loop
            %  assert(dist_squareds(li)==dist_squared);
            if i <= n_edges
                grad_coeff=(dist_squared > 0)*(ABNEG2*(dist_squared).^BNEG1)./ (a*dist_squared.^b + 1);
            else
                grad_coeff=(dist_squared > 0)*BG2S./(0.001 + dist_squared)./(a*dist_squared.^b + 1);
            end

            grad_coeff(isnan(grad_coeff)) = 0;

            %optimize clip()->grad = clip(grad_coeff .* (current - other));
            grad= max(-4, min(4, grad_coeff .* (current - other)));
            %grad = grad_coeff.* (current - other);
            current = current + grad * alpha;
            head_embedding(j,:) = current;
            if same_embedding
                tail_embedding(j,:) = current;
                other = other - grad * alpha;
                head_embedding(k,:) = other;
                tail_embedding(k,:) = other;
            end

            %optimize by vectorizing->epoch_of_next_sample(i) = epoch_of_next_sample(i) + epochs_per_sample(i);

            %optimize by vectorizing->n_neg_samples = floor((n - epoch_of_next_negative_sample(i)) / epochs_per_negative_sample(i));
            if i <= n_edges
                n_neg_samples=nNegSamples(lIdx);
                for p = 1:n_neg_samples
                    %k = randi(n_vertices);
                    k=randy(lIdx,p);
                    if same_embedding && j == k
                        continue
                    end

                    other = tail_embedding(k,:);

                    dist_squared = norm(current - other).^2;

                    %optimize->grad_coeff = (dist_squared > 0)*2*gamma* b*ones(1, dim)./(0.001 + dist_squared)./(a*dist_squared.^b + 1);
                    grad_coeff=(dist_squared > 0)*BG2S./(0.001 + dist_squared)./(a*dist_squared.^b + 1);
                    grad_coeff(isnan(grad_coeff)) = 0;
                    grad = grad_coeff.* (current - other);

                    %optimize-->grad = (grad_coeff > 0).*clip(grad_coeff .* (current - other)) + ~(grad_coeff > 0)*4.*ones(1, dim);
                    %grad=(grad_coeff > 0).*max(-4, min(4, grad_coeff .* (current - other))) + ~(grad_coeff > 0).*FOURS;

                    %the ABOVE statement remains faster than this optimization
                    %if any(grad_coeff <= 0)
                    %    grad=(grad_coeff > 0).*max(-4, min(4, grad_coeff .* (current - other))) + ~(grad_coeff > 0).*FOURS;
                    %else
                    %    grad=max(-4, min(4, grad_coeff .* (current - other)));
                    %end
                    current = current + grad * alpha;
                end

                head_embedding(j,:) = current;
                if same_embedding
                    tail_embedding(j,:) = current;
                end

                %optimize->epoch_of_next_negative_sample(i) = epoch_of_next_negative_sample(i)+(n_neg_samples * epochs_per_negative_sample(i));
            end
        end
        epoch_of_next_sample(l)=next;
        epoch_of_next_negative_sample(nl)=next_neg;
        alpha = initial_alpha * (1 - single(n)/single(n_epochs));

        progress_checkpoint = min(floor(n_epochs / 10), 50);
        if verbose && mod(n, progress_checkpoint) == 0
            fprintf('\t%d/%d epochs done\n', int32(n), int32(n_epochs));

            dists = sqrt(sum((head_embedding(head,:) - tail_embedding(tail,:)).^2, 2));
            %spread = 1; %Should pass in spread from UMAP object.

            result = cross_entropy(double(head_embedding), double(tail_embedding), double(head), double(tail), double(weights), a, b, same_embedding);

            disp(['The current cross entropy is ' num2str(result)]);
        end
    end

    embedding = head_embedding;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function f=loglikGaP(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

[Hkt, cx, cy, Wxk] = loglikGaPreadparam(x,varargin);
Hkt=max(Hkt,eps); % adjust small values to avoid undeflow

[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
P=Wxkbg*Hktbg; %current approximation

%Poisson contribution
t1=Vxt.*log(P) - P;
%Gamma contribution
t2=(alpha-1)*log(Hkt)-1/beta*Hkt-alpha*log(beta)-log(gamma(alpha));

f=sum(t1(:));%+sum(t2(:)));
f=-f; %conjugate gradient is mimimizing!
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function outModel=removeDrugReactions(inModel)
% outModel=removeDrugReactions(inModel)
%
% Removes reactions from Human-GEM related to drugs
%
% inModel       The model, typically Human-GEM or Mouse-GEM
% outModel      The resulting model.

drugNames = {'pravastatin';'Gliclazide';'atorvastatin';'fluvastatin';'fluvastain';'fluvstatin';'simvastatin';'cyclosporine';...
             'acetaminophen';'cerivastatin';'Tacrolimus';'ibuprofen';'lovastatin';'Losartan';'nifedipine';...
             'pitavastatin';'rosuvastatin';'Torasemide';'Midazolam'};
drugMets = contains(inModel.metNames,drugNames,'IgnoreCase',true);
drugRxns = (sum(inModel.S(drugMets,:) ~= 0,1) > 0).';
outModel = removeReactions(inModel,inModel.rxns(drugRxns));
end

author: SysBioChalmers

repo: Human-GEM


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  June 2013
% ------------------------------------------------------------------------

function [stats, extra] = join_masks(stats1, stats2)
stats.gt_set = stats1.gt_set;

% num_objects
stats.num_objects = stats1.num_objects;
if (stats2.num_objects~=stats1.num_objects)
    error('Incoherent number of number of objects')
end
stats.obj_classes = stats1.obj_classes;

% all_n_masks: [2913x1 double]
if size(stats1.all_n_masks,2)~=size(stats2.all_n_masks,2)
    error('Incoherent number of results')
end
n_results = size(stats2.all_n_masks,2);

stats.all_n_masks = stats1.all_n_masks + stats2.all_n_masks;

stats.max_J = zeros(stats.num_objects, n_results);
stats.max_fp = zeros(stats.num_objects, n_results);
stats.max_fn = zeros(stats.num_objects, n_results);
stats.max_inters = zeros(stats.num_objects, n_results);
extra.which = zeros(stats.num_objects, n_results);
for ii=1:n_results
    [stats.max_J(:,ii), which] = max([stats1.max_J(:,ii), stats2.max_J(:,ii)],[],2);
    extra.which(:,ii) = which;
    for jj=1:stats.num_objects
        if which(jj)==1
            stats.max_fp(jj,ii) = stats1.max_fp(jj,ii);
            stats.max_fn(jj,ii) = stats1.max_fn(jj,ii);
            stats.max_inters(jj,ii) = stats1.max_inters(jj,ii);
        else
            stats.max_fp(jj,ii) = stats2.max_fp(jj,ii);
            stats.max_fn(jj,ii) = stats2.max_fn(jj,ii);
            stats.max_inters(jj,ii) = stats2.max_inters(jj,ii);
        end
    end
end

% Recompute statistics
stats.jaccard_object = mean(stats.max_J);
stats.mean_n_masks   = mean(stats.all_n_masks);

class_ids = unique(stats.obj_classes);

% Compute per-class statistics and then mean
for ii=1:length(class_ids)
    curr_class = class_ids(ii);
    stats.per_class_results{ii}.num_objects = sum(stats.obj_classes==curr_class);

    stats.per_class_results{ii}.max_fp     = stats.max_fp(logical(stats.obj_classes==curr_class),:);
    stats.per_class_results{ii}.max_fn     = stats.max_fn(logical(stats.obj_classes==curr_class),:);
    stats.per_class_results{ii}.max_inters = stats.max_inters(logical(stats.obj_classes==curr_class),:);
    stats.per_class_results{ii}.max_J      = stats.max_J(logical(stats.obj_classes==curr_class),:);
    stats.per_class_results{ii}.meanmax    = mean(stats.per_class_results{ii}.max_J);

    stats.per_class_results{ii}.global_fp      = sum(stats.per_class_results{ii}.max_fp);
    stats.per_class_results{ii}.global_fn      = sum(stats.per_class_results{ii}.max_fn);
    stats.per_class_results{ii}.global_inters  = sum(stats.per_class_results{ii}.max_inters);

    % Compute per-class global J
    stats.per_class_results{ii}.global_J = ...
         stats.per_class_results{ii}.global_inters ./...
        (stats.per_class_results{ii}.global_inters+...
         stats.per_class_results{ii}.global_fp    +...
         stats.per_class_results{ii}.global_fn);

     % NaN means mask has no pixels --> P=1, R=0, F=0
     nan_mark = isnan(stats.per_class_results{ii}.global_J);
     stats.per_class_results{ii}.global_J(nan_mark) = 0;
end

% Compute global mean
tmp = [];
for ii=1:length(class_ids)
    tmp = [tmp; stats.per_class_results{ii}.global_J]; %#ok<AGROW>
end

stats.jaccard_class = mean(tmp);

author: s-gupta

repo: rcnn-depth

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


text:
function [lp,dlp] = priorTMulti(mu,s2,nu,x)

% Multivariate Laplace hyperparameter prior distribution.
% Compute log-likelihood and its derivative or draw a random sample.
% The prior distribution is parameterized as:
%
%   p(x) = ( 1 + r2/(nu-2) )^(-(nu+D)/2) / Z, where
%      Z = gamma(nu/2) * sqrt(det(pi*(nu-2)*s2)) / gamma((nu+D)/2), and
%  r2(x) = (x-mu)'*inv(s2)*(x-mu),
%
% further mu(Dx1) is the mean parameter, s2(Dx1) or s2(DxD) is the variance
% parameter, nu(1x1) > 2 is the degrees of freedom parameter and x(DxN) contains
% query hyperparameters for prior evaluation.
%
% For more help on design of priors, try "help priorDistributions".
%
% Copyright (c) by Hannes Nickisch, 2014-10-16.
%
% See also PRIORDISTRIBUTIONS.M, PRIORLAPLACE.M.

if nargin<3, error('mu, s2 and nu parameters need to be provided'), end
if ndims(mu)~=2 || size(mu,2)~=1, error('mu needs to be (Dx1)'), end
D = size(mu,1);
s2_ok = ndims(s2)==2 && all(size(s2)==[D,1] | size(s2)==[D,D]);
if ~s2_ok, error('s2 needs to be (DxD) or (Dx1)'), end
if size(s2,2)==D                                        % full multivariate case
  s = chol(s2)'; lds = sum(log(diag(s)));                 % lds = log(det(s2))/2
else                                                       % diagonal covariance
  s = sqrt(s2);  lds = sum(log(s));
end
if ~isscalar(nu), error('nu parameter needs to be a scalar'), end

if nargin<4                                             % return a random sample
  lp = randn(D,1);
  for d=1:D, lp(d) = sqrt((nu-2)/priorGamma(nu/2,2)) * lp(d); end  % unit sample
  if size(s,2)==D, lp = s*lp+mu; else lp = s.*lp+mu; end % affine transformation
  return
end
if D==1 && size(x,1)>1                            % both mu/s2 scalar => inflate
  D = size(x,1); mu = mu*ones(D,1); s = s*ones(D,1); lds = D*lds;
end
if ~(ndims(x)==2 && size(x,1)==D), error('x needs to be (Dxn)'), end

oN = ones(1,size(x,2)); oD = ones(D,1);
if size(s,2)==D
  xm = x-mu*oN; xs = s\xm;
else
  xm = x-mu*oN; xs = xm./(s*oN);
end

lZ = gammaln((nu+D)/2) - gammaln(nu/2) - D*log(pi*(nu-2))/2 - lds;
r2 = sum(xs.^2,1);
lp = -(nu+D)/2*log(1+r2/(nu-2)) + lZ;
dlp = -(nu+D)/(nu-2)*xs./(oD*(1+r2/(nu-2)));

if size(s,2)==D, dlp = s'\dlp; else dlp = dlp./(s*oN); end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ out datInfo] = opt_smt2cnt(Dat)
%OPT_SMT2CNT Summary of this function goes here
%   Detailed explanation goes here
if isstruct(Dat)
    out=Dat;
    if isfield(Dat, 'x')
        dat=Dat.x;
        if ndims(dat)==3  % smt2cnt
            [nD nTr nCH]=size(dat);
            dat=reshape(dat, [nD*nTr nCH]);
            out.x=dat;
            out.dimInfo=[nD nTr nCH];
        elseif ndims(dat)==2 %cnt2smt
            warning('OpenBMI: check the dimensions of Dat.x(DATA X TRIAL X CH)');
        end
    else
        warning('OpenBMI: check the field of x from input data: Dat.x');
    end
else
    dat=Dat;
    if ndims(dat)==3  % smt2cnt
        [nD nTr nCH]=size(dat);
        dat=reshape(dat, [nD*nTr nCH]);
        out=dat;
        datInfo=[nD nTr nCH];
    else
        warning('OpenBMI: check the dimensions of Dat.x, 2(DATA X CH) or 3(DATA X TRIAL X CH) are acceptable');
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function node2 = grOppositeNode(edges, node)
%GROPPOSITENODE Return opposite node in an edge
%
%   NODE2 = grOppositeNode(EDGE, NODE);
%   Return the index of the node opposite to NODE in EDGE.
%   If the edge does not contain node NODE, result is 0.
%
%   Works also if EDGE is a N-by-2 array of source and target vertex
%   indices, in this case the result NODE2 has as many rows as the number
%   of edges.
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-09-07
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY

% allocate memory for result
Ne = size(edges, 1);
node2 = zeros(Ne, 1);

% iterate on edges
for i = 1:Ne
    if edges(i,1) == node
        node2(i) = edges(i,2);
    elseif edges(i,2) == node
        node2(i) = edges(i,1);
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function IC = arrayToDims( I, dims, padEl )
% Pads or crops I appropriately so that size(IC)==dims.
%
% For each dimension d, if size(I,d) is larger then dims(d) then
% symmetrically crops along d (if cropping amount is odd crops one more
% unit from the start of the dimension).  If size(I,d) is smaller then
% dims(d) then symmetrically pads along d with padEl (if padding amount is
% even then pads one more unit along the start of the dimension).
%
% USAGE
%  IC = arrayToDims( I, dims, [padEl] )
%
% INPUTS
%  I         - n dim array to crop window from (for arrays can only crop)
%  dims      - dimensions to make I
%  padEl     - [0] element with which to pad
%
% OUTPUTS
%  IC        - cropped array
%
% EXAMPLE
%  I=randn(10); delta=1; IC=arrayToDims(I,size(I)-2*delta);
%
% See also ARRAYCROP, PADARRAY
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(padEl)); padEl=0; end
nd = ndims(I);  siz = size(I);
[dims,er] = checkNumArgs( dims, size(siz), 0, 1 ); error(er);
if(any(dims==0)); IC=[]; return; end

% get start and end locations for cropping
strLocs = ones( 1, nd );  endLocs = siz;
for d=1:nd
  delta = siz(d) - dims(d);
  if ( delta~=0 )
    deltaHalf = floor( delta / 2 );
    deltaRem = delta - 2*deltaHalf;
    strLocs(d) = 1 + (deltaHalf + deltaRem);
    endLocs(d) = siz(d) - deltaHalf;
  end
end

% call arrayCrop
IC = arrayCrop( I, strLocs, endLocs, padEl );

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
% hole_example1.m
%
% 2D example showing the "hole in sphere" effect for OSEM SPECT
% and various regularization methods that may help overcome it.
%
% Copyright 2004-5-13, Jeff Fessler, University of Michigan

if ~has_aspire, return, end

if ~isvar('A0'), printm 'generate system matrix'
	f.dir = test_dir;
	f.dsc0 = [f.dir 't0.dsc'];
	f.wtf0 = strrep(f.dsc0, 'dsc', 'wtf');
	if 1 % todo: replace with aspire_pair
		os_run(['wt -chat 0 dsc 12 fwhm_detector 5 >! ' f.dsc0])
		os_run(['echo y | wt -chat 0 gen ' f.dsc0 ' row'])
	end

	A0 = Gtomo2_wtmex(f.wtf0);
	ig = A0.arg.ig;
	sg = A0.arg.sg;
	nx = ig.nx; nb = sg.nb;
	ny = ig.ny; na = sg.na;
end

if ~isvar('Ab0'), printm 'Ab'
	f.nblock = 6;
	Ab0 = Gblock(A0, f.nblock);
	Ab20 = Gblock(A0, 20);
end

if ~isvar('yi'), printm 'projection data'
%	xtrue = [6.5 -0.5 5 5 0 100]; over = {'oversample', 4};
%	xtrue = ellipse_im(ig, xtrue, over{:});
	over = {'oversample', 1};
	xtrue1 = ellipse_im(ig, [18.5 -0.5 5 5 0 100], over{:});
	xtrue2 = ellipse_im(ig, [-3.5 -0.5 3 3 0 100], over{:});
	xtrue3 = ellipse_im(ig, [-18.5 -0.5 2 2 0 100], over{:});
	xtrue = xtrue1 + xtrue2 + xtrue3;
	im(xtrue, 'xtrue'), cbar
prompt

	ri = 1;
	yi = A0 * xtrue + ri;
	im(yi, 'yi'), cbar
prompt
end

% uniform initial image
xinit = ig.ones;

if 0 && ~isvar('xfbp'), printm 'FBP'
	tmp = fbp2(sg, ig);
	xfbp = fbp2(yi-ri, tmp, 'window', 'hann');
	im(xfbp, 'fbp'), cbar
prompt
end

if ~isvar('xe0'), printm 'EM'
	f.niter_em = 5;
%	xe0 = eml_em(xinit(ig.mask), A0, yi(:), 1, ri, [], f.niter_em);
	xe0 = eml_em(xinit(ig.mask), A0, yi(:), 1, ri, ...
		'isave', 'all', 'niter', f.niter_em);
	xe0 = ig.embed(xe0);
	im(xe0, 'EM'), cbar
prompt
end

xinit = xe0(:,:,end); % last EM iter

% OS-EM iterations
if ~isvar('xo0'), printm 'run os-em'
	f.niter = 10; % for testing
	if isempty(caller_name)
		f.niter = 200; % for full comparison
	end

	rri = ri * ones(nb,na);
	xo0 = eml_osem(xinit(ig.mask), Ab0, yi, [], rri, f.niter);
	xo0 = ig.embed(xo0);
	im(xo0(:,:,end), 'OSEM'), cbar
prompt
end

% regularization (edge preserving)
if ~isvar('R'), printm 'R'
	f.l2b = -9;
	f.l2b = -11; % for binary case
	f.potential = 'hyper3'; f.delta = 3;
%	f.potential = 'cauchy'; f.delta = 10; % this one used for grant
	f.potential = 'cauchy'; f.delta = 2;
%	R = Robject(ig.mask, 'type_denom', 'matlab', ...
%		'beta', 2^f.l2b, 'potential', f.potential, 'delta', f.delta);
	R = Reg1(ig.mask, 'type_denom', 'matlab', ...
		'beta', 2^f.l2b, 'pot_arg', {f.potential, f.delta});
	clear xp0
prompt
end

% Incremental EM-3, aka C-OS-3
if ~isvar('xp0'), printm 'C-OS-3'
%	xp0 = epl_os_emdp(xinit(ig.mask), Ab0, yi, ones(size(yi)), rri, R, ...
%		0+1*f.niter, 200, 1);
	xp0 = epl_inc(xinit(ig.mask), Ab20, yi, ones(size(yi)), rri, R, ...
		'niter', f.niter, 'hds', 3, 'os', 40, 'pixmax', 200, 'chat', 0);
	xp0 = ig.embed(xp0);
	im(xp0(:,:,end), 'C-OS-3')
prompt
end

% regularization with perfect side info
if ~isvar('Rs'), printm 'Rs'
	ug = ugibb_form(xtrue > 50, 'threshold', 0.5, 'offsets', R.offsets);
	im(ug)
	f.l2b = -9;
%	Rs = Robject(ig.mask, 'type_denom', 'matlab', ...
%		'beta', 2^f.l2b, 'user_wt', ug);
	Rs = Reg1(ig.mask, 'type_denom', 'matlab', ...
		'beta', 2^f.l2b, 'user_wt', ug);
	clear xs0
prompt
end

if ~isvar('xs0'), printm 'recon with side info'
	xs0 = epl_inc(xinit(ig.mask), Ab20, yi, ones(size(yi)), rri, Rs, ...
		'niter', f.niter, 'hds', 3, 'os', 40, 'pixmax', 200, 'chat', 0);
	xs0 = ig.embed(xs0);
	im(xs0(:,:,end), 'xs0')
prompt
end

if 1 && im, printm 'show "hole in sphere"'
	iy = ny/2;
	t = squeeze(xo0(:,iy,:));
	t = squeeze(xp0(:,iy,:));
	im(t)

	ix = 1:nx;
%	ix = 10:60;
	clf
%	subplot(211)
%	plot(ix, xtrue(ix,iy), 'y-o', 1:nx, t);

	plot(	ix, xtrue(ix,iy), 'y-x', ...
		ix, xo0(ix,iy,end), 'c-o', ...
		ix, xp0(ix,iy,end), 'g.-', ...
		ix, xs0(ix,iy,end), 'r.-')
	legend('True', 'ML: OS-EM', 'PL: C-OS-3', 'PL: side', 4)
	xlabel 'Horizontal Pixel Index'
	ylabel 'Activity [arbitrary units]'
	title 'Profile'
	axisy([0 160])
prompt
end

if im, printm 'show total activity per iteration'
	t = xtrue3 > 0;
	sum_true = sum(xtrue(t));
	t1 = reshape(xo0, [], f.niter+1); t1 = sum(t1(t,:)) / sum_true;
	t2 = reshape(xp0, [], f.niter+1); t2 = sum(t2(t,:)) / sum_true;
	t3 = reshape(xs0, [], f.niter+1); t3 = sum(t3(t,:)) / sum_true;
	ii = 0:f.niter;
	clf, % subplot(212)
	plot( ...
		ii(1), t3(1), 'y-s', ...
		ii(1), t2(1), 'g.-', ...
		ii(1), t1(1), 'c--o', ...
		ii(1:8:end), t3(1:8:end), 'ys', ...
		ii(1:8:end), t2(1:8:end), 'g.', ...
		ii(1:8:end), t1(1:8:end), 'co', ...
		ii, t3, 'y-', ...
		ii, t2, 'g-', ...
		ii, t1, 'c--', ...
		ii, 1, 'y:')
	legend('PL side', 'PL: C-OS-3', 'ML: OS-EM', 4)
%	legend('PL: C-OS-3 with side information', 'ML: OS-EM', 4)
	xlabel 'Iteration', ylabel 'Activity Recovery'
	axis([0 120 0.8 1.02])
	axis([0 120 0.6 1.02])
%	ir_savefig hole_example_recovery1 c
%	set(gcf, 'InvertHardCopy', 'off'); print -djpeg90 hole_example_recovery2
prompt
end

if 1 && im, printm 'examine "hole in sphere"'
	t = xtrue1(:) > 99;
	t1 = reshape(xo0, [nx*ny f.niter+1]);
	t1 = t1(t, :);
	t2 = reshape(xp0, [nx*ny f.niter+1]);
	t2 = t2(t, :);
	ii = 0:f.niter;
	clf, % subplot(212)
	plot(	...
		ii(1), max(t1(:,1), [], 1), 'c-o', ...
		ii(1), min(t1(:,1), [], 1), 'c--x', ...
		ii(1), max(t2(:,1), [], 1), 'g-', ...
		ii(1), min(t2(:,1), [], 1), 'g--', ...
		ii(1:8:end), max(t1(:,1:8:end), [], 1), 'co', ...
		ii(1:8:end), min(t1(:,1:8:end), [], 1), 'cx', ...
...%		ii(1:8:end), max(t2(:,1:8:end), [], 1), 'g-', ...
...%		ii(1:8:end), min(t2(:,1:8:end), [], 1), 'g--', ...
		ii, max(t1, [], 1), 'c-', ...
		ii, min(t1, [], 1), 'c--', ...
		ii, max(t2, [], 1), 'g-', ...
		ii, min(t2, [], 1), 'g--', ...
		ii, 100, 'y:')
	legend('ML: OSEM max', 'ML: OSEM min', 'PL: C-OS-3 max', 'PL: C-OS-3 min')
	xlabel 'Iteration'
	ylabel 'Activity [arbitrary units]'
	title 'Maximum and Minimum Values within ROI'
	axisy([0 165])
prompt
end

if 1, printm 'show reconstructions' % for Ken
	ix = 1:nx; iy=1:ny;
%	ix = 7+[nx/4:(3*nx/4)]; iy=ny/4:(3*ny/4);
	t = stackup(xtrue(ix,iy), xo0(ix,iy,end), xp0(ix,iy,end), xs0(ix,iy,end));
	im clf, im(t, [0 120]), cbar
%	t(t>120) = 120; t = uint8(255 * t/120);
%	imwrite(t', 'hole.tif')
end

% ir_savefig fig_hole_example1_new

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function str = showdetail(item)

% function str = showdetail(item)
% Display details for a cfg_files item.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showdetail.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

str = showdetail(item.cfg_item);
str{end+1} = 'Class  : cfg_files';
str{end+1} = 'A cellstr array of file names.';
str = [str; gencode(item.filter, 'filter:')];
str = [str; gencode(item.num,    'num   :')];

author: spm


text:
function [u_opt_total, x_pred, y_pred, x_pred_left, y_pred_left, x_pred_right, y_pred_right, ...
    vx_pred, d_pred, dot_d_pred, s_dot_pred, ...
    ax_pred, ax_tire_pred, ay_pred, TireUtilizationTarget, cost_values]  = ...
    transformMPCResult(status, u_opt_total, errorState, x_traj, y_traj, ...
    psi_traj, v_traj, ax_diff_traj, d_Target_m, ax_lim_mps2, ay_lim_mps2, ...
    vx_lin_mps, kappa_lin_radpm, UncertaintyTube, ...
    sys, drag_coefficient, roh_air, vehiclemass_kg, P_VDC_PositiveAxLimScale)
% function transformMPCResult
% Author:       Martin Euler
%               Alexander Wischnewski
% Description:
%   function used to process the output of the mpc solver
% Inputs/parameters:
%   status:             Solver status / > 0 -> ok / 0 -> not ok
%   u_opt_total:        optimal solution
%   error_state:        vehicle state in error coordinates
%   UncertaintyTube:    Matrix with all uncertainty matrices for the prediction horizon
%   x_traj:             x-East of target trajectory
%   y_traj:             y-North of target trajectory
%   psi_traj:           Heading of target trajectory
%   v_traj:             Velocity of target trajectory
%   ax_diff_traj:       Velocity derivative with respect to path coordinate of target trajectory
%   d_Target_m:         Target trajectory projected into curvilinear coordinate system
%   ax_lim_mps2:        Limit for longitudinal acceleration
%   ay_lim_mps2:        Limit for lateral acceleration
%   ax_dist_mps2:       Disturbance in longitudinal acceleration
%   ay_dist_mps2:       Disturbance in lateral acceleration
%   vx_lin_mps:         Linearization velocity
%   kappa_lin_radpm:    Linearization curvature
%   sys:                struct with necessary parameters and matrices
%   drag_coefficient:   Coefficient for air drag according to F = 1/2*roh_air*drag_coefficient
%   roh_air:            Air density
%   vehiclemass_kg:     Overall vehicle mass
%   P_VDC_PositiveAxLimScale Scaling factor for positive accelerations to consider RWD

% Outputs:
%   u_opt_total:        optimal solution
%   x_pred:             Predicted x-East positions (mean of uncertainty tube)
%   y_pred:             Predicted y-North positions (mean of uncertainty tube)
%   x_pred_left:        Predicted x-East positions (left side of uncertainty tube)
%   y_pred_left:        Predicted y-North positions (left side of uncertainty tube)
%   x_pred_right:       Predicted x-East positions (right side of uncertainty tube)
%   y_pred_right:       Predicted y-North positions (right side of uncertainty tube)
%   vx_pred:            Predicted velocity
%   d_pred:             Predicted lateral deviation
%   dot_d_pred:         Predicted lateral deviation derivative
%   s_dot_pred:         Predicted path progress
%   ax_pred:            Predicted longitudinal acceleration
%   ax_tire_pred:       Predicted longitudinal acceleration "at tire" (with driving resistances)
%   ay_pred:            Predicted lateral acceleration
%   TireUtilization:    Tire utilization (target)
%   cost_values:        debug information on separate cost terms

%% ------- define optimization problem variables ----------------------- %%
nx = 3; % number of state variables should equal sys.n_sys

% initialize output variables
cost_values = zeros(4, 1);

% check if solver solution is valid
if status > 0
    %% -------------------------- obtain predictions for nominal system ---------------------------- %%
    % calculate errorState prediction and improve computational efficiency by exploiting
    % knowledge that slacks do not contribute to this
    error_pred = sys.ABK_MPC(:, 1:sys.m_sys*(sys.N_hor+1))*u_opt_total(1:sys.m_sys*(sys.N_hor+1)) ...
        + sys.Ax0_MPC*errorState;
    % update velocity predictions and limit to ensure robust transformations
    vx_pred = v_traj - error_pred(1:nx:end);
    % update lateral error predictions
    d_pred = error_pred(2:nx:end);
    % calculate lateral uncertainties
    tight_d_m = sqrt(UncertaintyTube(2, 2:nx:end))';
    % update lateral error derivative predictions
    dot_d_pred = error_pred(3:nx:end);
    % update predictions for path coordinate
    d_psi_pred = asin(max(min(dot_d_pred./vx_pred, 0.999), -0.999));
    s_dot_pred = vx_pred.*cos(d_psi_pred)./(1-d_pred.*kappa_lin_radpm);
    % update position predictions (compensate for linearization error with respect to target path)
    x_pred = x_traj - sin(psi_traj+pi/2).*(d_pred-d_Target_m);
    y_pred = y_traj + cos(psi_traj+pi/2).*(d_pred-d_Target_m);
    % update uncertainty tube predictions
    % +1 for vehicle width
    x_pred_left = x_pred - sin(psi_traj+pi/2).*(tight_d_m+1);
    y_pred_left = y_pred + cos(psi_traj+pi/2).*(tight_d_m+1);
    x_pred_right = x_pred + sin(psi_traj+pi/2).*(tight_d_m+1);
    y_pred_right = y_pred - cos(psi_traj+pi/2).*(tight_d_m+1);
    % update longitudinal acceleration predictions
    ax_pred = ax_diff_traj.*(kappa_lin_radpm.*vx_lin_mps.*d_pred + ...
        2*vx_lin_mps-vx_pred) + ...
        roh_air*drag_coefficient.*vx_lin_mps.*error_pred(1:nx:end)./vehiclemass_kg +...
        - u_opt_total(1:2:2*(sys.N_hor+1));
    % include driving resistances
    ax_tire_pred = ax_pred + 0.5.*roh_air.*drag_coefficient.*vx_lin_mps.^2./vehiclemass_kg;
    % update lateral acceleration predictions
    ay_pred = kappa_lin_radpm.^2.*vx_lin_mps.^2.*d_pred - ax_diff_traj.*dot_d_pred ...
        - 2*kappa_lin_radpm.*vx_lin_mps.*error_pred(1:nx:end) ...
        + 2*kappa_lin_radpm.*vx_lin_mps.*(v_traj - vx_lin_mps) ...
        + kappa_lin_radpm.*vx_lin_mps.^2 + u_opt_total(2:2:2*(sys.N_hor+1));

    % tire utilization target
    if(ax_tire_pred(1) > 0)
        % consider RWD factor for positive accelerations
        TireUtilizationTarget = abs(ax_tire_pred(1)/(P_VDC_PositiveAxLimScale*ax_lim_mps2(1))) ...
            + abs(ay_pred(1)/ay_lim_mps2(1));
    else
        TireUtilizationTarget = abs(ax_tire_pred(1)/(ax_lim_mps2(1))) ...
            + abs(ay_pred(1)/ay_lim_mps2(1));
    end

    %% recalculate costs
    % calculations removed due to computation time restrictions
    cost_values(1) = 0;
    cost_values(2) = 0;
    cost_values(3) = 0;
    cost_values(4) = 0;

else
    % initialize all predictions
    s_dot_pred = 0.01*ones(41, 1);
    x_pred = zeros(41, 1);
    y_pred = zeros(41, 1);
    x_pred_left = zeros(41, 1);
    y_pred_left = zeros(41, 1);
    x_pred_right = zeros(41, 1);
    y_pred_right = zeros(41, 1);
    vx_pred = zeros(41, 1);
    d_pred = zeros(41, 1);
    dot_d_pred = zeros(41, 1);
    ax_pred = zeros(41, 1);
    ax_tire_pred = zeros(41, 1);
    ay_pred = zeros(41, 1);
    TireUtilizationTarget = 0;
end

author: TUMFTM


text:
function X = ipg(A,Y,X,maxiter)
%
% Interior-Point Method
%
% REFERENCE:
%   M. Merritt and Y. Zhang. An Interior-Point Gradient Method for
%   Large-Scale Totally Nonnegative Least Squares Problems. J. Optimization
%   Theory and Applications, Vol. 126. No. 1, pp. 191-202, July, 2005.
%

% initialize
B = A'*A;
BX = B*X;

for iter = 1:maxiter

   % step computation
    R = A*X - Y;
    S = A'*R;
    P = (X./max(BX,eps)).*S;
    BP = B*P;

    % steplength
     alpha1 = sum(P.*S,1)./max(sum(P.*BP,1),eps);
     alpha2 = -1./min(-eps,min(-P./max(X,eps),[],1) );
     alpha = min([.99*alpha2; alpha1],[],1);

    % update
    X = X - P.*repmat(alpha,size(B,2),1);
    BX = BX - BP.*repmat(alpha,size(B,2),1);

end


text:
function y=upsampleZero(x,P)
%%UPSAMPLEZERO Increase the sampling rate of a signal by inserting zeros
%               in the new samples. This is often used as a first step in
%               interpolatory filters. Adding zeros between samples causes
%               the frequency response of the fft to be aliased. One can
%               then apply a window in the frequency domain and take the
%               ifft to get an interpolated signal. Just adding zeros
%               results in an interpolation with a lot of Gibbs phenomenon.
%
%INPUTS: x An NX1 or 1XN vector or an NXnumVec matrix to upsample. If it is
%          a matrix, then the values are taken per column to upsample. If
%          it is a vector, then the output will have the same orientation
%          as the input (column or row).
%        P The positive integer factor to upsample. P=1 returns the original
%          signal.
%
%OUTPUTS: y The upsampled signal. If x is a vector, it will be (N*P)X1 or
%           1X(N*P), matcghing the input orientation. If a matrix, it will
%           be (N*P)XnumVec in size.
%
%EXAMPLE:
% x=1:5;
% P=3;%So there will be two zeros between each number.y=
% y=upsampleZero(x,P)
%
%December 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

isTranspose=false;
if(isvector(x))
    if(size(x,2)>size(x,1))
        isTranspose=true;
        x=x(:);
    end
end

N=size(x,1);
numVals=size(x,2);
y=zeros(N*P,numVals);
y(1:P:(N*P),:)=x;

if(isTranspose)
    y=y.';
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:
classdef (Abstract) Tx  < adi.common.RxTx & adi.common.DDS
    % Tx: Common shared functions between transmitter classes

    methods (Hidden, Access = protected)

        function numIn = getNumInputsImpl(obj)
            if strcmp(obj.DataSource,'DDS')
                numIn = 0;
            else
                numIn = 1;
            end
        end

    end

    methods (Hidden, Access = protected)

        function setupImpl(obj,data)
            if strcmp(obj.DataSource,'DMA')
                obj.SamplesPerFrame = size(data,1);
            end
            % Call the superclass method
            setupImpl@adi.common.RxTx(obj);
        end

        function valid = stepImpl(obj,dataIn)
            % valid = tx(data) returns a logical value that indicates
            % data was sent to the device correctly.
            %
            % When 'DataSource' is 'DMA' the input 'data' will be an [NxM]
            % vector where N is the length of the input data and M is the
            % number of elements in 'EnabledChannels'. 'data' should be
            % complex if the device assumes complex data.
            %
            % When 'DataSource' is 'DDS' the operator will take no inputs.
            % Running the operator simply will for a DDS settings update
            % and connect to the device if never run before.

            if strcmp(obj.DataSource,'DMA')
                if obj.ComplexData
                    % Interleave channels include I and Q
                    c = obj.channelCount;
                    s = size(dataIn);
                    assert(s(2)==c/2,sprintf('Data size must [Nx%d]\n',c/2));
                    index = 1;
                    %
                    outputData = complex(zeros(length(dataIn)*c,1));
                    for k = 1:2:c
                        outputData(k+0:c:end,1) = real(int16(dataIn(:,index).'));
                        outputData(k+1:c:end,1) = imag(int16(dataIn(:,index).'));
                        index = index + 1;
                    end

                else
                    c = obj.channelCount;
                    s = size(dataIn);
                    assert(s(2)==c,sprintf('Data size must [Nx%d]\n',c));
                    %
                    outputData = zeros(length(dataIn)*c,1);
                    for k = 1:c
                        outputData(k:c:end,1) = int16(dataIn(:,k).');
                    end
                end
                %%
                valid = sendData(obj,outputData);
            else
                valid = true;
                obj.DDSUpdate();
            end
        end

    end

end

author: analogdevicesinc


text:
function []=panel4hddisp(N,n,T,Units,endo,hd_estimates,stringdates1,decimaldates1,pref)

% plot the figure
if pref.plot
    hd=figure('Tag','BEARresults');
    set(hd,'Color',[0.9 0.9 0.9]);
    set(hd,'name','historical decomposition');
    % initiate the count
    count=0;
    % loop over units
    for ii=1:N
        % loop over endogenous variables
        for jj=1:n
            % loop over shocks
            for kk=1:n+1
                % increment count
                count=count+1;
                % then plot
                subplot(N*n,n+1,count)
                temp=hd_estimates{jj,kk,ii};
                hold on
                Xpatch=[decimaldates1' fliplr(decimaldates1')];
                Ypatch=[temp(1,:) fliplr(temp(3,:))];
                HDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
                set(HDpatch,'facealpha',0.5);
                set(HDpatch,'edgecolor','none');
                plot(decimaldates1,temp(2,:)','Color',[0.4 0.4 1],'LineWidth',2);
                plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
                hold off
                set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
                % top labels
                if (jj==1 && kk<=n)
                    title([Units{ii,1} '\_' endo{kk,1}],'FontWeight','normal');
                elseif (jj==1 && kk==n+1)
                    title('Exogenous','FontWeight','normal');
                end
                % side labels
                if kk==1
                    ylabel([Units{ii,1} '\_' endo{jj,1}],'FontWeight','normal');
                end
            end
        end
    end
    % top supertitle
    ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
    set(get(ax,'Title'),'Visible','on')
    title('Contribution of shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
    % side supertitle
    ylabel('Fluctuation of:','FontSize',12,'FontName','Times New Roman','FontWeight','normal');
    set(get(ax,'Ylabel'),'Visible','on')
end

% save on Excel

% create the cell that will be saved on excel
hdcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},2,5*(n+1));
% loop over units
for ii=1:N
    % loop over endogenous variables
    for jj=1:n
        % initiate the cell of results
        endocell={};
        % loop over shocks
        for kk=1:n
            % create a header
            header=[{['fluctuation of ' Units{ii,1} '_' endo{jj,1} ' due to ' Units{ii,1} '_' endo{kk,1} ' shocks']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
            % complete the cell
            tempcell=[[header;stringdates1 num2cell((hd_estimates{jj,kk,ii}(1,:))') num2cell((hd_estimates{jj,kk,ii}(2,:))') num2cell((hd_estimates{jj,kk,ii}(3,:))')] vertspace];
            % concatenate to the previous parts of unitcell
            endocell=[endocell tempcell];
        end
        % exogenous contribution
        header=[{['fluctuation of ' Units{ii,1} '_' endo{jj,1} ' due to exogenous']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
        % complete the cell
        tempcell=[[header;stringdates1 num2cell((hd_estimates{jj,n+1,ii}(1,:))') num2cell((hd_estimates{jj,n+1,ii}(2,:))') num2cell((hd_estimates{jj,n+1,ii}(3,:))')] vertspace];
        % concatenate to the previous parts of unitcell
        endocell=[endocell tempcell];
        % concatenate to the previous parts of irfcell
        hdcell=[hdcell;horzspace;endocell];
    end
    hdcell=[hdcell;horzspace];
end
% trim
hdcell=hdcell(3:end-2,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),hdcell,'hist decomposition','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
% display() - display an EEG data class underlying structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function display(obj);

    tmpMMO = memmapfile(obj.dataFile, 'writable', obj.writable, 'format', { 'single' obj.dimensions 'x' });
    if obj.transposed, disp('Warning: data does not display properly for memory mapped file which have been transposed'); end;
    disp(tmpMMO.data.x);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ihs, x_mu, y_mu, sig] = plotreswh(res, peval, dpixc, p, savethis, tilesonly, showbg, showpoints, shownumbers)
% function [ihs, x_mu, y_mu, sig] = plotreswh(res, peval, dpixc, p, savethis, tilesonly, showbg, showpoints, shownumbers)
% showbg = 0;
% showpoints = 1;
% savethis = 0;
% tilesonly = 1;
if ~(showbg == 1)
    res.w = res.w(:,1:end-1);
    res.h = res.h(1:end-1,:);
end

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

x_mu = [];
y_mu = [];
sig = [];
ca
% mh = mean (res.h,2);
% [mhs,ihs] = sort(mh,'descend');
% sm = sum(res.w.^2,1);
% [mhs, ihs] = sort(sm, 'descend');
[mhs, ihs] = sortcomponents(res.w);
% [mhs, ihs] = sortcomponents(res.h);
% ts = testimportance(reshape(dpixc, peval.nx*peval.ny, peval.nt), res.w, res.h);
% [mhs, ihs] = sort(ts, 'descend');

hsort=res.h(ihs,:);
wsort = res.w(:,ihs);
if showbg > 1
    hsort = hsort(1:showbg,:);
    wsort = wsort(:,1:showbg);
end

wr = reshape(wsort,peval.nx,peval.ny,size(wsort,2));

if savethis
    mkdir ('w')
end

sw = size(wr,3);
a = ceil(sqrt(sw));
b = ceil(sw/a);
maxd = max(dpixc,[],3);
maxd = maxd/max(maxd(:));
mwr = mean(squeeze(max(max(wr,[],1),[],2)));
for ii=1:sw
    [x_mu(ii), y_mu(ii), sig(ii), differ(ii)] = fitgauss2d(wr(:,:,ii));
end

if ~tilesonly
    for ii=1:sw
        %     dipshow(ii,wr(:,:,ii))
        %     colormap('jet')
        ims(wr(:,:,ii),'gray',1);
        switch showpoints
            case 0
                if savethis; SaveImageFULL(['w/w_' num2str(ii)]); end
            case 1
                hold on
                scatter(p.x_vec+.5, p.y_vec+.5,'r')
                hold off
                if savethis; SaveImageFULL(['w/w2_' num2str(ii)]); end
            case 2
                %                 mwr = max(max(wr(:,:,ii)));
                hold on
                contour(mwr*maxd)
                hold off
                if savethis; SaveImageFULL(['w/w2_' num2str(ii)]); end
            case 3
                %                 mwr = max(max(wr(:,:,ii)));
%                 [x_mu(ii), y_mu(ii), sig(ii), differ(ii)] = fitgauss2d(wr(:,:,ii));
                hold on
                if ~isempty(p)
                    scatter(p.x_vec+.5, p.y_vec+.5,'r')
                end
                scatter(x_mu(ii), y_mu(ii),sig(ii)*100, 'g')
                hold off
                if savethis; SaveImageFULL(['w/w2_' num2str(ii)]); end
        end
        %     figure(sw+1)
        %     subplot(a,b,ii)
        %     imagesc(wr(:,:,ihs(ii)))
        %     set(gca,'dataaspectratio',[1 1 1])
        %     set(gca,'xtick',[])
        %     set(gca,'ytick',[])
        %     xlabel(num2str(ii))
    end
end

figure ('name','w tiled');
imstiled(wr,[],'gray',1)
suptitle('w - ordered')
if savethis; SaveImageFULL('w_tiled'); end
% if showpoints
%     for ii=1:size(res.w,2)-1
%         dipshow(ii,wr(:,:,ihs(ii+1)))
%         hold on
%         scatter(p.x_vec-.5, p.y_vec-.5,'r')
%         hold off
%         SaveImageFULL(['w/w2_' num2str(ii)])
%     end
% end

figure ('name','H');
imagesc(hsort)
if ~showbg; ylim([0.5 peval.ncomp+0.5]); end
xlabel('slice # (time)')
ylabel('component')
title('H')
if savethis; SaveImageFULL(['Himage']); end

if tilesonly > 1
    cc = corrcoef(hsort');
    figure ('name','corrcoef(H'')');
    ims(cc)
    colorbar
    title('corrcoef(H'')')
    if savethis; SaveImageFULL('Hcorrel'); end

    figure ('name','abs(corrcoef(H''))');
    imagesc(abs(cc))
    set(gca, 'DataAspectRatio',[1 1 1])
    colorbar
    title('abs(corrcoef(H''))')
    if savethis; SaveImageFULL('Hcorrelabs'); end

    figure ('name','mean(H,2)');
    meanh = mean(hsort,2);
    bar(meanh)
    xlabel('component')
    ylabel('mean (H,2)')
    grid on
    title('sum(H,2)')
    % if ~showbg; xlim([1.5 peval.ncomp+1.5]); end
    if savethis; SaveImageFULL('Hmean_bar'); end

    figure ('name','corrcoef(W)');
    ims(corrcoef(wsort))
    colorbar
    title('corrcoef(W)')
    if savethis; SaveImageFULL('Wcorrel'); end
end
figure ('name','sum');
ims(sum(dpixc,3),'gray')
hold on
if ~isempty(p)
    scatter(p.x_vec+0.5, p.y_vec+0.5,'r')
end
scatter(x_mu, y_mu,sig.^2*100, 'b')
if shownumbers
    text(x_mu+.2, y_mu, num2cell([1:sw]) )
end
% scatter(x_mu, y_mu,sig.^2*100, [1:sw])
if savethis; SaveImageFULL('Sum'); end


text:
function M_hat = run_mc(params)
  M = params.M;
  Idx = params.Idx;
  warning('off','all');
  [numr,numc] = size(M);
  rank = 10;
  Known = find(Idx);
  data = M(Known);
  [U,S,V] = svp(Known,data,numr,numc,rank);
  M_hat = U*diag(S)*V';
  warning('on','all');
end

author: andrewssobral


text:
% <<<3step>>> Merge data from all subjects into one file
% -----------------------------------------------------------------------------------
% 1) loading files
% 2) merging
% 3) changing demensions
% -----------------------------------------------------------------------------------

fpath = ('D:\anedata_step\');

for i = 1:10

    eval(['load d2_MM_' num2str(i) '.mat']);

end

MM_c = cat(3, d2_MM_1, d2_MM_2, d2_MM_3, d2_MM_5, d2_MM_6, d2_MM_7, d2_MM_8, d2_MM_9, d2_MM_10);

MM1 = MM_c(:, :, 1);
MM2 = MM_c(:, :, 2);
MM3 = MM_c(:, :, 3);
MM4 = MM_c(:, :, 4);
MM5 = MM_c(:, :, 4);
MM6 = MM_c(:, :, 5);
MM7 = MM_c(:, :, 6);
MM8 = MM_c(:, :, 7);
MM9 = MM_c(:, :, 8);
MM10 = MM_c(:, :, 9);

MM_all = cat(2, MM1, MM2, MM3, MM5, MM6, MM7, MM8, MM9, MM10);

save MM_3step_exceptMM4


text:
function [c,ceq,dc,dceq]=fun11(x);
c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];
dc=[x(2)-1,-x(2);x(1)-1,-x(1)];
ceq=[];dceq=[];

author: Eurus-Holmes


text:
function b = r8sr_mxv ( n, nz, row, col, diag, off, x )

%*****************************************************************************80
%
%% R8SR_MXV multiplies a R8SR matrix times a vector.
%
%  Discussion:
%
%    The R8SR storage format stores the diagonal of a sparse matrix in DIAG.
%    The off-diagonal entries of row I are stored in entries ROW(I)
%    through ROW(I+1)-1 of OFF.  COL(J) records the column index
%    of the entry in A(J).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer NZ, the number of offdiagonal nonzero elements in A.
%
%    Input, integer ROW(N+1).  The nonzero offdiagonal elements of row I of A
%    are contained in A(ROW(I)) through A(ROW(I+1)-1).
%
%    Input, integer COL(NZ), contains the column index of the element
%    in the corresponding position in A.
%
%    Input, real DIAG(N), the diagonal elements of A.
%
%    Input, real OFF(NZ), the off-diagonal elements of A.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product A * X.
%
  b(1:n) = diag(1:n) .* x(1:n);

  for i = 1 : n
    for k = row(i) : row(i+1)-1
      j = col(k);
      b(i) = b(i) + off(k) * x(j);
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% TEST_MAXWELL_RING_G_NMNN: data function for Neumann boundary condition.

function g = test_maxwell_ring_g_nmnn (x, y, ind)

  [theta, r] = cart2pol (x, y);
  g = zeros ([2, size(x)]);
  switch (ind)
    case 1
      g(1,:,:) = sin(theta) .* (cos(x) - cos(y));
      g(2,:,:) = cos(theta) .* (cos(y) - cos(x));
    case 2
      g(1,:,:) = sin(theta) .* (cos(y) - cos(x));
      g(2,:,:) = cos(theta) .* (cos(x) - cos(y));
    case 3
      g(1,:,:) = cos(x) - cos(y);
    case 4
      g(2,:,:) = cos(y) - cos(x);
    otherwise
      error ('g_nmnn: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function a_out = r8vec_heap_d ( n, a )

%*****************************************************************************80
%
%% R8VEC_HEAP_D reorders an R8VEC into an descending heap.
%
%  Discussion:
%
%    A descending heap is an array A with the property that, for every index J,
%    A(J) >= A(2*J) and A(J) >= A(2*J+1), (as long as the indices
%    2*J and 2*J+1 are legal).
%
%  Diagram:
%
%                  A(1)
%                /      \
%            A(2)         A(3)
%          /     \        /  \
%      A(4)       A(5)  A(6) A(7)
%      /  \       /   \
%    A(8) A(9) A(10) A(11)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the size of the input array.
%
%    Input, real A(N), an unsorted array.
%
%    Output, real A_OUT(N), the array has been reordered into a heap.
%
  a_out(1:n) = a(1:n);
%
%  Only nodes N/2 down to 1 can be "parent" nodes.
%
  for i = floor ( n/2 ) : -1 : 1
%
%  Copy the value out of the parent node.
%  Position IFREE is now "open".
%
    key = a_out(i);
    ifree = i;

    while ( 1 )
%
%  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
%  IFREE.  (One or both may not exist because they exceed N.)
%
      m = 2 * ifree;
%
%  Does the first position exist?
%
      if ( n < m )
        break;
      end
%
%  Does the second position exist?
%
      if ( m + 1 <= n )
%
%  If both positions exist, take the larger of the two values,
%  and update M if necessary.
%
        if ( a_out(m) < a_out(m+1) )
          m = m + 1;
        end

      end
%
%  If the large descendant is larger than KEY, move it up,
%  and update IFREE, the location of the free position, and
%  consider the descendants of THIS position.
%
      if ( a_out(m) <= key )
        break;
      end

      a_out(ifree) = a_out(m);
      ifree = m;

    end
%
%  Once there is no more shifting to do, KEY moves into the free spot IFREE.
%
    a_out(ifree) = key;

  end

  return
end


text:
function Gsv = d2Imis_dVdSg(Cg, V, lam, vcart)
%D2IMIS_DVDSG   Computes 2nd derivatives of current balance w.r.t. V and Sg.
%
%   The derivatives can be take with respect to polar or cartesian coordinates
%   of voltage, depending on the 4th argument.
%
%   GSV = D2IMIS_DVDSG(CG, V, LAM)
%   GSV = D2IMIS_DVDSG(CG, V, LAM, 0)
%
%   Returns a matrix containing the partial derivatives w.r.t. voltage angle
%   and magnitude of the product of a vector LAM with the 1st partial
%   derivatives of the real and reactive power generation.
%
%   GSV = D2IMIS_DVDSG(CG, V, LAM, 1)
%
%   Returns a matrix containing the partial derivatives w.r.t. real and
%   imaginary parts of voltage of the product of a vector LAM with the 1st
%   partial derivatives of the real and reactive power generation.
%
%   Takes the generator connection matrix, complex voltage vector V and
%   nb x 1 vector of multipliers LAM. Output matrices are sparse.
%
%   Examples:
%       Cg = sparse(gen(:, GEN_BUS), 1:ng, -, nb, ng);
%       Gsv = d2Imis_dVdSg(Cg, V, lam);
%
%       Here the output matrix corresponds to:
%           Gsv = [ Gpa Gpv;
%                   Gqa Gqv ];
%           Gpa = d/dVa (dImis_dPg.' * lam)
%           Gpv = d/dVm (dImis_dPg.' * lam)
%           Gqa = d/dVa (dImis_dQg.' * lam)
%           Gqv = d/dVm (dImis_dQg.' * lam)
%
%       [Grr, Gri, Gir, Gii] = d2Imis_dVdSg(Cg, V, lam, 1);
%
%       Here the output matrices correspond to:
%           Gsv = [ Gpr Gpi;
%                   Gqr Gqi ];
%           Gpr = d/dVr (dImis_dPg.' * lam)
%           Gpi = d/dVi (dImis_dPg.' * lam)
%           Gqr = d/dVr (dImis_dQg.' * lam)
%           Gqi = d/dVi (dImis_dQg.' * lam)
%
%   For more details on the derivations behind the derivative code used
%   in MATPOWER, see:
%
%   [TN2]  R. D. Zimmerman, "AC Power Flows, Generalized OPF Costs and
%          their Derivatives using Complex Matrix Notation", MATPOWER
%          Technical Note 2, February 2010. [Online]. Available:
%          https://matpower.org/docs/TN2-OPF-Derivatives.pdf
%          doi: 10.5281/zenodo.3237866
%   [TN3]  B. Sereeter and R. D. Zimmerman, "Addendum to AC Power Flows and
%          their Derivatives using Complex Matrix Notation: Nodal Current
%          Balance," MATPOWER Technical Note 3, April 2018. [Online].
%          Available: https://matpower.org/docs/TN3-More-OPF-Derivatives.pdf
%          doi: 10.5281/zenodo.3237900
%   [TN4]  B. Sereeter and R. D. Zimmerman, "AC Power Flows and their
%          Derivatives using Complex Matrix Notation and Cartesian
%          Coordinate Voltages," MATPOWER Technical Note 4, April 2018.
%          [Online]. Available: https://matpower.org/docs/TN4-OPF-Derivatives-Cartesian.pdf
%          doi: 10.5281/zenodo.3237909

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

%% default input args
if nargin < 4
    vcart = 0;      %% default to polar coordinates
end

nb = length(V);

if vcart
    D = Cg' * sparse(1:nb,1:nb, lam./conj(V.^2), nb, nb);

    G_Pg_Vr = D;
    G_Pg_Vi = -1j * D;
    G_Qg_Vr = G_Pg_Vi;
    G_Qg_Vi = -D;

    Gsv = [ G_Pg_Vr G_Pg_Vi;
            G_Qg_Vr G_Qg_Vi ];
else
    D = sparse(1:nb, 1:nb, 1./abs(V), nb, nb);
    E = sparse(1:nb, 1:nb, lam./conj(V), nb, nb);
    K = Cg' * E;
    L = K * D;

    G_Pg_Va = -1j * K;
    G_Pg_Vm = L;
    G_Qg_Va = -K;
    G_Qg_Vm = -1j * L;

    Gsv = [ G_Pg_Va G_Pg_Vm;
            G_Qg_Va G_Qg_Vm ];
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower