text:
function [mg,vg,pg,pv]=gaussmixg(m,v,w,n)
%GAUSSMIXG global mean, variance and mode of a GMM
%
% Usage: (1) gaussmixg(m,v,w)               % plot the mean and mode positions of a GMM
%        (2) [mg,vg]=gaussmixg(m,v,w)       % find global mean and covariance of a GMM
%        (3) [mg,vg,pg]=gaussmixg(m,v,w)    % find global mean,covariance and mode of a GMM
%        (4) [mg,vg,pg,pv]=gaussmixg(m,v,w) % ... also find log probability of the peak
%
%  Inputs:  M(k,p) = mixture means for pg(p)
%           V(k,p) or V(p,p,k) variances (diagonal or full)
%           W(k,1) = mixture weights
%           N      = maximum number of modes to find [default 1]
%
% Outputs: MG(1,p) = global mean
%          VG(p,p) = global covariance
%          PG(N,p) = sorted list of N modes
%          PV(N,1) = log pdf at the modes PG(N,p) (in decreasing order)
%
%  This routine finds the global mean and covariance matrix of a Gaussian Mixture (GMM). It also
%  attempts to find up to N local maxima using a combination of the fixed point and quadratic
%  Newton-Raphson algorithms from [1]. Currently, N must be less than or equal to the number of
%  mixtures K. In general the PDF surface of a GMM can be very complicated with many local maxima [2]
%  and, as discussed in [1,2], this algorithm is not guaranteed to find the N highest. In [2], it is
%  conjectured that the number of local maxima is <=K for the following cases (a) P=1, (b) all
%  mixture covariance matrices are equal and (c) all mixture covariance matrices are multiples of
%  the identity.
%
% Refs:
%   [1]	M. Á. Carreira-Perpiñán. Mode-finding for mixtures of gaussian distributions.
%       IEEE Trans. Pattern Anal and Machine Intell, 22 (11): 13181323, 2000. doi: 10.1109/34.888716.
%   [2] M. Á. Carreira-Perpiñán and C. K. I. Williams. On the number of modes of a gaussian mixture.
%       In Proc Intl Conf on Scale Space Theories in Computer Vision, volume LNCS 2695, pages 625640,
%       Isle of Skye, June 2003. doi: 10.1007/3-540-44935-3_44.

% Bugs/Suggestions:
% (1) Sometimes the mode is not found, e.g. m=[0 1; 1 0];v=[.01 10; 10 .01];
%     has a true mode near (0,0). Could add to the list of mode candidates
%     all the pairwise intersections of the mixtures.
%     Another is: m=[0 0; 10 0.3]; v=[1 1; 1000 .001];
% (2) When merging candidates, we should keep the one with the highest probability
% (3) could preserve the fixed arrays between calls if p and/or k are unchanged
%
% See also: gaussmix, gaussmixd, gaussmixp, randvec

%      Copyright (C) Mike Brookes 2000-2012
%      Version: $Id: gaussmixg.m 3227 2013-07-04 15:42:04Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Algorithm parameters
nfp=2;          % number of fixed point iterations to do at start
maxloop=60;     % maximum number of iterations
ssf=0.1;        % factor when calculating minimum mode separation
%Sort out arguments
[k,p]=size(m);
if nargin<4
    n=1;
    if nargin<3
        w=ones(k,1);
        if nargin<2
            v=ones(k,p);
        end
    end
end
if ~nargout
    if nargin<4
        n=k;
    end
    nao=4;
else
    nao=nargout;
end

full=numel(size(v))>2 || k==1 && numel(v)>p; % test for full covariance matrices
if full && p==1   % if p=1 then we force diagonal covariance matrices
    v=reshape(v,1,k)';
    full=false;
end
w=w/sum(w);  % force w to sum to unity
% calculate the global mean and covariance
mg=w.'*m;
mz=m-mg(ones(k,1),:);   % means relative to global mean
if nao>2                % need to calculate the modes
    nx=k;               % number of pg values initially
    kk=reshape(repmat(1:nx,k,1),k*nx,1); % [1 1 1 2 2 2 ... nx nx nx]'
    km=repmat(1:k,1,nx)'; % [1:k 1:k ... 1:k]'
    % sort out indexing for all data value pairs; needed to eliminate duplicates
    nxp=nx*(nx-1)/2;
    ja=(1:nxp)';
    jb=floor((3+sqrt(8*ja-3))/2);       % [2 3 3 4 4 4 ... nx]'
    ja=ja-(jb.^2-3*jb+2)/2;             % [1 1:2 1:3 ... 1:nx-1]'
    jc=ones(nxp,1);
    % sort out indexing for vectorized upper triagular matrix
    npu=p*(p+1)/2;                      % number of distinct elements in a symmetrial (p,p) matrix
    kw=1:npu;
    ku=floor((1+sqrt(8*kw-3))/2);       % [1 2 2 3 3 3 ... p]
    kv=kw-(ku.^2-ku)/2;                 % [1 1:2 1:3 ... 1:p]
    zpp=zeros(p,p);
    zpp(kv+p*(ku-1))=kw;
    zpp(ku+p*(kv-1))=kw;
    kw=reshape(zpp,1,[]);               % maps vectorized upper triangular to vectorized full matrix
    kp=repmat(1:p,1,p);                 % row indices for a (p,p) matrix as a row vector
    kq=reshape(repmat(1:p,p,1),1,p^2);  % col indices for a (p,p) matrix as a row vector
    kr=p*kp-p+kq;                       % transpose indexing for a vectorized (p,p) matrix
    kd=1:p+1:p^2;                       % diagonal indices of a (p,p) matrix
    % unity vectors to make efficient replication
    wk=ones(k,1);
    wnx=ones(nx,1);
    if full
        vg=mz.'*(mz.*w(:,ones(1,p)))+reshape(reshape(v,p^2,k)*w,p,p);
        % now determine the mode
        vi=zeros(p*k,p);                    % stack of k inverse cov matrices each size p*p times -0.5
        vim=zeros(p*k,1);                   % stack of k vectors of the form -0.5*inv(vt)*m
        mtk=vim;                             % stack of k vectors of the form m
        lvm=zeros(k,1);
        wpk=repmat((1:p)',k,1);
        for i=1:k    % loop for each mixture
            [uvk,dvk]=eig(v(:,:,i));      % find eigenvalues
            dvk=diag(dvk);
            if(any(dvk<=0))
                error('Covariance matrix for mixture %d is not positive definite',i);
            end
            vik=-0.5*uvk*diag(dvk.^(-1))*uvk';   % calculate inverse including -0.5 factor
            vi((i-1)*p+(1:p),:)=vik;           % vi contains all mixture inverses stacked on top of each other
            vim((i-1)*p+(1:p))=vik*m(i,:)';   % vim contains vi*m for all mixtures stacked on top of each other
            mtk((i-1)*p+(1:p))=m(i,:)';       % mtk contains all mixture means stacked on top of each other
            lvm(i)=log(w(i))-0.5*sum(log(dvk));       % vm contains the weighted sqrt of det(vi) for each mixture
        end
        vif=reshape(permute(reshape(vi,p,k,p),[2 1 3]),k,p^2); % each covariance matrix as a vectorized row
        vimf=reshape(vim,p,k)'; % vi*m as a row for each mixture
        ss=sqrt(min(v(repmat(kd,k,1)+repmat(p^2*(0:k-1)',1,p)),[],1))*ssf/sqrt(p);   % minimum separation of modes [this is a conservative guess]
    else
        vg=mz.'*(mz.*w(:,ones(1,p)))+diag(w.'*v);
        % now determine the mode
        vi=-0.5*v.^(-1);                % vi(k,p) = data-independent scale factor in exponent
        vi2=vi(:,ku).*vi(:,kv);         % vi2(k,npu) = upper triangular Hessian data dependent term
        lvm=log(w)-0.5*sum(log(v),2);   % log of external scale factor (excluding -0.5*p*log(2pi) term)
        vim=vi.*m;
        vim2=vim(:,ku).*vim(:,kv);      % vim2(k,npu) = upper triangular Hessian data independent term
        vimvi=vim(:,kp).*vi(:,kq);      % vimvi(k,p^2) = vectorized Hessian term
        ss=sqrt(min(v,[],1))*ssf/sqrt(p);   % minimum separation of modes [this is a conservative guess]
    end
    pgf=zeros(nx,p);                    % space for fixed point update
    sv=0.01*ss;                         % convergence threshold
    pg=m;  % initialize mode candidates to mixture means
    i=1;   % loop counter
    %     gx=zeros(nx,p);   %%%%%%%%%%%% temp
    while i<=maxloop
        %         pg00=pg0;   %%%%%%%%%%%% temp
        pg0=pg;                         % save previous mode candidates pg(nx,p)
        if full
            py=reshape(sum(reshape((vi*pg'-vim(:,wnx)).*(pg(:,wpk)'-mtk(:,wnx)),p,nx*k),1),k,nx)+lvm(:,wnx);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=(px./ps(wk,:))';             % px(nx,k) = relative mixture probabilities for each data point (rows sum to 1)
            % calculate the fixed point update
            pxvif=px*vif;    % pxvif(nx,p^2)
            pxvimf=px*vimf;  % pxvimf(nx,p)
            for j=1:nx
                pgf(j,:)=pxvimf(j,:)/reshape(pxvif(j,:),p,p);
            end
        else
            py=reshape(sum((pg(kk,:)-m(km,:)).^2.*vi(km,:),2),k,nx)+lvm(:,wnx);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=(px./ps(wk,:))';             % px(nx,k) = relative mixture probabilities for each data point (rows sum to 1)
            % calculate the fixed point update
            pxvim=px*vim;
            pxvi=px*vi;
            pgf=pxvim./pxvi;       % fixed point update for all points
        end
        if i>nfp
            % calculate gradient and Hessian; see [1] equations (4) and (5)
            lp=log(ps)+mx;              % log prob of each data point
            if full
                %                 gx0=gx;   %%%%%%%%%%%% temp
                gx=pxvimf-reshape(sum(repmat(pg,p,1).*reshape(pxvif,[],p),2),nx,p);
                vimpg=repmat(vimf,nx,1)-reshape(permute(reshape(pg*vi',nx,p,k),[3 1 2]),[],p); % vimpg(k*nx,p)
                hx1=2*reshape(sum(reshape(repmat(reshape(px',[],1),1,npu).*vimpg(:,ku).*vimpg(:,kv),k,[]),1),nx,[]);
                hx=pxvif+hx1(:,kw);
            else
                gx=pxvim-pxvi.*pg;               % gradient for each data point (one row per point)
                hx1=px*vim2+(px*vi2).*pg(:,ku).*pg(:,kv);
                hx2=(px*(vimvi)).*pg(:,kq);
                hx=2*(hx1(:,kw)-hx2-hx2(:,kr));
                hx(:,kd)=hx(:,kd)+pxvi;
            end
            hx=reshape(hx',p,p,nx);
            for j=1:nx
                if all(eig(hx(:,:,j))<0)  % if positive definite
                    pg(j,:)=pg(j,:)+gx(j,:)/hx(:,:,j); % do a Newton-Raphson update
                    if full
                        pyj=sum(reshape((vi*pg(j,:)'-vim).*(pg(j,wpk)'-mtk),p,k),1)'+lvm;
                    else
                        pyj=sum((repmat(pg(j,:),k,1)-m).^2.*vi,2)+lvm;
                    end
                    mxj=max(pyj);                % find normalizing factor for each data point to prevent underflow when using exp()
                    pxj=exp(pyj-mxj);            % find normalized probability of each mixture for each datapoint
                    psj=sum(pxj,1);                   % total normalized likelihood of each data point
                    lpj=log(psj)+mxj;              % log prob of updated data point
                    if lpj<lp(j)       % check if the probability has decreased
                        pg(j,:)=pgf(j,:);   % if so, do fixed point update
                    end
                else
                    pg(j,:)=pgf(j,:);   % else do fixed point update
                end
            end
        else
            pg=pgf;       % fixed point update for all points
        end
        if all(all(abs(pg-pg0)<sv(wnx,:))) && i+2<maxloop
            maxloop=min(maxloop,i+2);        % two more loops if converged if converged
        end
        %         [all(pg==pgf,2) pg [i; repmat(NaN,nx-1,1)] pg-pg0]  %debug: [fixed-point x-mode iteration delta-x]
        jd=all(abs(pg(jb,:)-pg(ja,:))<ss(jc,:),2);   % find duplicate modes
        if any(jd)
            jx=sparse([(1:nx)';ja;jb],[(1:nx)';jb;ja],[wnx;jd;jd]);   % neighbour matrix
            kx=any((jx*jx)>0 & ~jx,2);  % find chains that  are not fully connected
            while any(kx)
                kx=any(jx(:,kx),2);     % destroy all links connected to these chains
                jx(kx,:)=0;
                jx(:,kx)=0;
                % jx(kx,kx)=1;
                kx=any((jx*jx)>0 & ~jx,2);
            end
            jx([1:nx+1:nx*nx (ja+(jb-1)*nx)'])=0; % reset the upper triangle + diagonal
            pg(any(jx,2),:)=[];   % delete the duplicates
            % update nx and anything that depends on it
            nx=size(pg,1);
            pgf=zeros(nx,p);                    % space for fixed point update
            wnx=ones(nx,1);
            nxp=nx*(nx-1)/2;
            ja=ja(1:nxp);
            jb=jb(1:nxp);
            kk=reshape(repmat(1:nx,k,1),k*nx,1); % [1 1 1 2 2 2 ... nx nx nx]'
            km=reshape(repmat(1:k,1,nx),k*nx,1); % [1:k 1:k ... 1:k]'
            jc=ones(nxp,1);
        end
        i=i+1;
    end
    %     calculate the log pdf at each mode
    if full
        py=reshape(sum(reshape((vi*pg'-vim(:,wnx)).*(pg(:,wpk)'-mtk(:,wnx)),p,nx*k),1),k,nx)+lvm(:,wnx);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
    else
        py=reshape(sum((pg(kk,:)-m(km,:)).^2.*vi(km,:),2),k,nx)+lvm(:,wnx);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
    end
    [pv,ix]=sort((log(ps)+mx)'-0.5*p*log(2*pi),'descend');
    pg=pg(ix,:);
    if n<numel(pv) % only keep the first n modes
        pg=pg(1:n,:);
        pv=pv(1:n);
    end
elseif nao>1
    if full
        vg=mz.'*(mz.*w(:,ones(1,p)))+reshape(reshape(v,p^2,k)*w,p,p);
    else
        vg=mz.'*(mz.*w(:,ones(1,p)))+diag(w.'*v);
    end
end

if ~nargout
    % now plot the result
    clf;
    pg1=pg(1,:);
    lpm=gaussmixp(mg,m,v,w);
    lpp=pv(1);
    switch p
        case 1
            gaussmixp([],m,v,w);
            hold on;
            ylim=get(gca,'ylim')';
            plot([mg mg]',ylim,'-k',[mg mg; mg mg]+[-1 1; -1 1]*sqrt(vg),ylim,':k');
            plot(pg(1),pv(1),'^k');
            if numel(pg)>1
                plot(pg(2:end),pv(2:end),'xk');
            end
            hold off;
            title(sprintf('Mean+-sd = %.3g+-%.3g LogP = %.3g, Mode\\Delta = %.3g LogP = %.3g',mg,sqrt(vg),lpm,pg1,lpp));
            xlabel('x');
        case 2
            gaussmixp([],m,v,w);
            hold on;
            t=linspace(0,2*pi,100);
            xysd=chol(vg)'*[cos(t); sin(t)]+repmat(mg',1,length(t));
            plot(xysd(1,:),xysd(2,:),':k',mg(1),mg(2),'ok');
            plot(pg(1,1),pg(1,2),'^k');
            if numel(pv)>1
                plot(pg(2:end,1),pg(2:end,2),'xk');
            end
            hold off;
            title(sprintf('Mean = (%.3g,%.3g) LogP = %.3g, Mode:\\Delta = (%.3g,%.3g) LogP = %.3g',mg,lpm,pg1,pv(1)));
            xlabel('x');
            ylabel('y');
        otherwise
            nx=200;
            nc=ceil(sqrt(p/2));
            nr=ceil(p/nc);
            sdx=sqrt(diag(vg))';  % std deviation
            minx=min([mg; pg],[],1)-1.5*sdx;
            maxx=max([mg; pg],[],1)+1.5*sdx;
            ix=2:p; % selected indices
            for i=1:p
                xi=linspace(minx(i),maxx(i),nx)';
                [mm,vm,wm]=gaussmixd(mg(ix),m,v,w,[],ix);
                ym=gaussmixp(xi,mm,vm,wm')+lpm-gaussmixp(mg(i),mm,vm,wm');
                [mp,vp,wp]=gaussmixd(pg(1,ix),m,v,w,[],ix);
                yp=gaussmixp(xi,mp,vp,wp')+lpp-gaussmixp(pg1(i),mp,vp,wp');
                subplot(nr,nc,i);
                plot(xi,ym,'-k',mg(i),lpm,'ok',xi,yp,':b',pg1(i),lpp,'^b');
                axisenlarge([-1 -1 -1 -1.05]);
                hold on
                plot([mg(i) mg(i)],get(gca,'ylim'),'-k');
                hold off
                xlabel(sprintf('x[%d], Mean = %.3g, Mode\\Delta = %.3g',i,mg(i),pg1(i)));
                if i==1
                    title(sprintf('Log Prob: Mean = %.3g, Mode\\Delta = %.3g',lpm,lpp));
                end
                ix(i)=i;
            end
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function theRank=rankTComposition(p,firstElMostSig)
%%RANKTCOMPOSITION Obtain the lexicographic order (the rank) of a particular
%                 composition starting at 0. A composition of n into t
%                 is a way of putting n unlabeled balls into t labeled
%                 slots. Compositions are tuples of t integers >0 that
%                 sum to n. Compositions are given in lexicographic order
%                 where the first element is the least significant unless
%                 otherwise specified.
%
%INPUTS: p An tX1 vector holding the current composition, whose elements
%          sum to n>=1.
% firstElMostSig An optional parameter specifying whether the first element
%          is the most or least significant. The default if omitted or an
%          empty matrix is passed is false (the first element is the least
%          significant).
%
%OUTPUTS: theRank The rank of the composition in a lexicographic ordering
%                 of compositions, counting from zero. There is a total of
%                 binomial(n-1,t-1) compositions. Lexicographic ordering is
%                 defined using p(1) as the least significant element.
%
%This is an implementation of the relation described in Chapter 7.2.1.3 of
%[1] for mapping combinations to compositions. However, a special case is
%if the input is one-dimensional. In such an instance, the rank is always
%0.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Special case, empty matrix.
if(isempty(p))
   theRank=[];
   return;
end

%Special case, given a scalar.
if(length(p)==1)
    theRank=0;
   return;
end

if(nargin<2||isempty(firstElMostSig))
    firstElMostSig=false;
end

if(firstElMostSig)
   p=flipud(p(:));
end

m=length(p);

%Turn the composition into the equivalent combination.
mCombo=m-1;

c=zeros(mCombo,1);
c(1)=p(1)-1;

for k=2:mCombo
    c(k)=p(k)+c(k-1);
end

%get the rank of the combination.
theRank=rankCombination(c);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y = LPFilter(x,fc)
%
% y = LPFilter(x,fc),
% Overview: Second order zero-phase Lowpass filter. The filter is used to
% estimate the baseline wander in the input ECG signal.
%
% inputs:
% x: vector or matrix of input ECG data (channels x samples)
% fc: -3dB cut-off frequency normalized by the sampling frequency
%
% output:
% y: vector or matrix of filtered ECG data (channels x samples)
%
%
% Open Source ECG Toolbox, version 1.0, November 2006
% Released under the GNU General Public License
% Copyright (C) 2006  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- GIPSA-Lab, INPG, Grenoble, France
% reza.sameni@gmail.com

% 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

k = .7;     % cut-off value
alpha = (1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc))-k^2*sin(2*pi*fc)^2))/(1-k);
y = zeros(size(x));
for i = 1:size(x,1)
    y(i,:) = filtfilt(1-alpha,[1 -alpha],x(i,:));
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [Offspring] = DEgenerator2(Problem,Population)

%------------------------------- 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".
%--------------------------------------------------------------------------

    cv = sum(max(0,Population.cons),2);

    FrontNo = NDSort(Population.objs,Population.cons,1);
    index1  = find(FrontNo==1);
    r       = floor(rand*length(index1))+1;
    best    = index1(r);

    [N,D] = size(Population(1).decs);
    trial = zeros(1*Problem.N,D);

    for i = 1 : Problem.N
        % DE/rand-to-bestf/1/bin--Convergence
        if rand > 0.5
            l = rand;
            if l <= 1/3
            	F = .6;
            elseif l <= 2/3
                F = 0.8;
            else
                F = 1.0;
            end
            l = rand;
            if l <= 1/3
                CR = .1;
            elseif l <= 2/3
                CR = 0.2;
            else
                CR = 1.0;
            end
            indexset    = 1 : Problem.N;
            indexset(i) = [];
            r1  = floor(rand*(Problem.N-1))+1;
            xr1 = indexset(r1);
            indexset(r1) = [];
            r2  = floor(rand*(Problem.N-2))+1;
            xr2 = indexset(r2)  ;
            r3  = floor(rand*(Problem.N-3))+1;
            xr3 = indexset(r3);
            Best_index = Population(best).decs;
            v      = Population(xr1).decs+rand*(Best_index-Population(xr1).decs)+F*(Population(xr2).decs-Population(xr3).decs);
            Lower  = repmat(Problem.lower,N,1);
            Upper  = repmat(Problem.upper,N,1);
            v      = min(max(v,Lower),Upper);
            Site   = rand(N,D) < CR;
            j_rand = floor(rand * D) + 1;
            Site(1, j_rand) = 1;
            Site_  = 1-Site;
            trial(i, :) = Site.*v+Site_.*Population(i).decs;
        % DE/current-to-rand/1--Diversity
        else
            l = rand;
            if l <= 1/3
                F = .6;
            elseif l <= 2/3
                F = 0.8;
            else
                F = 1.0;
            end
            indexset    = 1:Problem.N;
            indexset(i) = [];
            r1  = floor(rand*(Problem.N-1))+1;
            xr1 = indexset(r1);
            indexset(r1) = [];
            r2  = floor(rand*(Problem.N-2))+1;
            xr2 = indexset(r2);
            indexset(r2) = [];
            r3    = floor(rand*(Problem.N-3))+1;
            xr3   = indexset(r3);
            v     = Population(i).decs+rand*(Population(xr1).decs-Population(i).decs)+F*(Population(xr2).decs-Population(xr3).decs);
            Lower = repmat(Problem.lower,N,1);
            Upper = repmat(Problem.upper,N,1);
            trial(i, :) = min(max(v,Lower),Upper);
        end
    end
    Offspring = trial;
    Offspring = Problem.Evaluation(Offspring);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function y=v_teager(x,d,m)
%V_TEAGER calculate v_teager energy waveform Y=(X,D,M)
%
%  Inputs:  x         speech signal
%           d         dimension to apply filter along [default 1st non-singleton]
%           m         Normally Y has the same length as X and the first
%                     and last output samples are extrapolated. Setting m='x'
%                     supresses this extrapolation and Y will be two
%                     samples shorter than X
%
% Outputs:  Y         output signal: y(n)=abs(x(n))^2 - x(n+1)*conj(x(n-1))
%
% Calculates the Teager energy waveform [1]. The following waveforms give
% a constant output (independent of n) where A, B, C are real constants:
%  (a) x(n) = A*sin(B*n+C)    -->   y(n) = (A*sin(B))^2
%  (b) x(n) = A*n + B         -->   y(n) = A^2
%  (c) x(n) = A*exp(j(B*n+C)) -->   y(n) = A^2*(1-exp(2jB))
%  (d) x(n) = A*exp(B*n+C)    -->   y(n) = 0
%
% Reference:
%  [1]	J. Kaiser. On a simple algorithm to calculate the "energy" of a signal.
%       In Proc IEEE Intl Conf Acoustics, Speech and Signal Processing,
%       pages 381-384, vol.1, Apr. 1990. doi: 10.1109/ICASSP.1990.115702.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_teager.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e=size(x);
p=prod(e);
if nargin<2             % if no dimension given, find the first non-singleton
    d=find(e>1,1);
    if ~numel(d)
        d=1;
    end
end
k=e(d);                 % size of active dimension
q=p/k;                  % size of remainder
if d==1
    z=reshape(x,k,q);
else
    z=shiftdim(x,d-1);
    r=size(z);
    z=reshape(z,k,q);
end
if nargin>2 && any(m=='x')
    y=z(2:k-1,:).*conj(z(2:k-1,:))-z(3:k,:).*conj(z(1:k-2,:));
    k=k-2;              % we have lost two elements
elseif k>=4
    y=zeros(k,q);
    y(2:k-1,:)=z(2:k-1,:).*conj(z(2:k-1,:))-z(3:k,:).*conj(z(1:k-2,:));
    y(1,:)=2*y(2,:)-y(3,:);             % linearly interpolate the end points
    y(k,:)=2*y(k-1,:)-y(k-2,:);
elseif k==3
    y=repmat(x(2,:).*conj(x(2,:))-x(3,:).*conj(x(1,:)),3,1);
else
    y=zeros(k,q);
end
if d==1
    e(d)=k;
    y=reshape(y,e);
else
    r(1)=k;
    y=shiftdim(reshape(y,r),length(e)+1-d);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function lik = lik_loglogistic(varargin)
%LIK_LOGLOGISTIC Create a right censored log-logistic likelihood structure
%
%  Description
%    LIK = LIK_LOGLOGISTIC('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a likelihood structure for right censored log-logistic
%    survival model in which the named parameters have the
%    specified values. Any unspecified parameters are set to
%    default values.
%
%    LIK = LIK_LOGLOGISTIC(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood structure with the named parameters
%    altered with the specified values.
%
%    Parameters for loggaussian likelihood [default]
%      shape       - shape parameter r [1]
%      shape_prior - prior for shape [prior_logunif]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                  __ n
%      p(y|f, z) = || i=1 [ (r/exp(f_i)*(y_i/exp(f_i))^(r-1)/
%                               (1+(y_i/exp(f_i))^r))^(1-z_i)
%                               *(1+(y_i/exp(f_i))^r)^(-z_i) ]
%
%
%    where r is the shape parameter of log-logistic distribution.
%    z is a vector of censoring indicators with z = 0 for uncensored event
%    and z = 1 for right censored event.
%
%    When using the loggaussian likelihood you need to give the vector z
%    as an extra parameter to each function that requires also y.
%    For example, you should call gpla_e as follows: gpla_e(w, gp,
%    x, y, 'z', z)
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%

% Copyright (c) 2012 Ville Tolvanen

% 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.

  ip=inputParser;
  ip.FunctionName = 'LIK_LOGLOGISTIC';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('shape',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('shape_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Log-Logistic';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Log-Logistic')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('shape',ip.UsingDefaults)
    lik.shape = ip.Results.shape;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('shape_prior',ip.UsingDefaults)
    lik.p.shape=ip.Results.shape_prior;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_loglogistic_pak;
    lik.fh.unpak = @lik_loglogistic_unpak;
    lik.fh.lp = @lik_loglogistic_lp;
    lik.fh.lpg = @lik_loglogistic_lpg;
    lik.fh.ll = @lik_loglogistic_ll;
    lik.fh.llg = @lik_loglogistic_llg;
    lik.fh.llg2 = @lik_loglogistic_llg2;
    lik.fh.llg3 = @lik_loglogistic_llg3;
    lik.fh.tiltedMoments = @lik_loglogistic_tiltedMoments;
    lik.fh.siteDeriv = @lik_loglogistic_siteDeriv;
    lik.fh.invlink = @lik_loglogistic_invlink;
    lik.fh.predy = @lik_loglogistic_predy;
    lik.fh.recappend = @lik_loglogistic_recappend;
    lik.fh.predcdf=@lik_loglogistic_predcdf;
  end

end

function [w,s] = lik_loglogistic_pak(lik)
%LIK_LOGLOGISTIC_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_LOGLOGISTIC_PAK(LIK) takes a likelihood structure LIK and
%    combines the parameters into a single row vector W. This is a
%    mandatory subfunction used for example in energy and gradient
%    computations.
%
%       w = log(lik.shape)
%
%   See also
%   LIK_LOGLOGISTIC_UNPAK, GP_PAK

  w=[];s={};
  if ~isempty(lik.p.shape)
    w = log(lik.shape);
    s = [s; 'log(loglogistic.shape)'];
    [wh sh] = lik.p.shape.fh.pak(lik.p.shape);
    w = [w wh];
    s = [s; sh];
  end
end

function [lik, w] = lik_loglogistic_unpak(lik, w)
%LIK_LOGLOGISTIC_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    [LIK, W] = LIK_LOGLOGISTIC_UNPAK(W, LIK) takes a likelihood
%    structure LIK and extracts the parameters from the vector W
%    to the LIK structure. This is a mandatory subfunction used
%    for example in energy and gradient computations.
%
%   Assignment is inverse of
%       w = log(lik.shape)
%
%   See also
%   LIK_LOGLOGISTIC_PAK, GP_UNPAK

  if ~isempty(lik.p.shape)
    lik.shape = exp(w(1));
    w = w(2:end);
    [p, w] = lik.p.shape.fh.unpak(lik.p.shape, w);
    lik.p.shape = p;
  end
end

function lp = lik_loglogistic_lp(lik, varargin)
%LIK_LOGLOGISTIC_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_LOGLOGISTIC_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This subfunction
%    is needed when there are likelihood parameters.
%
%  See also
%    LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, LIK_LOGLOGISTIC_LLG2, GPLA_E

% If prior for shape parameter, add its contribution
  lp=0;
  if ~isempty(lik.p.shape)
    lp = lik.p.shape.fh.lp(lik.shape, lik.p.shape) +log(lik.shape);
  end

end

function lpg = lik_loglogistic_lpg(lik)
%LIK_LOGLOGISTIC_LPG  d log(prior)/dth of the likelihood
%                parameters th
%
%  Description
%    E = LIK_LOGLOGISTIC_LPG(LIK) takes a likelihood structure LIK and
%    returns d log(p(th))/dth, where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, LIK_LOGLOGISTIC_LLG2, GPLA_G

  lpg=[];
  if ~isempty(lik.p.shape)
    % Evaluate the gprior with respect to shape
    ggs = lik.p.shape.fh.lpg(lik.shape, lik.p.shape);
    lpg = ggs(1).*lik.shape + 1;
    if length(ggs) > 1
      lpg = [lpg ggs(2:end)];
    end
  end
end

function ll = lik_loglogistic_ll(lik, y, f, z)
%LIK_LOGLOGISTIC_LL  Log likelihood
%
%  Description
%    LL = LIK_LOGLOGISTIC_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z, and
%    latent values F. Returns the log likelihood, log p(y|f,z).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also
%    LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, LIK_LOGLOGISTIC_LLG2, GPLA_E

  if isempty(z)
    error(['lik_loglogistic -> lik_loglogistic_ll: missing z!    '...
           'loglogistic likelihood needs the censoring    '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loglogistic and gpla_e.               ']);
  end

  r = lik.shape;
  ll = sum((1-z).*(log(r)+(r-1).*log(y)-r.*f) +(z-2).*log(1+(y./exp(f)).^r));

end

function llg = lik_loglogistic_llg(lik, y, f, param, z)
%LIK_LOGLOGISTIC_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_LOGLOGISTIC_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z and
%    latent values F. Returns the gradient of the log likelihood
%    with respect to PARAM. At the moment PARAM can be 'param' or
%    'latent'. This subfunction is needed when using Laplace
%    approximation or MCMC for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_LLG2, LIK_LOGLOGISTIC_LLG3, GPLA_E

  if isempty(z)
    error(['lik_loglogistic -> lik_loglogistic_llg: missing z!    '...
           'loglogistic likelihood needs the censoring    '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loglogistic and gpla_e.               ']);
  end

  r = lik.shape;
  switch param
    case 'param'
      llg = sum((1-z).*(1/r+log(y)-f) + (z-2)./(1+(y./exp(f)).^r).* ...
             (y./exp(f)).^r.*(log(y)-f));
      % correction for the log transformation
      llg = llg.*lik.shape;
    case 'latent'
      llg = -r.*(1-z) - (z-2).*r.*(y./exp(f)).^r./(1+(y./exp(f)).^r);
  end
end

function llg2 = lik_loglogistic_llg2(lik, y, f, param, z)
%LIK_LOGLOGISTIC_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_LOGLOGISTIC_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z, and
%    latent values F. Returns the hessian of the log likelihood
%    with respect to PARAM. At the moment PARAM can be only
%    'latent'. LLG2 is a vector with diagonal elements of the
%    Hessian matrix (off diagonals are zero). This subfunction
%    is needed when using Laplace approximation or EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, GPLA_E

  if isempty(z)
    error(['lik_loglogistic -> lik_loglogistic_llg2: missing z!   '...
           'loglogistic likelihood needs the censoring   '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loglogistic and gpla_e.               ']);
  end

  r = lik.shape;
  switch param
    case 'param'

    case 'latent'
      llg2 = r.^2.*(z-2).*(y./exp(f)).^r./(1+(y./exp(f)).^r).^2;
    case 'latent+param'
      llg2 = (z-1) - (z-2).*(y./exp(f)).^r./(1+(y./exp(f)).^r) ...
               + (z-2).*r.*(y./exp(f)).^(2*r).*(log(y)-f)./(1+(y./exp(f)).^r).^2 ...
               - (z-2).*r.*(y./exp(f)).^r.*(log(y)-f)./(1+(y./exp(f)).^r);
      % correction due to the log transformation
      llg2 = llg2.*r;
  end
end

function llg3 = lik_loglogistic_llg3(lik, y, f, param, z)
%LIK_LOGLOGISTIC_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_LOGLOGISTIC_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z and
%    latent values F and returns the third gradients of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. LLG3 is a vector with third gradients. This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG2, GPLA_E, GPLA_G

  if isempty(z)
    error(['lik_loglogistic -> lik_loglogistic_llg3: missing z!   '...
           'loglogistic likelihood needs the censoring    '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loglogistic and gpla_e.               ']);
  end

  r = lik.shape;
  switch param
    case 'param'

    case 'latent'
      llg3 = r.^3.*(z-2).*(y./exp(f)).^r.*(-1+(y./exp(f)).^r)./(1+(y./exp(f)).^r).^3;
    case 'latent2+param'
      llg3 = -(r.*(z-2).*(y./exp(f)).^r.*(-2-2.*(y./exp(f)).^r + ...
              r.*(-1+(y./exp(f)).^r).*log(y./exp(f))))./(1+(y./exp(f)).^r).^3;
      % correction due to the log transformation
      llg3 = llg3.*lik.shape;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_loglogistic_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
%LIK_LOGLOGISTIC_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_LOGLOGISTIC_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, survival times
%    Y, censoring indicators Z, index I and cavity variance S2 and
%    mean MYY. Returns the zeroth moment M_0, mean M_1 and
%    variance M_2 of the posterior marginal (see Rasmussen and
%    Williams (2006): Gaussian processes for Machine Learning,
%    page 55). This subfunction is needed when using EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

 if isempty(z)
   error(['lik_loglogistic -> lik_loglogistic_tiltedMoments: missing z!'...
          'loglogistic likelihood needs the censoring            '...
          'indicators as an extra input z. See, for                 '...
          'example, lik_loglogistic and gpep_e.                       ']);
 end

  yy = y(i1);
  yc = 1-z(i1);
  r = lik.shape;
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    % get a function handle of an unnormalized tilted distribution
    % (likelihood * cavity = Negative-binomial * Gaussian)
    % and useful integration limits
    [tf,minf,maxf]=init_loglogistic_norm(yy(i),myy_i(i),sigm2_i(i),yc(i),r);

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigm2_i(i)
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      if sigm2hati1(i) >= sigm2_i(i)
        error('lik_loglogistic_tilted_moments: sigm2hati1 >= sigm2_i');
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [g_i] = lik_loglogistic_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_LOGLOGISTIC_SITEDERIV  Evaluate the expectation of the gradient
%                      of the log likelihood term with respect
%                      to the likelihood parameters for EP
%
%  Description [M_0, M_1, M2] =
%    LIK_LOGLOGISTIC_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIK, survival times Y, expected
%    counts Z, index I and cavity variance S2 and mean MYY.
%    Returns E_f [d log p(y_i|f_i) /d a], where a is the
%    likelihood parameter and the expectation is over the
%    marginal posterior. This term is needed when evaluating the
%    gradients of the marginal likelihood estimate Z_EP with
%    respect to the likelihood parameters (see Seeger (2008):
%    Expectation propagation for exponential families). This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods and there are likelihood parameters.
%
%  See also
%    GPEP_G

  if isempty(z)
    error(['lik_loglogistic -> lik_loglogistic_siteDeriv: missing z!'...
           'loglogistic likelihood needs the censoring        '...
           'indicators as an extra input z. See, for             '...
           'example, lik_loglogistic and gpla_e.                   ']);
  end

  yy = y(i1);
  yc = 1-z(i1);
  r = lik.shape;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Log-Gaussian * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_loglogistic_norm(yy,myy_i,sigm2_i,yc,r);
  % additionally get function handle for the derivative
  td = @deriv;

  % Integrate with quadgk
  [m_0, fhncnt] = quadgk(tf, minf, maxf);
  [g_i, fhncnt] = quadgk(@(f) td(f).*tf(f)./m_0, minf, maxf);
  g_i = g_i.*r;

  function g = deriv(f)
    g = yc.*(1/r+log(yy)-f) + (-1-yc)./(1+(yy./exp(f)).^r).* ...
             (yy./exp(f)).^r.*(log(yy)-f);
  end
end

function [lpy, Ey, Vary] = lik_loglogistic_predy(lik, Ef, Varf, yt, zt)
%LIK_LOGLOGISTIC_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_LOGLOGISTIC_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the survival times YT, censoring indicators ZT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_LOGLOGISTIC_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%
%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  if isempty(zt)
    error(['lik_loglogistic -> lik_loglogistic_predy: missing zt!'...
           'loglogistic likelihood needs the censoring    '...
           'indicators as an extra input zt. See, for         '...
           'example, lik_loglogistic and gpla_e.               ']);
  end

  yc = 1-zt;
  r = lik.shape;

  Ey=[];
  Vary=[];

  % Evaluate the posterior predictive densities of the given observations
  lpy = zeros(length(yt),1);
  for i1=1:length(yt)
    if abs(Ef(i1))>700
      lpy(i1) = NaN;
    else
      % get a function handle of the likelihood times posterior
      % (likelihood * posterior = Negative-binomial * Gaussian)
      % and useful integration limits
      [pdf,minf,maxf]=init_loglogistic_norm(...
        yt(i1),Ef(i1),Varf(i1),yc(i1),r);
      % integrate over the f to get posterior predictive distribution
      lpy(i1) = log(quadgk(pdf, minf, maxf));
    end
  end
end

function [df,minf,maxf] = init_loglogistic_norm(yy,myy_i,sigm2_i,yc,r)
%INIT_LOGLOGISTIC_NORM
%
%  Description
%    Return function handle to a function evaluating
%    loglogistic * Gaussian which is used for evaluating
%    (likelihood * cavity) or (likelihood * posterior) Return
%    also useful limits for integration. This is private function
%    for lik_loglogistic. This subfunction is needed by subfunctions
%    tiltedMoments, siteDeriv and predy.
%
%  See also
%    LIK_LOGLOGISTIC_TILTEDMOMENTS, LIK_LOGLOGISTIC_SITEDERIV,
%    LIK_LOGLOGISTIC_PREDY

% avoid repetitive evaluation of constant part
  ldconst =  yc.*log(r) + yc.*(r-1).*log(yy) ...
            - log(sigm2_i)/2 - log(2*pi)/2;

  % Create function handle for the function to be integrated
  df = @loglogistic_norm;
  % use log to avoid underflow, and derivates for faster search
  ld = @log_loglogistic_norm;
  ldg = @log_loglogistic_norm_g;
  ldg2 = @log_loglogistic_norm_g2;

  % Set the limits for integration
  if yc==0
    % with yy==0, the mode of the likelihood is not defined
    % use the mode of the Gaussian (cavity or posterior) as a first guess
    modef = myy_i;
  else
    % use precision weighted mean of the Gaussian approximation
    % of the loglogistic likelihood and Gaussian
    mu=-log(yc.^(1/r)./yy);
    s2=-r^2.*yc.*(-1-yc)./(1+yc).^2;
%     s2=1;
    modef = (myy_i/sigm2_i + mu/s2)/(1/sigm2_i + 1/s2);
  end
  % find the mode of the integrand using Newton iterations
  % few iterations is enough, since the first guess in the right direction
  niter=4;       % number of Newton iterations
  mindelta=1e-6; % tolerance in stopping Newton iterations
  for ni=1:niter
    g=ldg(modef);
    h=ldg2(modef);
    delta=-g/h;
    modef=modef+delta;
    if abs(delta)<mindelta
      break
    end
  end
  % integrand limits based on Gaussian approximation at mode
  modes=sqrt(-1/h);
  minf=modef-8*modes;
  maxf=modef+8*modes;
  modeld=ld(modef);
  iter=0;
  % check that density at end points is low enough
  lddiff=20; % min difference in log-density between mode and end-points
  minld=ld(minf);
  step=1;
  while minld>(modeld-lddiff)
    minf=minf-step*modes;
    minld=ld(minf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_loglogistic -> init_loglogistic_norm: ' ...
             'integration interval minimun not found ' ...
             'even after looking hard!'])
    end
  end
  maxld=ld(maxf);
  step=1;
  while maxld>(modeld-lddiff)
    maxf=maxf+step*modes;
    maxld=ld(maxf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_loglogistic -> init_loglogistic_norm: ' ...
             'integration interval maximun not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = loglogistic_norm(f)
  % loglogistic * Gaussian
    integrand = exp(ldconst ...
                     - yc.*r.*f +(-1-yc).*log(1+(yy./exp(f)).^r) ...
                    -0.5*(f-myy_i).^2./sigm2_i);
  end

  function log_int = log_loglogistic_norm(f)
  % log(loglogistic * Gaussian)
  % log_loglogistic_norm is used to avoid underflow when searching
  % integration interval
    log_int = ldconst ...
               -yc.*r.*f +(-1-yc).*log(1+(yy./exp(f)).^r) ...
              -0.5*(f-myy_i).^2./sigm2_i;
  end

  function g = log_loglogistic_norm_g(f)
  % d/df log(loglogistic * Gaussian)
  % derivative of log_loglogistic_norm
    g = -r.*yc - (-1-yc).*r.*(yy./exp(f)).^r./(1+(yy./exp(f)).^r) ...
        + (myy_i - f)./sigm2_i;
  end

  function g2 = log_loglogistic_norm_g2(f)
  % d^2/df^2 log(loglogistic * Gaussian)
  % second derivate of log_loglogistic_norm
    g2 =  r.^2.*(-1-yc).*(yy./exp(f)).^r./(1+(yy./exp(f)).^r).^2 ...
              -1/sigm2_i;
  end

end

function cdf = lik_loglogistic_predcdf(lik, Ef, Varf, yt)
%LIK_LOGLOGISTIC_PREDCDF  Returns the predictive cdf evaluated at yt
%
%  Description
%    CDF = LIK_LOGLOGISTIC_PREDCDF(LIK, EF, VARF, YT)
%    Returns the predictive cdf evaluated at YT given likelihood
%    structure LIK, posterior mean EF and posterior Variance VARF
%    of the latent variable. This subfunction is needed when using
%    functions gp_predcdf or gp_kfcv_cdf.
%
%  See also
%    GP_PREDCDF

  r = lik.shape;

  % Evaluate the posterior predictive densities of the given observations
  cdf = zeros(length(yt),1);
  for i1=1:length(yt)
    % Get a function handle of the likelihood times posterior
    % (likelihood * posterior = log-logistic * Gaussian)
    % and useful integration limits.
    % yc=0 when evaluating predictive cdf
    [pdf,minf,maxf]=init_loglogistic_norm(...
      yt(i1),Ef(i1),Varf(i1),0,r);
    % integrate over the f to get posterior predictive distribution
    cdf(i1) = 1-quadgk(pdf, minf, maxf);
  end
end

function p = lik_loglogistic_invlink(lik, f)
%LIK_loglogistic Returns values of inverse link function
%
%  Description
%    P = LIK_LOGLOGISTIC_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values of inverse link function P.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_PREDY

p = exp(f);
end

function reclik = lik_loglogistic_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = GPCF_LOGLOGISTIC_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    % Initialize the record
    reclik.type = 'Log-Logistic';

    % Initialize parameter
    reclik.shape = [];

    % Set the function handles
    reclik.fh.pak = @lik_loglogistic_pak;
    reclik.fh.unpak = @lik_loglogistic_unpak;
    reclik.fh.lp = @lik_loglogistic_lp;
    reclik.fh.lpg = @lik_loglogistic_lpg;
    reclik.fh.ll = @lik_loglogistic_ll;
    reclik.fh.llg = @lik_loglogistic_llg;
    reclik.fh.llg2 = @lik_loglogistic_llg2;
    reclik.fh.llg3 = @lik_loglogistic_llg3;
    reclik.fh.tiltedMoments = @lik_loglogistic_tiltedMoments;
    reclik.fh.invlink = @lik_loglogistic_invlink;
    reclik.fh.predy = @lik_loglogistic_predy;
    reclik.fh.predcdf=@lik_loglogistic_predcdf;
    reclik.fh.recappend = @lik_loglogistic_recappend;
    reclik.p=[];
    reclik.p.shape=[];
    if ~isempty(ri.p.shape)
      reclik.p.shape = ri.p.shape;
    end
  else
    % Append to the record
    reclik.shape(ri,:)=lik.shape;
    if ~isempty(lik.p)
      reclik.p.shape = lik.p.shape.fh.recappend(reclik.p.shape, ri, lik.p.shape);
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% NUMEL   Number of elements in an array or subscripted array expression.
%    N = NUMEL(A) returns the number of elements, N, in array A, equivalent
%    to PROD(SIZE(A)).
%
%    See also SIZE, PROD, SUBSREF.
%
%    Reference page in Doc Center
%       doc numel
%
%    Other functions named numel
%
%       categorical/numel      distributed/numel    tabular/numel
%       codistributed/numel    gpuArray/numel       tall/numel
%       Composite/numel        sym/numel            ts/numel
%       dataset/numel          symfun/numel
%

author: jmaih


text:
function [eventID] = ma_GetEventIDByEventNum(eventNum, script)
%getEventIDByEventNum Summary of this function goes here
%   Detailed explanation goes here
    try
        eventID = script{eventNum}.id;
    catch
        eventID = -1;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
close all
clear
clc

%% overall score
load('./AADBinfo.mat');
score = [testScore, trainScore];
figure;
hist(score, 10);
xlabel('score');
ylabel('#images');

%% attributes
attNameList = {'BalacingElements', 'ColorHarmony', 'Content', 'DoF', 'Light', 'MotionBlur', 'Object', 'Repetition', 'RuleOfThirds', 'Symmetry', 'VividColor' };
folderpath = './AADB_prepare';
attScore = cell(11,1);
for attID = 1:length(attNameList)
    fprintf('attribute-%s\n', attNameList{attID});
    attScore{attID} = [];
    filename = fullfile(folderpath, ['imgListTrainRegression_' attNameList{attID} '.txt']);
    numImage = numel(textread(filename,'%1c%*[^\n]'));
    scoreTMP = zeros(1,numImage);
    fn = fopen( filename, 'r' );
    for i = 1:numImage
        tline = fgets(fn);
        C = strsplit(tline, ' ');
%         imgName = C{1};
        imgLabel = str2double(C{2});
%         nameList{i} = imgName;
        scoreTMP(i) = imgLabel;
    end
    fclose(fn);
    attScore{attID} = [attScore{attID}, scoreTMP];

    filename = fullfile(folderpath, ['imgListTestRegression_' attNameList{attID} '.txt']);
    numImage = numel(textread(filename,'%1c%*[^\n]'));
    scoreTMP = zeros(1,numImage);
    fn = fopen( filename, 'r' );
    for i = 1:numImage
        tline = fgets(fn);
        C = strsplit(tline, ' ');
%         imgName = C{1};
        imgLabel = str2double(C{2});
%         nameList{i} = imgName;
        scoreTMP(i) = imgLabel;
    end
    fclose(fn);
    attScore{attID} = [attScore{attID}, scoreTMP];
end

%% show attribute distribution

A = cell2mat(attScore);
countA = zeros(11,3);
for attID = 1:length(attNameList)
    a = find(A(attID,:)==1);
    countA(attID, 1) = length(a);
    a = find(A(attID,:)==2);
    countA(attID, 2) = length(a);
    a = find(A(attID,:)==3);
    countA(attID, 3) = length(a);
end

figure;
h = bar(countA);
set(gca,'XTickLabel',{'BalacingElements'; 'ColorHarmony'; 'Content'; 'DoF'; 'Light'; 'MotionBlur'; 'Object'; 'Repetition'; 'RuleOfThirds'; 'Symmetry'; 'VividColor' })
l = {'negative', 'null', 'positive'};
legend(h,l);

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
%Q2TR	Convert unit-quaternion to homogeneous transform
%
%	T = q2tr(Q)
%
%	Return the rotational homogeneous transform corresponding to the unit
%	quaternion Q.
%
%	See also TR2Q

%	Copyright (C) 1993 Peter Corke
function t = q2tr(q)
	s = q(1);
	x = q(2);
	y = q(3);
	z = q(4);

	r = [	1-2*(y^2+z^2)	2*(x*y-s*z)	2*(x*z+s*y)
		2*(x*y+s*z)	1-2*(x^2+z^2)	2*(y*z-s*x)
		2*(x*z-s*y)	2*(y*z+s*x)	1-2*(x^2+y^2)	];
	t = eye(4,4);
	t(1:3,1:3) = r;
	t(4,4) = 1;

repo: ekfmonoslam


text:
function f = dubna(varargin)
% Construct the odf from the dubna data set as example for an SO3Fun.

f = mtexdata('dubnaODF');

if nargin > 0 && isa(varargin{1},'rotation')
  ori = varargin{1};
  f = f.eval(ori);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function gui_adv(x, y)
%
% This file is needed by FASTICAG

% This is the advanced options -dialog

% @(#)$Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables

% Handle to the window
global hf_FastICA_AdvOpt;

% Handles to some of the controls in window
global hpm_FastICA_finetune;
global he_FastICA_a1;
global he_FastICA_a2;
global he_FastICA_myy;
global he_FastICA_epsilon;
global he_FastICA_maxIterations;
global he_FastICA_maxFinetune;
global he_FastICA_sampleSize
global hpm_FastICA_initState;
global hb_FastICA_initGuess;
global ht_FastICA_initGuess;
global hpm_FastICA_displayMode;
global he_FastICA_displayInterval;
global hpm_FastICA_verbose;

% Some of the main variables needed
global g_FastICA_initGuess;
global g_FastICA_finetune;
global g_FastICA_a1;
global g_FastICA_a2;
global g_FastICA_myy;
global g_FastICA_epsilon;
global g_FastICA_maxNumIte;
global g_FastICA_maxFinetune;
global g_FastICA_initState;
global g_FastICA_sampleSize;
global g_FastICA_displayMo;
global g_FastICA_displayIn;
global g_FastICA_verbose;

global c_FastICA_appr_strD;
global c_FastICA_appr_strV;
global c_FastICA_finetune_strD;
global c_FastICA_finetune_strV;
global c_FastICA_iSta_strD;
global c_FastICA_iSta_strV;
global c_FastICA_dMod_strD;
global c_FastICA_dMod_strV;
global c_FastICA_verb_strD;
global c_FastICA_verb_strV;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Configuration options
FIGURENAME = 'FastICA: advanced options';
FIGURETAG = 'f_FastICAAdvOpt';
FIGURESIZE = [x y 450 280];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if this figure is already open - it should not!
% Can't have more than one copy - otherwise the global
% variables and handles can get mixed up.
if ~isempty(findobj('Tag',FIGURETAG))
  error('Error: advanced options dialog already open!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize some of the controls' values

% Did we already load some initial guess
pm_initState_Value = g_FastICA_initState;
if isempty(g_FastICA_initGuess) | (g_FastICA_initGuess == 1)
  t_initGuess_String = 'Not loaded';
else
  t_initGuess_String = 'Loaded';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the figure
a = figure('Color',[0.8 0.8 0.8], ...
	   'PaperType','a4letter', ...
	   'Name', FIGURENAME, ...
	   'NumberTitle', 'off', ...
	   'Tag', FIGURETAG, ...
	   'Position', FIGURESIZE, ...
	   'MenuBar', 'none');
set (a, 'Resize', 'off');

hf_FastICA_AdvOpt = a;

set(hf_FastICA_AdvOpt, 'HandleVisibility', 'callback');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From here on it get's ugly as I have not had time to clean it up

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the frames
pos_l=2;
pos_w=FIGURESIZE(3)-4;
pos_h=FIGURESIZE(4)-4;
pos_t=2;
h_f_adv_background = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_adv_background');

pos_w=120;
pos_l=FIGURESIZE(3)-(pos_w+2+2);
pos_h=FIGURESIZE(4)-2*4;
pos_t=4;
h_f_adv_side = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_adv_side');

pos_l=4;
pos_w=FIGURESIZE(3)-8-pos_w-2;
pos_h=FIGURESIZE(4)-8;
pos_t=4;
h_f_advopt = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_advopt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_advopt
bgc = get(h_f_advopt, 'BackgroundColor');

pos_w1=230;
pos_w2=70;

pos_frame=get(h_f_advopt, 'Position');
pos_h = 20;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - 6;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Fine-tune (g)', ...
  'Style','text', ...
  'Tag','t_727');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_finetune = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_finetune_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_finetune', ...
  'Value',g_FastICA_finetune);

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Parameter a1 (g = ''tanh'')', ...
  'Style','text', ...
  'Tag','t_22');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_a1 = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc Checka1', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_a1), ...
  'Style','edit', ...
  'Tag','e_a1');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Parameter a2 (g = ''gauss'')', ...
  'Style','text', ...
  'Tag','t_222');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_a2 = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc Checka2', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_a2), ...
  'Style','edit', ...
  'Tag','e_a2');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','mu (step size)', ...
  'Style','text', ...
  'Tag','t_223');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_myy = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc CheckMyy', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_myy), ...
  'Style','edit', ...
  'Tag','e_myy');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','epsilon (stopping criterion)', ...
  'Style','text', ...
  'Tag','t_23');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_epsilon = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_epsilon), ...
  'Style','edit', ...
  'Tag','e_epsilon');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Maximum number of iterations', ...
  'Style','text', ...
  'Tag','t_24');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_maxIterations = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_maxNumIte), ...
  'Style','edit', ...
  'Tag','e_maxIterations');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Maximum iterations in fine-tuning', ...
  'Style','text', ...
  'Tag','t_2412');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_maxFinetune = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_maxFinetune), ...
  'Style','edit', ...
  'Tag','e_maxFinetune');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Sample size (proportion)', ...
  'Style','text', ...
  'Tag','t_224');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_sampleSize = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc CheckSampleSize', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_sampleSize), ...
  'Style','edit', ...
  'Tag','e_sampleSize');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Initial state', ...
  'Style','text', ...
  'Tag','t_25');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_initState = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_iSta_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_initState', ...
  'Value',pm_initState_Value);

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
hb_FastICA_initGuess = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc loadGuess', ...
  'Position',[pos_l pos_t (pos_w1-60) pos_h], ...
  'String','Load initial guess', ...
  'UserData', g_FastICA_initGuess, ...
  'Tag','b_LoadGuess');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
ht_FastICA_initGuess = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[(pos_l) pos_t (pos_w2) pos_h], ...
  'String',t_initGuess_String, ...
  'Style','text', ...
  'Tag','t_initGuess');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Display mode', ...
  'Style','text', ...
  'Tag','t_27');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_displayMode = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_dMod_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_displayMode', ...
  'Value',g_FastICA_displayMo);

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Iterations between displays', ...
  'Style','text', ...
  'Tag','t_28');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_displayInterval = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_displayIn), ...
  'Style','edit', ...
  'Tag','e_displayInterval');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Verbose', ...
  'Style','text', ...
  'Tag','t_29');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_verbose = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_verb_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_verbose', ...
  'Value',g_FastICA_verbose);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_adv_side
pos_vspace = 6;
pos_hspace = 10;
pos_frame = get(h_f_adv_side, 'Position');
pos_w = 100;
pos_h = 30;
pos_l = pos_frame(1) + pos_hspace;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc OK', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','OK', ...
  'Tag','b_advOK');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Cancel', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Cancel', ...
  'Tag','b_advCancel');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Default', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Default', ...
  'Tag','b_advDefault');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Apply', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Apply', ...
  'Tag','b_advApply');

pos_t = pos_frame(2) + pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Help', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Help', ...
  'Tag','b_advHelp');


text:
%  FILE:   cnn_widerface.m
%
%    This function serves as the main function for training a model.
%
%  INPUT:  configuration     (see code for details)
%
%  OUTPUT: net               (trained network)
%          info              (training stats)

function [net, info] = cnn_widerface(varargin)
startup;

opts.keepDilatedZeros = false;
opts.inputSize = [500, 500];
opts.learningRate = 1e-4;

%% use customized training function ie. adam
opts.trainFn = '@cnn_train_dag';
opts.batchGetterFn = '@cnn_get_batch_logistic_zoom';
opts.freezeResNet = false;
opts.tag = '';
opts.clusterNum = 25;
opts.clusterName = '';
opts.bboxReg = true;
opts.skipLRMult = [1, 0.001, 0.0001, 0.00001];
opts.sampleSize = 256;
opts.posFraction = 0.5;
opts.posThresh = 0.7;
opts.negThresh = 0.3;
opts.border = [0, 0];
opts.pretrainModelPath = 'matconvnet/pascal-fcn8s-tvg-dag.mat';
opts.dataDir = fullfile('data','widerface') ;
opts.modelType = 'pascal-fcn8s-tvg-dag' ;
opts.networkType = 'dagnn' ;
opts.batchNormalization = true ;
opts.weightInitMethod = 'gaussian' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.minClusterSize = [2, 2];
opts.maxClusterSize = opts.inputSize;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.modelType ;
if opts.freezeResNet, sfx = ['freezed-' sfx] ; end
sfx = [sfx '-' 'sample' num2str(opts.sampleSize)] ;
sfx = [sfx '-' 'posfrac' num2str(opts.posFraction)] ;
sfx = [sfx '-' 'N' num2str(opts.clusterNum)];
if opts.bboxReg, sfx = [sfx '-' 'bboxreg']; end

if any(opts.inputSize~=500),
    sz = opts.inputSize;
    sfx = [sfx '-input' num2str(sz(1)) 'x' num2str(sz(2))];
end
if ~isempty(opts.clusterName)
    sfx = [sfx '-' 'cluster-' opts.clusterName];
end

opts.expDir = fullfile('models', ['widerface-' sfx]) ;
if ~isempty(opts.tag)
    opts.expDir = [opts.expDir '-' opts.tag];
end
[opts, varargin] = vl_argparse(opts, varargin) ;
opts.expDir

opts.batchSize = 10;
opts.numSubBatches = 1;
opts.numEpochs = 50;
opts.gpus = [1];
opts.numFetchThreads = 8;
opts.lite = false ;
opts.imdbPath = fullfile(opts.dataDir, 'imdb.mat');
opts = vl_argparse(opts, varargin) ;

opts.train = struct() ;
opts = vl_argparse(opts, varargin) ;

opts.train.gpus = opts.gpus;
opts.train.batchSize = opts.batchSize;
opts.train.numSubBatches = opts.numSubBatches;
opts.train.numEpochs = opts.numEpochs;
opts.train.learningRate = opts.learningRate;

opts.train.keepDilatedZeros = opts.keepDilatedZeros;

%% initialize model structure
fprintf('Trying to initialize the structure of %s\n', opts.modelType);
net = cnn_init('model', opts.modelType, ...
               'batchNormalization', opts.batchNormalization, ...
               'weightInitMethod', opts.weightInitMethod, ...
               'networkType', opts.networkType) ;

%% load pretrained weights
if ~isempty(opts.pretrainModelPath)
  fprintf('Loading pretrained weights from %s\n', opts.pretrainModelPath);
  if ~exist(opts.pretrainModelPath)
    [~,model_name,~] = fileparts(opts.pretrainModelPath);
    download_path = sprintf('./trained_models/%s.mat', model_name);
    download_url = sprintf('http://www.vlfeat.org/matconvnet/models/%s.mat', model_name);
    fprintf('Downloading model %s from %s to %s ... \n', model_name, download_url, download_path);
    cmd = sprintf('wget -O %s %s', download_path, download_url);
    system(cmd);
  end
  net = cnn_load_pretrain(net, opts.pretrainModelPath);
end

net.meta.inputSize = opts.inputSize;
net.meta.normalization.inputSize = opts.inputSize;
net.meta.normalization.border = opts.border;
net.meta.augmentation.transformation = 'none';
net.meta.augmentation.rgbVariance = [];

if ~exist(opts.expDir), mkdir(opts.expDir); end;

%% prepare data
if exist(opts.imdbPath)
    imdb = load(opts.imdbPath) ;
    fprintf('Loaded imdb from %s\n', opts.imdbPath);
else
    imdb = cnn_setup_imdb('dataDir', opts.dataDir);
    save(opts.imdbPath, '-struct', 'imdb') ;
    fprintf('Saved imdb to %s\n', opts.imdbPath);
end

%% save model options
optpath = fullfile(opts.expDir, 'opts.mat');
if ~exist(optpath), save(optpath, 'opts'); end

%% define batch getter function
%batchGetter = @cnn_get_batch_logistic_zoom;
if ~isempty(opts.batchGetterFn)
    batchGetter = str2func(opts.batchGetterFn);
end

%% compute image stats
imageStatsPath = fullfile(opts.dataDir, 'imageStats.mat') ;
if exist(imageStatsPath)
    load(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
else
    [averageImage, rgbMean, rgbCovariance] = getImageStats(batchGetter, ...
                                                      opts, net.meta, imdb) ;
    save(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
end
net.meta.augmentation.transformation = 'f5';
net.meta.normalization.averageImage = rgbMean ;
[v,d] = eig(rgbCovariance) ;
net.meta.augmentation.rgbVariance = 0.1*sqrt(d)*v' ;
clear v d ;

%% clustering
minh = opts.minClusterSize(1); minw = opts.minClusterSize(2);
maxh = opts.maxClusterSize(1); maxw = opts.maxClusterSize(2);
if isempty(opts.clusterName)
    tmp_str = ['%s_N%d_' sprintf('min%dx%d_max%dx%d.mat',minh,minw,maxh,maxw)];
else
    tmp_str = ['%s_N%d_' opts.clusterName '.mat'];
end
clusterName = @(name,nd)(sprintf(tmp_str, name, nd));

%clusterName = @(name,nd,minh,minw,maxh,maxw)(sprintf( ...
%    '%s_N%d_min%dx%d_max%dx%d.mat',name,nd,minh,minw,maxh,maxw));
%clusterPath = fullfile(opts.dataDir, clusterName(...
%    'RefBox',opts.clusterNum,minh,minw,maxh,maxw));
boxName = 'RefBox';
clusterPath = fullfile(opts.dataDir, clusterName(boxName,opts.clusterNum));

fprintf('cluster path: %s\n', clusterPath);

if ~exist(clusterPath)
  clusters = cluster_rects(imdb, opts.clusterNum, [minh minw], [maxh maxw]);
  save(clusterPath, 'clusters');
else
  load(clusterPath);
end
net.meta.clusters = clusters;

%% add predictors/losses
switch opts.modelType
  case 'resnet-101-simple'
    net = cnn_add_loss_fcn8s_resnet101_simple(opts, net);
  otherwise
    error(sprintf('Not Implemented: model type %s', opts.modelType));
end

%% compute receptive fields and canonical variable sizes
var2idx = containers.Map;
for i = 1:numel(net.vars)
    var2idx(net.vars(i).name) = i;
end
net.meta.var2idx = var2idx;

sz_ = opts.inputSize;
net.meta.varsizes = net.getVarSizes({'data',[sz_,3,1]});
net.meta.recfields = net.getVarReceptiveFields('data');

%% configure sampling hyperparameters
net.meta.sampleSize = opts.sampleSize;
net.meta.posFraction = opts.posFraction;
net.meta.posThresh = opts.posThresh;
net.meta.negThresh = opts.negThresh;

%% assign training function
trainFn = str2func(opts.trainFn);

%% print out training options
opts
opts.train

%% start training (no validation)
derOutputs = {'loss_cls', 1, 'loss_reg', 1};
[net, info] = trainFn(net, imdb, getBatchFn(batchGetter, opts, net.meta), ...
                      'expDir', opts.expDir,...
                      'derOutputs', derOutputs, ...
                      'val', nan,...
                      opts.train) ;

%% wrapper for batch getter function
function fn = getBatchFn(batchGetter, opts, meta)
useGpu = numel(opts.train.gpus) > 0 ;

bopts.numThreads = opts.numFetchThreads ;
bopts.inputSize = meta.normalization.inputSize ;
bopts.border = meta.normalization.border ;
bopts.averageImage = meta.normalization.averageImage ;
bopts.rgbVariance = meta.augmentation.rgbVariance ;
bopts.transformation = meta.augmentation.transformation ;

if isfield(meta, 'sampleSize'), bopts.sampleSize = meta.sampleSize; end;
if isfield(meta, 'posFraction'), bopts.posFraction = meta.posFraction; end;
if isfield(meta, 'posThresh'), bopts.posThresh = meta.posThresh; end;
if isfield(meta, 'negThresh'), bopts.negThresh = meta.negThresh; end;

if isfield(meta, 'clusters'), bopts.clusters = meta.clusters; end;
if isfield(meta, 'var2idx'), bopts.var2idx = meta.var2idx; end;
if isfield(meta, 'varsizes'), bopts.varsizes = meta.varsizes; end;
if isfield(meta, 'recfields'), bopts.recfields = meta.recfields; end;
switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(batchGetter, bopts,x,y) ;
  case 'dagnn'
    fn = @(x,y) getDagNNBatch(batchGetter, bopts,useGpu,x,y) ;
end

%% interface to batch getter function
function inputs = getDagNNBatch(batchGetter, opts, useGpu, imdb, batch)
imagePaths = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
imageSizes = imdb.images.size(batch, :);
%isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;
isVal = 0;
if isfield(opts, 'clusters') && ~isVal
    labelRects = imdb.labels.rects(batch);
else
    labelRects = [];
end

[images, clsmaps, regmaps] = batchGetter(imagePaths, imageSizes, labelRects, ...
                                       opts, 'prefetch', nargout == 0) ;

if nargout > 0
    % if we are training
    if useGpu && ~isempty(clsmaps) && ~isempty(regmaps)
        images = gpuArray(images) ;
    end
    inputs = {'data', images};

    if ~isempty(clsmaps)
        inputs(end+1:end+2) = {'label_cls', clsmaps};
    end
    if ~isempty(regmaps)
        inputs(end+1:end+2) = {'label_reg', regmaps};
    end
end

%% function to collect image statistics
function [averageImage, rgbMean, rgbCovariance] = ...
    getImageStats(batchGetter, opts, meta, imdb)
train = find(imdb.images.set == 1) ;
train = train(1: 101: end);
bs = 256 ;
opts.networkType = 'dagnn' ;
fn = getBatchFn(batchGetter, opts, meta) ;
avg = {}; rgbm1 = {}; rgbm2 = {};

for t=1:bs:numel(train)
    batch_time = tic ;
    batch = train(t:min(t+bs-1, numel(train))) ;
    fprintf('collecting image stats: batch starting with image %d ...', batch(1)) ;
    temp = fn(imdb, batch) ;
    temp = gather(temp{2});
    z = reshape(permute(temp,[3 1 2 4]),3,[]) ;
    n = size(z,2) ;
    avg{end+1} = mean(temp, 4) ;
    rgbm1{end+1} = sum(z,2)/n ;
    rgbm2{end+1} = z*z'/n ;
    batch_time = toc(batch_time) ;
    fprintf(' %.2f s (%.1f images/s)\n', batch_time, numel(batch)/ batch_time) ;
end
averageImage = mean(cat(4,avg{:}),4) ;
rgbm1 = mean(cat(2,rgbm1{:}),2) ;
rgbm2 = mean(cat(3,rgbm2{:}),3) ;
rgbMean = rgbm1 ;
rgbCovariance = rgbm2 - rgbm1*rgbm1' ;

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function sphere_grid_test10 ( )

%*****************************************************************************80
%
%% SPHERE_GRID_TEST10 tests SPHERE_CUBED_POINTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_GRID_TEST10\n' );
  fprintf ( 1, '  SPHERE_CUBED_POINTS computes points on a cubed sphere grid.\n' );

  n = 10;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of divisions on each face = %d\n', n );

  ns = sphere_cubed_point_num ( n );
  fprintf ( 1, '  Total number of points = %d\n', ns );

  xyz = sphere_cubed_points ( n, ns );

  r8mat_transpose_print_some ( 3, ns, xyz, 1, 1, 3, 20, ...
    '  Initial part of XYZ array:' );
%
%  Display the nodes on a sphere.
%
  figure ( )
  clf
  hold on
  [ x, y, z ] = sphere ( 20 );
  c = ones ( size ( z ) );
  surf ( x, y, z, c );
  plot3 ( xyz(1,:), xyz(2,:), xyz(3,:), 'b.', 'Markersize', 20 );
  axis equal
  grid on
  view ( 3 )
  xlabel ( '<--X-->' )
  ylabel ( '<--Y-->' )
  zlabel ( '<--Z-->' )
  title ( 'Cubed sphere points' );
  hold off
  filename = sprintf ( 'sphere_cubed_points_f%d.png', n );
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved to "%s".\n', filename );
%
%  Write the nodes to a file.
%
  if ( 1 )

    filename = sprintf ( 'sphere_cubed_grid_f%d.xyz', n );

    file_unit = fopen ( filename, 'wt' );
    for node = 1 : ns
      for i = 1 : 3
        fprintf ( file_unit, '  %f', xyz(i,node ) );
      end
      fprintf ( file_unit, '\n' );
    end
    fclose ( file_unit );

    fprintf ( 1, '  Wrote grid nodes to "%s".\n', filename );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Example 7.4: Binary hypothesis testing
% Figure 7.4
% Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Michael Grant, 2005-12-19

% Generate the data
P = [0.70  0.10
     0.20  0.10
     0.05  0.70
     0.05  0.10];
[n,m] = size(P);

% Construct the tradeoff curve by finding the
% the Pareto optimal deterministic detectors,
% which are the curve's vertices

nopts   = 1000;
weights = logspace(-5,5,nopts);
obj     = [0;1];
inds    = ones(n,1);

% minimize  -t1'*q1 - w*t2'*q2
% s.t.      t1+t2 = 1,  t1,t2 \geq 0

next = 2;
for i = 1 : nopts,
   PW = P * diag( [ 1 ; weights(i) ] );
   [ maxvals, maxinds ] = max( PW' );  % max elt in each row
   if (~isequal(maxinds', inds(:,next-1)))
       inds(:,next) = maxinds';
       T = zeros(m,n);
       for j=1:n
          T(maxinds(1,j),j) = 1;
       end;
       obj(:,next) = 1-diag(T*P);
       next = next+1;
   end;
end;
plot(obj(1,:), obj(2,:),[0 1], [0 1],'--');
grid on
for i=2:size(obj,2)-1
   text(obj(1,i),obj(2,i),['a', num2str(i-1)]);
end;

% Minimax detector: not deterministic

cvx_begin
    variables T( m, n ) D( m, m )
    minimize max( D(1,2), D(2,1) )
    subject to
        D == T * P;
        sum( T, 1 ) == 1;
        T >= 0;
cvx_end

objmp = 1 - diag( D );
text( objmp(1), objmp(2), 'b' );
xlabel('P_{fp}'); ylabel('P_{fn}');

%print -deps roc.eps

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function z=orthoAlt4Pres(P,P0,T0)
%%ORTHOALT4PRES Get the orthometric height (pressure altitude) for a
%               point up to 86km altitude using the U.S. Standard
%               Atmosphere 1976 model given a reading of the air pressure
%               at the point. The model is built into many altimeters.
%
%INPUTS: P A vector or matrix of the pressures in Pascals that are to be
%          converted into orthometric heights in meters.
%       P0 An optional parameter specifying the pressure at sea level at
%          the location in question in Pascals. This pressure might have
%          to be interpolated down below terrain given a measurement on
%          the terrain. If omitted, the value of
%          Constants.standardAtmosphericPressure is used for the pressure.
%          Except in certain instances, aircraft over 18,000ft altitude in
%          the United States have to have their altimeter calibrated to
%          the standard atmospheric pressure and not the true pressure on
%          the ground/ at sea level.
%       T0 An optional parameter specifying the temperature at sea level
%          in degrees Kelvin. The standard atmosphere only specifies one
%          value, which is in Constants.standardTemp and is the default
%          value used if this parameter is omitted.
%
%OUTPUTS: z The orthometric height of the point implied by the specified air
%           pressure.The standard does not specify a reference geoid/
%           permanent tide system.
%
%The atmospheric model is defined in [1]. and the full model goes up to
%1000km altitude. Up to an altitude of 32km, the U.S. standard atmosphere
%is identical to the International Civil Aviation Organization's (ICAO)
%standard atmosphere. Regulations related to how altimeters should be set
%with respect to a standard pressure are at
%https://www.faa.gov/air_traffic/publications/atpubs/aim/aim0702.html
%
%REFERENCES:
%[1] U.S. Standard Atmosphere, 1976, National Oceanic and Atmospheric
%   Administration Std., Oct. 1976.
%
%June 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    P0=Constants.standardAtmosphericPressure;%N/m
end

if(nargin<3)
    T0=Constants.standardTemp;
end

%%%Constants defined in the standard atmosphere.
%The radius of the spherical Earth used in the standard atmosphere. THis is
%not the same as the semi-major axis of the WGS-84 reference ellipsoid.
rE=6356766;%Meters
%The 1962 estimate of the UNiversal Gas Constant that is used in the
%standard.
R=8.31432;%J/(mole*K)
%The mean moleacular mass of air at sea level.
M0=28.9644e-3;%kg/mole
%Standard acceleration due to gravity.
g0=9.80665;%m/s

%Hb and LMb are from table 4 in the U.S. Standard Atmosphere 1976.
%Hb is the base geopotential height and LMb is the molecular scale
%temperature gradient from that base.
HbTable=[0;
    11000;
    20000;
    32000;
    47000;
    51000;
    71000;
    84852];
LMbTable=[-6.5e-3;
    0;
    1e-3;
    2.8e-3;
    0;
    -2.8e-3;
    -2e-3;
    0];

%TMb is computed using Equation 23 starting from T0. It isn't just
%tabulated in case the user wants to specify a non-standard T0.
TMbTable=zeros(8,1);
TMbTable(1)=T0;
for n=2:8
    TMbTable(n)=TMbTable(n-1)+LMbTable(n-1)*(HbTable(n)-HbTable(n-1));
end

numP=size(P(:));
z=zeros(size(P));
for curP=1:numP
    PCur=P(curP);

    %To determine the correct formula to use, we must first determine which
    %band of altitudes the pressure falls. Thus, we will first evaluate
    %pres4OrthoAlt for the pressure at all of the bands that separate, knowing
    %that the pressure monotonically decreases with altitude in the model.
    H=11000;
    zt=rE*(H/(rE-H));
    [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
    if(PCur<Ph)
        H=20000;
        zt=rE*(H/(rE-H));
        [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
        if(PCur<Ph)
            H=32000;
            zt=rE*(H/(rE-H));
            [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
            if(PCur<Ph)
                H=47000;
                zt=rE*(H/(rE-H));
                [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                if(PCur<Ph)
                    H=51000;
                    zt=rE*(H/(rE-H));
                    [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                    if(PCur<Ph)
                        H=71000;
                        zt=rE*(H/(rE-H));
                        [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                        if(PCur<Ph)
                            H=72000;
                            zt=rE*(H/(rE-H));
                            [~,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                        end
                    end
                end
            end
        end
    end

    %Solve the inverse of the appropriate air pressure equation to determine
    %the geopotential height.
    if(LMb==0)
        %The inverse of the exponential expression for air pressure as a
        %function of altitude in Equation 33b
        H=Hb-R*TMb*log(PCur/Pb)/(g0*M0);
    else
        %The inverse of air pressure according to equation 33a, which is the
        %linear formulation.
        H=Hb+(TMb/LMb)*((PCur/Pb)^(-LMb*R/(g0*M0))-1);
    end

    %Convert the geopotential height into an orthometric height using the
    %approximation in the standard.
    z(curP)=rE*(H/(rE-H));
end

end

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


text:
classdef nnbnorm < nntest
  properties (TestParameter)
    rows = {2 8 13}
    cols = {2 8 17}
    numDims = {1 3 4}
    batchSize = {2 7}
  end
  methods (Test)
    function regression(test, rows,cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs) ;
      g = test.randn(nd, 1) / test.range ;
      b = test.randn(nd, 1) / test.range ;
      epsilon = 0.001 ;

      [y,m] = vl_nnbnorm(x,g,b,'epsilon',epsilon) ;
      n = numel(x) / size(x,3) ;
      mu = sum(sum(sum(x,1),2),4) / n ;
      sigma2 = sum(sum(sum(bsxfun(@minus,x,mu).^2,1),2),4) / n + epsilon ;
      sigma = sqrt(sigma2) ;
      m_ = [mu(:),sigma(:)] ;
      test.eq(m,m_) ;

      g = reshape(g,1,1,[]) ;
      b = reshape(b,1,1,[]) ;
      y_ = bsxfun(@plus,b,bsxfun(@times,g,bsxfun(@rdivide,bsxfun(@minus,x,mu),sigma))) ;
      test.eq(y,y_) ;
    end

    function basic(test, rows, cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs) ;
      %g = test.randn(1, 1, nd, 1) ;
      %b = test.randn(1, 1, nd, 1) ;
      g = test.randn(nd, 1) / test.range ;
      b = test.randn(nd, 1) / test.range ;

      y = vl_nnbnorm(x,g,b) ;
      dzdy = test.randn(size(y)) ;
      [dzdx,dzdg,dzdb] = vl_nnbnorm(x,g,b,dzdy) ;

      test.der(@(x) vl_nnbnorm(x,g,b), x, dzdy, dzdx, test.range * 1e-3) ;
      test.der(@(g) vl_nnbnorm(x,g,b), g, dzdy, dzdg, 1e-2) ;
      test.der(@(b) vl_nnbnorm(x,g,b), b, dzdy, dzdb, 1e-3) ;
    end

    function givenMoments(test, rows, cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs) ;
      %g = test.randn(1, 1, nd, 1) ;
      %b = test.randn(1, 1, nd, 1) ;
      g = test.randn(nd, 1) / test.range ;
      b = test.randn(nd, 1) / test.range ;

      [y,m] = vl_nnbnorm(x,g,b) ;
      [y_,m_] = vl_nnbnorm(x,g,b,'moments',m) ;

      test.eq(y,y_) ;
      test.eq(m,m_) ;

      dzdy = test.randn(size(y)) ;
      [dzdx,dzdg,dzdb,m__] = vl_nnbnorm(x,g,b,dzdy) ;
      [dzdx_,dzdg_,dzdb_,m___] = vl_nnbnorm(x,g,b,dzdy,'moments',m) ;

      test.eq(m,m__)
      test.eq(m,m___)
      test.eq(dzdx,dzdx_) ;
      test.eq(dzdg,dzdg_) ;
      test.eq(dzdb,dzdb_) ;
    end

  end
end

repo: UDT


text:
function t = uminus(t)
%UMINUS Unary minus for ttensor.
%
%   See also TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

t.core = -t.core;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function test085 ( dim_num, n, z )

%*****************************************************************************80
%
%% TEST085 tests R0_MEASURE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST085\n' );
  fprintf ( 1, '  R0_MEASURE computes the R0 measure of quality.\n' );
  fprintf ( 1, '  The Riesz S = 0 energy measure  R0 = %14f\n', ...
    r0_measure ( dim_num, n, z ) );

  return
end


text:
% read_erplabevent() - import ERPLAB dataset events
%
% Usage:
%     >> event = read_erplabevent(filename, ...);
%
% Inputs:
%     filename - [string] file name
%
% Optional inputs:
%     'header' - FILEIO structure header
%
% Outputs:
%     event - FILEIO toolbox event structure
%
% Modified from read_eeglabevent

%123456789012345678901234567890123456789012345678901234567890123456789012
%
% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.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 event = read_erplabevent(filename, varargin)

if nargin < 1
  help read_erplabheader;
  return;
end;

hdr = ft_getopt(varargin, 'header');

if isempty(hdr)
  hdr = read_erplabheader(filename);
end

event    = [];                % these will be the output in FieldTrip format
oldevent = hdr.orig.bindescr;    % these are in ERPLAB format

for index = 1:length(oldevent)
    event(end+1).type     = 'trial';
    event(end  ).sample   = (index-1)*hdr.nSamples + 1;
    event(end  ).value    = oldevent{index};
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;
end;

author: spm


text:
function msm_to_mm_test16 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST16 tests MSM_TO_MM_COORDINATE_INTEGER_SKEW_SYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST16\n' );
  fprintf ( 1, '  Convert an MSM to MM coordinate integer skew-symmetric format.\n' );

  output_filename = 'msm_to_mm_test16.mm';

  a = i4mat_indicator ( 4, 4 );
  a = a - a';
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'coordinate', 'integer', 'skew-symmetric' );

  return
end


text:
function ellipsoid_monte_carlo_test03 ( )

%*****************************************************************************80
%
%% ELLIPSOID_MONTE_CARLO_TEST03 uses ELLIPSOID_SAMPLE on a 3D ellipse centered at (1,2,3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;

  a = [ ...
    9.0, 6.0, 3.0; ...
    6.0, 5.0, 4.0; ...
    3.0, 4.0, 9.0 ]';
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    0, 0, 3 ]';
  r = 0.5;
  v = [ 1.0, 2.0, 3.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSOID_MONTE_CARLO_TEST03\n' );
  fprintf ( 1, '  Use ELLIPSOID_SAMPLE to estimate integrals\n' );
  fprintf ( 1, '  in a 3D ellipse (x-v)'' * A * (x-v) <= r^2.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ellipsoid radius R = %g\n', r );
  r8vec_print ( m, v, '  Ellipsoid center V:' )
  r8mat_print ( m, m, a, '  Ellipsoid matrix A:' )

  volume = ellipsoid_volume ( m, a, v, r );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ellipsoid volume = %g\n', volume );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N        1              X               Y  ' );
  fprintf ( 1, '              Z                X^2            YZ              Z^3\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = ellipsoid_sample ( m, n, a, v, r, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e(1:m) = e_test(1:m,j);

      value = monomial_value ( m, n, e, x );

      result = volume * sum ( value(1:n) ) / n;

      fprintf ( 1, '  %14.6g', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  return
end


text:
function Clg=cct2clg(Cct,lat,lon)
% CCT2CLG  Convert CT covariance matrix to local geodetic.
%   Vectorized. Requires rotct2lg. See also CLG2CCT.
% Version: 2011-04-05
% Useage:  Clg=cct2clg(Cct,lat,lon)
% Input:   Cct - CT covariance matrix
%          lat - vector of station latitudes (rad)
%          lon - vector of station longitudes (rad)
% Output:  Clg - LG covariance matrix

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

n=length(lat);
if (n*3 ~= max(size(Cct)) )
  error('Size of lat,lon does not match size of Cct');
end

for i=1:n
  sinlat=sin(lat(i));
  coslat=cos(lat(i));
  sinlon=sin(lon(i));
  coslon=cos(lon(i));
  Ji=rotct2lg(lat(i),lon(i));
  indi=(i-1)*3+[1:3];
  for j=1:n
    sinlat=sin(lat(j));
    coslat=cos(lat(j));
    sinlon=sin(lon(j));
    coslon=cos(lon(j));
    Jj=rotct2lg(lat(j),lon(j));
    indj=(j-1)*3+[1:3];
    Clg(indi,indj)=Ji*Cct(indi,indj)*Jj';
  end
end

author: Sable

repo: mcbench-benchmarks