text:
function [sos,info] = sossolve(sos,solver_options)
% SOSSOLVE --- Solve a sum of squares program.
%
% SOSP = sossolve(SOSP)
%
% SOSP is the SOS program to be solved.
%
% Alternatively, SOSP = sossolve(SOSP,SOLVER_OPT) also defines the solver
% and/or the solver-specific options respectively by fields
%
% SOLVER_OPT.solver (name of the solver)
% SOLVER_OPT.params (a structure containing solver-specific parameters)
%
% The default values for solvers is 'SeDuMi' with parameter ALG = 2, which
% uses the xz-linearization in the corrector and parameter tol =1e-9. See
% SeDuMi help files or user manual for more detail.
%
% Using a second output argument such as [SOSP,INFO] = sossolve(SOSP) will
% return in INFO numerous information concerning feasibility  and CPU time
% that is generated by the SDP solver.
%

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.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, see <http://www.gnu.org/licenses/>.
%

% 12/25/01 - SP
% 01/05/02 - SP - primal
% 01/07/02 - SP - objective
% aug/13 - JA,GV - CDSP,SDPNAL,SDPA solvers and SOS matrix decomposition

if (nargin==1)
   %Default options from old sossolve
   solver_options.solver = 'sdpt3';
   solver_options.params.tol = 1e-9;
   solver_options.params.alg = 2;
elseif ((nargin==2) & ~isnumeric('solver_options') )%2 arguments given,
    if ~isfield(solver_options,'solver')
       solver_options.solver = 'sdpt3';
    end
    if ~isfield(solver_options,'params')
        solver_options.params.tol = 1e-9;%default values for SeDuMi
        solver_options.params.alg = 2;
    end
end

%whenever nargin>=2 options are overwritten
if (nargin==3)
    error('Current SOSTOOLS version does not support call to sossolve with 3 arguments, see manual.');
end;

if ~isempty(sos.solinfo.x)
    error('The SOS program is already solved.');
end;

% Adding slack variables to inequalities
sos.extravar.idx{1} = sos.var.idx{sos.var.num+1};
% SOS variables
I = [find(strcmp(sos.expr.type,'ineq')), find(strcmp(sos.expr.type,'sparse')), find(strcmp(sos.expr.type,'sparsemultipartite'))];
if ~isempty(I)
    sos = addextrasosvar(sos,I);
end;
% SOS variables type II (restricted on interval)
I = find(strcmp(sos.expr.type,'posint'));
if ~isempty(I)
    sos = addextrasosvar2(sos,I);
end;

% Processing all expressions

Atf = []; bf = [];
for i = 1:sos.expr.num
    Atf = [Atf, sos.expr.At{i}];
    bf = [bf; sos.expr.b{i}];
end;

% Processing all variables
[At,b,K,RR] = processvars(sos,Atf,bf);

% Objective function
c = sparse(size(At,1),1);

%% Added by PAP, for compatibility with MATLAB 6.5
if isempty(sos.objective);
    sos.objective = zeros(size(c(1:sos.var.idx{end}-1)));
end
%% End added stuff

c(1:sos.var.idx{end}-1) = c(1:sos.var.idx{end}-1) + sos.objective;       % 01/07/02
c = RR'*c;

pars = solver_options.params;

if strcmp(lower(solver_options.solver),'sedumi')
    % SeDuMi in action
    disp(['Size: ' num2str(size(At))]);
    disp([' ']);
    [x,y,info] = sedumi(At,b,c,K,pars);

elseif strcmp(lower(solver_options.solver),'sdpt3')
    % SDPT3 in action
    smallblkdim = 50;
    save sostoolsdata_forSDPT3 At b c K smallblkdim;
    [blk,At2,C2,b2] = read_sedumi('sostoolsdata_forSDPT3.mat');
    delete sostoolsdata_forSDPT3.mat;
    [obj,X,y,Z,infoSDPT] = sqlp(blk,At2,C2,b2,pars);
    x = zeros(length(c),1);
    cellidx = 1;
    if K.f ~= 0
        x(1:K.f) = X{1}(:);
        cellidx = 2;
    end;
    if K.s(1) ~= 0
        idxX = 1;
        idx = K.f+1;
        smblkdim = 100;
        deblkidx = find(K.s > smblkdim);
        spblkidx = find(K.s <= smblkdim);
        blknnz = [0 cumsum(K.s.*K.s)];
        for i = deblkidx
            dummy = X{cellidx};
            x(idx+blknnz(i):idx+blknnz(i+1)-1) = dummy(:);
            cellidx = cellidx+1;
        end;
        for i = spblkidx
            dummy = X{cellidx}(idxX:idxX+K.s(i)-1,idxX:idxX+K.s(i)-1);
            x(idx+blknnz(i):idx+blknnz(i+1)-1) = dummy(:);
            idxX = idxX+K.s(i);
        end;
    end;
	info.cpusec = infoSDPT.cputime;
	info.iter = infoSDPT.iter;
	if infoSDPT.termcode == 1
		info.pinf = 1;
	else
		info.pinf = (infoSDPT.pinfeas>0.1);
	end;
	if infoSDPT.termcode == 2
		info.dinf = 1;
	else
		info.dinf = (infoSDPT.dinfeas>0.1);
	end;
	if infoSDPT.termcode<= 0
		info.numerr = infoSDPT.termcode;
	else
		info.numerr = 0;
	end;

elseif strcmp(lower(solver_options.solver),'csdp') %6/6/13 JA CSDP interface
    %CSDP in action
    if exist('solver_options.params')
        pars = solver_options.params;
    else
        pars.objtol = 1e-9;
        pars.printlevel = 1;
    end
    if (isfield(K,'f')) %Convert free vars to non-negative LP vars
            n_free = K.f;
            [A,b,c,K] = convertf(At,b,c,K); %K.f set to zero
            At = A';
    end
    c = full(c);
    [x,y,z,info] = csdp(At,b,c,K,pars);
    c=sparse(c);
    % 7/6/13 JA Remove extra entries from x corresponding to LP vars
    if (isfield(K,'f')) %Convert free vars to non-negative LP vars
           index = [n_free+1:2*n_free];
           x(1:n_free) = x(1:n_free)-x(index);
           x(index) = [];
           At(index,:)=[];
           c(index) = [];
    end
    if info~=0
       info.dinf=1;
       info.pinf=1;
    else
       info.dinf=0;
       info.pinf=0;
    end;

elseif strcmp(lower(solver_options.solver),'sdpnal') %6/11/13 JA SDPNAL interface
    % SDPNAL in action
    smallblkdim = 50;
    save sostoolsdata_forSDPNAL At b c K smallblkdim;
    [blk,At2,C2,b2] = read_sedumi('sostoolsdata_forSDPNAL.mat');
    delete sostoolsdata_forSDPNAL.mat;
    try
        [obj,X,y,Z,info] = sdpnal(blk,At2,C2,b2,pars); %run history not returned;
    catch
        [obj,X,y,Z] = sdpnal(blk,At2,C2,b2,pars); %run history not returned;
        info = [];
    end

    x = zeros(length(c),1);
    cellidx = 1;
    if K.f ~= 0
        x(1:K.f) = X{1}(:);
        cellidx = 2;
    end;
    if K.s(1) ~= 0
        idxX = 1;
        idx = K.f+1;
        smblkdim = 100;
        deblkidx = find(K.s > smblkdim);
        spblkidx = find(K.s <= smblkdim);
        blknnz = [0 cumsum(K.s.*K.s)];
        for i = deblkidx
            dummy = X{cellidx};
            x(idx+blknnz(i):idx+blknnz(i+1)-1) = dummy(:);
            cellidx = cellidx+1;
        end;
        for i = spblkidx
            dummy = X{cellidx}(idxX:idxX+K.s(i)-1,idxX:idxX+K.s(i)-1);
            x(idx+blknnz(i):idx+blknnz(i+1)-1) = dummy(:);
            idxX = idxX+K.s(i);
        end;
    end;

    if ~isempty(info)
        if info.pinfeas>1e-6||info.pinfeas>1e-6
            info.dinf=1;
            info.pinf=1;
        else
            info.dinf=0;
            info.pinf=0;
        end;
    else
        info.dinf=1;
        info.pinf=1;
    end

    elseif strcmp(lower(solver_options.solver),'sdpa')
    % SDPA in action

    disp(['Size: ' num2str(size(At))]);
    disp([' ']);
    [x,y,info]=sedumiwrap(At',b,c,K,[],pars);

    info.primalObj
    info.primalError

    info.dualObj
    info.dualError

    if info~=0
       info.dinf=1;
       info.pinf=1;
    else
       info.dinf=0;
       info.pinf=0;
    end;

end;
disp([' ']);
disp(['Residual norm: ' num2str(norm(At'*x-b))]);
disp([' ']);
sos.solinfo.x = x;
sos.solinfo.y = y;
sos.solinfo.RRx = RR*x;
sos.solinfo.RRy = RR*(c-At*y);    % inv(RR') = RR
sos.solinfo.info = info;
sos.solinfo.solverOptions = solver_options;
disp(info)

%return;

% Constructing the (primal and dual) solution vectors and matrices
% If you want to have them, comment/delete the return command above.
% In the future version, these primal and dual solutions will be computed only
% when they are needed. We don't want to store redundant info.

for i = 1:sos.var.num
    switch sos.var.type{i}
    case 'poly'
        sos.solinfo.var.primal{i} = sos.solinfo.RRx(sos.var.idx{i}:sos.var.idx{i+1}-1);
        sos.solinfo.var.dual{i} = sos.solinfo.RRy(sos.var.idx{i}:sos.var.idx{i+1}-1);
    case 'sos'
        primaltemp = sos.solinfo.RRx(sos.var.idx{i}:sos.var.idx{i+1}-1);
        dualtemp = sos.solinfo.RRy(sos.var.idx{i}:sos.var.idx{i+1}-1);
        sos.solinfo.var.primal{i} = reshape(primaltemp,sqrt(length(primaltemp)),sqrt(length(primaltemp)));
        sos.solinfo.var.dual{i} = reshape(dualtemp,sqrt(length(dualtemp)),sqrt(length(dualtemp)));
    end;
end;

for i = 1:sos.extravar.num
    primaltemp = sos.solinfo.RRx(sos.extravar.idx{i}:sos.extravar.idx{i+1}-1);
    dualtemp = sos.solinfo.RRy(sos.extravar.idx{i}:sos.extravar.idx{i+1}-1);
    sos.solinfo.extravar.primal{i} = reshape(primaltemp,sqrt(length(primaltemp)),sqrt(length(primaltemp)));
    sos.solinfo.extravar.dual{i} = reshape(dualtemp,sqrt(length(dualtemp)),sqrt(length(dualtemp)));
end;

sos.solinfo.decvar.primal = sos.solinfo.RRx(1:sos.var.idx{1}-1);
sos.solinfo.decvar.dual = sos.solinfo.RRy(1:sos.var.idx{1}-1);

% ====================================================================================
function sos = addextrasosvar(sos,I)
% Adding slack SOS variables to inequalities

for i = I

    numstates = size(sos.expr.Z{i},2);%GV&JA 6/12/2013

    % Creating extra variable
    maxdeg = full(max(sum(sos.expr.Z{i},2)));
    mindeg = full(min(sum(sos.expr.Z{i},2)));
    Z = monomials(numstates,[floor(mindeg/2):ceil(maxdeg/2)]);
    %disp(['Original : ',num2str(size(Z,1))]);

    % Discarding unnecessary monomials
    maxdegree = max(sos.expr.Z{i},[],1)/2;
    mindegree = min(sos.expr.Z{i},[],1)/2;
    j = 1;
    while (j <= size(Z,1))
        Zdummy1 = maxdegree-Z(j,:);
        Zdummy2 = Z(j,:)-mindegree;
        idx = find([Zdummy1, Zdummy2]<0);
        if ~isempty(idx)
            Z = [Z(1:j-1,:); Z(j+1:end,:)];
        else
            j = j+1;
        end;
    end;
    %disp(['Optimized : ',num2str(size(Z,1))]);

    % Convex hull algorithm
    if strcmp(sos.expr.type{i},'sparse')
        Z2 = sos.expr.Z{i}/2;
        Z = inconvhull(full(Z),full(Z2));
        Z = makesparse(Z);
        %disp(['Optimized again : ',num2str(size(Z,1))]);
    end;

    if strcmp(sos.expr.type{i},'sparsemultipartite')
        Z2 = sos.expr.Z{i}/2;
        info2 = sos.expr.multipart{i};%the vectors of independent variables
        sizeinfo2m = length(info2);
        vecindex = [];
        for indm = 1:sizeinfo2m%for each set of independent variables
            sizeinfo2n(indm) = length(info2{indm});
            for indn = 1:sizeinfo2n(indm)
                if isfield(sos,'symvartable')%
                    varcheckindex = find(info2{indm}(indn)==sos.symvartable);
                    if ~isempty(varcheckindex)
                        vecindex{indm}(indn) = varcheckindex;
                    else
                        vecindex{indm}(indn) = length(info2{1})+find(info2{indm}(indn)==sos.varmat.symvartable);%GV&JA 6/12/2013
                    end

                else
                    % PJS 9/12/13: Update code to handle polynomial objects
                    var = info2{indm}(indn);
                    cvartable = char(sos.varmat.vartable);

                    if ispvar(var)
                        % Convert to string representation
                        var = var.varname;
                    end
                    varcheckindex = find(strcmp(var,sos.vartable));
                    if ~isempty(varcheckindex)
                        vecindex{indm}(indn) = varcheckindex;
                    else
                        vecindex{indm}(indn) = length(info2{1}) + find(strcmp(var,cvartable));
                    end

                    % PJS 9/12/13: Original Code to handle polynomial objects
                    %vecindex{indm}(indn) = find(strcmp(info2{indm}(indn).varname,sos.vartable));
                end;
            end
        end
        Z = sparsemultipart(full(Z),full(Z2),vecindex);
        Z = makesparse(Z);
    end;

    dimp = size(sos.expr.b{i},2);

    % Adding slack variables
    sos.extravar.num = sos.extravar.num + 1;
    var = sos.extravar.num;
    sos.extravar.Z{var} = makesparse(Z);
    [T,ZZ] = getconstraint(Z);
    sos.extravar.ZZ{var} = ZZ;
    sos.extravar.T{var} = T';
    %sos.extravar.idx{var+1} = sos.extravar.idx{var}+size(Z,1)^2;%GVcomment the next slack variable starts in the column i+dim(Z)^2 - the elements of the vectorized square matrix

    sos.extravar.idx{var+1} = sos.extravar.idx{var}+(size(Z,1)*dimp)^2;
    for j = 1:sos.expr.num
        sos.expr.At{j} = [sos.expr.At{j}; ...
                sparse(size(sos.extravar.T{var},1)*dimp^2,size(sos.expr.At{j},2))];
    end

    ZZ = flipud(ZZ);
    T = flipud(T);

    Zcheck = sos.expr.Z{i};
    %this is for the matrix case

    if dimp==1
        % JFS 6/3/2003: Ensure correct size:
        pc.Z = sos.extravar.ZZ{var};
        pc.F = -speye(size(pc.Z,1));
        [R1,R2,newZ] = findcommonZ(sos.expr.Z{i},pc.Z);
        % JFS 6/3/2003: Ensure correct size:

        if isempty(sos.expr.At{i})
            sos.expr.At{i} = sparse(size(sos.expr.At{i},1),size(R1,1));
        end
        %------------
        sos.expr.At{i} = sos.expr.At{i}*R1;
        lidx = sos.extravar.idx{var};
        uidx = sos.extravar.idx{var+1}-1;
        sos.expr.At{i}(lidx:uidx,:) = sos.expr.At{i}(lidx:uidx,:) - sos.extravar.T{var}*pc.F*R2;
        sos.expr.b{i} = R1'*sos.expr.b{i};
        sos.expr.Z{i} = newZ;

    else

        [R1,R2,Znew] = findcommonZ(Zcheck,ZZ);

        R1 = fliplr(R1);
        R2 = fliplr(R2);
        Znew = flipud(Znew);

        R1sum = sum(R1,1);
        T = R2'*T;

        ii = 1;
        sig_ZZ = size(ZZ,1);
        sig_Z = size(Z,1);
        sig_Znew = size(Znew,1);

        Tf = sparse(dimp^2*sig_Znew,(dimp*sig_Z)^2);
        Sv = sparse(sig_Znew*dimp^2,1);
        for j = 1:sig_Znew
            Mt0 = sparse(dimp,dimp*sig_Z^2);
            for k = 1:sig_Z
                Mt0(:, (dimp*sig_Z)*(k-1)+1:(dimp*sig_Z)*k) = kron(eye(dimp),T(j,(sig_Z)*(k-1)+1:(sig_Z)*k));
            end
            Tf((j-1)*dimp^2+1:j*dimp^2,:) = kron(eye(dimp),Mt0);

            if R1sum(j)==1
                Sv((j-1)*dimp^2+1:j*dimp^2)= reshape(sos.expr.b{i}( dimp*(ii-1)+1:dimp*ii,:)',dimp^2,1);
                if ii<size(Zcheck,1)
                    ii = ii+1;
                end
            else
                sos.expr.At{i} = [ sos.expr.At{i}(:,1:(j-1)*dimp^2) sparse(size(sos.expr.At{i},1),dimp^2) sos.expr.At{i}(:,(j-1)*dimp^2+1:end)];
            end
        end

        lidx = sos.extravar.idx{var};
        uidx = sos.extravar.idx{var+1}-1;
        sos.expr.At{i}(lidx:uidx,:) =  Tf';
        sos.expr.b{i} =  Sv;

    end
end;

% ====================================================================================
function sos = addextrasosvar2(sos,I)
% Adding slack SOS variables type II

numstates = size(sos.expr.Z{1},2);
for i = I
    % Creating extra variable
    maxdeg = full(max(sum(sos.expr.Z{i},2)));
    mindeg = full(min(sum(sos.expr.Z{i},2)));
    Z = monomials(numstates,[floor(mindeg/2):ceil(maxdeg/2)]);

    % Discarding unnecessary monomials
    maxdegree = max(sos.expr.Z{i},[],1)/2;
    mindegree = min(sos.expr.Z{i},[],1)/2;
    j = 1;
    while (j <= size(Z,1))
        Zdummy1 = maxdegree-Z(j,:);
        Zdummy2 = Z(j,:)-mindegree;
        idx = find([Zdummy1, Zdummy2]<0);
        if ~isempty(idx)
            Z = [Z(1:j-1,:); Z(j+1:end,:)];
        else
            j = j+1;
        end;
    end;

    % Add the variables

    for k = 0:1

    sos.extravar.num = sos.extravar.num + 1;
    var = sos.extravar.num;
    sos.extravar.Z{sos.extravar.num} = makesparse(Z);
    [T,ZZ] = getconstraint(Z);
    sos.extravar.ZZ{var} = ZZ;
    sos.extravar.T{var} = T';
    sos.extravar.idx{var+1} = sos.extravar.idx{var}+size(Z,1)^2;
    for j = 1:sos.expr.num
        sos.expr.At{j} = [sos.expr.At{j}; ...
                sparse(size(sos.extravar.T{var},1),size(sos.expr.At{j},2))];
    end;

    % Modifying expression
    degoffset = [k sparse(1,numstates-1)];
    pc.Z = sos.extravar.ZZ{var} + degoffset(ones(size(sos.extravar.ZZ{var},1),1),:);
    pc.F = -speye(size(pc.Z,1));
    [R1,R2,newZ] = findcommonZ(sos.expr.Z{i},pc.Z);
    % JFS 6/3/2003: Ensure correct size:
    if isempty(sos.expr.At{i})
       sos.expr.At{i} = sparse(size(sos.expr.At{i},1),size(R1,1));
    end
    %------------
    sos.expr.At{i} = sos.expr.At{i}*R1;
    lidx = sos.extravar.idx{var};
    uidx = sos.extravar.idx{var+1}-1;
    sos.expr.At{i}(lidx:uidx,:) = sos.expr.At{i}(lidx:uidx,:) - sos.extravar.T{var}*pc.F*R2;
    sos.expr.b{i} = R1'*sos.expr.b{i};
    sos.expr.Z{i} = newZ;

    Z = Z(find(Z<maxdegree));    % Discard the unnecessary monomial for the second variable

    end;

end;

% ====================================================================================
function [At,b,K,RR] = processvars(sos,Atf,bf)
% Processing all variables

% Decision variables
K.s = [];
K.f = sos.var.idx{1}-1;
RR = speye(K.f);

% Polynomial and SOS variables
for i = 1:sos.var.num
    switch sos.var.type{i}
    case 'poly'
        sizeX = sos.var.idx{i+1}-sos.var.idx{i};
        startidx = sos.var.idx{i};
        RR = spantiblkdiag(RR,speye(sizeX));
        K.f = sizeX+K.f;
    case 'sos'
        sizeX = sqrt(sos.var.idx{i+1}-sos.var.idx{i});
        startidx = sos.var.idx{i};
        RR = spblkdiag(RR,speye(sizeX^2));
        K.s = [K.s sizeX];
    end;
end;

for i = 1:sos.extravar.num
    sizeX = sqrt(sos.extravar.idx{i+1}-sos.extravar.idx{i});
    startidx = sos.extravar.idx{i};
    RR = spblkdiag(RR,speye(sizeX^2));
    K.s = [K.s sizeX];
end;

At = RR'*Atf;

b = bf;

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [success,msg,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  pxpermm] = ...
  ReadArenaParameters_Wrapper(InputDataType,varargin)

success = false;
msg = '';

[inmoviefile,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  pxpermm,leftovers] = myparse_nocheck(varargin,...
  'inmoviefile','',...
  'arenatype','None','arenacenterx',0,'arenacentery',0,...
  'arenaradius',123,'arenawidth',123,'arenaheight',123,...
  'pxpermm',1);

switch InputDataType,
  case {'Ctrax','CtraxPlusWings'},

    [success,msg,...
      arenatype,arenacenterx,arenacentery,...
      arenaradius,arenawidth,arenaheight,...
      pxpermm] = ...
      ReadArenaParameters_Ctrax(...
      'inmoviefile',inmoviefile,...
      leftovers{:},...
      'pxpermm',pxpermm,...
      'arenatype',arenatype,...
      'arenacenterx',arenacenterx,...
      'arenacentery',arenacentery,...
      'arenaradius',arenaradius,...
      'arenawidth',arenawidth,...
      'arenaheight',arenaheight);

  case 'LarvaeRiveraAlba',

    success = false;
    msg = 'Not implemented';

  case 'Motr',

    success = false;
    msg = 'No arena parameters stored in Motr files';

  case 'Qtrax',

    [success,msg,pxpermm] = ...
      ReadArenaParameters_Qtrax(...
      leftovers{:},...
      'pxpermm',pxpermm);

  case 'MAGATAnalyzer',

    [success,msg,pxpermm] = ...
      ReadArenaParameters_MAGATAnalyzer(...
      leftovers{:},...
      'pxpermm',pxpermm);

  case 'MWT',

    [success,msg,pxpermm] = ...
      ReadArenaParameters_MWT(...
      leftovers{:},...
      'pxpermm',pxpermm);

  case 'LarvaeReid',

    [success,msg,pxpermm] = ...
      ReadArenaParameters_LarvaeReid(...
      leftovers{:},...
      'pxpermm',pxpermm);

  case 'JCtrax',

    [success,msg,pxpermm] = ...
      ReadArenaParameters_JCtrax(...
      leftovers{:},...
      'pxpermm',pxpermm);

  otherwise
    success = false;
    msg = sprintf('Unknown data type %s',InputDataType);
end

author: kristinbranson

repo: JAABA


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
classdef SortGui < handle
    properties(Constant)
        PROP_ORDER='sortGui.sortOrder';
        PROP_ASCENDING='sortGui.sortAscending';
        PROP_SEARCH='sortGui.search';
        PROP_RESTRICT='sortGui.restrict';
    end

    properties
        fncRefresh;
        dlg;
    end

    properties(SetAccess=private)
        givenTip=false;
        path;
        strings;
        allChb;
        allMsg;
        sortIdxs;
        visibleIdxs;
        sortTypes;
        sortKeys;
        btnUp;
        btnDown;
        btnSortDir;
        btnSearch;
        options;
        sortAscending=true;
        nOptions;
        checkBoxes
        innerPnl;
        cb;
        map;
        propPrefix;
        allChbPnl;
        searchPnl;
        searchJtf;
        searching=false;
        found;
        foundLabel;
        foundIdx;
        btnRestrict;
        originalSize;
        originalHeight;
        firstCb;
        sortLabel;
    end

    methods
        function this=SortGui(dlg, allChb, allMsg, allChbPnl, options, ...
                checkBoxes, innerPnl)
            pnl=Gui.Panel;
            this.dlg=dlg;
            this.allChb=allChb;
            this.allMsg=allMsg;
            this.options=options;
            this.nOptions=length(options);
            this.sortIdxs=1:this.nOptions;
            this.checkBoxes=checkBoxes;
            this.innerPnl=innerPnl;
            [this.sortKeys, this.sortTypes]=Html.DecodeSort(options{1});
            this.sortLabel=Gui.Label('Sort:');
            if ~isempty(this.sortKeys)
                pnl.add(this.sortLabel);
                fl=javaObjectEDT('java.awt.FlowLayout', 0,0,0);
                pnl.setLayout(fl);
                labels=['sort by' this.sortKeys];
                cb=Gui.Combo(labels, 1,[],[],@(h,e)sortItems(this, h), ...
                    'Change order of list');
                pnl.add(cb);
                this.cb=cb;
                allChbPnl.add(Gui.Label('  '), 'Center');
            end
            this.path=BasicMap.Path;
            this.btnSortDir=Gui.ImageButton(...
                fullfile(this.path, 'sortDown.png'), ...
                'List order is ascending, click for DESCENDING?', ...
                @(h,e)toggleAscending(this));
            pnl.add(this.btnSortDir);
            this.btnSearch=Gui.ImageButton(...
                fullfile(this.path, 'find16.gif'), ...
                'Search items in this list', ...
                @(h,e)toggleSearch(this));
            pnl.add(this.btnSearch);
            this.initSearch;
            allChbPnl.add(pnl, 'East');
            this.allChbPnl=allChbPnl;
            this.seeSortCues;
        end

        function setSortDir(this)
            if this.sortAscending
                this.btnSortDir.setIcon(Gui.Icon('sortDown.png'));
                this.btnSortDir.setToolTipText(...
                    'Order is ascending, click for DESCENDING.');
            else
                this.btnSortDir.setIcon(Gui.Icon('sortUp.png'));
                this.btnSortDir.setToolTipText(...
                    'Order is descending, click for ASCENDING.');
            end
        end

        function toggleAscending(this)
            this.sortIdxs=flip(this.sortIdxs);
            this.refreshSearchAndSort(true);
            this.seeSortCues;
        end

        function seeSortCues(this)
            sorting=~isempty(this.cb) && this.cb.getSelectedIndex>0;
            this.btnSortDir.setVisible(sorting);
            this.sortLabel.setVisible(sorting);
        end

        function refreshSearchAndSort(this, sorting)
            if ~isempty(this.fncRefresh)
                if this.btnRestrict.isSelected
                    pnl=feval(this.fncRefresh);
                else
                    was=this.visibleIdxs;
                    this.visibleIdxs=[];
                    pnl=feval(this.fncRefresh);
                    this.visibleIdxs=was;
                end
                if ~isempty(this.dlg)
                    if isempty(this.originalSize)
                        sz=this.dlg.getSize;
                        if sz.height>0 && sz.width>0
                            this.originalSize=sz;
                            this.originalHeight=this.originalSize.height*.8;
                        end
                    end
                    this.dlg.pack;
                    if ~isempty(this.originalSize)
                        sz=this.dlg.getSize;
                        if ~sz.equals(this.originalSize)
                            this.dlg.setSize(this.originalSize);
                        end
                    end
                end
            else
                this.innerPnl.removeAll;
                if isempty(this.visibleIdxs) || ~this.btnRestrict.isSelected
                    for jj=1:this.nOptions
                        idx=this.sortIdxs(jj);
                        this.innerPnl.add(this.checkBoxes.get(idx-1));
                    end
                else
                    for jj=1:this.nOptions
                        idx=this.sortIdxs(jj);
                        if isempty(find(this.visibleIdxs==idx,1))
                            continue;
                        end
                        this.innerPnl.add(this.checkBoxes.get(idx-1));
                    end
                end
                if nargin>1 && sorting && ~isempty(this.dlg)
                    this.dlg.pack;
                else
                    this.innerPnl.repaint;
                end
            end
            this.sortAscending=~this.sortAscending;
            this.setAscending;
            this.setSortDir;
            this.handleFound;
            this.seeFirstCb;
            drawnow;
        end

        function setAscending(this)
            if ~isempty(this.map)
                if this.sortAscending
                    this.map.set(this.propAscending, 'true');
                else
                    this.map.set(this.propAscending, 'false');
                end
            end
        end

        function prop=propOrder(this)
            prop=[this.propPrefix SortGui.PROP_ORDER];
        end

        function prop=propAscending(this)
            prop=[this.propPrefix SortGui.PROP_ASCENDING];
        end

        function prop=propSearch(this)
            prop=[this.propPrefix SortGui.PROP_SEARCH];
        end

        function prop=propRestrict(this)
            prop=[this.propPrefix SortGui.PROP_RESTRICT];
        end

        function setProperties(this, map, propPrefix, searchDflt)
            if nargin<4
                searchDflt=false;
            end
            this.map=map;
            this.propPrefix=propPrefix;
            if map.has(this.propOrder)
                idx=map.getNumeric(this.propOrder, 0);
                if isnumeric(idx) && idx>0
                    this.cb.setSelectedIndex(idx-1);
                    %this.sort(idx);
                end
            end
            if ~map.is(this.propAscending, true)
                this.toggleAscending;
                this.setSortDir;
            end
            if map.is(this.propSearch, searchDflt)
                MatBasics.RunLater(@(h,e)toggleSearch(this), .12);
            end
            if map.is(this.propRestrict, true)
                this.btnRestrict.setSelected(true);
                this.restrict;
            end
        end

        function sortItems(this, h)
            idx=h.getSelectedIndex+1;
            if ~isempty(this.map)
                this.map.set(this.propOrder, idx);
            end
            this.sort(idx);
        end

        function sort(this, idx)
            if idx<=1
                this.sortIdxs=1:this.nOptions;
            else
                isNum=isequal('N', this.sortTypes{idx-1});
                if isNum
                    rm=realmin;
                    values=zeros(1, this.nOptions);
                else
                    values=cell(1, this.nOptions);
                end
                key=this.sortKeys{idx-1};
                for jj=1:this.nOptions
                    v=Html.DecodeSortValues(this.options{jj}, key);
                    if isNum
                        if isempty(v)
                            values(jj)=rm;
                        else
                            values(jj)=str2double(v{1});
                        end
                    else
                        if isempty(v)
                            values{jj}='';
                        else
                            values{jj}=v{1};
                        end
                    end
                end
                [~, this.sortIdxs]=sort(values);
            end
            if this.sortAscending
                this.sortIdxs=flip(this.sortIdxs);
            end
            this.sortAscending=~this.sortAscending;
            this.toggleAscending;
            this.setAscending;
        end

        function handleFound(this)
            if ~isempty(this.foundLabel)
                if isempty(this.found)
                    this.btnDown.setVisible(false);
                    this.btnUp.setVisible(false);
                    if this.searchJtf.getText.length==0
                        this.foundLabel.setText('');
                    else
                        this.foundLabel.setText(...
                            ['0/' ...
                            num2str(this.nOptions) '']);
                    end
                else
                    N=length(this.found);
                    if ~this.btnRestrict.isSelected
                        s=[num2str(this.foundIdx) '/' num2str(N)];
                        this.btnDown.setVisible(true);
                        this.btnUp.setVisible(true);
                    else
                        s=[num2str(N) '/' num2str(this.nOptions)];
                        this.btnDown.setVisible(false);
                        this.btnUp.setVisible(false);
                    end
                    this.foundLabel.setText(s);
                    if ispc
                        this.dlg.pack;
                    end
                end
            end
        end

        function initSearch(this)
            jtf=javaObjectEDT('javax.swing.JTextField');
            jtf.setColumns(15)
            %jtf.setHorizontalAlignment(jtf.RIGHT);
            jj=handle(jtf, 'CallbackProperties');
            set(jj, 'KeyTypedCallback', @(h,e)search(this, h));
            this.searchPnl=Gui.BorderPanel;
            pnl=Gui.Panel;
            fl=javaObjectEDT('java.awt.FlowLayout', 0,0,0);
            pnl.setLayout(fl);
            emptyBorder=javax.swing.BorderFactory.createEmptyBorder(0,0,0,0);
            pnl.setBorder(emptyBorder);
            pnl.add(javaObjectEDT('javax.swing.JLabel', ' '));
            this.btnRestrict=Gui.CheckBox(Html.Wrap(Html.Img(...
                'hide.png')), false, [],...
                [], @(h,e)restrict(this), ...
                Html.Wrap(['Select ' Html.Img('magnify.png') ...
                ' to restrict list to found items']));
            pnl.add(this.btnRestrict);

            this.btnUp=Gui.ImageButton(...
                fullfile(this.path, 'upArrow.png'), ...
                'Find next', ...
                @(h,e)nextSearch(this, -1));
            pnl.add(this.btnUp);
            this.btnDown=Gui.ImageButton(...
                fullfile(this.path, 'downArrow.png'), ...
                'Find previous', ...
                @(h,e)nextSearch(this, 1));
            pnl.add(this.btnDown);
            this.foundLabel=javaObjectEDT('javax.swing.JLabel');
            pnl.add(this.foundLabel);
            pnl.add(jtf);
            btnClear=Gui.ImageButton(...
                fullfile(BasicMap.Path, 'cancel.gif'), ...
                'Clear/close search', @(h,e)clearSearch(this));
            pnl.add(btnClear);
            pnl.add(javaObjectEDT('javax.swing.JLabel', '     '));
            this.searchPnl.add(pnl, 'East');
            this.searchJtf=jtf;
            this.btnDown.setPreferredSize(java.awt.Dimension(12, 16))
            this.btnUp.setPreferredSize(java.awt.Dimension(12, 16))
        end

        function restrict(this)
            if this.btnRestrict.isSelected
                this.btnUp.setVisible(false);
                this.btnDown.setVisible(false);
                if ~isempty(this.map)
                    this.map.set(this.propRestrict, 'true');
                end
            else
                this.btnUp.setVisible(true);
                this.btnDown.setVisible(true);
                if ~isempty(this.map)
                    this.map.set(this.propRestrict, 'false');
                end
            end
            this.sortAscending=~this.sortAscending;
            this.refreshSearchAndSort;
        end

        function toggleSearch(this)
            if this.searching
                this.allChbPnl.remove(this.searchPnl);
                if ~isempty(this.map)
                    this.map.set(this.propSearch, 'false');
                end
            else
                sz=this.innerPnl.getSize;
                if sz.width==0
                    sz=this.checkBoxes.get(0).getParent.getSize;
                end
                if sz.width>0 && sz.width<500
                    where='North';
                    if sz.width<200
                        this.searchJtf.setColumns(9);
                    end
                else
                    where='Center';
                end

                this.allChbPnl.add(this.searchPnl, where);
                this.searchJtf.requestFocus;
                if ~isempty(this.map)
                    this.map.set(this.propSearch, 'true');
                end

            end
            this.allChbPnl.repaint;
            this.searching=~this.searching;
            this.dlg.pack;
        end

        function nextSearch(this, i)
            N=length(this.found);
            this.foundIdx=this.foundIdx+i;
            if this.foundIdx>N
                this.foundIdx=1;
            elseif this.foundIdx<1
                this.foundIdx=N;
            end
            cb_=this.found{this.foundIdx};
            b=cb_.getBounds;
            b.height=b.height*2.5;
            cb_.getParent.scrollRectToVisible(b);
            this.handleFound;
        end

        function clearSearch(this)
            if this.searchJtf.getText.length==0
                this.toggleSearch;
            else
                this.searchJtf.setText('');
                this.found={};
                this.search(this.searchJtf);
                this.handleFound;
            end
        end

        function search(this, h)
            this.found={};
            this.visibleIdxs=[];
            N=this.nOptions-1;
            s=lower(char(h.getText));
            N=this.nOptions-1;
            if isempty(this.strings)
                this.strings={};
                for i=0:N
                    this.strings{end+1}=string(lower(...
                        char(edu.stanford.facs.swing.Basics.RemoveXml(...
                        this.checkBoxes.get(i).getText))));
                end
            end
            this.firstCb=this.checkBoxes.get(0);
            if isempty(s)
                for i=0:N
                    cb_=this.checkBoxes.get(i);
                    cb_.setEnabled(true);
                end
            else
                first=false;
                N=this.nOptions;
                for i=1:N
                    idx=this.sortIdxs(i);
                    cb_=this.checkBoxes.get(idx-1);
                    if this.strings{idx}.contains(s)
                        cb_.setEnabled(true);
                        if ~first
                            this.firstCb=cb_;
                            first=true;
                        end
                        this.found{end+1}=cb_;
                        this.visibleIdxs(end+1)=idx;
                    else
                        cb_.setEnabled(false);
                    end
                end
            end
            this.foundIdx=1;
            if this.btnRestrict.isSelected
                this.sortAscending=~this.sortAscending;
                this.refreshSearchAndSort;
            else
                this.handleFound;
                this.seeFirstCb;
                if ~this.givenTip && h.getText.length==1
                    this.givenTip=true;
                    MatBasics.RunLater(@(h,e)tip, .51);
                end
            end
            %this.dlg.pack;
            function tip
                try
                    BasicMap.Global.showToolTip(this.btnRestrict,...
                        [], 2, -50);
                catch ex
                    disp(ex);
                end
            end
        end

        function seeFirstCb(this)
            cb_=[];
            if ~isempty(this.firstCb)
                cb_=this.firstCb;
            else
                for jj=0:this.nOptions-1
                    if this.checkBoxes.get(jj).isSelected
                        cb_=this.checkBoxes.get(jj);
                    end
                end
            end
            if ~isempty(cb_)
                b=cb_.getBounds;
                b.height=b.height*2.5;
                cb_.getParent.scrollRectToVisible(b);
            end
        end

        function [idxs_, N]=setAllChbText(this)
            [idxs_, N]=Gui.SetAllChb(this.allChb, this.allMsg, this.checkBoxes);
        end

        function idxs=reorder(this, idxs)
            idxs=MatBasics.Sort1D(idxs, this.sortIdxs);
        end
    end

    methods(Static)
        function items=EncodeLabelsAndList(labels, list, gtp)
            if iscell(list)
                N=length(list);
            else
                N=list.size;
            end
            app=BasicMap.Global;
            items=cell(1, N);
            for i=1:N
                if iscell(list)
                    item=list{i};
                else
                    item=list.get(i-1);
                end
                frequency=sum(labels==i);
                cue1=Html.EncodeSort('frequency', frequency);
                cue2=Html.EncodeSort('name', item);
                if nargin>2
                    frequency=gtp.encodeCount(frequency) ;
                else
                    frequency=String.encodeCount(frequency);
                end

                items{i}=['' item app.smallStart...
                        ''...
                        ' <' frequency ...
                        ' events>'  app.smallEnd ...
                        cue1 cue2 ''];
            end
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%  logistic 1st-order logistic function
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function visualizeHOG(w)
% Visualize HOG features/weights.
%   visualizeHOG(w)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
% Copyright (C) 2007 Pedro Felzenszwalb, Deva Ramanan
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% Make pictures of positive and negative weights
bs = 20;
w = w(:,:,1:9);
scale = max(max(w(:)),max(-w(:)));
pos = HOGpicture(w, bs) * 255/scale;
neg = HOGpicture(-w, bs) * 255/scale;

% Put pictures together and draw
buff = 10;
pos = padarray(pos, [buff buff], 128, 'both');
if min(w(:)) < 0
  neg = padarray(neg, [buff buff], 128, 'both');
  im = uint8([pos; neg]);
else
  im = uint8(pos);
end
imagesc(im);
colormap gray;
axis equal;
axis off;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
% cart2topo() - convert xyz-cartesian channel coordinates
%               to polar topoplot() coordinates. Input data
%               are points on a sphere centered at (0,0,0)
%               or at optional input 'center'. This function
%               is now DEPRECATED! See Important warning below.
%
% Usage:  >> [th r] = cart2topo(xyz);   % 3-column [x y z] position data
%         >> [th r] = cart2topo(x,y,z); % separate x,y,z vectors
%         >> [th r x y z] = cart2topo(xyz,'key', 'val', ...);
%         >> [th r x y z] = cart2topo(x,y,z,'key', 'val', ...);
%
% Optional inputs:
%    'center'  = [X Y Z] use a known sphere center {Default: [0 0 0]}
%    'squeeze' = plotting squeeze factor (0[default]->1). -1 -> optimize
%                the squeeze factor automatically so that maximum radius is 0.5.
%    'optim'   = [0|1] find the best-fitting sphere center minimizing the
%                standard deviation of the radius values.
%    'gui'     = ['on'|'off'] pops up a gui for optional arguments.
%                Default is off.
%
% Example: >> [th r] = cart2topo(xyz,[1 0 4]);
%
% Notes: topoplot() does not plot channels with radius>0.5
%        Shrink radii to within this range to plot all channels.
%        [x y z] are returned after the optimization of the center
%        and optionally squeezing r towards it by factor 'squeeze'
%
% Important:
%   DEPRECATED: cart2topo() should NOT be used if elevation angle is less than 0
%   (for electrodes below zero plane) since then it returns INNACURATE results.
%   SUBSTITUTE: Use cart2topo = cart2sph() -> sph2topo().
%
% Authors: Scott Makeig, Luca Finelli & Arnaud Delorme SCCN/INC/UCSD,
%          La Jolla, 11/1999-03/2002
%
% See also: topo2sph(), sph2topo(), chancenter()

% Copyright (C) 11/1999 Scott Makeig, 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

% 3-16-00 improved help message -sm
% 1-25-02 put spherror subfunction inside cart2topo -ad
% 1-25-02 help pops-up if no arguments -ad
% 01-25-02 reformated help & license -ad
% 02-13-02 now center fitting works, need spherror outside cart2topo -lf
% 02-14-02 radii are squeezed of squeeze in to fit topoplot circle -lf
% 03-31-02 center fitting is optional
% 04-01-02 automatic squeeze calculation -ad & sm

function [th,r,xx,yy,zz] = cart2topo(x,varargin)

if nargin<1
    help cart2topo
    return;
end;
if nargin >= 2
	if ~ischar(varargin{1})
		y = varargin{1};
		z = varargin{2};
		varargin = varargin(3:end);
	end;
end;
if exist('y') ~= 1
	if size(x,2)==3 % separate 3-column data
		z = x(:,3);
		y = x(:,2);
		x = x(:,1);
	else
		error('Insufficient data in first argument');
	end
end;

g = [];
if ~isempty(varargin)
    try, g = struct(varargin{:});
    catch, error('Argument error in the {''param'', value} sequence'); end;
end;

try, g.optim;      catch, g.optim = 0; end;
try, g.squeeze;    catch, g.squeeze = 0; end;
try, g.center;     catch, g.center = [0 0 0]; end;
try, g.gui;        catch, g.gui = 'off'; end;

if g.squeeze>1
  fprintf('Warning: Squeeze must be less than 1.\n');
  return
end

if ~isempty(g.center) & size(g.center,2) ~= 3
  fprintf('Warning: Center must be [x y z].\n');
  return
end

% convert to columns
x = -x(:); % minus os for consistency between measures
y = -y(:);
z = z(:);

if any(z < 0)
    disp('WARNING: some electrodes lie below the z=0 plane, result may be innacurate')
    disp('         Instead use cart2sph() then sph2topo().')
end;
if strcmp(g.gui, 'on')
	[x y z newcenter] = chancenter(x, y, z, [], 1);
else
	if g.optim == 1
		[x y z newcenter] = chancenter(x, y, z, []);
	else
		[x y z newcenter] = chancenter(x, y, z, g.center);
	end;
end;
radius = (sqrt(x.^2+y.^2+z.^2));   % assume xyz values are on a sphere
xx=-x; yy=-y; zz=z;
x = x./radius; % make radius 1
y = y./radius;
z = z./radius;

r = x; th=x;

for n=1:size(x,1)
  if x(n)==0 & y(n)==0
    r(n) = 0;
  else
    r(n)  = pi/2-atan(z(n)./sqrt(x(n).^2+y(n).^2));
  end
end

r = r/pi;        % scale r to max 0.500
if g.squeeze < 0
    g.squeeze = 1 - 0.5/max(r); %(2*max(r)-1)/(2*rmax);
    fprintf('Electrodes will be squeezed together by %2.3g to show all\n', g.squeeze);
end;
r = r-g.squeeze*r; % squeeze electrodes in squeeze*100% to have all inside

for n=1:size(x,1)
  if abs(y(n))<1e-6
    if x(n)>0
      th(n) = -90;
    else % x(n) <= 0
      th(n) = 90;
    end
  else
    th(n) = atan(x(n)./y(n))*180/pi+90;
    if y(n)<0
      th(n) = th(n)+180;
    end
  end
  if th(n)>180
     th(n) = th(n)-360;
  end
  if th(n)<-180
     th(n) = th(n)+360;
  end
end

author: PatternRecognition

repo: OpenBMI


text:
%                    runBenchmarkTransmissionMatrix.m
%
% Benchmark algorthms using real/empirical data.  The measurement matrix is
% a transmission matrix, and the measurements were acquired using an
% optical aparatus. The data acquisition is described in "Coherent
% Inverse Scattering via Transmission Matrices: Efficient Phase Retrieval
% Algorithms and a Public Dataset."
%    This script will reconstruct images using various algorithms, and
% report the image quality of the results.
%
% See the header of transmissionMatrixExperiment.m, or the paper cited
% above for more details.
%
% This script does the following:
% 1. Set up parameters and create a list of
% algorithm structs.
% 2. Invoke the general benchmark function benchmarkTransmissionMatrix.
%
% The benchmark program compares the performance of specified algorithms on
% reconstructing an image.  By default, this script reconstructs a 16x16
% image.  However, it can also reconstruct a 40x40 or 64x64 image by
% changing the options below.  Note, the runtime and memory requirements
% are very high when reconstructing larger images.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

%clc
%clear
%close all

%% 1.Set up parameters
imageSize = 64;         % Side-length of image to reconstruct. Valid choices: 16,40,64
datasetSelection = 5;   % Which dataset to select measurements from.  Valid choices: 1-5
residualConstant = 0.3;

% Create a list of algorithms structs
wf = struct('initMethod','spectral','algorithm','wirtflow');
twf = struct('algorithm','twf');
rwf = struct('algorithm','rwf');
ampflow = struct('algorithm','amplitudeflow');
taf = struct('initMethod','orthogonal','algorithm','taf');
raf = struct('initMethod','weighted','algorithm','raf');
fienup = struct('algorithm','fienup');
gs = struct('algorithm','gerchbergsaxton');
cd = struct('algorithm','coordinatedescent','maxIters',50000);
kac = struct('algorithm','kaczmarz','maxIters',1000);
pmax = struct( 'algorithm','phasemax','maxIters',1000);
plamp = struct('algorithm', 'phaselamp','maxIters',1000);
scgm = struct('algorithm','sketchycgm');
plift = struct('algorithm','phaselift','initMethod','orthogonal');

% Grab your pick of algorithms.
algorithms = {plift};

%% 2. Run benchmark
benchmarkTransmissionMatrix(imageSize, datasetSelection, residualConstant, algorithms)

author: tomgoldstein

repo: phasepack-matlab


text:
% EEG_INSERTBOUND - insert boundary event in an EEG event structure.
%
% Usage:
%       >> [eventout newind] = eeg_insertbound( eventin, pnts, ...
%                                                abslatency, duration);
% Required Inputs:
%   eventin    - EEGLAB event structure (EEG.event)
%   pnts       - data points in EEG dataset (EEG.pnts * EEG.trials)
%   abslatency - absolute latency of regions in original dataset. Can
%                also be an array of [beg end] latencies with one row
%                per region removed. Then 'lengths' argument is ignored.
% Optional Inputs:
%   lengths    - lengths of removed regions
%
% Outputs:
%   eventout   - EEGLAB event output structure with added boundaries
%   newind     - Indices of the new events
%
% Notes:
%   This function performs the following:
%   1) add boundary events to the 'event' structure;
%        remove nested boundaries;
%        recompute the latencies of all events.
%   2) all latencies are given in (float) data points.
%        e.g., a latency of 2000.3 means 0.3 samples (at EEG.srate)
%              after the 2001st data frame (since first frame has latency 0).
%
% Author: Arnaud Delorme and Hilit Serby, SCCN, INC, UCSD, April, 19, 2004
%
% See also: EEG_EEGREJ, POP_MERGESET

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

function [eventin, newind] = eeg_insertbound( eventin, pnts, regions, lengths)

    if nargin < 3
        help eeg_insertbound;
        return;
    end

    eeglab_options;
    regions = round(regions);
    regions(regions < 1) = 1;
    regions(regions > pnts) = pnts;
    for i=2:size(regions,1)
        if regions(i-1,2) >= regions(i,1)
            regions(i,1) = regions(i-1,2)+1;
        end
    end

    if ~isempty(regions)
        fprintf('eeg_insertbound(): %d boundary (break) events added.\n', size(regions, 1));
    else
        return;
    end

    % recompute latencies of boundevents (in new dataset)
    % ---------------------------------------------------
    [~, tmpsort] = sort(regions(:,1));
    regions        = regions(tmpsort,:);
    lengths = regions(:,2)-regions(:,1)+1;

    if ~isempty(eventin) && isfield(eventin, 'latency')
         eventLatencies = [ eventin.latency ];
         if length(eventLatencies) ~= length(eventin)
             error('Some events do not have latencies')
         end
    else
        eventLatencies = [];
    end
    newEventLatencies = eventLatencies;
    oriLen            = length(eventin);
    rmEvent           = [];
	for iRegion = 1:size(regions,1) % sorted in decreasing order

        % find event succeeding boundary to insert event
        % at the correct location in the event structure
        % ----------------------------------------------
        tmpind    = find( eventLatencies - regions(iRegion,1) > 0 );
        if ~isempty(newEventLatencies)
            newEventLatencies(tmpind) = newEventLatencies(tmpind)-lengths(iRegion);
        end

        % insert event
        % ------------
        [tmpnest, addlength ]  = findnested(eventin, eventLatencies, regions(iRegion,:));
        rmEvent = [ rmEvent tmpnest ];
        %if regions(iRegion,1) % do not remove first event

            eventin(end+1).type   = eeg_boundarytype(eventin);
            eventin(end).latency  = regions(iRegion,1)-sum(lengths(1:iRegion-1))-0.5;
            eventin(end).duration = lengths(iRegion,1)+addlength;
        %end
    end

    % copy latencies
    % --------------
    if ~isempty(newEventLatencies)
        for iEvent = 1:oriLen
            eventin(iEvent).latency = newEventLatencies(iEvent);
        end
    end
    eventin(rmEvent) = [];

    % resort events
    % -------------
    if ~isempty(eventin) && isfield(eventin, 'latency')
%       eventin([ eventin.latency ] < 1) = [];
        eventin([ eventin.latency ] < 0) = [];
        alllatencies = [ eventin.latency ];
        [~, sortind] = sort(alllatencies);
        eventin = eventin(sortind);
        newind = sortind(oriLen+1:end);
    end

    if ~isempty(rmEvent)
        fprintf('eeg_insertbound(): boundary events inserted and %d events removed.\n', length(rmEvent));
    end

% look for nested events
% return indices of nested events and
% their total length
% -----------------------------------
function [ indEvents, addlen ] = findnested(event, eventlat, region)
    indEvents = find( eventlat > region(1) & eventlat < region(2));

    if ~isempty(event) && isfield(event,'type') && isfield(event, 'duration')
        boundaryInd = eeg_findboundaries( event(indEvents) );
        addlen      = sum( [ event(indEvents(boundaryInd)).duration ] );
    else
        addlen = 0;
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function reset_mem()
    global mem config;
    mem.layer_inputs = {};
    mem.activations = {};
    mem.orig_activation_size = {};
    mem.mask_inputs = {};
    mem.mask_activations = {};
    mem.mask_intermediate = {};
    mem.deltas = {};
    mem.conv2conv = {};
    mem.pool2conv = {};
    mem.gen_out_matrix = {};
    mem.convBpool = {};
    mem.convBconv = {};
    mem.pooling_matrix ={};
    mem.output = 0;
    mem.one_over_add_counts = 0;
    mem.grads = {};
    if strcmp(config.forward_pass_scheme{1}, 'conv_v_sr')
        % for SR only
        mem.grads{1} = {};
    end
end

author: jimmy-ren


text:
function o = spoptimget(options,name,default,flag)
% SPOPTIMGET   Get sparse grid optimization OPTIONS parameters.
%   VAL = SPOPTIMGET(OPTIONS,'NAME') extracts the value of the
%   named property from the sparse grid options structure OPTIONS,
%   returning an empty matrix if the property value is not
%   specified in OPTIONS. It is sufficient to type only the leading
%   characters that uniquely identify the property. Case is ignored
%   for property names. [] is a valid OPTIONS argument.
%
%   VAL = SPOPTIMGET(OPTIONS,'NAME',DEFAULT) extracts the named
%   property as above, but returns VAL = DEFAULT if the named
%   property is not specified in OPTIONS. For example
%
%       val = spoptimget(opts, 'Maximize', 'off');
%
%   returns val = 'off' if the Maximize property is not specified in opts.
%
%   See also SPOPTIMSET, SPCOMPSEARCH, SPFMINSEARCH, SPCGSEARCH,
%            SPMULTISTART

% Author : Andreas Klimke
% Version: 1.2
% Date   : November 13, 2006

% Change log:
% V1.2   : November 13, 2006
%          Cleaned up options.
% V1.1   : August 8, 2005
%          Added BisectGap and funTol parameters
% V1.0   : June 9, 2005
%          Initial version

% Note: SPOPTIMGET is similar in syntax and code to the options
% handling with ODEGET and ODESET of the MATLAB ODE suite by Marc
% Reichelt and Lawrence Shampine.

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

% undocumented usage for fast access with no error checking
if (nargin == 4) & isequal(flag,'fast')
   o = getknownfield(options,name,default);
   return
end

if nargin < 2
  error('Not enough input arguments.');
end
if nargin < 3
  default = [];
end

if ~isempty(options) & ~isa(options,'struct')
  error(['First argument must be an options structure created with' ...
	 ' SPOPTIMSET.']);
end

if isempty(options)
  o = default;
  return;
end

Names = {'Minimize', 'Maximize', 'TolFun', 'TolX', ...
	 'MaxIter', 'StartPoint', 'TestCorners', 'PrevResult', ...
	 'Method', 'NumStarts', 'OptimsetOptions', 'Display'};

m = length(Names);
for k = 1:m
  names{k} = lower(Names{k});
end

lowName = lower(name);
matched = strmatch(lowName,names);
if isempty(matched)               % if no matches
  error(sprintf(['Unrecognized property name ''%s''.  ' ...
                 'See SPOPTIMSET for possibilities.'], name));
elseif length(matched) > 1
  msg = sprintf('Ambiguous property name ''%s'' ', name);
  msg = [msg '(' Names{matched}];
  for k = j(2:length(matched))'
    msg = [msg ', ' Names{matched}];
  end
  msg = sprintf('%s).', msg);
  error(msg);
end

if any(strcmp(fieldnames(options),Names{matched}))
  o = options.(Names{matched});
  if isempty(o)
    o = default;
  end
else
  o = default;
end

%-----------------------------------------------------------------------
function v = getknownfield(s, f, d)
% GETKNOWNFIELD   Get field f from struct s, or else yield default d.

if isfield(s,f)   % s could be empty.
  v = subsref(s, struct('type','.','subs',f));
  if isempty(v)
    v = d;
  end
else
  v = d;
end

author: johannesgerer

repo: jburkardt-m


text:
function value = r4_sign ( x )

%*****************************************************************************80
%
%% R4_SIGN returns the sign of an R4.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( 0 <= x )
    value = +1.0;
  else
    value = -1.0;
  end

  return
end


text:
classdef OMAS_modelViewer
    properties
        figureProperties;
        targetObject;
        triadPosition = [-2;-2;-2];  % Position of the model
        triadScale = 0.5;
        R        = eye(3);      % The rotation of the body and triad together
        R_offset = eye(3);      % The offset between the body axis triad and geometry
        %triad_b;
        %triad_g;
    end
    methods
        % Constructor
        function obj = OMAS_modelViewer(varargin)
            % Input sanity check
            if numel(varargin) ~= 1
                error('Please provide a target object class to view');
            end
            % Assume input is an object label
            if ischar(varargin{1})
                obj.targetObject = eval(varargin{1});                      % Instantiate the object
            elseif ismember('objectDefinition',superclasses(varargin{1})) % Check for the object root class
                obj.targetObject = varargin{1};
            else
                error('The variable provided is not a valid OpenMAS object.');
            end
            % Get the default figure properties
            [obj.figureProperties] = obj.getFigureProperties();
        end
        % Set the model position
        function [obj] = setTriadPosition(obj,p)
            % Input sanity check
            assert(numel(p) == 3 && size(p,1) == 3,'Position must be a column vector [3x1]');
            % Set the triad's relative position
            obj.triadPosition = p;
        end
        % Set the model rotation
        function [obj] = setRotation(obj,eulers)
            % Input sanity check
            assert(numel(eulers) == 3,'Euler rotations must be a column vector [3x1]');
            % Get the rotation matrix
            obj.R = OMAS_geometry.eulersToRotationMatrix(eulers);
        end
        % Show scene
        function [figureHandle] = show(obj,objectClass)

            % Input checking
            if nargin > 1
                target = objectClass;
            else
                target = obj.targetObject;
            end
            % For clarity
            properties = obj.figureProperties;

            % Define plot
            figureHandle = figure('Name',target.name);
            ax = axes(figureHandle);
            hold on;

            % Other plot attributes
            title(target.name,'fontweight',properties.fontWeight,'fontsize',properties.titleFontSize);
            xlabel('x(m)','fontweight',properties.fontWeight,'fontSize',properties.axisFontSize);
            ylabel('y(m)','fontweight',properties.fontWeight,'fontSize',properties.axisFontSize);
            zlabel('z(m)','fontweight',properties.fontWeight,'fontSize',properties.axisFontSize);
            set(ax,'FontSize',properties.axisFontSize,'fontWeight',properties.fontWeight);
            set(ax,'Color',properties.axesColor);
            axis('tight','equal');
            grid on;  box on;
            % AXES-ON LOGICAL
            set(ax,'visible',properties.axesVisible)
            % SET ORIENTATION
            view(properties.orientation);
            % REPRESENT GEOMETRY AS A PATCH
            target.GEOMETRY = OMAS_graphics.normalise(target.GEOMETRY);
            patch(ax,...
                'Vertices',(target.GEOMETRY.vertices*obj.R_offset)*obj.R',...
                'Faces',target.GEOMETRY.faces,...
                'FaceColor',properties.faceColour,...
                'EdgeColor',properties.edgeColour,...
                'EdgeAlpha',properties.edgeAlpha,...
                'FaceAlpha',properties.faceAlpha,...
                'FaceLighting',properties.faceLighting,...
                'LineWidth',0.1);

            % ADD ALIGNED TRIAD
            triadHandle = OMAS_graphics.drawTriad(figureHandle,zeros(3,1),obj.R*obj.R_offset,obj.triadScale);
            set(triadHandle(1,:),'lineWidth',1.2);
            set(triadHandle(1,:),'lineStyle','-');
            % set(triadHandle(1,:),'Color','k');
            % ADD GLOBAL TRIAD
            triadHandle = OMAS_graphics.drawTriad(figureHandle,obj.triadPosition,eye(3),obj.triadScale);
            set(triadHandle(1,:),'lineWidth',1.2);
            set(triadHandle(1,:),'lineStyle','-');
            % set(triadHandle(1,:),'Color','k');
        end
    end
    methods (Static)
        % Get figure properties from 'OMAS_figureProperties'
        function [figureProperties] = getFigureProperties()
            % Get the figure properties structure common to OMAS
            figureProperties = OMAS_figureProperties();
            % Specific default parameters
            figureProperties.backgroundColor = 'w';
            figureProperties.orientation = [0 90]; %[75 44];
            figureProperties.axesVisible = 'off';
            figureProperties.axesColor = 'none';
            figureProperties.edgeAlpha = 0.02;
            figureProperties.edgeColour = 'k';
            figureProperties.edgeWidth = 0.1;
            figureProperties.faceAlpha = 0.05;
            figureProperties.faceColour = 'b';
        end
    end
end

repo: OpenMAS


text:
function [] = test_l1_logistic_regression()

    clc;
    clear;
    close all;

    %% Set algorithms
    if 0
        algorithms = gd_solver_list('ALL');
    else
        %algorithms = {'PG-BKT', 'PG-TFOCS-BKT', 'APG-BKT', 'APG-TFOCS-BKT', 'Newton-CHOLESKY', 'NCG-BKT','L-BFGS-TFOCS'};
        algorithms = {'APG-BKT', 'APG-TFOCS-BKT'};
    end

    %% prepare dataset
    if 1
        % generate synthtic data
        d = 100;
        n = 1000;
        data = logistic_regression_data_generator(n, d);
        x_train = data.x_train;
        y_train = data.y_train;
        x_test = data.x_test;
        y_test = data.y_test;
        d = size(x_train,1);
        w_opt = data.w_opt;
        lambda = 0.1;
    else
        % load pre-created synthetic data
        data = importdata('../data/logistic_regression/data_100d_10000.mat');
        x_train = data.x_train;
        y_train = data.y_train;
        x_test = data.x_test;
        y_test = data.y_test;
        d = size(x_train,1);
        n = length(y_train);
        w_opt = data.w_star;
        lambda = data.lambda;
    end

    %% define problem definitions
    problem = l1_logistic_regression(x_train, y_train, x_test, y_test, lambda);

    %% calculate solution
    if norm(w_opt)
    else
        % calculate solution
        w_opt = problem.calc_solution(1000, 0.05);
    end
    f_opt = problem.cost(w_opt);
    fprintf('f_opt: %.24e\n', f_opt);

    %% initialize
    w_init = rand(d,1);
    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_iter = 300;
        options.verbose = true;

        switch algorithms{alg_idx}
            case {'PG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'PG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'APG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'APG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-TFOCS'}

                options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'BFGS-TFOCS'}

                options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                %options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'NCG-BKT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                %options.step_alg = 'tfocs_backtracking';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    fprintf('\n\n');

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','cost', algorithms, w_list, info_list);
    % display iter vs. l1 norm, i.e. the toral number of non-zero elements
    display_graph('iter','reg', algorithms, w_list, info_list);

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
% Test file for @deltafun/feval.m

function pass = test_feval(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

seedRNG(1337)

%%
f = bndfun(@sin);
d = deltafun(f, struct('deltaMag', 1, 'deltaLoc', 0));
pass(1) = isinf(feval(d, 0)) && feval(d, 0) > 0;
pass(2) = isinf(feval(-d, 0)) && feval(-d, 0) <0;

%%
f = fun.constructor(@(x) sin(x));
d = deltafun(f, struct('deltaMag', [], 'deltaLoc', []));
x = rand(1, 4);
pass(2) = norm(feval(f, x) - feval(d, x), inf) == 0;

%%
x = rand(1,4);
d = deltafun(f, struct('deltaMag', rand(1, 4), 'deltaLoc', x));
pass(3) = all(isinf(feval(d, x)));

%%
x = chebfun('x');
d = dirac(x-1);
pass(4) = isinf(feval(d, 1));
val = feval(-d, 1);
pass(5) = isinf(val) && val < 0;
pass(6) = isinf(feval(d, 'right'));
pass(7) = feval(d, 'left') == 0;
pass(8) = feval(d, 1, 'left' ) == 0;
pass(9) = feval(d, 0, 'left' ) == 0;
pass(10) = feval(d, 0, 'right' ) == 0;
d = diff(heaviside(x));
pass(11) = isinf(feval(d, 0));
pass(12) = feval(d, 'right') == 0;
pass(13) = feval(d, 'left') == 0;
pass(14) = feval(d, 0, 'left' ) == 0;
pass(15) = feval(d, 0, 'right' ) == 0;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Y,XYZ] = spm_read_vols(V,mask)
% Read in entire image volumes
% FORMAT [Y,XYZ] = spm_read_vols(V,mask)
% V    - vector of mapped image volumes to read in (from spm_vol)
% mask - implicit zero mask?
% Y    - 4D matrix of image data, fourth dimension indexes images
% XYZ  - 3xn matrix of XYZ locations returned
%_______________________________________________________________________
%
% For image data types without a representation of NaN (see spm_type),
% implicit zero masking can be used. If mask is set, then zeros are
% treated as masked, and returned as NaN.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_read_vols.m 1143 2008-02-07 19:33:33Z spm $

%-Argument checks
%-----------------------------------------------------------------------
if nargin<2, mask = 0; end
if nargin<1, error('insufficient arguments'), end

spm_check_orientations(V);

%-Read in image data
%-----------------------------------------------------------------------
n = numel(V);         %-#images
Y = zeros([V(1).dim(1:3),n]);       %-image data matrix

for i=1:n,
  for p=1:V(1).dim(3)
    Y(:,:,p,i) = spm_slice_vol(V(i),spm_matrix([0 0 p]),V(i).dim(1:2),0);
  end
end

%-Apply implicit zero mask for image datatypes without a NaNrep
%-----------------------------------------------------------------------
if mask
    %-Work out images without NaNrep
    im = logical(zeros(n,1));
    for i=1:n, im(i)=~spm_type(V(i).dt(1),'NaNrep'); end
    %-Mask
    Y(Y(:,:,:,im)==0)=NaN;
end

%-Return as 3D matrix if single image
%-----------------------------------------------------------------------
if n==1; Y=Y(:,:,:,1); end

%-Compute XYZ co-ordinates (if required)
%-----------------------------------------------------------------------
if nargout>1
    [R,C,P]=ndgrid(1:V(1).dim(1),1:V(1).dim(2),1:V(1).dim(3));
    RCP = [R(:)';C(:)';P(:)'];
    clear R C P
    RCP(4,:)=1;
    XYZ = V(1).mat(1:3,:)*RCP;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function value = p32_exact ( dim_num )

%*****************************************************************************80
%
%% P32_EXACT returns the exact integral for problem 32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  c = [];
  c = p32_r8vec ( 'G', 'C', dim_num, c );

  z = [];
  z = p32_r8vec ( 'G', 'Z', dim_num, z );

  [ a, b ] = p32_lim ( dim_num );

  value = 1.0;

  for i = 1 : dim_num

    if ( z(i) <= a(i) )

      value = value * 0.0;

    elseif ( z(i) <= b(i) )

      if ( c(i) == 0.0 )
        value = value * ( z(i) - a(i) );
      else
        value = value * ( exp ( c(i) * z(i) ) - exp ( c(i) * a(i) ) ) / c(i);
      end

    else

      if ( c(i) == 0.0 )
        value = value * ( b(i) - a(i) );
      else
        value = value * ( exp ( c(i) * z(i) ) - exp ( c(i) * a(i) ) ) / c(i);
      end

    end

  end

  return
end