text:
function signal = set_selinterval(varargin)
% Selects a time interval from a data set.
% Signal = set_selinterval(Signal,Intervals,Unit)
%
% In:
%   Signal    :   Continuous EEGLAB data set from which an interval should be selected
%
%   Intervals :   Selection interval(s) formatted as [start end; start end; start end; ...] or
%                 a MATLAB index range for the samples (if IntervalUnit is set to 'range')
%
%   IntervalUnit : Unit of measurement for the interval. Either 'seconds', 'samples', 'fraction',
%                  or 'range' (default: seconds)
%
%   InsertBoundaryMarkers : Whether to insert boundary markers (default: false)
%
% Out:
%   Signal    :   data set restricted to the selected range. The following fields are updated:
%                 * .data and all other time series fields are reduced to the selected intervals
%                 * .event is reduced to selected events
%                 * .xmax/.pnts are updated
%                 * optionally boundary events are inserted
%
% Examples:
%   % for a continuous data set, retain only the data within 50s and 200s, as well as 1200s and 1500s
%   eeg = set_selinterval(eeg,[50 200; 1200 1500])
%
%   % for a continuous data set, retain only the data within the 1000's and the 10000's sample
%   eeg = set_selinterval(eeg,[1000 10000],'samples')
%
%   % for a continuous data set, retain only the last half of the data
%   eeg = set_selinterval(eeg,[0.5 1],'fraction')
%
%   % as before, but pass the arguments by name
%   eeg = set_selinterval('Signal',eeg, 'Intervals',[0.5 1], 'Unit','fraction')
%
% See also:
%   set_selepos, set_partition
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-01
dp;

% set_selinterval_version<2.0> -- for the cache

if ~exp_beginfun('editing') return; end

declare_properties('independent_channels',true,'independent_trials',false);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg({'intervals','Intervals'},[],[], 'Selection intervals. Array of the form [start end; start end; start end; ...].'), ...
    arg({'selunit','IntervalUnit','Unit','unit'},'seconds',{'seconds','samples','fraction','range'}, 'Interval unit. The unit of measurement for selection intervals (does not apply to SampleRange).'), ...
    arg({'insert_boundary_markers','InsertBoundaryMarkers'},false,[], 'Insert boundary markers. Whether to insert boundary markers (for EEGLAB compatibility).'));

utl_check_fields(signal,{'epoch','event','data','pnts','srate','xmin'},'signal','signal');
if signal.srate == 0
    error('Your signal needs to have a nonzero .srate value.'); end
if ~isempty(signal.epoch) || size(signal.data,3) > 1
    error('This function cannot be applied to epoched data. Use flt_window to perform interval selection on epoched data.'); end

switch selunit
    case 'fraction'
        if any(intervals<0 | intervals>1) %#ok<*NODEF>
            error('Fractional intervals must be in the [0,1] range.'); end
        intervals = 1+round(intervals*(size(signal.data,2)-1));
        samplerange = [];
    case 'seconds'
        intervals(intervals==Inf) = (size(signal.data,2)-1)/signal.srate;
        if any(intervals<0 | (1+round(intervals*signal.srate)) > size(signal.data,2))
            error('Your intervals exceed the data limit.'); end
        intervals = min(1+round(intervals*signal.srate),size(signal.data,2));
        samplerange = [];
    case 'range'
        samplerange = intervals(:)';
        intervals = [];
    case 'samples'
        samplerange = [];
    otherwise
        error(['Unsupported selection unit: ' selunit]);
end

% generate samplerange from intervals if not yet present
if isempty(samplerange)
    samplerange = cell(1,size(intervals,1));
    for k=1:length(samplerange)
        samplerange{k} = intervals(k,1):intervals(k,2); end
    samplerange = [samplerange{:}];
elseif isa(samplerange,'logical')
    samplerange = find(samplerange);
end
if any(round(samplerange) ~= samplerange)
    error('Your sample range contains non-integer indices: %s',hlp_tostring(samplerange)); end
if any(samplerange<1 | samplerange>size(signal.data,2))
    error('Your sample range exceeds the data bounds (%i); range is: %s.',size(signal.data,2),hlp_tostring(samplerange)); end

if isempty(samplerange)
    % clear data
    for field = utl_timeseries_fields(signal)
        signal.(field{1}) = []; end
    signal.event = [];
    signal.pnts = 0;
elseif ~isequal(samplerange,1:size(signal.data,2))
    % select range within the time series fields
    for field = utl_timeseries_fields(signal)
        if ~isempty(signal.(field{1}))
            try
                signal.(field{1}) = signal.(field{1})(:,samplerange,:,:,:,:,:,:);
            catch e
                error('The given sample range could not be extracted from the time-series field .%s with error: %s. Field size was: %i, epoch indices were: %s',field{1},e.message,size(signal.(field{1}),2),hlp_tostring(samplerange));
            end
        end
    end
    % select range within the events
    if ~isempty(signal.event)
        if ~isfield(signal.event,'latency')
            error('Your signal is lacking the required .event.latency field.'); end
        latency_numels = cellfun('prodofsize',{signal.event.latency});
        if any(latency_numels == 0)
            error('The given signal has one or more events with empty .latency field. This is not permitted.');
        elseif any(latency_numels ~= 1)
            error('The given signal has one or more events with a .latency value that is not a scalar. This is not permitted.');
        end
        % trim out-of-bounds events
        event_latency = [signal.event.latency];
        out_of_bounds = event_latency<1 | event_latency>signal.pnts;
        if any(out_of_bounds)
            disp_once('WARNING: your signal had %i events at out-of-bounds latencies that were removed.',nnz(out_of_bounds)); end
        signal.event(out_of_bounds) = [];
        % bin latencies into sparse array and apply selection
        [event_positions,residuals] = sparse_binning([signal.event.latency],[],signal.pnts);
        [ranks,sample_indices,event_indices] = find(event_positions(:,samplerange));
        residuals = residuals(:); sample_indices = sample_indices(:); event_indices = event_indices(:);
        % write event subset and update latencies
        if ~isempty(ranks)
            signal.event = signal.event(event_indices);
            [signal.event.latency] = arraydeal(sample_indices+residuals(event_indices));
        else
            signal.event = [];
        end
    end
    % optionally insert boundary markers
    if insert_boundary_markers
        % calculate gaps between intervals
        if isempty(intervals)
            tmp = diff(samplerange)-1;
            gap_offsets = find(tmp);
            gap_lengths = tmp(gap_offsets);
        else
            gap_offsets = intervals(1:end-1,2);
            gap_lengths = intervals(2:end,1) - intervals(1:end-1,2) - 1;
        end
        % append new events for each boundary between intervals
        insert_range = length(signal.event)+(1:length(gap_offsets));
        if insert_range
            [signal.event(insert_range).type] = deal('boundary');
            [signal.event(insert_range).latency] = arraydeal(gap_offsets+0.5);
            [signal.event(insert_range).duration] = arraydeal(gap_lengths);
            % resort events
            [dummy,order] = sort([signal.event.latency]); %#ok<ASGLU>
            signal.event = signal.event(order);
        end
    end
    % update misc fields
    signal.pnts = size(signal.data,2);
    signal.xmax = signal.xmin + (max(samplerange)-1)/signal.srate;
    signal.xmin = signal.xmin + (min(samplerange)-1)/signal.srate;
end

exp_endfun;

repo: REST


text:
classdef GaussianMultC
%%GAUSSIANMULTC A collection of functions for the Gaussian multivariate
%               copulae. Implemented functions are: PDF.
%
%REFERENCES:
%[1] H. Joe and D. Kurowicka, Dependence Modeling: Vine Copula Handbook.
%    World Scientific, 2011.
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    methods(Static)
        function jprob = PDF(u,R)
        %PDF Generates a joint probability for column vector u or a row
        %    vector of joint probabilities for matrix u with each column as
        %    the input vector.
        %
        %INPUTS: u: A d-by-n matrix where each column represents a
        %           d-dimensional random vector and n is the number of
        %           random vectors. Entries must be in [0,1].
        %        R: A d-by-d correlation matrix.
        %
        %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
        %
        if ~exist('u','var') || isempty(u)
            jprob = [];
            return
        end

        s2 = size(u,2);
        I = eye(size(R));
        jprob = zeros(1,s2);

        for col = 1:s2
            v = GaussianD.invCDF(u(:,col));
            jprob(col) = exp(-0.5*v'*(inv(R)-I)*v)/sqrt(det(R));
        end
        end
    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 vers = eegplugin_BCI2000import(fig, trystrs, catchstrs)
% (C) 2000-2009, BCI2000 Project
% http://www.bci2000.org

    vers = 'BCI2000import1.00';
    if nargin < 3
        error('eegplugin_BCI2000import requires 3 arguments');
    end;

    % add folder to path
    % ------------------
    if ~exist('BCI2000import')
        p = which('eegplugin_BCI2000import.m');
        p = p(1:findstr(p,'eegplugin_BCI2000import.m')-1);
        addpath([ p 'BCI2000import1.00' ] );
    end;

    % find import data menu
    % ---------------------
    menu = findobj(fig, 'tag', 'import data');

    % menu callbacks
    % --------------
    comcnt = [ trystrs.no_check '[EEG] = BCI2000import;'     catchstrs.new_and_hist ];

    % create menus
    % ------------
    uimenu( menu, 'label', 'From BCI2000 .DAT file', 'callback', comcnt, 'separator', 'on' );

author: PatternRecognition

repo: OpenBMI


text:
function [lb,ub,redundant,psstruct,infeasible] = milppresolve(A,b,lb,ub,integer_variables,binary_variables,changed_bounds);
%MILPPRESOLVE Internal function for presolving MILPs

% Simple bound pre-processing (paper by Savelsbergh)
% No code optimization at all

if nargin < 5
    integer_variables = [];
end
if nargin < 6
    binary_variables = [];
end
if nargin < 7
    changed_bounds = [];
end

ubnew = ub;
lbnew = lb;
goon = 1;

AL0A  = (A>0).*A;
AG0A  = (A<0).*A;

At = A';
use_indicies=ones(length(b),1);
used = full(any(A(:,find(changed_bounds)),2));
isbinary  = ismembc(1:length(lb),binary_variables);
isinteger = ismembc(1:length(lb),binary_variables);

goon = all(lb<=ub);
infeasible = ~goon;
if ~infeasible

    bi_up = AL0A*ub+AG0A*lb;
    bi_dn = AL0A*lb+AG0A*ub;

    while goon

        bminusbdn=b-bi_dn;
        for i = find(use_indicies & used)'

            [ii,jj,kk] = find(At(:,i));

            Cp = ii(kk>0);
            if ~isempty(Cp)
                new1=lbnew(Cp)+bminusbdn(i)./At(Cp,i);
                ubnew(Cp)=min(ubnew(Cp),new1);
            end

            Cm = ii(kk<0);
            if ~isempty(Cm)
                new2=ubnew(Cm)+bminusbdn(i)./At(Cm,i);
                lbnew(Cm)=max(lbnew(Cm),new2);
            end

            if any(lbnew>ubnew)
                infeasible = 1;
                break
            end

        end

        if ~infeasible

            lbnew(integer_variables) = round(lbnew(integer_variables)+0.49999);
            ubnew(integer_variables) = round(ubnew(integer_variables)-0.49999);

            lbnew(binary_variables) = round(lbnew(binary_variables)+0.49999);
            ubnew(binary_variables) = round(ubnew(binary_variables)-0.49999);

            goon = (~all((lb==lbnew) & (ub==ubnew))) & all(lbnew<=ubnew);

            used = (lb~=lbnew) | (ub~=ubnew);
            used = full(any(A(:,find(used)),2));

            lb = lbnew;
            ub = ubnew;

            bi_up = AL0A*ub+AG0A*lb;
            bi_dn = AL0A*lb+AG0A*ub;
            redundant = find(bi_up<=b);

            use_indicies = use_indicies & bi_up>b;
        else
            goon = 0;
        end
    end
    redundant = find(bi_up<=b);
else
    redundant=[];
end
psstruct.AL0A = AL0A;
psstruct.AG0A = AG0A;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Section 11.8.4: Network rate optimization
% Boyd & Vandenberghe "Convex Optimization"
% Argyrios Zymnis - 05/03/08
%
% We consider a network with n flows and L links. Each flow i,
% moves along a fixed predetermined path (i.e. a subset of the links)
% and has an associated rate x_i. Each link j has an associated capacity
% c_j. The total rate of all flows travelling along a link cannot exceed
% the link capacity. We can describe these link capacity limits using the
% flow-link incidence matrix A \in \reals^{L \times n}, where
% A_{ij} = 1, if flow j passes through link i and 0 otherwise.
% The link capacity constraints can be expressed as A*x <= c
% In the network rate problem the variables are the flow rates x. The
% objective is to choose the flow rates to maximize a separate utility
% function U, given by
%           U(x) = U_1(x_1)+U_2(x_2)+...+U_n(x_n)
% The network rate optimization problem is then
%           maximize    U(x)
%           subject to  A*x <= c
% Here we use U_i(x_i) = log x_i for all i

% Input data
rand('state',1)
L = 20;
n = 10;
k = 7; %average links per flow
A = double(rand(L,n) <= k/L);
c = 0.9*rand(L,1)+0.1;

% Solve network rate problem
cvx_begin
    variable x(n);
    maximize(sum(log(x)))
    subject to
        A*x <= c;  %#ok
cvx_end
primal_obj = cvx_optval;

% Solve dual problem to obtain link prices
cvx_begin
    variable lambda(L);
    minimize(c'*lambda-sum(log(A'*lambda))-n)
    subject to
        lambda >= 0;  %#ok
cvx_end
dual_obj = cvx_optval;

author: yu-jiang

repo: radpbook

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


text:
function varargout = setInputName(obj,varargin)
    nVarargs = length(varargin);
    name = varargin{1};
    if nVarargs == 2
        index = varargin{2};
    else
        index = 1:obj.dim.u;
    end
    varargout = cell(size(index));
    j = 1;
    for i=index
        obj.u(i) = sym(name{j});
        varargout{j} = sym(name{j});
        j = j+1;
    end
end

author: deng-haoyang

repo: ParNMPC


text:
%type3detransform.m
%written by: Duncan Po
%Date: June 27/2002
% utility file used in tree2contourlet.m
% usage [x,y] = type3detransform(z)
% decomposes z into x and y

function [x,y] = type3detransform(z)

row = size(z,1);
col = size(z,2);
z = z.';
z = reshape(z, col*2, row/2);
z = z.';
x = z(:, 1:col);
y = z(:,(col+1):(2*col));

author: Sable

repo: mcbench-benchmarks


text:
function value = i4_mach ( i )

