text:
function [D, X] = ODL(Y, k, lambda, opts, method)
% * Solving the following problem:
%  (D, X) = \arg\min_{D,X} 0.5||Y - DX||_F^2 + lambda||X||_1
% * Syntax: `[D, X] = ODL(Y, k, lambda, opts, sc_method)`
%   - INPUT:
%     + `Y`: collection of samples.4/7/2016 7:35:39 PM
%     + `k`: number of atoms in the desired dictionary.
%     + `lambda`: norm 1 regularization parameter.
%     + `opts`: option.
%     + `sc_method`: sparse coding method used in the sparse coefficient update. Possible values:
%       * `'fista'`: using FISTA algorithm. See also [`fista`](#fista).
%       * `'spams'`: using SPAMS toolbox [[12]](#fn_spams).
%   - OUTPUT:
%     + `D, X`: as in the problem.
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 4/7/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
	if nargin == 0
		addpath(fullfile('..', 'utils'));
        addpath(fullfile('..', 'build_spams'));
        d      = 50; % data dimension
        N      = 100; % number of samples
        k      = 50; % dictionary size
        lambda = 0.1;
        Y      = normc(rand(d, N));
        %
		opts.max_iter      = 500;
		opts.show_progress = 0;
		opts.check_grad    = false;
		opts.tol           = 1e-8;
		opts.verbose     = true;

	end
	%%
	opts = initOpts(opts);
	%%
	%% ========= initial D ==============================
	D = PickDfromY(Y, [0, size(Y,2)], k);
    X = zeros(size(D,2), size(Y,2));
    if opts.verbose
        fprintf('cost: %f', ODL_cost(Y, D, X, lambda));
    end
    optsX = opts;
	optsX.max_iter = 200;
	optsX.tol      = 1e-8;
    optsD = opts;
	optsD.max_iter = 200;
	optsD.tol      = 1e-8;
	iter = 0;
	while iter < opts.max_iter
		iter = iter + 1;
		%% ========= sparse coding step ==============================
		X = lasso_fista(Y, D, X, lambda, optsX);
       	if opts.verbose
			costX = ODL_cost(Y, D, X, lambda);
			fprintf('iter: %3d, costX = %5f\n', iter, costX)
		end
		%% ========= dictionary update step ==============================
		F = X*X'; E = Y*X';
		D = ODL_updateD(D, E, F, optsD);
		if opts.verbose
			costD = ODL_cost(Y, D, X, lambda);
			fprintf('iter: %3d, costD = %5f\n', iter, costD)
		end
	end
	%%
	if nargin == 0
		pause;
	end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function hf = symmetrize_filter(hf)

% Ensure hermetian symmetry.

for k = 1:length(hf)
    dc_ind = (size(hf{k},1) + 1) / 2;
    hf{k}(dc_ind+1:end,end,:) = conj(flipud(hf{k}(1:dc_ind-1,end,:)));
end

author: martin-danelljan

repo: ECO


text:
function z=convfft(x,h,d,m,h0,x1,x2)
%CONFFT 1-D convolution or correlation using FFT
%
%  Usage: (1) z=convfft(x,h,d1,'',1,1,size(x,d)+length(h)-1);            % equivalent to conv(x,h)
%
%         (2) hh=convfft(length(x),h,2,'z',1,1,length(x)+length(h)-1);   % precalculate fft(h) for multiple calls
%             z=convfft(x,hh);                                           % also equivalent to conv(x,h)
%
%         (3) z=convfft(x,h);                                            % equivalent to filter(h,1,x)
%             z=convfft(x,h,d1,'',1,1,size(x,d));                        % equivalent to filter(h,1,x)
%
%         (4) z=convfft(x,h,d,'X',1,2-max(size(x,d),length(h)),max(size(x,d),length(h)));  % equivalent to xcorr(x,h)
%
%         (5) z=convfft(x,h,d,'x',floor((1+length(h))/2),1,size(x,d));   % equivalent to imfilter(x,h)
%
%  Inputs: x     input vector or array (or size(x,d) for the 'z' option)
%          h     impulse response (or z output from previous call with the 'z' option)
%          d     dimension of x to do convolution along [first non-singleton]
%          m     mode options (see below)
%          h0    origin sample number in h (can be outside the range 1:length(h)) [default: 1]
%          x1,x2 range of rows/columns in x to align with origin of h (can be outside the range 1:size(x,d)) [default: 1,size(x,d)]
%
% Outputs: z     output from convolution/correlation. The same size and shape as x except that size(z,d)=x2-x1+1
%                If m='z' is specified, then z is a structure that can be used as h in a subsequent call
%
% Mode is any sensible combination of the following
%          'x'   perform real correlation rather than convolution (reflects h around sample h0)
%          'X'   perform complex correlation rather than convolution (reflects and conjugates h around sample h0)
%          'z'   Precalculate fft(h) for efficiency. Input d must be given explicitly and input x equals size(x,d).

%      Copyright (C) Mike Brookes 2016-2017
%      Version: $Id: convfft.m 10118 2017-09-17 19:45:51Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isstruct(h) % normal input calling sequence
    if nargin<4 || isempty(m)
        m='';                       % set default mode string
    end
    if nargin<5 || isempty(h0)
        h0=1;                       % set default h origin
    end
    s=size(x);                      % get the structure of x
    ps=numel(x);                  	% number of elements in x
    ns=length(s);                   % number of dimensions in x
    if any(m=='z')      % output pre-computed structure instead of convolution
        if nargin<3 || isempty(d)   % if d unspecified
            error('d must be specified explicitly');
        end
        if numel(x)~=1              % x contains nx for a future call
            error('x must equal size(*,d)');
        end
        nx=x;                       % save x as the value of nx
    else
        if nargin<3 || isempty(d)       % if d unspecified
            if ps<2
                d=1;                    % if d is a singleton or is empty
            else
                d=find(s>1,1);          % d = first nonsingleton dimension
            end
        end
        nx=s(d);                        % length in correlation dimension
    end
    k=ps/nx;                        % total size of all other dimensions
    if nargin<6 || isempty(x1)
        x1=1;                       % default initial lag
    end
    if nargin<7 || isempty(x2)
        x2=nx;                      % default final lag
    end
    nz=x2-x1+1;                     % number of output lags required
    h=h(:);                         % force h to be a column
    nh=length(h);
    if any(m=='X')                  % do complex correlation rather than convolution
        h=conj(h(nh:-1:1));         % reflect and conjugate h
        h0=nh+1-h0;                 % reflect the position of h0
    elseif any(m=='x')              % do real correlation
        h=h(nh:-1:1);               % reflect h
        h0=nh+1-h0;                 % reflect the position of h0
    end
    hmin=h0+x1-nx;                  % smallest h index ever used
    hmax=h0+x2-1;                   % largest h index ever used
    xmin=x1+h0-nh;                  % smallest x index ever used
    xmax=x2+h0-1;                   % largest x index ever used
    if hmin>1 || hmax<nh            % we can delete some unused h values
        hmin=max(hmin,1);
        h=h(hmin:min(hmax,nh));     % trim h if possible
        nh=length(h);
        h0=h0-hmin+1;               % update h0 to new position
    end
    if xmin>1 || xmax<nx            % we can delete some unused x values
        vmin=max(xmin,1);           % we will trim v to v(vmin:vmax)
        vmax=min(xmax,nx);
        x1=x1-vmin+h0;              % update x1,x2 to new positions assuming h0=0
        x2=x2-vmin+h0;
    else
        vmin=1;
        vmax=nx;
        x1=x1+h0-1;                 % update x1,x2 to new positions assuming h0=0
        x2=x2+h0-1;
    end
    nv=vmax-vmin+1;                 % number of elements of v to retain
    nxz=min(max(max(nh-x1,0),max(x2-nv,0)),nh-1); % number of zeros to add to v is the larger of the number needed at each end
    [fnx,enx]=log2(nv+nxz);         % round up length of zero-padded v to next power of 2
    nf=pow2(1,enx-(fnx==0.5));      % actual length of dft is a power of 2
    fmin=max(x1,1);                 % range of indices to extract from circular convolution
    fmax=min(x2,min(nf,nx+nh-1));
    zmin=max(1,2-x1);               % range of indices for non-zero output values
    zmax=zmin+fmax-fmin;
    fh=fft(h,nf,1);
else                % h is the z output of a previous call with the 'z' option
    d=h.d;          % x dimension to convolve along
    nx=h.nx;        % original size(x,d)
    ns=h.ns;        % number of dimensions of x
    nv=h.nv;        % number x values needed (might be < nx)
    vmin=h.vmin;    % first x value needed (might be > 1)
    vmax=h.vmax;    % last x value needed (might be < nx)
    nf=h.nf;        % dft length
    fh=h.fh;        % dft of impulse response h
    fmin=h.fmin;    % first fh value needed
    fmax=h.fmax;    % last fh value needed
    nz=h.nz;        % output size(z,d)
    zmin=h.zmin;    % first non-zero z value
    zmax=h.zmax;    % last non-zero z value
    s=size(x);      % get the structure of x
    k=numel(x)/nx;  % total size of the rest of x
    if size(x,d)~=nx || length(s)~=ns
        error('input x has incompatible dimensions');
    end
end
if  ~isstruct(h) && any(m=='z')    % save required information as a structure for next call
    z.d=d;          % x dimension to convolve along
    z.nx=nx;        % original size(x,d)
    z.ns=ns;        % number of dimensions of x
    z.nv=nv;        % number x values needed (might be < nx)
    z.vmin=vmin;    % first x value needed (might be > 1)
    z.vmax=vmax;    % last x value needed (might be < nx)
    z.nf=nf;        % dft length
    z.fh=fh;        % dft of impulse response h
    z.fmin=fmin;    % first fh value needed
    z.fmax=fmax;    % last fh value needed
    z.nz=nz;        % output size(z,d)
    z.zmin=zmin;    % first non-zero z value
    z.zmax=zmax;    % last non-zero z value
elseif ~isempty(x)  % there is data to convolve
    if d==1         % reshape x if necessary
        v=reshape(x,nx,k);
    else
        v=reshape(permute(x,[d:ns 1:d-1]),nx,k);
    end
    if nv<nx        % we can delete some unused v values
        v=v(vmin:vmax,:);
    end
    zz=ifft(fft(v,nf,1).*repmat(fh,1,k));    % do the convolution
    z=zeros(nz,k);                  % reserve space for the output
    z(zmin:zmax,:)=zz(fmin:fmax,:); % extract values from the convolution
    s(d)=nz;                        % update the size of the output
    if d==1
        z=reshape(z,s);
    else
        z=permute(reshape(z,s([d:ns 1:d-1])),[ns+2-d:ns 1:ns+1-d]);
    end
else
    z=[];           % no input data
end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [Cl, Cd] = a2clcd(alfa)

% [Cl, Cd] = a2clcd(alfa) computes hydrodynamic coefficients Cl and Cd
% for the fins

% Costants
CLa = 2.865;               % CLalfa [rad^-1]
CDmin = 0.0115;            % CDmin
K = 0.1309;                % K
ALFA1 = 0.419;             % alfa_stall [rad]
ALFA2 = 0.7854;            % alfa45 [rad]

C1 = -1.0572;              % interpolating coefficients
C2 = 1.6434;               % between zone 1 and 3
C3 = 1.6759;
C4 = -0.5021;

CT = 1.15;

% sign correction
mod_alfa=abs(alfa-sign(alfa)*(abs(alfa)>pi/2)*pi);

if mod_alfa < ALFA1
   % zone 1
   Cl = CLa * mod_alfa ;
   Cd = CDmin + K * Cl^2 ;
elseif  mod_alfa < ALFA2
   % zone2
   Cl = C1 * mod_alfa + C2;
   Cd = C3 * mod_alfa + C4;
else
  % zone 3 , piastra
   Cl = CT*cos(mod_alfa);
   Cd = CT*sin(mod_alfa);
end;

% sign correction
Cl = Cl*sign(sin(2*alfa));

author: Sable

repo: mcbench-benchmarks


text:
function hist = read_ctf_hist(filename)

% READ_CTF_HIST

% Copyright (C) 2010, Arjen Stolk

[p, f, x] = fileparts(filename);

if strcmp(x, '.ds')
  % the filename should point to the hist file, not the ds
  filename = fullfile(p, [f x], [f '.hist']);
end

fileline = 0;
fid      = fopen(filename,'r');

% the hist file is empty sometimes, in which case the actual reading does not return any useful information
hist = [];

while fileline >= 0
  fileline = fgets(fid);
  if ~isempty(findstr(fileline,'Collection started'))
    startdate = sscanf(fileline(findstr(fileline,'Collection started:'):end),'Collection started: %s');
    hist.starttime = sscanf(fileline(findstr(fileline,startdate):end),strcat(startdate, '%s'));
    hist.startdate = startdate;
  end
  if ~isempty(findstr(fileline,'Collection stopped'))
    stopdate = sscanf(fileline(findstr(fileline,'Collection stopped:'):end),'Collection stopped: %s');
    hist.stoptime = sscanf(fileline(findstr(fileline,stopdate):end),strcat(stopdate, '%s'));
    hist.stopdate = stopdate;
  end
  if ~isempty(findstr(fileline,'Dataset name'))
    hist.datasetname = sscanf(fileline(findstr(fileline,'Dataset name'):end),'Dataset name %s');
  end
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function w = dwt(x, J, af)

% Discrete 1-D Wavelet Transform
%
% USAGE:
%    w = dwt(x, J, af)
% INPUT:
%    x - N-point vector, where
%            1) N is divisible by 2^J
%            2) N >= 2^(J-1)*length(af)
%    J - number of stages
%    af - analysis filters
%    af(:, 1) - lowpass filter (even length)
%    af(:, 2) - highpass filter (evenlength)
% OUTPUT:
%    w{j}, j = 1..J+1 - DWT coefficients
% EXAMPLE:
%    [af, sf] = farras;
%    x = rand(1,64);
%    w = dwt(x,3,af);
%    y = idwt(w,3,sf);
%    err = x - y;
%    max(abs(err))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

for k = 1:J
    [x w{k}] = afb(x, af);
end
w{J+1} = x;

author: thomaskuestner


text:
%==================================================================
%                        General Data File
% Title: TETREAHEDRA
% Units: SI
% Dimensions: 3D
% Type of problem: Plane_Stress
% Type of Phisics: ELASTIC
% Micro/Macro: MACRO
%
%==================================================================

%% Data

Data_prb = {
'TETRAHEDRA';
'SI';
'3D';
'Plane_Stress';
'ELASTIC';
'MACRO';
};

%% Coordinates
% Node                X                Y                Z

coord = [
1            2            0            1
2            2            0          0.5
3            2          0.5            1
4            2          0.5          0.5
5            2            0            0
6            2            1            1
7            1            0            1
8            1            0          0.5
9            1          0.5            1
10            2            1          0.5
11            2          0.5            0
12            1          0.5          0.5
13            2            1            0
14            1            0            0
15            1            1            1
16            1          0.5            0
17            1            1          0.5
18            1            1            0
19            0            0            1
20            0            0          0.5
21            0          0.5            1
22            0          0.5          0.5
23            0            1            1
24            0            0            0
25            0          0.5            0
26            0            1          0.5
27            0            1            0
];

%% Conectivities
% Element        Node(1)                Node(2)                Node(3)                Node(4)                Material

connec = [
1 10 17 12 6 0
2 6 10 4 12 0
3 15 9 12 6 0
4 6 3 9 12 0
5 6 15 17 12 0
6 3 4 12 6 0
7 11 16 12 13 0
8 13 11 4 12 0
9 18 17 12 13 0
10 13 10 17 12 0
11 13 18 16 12 0
12 10 4 12 13 0
13 9 21 23 12 0
14 12 9 15 23 0
15 22 26 23 12 0
16 12 17 26 23 0
17 12 22 21 23 0
18 17 15 23 12 0
19 17 26 27 12 0
20 12 17 18 27 0
21 22 25 27 12 0
22 12 16 25 27 0
23 12 22 26 27 0
24 16 18 27 12 0
25 3 9 12 1 0
26 1 3 4 12 0
27 7 8 12 1 0
28 1 2 8 12 0
29 1 7 9 12 0
30 2 4 12 1 0
31 2 8 12 5 0
32 5 2 4 12 0
33 14 16 12 5 0
34 5 11 16 12 0
35 5 14 8 12 0
36 11 4 12 5 0
37 8 20 19 12 0
38 12 8 7 19 0
39 22 21 19 12 0
40 12 9 21 19 0
41 12 22 20 19 0
42 9 7 19 12 0
43 16 25 24 12 0
44 12 16 14 24 0
45 22 20 24 12 0
46 12 8 20 24 0
47 12 22 25 24 0
48 8 14 24 12 0
];

%% Variable Prescribed
% Node            Dimension                Value

dirichlet_data = [
19 1 0
19 2 0
19 3 0
20 1 0
20 2 0
20 3 0
21 1 0
21 2 0
21 3 0
22 1 0
22 2 0
22 3 0
23 1 0
23 2 0
23 3 0
24 1 0
24 2 0
24 3 0
25 1 0
25 2 0
25 3 0
26 1 0
26 2 0
26 3 0
27 1 0
27 2 0
27 3 0
];

%% Force Prescribed
% Node                Dimension                Value

pointload_complete = [
4 1 0
4 2 -1
];

%% Volumetric Force
% Element        Dim                Force_Dim

Vol_force = [
];

%% Group Elements
% Element        Group_num

Group = [
];

%% Initial Holes
% Elements that are considered holes initially
% Element

Initial_holes = [
];

%% Boundary Elements
% Elements that can not be removed
% Element

Boundary_elements = [
];

%% Micro gauss post
%
% Element

Micro_gauss_post = [
];

%% Micro Slave-Master
% Nodes that are Slaves
% Nodes             Value (1-Slave,0-Master)

Micro_slave = [
];

%% Nodes solid
% Nodes that must remain
% Nodes

nodesolid = unique(pointload_complete(:,1));

%% External border Elements
% Detect the elements that define the edge of the domain
% Element               Node(1)           Node(2)

External_border_elements = [
];

%% External border Nodes
% Detect the nodes that define the edge of the domain
% Node

External_border_nodes = [
];

%% Materials
% Materials that have been used
% Material_Num              Mat_density        Young_Modulus        Poisson

Materials = [
];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function net = init_resnet_toyDigitV2(varargin)
%CNN_IMAGENET_INIT_RESNET  Initialize the ResNet-50 model for ImageNet classification

opts.classNames = {} ;
opts.classDescriptions = {} ;
opts.averageImage = zeros(3,1) ;
opts.colorDeviation = zeros(3) ;
opts.cudnnWorkspaceLimit = 1024*1024*1204 ; % 1GB
opts = vl_argparse(opts, varargin) ;

net = dagnn.DagNN() ;

lastAdded.var = 'input' ;
lastAdded.depth = 3 ;

    function Conv(name, ksize, depth, varargin)
        % Helper function to add a Convolutional + BatchNorm + ReLU
        % sequence to the network.
        args.relu = true ;
        args.downsample = false ;
        args.bias = false ;
        args = vl_argparse(args, varargin) ;
        if args.downsample, stride = 2 ; else stride = 1 ; end
        if args.bias, pars = {[name '_f'], [name '_b']} ; else pars = {[name '_f']} ; end
        net.addLayer([name  '_conv'], ...
            dagnn.Conv('size', [ksize ksize lastAdded.depth depth], ...
            'stride', stride, ....
            'pad', (ksize - 1) / 2, ...
            'hasBias', args.bias, ...
            'opts', {'cudnnworkspacelimit', opts.cudnnWorkspaceLimit}), ...
            lastAdded.var, ...
            [name '_conv'], ...
            pars) ;
        net.addLayer([name '_bn'], ...
            dagnn.BatchNorm('numChannels', depth, 'epsilon', 1e-5), ...
            [name '_conv'], ...
            [name '_bn'], ...
            {[name '_bn_w'], [name '_bn_b'], [name '_bn_m']}) ;
        lastAdded.depth = depth ;
        lastAdded.var = [name '_bn'] ;
        if args.relu
            net.addLayer([name '_relu'] , ...
                dagnn.ReLU(), ...
                lastAdded.var, ...
                [name '_relu']) ;
            lastAdded.var = [name '_relu'] ;
        end
    end

% -------------------------------------------------------------------------
% Add input section
% -------------------------------------------------------------------------

Conv('conv1', 3, 32, ...
    'relu', true, ...
    'bias', false, ...
    'downsample', 0) ;

% net.addLayer(...
%   'conv1_pool' , ...
%   dagnn.Pooling('poolSize', [3 3], ...
%                 'stride', 2, ...
%                 'pad', 1,  ...
%                 'method', 'max'), ...
%   lastAdded.var, ...
%   'conv1') ;
lastAdded.var = 'conv1_relu' ;

% -------------------------------------------------------------------------
% Add intermediate sections
% -------------------------------------------------------------------------

for s = 2:5

    switch s
        case 2, sectionLen = 3 ;
        case 3, sectionLen = 3 ; % 4 8 ;
        case 4, sectionLen = 3 ; % 6 23 ; % 36 ;
        case 5, sectionLen = 3 ;
    end

    % -----------------------------------------------------------------------
    % Add intermediate segments for each section
    for l = 1:sectionLen
        depth = 2^(s+2) ;
        sectionInput = lastAdded ;
        name = sprintf('conv%d_%d', s, l)  ;

        % Optional adapter layer
        if l == 1
%             Conv([name '_adapt_conv'], 1, 2^(s+6), 'downsample', s >= 3, 'relu', false) ;
            Conv([name '_adapt_conv'], 1, 2^(s+4), 'downsample', 0, 'relu', false) ;
        end
        sumInput = lastAdded ;

        % ABC: 1x1, 3x3, 1x1; downsample if first segment in section from
        % section 2 onwards.
        lastAdded = sectionInput ;
        %Conv([name 'a'], 1, 2^(s+4), 'downsample', (s >= 3) & l == 1) ;
        %Conv([name 'b'], 3, 2^(s+4)) ;
        Conv([name 'a'], 1, 2^(s+2)) ;
%         Conv([name 'b'], 3, 2^(s+4), 'downsample', (s >= 3) & l == 1) ;
        Conv([name 'b'], 3, 2^(s+2), 'downsample', 0 ) ;
        Conv([name 'c'], 1, 2^(s+4), 'relu', false) ;

        % Sum layer
        net.addLayer([name '_sum'] , ...
            dagnn.Sum(), ...
            {sumInput.var, lastAdded.var}, ...
            [name '_sum']) ;
        net.addLayer([name '_relu'] , ...
            dagnn.ReLU(), ...
            [name '_sum'], ...
            name) ;
        lastAdded.var = name ;
    end
end

% net.addLayer('prediction_avg' , ...
%     dagnn.Pooling('poolSize', [7 7], 'method', 'avg'), ...
%     lastAdded.var, ...
%     'prediction_avg') ;
%
% net.addLayer('prediction' , ...
%     dagnn.Conv('size', [1 1 2048 1000]), ...
%     'prediction_avg', ...
%     'prediction', ...
%     {'prediction_f', 'prediction_b'}) ;
%
% net.addLayer('loss', ...
%     dagnn.Loss('loss', 'softmaxlog') ,...
%     {'prediction', 'label'}, ...
%     'objective') ;
%
% net.addLayer('top1error', ...
%     dagnn.Loss('loss', 'classerror'), ...
%     {'prediction', 'label'}, ...
%     'top1error') ;
%
% net.addLayer('top5error', ...
%     dagnn.Loss('loss', 'topkerror', 'opts', {'topK', 5}), ...
%     {'prediction', 'label'}, ...
%     'top5error') ;

% -------------------------------------------------------------------------
%                                                           Meta parameters
% -------------------------------------------------------------------------

% net.meta.normalization.imageSize = [224 224 3] ;
% net.meta.inputSize = [net.meta.normalization.imageSize, 32] ;
% net.meta.normalization.cropSize = net.meta.normalization.imageSize(1) / 256 ;
% net.meta.normalization.averageImage = opts.averageImage ;
%
% net.meta.classes.name = opts.classNames ;
% net.meta.classes.description = opts.classDescriptions ;
%
% net.meta.augmentation.jitterLocation = true ;
% net.meta.augmentation.jitterFlip = true ;
% net.meta.augmentation.jitterBrightness = double(0.1 * opts.colorDeviation) ;
% net.meta.augmentation.jitterAspect = [3/4, 4/3] ;
% net.meta.augmentation.jitterScale  = [0.4, 1.1] ;
% %net.meta.augmentation.jitterSaturation = 0.4 ;
% %net.meta.augmentation.jitterContrast = 0.4 ;
%
% net.meta.inputSize = {'input', [net.meta.normalization.imageSize 32]} ;
%
% %lr = logspace(-1, -3, 60) ;
% lr = [0.1 * ones(1,30), 0.01*ones(1,30), 0.001*ones(1,30)] ;
% net.meta.trainOpts.learningRate = lr ;
% net.meta.trainOpts.numEpochs = numel(lr) ;
% net.meta.trainOpts.momentum = 0.9 ;
% net.meta.trainOpts.batchSize = 256 ;
% net.meta.trainOpts.numSubBatches = 4 ;
% net.meta.trainOpts.weightDecay = 0.0001 ;

% Init parameters randomly
net.initParams() ;

% For uniformity with the other ImageNet networks, t
% the input data is *not* normalized to have unit standard deviation,
% whereas this is enforced by batch normalization deeper down.
% The ImageNet standard deviation (for each of R, G, and B) is about 60, so
% we adjust the weights and learing rate accordingly in the first layer.
%
% This simple change improves performance almost +1% top 1 error.
% p = net.getParamIndex('conv1_f') ;
% net.params(p).value = net.params(p).value / 100 ;
% net.params(p).learningRate = net.params(p).learningRate / 100^2 ;
%
% for l = 1:numel(net.layers)
%     if isa(net.layers(l).block, 'dagnn.BatchNorm')
%         k = net.getParamIndex(net.layers(l).params{3}) ;
%         net.params(k).learningRate = 0.3 ;
%     end
% end

end

author: aimerykong


text:
function x = normal_ms_cdf_inv ( cdf, mu, sigma )

%*****************************************************************************80
%
%% NORMAL_MS_CDF_INV inverts the Normal CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Input, real MU, SIGMA, the parameters of the PDF.
%    0.0 < SIGMA.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NORMAL_MS_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'NORMAL_MS_CDF_INV - Fatal error!' );
  end

  x2 = normal_01_cdf_inv ( cdf );

  x = mu + sigma * x2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function smallpot = marginalize_pot(bigpot, keepdom, sumoverdom, nodesizes)
% MARGINALIZE_POT Marginalize a mpot onto a smaller domain.
% smallpot = marginalize_pot(bigpot, keep)

keepsize = sum(nodesizes(keepdom));
[A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(bigpot.A, bigpot.B, bigpot.C, keepdom, sumoverdom, nodesizes);
smallpot = scgcpot(keepsize, bigpot.ctailsize, bigpot.p, A1, B1, C11);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clear all; close all; clc;

load catData_w.mat
load dogData_w.mat
CD=[dog_wave cat_wave];
[u,s,v]=svd(CD-mean(CD(:)));

xtrain=[v(1:60,2:2:4); v(81:140,2:2:4)];
label=[ones(60,1); -1*ones(60,1)];
test=[v(61:80,2:2:4); v(141:160,2:2:4)];

class=classify(test,xtrain,label);
subplot(4,1,3), bar(class,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis off

truth=[ones(20,1); -1*ones(20,1)];
E=100-sum(0.5*abs(class-truth))/40*100

plot(v(1:80,2),v(1:80,4),'ro','Linewidth',[1],'MarkerEdgeColor','k','MarkerFaceColor',[0 1 0.2],'MarkerSize',8), hold on
plot(v(81:end,2),v(81:end,4),'bo','Linewidth',[1],'MarkerEdgeColor','k','MarkerFaceColor',[0.9 0 1],'MarkerSize',8), set(gca,'Fontsize',[15])

figure(2)
for j=1:2
    subplot(2,2,j)
    U3=flipud(reshape(u(:,2*j),32,32));
    pcolor(U3), colormap(hot), axis off
end

%%
load dogData.mat
load catData.mat
CD=double([dog cat]);
[u,s,v]=svd(CD-mean(CD(:)));

xtrain=[v(1:60,2:2:4); v(81:140,2:2:4)];
label=[ones(60,1); -1*ones(60,1)];
test=[v(61:80,2:2:4); v(141:160,2:2:4)];

class=classify(test,xtrain,label);
subplot(4,1,4), bar(class,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis off

%% cross-validate

for jj=1:100;
    r1=randperm(80); r2=randperm(80);
    ind1=r1(1:60); ind2=r2(1:60)+60;
    ind1t=r1(61:80); ind2t=r2(61:80)+60;

    xtrain=[v(ind1,2:2:4); v(ind2,2:2:4)];
    test=[v(ind1t,2:2:4); v(ind2t,2:2:4)];

    label=[ones(60,1); -1*ones(60,1)];
    truth=[ones(20,1); -1*ones(20,1)];
    class=classify(test,xtrain,label);
    E(jj)=sum(abs(class-truth))/40*100;
end

figure(4)
bar(E,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis([0 100 0 100])
hold on
plot([0 100],[mean(E) mean(E)],'r:','Linewidth',[2])
set(gca,'Fontsize',[15])

%%
close all
figure(1)
load catData_w.mat
load dogData_w.mat
CD=[dog_wave cat_wave];
[u,s,v]=svd(CD-mean(CD(:)));

for j=1:2
    subplot(2,2,j)
    plot(v(1:80,2),v(1:80,4),'ro','Linewidth',[1],'MarkerEdgeColor','k',...
        'MarkerFaceColor',[0 1 0.2],...
        'MarkerSize',8), hold on   % [0.49 1 .63]
    plot(v(81:end,2),v(81:end,4),'bo','Linewidth',[1],'MarkerEdgeColor','k',...
        'MarkerFaceColor',[0.9 0 1],...
        'MarkerSize',8)
    set(gca,'Fontsize',[15]), hold on
end

xtrain=[v(1:60,2:2:4); v(81:140,2:2:4)];
label=[ones(60,1); -1*ones(60,1)];
test=[v(61:80,2:2:4); v(141:160,2:2:4)];

subplot(2,2,1)
[class,~,~,~,coeff]=classify(test,xtrain,label);
K = coeff(1,2).const;
L = coeff(1,2).linear;
f = @(x,y) K + [x y]*L;
h2 = ezplot(f,[-.15 0.25 -.3 0.2]);
set(h2,'Linecolor','k','Linewidth',[2])
xlabel(''),ylabel(''),title('')

subplot(2,2,2)
[class,~,~,~,coeff]=classify(test,xtrain,label,'quadratic');
K = coeff(1,2).const;
L = coeff(1,2).linear;
Q = coeff(1,2).quadratic;
f = @(x,y) K + [x y]*L + sum(([x y]*Q) .* [x y], 2);
h2 = ezplot(f,[-.15 0.25 -.3 0.2]);
set(h2,'Linecolor','k','Linewidth',[2])
xlabel(''),ylabel(''),title('')

%set(gca,'Linewidth',[2])
%set(gcf,'LineStyle','k-')

%% figure 0
figure(5), subplot(2,2,1)

load catData_w.mat
load dogData_w.mat
CD=[dog_wave cat_wave];
[u,s,v]=svd(CD-mean(CD(:)));

plot(v(1:80,2),v(1:80,4),'ro','Linewidth',[1],'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0 1 0.2],...
    'MarkerSize',8), hold on   % [0.49 1 .63]
plot(v(81:end,2),v(81:end,4),'bo','Linewidth',[1],'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.9 0 1],...
    'MarkerSize',8)
axis off

x=linspace(-0.1,-.4,200); %x2=cos(theta)*x + sin(theta)*y;
y=0.2*exp(-1000*(x2+0.25).^2);
theta=-20; Ar=[cos(theta) sin(theta); -sin(theta) cos(theta)];
%x2=cos(theta)*x + sin(theta)*y;
%y2=-sin(theta)*x + cos(theta)*y;
hold on
fill(x,y,[240 194 224]/255,'MarkerEdgeColor',[0.9 0 1],'Linewidth',[2])
set(gca,'Fontsize',[15])

plot([-0.05  0.25],[-.5 -0.25],'k')

clear x
clear y
figure(6)
x=-7:0.01:7;
y=exp(-(x-2).^2);
yn=y/trapz(x,y);
y2=exp(-0.5*(x+2).^2);
yn2=y2/trapz(x,y2);
fill(x,yn,[240 194 224]/255,'Linewidth',[2]), axis off
hold on
fill(x,yn2,[179 255 179]/255,'Linewidth',[2]),
alpha(0.5)
set(gca,'Linewidth',[2])

figure(7)
y=exp(-0.3*(x-0.6).^2);
yn=y/trapz(x,y);
y2=exp(-0.2*(x+0.5).^2);
yn2=y2/trapz(x,y2);
fill(x,yn,[240 194 224]/255,'Linewidth',[2]), axis off
hold on
fill(x,yn2,[179 255 179]/255,'Linewidth',[2]),
alpha(0.5)
set(gca,'Linewidth',[2])

author: dynamicslab


text:
function H = entropy(v, scale)
% ENTROPY Entropy log base 2
% H = entropy(v)
% If v is a matrix, we compute the entropy of each column
%
% % H = entropy(v,1) means we scale the result so that it lies in [0,1]

if nargin < 2, scale = 0; end

v = v + (v==0);
H = -1 * sum(v .* log2(v), 1); % sum the rows

if scale
  n = size(v, 1);
  unif = normalise(ones(n,1));
  H = H / entropy(unif);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% The COBRAToolbox: testConvertCNAModelToCbModel.m
%
% Purpose:
%     - convert a CNA model to a COBRA model and vice versa
%

% save the current path
currentDir = pwd;

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

% check if CNA is properly installed
statusCNA = checkCNAinstallation();

if statusCNA
    model = getDistributedModel('ecoli_core_model.mat');

    % load reference data
    load('refData_cnaModel.mat')

    % convert a COBRA model to a CNA model
    cnaModel_new = convertCbModelToCNAModel(model);

    % assert if the cnaModel_new is equal to the reference CNA model
    if 0
        assert(isequaln(cnaModel_new, cnaModel));
    else
        % function [cnap, errval]= CNAgenerateMFNetwork(cnap,nodisp) in CellNetAnalyzer: API function CNAgenerateMFNetwork
        % has been updated to include new default fields, that do not
        % match, so remove them. Probably cnaModel_new needs to be updated.
        if isfield(cnaModel_new,'color1')
            cnaModel_new= rmfield(cnaModel_new,{'color1','color2','color3','color4','textColor','macroSynthColor','specBoxColor','unsaved_changes'});
        end
        assert(isequaln(cnaModel_new, cnaModel));
    end

    % convert the reference CNA model to the COBRA model
    model_new = convertCNAModelToCbModel(cnaModel);

    assert(isequal(model_new.rxns, model.rxns));
    assert(isequal(model_new.mets, model.mets));
    assert(isequal(model_new.S, model.S));
    assert(isequal(model_new.lb, model.lb));
    assert(isequal(model_new.ub, model.ub));
    %Osense is no longer a model field, so we need to check, whether the
    %generated LPproblems are the same (or we have to make a lot of if
    %statements...)
    LPproblemA = buildLPproblemFromModel(model_new);
    LPproblemB = buildLPproblemFromModel(model);
    assert(isequal(LPproblemA.osense*LPproblemA.c, LPproblemB.osense*LPproblemB.c));
    assert(isequal(model_new.metNames, model.metNames));
    assert(isequal(model_new.b, model.b));
    % Note: rxnNames is different
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function rd = gsp_resistance_distance(G,param)
%GSP_RESISTANCE_DISTANCE see gsp_resistance_distances

if nargin<2
    param = struct;
end

rd = gsp_resistance_distances(G,param);

end

repo: gspbox


text:
function [pos] = moveinward(pos, move)
%This functions moves 'pos' inward according to their normals by 'move'
%units
propos = elproj(pos); % projection to 2D
tri = delaunay(propos); %creates delaunay triangulation of 2D plane, which will be used for the the 3D case
nor = normals(pos,tri); %compute normals of surface
ori = surfaceorientation(pos, tri, nor);
if ori==1
  % the normals are outward oriented
elseif ori==-1
  % the normals are inward oriented
  nor = -nor;
else
  ft_warning('cannot determine the orientation of the vertex normals');
end
pos = pos-move*nor; % moves pos inwards according to their normals
end

author: spm


text:
% Display the contents of a pyramid, as returned by functions
% 'laplacian_pyramid' or 'gaussian pyramid'
%
% tom.mertens@gmail.com, August 2007
%

function display_pyramid(pyr)

L = length(pyr);
r = size(pyr{1},1);
c = size(pyr{1},2);
k = size(pyr{1},3);
R = zeros(r,2*c,k);

offset = 1;
for l = 1:L
    I = pyr{l};
    r = size(I,1);
    c = size(I,2);
    R(1:r, offset:offset-1+c, :) = I;
    offset = offset + c;
end

if (min(R(:)) < 1e-5)
    %make negative values displayable
    a = min(R(:));
    b = max(R(:));
    R = (R - a) / (b - a);
end

figure; imshow(R);

author: mahmoudnafifi


text:
function [x, f, eflag, output, lambda] = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
%MIQPS_GLPK  Mixed Integer Linear Program Solver based on GLPK - GNU Linear Programming Kit.
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
%       MIQPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GLPK(PROBLEM)
%   A wrapper function providing a standardized interface for using
%   GLKP to solve the following LP (linear programming) problem:
%
%       min C'*X
%        X
%
%   subject to
%
%       L <= A*X <= U       (linear constraints)
%       XMIN <= X <= XMAX   (variable bounds)
%
%   Inputs (all optional except H, C, A and L):
%       H : dummy matrix (possibly sparse) of quadratic cost coefficients
%           for QP problems, which GLPK does not handle
%       C : vector of linear cost coefficients
%       A, L, U : define the optional linear constraints. Default
%           values for the elements of L and U are -Inf and Inf,
%           respectively.
%       XMIN, XMAX : optional lower and upper bounds on the
%           X variables, defaults are -Inf and Inf, respectively.
%       X0 : optional starting value of optimization vector X (NOT USED)
%       VTYPE : character string of length NX (number of elements in X),
%               or 1 (value applies to all variables in x),
%               allowed values are 'C' (continuous), 'B' (binary) or
%               'I' (integer).
%       OPT : optional options structure with the following fields,
%           all of which are also optional (default values shown in
%           parentheses)
%           verbose (0) - controls level of progress output displayed
%               0 = no progress output
%               1 = some progress output
%               2 = verbose progress output
%           skip_prices (0) - flag that specifies whether or not to
%               skip the price computation stage, in which the problem
%               is re-solved for only the continuous variables, with all
%               others being constrained to their solved values
%           price_stage_warn_tol (1e-7) - tolerance on the objective fcn
%               value and primal variable relative match required to avoid
%               mis-match warning message
%           glpk_opt - options struct for GLPK, value in verbose
%                   overrides these options
%       PROBLEM : The inputs can alternatively be supplied in a single
%           PROBLEM struct with fields corresponding to the input arguments
%           described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt
%
%   Outputs:
%       X : solution vector
%       F : final objective function value
%       EXITFLAG : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other
%       OUTPUT : struct with errnum and status fields from GLPK output args
%       LAMBDA : struct containing the Langrange and Kuhn-Tucker
%           multipliers on the constraints, with fields:
%           mu_l - lower (left-hand) limit on linear constraints
%           mu_u - upper (right-hand) limit on linear constraints
%           lower - lower bound on optimization variables
%           upper - upper bound on optimization variables
%
%   Note the calling syntax is almost identical to that of GLPK. The main
%   difference is that the linear constraints are specified with A, L, U
%   instead of A, B, Aeq, Beq.
%
%   Calling syntax options:
%       [x, f, exitflag, output, lambda] = ...
%           miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
%
%       x = miqps_glpk(H, c, A, l, u)
%       x = miqps_glpk(H, c, A, l, u, xmin, xmax)
%       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0)
%       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype)
%       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
%       x = miqps_glpk(problem), where problem is a struct with fields:
%                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
%                       all fields except 'c', 'A' and 'l' or 'u' are optional
%       x = miqps_glpk(...)
%       [x, f] = miqps_glpk(...)
%       [x, f, exitflag] = miqps_glpk(...)
%       [x, f, exitflag, output] = miqps_glpk(...)
%       [x, f, exitflag, output, lambda] = miqps_glpk(...)
%
%
%   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
%       H = [   1003.1  4.3     6.3     5.9;
%               4.3     2.2     2.1     3.9;
%               6.3     2.1     3.5     4.8;
%               5.9     3.9     4.8     10  ];
%       c = zeros(4,1);
%       A = [   1       1       1       1;
%               0.17    0.11    0.10    0.18    ];
%       l = [1; 0.10];
%       u = [1; Inf];
%       xmin = zeros(4,1);
%       x0 = [1; 0; 0; 1];
%       opt = struct('verbose', 2);
%       [x, f, s, out, lambda] = miqps_glpk(H, c, A, l, u, xmin, [], x0, vtype, opt);
%
%   See also MIQPS_MASTER, GLPK.

%   MP-Opt-Model
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% check for Optimization Toolbox
% if ~have_feature('quadprog')
%     error('miqps_glpk: requires the MEX interface to GLPK');
% end

%%----- input argument handling  -----
%% gather inputs
if nargin == 1 && isstruct(H)       %% problem struct
    p = H;
    if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
    if isfield(p, 'vtype'), vtype = p.vtype;else,   vtype = []; end
    if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
    if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
    if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
    if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
    if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
    if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
    if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
    if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
else                                %% individual args
    if nargin < 10
        opt = [];
        if nargin < 9
            vtype = [];
            if nargin < 8
                x0 = [];
                if nargin < 7
                    xmax = [];
                    if nargin < 6
                        xmin = [];
                    end
                end
            end
        end
    end
end

%% define nx, set default values for missing optional inputs
if isempty(H) || ~any(any(H))
    if isempty(A) && isempty(xmin) && isempty(xmax)
        error('miqps_glpk: LP problem must include constraints or variable bounds');
    else
        if ~isempty(A)
            nx = size(A, 2);
        elseif ~isempty(xmin)
            nx = length(xmin);
        else    % if ~isempty(xmax)
            nx = length(xmax);
        end
    end
else
    error('miqps_glpk: GLPK handles only LP problems, not QP problems');
    nx = size(H, 1);
end
if isempty(c)
    c = zeros(nx, 1);
end
if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
                                 (isempty(u) || all(u == Inf)))
    A = sparse(0,nx);           %% no limits => no linear constraints
end
nA = size(A, 1);                %% number of original linear constraints
if isempty(u)                   %% By default, linear inequalities are ...
    u = Inf(nA, 1);             %% ... unbounded above and ...
end
if isempty(l)
    l = -Inf(nA, 1);            %% ... unbounded below.
end
if isempty(xmin)                %% By default, optimization variables are ...
    xmin = -Inf(nx, 1);         %% ... unbounded below and ...
end
if isempty(xmax)
    xmax = Inf(nx, 1);          %% ... unbounded above.
end
if isempty(x0)
    x0 = zeros(nx, 1);
end

%% default options
if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
    verbose = opt.verbose;
else
    verbose = 0;
end

%% split up linear constraints
ieq = find( abs(u-l) <= eps );          %% equality
igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
bb = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)];

%% grab some dimensions
nlt = length(ilt);      %% number of upper bounded linear inequalities
ngt = length(igt);      %% number of lower bounded linear inequalities
nbx = length(ibx);      %% number of doubly bounded linear inequalities
neq = length(ieq);      %% number of equalities
nie = nlt+ngt+2*nbx;    %% number of inequalities

ctype = [repmat('S', neq, 1); repmat('U', nlt+ngt+2*nbx, 1)];

if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I'))
    mi = 0;
    vtype = repmat('C', nx, 1);
else
    mi = 1;
    %% expand vtype to nx elements if necessary
    if length(vtype) == 1 && nx > 1
        vtype = char(vtype * ones(nx, 1));
    elseif size(vtype, 2) > 1   %% make sure it's a col vector
        vtype = vtype';
    end
end
%% convert 'B' variables to 'I' and clip bounds to [0, 1]
k = find(vtype == 'B');
if ~isempty(k)
    kk = find(xmax(k) > 1);
    xmax(k(kk)) = 1;
    kk = find(xmin(k) < 0);
    xmin(k(kk)) = 0;
    vtype(k) = 'I';
end

%% set options struct for GLPK
if ~isempty(opt) && isfield(opt, 'glpk_opt') && ~isempty(opt.glpk_opt)
    glpk_opt = glpk_options(opt.glpk_opt);
else
    glpk_opt = glpk_options;
end
glpk_opt.msglev = verbose;

%% call the solver
[x, f, errnum, extra] = ...
    glpk(c, AA, bb, xmin, xmax, ctype, vtype, 1, glpk_opt);

%% set exit flag
if isfield(extra, 'status')             %% status found in extra.status
    output.errnum = errnum;
    output.status = extra.status;
    eflag = -errnum;
    if eflag == 0 && extra.status == 5
        eflag = 1;
    end
else                                    %% status found in errnum
    output.errnum = [];
    output.status = errnum;
    if have_feature('octave')
        if errnum == 180 || errnum == 151 || errnum == 171
            eflag = 1;
        else
            eflag = 0;
        end
    else
        if errnum == 5
            eflag = 1;
        else
            eflag = 0;
        end
    end
end

%% repackage lambdas
if isempty(extra) || ~isfield(extra, 'lambda') || isempty(extra.lambda)
    lambda = struct( ...
        'mu_l', zeros(nA, 1), ...
        'mu_u', zeros(nA, 1), ...
        'lower', zeros(nx, 1), ...
        'upper', zeros(nx, 1) ...
    );
else
    lam.eqlin = extra.lambda(1:neq);
    lam.ineqlin = extra.lambda(neq+(1:nie));
    lam.lower = extra.redcosts;
    lam.upper = -extra.redcosts;
    lam.lower(lam.lower < 0) = 0;
    lam.upper(lam.upper < 0) = 0;
    kl = find(lam.eqlin > 0);   %% lower bound binding
    ku = find(lam.eqlin < 0);   %% upper bound binding

    mu_l = zeros(nA, 1);
    mu_l(ieq(kl)) = lam.eqlin(kl);
    mu_l(igt) = -lam.ineqlin(nlt+(1:ngt));
    mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx));

    mu_u = zeros(nA, 1);
    mu_u(ieq(ku)) = -lam.eqlin(ku);
    mu_u(ilt) = -lam.ineqlin(1:nlt);
    mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));

    lambda = struct( ...
        'mu_l', mu_l, ...
        'mu_u', mu_u, ...
        'lower', lam.lower(1:nx), ...
        'upper', lam.upper(1:nx) ...
    );
end

if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
    if verbose
        fprintf('--- Integer stage complete, starting price computation stage ---\n');
    end
    if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol)
        tol = opt.price_stage_warn_tol;
    else
        tol = 1e-7;
    end
    k = find(vtype == 'I' | vtype == 'B');
    x(k) = round(x(k));
    xmin(k) = x(k);
    xmax(k) = x(k);
    x0 = x;
    opt.glpk_opt.lpsolver = 1;      %% simplex
    opt.glpk_opt.dual = 0;          %% primal simplex
    if have_feature('octave') && have_feature('octave', 'vnum') >= 3.007
        opt.glpk_opt.dual = 1;      %% primal simplex
    end

    [x_, f_, eflag_, output_, lambda] = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt);
    if eflag ~= eflag_
        error('miqps_glpk: EXITFLAG from price computation stage = %d', eflag_);
    end
    if abs(f - f_)/max(abs(f), 1) > tol
        warning('miqps_glpk: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
    end
    xn = x;
    xn(abs(xn)<1) = 1;
    [mx, k] = max(abs(x - x_) ./ xn);
    if mx > tol
        warning('miqps_glpk: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
    end
    output.price_stage = output_;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
    %% Regional Channel Selection
    % 5 Region
    channel_f=epoch.x(:,1:7,:);
    channel_o=epoch.x(:,28:32,:);
    channel_p=epoch.x(:,23:27,:);
    channel_c=epoch.x(:,[8:11,13:15,18:21],:);
    channel_t=epoch.x(:,[12,16:17,22],:);

    % 5 Region + Vertical
    channel_f1=epoch.x(:,[1,3:4],:);
    channel_f2=epoch.x(:,[2,6:7],:);
    channel_f3=epoch.x(:,5,:);
    channel_o1=epoch.x(:,29,:);
    channel_o2=epoch.x(:,31,:);
    channel_o3=epoch.x(:,30,:);
    channel_p1=epoch.x(:,[18:19,23:24,28],:);
    channel_p2=epoch.x(:,[20:21,26:27,32],:);
    channel_p3=epoch.x(:,25,:);
    channel_c1=epoch.x(:,[8:9,13],:);
    channel_c2=epoch.x(:,[10:11,15],:);
    channel_c3=epoch.x(:,14,:);
    channel_t1=epoch.x(:,[12,17],:);
    channel_t2=epoch.x(:,[16,22],:);


text:
function boundary_segment_num = p07_boundary_segment_num ( )

%*****************************************************************************80
%
%% P07_BOUNDARY_SEGMENT_NUM counts the boundary segments in problem 07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
%
  boundary_segment_num = 1;

  return
end


text:
function r8gb_print ( m, n, ml, mu, a, title )

%*****************************************************************************80
%
%% R8GB_PRINT prints a banded matrix.
%
%  Discussion:
%
%    An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU
%    is assumed to be entirely zero, except for the main diagonal, and
%    entries in the ML nearest subdiagonals, and MU nearest superdiagonals.
%
%    LINPACK and LAPACK "R8GB" storage for such a matrix generally includes
%    room for ML extra superdiagonals, which may be required to store
%    nonzero entries generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than min(M,N)-1..
%
%    Input, real A(2*ML+MU+1,N), the M by N band matrix, stored in LINPACK
%    or LAPACK general band storage mode.
%
%    Input, string TITLE, a title to be printed.
%
  r8gb_print_some ( m, n, ml, mu, a, 1, 1, m, n, title );

  return
end


text:
function assign (varargin)
% assign - make a list of assignments (matlab 5 or higher)
%
%	ASSIGN('VAR1', VAL1, 'VAR2', VAL2, ...) makes the assignments
%	VAR1 = VAL1; VAR2 = VAL2; ... in the caller's workspace.
%
%	This is most useful when passing in an option list to a
%	function.  Thus in the function which starts:
%		function foo(x,y,varargin)
%		z = 0;
%		assign(varargin{:});
%	the variable z can be given a non-default value by calling the
%	function like so: foo(x,y,'z',4);

vars = {varargin{1:2:end}};
vals = {varargin{2:2:end}};

% use length(vals) not length(vars) in case of odd number of arguments
for i = 1:length(vals), % changed 1 to a 2
  assignin('caller', vars{i}, vals{i});
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function y = base64encode(x, alg, isChunked, url_safe)
%BASE64ENCODE Perform base64 encoding on a string.
% INPUT:
%   x    - block of data to be encoded.  Can be a string or a numeric
%          vector containing integers in the range 0-255.
%   alg  - Algorithm to use: can take values 'java' or 'matlab'. Optional
%          variable defaulting to 'java' which is a little faster. If
%          'java' is chosen than core of the code is performed by a call to
%          a java library. Optionally all operations can be performed using
%          matleb code.
%   isChunked - encode output into 76 character blocks. The returned
%          encoded string is broken into lines of no more than
%          76 characters each, and each line will end with EOL. Notice that
%          if resulting string is saved as part of an xml file, those EOL's
%          are often stripped by xmlwrite funtrion prior to saving.
%   url_safe - use Modified Base64 for URL applications ('base64url'
%   encoding) "Base64 alphabet" ([A-Za-z0-9-_=]).
%
%
% OUTPUT:
%   y    - character array using only "Base64 alphabet" characters
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).
%   The Base64 encoding is designed to represent arbitrary sequences of
%   octets in a form that need not be humanly readable.  A 65-character
%   subset ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be
%   represented per printable character.
%
%   See also BASE64DECODE.
%
%   Written by Jarek Tuszynski, SAIC, jaroslaw.w.tuszynski_at_saic.com
%
%   Matlab version based on 2004 code by Peter J. Acklam
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam
%   http://home.online.no/~pjacklam/matlab/software/util/datautil/base64encode.m

if nargin<2, alg='java';      end
if nargin<3, isChunked=false; end
if ~islogical(isChunked)
  if isnumeric(isChunked)
    isChunked=(isChunked>0);
  else
    isChunked=false;
  end
end
if nargin<4, url_safe=false; end
if ~islogical(url_safe)
  if isnumeric(url_safe)
    url_safe=(url_safe>0);
  else
    url_safe=false;
  end
end

%% if x happen to be a filename than read the file
if (numel(x)<256)
  if (exist(x, 'file')==2)
    fid = fopen(x,'rb');
    x = fread(fid, 'uint8');             % read image file as a raw binary
    fclose(fid);
  end
end

%% Perform conversion
switch (alg)
  case 'java'
    base64 = org.apache.commons.codec.binary.Base64;
    y = base64.encodeBase64(x, isChunked);
    if url_safe
      y = strrep(y,'=','-');
      y = strrep(y,'/','_');
    end

  case 'matlab'

    %% add padding if necessary, to make the length of x a multiple of 3
    x   = uint8(x(:));
    ndbytes = length(x);                 % number of decoded bytes
    nchunks = ceil(ndbytes / 3);         % number of chunks/groups
    if rem(ndbytes, 3)>0
      x(end+1 : 3*nchunks) = 0;          % add padding
    end
    x = reshape(x, [3, nchunks]);        % reshape the data
    y = repmat(uint8(0), 4, nchunks);    % for the encoded data

    %% Split up every 3 bytes into 4 pieces
    %    aaaaaabb bbbbcccc ccdddddd
    % to form
    %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
    y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)
    y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest  bits of x(1,:)
    y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)
    y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest  bits of x(2,:)
    y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)
    y(4,:) = bitand(x(3,:), 63);                    % 6 lowest  bits of x(3,:)

    %% Perform the mapping
    %   0  - 25  ->  A-Z
    %   26 - 51  ->  a-z
    %   52 - 61  ->  0-9
    %   62       ->  +
    %   63       ->  /
    map = ['A':'Z', 'a':'z', '0':'9', '+/'];
    if (url_safe), map(63:64)='-_'; end
    y = map(y(:)+1);

    %% Add padding if necessary.
    npbytes = 3 * nchunks - ndbytes;    % number of padding bytes
    if npbytes>0
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
    end

    %% break into lines with length LineLength
    if (isChunked)
      eol = sprintf('\n');
      nebytes = numel(y);
      nlines  = ceil(nebytes / 76);     % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = [];
    end
end

%% reshape to a row vector and make it a character array
y = char(reshape(y, 1, numel(y)));

author: TadasBaltrusaitis

repo: OpenFace


text:
function gamblers_ruin_plot ( a_stakes, b_stakes )

%*****************************************************************************80
%
%% GAMBLERS_RUIN_PLOT plots a game of gambler's ruin.
%
%  Discussion:
%
%    Two players, A and B, repeatedly toss a coin.
%    For heads, A wins one dollar from B;
%    For tails, B wins one dollar from A.
%    Play continues until one player is bankrupt.
%
%    The program simulates the game, and then produces a plot of the
%    amount of money that A has at each stage of the game.  At the end,
%    A must have either nothing or all the money.
%
%    The program produces a plot of A's money.  It also displays the
%    initial stakes, the number of steps, and the number of times the
%    lead changed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A_STAKES, B_STAKES, the number of dollars that A and
%    B have initially.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAMBLERS_RUIN_PLOT\n' );
  fprintf ( 1, '  MATLAB version\n' );

  step_num = 0;
  leader = '0';
  flip_num = -1;
  a = a_stakes;
  b = b_stakes;
  value(1) = a;

  while ( 0 < a & 0 < b )

    step_num = step_num + 1;

    r = rand ( );

    if ( r <= 0.5 )
      a = a + 1;
      b = b - 1;
    else
      a = a - 1;
      b = b + 1;
    end

    fprintf ( 1, ' %d %d\n', a, b );
    value(step_num+1) = a;

    if ( a_stakes < a & leader ~= 'A' )
      leader = 'A';
      flip_num = flip_num + 1;
    elseif ( a < a_stakes & leader ~= 'B' )
      leader = 'B';
      flip_num = flip_num + 1;
    end

  end
%
%  Plot the results.
%
  clf
  hold on
  plot ( [ 0, step_num], [ a_stakes, a_stakes ], 'r-', 'LineWidth', 2 );
  plot ( [ 0, step_num], [ 0, 0 ], 'r-', 'LineWidth', 2 );
  plot ( [ 0, step_num], [ a_stakes + b_stakes, a_stakes + b_stakes ], 'r-', 'LineWidth', 2 );

  steps = 0 : step_num;
  plot ( steps, value, 'b-', 'LineWidth', 2 );

  title_string = sprintf ( 'Gambler''s Ruin - A = %d, B = %d, Steps = %d, Flips = %d', ...
    a_stakes, b_stakes, step_num, flip_num );

  title ( title_string );
  xlabel ( 'Coin Tosses' )
  ylabel ( 'A''s Money' );
  axis tight

  hold off

  return
end


text:
function varargout = sigmoid(varargin)
% VL_SIGMOID Sigmoid function
%   Y = VL_SIGMOID(X) returns
%
%    Y = 1 ./ (1 + EXP(X)) ;
%
%   Remark::
%     Useful properties of the sigmoid function are:
%
%     -  1 - VL_SIGMOID(X) = VL_SIGMOID(-X)
%     -  Centered sigmoid: 2 * VL_SIGMOID(X) - 1 ;
%     -  VL_SIGMOID(X) = (EXP(X/2) - EXP(X/2)) / (EXP(X/2) + EXP(X/2))
%
%   See also: VL_DSIGMOID(), VL_HELP().
[varargout{1:nargout}] = vl_sigmoid(varargin{:});

author: yihui-he

repo: panorama

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


text:
%compute dcentralheadmag
function [data,units]=compute_dcentralheadmag(trx,n)

% inputfilename=[outputfolder,'centralheadmag.mat'];
% if ~exist(inputfilename,'file')
%     [trx]=compute_centralheadmag(trx,outputfolder);
% end
% load([outputfolder,'centralheadmag.mat'], 'data')
% centralheadmag=data;
%

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dcentralheadmag=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dcentralheadmag{1,i}=(trx(larva).centralheadmag(2:end)-trx(larva).centralheadmag(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dcentralheadmag;
% filename=[outputfolder, 'dcentralheadmag.mat'];
% save(filename, 'data', 'units')

author: kristinbranson

repo: JAABA


text:
function [cvx_optval,P,q,r,X,lambda] = cheb(A,b,Sigma);

% Computes Chebyshev lower bounds on probability vectors
%
% Calculates a lower bound on the probability that a random vector
% x with mean zero and covariance Sigma satisfies A x <= b
%
% Sigma must be positive definite
%
% output arguments:
% - prob: lower bound on probability
% - P,q,r: x'*P*x + 2*q'*x + r is a quadratic function
%   that majorizes the 0-1 indicator function of the complement
%   of the polyhedron,
% - X, lambda:  a discrete distribution with mean zero, covariance
%   Sigma and Prob(X not in C)  >= 1-prob

%
% maximize  1 - Tr Sigma*P - r
% s.t.      [ P  q     ]             [ 0      a_i/2 ]
%           [ q' r - 1 ] >= tau(i) * [ a_i'/2  -b_i ], i=1,...,m
%           taui >= 0
%           [ P q  ]
%           [ q' r ] >= 0
%
% variables P in Sn, q in Rn, r in R
%

[ m, n ] = size( A );
cvx_begin sdp quiet
    variable P(n,n) symmetric
    variables q(n) r tau(m)
    dual variables Z{m}
    maximize( 1 - trace( Sigma * P ) - r )
    subject to
        for i = 1 : m,
            qadj = q - 0.5 * tau(i) * A(i,:)';
            radj = r - 1 + tau(i) * b(i);
            [ P, qadj ; qadj', radj ] >= 0 : Z{i};
        end
        [ P, q ; q', r ] >= 0;
        tau >= 0;
cvx_end

if nargout < 4,
    return
end

X = [];
lambda = [];
for i=1:m
   Zi = Z{i};
   if (abs(Zi(3,3)) > 1e-4)
      lambda = [lambda; Zi(3,3)];
      X = [X Zi(1:2,3)/Zi(3,3)];
   end;
end;
mu = 1-sum(lambda);
if (mu>1e-5)
   w = (-X*lambda)/mu;
   W = (Sigma - X*diag(lambda)*X')/mu;
   [v,d] = eig(W-w*w');
   d = diag(d);
   s = sum(d>1e-5);
   if (d(1) > 1e-5)
      X = [X w+sqrt(s)*sqrt(d(1))*v(:,1) ...
            w-sqrt(s)*sqrt(d(1))*v(:,1)];
      lambda = [lambda; mu/(2*s); mu/(2*s)];
   elseif (d(2) > 1e-5)
      X = [X w+sqrt(s)*sqrt(d(2))*v(:,2) ...
            w-sqrt(s)*sqrt(d(2))*v(:,2)];
      lambda = [lambda; mu/(2*s); mu/(2*s)];
   else
      X = [X w];
      lambda = [lambda; mu];
   end;
end;

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [f,Ls]=ifilterbank(c,g,a,varargin)
%IFILTERBANK  Filter bank inversion
%   Usage:  f=ifilterbank(c,g,a);
%
%   `ifilterbank(c,g,a)` synthesizes a signal *f* from the coefficients *c*
%   using the filters stored in *g* for a channel subsampling rate of *a* (the
%   hop-size). The coefficients has to be in the format returned by
%   either |filterbank| or |ufilterbank|.
%
%   The filter format for *g* is the same as for |filterbank|.
%
%   If perfect reconstruction is desired, the filters must be the duals
%   of the filters used to generate the coefficients. See the help on
%   |filterbankdual|.
%
%   Additional parameters
%   ---------------------
%
%   'complex' (default), 'real'
%       The 'real' flag indicates that the filters *g* cover only the positive
%       frequencies and does `2*real(f)` to effectivelly mirror the filters to
%       cover also the negative frequencies.
%
%   See also: filterbank, ufilterbank, filterbankdual
%
%   References: bohlfe02

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'pfilt'};
definput.keyvals.Ls=[];
definput.flags.complex={'complex','real'};
[flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);

L=filterbanklengthcoef(c,a);

if iscell(c)
  M=numel(c);
else
  M=size(c,2);
end;

[g,asan]=filterbankwin(g,a,L,'normal');

if numel(g)~=M
  error(['%s: Number of filters must be equal to the number of channels ' ...
            'of coefficients.'],upper(mfilename));
end

 if size(a,1)>1
   if  size(a,1)~=M
     error(['%s: The number of entries in "a" must match the number of ' ...
            'filters.'],upper(mfilename));
   end;
 end;

g = comp_filterbank_pre(g,asan,L,kv.crossover);

% Handle ufilterbank output format here
if isnumeric(c)
   ctmp = c;
   c = cell(M,1);
   for m=1:M
      c{m}=squeeze(ctmp(:,m,:));
   end;
end

f = comp_ifilterbank(c,g,asan,L);

if flags.do_real
    f = 2*real(f);
end

% Cut or extend f to the correct length, if desired.
if ~isempty(Ls)
  f=postpad(f,Ls);
else
  Ls=L;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function sqnorm = multisqnorm(A)
% Returns the squared Frobenius norms of the slices of a 3D matrix.
%
% function sqnorm = multisqnorm(A)
%
% Given a 3-dimensional matrix A of size n-by-m-by-N, returns a column
% vector of length N such that sqnorm(i) = norm(A(:, :, i), 'fro')^2.
%
% See also: multiprod multitransp multitrace norms

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 17, 2015.
% Contributors:
% Change log:

    assert(ndims(A) <= 3, ...
           ['multisqnorm is only well defined for matrix arrays of 3 ' ...
            'or less dimensions.']);
    [n, m, N] = size(A);

    % This is equivalent to squeeze(sum(norms(A, 2, 1).^2)), but faster.
    sqnorm = sum(reshape(A, n*m, N).^2, 1)';

end

author: NicolasBoumal

repo: manopt


text:
function [psfMatData, p] = constructMatrix_new( PSF, center, boundary, imsize )
%
%       psfMatData = constructMatrix_new( PSF, center, boundary, imsize );
%
%  Construct psfMatrix data.
%
%  Given a PSF and the locatio of the corresponding point source,
%  this function sets up the data needed to do efficient matrix-vector
%  multiplication.
%
%  This is meant only for spatially invariant PSFs.
%
%  Input:
%         PSF  -  PSF image array
%      center  -  location of point source
%    boundary  -  desired boundary condition
%      imsize  -  size of restored image
%
%  Output:
%   psfMatData -  array containing the (complex) data needed to
%                 do efficient matrix-vector multiplications.
%
switch boundary
    case 'periodic'
        padSize = imsize - size(PSF);
        PSF = padarray(PSF, padSize, 'post');
        psfMatData = fft2(circshift(PSF, 1-center));
        p = [];
    otherwise
        [U, S, V] = svd(PSF);
        [m, n] = size(PSF);
        imsize_pad = 2.^nextpow2(2*imsize);
        minU = abs(min(min(U(:,1))));
        maxU = max(max(abs(U(:,1))));
       if minU == maxU
           U = -U;
           V = -V;
       end
       %
       %  Find number of Kronecker terms to represent the matrix.
       %
       Rtol = 10*eps;
       MaxTerms = sum(diag(S)/S(1,1) >= Rtol);

       %
       %  Get quantities needed to multiply by the Toeplitz pieces:
       %
       bidx_top1 = 1:n-center(2)+1;
       bidx_top2 = center(2):n;
       bidx_bot1 = imsize_pad(2)-center(2)+2:imsize_pad(2);
       bidx_bot2 = 1:center(2)-1;
       cidx_top1 = 1:m-center(1)+1;
       cidx_top2 = center(1):m;
       cidx_bot1 = imsize_pad(1)-center(1)+2:imsize_pad(1);
       cidx_bot2 = 1:center(1)-1;

       BT = V(:,1:MaxTerms)*sqrt(S(1:MaxTerms,1:MaxTerms));
       CT = U(:,1:MaxTerms)*sqrt(S(1:MaxTerms,1:MaxTerms));

       BT_hat = zeros(imsize_pad(2), MaxTerms);
       CT_hat = zeros(imsize_pad(1), MaxTerms);

       BT_hat(bidx_top1,:) = BT(bidx_top2,:);
       BT_hat(bidx_bot1,:) = BT(bidx_bot2,:);
       CT_hat(cidx_top1,:) = CT(cidx_top2,:);
       CT_hat(cidx_bot1,:) = CT(cidx_bot2,:);
       lambdaT = fft(CT_hat);
       deltaT = fft(BT_hat);

       %
       %  Get quantities needed to multiply by the Hankel pieces:
       %
       bidx_top1 = 1:n-center(2);
       bidx_top2 = center(2)+1:n;
       bidx_bot1 = imsize_pad(2)-center(2)+1:imsize_pad(2)-1;
       bidx_bot2 = 1:center(2)-1;
       cidx_top1 = 1:m-center(1);
       cidx_top2 = center(1)+1:m;
       cidx_bot1 = imsize_pad(1)-center(1)+1:imsize_pad(1)-1;
       cidx_bot2 = 1:center(1)-1;

       BH_hat = zeros(imsize_pad(2), MaxTerms);
       CH_hat = zeros(imsize_pad(1), MaxTerms);
       BH_hat(bidx_top1,:) = BT(bidx_top2,:);
       BH_hat(bidx_bot1,:) = BT(bidx_bot2,:);
       CH_hat(cidx_top1,:) = CT(cidx_top2,:);
       CH_hat(cidx_bot1,:) = CT(cidx_bot2,:);
       lambdaH = fft(CH_hat);
       deltaH = fft(BH_hat);

       %
       %  Need to permute these with a shift matrix:
       %
       pB = [1;(imsize_pad(2):-1:2)'];
       pC = [1;(imsize_pad(1):-1:2)'];

       lambdaH = lambdaH(pC,:);
       deltaH = deltaH(pB,:);

       psfMatData = zeros(imsize_pad(1), imsize_pad(2), 4);
       psfMatData(:,:,1) = lambdaT*deltaT';
       psfMatData(:,:,2) = lambdaT*deltaH';
       psfMatData(:,:,3) = lambdaH*deltaT';
       psfMatData(:,:,4) = lambdaH*deltaH';
       p = [pB, pC];

end

repo: IRtools


text:
function update_visualization_func = show_video(img_files, video_path, resize_image)
%SHOW_VIDEO
%   Visualizes a tracker in an interactive figure, given a cell array of
%   image file names, their path, and whether to resize the images to
%   half size or not.
%
%   This function returns an UPDATE_VISUALIZATION function handle, that
%   can be called with a frame number and a bounding box [x, y, width,
%   height], as soon as the results for a new frame have been calculated.
%   This way, your results are shown in real-time, but they are also
%   remembered so you can navigate and inspect the video afterwards.
%   Press 'Esc' to send a stop signal (returned by UPDATE_VISUALIZATION).
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/

	%store one instance per frame
	num_frames = numel(img_files);
	boxes = cell(num_frames,1);

	%create window
	[fig_h, axes_h, unused, scroll] = videofig(num_frames, @redraw, [], [], @on_key_press);  %#ok, unused outputs
	set(fig_h, 'Number','off', 'Name', ['Tracker - ' video_path])
	axis off;

	%image and rectangle handles start empty, they are initialized later
	im_h = [];
	rect_h = [];

	update_visualization_func = @update_visualization;
	stop_tracker = false;

	function stop = update_visualization(frame, box)
		%store the tracker instance for one frame, and show it. returns
		%true if processing should stop (user pressed 'Esc').
		boxes{frame} = box;
		scroll(frame);
		stop = stop_tracker;
	end

	function redraw(frame)
		%render main image
		im = imread([video_path img_files{frame}]);
		if size(im,3) > 1,
			im = rgb2gray(im);
		end
		if resize_image,
			im = imresize(im, 0.5);
		end

		if isempty(im_h),  %create image
			im_h = imshow(im, 'Border','tight', 'InitialMag',200, 'Parent',axes_h);
		else  %just update it
			set(im_h, 'CData', im)
		end

		%render target bounding box for this frame
		if isempty(rect_h),  %create it for the first time
			rect_h = rectangle('Position',[0,0,1,1], 'EdgeColor','g', 'Parent',axes_h);
		end
		if ~isempty(boxes{frame}),
			set(rect_h, 'Visible', 'on', 'Position', boxes{frame});
		else
			set(rect_h, 'Visible', 'off');
		end
	end

	function on_key_press(key)
		if strcmp(key, 'escape'),  %stop on 'Esc'
			stop_tracker = true;
		end
	end

end

repo: KCF


text:
function planC = createUniformlySlicedPlanC(planC)
%
% function planC = createUniformlySlicedPlanC(planC)
%
% This function creates uniformly sliced planC
%
% APA, 08/19/2011

global planC

indexS = planC{end};

% Get scan grid
zValuesV = [planC{indexS.scan}.scanInfo(:).zValue];
zDiffV = diff(zValuesV);
if ~all(zDiffV == min(zDiffV))
   % Create new z-vector
   minSpac = min(zDiffV);
   zValNewV = zValuesV(1):minSpac:zValuesV(end);

else
    % Already uniform plan
    return

end

% Loop over new z-values and create new scanArray
numSlices = length(zValNewV);

scanInfoS = planC{indexS.scan}.scanInfo(1);
scanInfoS.sliceThickness = minSpac;
scanInfoS.voxelThickness = minSpac;

for i = 1:numSlices

    % Find the nearest non-uniform slice
    nearestZindexV(i) = findnearest(zValuesV,zValNewV(i));

    scanArrayNewM(:,:,i) = planC{indexS.scan}.scanArray(:,:,nearestZindexV(i));

    scanInfoS.zValue = zValNewV(i);

    newScanInfoS(i) = scanInfoS;

end

planC{indexS.scan}.scanArray = scanArrayNewM;

clear scanArrayNewM;

planC{indexS.scan}.scanInfo = newScanInfoS;

% Change structure z-cords
numStructures = planC{indexS.structures};
for strNum = 1:length(numStructures)

    newStructureS = planC{indexS.structures}(strNum);

    newStructureS.numberOfScans = numSlices;

    newStructureS.contour(:) = [];
    newStructureS.contour
    newStructureS.contour(1:numSlices) = planC{indexS.structures}(strNum).contour(nearestZindexV);

    for slcNum = 1:length(newStructureS.contour)

        for segNum = 1:length(newStructureS.contour(slcNum).segments)

            if ~isempty(newStructureS.contour(slcNum).segments(segNum).points)

                newStructureS.contour(slcNum).segments(segNum).points(:,3) = zValNewV(slcNum);

            end

        end

    end

    planC{indexS.structures}(strNum) = newStructureS;

end

% ReRaster and Uniformize

reRasterAndUniformize;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function x=mdl(m)
%MDL  Returns Rissanen's Minimum Description Length.
% m=model that has been estimated using System Identication toolbox.
% This function requires System Identification toolbox to work.
% Plug-compatible with built-in functions aic(m) and fpe(m).
% MDL can be used like AIC or FPE to compare models of different
% complexities.  Choose model with lowest MDL or AIC or FPE.
% Pintelon & Schoukens (2001) pp. 329,550 say MDL is better
% than AIC; AIC tends to select a too-complex model.
% Example: Compute & print MDL and AIC for an AR model of order 10.
%   Data=iddata(y,[],1/Fs);
%   m_fb=ar(Data,10,'fb');
%   fprintf('MDL=%.3d; AIC=%.3f\n',mdl(m_fb),aic(m_fb));
% William C Rose 2007-06-05.

d=size(m,'Npar');   % d=number of model parameters
N=m.es.DataLength;  % N=number of data points fitted
V=m.es.LossFcn;     % V=loss function;
x=V*(1+d*log(N)/N); % Rissanen's Minimum Description Length

author: Sable

repo: mcbench-benchmarks


text:
function res = arg_nogui(varargin)
% Like arg(), but not displayed by GUIs.
% Spec = arg_nogui(Names,Default,Range,Help,Options...)
%
% This type of function argument specifier behaves like arg(), except that it will not be displayed
% in GUIs that are generated for the function. This is mainly used for optional arguments that have
% a format that is too complex to be meaningfully edited by the user (e.g. a large matrix, or a
% long cell array), or undocumented or possibly confusing (e.g. some functions have special internal
% or reserved arguments).
%
% In:
%   Names : The name(s) of the argument. At least one must be specified, and if multiple are
%           specified, they must be passed in a cell array.
%           * The first name specified is the argument's "code" name, as it should appear in the
%             function's code (= the name under which arg_define() returns it to the function).
%           * The second name, if specified, is the "Human-readable" name, which is exposed in the
%             GUIs (if omitted, the code name is displayed). For consistency with other MATLAB
%             functions it should be in CamelCase.
%           * Further specified names are aliases for the argument (e.g., for backwards
%             compatibility with older function syntaxes/parameter names).
%
%   Default : Optionally the default value of the argument; can be any data structure (default: []).
%
%   Range : Optionally a range of admissible values (default: []).
%           * If empty, no range is enforced.
%           * If a cell array, each cell is considered one of the allowed values.
%           * If a 2-element numeric vector, the two values are considered the numeric range of the
%             data (inclusive).
%
%   Help : The help text for this argument, optional. (default: '').
%
%   Options... : Optional name-value pairs to denote additional properties, same as in arg().
%
% Out:
%   Spec : A cell array, that, when called as feval(Spec{1},reptype,Spec{2}{:}), yields a
%          specification of the argument, for use by arg_define.
%
% Examples:
%   function myfunction(varargin)
%   arg_define(varargin, ...
%       arg('arg1',10,[],'Some argument.'), ...
%       arg_nogui('arg_hidden',1001,[],'Hidden argument. This one will not be displayed in a GUI; reserved for special purposes.'));
%
% See also:
%   arg, arg_norep, arg_sub, arg_subswitch, arg_subtoggle, arg_define
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-09-24

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

if nargin == 1
    res = {'expand_arg',[varargin {[],[],[],'displayable',false}]};
elseif nargin >= 4
    res = {'expand_arg',[varargin(1:4) {'displayable',false} varargin(5:end)]};
elseif nargin == 2
    res = {'expand_arg',[varargin {[],[],'displayable',false}]};
else
    res = {'expand_arg',[varargin {[],'displayable',false}]};
end

repo: REST


text:
function X = rocholForeSub(ch, Y);

% ROCHOLFORESUB Foreward substitute the representation of the rank one Cholesky.

% ROCHOL

%/~
% This would be the long way of doing it.
%L = rocholExtract(ch);
%X2 = L\Y;
%~/

X = zeros(size(Y));
t = zeros(1, size(Y, 2));
X(1, :) = Y(1, :)/ch.s(1);
for i = 2:ch.n
  if i == 2
    t = Y(1, :)*ch.u(1);
  else
    t = t + ch.s(i-1)*ch.u(i-1)*X(i-1, :);
  end
  X(i, :) = (Y(i, :)-ch.v(i)*t)/ch.s(i);
end

%/~
%disp(max(max(X - X2)));
%~/

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function X = zscore(X)

%ZSCORE standardized z-scores for GRIDobj
%
% Syntax
%
%     Z = zscore(X)
%
% Description
%
%     zscore returns the z-score for each element of GRIDobj X such that
%     all values of X are centered to have mean 0 and scaled to have
%     standard deviation 1.
%
% Input arguments
%
%     X     GRIDobj
%
% Output arguments
%
%     Y     GRIDobj
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     Z = zscore(DEM);
%     imagesc(Z); colorbar
%
% See also: GRIDobj, zscore
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 15. March, 2016

I = ~isnan(X.Z);
X.Z(I) = (X.Z(I(:)) - mean(X.Z(I(:))))./std(X.Z(I(:)));

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function [ t_start, y_start ] = p18_start ( neqn )

%*****************************************************************************80
%
%% P18_START returns the starting point for problem p18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  delta = p18_param ( 'GET', 'DELTA', [] );
  y_start(1) = 1.0 - delta;
  y_start(2) = 0.0;
  y_start(3) = 0.0;
  y_start(4) = sqrt ( ( 1.0 + delta ) / ( 1.0 - delta ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ nr, nt, nc ] = circle_rt_size ( rule )

%*****************************************************************************80
%
%% CIRCLE_RT_SIZE sizes an R, THETA product quadrature rule in the unit circle.
%
%  Discussion:
%
%    For a given value of RULE, here are the number of points used at the
%    center (NC), the number of points along the radial direction (NR) and
%    the number of points along the theta direction (NT).  The total number
%    of points in the rule will be
%
%      Total = NC + NR * NT.
%
%    The user, when choosing RULE, must allocate enough space in the arrays
%    RA, RW, TA and TW for the resulting values of NR and NT.
%
%    RULE  NC  NR  NT  Total
%    ----  --  --  --  -----
%       1   1   0   0      1
%       2   0   1   4      4
%       3   1   1   4      5
%       4   1   1   6      7
%       5   1   2   4      9
%       6   0   3   4     12
%       7   1   2  10     21
%       8   0   4  16     64
%       9   0   5  20    120
%
%    The integral of F(X,Y) over the unit circle is approximated by
%
%      Integral ( X*X + Y*Y <= 1 ) F(X,Y) dx dy
%      = Integral ( 0 <= R <= 1, 0 <= T <= 2PI ) F(R*cos(T),R*sin(T)) r dr dt
%      = approximately
%        ZW * F(0,0)
%        + sum ( 1 <= I <= NR ) Sum ( 1 <= J <= NT )
%        RW(I) * TW(J) * F ( R(I) * cos ( TA(J) ), R(I) * sin ( TA(J) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer RULE, the rule desired.
%
%    Output, integer NR, the number of R abscissas.
%
%    Output, integer NT, the number of Theta abscissas.
%
%    Output, integer NC, the number of center abscissas (0 or 1).
%
  if ( rule == 1 )

    nr = 0;
    nt = 0;
    nc = 1;

  elseif ( rule == 2 )

    nr = 1;
    nt = 4;
    nc = 0;

  elseif ( rule == 3 )

    nr = 1;
    nt = 4;
    nc = 1;

  elseif ( rule == 4 )

    nr = 1;
    nt = 6;
    nc = 1;

  elseif ( rule == 5 )

    nr = 2;
    nt = 4;
    nc = 1;

  elseif ( rule == 6 )

    nr = 3;
    nt = 4;
    nc = 0;

  elseif ( rule == 7 )

    nr = 2;
    nt = 10;
    nc = 1;

  elseif ( rule == 8 )

    nr = 4;
    nt = 16;
    nc = 0;

  elseif ( rule == 9 )

    nr = 5;
    nt = 20;
    nc = 0;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE_RT_SIZE - Fatal error!\n' );
    fprintf ( 1, '  There is no rule of index %d\n', rule );
    error ( 'CIRCLE_RT_SIZE - Fatal error!' );

  end

  return
end


text:
function infsup(c)
%INFSUP       Display of interval taylor in infsup notation
%
%   infsup(c)
%

% written  05/21/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/07/12     S.M. Rump  complete redesign
%

  INTLAB_TAYLOR_ORDER = getappdata(0,'INTLAB_TAYLOR_ORDER');

  loose = strcmp(get(0,'FormatSpacing'),'loose');

  name = inputname(1);
  if isempty(name)                    % happens for display(taylorinit(random))
    name = 'ans';
  end

  numvar = size(c.t,1)-1;
  if numvar~=INTLAB_TAYLOR_ORDER
    warning('**** number of dependent variables and partial derivatives do not coincide')
  end

  INTLAB_INTVAL_DISPLAY = getappdata(0,'INTLAB_INTVAL_DISPLAY');
  setappdata(0,'INTLAB_INTVAL_DISPLAY','DisplayInfsup');
  display(c,name)
  setappdata(0,'INTLAB_INTVAL_DISPLAY',INTLAB_INTVAL_DISPLAY);

repo: OpenMAS


text:
% -----------------------------------------------------------------  %
% Matlab Programs included the Appendix B in the book:               %
%  Xin-She Yang, Engineering Optimization: An Introduction           %
%                with Metaheuristic Applications                     %
%  Published by John Wiley & Sons, USA, July 2010                    %
%  ISBN: 978-0-470-58246-6,   Hardcover, 347 pages                   %
% -----------------------------------------------------------------  %
% Citation detail:                                                   %
% X.-S. Yang, Engineering Optimization: An Introduction with         %
% Metaheuristic Application, Wiley, USA, (2010).                     %
%                                                                    %
% http://www.wiley.com/WileyCDA/WileyTitle/productCd-0470582464.html %
% http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0470582464.html  %
% -----------------------------------------------------------------  %
% ===== ftp://  ===== ftp://   ===== ftp:// =======================  %
% Matlab files ftp site at Wiley                                     %
% ftp://ftp.wiley.com/public/sci_tech_med/engineering_optimization   %
% ----------------------------------------------------------------   %

% Design Optimization of a Pressure Vessel using fmincon             %
% This is essentially a different problem as integer multiples       %
% are not applied.  If leaves an exercise to modify the program that %
% the first two design variables are the integer multiples of 0.0625 %
% ------------------------------------------------------------------ %

function B7_pressure

  d=0.0625;
  % Linear inequality constraints
  A=[-1 0 0.0193 0;0 -1 0.00954 0;0 0 0 1];
  b=[0; 0; 240];
  % Simple bounds
  Lb=[d; d; 10; 10];
  Ub=[99*d; 99*d; 200; 200];
  x0=Lb+(Ub-Lb).*rand(size(Lb));
  options=optimset('Display','iter','TolFun',1e-08);
  [x,fval]=fmincon(@objfun,x0,A,b,[],[],Lb,Ub,@nonfun,options)

% The objective function
function f=objfun(x)
f=0.6224*x(1)*x(3)*x(4)+1.7781*x(2)*x(3)^2 ...
  +3.1661*x(1)^2*x(4)+19.84*x(1)^2*x(3)

% Nonlinear constraints
function [g,geq]=nonfun(x)
% Nonlinear inequality
  g=-pi*x(3)^2*x(4)-4*pi/3*x(3)^3+1296000;
% Equality constraint [none]
  geq=[];


text:
function q=rotMat2Quat(R,handed)
%%ROTMAT2QUAT Get a unit quaternion corresponding to a particular rotation
%             matrix. The quaternion can be chosen to support standard
%             right-handed quaternion multiplication rules, or non-standard
%             left-handed rules that some authors choose to use.
%
%INPUTS:  M A 3X3 orthonormal real rotation matrix.
%    handed The handedness of the quaternion. If omitted, it is assumed
%           that the quaternion is right-handed (the standard). Possible
%           values are:
%           'right' The default if omitted. The quaternion multiplication
%                   is assumed right-handed (standard).
%           'left'  The quaternion multiplication is assumed left-handed.
%                   This is used in someplaces, including the reference
%                   from Shuster, below.
%
%OUTPUTS: q A 4X1 unit quaternion corresponding to the rotation matrix. The
%           quaternion is ordered [cos(theta/2);sin(theta/2)u'] where u is
%           a unit vector for the axis of rotation and theta is the
%           counterclockwise (right-handed) or clockwise (left-handed)
%           rotation angle about that unit vector.
%
%The formulae for converting a rotation matrix are from [1], where a minor
%change has been performed to support both right and left-handed
%quaternions. Both q and -q represent the same rotations. Here, the
%solution is chosen so that when considering a left-handed rotation, the
%sign of the largest magnitude element is positive. When considering a
%right handed rotation, the sign of the largest element may or may not be
%positive.
%
%A quaternion of form q(1)+i*q(2)+j*q(3)+k*q(4) that obeys right-handed
%multiplication rules supports the following rules for multiplication of i,
%j, and k, where an item in a row is multiplied by an item in the column to
%get the result:
%  i,  j, k
%i -1, k,-j
%j -k,-1, i
%k  j,-i,-1
%On the other hand, left-handed multiplication rules flip the signs of the
%off-diagonal terms:
%  i,  j, k
%i -1,-k, j
%j  k,-1,-i
%k -j, i,-1
%
%REFERENCES:
%[1] W. F. Phillips, C. E. Hailey, and G. A. Gebert, "Review of attitude
%    representations used for aircraft kinematics," Journal of Aircraft,
%    vol. 38, no. 4, pp. 718-737, Jul. - Aug. 2001.
%
%August 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(handed))
    handed='right';
end

R11=R(1,1);
R22=R(2,2);
R33=R(3,3);
R12=R(1,2);
R21=R(2,1);
R13=R(1,3);
R31=R(3,1);
R23=R(2,3);
R32=R(3,2);

%The square of the quaternion...
q2=(1/4)*[1+R11+R22+R33;
          1+R11-R22-R33;
          1-R11+R22-R33;
          1-R11-R22+R33];

q2Max=max(q2);

if(q2Max==q2(1))
    q0=0.5*sqrt(1+R11+R22+R33);
    q1=(1/(4*q0))*(R23-R32);
    q2=(1/(4*q0))*(R31-R13);
    q3=(1/(4*q0))*(R12-R21);
elseif(q2Max==q2(2))
    q1=0.5*sqrt(1+R11-R22-R33);
    q0=(1/(4*q1))*(R23-R32);
    q2=(1/(4*q1))*(R12+R21);
    q3=(1/(4*q1))*(R31+R13);
elseif(q2Max==q2(3))
    q2=0.5*sqrt(1-R11+R22-R33);
    q0=(1/(4*q2))*(R31-R13);
    q1=(1/(4*q2))*(R12+R21);
    q3=(1/(4*q2))*(R23+R32);
else
    q3=0.5*sqrt(1-R11-R22+R33);
    q0=(1/(4*q3))*(R12-R21);
    q1=(1/(4*q3))*(R31+R13);
    q2=(1/(4*q3))*(R23+R32);
end

%The above implementation provides a left-handed quaternion, so the sign of
%the vector part needs to be flipped if the standard right-handed system is
%used.
switch(handed)
    case 'right'
        q=[q0;-q1;-q2;-q3];
    case 'left'
        q=[q0;q1;q2;q3];
    otherwise
        error('Invalid handedness provided.')
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Plots the Lagrangian structure with:
%           (1): the background velocity field
%           (2): the Lagrangian structure itself
%           (3): the vorticity field in a colormap
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [loc, diffy] =  please_Plot_Results(ds,X,Y,U,V,vort,uMag,p,C,chiX,chiY,lagPlot,velPlot,vortPlot,pressPlot,uMagPlot,conPlot,firstPrint,loc,diffy,spacing)

%X,Y:  (x,y) values
%U,V:  x-directed, y-directed velocities respectively
%vort: vorticity
%uMag: magnitude of velocity
%p:    pressure

%FLAGS FOR PLOTTING:
%pressPlot - if you want pressure plot = 1
%uMagPlot  - if you want mag. velocity plot = 1
%vortPlot  - if you want vorticity plot = 1
%velPlot   - if you want velocity plot = 1
%lagPlot   - if you want lag. point ONLY plot = 1
%conPlot   - if you want concentration plot =1
%
% Assumption: Assuming chiX and chiY are column vectors
% Assumption: Assuming chiX(i+1)-chiX(i) < .5 and chiY(i+1)-chiY(i) < .5, for all points that don't cross the boundary
%
% This script has been adapted from the Mat-IB code by Hiens and Stockie, 2014.
%

Lx = X(1,end)+X(1,2);
Ly = Y(end,1)+Y(2,1);

%Shift so inside of interval [0,Lx] or [0,Ly]
chiX = mod(chiX,Lx);
chiY = mod(chiY,Ly);

%Connect Geometry Appropriately At Beginning Of Simulation then Store it!
if firstPrint == 1
    diffX = abs([chiX(2:end);chiX(1)]-chiX);
    diffY = abs([chiY(2:end);chiY(1)]-chiY);

    locX = find(diffX > spacing ); % assuming x value can't change more than Lx/2 (without crossing boundary)
    locY = find(diffY > spacing ); % assuming y value can't change more than Ly/2 (without crossing boundary)
    loc =  sort(unique([locX;locY]));

    diffy = sqrt( (chiX(end)-chiX(1) )^2 + ( chiY(end)-chiY(1) )^2 );
end

clf; %Clear previous plots :)

[N1,N2,num_con]=size(C);
figure(1)
numPlots = lagPlot+velPlot+vortPlot+pressPlot+uMagPlot+conPlot+(num_con-1);

ct = 1;

% % % % % PLOTS LAGRANGIAN POINTS ONLY (if selected) % % % % %

if lagPlot == 1
    subplot(1,numPlots,ct)
    axis([0 Lx 0 Ly]);
    title('LAGRANGIAN PTS');
    xlabel('x'); ylabel('y');
    hold all;

    loc = [0;loc;length(chiX)];
    for i=2:length(loc)
       plot(chiX(loc(i-1)+1:loc(i)),chiY(loc(i-1)+1:loc(i)),'m','LineWidth',3);
    end
    if diffy < 5*ds
        xTemp = [chiX(1) chiX(end)];
        yTemp = [chiY(1) chiY(end)];
        plot(xTemp(1:2),yTemp(1:2),'m','LineWidth',3);
    end

    axis square;

    ct=ct+1;
end

% % % % % PLOTS VORTICITY + LAG Pts. (if selected) % % % % %

if vortPlot == 1

    subplot(1,numPlots,ct)
    %
    axis([0 Lx 0 Ly]);
    title('VORTICITY');
    xlabel('x'); ylabel('y');
    hold all;

    %Compute Vorticity and Plot It against Lagrangian Grid!
    x = X(1,:); y = Y(:,1);
    contourf(x,y,flipud(rot90(vort)),10); hold on;

    loc = [0;loc;length(chiX)];
    for i=2:length(loc)
        plot(chiX(loc(i-1)+1:loc(i)),chiY(loc(i-1)+1:loc(i)),'m','LineWidth',3);
    end
    if diffy < 5*ds
        xTemp = [chiX(1) chiX(end)];
        yTemp = [chiY(1) chiY(end)];
        plot(xTemp(1:2),yTemp(1:2),'m','LineWidth',3);
    end

    axis square;

    ct=ct+1;
end

% % % % % PLOTS PRESSURE + LAG Pts. (if selected) % % % % %

if pressPlot == 1

    subplot(1,numPlots,ct)
    %
    axis([0 Lx 0 Ly]);
    title('PRESSURE');
    xlabel('x'); ylabel('y');
    hold all;

    %Use Pressure and Plot It against Lagrangian Grid!
    x = X(1,:); y = Y(:,1);
    contourf(x,y,p,6); hold on;

    loc = [0;loc;length(chiX)];
    for i=2:length(loc)
        plot(chiX(loc(i-1)+1:loc(i)),chiY(loc(i-1)+1:loc(i)),'m','LineWidth',3);
    end
    if diffy < 5*ds
        xTemp = [chiX(1) chiX(end)];
        yTemp = [chiY(1) chiY(end)];
        plot(xTemp(1:2),yTemp(1:2),'m','LineWidth',3);
    end

    axis square;

    ct=ct+1;
end

% % % % % PLOTS MAGNITUDE OF VELOCITY + LAG Pts. (if selected) % % % % %

if uMagPlot == 1

    subplot(1,numPlots,ct)
    %
    axis([0 Lx 0 Ly]);
    title('MAGNITUDE OF VELOCITY');
    xlabel('x'); ylabel('y');
    hold all;

    %Use Mag. Velocity and Plot It against Lagrangian Grid!
    x = X(1,:); y = Y(:,1);
    contourf(x,y,uMag,6); hold on;

    loc = [0;loc;length(chiX)];
    for i=2:length(loc)
        plot(chiX(loc(i-1)+1:loc(i)),chiY(loc(i-1)+1:loc(i)),'m','LineWidth',3);
    end
    if diffy < 5*ds
        xTemp = [chiX(1) chiX(end)];
        yTemp = [chiY(1) chiY(end)];
        plot(xTemp(1:2),yTemp(1:2),'m','LineWidth',3);
    end

    axis square;

    ct=ct+1;
end

% % % % % PLOTS VELOCITY FIELD + LAG Pts. (if selected) % % % % %

if velPlot == 1

    subplot(1,numPlots,ct)
    %
    axis([0 Lx 0 Ly]);
    title('VELOCITY');
    xlabel('x'); ylabel('y');
    hold all;

    quiver(X,Y,U,V); %Print Velocity Field

    loc = [0;loc;length(chiX)];
    for i=2:length(loc)
       plot(chiX(loc(i-1)+1:loc(i)),chiY(loc(i-1)+1:loc(i)),'m','LineWidth',3);
    end
    if diffy < 5*ds
        xTemp = [chiX(1) chiX(end)];
        yTemp = [chiY(1) chiY(end)];
        plot(xTemp(1:2),yTemp(1:2),'m','LineWidth',3);
    end

    axis square;

    ct=ct+1;
end

if conPlot == 1
    % plotting all concentrations
    for ic=1:num_con
        subplot(1,numPlots,ct)
    %
        axis([0 Lx 0 Ly]);
        title(sprintf('Concentration%g',ic));
        xlabel('x'); ylabel('y');
        hold all;

        x = X(1,:); y = Y(:,1);
        pcolor(X,Y,C(:,:,ic));
        shading interp;

        loc = [0;loc;length(chiX)];
        for i=2:length(loc)
            plot(chiX(loc(i-1)+1:loc(i)),chiY(loc(i-1)+1:loc(i)),'m','LineWidth',3);
        end
        if diffy < 5*ds
            xTemp = [chiX(1) chiX(end)];
            yTemp = [chiY(1) chiY(end)];
            plot(xTemp(1:2),yTemp(1:2),'m','LineWidth',3);
        end

        axis square;

        ct=ct+1;
    end
end

drawnow;

hold off;
set(gca,'Box','on');

author: nickabattista