text:
function varargout = ssofm(varargin)
% SSOFM M-file for ssofm.fig
%      SSOFM, by itself, creates a new SSOFM or raises the existing
%      singleton*.
%
%      H = SSOFM returns the handle to a new SSOFM or the handle to
%      the existing singleton*.
%
%      SSOFM('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SSOFM.M with the given input arguments.
%
%      SSOFM('Property','Value',...) creates a new SSOFM or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ssofm_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ssofm_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help ssofm

% Last Modified by GUIDE v2.5 10-Jun-2006 02:13:27

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ssofm_OpeningFcn, ...
                   'gui_OutputFcn',  @ssofm_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ssofm is made visible.
function ssofm_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ssofm (see VARARGIN)

% Choose default command line output for ssofm
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes ssofm wait for user response (see UIRESUME)
% uiwait(handles.ssofmFigure);

% --- Outputs from this function are returned to the command line.
function varargout = ssofm_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --------------------------------------------------------------------
function exit_Callback(hObject, eventdata, handles)
% hObject    handle to exit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close;

% --------------------------------------------------------------------
function file_Callback(hObject, eventdata, handles)
% hObject    handle to file (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function loaddata_Callback(hObject, eventdata, handles)
% hObject    handle to load (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filenameData, pathnameData]=uigetfile('*.mat','Load data...','MultiSelect', 'off');
set(handles.data2Text,'String',filenameData);
handles.Data = [pathnameData filenameData];
guidata(ssofm,handles)
a = get(handles.data2Text,'String');
b = get(handles.ssofm2Text,'String');
if a(1:3)~='...' & b(1:3)~='...'
    set(handles.trainButton,'Enable','on')
end

function epochsEdit_Callback(hObject, eventdata, handles)
% hObject    handle to epochsEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of epochsEdit as text
%        str2double(get(hObject,'String')) returns contents of epochsEdit as a double

% --- Executes during object creation, after setting all properties.
function epochsEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to epochsEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on button press in trainButton.
function trainButton_Callback(hObject, eventdata, handles)
% hObject    handle to trainButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(hObject,'Enable','off');
set(handles.plotglyphButton,'Enable','off');
set(handles.retrainButton,'Enable','off');
s = str2num(get(handles.sizeEdit,'String'));
epochs = str2num(get(handles.epochsEdit,'String'));
if get(handles.plotRadio,'Value') == 1
    showglyph = 'plot';
end
if get(handles.glyphRadio,'Value') == 1
    showglyph = 'glyph';
end
if get(handles.aviRadio,'Value') == 1
    showglyph = get(handles.aviText,'String');
    if isempty(showglyph)
        showglyph = 'test.avi';
    else
        showglyph = [get(handles.aviText,'String') '.avi'];
    end
end
load(handles.Data);
load(handles.Ssofm);
[w,g,r,c,freq]=trainGlyph(P,X,C,s,epochs,showglyph);
handles.w = w;
handles.g = g;
handles.r = r;
handles.c = c;
handles.freq = freq;
handles.X = X;
guidata(ssofm,handles);
set(hObject,'Enable','on');
set(handles.plotglyphButton,'Enable','on');
set(handles.retrainButton,'Enable','on');

function sizeEdit_Callback(hObject, eventdata, handles)
% hObject    handle to sizeEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of sizeEdit as text
%        str2double(get(hObject,'String')) returns contents of sizeEdit as a double

% --- Executes during object creation, after setting all properties.
function sizeEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sizeEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --------------------------------------------------------------------
function loadssofm_Callback(hObject, eventdata, handles)
% hObject    handle to loadssofm (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filenameSsofm, pathnameSsofm]=uigetfile('*.mat','Load S-SOFM...','MultiSelect', 'off');
set(handles.ssofm2Text,'String',filenameSsofm);
handles.Ssofm = [pathnameSsofm filenameSsofm];
guidata(ssofm,handles);
a = get(handles.data2Text,'String');
b = get(handles.ssofm2Text,'String');
if a(1:3)~='...' & b(1:3)~='...'
    set(handles.trainButton,'Enable','on')
end

% --------------------------------------------------------------------
function about_Callback(hObject, eventdata, handles)
% hObject    handle to about (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
ssofmAbout;

% --------------------------------------------------------------------
function help_Callback(hObject, eventdata, handles)
% hObject    handle to help (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function aviText_Callback(hObject, eventdata, handles)
% hObject    handle to aviText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of aviText as text
%        str2double(get(hObject,'String')) returns contents of aviText as a double

% --- Executes during object creation, after setting all properties.
function aviText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to aviText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on button press in plotglyphButton.
function plotglyphButton_Callback(hObject, eventdata, handles)
% hObject    handle to plotglyphButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure;
glyph(handles.X,handles.r);

% --- Executes on button press in aviRadio.
function aviRadio_Callback(hObject, eventdata, handles)
% hObject    handle to aviRadio (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of aviRadio
if get(hObject,'Value')==1
    set(handles.aviText,'Enable','on');
else
    set(handles.aviText,'Enable','off');
end

% --- Executes on button press in retrainButton.
function retrainButton_Callback(hObject, eventdata, handles)
% hObject    handle to retrainButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

w = handles.w;
freq = handles.freq;
X = handles.X;
load(handles.Data);
load(handles.Ssofm);
s = str2num(get(handles.sizeEdit,'String'));
if get(handles.plotRadio,'Value') == 1
    showglyph = 'plot';
end
if get(handles.glyphRadio,'Value') == 1
    showglyph = 'glyph';
end
if get(handles.aviRadio,'Value') == 1
    showglyph = get(handles.aviText,'String');
    if isempty(showglyph)
        showglyph = 'test.avi';
    else
        showglyph = [get(handles.aviText,'String') '.avi'];
    end
end

j=1;
% For each epoch
epochs = str2num(get(handles.epochsEdit,'String'));
for i=1:max(epochs)
    % Define the learning rate
    a=0.1*exp(-i/max(epochs));
    % Update the weights
    [w,freq]=updateGlyph(P,X,C,w,a,s,freq);
    % Update g
    g(i+1,:)=[mean(stdrc(w)) 1.96*std(stdrc(w))/sqrt(size(X,1))];
    if strcmp(showglyph,{'none','plot'})==0
        % Standardize the range
        r=stdrc(w);
        % Draw the glyph
        glyph(X,r);
        drawnow
        if strcmp(showglyph,'glyph')==0
            % Add another frame
            aviobj = addframe(aviobj,getframe(gcf));
        end
    elseif showglyph=='plot'
        %plot(0:i,g(1:i+1),'.-','LineWidth',2);
        errorbar(0:i,g(1:i+1,1),g(1:i+1,2),'.-','LineWidth',2);
        axis([0 max(1,i) 1 1.2])
        xlabel('epoch #');
        ylabel('convergence');
        if i<=10
            if i==1
                set(gca,'XTick',0:2);
            else
                set(gca,'XTick',1:i);
            end
        else
            set(gca,'XTickMode','auto');
        end
        drawnow
    else
        % Show progress
        if i==epochs(j)
            fprintf('Epoch %4.0f, log10(convergence) %9.4f\n',i,log10(g(i)));
            j=j+1;
        end
    end
end

[r,c]=stdrc(w);

if strcmp(showglyph,{'none','plot','glyph'})==0
    % Close the avi file
    aviobj = close(aviobj);
end

handles.w = w;
handles.freq = freq;
handles.X = X;


text:
function[image] = labelMapToColorImage(labelMap, labelCount, cmap)
% [image] = labelMapToColorImage(labelMap, [labelCount], [cmap])
%
% Convert an mxn labelMap to a color image.
% Each index is replaced by a certain color from a colormap.
% If labelCount is specified, we use absolute colors, otherwise the colors
% are sorted by label indices.
%
% Copyright by Holger Caesar, 2014

% Default arguments
if ~exist('labelCount', 'var'),
    labelCount = [];
end;
if ~exist('cmap', 'var'),
    cmap = @jet;
end;

% Get a unique list of all labels
labelList = unique(labelMap(:));
labelList(labelList == 0) = [];
labelListCount = numel(labelList);

% Initialize result
image = zeros(size(labelMap, 1), size(labelMap, 2), 3);
if isempty(labelCount),
    % Take a variable color scheme
    colormap = cmap(labelListCount);
else
    % Take a fixed color scheme relative to the number of labels in the
    % dataset
    colormap = cmap(labelCount);
    colormap = colormap(labelList, :);
end;

% Go through each label and replace its pixels by a color
for labelListIdx = 1 : labelListCount,
    labelIdx = labelList(labelListIdx);
    indices = labelMap == labelIdx;

    curImage = cat(3, ...
        indices .* colormap(labelListIdx, 1), ...
        indices .* colormap(labelListIdx, 2), ...
        indices .* colormap(labelListIdx, 3));
    image = image + curImage;
end;

author: nightrome

repo: matconvnet-calvin


text:
function metID = findMetIDs(model, metList)
% Finds metabolite numbers in a model
%
% USAGE:
%
%    metID = findMetIds(model, metList)
%
% INPUTS:
%    model:      COBRA model structure
%    metList:    List of metabolites
%
% OUTPUT:
%    metID:      List of metabolite IDs corresponding to `metList`
%
% .. Author: - Jan Schellenberger 8/15/08

if (iscell(metList))
    [tmp,metID] = ismember(metList,model.mets);
else
    metID = find(strcmp(model.mets,metList));
    if (isempty(metID))
        metID = 0;
    end
    if (length(metID) > 1)
        metID = metID(1);
    end
end

author: opencobra

repo: cobratoolbox


text:
function [CMerged,omega,xMerged]=covarianceIntersect(C1,C2,optCrit,x1,x2)
%%COVARIANCEINTERSECT Perform covariance intersection. This is a method of
%                     fusing the first two moments of estimates when the
%                     correlation between the estimates is unknown. Given
%                     the covariance matrices of the estimates, this
%                     function returns a covariance matrix and a scaling
%                     factor, which can be used to fuse estimates.
%                     Alternatively, if the estimates themselves are given,
%                     this function can also return the merged estimate.
%                     Compare this function to ellipsoidIntersect.
%
%INPUTS: C1 The nXn positive definite covariance matrix of the first
%           estimate.
%        C2 The nXn positive definite covariance matrix of the second
%           estimate.
%   optCrit An optional  parameter specifying the optimization criterion
%           for the fusion. Possible values are
%           'det' (The default if omitted or an empty matrix is passed) The
%                 fusion is performed to optimize the determinant of the
%                 matrix on the output.
%           'trace' The fusion is performed to minimize the trace of the
%                 output matrix.
%    x1,x2 The optional nX1 estimate vectors to be merged. These are only
%          used if xMerged is requested on the output.
%
%OUTPUTS: CMerged The nXn fused covariance matrix.
%           omega The weighting that played a role in the fusion of the
%                 covariance matrix/ state estimates. Specifically,
%                 CMerged=inv(omega*inv(C1)+(1-omega)*inv(C2))
%                 and the merged estimates are
%                 xMerged=CMerged*(omega*inv(C1)*x1+(1-omega)*inv(C2)*x2)
%         xMerged The merged estimate. This requires x1,x2 to be given on
%                 the input. The covariance matrix associated with the
%                 merged estimate is CMerged.
%
%The first (non-dissertation) publication of covariance intersection for
%fusing measurements having unknown correlations is [2]. Various optimality
%criteria are derived in [3]. However, the most practical approach is that
%of [1], which clearly expresses the solution in terms of the solution of
%polynomials, with some solutions given in closed form. Thus, this function
%uses the closed form solutions, when available, and the general polynomial
%solutions otherwise.
%
%Note that covariance intersection is overly convervative in its covariance
%estimates as mentioned in Chapter 9.3.7 of [5].
%
%EXAMPLE:
%Here we use the numerical values form the example for ellipsoidal
%intersection in [1]. The plots in the paper appear to be incorrect.
% xi=[1;-2];
% Pi=[3,0;0,0.4];
% xj=[-2;-1];
% Pj=[2,-0.8;-0.8,1];
%
% [P,~,x]=covarianceIntersect(Pi,Pj,'det',xi,xj);
% figure()
% clf
% hold on
% drawEllipse(xi,inv(Pi),[],'--r')
% drawEllipse(xj,inv(Pj),[],'--g')
% drawEllipse(x,inv(P),[],'-b')
%
%REFERENCES:
%[1] M. Reinhardt, B. Noack, and U. D. Hanebeck, "Closed-form optimization
%    of covariance intersection for low-dimensional matrices," in
%    Proceedings of the 15th International Conference on Information
%    Fusion, Singapore, 9-12 Jun. 2012, pp. 1891-1896.
%[2] S. J. Julier and J. K. Uhlmann, "A non-divergent estimation algorithm
%    in the presence of unknown correlations," in Proceedings of the
%    American Control Conference, vol. 4, Albuquerque, NM, 4-6 Jun. 1997,
%    pp. 2369-2373.
%[3] L. Chen, P. O. Arambel, and R. K. Mehra, "Estimation under unknown
%    correlation: Covariance intersection revisited," IEEE Transactions on
%    Automatic Control, vol. 47, no. 11, pp. 1879-1882, Nov. 2002.
%[4] J. Sijs, M. Lazar, and P. Bosch, "State fusion with unknown
%    correlation: Ellipsoidal intersection," in Proceedings of the 2010
%    American Control Conference, Baltimore, MD, 30 Jun. - 2 Jul. 2010.
%[5] Y. Bar-Shalom, P. K. Willett, and X. Tian, Tracking and Data Fusion.
%    Storrs, CT: YBS Publishing, 2011.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(optCrit))
    optCrit='det';
end

n=length(C1);

%The transformation matrix in Equation 2 of [1], but here we use an
%SVD-based algorithm instead of an eigenvalue-based algorithm, because it
%is numerically more stable.
A=twoMatDiag(C1,C2);

%The positive scalars d.
d=diag(A*C2*A');
dBar=1./d;
dTilde=dBar./(1-dBar);

invC1=inv(C1);
invC2=inv(C2);

%Determinant minimization
switch(optCrit)
    case 'det'
        if(n==1)
            if(C1<C2)
                omega=1;
            else
                omega=0;
            end
        elseif(n==2)
            %The solution is Equation 11 in [1].
            omega=-(1/2)*(dTilde(1)+dTilde(2));
        elseif(n==3)
            sqrtTerm=sqrt(sum(dTilde.^2)-dTilde(1)*dTilde(2)-dTilde(1)*dTilde(3)-dTilde(2)*dTilde(3));
            %Equation 12 in [1].
            omega(1,1)=-(1/3)*(sum(dTilde)+sqrtTerm);
            omega(2,1)=-(1/3)*(sum(dTilde)-sqrtTerm);
        elseif(n==4)
            b=(3/4)*sum(dTilde);

            c=0;
            for i=1:n
               for j=(i+1):n
                   c=c+dTilde(i)*dTilde(j);
               end
            end
            c=(1/2)*c;

            d=0;
            for i=1:n
                for j=(i+1):n
                    for k=(j+1):n
                        d=d+dTilde(i)*dTilde(j)*dTilde(k);
                    end
                end
            end
            d=(1/4)*d;

            Q=sqrt((2*b^3-9*b*c+27*d)^2-4*(b^2-3*c)^3);
            C=((1/2)*(Q+2*b^3-9*b*c+27*d))^(1/3);

            %Equation 20
            omega(1,1)=-(b/3)-C/3-(b^2-3*c)/(3*C);
            omega(2,1)=-(b/3)+C*(1+1j*sqrt(3))/6+(1-1j*sqrt(3))*(b^2-3*c)/(6*C);
            omega(3,1)=-(b/3)+C*(1-1j*sqrt(3))/6+(1+1j*sqrt(3))*(b^2-3*c)/(6*C);
        else%For n>4
            %The pi values are given in the unnumbered equation after
            %Equation 10 in [1].
            piVals=zeros(1,n-1);
            for curPi=1:(n-1)
                piIdx=1:curPi;
                while(1)
                    %Add the current  product
                    piVals(curPi)=piVals(curPi)+prod(dTilde(piIdx));

                    %Increment the value in the innermost sum.
                    curLevel=curPi;
                    allDone=false;
                    while(piIdx(curLevel)+1>(n-curPi+curLevel))
                        curLevel=curLevel-1;
                        if(curLevel<1)
                            allDone=true;
                            break;
                        end
                    end
                    if(allDone)
                        break;
                    end
                    %Increment the index
                    piIdx(curLevel)=piIdx(curLevel)+1;

                    %Fill in the minimum values in further nested sums.
                    for k=(curLevel+1):curPi
                        piIdx(k)=piIdx(k-1)+1;
                    end
                end
            end

            %Given the pi values, build the polynomial in omega.
            poly2Solve=(n:-1:1).*[1,piVals];

            %Omega hypotheses are the roots of an order n-1 polynomial in
            %Equation 10 of [1].
            omega=roots(poly2Solve);
        end

        if(n>1)
        %Given a set of omega values, we have to choose the best solution.
        %We also have to include omega=0 and omega=1 as possibilities as
        %per Algorithm 1 in [1]. Some of the solutions might be slightly
        %imaginary due to finite precision errors. We will just discard the
        %imaginary parts. This does not affect the choice of the optimal
        %solution.
            omega=[real(omega);0;1];
            numHyp=length(omega);

            minCost=Inf;
            minOmega=[];
            for curHyp=1:numHyp
                if(omega(curHyp)>=0&&omega(curHyp)<=1)
                    %The determinant of a matrix inverse is the inverse of
                    %the determinant.
                    cost=1/det(omega(curHyp)*invC1+(1-omega(curHyp))*invC2);

                    if(cost<minCost)
                       minCost=cost;
                       minOmega=omega(curHyp);
                    end
                end
            end

            omega=minOmega;
        end
    case 'trace'
        if(n==1)
            %The solution is the same as for the determinant minimization
            %criterion.
            if(C1<C2)
                omega=1;
            else
                omega=0;
            end
        elseif(n==2)
            invA=inv(A);
            a=diag(invA'*invA);

            denom=a(1)*(1+dTilde(1))+a(2)*(1+dTilde(2));
            p=(a(1)*dTilde(2)*(1+dTilde(1))+a(2)*dTilde(1)*(1+dTilde(2)))/denom;
            q=(a(1)*dTilde(2)^2*(1+dTilde(1))+a(2)*dTilde(1)^2*(1+dTilde(2)))/denom;

            %Equation 19
            omega(1,1)=-p+sqrt(p^2-q);
            omega(2,1)=-p-sqrt(p^2-q);
        else%for n>=3
            %Use convolutions to build up the polynomial in Equation 17 in
            %[1].
            invA=inv(A);
            a=diag(invA'*invA);

            %The polynomial is degree 2*(n-1).
            numPolyEls=2*(n-1)+1;
            poly2Solve=zeros(numPolyEls,1);
            for i=1:n
                coeff=a(i)*(1+dTilde(i));

                jPoly=zeros(numPolyEls,1);
                jPoly(end)=1;
                for j=1:n
                    if(j==i)
                        continue;
                    end

                    jPoly=conv(jPoly,[1;2*dTilde(j);dTilde(j)^2]);
                    %Get rid of zero-padding at the beginning.
                    jPoly=jPoly((end-numPolyEls+1):end);
                end
                poly2Solve=poly2Solve+coeff*jPoly;
            end

            %omega hypotheses are the roots of the equation.
            omega=roots(poly2Solve);
        end

        %Given a set of omega values, we have to choose the best solution.
        %We also have to include omega=0 and omega=1 as possibilities as
        %per Algorithm 1 in [1]. Some of the solutions might be slightly
        %imaginary due to finite precision errors. We will just discard the
        %imaginary parts. This does not affect the choice of the optimal
        %solution.
        if(n>1)
            omega=[real(omega);0;1];
            numHyp=length(omega);

            minCost=Inf;
            minOmega=[];
            for curHyp=1:numHyp
                if(omega(curHyp)>=0&&omega(curHyp)<=1)
                    cost=trace(inv(omega(curHyp)*invC1+(1-omega(curHyp))*invC2));

                    if(cost<minCost)
                       minCost=cost;
                       minOmega=omega(curHyp);
                    end
                end
            end

            omega=minOmega;
        end
    otherwise
        error('Invalid optimality criterion chosen.')
end

CMergedInv=omega*invC1+(1-omega)*invC2;
CMerged=inv(CMergedInv);

if(nargout>2)
    xMerged=CMergedInv\(omega*invC1*x1+(1-omega)*invC2*x2);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [A,b,c,lb,ub] = ls_preprocess(A,b,c,lb,ub,opts)
% PREPROCESS  - Preprocessing input data.
% Usage:  [A,b,c,lb,ub,probData] = preprocess(A,b,c,lb,ub,probData)

% Yin Zhang, last updated April, 1995
% Department of Mathematics and Statistics
% University of Maryland  Baltimore County
% Modified J.Currie AUT May 2013

global probData

verb = opts.verb;

if(verb), fprintf('LIPSOL Preprocessing ...\n'); end

%Ensure columns
if(size(b,2) > 1), b = b'; end
if(size(c,2) > 1), c = c'; end
if(size(lb,2) > 1), lb = lb'; end
if(size(ub,2) > 1), ub = ub'; end

%Check problem sizes
[m,n] = size(A);
if(m ~= size(b,1)), error('A does not have the same number of rows as b'); end
if(n ~= size(c,1)), error('A does not have the same number of colums as c'); end
if(n ~= size(lb,1)), error('lb does not have the same number of elements as c'); end
if(n ~= size(ub,1)), error('ub does not have the same number of elements as c'); end

if any(lb > ub)
   if(verb), fprintf('\nPreprocessor: Lower bound exceeds upper bound\n'); end
   probData.isFeasible = 0;
   return;
end

%Convert Infinite Bounds to "BIG"
ub(isinf(ub)) = opts.big;
%Check for infinite lower bounds (should be allowed, but no goood...)
if(any(isinf(lb))), throwAsCaller(MException('OPTI:LIPSOLINFLB','LIPSOL only solves problems with Finite Lower Bounds')); end

if(~issparse(A)), A = sparse(A); end;
b = sparse(b); c = sparse(c);
lb = sparse(lb);

%----- delete fixed variables -----
fixed = lb == ub;
probData.Fixed_exist = any(fixed);
if (probData.Fixed_exist)
   ifix = find(fixed);
   infx = find(1 - fixed);
   xfix = lb(ifix);
   c    = c(infx);
   b    = b - A(:,ifix)*sparse(xfix);
   A    = A(:,infx);
   lb = lb(infx);
   ub = ub(infx);
   if(verb), fprintf(' - %i fixed vars\n', length(ifix)); end
   probData.data_changed = 1;
   probData.ifix = ifix;
   probData.infx = infx;
   probData.xfix = xfix;
end

%----- delete zero rows -----
rnnzct = sum(spones(A'));
if any(rnnzct == 0)
   izrows = find(rnnzct == 0);
   if any(b(izrows) ~= 0)
      if(verb), fprintf('\nPreprocessor: problem infeasible\n'); end
      probData.isFeasible = 0;
      return;
   end;
   inzrows = find(rnnzct > 0);
   A = A(inzrows,:); b = b(inzrows);
   rnnzct = rnnzct(inzrows);
   if(verb), fprintf(' - %i 0-rows', length(izrows)); end
   probData.data_changed = 1; m = size(A,1);
end

%----- make A structurally "full rank" -----
sprk = sprank(A');
if (sprk < m)
   [dmp, ~] = dmperm(A);
   irow = dmp(1:sprk);
   A = A(irow,:); b = b(irow);
   rnnzct = rnnzct(irow);
   if(verb), fprintf(' - %i dep-rows\n', m-sprk); end
   probData.data_changed = 1;
end

%----- delete zero columns -----
zrcol = (max(abs(A)) == 0)';
if any(zrcol == 1)
   probData.Zrcols_exist = 1;
   izrcol = find(zrcol);
   if any(c(izrcol) < 0 & ub(izrcol) > opts.big-1)
      if(verb), fprintf('\nPreprocessor: problem unbounded below\n'); end
      probData.isFeasible = 0;
      return;
   end
   xzrcol = zeros(size(izrcol))...
          + (c(izrcol) < 0).*ub(izrcol)...
          + (c(izrcol) > 0).*lb(izrcol);
   inzcol = find(1 - zrcol);
   A = A(:,inzcol);
   c = c(inzcol);
   lb = lb(inzcol);
   ub = ub(inzcol);
   if(verb), fprintf(' - %i 0-columns\n', nnz(zrcol)); end
   probData.data_changed = 1;
   probData.izrcol = izrcol;
   probData.xzrcol = xzrcol;
   probData.inzcol = inzcol;
end

%----- solve singleton rows -----
singleton = (rnnzct == 1);
nsgrows = nnz(singleton);
if nsgrows >= max(1, .01*size(A,1))
   probData.Sgtons_exist = 1;
   isgrows = find(singleton);
   iothers = find(1 - singleton);
   if(verb), fprintf(' - %i singletons\n',nsgrows); end

   Atmp = A(isgrows,:); Atmp1 = spones(Atmp); btmp = b(isgrows);
   if nsgrows == 1
      isolved  = find(Atmp1);
      insolved = find(Atmp1 == 0);
      xsolved  = b(isgrows)/Atmp(isolved);
   else
      colnnzct = sum(Atmp1);
      isolved  = find(colnnzct);
      insolved = find(colnnzct == 0);
      [ii, jj] = find(Atmp);
      Atmp = Atmp(ii,jj); btmp = btmp(ii);
      xsolved  = btmp./diag(Atmp);
      if any(colnnzct >  1)
         repeat = diff([0; jj]) == 0;
         for i = 1:length(xsolved) - 1
             if repeat(i+1) && xsolved(i+1) ~= xsolved(i)
                if(verb), fprintf('\nPreprocessor: problem infeasible\n'); end
                probData.isFeasible = 0;
                return;
             end;
         end;
         ii = find(~repeat); jj = ii;
         Atmp = Atmp(ii,jj); btmp = btmp(ii);
         xsolved  = btmp./diag(Atmp);
      end;
   end;

   if (any(xsolved < lb(isolved)) || any(xsolved > ub(isolved)))
      if(verb), fprintf('\nPreprocessor: problem infeasible\n'); end
      probData.isFeasible = 0;
      return;
   end;

   b = b(iothers) - A(iothers,isolved)*xsolved;
   A = A(iothers, insolved);
   c = c(insolved);
   lb = lb(insolved);
   ub = ub(insolved);
   probData.data_changed = 1;
   probData.isolved = isolved;
   probData.insolved = insolved;
   probData.xsolved = xsolved;
end

%----- shift nonzero lower bounds -----
probData.Lbounds_non0 = any(lb ~= 0);
if (probData.Lbounds_non0)
   b = b - A*lb;
   probData.data_changed = 1;
end

%----- find upper bounds -----
iubounds = ub < opts.big - 1;
probData.Ubounds_exist = full(any(iubounds));
if (probData.Ubounds_exist)
   ub = sparse(iubounds.*(ub-lb));
   probData.nub = nnz(ub);
end

[m, n] = size(A); probData.NNZA = nnz(A);
if(verb), fprintf('[m n] = [%d %d]\n',m,n); end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% plotdata() - plot concatenated multichannel data epochs in two-column format
%
% Usage:   >> plotdata(data)
%          >> plotdata(data,frames)
%          >> plotdata(data,frames,limits,title,channames,colors,rtitle,ydir)
%
% Necessary input:
%   data       = data consisting of consecutive epochs of (chans,frames)
%
% Optional inputs:
%   frames     = time frames/points per epoch {default: 0 -> data length}
%  [limits]    = [xmin xmax ymin ymax]  (x's in ms)
%                {default|0 (or both y's 0) -> use data limits)
%  'title'     = plot title {default|0 -> none}
%  'channames' = channel location file or structure (see readlocs())
%                {default|0 -> channel numbers}
%  'colors'    = file of color codes, 3 chars per line
%                ( '.' = space) {default|0 -> default color order}
%  'rtitle'    = right-side plot title {default|0 -> none}
%   ydir       = y-axis polarity (1 -> pos-up; -1 -> neg-up)
%                {default|0 -> set from default YDIR in 'icadefs.m'}
%
% Authors: Scott Makeig, Arnaud Delorme, Tzyy-Ping Jung,
%          SCCN/INC/UCSD, La Jolla, 05-01-96
%
% See also: plottopo(), timtopo(), envtopo(), headplot(), compmap(), eegmovie()

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

% 5-1-96 from showerps.m  -sm from showerp.m -tpj
% 5-3-96 added default channel numbering, frames & title -sm
% 5-17-96 added nargin tests below -sm
% 5-21-96 added right titles -sm
% 6-29-96 removed Postscript file query -sm
% 7-22-96 restored lines to fill printed page with figure -sm
% 7-29-96 added [channumbers] option for channames argument. -sm
%         changed "channels" argument to "channames" in help statement above -sm
% 1-6-97  debugged min/max time and +/- printing -tpj & sm
% 3-3-97  restored previous Default axis parameters at end -sm
% 4-2-97  debugged 32-epoch plotting -sm
% 5-10-97 added no-args check -sm
% 5-20-97 read icadefs.m for MAXPLOTDATACHANS and MAXPLOTDATAEPOCHS -sm
% 6-23-97 use returns instead of errorcodes -sm
% 10-31-97 use normalized PaperUnits for US/A4 compatibility,
%          fixed [xy]{min,max} printing, added clf, adding Clipping off,
%          fixed scaling, added limits tests -sm & ch
% 11-19-97 removed an 'orient' command that caused problems printing -sm
% 07-15-98 added 'ydir' arg, made pos-up the default -sm
% 07-24-98 fixed 'ydir' arg, and pos-up default -sm
% 01-02-99 added warning about frames not dividing data length -sm
% 02-19-99 debugged axis limits -sm
% 11-21-01 add compatibility to load .loc files for channames -ad
% 01-25-02 reformated help & license, added links -ad
% 02-16-02 added axcopy -ad & sm
% 03-15-02 added readlocs and the use of eloc input structure -ad
% 03-15-02 added smart axcopy -ad

function plotdata(data,frames,limits,plottitle,channels,colors,righttitle,ydr)

if nargin < 1,
    help plotdata
    return
end

%
% Set defaults
%
FONTSIZE = 12;     % font size to use for labels
TICKFONTSIZE=12;   % font size to use for axis labels
icadefs;           % read MAXPLOTDATACHANS constant from icadefs.m
                   % read YDIR
if ~exist('YDIR')
    fprintf('YDIR not read from ''icadefs.m'' - plotting positive-up.\n');
    YDIR = 1;
end
DEFAULT_SIGN = YDIR;  % default to plotting positive-up (1) or negative-up (-1)
                      % Now, read sign from YDIR in icadefs.m
ISSPEC = 0;           % default - 0 = not spectral data, 1 = pos-only spectra

axcolor= get(0,'DefaultAxesXcolor'); % find what the default x-axis color is
plotfile = 'plotdata.ps';
ls_plotfile = 'ls -l plotdata.ps';

%
%%%%%%%%%%%%%%%%%%%%%%%%%% Substitute defaults for missing parameters %%%%%
%
SIGN = DEFAULT_SIGN; % set default ydir from YDIR in icadefs.m
if nargin < 8
   ydr = [];
end
if ~isempty(ydr)& ydr ~=0  % override default from commandline
 if ydr == 1
   SIGN = 1;
 else
   SIGN = -1;
 end
end

if nargin < 7,
    righttitle = 0;
end;
if nargin < 6
    colors = 0;
end
if nargin < 5
    channels = [1:size(data,1)]; % default channames = 1:chans
end
if nargin < 4
    plottitle = 0; %CJH
end
limitset = 0;
if nargin < 3,
    limits = 0;
elseif length(limits)>1
    limitset = 1;
end
if nargin < 2,
    frames = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%% Test parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  [chans,framestotal]=size(data);             % data size
  if frames <=0,
    frames = framestotal;    % default
    datasets=1;
  elseif frames==1,
    fprintf('plotdata: cannot plot less than 2 frames per trace.\n');
    return
    datasets=1;
  else
    datasets = fix(framestotal/frames);        % number of traces to overplot
    if datasets*frames < framestotal
         fprintf('\nWARNING: %d points at end of data will not be plotted.\n',...
                    framestotal-datasets*frames);
    end
  end;

  if chans>MAXPLOTDATACHANS,
    fprintf('plotdata: not set up to plot more than %d channels.\n',...
                       MAXPLOTDATACHANS);
    return
  end;
  if datasets>MAXPLOTDATAEPOCHS
      fprintf('plotdata: not set up to plot more than %d epochs.\n',...
                       MAXPLOTDATAEPOCHS);
    return
  end;
  if datasets<1
      fprintf('plotdata: cannot plot less than 1 epoch!\n');
      return
  end;
%
%%%%%%%%%%%%% Extend the size of the plotting area in the window %%%%%%%%%%%%
%
  curfig = gcf;
  h=figure(curfig);
  set(h,'Color',BACKCOLOR); % set the background color
  set(h,'PaperUnits','normalized'); % use percentages to avoid US/A4 difference
  set(h,'PaperPosition',[0.0235308 0.0272775 0.894169 0.909249]); % equivalent
  % orient portrait
  axis('normal');
%
%%%%%%%%%%%%%%%%%%%% Read the channel names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if isnumeric(channels) & channels(1)==0,
    channels = [1:size(data,1)];
  end;
  if isnumeric(channels),
      channames = num2str(channels(:));                   %%CJH
  else
      [tmp channames] = readlocs(channels);
      channames = strvcat(channames{:});
  end;
%    chid = fopen(channels,'r');
%    if chid <3,
%        fprintf('plotdata(): cannot open file %s.\n',channels);
%        return
%     end;
%     if isempty( findstr( lower(channels), '.loc') )
%    		channames = fscanf(chid,'%s',[4 Inf]);
%    		channames = channames';
%     else
%        	channames = fscanf(chid,'%d %f %f  %s',[7 128]);
%    		channames = char(channames(4:7,:)');
%     end;
%     ii = find(channames == '.');
%     channames(ii) = ' ';

    %channames = channames';
    %   [r c] = size(channames);
    %for i=1:r
    %    for j=1:c
    %        if channames(i,j)=='.',
    %            channames(i,j)=' ';
    %        end;
    %    end;
    %end;
%    end; % setting channames

%
%%%%%%%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if colors ~=0,
    if ~isstr(colors)
       fprintf('plotdata(): color file name must be a string.\n');
       return
    end
    cid = fopen(colors,'r');
    % fprintf('cid = %d\n',cid);
    if cid <3,
        fprintf('plotdata: cannot open file %s.\n',colors);
        return
    end;
    colors = fscanf(cid,'%s',[3 Inf]);
    colors = colors';
       [r c] = size(colors);
    for i=1:r
        for j=1:c
            if colors(i,j)=='.',
                colors(i,j)=' ';
            end;
        end;
    end;
  else % use default color order (no yellow!)
     colors =['r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  '];
     colors = [colors; colors];  % make > 64 available
  end;
  for c=1:size(colors,1)   % make white traces black unless axis color is white
    if colors(c,1)=='w' & axcolor~=[1 1 1]
         colors(c,1)='k';
    end
  end

%
%%%%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if limits==0,      % == 0 or [0 0 0 0]
    xmin=0;
    xmax=frames-1;
    ymin=min(min(data));
    ymax=max(max(data));
    yrange = ymax-ymin;
    ymin = ymin - 0.00*yrange;
    ymax = ymax + 0.00*yrange;
  else
    if length(limits)~=4,
      fprintf( ...
       'plotdata: limits should be 0 or an array [xmin xmax ymin ymax].\n');
      return
    end;
    xmin = limits(1);
    xmax = limits(2);
    ymin = limits(3);
    ymax = limits(4);
  end;

  if xmax == 0 & xmin == 0,
    x = (0:1:frames-1);
    xmin = 0;
    xmax = frames-1;
  else
    dx = (xmax-xmin)/(frames-1);
    x=xmin*ones(1,frames)+dx*(0:frames-1); % compute x-values
  end;
  if xmax<=xmin,
      fprintf('plotdata() - xmax must be > xmin.\n')
      return
  end

  if ymax == 0 & ymin == 0,
      ymax=double(max(max(data)));
      ymin=double(min(min(data)));
      yrange = ymax-ymin;
      % ymin = ymin - 0.00*yrange;
      % ymax = ymax + 0.00*yrange;
  end
  if ymax<=ymin,
      fprintf('plotdata() - ymax must be > ymin.\n')
      return
  end

  xlabel = 'Time (ms)';
  if ymin >= 0 & xmin >= 0,          % For all-positive (spectral) data
    ISSPEC = 1;
    SIGN = 1;
    fprintf('\nPlotting positive up. Assuming data are spectra.\n');
    xlabel = 'Freq (Hz)';
    ymin = 0;                        % plot positive-up
  end;

%
%%%%%%%%%%%%%%%%%%%%%%%% Set up plotting environment %%%%%%%%%%%%%%%%%%%%%%%%%
%
  h = gcf;
  % set(h,'YLim',[ymin ymax]);       % set default plotting parameters
  % set(h,'XLim',[xmin xmax]);
  % set(h,'FontSize',18);
  % set(h,'DefaultLineLineWidth',1); % for thinner postscript lines
%
%%%%%%%%%%%%%%%%%%%%%%%%%% Print plot info %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  clf;   % clear the current figure

  % print plottitle over (left) sbplot 1
  if plottitle==0,
      plottitle = '';
  end
  if righttitle==0,
      righttitle = '';
  end

  if ~isempty(righttitle)
      sbplot(ceil(chans/2),2,2), h=gca;%title([righttitle], 'FontSize',FONTSIZE); % title plot and
      set(h,'FontSize',FONTSIZE);            % choose font size
      set(h,'YLim',[ymin ymax]);            % set default plotting parameters
      set(h,'XLim',[xmin xmax]);
  end;

  msg = ['\nPlotting %d traces of %d frames with colors: '];
  for c=1:datasets
    msg = [msg  colors(mod(c-1,length(colors))+1,:)];
  end
  msg = [msg ' -> \n'];    % print starting info on screen . . .
  fprintf(...
    '\n  limits: [xmin,xmax,ymin,ymax] = [%4.1f %4.1f %4.2f %4.2f]\n',...
                xmin,xmax,ymin,ymax);
  fprintf(msg,datasets,frames);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

  xdiff=xmax-xmin;
  ydiff=ymax-ymin;

  for P=0:datasets-1, %  for each data epoch
      fprintf('\ntrace %d: ',P+1);

    for I=1:chans,        % for each data channel

        index=(2*((rem(I-1,ceil(chans/2))+1)))-1+floor(2*(I-1)/chans);
                                      % = 1 3 5 .. 2 4 6 ..
                                      % plot down left side of page first
        h=sbplot(ceil(chans/2),2,index); h=gca;
        pos = get(h,'position');
        set(h,'position',[pos(1)-pos(4)*0.5 pos(2), pos(3),pos(4)*2]); % increase sbplot-height

        hold on;
        set(h,'YLim',[ymin ymax]);    % set default plotting parameters
        set(h,'XLim',[xmin xmax]);

        axislcolor = get(gca,'Xcolor');   %%CJH

        %
        %%%%%%%%%%%%%%%%%%%%% Plot two-sided time-series data %%%%%%%%%%%%%%%%%%%
        %
        if ~ISSPEC % not spectral data

            ymin = double(ymin);
            ymax = double(ymax);

            if ymin == ymax, ymin = ymin-1; ymax = ymax+1; end;
            plot(x,SIGN*data(I,1+P*frames:1+P*frames+frames-1),colors(mod(P,length(colors))+1));

            if SIGN > 0
                axis([xmin xmax ymin ymax]);          % set axis bounds (pos up)
            else
                axis([xmin xmax -1*ymax -1*ymin]);     % set axis bounds (neg up)
            end

            if P==datasets-1,            % on last traces
                if I==floor((chans+1)/2),   % draw +/0 on lowest left plot
                    signx = double(xmin-0.04*xdiff);

                    if SIGN > 0  % pos up
                        axis('off');hl=text(signx,ymin,num2str(ymin,3));        % text ymin
                        axis('off');hi=text(signx,ymax,['+' num2str(ymax,3)]);  % text +ymax
                    else         % neg up
                        axis('off');hl=text(signx,-1*ymin,num2str(ymin,3));        % text ymin
                        axis('off');hi=text(signx,-1*ymax,['+' num2str(ymax,3)]);  % text +ymax
                    end
                    set(hl,'FontSize',TICKFONTSIZE);         % choose font size
                    set(hl,'HorizontalAlignment','right','Clipping','off');
                    set(hi,'FontSize',TICKFONTSIZE);         % choose font size
                    set(hi,'HorizontalAlignment','right','Clipping','off');
                end

                if I==chans & limitset,    % draw timescale on lowest right plot
                    ytick = double(-ymax-0.25*ydiff);

                    tick = [int2str(xmin)]; h=text(xmin,ytick,tick); % min time
                    set(h,'FontSize',TICKFONTSIZE);         % choose font size
                    set(h,'HorizontalAlignment','center',...
                          'Clipping','off');  % center text

                    tick = [xlabel]; h=text(xmin+xdiff/2,ytick,tick); % xlabel
                    set(h,'FontSize',TICKFONTSIZE);         % choose font size
                    set(h,'HorizontalAlignment','center',...
                          'Clipping','off');  % center text

                    tick = [int2str(xmax)]; h=text(xmax,ytick,tick); % max time
                    set(h,'FontSize',TICKFONTSIZE);         % choose font size
                    set(h,'HorizontalAlignment','center',...
                      'Clipping','off');  % center text
                end;
        end;
      %
      %%%%%%%%%%%%%%%%%%%%% Plot spectral data positive-up [0,ymax] %%%%%%%%%%%%%%%%%%%%%%%%
      %
      else % ISSPEC
        ymin=0;
        plot(x,SIGN*data(I,1+P*frames:1+P*frames+frames-1),colors(mod(P,length(colors))+1));
        ymaxm = ymax;

        % ymin = 0.01;
        % ymaxm = 10.^ceil(log(ymax)/log(10.));
        % if ymaxm/2. > ymax,
        %    ymaxm = ymaxm/2.;
        % end;

        axis([xmin xmax ymin ymaxm]);      % set axis values

        if P==datasets-1,                  % on last trace
         if I==floor((chans+1)/2),         % draw +/0 on lowest left plot
          signx = xmin-0.04*xdiff;

          axis('off');h=text(signx,ymax,['+' num2str(ymax,3)]);
            set(h,'FontSize',TICKFONTSIZE);
            set(h,'HorizontalAlignment','right','Clipping','off');

          axis('off');h=text(signx,0,'0');
            set(h,'FontSize',TICKFONTSIZE);
            set(h,'HorizontalAlignment','right','Clipping','off');
         end;

         if I==chans,                    % draw freq scale on lowest right plot
            ytick = -0.25*ymax;

          tick = [num2str(round(10*xmin)/10) ]; h=text(xmin,ytick,tick);
            set(h,'FontSize',TICKFONTSIZE);
            set(h,'HorizontalAlignment','center','Clipping','off');

          tick = [xlabel]; h=text(xmin+xdiff/2,ytick,tick);
            set(h,'FontSize',TICKFONTSIZE);
            set(h,'HorizontalAlignment','center','Clipping','off');

          tick = [num2str(round(10*xmax)/10) ]; h=text(xmax,ytick,tick);
            set(h,'FontSize',TICKFONTSIZE);
            set(h,'HorizontalAlignment','center','Clipping','off');

        end; % if last chan
      end % if last data
     end; % if ~ISSPEC

      %%%%%%%%%%%%%%%%%%%%%%% Print channel names and lines %%%%%%%%%%%%%%%%%%%%%%%%%%

      if P==datasets-1
          if ~ISSPEC
              axis('off');
              plot([0 0],[ymin ymax],'color',axislcolor); % draw vert %%CJH
                                                          % axis at time 0
          else  % ISSPEC
              axis('off');plot([xmin xmin],[0 ymax],'color',axislcolor);
          end

          % secondx = 200;                                % draw second vert axis
          % axis('off');plot([secondx secondx],[ymin ymax],'color',axislcolor);

          axis('off');
          plot([xmin xmax],[0 0],'color',axislcolor);     % draw horizontal axis

          if ~isempty(channels),                          % print channames
              if ~ISSPEC
                  if ymin <= 0 & ymax >= 0,
                  yht = 0;
                  else
                      yht = nan_mean(SIGN*data(I,1+P*frames:1+P*frames+frames-1));
                  end
                  axis('off'),h=text(xmin-0.04*xdiff,double(yht),[channames(I,:)]);
                  set(h,'HorizontalAlignment','right');      % print before traces
                  set(h,'FontSize',FONTSIZE);                % choose font size

                  % axis('off'),h=text(xmax+0.10*xdiff,yht,[channames(I,:)]);
                  %    set(h,'HorizontalAlignment','left');      % print after traces

              else % ISSPEC
                  axis('off'),h=text(xmin-0.04*xdiff,ymax/2,[channames(I,:)]);
                  set(h,'HorizontalAlignment','right');      % print before traces
                  set(h,'FontSize',FONTSIZE);                % choose font size

                  % axis('off'),h=text(xmax+0.10*xdiff,ymax/2,[channames(I,:)]);
                  %    set(h,'HorizontalAlignment','left');      % print after traces

              end;
          end;
      end;

     fprintf(' %d',I);
    end; % subplot
  end; % dataset
  fprintf('\n');
%
%%%%%%%%%%%%%%%%%% Make printed figure fill page %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  curfig = gcf;
  h=figure(curfig);
  % set(h,'PaperPosition',[0.2 0.3 7.6 10]); % stretch out the plot on the page
%
%%%%%%%%%%%%%%%%%% Restore plot environment %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  % set(h,'DefaultAxesYLim',aylim);      % restore previous plotting parameters
  % set(h,'DefaultAxesXLim',axlim);
  % set(h,'DefaultAxesFontSize',axfont);

%%%%%%%%%%%%%%%%%% Add axcopy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if plottitle
	h = textsc(plottitle, 'title');
	set(h, 'fontsize', FONTSIZE);
end;
axcopy(gcf, 'axis on');

if 0,    % START DETOUR XXXXXXXXXXXXX
%
%%%%%%%%%%%%%%%%%% Save plot to disk if asked %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
   if  plotfile ~= '', %
     n=0; y=1;
     answer = input('plotdata: Save plot as Postscript file? (y/n) ');
     if answer==1,
         fprintf('\nSaving figure as %s ... ',plotfile);
         curfig = gcf;
         h=figure(curfig);
         % set(h,'PaperPosition',[0.2 0.3 7.6 10]);
                                     % stretch out the plot on the page
         eval (['print -dpsc ' plotfile]);
         fprintf('saved. Move or remove file!\n');
         unix(ls_plotfile);
     end
   end
end       % END DETOUR XXXXXXXXXXXXX

function out = nan_mean(in)

    nans = find(isnan(in));
    in(nans) = 0;
    sums = sum(in);
    nonnans = ones(size(in));
    nonnans(nans) = 0;
    nonnans = sum(nonnans);
    nononnans = find(nonnans==0);
    nonnans(nononnans) = 1;
    out = sum(in)./nonnans;
    out(nononnans) = NaN;

author: PatternRecognition

repo: OpenBMI


text:
function training(video, method, frames)

opts.expDir = ['net/' method '/' num2str(frames) '/' video];

opts.train.batchSize = 5 ;
opts.train.numEpochs = 20;
opts.train.continue = true ;
opts.train.useGpu = true ;
opts.train.learningRate = 1e-3;
opts.train.expDir = opts.expDir;

% --------------------------------------------------------------------
% Prepare data
% --------------------------------------------------------------------
imgDir = ['../' video '/input'];
labelDir = ['../' video '/GT'];

grayDir = fullfile('../result/', method, '/', num2str(num_frames), video);

imdb = getImdb_new(imgDir, labelDir, grayDir);

mask = imread(['../' video '/ROI.bmp']);
mask = mask(:, :, 1);
A = max(max(mask));
mask(mask == A) = 1;
if size(mask, 1) > 400 || size(mask,2) > 400
    mask = imresize(mask, 0.5, 'nearest');
end
imdb.mask = single(double(mask));
imdb.half_size = 15;

%%%%%%Yi%%%%%% redefined the net
load('net');

net.layers{1} = struct('type', 'conv', ...
    'filters', 0.01 * randn(7, 7, 4, 32, 'single'), ...
    'biases', zeros(1, 32, 'single'), ...
    'stride', 1, ...
    'pad', 0) ;

net.layers{end-1} = struct('type', 'conv', ...
    'filters', 0.1*randn(1,1,64,1, 'single'), ...
    'biases', zeros(1, 1, 'single'), ...
    'stride', 1, ...
    'pad', 0) ;
net.layers{end} = struct('type', 'sigmoidcrossentropyloss');

load('meanPixel.mat');
imdb.meanPixel = meanPixel;

[net,info] = cnn_train_adagrad(net, imdb, @getBatch,...
    opts.train, 'errorType', 'euclideanloss', ...
    'conserveMemory', true);
end

function [im, labels, mask] = getBatch(imdb, batch)

half_size = imdb.half_size;
meanPixel = imdb.meanPixel;
meanPixel(:,:,4) =0;

for ii = 1 : numel(batch)
    imagename = imdb.images.name{batch(ii)};
    im_ii = single(imread(imagename));

    labelname = imdb.images.labels{batch(ii)};
    roi = imread(labelname);
    labels_ii = zeros(size(roi, 1), size(roi, 2));
    labels_ii( roi == 50 )  = 0.25;       %shade
    labels_ii( roi == 170 ) = 0.75;       %object boundary
    labels_ii( roi == 255 ) = 1;          %foreground

    % resize the image to half size
    if size(im_ii, 1) > 400 || size(im_ii, 2) >400
        im_ii = imresize(im_ii, 0.5, 'nearest');
        labels_ii = imresize(labels_ii, 0.5, 'nearest');
    end

    grayname =imdb.images.gray_name{batch(ii)};

    im_ii(:,:,4) = single(imread(grayname));

    im_large = padarray(im_ii, [half_size, half_size], 'symmetric');
    im_ii = bsxfun(@minus, im_large, meanPixel);

    im(:, :, :, ii) = im_ii;
    labels(:, :, 1, ii) = labels_ii;
    labels(:, :, 2, ii) = double(imdb.mask);
end
end

function imdb = getImdb_new(imgDir, labelDir, grayDir)
files = dir([imgDir '/*.jpg']);
label_files = dir([labelDir '/*.png']);

names = {};
labels = {};
gray_names = {};

for ii = 1:numel(files)
    names{end+1} = [imgDir '/' files(ii).name];
    labels{end+1} = [labelDir '/' label_files(ii).name];

    %prob_name = strrep(label_files(ii).name, 'gt', 'in');

    prob_name = strrep(label_files(ii).name, 'gt', 'in');

    gray_names{end+1} = [grayDir '/' prob_name];

end

imdb.images.set = ones(1,numel(names));
imdb.images.name = names ;
imdb.images.gray_name = gray_names;
imdb.images.labels = labels;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
A=sprand(200,200000,0.05);

tic
AAt=mexCalcAAt(A);
t=toc;
fprintf('mex-file time: %fs\n',t);

tic
AAt2=A*A';
t=toc;
fprintf('matlab time: %fs\n',t);

sum((AAt(:)-AAt2(:)).^2)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab