text:
function ticklabelformat(hAxes,axName,format)
%TICKLABELFORMAT  Sets axis tick labels format
%
% TICKLABELFORMAT enables setting the format of the tick labels
%
% Syntax:
%    ticklabelformat(hAxes,axName,format)
%
% Input Parameters:
%    hAxes  - handle to the modified axes, such as returned by the gca function
%    axName - name(s) of axles to modify: 'x','y','z' or combination (e.g. 'xy')
%    format - format of the tick labels in sprintf format (e.g. '%.1f V') or a
%             function handle that will be called whenever labels need to be updated
%
%    Note: Calling TICKLABELFORMAT again with an empty ([] or '') format will revert
%    ^^^^  to Matlab's normal tick labels display behavior
%
% Examples:
%    ticklabelformat(gca,'y','%.6g V') - sets y axis on current axes to display 6 significant digits
%    ticklabelformat(gca,'xy','%.2f')  - sets x & y axes on current axes to display 2 decimal digits
%    ticklabelformat(gca,'z',@myCbFcn) - sets a function to update the Z tick labels on current axes
%    ticklabelformat(gca,'z',{@myCbFcn,extraData}) - sets an update function as above, with extra data
%
% Warning:
%    This code heavily relies on undocumented and unsupported Matlab functionality.
%    It works on Matlab 7+, but use at your own risk!
%
% Technical description and more details:
%    http://UndocumentedMatlab.com/blog/setting-axes-tick-labels-format/
%
% Bugs and suggestions:
%    Please send to Yair Altman (altmany at gmail dot com)
%
% Change log:
%    2012-04-18: first version
%
% See also: sprintf, gca

    % Check # of args (we now have narginchk but this is not available on older Matlab releases)
    if nargin < 1
        help(mfilename)
        return;
    elseif nargin < 3
        error('YMA:TICKLABELFORMAT:ARGS','Not enough input arguments');
    end

    % Check input args
    if ~ishandle(hAxes) || ~isa(handle(hAxes),'axes')
        error('YMA:TICKLABELFORMAT:hAxes','hAxes input argument must be a valid axes handle');
    elseif ~ischar(axName)
        error('YMA:TICKLABELFORMAT:axName','axName input argument must be a string');
    elseif ~isempty(format) && ~ischar(format) && ~isa(format,'function_handle') && ~iscell(format)
        error('YMA:TICKLABELFORMAT:format','format input argument must be a string or function handle');
    end

    % normalize axes name(s) to lowercase
    axName = lower(axName);

    if strfind(axName,'x')
        install_adjust_ticklbl(hAxes,'X',format)
    end
    if strfind(axName,'y')
        install_adjust_ticklbl(hAxes,'Y',format)
    end
    if strfind(axName,'z')
        install_adjust_ticklbl(hAxes,'Z',format)
    end

% Install the new tick labels for the specified axes
function install_adjust_ticklbl(hAxes,axName,format)

    % If empty format was specified
    if isempty(format)
        % Remove the current format (revert to default Matlab format)
        set(hAxes,[axName 'TickLabelMode'],'auto')
        setappdata(hAxes,[axName 'TickListener'],[])
        return
    end

    % Determine whether to use the specified format as a
    % sprintf format or a user-specified callback function
    if ischar(format)
        cb = {@adjust_ticklbl axName format};
    else
        cb = format;
    end

    % Now install axis tick listeners to adjust tick labels
    % (use undocumented feature for adjustments)
    ha = handle(hAxes);
    hp = findprop(ha,[axName 'Tick']);
    hl = handle.listener(ha,hp,'PropertyPostSet',cb);
    setappdata(hAxes,[axName 'TickListener'],hl)

    % Adjust tick labels now
    %eventData.AffectedObject = hAxes;
    %adjust_ticklbl([],eventData,axName,format)
    set(hAxes,[axName 'TickLabelMode'],'manual')
    set(hAxes,[axName 'TickLabelMode'],'auto')

% Default tick labels update callback function (used if user did not specify their own function)
function adjust_ticklbl(hProp,eventData,axName,format)	%#ok<INUSL>
    hAxes = eventData.AffectedObject;
    tickValues = get(hAxes,[axName 'Tick']);
    tickLabels = arrayfun(@(x)(sprintf(format,x)),tickValues,'UniformOutput',false);
    set(hAxes,[axName 'TickLabel'],tickLabels)


text:
function studies = studyLoad;
% Load a studies struct, specifying saved groupings of sessions for
% mrVista 2.
%
%  studies = studyLoad;
%
% STRUCTURE OF STUDIES VARIABLE:
%   studies is a struct array that is saved in the following path relative
%   to the root mrVista2 repository:
%       mrVista2/study/studies.mat
%   If this is not found, it creates and saves it, with a default "(Recent
%   Sessions)" study which keeps track of sessions viewed with the mrVista
%   2 session GUI.
%
%   Each entry in studies has the following fields
%       name: brief name of study.
%       sessions: set of paths to sessions contained within the study. Can
%           be absolute paths, or relative to a parent directory (see
%           studyDir below). A cell array of strings.
%       studyDir: provides a parent directory for the study. If non-empty,
%           this could be a string that is evaluated. That is, it could be
%           a regular string, e.g., 'X:\myStudies\study1\', or a MATLAB
%           command that will evaluate to the path on the local machine,
%           e.g. 'fullfile(RAID, 'myStudies', 'study1');'. This allows
%           sessions to be accessed from multiple machines. If empty,
%           will assume the sessions field contains absolute paths.
%           The studyDir field may also be useful for saving the results of
%           study-level analyses, across sessions.
%       params: a struct of parameters relevant to analyses specific to
%       this study. The details can vary depending on the study.
%       comments: string containing any comments about the study.
%
%
% ras 07/03/06.
parentDir = fileparts(which(mfilename)); % parent directory of this M-file.
studyPath = fullfile(parentDir, 'studies.mat');

if ~exist(studyPath, 'file')
    fprintf('Creating file %s which will contain saved study information...', ...
            studyPath);

    studies(1) = studyCreate('(Recent Sessions)');
    studies(1).sessions = {pwd};

    save(studyPath, 'studies');
    disp('Done.')
else
    load(studyPath, 'studies');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function btf_make
%BTF_MAKE compile BTF for use in MATLAB
% Your current working directory must be BTF/MATLAB for this function to work.
%
% Example:
%   btf_make
%
% See also btf, maxtrans, stroncomp, dmperm.

% Copyright 2004-2007, Tim Davis, University of Florida

details = 0 ;       % if 1, print details of each command

mexcmd = 'mex -O -DDLONG -I../Include -I../../UFconfig ' ;
if (~isempty (strfind (computer, '64')))
    mexcmd = [mexcmd '-largeArrayDims '] ;
end

s = [mexcmd 'maxtrans.c ../Source/btf_maxtrans.c'] ;
if (details)
    fprintf ('%s\n', s) ;
end
eval (s) ;

s = [mexcmd 'strongcomp.c ../Source/btf_strongcomp.c'] ;
if (details)
    fprintf ('%s\n', s) ;
end
eval (s) ;

s = [mexcmd 'btf.c ../Source/btf_maxtrans.c ../Source/btf_strongcomp.c ../Source/btf_order.c'] ;
if (details)
    fprintf ('%s\n', s) ;
end
eval (s) ;

fprintf ('BTF successfully compiled.\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function latLonPts=polarAzEquidistProj2Ellipse(pAzPts,latLonRef,a,f)
%%POLARAZEQUIDISTANTPROJ2ELLIPSE Given points as polar coordinates of an
%       azimuthal equidistant projection, convert the points to latitudes
%       and longitudes on a reference ellipsoid (or sphere). A third height
%       coordinate can also be provided, which doesn't change in the
%       conversion.
%
%INPUTS: pAzPts A 2XN set of the [ground distance; heading] points, with the
%              heading given in radians East of North, to convert.
%              Alternatively, if heights are also given, this can be a 3XN
%              set of points with the height being the third dimension.
%    latLonRef A 2X1 [latitude;longitude] reference point in radians about
%              which the projection is taken.
%            a The semi-major axis of the reference ellipsoid (in meters).
%              If this argument is omitted or an empty matrix is passed,
%              the value in Constants.WGS84SemiMajorAxis is used.
%            f The flattening factor of the reference ellipsoid. If this
%              argument is omitted or an empty matrix is passed, the value
%              in Constants.WGS84Flattening is used.
%
%OUTPUTS: latLonPts The 2XN set of converted [latitude;longitude] points in
%                   radians. If heights were given in xyPts, then this is a
%                   3XN set of converted [latitude;longitude;height] points
%                   with the third row the same as in xyPts.
%
%The conversion is described in Chapter 25 of [1], where expressions for
%a spherical Earth are given. However, we do not use those formulae. The
%latitude and longitude of the point can be obtained using the
%directGreatCircleProb function for a spherical Earth or the
%directGeodeticProb for an ellipsoidal Earth, which is just what is done
%here.
%
%December 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(f))
    f=Constants.WGS84Flattening;
end

if(nargin<3||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

hasHeight=(size(pAzPts,1)==3);
numPts=size(pAzPts,2);
if(hasHeight)
    latLonPts=zeros(3,numPts);
else
    latLonPts=zeros(2,numPts);
end
if(f==0)
    %Under a spherical Earth approximation.
    for k=1:numPts
        distVal=pAzPts(1,k);
        az=pAzPts(2,k);

        latLonPts(1:2,k)=directGreatCircleProb(latLonRef,az,distVal,a);
    end
else
    %Under an ellipsoidal Earth approximation.
    for k=1:numPts
        distVal=pAzPts(1,k);
        az=pAzPts(2,k);

        latLonPts(1:2,k)=directGeodeticProb(latLonRef,az,distVal,a,f);
    end
end

if(hasHeight)
    latLonPts(3,:)=pAzPts(3,:);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function update_from_affine_matrix(this, ...
    affineMatrix)
% Updates properties of MrImageGeometry from affine 4x4 transformation
% matrix
%
%   Y = MrImageGeometry()
%   Y.update_from_affine_matrix(affineMatrix)
%
% This is a method of class MrImageGeometry.
%
% IN
%
% OUT
%
% EXAMPLE
%   update_from_affine_matrix
%
%   See also MrImageGeometry tapas_uniqc_spm_matrix, tapas_uniqc_spm_imatrix

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-07-27
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% round to N decimals for specified precision, to avoid small numbers < double precision
N = floor(abs(log10(eps('double'))));
P = round(tapas_uniqc_spm_imatrix(affineMatrix),N);

this.offcenter_mm       = P(1:3);
this.rotation_deg       = P(4:6)/pi*180;
this.resolution_mm      = P(7:9);
this.shear              = P(10:12);
this.FOV_mm             = this.resolution_mm.*...
    this.nVoxels(1:3);

author: translationalneuromodeling

repo: tapas


text:
function distance = word_distance(word1, word2, model)
% Shows the L2 distance between word1 and word2 in the word_embedding_weights.
% Inputs:
%   word1: The first word as a string.
%   word2: The second word as a string.
%   model: Model returned by the training script.
% Example usage:
%   word_distance('school', 'university', model);

word_embedding_weights = model.word_embedding_weights;
vocab = model.vocab;
id1 = strmatch(word1, vocab, 'exact');
id2 = strmatch(word2, vocab, 'exact');
if ~any(id1)
  fprintf(1, 'Word ''%s\'' not in vocabulary.\n', word1);
  return;
end
if ~any(id2)
  fprintf(1, 'Word ''%s\'' not in vocabulary.\n', word2);
  return;
end
word_rep1 = word_embedding_weights(id1, :);
word_rep2 = word_embedding_weights(id2, :);
diff = word_rep1 - word_rep2;
distance = sqrt(sum(diff .* diff));

repo: neural-nets


text:
function obj = tidal_data_to_ob(obj,tidal_database,const)
% obj = tidal_data_to_ob(obj,tidal_database,const)
% Input a msh class obj with open boundary locations and tidal
% constituents, and interpolate the tidal database constituents onto it.
% Put the result into the f15 struct of the msh obj.
%
% Requires: m_map for projection
%
% Created by William Pringle March 15 2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check entry
% if we are using wildcard for each constituent in list of files
cidx = strfind(tidal_database,'**') ;
if ~isempty(cidx)
    tidal_database(cidx:cidx+1) = lower(const{1});
end
if ~exist(tidal_database,'file')
   error(['tidal database file does not exist: ' tidal_database])
end
if isempty(obj.op)
   error('No open boundaries in the mesh to put BCs on')
end

% Select desired projection (using m_map)
proj = 'Mercator';

% Specify limits of grid
lat_min = min(obj.p(:,2)) - 1; lat_max = max(obj.p(:,2)) + 1;
lon_min = min(obj.p(:,1)) - 1; lon_max = max(obj.p(:,1)) + 1;

% doing the projection
m_proj(proj,'lon',[ lon_min lon_max],...
            'lat',[ lat_min lat_max])

%% Get boundary info
opedat = obj.op;
b_lon = zeros(opedat.neta,1);
b_lat = zeros(opedat.neta,1);
node_num = zeros(opedat.neta,1);
ns = 1; ne = 0;
for n = 1:opedat.nope
    ne = ne + opedat.nvdll(n);
    node_num(ns:ne) = opedat.nbdv(1:opedat.nvdll(n),n);
    b_lon(ns:ne)    = obj.p(node_num(ns:ne),1);
    b_lat(ns:ne)    = obj.p(node_num(ns:ne),2);
    ns = ne + 1;
end

% Do projection
[b_x,b_y] = m_ll2xy(b_lon,b_lat);

%% Load tide data and make vectors
lon = ncread(tidal_database,'lon_z');
lat = ncread(tidal_database,'lat_z');
const_t = ncread(tidal_database,'con');
if min(size(lon)) == 1
    [lon,lat] = meshgrid(lon,lat);
end
if max(obj.p(:,1)) <= 180
    % change to -180 to 180 format if msh is in that format
    lon(lon > 180) = lon(lon > 180) - 360;
end
lon_x = reshape(lon,[],1);
lat_y = reshape(lat,[],1);

% Delete uncessecary portions
% First delete by square
I = find(lon_x < lon_min | lon_x > lon_max | ...
         lat_y < lat_min | lat_y > lat_max);
lon_x(I) = []; lat_y(I) = [];

% Get 20 nearest neighbours
Kd = ourKNNsearch([lon_x,lat_y]',[b_lon, b_lat]',20);
Kd = unique(Kd);

% The new lon and lat vectors of data
lon_x = lon_x(Kd); lat_y = lat_y(Kd);

% Do the projection
[x,y] = m_ll2xy(lon_x,lat_y);

%% Now interpolate into f15 struct
keep = true(obj.f15.nbfr,1);
for j = 1:obj.f15.nbfr
    % Read the current consituent
    if ~isempty(cidx)
        tidal_database(cidx:cidx+1) = lower(const{j});
        if exist(tidal_database,'file')
            const_t = ncread(tidal_database,'con');
        end
    end
    k = find(startsWith(string(const_t'),lower(const{j})));
    if isempty(k)
        disp(['No tidal data in file for constituent ' const{j}])
        obj.f15.opealpha(j).name = const{j};
        obj.f15.opealpha(j).val = 0*[b_x b_x];
        continue
    end
        % For real part
    if min(size(const_t)) > 1
        Re_now = ncread(tidal_database,'hRe',[1 1 k],[size(lon) 1]);
    else
        Re_now = double(ncread(tidal_database,'hRe'))/1000;
    end
    % reshape to vector
    Re_now = reshape(Re_now,[],1);
    % For imaginary part
    if min(size(const_t)) > 1
        Im_now = ncread(tidal_database,'hIm',[1 1 k],[size(lon) 1]);
    else
        Im_now = double(ncread(tidal_database,'hIm'))/1000;
    end
    % reshape to vector
    Im_now = reshape(Im_now,[],1);
    % Eliminate regions outside of search area and on land
    % Linear extrapolation of ocean values will be conducted where
    % boundary nodes fall inside a land cell of the tidal data.
    Re_now(I) = []; Re_now = Re_now(Kd);
    K = find(Re_now == 0); Re_now(K) = [];
    Im_now(I) = []; Im_now = Im_now(Kd); Im_now(K) = [];
    xx = x; yy = y; xx(K) = []; yy(K) = [];
    % Make into complex number
    Z = Re_now - Im_now*1i;
    % Do the scattered Interpolation
    F = scatteredInterpolant(xx,yy,Z,'natural');
    BZ = F(b_x,b_y);
        %
    % Convert real and imaginary parts to amplitude and phase
    amp_b = abs(BZ);
    phs_b = rad2deg(angle(BZ));
    % Convert to 0 to 360;
    % phs_b that is positive 0 - 180 stays same.
    % phs_b that is negative -180 - 0 becomes 180 - 360
    phs_b(phs_b < 0) = phs_b(phs_b < 0) + 360;

    obj.f15.opealpha(j).name = const{j};
    obj.f15.opealpha(j).val = [amp_b phs_b];

end
obj.f15.nbfr = length(find(keep));
obj.f15.opealpha = obj.f15.opealpha(keep);
obj.f15.bountag = obj.f15.bountag(keep);
%EOF
end

author: CHLNDDEV


text:
function r = r8pbu_resid ( m, n, mu, a, x, b )

%*****************************************************************************80
%
%% R8PBU_RESID computes the residual R = B-A*X for R8PBU matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer MU, the number of superdiagonals in the matrix.
%    MU must be at least 0 and no more than N-1.
%
%    Input, real A(MU+1,N), the matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Input, real B(M), the desired result A * x.
%
%    Output, real R(M), the residual R = B - A * X.
%
  r = r8pbu_mv ( m, n, mu, a, x );

  r(1:m) = b(1:m) - r(1:m);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function region_bb = get_region_bbox(regions, sp_bb)

for obj = 1:numel(regions)
   sub_bb = sp_bb(regions{obj},:);
   if(numel(sub_bb) > 0)
      region_bb(obj, :) = [min(sub_bb(:, 1)), min(sub_bb(:, 2)), max(sub_bb(:,3)), max(sub_bb(:, 4))];
   else
      region_bb(obj, :) = [-1 -1 -1 -1];
   end
end

author: Cloud-CV

repo: object-proposals


text:
function varargout = randsample(varargin)
% - Produces a column, like Matlab does (instead of Octave, which produces a row), when the first input is a scalar
% - Allows first input to be char. It is then treated as the population to sample from (not as a number)
if ischar(varargin{1})
    varargout{1} = varargin{1}(randsample(numel(varargin{1}), varargin{2:end}));
else
    varargout{1} = builtin('randsample', varargin{:});
    if numel(varargin{1})==1
        varargout{1} = varargout{1}(:);
    end
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function lung_roi = PTKGetLeftLungROIFromLeftAndRightLungs(lung_image, left_and_right_lung_mask, reporting)
    % PTKGetLeftLungROIFromLeftAndRightLungs. Extracts a region of interest for the left lung given the original image and the left and right lung mask.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    mask = left_and_right_lung_mask.Copy;
    mask.ChangeRawImage(uint8(mask.RawImage == 2));
    mask.CropToFitWithBorder(5);
    mask.AddBorder(5);
    lung_roi = lung_image.Copy;
    lung_roi.ResizeToMatch(mask);
end

author: tomdoel

repo: pulmonarytoolkit


text:
%%*************************************************************************
%% HSDbicgstab
%%
%% [xx,resnrm,flag] = HSDbicgstab(A,b,M1,tol,maxit)
%%
%% iterate on  bb - (M1)*AA*x
%%
%% r = b-A*xtrue;
%%
%%*************************************************************************

function [xx,resnrm,flag] = HSDbicgstab(A,b,M1,tol,maxit,printlevel)

N = length(b);
if (nargin < 6); printlevel = 1; end
if (nargin < 5) || isempty(maxit); maxit = max(20,length(A.mat22)); end;
if (nargin < 4) || isempty(tol); tol = 1e-8; end;
tolb = min(1e-4,tol*norm(b));
flag = 1;

x = zeros(N,1);
if isstruct(A); r = b-matvec(A,x); else r = b-mexMatvec(A,x); end;
err = norm(r); resnrm(1) = err;  minresnrm = err; xx = x;
%%if (err < tolb); return; end

omega = 1.0;
r_tld = r;
%%
%%
%%
smtol = 1e-40;
for iter = 1:maxit,

    rho   = (r_tld'*r);
    if (abs(rho) < smtol)
        flag = 2;
        if (printlevel); fprintf('*'); end;
        break;
    end
    if (iter > 1)
        beta  = (rho/rho_1)* (alp/omega);
        p = r + beta*(p - omega*v);
    else
        p = r;
    end
    p_hat = precond(A,M1,p);
    if isstruct(A); v = matvec(A,p_hat); else v = mexMatvec(A,p_hat); end;
    alp = rho / (r_tld'*v);
    s = r - alp*v;
    %%
    s_hat = precond(A,M1,s);
    if isstruct(A); t = matvec(A,s_hat); else t = mexMatvec(A,s_hat); end;
    omega = (t'*s) / (t'*t);
    x = x + alp*p_hat + omega*s_hat;
    r = s - omega*t;
    rho_1 = rho;
    %%
    %% check convergence
    %%
    err = norm(r); resnrm(iter+1) = err; %#ok
    if (err < minresnrm);
        xx = x; minresnrm = err;
    end
    if (err < tolb)
        break;
    end
    if (err > 10*minresnrm)
        if (printlevel); fprintf('^'); end
        break;
    end
    if (abs(omega) < smtol)
        flag = 2;
        if (printlevel); fprintf('*'); end;
        break;
    end
end
%%
%%*************************************************************************
%%*************************************************************************
%% matvec: matrix-vector multiply.
%% matrix = [A.mat11, A.mat12; A.mat12', A.mat22]
%%*************************************************************************

function Ax = matvec(A,x)

m = length(A.mat11); m2 = length(x)-m;
if (m2 > 0)
    x1 = full(x(1:m));
else
    x1 = full(x);
end
Ax = mexMatvec(A.mat11,x1);
if (m2 > 0)
    x2 = full(x(m+1:m+m2));
    Ax = Ax + mexMatvec(A.mat12,x2);
    Ax2 = mexMatvec(A.mat12,x1,1) + mexMatvec(A.mat22,x2);
    Ax = [Ax; Ax2];
end
%%*************************************************************************
%% precond:
%%*************************************************************************

function Mx = precond(A,L,x)

m = L.matdim; m2 = length(x)-m;
if (m2 > 0)
    x1 = full(x(1:m));
else
    x1 = full(x);
end
if (m2 > 0)
    x2 = x(m+1:m+m2);
    w = linsysolvefun(L,x1);
    z = mexMatvec(A.mat12,w,1) -x2;
    z = L.Mu \ (L.Ml \ (L.Mp*z));
    x1 = x1 - mexMatvec(A.mat12,z);
end
%%
Mx = linsysolvefun(L,x1);
%%
if (m2 > 0)
    Mx = [Mx; z];
end
%%*************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [p]=rmEvaluate2Fits(rss1,rss2,df1,df2,n,units,wtest)
% rmEvaluate2Fits - probability that fit1 == fit2
%
% output = rmEvaluate2Fits(rss1,rss2,df1,df2,n,units,wtest);
%
% Performs f-test and general likelihood ratio (Wilks likelihood)
% tests. Outputs the probability of Ho: data1 == data2.
% Assumes df1>=df2.
%
% rss  : residual sum of squares (model 1 and 2)
% df   : degrees of freedom      (model 1 and 2)
% n    : number of data points
% units: 'p', 'log10p' ['p']
% wtest: 'ftest','glr'
%
% 2006/02 SOD: wrote it.

if ieNotDefined('rss1'), error('Need rss1.'); end;
if ieNotDefined('rss2'), error('Need rss2.'); end;
if ieNotDefined('df1'),  error('Need df1.');  end;
if ieNotDefined('df2'),  error('Need df2.');  end;
if ieNotDefined('n'),    error('Need n.');    end;
if ieNotDefined('units'),wtest = 'p';         end;
if ieNotDefined('wtest'),wtest = 'glr';       end;

% some sanity checks
if df1 == df2,
  disp(sprintf('[%s]:error: df1 == df2.',mfilename));
  return;
end;
if df1<df2,
  disp(sprintf('[%s]:warning: df1>df2 flipping datasets!',mfilename));
  dftmp = df1;
  df1   = df2;
  df2   = dftmp;
  rsstmp = rss1;
  rss1   = rss2;
  rss2   = rsstmp;
end;

% now do tests
switch lower(wtest),
 case {'f','ftest'}
  % ftest
  if df1~=df2,
    f = ((rss1-rss2)./rss2) ./ ((df1-df2)./df2);
    p = 1-fcdf(f,df1-df2,df2);
  else
    f = rss1./rss2;
    p = 1-fcdf(f,df2,df1);
  end;

 case {'g','glr','wilks'}
  % From Fan et al (2001) Generalized likelihood ratio statistics and Wilks
  % phenomenon. The Annals of Statistics. 29(1): 153-193.
  % calculate lambda (glr computation)
  lambda = (n./2).*log(rss1./rss2);

  % probability function: Wilks phenomenon lambda approximates
  % chi-squared distribution.
  df = df1-df2;
  p = 1-chi2cdf(lambda,df);
end;

% convert units if requested
switch lower(units),
 case 'log10p',
  sp         = sign(p);
  sp(sp==0)  = 1;
  p          = abs(p);
  p(p<1e-50) = 10^-50; % remove 0
  p          = -log10(p) .* sp;
 case 'p',
  % do nothing
 otherwise,
  error('[%s]:Unknown unit: %s',mfilename,units);
end;

return;


text:
function sftpack_test06 ( )

%*****************************************************************************80
%
%% TEST06 tests R8VEC_SST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
  n = 256;
  alo = 0.0;
  ahi = 5.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  For slow sine transforms,\n' );
  fprintf ( 1, '  R8VEC_SST does a forward or backward transform.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data items is N = %d\n', n );
%
%  Set the data values.
%
  seed = 123456789;

  [ c, seed ] = r8vec_uniform ( n, alo, ahi, seed );

  r8vec_print_part ( n, c, 10, '  The original data:' );
%
%  Compute the coefficients.
%
  d = r8vec_sst ( n, c );

  r8vec_print_part ( n, d, 10, '  The sine coefficients:' );
%
%  Now compute inverse transform of coefficients.  Should get back the
%  original data.

  e = r8vec_sst ( n, d );

  e(1:n) = e(1:n) / 2 / ( n + 1 );

  r8vec_print_part ( n, e, 10, '  The retrieved data:' );

  return
end


text:
function [I1,I2] = hmxSubdivide(X,dim)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxSubdivide.m                                |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Subdivide particles with median repartition   |
%|  `---'  |                                                              |
%+========================================================================+

% Sort data on largest dimension
[~,I] = sort(X(:,dim));

% Equal repartition
I1 = I(1:floor(end/2));
I2 = I(floor(end/2)+1:end);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [FrEn,avLL] = evalfreeenergy_ehmm(T,Gamma,Xi,ehmm,residuals,XX,todo)
% Computes the Free Energy of an HMM depending on observation model
%
% INPUT
% T            length of series
% Gamma        probability of states conditioned on data
% Xi           joint probability of past and future states conditioned on data
% ehmm          bsmm structure
% residuals    in case we train on residuals, the value of those.
%
% OUTPUT
% FrEn         the variational free energy, separated in different terms:
%                   element 1: Gamma Entropy
%                   element 2: data negative loglikelihood
%                   element 3: Gamma negative loglikelihood
%                   element 4: KL for initial and transition probabilities
%                   elements 5-: KL for the state parameters
% avLL         log likelihood of the observed data, per trial
%
% Author: Diego Vidaurre, OHBA, University of Oxford

if nargin<7, todo = ones(1,5); end

K = ehmm.K;

if nargin>=6, Tres = size(residuals,1);
else, Tres = sum(T) - length(T)*ehmm.train.maxorder;
end
setstateoptions;
ltpi = sum(regressed)/2 * log(2*pi);
np = size(XX,2); ndim = size(residuals,2);

% Entropy of Gamma
Entr = [];
if todo(1)==1
    Entr = GammaEntropy_ehmm(Gamma,Xi,T);
end

% loglikelihood of Gamma
avLLGamma = [];
if todo(3)==1
    avLLGamma = GammaavLL_ehmm(ehmm,Xi);
end

% P and Pi KL
KLdivTran = [];
if todo(4)==1
    KLdivTran = KLtransition_ehmm(ehmm);
end

% state KL
KLdiv = [];
if todo(5)==1

    OmegaKL = 0;
    for n = 1:ndim
        if ~regressed(n), continue; end
        OmegaKL = OmegaKL + gamma_kl(ehmm.Omega.Gam_shape,ehmm.prior.Omega.Gam_shape, ...
            ehmm.Omega.Gam_rate(n),ehmm.prior.Omega.Gam_rate(n));
    end
    KLdiv = [KLdiv OmegaKL];

    WKL = 0;
    Sind_all = [repmat(Sind(:,n),K,1) == 1; false(np,1)];
    for n = 1:ndim
        prior_prec = [];
        for k = 1:K
            hs = ehmm.state(k);
            pr = ehmm.state(k).prior;
            if ~regressed(n), continue; end
            if ~train.zeromean, prior_prec = hs.prior.Mean.iS(n); end
            if train.order > 0
                ndim_n = sum(S(:,n));
                alphaterm = repmat( (hs.alpha.Gam_shape ./  hs.alpha.Gam_rate), ndim_n, 1);
                if ndim==1
                    prior_prec = [prior_prec; alphaterm(:)] ;
                else
                    sigmaterm = repmat(hs.sigma.Gam_shape(S(:,n),n) ./ ...
                        hs.sigma.Gam_rate(S(:,n),n), length(orders), 1);
                    prior_prec = [prior_prec; sigmaterm .* alphaterm(:)] ;
                end
            end
        end
        prior_var = diag(1 ./ prior_prec);
        WKL = WKL + gauss_kl(ehmm.state_shared(n).Mu_W(Sind_all),zeros(np*K,1), ...
            ehmm.state_shared(n).S_W(Sind_all,Sind_all), prior_var);
    end
    KLdiv = [KLdiv WKL];

    sigmaKL = 0; alphaKL = 0;
    for k = 1:K+1
        if isempty(train.prior) && ~isempty(orders) && ~train.uniqueAR && ndim>1
            for n1 = 1:ndim
                for n2 = 1:ndim
                    if (train.symmetricprior && n2<n1) || S(n1,n2)==0, continue; end
                    sigmaKL = sigmaKL + gamma_kl(hs.sigma.Gam_shape(n1,n2),pr.sigma.Gam_shape(n1,n2), ...
                        hs.sigma.Gam_rate(n1,n2),pr.sigma.Gam_rate(n1,n2));
                end
            end
        end
        if isempty(train.prior) &&  ~isempty(orders)
            for i = 1:length(orders)
                alphaKL = alphaKL + gamma_kl(hs.alpha.Gam_shape,pr.alpha.Gam_shape, ...
                    hs.alpha.Gam_rate(i),pr.alpha.Gam_rate(i));
            end
        end
    end
    KLdiv = [KLdiv sigmaKL alphaKL];

end

% data log likelihood
savLL = [];
if todo(2)==1

    Gamma = [Gamma prod(1-Gamma,2) ];
    Gamma = rdiv(Gamma,sum(Gamma,2));

    avLL = zeros(Tres,1);

    ldetWishB = 0;
    PsiWish_alphasum = 0;
    for n = 1:ndim
        if ~regressed(n), continue; end
        ldetWishB = ldetWishB+0.5*log(ehmm.Omega.Gam_rate(n));
        PsiWish_alphasum = PsiWish_alphasum+0.5*psi(ehmm.Omega.Gam_shape);
    end
    C = ehmm.Omega.Gam_shape ./ ehmm.Omega.Gam_rate;
    avLL = avLL + (-ltpi-ldetWishB+PsiWish_alphasum);

    for k = 1:K+1

        % error
        meand = XX * ehmm.state(n).W.Mu_W(:,regressed);
        d = residuals(:,regressed) - meand;
        Cd = bsxfun(@times, C(regressed), d)';
        dist = zeros(Tres,1);
        for n = 1:sum(regressed)
            dist = dist - 0.5 * (d(:,n).*Cd(n,:)');
        end

        % Covariance of the distance
        NormWishtrace = zeros(Tres,1);
%         if ndim == 1
%             NormWishtrace = NormWishtrace + 0.5 * C * ...
%                 sum( (XX * ehmm.state(k).W.S_W) .* XX, 2);
%         else
%             for n = 1:ndim
%                 if ~regressed(n), continue; end
%                 NormWishtrace = NormWishtrace + 0.5 * C(n) * ...
%                     sum( (XX(:,Sind(:,n)) * ...
%                     permute(ehmm.state(k).W.S_W(n,Sind(:,n),Sind(:,n)),[2 3 1])) ...
%                     .* XX(:,Sind(:,n)), 2);
%             end
%         end

        avLL = avLL + Gamma(:,k) .* (dist - NormWishtrace);

    end
    savLL = sum(avLL);
end

FrEn = [-Entr -savLL -avLLGamma +KLdivTran +KLdiv];

% fprintf(' %.10g + %.10g + %.10g + %.10g + %.10g + %.10g =  %.10g \n',...
%     sum(-Entr) ,sum(-(-ltpi-ldetWishB+PsiWish_alphasum)*Tres) , sum(-(dist - NormWishtrace)), ...
%     sum(-avLLGamma) , sum(+KLdivTran) ,sum(+KLdiv) ,sum(FrEn));

end

author: OHBA-analysis

repo: HMM-MAR


text:
function TF = isempty(TC)

%ISEMPTY Display threecomp object
%   TF = ISEMPTY(TC) returns an array the same size as TC, containing
%   logical 1 (true) where the elements of TC are empty, and logical 0
%   (false) elsewhere. Empty is defined as containging no waveform data.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% CHECK INPUTS
if ~isa(TC,'threecomp')
    error('First argment must be a threecomp object');
end
if nargin~=1
    error('Incorrect number of arguments');
end

% CHECK EACH ELEMENT TO SEE IF IT IS EMPTY
[objSize1,objSize2] = size(TC);
numObj = numel(TC);
TF = zeros(size(TC));
for n = 1:numObj
    numSamples = min(get(TC(n).traces,'DATA_LENGTH'));
    if numSamples == 0
        TF(n) = 1;
    else
        TF(n) = 0;
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function [ y, next ] = p08_equil ( neqn, next )

%*****************************************************************************80
%
%% P08_EQUIL returns equilibrium solutions of problem p08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%

  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = [ 0.0; 0.0; 0.0 ];
  else
    next = 0;
    y = [];
  end

  return
end


text:
function element_size = p05_element_size ( )

%*****************************************************************************80
%
%% P05_ELEMENT_SIZE returns a typical element size for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real ELEMENT_SIZE, a typical element size.
%
  element_size = 0.015;

  return
end


text:
function f = logsigmoid_fh(w)
% This is a factory for a function handle to an MV2DF, which represents
% the vectorization of the logsigmoid function. The mapping is, in
% MATLAB-style code:
%
%   y = log(sigmoid(w)) = log(1./1+exp(-w)) = -log(1+exp(-w))
%
% Inputs:
%   m: the number of inputs to each individual logsumexp calculation.
%   direction: 1 sums down columns, or 2 sums accross rows.
%   w: optional, if ssupplied
%
% Outputs:
%   f: a function handle to the MV2DF described above.
%
% see: MV2DF_API_DEFINITION.readme

if nargin==0
    test_this();
    return;
end

f = vectorized_function([],@(x)F0(x));

if exist('w','var') && ~isempty(w)
    f = f(w);
end

end

function [y,f1] = F0(x)
logp1 = -neglogsigmoid(x);
logp2 = -neglogsigmoid(-x);
y = logp1;
f1 = @() F1(logp1,logp2);
end

function [J,f2,linear] = F1(logp1,logp2)
linear = false;
J = exp(logp2);
f2 = @(dx) F2(dx,logp1,logp2);
end

function h = F2(dx,logp1,logp2)
h = -dx.*exp(logp1+logp2);
end

function test_this()
n = 10;
f = logsigmoid_fh([]);
x = randn(n,1);
test_MV2DF(f,x);
end

author: nesl


text:
function [trans,rot,scale,skew] = affineDecompose(A)
%
% [trans,rot,scale,skew] = affineDecompose(A)
%
% Decomposes an affine (orthogonal linear) transformation matrix into its
% four componets: translations, rotations scales and skews.
%
% (See affineBuild for more info.)
%
% HISTORY:
%   2004.03.12 RFD (bob@white.stanford.edu) shamelessly copied the core
%   algorithm from spm99's spm_imatrix.m.
%   2005.07.08 ras (sayres at stanford edu) imported into mrVista 2.0

% Translations and zooms
%-----------------------------------------------------------------------
R         = A(1:3,1:3);
C         = chol(R'*R);
if(size(A,2)>3)
    trans = A(1:3,4)';
else
    trans = [0 0 0];
end
scale = diag(C)';
if det(R)<0, scale(1)=-scale(1);end % Fix for -ve determinants

% Shears
%-----------------------------------------------------------------------
C         = diag(diag(C))\C;
skew      = C([4 7 8]);
R0        = affineBuild([0 0 0], [0 0 0], scale, skew);
R0        = R0(1:3,1:3);
R1        = R/R0;

% This just leaves rotations in matrix R1
%-----------------------------------------------------------------------
%[          c5*c6,           c5*s6, s5]
%[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
%[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]

rot(2) = asin(rang(R1(1,3)));
if (abs(rot(2))-pi/2).^2 < 1e-9,
	rot(1) = 0;
	rot(3) = atan2(-rang(R1(2,1)), rang(-R1(3,1)/R1(1,3)));
else,
	c    = cos(rot(2));
	rot(1) = atan2(rang(R1(2,3)/c), rang(R1(3,3)/c));
	rot(3) = atan2(rang(R1(1,2)/c), rang(R1(1,1)/c));
end;

return;

% There may be slight rounding errors making b>1 or b<-1.
function a = rang(b)
a = min(max(b, -1), 1);
return;


text:
function rx = rxOut(rx);
%
% rx = rxOut(rx);
%
% For mrRx, nudge the prescription 'out' (along the vector pointing from
% the middle of the first slice toward the middle of the last slice) a
% certain delta according to the value of the nudge slider.
%
% ras, 08/05.
if ieNotDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

%%%%%get params
nudgeVal = get(rx.ui.nudge.sliderHandle,'Value')/6; % make a small nudge

%%%%%compute a vector pointing 'in'
% middle X, Y coords
midX = rx.rxDims(2)/2;
midY = rx.rxDims(1)/2;

% cols of coords are: middle of last slice, middle of first slice
coords(:,1) = [midX; midY; 1];
coords(:,2) = [midX; midY; rx.rxDims(3)];

% convert coords into volume coordinate space
volCoords = rx2vol(rx,coords);

% get vector as diff. b/w these two points
vec = diff(volCoords');

%%%%%modify translation accordingly
[trans rot scale skew] = affineDecompose(rx.xform);
trans = trans + nudgeVal.*vec;
newXform = affineBuild(trans,rot,scale,skew);
rx = rxSetXform(rx,newXform);

return


text:
filename='test2d_micro';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'horizontal';
cost={'chomog_alphabeta'};
weights=[1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 1 0]';
micro.beta =[1 1 0]';

% For all tests
plotting = false;
printing = false;
monitoring = false;
maxiter = 3;


text:
function [jiDim,jiUL,resolu,ZC,bands]=envihdrread(filename)
% any ENVI file uint 8, int16, and uint16
% [jiDim,jiUL,resolu,ZC]=envihdrread(filename)
% Input file directory
% Output 1) im = data
% Output 2) jiDim dimension cols and rows
% Output 3)jiUL UpperLeft coner of the pixel (x,y)
% Output  4) resolution [x,y]
% Output 5) ZC zone code

filename_HDR=[filename,'.HDR'];
filename_hdr=[filename,'.hdr'];

fid_in1=fopen(filename_hdr,'r');
fid_in2=fopen(filename_HDR,'r');

if fid_in1~=-1
    fid_in=fid_in1;
elseif fid_in2~=-1
    fid_in=fid_in2;
else
    fprintf('Wrong ENVI header file!\n');
    fprintf('%s\n',filename); % show envi hdr file
    return;
end

geo_char=fscanf(fid_in,'%c',inf);
fclose(fid_in);

geo_char=geo_char';
geo_str=strread(geo_char,'%s');

indx_samples = strmatch('samples',geo_str)+2;
indx_lines = strmatch('lines',geo_str)+2;
indx_bands = strmatch('bands',geo_str)+2;
indx_datatype = strmatch('data',geo_str)+3;
indx_interleave = strmatch('interleave',geo_str)+2;
indx_xUL = strmatch('map',geo_str)+6;
indx_yUL = strmatch('map',geo_str)+7;
indx_xreso = strmatch('map',geo_str)+8;
indx_yreso = strmatch('map',geo_str)+9;
indx_zc = strmatch('map',geo_str)+10;

% read input image hdr
cols = str2double(geo_str(indx_samples));
rows = str2double(geo_str(indx_lines));
jiDim = [cols,rows];

bands = str2double(geo_str(indx_bands));
datatype = str2double(geo_str(indx_datatype));
interleave = char(geo_str(indx_interleave));
jiUL(1)=str2double(geo_str(indx_xUL));
jiUL(2)=str2double(geo_str(indx_yUL));
resolu(1)=str2double(geo_str(indx_xreso));
resolu(2)=str2double(geo_str(indx_yreso));
ZC=str2double(geo_str(indx_zc));

if datatype == 1
    in_type = 'uint8';
elseif datatype == 2
    in_type = 'int16';
elseif datatype == 3
    in_type = 'int32';
elseif datatype == 4
    in_type = 'single';
elseif datatype == 12
    in_type = 'uint16';
else
    error('Invalid read data type!');
end

% % Define the data set.
% % Read every other band of the data using the Band-Sequential format.
% im = multibandread(filename, [rows cols bands], ...
%                     [in_type,'=>',in_type], 0, interleave, 'ieee-le');
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
filename='Gripping_tetrahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
classdef MimToolList < handle
    % MimToolList. Stores a list of MimTools
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        Tools % Map of tools
    end

    properties (Access = private)
        ToolCallback
        ViewerPanel

        CineTool
        WindowLevelTool
        ZoomTool
        PanTool
        PanMatlabTool
        Rotate3dMatlabTool
        ZoomMatlabTool
        EditTool
        ReplaceColourTool
        MapColourTool
        MarkerPointTool
    end

    methods
        function obj = MimToolList(marker_manager, tool_callback, viewer_panel, image_parameters, image_display_parameters)
            obj.ToolCallback = tool_callback;
            obj.ViewerPanel = viewer_panel;

            obj.CineTool = MimCineTool(image_parameters, tool_callback);
            obj.WindowLevelTool = MimWindowLevelTool(image_display_parameters, tool_callback, viewer_panel);
            obj.ZoomTool = MimZoomTool(image_parameters, tool_callback);
            obj.PanTool = MimPanTool(tool_callback);
            obj.PanMatlabTool = MimPanMatlabTool(tool_callback);
            obj.Rotate3dMatlabTool = MimRotate3dMatlabTool(tool_callback);
            obj.ZoomMatlabTool = MimZoomMatlabTool(tool_callback);
            obj.EditTool = MimEditManager(obj.ViewerPanel);
            obj.ReplaceColourTool = MimReplaceColourTool(obj.ViewerPanel);
            obj.MapColourTool = MimMapColourTool(obj.ViewerPanel);
            obj.MarkerPointTool = MimMarkerPointTool(marker_manager, obj.ViewerPanel);

            tool_list = {obj.ZoomMatlabTool, obj.PanMatlabTool, obj.Rotate3dMatlabTool, obj.MarkerPointTool, obj.WindowLevelTool, obj.CineTool, obj.EditTool, obj.ReplaceColourTool, obj.MapColourTool};

            obj.Tools = containers.Map;
            for tool_set = tool_list
                tool = tool_set{1};
                tag = tool.Tag;
                obj.Tools(tag) = tool;
            end
        end

        function tool = GetCurrentToolForSelectedControl(obj, mouse_is_down, keyboard_modifier)

            tool = obj.GetCurrentTool(mouse_is_down, keyboard_modifier, obj.ViewerPanel.SelectedControl);
        end

        function tool = GetCurrentTool(obj, mouse_is_down, keyboard_modifier, selected_control)
            if ~isempty(keyboard_modifier) && ~isempty(mouse_is_down) && mouse_is_down
                if strcmp(keyboard_modifier, 'extend')
                    tool = obj.PanTool;
                    return;
                elseif strcmp(keyboard_modifier, 'alt')
                    tool = obj.ZoomTool;
                    return;
                end
            end
            tool = obj.Tools(selected_control);
        end

        function tool = GetTool(obj, tag)
            tool = obj.Tools(tag);
        end

        function UpdateTools(obj)
            tool_list = obj.Tools.values;

            % Disable inactive tools before enabling new tools
            for tool_set = tool_list
                tool = tool_set{1};
                if ~strcmp(obj.ViewerPanel.SelectedControl, tool.Tag);
                    tool.Enable(false);
                end
            end

            % Enable tools
            for tool_set = tool_list
                tool = tool_set{1};
                if strcmp(obj.ViewerPanel.SelectedControl, tool.Tag)
                    tool.Enable(true);
                    % Set the context menu on the images. Note it would be hidden by
                    % the images if it were set on the axes
                    obj.ViewerPanel.BackgroundLayer.SetContextMenu(tool.GetContextMenu);
                    obj.ViewerPanel.SegmentationLayer.SetContextMenu(tool.GetContextMenu);
                end
            end
        end

        function input_has_been_processed = ShortcutKeys(obj, key, selected_control)
            % Each tool has a shortcut key to select it
            for tool = obj.Tools.values
                if strcmpi(key, tool{1}.ShortcutKey)
                    obj.ViewerPanel.SetControl(tool{1}.Tag);
                    input_has_been_processed = true;
                    return
                end
            end

            % Otherwise let the currently selected tool process shortcuts
            input_has_been_processed = obj.Tools(selected_control).Keypressed(key);
        end

        function ImageChanged(obj)
            tool = obj.GetCurrentPrimaryTool;
            if ~isempty(tool) && tool.Enabled
                tool.ImageChanged;
            end
        end

        function NewOrientation(obj)
            tool = obj.GetCurrentPrimaryTool;
            if ~isempty(tool) && tool.Enabled
                tool.NewOrientation;
            end
        end

        function OverlayImageChanged(obj)
            tool = obj.GetCurrentPrimaryTool;
            if ~isempty(tool) && tool.Enabled
                tool.OverlayImageChanged;
            end
        end

        function NewSlice(obj)
            tool = obj.GetCurrentPrimaryTool;
            if ~isempty(tool) && tool.Enabled
                tool.NewSlice;
            end
        end

    end

    methods (Access = private)
        function tool = GetCurrentPrimaryTool(obj)
            selected_control = obj.ViewerPanel.SelectedControl;
            if obj.Tools.isKey(selected_control)
                tool = obj.Tools(selected_control);
            else
                tool = [];
            end
        end
    end
end


text:
function res = penalized_coefficients(lasso_object, indices, alignment)
% -------------------------------------------------------------------------
% function res = penalized_coefficients(penalized_object, indices, alignment)
% -------------------------------------------------------------------------
% PURPOSE: This function returns the LASSO coefficients for which the
%          chosen index equals those listed in indices.
%          The alignment index can be chosen to be:
%          lambdas, normalized L1 norms and L1 norms;
% -------------------------------------------------------------------------
% INPUTS:
% lasso_object: a structure as the one returned by the lasso routine:
% indices:      index of the points at which interpolation is wanted
% alignment:    one string of ('normalized_l1', 'l1'm 'lambdas')
%               indicating which index to use for alignment
% -------------------------------------------------------------------------
% OUTPUTS:
% res:          a structure containing  the coefficients and the
%               regularization parameter at the interpolated points
%               intercept:
%               nbeta:
%               beta:
%               lambda:
% -------------------------------------------------------------------------
% Author: Guilherme V. Rocha
%         Department of Statistics
%         University of California, Berkeley
%         gvrocha@stat.berkeley.edu, gvrocha@gmail.com
% 2006/09
% -------------------------------------------------------------------------
% See also: LASSO, BLASSOL2

% 0. Checking input parameters:
%==========================================================================
if nargin < 3
  alignment = 'normalized_penalty';
end;

% 1. Interpolates according to the chosen L1 norm (normalized x non normalized)
%==========================================================================
switch(lower(alignment))
  case{'normalized_penalty'}
    sizes = lasso_object.npenalty;
  case{'penalty'}
    sizes = lasso_object.penalty;
  case{'lambda'}
    sizes = lasso_object.lambda;
  otherwise
    error('Unrecognized alignment string');
end;

% Takes care of the case when there are entries with repeated L1s:
[unique_sizes, I, J]   = unique(sizes);

% If one of the points to be calculated is beyond the end of the path, use
% the end of path coeffcients for it:
indices       = min(indices, max(sizes));
indices       = reshape(indices, prod(size(indices)), 1);
res.intercept = interp1(unique_sizes, lasso_object.intercept(I), indices);
res.beta      = interp1(unique_sizes, lasso_object.beta(I,:), indices);
res.nbeta     = interp1(unique_sizes, lasso_object.nbeta(I,:), indices);
if(~isempty(lasso_object.lambda))
  res.lambda    = interp1(unique_sizes, lasso_object.lambda(I,:), indices);
end;
res.indices   = indices';

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function s=minmod(a,b)
% The minmod limiter
s=0.5*(sign(a)+sign(b)).*min(abs(a),abs(b));

repo: Aero-matlab


text:
% rmart() - Remove eye artifacts from EEG data using regression with
%           multiple time lags. Each channel is first made mean-zero.
%           After JL Kenemans et al., Psychophysiology 28:114-21, 1991.
%
%   Usage: >> rmart('datafile','outfile',nchans,chanlist,eogchan,[threshold])
% Example: >> rmart('noisy.floats','clean.floats',31,[2:31],7)
%
%   Input:   datafile - input float data file, multiplexed by channel
%            outfile  - name of output float data file
%            nchans   - number of channels in datafile
%            chanlist - indices of EEG channel(s) to process (1,...,nchans)
%            eogchan  - regressing channel indices(s) (1,...,nchans)
%            threshold- abs threshold value to trigger regression {def|0 -> 80}
%
% Output: Writes [length(chanlist),size(data,2)] floats to 'outfile'
%
% Note: Regression epoch length and number of lags are set in the script.
%       Some machines may require a new byte_order value in the script.
%       note that runica() -> icaproj() should give better results! See
%       Jung et al., Psychophysiology 111:1745-58, 2000.
%
% Author: Tzyy-Ping Jung, SCCN/INC/UCSD, La Jolla, 1997

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

% 2-22-97  Tzyy-Ping Jung  CNL/Salk Institute, La Jolla, CA
% 2-24-97  Formatted for ICA package release -Scott Makeig
% 12-10-97 Changed name from rmartifact to rmart for toolbox inclusion -sm
% 12-11-97 Adapted to read/write a float matrix -sm & sw
% 09-14-00 Added comments and help -sm
% 01-25-02 reformated help & license -ad

function rmart(datafile,outfile,nchans,chanlist,eogchan,threshold)

if nargin < 5
   help rmart
   return
end

%
% The following parameters may be fine-tuned for a data set
%
DEF_THRESHOLD = 80; % trigger regression on blocks exceeding this (default)
epoch         = 80; % remove artifacts in successive blocks of this length
nlags         = 40; % perform multiple regression filtering of this length

byte_order    = 'b';% (machine-dependent) byte order code for fopen();
MAKE_MEAN_ZERO = 1  % 1/0 flag removing mean offset from each channel

fprintf('Performing artifact regression on data in %s.\n',datafile);

if nargin<6
   threshold = 0;
end
if threshold == 0,
   threshold = DEF_THRESHOLD;
end
fprintf('Regression threshold %g.\n',threshold);

%
% Read the input data
%
[fid,msg]=fopen(datafile,'r',byte_order); % open datafile
if fid < 3,
  fprintf('rmart() - could not open data file: %s\n',msg);
  exit 1
end
data=(fread(fid,'float'))';
status=fclose('all');
if rem(length(data),nchans) == 0 % check length
   fprintf('rmart() - data length not divisible by %d chans.\n',nchans);
    return
end

data = reshape(data,nchans,length(data)/nchans);
[chans,frames] = size(data);
fprintf('Data of size [%d,%d] read.\n',chans,frames);
eog = data(eogchan,:);
data = data(chanlist,:);
procchans = length(chanlist);

fprintf('Regression epoch length %d frames.\n',epoch);
fprintf('Using %d regression lags.\n',nlags);
if length(eogchan)> 1
  fprintf('Processing %d of %d channels using %d EOG channels.\n',...
                        procchans,chans,length(eogchan));
else
  fprintf('Processing %d of %d channels using EOG channel %d.\n',...
                        procchans,chans,eogchan);
end

%
% Process the data
%
for i=1:procchans
  chan = chanlist(i);
  idx=[];
  frame=1+epoch/2+nlags/2;
  if MAKE_MEAN_ZERO
    data(chan,:) = data(chan,:) - mean(data(chan,:)); % make mean-zero
  end

  % Search the EOG & EEG records for values above threshold,
  % Selected frame numbers are registered in the variable "idx".
  % The entries in "idx" are at least epoch apart to avoid double
  % compensation (regression) on the same portion of the EEG data.

  while frame <= length(eog)-epoch/2-nlags/2,  % foreach epoch in channel
    stop = min(frame+epoch-1,eogframes);
    tmp= ...
        find( abs(eog(frame:stop)) >= threshold ...
             | abs(data(chan,frame:stop)) >= threshold);
                                           % find beyond-threshold values
    if length(tmp) ~= 0
      mark = tmp(1)+frame-1;
      if  length(idx) ~= 0
        if mark-idx(length(idx))  < epoch,
           idx=[idx idx(length(idx))+epoch]; % To guarantee idx(i) & idx(i-1)
                                             % are at least EPOCH points apart
           frame = idx(length(idx))+epoch/2;
        else
           idx=[idx mark];
           frame = mark + epoch/2;
        end
      else
        idx=[idx mark];
        frame = mark + epoch/2;
      end
    else
      frame=frame+epoch;
    end
  end % while

  % For each registered frame, take "epoch" points
  % surrounding it from the EEG, and "epoch + lag" points
  % from the EOG channel. Then perform multivariate
  % linear regression on EEG channel.

  for j=1:length(idx);
     art=ones(1,epoch);
     eogtmp=eog(idx(j)-epoch/2-nlags/2:idx(j)+epoch/2-1+nlags/2);

     % Collect EOG data from lag/2 points before to lag/2 points
     % after the regression window.
     for J=nlags:-1:1,
       art=[art ; eogtmp(J:J+epoch-1)];
     end
     eegtmp=data(chan,idx(j)-epoch/2:idx(j)+epoch/2-1);

     eegeog=eegtmp*art';        % perform the regression here
     eogeog=art*art';
     b=eegeog/eogeog;
     eegtmp=eegtmp-b*art;
     data(chan,idx(j)-epoch/2:idx(j)+epoch/2-1)=eegtmp;
   end % j
end % i

%
% Write output file
%
[fid,msg]=fopen(outfile,'w',byte_order);
if fid < 3
  fprintf('rmart() - could not open output file: %s\n',msg);
  return
end
count = fwrite(fid,data,'float');
if count == procchans*frames,
  fprintf('Output file "%s" written, size = [%d,%d] \n\n',...
             outfile,procchans,frames);
else
  fprintf('rmart(): Output file "%s" written, SIZE ONLY [%d,%g]\n',...
             outfile,procchans,count/procchans);
end
fclose('all');

author: PatternRecognition

repo: OpenBMI


text:
function [theta,b,QMin]=LTSRegression(x,y,h)
%%LTSREGRESSION Perform least trimmed squares (LTS) regression using the
%               fast, approximate algorithm of [1] that is identified as
%               best for <600 data points. This means find theta and b that
%               minimizes
%               sum_{i=1}^h r_i^2
%               where r_i is the ith largest value of
%               abs(y(i)-theta*x(:,i)-b)
%               This type of regression is robust to outliers.
%
%INPUTS: x,y These are the real samples of data. x is an xDimXN set of N
%            vectors and y is a set of corresponding scalar values. The
%            regression searches for theta and b to fit a line such that
%            y=theta'*x(:,i)+b.
%          h The trim value of the estimator. This must be between
%            fix((N+xDim+2)/2) and n. The default if omitted or an empty
%            matrix is passed is the value that allows for the highest
%            breakdown value, which is fix((N+xDim+2)/2).
%
%OUTPUTS: theta An xDimX1 vector. The multivariate slope of the line.
%             b The scalar y intercept of the data.
%          QMin The residual of the data.
%
%The algorithm described in [1] without the y intercept adjustment is used.
%
%EXAMPLE:
%This is similar to an example in [1]. There is a distribution that one
%might expect to be well approximated with a line corrupted by a minority
%of samples of a different distribution at the side. The plot is zoomed in
%so not all samples are shown. The standard linear regression result is
%shown for comparison.
% x=100*randn(1,800/2);
% y=x+1+randn(1,800/2);
%
% xy=GaussianD.rand(200/2,[50;0],25*eye(2,2));
% x=[x,xy(1,:)];
% y=[y,xy(2,:)];
%
% figure(1)
% clf
% hold on
% scatter(x,y,'b','filled')
%
% coeffs=multilinRegress([x;y]);
% numPoints=100;
% xEst=linspace(-200,200,numPoints);
% yEst=xEst*coeffs(1)+coeffs(2);
% plot(xEst,yEst,'--k','linewidth',2)
% [theta,b]=LTSRegression(x,y);
% yEst=theta*xEst+b;
% plot(xEst,yEst,'-r','linewidth',2)
% axis([-30-20,70+20,-30-20,50+20])%Zoom in.
% legend('Samples','Least Squares Fit','LTS Fit','location','NorthWest')
% h1=xlabel('x');
% h2=ylabel('y');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] P. J. Rousseeuw and V. Katrien, "Computing LTS regression for large
%    data sets," Data Mining and Knowledge Discovery, vol. 12, no. 1, pp.
%    26-45, Jan. 2006.
%
%August 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(x,2);

x=[x;ones(1,n)];

y=y(:);

p=size(x,1);%p>=2

if(nargin<3||isempty(h))
    h=fix((n+p+1)/2);
end

if(h<fix((n+p+1)/2)||h>n)
   error('Invalid value of h provided.')
end

totalStarts=binomial(n,p);
Q3Heap=BinaryHeap(10,false);
for k=1:min(totalStarts,500)
    if(totalStarts<500)
        %Do all possible initializations.
        idx=unrankCombination(k-1,n,p)+1;
    else
        %Get a random subset of points in x.
        idx=randCombination(n,p)+1;
    end

    xCur=x(:,idx);
    yCur=y(idx);
    %If what we get is not full rank, then we just use a
    %pseudoinverse. This skips issues of determining which and how
    %many additional columns to add.
    if(rank(xCur)<p)
        theta0=pinv(xCur)'*yCur;
    else
        theta0=xCur'\yCur;
    end

    for i=1:3
        r0=y'-sum(bsxfun(@times,x,theta0),1);
        [~,idx]=sort(abs(r0),'ascend');
        idx=idx(1:h);

        xCur=x(:,idx);
        yCur=y(idx);
        theta0=pinv(xCur)'*yCur;
    end

    r2=(y'-sum(bsxfun(@times,x,theta0),1)).^2;
    Q3=sum(mink(r2,h));

    if(Q3Heap.heapSize()<10)
        Q3Heap.insert(Q3,idx);
    else
        topEl=Q3Heap.getTop();
        if(topEl.key>Q3)
            Q3Heap.deleteTop();
            Q3Heap.insert(Q3,idx);
        end
    end
end

%For the 10 Q3 values in the heap, iterate the C-step in Section 2
%until convergence.
theta0List=zeros(p,10);
QList=zeros(10,1);
for k=1:10
    topEl=Q3Heap.getTop();
    Q3Heap.deleteTop();

    idx=topEl.value;
    xCur=x(:,idx);
    yCur=y(idx);
    theta0=pinv(xCur)'*yCur;
    idxOld=zeros(h,1);

    %It was noted that convergence should be for fewer than 10
    %iterations, so we make the maximum 10.
    for curIter=1:10
        if(all(idx(:)==idxOld(:)))
            break;
        end
        r0=y'-sum(bsxfun(@times,x,theta0),1);
        [~,idx]=sort(abs(r0),'ascend');
        idx=idx(1:h);

        xCur=x(:,idx);
        yCur=y(idx);
        theta0=pinv(xCur)'*yCur;
    end

    r2=(y'-sum(bsxfun(@times,x,theta0),1)).^2;
    Q=sum(mink(r2,h));

    QList(k)=Q;
    theta0List(:,k)=theta0;
end

[QMin,idx]=min(QList);
theta0Min=theta0List(:,idx);

theta=theta0Min(1:(p-1));
b=theta0Min(p);
end

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


text:
function [Torr] = mbar2Torr(mbar)
% Convert pressure from millibars to torr.
% Chad Greene 2012
Torr = mbar*0.750062;

author: Sable

repo: mcbench-benchmarks


text:
%% RECIPROCAL VO/FORMATION CONTROL AGENT (agent_formation_RVO.m) %%%%%%%%%%
% This agent is a hybrid of the formation control element agent class and
% the Reciprocal Velocity Obstacle (VO) class.

% Author: James A. Douthwaite

classdef agent_formation_RVO < agent_formation_VO & agent_RVO
    %% INITIALISE THE AGENT SPECIFIC PARAMETERS
    properties

    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function obj = agent_formation_RVO(varargin)
            % CALL THE SUPERCLASS CONSTRUCTOR
            obj = obj@agent_formation_VO(varargin);

            % //////////////// Check for user overrides ///////////////////
            [obj] = obj.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % NOTES:
        % This class simply overrides the 'process_timeCycle' function of
        % the superclass 'agent_formation_VO' such that it considers the
        % 'avoidanceCorrection' of the second super class 'agent_RVO'.
    end
end
% AGENT STATE VECTOR [x;y;z;phi;theta;psi;xdot;ydot;zdot;phidot;thetadot;psidot]

repo: OpenMAS


text:
function J=vl_imwhiten(I,alpha,cutoff)
% VL_IMWHITEN  Whiten an image
%   J = VL_IMWHITEN(I,ALPHA) approximatively whitens the power spectrum
%   of the natural image I. The algorithm assumes that the modulus of
%   the spectrum decays as 1/f^ALPHA (f is the frequency).
%
%   VL_IMWHITEN(I) uses ALPHA=1 (a typical value for natural images).
%
%   VL_IMWHITEN(I,ALPHA,CUTOFF) also applies a low-pass filter with
%   cutoff frequency equal to CUTOFF x FN, where FN is the Nyquist
%   frequency (half of the sampling frequency).
%
%   See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if ~exist('alpha','var'),  alpha = 1 ; end
if ~exist('cutoff','var'), cutoff = [] ; end

[M,N]=size(I) ;

% Frequency domain
fn = 0.5 ; % Nyquist freq (=1/2T, T=1)
fx_range=linspace(-fn, fn, N) ;
fy_range=linspace(-fn, fn, M) ;
[fx fy]=meshgrid(fx_range, fy_range) ;

% Whitening filter
rho=sqrt(fx.*fx+fy.*fy);
filt=rho.^alpha ;

% Low-pass filter
if ~isempty(cutoff)
  fcut = cutoff * fn ;
  filt = filt .* exp(-(rho/fcut).^4);
  %filt = filt .* exp( - 0.5 * (rho / fcut) .^ 2);
end

% Apply filter
J = real(ifft2(fft2(I).*fftshift(filt))) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
%% data association for every kalman filter (every object from previous frame)
function [k, ocn, szn] = asc(k, ocn, szn)
for i             = 1 : size(k, 2)                      % for every kalman filter
[k(i).s, idx]     = assoc(k(i).s, ocn);                 % apply data association
k(i).sz           = szn(:, idx);
szn(:, idx)       = [];
ocn(:, idx)       = [];                                 % eliminated checked objects from the new frame and go for
end                                                     % associating remained objects to the next kalman filter
end
%% data association: 1. no object, 2. one object, 3. multiple objects
function [s, idx] = assoc(s, ocn)
can    = ocn(1 : 2, :);                                 % candidates
pe     = [s.x(1); s.x(3)];                              % previous estimate (to do: replace with last prediction!)
gate   = 4;                                             % gating: 7
idx    = ((sum((can - repmat(pe, 1, size(can, 2)))...   % indexes of objects inside the gate
         .^2)) .^0.5) < gate;
if     sum(idx) == 0                                    % 1. no object
s.z    = pe;                                            % previous estimate/prediction!
elseif sum(idx) == 1                                    % 2. one object
s.z    = can(:, idx);                                   % associated object
elseif sum(idx) > 1                                     % 3. multiple objects (take nearest object)
idx    = ((sum((can - repmat(pe, 1, size(can, 2)))...   % index of the nearest object
         .^2)) .^0.5) == min((sum((can - repmat(pe...
         , 1, size(can, 2))) .^2)) .^0.5);
s.z    = can(:, idx);                                   % associated object
end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function varargout = plot(f, varargin)
%PLOT   Basic linear plot for TRIGTECH objects.
%   PLOT(F) plots the TRIGTECH object F.
%
%   PLOT(F, S) allows various line types, plot symbols, and colors to be used
%   when S is a character string made from one element from any or all the
%   following 3 columns:
%
%            b     blue          .     point              -     solid
%            g     green         o     circle             :     dotted
%            r     red           x     x-mark             -.    dashdot
%            c     cyan          +     plus               --    dashed
%            m     magenta       *     star             (none)  no line
%            y     yellow        s     square
%            k     black         d     diamond
%            w     white         v     triangle (down)
%                                ^     triangle (up)
%                                <     triangle (left)
%                                >     triangle (right)
%                                p     pentagram
%                                h     hexagram
%
%   The entries from the centre columns are plotted at TRIGPTS being
%   used to represent F. If no options from this column are chosen, 'o' is
%   chosen by default if LENGTH(F) < 256.
%
%   The X,Y pairs, or X,Y,S triples, can be followed by parameter/value pairs to
%   specify additional properties of the lines. For example,
%            f = trigtech.constructor(@sin);
%            plot(f, 'LineWidth', 2, 'Color', [.6 0 0])
%   will create a plot with a dark red line width of 2 points.
%
%   H1 = PLOT(F) returns a column vector of handles to lineseries objects, one
%   handle per plotted line (in the case of array-valued TRIGTECH objects).
%   [H1, H2] returns a second vector of column handles, this time for each of
%   the marker plots.
%
% See also PLOT3, PLOTDATA.

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

% Deal with an empty input:
if ( isempty(f) )
    if ( nargout == 1 )
        varargout{1} = plot([]);
    end
    return
end

% Store the hold state of the current axis:
holdState = ishold;

%%
% Get the data for plotting from PLOTDATA():
if ( nargin > 1 && isa(varargin{1}, 'trigtech') )
    % Deal with plot(f, g);
    g = varargin{1};
    varargin(1) = [];
    % We can only plot real against real:
    f = real(f);
    g = real(g);
    % Call PLOTDATA():
    data = plotData(f, g);
else
    g = [];
    data = plotData(f);
end

%%
% Plot the curve:
if ( isreal(f) )
    h1 = plot(real(data.xLine), real(data.yLine), varargin{:});
else
    h1 = plot(data.yLine, varargin{:});
end
set(h1, 'Marker', 'none')
hold on

%%
% Plot the points:
if ( isreal(f) )
    h2 = plot(real(data.xPoints), real(data.yPoints), varargin{:});
else
    h2 = plot(data.yPoints, varargin{:});
end

% Change the style accordingly:
set(h2,'LineStyle', 'none')
if ( all(strcmp(get(h2, 'Marker'), 'none')) && (length(f) < 257) )
    set(h2, 'Marker', 'o')
end

%%

% Return hold state to what it was before:
if ( ~holdState )
    hold off
end

% Give an output if one was requested:
if ( nargout > 0 )
    varargout{1} = h1;
    varargout{2} = h2;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [yc] = checkCosts(obj, x, output_file,permission)
    % Check the value of const function

    if nargin > 2
        % print to the file
        if nargin < 4
            permission = 'w';
        else
            validatestring(permission, {'a','w'});
        end
        f_id = fopen(output_file, permission);
    else
        % print on the screen
        f_id = 1;
    end

    fprintf(f_id, '**************************************************\n');
    fprintf(f_id, 'Checking cost value of %s:\n', obj.Name);

    cost_table = obj.CostTable;
    [n_node, n_cost] = size(cost_table);
    yc = zeros(1,n_cost);
    fprintf(f_id, '**************************************************\n');
    fprintf(f_id,'%12s \t %12s\n','Cost','Value');

    for j=1:n_cost
        cost_name = obj.CostTable.Properties.VariableNames{j};
        cost_array = obj.CostTable.(cost_name);
        fprintf(f_id, 'Cost: %s \n', cost_name);
        for k=1:n_node
            cost = cost_array(k);
            if cost.Dimension ~=0
                dep_constr = getSummands(cost);
                for ll = 1:numel(dep_constr)
                    dep_var = dep_constr(ll).DepVariables;
                    var = arrayfun(@(v)x(v.Indices(:)),dep_var,'UniformOutput',false); % dependent variables
                    if isempty(dep_constr(ll).AuxData)
                        yc(j) = yc(j) + feval(dep_constr(ll).Funcs.Func, var{:});
                    else
                        yc(j) = yc(j) + feval(dep_constr(ll).Funcs.Func, var{:}, dep_constr(ll).AuxData{:});
                    end

                end
            end
        end
        fprintf(f_id,'%12s \t %12.8E\n',cost_name,yc(j));
    end
    fprintf(f_id, '**************************************************\n');

    if f_id ~= 1
        fclose(f_id);
    end
end

author: ayonga

repo: frost-dev

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


text:
classdef (Abstract) AbstractSimulator < handle
% Abstract class for movement simulation.
% Author: Yu Jiang
% Contact: yu.jiang@nyu.edu
% Copyright 2015 Yu Jiang

    properties (Constant)
        % Constant simulation parameters
        tau = 0.05;   % Time constant
        m1 = 2;       % mass on x direction
        m2 = 2;       % mass on y direction

        c1 = 0.15/2;  % noise scale
        c2 = 0.05/2;  % noise scale
        dt_ = 0.005;   % sample time for learning

        Q0 = [500 0; 0 1000];   % Initial weighting matrices
        R = diag([0.01,0.01]);  % Initial weighting matrices

        % The null filed dynamics
        A0 = [zeros(2) eye(2) zeros(2);
              zeros(2,4) diag([1/AbstractSimulator.m1 1/AbstractSimulator.m2])
              zeros(2,4) -diag([1/AbstractSimulator.tau 1/AbstractSimulator.tau])];
        B =  [zeros(4,2); diag([1/AbstractSimulator.tau 1/AbstractSimulator.tau])];
    end

    properties
        % Variables
        dt % Sample time for simulation
        A  % Dynamics with Force field
    end

    properties
        % Others
        fig1;
        fig2;
        fig3;
    end
end

author: yu-jiang

repo: radpbook

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


text:
function lrdata = getlrdata(F)
F = flatten(F);
k = 1;
lrdata = [];
for i = 1:length(F.clauses)
    if F.clauses{i}.type == 14
        lrdata{end+1} =  F.clauses{i}.data;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
clear,clc
% compute the background image
Imzero = zeros(240,320,3);
for i = 1:5
Im{i} = double(imread(['DATA/',int2str(i),'.jpg']));
Imzero = Im{i}+Imzero;
end
Imback = Imzero/5;
[MR,MC,Dim] = size(Imback);

% Kalman filter initialization
R=[[0.2845,0.0045]',[0.0045,0.0455]'];
H=[[1,0]',[0,1]',[0,0]',[0,0]'];
Q=0.01*eye(4);
P = 100*eye(4);
dt=1;
A=[[1,0,0,0]',[0,1,0,0]',[dt,0,1,0]',[0,dt,0,1]'];
g = 6; % pixels^2/time step
Bu = [0,0,0,g]';
kfinit=0;
x=zeros(100,4);

% loop over all images
for i = 1 : 60
  % load image
  Im = (imread(['DATA/',int2str(i), '.jpg']));
  imshow(Im)
  imshow(Im)
  Imwork = double(Im);

  %extract ball
  [cc(i),cr(i),radius,flag] = extractball(Imwork,Imback,i);
  if flag==0
    continue
  end

  hold on
    for c = -1*radius: radius/20 : 1*radius
      r = sqrt(radius^2-c^2);
      plot(cc(i)+c,cr(i)+r,'g.')
      plot(cc(i)+c,cr(i)-r,'g.')
    end
  % Kalman update
i
  if kfinit==0
    xp = [MC/2,MR/2,0,0]'
  else
    xp=A*x(i-1,:)' + Bu
  end
  kfinit=1;
  PP = A*P*A' + Q
  K = PP*H'*inv(H*PP*H'+R)
  x(i,:) = (xp + K*([cc(i),cr(i)]' - H*xp))';
  x(i,:)
  [cc(i),cr(i)]
  P = (eye(4)-K*H)*PP

  hold on
    for c = -1*radius: radius/20 : 1*radius
      r = sqrt(radius^2-c^2);
      plot(x(i,1)+c,x(i,2)+r,'r.')
      plot(x(i,1)+c,x(i,2)-r,'r.')
    end
      pause(0.3)
end

% show positions
  figure
  plot(cc,'r*')
  hold on
  plot(cr,'g*')
%end

%estimate image noise (R) from stationary ball
  posn = [cc(55:60)',cr(55:60)'];
  mp = mean(posn);
  diffp = posn - ones(6,1)*mp;
  Rnew = (diffp'*diffp)/5;


text:
%MDL_STANFORD Create model of Stanford arm
%
%      mdl_stanford
%
% Script creates the workspace variable stanf which describes the
% kinematic and dynamic characteristics of the Stanford (Scheinman) arm.
%
% Also defines the vectors:
%   qz   zero joint angle configuration.
%
% Note::
% - SI units are used.
% - Gear ratios not currently known, though reflected armature inertia
%   is known, so gear ratios are set to 1.
%
% References::
% - Kinematic data from "Modelling, Trajectory calculation and Servoing of
%   a computer controlled arm".  Stanford AIM-177.  Figure 2.3
% - Dynamic data from "Robot manipulators: mathematics, programming and control"
%   Paul 1981, Tables 6.5, 6.6
% - Dobrotin & Scheinman, "Design of a computer controlled manipulator for
%   robot research", IJCAI, 1973.
%
% See also SerialLink, mdl_puma560, mdl_puma560akb.

% MODEL: Stanford, Stanford Arm, prismatic, 6DOF, standard_DH

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

clear L
%             th    d       a    alpha
L(1) = Link([ 0     0.412   0   -pi/2     0]);
L(2) = Link([ 0     0.154   0    pi/2     0]);
L(3) = Link([ -pi/2 0       0    0        1]);  % PRISMATIC link
L(4) = Link([ 0     0       0   -pi/2     0]);
L(5) = Link([ 0     0       0    pi/2     0]);
L(6) = Link([ 0     0.263   0    0        0]);

% guestimates of some parameters
%
% According to the IJCAI paper the rack is 38in and the maximum reach is 52in
% From the image http://infolab.stanford.edu/pub/voy/museum/pictures/display/robots/IMG_2408ArmCenter.JPG
% and scaled by the rack length (38in) it looks like the minimum stroke is 12in.
%
L(3).qlim = [12 12+38] * 0.0254;

% According to the IJCAI paper
L(1).qlim = [-170 170]*pi/180;
L(2).qlim = [-170 170]*pi/180;
L(4).qlim = [-170 170]*pi/180;
L(5).qlim = [-90 90]*pi/180;
L(6).qlim = [-170 170]*pi/180;

L(1).m = 9.29;
L(2).m = 5.01;
L(3).m = 4.25;
L(4).m = 1.08;
L(5).m = 0.630;
L(6).m = 0.51;

L(1).Jm = 0.953;
L(2).Jm = 2.193;
L(3).Jm = 0.782;
L(4).Jm = 0.106;
L(5).Jm = 0.097;
L(6).Jm = 0.020;

L(1).G = 1;
L(2).G = 1;
L(3).G = 1;
L(4).G = 1;
L(5).G = 1;
L(6).G = 1;

L(1).I = [0.276   0.255   0.071   0   0   0];
L(2).I = [0.108   0.018   0.100   0   0   0];
L(3).I = [2.51    2.51    0.006   0   0   0 ];
L(4).I = [0.002   0.001   0.001   0   0   0 ];
L(5).I = [0.003   0.0004  0       0   0   0];
L(6).I = [0.013   0.013   0.0003  0   0   0];

L(1).r = [0    0.0175 -0.1105];
L(2).r = [0   -1.054  0];
L(3).r = [0    0      -6.447];
L(4).r = [0    0.092  -0.054];
L(5).r = [0    0.566   0.003];
L(6).r = [0    0       1.554];

qz = [0 0 0 0 0 0];

stanf = SerialLink(L, 'name', 'Stanford arm');
stanf.plotopt = {'workspace', [-2 2 -2 2 -2 2]};
stanf.model3d = 'example/stanford';

author: Allopart

repo: rbpf-gmapping


text:
function res=cosmo_distatis(ds, varargin)
% apply DISTATIS measure to each feature
%
% res=cosmo_statis_measure(ds, opt)
%
% Inputs:
%    ds               dataset struct with dissimilarity values; usually
%                     the output from @cosmo_dissimilarity_matrix_measure
%                     applied to each subject followed by cosmo_stack. It
%                     can also be a cell with datasets (one per subject).
%    'return', d      d can be 'distance' (default) or 'crossproduct'.
%                     'distance' returns a distance matrix, whereas
%                     'crossproduct' returns a crossproduct matrix
%    'split_by', s    sample attribute that discriminates chunks
%                     (participants) (default: 'chunks')
%    'shape', sh      shape of output if it were unflattened using
%                     cosmo_unflatten, either 'square' (default) or
%                     'triangle' (which gives the lower diagonal of the
%                     distance matrix)
%
% Returns:
%    res              result dataset struct with feature-wise optimal
%                     compromise distance matrix across subjects
%      .samples
%
%
% Example:
%     % (This example cannot be documentation tested using Octave,
%     %  since Octave does not allow for-loops with evalc)
%     cosmo_skip_test_if_no_external('matlab');
%     %
%     ds=cosmo_synthetic_dataset('nsubjects',5,'nchunks',1,'ntargets',4);
%     %
%     % define neighborhood (here a searchlight with radius of 1 voxel)
%     nbrhood=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
%     %
%     % define measure
%     measure=@cosmo_dissimilarity_matrix_measure;
%     % each subject is a chunk
%     ds.sa.chunks=ds.sa.subject;
%     % compute DSM for each subject
%     sp=cosmo_split(ds,'chunks');
%     for k=1:numel(sp)
%         sp{k}=cosmo_searchlight(sp{k},nbrhood,measure,'progress',false);
%         sp{k}.sa.chunks=ones(6,1)*k;
%     end
%     % merge results
%     dsms=cosmo_stack(sp);
%     %
%     r=cosmo_distatis(dsms,'return','distance','progress',false);
%     cosmo_disp(r);
%     %|| .samples
%     %||   [     0         0         0         0         0         0
%     %||     0.818      1.09      0.77     0.653      1.03     0.421
%     %||     0.869       1.3      1.06      1.04     0.932      1.07
%     %||       :         :         :         :         :         :
%     %||      1.16     0.889      0.99     0.631      1.48     0.621
%     %||     0.268     0.952     0.965     0.462     0.943      1.04
%     %||         0         0         0         0         0         0 ]@16x6
%     %|| .fa
%     %||   .center_ids
%     %||     [ 1         2         3         4         5         6 ]
%     %||   .i
%     %||     [ 1         2         3         1         2         3 ]
%     %||   .j
%     %||     [ 1         1         1         2         2         2 ]
%     %||   .k
%     %||     [ 1         1         1         1         1         1 ]
%     %||   .nvoxels
%     %||     [ 3         4         3         3         4         3 ]
%     %||   .radius
%     %||     [ 1         1         1         1         1         1 ]
%     %||   .quality
%     %||     [ 0.685     0.742     0.617     0.648     0.757     0.591 ]
%     %||   .nchunks
%     %||     [ 5         5         5         5         5         5 ]
%     %|| .a
%     %||   .fdim
%     %||     .labels
%     %||       { 'i'  'j'  'k' }
%     %||     .values
%     %||       { [ 1         2         3 ]  [ 1         2 ]  [ 1 ] }
%     %||   .sdim
%     %||     .labels
%     %||       { 'targets1'  'targets2' }
%     %||     .values
%     %||       { [ 1    [ 1
%     %||           2      2
%     %||           3      3
%     %||           4 ]    4 ] }
%     %||   .vol
%     %||     .mat
%     %||       [ 2         0         0        -3
%     %||         0         2         0        -3
%     %||         0         0         2        -3
%     %||         0         0         0         1 ]
%     %||     .dim
%     %||       [ 3         2         1 ]
%     %||     .xform
%     %||       'scanner_anat'
%     %|| .sa
%     %||   .targets1
%     %||     [ 1
%     %||       2
%     %||       3
%     %||       :
%     %||       2
%     %||       3
%     %||       4 ]@16x1
%     %||   .targets2
%     %||     [ 1
%     %||       1
%     %||       1
%     %||       :
%     %||       4
%     %||       4
%     %||       4 ]@16x1
%
% Reference:
%   - Abdi, H., Valentin, D., O?Toole, A. J., & Edelman, B. (2005).
%     DISTATIS: The analysis of multiple distance matrices. In
%     Proceedings of the IEEE Computer Society: International conference
%     on computer vision and pattern recognition, San Diego, CA, USA
%     (pp. 42?47).
%
% Notes:
%   - DISTATIS tries to find an optimal compromise distance matrix across
%     the different samples (participants)
%   - Output can be reshape to matrix or array form using
%     cosmo_unflatten(res,1)
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    cosmo_check_external('distatis');

    defaults.return='distance';
    defaults.split_by='chunks';
    defaults.shape='square';
    defaults.mask_output=[];
    defaults.progress=100;
    defaults.feature_ids=[];
    defaults.autoscale=true;
    defaults.abs_correlation=false;
    defaults.weights='eig';

    opt=cosmo_structjoin(defaults,varargin);

    subject_cell=get_subject_data(ds,opt);
    nsubj=numel(subject_cell);

    [dsms,nclasses,dim_labels,dim_values]=get_dsms(subject_cell);

    feature_ids=get_feature_ids(size(dsms{1},3),opt);
    nfeatures=numel(feature_ids);

    quality=zeros(1,nfeatures);
    nobservations=zeros(1,nfeatures);

    prev_msg='';
    clock_start=clock();
    show_progress=nfeatures>1 && opt.progress;

    for k=1:nfeatures
        feature_id=feature_ids(k);
        x=zeros(nclasses*nclasses,nsubj);

        for j=1:nsubj
            dsm=dsms{j}(:,:,feature_id);
            x(:,j)=distance2crossproduct(dsm, opt.autoscale);
        end

        [x,subj_msk]=cosmo_remove_useless_data(x);
        nkeep=sum(subj_msk);

        % equivalent, but slower:
        % [e,v]=eigs(c,1);

        [ew,v]=get_weights(x, feature_id, nkeep, opt);

        % compute compromise
        compromise=x*ew;

        result=convert_compromise(compromise, opt);

        if feature_id==1
            % allocate space
            samples=zeros(numel(result),nfeatures);
        end

        samples(:,k)=result;

        quality(:,k)=v/nkeep;
        nobservations(:,k)=nkeep;

        if show_progress && (k<10 || ...
                                mod(k, opt.progress)==0 || ...
                                k==nfeatures)
            status=sprintf('quality=%.3f%% (avg)',mean(quality(1:k)));
            prev_msg=cosmo_show_progress(clock_start,k/nfeatures,...
                                                        status,prev_msg);
        end
    end

    % set output in either triangular or square shape
    [res,i,j]=get_samples_in_shape(samples,nclasses,opt.shape);
    res=copy_fields(ds,res,{'fa','a'});

    % add attributes
    res.fa.quality=quality;
    res.fa.nchunks=nobservations;
    res.a.sdim=struct();
    res.a.sdim.labels=dim_labels;
    res.a.sdim.values=dim_values;

    res.sa.(dim_labels{1})=i(:);
    res.sa.(dim_labels{2})=j(:);

    cosmo_check_dataset(res);

function [res,i,j]=get_samples_in_shape(samples,nclasses,shape)
    res=struct();
    switch shape
        case 'triangle'
            [msk,i,j]=distance_matrix_mask(nclasses);
            res.samples=samples(msk(:),:);
        case 'square'
            res.samples=samples;
            [i,j]=find(ones(nclasses));
        otherwise
            error('unsupported direction %s', shape);
    end

function dst=copy_fields(src,dst,keys)
    for k=1:numel(keys)
        key=keys{k};
        if isfield(src,key)
            dst.(key)=src.(key);
        end
    end

function feature_ids=get_feature_ids(nfeatures, opt)
    feature_ids=opt.feature_ids;
    if isempty(feature_ids);
        feature_ids=1:nfeatures;
    end

function [ew,v]=get_weights(x, feature_id, nkeep, opt)
    switch opt.weights
        case 'eig'
            [ew,v]=eigen_weights(x, feature_id);

        case 'uniform'
            % all the same (allowing for comparison with 'eig')
            ew=ones(nkeep,1)/nkeep;
            v=0;

        otherwise
            error('illegal weight %s', opt.weights);
    end

function subject_cell=get_subject_data(ds,opt)
    if isstruct(ds)
        subject_cell=cosmo_split(ds,opt.split_by);
    else
        subject_cell=ds;
    end

    if numel(subject_cell)==0
        error('empty input');
    end

function [ew,v]=eigen_weights(x, feature_id)

    c=cosmo_corr(x);

    negative_c=c<0;

    if any(negative_c(:))

        [i,j]=find(negative_c);
        error(['feature %d has negative correlation between '...
                'sample %d and %d, which is not supported by '...
                'distatis. DISTATIS assumes that the similarity '...
                'data from all samples (typically: participants) '...
                'correlate positively. Because that is not the '...
                'case, you cannot use DISTATIS analysis on this '...
                'data. '],...
                feature_id,i(1),j(1));
    end

    [v,e]=fast_eig1(c);

    if all(e<0)
        e=-e;
    end

    assert(all(e>0));
    assert(v>0);

    % normalize first eigenvector
    ew=e/sum(e);

function result=convert_compromise(compromise, opt)
    switch opt.return
        case 'crossproduct'
            result=compromise;
        case 'distance'
            result=crossproduct2distance(compromise);
        otherwise
            error('illegal opt.return');
    end

function z=crossproduct2distance(x)
    n=sqrt(numel(x));
    e=ones(n,1);
    d=x(1:(n+1):end);
    dd=d*e';
    ddt=dd';
    y=dd(:)+ddt(:)-2*x;
    z=ensure_distance_vector(y);

function assert_symmetric(x, tolerance)
    if nargin<2, tolerance=1e-8; end

    % assert x is a square matrix
    sz=size(x);
    assert(isequal(sz,sz([2 1])));

    xx=x'-x;

    msk=xx>tolerance;
    if any(msk)
        [i,j]=find(msk,1);
        error('not symmetric: x(%d,%d)=%d ~= %d=x(%d,%d)',...
                i,j,x(i,j),x(j,i),j,i);
    end

function z_vec=distance2crossproduct(x, autoscale)

    n=size(x,1);
    e=ones(n,1);
    m=e*(1/n);
    ee=eye(n)-e*m';
    y=-.5*ee*(x+x')*ee';
    if autoscale
        z=(1/fast_eig1(y))*y;
    else
        z=y;
    end
    assert_symmetric(z);
    % equivalent, but slower:
    % z=(1/eigs(y,1))*y(:);

    z_vec=z(:);

function [lambda,pivot]=fast_eig1(x)
    % returns the first eigenvalue in lambda, and the corresponding
    % eigenvector in pivot
    if cosmo_wtf('is_matlab')
        [pivot,lambda]=eigs(x,1);
    else
        % There seems a bug in Octave for 'eigs',
        % so use 'eig' instead.
        % http://savannah.gnu.org/bugs/?44004
        [e,v]=eig(x);
        diag_v=diag(v);

        % find largest eigenvalue and eigenvector
        [lambda,i]=max(diag_v);
        pivot=e(:,i);
    end

    % The code below is disabled because under certain circumstances
    % it would return a near-zero eigenvalue if indeed one eigenvalue (but
    % not the largest one) is zero.
    % % compute first (largest) eigenvalue and corresponding eigenvector
    % % using power iteration method; benchmarking suggests this can be up to
    % % five times as fast as using eigs(x,1)
    % n=size(x,1);
    % pivot=ones(n,1);
    % tolerance=1e-8;
    % max_iter=1000;
    %
    % old_lambda=NaN;
    % for k=1:max_iter
    %     z=x*pivot;
    %     pivot=z / norm(z);
    %
    %     lambda=pivot'*z;
    %     if abs(lambda-old_lambda)/lambda<tolerance
    %         z=x*pivot;
    %         pivot=z / sqrt(sum(z.^2));
    %
    %         lambda=pivot'*z;
    %         return
    %     end
    %     old_lambda=lambda;
    % end
    %
    % % matlab fallback
    % [pivot,lambda]=eigs(x,1);

function y=ensure_distance_vector(x)
    tolerance=1e-8;

    n=sqrt(numel(x));
    xsq=reshape(x,n,n);

    dx=diag(xsq);
    assert(all(dx<tolerance));

    xsq=xsq-diag(dx);

    delta=xsq-xsq';
    assert(all(delta(:)<tolerance));

    xsq=.5*(xsq+xsq');
    y=xsq(:);

function [dsms,nclasses,dim_labels,dim_values]=get_dsms(data_cell)
    nsubj=numel(data_cell);

    % allocate
    dsms=cell(nsubj,1);
    for k=1:numel(data_cell)
        data=data_cell{k};

        % get data
        [dsm,dim_labels,dim_values,is_ds]=get_dsm(data);

        % store data
        dsms{k}=dsm;

        if k==1
            nclasses=size(dsm,1);
            first_dim_labels=dim_labels;
            first_dim_values=dim_values;

            data_first=data;
        else

            if ~isequal(first_dim_labels,dim_labels)
                error('dim label mismatch between subject 1 and %d',k);
            end
            if ~isequal(first_dim_values,dim_values)
                error('dim label mismatch between subject 1 and %d',k);
            end

            % check for compatibility over subjects, raises an error if not
            % kosher
            if is_ds
                cosmo_stack({cosmo_slice(data,1),...
                                cosmo_slice(data_first,1)},1,'unique');
            end
        end
    end

function [msk,i,j]=distance_matrix_mask(nclasses)
    msk=triu(repmat(1:nclasses,nclasses,1),1)'>0;
    [i,j]=find(msk);

function [dsm, dim_labels, dim_values, is_ds]=get_dsm(data)
    is_ds=isstruct(data);
    if is_ds
        [dsm,dim_labels,dim_values]=cosmo_unflatten(data,1);
    elseif isnumeric(data)
        sz=size(data);
        if numel(sz)~=2
            error('only vectorized distance matrices are supported');
        end
        [n,nfeatures]=size(data);

        side=(1+sqrt(1+8*n))/2; % so that side*(side-1)/2==n
        if ~isequal(side, round(side))
            error(['size %d of input vector is not correct for '...
                    'the number of elements below the diagonal of a '...
                    'square (distance) matrix'], n);
        end

        [msk,i,j]=distance_matrix_mask(side);
        dsm=zeros([side,side,nfeatures]);

        assert(numel(i)==n);
        for pos=1:n
            dsm(i(pos),j(pos),:)=data(pos,:);
        end

        sq1=cosmo_squareform(data(:,1));
        dsm1=dsm(:,:,1);
        assert(isequal(sq1,dsm1+dsm1'));

        dim_labels={'targets1','targets2'};
        dim_values={(1:side)',(1:side)'};
    else
        error('illegal input: expect dataset struct, or cell with arrays');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function hdr = mayo2nifti1(ohdr,mat)
% Convert from an ANALYZE to a NIFTI-1 header
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

if isfield(ohdr,'magic'),
    hdr = ohdr;
    return;
end;
hdr            = empty_hdr;
hdr.dim        = ohdr.dim;
hdr.datatype   = ohdr.datatype;
hdr.bitpix     = ohdr.bitpix;
hdr.pixdim     = ohdr.pixdim;
hdr.vox_offset = ohdr.vox_offset;
hdr.scl_slope  = ohdr.roi_scale;
hdr.scl_inter  = ohdr.funused1;
hdr.descrip    = ohdr.descrip;
hdr.aux_file   = ohdr.aux_file;
hdr.glmax      = ohdr.glmax;
hdr.glmin      = ohdr.glmin;
hdr.cal_max    = ohdr.cal_max;
hdr.cal_min    = ohdr.cal_min;
hdr.magic      = 'ni1';

switch hdr.datatype,
case 130, hdr.datatype = 256; %  int8
case 132, hdr.datatype = 512; % uint16
case 136, hdr.datatype = 768; % uint32
end;

if nargin<2,
    % No mat, so create the equivalent from the hdr...
    if any(ohdr.origin(1:3)), origin = double(ohdr.origin(1:3));
    else                      origin = (double(ohdr.dim(2:4))+1)/2; end;
    vox    = double(ohdr.pixdim(2:4));
    if vox(1)<0,
        % Assume FSL orientation
        flp    = 0;
    else
        % Assume SPM or proper Analyze
        flp    = spm_flip_analyze_images;
    end;
    if all(vox == 0), vox = [1 1 1]; end;
    off    = -vox.*origin;
    mat    = [vox(1) 0 0 off(1) ; 0 vox(2) 0 off(2) ; 0 0 vox(3) off(3) ; 0 0 0 1];
    if flp,
        %disp(['Assuming that image is stored left-handed']);
        mat = diag([-1 1 1 1])*mat;
    else
        %disp(['Assuming that image is stored right-handed']);
    end;
end;

hdr            = encode_qform0(mat,hdr);
mat            = mat*[eye(4,3) [1 1 1 1]'];
hdr.srow_x     = mat(1,:);
hdr.srow_y     = mat(2,:);
hdr.srow_z     = mat(3,:);
hdr.qform_code = 2;
hdr.sform_code = 2;
return;

author: fieldtrip

repo: fieldtrip