text:
function [displacementField, certainty] = demons2d(moving, fixed, varargin)
% DEMONS2D Estimates a displacement field using the original demons scheme
%
% [displacementField, certainty] = demons2d(moving, fixed)
%
% INPUT ARGUMENTS
% moving            - Moving iamge
% fixed             - Fixed image
%
% OPTIONAL INPUT ARGUMENTS
% 'method'                  - Use the gradient of one of the images or both
%                             'fixed', 'moving', 'symmetric'
%
% 'maxDisplacemet'          - Max displacement of the estimated displacement field
%                             2 (default)
%
% 'multiModal'              - Set whether to perform multi-modal or
%                             uni-modal image registration
%                             false (default), true
%                             Setting to true will force method to 'fixed'
%
% 'numberOfChannels'        - Number of channels to use in when computing
%                             the entropy (based on channel coding). This
%                             is only relevant if multiModal is set to
%                             true.
%                             Default value is 8
%
% OUTPUT ARGUMENTS
% displacementField - Estimated displacement field
% certainty         - Certainty map related to the estimated displacement
%                     field

% Copyright (c) 2011
% danne.forsberg@outlook.com
%
% 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/>.

% Default parameters
method = 'symmetric';
maxDisplacement = 2;

% multi-modal
multiModal = false;

% Only valid for multi-modal registration
numberOfChannels = 8;

% Overwrites default parameters
for k=1:2:length(varargin),
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end

if (multiModal && ~strcmp(method,'fixed'))
    method = 'fixed';
end

displacementField = cell(2,1);
displacementField{1} = zeros(size(moving));
displacementField{2} = zeros(size(moving));

switch method
    case 'moving'
        % Estimate gradient of moving image
        [movingGradX movingGradY] = gradient(fixed);
        gradX = (movingGradX);
        gradY = (movingGradY);
    case 'fixed'
        % Estimate gradient of fixed image
        [fixedGradX fixedGradY] = gradient(fixed);
        gradX = (fixedGradX);
        gradY = (fixedGradY);
    case 'symmetric'
        % Estimate gradient of fixed image
        [fixedGradX fixedGradY] = gradient(fixed);
        % Estimate gradient of moving image
        [movingGradX movingGradY] = gradient(fixed);
        gradX = 0.5*(fixedGradX + movingGradX);
        gradY = 0.5*(fixedGradY + movingGradY);
end

if multiModal
    numerator = estimate_delta_c(fixed,moving,numberOfChannels);
else
    numerator = (fixed - moving);
end

denominator = numerator.^2 + (gradX.^2 +  gradY.^2) + eps;

displacementField{1} = (numerator.*gradX)./denominator;
displacementField{2} = (numerator.*gradY)./denominator;

currentMaxDisplacement = max(sqrt(displacementField{1}(:).^2 + displacementField{2}(:).^2));

if currentMaxDisplacement > maxDisplacement
    normalizingFactor = currentMaxDisplacement / maxDisplacement;
    displacementField{1} = displacementField{1}/normalizingFactor;
    displacementField{2} = displacementField{2}/normalizingFactor;
end

certainty = sqrt(gradX.^2 +  gradY.^2);

author: fordanic

repo: image-registration


text:
function [ops, stat, Fcell, FcellNeu]      = extractSignalsNoOverlaps2(ops, m, stat)

ops.saveNeuropil = getOr(ops, 'saveNeuropil', 0);

Nk = numel(stat);

Ny = numel(ops.yrange);
Nx = numel(ops.xrange);

% make new set of basis functions (larger coverage)
ops.neuropilRange = 10;

S = getNeuropilBasis(ops, Ny, Nx, 'raisedcosyne'); % 'raisedcosyne', 'Fourier'
S = normc(S);

% S = m.S;

S = reshape(S, [], size(S, ndims(S)));
nBasis = size(S,2);
S = double(S);

% initialize mask
maskNeu = ones(size(S,1), 1);

stat = getNonOverlapROIs(stat, Ny, Nx);

% create cell masks and cell exclusion areas
[stat, cellPix, cellMasks] = createCellMasks(stat, Ny, Nx);
cellMasks   = sparse(double(cellMasks(:,:)));

LtS = cellMasks * S;

% LtS = zeros(Nk, size(S,2));
% lam = cell(Nk,1);
% for k = 1:Nk
%     ix = stat(k).ipix(~stat(k).isoverlap);
%
%
%     lam{k} = stat(k).lam(~stat(k).isoverlap)';
%     lam{k} = lam{k}/sum(lam{k}(:));
%
%     maskNeu(stat(k).ipix)= 0;
%     if numel(ix)==0 || sum(~stat(k).isoverlap)==0
%         LtS(k,:) = 0;
%     else
%         LtS(k,:) = lam{k} * S(ix, :);
%     end
% end

% add all pixels within X um
if isfield(ops, 'exclFracCell') && ops.exclFracCell>0
    H       = fspecial('disk', round(ops.diameter * ops.exclFracCell));
    maskNeu = reshape(maskNeu, Ny, Nx);
    maskNeu = imfilter(maskNeu, H, 'replicate');
    maskNeu = single(maskNeu(:) > 1-1e-3);
end
%% get signals
S    = bsxfun(@times, S, maskNeu(:));
StS = S' * S;
StS = StS + 1e-2 * eye(size(StS));
iLtS = LtS * (StS\S');

nimgbatch = 2000;

ix = 0;
fclose all;
fid = fopen(ops.RegFile, 'r');

tic
F        = zeros(Nk, sum(ops.Nframes), 'single');
Fneu    = zeros(Nk, sum(ops.Nframes), 'single');
if ops.saveNeuropil
    Ntraces = zeros(nBasis, sum(ops.Nframes), 'single');
end
% S    = bsxfun(@times, S, maskNeu(:));

[Ly Lx] = size(ops.mimg1);

mimg1 = ops.mimg1(ops.yrange, ops.xrange);
% find indices of good clusters
while 1
    data = fread(fid,  Ly*Lx*nimgbatch, '*int16');
    if isempty(data)
       break;
    end
    data = reshape(data, Ly, Lx, []);

    data = data(ops.yrange, ops.xrange, :);
    data = double(data);

%     data = bsxfun(@minus, data, mimg1);

    NT   = size(data,3);
    data = reshape(data, [], NT);

    % process the data

%     Ftemp = zeros(Nk, NT, 'single');
%     for k = 1:Nk
%        ipix = stat(k).ipix(~stat(k).isoverlap)';
%        if ~isempty(ipix)
%            Ftemp(k,:) = lam{k} * data(ipix,:);
%        end
%     end
    Ftemp               = cellMasks * data;
    F(:,ix + (1:NT))    = Ftemp;

%     Tneu                = StS\(S' * data);
    Ftemp2              = iLtS * data; %LtS * Tneu;
    Fneu(:,ix + (1:NT)) = Ftemp2;

%     Fneu(:,ix + (1:NT))     = m.LtS * Fdeconv(1+Nk:end, :); % estimated neuropil
%     F(:,ix + (1:NT))        = Fneu(:,ix + (1:NT)) + Fdeconv(1:Nk, :); % estimated ROI signal

    if ops.saveNeuropil
        Tneu =  StS\(S' * data);
        Ntraces(:,ix + (1:NT)) = Tneu;
    end

    ix = ix + NT;
    if rem(ix, 3*NT)==0
        fprintf('Frame %d done in time %2.2f \n', ix, toc)
    end
end
fclose(fid);

%%
% get activity stats
[stat, F, Fneu] = getActivityStats(ops, stat, F, Fneu);

%
csumNframes = [0 cumsum(ops.Nframes)];
Fcell       = cell(1, length(ops.Nframes));
FcellNeu    = cell(1, length(ops.Nframes));
for i = 1:length(ops.Nframes)
    Fcell{i}     = F(:, csumNframes(i) + (1:ops.Nframes(i)));
    FcellNeu{i}  = Fneu(:, csumNframes(i) + (1:ops.Nframes(i)));
end

if getOr(ops, 'saveNeuropil', 0)
    S = reshape(S, numel(ops.yrange), numel(ops.xrange), nBasis);
    save(sprintf('%s/NEU_%s_%s_plane%d.mat', ops.ResultsSavePath, ...
        ops.mouse_name, ops.date, ops.iplane),  'ops', 'S', 'Ntraces', '-v7.3')
end

author: cortex-lab


text:
function rotated=uniform_rotation_3D(v)
% This function perform a unformly random rotation 3D.
% It is clear that in 2 dimensions that this means the
% rotation angle is uniformly distributed between (0, 2pi).
% However, it is NOT correct and does NOT carry over to
% higher dimensions.
%
% More information see Le?n, Carlos A.; Mass?, Jean-Claude;
% Rivest, Louis-Paul (February 2006), "A statistical model
% for random rotations", Journal of Multivariate Analysis 97
% (2): 412?430, doi:10.1016/j.jmva.2005.03.009, ISSN 0047-259X

% Arguments (input):
%   v - 3*1 vector needs to be rotated
%
%
% arguments (output):
%   rotated - 3*1 vector after unformly random rotation
%
%
% References: J. Arvo (1992), Fast Random Rotation Matrices
%
% Author: Yaming Wang. ETH Zurich.
% email: yaming.wang@sed.ethz.ch
% Release: 1.00
% Release data: 18.09.2012

rot_1=2*pi*rand(1);
matrix_1=[cos(rot_1),sin(rot_1),0;...
    -sin(rot_1),cos(rot_1),0;...
    0,0,1];

rot_2a=2*pi*rand(1);
rot_2b=rand(1);

Householder_v=[cos(rot_2a)*sqrt(rot_2b)
    sin(rot_2a)*sqrt(rot_2b)
    sqrt(1-rot_2b)];

% Householder matrix
matrix_2=eye(3)-2*Householder_v*Householder_v';

rotated=-matrix_2*matrix_1*v;

fprintf(1,'rotation finished.\n');

% ====================================================
%      end of function
% ====================================================

author: Sable

repo: mcbench-benchmarks


text:
function acpc_aligned_nifti = mrAnatAutoAlignAcpcNifti(anatomical_in, anatomical_out)
%
%  acpc_aligned_nifti = mrAnatAutoAlignAcpcNifti(anatomical_in, [anatomical_out])
%
% Generate a ac-pc aligned nifti image using Talairach coordinates
% converted to MNI space.
%
% INPUTS:
%   anatomical_in:	Full path to raw (unaligned) anatomical image
%   anatomical_out:	(optional) Full path for output file. Defaults to
%                   <anatomical_in> + '_acpc.nii.gz'
%
% OUTPUTS:
%   anatomical_out: Full path to the aligned nifti file.
%
% (C) Stanford Vista Lab, 2016
%

%% Handle I/O

if ~exist('anatomical_in','var') || isempty(anatomical_in) || ~exist(anatomical_in, 'file')
    error('No input file was defined.');
end

% If no output file was defined, set to input + '_acpc.nii.gz'
if ~exist('anatomical_out','var') || isempty(anatomical_out)
    base = strsplit(anatomical_in, '.nii');
    [p, f] = fileparts(base{1});
    anatomical_out = fullfile(p, [f, '_acpc.nii.gz']);
end

%% Perform acpc alignment

% The ACPC coordinates in Talairach
ACPC_COORDS = [0,0,0; 0,-16,0; 0,-8,40];

% Read in the file
ni = niftiRead(anatomical_in);

% Apply the cannonical transform
ni = niftiApplyCannonicalXform(ni);

% Use the MNI_T1 template
template =  fullfile(mrDiffusionDir, 'templates', 'MNI_T1.nii.gz');

% Compute the spatial normalization
sn = mrAnatComputeSpmSpatialNorm(ni.data, ni.qto_xyz, template);

% Get the iamge coordinates for the mid-line, ac, and pc
coords = mrAnatGetImageCoordsFromSn(sn, tal2mni(ACPC_COORDS)', true)';

% Use the coords to generate acpc-aligned image
mrAnatAverageAcpcNifti(ni, anatomical_out, coords, [], [], [], false);

%% Check for anatomical_out file

if ~exist(anatomical_out,'file')
    warning('Auto AC-PC Alignment failed.');
else
    acpc_aligned_nifti = anatomical_out;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]};

%% AR(1)
g = .95;
[Y, trueC, trueSpikes] = gen_data();
N = size(Y, 1);

% run deconvolution
lam = 2.4;
[c_oasis, s_oasis] = oasisAR1(Y(1,:), g, lam);
[c_foopsi, s_foopsi] = foopsi(Y(1,:), g, lam);

% plot results
figure('papersize', [15, 4]);
init_fig;

% c
axes('position', [.05, .57, .95, .37]);
hold on;
plot(Y(1,:), 'color', col{8}/255);
alpha(.7);
plot(trueC(1,:), 'color', col{3}/255, 'linewidth', 1.5);
plot(c_oasis, 'color', col{1}/255);
plot(c_foopsi, '--', 'color', col{7}/255);
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', []);
set(gca, 'ytick', 0:2);
ylabel('Fluor.');
box off;
legend('Data', 'Truth', 'OASIS', 'CVX', 'location', 'northeast', 'orientation', 'horizontal');
% s
axes('position', [.05, .18, .95, .37]);
hold on;
plot(trueSpikes(1,:), 'color', col{3}/255, 'linewidth', 1.5);
plot(s_oasis, 'color', col{1}/255);
plot(s_foopsi, '--', 'color', col{7}/255);
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', get(gca, 'xtick')/30);
set(gca, 'ytick', [0,1]);
xlabel('Time [s]');
ylabel('Activity.');

saveas(gcf, 'fig/traceAR1.pdf');

%% time it
fprintf('\n**************AR 1**************\n');
tic;
for m=1:N
    [c_oasis, s_oasis] = oasisAR1(Y(m,:), g, lam);
end
fprintf('OASIS: %.3f seconds\n', toc);

tic;
for m=1:N
    [c_foopsi, s_foopsi] = foopsi(Y(1,:), g, lam); %#ok<ASGLU>
end
fprintf('FOOPSI: %.3f seconds\n', toc);
fprintf('\n**************AR 1**************\n');

%% AR(2)
% simulation
g = [1.7, -0.712];
sn = 1;
seed = 3;
[Y, trueC, trueSpikes] = gen_data(g, sn, [], [], [], [], [], seed);
[N, T] = size(Y);

% deconvolution
lam = 25;
[c_onnls, s_onnls] = onnls(Y(1,:), g, lam);
[c_foopsi, s_foopsi] = foopsi(Y(1,:), g, lam);
figure('papersize', [15, 4]);
init_fig;

% plot results
figure('papersize', [15, 4]);
init_fig;

% c
axes('position', [.05, .57, .46, .37]);
hold on;
plot(Y(1,:), 'color', col{8}/255, 'linewidth', 0.5);
alpha(.7);
plot(trueC(1,:), 'color', col{3}/255, 'linewidth', 1.5);
plot(c_onnls, 'color', col{1}/255);
plot(c_foopsi, '--', 'color', col{7}/255);
axis tight;
xlim([0, 1200]);
set(gca, 'xtick', 0:300:1500);
set(gca, 'xticklabel', []);
ylim(round([1+min(Y(1,:)), max(Y(1,:))-0.5]));
ylabel('Fluor.');
box off;
% s
axes('position', [.05, .18, .46, .37]);
hold on;
plot(trueSpikes(1,:), 'color', col{3}/255, 'linewidth', 1.5);
plot(s_onnls, 'color', col{1}/255);
plot(s_foopsi, '--', 'color', col{7}/255);
axis tight;
xlim([0, 1200]);
set(gca, 'xtick', 0:300:1500);
set(gca, 'xticklabel', get(gca, 'xtick')/30);
set(gca, 'ytick', [0,1]);
xlabel('Time [s]');
ylabel('Activity.');

%% timeit
fprintf('\n**************AR 2**************\n');
tic;
for m=1:N
    [c_onnls, s_onnls] = onnls(Y(m,:), g, lam);
end
fprintf('online NNLS: %.3f seconds\n', toc);

tic;
for m=1:N
    [c_foopsi, s_foopsi] = foopsi(Y(m,:), g, lam);
end
fprintf('FOOPSI: %.3f seconds\n', toc);
fprintf('\n**************AR 2**************\n');

author: zhoupc


text:
function [pts] = helen_annotation_load( fpath )

%fpath = ['/home/tranngoctrung/Workspace/These/3dpose/data/input/HELEN/points/' '232194_1.txt'];

pts = [];

if exist(fpath,'file') ~= 0

fp = fopen(fpath,'r');

npts = 194; %predefined
pts  = zeros(npts,2);

count = 0;

while ~feof(fp)

    count = count + 1;

    line = fgetl(fp);

    if count > 1
        terms = strsplit(line,',');
        pts(count-1,1) = str2num(terms{1});
        pts(count-1,2) = str2num(terms{2});
    end

end

end

fclose(fp);

end

author: tntrung


text:
classdef linop
%linop A class for linear operators

% written by M. Storath
% $Date: 2015-10-15 15:07:43 +0200 (Do, 15. Okt 2015) $	$Revision: 132 $

    properties
        % function handle for evaluation A * x
        eval;
        % function handle for evaluation A' * x
        ctrans;
        % A'A
        normalOp;
        % true if A'A positive definite
        posdef;
        % true if A'A positive definite
        lseSolver;
    end

    methods
        % constructor
        function A = linop(eval, ctrans, varargin)
            % if it is a matrix
            if ~isa(eval, 'function_handle')
                A.eval = @(x) eval * x;
                A.ctrans = @(x) eval' * x;
                M = eval' * eval;
                A.normalOp = @(x) M * x;
            else
                A.eval = eval;
                A.ctrans = ctrans;
                A.normalOp = @(x) A.ctrans(A.eval(x));
            end
            ip = inputParser;
            addParamValue(ip, 'posdef', false);
            addParamValue(ip, 'lseSolver', []);
            parse(ip, varargin{:});
            par = ip.Results;
            A.posdef = par.posdef;
            A.lseSolver = par.lseSolver;
        end

        % ctranspose (')
        function C = ctranspose(A)
            C =  linop(A.ctrans, A.eval);
        end

        % mtimes (*)
        function C = mtimes( A, B )
            if isa(B, 'linop')
                C = linop( @(x) A.eval(B.eval(x)), @(x) B.ctrans(A.ctrans(x)));
            else
                C = A.eval(B);
            end
        end

        % size
        function s = size(A)
            %warning('Size is deprecated for class linop');
            s = [1 1];
        end

    end
end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
classdef (Abstract) Tire
    % Tire Tire abstract class
    %
    % This class must be inherited by tire models that will be used with the simulator. This class has no properties, since parameters of tire models differ significantly from one another.

    methods(Abstract)
        Characteristic(self, alpha, varargin)
    end
end

%% See Also
%
% <../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function [fg, OK] = dtiModifyFiberGroup(fg)
% Window popup to allow user to set properties of a fiber group in the UI
%
%  [fg, OK] = dtiModifyFiberGroup(fg)
%
% Fiber group name, thickness, visibility, and color can be changed.
%
% HISTORY:
%   2003.10.29 RFD (bob@white.stanford.edu) wrote it.
%
% (c) Stanford VISTA Team

if(~exist('fg','var') || isempty(fg)), fg = dtiNewFiberGroup; end

ans = inputdlg({'name:','thickness:','visible:','color:'},...
    'Fiber Group Info',1,...
    {fg.name, num2str(fg.thickness), num2str(fg.visible), num2str(fg.colorRgb)});

if(~isempty(ans))
    fg.name = ans{1};
    fg.thickness = str2num(ans{2});
    fg.visible = str2num(ans{3});
    fg.colorRgb = str2num(ans{4});
    OK = true;
else
    OK = false;
end

return;


text:
function f = comp_ifilterbank(c,g,a,L)
%COMP_IFILTERBANK Compute inverse filterbank

%   called by comp_ifilterbank, performs frequency domain
%   filtering for filters with finite impulse responses

M = numel(g);
classname = assert_classname(c{1});

% Divide filters into time domain and frequency domain groups
mFreq = 1:M;
mTime = mFreq(cellfun(@(gEl) isfield(gEl,'h') ,g)>0);
mFreq(mTime) = [];

f = [];

if ~isempty(mTime)
   % Pick imp. resp.
   gtime = cellfun(@(gEl) gEl.h, g(mTime),'UniformOutput',0);

   % Call the routine
   gskip = cellfun(@(gEl) gEl.offset ,g(mTime));
   f = comp_ifilterbank_td(c(mTime),gtime,a(mTime),L,gskip,'per');
end

if ~isempty(mFreq)
   % Pick frequency domain filters
   gfreq = g(mFreq);
   % Divide filters into the full-length and band-limited groups
   mFreqFullL = 1:numel(gfreq);
   amFreqCell = mat2cell(a(mFreq,:).',size(a,2),ones(1,numel(mFreq)));
   mFreqBL = mFreqFullL(cellfun(@(gEl,aEl) numel(gEl.H)~=L || (numel(aEl)>1 && aEl(2) ~=1), gfreq(:),amFreqCell(:))>0);
   mFreqFullL(mFreqBL) = [];

   mFreqFullL = mFreq(mFreqFullL);
   mFreqBL = mFreq(mFreqBL);

   F = [];
   if ~isempty(mFreqBL)
      conjG = cellfun(@(gEl) cast(gEl.H,classname), g(mFreqBL),'UniformOutput',0);
      foff = cellfun(@(gEl) gEl.foff, g(mFreqBL));
      % Cast from logical to double.
      realonly = cellfun(@(gEl) cast(isfield(gEl,'realonly') && gEl.realonly,'double'), g(mFreqBL));
      F = comp_ifilterbank_fftbl(c(mFreqBL),conjG,foff,a(mFreqBL,:),realonly);
   end

   if ~isempty(mFreqFullL)
      conjG = cellfun(@(gEl) cast(gEl.H,classname), g(mFreqFullL),'UniformOutput',0);

      % In case some of the filters were BL
      if isempty(F)
         F = comp_ifilterbank_fft(c(mFreqFullL),conjG,a(mFreqFullL));
      else
         F = F + comp_ifilterbank_fft(c(mFreqFullL),conjG,a(mFreqFullL));
      end
   end

   % In case some of the filters were TD
   if isempty(f)
      f = ifft(F);
   else
      f = f + ifft(F);
   end
end

% W = size(c{1},2);
% M = numel(g);
% classname = assert_classname(c{1});
%
% f=zeros(L,W,classname);
%
% % This routine must handle the following cases
% %
% %   * Time-side or frequency-side filters (test for  isfield(g,'H'))
% %
% %   * Cell array or matrix input (test for iscell(c))
% %
% %   * Regular or fractional subsampling (test for info.isfractional)
%
%
% for m=1:M
%     conjG=conj(comp_transferfunction(g{m},L));
%
%     % For Octave 3.6 compatibility
%     conjG=cast(conjG,classname);
%
%     % Handle fractional subsampling (this implies frequency side filters)
%     if isfield(g{m},'H') && numel(g{m}.H)~=L
%         N=size(c{m},1);
%         Llarge=ceil(L/N)*N;
%         amod=Llarge/N;
%
%         for w=1:W
%             % This repmat cannot be replaced by bsxfun
%             innerstuff=middlepad(circshift(repmat(fft(c{m}(:,w)),amod,1),-g{m}.foff),L);
%             innerstuff(numel(g{m}.H)+1:end) = 0;
%             f(:,w)=f(:,w)+(circshift(innerstuff.*circshift(conjG,-g{m}.foff),g{m}.foff));
%         end;
%     else
%         if iscell(c)
%             for w=1:W
%                 % This repmat cannot be replaced by bsxfun
%                 f(:,w)=f(:,w)+(repmat(fft(c{m}(:,w)),a(m),1).*conjG);
%             end;
%         else
%             for w=1:W
%                 % This repmat cannot be replaced by bsxfun
%                 f(:,w)=f(:,w)+(repmat(fft(c(:,m,w)),a(m),1).*conjG);
%             end;
%         end;
%     end;
% end;
%
% f = ifft(f);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
clc
close all
clear
% This file tests safety performance for various approaches

%% Setup and simulations
timestep = 0.1;
system_param.A = [[1 ,timestep, 0]; [0, 1, timestep]; [0, 0, 1]];
system_param.B = [0; 0; timestep];
system_param.ul = -1;
system_param.uu = 1;
system_param.timestep = timestep;
x0 = [-2.0;0.0;1.0];
t0 = 0.0;

gammalist = [0.05, 0.10, 0.15, 0.2];

%% MPC-DCBF simulator
mpcdcbf_simulators = {};
for index = 1:length(gammalist)
    mpcdcbf_simulators{index} = CBFDT(system_param, x0, t0);
    param_mpcdcbf = ParamMPCDCBF(8, gammalist(index), 10.0*eye(3), 10.0*eye(3), 1.0);
    mpcdcbf_simulators{index}.setOpt('mpcdcbf', param_mpcdcbf);
    mpcdcbf_simulators{index}.sim(10.0);
end

%% MPC-GCBF simulator
mpcgcbf_simulators = {};
for index = 1:length(gammalist)
    mpcgcbf_simulators{index} = CBFDT(system_param, x0, t0);
    param_mpcgcbf = ParamMPCGCBF(8, gammalist(index), 10.0*eye(3), 10.0*eye(3), 1.0);
    mpcgcbf_simulators{index}.setOpt('mpcgcbf', param_mpcgcbf);
    mpcgcbf_simulators{index}.sim(10.0);
end

%% NMPC-DCBF simulator
nmpcdcbf_simulators = {};
for index = 1:length(gammalist)
    nmpcdcbf_simulators{index} = CBFDT(system_param, x0, t0);
    param_nmpcdcbf = ParamNMPCDCBF(8, 8, gammalist(index), 10.0*eye(3), 10.0*eye(3), 1.0, 10.0);
    nmpcdcbf_simulators{index}.setOpt('nmpcdcbf', param_nmpcdcbf);
    nmpcdcbf_simulators{index}.sim(10.0);
end

%% Plotting
figure('Renderer', 'painters', 'Position', [0 0 500 400]);
color1 = '[0, 0.4470, 0.7410]';
color2 = '[0.8500, 0.3250, 0.0980]';
color3 = '[0.9290, 0.6940, 0.1250]';
color4 = '[0.4940, 0.1840, 0.5560]';
set(gca,'LineWidth', 0.2, 'FontSize', 20);
hold on;
grid on;
set(gca, 'YScale', 'log');
plot(mpcdcbf_simulators{2}.tlog, -mpcdcbf_simulators{2}.xlog(1, :), 'Color', color2, 'LineWidth', 1.0);
plot(mpcdcbf_simulators{3}.tlog, -mpcdcbf_simulators{3}.xlog(1, :), 'Color', color3, 'LineWidth', 1.0);
plot(mpcdcbf_simulators{4}.tlog, -mpcdcbf_simulators{4}.xlog(1, :), 'Color', color4, 'LineWidth', 1.0);
h_legend = legend('MPC-DCBF ($\gamma=0.10$)', 'MPC-DCBF ($\gamma=0.15$)', 'MPC-DCBF ($\gamma=0.20$)');
set(h_legend, 'Interpreter','latex', 'Location', 'SouthWest');
xlim([0, 10]);
% save data and generate figures
print(gcf, 'figures/safety-mpcdcbf.eps', '-depsc');
print(gcf, 'figures/safety-mpcdcbf.png', '-dpng', '-r800');
save('data/safety-mpcdcbf.mat');

figure('Renderer', 'painters', 'Position', [0 0 500 400]);
set(gca,'LineWidth', 0.2, 'FontSize', 20);
hold on;
grid on;
set(gca, 'YScale', 'log');
plot(mpcgcbf_simulators{2}.tlog, -mpcgcbf_simulators{2}.xlog(1, :), 'Color', color2, 'LineWidth', 1.0);
plot(mpcgcbf_simulators{3}.tlog, -mpcgcbf_simulators{3}.xlog(1, :), 'Color', color3, 'LineWidth', 1.0);
plot(mpcgcbf_simulators{4}.tlog, -mpcgcbf_simulators{4}.xlog(1, :), 'Color', color4, 'LineWidth', 1.0);
h_legend = legend('MPC-GCBF ($\gamma=0.10$)', 'MPC-GCBF ($\gamma=0.15$)', 'MPC-GCBF ($\gamma=0.20$)');
set(h_legend, 'Interpreter','latex', 'Location', 'SouthWest');
xlim([0, 10]);
ylim([0.01, 2]);
% save data and generate figures
print(gcf, 'figures/safety-mpcgcbf.eps', '-depsc');
print(gcf, 'figures/safety-mpcgcbf.png', '-dpng', '-r800');
save('data/safety-mpcgcbf.mat');

figure('Renderer', 'painters', 'Position', [0 0 500 400]);
set(gca,'LineWidth', 0.2, 'FontSize', 20);
hold on;
grid on;
set(gca, 'YScale', 'log');
plot(nmpcdcbf_simulators{1}.tlog, -nmpcdcbf_simulators{1}.xlog(1, :), 'Color', color1, 'LineWidth', 1.0);
plot(nmpcdcbf_simulators{2}.tlog, -nmpcdcbf_simulators{2}.xlog(1, :), 'Color', color2, 'LineWidth', 1.0);
plot(nmpcdcbf_simulators{3}.tlog, -nmpcdcbf_simulators{3}.xlog(1, :), 'Color', color3, 'LineWidth', 1.0);
plot(nmpcdcbf_simulators{4}.tlog, -nmpcdcbf_simulators{4}.xlog(1, :), 'Color', color4, 'LineWidth', 1.0);
h_legend = legend('NMPC-DCBF ($\gamma=0.05$)', 'NMPC-DCBF ($\gamma=0.10$)', 'NMPC-DCBF ($\gamma=0.15$)', 'NMPC-DCBF ($\gamma=0.20$)');
set(h_legend, 'Interpreter','latex', 'Location', 'SouthWest');
xlim([0, 10]);
% save data and generate figures
print(gcf, 'figures/safety-nmpcdcbf.eps', '-depsc');
print(gcf, 'figures/safety-nmpcdcbf.png', '-dpng', '-r800');
save('data/safety-nmpccbf.mat');

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function [ps,params] = fused_sigmoid(w,input_data)
%
% Algorithm: ps = sigmoid( alpha'*input_data +beta)
%
%
% Inputs:
%  w: is [alpha; beta], where alpha is D-by-1 and beta is scalar.
%     Use w=[] to let output ps be an MV2DF function handle.
%     If w is a function handle to an MV2DF then ps is the function handle
%     to the composition of w and this function.
%
%  input_data: D-by-T matrix
%
%
% Outputs:
%   ps: function handle (if w=[], or w is handle), or numeric T-by-1
%   params.get_w0(ssat): returns w0 for optimization initialization,
%                        0<ssat<1 is required average sigmoid output.
%   params.tail: is tail of parameter w, which is not consumed by this
%                function.

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

[dim,n] = size(input_data);
wsz = dim+1;
[whead,wtail] = splitvec_fh(wsz,w);
params.get_w0 = @(ssat) init_w0(ssat,dim);
params.tail = wtail;

y = fusion_mv2df(whead,input_data);
ps = sigmoid_mv2df(y);

    function w0 = init_w0(ssat,dim)
    alpha = zeros(dim,1);
    beta = logit(ssat);
    w0 = [alpha;beta];
    end

end

function test_this()

K = 5;
T = 10;
data = randn(K,T);

ssat = 0.99;
[sys,params] = fused_sigmoid([],data);

w0 = params.get_w0(ssat);
test_MV2DF(sys,w0);

ps = fused_sigmoid(w0,data),

end

author: nesl


text:
% SERVER.SENDREQUEST sends an RPC request to a slave.
function sendRequest(server, slaveIdx, numArgsOut, reqHook, varargin)
  if (slaveIdx > server.slaveCount)
    error('Slave index is out of bounds.');
  end

  % set up output stream
  req.hook = reqHook;
  req.vars = varargin;
  req.nargout = numArgsOut;

  % set up output stream, send object
  if (server.options.emulate)
    server.savedObjects{slaveIdx} = matlab2java(req);
  else
    if (slaveIdx > length(server.slave_sock))
      error('No socket for slave.');
    end

    os=server.slave_sock{slaveIdx}.getOutputStream();
    oos=java.io.ObjectOutputStream(os);
    oos.writeObject(matlab2java(req));
  end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function SNR = svargplvmSNR(model)

% SVARGPLVMSNR A small utility to display the Signal to Noise ratio of an optimised svargplvm model
% DESC The signal to noise ratio is defined for every modality, and it is shows the amount of  true signal
% vs the amount of "noise" which is assumed to have generated the data.
% SEEALSO : svargplvmOptimiseModel
%
% COPYRIGHT : Andreas C. Damianou, 2012

% VARGPLVM

for i=1:model.numModels
    if model.comp{i}.DgtN && isfield(model.comp{i}, 'mOrig') && ~isempty(model.comp{i}.mOrig)
        varData = var(model.comp{i}.mOrig(:));
    else
        varData = var(model.comp{i}.m(:));
    end
    SNR{i} = varData/(1/model.comp{i}.beta);
end

author: SheffieldML

repo: vargplvm


text:
function plotsig(C,sig,t,f,c)
% Function to plot C where it is higher than a threshold sig
% useful for plotting coherence
% Usage: plotsig(C,sig,t,f)
% Inputs:
% C: input array t x f - also works for a single vector
% sig: significance level
% t: t axis grid for plot
% f: f axis grid for plot.
% c: color to use (default blue)-only meaningful for a line plot
if nargin < 4; error('Need at least 4 arguments'); end;
if nargin < 5 | isempty(c); c='b'; end;
[T,F]=size(C);
if F==1; C=C'; [T,F]=size(C);end;
if T~=length(t) | F~=length(f);
    error('frequency and/or time axes are incompatible with data');
end;
if T==1;
    dim=max(T,F);
    C=C(:);
    indx=find(C>sig);
    plot(f,C,c);
%     plot(f,C,f,mask.*C)
    line(get(gca,'xlim'),[sig sig]);
    xlabel('f'); ylabel('|C|');
else
    mask=zeros(T,F);
    for n=1:length(t);
        for m=1:length(f);
           if C(n,m)>sig
              mask(n,m)=1;
           end;
        end;
    end;
    imagesc(t,f,(mask.*C)'); axis xy; colorbar
    xlabel('t'); ylabel('f');
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function test_bug2163

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_spike read_neuralynx_nse

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2163/test.nse');
spike = ft_read_spike(filename);
assert(size(spike.waveform{1},2)==64); % there should be 64 samples per spike snippet in this test file

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2163/Sc1.nse');
spike = ft_read_spike(filename);
assert(size(spike.waveform{1},2)==32); % this is an original file with 32 samples

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% NMF-ALS: NMF solved by Alternating Least Squares
% process_video('NMF', 'NMF-ALS', 'dataset/demo.avi', 'output/demo_NMF-ALS.avi');
alg_path_aux = fullfile(lrs_conf.nmf_path,'NMF-DTU-Toolbox');
addpath(genpath(alg_path_aux));

M = sparse(M);

% als: Alternating least squares.
[W, H] = nmf(M,1,'als');

L = W * H;
S = M - L;

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function f = circconv(f, g)
%CONV   Circular convolution of TRIGTECH objects.
%   H = CIRCCONV(F, G) produces the convolution of TRIGTECH objects F and G:
%                     -
%                    /
%           H(x) =   |    F(t) G(x-t) dt,  x in [-pi, pi]
%                    /
%                   -
%   Note that CIRCCONV only supports smooth periodic functions on [-pi,pi].
%
%   Example:
%     f = trigtech(@(x) exp(cos(40*pi*x)));
%     g = trigtech(@(x) exp(-(20*x).^2);
%     h = circconv(f,g);
%     plot(h);

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

% Return empty for an empty input:
if ( isempty(f) || isempty(g) )
    f = trigtech();
    return
end

% No support for array-valued trigtech objects:
if ( (size(f, 2) > 1) || (size(g, 2) > 1) )
    error('CHEBFUN:TRIGTECH:conv:array', ...
        'No support for array-valued TRIGTECH objects.');
end

% Get the sizes of the TRIGTECH objects
nf = size(f.coeffs, 1);
ng = size(g.coeffs, 1);

% Make the TRIGTECH objects the same length.
if ( nf > ng )
    % Increase the length of g (via PROLONG):
    g = prolong(g, nf);
elseif ( nf < ng )
    % Increase the length of f (via PROLONG):
    f = prolong(f, ng);
end
n = size(f.coeffs,1);

% Convolution is just multiplication of the Fourier coefficients.
% Shift g horizontally to -1.
g = circshift(g,-1);
f.values = 2/n*ifft(fft(f.values).*fft(g.values));
f.coeffs = f.vals2coeffs(f.values);

% TODO:  Why do we simplify twice?  (Once here and once below.)
f = simplify(f);

f.ishappy = f.ishappy && g.ishappy;
f.isReal = f.isReal && g.isReal;  % Are you real happy though?

f.values(:,f.isReal) = real(f.values(:,f.isReal));

if ( f.ishappy )
    f = simplify(f);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [p,ll]=parametersfbs(xf,xs,y,wstilda,wstilda2,wstilda3,p)
T=length(xf(1,:));

for t=2:T
ep1(t)=p(7)*normpdf(xs(t),p(1)*xs(t-1)+p(5),sqrt((1+p(6)))*p(2))/(p(7)*normpdf(xs(t),p(1)*xs(t-1)+p(5),sqrt((1+p(6)))*p(2))+(1-p(7))*normpdf(xs(t),p(1)*xs(t-1),p(2)));
ep2(t)=1-ep1(t);
end

e11(T)=(xf(:,T).*xf(:,T))'*wstilda3(:,T);
for t=1:T-1
e21(t)=(xf(:,t+1).*xf(:,t))'*wstilda2(:,t);
e11(t)=(xf(:,t).*xf(:,t))'*wstilda3(:,t);
eyf(t+1)=exp(xf(:,t+1))'*wstilda(:,t+1);
eff(t+1)=exp(2*xf(:,t+1))'*wstilda(:,t+1);
end

esumx2x2p1=sum(e11(2:end).*ep1(2:end));
esumx2x1p1=sum(e21.*ep1(2:end));
esumx1x1p1=sum(e11(1:end-1).*ep1(2:end));
esumx1p1=sum(xs(1:end-1).*ep1(2:end));
esumx2p1=sum(xs(2:end).*ep1(2:end));
esump1=sum(ep1(2:end));

esumx2x2p2=sum(e11(2:end).*ep2(2:end));
esumx2x1p2=sum(e21.*ep2(2:end));
esumx1x1p2=sum(e11(1:end-1).*ep2(2:end));

top=(1/((p(6)+1)*p(2)^2))*(esumx2x1p1-p(5)*esumx1p1)+(1/p(2)^2)*esumx2x1p2;
bot=(1/((p(6)+1)*p(2)^2))*esumx1x1p1+(1/p(2)^2)*esumx1x1p2;
p(1)=top/bot;

first=esumx2x2p1-2*p(5)*esumx2p1+sum(p(5)^2*ep1(2:end))-2*p(1)*esumx2x1p1+2*p(1)*p(5)*esumx1p1+p(1)^2*esumx1x1p1;
second=esumx2x2p2-2*p(1)*esumx2x1p2+p(1)^2*esumx1x1p2;
p(2)=abs(((1/(p(6)+1))*first+second)/(T-1))^.5;

sumyy=sum(y.*y);
sumyfx=sum(y.*eyf);
sumfxfx=sum(eff(2:end));
p(3)=sumyfx/sumfxfx;

p(4)=sqrt(abs(sumyy-2*p(3)*sumyfx+p(3)^2*sumfxfx)/T);

p(5)=abs((esumx2p1-p(1)*esumx1p1)/sum(ep1(2:end)));

p(6)=abs(first/(p(2)^2*esump1)-1);
p(7)=esump1/T;

ll=likelihoodj(y,xs,p(1),p(2),p(3),p(4),p(5),p(6),p(7));


text:
function cS =  biotite
% Mica

cs_Bt = crystalSymmetry('2/m',[0.577 1 1.105], [90, 100.2, 90]*degree);
N = Miller({0,0,1},{0,1,0},{3,3,1},cs_Bt);
dist = [0.2, 1, 5];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function x = sample_paths_eigen ( n, n2, rhomax, rho0, correlation )

%*****************************************************************************80
%
%% SAMPLE_PATHS_EIGEN: sample paths for stationary correlation functions.
%
%  Discussion:
%
%    This method uses the eigen-decomposition of the correlation matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points on each path.
%
%    Input, integer N2, the number of paths.
%
%    Input, real RHOMAX, the maximum value of RHO.
%
%    Input, real RHO0, the correlation length.
%
%    Input, @CORRELATION, a handle for a correlation function.
%
%    Output, real X(N,N2), the sample paths.
%

%
%  Choose N equally spaced sample points from 0 to RHOMAX.
%
  rho_vec = linspace ( 0.0, rhomax, n );
%
%  Evaluate the correlation function.
%
  cor_vec = correlation ( n, rho_vec, rho0 );

  cor = zeros ( n, n );
%
%  Construct the correlation matrix;
%
%  From the vector
%    [ C(0), C(1), C(2), ... C(N-1) ]
%  construct the vector
%    [ C(N-1), ..., C(2), C(1), C(0), C(1), C(2), ...  C(N-1) ]
%  Every row of the correlation matrix can be constructed by a subvector
%  of this vector.
%
  cor_vec = [ cor_vec(n:-1:2)', cor_vec(1:n)' ];
  for i = 1 : n
    cor(i,1:n) = cor_vec(n+1-i:2*n-i);
  end
%
%  Get the eigendecomposition of COR:
%
%    COR = V * D * V'.
%
%  Because COR is symmetric, V is orthogonal.
%
  [ v, d ] = eig ( cor );
%
%  We assume COR is non-negative definite, and hence that there
%  are no negative eigenvalues.  If this is not the case,
%  warn the user, hope the numbers are only slightly negative,
%  and reset them to 0.
%
  dmin = min ( min ( d ) );

  if ( dmin < - sqrt ( eps ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SAMPLE_PATHS_EIGEN - Warning!\n' );
    fprintf ( 1, '  Negative eigenvalues observed as low as %g\n', dmin );
  end

  d = max ( d, 0.0 );
%
%  Compute the eigenvalues of the factor C.
%
  sqrt_d = sqrt ( d );
%
%  Compute C, such that C' * C = COR.
%
  c = v * sqrt_d * v';
%
%  Compute N independent random normal values.
%
  r(1:n,1:n2) = randn ( n, n2 );
%
%  Get the variables X which have correlation COR.
%
  x(1:n,1:n2) = c(1:n,1:n) * r(1:n,1:n2);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = cumsum( x, dim )

%Disciplined convex/geometric programming information for SUM:
%   CUMSUM(X) and CUMSUM(X,DIM) are vectorized forms of addition. So
%   when CUMSUM is used in a DCP or DGP, elements in each subvector
%   must satisfy the corresponding combination rules for addition (see
%   PLUS). For example, suppose that X looks like this:
%      X = [ convex concave affine  ;
%            affine concave concave ]
%   Then CUMSUM(X,1) would be permittted, but CUMSUM(X,2) would not,
%   because the top row contains the sum of convex and concave terms, in
%   violation of the DCP ruleset. For DGPs, addition rules dictate that
%   the elements of X must be log-convex or log-affine.

s = x.size_;
switch nargin,
    case 0,
        error( 'Not enough input arguments.' );
    case 1,
        dim = cvx_default_dimension( s );
    case 2,
        if ~cvx_check_dimension( dim, false ),
            error( 'Second argument must be a dimension.' );
        end
end

if dim > length( s ) || s( dim ) <= 1,

    y = x;

else

    b = x.basis_;
    sb = size( b );
    need_perm = any( s( 1 : dim - 1 ) > 1 );
    if need_perm,
        ndxs = reshape( 1 : prod( s ), s );
        ndxs = permute( ndxs, [ dim, 1 : dim - 1, dim + 1 : length( s ) ] );
        b = b( :, ndxs );
    end
    b = reshape( b, prod( sb ) / s( dim ), s( dim ) );
    b = cumsum( b, 2 );
    b = reshape( b, sb );
    if need_perm,
        b( :, ndxs ) = b;
    end
    y = cvx( s, b );
    v = cvx_vexity( y );
    if any( isnan( v( : ) ) ),
        error( 'Disciplined convex programming error:\n   Illegal addition encountered (e.g., {convex} + {concave}).', 1 ); %#ok
    end

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:
classdef CalBatteryAndThermometerDataClass < BatteryAndThermometerDataBaseClass

    %% Public 'read-only' properties
    properties (SetAccess = private)
        FileNameAppendage = '_CalBattAndTherm.csv';
    end

    %% Public methods
    methods (Access = public)
        function obj = CalBatteryAndThermometerDataClass(varargin)
            fileNamePrefix = varargin{1};
            for i = 2:2:nargin
                if  strcmp(varargin{i}, 'SampleRate'), obj.SampleRate = varargin{i+1};
                else error('Invalid argument.');
                end
            end
            obj.Import(fileNamePrefix);

            % Set protected parent class variables
            obj.ThermometerUnits = '^\circC';
            obj.BatteryUnits = 'G';
        end
    end
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
function z = rk4_cannon(t,z0,P)
% z = rk4_cannon(dyn,t,z0,P)
%
% This function is used to perform a 4th-order Runge-Kutta
% integration of the cannon dynamical system
%
% INPUTS:
%   t = [1 x nTime] vector of times, created by linspace
%   z0 = [nState (x nSim)] matrix of initial states
%   P = parameter to pass to dynamics
%
% OUTPUTS:
%   z = [nState (x nSim) x nTime ] matrix of trajectories
%
% NOTES:
%   1) See rk4.m for the generic version of this function
%   2) I've hard-coded the dynamics function into this version because it
%      is faster than calling the anonymous function with an embedded
%      parameter.
%

[nState,nSim] = size(z0);
nTime = length(t);

if nSim==1
    z = zeros(nState,nTime);
    z(:,1) = z0;
    for i=1:(nTime-1)
        dt = t(i+1)-t(i);
        k1 = cannonDynamics(t(i),  z(:,i),P);
        k2 = cannonDynamics(t(i)+0.5*dt,  z(:,i) + 0.5*dt*k1,P);
        k3 = cannonDynamics(t(i)+0.5*dt,  z(:,i) + 0.5*dt*k2,P);
        k4 = cannonDynamics(t(i)+dt,  z(:,i) + dt*k3,P);
        z(:,i+1) = z(:,i) + (dt/6)*(k1+2*k2+2*k3+k4);
    end
else
    z = zeros(nState,nSim,nTime);
    z(:,:,1) = z0;
    for i=1:(nTime-1)
        dt = t(i+1)-t(i);
        k1 = cannonDynamics(t(i),  z(:,:,i),P);
        k2 = cannonDynamics(t(i)+0.5*dt,  z(:,:,i) + 0.5*dt*k1,P);
        k3 = cannonDynamics(t(i)+0.5*dt,  z(:,:,i) + 0.5*dt*k2,P);
        k4 = cannonDynamics(t(i)+dt,  z(:,:,i) + dt*k3,P);
        z(:,:,i+1) = z(:,:,i) + (dt/6)*(k1+2*k2+2*k3+k4);
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%-------------------------------------------------------------------------
% Data_Analysis.m
% This file is about sample codes.
% If you want to execute this sample codes, you downloaded sampleData.zip.
% Go to website below and download it.
% (https://osf.io/kprvn/?view_only=0bfb915de5bf4927baad3f92125feee3)
%
% Each paragraph is followed by a detailed description.
% All codes can
%-------------------------------------------------------------------------

%% Convert_sample code
% This code coverts raw file to mat file
edit 60EEG_4EOG_7EMG_converter

%% FBCSP_RLDA_sample codes (reaching, grasp and twist)
% This code is about to filter bank common spatial pattern (FBCSP) with
% regularized linear discriminant analysis (RLDA).
% Each line is about reaching, grasp and twist.
edit FBCSP_RLDA_reaching
edit FBCSP_RLDA_grasp
edit FBCSP_RLDA_twist

%% Signal figure_sample codes
% These codes show each of raw signals of electroencephalogram (EEG), electrooculography (EOG) and electromyography (EMG) figure.
edit figureEEG
edit figureEOG
edit figureEMG

%% Scalp plot_sample codes
%
edit plotScalp_session1_SampleCode
edit plotScalp_session2_SampleCode
edit plotScalp_session3_SampleCode

%% ERSP_sample codes
% These codes show event-related spectral dynamics (ERSP) of each sessions.
edit ERSP_session1_SampleCode
edit ERSP_session2_SampleCode
edit ERSP_session3_SampleCode

author: PatternRecognition

repo: OpenBMI


text:
%  crs_linear_filter_update_cell_estim_K -- filter with update of K. gain
%
%  ::
%
%
%    [loglik,Incr,retcode,Filters]=crs_linear_filter_update_cell_estim_K(...
%     syst,y,U,z,options)
%
%  Args:
%
%     - **syst** [struct]: structure containing:
%
%           - **PAI00** [vector]: initial probability distributions of regimes
%
%           - **a** [cell]: initial conditions in each regime
%
%           - **Qfunc** [function handle]: transition matrix generator
%
%           - **ff** [function handle]: ft=ff(rt,xt,et), where rt is the
%           regime, xt is the vector of state variables and et the vector of
%           shocks
%
%           - **P** [cell]: initial covariance matrix of the states in each
%           regime
%
%           - **H** [cell]: Measurement error covariance matrices in each regime
%
%           - **SIGeta** [cell]: Covariance matrix of structural shocks.
%
%     - **y** [matrix]: ny x T x npages matrix of data
%
%     - **U** [[]|matrix]: ndx x T matrix of exogenous data
%
%     - **z** [function handle|logical|vector]: linear connection of the
%     observables to the state.
%
%     - **include_in_likelihood** [logical]: selector of increments to include
%     in the likelihood calculation
%
%     - **options** [struct]: structure with various options
%
%  Returns:
%     :
%
%     - **loglik** [scalar]: log likelihood
%
%     - **Incr** [vector]: increments of elements going into the likelihood
%
%     - **retcode** [{0}|integer]: flag for problems.
%
%     - **Filters** [struct]: Filtered, updated and smoothed variables
%
%  Note:
%
%     - The filter checks for violations of constraints and uses the
%     information in the database to reset the offending variables to their
%     values in the database. When this occurs, the Kalman gain is recomputed
%     so that the traditional updating equation holds.
%
%     - This strategy is adopted so as to permit the use of the efficient
%     smoothing algorithm of Durbin and Koopman instead of the classical
%     smoothing algorithm which requires multiple inversions of a potentially
%     singular covariance matrix.
%
%  Example:
%
%     See also:
%

author: jmaih


text:
% Version 1.000
%
% Code provided by Ruslan Salakhutdinov and Geoff Hinton
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

function [f, df] = CG_MNIST(VV,Dim,XX);

l1 = Dim(1);
l2 = Dim(2);
l3 = Dim(3);
l4= Dim(4);
l5= Dim(5);
l6= Dim(6);
l7= Dim(7);
l8= Dim(8);
l9= Dim(9);
N = size(XX,1);

% Do decomversion.
 w1 = reshape(VV(1:(l1+1)*l2),l1+1,l2);
 xxx = (l1+1)*l2;
 w2 = reshape(VV(xxx+1:xxx+(l2+1)*l3),l2+1,l3);
 xxx = xxx+(l2+1)*l3;
 w3 = reshape(VV(xxx+1:xxx+(l3+1)*l4),l3+1,l4);
 xxx = xxx+(l3+1)*l4;
 w4 = reshape(VV(xxx+1:xxx+(l4+1)*l5),l4+1,l5);
 xxx = xxx+(l4+1)*l5;
 w5 = reshape(VV(xxx+1:xxx+(l5+1)*l6),l5+1,l6);
 xxx = xxx+(l5+1)*l6;
 w6 = reshape(VV(xxx+1:xxx+(l6+1)*l7),l6+1,l7);
 xxx = xxx+(l6+1)*l7;
 w7 = reshape(VV(xxx+1:xxx+(l7+1)*l8),l7+1,l8);
 xxx = xxx+(l7+1)*l8;
 w8 = reshape(VV(xxx+1:xxx+(l8+1)*l9),l8+1,l9);

  XX = [XX ones(N,1)];
  w1probs = 1./(1 + exp(-XX*w1)); w1probs = [w1probs  ones(N,1)];
  w2probs = 1./(1 + exp(-w1probs*w2)); w2probs = [w2probs ones(N,1)];
  w3probs = 1./(1 + exp(-w2probs*w3)); w3probs = [w3probs  ones(N,1)];
  w4probs = w3probs*w4; w4probs = [w4probs  ones(N,1)];
  w5probs = 1./(1 + exp(-w4probs*w5)); w5probs = [w5probs  ones(N,1)];
  w6probs = 1./(1 + exp(-w5probs*w6)); w6probs = [w6probs  ones(N,1)];
  w7probs = 1./(1 + exp(-w6probs*w7)); w7probs = [w7probs  ones(N,1)];
  XXout = 1./(1 + exp(-w7probs*w8));

f = -1/N*sum(sum( XX(:,1:end-1).*log(XXout) + (1-XX(:,1:end-1)).*log(1-XXout)));
IO = 1/N*(XXout-XX(:,1:end-1));
Ix8=IO;
dw8 =  w7probs'*Ix8;

Ix7 = (Ix8*w8').*w7probs.*(1-w7probs);
Ix7 = Ix7(:,1:end-1);
dw7 =  w6probs'*Ix7;

Ix6 = (Ix7*w7').*w6probs.*(1-w6probs);
Ix6 = Ix6(:,1:end-1);
dw6 =  w5probs'*Ix6;

Ix5 = (Ix6*w6').*w5probs.*(1-w5probs);
Ix5 = Ix5(:,1:end-1);
dw5 =  w4probs'*Ix5;

Ix4 = (Ix5*w5');
Ix4 = Ix4(:,1:end-1);
dw4 =  w3probs'*Ix4;

Ix3 = (Ix4*w4').*w3probs.*(1-w3probs);
Ix3 = Ix3(:,1:end-1);
dw3 =  w2probs'*Ix3;

Ix2 = (Ix3*w3').*w2probs.*(1-w2probs);
Ix2 = Ix2(:,1:end-1);
dw2 =  w1probs'*Ix2;

Ix1 = (Ix2*w2').*w1probs.*(1-w1probs);
Ix1 = Ix1(:,1:end-1);
dw1 =  XX'*Ix1;

df = [dw1(:)' dw2(:)' dw3(:)' dw4(:)' dw5(:)' dw6(:)'  dw7(:)'  dw8(:)'  ]';

author: qiuwch

repo: DeepLearning

save_path: github-repos/MATLAB/qiuwch-DeepLearning


text:
function [beta_median,beta_std,beta_lbound,beta_ubound,sigma_median,beta_mean_median,beta_mean_lbound,beta_mean_ubound,sigma_mean_median]...
    =panel4estimates(N,n,q,beta_gibbs,sigma_gibbs,cband, beta_mean,sigma_mean)

% compute first the  percentiles for the mean model values, i.e. beta_mean
% and sigma_man
%
% for jj=1:n^2
%  sigma_mean_median(jj,1)=[quantile(sigma_mean(jj,:),0.5)];
% end
%    for jj=1:q
%    beta_mean_median(jj,1)=[quantile(beta_mean(jj,:),0.5)];
%    beta_mean_lbound(jj,1)=[quantile(beta_gibbs(jj,1),(1-cband)/2)];
%    beta_mean_ubound(jj,1)=[quantile(beta_gibbs(jj,1),1-(1-cband)/2)];
%    end
% as the VAR coefficients are estimated for each unit, loop over units
% move to the country specific coefficients
for ii=1:N
   % loop over VAR coefficients
   for jj=1:q
   beta_median(jj,1,ii)=[quantile(beta_gibbs(jj,:,ii),0.5)];
   beta_std(jj,1,ii)=std(beta_gibbs(jj,:,ii));
   beta_lbound(jj,1,ii)=[quantile(beta_gibbs(jj,1,ii),(1-cband)/2)];
   beta_ubound(jj,1,ii)=[quantile(beta_gibbs(jj,1,ii),1-(1-cband)/2)];
   end
   % loop over sigma entries
   for jj=1:n^2
   sigma_median(jj,1,ii)=[quantile(sigma_gibbs(jj,:,ii),0.5)];
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function setOptsDoubleValueStrInGUI(handles, options, optPropName, handlesFieldName)
    value = options.(optPropName);

    if(isnan(value))
        str = '';
    else
        str = fullAccNum2Str(double(value));
    end

    handles.(handlesFieldName).String = str;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function reg_filter = get_reg_filter(sz, target_sz, params, reg_window_edge)

% Compute the spatial regularization function and derive the corresponding
% filter operation used for the optimization

if nargin < 3 || isempty(reg_window_edge)
    reg_window_edge = params.reg_window_edge;
end

if params.use_reg_window
    % create weight window
    reg_window_power = params.reg_window_power;

    % normalization factor
    reg_scale = 0.5 * target_sz;

    % construct grid
    wrg = -(sz(1)-1)/2:(sz(1)-1)/2;
    wcg = -(sz(2)-1)/2:(sz(2)-1)/2;
    [wrs, wcs] = ndgrid(wrg, wcg);

    % construct the regukarization window
    reg_window = (reg_window_edge - params.reg_window_min) * (abs(wrs/reg_scale(1)).^reg_window_power + abs(wcs/reg_scale(2)).^reg_window_power) + params.reg_window_min;

    % compute the DFT and enforce sparsity
    reg_window_dft = fft2(reg_window) / prod(sz);
    reg_window_dft(abs(reg_window_dft) < params.reg_sparsity_threshold * max(abs(reg_window_dft(:)))) = 0;

    % do the inverse transform, correct window minimum
    reg_window_sparse = real(ifft2(reg_window_dft));
    reg_window_dft(1,1) = reg_window_dft(1,1) - prod(sz) * min(reg_window_sparse(:)) + params.reg_window_min;
    reg_window_dft = fftshift(reg_window_dft);

    % find the regularization filter by removing the zeros
    reg_filter = single(real(reg_window_dft(~all(reg_window_dft==0,2), ~all(reg_window_dft==0,1))));
else
    % else use a scaled identity matrix
    reg_filter = single(params.reg_window_min);
end

author: martin-danelljan

repo: Continuous-ConvOp


text:
function poseEst = refinePose2D(edges, thEst, posesInit)

isDebug = 0;

m = size(edges,1);
nrNodes = length(thEst); % thEst is the orientation estimate
n = nrNodes-1; % number of observable nodes

J = sparse(spalloc(4*m, 3*nrNodes, 6*4*m)); % f(x) in (4*m) x 3*nrNodes
% since each measurement includes a Cartesian measurement in R^2 and
% a rotation measurement (converted to R^2 vector). The unknown x contains
% 3 variables per node (x,y,th). Each row of the Jacobian involves 2 nodes
% and can have at most 6 nonzero quantities
fhat = sparse(zeros(4*m,1));
Mij = [-1 0; 0 -1];
cost = 1e+20; % initial cost set to infinite
% Initial guess
if  nargin < 3
  poseEst = [zeros(2*nrNodes,1); thEst];
else
  poseEst = [reshape(posesInit(:,1:2)',2*nrNodes,1); posesInit(:,3)]; % written as a vector, with xy positions first, then theta
end

nonAnchorInds = [3:2*nrNodes,2*nrNodes+2:3*nrNodes];
JR = spalloc(4*m, nrNodes, 2*4*m); % 2 nonzero elements per row
JC = spalloc(4*m, 2*nrNodes, 4*4*m); % 4 nonzero elements per row

for iter=1:10
  %% linerize
  costPrec = cost;
  for k=1:m
    tic
    id1 = edges(k,1);
    id2 = edges(k,2);
    Ri = rot2D(poseEst(2*nrNodes+id1));
    Rj = rot2D(poseEst(2*nrNodes+id2));
    Rij = rot2D(edges(k,5));
    Deltaij = edges(k,3:4)';
    DeltaijMat = [Deltaij(1)  -Deltaij(2); Deltaij(2)  Deltaij(1)];

    id1c = blockToMatIndices(id1,2);
    p1 = poseEst(id1c);
    id2c = blockToMatIndices(id2,2);
    p2 = poseEst(id2c);

    fcart = p2 - p1 - DeltaijMat * Ri(:,1);
    frot  = Rij * Ri(:,1) + Mij * Rj(:,1);

    %% here we do not need the 1/2, since is [cos sin] instead of rotation
    inf_th = sqrt ( edges(k,11) ); % inverse of std for rotation measurements
    inf_cart = sqrt ( edges(k,6) ); % inverse of std for position measurements

    % each measurement f in R^4 is f=[fcart frot]
    rowInd = blockToMatIndices(k,4);
    fhat(rowInd) = [inf_cart * fcart; inf_th * frot];

    JR(rowInd,id1) = [ - inf_cart * DeltaijMat * Ri(:,2);
                         inf_th * Rij * Ri(:,2) ];
    JR(rowInd,id2) = [ zeros(2,1)          ;
                       inf_th * Mij * Rj(:,2) ]; % thj not in cartesian measurements

    JC(rowInd,id1c) = [-inf_cart * eye(2) ;
                        zeros(2,2)];
    JC(rowInd,id2c) = [ inf_cart * eye(2) ;
                        zeros(2,2)];
  end
  %% current cost
  cost = norm(fhat)^2;

  %% delete anchor
  Jtot = [JC(:,3:end) JR(:,2:end)]; % cartesian part first
  corrEst = - (Jtot'*Jtot) \ (Jtot' * fhat);% minus because of fhat

  %% update estimate
  poseEst(nonAnchorInds) = poseEst(nonAnchorInds) + corrEst; % update estimate, preserving anchors

  %% check stopping condition
  relCostChange = abs((cost - costPrec)/costPrec);
  if (isDebug>0)
    fprintf('Current cost: %f, norm of the correction: %f, relative decrease: %.10f \n', cost, norm(corrEst), relCostChange);
  end
  if norm(corrEst) < 1e-4 || relCostChange < 1e-5
    break;
  end
end

if (isDebug==2)
  figure
  plot(poseEst(1:2:2*nrNodes),poseEst(2:2:2*nrNodes),'-k');
  title('Estimate from refinePose2D')
end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
% example8_10.m
net1=competlayer(40);
load simpleclass_dataset
net2=newsom(simpleclassInputs,[5 8]);
view(net1)
view(net2)

web -broswer http://www.ilovematlab.cn/forum-222-1.html

author: HuangCongQing


text:
function [boutmat,line_names] = ethogram_plot_core(hfig,labels,T0,T1,expdirs,behaviornames,jabfiles,varargin)
% [boutmat,line_names] = ethogram_plot_core(hfig,labels,expdirs,T0,T1,behaviornames,jabfiles,p1,v1,...)
%
% labels: nexps x nsamps x nbehs label array. 0s and 1s are plotted as
% none and behavior-present, respectively. -1's may also be supplied, which
% are plotted as negative-behavior (faded color).
%
% Optional PVs:
%
% Input arg jabfiles are used only to enable click-on-bout-to-open-in-JAABA.

[trialdividers,pulldownselect,behaviorcolorslines] = ...
  myparse(varargin,...
  'trialdividers',false,...
  'pulldownselect',false,...
  'behaviorcolorslines',false);

figure(hfig);
clf;
ax = axes;

nexps = numel(expdirs);
nT = T1-T0+1;
nbehaviors = numel(behaviornames);
assert(size(labels,1)==nexps);
assert(size(labels,2)==nT);
assert(size(labels,3)==nbehaviors);
assert(numel(jabfiles)==nbehaviors,...
  'Input arg ''jabfiles'' must be a cellstr specifying the containing jab for each behavior');

if behaviorcolorslines || nbehaviors>7
  behaviorcolors = lines(nbehaviors);
else
  BEH_COLORS = [
         0         0    1.0000; ...
         0    0.5000         0; ...
    1.0000         0         0; ...
         0    0.7500    0.7500; ...
    0.7500         0    0.7500; ...
    0.7500    0.6000         0; ...
    0.2500    0.2500    0.2500];
  behaviorcolors = BEH_COLORS;
end
nobehaviorcolors = behaviorcolors + 0.67*(1-behaviorcolors); % gray transform
yrad = 1/nbehaviors;
xl = [T0-10,T1+10];

h = cell(nexps,nbehaviors); % patch handles; we keep track of these but ends up unused
hNo = cell(nexps,nbehaviors);
hlegend = nan(1,nbehaviors);
boutmat = cell(nexps,nbehaviors);
isfirst = true;
for j = 1:nbehaviors,
  behclr = behaviorcolors(j,:);
  nobehclr = nobehaviorcolors(j,:);
  for i = 1:nexps,

    if trialdividers && j==1
      % plot dividers for each exp
      plot(xl,repmat(i-yrad/2,1,2),'k:');
    end
    if isfirst,
      isfirst = false;
      hold on;
    end

    [starts,ends] = get_interval_ends(labels(i,:,j)>0);
    boutmat{i,j} = [starts(:) ends(:)];

    % pos behaviors
    if numel(starts)>0
      h{i,j} = hlpPatch(starts,ends,i,j,yrad,T0,behclr,{'FaceAlpha',0.4});
      if isnan(hlegend(j)),
        hlegend(j) = h{i,j}(1);
      end
    end

    % negative behaviors
    [starts,ends] = get_interval_ends(labels(i,:,j)<0);
    if numel(starts)>0
      hNo{i,j} = hlpPatch(starts,ends,i,j,yrad,T0,nobehclr,{'FaceAlpha',0.4});
    end

  end
end

line_names = cell(1,nexps);
for i = 1:nexps,
  [~,line_names{i}] = myfileparts(expdirs{i});
end
legend(hlegend(~isnan(hlegend)),behaviornames(~isnan(hlegend)));
yoffset = (nbehaviors-1)/nbehaviors/2;
set(ax,'YTick',yoffset+(1:nexps),'YTickLabel',line_names,'XLim',xl,'YLim',[1,nexps+1],'YDir','reverse','Box','off','TickDir','out');
if isprop(ax,'TickLabelInterpreter')
  ax.TickLabelInterpreter = 'none';
end

% set up axis guidata; data used in callbacks
gdata = struct();
gdata.T0 = T0;
gdata.T1 = T1;
gdata.labels = labels;
gdata.jabfiles = jabfiles;
gdata.expdirs = expdirs;
gdata.Nexp = numel(expdirs);
gdata.behaviornames = behaviornames;
gdata.Nbeh = numel(behaviornames);
gdata.behaviorcolors = behaviorcolors;
gdata.hAutoMarks = [];

% gdata.EXP_SELECTED_COLOR = [1 1 0];
% gdata.EXP_SELECTED_ALPHA = 1;
gdata.EXPZOOMPLUSMINUS = 1; % when zooming on exp, include +/- this many neighboring exps
% gdata.tfExpSelInProgress = false;
% gdata.tfExpIsUserSelected = false(numel(expdirs),1);
% gdata.hUserSelectionPatches = nan(numel(expdirs),1);
% gdata.customGroupVizInfo = cell(0,1); % customGroupVizInfo{i}.name % name of grp
% customGroupVizInfo{i}.tf  % indicator vec for ith viz group
% customGroupVizInfo{i}.h   % plot handles for ith viz group
% Note: due to arbitrary removal, customGroupVizInfo{i} may be [] for
% arbitrary i.
% gdata.CUSTOMGROUPMODN = 7; % this the "mod" for lines()

% gdata.axGroupStats = []; % second/overlay axis for group stats
guidata(ax,gdata);

if pulldownselect
  set(hfig,'units','pixels');
  set(hfig,'position',[100 100 1000 800],'toolbar','figure'); % adding uicontrols to figure below

  tmpunits = get(ax,'Units');
  set(ax,'Units','normalized');
  axpos = get(ax,'Position');
  set(ax,'Units',tmpunits);

  selections = [{'<all exps>'};expdirs(:)];
  uicontrol('style','popupmenu','String',selections,'fontsize',7,'units','normalized',...
    'position',[axpos(1) 0.955 axpos(3) 0.02],'callback',@(zS,zE)cbkExpZoom(zS,ax));
  uicontrol('style','text','String','Experiment zoom: ','fontweight','bold','horizontalalignment','right',...
    'units','normalized','position',[axpos(1)-0.2 0.95 0.2 0.02],'backgroundcolor',get(hfig,'color'));
end

set(ax,'ButtonDownFcn',@(src,evt)cbkAxBDF(src));

function hPatch = hlpPatch(starts,ends,i,j,yrad,T0,clr,patchArgs)
nbouts = numel(starts);
ends = ends-0.5;
starts = starts-0.5;
x = [starts;starts;ends;ends;starts]+T0-1;
y = i+yrad*[0;1;1;0;0]+yrad*(j-1);
xx = x;
yy = repmat(y,[1 nbouts]);
hPatch = patch(xx,yy,-1*ones(size(xx)),clr,...
  'MarkerFaceColor',clr,'HitTest','off','EdgeColor','none',patchArgs{:});

function cbkExpZoom(src,ax)
gdata = guidata(ax);

val = get(src,'value');
str = get(src,'String');
str = str{val};

expnum = find(strcmp(str,gdata.expdirs));
if isempty(expnum)
  ylim(ax,'auto');
else
  ylim(ax,[expnum-gdata.EXPZOOMPLUSMINUS expnum+gdata.EXPZOOMPLUSMINUS+1]);
end

function cbkAxBDF(ax)

gdata = guidata(ax);
ty = get(ancestor(ax,'figure'),'SelectionType');
switch ty
  case 'extend' % Shift-click
    yl = ylim(ax);
    if diff(yl)<=1.5
      % fully zoomed; zoom out
      ylim(ax,'auto');
    elseif diff(yl)<=gdata.EXPZOOMPLUSMINUS*2+1
      % medium zoom; go to full zoom (single exp)
      pt = get(ax,'currentpoint');
      expnum = floor(pt(1,2));
      ylim(ax,[expnum expnum+1]);
    else
      % go to medium zoom
      pt = get(ax,'currentpoint');
      expnum = floor(pt(1,2));
      ylim(ax,[expnum-gdata.EXPZOOMPLUSMINUS expnum+gdata.EXPZOOMPLUSMINUS+1]);
    end
  case 'alt' % Ctrl-click
%     pt = get(ax,'currentpoint');
%     expnum = floor(pt(1,2));
%     gdata = guidata(ax);
%     if expnum > gdata.Nexp % Can occur if ctrl-click slightly 'below' axis
%       warning('ethogram_plot:noSuchExp','Clicked point beyond range of experiments.');
%     else
%       gdata = lclToggleExpSelection(gdata,expnum);
%     end
%     gdata.tfExpSelInProgress = true;
%     guidata(ax,gdata);
  otherwise
    patchClickCallback(ax);
end

function patchClickCallback(h)
% Check for double click

% AL 20150113: This is not working

persistent chk

dothis = true;
if isempty(chk)
  chk = 1;
  pause(0.5); %Add a delay to distinguish single click from a double click
  if chk == 1
    chk = [];
    return;
  end
else
  dothis = false;
  chk = [];
end

if ~dothis, return; end
% Below executes only on double click

gdata = guidata(h);

% if gdata.scoresnotjabs
%   warning('ethogram_plot_core:nojabs','No jab files specified; clicking on bouts is disabled.');
%   return;
% end

pt = get(h,'currentpoint');
frameno = round(pt(1,1));
if frameno<0 || frameno>size(gdata.labels,2),
  return;
end
expnum = floor(pt(1,2));

% determine behavior number/index so we know which jab to open
behnum = find(gdata.labels(expnum,frameno,:));
if numel(behnum)==0, behnum = 1:numel(gdata.behaviornames); end
if numel(behnum)>1
  [sel,ok] = listdlg('ListSize',[650 220],'Name','Select a Behavior to see in JAABA',...
    'ListString',gdata.behaviornames(behnum),'SelectionMode','single','OKString','Select',...
    'InitialValue',1);
  if ok == 0, return; end
  behnum = behnum(sel);
end
assert(any(behnum==(1:gdata.Nbeh)));
selectedJab = gdata.jabfiles{behnum};

JLabelInterf.openJabExpFrame(selectedJab,gdata.expdirs{expnum},frameno);

author: kristinbranson

repo: JAABA


text:
function out = proj_spectahedron_box(X,r,l,u,eq_flag)
%PROJ_SPECTAHEDRON_BOX computes the orthogonal projection onto the spectahedron {X : Tr(X) = r, X psd, l<=eig(X)<=u} or
%                                                onto the full spectahedron {X : Tr(X)<= r, X psd, l<=eig(X)<=u}
%
%  Usage:
%  out = PROJ_SPECTAHEDRON_BOX(X,[r],[l],[u],[eq_flag])
%  ===========================================
%  Input:
%  X - matrix to be projected
%  r - positive scalar [default: 1]
%  l - lower bound (scalar) [default: 0]
%  u - upper bound (scalar) [default: inf]
%  eq_flag - a flag that determines whether the projection is onto the spectahedron ('eq') or
%                 the full spectahedron ('ineq') [defualt: 'eq']
%  ===========================================
%  Assumptions:
%  X symmetric
%  r > 0
%  ===========================================
%  Output:
%  out - projection matrix

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

%reading the user X and setting defalut values when required.
if (nargin < 1)
    error ('usage: spectahedron_box(X,[r],[l],[u],[eq_flag])') ;
end

if (( nargin < 5) || (isempty(eq_flag)) )
    %eq_flag is not given, setting to default value: true
    eq_flag = 'eq' ;
end

if ((nargin < 4) || (isempty(u)))
    %u is not given, setting to default value: inf
    u=inf ;
end

if ((nargin < 3) || (isempty(l)))
    %l is not given, setting to default value: 0
    l=0 ;
end

if ((nargin < 2) || (isempty(r)))
    %r is not given, setting to default value: 1
    r = 1 ;
end

if ( r <= 0 )
   error('Set is infeasible') ;
end

eps = 1e-10 ; % defalut value for eps : 1e-10
if ((size(X,1) ~= size(X,2)) || (norm( X - X') > eps))
    error('usage: spectahedron_box(X,[r],[eq_flag]) - X should be a symmetric matrix') ;
end

X = 0.5 * (X + X');

[V,D] = eig(X) ;

if (strcmp(eq_flag,'eq') == 1)
    %call proj_simplex with eq
   out = V * diag(proj_hyperplane_box(diag(D),ones(length(X),1),r,l,u)) * V' ;
   out=(out+out')/2;
else
    if (strcmp(eq_flag,'ineq') == 1)
        %call proj_simplex with ineq
        out = V * diag((proj_halfspace_box(diag(D),ones(length(X),1),r,l,u))) * V' ;
        out=(out+out')/2;
    else
        error('usage: spectahedron(X,[r],[eq_flag]) - eq_flag should be either eq or ineq') ;
    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function convertDicom2Vista_UU_UMC(subjectName,SkipVols,TR,SessionName)
% convertDicom2Vista_UU_UMC(subjectName,SkipVols,TR,SessionName)
%
%
% Input:
%   subjectName: subject last name, no capitals
%   SkipVols: number of initial time-frames to remove from the data
%   TR: repetition time of fMRI scanner (should really get from the data)
%   SessionName: name of the experiment
%
% The script takes (3D) nifti files converted usign r2agui and converts
% them to mrVista format.
%
% It is adapted from prepareUCSFParnassusUniversal and relies on spm(5) to
% read the nifti files.
%

%% Some checks

if(isempty(which('spm_dicom_headers')))
    error('Requires spm2+ tools! I recommend spm5');
end
if ~exist('subjectName','var') || isempty(subjectName)
    error('Need the lastname of subject as input');
end
if ~exist('SessionName','var') || isempty(SessionName)
    error('Need SessionName, i.e. what is the experiment called');
end
if ~exist('SkipVols','var') || isempty(SkipVols);
    error('Need SkipVols, i.e. the number of initial time frames to remove from the data');
end
if ~exist('TR','var') || isempty(TR);
    error('Need TR, FIX ME I need to be able to read this from the data');
end

%% Clean up old mrVista structure
fprintf('[%s]:WARNING:Removing old mrVista structures!',mfilename);
if exist('mrSESSION_backup.mat','file')
    delete mrSESSION_backup
end
if exist('mrSESSION.mat','file')
    delete mrSESSION
end
if exist('prepareUUvars.mat','file')
    delete prepareUUvars.mat
end

%% Some defaults
if ~exist('Raw','dir')
    error('Data needs to be stored in directories under Raw directory');
end
RawDir = 'Raw';

%% find scans
d = dir(RawDir);
% find directory names
d = d([d.isdir]);
% remove '.' and '..'
remove = false(numel(d),1);
for n=1:numel(d)
    if strcmp(d(n).name,'.') || strcmp(d(n).name,'..')
        remove(n) = true;
    end
end
dirs = d(~remove);

if isempty(dirs)
    fprintf(1,'[%s]:WARNING:No directories found in %sRaw\n',mfilename,filesep);
    if ~isempty(dir(fullfile(RawDir,'*nii*')))
        dirs(1).name = [];
        fprintf(1,'[%s]:WARNING:But found nifti files continuing...\n',mfilename);
    else
        fprintf(1,'[%s]:WARNING:EXITING.\n',mfilename);
        return
    end
end

%% find functionals
c = 0;
for n=1:numel(dirs)
    % find files
    f = dir(fullfile(RawDir,dirs(n).name,'*.nii'));

    if numel(f)>10
        c = c+1;
        foundScans(c).Action = 'Functional'; %#ok<AGROW>

        hdr = spm_vol(['Raw/' dirs(n).name filesep f(1).name]);
        if c==1
            roughAlignment = hdr.mat;
        end

        % sort files in order
        fileNameOrder=zeros(numel(f),1);
        % hack time number is after (-)
        for ii=1:numel(f)
            sep1 = findstr(f(ii).name,'-');
            sep2 = findstr(f(ii).name,'.');
            fileNameOrder(ii) = str2double(f(ii).name(sep1+1:sep2-1));
        end
        [junk,sIndex] = sort(fileNameOrder); % Because we don't know how DIR might order the output
        f             = f(sIndex); % Apply the sort

        foundScans(c).FileNames = cell(numel(f),1); %#ok<AGROW>
        for ii=1:numel(f);
            foundScans(c).FileNames{ii} = f(ii).name; %#ok<AGROW>
        end
        foundScans(c).DirName   = dirs(n).name;%#ok<AGROW>
        foundScans(c).Sequence  = hdr.descrip;%#ok<AGROW>
        foundScans(c).SkipVols  = SkipVols;%#ok<AGROW>
        foundScans(c).Cycles    = 8;%#ok<AGROW>

        foundScans(c).TR        = TR;%#ok<AGROW>
        foundScans(c).RawDir    = RawDir;%#ok<AGROW>
        foundScans(c).Files     = numel(f);%#ok<AGROW>
        foundScans(c).Slices    = hdr.dim(3);%#ok<AGROW>
        foundScans(c).Volumes   = numel(f);%#ok<AGROW>
    else
        c = c+1;
        foundScans(c).Action = 'Anatomical'; %#ok<AGROW>

        hdr = spm_vol(['Raw/' dirs(n).name filesep f(1).name]);
        if c==1
            roughAlignment = hdr.mat;
        end

        % sort files in order
        fileNameOrder=zeros(numel(f),1);
        % hack time number is after (-)
        for ii=1:numel(f)
            sep1 = findstr(f(ii).name,'-');
            sep2 = findstr(f(ii).name,'.');
            fileNameOrder(ii) = str2double(f(ii).name(sep1+1:sep2-1));
        end
        [junk,sIndex] = sort(fileNameOrder); % Because we don't know how DIR might order the output
        f             = f(sIndex); % Apply the sort

        foundScans(c).FileNames = cell(numel(f),1); %#ok<AGROW>
        for ii=1:numel(f);
            foundScans(c).FileNames{ii} = f(ii).name; %#ok<AGROW>
        end
        foundScans(c).DirName   = dirs(n).name;%#ok<AGROW>
        foundScans(c).Sequence  = hdr.descrip;%#ok<AGROW>
        foundScans(c).SkipVols  = SkipVols;%#ok<AGROW>
        foundScans(c).Cycles    = 8;%#ok<AGROW>

        foundScans(c).TR        = TR;%#ok<AGROW>
        foundScans(c).RawDir    = RawDir;%#ok<AGROW>
        foundScans(c).Files     = numel(f);%#ok<AGROW>
        foundScans(c).Slices    = hdr.dim(3);%#ok<AGROW>
        foundScans(c).Volumes   = numel(f);%#ok<AGROW>
    end
end

disp('Making directories');
mkdir('Inplane');
mkdir('Volume');
mkdir('Gray');

save prepareUUvars

status=mrInitRetFromNifti(foundScans,lower(subjectName),SessionName);

disp('Result of mrInitRetFromNifti');
disp(status);

loadSession;
mrSESSION.alignment = roughAlignment;
saveSession;

return

author: vistalab

repo: vistasoft


text:
function [conflicts] = verify(TC)

%VERIFY Check for property conflicts in the threecomp object.
% [CONFLICTS] = VERIFY(TC) Check property conflicts of a threecomp object. Runs a
% series of tests to ensure that the threecomp object meets minimum
% standards and is suitable for further processing.
%
% The output term CONFLICTS is a list of booleans showing the result of
% various consistency tests. A conflict string of '00000' shows that no
% conflicts were found. The presence of 1's denote failed tests. Specifically:
%   position 1:  start times don't match
%   position 2:  end times don't match
%   position 3:  sample rates don't match
%   position 4:  units don't match
%   position 5:  station names don't match
%   position 6:  insufficient data length. Minimum 25 samples required.
%   position 7:  component tuple not recognized (ZNE, ZRT, Z21)
%
% This function is mostly likely to be called from within other
% threecomp functions.

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

for n = 1:length(TC)
    conflicts{n} = do_one(TC(n),num2str(n));
end
conflicts = reshape(conflicts,size(TC));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process a single threecomp trace tuple
function conflict = do_one(TC,num)

conflict = '0000000';

if any(isempty(TC.traces))
    disp(['problem in element ' num ': one or more empty waveforms']);
    conflict(6) = '1';
    return
end

% CHECK START TIMES
t = get(TC.traces,'START_MATLAB');
if ( t(2)~=t(1) ) || ( t(3)~=t(1) )
    disp(['problem in element ' num ': traces have different start times']);
    conflict(1) = '1';
end;

% CHECK END TIMES
t = get(TC.traces,'END_MATLAB');
if ( t(2)~=t(1) ) || ( t(3)~=t(1) )
    disp(['problem in element ' num ': traces have different end times']);
    conflict(2) = '1';
end;

% CHECK FREQUENCIES
freq = get(TC.traces,'FREQ');
if ( freq(2)~=freq(1) ) || ( freq(3)~=freq(1) )
    disp(['problem in element ' num ': traces have different frequencies']);
    conflict(3) = '1';
end;

% CHECK FOR COMMON UNITS
units = get(TC.traces,'UNITS');
tf = strcmpi(units,units(1));
if ~all(tf)
    disp(['problem in element ' num ': traces have different units']);
    conflict(4) = '1';
end

% CHECK FOR COMMON STATION, NETWORK, LOCATION NAMES
sta = get(TC.traces,'STATION');
tf = strcmpi(sta,sta(1));
if ~all(tf)
    disp(['problem in element ' num ': traces have different station codes.']);
    conflict(5) = '1';
end

% CHECK FOR SUFFICIENT DATA LENGTH
dataLength = get(TC.traces,'DATA_LENGTH');
if min(dataLength)<25
    disp(['problem in element ' num ': traces are too short. Not enough samples.']);
    conflict(6) = '1';
end

% CHECK ARRANGEMENT OF CHANNELS IN WAVEFORM OBJECT
channels = get(TC.traces,'CHANNEL');
if any([ length(channels{1})<3 length(channels{2})<3 length(channels{3})<3 ])
    disp(['problem in element ' num ': channel codes have less than three characters']);
    conflict(7) = '1';
elseif ~strcmpi(channels{1}(3),'Z')
    disp(['problem in element ' num ': First channel character is not Z component']);
    conflict(7) = '1';
elseif ~strcmpi(channels{2}(3),'N') && ~strcmpi(channels{1,2}(end),'R') && ~strcmpi(channels{1,2}(end),'1')
    disp(['problem in element ' num ': channel character is not N, R, or 2 component']);
    conflict(7) = '1';
elseif ~strcmpi(channels{3}(3),'E') && ~strcmpi(channels{1,3}(end),'T') && ~strcmpi(channels{1,3}(end),'2')
    disp(['problem in element ' num ': Third channel character is not E, T, or 1 component']);
    conflict(7) = '1';
end

author: geoscience-community-codes

repo: GISMO


text:
function script_rfcn_VOC0712_ResNet101_OHEM_ss()
% script_rfcn_VOC0712_ResNet101_OHEM_ss()
% RFCN training and testing with OHEM using ResNet101 model and selective
% search proposals
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

clc;
clear mex;
clear is_valid_handle; % to clear init_key
run(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'startup'));
%% -------------------- CONFIG --------------------
opts.caffe_version          = 'caffe_rfcn';
opts.gpu_id                 = auto_select_gpu;
active_caffe_mex(opts.gpu_id, opts.caffe_version);

% model
model                       = Model.ResNet101_for_RFCN_VOC0712_OHEM();
% cache name
opts.cache_name             = 'rfcn_VOC0712_ResNet101_OHEM_ss';
% config
conf                        = rfcn_config_ohem('image_means', model.mean_image);
% train/test data
fprintf('Loading dataset...')
dataset                     = [];
dataset                     = Dataset.voc0712_trainval_ss(dataset, 'train', conf.use_flipped);
dataset                     = Dataset.voc2007_test_ss(dataset, 'test', false);
fprintf('Done.\n');

% do validation, or not
opts.do_val                 = true;

%% -------------------- TRAINING --------------------

opts.rfcn_model        = rfcn_train(conf, dataset.imdb_train, dataset.roidb_train, ...
                                'do_val',           opts.do_val, ...
                                'imdb_val',         dataset.imdb_test, ...
                                'roidb_val',        dataset.roidb_test, ...
                                'solver_def_file',  model.solver_def_file, ...
                                'net_file',         model.net_file, ...
                                'cache_name',       opts.cache_name, ...
                                'caffe_version',    opts.caffe_version);
assert(exist(opts.rfcn_model, 'file') ~= 0, 'not found trained model');

%% -------------------- TESTING --------------------
                          rfcn_test(conf, dataset.imdb_test, dataset.roidb_test, ...
                                'net_def_file',     model.test_net_def_file, ...
                                'net_file',         opts.rfcn_model, ...
                                'cache_name',       opts.cache_name,...
                                'ignore_cache',     true);

end

repo: R-FCN


text:
classdef ProductComputerOutputCreator < handle

    properties (Access = private)
        fourthOrder
        secondOrder
        secondOrderOut
    end

    methods (Access = public)

        function obj = ProductComputerOutputCreator(C,e)
            obj.init(C,e)
            obj.createSecondOrderOut()
        end

        function s = getSecondOrderOut(obj)
            s = obj.secondOrderOut;
        end
    end

    methods (Access = private)

        function init(obj,C,e)
            obj.fourthOrder = C;
            obj.secondOrder = e;
        end

        function createSecondOrderOut(obj)

            if obj.isTensor()
                if obj.is3D()
                    if obj.isStress()
                        s = Stress3DTensor;
                    elseif obj.isStrain()
                        s = Strain3DTensor;
                    end
                elseif obj.isPlaneStress()
                    if obj.isStress()
                        s = StressPlaneStressTensor;
                    elseif obj.isStrain()
                        s = StrainPlaneStressTensor;
                    end
                end
            elseif obj.isVoigt()
                if obj.is3D()
                    if obj.isStress()
                        s = Stress3DVoigtTensor();
                    elseif obj.isStrain()
                        s = Strain3DTensor;
                    end
                elseif obj.isPlaneStress()
                    if obj.isStress()
                        s = StressPlaneStressVoigtTensor;
                    elseif obj.isStrain()
                        s = StrainPlaneStressVoigtTensor;
                    end
                end
            end
            obj.secondOrderOut = s;

        end

        function itIs = isVoigt(obj)
            isFourthVoigt = strcmp(obj.fourthOrder.getRepresentation(),'voigt');
            isSecondVoigt = strcmp(obj.secondOrder.getRepresentation(),'voigt');
            itIs = isFourthVoigt && isSecondVoigt;
        end

        function itIs = isTensor(obj)
            isFourthVoigt = strcmp(obj.fourthOrder.getRepresentation(),'tensor');
            isSecondVoigt = strcmp(obj.secondOrder.getRepresentation(),'tensor');
            itIs = isFourthVoigt && isSecondVoigt;
        end

        function itIs = isStress(obj)
            isConstiutive = strcmp(obj.fourthOrder.getOrder(),'fourth');
            isStrian = strcmp(obj.secondOrder.getFieldName(),'strain');
            itIs = isConstiutive && isStrian;
        end

        function itIs = isStrain(obj)
            isConstiutive = strcmp(obj.fourthOrder.getOrder(),'fourth');
            isStress = strcmp(obj.secondOrder.getFieldName(),'stress');
            itIs = isConstiutive && isStress;
        end

        function itIs = is3D(obj)
            isFourth3D = strcmp(obj.fourthOrder.getElasticityCase(),'3D');
            isSecond3D = strcmp(obj.secondOrder.getElasticityCase(),'3D');
            itIs = isFourth3D && isSecond3D;
        end

        function itIs = isPlaneStress(obj)
            isFourthPS = strcmp(obj.fourthOrder.getElasticityCase(),'planeStress');
            isSecondPS = strcmp(obj.secondOrder.getElasticityCase(),'planeStress');
            itIs = isFourthPS && isSecondPS;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function fx = p33_fun ( n, x )

%*****************************************************************************80
%
%% P33_FUN evaluates the integrand for problem 33.
%
%  Discussion:
%
%    The integrand is singular at both endpoints of the interval.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    sqrt ( - ln ( x ) )
%
%  Exact Integral:
%
%    sqrt ( pi ) / 2
%
%  Approximate Integral (20 digits):
%
%    0.88622692545275801365...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kendall Atkinson,
%    An Introduction to Numerical Analysis,
%    Prentice Hall, 1984, page 307.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = sqrt ( - log ( x ) );

  i = find ( x == 0.0 );
  fx(i) = 0.0;

  return
end


text:
function Ma = timesv(M,a);
%function Ma = timesv(M,a);
%
% Multiplies a by M:
% a(:,:,i) = M*a(:,:,i)
%
% See also FILTER.

dima = size(a,1);
I = eye(dima);
Ma = armafilterv(a,I,M);

author: Sable

repo: mcbench-benchmarks


text:
function [ux,ispos] = psdfactor(x,K)
% [ux,ispos] = psdfactor(x,K)
%
% PSDFACTOR  UX'*UX Cholesky factorization
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%
% disp('The SeDuMi binaries are not installed.')
% disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
% disp('For more information see the file Install.txt.')
% error(' ')

Ks = K.s;
if isempty(Ks),
    ux = zeros(0,1);
    ispos = true;
    return
end
Kq = Ks .* Ks;
nr = K.rsdpN;
nc = length(Ks);
N  = sum(Kq) + sum(Kq(nr+1:end));
ux = zeros(N,1);
xi = length(x) - N;
ui = 0;
for i = 1 : nc,
    ki = Ks(i);
    qi = Kq(i);
    XX = x(xi+1:xi+qi);
    xi = xi + qi;
    if i > nr,
        XX = XX + 1j*x(xi+1:xi+qi);
        xi = xi + qi;
    end
    XX = reshape(XX,ki,ki);
    [XX,flag]=chol(XX,'lower');
    if flag,
        ispos = false;
        return
    end
    XX = XX + tril(XX,-1)';
    ux(ui+1:ui+qi) = real(XX);
    ui = ui + qi;
    if i > nr,
        ux(ui+1:ui+qi) = imag(XX);
        ui = ui + qi;
    end
end
ispos = true;

author: yu-jiang

repo: radpbook

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


text:
function [labels adjlist pE] = trimapGenerateMultipleSegmentations2( image, imsegs, edgeClassifier, ecal, t, size_array )
    adjmat = imsegs.adjmat;
    segimage = imsegs.segimage;

    spFeat = getSuperpixelData_ver2(image, imsegs);
    [edgeFeat adjlist] = getEdgeData_ver2( imsegs, spFeat );
%     spFeat = mcmcGetSuperpixelData( image, imsegs );
%     [edgeFeat adjlist] = mcmcGetEdgeData( imsegs, spFeat );
    pE=test_boosted_dt_mc(edgeClassifier,edgeFeat);
    pE = 1 ./ (1+exp(ecal(1)*pE+ecal(2)));

    nSuperpixel = max(segimage(:));
    % labels = mexMergeAdjacentRegions2( adjlist, pE, nSuperpixel, t );
    labels = mexMergeAdjRegs_Felzenszwalb( adjlist, pE, nSuperpixel, t, size_array );

%     image_lab = rgb2lab( image );
%     bins = [8 16 16];
%     Q = computeQuantMatrix( image_lab, bins );
%     region_hist = computeRegionHist(Q, bins, segimage);
%
%     num_region = max(segimage(:));
%     region_dist = zeros(num_region, num_region);
%
%     ind = find(adjmat);
%     for ix = 1 : length(ind)
%         [x y] = ind2sub([num_region, num_region], ind(ix));
%         region_dist(x, y) = histDist(region_hist(x,:), region_hist(y,:));
%     end
%
%     t = 0.2:0.05:0.8;
%     labels = mexMergeAdjacentRegions(region_dist, t);

    for jx = 1 : size(labels, 2)
        L = labels(:, jx);
        temp_label = unique(L);
        for ix = 1 : length(temp_label)
            idx = find( L == temp_label(ix) );
            labels(idx, jx) = ix;
        end
    end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function DEM_demo_fMRI_PEB
% Demonstration of PEB for multisession spectral DCM studies
%__________________________________________________________________________
% This demonstration routine illustrates the analysis of a multisession
% fMRI study using spectral DCM. Crucially, the between session effects are
% characterised using empirical Bayes and Bayesian model reduction. This
% means that the original session data are only inverted once (at the
% within session level). The resulting posterior estimates and then used to
% make inferences about between session effects (e.g., time or drug
% effects). The basic question addressed in this sort of analysis is where
% between session effects are expressed in terms of connectivity or
% parameters of neuronal fluctuations. These sorts of effects are specified
% in a second level design matrix in the usual way and can be identified
% using Bayesian model reduction.
%
% in this example, we analyse three sessions with a monotonic change in the
% intrinsic (self) connectivity over three sessions. This involves
% decreases in diagonal A parameters at the first two levels of a simple
% three node hierarchy - and an increase at the highest (third) level.
% Physiologically, this corresponds to a decrease in self-inhibition (or
% increase in excitability) in the lower notes for regions, as time goes
% on.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_fMRI_PEB.m 7679 2019-10-24 15:54:07Z spm $

% Simulate timeseries
%==========================================================================
rng('default')

% assume we have data from eight subjects with a TR of two seconds
% -------------------------------------------------------------------------
T  = 256*8;                           % number of observations (scans)
TR = 2;                               % repetition time or timing
t  = (1:T)*TR;                        % observation times
n  = 3;                               % number of regions or nodes

% setup model
% -------------------------------------------------------------------------
options.nonlinear  = 0;
options.two_state  = 0;
options.stochastic = 0;
options.centre     = 1;
options.induced    = 1;

A   = ones(n,n);
B   = zeros(n,n,0);
C   = zeros(n,n);
D   = zeros(n,n,0);
pP  = spm_dcm_fmri_priors(A,B,C,D,options);

% true parameters - a simple hierarchy of three nodes
% -------------------------------------------------------------------------
pP.A = [  0  -.2    0;
         .2    0  -.1
          0   .3    0];
pP.C = eye(n,n);
pP.transit = randn(n,1)/16;

% spectral density of neuronal fluctuations and observation noise (eg)
% -------------------------------------------------------------------------
[Gu,Gn,Hz,dt] = spm_csd_fmri_gu(pP,TR);
Gu            = Gu(:,1,1)*ones(1,n);
Gn            = Gn(:,1,1)*ones(1,n);

% session specific parameters: here, intrinsic connectivity
% -------------------------------------------------------------------------
X    = exp(-(0:2))';                    % between session decay
B    = [  1    0    0;                  % of intrinsic connections
          0    1    0
          0    0   -1];
X    = X - mean(X);                     % detrended session-effects

% simulate response to endogenous fluctuations
%==========================================================================
for s = 1:numel(X)

    % parameters in this session
    % ---------------------------------------------------------------------
    P    = pP;
    P.A  = P.A + B*X(s);

    % integrate states with endogenous fluctuations (u)
    % ---------------------------------------------------------------------
    M.f  = 'spm_fx_fmri';
    M.x  = sparse(n,5);
    U.u  = spm_rand_power_law(Gu,Hz,dt,T);
    U.dt = TR;
    x    = spm_int_J(P,M,U);

    % haemodynamic observer
    % ---------------------------------------------------------------------
    for i = 1:T
        y(i,:) = spm_gx_fmri(spm_unvec(x(i,:),M.x),[],P)';
    end

    % response with observation noise (e)
    % ---------------------------------------------------------------------
    e       = spm_rand_power_law(Gn,Hz,dt,T);
    Y(s).y  = y + e;
    Y(s).dt = TR;
    PP(s)   = P;

end

% show simulated response (last session)
%--------------------------------------------------------------------------
i = 1:128;
spm_figure('Getwin','Figure 1'); clf
subplot(2,2,1)
plot(t(i),U.u(i,:))
title('Endogenous fluctuations','FontSize',16)
xlabel('Time (seconds)')
ylabel('Amplitude')
axis square

subplot(2,2,2), hold off
plot(t(i),x(i,n + 1:end),'c'), hold on
plot(t(i),x(i,1:n)), hold off
title('Hidden states','FontSize',16)
xlabel('Time (seconds)')
ylabel('Amplitude')
axis square

subplot(2,2,3)
plot(t(i),y(i,:),t(i),e(i,:),':')
title('Hemodynamic response and noise','FontSize',16)
xlabel('Time (seconds)')
ylabel('Amplitude')
axis square

% nonlinear system identification (DCM for CSD)
%==========================================================================
dcm.options = options;

dcm.a    = logical(pP.A);
dcm.b    = zeros(n,n,1);
dcm.c    = zeros(n,1);
dcm.d    = zeros(n,n,0);

% add responsse to each session
% -------------------------------------------------------------------------
for s = 1:numel(X)
    DCM{s,1}   = dcm;
    DCM{s,1}.Y = Y(s);
end

% first level inversion - spectral DCM
% =========================================================================
CSD = spm_dcm_peb_fit(DCM);

% show estimates for a single session
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 1');

j  = 1:n^2;
d  = 3;
subplot(2,2,4); hold off
spm_plot_ci(CSD{d}.Ep.A(:),CSD{d}.Cp(j,j)), hold on
bar(PP(d).A(:),1/4), hold off
title('True and MAP connections (Deterministic)','FontSize',16)
axis square

% Parametric empirical Bayes
%==========================================================================
% having inverted every session, we can construct a between session model
% at the second level. We are interested in identifying which parameters
% change according to the session specific effects encoded in the design
% matrix. This design matrix is supplemented with a constant term, creating
% two times the number of parameters at the second level. Using Bayesian
% model reduction, we can then examine models with and without each
% parameter for the  constant (first) and hypothesised (second) explanatory
% variables in the design matrix. Crucially, we want to explain as much  as
% possible using just between session differences. This can be implemented
% by setting the between session covariance of random effects to a
% relatively small value; here the original prior variance divided by 32.
%--------------------------------------------------------------------------
clear M;
beta  = 32;

M.X   = [X.^0 X];                  % between session explanatory variables
M.hE  = 0;                         % prior expectation of log precision
M.hC  = 1/16;                      % prior covariance of precision
M.bE  = CSD{1}.M.pE;               % prior expectations over sessions
M.bC  = CSD{1}.M.pC;               % prior covariance over sessions
M.pC  = CSD{1}.M.pC/beta;          % prior covariance between sessions

field = {'A','a'};                 % parameters of interest
PEB   = spm_dcm_peb(CSD,M,field);  % empirical Bayesian inversion
BMA   = spm_dcm_peb_bmc(PEB);      % Bayesian model reduction and averaging

% overlay true between session effects over BMA density over parameters
%--------------------------------------------------------------------------
j     = spm_find_pC(CSD{1},'A');
p     = spm_vec(B); p = p(j);
subplot(3,2,4), hold on, bar(p,1/4), hold off

return

author: spm


text:
function b = reduceboxes(model, bs)
% Eliminate columns for filters that are not used.
%   b = reduceboxes(model, bs)
%
%   E.g., [0 0 0 0 10 20 110 120] -> [10 20 110 120]
%   Index end-1 is the component label and index end is the
%   detection score.
%
%   This function assumes that model is a mixture model where
%   each component always places exactly the same number of filters.
%
% Return value
%   b       Filter bounding boxes with unused filter columns removed
% Arguments
%   model   Object model
%   bs      Filter bounding boxes

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

% Only reduce boxes for mixtures of star models
if model.type ~= model_types.MixStar
  b = bs;
  return;
end

% n = #filters per component (assuming all components have
% the same number of parts)
n = length(model.rules{model.start}(1).rhs);
% n*4+2 := 4 coordinates per boxes plus the component index
% and score
b = zeros(size(bs, 1), n*4+2);
maxc = max(bs(:,end-1));
for i = 1:maxc
  % process boxes for component i
  I = find(bs(:,end-1) == i);
  tmp = bs(I,:);
  del = [];
  % find unused filters
  for j = 1:4:size(bs, 2)-2
    % count # of non-zero coordinates
    s = sum(sum(tmp(:,j:j+3)~=0));
    % the filter was not used if all coordinates are zero
    if s == 0
      del = [del j:j+3];
    end
  end
  % remove all unused filters
  tmp(:,del) = [];
  b(I,:) = tmp;
end

author: rbgirshick

repo: voc-dpm

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