text:
function [sig,delay_offset] = delayline(sig,dt,weight,conf)
%DELAYLINE (fractional) delay line with weights
%
%   Usage: [sig,delay_offset] = delayline(sig,dt,weight,conf)
%
%   Input parameter:
%       sig     - input signal (vector), can be in the form of [N C], or
%                 [M C N], where
%                     N ... samples
%                     C ... channels (most probably 2)
%                     M ... number of measurements
%                 If the input is [M C N], the length of dt and weight has to be
%                 1 or M*C. In the last case the first M entries in dt are
%                 applied to the first channel and so on.
%       dt      - delay / s
%       weight  - amplitude weighting factor
%       conf    - configuration struct (see SFS_config).
%
%   Output parameter:
%       sig             - delayed signal
%       delay_offset    - additional delay / s
%                         This is added by the fractional delayline filters to
%                         all channels. For integer delays this is 0.
%
%   DELAYLINE(sig,dt,weight,conf) implementes a delayline, that delays the given
%   signal by dt samples and applies an amplitude weighting factor. The delay is
%   implemented as integer delay or fractional delay filter, see delayline
%   section in SFS_config for possible settings. As default setting an integer
%   delayline is used.
%
%   See also: get_ir, driving_function_imp_wfs

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Configuration ===================================================
% Check for old configuration
if isfield(conf, 'usefracdelay')
    error(['%s: conf.usefracdelay is deprecated, please use conf.delayline', ...
           ' instead. See SFS_config for details.'],upper(mfilename));
end
fs = conf.fs;
delay = conf.delayline;

%% ===== Preparation =====================================================
% --- Reshape signals ---
% Check if the signal is an impulse response given in SOFA conventions [M C N],
% or in usual [N C] convention, where
% M ... number of measurements
% C ... number of channels
% N ... number of samples
if ndims(sig)==3
    [M,C,samples] = size(sig);
    channels = M * C;
    % Reshape [M C N] => [N C*M], this will be redone at the end of the function
    sig = reshape(sig,[channels,samples])';
    reshaped = true;
else
    % Assume standard format [N C]
    [samples,channels] = size(sig);
    reshaped = false;
end

%% ===== Resampling ======================================================
% The resampling is applied independently from the actual fractional/integer
% delay handling performed in the next step. The resampling is redone at the end
% of the file.
% If resampling is used together with the integer delay filter, this is already
% a usage of fractional delay due to the upsampling.
%
switch delay.resampling
case 'none'
    rfactor = 1.0;
    delay_offset = 0;
case 'matlab'
    rfactor = delay.resamplingfactor;
    delay_offset = 0;
    sig = resample(sig,rfactor,1);
case 'pm'
    % === Parks-McClellan linear phase FIR filter ===
    rfactor = delay.resamplingfactor;
    rfilt = pm_filter(rfactor*delay.resamplingorder, 0.9/rfactor, ...
      1/rfactor);
    delay_offset = delay.resamplingorder*rfactor/2;

    sig = reshape(sig,1,channels*samples);
    sig = [sig; zeros(rfactor-1,channels*samples)];
    sig = reshape(sig,rfactor*samples,channels);

    sig = convolution(rfactor*rfilt, sig);
otherwise
    error('%s: "%s": unknown resampling method',upper(mfilename), ...
        delay.resampling);
end

%% ===== Expansion of signals, delays or weights =========================
% --- Expand channels
if channels==1
    channels = max(length(dt),length(weight));
    sig = repmat(sig,[1 channels]);
end

% --- Expand dt and weight ---
% If only single valued time delay and weight is given, create vectors
if channels>1 && length(dt)==1, dt=repmat(dt,[1 channels]); end
if channels>1 && length(weight)==1, weight=repmat(weight,[1 channels]); end

%% ===== Conversion to integer delay =====================================
dt = dt.*rfactor.*fs;  % resampled delay / samples
samples = rfactor.*samples;  % length of resampled signals
switch delay.filter
case 'integer'
    % === Integer delays ===
    idt = round(dt);  % round to nearest integer delay
    delay_offset = delay_offset + 0;
case 'zoh'
    % === Zero-order hold ===
    idt = ceil(dt);  % round to next larger integer delay
    delay_offset = delay_offset + 0;
case 'lagrange'
    % === Lagrange polynomial interpolator ===
    if iseven(delay.filterorder)
        idt = round(dt);  % round delay for even order
    else
        idt = floor(dt);  % floor delay for odd order
    end
    fdt = dt - idt;  % fractional part of delays
    b = lagrange_filter(delay.filterorder,fdt);
    a = ones(1,channels);
    delay_offset = delay_offset + floor(delay.filterorder/2);
case 'thiran'
    % === Thiran's allpass filter for maximally flat group delay ===
    idt = round(dt);  % integer part of delays
    fdt = dt - idt;  % fractional part of delays
    [b,a] = thiran_filter(delay.filterorder,fdt);
    delay_offset = delay_offset + delay.filterorder;
case 'least_squares'
    % ==== Least squares interpolation filter ===
    idt = floor(dt);  % integer part of delays
    fdt = dt - idt;  % fractional part of delays
    b = zeros(delay.filterorder+1,channels);
    for ii=1:channels
        b(:,ii) = general_least_squares(delay.filterorder+1,fdt(ii),0.90);
    end
    a = ones(1,channels);
    delay_offset = delay_offset + floor(delay.filterorder/2);
case 'farrow'
    % === Farrow-structure ===
    % Based on the assumption, that each coefficient h(n) of the fractional
    % delay filter can be expressed as a polynomial in d (frac. delay), i.e.
    %            __
    %           \  NPol
    % h_d(n) ~=  >      c_m(n) d^m
    %           /__m=0
    %
    % For some Filter design methods, e.g. Lagrange Interpolators, this is
    % perfectly possible. For other, a uniform grid of test delays d_q is
    % used to fit the polynomials to the desired coefficient(n) find a set
    % polynomial which approximates each coefficient of the desired filter.
    % This structure allows to perform the convolution independently from
    % the delay and reuse the results of the filter for different delays.
    %                           __
    %                          \  NPol
    % y(n) = h_d(n) * x(n) ~=   >      ( c_m(n)*x(n) ) d^m
    %                          /__m=0
    %
    % The above representation shows that the convolution of the input
    % signal x can be performed by first convolving c_m and x and
    % incorporating the delay d afterwards.
    %
    % number of parallel filters, i.e. order of polynomial + 1
    % Nfilter = delay.filternumber;
    to_be_implemented(mfilename);
otherwise
    error('%s: \"%s\" is an unknown delayline filter', ...
        upper(mfilename),delay.filter);
end
% Apply filter if needed
if exist('a','var') && exist('b','var')
    for ii=1:channels
        sig(:,ii) = filter(b(:,ii),a(:,ii),sig(:,ii));
    end
end

%% ===== Integer delayline ===============================================
% Handling of too long delay values (returns vector of zeros)
idt(abs(idt)>samples) = samples;
% Handle positive or negative delays
for ii=1:channels
    if idt(ii)>=0
        sig(:,ii) = [zeros(idt(ii),1); weight(ii)*sig(1:end-idt(ii),ii)];
    else
        sig(:,ii) = [weight(ii)*sig(-idt(ii)+1:end,ii); zeros(-idt(ii),1)];
    end
end

%% ===== Postprocessing ==================================================
% --- Downsampling ---
if rfactor~=1
    sig = sig(1:rfactor:samples,:);
    delay_offset = delay_offset ./ rfactor;
end
% --- Undo reshape ---
% [N M*C] => [M C N]
if reshaped
    % C might have changed due to replication of single-channel input
    sig = reshape(sig',M,[],size(sig,1));
end
% --- delay_offset in seconds ---
delay_offset = delay_offset / fs;

author: sfstoolbox

repo: sfs-matlab


text:
function linplus_test527 ( )

%*****************************************************************************80
%
%% TEST527 tests R8PP_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST527\n' );
  fprintf ( 1, '  R8PP_INDICATOR sets up a R8PP indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  a = r8pp_indicator ( n );

  r8pp_print ( n, a, '  The R8PP indicator matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = tangentspace(A, B)

% Get tangent space between A and B

for i = 1:size(A, 3)
    out(:, :, i) = sqrtm(B) * logm(B ^ (-1 / 2) * A(:, :, i) * B ^ (-1 / 2)) * sqrtm(B);
end

author: PatternRecognition

repo: OpenBMI


text:
function varargout = sound_field_mono_plane_wave(X,Y,Z,xs,f,conf)
%SOUND_FIELD_MONO_PLANE_WAVE sound field of a plane wave
%
%   Usage: [P,x,y,z] = sound_field_mono_plane_wave(X,Y,Z,xs,f,conf)
%
%   Input parameters:
%       X           - x-axis / m; single value or [xmin,xmax] or nD-array
%       Y           - y-axis / m; single value or [ymin,ymax] or nD-array
%       Z           - z-axis / m; single value or [zmin,zmax] or nD-array
%       xs          - direction of the plane wave
%       f           - monochromatic frequency / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       P           - Simulated sound field
%       x           - corresponding x values / m
%       y           - corresponding y values / m
%       z           - corresponding z values / m
%
%   SOUND_FIELD_MONO_PLANE_WAVE(X,Y,Z,xs,f,conf) simulates a monochromatic sound
%   field of a plane wave going in the direction xs for the frequency f.
%
%   To plot the result use:
%   plot_sound_field(P,X,Y,Z,conf);
%   or simple call the function without output argument:
%   sound_field_mono_plane_wave(X,Y,Z,xs,f,conf)
%
%   See also: sound_field_mono, plot_sound_field, sound_field_mono_point_source

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 6;
nargmax = 6;
narginchk(nargmin,nargmax);
isargxs(xs);
isargstruct(conf);

%% ===== Computation ====================================================
% Disable the plotting of a source, because we have a plane wave
conf.plot.loudspeakers = 0;
[varargout{1:nargout}] = sound_field_mono(X,Y,Z,[xs 0 1 0 1],'pw',1,f,conf);


text:
function x = circle01_to_circle ( alpha, r, c, n, x )

%*****************************************************************************80
%
%% CIRCLE01_TO_CIRCLE maps points from the unit circle to a general one.
%
%  Discussion:
%
%    To map data, defined in the unit circle, to a circle with radius R
%    and center C, and to rotate the original data by an angle of ALPHA first,
%
%      X = R * ( cos(alpha) -sin(alpha) ) * X + CX
%      Y       ( sin(alpha)  cos(alpha) )   Y   CY
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the angle, in radians, by which the unit circle
%    data should be rotated.
%
%    Input, real R, the scale factor by which the unit circle data
%    should be stretched.
%
%    Input, real C(2,1), the translation to be applied.
%
%    Input, integer N, the number of points to transform.
%
%    Input/output, real X(2,N), the points to be transformed.
%

%
%  Rotation:
%
%  Don't write "cos ( alpha )" instead of "cos(alpha)" because Matlab's
%  bebuggered automatic array creator thinks blanks separate matrix elements!
%
  rot = [ cos(alpha), -sin(alpha); ...
          sin(alpha),  cos(alpha) ];

  x(1:2,1:n) = c(1:2,1:2) * x(1:2,1:n);
%
%  Dilation:
%
  x(1:2,1:n) = r * x(1:2,1:n);
%
%  Translation:
%
  x(1:2,1:n) = x(1:2,1:n) + repmat ( c(1:2,1), 1, n );

  return
end


text:
function N = setupLaplace( dom )
%SETUPLAPLACE( DOM )  Construct a chebop2 object for Laplace operator.
% A small piece of code that is faster than calling the chebop2
% constructor for forming the Laplace operator on DOM.

    N = chebop2();
    N.domain = dom;
    N.op = @(u) lap(u);
    N.coeffs = [ 0 0 1 ; 0 0 0 ; 1 0 0 ];
    N.xorder = 2;
    N.yorder = 2;
    N.rhs = 0;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

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

l = r.c_prc.n_levels;

pvec(1:l)         = ptrans(1:l);                           % mu_0
pstruct.mu_0      = pvec(1:l);
pvec(l+1:2*l)     = exp(ptrans(l+1:2*l));                  % sa_0
pstruct.sa_0      = pvec(l+1:2*l);
pvec(2*l+1:3*l)   = ptrans(2*l+1:3*l);                     % rho
pstruct.rho       = pvec(2*l+1:3*l);
pvec(3*l+1:4*l-1) = exp(ptrans(3*l+1:4*l-1));              % ka
pstruct.ka        = pvec(3*l+1:4*l-1);
pvec(4*l:5*l-1)   = ptrans(4*l:5*l-1);                     % om
pstruct.om        = pvec(4*l:5*l-1);

return;

author: translationalneuromodeling

repo: tapas


text:
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
im = mksq(64, 1);
%im = zeros(64,64);
%im(32,:) = 1;
f1
idisp(im)
disp('create a square')
pause

f2
edge = icanny(im);
%edge = zeros(64,64);
%edge(32,:) = 1;
idisp(edge);
disp('take Laplacian to get 4 line segments')
pause

f3
h = Hough(edge)
disp('now show the Hough accumulator in which 4 peaks can be seen')
h.show()

f2
h.plot(4,'g')
h.peaks(4)
f1
h.plot(4,'g')

author: petercorke

repo: machinevision-toolbox-matlab


text:
classdef EventWithEncoderSpikes < Simulink.IntEnumType
    %EVENTWITHENCODERSPIKES Class for selecting a None, Warning or EmergencyStop signal
    %
    %   Class for selecting a None, Warning or EmergencyStop signal triggered by an encoder spike
    %   event.

    enumeration
        None(0)
        Warning(1)
        EmergencyStop(2)
    end
end

author: robotology

repo: whole-body-controllers


text:
%kkshad 'Compute the object surface shading giving a 3D cue '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kshad.pane file
%
% Parameters:
% InputFile: i1 'Z-buffer', required: 'Distance between view plane and object surface'
% OutputFile: o 'Shading', required: 'Shading information using Phong s model'
% Integer: background 'background', default: 255: 'Image background color (gray level)'
% Integer: Imin ' minimum ', default: 100: 'Minimum intensity of light that reaches the object surface '
% Integer: Imax ' maximum ', default: 244: 'Maximum intensity of light that reaches the object surface '
% InputFile: i2 'Normal Vectors', optional: 'Surface normal vectors'
% Double: Ka 'ambient ', default: 20: 'Coeficient of reflection for the ambient light'
% Double: Kd 'difuse ', default: 80: 'Coeficient of reflection for the difuse light'
% Double: Ks 'specular', default: 20: 'Coficient of reflection for the specular light'
% Integer: n 'n factor', default: 7: 'Falling factor for the specular light'
% Double: alpha 'alpha', default: 0: 'View plane rotation angle around Z axis'
% Double: beta 'beta ', default: 0: 'View plane rotation angle around X axis'
% InputFile: i3 '"Texture"', optional: '"Texture" information to be incorporated on shading'
% Double: Kt '"texture"', default: 0: 'Coficient of texture influence'
%
% Example: o = kkshad({i1, i2, i3}, {'i1','';'o','';'background',255;'Imin',100;'Imax',244;'i2','';'Ka',20;'Kd',80;'Ks',20;'n',7;'alpha',0;'beta',0;'i3','';'Kt',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kshad - Compute shading
%
%  DESCRIPTION
% Create a shaded image file based on z-buffer and/or normal vectors.
%
% The operator "kshad"
% creates a shaded image file of type KUBYTE from the z-buffer of a volume and
% optionally, also from the normal vectors of the volume.
%
% There are two
% shading possibilities: depth shading and gradient shading. The depth shading
% uses only the z-buffer information from the input file. The gradient shading
% uses the z-buffer and the angles between the normals and the light source
% direction (given by angles alpha and beta) to create the shaded image.
%
% The normal vectors can be computed in three different ways: in the image
% space, in the object space and in the voxel space. Each one will possibly
% give a different shaded image. These vectors must have their directions
% inverted (pointing towards the interior of the object).
%
% The input parameters Imax, Imin, Ka, Kd, Ks, and n, are chosen to
% adjust the gradient shading quality according to the Phong equation (see
% equation below).
% However, the depth shading only needs the Imax and Imin parameters. The
% image background color (gray level) also can be chosen by the user. If the
% texture information is to be added, then the texture file as well as the
% texture coeficient Kt must be given.
%
% Angles alpha and beta estabilish rotations of the light source around z and x
% axes respectively (see figure in man page of kzbuff), and they are
% recommended to be
% the same used to generate the z-buffer image (so that observer and light
% source have the same position), avoiding problems with shadowing, not
% implemented. Different values of alpha and beta can be used without
% restrictions when there is only one object and its visible surface is
% convex. The rotations are made firstly around x and then around z.
%
% Depth Shading Intensity (Id):
% Id(u,v) = [(Imax-Imin)/(dmin-dmax)] * [d(u,v)-dmax] + Imin
% where
% d corresponds to the distance in the z-buffer image
%
% Gradient Shading Intensity (Ig):
% Ig(u,v) = Imax*Ka + Id(u,v)*{kd*cos[theta(u,v)] + Ks*[cos(2*theta(u,v))^n] +
% Kt*T(u,v)}
% where
% theta is the angle between the surface nomal and the ray of light;
% Ka is the coeficient of reflection for the ambient light;
% Kd is the coeficient of diffuse reflection;
% Ks is the coeficient of specular reflection;
% Kt is the coeficient of texture;
% T is the optional texture information;
%
%
%
%  EXAMPLES
% kshad -i1 zbuffer.viff -i2 normals.viff -o image.viff -Ka 30 -Kd 80 -Ks 20
%       -n 7 -Imax 220 -Imin 100 -background 255 -alpha 0 -beta 0
%
% The program will create from the zbuffer.viff and normals.viff files an image
% represented by the
% image.viff file using the gradient shading. The gradient
% shading will use 30% (Ka) of reflection coeficient for the ambient light,
% 80% (Kd) of difuse reflection coeficient and 20% (Ks) of specular reflection
% coeficient for the light source. The maximum (Imax) and minimum (Imin)
% intensity of the light source reflected on the object surface will be
% respectively 220 and 100. The falling factor (n) for the specular reflection
% light component will be 7 and the image background color (gray level) will
% be 255. The zbuffer.viff image was generated by a view plane rotated zero
% degree around x axis and zero degree around z axis, therefore the position
% of the light source is given by setting angles alpha and beta to zero.
%
%  "SEE ALSO"
% kzbuff, kisnorm, kvsnorm, ktextu, kvoxext
%
%  RESTRICTIONS
% The input objects must have only the value segment.
%
% The z-buffer and texture
% input objects can not have dimention e > 1. The input object normal must have
% dimention e=3.
%
% The input objects must match in their dimentions w and h.
%
% The input objects can not be of data types KBIT and KCOMPLEX.
%
% In case of t > 1 in the input objects, the operator will be applied to the time
% t=0 only.
%
% None of the input and output objects are referenced, therefore some attributes
% may change, as the VALUE_POSITION, for example.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993, 1994, 1995 UNICAMP, R A Lotufo,  All rights reserved.
%

function varargout = kkshad(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkshad(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'background', 255;'Imin', 100;'Imax', 244;'i2', '__input';'Ka', 20;'Kd', 80;'Ks', 20;'n', 7;'alpha', 0;'beta', 0;'i3', '__input';'Kt', 0};
maxval={0,0,255,255,255,1,100,100,100,200,0,0,1,100};
minval={0,0,0,0,0,1,0,0,0,0,0,0,1,0};
istoggle=[0,0,1,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Integer','Integer','Integer','InputFile','Double','Double','Double','Integer','Double','Double','InputFile','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kshad"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef tflowmex < matlab.unittest.TestCase

  properties
    numLabs
  end

  properties (TestParameter)
    device = {'cpu', 'gpu'}
  end

  methods (TestClassSetup)
    function openParpool(test)
      test.numLabs = 4 ;
      if ~isempty(gcp('nocreate'))
        delete(gcp) ;
      end
      cl = parcluster('local');
      cl.NumWorkers = test.numLabs ;
      parpool(cl,test.numLabs) ;
    end
  end

  methods (TestClassTeardown)
    function closeParpool(test)
      if ~isempty(gcp('nocreate'))
        delete(gcp) ;
      end
    end
  end

  methods (TestMethodTeardown)
    function reset(test)
      if ~isempty(gcp('nocreate'))
        spmd
          vl_tflow('reset') ;
        end
      end
    end
  end

  methods (Test)
    function test_basic(test, device)
      haveGpu = strcmp(device,'gpu') ;
      T = 10 ;
      finalValue = (test.numLabs * (test.numLabs + 1)) / 2 * ...
          test.numLabs^(T-1) ;
      format = {'single', [10 10],     'x0', 'cpu', ;
                'single', [10 15 3 2], 'x1', 'cpu'  ;
                'double', [3 4 5 6],   'x2', 'cpu' } ;
      if haveGpu
        format = vertcat(format, ...
                         {'single', [10 10],     'x0g', 'gpu', ;
                          'single', [10 15 3 2], 'x1g', 'gpu'  ;
                          'double', [3 4 5 6],   'x2g', 'gpu' }) ;
      end
      spmd
        if haveGpu
          gpuDevice(1) ;
        end
        for i = 1:size(format,1)
          x{i} = tflowmex.makeArray(format(i,:)) + labindex ;
        end
        labBarrier() ;
        vl_tflow('init',format,labindex,numlabs) ;
        for t = 1:T
          for i = 1:size(format,1),
            vl_tflow('push',format{i,3},x{i}) ;
          end
          for i = 1:size(format,1),
            x{i} = vl_tflow('pull',format{i,3}) ;
          end
        end
        for i = 1:size(format,1)
          assert(x{i}(1) == finalValue) ;
        end
      end
    end

    function test_inplace(test, device)
      haveGpu = strcmp(device,'gpu') ;
      if ~haveGpu, return ; end
      T = 10 ;
      finalValue = (test.numLabs * (test.numLabs + 1)) / 2 * ...
          test.numLabs^(T-1) ;
      format = {'single', [10 10],     'x0', 'gpu', ;
                'single', [10 15 3 2], 'x1', 'gpu'  ;
                'double', [3 4 5 6],   'x2', 'gpu' } ;
      gpuDevice(1) ;
      spmd
        for i = 1:size(format,1)
          x{i} = tflowmex.makeArray(format(i,:)) + labindex ;
        end
        vl_tflow('init',format,labindex,numlabs) ;
        for t = 1:T
          for i = 1:size(format,1),
            vl_tflow('push',format{i,3},x{i},'inplace') ;
          end
          for i = 1:size(format,1),
            vl_tflow('pull',format{i,3},'inplace') ;
          end
        end
        for i = 1:size(format,1)
          assert(x{i}(1) == finalValue) ;
        end
      end
    end

  end

  methods (Static)
    function x = makeArray(format)
      x = zeros(format{1,2},format{1,1}) ;
      if numel(format) > 3 && strcmp(format{4},'gpu')
        x = gpuArray(x) ;
      end
    end
  end

end

repo: CREST-Release


text:
function channelClose(sockChannel)
	%
	try
		pause(1);
		sockChannel.close();
	catch
		%do nothing
		s	= lasterror;
		channelErrMsg('channelClose:badClose','could not close channel');
		channelErrMsg(s.identifier, s.message);
	end

	return;
end

author: Sable

repo: mcbench-benchmarks


text:
function [N,cN,cF] = per_corner_normals(V,F,varargin)
  % N = per_corner_normals(V,F,varargin)
  %
  % Inputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of mesh triangle indices into V
  %   Optional:
  %     'Angle'  followed by dihedral angle considered sharp. {pi*0.11}
  % Outputs:
  %   N  3*#F by 3 list of per-corner normals
  %   cN  #cN by 3 list of per-corner normals before expanding to corners
  %   cF  #F by 3 list of indices into cN
  %
  % Example:
  %   [V,F] = cylinder_mesh(1,20,'Caps',true);
  %   N = per_corner_normals(V,F);
  %   tsurf(F,V,'FaceColor','w',falpha(1,0.5));
  %   hold on;
  %   qvr(V(F,:),N);
  %   hold off;
  %   axis equal;
  %
  % Example:
  %   [V,F] = cylinder_mesh(1,20,'Caps',true);
  %   N = per_corner_normals(V,F);
  %   % zip into gltf style mesh (shared F, per-vertex attributes)
  %   [uVN,~,I] = unique([V(F,:) N],'rows');
  %   uV = uVN(:,1:3);
  %   uN = uVN(:,3+(1:3));
  %   uF = reshape(1:size(F,1)*3,size(F));
  %   uF = I(uF);
  %
  % Example:
  %   [V,F] = cylinder_mesh(1,20,'Caps',true);
  %   % zip into an obj style mesh (indexed V,N with separate faces)
  %   [~,cN,cF] = per_corner_normals(V,F);
  %   writeOBJ('per_corner_normals-test.obj',V,F,[],[],cN,cF);
  %
  %

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Angle'}, ...
    {'angle'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  E = sharp_edges(V,F);
  [cF,I] = cut_edges(F,E);
  cV = V(I,:);
  cN = per_vertex_normals(cV,cF);
  N = cN(cF,:);
end

author: alecjacobson

repo: gptoolbox


text:
function  gradH = lfmwhiteComputeGradThetaH2(gamma1, gamma2, t1, t2, ...
    gradTheta, isStationary1, isStationary2)

% LFMWHITECOMPUTEGRADTHETAH2 computes a portion of the LFM-WHITE kernel's gradient w.r.t. theta.
% FORMAT
% DESC Helper function for computing part of the gradient of
% the LFM-WHITE kernel w.r.t. to a generic parameter theta (mass, spring or
% damper). Used for obtaining the gradients w.r.t. parameters related to
% the second argument (gamma2).
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG gradTheta : gradient of gamma w.r.t. the generic parameter theta.
% ARG isStationary1: indicates whether the stationary version of the first
% kernel is used (TRUE) or not (FALSE). Set to FALSE by default.
% ARG isStationary2: indicates whether the stationary version of the second
% kernel is used (TRUE) or not (FALSE). Set to FALSE by default.
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : David Luengo, 2009
%
% SEEALSO : lfmwhiteKernParamInit, lfmwhiteXlfmwhiteKernGradient,
% lfmwhiteComputeH, lfmwhiteComputeGradThetaH1

% KERN

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end

if nargin < 6
    isStationary1 = false;
end
if nargin < 7
    isStationary2 = false;
end

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
deltaT = T1 - T2;
indT = deltaT >= 0;
psi = gamma1 * (1-indT) + gamma2 * indT;
gradH = -(lfmwhiteComputeH(gamma1, gamma2, t1, t2, isStationary1, isStationary2) ...
    + abs(deltaT) .* exp(-psi .* abs(deltaT)) .* indT);
if ((isStationary1 == false) | (isStationary2 == false))
    gradH = gradH + T1 .* exp(-(gamma2 * T1 * double(isStationary2 == false) ...
        + gamma1 * T2 * double(isStationary1 == false)));
end
gradH = gradH * gradTheta / (gamma1 + gamma2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function data = parse_json(string)
% DATA = PARSE_JSON(string)
% This function parses a JSON string and returns a cell array with the
% parsed data. JSON objects are converted to structures and JSON arrays are
% converted to cell arrays.

% F. Glineur, 2009
% (inspired by the JSON parser by Joel Feenstra on MATLAB File Exchange
% (http://www.mathworks.com/matlabcentral/fileexchange/20565) but with
% faster handling of strings)

% Copyright (c) 2011, François Glineur
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

assert(ischar(string), 'Input must be a string or a filename.');
if exist(string, 'file')
  string = fileread(string);
end

pos = 1;
len = length(string);
% String delimiters and escape characters are identified beforehand to improve speed
esc = regexp(string, '["\\]'); index_esc = 1; len_esc = length(esc);

if pos <= len
  switch(next_char)
    case '{'
      data = parse_object;
    case '['
      data = parse_array;
    otherwise
      error_pos('Outer level structure must be an object or an array');
  end
end

  function object = parse_object
    parse_char('{');
    object = [];
    if next_char ~= '}'
      while 1
        str = parse_string;
        if isempty(str)
          error_pos('Name of value at position %d cannot be empty');
        end
        parse_char(':');
        val = parse_value;
        object.(valid_field(str)) = val;
        if next_char == '}'
          break;
        end
        parse_char(',');
      end
    end
    parse_char('}');
  end

  function object = parse_array
    parse_char('[');
    object = cell(0, 1);
    if next_char ~= ']'
      while 1
        val = parse_value;
        object{end+1} = val;
        if next_char == ']'
          break;
        end
        parse_char(',');
      end
    end
    parse_char(']');
  end

  function parse_char(c)
    skip_whitespace;
    if pos > len || string(pos) ~= c
      error_pos(sprintf('Expected %c at position %%d', c));
    else
      pos = pos + 1;
      skip_whitespace;
    end
  end

  function c = next_char
    skip_whitespace;
    if pos > len
      c = [];
    else
      c = string(pos);
    end
  end

  function skip_whitespace
    while pos <= len && isspace(string(pos))
      pos = pos + 1;
    end
  end

  function str = parse_string
    if string(pos) ~= '"'
      error_pos('String starting with " expected at position %d');
    else
      pos = pos + 1;
    end
    str = '';
    while pos <= len
      while index_esc <= len_esc && esc(index_esc) < pos
        index_esc = index_esc + 1;
      end
      if index_esc > len_esc
        str = [str string(pos:end)];
        pos = len + 1;
        break;
      else
        str = [str string(pos:esc(index_esc)-1)];
        pos = esc(index_esc);
      end
      switch string(pos)
        case '"'
          pos = pos + 1;
          return;
        case '\'
          if pos+1 > len
            error_pos('End of file reached right after escape character');
          end
          pos = pos + 1;
          switch string(pos)
            case {'"' '\' '/'}
              str(end+1) = string(pos);
              pos = pos + 1;
            case {'b' 'f' 'n' 'r' 't'}
              str(end+1) = sprintf(['\' string(pos)]);
              pos = pos + 1;
            case 'u'
              if pos+4 > len
                error_pos('End of file reached in escaped unicode character');
              end
              str(end+1:end+6) = string(pos-1:pos+4);
              pos = pos + 5;
          end
        otherwise % should never happen
          str(end+1) = string(pos);
          pos = pos + 1;
      end
    end
    error_pos('End of file while expecting end of string');
  end

  function num = parse_number
    [num, one, err, delta] = sscanf(string(pos:min(len,pos+20)), '%f', 1); % TODO : compare with json(pos:end)
    if ~isempty(err)
      error_pos('Error reading number at position %d');
    end
    pos = pos + delta-1;
  end

  function val = parse_value
    switch(string(pos))
      case '"'
        val = parse_string;
        return;
      case '['
        val = parse_array;
        return;
      case '{'
        val = parse_object;
        return;
      case {'-','0','1','2','3','4','5','6','7','8','9'}
        val = parse_number;
        return;
      case 't'
        if pos+3 <= len && strcmpi(string(pos:pos+3), 'true')
          val = true;
          pos = pos + 4;
          return;
        end
      case 'f'
        if pos+4 <= len && strcmpi(string(pos:pos+4), 'false')
          val = false;
          pos = pos + 5;
          return;
        end
      case 'n'
        if pos+3 <= len && strcmpi(string(pos:pos+3), 'null')
          val = [];
          pos = pos + 4;
          return;
        end
    end
    error_pos('Value expected at position %d');
  end

  function error_pos(msg)
    poss = max(min([pos-15 pos-1 pos pos+20],len),1);
    if poss(3) == poss(2)
      poss(3:4) = poss(2)+[0 -1];         % display nothing after
    end
    msg = [sprintf(msg, pos) ' : ... ' string(poss(1):poss(2)) '<error>' string(poss(3):poss(4)) ' ... '];
    ME = MException('JSONparser:invalidFormat', msg);
    throw(ME);
  end

  function str = valid_field(str)
    % From MATLAB doc: field names must begin with a letter, which may be
    % followed by any combination of letters, digits, and underscores.
    % Invalid characters will be converted to underscores, and the prefix
    % "alpha_" will be added if first character is not a letter.
    if ~isletter(str(1))
      str = ['alpha_' str];
    end
    str(~isletter(str) & ~('0' <= str & str <= '9')) = '_';
  end

end

author: hpatches

repo: hpatches-benchmark


text:
% Fig. 9.12  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
% script to generate figure 9.12
% using the general simulation nonlin
N=1;
a=2;
r=0;
num=[1 2 1];
den=[1 0 0 0];
rlocus(num,den)
axis([-6 2 -3 3])
title('Root locus for the system of Figure 9.12')
 z=0:.1:.9;
 wn= 1:6;
 sgrid(z, wn)
 hold off


text:
function DataSet = prtDataGenMary
% prtDataGenMary  Generate unimodal M-ary example data
%
%  DATASET = prtDataGenMary returns a prtDataSetClass with randomly
%  generated data according to the following distribution.
%
%       H1: N([0 0],0.5*eye(2))
%       H2: N([0.5 0.5],0.1*eye(2))
%       H3: N([-2 -2],eye(2))
%
%  Example:
%
%  ds = prtDataGenMary;
%  plot(ds)
%
%   See also: prtDataSetClass, prtDataGenBiModal, prtDataGenIris,
%   prtDataGenMary, prtDataGenNoisySinc, prtDataGenOldFaithful,
%   prtDataGenSpiral, prtDataGenUnimodal, prtDataGenUnimodal, prtDataGenXor

rvH1 = prtRvMvn('mu',[0 0],'sigma',0.5*eye(2));
rvH2 = prtRvMvn('mu',[0.5 0.5],'sigma',0.1*eye(2));
rvH3 = prtRvMvn('mu',[-2 -2],'sigma',eye(2));
X = cat(1,draw(rvH1,100),draw(rvH2,100),draw(rvH3,100));
Y = prtUtilY(0,100,100,100);

DataSet = prtDataSetClass(X,Y,'name','prtDataGenMary');

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [ image_paths, image_set_name ] = get_image_paths_from_voc( voc_path, image_set, voc_year )

VOCopts         = initVOCOpts( voc_path, voc_year );
VOCopts.testset = image_set;
image_set_name  = ['voc_', voc_year, '_' image_set];

image_ext       = '.jpg';
image_dir       = fileparts(VOCopts.imgpath);
image_ids       = textread(sprintf(VOCopts.imgsetpath, image_set), '%s');
image_paths     = strcat([image_dir, filesep], image_ids, image_ext);
end

author: gidariss

repo: mrcnn-object-detection


text: %NOT Logical NOT. Datafile overload

author: marianux

repo: ecg-kit

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


text:
clc
clear

I=imread('peppers.png');
I=double(I);
class_number=3;
potential=0.5;
maxIter=30;
seg=ICM(I,class_number,potential,maxIter);
figure;
imshow(I);


text:
function [ml] = in32ml(in3)
% Convert volume from cubic inches to milliliters.
% Chad Greene 2012
ml = in3*16.387064;


text:
function S2G = HEALPixS2Grid(varargin)
% defines an equispaced spherical grid
%
% Syntax
%   equispacedS2Grid('points',300)
%   equispacedS2Grid('resolution',5*degree)
%
% Options
%  points     - number of points to be generated
%  resolution - resolution of the grid
%  hemisphere - 'lower', 'uper', 'complete', 'sphere', 'identified'
%  minRho     - starting rho angle (default 0)
%  maxRho     - maximum rho angle (default 2*pi)
%  minTheta   - starting theta angle (default 0)
%  maxTheta   - maximum theta angle (default pi)
%
% Flags
%  antipodal  - include <VectorsAxes.html antipodal symmetry>
%  restrict2MinMax - restrict margins to min / max
%  no_center  - ommit point at center
%
% See also
% regularS2Grid plotS2Grid

%  ind = rhoInside(rho,minrho,maxrho) & theta >= mintheta;
%  theta = theta(ind);
%  rho = rho(ind);
%
%  if isnumeric(maxtheta)
%    ind = theta <= maxtheta;
%  else
%    ind = theta <= maxtheta(rho);
%  end
%  theta = theta(ind);
%  rho = rho(ind);

res = get_option(varargin,'resolution',2.5*degree);

n = round(log(1/res)/log(2));

polar = HealpixGenerateSampling(2^n,1);

S2G = vector3d.byPolar(polar(:,1),polar(:,2));

return

% extract options
bounds = getPolarRange(varargin{:});

% get number of points
if check_option(varargin,'points') % calculate resolution
  ntheta = N2ntheta(fix(get_option(varargin,'points')),...
    bounds.VR{2},bounds.VR{4}-bounds.VR{3});%Check this
  res =  (bounds.VR{2}-bounds.VR{1}) / ntheta;
else
  res = get_option(varargin,'resolution',2.5*degree);
  res =  2* bounds.VR{2} / round(2 * bounds.VR{2} / res);
  ntheta = fix(round(2 * (bounds.VR{2}-bounds.VR{1}) / res + ...
    check_option(varargin,'no_center') )/2);
end

% define polar angle
if check_option(varargin,'no_center')
  theta =  (0.5:ntheta-0.5)*res;
else
  theta = (0:ntheta)*res;
end
theta = bounds.VR{1} + theta;

% define azimuth angles
identified = check_option(varargin,'antipodal');

rhGrid = repmat(S1Grid([],bounds.FR{3},bounds.FR{4} + pi,'periodic'),...
  1,length(theta));
for j = 1:length(theta)

  th = theta(j);
  if isappr(th,pi/2) && isappr(bounds.drho,2*pi) && identified
    rhGrid(j).max = bounds.FR{3} + pi;
    rhGrid(j).points = bounds.VR{3} + res*(0.5*mod(j,2)+(0:2*ntheta-1));
  else
    steps = max(round(sin(th) * bounds.drho / bounds.dtheta * ntheta),1);
    rhGrid(j).points = bounds.VR{3} + (0:steps-1 )* bounds.drho /steps + ...
      mod(j,2) * bounds.drho/steps/2;
  end
end

theta = S1Grid(theta,bounds.FR{1},bounds.FR{2});

if identified, opt = {'antipodal'}; else, opt = {}; end
S2G = S2Grid(theta,rhGrid,opt{:});
S2G = S2G.setOption('resolution',res);

% restrict to spherical region if specified
sR = getClass(varargin,'sphericalRegion');
if ~isempty(sR), S2G = S2G.subGrid(sR.checkInside(S2G,'noAntipodal')); end

end

% ---------------------------------------------------------
function ntheta = N2ntheta(N,maxtheta,maxrho)
ntheta = 1;
while calcAnz(ntheta,0,maxtheta,maxrho) < N
  ntheta = ntheta + 1;
end
if (calcAnz(ntheta,0,maxtheta,maxrho) - N) > (N-calcAnz(ntheta-1,0,maxtheta,maxrho))
  ntheta = ntheta-1;
end

end

function c = calcAnz(N,tmin,dt,dr)
c = sum(round(sin(tmin+dt/N*(1:N)) * dr/dt * N));
end

author: mtex-toolbox

repo: mtex

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


text:
function [fAValue] = calc_MaxLikelihoodACombined(mCatalog, mControl, fBValue, bReBin)
    % Maximum likelihood a-values for a combined catalog containing different periods while using a fixed b-value
    % [fAValue] = calc_MaxLikelihoodACombined(mCatalog, mControl, fBValue, bReBin)
    %
    %
    % Input parameters:
    %   mCatalog        Catalog containing different periods with varying
    %                   magnitude of completeness
    %   mControl        Controlmatrix containing informations about the single catalogs
    %                   mControl(n,:) contains information about caCatalogs{n}
    %                   Column 1: Starting time of catalog
    %                   Column 2: Magnitude of completeness
    %                   Column 3: Starting magnitude bin
    %                   Column 4: Magnitude bin stepsize (must be 0.1)
    %   fBValue         Fixed b-value
    %   bReBin          Rebin the catalogs to the given binning in mControl
    %                   (must be 0)
    %
    % Output parameters:
    %   fAValue         Maximum likelihood a-value
    %
    % Danijel Schorlemmer
    % July 17, 2002

    % Get the number of different periods in the catalog
    [nRow_, nColumn_] = size(mControl);
    % Init vector with ending times for each period
    vMaxTime_ = zeros(nRow_, 1);
    % Loop over the control matrix
    for nCnt_ = 1:nRow_
        % Determine starting time of period
        fMinTime_ = mControl(nCnt_, 1);
        % Determine ending time of period
        if nCnt_ < nRow_
            vMaxTime_(nCnt_) = mControl(nCnt_+1, 1);
        else
            vMaxTime_(nCnt_) = max(mCatalog.Date);
        end
        % Create subcatalog for period
        vSel_ = (mCatalog.Date >= fMinTime_) & (mCatalog.Date < vMaxTime_(nCnt_));
        mTmpCatalog_ = mCatalog.subset(vSel_);
        % Rebin the catalog (should not be used)
        if bReBin
            fMaxMag_ = max(mTmpCatalog_(:,6));
            for nBin_ = mControl(nCnt_,3):mControl(nCnt_,4):(fMaxMag_ + mControl(nCnt_,4))
                vSel_ = (mTmpCatalog_(:,6) >= (nBin_ - (mControl(nCnt_,4)/2))) & (mTmpCatalog_(:,6) < (nBin_ + (mControl(nCnt_,4)/2)));
                mTmpCatalog_(vSel_,6) = nBin_;
            end
        end
        % Cut the subcatalog at magnitude of completeness
        vSel_ = mTmpCatalog_(:,6) >= mControl(nCnt_,2);
        mTmpCatalog_ = mTmpCatalog_(vSel_,:);
        % Store the subcatalog
        caCatalogs_{nCnt_} = mTmpCatalog_;
    end
    % Add the ending times to the control matrix
    mControl = [mControl vMaxTime_];
    % Set the callback starting values
    vStartValue = 1;
    % Find the maximum likelihood solution
    [fAValue, ~, bExitFlag_] = fminsearch(@callback_LogLikelihoodACombined, vStartValue, [], caCatalogs_, mControl, fBValue);
end

author: CelsoReyes


text:
function f = filterobject(anything, cutoff, poles)
% FILTEROBJECT constructor for a filter object
%   f = filterobject() creates a new default filterobject
%   f = filterobject(filterobject) duplicates a filterobject
%   f = filterobject(type, cutoff, poles) create user-defined filterobject
%
%   TYPE:  'B' : Bandpass, 'H' : Highpass, 'L' : Lowpass
%   CUTOFF: [low, high] for bandpass, single value for others
%   POLES: number of poles used in the filter
%
%   default new object: Bandpass from 0.8 Hz to 5Hz, 2 poles (as dictated
%   by this constructor file

% VERSION: 1.0 of filter objects
% AUTHOR: Celso Reyes
% LASTUPDATE: 1/30/2007

load_global_namespace;

switch nargin
    case 0
            %create a fresh filterobject
            f.type = 'B';
            f.cutoff = [0.8 5];
            f.poles = 2;
            f = class(f, 'filterobject');

    case 1
        if isa(anything, 'filterobject')
            f = anything;
        else
            error(['trying to use filterobject(' class(anything) ')']);
        end

    case 3
        f = filterobject;
        f = set(f,'type',anything, 'cutoff', cutoff,'poles',poles);

    otherwise
        disp('Invalid arguments in filterobject constructor');
end;

%% LOAD filterobject's global namespace
% replaces SUITE_STUFF
function load_global_namespace()

persistent FILTER_NAMESPACE

if FILTER_NAMESPACE
    return
else
    FILTER_NAMESPACE = true;
end
%
% global FILTER_LONG
% global FILTER_STANDARD
%
% FILTER_LONG = filterobject('L',0.999,2);
% FILTER_STANDARD = filterobject('B', [0.8 5.0] ,2);

author: geoscience-community-codes

repo: GISMO


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

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

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = (((X*theta)-y)' * ((X*theta)-y))/(2*m);

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

end

author: UtkarshPathrabe


text:
function Goal_Idx = Select_Goal_Point(State,Params,Trajectory)

th = State(3);

if strcmp(Params.Sim.Direction,'reverse')
    Current_Direction = -[-sin(th); cos(th)];
    Dir = 2;
else
    Current_Direction = [-sin(th); cos(th)];
    Dir = 1;
end

Npts = Params.Traj.Npts;
Idx = 1:Npts;

%Determine the relative position vector for all points in Trajectory
Rel_Pos = zeros(2,Npts);
for i=1:Npts
    Rel_Pos(1,i) = Trajectory.States{Dir,i}(1) - State(1); %Relative X position
    Rel_Pos(2,i) = Trajectory.States{Dir,i}(2) - State(2); %Relative Y position
end

%[Index,Distance,Dot_Product]
Data = zeros(Npts,3);

Data(:,1) = Idx';
Data(:,2) = sqrt(Rel_Pos(1,:).^2 + Rel_Pos(2,:).^2);
Data(:,3) = Dot(Current_Direction,Rel_Pos)./Data(:,2)';   %==Cos(heading)

Wrong_Way = Data(:,3) < cos(Params.Ctl.Max_Target_Angle);
Too_Close = Data(:,2) < Params.Ctl.Minimum_Goal_Distance;
End_of_Road = Data(:,1) > (Params.Traj.Npts - Params.Traj.End_Buffer);

N_remove = sum(Wrong_Way | Too_Close | End_of_Road);

if N_remove < Npts   %Then at least one valid point remains
    Data(Wrong_Way | Too_Close | End_of_Road,:) = [];   %Remove invalid points
    Data_Sorted = sortrows(Data,2);
    %Select the closest valid goal point:
    Goal_Idx = Data_Sorted(1,1);
else
    N_backwards = sum(Wrong_Way);
    if N_backwards < Npts   %Then we can just pick a point that is too close
        %Data(Wrong_Way,:) = [];   %Remove invalid points
        %Data_Sorted = sortrows(Data,2);
        %Select the closest valid goal point:
        %Goal_Idx = Data_Sorted(1,1);
        Goal_Idx = -1;   %Delimiter to show that the road is ending
    else  %No valid points. Give up and select the closest
        %Data_Sorted = sortrows(Data,2);
        %Select the closest valid goal point:
        Goal_Idx = -1;   %Delimiter to show that the road is ended
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function test_failed=test_dgt
%TEST_DGT  Test DGT
%
%  This script runs a throrough test of the DGT routine,
%  testing it on a range of input parameters.
%
%  The computational backend is tested this way, but the
%  interface is not.
%
%  The script tests dgt, idgt, gabdual and gabtight.
%
%  Use TEST_WFAC and TEST_DGT_FAC for more specific testing
%  of the DGT backend.

Lr=[24,16,144,108,144,24,135,35,77,20];
ar=[ 4, 4,  9,  9, 12, 6,  9, 5, 7, 1];
Mr=[ 6, 8, 16, 12, 24, 8,  9, 7,11,20];

test_failed=0;

disp(' ===============  TEST_DGT ================');

disp('--- Used subroutines ---');

which comp_wfac
which comp_iwfac
which comp_sepdgt
which comp_isepdgt
which comp_sepdgtreal
which comp_isepdgtreal
which comp_gabdual_long
which comp_gabtight_long

for ii=1:length(Lr);

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  b=L/M;
  N=L/a;
  c=gcd(a,M);
  d=gcd(b,N);
  p=a/c;
  q=M/c;

  for rtype=1:2

    if rtype==1
      rname='REAL ';
      g=tester_rand(L,1);
    else
      rname='CMPLX';
      g=tester_crand(L,1);
    end;

    global LTFAT_TEST_TYPE;
    if strcmpi(LTFAT_TEST_TYPE,'single')
        C = gabframebounds(g,a,M);
        while C>1e3
%             warning(sprintf(['The frame is too badly conditioned '...
%                              'for single precision. Cond. num. %d. '...
%                              ' Trying again.'],C));

                         if rtype==1
                             rname='REAL ';
                             g=tester_rand(L,1);
                         else
                             rname='CMPLX';
                             g=tester_crand(L,1);
                         end;
                         C = gabframebounds(g,a,M);
        end
    end

    gd=gabdual(g,a,M);
    gt=gabtight(g,a,M);

    % --- Test windows against their reference implementations. ---

    ref_gd=ref_gabdual(g,a,M);
    res=norm(ref_gd-gd);
    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    fprintf(['REFDUAL %s L:%3i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i '...
               '%0.5g %s\n'],rname,L,a,b,c,d,p,q,res,fail);

    ref_gt=ref_gabtight(g,a,M);
    res=norm(ref_gt-gt);
    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['REFTIGHT %s L:%3i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i '...
               '%0.5g %s'],rname,L,a,b,c,d,p,q,res,fail);
    disp(s);

    % ---- Test gabdualnorm --------------------------------------
    res=gabdualnorm(g,gd,a,M);
    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['DUALNORM1 %s L:%3i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i '...
               '%0.5g %s'],rname,L,a,b,c,d,p,q,res,fail);
    disp(s);

    [o1,o2]=gabdualnorm(g,gd,a,M);
    res=o1-1+o2;
    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['DUALNORM2 %s L:%3i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i '...
               '%0.5g %s'],rname,L,a,b,c,d,p,q,res,fail);
    disp(s);

    for W=1:3

      if rtype==1
        f=tester_rand(L,W);
      else
        f=tester_crand(L,W);
      end;

      % --- Test DGT against its reference implementation. ---

      cc=dgt(f,g,a,M);
      cc2=ref_dgt(f,g,a,M);

      cdiff=cc-cast(cc2,class(cc));
      res=norm(cdiff(:));
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf(['REF %s L:%3i W:%2i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i '...
                 '%0.5g %s'],rname,L,W,a,b,c,d,p,q,res,fail);
      disp(s)

      % --- Test reconstruction of IDGT using a canonical dual window. ---

      r=idgt(cc,gd,a);
      res=norm(f-r,'fro');
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf(['REC %s L:%3i W:%2i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i ' ...
                 '%0.5g %s'],rname,L,W,a,b,c,d,p,q,res,fail);
      disp(s)

      % --- Test reconstruction of IDGT using a canonical tight window. ---

      res=norm(f-idgt(dgt(f,gt,a,M),gt,a),'fro');
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf(['TIG %s L:%3i W:%2i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i ' ...
                 '%0.5g %s'],rname,L,W,a,b,c,d,p,q,res,fail);
      disp(s);

      % Test the real valued transform
      if rtype==1

        % --- Reference test ---
        ccreal=dgtreal(f,g,a,M);
        M2=floor(M/2)+1;

        cdiff=cc(1:M2,:,:)-ccreal;
        res=norm(cdiff(:));
        [test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=sprintf(['REFREAL   L:%3i W:%2i a:%3i b:%3i c:%3i d:%3i p:%3i ' ...
                   'q:%3i %0.5g %s'],L,W,a,b,c,d,p,q,res,fail);
        disp(s);

        % --- Reconstruction test ---

        rreal=idgtreal(ccreal,gd,a,M);

        res=norm(f-rreal,'fro');
        [test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=sprintf(['RECREAL   L:%3i W:%2i a:%3i b:%3i c:%3i d:%3i p:%3i ' ...
                   'q:%3i %0.5g %s'],L,W,a,b,c,d,p,q,res,fail);
        disp(s)

      end;
    end;

  end;

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function kf = polynomial_correlation(xf, yf, a, b)
%POLYNOMIAL_CORRELATION Polynomial Kernel at all shifts, i.e. kernel correlation.
%   Evaluates a polynomial kernel with constant A and exponent B, for all
%   relative shifts between input images XF and YF, which must both be MxN.
%   They must also be periodic (ie., pre-processed with a cosine window).
%   The result is an MxN map of responses.
%
%   Inputs and output are all in the Fourier domain.
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/

	%cross-correlation term in Fourier domain
	xyf = xf .* conj(yf);
	xy = sum(real(ifft2(xyf)), 3);  %to spatial domain

	%calculate polynomial response for all positions, then go back to the
	%Fourier domain
	kf = fft2((xy / numel(xf) + a) .^ b);

end

repo: Context-Aware-CF-Tracking


text:
function value = p21_autonomous ( )

%*****************************************************************************80
%
%% P21_AUTONOMOUS reports whether problem P21 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [out, revertClass] = tofloat(in)
%TOFLOAT Converts an image to floating point.
%   [OUT, REVERTCLASS] = TOFLOAT(IN) converts the input image IN to
%   floating-point. If IN is a double or single image, then OUT
%   equals IN.  Otherwise, OUT equals IM2DOUBLE(IN).  REVERTCLASS is
%   a function handle that can be used to convert back to the class
%   of IN.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

identity = @(x) x;

todouble = @im2double;

converter_table = {'uint8',   todouble, @im2uint8
   'uint16',  todouble, @im2uint16
   'int16',   todouble, @im2int16
   'logical', todouble, @logical
   'double',  identity, identity
   'single',  identity, identity};

classIndex = find(strcmp(class(in), converter_table(:, 1)));

if isempty(classIndex)
   error('Unsupported input image class.');
end

out = converter_table{classIndex, 2}(in);

revertClass = converter_table{classIndex, 3};

author: dipum

repo: dipum-toolbox

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


text:
function [acc pred softpred] = testbisvm(model, X, labels, params)
vecstest = sparse(params.dictsize, length(X));
pred = zeros(1, length(X));
for i = 1:length(X)
    nbsentraw = X{i};
    nbsentraw=nbsentraw(nbsentraw>0);
    updates = unique(nbsentraw);
    vecstest(updates,i) = 1;
end

[pred, acc, softpred] = predict(double(labels'), vecstest', model);
%[predicted_labelf, accuracyf, decision_valuesf] = svmpredict(double(labels'), vecstest', modelf, '-t 2');
pred = pred';
softpred = softpred';

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
function g = std( f, varargin )
%STD   Standard deviation of a SEPARABLEAPPROX along one variable.
%   G = STD(F) returns the standard deviation of F in the y-variable (default).
%   That is, if F is defined on the rectangle [a,b] x [c,d] then
%
%                         d
%                        /
%     std(F)^2 = 1/(d-c) | ( F(x,y) - mean(F,1) )^2 dy
%                        /
%                        c
%
%   G = STD(F, FLAG, DIM) takes the standard deviation along the y-variable if
%   DIM = 1 and along the x-variable if DIM = 2. The FLAG is ignored and kept in
%   this function so the syntax agrees with the Matlab STD command.
%
% See also CHEBFUN/STD, SEPARABLEAPPROX/MEAN.

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

% Empty check:
if ( isempty( f ) )
    g = chebfun();
    return
end

dom = f.domain;
if ( nargin < 3 )
    dim = 1;   % default to std over the y-variable.
elseif ( nargin == 3 )
    dim = varargin{ 2 };
else
    error( 'CHEBFUN:SEPARABLEAPPROX:std:nargin', 'Too many input arguments.' );
end

if ( dim == 1 )          % y-variable.
    mx = chebfun2(@(x,y) feval( mean(f, 1).', x), dom);
    g = sqrt(1/(diff(dom(3:4))) * sum((f - mx).^2, 1)).' ;
elseif ( dim == 2 )      %  x-variable.
    my = chebfun2(@(x,y) feval( mean(f, 2), y), dom);
    g = sqrt(1/(diff(dom(1:2))) * sum((f - my).^2, 2));
else
    error('CHEBFUN:SEPARABLEAPPROX:std:dim', ...
        'Third argument should have value 1 or 2.');
end

end

author: chebfun

repo: chebfun


text:
% GAUSS3D - generate a 3-dimensional gaussian matrix
%
% Usage:
%   >> [ gaussmatrix ] = gauss2d( nX, nY, nZ);
%   >> [ gaussmatrix ] = gauss2d( nX, nY, nZ, ...
%                              sigmaX, sigmaY, sigmaZ, ...
%                              centerX, centerY, centerZ, mask)
%
% Example:
%   >> gauss3d(3,3,3); % generate a 3x3x3 gaussian matrix
%
% Inputs:
%   nX    - number of values in first dimension
%   nY    - number of values in second dimension
%   nZ    - number of values in third dimension
%   sigmaX  - width of standard deviation in first dim (default: nX/5)
%   sigmaY  - width of standard deviation in second dim (default: nY/5)
%   sigmaZ  - width of standard deviation in third dim (default: nZ/5)
%   centerX - location of center (default: nX/2)
%   centerY - location of center (default: nY/2)
%   centerZ - location of center (default: nZ/2)
%   mask    - (0->1) percentage of low values in the matrix to mask
%             with zeros (default: 0 or none)
%
% Output:
%   gaussmatrix - 3-D gaussian matrix
%
% Author: Arnaud Delorme, 2009

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

function mat = gauss3d( sizeX, sizeY, sizeZ, sigmaX, sigmaY, sigmaZ, meanX, meanY, meanZ, cut);

if nargin < 2
	help gauss2d
	return;
end
if nargin < 4
	sigmaX = sizeX/5;
end
if nargin < 5
	sigmaY = sizeY/5;
end
if nargin < 6
	sigmaZ = sizeZ/5;
end
if nargin < 7
	meanX = (sizeX+1)/2;
end
if nargin < 8
	meanY = (sizeY+1)/2;
end
if nargin < 9
	meanZ = (sizeZ+1)/2;
end
if nargin < 10
	cut = 0;
end

[X,Y,Z] = ndgrid(1:sizeX,1:sizeY,1:sizeZ);

mat = exp(-0.5*(  ((X-meanX)/sigmaX).*((X-meanX)/sigmaX)...
				 +((Y-meanY)/sigmaY).*((Y-meanY)/sigmaY)...
				 +((Z-meanZ)/sigmaZ).*((Z-meanZ)/sigmaZ)))...
            			/((sigmaX*sigmaY*sigmaZ)^(0.5)*pi);

if cut > 0
	maximum = max(mat(:))*cut;
	I = find(mat < maximum);
	mat(I) = 0;
end

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function xx = to_gpu(x)
    xx = gpuArray(single(x));
end

author: jimmy-ren


text:
function out=cot(x)

out=cos(x)./sin(x);

author: opencobra

repo: cobratoolbox


text:
%% COLLATZ_POOL uses the MATLABPOOL command to run the COLLATZ code.
%
%  Discussion:
%
%    Output printed by the function appears directly on the screen.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2010
%
%  Author:
%
%    John Burkardt
%
  matlabpool open local 4

  n = 10000000;

  tic
  j_max = collatz_fun ( n );
  toc

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Searched integers 1 through %d\n', n );
  fprintf ( 1, '  Longest Collatz sequence had length %d\n', j_max );

  matlabpool close


text:
function obj = updateQuiver(obj, quiverIndex)

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

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(quiverIndex).AssociatedAxis);

%-QUIVER DATA STRUCTURE- %
quiver_data = get(obj.State.Plot(quiverIndex).Handle);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

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

%-quiver xaxis-%
obj.data{quiverIndex}.xaxis = ['x' num2str(xsource)];

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

%-quiver yaxis-%
obj.data{quiverIndex}.yaxis = ['y' num2str(ysource)];

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

%-quiver type-%
obj.data{quiverIndex}.type = 'scatter';

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

%-quiver visible-%
obj.data{quiverIndex}.visible = strcmp(quiver_data.Visible,'on');

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

%-quiver mode-%
obj.data{quiverIndex}.mode = 'lines';

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

%-scatter name-%
obj.data{quiverIndex}.name = quiver_data.DisplayName;

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

%-quiver line color-%
col = 255*quiver_data.Color;
obj.data{quiverIndex}.line.color = ['rgb(' col(1) ',' col(2) ',' col(3) ')'];

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

% check for x/y vectors
if isvector(quiver_data.XData)
    [quiver_data.XData, quiver_data.YData] = meshgrid(quiver_data.XData,quiver_data.YData);
end

%-format data-%
xdata = reshape(quiver_data.XData,1,size(quiver_data.XData,1)*size(quiver_data.XData,1));
ydata = reshape(quiver_data.YData,1,size(quiver_data.YData,1)*size(quiver_data.YData,1));
udata = reshape(quiver_data.UData,1,size(quiver_data.UData,1)*size(quiver_data.UData,1));
vdata = reshape(quiver_data.VData,1,size(quiver_data.VData,1)*size(quiver_data.VData,1));

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

%-quiver x-%
m = 1;
for n = 1:length(xdata)
obj.data{quiverIndex}.x(m) = xdata(n);
obj.data{quiverIndex}.x(m+1) = xdata(n) + 0.1*udata(n);
obj.data{quiverIndex}.x(m+2) = nan;
m = m + 3;
end

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

%-quiver y-%
m = 1;
for n = 1:length(ydata)
obj.data{quiverIndex}.y(m) = ydata(n);
obj.data{quiverIndex}.y(m+1) = ydata(n) + 0.1*vdata(n);
obj.data{quiverIndex}.y(m+2) = nan;
m = m + 3;
end

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

%-quiver barbs-%
if isHG2()
    % 'MaxHeadSize' scalar, matlab clips to 0.2 in r2014b
    maxheadsize = quiver_data.MaxHeadSize;
    % barb angular width, not supported by matlab
    head_width = deg2rad(17.5);
    for n = 1:length(xdata)
        % length of arrow
        l = norm([0.1*udata(n), 0.1*vdata(n)]);

        % angle of arrow
        phi = atan2(vdata(n),udata(n));

        % make barb with specified angular width, length is prop. to arrow
        barb = [...
            [-maxheadsize*l*cos(head_width), maxheadsize*l*sin(head_width)]; ...
            [0, 0]; ...
            [-maxheadsize*l*cos(head_width), -maxheadsize*l*sin(head_width)]; ...
            [nan, nan]; ...
            ]';

        % affine matrix: rotate by arrow angle and translate to end of arrow
        barb_transformation = affine2d([...
            [cos(phi), sin(phi), 0]; ...
            [-sin(phi), cos(phi), 0]; ...
            [xdata(n) + 0.1*udata(n), ydata(n) + 0.1*vdata(n), 1];
            ]);

        % place barb at end of arrow
        barb = transformPointsForward(barb_transformation, barb')';

        % add barb to plot data, inserting is optimized in matlab >2010
        for col = 1:4
            obj.data{quiverIndex}.x(end+1) = barb(1,col); % point 1
            obj.data{quiverIndex}.y(end+1) = barb(2,col);
        end
    end
end

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

%-scatter showlegend-%
leg = get(quiver_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{quiverIndex}.showlegend = showleg;

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

end

author: goGPS-Project


function plotline(varargin)
%PLOTLINE Plot points specified by [x,y] pairs.
% Function PLOTLINE accepts an arbitrary number of
% [x,y] points and plots a line connecting them.
% In addition, it can accept a line specification
% string, and pass that string on to function plot.

% Define variables:
%   ii        -- Index variable
%   jj        -- Index variable
%   linespec  -- String defining plot characteristics
%   msg       -- Error message
%   varargin  -- Cell array containing input arguments
%   x         -- x values to plot
%   y         -- y values to plot
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/05/07    S. J. Chapman        Original code

% Check for a legal number of input arguments.
% We need at least 2 points to plot a line...
msg = nargchk(2,Inf,nargin);
error(msg);

% Initialize values
jj = 0;
linespec = '';

% Get the x and y values, making sure to save the line
% specification string, if one exists.
for ii = 1:nargin

   % Is this argument an [x,y] pair or the line
   % specification?
   if ischar(varargin{ii})

      % Save line specification
      linespec = varargin{ii};

   else

      % This is an [x,y] pair.  Recover the values.
      jj = jj + 1;
      x(jj) = varargin{ii}(1);
      y(jj) = varargin{ii}(2);

   end
end

% Plot function.
if isempty(linespec)
   plot(x,y);
else
   plot(x,y,linespec);
end


text:
function [percentAUCdecrease,varNames] = calcFeatureImportanceRF_HN_Old(pathRF,nPerms,nameOutcome,fSetName,seed)

startpath = pwd;

cd(pathRF), load('training'), load('testing')
if strcmp(nameOutcome,'DeathSign')
    outcome = testing.outcomes.Death;
else
    outcome = testing.outcomes.(nameOutcome);
end
nInst = numel(outcome); % Number of testing patients
perms = zeros(nInst,nPerms); rng(seed)
for p = 1:nPerms
    perms(:,p) = datasample(1:nInst,nInst,'Replace',false)';
end

results = load(['testResultsRF_',fSetName,'_',nameOutcome]); results = struct2cell(results); results = results{1};
AUC_baseline = results.AUC;
RF = load(['RF_',fSetName,'_',nameOutcome]); RF = struct2cell(RF); RF = RF{1};
indClinic = training.clinical.bestAdd.(nameOutcome).(fSetName);
text = testing.textures.(nameOutcome).(fSetName); nText = size(text,2);
tableTest = [text,testing.clinical.table(:,indClinic)];
varNames = tableTest.Properties.VariableNames'; nVar = numel(varNames);
percentAUCdecrease = zeros(nVar,1);
for v = 1:nVar
    tableTemp = tableTest;
    percentAUCdecrease_temp = 0;
    for p = 1:nPerms
        tableTemp(:,v) = tableTest(perms(:,p),v);
        [prob] = predictRF(tableTemp,RF);
        [AUC,~,~,~] = calcPerformMetrics(prob,outcome,0.5);
        percentAUCdecrease_temp = percentAUCdecrease_temp + (AUC - AUC_baseline)/AUC_baseline;
    end
    if percentAUCdecrease_temp > 0, percentAUCdecrease_temp = percentAUCdecrease_temp * -1; end
    percentAUCdecrease_temp = percentAUCdecrease_temp/nPerms;
    percentAUCdecrease(v) = percentAUCdecrease_temp;
end
[percentAUCdecrease,ind] = sort(percentAUCdecrease,'descend');
varNames = varNames(ind);

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [J,alpha,beta,T] = variationalDynLoreta(Ut,Y,s2,iLV,L,alpha,beta,options)

%[J,alpha,beta,T] = variationalDynLoreta(Ut,Y,s2,iLV,L,alpha,beta,options)
%
% Computes the posterior distribution of the parameters J given some data V.
% The program solves levels of inference: 1) optimization of parameters J, and
% 2) optimization of hyperparameters alpha and beta. See Trujillo-Barreto
% et. al. (2004) for details.
%
% Ut,s2, and iLV are defined as follows:
%     Y: Nsensors x time points data matrix
%     K: N x P predictor matrix
%     L: sparse P x P square root of the precision matrix
%     [U,s,V] = svd( K*inv(L) )
%     iLV = inv(L)*V
%     s2  = s.^2
%
% alpha, beta: hyperparameters
% J: estimated parapeters
%
%                     P(V|J,alpha)*P(J|beta)
% P(J|V,alpha,beta) = ----------------------
%                        P(V|alpha,beta)
%
% Author: Alejandro Ojeda, SCCN/INC/UCSD, Jan-2013
%
% References:
%   Trujillo-Barreto, N., Aubert-Vazquez, E., Valdes-Sosa, P.A., 2004.
%     Bayesian model averaging in EEG/MEG imaging. NeuroImage 21, 1300???1319

if nargin < 5, error('Not enough input arguments.');end
if nargin < 8
    options.maxTol = 1e-3;
    options.maxIter = 100;
    options.verbose = true;
    options.gridSize = 100;
end

s = s2.^(0.5);
n = length(s);
p = size(L,1);
ntp = size(Y,2);

% Initialize hyperparameters
if nargin < 8
    UtY = Ut*Y(:,1);
    tol = max([n p])*eps(max(s));
    lambda2 = logspace(log10(tol),log10(max(s)),options.gridSize);
    gcv = zeros(options.gridSize,1);
    parfor k=1:options.gridSize
        d = lambda2(k)./(s2+lambda2(k));
        f = diag(d)*UtY;
        gcv(k) = dot(f,f,1)/sum(d)^2;
    end
    loc = getMinima(gcv);
    if isempty(loc), loc = 1;end
    loc = loc(end);
    lambda2 = lambda2(loc);

    alpha = 0.001*(Y(:)'*Y(:))/n;
    beta = alpha*lambda2;
end
err = inf;

for it=1:options.maxIter
    if err < options.maxTol, break;end

    % computing statistics
    H = Ut'*diag(alpha.*s2./(alpha.*s2+beta))*Ut;
    SSE = mean( (Y - H*Y).^2 ,2);
    sigma2 = mean(SSE);
    kSjk = diag(H);

    q = diag(alpha.*s./(alpha.*s2+beta))*Ut;%*Y;
    Sj = mean((iLV*q).^2,2);
    aic = -2*log(sigma2) + 2*p;
    alpha_old = alpha;
    beta_old = beta;

    % updating hyperparameters
    alpha = updateAlpha(SSE,kSjk,ntp);
    beta  = updateBeta(Sj,p);

    err = 0.5*abs(sum(alpha_old-alpha)) + 0.5*abs(beta_old-beta);
    if options.verbose
        disp([num2str(it) ' => alpha: ' num2str(alpha_old) '  beta: ' num2str(beta_old) ' sse: ' num2str(mean(SSE)) ' hyrp. error: ' num2str(err) ' aic: ' num2str(aic)]);
    end
end
if it == options.maxIter, warning('Maximum iteration reached. Failed to converge.');end

% parameters's estimation
T = iLV*diag(alpha.*s./(alpha.*s2+beta))*Ut;
J = T*Y;

% standardized Loreta
E = sum(Y-H*Y,2);
sigma = E'*E/(n-trace(H));
dT = 1./sqrt(dot(T,T,2));
S = 1./sigma*dT;
S = S./std(eps+S);
T = bsxfun(@times,T,S);%sqrt(p)*
J = bsxfun(@times,J,S);%sqrt(p)*
end

%---
function indmin = getMinima(x)
fminor = diff(x)>=0;
fminor = ~fminor(1:end-1, :) & fminor(2:end, :);
fminor = [0; fminor; 0];
indmin = find(fminor);
end

function alpha = updateAlpha(SSE,kSjk,ntp,bp,cp)
if nargin < 4, bp = 3;end
if nargin < 5, cp = 3;end
b = 1./( (1/bp) + 0.5*SSE + 0.5*kSjk );
c = cp + ntp/2;
alpha = median(1./(b*c));
end

function beta = updateBeta(r,p,bp,cp)
if nargin < 3, bp = 3;end
if nargin < 4, cp = 3;end
b = 1./( 1/bp + sum(r) );
c = cp+p/2;
beta = 1/(b*c);
end

repo: REST


text:
function w = addhistory(w, whathappened,varargin)
   %ADDHISTORY function in charge of adding history to a waveform
   %   waveform = addhistory(waveform, whathappened);
   %   waveform = addhistory(waveform, formatString, [variables...])
   %
   %   Mostly, addhistory is internal to the waveform methods, and will be
   %   updated automatically when you change the waveform.  Things like
   %   filtering, taking a hilbert envelope, or detrending will automatically
   %   drop notes into the waveform's HISTORY field.
   %
   %   The second way of using addhistory follows the syntax of sprintf
   %
   %   Input Arguments
   %       WAVEFORM: a waveform object   N-DIMENSIONAL
   %       WHATHAPPENED: absolutely anything.  Really.
   %
   %   AddHistory appends not only what happened, but also keeps track of WHEN
   %   it happened.
   %
   %   Control of whether or not history is added automatically lies within
   %   the waveform constructor (in a global variable called WAVEFORM_HISTORY)
   %
   %   example
   %       w = waveform; %create a blank waveform
   %       w = addhistory(w,'the following procedures done by DoIt.m');
   %       N = 1; M = 'Today'
   %       w = addhistory(w,'this is sample #%d date:%s',N,M);
   %       % output: "this is sample #1 date:Today"
   %
   %
   % See also WAVEFORM/ADDFIELD, WAVEFORM/GET, WAVEFORM/HISTORY, SPRINTF

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   global WAVEFORM_HISTORY
   if nargin > 2,
      whathappened = sprintf(whathappened, varargin{:});
   end
   if WAVEFORM_HISTORY
      modtime = now;
      for N = 1 : numel(w);
         %History is stored in a cell, the format of which is [WHAT, WHEN]
         if isempty(w(N).history)
            w(N).history = {whathappened, modtime};
         else
            w(N).history = [w(N).history; {whathappened, modtime}];
         end
      end
   end
end


text:
% The COBRAToolbox: testSolveCobraMILP.m
%
% Purpose:
%     - testSolveCobraMILP tests the SolveCobraMILP function and its different methods
%
% Author:
%     - Original file: Joseph Kang 11/16/09
%     - CI integration: Laurent Heirendt, March 2017
%
% Note:
%       test is performed on objective as solution can vary between machines, solver version etc..

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSolveCobraMILP'));
cd(fileDir);

% save original solve
global CBT_MILP_SOLVER;
orig_solver = CBT_MILP_SOLVER;

% test solver packages

global SOLVERS

% Do this test for all available MIQP solvers
UseIfAvailable = fieldnames(SOLVERS);  % We will simply use all available solvers that are MIQP solvers.
solverPkgs = prepareTest('needsMILP', true, 'useSolversIfAvailable', UseIfAvailable);

% set the tolerance
tol = 1e-8;

for k = 1:length(solverPkgs.MILP)
    fprintf('   Running solveCobraMILP using %s ... ', solverPkgs.MILP{k});

    % change the COBRA solver (LP)
    solverOK = changeCobraSolver(solverPkgs.MILP{k}, 'MILP', 0);

    if solverOK
        % MILP Solver test: chemeng.ed.ac.uk/~jwp/MSO/section5/milp.html

        % set up MILP problem.
        MILPproblem.c = [20; 6; 8];
        MILPproblem.A = [0.8, 0.2, 0.3;
                         0.4, 0.3, 0;
                         0.2, 0, 0.1];
        MILPproblem.b = [20; 10; 5];
        MILPproblem.lb = [0; 0; 0];
        MILPproblem.ub = [1000; 1000; 1000];
        MILPproblem.osense = -1;
        MILPproblem.csense = ['L'; 'L'; 'L'];
        MILPproblem.vartype = ['I'; 'I'; 'I'];
        MILPproblem.x0 = [0, 0, 0];
        pass = 1;

        % solve MILP problem setting the relative MIP gap tolerance and integrality tolerance to 1e-12 using parameters structure.
        if strcmp(solverPkgs.MILP{k}, 'cplex_direct') || strcmp(solverPkgs.MILP{k}, 'tomlab_cplex')
            parameters.relMipGapTol = 1e-12;
            parameters.intTol = 1e-12;
            MILPsolution = solveCobraMILP(MILPproblem, parameters);
            % check if MILP can be solved without x0 supplied
            MILPsolution2 = solveCobraMILP(rmfield(MILPproblem, 'x0'), parameters);
        else
            MILPsolution = solveCobraMILP(MILPproblem);
            % check if MILP can be solved without x0 supplied
            MILPsolution2 = solveCobraMILP(rmfield(MILPproblem, 'x0'));
        end

        % check results with expected answer.
        assert(all(abs(MILPsolution.int - [0; 31; 46]) < tol))
        assert(abs(MILPsolution.obj - 554) < tol)
        assert(abs(MILPsolution2.obj - 554) < tol)

        if strcmp(solverPkgs.MILP{k}, 'ibm_cplex')
            % test IBM-Cplex-specific parameters. Solve with the below parameters changed
            cplexParams = struct();
            cplexParams.emphasis.mip = 0;  % MIP emphasis: balance optimality and integer feasibility
            cplexParams.mip.strategy.search = 2;  % MIP search method: dynamic search
            MILPsolution = solveCobraMILP(MILPproblem, cplexParams, 'logFile', 'testIBMcplexMILPparam1.log');
            % check expected answer
            assert(all(abs(MILPsolution.int - [0; 31; 46]) < tol))
            assert(abs(MILPsolution.obj - 554) < tol)

            % solve with the parameters changed to other values
            cplexParams.emphasis.mip = 1;  % MIP emphasis: integer feasibility.
            cplexParams.mip.strategy.search = 1;  % MIP search method: traditional branch-and-cut search
            MILPsolution = solveCobraMILP(MILPproblem, cplexParams, 'logFile', 'testIBMcplexMILPparam2.log');
            % check expected answer
            assert(all(abs(MILPsolution.int - [0; 31; 46]) < tol))
            assert(abs(MILPsolution.obj - 554) < tol)

            % compare the log files to see whether the parameter changes are implemented
            testLog = {''; ''};
            paramsInLog = cell(2, 1);
            % text that should be found during the first test
            paramsInLog{1} = {'balance optimality and feasibility'; 'dynamic search'};
            % text that should be found during the second test
            paramsInLog{2} = {'integer feasibility'; 'branch-and-cut'};
            for jTest = 1:2
                % read the log files
                f = fopen(['testIBMcplexMILPparam' num2str(jTest) '.log'], 'r');
                l = fgets(f);
                while ~isequal(l, -1)
                    testLog{jTest} = [testLog{jTest}, l];
                    l = fgets(f);
                end
                fclose(f);
                % check that the expected parameter values are set, and the unexpected are not set.
                for p = 1:2
                    assert(~isempty(strfind(testLog{jTest}, paramsInLog{jTest}{p})));
                    assert(isempty(strfind(testLog{jTest}, paramsInLog{setdiff(1:2, jTest)}{p})));
                end
                % delete the log files
                delete(['testIBMcplexMILPparam' num2str(jTest) '.log']);
            end
            fprintf('Test ibm_cplex output to command window ...\n')
            % solve without logToFile = 1
            diary test_ibm_cplex_output_to_console1.txt
            sol = solveCobraMILP(MILPproblem);
            diary off
            % read the diary, which should be empty
            f = fopen('test_ibm_cplex_output_to_console1.txt', 'r');
            l = fgets(f);
            assert(isequal(l, -1))
            fclose(f);
            delete('test_ibm_cplex_output_to_console1.txt')

            % solve wit logToFile = 1
            diary test_ibm_cplex_output_to_console2.txt
            sol = solveCobraMILP(MILPproblem, 'logFile', 1);
            diary off
            % read the diary, which should be non-empty
            f = fopen('test_ibm_cplex_output_to_console2.txt', 'r');
            l = fgets(f);
            line = 0;
            while ~isequal(l, -1)
                line = line + 1;
                l = fgets(f);
            end
            fclose(f);
            assert(line > 3)
            delete('test_ibm_cplex_output_to_console2.txt')
            fprintf('Test ibm_cplex output to command window ... Done\n')

        end

        if strcmp(solverPkgs.MILP{k}, 'gurobi')
            % check additional parameters for Gurobi
            % temporarily shut down warning
            warning_stat = warning;
            warning off
            MILPproblem = struct();
            MILPproblem.A = [speye(10, 20), -3 * rand(10, 30)];
            MILPproblem.b = zeros(10, 1);
            MILPproblem.c = ones(50, 1);
            MILPproblem.lb = [-1000 * ones(35, 1); zeros(15, 1)];
            MILPproblem.ub = [1000 * ones(35, 1); ones(15, 1)];
            MILPproblem.vartype = char(['C' * ones(1, 20), 'I' * ones(1, 15), 'B' * ones(1, 15)]);
            MILPproblem.csense = char('E' * ones(1, 10));
            MILPproblem.osense = -1;
            % test TimeLimit as a gurobi-specific parameter
            sol = solveCobraMILP(MILPproblem, struct('TimeLimit', 0));
            assert(strcmp(sol.origStat, 'TIME_LIMIT'))
            % restore previous warning state
            warning(warning_stat)

            % check user-supplied x0
            MILPproblem.A = rand(10, 20);
            MILPproblem.b = 1000 * ones(10, 1);
            MILPproblem.c = zeros(20, 1);
            MILPproblem.lb = zeros(20, 1);
            MILPproblem.ub = ones(20, 1);
            MILPproblem.vartype = char(['C' * ones(1, 10), 'B' * ones(1, 10)]);
            MILPproblem.csense = char('L' * ones(1, 10));

            % no objective function. The supplied should be the returned
            % (if not everything becomes zero after presolve)
            MILPproblem.x0 = zeros(20, 1);
            sol = solveCobraMILP(MILPproblem);
            assert(isequal(sol.full, MILPproblem.x0));

            MILPproblem.x0 = ones(20, 1);
            sol = solveCobraMILP(MILPproblem);
            assert(isequal(sol.full, MILPproblem.x0));

        end
    end

    % output a success message
    fprintf('Done.\n');
end

% test ibm_cplex output to command window
solverOK = changeCobraSolver('ibm_cplex', 'MILP', 0);
if solverOK
    fprintf('Test ibm_cplex output to command window ...\n')
    % solve without logToFile = 1
    diary test_ibm_cplex_output_to_console1.txt
    sol = solveCobraMILP(MILPproblem);
    diary off
    % read the diary, which should be empty
    f = fopen('test_ibm_cplex_output_to_console1.txt', 'r');
    l = fgets(f);
    assert(isequal(l, -1))
    fclose(f);
    delete('test_ibm_cplex_output_to_console1.txt')

    % solve wit logToFile = 1
    diary test_ibm_cplex_output_to_console2.txt
    sol = solveCobraMILP(MILPproblem, 'logFile', 1);
    diary off
    % read the diary, which should be non-empty
    f = fopen('test_ibm_cplex_output_to_console2.txt', 'r');
    l = fgets(f);
    line = 0;
    while ~isequal(l, -1)
        line = line + 1;
        l = fgets(f);
    end
    fclose(f);
    assert(line > 3)
    delete('test_ibm_cplex_output_to_console2.txt')
    fprintf('Test ibm_cplex output to command window ... Done\n')
end

% remove the generated file
fullFileNamePath = [fileparts(which(mfilename)), filesep, 'MILPProblem.mat'];
if exist(fullFileNamePath, 'file') == 2
    delete(fullFileNamePath);
end

% change back to the original solver
if ~isempty(orig_solver)
    changeCobraSolver(orig_solver, 'MILP', 0);
end

% change the directory
cd(currentDir)


text:
function [ref_signal, refsig_fs, ref_signal_ds] = ...
    read_signal(sample_name, num_devices)

    global GYRO_FS;
    [ref_signal, refsig_fs] = audioread(['samples/' sample_name '.wav']);
    ref_signal_ds = resample(ref_signal, GYRO_FS * num_devices, refsig_fs);
    ref_signal_ds = normalization(ref_signal_ds);
    display(['Samples in ' sample_name ' (resampled): ' ...
        num2str(length(ref_signal_ds))]);
    figure;
    subplot(121);
    fft_plot(ref_signal_ds, GYRO_FS * num_devices);
    title([sample_name ' FFT']);

end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
names_list = { ...
%    {'red-9' 'blue-7'}, ...
%    {'blue-1' 'red-5'}, ...
%    {'blue-5' 'red-3'}, ...
%    {'red-7' 'yellow-1'}, ...
%    {'red-5' 'blue-7'}, ...
%    {'yellow-2' 'purple-1'}, ...
%    {'green-2' 'red-5'}, ...  % PBM
%    {'green-2' 'purple-1'}, ...
% {'street-3', 'street-5'}, ...
% {'nature-3' 'art'}, ...
% {'street-4' 'street-6'}, ...
% {'street-2' 'street-7'}, ...
% {'street-8' 'street-1'}, ...
% {'nature-7' 'street-4'}, ...
% {'nature-7' 'street-2'}, ...
% {'street-11' 'street-10'}, ...
% {'colors1' 'colors2'}, ...
% {'colors3' 'colors4'}, ...
% {'colors5' 'colors6'}, ...
% {'colors7' 'colors8'}, ...
% {'colors9' 'colors8'}, ...
% {'nature-7' 'street-1'}, ...
% {'street-2' 'street-4'}, ...
% {'colors1' 'street-1'}, ...
% {'colors2' 'street-2'}, ...
% {'colors3' 'street-3'}, ...
% {'colors4' 'nature-7'}, ...
% {'colors5' 'street-4'}, ...
% {'colors6' 'street-5'}, ...
% {'colors7' 'street-7'}, ...
% {'colors8' 'street-10'}, ...
% {'colors9' 'street-11'}, ...
% {'colors10' 'vegetable-9'}, ...
% {'nature-2' 'portrait-2'}, ...
% {'nature-3' 'vegetable-1'}, ...
% {'nature-4' 'vegetable-6'}, ...
% {'nature-5' 'vegetable-11'}, ...
% {'nature-6' 'vegetable-12'}, ...
 {'green-5' 'red-5'}, ...
 {'green-6' 'red-5'}, ...
 {'green-7' 'red-5'}, ...
 {'landscape-1' 'red-5'}, ...
    };

for it=1:length(names_list)

    names = names_list{it};
    disp(['------ ' names{1} '/' names{2} ' ------']);
    generateColorFigure;
end

author: gpeyre