text:
%%***************************************************************
%% linsysolve: solve linear system to get dy, and direction
%%             corresponding to unrestricted variables.
%%
%% [xx,coeff,L,resnrm] = linsysolve(schur,UU,Afree,EE,rhs);
%%
%% child functions: symqmr.m, mybicgstable.m, linsysolvefun.m
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***************************************************************

   function [xx,coeff,L,resnrm] = linsysolve(par,schur,UU,Afree,EE,rhs);

    global solve_ok existspcholsymb exist_analytic_term
    global nnzmat nnzmatold matfct_options matfct_options_old use_LU
    global Lsymb msg diagR diagRold

    spdensity  = par.spdensity;
    printlevel = par.printlevel;
    iter       = par.iter;
    err        = max([par.relgap,par.pinfeas,par.dinfeas]);

    m = length(schur);
    if (iter==1);
       use_LU = 0; matfct_options_old = '';
       diagR = ones(m,1);
    end
    if isempty(nnzmatold); nnzmatold = 0; end
    diagRold = diagR;
%%
%% schur = schur + rho*diagschur + lam*AAt
%%
    diagschur = abs(full(diag(schur)));
    if (par.ublksize)
       minrho = 1e-15;
    else
       minrho = 1e-17;
    end
    const = 1e-2;
    ratio = max(diagR)/min(diagR);
    if (par.depconstr) | (ratio > 1e10) | (iter < 5)
       %% important: do not perturb beyond certain threshold
       %% since it will adversely affect prim_infeas of fp43
       %%
       rho = max(minrho,min(1e-10,const*norm(par.rp)/(1+norm(diagschur.*par.dy))));
       if (exist_analytic_term); rho = 0; end
       pertdiag = rho*max(1e-8,diagschur);
       mexschurfun(schur,pertdiag);
       %%if (printlevel>2); fprintf(' %2.1e',rho); end
    end
    if (par.depconstr) | (par.ZpATynorm > 1e10) | (par.ublksize) | (iter < 10)
       %% Note: do not add this perturbation even if ratio is large.
       %% It adversely affects hinf15.
       %%
       lam = 0.1*min(1e-14,const*norm(par.rp)/(1+norm(par.AAt*par.dy)));
       if (exist_analytic_term); lam = 0; end
       mexschurfun(schur,lam*par.AAt);
       %%if (printlevel>2); fprintf(' *%2.1e ',lam); end
    end
    if (max(diagschur)/min(diagschur) > 1e14) & (par.blkdim(2) == 0) ...
       & (iter > 10)
       tol = 1e-8;
       idx = find(diagschur < tol); len = length(idx);
       pertdiagschur = zeros(m,1);
       if (len > 0 & len < 5) & (norm(rhs(idx)) < tol)
          pertdiagschur(idx) = 1*ones(length(idx),1);
          mexschurfun(schur,pertdiagschur);
          if (printlevel>2); fprintf('#'); end
       end
    end
%%
%% assemble coefficient matrix
%%
    len = size(Afree,2);
    if ~isempty(EE)
       EE(:,[1 2]) = len + EE(:,[1 2]); %% adjust for ublk
    end
    EE = [(1:len)' (1:len)' zeros(len,1); EE];
    if isempty(EE)
       coeff.mat22 = [];
    else
       coeff.mat22 = spconvert(EE);
    end
    if (size(Afree,2) | size(UU,2))
       coeff.mat12 = [Afree, UU];
    else
       coeff.mat12 = [];
    end
    coeff.mat11 = schur; %% important to use perturbed schur matrix
    ncolU = size(coeff.mat12,2);
%%
%% pad rhs with zero vector
%% decide which solution methods to use
%%
    rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
    if (ncolU > 300); use_LU = 1; end
%%
%% Cholesky factorization
%%
    L = []; resnrm = norm(rhs); xx = inf*ones(m,1);
    if (~use_LU)
       solve_ok = 1;  solvesys = 1;
       nnzmat = mexnnz(coeff.mat11);
       nnzmatdiff = (nnzmat ~= nnzmatold);
       if (nnzmat > spdensity*m^2) | (m < 500)
          matfct_options = 'chol';
       else
          if (par.matlabversion >= 7.3) %% & (par.computer == 64)
             matfct_options = 'spchol';
	  else
             matfct_options = 'myspchol';
          end
       end
       if (printlevel>2); fprintf(' %s ',matfct_options); end
       if strcmp(matfct_options,'chol')
          if issparse(schur); schur = full(schur); end;
          if (iter<=5); %% to fix strange anonmaly in Matlab
             mexschurfun(schur,1e-20,2);
          end
          L.matfct_options = 'chol';
          L.perm = [1:m];
          [L.R,indef] = chol(schur);
          diagR = diag(L.R).^2;
          if (indef)
             diagR = diagRold;
 	     solve_ok = -2; solvesys = 0;
             msg = 'linsysolve: Schur complement matrix not pos. def.';
             if (printlevel); fprintf('\n  %s',msg); end
          end
       elseif strcmp(matfct_options,'spchol')
          if ~issparse(schur); schur = sparse(schur); end;
          L.matfct_options = 'spchol';
          [L.R,indef,L.perm] = chol(schur,'vector');
          L.Rt = L.R';
          diagR = full(diag(L.R)).^2;
          if (indef)
             diagR = diagRold;
 	     solve_ok = -2; solvesys = 0;
             msg = 'linsysolve: Schur complement matrix not pos. def.';
             if (printlevel); fprintf('\n  %s',msg); end
          end
       elseif strcmp(matfct_options,'myspchol')
          if ~issparse(schur), schur = sparse(schur); end;
          if (nnzmatdiff | ~strcmp(matfct_options,matfct_options_old))
             [Lsymb,flag] = symbcholfun(schur,par.cachesize);
             if (flag)
                solve_ok = -2;  solvesys = 0;
                existspcholsymb = 0;
                msg = 'linsysolve: symbolic factorization fails';
                if (printlevel); fprintf('\n  %s',msg); end
                use_LU = 1;
             else
                existspcholsymb = 1;
             end
          end
          if (existspcholsymb)
             L = sparcholfun(Lsymb,schur);
             L.matfct_options  = 'myspchol';
             L.d(find(L.skip)) = 1e20;
             if any(L.skip) & (ncolU)
                solve_ok = -3; solvesys = 0;
                existspcholsymb = 0;
                use_LU = 1;
                if (printlevel)
                   fprintf('\n  L.skip exists but ncolU > 0.');
                   fprintf('\n  switch to LU factor.');
                end
             end
          end
       end
       if (solvesys)
          if (ncolU)
             tmp = coeff.mat12'*linsysolvefun(L,coeff.mat12)-coeff.mat22;
	     if issparse(tmp); tmp = full(tmp); end
             tmp = 0.5*(tmp + tmp');
             [L.Ml,L.Mu,L.Mp] = lu(tmp);
             tol = 1e-16;
             condest = max(abs(diag(L.Mu)))/min(abs(diag(L.Mu)));
             idx = find(abs(diag(L.Mu)) < tol);
             if ~isempty(idx) | (condest > 1e18);
                solvesys = 0; solve_ok = -4;
                use_LU = 1;
                msg = 'SMW too ill-conditioned, switch to LU factor';
                if (printlevel); fprintf('\n  %s.',msg); end
             end
          end
          [xx,resnrm,solve_ok] = symqmr(coeff,rhs,L,[],[],printlevel);
          if (solve_ok <= 0.3) & (printlevel)
             fprintf('\n  warning: symqmr failed: %3.1f ',solve_ok);
          end
       end
       if (solve_ok <= 0.3)
          tol = 1e-10;
          if (m < 7e3 & strcmp(matfct_options,'chol') & (err > tol)) ...
             | (m < 2e5 & strcmp(matfct_options,'spchol') & (err > tol)) ...
             | (m < 2e5 & strcmp(matfct_options,'myspchol') & (err > tol))
             use_LU = 1;
             if (printlevel); fprintf('\n  switch to LU factor.'); end
          end
       end
    end
%%
%% LU factorization
%%
    if (use_LU)
       nnzmat = mexnnz(coeff.mat11)+mexnnz(coeff.mat12);
       nnzmatdiff = (nnzmat ~= nnzmatold);
       solve_ok = 1; solvesys = 1;
       if ~isempty(coeff.mat22)
          raugmat = [coeff.mat11, coeff.mat12; coeff.mat12', coeff.mat22];
       else
          raugmat = coeff.mat11;
       end
       if (nnzmat > spdensity*m^2) | (m+ncolU < 500)
          matfct_options = 'lu';
       else
          matfct_options = 'splu';
       end
       if (printlevel>2); fprintf(' %s ',matfct_options); end
       if strcmp(matfct_options,'lu')
          if issparse(raugmat); raugmat = full(raugmat); end
          [L.l,L.u,L.p] = lu(raugmat);
          L.matfct_options = 'lu';
          L.p = sparse(L.p);
          idx = find(abs(diag(L.u)) < 1e-30);
          if ~isempty(idx)
             msg = 'linsysolve: matrix is singular';
             if (printlevel); fprintf('\n  %s',msg); end
             solvesys = 0;
          end
          [ii,jj] = find(L.p); [dummy,idx] = sort(ii); L.perm = jj(idx);
       end
       if strcmp(matfct_options,'splu')
          if ~issparse(raugmat); raugmat = sparse(raugmat); end
          L.perm = [1:length(raugmat)];
          L.matfct_options = 'splu';
          L.symmatrix = 1;
          [L.l,L.u,L.p,L.q] = lu(raugmat);
          [ii,jj] = find(L.p); L.p = ii;
          [ii,jj] = find(L.q); L.q = ii;
       end
       if (solvesys)
          [xx,resnrm,solve_ok] = mybicgstab(coeff,rhs,L,[],[],printlevel);
          if (solve_ok<=0) & (printlevel)
             fprintf('\n  warning: bicgstab fails: %3.1f,',solve_ok);
          end
       end
    end
    if (printlevel>2); fprintf('%2.0d ',length(resnrm)-1); end
%%
    nnzmatold = nnzmat; matfct_options_old = matfct_options;
%%***************************************************************


text:
function munu_m_i(M,ML,MU,N)
%
%  Generates of the data used in 'munu_m'. Data are stored in global
%  variables.
%
%  ML, MU are generated by the routine 'munu_pre'
%
%  Calling sequence:
%
%    munu_m_i(M,ML,MU,N)
%
%  Input:
%
%    M         real matrix M
%    ML, MU    real LU factors,
%    N         real matrix N.
%
%  After calling this routine ML and MU can be deleted to save memory.
%
%
%  LYAPACK 1.0 (Thilo Penzl, September 1999)

if nargin~=4
  error('Wrong number of input arguments.');
end

global LP_M LP_ML LP_MU LP_N

LP_ML = ML;
LP_MU = MU;
LP_M = M;
LP_N = N;

author: Sable

repo: mcbench-benchmarks


text:
function [af, sf] = farras

% Farras nearly symmetric filters for orthogonal
% 2-channel perfect reconstruction filter bank
%
% USAGE:
%    [af, sf] = farras
% OUTPUT:
%    af - analysis filters
%    sf - synthesis filters
% REFERENCE:
%    A. F. Abdelnour and I. W. Selesnick.
%    "Nearly symmetric orthogonal wavelet bases",
%    Proc. IEEE Int. Conf. Acoust., Speech,
%    Signal Processing (ICASSP), May 2001.
% See afb, dwt.
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

af = [
                  0  -0.01122679215254
                  0   0.01122679215254
  -0.08838834764832   0.08838834764832
   0.08838834764832   0.08838834764832
   0.69587998903400  -0.69587998903400
   0.69587998903400   0.69587998903400
   0.08838834764832  -0.08838834764832
  -0.08838834764832  -0.08838834764832
   0.01122679215254                  0
   0.01122679215254                  0
];

sf = af(end:-1:1, :);

author: thomaskuestner


text:
        % Deschiderea unei ferestre de figuri
figure(1)
t = -10*pi:pi/250:10*pi;
        % Desenarea cometei tridimensionale
comet3((cos(2*t).^2).*sin(t),(sin(2*t).^2).*cos(t),t);


text:
function e = mohsst5_rmsew(Y, Yh)

Z = Y-Yh;
I = ~isnan(Z);

% Weight proportional to area size
load('mohsst5_data.mat', 'latitude', 'longitude');
%[lat,lon] = meshgrid(latitude, longitude);
[lon,lat] = meshgrid(longitude, latitude);
W = repmat(cosd(lat(:)), [1, size(Y,2)]);

e = rmsew(Z(I), W(I));

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function s = paramvalue2struct(varargin)

for i = 1:2:numel(varargin)

    s.(varargin{i}) = varargin{i+1};

end

end

author: pglira


text:
function [w,W]=morlwave(N,fo,fs,str)
%MORLWAVE  Morlet wavelet.
%
%   The Morlet wavelet is not recommended.  Instead use a generalized Morse
%   wavelet with GAMMA set to 3.  This function is included for testing and
%   comparison purposes.  See Lilly and Olhede (2009, 2012) for details.
%
%   PSI=MORLWAVE(N,FMAX,F) returns the complex-valued Morlet wavelet of
%   length N, with a peak in the frequency domain at *radian* F, based on a
%   mother wavelet having a peak frequency FMAX>1.
%
%   If F is a scalar, PSI is a column vector of length M.  If F is an
%   array, PSI is a matrix of size M x LENGTH(F), with the columns in order
%   of decreasing frequency.
%
%   The peak frequency FMAX should be greater than one.  As FMAX increases
%   from unity, the number of oscillations in the wavelet increases.
%
%   Note that the wavelets are centered at the midpoint in time, row
%   number ROUND(SIZE(PSI,1)/2).
%
%   [PSI,PSIF]=MORLWAVE(N,FMAX,F) also returns a size SIZE(PSI) matrix
%   PSIF which is the frequency-domain version of the wavelets.
%
%   For further details on this representation of the Morlet wavelet, see
%   Appendix A of
%
%       Lilly and Olhede (2009).  Higher-order properties of analytic
%           wavelets.  IEEE Trans. Sig. Proc., 57 (1), 146--160.
%   _________________________________________________________________
%
%   Normalization
%
%   MORLWAVE supports two kinds of normalization for the wavelets.
%
%   MORLWAVE(..., 'bandpass') uses "bandpass normalization", the
%   default.  See Lilly and Gascard (2006).  This implies that the
%   FFT of wavelet has a maximum value of 2.
%
%   MORLWAVE(..., 'energy') uses the unit energy normlization.  The
%   time-domain wavelet energy SUM(ABS(PSI).^2) is then unity.
%
%   Note MORLWAVE now uses bandpass-normalization by default.
%   _________________________________________________________________
%
%   See also MORLFREQ, MORLPROPS, MORSEWAVE.
%
%   Usage: psi=morlwave(n,fmax,f);
%          [psi,psif]=morlwave(n,fmax,f);
%
%   'morlwave --t' run a test.
%   'morlwave --f' generates a sample figure.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 F. Rekibi and J. M. Lilly
%                         --- type 'help jlab_license' for details

if strcmpi(N,'--t')
    morlwave_test;return
end

if strcmpi(N,'--f')
    type makefigs_morlwave
    makefigs_morlwave;
    return
end

if nargin<4
    str='bandpass';
end
%/********************************************************
%Enforce convention of high frequencies first
fs=fs(:);
if length(fs)>1
  if fs(2)-fs(1)>0
    fs=flipud(fs);
  end
end
%\********************************************************

for i=1:length(fs)
    w(:,i)= morlwave1(N,abs(fs(i)),abs(fo),str);
    if fs(i)<0
        w(:,i)=conj(w(:,i));
    end
end

wshifted=0*w;
for i=1:size(w,2)
  wshifted(:,i)=fftshift(w(:,i));
end
W=fft(wshifted);

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

function  [w]=morlwave1(N,fs,fo,str)

if fo <=1
    error('MORLWAVE frequency FMAX should be > 1')
end

s=frac(fo,fs);
a=morla(fo,str);

if ~isempty(strfind(str,'ene'))
    a=a.*frac(1,sqrt(s));
elseif ~isempty(strfind(str,'ban'))
    a=a.*frac(1,s);
end

t=(1:N)';   % Choose dt = 1
t=t-mean(t);
t=t./s;

nu=morlfreq(fo);
w=a.*exp(-frac(1,2).*t.^2).*(rot(nu.*t)-exp(-frac(1,2).*nu.^2));

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

function[a]=morla(fm,str)
%MORLA  Returns the Morlet wavelet amplitude coefficient "a".
%
%   A=MORLA(FMAX), given a peak frequency FMAX, returns the Morlet
%   wavelet amplitude called "A_{NU}" by Lilly and Olhede (2010).
%
%   By default, A is chosen such that the maximum of the frequency-
%   domain wavelet is equal to 2, the ``bandpass normalization.''
%
%   A=MORLA(FMAX,'energy') instead returns the coefficient giving
%   the wavelet unit energy.
%
%   Usage: a=morla(fmax);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008 J.M. Lilly --- type 'help jlab_license' for details

fnu=morlfreq(fm);

if strcmpi(str(1:3),'ban')
    a=2*frac(fm,fnu).*exp(frac(1,2).*(fm-fnu).^2)*frac(1,sqrt(2*pi));
elseif strcmpi(str(1:3),'ene')
    a=frac(1,sqrt(sqrt(pi))).*frac(1,sqrt(1+exp(-(fnu).^2)-2.*exp(-frac(3,4).*(fnu).^2)));
end

%\**************************************************************************

function []= morlwave_test

N=1028; % nombre de points
fs=logspace(log10(1/4),log10(10*1/4),10)';   %change to log space
[w,W]=morlwave(N,10/4,fs,'energy');
t=1:length(w);
t=t-mean(t);

tol=5e-3;

b=aresame(mean(w,1),0*w(1,:),tol);
reporttest('MORLWAVE zero mean',b)

b=aresame(sum(abs(w).^2,1),1+0*w(1,:),tol);
reporttest('MORLWAVE unit energy for energy normalization',b)

[w,W]=morlwave(N,10/4,fs(end),'bandpass');
b=aresame(maxmax(abs(W)),2,tol);
reporttest('MORLWAVE maximum of 2 for bandpass normalization',b)

fs=1./logspace(log10(5),log10(40))';
N=1023;
[psi,Psi]=morlwave(N,2*pi/4,fs,'bandpass');
reporttest('MORLWAVE Fourier transform, N odd, FS positive',aresame(abs(fft(psi,[],1)),abs(Psi),1e-8))
[psi,Psi]=morlwave(N+1,2*pi/4,fs,'bandpass');
reporttest('MORLWAVE Fourier transform, N even, FS positive',aresame(abs(fft(psi,[],1)),abs(Psi),1e-8))
[psi,Psi]=morlwave(N,2*pi/4,-fs,'bandpass');
reporttest('MORLWAVE Fourier transform, N odd, FS negative',aresame(abs(fft(psi,[],1)),abs(Psi),1e-8))
[psi,Psi]=morlwave(N+1,2*pi/4,-fs,'bandpass');
reporttest('MORLWAVE Fourier transform, N even, FS negative',aresame(abs(fft(psi,[],1)),abs(Psi),1e-8))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function distr = Create1DDistribution(values)
%
%
%        distr = Create1DDistribution(values)
%
%
%        Input:
%           -values: values
%
%        Output:
%           -distr: a 1D distribution
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

sumPDF = sum(values);

if(sumPDF <= 0)
    PDF = zeros(size(values));
    sumPDF = 1.0;
else
    PDF = values;
end

CDF = cumsum(PDF);
maxCDF = max(CDF);

if(maxCDF > 0)
    CDF = CDF / maxCDF;
end

distr = struct('PDF', PDF / sumPDF, 'CDF', CDF, 'maxCDF', maxCDF);

end

author: banterle


text:
function numgrad = computeNumericalGradient(J, theta)
% numgrad = computeNumericalGradient(J, theta)
% theta: a vector of parameters
% J: a function that outputs a real-number. Calling y = J(theta) will return the
% function value at theta.

% Initialize numgrad with zeros
numgrad = zeros(size(theta));

%% ---------- YOUR CODE HERE --------------------------------------
% Instructions:
% Implement numerical gradient checking, and return the result in numgrad.
% (See Section 2.3 of the lecture notes.)
% You should write code so that numgrad(i) is (the numerical approximation to) the
% partial derivative of J with respect to the i-th input argument, evaluated at theta.
% I.e., numgrad(i) should be the (approximately) the partial derivative of J with
% respect to theta(i).
%
% Hint: You will probably want to compute the elements of numgrad one at a time.

epsilon = 1e-4;

for i =1:length(numgrad)
    oldT = theta(i);
    theta(i)=oldT+epsilon;
    pos = J(theta);
    theta(i)=oldT-epsilon;
    neg = J(theta);
    numgrad(i) = (pos-neg)/(2*epsilon);
    theta(i)=oldT;
    if mod(i,100)==0
       fprintf('Done with %d\n',i);
    end;
end;

%% ---------------------------------------------------------------
end

author: xuzhenqi

repo: cnn

save_path: github-repos/MATLAB/xuzhenqi-cnn


text:
% Convert an image to polar coordinates
%
% This function converts an image to polar coordinates. The function is useful to represent
% circular objects in image as though they are rectangular. One example is a transformation
% of a firing rate map recorded in a circular environment. Such image contains a circle of firing
% rates. Polar image will have degrees along y-axis and distance from the middle point of arena
% towards it's end along y-axis. If you map is [20 20]. Middle point is [10 10]. A pixel with
% coordinates [10 15] will have polar coordinates [0 5], i.e. 0 degrees, 5 is the distance from [10 10].
% There is no exact transformation and some values are linearly interpolated.
% The function creates an image that have 360 values along y-axis and appr. half of the width/height of im
% along the x-axis. If im is [36 36], then polar image size is [18 360].
%
%  USAGE
%    pim = general.im2polar(im)
%    im     2D matrix with the image to convert, usually a square matrix of size [NxN].
%    pim    Polar version of im. 2D matrix of size [N/2 360].
%
function pim = im2polar(im)
    [r, c] = size(im);
    cx = c/2 + 0.5;
    cy = r/2 + 0.5;

    rmax = min([cx-1, c-cx, cy-1, r-cy]);
    [theta, radius] = meshgrid(1:360, 0:rmax);

    xi = radius .* cosd(theta) + cx;  % Locations in image to interpolate data
    yi = radius .* sind(theta) + cy;  % from.

    [x, y] = meshgrid(1:c, 1:r);
    pim = interp2(x, y, im, xi, yi);
end

author: kavli-ntnu


text:
function [D,X,Err] = simplestNMF(Y, opts)
    [n,T] = size(Y);

    if nargin<2
        opts = struct();
    end
    opts = initialize_opts(Y, opts);
    D = opts.D0;
    X = opts.X0;

    iter = 0;
    Err = zeros(opts.max_iter, 1);
    lambda_factor = repmat(opts.lambda.*ones(opts.m,1), 1, T);

    %% ----- Preprocess [start] -----
    nmflib_in_options.metric_type = opts.metric_type;
    nmflib_in_options.verbose = opts.verbose;
    nmflib_in_options.x_init.W = D;
    nmflib_in_options.x_init.H = X;
    nmflib_in_options.d_beta = opts.d_beta;
    [D, X, ~, nmflib_infos, nmflib_options] = pre_process('simplestNMF', Y, nmflib_in_options);
    % other necessary processes should be inserted here.
    %% ----- Preprocess [end] -----

    while iter < opts.max_iter
        iter = iter + 1;

        %update D
        if ~isempty(opts.updateD)
            %DX = D(:, opts.updateD)*X(opts.updateD, :);
            DX = D*X;
            if opts.d_beta<2
                DX(DX==0) = eps;
            end
            pgradD = (DX.^(opts.d_beta-1))*X(opts.updateD, :)';
            pgradD(pgradD==0)=eps;
            ngradD = ((DX.^(opts.d_beta-2)).*Y)*X(opts.updateD, :)';
            D(:, opts.updateD) = D(:, opts.updateD) .* ngradD ./ pgradD;
        end

        %update X
        DX = D*X;
        if opts.d_beta<2
            DX(DX==0) = eps;
        end
        pgradX = D'*(DX.^(opts.d_beta-1));
        ngradX = D'*((DX.^(opts.d_beta-2)).*Y);
        X = X .* ngradX ./ (pgradX + lambda_factor);

        DX = D*X;
        DX(DX==0) = eps;
        Err(iter) = beta_divergence(Y, DX, opts.d_beta);
        %sparsy = sum(sum(lambda_factor.*X));
        delta = inf;
        if iter>1
            delta = (Err(iter-1)-Err(iter))/Err(iter-1);
        end

        %% ----- Innerprocess [start] -----
        nmflib_infos = inner_process('simplestNMF', Y, D, X, [], nmflib_options, nmflib_infos, iter);
        %% ----- Innerprocess [end] -----

        %fprintf('iter = %d, Err = %f, delta = %f, sparsy = %f\n', iter, Err(iter), delta, sparsy);
        %if delta < opts.conv_value
        %    break;
        %end
    end
    Err(iter+1:end) = [];
end

function opts = initialize_opts(Y, opts)
    [n, T] = size(Y);
    if ~isfield(opts, 'm') opts.m = 1; end
    if ~isfield(opts, 'conv_value') opts.conv_value = 1e-3; end
    if ~isfield(opts, 'max_iter') opts.max_iter = 1000; end
    if ~isfield(opts, 'lambda') opts.lambda = eps; end
    if ~isfield(opts, 'beta') opts.d_beta = 2; end
    if ~isfield(opts, 'updateD') opts.updateD = 1:opts.m; end
    if ~isfield(opts, 'D0') opts.D0 = rand(n, opts.m); end
    if ~isfield(opts, 'X0') opts.X0 = rand(opts.m, T); end
end

function r = beta_divergence(A, B, beta)
    switch beta
        case 0
            r = sum(sum (A./B - log(A./B+eps) - 1));
        case 1
            r = sum(sum(A.*log(A./B+eps) - A + B));
        case 2
            r = .5 * norm(A - B, 'fro').^2;
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% std_checkset() - check STUDY set consistency
%
% Usage: >> [STUDY, ALLEEG] = std_checkset(STUDY, ALLEEG);
%
% Input:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of EEG datasets included in the STUDY structure
%
% Output:
%   STUDY      - a new STUDY set containing some or all of the datasets in ALLEEG,
%                plus additional information from the optional inputs above.
%   ALLEEG     - an EEGLAB vector of EEG sets included in the STUDY structure
%
% Authors:  Arnaud Delorme & Hilit Serby, SCCN, INC, UCSD, November 2005

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, arno@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

function [STUDY, ALLEEG, command] = std_checkset(STUDY, ALLEEG, option);

if nargin < 2
    help std_checkset;
    return;
end;
command  = '';

if isempty(STUDY), return; end;
studywasempty = 0;

modif = 0;
if ~isfield(STUDY, 'name'),  STUDY.name  = ''; modif = 1; end;
if ~isfield(STUDY, 'task'),  STUDY.task  = ''; modif = 1; end;
if ~isfield(STUDY, 'notes'), STUDY.notes = ''; modif = 1; end;
if ~isfield(STUDY, 'filename'),  STUDY.filename  = ''; modif = 1; end;
if ~isfield(STUDY, 'filepath'),  STUDY.filepath  = ''; modif = 1; end;
if ~isfield(STUDY, 'history'),   STUDY.history   = ''; modif = 1; end;
if ~isfield(STUDY, 'subject'),   STUDY.subject   = {}; modif = 1; end;
if ~isfield(STUDY, 'group'),     STUDY.group     = {}; modif = 1; end;
if ~isfield(STUDY, 'session'),   STUDY.session   = {}; modif = 1; end;
if ~isfield(STUDY, 'condition'), STUDY.condition = {}; modif = 1; end;
if ~isfield(STUDY, 'setind'),    STUDY.setind    = {}; modif = 1; end;
if ~isfield(STUDY, 'etc'),       STUDY.etc       = []; modif = 1; end;
if ~isfield(STUDY, 'etc.warnmemory'), STUDY.etc.warnmemory = 1; modif = 1; end;
if ~isfield(STUDY, 'preclust'),  STUDY.preclust  = []; modif = 1; end;
if ~isfield(STUDY, 'datasetinfo'), STUDY.datasetinfo = []; modif = 1; end;
if ~isfield(STUDY.etc, 'version'), STUDY.etc.version = []; modif = 1; end;
if ~isfield(STUDY.preclust, 'erpclusttimes' ),  STUDY.preclust.erpclusttimes = []; modif = 1; end;
if ~isfield(STUDY.preclust, 'specclustfreqs' ), STUDY.preclust.specclustfreqs = []; modif = 1; end;
if ~isfield(STUDY.preclust, 'erspclustfreqs' ), STUDY.preclust.erspclustfreqs = []; modif = 1; end;
if ~isfield(STUDY.preclust, 'erspclusttimes' ), STUDY.preclust.erspclusttimes = []; modif = 1; end;
if ~isfield(STUDY.datasetinfo, 'comps') & ~isempty(STUDY.datasetinfo), STUDY.datasetinfo(1).comps = []; modif = 1; end;
if ~isfield(STUDY.datasetinfo, 'index') & ~isempty(STUDY.datasetinfo), STUDY.datasetinfo(1).index = []; modif = 1; end;

% all summary fields
% ------------------
try, subject = unique_bc({ STUDY.datasetinfo.subject });
catch,
    subject = '';
    disp('Important warning: some datasets do not have subject codes; some functions may crash!');
end;
try, group = unique_bc({ STUDY.datasetinfo.group });
catch,
    group = {};
    % disp('Important warning: some datasets do not have group codes; some functions may crash!');
end;
try, condition = unique_bc({ STUDY.datasetinfo.condition });
catch,
    condition = {};
    disp('Important warning: some datasets do not have condition codes; some functions may crash!');
end;
try, session = unique_bc([STUDY.datasetinfo.session]);
catch,
    session = '';
    % disp('Important warning: some datasets do not have session numbers; some functions may crash!');
end;
if ~isequal(STUDY.subject,   subject  ), STUDY.subject   = subject;   modif = 1; end;
if ~isequal(STUDY.group,     group    ), STUDY.group     = group;     modif = 1; end;
if ~isequal(STUDY.condition, condition), STUDY.condition = condition; modif = 1; end;
if ~isequal(STUDY.session,   session  ), STUDY.session   = session;   modif = 1; end;

% check dataset info consistency
% ------------------------------
for k = 1:length(STUDY.datasetinfo)
    if ~strcmpi(STUDY.datasetinfo(k).filename, ALLEEG(k).filename)
        STUDY.datasetinfo(k).filename = ALLEEG(k).filename; modif = 1;
        fprintf('Warning: file name has changed for dataset %d and the study has been updated\n', k);
        fprintf('         to discard this change in the study, reload it from disk\n');
    end;
end;

% recompute setind array (setind is deprecated but we keep it anyway)
% -------------------------------------------------------------------
setind  = [];
sameica = std_findsameica(ALLEEG);
for index = 1:length(sameica)
    setind(length(sameica{index}):-1:1,index) = sameica{index}';
end;
setind(find(setind == 0)) = NaN;
if any(isnan(setind))
    warndlg('Warning: non-uniform set of dataset, some function might not work');
end
if ~isequal(setind, STUDY.setind), STUDY.setind = setind; modif = 1; end;

% check that dipfit is present in all datasets
% --------------------------------------------
for cc = 1:length(sameica)
    idat = [];
    for tmpi = 1:length(sameica{cc})
        if isfield(ALLEEG(sameica{cc}(tmpi)).dipfit, 'model')
            idat = sameica{cc}(tmpi);
        end;
    end;
    if ~isempty(idat)
        for tmpi = 1:length(sameica{cc})
            if ~isfield(ALLEEG(sameica{cc}(tmpi)).dipfit, 'model')
                ALLEEG(sameica{cc}(tmpi)).dipfit = ALLEEG(idat).dipfit;
                ALLEEG(sameica{cc}(tmpi)).saved  = 'no';
                fprintf('Warning: no ICA dipoles for dataset %d, using dipoles from dataset %d (same ICA)\n', sameica{cc}(tmpi), idat);
            end;
        end;
    end;
end;

% put in fake channels if channel labels are missing
% --------------------------------------------------
chanlabels = { ALLEEG.chanlocs };
if any(cellfun(@isempty, chanlabels))
    if any(~cellfun(@isempty, chanlabels))
        disp('********************************************************************');
        disp(' IMPORTANT WARNING: SOME DATASETS DO NOT HAVE CHANNEL LABELS AND ');
        disp(' SOME OTHERs HAVE CHANNEL LABELS. GENERATING CHANNEL LABELS FOR ');
        disp(' THE FORMER DATASETS (THIS SHOULD PROBABLY BE FIXED BY THE USER).');
        disp('********************************************************************');
    end;
    disp('Generating channel labels for all datasets...');
    for currentind = 1:length(ALLEEG)
        for ind = 1:ALLEEG(currentind).nbchan
            ALLEEG(currentind).chanlocs(ind).labels = int2str(ind);
        end;
    end;
    ALLEEG(currentind).saved = 'no';
end;

if length( unique( [ ALLEEG.srate ] )) > 1
    disp('********************************************************************');
    disp(' IMPORTANT WARNING: SOME DATASETS DO NOT HAVE THE SAME SAMPLING ');
    disp(' RATE AND THIS WILL MAKE MOST OF THE STUDY FUNCTIONS CRASH. THIS');
    disp(' SHOULD PROBABLY BE FIXED BY THE USER.');
    disp('********************************************************************');
end;

% check cluster array
% -------------------
rebuild_design = 0;
if ~isfield(STUDY, 'cluster'), STUDY.cluster = []; modif = 1; end;
if ~isfield(STUDY, 'changrp'), STUDY.changrp = []; modif = 1; end;
if isempty(STUDY.changrp) && isempty(STUDY.cluster)
    rebuild_design = 1;
end;
if isfield(STUDY.cluster, 'sets'),
    if max(STUDY.cluster(1).sets(:)) > length(STUDY.datasetinfo)
        disp('Warning: Some datasets had been removed from the STUDY, clusters have been reinitialized');
        STUDY.cluster = [];
    end;
end;
if ~studywasempty
    if isempty(STUDY.cluster)
        modif = 1;
        [STUDY] = std_createclust(STUDY, ALLEEG, 'parentcluster', 'on');
    end;
    if length(STUDY.cluster(1).child) == length(STUDY.cluster)-1 && length(STUDY.cluster) > 1
        newchild = { STUDY.cluster(2:end).name };
        if ~isequal(STUDY.cluster(1).child, newchild)
            STUDY.cluster(1).child = newchild;
        end;
    end;
end;

% create STUDY design if it is not present
% ----------------------------------------
if ~studywasempty
    if isfield(STUDY.datasetinfo, 'trialinfo')
        alltrialinfo = { STUDY.datasetinfo.trialinfo };
        if any(cellfun(@isempty, alltrialinfo)) && any(~cellfun(@isempty, alltrialinfo))
            disp('Rebuilding trial information structure for STUDY');
            STUDY  = std_maketrialinfo(STUDY, ALLEEG); % some dataset do not have trialinfo and
            % some other have it, remake it for everybody
        end;
    end;
    if ~isfield(STUDY, 'design') || isempty(STUDY.design) || ~isfield(STUDY.design, 'name')
        STUDY  = std_maketrialinfo(STUDY, ALLEEG);
        STUDY  = std_makedesign(STUDY, ALLEEG);
        STUDY  = std_selectdesign(STUDY, ALLEEG,1);
        rebuild_design = 0;
    else
        if isfield(STUDY.design, 'indvar1')
            STUDY  = std_convertdesign(STUDY, ALLEEG);
        end;

        % convert combined independent variable values
        % between dash to cell array of strings
        % -------------------------------------
        for inddes = 1:length(STUDY.design)
            if length(STUDY.design(inddes).variable) == 0
                STUDY.design(inddes).variable(1).label = '';
                STUDY.design(inddes).variable(1).value = [];
            end;
            if length(STUDY.design(inddes).variable) == 1
                STUDY.design(inddes).variable(2).label = '';
                STUDY.design(inddes).variable(2).value = [];
            end;
            if ~isfield(STUDY.design(inddes).variable, 'pairing')
                STUDY.design(inddes).variable(1).pairing = 'on';
                STUDY.design(inddes).variable(2).pairing = 'on';
            end;
            for indvar = 1:length(STUDY.design(inddes).variable)
                for indval = 1:length(STUDY.design(inddes).variable(indvar).value)
                    STUDY.design(inddes).variable(indvar).value{indval} = convertindvarval(STUDY.design(inddes).variable(indvar).value{indval});
                end;
            end;
        end;

       if ~isfield(STUDY.design(1), 'cell') || isempty(STUDY.design(1).cell)
           fprintf('Warning: Importing STUDY from a newer version of EEGLAB - some information will be lost\n');
           STUDY = std_makedesign(STUDY, ALLEEG, 1, STUDY.design(1), 'defaultdesign', 'forceoff');
       end;
       for inddes = 1:length(STUDY.design)
             for indcell = 1:length(STUDY.design(inddes).cell)
                for indval = 1:length(STUDY.design(inddes).cell(indcell).value)
                    STUDY.design(inddes).cell(indcell).value{indval} = convertindvarval(STUDY.design(inddes).cell(indcell).value{indval});
                end;
            end;
            for indinclude = 1:length(STUDY.design(inddes).include)
                if iscell(STUDY.design(inddes).include{indinclude})
                    for indval = 1:length(STUDY.design(inddes).include{indinclude})
                        STUDY.design(inddes).include{indinclude}{indval} = convertindvarval(STUDY.design(inddes).include{indinclude}{indval});
                    end;
                end;
            end;

            % check for duplicate entries in filebase
            % ---------------------------------------
            if length( { STUDY.design(inddes).cell.filebase } ) > length(unique({ STUDY.design(inddes).cell.filebase }))
                if ~isempty(findstr('design_', STUDY.design(inddes).cell(1).filebase))
                    error('There is a problem with your STUDY, contact EEGLAB support');
                else
                    fprintf('Duplicate entry detected in Design %d, reinitializing design\n', inddes);
                    [STUDY com] = std_makedesign(STUDY, ALLEEG, inddes, STUDY.design(inddes), 'defaultdesign', 'forceoff');
                end
            end;
        end;
    end;

    if rebuild_design % in case datasets have been added or removed
        STUDY = std_rebuilddesign(STUDY, ALLEEG);
    end;

    % scan design to fix old paring format
    % ------------------------------------
    for design = 1:length(STUDY.design)
        for var = 1:length(STUDY.design(design).variable)
            if isstr(STUDY.design(design).variable(1).pairing)
                if strcmpi(STUDY.design(design).variable(1).pairing, 'paired')
                    STUDY.design(design).variable(1).pairing = 'on';
                elseif strcmpi(STUDY.design(design).variable(1).pairing, 'unpaired')
                    STUDY.design(design).variable(1).pairing = 'off';
                end;
            end;
            if isstr(STUDY.design(design).variable(2).pairing)
                if strcmpi(STUDY.design(design).variable(2).pairing, 'paired')
                    STUDY.design(design).variable(2).pairing = 'on';
                elseif strcmpi(STUDY.design(design).variable(2).pairing, 'unpaired')
                    STUDY.design(design).variable(2).pairing = 'off';
                end;
            end;
        end;
    end;

    % add filepath field if absent
    for ind = 1:length(STUDY.design)
        if ~isfield(STUDY.design, 'filepath') || (isnumeric(STUDY.design(ind).filepath) && isempty(STUDY.design(ind).filepath))
            STUDY.design(ind).filepath = '';
            STUDY.saved = 'no';
            modif = 1;
        end;
    end;

    % check that ICA is present and if it is update STUDY.datasetinfo
    allcompsSTUDY  = { STUDY.datasetinfo.comps };
    allcompsALLEEG = { ALLEEG.icaweights };
    if all(cellfun(@isempty, allcompsSTUDY)) && ~all(cellfun(@isempty, allcompsALLEEG))
        for index = 1:length(STUDY.datasetinfo)
            STUDY.datasetinfo(index).comps = [1:size(ALLEEG(index).icaweights,1)];
        end;
    end;

    % make channel groups
    % -------------------
    if ~isfield(STUDY, 'changrp') || isempty(STUDY.changrp)
        STUDY = std_changroup(STUDY, ALLEEG);
        modif = 1;
    end;
end;

% determine if there has been any change
% --------------------------------------
if modif;
    STUDY.saved = 'no';
    command = '[STUDY ALLEEG] = std_checkset(STUDY, ALLEEG);';
    addToHistory = true;
    % check duplicate
    if length(STUDY.history) >= length(command) && strcmpi(STUDY.history(end-length(command)+1:end), command)
        addToHistory = false;
    end;
    if addToHistory
        STUDY.history =  sprintf('%s\n%s',  STUDY.history, command);
    end;
end;

% convert combined independent variables
% --------------------------------------
function val = convertindvarval(val);
if isstr(val)
    inddash = findstr(' - ', val);
    if isempty(inddash), return; end;
    inddash = [ -2 inddash length(val)+1];
    for ind = 1:length(inddash)-1
        newval{ind} = val(inddash(ind)+3:inddash(ind+1)-1);
    end;
    val = newval;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [v, V_p, V_k, V_d] = pinHoleDepth(p,k,d)

% PINHOLEDEPTH Pin hole projection with distance measurement
%   [v, V_p, V_k, V_d] = PINHOLEDEPTH(l, k, d) projects the Eucliden point
%   l into a pinhole+depth camera. k is the intrinsic vector of the
%   subjacent pinhole model, and d the distortion vector.
%
%   See also PINHOLE, INVPINHOLEDEPTH.

%   Copyright 2015-     Joan Sola @ IRI-UPC-CSIC.

if nargout == 1
    [v(1:2,:),v(3,:)] = pinHole(p,k,d);
else
    [u,s,U_p,U_k,U_d] = pinHole(p,k,d);
    S_p = [0 0 1]; % s = p(3,1) --> ds/dp = [0 0 1];
    v = [u;s];
    V_p = [U_p ; S_p];
    V_k = [U_k ; zeros(1, length(k))];
    V_d = [U_d ; zeros(1, length(d))];
end

return

%%
syms x y z u0 v0 au av real
p = [x;y;z];
k = [u0;v0;au;av];
d = [];

[v, V_p, V_k, V_d] = pinHoleDepth(p,k,d);

simplify(V_p - jacobian(v,p))
simplify(V_k - jacobian(v,k))

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [ stats ] = ml_get_cv_grid_states_regression(test,train)
%ML_GET_CV_GRID_STATES_REGRESSION
%
%   input -----------------------------------------------------------------
%
%       o test  : cell array (M x 1), M is the number of parameters for
%                                     which K-fold CV has been evaluated on.
%
%           test{i}.mse  : (1 x K), K number of
%
%   output ----------------------------------------------------------------
%
%       o stats : struct , returns statistics on train and test sets
%
%           stats.train  : (statistics on the training data set)
%
%               - stats.train.acc : (statistics on the accuracy for the training set)
%
%                   - stats.train.mse.mean
%                   - stats.train.mse.std
%
%       the structure is the same for the test statistics
%
%

%% Check if grid search was done on > 1 parameters

N = size(test);

if length(N) > 1
    test = test(:);
    train = train(:);
end

%% Prepare train statistics struct
M                    = length(test);

%% Get statistics

for m=1:M

    %% Get Mean for [Train]
    stats.train.mse.mean(m)     = mean(train{m}.mse);
    stats.train.nmse.mean(m)    = mean(train{m}.nmse);
    stats.train.rmse.mean(m)    = mean(train{m}.rmse);
    stats.train.nrmse.mean(m)   = mean(train{m}.nrmse);
    stats.train.mae.mean(m)     = mean(train{m}.mae);
    stats.train.mare.mean(m)    = mean(train{m}.mare);
    stats.train.r.mean(m)       = mean(train{m}.r);
    stats.train.d.mean(m)       = mean(train{m}.d);
    stats.train.e.mean(m)       = mean(train{m}.e);
    stats.train.me.mean(m)      = mean(train{m}.me);
    stats.train.mre.mean(m)     = mean(train{m}.mre);

    %% Get Std for [Train]

    stats.train.mse.std(m)     = std(train{m}.mse);
    stats.train.nmse.std(m)    = std(train{m}.nmse);
    stats.train.rmse.std(m)    = std(train{m}.rmse);
    stats.train.nrmse.std(m)   = std(train{m}.nrmse);
    stats.train.mae.std(m)     = std(train{m}.mae);
    stats.train.mare.std(m)    = std(train{m}.mare);
    stats.train.r.std(m)       = std(train{m}.r);
    stats.train.d.std(m)       = std(train{m}.d);
    stats.train.e.std(m)       = std(train{m}.e);
    stats.train.me.std(m)      = std(train{m}.me);
    stats.train.mre.std(m)     = std(train{m}.mre);

    %% Get Mean for [Test]
    stats.test.mse.mean(m)     = mean(test{m}.mse);
    stats.test.nmse.mean(m)    = mean(test{m}.nmse);
    stats.test.rmse.mean(m)    = mean(test{m}.rmse);
    stats.test.nrmse.mean(m)   = mean(test{m}.nrmse);
    stats.test.mae.mean(m)     = mean(test{m}.mae);
    stats.test.mare.mean(m)    = mean(test{m}.mare);
    stats.test.r.mean(m)       = mean(test{m}.r);
    stats.test.d.mean(m)       = mean(test{m}.d);
    stats.test.e.mean(m)       = mean(test{m}.e);
    stats.test.me.mean(m)      = mean(test{m}.me);
    stats.test.mre.mean(m)     = mean(test{m}.mre);

    %% Get Std for [Test]

    stats.test.mse.std(m)     = std(test{m}.mse);
    stats.test.nmse.std(m)    = std(test{m}.nmse);
    stats.test.rmse.std(m)    = std(test{m}.rmse);
    stats.test.nrmse.std(m)   = std(test{m}.nrmse);
    stats.test.mae.std(m)     = std(test{m}.mae);
    stats.test.mare.std(m)    = std(test{m}.mare);
    stats.test.r.std(m)       = std(test{m}.r);
    stats.test.d.std(m)       = std(test{m}.d);
    stats.test.e.std(m)       = std(test{m}.e);
    stats.test.me.std(m)      = std(test{m}.me);
    stats.test.mre.std(m)     = std(test{m}.mre);

end

if length(N) > 1

    stats.train.mse.mean     = reshape(stats.train.mse.mean,N);
    stats.train.nmse.mean    = reshape(stats.train.nmse.mean,N);
    stats.train.rmse.mean    = reshape(stats.train.rmse.mean,N);
    stats.train.nrmse.mean   = reshape(stats.train.nrmse.mean,N);
    stats.train.mae.mean     = reshape(stats.train.mae.mean,N);
    stats.train.mare.mean    = reshape(stats.train.mare.mean,N);
    stats.train.r.mean       = reshape(stats.train.r.mean,N);
    stats.train.d.mean       = reshape(stats.train.d.mean,N);
    stats.train.e.mean       = reshape(stats.train.e.mean,N);
    stats.train.me.mean      = reshape(stats.train.me.mean,N);
    stats.train.mre.mean     = reshape(stats.train.mre.mean,N);

    stats.train.mse.std     = reshape(stats.train.mse.std,N);
    stats.train.nmse.std    = reshape(stats.train.nmse.std,N);
    stats.train.rmse.std    = reshape(stats.train.rmse.std,N);
    stats.train.nrmse.std   = reshape(stats.train.nrmse.std,N);
    stats.train.mae.std     = reshape(stats.train.mae.std,N);
    stats.train.mare.std    = reshape(stats.train.mare.std,N);
    stats.train.r.std       = reshape(stats.train.r.std,N);
    stats.train.d.std       = reshape(stats.train.d.std,N);
    stats.train.e.std       = reshape(stats.train.e.std,N);
    stats.train.me.std      = reshape(stats.train.me.std,N);
    stats.train.mre.std     = reshape(stats.train.mre.std,N);

    stats.test.mse.mean     = reshape(stats.test.mse.mean,N);
    stats.test.nmse.mean    = reshape(stats.test.nmse.mean,N);
    stats.test.rmse.mean    = reshape(stats.test.rmse.mean,N);
    stats.test.nrmse.mean   = reshape(stats.test.nrmse.mean,N);
    stats.test.mae.mean     = reshape(stats.test.mae.mean,N);
    stats.test.mare.mean    = reshape(stats.test.mare.mean,N);
    stats.test.r.mean       = reshape(stats.test.r.mean,N);
    stats.test.d.mean       = reshape(stats.test.d.mean,N);
    stats.test.e.mean       = reshape(stats.test.e.mean,N);
    stats.test.me.mean      = reshape(stats.test.me.mean,N);
    stats.test.mre.mean     = reshape(stats.test.mre.mean,N);

    stats.test.mse.std     = reshape(stats.test.mse.std,N);
    stats.test.nmse.std    = reshape(stats.test.nmse.std,N);
    stats.test.rmse.std    = reshape(stats.test.rmse.std,N);
    stats.test.nrmse.std   = reshape(stats.test.nrmse.std,N);
    stats.test.mae.std     = reshape(stats.test.mae.std,N);
    stats.test.mare.std    = reshape(stats.test.mare.std,N);
    stats.test.r.std       = reshape(stats.test.r.std,N);
    stats.test.d.std       = reshape(stats.test.d.std,N);
    stats.test.e.std       = reshape(stats.test.e.std,N);
    stats.test.me.std      = reshape(stats.test.me.std,N);
    stats.test.mre.std     = reshape(stats.test.mre.std,N);

end

end

author: epfl-lasa


text:
function mfunc = minmodB(v,M,h)

% function mfunc = minmodB(v,M,h)
% Purpose: Implement the TVB modified midmod function. v is a vector

mfunc = v(1,:);
ids = find(abs(mfunc) > M*h.^2);

if(size(ids,2)>0)
  mfunc(ids) = minmod(v(:,ids));
end
return

repo: Aero-matlab


text:
function [axisLabelCell,orientationStr,iop] = returnViewerAxisLabels(planC,scanNum)

indexS = planC{end};

if ~exist('scanNum','var') || isempty(scanNum)
    scanNum = 1;
end

try
    iop = planC{indexS.scan}(scanNum).scanInfo(1).imageOrientationPatient;

    if isempty(iop)
        iop = [1,0,0,0,1,0]'; % temporary. assign this during loading planC
    end

    [~,idx2] = max(abs(iop(1:3)));
    [~,idx1] = max(abs(iop(4:end)));

    iopProj = zeros(size(iop));
    iopProj(idx1 + 3) = iop(idx1 + 3);
    iopProj(idx2) = iop(idx2);

    iopProjSign = iopProj ./ abs(iop);

    iopProjSign(isnan(iopProjSign)) = 0;

    if isequal(abs(iopProjSign),[1; 0; 0; 0; 1; 0])
        sliceType = 'axial';
    elseif isequal(abs(iopProjSign),[1; 0; 0; 0; 0; 1])
        sliceType = 'coronal';
    else
        sliceType = 'sagittal';
    end

    switch sliceType

        case 'axial'
            % definitions
            HFS = [1; 0; 0; 0; 1; 0];
            HFP = [-1; 0; 0; 0; -1; 0];
            FFS = [-1; 0; 0; 0; 1; 0];
            FFP = [1; 0; 0; 0; -1; 0];

            if isequal(iopProjSign(:), HFS)
                axisLabelCell = {'A','P';'R','L';'S','I'};
                orientationStr = 'HFS';
            elseif isequal(iopProjSign(:), HFP)
                axisLabelCell = {'P','A';'L','R';'S','I'};
                orientationStr = 'HFP';
            elseif isequal(iopProjSign(:), FFS)
                axisLabelCell = {'A','P';'L','R';'S','I'};
                orientationStr = 'FFS';
            else %FFP
                axisLabelCell = {'P','A';'R','L';'S','I'};
                orientationStr = 'FFP';
            end
        case 'coronal'
            HFS = [1; 0; 0; 0; 0; 1];
            HFP = [-1; 0; 0; 0; 0; 1];
            FFS = [-1; 0; 0; 0; 0; -1];
            FFP = [1; 0; 0; 0; 0; -1];

            if isequal(iopProjSign(:), HFS)
                axisLabelCell = {'S','I';'R','L';'A','P'};
                orientationStr = 'HFS';
            elseif isequal(iopProjSign(:), HFP)
                axisLabelCell = {'S','I';'R','L';'A','P'};
                orientationStr = 'HFP';
            elseif isequal(iopProjSign(:), FFS)
                axisLabelCell = {'S','I';'R','L';'A','P'};
                orientationStr = 'FFS';
            else %FFP
                axisLabelCell = {'S','I';'R','L';'A','P'};
                orientationStr = 'FFP';
            end

        case 'sagittal' % eg [0; 1; 0; 0; 0; -1]
            HFS = [0; 1; 0; 0; 0; 1];
            HFP = [0; -1; 0; 0; 0; 1];
            FFS = [0; -1; 0; 0; 0; -1];
            FFP = [0; 1; 0; 0; 0; -1];

            if isequal(iopProjSign(:), HFS)
                axisLabelCell = {'A','P';'S','I';'R','L'};
                orientationStr = 'HFS';
            elseif isequal(iopProjSign(:), HFP)
                axisLabelCell = {'A','P';'I','S';'R','L'};
                orientationStr = 'HFP';
            elseif isequal(iopProjSign(:), FFS)
                axisLabelCell = {'P','A';'S','I';'R','L'};
                orientationStr = 'FFS';
            else %FFP
                axisLabelCell = {'P','A';'I','S';'R','L'};
                orientationStr = 'FFP';
            end
    end
catch err
    %default HFS
    disp(err);
    disp('Defaulting to HFS orientation');
    axisLabelCell = {'A','P';'R','L';'S','I'};
    orientationStr = 'HFS';
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef TriP < ALGORITHM
% <multi/many> <real/integer> <constrained>
% Tri-population based coevolutionary algorithm
% delta --- 0.9 --- The probability of choosing parents locally
% nr    ---   2 --- Maximum number of solutions replaced by each offspring

%------------------------------- Reference --------------------------------
% F. Ming, W. Gong, L. Wang, and C. Lu, A tri-population based
% co-evolutionary framework for constrained multi-objective optimization
% problems, Swarm and Evolutionary Computation, 2022, 70: 101055.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)

            %% Parameter setting
            [delta,nr] = Algorithm.ParameterSet(0.9,2);

            %% Generate the weight vectors
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            T = ceil(Problem.N/10);

            %% Detect the neighbours of each solution
            B = pdist2(W,W);
            [~,B] = sort(B,2);
            B = B(:,1:T);

            %% Generate random population
            Population  = Problem.Initialization();
            Z           = min(Population.objs,[],1);
            Population1 = Problem.Initialization();
            Population2 = Problem.Initialization();
            alpha_c     = 2./(1+exp(1).^(-Problem.FE*10/(3*Problem.maxFE)))-1;
            para        = ceil(Problem.maxFE/Problem.N)/2 - ceil(Problem.FE/(3*Problem.N));

            %% Evaluate the Population
            Tc               = 0.9 * ceil(Problem.maxFE/Problem.N);
            last_gen         = 20;
            change_threshold = 1e-1;
            search_stage     = 1; % 1 for push stage,otherwise,it is in pull stage.
            max_change       = 1;
            epsilon_k        = 0;
            epsilon_0        = 0;
            cp               = 2;
            alpha            = 0.95;
            tao              = 0.05;
            ideal_points     = zeros(ceil(Problem.maxFE/Problem.N),Problem.M);
            nadir_points     = zeros(ceil(Problem.maxFE/Problem.N),Problem.M);
            arch             = archive(Population,Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                gen        = ceil(Problem.FE/(2*Problem.N));
                pop_cons   = Population.cons;
                cv         = overall_cv(pop_cons);
                population = [Population.decs,Population.objs,cv];
                rf         = sum(cv <= 1e-6) / Problem.N;
                ideal_points(gen,:) = Z;
                nadir_points(gen,:) = max(population(:,Problem.D + 1 : Problem.D + Problem.M),[],1);

                % The maximumrate of change of ideal and nadir points rk is calculated.
                if gen >= last_gen
                    max_change = calc_maxchange(ideal_points,nadir_points,gen,last_gen);
                end

                % The value of e(k) and the search strategy are set.
                if gen < Tc
                    if max_change <= change_threshold && search_stage == 1
                        search_stage = -1;
                        epsilon_0 = max(population(:,end),[],1);
                        epsilon_k = epsilon_0;
                    end
                    if search_stage == -1
                        epsilon_k =  update_epsilon(tao,epsilon_k,epsilon_0,rf,alpha,gen,Tc,cp);
                    end
                else
                    epsilon_k = 0;
                end

                % For each solution
                for i = 1 : Problem.N
                    % Choose the parents
                    if rand < delta
                        P = B(i,randperm(size(B,2)));
                    else
                        P = randperm(Problem.N);
                    end

                    % Generate an offspring
                    Offspring = OperatorDE(Problem,Population(i),Population(P(1)),Population(P(2)));

                    % Update the ideal point
                    Z = min(Z,Offspring.obj);

                    g_old  = max(abs(Population(P).objs-repmat(Z,length(P),1)).*W(P,:),[],2);
                    g_new  = max(repmat(abs(Offspring.obj-Z),length(P),1).*W(P,:),[],2);
                    cv_old = overall_cv(Population(P).cons);
                    cv_new = overall_cv(Offspring.con) * ones(length(P),1);

                    if search_stage == 1 % Push Stage
                        Population(P(find(g_old>=g_new,nr))) = Offspring;
                    else  % Pull Stage  &&  An improved epsilon constraint-handling is employed to deal with constraints
                        Population(P(find(((g_old >= g_new) & (((cv_old <= epsilon_k) & (cv_new <= epsilon_k)) | (cv_old == cv_new)) | (cv_new < cv_old) ), nr))) = Offspring;
                    end
                end

                % Re-rank
                Population1 = Population1(randperm(Problem.N));
                Population2 = Population2(randperm(Problem.N));
                Lia   = ismember(Population2.objs,Population1.objs, 'rows');
                gamma = 1-sum(Lia)/Problem.N;
                [Population1,FrontNo1] = EnvironmentalSelection(Population1,Problem.N,alpha_c);
                [Population2,FrontNo2] = EnvironmentalSelection_noCon(Population2,Problem.N,alpha_c,gamma,para);

                % Offspring Reproduction
                Population_all   = [Population1,Population2];
                RankSolution_all = [FrontNo1,FrontNo2];
                MatingPool = TournamentSelection(2,2*Problem.N,RankSolution_all);
                Offspring  = OperatorGAhalf(Problem,Population_all(MatingPool));

                % Environmental Selection
                alpha_c = 2./(1+exp(1).^(-Problem.FE*5/Problem.maxFE)) - 1;
                para  = ceil(Problem.maxFE/Problem.N)/2 - ceil(Problem.FE/(2*Problem.N));
                [Population1,~] = EnvironmentalSelection([Population1,Offspring],Problem.N,alpha_c);
                [Population2,~] = EnvironmentalSelection_noCon([Population2,Offspring],Problem.N,alpha_c,gamma,para);

                % Output the non-dominated and feasible solutions.
                arch = archive([arch,Population],Problem.N);
                if Problem.FE >= Problem.maxFE
                    arch = archive([arch,Population1],Problem.N);
                    Population = arch;
                end
            end
        end
    end
end

function result = overall_cv(cv)
% The Overall Constraint Violation

    cv(cv <= 0) = 0;cv = abs(cv);
    result = sum(cv,2);
end

function max_change = calc_maxchange(ideal_points,nadir_points,gen,last_gen)
% Calculate the Maximum Rate of Change

    delta_value = 1e-6 * ones(1,size(ideal_points,2));
    rz = abs((ideal_points(gen,:) - ideal_points(gen - last_gen + 1,:)) ./ max(ideal_points(gen - last_gen + 1,:),delta_value));
    nrz = abs((nadir_points(gen,:) - nadir_points(gen - last_gen + 1,:)) ./ max(nadir_points(gen - last_gen + 1,:),delta_value));
    max_change = max([rz, nrz]);
end

function result = update_epsilon(tao,epsilon_k,epsilon_0,rf,alpha,gen,Tc,cp)
    if rf < alpha
        result = (1 - tao) * epsilon_k;
    else
        result = epsilon_0 * ((1 - (gen / Tc)) ^ cp);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function grid_index = sparse_grid_cc_index ( dim_num, level_max, point_num )

%*****************************************************************************80
%
%% SPARSE_GRID_CC_INDEX indexes the points forming a sparse grid.
%
%  Discussion:
%
%    The points forming the sparse grid are guaranteed to be a subset
%    of a certain product grid.  The product grid is formed by DIM_NUM
%    copies of a 1D rule of fixed order.  The orders of the 1D rule,
%    (called ORDER_1D) and the order of the product grid, (called ORDER)
%    are determined from the value LEVEL_MAX.
%
%    Thus, any point in the product grid can be identified by its grid index,
%    a set of DIM_NUM indices, each between 1 and ORDER_1D.
%
%    This routine creates the GRID_INDEX array, listing (uniquely) the
%    points of the sparse grid.
%
%    An assumption has been made that the 1D rule is closed (includes
%    the interval endpoints) and nested (points that are part of a rule
%    of a given level will be part of every rule of higher level).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer POINT_NUM, the total number of points in the grids.
%
%    Output, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of point indices,
%    representing a subset of the product grid of level LEVEL_MAX,
%    representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%

%
%  The outer loop generates LEVELs from 0 to LEVEL_MAX.
%
  point_num2 = 0;

  for level = 0 : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_closed ( dim_num, level_1d );
%
%  The product of the 1D orders gives us the number of points in this grid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  The inner (hidden) loop generates all points corresponding to given grid.
%
      grid_index2 = multigrid_index0 ( dim_num, order_1d, order_nd );
%
%  Adjust these grid indices to reflect LEVEL_MAX.
%
      grid_index2 = multigrid_scale_closed ( dim_num, order_nd, level_max, level_1d, ...
        grid_index2 );
%
%  Determine the first level of appearance of each of the points.
%
      grid_level = abscissa_level_closed_nd ( level_max, dim_num, order_nd, ....
        grid_index2 );
%
%  Only keep those points which first appear on this level.
%
      for point = 1 : order_nd

        if ( grid_level(point) == level )

          point_num2 = point_num2 + 1;

          grid_index(1:dim_num,point_num2) = grid_index2(1:dim_num,point);

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function binary_image = PTKSkeletonise(binary_image, fixed_points_global, reporting)
    % PTKSkeletonise. Performs a skeletonisation on a segmented airway tree.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if ~isa(binary_image, 'PTKImage')
        error('Requires a PTKImage as input');
    end

    if ~exist('reporting', 'var')
        reporting = CoreReportingDefault;
    end

    if isdeployed || exist('PTKFastIsSimplePoint') == 3 %#ok<EXIST>
        use_mex_simple_point = true;
    else
        use_mex_simple_point = false;
        warning_message = 'Could not find compiled mex file PTKFastIsSimplePoint. Using a slower version PTKIsSimplePoint instead. Increase program speed by running mex PTKFastIsSimplePoint.cpp in the mex folder.';
        reporting.ShowWarning('PTKSkeletonise:PTKFastIsSimplePointnotFound', warning_message, []);
    end

    fixed_points = binary_image.GlobalToLocalIndices(fixed_points_global);

    % Marks endpoints with 3
    binary_image = binary_image.Copy;
    binary_image.ChangeRawImage(MarkEndpoints(binary_image.RawImage, fixed_points));
    total_number_of_points = sum(binary_image.RawImage(:) > 0);

    binary_image.AddBorder(2);
    direction_vectors = CalculateDirectionVectors;

    raw_image = binary_image.RawImage;

    previous_image = zeros(size(raw_image), 'uint8');

    iteration = 0;

    while ~isequal(previous_image, raw_image)
        previous_image = raw_image;

        iteration = iteration + 1;
        if (iteration > 20)
            if isempty(reporting)
                error('Maximum number of iterations exceeded. This can occur if not all the airway endpoints have been specified correctly.');
            else
                reporting.Error('PTKSkeletonise:MaximumIterationsExceeded', 'Maximum number of iterations exceeded. This can occur if not all the airway endpoints have been specified correctly.');
            end
        end

        % For each of the 6 principal directions
        for direction = [5, 23, 11, 17, 13, 15]

            if ~isempty(reporting)
                number_remaining_points = sum(raw_image(:) > 0);
                progress_value = round(100*(1-number_remaining_points/total_number_of_points));
                reporting.UpdateProgressAndMessage(progress_value, ['Skeletonisation: Iteration ' int2str(iteration)]);
            end

            if ~isempty(reporting)
                if reporting.HasBeenCancelled
                    error('User cancelled');
                end
            end

            direction_vector = direction_vectors(direction,:);
            i = direction_vector(1);
            j = direction_vector(2);
            k = direction_vector(3);

            % Detect border points and get their indices
            [b_i, b_j, b_k] = ind2sub(size(raw_image) - [2 2 2], ...
                find((1 == raw_image(2:end-1, 2:end-1, 2:end-1)) & (0 == raw_image(2+i:end-1+i,2+j:end-1+j,2+k:end-1+k))));
            b_i = b_i + 1; b_j = b_j + 1; b_k = b_k + 1;

            % Iterate through each border point and delete (set to zero) if
            % it is a simple point
            for i = 1 : length(b_i)
                raw_image(b_i(i), b_j(i), b_k(i)) = ~IsPointSimple(raw_image, b_i(i), b_j(i), b_k(i), use_mex_simple_point);
            end
        end
    end
    binary_image.ChangeRawImage(raw_image);
    binary_image.RemoveBorder(2);
end

function is_simple = IsPointSimple(binary_image, i, j, k, use_mex_simple_point)
    if use_mex_simple_point
        % MEX function (fast)
        is_simple = PTKFastIsSimplePoint((binary_image(i-1:i+1, j-1:j+1, k-1:k+1)));
    else
        % Matlab function (slow)
        is_simple = PTKIsSimplePoint(binary_image(i-1:i+1, j-1:j+1, k-1:k+1));
    end
end

function binary_image = MarkEndpoints(binary_image, fixed_points)
    binary_image = int8(binary_image ~= 0);
    binary_image(fixed_points) = 3;
end

function direction_vectors = CalculateDirectionVectors
    [i, j, k] = ind2sub([3 3 3], 1 : 27);
    direction_vectors = [i' - 2, j' - 2, k' - 2];
end

author: tomdoel

repo: pulmonarytoolkit


text:
function gsp()
% GSP General Piecewise Spline Interpolation
%
% GSP is a graphical user interface which takes x and y points as inputs, along with
% the order m, then outputs mth order splines between each of the x points. In most
% cases additional constraints are to be entered as well.
%
% Example:
%
% Input:
%
% x points : [1 2 3 4]
% y points : [2 7 -1 5]
% Order    : 2
% *One Additional Constraint is Required*
% f''(1) = 0
%
% Output:
% Polynomials:
%   0             5            -3
% -13            57           -55
%  27          -183           305
%
% In the Range:
% 1  2
% 2  3
% 3  4
%
% This means the following:
%
% for 1 < x < 2   y = 0*x^2 + 5*x - 3
% for 2 < x < 3   y = -13*x^2 + 57*x - 55
% for 3 < x < 4   y = 27*x^2 - 183*x + 305
%
% You can obtain the mentioned forms using poly2sym (requires Symbolic Math Toolbox).
%
% If you want to delete constraints, put its values as '*'. For more help and examples,
% see the pictures accompanied in the zip file. An optional feature can be
% accesessed if the function MQUAKE is present. mquake.m can be found here:
% <http://www.mathworks.com/matlabcentral/fileexchange/22816>
%
% To launch the GUI, type gsp in the command window with this file in the current directory.
% Alternatively, you can choose Debug -> Run from this editor window, or press F5.
%
%
% Husam Aldahiyat, 2009
% numandina@gmail.com
%
global aA

%% figure and uicontrols
figure('units','normalized','position',[.2 .2 .65 .65],'menubar','none','numbertitle','off','color','w','name','General Splines')
axes('position',[.35 .1 .5 .5])
ed1=uicontrol('style','edit','units','normalized','position',[.025 .875 .1 .05],'backgroundcolor','w');
uicontrol('style','text','units','normalized','position',[.025 .94 .1 .025],'backgroundcolor','w','string','x Points');
ed2=uicontrol('style','edit','units','normalized','position',[.025 .75 .1 .05],'backgroundcolor','w');
uicontrol('style','text','units','normalized','position',[.025 .81 .1 .025],'backgroundcolor','w','string','Y Points');
px=uicontrol('style','listbox','units','normalized','position',[.2 .75 .1 .2],'max',2,'backgroundcolor','w',...
	'string',{'Example 1';'Example 2';'Example 3';'Example 4';'Example 5';'Example 6'},'callback',@pn);
ed3=uicontrol('style','edit','units','normalized','position',[.045 .63 .05 .05],'backgroundcolor','w','callback',@cond);
uicontrol('style','text','units','normalized','position',[.025 .69 .1 .025],'backgroundcolor','w','string','Spline Order');
err=uicontrol('style','text','foregroundcolor','r','units','normalized','position',[.025 .55 .22 .05],'backgroundcolor','w',...
	'horizontalalignment','left');
uicontrol('style','pushbutton','units','normalized','position',[.15 .63 .05 .05],'backgroundcolor','w','callback',@go,'string','Solve');
pv=uicontrol('style','pushbutton','units','normalized','position',[.22 .63 .075 .05],'backgroundcolor','w','callback',@dv,...
	'string','Continuity','visible','off');
thih=uicontrol('style','edit','units','normalized','position',[.025 .025 .125 .525],'backgroundcolor','w','max',2,'string','',...
	'horizontalalignment','right','fontsize',13,'fontname','calibri');
edc=uicontrol('style','edit','units','normalized','position',[.15 .025 .075 .525],'backgroundcolor','w','max',2,...
	'horizontalalignment','left','fontsize',13,'fontname','calibri','string','');
Ap=uicontrol('style','listbox','units','normalized','position',[.35 .675 .5 .25],'backgroundcolor','w','max',2,'fontname','courier',...
	'callback',@chv);
Ac=uicontrol('style','listbox','units','normalized','position',[.85 .675 .125 .25],'backgroundcolor','w','max',2,'fontname','courier');
uicontrol('style','text','units','normalized','position',[.35 .935 .5 .025],'backgroundcolor','w','max',2,'string','Spline Polynomials');
uicontrol('style','text','units','normalized','position',[.85 .935 .12 .025],'backgroundcolor','w','max',2,'string','For x in the range');

	function cond(varargin)
		m=str2double(get(ed3,'string'));
		if m>1
			x=str2num(get(ed1,'string')); %#ok
			if isempty(x)
				set(err,'string','Please give values for x')
				return
			end
			set(err,'string',sprintf('%d Additional Condition(s) Needed',m-1))
			trip=cell(m*2,1);
			for k1=1:m
				koi=repmat('''',1,k1);
				trip{k1}=['f ',koi,'(',num2str(x(1)),') = '];
				trip{k1+m}=['f ',koi,'(',num2str(x(end)),') = '];
			end
			set(thih,'string',trip);
		else
			set(err,'string','No Additional Constraints Required');
			set(edc,'string','')
			set(thih,'string','')
		end
		if m>1
			set(edc,'string',repmat('*',m*2,1))
		end
	end

	function pn(varargin)
		switch get(px,'value')
			case 1
				 set(ed1,'string','0 10')
				 set(ed2,'string','0 0')
				 set(ed3,'string','2')
				 cond
				 set(err,'string','')
				 set(edc,'string',['5';'*';'*';'*'])
			case 2
				 set(ed1,'string','0 1 2')
				 set(ed2,'string','0 1 0')
				 set(ed3,'string','3')
				 cond
				 set(err,'string','')
				 set(edc,'string',['0';'*';'*';'0';'*';'*'])
			case 3
				 set(ed1,'string','0 1 2')
				 set(ed2,'string','0 1 0')
				 set(ed3,'string','5')
				 cond
				 set(err,'string','')
				 set(edc,'string',['0';'0';'*';'*';'*';'0';'0';'*';'*';'*'])
			case 4
				 set(ed1,'string','1 2 3 4 7')
				 set(ed2,'string','2 -4 0 3 0')
				 set(ed3,'string','3')
				 cond
				 set(err,'string','')
				 set(edc,'string',['*';'0';'*';'*';'0';'*'])
			case 5
				 set(ed1,'string','1 2 4 6 8 10')
				 set(ed2,'string','2 4 -2 3 1 0')
				 set(ed3,'string','5')
				 cond
				 set(err,'string','')
				 set(edc,'string',strvcat('5','0','*','*','*','-4','3','*','*','*')) %#ok
			case 6
				 set(ed1,'string','1 2 4 6 8 10')
				 set(ed2,'string','2 4 -2 3 1 0')
				 set(ed3,'string','5')
				 cond
				 set(err,'string','')
				 set(edc,'string',strvcat('-5','0','*','*','*','4','3','*','*','*')) %#ok
		end
	end

%% Solution

	function go(varargin)
		set(pv,'visible','off')
		x=str2num(get(ed1,'string')); %#ok
		y=str2num(get(ed2,'string')); %#ok
		if length(y)~=length(x)
			set(err,'string','Points are not of equal lengths')
			return
		end
		if any(find(abs(sort(x)-x)>1e-10))
			set(err,'string','x Points need to be ascending')
			return
		end
		if length(unique(x))<length(x)
			set(err,'string','Repeated x points found')
			return
		end
		m=str2double(get(ed3,'string'));
		n=length(x)-1; % number of splines
		g=fliplr(0:m);
		DAT=zeros((m+1)*n);
		B=zeros(n*(1+m),1);
		B(1:n+1)=y;
		S=zeros(n,m+1);
		for k=1:n
			DAT(k,k+m*k-m:k+m*k)=x(k).^g;
		end
		DAT(k+1,k+m*k-m:k+m*k)=x(k+1).^g;
		k=k+2;
		poln=ones(1,m+1);
		for k2=1:m
			pol=[poln,zeros(1,(k2-1))];
			for p=1:n-1
				DAT(k+p+k2*n-k2-n,p*m+p-m:p*m+p+m+1)=...
					repmat(pol,1,2).*repmat(x(p+1),1,(m+1)*2).^repmat(g-(k2-1),1,2).*[ones(1,m+1),-ones(1,m+1)];
			end
			poln=polyder(poln);
		end
		aa=(get(edc,'string'));
		if m>1
			con=find(~ismember(1:m*2,findstr(aa(:,1)','*'))==1);
			if length(con)>m-1
				set(err,'string','Too many constraints entered')
				return
			end
			if isempty(con)
				set(err,'string','Too few constraints entered')
				return
			end
			for lp=1:m-1
				poln=ones(1,m+1);
				wo=con(1);
				if wo>m
					con(1)=con(1)-m;
				end
				for h1=1:con(1)
					poln=polyder(poln);
				end
				if wo>m
					DAT(end-(m-1-lp),end-m:end)=repmat(x(end),1,m+1).^(g-con(1)).*[poln,zeros(1,m+1-length(poln))];
					B(end-(m-1-lp))=str2double(aa(wo,:));
				end
				if wo<=m
					DAT(end-(m-1-lp),1:m+1)=[poln,zeros(1,m+1-length(poln))].*repmat(x(1),1,m+1).^(g-con(1));
					B(end-(m-1-lp))=str2double(aa(con(1),:));
				end
				con(1)=[];
			end
		end
		DAT(isnan(DAT))=0;
		s=DAT\B;
		C=zeros(n,2);
		for k=1:n
			S(k,:)=s(k*(m+1)-m:k*(m+1));
			C(k,:)=x(k:k+1);
		end
		set(err,'string','')

%% Displaying Results

		set(Ap,'string',num2str(S))
		set(Ac,'string',num2str(C))
		cla reset
		hold on
		st=[1 0 0;0 1 0];
		for k=1:n
			t=linspace(C(k,1),C(k,2),100);
			f=polyval(S(k,:),t);
			aA(k)=plot(t,f);
			set(findobj('color','b'),'color',st(mod(k,2)+1,:),'linewidth',2)
		end
		plot(x,y,'.')
		set(findobj('marker','.'),'markersize',10,'color','k')
		xlabel('x')
		ylabel('y')
		xlim([x(1)-1 x(end)+1])
		h=axis;
		if h(3)>min(y)-1
			h(3)=min(y)-1;
		end
		if h(4)<max(y)+1
			h(4)=max(y)+1;
		end
		axis(h)
		set(pv,'visible','on')
	end

	function dv(varargin)
		S=str2num(get(Ap,'string')); %#ok
		C=str2num(get(Ac,'string')); %#ok
		st=[1 0 0;0 1 0];
		m=str2double(get(ed3,'string'));
		n=length(str2num(get(ed1,'string')))-1; %#ok
		figure('color','w','menubar','none','numbertitle','off','name','Continuity')
		for k2=1:m
			for g=1:n
				S(g,:)=[zeros(m+1-length([polyder(S(g,1:end-k2+1)),zeros(1,k2)])),polyder(S(g,1:end-k2+1)),zeros(1,k2)];
			end
			subplot(ceil(m/2),ceil(m/ceil(m/2)),k2)
			hold on
			for k=1:n
				t=linspace(C(k,1),C(k,2),100);
				f=polyval(S(k,1:end-k2),t);
				plot(t,f);
				set(findobj('color','b'),'color',st(mod(k,2)+1,:),'linewidth',2)
			end
			title(sprintf('Derivate Order %d',k2))
		end
	end

	function chv(varargin)
		h=axis;
		try
			mquake(aA(get(Ap,'value')),(h(2)+h(4))/4,.5,1)
		catch
		end
	end

end


text:
function y = subsasgn(X,I,Y)
%SUBASGN (overloaded)

try
    if strcmp('()',I.type)
        X_is_spdvar = isa(X,'sdpvar') |  isa(X,'ndsdpvar');
        Y_is_spdvar = isa(Y,'sdpvar') |  isa(Y,'ndsdpvar');
        if islogical(I.subs{1})
            I.subs{1} = double(find(I.subs{1}));
        end
        if any(I.subs{1} <=0)
            error('Index into matrix is negative or zero.');
        end

        if length(I.subs)>2 && X_is_spdvar
            y = subsasgn(ndsdpvar(X),I,Y);
            return
        end

        switch 2*X_is_spdvar+Y_is_spdvar
            case 1
                % This code does not work properly
                % Only work if b is undefined!!?!!
                % generally ugly code...
                y = Y;
                [n_y,m_y] = size(Y);
                y_lmi_variables = y.lmi_variables;
                try
                    X0 = subsasgn(full(X),I,full(reshape(Y.basis(:,1),n_y,m_y)));
                    dim = size(X0);
                    y.basis = reshape(X0,prod(dim),1);
                    X = full(X)*0;
                    for i = 1:length(y_lmi_variables)
                        X0 = subsasgn(X,I,full(reshape(Y.basis(:,i+1),n_y,m_y)));
                        y.basis(:,i+1) = reshape(X0,prod(dim),1);
                    end
                    y.dim = dim;
                    % Reset info about conic terms
                    y.conicinfo = [0 0];
                    y.basis = sparse(y.basis);
                    if length(dim)>2
                        y = ndsdpvar(y);
                    end
                catch
                    error(lasterr)
                end
            case 2
                if ~isempty(Y)
                    if isa(Y,'uint8') || isa(Y,'uint16') || isa(Y,'uint32') || isa(Y,'uint64')
                        Y = sparse(double(Y));
                    elseif isnumeric(Y)
                        Y = sparse(Y);
                    else
                        Y = sparse(double(Y));
                    end
                end
                y = X;

                % Special code for speed
                % elements in vector replaced with constants
                if min(X.dim(1),X.dim(2))==1 & (length(I.subs)==1)
                     y = X;
                     if isempty(Y)
                         y.basis(I.subs{1},:) = [];
                         if X.dim(1) == 1
                             y.dim(2) = y.dim(2) - length(unique(I.subs{1}));
                         else
                             y.dim(1) = y.dim(1) - length(unique(I.subs{1}));
                         end
                     else
                         y.basis(I.subs{1},1) = Y;
                         y.basis(I.subs{1},2:end) = 0;
                     end
                     if prod(y.dim)~=size(y.basis,1)
                         % Ah bugger, the dimension of the object was changed)
                         aux = X.basis(:,1);
                         aux = reshape(aux,X.dim);
                         aux(I.subs{1})=Y;
                         y.dim = size(aux);
                     end
                     y = clean(y);
                     % Reset info about conic terms
                     if isa(y,'sdpvar')
                         y.conicinfo = [0 0];
                     end
                     return;
                end

                x_lmi_variables = X.lmi_variables;
                lmi_variables = [];

                n = y.dim(1);
                m = y.dim(2);
                subX = sparse(subsasgn(full(reshape(X.basis(:,1),n,m)),I,Y));
                y.basis = subX(:);
                if isa(I.subs{1},'char')
                    I.subs{1} = 1:n;
                end
                if length(I.subs)>1
                    if isa(I.subs{2},'char')
                        I.subs{2} = 1:m;
                    end
                end
                if length(I.subs)>1
                    if length(I.subs{1})==1 & length(I.subs{2})~=1
                        I.subs{1} = repmat(I.subs{1},size(I.subs{2},1),size(I.subs{2},2));
                    elseif length(I.subs{2})==1 & length(I.subs{1})~=1
                        I.subs{2} = repmat(I.subs{2},size(I.subs{1},1),size(I.subs{1},2));
                    end
                end

                if length(I.subs)>1
                    ii = kron(I.subs{1}(:),ones(length(I.subs{2}),1));
                    jj = kron(ones(length(I.subs{1}),1),I.subs{2}(:));
                    LinearIndex = sub2ind([n m],ii,jj);
                else
                    LinearIndex = I.subs{1};
                end

                if isempty(Y)
                    X.basis = X.basis(:,2:end);
                    X.basis(LinearIndex,:) = [];
                    y.basis = [y.basis(:,1) X.basis];
                else
                    X.basis(LinearIndex,2:end)=sparse(0);
                    y.basis = [y.basis(:,1) X.basis(:,2:end)];
                end

                y.dim(1) = size(subX,1);
                y.dim(2) = size(subX,2);
                if ~isempty(x_lmi_variables)
                    y = clean(y);
                end
                if isa(y,'sdpvar')
                    % Reset info about conic terms
                    y.conicinfo = [0 0];
                end

            case 3
                z = X;

                x_lmi_variables = X.lmi_variables;
                y_lmi_variables = Y.lmi_variables;

                % In a first run, we fix the constant term and null terms in the X basis
                lmi_variables = [];
                nx = X.dim(1);
                mx = X.dim(2);
                ny = Y.dim(1);
                my = Y.dim(2);

                if (mx==1) & (my == 1) & isempty(setdiff(y_lmi_variables,x_lmi_variables)) & (max(I.subs{1}) < nx) & length(I.subs)==1 & length(unique(I.subs{1}))==length(I.subs{1}) ;
                    % Fast specialized code for Didier
                     y = specialcode(X,Y,I);
                     return
                end

                subX = subsasgn(reshape(X.basis(:,1),nx,mx),I,reshape(Y.basis(:,1),ny,my));
                [newnx, newmx] = size(subX);

                j = 1;

                yz = reshape(1:ny*my,ny,my);
                subX2 = subsasgn(reshape(zeros(nx*mx,1),nx,mx),I,yz);
                subX2 = subX2(:);
                [ix,jx,sx] = find(subX2);
                yz = 0*reshape(Y.basis(:,1),ny,my);
                lmi_variables = zeros(1,length(x_lmi_variables));

                A = reshape(1:nx*mx,nx,mx);
                B = reshape(1:newnx*newmx,newnx,newmx);

                rm = B(1:nx,1:mx);rm = rm(:);
                [iix,jjx,ssx] = find(X.basis(:,2:end));
                z.basis = [subX(:) sparse(rm(iix),jjx,ssx,newnx*newmx,size(X.basis,2)-1)];
                z.basis(ix,2:end) = 0;

                keep = find(any(z.basis(:,2:end),1));
                z.basis = z.basis(:,[1 1+keep]);
                lmi_variables2 = x_lmi_variables(keep);

                z.lmi_variables = lmi_variables2;
                lmi_variables = lmi_variables2;

                all_lmi_variables = union(lmi_variables,y_lmi_variables);
                in_z = ismembcYALMIP(all_lmi_variables,lmi_variables);
                in_y = ismembcYALMIP(all_lmi_variables,y_lmi_variables);
                z_ind = 2;
                y_ind = 2;
                basis = spalloc(size(z.basis,1),1+length(all_lmi_variables),0);
                basis(:,1) = z.basis(:,1);
                nz = size(subX,1);
                mz = size(subX,2);
                template = full(0*reshape(X.basis(:,1),nx,mx));
                in_yin_z = 2*in_y + in_z;
                if all(in_yin_z<3)
                    case1 = find(in_yin_z==1);
                    if ~isempty(case1)
                        basis(:,case1+1) = z.basis(:,2:1+length(case1));
                        in_yin_z(case1) = 0;
                    end
                end
                    % Let's identify the indices at which we need to intervene
                    case1 = find(in_yin_z==1);
                    checkI = union(find(in_yin_z >= 2), setdiff(case1,case1+1));
                    if size(checkI,1) > size(checkI,2)
                        % Sometimes the in_yin_z vector is vertical (not
                        % always though), so we make sure that checkI, on
                        % which we'll iterate, is horizontal.
                        checkI = checkI';
                    end
                    for i = checkI
                        switch in_yin_z(i)
                            case 1
                                % We look for the end of the block of ones
                                % starting at i
                                iend = i + find([in_yin_z(i+1:end) 0] ~= 1, 1, 'first')-1;
                                basis(:,i+1:iend+1) = z.basis(:,z_ind:z_ind+(iend-i));z_ind = z_ind+(iend-i)+1;
                        case 2
                            temp = sparse(subsasgn(template,I,full(reshape(Y.basis(:,y_ind),ny,my))));
                            basis(:,i+1) = temp(:);
                            y_ind = y_ind+1;
                        case 3
                            Z1 = z.basis(:,z_ind);
                            Z4 = Y.basis(:,y_ind);
                            Z3 = reshape(Z4,ny,my);
                            Z2 = sparse(subsasgn(0*reshape(full(X.basis(:,1)),nx,mx),I,Z3));
                            temp = reshape(Z1,nz,mz)+Z2;
                            basis(:,i+1) = temp(:);
                            z_ind = z_ind+1;
                            y_ind = y_ind+1;
                        otherwise
                    end
                end;
                z.dim(1) = nz;
                z.dim(2) = mz;
                z.basis = basis;
                z.lmi_variables = all_lmi_variables(:)';
                y = z;
                % Reset info about conic terms
                y.conicinfo = [0 0];
            otherwise
        end
    else
        error('Reference type not supported');
    end

catch
    error(lasterr)
end

function y = specialcode(X,Y,I)

y = X;
X_basis = X.basis;
Y_basis = Y.basis;
ind = I.subs{1};ind = ind(:);
yvar_in_xvar = zeros(length(Y.lmi_variables),1);
for i = 1:length(Y.lmi_variables);
    yvar_in_xvar(i) = find(X.lmi_variables==Y.lmi_variables(i));
end
y.basis(ind,:) = 0;
mapper = [1 1+yvar_in_xvar(:)'];mapper = mapper(:);
[i,j,k] = find(y.basis);
[ib,jb,kb] = find(Y_basis);
i = [i(:);ind(ib(:))];
j = [j(:);mapper(jb(:))];
k = [k(:);kb(:)];
y.basis = sparse(i,j,k,size(y.basis,1),size(y.basis,2));
y = clean(y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% interp_table_test.m

hr = inline('(1-abs(t/(J/2))) .* (abs(t) <= J/2)', 't', 'J');
hi = inline('(1-abs(t/(J/2))).^2 .* (abs(t) <= J/2)', 't', 'J');

% 1D
if 0
	L = 10;
	J = 6;
	s = [-J/2*L:J/2*L]'/L;
	hs = hr(s, J) + 1i * hi(s, J);
	if length(hs) ~= J*L+1, error 'size', end

	clf, subplot(211)
	plot(s, real(hs), 'c.-', s, imag(hs), 'y.-')

	K = 20;
	ck = zeros(K,1);
	ck(2+1) = 1;
	%ck = complexify(ck);
	tm = linspace(-2*K, 2*K, 2001)';
	fm = interp1_table_mex(ck, hs, int32(J), int32(L), tm);

	subplot(212)
	plot(tm, real(fm), 'c.-', tm, imag(fm), 'y.-')
end

% 2D
if 0
	L = [2^5 2^4];
	J = [6 4];
	s1 = [-J(1)/2*L(1):J(1)/2*L(1)]'/L(1);
	s2 = [-J(2)/2*L(2):J(2)/2*L(2)]'/L(2);
	h1 = 0*hr(s1, J(1)) + 1i * hi(s1, J(1));
	h2 = 0*hr(s2, J(2)) + 1i * hi(s2, J(2));
%	h1 = complexify(h1);
%	h2 = complexify(h2);

	clf, subplot(211)
	plot(	s1, real(h1), 'c.-', s1, imag(h1), 'y.-', ...
		s2, real(h2), 'g.-', s2, imag(h2), 'm.-')

	K = [8 12];
	ck = zeros(K);
	ck(0+1, 0+1) = -1i;
%	ck = complexify(ck);
	t1 = linspace(-2*K(1), 2*K(1), 201)';
	t2 = linspace(-2*K(2), 2*K(2), 199)';
	t1 = linspace(0, K(1), 201)';
	t2 = linspace(0, K(2), 199)';
	[tt1 tt2] = ndgrid(t1, t2);
	tm = [tt1(:) tt2(:)];
%	tic
	fm = interp2_table_mex(ck, h1, h2, int32(J), int32(L), tm);
%	toc
	fm = reshape(fm, size(tt1));

	im(121, t1, t2, real(fm), 'real'), cbar
	im(122, t1, t2, imag(fm), 'imag'), cbar
end

% 3D
if 1
	L = [2^5 2^4 2^6];
	J = [5 3 4];
	s1 = [-J(1)/2*L(1):J(1)/2*L(1)]'/L(1);
	s2 = [-J(2)/2*L(2):J(2)/2*L(2)]'/L(2);
	s3 = [-J(3)/2*L(3):J(3)/2*L(3)]'/L(3);
	h1 = 1 * hr(s1, J(1)) + 0i * hi(s1, J(1));
	h2 = 1 * hr(s2, J(2)) + 0i * hi(s2, J(2));
	h3 = 1 * hr(s3, J(3)) + 0i * hi(s3, J(3));
	h1 = complexify(h1);
	h2 = complexify(h2);
	h3 = complexify(h3);

	clf, subplot(211)
	plot(	s1, real(h1), 'c.-', s1, imag(h1), 'y.-', ...
		s2, real(h2), 'g.-', s2, imag(h2), 'm.-', ...
		s3, real(h3), 'y.-', s3, imag(h3), 'w.-')

	K = [8 12 4];
	ck = zeros(K);
	ck(0+1, 0+1) = 1i;
	ck = complexify(ck);
	t1 = linspace(-2*K(1), 2*K(1), 69)';
	t2 = linspace(-2*K(2), 2*K(2), 89)';
	t3 = linspace(-2*K(3), 2*K(3), 9)';
%	t1 = linspace(0, K(1), 201)';
%	t2 = linspace(0, K(2), 199)';
	t3 = [-1 0 1];
	[tt1 tt2 tt3] = ndgrid(t1, t2, t3);
	tm = [tt1(:) tt2(:) tt3(:)];
	tic
	fm = interp3_table_mex(ck, h1, h2, h3, int32(J), int32(L), tm);
	toc
	fm = reshape(fm, size(tt1));

	im(121, t1, t2, real(fm), 'real'), cbar
	im(122, t1, t2, imag(fm), 'imag'), cbar
	if length(t3) == 1 % compare to 2D
		f2 = interp2_table_mex(ck, h1, h2, ...
			int32(J(1:2)), int32(L(1:2)), tm(:,1:2));
		f2 = reshape(f2, size(tt1));
		max_percent_diff(f2, fm)
	end
end


text:
function R = gamrand1(A, B)
%GAMRAND1 Random matrices from gamma distribution.
%
%   R = GAMRAND1(A,B) returns a matrix of random numbers chosen
%   from the gamma distribution with parameters A and B.
%   Both parameters have to be scalar.
%
%   Note: Parameterization as in (Neal, 1996).
%      A is mean of the distribution
%      B is degrees of freedom
%
%	See also INVGAMRAND
%
% Copyright (c) 1999 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

error('No mex-file for this architecture. See Matlab help and convert.m in ./linuxCsource or ./winCsource for help.')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function fx2 = p15_fx2 ( x )

%*****************************************************************************80
%
%% P15_FX2 evaluates the second derivative of the function for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  fx2 = ( - 2.0 - 30.0 * x.^2 + 36.0 * x.^4 ) * r8_cube_root ( x ) ...
    * exp ( - x.^2 ) / ( 9.0 * x.^2 );

  return
end


text:
% This function compute the instantaneous frequency from phase spectrogram
% Input:
%   phase: D x T matrix of phase spectrogram, where D is FFT-bin number and
%       T is number of frames.
% Output:
%   instantaneous frequency of the size size as phase.
% The implementation is based on equation (3) of
%   Krawczyk, Martin, and Timo Gerkmann. "STFT Phase Reconstruction in
%   Voiced Speech for an Improved Single-Channel Speech Enhancement." (2014).
%
% Author: Xiong Xiao, Nanyang Tech. Univ., Singapore.
% Created: 06 Jan 2015
% Last modified: 06 Jan 2015
%
function [instan_freq] = comp_instan_freq(phase)

phase_delta = phase(:,2:end) - phase(:,1:end-1);
phase_delta = [phase(:,1) phase_delta];

[dim, nFr] = size(phase_delta);

phase_delta_vec = reshape(phase_delta, dim*nFr, 1);
phase_delta_vec = get_principal_value(phase_delta_vec);
instan_freq = reshape(phase_delta_vec, dim, nFr);

if 0
    subplot(3,1,1); imagesc(phase); title('Phase'); colorbar;
    subplot(3,1,2); imagesc(phase_delta); title('Phase delta');colorbar;
    subplot(3,1,3); imagesc(instan_freq); title('Phase delta principal value - Instantaneous frequency');colorbar;
end

end

function phase_delta_vec = get_principal_value(phase_delta_vec)
idx = find(phase_delta_vec>pi);
phase_delta_vec(idx) = phase_delta_vec(idx) - 2*pi;

idx = find(phase_delta_vec<-pi);
phase_delta_vec(idx) = phase_delta_vec(idx) + 2*pi;

end

author: singaxiong

repo: SignalGraph


text:
function Lms = gallager_ach(Ns, delta, epsil)
% computes gallager's achievability bound for BSC
%
% We use it in the form
%
%     P_e \le M^s  (   Sum_y [ Sum_x Q(x) W(y|x)^{1/1+s} ]^(1+s)   )^n
%
% Note: we then conclude that from (M, P_e) code with AVG P_e
% we can get (M/2, 2 P_e) code with MAXIMAL. Thus we must take P_e = epsil/2 and then
% subtract 1 bit from gallager's result
%

Lms = [];

for n = Ns;

	f = @(x) -( log2(epsil/2) ./ x + n - n .* (1+x) ./ x .* log2( delta.^(1./(1+x)) + (1-delta).^(1./(1+x))));

	[lambda f] = fminbnd(f, 0, 1);

	% Cut half of the codewords

	logm = -f-1;

	disp(sprintf('--- Gallager: n=%d, best_lambda = %g, log M >= %g', n, lambda, logm));

	if ( lambda <= 0 ) || (lambda >= 1)
		disp(sprintf('------- WARNING: lambda on the boundary! params: n=%d, delta=%g, epsil=%g', n, delta, epsil));
	end
	Lms = [Lms logm];
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function closeInplaneWindow()
%
% closeInplaneWindow()
%
% clears INPLANE while closing the corresponding window.
%
% djh, 1/26/98
mrGlobals

% Get the current figure
currentFigure = get(0,'CurrentFigure');

% Find and clear the corresponding INPLANE structure
try
    for s = 1:length(INPLANE)
        if checkfields(INPLANE{s}, 'ui', 'figNum')
            if INPLANE{s}.ui.figNum==currentFigure
                % let's try saving the prefs as we go,
                % so we pick up where we left off... (ras 01/06)
                if ispref('VISTA', 'savePrefs')
                    flag = getpref('VISTA', 'savePrefs');
                    if flag==1
                        savePrefs(INPLANE{s});
                    end
                end

                % if this was the selected inplane, well, it ain't no more.
                if selectedINPLANE == s
                    selectedINPLANE = [];
                end

                % finally, clear the global variable
                INPLANE{s} = [];
            end
        end
    end
end

% Delete the window
delete(currentFigure);

% If there are no open Inplane windows left, re-initialize it as empty
if isempty(cellfind(INPLANE)), INPLANE = {}; end

% Check if there are no more open views: if not, clean the workspace
if isempty(cellfind(INPLANE)) & isempty(cellfind(VOLUME)) & ...
   isempty(cellfind(FLAT))
    mrvCleanWorkspace;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft