text:
function copyaxes(source, destination, varargin)
%
% COPYAXES - copy a handle object axes inpt
%
% COPYAXES(SOURCE, DESTINATION) - copy axes from SOURCE to DESTINATION
%
% COPYAXES(..., isInSubplot) - if the destination is in a subplot figure
%       (default false).
%
% COPYAXES(..., isLegend) - if the axes is a legend (only if isInSubplot is
%       true) (default false)
%
%
% EXAMPLE: Copy a axes with plot
%
%   plot([1:0.1:2*pi], sin([1:0.1:2*pi]));
%   title('sin function')
%   xlabel('x')
%   ylabel('sin(x)')
%   ax = gca;
%
%   figure;
%   ax_new = axes;
%   copyaxes(ax, ax_new)
%
%
% EXAMPLE: Copy a axes with bar
%
%   bar(rand(10,5),'stacked');
%   title('bar stacket function')
%   xlabel('x label')
%   ylabel('y label')
%   ax = gca;
%
%   figure;
%   ax_new = axes;
%   copyaxes(ax, ax_new)
%
%
% EXAMPLE: Copy a axes and legend (the legend is an axes object)
%
%   plot([1:0.1:2*pi], sin([1:0.1:2*pi]));
%   title('sin function')
%   xlabel('x')
%   ylabel('sin(x)')
%   lg = legend('sin');
%   ax = gca;
%
%   figure;
%   ax_new = axes;
%   lg_new = axes;
%   copyaxes(ax, ax_new)
%   copyaxes(lg, lg_new)
%
% EXAMPLE: Copy an axes with a surface in a subplot.
%   Colormap is a figure property, not axes property.
%
%   k = 5;
%   n = 2^k-1;
%   [x,y,z] = sphere(n);
%   c = hadamard(2^k);
%   surf(x,y,z,c);
%   axis equal
%   title('sphere')
%   xlabel('x label')
%   ylabel('y label')
%   zlabel('z label')
%   ax = gca;
%   colormap([1  1  0; 0  1  1])
%
%   figure;
%   ax_new = subplot(2,2,1);
%   copyaxes(ax, ax_new, true)
%   colormap([1  1  0; 0  1  1])
%
% EXAMPLE:  Copy a surface in a subplot
%
%   [x,y] = meshgrid([-2:.4:2]);
%   Z = x.*exp(-x.^2-y.^2);
%   fh = figure('Position',[350 275 400 300],'Color','w');
%   ax = subplot(1,2,1);
%     set(ax, 'Color',[.8 .8 .8],'XTick',[-2 -1 0 1 2],...
%           'YTick',[-2 -1 0 1 2]);
%   sh = surface('XData',x,'YData',y,'ZData',Z,...
%              'FaceColor',get(ax,'Color')+.1,...
%              'EdgeColor','k','Marker','o',...
%              'MarkerFaceColor',[.5 1 .85]);
%   view(3)
%
%   ax_new = subplot(1,2,2);
%   copyaxes(ax, ax_new, true)
%
% EXAMPLE:  Copy an hist in a subplot
%
%     x = -4:0.1:4;
%     y = randn(500,1);
%     figure;
%     ax = subplot(1,2,1);
%     hist(y,x)
%
%     ax_new = subplot(1,2,2);
%     copyaxes(ax, ax_new, true)
%
%
% EXAMPLE:  Copy an pie in a subplot
%
%     x = [1 3 0.5 2.5 2];
%     explode = [0 1 0 0 0];
%     figure;
%     colormap jet
%     ax =  subplot(1,2,1);
%     pie(x,explode)
%
%     ax_new = subplot(1,2,2);
%     copyaxes(ax, ax_new, true)
%
%
% EXAMPLE:  Copy an figure in a subplot
%
%     load earth
%     figure
%     ax = subplot(1,2,1);
%     image(X); colormap(map)
%
%     ax_new = subplot(1,2,2);
%     copyaxes(ax, ax_new, true)
%
%
% tags: figure, copy, axes, line, pie, bar, hist, surface,
%
%
% author: Mar Callau-Zori
% PhD student - Universidad Politecnica de Madrid
%
% version: 1.2, December 2011
%

    isInSubplot = false;
    isLegend    = false;

    if nargin>2
        switch nargin

            case 3
                isInSubplot = varargin{1};
            case 4
                isInSubplot = varargin{1};
                isLegend    = varargin{2};
            otherwise
                error('ERROR:copyaxes:argChk', 'Wrong input arguments');
        end
    end

    if ~isLegend
        copyobj(get(source, 'Children'), destination);

        title = copyobj(get(source, 'Title'), destination);
        set(destination, 'Title', title);

        xlabel = copyobj(get(source, 'XLabel'), destination);
        set(destination, 'XLabel', xlabel);

        ylabel = copyobj(get(source, 'YLabel'), destination);
        set(destination, 'YLabel', ylabel);

        zlabel = copyobj(get(source, 'ZLabel'), destination);
        set(destination, 'ZLabel', zlabel);
    else
        set(destination, 'string', get(source, 'string'))
    end

    properties_str = {  'Units'
                        'ActivePositionProperty'
                        'ALim'
                        'ALimMode'
                        'AmbientLightColor'
                        'Box'
                        'CameraPosition'
                        'CameraPositionMode'
                        'CameraTarget'
                        'CameraTargetMode'
                        'CameraUpVector'
                        'CameraUpVectorMode'
                        'CameraViewAngle'
                        'CameraViewAngleMode'
                        'CLim'
                        'CLimMode'
                        'Color'
                        'CurrentPoint'
                        'ColorOrder'
                        'DataAspectRatio'
                        'DataAspectRatioMode'
                        'DrawMode'
                        'FontAngle'
                        'FontName'
                        'FontSize'
                        'FontUnits'
                        'FontWeight'
                        'GridLineStyle'
                        'Layer'
                        'LineStyleOrder'
                        'LineWidth'
                        'MinorGridLineStyle'
                        'NextPlot'
                        'OuterPosition'
                        'PlotBoxAspectRatio'
                        'PlotBoxAspectRatioMode'
                        'Projection'
                        'Position'
                        'TickLength'
                        'TickDir'
                        'TickDirMode'
                        'TightInset'
                        'View'
                        'XColor'
                        'XDir'
                        'XGrid'
                        'XAxisLocation'
                        'XLim'
                        'XLimMode'
                        'XMinorGrid'
                        'XMinorTick'
                        'XScale'
                        'XTick'
                        'XTickLabel'
                        'XTickLabelMode'
                        'XTickMode'
                        'YColor'
                        'YDir'
                        'YGrid'
                        'YAxisLocation'
                        'YLim'
                        'YLimMode'
                        'YMinorGrid'
                        'YMinorTick'
                        'YScale'
                        'YTick'
                        'YTickLabel'
                        'YTickLabelMode'
                        'YTickMode'
                        'ZColor'
                        'ZDir'
                        'ZGrid'
                        'ZLim'
                        'ZLimMode'
                        'ZMinorGrid'
                        'ZMinorTick'
                        'ZScale'
                        'ZTick'
                        'ZTickLabel'
                        'ZTickLabelMode'
                        'ZTickMode'
                        'BeingDeleted'
                        'ButtonDownFc'
                        'Clipping'
                        'CreateFcn'
                        'DeleteFcn'
                        'BusyAction'
                        'HandleVisibility'
                        'HitTest'
                        'Interruptible'
                        'Selected'
                        'SelectionHighlight'
                        'Tag'
                        'Type'
                        'UIContextMenu'
                        'UserData'
                        'Visible'};

    for i=1:length(properties_str)

        if (strcmpi(properties_str{i}, 'position') || strcmpi(properties_str{i}, 'outerposition'))

            if ~isInSubplot
                try
                    set(destination, properties_str{i}, get(source, properties_str{i}));
                catch e
                    if ~strcmpi(e.identifier, 'MATLAB:hg:propswch:FindObjFailed') && ...
                       ~strcmpi(e.identifier, 'MATLAB:hg:g_object:MustBeInSameFigure')
                            rethrow(e);
                    end
                end

            elseif isLegend
                % change the reference
                    try
                        units_dst = get(destination, 'units');
                        units_src = get(source, 'units');
                        set(destination, 'units', 'pixels');
                        set(source, 'units', 'pixels');
                        aux_dst = get(destination, properties_str{i});
                        aux_src = get(source, properties_str{i});
                        aux_dst(3) =  aux_src(3);
                        aux_dst(4) =  aux_src(4);
                        set(destination, properties_str{i}, aux_dst);
                        set(destination, 'units', units_dst);
                        set(source, 'units', units_src);
                    catch e
                        if ~strcmpi(e.identifier, 'MATLAB:hg:propswch:FindObjFailed') && ...
                           ~strcmpi(e.identifier, 'MATLAB:hg:g_object:MustBeInSameFigure')
                                rethrow(e);
                        end
                    end
            end
        else
            try
                set(destination, properties_str{i}, get(source, properties_str{i}));
            catch e
                if ~strcmpi(e.identifier, 'MATLAB:hg:propswch:FindObjFailed') && ...
                   ~strcmpi(e.identifier, 'MATLAB:hg:g_object:MustBeInSameFigure')
                         rethrow(e);
                end
            end
        end
    end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function upsilonjv = lfmjvComputeUpsilonMatrix(gamma, sigma2, t1, t2, mode)

% LFMJVCOMPUTEUPSILONMATRIX Upsilon matrix jolt. vel. with t1, t2 limits
% FORMAT
% DESC computes a portion of the LFMJV kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010
%
% SEEALSO : lfmComputeUpsilonMatrix.F, lfmvpComputeUpsilonMatrix.m

% KERN

sigma = sqrt(sigma2);
gridt1 = repmat(t1, 1, length(t2));
gridt2 = repmat(t2', length(t1), 1);
timeGrid = gridt1 - gridt2;

if mode==0
    upsilonjv = gamma^2*lfmvvComputeUpsilonMatrix(gamma, sigma2, t1, t2, mode) ...
        - (4/(sqrt(pi)*sigma^3))*exp(-(timeGrid.^2)./sigma2).* ...
        ((gamma + (2*timeGrid)/sigma2).*(1-(2*timeGrid.^2)/sigma2) + 4*timeGrid/sigma2);
else
    upsilonjv = gamma^2*lfmvvComputeUpsilonMatrix(gamma, sigma2, t1, t2, mode) ...
        - (4/(sqrt(pi)*sigma^3))*exp(-(timeGrid.^2)./sigma2).* ...
        ((gamma + (2*timeGrid)/sigma2).*(1-(2*timeGrid.^2)/sigma2) + 4*timeGrid/sigma2) ...
        + ((4*gamma)/(sqrt(pi)*sigma^3))*exp(-gamma*t1)*...
        ((t2.*(gamma - 2*t2/sigma2) + 1).*exp(-t2.^2/sigma2)).';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function g = whiteblockKernDiagGradient(kern, x, covDiag)

% WHITEBLOCKKERNDIAGGRADIENT WHITEBLOCK kernel's diagonal gradient wrt par.
% FORMAT
% DESC computes the gradient of the diagonal of the white noise kernel
% block matrix with respect to the parameters of the kernel. The
% parameters' gradients are returned in the order given by the
% whiteblockKernExtractParam command.
% ARG kern : the kernel structure for which the gradients are
% computed.
% ARG x : the input data for which the gradient is being computed.
% ARG covDiag : partial derivatives of the function of interest with
% respect to the diagonal elements of the kernel.
% RETURN g : gradients of the relevant function with respect to each
% of the parameters. Ordering should match the ordering given in
% whiteblockKernExtractParam.
%
% SEEALSO : whiteblockKernParamInit, kernDiagGradient
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

g = zeros(1, kern.nout);
startOne = 1;
endOne = 0;
for i=1:kern.nout
    endOne = endOne + size(x,1);
    g(1, i) = sum(covDiag(startOne:endOne));
    startOne = endOne + 1;
end


text:
function a = band(a,p,q)
%BAND         Extract band from matrix a, lower bandwidth p, upper bandwidth q
%   if parameter q is omitted, q:=p
%
%   c = band(a,p,q)
%

% written  09/28/01     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    improved performance
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin<3
    q = p;
  end

  a = tril(triu(a,-p),q);

repo: OpenMAS


text:
function y=rfft(x,n,d)
%RFFT     Calculate the DFT of real data Y=(X,N,D)
% Data is truncated/padded to length N if specified.
%   N even:	(N+2)/2 points are returned with
% 			the first and last being real
%   N odd:	(N+1)/2 points are returned with the
% 			first being real
% In all cases fix(1+N/2) points are returned
% D is the dimension along which to do the DFT

s=size(x);
if prod(s)==1
    y=x
else
    if nargin <3 || isempty(d)
        d=find(s>1);
        d=d(1);
        if nargin<2
            n=s(d);
        end
    end
    if isempty(n)
        n=s(d);
    end
    y=fft(x,n,d);
    y=reshape(y,prod(s(1:d-1)),n,prod(s(d+1:end)));
    s(d)=1+fix(n/2);
    y(:,s(d)+1:end,:)=[];
    y=reshape(y,s);
end

author: bastamon


text:
function [Estimation_X] = REKF_propagate(Estimation_X, OdometryFromThis2Next, Sigma_ODO )

% retrieve v and w
v = OdometryFromThis2Next(1:3);
w = OdometryFromThis2Next(4:6);

NumberOfLandmarks = size(Estimation_X.landmarks, 2);
%Jrw = jaco_r(-w);
Jrw=eye(3);
ExpMinusM = so3_exp(-w);

temp = repmat({Estimation_X.orientation}, NumberOfLandmarks+2,1 );
A = blkdiag(temp{:});
 A(4:6,1:3)=skew(Estimation_X.position)*Estimation_X.orientation;

if NumberOfLandmarks>0
   for i=1:NumberOfLandmarks
    A(6+3*i-2:6+3*i,1:3)=skew(Estimation_X.landmarks(1:3,i))*Estimation_X.orientation;
   end
end

%B1=[-Jrw  zeros(3,3); -skew(v)*Jrw -eye(3)];

B1=[eye(3)  zeros(3,3); zeros(3,3) eye(3)];

B=[B1; sparse(3*NumberOfLandmarks,6)];
B=sparse(B);

adA=A*B;

odoCov=diag([w.^2;v.^2])*Sigma_ODO^2;

% final update the covariance
Estimation_X.cov = Estimation_X.cov+ adA*odoCov*adA';

% update position and orientation
Estimation_X.position = Estimation_X.position+Estimation_X.orientation*v;
Estimation_X.orientation = Estimation_X.orientation*so3_exp(w);

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function [sigma_points, w_m, w_c] = compute_sigma_points(mu, sigma, lambda, alpha, beta)
% This function samples 2n+1 sigma points from the distribution given by mu and sigma
% according to the unscented transform, where n is the dimensionality of mu.
% Each column of sigma_points should represent one sigma point
% i.e. sigma_points has a dimensionality of nx2n+1.
% The corresponding weights w_m and w_c of the points are computed using lambda, alpha, and beta:
% w_m = [w_m_0, ..., w_m_2n], w_c = [w_c_0, ..., w_c_2n] (i.e. each of size 1x2n+1)
% They are later used to recover the mean and covariance respectively.

n = length(mu);
sigma_points = zeros(n,2*n+1);
w_m = zeros(1,2*n+1);
w_c = zeros(1,2*n+1);

% TODO: compute all sigma points
sigma_points(:,1) = mu;
addTerm = sqrtm((n+lambda)*sigma);
for i=1:n
    sigma_points(:,i+1) = mu + addTerm(:,i);
endfor
for i=n+1:2*n
    sigma_points(:,i+1) = mu - addTerm(:,i-n);
endfor

% TODO compute weight vectors w_m and w_c
w_m(1,1) = lambda/(n+lambda);
w_c(1,1) = w_m(1) + (1 - alpha*alpha + beta);
w_m(1,2:end) = 1/(2*(n+lambda));
w_c(1,2:end) = 1/(2*(n+lambda));
end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function [scd] = ft_scalpcurrentdensity(cfg, data)

% FT_SCALPCURRENTDENSITY computes an estimate of the SCD using the
% second-order derivative (the surface Laplacian) of the EEG potential
% distribution
%
% The relation between the surface Laplacian and the SCD is explained
% in more detail on http://tinyurl.com/ptovowl.
%
% Use as
%   [data] = ft_scalpcurrentdensity(cfg, data)
% or
%   [timelock] = ft_scalpcurrentdensity(cfg, timelock)
% where the input data is obtained from FT_PREPROCESSING or from
% FT_TIMELOCKANALYSIS. The output data has the same format as the input
% and can be used in combination with most other FieldTrip functions
% such as FT_FREQNALYSIS or FT_TOPOPLOTER.
%
% The configuration should contain
%   cfg.method       = 'finite' for finite-difference method or
%                      'spline' for spherical spline method
%                      'hjorth' for Hjorth approximation method
%   cfg.elec         = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.trials       = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.feedback     = string, 'no', 'text', 'textbar', 'gui' (default = 'text')
%
% The finite method require the following
%   cfg.conductivity = conductivity of the skin (default = 0.33 S/m)
%
% The spline and finite method require the following
%   cfg.conductivity = conductivity of the skin (default = 0.33 S/m)
%   cfg.lambda       = regularization parameter (default = 1e-05)
%   cfg.order        = order of the splines (default = 4)
%   cfg.degree       = degree of legendre polynomials (default for
%                       <=32 electrodes  = 9,
%                       <=64 electrodes  = 14,
%                       <=128 electrodes = 20,
%                       else             = 32
%
% The hjorth method requires the following
%   cfg.neighbours   = neighbourhood structure, see FT_PREPARE_NEIGHBOURS
%
% For the spline method you can specify the following
%   cfg.badchannel      = cell-array, see FT_CHANNELSELECTION for details (default = [])
%
% Note that the skin conductivity, electrode dimensions and the potential
% all have to be expressed in the same SI units, otherwise the units of
% the SCD values are not scaled correctly. The spatial distribution still
% will be correct.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% The 'finite' method implements
%   TF Oostendorp, A van Oosterom; The surface Laplacian of the potential:
%   theory and application. IEEE Trans Biomed Eng, 43(4): 394-405, 1996.
%   G Huiskamp; Difference formulas for the surface Laplacian on a
%   triangulated sphere. Journal of Computational Physics, 2(95): 477-496,
%   1991.
%
% The 'spline' method implements
%   F. Perrin, J. Pernier, O. Bertrand, and J. F. Echallier.
%   Spherical splines for scalp potential and curernt density mapping.
%   Electroencephalogr Clin Neurophysiol, 72:184-187, 1989
% including their corrections in
%   F. Perrin, J. Pernier, O. Bertrand, and J. F. Echallier.
%   Corrigenda: EEG 02274, Electroencephalography and Clinical
%   Neurophysiology 76:565.
%
% The 'hjorth' method implements
%   B. Hjort; An on-line transformation of EEG scalp potentials into
%   orthogonal source derivation. Electroencephalography and Clinical
%   Neurophysiology 39:526-530, 1975.
%
% See also FT_PREPROCESSING, FT_TIMELOCKANALYSIS, FT_FREQNALYSIS, FT_TOPOPLOTER.

% Copyright (C) 2004-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.method          = ft_getopt(cfg, 'method',       'spline');
cfg.conductivity    = ft_getopt(cfg, 'conductivity', 0.33); % in S/m
cfg.trials          = ft_getopt(cfg, 'trials',       'all', 1);
cfg.feedback        = ft_getopt(cfg, 'feedback',     'text');
cfg.badchannel      = ft_getopt(cfg, 'badchannel',     {});

switch cfg.method
  case 'hjorth'
    cfg = ft_checkconfig(cfg, 'required', {'neighbours'});
  case 'spline'
    cfg.lambda  = ft_getopt(cfg, 'lambda', 1e-5);
    cfg.order   = ft_getopt(cfg, 'order', 4);
    cfg.degree  = ft_getopt(cfg, 'degree', []);

    if isempty(cfg.degree) % determines degree of Legendre polynomials bases on number of electrodes
      nchan = numel(data.label);
      if nchan<=32
        cfg.degree = 9;
      elseif nchan<=64
        cfg.degree = 14;
      elseif nchan<=128
        cfg.degree = 20;
      else
        cfg.degree = 32;
      end
    end
  otherwise
    cfg = ft_checkconfig(cfg); % perform a simple consistency check
end

% store original datatype
dtype = ft_datatype(data);

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes', 'ismeg', []);

% get the electrode positions
tmpcfg = cfg;
tmpcfg.senstype = 'EEG';
elec = ft_fetch_sens(tmpcfg, data);

% select channels and trials of interest
tmpcfg = keepfields(cfg, {'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
tmpcfg.channel = elec.label;
data = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

Ntrials = numel(data.trial);

if isempty(cfg.badchannel)
  % check if the first sample of the first trial contains NaNs; if so treat it as a bad channel
  cfg.badchannel = ft_channelselection(find(isnan(data.trial{1}(:,1))), data.label);
end

% match the order of the data channels with the channel positions, order them according to the data
[datindx, elecindx] = match_str(data.label, elec.label);
[goodindx, tmp]     = match_str(data.label, setdiff(data.label, cfg.badchannel, 'stable'));

if ~isempty(cfg.badchannel)
  ft_info('detected channel %s as bad\n', cfg.badchannel{:});
  tmpcfg         = [];
  tmpcfg.channel = data.label(goodindx);
  data           = ft_selectdata(tmpcfg, data);
end

allchanpos  = elec.chanpos(elecindx,:); % the position of all channels, ordered according to the data
goodchanpos = allchanpos(goodindx,:);   % the position of good channels

% compute SCD for each trial
if strcmp(cfg.method, 'spline')
  fprintf('Checking spherical fit... ');
  [c, r] = fitsphere(allchanpos);
  d = allchanpos - repmat(c, size(allchanpos,1), 1);
  d = sqrt(sum(d.^2, 2));
  d = mean(abs(d) / r);
  if abs(d-1) > 0.1
    ft_warning('bad spherical fit (residual: %.2f%%). The interpolation will be inaccurate.', 100*(d-1));
  elseif abs(d-1) < 0.01
    fprintf('perfect spherical fit (residual: %.1f%%)\n', 100*(d-1));
  else
    fprintf('good spherical fit (residual: %.1f%%)\n', 100*(d-1));
  end
  % Builds the spatial filter only once.
  fprintf('Calculating the filter to build the SCD.\n');
  [WVo, WLo] = sphsplint(goodchanpos, allchanpos, cfg.order, cfg.degree, cfg.lambda);
  % Creates a montage to apply the spatial filter.
  montage.tra      = WLo;
  montage.labelold = elec.label(elecindx(goodindx));
  montage.labelnew = elec.label(elecindx);
  % Applies the montage to both the data and electrode definition
  scd  = ft_apply_montage(data, montage);
  elec = ft_apply_montage(elec, montage);

elseif strcmp(cfg.method, 'finite')
  if ~isempty(cfg.badchannel)
    ft_error('the method "%s" does not support the specification of bad channels', cfg.method);
  end
  % the finite difference approach requires a triangulation
  prj = elproj(allchanpos);
  tri = delaunay(prj(:,1), prj(:,2));
  % the new electrode montage only needs to be computed once for all trials
  montage.tra = lapcal(allchanpos, tri);
  montage.labelold = data.label;
  montage.labelnew = data.label;
  % apply the montage to the data, also update the electrode definition
  scd  = ft_apply_montage(data, montage);
  elec = ft_apply_montage(elec, montage);

elseif strcmp(cfg.method, 'hjorth')
  if ~isempty(cfg.badchannel)
    ft_error('the method "%s" does not support the specification of bad channels', cfg.method);
  end
  % convert the neighbourhood structure into a montage
  labelnew = {};
  labelold = {};
  for i=1:length(cfg.neighbours)
    labelnew = cat(2, labelnew, cfg.neighbours(i).label);
    labelold = cat(2, labelold, cfg.neighbours(i).neighblabel(:)');
  end
  labelold = cat(2, labelnew, labelold);
  labelold = unique(labelold);
  tra = zeros(length(labelnew), length(labelold));
  for i=1:length(cfg.neighbours)
    thischan   = match_str(labelold, cfg.neighbours(i).label);
    thisneighb = match_str(labelold, cfg.neighbours(i).neighblabel);
    tra(i, thischan) = 1;
    tra(i, thisneighb) = -1/length(thisneighb);
  end
  % combine it in a montage
  montage.tra = tra;
  montage.labelold = labelold;
  montage.labelnew = labelnew;
  % apply the montage to the data, also update the electrode definition
  scd  = ft_apply_montage(data, montage);
  elec = ft_apply_montage(elec, montage);

else
  ft_error('unknown method "%s"', cfg.method);
end

if strcmp(cfg.method, 'spline') || strcmp(cfg.method, 'finite')
  % correct the units
  ft_warning('trying to correct the units, assuming uV and mm');
  for trlop=1:Ntrials
    % The surface laplacian is proportional to potential divided by squared distance which means that, if
    % - input potential is in uV, which is 10^6 too large
    % - units of electrode positions are in mm, which is 10^3 too large
    % these two cancel out against each other. Hence the computed laplacian
    % is in SI units (MKS).
    scd.trial{trlop} = cfg.conductivity * -1 * scd.trial{trlop};
  end
  fprintf('output surface laplacian is in V/m^2\n');
else
  fprintf('output Hjorth filtered potential is in uV\n');
end

% Adds the electrode definition to the data.
scd.elec = elec;

% convert back to input type if necessary
switch dtype
  case 'timelock'
    scd = ft_checkdata(scd, 'datatype', 'timelock');
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data

% rename the output variable to accomodate the savevar postamble
data = scd;

ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [psd, ft] = gsp_jtv_estimate_psd(G, X, param)
%GSP_JTV_ESTIMATE_PSD Estimate the PSD of a time vertex process
%   Usage:  psd = gsp_jtv_estimate_psd(G, X)
%           psd = gsp_jtv_estimate_psd(G, X, param)
%
%   Input parameters:
%         G          : Graph
%         X          : Signal(s) (matrix of NxTxR)
%         param      : Optional parameters
%   Output parameters:
%         psd        : PSD matrix
%         ft         : Filter type
%
%   This function estimates the PSD for vertex time processes. It is a
%   generalization of the Bartlett method.
%
%   Additional parameters
%   ---------------------
%
%   * *param.estimator*  : Method used for the estimation. By default 'TA'
%     - 'TA': Average over the time domain.
%     - 'FC': Fourier convolution a trick used to reduce the computation.
%     - 'VA': Average over the vertex domain
%     - 'TVA': Average over the time and the vertex domain
%     - 'FA' : Average on multiple realization in the joint spectral domain
%   * *param.L*       : Length of the window for the 'GSTFT' method
%   * *param.kernel*  : Kernel for the convolution for 'FC' method.
%
%
%   References: perraudin2016stationary

% Author : Nathanael Perraudin
% Date: 6 January 2016
% Testing: test_gsp_estimate_vertex_time_psd

if nargin<3
    param = struct;
end
if ~isfield(param,'estimator')
    param.estimator = 'FA';
end

% number of realizations
[N, T, R] = size(X);

switch param.estimator

    % This method compute the Fourier transform of the signals and average
    % the modulus
    case 'FA'
        if ~isfield(param, 'L')
            param.L = compute_default_L(X,T,R);
        end
        Xhat = gsp_jft(G,reshape(X,N,param.L,R*T/param.L));
        psd = mean(abs(Xhat).^2,3);
        ft = 'js';

    % This method splits time in windows, and takes an average of the psd
    % over each time window. The size of each window (L) should be
    % equal to the maximum correlation distance in the process
    case 'TA'

        if ~gsp_check_fourier(G)
            error('The Fourier basis is needed the Fourier basis')
            %G = gsp_compute_fourier_basis(G);
        end

        if ~isfield(param, 'L')
            param.L = compute_default_L(X,T,R);
        end
        if ~isfield(param, 'a'),       param.a = round(param.L/2); end
        if ~isfield(param, 'M'),       param.M = G.jtv.T; end
        if ~isfield(param,'win_type'), param.win_type = 'itersine'; end

        if round(T/param.L) - (T/param.L) ~= 0
            warning('Optimally the length of the signal should be divisible by param.L.')
            while round(T/param.L) - (T/param.L) ~= 0
                param.L = param.L + 1;
            end
        end

        % Compute a Gabor window
        w = gabwin(param.win_type, param.a, param.L);

        % estimate the psd as the average over each realization
        psd = zeros(G.N, G.jtv.T);
        for r = 1:R

            % do a windowed JFT
            coeff = gsp_jtwgft(G, w, X(:,:,r), param);

            % compute the psd for each window
            psd_win = abs(coeff(:, 1:end-1, :)).^2;

            % average all windows to obtain an estimate of the psd
            psd_est = squeeze(mean(psd_win, 2));

            % average over realizations
            psd = psd + transpose(psd_est);
        end

        psd = psd / R;

        % normalize energy
        psd = psd * (norm(X(:), 'fro')^2/R / sum(psd(:)));

    % Go to the joint frequency domain and smooth out the PSD with a kernel.
        ft = 'js';

    case 'FC'
        if ~isfield(param, 'kernel')
            param.kernel = exp(-(-20:20).^2/3);
        end
        h = param.kernel;

        if ~gsp_check_fourier(G)
            G = gsp_compute_fourier_basis(G);
        end

        Xhat2 = abs(gsp_jft(G,X)).^2;

        psd = conv2(h', h', mean(Xhat2,3), 'same') / norm(h, 1)^2;
%         psd = conv2(mean(Xhat2,3), param.kernel, 'same') / norm(param.kernel, 1);

        % normalize energy
        psd = psd * (norm(X(:), 'fro')^2/R / sum(psd(:)));
        ft = 'js';

    % For scalability,
    case 'VA'

        if ~isfield(G,'boundary'); param.boundary = 'periodic'; end

        switch param.boundary

            case 'periodic'
                Xhat = fft(X,[],2)/sqrt(T);

            case 'reflecting'
                error('Sorry. Not implemented yet...')

            otherwise
                error('Unknown boundary condition');
        end

        if ~isfield(G,'lmax')
            G = gsp_estimate_lmax(G);

            warning(['GSP_PSD_ESTIMATION: The variable lmax is not ',...
                'available. The function will compute it for you. ',...
                'However, if you apply many time this function, you ',...
                'should precompute it using the function: ',...
                'gsp_estimate_lmax']);
        end

        if ~isfield(param,'Nfilt'),   param.Nfilt = 50; end
        if ~isfield(param,'Nrandom'), param.Nrandom = max(10, R); end
        if ~isfield(param,'g0')
            sigma = sqrt(2*G.lmax/param.Nfilt^2 * (param.Nfilt + 1));
            param.g0 = @(x) exp(-x.^2/sigma^2);
        end

        % Design the frame
        [ g , mu ] = gsp_design_translates(G, param.g0, param.Nfilt);
        %[ g , mu ] = gsp_design_itersine(G,Nfilt );

        % Estimate the energy of the window
        if gsp_check_fourier(G)
            mu_y2 = sum(abs(gsp_filter_evaluate(g,G.e)).^2,1)';
        else
            w = randn(N, param.Nrandom);
            x_filt2    = gsp_vec2mat( gsp_filter(G,g,w,param), param.Nfilt);
            n2_x_filt2 = sum(abs(x_filt2).^2,1);
            mu_y2      = reshape(mean(n2_x_filt2,3),[],1);
        end

        psd = cell(1,T);

        for ii = 1:T
            % Perform the filtering
            x_filt = gsp_vec2mat( ...
                gsp_filter(G, g, squeeze(Xhat(:,ii,:)), param), ...
                param.Nfilt);
            % estimate the points
            n2_x_filt = sum(abs(x_filt).^2,1);
            mu_y = reshape(mean(n2_x_filt,3),[],1);

            % Interpolate to obtain a nice filter.
            psd{ii} = @(s) max(spline(mu,mu_y./mu_y2,s), 0);
        end
        ft = 'js-array';

   case 'TVA'

        if ~isfield(param, 'L')

            threshold = 0.05;

            % compute the average correlation distance
            C_T = zeros(T); Toep = toeplitz(0:(T-1)); Tmax = round(0.5*T);

            for r = 1:R, C_T = C_T + abs(X(:,:,r)'*X(:,:,r)) / R; end
            cor = zeros(Tmax,1);
            for t = 1:Tmax, cor(t) = mean( vec(C_T(Toep == t-1)) ); end
            cor = normalize_data(cor, 0, 1);

            % select the window length conservatively
            param.L = min(4*find(cor>=threshold, 1, 'last' ), T);

            % visualize the correlation
            % figure; plot(0:(Tmax-1), cor, '-o', [1 1]*param.L, [0 1],
            % 'r-'); xlabel('correlation distance');

            % param.L = round(T/4);
        end

       % Make sure that the length of the signal should be divisible by param.L
       if round(T/param.L) - (T/param.L) ~= 0
           % warning('Ideally the length of the signal should be divisible by param.L.')
           while (round(T/param.L) - (T/param.L) ~= 0) || (round(param.L/2) - (param.L/2) ~= 0)
              param.L = param.L + 1;
           end
       end

       if ~isfield(param, 'a'),       param.a = round(param.L/2);  end
       if ~isfield(param, 'M'),       param.M = G.jtv.T;           end
       if ~isfield(param,'win_type'), param.win_type = 'itersine'; end
       if ~isfield(G,'boundary');     param.boundary = 'periodic'; end

       % Compute a Gabor window
       w = gabwin(param.win_type, param.a, param.L);

        switch param.boundary

            case 'periodic'
                S = zeros(N, param.M, ceil(T/param.a), R);
                for r = 1:R
                    % do a discrete gabore transform
                    X_gabore = dgt(transpose(X(:,:,r)), w, param.a, param.M);
                    S(:,:,:,r) = permute( X_gabore, [3,1,2]);
                end
            case 'reflecting'
                error('Sorry. Not implemented yet...')
            otherwise
                error('Unknown boundary condition');
        end

        if ~isfield(G,'lmax')
            G = gsp_estimate_lmax(G);

            warning(['GSP_PSD_ESTIMATION: The variable lmax is not ',...
                'available. The function will compute it for you. ',...
                'However, if you apply many time this function, you ',...
                'should precompute it using the function: ',...
                'gsp_estimate_lmax']);
        end

        if ~isfield(param,'Nfilt'),   param.Nfilt = 50;          end
        if ~isfield(param,'Nrandom'), param.Nrandom = max(10,R); end
        if ~isfield(param,'g0')
            sigma = sqrt(2*G.lmax/param.Nfilt^2 * (param.Nfilt + 1));
            param.g0 = @(x) exp(-x.^2/sigma^2);
        end

        % Design the frame
        [ g , mu ] = gsp_design_translates(G, param.g0, param.Nfilt);
        %[ g , mu ] = gsp_design_itersine(G,Nfilt );

        % Estimate the energy of the window
        if gsp_check_fourier(G)
            mu_y2 = sum(abs(gsp_filter_evaluate(g, G.e)).^2, 1)';
        else
            w = randn(N, param.Nrandom);
            x_filt2 = gsp_vec2mat( gsp_filter_analysis(G,g,w,param), param.Nfilt);
            n2_x_filt2 = sum(abs(x_filt2).^2, 1);
            mu_y2 = reshape(mean(n2_x_filt2,3), [], 1);
        end

        psd = cell(1, param.M);

%         if gsp_check_fourier(G) && param.use_fast
%             warning('This may use a lot of ram!')
%             [N, T, Nt, Ns ] = size(S);
%             Nf = numel(g);
%             ShatG = gsp_gft(G,S);
%             filt_eval = gsp_filter_evaluate(g, G.e);
%             filt_eval = repmat(reshape(filt_eval,[N,1,1,1,Nf]), 1,T,Nt,Ns ,1);
%             % Perform the filtering
%             S_filt = gsp_igft(G, repmat(ShatG,[1,1,1,1,Nf]) .* filt_eval);
%             n2_S_filt = squeeze(mean(mean(sum(abs(S_filt).^2,1),3),4));
%             for ii = 1:param.M
%
%                 % Interpolate to obtain a nice filter.
%                 psd{ii} = @(s) max(spline(mu,reshape(n2_S_filt(ii,:),[],1)./mu_y2,s), 0);
%             end
%         else
        for ii = 1:param.M

            % Perform the filtering
            x_filt = gsp_vec2mat( ...
                gsp_filter(G, g, reshape(S(:,ii,:,:),N,[]), param), ...
                param.Nfilt);

            % estimate the points
            n2_x_filt = sum(abs(x_filt).^2,1);
            mu_y = reshape(mean(n2_x_filt,3),[],1);

            % Interpolate to obtain a nice filter.
            psd{ii} = @(s) max(spline(mu,mu_y./mu_y2,s), 0);
        end
%         end

        ft = 'js-array';

    otherwise
        error('Unknown method')
end

end

function [x] = normalize_data(x, xmin, xmax)
%NORMALIZE_VECTOR Normalize x between xmin and xmax
%   x might be a scalar, vector, or matrix.

%normalize to [0,1]
x = (x - min(min(x))) ./ (max(max(x)) - min(min(x)));

if exist('xmax', 'var') && exist('xmin', 'var')
   x = x .* (xmax - xmin) + xmin;
end

end

function L =  compute_default_L(X,T,R)
% compute the correlation distance
C_T = zeros(T);
for r = 1:R
    C_T = C_T + abs(X(:,:,r)'*X(:,:,r)) / R;
end
Toep = toeplitz(0:(T-1));
Tmax = round(0.5*T);
cor = zeros(Tmax,1);
for t = 1:Tmax
    cor(t) = mean( vec(C_T(Toep == t-1)) );
end
cor = normalize_data(cor, 0, 1);
threshold = 0.05;
L = min(2*find(cor>=threshold, 1, 'last'), T);

%L = round(T/4);
%figure; plot(0:(Tmax-1), cor, '-o', [1 1]*param.L, [0 1], 'r-');

end

repo: gspbox


text:
function test_failed = test_libltfat_ifftreal(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
phaseconv = enuminfo.ltfat_phaseconvention;

fftwflags = struct('FFTW_MEASURE',0,'FFTW_ESTIMATE',64,'FFTW_PATIENT',32,'FFTW_DESTROY_INPUT',1,...
    'FFTW_UNALIGNED',2,'FFTW_EXHAUSTIVE',8,'FFTW_PRESERVE_INPUT',16);

Larr  = [351];
Warr  = [6];

for idx = 1:numel(Larr)
    L = Larr(idx);
    W = Warr(idx);
    M2 = floor(L/2) + 1;

    f = randn(M2,W,flags.complexity) + 1i*randn(M2,W,flags.complexity);
    fi = complex2interleaved(f);
    fPtr = libpointer(dataPtr,fi);

    c = cast(randn(L,W),flags.complexity);
    coutPtr = libpointer(dataPtr,c);

    truec = ifftreal(f,L)*L;

    funname = makelibraryname('ifftreal',flags.complexity,0);
    status = calllib('libltfat',funname,fPtr,L,W,coutPtr);

    res = norm(truec - coutPtr.Value,'fro');
    [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
    fprintf(['FFT   L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    % With plan
    c = cast(randn(L,W),flags.complexity);
    coutPtr = libpointer(dataPtr,c);

    plan = libpointer();
    funname = makelibraryname('ifftreal_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,L,W,fPtr,coutPtr, fftwflags.FFTW_ESTIMATE, plan);

    funname = makelibraryname('ifftreal_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan);

    res = norm(truec - coutPtr.Value,'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    c = cast(randn(L,W),flags.complexity);
    coutPtr = libpointer(dataPtr,c);

    funname = makelibraryname('ifftreal_execute_newarray',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan,fPtr,coutPtr);

    res = norm(truec - coutPtr.Value,'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    funname = makelibraryname('ifftreal_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    %%%%%% Inplace
    c = f;
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('ifftreal_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,L,W, coutPtr, coutPtr, fftwflags.FFTW_MEASURE, plan);

    funname = makelibraryname('ifftreal_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan);

    ctmp = coutPtr.Value;
    res = norm(truec - ctmp(1:L,:),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    c = f;
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    funname = makelibraryname('ifftreal_execute_newarray',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan,coutPtr,coutPtr);

    ctmp = (coutPtr.Value);
    res = norm(truec - ctmp(1:L,:),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    funname = makelibraryname('ifftreal_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function marginal = marginal_nodes(engine, nodes, t, add_ev)
% MARGINAL_NODES Compute the marginal on the specified query nodes (loopy_unrolled_dbn)
% marginal = marginal_nodes(engine, nodes, t)
%
% 't' specifies the time slice of the earliest node in 'nodes'.
% 'nodes' must occur in some clique.
%
% Example:
% Consider a DBN with 2 nodes per slice.
% Then t=2, nodes=[1 3] refers to node 1 in slice 2 and node 1 in slice 3,
% i.e., nodes 3 and 5 in the unrolled network,

if nargin < 3, t = 1; end
if nargin < 4, add_ev = 0; end

marginal = marginal_nodes(engine.unrolled_engine, nodes + (t-1)*engine.ss, add_ev);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [] = draw_funciton(f, x_range_min, x_range_max, y_range_min, y_range_max, unit_len)

    xCoarse = x_range_min:unit_len:x_range_max;
    yCoarse = y_range_min:unit_len:y_range_max;
    [XX,YY] = meshgrid(xCoarse,yCoarse);
    row_size = size(XX,1);
    col_size = size(XX,2);
    for j=1:col_size
        for i=1:row_size
           w = [XX(i,j); YY(i,j)];
           ZZ(i,j) = f(w);
        end
    end

    figure

    % draw contour
    contour(XX,YY,ZZ,50)

    % draw x-axis
    x = [x_range_min x_range_max];
    y = [0 0];
    line(x,y,'Color','black','LineStyle','-')

    % draw y-axis
    x = [0 0];
    y = [y_range_min y_range_max];
    line(x,y,'Color','black','LineStyle','-')

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function output = callmosek(model)

% Retrieve needed data
options = model.options;
F_struc = model.F_struc;
c       = model.c;
Q       = model.Q;
K       = model.K;
x0      = model.x0;
integer_variables = model.integer_variables;
binary_variables = model.binary_variables;
extended_variables = model.extended_variables;
ub      = model.ub;
lb      = model.lb;
mt      = model.monomtable;

% *********************************
% What type of variables do we have
% *********************************
model.linear_variables = find((sum(abs(mt),2)==1) & (any(mt==1,2)));
model.nonlinear_variables = setdiff((1:size(mt,1))',model.linear_variables);
model.sigmonial_variables = find(any(0>mt,2) | any(mt-fix(mt),2));

% Some meta-solver thinks we handle binaries
if ~isempty(model.binary_variables)
    integer_variables = union(model.integer_variables, model.binary_variables);
    if isempty(lb)
        lb = repmat(-inf,1,length(model.c));
    end
    lb(model.binary_variables) = max(lb(model.binary_variables),0);
    if isempty(ub)
        ub = repmat(inf,1,length(model.c));
    end
    ub(model.binary_variables) = min(ub(model.binary_variables),1);
    model.lb = lb;
    model.ub = ub;
    model.integer_variables = integer_variables;
end

% Some meta solvers might construct model with empty cones
if any(model.K.s) && any(model.K.s == 0)
    model.K.s(model.K.s==0)=[];
end
if any(model.K.q) && any(model.K.q == 0)
    model.K.q(model.K.q==0)=[];
end

if ~isempty(model.sigmonial_variables) || isequal(model.solver.version,'GEOMETRIC')
    [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_geometric(model);
else
    [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_lpqpsocpsdp(model);
end

% YALMIP has introduced internal variables for socp/exp cones etc
if length(x) > 0 && length(x) ~= length(model.c)
    x = x(1:length(model.c));
end

% Save all data sent to solver?
if options.savesolverinput
    solverinput.prob = prob;
    solverinput.param = options.mosek;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.r = r;
    solveroutput.res = res;

else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,model.solver.tag,solverinput,solveroutput,solvertime);

function [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_lpqpsocpsdp(model);

[model,output] = normalizeExponentialCone(model);
if output.problem == 23
    % Cases such as 1-exp(t) >= x'*x
    [model,problem] = splitCones(model);
    if problem
        problem = 19;
        x = [];
        D_struc = [];
        r = [];
        res = [];
        solvertime = 0;
        prob = [];
    else
        [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_lpqpsocpsdp(model);
    end
elseif output.problem
    problem = output.problem;
    x = [];
    D_struc = [];
    r = [];
    res = [];
    solvertime = 0;
    prob = [];
else
    if isempty(model.integer_variables)
        % Mosek cannot exploit this for continuous models (and if we keep
        % it from e.g. bnb, then the prmal version will be called, which
        % will run much slower for large SDPs)
        model.x0 = [];
    end
    if nnz(model.Q)==0 && isempty(model.integer_variables) && isempty(model.x0)
        % Standard cone problem which we can model by sending our standard dual
        % and then recover solution via Moseks dual
        [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_dual(model);
    else
        % Integer conic program
        % Quadratic objective
        % Exponential cones
        [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_primal(model);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function view = er_voxRMap(view, roi, scans, dt, conditions, subsets);
%
% view = er_voxRMap(<view=cur view>,  <roi=selected>, <scans=scan group>, ...
%                          <dt=cur dt>, <conditions=all>, <subsets=even odd>);
%
% Export a map of "voxel reliability": for each voxel, the correlation
% coefficient R between response amplitudes across conditions for two
% subsets. By default, these subsets are even and odd runs in the scan
% group.
%
% If 'gray' is provided as the ROI, will step through gray matter computing
% the map in a memory-efficient manner.
%
% subsets: 2x1 cell, e.g. {[1 3 5 7] [2 4 6 8]}, showing which subsets of
% runs to compare. By default, will compare even and odd runs across the
% selected scans.
%
% ras, 05/01/06.
if notDefined('view'), view = getCurView;               end
if notDefined('roi'),  roi = view.selectedROI;          end
if notDefined('scans'), scans = er_getScanGroup(view);  end
if notDefined('dt'), dt = view.curDataType;              end
if notDefined('conditions'),
    trials = er_concatParfiles(view);
    conditions = trials.condNums(trials.condNums>0);
end
if notDefined('subsets'), % {odd even}
    subsets = {scans(1:2:end) scans(2:2:end)};
end

if ~isequal(roi, 'gray')
    % for an ROI: call mv_exportMap
    mv = mv_init(view, roi, scans, dt);
    mv_exportMap(mv, 'voxel reliability', 1);
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If we got here, we're doing a memory-efficient gray matter computation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roi = getGrayRoi(view);
roi.coords = roiSubCoords(view, roi.coords); % remove redundant voxels
nVoxels = size(roi.coords, 2);

hwait = mrvWaitbar(0, 'Computing Voxel Reliability Across Gray Matter...');
set(hwait, 'Position', get(hwait,'Position')+[0 100 0 0]);
for a = 1:2000:nVoxels
    b = min(nVoxels, a+1999);
    subRoi = roi;
    subRoi.coords = roi.coords(:,a:b);
    subRoi.name = 'ROI1'; % will prevent lengthy caching of each sub-ROI

    mv = mv_init(view, subRoi);

    mv.params.selConds = conditions;
    mv = mv_reliability(mv, 'plotFlag', 0);

    voxR(a:b) = mv.wta.voxR;

    mrvWaitbar(b/nVoxels, hwait);
end
close(hwait);

% plug in the values to the map volume:
mapdims = viewGet(view,'dataSize');
nScans = viewGet(view,'numScans');
mapvol = zeros(mapdims);
scan = scans(1);
ind = roiIndices(view, roi.coords);
mapvol(ind) = voxR;

%%%%%set in view
map = cell(1,nScans);
mapName = sprintf('VoxR_%s', roi.name);
map{scan} = mapvol;
view = setParameterMap(view, map, mapName);
saveParameterMap(view, [], 1);
disp('Finished computing gray matter voxel reliability map.')

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function roi = getGrayRoi(view);
% roi = getGrayRoi(view);
% find, load, or make a gray ROI for the view.

% try finding a loaded ROI
if ~isempty(view.ROIs)
    existingRois = {view.ROIs.name};
    N = cellfind(existingRois, 'gray');
    if ~isempty(N), roi = view.ROIs(N); return; end
end

% try loading a saved ROI
w = what(roiDir(view));
savedRois = w.mat';
if ~isempty(savedRois)
    N = cellfind(savedRois, 'gray.mat');
    if ~isempty(N)
        load(fullfile(roiDir(view), 'gray.mat'), 'ROI');
        roi = ROI;
        return
    end
end

% got here: need to make it
try
    view = makeGrayROI(view);
    roi = view.ROIs(end);
    saveROI(view, roi);
catch
    disp('Couldn''t make a gray ROI. Segmentation installed?')
    roi = [];
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function y = geo_mean( x, dim, w )

%GEO_MEAN   Geometric mean.
%   Y=GEO_MEAN(X), where X is a vector, computes the geometrix mean of X. If any
%   of the elements of X are negative, then Y=-Inf. Otherwise, it is equivalent
%   to Y=PROD(X).^(1/LENGTH(X)). All elements must be real.
%
%   For matrices, GEO_MEAN(X) is a row vector containing the geometric means of
%   the columns. For N-D arrays, GEO_MEAN(X) is an array of the geometric means
%   taken along the first non-singleton dimension of X.
%
%   GEO_MEAN(X,DIM) takes the geometric mean along the dimension DIM of X.
%
%   GEO_MEAN(X,DIM,W), where W is a vector of nonnegative integers, computes a
%   weighted geometric mean Y = PROD(X.^W)^(1/SUM(W)). This is more efficient
%   than replicating the values of X W times. Note that W must be a vector,
%   even if X is a matrix, and its length must be the same as SIZE(X,DIM).
%
%   Disciplined convex programming information:
%       GEO_MEAN is concave  and nondecreasing; therefore, when used in CVX
%       specifications, its argument must be concave.

%
% Check arguments
%

error( nargchk( 1, 3, nargin ) );
if ~isreal( x ),
    error( 'First argument must be real.' );
elseif nargin < 2,
    dim = cvx_default_dimension( size( x ) );
elseif ~cvx_check_dimension( dim ),
    error( 'Second argument must be a positive integer.' );
end
sx = size( x );
nx = sx( dim );

%
% Third argument check
%

if nargin < 3 || isempty( w ),
    w = [];
elseif numel( w ) ~= length( w ) || ~isnumeric( w ) || ~isreal( w ) || any( w < 0 ) || any( w ~= floor( w ) ),
    error( 'Third argument must be a vector of nonnegative integers.' );
elseif length( w ) ~= nx,
    error( 'Third argument must be a vector of length %d.', nx );
else
    w = reshape( w, 1, nx );
end

if nx == 0,
    sx( dim ) = 1;
    y = ones( sx );
else
    if nx == 1,
        y = x;
    elseif isempty( w ) || ~any( diff( w ) ),
        y = exp( sum( log( max( x, realmin ) ), dim ) * ( 1 / nx ) );
    elseif dim == 1,
        y = exp( w * log( max( x, realmin ) ) * ( 1 / sum( w ) ) );
    else
        pvec = [ dim, 1 : dim - 1, dim + 1 : ndims( x ) ];
        y = ipermute( exp( w * log( max( permute( x, pvec ), realmin ) ) * ( 1 / sum( w ) ) ), pvec );
    end
    xmin = min( x, [], dim );
    y( xmin <  0 ) = -Inf;
    y( xmin == 0 ) = 0;
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
clear; close all; clc; dbstop if error;

%--------------------------------------------------------------------------
% This demo process the data from the RIT Target Detection Blind Test
% contest which is located at: http://dirsapps.cis.rit.edu/blindtest/
% To use this file, you select a target detection algorithm and a target to
% find and then the script runs the algorithm and outputs the data into the
% outputDir.  Two files are outputted, a .img and a .hdr. You upload these
% files to the RIT website and they are automatically scored.
%--------------------------------------------------------------------------
% Parameters
inputFilename = 'data\blind_test\HyMap\blind_test_refl.img';
fasticaToolboxPath = '..\matlab_hyperspectral_toolbox\trunk\FastICA_25';
targetFilenames = {'data\blind_test\SPL\F5\F5_f.txt'};
outputDir = 'RIT Data Results';
% See switch statement for algorithm choices
algorithm = 'ace'
%algorithm = 'rmf-sum';
%algorithm = 'plmf'
%algorithm = 'matchedFilter';
%algorithm = 'sam'
%--------------------------------------------------------------------------

addpath('gmm');

addpath(fasticaToolboxPath);
mkdir(outputDir);

% Read in the data
w = 280;
h = 800;
p = 126;
M = multibandread(inputFilename, [w h p], 'int16', 0, 'bil', 'ieee-le')/1e4;
lData = hyperGetHymapWavelengthsNm();

% Read in target signatures
[sig1, lSig] = hyperGetEnviSignature(targetFilenames{1});

% Get signature from data for comparison
fsig1 = squeeze(M(122,495,:));
%sig1 = fsig1;

% Resample data to commone wavelength set
desiredLambdas = lData;
sig1 = squeeze(hyperResample(sig1, lSig, desiredLambdas));
figure; plot(sig1); grid on; title('Signature 1');
    xlabel('Wavelength [nm]'); ylabel('Reflectance [%]');
    hold on; plot(fsig1, '--');
    legend('Recorded', 'From Image');

goodBands = 1:p; %[3:63 69:93 98:123];

% Image sharpening
if 0
    ff = fspecial('unsharp',0.2);
    for k=1:p
        M(:,:,k) = imfilter(M(:,:,k),ff,'same');
        M(:,:,k) = imfilter(M(:,:,k),ff,'same');
        %M(:,:,k) = imfilter(M(:,:,k),ff,'same');
    end
end

figure; imagesc(M(:,:,40)); axis image; colormap(gray);

% Try to discover in-situ to lab kernel.
% TODO
% sub(:,1) = M(144,515,:);
% sub(:,2) = M(144,516,:);
% sub(:,3) = M(144,517,:);
% sub(:,4) = M(145,515,:);
% sub(:,5) = M(145,516,:);
% sub(:,6) = M(145,517,:);
% sub(:,7) = M(146,515,:);
% sub(:,8) = M(146,516,:);
% sub(:,9) = M(146,517,:);
%
% alpha = pinv(sub)*sig1; %alpha = alpha ./ sum(alpha(:));
% err = sub*alpha - sig1; err = err - mean(err); badBands = find(abs(err)>0.02);
% goodBands = setxor(1:p,badBands);
% figure; plot(err); hold on; plot(sig1,'.'); plot(fsig1,'.-'); hold off; grid on;
%     legend({'err','lab sig','in situ sig'})
% alpha = reshape(alpha,3,3);
% figure; imagesc(alpha);
%
% for k=1:p
%     %M(:,:,p) = conv2(M(:,:,p),alpha,'same');
% end

% Emperical dervied
goodBands = [3     4     5     6     7     8     9    10    11    12    14    15    16    17    18 ...
    19    20    22    23    24    26    28    29    31    32    33    34    35    36    37 ...
    38    39    40    41    42    43    44    45    46    49    51    52    53    54    55 ...
    56    57    58    59    60    61    62    66    69    70    71    72    86    87    88 ...
    89    90    91    92    93    96    97    98    99   100   101   102   103   104   105 ...
   106   107   108   109   110   111   112   113   115   116   117   119   120   121   122];
%sig1 = squeeze(hyperResample(sig1, lSig, desiredLambdas));
figure; plot(sig1(goodBands)); grid on; title('Signature 1 - good bands only');
    xlabel('Wavelength [nm]'); ylabel('Reflectance [%]');
    hold on; plot(fsig1(goodBands), '--');
    legend('Recorded', 'From Image');

% Display data
M = hyperConvert2d(M);
%[M, H, snr] = hyperMnf(M, w, h);
M_pct = hyperPct(M, 3);
M_pct = hyperNormalize(hyperConvert3d(M_pct,w,h,3));
figure; imagesc(M_pct); axis image; title('Scene');

% Data conditioning
M = M(goodBands, :);
sig1 = sig1(goodBands);
%fsig1 = fsig1(goodBands);
%sig1 = fsig1;

%q = hyperHfcVd(M);

% Do PCT
if 0
    M = [M sig1];
    %[M,V] = hyperPct(M,size(M,1));
    [M,V] = hyperPct(M,55);
    sig1 = M(:,end);
    M = M(:,1:end-1);
    p = size(M,1);
    goodBands = 1:p;
end

%q = hyperHfcVd(M);
q = 39;

algorithm = lower(algorithm);
tic
switch algorithm
    case 'ica-eea'
        [U, X] = hyperIcaEea(M, 50, sig1);
        r = X(1,:);
        r = hyperConvert3d(r, w, h, 1);
    case 'rx'
        r = hyperConvert3d(hyperRxDetector(M), w, h, 1);
    case 'matchedfilter'
        r = hyperConvert3d(hyperMatchedFilter(M, sig1), w, h, 1);
    case 'ace'
        r = hyperConvert3d(hyperAce(M, sig1), w, h, 1);
    case 'mace'
        r = hyperConvert3d(hyperMace(M, sig1), w, h, 1);
    case 'sid'
        r = hyperConvert3d(hyperSid(M, sig1), w, h, 1);
    case 'cem'
        r = hyperConvert3d(hyperCem(M, sig1), w, h, 1);
    case 'plmf'
        r = hyperPlmf(hyperConvert3d(M,w,h,p),sig1,9);
    case 'rmf-sum'
        r = hyperRmf(hyperConvert3d(M,w,h,p),sig1,11,'sum');
    case 'rmf-meanlocal'
        r = hyperRmf(hyperConvert3d(M,w,h,p),sig1,11,'meanLocal');
    case 'rmf-meangloballocal'
        r = hyperRmf(hyperConvert3d(M,w,h,p),sig1,11,'meanGlobalLocal');
    case 'glrt'
        r = hyperConvert3d(hyperGlrt(M, sig1), w, h, 1);
    case 'osp'
        U = hyperAtgp(M, q, sig1);
        r = hyperConvert3d(hyperOsp(M, U, sig1), w, h, 1);
    case 'amsd'
        r = hyperConvert3d(hyperAmsd(M, U, sig1), w, h, 1);
    case 'hud'
        U = hyperAtgp(M, q, sig1);
        r = hyperConvert3d(hyperHud(M, U, sig1), w, h, 1);
    case 'nnls'
        U = hyperAtgp(M, q, sig1);
        r = hyperConvert3d(hyperNnls(M,U),w,h,q);
        r = r(:,:,1);
    case 'fcls'
        U = hyperAtgp(M, q, sig1);
        r = hyperConvert3d(hyperFcls(M,U),w,h,q);
        r = r(:,:,1);
    case 'ucls'
        U = hyperAtgp(M, q, sig1);
        r = hyperConvert3d(hyperUcls(M,U),w,h,q);
        r = r(:,:,1);
    case 'sam'
        r = (1./(eps+hyperConvert3d(hyperSam(M, sig1), w, h, 1)));
    otherwise
        error('Incorrect algorithm name specified!\n');
end
toc

% Display results and write to file
figure; imagesc(r); axis image; colorbar;
    title(algorithm);

[a,b]=sort(r(:),'descend');
tmp = a(1:20);
figure; plot(tmp./tmp(1)); grid on;
[x, y, val] = hyperMax2d(r);

% d1 = r(122,494)
% d2 = r(127,490)
% N = prod(size(r));
% [v] = sort(r(:),'ascend');
% idx = find(v==d2);
% N-idx
% figure; hist(r(:),100);

tmp = (hyperNormalize(r)*2^10);
multibandwrite(tmp, sprintf('%s\\results.img', outputDir), 'bil', 'PRECISION', 'int16', 'MACHFMT', 'ieee-le');

[pd,fa] = hyperRoc(r);
figure; plot(fa,pd,'.'); grid on; title(sprintf('%s\n%s',algorithm, targetFilenames{1}));

author: davidkun

repo: HyperSpectralToolbox


text:
function varargout = boundingbox(F,ops,x,fun_eval)
%BOUNDINGBOX Computes bounding box of a constraint
%
% If only one output is requested, only the symbolic model is returned
%  B = boundingbox(F)
%
% If three outputs are requested, numerical values are returned too
%  [B,L,U] = boundingbox(F)
%
% A second argument can be used to specificy solver settings
%  B = boundingbox(F,sdpsettings('solver','cplex'))
%
% A third argument can (should!) be used to obtain a bounding box in a
% particular set of variables, i.e. the bounding box of a projection.
% Unless you specify which variables you are computing the bounding box
% w.r.t, it will be very hard for you to understand which variables the
% values in L and U relate to
%  [B,L,U] = boundingbox(F,[],x)
% B will now be the box [L <= x <= U] (infinite bounds not included)

if nargin < 3 | isempty(x)
    x = recover(depends(F));
else
    x = x(:);
    if ~isa(x,'sdpvar')
        error('The third argument should be an SDPVAR obeject');
    end
end

if nargin < 2
    ops = sdpsettings('verbose',0);
end
if isempty(ops)
    ops = sdpsettings('verbose',0);
end
if ~isa(ops,'struct')
    error('The second argument should be an SDPSETTINGS struct (or empty)');
end

sol = solvesdp(F,[x;-x],ops);
n = length(x);
sols = {};
vals = {};
for i = 1:n
    xi = x(i);
    if isa(xi,'sdpvar')
        if sol.problem(i)==0
            L(i,1) = double(xi,i);
            if nargin > 3
              sols{end+1} = double(x,i);
              vals{end+1} = double(fun_eval,i);
            end
        else
            L(i,1) = -inf;
        end
        if sol.problem(n+i)==0
            U(i,1) = double(xi,n+i);
            if nargin > 3
              sols{end+1} = double(x,n+i);
              vals{end+1} = double(fun_eval,n+i);
            end
        else
            U(i,1) = inf;
        end
    else
        L(i,1) = xi;
        U(i,1) = xi;
    end
end

% Only add finite bounds
Lf = find(~isinf(L));
Uf = find(~isinf(U));
B = [];
if ~isempty(Lf)
    xLf = x(Lf);
    if isa(xLf,'sdpvar')
        B = [B, (xLf >= L(Lf)):'Finite lower bounds'];
    end
end
if ~isempty(Uf)
    xUf = x(Uf);
    if isa(xUf,'sdpvar')
        B = [B, (xUf <= U(Uf)):'Finite upper bounds'];
    end
end
varargout{1} = B;
varargout{2} = L;
varargout{3} = U;
varargout{4} = sols;
varargout{5} = vals;


text:
%% Copyright (C) 2014, 2016, 2018-2019 Colin B. Macdonald
%% Copyright (C) 2016 Lagu
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym imag (@var{z})
%% Imaginary part of a symbolic expression.
%%
%% Examples:
%% @example
%% @group
%% syms z
%% imag(z)
%%   @result{} ans = (sym) im(z)
%% @end group
%%
%% @group
%% syms x real
%% imag(x)
%%   @result{} ans = (sym) 0
%% imag(1i*x)
%%   @result{} ans = (sym) x
%% @end group
%%
%% @group
%% imag([x  sym(pi) + 6i  7  3i])
%%   @result{} ans = (sym) [0  6  0  3]  (1×4 matrix)
%% @end group
%% @end example
%%
%% @seealso{@@sym/real, @@sym/conj, @@sym/ctranspose}
%% @end defmethod

function y = imag(z)
  if (nargin ~= 1)
    print_usage ();
  end

  y = elementwise_op ('im', z);

end

%!assert (isequal (imag (sym (4) + 3i),3))

%!test
%! syms x y real
%! z = x + 1i*y;
%! assert (isequal (imag (z),y))

%!test
%! syms x y real
%! Z = [4  x + 1i*y; 1i*y 4 + 3i];
%! assert (isequal (imag (Z),[0 y; y 3]))

%!test
%! syms x real
%! d = exp (x*i);
%! assert (isequal (imag (d), sin (x)))

%!test
%! % round trip
%! syms x
%! d = 3 - 5i;
%! f = imag (x);
%! A = imag (d);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B)

repo: octsympy


text:
function right_impact_constraints(nlp, src, tar, bounds, varargin)

    plant = nlp.Plant;

    % fist call the class method
    plant.rigidImpactConstraint(nlp, src, tar, bounds, varargin{:});

    % Don't need time continuity constraint
    removeConstraint(nlp,'tContDomain');

    % the relabeling of joint coordiante is no longer valid
    removeConstraint(nlp,'xDiscreteMapRightImpact');

    % the configuration only depends on the relabeling matrix
    R = plant.R;
    x = plant.States.x;
    xn = plant.States.xn;
    x_diff = R*x-xn;
    x_map = SymFunction(['xDiscreteMap' plant.Name],x_diff(3:end),{x,xn});

    addNodeConstraint(nlp, x_map, {'x','xn'}, 'first', 0, 0, 'Linear');

end

author: ayonga

repo: frost-dev

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


text:
function item = update_deps(item, oid, nid)

% function item = update_deps(item, varargin)
% This function will run cfg_dep/update_deps in all leaf (cfg_entry,
% cfg_menu, cfg_files) nodes of a configuration tree and update their
% dependency information (mod_job_ids) if necessary. It will also update
% cfg_exbranch mod_job_ids (item.id and ids in item.tdeps and item.sdeps).
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: update_deps.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

id = subsref(item, substruct('.', 'id'));
for k = 1:numel(oid)
    if isequal(id, oid{k})
        item = subsasgn(item, substruct('.', 'id'), nid{k});
        for cs = 1:numel(item.sout)
            item.sout(cs).src_exbranch = nid{k};
        end;
        break;
    end;
end;
if ~isempty(item.tdeps)
    item.tdeps = update_deps(item.tdeps, oid, nid);
end;
if ~isempty(item.sdeps)
    item.sdeps = update_deps(item.sdeps, oid, nid);
end;

val = subsref(item, substruct('.','val'));
for k = 1:numel(val)
    val{k} = update_deps(val{k}, oid, nid);
end;
item = subsasgn(item, substruct('.','val'), val);

author: spm


text:
function sitem = cfg2struct(item)

% function sitem = cfg2struct(item)
% Return a struct containing all fields of item plus a field type. This is
% the method suitable for cfg_exbranch.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg2struct.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

% Get parent struct, re-classify as field 'type'
sitem = cfg2struct(item.cfg_branch);
sitem.type = class(item);

% Need to cycle through added fields
fn = mysubs_fields;
for k = 1:numel(fn)
    sitem.(fn{k}) = item.(fn{k});
end;


text:
function obj=geom_vline(obj,varargin)
% geom_abline Display vertical reference lines in each facet
%
% Example syntax: gramm_object.geom_abline('xintercept',1,'style','k--')
% See geom_abline for details

p=inputParser;
my_addParameter(p,'xintercept',0);
my_addParameter(p,'style','k--');
my_addParameter(p,'extent',2);
my_addParameter(p,'linewidth',1);
parse(p,varargin{:});

for obj_ind=1:numel(obj)
    obj(obj_ind).abline=fill_abline(obj(obj_ind).abline,NaN,NaN,p.Results.xintercept,NaN,@(x)x,p.Results.style,p.Results.extent,p.Results.linewidth);
end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
%% import_FBG
% Below is a demonstration of the features of the |import_FBG| function

%%
clear; close all; clc;

%% Syntax
% |D=import_FBG(loadName);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% void show_4dvideo(4dmatrix)
%
function [] = show_4dvideo(video)
  for i = 1 : size(video,4)
    frame = video(:,:,:,i);
    imshow(frame,[],'InitialMagnification','fit');
    disp(i);
    pause(0.01);
  end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%--------------------------------------------------------------------------
%   [vpp] = dbm2vpp(dbm,R)
%--------------------------------------------------------------------------
%   功能:
%   dbm 单位转化为电压值
%--------------------------------------------------------------------------
%   输入:
%           dbm                 signal peak to peak volt
%           R                   阻抗匹配(默认为50Ω)
%   输出:
%           vpp
%--------------------------------------------------------------------------
%   例子:
%   dBm2vpp(0)
%   ans = 0.6325
%   dbm2vpp(1,50)
%   ans = 0.7096
%--------------------------------------------------------------------------
function vpp = dbm2vpp(dbm,R)
if nargin==1
    R=50;
end
vpp = 2*10.^((dbm-30+10*log10(2.*R))./20);
end

repo: SP


text:
function i4vec_amax_test ( )

%*****************************************************************************80
%
%% I4VEC_AMAX_TEST tests I4VEC_AMAX;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_AMAX_TEST\n' );
  fprintf ( 1, '  For an integer vector:\n' );
  fprintf ( 1, '  I4VEC_AMAX:   maximum absolute entry;\n' );

  seed = 123456789;
  b = -n;
  c = n;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Input vector:' );

  aval = i4vec_amax ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum absolute value: %f\n', aval );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lp=dir_lpdf(x,a)
%DIR_LPDF   Log probability density function of uniform Dirichlet
%           distribution
%
%       Description:
%       LP = DIR_LPDF(X, A) returns the lpdf of Dirichlet distribution
%       with A at X

% Copyright (c) 2000 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

lp=gammaln(sum(a))-sum(gammaln)+sum(log(x).^(a-1));


text:
%CVX_END_SET  Completes a cvx set specification.

cvx_end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
% DYNAMICPLOT
%   obj = dynamicplot(genFcn, T)
%    creates a figure that is automatically updated (by calling a
%    generator function genFcn) every T seconds. genFcn does not take
%    any arguments, and returns a vector of values to be plotted.
%    obj is a handle to the created object.
%
% Example:
%    d = dynamicplot(@() randn(1,100), 0.5);
%    methods(d)  % shows all the methods
%    help dynamicplot/delete % shows help on an individual method
%
% -------------------------------------------------------------------
%
% The DYNAMICPLOT class demonstrates how to get reference (aka. handle)
% semantics using MATLAB classes. Typically when one modifies a MATLAB
% object, it needs be returned to the caller:
%       obj = set(obj, 'myParameter', 10.0);
% But with handle semantics, you can simply do this:
%      set(obj, 'myParameter', 10.0);
% This is especially useful with objects that need to be updated
% autonomously (e.g., acquiring data from a data aquisition device
% or an instrument).
%
% The handle semantics were accomplished as follows. (1) the constructor
% defines the instance variables and (2) returns a set of function handles
% to access these instance variables. Other methods (e.g., set(obj, ...);
% get(obj, ...)), use the function handles to access and modify the data.
% This technique also enforces data-hiding with the object.
%
% Detailed example of handle semantics:
%
%   d = dynamicplot(@() randn(1,100), 0.5);
%   struct(d)
%   ans =
%          getLineHandle: @dynamicplot/getLineHandle
%        setUpdatePeriod: @dynamicplot/setUpdatePeriod
%        getUpdatePeriod: @dynamicplot/getUpdatePeriod
%      getMostRecentData: @dynamicplot/getMostRecentData
%         getTimerObject: @dynamicplot/getTimerObject
%     getNumAcquisitions: @dynamicplot/getNumAcquisitions
%        getFigureHandle: @dynamicplot/getFigureHandle
%          checkValidity: @dynamicplot/checkValidity
%            setValidity: @dynamicplot/setValidity
%
%
%   getdata(d) % issue this command several times
%   setUpdatePeriod(d, 0.2)
%   getNumUpdates(d) % issue this several times
%   setLineStyle(d,'color','r');
%
%   % Object assignment respects the handle semantics
%   d2 = d;  % both d and d2 refer to the same objects
%   pauseDisplay(d);
%   resumeDisplay(d2);
%
%   % Validity checks can be built into the object itself
%   delete(d);
%   getdata(d) % will generate an error
%   getdata(d2) % this generates error also
%
%   % It is possible to have multiple objects without any clashes
%   d1 = dynamicplot(@() randn(1,100), 0.1);
%   d2 = dynamicplot(@() rand(1,100), 0.1);
%   pauseDisplay(d1);
%   getNumUpdates(d1) % issue these a few times
%   getNumUpdates(d2)
%   resumeDisplay(d1);
%   delete(d1);
%
%
%  Gautam Vallabha, Oct-30-2007, Gautam.Vallabha@mathworks.com

%  Revision 1.0, Oct-30-2007
%    - File created

function s = dynamicplot(genFcn, period)

% class variables
isValid = true;

hFigure = figure('name', 'Dynamic Plot');
hAxes = axes('position', [0.13 0.11 0.775 0.815]);
hLine = plot(0,0, 'parent', hAxes);

generatorFcn = genFcn;
mostRecentData = [];
numberOfAcquisitions = 0;
updatePeriod = period;

timerObj = timer('period', updatePeriod, ...
                 'executionMode', 'fixedRate', ...
                 'timerFcn', @timerCallback);
start(timerObj);

% The class struct contains only function handles

s.getLineHandle = @getLineHandle;
s.setUpdatePeriod = @setUpdatePeriod;
s.getUpdatePeriod = @getUpdatePeriod;
s.getMostRecentData = @getMostRecentData;
s.getTimerObject = @getTimerObject;
s.getNumAcquisitions = @getNumAcquisitions;
s.getFigureHandle = @getFigureHandle;
s.checkValidity = @checkValidity;
s.setValidity = @setValidity;

s = class(s, 'dynamicplot');

%% -------------------
% These are essentially private get/set methods

    function out = checkValidity()
        out = isValid;
        if ~isValid
            error('Invalid DYNAMICPLOT object');
        end
    end

    function setValidity(newState)
        isValid = logical(newState);
    end

    function out = getFigureHandle()
       out = hFigure;
    end

    function out = getNumAcquisitions()
        out = numberOfAcquisitions;
    end

    function out = getTimerObject()
        out = timerObj;
    end

    function out = getLineHandle()
        out = hLine;
    end

    function setUpdatePeriod(newPeriod)
        updatePeriod = newPeriod;
    end

    function out = getUpdatePeriod()
        out = updatePeriod;
    end

    function out = getMostRecentData()
        out = mostRecentData;
    end

%% -------------------
    function timerCallback(hobj, eventinfo) %#ok<INUSD>
        mostRecentData = generatorFcn();
        if ~ishandle(hLine)
           % someone deleted the figure, create a new one
           hFigure = figure('name', 'Dynamic Plot');
           hAxes = axes('position', [0.13 0.11 0.775 0.815]);
           hLine = plot(0,0, 'parent', hAxes);
        end

        set(hLine, 'xdata', 1:length(mostRecentData), 'ydata', mostRecentData);
        drawnow;
        numberOfAcquisitions = numberOfAcquisitions + 1;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function rouderFigure2
% Runs the simulations of Rouder et al. 2009 in Figure 2,
% Shows the relationship between T test results, the Bayes Factor and
% sample size.
T = 1:0.05:10;
N = round(logspace(log10(5),log10(5000),100));
criticalT = nan(numel(N),3);
cntr =0;
for n=N
    cntr= cntr+1;
    for t= T
        % We call the bf.ttest function directly with
        % the results of a t-test (e.g. as returned by ttest.m)
        stats.tstat = t;
        stats.df    = n-1;
        stats.N     = n;
        stats.tail  = 'both';
        stats.p     = NaN; % no tused
        % Call the class member
        bf10 = bf.ttest([],'stats',stats);
        % Check whether we;ve reach any of the thresholds.
        % If so, store.
        if bf10>3 && isnan(criticalT(cntr,1))
            criticalT(cntr,1) = t;
        end
         if bf10>10 && isnan(criticalT(cntr,2))
            criticalT(cntr,2) = t;
         end
         if bf10>30 && isnan(criticalT(cntr,3))
            criticalT(cntr,3) = t;
            break; % Goto next n
         end
    end
end

%%
figure(2);
clf;
plot(N,criticalT)
xlabel 'Sample Size'
ylabel 'Critical t-value'
legend('BF=3','BF=10','BF=30')
set(gca,'XScale','Log','YLim',[1 6],'YTick',2:6,'XTick',[5 20 50 200 1000 5000],'XLim',[4 6000])
end

repo: fieldtrip


text:
function testspecgram(data)

%cd 'C:\Documents and Settings\Admin\Desktop\';
%data=wavread('bird109_26519_on_Aug_19_16_33.wav');

params.tapers = [3 5];
params.fpass = [100 20000];
params.Fs = 44100;
params.pad = 2;
max_time =10; % seconds per run
max_tapers = 150;
increment = 1.5;
profile on

nsamples = 1000;
if 1
slow_results = [];

while 1
    tic
    [S,t,f] = mtspecgramc_slow( data(1:nsamples), [0.01 0.001], params );
    time = toc;
    result = [nsamples time];
    fprintf( 'ran %d samples in %d seconds\n',nsamples, time );
    slow_results = [slow_results ;result];
    if time > max_time
        break
    end
    nsamples = round(nsamples * increment);
end
slow_results
fig=figure();
ax=axes('XScale','log','YScale','log');
axes(ax);
h=line( 'Xdata',slow_results(:,1),'Ydata',slow_results(:,2),'Marker','*');
xlabel('number of samples')
ylabel('time');
title('Original mtspecgramc');
grid on
drawnow;
saveas(fig,'datalength_slow.png');

nsamples = 1000;
fast_results = [];
while 1
    tic
    [S,t,f] = mtspecgramc( data(1:nsamples), [0.01 0.001], params );
    time = toc;
    result = [nsamples time];
    fprintf( 'ran %d samples in %d seconds\n',nsamples, time );
    fast_results = [fast_results ;result];
    if time > max_time
        break
    end
    nsamples = round(nsamples * increment);
end
fast_results
fig=figure();
ax=axes('XScale','log','YScale','log');
axes(ax);
h=line( 'Xdata',fast_results(:,1),'Ydata',fast_results(:,2),'Marker','*');
xlabel('number of samples')
ylabel('time');
title('Modified mtspecgramc - preallocate space');
grid on
drawnow;
saveas(fig,'datalength_fast.png');

compare = [];
n = 1;
while n <= min(length(slow_results(:,1)),length(fast_results(:,1)))
    compare_one = [slow_results(n,1) slow_results(n,2)/fast_results(n,2)];
    compare = [compare ;compare_one];
    n = n + 1;
end
compare
fig=figure();
ax=axes('XScale','log','YScale','lin');
axes(ax);
h=line( 'Xdata',compare(:,1),'Ydata',compare(:,2),'Marker','*');
title('Preallocation slowdown/speedup of mtspecgramc');
xlabel('number of samples')
ylabel('speedup');
grid on
drawnow;
saveas(fig,'speedup.png');

end;

nsamples=10000;
results = [];
n = 1;
while 1
    tic
    params.tapers = [n (2*n-1)];
    [S,t,f] = mtspecgramc( data(1:nsamples), [0.01 0.001], params );
    time = toc;
    result = [params.tapers(2) time];
    fprintf( 'ran %d samples in %d seconds with tapers %d %d\n',nsamples, time,params.tapers(1),params.tapers(2) );
   results = [results ;result];
    if time > max_time || params.tapers(2) > max_tapers
        break
    end
    n = round(n * increment);
end
fig=figure();
ax=axes('XScale','log','YScale','log');
axes(ax);
h=line( 'Xdata',results(:,1),'Ydata',results(:,2),'Marker','*');
xlabel('tapers')
ylabel('time');

drawnow;
saveas(fig,'tapers.png');

stats = profile('info')

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.01 0.05], [0.1 0.01], [0.1 0.01]);
if ~make_it_tight,  clear subplot;  end

%% Upper and Lower Subplots with Titles
income = [3.2,4.1,5.0,5.6];
outgo = [2.5,4.0,3.35,4.9];
subplot(2,1,1); plot(income)
title('Income')
subplot(2,1,2); plot(outgo)
title('Outgo')

%% Subplots in Quadrants
figure
subplot(2,2,1)
text(.5,.5,{'subplot(2,2,1)';'or subplot 221'},...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,2)
text(.5,.5,{'subplot(2,2,2)';'or subplot 222'},...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,3)
text(.5,.5,{'subplot(2,2,3)';'or subplot 223'},...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,4)
text(.5,.5,{'subplot(2,2,4)';'or subplot 224'},...
    'FontSize',14,'HorizontalAlignment','center')

%% Asymmetrical Subplots
figure
subplot(2,2,[1 3])
text(.5,.5,'subplot(2,2,[1 3])',...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,2)
text(.5,.5,'subplot(2,2,2)',...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,4)
text(.5,.5,'subplot(2,2,4)',...
    'FontSize',14,'HorizontalAlignment','center')

%%
figure
subplot(2,2,1:2)
text(.5,.5,'subplot(2,2,1:2)',...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,3)
text(.5,.5,'subplot(2,2,3)',...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,4)
text(.5,.5,'subplot(2,2,4)',...
    'FontSize',14,'HorizontalAlignment','center')

%% Plotting Axes Over Subplots
figure
y = zeros(4,15);
for k = 1:4
    y(k,:) = rand(1,15);
    subplot(2, 2, k)
    plot(y(k,:));
end
hax = axes('Position', [.35, .35, .3, .3]);
bar(hax,y,'EdgeColor','none')
set(hax,'XTick',[])

author: epfl-lasa


% SubgradientOpt.m // Jon Lee
%
% Apply Subgradient Optimization to
%   z = min c'x
%        Ex  = h
%        Ax  = b
%         x >= 0
%
% NOTE: WE ONLY HANDLE PROBLEMS WITH BOUNDED SUBPROBLEM LPs.
%       IF WE ENCOUNTER AN UNBOUNDED SUBPROBLEM, THEN WE QUIT.

% generate a random example
n = 50; % number of variables
m1 = 15; % number of equations to relax
m2 = 10; % number of equations to keep

rng('default');
rng(1);  % set seed
E = rand(m1,n);  % random m1-by-n matrix
A = rand(m2,n);  % random m2-by-n matrix

w = rand(n,1);     % generate rhs's so problem is feasible
h = E*w;
b = A*w;

c = rand(n,1);     % generate random objective
% first we will calculate the optimal value z of the original LP, just to
% compare later
disp('Optimizing the original LP without Lagrangian relaxation');
[x,z,exitflag] = linprog(c,[],[],[E; A],[h; b],zeros(n,1),[]);
if (exitflag < 1)
     disp('fail 1: original LP (without using LR) did not solve correctly');
     return;
end;

disp('Optimal value of the original LP:');
z

disp('Starting Subgradient Optimization');

MAXIT = 100; % number of iterations

results1 = zeros(MAXIT,1);
results2 = zeros(MAXIT,1);

k = 0; % iteration counter
y = zeros(m1,1); % initialize y as you like
g = zeros(m1,1); % initialize g arbitrarily --- it will essentially be ignored
stepsize = 0;    % just to have us really start at the initialized y
bestlb = -Inf;
while k < MAXIT
   k = k + 1; % increment the iteration counter
   y = y + stepsize*g; % take a step in the direction of the subgradient
   % solve the Lagrangian
   [x,subval,exitflag,output,dualsol] = linprog((c'-y'*E)',[],[],A,b,zeros(n,1),[]);
   v = h'*y + (c'-y'*E)*x;  % calculate the value of the Lagrangian
   disp(v);
   bestlb = max(bestlb,v);
   results1(k)=k;
   results2(k)=v;
   g = h - E*x; % calculate the subgradient
   stepsize = 1/k; % calculate the step size
end

bestlb

z

clf;  % clear figure
% plot the sequence of MAXIT lower bounds
plot(results1,results2,'k--.','MarkerSize',15);
hold on;
% plot a horizontal line at height z
plot([1,MAXIT],[ z z ], 'r-','LineWidth',2.5)
axis tight;
xlabel('Iteration');
ylabel('Lagrangian lower bound');
legend('Lagrangian LB','z','Location','SouthEast');

print -djpeg SubgradientOpt.jpeg;

% Next, if y really solves the Lagrangian dual, and pi is the
% optimal solution to the Lagrangian subproblem corresponding to y,
% then y and pi should be an optimal dual solution to the given problem.
%
% Let's see if y and pi are nearly dual feasible.

pi = - dualsol.eqlin;
Total_Dual_Infeasibility = norm(min(c' - y'*E - pi'*A,zeros(1,n)))


text:
function [cl,clpos,clneg,clrpos,clrneg] = cluster_sigregions(cl,pthr,varargin)
% [cl,clpos,clneg,clrpos,clrneg] = cluster_ttest(clusters,pthr,[behavioral regressor])
% tor wager  Oct 9, 2003
%
% does t-tests on cluster timeseries and all voxels
% which should contain group data (e.g., con* data)
%
% if optional behavioral regressor is entered,
% computes cluster extent nonparametric p-values, and reports
% extent npm p-values for intercept and behavioral reg.
%
% writes output to cluster_ttest_output.txt

fid = fopen('cluster_ttest_output.txt','a');

k=1;
if length(varargin) > 0, beh = varargin{1};,k=size(beh,2);,end

outstr1 = []; outstr2 = []; outstr3 = []; outstr4 = [];
try,outstr = sprintf('\n%s\n',[cl(1).title ' ' cl(1).imP(1,:)]);
catch, outstr = sprintf('\n%s\n',[cl(1).title ' ' cl(1).imnames(1,:)]);,
end
outstr = [outstr sprintf('\nCl.\tVoxels\tThreshold\tMask\tData\t\n')];
outstr1 = [outstr1 sprintf('\tx\ty\tz\tVoxels\tCluster corr. p.\tMax. t\tCorrected p\tDirection of effect\tr\n')];
outstr2 = outstr1; outstr3 = outstr1; outstr4 = outstr1;

for i = 1:length(cl)

% df = n - k    primary threshold
tthr1 = abs(tinv(pthr,size(cl(i).all_data,1) - k));

% ----------------------------------------
% * do nonparametric permutation
% ----------------------------------------
if ~isfield(cl(i),'npm_ttest'),cl(i).npm_ttest = [];,end
if isfield(cl(i).npm_ttest,'tmax') & isfield(cl(i).npm_ttest,'textmax')
    tthr = cl(i).npm_ttest.tthr ;
    rthr = cl(i).npm_ttest.rthr ;
    tnegthr = cl(i).npm_ttest.tnegthr ;
    rnegthr = cl(i).npm_ttest.rnegthr ;
    tmax = cl(i).npm_ttest.tmax ;
    tmin = cl(i).npm_ttest.tmin ;
    rmax = cl(i).npm_ttest.rmax ;
    rmin = cl(i).npm_ttest.rmin ;
    textmax = cl(i).npm_ttest.textmax ;
    rextmax = cl(i).npm_ttest.rextmax ;
else

[tthr,rthr,tnegthr,rnegthr,tmax,rmax,tmin,rmin,textmax,rextmax] = ...
    npm_ttest(cl(i).all_data,1000,beh,pthr,cl(i).XYZ);
close

cl(i).npm_ttest.tthr = tthr;
cl(i).npm_ttest.rthr = rthr;
cl(i).npm_ttest.tnegthr = tnegthr;
cl(i).npm_ttest.rnegthr = rnegthr;
cl(i).npm_ttest.tmax = tmax;
cl(i).npm_ttest.rmax = rmax;
cl(i).npm_ttest.textmax = textmax;
cl(i).npm_ttest.rextmax = rextmax;
cl(i).npm_ttest.tmin = tmin;
cl(i).npm_ttest.rmin = rmin;

end

%outstr = [outstr sprintf('%3.0f\t%3.0f\t%3.4f\t%s\t%s\t\n',i,cl(i).numVox,pthr,cl(i).P,cl(i).imP(1,:))];

% ----------------------------------------
% * get sig clusters
% ----------------------------------------

clpos(i).from_cluster = i;
clpos(i).M = cl(i).M;
clpos(i).voxSize = cl(i).voxSize;
clpos(i).threshold = tthr1;
clpos(i).title = 'Positive effects';

clrpos(i).from_cluster = i;
clrpos(i).M = cl(i).M;
clrpos(i).voxSize = cl(i).voxSize;
clrpos(i).threshold = tthr1;
clrpos(i).title = 'Positive correlations';

if isfield(cl(i).ttest,'t'), sig = cl(i).ttest.t > tthr1;,
else, warning(['No t field for cl(' num2str(i) ]), sig = [0;0];
end

if any(sig(1,:))
    clpos(i).XYZ = cl(i).XYZ(:,find(sig(1,:)));
    clpos(i).XYZmm = cl(i).XYZmm(:,find(sig(1,:)));
    clpos(i).Z = cl(i).Z(:,find(sig(1,:)));
    clpos(i).t = cl(i).ttest.t(:,find(sig(1,:)));
    gopos(i) = 1;

    % print table entry
    wcl = spm_clusters(clpos(i).XYZ);
    clpos(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clpos(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clpos(i).XYZmm(:,wh)';,end
        nv = length(wh);

        if length(wh)>1,max_t = max(clpos(i).t(:,wh)'); max_t = max_t(1);,
        else, max_t = clpos(i).t(1,wh);
        end

        clpos(i).center = xyz;

        extcor_p = 1 - (sum(textmax <= nv) ./ length(textmax));
        max_p = 1 - (sum(tmax <= max_t) ./ length(tmax));

        outstr1 = [outstr1 sprintf('%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t+\n', ...
            clpos(i).from_cluster,xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p)];
    end

else
    gopos(i) = 0;
end

if any(sig(2,:))
    clrpos(i).XYZ = cl(i).XYZ(:,find(sig(2,:)));
    clrpos(i).XYZmm = cl(i).XYZmm(:,find(sig(2,:)));
    clrpos(i).Z = cl(i).Z(:,find(sig(2,:)));
    clrpos(i).t = cl(i).ttest.t(:,find(sig(2,:)));
    gorpos(i) = 1;

    % print table entry
    wcl = spm_clusters(clrpos(i).XYZ);
    clrpos(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clrpos(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clrpos(i).XYZmm(:,wh)';,end
        nv = length(wh);

        tmp = cl(i).all_data(:,find(sig(2,:)));
        tmp = tmp(:,wh);
        tmp = corrcoef([tmp beh]);
        tmp = tmp(end,1:end-1);
        max_r = max(tmp);
        if max_r > .99, warning('problem?'),keyboard,end

        if length(wh)>1,max_t = max(clrpos(i).t(:,wh)'); max_t = max_t(2);,
        else, max_t = clrpos(i).t(2,wh);
        end
        clrpos(i).center = xyz;

        extcor_p = 1 - sum(rextmax <= nv) ./ length(rextmax);
        max_p = 1 - sum(rmax <= max_t) ./ length(rmax);

        outstr2 = [outstr2 sprintf('%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t+\t%3.2f\n', ...
            clrpos(i).from_cluster,xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p,max_r)];
    end

else
    gorpos(i) = 0;
end

% negative

clneg(i).from_cluster = i;
clneg(i).M = cl(i).M;
clneg(i).voxSize = cl(i).voxSize;
clneg(i).threshold = tthr1;
clneg(i).title = 'Negative effects';

clrneg(i).from_cluster = i;
clrneg(i).M = cl(i).M;
clrneg(i).voxSize = cl(i).voxSize;
clrneg(i).threshold = tthr1;
clrneg(i).title = 'Negative correlations';

if isfield(cl(i).ttest,'t'), sig = cl(i).ttest.t < -tthr1;,
else, warning(['No t field for cl(' num2str(i) ]), sig = [0;0];
end

if any(sig(1,:))
    clneg(i).XYZ = cl(i).XYZ(:,find(sig(1,:)));
    clneg(i).XYZmm = cl(i).XYZmm(:,find(sig(1,:)));
    clneg(i).Z = cl(i).Z(:,find(sig(1,:)));
    clneg(i).t = cl(i).ttest.t(:,find(sig(1,:)));
    goneg(i) = 1;

    % print table entry
    wcl = spm_clusters(clneg(i).XYZ);
    clneg(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clneg(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clpos(i).XYZmm(:,wh)';,end
        nv = length(wh);

        if length(wh)>1,max_t = min(clneg(i).t(:,wh)'); max_t = max_t(1);,
        else, max_t = clneg(i).t(1,wh);
        end

        clneg(i).center = xyz;

        extcor_p = 1 - sum(textmax <= nv) ./ length(textmax);
        max_p = 1 - sum(tmin >= max_t) ./ length(tmin);

        outstr3 = [outstr3 sprintf('%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t-\n', ...
            clneg(i).from_cluster,xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p)];
    end

else
    goneg(i) = 0;
end

if any(sig(2,:))
    clrneg(i).XYZ = cl(i).XYZ(:,find(sig(2,:)));
    clrneg(i).XYZmm = cl(i).XYZmm(:,find(sig(2,:)));
    clrneg(i).Z = cl(i).Z(:,find(sig(2,:)));
    clrneg(i).t = cl(i).ttest.t(:,find(sig(2,:)));
    gorneg(i) = 1;

    % print table entry
    wcl = spm_clusters(clrneg(i).XYZ);
    clrneg(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clrneg(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clrneg(i).XYZmm(:,wh)';,end
        nv = length(wh);

        if length(wh)>1,max_t = min(clrneg(i).t(:,wh)'); max_t = max_t(2);,
        else, max_t = clrneg(i).t(2,wh);
        end

        tmp = cl(i).all_data(:,find(sig(2,:)));
        tmp = tmp(:,wh);
        tmp = corrcoef([tmp beh]);
        tmp = tmp(end,1:end-1);
        max_r = min(tmp);
        if max_r > .99, warning('problem?'),keyboard,end

        clrneg(i).center = xyz;

        extcor_p = 1 - sum(rextmax <= nv) ./ length(rextmax);
        max_p = 1 - sum(rmin >= max_t) ./ length(rmin);

        outstr4 = [outstr4 sprintf('%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t-\t%3.2f\t\n', ...
            clrneg(i).from_cluster,xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p,max_r)];
    end

else
    gorneg(i) = 0;
end

end % loop through clusters

outstr = [outstr sprintf('\n')];
outstr1 = [outstr1 sprintf('\n')];
outstr2 = [outstr2 sprintf('\n')];
outstr3 = [outstr3 sprintf('\n')];
outstr4 = [outstr4 sprintf('\n')];

disp([outstr outstr1 outstr2 outstr3 outstr4])
outstr(outstr == '\') = '_';
fprintf(fid,[outstr outstr1 outstr2 outstr3 outstr4]);
fclose(fid);

% Make montage
gomontage = 0;

if gomontage
try

if ~isempty(clpos) & ~isempty(clneg)
    montage_clusters([],cl(unique(cat(2,[clpos.from_cluster clneg.from_cluster]))),clpos,clneg,{'k' 'r' 'b'},'nooverlap')
    %hh=get(gcf,'Children');
    %for i = 1:length(hh),axes(hh(i));, h = findobj(gca,'FaceColor',[0 0 0]);, set(h,'FaceAlpha',.5) ,end
elseif ~isempty(clpos)
    montage_clusters([],cl(cat(2,clpos.from_cluster)),clpos,{'k' 'r'},'nooverlap')
elseif ~isempty(clneg)
    montage_clusters([],cl,clneg,{'k' 'b'},'nooverlap')
else
    disp('No main contrast effects')
end

if ~isempty(clrpos) & ~isempty(clrneg)
    montage_clusters([],cl(unique(cat(2,[clrpos.from_cluster clrneg.from_cluster]))),clrpos,clrneg,{'k' 'r' 'b'},'nooverlap')
elseif ~isempty(clrpos)
    montage_clusters([],cl(cat(2,clrpos.from_cluster)),clrpos,{'k' 'r'},'nooverlap')
elseif ~isempty(clrneg)
    montage_clusters([],cl(cat(2,clrneg.from_cluster)),clrneg,{'k' 'b'},'nooverlap')
else
    disp('No correlation effects')
end

catch
    disp('error with montage')
end
end

try, clpos = redefine_clusters(clpos);,catch,end
try, clrpos = redefine_clusters(clrpos);,catch,end
try, clneg = redefine_clusters(clneg);,catch,end
try, clrneg = redefine_clusters(clrneg);,catch,end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [TorF, vstr, rdate] = have_feature_e4st()
%HAVE_FEATURE_E4ST  Detect availability/version info for E4ST
%
%   Feature detection function implementing 'e4st' tag for HAVE_FEATURE
%   to detect availability/version of E4ST, the Engineering, Economic, and
%   Environmental Electricity Simulation Tool (https://e4st.com).
%
%   See also HAVE_FEATURE.

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

TorF = exist('e4st_ver', 'file') == 2;
if TorF
    v = e4st_ver('all');
    vstr = v.Version;
    rdate = v.Date;
else
    vstr = '';
    rdate = '';
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower