text:
%PROBLEM SPECIFICATION
% Specification of parameters, grid and boundary conditions

global Re

tic
if geval == 1	% Horizontal Poiseuille flow to the right
  dt = 0.1;	% Time step
  tend = 10;	% End time
  Re = 100;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 1/2;  % Omega-scheme for time-stepping
  umax = 0.75;	% Estimates of max. velocity components required for stability
  vmax = 0;	% condition for Adams-Bashforth-crank-Nicolson scheme

  xseglen = [1,1];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1,1];	% Similar to xseglen in y-direction.
  nx = [4,4];	% Number of cells along x-segments.
  ny = [4,4];	% Number of cells along y-segments.

  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [1,1; 1,1];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [2,2; 3,3];

elseif geval == 2	% Vertical Poiseuille flow upward
  dt = 0.5;	% Time step
  tend = 10;	% End time
  Re = 100;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 1/2;  % Omega-scheme for time-stepping
  umax = 0;	% Estimates of max. velocity components required for stability
  vmax = 0.75;	% condition for Adams-Bashforth-crank-Nicolson scheme

  xseglen = [1,1];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1,1];	% Similar to xseglen in y-direction.
  nx = [8,16];	% Number of cells along x-segments.
  ny = [8,8];	% Number of cells along y-segments.

  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [2,2; 3,3];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [1,1; 1,1];

elseif geval == 3	% Backward facing step
  dt = 3;	% Time step
  tend = 140;	% End time
  Re = 100;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 1/2;  % Omega-scheme for time-stepping
  umax = 1.25;	% Estimates of max. velocity components required for stability
  vmax = 0.25;	% condition for Adams-Bashforth-Crank-Nicolson scheme
  		% (not used in program ns1)

  xseglen = [5,5];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1,1];	% Similar to xseglen in y-direction.
  nx = [20,10];	% Number of cells along x-segments.
  ny = [25,25];	% Number of cells along y-segments.

  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [1,1; 1,1];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [1,2; 3,3];

elseif geval == 4	% Driven cavity
  dt = 0.5;	% Time step
  tend = 25;	% End time
  Re = 200;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 0.55;  % Omega-scheme for time-stepping
  umax = 1;	% Estimates of max. velocity components required for stability
  vmax = 0;	% condition for Adams-Bashforth-crank-Nicolson scheme

  xseglen = [1,1];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1,1];	% Similar to xseglen in y-direction.
  nx = [15,15];	% Number of cells along x-segments.
  ny = [15,15];	% Number of cells along y-segments.
  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [1,1; 2,2];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [1,1; 1,1];

elseif geval == 5	% Uniform flow under angle alpha in [0,pi/2] without
			% segmentation
  alpha = pi/4;
  dt = 3;	% Time step
  tend = 60;	% End time
  Re = 100;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 1/2;  % Omega-scheme for time-stepping
  umax = abs(cos(alpha));% Estimates of max. velocity components required for stability
  vmax = abs(sin(alpha));% condition for Adams-Bashforth-crank-Nicolson scheme

  xseglen = [1];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1];	% Similar to xseglen in y-direction.
  nx = [2];	% Number of cells along x-segments.
  ny = [2];	% Number of cells along y-segments.

  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [2; 3];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [2; 3];

elseif geval == 6	% Uniform flow under angle alpha in [0,pi/2] with
			% segmentation
  alpha = pi/4;
  dt = 3;	% Time step
  tend = 60;	% End time
  Re = 100;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 1/2;  % Omega-scheme for time-stepping
  umax = abs(cos(alpha));% Estimates of max. velocity components required for stability
  vmax = abs(sin(alpha));% condition for Adams-Bashforth-crank-Nicolson scheme

  xseglen = [1,1];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1,1];	% Similar to xseglen in y-direction.
  nx = [2,6];	% Number of cells along x-segments.
  ny = [2,3];	% Number of cells along y-segments.

  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [2,2; 3,3];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [2,2; 3,3];

elseif geval == 7	% Horizontal Poiseuille flow to the left
  dt = 0.5;	% Time step
  tend = 60;	% End time
  Re = 100;	% Reynolds number based on unit length and unit velocity
  central = 1;	% Enter 1 for central scheme or something else for upwind scheme
  omega = 1/2;  % Omega-scheme for time-stepping
  umax = 0.75;	% Estimates of max. velocity components required for stability
  vmax = 0;	% condition for Adams-Bashforth-crank-Nicolson scheme

  xseglen = [1,1];	% Enter segment lengths of horizontal boundary in
		% order of increasing x. Number of segments is arbitrary.
		% These segments are used for grid generation and boundary
		% conditions.
  yseglen = [1,1];	% Similar to xseglen in y-direction.
  nx = [4,4];	% Number of cells along x-segments.
  ny = [4,4];	% Number of cells along y-segments.

  % Type of boundary condition: 1: no-slip
  %				2: inflow
  %				3: outflow
  % Give type of boundary condition along lower (first row) and upper
  %(second row) horizontal segments in order of increasing x:

  xbc = [1,1; 1,1];

  % Give type of boundary condition along left (first row) and right
  % (second row) vertical segments in order of increasing y:

  ybc = [3,3; 2,2];

else
  error('Wrong value in input for parameter geval')
end

if (size(xseglen)~=size(nx))|(size(yseglen)~=size(ny))
  error('Wrong correspondence between xseglen,nx or yseglen,ny in problem_specification')
end
if (size(xbc(1,:))~=size(nx))|(size(ybc(1,:))~=size(nx))
  error('Wrong correspondence between nx,xbc or ny,ybc in problem_specification')
end

tijd = toc; disp(['problem_specification time = ',num2str(tijd)])


text:
function [V,F,I] = capsule(n,a)
  % [V,F,I] = capsule(n,a)
  %
  % Inputs:
  %   n  number of faces along latitude of each hemisphere
  %   a  length of cylindrical portion
  % Outputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices into rows of V
  %   I  #F list I(f) = 1,2 or 3 indicating bottom hemisphere, cylindrical part,
  %     top hemisphere, respectively.
  %
  % See also: stadium
  R = 1;
  s = 1;
  s = ceil( (2*n)/(2*pi*R)* a);
  [CV,CF] = cylinder_mesh(R,2*n,'Stacks',s);
  CV = CV.*[1 1 a];
  [X,Y,Z] = sphere(2*n);
  [SF,SV] = surf2patch(X,Y,Z,'triangles');
  [SV,~,J] = remove_duplicate_vertices(SV,1e-7);
  SF=J(SF);
  SF = SF(SF(:,1)~=SF(:,2) & SF(:,2)~=SF(:,3) & SF(:,3)~=SF(:,1),:);
  SBC = barycenter(SV,SF);
  T = SBC(:,3)>0;
  TV = SV+[0 0 a];
  TF = SF(T,:);
  BV = SV;
  BF = SF(~T,:);
  V = [BV;CV;TV];
  F = [BF;size(BV,1)+[CF;size(CV,1)+TF]];
  [V,~,J] = remove_duplicate_vertices(V,1e-7);
  F=J(F);
  [V,~,~,F] = remove_unreferenced(V,F);
  I = [repmat(1,size(BF,1),1);repmat(2,size(CF,1),1);repmat(3,size(TF,1),1)];
end

author: alecjacobson

repo: gptoolbox


text:
% Function to derive residual peak prominence parameter
%
% Description
% Function to derive residual peak prominence parameter
%
%
% Inputs
%  rep      : [samples] [Nx1]  Resonator output of LP-residual of speech
%                              signal
%  fs       : [Hz]      [1x1]  Sampling frequency
%  winLen   : [samples] [1x1]  Window length
%  x        : [samples] [Nx1]  Speech signal
%  Es       : [dB]      [Mx1]  Energy contour
%
% Outputs
%  peak_prom: [samples] [Px1]  Peak prominence contour
%  peak_t   : [samples] [Px1]  Time locations of Peak prominence contour
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [2] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [3] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
% Copyright (c) 2013 University of Mons, FNRS & 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>

function [peak_prom,peak_t,peak_idx] = get_res_peak_prom(rep,fs,winLen,x,Es)

%% Initial settings
Es=interp1(linspace(1,length(x),length(Es)),Es,1:length(x)); % Interpolated energy contour

ener_tol = 400/1000*fs;
ener_thresh=4;

peak_prom=zeros(1,round(length(rep)/winLen));
peak_t=zeros(1,round(length(rep)/winLen));
peak_idx=zeros(1,round(length(rep)/winLen));
winLenRem=winLen*.2; % just under half so periodic signals will have space for another peak

start=1;
stop=start+winLen-1;
n=1;
shift=20/1000*fs;

%% Do processing
while stop <= length(rep)

    frame_init = rep(start:stop);
    peak_t(n)=round(mean([start stop]));
    [~,maxIdx]=max(abs(frame_init));

    if maxIdx+start-round(winLen/2) < 1 || start+maxIdx+round(winLen/2) > length(rep)
        peak_prom(n)=0;
    else
        start2=maxIdx+start-round(winLen/2);
        stop2=start+maxIdx+round(winLen/2);
        frame2=rep(start2:stop2);
        [~,maxIdx]=max(abs(frame2));
        if frame2(maxIdx)< 0
            frame2=frame2*-1;
        end

        peak_idx(n) = maxIdx+start2-1;
        midpoint = round(winLen/2);
        frame2=frame2/max(frame2); % Normalise amplitude
        frame2(round(midpoint-winLenRem):round(midpoint+winLenRem))=0;
        frame2(frame2<0)=0;

        peak_prom(n)=1-max(frame2);

        if maxIdx+start-ener_tol < 1
            start_ener=1;
        else start_ener = maxIdx+start-ener_tol;
        end
        if maxIdx+start+ener_tol > length(x)
            stop_ener = length(x);
        else stop_ener = maxIdx+start+ener_tol;
        end

        max_ener=max(Es(start_ener:stop_ener));

        if Es(maxIdx+start) < max_ener-ener_thresh
            peak_prom(n)=0;
        end
    end

    start=start+shift;
    stop=start+winLen-1;
    n=n+1;
end

% Median filtering to remove outliers
peak_prom=medfilt1(peak_prom,5);

peak_prom(peak_t==0)=[];
peak_idx(peak_t==0)=[];
peak_t(peak_t==0)=[];

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function[b]=maxmax(x)
%MAXMAX  MAXMAX(X)=MAX(X(ISFINITE(X)))
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2012 J.M. Lilly --- type 'help jlab_license' for details
b=max(x(isfinite(x(:))));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [erg,arg] = getopt(str,default,varargin)
%JM-16-mar-00
% [erg,arg] = getoption('dummy',default,'a',a,...,'dummy',dummyval,...,'z',z)
% erg = dummyval,
% arg = getoption('dummy','a',a,...,'z',z)
%
% (c) 2000 - Jan Modersitzki

  j   = max(find(strcmp(varargin,str)));
  if isempty(j),
    erg = default;
    arg = varargin;
  else
    erg = varargin{j+1};
    if j == length(varargin),
      arg = {varargin{1:j-1}};
    else
      arg = {varargin{1:j-1},varargin{j+2:end}};
    end;
  end;

  if ~strcmp(str,'getoptout'),
    out = getopt('getoptout',0,arg{:});
    if out, showopt(str,erg); end;
  end;

function showopt(str,erg)

  fprintf('%-20s ',str);

  if isnumeric(erg),

    if size(erg) == [1,1],
      if round(erg)==erg,
	fprintf('%d\n',erg);
      else
	fprintf('%12.6e\n',erg);
      end;

    else
      fprintf('is %dx%d\n',size(erg));
    end;

  elseif isstr(erg),

    fprintf('%s\n',erg);

  end;

return;
%----------------------------------------------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [f] = spm_cost_SHC_fxa(x,v,a,P)
% equations of motion for a foraging problem
% FORMAT [f] = spm_cost_SHC_fxa(x,v,a,P)
%
% x   - hidden states
% v   - exogenous inputs
% a   - action
% P   - parameters for mountain car
%
% returns f = dx/dt (see spm_cost_SHC_fx)
% These equations of motion model dissipative flow x.x and x.v on a flat
% potential and increases in physiological states x.q as radial basis
% functions of secrete locations. The agent has to discover these
% locations % using an appropriate policy. This generative process would
% also substitute for Morris water-maze simulations or unbounded saccades.
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_cost_SHC_fxa.m 3757 2010-03-08 11:41:53Z guillaume $

% location and radius of attractors A (only A.q attractors deliver reward)
%--------------------------------------------------------------------------
global A; X   = A.x(:,A.q);

% physical flow
%--------------------------------------------------------------------------
f   = x;
f.x = x.v;
f.v = a - x.x*2 - x.v*4;

% physiological flow
%--------------------------------------------------------------------------
for i = 1:size(X,2)
    f.q(i) = (norm(x.x - X(:,i)) < A.d) - x.q(i)/2;
end

% flow
%--------------------------------------------------------------------------
dt  = 1/8;
f.x = f.x*dt;
f.v = f.v*dt;
f.q = f.q*dt;

author: spm


text:
function [ a, seed ] = r82_uniform_ab ( b, c, seed )

%*****************************************************************************80
%
%% R82_UNIFORM_AB returns a random R82 in a given range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real B, C, the minimum and maximum values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(2), the randomly chosen values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ t(1:2), seed ] = r8vec_uniform_01 ( 2, seed );

  a(1:2) = ( 1.0 - t(1:2) ) * b + t(1:2) * c;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [pptPath status] = prfMakePPTFromImages(sessDir, dtNames);
% Using images created from the meshscript_[date] script for each object
% pRF session, create a power point deck with images from both hemispheres,
% for easy comparison.
%
%   [pptPath status] = prfMakePPTFromImages([sessDir=pwd], [dtNames]);
%
% ras, 08/2008.
if notDefined('sessDir'),	sessDir = pwd;			end
if notDefined('dtNames'),
	dtNames = {'Averages' 'FaceAverages' 'HouseAverages' 'CheckerboardAverages'};

	% alternate: include odd/even subsets
% 	dtNames = {'Averages' 'FaceAverages' 'HouseAverages' 'CheckerboardAverages' ...
%			   'Odd' 'Even' 'FaceOdd' 'FaceEven' 'HouseOdd' 'HouseEven' ...
% 				'CheckerboardOdd' 'ChecerkboardEven'};
end

%% get params
cd(sessDir);
load mrSESSION
sessName = mrSESSION.sessionCode;
fields = {'Variance Explained' 'Eccentricity' 'Polar Angle' 'pRF Diameter'};
viewAngles = {'Medial' 'LateralOccipital' 'VentralZoom'};
hemiOrder = {[2 1] [1 2] [2 1]}; % L=1, R=2: order on the slide

%% open the PPT file
pptFile = sprintf('pRF Mesh Images %s.ppt', sessName);
pptPath = fullfile(prfDir, 'results', pptFile);
[ppt op] = pptOpen(pptPath);

%% main loop
for dt = dtNames
	% make a title slide for this analysis
	hFig = figure('Color', 'w');
	axes; axis off;
	pptPaste(op, hFig, 'meta', dt{1});
	close(hFig);

	for v = viewAngles
		for f = fields
			% load the images
			lh = imread( sprintf('Images/%s LH %s %s.png', v{1}, f{1}, dt{1}) );
			rh = imread( sprintf('Images/%s RH %s %s.png', v{1}, f{1}, dt{1}) );

			% display in a figure
			hFig = figure('Color', 'w', 'Units', 'norm', 'Position', [.1 .1 .8 .8]);
			hAx(1) = subplot('Position', [0 0 .5 1]);
			hAx(2) = subplot('Position', [.5 0 .5 1]);

			if isequal( hemiOrder{cellfind(viewAngles, v{1})}, [1 2] );
				% display left, then right
				axes(hAx(1));  imagesc(lh);  axis image;  axis off;  title('LH')
				axes(hAx(2));  imagesc(rh);  axis image;  axis off;  title('RH')
			else
				% display right, then left
				axes(hAx(1));  imagesc(rh);  axis image;  axis off;  title('RH')
				axes(hAx(2));  imagesc(lh);  axis image;  axis off;	 title('LH')
			end

			% add a colorbar
			switch f{1}
				case 'Variance Explained',
					cmap = mrvColorMaps('blueredyellow', 256);
% 					cmap = hot(256);
					clim = [0 .8];
					colorWheel = 0;
				case 'Eccentricity',
					cmap = mrvColorMaps('hsvshort', 256);
					clim = [0 30];
					colorWheel = 0;
				case 'Polar Angle',
					cmap = hsv(256);
					clim = [0 360];
					colorWheel = 1;
				case 'pRF Diameter',
					cmap = flipud( mrvColorMaps('hsvshort', 256) );
% 					cmap = jet(256);
					clim = [0 15];
					colorWheel = 0;
			end
			cbar = cbarCreate(cmap, f{1}, 'Clim', clim, ...
							  'ColorWheel', colorWheel, ...
							  'ColorWheelStart', 180);
			hLeg = mrvPanel('below', .15);
			axes('Parent', hLeg, 'Units', 'norm', 'Position', [.3 .3 .4 .4]);
			cbarDraw(cbar, gca);

			% copy into PPT
			ttl = sprintf('%s %s', dt{1}, f{1});
			pptPaste(op, hFig, 'meta', ttl);

			% close figure
			close(hFig);
		end
	end
end

pptClose(op, ppt, pptPath);

fprintf('\n\n\t*** [%s]: All done! ***\n\n', mfilename);

return


text:
function M1 = vox2ras_0to1(M0)
% M1 = vox2ras_0to1(M0)
%
% Converts a 0-based vox2ras matrix to 1-based, ie,
% Pxyz = M0*[c r s 1]' = M1*[(c+1) (r+1) (s+1) 1]'
%

%
% vox2ras_0to1.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

M1 = [];

if(nargin ~= 1)
  fprintf('M1 = vox2ras_0to1(M0)\n');
  return;
end

Q = zeros(4);
Q(1:3,4) = ones(3,1);

M1 = inv(inv(M0)+Q);

return;


text:
function pyra = projectpyramid(model, pyra)

% pyra = projectpyramid(model, pyra)
%
% Project feature pyramid pyra onto PCA eigenvectors stored
% in model.coeff.

for i = 1:length(pyra.feat)
  pyra.feat{i} = project(pyra.feat{i}, model.coeff);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% FT_PREAMBLE_LOADVAR is a helper script that optionally loads one or
% multiple FieldTrip data structures from mat files on disk, as an
% alternative to the user specifying the data structures as input variables
% to the calling function. This makes use of the cfg.inputfile variable.
%
% Furthermore, this function computes the MD5 hash of the input data structures for
% provenance.
%
% Use as
%   ft_preamble loadvar data
%   ft_preamble loadvar source mri
%   ft_preamble loadvar varargin
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_POSTAMBLE_SAVEVAR, FT_PREAMBLE_PROVENANCE

% Copyright (C) 2011-2019, Robert Oostenveld, DCCN
%
% 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$

% use an anonymous function
assign = @(var, val) assignin('caller', var, val);

if (isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)) || exist('Fief7bee_reproducescript', 'var')
  % the input data should be read from file

  if exist('Fief7bee_reproducescript', 'var')
    % the script, input and output files are written to a directory

    % write the function input variables to a MATLAB file
    if isequal(preamble_argin, {'varargin'})
      cfg.inputfile = {};
      iW1aenge_now = datestr(now, 30);
      for i=1:(ft_nargin-1)
        cfg.inputfile{i} = make_or_fetch_inputfile(Fief7bee_reproducescript,...
          sprintf('%s_%s_input_varargin_%d.mat', iW1aenge_now, FjmoT6aA_highest_ft, i),...
          'data', varargin{i});
      end
    else
      cfg.inputfile = {};
      iW1aenge_now = datestr(now, 30);
      for i=1:(ft_nargin-1)
        cfg.inputfile{i} = make_or_fetch_inputfile(Fief7bee_reproducescript,...
          sprintf('%s_%s_input_%s.mat', iW1aenge_now, FjmoT6aA_highest_ft, preamble_argin{i}),...
          preamble_argin{i}, eval(preamble_argin{i}));
      end
    end

  elseif ft_nargin>1
    ft_error('cfg.inputfile should not be used in conjunction with giving input data to this function');
  end

  if isfield(cfg, 'inputlock') && ~isempty(cfg.inputlock)
    mutexlock(cfg.inputlock);
  end

  if isequal(preamble_argin, {'varargin'}) && ~iscell(cfg.inputfile)
    % this should be a cell-array, otherwise it cannot be assigned to varargin
    cfg.inputfile = {cfg.inputfile};
  end

  if iscell(cfg.inputfile)
    if isequal(preamble_argin, {'varargin'})
      % read multiple inputs and copy them into varargin
      tmp = {};
      for i=1:length(cfg.inputfile)
        tmp{i} = loadvar(cfg.inputfile{i}, 'data');
      end % for
      assign('varargin', tmp);
      clear i tmp
    else
      % preamble_argin is a cell-array containing the variable names
      for i=1:length(cfg.inputfile)
        assign(preamble_argin{i}, loadvar(cfg.inputfile{i}, preamble_argin{i}));
      end % for
      clear i
    end
  else
    % preamble_argin{1} contains the name of the only variable
    assign(preamble_argin{1}, loadvar(cfg.inputfile, preamble_argin{1}));
  end

  if isfield(cfg, 'inputlock') && ~isempty(cfg.inputlock)
    mutexunlock(cfg.inputlock);
  end

end % if inputfile

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/mmv_phase_transition_snr_30db_s_8.mat';
N = 256;
S = 8;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 30;
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function [output] = cellall(x, val)

if nargin==1, val = 1; end
for k = 1:prod(size(x))
  tmp(k) = all(x{k}==val);
end
output = all(tmp==1);


text:
function [ nuniq, b1, b2 ] = i4vec2_sorted_unique ( n, a1, a2 )

%*****************************************************************************80
%
%% I4VEC2_SORTED_UNIQUE keeps unique elements of a sorted I4VEC2.
%
%  Discussion:
%
%    Item I is stored as the pair A1(I), A2(I).
%
%    The items must have been sorted, or at least it must be the
%    case that equal items are stored in adjacent vector locations.
%
%    If the items were not sorted, then this routine will only
%    replace a string of equal values by a single representative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of items.
%
%    Input, integer A1(N), A2(N), the array of N items.
%
%    Output, integer NUNIQ, the number of unique items.
%
%    Output, integer B1(N), B2(N), an array of NUNIQ unique items.
%
  nuniq = 0;

  if ( n <= 0 )
    return
  end

  itest = 1;
  nuniq = 1;
  b1(nuniq) = a1(itest);
  b2(nuniq) = a2(itest);

  for itest = 2 : n

    if ( a1(itest) ~= b1(nuniq) | a2(itest) ~= b2(nuniq) )

      nuniq = nuniq + 1;

      b1(nuniq) = a1(itest);
      b2(nuniq) = a2(itest);

    end

  end

  return
end


text:
%GETDERIVKERNELS  Returns filter coefficients for computing spatial image derivatives
%
%     [kx, ky] = cv.getDerivKernels('OptionName', optionValue, ...)
%
% ## Output
% * __kx__ Output matrix of row filter coefficients. It has the type `KType`.
% * __ky__ Output matrix of column filter coefficients. It has the type
%   `KType`.
%
% ## Options
% * __Dx__ Derivative order in respect of x. default 1
% * __Dy__ Derivative order in respect of y. default 1
% * __KSize__ Aperture size. It can be 'Scharr', 1, 3, 5, or 7. default 3.
% * __Normalize__ Flag indicating whether to normalize (scale down) the filter
%   coefficients or not. Theoretically, the coefficients should have the
%   `denominator = 2^(KSize*2-Dx-Dy-2)`. If you are going to filter
%   floating-point images, you are likely to use the normalized kernels. But
%   if you compute derivatives of an 8-bit image, store the results in a
%   16-bit image, and wish to preserve all the fractional bits, you may want
%   to set `Normalize=false`. default false
% * __KType__ Type of filter coefficients. It can be `single` or `double`.
%   default `single`
%
% The function computes and returns the filter coefficients for spatial image
% derivatives. When `KSize='Scharr'`, the Scharr 3x3 kernels are generated
% (see cv.Scharr). Otherwise, Sobel kernels are generated (see cv.Sobel). The
% filters are normally passed to cv.sepFilter2D.
%
% See also: cv.sepFilter2D, cv.getDerivKernels, cv.getStructuringElement
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function savestl(node,elem,fname,solidname)
%
% savestl(node,elem,fname,solidname)
%
% save a tetrahedral mesh to an STL (Standard Tessellation Language) file
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2012/02/14
%
% input:
%      node: input, surface node list, dimension Nx3
%      elem: input, tetrahedral element list; if size(elem,2)==3, it is a surface
%      fname: output file name
%      solidname: an optional string for the name of the object
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fid=fopen(fname,'wt');
if(fid==-1)
    error('You do not have permission to save mesh files.');
end

if(nargin<4)
    solidname='';
end
fprintf(fid,'solid %s\n', solidname);

if(nargin<3)
    error('incomplete input');
end
if(isempty(node) || size(node,2)<3)
    error('invalid node input');
end
if(~isempty(elem))
  if(size(elem,2)>=5)
	elem(:,5:end)=[];
  end
  if(size(elem,2)==4)
	elem=meshreorient(node,elem);
	elem=volface(elem);
  end
  ev=surfplane(node,elem);
  ev=ev(:,1:3)./repmat(sqrt(sum(ev(:,1:3).*ev(:,1:3),2)),1,3);
  len=size(elem,1);
  for i=1:len
      fprintf(fid,['facet normal %e %e %e\n' ...
                   ' outer loop\n' ...
                   '  vertex %e %e %e\n' ...
                   '  vertex %e %e %e\n' ...
                   '  vertex %e %e %e\n' ...
                   ' endloop\n' ...
                   'endfacet\n'], [ev(i,:)',node(elem(i,:),1:3)']);
  end
end

fprintf(fid,'endsolid %s\n', solidname);

fclose(fid);


text:
function [tf, bool,ind] = isvar(T,varName)
%checks if varName is a column of table T

if isstruct(T)
    tf = isfield(T,varName);
    bool = tf;
    ind = 1;
else
    bool = ismember(T.Properties.VariableNames,varName);
    tf = any(bool);
    ind = find(bool);
end

end

author: opencobra

repo: cobratoolbox


text:
function acv= autocov(x,maxlag)
%
% autocov computes the sample autocovariance of a time series x for lags
% from 0 to maxlag, returning a column vector of length maxlag+1.  x must
% be a column vector having length m not less than maxlag+1.  If no value
% is supplied for maxlag, the default is the minimum of m-1 and 100.
%
% Dr. Phillip M. Feldman
% Last update 21 June 2008
%
% Based on equations on p. 19 of "Introduction to Time Series and
% Forecasting" 2nd Edition by Brockwell and Davis.

% Section 1: Check input arguments and supply default value for maxlag
% if needed.

if nargin < 1, error('Missing input vector.'); end

[m n]= size(x);
if (n ~= 1)
   error('x must be a column vector.')
end
if (m <= maxlag)
   error('The length of the input vector x must be at least maxlag+1.');
end

if nargin < 2, maxlag= min(m-1,100); end

% Section 2: Compute autocovariance.

% Remove mean from x:
x= x - mean(x);

% For faster running time, we pre-allocate the output array:
acv= zeros(maxlag+1,1);

% Compute autocovariance:
for h= 0 : maxlag
   % Take matrix product of row vector and column vector to obtain a
   % scalar.  The row vector contains the first n-h elements of x; the
   % column vector contains the last n-h elements of x.
   acv(h+1)= x(1:m-h)' * x(1+h:m);
end

acv= acv / m;

author: Sable

repo: mcbench-benchmarks


text:
%Copyright (c) October,15 2008 by Varsha Hedau, UIUC.  All rights reserved.
function [Xpnts2,Vote2,VoteArr2] = RemoveRedundantPoints(Xpnts,Vote,VoteArr,w,h)

%%%%% remove redundant points %%%%%%
currid=1;
done=0;
Xpnts2=Xpnts;
Vote2=Vote;
VoteArr2=VoteArr;

while size(Xpnts,1)>0
    Xpnts2(currid,:) = Xpnts(1,:);
    Vote2(currid) = Vote(1);
    VoteArr2(:,currid) = VoteArr(:,1);
    currid=currid+1;
    dists = (Xpnts(1,1)-Xpnts(:,1)).^2 + ...
        (Xpnts(1,2)-Xpnts(:,2)).^2;
    if sqrt((Xpnts(1,1)-w)^2+(Xpnts(1,2)-h)^2)/sqrt((w/2)^2+(h/2)^2) < 1
        thres=10;
    else
        thres=20*(sqrt((Xpnts(1,1)-w)^2+(Xpnts(1,2)-h)^2)/sqrt((w/2)^2+(h/2)^2));
    end
    inds=find(dists>thres);
    Xpnts = Xpnts(inds,:);
    Vote = Vote(inds,:);
    VoteArr = VoteArr(:,inds);
end
Xpnts2 = Xpnts2(1:currid-1,:);
Vote2 = Vote2(1:currid-1);
VoteArr2 = VoteArr2(:,1:currid-1);

return;

author: zouchuhang

repo: LayoutNet


text:
function data = rawsignal2power(data,T)
% Gets *square of* the power time series using the Hilbert transform,
% as done in Baker et al (2014), eLife.
%
% Diego Vidaurre, OHBA, University of Oxford (2017)

if isstruct(data), ndim = size(data.X,2);
else, ndim = size(data,2);
end

for n = 1:length(T)
    ind = sum(T(1:n-1))+1:sum(T(1:n));
    if isstruct(data)
        for j = 1:ndim
            data.X(ind,j) = abs(hilbert(data.X(ind,j)));
        end
    else
        for j = 1:ndim
            data(ind,j) = abs(hilbert(data(ind,j)));
        end
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
clear all;
randn('seed',0);
% Data are generated
X=randn(20000,100);
X=X./repmat(sqrt(sum(X.^2)),[size(X,1) 1]);

% parameter of the optimization procedure are chosen
param.numThreads=-1; % number of processors/cores to use; the default choice is -1
                    % and uses all the cores of the machine

param.pos=0;
param.mode=1;       % projection on the l1 ball
param.thrs=2;
tic
X1=mexSparseProject(X,param);
t=toc;
toc
fprintf('%f signals of size %d projected per second\n',size(X,2)/t,size(X,1));
fprintf('Checking constraint: %f, %f\n',min(sum(abs(X1))),max(sum(abs(X1))));

param.mode=2;       % projection on the Elastic-Net
param.lambda1=0.15;

tic
X1=mexSparseProject(X,param);
t=toc;
toc
fprintf('%f signals of size %d projected per second\n',size(X,2)/t,size(X,1));
constraints=sum((X1.^2))+param.lambda1*sum(abs(X1));
fprintf('Checking constraint: %f, %f\n',min(constraints),max(constraints));

param.mode=6;       % projection on the FLSA
param.lambda1=0.7;
param.lambda2=0.7;
param.lambda3=1.0;

X=rand(2000,100);
X=X./repmat(sqrt(sum(X.^2)),[size(X,1) 1]);

tic
X1=mexSparseProject(X,param);
t=toc;
toc
fprintf('%f signals of size %d projected per second\n',size(X,2)/t,size(X,1));
constraints=0.5*param.lambda3*sum(X1.^2)+param.lambda1*sum(abs(X1))+param.lambda2*sum(abs(X1(2:end,:)-X1(1:end-1,:)));
fprintf('Checking constraint: %f, %f\n',mean(constraints),max(constraints));
fprintf('Projection is approximate (stops at a kink)\n',mean(constraints),max(constraints));

param.mode=6;       % projection on the FLSA
param.lambda1=0.7;
param.lambda2=0.7;
param.lambda3=1.0;

X=rand(2000,100);
X=X./repmat(sqrt(sum(X.^2)),[size(X,1) 1]);

tic
X1=mexSparseProject(X,param);
t=toc;
toc
fprintf('%f signals of size %d projected per second\n',size(X,2)/t,size(X,1));
constraints=0.5*param.lambda3*sum(X1.^2)+param.lambda1*sum(abs(X1))+param.lambda2*sum(abs(X1(2:end,:)-X1(1:end-1,:)));
fprintf('Checking constraint: %f, %f\n',mean(constraints),max(constraints));
fprintf('Projection is approximate (stops at a kink)\n',mean(constraints),max(constraints));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% Demo of using different line codings

bits = [1 0 1 0 0 0 1 1 0];
bitrate = 1; % bits per second

figure;
[t,s] = unrz(bits,bitrate);
plot(t,s,'LineWidth',3);
axis([0 t(end) -0.1 1.1])
grid on;
title(['Unipolar NRZ: [' num2str(bits) ']']);

figure;
[t,s] = urz(bits,bitrate);
plot(t,s,'LineWidth',3);
axis([0 t(end) -0.1 1.1])
grid on;
title(['Unipolar RZ: [' num2str(bits) ']']);

figure;
[t,s] = prz(bits,bitrate);
plot(t,s,'LineWidth',3);
axis([0 t(end) -1.1 1.1])
grid on;
title(['Polar RZ: [' num2str(bits) ']']);

figure;
[t,s] = manchester(bits,bitrate);
plot(t,s,'LineWidth',3);
axis([0 t(end) -1.1 1.1])
grid on;
title(['Manchester: [' num2str(bits) ']']);


text:
% SYNTAX:
% dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )
%
% UI NAME:
% OD_to_Conc
%
% DESCRIPTION:
% Convert OD to concentrations
%
% INPUTS:
% dod: the change in OD (#time points x #channels)
% SD:  the SD structure
% ppf: Partial path length factors for each wavelength. If there are 2 wavelengths
%      of data, then this is a vector of 2 elements.  Typical value is ~6 for each
%      wavelength if the absorption change is uniform over the volume of tissue measured.
%      To approximate the partial volume effect of a small localized absorption change
%      within an adult human head, this value could be as small as 0.1. Convention is
%      becoming to set ppf=1 and to not divide by the source-detector separation such that
%      the resultant "concentration" is in units of Molar mm (or Molar cm if those are the
%      spatial units). This is becoming wide spread in the literature but there is no
%      fixed citation. Use a value of 1 to choose this option.
%
% OUTPUTS:
% dc: the concentration data (#time points x 3 x #SD pairs
%     3 concentrations are returned (HbO, HbR, HbT)
%
% USAGE OPTIONS:
% Delta_OD_to_Conc: dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )
%
% PARAMETERS:
% ppf: [1.0, 1.0, 1.0]
%
function dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )

dc = [];
nWav = length(SD.Lambda);
ml = SD.MeasList;

if length(ppf) < nWav
    errordlg('WARNING: Data contains more than 3 wavelengths. Using PPF value of 1 for all wavelengths.');
    ppf = ones(1, nWav);
elseif length(ppf) > nWav
    d = length(ppf)-nWav;
    ppf(end-d+1:end) = [];
end

if ~isempty(find(ppf==1))
    ppf = ones(size(ppf));
end

nTpts = size(dod,1);

e = GetExtinctions( SD.Lambda );
if ~isfield(SD,'SpatialUnit')
    e = e(:,1:2) / 10; % convert from /cm to /mm
elseif strcmpi(SD.SpatialUnit,'mm')
    e = e(:,1:2) / 10; % convert from /cm to /mm
elseif strcmpi(SD.SpatialUnit,'cm')
    e = e(:,1:2) ;
end
einv = inv( e'*e )*e';

lst = find( ml(:,4)==1 );
for idx=1:length(lst)
    idx1 = lst(idx);
    idx2 = find( ml(:,4)>1 & ml(:,1)==ml(idx1,1) & ml(:,2)==ml(idx1,2) );
    rho = norm(SD.SrcPos(ml(idx1,1),:)-SD.DetPos(ml(idx1,2),:));
    if ppf(1)~=1
        dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)*rho*ppf))' )';
    else
        dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)))' )';
    end
end
dc(:,3,:) = dc(:,1,:) + dc(:,2,:);

author: BUNPC


text:
function parameters = connfreq_estimation(ns_ps,parameters)

% parameters = struct('n',2,'len',len_val,'ad',0.95,'as',0.8,'ap',0.2,'beta',0.8,'beta1',0.98,'gamma',0.998,'alpha',0.7,...
%             'pk',zeros(len_val,1),'noise_ps',ns_ps,'pxk_old',ns_ps,'pxk',ns_ps,'pnk_old',ns_ps,'pnk',ns_ps);

D = parameters.D;
len = parameters.len;
b = parameters.b;
V = parameters.V;
U = parameters.U;
gamma1 = parameters.gamma1;
gamma2 = parameters.gamma2;
alpha_max = parameters.alpha_max;
beta_min = parameters.beta_min;

alpha_c = parameters.alpha_c;
noise_ps = parameters.noise_ps;
Rmin_old = parameters.Rmin_old;
Pmin = parameters.Pmin;
Pmin_sw = parameters.Pmin_sw;
P = parameters.SmthdP;
u1 = parameters.u1;
j = parameters.j;
stored_min = parameters.stored_min;

P_noise_est = zeros(len,1);

%Spectral smoothing
%by equation(4)
% P_vector = zeros(2*D+1,1);
% P_y = ns_ps;
% for k = 1:len
%     if k>D&k<=(len - D)
%         for i = 1:2*D+1
%             P_vector(i) = ns_ps(k-D+i-1);%.^2;
%         end
%         P_y(k) = b*P_vector;
%     end
% end

P_y=smoothing(ns_ps,b,D);   % spectral smoothing according to Eq. 4

R = sum(P)/sum(ns_ps);
alpha_c_tild = 1/(1+(R-1)^2);
alpha_c = alpha_c*0.7 + 0.3*max(alpha_c_tild, 0.7);
alpha = (alpha_max*alpha_c)./((P./(noise_ps+eps) -1).^2 +1);

%
%Temporal smoothing
power_min = sum(Pmin);
power_noise_ps = sum(noise_ps); %Eq.(8)

%Pold=P; Pminold=Pmin; Pmin_swold=Pmin_sw;
%Decision=zeros(len,1);
% for k = 1:len
%     P(k) = alpha(k)*P(k) + (1-alpha(k))*P_y(k); %Eq.(5)
%     %Speech presence decision
%     if P(k)<Pmin(k)
%         Pmin(k) = P(k);
%     end
%     if P(k)<Pmin_sw(k)
%         Pmin_sw(k) = P(k);
%     end
%     if P(k)>gamma1*Pmin(k)
%         D_1(k) = 1;
%     else
%         D_1(k) = 0;
%     end
%
%     if P(k)>(Pmin(k) + gamma2*power_min/len)
%         D_2(k) = 1;
%     else
%         D_2(k) = 0;
%     end
%     Decision(k) = D_1(k)*D_2(k);
% end

Decision=zeros(len,1);
P=alpha.*P+(1-alpha).*P_y;
%Pmin2=Pminold;
Pmin=min(Pmin,P);
%Pmin_sw2=Pmin_swold;
Pmin_sw=min(Pmin_sw,P);
D_1a=zeros(len,1);
indx=find(P>gamma1*Pmin);
if ~isempty(indx),D_1a(indx)=1; end;
D_2a=zeros(len,1);
indx2=find(P>(Pmin+gamma2*power_min/len));
if ~isempty(indx2), D_2a(indx2)=1; end;
Decision=D_1a.*D_2a;

%Noise periodogram estimation
Rmin_tild = power_noise_ps/(power_min+eps); % Bias factor

if sum(Decision)>0
    Rmin = Rmin_old;
else
    Rmin = beta_min*Rmin_old + (1-beta_min)*Rmin_tild;  %Eq.(18)
end

% for k = 1:len
%     if Decision(k)==1
%         noise_ps(k) = Rmin*Pmin(k);
%     else
%         noise_ps(k) = ns_ps(k);%
%     end
% end

noise_ps=ns_ps;
indd=find(Decision==1);
if ~isempty(indd), noise_ps(indd)=Rmin*Pmin(indd);, end;

%Temporal minimum tracking
%use window to find the minimum
j = j+1;
if j==V
    stored_min(:,u1) = Pmin_sw;
    u1 = u1+1;
    if u1==U+1;
        u1=1;
    end
    Pmin = min(stored_min,[],2);
    Pmin_sw = P;
    j = 0;
end
%

parameters.alpha_c = alpha_c;
parameters.noise_ps = noise_ps;
parameters.Rmin_old = Rmin;
parameters.Pmin = Pmin;
parameters.Pmin_sw = Pmin_sw;
parameters.SmthdP = P;
parameters.u1 = u1;
parameters.j = j;
parameters.alpha = alpha;
parameters.stored_min = stored_min;
parameters.Decision = Decision;

% ----------------------------------------------
function y=smoothing (x,win,N);

len=length(x);
win1=win(1:N+1);
win2=win(N+2:2*N+1);
y1=filter(fliplr(win1),[1],x);

x2=zeros(len,1);
x2(1:len-N)=x(N+1:len);

y2=filter(fliplr(win2),[1],x2);

y=(y1+y2);

author: jtkim-kaist

repo: Speech-enhancement


text:
function F_enhance= RSWHEg_gray( F )
% Bi Bi histogram equalization with recursively seperated weighted histgram
% equalization with mean and mean
% clear all;
% close all;
% clc;
%close figure;
% x0=imread('pout.tif');
%x0=imread('eswar.jpg');
% x0=imread('54.tif');
% x0=imread('source24_1.tif');
% figure;
% imshow(x0);
% x0=imread('untitled.bmp');
 %x0=imread('plane.bmp');
%x0=uint8(x0);
%
%  %x0=imread('12.bmp');
% %u=double(i);
% h = [0.0086,0.0856,0.0086;
%      0.0856,0.8561,0.0856;
%      0.0086,0.0856,0.0086];
% h = h/sum1(h(:));  % Normalize the filter
% uSmooth = conv2padded(double(x0),double(h));
%figure(1);
%imshow(uint8(uSmooth));
 %luma=uint8(uSmooth(:,:,1));
%  XM=0;
% y0=rgb2ycbcr(x0);
luma=F;%(y0(:,:,1));
%luma=wiener2(luma);
%luma=luma-60;
[m,n]=size(luma);
%luma=luma;
% cb=y0(:,:,2);
% cr=y0(:,:,3);
k=mean(mean(luma));
XM=k;
XG=round(255/2);
[pixelCounts bins] = imhist(luma, 256);
[m12,xmax]=max(imhist(luma,256));
[m21,xmin]=min(imhist(luma,256));
Pmax=(max(imhist(luma,256))/(m*n));
Pmin=(min(imhist(luma,256))/(m*n));
beta=Pmax*((abs(XM-XG))/(xmax-xmin));
r=round(k);
l=length(luma(luma<=r));
         listindex=find(luma<=r);
 k1=reshape(luma(listindex),l,1);%obtained below fist mean values put it into an arry
 k2=mean(k1);%find mean of mean in first half
 k3=round(k2);%round to nearest integer to find pdf and cdf till first quarter
  sum1=0;
 l1=length(luma(luma<=k3));%find all  values in luma falls below fist mean
 listindex1=find(luma<=k3);%all  values in luma falls below fist mean given in their locations
 k4=reshape(luma(listindex1),l1,1);%arrange all fist mean values in a vector
     xpdf=hist(k4,[0:k3]);%pdf from 0:r
     xpdf=xpdf/(m*n);
     %alpa0=max(cumsum(xpdf/(m*n)));%from line 39 to 48 new rswhe code add this to other quarters
     alpa=(cumsum(xpdf));
     alpa0=alpa(end);
      %xpdf=xpdf/l1;%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form 0 to r.
      sum1=0;
      for l2=0:k3

          pw(l2+1)=(Pmax*((xpdf(l2+1)-Pmin)/(Pmax-Pmin))).^alpa0+beta;
          sum1=sum1+pw(l2+1);
      end
       for l2=0:k3
          pwn(l2+1)=pw(l2+1)/(sum1);
      end
      %plot(pwn);
      %(imhist(luma));
      %xlabel('gray levels up to mean');
      %ylabel('pdf up to mean');
      %title('histogram for half an image up to mean');
      %sk=xpdf*triu(ones(k3+1));
      sk=pwn*triu(ones(k3+1));
%       figure(2);
%       plot(sk);
%       xlabel('gray levels upto 1st mean');
%       ylabel('cdf upto mean');
%       title('cdf for half  of an image up to ist mean');
     % alpha=0.7;
      for l2=0:k3
          if(xpdf(l2+1)>0)
              list1=find(k4==l2);%find value in an vector i.e converted from matrix
              %list(list1)=alpha*sk(l2+1)*(k3+1)+(1-alpha)*(k3+1);
              list(list1)=(sk(l2+1)*(k3+1));%map dont disturb to get bhe as
              %list(list1)=alpha*sk(l2+1)*(k3+1)+(1-alpha)*(k3+1);
              %it is 13/3/2011
              ert(l2+1)=(sk(l2+1)*(k3+1));
          end
      end
          p=zeros(m,n);
      p(listindex1)= list; %listindex1 mxn locations list transformed values.
%      figure(3);
%      imshow(p);
%      xlabel('gray levels up to first mean');
%       ylabel('luma component equilized image up to first mean');
%      title('processed luma image up to first mean');
     k=mean(mean(luma));
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2nd quarter
r=round(k);
l=length(luma(luma<=r));
        listindex=find(luma<=r);
 k1=reshape(luma(listindex),l,1);%obtained below fist mean values put it into an arry
 k2=mean(k1);%find mean of mean in first half
 k3=round(k2);%round to nearest integer to find pdf and cdf till first quarter
%   sum1=0;
   b=k3;
%  l2=length(luma(luma>k3));%find all  values in luma falls below fist mean
  listindex2=find((luma>k3)&(luma<=r));%all  values in luma falls below fist mean given in their locations
  l2=length(listindex2);
  k5=reshape(luma(listindex2),l2,1);%arrange all fist mean values in a vector
      x2pdf=hist(k5,[k3+1:r]);%pdf from 0:r
       x2pdf=x2pdf/(m*n);%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form 0 to r.
       %alpa1=max(cumsum(x2pdf/(m*n)));%from line 39 to 48 new rswhe code add this to other quarters
       alpa=(cumsum(x2pdf));%/(m*n)))
       alpa1=alpa(end);
      %x2pdf=x2pdf/l2;%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form 0 to r.
      sum1=0;
      for l2k=0:r-(k3+1)
          pw1(l2k+1)=(Pmax*((x2pdf(l2k+1)-Pmin)/(Pmax-Pmin))).^alpa1+beta;
          sum1=sum1+pw1(l2k+1);
      end
       for l2k=0:r-(k3+1)
          pw1n(l2k+1)=pw1(l2k+1)/(sum1);
      end
%        figure(4);
%        plot(x2pdf);
%        xlabel('gray levels 2nd mean');
%        ylabel('pdf of 2nd mean');
%        title('histogram for 2nd mean');
       sk2=pw1n*triu(ones(r-k3));
%        figure(5);
%        plot(sk2);
%        xlabel('gray levels of 2nd mean');
%        ylabel('cdf upto mean');
%        title('cdf for half  of an image of 2nd mean');
       k2u=1;
       for l3=k3+1:r
           if(pw1n(k2u)>0)
               list2=find(k5==l3);%find value in an vector i.e converted from matrix
               %list1(list2)=alpha*(sk2(k2u)*(r-k3)+(k3+1))+(1-alpha)*(k3+1);
               list1(list2)=((k3+1)+(sk2(k2u))*(r-k3-1));%map dont disturb to
               %list1(list2)=alpha*(sk2(k2u)*(r-k3)+(k3+1))+(1-alpha)*(k3+1);
               %get BHE 13/3/2011
               ert(l3)=((k3+1)+(sk2(k2u))*(r-k3-1));
           end
               k2u=k2u+1;

       end

       p(listindex2)= list1;
%       figure(6);
%       imshow(p);
%      xlabel('gray levels up to first 2nd mean');
%       ylabel('luma component equilized image 2nd mean');
%      title('processed luma image up to 2nd mean');
 %lupper30=length(luma(luma>r);
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 lupper30=length(luma(luma>r));
%for i=0:r
    %if(luma(luma<=r))
        listindexupper3=find(luma>r);
   % end
%end
 k1upper=reshape(luma(listindexupper3),lupper30,1);
 mean3=mean(k1upper);
  r3=round(mean3);
 %length30=length((luma<=r3)&(luma>r));
  listindexupper30=find((luma>r)&(luma<=r3));
  length30=length(listindexupper30);
  k30upper=reshape(luma(listindexupper30),length30,1);

  %length30=length((luma<=r3)&(luma>r));
%  sum1=0;
alpa=0;
      xpdfupper30=hist(k30upper,[r+1:r3]);%pdf from r+1:r3
      xpdfupper30=xpdfupper30/(m*n);%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form r+1 to 255.
%alpa2=max(cumsum(xpdfupper30/(m*n)));%from line 39 to 48 new rswhe code add this to other quarters
alpa=(cumsum(xpdfupper30));%/(m*n)))
alpa2=alpa(end);
      %x3pdf=xpdfupper30/length30;%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form 0 to r.
sum1=0;
      for l2k=0:r3-(r+1)
          %pw2(l2k+1)=(Pmax*((x3pdf(l2k+1)-Pmin)/(Pmax-Pmin))).^alpa2+beta;
          pw2(l2k+1)=(Pmax*((xpdfupper30(l2k+1)-Pmin)/(Pmax-Pmin))).^alpa2+beta;
          sum1=sum1+pw2(l2k+1);
      end
       for l2k=0:r3-(r+1)
          pw2n(l2k+1)=pw2(l2k+1)/(sum1);
      end
% figure(7);
%      plot(pw2n);
%      xlabel('gray levels   3rd mean');
%      ylabel('pdf of 3rd mean');
%      title('histogram for upper half an image  3rd mean');
      skupper30=pw2n*triu(ones(r3-r));
%       figure(8);
%       plot(skupper30);
%       xlabel('gray levels after mean');
%       ylabel('cdf after mean');
%       title('cdf for upper half  of an image after mean');
      k3u=1;
      for k3upper=(r+1):r3
          if(pw2n(k3u)>0)
              list1upper30=find(k30upper==k3upper);%find value in an vector i.e converted from matrix
              listnew(list1upper30)=((r+1)+skupper30(k3u)*(r3-r-1));%map dont
              %listnew(list1upper30)=alpha*(skupper30(k3u)*(r3-r)+(r+1))+(1-alpha)*(r+1);
              %disturg to get original heq 14/3/2011
              %listnew(list1upper30)=alpha*(skupper30(k3u)*(r3-r)+(r+1))+(1-alpha)*(r+1);
              ert(k3upper)=((r+1)+skupper30(k3u)*(r3-r-1));
          end
         k3u=k3u+1;

     end

%           p2=zeros(m,n);

     p(listindexupper30)= listnew;
%     figure(9);
%      imshow(p);
%      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     listindexupper40=find(luma>r3);
     lupper40=length(listindexupper40);
   % end
%end
 k1upper40=reshape(luma(listindexupper40),lupper40,1);
 %sum1=0;
 %for i=0:r
 alpa=0;
     xpdfupper40=hist(k1upper40,[r3+1:255]);%pdf from r+1:255
     xpdfupper40=xpdfupper40/(m*n);%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form r+1 to 255.
alpa=(cumsum(xpdfupper40));%/(m*n)));%from line 39 to 48 new rswhe code add this to other quarters
alpa3=alpa(end);
     % x4pdf=xpdfupper40/lupper40;%normalized pdf to get nk/n,l=sum1 of xpdf,total no of pixels form 0 to r.
sum1=0;
      for l2k=0:(255-(r3+1))
          pw3(l2k+1)=(Pmax*((xpdfupper40(l2k+1)-Pmin)/(Pmax-Pmin))).^alpa3+beta;
          sum1=sum1+pw3(l2k+1);
      end
       for l2k=0:(255-(r3+1))
          pw3n(l2k+1)=pw3(l2k+1)/(sum1);
      end
% figure(10);
%      plot(xpdfupper40);
%      xlabel('gray levels after 4mean');
%      ylabel('pdf after 4mean');
%      title('histogram for upper half an image after 4mean');
     skupper40=pw3n*triu(ones(255-r3));
%      figure(11);
%      plot(skupper40);
%      xlabel('gray levels after 4mean');
%      ylabel('cdf after 4mean');
%      title('cdf for upper half  of an image after 4mean');
     k4u=1;
     for k4upper=(r3+1):255
         if(pw3n(k4u)>0)
             list1upper40=find(k1upper40==k4upper);%find value in an vector i.e converted from matrix
             %for k2u=1:58
             %listnew4(list1upper40)=alpha*(198+skupper40(k4u))*(255-r3)+(1-alpha)*(r3+1);
             listnew4(list1upper40)=(r3+1)+skupper40(k4u)*(255-r3-1);%map
             %listnew4(list1upper40)=alpha*((r3+1+skupper40(k4u))*(255-r3))+(1-alpha)*(r3+1);
                    % listnew4(list1upper40)=alpha*(198+skupper40(k4u)*(255-r3))+(1-alpha)*(r3+1);
             %dont disturb to get original bbhe 14/3/2011
         %end
         ert(k4upper)=(r3+1)+skupper40(k4u)*(255-r3-1);
         end
         k4u=k4u+1;
     end

%      for i=0:l-1
          %p3=zeros(m,n);
%          if (p(listindex))
%              p(:)=list;
%          end
%      end

     p(listindexupper40)= listnew4;
%     figure(12);
%     imshow(p);
%     xlabel('gray levels after 4mean');
%      ylabel('luma component equilized image after 4mean');
%      title('processed luma image after 4mean');
     ommmmm=p;
     F_enhance=ommmmm;
     end
     %ommmmm=p1+p;
%     figure(13);
%     imshow(ommmmm);
%      colormap('gray');
%       xlabel('gray level');
%      ylabel('combined lower and upper half luma component equilized image');
%      title('combined luma image');
% %     figure(8);
%     image(om);
% %     for j1=0:255
% %         count=0;
% %         for i1=0:m*n-1
% %             if om(i1+1)==j1
% %                 count=count+1;
% %             end
% %         end
% %             prob(j1+1)=count/m*n;
% %     end
% %     figure(16);
% %     plot(prob);
% %
% %    for j2=0:255
% %         count1=0;
% %         for i2=0:m*n-1
% %             if luma(i2+1)==j2
% %                 count1=count1+1;
% %             end
% %         end
% %             prob2(j2+1)=count1/m*n;
% %     end
% %     figure(17);
% %     plot(prob2);
% %     xlabel('gray levels after mean');
% %      ylabel('luma component equilized image after mean');
% %      title('processed luma image after mean');
% %     ommmmm=p1+p;
% %     figure(7);
% %      colormap('gray');
% %       xlabel('gray level');
% %      ylabel('combined lower and upper half luma component equilized image');
% %      title('combined luma image');
% %     image(ommmmm);
%k=imshow(ommmmm);
%      cat1=cat(3,ommmmm,cb,cr);
% %     figure(14);
% %     imshow(cat1);
% % %     xlabel('gray level(ycbcr)');
% %      ylabel('combined lower and upper half luma,cromablue,croma red component equilized image');
% %      title('luma croma b and r color processed image');
%     catconversion=ycbcr2rgb(uint8(cat1));
%     figure(15);
%     imshow(uint8(ommmmm));
%     xlabel('gray level(rgb)');
%      ylabel('combined lower and upper half  RGB component equilized image');
%      title('converted from ycbcr2rgb color(RGB) processed image');
% YM=0;
%      YM=mean(mean(ommmmm));
%      AMBE=abs(XM-YM);
%      disp('Absolute mean brightness error=');
%      disp(AMBE);
     %MSE = meanSquareError(luma, ommmmm);
    % figure(16);
%disp('mean Square Error = ');
%disp(MSE);
%PSNR1 = PSNR(luma, ommmmm);
%figure(17);
% disp('Peak Signal to Noise Ratio = ');
% disp(PSNR1);
% E=entropy(uint8(ommmmm));
% disp('Entropy=');
% disp(E);
%      figure(16);
%      (imhist(uint8(ommmmm)));

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function mv_plotScans(view,groupFlag);
% mv_plotScans(view,[groupFlag]);
%
% Shell/dialog for calling multi voxel UI for multiple scans.
%
% The idea here is to concatenate similar scans together and view the
% concatenated time course (mainly event-related or non-AB block designs,
% but can also be used if you have many cyclic scans and want to see all
% the cycles w/o averaging). If the scans have parfiles assigned, ras_tc will
% concatenate them together also.
%
% groupFlag, if passed as 1, rather than prompting for the scans to plot
% the code will look for a 'scanGroup' field in the dataTYPES.scanParams
% struct, which specifies which scans to use. If it can't find it, it
% will create it calling er_groupScans.
%
% For simplicity, assumes the data will be from the currently-selected ROI.
%
% 04/12/05 ras: adapted from tc_plotScans.
global dataTYPES HOMEDIR;

if nargin < 1
    help mv_plotScans;
    return
end

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

cdt = view.curDataType;

if groupFlag==1     % use pre-assigned group of scans
	[scans,dt] = er_getScanGroup(view);
else
    [scans ok] = er_selectScans(view,'View voxel data from which scans?');
    if ~ok  return;  end
    dt = cdt; % data type is current data type
end

view.curDataType = dt;

multiVoxelUI(view,[],scans);

view.curDataType = cdt; % switch back to orig. data type

return


text:
%% (Internal) Find unique elements in a vector within a win_size tolerance window
%
%     [unique_elements idx ] = unique_w_tolerance( val1, win_size )
%
%
% Arguments:
%
%   + val1: data elements,
%
%   + win_size: tolerance to consider val1(i) == val1(i+1).
%
% Output:
%
%   + sft_intersect1: the elements in the soft intersection
%
%   + idx1, idx2: the indexes of val1(idx1) which are in the soft
%   intersection.
%
% Example:
%
%
% See also soft_set_intersect
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate: 14/05/2016
% Last update: 14/05/2016
% Copyright 2008-2016
%
function [unique_elements, idx ] = unique_w_tolerance( val1, win_size )

if( nargin < 1 || isempty(val1) )
    unique_elements = [];
    idx = [];
    return;
end

if( nargin < 2 || isempty(win_size) )
    win_size = 1;
end

[val1, idx] = unique(val1);

if( length( val1) == 1 )
    unique_elements = val1;
    return;
end

bAux = true;
while( any(bAux) )
    bAux = abs(val1(1:end-1) - val1(2:end)) <= win_size;
    bAux = [bAux; bAux(end)];
    % keep first occurrence
    % avoid removing consecutive occurrences, keep only the first.
    % Use a logic function to do this.
    bAux1 = toeplitz(double([ false(3,1); bAux ]));
    bAux1 = logical(bAux1(3:end-1,1:4));
    bAux1 = bAux1(:,1) & (~bAux1(:,2)) | bAux1(:,1) & bAux1(:,3) & (~bAux1(:,4));

    val1(bAux1) = [];
    idx(bAux1) = [];

end

unique_elements = val1;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%demo_SOFAplotGeometry - Script demonstrating the usage of SOFAplotGeometry.

% #Author: Piotr Majdak
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or – as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% load a SOFA file in SimpleFreeFieldHRIR

SOFAfile=fullfile(SOFAdbPath,'database','widespread','ICO1m_00139.sofa');
disp(['Loading: ' SOFAfile]);
Obj=SOFAload(SOFAfile);

% plot all measurements
SOFAplotGeometry(Obj);
title(['Geometry SimpleFreeFieldHRIR, ' num2str(Obj.API.M) ' position(s)'])
set(gcf, 'Name', mfilename);

% % only show every 45th measurement
index = 1:45:Obj.API.M;
SOFAplotGeometry(Obj,index);
title(['Geometry SimpleFreeFieldHRIR, reduced to ' num2str(size(index,2)) ' position(s)'])
set(gcf, 'Name', mfilename);

% % %% load a SingleRoomDRIR SOFA file (outdated)
% disp(['Loading: ' 'db://' fullfile('database','thk','DRIR_LBS_VSA_1202RS_SBL.sofa')]);
% Obj=SOFAload(['db://' ...
%   fullfile('database','thk','DRIR_LBS_VSA_1202RS_SBL.sofa')]);
%
% % plot SOFA Object with 1202 Receivers
% SOFAplotGeometry(Obj);
% set(gcf, 'Name', mfilename)

% % remove all but one Receiver
Obj.ReceiverPosition = [0 0.09 0];
Obj.ReceiverPosition_Type = 'cartesian';
Obj.Data.IR = Obj.Data.IR(:,1,:);
Obj.Data.Delay = Obj.Data.Delay(:,1,:);
Obj = SOFAupdateDimensions(Obj);

SOFAplotGeometry(Obj);
title(['Geometry SimpleFreeFieldHRIR, ' num2str(Obj.API.R) ' receiver(s), ' num2str(Obj.API.M) ' position(s)'])
set(gcf, 'Name', mfilename);

% %% load a GeneralFIR SOFA file
SOFAfile=fullfile(SOFAdbPath,'database', 'tu-berlin','FABIAN_CTF_modeled.sofa');
Obj = SOFAload(SOFAfile);

SOFAplotGeometry(Obj);
title(['Geometry GeneralFIR, ' num2str(Obj.API.R) ' receiver(s), ' num2str(Obj.API.M) ' position(s)'])
set(gcf, 'Name', mfilename);

% % %% load example with room geometry (outdated)
% disp(['Loading: ' SOFAfile]);
% SOFAfile = fullfile(SOFAdbPath,'sofatoolbox_test', 'Oldenburg_OfficeII.sofa');
% Obj = SOFAload(SOFAfile);
%
% SOFAplotGeometry(Obj);
% set(gcf, 'Name', mfilename)

% %% if exists try plotting SOFA file containing spherical harmonic emitter
% if exist(fullfile(SOFAdbPath,'demo_SHforHRTFs_SH.sofa'))
%    Obj =  SOFAload(fullfile(SOFAdbPath,'demo_SHforHRTFs_SH.sofa'));
%     SOFAplotGeometry(Obj);
% else
%     error('Run demoSHforHRTFS.m first.')
% end

author: sofacoustics

repo: SOFAtoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_noiseless_s_32.mat';
N = 256;
S = 32;
pta = spx.pursuit.PhaseTransitionAnalysis(N);

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);


text:
function [X,Y,Z] = adjacency_plot_und(aij,coor)
%ADJACENCY_PLOT_UND     Quick visualization tool
%
%   [X,Y,Z] = ADJACENCY_PLOT(AIJ,COOR) takes adjacency matrix AIJ and node
%   spatial coordinates COOR and generates three vectors that can be used
%   for quickly plotting the edges in AIJ. If no coordinates are specified,
%   then each node is assigned a position on a circle. COOR can, in
%   general, be 2D or 3D.
%
%   Example:
%
%   >> load AIJ;                                % load your adjacency matrix
%   >> load COOR;                               % load 3D coordinates for each node
%   >> [x,y,z] = adjacency_plot_und(AIJ,COOR);  % call function
%   >> plot3(x,y,z);                            % plots network as a single line object
%
%   If COOR were 2D, the PLOT3 command changes to a PLOT command.
%
%   NOTE: This function is similar to MATLAB's GPLOT command.
%
%   Richard Betzel, Indiana University, 2013

n = length(aij);
if nargin < 2
    coor = zeros(n,2);
    for i = 1:n
        coor(i,:) = [cos(2*pi*(i - 1)./n), sin(2*pi*(i - 1)./n)];
    end
end

[i,j] = find(triu(aij,1));
[~, p] = sort(max(i,j));
i = i(p);
j = j(p);

X = [ coor(i,1) coor(j,1)]';
Y = [ coor(i,2) coor(j,2)]';
if size(coor,2) == 3
    Z = [ coor(i,3) coor(j,3)]';
end
if isfloat(coor) || nargout ~= 0
    X = [X; NaN(size(i))'];
    Y = [Y; NaN(size(i))'];
    if size(coor,2) == 3
        Z = [Z; NaN(size(i))'];
    end
end

X = X(:);
Y = Y(:);
if size(coor,2) == 3
    Z = Z(:);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function x = sample(p, n)
% SAMPLE    Sample from categorical distribution.
% Returns a row vector of integers, sampled according to the probability
% distribution p.
% Uses the stick-breaking algorithm.
% Much faster algorithms are also possible.

if nargin < 2
  n = 1;
end

cdf = cumsum(p(:));
for i = 1:n
  x(i) = sum(cdf < rand) + 1;
end


text:
function Population = EnvironmentalSelection(Population,W,N)
% The environmental selection of MOEA/IGD-NS

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,N);
    Next = FrontNo < MaxFNo;

    %% Select the solutions in the last front
    Last   = find(FrontNo==MaxFNo);
    Choose = LastSelection(Population(Last).objs,W,N-sum(Next));
    Next(Last(Choose)) = true;
    % Population for next generation
    Population = Population(Next);
end

function Remain = LastSelection(PopObj,W,K)
% Select part of the solutions in the last front

    N  = size(PopObj,1);
    NW = size(W,1);

    %% Calculate the distance between each solution and point
    Distance = pdist2(PopObj,W);
    Con      = min(Distance,[],2);

    %% Delete the solution which has the smallest metric contribution one by one
    [dis,rank] = sort(Distance,1);
    Remain     = true(1,N);
    while sum(Remain) > K
        % Calculate the fitness of outliers
        Outliers = Remain;
        Outliers(rank(1,:)) = false;
        METRIC   = sum(dis(1,:)) + sum(Con(Outliers));
        Metrics  = inf(1,N);
        Metrics(Outliers) = METRIC - Con(Outliers);
        % Calculate the fitness of other solutions
        for p = find(Remain & ~Outliers)
            temp = rank(1,:) == p;
            outliers = false(1,N);
            outliers(rank(2,temp)) = true;
            outliers = outliers & Outliers;
            Metrics(p) = METRIC - sum(dis(1,temp)) + sum(dis(2,temp)) - sum(Con(outliers));
        end
        % Delete the worst solution and update the variables
        [~,del] = min(Metrics);
        temp = rank ~= del;
        dis  = reshape(dis(temp),sum(Remain)-1,NW);
        rank = reshape(rank(temp),sum(Remain)-1,NW);
        Remain(del) = false;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [s, S_l] = idpLin2seg(l)

% IDPLIN2SEG IDP line to segment conversion
%   IDPLIN2SEG  returns a 3d segment with the two suppor points of the IDP
%   line L.
%
%   [s, S_l] = IDPLIN2SEG(L) returns the Jacobian wrt L.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1

    [e1, e2] = idpLin2idpPnts(l);
    p1       = idp2euc(e1);
    p2       = idp2euc(e2);

    s        = [p1;p2];

else

    [e1, e2, E1_l, E2_l] = idpLin2idpPnts(l);
    [p1, P1_e1]          = idp2euc(e1);
    [p2, P2_e2]          = idp2euc(e2);

    s   = [p1;p2];
    S_l = [P1_e1*E1_l ; P2_e2*E2_l];

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function dt =  testing(algo,dt)
%
%   Generates labels using knn.
%
  [numOfSam,vecDim,outDim]=get_dim(dt);

  % kernel or distance - that is the question ?!
  if isa( algo.child, 'distance')
    distflag = 1;
  else
    distflag = 0;
    trainXP=get_norm(algo.child,algo.dat).^2;
    testXP=get_norm(algo.child,dt).^2;
  end

%  flag=0;
%  if isa(algo.child,'distance')
%    if strcmp(algo.child.dist,'custom')
%      flag=1;
%    end
%    if strcmp(algo.child.dist,'custom_fast')
%      flag=1;
%    end
%  end
%  if flag==0 %% ok to compute norm
%   trainXP=get_norm(algo.child,algo.dat).^2;
%   testXP=get_norm(algo.child,dt).^2;
%  end

  y=get_y(algo.dat);
  yTemp=get_y(dt);
  yEst=yTemp*0-1;

% Possible to calc in batch or as for loop, whereas the batch way requires
% the whole kernel matrix to fit in memory.
if algo.batch==1,

  if distflag==0
    dist = trainXP*ones(1,length(testXP)) - 2*calc(algo.child,dt,algo.dat) + ones(length(trainXP),1)*testXP';
  else
    dist = calc(algo.child,dt,algo.dat);
  end

  [val,valInd]=sort(dist);

  if(~isempty(y))
      for i=1:numOfSam,
          yTmp(i,:) = mean(y(valInd(1:algo.k,i),:),1);
      end;
  else
   yTmp=[];
  end
 if algo.algorithm.use_signed_output==1,
   if size(y,2)>1  %% <---multi-class pattern recognition
     [maxComp maxCompInd]=max(yTmp,[],2);
     for i = 1:numOfSam,
       yEst(i,maxCompInd(i))=1;
     end;
   else    %%<--- binary pattern recognition
     yEst=sign(yTmp);
   end
 else
   yEst=yTmp;
 end

 if algo.output_preimage
   yEst=[]; ind=algo.dat.index;
   for i=1:algo.k
     yEst=[yEst ind(valInd(i,:))'];
   end
 end

else
%<<-----online calculations, instead of batch (does not need as much memory
%           as batch)----->

  for i=1:numOfSam
    dist= (trainXP - 2*get_kernel(algo.child,dt,algo.dat,i) ) + testXP(i);

    [val,valInd]=sort(dist);

   if algo.algorithm.use_signed_output==1,
    if size(y,2)>1  %% <---multi-class pattern recognition
      tmp=sum(y(valInd(1:algo.k),:),1);
      [maxComp maxCompInd]=max(tmp);
      yEst(i,maxCompInd)=1;
    else            %% binary pattern recognition
      yEst(i)=sign(mean(y(valInd(1:algo.k),:)));
    end
   else
     yEst(i,:)=mean(y(valInd(1:algo.k),:),1);
   end
  end

  if algo.output_preimage
    ind=algo.dat.index;
    for j=1:algo.k
      yEst(i,j)=[yEst ind(valInd(j,:))'];
    end
   end;

%<<-----------end online------------------>>
end;

  if algo.algorithm.use_signed_output==1 %% <--- break deadlocks
    yEst(yEst==0)=-1;
  end

  dt=set_x(dt,yEst);
  dt=set_name(dt,[get_name(dt) ' -> ' get_name(algo)]);


text:
D = spm_eeg_load('dspm8_faces_run1.mat');

montage.labelorg = D.chanlabels;

montage.labelnew = [montage.labelorg(1:128), 'HEOG', 'VEOG'];

tra = eye(D.nchannels);
tra(129:end, :) = [];
tra = detrend(tra, 'constant');

% HEOG
tra(129, [131 132]) = [1 -1];

% VEOG
tra(130, [135 136]) = [1 -1];

montage.tra = tra;

save faces_eeg_montage.mat montage

author: spm


text:
function h= rbmVtoH(model, v)
%go from visible to hidden based on type

if isequal(model.type, 'BB')
    h= logistic(v*model.W + repmat(model.b,size(v,1),1));
elseif isequal(model.type, 'BG')
    h= v*model.W + repmat(model.b,size(v,1),1);
elseif isequal(model.type, 'GB')
    h= logistic(v*model.W/model.sigma + repmat(model.b,size(v,1),1));
end


text:
%% test2.m - exp03
% This script shows histograms for errors computed in test1.m
%
% 02-08-10 Michal Uricar
% 07-03-11 Michal Uricar

% clc;
close all; clearvars;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Add path

addpath('./Functions/');

%% Load errors

load('./results/errors.mat');

%% Options

% show images
% opt.verbose = true;
opt.verbose = false;
% save images
opt.save = true;
% opt.save = false;

%% Show histograms

if (~exist('./img/', 'dir'))
    mkdir('./img/');
end;

% get screen size
scrsz = get(0,'ScreenSize');

% histogram of mean errors
x = 0:options.bw(1); y = round(L);
hist_mean = hist(y, x);
% cumulative histogram of mean errors
n_elements = histc(y, x); c_elements = cumsum(n_elements);
if (opt.verbose)
    figure; hist(y, x);
    set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
    set(gca, 'XTick', x); set(gca, 'xlim', [-0.5 options.bw(1)+0.5]);
    xlabel('Distance in pixels'); ylabel('Count of occurrences');
    set(gca, 'YTick', 0:50:2*c_elements(end));
    title('Histogram of mean errors');
    if (opt.save)
        saveas(gcf, './img/histogram_mean.png');
        saveas(gcf, './img/histogram_mean.fig');
    end;

    figure; bar(x, c_elements / max(c_elements));
    xlabel('RMS [px]'); ylabel('Count of occurrences [%]');
    set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
    set(gca, 'XTick', x); set(gca, 'xlim', [-0.5 options.bw(1)+0.5]);
    set(gca, 'YTick', 0:0.05:1);
    title('Cumulative historgam of mean errors');
    if (opt.save)
        saveas(gcf, './img/cumhist_mean.png');
        saveas(gcf, './img/cumhist_mean.fig');
    end;
end;

% histogram of max errors
x = 0:options.bw(1); y = round(err_maxdist);
hist_max = hist(y, x);
% cumulative histogram of max errors
n_elements = histc(y, x); c_elements = cumsum(n_elements);
if (opt.verbose)
    figure; hist(y, x);
    set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
    set(gca, 'XTick', x); set(gca, 'xlim', [-0.5 options.bw(1)+0.5]);
    xlabel('Distance in pixels'); ylabel('Count of occurrences');
    set(gca, 'YTick', 0:50:2*c_elements(end));
    title('Histogram of max errors');
    if (opt.save)
        saveas(gcf, './img/histogram_max.png');
        saveas(gcf, './img/histogram_max.fig');
    end;

    figure; bar(x, c_elements / max(c_elements));
    xlabel('RMS [px]'); ylabel('Count of occurrences [%]');
    set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
    set(gca, 'XTick', x); set(gca, 'xlim', [-0.5 options.bw(1)+0.5]);
    set(gca, 'YTick', 0:0.05:1);
    title('Cumulative historgam of max errors');
    if (opt.save)
        saveas(gcf, './img/cumhist_max.png');
        saveas(gcf, './img/cumhist_max.fig');
    end;
end;
% save histograms
save('./results/histograms.mat', 'hist_mean', 'hist_max');

%% Show mean and max error for each point from ground truth

mean_nose = mean(err_nose);
mean_canthus_rl = mean(err_canthus_rl);
mean_canthus_lr = mean(err_canthus_lr);
mean_mouth_corner_r = mean(err_mouth_corner_r);
mean_mouth_corner_l = mean(err_mouth_corner_l);
mean_canthus_rr = mean(err_canthus_rr);
mean_canthus_ll = mean(err_canthus_ll);

max_nose = max(err_nose);
max_canthus_rl = max(err_canthus_rl);
max_canthus_lr = max(err_canthus_lr);
max_mouth_corner_r = max(err_mouth_corner_r);
max_mouth_corner_l = max(err_mouth_corner_l);
max_canthus_rr = max(err_canthus_rr);
max_canthus_ll = max(err_canthus_ll);

fprintf('\n\n__________________________________________________________\n');
fprintf('\nResults:\n');
fprintf('Mean error on nose: \t\t\t\t%.4f\t\n', mean_nose);
fprintf('Mean error on mean_canthus_rl: \t\t%.4f\t\n', mean_canthus_rl);
fprintf('Mean error on mean_canthus_lr: \t\t%.4f\t\n', mean_canthus_lr);
fprintf('Mean error on mean_mouth_corner_r: \t%.4f\t\n', mean_mouth_corner_r);
fprintf('Mean error on mean_mouth_corner_l: \t%.4f\t\n', mean_mouth_corner_l);
fprintf('Mean error on mean_canthus_rr: \t\t%.4f\t\n', mean_canthus_rr);
fprintf('Mean error on mean_canthus_ll: \t\t%.4f\t\n', mean_canthus_ll);
fprintf('__________________________________________________________\n');

fprintf('\nR_mean_tst: \t\t\t\t\t%.5f\n', Rtst);
fprintf('R_max_tst: \t\t\t\t\t\t%.5f\n', sum(err_maxdist)/length(err_maxdist));

fprintf('__________________________________________________________\n');

%%

xSC = 0:ceil(max(L));
ySC = round(L);
hist_meanSC = hist(ySC, xSC);
% cumulative histogram of mean errors
n_elementsSC = histc(ySC, xSC);
c_elementsSC = cumsum(n_elementsSC);
% max
xSCmax = 0:ceil(max(L));
ySCmax = round(err_maxdist);
hist_maxSC = hist(ySCmax, xSCmax);
% cumulative histogram of max errors
n_elementsSCmax = histc(ySCmax, xSCmax);
c_elementsSCmax = cumsum(n_elementsSCmax);

fig = figure;
y1 = c_elementsSC / max(c_elementsSC) * 100;
plot(xSC(1:end), y1(1:end), 'b', 'LineWidth', 3); hold on; grid on;
set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
set(gca, 'XTick', xSC); set(gca, 'YTick', 0:5:100, 'FontSize', 15); set(gca, 'ylim', [0 105]);
l1 = legend('corners dataset, normalized trn loss', 'Location', 'SouthEast');
xlabel(['relative error [%]'], 'FontSize', 20);
ylabel('Count of occurrences [%]', 'FontSize', 20);
title('Cumulative historgam of mean errors', 'FontSize', 20);

% r = 781; c = 731;
r = 801; c = 801;
% Sets position and size of figure on the screen
set(fig, 'Units', 'pixels', 'position', [0 0 c r] );
% Sets axes to fill the figure space
% set(gca, 'Units', 'pixels', 'position', [0 0 c+1 r+1 ]);
% Sets print properties; Looks like 1 pixel = (3/4)th of a point
set(fig, 'paperunits', 'points', 'papersize', [fix((c-1)*(3/4))+1 fix((r-1)*(3/4))+1]);
set(fig, 'paperunits', 'normalized', 'paperposition', [0 0 1 1]);
print( fig, sprintf('-r%d', ceil(72*(4/3))), '-dpng', './img/cumhist-mean.png');

fig = figure;
y2 = c_elementsSCmax / max(c_elementsSCmax) * 100;
plot(xSCmax(1:end), y2(1:end), 'b', 'LineWidth', 3); hold on; grid on;
set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
set(gca, 'XTick', xSCmax); set(gca, 'YTick', 0:5:100, 'FontSize', 15);
set(gca, 'ylim', [0 105]);
l1 = legend('corners dataset, normalized trn loss', 'Location', 'SouthEast');
xlabel(['relative error [%]'], 'FontSize', 20);
ylabel('Count of occurrences [%]', 'FontSize', 20);
title('Cumulative historgam of max errors', 'FontSize', 20);

r = 801; c = 801;
% Sets position and size of figure on the screen
set(fig, 'Units', 'pixels', 'position', [0 0 c r] );
% Sets axes to fill the figure space
% set(gca, 'Units', 'pixels', 'position', [0 0 c+1 r+1 ]);
% Sets print properties; Looks like 1 pixel = (3/4)th of a point
set(fig, 'paperunits', 'points', 'papersize', [fix((c-1)*(3/4))+1 fix((r-1)*(3/4))+1]);
set(fig, 'paperunits', 'normalized', 'paperposition', [0 0 1 1]);
print( fig, sprintf('-r%d', ceil(72*(4/3))), '-dpng', './img/cumhist-max.png');

%% Timestamp

fprintf(1,'Finished on %s\n\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function [map,data] = PhaseMango(spikes,phases,varargin)

%PhaseMango - Compute phase as a function of spike rate and acceleration.
%
% Compute spike phase as a function of spike rate and acceleration, as in
% Harris et al. (2002). The name of the function refers to the aspect of
% the resulting phase plot.
%
%  USAGE
%
%    [map,data] = PhaseMango(spikes,phases,<options>)
%
%    spikes         spike timestamps
%    phases         phase samples (see Phase)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'smooth'      smoothing size in bins (0 = no smoothing, default = 2)
%     'nBins'       number of horizontal and vertical bins (default = [50 50])
%     'limits'      [maxRate minAccel maxAccel] (default = [8 -2 2])
%     'maxGap'      time gaps between successive position samples exceeding
%                   this threshold (e.g. undetects) will not be interpolated
%                   (default = 100 ms)
%     'boundaries'  onset and offset for single-lap phase precession can be
%                   determined either automatically based on spike count
%                   ('count', default) or using explicit firing field
%                   boundaries ([Xstart Xstop], where each X is in [0..1])
%    =========================================================================
%
%  OUTPUT
%
%    map.x                x bins
%    map.y                y bins
%    map.phase            phase map
%    map.count            count map
%
%    data.rate            spike rate for each cycle
%    data.acceleration    spike acceleration for each cycle
%    data.phase           mean spike phase (in radians) for each cycle
%
%  SEE
%
%    See also Phase, PhasePrecession, PlotPhaseMango.

% Copyright (C) 2004-2012 by Michaël Zugaro
%
% 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.

% Default values
maxGap = 0.1;
smooth = 2;
nBins = 50;
limits = [8 -2 2];

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help PhaseMango'' for details).');
end

% Check parameter sizes
if ~isdvector(spikes),
	error('Parameter ''spikes'' is not a vector (type ''help PhaseMango'' for details).');
end
if size(phases,2) ~= 2,
	error('Parameter ''phases'' is not a Nx2 matrix (type ''help PhaseMango'' for details).');
end
isradians(phases(:,2));

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help PhaseMango'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'smooth',
			smooth = varargin{i+1};
			if ~isdvector(smooth,'>=0') || length(smooth) > 2,
				error('Incorrect value for property ''smooth'' (type ''help PhaseMango'' for details).');
			end
		case 'nbins',
			nBins = varargin{i+1};
			if ~isivector(nBins,'>0') || length(nBins) > 2,
				error('Incorrect value for property ''nBins'' (type ''help PhaseMango'' for details).');
			end
		case 'limits',
			limits = varargin{i+1};
			if ~isivector(limits,'#3') || limits(1) <0 || limits(3) < limits(2),
				error('Incorrect value for property ''limits'' (type ''help PhaseMango'' for details).');
			end
		case 'maxgap',
			maxGap = varargin{i+1};
			if ~isdscalar(maxGap,'>0'),
				error('Incorrect value for property ''maxGap'' (type ''help PhaseMango'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help PhaseMango'' for details).']);
	end
end

% Default values
data.time = [];
data.rate = [];
data.acceleration = [];
data.phase = [];
map.x = [];
map.y = [];
map.phase = [];
map.count = [];

% Compute spike phases
if isempty(spikes), return; end
spikePhases = Interpolate(phases,spikes,'trim','off','type','circular');
spikePhases = exp(j*spikePhases(:,2));
if isempty(spikePhases), return; end

% Determine theta cycles (start/stop)
phases(:,2) = wrap(phases(:,2),1); % in [-pi,pi]
start = phases(ZeroCrossings(phases),1);
stop = start(2:end);
start = start(1:end-1);
data.time = start;

% Determine to which theta cycle each spike belongs
[~,spikeCycles] = InIntervals(spikes,[start stop]);
partialCycles = spikeCycles==0;
spikeCycles(partialCycles) = [];

% Compute number of spikes in each cycle
nSpikes = Accumulate(spikeCycles,1,size(start));

% Compute average phase in each cycle
spikePhases(partialCycles) = [];
data.phase = angle(Accumulate(spikeCycles,spikePhases,size(start))./nSpikes);

phase = Accumulate(spikeCycles,spikePhases,size(start))./nSpikes;
phase(nSpikes==0) = 0;
kernel = ones(7,1)/7;
data.phase = angle(conv(phase,kernel,'same'));

% Estimate spike rate and acceleration in each cycle

% Rate = average number of spikes over 7 surrounding cycles, which can be
% efficiently computed as the convolution of nSpikes with the kernel (1/7,..,1/7)
kernel = ones(7,1)/7;
data.rate = conv(nSpikes,kernel,'same');
% Acceleration = slope of linear regression over 7 surrounding cycles
% The general formula is b = (E[xy]-E[x]E[y]) / (E[x²]-E[x]²)
% To compute this efficiently, we first note that the slope is unchanged
% by translations along x the axis: we choose x in {-3..3}, so that E[x]=0.
% We then compute E[x²] = 4, and end up with b = 1/4 E[xy], which can be
% computed as the convolution of y with the kernel (3,2,1,0,-1,-2,-3)/7.
kernel = (3:-1:-3)/7;
data.acceleration = conv(nSpikes,kernel,'same')/4;

% Number of bins for x and y
nBinsX = nBins(1);
if length(nBins) == 1,
	nBinsY = nBinsX;
	nBins(2) = nBins;
else
	nBinsY = nBins(2);
end

% Min and max speed and acceleration
xLims = [0 limits(1)];
yLims = limits(2:3);

% Bin rate and acceleration
rate = Bin(data.rate,xLims,nBinsX);
acceleration = Bin(data.acceleration,yLims,nBinsY);

% Compute phase map
map.phase = angle(Smooth(Accumulate([rate acceleration],exp(j*data.phase),[nBinsX nBinsY]),smooth))';
map.count = Smooth(Accumulate([rate acceleration],1,[nBinsX nBinsY]),smooth)';

% Bins
map.x = linspace(xLims(1),xLims(2),nBinsX);
map.y = linspace(yLims(1),yLims(2),nBinsY);

%  figure;
%  hold on;
%  PlotXY(phases);
%  PlotTicks(spikes,'size',1,'k');
%  plot(start,nSpikes,'b*');
%  plot(start,data.rate,'r*');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the logistic
%  regression exercise. You will need to complete the following functions
%  in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the exam scores and the third column
%  contains the label.

data = load('ex2data1.txt');
X = data(:, [1, 2]); y = data(:, 3);

%% ==================== Part 1: Plotting ====================
%  We start the exercise by first plotting the data to understand the
%  the problem we are working with.

fprintf(['Plotting data with + indicating (y = 1) examples and o ' ...
         'indicating (y = 0) examples.\n']);

plotData(X, y);

% Put some labels
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============ Part 2: Compute Cost and Gradient ============
%  In this part of the exercise, you will implement the cost and gradient
%  for logistic regression. You neeed to complete the code in
%  costFunction.m

%  Setup the data matrix appropriately, and add ones for the intercept term
[m, n] = size(X);

% Add intercept term to x and X_test
X = [ones(m, 1) X];

% Initialize fitting parameters
initial_theta = zeros(n + 1, 1);

% Compute and display initial cost and gradient
[cost, grad] = costFunction(initial_theta, X, y);

fprintf('Cost at initial theta (zeros): %f\n', cost);
fprintf('Expected cost (approx): 0.693\n');
fprintf('Gradient at initial theta (zeros): \n');
fprintf(' %f \n', grad);
fprintf('Expected gradients (approx):\n -0.1000\n -12.0092\n -11.2628\n');

% Compute and display cost and gradient with non-zero theta
test_theta = [-24; 0.2; 0.2];
[cost, grad] = costFunction(test_theta, X, y);

fprintf('\nCost at test theta: %f\n', cost);
fprintf('Expected cost (approx): 0.218\n');
fprintf('Gradient at test theta: \n');
fprintf(' %f \n', grad);
fprintf('Expected gradients (approx):\n 0.043\n 2.566\n 2.647\n');

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 3: Optimizing using fminunc  =============
%  In this exercise, you will use a built-in function (fminunc) to find the
%  optimal parameters theta.

%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);

%  Run fminunc to obtain the optimal theta
%  This function will return theta and the cost
[theta, cost] = ...
	fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);

% Print theta to screen
fprintf('Cost at theta found by fminunc: %f\n', cost);
fprintf('Expected cost (approx): 0.203\n');
fprintf('theta: \n');
fprintf(' %f \n', theta);
fprintf('Expected theta (approx):\n');
fprintf(' -25.161\n 0.206\n 0.201\n');

% Plot Boundary
plotDecisionBoundary(theta, X, y);

% Put some labels
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============== Part 4: Predict and Accuracies ==============
%  After learning the parameters, you'll like to use it to predict the outcomes
%  on unseen data. In this part, you will use the logistic regression model
%  to predict the probability that a student with score 45 on exam 1 and
%  score 85 on exam 2 will be admitted.
%
%  Furthermore, you will compute the training and test set accuracies of
%  our model.
%
%  Your task is to complete the code in predict.m

%  Predict probability for a student with score 45 on exam 1
%  and score 85 on exam 2

prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n'], prob);
fprintf('Expected value: 0.775 +/- 0.002\n\n');

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);
fprintf('Expected accuracy (approx): 89.0\n');
fprintf('\n');

author: Ayatans

repo: Machine-Learning-homework


text:
function angle = Compute_Angle(xy_long)

x1 = xy_long(:, 1);
y1 = xy_long(:, 2);
K1 = -(y1(2)-y1(1))/(x1(2)-x1(1));
angle = atan(K1)*180/pi;

author: decouples


text:
%% subtri
% Below is a demonstration of the features of the |subtri| function

%% Syntax
% |[Fs,Vs]=subtri(F,V,n,uniqueOpt);|

%% Description
% The |subtri| function enables refinement of triangulated data

%% Examples

clear; close all; clc;

%%
% Plot Settings
fontSize=15;
faceAlpha=1;
edgeColor=0.2*ones(1,3);
edgeWidth=1.5;
markerSize=35;
markerSize2=20;

%% Refining a triangle

V=[0 0 0; 1 0 0; 0.5 sqrt(3)/2 0];
F=[1 2 3];

n=0:1:3; %Number of added edge nodes
pColors=gjet(numel(n));

cFigure;
for q=1:1:numel(n)
    [Fs,Vs]=subtri(F,V,n(q));
    subplot(2,2,q); hold on;
    title([num2str(n(q)),' added edge nodes'],'FontSize',fontSize);
    gpatch(Fs,Vs,pColors(q,:),'k');
    plotV(Vs,'k.','markerSize',markerSize2);
    plotV(V,'k.','markerSize',markerSize);

    axis equal; axis tight; view(2);
end
drawnow;

%% Refining a tetrahedron

[V,F]=platonic_solid(1,1);

n=0:1:3; %Number of added edge nodes
pColors=gjet(numel(n));

cFigure;
for q=1:1:numel(n)
    [Fs,Vs]=subtri(F,V,n(q));
    subplot(2,2,q); hold on;
    title([num2str(n(q)),' added edge nodes'],'FontSize',fontSize);
    gpatch(Fs,Vs,pColors(q,:),'k');
    plotV(Vs,'k.','markerSize',markerSize2);
    plotV(V,'k.','markerSize',markerSize);
    axisGeom(gca,fontSize);
    camlight headlight;
end
drawnow;

%% Refining triangulated surfaces in general

[F,V]=geoSphere(1,1);

n=[0 1 2 3]; %Number of added edge nodes
pColors=gjet(numel(n));

cFigure;
for q=1:1:numel(n)
    [Fs,Vs]=subtri(F,V,n(q));
    subplot(2,2,q); hold on;
    title([num2str(n(q)),' added edge nodes'],'FontSize',fontSize);
    gpatch(Fs,Vs,pColors(q,:),'k');
    axisGeom(gca,fontSize);
    camlight headlight;
end
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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:
function [foodMenu] = extractDietFromModel(model)
%This function analysis an VMH model and extracts the diet based on the
%lower bound constraints.
%
% USAGE:
%   [foodMenu] = extractDietFromModel(model)
% INPUT:
%   model:     A human model (WBM)
% OUTPUT:
%   foodMenu:   A nx2 cell array containing n dietary metabolites and their
%               corresponding flux
% Authors:
%       Bronson R. Weston, 2022

%Identify food reactions
foodRxns=find(contains(model.rxns,'Food_EX_'));
foodRxns=foodRxns(model.lb(foodRxns)<0);

%Identify dietary metabolite reactions
metRxns=find(contains(model.rxns,'Diet_EX_'));
metRxns=metRxns(model.lb(metRxns)<0);
foodFlux=[];
metFlux=[];
metItems={};
foodItems={};

%Excract lower bound constraints of food and diet reactions
if length(foodRxns)>0
    foodItems = model.rxns(foodRxns);
    foodFlux = -1*model.lb(foodRxns);
end
if length(metRxns)>0
    metItems = model.rxns(metRxns);
    metFlux = -1*model.lb(metRxns);
end

%Specify foodMenu
foodMenu=[foodItems,num2cell(foodFlux);metItems,num2cell(metFlux)];
end

author: opencobra


text:
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: copyright_cfg.m 3944 2010-06-23 08:53:40Z volkmar $

rev = '$Rev: 3944 $';