text:
function varargout = gui_gridOmori(varargin)
    %    Application M-file for gui_gridOmori.fig
    %    FIG = SV_GRID launch sv_grid GUI.
    %    SV_GRIDMc('callback_name', ...) invoke the named callback.

    % Last Modified by GUIDE v2.0 22-Aug-2002 10:30:09

    if nargin == 1  % LAUNCH GUI

        fig = openfig(mfilename,'reuse');

        % Use system color scheme for figure:
        set(fig,'Color',get(groot,'defaultUicontrolBackgroundColor'));

        % Generate a structure of handles to pass to callbacks, and store it.
        handles = guihandles(fig);
        guidata(fig, handles);

        if varargin{1} == 1
            set(handles.lblSpacingDepth, 'String', 'Latitude spacing [deg]:');
            set(handles.txtSpacingDepth, 'String', '0.1');
            set(handles.lblSpacingHorizontal, 'String', 'Longitude spacing [deg]:');
            set(handles.txtSpacingHorizontal, 'String', '0.1');
        else
            set(handles.lblSpacingDepth, 'String', 'Depth spacing [km]:');
            set(handles.txtSpacingDepth, 'String', '1');
            set(handles.lblSpacingHorizontal, 'String', 'Horizontal spacing [km]:');
            set(handles.txtSpacingHorizontal, 'String', '1');
        end
        if get(handles.radNumber, 'Value') == 1
            set(handles.txtNumber, 'Enable', 'on');
            set(handles.radRadius, 'Value', 0);
            set(handles.txtRadius, 'Enable', 'off');
        else
            set(handles.txtNumber, 'Enable', 'off');
            set(handles.radRadius, 'Value', 1);
            set(handles.txtRadius, 'Enable', 'on');
        end

        % if get(handles.radSplitTime, 'Value') == 1
        %     set(handles.txtSplitTime, 'Enable', 'on');
        %     set(handles.txtTimePeriod, 'Enable', 'on');
        % else
        %     set(handles.radSplitTime, 'Value', 1);
        %     set(handles.txtSplitTime, 'Enable', 'on');
        %     set(handles.txtTimePeriod, 'Enable', 'on');
        % end
        % Wait for callbacks to run and window to be dismissed:
        uiwait(fig);

        if nargout > 0
            varargout{1} = fig;
        end

    elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

        try
            if (nargout)
                [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
            else
                feval(varargin{:}); % FEVAL switchyard
            end
        catch
            disp(lasterr);
        end

    end

    %| ABOUT CALLBACKS:
    %| GUIDE automatically appends subfunction prototypes to this file, and
    %| sets objects' callback properties to call them through the FEVAL
    %| switchyard above. This comment describes that mechanism.
    %|
    %| Each callback subfunction declaration has the following form:
    %| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
    %|
    %| The subfunction name is composed using the object's Tag and the
    %| callback type separated by '_', e.g. 'slider2_Callback',
    %| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
    %|
    %| H is the callback object's handle (obtained using GCBO).
    %|
    %| EVENTDATA is empty, but reserved for future use.
    %|
    %| HANDLES is a structure containing handles of components in GUI using
    %| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
    %| structure is created at GUI startup using GUIHANDLES and stored in
    %| the figure's application data using GUIDATA. A copy of the structure
    %| is passed to each callback.  You can store additional information in
    %| this structure at GUI startup, and you can change the structure
    %| during callbacks.  Call guidata(h, handles) after changing your
    %| copy to replace the stored original so that subsequent callbacks see
    %| the updates. Type "help guihandles" and "help guidata" for more
    %| information.
    %|
    %| VARARGIN contains any extra arguments you have passed to the
    %| callback. Specify the extra arguments by editing the callback
    %| property in the inspector. By default, GUIDE sets the property to:
    %| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
    %| Add any extra arguments after the last argument, before the final
    %| closing parenthesis.

    % --------------------------------------------------------------------
    % Radiobutton: Number of events
function varargout = radNumber_Callback(h, eventdata, handles, varargin)

    if get(handles.radNumber, 'Value') == 1
        set(handles.txtNumber, 'Enable', 'on');
        set(handles.txtMaxRadius, 'Enable', 'on');
        set(handles.radRadius, 'Value', 0);
        set(handles.txtRadius, 'Enable', 'off');
        set(handles.radRectangle, 'Value', 0);
        set(handles.txtSizeRectHorizontal, 'Enable', 'off');
        set(handles.txtSizeRectDepth, 'Enable', 'off');
    else
        set(handles.txtNumber, 'Enable', 'off');
        set(handles.txtMaxRadius, 'Enable', 'off');
        set(handles.radRadius, 'Value', 1);
        set(handles.txtRadius, 'Enable', 'on');
        set(handles.radRectangle, 'Value', 0);
        set(handles.txtSizeRectHorizontal, 'Enable', 'off');
        set(handles.txtSizeRectDepth, 'Enable', 'off');
    end

    % --------------------------------------------------------------------
    % Radiobutton: Constant radius
function varargout = radRadius_Callback(h, eventdata, handles, varargin)

    if get(handles.radRadius, 'Value') == 1
        set(handles.txtRadius, 'Enable', 'on');
        set(handles.radNumber, 'Value', 0);
        set(handles.txtMaxRadius, 'Enable', 'off');
        set(handles.txtNumber, 'Enable', 'off');
        set(handles.radRectangle, 'Value', 0);
        set(handles.txtSizeRectHorizontal, 'Enable', 'off');
        set(handles.txtSizeRectDepth, 'Enable', 'off');
    else
        set(handles.txtRadius, 'Enable', 'off');
        set(handles.txtMaxRadius, 'Enable', 'on');
        set(handles.radNumber, 'Value', 1);
        set(handles.txtNumber, 'Enable', 'on');
        set(handles.radRectangle, 'Value', 0);
        set(handles.txtSizeRectHorizontal, 'Enable', 'off');
        set(handles.txtSizeRectDepth, 'Enable', 'off');
    end

    % --------------------------------------------------------------------
    % Radiobutton: Rectangle size
function varargout = radRectangle_Callback(h, eventdata, handles, varargin)

    if get(handles.radRectangle, 'Value') == 1
        set(handles.txtSizeRectHorizontal, 'Enable', 'on');
        set(handles.txtSizeRectDepth, 'Enable', 'on');
        set(handles.radRadius, 'Value', 0);
        set(handles.txtRadius, 'Enable', 'off');
        set(handles.radNumber, 'Value', 0);
        set(handles.txtNumber, 'Enable', 'off');
        set(handles.txtMaxRadius, 'Enable', 'off');
    else
        set(handles.radRadius, 'Value', 0);
        set(handles.txtRadius, 'Enable', 'off');
        set(handles.radNumber, 'Value', 1);
        set(handles.txtNumber, 'Enable', 'on');
        set(handles.txtMaxRadius, 'Enable', 'on');
    end

    % --------------------------------------------------------------------
    % Radiobutton: Split Time
function varargout = radSplitTime_Callback(h, eventdata, handles, varargin)

    if get(handles.radSplitTime, 'Value') == 1
        set(handles.txtSplitTime, 'Enable', 'on');
        set(handles.txtTimePeriod, 'Enable', 'on');
    else
        set(handles.txtTimePeriod, 'Enable', 'on');
        set(handles.txtSplitTime, 'Enable', 'off');
    end

    % --------------------------------------------------------------------
function varargout = chkTimePeriod_Callback(h, eventdata, handles, varargin)

    if get(handles.chkTimeperiod, 'Value') == 1
        set(handles.txtTimeperiod, 'Enable', 'on');
    else
        set(handles.txtTimeperiod, 'Enable', 'off');
    end

    % --------------------------------------------------------------------
function varargout = chkSaveParameter_Callback(h, eventdata, handles, varargin)

    if get(handles.chkSaveParameter, 'Value') == 1
        set(handles.btnSaveParameter, 'Enable', 'on');
    else
        set(handles.btnSaveParameter, 'Enable', 'off');
    end

    % --------------------------------------------------------------------
function varargout = btnSaveParameter_Callback(h, eventdata, handles, varargin)

    [newfile, newpath] = uiputfile('*.mat','Choose parameter file');
    % Cancel pressed?
    if isequal(newfile, 0)  ||  isequal(newpath, 0)
        return;
    end
    % Everything ok?
    newfile = [newpath newfile];
    if length(newfile) > 1
        set(handles.lblSaveParameter, 'String', newfile);
    end

    % --------------------------------------------------------------------
function varargout = btnOK_Callback(h, eventdata, handles, varargin)

    handles.answer = 1;
    guidata(h, handles);
    uiresume(handles.figGrid);

    % --------------------------------------------------------------------
function varargout = btnCancel_Callback(h, eventdata, handles, varargin)

    handles.answer = 0;
    guidata(h, handles);
    uiresume(handles.figGrid);

author: CelsoReyes


text:
% DEMBARENCOMAP3 Optimise model using MAP approximation with MLP kernel and linear response, using SCG optimisation and new PUMA processed data with noise variance estimation.

% SHEFFIELDML
clear,clc
expNo = 1;
type = 'mapMLPwithNoise';
rand('seed',1);

load demBarenco1;
origModel = model;
clear model
saveFigures = true;

colordef white
[y, yvar, gene, times, scale, rawExp, rawVar] = gpsimLoadBarencoPUMAData;

numGenes = size(gene, 1);
% Get the default options structure.
options = gpsimMapOptions(numGenes);
options.kern = 'mlp';
options.nonLinearity = 'linear';
options.includeNoise = 0;
options.startPoint = 0;                 % TF starts from 0;
options.endPoint = 12;
options.priorProteinTimes = [0]';
options.priorProtein = [0]';
options.includeNoise = 1;
% Modify according to the range of the TF
options.intPoints = (options.endPoint-options.startPoint)/0.1 + 1 ;
options.gParam = [];
options.ngParam = 0;
options.optimiser = 'scg';

options.fix(1).index = 14;         % S of p21
options.fix(1).value = expTransform(1, 'xtoa');
options.fix(2).index = 15;          % D of p21
options.fix(2).value = expTransform(0.8, 'xtoa');
options.fix(3).index = 3;
options.fix(3).value =  expTransform(1, 'xtoa');
% options.fix(2).index = 1;
% options.fix(2).value =  expTransform(1, 'xtoa');;

for i =1:3
    times = times;
%    options.B = origModel.comp{1}.B;
%    options.D = origModel.comp{1}.D; % [0.6417 0.4094 0.4894 0.8000 0.4992];
%    options.S = origModel.comp{1}.S; % [0.4206 0.3255 0.1486 1 0.1870];
    options.S = ones(1, 5);
    options.D = rand(1, 5);
    mu = mean(y{i}, 1);
    options.B = options.D.*mu;
    model.comp{i} = gpsimMapCreate(numGenes, 1, times, y{i}, yvar{i}, options);
    if strcmp(options.kern, 'mlp')
        model.comp{i}.kern.weightVariance = 30;
        model.comp{i}.kern.biasVariance = 1000;
        % This forces kernel recompute.
        params = gpsimMapExtractParam(model.comp{i});
        model.comp{i} = gpsimMapExpandParam(model.comp{i}, params);
    end

end

param = 0;
Nrep = length(model.comp);

for i = 1:Nrep
    paramvec{i} = gpsimMapExtractParam(model.comp{i}); %vector of gamma estimates
    param = param + paramvec{i};
end

param = param/Nrep;

iters = 300;

options = optOptions;
options(14) = iters;
options(9) = 0;
options(1) = 1;

newparam = scg('gpsimMapObjective', param,  options, ...
               'gpsimMapGradients', model);

% [newparam, ll, index] = minimize(param', 'gpsimMapGradFuncWrapper', iters, model);

for rep = 1:Nrep
    options = defaultOptions;
    model.comp{rep} = gpsimMapExpandParam(model.comp{rep}, newparam);
    model.comp{rep} = gpsimMapUpdateF(model.comp{rep}, options);
    model.comp{rep} = gpsimMapUpdateYpredVar(model.comp{rep});
end

type(1) = upper(type(1));
save(['demBarenco' type num2str(expNo)], 'model', 'type', 'expNo', 'scale')
gpsimMapBarencoResults(model, type, expNo, saveFigures, scale)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%  Figure 3.44      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  Example 3.35
clf;
t=[0 .1 .2 .3 .4 .5 1 1.5 2 2.5 3 4 10];
y=[0 .005 .034 .085 .140 .215 .510 .7 .817 .890 .932 .975 .9999];
dy=ones(1,13)-y;
A=-1.37;     % adjust A & alpha to get best fit
alpha=1.;     % to the data by trial and error
axis([0 6 -2 1]);
fity=-A*exp(-alpha*t);
dyl=log10(dy);
fityl=log10(fity);
plot(t,dyl,'o',t,fityl,'-');
title('Fig. 3.44 Plot of log10[y(\infty) - y(t)] versus t');
ylabel('log|1 - y(t)|');
xlabel('Time (sec)');
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
%ADAPTIVETHRESHOLD  Applies an adaptive threshold to an array
%
%     dst = cv.adaptiveThreshold(src)
%     dst = cv.adaptiveThreshold(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Source 8-bit single-channel `uint8` image.
%
% ## Output
% * __dst__ Destination image of the same size and the same type as `src`.
%
% ## Options
% * __MaxValue__ Non-zero value assigned to the pixels for which the condition
%   is satisfied. See the details below. default 255
% * __Method__ Adaptive thresholding algorithm to use, default 'Mean'. One of:
%   * __Mean__ the threshold value `T(x,y)` is a mean of the
%     `BlockSize x BlockSize` neighborhood of `(x,y)` minus `C`
%   * __Gaussian__ the threshold value `T(x,y)` is a weighted sum
%     (cross-correlation with a Gaussian window) of the `BlockSize x BlockSize`
%     neighborhood of `(x,y)` minus `C`. The default sigma (standard deviation)
%     is used for the specified `BlockSize`. See cv.getGaussianKernel
% * __Type__ Thresholding type, default 'Binary'. One of:
%   * __Binary__    `dst(x,y) = (src(x,y) > thresh) ? maxValue : 0`
%   * __BinaryInv__ `dst(x,y) = (src(x,y) > thresh) ? 0 : maxValue`
% * __BlockSize__ Size of a pixel neighborhood that is used to calculate a
%   threshold value for the pixel: 3, 5, 7, and so on. Default 3
% * __C__ Constant subtracted from the mean or weighted mean. Normally, it is
%   positive but may be zero or negative as well. Default 5
%
% The function transforms a grayscale image to a binary image according to
% the formulae:
%
% * __Binary__
%
%                  | maxValue, if src(x,y) > T(x,y)
%       dst(x,y) = |
%                  | 0, otherwise
%
% * __BinaryInv__
%
%                  | 0, if src(x,y) > T(x,y)
%       dst(x,y) = |
%                  | maxValue, otherwise
%
% where `T(x,y)` is a threshold calculated individually for each pixel (see
% `Method` parameter).
%
% Note: Internally the options `'BorderType','Replicate', 'Isolated',true` are
% used to process boundaries (see cv.copyMakeBorder).
%
% See also: cv.threshold, cv.blur, cv.GaussianBlur, adaptthresh, imbinarize
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%
%      JAWAHARLAL NEHRU TECHNOLOGICAL UNIVERSITY
%
%
%
%

clear;
clc;
clear all;
close all;
display('           ');
display('           ');

display('           ');
display('               SOME EXPERIMENTS ON IMAGE DENOISING USING WAVELETS ');

display('           ');
display('           ');
display('         RAJA RAO           ');

display('           ');
display('           ');

display('select the image');

display('           1:lena.png');
display('           2:barbara.png');
display('           3:boat.png');
display('           4:house.png');
display('           5:peppers256.png');
display('           6:cameraman.jpg');
display('           ');
display('           7:hyderabad.png');
display('           8:friendgray.jpg');
display('           ');

ss1=input('enter your choice:   ');
switch ss1
    case 1
       f=imread('lena.png');
       %f=imread('babu.jpg');
    case 2
        f=imread('barbara.png');
    case 3
        f=imread('boat.png');
    case 4
        f=imread('house.png');
        case 5
        f=imread('peppers256.png');
    case 6
        f=imread('cameraman.jpg');
    case 7
        f=imread('hyderabad512.png');
    case 8
        f=imread('friendgray.jpg');
end

subplot(2,2,1), imshow(f);title('original image');

display('enter the type of noise:');
display('    1    for salt & pepper');
display('    2    for gaussian');
display('    3    for poisson');
display('    4    for speckle');

ud=input('enter the value:');

switch ud
    case 1
        display('enter the % of noise(Ex:0.2)');
        ud1=input('pls enter:     ');
        g=imnoise(f,'salt & pepper',ud1);
    case 2

%f=imread('peppers256.png');
%subplot(2,2,1),imshow(f);
display('enter the noise varience:  ');
va=input('enter between 0.01 to 0.09:   ');
g=imnoise(f,'gaussian',0,va);
    case 3
       % display('enter the % of noise(Ex:0.2)');
        %ud1=input('pls enter:     ');
        g=imnoise(f,'poisson');
        case 4
        display('enter the varience of noise(Ex:0.02)');
        ud1=input('pls enter:     ');
        g=imnoise(f,'speckle',ud1);

end
%g=imnoise(f,'salt & pepper',01);
subplot(2,2,2),imshow(g);title('noisy image');

%[ca,ch,cv,cd] = dwt2(g,'db2');
%c=[ca ch;cv cd];
%subplot(2,2,3),imshow(uint8(c));

x=g;
% Use wdencmp for image de-noising.
% find default values (see ddencmp).
[thr,sorh,keepapp] = ddencmp('den','wv',x);
display('');
display('select wavelet');
display('enter 1 for haar wavelet');
display('enter 2 for db2 wavelet');
display('enter 3 for db4 wavelet');
display('enter 4 for sym wavelet');
display('enter 5 for sym wavelet');
display('enter 6 for bior wavelet');
display('enter 7 for bior wavelet');
display('enter 8 for mexh wavelet');
display('enter 9 for coif  wavelet');
display('enter 10 for meyr wavelet');
display('enter 11 for morl wavelet');
display('enter 12 for  rbio wavelet');
display('press any key to quit');
display('');

ww=input('enter your choice:    ');
switch ww
    case 1
        wv='haar';
    case 2
        wv='db2';
    case 3
        wv='db4' ;
    case 4
        wv='sym2'
    case 5
        wv='sym4';
    case 6
        wv='bior1.1';
    case 7
       wv='bior6.8';
    case 8
        wv='mexh';
    case 9
        wv='coif5';
    case 10
        wv='dmey';
    case 11
        wv='mor1';
    case 12
        wv='jpeg9.7';
    otherwise
        quit;
end
display('');
display('enter 1 for soft thresholding');
display('enter 2 for hard thresholding');
display('enter 3 for bayes soft thresholding');
sorh=input('sorh:   ');

display('enter the level of decomposition');
level=input(' enter 1 or 2 :    ');

switch sorh
    case 1
        sorh='s';
        xd = wdencmp('gbl',x,wv,level,thr,sorh,keepapp);
    case 2
        sorh='h';
        xd = wdencmp('gbl',x,wv,level,thr,sorh,keepapp);
    case 3
        %%%%%%%%%%%%%%%%%%%%%
       % clear all;
%close all;
%clc;

%Denoising using Bayes soft thresholding

%Note: Figure window 1 displays the original image, fig 2 the noisy img
%fig 3 denoised img by bayes soft thresholding

%Reading the image
%pic=imread('elaine','png');
pic=f;
%figure, imagesc(pic);colormap(gray);

%Define the Noise Variance and adding Gaussian noise
%While using 'imnoise' the pixel values(0 to 255) are converted to double in the range 0 to 1
%So variance also has to be suitably converted
sig=15;
V=(sig/256)^2;
npic=g;
%npic=imnoise(pic,'gaussian',0,V);
%figure, imagesc(npic);colormap(gray);

%Define the type of wavelet(filterbank) used and the number of scales in the wavelet decomp
filtertype=wv;
levels=level;

%Doing the wavelet decomposition
[C,S]=wavedec2(npic,levels,filtertype);

st=(S(1,1)^2)+1;
bayesC=[C(1:st-1),zeros(1,length(st:1:length(C)))];
var=length(C)-S(size(S,1)-1,1)^2+1;

%Calculating sigmahat
sigmahat=median(abs(C(var:length(C))))/0.6745;

for jj=2:size(S,1)-1
    %for the H detail coefficients
    coefh=C(st:st+S(jj,1)^2-1);
    thr=bayes(coefh,sigmahat);
    bayesC(st:st+S(jj,1)^2-1)=sthresh(coefh,thr);
    st=st+S(jj,1)^2;

    % for the V detail coefficients
    coefv=C(st:st+S(jj,1)^2-1);
    thr=bayes(coefv,sigmahat);
    bayesC(st:st+S(jj,1)^2-1)=sthresh(coefv,thr);
    st=st+S(jj,1)^2;

    %for Diag detail coefficients
    coefd=C(st:st+S(jj,1)^2-1);
    thr=bayes(coefd,sigmahat);
    bayesC(st:st+S(jj,1)^2-1)=sthresh(coefd,thr);
    st=st+S(jj,1)^2;
end

%Reconstructing the image from the Bayes-thresholded wavelet coefficients
bayespic=waverec2(bayesC,S,filtertype);
xd=bayespic;
%Displaying the Bayes-denoised image
%figure, imagesc(uint8(bayespic));colormap(gray);
display('IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 9, NO. 9, SEPTEMBER 2000');

display('IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 9, NO. 9, SEPTEMBER 2000');
display('Adaptive Wavelet Thresholding for Image Denoising and Compression');
display('S. Grace Chang, Student Member, IEEE, Bin Yu, Senior Member, IEEE, and Martin Vetterli, Fellow, IEEE');

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

end

%sorh=sorh;
% de-noise image using global thresholding option.

%f=imread('peppers256.png');
[c,s]=wavefast(g,level,wv);
subplot(2,2,3),wave2gray(c,s,8);title('decomposed structure');

subplot(2,2,4),xd=uint8(xd);
imshow(xd);title('denoised image');
%subplot(2,2,4),sub=f-xd;
%sub=abs(1.2*sub);
%imshow(im2uint8(sub));title('difference image');
ff=im2double(f);xdd=im2double(xd);
display('      ');
display('      ');
display('reference: To calcullate signal to noise ratio');
display('Makoto Miyahara');
display('"Objective Picture Quality Scale (PQS) for Image Coding"');
display('IEEE Trans. on Comm., Vol 46, No.9, 1998.');
display('      ');
display('      ');
snr=wpsnr(ff,xdd)

display('      ');
display('      ');
mse=compare11(ff,xdd)


text:
function varargout = process_psth_per_channel( varargin )
% PROCESS_PSTH_PER_CHANNEL: Computes the PSTH per channel.

% It displays the binned firing rate on each channel (of only the first
% neuron on each channel if multiple have been detected). This can be nicely
% visualized on the cortical surface if the positions of the electrodes
% have been set, and show real time firing rate.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription()
    % Description the process
    sProcess.Comment     = 'PSTH per channel';
    sProcess.FileTag     = 'raster';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Electrophysiology';
    sProcess.Index       = 1229;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/e-phys/functions';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data'};
    sProcess.OutputTypes = {'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Options: Bin size
    sProcess.options.binsize.Comment = 'Bin size: ';
    sProcess.options.binsize.Type    = 'value';
    sProcess.options.binsize.Value   = {0.05, 'ms', 1};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput)
    % Initialize returned values
    OutputFiles = {};

    % ==== OPTIONS =====
    % Bin size
    if isfield(sProcess.options, 'binsize') && ~isempty(sProcess.options.binsize) && ~isempty(sProcess.options.binsize.Value) && iscell(sProcess.options.binsize.Value) && sProcess.options.binsize.Value{1} > 0
        bin_size = sProcess.options.binsize.Value{1};
    else
        bst_report('Error', sProcess, sInput, 'Positive bin size required.');
        return;
    end

    % ===== LOAD INPUT FILES =====
    % Load data file
    DataMat = in_bst_data(sInput.FileName, 'Time', 'Events', 'Comment', 'Device', 'ChannelFlag', 'History');
    sampling_rate = round(abs(1. / (DataMat.Time(2) - DataMat.Time(1))));
    % Load channel file
    ChannelMat = in_bst_channel(sInput.ChannelFile);

    % ===== COMPUTE BINNING =====
    % Define bins
    nBins = floor(length(DataMat.Time) / (bin_size * sampling_rate));
    bins = linspace(DataMat.Time(1), DataMat.Time(end), nBins+1);
    single_file_binning = zeros(length(ChannelMat.Channel), nBins);
    % Process channel by channel
    for iChan = 1:length(ChannelMat.Channel)
        for iEvent = 1:size(DataMat.Events,2)

            % Bin ONLY THE FIRST NEURON'S SPIKES if there are multiple neurons!
            if panel_spikes('IsSpikeEvent', DataMat.Events(iEvent).label) ...
                    && panel_spikes('IsFirstNeuron', DataMat.Events(iEvent).label) ...
                    && strcmp(ChannelMat.Channel(iChan).Name, panel_spikes('GetChannelOfSpikeEvent', DataMat.Events(iEvent).label))

                outside_up = DataMat.Events(iEvent).times >= bins(end); % This snippet takes care of some spikes that occur outside of the window of Time due to precision incompatibility.
                DataMat.Events(iEvent).times(outside_up) = bins(end) - 0.001; % Make sure it is inside the bin. Add 1ms offset
                outside_down = DataMat.Events(iEvent).times <= bins(1);
                DataMat.Events(iEvent).times(outside_down) = bins(1) + 0.001; % Make sure it is inside the bin. Add 1ms offset

                [tmp, bin_it_belongs_to] = histc(DataMat.Events(iEvent).times, bins);

                unique_bin = unique(bin_it_belongs_to);
                occurences = [unique_bin; histc(bin_it_belongs_to, unique_bin)];

                single_file_binning(iChan,occurences(1,:)) = occurences(2,:)/bin_size; % The division by the bin_size gives the Firing Rate
                break
            end
        end

    end

    % Events have to be converted to the sampling rate of the binning
    convertedEvents = DataMat.Events;
    for iEvent = 1:length(DataMat.Events)
        [tmp, bin_it_belongs_to] = histc(DataMat.Events(iEvent).times, bins);
        bin_it_belongs_to(bin_it_belongs_to==0) = 1;
        convertedEvents(iEvent).times   = bins(bin_it_belongs_to);
    end
    Events = convertedEvents;

    % ===== SAVE RESULTS =====
    % Prepare output file structure
    FileMat = db_template('datamat');
    FileMat.F           = single_file_binning;
    FileMat.Time        = diff(bins(1:2))/2+bins(1:end-1);
    FileMat.Comment     = ['PSTH: ' DataMat.Comment];
    FileMat.DataType    = 'recordings';
    FileMat.ChannelFlag = DataMat.ChannelFlag;
    FileMat.Device      = DataMat.Device;
    FileMat.Events      = Events;
    FileMat.nAvg        = 1;
    FileMat.History     = DataMat.History;

    % Add history field
    FileMat = bst_history('add', FileMat, 'ptsh', ['PSTH per electrode: ' num2str(bin_size) ' ms']);
    FileMat = bst_history('add', FileMat, 'ptsh', ['Input file: ' sInput.FileName]);
    % Output filename
    FileName = bst_process('GetNewFilename', bst_fileparts(sInput.FileName), 'data_psth');
    OutputFiles = {FileName};
    % Save output file and add to database
    bst_save(FileName, FileMat, 'v6');
    db_add_data(sInput.iStudy, FileName, FileMat);
end

author: brainstorm-tools


text:
function sl_customization(cm)
% Customizes Simulink context menu with all blocks from a model file
% custom_context_menu_blocks.mdl file. These blocks are then accessible
% through context menu: Add custom block / List of blocks from library.
% The chosen block is added at the position of mouse context menu click.
%
% Please feel free to modify custom_context_menu_blocks.mdl to meet your
% expectations.
%
% OUTPUTS:  Customised Simulink context
%
% Author: Paul Proteus (e-mail: proteus.paul (at) yahoo (dot) com)
% Version: 1.0
% Changes tracker:  31.01.2011  - First version
%
% Note:  It might be neccessery to run sl_refresh_customizations to refresh
%        Simulink context menu or set the Matlab path correctly.
%
% Inspired by and parts of code taken from:
%   1. Alex - Add Simulink Commonly Used Blocks with context menu (#22067)
%   2. Giacomo Faggiani - Create From Blocks (#25736)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Register custom menu function.
cm.addCustomMenuFcn('Simulink:PreContextMenu', @getMyMenuItems);
end

%% Define the custom menu function.
function schemaFcns = getMyMenuItems(callbackInfo)
  schemaFcns = {@addSimulinkBlock};
end

%% Define the schema function for first menu item.
function schema = addSimulinkBlock(callbackInfo)
  schema = sl_container_schema;
  schema.label  = 'Add custom block';
  % Get mouse position
  mousePosition = get(0,'PointerLocation');
  % Get available block from the library
  customMDL = 'custom_context_menu_blocks';
  load_system(customMDL);
  availableBlocks = find_system(customMDL,'SearchDepth',1);
  for iBlock=2:length(availableBlocks),
     blockName = availableBlocks{iBlock};
     blockName = blockName(max(strfind(blockName,'/'))+1:end);
     childFunc(iBlock-1) = {{@addLibraryBlock,{customMDL,blockName,mousePosition}}};
  end
  schema.childrenFcns = childFunc;
end

function schema = addLibraryBlock(callbackInfo)
  schema = sl_action_schema;
  schema.label = callbackInfo.userdata{2};
  schema.userdata = callbackInfo.userdata;
  schema.callback = @slAddLibraryBlock;
end

% SUB-FUNCTIONS
function slAddLibraryBlock(inArgs)

% Define initial information
modelName     = gcs;
customMDL     = inArgs.Userdata{1};
blockName     = inArgs.Userdata{2};
mousePosition = inArgs.Userdata{3};

% Source and destination path
dstStr = [modelName '/' blockName];
srcStr = [customMDL '/' blockName];

% Get more information for position calculation
blockSizeRef = get_param(srcStr, 'Position');
locationBase = get_param(modelName,'Location');
scrollOffset = get_param(modelName,'ScrollbarOffset');
screenSize   = get(0,'ScreenSize');

% Calculate new position
p_X = mousePosition(1) - locationBase(1) + scrollOffset(1);
p_Y = screenSize(4) - mousePosition(2) - locationBase(2) + scrollOffset(2) ;
width = blockSizeRef(3) - blockSizeRef(1);
height = blockSizeRef(4) - blockSizeRef(2);
location = [p_X p_Y p_X+width p_Y+height];

% Add block
block = add_block(srcStr, dstStr, 'MakeNameUnique', 'on');
set(block, 'Position', location);

close_system(customMDL);

end


text:
function flag = isodd(n)
%ISODD returns true if the input is odd.

flag = mod(n,2);

% Copyright 2002 - 2009 The MathWorks, Inc.


text:
% Copyright (C) 1993-2013, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for MATLAB (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% We will show how to estimate the background of a scene, those pixels that
% are constant from frame to frame, and to estimate those parts of the scene
% that are changing.  We will demonstrate this on a live video feed.

% Create an object which represents a video camera connected to your computer
% If you have a laptop you should see the video recording light come on
%
% Type the q-key to stop this running

camera = VideoCamera('grey', 'double');

% take a frame from the camera as initial estimate of background
bg = camera.grab();

% rate of change threshold
sigma = 5/255.0;

% now we loop over all frames
while true
    % grab the image
    im = camera.grab();
    % compute the difference between the new frame and background estimate
    d = im-bg;
    % clip it to the range -sigma to sigma
    d = max(min(d, sigma), -sigma);
    % update the background, adjust background pixel values towards the new image
    bg = bg + d;

    % display live, estimated background and motion pixels
    idisp(im, 'nogui', 'figure', 1, 'title', 'live video');
    idisp(bg, 'nogui', 'figure', 2, 'title', 'background image');
    idisp(abs(im-bg), 'nogui', 'figure', 3, 'title', 'motion image');

    drawnow

    % test whether to quit, hit 'q'
    if get(gcf,'CurrentCharacter') == 'q'
        set(gcf,'CurrentCharacter', '-')
        break;
    end;
end
clear camera   % close the video source

author: petercorke

repo: machinevision-toolbox-matlab


text:
%VGG_SINGF_FROM_FF  Linearly combines two 3x3 matrices to a singular one.
%
%   a = vgg_singF_from_FF(F)  computes scalar(s) a such that given two 3x3 matrices F{1} and F{2},
%   it is det( a*F{1} + (1-a)*F{2} ) == 0.

function a = vgg_singF_from_FF(F)

% precompute determinants made from columns of F{1}, F{2}
for i1 = 1:2
  for i2 = 1:2
    for i3 = 1:2
      D(i1,i2,i3) = det([F{i1}(:,1) F{i2}(:,2) F{i3}(:,3)]);
    end
  end
end

% Solve The cubic equation for a
a = roots([-D(2,1,1)+D(1,2,2)+D(1,1,1)+D(2,2,1)+D(2,1,2)-D(1,2,1)-D(1,1,2)-D(2,2,2)
            D(1,1,2)-2*D(1,2,2)-2*D(2,1,2)+D(2,1,1)-2*D(2,2,1)+D(1,2,1)+3*D(2,2,2)
            D(2,2,1)+D(1,2,2)+D(2,1,2)-3*D(2,2,2)
            D(2,2,2)]);
a = a(abs(imag(a))<10*eps);

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function [ x, odata, opts ] = solver_SLOPE( A, b, lambda, x0, opts )
% SOLVER_SLOPE Sorted l1-regularized least squares problem,
% [ beta, odata, opts ] = solver_SLOPE( X, y, lambda, beta0, opts )
%    Solves the l1-regularized least squares problem, using the sorted/ordered l1 norm,
%        minimize (1/2)*norm( A * x - b )^2 + norm( lasso.*sort(abs(x),'descend'), 1 )
%    using the Auslender/Teboulle variant with restart. X must be a matrix
%    or a linear operator, y must be a vector, and lambda must be a real
%    positive vector in decreasing order.
%    The initial point beta0 and option structure opts are both optional.
%
% SLOPE stands for Sorted L-One Penalized Estimation
%
% Reference:
%   "Statistical Estimation and Testing via the Ordered l1 Norm"
%   by M. Bogdan, E. van den Berg, W. Su, and E. J. Candès, 2013
%   http://www-stat.stanford.edu/~candes/OrderedL1/
%
%   See also solver_L1RLS.m, solver_LASSO.m, prox_Sl1.m

error(nargchk(3,5,nargin));
if nargin < 4, x0 = []; end
if nargin < 5, opts = []; end
if ~isfield( opts, 'restart' ),
    opts.restart = 100;
end

[x,odata,opts] = tfocs( smooth_quad, { A, -b }, prox_Sl1( lambda ), x0, opts );

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [C, L, U] = SpectralClustering(W, k, Type)
%SPECTRALCLUSTERING Executes spectral clustering algorithm
%   Executes the spectral clustering algorithm defined by
%   Type on the adjacency matrix W and returns the k cluster
%   indicator vectors as columns in C.
%   If L and U are also called, the (normalized) Laplacian and
%   eigenvectors will also be returned.
%
%   'W' - Adjacency matrix, needs to be square
%   'k' - Number of clusters to look for
%   'Type' - Defines the type of spectral clustering algorithm
%            that should be used. Choices are:
%      1 - Unnormalized
%      2 - Normalized according to Shi and Malik (2000)
%      3 - Normalized according to Jordan and Weiss (2002)
%
%   References:
%   - Ulrike von Luxburg, "A Tutorial on Spectral Clustering",
%     Statistics and Computing 17 (4), 2007
%
%   Author: Ingo Buerk
%   Year  : 2011/2012
%   Bachelor Thesis

% calculate degree matrix
degs = sum(W, 2);
D    = sparse(1:size(W, 1), 1:size(W, 2), degs);

% compute unnormalized Laplacian
L = D - W;

% compute normalized Laplacian if needed
switch Type
    case 2
        % avoid dividing by zero
        degs(degs == 0) = eps;
        % calculate inverse of D
        D = spdiags(1./degs, 0, size(D, 1), size(D, 2));

        % calculate normalized Laplacian
        L = D * L;
    case 3
        % avoid dividing by zero
        degs(degs == 0) = eps;
        % calculate D^(-1/2)
        D = spdiags(1./(degs.^0.5), 0, size(D, 1), size(D, 2));

        % calculate normalized Laplacian
        L = D * L * D;
end

% compute the eigenvectors corresponding to the k smallest
% eigenvalues
diff   = eps;
[U, ~] = eigs(L, k, diff);

% in case of the Jordan-Weiss algorithm, we need to normalize
% the eigenvectors row-wise
if Type == 3
    U = bsxfun(@rdivide, U, sqrt(sum(U.^2, 2)));
end

% now use the k-means algorithm to cluster U row-wise
% C will be a n-by-1 matrix containing the cluster number for
% each data point
C = kmeans(U, k, 'start', 'sample', ...
                 'EmptyAction', 'singleton');

% now convert C to a n-by-k matrix containing the k indicator
% vectors as columns
C = sparse(1:size(D, 1), C, 1);

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
%% Dislocation Density Tensor
%
%% PLEASE HELP AND ADD CONTENT HERE
%

author: mtex-toolbox

repo: mtex

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


text:
function kernWriteParamsToFID(kern, FID)

% KERNWRITEPARAMSTOFID Write the kernel parameters to a stream.
% FORMAT
% DESC writes kernel parameters to a stream.
% ARG kern : the kernel that the parameters are in.
% ARG FID : the stream to which parameters are written.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2008
%
% SEEALSO : kernWriteToFID, modelWriteToFID

% KERN

if strcmp(kern.type, 'cmpnd') || strcmp(kern.type, 'tensor')
  componentKernWriteParamsToFID(kern, FID);
else
  writeIntToFID(FID, 'inputDim', kern.inputDimension);
  writeIntToFID(FID, 'numParams', kern.nParams);

  if strcmp(kern.type, 'poly') | strcmp(kern.type, 'polyard')
    writeIntToFID(FID, 'numParams', kern.degree);
  end

  fhandle = str2func([kern.type 'KernExtractParam']);
  params = fhandle(kern);
  doubleMatrixWriteToFID(params, FID);
  if isfield(kern, 'priors')
    writeIntToFID(FID, 'numPriors', length(kern.priors));
    for j=1:length(kern.priors)
      writeIntToFID(FID, 'priorIndex', kern.priors(j).index);
      priorWriteToFID(FID, prior);
    end
  else
    writeIntToFID(FID, 'numPriors', 0);
  end
end


text:
function C=covdelays(A, tdelay)
% C=covdelays(A, tdelays)
% Computes (unbiased) correlation matrices of MxN matrix A (rows are time series) for
% time delayes specified in the 1xT vector TDELAY

[m,n]=size(A);
t=length(tdelay);
C=zeros(m,m,t);
for indext = 1:t
    Ashift = (circshift(A',tdelay(indext)))';
    C(:,:,indext)=1/n*A*Ashift';

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function v1 = ctranspose(v1)
% transpose vector

v1.x = v1.x';
v1.y = v1.y';
v1.z = v1.z';


text:
function [xroot, froot] = brent (f, x1, x2, rtol)

% solve for a single real root of a nonlinear equation

% Brent's method

% input

%  f    = objective function coded as y = f(x)
%  x1   = lower bound of search interval
%  x2   = upper bound of search interval
%  rtol = algorithm convergence criterion

% output

%  xroot  = real root of f(x) = 0
%  froot  = function value at f(x) = 0

% Orbital Mechanics with Matlab

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

global iter;

% machine epsilon

eps = 2.23e-16;

e = 0;

a = x1;
b = x2;

fa = feval(f, a);

fb = feval(f, b);

fc = fb;

for iter = 1:1:50

    if (fb * fc > 0)
       c = a;
       fc = fa;
       d = b - a;
       e = d;
    end

    if (abs(fc) < abs(fb))
       a = b;
       b = c;
       c = a;
       fa = fb;
       fb = fc;
       fc = fa;
    end

    tol1 = 2 * eps * abs(b) + 0.5 * rtol;

    xm = 0.5 * (c - b);

    if (abs(xm) <= tol1 || fb == 0)
       break;
    end

    if (abs(e) >= tol1 && abs(fa) > abs(fb))
       s = fb / fa;

       if (a == c)
          p = 2 * xm * s;
          q = 1 - s;
       else
          q = fa / fc;
          r = fb / fc;
          p = s * (2 * xm * q * (q - r) - (b - a) * (r - 1));
          q = (q - 1) * (r - 1) * (s - 1);
       end

       if (p > 0)
          q = -q;
       end

       p = abs(p);

       min = abs(e * q);

       tmp = 3 * xm * q - abs(tol1 * q);

       if (min < tmp)
          min = tmp;
       end

       if (2 * p < min)
          e = d;
          d = p / q;
       else
          d = xm;
          e = d;
       end
    else
       d = xm;
       e = d;
    end

    a = b;
    fa = fb;

    if (abs(d) > tol1)
       b = b + d;
    else
       b = b + sign(xm) * tol1;
    end

    fb = feval(f, b);

end

xroot = b;
froot = fb;


text:
function tc = tc_sparklineWholeTc(tc);
%
% tc = tc_sparklineWholeTc(tc);
%
% An attempt to visualize a whole time course using
% Edward Tufte's sparkline method.
%
% ras, 01/2007.
if notDefined('tc'),    tc = get(gcf, 'UserData');      end

%% clear previous objects in figure
delete( findobj('Parent', tc.ui.plot) );

%% get time vector t in seconds
nFrames = length(tc.wholeTc);
t = [0:nFrames-1] .* tc.params.framePeriod;

%% plot time course in small axes
axes('Parent', tc.ui.plot, 'Units', 'norm', 'Position', [.1 .6 .65 .3]);
plot(t, tc.wholeTc, 'k');
hold on, axis off;
AX = axis;
w = AX(2) - AX(1);  % width
h = AX(4) - AX(3);  % height

%% add time point annotation underneath the time course
% mark only a few salient time points
nPoints = 8;
pts = [t(1) t(mod(t,100)==0) t(end)];
pts = pts( round( linspace(1, length(pts), nPoints) ) );
pts = unique(pts);

for p = pts
    f = find(t==p);
    plot(p, tc.wholeTc(f), 'r.');
    if p==pts(1)
        txt = sprintf('t = %i s', p);
    else
        txt = num2str(p);
    end
    text(p, -0.33*h, txt, 'FontName', 'Helvetica', ...
         'FontSize', 9, 'HorizontalAlignment', 'center');
end

%% add color bars indicating condition underneath the text annotation
cond = er_resample(tc.trials.onsetSecs, tc.trials.cond, t);
for i = 1:length(tc.trials.condColors)
    cmap(i,:) = colorLookup(tc.trials.condColors{i});
end
cmap(1,:) = [1 1 1];
condImg = ind2rgb(repmat(cond(:)'+1, 2, 1), cmap);
image(t, min(tc.wholeTc)-.2*h, condImg);
% axis([AX(1:2) AX(3)-h AX(4)])

% % add a scrollbar
% scrollbar(gca, tc.wholeTc);

return

% The # of samples we show per row depends on the current
% aspect ratio of the figure. (Although the user can
% resize it, we want it to look nice in the figure as given.)
pos = get(gcf, 'Position');
rFig = pos(3) / pos(4);  % xSize / ySize

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [I3, I2]=Cell_centri(I0,  coeff, latent, area_real, boundary, t5)
if(~exist('t5','var'))
    t5=.2;
end
area_right=area_real+2*t5*sqrt(latent(1))*coeff(:,1);
area_left=area_real-2*t5*sqrt(latent(1))*coeff(:,1);

[inner(:,2),inner(:,1)]=find(I0>=0);
centerlized1=inner-repmat(area_right', size(inner,1),1);
score1=centerlized1*coeff(:,1);

centerlized2=inner-repmat(area_left', size(inner,1),1);
score2=centerlized2*coeff(:,1);
I1=I0;
I1((score1.*score2)<0)=0;
I1((score1.*score2)>0)=1;
I2=I1;
figure;
imshow(I1)
hold on
h=imshow(I0);
alpha(h, 0.5)
I3=false(size(I0));

%% get rid of the part left
stats = regionprops(I0, 'Area','FilledArea','MajorAxisLength',...
    'MinorAxisLength', 'Orientation', 'Eccentricity','PixelIdxList');
AA=struct2cell(stats);
% area=zeros(1,size(AA,2));
% ratio=zeros(1,size(AA,2));
plist=cell(1,size(AA,2));
for i=1:size(AA,2)
%     area(i)=AA{1,i};
%     ratio(i)=AA{2,i}./AA{3,i};
    plist{i}=AA{7,i};
    if all(I1(plist{i})==0)
        I3(plist{i})=1;
    end

end
[A, B]=find(I3==1);
[IDX, D]=knnsearch(boundary, [B,A]);
index1=find(I3==1);
I3=false(size(I0));
I3(index1(D<1.5*sqrt(latent(2))))=1;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


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

[nx,p]=size(x); ny=size(y,1);
if nargin<3 | isempty(mode) mode='0'; end
if any(mode=='d') | (mode~='x' & nx==ny)
   nx=min(nx,ny);
   z=x(1:nx,:)-y(1:nx,:);
   if nargin<4
      d=sum(z.*conj(z),2);
   elseif min(size(w))==1
      wv=w(:).';
      d=sum(z.*wv(ones(size(z,1),1),:).*conj(z),2);
   else
      d=sum(z*w.*conj(z),2);
   end
else
   if p>1
      if nargin<4
         z=permute(x(:,:,ones(1,ny)),[1 3 2])-permute(y(:,:,ones(1,nx)),[3 1 2]);
         d=sum(z.*conj(z),3);
      else
         nxy=nx*ny;
         z=reshape(permute(x(:,:,ones(1,ny)),[1 3 2])-permute(y(:,:,ones(1,nx)),[3 1 2]),nxy,p);
         if min(size(w))==1
            wv=w(:).';
            d=reshape(sum(z.*wv(ones(nxy,1),:).*conj(z),2),nx,ny);
         else
            d=reshape(sum(z*w.*conj(z),2),nx,ny);
         end
      end
   else
      z=x(:,ones(1,ny))-y(:,ones(1,nx)).';
      if nargin<4
         d=z.*conj(z);
      else
         d=w*z.*conj(z);
      end
   end
end
if any(mode=='s')
   d=sqrt(d);
end

author: bastamon


text:
function [ x, seed ] = disk01_sample ( n, seed )

%*****************************************************************************80
%
%% DISK01_SAMPLE uniformly samples the unit disk.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real X(2,N), the points.
%
  x = randn ( 2, n );
  norm = ones ( 1, 2 ) * ( x.^2 );
  norm = sqrt ( norm );
  for i = 1 : 2
    x(i,1:n) = x(i,1:n) ./ norm(1:n);
  end

  for j = 1 : n
    r = rand ( 1, 1 );
    x(1:2,j) = sqrt ( r  ) * x(1:2,j);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function obj = t3_pTop_func_new(in1,in2,in3)
coefs_tq1_1 = in2(1);
coefs_tq1_2 = in2(4);
coefs_tq1_3 = in2(7);
coefs_tq1_4 = in2(10);
coefs_tq1_5 = in2(13);
coefs_tq1_6 = in2(16);
coefs_tq1_7 = in2(19);
coefs_tq1_8 = in2(22);
coefs_tq1_9 = in2(25);
coefs_tq2_1 = in2(2);
coefs_tq2_2 = in2(5);
coefs_tq2_3 = in2(8);
coefs_tq2_4 = in2(11);
coefs_tq2_5 = in2(14);
coefs_tq2_6 = in2(17);
coefs_tq2_7 = in2(20);
coefs_tq2_8 = in2(23);
coefs_tq2_9 = in2(26);
coefs_tq3_1 = in2(3);
coefs_tq3_2 = in2(6);
coefs_tq3_3 = in2(9);
coefs_tq3_4 = in2(12);
coefs_tq3_5 = in2(15);
coefs_tq3_6 = in2(18);
coefs_tq3_7 = in2(21);
coefs_tq3_8 = in2(24);
coefs_tq3_9 = in2(27);
coefs_tq1_10 = in2(28);
coefs_tq1_11 = in2(31);
coefs_tq2_10 = in2(29);
coefs_tq2_11 = in2(32);
coefs_tq3_10 = in2(30);
coefs_tq3_11 = in2(33);
pinvG3_1 = in1(3);
pinvG3_2 = in1(6);
pinvG3_3 = in1(9);
q0 = in3(1,:);
q1 = in3(2,:);
q2 = in3(3,:);
q3 = in3(4,:);
obj = q0.^2.*(coefs_tq1_1.*pinvG3_1+coefs_tq2_1.*pinvG3_2+coefs_tq3_1.*pinvG3_3)+q1.^2.*(coefs_tq1_5.*pinvG3_1+coefs_tq2_5.*pinvG3_2+coefs_tq3_5.*pinvG3_3)+q2.^2.*(coefs_tq1_8.*pinvG3_1+coefs_tq2_8.*pinvG3_2+coefs_tq3_8.*pinvG3_3)+q3.^2.*(coefs_tq1_10.*pinvG3_1+coefs_tq2_10.*pinvG3_2+coefs_tq3_10.*pinvG3_3)+coefs_tq1_11.*pinvG3_1+coefs_tq2_11.*pinvG3_2+coefs_tq3_11.*pinvG3_3+q0.*q1.*(coefs_tq1_2.*pinvG3_1+coefs_tq2_2.*pinvG3_2+coefs_tq3_2.*pinvG3_3)+q0.*q2.*(coefs_tq1_3.*pinvG3_1+coefs_tq2_3.*pinvG3_2+coefs_tq3_3.*pinvG3_3)+q0.*q3.*(coefs_tq1_4.*pinvG3_1+coefs_tq2_4.*pinvG3_2+coefs_tq3_4.*pinvG3_3)+q1.*q2.*(coefs_tq1_6.*pinvG3_1+coefs_tq2_6.*pinvG3_2+coefs_tq3_6.*pinvG3_3)+q1.*q3.*(coefs_tq1_7.*pinvG3_1+coefs_tq2_7.*pinvG3_2+coefs_tq3_7.*pinvG3_3)+q2.*q3.*(coefs_tq1_9.*pinvG3_1+coefs_tq2_9.*pinvG3_2+coefs_tq3_9.*pinvG3_3);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function f = ddivE(Ekt, varargin)
% function f = ddivE(Ekt, varargin)
% Vxt = varargin{1};  %data
% Wxk = varargin{2};  %estimated H from NMF update
% alpha = varargin{3}; %rate

Vxt = varargin{1};  %data
Wxk = varargin{2};  %estimated values of W
alpha = varargin{3};

Hkt=exp(alpha*Ekt);
f = ddivergence(Vxt,Wxk*Hkt);


text:
data_masterdir = GetCurrentDataDir();
range_fish = [2,3,4,5,6];
%% a) Histogram of Pairwise Distances for One Fish

%
clear T N
figure;
hold on;
for i = 1:length(range_fish),
    i_fish = range_fish(i);
    [f_cells,M,N(i),T(i),hfig] = getFullFishData(hfig,i_fish,VAR);

    num_select = 1000;%num_cells;
    idx_select = randperm(N(i), num_select);

    f_cells_select = f_cells(idx_select,:);

    %% Correlation
    tic
    Dpairs = 1-pdist(f_cells_select,'correlation');
    % Dpairs is in correlation, not corr dist
    toc

    %% Plot Histogram
    binwidth = .005;
    edges = -1:binwidth:1;

    [counts,edges] = histcounts(Dpairs,edges);
    allBins = edges(1:end-1) + binwidth/2;

    allCorrVals = counts;
    leg{i} = ['Fish ',num2str(i_fish)];

    plot(allBins,allCorrVals/trapz(allBins,allCorrVals));
    xlabel('Pairwise Correlations');
    ylabel('Normalized Frequency');
    xlim([-0.2,0.2]);

end

% Compare to Random Data
T_control = T(end);
f_cells_control = normrnd(0,1,num_select,T_control);
Dpairs_control = 1-pdist(f_cells_control,'correlation');

[ctrlCounts,edgesCtrl] = histcounts(Dpairs_control,edges);
control_allBins = edgesCtrl(1:end-1) + binwidth/2;
control_allCorrVals = ctrlCounts;

plot(control_allBins,control_allCorrVals/trapz(control_allBins,control_allCorrVals)...
    ,'k');

leg{i+1} = 'Control';
legend(leg)

%% b) Percentile Threshold Plot

% Percentile Thresh
threshPerc = .99;
figure;
xlabel('Number of Frames');
ylabel('99th Percentile Correlation');

hold on;
for num_fish = 1:length(range_fish),

    i_fish = range_fish(num_fish);
    [f_cells,M,N(num_fish),T(num_fish),hfig] = getFullFishData(hfig,i_fish,VAR);

    num_select = 1000;%num_cells;
    idx_select = randperm(N(num_fish), num_select);

    f_cells_select = f_cells(idx_select,:);

    T_trunc = T(num_fish)/10:T(num_fish)/10:T(num_fish); %hardcode for now
    clear threshCorr;
    for i = 1:length(T_trunc)
        f_trunc = f_cells_select(:,1:T_trunc(i));

        % Correlation Distance
        tic
        Dpairs = pdist(f_trunc,'correlation');
        % 138 sec for 86k cells
        % 1.8 sec for 10k cells
        toc
        %figure;
        %h = histogram(1-Dpairs);%Plot corr, not corr dist

        sortCorr = sort(1-Dpairs);
        threshIDX = round(length(sortCorr)*threshPerc);
        threshCorr(i) = sortCorr(threshIDX);
    end
    plot(T_trunc,threshCorr,'o-')
    leg2{num_fish} = ['Fish ',num2str(i_fish)];
end

% Compare to Random Data
T_control = T(end);
f_cells_control = normrnd(0,1,num_select,T_control);
Dpairs_control = 1-pdist(f_cells_control,'correlation');

T_trunc = T_control/10:T_control/10:T_control; %hardcode for now
clear threshCorr;
clear threshCorr;
for i = 1:length(T_trunc)
    f_trunc = f_cells_control(:,1:T_trunc(i));

    % Correlation Distance
    tic
    Dpairs = pdist(f_trunc,'correlation');
    % 138 sec for 86k cells
    % 1.8 sec for 10k cells
    toc
    %figure;
    %h = histogram(1-Dpairs);%Plot corr, not corr dist
    sortCorr = sort(1-Dpairs);
    threshIDX = round(length(sortCorr)*threshPerc);
    threshCorr(i) = sortCorr(threshIDX);
end
plot(T_trunc,threshCorr,'ko-')
leg2{num_fish+1} = 'control';
ylim([0 .5]);
legend(leg2);
grid on;

author: xiuyechen

repo: FishExplorer


text:
function printx(x,rtk)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
opt=rtk.opt; nf=rtk.NF;

fprintf('pos   = ');fprintf('%14.5f %14.5f %14.5f',x(1),x(2),x(3));fprintf('\n');

if rtk.NP==9
    fprintf('vel   = ');fprintf('%9.5f%9.5f%9.5f',x(4),x(5),x(6));fprintf('\n');
    fprintf('acc   = ');fprintf('%9.5f%9.5f%9.5f',x(7),x(8),x(9));fprintf('\n');
end

if rtk.NI>0
    fprintf('iono  = ');
    for i=1:glc.MAXSAT
        if abs(x(rtk.ii+i))>=1e-6,fprintf('%9.5f',x(rtk.ii+i));end
    end
    fprintf('\n');
end

if rtk.NT>0
    fprintf('trop  = ');
    if opt.tropopt==glc.TROPOPT_EST
        fprintf('%9.5f',x(rtk.itr+1));
        fprintf('%9.5f',x(rtk.itb+1));
    else
        fprintf('%9.5f',x(rtk.itr+1));
        fprintf('%9.5f',x(rtk.itr+2));
        fprintf('%9.5f',x(rtk.itr+3));
        fprintf('%9.5f',x(rtk.itb+1));
        fprintf('%9.5f',x(rtk.itb+2));
        fprintf('%9.5f',x(rtk.itb+3));
    end
    fprintf('\n');
end

if rtk.NL>0
    fprintf('glo   = ');
    fprintf('%9.5f',x(rtk.il+1)); fprintf('%9.5f',x(rtk.il+2));
    fprintf('\n');
end

if rtk.NB>0
    fprintf('bias1 = ');
    for i=1:glc.MAXSAT
        if abs(x(rtk.ib+i))>=1e-6,fprintf('%10.5f',x(rtk.ib+i));end
    end
    fprintf('\n');
    if nf>=2
        fprintf('bias2 = ');
        for i=1:glc.MAXSAT
            if abs(x(rtk.ib+glc.MAXSAT+i))>=1e-6
                fprintf('%10.5f',x(rtk.ib+glc.MAXSAT+i));
            end
        end
        fprintf('\n');
    end
    if nf>=3
        fprintf('bias3 = ');
        for i=1:glc.MAXSAT
            if abs(x(rtk.ib+glc.MAXSAT*2+i))>=1e-6
                fprintf('%10.5f',x(rtk.ib+glc.MAXSAT*2+i));
            end
        end
        fprintf('\n');
    end
end

return

repo: GINav


text:
function disp(f)
% DISP - Filterobject disp overloaded operator

% VERSION: 1.0 of filter objects
% AUTHOR: Celso Reyes
% LASTUPDATE: 1/30/2007

if numel(f) > 1;
    disp(' ');
    DispStr = sprintf('%d',size(f,1));
    for n = 2 : numel(size(f))
        DispStr = sprintf('%sx%d', DispStr, size(f,n));
    end
    disp(sprintf('%s %s object with fields:', DispStr, class(f)));
    disp('    type');
    disp('    cutoff');
    disp('    poles');
else
    switch f.type
        case 'B'
            typename = 'Bandpass';
        case 'L'
            typename = 'Low-pass';
        case 'H'
            typename = 'High-pass';
    end;
    disp(['      type: ' f.type ' (' typename ')']);
    disp(['    cutoff: ' num2str(f.cutoff) ' Hz']);
    disp(['     poles: ' num2str(f.poles)]);
end;

author: geoscience-community-codes

repo: GISMO


text:
function pool_forward()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.activations{curr_layer_idx} = bsxfun(@plus, bsxfun(@times, mem.layer_inputs{curr_layer_idx}, config.weights{curr_layer_idx}'), ...
                                                      config.weights{config.layer_num+curr_layer_idx}');
    config.misc.current_layer = curr_layer_idx + 1;
end

author: jimmy-ren


text:
function failed_tutorial_spike

% MEM 8gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_spike ft_spike_select ft_spike_waveform ft_spike_maketrials ft_spike_isi ft_spike_plot_isireturn ft_spike_psth ft_spikedensity ft_spike_plot_raster ft_spike_rate ft_spike_select ft_spike_xcorr ft_spike_jpsth ft_spike_plot_jpsth

filenex = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/spikefield/p029_sort_final_01.nex');
spike   = ft_read_spike(filenex);

cfg              = [];
cfg.spikechannel = {'sig002a_wf', 'sig003a_wf'}; % select only the two single units
spike = ft_spike_select(cfg, spike);

cfg             = [];
cfg.fsample     = 40000;
cfg.interpolate = 1; % keep the density of samples as is
[wave, spikeCleaned] = ft_spike_waveform(cfg,spike);

for k = [1 2]
  figure,
  sl = squeeze(wave.dof(k,:,:))>1000; % only keep samples with enough spikes
  plot(wave.time(sl), squeeze(wave.avg(k,:,sl)),'k') % factor 10^6 to get microseconds
  hold on

  % plot the standard deviation
  plot(wave.time(sl), squeeze(wave.avg(k,:,sl))+sqrt(squeeze(wave.var(k,:,sl))),'k--')
  plot(wave.time(sl), squeeze(wave.avg(k,:,sl))-sqrt(squeeze(wave.var(k,:,sl))),'k--')

  axis tight
  set(gca,'Box', 'off')
  xlabel('time')
  ylabel('normalized voltage')
end

event = ft_read_event(filenex);

cfg          = [];
cfg.dataset  = filenex;
cfg.trialfun = 'trialfun_stimon';
cfg = ft_definetrial(cfg);

cfg.timestampspersecond =  spike.hdr.FileHeader.Frequency; % 40000
spikeTrials = ft_spike_maketrials(cfg,spike);

spikeTrials.trialtime(1:5,:)
spikeTrials.cfg.trl(1:5,:)

cfg                     = [];
hdr                     = ft_read_header(filenex);
cfg.trl                 = [0 hdr.nSamples*hdr.TimeStampPerSample 0];
cfg.timestampspersecond =  spike.hdr.FileHeader.Frequency; % 40000
spike_notrials   = ft_spike_maketrials(cfg,spike);

dat = ft_checkdata(spikeTrials,'datatype', 'raw', 'fsample', 1000);

dat.trial{1}(:,4000:4004)

spike_converted = ft_checkdata(dat,'datatype', 'spike');

cfg       = [];
cfg.bins  = [0:0.0005:0.1]; % use bins of 0.5 milliseconds
cfg.param = 'coeffvar'; % compute the coefficient of variation (sd/mn of isis)
isih = ft_spike_isi(cfg,spikeTrials);

for k = [1 2] % only do for the single units
  cfg              = [];
  cfg.spikechannel = isih.label{k};
  cfg.interpolate  = 5; % interpolate at 5 times the original density
  cfg.window       = 'gausswin'; % use a gaussian window to smooth
  cfg.winlen       = 0.004; % the window by which we smooth has size 4 by 4 ms
  cfg.colormap     = jet(300); % colormap
  cfg.scatter      = 'no'; % do not plot the individual isis per spike as scatters
  figure, ft_spike_plot_isireturn(cfg,isih)
end

% read in the .t file
filet = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/spike/tt6_7.t');

cfg         = [];
cfg.dataset = filet;
spike2 = ft_read_spike(cfg.dataset);

% convert timestamps to seconds
cfg                     = [];
cfg.trl                 = [0 max(spike2.timestamp{1})+1 0];
cfg.timestampspersecond = 10^6;
spike2Trial = ft_spike_maketrials(cfg,spike2);

% run the isi histogram
cfg      = [];
cfg.bins = [0:0.001:0.2];
isih = ft_spike_isi(cfg,spike2Trial);

% plot the isi histogram with the Poincare return map
cfg             = [];
cfg.interpolate = 5;
cfg.window      = 'gausswin';
cfg.winlen      = 0.005;
cfg.scatter     = 'no';
cfg.colormap    = jet(300);
figure, ft_spike_plot_isireturn(cfg,isih)

cfg             = [];
cfg.binsize     =  0.1; % if cfgPsth.binsize = 'scott' or 'sqrt', we estimate the optimal bin size from the data itself
cfg.outputunit  = 'rate'; % give as an output the firing rate
cfg.latency     = [-1 3]; % between -1 and 3 sec.
cfg.vartriallen = 'yes'; % variable trial lengths are accepted
cfg.keeptrials  = 'yes'; % keep the psth per trial in the output
psth = ft_spike_psth(cfg,spikeTrials);

cfg         = [];
cfg.binsize =  [0.05];
cfg.latency = [-1 3];
psth = ft_spike_psth(cfg,spikeTrials);

cfg              = [];
cfg.topplotfunc  = 'line'; % plot as a line
cfg.spikechannel = spikeTrials.label([1 2]);
cfg.latency      = [-1 3];
cfg.errorbars    = 'std'; % plot with the standard deviation
cfg.interactive  = 'no'; % toggle off interactive mode
figure, ft_spike_plot_raster(cfg,spikeTrials, psth)

cfg         = [];
cfg.latency = [-1 3];
cfg.timwin  = [-0.025 0.025];
cfg.fsample = 1000; % sample at 1000 hz
sdf = ft_spikedensity(cfg,spikeTrials);

cfg              = [];
cfg.topplotfunc  = 'line'; % plot as a line plot
cfg.spikechannel = spikeTrials.label([1 2]); % can also select one unit here
cfg.latency      = [-1 3];
cfg.errorbars    = 'std'; % plot with standard deviation
cfg.interactive  = 'no'; % toggle off interactive mode
figure, ft_spike_plot_raster(cfg,spikeTrials, sdf)

cfg         = [];
cfg.latency = [-1 3];
cfg.timwin  = [-0.1 0.1];
cfg.fsample = 1000; % sample at 1000 hz
[sdf, sdfdata] = ft_spikedensity(cfg,spikeTrials);

sdfdata.trial{1}(:,3000:3005)

cfg            = [];
cfg.latency    = [0.3 max(spikeTrials.trialtime(:))]; % sustained response period
cfg.keeptrials = 'yes';
rate = ft_spike_rate(cfg,spikeTrials);

[R,P] = corrcoef(rate.trial);

% read in the data, select the channels and define the trials
spike = ft_read_spike(filenex);

cfg              = [];
cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'}; % select only the MUA
spike = ft_spike_select(cfg, spike);

cfg          = [];
cfg.dataset  = filenex;
cfg.trialfun = 'trialfun_stimon';
cfg = ft_definetrial(cfg);
cfg.timestampspersecond =  spike.hdr.FileHeader.Frequency; % 40000
spikeTrials = ft_spike_maketrials(cfg,spike);

cfg             = [];
cfg.maxlag      = 0.2; % maximum 200 ms
cfg.binsize     = 0.001; % bins of 1 ms
cfg.outputunit  = 'proportion'; % make unit area
cfg.latency     = [-2.5 0];
cfg.vartriallen = 'no'; % do not allow variable trial lengths
cfg.method      = 'xcorr'; % compute the normal cross-correlogram
Xc = ft_spike_xcorr(cfg,spikeTrials);

% compute the shuffled correlogram
cfg.method      = 'shiftpredictor'; % compute the shift predictor
Xshuff = ft_spike_xcorr(cfg,spikeTrials);

iCmb = 3;
jCmb = 4;
figure
xcSmoothed = conv(squeeze(Xc.xcorr(iCmb,jCmb,:)),ones(1,5)./5,'same'); % do some smoothing
hd = plot(Xc.time(3:end-2),xcSmoothed(3:end-2),'k'); % leave out borders (because of smoothing)
hold on
xcSmoothed = conv(squeeze(Xshuff.shiftpredictor(iCmb,jCmb,:)),ones(1,5)./5,'same');
plot(Xc.time(3:end-2),xcSmoothed(3:end-2),'r')
hold on
xlabel('delay')
ylabel('proportion of coincidences')
title([Xc.label{iCmb} Xc.label{jCmb}])
axis tight

% compute the spike densities
cfg         = [];
cfg.latency = [-1 3];
cfg.timwin  = [-0.025 0.025];
cfg.fsample = 200;
[sdf] = ft_spikedensity(cfg,spikeTrials);

% compute the joint psth
cfg               = [];
cfg.normalization = 'no';
cfg.channelcmb    = spikeTrials.label(3:4);
cfg.method        = 'jpsth';
jpsth = ft_spike_jpsth(cfg,sdf);

cfg.method = 'shiftpredictor';
jpsthShuff = ft_spike_jpsth(cfg,sdf);

% subtract the predictor
jpsthSubtr = jpsth;
jpsthSubtr.jpsth = jpsth.jpsth-jpsthShuff.shiftpredictor;

cfg        = [];
figure
ft_spike_plot_jpsth(cfg,jpsth)
figure
ft_spike_plot_jpsth(cfg,jpsthShuff)
figure
ft_spike_plot_jpsth(cfg,jpsthSubtr)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [A, names] = mk_adj_mat(connections, names, topological)
% MK_ADJ_MAT Make a directed adjacency matrix from a list of connections between named nodes.
%
% A = mk_adj_mat(connections, name)
% This is best explaine by an example:
%   names = {'WetGrass', 'Sprinkler', 'Cloudy', 'Rain'};
%   connections = {'Cloudy', 'Sprinkler'; 'Cloudy', 'Rain'; 'Sprinkler', 'WetGrass'; 'Rain', 'WetGrass'};
% adds the arcs C -> S, C -> R, S -> W, R -> W. Node 1 is W, 2 is S, 3 is C, 4 is R.
%
% [A, names] = mk_adj_mat(connections, name, 1)
% The last argument of 1 indicates that we should topologically sort the nodes (parents before children).
% In the example, the numbering becomes: node 1 is C, 2 is R, 3 is S, 4 is W
% and the return value of names gets permuted to {'Cloudy', 'Rain', 'Sprinkler', 'WetGrass'}.
% Note that topological sorting the graph is only possible if it has no directed cycles.

if nargin < 3, topological = 0; end

n=length(names);
A=zeros(n);
[nr nc] = size(connections);
for r=1:nr
  from = strmatch(connections{r,1}, names, 'exact');
  assert(~isempty(from));
  to = strmatch(connections{r,2}, names, 'exact');
  assert(~isempty(to));
  %fprintf(1, 'from %s %d to %s %d\n', connections{r,1}, from, connections{r,2}, to);
  A(from,to) = 1;
end

if topological
  order = topological_sort(A);
  A = A(order, order);
  names = names(order);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
 function run_mfile_local(arg, varargin)
%function run_mfile_local(arg)
%|
%| run an mfile in a local environment so that workspace variables
%| are untouched.  useful for tests.
%| arg can be a cell with many mfiles to run
%| options
%|	'draw'	1|0	1 to draw after each test (default: 0)
%|	'pause'	1|0	1 to pause after each test (default: 0)
%|	'abort'	1|0	1 to abort if any test fails (default: 0)
%|
%| Copyright 2005-6-21, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

opt.draw = false;
opt.pause = false;
opt.abort = false;
opt = vararg_pair(opt, varargin);

% track which tests open a figure even though they should not if im disabled
check_fig = ~figure_opened && ~im;
bad_fig = {};

test_bad = {};
test_good = {};
if iscell(arg)

	tmp = cpu('etic');

	for ii=1:length(arg)
		printf('\n\nTesting: %s\n\n', arg{ii})
		try
			time1 = cpu('etic');
			run_mfile_local(arg{ii})
			time1 = cpu(time1);
			test_good{end+1} = sprintf('%s (%.2f)', arg{ii}, time1);
		catch
			time1 = cpu(time1);
			test_bad{end+1} = sprintf('%s (%.2f)', arg{ii}, time1);
			if opt.abort
				fail('test "%s" failed', arg{ii})
			end
		end
		if opt.draw, drawnow, end
		if opt.pause, printm 'pausing', pause, end

		drawnow;
		if check_fig && figure_opened
			bad_fig{end+1} = arg{ii};
			close all
		end
	end

	if length(test_good)
		printf('\n----\n')
		printm('in "%s" the following test(s) passed:\n', caller_name)
		disp(char(test_good))
	end

	if length(bad_fig)
		warn('in %s: the following test(s) had figure issues:', caller_name)
		disp(char(bad_fig))
	end

	cpu(tmp, [caller_name ' time'])

	if length(test_bad)
		warn('in "%s" the following test(s) failed:', caller_name)
		disp(char(test_bad))
		error 'the above tests failed'
	else
		printm('in "%s" all %d tests passed!', caller_name, length(arg))
	end

else
	eval(arg)
end

% try to determine if any figure window is opened
function out = figure_opened
if isfreemat
	out = true; % freemat will not tell, so assume yes
else
	out = ~isempty(get(0, 'children')); % matlab knows
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function plotem()

figure(1); clf; hold on;
pr = load('pb/bgtg_0.01_0.02_64/pr.txt');
plot(pr(:,2),pr(:,3),'co-');
pr = load('pb/bgtg/pr.txt');
plot(pr(:,2),pr(:,3),'bx-');
pr = load('pb/cgtg/pr.txt');
plot(pr(:,2),pr(:,3),'rx-');
prplot('');
legend('bgtg zone=2','bgtg zone=10','cgtg');
return

figure(1); clf; hold on;
prplot('all');
pr = load('pb/bgtg_0.01_0.02_64/pr.txt');
plot(pr(:,2),pr(:,3),'c-');
pr = load('pb/bgtg_hys_0.1/pr.txt');
plot(pr(:,2),pr(:,3),'rx-');
pr = load('pb/bgtg_hys_0.3/pr.txt');
plot(pr(:,2),pr(:,3),'ro-');
pr = load('pb/bgtg_hys/pr.txt');
plot(pr(:,2),pr(:,3),'mo-');
pr = load('pb/bgtg_hys_0.5/pr.txt');
plot(pr(:,2),pr(:,3),'rv-');
pr = load('pb/bgtg_hys_0.7/pr.txt');
plot(pr(:,2),pr(:,3),'bx-');
pr = load('pb/bgtg_hys_0.9/pr.txt');
plot(pr(:,2),pr(:,3),'bo-');
legend('base','0.1','0.3','0.4','0.5','0.7','0.9');
return

figure(1); clf; hold on;
pr = load('pb/random/pr.txt');
plot(pr(:,2),pr(:,3),'cx-');
pr = load('pb/gm_2/pr.txt');
plot(pr(:,2),pr(:,3),'mx-');
pr = load('pb/canny_2/pr.txt');
plot(pr(:,2),pr(:,3),'mo-');
pr = load('pb/gm_2_4/pr.txt');
plot(pr(:,2),pr(:,3),'mv-');
pr = load('pb/2mm_2/pr.txt');
plot(pr(:,2),pr(:,3),'gx-');
pr = load('pb/bgtg_0.01_0.02_64/pr.txt');
plot(pr(:,2),pr(:,3),'bx-');
pr = load('pb/bgtg_hys/pr.txt');
plot(pr(:,2),pr(:,3),'bo-');
prplot('all');
legend('random','gm(2)','canny(2)','gm(2,4)','2mm(2)','bgtg','bgtg(hys)');
return

figure(2); clf; hold on;
pr = load('pb/gm_1/pr.txt');
plot(pr(:,2),pr(:,3),'bx-');
pr = load('pb/gm_2/pr.txt');
plot(pr(:,2),pr(:,3),'bv-');
pr = load('pb/gm_4/pr.txt');
plot(pr(:,2),pr(:,3),'b^-');
pr = load('pb/gm_8/pr.txt');
plot(pr(:,2),pr(:,3),'bo-');
pr = load('pb/gm_16/pr.txt');
plot(pr(:,2),pr(:,3),'bs-');
prplot('GM, Single Scale');
legend('gm-1','gm-2','gm-4','gm-8','gm-16');
return

figure(3); clf; hold on;
pr = load('pb/2mm_1/pr.txt');
plot(pr(:,2),pr(:,3),'rx-');
pr = load('pb/2mm_2/pr.txt');
plot(pr(:,2),pr(:,3),'rv-');
pr = load('pb/2mm_4/pr.txt');
plot(pr(:,2),pr(:,3),'r^-');
pr = load('pb/2mm_8/pr.txt');
plot(pr(:,2),pr(:,3),'ro-');
pr = load('pb/2mm_16/pr.txt');
plot(pr(:,2),pr(:,3),'rs-');
prplot('2MM, Single Scale');
legend('2mm-1','2mm-2','2mm-4','2mm-8','2mm-16');

figure(4); clf; hold on;
pr = load('pb/gm_4/pr.txt');
plot(pr(:,2),pr(:,3),'rx-');
pr = load('pb/gm_1_2/pr.txt');
plot(pr(:,2),pr(:,3),'bx-');
pr = load('pb/gm_2_4/pr.txt');
plot(pr(:,2),pr(:,3),'bv-');
pr = load('pb/gm_4_8/pr.txt');
plot(pr(:,2),pr(:,3),'b^-');
prplot('GM, Multi-Scale');
legend('gm-4','gm-1-2','gm-2-4','gm-4-8');

figure(5); clf; hold on;
pr = load('pb/2mm_2/pr.txt');
plot(pr(:,2),pr(:,3),'rx-');
pr = load('pb/2mm_1_2/pr.txt');
plot(pr(:,2),pr(:,3),'bx-');
pr = load('pb/2mm_2_4/pr.txt');
plot(pr(:,2),pr(:,3),'bv-');
prplot('2MM, Multi-Scale');
legend('2mm-2','2mm-1-2','2mm-2-4');

function prplot(ti)
axis([0 1 0 1]);
axis square;
box on;
xlabel('Recall');
ylabel('Precision');
title(ti);

author: Cloud-CV

repo: object-proposals


text:
function string = ASAversion2numstr(version)
%ASAVERSION2NUMSTR Version identifier to number as string type
%   NUMSTRING = ASAVERSION2NUMSTR(VERSION) converts a numeric ARMASA
%   version identifier into a character string NUMSTRING containing a
%   numeric conversion of the version identifier.
%
%   ARMASA function version identifiers are implemented as date vectors
%   according to the Matlab format, see DATEVEC. Each date vector can be
%   associated with an equivalent serial date number, see DATENUM.
%   ASAVERSION2NUMSTR accepts both numeric formats and converts them into
%   the alternative numeric format which is subsequently converted into a
%   character string.
%
%   Example:
%
%   This function can be used to quickly generate a version identifier
%   for an ARMASA main function that has just been created or modified.
%   Typing,
%
%     ASAversion2numstr(now)
%
%   returns the current date and time in a date vector string that can
%   directly be pasted at the appropriate location in the ARMASA function
%   m-file.
%
%   See also: ASACONTROL, ASAVERSIONCHK, DATENUM, DATESTR, DATEVEC, NOW.

if ~isnumeric(version) |  ~(isequal(length(version),6) | isequal(length(version),1))
   error(ASAerr(10,{'version';mfilename}))
end

if isequal(length(version),6)
   v = version;
   string = sprintf('%0.5f',datenum(v(1),v(2),v(3),v(4),v(5),v(6)));
else
   string = sprintf('%0.0f %0.0f %0.0f %0.0f %0.0f %0.0f',datevec(version));
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% [2000 12 30 20 0 0]    W. Wunderink           wwunderink01@freeler.nl

author: Sable

repo: mcbench-benchmarks


text:
%--------------------------------------------------------------------------
%	SurfBox-MATLAB (c)
%--------------------------------------------------------------------------
%
%	Yue M. Lu and Minh N. Do
%
%--------------------------------------------------------------------------
%
%	SurfEP.m
%
%	First created: 12-12-08
%	Last modified: 12-12-08
%
%--------------------------------------------------------------------------

function Ep = SurfEP(sz, Lev_array, Pyr_mode, bo, nRepeat)

%   Estimating the subband scaling coefficients of the surfacelet transform
%   using a Monte Carlo method.
%
%   Input:
%
%   sz: the size of a N-dimensional array (N >= 2)
%
%   Lev_array: an L by 1 cell array, with each cell being an N by N matrix
%   for NDFB decomposition, or a scaler (either 1 or 2) for dual-tree
%   wavelet decomposition.
%
%   Pyr_mode: type of multiscale pyramid, corresponding to different levels
%   of redundancy.
%   1:       ~ 6.4
%   1.5:     ~ 4.0
%   2:       ~ 3.4
%
%   bo: the order of the checkerboard filters. Default: bo = 12
%
%   nRepeat: the number of trials in the Monte Carlo method.
%
%   Output:
%
%   Ep: the estimated subband scaling coefficients

h = waitbar(0, ['Estimating the subband scaling coefficients.' ...
        'This can be a slow process, but it only needs to be done once for a given configuration.']);

for n = 1 : nRepeat
    X = randn(sz);
    [Y, Recinfo] = Surfdec(X, Pyr_mode, Lev_array, 'ritf', 'bo', bo);
    clear X; % to save some memory

    if n == 1
        Ep = Y;

        for j = 1 : length(Y) - 1
            for dd = 1 : length(Y{j})
                for sd = 1 : length(Y{j}{dd})
                    Ep{j}{dd}{sd} = Ep{j}{dd}{sd} .^ 2;
                end
            end
        end

    else
        for j = 1 : length(Y) - 1
            for dd = 1 : length(Y{j})
                for sd = 1 : length(Y{j}{dd})
                    Ep{j}{dd}{sd} = Ep{j}{dd}{sd} + Y{j}{dd}{sd} .^ 2;
                end
            end
        end
    end

    clear Y;

    waitbar(n / nRepeat, h);
end
close(h);

for j = 1 : length(Ep) - 1
    for dd = 1 : length(Ep{j})
        for sd = 1 : length(Ep{j}{dd})
            Ep{j}{dd}{sd} = sqrt(Ep{j}{dd}{sd} / (nRepeat - 1));
        end
    end
end

author: thomaskuestner


text:
function [ quasi, seed_new ] = sobol ( dim_num, seed )

%% SOBOL generates a new quasirandom Sobol vector with each call.
%
%  Discussion:
%
%    The routine adapts the ideas of Antonov and Saleev.
%
%  Modified:
%
%    30 March 2003
%
%  Reference:
%
%    Antonov and Saleev,
%    USSR Computational Mathematics and Mathematical Physics,
%    Volume 19, 1980, pages 252 - 256.
%
%    Paul Bratley and Bennett Fox,
%    Algorithm 659:
%    Implementing Sobol's Quasirandom Sequence Generator,
%    ACM Transactions on Mathematical Software,
%    Volume 14, Number 1, pages 88-100, 1988.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    I Sobol,
%    USSR Computational Mathematics and Mathematical Physics,
%    Volume 16, pages 236-242, 1977.
%
%    I Sobol and Levitan,
%    The Production of Points Uniformly Distributed in a Multidimensional
%    Cube (in Russian),
%    Preprint IPM Akad. Nauk SSSR,
%    Number 40, Moscow 1976.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%    DIM_NUM must satisfy 2 <= DIM_NUM <= 40.
%
%    Input/output, integer SEED, the "seed" for the sequence.
%    This is essentially the index in the sequence of the quasirandom
%    value to be generated.  On output, SEED has been set to the
%    appropriate next value, usually simply SEED+1.
%    If SEED is less than 0 on input, it is treated as though it were 0.
%    An input value of 0 requests the first (0-th) element of the sequence.
%
%    Output, real QUASI(DIM_NUM), the next quasirandom vector.
%
  global SOBOL_lastq;
  global SOBOL_seed;

  dim_max = 40;
%
%  Initialize (part of) V.
%
    v(1:40,1:30) = zeros(40,30);

    v(1:40,1) = [ ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]';

    v(3:40,2) = [ ...
            1, 3, 1, 3, 1, 3, 3, 1, ...
      3, 1, 3, 1, 3, 1, 1, 3, 1, 3, ...
      1, 3, 1, 3, 3, 1, 3, 1, 3, 1, ...
      3, 1, 1, 3, 1, 3, 1, 3, 1, 3 ]';

    v(4:40,3) = [ ...
               7, 5, 1, 3, 3, 7, 5, ...
      5, 7, 7, 1, 3, 3, 7, 5, 1, 1, ...
      5, 3, 3, 1, 7, 5, 1, 3, 3, 7, ...
      5, 1, 1, 5, 7, 7, 5, 1, 3, 3 ]';

    v(6:40,4) = [ ...
                     1, 7, 9,13,11, ...
      1, 3, 7, 9, 5,13,13,11, 3,15, ...
      5, 3,15, 7, 9,13, 9, 1,11, 7, ...
      5,15, 1,15,11, 5, 3, 1, 7, 9 ]';

    v(8:40,5) = [ ...
                           9, 3,27, ...
     15,29,21,23,19,11,25, 7,13,17, ...
      1,25,29, 3,31,11, 5,23,27,19, ...
     21, 5, 1,17,13, 7,15, 9,31, 9 ]';

    v(14:40,6) = [ ...
              37,33, 7, 5,11,39,63, ...
     27,17,15,23,29, 3,21,13,31,25, ...
      9,49,33,19,29,11,19,27,15,25 ]';

    v(20:40,7) = [ ...
                                         13, ...
     33,115, 41, 79, 17, 29,119, 75, 73,105, ...
      7, 59, 65, 21,  3,113, 61, 89, 45,107 ]';

    v(38:40,8) = [ ...
                                  7, 23, 39 ]';
%
%  Set POLY.
%
    poly(1:40)= [ ...
        1,   3,   7,  11,  13,  19,  25,  37,  59,  47, ...
       61,  55,  41,  67,  97,  91, 109, 103, 115, 131, ...
      193, 137, 145, 143, 241, 157, 185, 167, 229, 171, ...
      213, 191, 253, 203, 211, 239, 247, 285, 369, 299 ];
%
%  Check parameters.
%
    if ( dim_num < 2 | dim_max < dim_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SOBOL - Fatal error!\n' );
      fprintf ( 1, '  The spatial dimension DIM_NUM should satisfy:\n' );
      fprintf ( 1, '    2 <= DIM_NUM <= %d\n', dim_max );
      fprintf ( 1, '  But this input value is DIM_NUM = %d\n', dim_num );
      return
    end

    atmost = 2^30 - 1;
%
%  Find the number of bits in ATMOST.
%
    maxcol = bit_hi1_base_2 ( atmost );
%
%  Initialize row 1 of V.
%
    v(1,1:maxcol) = 1;
%
%  Initialize the remaining rows of V.
%
    for ( i = 2 : dim_num )
%
%  The bit pattern of the integer POLY(I) gives the form
%  of polynomial I.
%
%  Find the degree of polynomial I from binary encoding.
%
      j = poly(i);
      m = 0;

      while ( 1 )

        j = floor ( j / 2 );

        if ( j <= 0 )
          break;
        end

        m = m + 1;

      end
%
%  We expand this bit pattern to separate components of the logical array INCLUD.
%
      j = poly(i);
      for ( k = m : -1 : 1 )
        j2 = floor ( j / 2 );
        includ(k) = ( j ~= 2 * j2 );
        j = j2;
      end
%
%  Calculate the remaining elements of row I as explained
%  in Bratley and Fox, section 2.
%
      for ( j = m + 1 : maxcol )

        newv = v(i,j-m);
        l = 1;

        for ( k = 1 : m )

          l = 2 * l;

          if ( includ(k) )
            newv = exor ( newv, l * v(i,j-k) );
          end

        end

        v(i,j) = newv;

      end

    end
%
%  Multiply columns of V by appropriate power of 2.
%
    l = 1;
    for ( j = maxcol-1 : -1 : 1 )
      l = 2 * l;
      v(1:dim_num,j) = v(1:dim_num,j) * l;
    end
%
%  RECIPD is 1/(common denominator of the elements in V).
%
    recipd = 1.0E+00 / ( 2 * l );

  seed = floor ( seed );

  if ( seed < 0 )
    seed = 0;
  end

  if ( seed == 0 )

    l = 1;
    SOBOL_lastq(1:dim_num) = 0;

  elseif ( seed == SOBOL_seed + 1 )
%
%  Find the position of the right-hand zero in SEED.
%
    l = bit_lo0_base_2 ( seed );

  elseif ( seed <= SOBOL_seed )

    SOBOL_seed = 0;
    l = 1;
    SOBOL_lastq(1:dim_num) = 0;

    for ( seed_temp = SOBOL_seed : seed-1 )

      l = bit_lo0_base_2 ( seed_temp );

      for ( i = 1 : dim_num )
        SOBOL_lastq(i) = exor ( SOBOL_lastq(i), v(i,l) );
      end

    end

    l = bit_lo0_base_2 ( seed );

  elseif ( SOBOL_seed+1 < seed )

    for ( seed_temp = SOBOL_seed+1 : seed-1 )

      l = bit_lo0_base_2 ( seed_temp );

      for ( i = 1 : dim_num )
        SOBOL_lastq(i) = exor ( SOBOL_lastq(i), v(i,l) );
      end

    end

    l = bit_lo0_base_2 ( seed );

  end
%
%  Check that the user is not calling too many times!
%
  if ( maxcol < l )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SOBOL - Fatal error!\n' );
    fprintf ( 1, '  Too many calls!\n' );
    fprintf ( 1, '  MAXCOL = %d\n', maxcol );
    fprintf ( 1, '  L =      %d\n', l );
    return
  end
%
%  Calculate the new components of QUASI.
%
  for ( i = 1 : dim_num )

    quasi(i) = SOBOL_lastq(i) * recipd;

    SOBOL_lastq(i) = exor ( SOBOL_lastq(i), v(i,l) );

  end

  SOBOL_seed = seed;
  seed_new = seed + 1;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function DetectRedCells(db, ops0)

ops = build_ops3(db, ops0);
redcells = [];
for i = 1:length(ops.planesToProcess)
    iplane  = ops.planesToProcess(i);
    try
        try
            fname = sprintf('%s/F_%s_%s_plane%d.mat', ops.ResultsSavePath, ops.mouse_name, ops.date, iplane);
            dat = load(fname);
            while isfield(dat, 'dat')
                dat = dat.dat;
            end
        catch
            fname = sprintf('%s/F_%s_%s_plane%d_Nk%d_proc.mat', ops.ResultsSavePath, ops.mouse_name, ops.date, iplane, ops.Nk);
            dat = load(fname);
            while isfield(dat, 'dat')
                dat = dat.dat;
            end
        end
    catch
        fname = sprintf('%s/F_%s_%s_plane%d_Nk%d.mat', ops.ResultsSavePath, ops.mouse_name, ops.date, iplane, ops.Nk);
        dat = load(fname);
        while isfield(dat, 'dat')
            dat = dat.dat;
        end
    end

    % subtract bleedthrough of green into red channel
    % if there's a specific green channel (from a red+green short session)
    if isfield(dat.ops,'mimgGREEN')
        mimgG = dat.ops.mimgGREEN;
    else
        mimgG = dat.ops.mimg;
    end
    mimgR = dat.ops.mimgRED;

    % subsample mimgR and mimgG
    mimgR = mimgR(dat.ops.yrange,dat.ops.xrange);
    mimgG = mimgG(dat.ops.yrange,dat.ops.xrange);
    [nY nX] = size(mimgR);
    %
    % pixels with ROIs
    iscell = [dat.stat.mrs]./[dat.stat.mrs0]<1.2 & [dat.stat.npix]<150 & [dat.stat.npix]>20;
    icell  = find(iscell);
    cellpix = false(numel(mimgR),1);
    for j = icell
        cellpix(dat.stat(j).ipix) = 1;
    end

    nonrig = 1;
    if nonrig
        % non-rigid linear regression
        nblks               = 3;
        yB                  = round(linspace(1,nY,nblks+1));
        xB                  = round(linspace(1,nX,nblks+1));
        [yBL,xBL,numBlocks] = MakeQuadrants(yB,xB);
        msk                 = zeros(nY,nX,length(yB));
        for j = 1:length(yBL)
            xg              = mimgG(yBL{j},xBL{j});
            xr              = mimgR(yBL{j},xBL{j});
            A               = polyfit(xg,xr,1);
            msk(:,:,j)      = QuadrantMask(yBL,xBL,nY,nX,j);
            gw(j)           = A(1);
        end
        msk                 = bsxfun(@times,msk,1./sum(msk,3));
        pixweight           = zeros(nY,nX);
        for j = 1:length(yBL)
            pixweight = pixweight + msk(:,:,j)*gw(j);
        end
        mimgR0              = mimgR - pixweight.*mimgG;
    else
        % fit to non-ROI pixels
        mimgGn              = mimgG;
        mimgGn(cellpix)     = 0;
        mimgRn              = mimgR;
        mimgRn(cellpix)     = 0;

        xg                  = mimgGn;%my_conv2(mimgGn,[1 1],[1 2]);
        xr                  = mimgRn;%my_conv2(mimgRn,[1 1],[1 2]);
        A                   = polyfit(xg(:),xr(:),1);
        plot(xg(1:100:end)*A(1)+A(2),xr(1:100:end),'.')
        mimgR0              = mimgR - (A(1)*mimgG);
    end
    dat.ops.mimgREDcorrected = mimgR0;

    %%
    mimgR0 = normalize_image(mimgR0);
    %%%%% compute overlap with pixel map
    % (exclude pixels from other cells (cellpix))

    redpix = zeros(length(dat.stat),2);
    [xx,yy] = meshgrid([1:size(mimgR0,1)],[1:size(mimgR0,2)]);
    xx=xx(:); yy=yy(:);
    xy2 = xx.^2+yy.^2;

    for j = 1:numel(dat.stat)
        ipix                 = dat.stat(j).ipix;
        imgpix               = false(size(mimgR0));
        imgpix(ipix)         = 1;
        rpix                 = mimgR0(ipix);%.*(dat.stat(j).lambda'/sum(dat.stat(j).lambda));
        [ix,iy]              = ind2sub(size(mimgR0),ipix);
        params               = FitCircle(ix,iy);
        cellradius           = params.ra;
        yc                   = params.xc;
        xc                   = params.yc;
        icircle              = (xy2 - 2*xc*xx - 2*yc*yy + xc^2 + yc^2) < (cellradius+15)^2;
        extpix               = icircle & ~imgpix(:) & ~cellpix;
        ext_rpix             = mimgR0(extpix);%/length(extpix);
        redpix(j,:)          = [mean(rpix) mean(ext_rpix)];
    end
    redpix = redpix - min(redpix(:));
%     redpix(~dat.cl.isroi,:) = NaN;

    % set threshold for redpix
    if isfield(ops,'redthres')
        redthres = ops.redthres;
    else
        redthres = 1.5;
    end
    if isfield(ops,'redmax')
        redmax = ops.redmax;
    else
        redmax = 1;
    end
    rrat = redpix(:,1)./(redpix(:,2)+redpix(:,1));
    redcell  = rrat > nanmean(rrat)+redthres*nanstd(rrat);
    notred   = rrat <= nanmean(rrat) + redmax*nanstd(rrat);
    %%
    fprintf('plane %d  reds %d\n',iplane,sum(redcell(:)&iscell(:)));

%     dat.cl.redcell = redcell(:);
%     dat.cl.notred  = notred(:);

    for j = 1:length(dat.stat)
        dat.stat(j).redcell = redcell(j);
        dat.stat(j).redprob = rrat(j);
    end

    save(fname, '-struct', 'dat')
end

author: cortex-lab


text:
function [ res, jac ] = opt05_rj ( x, flag, a )

%*****************************************************************************80
%
%% OPT05_RJ evaluates RES and JAC for test case #5.
%
%  Discussion:
%
%    This example is known as the Rosenbrock "banana" function.
%
%    This example is discussed in Dennis and Schnabel, page 157.
%    (There is a misprint in Dennis and Schabel.  When they print the
%    rescaled function, they include an incorrect sign on the second term.)
%
%    This example tests the scaling matrices.
%
%    Suggested starting points are
%
%      X(init) = ( -1.2/A, A)
%      or
%      X(init) = ( 6.39/A, -0.221*A )
%
%    The optimizing value is
%
%      X* = (1/A, A)
%
%    for which
%
%      RES(X*) = (0,0).
%
%    A typical value of A is 1.  For values of A that are greater or lesser
%    than 1, the optimization can take significantly longer.
%
%  Modified:
%
%    02 January 2008
%
%  Author:
%
%    Jeff Borggaard,
%    Gene Cliff,
%    Virginia Tech.
%
%  Reference:
%
%    John Dennis, Robert Schnabel,
%    Numerical Methods for Unconstrained Optimization
%    and Nonlinear Equations,
%    SIAM, 1996,
%    ISBN13: 978-0-898713-64-0,
%    LC: QA402.5.D44.
%
%  Parameters:
%
%    Input, real X(2), the evaluation point.
%
%    Input, string FLAG, indicates what must be computed.
%    'f' means only the value of RES is needed,
%    'g' means only the value of JAC is needed,
%    'all' means RES and JAC are needed.
%    It is acceptable to behave as though FLAG was 'all'
%    on every call.
%
%    Input, real A, the scale factor.  A typical value is 1.
%    A should not be 0.
%
%    Output, real RES(2,1), the residual column vector.
%
%    Output, real J(2,2), the Jacobian matrix.
%
  n = length ( x );

  if ( n ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'OPT05_RJ - Fatal error!\n' );
    fprintf ( '  The input vector X should have length 2.\n'),
    fprintf ( '  Instead, it has length = %d.\n', n );
    keyboard
  end

  if ( a == 0 )
    fprintf ( '\n' );
    fprintf ( 'OPT05_RJ - Fatal error!\n' );
    fprintf ( '  The scale factor A should be nonzero.\n'),
    keyboard
  end

  res(1,1) = 10 * ( ( a * x(1) )^2 - x(2) / a );
  res(2,1) = 1 - a * x(1);

  jac(1,1) = 20 * a * a * x(1);
  jac(1,2) = - 10 / a;

  jac(2,1) = - a;
  jac(2,2) = 0.0;


text:
classdef dme_xfmr3p_opf < mp.dme_xfmr3p & mp.dme_shared_opf

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

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% effOFexpDesignsNoOverlap2
% ER-fMRI data analysis
% script:
% throughout we assume TR=2.

loadDataYes=1; % =1 -> loads existing data file

%cd /home/giedrius/giedrius/projects/matlab/erfmri
pluginCorrYes=0;
overlapYes=0;
nReps=10000; %
nVals=3;
n2=12; % assumed # TRs to cover HDR fn
pwRange=6:10;

daleYes=1; % 1-usual Dale efficiency. 0-Fisher
% event matrices are defined based on vectors with the average
% energy removed

cutEndPad=1; %cutting reduces efficiency!

if ~loadDataYes,
rndEffAll=[];
mEffAll=[];

tic
for pwr=pwRange,
% creating event vectors
ms=m2bin([mseq(2,pwr,0,1)])'; % fixing to the length of 2^n
n=length(ms)  % scan duration in TRs
if pluginCorrYes
  b=[0.406;0.8825]; % parameters from SPG fMRI data
  fittedACorr=autocorrFnct(b,1:n+(n2-1)*rem(cutEndPad+1,2));
  Cninv=inv(toeplitz(fittedACorr));
else
  Cninv=eye(n+(n2-1)*rem(cutEndPad+1,2));
end;

mRange=(2^pwRange(1)-2);
mEff=zeros(1,mRange);

%shift1=ceil(2^pwRange(1)/2);
%XXXXXXXXXXXXXXXXXXXXX
for shift=1:mRange,
	%shift=2^(pwr-1)+5;
	foo=[ms(shift+1:end) ms(1:shift)];%shifted version
	%foo1=[ms(shift1+1:end) ms(1:shift1)];%shifted version
	foo=ms+foo*2; %+foo1*4;

	mEvent=zeros(nVals,n);
	for k=1:nVals,
   	  ind=find(foo==k);
   	  mEvent(k,ind)=1;
	end;

	% defining event matrix as convolution matrix
	%eventMatrix=makeEventMtrx(mEvent-(ones(length(mEvent),1)*...
	%sum(mEvent'))'/length(mEvent),n2);

	eventMatrix=makeEventMtrx(mEvent,n2);

        if cutEndPad, eventMatrix=eventMatrix(1:n,:); end;
	eventMatrix=eventMatrix-ones(size(eventMatrix,1),1)* ...
	    sum(eventMatrix)/size(eventMatrix,1);

   if daleYes,
	  designEff=1/trace(inv(eventMatrix'*Cninv*eventMatrix));
	else
	  designEff=trace(eventMatrix'*Cninv*eventMatrix);
	end;

	mEff(shift)=designEff;
end;

%mProbab=sum(mEvent'); mProbab=mean(mProbab)/n;
plot(mEff,'k','LineWidth',2);
set(gca,'FontSize',12);
xlabel('Cyclical shift of event vector #3','FontSize',16)
ylabel('Efficiency','FontSize',16)
set(gca,'Position',[.2,.15,.7,.7]);
set(gcf,'PaperPosition',[1,1,5,4]);
%axis([0 n 0 .5])
drawnow;
%print -dpsc2 cyclingNonOverlapping
pause

mEffAll=[mEffAll, max(mEff)];

% for simply randomized designs:
nVals=size(mEvent,1);
%n=n+1
if pluginCorrYes
  b=[0.406;0.8825]; % parameters from SPG fMRI data
  fittedACorr=autocorrFnct(b,1:n+(n2-1)*rem(cutEndPad+1,2));
  Cninv=inv(toeplitz(fittedACorr));
else
  Cninv=eye(n+(n2-1)*rem(cutEndPad+1,2));
end;

if pwr<8, correction=0.02;
elseif pwr<9, correction=0.005;
elseif pwr==10, correction=0.001;
end;
correction=0;

pEv=nVals/(nVals+1)+correction; %

rndEff=zeros(1,nReps);
for k=1:nReps;
   rndEvent=balancedRnd(n,nVals,pEv,overlapYes);%sum(rndEvent')

   % defining event matrix as convolution matrix
   %eventMatrix=makeEventMtrx(rndEvent-(ones(length(rndEvent),1)*...
   %sum(rndEvent'))'/length(rndEvent),n2);

   eventMatrix=makeEventMtrx(rndEvent,n2);

   if cutEndPad, eventMatrix=eventMatrix(1:n,:); end;
   eventMatrix=eventMatrix-ones(size(eventMatrix,1),1)* ...
       sum(eventMatrix)/size(eventMatrix,1);

   if daleYes,
	  designEff=1/trace(inv(eventMatrix'*Cninv*eventMatrix));
   else
	  designEff=trace(eventMatrix'*Cninv*eventMatrix);
   end;
   rndEff(k)=designEff;
end;

rndEffAll=[rndEffAll;rndEff];
end;
toc;

else
  eval(['load dataEffNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh']);
end %fi loadDataYes

nn=2.^[pwRange]-1;

%[b,a]=hist(rndEff,20);
%bar(a,b/nReps); hold on;
%plot(max(mEff),0,'k*');

figure(1); clf;
loglog(nn,mEffAll,'r*');hold on;

% theoretical max -- only for one event type:
%theoMax=(nn+1)/n2/4;
%loglog(nn,theoMax,'g+');hold on;

p99=prctile(rndEffAll',99.9);
p00=prctile(rndEffAll',0.1);
med=median(rndEffAll');
loglog(nn,med,'k.-','LineWidth',2);
loglog([nn;nn],[p00;p99],'k')
loglog([nn*.93;nn*1.1],[p00;p00],'k')
loglog([nn*.93;nn*1.1],[p99;p99],'k')

axis([50 10^3*1.2 10^-2 10^1]);
%set(gca,'YTick',[0.003 .01 .1 .3 1 3 10 30]);
set(gca,'XTick',nn);
set(gca,'LineWidth',2,'FontSize',12);

xlabel('Sequence length','FontSize',12);
ylabel('Efficiency','FontSize',12);
title(['No Overlap: n_e=',num2str(nVals),' p=',num2str(pEv),' n_h=',num2str(n2)]);

set(gca,'Position',[.2,.15,.7,.7]);
set(gcf,'PaperPosition',[1,1,4,3]);
eval(['print -dpsc2 effNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh']);
disp(['print -dpsc2 effNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh']);

figure(2); clf;
semilogx(nn,mEffAll./max(rndEffAll'),'-ok','LineWidth',2); hold on;
semilogx(nn,mEffAll./med,':+k','LineWidth',2)
set(gca,'Position',[.2,.15,.7,.7]);
xlabel('Sequence length','FontSize',12);
ylabel('m-seq/random efficiency ratio','FontSize',12);
title(['No Overlap: n_e=',num2str(nVals),' p=',num2str(pEv),' n_h=',num2str(n2)]);
set(gca,'XTick',nn)
set(gca,'LineWidth',2,'FontSize',12);
set(gca,'Position',[.2,.15,.7,.7]);
set(gcf,'PaperPosition',[1,1,4,3]);
axis([50 10^3*1.2 .9 2.5]);
eval(['print -dpsc2 effNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nhRATIO']);
disp(['print -dpsc2 effNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nhRATIO']);

eval(['save dataEffNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh rndEffAll p99 p00 med mEffAll']);
disp(['save dataEffNOV',num2str(nVals),'ev',num2str(nReps),'reps',num2str(n2),'nh rndEffAll p99 p00 med mEffAll']);


text:
classdef VehicleSimpleLinear < VehicleDynamicsLateral.VehicleSimple
    % VehicleSimpleLinear Linear simple vehicle model.
    %
    % It inherits properties from VehicleSimple.

    methods
        function self = VehicleSimpleLinear()
            % Constructor for the vehicle
            self.mF0 = 700;
            self.mR0 = 600;
            self.IT = 10000;
            self.lT = 3.5;
            self.nF = 2;
            self.nR = 2;
            self.wT = 2;
            self.muy = .8;
            self.deltaf = 0;
            self.Fxf = 0;
            self.Fxr = 0;
        end

        %% Model
        % Function with the model
        function dx = Model(self, t, states,tspan)
            % Data
            mT = self.mT;
            IT = self.IT;
            a = self.a;
            b = self.b;
            nF = self.nF;
            nR = self.nR;
            muy = self.muy;

            g = 9.81;                 % Gravity [m/s^2]

            FzF = self.mF0 * g;       % Vertical load @ F [N]
            FzR = self.mR0 * g;       % Vertical load @ R [N]

            v0 = 20;                  % [m/s]

            % State variables
            X = states(1,1);         % Not used
            Y = states(2,1);         % Not used
            PSI     = states(3,1);
            VT       = states(4,1);
            ALPHAT  = states(5,1);
            dPSI    = states(6,1);

            if isa(self.deltaf,'function_handle')
                deltaf = self.deltaf([X;Y;PSI;VT;ALPHAT;dPSI],t);
            elseif length(self.deltaf)>1
                deltaf = interp1(tspan,self.deltaf,t);
            else
                deltaf = self.deltaf;
            end

            % Slip angles
            ALPHAF = ALPHAT + a/v0*dPSI - deltaf;
            ALPHAR = ALPHAT - b/v0*dPSI;

            % Longitudinal forces
            if isa(self.Fxf,'function_handle')
                FxF = self.Fxf([X;Y;PSI;VT;ALPHAT;dPSI],t);
            elseif length(self.Fxf)>1
                FxF = interp1(tspan,self.Fxf,t);
            else
                FxF = self.Fxf;
            end

            if isa(self.Fxr,'function_handle')
                FxR = self.Fxr([X;Y;PSI;VT;ALPHAT;dPSI],t);
            elseif length(self.Fxr)>1
                FxR = interp1(tspan,self.Fxr,t);
            else
                FxR = self.Fxr;
            end

            % Lateral force
            FyF = nF * self.tire.Characteristic(ALPHAF, FzF / nF, muy);
            FyR = nR * self.tire.Characteristic(ALPHAR, FzR / nR, muy);

            % State equations
            dx(1,1) = VT;
            dx(2,1) = v0*(PSI + ALPHAT);
            dx(3,1) = dPSI;
            dx(4,1) = (FxF + FxR)/mT;
            dx(5,1) = (FyF + FyR)/(mT*v0) - dPSI;
            dx(6,1) = (a*FyF - b*FyR)/IT;

        end
    end
end

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

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function z = r82_to_c8 ( a, b )

%*****************************************************************************80
%
%% R82_TO_C8 makes a complex number from two R8's.
%
%  Discussion:
%
%    It is extremely unnatural for a programmer to find that
%    "i" is, by default, the imaginary unit in MATLAB.
%
%    So the only safe thing to do, in my opinion, is to hide
%    it somewhere, as in this function!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the real and imaginary parts.
%
%    Output, complex Z, the value of ( A + B * i ).
%
  z = a + b * i;

  return
end


text:
function seq = report_tracking_result(seq, result)

if strcmpi(seq.format, 'otb')
    seq.rect_position(seq.frame,:) = round([result.center_pos([2,1]) - (result.target_size([2,1]) - 1)/2, result.target_size([2,1])]);
elseif strcmpi(seq.format, 'vot')
    if seq.frame > 1
        bb_scale = 1;
        sz = result.target_size / bb_scale;
        tl = result.center_pos - (sz - 1)/2;
        br = result.center_pos + (sz - 1)/2;
        x1 = tl(2); y1 = tl(1);
        x2 = br(2); y2 = br(1);
        result_box = round(double([x1 y1 x2 y1 x2 y2 x1 y2]));
        if any(isnan(result_box) | isinf(result_box))
            error('Illegal values in the result.')
        end
%         if any(result_box < 0)
%             error('Negative values')
%         end
        seq.handle = seq.handle.report(seq.handle, result_box);
    end
else
    error('Uknown sequence format');
end

author: martin-danelljan

repo: ECO


text:
function test_ft_volumedownsample

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_volumedownsample SPM

mri = [];
mri.anatomy = randn(181,217,181);
mri.dim = [181 217 181];
mri.transform = eye(4);
mri.unit = 'cm';
mri.coordsys = 'ctf';

cfg = [];
cfg.downsample = 2;
cfg.parameter = 'anatomy';
mriout = ft_volumedownsample(cfg, mri);