%*****************************************************************************80
%
%% I4_MACH returns integer machine constants.
%
%  Discussion:
%
%    Input/output unit numbers.
%
%      I1MACH(1) = the standard input unit.
%      I1MACH(2) = the standard output unit.
%      I1MACH(3) = the standard punch unit.
%      I1MACH(4) = the standard error message unit.
%
%    Words.
%
%      I1MACH(5) = the number of bits per integer storage unit.
%      I1MACH(6) = the number of characters per integer storage unit.
%
%    Integers.
%
%    Assume integers are represented in the S digit base A form:
%
%      Sign * (X(S-1)*A**(S-1) + ... + X(1)*A + X(0))
%
%    where 0 <= X(1:S-1) < A.
%
%      I1MACH(7) = A, the base.
%      I1MACH(8) = S, the number of base A digits.
%      I1MACH(9) = A**S-1, the largest integer.
%
%    Floating point numbers
%
%    Assume floating point numbers are represented in the T digit
%    base B form:
%
%      Sign * (B**E) * ((X(1)/B) + ... + (X(T)/B**T) )
%
%    where 0 <= X(I) < B for I=1 to T, 0 < X(1) and EMIN <= E <= EMAX.
%
%      I1MACH(10) = B, the base.
%
%    Single precision
%
%      I1MACH(11) = T, the number of base B digits.
%      I1MACH(12) = EMIN, the smallest exponent E.
%      I1MACH(13) = EMAX, the largest exponent E.
%
%    Double precision
%
%      I1MACH(14) = T, the number of base B digits.
%      I1MACH(15) = EMIN, the smallest exponent E.
%      I1MACH(16) = EMAX, the largest exponent E.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2007
%
%  Author:
%
%    Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Phyllis Fox, Andrew Hall, Norman Schryer,
%    Algorithm 528,
%    Framework for a Portable Library,
%    ACM Transactions on Mathematical Software,
%    Volume 4, Number 2, June 1978, page 176-188.
%
%  Parameters:
%
%    Input, integer I, chooses the parameter to be returned.
%    1 <= I <= 16.
%
%    Output, integer VALUE, the value of the chosen parameter.
%
  if ( i < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 16.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0;
    error ( 'I4_MACH - Fatal error!' );
  elseif ( i == 1 )
    value = 5;
  elseif ( i == 2 )
    value = 6;
  elseif ( i == 3 )
    value = 7;
  elseif ( i == 4 )
    value = 6;
  elseif ( i == 5 )
    value = 32;
  elseif ( i == 6 )
    value = 4;
  elseif ( i == 7 )
    value = 2;
  elseif ( i == 8 )
    value = 31;
  elseif ( i == 9 )
    value = 2147483647;
  elseif ( i == 10 )
    value = 2;
  elseif ( i == 11 )
    value = 24;
  elseif ( i == 12 )
    value = -125;
  elseif ( i == 13 )
    value = 128;
  elseif ( i == 14 )
    value = 53;
  elseif ( i == 15 )
    value = -1021;
  elseif ( i == 16 )
    value = 1024;
  elseif ( 16 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 16.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0;
    error ( 'I4_MACH - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fimu, imudata, previousImuData]=readimuheader(imufile, ...
    previousImuData, startTime, numPrevImuDataToKeep, imuFileType)
% input
% imuFileType
%   0: plain text 3dm gx3-35 data,
%   1: H764G 1 comprehensive csv
%   2: steval mki 062v2 from Yujia file
%   3: steval iNemo suite output data
%   4: microstrain csv
%   5: epson csv

% imufile: imu file name

% To load IMU data from the very beginning, set startTime <= 0.

% output:
% fimu the file pointer after opening the file
% imudata: the entry of imu data that has timestamp no less than startTime
% of format [time(sec), ax, ay, az(,/s^2), wx, wy, wz(rad/s)]
% previousImuData: stores numPrevImuDataToKeep inertial data that are just
% before the retrieved imu data

if(nargin<5)
    imuFileType=0;
end
switch(imuFileType)
    case 0% plain text 3dm gx3-35 data
        % make sure previousImuData is transferred correctly, no memory leaking
        fimu=fopen(imufile,'r');
        fgetl(fimu);% remove the header
        hstream= fgetl(fimu);
        mass=textscan(hstream,'%f','delimiter',',');
        imudata=mass{1};
        while(imudata(7,1)==0||imudata(1,1)<startTime)
            if(length(previousImuData)==numPrevImuDataToKeep)
                previousImuData(1) = [];
            end
            previousImuData{end + 1} = imudata(:,end);
            hstream= fgetl(fimu);
            mass=textscan(hstream,'%f','delimiter',',');
            imudata=mass{1};
        end
    case 1 % H764G-1 comprehensive csv data
        %open files
        fimu=fopen(imufile,'r');
        h= fgetl(fimu);
        while(true)
            if(contains(h,'RAW 01'))
                break;
            else h= fgetl(fimu);
            end
        end
        h= fgetl(fimu);
        mass=textscan(h,'%f','delimiter',',');
        imudata=mass{1};
        while(imudata(9,1)==0||imudata(2,1)<startTime)
            imudata=imudata(2:8,1);% gps time, xyz delta v, xyz delta theta
            imudata(2:4,1)=imudata(2:4,1)*.3048;% convert to metric unit meter
            if(length(previousImuData)==numPrevImuDataToKeep)
                previousImuData(1) = [];
            end
            previousImuData{end + 1} = imudata(:,end);
            h= fgetl(fimu);
            mass=textscan(h,'%f','delimiter',',');
            imudata=mass{1};
        end
        imudata=imudata(2:8,1);% gps time, xyz delta v, xyz delta theta
        imudata(2:4,1)=imudata(2:4,1)*.3048;% convert to metric unit meter
    case 2 % steval mki 062v2 from Yujia file
        % make sure previousImuData is transferred correctly, no memory leaking
        fimu=fopen(imufile,'r');
        hstream= fgetl(fimu);
        if(isempty(hstream))
            hstream= fgetl(fimu);
        end
        hstream=fgetl(fimu);% discard first two lines
        hstream=fgetl(fimu);
        mass=textscan(hstream,'%f','delimiter',' ');
        imudata=mass{1};
        hstream= fgetl(fimu);
        imu_scalefactor = 9.78;
        while(imudata(1,1)<startTime)
            imudata(2:4,1) = imudata(2:4,1)/1000.0 * imu_scalefactor;
            imudata(5:7,1) = imudata(5:7,1)*pi/180;
            imudata(6,1) = -imudata(6,1); % wrong data format
            if(length(previousImuData)==numPrevImuDataToKeep)
                previousImuData(1) = [];
            end
            previousImuData{end + 1} = imudata(:,end);
            hstream= fgetl(fimu);
            mass=textscan(hstream,'%f','delimiter',' ');
            imudata=mass{1};
            hstream= fgetl(fimu);
        end
        imudata(2:4,1) = imudata(2:4,1)/1000.0 * imu_scalefactor;
        imudata(5:7,1) = imudata(5:7,1)*pi/180;
        imudata(6,1) = -imudata(6,1);
    case 3 % steval iNemo suite output data
        fimu=fopen(imufile,'r');
        fgetl(fimu); % discard first line
        hstream=fgetl(fimu);
        mass=textscan(hstream,'%f','delimiter',',');
        imudata=mass{1};
        imudata(1)=imudata(1)/1000;
        imu_scalefactor = 9.8;
        while(imudata(1,1)<startTime)
            imudata(2:4,1) = imudata(2:4,1)/1000.0 * imu_scalefactor;
            imudata(5:7,1) = imudata(5:7,1)*pi/180;
            imudata(6,1) = -imudata(6,1); % wrong data format
            if(length(previousImuData)==numPrevImuDataToKeep)
                previousImuData(1) = [];
            end
            previousImuData{end + 1} = imudata(:,end);
            hstream= fgetl(fimu);
            mass=textscan(hstream,'%f','delimiter',',');
            imudata=mass{1};
            imudata(1)=imudata(1)/1000;
        end
        imudata(2:4,1) = imudata(2:4,1)/1000.0 * imu_scalefactor;
        imudata(5:7,1) = imudata(5:7,1)*pi/180;
        imudata(6,1) = -imudata(6,1);
    case 4 % microstrain csv data
        fimu=fopen(imufile,'r');
        % remove first 16 lines
        headerLines=16;
        for i=1:headerLines
            fgetl(fimu); % discard first lines
        end
        hstream=fgetl(fimu);
        mass=textscan(hstream,'%f','delimiter',',');
        imudata=mass{1};
        imudata=imudata([3,8:13]);

        while(imudata(1)<startTime || sum(isnan(imudata)))
            if(~sum(isnan(imudata)))
                imudata(2:4)=imudata(2:4)*9.80665;
                if(length(previousImuData)==numPrevImuDataToKeep)
                    previousImuData(1) = [];
                end
                previousImuData{end + 1} = imudata(:,end);
            end
            hstream=fgetl(fimu);
            if (~ischar(hstream))
                imudata=[];
                return;
            end
            mass=textscan(hstream,'%f','delimiter',',');
            imudata=mass{1};
            imudata=imudata([3,8:13]);
        end
        imudata(2:4)=imudata(2:4)*9.80665;
    case 5 % read the header of epson csv data
        fimu=fopen(imufile,'r');
        % remove first 4 lines
        headerLines=4;
        for i=1:headerLines
            fgetl(fimu); % discard header
        end
        hstream=fgetl(fimu);
        mass=textscan(hstream,'%f','delimiter',',');
        imudata=mass{1};
        imudata=imudata([2,[8:10, 5:7]]);
        while(imudata(1)<startTime)
            imudata(2:4)=imudata(2:4)*9.8/1000;
            imudata(5:7)=imudata(5:7)*pi/180;
            if(length(previousImuData)==numPrevImuDataToKeep)
                previousImuData(1) = [];
            end
            previousImuData{end + 1} = imudata(:,end);
            hstream=fgetl(fimu);
            if (~ischar(hstream))
                imudata=[];
                return;
            end
            mass=textscan(hstream,'%f','delimiter',',');
            imudata=mass{1};
            imudata=imudata([2,[8:10, 5:7]]);
        end
        imudata(2:4)=imudata(2:4)*9.80665/1000;
        imudata(5:7)=imudata(5:7)*pi/180;
    otherwise
        fprintf('IMU file type %d is unsupported!\n', imuFileType);
        imudata = [];
end

repo: ekfmonoslam


text:
% Demo for phase-space reduction (cf. Centre-Manifold theory): In this
% example we show that Generalised filtering can be used to approximate the
% dynamics on the Centre-Manifold of a chaotic system. This example uses a
% nominal polynomial for the reduced dynamics (cf, order parameters) and
% exploits parameter optimisation to learn their coefficients. The approach
% appeals to Centre Manifold Theory and the enslaving principle by assuming
% that stable modes collapse sufficiently quickly to be treated as random
% fluctuations. This means that the unstable (slow modes or patterns)
% dynamics are sufficient to characterise the emergent order (e.g. pattern
% formation) and can be accessed using a generative model that tries to
% explain observed motion in terms of a few hidden states (i.e. order
% parameters) that implicitly enslave random fluctuations.
%
% THIS DEMO IS UNDER CONSTRUCTION

% non-hierarchical non-linear generative model (dynamic & chaotic)
%==========================================================================
spm_figure('GetWin','Figure 1');

% get model parameters and initial states
%--------------------------------------------------------------------------
x      = [4 4 28]';
P.A    = [-10  10  0;
           32 -1   0;
           0   0  -8/3];
P.B{1} = [ 0  0  0;
           0  0 -1;
           0  1  0];
P.B{2} = [ 0  0  0;
           0  0  0;
           0  0  0];
P.B{3} = [ 0  0  0;
           0  0  0;
           0  0  0];
P.C    = {};

V.A    = [ 1  1  1;
           1  1  1;
           1  1  1];
V.B{1} = [ 1  1  1;
           1  1  1;
           1  1  1];
V.B{2} = [ 0  1  1;
           0  1  1;
           0  1  1];
V.B{3} = [ 0  0  1;
           0  0  1;
           0  0  1];
V.C    = {};
pE     = spm_unvec(spm_vec(P)/2,P);
pC     = diag(spm_vec(V));

% level 1
%--------------------------------------------------------------------------
M(1).f  = inline('spm_fx_poly(x,v,P)/32','x','v','P');
M(1).g  = inline('x','x','v','P');
M(1).x  = x;
M(1).pE = P;
M(1).V  = exp(8);
M(1).W  = exp(8);

% level 2
%--------------------------------------------------------------------------
M(2).v  = 0;
M(2).V  = exp(16);

% create data
%==========================================================================

% create innovations & add causes
%--------------------------------------------------------------------------
N       = 128;
U       = sparse(1,N);
DEM     = spm_DEM_generate(M,U);
spm_DEM_qU(DEM.pU)

% DEM estimation
%==========================================================================
DEM.M(1).E.n  = 4;
DEM.M(1).E.s  = 1/8;
DEM.M(1).pE   = pE;
DEM.M(1).pC   = pC*exp(0);
DEM           = spm_DEM(DEM);

M(1).pE = DEM.qP.P{1};
SIM     = spm_DEM_generate(M,U);

% graphics
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');
spm_DEM_qU(SIM.pU)

subplot(2,2,3)
x      = DEM.pU.x{1};
plot(x(1,:),x(2,:))
title('true manifold dynamics')
xlabel('time')
axis square

subplot(2,2,4)
x      = SIM.pU.x{1};
plot(x(1,:),x(2,:))
title('estimated dynamics')
xlabel('time')
axis square

return

% Notes for creation of a globally coupled map
%==========================================================================

% Demo of synchronization manifold using coupled Lorenz attractors
%--------------------------------------------------------------------------
W    = 2;                            % amplitude of random fluctuations
N    = 16;                           % number of (Lorenz) oscillators
T    = 512;                          % number of time bins
dt   = 1/32;                         % time interval

% parameters (set s.d. P.t to 1/4 to see oscillator death)
%--------------------------------------------------------------------------
P.t  = randn(N,1)/16;                 % variations in log-rate constants
P.k  = 2;                            % global coupling parameter

% states
%--------------------------------------------------------------------------
x      = randn(3,N)*8;               % microstates
x(1,:) = x(1,:) + 0;
x(2,:) = x(2,:) + 0;
x(3,:) = x(3,:) + 28;
v      = 0;

% integrate
%--------------------------------------------------------------------------
for i = 1:T

    [dfdx f]  = spm_diff('spm_lorenz_k',x,v,P,1);
    f         = f + randn(N*3,1)*W;
    dx        = spm_dx(dfdx,f,dt);
    x         = x + spm_unvec(dx,x);
    y(:,i)    = x(:);               % microstates
    X(:,i)    = mean(x,2);          % macrostates

end

% plot
%--------------------------------------------------------------------------
subplot(2,1,1)
plot(y',':'), hold on
plot(X','k'), hold off
axis tight

t = T/4:T;
subplot(2,1,2)
plot(X(1,t),X(2,t),'k')
axis square
axis([-16 16 -16 16])

return

% plot synchronisation manifold
%--------------------------------------------------------------------------
clf
subplot(2,1,1)
for i = 1:12
    plot(y(i*3 - 2,t),y(i*3 + 1,t),'.','Color',[1 1/2 1/2]),hold on
end
plot(X(1,t),X(1,t),'.k'),hold off
axis square
axis([-16 16 -16 16])

author: spm


function [V_RF, V_U, iter] = mo_algorithm(V_RF, Vn, H1)

global manifold;
[Nt, Nrf] = size(V_RF);
problem.M = manifold;

problem.cost = @(x)MMSE_cost(x,H1,Vn);
problem.egrad = @(x)MMSE_egrad(x,H1,Vn);

[x,iter] = conjugategradient(problem,V_RF(:));

V_RF = reshape(x,Nt,Nrf);
V_U = inv(V_RF'*H1 * H1'* V_RF+ 1 * Vn *(V_RF)'*V_RF)*V_RF'*H1;


text:
function fileName = dnetWriteResult(model, dataSet, number)

% DNETWRITERESULT Write a DNET result.
% FORMAT
% DESC writes a DNET result.
% ARG model : the model to write.
% ARG dataSet : the name of the data set to write.
% ARG number : the number of the DNET run to write.
% RETURN fileName : the file name used to write.
%
% SEEALSO : dnetLoadResult
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS

dataSet(1) = upper(dataSet(1));
type = model.type;
type(1) = upper(type(1));
fileName = ['dem' dataSet type num2str(number)];

[mapping, dnetInfo] = dnetDeconstruct(model);

save(fileName, 'mapping', 'dnetInfo');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% REALPROBA - compute the effective probability of the value
%               in the sample.
%
% Usage:
%   >> [probaMap, probaDist ] = realproba( data, discret);
%
% Inputs:
%   data       - the data onto which compute the probability
%   discret    - discretisation factor (default: (size of data)/5)
%                if 0 base the computation on a Gaussian
%                approximation of the data
%
% Outputs:
%   probaMap   - the probabilities associated with the values
%   probaDist  - the probabilities distribution
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001

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

function [ probaMap, sortbox ] = realproba( data, bins );

if nargin < 1
	help realproba;
	return;
end
if nargin < 2
	bins = round(size(data,1)*size(data,2)/5);
end;

if bins > 0
	% COMPUTE THE DENSITY FUNCTION
	% ----------------------------
	SIZE = size(data,1)*size(data,2);
	sortbox = zeros(1,bins);
	minimum =  min(data(:));
	maximum =  max(data(:));
	data = floor((data - minimum )/(maximum - minimum)*(bins-1))+1;
    if any(any(isnan(data))), warning('Binning failed - could be due to zeroed out channel'); end
	for index=1:SIZE
		sortbox(data(index)) = sortbox(data(index))+1;
	end
	probaMap = sortbox(data) / SIZE;
	sortbox  = sortbox / SIZE;
else
	% BASE OVER ERROR FUNCTION
	% ------------------------
	data     = (data-mean(data(:)))./std(data(:));
	probaMap = exp(-0.5*( data.*data ))/(2*pi);
	probaMap = probaMap/sum(probaMap); % because the total surface under a normalized Gaussian is 2
	sortbox  = probaMap/sum(probaMap);
end
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function createSBMLFiles(refinedFolder, sbmlFolder)
% Creates SBML files for the created refined reconstructions. This may be
% time-consuming.
%
% USAGE
%   createTestSBMLFiles(refinedFolder, sbmlFolder)
%
% INPUTS
% refinedFolder            Folder with refined COBRA models generated by
%                          the refinement pipeline
% sbmlFolder               Folder where SBML files, if desired, will be saved
%
% AUTHOR:
%       - Almut Heinken, 03/2020

dInfo = dir(sbmlFolder);
alreadyCreated={dInfo.name};
alreadyCreated=alreadyCreated';
alreadyCreated(~contains(alreadyCreated(:,1),'.xml'),:)=[];
alreadyCreated=strrep(alreadyCreated,'.xml','.mat');

dInfo = dir(refinedFolder);
modelList={dInfo.name};
modelList=modelList';
modelList(find(strcmp(modelList(:,1),'.')),:)=[];
modelList(find(strcmp(modelList(:,1),'..')),:)=[];

for i=1:length(modelList)
    if ~any(strcmp(modelList{i},alreadyCreated))
    model = readCbModel([refinedFolder filesep modelList{i}]);
    writeCbModel(model, 'format', 'sbml', 'fileName', [sbmlFolder filesep strrep(modelList{i},'.mat','')]);
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function [varargout] = iscolumn(varargin)

% ISCOLUMNN is a drop-in replacement for the same function that was
% introduced in MATLAB R2010b.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to automatically remove incorrect compat folders from the path, see https://github.com/fieldtrip/fieldtrip/issues/899

if isempty(strfind(mfilename('fullpath'), matlabroot))
  % automatic cleanup does not work when FieldTrip is installed inside the MATLAB path, see https://github.com/fieldtrip/fieldtrip/issues/1527

  alternatives = which(mfilename, '-all');
  if ~iscell(alternatives)
    % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
    alternatives = {alternatives};
  end

  keep = true(size(alternatives));
  for i=1:numel(alternatives)
    keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
    keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
  end
  alternatives = alternatives(keep);

  if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
    % remove this directory from the path
    p = fileparts(mfilename('fullpath'));
    warning('removing "%s" from your path, see http://bit.ly/2SPPjUS', p);
    rmpath(p);
    % call the original MATLAB function
    if exist(mfilename, 'builtin')
      [varargout{1:nargout}] = builtin(mfilename, varargin{:});
    else
      [varargout{1:nargout}] = feval(mfilename, varargin{:});
    end
    return
  end

end % automatic cleanup of compat directories

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = iscolumn(x)

% deal with the input arguments
if nargin==1
  [x] = deal(varargin{1:1});
else
  error('incorrect number of input arguments')
end

tf = length(size(x))==2 && size(x,2)==1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [m,v,w,g,f,pp,gg]=gaussmix(x,c,l,m0,v0,w0)
%GAUSSMIX fits a gaussian mixture pdf to a set of data observations [m,v,w,g,f]=(x,c,l,m0,v0,w0)
%
% Inputs: n data values, k mixtures, p parameters, l loops
%
%     X(n,p)   Input data vectors, one per row.
%     c(1)     Minimum variance of normalized data (Use [] to take default value of 1/n^2)
%     L        The integer portion of l gives a maximum loop count. The fractional portion gives
%              an optional stopping threshold. Iteration will cease if the increase in
%              log likelihood density per data point is less than this value. Thus l=10.001 will
%              stop after 10 iterations or when the increase in log likelihood falls below
%              0.001.
%              As a special case, if L=0, then the first three outputs are omitted.
%              Use [] to take default value of 100.0001
%     M0(k,p)  Initial mixture means, one row per mixture.
%     V0(k,p)  Initial mixture variances, one row per mixture.
%      or V0(p,p,k)  one full-covariance matrix per mixture
%     W0(k,1)  Initial mixture weights, one per mixture. The weights should sum to unity.
%
%     Alternatively, if initial values for M0, V0 and W0 are not given explicitly:
%
%     M0       Number of mixtures required
%     V0       Initialization mode:
%                'f'    Initialize with K randomly selected data points [default]
%                'p'    Initialize with centroids and variances of random partitions
%                'k'    k-means algorithm ('kf' and 'kp' determine initialization of kmeans)
%                'h'    k-harmonic means algorithm ('hf' and 'hp' determine initialization of kmeans)
%                's'    do not scale data during initialization to have equal variances
%                'm'    M0 contains the initial centres
%                'v'    full covariance matrices
%              Mode 'hf' [the default] generally gives the best results but 'f' is faster and often OK
%
% Outputs: (Note that M, V and W are omitted if L==0)
%
%     M(k,p)   Mixture means, one row per mixture. (omitted if L==0)
%     V(k,p)   Mixture variances, one row per mixture. (omitted if L==0)
%       or V(p,p,k) if full covariance matrices in use (i.e. either 'v' option or V0(p,p,k) specified)
%     W(k,1)   Mixture weights, one per mixture. The weights will sum to unity. (omitted if L==0)
%     G       Average log probability of the input data points.
%     F        Fisher's Discriminant measures how well the data divides into classes.
%              It is the ratio of the between-mixture variance to the average mixture variance: a
%              high value means the classes (mixtures) are well separated.
%     PP(n,1)  Log probability of each data point
%     GG(l+1,1) Average log probabilities at the beginning of each iteration and at the end

%  Bugs/Suggestions
%     (2) Allow processing in chunks by outputting/reinputting an array of sufficient statistics
%     (6) Other initialization options:
%              'l'    LBG algorithm
%              'm'    Move-means (dog-rabbit) algorithm

%      Copyright (C) Mike Brookes 2000-2009
%      Version: $Id: gaussmix.m,v 1.22 2009/09/21 14:03:13 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,p]=size(x);
mx0=sum(x,1)/n;         % calculate mean and variance of input data
vx0=sum(x.^2,1)/n-mx0.^2;
sx0=sqrt(vx0);
sx0(sx0==0)=1;      % do not divide by zero when scaling
scaled=0;           % data is not yet scaled
memsize=voicebox('memsize');    % set memory size to use
if isempty(c)
    c=1/n^2;
else
    c=c(1);         % just to prevent legacy code failing
end
fulliv=0;           % initial variance is not full
if isempty(l)
    l=100+1e-4;         % max loop count + stopping threshold
end
if nargin<6             % no initial values specified for m0, v0, w0
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  No initialvalues given, so we must use k-means or equivalent
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if nargin<5
        v0='hf';                 % default initialization mode: hf
    end
    if any(v0=='m')
        k=size(m0,1);
    else
        k=m0;
    end
    fv=any(v0=='v');                % full covariance matrices requested
    if n<=k                         % each data point can have its own mixture
        xs=(x-mx0(ones(n,1),:))./sx0(ones(n,1),:);          % scale the data
        m=xs(mod((1:k)-1,n)+1,:);    % just include all points several times
        v=zeros(k,p);               % will be set to floor later
        w=zeros(k,1);
        w(1:n)=1/n;
        if l>0
            l=0.1;                  % no point in iterating
        end
    else                            % more points than mixtures
        if any(v0=='s')
            xs=x;                   % do not scale data during initialization
        else
            xs=(x-mx0(ones(n,1),:))./sx0(ones(n,1),:);  % else scale now
            if any(v0=='m')
                m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:);  % scale specified means as well
            end
        end
        w=repmat(1/k,k,1);                      % all mixtures equally likely
        if any(v0=='k')                         % k-means initialization
            if any(v0=='m')
                [m,e,j]=kmeans(xs,k,m);
            elseif any(v0=='p')
                [m,e,j]=kmeans(xs,k,'p');
            else
                [m,e,j]=kmeans(xs,k,'f');
            end
        elseif any(v0=='h')                     % k-harmonic means initialization
            if any(v0=='m')
                [m,e,j]=kmeanhar(xs,k,[],4,m);
            else
                if any(v0=='p')
                    [m,e,j]=kmeanhar(xs,k,[],4,'p');
                else
                    [m,e,j]=kmeanhar(xs,k,[],4,'f');
                end
            end
        elseif any(v0=='p')                     % Initialize using a random partition
            j=ceil(rand(n,1)*k);                % allocate to random clusters
            j(rnsubset(k,n))=1:k;               % but force at least one point per cluster
            for i=1:k
                m(i,:)=mean(xs(j==i,:),1);
            end
        else
            if any(v0=='m')
                m=m0;                           % use specified centres
            else
                m=xs(rnsubset(k,n),:);          % Forgy initialization: sample k centres without replacement [default]
            end
            [e,j]=kmeans(xs,k,m,0);             % find out the cluster allocation
        end
        if any(v0=='s')
            xs=(x-mx0(ones(n,1),:))./sx0(ones(n,1),:);      % scale data now if not done previously
        end
        v=zeros(k,p);                   % diagonal covariances
        w=zeros(k,1);
        for i=1:k
            ni=sum(j==i);               % number assigned to this centre
            w(i)=(ni+1)/(n+k);          % weight of this mixture
            if ni
                v(i,:)=sum((xs(j==i,:)-repmat(m(i,:),ni,1)).^2,1)/ni;
            else
                v(i,:)=zeros(1,p);
            end
        end
    end
else
    %%%%%%%%%%%%%%%%%%%%%%%%
    % use initial values given as input parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [k,p]=size(m0);
    xs=(x-mx0(ones(n,1),:))./sx0(ones(n,1),:);          % scale the data
    m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:);          % and the means
    v=v0;
    w=w0;
    fv=ndims(v)>2 || size(v,1)>k;                       % full covariance matrix is supplied
    if fv
        mk=eye(p)==0;                                    % off-diagonal elements
        fulliv=any(v(repmat(mk,[1 1 k]))~=0);            % check if any are non-zero
        if ~fulliv
            v=reshape(v(repmat(~mk,[1 1 k])),p,k)'./repmat(sx0.^2,k,1);   % just pick out and scale the diagonal elements for now
        else
            v=v./repmat(sx0'*sx0,[1 1 k]);              % scale the full covariance matrix
        end
    end
end
lsx=sum(log(sx0));
if ~fulliv          % initializing with diagonal covariance
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Diagonal Covariance matrices  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    v=max(v,c);         % apply the lower bound
    xs2=xs.^2;          % square the data for variance calculations

    % If data size is large then do calculations in chunks

    nb=min(n,max(1,floor(memsize/(8*p*k))));    % chunk size for testing data points
    nl=ceil(n/nb);                  % number of chunks
    jx0=n-(nl-1)*nb;                % size of first chunk

    im=repmat(1:k,1,nb); im=im(:);
    th=(l-floor(l))*n;
    sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values
    lp=floor(l)+sd;   % extra loop needed to calculate final G value

    lpx=zeros(1,n);             % log probability of each data point
    wk=ones(k,1);
    wp=ones(1,p);
    wnb=ones(1,nb);
    wnj=ones(1,jx0);

    % EM loop

    g=0;                           % dummy initial value for comparison
    gg=zeros(lp+1,1);
    ss=sd;                       % initialize stopping count (0 or 1)
    for j=1:lp
        g1=g;                    % save previous log likelihood (2*pi factor omitted)
        m1=m;                       % save previous means, variances and weights
        v1=v;
        w1=w;
        vi=-0.5*v.^(-1);                % data-independent scale factor in exponent
        lvm=log(w)-0.5*sum(log(v),2);   % log of external scale factor (excluding -0.5*p*log(2pi) term)

        % first do partial chunk

        jx=jx0;
        ii=1:jx;
        kk=repmat(ii,k,1);
        km=repmat(1:k,1,jx);
        py=reshape(sum((xs(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx)+lvm(:,wnj);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
        px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
        lpx(ii)=log(ps)+mx;
        pk=sum(px,2);                   % effective number of data points for each mixture (could be zero due to underflow)
        sx=px*xs(ii,:);
        sx2=px*xs2(ii,:);

        for il=2:nl
            ix=jx+1;
            jx=jx+nb;                    % increment upper limit
            ii=ix:jx;
            kk=repmat(ii,k,1);
            py=reshape(sum((xs(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb)+lvm(:,wnb);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
            lpx(ii)=log(ps)+mx;
            pk=pk+sum(px,2);                   % effective number of data points for each mixture (could be zero due to underflow)
            sx=sx+px*xs(ii,:);
            sx2=sx2+px*xs2(ii,:);
        end
        g=sum(lpx);                    % total log probability summed over all data points
        gg(j)=g;
        w=pk/n;                         % normalize to get the weights
        if pk                       % if all elements of pk are non-zero
            m=sx./pk(:,wp);
            v=sx2./pk(:,wp);
        else
            wm=pk==0;                       % mask indicating mixtures with zero weights
            nz=sum(wm);                  % number of zero-weight mixtures
            [vv,mk]=sort(lpx);             % find the lowest probability data points
            m=zeros(k,p);                   % initialize means and variances to zero (variances are floored later)
            v=m;
            m(wm,:)=xs(mk(1:nz),:);                % set zero-weight mixture means to worst-fitted data points
            w(wm)=1/n;                      % set these weights non-zero
            w=w*n/(n+nz);                   % normalize so the weights sum to unity
            wm=~wm;                         % mask for non-zero weights
            m(wm,:)=sx(wm,:)./pk(wm,wp);  % recalculate means and variances for mixtures with a non-zero weight
            v(wm,:)=sx2(wm,:)./pk(wm,wp);
        end
        v=max(v-m.^2,c);                % apply floor to variances

        if g-g1<=th && j>1
            if ~ss, break; end  %  stop
            ss=ss-1;       % stop next time
        end

    end
    if sd && ~fv  % we need to calculate the final probabilities
        pp=lpx'-0.5*p*log(2*pi)-lsx;   % log of total probability of each data point
        gg=gg(1:j)/n-0.5*p*log(2*pi)-lsx;    % average log prob at each iteration
        g=gg(end);
        %     gg' % *** DEBUG ***
        m=m1;       % back up to previous iteration
        v=v1;
        w=w1;
        mm=sum(m,1)/k;
        f=(m(:)'*m(:)-k*mm(:)'*mm(:))/sum(v(:));
    end
    if ~fv
        m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:);  % unscale means
        v=v.*repmat(sx0.^2,k,1);                % and variances
    else
        v1=v;
        v=zeros(p,p,k);
        mk=eye(p)==1;                                    % mask for diagonal elements
        v(repmat(mk,[1 1 k]))=v1';            % set from v1
    end
end
if fv              % check if full covariance matrices were requested
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Full Covariance matrices  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    pl=p*(p+1)/2;
    lix=1:p^2;
    cix=repmat(1:p,p,1);
    rix=cix';
    lix(cix>rix)=[];                                        % index of lower triangular elements
    cix=cix(lix);                                           % index of lower triangular columns
    rix=rix(lix);                                           % index of lower triangular rows
    dix=find(rix==cix);
    lixi=zeros(p,p);
    lixi(lix)=1:pl;
    lixi=lixi';
    lixi(lix)=1:pl;                                        % reverse index to build full matrices
    v=reshape(v,p^2,k);
    v=v(lix,:)';                                            % lower triangular in rows

    % If data size is large then do calculations in chunks

    nb=min(n,max(1,floor(memsize/(24*p*k))));    % chunk size for testing data points
    nl=ceil(n/nb);                  % number of chunks
    jx0=n-(nl-1)*nb;                % size of first chunk
    %
    th=(l-floor(l))*n;
    sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values
    lp=floor(l)+sd;   % extra loop needed to calculate final G value
    %
    lpx=zeros(1,n);             % log probability of each data point
    wk=ones(k,1);
    wp=ones(1,p);

    wpl=ones(1,pl);       % 1 index for lower triangular matrix
    wnb=ones(1,nb);
    wnj=ones(1,jx0);

    % EM loop

    g=0;                           % dummy initial value for comparison
    gg=zeros(lp+1,1);
    ss=sd;                       % initialize stopping count (0 or 1)
    vi=zeros(p*k,p);                    % stack of k inverse cov matrices each size p*p
    vim=zeros(p*k,1);                   % stack of k vectors of the form inv(v)*m
    mtk=vim;                             % stack of k vectors of the form m
    lvm=zeros(k,1);
    wpk=repmat((1:p)',k,1);
    for j=1:lp
        g1=g;                    % save previous log likelihood (2*pi factor omitted)
        m1=m;                       % save previous means, variances and weights
        v1=v;
        w1=w;

        for ik=1:k

            % these lines added for debugging only
            %             vk=reshape(v(k,lixi),p,p);
            %             condk(ik)=cond(vk);
            %%%%%%%%%%%%%%%%%%%%
            [uvk,dvk]=eig(reshape(v(ik,lixi),p,p));      % convert lower triangular to full and find eigenvalues
            dvk=max(diag(dvk),c);                           % apply variance floor to eigenvalues
            vik=-0.5*uvk*diag(dvk.^(-1))*uvk';   % calculate inverse
            vi((ik-1)*p+(1:p),:)=vik;           % vi contains all mixture inverses stacked on top of each other
            vim((ik-1)*p+(1:p))=vik*m(ik,:)';   % vim contains vi*m for all mixtures stacked on top of each other
            mtk((ik-1)*p+(1:p))=m(ik,:)';       % mtk contains all mixture means stacked on top of each other
            lvm(ik)=log(w(ik))-0.5*sum(log(dvk));       % vm contains the weighted sqrt of det(vi) for each mixture
        end
        %
        %         % first do partial chunk
        %
        jx=jx0;
        ii=1:jx;
        xii=xs(ii,:).';
        py=reshape(sum(reshape((vi*xii-vim(:,wnj)).*(xii(wpk,:)-mtk(:,wnj)),p,jx*k),1),k,jx)+lvm(:,wnj);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));  % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
        px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
        lpx(ii)=log(ps)+mx;
        pk=sum(px,2);                   % effective number of data points for each mixture (could be zero due to underflow)
        sx=px*xs(ii,:);
        sx2=px*(xs(ii,rix).*xs(ii,cix));            % accumulator for variance calculation (lower tri cov matrix as a row)

        for il=2:nl
            ix=jx+1;
            jx=jx+nb;        % increment upper limit
            ii=ix:jx;
            xii=xs(ii,:).';
            py=reshape(sum(reshape((vi*xii-vim(:,wnb)).*(xii(wpk,:)-mtk(:,wnb)),p,nb*k),1),k,nb)+lvm(:,wnb);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));  % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
            lpx(ii)=log(ps)+mx;
            pk=pk+sum(px,2);                % effective number of data points for each mixture (could be zero due to underflow)
            sx=sx+px*xs(ii,:);               % accumulator for mean calculation
            sx2=sx2+px*(xs(ii,rix).*xs(ii,cix));            % accumulator for variance calculation
        end
        g=sum(lpx);                    % total log probability summed over all data points
        gg(j)=g;                        % save convergence history
        w=pk/n;                         % normalize to get the column of weights
        if pk                       % if all elements of pk are non-zero
            m=sx./pk(:,wp);         % find mean and mean square
            v=sx2./pk(:,wpl);
        else
            wm=pk==0;                       % mask indicating mixtures with zero weights
            nz=sum(wm);                  % number of zero-weight mixtures
            [vv,mk]=sort(lpx);             % find the lowest probability data points
            m=zeros(k,p);                   % initialize means and variances to zero (variances are floored later)
            v=zeros(k,pl);
            m(wm,:)=xs(mk(1:nz),:);                % set zero-weight mixture means to worst-fitted data points
            w(wm)=1/n;                      % set these weights non-zero
            w=w*n/(n+nz);                   % normalize so the weights sum to unity
            wm=~wm;                         % mask for non-zero weights
            m(wm,:)=sx(wm,:)./pk(wm,wp);  % recalculate means and variances for mixtures with a non-zero weight
            v(wm,:)=sx2(wm,:)./pk(wm,wpl);
        end
        v=v-m(:,cix).*m(:,rix);                 % subtract off mean squared
        if g-g1<=th && j>1
            if ~ss, break; end  %  stop
            ss=ss-1;       % stop next time
        end
    end
    if sd  % we need to calculate the final probabilities
        pp=lpx'-0.5*p*log(2*pi)-lsx;   % log of total probability of each data point
        gg=gg(1:j)/n-0.5*p*log(2*pi)-lsx;    % average log prob at each iteration
        g=gg(end);
        %             gg' % *** DEBUG ONLY ***
        m=m1;       % back up to previous iteration
        v=zeros(p,p,k);
        trv=0;      % sum of variance matrix traces
        for ik=1:k
            [uvk,dvk]=eig(reshape(v1(ik,lixi),p,p));      % convert lower triangular to full and find eigenvectors
            dvk=max(diag(dvk),c);                       % apply variance floor
            v(:,:,ik)=uvk*diag(dvk)*uvk';          % reconstitute full matrix
            trv=trv+sum(dvk);
        end
        w=w1;
        mm=sum(m,1)/k;
        f=(m(:)'*m(:)-k*mm(:)'*mm(:))/trv;
    else
        v1=v;               % lower triangular form
        v=zeros(p,p,k);
        for ik=1:k
            [uvk,dvk,]=eig(reshape(v1(ik,lixi),p,p));      % convert lower triangular to full and find eigenvectors
            dvk=max(diag(dvk),c);                       % apply variance floor
            v(:,:,ik)=uvk*diag(dvk)*uvk';          % reconstitute full matrix
        end
    end
    m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:);  % unscale means
    v=v.*repmat(sx0'*sx0,[1 1 k]);
end
if l==0         % suppress the first three output arguments if l==0
    m=g;
    v=f;
    w=pp;
end

author: decouples


text:
function d=disteusq(x,y,mode,w)
%DISTEUSQ calculate euclidean, squared euclidean or mahanalobis distance D=(X,Y,MODE,W)
%
% Inputs: X,Y         Vector sets to be compared. Each row contains a data vector.
%                     X and Y must have the same number of columns.
%
%         MODE        Character string selecting the following options:
%                         'x'  Calculate the full distance matrix from every row of X to every row of Y
%                         'd'  Calculate only the distance between corresponding rows of X and Y
%                              The default is 'd' if X and Y have the same number of rows otherwise 'x'.
%                         's'  take the square-root of the result to give the euclidean distance.
%
%         W           Optional weighting matrix: the distance calculated is (x-y)*W*(x-y)'
%                     If W is a vector, then the matrix diag(W) is used.
%
% Output: D           If MODE='d' then D is a column vector with the same number of rows as the shorter of X and Y.
%                     If MODE='x' then D is a matrix with the same number of rows as X and the same number of columns as Y'.
%

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: disteusq.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nx,p]=size(x); ny=size(y,1);
if nargin<3 | isempty(mode) mode='0'; end
if any(mode=='d') | (mode~='x' & nx==ny)

    % Do pairwise distance calculation

    nx=min(nx,ny);
    z=double(x(1:nx,:))-double(y(1:nx,:));
    if nargin<4
        d=sum(z.*conj(z),2);
    elseif min(size(w))==1
        wv=w(:).';
        d=sum(z.*wv(ones(size(z,1),1),:).*conj(z),2);
    else
        d=sum(z*w.*conj(z),2);
    end
else

    % Calculate full distance matrix

    if p>1

        % x and y are matrices

        if nargin<4
            z=permute(double(x(:,:,ones(1,ny))),[1 3 2])-permute(double(y(:,:,ones(1,nx))),[3 1 2]);
            d=sum(z.*conj(z),3);
        else
            nxy=nx*ny;
            z=reshape(permute(double(x(:,:,ones(1,ny))),[1 3 2])-permute(double(y(:,:,ones(1,nx))),[3 1 2]),nxy,p);
            if min(size(w))==1
                wv=w(:).';
                d=reshape(sum(z.*wv(ones(nxy,1),:).*conj(z),2),nx,ny);
            else
                d=reshape(sum(z*w.*conj(z),2),nx,ny);
            end
        end
    else

        % x and y are vectors

        z=double(x(:,ones(1,ny)))-double(y(:,ones(1,nx))).';
        if nargin<4
            d=z.*conj(z);
        else
            d=w*z.*conj(z);
        end
    end
end
if any(mode=='s')
    d=sqrt(d);
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function news ( filename, thresh )

%*****************************************************************************80
%
%% NEWS demonstrates the NEWS stencil for edge detection.
%
%  Discussion:
%
%    Given a black and white image A, which we regard as an M by N array
%    of pixels, we want to produce an array E of the same shape, which
%    contains information describing the location of edges.
%
%    A simple algorithm for trying to detect edges in an array that
%    represents an image is the NEWS scheme.  For each pixel A(C),
%    we consider its North, East, West, and South pixel neighbors.  The
%    indexing of arrays and images do not correspond, so we will use
%    these directions instead:
%
%             A(N)
%              |
%              |
%      A(W)---A(C)---A(E)
%              |
%              |
%             A(S)
%
%    Entry E(C) of the edge array will be computed by
%
%      E(C) = abs ( A(N) - A(S) ) + abs ( A(E) - A(W) )
%
%    Pixels of A that represent edges will tend to have high values
%    of E, while pixels that are interior to a region of roughly the
%    same shade will tend to have low values.
%
%    Thus, an edge detection scheme would use the NEWS stencil to
%    compute the E array, determine E_MAX, the maximum entry in E,
%    choose some threshold value E_THRESH, and declare pixel A(I,J)
%    to be associated with an edge whenever E(I,J) is greater than E_THRESH.
%
%    In this program, we demonstrate the NEWS stencil using a PGM
%    grayscale image of coins.  At the end, we use the edge information
%    to produce a color image in which the edges of the coins have been
%    outlined in red.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the name of the file containing the image.
%    The image should be a grayscale image, not color!
%
%    Input, integer THRESH, the threshhold for the edge detection.
%    0 <= THRESH <= 255 is required.  A value of 50 is the default.
%    Higher values are more selective.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NEWS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate the NEWS stencil for edge detection\n' );
  fprintf ( 1, '  in images.\n' );

  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    filename = input ( 'Enter the name of the image file:  ' );
  end

  if ( nargin < 2 )
    thresh = 50;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reading "%s".\n', filename );
%
%  Read the file into "A".
%
  a = imread ( filename );
%
%  Get the size of the array.
%  We'll use it later when we add a border to the data.
%
  [ m, n ] = size ( a );
%
%  Display the input image.
%
  figure ( 1 )
  imshow ( a );
  title ( 'Initial image' )

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Figure 1:\n' );
  fprintf ( 1, ' This is the original gray scale image.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Press return.\n' );
  pause
%
%  For neatness, we add a border of zeros to the image,
%  then fill in the border by copying the nearby original values.
%  This will be our M+2 by N+2 data array B.
%
  b = zeros ( m + 2, n + 2 );
  b(2:m+1,2:n+1) = double ( a );

  b(1,  2:n+1) = b(2,2:n+1);
  b(m+2,2:n+1) = b(m+1,2:n+1);

  b(2:m+1,1)   = b(2:m+1,2);
  b(2:m+1,n+2) = b(2:m+1,n+1);

  b(1,1)       = ( b(1,2)     + b(2,1)     ) / 2.0;
  b(m+2,1)     = ( b(m+2,2)   + b(m+1,1)   ) / 2.0;
  b(1,n+2)     = ( b(1,n+1)   + b(2,n+2)   ) / 2.0;
  b(m+2,n+2)   = ( b(m+2,n+1) + b(m+1,n+2) ) / 2.0;

  figure ( 2 )
  imshow ( uint8 ( b ) );
  title ( 'Same image, with a 1-pixel border.' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Figure 2:\n' );
  fprintf ( 1, '  We have added a border of 1-pixel thickness.\n' );
  fprintf ( 1, '  Now all the pixels in the original picture \n' );
  fprintf ( 1, '  have a 3x3 pixel neighborhood.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Press return.\n' );
  pause
%
%  Apply the NEWS Operator.  We don't process the boundary pixels.
%
%  The stencil we use is:
%
%   |  0 +1  0 |     |  0  0   0 |
%   |  0  0  0 |  +  | -1  0  +1 |
%   |  0 -1  0 |     |  0  0   0 |
%
  e(2:m+1,2:n+1) = abs ( - b(1:m,2:n+1) + b(3:m+2,2:n+1) ) ...
                 + abs ( - b(2:m+1,1:n) + b(2:m+1,3:n+2) );
%
%  The values in E must be rescaled to run from 0 to 255
%  if we are going to display them.
%
  emin = min ( min ( e ) );
  emax = max ( max ( e ) );
  e = round ( 255 * ( e - emin ) / ( emax - emin ) );

  figure ( 3 )
  imshow ( uint8 ( e ) );
  title ( 'All the E data.' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Figure 3:\n' );
  fprintf ( 1, '  We computed the value of E for each pixel,\n' );
  fprintf ( 1, '  and scaled it to [0,255].\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Press return.\n' );
  pause
%
%  Threshold the data.
%
  e = 255 * ( thresh < e );

  figure ( 4 )
  imshow ( uint8 ( e ) );
  title ( 'E data above the threshold.' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Figure 4:\n' );
  fprintf ( 1, '  We zeroed all the pixels with a value of E\n' );
  fprintf ( 1, '  less than the threshold of %d\n', thresh )
  fprintf ( 1, '  so the edges show up better.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Press return.\n' );
  pause
%
%  Reverse the image.
%
  e_reverse = 255 - e;

  figure ( 5 )
  imshow ( uint8 ( e_reverse ) );
  title ( 'E data above the threshold (reverse video).' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Figure 5:\n' );
  fprintf ( 1, '  Black lines on white background are MUCH easier\n' );
  fprintf ( 1, '  to read.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Press return.\n' );
  pause
%
%  Make an RGB image of the gray picture, with red highlights.
%
  e2 = max ( e(2:m+1,2:n+1), double ( a ) );
  a2 = uint8 ( e2 );

  r = a2;
  g = a;
  b = a;

  rgb = cat ( 3, r, g, b );

  figure ( 6 )
  imshow ( rgb );
  title ( 'Original gray data, with edges in red.' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Figure 6:\n' );
  fprintf ( 1, '  This is actually an RGB "color" image.\n' );
  fprintf ( 1, '  This way, we can show the detected edges in red.\n' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NEWS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function [bandsLin2D, bandsLinAv, ...
    mB2D, mBAv, ...
    names2D, namesAv, ...
    namesmB2D, namesmBAv] = extractBandsLin(data, params)
% Do FFT over whole segement and return average bin power
% Later maybe add time divisions

% Delta: 1-3Hz
% Theta: 4-7Hz
% Alpha1: 8-9Hz
% Alpha2: 10-12Hz
% Beta1: 13-17Hz
% Beta2: 18-30Hz
% Gamma1: 31-40Hz
% Gamma2: 41-50Hz
% Higher: 51->200 Hz
bLims = [[1;3], [4;7], [8;9], [10;12], [13;17], [18;30], [31;40], ...
    [41;50], [51;70], [71;150], [151;250]]; % Last 3 mod M37
nBands = size(bLims,2);
nChans = 16;

bands2D = NaN(nBands, 16, 'single');

mB2D =  NaN(1, 16, 'single');

for c = 1:nChans

    Fs = 400;
    T = 1/Fs;
    L = size(data,1);

    Y = fft(data(:,c));

    P2 = abs(Y/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);

    f = Fs*(0:(L/2))/L;

    if params.plotOn
        plot(f,P1)
        title('Single-Sided Amplitude Spectrum of X(t)')
        xlabel('f (Hz)')
        ylabel('|P1(f)|')
    end

    for b = 1:nBands
        bIdx = f>=bLims(1,b) & f<=bLims(2,b);

        mPower = mean(P1(bIdx));

        bands2D(b, c) = mPower;
    end

    [~, mB2D(1,c)] = max(bands2D(:,c));
end

% Max bands
mBAv = single(mean(mB2D));
namesmB2D = (string('maxBand_c') + (1:16)')';
namesmBAv = 'maxBandAv';

% Bands lin
bandsLinAv = mean(bands2D,2)';
bandsLin2D = reshape(bands2D, 1, nChans*nBands);

names2D = cellstr([repmat('bandsLin2D_b', nChans*nBands,1), ...
    num2str(repmat((1:nBands)', nChans, 1)), ...
    repmat('_c', nChans*nBands, 1), ...
    num2str(reshape(repmat((1:16), nBands, 1), nBands*nChans,1))]);
names2D = strrep(names2D, ' ' , '')';

namesAv = cellstr([repmat('bandsLinAv_b', nBands,1), ...
    num2str((1:nBands)')])';

author: garethjns

repo: Kaggle-EEG


text:
% PURPOSE: demo of chowlin()
%          Temporal disaggregation with indicators.
% 			  Chow-Lin method
%---------------------------------------------------
% USAGE: chowlin_d
%---------------------------------------------------

close all; clear all; clc;

% Low-frequency data: Spain's Exports of Goods. 1995 prices

Y=[  20499
     23477
     25058
     27708
     31584
     31898
     30233
     32235
     34049
     36035
     39795
     44299
     47426
     52339
     62949
     69885
     77174
     90133
     96496
    102776
    113026
    115573 ];

% High-frequency data: Spain's Registered exports of goods deflated by
%                      unit value index.

x=[   5162
      5054
      4049
      5196
      4972
      5606
      5844
      6196
      6526
      5671
      5631
      6510
      6575
      6797
      5973
      6796
      8404
      8260
      7058
      7403
      7934
      7762
      7087
      8659
      7471
      8082
      6700
      8117
      8271
      8336
      7698
      8372
      9120
      8911
      8035
      8613
      9725
      9529
      7774
      9295
     10357
     10372
      9056
     10812
     11989
     11839
      9686
     11736
     12878
     12211
     10278
     12321
     13267
     12973
     11268
     15008
     16565
     15641
     13684
     17254
     18613
     17774
     14966
     18543
     19287
     19399
     17299
     21065
     20687
     23215
     21382
     24935
     24256
     25558
     21680
     24951
     25284
     26149
     23344
     27754
     28271
     29835
     26148
     30917
     30494
     30486
     26153
     29930 ];

% ---------------------------------------------
% Inputs for td library

% Type of aggregation
ta=1;
% Frequency conversion
sc=4;
% Method of estimation
type=1;
% Intercept
opC = -1;
% Interval of rho for grid search
% rl = [-.33 .80];
rl = 0.57;
% rl = [];
% Name of ASCII file for output
file_sal='td.sal';
% Calling the function: output is loaded in a structure called res
res=chowlin(Y,x,ta,sc,type,opC,rl);
% Calling printing function
tdprint(res,file_sal);
edit td.sal;
% Calling graph function
tdplot(res);


text:
% VTKQUIVER creates a vtk-file containing a glyph plot similar to Matlab's
%    quiver3 which can be called by a vtk-visualization tool.
%
%    VTKQUIVER(X,Y,Z,U,V,W,VARNAME,FILENAME) plots vectors as glyphs with
%    components (u,v,w) at the points (x,y,z).  The matrices X,Y,Z,U,V,W
%    must all be the same size and contain the corresponding position and
%    velocity components.  VARNAME should contain the description of the
%    visualized variable which is required by the vtk file format. FILENAME
%    is the name of the .vtu-file to store the data.  If the string does
%    not end with '.vtu', this extension is attached.
%
%    VTKQUIVER(X,Y,[],U,V,[],VARNAME,FILENAME) can be used to abbreviate
%    VTKQUIVER(X,Y,zeros(size(X)),U,V,zeros(size(X)),VARNAME,FILENAME), ie.
%    for 2d plots (vtk expects 3d data).
%
%    To visualize the glyphs you have to load the created .vtu-file via
%    paraview and hit the big `Glyph'-button.
%
% Example
%    [X,Y] = meshgrid(-2:0.25:2,-1:0.2:1);
%    Z = X.*exp(-X.^2-Y.^2);
%    [U,V,W] = surfnorm(X,Y,Z);
%    vtkquiver(X,Y,Z,U,V,W,'velocity','vtkquiver.vtu');
%    % open `vtkquiver.vtu' via Paraview or Mayavi2 and choose `Glyph'
%
% See also quiver, quiver3, vtktrisurf
%
% Copyright see license.txt
%
% Author:     Florian Frank
% eMail:      snflfran@gmx.net
% Version:    1.00
% Date:       Jun 16th, 2010

function vtkquiver(x, y, z, u, v, w, varname, filename)

% ASSERTIONS
assert(nargin == 8)
assert(ischar(varname) && ischar(filename))
dim = length(x(:));
if isempty(z); z = zeros(size(x)); end
if isempty(w); w = zeros(size(x)); end
assert(dim == length(y(:)) && dim == length(z(:)) && ...
  dim == length(u(:)) && dim == length(v(:)) && dim == length(w(:)))

% OPEN FILE
if ~strcmp(filename(end-3:end), '.vtu') % append file extension if not specified yet
  filename = [filename '.vtu'];
end

file = fopen(filename, 'wt');

% HEADER
fprintf(file, '<?xml version="1.0"?>\n');
fprintf(file, '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" compressor="vtkZLibDataCompressor">\n');
fprintf(file, '  <UnstructuredGrid>\n');
fprintf(file, '    <Piece NumberOfPoints="%d" NumberOfCells="%d">\n', dim, 0);

% POINTS
fprintf(file, '      <Points>\n');
fprintf(file, '        <DataArray type="Float32" NumberOfComponents="3" format="ascii">\n');
for k = 1 : dim
  fprintf(file, '          %.3e %.3e %.3e\n', x(k), y(k), z(k));
end
fprintf(file, '        </DataArray>\n');
fprintf(file, '      </Points>\n');

% CELLS (empty but required for paraview)
fprintf(file, '      <Cells>\n');
fprintf(file, '        <DataArray type="Int32" Name="connectivity" format="ascii">\n');
fprintf(file, '        </DataArray>\n');
fprintf(file, '        <DataArray type="Int32" Name="offsets" format="ascii">\n');
fprintf(file, '        </DataArray>\n');
fprintf(file, '        <DataArray type="UInt8" Name="types" format="ascii">\n');
fprintf(file, '        </DataArray>\n');
fprintf(file, '      </Cells>\n');

% % CELL DATA (empty but required for paraview)
% fprintf(file, '      <CellData>\n');
%   fprintf(file, '        <DataArray type="Float32" Name="void" format="ascii">\n');
%   fprintf(file, '        </DataArray>\n');
% fprintf(file, '      </CellData>\n');

% POINT DATA
fprintf(file, '      <PointData Vectors="%s">\n', varname); % def of std value
fprintf(file, '        <DataArray type="Float32" Name="%s" NumberOfComponents="3" format="ascii">\n', varname);
for k = 1 : dim
  fprintf(file, '          %.3e %.3e %.3e\n', u(k), v(k), w(k));
end
fprintf(file, '        </DataArray>\n');
fprintf(file, '      </PointData>\n');

% FOOTER
fprintf(file, '    </Piece>\n');
fprintf(file, '  </UnstructuredGrid>\n');
fprintf(file, '</VTKFile>\n');

fclose(file);

return

end


text:
function M = hmxSpy(varargin)
%+========================================================================+
%|                                                                        |
%|         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       : hmxSpy.m                                      |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Spy H-Matrix architecture                     |
%|  `---'  |                                                              |
%+========================================================================+

% Input data
Mh = varargin{1};

% Initialize visu block
if (Mh.typ > 0)
    M        = sparse(size(Mh,1),size(Mh,2));
    M(:,1)   = 1;
    M(:,end) = 1;
    M(1,:)   = 1;
    M(end,:) = 1;
end

%%%% H-Matrix (recursion)
if (Mh.typ == 0)
    A = hmxSpy(Mh.chd{1},[]);
    B = hmxSpy(Mh.chd{2},[]);
    C = hmxSpy(Mh.chd{3},[]);
    D = hmxSpy(Mh.chd{4},[]);
    M = [A,B;C,D];

%%% Compressed leaf
elseif (Mh.typ == 1)
    rk = size(Mh.dat{1},2);
    if (rk > 0)
        M(:,rk) = 1;
        M(rk,:) = 1;
    else
        M = 3 .* M;
    end

%%% Full leaf
elseif (Mh.typ == 2)
    if issparse(Mh.dat)
        M = 4 .* (M + (abs(Mh.dat)>0));
    else
        M = 2 .* M;
    end

%%% Unknown type
else
    error('hmxSpy.m : unavailable case')
end

% Graphical representation
if (length(varargin) == 1)
    spy(M==1,'b')
    hold on
    spy(M==2,'r')
    spy(M==3,'g')
    spy(M==4,'m')
    hold off
    title('openHmx : H-Matrix structure');
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function plotSigmoidalResponse_LGG(resp,Y,modelCI,nameOutcome,pathFig)
% -------------------------------------------------------------------------
% function plotSigmoidalResponse_LGG(resp,Y,modelCI,nameOutcome,pathFig)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function plots the sigmoidal response (P(y=1|x)) of a linear
% multivariable model response. See ref. [1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. resp: Column vector of size [nInst X 1] specifying the linear
%          multivariable model response, where 'nInst' refers to the number
%          of instances.
% 2. Y: Column vector of size [nInst X 1] specifying the outcome status
%       (1 or 0) for all instances.
% 3. modelCI: Column vector of size [nInst X 2] specifying the 95%
%             confidence interval on the multivariable model response as
%             defined by the 2.5 (modelCI(i,1)) and the 97.5 (modelCI(i,2))
%             percentiles, for the ith instance. To plot without CIs, use []
%             for the 3rd argument.See ref. [1] for more
%             details.
% 4. nameOutcome: String specifying the name of the
%                 modeled outcome.
%                 --> Ex: 'progression'
% 5. pathFig: (optional).  Full path to where figure is saved without
%             displaying it. Put '' for displaying the figure and not
%             saving it to 'pathFig' (default).
%             --> Ex: ''
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

if nargin < 5
    pathFig = '';
end

prob = 1./(1 + exp(-resp)); % Sigmoidal response, or P(y=1|x)
respPos = resp(Y==1);
respNeg = resp(Y==0);
probPos = prob(Y==1);
probNeg = prob(Y==0);
name = nameOutcome;
symbols = {'ob','xr'}; % First entry for positive instances, second entry for negative instances

if isempty(pathFig)
    figure
else
    h = figure('visible','off');
end
[sortResp,~] = sort(resp);
sigX = (sortResp(1)-5):0.1:sortResp(end)+5;
sigY = 1./(1 + exp(-sigX));
plot(sigX,sigY,'-k','LineWidth',3);
hold on
if ~isempty(modelCI)
    lowXpos = modelCI(Y==1,1); highXpos = modelCI(Y==1,2);
    lowXneg = modelCI(Y==0,1); highXneg = modelCI(Y==0,2);
    [~] = herrorbar(respPos,probPos,respPos-lowXpos,highXpos-respPos,'ob');
    hold on
    [~] = herrorbar(respNeg,probNeg,respNeg-lowXneg,highXneg-respNeg,'xr');
    legend('Sigmoidal response',['95% CI: ',name,' = 1'],['Status: ',name,' = 1'],['95% CI: ',name,' = 0'],['Status: ',name,' = 0'],'Location','NorthWest')
else
    plot(respPos,probPos,symbols{1},'LineWidth',6,'MarkerSize',18,'MarkerFaceColor',symbols{1}(end),'MarkerEdgeColor',symbols{1}(end));
    hold on
    plot(respNeg,probNeg,symbols{2},'LineWidth',6,'MarkerSize',20,'MarkerFaceColor',symbols{2}(end),'MarkerEdgeColor',symbols{2}(end));
    legend('Sigmoidal response',['Status: ',name,' = 1'],['Status: ',name,' = 0'],'Location','NorthWest')
end
xlabel('Multivariable model response','FontSize',30)
ylabel(['Probability that ',name,' = 1'],'FontSize',30)
axis([sortResp(1)-0.33*abs(sortResp(1)),sortResp(end)+0.33*abs(sortResp(end)),-0.05,1.05])
set(gca,'FontSize',24)

if ~isempty(pathFig)
    cd(pathFig)
    saveas(h,['sigPlotResponse_',name],'fig')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%kImportAny 'Generic Routine Interface (glue)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ImportAny.pane file
%
% Parameters:
% InputFile: i1 'Name of Input Image', required: 'input filename #1'
% OutputFile: o1 'Output Filename ', required: 'output filename #1'
%
% Example: o1 = kImportAny(i1, {'i1','';'o1',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ImportAny - Imports any Imageformat using convert
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kImportAny(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kImportAny(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o1', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'conviff"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ ret ] = setPin12Off( t_Kuka )
%% About:
% This function is used to set off the Pin 12 on the KUKA iiwa 7 R 800.

%% Syntax:
% [ ret ] = setPin12Off( t_Kuka )

%% Arreguments:
% t_Kuka: is the TCP/IP connection

%% Return value:
% ret: a boolean variable
% true: if the message is treated successfully
% false: if an error has happened

% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='pin12off';
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
ret=checkAcknowledgment(message);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function erase_and_display(old_str, new_str)
% Backspaces the number of chars of the length of old_str, then prints new_str

    fprintf([repmat('\b', 1, length(old_str)) new_str]);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef DataPuller < handle
%DataPuller Pull olympiad data from SAGE.
%   A DataPuller knows how, for a particular assay, to pull data from SAGE
%   based on a user-specified query, and merge/clean the results into a
%   single struct array. This is an assay-dependent operation; each assay
%   has differerent SAGE views and different requirements for merging and
%   cleaning.
%
%   DataPuller is meant to be used in conjunction with the
%   OlyDat.DataSelector.
%
%   See also OlyDat.DataSelector, OlyDat.Browser.

    properties (Abstract)
        % if true, pull from FS by default.
        tfPullFSDefault;

        % cellstr of options displayed in DataSelector pull-down and passed
        % to pullData function. can be empty.
        options;
    end

    methods (Abstract)
        %pullData/pullDataFS Pull Data from SAGE
        %   pullData is a SAGE data-pulling function for use with
        %   OlyDat.DataSelector. queryObjs is a cell array of query objects
        %   created by the user using the DataSelector. datapath is the
        %   user-specified folder containing experimental directories.
        %   option is a string, selected by the user in the DataSelector;
        %   it is one of the options specified in the DataPuller .options
        %   property.
        %
        %   data should be a struct array that results from pulling data
        %   from the appropriate dataSets and cleaning and merging that
        %   data. This struct array should have one element per experiment
        %   for use with the OlyDat.Browser.
        data = pullData(obj,queryObjs,option)
        data = pullDataFS(obj,queryObjs,datapath,option)
    end

    methods % hook methods

        function s = getCustomFieldPulldownInfo(obj) %#ok<MANU>
            % s = getCustomFieldPulldownInfo(obj)
            % Return info to resolve ambiguities when two or more datasets
            % share fields with the same name. Information returned by this
            % method is used by the DataSelector when populating pulldowns for
            % specifying the field in a query.
            %
            % This is a hook method intended to be overridden.
            %
            % s: struct array with fields
            %   * .Name: name of field
            %   * .DataSetName: name of DataSet whose DataSetField object
            %   should be used
            %   * .PrettyName: (optional). Currently not used/implemented.

            % Default implementation: none
            s = [];
        end

    end

    methods (Static) % utilities

        function fsp = reseatFileSystemPath(fsp,datapath)
            % fsp = reseatFileSystemPath(fsp,datapath)
            [~,fsp] = fileparts(fsp);
            fsp = fullfile(datapath,fsp);
        end

    end
end

author: kristinbranson

repo: JAABA


text:
function obj = getColorPrior(image,k)
    x = reshape(permute(image,[3 1 2]),[3 size(image,1) * size(image,2)]);
    vgg_mean =  [103.939, 116.779, 123.68];
    x = double(x([3 2 1],:));
    x = bsxfun(@minus,x,vgg_mean');
    options = statset('MaxIter',1000);
    while k>0
        try
            fprintf('color GMM, try clusters=%d...\n',k);
            obj = gmdistribution.fit(x',k,'Options',options);
            break;
        catch
            k = k-1;
        end;
    end;
end

author: happynear

repo: DeepVisualization


text:
function [ sicd_meta ] = meta2sicd_palsar2vol( native_meta )
%META2SICD_PALSAR2VOL Converts ALOS PALSAR 2 volume file into a SICD-style metadata structure
%
% Takes as input a metadata structure from read_ceos_vol_meta.
%
% Written by: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% CollectionInfo
if strncmpi(native_meta.vol_set_id, 'ALOS2', 5)
    sicd_meta.CollectionInfo.CollectorName = 'ALOS2';
    sicd_meta.CollectionInfo.RadarMode.ModeID=native_meta.text.prod_id(9:11);
    % Perhaps ScanSAR modes should be DYNAMIC STRIPMAP?
    if strncmpi('SBS',sicd_meta.CollectionInfo.RadarMode.ModeID,3)
        sicd_meta.CollectionInfo.RadarMode.ModeType='SPOTLIGHT';
    else
        sicd_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
    end
elseif strncmpi(native_meta.vol_set_id, 'STRIX', 5)
    sicd_meta.CollectionInfo.CollectorName = deblank(native_meta.vol_set_id);
    sicd_meta.CollectionInfo.RadarMode.ModeID=native_meta.text.prod_id(9:10);
    if strncmpi('SM',sicd_meta.CollectionInfo.RadarMode.ModeID,2)
        sicd_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
    elseif strncmpi('SL',sicd_meta.CollectionInfo.RadarMode.ModeID,2)
        sicd_meta.CollectionInfo.RadarMode.ModeType='DYNAMIC STRIPMAP';
    end
end
% Do we want to convert CoreName to NGA-style pattern?
sicd_meta.CollectionInfo.CoreName = native_meta.text.scene_id(8:end);
sicd_meta.CollectionInfo.CollectType='MONOSTATIC';
sicd_meta.CollectionInfo.Classification='UNCLASSIFIED';

%% ImageCreation
sicd_meta.ImageCreation.DateTime = ...
    datenum([native_meta.log_vol_create_date ...
             native_meta.log_vol_create_time(1:6)] ,'yyyymmddHHMMSS');
sicd_meta.ImageCreation.Profile='Prototype';

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint