text:
function EditDataType_CopyFields(topH);
%
% For the Edit Data Type dialog, copy fields to user-selected scans.
%
% EditDataType_CopyFields(topH);
%
%
% ras, 02/06.
% ras, 07/06: streamlined and simplified.

% Unpack the user data from the edit figure
uiData = get(topH, 'UserData');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% have the user select the target scans %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
scanList = {uiData.dataType.scanParams.annotation};
for i=1:length(scanList)
    scanList{i} = sprintf('Scan %i: %s', i, scanList{i});
end
title = 'Copy Current Params to which scans?';
[tgtScans, ok] = listdlg('PromptString', title,...
    'ListSize', [300 400], 'ListString',scanList, ...
    'InitialValue', 1 ,'OKString', 'OK');
if ~ok, return; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% copy over the fields, one by one %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scan Params Fields
nFields = length(uiData.scanData); % Number of scan edit fields
for iField=1:nFields
    editFlag = uiData.scanData(iField).edit;
    if editFlag
        hContent = uiData.scanData(iField).handle;
        fName = uiData.scanData(iField).field;
        value = get(hContent, 'String');

        % if the value is normally numeric, run num2str
        oldVal = uiData.dataType.scanParams(tgtScans(1)).(fName);
        if isnumeric(oldVal)
            value = str2num(value);
        end

        for tgt = tgtScans
            uiData.dataType.scanParams(tgt).(fName) = value;
        end
    end
end

% Blocked (Traveling-Wave) Params Fields
nFields = length(uiData.blockData); % Number of scan edit fields
for iField=1:nFields
    editFlag = uiData.blockData(iField).edit;
    if editFlag
        hContent = uiData.blockData(iField).handle;
        fName = uiData.blockData(iField).field;
        value = get(hContent, 'String');

        % if the value is normally numeric, run num2str
        oldVal = uiData.dataType.blockedAnalysisParams(tgtScans(1)).(fName);
        if isnumeric(oldVal)
            value = str2num(value);
        end

        for tgt = tgtScans
            uiData.dataType.blockedAnalysisParams(tgt).(fName) = value;
        end
    end
end

% Scan Params Fields
nFields = length(uiData.eventData); % Number of scan edit fields
for iField=1:nFields
    editFlag = uiData.eventData(iField).edit;
    if editFlag
        hContent = uiData.eventData(iField).handle;
        fName = uiData.eventData(iField).field;
        value = get(hContent, 'String');

        % if the value is normally numeric, run num2str
        oldVal = uiData.dataType.eventAnalysisParams(tgtScans(1)).(fName);
        if isnumeric(oldVal)
            value = str2num(value);
        end

        for tgt = tgtScans
            uiData.dataType.eventAnalysisParams(tgt).(fName) = value;
        end
    end
end

set(topH, 'UserData', uiData);

return


text:
function result = p00_simpson ( prob, int_num )

%*****************************************************************************80
%
%% P00_SIMPSON applies the composite Simpson rule to integrate a function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem index.
%
%    Input, integer INT_NUM, the number of subintervals.
%
%    Output, real RESULT, the approximate integral.
%
  [ a, b ] = p00_lim ( prob );

  h = ( b - a ) / ( int_num );
%
%  This section of the code was rewritten to vectorize.
%
  x(1:2*int_num+1) = ( ( 2*int_num - ( 0 : 2*int_num ) ) * a   ...
                   + (               ( 0 : 2*int_num ) ) * b ) ...
                     / ( 2*int_num                         );

  if ( 1 == int_num )
    result = ...
                    p00_fun ( prob, 1,           x(1) ) ...
      + 4.0 *       p00_fun ( prob, int_num,     x(2) ) ...
      +             p00_fun ( prob, 1,           x(3) );
  else
    result = ...
                    p00_fun ( prob, 1,           x(1)               ) ...
      + 2.0 * sum ( p00_fun ( prob, int_num - 1, x(3:2:2*int_num-1) ) ) ...
      + 4.0 * sum ( p00_fun ( prob, int_num,     x(2:2:2*int_num)   ) )...
      +             p00_fun ( prob, 1,           x(2*int_num+1)     );
   end

  result = result * h / 6.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [res,resL2,qualMeasOut]=SART_TV(proj,geo,angles,niter,varargin)
% SART_TV solves Cone Beam CT image reconstruction using Oriented Subsets
%              Simultaneous Algebraic Reconstruction Technique algorithm
%
%   SART_TV(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem
%   using the projection data PROJ taken over ALPHA angles, corresponding
%   to the geometry described in GEO, using NITER iterations.
%
%   SART_TV(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter. Default is 1
%
%   'lambda_red':  Reduction of lambda. Every iteration
%                  lambda=lambdared*lambda. Default is 0.99
%
%   'Init':        Describes different initialization techniques.
%                  'none'     : Initializes the image to zeros (default)
%                  'FDK'      : Initializes image to FDK reconstruction
%                  'multigrid': Initializes image by solving the problem in
%                               small scale and increasing it when relative
%                               convergence is reached.
%                  'image'    : Initialization using a user specified
%                               image. Not recommended unless you really
%                               know what you are doing.
%   'InitImg'      an image for the 'image' initialization. Avoid.
%
%   'TViter'       number of iterations in the TV step. Default 50
%
%   'TVlambda'     hyperparameter in TV iteration. It gives the ratio of
%                  importance of the image vs the minimum total variation.
%                  default is 15. Lower means more TV denoising.
%
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
%
%   'QualMeas'     Asks the algorithm for a set of quality measurement
%                  parameters. Input should contain a cell array of desired
%                  quality measurement names. Example: {'CC','RMSE','MSSIM'}
%                  These will be computed in each iteration.
% 'OrderStrategy'  Chooses the subset ordering strategy. Options are
%                  'ordered' : uses them in the input order, but divided
%                  'random'  : orders them randomly
%                  'angularDistance': chooses the next subset with the
%                                     biggest angular distance with the ones used.
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------

%% Deal with input parameters
[lambda,res,lamdbared,verbose,QualMeasOpts,TViter,TVlambda,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);
measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<3 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);

resL2=zeros(1,niter);
if nargout>1
    computeL2=true;
else
    computeL2=false;
end

blocksize=1;
[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy);

angles_reorder=cell2mat(alphablocks);
index_angles=cell2mat(orig_index);
% does detector rotation exist?
if ~isfield(geo,'rotDetector')
    geo.rotDetector=[0;0;0];
end
%% Create weighting matrices

% Projection weight, W
W=computeW(geo,angles,gpuids);

% Back-Projection weight, V
V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids);

if redundancy_weights
    % Data redundancy weighting, W_r implemented using Wang weighting
    % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc

    num_frames = size(proj,3);
    W_r = redundancy_weighting(geo);
    W_r = repmat(W_r,[1,1,num_frames]);
    % disp('Size of redundancy weighting matrix');
    % disp(size(W_r));
    W = W.*W_r; % include redundancy weighting in W
end

%% Iterate
offOrigin=geo.offOrigin;
offDetector=geo.offDetector;
rotDetector=geo.rotDetector;
DSD=geo.DSD;
DSO=geo.DSO;
% TODO : Add options for Stopping criteria
for ii=1:niter
    if (ii==1 && verbose==1);tic;end
    % If quality is going to be measured, then we need to save previous image
    % THIS TAKES MEMORY!
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = res; % only store if necesary
    end

    for jj=1:size(angles,2)
        if size(offOrigin,2)==size(angles,2)
            geo.offOrigin=offOrigin(:,index_angles(:,jj));
        end
        if size(offDetector,2)==size(angles,2)
            geo.offDetector=offDetector(:,index_angles(:,jj));
        end
        if size(rotDetector,2)==size(angles,2)
            geo.rotDetector=rotDetector(:,index_angles(:,jj));
        end
        if size(DSD,2)==size(angles,2)
            geo.DSD=DSD(jj);
        end
        if size(DSO,2)==size(angles,2)
            geo.DSO=DSO(jj);
        end
        %         proj_err=proj(:,:,jj)-Ax(res,geo,angles(:,jj));     %                                 (b-Ax)
        %         weighted_err=W(:,:,jj).*proj_err;                   %                          W^-1 * (b-Ax)
        %         backprj=Atb(weighted_err,geo,angles(:,jj));         %                     At * W^-1 * (b-Ax)
        %         weigth_backprj=bsxfun(@times,1./V(:,:,jj),backprj); %                 V * At * W^-1 * (b-Ax)
        %         res=res+lambda*weigth_backprj;                      % x= x + lambda * V * At * W^-1 * (b-Ax)
        res=res+lambda* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,jj).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids));
        if nonneg
            res=max(res,0);
        end
    end

    % If quality is being measured
    if measurequality
        qualMeasOut(:,ii)=Measure_Quality(res,res_prev,QualMeasOpts);
    end

    lambda=lambda*lamdbared;
    % TV denoising
    res=im3DDenoise(res,'TV',TViter,TVlambda,'gpuids',gpuids);

    if computeL2
        geo.offOrigin=offOrigin;
        geo.offDetector=offDetector;
        geo.DSD=DSD;
        geo.rotDetector=rotDetector;
        resL2(ii)=im3Dnorm(proj(:,:,index_angles)-Ax(res,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax
        % If the error is not minimized.
        if  ii~=1 && resL2(ii)>resL2(ii-1)
            if verbose
                disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]);
            end
            return
        end
    end

    if (ii==1 && verbose==1)
        expected_time=toc*niter;
        disp('SART_TV');
        disp(['Expected duration   :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end
end

end

function initres=init_multigrid(proj,geo,alpha,TViter,TVlambda,gpuids)

finalsize=geo.nVoxel;
% start with 64
geo.nVoxel=[64;64;64];
geo.dVoxel=geo.sVoxel./geo.nVoxel;
if any(finalsize<geo.nVoxel)
    initres=zeros(finalsize');
    return;
end
niter=100;
initres=zeros(geo.nVoxel','single');
while ~isequal(geo.nVoxel,finalsize)

    % solve subsampled grid
    initres=SART_TV(proj,geo,alpha,niter,'Init','image','InitImg',initres,'Verbose',0,'TViter',TViter,'TVlambda',TVlambda,'gpuids',gpuids);

    % Get new dims.
    geo.nVoxel=geo.nVoxel*2;
    geo.nVoxel(geo.nVoxel>finalsize)=finalsize(geo.nVoxel>finalsize);
    geo.dVoxel=geo.sVoxel./geo.nVoxel;
    % Upsample!
    % (hopefully computer has enough memory............)
    [y, x, z]=ndgrid(linspace(1,size(initres,1),geo.nVoxel(1)),...
        linspace(1,size(initres,2),geo.nVoxel(2)),...
        linspace(1,size(initres,3),geo.nVoxel(3)));
    initres=interp3(initres,x,y,z);
    clear x y z
end
end

function [lambda,res,lamdbared,verbose,QualMeasOpts,TViter,TVlambda,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin)
opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','tviter','tvlambda','orderstrategy','nonneg','gpuids','redundancy_weighting','groundtruth'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:SART_TV:InvalidInput','Invalid number of inputs')
end
multigrid=false;
% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:SART_TV:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option is not default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:SART_TV:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        % % % % % % % Verbose
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
            % % % % % % % hyperparameter, LAMBDA
        case 'lambda'
            if default
                lambda=1;
            else
                if length(val)>1 || ~isnumeric(val)
                    error('TIGRE:SART_TV:InvalidInput','Invalid lambda')
                end
                lambda=val;
            end
        case 'lambda_red'
            if default
                lamdbared=0.99;
            else
                if length(val)>1 || ~isnumeric(val)
                    error('TIGRE:SART_TV:InvalidInput','Invalid lambda')
                end
                lamdbared=val;
            end
        case 'init'
            res=[];
            if default || strcmp(val,'none')
                res=zeros(geo.nVoxel','single');
                continue
            end
            if strcmp(val,'FDK')
                res=FDK(proj,geo,alpha);
                continue
            end
            if strcmp(val,'multigrid')
                multigrid=true;
                continue
            end
            if strcmp(val,'image')
                initwithimage=1;
                continue
            end
            if isempty(res)
                error('TIGRE:SART_TV:InvalidInput','Invalid Init option')
            end
            % % % % % % % ERROR
        case 'initimg'
            if default
                continue
            end
            if exist('initwithimage','var')
                if isequal(size(val),geo.nVoxel')
                    res=single(val);
                else
                    error('TIGRE:SART_TV:InvalidInput','Invalid image for initialization');
                end
            end
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('CBCT:SART_TV:InvalidInput','Invalid quality measurement parameters');
                end
            end
        case 'tviter'
            if default
                TViter=50;
            else
                TViter=val;
            end
        case 'tvlambda'
            if default
                TVlambda=50;
            else
                TVlambda=val;
            end
        case 'orderstrategy'
            if default
                OrderStrategy='random';
            else
                OrderStrategy=val;
            end

        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:SART_TV:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in SART()']);
    end
end
if multigrid
    res=init_multigrid(proj,geo,alpha,TViter,TVlambda,gpuids);
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function k = outerProd(f, g, h)
%OUTERPRODUCT    The outer product of three CHEBFUN objects.
%
% This should be moved to chebfun/outerProduct

%   K = OUTERPRODUCT(F, G, H) returns the CHEBFUN3 representing
%   H(x,y) = F(x)G(y)H(z),
%   where F, G and H are three CHEBFUN objects.

% TODO: This command could use a compression-like algorithm, but instead
%   here we are just concatenating cols, rows and tubes.

% Empty check:
if ( isempty(f) || isempty(g) || isempty(h) )
    k = chebfun3();
    return
end

% Check they are all chebfun objects
if ( ~isa(f, 'chebfun') || ~isa(g, 'chebfun') || ~isa(h, 'chebfun'))
   error('CHEBFUN:CHEBFUN3:outerProd:badInputs', ...
       'Outer product must involve three chebfun objects.');
end

% Extract out information:
fdom = domain(f);
gdom = domain(g);
hdom = domain(h);
dom = [fdom, gdom, hdom];

k = chebfun3();
% Form outerproduct:
if ( size(f, 2) == size(g, 2) && size(f, 2) == size(h, 2) )
    k.cols = f;
    k.rows = g;
    k.tubes = h;
    k.core = zeros(size(f, 2), size(g, 2), size(h, 2));
    for i=1:size(f, 2)
        k.core(i,i,i) = 1;
    end
    k.domain = dom;
else
    error('CHEBFUN:CHEBFUN3:outerProd:sizes', ...
        'Sizes not consistent for outer product.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function obj = addJumpConstraint(obj, edge, src, tar, bounds, varargin)
    % Add jump constraints between the edge and neighboring two nodes
    % (source and target)
    %
    % These constraints include the continuity of the states/time
    % variables, as well as system-specific discrete map constraint of the
    % guard dynamics, and user-specific constraints.
    %
    % Parameters:
    % src: the source node NLP @type TrajectoryOptimization
    % edge: the edge NLP @type TrajectoryOptimization
    % tar: the target node NLP @type TrajectoryOptimization
    % bounds: the boundary values @type struct
    % varargin: extra argument @type varargin

    %% continuity of time
    t_s = SymVariable('ts',[2,1]);
    t_n = SymVariable('tn',[2,1]);
    t_cont = SymFunction('tContDomain',flatten(t_s(2)-t_n(1)),{t_s,t_n});

    % create a NlpFunction for 'edge' NLP, but use the NlpVariables from
    % 'src' and 'tar' NLPs.
    if src.Options.DistributeTimeVariable
        src_time_node = src.NumNode;
    else
        src_time_node = 1;
    end
    t_cstr = NlpFunction('Name','tContDomain',...
        'Dimension',1,...
        'lb', 0,...
        'ub', 0,...
        'Type','Linear',...
        'SymFun',t_cont,...
        'DepVariables',[src.OptVarTable.T(src_time_node);tar.OptVarTable.T(1)]);
    edge.addConstraint('tContDomain','first',t_cstr);

    %% state continuity (src <-> edge)
    x_s = src.Plant.States.x;
    x_e = SymVariable('xp',size(x_s));
    x_cont_src = SymFunction(['xMinusCont_' edge.Name],x_s-x_e,{x_s,x_e});
    x_src_cstr = NlpFunction('Name',['xMinusCont_' edge.Name],...
        'Dimension',src.Plant.numState,...
        'lb', 0,...
        'ub', 0,...
        'Type','Linear',...
        'SymFun',x_cont_src,...
        'DepVariables',[src.OptVarTable.x(end);edge.OptVarTable.x(1)]);
    edge.addConstraint('xMinusCont','first',x_src_cstr);

    %% state continuity (edge <-> tar)
    x_t = tar.Plant.States.x;
    x_e = edge.Plant.States.xn;
    x_cont_tar = SymFunction(['xPlusCont_' edge.Name],x_e-x_t,{x_e,x_t});
    x_tar_cstr = NlpFunction('Name',['xPlusCont_' edge.Name],...
        'Dimension',tar.Plant.numState,...
        'lb', 0,...
        'ub', 0,...
        'Type','Linear',...
        'SymFun',x_cont_tar,...
        'DepVariables',[edge.OptVarTable.xn(1);tar.OptVarTable.x(1)]);
    edge.addConstraint('xPlusCont','first',x_tar_cstr);

    if strcmp(edge.Plant.Type,'SecondOrder')
        %% state derivative continuity (src <-> edge)
        dx_s = src.Plant.States.dx;
        dx_e = SymVariable('xp',size(x_s));
        dx_cont_src = SymFunction(['dxMinusCont_' edge.Name],dx_s-dx_e,{dx_s,dx_e});
        dx_src_cstr = NlpFunction('Name',['dxMinusCont_' edge.Name],...
            'Dimension',src.Plant.numState,...
            'lb', 0,...
            'ub', 0,...
            'Type','Linear',...
            'SymFun',dx_cont_src,...
            'DepVariables',[src.OptVarTable.dx(end);edge.OptVarTable.dx(1)]);
        edge.addConstraint('dxMinusCont','first',dx_src_cstr);

        %% state derivative continuity (edge <-> tar)
        dx_t = tar.Plant.States.dx;
        dx_e = edge.Plant.States.dxn;
        dx_cont_tar = SymFunction(['dxPlusCont_' edge.Name],dx_e-dx_t,{dx_e,dx_t});
        dx_tar_cstr = NlpFunction('Name',['dxPlusCont_' edge.Name],...
            'Dimension',tar.Plant.numState,...
            'lb', 0,...
            'ub', 0,...
            'Type','Linear',...
            'SymFun',dx_cont_tar,...
            'DepVariables',[edge.OptVarTable.dxn(1);tar.OptVarTable.dx(1)]);
        edge.addConstraint('dxPlusCont','first',dx_tar_cstr);
    end

    %% the event function constraint for the source domain
    event_list = fieldnames(src.Plant.EventFuncs);  % all events
    if ~isempty(event_list)
        % find the index of the event associated with the edge
        event_idx = str_index(edge.Plant.EventName,event_list);
        % extract the event function object using the index
        event_obj = src.Plant.EventFuncs.(event_list{event_idx});
        % impose the NLP constraints (unilateral constraints)
        event_obj.imposeNLPConstraint(src);
        % update the upper bound at the last node to be zero (to ensure equality)
        event_cstr_name = event_obj.ConstrExpr.Name;
        updateConstrProp(src,event_cstr_name,'last','ub',0);
    end

    %% call the system constraint callback method for the discrete dyamics
    plant = edge.Plant;
    plant.UserNlpConstraint(edge, src, tar, bounds, varargin{:});
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
% vgg_selfcalib_qaffine  Upgrading projective to quasi-affine reconstruction.
%
% Given projective reconstruction [P,X] with correct signs of P and X
% (the output of vgg_signsPX_from_x), it finds homography H transforming
% [P,X] to quasi-affine reconstruction [Pq,Xq] = [P*inv(H),H*X].
% Let Ainf=[0 0 0 1] be plane at infinity, then [Pq,Xq] has the property that :-
%
%   Ainf * Xq > 0                (all scene points in front of plane at infty)
%   Ainf * vgg_wedge(Pq{k}) > 0  (all camera centers in front of plane at infty)
%
% H = vgg_selfcalib_qaffine(P,X), where
%   P ... cell(K) of double(3,4), camera matrices. K is number of cameras.
%     P also can be 3x4xK array.
%   X ... double(4,N), scene points in homog. coordinates.
%   H ... cell{I} of double(4,4), homographies upgrading [P,X] to quasi-affine reconstruction.
%     There can be 0, 1, or 2 solution classes (corresponding to I=0,1,2) :-
%       - I==0 ... no solution, ie [P,X] cannot be transformed to any affine scene.
%       - I==1 ... 1 solution, ie camera centers and scene points
%           are not separable by a plane in the true scene.
%       - I==2 ... 2 solutions, ie camera centers and scene points are separable
%           by a plane in the true scene. Then there are two solutions for plane at infinity,
%           differing by sign(det(H{i})). The two reconstruction corresponding to H{1} and H{2}
%           have oppposite handedness and we cannot say which handedness is that of the true scene.
% (Note: by 'solution' we mean rather 'class of solutions' - indeed there are infinitely many
% solutions if I>0, and linear programming chooses a single solution out of them.)
%
% EXAMPLE: Let [P,X] be a projective reconstruciton from homogeneous image points x.
% Upgrade to quasi-affine reconstruction is done as follows:
%   [P,X] = vgg_signsPX_from_x(P,X,x);
%   H = vgg_selfcalib_qaffine(P,X);
%   H = H{1}; % single solution assumed
%   P = P*inv(H);
%   X = H*X;
% If either of rows 1 and 2 returns no solution, there's something wrong with
% the reconstruction, eg an outlier.

% T.Werner, Feb 2002, werner@robots.ox.ac.uk

function H = vgg_selfcalib_qaffine(P,X)

if ndims(P)==3
  for k = 1:size(P,3)
    Q{k} = P(:,:,k);
  end
  P = Q;
end

[D N] = size(X);
K = length(P);

for k = 1:K
  C(:,k) = vgg_wedge(P{k}); % oriented camera centers
end

% Solve chiral equalities:
%
% A := found plane at infinity
% detH := required det(H)
% (A and detH can be none, one, or two according to the number of solution classes)
detH = [];
A = [];
for detHa = [-1 1]
  Aa = sephplane([X detHa*C]);
  if ~isempty(Aa)
    A = [A; Aa];
    detH = [detH; detHa];
  end
end

if isempty(A)
  H = {};
  return
end

% compose final homography H
for i = 1:size(A,1)

  % find H{i} such that H{i}(4,:)==A
  [dummy,dummy,H{i}] = svd(A(i,:),0);
  H{i} = H{i}(end:-1:1,:);

  % make det(H{i}) the same sign as detH(i)
  if det(H{i})*detH(i) < 0
    H{i} = H{i}([2 1 3 4],:);
  end

  % 'beautifier' of X:
  % Do singular value equalization on the set X,
  % i.e., make mean(nhom(H{i}*X),2)==[0;0;0] and svd(nhom(H{i}*X))==[1;1;1].
  Xi = vgg_get_nonhomg(H{i}*X);
  c = mean(Xi,2); % centroid
  Xi = Xi - c*ones(1,N);
  [U,S] = eig(Xi*Xi'); % sv equalization
  S = diag(1./sqrt(diag(S)));
  K = S*U';
  if det(K) < 0 % we want the sv equalization to be parity-preserving
    K = -K;
  end
  H{i} = [ K -K*c; 0 0 0 1 ]*H{i};

end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% A = sephplane(X)  Finds separating hyperplane A such that all(A*X)>0.
% If no solution exists, A = [].
% Works for any dimension of X.
function A = sephplane(X)

[D,N] = size(X);

X = X ./ (ones(D,1)*sqrt(sum(X.*X)));
A = [-X' ones(N,1)];
b = zeros(size(A,1),1);
f = [zeros(1,D) -1]';
LB = [-ones(1,D) 0];
UB = [ones(1,D) Inf];
fprintf('vgg selfcalib_qaffine: linprog for %d %dd pts ... ', size(X,2), D);
options = optimset('linprog');
options.Display = 'off';
[res,FVAL,EXITFLAG] = linprog(f,A,b,[],[],LB,UB, [], options);
if isempty(res)
  fprintf('no feasible plane\n');
  A = [];
  return
end
A = res(1:D)';

if ~all(A*X > 0)
  fprintf('feasible plane returned, but is not in fact feasible\n');
  A = [];
end

fprintf('Got plane [%.2f %.2f %.2f %.2f]\n', A);

return

%i = k2i(k)
% Computes indices of joint point matrix rows corresponding to views k.
%% function i = k2i(k,step)
%% k = k(:)';
%% i = [1:3]'*ones(size(k)) + 3*(ones(3,1)*k-1);
%% i = i(:);

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function [B, A] = swap(A, B)
%SWAP Swap two variables without an intermediate copy.
% Usage:
%   [B, A] = swap(A, B)

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
%  Copyright (C) 2004 Yokogawa Electric Corporation, All Rights Reserved.
%
%  usage:
%    [trigger_event] = GetMeg160TriggerEventM( fid )
%
%  arguments:
%    fid             : file ID
%
%  return values:
%    trigger_event   : n x 1 vector which include multi trigger code (n: epoch count of evoked raw data)
%
%  confirmation of revision:
%   GetMeg160TriggerEventM( Inf ) will show and return revision of this function.
%

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout=readflac(varargin)
%
% For calling details please see v_readflac.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_readflac.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: readflac.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_readflac(varargin{:});
else
    v_readflac(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function mean = extreme_values_mean ( a, b )

%*****************************************************************************80
%
%% EXTREME_VALUES_MEAN returns the mean of the Extreme Values PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a + b * euler_constant ( );

  return
end


text:
function errorbar_width(h, x, interval)
% Work with errorbar.m: Adjust the width of errorbar
%
% :Usage:
% ::
%
%    errorbar_width(h, x, interval)
%
% :Inputs:
%
%   **h:**
%        errorbar graphic handle
%
%   **x:**
%        vector x, which is used in errorbar
%
%   **interval:**
%        e.g., [-.1 .1] or [0 0]
%
% :Examples: you can see this output in
% http://wagerlab.colorado.edu/wiki/doku.php/help/core/figure_gallery
% ::
%
%    x = 1:5; % x values
%    y = [32 40 55 84 130]; % mean
%    e = [6 6 6 6 6]; % standard error of the mean
%
%    create_figure(y_axis);
%    set(gcf, 'Position', [1   512   268   194]);
%    col = [0.3333    0.6588    1.0000];
%    markercol = col-.2;
%
%    h = errorbar(x, y, e, 'o', 'color', 'k', 'linewidth', 1.5, 'markersize', 7, 'markerfacecolor', col);
%    hold on;
%    sepplot(x, y, .75, 'color', col, 'linewidth', 2);
%    errorbar_width(h, x, [0 0]); % here
%
%    set(gca, 'xlim', [.5 5.5], 'linewidth', 1.5);
%
%    try
%       pagesetup(gcf);
%       saveas(gcf, 'example.pdf');
%    catch
%       pagesetup(gcf);
%       saveas(gcf, 'example.pdf');
%    end
%
% ..
%    Copyright (C) 2014  Wani Woo
% ..

xdata = [];
for i = x
    xdata = [xdata repmat(i,1,2) NaN (repmat(i,1,2) + interval) NaN (repmat(i,1,2) + interval) NaN];
end

hh = get(h, 'children'); set(hh(2), 'XData', xdata);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function result = CDPLSA_enterFunc(numCluster,numIter,inputPath)
    %% ======================================================================
    %%STEP 1: load the data
    fprintf('start load the data...\n');
     [TrainData, TestData, TrainLabel, TestLabel, numSource, numTrain, numTarget, numTest] = CDPLSA_loadData(inputPath);
%     [TrainData, TestData, TrainLabel, TestLabel, numSource, numTrain, numTarget, numTest] = original_loadData(inputPath);
     TrainData = double(TrainData);
      TestData = double(TestData);
      TrainLabel = double(TrainLabel);
      TestLabel = double(TestLabel);
    TrainData = sparse(TrainData);
     TestData = sparse(TestData);
    fprintf('training the model and testing...\n');
    [Results, pz_d] = CD_PLSA(TrainData, TestData, TrainLabel, TestLabel, numCluster, numIter, numSource, numTrain, numTarget, numTest);

    result = Results(end);
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%%  Introduction to *image_scat*
% *image_scat* is a function to visualize scattering coefficients. It will
% display the output of every path of the scattering.
%
%% Usage
% output_image = *image_scat*(S); where *S* is the output of *scat*.
%
%% Example
% The following code will produce 3 images.
x = mandrill;
[Wop,filters] = wavelet_factory_2d(size(x));
[S,U] = scat(x,Wop);
image_scat(S);

% They correspond to the resulting scattering coefficients of each paths
% for each layer and each rotation. The legend indicates the path and the resolution.
% These filters are normalized but the normalization can be removed by
% modifying the option of image_scat. The caption can also be removed with
% a specific option.
%
% To improve legibility, the image background behind the caption is
% artificially enlightened.
%
% On the Figure 2 for instance, there is 2 different j. The first one
% corresponds to the scale of the used wavelets whereas the second one
% corresponds to the scale of the low-pass filter which is thus constant.

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
classdef form_acc < mp.form_ac
%MP.FORM_ACC  MATPOWER Formulation class for AC cartesian voltage formulations
%   Each concrete Network Model Element class must inherit, at least
%   indirectly, from both MP.NM_ELEMENT and MP.FORM.
%
%   Subclass of MP.FORM_AC.
%   MP.FORM provides properties and methods related to the specific
%   formulation (e.g. DC version, AC polar power version, etc.)
%
%   Properties
%       (model parameters inherited from MP.FORM_AC)
%
%   Methods
%       form_name() - returns string w/name of formulation ('AC-cartesian formulation')
%       form_tag() - returns string w/short label for formulation ('acc')

%   MATPOWER
%   Copyright (c) 2019-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function name = form_name(obj)
            name = 'AC-cartesian';
        end

        function tag = form_tag(obj)
            tag = 'acc';
        end

        function vtypes = model_vvars(obj)
            vtypes = {'vr', 'vi'};
        end

        function [Iu, Iw] = port_inj_current_jac(obj, ...
                n, v_, Y, M, invdiagvic, diagSlincJ)
            % [Iu, Iw] = obj.port_inj_current_jac(...)

            %% intermediate terms
            E = invdiagvic * (conj(M) - invdiagvic * diagSlincJ);

%             %% linear current term
%             Iu = Y;
%             Iw = 1j * Y;
%
%             %% + current from linear power term
%             Iu = Iu + E;
%             Iw = Iw - 1j * E;

            Iu = Y + E;
            Iw = 1j * (Y - E);
        end

        function [Iuu, Iuw, Iww] = port_inj_current_hess_v(obj, x_, lam, v_, z_, diaginvic, Y, M, diagSlincJ, dlamJ)
            % [Iuu, Iuw, Iww] = obj.port_inj_current_hess_v(x_, lam)
            % [Iuu, Iuw, Iww] = obj.port_inj_current_hess_v(x_, lam, sysx)
            % [Iuu, Iuw, Iww] = obj.port_inj_current_hess_v(x_, lam, sysx, idx)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, v_, z_, diaginvic, Y, M, diagSlincJ, dlamJ)

            if nargin < 10
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                [Y, M, N, s] = obj.get_params(idx, {'Y', 'M', 'N', 's'});
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                %% compute linear power injections
                if isempty(z_)
                    Slin = M*v_ + s;
                else
                    Slin = M*v_ + N*z_ + s;
                end

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diaginvic = sparse(1:ni, 1:ni, 1 ./ conj(vi_), ni, ni);
                if isempty(idx)     %% all ports
                    diagSlincJ = sparse(1:n, 1:n, conj(Slin), n, n);
                    dlamJ      = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    diagSlincJ = sparse(1:ni, idx, conj(Slin), ni, n);
                    dlamJ      = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            % A = diaginvic * conj(M);
            % B = diaginvic * conj(N);
            % C = diaginvic * diaginvic * diagSlincJ;
            D = (diaginvic * dlamJ).';
            E = diaginvic * (conj(M) - diaginvic * diagSlincJ);
            % E = A - C;
            F = D * E;
            G = -(F.' + F);
            % H = -D * B;

            %% linear current term
            %% second derivatives all zero

            %% current from linear power term
            Iuu = G;
            Iuw = -1j * G;
            Iww = -G;
        end

        function [Iuzr, Iuzi, Iwzr, Iwzi] = port_inj_current_hess_vz(obj, x_, lam, v_, z_, diaginvic, N, dlamJ)
            % [Iuzr, Iuzi, Iwzr, Iwzi] = obj.port_inj_current_hess_vz(x_, lam)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, sysx)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, sysx, idx)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, v_, z_, diaginvic, N, dlamJ)

            if nargin < 8
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                N = obj.get_params(idx, 'N');
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diaginvic = sparse(1:ni, 1:ni, 1 ./ conj(vi_), ni, ni);
                if isempty(idx)     %% all ports
                    dlamJ = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    dlamJ = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            % A = diaginvic * conj(M);
            B = diaginvic * conj(N);
            % C = diaginvic * diaginvic * diagSlincJ;
            D = (diaginvic * dlamJ).';
            % E = diaginvic * (conj(M) - diaginvic * diagSlincJ);
            % E = A - C;
            % F = D * E;
            % G = -(F.' + F);
            H = -D * B;

            %% current from linear power term
            Iuzr = H;
            Iuzi = -1j * H;
            Iwzr = Iuzi;
            Iwzi = -H;
        end

        function [Su, Sw] = port_inj_power_jac(obj, ...
                n, v_, Y, M, diagv, diagvi, diagIlincJ)
            % [Su, Sw] = obj.port_inj_power_jac(...)

            %% intermediate terms
%             A = diagIlincJ;
            B = diagvi * conj(Y);

%             %% linear power term
%             Su = M;
%             Sw = 1j * M;
%
%             %% + power from linear current term
%             Su = Su + A + B;
%             Sw = Sw + 1j * (A - B);

            A = M + diagIlincJ;
            Su = A + B;
            Sw = 1j * (A - B);
        end

        function [Suu, Suw, Sww] = port_inj_power_hess_v(obj, x_, lam, v_, z_, diagvi, Y, M, diagIlincJ, dlamJ)
            % [Suu, Suw, Sww] = obj.port_inj_power_hess_v(x_, lam)
            % [Suu, Suw, Sww] = obj.port_inj_power_hess_v(x_, lam, sysx)
            % [Suu, Suw, Sww] = obj.port_inj_power_hess_v(x_, lam, sysx, idx)
            % [...] = obj.port_inj_power_hess_v(x_, lam, v_, z_, diagvi, Y, M, diagIlincJ, dlamJ)

            if nargin < 10
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                [Y, L, M, i] = obj.get_params(idx, {'Y', 'L', 'M', 'i'});
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                %% compute linear current injections
                if isempty(z_)
                    Ilin = Y*v_ + i;
                else
                    Ilin = Y*v_ + L*z_ + i;
                end

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diagvi = sparse(1:ni, 1:ni, vi_, ni, ni);
                if isempty(idx)     %% all ports
                    diagIlincJ = sparse(1:n, 1:n, conj(Ilin), n, n);
                    dlamJ      = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    diagIlincJ = sparse(1:ni, idx, conj(Ilin), ni, n);
                    dlamJ      = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            % D = dlamJ.';
            E = dlamJ.' * conj(Y);  %% D * conj(Y);
            % F = E + E.';
            % G = 1j * (E - E.');

            %% linear power term
            %% second derivatives all zero

            %% power from linear current term
            Suu = E + E.';          %% F
            Suw = 1j * (E.' - E);   %% G.'
            Sww = Suu;
        end

        function [Suzr, Suzi, Swzr, Swzi] = port_inj_power_hess_vz(obj, x_, lam, v_, z_, diagvi, L, dlamJ)
            % [Suzr, Suzi, Swzr, Swzi] = obj.port_inj_power_hess_vz(x_, lam)
            % [...] = obj.port_inj_power_hess_vz(x_, lam, sysx)
            % [...] = obj.port_inj_power_hess_vz(x_, lam, sysx, idx)
            % [...] = obj.port_inj_power_hess_vz(x_, lam, v_, z_, diagvi, L, dlamJ)

            if nargin < 8
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                L = obj.get_params(idx, 'L');
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diagvi = sparse(1:ni, 1:ni, vi_, ni, ni);
                if isempty(idx)     %% all ports
                    dlamJ = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    dlamJ = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            D = dlamJ.';
            H = D * conj(L);

            %% power from linear current term
            Suzr = H;
            Suzi = -1j * H;
            Swzr = 1j * H;
            Swzi = H;
        end

        function [va, vm] = aux_data_va_vm(obj, ad)
            v_ = ad.vr + 1j * ad.vi;
            va = angle(v_);
            vm = abs(v_);
        end

        function [g, dg] = va_fcn(obj, xx, idx, lim)
            %% lim can be a vector value for equality constraint or
            %% upper bound on va, or a cell array with {vamin, vamax}
            %% for double bounds

            %% unpack data
            [vr, vi] = deal(xx{:});

            %% compute voltage angle mismatch
            if isempty(idx)
                va = angle(vr + 1j * vi);
            else
                va = angle(vr(idx) + 1j * vi(idx));
            end
            if iscell(lim)
                g = [ lim{1} - va;
                      va - lim{2}   ];
            else
                g = va - lim;
            end

            if nargout > 1
                %% compute partials of voltage angle w.r.t vr and vi
                nn = length(vr);
                if isempty(idx)
                    idx = 1:nn;
                end
                n = length(idx);
                vm2 = vr(idx).^2 + vi(idx).^2;
                dva_dvr = sparse(1:n, idx, -vi(idx) ./ vm2, n, nn);
                dva_dvi = sparse(1:n, idx,  vr(idx) ./ vm2, n, nn);
                if iscell(lim)
                    dg = [ -dva_dvr -dva_dvi;       %% va w.r.t vr, vi
                            dva_dvr  dva_dvi  ];    %% va w.r.t vr, vi
                else
                    dg = [dva_dvr dva_dvi];         %% va w.r.t vr, vi
                end
            end
        end

        function d2G = va_hess(obj, xx, lam, idx)
            %% unpack data
            [vr, vi] = deal(xx{:});
            nn = length(vr);

            %% evaluate Hessian of voltage angle function
            if isempty(idx)
                vvr = vr;
                vvi = vi;
                idx = 1:nn;
            else
                vvr = vr(idx);
                vvi = vi(idx);
            end
            vvr2 = vvr.^2;
            vvi2 = vvi.^2;
            vvm4 = (vvr2 + vvi2).^2;
            n = length(idx);
            if length(lam) == n     %% upper bound or equality
                lamvm4 = lam ./ vvm4;
            else                    %% doubly bounded (use lam_ub-lam_lb)
                lamvm4 = (lam(n+1:2*n) - lam(1:n)) ./ vvm4;
            end
            d2vref_rr = sparse(idx, idx, 2 * lamvm4 .*  vvr .* vvi,   nn, nn);
            d2vref_ri = sparse(idx, idx,     lamvm4 .* (vvi2 - vvr2), nn, nn);

            %% construct Hessian
            d2G = [ d2vref_rr  d2vref_ri;
                    d2vref_ri -d2vref_rr ];
        end

        function [g, dg] = vm2_fcn(obj, xx, idx, lim)
            %% lim can be a scalar value for equality constraint or
            %% upper bound on vm^2, or a cell array with {vm2min, vm2max}
            %% for double bounds

            %% unpack data
            [vr, vi] = deal(xx{:});

            %% compute voltage magnitude^2 mismatch
            if isempty(idx)
                vm2 = vr.^2 + vi.^2;
            else
                vm2 = vr(idx).^2 + vi(idx).^2;
            end
            if iscell(lim)
                g = [ lim{1} - vm2;
                      vm2 - lim{2}   ];
            else
                g = vm2 - lim;
            end

            if nargout > 1
                %% compute partials of voltage magnitude^2 w.r.t vr and vi
                nn = length(vr);
                if isempty(idx)
                    idx = 1:nn;
                end
                n = length(idx);
                dvm_dvr = sparse(1:n, idx, 2 * vr(idx), n, nn);
                dvm_dvi = sparse(1:n, idx, 2 * vi(idx), n, nn);
                if iscell(lim)
                    dg = [ -dvm_dvr -dvm_dvi;
                            dvm_dvr  dvm_dvi  ];    %% vm2 w.r.t vr, vi
                else
                    dg = [ dvm_dvr dvm_dvi ];       %% vm2 w.r.t vr, vi
                end
            end
        end

        function d2G = vm2_hess(obj, xx, lam, idx)
            %% unpack data
            [vr, vi] = deal(xx{:});
            nn = length(vr);

            %% evaluate Hessian of voltage magnitude^2 function
            if isempty(idx)
                idx = 1:nn;
            end
            n = length(idx);
            if length(lam) == n     %% upper bound or equality
                dlam = sparse(idx, idx, 2*lam, nn, nn);
            else                    %% doubly bounded (use lam_ub-lam_lb)
                dlam = sparse(idx, idx, 2*(lam(n+1:2*n) - lam(1:n)), nn, nn);
            end

            %% construct Hessian
            zz = sparse(nn, nn);
            d2G = [dlam zz; zz dlam];
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef humanize
    % humanize functions
    % collection of static functions to convert si units into 'human readable'
    % strings. Currently implemented functions are
    %
    % humanize.seconds      % creates a time string
    % humanize.bytes        % creates a data size string
    % humanize.clock
    %
    % written by T.Sumpf (tsumpf@gwdg.de) June 2010

    methods (Static)
        function str = bytes(nBytes, precision)
            if nargin < 2
                precision = '%1.3f';
            end
            order = log2(nBytes);
            order = floor(order/10) * 10;

            unitValue = nBytes / 2^order;

            switch(order)
                case 0
                    str = [num2str(nBytes),' byte'];
                case 10
                    str = [num2str(unitValue,precision),' kB'];
                case 20
                    str = [num2str(unitValue,precision),' MB'];
                case 30
                    str = [num2str(unitValue,precision),' GB'];
                otherwise
                    str = [num2str( nBytes / 2^40,precision),' TB'];
            end

        end

        function str = seconds(secs, precision)
            if nargin < 2
                precision = '%1.3f';
            end

            % decimal system
            if secs < 1e-3
                str = [num2str(secs*1e6,precision),' us'];
                return
            end

            if secs < 1
                str = [num2str(secs*1e3,precision),' ms'];
                return
            end

            if secs < 60
                str = [num2str(secs,precision),'s'];
                return
            end

            % ...else
            % "date system"

            remSecs = floor(mod(secs,60));
            mins    = mod(floor(secs/60),60);
            hours   = mod(floor(secs/(60*60)),24);
            days    = floor(secs/(60*60*24));

            str = [num2str(remSecs),'s'];

            if secs >= 60
                str = [num2str(mins),'min : ',str];
            end

            if secs >= (60 * 60)
                str = [num2str(hours),'h : ',str];
            end

            if secs >= (60 * 60 * 24)
                str = [num2str(days),'d : ',str];
            end

        end

        function str = clock(clk, date_time_or_both)
            if nargin < 1
                date_time_or_both = 'both';
            end

            c = fix(clk);

            d = sprintf('%d.%d.%d',c(3),c(2),c(1));
            t = sprintf('%d:%d:%d',c(4),c(5),c(6));

            switch lower(date_time_or_both)
                case 'date'
                    str = d;
                case 'time'
                    str = t;
                case 'both'
                    str = [d, ' ' , t];
            end

        end
    end
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function mask = maskSphere( d, r, show )
% Creates an 'image' of a n-dimensional hypersphere.
%
% Useful for testing visualization procedures.
%
% Create a d-dimensional matrix mask of dimensions [s x s x s x ... x s]
% where s=2r+1.  Each element mask(x1,x2...,xd) is 1 if (x1,...xd) falls
% inside the centered hypersphere of rad r.
%
% In 1d this corresponds to a vector of the form [.. 0 1 1 1 1 1 0 ..]
% In 2d this corresponds to an image of a white circle.
% In 3d and 4d, well try it and 'see' what it looks like.
%
% USAGE
%  mask = maskSphere( d, r, [show] )
%
% INPUTS
%  d         - dimension (any positive integer)
%  r         - sphere integer radius
%  show      - [] figure in which to display results
%
% OUTPUTS
%  mask      - [s x s x s x ... x s] hypersphere image
%
% EXAMPLE
%  mcircle = maskSphere( 2, 20, 1 );
%  msphere = maskSphere( 3, 10, 2 );
%  msphere = maskSphere( 4, 10, 3 );
%
% See also MASKCIRCLE, MONTAGE2
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(show) ); show = 1; end

xs=cell(1,d);
for i=1:d; xs{i}=-r:r; end
if( d>1 ); [xs{:}] = ndgrid(xs{:}); else xs{1}=xs{1}'; end
mask=xs{1}.^2;
for i=2:d; mask=mask+xs{i}.^2; end
mask = double( mask < (r+1)^2 );

if( show )
  figure(show); clf;
  if( d<=2 )
    im( mask );
  elseif any(d==3:4)
    montage2( mask );
  else
    disp('no visualization available for d>4');
  end
end

author: zouchuhang

repo: LayoutNet


text:
function [C,P]=knn(d, Cp, K)

%KNN K-Nearest Neighbor classifier using an arbitrary distance matrix
%
%  [C,P]=knn(d, Cp, [K])
%
%  Input and output arguments ([]'s are optional):
%   d     (matrix) of size NxP: This is a precalculated dissimilarity (distance matrix).
%           P is the number of prototype vectors and N is the number of data vectors
%           That is, d(i,j) is the distance between data item i and prototype j.
%   Cp    (vector) of size Px1 that contains integer class labels. Cp(j) is the class of
%            jth prototype.
%   [K]   (scalar) the maximum K in K-NN classifier, default is 1
%   C     (matrix) of size NxK: integers indicating the class
%           decision for data items according to the K-NN rule for each K.
%           C(i,K) is the classification for data item i using the K-NN rule
%   P     (matrix) of size NxkxK: the relative amount of prototypes of
%           each class among the K closest prototypes for each classifiee.
%           That is, P(i,j,K) is the relative amount of prototypes of class j
%           among K nearest prototypes for data item i.
%
% If there is a tie between representatives of two or more classes
% among the K closest neighbors to the classifiee, the class i selected randomly
% among these candidates.
%
% IMPORTANT  If K>1 this function uses 'sort' which is considerably slower than
%            'max' which is used for K=1. If K>1 the knn always calculates
%            results for all K-NN models from 1-NN up to K-NN.
%
% EXAMPLE 1
%
% sP;                           % a SOM Toolbox data struct containing labeled prototype vectors
% [Cp,label]=som_label2num(sP); % get integer class labels for prototype vectors
% sD;                           % a SOM Toolbox data struct containing vectors to be classified
% d=som_eucdist2(sD,sP);        % calculate euclidean distance matrix
% class=knn(d,Cp,10);           % classify using 1,2,...,10-rules
% class(:,5);                   % includes results for 5NN
% label(class(:,5))             % original class labels for 5NN
%
% EXAMPLE 2 (leave-one-out-crossvalidate KNN for selection of proper K)
%
% P;                          % a data matrix of prototype vectors (rows)
% Cp;                         % column vector of integer class labels for vectors in P
% d=som_eucdist2(P,P);        % calculate euclidean distance matrix PxP
% d(eye(size(d))==1)=NaN;     % set self-dissimilarity to NaN:
%                             % this drops the prototype itself away from its neighborhood
%                             % leave-one-out-crossvalidation (LOOCV)
% class=knn(d,Cp,size(P,1));  % classify using all possible K
%                             % calculate and plot LOOC-validated errors for all K
% failratep = ...
%  100*sum((class~=repmat(Cp,1,size(P,1))))./size(P,1); plot(1:size(P,1),failratep)

% See also SOM_LABEL2NUM, SOM_EUCDIST2, PDIST.
%
% Contributed to SOM Toolbox 2.0, October 29th, 2000 by Johan Himberg
% Copyright (c) by Johan Himberg
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta Johan 291000

%% Init %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check K
if nargin<3 || isempty(K),
  K=1;
end

if ~vis_valuetype(K,{'1x1'})
  error('Value for K must be a scalar');
end

% Check that dist is a matrix
if ~vis_valuetype(d,{'nxm'}),
  error('Distance matrix not valid.')
end

[N_data N_proto]=size(d);

% Check class label vector: must be numerical and of integers
if ~vis_valuetype(Cp,{[N_proto 1]});
  error(['Class vector is invalid: has to be a N-of-data_rows x 1' ...
	 ' vector of integers']);
elseif sum(fix(Cp)-Cp)~=0
  error('Class labels in vector ''Cp'' must be integers.');
end

if size(d,2) ~= length(Cp),
  error('Distance matrix and prototype class vector dimensions do not match.');
end

% Check if the classes are given as labels (no class input arg.)
% if they are take them from prototype struct

% Find all class labels
ClassIndex=unique(Cp);
N_class=length(ClassIndex); % number of different classes

%%%% Classification %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if K==1,   % sort distances only if K>1

  % 1NN
  % Select the closest prototype
  [~,proto_index]=min(d,[],2);
  C=Cp(proto_index);

else

  % Sort the prototypes for each classifiee according to distance
  [~, proto_index]=sort(d);

  %% Select up to K closest prototypes
  proto_index=proto_index(1:K,:);
  knn_class=Cp(proto_index);
  for i=1:N_class,
    classcounter(:,:,i)=cumsum(knn_class==ClassIndex(i));
  end

  %% Vote between classes of K neighbors
  [winner,vote_index]=max(classcounter,[],3);

  %%% Handle ties

  % Set index to classes that got as much votes as winner

  equal_to_winner=(repmat(winner,[1 1 N_class])==classcounter);

  % set index to ties
  [tie_indexi,tie_indexj]=find(sum(equal_to_winner,3)>1); % drop the winner from counter

  % Go through tie cases and reset vote_index randomly to one
  % of them

  for i=1:length(tie_indexi),
    tie_class_index=find(squeeze(equal_to_winner(tie_indexi(i),tie_indexj(i),:)));
    fortuna=randperm(length(tie_class_index));
    vote_index(tie_indexi(i),tie_indexj(i))=tie_class_index(fortuna(1));
  end

  C=ClassIndex(vote_index)';
end

%% Build output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Relative amount of classes in K neighbors for each classifiee

if K==1,
  P=zeros(N_data,N_class);
  if nargout>1,
    for i=1:N_data,
      P(i,ClassIndex==C(i))=1;
    end
  end
else
  P=shiftdim(classcounter,1)./repmat(shiftdim(1:K,-1), [N_data N_class 1]);
end

author: ilarinieminen

repo: SOM-Toolbox


text:
clear all; close all;
I=imread('peppers.png');
J=rgb2gray(I);
K=fft2(J);
K=fftshift(K);
L=abs(K/256);
figure;
subplot(121);
imshow(J);
subplot(122);
imshow(uint8(L));

author: HuangCongQing


text:
function value = monomial_value ( dim_num, point_num, x, expon )

%*****************************************************************************80
%
%% MONOMIAL_VALUE evaluates a monomial.
%
%  Discussion:
%
%    This routine evaluates a monomial of the form
%
%      product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
%
%    where the exponents are nonnegative integers.  Note that
%    if the combination 0^0 is encountered, it should be treated
%    as 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points at which the
%    monomial is to be evaluated.
%
%    Input, real X(DIM_NUM,POINT_NUM), the point coordinates.
%
%    Input, integer EXPON(DIM_NUM), the exponents.
%
%    Output, real VALUE(POINT_NUM), the value of the monomial.
%
  value(1,1:point_num) = 1.0;

  for dim = 1 : dim_num
    if ( 0 ~= expon(dim) )
      value(1,1:point_num) = value(1,1:point_num) .* ( x(dim,1:point_num).^expon(dim) );
    end
  end

  return
end


text:
function [Model, Info] = linear_sparse_stepwise_vec(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     Automatic Relevance Prior for each input dimension
%     is imposed to get sparse weight matrix
%
% Notice !!
%     Delay embedding is done in this module,
%     then input vector should be original input data without embedding
%
%   [Model, Info] = linear_sparse_stepwise(X,Y,Model,parm)
%     Scalar iteration version
% --- Input
%  Y  : Output data ( N x T x Ntrial )
%  X  : Input data  ( M x (T + (Dtau-1)*Tau)  x Ntrial)
%  N  =  # of output
%  M  =  # of input (original input space dimension)
%  T  =  # of time sample
%
%  Estimate the following model
%    Y(t) = W * [X(:, t + (Dtau-1)*Tau ); ...; X(:, t)]
%
%  Model : Structure for estimated model
%  if Model is empty, initialization is done before training
%  if Model is previous training result, re-training us done
%
%  parm  : Structure for learning parameter
%  parm.Npre_train :  # of VB-update in initial training
%  parm.Ntrain :  # of training
%  parm.Nskip  :  skip # for print
%  parm.a_min  :  Min value for pruning small variance component
%  parm.Prune  :  = 1 : Prune small variance & irrelevant input dimension
%
%  parm.Tau       = Lag time
%  parm.Dtau      = Number of embedding dimension
%    Total input dimension in embedding space is (M * parm.Dtau)
% --- Output
%  Model : Structure for estimated model
%  Model.SY  :  Noise variance         ( 1 x 1 )
%  Model.W   :  Weight matrix          ( N x M*D ) , D = parm.Dtau
%  Model.A   :  Prior weight variance  ( 1 x M*D )
%  Model.ix_act : Active index for W after pruning
%
%  Info  : Structure for learning process history
%  Info.FE  = LP + H : Free energy
%  Info.LP  = Log likelihood
%  Info.H   = - Model entropy
%
% 2009-10-18 Made by M. Sato

MINVAL = 1.0e-15;

fprintf('linear_sparse_stepwise start\n')

% Dimension
[N ,T  ,Ntrial]  = size(Y); % N = # of output
[M ,Tx ,Ntrialx] = size(X); % M = # of input without embedding

Tall = T * Ntrial;

if Ntrial~=Ntrialx, error('# of trial is different in X and Y'); end;

Nskip  = 100;   % skip steps for display info
a_min = 1e-14;	% Minimum value for weight pruning
Fdiff = 1e-12; % Threshold for convergence
Ncheck = 100;   % Minimum number of training iteration
Fstep  = 5;     % Free energy convergence check step
Prune  = 1;     % Prune mode
space_ARD=1;	% if space_ARD=1, ARD is done only for space dimension

if isfield(parm,'Nskip'), Nskip  = parm.Nskip; end;
if isfield(parm,'Fdiff'), Fdiff   = parm.Fdiff; end;
if isfield(parm,'a_min'), a_min   = parm.a_min ; end;
if isfield(parm,'Prune'), Prune = parm.Prune; end;
if isfield(parm,'Ncheck'), Ncheck = parm.Ncheck; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'space_ARD'),space_ARD  = parm.space_ARD ; end

% # of embedding dimension
if isfield(parm,'Dtau')
	D    = parm.Dtau;
	tau  = parm.Tau;
else
	D    = 1;
	tau  = 1;
end
% Parameters
Ntrain = parm.Ntrain;

if Ntrain < 1, Info = []; return; end;

if isfield(parm,'Npre_train')
	Npre_train = parm.Npre_train;
else
%	Npre_train = Ntrain;
	if Tall >= 2*M*D
		Npre_train = 0;
	elseif Tall >= M*D
		Npre_train = fix(Ntrain/2);
	else
		Npre_train = Ntrain;
	end
end

if Npre_train > Ntrain, Npre_train = Ntrain; end;

if isfield(parm,'Nupdate')
	Nupdate = parm.Nupdate;
else
	Nupdate = 1;
end

fprintf('--- Output Dimension  = %d\n',N)
fprintf('--- Input  Dimension  = %d\n',M)
fprintf('--- Embedding  Dimension  = %d\n',D)
fprintf('--- Number of trials  = %d\n',Ntrial)
fprintf('--- Number of training sample = %d\n',Tall)
fprintf('--- Total update iteration    = %d (%d)\n',Ntrain,Npre_train)

% original input dimension
Xdim  = M;
M = Xdim*D;
M_ALL = M;

%
% --- Initialization
%

% Input/Output variance
sx = mean((X(:) - mean(X(:))).^2);
sy = mean((Y(:) - mean(Y(:))).^2);

SY0 = mean(sy);
A0  = 1./mean(sx);

if isfield(parm,'Ta0') && parm.Ta0 > 0,
	Ta0 = parm.Ta0;
	a0  = parm.a0 * A0;
else
	Ta0 = 0;
	a0  = 1;
end

% Input covariance
XX  = sum(sum(X.^2,3),2)/(Tx*Ntrial);
XX  = repmat(XX', [1 D]);% 1 x M

if isfield(Model,'ix_act')
	SY  = Model.SY;  % 1 x 1

	% Active index
	ix_act = Model.ix_act;
	A = zeros(1,M_ALL);
	W = zeros(N,M_ALL);

	A(ix_act)   = sum(Model.A,1);
	W(:,ix_act) = Model.W;

	if M_ALL ~= Xdim*D,
		fprintf('M_ALL=%d,Xdim=%d,D=%d\n',M_ALL ,Xdim ,D)
		error('M_ALL ~= Xdim*D')
	end

	% Active index in input space without embedding
	IX_dim = find( sum(reshape(A,[Xdim,D]),2) > 0);

	IX_act = repmat( (0:(D-1))* Xdim ,[length(IX_dim) 1]) ...
	       + repmat(IX_dim, [1 D]);
	IX_act = IX_act(:);

	W  = W(:,IX_act) ;  % N x (M*D)
	A  = A(IX_act) ;    % 1 x (M*D)
	X  = X(IX_dim,:,:);
	XX = XX(IX_act);

	M     = length(IX_act);
	Xdim  = length(IX_dim);
else
	A = repmat(A0, [1,M_ALL]);
	W = zeros(N,M_ALL);
	SY = SY0;
	% Save original input for pruning
	IX_act  = 1:M_ALL;
	IX_dim  = 1:Xdim;
end

ix_act  = 1:M;
ix_dim  = 1:Xdim;
M_all   = M;

% Delay embedding index for W
Wid = [(0:(D-1))'* Xdim + 1 , (1:D)'* Xdim];

% Working variable
G_A = zeros(1,M);       % N x M

if size(A,1)==Xdim && size(A,2)==1
	A   = repmat(A',[1 D]);
elseif size(A,1)==1 && size(A,2)==Xdim
	A   = repmat(A,[1 D]);
end

if length(SY)==1, SY= repmat(SY,N,1); end;

SW = 1./(Tall*XX + 1./A );

fprintf('a_min = %g\n', a_min)
fprintf('SY0 = %g\n', SY0)
fprintf('SY  = %g\n', mean(SY))

% Free energy histry
FE  = zeros(Ntrain,1);
LP  = zeros(Ntrain,1);
H   = zeros(Ntrain,1);
Err = zeros(Ntrain,1);
Mhist = zeros(Ntrain,1);

% ARD hyper param. history
if isfield(parm,'Debug') && ~isempty(parm.Debug) && parm.Debug > 0
	Debug = 1;
	A_tmp = zeros(N*M_ALL, ceil(Ntrain/Nskip));
else
	Debug = 0;
end

% recover all component
A_all  = zeros(1,M_all);       % N x M

[dY] = error_delay_time(X,Y,W,tau);

k_save  = 0;

%%%%%% Learning Loop %%%%%%
for k=1:Ntrain
	% Ainv = alpha , A = 1/alpha
	Ainv = 1./A;

	% E = ( (Y-W*X)^2 +  W^2 * Ainv )/SY
    [W]  = weight_update_embed(X, dY, W, Tall*XX, Ainv, tau);
	[dY] = error_delay_time(X, Y, W, tau);

    dYY = sum(dY.^2,2)/(Tall); % N x 1
	WW  = W.^2;

    % Noise variance update
    SY  = dYY + WW * Ainv'/(Tall);
    % Prevent zero variance
    SY  = max( SY, MINVAL);

    % Weight variance
	SW = 1./(Tall*XX + Ainv );

    % Log variance
    SWA     = max( SW .* Ainv , MINVAL);
    log_sw  = N*(sum( log(SWA) - SWA + 1 ));
    log_sy  = sum( log(SY) );
    log_a   = Ta0*sum( log(Ainv) - a0.*Ainv + 1 );

    % Free energy
    LP(k)  = - (0.5) * (log_sy + N*sum(SW.*XX));
    H(k)   = 0.5*( log_sw + log_a  );
    FE(k)  = LP(k) + H(k)/Tall;
    Err(k) = sum(dYY)/(N*SY0);

	% E = ( (Y-W*X)^2 +  W^2 * Ainv )/SY

    % Hyper parameter for weight variance (ARD)
	if mod(k,Nupdate) == 0,
		% Ainv = 1./A;
	    % SW = 1./( Tall*XX + Ainv );

	    % G_A = 1 - (SW)./A;
	    %     = (Tall*XX + Ainv - 1/A) * SW
	    G_A = Tall.* ( SW ) .* XX;
	    G_A = max((G_A), MINVAL);

		% N*A = ( (W.^2) + SW )/SY ;
		if k <= Npre_train,
			% VB update rule (Stable)
			%  ARD for each weight
			A  = ((1./SY)'*WW + N*SW + 2*Ta0*a0)/( N + 2*Ta0 );
		else
			% Accelerated update rule (Unstable)
		    A  = sqrt( A.* ((1./SY)'*WW) ./(G_A * N) );
		    %A  = (WW + 2*Ta0*a0)./(G_A * N + 2*Ta0);
		end

		if space_ARD==1
			Am = mean(reshape(A,[M/D,D]),2);
			A  = repmat(Am', 1,D);
		end

	    % Prune small variance
	    if Prune == 1
		    % Find active input dimension
		    ix_act_old = ix_act;
		    ix_dim_old = ix_dim;

		    % recover all component
			A_all  = zeros(1,M_all);       % 1 x M
		    %  A_all(:,ix_act)  = A;
		    WWW = sum(WW,1);
		    A_all(ix_act) = WWW/max(WWW);

		    % find active input dimension (absolute index)
		    ix_dim = find( sum(reshape(A_all,[Xdim,D]),2) > a_min );
		    ix_act = repmat(ix_dim, [1 D]) ...
		           + repmat([0:D-1]*Xdim, [length(ix_dim) 1]);
		    ix_act = ix_act(:);

		    Mnew   = length(ix_act);  		% # of effective input

		    if Mnew < M,
			    % convert to relative index
			    jx_act = trans_index(ix_act,ix_act_old,M_all);
			    jx_dim = trans_index(ix_dim,ix_dim_old,Xdim);

			    M   = Mnew;
			    A   = A(jx_act) ;    % 1 x M
			    W   = W(:,jx_act) ;  % N x M
			    SW  = SW(:,jx_act);  % N x M

			    X	 = X(jx_dim,:,:);  	 	% Xdim x T
			    XX	 = XX(jx_act);  	 	% 1 x M
			end
	    end
	    % END of if Prune == 1

		A = max(A,MINVAL);

	end
	% END of if mod(k,Nupdate) == 0

	Mhist(k) = M;

    if mod(k, Nskip)==0
        % Save history
		if Debug == 1
        	k_save = k_save + 1;
        	A_tmp(:,k_save) = A_all(:);
		end

        fprintf('Iter = %4d, M = %4d, err = %g, F = %g, H = %g\n', ...
               k, length(ix_act), Err(k), FE(k), - 2*H(k)/log(T));
    end

    % Convergence check
	if k > Ncheck,
		Fdif = (FE(k) - FE(k-Fstep))/(abs(FE(k))+eps);
	else
		Fdif = Fdiff + 1;
	end

	if (Fdiff > abs(Fdif)),
		fprintf('Converged : Free energy change = %g\n',Fdif)
		break;
	end;
end

ix_act = IX_act(ix_act);

% Active index
Model.ix_act = ix_act;
Model.M_all  = M_ALL ;

% Save output variable
Model.A    = A ;
Model.W    = W ;
Model.SY   = SY;
Model.SW   = SW; % = 1./(Tall*XX + diag(Ainv))

Model.method = 'linear_sparse_stepwise';
Model.mode   = 'scalar';
Model.sparse = 'sparse';

% Save history
Info.FE  = FE(1:k);
Info.LP  = LP(1:k);
Info.H   = H(1:k) ;
Info.Err = Err(1:k);
Info.M   = Mhist(1:k);

if exist('A_tmp','var')
	Info.A   = A_tmp(:,1:k_save) ;
end

% recover all component
%W_all   = zeros(N,M_all);
%SW_all  = zeros(N,M_all);
%
%A_all(:,ix_act)  = A;
%W_all(:,ix_act)  = W ;  % N x M
%SW_all(:,ix_act) = SW;  % N x M

%%% ---- Index transformation from old active_index to current active_index
function	jx = trans_index(ix,ix_old,M)

N = length(ix_old);
Itrans = zeros(M,1);
Itrans(ix_old) = 1:N;

jx = Itrans(ix);

author: KamitaniLab

repo: GenericObjectDecoding


text:
% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix

try
	Type = ver;
	% This part is for OCTAVE
	if(strcmp(Type(1).Name, 'Octave') == 1)
		mex libsvmread.c
		mex libsvmwrite.c
		mex svmtrain.c ../svm.cpp svm_model_matlab.c
		mex svmpredict.c ../svm.cpp svm_model_matlab.c
	% This part is for MATLAB
	% Add -largeArrayDims on 64-bit machines of MATLAB
	else
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims svmtrain.c ../svm.cpp svm_model_matlab.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims svmpredict.c ../svm.cpp svm_model_matlab.c
	end
catch
	fprintf('If make.m failes, please check README about detailed instructions.\n');
end

author: rbgirshick

repo: voc-dpm

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


text:
function [sFile, nEpochs] = in_epochs_nwb(sFile, nwb2)
% IN_TRIALS_NWB Read trials from a Neurodata Without Borders .nwb file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Konstantinos Nasiotis, 2019,2020

% Check if there are trials present
if isempty(nwb2.intervals_epochs)
    nEpochs = 1;
    sFile.prop.nAvg = 1;
    sFile.format    = 'NWB-CONTINUOUS';
    return

else
    % Get epochs
    timeBoundsTrials = double([nwb2.intervals_epochs.start_time.data.load nwb2.intervals_epochs.stop_time.data.load]);

    % Get number of epochs
    nEpochs = length(nwb2.intervals_epochs.start_time.data.load);

    % === EPOCHS FILE ===
    % Build epochs structure
    for iEpoch = 1:nEpochs
        try
            epochLabel = nwb2.intervals_epochs.tags.data.load(iEpoch);
            sFile.epochs(iEpoch).label       = ['Epoch #(' num2str(iEpoch) ') - ' epochLabel{1}];
        catch
            sFile.epochs(iEpoch).label       = ['Epoch #(' num2str(iEpoch) ')'];
        end
        sFile.epochs(iEpoch).times       = timeBoundsTrials(iEpoch,:);
%         sFile.epochs(iEpoch).samples     = round(sFile.epochs(iEpoch).times * sFile.prop.sfreq);
        sFile.epochs(iEpoch).nAvg        = 1;
        sFile.epochs(iEpoch).select      = 1;
        sFile.epochs(iEpoch).bad         = 0;
%         sFile.epochs(iEpoch).bad         = badTrials(iEpoch);
        sFile.epochs(iEpoch).channelflag = [];
    end

    sFile.format    = 'NWB';
end

end

author: brainstorm-tools


text:
function [ normalisation_options ] = Parse_settings( sigma, patch_type, ratio_neg, num_samples, varargin)
%PARSE_SETTINGS Summary of this function goes here
%   Detailed explanation goes here

    % creating the parameters to use when training colour (intensity) patches
    normalisation_options = struct;

    % this is what currently is expected (although could potentially have
    % bigger or smaller support regions
    normalisation_options.patchSize = [11 11];

    % The region size of a region that is taken for training around an
    % aligned or misaligned landmark
    if(sum(strcmp(varargin,'normalisation_size')))
        ind = find(strcmp(varargin,'normalisation_size')) + 1;
        normalisation_options.normalisationRegion = [varargin{ind}, varargin{ind}];
    else
       normalisation_options.normalisationRegion = [21 21];
    end

    % This specifies the split of data ratios
    normalisation_options.svmRatio = 0.8; % proportion of data used for training SVR
    normalisation_options.logitRatio = 0.1; % proportion of data for training logistic regressors
    % the rest is used for testing and provides the correlation and rms scores

    % should normalised cross correlation or just cross correlation should
    % be used on the patch as an SVR
    normalisation_options.useNormalisedCrossCorr = 1;

    % the patch types to be used (for now 'reg' (raw pixel values), and
    % 'grad' (gradient intensity values)
    normalisation_options.patch_type = patch_type;

    % number of training samples to use
    normalisation_options.numSamples = num_samples;

    normalisation_options.sigma = sigma;

    normalisation_options.rate_negative = ratio_neg;

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [x, infos] = sep_symm_nmtf(V, rank, in_options)
% Separable symmetric nonnegative matrix tri-factorization (Sep-Symm-NMTF)
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       Arora, Ge, Halpern, Mimno, Moitra, Sontag, Wu, Zhu,
%       "A practical algorithm for topic modeling with provable guarantees,"
%       International Conference on Machine Learning (ICML), pp. 280-288,
%       2013.
%
%
% This file is part of NMFLibrary.
%
% This file has been ported from
%       septrisymNMF.m at https://gitlab.com/ngillis/nmfbook/-/tree/master/algorithms
%       by Nicolas Gillis (nicolas.gillis@umons.ac.be)
%
% Change log:
%
%       June 14, 2022 (Hiroyuki Kasai): Ported initial version
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.delta = 1e-6;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    method_name = 'Sep-Symm-NMTF';
    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % idendity K such that V(K,K) = W(K,:) S W(K,:)^T
    spa_options.normalize = 1;
    spa_sol = spa(V, rank, spa_options);

    % solve V(K,K) z = q = V(K,:)
    q = V(:, spa_sol.K)' * ones(m,1);
    % y = A(K,K)\q; % This works in noiseless conditions
    options.alg_name = method_name;
    [y, WtW, WtV, infos] = nnls_solver(q, V(spa_sol.K, spa_sol.K), options);
    fprintf('\n');

    % recover S and W
    S = diag(y) * V(spa_sol.K, spa_sol.K) * diag(y);
    % W = V(:,K)/( diag(z)*V(K,K) ); % This works in noiseless conditions
    [W, WtW, WtV, infos] = nnls_solver(V(:,spa_sol.K)', (diag(y) * V(spa_sol.K, spa_sol.K))', options);
    fprintf('\n');

    if options.verbose > 0
        fprintf('# %s: finished.\n', options.alg_name);
    end

    %% store results
    x.S = S;
    x.W = W';

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function x = mono_next_grlex ( m, x )

%*****************************************************************************80
%
%% MONO_NEXT_GRLEX: grlex next monomial.
%
%  Discussion:
%
%    Example:
%
%    M = 3
%
%    #  X(1)  X(2)  X(3)  Degree
%      +------------------------
%    1 |  0     0     0        0
%      |
%    2 |  0     0     1        1
%    3 |  0     1     0        1
%    4 |  1     0     0        1
%      |
%    5 |  0     0     2        2
%    6 |  0     1     1        2
%    7 |  0     2     0        2
%    8 |  1     0     1        2
%    9 |  1     1     0        2
%   10 |  2     0     0        2
%      |
%   11 |  0     0     3        3
%   12 |  0     1     2        3
%   13 |  0     2     1        3
%   14 |  0     3     0        3
%   15 |  1     0     2        3
%   16 |  1     1     1        3
%   17 |  1     2     0        3
%   18 |  2     0     1        3
%   19 |  2     1     0        3
%   20 |  3     0     0        3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer X(M), the current monomial.
%    The first item is X = [ 0, 0, ..., 0, 0 ].
%
%    Output, integer X(M), the next monomial.
%

%
%  Ensure that 1 <= M.
%
  if ( m < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MONO_NEXT_GRLEX - Fatal error!' );
    fprintf ( 1, '  M < 1\n' );
    error ( 'MONO_NEXT_GRLEX - Fatal error!' );
  end
%
%  Ensure that 0 <= XC(I).
%
  for i = 1 : m
    if ( x(i) < 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'MONO_NEXT_GRLEX - Fatal error!' );
      fprintf ( 1, '  X(I) < 0\n' );
      error ( 'MONO_NEXT_GRLEX - Fatal error!' );
    end
  end
%
%  Find I, the index of the rightmost nonzero entry of X.
%
  i = 0;
  for j = m : -1 : 1
    if ( 0 < x(j) )
      i = j;
      break
    end
  end
%
%  set T = X(I)
%  set X(I) to zero,
%  increase X(I-1) by 1,
%  increment X(M) by T-1.
%
  if ( i == 0 )
    x(m) = 1;
    return
  elseif ( i == 1 )
    t = x(1) + 1;
    im1 = m;
  elseif ( 1 < i )
    t = x(i);
    im1 = i - 1;
  end

  x(i) = 0;
  x(im1) = x(im1) + 1;
  x(m) = x(m) + t - 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = eegsynth_tsv(filename, hdr)

% EEGSYNTH_TSV is called from FT_READ_EVENT to read the events from a tsv file
% written by the recordtrigger module. The .tsv file should also contain a
% synchronization trigger from the recordsignal module.
%
% Use as
%   hdr = events_tsv(filename)
%   evt = events_tsv(filename, hdr)
% to read the header or the event information. Note that when reading the header, the
% number of channels in the actual data is unknown.
%
% See https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/05-task-events.html
%
% See also FT_FILETYPE, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, QUALISYS_TSV, EVENTS_TSV, BIDS_TSV

% FIXME it might be that the events are one sample off, but I cannot be bothered
% checking that precisely at the moment of the initial implementation.

% Copyright (C) 2022, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

needhdr = (nargin==1);
needevt = (nargin==2);

% read the table from file
t = readtable(filename, 'filetype', 'text', 'delimiter', '\t');

type = unique(t.event);
sync = type(endsWith(type, 'synchronize'));
assert(length(sync)==1, 'there should be exactly one type of trigger that ends with "synchronize"')
sync = sync{1};

% find the mapping between timestamps and samples
sel = strcmp(t.event, sync);
firstTimeStamp = t.timestamp(1);
sec = seconds(t.timestamp(sel) - firstTimeStamp);
smp = t.value(sel); % the "synchronize" triggers contain the sample number of the corresponding EEG file

% figure
% plot(sec, smp, '.')
% xlabel('seconds')
% ylabel('samples')

if needhdr
  % fit a first-order polynomial that expresses samples per timestamp
  p = polyfit(sec, smp, 1);
  samplesPerTimeStamp = p(1);

  % fit a first-order polynomial that expresses timestamps per sample
  p = polyfit(smp, sec, 1);

  % construct a dummy header
  hdr = [];
  hdr.nChans = 1;
  hdr.label = {'unknown'};
  hdr.nSamples = inf;
  hdr.nSamplesPre = 0;
  hdr.nTrials = 1;
  hdr.Fs = samplesPerTimeStamp;
  hdr.FirstTimeStamp = firstTimeStamp + seconds(polyval(p, 1));
  hdr.TimeStampPerSample  = samplesPerTimeStamp;

  % return the header
  varargout = {hdr};
end % if needhdr

if needevt
  % fit a first-order polynomial that expresses samples per timestamp
  p = polyfit(sec, smp, 1);

  % estimate the sample number from the timestamp for all other events
  sec = seconds(t.timestamp - firstTimeStamp);
  smp = round(polyval(p, sec)) - 1;

  event = struct;
  for i=1:numel(smp)
    event(i).type = t.event{i};
    event(i).value = t.value(i);
    event(i).sample = smp(i);
    event(i).duration = [];
    event(i).offset = [];
  end

  % return the events
  varargout = {event};
end % if needevt


text:
function out = isnan(f)
%ISNAN   Test if a TRIGTECH has any NaN values.
%   ISNAN(F) returns TRUE if F has any NaN values and FALSE otherwise.
%
% See also ISFINITE, ISINF.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Check if any values are NaN:
out = any(isnan(f.values(:)));

end

author: chebfun

repo: chebfun


text:
function varargout = about_dlg(varargin)
% ABOUT_DLG M-file for about_dlg.fig
%      ABOUT_DLG, by itself, creates a new ABOUT_DLG or raises the existing
%      singleton*.
%
%      H = ABOUT_DLG returns the handle to a new ABOUT_DLG or the handle to
%      the existing singleton*.
%
%      ABOUT_DLG('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ABOUT_DLG.M with the given input arguments.
%
%      ABOUT_DLG('Property','Value',...) creates a new ABOUT_DLG or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before about_dlg_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to about_dlg_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 about_dlg

% Last Modified by GUIDE v2.5 14-Nov-2005 01:26:22

%* Author: Ranjith Unnikrishnan                                          *
%* Carnegie Mellon University, Vision and Mobile Robotics Laboratory     *
%* THE MATERIAL EMBODIED IN THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"     *
%* AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,      *
%* INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR      *
%* FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL CARNEGIE MELLON  *
%* UNIVERSITY BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,            *
%* SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY         *
%* KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,        *
%* LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF     *
%* THIRD PARTIES, WHETHER OR NOT CARNEGIE MELLON UNIVERSITY HAS BEEN     *
%* ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON        *
%* ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE     *
%* POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.                      *
%

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @about_dlg_OpeningFcn, ...
                   'gui_OutputFcn',  @about_dlg_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 about_dlg is made visible.
function about_dlg_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 about_dlg (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes about_dlg wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = about_dlg_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;

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

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
function M = shapefitfactory(VJt)
% Linear manifold structure for optimization over the ShapeFit search space
%
% function M = shapefitfactory(VJt)
%
% Input: VJt is a matrix of size dxn, such that VJt * ones(n, 1) = 0.
%
% Returns M, a structure describing the Euclidean space of d-by-n matrices
% equipped with the standard Frobenius distance and associated trace inner
% product, as a manifold for Manopt. Matrices on M, denoted by T, have size
% dxn and obey T*ones(n, 1) = 0 (centered columns) and <VJt, T> = 1, where
% <A, B> = Trace(A' * B).
%
% See this paper: http://arxiv.org/abs/1506.01437
% ShapeFit: Exact location recovery from corrupted pairwise directions, 2015
% Paul Hand, Choongbum Lee, Vladislav Voroninski
%
% See also: shapefit_smoothed

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

    [d, n] = size(VJt);

    M.name = @() sprintf('ShapeFit space of size %d x %d', d, n);

    M.dim = @() d*n - d - 1;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) norm(x-y, 'fro');

    M.typicaldist = @() sqrt(d*n);

    M.proj = @(T, U) projection(U);
    VJt_normed = VJt / norm(VJt, 'fro');
    function PU = projection(U)
        % Center the columns
        PU = bsxfun(@minus, U, mean(U, 2));
        % Remove component along VJt
        % Note: these two actions can be executed separately, without
        % interference, owing to VJt having centered columns itself.
        PU = PU - (VJt_normed(:)'*U(:))*VJt_normed;
    end

    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @(x, eg, eh, d) projection(eh);

    M.tangent = @(x, d) d;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.randvec = @(x) randvec();
    function u = randvec()
        u = projection(randn(d, n));
        u = u / norm(u, 'fro');
    end

    % We exploit the fact that VJt_normed belongs to the manifold
    M.rand = @() VJt_normed + randn(1) * randvec();

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(d, n);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [d, n]);
    M.vecmatareisometries = @() true;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function varargout=wfdbexec(varargin)
%
% [output]=wfdbexec(commandName,inputArguments,logLevel)
% [nativeCommands]=wfdbexec()
%
% Executes a WFDB native command ('commandName'), with input arguments
% given by 'inputArguments', and outputs everything into a cell array
% called 'output'. The execution is done through a system call outside of
% MATLAB. If no input arguments are provided to the function, the
% output is an array list of all available native WFDB commands in this OS.
% Substantial, non-trivial text parsing maybe required for processing the
% output into a useful MATLAB numerical variable.
%
% WFDBEXEC is useful in cases where a specific command, or feature for a
% command, is not yet in implemented in the MATLAB wrapper. This function
% may be useful for those wishing to do some debugging or performance comparisons.
%
% The user should be very carefull when using this command. The user should
% be familiar with the input and output arguments of the native command
% that he/she is using. In some cases, information about the command can be
% obtained by running the command by itself, of with either '-h',
% '-help','--help' as input arguments in order to access the command help
% information. Offcourse, you can also look at the man page or source code
% for the command at http://www.physionet.org/.
%
%
%
%Required Parameters:
%
% commandName
%       String specifying the command to be called. To ge a list of
%       available WFDB commands that may be runnable by system calls run
%       this function by itself:
%               [nativeCommands]=wfdbexec();
%
%inputArguments
%       Cell array of strings specifying the inputArguments to the command.
%       Each command flag is an element on the cell array. Command flags
%       that require additional parameters should be followed by another
%       String element with the requied parameter(s). In cases you are
%       usign a command that is already implemente in a MATLAB wrapper, it
%       maybe helpful to look at that command's MATLAB code.
%
%Optional Parameters:
%
%       logLevel
%       1x1 integer that specifies the logleve (verbosity) of the execution process.
%       Options are:
%                     0 OFF (Default)
%                     1 SEVERE
%                     2 WARNING
%                     3 INFO
%                     4 FINEST
%                     5 ALL
%
%Output:
%
%output
%       Nx1 Cell array list of Strings.
%
%
% Written by Ikaro Silva, 2013
% Last Modified: January 7, 2014
% Version 0.0.1
% Since 0.9.5
%
% %Get a  list of all WFDB native executables
% nativeComands=wfdbexec()
%
% %Get help on RDSAMP
% out=wfdbexec('rdsamp',{'-h'})
%
% %Read the first five samples of a record
% %Note: This will be very inefficient because we are not buffering
% %and we are returning a list of Strings (instead of doubles and ints).
% %This is provided just as an example on how WFDBEXEC works.
% out=wfdbexec('rdsamp',{'-r','mitdb/100','-t','s5'})
%
% See also WFDB

%endOfHelp
logLevel=[];

[~,config]=wfdbloadlib;

if(nargin==0)
    %With no arguments passed in, we provide the user a list of native
    %commands available for this OS.
    del=findstr(' ',config.WFDB_NATIVE_BIN);
    if(~isempty(del) && del(1)==1)
        config.WFDB_NATIVE_BIN(1)=[];
    end
    commandList={};
    if(ispc)
        [~,message,~] = fileattrib([config.WFDB_NATIVE_BIN '/bin/*.exec']);
        commandList=cell(length(message),1);
        for n=1:length(message)
            str=message(n).Name;
            st=strfind(str,filesep);
            if(~isempty(st))
                str=str(st(end)+1:end);
            end
            commandList(n)={str};
        end
    else
        [~,message,~] = fileattrib([config.WFDB_NATIVE_BIN '/bin/*']);
        for n=1:length(message)
            if(~message(n).directory && message(n).UserExecute)
                str=message(n).Name;
                st=strfind(str,filesep);
                if(~isempty(st))
                    str=str(st(end)+1:end);
                end
                if(isempty(strfind(str,'lib')))
                    commandList(end+1)={str};
                end
            end
        end
    end
    varargout{1}=commandList(:);
else
    if(nargin>2)
        logLevel=varargin{3};
    end
    if(~wfdbloadlib)
        %Add classes to dynamic path
        wfdbloadlib;
    end
    javaWfdbExec=org.physionet.wfdb.Wfdbexec(varargin(1),strcmp(config.customArchFlag,'true'));
    if(~isempty(logLevel))
        javaWfdbExec.setLogLevel(logLevel);
    end
    try
        %System call
        varargout{1}=cell(javaWfdbExec.execToStringList(varargin{2}).toArray);
        clear javaWfdbExec; %Clean-up to avoid leaked Java classes
    catch exception
        clear javaWfdbExec %Clean-up to avoid leaked Java classes
        rethrow(exception)
    end

end

author: ikarosilva

repo: wfdb-app-toolbox


text:
function h = bc_test01 ( x_num, x, t, h )

%*****************************************************************************80
%
%% BC_TEST01 evaluates the boundary conditions for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM,1), the node coordinates.
%
%    Input, real T, the current time.
%
%    Input, real H(X_NUM,1), the current heat values.
%
%    Output, real H(X_NUM,1), the current heat values, after boundary
%    conditions have been imposed.
%
  h(1,1)     = 90.0;
  h(x_num,1) = 70.0;

  return
end


text:
view=getSelectedVolume;
classFileName = [view.leftPath(1:end-5) '.class'];
c = readClassFile(classFileName);
anatSize = size(view.anat);

voi = c.header.voi;
anat = uint8(zeros([anatSize, 3]));
wm = repmat(logical(0), anatSize);
wm([voi(3):voi(4)], [voi(1):voi(2)], [voi(5):voi(6)]) = permute(c.data, [2,1,3])==c.type.white;
for(ii=1:3)
    anat(:,:,:,ii) = view.anat;
end
tmp = anat(:,:,:,1);
tmp(wm) = uint8(double(tmp(wm))+20);
anat(:,:,:,1) = tmp;

mHires = makeMontage3(anat(:,:,:,1), anat(:,:,:,2), anat(:,:,:,3), [voi(5):voi(6)], view.mmPerVox(1), 0);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = file_char_count ( file_name )

%*****************************************************************************80
%
%% FILE_CHAR_COUNT counts the number of characters in a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the input file.
%
%    Output, integer VALUE, the number of characters in the file.
%
  value = -1;

  input = fopen ( file_name );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_CHAR_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', file_name );
    error ( 'FILE_CHAR_COUNT - Error!' );
    return;
  end

  value = 0;

  while ( 1 )

    [ c, count ] = fread ( input, 1 );

    if ( count < 1 )
      break;
    end

    value = value + 1;

  end

  fclose ( input );

  return
end