text:
function [Iq,in_avg,ext_avg,in_min,ext_max]=icassoStability(sR,L,graphmode)
%function Iq=icassoStability(sR,L,graphmode)
%
%PURPOSE
%
%To compute and/or plot the stability (quality) indices of the ICA
%estimate-clusters.
%
% Iq=avg(intra-cluster similarity) - avg(extra-cluster similarity)
%
%See publication Himberg et al. (2004), "Validating the independent
%components of neuroimaging time-series via clustering and
%visualization". NeuroImage, 22:3(1214-1222). Ideally, each ICA
%estimate-cluster should have Iq=1. The smaller the value is, the
%less stable (compact and isolated) the estimate-cluster is.
%
%EXAMPLES OF BASIC USAGE
%
%   Iq=icassoStability(sR);
%
%returns the stability index for each estimate (centrotype) using
%the default number of estimate-clusters.
%
%   Iq=icassoStability(sR,13,'plotindex');
%
%...same as previous but also plots the values (in rank order) and
%uses 13 estimate clusters instead of default.
%
%INPUT
%
%[An argument in brackets is optional. If it isn't  given or it's
% an empty matrix/string, the function will use a default value.]
%
% sR          (struct) Icasso result data structure; clustering
%              must be readily computed.
% [L]         (scalar) number of clusters, default: (reduced) data
%              dimension
% [graphmode] (string) 'none' (default) | 'plotindex' | 'plotstat'
%
%OUTPUTS
%
% Iq  (Lx1 matrix) Iq(C) contains the index value for estimate C
% (actually estimate-cluster C)
%
%DETAILS
%
%The intra-cluster similarities for cluster C mean the mutual
%similarities between the estimates in C, and the extra-cluster
%similarities for C are the similarities between the estimates in C
%and the estimates not in C. The index is the average intra-cluster
%similarity subtracted by the average intra-cluster similarity.
%
%In default graph mode 'none'
% The function computes the index and returns it. No graphical
% output. Iq(C) contains the index for cluster C (That is, for
% estimates sR.cluster.partition(L,:)==C.
%
%In graph mode 'plotindex'
% The function computes the quality index as in mode 'none'
% and plots the index in the active axis (gca) as follows:
% The X-axis shows the value of the index for clusters. The clusters
% are ranked on the Y-axis in descending order according to the
% index. The Y-axis legend shows the cluster number (the same as in
% sR.cluster.partition(L,:)). The values are indicated by black
% dots connected into each other by a dotted line.
%
%In graph mode 'plotstat'
% As previous one but instead of the index,
% some cluster statistics are shown:
%  Light red  indicates the average intra-cluster similarity for a cluster
%  Light blue               average extra-cluster
%  Clear red                minimum intra-cluster
%  Clear blue               maximum extra-cluster
%
%SEE ALSO
% clusterquality
% icassoViz
% icassoShow

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 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.

% ver 1.2 100501 johan

if nargin<2|isempty(L),
  L=icassoGet(sR,'rdim');
  disp([sprintf('\n') 'Number of estimate-clusters not given: using (reduced)' ...
	   ' data dimension.']);
end

if nargin<3 | isempty(graphmode),
  graphmode='none';
end

clusternumber=1:L;
partition=sR.cluster.partition(L,:);

for i=1:L,
  NofEstimates(i)=sum(partition==i);
end

% compute cluster quality index

[Iq,in_avg,ext_avg]=clusterquality('mean',sR.cluster.similarity,partition);
[tmp,in_min,ext_max]=clusterquality('minmax',sR.cluster.similarity,partition);

% Sort entries according to score

[tmp,order]=sort(-Iq);
score=Iq(order);
in_avg=in_avg(order);
ext_avg=ext_avg(order);
in_min=in_min(order);
ext_max=ext_max(order);
clusternumber=clusternumber(order);
NofEstimates=NofEstimates(order);

% make tick mark labels
for i=1:L,
  ticklabel{i}=num2str(clusternumber(i));
end

% Plot
switch lower(graphmode)
 case 'none'
  return;
 case 'plotindex'
  cla reset;
  h_score=plot(score,1:L,'o:');
  set(h_score,'color','k','markersize',8,'markerfacecolor','k');
  set(gca,'ytick',1:L,'yticklabel',ticklabel,'ydir','reve');
  ax=axis;

  if ax(1)<0,
    axis([ax(1) 1  0.5 L+.5]);
  else
    axis([0 1  0.5 L+.5]);
  end

  grid on;
  xlabel('I_q=avg(S(i)_{int})-avg(S(i)_{ext})');
  title('Stability index (I_q) for ICA estimate-clusters');
 case 'plotstat'
  cla reset;

  h_in_avg=barh(in_avg,0.7); hold on
  h_ext_max=barh(-ext_max,0.35);
  h_in_min=barh(in_min,0.35);
  h_ext_avg=barh(-ext_avg,0.7);

  axis on;
  set(h_in_avg,'facecolor',[1 0.7 0.7]); hold on;
  set(h_in_min,'facecolor',[1 0 0]);
  set(h_ext_avg,'facecolor',[0.7 0.7 1]);
  set(h_ext_max,'facecolor',[0 0 1]);
  axis([-1 1 0.5 L+.5]);

  set(gca,'ytick',1:L,'yticklabel',ticklabel,...
	  'ydir','reve','xtick',[-1 -.5 0 .5 1],...
	  'xgrid','on',...
	  'xticklabel',{'1' '0.5 ' '0' '0.5' '1'});

  title('Statistics on within- and between-cluster similarities');
  xlabel('clusters are ordered according to I_q=avg(S(i)_{int})-avg(S(i)_{ext})');
  legend([h_in_avg(1),h_ext_avg(1),h_in_min(1),h_ext_max(1)],...
	 {'mean S_{in}' 'mean S_{ex}' 'min S_{in}' 'max S_{ex}'},-1);
 otherwise
  error('Graphmode must be ''none'',''plotindex'' or ''plotstat''.');
end

ylabel('Label')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function value = r4_gamma ( x )

%*****************************************************************************80
%
%% R4_GAMMA evaluates the gamma function of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the gamma function of X.
%
  persistent dxrel
  persistent gcs
  persistent ngcs
  persistent sq2pil
  persistent xmax
  persistent xmin
  persistent xsml

  sq2pil = 0.91893853320467274;

  if ( isempty ( ngcs ) )

    gcs = [ ...
      0.008571195590989331, ...
      0.004415381324841007, ...
      0.05685043681599363, ...
     -0.004219835396418561, ...
      0.001326808181212460, ...
     -0.0001893024529798880, ...
      0.0000360692532744124, ...
     -0.0000060567619044608, ...
      0.0000010558295463022, ...
     -0.0000001811967365542, ...
      0.0000000311772496471, ...
     -0.0000000053542196390, ...
      0.0000000009193275519, ...
     -0.0000000001577941280, ...
      0.0000000000270798062, ...
     -0.0000000000046468186, ...
      0.0000000000007973350, ...
     -0.0000000000001368078, ...
      0.0000000000000234731, ...
     -0.0000000000000040274, ...
      0.0000000000000006910, ...
     -0.0000000000000001185, ...
      0.0000000000000000203 ]';

    ngcs = r4_inits ( gcs, 23, 0.1 * r4_mach ( 3 ) );
    [ xmin, xmax ] = r4_gaml ( );
    xsml = exp ( max ( log ( r4_mach ( 1 ) ), ...
      - log ( r4_mach ( 2 ) ) ) + 0.01 );
    dxrel = sqrt ( r4_mach ( 4 ) );

  end

  y = abs ( x );

  if ( y <= 10.0 )

    n = r4_aint ( x );
    if ( x < 0.0 )
      n = n - 1;
    end
    y = x - n;
    n = n - 1;
    value = 0.9375 + r4_csevl ( 2.0 * y - 1.0, gcs, ngcs );

    if ( n == 0 )

      return

    elseif ( n < 0 )

      n = - n;

      if ( x == 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R4_GAMMA - Fatal error!\n' );
        fprintf ( 1, '  X is 0.\n' );
        error ( 'R4_GAMMA - Fatal error!' )
      end

      if ( x < 0.0 && x + n - 2 == 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R4_GAMMA - Fatal error!\n' );
        fprintf ( 1, '  X is a negative integer.\n' );
        error ( 'R4_GAMMA - Fatal error!' )
      end

      if ( x < - 0.5 && ...
        abs ( ( x - r4_aint ( x - 0.5 ) ) / x ) < dxrel )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R4_GAMMA - Warning!\n' );
        fprintf ( 1, '  X too near a negative integer,\n' );
        fprintf ( 1, '  answer is half precision.\n' );
      end

      if ( y < xsml )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R4_GAMMA - Fatal error!\n' );
        fprintf ( 1, '  X is so close to zero that Gamma overflows.\n' );
        error ( 'R4_GAMMA - Fatal error!' )
      end

      for i = 1 : n
        value = value / ( x + i - 1 );
      end

    elseif ( n == 0 )

    else

      for i = 1 : n
        value = ( y + i ) * value;
      end

    end

  else

    if ( xmax < x )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R4_GAMMA - Fatal error!\n' );
      fprintf ( 1, '  X so big that Gamma overflows.\n' );
      error ( 'R4_GAMMA - Fatal error!' )
    end
%
%  Underflow.
%
    if ( x < xmin )
      value = 0.0;
      return
    end

    value = exp ( ( y - 0.5 ) * log ( y ) - y + sq2pil + r4_lgmc ( y ) );

    if ( 0.0 < x )
      return
    end

    if ( abs ( ( x - r4_aint ( x - 0.5 ) ) / x ) < dxrel )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R4_GAMMA - Warning!\n' );
      fprintf ( 1, '  X too near a negative integer,\n' );
      fprintf ( 1, '  answer is half precision.\n' );
    end

    sinpiy = sin ( pi * y );

    if ( sinpiy == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R4_GAMMA - Fatal error!\n' );
      fprintf ( 1, '  X is a negative integer.\n' );
      error ( 'R4_GAMMA - Fatal error!' )
    end

    value = - pi / ( y * sinpiy * value );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Load, modify and export a fig file as an eps file.

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;  % frequency
Vm = 10; % peak
phi = 0; % phase

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v = Vm*sin(th+phi);

% plot it
figure;
plot(t*1E3, v);

% change settings
opt = [];
opt.XLabel = 'Time, t (ms)'; % xlabel
opt.YLabel = 'Voltage, V (V)'; %ylabel
opt.Title = 'Voltage as a function of time';

% Save? comment the following line if you do not want to save
opt.FileName = 'plotSimple1.eps';

% apply the settings
setPlotProp(opt);

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
function [centroid, dimension_weight, class] = Subspace_Kmeans(data, iteration, K, beta, verbose)

% Pre-allocate weights
centroid = zeros(K,size(data,2)); % centroid for each class
class = zeros(size(data,1),1); % classification result for each observation
dimension_weight = ones(K,size(data,2)); % weights for each dimension in each class

% Initialization
n = size(data,1);
J = 0; % used to record objective function value

% First round
centroid = data(unidrnd(n,K,1),:); % randomly choose K initial centroids
dimension_weight = 1/size(data,2) * dimension_weight; % initialize weights for each dimension using uniform distribution

% Partially optimization
for i = 1:iteration

    J0 = J;

    % body part
    class = classify(centroid,dimension_weight,data);
    centroid = centroid_update(centroid,data,class);
    [dimension_weight,J] = dimension_weight_update(K,centroid,class,dimension_weight,data,beta);

    % whether display intermediate information
    if(verbose == 1)
        disp(['Objective Function(J) Value: ',num2str(J)]);
    end

    % early stop condition
    if(abs(J-J0)<1e-9)
        fprintf('*** Clustering terminates after %i iterations ***\n',i);
        break;
    end
end
end

% Details on the scheme of updating dimension weights
function [alpha,J] = dimension_weight_update(K,m,c,alpha,data,beta)
all_matrix = zeros(size(data,1),size(data,2));
for i = 1:size(data,1)
    all_matrix(i,:) = (m(c(i,1),:)-data(i,:)) .* (m(c(i,1),:)-data(i,:));
end
J = 0;

for i = 1:K
    class_matrix = all_matrix(c==i,:); % pick out data in all_matrix that corresponds to class i
    t_sum = sum(class_matrix,1) + 1e-9; % adding a small positive constant to make weights computatable
    J = J + (alpha(i,:).^beta) * t_sum';
    if(size(class_matrix,1)>0)
        for j = 1:size(alpha,2)
            % variables used to avoid redundant calculation
            tt_numerator = t_sum(1,j);
            t_denominator = sum(tt_numerator./t_sum,2);
            alpha(i,j) = 1/((t_denominator+1e-9)^(1/(beta-1)));
        end
    else
        continue;
    end
end

% Normalize
for i = 1:size(alpha,1)
    alpha(i,:) = alpha(i,:)./(sum(alpha(i,:),2));
end

end

function [result] = classify(m,alpha,data)

result = zeros(size(data,1),1);

% Construct temporary matrix for efficiently computing dimention-weighted distance
matrix = zeros(size(m,1),size(data,2));

for i = 1:size(result,1)
    for j = 1:size(m,1)
        matrix(j,:) = (m(j,:)-data(i,:)) .* (m(j,:)-data(i,:));
    end
    temp = sum(alpha .* matrix,2);

    % To avoid more than one class having the minimum distance.
    t_index = find(temp == min(temp));
    result(i,1) = t_index(1,1);
end
end

function [m] = centroid_update(m,data,c)
for i = 1:size(m,1)
    if(~isempty(data(c==i,:)))
        m(i,:) = mean(data(c==i,:),1);
    else
        continue;
    end
end
end

author: xuyxu

repo: Clustering

save_path: github-repos/MATLAB/xuyxu-Clustering


text:
function [errorcode, varargout] = common_size (varargin)

% COMMON_SIZE  Checks that all inputs are either scalar or of common size
%  [ERR, Y1, ...] = common_size(X1, ...)
%  Determine if all input arguments are either scalar or of common
%  size.  If so, ERR is zero, and YI is a matrix of the
%  common size with all entries equal to XI if this is a scalar or
%  XI otherwise.  If the inputs cannot be brought to a common size,
%  errorcode is 1, and YI is XI.
%
%  Example:
%   [errorcode, a, b] = common_size([1 2; 3 4], 5)
%      >> errorcode = 0
%      >> a = [ 1, 2; 3, 4 ]
%      >> b = [ 5, 5; 5, 5 ]

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
% Copyright (C) 1995-2007 Kurt Hornik
% Copyright (C) 2008 Dynare Team
% Copyright (C) 2021 Geomatics Research & Development srl (GReD)
% Adapted from Octave
%  Written by:       KH <Kurt.Hornik@wu-wien.ac.at>
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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 3 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, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

if (nargin < 2)
    error ('common_size: only makes sense if nargin >= 2');
end

len = 2;
for i = 1 : nargin
    sz =  size (varargin{i});
    if (length (sz) < len)
        s(i,:) = [sz, ones(1,len - length(sz))];
    else
        if (length (sz) > len)
            if (i > 1)
                s = [s, ones(size(s,1), length(sz) - len)];
            end
            len = length (sz);
        end
        s(i,:) = sz;
    end
end

m = max (s);
if (any (any ((s ~= 1)') & any ((s ~= ones (nargin, 1) * m)')))
    errorcode = 1;
    varargout = varargin;
else
    errorcode = 0;
    for i = 1 : nargin
        varargout{i} = varargin{i};
        if (prod (s(i,:)) == 1)
            varargout{i} = varargout{i} * ones (m);
        end
    end
end

author: goGPS-Project


text:
function [ res ] = point_in_rect( pt, rect )
% Tests if the point pt is within rectangle rect
%  rect is [left_x top_y right_x bottom_y]

res = (pt(:,1) >= rect(1)) & (pt(:,2) >= rect(2)) & (pt(:,1) <= rect(3)) & (pt(:,2) <= rect(4));

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% LFDispLawnmower - visualize a 4D light field animating a lawnmower path through two dimensions
%
% Usage:
%     FigureHandle = LFDispLawnmower( LF )
%     FigureHandle = LFDispLawnmower( LF, [RenderOptions], [ScaleFactor] )
%
%
% Inputs:
%
%     LF : a colour or single-channel light field, and can a floating point or integer format. For
%          display, it is converted to 8 bits per channel. If LF contains more than three colour
%          channels, as is the case when a weight channel is present, only the first three are used.
%
% Optional Inputs:
%
%     RenderOptions
%        .FrameDelay         : delay between frames, in seconds. Default 1/30.
%        .NumCycles          : how many times the path should repeat. Default inf.
%
%     ScaleFactor : Adjusts the size of the display -- 1 means no change, 2 means twice as
%                   big, etc. Integer values are recommended to avoid scaling artifacts.
%
% Outputs:
%
%     FigureHandle
%
%
% User guide: LFToolbox.pdf
% See also: LFDisp, LFDispMousePan, LFDispVidCirc, LFDispProj, LFDispTiles

% Copyright (c) 2013-2020 Donald G. Dansereau

% todo: smooth transition from vertical to horizontal
% todo: optional horz/vert only

function FigureHandle = LFDispLawnmower( LF, RenderOptions, ScaleFactor )

%---Defaults---
ScaleFactor = LFDefaultVal('ScaleFactor', 1);
RenderOptions = LFDefaultField( 'RenderOptions', 'FrameDelay', 1/30 );
RenderOptions = LFDefaultField( 'RenderOptions', 'NumCycles', inf );

%---Check for a light field with only one spatial dim---
if( ndims(LF) == 3 || (ndims(LF)==4 && (size(LF,4)==3 || size(LF,4)==4)) )
	LF = reshape(LF, [1, size(LF)]);
end
LFSize = size(LF);

%---Check for mono and clip off the weight channel if present---
HasWeight = (ndims(LF)>2 && (LFSize(end)==2 || LFSize(end)==4));
if( HasWeight )
	LF = LF(:,:,:,:,1:end-1);
end

%---Rescale for 8-bit display---
if( isfloat(LF) )
	LF = uint8(LF ./ max(LF(:)) .* 255);
else
	LF = uint8(LF.*(255 / double(intmax(class(LF)))));
end

%---Setup the motion path---
[TSize,SSize, ~,~] = size(LF(:,:,:,:,1));
[tt,ss] = ndgrid(1:TSize, 1:SSize);

%---Setup the display---
TIdx = 1;
SIdx = 1;
[ImageHandle,FigureHandle] = LFDispSetup( squeeze(LF(TIdx,SIdx,:,:,:)), ScaleFactor );

while( 1 )
	CurRep = 0;

	% left/right
	TVec = tt;
	SVec = ss;
	SVec(2:2:end,:) = SVec(2:2:end,end:-1:1);
	SVec = SVec';
	TVec = TVec';

	for( PoseIdx = 1:numel(tt) )

		TIdx = TVec(PoseIdx);
		SIdx = SVec(PoseIdx);

		CurFrame =  squeeze(LF( TIdx, SIdx, :,:,: ));
		set(ImageHandle,'cdata', CurFrame );

		pause(RenderOptions.FrameDelay)

	end

	% Up/down
	TVec = tt;
	TVec(:,2:2:end) = TVec(end:-1:1, 2:2:end);
	SVec = ss;

	for( PoseIdx = 1:numel(tt) )

		TIdx = TVec(PoseIdx);
		SIdx = SVec(PoseIdx);

		CurFrame =  squeeze(LF( TIdx, SIdx, :,:,: ));
		set(ImageHandle,'cdata', CurFrame );

		pause(RenderOptions.FrameDelay)

	end
	CurRep = CurRep + 1;

	if( CurRep >= RenderOptions.NumCycles )
		break;
	end
end

repo: LFToolbox


text:
function density = p14_density ( m, n, point )

%*****************************************************************************80
%
%% P14_DENSITY returns the density for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real DENSITY(N), the mesh density at each point.
%
  density = ones ( 1, n );

  return
end


text:
clear all; close all; clc

%% folder setup
isSaveFig = 1;
isPlotFig = 1;

outputDir = GetOutputDataDir;

ClusterIDs = [2,1]; % init; can overrride
prct_const = 2; % init; can overrride

%% init

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

setappdata(hfig,'isMotorseed',1);

%%
IM_scatter = cell(2,18); % L and R

IM_maps = cell(5,18);
% IM_Xonly = cell(n_reg,18);
% IM_Yonly = cell(n_reg,18);
% IM_XY = cell(n_reg,18);
% IM_X = cell(n_reg,18);
% IM_Y = cell(n_reg,18);

caseID = 1;
switch caseID
    case 1
        load(fullfile(outputDir,'4D_SM_betas.mat'));
        M_reg_name = {'2x2motormaps'};
        range_fish = GetFishRange;% init; can overrride
    case 2
        load(fullfile(outputDir,'4D_SM_stimrangePT_betas.mat'));
        M_reg_name = {'2x2motormaps_PT'};
        range_fish = 6:18;
    case 3
        load(fullfile(outputDir,'4D_SM_stimrangeOMR_betas.mat'));
        M_reg_name = {'2x2motormaps_OMR'};
        range_fish = 8:18;
    case 4
        load(fullfile(outputDir,'4D_SM_stimrangelooming_betas.mat'));
        M_reg_name = {'2x2motormaps_looming'};
        range_fish = [9:15,17:18];
    case 5
        load(fullfile(outputDir,'4D_SM_stimrangeDF_betas.mat'));
        M_reg_name = {'2x2motormaps_DF'};
        range_fish = [12:15,17:18];
end
%% run fish
for i_fish = range_fish

    cIX_all = LoadSingleFishDefault(i_fish,hfig,ClusterIDs,[],0);

    %% main loop for left/right motor
    M_pass = cell(5,2);

    for i_lr = 1:2
        %% scatter plot with 4D components
        % loaded betas for this fish for the combo data (including both stim)
        betas = Betas{i_lr,i_fish};
        % set up plot dimensions
        X = betas(:,1);%b3;%b1;
        Y = betas(:,2);
        %                 Y = sqrt(b2.^2+b3.^2);%b2;
        Xname = 'motor res.';
        Yname = 'motor avr.';

        numcell = length(X);

        A = X;
        topN = round(prct_const/100*numcell); % top x% cutoff
        [~,IX] = sort(A,'descend');
        thresA = A(IX(topN));

        B = Y;
        topN = round(prct_const/100*numcell);
        [~,IX] = sort(B,'descend');
        thresB = B(IX(topN));

        IX_passX = find(A>=thresA);
        IX_passY = find(B>=thresB);
        IX_passXonly = setdiff(IX_passX,IX_passY);
        IX_passYonly = setdiff(IX_passY,IX_passX);%find(B>=thresB);%
        IX_passXY = intersect(IX_passX,IX_passY);
        IX_pass = union(IX_passX,IX_passY);
        IX_fail = intersect(find(A<thresA),find(B<thresB));%find(A<thresA);

        % get min/max
        x0 = min(X(IX_pass));
        x1 = max(X(IX_pass));
        y0 = min(Y(IX_pass));
        y1 = max(Y(IX_pass));

        gIX_in = (1:length(X))';

        M_pass{1,i_lr} = IX_passX;
        M_pass{2,i_lr} = IX_passY;
        M_pass{3,i_lr} = IX_passXonly;
        M_pass{4,i_lr} = IX_passYonly;
        M_pass{5,i_lr} = IX_passXY;
        %         PassX_2{i_lr} = IX_passX;
        %         PassY_2{i_lr} = IX_passY;
        %         PassXonly_2{i_lr} = IX_passXonly;
        %         PassYonly_2{i_lr} = IX_passYonly;
        %         PassXY_2{i_lr} = IX_passXY;

        %% scatter plot
        clrX = [0.3,0.8,0];
        clrY = [0.9,0.2,0.9];
        clrXY = [0,0.3,0.9];
        clr_fail = [0.5,0.5,0.5];

        h = figure('Position',[500,100,300,250]); hold on
        scatter(X,Y,1,clr_fail,'filled')

        scatter(X(IX_passXonly),Y(IX_passXonly),1,clrX);%,'filled');
        scatter(X(IX_passYonly),Y(IX_passYonly),1,clrY);%,'filled');
        scatter(X(IX_passXY),Y(IX_passXY),1,clrXY);%[1,0.5,0.5]);%,'filled');

        plot([x0,x1],[thresB,thresB],'k--');
        plot([thresA,thresA],[y0,y1],'k--');

        xlabel(Xname);ylabel(Yname);
        axis equal
        %         set(gca,'XTick',0:0.5:1,'YTick',-0.5:0.5:1);
        axis tight
        %% save bubble plot
        IM_scatter{i_lr,i_fish} = print('-RGBImage');
        close(h)

    end

    %% component anat map
    setappdata(hfig,'clrmap_name','hsv_old');
    for i_map = 1:5
        cIX12 = M_pass(i_map,:);
        cIX = vertcat(cIX12{1},cIX12{2});
        gIX = vertcat(ones(size(cIX12{1})),2*ones(size(cIX12{2})));
        I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
        [h,~,im] = DrawCellsOnAnat(I);
        close(h);
        IM_maps{i_map,i_fish} = im;

        %export code - added by Joe and Misha
        cIX_abs = I.absIX(I.cIX); %is referencing through absIX correct?
        M_xyz = I.CellXYZ(cIX_abs,:);
        % Determine size from anat sizes, apparently yxz
        siz_ = [size(I.anat_yx,1), size(I.anat_yx,2), size(I.anat_zx,1)];
        func_map=zeros(siz_(1),siz_(2),siz_(3));
        anat_map=zeros(siz_(1),siz_(2),siz_(3));
        for i=1:size(M_xyz,1),func_map(M_xyz(i,1),M_xyz(i,2),M_xyz(i,3))=1;end
        for i=1:size(I.CellXYZ,1),anat_map(I.CellXYZ(i,1),I.CellXYZ(i,2),I.CellXYZ(i,3))=1;end
        eval(['save export_fig3_fish' num2str(i_fish) '_case' num2str(caseID) '_mapnum' num2str(i_map) ' func_map anat_map'])
        %mapnums are above threshold for  X, Y, onlyX, onlyY, X and Y

    end

end

%% save as tiff stack
M_lr = {'-L','-R'};

M_comp_names = {'passX','passY','passXonly','passYonly','passXY'};
n_reg = 1;
for i_set = 1:n_reg
    range_im = range_fish;%M_fishrange_im{i_set};

    for i_lr = 1:2
        tiffdir = fullfile(outputDir,[M_reg_name{i_set},M_lr{i_lr},'_scatter_allfish.tiff']);
        IM = IM_scatter(i_lr,range_im);
        SaveImToTiffStack(IM,tiffdir);
    end

    for i_map = 1:5
        tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_',M_comp_names{i_map},'_anat_allfish.tiff']);
        IM = IM_maps(i_map,range_im);
        SaveImToTiffStack(IM,tiffdir);
    end
end

%% Average Plot
% M_k_scale = {1,1.5,1};
% M_k_contrast = {1.2,1.5,1.2};

for i_set = 1:n_reg
    range_im = range_fish;%M_fishrange_im{i_set};%[1:3,5:7];%[1:3,5:18];

    for i_map = 1:5
        %%
        IM = IM_maps(i_map,range_im);%IM_2(i_lr,range_im);

        % adjust params for visualization
        k_scale = 0.5;%1/1.5;%M_k_scale{i_set};
        k_contrast = 1;%M_k_contrast{i_set};

        [h_anat,im_avr] = AverageAnatPlot(IM,k_contrast,k_scale);

        tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_',M_comp_names{i_map},'_anat_avr.tiff']);
        imwrite(im_avr, tiffdir, 'compression','none','writemode','overwrite');
    end
end

author: xiuyechen

repo: FishExplorer


text:
function runGmsh(geo_name,gmsh_path)
    if isempty(gmsh_path)
        gmsh_path='C:\Users\icberg\OneDrive - UAB - The University of Alabama at Birmingham\Documents\FEM Model Heart Chip';
    end

    oldFolder = cd(gmsh_path);
    system(['gmsh.exe ' geo_name ' -0']);
    cd(oldFolder)
end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%-------------------------------------------------------------------------
% Coupling Kinetic and Fluid System of Euler Equations
% To solve Shock Tube Problem
%
% Kinetic equation use SBBGK and Hydrodynamic equation use Roe Euler
%
% Based on:
%  [1] Pierre Degond, Giacomo Dimarco and Luc Mieussens
%      A moving interface method for dynamic kinetic-fluid coupling
%      Journal of Computation Physics 227(2007)1176-1208
%
% By Manuel Diaz and Yun-Da Tsai
% 007@IAM 25.01.2013
%-------------------------------------------------------------------------

clear all; clc; close all;

%% Global Variables
global CFL r_time theta dt dtdx nx
global w k nv
global gamma etpfix

%% Controling Parameters
name        ='SBBGK1d'; % Simulation Name
CFL         = 0.05;     % CFL condition
r_time      = 1/10000;  % Relaxation time
tEnd        = 0.1;      % End time
theta       = 0;        % {-1} BE, {0} MB, {1} FD.
quad        = 2;        % for NC = 1 , GH = 2
method      = 1;        % for TVD = 1, WENO3 = 2, WENO5 = 3
IC_case     = 1;        % IC: {1}Sod's, {2}LE, {3}RE, {4}DS, {5}SS, {6}Cavitation
plot_figs   = 1;        % 0: no, 1: yes please!
write_ans   = 0;        % 0: no, 1: yes please!
gamma       = 2.4;      % Ratio of specific heats
flxtype     = 2;        % {1} Roe, {2} LF, {3} LLF, {4} Upwind <-non-conservative!
etpfix      = 0.90;     % {#} Harten's sonic entropy fix value, {0} no entropy fix

%% Space Discretization
nx  = 60;                      % number of cells
x   = linspace(0,1,nx);         % Physical domain -x
dx  = max(x(2:end)-x(1:end-1)); % delta x

%% Load Initial Condition
[z0,ux0,t0,p0,rho0,E0] = SSBGK_IC1d(x,IC_case);

%% Load Initial Cut Function
xa = 0;   % buffer left boundary
xb = 1;   % buffer right boundary
%h0  = cutfunc(x,xa,xb);  % Physical Cut Function
%h0 = ones(size(x));
h0 = zeros(size(x));

%% Discretization of the Velocity Space
% Microscopic Velocity Discretization (using Discrete Ordinate Method)
% that is to make coincide discrete values of microscopic velocities with
% values as the value points for using a quadrature method, so that we can
% integrate the velocity probability distribution to recover our
% macroscopics properties.
switch quad

    case{1} % Newton Cotes Quadrature:
    V  = [-20,20];  % range: a to b
    nv = 200;       % nodes desired (may not the actual value)
    [c,w,k] = cotes_xw(V(1),V(2),nv,5); % Using Netwon Cotes Degree 5

    case{2} % Gauss Hermite Quadrature:
    nv = 60;          % nodes desired (the actual value)
    [c,w] = GaussHermite(nv); % for integrating range: -inf to inf
    k = 1;            % quadrature constant.
    w = w.*exp(c.^2); % weighting function of the Gauss-Hermite quadrature

    otherwise
        error('Order must be between 1 and 2');
end

%% Applying DOM
% The actual nv value will be computed using 'lenght' vector function:
nv = length(c);
% Remap velocity points
    c = repmat(c,1,nx);     w = repmat(w,1,nx);
% Remap classical IC
    [rho0,ux0,p0] = apply_DOM(rho0,ux0,p0,nv);
% Remap Semiclassical IC
    [z0,t0,E0] = apply_DOM(z0,t0,E0,nv);
% Remap h coeficient
  h0 = repmat(h0,nv,1);

%% Semi0-classical Equilibrium Distribution Function
M0 = f_equilibrium_1d(z0,ux0,c,t0,theta);

%% Load initial Conditions and Spliting of information
M = M0; rho = rho0; ux = ux0; t = t0; p = p0; z = z0; h = h0; E = E0;

%% Split ICs in R and L
rhor = h.*rho0;     rhol = (1-h).*rho0;
uxr = h.*ux0;       uxl = (1-h).*ux0;
pr = h.*p0;         pl = (1-h).*p0;

[zr,~,tr,~] = macroproperties1d(rhor,rhor.*uxr,pr+0.5*rhor.*uxr.^2,nx,nv,theta);
[zl,~,tl,~] = macroproperties1d(rhol,rhol.*uxl,pl+0.5*rhol.*uxl.^2,nx,nv,theta);

fr = f_equilibrium_1d(zr,uxr,c,tr,theta);

%% Main Loop
% Compute next time step
dt = dx*CFL/max(c(:,1));
time  = 0:dt:tEnd;
dtdx = dt/dx;

Ml = f_equilibrium_1d(zl,uxl,c,tl,theta) ;
Ml(isnan(Ml)) = 0;     fr(isnan(fr)) = 0;
f = fr + Ml; % computed here for ploting purposes

count = 1; % iteration counter
tic;
for tsteps = time

    % Compute M_eq for the entire domain
    M_eq = f_equilibrium_1d(z,ux,c,t,theta);     % total % theta = 0;

    % Break total information into r and l
    if (tsteps>1)
        rhol = (1-h).*rho;
        uxl = (1-h).*ux;
        pl = (1-h).*p;
    end

    % Plot IC
    %if plot_figs == 1; surf(f); end;

    % Compute vector 'q'
    q = [rhol(1,:) ; rhol(1,:).*uxl(1,:) ; pl(1,:)+0.5*rhol(1,:).*uxl(1,:).^2];

    % Update Physical cut function 'h'
    h_next = h; % this means: fixed buffer assumption

% uncomment for BGK ******************************************************
%     % Evaluate Modified Boltzmann BGK
%     [fr_next] = ModSBBGK(h,h_next,M_eq,fr,Ml,c,flxtype);
%
%     % Compute macroscopic moments
%     [rhor,rhour,Er] = macromoments1d(k,w,fr_next,c);
%
%     % macroscopic properties
%     [zr,uxr,tr,pr] = macroproperties1d(rhor,rhour,Er,nx,nv,theta);
%
%         rhol = zeros(size(x));
%         rhoul = zeros(size(x));
%         El = zeros(size(x));
%         Ml_next = zeros(size(fr));
% uncomment for BGK ******************************************************

% uncomment for EULER*****************************************************
    rhor = zeros(size(x));
    rhour = zeros(size(x));
    Er = zeros(size(x));
    fr_next = zeros(size(Ml));

    % Evaluate Modificed Roe Euler Solver
    [rhol,rhoul,El] = ModEuler(h,h_next,q,fr_next,c);

    %plot partial result
    if plot_figs == 1;
        subplot(1,3,1); plot(x,rhol,'o'); title('Density');
        subplot(1,3,2); plot(x,rhoul./rhol,'o'); title('Velocity');
        subplot(1,3,3); plot(x,El,'o'); title('Energy');
    end

    % macroscopic properties
    [zl,uxl,tl,pl] = macroproperties1d(rhol,rhoul,El,nx,nv,theta);

    % Apply DOM
    [tl,zl,uxl] = apply_DOM(tl,zl,uxl,nv);

    % Compute Ml_next
    Ml_next = f_equilibrium_1d(zl,uxl,c,tl,theta) ;

% uncomment for EULER*****************************************************

    % New time step info: sum left and righ values with 'NaN' filter sum
    % function.
    Ml_next(isnan(Ml)) = 0;     fr_next(isnan(fr_next)) = 0;
    f   = fr_next + Ml_next;    % total f
    rho = nansum([rhol;rhor]);  % total density
    rhou= nansum([rhoul;rhour]);% total velocity in x
    E   = nansum([El;Er]);      % total persure

    %Apply Neumann BC's in total variables
    f(:,1)  = f(:,2);           f(:,end)  = f(:,end-1);
    rho(:,1)= rho(:,2);         rho(:,end)= rho(:,end-1);
    rhou(:,1) = rhou(:,2);      rhou(:,end) = rhou(:,end-1);
    E(:,1)  = E(:,2);           E(:,end)  = E(:,end-1);

    % Recover Semiclassical Conditions for next time step
    [z,ux,t,p] = macroproperties1d(rho,rhou,E,nx,nv,theta);

    %Apply DOM
    [rho,ux,p] = apply_DOM(rho,ux,p,nv);
    [t,z,E] = apply_DOM(t,z,E,nv);

    % Update information
    fr = fr_next;

    % update counter
    count = count+1;

    % update plot
    drawnow
end
toc;
% write/plot final output
%plot(x,r_total,'o');

repo: Aero-matlab


text:
function [ret,algo] = training(algo,dat)

  ret=feval(algo.type,algo,dat);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [fx,dF_dX,dF_dTheta] = f_lin1D(Xt,Theta,ut,inF)
% dummy 2D linear evolution function

deltat = inF.delta_t;

try
    a = inF.a.^-1;
catch
    a = 1;
end
try
    [uu,dudtheta] = feval(inF.u_fname,Theta(2:end),ut(2:end),inF);
catch
    uu = 0;
    if size(Theta,1) > 1
        dudtheta = zeros(size(Theta,1)-1,1);
    else
        dudtheta = [];
    end
end

A = -a.*exp(Theta(1));
fx = Xt + deltat.*(A*Xt + ut(1) + uu);
dF_dX = 1 + deltat.*A;
dF_dTheta = deltat.*([Xt*A;zeros(size(Theta,1)-1,1)] + [0;dudtheta]);

author: MBB-team

repo: VBA-toolbox


text:
function tests = test_spaces
  tests = functiontests(localfunctions);
end

function test_orth_complement(testCase)
    M = 100;
    N1 = 10;
    N2 = 10;
    N = N1 + N2;
    X = randn(M, N) ./ sqrt(M);
    A = X(:, 1:N1);
    B = X(:, N1+1:end);
    C = spx.la.spaces.orth_complement(A, B);
    D = A' * C;
    verifyEqual(testCase, D, zeros(size(D)), 'AbsTol', 1e-6);
end

function test_insersection_space(testCase)
    M = 100;
    N = 10;
    X =randn(M, N);
    A = [X(:, 1:4) X(:, 5:6)];
    B = [X(:, 1:4) X(:, 7:8)];
    C = spx.la.spaces.insersection_space(A, B);
    verifyEqual(testCase, rank(C), 4);
end

function test_principal_angles_orth_cos(testCase)
    M = 10;
    N = 10;
    K = N /2 ;
    X = randn(M, N);
    X = orth(X);
    A = X(:, 1:K);
    B = X(:, K+1:N);
    angles = spx.la.spaces.principal_angles_orth_cos(A, B);
    verifyEqual(testCase, angles, zeros(K, 1), 'AbsTol', 1e-6);
end

function test_two_spaces_at_angle(testCase)
    N = 50;
    theta = pi/4;
    [A, B] = spx.data.synthetic.subspaces.two_spaces_at_angle(N, theta);
    verifyTrue(testCase, spx.matrix.is_orthonormal(A));
    verifyTrue(testCase, spx.matrix.is_orthonormal(B));
    phi = spx.la.spaces.smallest_angle_rad(A, B);
    verifyEqual (testCase, phi, theta, 'AbsTol', 1e-12);
    thetas = 0.1:0.1:pi/2;
    for i=1:numel(thetas)
        theta = thetas(i);
        [A, B] = spx.data.synthetic.subspaces.two_spaces_at_angle(N, theta);
        verifyTrue(testCase, spx.matrix.is_orthonormal(A));
        verifyTrue(testCase, spx.matrix.is_orthonormal(B));
        phi = spx.la.spaces.smallest_angle_rad(A, B);
        verifyEqual (testCase, phi, theta, 'AbsTol', 1e-12);
    end
end

function test_have_same_column_spans(testCase)
    A = randn(6, 3);
    B = orth(A);
    verifyTrue(testCase, spx.la.spaces.have_same_column_spans(A, B));
    C = randn(6, 3);
    verifyFalse(testCase, spx.la.spaces.have_same_column_spans(A, C));
end

function test_three_disjoint_spaces_at_angle(testCase)
    N = 50;
    theta = pi/6;
    [A, B, C] = spx.la.spaces.three_disjoint_spaces_at_angle(N, theta);
    spx.la.spaces.describe_three_spaces(A, B, C);
    verifyTrue(testCase, spx.matrix.is_orthonormal(A));
    verifyTrue(testCase, spx.matrix.is_orthonormal(B));
    verifyTrue(testCase, spx.matrix.is_orthonormal(C));
    phi = spx.la.spaces.smallest_angle_rad(A, B);
    verifyEqual (testCase, phi, theta, 'AbsTol', 1e-12);
    phi = spx.la.spaces.smallest_angle_rad(B, C);
    verifyEqual (testCase, phi, theta, 'AbsTol', 1e-12);
    phi = spx.la.spaces.smallest_angle_rad(A, C);
    verifyEqual (testCase, phi,2*theta, 'AbsTol', 1e-12);
end

author: indigits

repo: sparse-plex


text:
function ChanceConstraint = horzcat(varargin)

if nargin == 2
    ChanceConstraint = varargin{1};
    for i = 1:length(varargin{2}.P)
        ChanceConstraint.P{end+1} = varargin{2}.P{i};
        ChanceConstraint.level{end+1} = varargin{2}.level{i};
    end
else
    ChanceConstraint = horzcat(varargin{1},horzcat(varargin{2:end}));
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function ffwd(H)
%FFWD Fast forward movie in movie player

feval(H.fcns.ffwd,[],[],H.hfig);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function measurements = posegraphMeasurementMatrix(graph)
if(strcmpi(graph.format,'2d'))
    n_edges = size(graph.edges,1);
    measurements = arrayfun(@(x) struct('i',[],'j',[], ...
        'R', [], 't', [], 'tau', [], 'kappa', []), 1:n_edges);
    for i = 1:n_edges
        edge = graph.edges(i, :);
        R_ij = rot2D(edge(5));

        Omega = [edge(6) edge(7) edge(8); ...
            edge(7) edge(9) edge(10); ...
            edge(8) edge(10) edge(11)];
        Omega = [R_ij zeros(2,1); 0 0 1]' * Omega * [R_ij zeros(2,1); 0 0 1];
        if(min(eig(Omega)) < 0)
            warning('Error in full covariance construction')
            disp(min(eig(Omega)));
        end

        measurements(i).i = edge(1);
        measurements(i).j = edge(2);
        measurements(i).t = edge(3:4)';
        measurements(i).R = R_ij;
        measurements(i).Omega = Omega;
        measurements(i).tau = 2/trace(inv(Omega(1:2,1:2)));
        measurements(i).kappa = Omega(3,3);
        measurements(i).weight = 1;
    end
elseif(strcmpi(graph.format,'3d'))
    n_edges = size(graph.edges,2);
    measurements = arrayfun(@(x) struct('i',[],'j',[], ...
        'R', [], 't', [], 'tau', [], 'kappa', []), 1:n_edges);
    for i = 1:n_edges
        edge = graph.edges(i);
        measurements(i).i = edge.i;
        measurements(i).j = edge.j;
        measurements(i).t = edge.t;
        measurements(i).R = edge.R;
        measurements(i).Omega = edge.Info;
        measurements(i).tau = 3/trace(inv(edge.Info(1:3,1:3)));
        measurements(i).kappa = 3/(2*trace(inv(edge.Info(end-2:end,end-2:end))));
        measurements(i).weight = 1;
    end
else
    error('Unknown format...')
end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
% The mapping from 61 phones set to 48 phonse set for TIMIT. The mapping
% table is as follows:
%
% cl vcl epi --> sil
% el --> l
% en --> n
% zh --> sh
% aa --> ao
% ix --> ih
% ax --> ah
%
% Author: Xiong Xiao, NTU, Singapore
% Date: 1 Jun 2016

function phone_seq_39 = TIMIT_map48to39(phone_seq_48)

mapping_table = {
    'cl'  'sil';
    'vcl' 'sil';
    'epi' 'sil';
    'el'  'l';
    'en'  'n';
    'zh'  'sh';
    'aa'  'ao';
    'ix'  'ih';
    'ax'  'ah';};

phone_seq_39 = phone_seq_48;
for i=1:size(mapping_table,1)
    found = strcmpi(phone_seq_39, mapping_table{i,1});
    if sum(found)>0
        idx = find(found==1);
        for j=1:length(idx)
            phone_seq_39{idx(j)} = mapping_table{i,2};
        end
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function varargout = wave_clus_aux4(varargin)
% WAVE_CLUS_AUX1 M-file for wave_clus_aux.fig
%      WAVE_CLUS_AUX, by itself, creates a new WAVE_CLUS_AUX or raises the existing
%      singleton*.
%
%      H = WAVE_CLUS_AUX returns the handle to a new WAVE_CLUS_AUX or the handle to
%      the existing singleton*.
%
%      WAVE_CLUS_AUX('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in WAVE_CLUS_AUX.M with the given input arguments.
%
%      WAVE_CLUS_AUX('Property','Value',...) creates a new WAVE_CLUS_AUX or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before wave_clus_aux_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to wave_clus_aux_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help wave_clus_aux

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @wave_clus_aux4_OpeningFcn, ...
                   'gui_OutputFcn',  @wave_clus_aux4_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before wave_clus_aux is made visible.
function wave_clus_aux4_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to wave_clus_aux (see VARARGIN)

% Choose default command line output for wave_clus_aux
handles.output = hObject;

set(handles.isi24_accept_button,'value',1);
set(handles.isi25_accept_button,'value',1);
set(handles.isi26_accept_button,'value',1);
set(handles.isi27_accept_button,'value',1);
set(handles.isi28_accept_button,'value',1);
set(handles.fix24_button,'value',0);
set(handles.fix25_button,'value',0);
set(handles.fix26_button,'value',0);
set(handles.fix27_button,'value',0);
set(handles.fix28_button,'value',0);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes wave_clus_aux wait for user response (see UIRESUME)
% uiwait(handles.wave_clus_aux);

% --- Outputs from this function are returned to the command line.
function varargout = wave_clus_aux4_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

h_figs=get(0,'children');
h_fig = findobj(h_figs,'tag','wave_clus_figure');
USER_DATA = get(h_fig,'UserData');
par = USER_DATA{1};

for i = 24:28
	si = num2str(i);
	set(eval(['handles.isi' si '_accept_button']),'value',1);
	set(eval(['handles.isi' si '_reject_button']),'value',0);
	set(eval(['handles.fix' si '_button']),'value',0);

	eval(['set(handles.isi' si '_nbins,''string'',par.nbins' si ');']);
	eval(['set(handles.isi' si '_bin_step,''string'',par.bin_step' si ');']);

end

plot_spikes_aux(handles,USER_DATA,4);

author: csn-le


text:
function [roc_avg, pr_avg] = DTINet(seed, nFold, interaction, drug_feat, prot_feat, dim_imc)
	rng(seed);
	Pint = find(interaction); % pair of interaction
	Nint = length(Pint);
	Pnoint = find(~interaction);
	Pnoint = Pnoint(randperm(length(Pnoint), Nint * 1));
	Nnoint = length(Pnoint);

	posFilt = crossvalind('Kfold', Nint, nFold);
	negFilt = crossvalind('Kfold', Nnoint, nFold);

	AUROC = zeros(nFold, 1);
	AUPRC = zeros(nFold, 1);

	for foldID = 1 : nFold
		train_posIdx = Pint(posFilt ~= foldID);
		train_negIdx = Pnoint(negFilt ~= foldID);
		train_idx = [train_posIdx; train_negIdx];
		Ytrain = [ones(length(train_posIdx), 1); zeros(length(train_negIdx), 1)];
		fprintf('Train data: %d positives, %d negatives\n', sum(Ytrain == 1), sum(Ytrain == 0));

		test_posIdx = Pint(posFilt == foldID);
		test_negIdx = Pnoint(negFilt == foldID);
		test_idx = [test_posIdx; test_negIdx];
		Ytest = [ones(length(test_posIdx), 1); zeros(length(test_negIdx), 1)];
		fprintf('Test data: %d positives, %d negatives\n', sum(Ytest == 1), sum(Ytest == 0));

		[I, J] = ind2sub(size(interaction), train_idx);
		Xtrain = sparse(I, J, Ytrain, size(interaction, 1), size(interaction, 2));

		[W, H, ~] = train_mf(Xtrain, sparse(drug_feat), sparse(prot_feat), ...
						[' -l ' num2str(1) ' -k ' num2str(dim_imc) ' -t 10' ' -s ' num2str(10)]);
		Zscore = drug_feat * W' * H * prot_feat';
		Ypred = Zscore(test_idx);

		[trainroc, trainpr] = auc(Ytrain, Zscore(train_idx), 1e-6);
		[testroc, testpr] = auc(Ytest, Ypred, 1e-6);
		AUROC(foldID) = testroc;
		AUPRC(foldID) = testpr;
		fprintf('Fold %d, Train: AUROC=%f AUPR=%f; Test: AUROC=%f, AUPR=%f\n', foldID, trainroc, trainpr, testroc, testpr);
	end
	roc_avg = mean(AUROC);
    pr_avg = mean(AUPRC);
end

author: luoyunan

repo: DTINet

save_path: github-repos/MATLAB/luoyunan-DTINet


text:
function accuracy = ldrc(TrainSet, TestSet, train_num, test_num, class_num, reduce_dimension, options)
% Linear discriminant regression classificatoin (LDRC) algorithm
%
% Inputs:
%       TrainSet            train sets of size dxn, where d is dimension and n is number of sets
%       TestSet             test sets of size dxn, where d is dimension and n is number of sets
%       test_num            numner of test sets
%       class_num           numner of classes
%       options             options
% Output:
%       accuracy            classification accurary
%
% References:
%       S.-M. Huang and J.-F. Yang,
%       "Linear discriminant regression classification for face recognition,"
%       IEEE Signal Processing Letters, vol.20, no.1, pp.91-94, 2013.
%
%
% Created by H.Kasai on July 05, 2017

    % extract options
    if ~isfield(options, 'verbose')
        verbose = false;
    else
        verbose = options.verbose;
    end

    if verbose
        fprintf('# LDRC: calculate projection matrix U ...');
    end
    Eb = 0;
    Ew = 0;
    for j = 1 : train_num

        j_class = TrainSet.y(1, j);

        % calcuate Ew
        intra_class_index = find(TrainSet.y == j_class);
        intra_class_index(intra_class_index == j) = [];
        X_intra_wo_j = TrainSet.X(:, intra_class_index);        % X^{intra}
        for k = 1 : size(X_intra_wo_j, 2)
            error_intra = TrainSet.X(:,j) - X_intra_wo_j(:,k);
            Ew = Ew + error_intra * error_intra';               % Eq.(14)
        end

        % calcuate Eb
        inter_class_index = find(TrainSet.y ~= j_class);
        X_inter = TrainSet.X(:, inter_class_index);             % X^{inter}

        for k = 1 : size(X_inter, 2)
            error_inter = TrainSet.X(:,j) - X_inter(:,k);
        end
        Eb = Eb + error_inter * error_inter';                   % Eq.(13)
    end

    % calc
    Ew = Ew/train_num;                                          % Eq.(14)
    Eb = Eb/(train_num * (class_num-1));                        % Eq.(13)
    epsilon = 0.001;
    Ew = Ew + epsilon*eye(size(Ew,2));

    [V, ~] = eig(Eb, Ew);                                       % Eq.(18)
    for i = 1:reduce_dimension
        U(:,i) = V(:,size(V,2)+1-i);
    end
    if verbose
        fprintf('done\n');
    end

    % reduce dimention
    TrainSet.X_red  = U' * TrainSet.X;
    TestSet.X_red   = U' * TestSet.X;

    % prepare projection matrix (hat matrix)
    H = cell(1, class_num);
    for i = 1 : class_num
        class_index = find(TrainSet.y == i);

        X_i = TrainSet.X_red(:, class_index);
        alpha_i = pinv(X_i' * X_i) *  X_i';
        H{i} = X_i * alpha_i;

        if verbose
            fprintf('# LDRC: calc Hi for class : %03d/%03d (samples: %03d)\n', i, class_num, length(class_index));
        end
    end

    % prepare predicted label array
    identity = zeros(1, test_num);

    % predict the class
    for j = 1 : test_num

        dis = zeros(1, class_num);
        y = TestSet.X_red(:, j);
        for i = 1 : class_num
            y_pred = H{i} * y;
            dis(1, i) = norm(y - y_pred);
        end
        [~, label] = min(dis);
        identity(j) = label;

        if verbose
            correct = (label == TestSet.y(1, j));
            fprintf('# LDRC: test:%03d, predict class: %03d --> ground truth :%03d (%d)\n', j, label, TestSet.y(1, j), correct);
        end
    end

    % calculate accuracy
    correct_num = sum(identity == TestSet.y);
    accuracy = correct_num / test_num;
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function linpack_z_test14 ( )

%*****************************************************************************80
%
%% TEST14 tests ZHPCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST14\n' );
  fprintf ( 1, '  For a double precision complex (C)\n' );
  fprintf ( 1, '  Hermitian matrix using packed storage (HP),\n' );
  fprintf ( 1, '  ZHPCO factors the matrix and estimates\n' );
  fprintf ( 1, '  the reciprocal condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  k = 0;
  seed = 123456789;

  for j = 1 : n

    for i = 1 : j-1
      k = k + 1;
      [ a(k), seed ] = c8_uniform_01 ( seed );
      a_save(i,j) = a(k);
      a_save(j,i) = conj ( a(k) );
    end

    k = k + 1;
    [ a(k), seed ] = r8_uniform_01 ( seed );
    a_save(j,j) = a(k);

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a_save(i,j) ), imag ( a_save(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, rcond ] = zhpco ( a, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition RCOND = %f\n', rcond );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [outputs] = ml_clustering(X,options,varargin)
%ML_CLUSTERING Applies clustering technique to the data
%
%   input -----------------------------------------------------------------
%
%       o X : (N x D), a data set with N samples each being of dimension D.
%
%       o options : structure
%
%               options.method_name  = 'kmeans' : this options should be given
%
%  output ----------------------------------------------------------------
%
%

%% Default options

outputs = [];
K       = 1;  % default number of clusters is assumed to be 1.
kernel  = 'gauss';
kpar    = 1;

%% Check input

if ~isfield(options,'method_name'),  error('field of method_name of structure options was not defined!'); end
if isfield(options,'K'),             K = options.K;end
if isfield(options,'kernel')         [kernel, kpar] = sanitise_kernel_input(options);end

switch options.method_name

    case 'kmeans'

        if isempty(varargin)
            [idx, centroids]    = kmeans(X, K);
        else
            [idx, centroids]    = kmeans(X, K,varargin{:});
        end

        pnames                  =  {   'distance'  'start' 'replicates' 'emptyaction' 'onlinephase' 'options' 'maxiter' 'display'};
        dflts                   =  {'sqeuclidean' 'plus'          []  'singleton'         'on'        []        []        []};
        distance                =  internal.stats.parseArgs(pnames, dflts, varargin{:});
        distNames               = {'sqeuclidean','cityblock','cosine','correlation','hamming'};
        distance                =  internal.stats.getParamVal(distance,distNames,'''Distance''');

        outputs.distance        = distance;
        outputs.K               = K;
        outputs.method_name     = 'kmeans';
        outputs.labels          = idx;
        outputs.centroids       = centroids;

    case 'kernel-kmeans'

        % Old implementation, not working properly
%        Kn                       = gram(X, X, kernel,kpar(1),kpar(2));
%        [idx,centroids,eigens]   = kernelkmeans(Kn, K);
%        outputs.K               = K;
%        outputs.method_name     = 'kernel-kmeans';
%        outputs.labels          = idx;
%        outputs.centroids       = centroids;
%        outputs.eigens          = eigens;
%        outputs.kernel          = kernel;
%        outputs.kpar            = kpar;

       % New implementation, working! But no decision boundaries or
       % isolines .. not yet at least
       [M, ~] = size(X);
       init = ceil(K*rand(1,M));
       % Using new implementation
       if strcmp(kernel,'gauss')
            % Gaussian Kernel
            [y,model,mse] = knKmeans(X',init,@knGauss, kpar(1));
       elseif strcmp(kernel,'poly')
            % Gaussian Kernel
            [y,model,mse] = knKmeans(X',init,@knPoly, kpar(2), kpar(1));
       end
       outputs.K               = K;
       outputs.eigens          = model.eigens;
       outputs.lambda          = diag(model.lambda);
       outputs.centroids       = model.centroids;
       outputs.method_name     = 'kernel-kmeans';
       outputs.mse             = mse;
       outputs.labels          = y;
       outputs.kernel          = kernel;
       outputs.kpar            = kpar;

end

author: epfl-lasa


text:
function tests = test_stackSlice(varargin)
%TEST_STACKSLICE  One-line description here, please.
%
%   output = test_stackSlice(input)
%
%   Example
%   test_stackSlice
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-12-02,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function test_sliceX_gray(testCase) %#ok<*DEFNU>

img = createTestImage;
dim = stackSize(img);

sliceYZ = stackSlice(img, 1, 5);
assertEqual(testCase, [dim(3) dim(2)], size(sliceYZ));

sliceYZ = stackSlice(img, 'x', 5);
assertEqual(testCase, [dim(3) dim(2)], size(sliceYZ));

function test_sliceY_gray(testCase) %#ok<*DEFNU>

img = createTestImage;
dim = stackSize(img);

sliceZX = stackSlice(img, 2, 5);
assertEqual(testCase, [dim(1) dim(3)], size(sliceZX));

sliceZX = stackSlice(img, 'y', 5);
assertEqual(testCase, [dim(1) dim(3)], size(sliceZX));

function test_sliceZ_gray(testCase) %#ok<*DEFNU>

img = createTestImage;
dim = stackSize(img);

sliceXY = stackSlice(img, 3, 5);
assertEqual(testCase, [dim(2) dim(1)], size(sliceXY));

sliceXY = stackSlice(img, 'z', 5);
assertEqual(testCase, [dim(2) dim(1)], size(sliceXY));

function img = createTestImage

[x, y, z] = meshgrid(1:10, 1:15, 1:20);
img = 5*x + 4*y + 3*z;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [varargout] = rfbevent(varargin)

% RFBEVENT sends a keyboard or mouse event to a VNC server
%
% RFB ("remote frame buffer") is a simple protocol for remote access to
% graphical user interfaces. Because it works at the framebuffer level it
% is applicable to all windowing systems and applications, including X11,
% Windows and Macintosh. RFB is the protocol used in VNC (Virtual Network
% Computing).
%
% The remote endpoint where the user sits (i.e. the display plus keyboard
% and/or pointer) is called the RFB client or viewer. The endpoint where
% changes to the framebuffer originate (i.e. the windowing system and
% applications) is known as the RFB server.
%
% Use as
%   rfbevent(display, passwd, eventtype, eventvalue, ...)
%
% Some examples
%   rfbevent('vncserver:5901', 'yourpasswd', 'Text',   'xclock')          % type multiple characters
%   rfbevent('vncserver:5901', 'yourpasswd', 'Button', 'Return')          % single key event, press and release
%   rfbevent('vncserver:5901', 'yourpasswd', 'Button', 'Return',  0)      % single key event, press and release
%   rfbevent('vncserver:5901', 'yourpasswd', 'Button', 'Return',  1)      % single key event, press only
%   rfbevent('vncserver:5901', 'yourpasswd', 'Button', 'Return', -1)      % single key event, release only
%   rfbevent('vncserver:5901', 'yourpasswd', 'Pointer', [20 100])         % only mouse position
%   rfbevent('vncserver:5901', 'yourpasswd', 'Pointer', [20 100 1])       % mouse position and button 1, press and release
%   rfbevent('vncserver:5901', 'yourpasswd', 'Pointer', [20 100 1],  0)   % mouse position and button 1, press and release
%   rfbevent('vncserver:5901', 'yourpasswd', 'Pointer', [20 100 1],  1)   % mouse position and button 1, press only
%   rfbevent('vncserver:5901', 'yourpasswd', 'Pointer', [20 100 1], -1)   % mouse position and button 1, release only
%
% Note that the password has to be represented as plain text in the matlab
% script/function that is using RFBEVENT, which poses a potential security
% problem. The password is sent over the network to the VNC server after
% being encrypted.
%
% This implements the KeyEvent and PointerEvent messages according to
% "The RFB Protocol" by Tristan Richardson (RealVNC Ltd, formerly of
% Olivetti Research Ltd / AT&T Labs Cambridge) Version 3.8 (Last updated 8
% June 2007), http://www.realvnc.com/docs/rfbproto.pdf

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

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);

  % the following is specific for this particular mex file
  mex rfbevent.c d3des.c

  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was just created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end

author: fieldtrip

repo: fieldtrip


text:
function [score,in,out]=clusterquality(f,S,partition)
%function [score,in,out]=clusterquality(f,S,partition)
%
%PURPOSE
%
%To compute a simple quality (compactness) index for the given
%clusters. The more dense and isolated a cluster is the bigger
%value the index gets; however, if there is only one item in a
%cluster the index is not computed (a NaN is returned for that
%cluster)
%
%INPUT
%
% f          (string)     cluster score function: 'mean' or 'minmax'
% S          (NxN matrix) similarity matrix
% partition  (Nx1 vector) partition vector (see explanation for
%              'partition vector' in function hcluster)
%
%OUTPUT
%
%score       (Kx1 vector) score(k) contains the index value for
%             cluster k: score(k)=in(k)-out(k)
%in          (Kx1 vector)  see above
%out         (Kx1 vector)  see above
%
%DETAILS
%
%For f='mean' the average extra-cluster similarity is subtracted
%from the mean intra-cluster similarity.
%
%   Iq=avg(intra-cluster similarity) - avg(extra-cluster similarity)
%
%The index is described in publication Himberg et al. (2004),
%"Validating the independent components of neuroimaging time-series
%via clustering and visualization". NeuroImage, 22:3(1214-1222). It
%is used by default in Icasso functions.
%
%(The intra-cluster similarities for cluster C mean the mutual
%similarities between the estimates in C, and the extra-cluster
%similarities for C are the similarities between the estimates in C
%and the estimates not in C.)
%
%Note: if there is only one item in a cluster, index
%value NaN is returned for that cluster.
%
%You can also specify a more conservative but less robust criteria
%('minmax') where the maximum extra-cluster similarity is
%subtracted from the minimum intra-cluster similarity: if
%f='minmax' the maximum extra-cluster similarity is subtracted from
%the minimum intra-cluster similarity.
%
%SEE ALSO
% clusterstat
% icassoStability

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 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.

% ver 1.2 100105

Ncluster=max(partition);
s=clusterstat(S,partition);

% Compute score
switch lower(f)
 case 'mean'
  in=s.internal.avg;
  out=s.external.avg;
 case 'minmax'
  in=s.internal.min;
  out=s.external.max;
 otherwise
  error('Unrecognized score function.');
end

score=in-out;


text:
% EEG_INSERTBOUND - insert boundary event in an EEG event structure.
%
% Usage:
%       >> [eventout indnew] = eeg_insertbound( eventin, pnts, ...
%                                                abslatency, duration);
% Required Inputs:
%   eventin    - EEGLAB event structure (EEG.event)
%   pnts       - data points in EEG dataset (EEG.pnts * EEG.trials)
%   abslatency - absolute latency of regions in original dataset. Can
%                also be an array of [beg end] latencies with one row
%                per region removed. Then 'lengths' argument is ignored.
% Optional Inputs:
%   lengths    - lengths of removed regions
%
% Outputs:
%   eventout   - EEGLAB event output structure with added boundaries
%   indnew     - array of indices returning new event index for any old
%                (input eventin) event index
% Notes:
%   This function performs the following:
%   1) add boundary events to the 'event' structure;
%        remove nested boundaries;
%        recompute the latencies of all events.
%   2) all latencies are given in (float) data points.
%        e.g., a latency of 2000.3 means 0.3 samples (at EEG.srate)
%              after the 2001st data frame (since first frame has latency 0).
%
% Author: Arnaud Delorme and Hilit Serby, SCCN, INC, UCSD, April, 19, 2004
%
% See also: EEG_EEGREJ, POP_MERGESET

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

function [eventout,indnew] = eeg_insertboundold( eventin, pnts, regions, lengths)

    if nargin < 3
        help eeg_insertbound;
        return;
    end

    if size(regions,2) ~= 1 && exist('lengths') ~= 1
        lengths = regions(:,2)-regions(:,1)+1;
        regions = regions(:,1);
    end
    if exist('lengths') ~= 1
        lengths = zeros(size(regions));
    end

    % resort events
    if ~isempty(eventin)
        [~,sortedEvents] = sort([eventin.latency]);
        eventin = eventin(sortedEvents);
    end

    if ~isempty(regions)
        fprintf('eeg_insertbound(): %d boundary (break) events added.\n', size(regions, 1));
    else
        return;
    end

    % recompute latencies of boundevents (in new dataset)
    % ---------------------------------------------------
    [regions, tmpsort] = sort(regions);
    lengths            = lengths(tmpsort);
    boundevents        = regions(:,1)-0.5;

    % sort boundevents by decreasing order (otherwise bug in new event index)
    % ------------------------------------
    boundevents = boundevents(end:-1:1);
    lengths     = lengths    (end:-1:1);
    eventout    = eventin;
    indnew      = 1:length(eventin);
    countrm     = 0;
	for tmpindex = 1:length(boundevents) % sorted in decreasing order
        if boundevents(tmpindex) >= 0.5 && boundevents(tmpindex) <= pnts

            % find event succeeding boundary to insert event
            % at the correct location in the event structure
            % ----------------------------------------------
            if ~isempty(eventout) && isfield(eventout, 'latency')
                alllats   = [ eventout.latency ] - boundevents(tmpindex);
                tmpind    = find( alllats >= 0 );
                [~, tmpind2 ] = min(alllats(tmpind));
                tmpind2       = tmpind(tmpind2);
            else
                tmpind2 = [];
            end

            % insert event at tmpind2
            % -----------------------
            if ~isempty(tmpind2)
                eventout(end+1).type = eeg_boundarytype(eventout);
                tmp = eventout(end);
                eventout(tmpind2+1:end) = eventout(tmpind2:end-1);
                eventout(tmpind2) = tmp;
                indnew(tmpind2:end) = indnew(tmpind2:end)+1;
            else
                tmpind2 = length(eventout)+1;
                eventout(tmpind2).type = eeg_boundarytype(eventout);
            end
            eventout(tmpind2).latency  = boundevents(tmpindex);
            eventout(tmpind2).duration = lengths(tmpindex); % just to create field

            [ tmpnest, addlength ] = findnested(eventout, tmpind2);

            % recompute latencies and remove events in the rejected region
            % ------------------------------------------------------------
            eventout(tmpnest) = [];
            countrm           = countrm+length(tmpnest);
            for latind = tmpind2+1:length(eventout)
                eventout(latind).latency = eventout(latind).latency-lengths(tmpindex);
            end

            % add lengths of previous events (must be done after above)
            % ---------------------------------------------------------
            eventout(tmpind2).duration = lengths(tmpindex)+addlength;
            if eventout(tmpind2).duration == 0, eventout(tmpind2).duration=NaN; end

        end
	end

    if countrm > 0
        fprintf('eeg_insertbound(): event latencies recomputed and %d events removed.\n', countrm);
    end

% look for nested events
% return indices of nested events and
% their total length
% -----------------------------------
function [ indnested, addlen ] = findnested(event, ind)
    indnested = [];
    addlen = 0;
    tmpind = ind+1;

    while tmpind <= length(event) && ...
        event(tmpind).latency < event(ind).latency+event(ind).duration
        if eeg_isboundary(event(tmpind))
            if ~isempty( event(tmpind).duration )
                addlen    = addlen + event(tmpind).duration;
                % otherwise old event duration or merge data discontinuity
            end
        end
        indnested = [ indnested tmpind ];
        tmpind = tmpind+1;
    end

% remove urevent and recompute indices
% THIS FUNCTION IS DEPRECATED
% ------------------------------------
function [event, urevent] = removenested(event, urevent, nestind)

    if ~isempty(nestind)
        fprintf('eeg_insertbound() debug msg: removing %d nested urevents\n', length(nestind));
        nestind = sort(nestind);
        urind = [ event.urevent ]; % this must not be done in the loop
                                   % since the indices are dynamically updated
    end

    for ind = 1:length(nestind)
        % find event urindices higher than the urevent to suppress
        % --------------------------------------------------------
        tmpind = find( urind > nestind(ind) );
        for indevent = tmpind
            event(indevent).urevent = event(indevent).urevent-1;
        end
    end

    urevent(nestind) = [];

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
  function [wt,info] = wfbtinit(wtdef,varargin)
%WFBTINIT Initialize Filterbank Tree
%   Usage:  wt = wfbtinit(wtdef);
%
%   Input parameters:
%         wtdef : Filterbank tree definition.
%
%   Output parameters:
%         wt    : Structure describing the filter tree.
%
%   `wfbtinit({w,J,flag})` creates a filterbank tree of depth *J*. The
%   parameter *w* defines a basic wavelet filterbank. For all possible
%   formats see |fwt|.  The following optional flags (still inside of the
%   cell-array) are recognized:
%
%   'dwt','full','doubleband','quadband','octaband'
%     Type of the tree to be created.
%
%   `wfbtinit({w,J,flag,'mod',mod})` creates a filterbank tree as before,
%   but modified according to the value of `mod`.
%   Recognized options:
%
%   `'powshiftable'`
%      Changes subsampling factors of the root to 1. This results in redundant
%      near-shift invariant representation.
%
%   The returned structure *wt* has the following fields:
%
%   .nodes
%      Cell-array of structures obtained from |fwtinit|. Each element
%      define a basic wavelet filterbank.
%
%   .children
%      Indexes of children nodes
%
%   .parents
%      Indexes of a parent node
%
%   .forder
%      Frequency ordering of the resultant frequency bands.
%
%   The structure together with functions from the `wfbtmanip`
%   subdirectory acts as an abstract data structure tree.
%
%   Regular `wfbtinit` flags:
%
%   'freq','nat'
%     Frequency or natural ordering of the coefficient subbands. The direct
%     usage of the wavelet tree (`'nat'` option) does not produce coefficient
%     subbans ordered according to the frequency. To achieve that, some
%     filter shuffling has to be done (`'freq'` option).
%
%   See also: wfbtput, wfbtremove

% AUTHOR: Zdenek Prusa

% TO DO: Do some caching

% Output structure definition.
% Effectively, it describes a ADT tree.
% .nodes, .children, .parents ale all arrays of the same length and the j-th
% node in the tree is desribed by wt.nodes{j}, wt.children{j} and
% wt.parents(j). wt.nodes{j} is the actual data stored in the tree node
% (a structure returned form fwtinit) and wt.parents(j) and wt.children{j}
% define relationship to other nodes. wt.parents(j) is an index of the
% parent in the arrays. wt.children{j} is an array of indexes of the
% children nodes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wt.nodes = {};
wt.children = {};
wt.parents = [];
wt.freqOrder = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
info.istight = 0;

% return empty struct if no argument was passed
if(nargin<1)
    return;
end

do_strict = 0;
do_dual = 0;

% Check 'strict'
if iscell(wtdef) && ischar(wtdef{1}) && strcmpi(wtdef{1},'strict')
   do_strict = 1;
   wtdef = wtdef{2:end};
end
% Check 'dual'
if iscell(wtdef) && ischar(wtdef{1}) && strcmpi(wtdef{1},'dual')
   do_dual = 1;
   wtdef = wtdef{2:end};
end

definput.import = {'wfbtcommon'};
[flags,kv]=ltfatarghelper({},definput,varargin);

% If wtdef is already a structure
if (isstruct(wtdef)&&isfield(wtdef,'nodes'))
    if isempty(wtdef.nodes)
        error('%s: The tree struct is empty.',upper(mfilename));
    end

    wt = wtdef;

    if do_dual || do_strict
       nodesArg = wt.nodes;
       if do_dual
          nodesArg = cellfun(@(nEl) {'dual',nEl},nodesArg,'UniformOutput',0);
       end
       if do_strict
          nodesArg = cellfun(@(nEl) {'strict',nEl},nodesArg,'UniformOutput',0);
       end
       info.istight = 1;
       wt.nodes = {};
       for ii=1:numel(nodesArg)
        % wt.nodes = cellfun(@(nEl) fwtinit(nEl),nodesArg,'UniformOutput',0);
          [wt.nodes{ii},infotmp] = fwtinit(nodesArg{ii});
          if info.istight
              info.istight = infotmp.istight;
          end
       end
       % Do the filter frequency shuffling again, since the filters were
       % overwritten in fwtinit.
       if wt.freqOrder
          wt = nat2freqOrder(wt);
       end
    end

    % Do filter shuffling if flags.do_freq differs from the wt.freqOrder.
    % Frequency and natural oreding coincide for DWT.
    if wt.freqOrder ~= flags.do_freq
       wt = nat2freqOrder(wt);
       wt.freqOrder = ~wt.freqOrder;
    end
    return;
end

% break if the input parameter is not in the correct format
if ~(iscell(wtdef)) || isempty(wtdef)
    error('%s: Unsupported filterbank tree definition.',upper(mfilename));
end

% Creating new tree
% Now wtdef is this {w,J,flag}
wdef = wtdef{1};
definput = [];
definput.flags.treetype = {'full','dwt','doubleband','quadband',...
                           'octaband','root'};
definput.keyvals.mod = [];
definput.keyvals.overcomplete = [];
definput.keyvals.J = [];
[flags2,kv2,J]=ltfatarghelper({'J'},definput,wtdef(2:end));

complainif_notposint(J,'J');

if do_dual
   wdef = {'dual',wdef};
end

if do_strict
   wdef = {'strict',wdef};
end

do_powshiftable = 0;
% Is the first level filterbank different?
if ~isempty(kv2.mod)
    if ischar(kv2.mod)
        if strcmpi(kv2.mod,'powshiftable')
           if ~flags2.do_dwt
              error('%s: powshiftable is only valid with the dwt flag.',...
                    upper(mfilename));
           end
           do_powshiftable = 1;
        else
           error('%s: Not recognized value for the mod key.',upper(mfilename));
        end
    else
        error('%s: Not recognized value for the first key.',upper(mfilename));
    end
end

if ~isempty(kv2.overcomplete)
     error('%s: TO DO: overcomplete.',upper(mfilename));
end

[w, info] = fwtinit(wdef);

% Doing one-node tree
if flags2.do_root
   J = 1;
end

if flags2.do_dwt
   % fill the structure to represent a DWT tree
   for jj=0:J-1
      wt = wfbtput(jj,0,w,wt);
   end
elseif flags2.do_full
   % fill the structure to represent a full wavelet tree
   for jj=0:J-1
     % for ii=0:numel(w.g)^(jj)-1
         wt = wfbtput(jj,0:numel(w.g)^(jj)-1,w,wt);
     % end
   end
elseif flags2.do_doubleband
   % fill the structure to represent a double band tree
   for jj=0:J-1
     % for ii=0:numel(w.g)^(jj)-1
         wt = wfbtput(2*jj,0,w,wt);
         wt = wfbtput(2*jj+1,0:1,w,wt);
     % end
   end

elseif flags2.do_quadband
   % fill the structure to represent a quad band tree
   for jj=0:J-1
     % for ii=0:numel(w.g)^(jj)-1
         wt = wfbtput(3*jj,0,w,wt);
         wt = wfbtput(3*jj+1,0:1,w,wt);
         wt = wfbtput(3*jj+2,0:3,w,wt);
     % end
   end
elseif flags2.do_octaband
   % fill the structure to represent a octa band tree
   for jj=0:J-1
     % for ii=0:numel(w.g)^(jj)-1
         wt = wfbtput(4*jj,0,w,wt);
         wt = wfbtput(4*jj+1,0:1,w,wt);
         wt = wfbtput(4*jj+2,0:3,w,wt);
         wt = wfbtput(4*jj+3,0:7,w,wt);
     % end
   end
end

% Do filter shuffling if frequency ordering is required,
wt.freqOrder = flags.do_freq;
if flags.do_freq
   wt = nat2freqOrder(wt);
end

if do_powshiftable
    % Change subsampling factors of the root to 1
    wt.nodes{wt.parents==0}.a(:) = 1;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function xfinal = my_ola(xi, len, len1)
len2 = len - len1;

nFr = size(xi,2);
k=1;
xfinal=zeros(nFr*len2,1);
x_old = zeros(len1,1);
for j=1:nFr
%     plot(x_old); hold on;
%     plot(xi(1:len1),'r'); hold off;
%     legend('x_old', 'xi'); pause;
    xfinal(k:k+len1-1) = x_old + xi(1:len1,j);
%     xfinal(k:k+len1-1) = xi(1:len1,j);
    x_old = xi(1+len2:len,j);
    k = k + len2;

%     plot(xfinal(1:k)); pause;
end


text:
function GraphicGUI(state);

warning off MATLAB:divideByZero;

persistent Fs Frame Owner NdB

% All inputs to the blocks are either in dBs or rads/seconds. Hence any
%value sent as inputs to the parEq blocks in the model should be converted
%to these or created as rads/sec.

% When calling PEQ.M, all gain values can be sent as absolute dB and freq
%values sent as rads/seconds; This way the program is kept independent of
%the incoming audio's sampling frequency. Only for the display absolute
%frequencies are displayed; but graphically chosen values are kept as
%rads/sec

% This gets the Frame variable and defines it locally;

color1 = [255 0 0]/255;     % Red
color2 = [0 0 255]/255;     % Blue
color3 = [0 0 0]/255;       % Black
axisColor = [159 188 191]/255;
figColor = [211 226 226]/255;

if nargin==0 | strcmp(state,'reset')

    Owner=gcs;
    NdB = str2num(get_param([Owner '/band1'],'NdB'));
    %Get NdB value from the model
   % This sets the Frame persistent variable
    eval(get_param(Owner, 'PreLoadFcn'));
    % This sets the Fs sampling rate variable persistent
    [y,Fs,bits] = wavread(get_param([Owner,'/From Wave File'],'FileName'),[1 100]);
    load GraphicEQ_state

    bandw1 = eq_state.BandWidth1;
    centerfreq1 =eq_state.CenterFreq1;
    pk1= eq_state.peakgain1;

    bandw2= eq_state.BandWidth2;
    centerfreq2 = eq_state.CenterFreq2;
    pk2=eq_state.peakgain2;

    bandw3 =eq_state.BandWidth3;
    centerfreq3 = eq_state.CenterFreq3;
    pk3=eq_state.peakgain3;

    %Initial Creation of the filters
    %Filter1
     [b1 a1] = peq(0,pk1,centerfreq1,bandw1,NdB);
   % [b1 a1] = peq(0,14.5,centerfreq1,bandw1,NdB);

    [h1,w1] = freqz(b1,a1);
    %Filter2
    [b2 a2] = peq(0,pk2,centerfreq2,bandw2,NdB);
    [h2,w2] = freqz(b2,a2);
    %Filter3
    [b3 a3] = peq(0,pk3,centerfreq3,bandw3,NdB);
    [h3,w3] = freqz(b3,a3);

    %Convert back to dB
    h1 = 20*log10(abs(h1));
    h2 = 20*log10(abs(h2));
    h3 = 20*log10(abs(h3));

    %-------------------BAND ONE-----------------------------------------
    %Create the plot
    if 1
        line1 = semilogx(w1,h1,'tag','line1'); hold on;
    else
        line1 = plot(w1,h1,'tag','line1'); hold on;  %  ENABLE THIS to
    % compare with xfer function in model
    end;

    set(line1,'color',color1);
    %axis off;
    grid on;
    xlim(([0 1])*pi);
    ylim([-15 15]);

    %Customize the looks of the Figure Window
    set(gcf,'menubar','none');

    % and the axis
    xlabels=[.01  .02  .05  .1 .2 .5 ];   % setup 1 2 5 sequence for xlabels
    set(gca,'xtick',xlabels*pi);

    for k=1:length(xlabels)
        xstring{k} = sprintf('%3.0f',(xlabels(k)*Fs/2));
    end;
    set(gca,'xticklabel',xstring)

    set(gcf,'color',figColor);
    set(gca,'color',axisColor);
    ScreenSize = get(0,'screensize');
    if ( exist('state') & ~strcmp(state,'reset')) | (nargin==0)
        set(gcf,'position',[ScreenSize(3)/3 ScreenSize(4)/2 516 350]);
    end
    set(gca,'position',[0.1000    0.1500    0.8150    0.6150]);
    set(gcf,'resize','off');
    set(gca,'xcolor','k','ycolor','k');

    %set(get(gca,'xlabel'),'String','Frequency in Hz','color','k','position',[.5,-2.2],'fontangle','italic');
    set(get(gca,'ylabel'),'String','Gain in dB','color','k','fontangle','italic');

    %Create Help Text at the top of the Figure Window
    htext = text(0.5,.5,'Click and drag on the Round Markers or the Lines');
    set(htext,'color','k','fontangle','oblique','fontweight','bold','tag','text1','position',[.0055 23.4028],...
        'EdgeColor','black','BackgroundColor','y');

    %Create Text box for Parameter Value display
    freqText = text(0.0057,19.6528,'Center Frequency:','color',figColor,'tag','freqText');

    peakText = text(0.0856,19.7917,'Peak Value:','color',figColor,'tag','peakText');

    bandwText = text(0.6207,19.9306,'Band Width:','color',figColor,'tag','bandwText');

    %Create the Note Text
    noteText = text(0.0036,-20.2083,'The overall response (yellow dashed line) is centered at 0dB.',...
        'color','k','tag','noteText','fontweight','demi','fontangle','italic',...
        'EdgeColor','black','BackgroundColor',[.7 .9 .7]);

    %Create a Reset Button
    reset = uicontrol(gcf);
    set(reset,'position',[460 5 50 15],'string','Reset');
    set(reset,'callback','GraphicGUI(''reset'')');
    %% ------------------ BAND 1 -------------------------------------------
    %Set parameters in the Simulink Model
    set_param([Owner,'/PeakGain1'],'value',num2str(pk1) );
    set_param([Owner,'/CenterFreq1'],'value',num2str(centerfreq1));
    set_param([Owner,'/BandWidth1'],'value',num2str(bandw1) );

    %Add Markers that will be used for manipulating the parameters
    if pk1>0
       [h1max,indexh1] = max(h1);
   else
       [h1max,indexh1] = min(h1);
   end

    w1max = w1(indexh1);
    % the marker is plotted as a line object
    gain1 = line(w1max,h1max,'Marker','o','MarkerEdgeColor','k', ...
        'MarkerFaceColor',color1,'tag','peakgain1','LineStyle','none','markersize',6);
    %-------------------END CODE FOR BAND ONE--------------------------------------

    %-------------------BAND TWO--------------------------------------------------

    %Create the plot
    %line2 = plot(w2*fs/(2*pi),h2,'b','tag','line2');
    line2 = semilogx(w2,h2,'tag','line2');
    set(line2,'color',color2);

    %Set parameters in the Simulink Model
    set_param([Owner,'/PeakGain2'],'value',num2str(pk2) );
    set_param([Owner,'/CenterFreq2'],'value',num2str(centerfreq2) );
    set_param([Owner,'/BandWidth2'],'value',num2str(bandw2) );

    %Add Markers that will be used for manipulating the parameters
    if pk2>0
        [h2max,indexh2] = max(h2);    % this does not work, it can be a max or a min!
    else
        [h2max,indexh2] = min(h2);    % this does not work, it can be a max or a min!
    end;

    w2max = w2(indexh2);
    % the marker is plotted as a line object
    gain2 = line(w2max,h2max,'Marker','o','MarkerEdgeColor','k', ...
        'MarkerFaceColor',color2,'tag','peakgain2','LineStyle','none','markersize',6);
    %-------------------END CODE FOR BAND TWO--------------------------------------

    %--------------------------BAND THREE------------------------------------------
    %Create the plot
    %line3 = plot(w3*fs/(2*pi),h3,'m','tag','line3');
    line3 = semilogx(w3,h3,'tag','line3');
    set(line3,'color',color3);
    %Set parameters in the Simulink Model
    set_param([Owner,'/PeakGain3'],'value',num2str(pk3));
    set_param([Owner,'/CenterFreq3'],'value',num2str(centerfreq3) );
    set_param([Owner,'/BandWidth3'],'value',num2str(bandw3));

    %Add Markers that will be used for manipulating the parameters
    if pk3>0
       [h3max,indexh3] = max(h3);
    else
       [h3max,indexh3] = min(h3);
    end
    w3max = w3(indexh3);
    % the marker is plotted as a line object
    gain3 = line(w3max,h3max,'Marker','o','MarkerEdgeColor','k', ...
        'MarkerFaceColor',color3,'tag','peakgain3','LineStyle','none','markersize',6);
    %-------------------END CODE FOR BAND THREE--------------------------------------

    sumH = h1+h2+h3;
    %Normalize the sum of responses to value bet -15 & +15
    normSumH = sumNorm(sumH);
    %Plot sum of all responses
    line4 = plot(w2,normSumH,'y--','tag','line4','linewidth',2,'hittest','off');
    hold off;

    %Store sum of all response in GCF to display sum of responses
    allH = normSumH;
    setappdata(gcf,'allH', allH);

    %Also save the equalizer responses and current parameters in a
    %structure: eqData.H ,eqData.centerFreq, eqData.bandW
    eqData.H = [h1 h2 h3];
    eqData.centerFreq = [centerfreq1 centerfreq2 centerfreq3];
    eqData.bandW = [bandw1 bandw2 bandw3];
    setappdata(gcf,'eqData', eqData);

    %Set the button down function
    set(gcf,'WindowButtonDownFcn','GraphicGUI(''down'')');

    %Eliminate flicker
    set(gcf,'DoubleBuffer','on');

    %Set the mouse move while button down function
    set(gcf,'WindowButtonMotionFcn','','WindowButtonUpFcn','');

    % Execute the WindowButtonDownFcn
elseif strcmp(state,'down')

    %Get the current complete info on filter responses
    EqData = getappdata(gcf,'eqData');
    cfreq = EqData.centerFreq;

    %Identify the Band that was clicked
    htype = get(gco,'type');

    %If Line is clicked, then set the Point Down information in the Figure
    if strcmp(htype,'line')
        tag = get(gco,'tag');
        tagIndex = eval(tag(end));

        set(gcf,'WindowButtonMotionFcn','GraphicGUI(''move'')', ...
            'WindowButtonUpFcn','GraphicGUI(''up'')');

        cp = get(gca,'CurrentPoint');
        xDown = cp(1,1);
        yDown = cp(1,2);

        setappdata(gcf,'pointDown',[xDown cfreq(tagIndex)]);

        text1 = findobj(gcf,'tag','text1');
        line4 = findobj(gcf,'tag','line4');
        if strcmp(tag,'peakgain1')|strcmp(tag,'peakgain2')|strcmp(tag,'peakgain3')
            set(text1,'string','Drag to move the peak around');
        elseif strcmp(tag,'line1')|strcmp(tag,'line2')|strcmp(tag,'line3')
            set(text1,'string','Drag to change the bandwidth');
        elseif strcmp(tag,'line4')
            set(text1,'string','This is the actual response. Change this by moving individual bands');
        end
    end

% Execute the WindowButtonMotionFcn
elseif strcmp(state,'move')

    %Get the current complete info on filter responses
    EqData = getappdata(gcf,'eqData');
    cfreq = EqData.centerFreq;
    H = EqData.H;
    bandw = EqData.bandW;

    %Find handles of blocks in Simulink Model whose value is set here
    text1 = findobj(gcf,'tag','text1');

    line1 = findobj(gcf,'tag','line1');
    gain1 = findobj(gcf,'tag','peakgain1');

    line2 = findobj(gcf,'tag','line2');
    gain2 = findobj(gcf,'tag','peakgain2');

    line3 = findobj(gcf,'tag','line3');
    gain3 = findobj(gcf,'tag','peakgain3');

    line4 = findobj(gcf,'tag','line4');

    freqText = findobj(gcf,'tag','freqText');
    peakText = findobj(gcf,'tag','peakText');
    bandwText = findobj(gcf,'tag','bandwText');

    cp = get(gca,'CurrentPoint');

    %-----------------------------------
    %If the user drags the point out of the axis reset
    %the corresponding point to the axis limits
    x = cp(1,1); xlims = get(gca,'xlim');
    %if x<xlims(1), x = xlims(1);end;
    if x<(45*2*pi/Fs), x = (45*2*pi/Fs) ;end; % Center frequency cannot be lower than 45 Hz
    if x>xlims(2), x = xlims(2);end;

    y = cp(1,2); ylims = get(gca,'ylim');
    if y<ylims(1), y = ylims(1);end;
    if y>ylims(2), y = ylims(2);end;
    %------------------------------------

    tag = get(gco,'tag');
    tagIndex = eval(tag(end));
    switch tag(end)
        case '1'
            col = color1;
        case '2'
            col = color2;
        case '3'
            col = color3;
    end
    %Get the original ButtonDown x,y coordinates
    xyDown = getappdata(gcf,'pointDown');
    allH = getappdata(gcf,'allH');

    %IF THE GAIN POINT IS MOVED ====>>>
    if strcmp(tag,'peakgain1')|strcmp(tag,'peakgain2')|strcmp(tag,'peakgain3')
        set(text1,'string','Change the Center Frequency and Peak Gain of the frequency band');
        myH = H(:,tagIndex);

        %Change the Shape and color of Marker while it's moving
        set(eval(['gain' tag(end)]),'xdata',x,'ydata',y,'marker','diamond','markersize',12,'markerfacecolor','y');
        [newb, newa] = peq(0,y,x,bandw(tagIndex),NdB);
        [newh,neww] = freqz(newb,newa);
        newh = 20*log10(abs(newh));
        [watever maxW] = max(newh);

        set(eval(['line' tag(end)]),'ydata',newh,'xdata',neww,'linewidth',2);
        H(:,tagIndex) = newh;

        %Update the sum response plot to improve readability
        allH = sumNorm(sum(H,2));
        set(line4,'ydata',allH,'xdata',neww);

        %Set the new parameters for the modified band
        EqData.H = H;
        cfreq(tagIndex) = x;
        EqData.centerFreq = cfreq;
        EqData.bandW = bandw;
        setappdata(gcf,'eqData',EqData);
        setappdata(gcf,'allH',allH);

        %Display the parameters as the user moves the UIobjects
        numFreq = numFormat(x,Fs);
        numBandw = numFormat(bandw(tagIndex),Fs);
        peakVal = sprintf('%0.1f',y);
        set(freqText,'string',['Center Frequency: ',numFreq],'color',col,'fontangle','italic');
        set(peakText,'string',['Peak Value: ',peakVal,'dB'],'color',col,'fontangle','italic');
        set(bandwText,'string',['Bandwidth: ',numBandw],'color',col,'fontangle','italic');

        drawnow

    %IF THE LINE (Bandwidth) IS MOVED ====>>>
    elseif strcmp(tag,'line1')|strcmp(tag,'line2')|strcmp(tag,'line3')
        set(text1,'string','Dragging the colored lines changes the bandwidth of the frequency band ');
        %Get the response of the current band
        myH = H(:,tagIndex);
        %Change the Shape and color of line while it's moving
        set(eval(['line' tag(end)]),'color','c','linewidth',3);

        %Get the value of the peak gain
        peakG = get(eval(['gain' tag(end)]),'ydata');

        %Compute new Bandwidth
        %bwDiff = sign(xyDown(1)-xyDown(2)) * (x-xyDown(1));
        %Difference between current pressed x and centerFreq of that band

        bwDiff = (x-xyDown(2));
        bandww = 2*abs(bwDiff);
        if bandww<50*2*pi/Fs
            bandww = 50*2*pi/Fs; %Minimum Bandwidth allowed is 50 hertz
        elseif bandww>(20000*2*pi/Fs) %Maximum Bandwidth allowed is 20000 hertz
            bandww = 20000*2*pi/Fs;
        end
        bandw(tagIndex) = bandww;

        %Center Frequency cfreq doesnt change

        %Compute New Filter parameters
        [newb, newa] = peq(0,peakG,cfreq(tagIndex),bandww,NdB);
        [newh,neww] = freqz(newb,newa);
        newh = 20*log10(abs(newh));

        %Update Line for the chosen band
        set(eval(['line' tag(end)]),'ydata',newh,'xdata',neww);
        set(eval(['gain' tag(end)]),'Marker','o','MarkerEdgeColor','k',...
             'MarkerFaceColor',col,'LineStyle','none','markersize',6);

        %Update the sum response plot along with a bias to improve
        %readability

        H(:,tagIndex) = newh;
        EqData.H = H;
        %EqData.centerFreq = cfreq;
        EqData.bandW = bandw;
        setappdata(gcf,'eqData',EqData);

        allH = sumNorm(sum(H,2));
        setappdata(gcf,'allH',allH);
        set(line4,'ydata',allH,'xdata',neww);

        peakVal = sprintf('%.1f',peakG);
        numFreq = numFormat(cfreq(tagIndex),Fs);
        numBandw = numFormat(bandw(tagIndex),Fs);

        set(freqText,'string',['Center Frequency: ',numFreq],'color',col,'fontangle','italic');
        set(peakText,'string',['Peak Value: ',peakVal,' dB'],'color',col,'fontangle','italic');
        set(bandwText,'string',['Bandwidth: ',numBandw],'color',col,'fontangle','italic');

        drawnow
    end;

% Execute the WindowButtonUpFcn
elseif strcmp(state,'up')

    tag = get(gco,'Tag');
    tagIndex = eval(tag(end));
    switch tag(end)
        case '1', col = color1;
        case '2', col = color2;
        case '3', col = color3;
    end

    eqData = getappdata(gcf,'eqData');
    H = eqData.H;
    myH = H(:,tagIndex);
    hmax = max(myH);
    hmin = min(myH);

    if abs(hmax) < 0.00001
        hmax = hmin;
    end

    cf = eqData.centerFreq;
    cfreq = cf(tagIndex);
    bw = eqData.bandW;
    bandw = bw(tagIndex);

    text1 = findobj(gcf,'tag','text1');
    gain1 = findobj(gcf,'tag','peakgain1');
    line1 = findobj(gcf,'tag','line1');

    gain2 = findobj(gcf,'tag','peakgain2');
    line2 = findobj(gcf,'tag','line2');

    gain3 = findobj(gcf,'tag','peakgain3');
    line3 = findobj(gcf,'tag','line3');

    %Set the Change in parameters in the Simulink Model
    set_param([Owner,'/PeakGain',tag(end)],'value',num2str(hmax)) ;
    set_param([Owner,'/CenterFreq' tag(end)] ,'value',num2str(cfreq));
    set_param([Owner,'/BandWidth' tag(end)],'value',num2str(bandw));

    allH = getappdata(gcf,'allH');

%     %Compute the total decibel energy in the final response, to be used for
%     %normalizing the output power
%     sumDB = trapz(allH); %Perform trapezoidal integration
%     if ~isnan(sumDB)|~isinf(sumDB)
%         gainDB = 25 + sign((4000 - sumDB))*0.7*log10(abs(4000 - sumDB)); %Compute the gain by which the output is to be muliplied
%         set_param([gcs,'/outputGain'],'db',num2str(gainDB));
%     end

    set(text1,'string','Click and drag on the Colored Round Markers or the Colored Lines');
    set(eval(['gain' tag(end)]),'Marker','o','MarkerEdgeColor','k',...
            'MarkerFaceColor',col,'LineStyle','none','markersize',6);
    set(eval(['line' tag(end)]),'color',col,'linewidth',1);

    set(gcf,'WindowButtonMotionFcn','');
    set(gcf,'WindowButtonUpFcn','');

elseif strcmp(state,'save_state')
    eq_state.peakgain1 =  str2num(get_param([Owner,'/PeakGain1'],'value'));
    eq_state.CenterFreq1= str2num(get_param([Owner,'/CenterFreq1'],'value'));
    eq_state.BandWidth1=  str2num(get_param([Owner,'/BandWidth1'],'value'));

    eq_state.peakgain2   =  str2num(get_param([Owner,'/PeakGain2'],'value'));
    eq_state.CenterFreq2 =  str2num(get_param([Owner,'/CenterFreq2'],'value'));
    eq_state.BandWidth2  =  str2num(get_param([Owner,'/BandWidth2'],'value'));

    eq_state.peakgain3   =  str2num(get_param([Owner,'/PeakGain3'],'value'));
    eq_state.CenterFreq3 =  str2num(get_param([Owner,'/CenterFreq3'],'value'));
    eq_state.BandWidth3  =  str2num(get_param([Owner,'/BandWidth3'],'value'));

    eq_state.AF_Gain   =    str2num(get_param([Owner,'/AF_Gain'],'gain'));

   save GraphicEQ_State  eq_state
   %%save junk eq_state

end

function out = numFormat(x,Fs)
%This function is for formatting the display of the frequency values
x = x*Fs/(2*pi); %Convert from radians/sec to absolute frequency in Hz
if (x/1000)<1
    out = [num2str(round(x)),'Hz'];
else
    s = sprintf('%0.3g',x/1000);
    out = [s,'kHz'];
end

function normSumH = sumNorm(x)
% Look for the center of the freq response range, and force this to be 0 dB.
normSumH= x -(max(x)+min(x))/2;

author: Sable

repo: mcbench-benchmarks


text:
clear
close all
clc;

addpath('./local_functions_demo1')
addpath '../libs/exportFig/';
addpath './fun4MeanShift';

path_to_matconvnet = '../libs/matconvnet-1.0-beta23_modifiedDagnn';
run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));

dataDir = './dataset';
%% read matconvnet model
% set GPU
flagSaveFig = true; % {true false} whether to store the result

h = 28;
w = 28;
panelSZ = round(64/2)*2;
stepSize = (panelSZ-h)/2;
posList = [1, stepSize, 2*stepSize];
posList = repmat(posList, 3, 1);
posList = cat(3, posList, posList');
posList = reshape(posList, [], 2);
N = size(posList,1)+1;
%% prepare imdb
imdb.path_to_dataset = './toydata_v3';

fileList = dir( fullfile(imdb.path_to_dataset, 'batch*mat') );
imdb.imgList = {fileList.name};
imdb.set = ones(1, length(imdb.imgList));

fileList = dir( fullfile(imdb.path_to_dataset, 'test*mat') );
imdb.imgList(end+1:end+length({fileList.name})) = {fileList.name};
imdb.set(end+1:end+length({fileList.name})) = 2;

imdb.sets.name = {'train', 'val'};

% meanValue = 0;
% for i = 1:sum(imdb.set==1)
%     if mod(i,100) == 0
%         fprintf('\t%d/%d\n', i, sum(imdb.set==1));
%     end
%     tmpMat = load(fullfile(imdb.path_to_dataset, imdb.imgList{i}));
%     meanValue = meanValue + mean(tmpMat.imgMat(:));
% end
% meanValue = meanValue / sum(imdb.set==1);

meanValue = reshape([0.0717, 0.0820, 0.0715], [1 1 3]);
imdb.meta.meanvalue = meanValue; % 0.0981
imdb.meta.className = {};
for i = 1:9
    imdb.meta.className{end+1} = int2str(i);
end
imdb.meta.className{end+1} = '10';
imdb.meta.className{end+1} = 'void';
imdb.meta.classNum = length(imdb.meta.className);
imdb.meta.height = panelSZ;
imdb.meta.width = panelSZ;

save('imdb_toydata_v3_from_mnist.mat', 'imdb')

author: aimerykong