text:
function bd = mne_block_diag(A,n);
%
%   function bd = mne_block_diag(A,n)
%
%   Make or extract a sparse block diagonal matrix
%
%   If A is not sparse, then returns a sparse block diagonal "bd", diagonalized from the
%   elements in "A".
%   "A" is ma x na, comprising bdn=(na/"n") blocks of submatrices.
%   Each submatrix is ma x "n", and these submatrices are
%   placed down the diagonal of the matrix.
%
%   If A is already sparse, then the operation is reversed, yielding a block
%   row matrix, where each set of n columns corresponds to a block element
%   from the block diagonal.
%
%   Routine uses NO for-loops for speed considerations.

%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
%
% Copyright (c) 2005 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% Author: John C. Mosher 1993 - 2004
%
%
% Modifications for mne Matlab toolbox
%
%   Matti Hamalainen
%   2006
%   Revision 1.2  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%

if(~issparse(A)),        % then make block sparse
    [ma,na] = size(A);
    bdn = na/n;             % number of submatrices

    if(bdn - fix(bdn)),
        error('Width of matrix must be even multiple of n');
    end

    tmp = reshape([1:(ma*bdn)]',ma,bdn);
    i = zeros(ma*n,bdn);
    for iblock = 1:n,
        i((iblock-1)*ma+[1:ma],:) = tmp;
    end

    i = i(:);             % row indices foreach sparse bd

    j = [1:na];
    j = j(ones(ma,1),:);
    j = j(:);             % column indices foreach sparse bd

    bd = sparse(i,j,A(:));

else                 % already is sparse, unblock it

    [mA,na] = size(A);        % matrix always has na columns
    % how many entries in the first column?
    bdn = na/n;            % number of blocks
    ma = mA/bdn;            % rows in first block

    % blocks may themselves contain zero entries.  Build indexing as above
    tmp = reshape([1:(ma*bdn)]',ma,bdn);
    i = zeros(ma*n,bdn);
    for iblock = 1:n,
        i((iblock-1)*ma+[1:ma],:) = tmp;
    end

    i = i(:);             % row indices foreach sparse bd

    j = [0:mA:(mA*(na-1))];
    j = j(ones(ma,1),:);
    j = j(:);

    i = i + j;

    bd = full(A(i));     % column vector
    bd = reshape(bd,ma,na);    % full matrix
end

return

author: fieldtrip


text:
function P = rosa_lineextract(P, RADIUS, bUpdateConnectivity)
%
% we RECOVER CONNECTIVITY after each edge collapse, which is slow, but the
% result sampling is more uniform.
% bUpdateConnectivity: 0 for RECOVER CONNECTIVITY after each edge collapse, which is slow, but the
% result sampling is more uniform.
% bUpdateConnectivity: 1 for update CONNECTIVITY after each edge collapse, which is faster, and the
% result sampling is almost as same as that of the original method (bUpdateConnectivity=0).
%
% @author: Andrea Tagliasacchi
% reformed by jjcao
% @reform-data:     2009-9-3   using 1-ring to construct the connectivity
% graph
% @reform-data:     2010-8-20  decompose the original function into several
% small functions.
SHOW_RESULTS = true;

options.collapse_order = 1;
if (bUpdateConnectivity)
    [P.spls,P.corresp] = farthest_sampling_by_sphere(P.cpts, P.sample_radius);
    P.spls_adj = connect_by_inherit_neigh(P.cpts, P.spls, P.corresp, P.rings);
    [P.spls, P.spls_adj,P.corresp] = edge_collapse_update(P.cpts, P.spls, P.corresp, P.rings, P.spls_adj, options);
else
    [P.spls,P.corresp] = farthest_sampling_by_sphere(P.cpts, P.sample_radius);
    P.spls_adj = connect_by_inherit_neigh(P.cpts, P.spls, P.corresp, P.rings);
    [P.spls, P.spls_adj,P.corresp] = edge_collapse(P.cpts, P.spls, P.corresp, P.rings, P.spls_adj, options);
end

%%
if (SHOW_RESULTS)
    figure; movegui('northeast');set(gcf,'color','white');hold on;
    plot3( P.spls(:,1), P.spls(:,2), P.spls(:,3), '.r', 'markersize', 5);
    axis off; axis equal;set(gcf,'Renderer','OpenGL');
    plot_connectivity(P.spls, P.spls_adj);
    view3d zoom;
end

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function p = patchline(xs,ys,varargin)
% Plot lines as patches (efficiently)
%
% SYNTAX:
%     patchline(xs,ys)
%     patchline(xs,ys,zs,...)
%     patchline(xs,ys,zs,'PropertyName',propertyvalue,...)
%     p = patchline(...)
%
% PROPERTIES:
%     Accepts all parameter-values accepted by PATCH.
%
% DESCRIPTION:
%     p = patchline(xs,ys,zs,'PropertyName',propertyvalue,...)
%         Takes a vector of x-values (xs) and a same-sized
%         vector of y-values (ys). z-values (zs) are
%         supported, but optional; if specified, zs must
%         occupy the third input position. Takes all P-V
%         pairs supported by PATCH. Returns in p the handle
%         to the resulting patch object.
%
% NOTES:
%     Note that we are drawing 0-thickness patches here,
%     represented only by their edges. FACE PROPERTIES WILL
%     NOT NOTICEABLY AFFECT THESE OBJECTS! (Modify the
%     properties of the edges instead.)
%
%     LINUX (UNIX) USERS: One test-user found that this code
%     worked well on his Windows machine, but crashed his
%     Linux box. We traced the problem to an openGL issue;
%     the problem can be fixed by calling 'opengl software'
%     in your <http://www.mathworks.com/help/techdoc/ref/startup.html startup.m>.
%     (That command is valid at startup, but not at runtime,
%     on a unix machine.)
%
% EXAMPLES:
%%% Example 1:
%
% n = 10;
% xs = rand(n,1);
% ys = rand(n,1);
% zs = rand(n,1)*3;
% plot3(xs,ys,zs,'r.')
% xlabel('x');ylabel('y');zlabel('z');
% p  = patchline(xs,ys,zs,'linestyle','--','edgecolor','g',...
%     'linewidth',3,'edgealpha',0.2);
%
%%% Example 2: (Note "hold on" not necessary here!)
%
% t = 0:pi/64:4*pi;
% p(1) = patchline(t,sin(t),'edgecolor','b','linewidth',2,'edgealpha',0.5);
% p(2) = patchline(t,cos(t),'edgecolor','r','linewidth',2,'edgealpha',0.5);
% l = legend('sine(t)','cosine(t)');
% tmp = sort(findobj(l,'type','patch'));
% for ii = 1:numel(tmp)
%     set(tmp(ii),'facecolor',get(p(ii),'edgecolor'),'facealpha',get(p(ii),'edgealpha'),'edgecolor','none')
% end
%
%%% Example 3 (requires Image Processing Toolbox):
%%%   (NOTE that this is NOT the same as showing a transparent image on
%%%         of the existing image. (That functionality is
%%%         available using showMaskAsOverlay or imoverlay).
%%%         Instead, patchline plots transparent lines over
%%%         the image.)
%
% img = imread('rice.png');
% imshow(img)
% img = imtophat(img,strel('disk',15));
% grains = im2bw(img,graythresh(img));
% grains = bwareaopen(grains,10);
% edges = edge(grains,'canny');
% boundaries = bwboundaries(edges,'noholes');
% cmap = jet(numel(boundaries));
% ind = randperm(numel(boundaries));
% for ii = 1:numel(boundaries)
% patchline(boundaries{ii}(:,2),boundaries{ii}(:,1),...
%     'edgealpha',0.2,'edgecolor',cmap(ind(ii),:),'linewidth',3);
% end
%
% Written by Brett Shoelson, PhD
% brett.shoelson@mathworks.com
% 5/31/2012
%
% Revisions:
% 6/26 Improved rice.png example, modified FEX image.
%
% Copyright 2012 MathWorks, Inc.
%
% See also: patch, line, plot

[zs,PVs] = parseInputs(varargin{:});
if rem(numel(PVs),2) ~= 0
    % Odd number of inputs!
    error('patchline: Parameter-Values must be entered in valid pairs')
end

% Facecolor = 'k' is (essentially) ignored here, but syntactically necessary
if isempty(zs)
    p = patch([xs(:);NaN],[ys(:);NaN],'k');
else
    p = patch([xs(:);NaN],[ys(:);NaN],[zs(:);NaN],'k');
end

% Apply PV pairs
for ii = 1:2:numel(PVs)
    set(p,PVs{ii},PVs{ii+1})
end
if nargout == 0
    clear p
end

function [zs,PVs] = parseInputs(varargin)
if isnumeric(varargin{1})
    zs = varargin{1};
    PVs = varargin(2:end);
else
    PVs = varargin;
    zs = [];
end

author: yuxng


text:
function test_ft_topoplotER

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_topoplotER ft_topoplotTFR ft_topoplotIC

% this script tests the functionality of ft_topoplotER with respect to the
% different input datatypes. no other functionality is tested.
% the script has been written in order to test a clean up of the code

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/eeg/preproc_neuroscan16');
load(filename)

%there's an unresolved issue with duplicate labels 'FREE'
%FIXME
data.label{1} = 'FREE1';
data.label{2} = 'FREE2';
data.label{3} = 'FREE3';
data.label{4} = 'FREE4';

cfg = [];
cfg.preproc.demean = 'yes';
tlck = ft_timelockanalysis(cfg, data);

cfg = [];
cfg.layout = 'biosemi64.lay';
lay = ft_prepare_layout(cfg);

cfg        = [];
cfg.layout = lay;
cfg.interactive = 'yes';

%plot timelocked-data
figure; ft_topoplotER(cfg, tlck);drawnow

%reproduce bug 1239
cfgtmp = cfg;
cfgtmp.highlight = 'on';
cfgtmp.highlightchannel = {};
figure; ft_topoplotER(cfgtmp, tlck);drawnow

%plot subplots
cfg.xlim = (-1:0.25:1);
figure; ft_topoplotER(cfg, tlck);drawnow
cfg = rmfield(cfg, 'xlim');

%create component-data
cfgc = [];
cfgc.method = 'pca';
comp = ft_componentanalysis(cfgc, data);

%plot component-data
cfg.interactive = 'no';
cfg.component   = 1:10;
figure; ft_topoplotIC(cfg, comp);drawnow

cfg.component   = 5;
figure; ft_topoplotIC(cfg, comp);drawnow

%create frequency-data
cfgf = [];
cfgf.method = 'mtmfft';
cfgf.foilim = [0 100];
cfgf.taper  = 'hanning';
cfgf.output = 'pow';
freq = ft_freqanalysis(cfgf, data);

%plot frequency-data
cfg.interactive = 'yes';
%cfg = rmfield(cfg, 'component');
figure; ft_topoplotER(cfg, freq);drawnow

%create connectivity-data
cfgf.output = 'fourier';
freq2 = ft_freqanalysis(cfgf, data);

cfgc2 = [];
cfgc2.method = 'coh';
cfgc2.channel = freq2.label(5:end);
coh   = ft_connectivityanalysis(cfgc2, freq2);

%plot connectivity-data
cfg.refchannel = 'gui';
cfg.parameter = 'cohspctrm';
figure; ft_topoplotER(cfg, coh);drawnow % FIXME this causes a crash when a new reference is selected and the old one is not unselected
% FIXME it also crashes when more than one ref is selected

%create connectivity-data with sparse linear indexing
cfgc2.channelcmb = [repmat(freq2.label(5),[numel(freq2.label)-1 1]) freq2.label([1:4 6:end])];
coh2  = ft_connectivityanalysis(cfgc2, freq2);

%plot
cfg = [];
cfg.layout = 'biosemi64.lay';
cfg.parameter = 'cohspctrm';
% cfg.refchannel = 'gui';
cfg.refchannel = freq2.label{5};

coh2 = ft_checkdata(coh2, 'cmbstyle', 'full');
figure; ft_topoplotER(cfg, coh2);drawnow

coh2 = ft_checkdata(coh2, 'cmbstyle', 'sparse');
figure; ft_topoplotER(cfg, coh2);drawnow

%create connectivity-data with very sparse linear indexing
cfgc2.channelcmb = cfgc2.channelcmb(1:25,:);
coh3   = ft_connectivityanalysis(cfgc2, freq2);

%plot
figure; ft_topoplotER(cfg, coh3);drawnow

%create connectivity-data with even sparser linear indexing
cfgc2.channelcmb = [repmat(freq2.label(5),[10 1]) freq2.label(21:30);repmat(freq2.label(10),[10 1]) freq2.label(21:30)];
coh4 = ft_connectivityanalysis(cfgc2, freq2);

%plot
cfg.refchannel = freq2.label{5};

coh2 = ft_checkdata(coh2, 'cmbstyle', 'full');
figure; ft_topoplotER(cfg, coh4);drawnow

%plot: this works
coh2 = ft_checkdata(coh2, 'cmbstyle', 'sparse');
figure; ft_topoplotER(cfg, coh4);drawnow

%create connectivity-data with asymmetry
%the data are probably not full-rank creating a problem for the sf
%freq3 = ft_selectdata(freq2, 'channel', freq2.label(1:30));
%%subselection of channels does not help

cfgc2 = [];
cfgc2.method = 'granger';
granger = ft_connectivityanalysis(cfgc2, freq2);

cfg.refchannel = 'gui';
cfg.parameter = 'grangerspctrm';
figure; ft_topoplotER(cfg, granger);drawnow

%plot a stat structure, containing only 1 freq bin (so without freq-field)
stat = freq;
stat.stat = freq.powspctrm(:,10);
stat = rmfield(stat, 'freq');
if isfield(stat, 'cumtapcnt')
  stat = rmfield(stat, 'cumtapcnt');
end
stat = rmfield(stat, 'powspctrm');
stat.dimord = 'chan';

cfg = rmfield(cfg, 'refchannel');
cfg.parameter = 'stat';
cfg.interactive = 'no';
figure; ft_topoplotER(cfg, stat);drawnow

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ T, W ] = rd_lin_spline ( w0, t_array, n, c_array )

%*****************************************************************************80
%
%% RD_LIN_SPLINE discretizes and solves a 1D reaction diffusion problem.
%
%  Discussion:
%
%    The problem includes homogeneous Neumann boundary conditions at both ends.
%
%    The dynamics are given by the following reaction/diffusion equation
%    for the function W(T,X):
%
%      W_t = W_xx + NL(W,c),
%
%    where NL(W,c) is a polynomial in W:
%
%       NL(W,c) = c(1) + c(2) * W + c(3) * W^2 + c(4) * W^3
%
%    with Neumann boundary conditions at X = 0.0 and X = 1.0:
%
%      W_x(T,0.0) = 0.0
%      W_x(T,1.0) = 0.0
%
%    and initial condition at T = 0.0:
%
%      W(0,X) = sin ( pi * X ).
%
%    The problem is to be solved for 0.0 <= T <= 4.0, 0.0 <= X <= 1.0.
%
%    We use a finite element approximation with piecewise linear "hat"
%    functions. The resulting ODE model is:
%
%      M * w-dot(t) = - K * w(t) + NL(w, c)
%
%    where M is the mass matrix, K the stiffness matrix, and NL the nonlinear
%    term.
%
%      NL(w, c_array)=c_array(1)+ c_array(2)*w + c_array(3)*w^2 + c_array(4)*w^3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2011
%
%  Author:
%
%   Eugene Cliff
%
%  Reference:
%
%    Jeffrey Borggaard, John Burkardt, John Burns, Eugene Cliff,
%    Working Notes on a Reaction Diffusion Model: a Finite Element Formulation.
%
%  Parameters:
%
%    w0     - handle to initial function (on [0, 1])
%
%    t_array- array of times for output
%
%    n      - grid  parameter ( > 2)
%
%    c_array- coefficients in nonlinear term
%

%
%  Assemble the mass matrix.
%
    M = spdiags([ones(n+1,1) [2;4*ones(n-1,1);2] ones(n+1,1)], ...
                                                 [-1 0 1], n+1, n+1)/(6*n);
%
%  Assemble the stiff matrix.
%
    K = n*spdiags([-ones(n+1,1) [1;2*ones(n-1,1);1] -ones(n+1,1)], ...
                                                 [-1 0 1], n+1, n+1);
%
%  Set the initial condition.
%
   wr  = zeros(n+1,1);
   h_wr = @(x) w0(x).*basic_hat(n*x );
   wr(1) = quad(h_wr, 0, 1/n);                    % < w_0, \phi_1 >
   h_wr = @(x) w0(x).*basic_hat(n*x - n);
   wr(n+1) = quad(h_wr, (n-1)/n, 1);              % < w_0, \phi_n+1 >
   for jj=2:n
      h_wr = @(x) w0(x).*basic_hat(n*x + 1 - jj);
      wr(jj) = quad(h_wr, (jj-2)/n, jj/n);        % < w_0, \phi_i >
   end

   w_0 = M \ wr;

   ODE_opt = odeset('Mass', M);                  % constant mass matrix
   h_rhs = @(t, w) -K*w  + NL(w, c_array, n, M); % handle for rhs function

   [T, W] = ode15s(h_rhs, t_array, w_0, ODE_opt);% invoke ODE solver

  return
end
function val = NL (w, c, n, M )

%*****************************************************************************80
%
%% NL evaluates the nonlinear term in the reaction-diffusion equation.
%
      val = (c(1)/n)* [1; 2*ones(n-1,1); 1] ... % constant
          +  c(2) * M * w     ...               % linear
          +  c(3) * Nq(w, n)  ...               % quadratic
          +  c(4) * Nc(w, n) ;                  % cubic

  return
end
function val = Nq(w, n)

%*****************************************************************************80
%
%% NQ evaluates a quadratic nonlinear finite element function.
%
     w2 = w(:).^2; wx = (w(1:end-1,1)+ w(2:end,1)).^2;

     val = [2*w2(1)+wx(1); wx(1:end-1)+4*w2(2:end-1)+wx(2:end) ; ...
                                                 wx(end)+2*w2(end)]/(12*n);

  return
end
function val = Nc(w, n)

%*****************************************************************************80
%
%% NC evaluates a cubiic nonlinear finite element function.
%
     w2 = w(:).*w(:); w3 = w(:).*w2(:); wx = (w(1:end-1,1)+w(2:end,1)).^3;

     val = [3*w3(1)+wx(1)-w(1)*w(2)^2;
   wx(1:end-1)+6*w3(2:end-1)+wx(2:end)-w(2:end-1).*(w2(1:end-2)+w2(3:end));
                               wx(end)+3*w3(end)-w(end)*w(end-1)^2]/(20*n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Introduction to SURF (Speeded-Up Robust Features)
% SURF keypoint detection + keypoint drawing with OpenCV functions
%
% In this sample you will learn:
%
% * The basics of SURF.
% * How to use <matlab:doc('cv.SURF') cv.SURF> and its function
%   |cv.SURF.detect| to perform the detection process in order to find
%   interest points.
% * How to use the function <matlab:doc('cv.drawKeypoints') cv.drawKeypoints>
%   to draw the detected keypoints.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d7/d66/tutorial_feature_detection.html>
% * <https://docs.opencv.org/3.2.0/df/dd2/tutorial_py_surf_intro.html>
%

%% Theory
%
% We previously saw SIFT for keypoint detection and description. But it was
% comparatively slow and people needed more speeded-up version. In 2006, three
% people, Bay, H., Tuytelaars, T. and Van Gool, L, published another paper,
% "SURF: Speeded Up Robust Features" which introduced a new algorithm called
% SURF. As name suggests, it is a speeded-up version of SIFT.
%
% In SIFT, Lowe approximated Laplacian of Gaussian with Difference of Gaussian
% for finding scale-space. SURF goes a little further and approximates LoG
% with Box Filter. Below image shows a demonstration of such an approximation.
% One big advantage of this approximation is that, convolution with box filter
% can be easily calculated with the help of integral images. And it can be
% done in parallel for different scales. Also the SURF rely on determinant of
% Hessian matrix for both scale and location.
%
% <<https://docs.opencv.org/3.2.0/surf_boxfilter.jpg>>
%
% For orientation assignment, SURF uses wavelet responses in horizontal and
% vertical direction for a neighbourhood of size 6s. Adequate Gaussian weights
% are also applied to it. Then they are plotted in a space as given in below
% image. The dominant orientation is estimated by calculating the sum of all
% responses within a sliding orientation window of angle 60 degrees.
% Interesting thing is that, wavelet response can be found out using integral
% images very easily at any scale. For many applications, rotation invariance
% is not required, so no need of finding this orientation, which speeds up the
% process. SURF provides such a functionality called Upright-SURF or U-SURF.
% It improves speed and is robust upto $\pm 15^{\circ}$. OpenCV supports both,
% depending upon the flag, |Upright|. If it is 0, orientation is calculated.
% If it is 1, orientation is not calculated and it is faster.
%
% <<https://docs.opencv.org/3.2.0/surf_orientation.jpg>>
%
% For feature description, SURF uses Wavelet responses in horizontal and
% vertical direction (again, use of integral images makes things easier). A
% neighbourhood of size 20sX20s is taken around the keypoint where |s| is the
% size. It is divided into 4x4 subregions. For each subregion, horizontal and
% vertical wavelet responses are taken and a vector is formed like this,
% $v=( \sum{d_x}, \sum{d_y}, \sum{|d_x|}, \sum{|d_y|})$. This when represented
% as a vector gives SURF feature descriptor with total 64 dimensions. Lower
% the dimension, higher the speed of computation and matching, but provide
% better distinctiveness of features.
%
% For more distinctiveness, SURF feature descriptor has an extended 128
% dimension version. The sums of $d_x$ and $|d_x|$ are computed separately for
% $d_y < 0$ and $d_y \geq 0$. Similarly, the sums of $d_y$ and $|d_y|$ are
% split up according to the sign of $d_x$ , thereby doubling the number of
% features. It doesn't add much computation complexity. OpenCV supports both
% by setting the value of flag |Extended| with 0 and 1 for 64-dim and 128-dim
% respectively (default is 128-dim).
%
% Another important improvement is the use of sign of Laplacian (trace of
% Hessian Matrix) for underlying interest point. It adds no computation cost
% since it is already computed during detection. The sign of the Laplacian
% distinguishes bright blobs on dark backgrounds from the reverse situation.
% In the matching stage, we only compare features if they have the same type
% of contrast (as shown in image below). This minimal information allows for
% faster matching, without reducing the descriptor's performance.
%
% <<https://docs.opencv.org/3.2.0/surf_matching.jpg>>
%
% In short, SURF adds a lot of features to improve the speed in every step.
% Analysis shows it is 3 times faster than SIFT while performance is
% comparable to SIFT. SURF is good at handling images with blurring and
% rotation, but not good at handling viewpoint change and illumination change.
%
% OpenCV provides SURF functionalities just like SIFT. You initiate a SURF
% object with some optional conditions like 64/128-dim descriptors,
% Upright/Normal SURF, etc. All the details are well explained in docs. Then
% as we did in SIFT, we can use |cv.SURF.detect, |cv.SURF.compute|, etc. for
% finding keypoints and descriptors.
%

%% Code

%%
% Read image as grayscale
img = cv.imread(fullfile(mexopencv.root(),'test','butterfly.jpg'), ...
    'Grayscale',true);

%%
% Detect keypoints using SURF Detector
detector = cv.SURF();
detector.HessianThreshold = 400;
tic
keypoints = detector.detect(img);
toc

%%
whos keypoints
disp(keypoints(1))

%%
% Draw keypoints
out = cv.drawKeypoints(img, keypoints);
imshow(out);

%%
% We increase the Hessian Threshold to some large value.
% This is just to avoid drawing too many keypoints.
% In real applications, it is better to have a value between 300 and 500,
% as we may need all those keypoints when matching.
detector.HessianThreshold = 50000;
keypoints = detector.detect(img);
fprintf('%d keypoints\n', numel(keypoints));

%%
out = cv.drawKeypoints(img, keypoints, ...
    'Color',[255 0 0], 'DrawRichKeypoints',true);
imshow(out);

%%
% Now we apply U-SURF, so that it won't find the orientation.
detector.Upright = true;
keypoints = detector.detect(img);
out = cv.drawKeypoints(img, keypoints, ...
    'Color',[255 0 0], 'DrawRichKeypoints',true);
imshow(out);

%%
% Finally we change the descriptor size to 128.
disp(detector.descriptorSize)
detector.Extended = true;
disp(detector.descriptorSize)
descriptors = detector.compute(img, keypoints);
whos descriptors

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
  function [sino, hn, nn, Hk] = fbp2_sino_hilbert(sino, varargin)
%|function [sino, hn, nn, Hk] = fbp2_sino_hilbert(sino, [options])
%|
%| Apply band-limited Hilbert-transform filter to 2D sinogram.
%| Frequency response: H(u) = -1i * sign(u) * rect(u/2/umax)
%|
%| in
%|	sino	[nb (L)] sinogram(s)
%|
%| options
%|	dr | ds	(real)	sample spacing (in distance units, e.g., cm) (default 1)
%|	npad		# of padded samples. (default: 0, means next power of 2)
%|	decon1		deconvolve effect of linear interpolator? (default: 0)
%|	window	[npad]	samples of apodization window function
%|			for [-np/2,...,np/2-1]. (default: '' = plain ramp)
%|			or a string like 'hann' for some predefined windows
%| out
%|	sino	[nb (L)] filtered sinogram rows
%|	hn	[npad]	samples of band-limited Hilbert transform filter
%|	nn	[npad]	[-np/2,...,np/2-1] vector for convenience
%|	Hk	[npad]	spectral samples on [0 ... np-1]
%|
%| Copyright 2011-07-16, Jeff Fessler, University of Michigan

if nargin == 1 && streq(sino, 'test'), fbp2_sino_hilbert_test, clear, return, end
if nargin < 1, ir_usage, end

arg.dr = 1;
arg.window = '';
arg.npad = 0;
arg.decon1 = false;
arg = vararg_pair(arg, varargin, 'subs', {'ds', 'dr'});

dims = size(sino);
sino = reshape(sino, dims(1), []);
[sino hn nn Hk] = fbp2_sino_hilbert_do(sino, ...
	arg.dr, arg.window, arg.npad, arg.decon1);
sino = reshape(sino, [size(sino, 1) dims(2:end)]);

% fbp2_sino_hilbert_do()
function [sino, hn, nn, Hk] ...
	= fbp2_sino_hilbert_do(sino, dr, window, npad, decon1);

[nb na] = size(sino);
if ~npad
	npad = 2^ceil(log2(2*nb-1)); % padded size
end
sino = [sino; zeros(npad-nb,na)]; % padded sinogram

[hn nn] = fbp2_filter_hilbert_make(npad, dr);

Hk = 1i * imag_check(fft(fftshift(hn))); % trick: pure imaginary!

Hk = Hk .* fbp2_window(npad, window);

Hk = dr * Hk; % differential for discrete-space convolution vs integral

% linear interpolation is like blur with a triangular response,
% so we can compensate for this approximately in frequency domain
if decon1
	Hk = Hk ./ fftshift(nufft_sinc(nn / npad).^2);
end

sino = ifft_sym( fft(sino, [], 1) .* repmat(Hk, [1 na]), [], 1); % apply filter

% fbp2_filter_hilbert_make()
function [hn, nn] = fbp2_filter_hilbert_make(n, dr)
nn = [-(n/2):(n/2-1)]';
u0 = 1/2/dr;
hn = 2 * u0 * nufft_sinc(nn/2) .* sin(pi/2 * nn);

% imag_check()
% take imaginary part but check that real part is negligible
function out = imag_check(in, varargin)
out = reale(-1i * in, varargin{:});

% fbp2_sino_hilbert_test()
function fbp2_sino_hilbert_test
nb = 2^5;
dr = 0.1;
[sino1 h1 nn H1] = fbp2_sino_hilbert(zeros(nb,2), 'dr', dr);
h2 = 1 ./ (pi * nn * dr); % samples of ideal non-bandlimited

if im
	clf, subplot(121)
	plot(nn, h1, '.-', nn, h2, '-')
	xlabel 'n', ylabel 'h[n]'
	legend('band-limited', 'ideal')

	subplot(122)
	plot(nn, imag_check(fftshift(H1)), '.-', nn, -sign(nn), '-')
	xlabel 'k', ylabel 'imag(H[k])'
	legend('band-limited', 'ideal', 'location', 'north')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function c = integrate(c)

%INTEGRATE   integrates each trace.
%
% C = INTEGRATE(C) integrates each trace. Prior to integrating, the trend
% and mean of each trace is removed. Depending on the application, it is
% recommended that users also consider applying a gentle high pass filter
% to remove unwanted low frequencies that can dominate the integrated
% waveforms.

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

c.W = detrend(c.W);
c.W = demean(c.W);
c.W = integrate(c.W);

author: geoscience-community-codes

repo: GISMO


text:
function rez = runTemplates(rez)

% this function will run Kilosort2 initializing at some previously found
% templates. These must be specified in rez.W, rez.U and rez.mu. The batch
% number at which the new template run is started is by default 1 (modify it by changing rez.istart).

% If you don't mind changes in template throughout the recording, but want
% to use the results from a previous session after splitting and merging, just pass
% the rez output from this previous session.

% Keep in mind that these templates will nonetheless
% change substantially as they track neurons throughout the new recording.
% As such, it is best to use the templates obtained at the very end of the
% previous recording, which are found in rez.WA(:,:,:,end),
% rez.UA(:,:,:,end), rez.muA(:,end). Keep in mind that any merging and
% splitting steps are done after running the templates, and would result in
% difference between sessions if done separately on each session.

% update: these recommendations no longer apply for the datashift version!

if sum(isfield(rez, {'istart'}))<1
    warning('if using pre-loaded templates, please specify istart, defaulting to 1');
    rez.istart = 1;
end

if sum(isfield(rez, {'W', 'U', 'mu'}))<3
    error('missing at least one field: W, U, mu');
end

Nbatches = rez.ops.Nbatch;

iorder = 1:Nbatches;

[rez, st3, fW,fWpc] = trackAndSort(rez, iorder);

% sort all spikes by batch -- to keep similar batches together,
% which avoids false splits in splitAllClusters. Break ties
% [~, isort] = sortrows(st3,[5,1,2,3,4]);
% st3 = st3(isort, :);
% fW = fW(:, isort);
% fWpc = fWpc(:, :, isort);

% just display the total number of spikes
fprintf( 'Number of spikes before applying cutoff: %d\n', size(st3,1));

rez.st3 = st3;
rez.st2 = st3; % keep also an st2 copy, because st3 will be over-written by one of the post-processing steps

% the template features are stored in cProj, like in Kilosort1
rez.cProj    = fW';

%  permute the PC projections in the right order
rez.cProjPC     = permute(fWpc, [3 2 1]); %zeros(size(st3,1), 3, nNeighPC, 'single');
% iNeighPC keeps the indices of the channels corresponding to the PC features

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [ X ] = hyperFcls( M, U )
%HYPERFCLS Performs fully constrained least squares on pixels of M.
% hyperFcls performs fully constrained least squares of each pixel in M
% using the endmember signatures of U.  Fully constrained least squares
% is least squares with the abundance sum-to-one constraint (ASC) and the
% abundance nonnegative constraint (ANC).
%
% Usage
%   [ X ] = hyperFcls( M, U )
% Inputs
%   M - HSI data matrix (p x N)
%   U - Matrix of endmembers (p x q)
% Outputs
%   X - Abundance maps (q x N)
%
% References
%   "Fully Constrained Least-Squares Based Linear Unmixing." Daniel Heinz,
% Chein-I Chang, and Mark L.G. Althouse. IEEE. 1999.
%
% Written by
%    Matlab Hyperspectral Toolbox - Toolbox of advanced algorithms for
%    hyperspectral processing and exploitation.
%    http://sourceforge.net/apps/mediawiki/matlabhyperspec/index.php?title=Main_Page.
%    Toolbox Available at http://sourceforge.net/projects/matlabhyperspec/.

if (ndims(U) ~= 2)
    error('M must be a p x q matrix.');
end

[p1, N] = size(M);
[p2, q] = size(U);
if (p1 ~= p2)
    error('M and U must have the same number of spectral bands.');
end

p = p1;
X = zeros(q, N);
Mbckp = U;
for n1 = 1:N
    count = q;
    done = 0;
    ref = 1:q;
    r = M(:, n1);
    U = Mbckp;
    while not(done)
        als_hat = inv(U.'*U)*U.'*r;
        s = inv(U.'*U)*ones(count, 1);

        % IEEE Magazine method (http://www.planetary.brown.edu/pdfs/3096.pdf)
        % Contains correction to sign.  Error in original paper.
        afcls_hat = als_hat - inv(U.'*U)*ones(count, 1)*inv(ones(1, count)*inv(U.'*U)*ones(count, 1))*(ones(1, count)*als_hat-1);

        % See if all components are positive.  If so, then stop.
        if (sum(afcls_hat>0) == count)
            alpha = zeros(q, 1);
            alpha(ref) = afcls_hat;
            break;
        end
        % Multiply negative elements by their counterpart in the s vector.
        % Find largest abs(a_ij, s_ij) and remove entry from alpha.
        idx = find(afcls_hat<0);
        afcls_hat(idx) = afcls_hat(idx) ./ s(idx);
        [val, maxIdx] = max(abs(afcls_hat(idx)));
        maxIdx = idx(maxIdx);
        alpha(maxIdx) = 0;
        keep = setdiff(1:size(U, 2), maxIdx);
        U = U(:, keep);
        count = count - 1;
        ref = ref(keep);
    end
    X(:, n1) = alpha;
end

return;

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [ Acc ] = ssvep_performance(CNT , params)
opt = opt_cellToStruct(params);
%% CCA - Anaysis
SMT=[];
for onoff=1:2
    cnt = CNT{onoff};
    cnt=prep_filter(cnt, {'frequency', opt.band});
    CNTch = prep_selectChannels(cnt, {'Index', opt.channel_index});
    SMT_iter = prep_segmentation(CNTch, {'interval', opt.time_interval});
    if onoff==1
        SMT= SMT_iter;
        clear SMT_iter
    else
        SMT = prep_addTrials(SMT, SMT_iter);
        clear SMT_iter
    end
end

tot = size(SMT.x, 2);
count1= tot;
for i = 1: size(SMT.x, 2)
    res_cca = ssvep_cca_analysis(squeeze(SMT.x(:,i,:)),{'marker',opt.marker;'freq', opt.freq;'fs', opt.fs;'time',opt.time});
    [~, ind] = max(res_cca);
    if SMT.y_dec(i) ~= ind
        count1 = count1 -1;
    end
end
Acc =count1/tot;
clear CNTch SMT tot count1 i res_cca in
end

author: PatternRecognition

repo: OpenBMI


text:
function a_n = extractAxisCoeffs(a_nm)
%EXTRACTAXISCOEFFS Exctract the N+1 coeffs of m=0 from the full (order+1)^2
% coefficient vector
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% EXTRACTAXISCOEFFS.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = sqrt(size(a_nm,1))-1;
nSets = size(a_nm,2);
a_n = zeros(N+1,nSets);
for ns=1:nSets
    for n=0:N
        a_n(n+1,ns) = a_nm(n^2 + n + 1);
    end
end

end

author: polarch

repo: Spherical-Array-Processing


text:
function lo=lpcrf2lo(rf)
%LPCRF2LO Convert reflection coefficients to log area ratios LO=(RF)
%the output values are limited to about +-14.5

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcrf2lo.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r=max(min(rf,1-1E-6),1E-6-1);
lo=log((1-r)./(1+r));

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function F=makeRFSfilters
% Returns the RFS filter bank of size 49x49x38 in F. The MR8, MR4 and
% MRS4 sets are all derived from this filter bank. To convolve an
% image I with the filter bank you can either use the matlab function
% conv2, i.e. responses(:,:,i)=conv2(I,F(:,:,i),'valid'), or use the
% Fourier transform.

  SUP=49;                 % Support of the largest filter (must be odd)
  SCALEX=[1,2,4];         % Sigma_{x} for the oriented filters
  NORIENT=6;              % Number of orientations

  NROTINV=2;
  NBAR=length(SCALEX)*NORIENT;
  NEDGE=length(SCALEX)*NORIENT;
  NF=NBAR+NEDGE+NROTINV;
  F=zeros(SUP,SUP,NF);
  hsup=(SUP-1)/2;
  [x,y]=meshgrid([-hsup:hsup],[hsup:-1:-hsup]);
  orgpts=[x(:) y(:)]';

  count=1;
  for scale=1:length(SCALEX),
    for orient=0:NORIENT-1,
      angle=pi*orient/NORIENT;  % Not 2pi as filters have symmetry
      c=cos(angle);s=sin(angle);
      rotpts=[c -s;s c]*orgpts;
      F(:,:,count)=makefilter(SCALEX(scale),0,1,rotpts,SUP);
      F(:,:,count+NEDGE)=makefilter(SCALEX(scale),0,2,rotpts,SUP);
      count=count+1;
    end;
  end;
  F(:,:,NBAR+NEDGE+1)=normalise(fspecial('gaussian',SUP,10));
  F(:,:,NBAR+NEDGE+2)=normalise(fspecial('log',SUP,10));
return

function f=makefilter(scale,phasex,phasey,pts,sup)
  gx=gauss1d(3*scale,0,pts(1,:),phasex);
  gy=gauss1d(scale,0,pts(2,:),phasey);
  f=normalise(reshape(gx.*gy,sup,sup));
return

function g=gauss1d(sigma,mean,x,ord)
% Function to compute gaussian derivatives of order 0 <= ord < 3
% evaluated at x.

  x=x-mean;num=x.*x;
  variance=sigma^2;
  denom=2*variance;
  g=exp(-num/denom)/sqrt(pi*denom);
  switch ord,
    case 1, g=-g.*(x/variance);
    case 2, g=g.*((num-variance)/(variance^2));
  end;
return

function f=normalise(f), f=f-mean(f(:)); f=f/sum(abs(f(:))); return

author: Cloud-CV

repo: object-proposals


text:
function rs = rang(s)

%tstoolbox/@signal/rang
%   Syntax:
%     * rs = rang(s)
%
%   Transform scalar time series to rang values.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

c = rang(s.core);
rs = signal(c, s);				% special constructor calling syntax for working routines

rs = setyunit(rs, unit);
rs = setyname(rs, 'Rang');
rs = addhistory(rs,  ['Transformed to rang values'] );
rs = addcommandlines(rs, 's = rang(s');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function outcome = Convertion(Average,G,g,last_gen,change_threshold)
% Conversion condition for enter next stage

%------------------------------- 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".
%--------------------------------------------------------------------------

    if (G-g > last_gen) && (max(abs(Average(G,:)-Average(G-last_gen,:)))<=change_threshold)
        outcome = true;
    else
        outcome = false;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function F = stream2d(A,vector)
%% Stream function for any D2Qn Lattice
% There are only 4 basic movements and any other is a combination of the
% basic ones:

% Matlab equivalent solution:
% F = circshift(A,[y,x])

%% Evaluate input data

v = length(vector);

    y = vector(1);
    x = vector(2);

[m,n] = size(A);

    b = zeros(1,n);
    c = zeros(m,1);

%% Stream in XY plane

if x == 1       % x+ stream
    b = A(:,n);
    for i = n:-1:2;
        A(:,i) = A(:,i-1);
    end
    A(:,1) = b;
elseif x == -1  % x- stream
    b = A(:,1);
    for i = 1:n-1;
        A(:,i) = A(:,i+1);
    end
    A(:,n) = b;
else            % x=0 stream
    %do nothing
end

if y == 1       % y+ stream
    c = A(m,:);
    for j = m:-1:2;
        A(j,:) = A(j-1,:);
    end
    A(1,:) = c;
elseif y == -1  % y- stream
    c = A(1,:);
    for j = 1:m-1;
        A(j,:) = A(j+1,:);
    end
    A(m,:) = c;
else            % y=0 stream
    %do nothing
end

F = A;

repo: Aero-matlab


text:
function obj = coef_Jacob2_qt_pnp_func_new(in1)
coef_J2 = in1(:,2);
coef_J5 = in1(:,5);
coef_J6 = in1(:,6);
coef_J7 = in1(:,7);
coef_J20 = in1(:,20);
coef_J21 = in1(:,21);
coef_J22 = in1(:,22);
coef_J40 = in1(:,40);
coef_J14 = in1(:,14);
coef_J41 = in1(:,41);
coef_J50 = in1(:,50);
coef_J15 = in1(:,15);
coef_J33 = in1(:,33);
coef_J42 = in1(:,42);
coef_J51 = in1(:,51);
coef_J16 = in1(:,16);
coef_J34 = in1(:,34);
coef_J43 = in1(:,43);
coef_J17 = in1(:,17);
coef_J35 = in1(:,35);
coef_J44 = in1(:,44);
coef_J18 = in1(:,18);
coef_J36 = in1(:,36);
coef_J45 = in1(:,45);
coef_J19 = in1(:,19);
coef_J37 = in1(:,37);
coef_J46 = in1(:,46);
coef_J38 = in1(:,38);
coef_J47 = in1(:,47);
coef_J39 = in1(:,39);
coef_J48 = in1(:,48);
coef_J49 = in1(:,49);
obj = [-coef_J2+coef_J14.*3.0,coef_J15.*2.0,coef_J16.*2.0,-coef_J2+coef_J17,coef_J18,-coef_J2+coef_J19,coef_J20,coef_J21,coef_J22,coef_J2,coef_J5.*-2.0+coef_J33.*4.0,-coef_J6+coef_J34.*3.0,-coef_J7+coef_J35.*3.0,coef_J5.*-2.0+coef_J36.*2.0,coef_J37.*2.0,coef_J5.*-2.0+coef_J38.*2.0,coef_J39.*2.0,coef_J40.*2.0,coef_J41.*2.0,coef_J5.*2.0,-coef_J6+coef_J42,-coef_J7+coef_J43,-coef_J6+coef_J44,coef_J45,coef_J46,coef_J47,coef_J6,-coef_J7+coef_J48,coef_J49,coef_J50,coef_J51,coef_J7];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% [INPUT]
% r = A vector of floats (-Inf,Inf) of length t representing the logarithmic returns.
% v = A vector of floats [0,Inf) of length t representing the trading volumes.
% sv = A float t-by-k matrix (-Inf,Inf) representing the state variables.
% bw = An integer [21,252] representing the dimension of each rolling window.
% mem = A string representing the MEM type:
%   - 'B' for Baseline MEM;
%   - 'A' for Asymmetric MEM;
%   - 'P' for Asymmetric Power MEM;
%   - 'S' for Spline MEM.
% mag = An integer [1,Inf) obtained as 10^x representing the magnitude of logarithmic returns and trading volumes (optional, default=[]).
%
% [OUTPUT]
% illiq = A column vector of floats [0,Inf) of length t representing the ILLIQ indicator.
% illiqc = A column vector of floats [0,Inf) of length t representing the ILLIQ indicator with covariates if state variables are provided, an empty array otherwise.
% knots = A row vector of floats [0,Inf) containing the optimal number of knots if Spline MEM is used, an empty array otherwise.

function [illiq,illiqc,knots] = illiq_indicator(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('r',@(x)validateattributes(x,{'double'},{'real' 'finite' 'vector' 'nonempty'}));
        ip.addRequired('v',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('sv',@(x)validateattributes(x,{'double'},{'real' 'finite'}));
        ip.addRequired('bw',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 21 '<=' 252 'scalar'}));
        ip.addRequired('mem',@(x)any(validatestring(x,{'A' 'B' 'P' 'S'})));
        ip.addOptional('mag',[],@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [r,v,sv,mag] = validate_input(ipr.r,ipr.v,ipr.sv,ipr.mag);
    bw = ipr.bw;
    mem = ipr.mem;

    nargoutchk(2,3);

    [illiq,illiqc,knots] = illiq_indicator_internal(r,v,sv,bw,mem,mag);

end

function [illiq,illiqc,knots] = illiq_indicator_internal(r,v,sv,bw,mem,mag)

    alpha = 2 / (bw + 1);

    input = mag .* (abs(r) ./ v);
    input(~isfinite(input) | (input == 0)) = NaN;
    input(isnan(input)) = mean(input,'omitnan');

    if (any(strcmp(mem,{'A' 'P'})))
        input = [input r];
    end

    knots = [];

    [illiq,~,mem_params] = multiplicative_error(input,mem);
    illiq = [illiq(1); filter(alpha,[1 (alpha - 1)],illiq(2:end),(1 - alpha) * illiq(1))];
    illiq = (illiq - min(illiq)) ./ (max(illiq) - min(illiq));

    if (strcmp(mem,'S'))
        knots(1) = mem_params(1);
    end

    if (isempty(sv))
        illiqc = [];
    else
        [illiqc,~,mem_params] = multiplicative_error([input sv],mem);
        illiqc = [illiqc(1); filter(alpha,[1 (alpha - 1)],illiqc(2:end),(1 - alpha) * illiqc(1))];
        illiqc = (illiqc - min(illiqc)) ./ (max(illiqc) - min(illiqc));

        if (strcmp(mem,'S'))
            knots(2) = mem_params(1);
        end
    end

end

function [r,v,sv,mag] = validate_input(r,v,sv,mag)

    data = {r(:) v(:)};

    l = unique(cellfun(@numel,data));

    if (numel(l) ~= 1)
        error('The number of elements of ''r'' and ''v'' must be equal.');
    end

    if (l < 5)
        error('The value of ''r'' and ''v'' is invalid. Expected inputs to be vectors containing at least 5 elements.');
    end

    [r,v] = deal(data{:});

    if (~isempty(sv))
        if (size(sv,1) ~= l)
            error(['The value of ''sv'' is invalid. Expected input to be a matrix with ' num2str(l) ' rows.']);
        end
    end

    if (isempty(mag))
        mag_r = floor(round((log(abs(r)) ./ log(10)),15));
        mag_r(~isfinite(mag_r)) = [];
        mag_r = round(abs(mean(mag_r)),0);

        mag_v = floor(round((log(abs(v)) ./ log(10)),15));
        mag_v(~isfinite(mag_v)) = [];
        mag_v = round(mean(mag_v),0);

        mag = 10^(mag_r + mag_v);
    else
        if ((mag ~= 1) && (rem(mag,10) ~= 0))
            error('The value of ''mag'' is invalid. Expected input to be an integer obtained as 10^x.');
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [x,r,g,info] = spg_lasso(A,b,tau,varargin )
%SPG_LASSO  Solve the LASSO problem
%
%   SPG_LASSO is designed to solve the LASSO problem
%
%   (LASSO)  minimize  ||AX - B||_2  subject to  ||X||_1 <= tau,
%
%   where A is an M-by-N matrix, B is an M-vector, and TAU is a
%   nonnegative scalar.  In all cases below, A can be an explicit M-by-N
%   matrix or matrix-like object for which the operations  A*x  and  A'*y
%   are defined (i.e., matrix-vector multiplication with A and its
%   adjoint.)
%
%   Also, A can be a function handle that points to a function with the
%   signature
%
%   v = A(w,mode)   which returns  v = A *w  if mode == 1;
%                                  v = A'*w  if mode == 2.
%
%   X = SPG_LASSO(A,B,TAU) solves the LASSO problem.
%
%   X = SPG_LASSO(A,B,TAU,OPTIONS) specifies options that are set using
%   SPGSETPARMS.
%
%   [X,R,G,INFO] = SPG_LASSO(A,B,TAU,OPTIONS) additionally returns the
%   residual R = B - A*X, the objective gradient G = A'*R, and an INFO
%   structure.  (See SPGL1 for a description of this last output argument.)
%
%   See also spgl1, spgSetParms, spg_bp, spg_bpdn.

%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
%   http://www.cs.ubc.ca/labs/scl/spgl1
%   $Id: spg_lasso.m 1074 2008-08-19 05:24:28Z ewout78 $

if ~exist('tau','var'), tau = []; end
if ~exist('b','var') || isempty(b)
    error('Second argument cannot be empty.');
end
if ~exist('A','var') || isempty(A)
    error('First argument cannot be empty.');
end

sigma = [];
x0  = [];
[x,r,g,info] = spgl1(A,b,tau,sigma,x0,varargin{:});

author: mpf


text:
function [ I ] = kernelcmi( x, y, z, h, ind )
% Kernel-based estimate for conditional mutual information I(X, Y|Z)
% h - kernel width; ind - subset of data on which to estimate MI

[Nx, Mx]=size(x);
[Ny, My]=size(y);
[Nz, Mz]=size(z);

if any([Nx Ny Nz My Mz] ~= [1 1 1 Mx Mx])
    error('Bad sizes of arguments');
end

if nargin < 4
    % Yields unbiased estiamte when Mx->inf
    % and low MSE for two joint gaussian variables
    alpha = 0.25;
    h = (Mx + 1) / sqrt(12) / Mx ^ (1 + alpha);
end

if nargin < 5
    ind = 1:Mx;
end

% Copula-transform variables
x = ctransform(x);
y = ctransform(y);
z = ctransform(z);

h2 = 2*h^2;

% Pointwise values for kernels
Kx = squareform(exp(-ssqd([x;x])/h2))+eye(Mx);
Ky = squareform(exp(-ssqd([y;y])/h2))+eye(Mx);
Kz = squareform(exp(-ssqd([z;z])/h2))+eye(Mx);

% Kernel sums for marginal probabilities
Cx = sum(Kx);
Cy = sum(Ky);

% Kernel products for joint probabilities
Kxz = Kx.*Kz;
Kyz = Ky.*Kz;
Kxyz = Kx.*Ky.*Kz;

f = ((Cx.*Cy)*Kz).*sum(Kxyz)./(Cx*Kyz)./(Cy*Kxz);
I = mean(log(f(ind)));

end

author: Sable

repo: mcbench-benchmarks


text:
function[]=makefigs_patchcontourf
%MAKEFIGS_PATCHCONTOURF  Makes a sample figure for PATCHCONTOURF.

load jtopo
use jtopo

figure
h=patchcontourf(lon,lat,topo,0,'k');latratio(30)
axis([-180 180 -75 75]),boxon,axis tight
title('Continents filled with PATCHCONTOURF')

%Make a second figure if m_map is installed
if exist('m_map')==7
    figure
    m_proj('miller','lat',75);
    m_grid('linestyle','none')
    patchcontourf(lon,lat,topo,0,'k','m_map');
    title('Continents filled with PATCHCONTOURF, M-MAP option')
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function addAllenCtxOutlines(bregma, lambda, lineColor, pixSize)
% to get bregma/lambda, do:
% >> figure; imagesc(myImg); axis image;
% then use the cursor to click on bregma, then swap the coords (if cursor
% says X=xx and Y=yy, then bregma is [yy, xx]).
% or leave blank to get prompted with ginput.
% lambda just needs to be on the midline posterior to bregma.
if isempty(bregma)
    fprintf(1, 'click bregma\n');
    [x,y] = ginput(1);
    bregma = [y,x]
    fprintf(1, 'click lambda\n');
    [x,y] = ginput(1);
    lambda = [y,x]
end

apDir = lambda-bregma; apDir = apDir./norm(apDir);

if nargin<4
    pixSize = 0.0217; % mm/pix. This is for PCO edge 5.5 with 0.6x mag (as kilotrode)
end

ccfbregma = allenCCFbregma()/100/pixSize;

load(fullfile(fileparts(mfilename('fullpath')), 'ctxOutlines.mat'));

hold on;
for q = 1:numel(coords) % coords is from ctxOutlines.mat

    % these are in 10um voxel coordinates, so first convert to mm, then to
    % pixels
    cx = coords(q).x/100/pixSize;
    cy = coords(q).y/100/pixSize;

    % to do this transformation, first subtract bregma to zero, then
    % rotate, then add back the other bregma
    cx = cx-ccfbregma(3); cy = cy-ccfbregma(1);

    T = affineMat.rot(-atan(apDir(2)/apDir(1)));

    newc = T*[cx cy ones(size(cx))]';
    cx = newc(1,:)'; cy = newc(2,:)';

    cx = cx+bregma(2); cy = cy+bregma(1);

    plot(cx,cy, 'LineWidth', 1.0, 'Color', lineColor);
    hold on;
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function bb2()
% BB2 Bouncing Smiley Ball Physics
% BB2 is a graphical User Interface that creates a simple simulation of a ball bouncing
% off the ground.
%
% Ball Parameters are:
%
% Height from the center of the ball to the ground is 10. This is constant and cannot
% change. Instead, change the other parameters since they're all relative.
%
% Radius: Radius of the ball. It's best to put it between 1 and 9. Any value equal to
% or above 10 will present problems as 10 is the distance from the centre of the ball
% to the ground.
%
% Gravity: Affect the ratio in which the ball gains speed in the negative Y direction.
%
% Initial Yv: Initial velocity in the Y direction. Can be positive, negative, or anything.
%
% Initial Xv: Initial velocity in the X direction. Can be positive, negative, or anything.
%
% Vertical Speed Conservation (%): Affects how much energy (velocity) is retained in the
% Y direction after hitting the ground. A value of 0 to 100 is normal. A value of 0 means
% no energy is retained and the ball comes to a halt after hitting the ground. A value
% of 100 means all energy is retained and the ball will bounce indefinitely. A value of
% over 100 will cause energy to increase after hitting the ground.
%
% Horizontal Speed Conservation (%): Same as Vertical Speed Conservation but in the
% X direction.
%
% Initial Angular Velocity (CW): The initial rotational velocity of the ball (spin).
%
% Plotting a path enables you to see the x-y diagram of the ball's movement and gives
% you the ability to pan the axes with the mouse.
%
%
% Please note that the ball may stop moving at different locations for different
% trials, even with the same input parameters. This is due to the slight
% randomization of the deformations of the ball and the ground, which could lead
% to the ball getting "stuck" to the ground at times. This is normal and
% deliberate.
%
% Set the vertical speed conservation as a value larger than 100 (e.g. 170) for
% a small bonus ;)
%
% I'd like to thank John D'Errico and Kennith Eaton whose ideas helped me a lot.
%
%
% To launch the GUI, type bb2 in the command window with this file in the current
% directory. Alternatively, you can choose Debug -> Run from this editor window, or
% press F5.
%
% Tested on MATLAB 7.4 (R2007a)
%
% Husam Aldahiyat, 2009
% numandina@gmail.com
%

figure('units','normalized','position',[.2 .2 .65 .65],'menubar','none','numbertitle','off','color','w','name','Bouncing Smiley')
axes('position',[.25 .05 .65 .75])
ed1=uicontrol('style','edit','units','normalized','position',[.025 .895 .1 .05],'backgroundcolor','w','string','1','callback',@init);
uicontrol('style','text','units','normalized','position',[.025 .96 .1 .025],'backgroundcolor','w','string','Radius');
ed2=uicontrol('style','edit','units','normalized','position',[.025 .775 .1 .05],'backgroundcolor','w','string','9.81');
uicontrol('style','text','units','normalized','position',[.025 .84 .1 .025],'backgroundcolor','w','string','Gravity');
ed3=uicontrol('style','edit','units','normalized','position',[.025 .65 .1 .05],'backgroundcolor','w','string','0');
uicontrol('style','text','units','normalized','position',[.025 .72 .1 .025],'backgroundcolor','w','string','Initial Yv');
ed4=uicontrol('style','edit','units','normalized','position',[.025 .535 .1 .05],'backgroundcolor','w','string','500');
uicontrol('style','text','units','normalized','position',[.025 .6 .1 .025],'backgroundcolor','w','string','Initial Xv');
ed5=uicontrol('style','edit','units','normalized','position',[.025 .415 .1 .05],'backgroundcolor','w','string','70');
uicontrol('style','text','units','normalized','position',[.005 .48 .21 .025],'backgroundcolor','w','string',...
	'Vertical Speed Conservation (%)','horizontalalignment','left');
ed6=uicontrol('style','edit','units','normalized','position',[.025 .295 .1 .05],'backgroundcolor','w','string','50');
uicontrol('style','text','units','normalized','position',[.005 .36 .21 .025],'backgroundcolor','w','string',...
	'Horizontal Speed Conservation (%)','horizontalalignment','left');
ed7=uicontrol('style','edit','units','normalized','position',[.025 .175 .1 .05],'backgroundcolor','w','string','25');
uicontrol('style','text','units','normalized','position',[.005 .24 .21 .025],'backgroundcolor','w','string',...
	'Initial Angular Velocity (CW)','horizontalalignment','left');

tb=uicontrol('style','togglebutton','string','Start','callback',@go,'units','normalized','position',[.025 .05 .1 .05],...
	'backgroundcolor','g');
chk=uicontrol('style','checkbox','string','Plot Path','units','normalized','position',[.025 .125 .1 .025],'backgroundcolor','w');
dummy1=uicontrol('style','text','max',2,'visible','off');
init;

	function [hb,h2,hx,h1,r,t,KE,hL]=init(varargin)
		r=str2double(get(ed1,'string'));
		t=linspace(0,2*pi,100);
		cla reset
		hold on
		h1=fill((5)+r*cos(t),(10)+r*sin(t),'y');
		am =[
			0    0    0    0    0    1    1    0    0    0    0    0    0    0    0    0    0    1    1    0    0    0    0    0    0
			0    0    0    0    1    1    1    1    0    0    0    0    0    0    0    0    1    1    1    1    0    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    0    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
			0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
			0    0    0    0    1    1    1    1    0    0    0    0    0    0    0    0    1    1    1    1    0    0    0    0    0
			0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
			1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1
			1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1
			1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1
			1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1
			1    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1
			1    1    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1    1
			0    0    1    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1    0    0
			0    0    0    0    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1    0    0    0
			0    0    0    0    0    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    0    0    0    0    0
			0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    1    1    1    1    0    0    0    0    0    0
			];
		am=[zeros(2,25);am;zeros(2,25)];
		colormap([1 1 0;0 0 0]);
		axis([0 10 -1 10+r+1])
		axis equal
		hx=axis;
		hL=line(1000.*[hx(1),hx(2)],[10 10]);
		axis(hx)
		set(hL,'color','c','linestyle','--')
		set(h1,'linewidth',3)
		hL=line(1000.*[hx(1) hx(2)],[0 0]);
		set(hL,'color','k','linestyle','-','linewidth',3)
		h2=plot(5,10,'k+');
		hb=[h1,h2];
		KE=surface(flipud(am),'facecolor','texturemap','xdata',linspace(get(h2,'xdata')-r/2,get(h2,'xdata')+r/2,...
				25),'ydata',linspace(get(h2,'ydata')-r/2,get(h2,'ydata')+r/2,25),'edgecolor','none');
	end

	function go(varargin)
		set(tb,'value',0)
		if strcmp(get(tb,'string'),'Start')
			set(tb,'string','Stop','backgroundcolor','r')
		else
			set(tb,'string','Start','backgroundcolor','g')
			return
		end
		set(dummy1,'string','')
		dt=1;
		[hb,h2,hx,h1,r,t,KE]=init;
		rot=str2double(get(ed7,'string'));
		a=str2double(get(ed2,'string'))/2000;
		Yv=-1*str2double(get(ed3,'string'))/2000;
		Xv=str2double(get(ed4,'string'))/2000;
		e=str2double(get(ed5,'string'))/100;
		f=str2double(get(ed6,'string'))*.7/100;
		pc=0;
		c=0;
		set(h2,'visible','off')
		ll=0;
		while strcmp(get(tb,'string'),'Stop')
			if c<5
				am =[
				0    0    0    0    0    1    1    0    0    0    0    0    0    0    0    0    0    1    1    0    0    0    0    0    0
				0    0    0    0    1    1    1    1    0    0    0    0    0    0    0    0    1    1    1    1    0    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    0    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
				0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
				0    0    0    0    1    1    1    1    0    0    0    0    0    0    0    0    1    1    1    1    0    0    0    0    0
				0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
				1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1
				1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1
				1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1
				1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1
				1    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1
				1    0    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    0    1
				0    0    0    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1    0    0
				0    0    0    0    1    1    1    0    0    0    0    0    0    0    0    0    0    0    0    1    1    1    0    0    0
				0    0    0    0    0    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    0    0    0    0    0
				0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    1    1    1    1    0    0    0    0    0    0
				];
			am=[zeros(2,25);am;zeros(2,25)]; %#ok
			colormap([1 1 0;0 0 0])
			set(h1,'facecolor','y')
			p=0;
			else
				c=c-1;
				p=p+1;
			end
			if p>30 && rot<.1
				set(tb,'string','Start','backgroundcolor','g');
				break
			end
			pc=pc+1;
			yc=get(h2,'ydata');
			if get(chk,'value')
				plot(get(h2,'xdata'),yc,'k.')
				pan on
			end
			Yv=Yv+a*dt;
			if yc-(Yv+a*dt)*dt<=r-min(get(h1,'ydata'))
				if Yv>1
					figure('units','normalized','menubar','none','numbertitle','off','color','w')
					set(gcf,'position',[0 -.3 1 1.27])
					axes('position',[0 0 1 1 ])
					t=-4:.01:4;
					sound(cos(300.*t.^2));
					a=imread('face.jpg');
					image(a)
					return
				end
				Yv=Yv+a*dt/(Yv+a*dt)*(yc-r);
				Yv=-Yv*e;
				yt=get(h1,'ydata');
				xxt=get(h1,'xdata');
				ys=sort(yt);
				xr=xxt((ismember(yt,ys(1:10)))==1);
				xrp=str2num(get(dummy1,'string')); %#ok
				set(dummy1,'string',mat2str(xr));
				Xv=Xv*f;
				rot2=rot/10+Xv/r*100;
				Xv=Xv+rot/250;
				rot=rot2;
				% rot=sqrt(10*a*hg/7/r^2)*10;
				if ~isempty(xrp) && sign(Xv+~sign(Xv))*max(xrp)>min(xr)*sign(Xv+~sign(Xv))
					newy=get(pl,'ydata')-rand(size(get(pl,'ydata')))/10*abs(Yv);
					newy(1)=0;
					newy(end)=0;
					set(pl,'ydata',newy);
				else
					sar=zeros(1,8)-(rand(1,8).*(1-f).*abs(Yv));
					ltt=line([min(xr) max(xr)],[0 0]);
					set(ltt,'linewidth',3,'color','w')
					pl=plot(sort(xr),[zeros(1,length(xr)-9),sar,0]);
					set(pl,'linewidth',2.5,'color','k')
				end
 				yt((ismember(yt,ys(1:10))))=yt((ismember(yt,ys(1:10))))+...
					transpose(rand(1,length(yt((ismember(yt,ys(1:10)))))).*(1-f).*abs(Yv));
				set(h2,'ydata',r)
				set(h1,'ydata',yt)
				am =[
					0    0    0    0    0    1    1    0    0    0    0    0    0    0    0    0    0    1    1    0    0    0    0    0    0
					0    0    0    0    1    1    1    1    0    0    0    0    0    0    0    0    1    1    1    1    0    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    1    1    1    1    1    1    0    0    0    0    0    0    1    1    1    1    1    1    0    0    0    0
					0    0    0    0    1    1    1    1    0    0    0    0    0    0    0    0    1    1    1    1    0    0    0    0    0
					0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    0    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    0    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    1    0    0    0    0    1    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    1    0    0    0    0    1    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    0    0    0    0    0    0    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    0    0    0    0    0    0    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    0    0    0    0    0    0    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    1    1    0    0    0    0    0    0    0    0
					0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
					];
				am=[zeros(2,25);am;zeros(2,25)]; %#ok
				colormap([1 0 0;0 0 0])
				set(h1,'facecolor','r')
				c=15;
			else
				for k=1:2
					set(hb(k),'ydata',get(hb(k),'ydata')-Yv*dt)
				end
			end
			set(KE,'xdata',get(KE,'xdata')+Xv*dt)
			set(hb(1),'xdata',get(hb(1),'xdata')+Xv*dt)
			set(hb(2),'xdata',get(hb(2),'xdata')+Xv*dt)
			l=axis;
			if get(h2,'xdata')>l(2)-r && Xv>0
				xlim([hx(1),hx(2)]+l(2))
			end
			if get(h2,'ydata')>l(4)-r && Yv<0
				ylim([hx(3),hx(4)]+l(4))
			end
			if get(h2,'xdata')<l(1)+r && Xv<0
				xlim([l(1),l(2)]-(hx(2)-hx(1)))
			end
			if get(h2,'ydata')<l(3)+r && Yv>0
				ylim([l(3),l(4)]-(hx(4)-hx(3)))
			end
			ccc=(max(get(h1,'ydata'))+min(get(h1,'ydata')))/2;
			ccc2=(max(get(h1,'xdata'))+min(get(h1,'xdata')))/2;
			ll=ll+rot;
			delete(KE)
			KE=surface(flipud(am),'facecolor','texturemap','xdata',linspace(ccc2-r/2,ccc2+r/2,...
				25),'ydata',linspace(ccc-r/2,ccc+r/2,25),'edgecolor','none');
			rotate(KE,[0 0 1],-ll,[ccc2 ccc 1])
			rotate(h1,[0 0 1],-rot,[ccc2,ccc 1])
			drawnow()
		end
	end
end


text:
function data=createdatamatpt(data,E,win)
% Helper function to create an event triggered matrix from a single
% channel of spike times.
% Usage:  data=createdatamatpt(data,E,win)
% Inputs:
% data   (input spike times as a structural array or as a column vector) - required
% E      (events to use as triggers) - required
% win    (window around triggers to use data matrix -[winl winr]) - required
%          e.g [1 1] uses a window starting 1 sec before E and
%              ending 1 sec after E if E and data are in secs.
% Note that E, win and data must have consistent units
% Outputs:
% data      (event triggered data as a structural array - times are stored
% relative to the E-winl
%
if nargin < 3; error('Need all input arguments'); end;
if isstruct(data);
   fnames=fieldnames(data);
   eval(['dtmp=data.' fnames{1} ';'])
else
   dtmp=data(:);
end;
NE=length(E);
winl=win(1);
winr=win(2);
data2(1:NE)=struct('times',[]);
for n=1:NE,
    indx=find(dtmp > E(n)-winl & dtmp<= E(n)+winr);
    if ~isempty(indx)
       data2(n).times=dtmp(indx)-E(n)+winl;
    else
       data2(n).times=[];
    end
end
data=data2;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
fid = fopen('list.txt');
C = textscan(fid, '%s %d');
all_class = 10575;
image_num_in_class = zeros(all_class,1);
image_in_class = cell(all_class,1);
for i = 1:all_class
    image_in_class{i} = find(C{2} == i -1);
    image_num_in_class(i) = length(image_in_class{i});
end;
fclose(fid);
% clear all;
% load('temp1.mat');

test_class = 2000-1;
test_num = sum(C{2}>all_class-test_class);
total = length(C{2});
train_num = total - test_num;
same_r = 0.5;

same_pair_l = floor(test_num / 2 * same_r);
same_in_class = randi(test_class,floor(same_pair_l * 2),1) + all_class-test_class;
same_in_class_unique = unique(same_in_class);
got_in_same = ones(total,1);
got_in_same(train_num+1:end) = 0;
pair1_same = zeros(same_pair_l * 10,1);
pair2_same = zeros(same_pair_l * 10,1);
p=1;
for i = 1: length(same_in_class_unique)
    if mod(p,same_pair_l/100) == 0
        disp([p same_pair_l]);
    end;
    c = same_in_class_unique(i);
    n = sum(same_in_class == same_in_class_unique(i)) * 2;
    s = image_num_in_class(c);
    if n > s
        n = floor(s / 2) * 2;
    end;
    idx_same = randperm(s);
    for j = 1 : n / 2
        pair1_same(p) = image_in_class{c}(idx_same(j*2-1));
        pair2_same(p) = image_in_class{c}(idx_same(j*2));
%         got_in_same(image_in_class{c}(idx_same(j*2-1))) = 1;
%         got_in_same(image_in_class{c}(idx_same(j*2))) = 1;
%         if p==same_pair_l
%             break;
%         end;
        p = p + 1;
    end;
end;
assert(p>same_pair_l);
idx = randperm(p-1);
pair1_same = pair1_same(idx(1:same_pair_l));
pair2_same = pair2_same(idx(1:same_pair_l));
assert(sum(C{2}(pair1_same) == C{2}(pair2_same)) == same_pair_l);
got_in_same(pair1_same) = 1;
got_in_same(pair2_same) = 1;
diff_num = test_num - length(pair1_same) * 2;
diff_sample = find(got_in_same==0);
idx = randperm(diff_num);
diff_sample = diff_sample(idx);
pair1_diff = diff_sample(1:(diff_num/2));
pair2_diff = diff_sample(diff_num/2 + 1:end);
err = find(C{2}(pair1_diff) == C{2}(pair2_diff));
% for i = 1 : length(err)
    ridx = randperm(length(err));
    pair1_diff(err) = pair1_diff(err(ridx));
% end;
assert(sum(C{2}(pair1_diff) == C{2}(pair2_diff)) == 0);
pair1 = [pair1_same;pair1_diff];
pair2 = [pair2_same;pair2_diff];
label_sim = [ones(length(pair1_same),1);zeros(length(pair1_diff),1)];
idx = randperm(length(pair1));
pair1 = pair1(idx);
pair2 = pair2(idx);
label_sim = label_sim(idx);
assert(sum((C{2}(pair1) == C{2}(pair2)) ~= label_sim) == 0);
fid1 = fopen('pair_data1_test.txt','w');
fid2 = fopen('pair_data2_test.txt','w');
fid3 = fopen('sim_label_test.txt','w');
for i = 1 : length(pair1)
    fprintf(fid1,'%s %d\r\n',C{1}{pair1(i)},C{2}(pair1(i)));
    fprintf(fid2,'%s %d\r\n',C{1}{pair2(i)},C{2}(pair2(i)));
    fprintf(fid3,'%d\r\n',label_sim(i));
end;
fclose(fid1);
fclose(fid2);
fclose(fid3);

author: happynear

repo: FaceVerification


text:
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrwvt
%TFRWVT	Unit test for the function TFRWV.

%	F. Auger, December 1995 - O. Lemoine, March 1996.

N=128;

% Covariance by translation in time
t1=60; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrwv(sig1);
tfr2=tfrwv(sig2);
[tr,tc]=size(tfr1);
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrwv test 1 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrwv(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrwv test 2 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrwv(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrwv test 3 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrwv(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrwv test 4 failed');
end

% Frequency-marginal
sig=noisecg(N);
tfr=tfrwv(sig);
FFT=fft(sig);
psd1=abs(FFT(2:N/2)).^2/(2*N);
psd2=mean(tfr(1:2:N,:)')';
if any(abs(psd1-psd2(2:N/2))>sqrt(eps)),
 error('tfrwv test 5 failed');
end

% Compatibility with filtering
h=amgauss(N,N/2,sqrt(N)).*fmlin(N);
x=amgauss(N).*fmconst(N);
y=conv(x,h); y=y(:);  y=y((N/2+1):3*N/2);
tfrx=tfrwv(x);
tfry=tfrwv(y)/(2*N);
tfrh=tfrwv(h);
tfr=zeros(N);
for f=1:N,
 tmp=conv(tfrx(f,:),tfrh(f,:));
 tfr(f,:)=tmp((N/2+1):3*N/2)/N;
end
if any(any(abs(tfr-tfry)>sqrt(eps))),
 error('tfrwv test 6 failed');
end

% Unitarity
x1=amgauss(N).*fmlin(N);
x2=amgauss(N).*fmsin(N);
tfr1=tfrwv(x1);
tfr2=tfrwv(x2);
cor1=abs(x1'*x2)^2/2;
cor2=sum(sum(tfr1.*conj(tfr2)))/N;
if abs(cor1-cor2)>sqrt(eps),
 error('tfrwv test 7 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(N/4,1);fmlin(N/2);zeros(N/4,1)];
tfr=tfrwv(sig);
if sum(any(abs(tfr(:,1:N/4-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(3*N/4+1):N))>sqrt(eps))),
 error('tfrwv test 8 failed');
end

% time localization
t0=30; sig= ((1:N)'==t0);
tfr=tfrwv(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrwv test 9 failed');
end;

% frequency localization
f0=10;
sig=fmconst(N+6,f0/N);
tfr=tfrwv(sig,N/2+2,N);
if (find(tfr>1/N)~=2*f0+1)|(abs(mean(tfr)-1.0)>2.0*eps),
 error('tfrwv test 10 failed');
end;

N=127;

% Covariance by translation in time
t1=60; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrwv(sig1);
tfr2=tfrwv(sig2);
[tr,tc]=size(tfr1);
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrwv test 11 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrwv(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrwv test 12 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrwv(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrwv test 13 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrwv(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrwv test 14 failed');
end

% Frequency-marginal
sig=noisecg(N);
tfr=tfrwv(sig);
FFT=fft(sig);
psd1=abs(FFT(2:fix(N/2)-3)).^2/(2*N);
psd2=mean(tfr(1:2:N,:)')';
if any(abs(psd1-psd2(2:fix(N/2)-3))>5e-2),
 error('tfrwv test 15 failed');
end

% Compatibility with filtering
h=amgauss(N,N/2,sqrt(N)).*fmlin(N);
x=amgauss(N).*fmconst(N);
y=conv(x,h);  y=y(:); y=y((round(N/2)+1):round(3*N/2));
tfrx=tfrwv(x);
tfry=tfrwv(y)/(2*N);
tfrh=tfrwv(h);
tfr=zeros(N);
for f=1:N,
 tmp=conv(tfrx(f,:),tfrh(f,:));
 tfr(f,:)=tmp((round(N/2)+1):round(3*N/2))/N;
end
if any(any(abs(tfr-tfry)>sqrt(eps))),
 error('tfrwv test 16 failed');
end

% Unitarity
x1=amgauss(N).*fmlin(N);
x2=amgauss(N).*fmsin(N);
tfr1=tfrwv(x1);
tfr2=tfrwv(x2);
cor1=abs(x1'*x2)^2/2;
cor2=sum(sum(tfr1.*conj(tfr2)))/N;
if abs(cor1-cor2)>sqrt(eps),
 error('tfrwv test 17 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(round(N/4),1);fmlin(round(N/2));zeros(round(N/4),1)];
tfr=tfrwv(sig);
if sum(any(abs(tfr(:,1:round(N/4)-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(round(3*N/4)+2):N))>sqrt(eps))),
 error('tfrwv test 18 failed');
end

% time localization
t0=30; sig= ((1:N)'==t0);
tfr=tfrwv(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrwv test 19 failed');
end;

% frequency localization
f0=10;
sig=fmconst(N+6,f0/N);
tfr=tfrwv(sig,round(N/2)+2,N);
if (find(tfr>1/N)~=2*f0+1)|(abs(mean(tfr)-1.0)>sqrt(eps)),
 error('tfrwv test 20 failed');
end;

author: HeLiangHIT


text:
function [dataout,condidout] = fl_condidselecthed(data,condid,hedin,hedout);
% Uses hierarchical event description to select and rename conditions

%data is nvar x ntimes x ntrials

% hedin = '*/*/*/13,14,15,16,54,55,60,64,65,69'; %keep only conditions 13,14, etc
% hedout = '*/-/-/-'; % condidout should only have the first descriptor

%     modifier = 'animate/*/*'
%     modifier = 'animate/*/1,2,3,4,5,6'
%     modifier = 'animate/*/7,8'
%     modifier = 'inanimate/*/*/54,55,60,64,65,69'

%% split hierarchical event descriptor (hed) tag
hedinsplit = strsplit(hedin,'/');
hedoutsplit = strsplit(hedout,'/');
m = length(hedinsplit);

%% split condid tags
n = length(condid);
for i = 1:n
    condidsplit{i} = strsplit(condid{i},'/');
end
keep = ones(n,1);

%% modify condid based on hedin
for i = 1:m

    if ~strcmp(hedinsplit{i},'*')% if not '*', then a specific tag selection is given

        hedinspliti = strsplit(hedinsplit{i},','); %tags to keep (separated by comma)
        for j = 1:n
            if keep(j)
                try %just in case some events's have less elements than hed modifier
                    if ~nnz(strcmp(condidsplit{j}{i},hedinspliti)) %if requested tag does not exist
                        keep(j) = 0; %discard this item
                    end
                end
            end
        end
    end

end

%% modify condid based on hedout
for i = 1:m

    if strcmp(hedoutsplit{i},'*') %pass input to output
        continue
    elseif strcmp(hedoutsplit{i},'-') %if remove (tag '-')
        for j = 1:n
            condidsplit{j}{i} = '-';
        end
    else %if an output name is given
        for j = 1:n
            condidsplit{j}{i} = hedoutsplit{i};
        end
    end

end

%% join condid tags
for i = 1:n
    condidout{i} = strjoin(condidsplit{i},'/');
end

%% keep only data with tags sharing descriptors with hed
dataout = data(:,:,keep==1);
condidout = condidout(keep==1);

return

author: brainstorm-tools


text:
classdef TestCalcHist
    %TestCalcHist

    properties (Constant)
        img = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_1
            % 2D histogram
            im = imread(TestCalcHist.img);
            histSize = [30, 32];
            edges1 = linspace(0, 256, histSize(1)+1);
            edges2 = linspace(0, 256, histSize(2)+1);
            edges = {edges1, edges2};

            H = cv.calcHist(im(:,:,[1 2]), edges);
            validateattributes(H, {'single'}, {'nonsparse', 'size',histSize});

            HH = cv.calcHist({im(:,:,1), im(:,:,2)}, edges);
            validateattributes(HH, {'single'}, {'nonsparse', 'size',histSize});
            assert(isequal(H,HH));

            HH = cv.calcHist(im, edges, 'Channels',[1 2]-1);
            validateattributes(HH, {'single'}, {'nonsparse', 'size',histSize});
            assert(isequal(H,HH));

            HH = cv.calcHist(im, edges, 'Channels',[1 2]-1, 'HistSize',histSize);
            validateattributes(HH, {'single'}, {'nonsparse', 'size',histSize});
            assert(isequal(H,HH));

            HH = cv.calcHist(im, {edges1([1 end]), edges2([1 end])}, ...
                'Uniform',true, 'Channels',[1 2]-1, 'HistSize',histSize);
            validateattributes(HH, {'single'}, {'nonsparse', 'size',histSize});
            assert(isequal(H,HH));

            HH(:) = 0;
            HH = cv.calcHist(im(:,:,[1 2]), edges, 'Hist',HH);
            validateattributes(HH, {'single'}, {'nonsparse', 'size',histSize});
            assert(isequal(H,HH));

            HH = cv.calcHist(im(:,:,[1 2]), edges, 'Sparse',true);
            validateattributes(HH, {'double'}, {'size',histSize});
            assert(issparse(HH) && isequal(H,full(HH)));
        end

        function test_histc
            % compare against HISTC
            im = cv.imread(TestCalcHist.img, 'Grayscale',true);  % uint8 grayscale
            edges = [0 50 100 150 200 256];                      % 1D histogram
            H1 = cv.calcHist(im, edges);
            H2  = histc(im(:), edges);
            assert(isequal(H1, H2(1:end-1)));

            % with a mask
            mask = false(size(im));
            mask(100:300,100:300) = true;
            H1 = cv.calcHist(im, edges, 'Mask',mask);
            H2  = histc(im(mask), edges);
            assert(isequal(H1, H2(1:end-1)));
        end

        function test_histcounts
            % compare against the new HISTCOUNTS function
            if mexopencv.isOctave() || verLessThan('matlab','8.4')
                error('mexopencv:testskip', 'toolbox');
            end
            im = cv.imread(TestCalcHist.img, 'Grayscale',true);  % uint8 grayscale
            edges = [0 50 100 150 200 256];                      % 1D histogram
            H1 = cv.calcHist(im, edges);
            H2  = histcounts(im, [edges(1:end-1) 255]);
            assert(isequal(H1(:), H2(:)));

            % with a mask
            mask = false(size(im));
            mask(100:300,100:300) = true;
            H1 = cv.calcHist(im, edges, 'Mask',mask);
            H2  = histcounts(im(mask), [edges(1:end-1) 255]);
            assert(isequal(H1(:), H2(:)));
        end

        function test_nd_hist
            % multi-dimensional histogram (4-D)
            X = rand([400,500,4], 'single');
            [h,w,~] = size(X);
            edges = {[0 1], [0 1], [0 1], [0 1]};
            histSize = [7 8 9 10];
            H = cv.calcHist(X, edges, 'HistSize',histSize, 'Uniform',true);
            validateattributes(H, {'single'}, {'ndims',4, 'size',histSize});
            assert(isequal(sum(H(:)), h*w));
        end

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

end

author: kyamagu

repo: mexopencv


text:
%% PRIME_BATCH_LOCAL uses the BATCH command to run the PRIME code locally.
%
%  Discussion:
%
%    The PRIME code is a function, so first we must write a script
%    called PRIME_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_BATCH_LOCAL:\n' );
  fprintf ( 1, '  Run PRIME_SCRIPT locally.\n' );
%
%  BATCH sends the script for execution.
%
  my_job = batch ( 'prime_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'prime_fun' }, ...
    'matlabpool', 4 );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  LOAD makes the script's workspace available.
%
%  total = total number of primes.
%
  load ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Total number of primes = %d\n', total );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_BATCH_LOCAL:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
%{
    Tests the extended "analysis" form of basis pursuit de-noising

    min_1 alpha*||W_1 x||_1 + beta*|| W_2 x ||_1
s.t.
    || A(x) - b || <= eps

The solvers solve a regularized version

see also test_sBP.m and test_sBPDN.m and test_sBPDN_W.m

%}

% Before running this, please add the TFOCS base directory to your path
myAwgn = @(x,snr) x + ...
        10^( (10*log10(sum(abs(x(:)).^2)/length(x(:))) - snr)/20 )*randn(size(x));

% Try to load the problem from disk
fileName = fullfile('reference_solutions','basispursuit_WW_problem1_smoothed_noisy');
randn('state',34324);
rand('state',34324);

N = 512;
M = round(N/2);
K = round(M/5);

A = randn(M,N);
x = zeros(N,1);

% introduce a sparsifying transform "W"
d  = round(4*N);     % redundant
Wf = @(x) dct(x,d);  % zero-padded DCT
downsample  = @(x) x(1:N,:);
Wt = @(y) downsample( idct(y) );    % transpose of W
W  = Wf(eye(N));     % make an explicit matrix for CVX

% and add another transform:
d2 = round(2*N);
W2 = randn(d2,N);
if exist([fileName,'.mat'],'file')
    load(fileName);
    fprintf('Loaded problem from %s\n', fileName );
else

    % Generate a new problem

    alpha = 1;
    beta = 2;

    % the signal x consists of several pure tones at random frequencies
    for k = 1:K
        x = x + randn()*sin( rand()*pi*(1:N) + 2*pi*rand() ).';
    end

    b_original = A*x;
    snr = 40;  % SNR in dB
    b   = myAwgn(b_original,snr);
    EPS = norm(b-b_original);
    x_original = x;

    mu = .01*norm(Wf(x),Inf);
    x0 = zeros(N,1);

    % get reference via CVX
    tic
    cvx_begin
        cvx_precision high
        variable xcvx(N,1)
        minimize alpha*norm(W*xcvx,1) + beta*norm(W2*xcvx,1) + ...
            mu/2*sum_square(xcvx-x0)
        subject to
            norm(A*xcvx - b ) <= EPS
    cvx_end
    time_IPM = toc;
    x_ref = xcvx;
    objF        = @(x)alpha*norm(W*x,1) +beta*norm(W2*x,1)+ mu/2*norm(x-x0).^2;
    obj_ref     = objF(x_ref);

    save(fileName,'x_ref','b','x_original','mu',...
        'EPS','b_original','obj_ref','x0','time_IPM','snr',...
        'd','d2','alpha','beta');
    fprintf('Saved data to file %s\n', fileName);

end

[M,N]           = size(A);
norm_x_ref      = norm(x_ref);
norm_x_orig     = norm(x_original);
er_ref          = @(x) norm(x-x_ref)/norm_x_ref;
er_signal       = @(x) norm(x-x_original)/norm_x_orig;
resid           = @(x) norm(A*x-b)/norm(b);  % change if b is noisy

fprintf('\tA is %d x %d\n', M, N );
fprintf('\tl1-norm solution and original signal differ by %.2e (mu = %.2e)\n', ...
    norm(x_ref - x_original)/norm(x_original),mu );

%% Call the TFOCS solver
objF            = @(x)alpha*norm(W*x,1) +beta*norm(W2*x,1)+ mu/2*norm(x-x0).^2;
infeasF         = @(x)norm(A*x-b) - EPS;
er              = er_ref;  % error with reference solution (from IPM)
opts            = [];
opts.errFcn     = { @(f,dual,primal) er(primal), ...
                    @(f,dual,primal) obj_ref - f, ...
                    @(f,dual,primal) infeasF(primal) };
opts.maxIts     = 2000;
opts.tol        = 1e-10;
% opts.normA2     = norm(A*A');
% opts.normW2     = norm(W'*W);
z0  = [];   % we don't have a good guess for the dual
tic;
[ x, out, optsOut ] = solver_sBPDN_WW( A, alpha,W,beta,W2,b, EPS, mu, x0, z0, opts );
time_TFOCS = toc;
fprintf('x is sub-optimal by %.2e, and infeasible by %.2e\n',...
    objF(x) - obj_ref, infeasF(x) );

fprintf('Solution has %d nonzeros.  Error vs. IPM solution is %.2e\n',...
    nnz(x), er(x) );

% Check that we are within allowable bounds
if out.err(end,1) < 1e-3
    disp('Everything is working');
else
    error('Failed the test');
end

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

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS