text:
function h = labelpoints (xpos, ypos, labels, varargin)
%  h = labelpoints (xpos, ypos, labels, position, buffer, adjust_axes, {parameters})
%
%   Given x and y position vectors (xpos, ypos) and given
%     a vector of labels this script will label all data points
%     and output the label handles in vector, h.
%
%   'xpos' and 'ypos' are required to be the same length. However, if all labels
%     fall along the same horizontal or vertical line, the function will accept
%     a single number for that xpos or ypos (see examples).
%
%   'labels' should be a cell or numerical array of the
%     same length as xpos or ypos.  Alternatively it can be a singleton
%     that will be replicated for all labels (see examples).
%
%   'position' (optional) describes the position of the labels
%     relative to their locations by entering one of the following
%     abbreviated compass directions in single quotes (default is 'NW').
%     N, S, E, W, NE, NW, SE, SW, center
%
%   'buffer' (optional) is a number between 0:1 that adds distance between
%     the label and the plotted point where 0 (default) is none and 1 is
%     1/10th of the axis.  Ignored for 'center' position.
%
%   'adjust_axes' (optional, default=0): depending on the positioning of labels,
%     some may fall beyond the axis limits. adjust_axes = 1 will readjust xlim
%     & ylim slightly so labels are not beyond axis limits.
%
%   This function includes three optional parameters to label only the outliers.
%     This can come in handy when there are many points but only the outliers
%     should be labeled. (may require stats toolbox)  (see examples)
%       'outliers_SD', N   -  will only label points that are greater than N
%                             standard deviations from the median of xpos or ypos.
%       'outliers_Q', N    -  will only label points that are greater than N
%                             times the interquartile range of xpos or ypos.
%       'outliers_N', N    -  will calculate the distance of each point from the median
%                             point(xpos,ypos) and will label the N-furthest points.
%                             Alternatively, N can be a decimal to label N% of the points.
%
%   The following parameters may also be entered in any order
%     'FontSize', N -   font size of all labels
%     'Color', S    -   font color of all labels ('y' 'm' 'c' 'r' 'g' 'b' 'w' 'k')
%     'rotation', N -   will rotate the labels N degrees about label-center
%                       (positive is counterclockwise).
%
%
%   Examples:
%      Fake Data:
%       x = 1:10;  y=rand(1,10);
%       scatter(x,y)
%       labs = {'a' 'b' 'c' 'd' 'middle' 'f' 'g' 'h' 'i' 'last'};
%
%      Label Examples
%       txt_h = labelpoints(x, y, labs);
%       txt_h = labelpoints(x, y, labs, 'E', 0.15);
%       txt_h = labelpoints(x, y, labs, 'E', 0.15, 1);
%       txt_h = labelpoints(x, y, labs, 'W', 0.15, 1, 'FontSize', 14, 'Color', 'r');
%       txt_h = labelpoints(x, y, labs, 'W', 0.15, 1, 'FontSize', 12, 'Color', 'm', 'rotation', 45);
%
%      Also works for
%       labs = [1:1:10];
%       labs = {'Sofia' '' '' '' 'Bucharest' '' '' 'Belgrade' '' 'Ankara'}
%       labs = '*';
%       labs = 'string';
%
%      When all labels share same xpos or ypos
%         boxplot(1:10)
%         labelpoints(0.8, [3, 5.5, 8], {'25%' '50%' '75%'}, 'center');
%
%      Outlier Examples
%       Fake Data:
%         x = [rand(1,30), rand(1,8)*2];
%         y = [rand(1,30), rand(1,8)*2];
%         scatter(x, y)
%         labs = 1:38;
%        labelpoints(x, y, labs, 'N', 0.1, 1, 'outliers_N', 5);                  %will label 5 furthest points from median
%        labelpoints(x, y, labs, 'N', 0.1, 1, 'outliers_N', 0.1, 'Color', 'r');  %will label 10% of furthest points from median
%        labelpoints(x, y, labs, 'N', 0.1, 1, 'outliers_SD', 2);                 %will label all points > 2 SD from median
%        labelpoints(x, y, labs, 'N', 0.1, 1, 'outliers_Q', 1.5);                %will label points greater that 1.5 x IQR
%
%   Alternative use:
%     Density Distributions:
%       x = randn(1,100); y = 1:100;
%       scatter(x,y)
%       labs = '|';
%       txt_h = labelpoints(x, 8, labs, 'center');
%
%     Single Labeling
%       x = 2004:2013;  y=rand(1,10);
%       plot(x,y,'-o')
%       labs = 'acquisition';
%       labelpoints(x(3), y(3), labs, 'N', 0.2, 1);
%       labelpoints(2008.5, min(ylim), {['labelpoints.m   ', datestr(now, 'mm/dd/yy')]}, 'N', 0.3, 0, 'fontsize', 12, 'color', 'm');
%
%     Use labels instead of markers
%       x = randn(1,15); y = randn(1,15);
%       labs = char('a'+(1:15)-1)';
%       labelpoints(x, y, labs, 'center', 0, 1, 'color', 'b');
%
%  140331 v.1
%  141115 v.2
% Copyright (c) 2014, Adam Danz
%All rights reserved

% source: http://www.mathworks.com/matlabcentral/fileexchange/46891-labelpoints

% Changes history
%   11/02/14    if labels are entered as char, they are convered to cell
%   11/03/14    varargin added to accomodate auto_outlier feature
%   11/03/14    convert inputs to row vector if they are column vectors
%   11/04/14    now position in put is not case sensitive
%   11/04/14    now 1 label can be entered to label multiple points
%   11/04/14    now fontsize and color can be specified by params
%   11/05/14    changed 'outlier' to 'outlier_SD' and added 'outlier_Q'
%   11/07/14    added option to rotate test
%   11/15/14    added 'outliers_N' input option and cleaned up text rotation section.
%   11/19/14    curr_extent is not always a cell.  fixed.
%   11/20/14    when outliers_N is selected N is min(N, lenght(xpos));
%   11/21/14    removes entire point and label when there is an 'inf' value

%%
% Check Class of 'labels'
    %If 'labels' are numberical, convert to cell
    if isnumeric(labels) == 1
        labels = num2cell(labels);
    end

    % if 'labels' are char, convert to cell
    if ischar(labels)
        labels = cellstr(labels);
    end

% if all labels share the same xpos or ypos (only 1 value entered in 1 of the position vectors)
    if length(xpos)==1 && length(ypos)>1
        xpos = repmat(xpos, size(ypos));
    elseif length(ypos)==1 && length(xpos)>1
        ypos = repmat(ypos, size(xpos));
    end

% if only one lable is entered for all points, replicate it
    if length(labels)==1 && length(xpos) > 1
        labels = repmat(labels, [1, length(xpos)]);
    end

% ensures xpos, ypos, and labels are all row vectors
    if iscolumn(xpos);      xpos = xpos';       end
    if iscolumn(ypos);      ypos = ypos';       end
    if iscolumn(labels);    labels = labels';   end

%check that x, y, and labels are same length
    if isequal(length(xpos), length(ypos), length(labels)) == 0
        error('xpos, ypos, and labels must all be the same length unless using one input for labels.')
    end

%if an 'inf' value is entered, this will remove that entire point and label
    xinf = find(xpos==inf);
    yinf = find(ypos==inf);
    findinf = [xinf yinf];
    if ~isempty(findinf)
        xpos(findinf)=[];
        ypos(findinf)=[];
        labels(findinf) = [];
    end

%Validate inputs and optional parameters
    validPositions = {'N' 'NE' 'E' 'SE' 'S' 'SW' 'W' 'NW' 'center'};
    checkPosition = @(x) any(validatestring(x, validPositions));

    validColors = {'y' 'm' 'c' 'r' 'g' 'b' 'w' 'k'};
    checkColors = @(x) any(validatestring(x, validColors));

    p = inputParser;
    p.FunctionName = mfilename;
    addRequired(p, 'xpos', @isnumeric);
    addRequired(p, 'ypos', @isnumeric);
    addRequired(p, 'labels');
    addOptional(p, 'position', 'NW', checkPosition);
    addOptional(p, 'buffer', 0, @isnumeric);
    addOptional(p, 'adjust_axes', 0, @isnumeric);

    addParamValue(p, 'outliers_SD', 3, @isnumeric);
    addParamValue(p, 'outliers_Q', 1.5, @isnumeric);
    addParamValue(p, 'outliers_N', 1, @isnumeric);

    addParamValue(p, 'FontSize', 10, @isnumeric);
    addParamValue(p, 'Color', 'k', checkColors);
    addParamValue(p, 'rotation', 0, @isnumeric);
    parse(p, xpos, ypos, labels, varargin{:})

mfile = [mfilename,'.m'];

%calculate buffer
    a = axis/10;% I've somewhat arbitrarily divided by 10 to make 'buffer' more sensitive
    u1 = 0;     %x offset
    u2 = 0;     %y offset

%assign position
    switch upper(p.Results.position)
        case 'E',       va = 'middle'; ha = 'left';         u1 = a(2)-a(1);
        case 'W',       va = 'middle'; ha = 'right';        u1 = (a(2)-a(1))*-1;
        case 'N',       va = 'bottom'; ha = 'center';       u2 = a(4)-a(3);
        case 'S',       va = 'top'; ha = 'center';          u2 = (a(4)-a(3))*-1;
        case 'NE',      va = 'bottom'; ha = 'left';         u1 = (a(2)-a(1))/2;     u2 = (a(4)-a(3))/2;
        case 'NW',      va = 'bottom'; ha = 'right';        u1 = (a(2)-a(1))*-0.5;  u2 = (a(4)-a(3))/2;
        case 'SE',      va = 'top'; ha = 'left';            u1 = (a(2)-a(1))/2;     u2 = (a(4)-a(3))*-0.5;
        case 'SW',      va = 'top'; ha = 'right';           u1 = (a(2)-a(1))*-0.5;  u2 = (a(4)-a(3))*-0.5;
        case 'CENTER',  va = 'middle'; ha = 'center';
    end

%Factor in buffer
    u1 = u1*p.Results.buffer;
    u2 = u2*p.Results.buffer;

%If outliers parameters are selected
    if sum(strcmp(varargin, 'outliers_SD')) == 1 || sum(strcmp(varargin, 'outliers_Q')) == 1 || sum(strcmp(varargin, 'outliers_N')) == 1
        if sum(strcmp(varargin, 'outliers_SD')) == 1
            outlier_idx = logical(abs(xpos - median(xpos)) > p.Results.outliers_SD*std(xpos)  |  abs(ypos - median(ypos)) > p.Results.outliers_SD*std(ypos)); %index of outliers

        elseif sum(strcmp(varargin, 'outliers_Q')) == 1
            xbounds = [prctile(xpos,25) - p.Results.outliers_Q * iqr(xpos) , prctile(xpos, 75) + p.Results.outliers_Q * iqr(xpos)];   %[lower upper] bounds of outliers
            ybounds = [prctile(ypos,25) - p.Results.outliers_Q * iqr(ypos) , prctile(ypos, 75) + p.Results.outliers_Q * iqr(ypos)];   %[lower upper] bounds of outliers
            outlier_idx = logical(ypos<ybounds(1) | ypos>ybounds(2) |  xpos<xbounds(1) | xpos>xbounds(2));

        elseif sum(strcmp(varargin, 'outliers_N')) == 1
            if p.Results.outliers_N<1;
                N = round(length(xpos) * p.Results.outliers_N);
            else
                N = min(p.Results.outliers_N, length(xpos));        %ensures that user cannot label more outliers than coordinates.
            end
             medianpoint = repmat([median(xpos) median(ypos)], [length(xpos),1]);
             paired = horzcat(xpos', ypos');
             distances = (((medianpoint(:,1)-paired(:,1)).^2)  +  ((medianpoint(:,2)-paired(:,2)).^2)).^(1/2);       %all distances from median
             [~, idx] = sort(distances, 'descend');
             outlier_idx = false(1,length(xpos));
             outlier_idx(idx(1:N))=1;
        end

        xpos = xpos(outlier_idx);
        ypos = ypos(outlier_idx);
        labels = labels(outlier_idx);

        if any(outlier_idx) == 0;           %dispay msg if there are no outliers to label
            disp(['There are no outliers to label in ', mfile,'.'])
            disp('Change outlier value for less sensitivity; See help file.');
        end
    end

%Label points
    h = text(xpos+u1 , ypos+u2, labels, 'VerticalAlignment',va, 'HorizontalAlignment',ha, 'FontSize', p.Results.FontSize, 'color', p.Results.Color);

%Rotate text if specified
    if sum(strcmp(varargin, 'rotation')) == 1               %if rotation parameter is called in inputs
        xl = xlim;      yl = ylim;                          %In case text rotation auto adjusts axes.
        curr_extent = get(h, 'extent');                     %Need to store current center point of all labels since text rotation relocates position
        if iscell(curr_extent); cell2mat(curr_extent); end
        hold on
        curr_position = [curr_extent(:,1)+(curr_extent(:,3)/2),curr_extent(:,2)+(curr_extent(:,4)/2)];          %uses extent to locate center of label
        set(h, 'rotation', p.Results.rotation, 'VerticalAlignment','middle', 'HorizontalAlignment','center');  	%note: text rotation changes alignment which is why they need to be centered back to specifications.
        for i = 1:length(h)                                 %after rotation, reposition lables back to desired location
            set(h(i), 'position', curr_position(i,:))
        end
        set(gca, 'xlim', xl); set(gca, 'ylim', yl);         %In case text rotation auto adjusts axes.
    end

%Determine if any labels go beyond axis limits and adjust if desired  (adjust_axes = 0 or 1)
    if p.Results.adjust_axes == 1   &&   ~isempty(h)
        x_adj = sign(u1+0.0000001);                 %the addition is to avoid '0'
        y_adj = sign(u2+0.0000001);                 %the addition is to avoid '0'

        labelextent = get(h, 'extent');
        if isequal(class(labelextent),'cell')
           labelextent = cat(1, labelextent{:});
        end
        xl = xlim;      yl = ylim;
        lablimX = [min(labelextent(:,1)), max(labelextent(:,1)+(labelextent(:,3).*x_adj))] +u1;
        lablimY = [min(labelextent(:,2)), max(labelextent(:,2)+(labelextent(:,4).*y_adj))] +u2;

        xlim([min(min(xl), min(lablimX)), max(max(xl), max(lablimX))])
        ylim([min(min(yl), min(lablimY)), max(max(yl), max(lablimY))])
    end

end

%% Notes
% a video (not mine) explaining this method:  http://blogs.mathworks.com/videos/2012/05/30/how-to-label-a-series-of-points-on-a-plot-in-matlab/
% Text properties :  http://www.mathworks.com/help/matlab/ref/text_props.html
% info on input parsing:  http://www.mathworks.com/help/matlab/matlab_prog/parse-function-inputs.html
%                   and   http://www.mathworks.com/help/matlab/ref/inputparser-class.html
% Outlier info : https://docs.oracle.com/cd/E17236_01/epm.1112/cb_statistical/frameset.htm?ch07s02s10s01.html

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [Hdraw,HvarsDraw,phi_Hdraw,Adraw,h0]=sampleH_MHstep(yData,Psi,A,B,phi_H,HvarsOld,priorValues,dataValues,n,h0)

%yData=YData
%Psi=PsiDraw_prop
%A=Adraw
%B=Bdraw
%phi_H=phi_Hdraw'
%HvarsOld=HvarsDraw
%% Initialize
offset_c=priorValues.offset_c; %constant for log transformation

[T,M]=size(yData);
p = size(B,1)/M; %lags of the B vector that has dimensions MxM*p (no constant), each column is one regression

% obtain prior data
startMeanVector=priorValues.mean_ln_h0;  %log mean of the initial state (variance scaling parameters) as residuals from an AR(4) in the training sample
startVarVector=priorValues.var_ln_h0;    %variance of the initial state for the elements of lambda

% obtain prior data
b0=priorValues.mean_ln_h0;  %log mean of the initial state (variance scaling parameters) as residuals from an AR(4) in the training sample
a0=eye(M)/priorValues.var_ln_h0(1,1);    %variance of the initial state for the elements of lambda

priorPhi_H=priorValues.phi_h;            %centering parameter for the inverse gamma of the variance of the innovations governing the random walk for lambda
priorD_H  =priorValues.d_h;              %scaling parameter for the inverse gamma distribution for the variance of the innovations governing the random walk for lambda

kappa= priorValues.kappa;
gamma= 1; %priorValues.gamma;

S_h = priorPhi_H*ones(n,1);
nu_h = priorD_H*ones(n,1);

Y_Psi=yData-Psi;  %subtract the local mean
%Y_Psi(1:p,:)=yData(1:p,:)-ones(p,1)*mean(yData(1:p,:)); %also generate initial conditions for the construction of the lagmatrix
Y_Psi(1:p,:)=yData(1:p,:)-ones(p,1)*(Psi(p+1,:)); %also generate initial conditions for the construction of the lagmatrix

% X_Psi = lagmatrix(Y_Psi,1:p); %create RHS of the VAR part
% X_Psi = X_Psi(p+1:end,:);     %remove the first p rows of RHS
X_Psi = bear.lagx(Y_Psi,p-1);
X_Psi = X_Psi(1:end-1,:);
Y_Psi=Y_Psi(p+1:end,:);       %and do so for LHS

E=Y_Psi-X_Psi*B;              %VAR residuals
Escaled=E*A';                 %transform such that VAR residuals have variance Lambda

   for jj=1:T-p
   epst(:,:,jj)=E(jj,:)'; %generate period specific errors
   end

   scaling = ones(n,1);
   Abelowdiag = cell(n,1);
   for kk=2:size(E,2)
       Abelowdiag{kk,1} = A(kk,1:kk-1)';
   end

  %logHvarsOld = HvarsOld;
   logHvarsOld = log(HvarsOld);
   logHdraw = logHvarsOld(p+1:end,:);
%% draw the series lambda_i,t from their conditional posteriors, i=1,...,n and t=1,...,T

% consider variables in turn
   for jj=1:n
      % consider periods in turn
      for kk=p+1:T
      % a candidate value will be drawn from N(lambdabar,phibar)
      % the definitions of lambdabar and phibar varies with the period, thus define them first
         % if the period is the first period
         if kk==p+1
         %lambdabar=(gamma/(1+gamma^2))*((startMeanVector(jj,1))+logHvarsOld(kk+1,jj));
         lambdabar=(gamma*(h0(jj,1)+logHvarsOld(kk+1,jj)))/(1/kappa+gamma^2);
         phibar=phi_H(1,jj)/(1/kappa+gamma^2);
         %phibar=startVarVector(jj,1)+phi_H(1,jj)/(1/kappa+gamma^2);
         %if the period is the final period
         elseif kk==T
         lambdabar=gamma*logHdraw(end-1,jj);
         phibar=phi_H(1,jj);
         % if the period is any period in-between
         else
         lambdabar=(gamma/(1+gamma^2))*(logHdraw(kk-p-1,jj)+logHvarsOld(kk+1,jj));
         phibar=phi_H(1,jj)/(1+gamma^2);
         end
      % now draw the candidate
      cand=lambdabar+phibar^0.5*randn;

%       end
      % compute the acceptance probability
      prob=bear.mhprob2(jj,cand,logHvarsOld(kk,jj),scaling(jj,1),epst(:,1,kk-p),Abelowdiag{jj,1});
      % draw a uniform random number
      draw=rand;
         % keep the candidate if the draw value is lower than the prob
         if draw<=prob
         logHdraw(kk-p,jj)=cand;
%          else
%          logHdraw(kk-p,jj)=logHvarsOld(kk,jj);
%          % if not, just keep the former value
          end
      end
   end

G=speye(T-p)-sparse(diag(gamma*ones(T-p-1,1),-1));
% set the value for omega
omega=10;
% compute I_omega
I_o=sparse(diag([1/kappa;ones(T-p-1,1)]));
GIG=G'*I_o*G;
alpha0=0.01;
alphabar=T-p+alpha0;
delta0=0.01;

% draw the parameters in turn
   for jj=1:n
   % estimate deltabar
   deltabar=logHdraw(1:end,jj)'*GIG*logHdraw(1:end,jj)+delta0;
   % draw the value phi_i
   phi_Hdraw(jj)=bear.igrandn(alphabar/2,deltabar/2);
   end

   %% draw A
   Hvars=exp([zeros(n,p) logHdraw']');

   [Adraw]=bear.sampleA_H(yData,Psi,B,Hvars,T,priorValues,dataValues);
   AdrawInv=Adraw\eye(n);

   %% sample the initial value
Kh0 = a0 + sparse(1:n,1:n,1./phi_Hdraw');
h0_hat = Kh0\(a0*b0 + log(Hvars(p+1,:)')./phi_Hdraw');
h0 = h0_hat + chol(Kh0,'lower')'\randn(n,1);

   %% Construct H
Hdraw=zeros(n,n,T);

for t=1:T
    Hdraw(:,:,t)=AdrawInv*diag(Hvars(t,:))*AdrawInv';
end

HvarsDraw=Hvars;

author: european-central-bank

repo: BEAR-toolbox


text:
function plot_fourier_approx(ai, n, m, normalized, color, line_width)

% This function will plot the fourier approximation, given a chain code (ai),
% number of harmonic elements (n), and number of points for reconstruction (m).
% Normalization can be applied by setting "normalized = 1".

    if (nargin < 5)
        color = 'b';
        line_width = 2;
    end

    if (nargin < 6)
        line_width = 2;
    end

    % Do Fourier approximatoin
    k = size(ai, 2);
    x_ = fourier_approx(ai, n, m, normalized);

    % Make it closed contour
    x = [x_; x_(1,1) x_(1,2)];

    plot(x(:,1), x(:,2), color, 'linewidth', line_width);

end

author: Sable

repo: mcbench-benchmarks


text:
% make movie of whole-brain activity, with stimulus and motor icons.

%% (manually run twice) 2 stimulus types: Phototaxis and OMR
isPTnotOMR = 0; % 1
choose_xyz = 'x';

%% load files
if isPTnotOMR
    filedir = 'C:\Users\Xiu\Downloads\chuckFish8\phototaxis\max_z';
    load('C:\Users\Xiu\Downloads\chuckFish8\phototaxis\PT_frame.mat');
    swim = PT_frame;
else
    filedir = ['C:\Users\Xiu\Downloads\chuckFish8\omr\max_',choose_xyz];
    load('C:\Users\Xiu\Downloads\chuckFish8\omr\OMR_swim_frame.mat');
    swim = OMR_swim_frame;
end

% load images
dir1 = dir(fullfile(filedir, '*.tif'));
nFrames = length(dir1);

% load first image
im1 = imread(fullfile(filedir,dir1(1).name));
[s1,s2] = size(im1);

% load all images
IM = zeros(s1,s2,nFrames);
for i_frame = 1:nFrames
    IM(:,:,i_frame) = imread(fullfile(filedir,dir1(i_frame).name));
end

%% subtract average (dFF) for this chunk of data (output in correct scale)
imavr = mean(IM,3);

IMavr = repmat(imavr,1,1,nFrames);
dIM0 = (IM-IMavr)./IMavr;

%% add fish outline
if choose_xyz == 'x'
    AI_fishoutline = imread('C:\Users\Xiu\Dropbox (Personal)\!Proj FishExplorer\presentation\Media\AI_fishoutline.tif');
    imfish = AI_fishoutline(:,:,1);
    IM_fish = double(repmat(imfish,1,1,nFrames));
    dIM = dIM0 + IM_fish;
else
    dIM = dIM0;

end
%% draw movie
if false
    %%
    figure; colormap(hot);
    clim = [0.1,1];
    for i_frame = 1:nFrames
        imagesc(dIM(:,:,i_frame),clim);axis equal; axis off
        drawnow;
        F(i_frame) = getframe;%im2frame(C);
        %    cdata = print('-RGBImage','-r300');
    end
    colorbar
    %%
    im = dIM(:,:,i_frame);
    k_zres_ratio = 20; % 19.7
    im = imresize(im, [s1*k_zres_ratio,s2],'nearest');
    figure;
    imagesc(im)
end

%% add stimulus
% manual
if isPTnotOMR % 0:white. -1:right; 1:left
    stim_1rep = horzcat(0*ones(1,20),-1*ones(1,40),0*ones(1,20),1*ones(1,40));
    stim = repmat(stim_1rep,1,4);
else % -2:baseline. 0:forward. -1:right. 1:left.
    stim_1rep = horzcat(-2*ones(1,20),0*ones(1,30),-2*ones(1,20),-1*ones(1,30),-2*ones(1,20),1*ones(1,30));
    stim = repmat(stim_1rep,1,3);
end

figure;hold on;
plot(stim,'k');
plot(swim(:,1),'r');

turns = swim(:,1);
turns = turns/max(abs(turns));
turns = smooth(turns,10);
fw = swim(:,2); % not used...

%% draw movie with stim icon
% pad with extra space for stim icon
IM2 = horzcat(zeros(568,20,nFrames),dIM);
[s1,s2,~] = size(IM2);

r1 = 20;
r2 = 40;
circlemaskIX = MakeCircularMask(r1,[s1,s2]);
center_ix = sub2ind([s1,s2],s1-r2,r2);
stimIX_left = circlemaskIX+center_ix;
center_ix = sub2ind([s1,s2],r2,r2);
stimIX_right = circlemaskIX+center_ix;
stimIX_both = union(stimIX_left,stimIX_right);

figure; colormap(hot);
clim = [0.1,1];

for i_frame = 1:nFrames
    im = IM2(:,:,i_frame);
    % add stim
    switch stim(i_frame)
        case 1 % left
            im(stimIX_left) = 0.5;
        case 0 % forward
            if ~isPTnotOMR
                im(stimIX_both) = 0.5;
            end
        case -1 % right
            im(stimIX_right) = 0.5;
    end

    % add motor
    if turns(i_frame)>0
        % (left)
        radius = round(sqrt(abs(turns(i_frame)))*r2);
        circlemaskIX = MakeCircularMask(radius,[s1,s2]);
        center_ix = sub2ind([s1,s2],s1-r2,s2-r2);
        motor_left = circlemaskIX+center_ix;
        im(motor_left) = 1;
    else
        % (right)
        radius = round(sqrt(abs(turns(i_frame)))*r2);
        circlemaskIX = MakeCircularMask(radius,[s1,s2]);
        center_ix = sub2ind([s1,s2],r2,s2-r2);
        motor_right = circlemaskIX+center_ix;
        im(motor_right) = 1;
    end

    % draw
    imagesc(im,clim);axis equal; axis off
%     drawnow;
    F(i_frame) = getframe;
end

%% show movie again, if desired
figure
nloops = 3;
fps = 20;
movie(F,nloops,fps);
axis off

%% write video to file
if isPTnotOMR
    videoname = 'PT_dFF_video_Fish8.avi';
else
    videoname = 'OMR_dFF_video_Fish8.avi';
end
myVideo = VideoWriter(videoname);
myVideo.FrameRate = 20;  % Default 30
myVideo.Quality = 50;    % Default 75
open(myVideo);
writeVideo(myVideo, F);
close(myVideo);

author: xiuyechen

repo: FishExplorer


text:
function psi = psifunc(centers,x,dis)
%This function is the psi second order spline function
% If other orders are desired, this function can be replaced
% centers are the knot locations
%x is a nxd matrix where d is the dimension and n are the number of points
%dis is an optional distance matrix. If it is not supplied, simple
% Euclidean distance will be calculated. This can also be used to save
% computational time if the matrix is calculated elsewhere

if nargin<3
    dis = pdist2(centers,x);
end
r=dis.*dis;
r(r==0)=1;
psi=r.*log(r);
end


text:
function partition=new_partitions(zones,sigma,r,N,x)
%sigma=1;
%r=8.4;
T=r/2/sigma;
kappa=size(zones);
kappa=kappa(1);
%width=sigma*(2*N/r-1);
%x=-width:1/200:width;

partition=zeros(kappa,length(x));

if kappa==1
  k=kappa;
  for j=1:length(x)
      if x(j)<=(zones(k,1)-1)/T+sigma
        partition(k,j)=0;
      elseif (zones(k,1)-1)/T+sigma<x(j) & x(j)<-sigma
        partition(k,j)=rho((-x(j)-sigma)/(-2*sigma-(zones(k,1)-1)/T));
      elseif -sigma<=x(j) & x(j)<=sigma
        partition(k,j)=1;
      elseif sigma<x(j) & x(j)<(zones(k,2)+1)/T-sigma
        partition(k,j)=rho((x(j)-sigma)/((zones(k,2)+1)/T-2*sigma));
      else
        partition(k,j)=0;
      end
  end
else
for k=1:kappa
  if k==1
    for j=1:length(x)
      if x(j)<=(zones(k,1)-1)/T+sigma
        partition(k,j)=0;
      elseif (zones(k,1)-1)/T+sigma<x(j) & x(j)<-sigma
        partition(k,j)=rho((-x(j)-sigma)/(-2*sigma-(zones(k,1)-1)/T));
      elseif -sigma<=x(j) & x(j)<=(zones(k+1,1)-1)/T+sigma
        partition(k,j)=1;
      elseif (zones(k+1,1)-1)/T+sigma<x(j) & x(j)<(zones(k,2)+1)/T-sigma
	partition(k,j)=rho((x(j)-((zones(k+1,1)-1)/T+sigma))/((zones(k,2)-zones(k+1,1)+2)/T-2*sigma));
%      elseif (1-N)/T+sigma<x(j) & x(j)<(1-N/2)/T
%        partition(k,j)=1-rho((-x(j)+(1-N/2)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(1-N/2)/T
%        partition(k,j)=1/2;
%      elseif (1-N/2)/T<x(j) & x(j)<1/T-sigma
%        partition(k,j)=rho((x(j)-(1-N/2)/T)/(N/T/2-sigma))/2;
      else
        partition(k,j)=0;
      end
    end
  elseif k==kappa
    for j=1:length(x)
      if x(j)<=(zones(k,1)-1)/T+sigma
        partition(k,j)=0;
      elseif (zones(k,1)-1)/T+sigma<x(j) & x(j)<(zones(k-1,2)+1)/T-sigma
        partition(k,j)=1-rho((x(j)-((zones(k,1)-1)/T+sigma))/((zones(k-1,2)-zones(k,1)+2)/T-2*sigma));
%      elseif sigma-1/T<x(j) & x(j)<(N/2-1)/T
%        partition(k,j)=rho((-x(j)+(N/2-1)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(N/2-1)/T
%        partition(k,j)=1/2;
%      elseif (N/2-1)/T<x(j) & x(j)<(N-1)/T-sigma
%        partition(k,j)=1-rho((x(j)-(N/2-1)/T)/(N/T/2-sigma))/2;
      elseif (zones(k-1,2)+1)/T-sigma<=x(j) & x(j)<=sigma
        partition(k,j)=1;
      elseif sigma<x(j) & x(j)<(zones(k,2)+1)/T-sigma
        partition(k,j)=rho((x(j)-sigma)/((zones(k,2)+1)/T-2*sigma));
      else
        partition(k,j)=0;
      end
    end
  else
    for j=1:length(x)
      if x(j)<=(zones(k,1)-1)/T+sigma
        partition(k,j)=0;
      elseif (zones(k,1)-1)/T+sigma<x(j) & x(j)<(zones(k-1,2)+1)/T-sigma
        partition(k,j)=1-rho((x(j)-((zones(k,1)-1)/T+sigma))/((zones(k-1,2)-zones(k,1)+2)/T-2*sigma));
%      elseif (-N+k-1)/T+sigma<x(j) & x(j)<(k-1-N/2)/T
%        partition(k,j)=rho((-x(j)+(k-1-N/2)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(k-1-N/2)/T
%        partition(k,j)=1/2;
%      elseif (k-1-N/2)/T<x(j) & x(j)<(k-1)/T-sigma
%        partition(k,j)=1-rho((x(j)-(k-1-N/2)/T)/(N/T/2-sigma))/2;
      elseif (zones(k-1,2)+1)/T-sigma<=x(j) & x(j)<=(zones(k+1,1)-1)/T+sigma
        partition(k,j)=1;
      elseif (zones(k+1,1)-1)/T+sigma<x(j) & x(j)<(zones(k,2)+1)/T-sigma
        partition(k,j)=rho((x(j)-((zones(k+1,1)-1)/T+sigma))/((zones(k,2)-zones(k+1,1)+2)/T-2*sigma));
%      elseif (-N+k)/T+sigma<x(j) & x(j)<(k-N/2)/T
%        partition(k,j)=1-rho((-x(j)+(k-N/2)/T)/(N/T/2-sigma))/2;
%      elseif x(j)==(k-N/2)/T
%        partition(k,j)=1/2;
%      elseif (k-N/2)/T<x(j) & x(j)<k/T-sigma
%        partition(k,j)=rho((x(j)-(k-N/2)/T)/(N/T/2-sigma))/2;
      else
        partition(k,j)=0;
      end
    end
  end
%  plot(x,partition(k,:))
%  axis([min(x) max(x) 0 1])
%  [k kappa]
%  zones(k,:)
%  pause
end
end
%plot(x,partition'); %'
%pause

%total=zeros(size(x));
%for j=1:N
%total=total+partition(j,:);
%end
%spy(1-total);

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% GEOMETRIC ASIAN OPTION PRICER  (RUN SCRIPT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Descritpion: Script to Price Geometric Asian options in Levy Models using the PROJ method
% Author:      Justin Kirkby
%
% Reference:   (1) An Efficient Transform Method For Asian Option Pricing, SIAM J. Financial Math., 2016
%              (2) Efficient Option Pricing by Frame Duality with the Fast Fourier Transform.
%                  SIAM J. Financial Math (2015), Kirkby, J.L
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[folder, name, ext] = fileparts(which( mfilename('fullpath')));
cd(folder);
addpath('../RN_CHF')
addpath('../Helper_Functions')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 1) CHOOSE CONTRACT/GENERAL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
call = 1;    %For call use 1 (else, its a put)
S_0  = 100;  %Initial price
W    = 100;  %Strike            %NOTE: no error handling in place for extreme values of W (increase grid if strike falls outside)
r    = 0.05;  %Interest rate
q    = 0.00;  %dividend yield
T    = 1;    %Time (in years)
M    = 52;   %Number of monitoring points

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 2) CHOOSE MODEL PARAMETERS (Levy Models)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model = 1;   %See Models Below (e.g. model 1 is Black Scholes), and choose specific params
params = {};

if model == 1 %BSM (Black Scholes Merton)
    params.sigmaBSM = 0.15;    %CHOOSE

elseif model == 2 %CGMY
    params.C  = 0.02;
    params.G  = 5;
    params.MM = 15;
    params.Y  = 1.2;

elseif model == 3 %NIG
    params.alpha = 15;
    params.beta  = -5;
    params.delta = 0.5;

elseif model == 4 %MJD (Merton Jump Diffusion)
    params.sigma  = 0.12;
    params.lam    = 0.4;
    params.muj    = -0.12;
    params.sigmaj = 0.18;

elseif model == 5 %Kou Double Expo
    params.sigma = 0.15;
    params.lam   = 3;
    params.p_up  = 0.2;
    params.eta1  = 25;
    params.eta2  = 10;

elseif model == 8 % Variance Gamma
    params.sigma = 0.2;
    params.nu = 0.85;
    params.theta = 0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 3) CHOOSE PROJ PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UseCumulant = 1;  %Set to 1 to use the cumulant base rule (Approach 1) to determine gridwidth, else used fixed witdth (Approach 2)

%---------------------
% APPROACH 1: Cumulant Based approach for grid width
% (see "Robust Option Pricing with Characteritics Functions and the BSpline Order of Density Projection")
%---------------------
if UseCumulant ==1  %With cumulant based rule, choose N and Alpha (N = 2^(P+Pbar) based on second approach)
    logN  = 8;   %Uses N = 2^logN  gridpoint
    L1 = 10;  % determines grid witdth (usually set L1 = 8 to 15 for Levy)
%---------------------
% APPROACH 2: Manual GridWidth approach
%---------------------
else %Manually specify resolution and Pbar
    P     = 6;  % resolution is 2^P
    Pbar  = 3;  % Determines density truncation grid with, 2^Pbar
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PRICE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Note: rnCHF is the risk netural CHF, c1,c2,c4 are the cumulants
dt = T/M;
modelInput = getModelInput(model, dt, r, q, params);

if UseCumulant ==1  % Choose density truncation width based on cumulants
    alpha = getTruncationAlpha(T, L1, modelInput, model);
else    % Manually supply density truncation width above
    logN = P + Pbar;
    alpha = 2^Pbar/2;
end
N = 2^logN;    % grid roughly centered on [c1 - alph, c1 + alph]

tic
price = PROJ_Geometric_Asian(N, alpha, S_0, M, W, call, T, r, q, modelInput.rnSYMB);
toc

fprintf('Geometric Price: %.8f \n', price)

compare_arithmetic = 1;

if compare_arithmetic
    addpath('../Asian_Options')
    price_arith = PROJ_Asian( N,alpha,S_0,M,W,call,T,r,q, modelInput.rnCHF, modelInput.RNmu*dt);
    fprintf('Arithmetic Price: %.8f \n', price_arith)
end


text:
% Locate place fields in a firing map.
%
% Identifies the place fields in 2D firing map. First map is converted to binary image
% by applying a global threshold. During the second step connected regions are labeled
% by function bwlabel. Final step includes filtering of identified regions and verification
% that they are valid place fields.
%
%  USAGE
%   [fieldsMap, fields] = analyses.palcefield(map, <options>)
%   map           Firing rate map either structure obtained using analyses.map
%                 or a matrix that represents firing map.
%   <options>     optional list of property-value pairs (see table below)
%
%   =========================================================================
%    Properties     Values
%   -------------------------------------------------------------------------
%    'threshold'    value above threshold*peak belong to a field (default = 0.2).
%    'minBins'      Minimum number of bins in a place field, i.e. total number
%                   of bins of place field or an area of place field. Fields with
%                   fewer bins are not considered as place fields. Remember to
%                   adjust this value when you change the bin width.
%                   (default = 9).
%    'minPeak'      peaks smaller than this value are considered spurious
%                   and ignored (default = 1). Peak is normally a rate, however
%                   it's units not necessary are Hz.
%    'binWidth'     Width of the bins in cm. It is used to calculate field size.
%                   (default = 1).
%    'pos'          Position samples. Used to calculate posInd. If not provided,
%                   then posInd will be an empty matrix.
%    'minMeanRate'  Fields with mean rate smaller that this value are ignored (default = 0).
%   =========================================================================
%
%  OUTPUT
%
%   fieldsMap       Matrix of the same size as firing map. The elements of fieldsMap
%                   are integer values greater or equal to 0. The elements labeled 0
%                   are the background (not a place field). The pixels labeled 1 make up
%                   first field; the pixels labeled 2 make up a second object; and so on.
%
%   fields          Structure with information about each field. Structure fields are:
%       row         vector of rows that constitute this field;
%       col         vector of columns that constitute this field;
%       area        area of field measured by regionprops function, i.e. number of bins in the field;
%       bbox        bounding box of the field;
%       peak        field peak value;
%       size        field size in cm; Calculated with help of binWidth argument.
%       x, y        field centre of mass point;
%       meanRate    mean firing rate;
%       PixelIdxList    linear list of indices that constitute the field.
%       peakX       X-coordinate of field's peak
%       peakY       Y-coordinate of field's peak
%       map         Binary map of the field. It has the same size as firing map. The elements of 'map'
%                   equal to one if they belong to the field and are zeros otherwise.
%       posInd      Indices of position samples that correspond to this field. In case position sample
%                   matrix is of size Nx5 ([t x y x1 y1]), posInd corresponds to the left most position
%                   columns ([x y]). If pos argument is not provided, then posInd will be an empty matrix.
%
%  EXAMPLES
%
%   1. To obtain positions of the first field run:
%       [~, fields] = analyses.placefield(map, 'threshold', 0.3, 'minBins', 5, 'minPeak', 0.1, 'pos', pos);
%       fieldPos = pos(fields(1).posInd, :);
%

% Note to developers:
% In order to improve speed and capture all fields the map is processed in two passes. First
% we extract fields based solely on global peak value. This should be enough for most of the cases.
% Secondly we check if extracted fields maxima is the same as results of imregionalmax. If they are
% different, then we perform another extraction, which is slower.
%
function [fieldsMap, fields] = placefield(map, varargin)
    inp = inputParser;
    defaultMinPeak = 1.;
    defaultThreshold = 0.2;
    defaultMinBins = 9;
    defaultBinWidth = 1;
    defaultPos = [];
    defaultMinMeanRate = 0;

    % input argument check functions
    checkThreshold = @(x) helpers.isdscalar(x, '>=0', '<=1');
    checkScalarZero = @(x) helpers.isiscalar(x, '>=0');
    checkDScalar = @(x) helpers.isdscalar(x, '>0');
    checkDScalarZero = @(x) helpers.isdscalar(x, '>=0');

    % fill input parser object
    addRequired(inp, 'map');
    addParameter(inp, 'threshold', defaultThreshold, checkThreshold);
    addParameter(inp, 'minBins', defaultMinBins, checkScalarZero);
    addParameter(inp, 'minPeak', defaultMinPeak, checkDScalarZero);
    addParameter(inp, 'binWidth', defaultBinWidth, checkDScalar);
    addParameter(inp, 'minMeanRate', defaultMinMeanRate, checkScalarZero);
    addParameter(inp, 'pos', defaultPos, @(x) ismatrix(x) && size(x, 2) >= 3);

    parse(inp, map, varargin{:});

    % get parsed arguments
    minPeak = inp.Results.minPeak;
    threshold = inp.Results.threshold;
    minBins = inp.Results.minBins;
    binWidth = inp.Results.binWidth;
    pos = inp.Results.pos;
    minMeanRate = inp.Results.minMeanRate;

    originalMap = [];
    if isstruct(map)
        originalMap = map;
        map = map.z;
    end

    fieldsMap = zeros(size(map));
    fields = struct('row', {}, 'col', {}, ...
        'size', {}, 'peak', {}, 'peakX', {}, 'peakY', {}, ...
        'area', {}, 'bbox', {}, ...
        'x', {}, 'y', {}, ...
        'meanRate', {}, 'PixelIdxList', {}, ...
        'map', {}, 'posInd', {} ...
        );

    if isempty(map)
        return
    end

    globalPeak = nanmax(nanmax(map));
    if isnan(globalPeak) || globalPeak == 0
        return;
    end

    mapNans = isnan(map);
    map(mapNans) = 0;
    regionalMaxMap = imregionalmax(map, 4); % obtain all local maxima
    testRegionalMx = zeros(size(map)); % map that will contain located fields maxima

    [ir, ic] = find(regionalMaxMap > 0);            % get peak locations
    foundPeaks = map(sub2ind(size(map), ir ,ic));   % obtain peaks value

    % remove peaks that have smaller rate than minPeak
    selected = foundPeaks < minPeak;
    regionalMaxMap(ir(selected), ic(selected)) = 0;

    % Counter for the number of fields
    nFields = 0;

    binMap = ones(size(map));

    binMap(map < globalPeak*threshold) = 0;
    binMap(mapNans) = 0;

    cc = bwconncomp(binMap, 4);
    stats = regionprops(cc, 'Area', 'BoundingBox', 'Centroid');

    for i = 1:cc.NumObjects
        linInd = cc.PixelIdxList{i};
        [fieldPeak, peakLinearInd] = max(map(linInd));
        [r, c] = ind2sub(size(map), linInd);
        meanRate = nanmean(map(linInd));

        [pr, pc] = ind2sub(size(map), linInd(peakLinearInd));
        % mark this field as visited even if we reject if further down
        testRegionalMx(pr, pc) = 1;

        if fieldPeak < minPeak
            continue;
        end
        if meanRate < minMeanRate
            continue;
        end

        if length(r) >= minBins
            nFields = nFields + 1;

            fields(nFields).row = r;
            fields(nFields).col = c;
            fields(nFields).size = length(r) * binWidth^2;
            fields(nFields).peak = fieldPeak;
            fields(nFields).peakX = pc;
            fields(nFields).peakY = pr;
            fields(nFields).area = stats(i).Area;
            fields(nFields).bbox = stats(i).BoundingBox;
            fields(nFields).PixelIdxList = linInd;

            fields(nFields).x = stats(i).Centroid(1);
            fields(nFields).y = stats(i).Centroid(2);

            fields(nFields).meanRate = meanRate;
            fields(nFields).map = zeros(size(map));
            fields(nFields).map(linInd) = 1;
            fields(nFields).map(mapNans) = nan;

            % fields(nFields).Extent = stats(i).Extent;
            % fields(nFields).Orientation = stats(i).Orientation;
            % fields(nFields).Eccentricity = stats(i).Eccentricity;
            % fields(nFields).MinorAxisLength = stats(i).MinorAxisLength;
            % fields(nFields).MajorAxisLength = stats(i).MajorAxisLength;
            % fields(nFields).Perimeter = stats(i).Perimeter;

            if ~isempty(pos)
                intRect = ceil(fields(nFields).bbox); % bounding box with integers
                if isempty(originalMap)
                    % we do not have position distribution from the map, so assume that it's based
                    % on data min/max values
                    nBins = size(map);
                    limitsX = [nanmin(pos(:, bntConstants.PosX)) nanmax(pos(:, bntConstants.PosX))];
                    limitsY = [nanmin(pos(:, bntConstants.PosY)) nanmax(pos(:, bntConstants.PosY))];
                    xSpace = linspace(limitsX(1), limitsX(2), nBins(2));
                    ySpace = linspace(limitsY(1), limitsY(2), nBins(1));
                else
                    xSpace = originalMap.x;
                    ySpace = originalMap.y;
                end

                xMin = xSpace(intRect(1));
                xMax = xSpace(intRect(1) + intRect(3) - 1);

                yMin = ySpace(intRect(2));
                yMax = ySpace(intRect(2) + intRect(4) - 1);

                posIndX = pos(:, bntConstants.PosX) >= xMin & pos(:, bntConstants.PosX) <= xMax;
                posIndY = pos(:, bntConstants.PosY) >= yMin & pos(:, bntConstants.PosY) <= yMax;
                fields(nFields).posInd = find(posIndX & posIndY);
            else
                fields(nFields).posInd = [];
            end

            fieldsMap(linInd) = nFields;
        end
    end

    if ~isequaln(testRegionalMx, regionalMaxMap)
        map(fieldsMap > 0) = 0; % turn off map values for known fields, prevent field duplicates

        % we have some uncounted fields
        leftPeaksMap = regionalMaxMap - testRegionalMx;
        [ir, ic] = find(leftPeaksMap > 0);                    % get peak locations
        foundPeaks = map(sub2ind(size(leftPeaksMap), ir ,ic));   % obtain peaks value
        mapThresholds = foundPeaks * threshold;
        if sum(foundPeaks) == 0
            return;
        end

        finalMap = zeros(size(map));
        for i = 1:length(foundPeaks)
            binMap = zeros(size(map));
            binMap( (map > mapThresholds(i)) & (map <= foundPeaks(i)) ) = 1;
            binMap(mapNans) = 0;

            [binMap, linInd] = bwselect(binMap, ic(i), ir(i), 4); % leave only region that relates to current field
            if isempty(linInd)
                continue;
            end

            % check for minimum number of bins
            [r, ~] = ind2sub(size(map), linInd);
            if length(r) < minBins
                continue;
            end

            stats = regionprops(binMap, 'Centroid', 'EulerNumber'); % find statistics on field candidates
            distToPeak = sqrt((stats.Centroid(1) - ic(i))^2 + (stats.Centroid(2) - ir(i))^2);
    %         distToPeak = pdist2(stats.Centroid, [ic(i) ir(i)]);
            if distToPeak > 4 || stats.EulerNumber < 1 % we want object without holes (euler number)
                continue;
            end
            finalMap(linInd) = 1;
        end

        cc = bwconncomp(finalMap, 4); % somehow regionprops is not working if finalMap is passed directly
        stats = regionprops(cc, 'Area', 'BoundingBox', 'Centroid');

        for fieldInd = 1:length(stats)
            linInd = cc.PixelIdxList{fieldInd};
            [r, c] = ind2sub(size(map), linInd);
            meanRate = nanmean(map(linInd));
            [peakRate, peakInd] = nanmax(map(linInd));
            [pr, pc] = ind2sub(size(map), linInd(peakInd));

            if peakRate < minPeak
                continue;
            end
            if meanRate < minMeanRate
                continue;
            end

            nFields = nFields + 1;

            fields(nFields).row = r;
            fields(nFields).col = c;
            fields(nFields).size = length(r) * binWidth^2;
            fields(nFields).peak = peakRate;
            fields(nFields).peakX = pc;
            fields(nFields).peakY = pr;
            fields(nFields).area = stats(fieldInd).Area;
            fields(nFields).bbox = stats(fieldInd).BoundingBox;
            fields(nFields).PixelIdxList = linInd;

            fields(nFields).x = stats(fieldInd).Centroid(1);
            fields(nFields).y = stats(fieldInd).Centroid(2);

            fields(nFields).meanRate = meanRate;
            fields(nFields).map = zeros(size(map));
            fields(nFields).map(linInd) = 1;
            fields(nFields).map(mapNans) = nan;

            if ~isempty(pos)
                intRect = ceil(fields(nFields).bbox); % bounding box with integers
                if isempty(originalMap)
                    % we do not have position distribution from the map, so assume that it's based
                    % on data min/max values
                    nBins = size(map);
                    limitsX = [nanmin(pos(:, bntConstants.PosX)) nanmax(pos(:, bntConstants.PosX))];
                    limitsY = [nanmin(pos(:, bntConstants.PosY)) nanmax(pos(:, bntConstants.PosY))];
                    xSpace = linspace(limitsX(1), limitsX(2), nBins(2));
                    ySpace = linspace(limitsY(1), limitsY(2), nBins(1));
                else
                    xSpace = originalMap.x;
                    ySpace = originalMap.y;
                end

                xMin = xSpace(intRect(1));
                xMax = xSpace(intRect(1) + intRect(3) - 1);

                yMin = ySpace(intRect(2));
                yMax = ySpace(intRect(2) + intRect(4) - 1);

                posIndX = pos(:, bntConstants.PosX) >= xMin & pos(:, bntConstants.PosX) <= xMax;
                posIndY = pos(:, bntConstants.PosY) >= yMin & pos(:, bntConstants.PosY) <= yMax;
                fields(nFields).posInd = find(posIndX & posIndY);
            else
                fields(nFields).posInd = [];
            end

            fieldsMap(linInd) = nFields;
        end
    end
end

author: kavli-ntnu


text:
function sconfounds = spm_cfg_eeg_spatial_confounds
% configuration file for reading montage files
%__________________________________________________________________________
% Copyright (C) 2014-2016 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_spatial_confounds.m 6926 2016-11-09 22:13:19Z guillaume $

D = cfg_files;
D.tag = 'D';
D.name = 'File Name';
D.filter = 'mat';
D.num = [1 1];
D.help = {'Select the EEG mat file.'};

timewin         = cfg_entry;
timewin.tag     = 'timewin';
timewin.name    = 'Time window';
timewin.strtype = 'r';
timewin.num     = [1 2];
timewin.val     = {[-Inf Inf]};
timewin.help    = {'Time window (ms)'};

condlabel = cfg_entry;
condlabel.tag = 'conditions';
condlabel.name = 'Condition label';
condlabel.strtype = 's';
condlabel.help = {''};

conditions = cfg_repeat;
conditions.tag = 'condrepeat';
conditions.name = 'Conditions';
conditions.help = {'Specify the labels of the conditions to include in the SVD.'};
conditions.num  = [0 Inf];
conditions.values  = {condlabel};
conditions.val = {};

ncomp = cfg_entry;
ncomp.tag = 'ncomp';
ncomp.name = 'Number of components';
ncomp.strtype = 'n';
ncomp.num = [1 1];
ncomp.help = {'Number of confound components to keep.'};

threshold = cfg_entry;
threshold.tag = 'threshold';
threshold.name = 'Threshold';
threshold.strtype = 'r';
threshold.val = {NaN};
threshold.num = [1 1];
threshold.help = {'Threshold for data amplitude after correction.',...
    'If defined this overrides number of components.'};

svd = cfg_branch;
svd.tag = 'svd';
svd.name = 'SVD';
svd.val = {timewin, conditions, threshold, ncomp};
svd.help = {'Define confounds from SVD of artefact samples.'};

conffile = cfg_files;
conffile.tag = 'conffile';
conffile.name = 'File name';
conffile.filter = 'mat';
conffile.num = [1 1];
conffile.help = {'Select the EEG mat file.'};

spmeeg = cfg_branch;
spmeeg.tag = 'spmeeg';
spmeeg.name = 'SPM M/EEG dataset';
spmeeg.val = {conffile};
spmeeg.help = {'Load confounds defined in a different dataset.'};

conffile = cfg_files;
conffile.tag = 'conffile';
conffile.name = 'BESA confounds file';
conffile.filter = '.*.bsa';
conffile.num = [1 1];
conffile.help = {'Select the EEG mat file.'};

besa = cfg_branch;
besa.tag = 'besa';
besa.name = 'BESA';
besa.val = {conffile};
besa.help = {'Load confounds defined in BESA from a .bsa file.'};

eyes = cfg_const;
eyes.tag = 'eyes';
eyes.name = 'Eyes';
eyes.val  = {1};
eyes.help = {'Defined confounds based on standard locations for the eyes',...
    'This is a very crude method, a data-driven method is preferable'};

clr = cfg_const;
clr.tag = 'clear';
clr.name = 'Clear';
clr.val  = {1};
clr.help = {'Clear previously defined spatial confounds.'};

mode = cfg_repeat;
mode.tag = 'mode';
mode.name = 'Mode';
mode.values = {svd, spmeeg, besa, eyes, clr};
mode.help = {'Select methods for defining spatial confounds.'};

sconfounds = cfg_exbranch;
sconfounds.tag = 'sconfounds';
sconfounds.name = 'Define spatial confounds';
sconfounds.val = {D, mode};
sconfounds.help = {'Define spatial confounds for topography-based correction of artefacts.'};
sconfounds.prog = @eeg_sconfounds;
sconfounds.vout = @vout_eeg_sconfounds;
sconfounds.modality = {'EEG'};

function out = eeg_sconfounds(job)

D = char(job.D{1});

for i = 1:numel(job.mode)
    mode = job.mode{i};
    S = [];
    S.D = D;
    S.mode = char(fieldnames(mode));
    switch S.mode
        case 'svd'
            if ~isnan(mode.svd.threshold)
                S.threshold = mode.svd.threshold;
            else
                S.ncomp = mode.svd.ncomp;
            end
            S.timewin = mode.svd.timewin;
        case {'besa', 'spmeeg'}
            S.conffile = char(mode.(S.mode).conffile);
        otherwise
            % do nothing
    end

    D = spm_eeg_spatial_confounds(S);
end

out.D = {fullfile(D)};

function dep = vout_eeg_sconfounds(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Dataset with spatial confounds';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function sparse_grid_hw_test ( )

%*****************************************************************************80
%
%% SPARSE_GRID_HW_TEST tests the SPARSE_GRID_HW library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_HW_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SPARSE_GRID_HW library.\n' );

  cce_test ( );
  cce_sparse_test ( );

  ccl_test ( );
  ccl_sparse_test ( );

  ccs_test ( );
  ccs_sparse_test ( );

  get_seq_test ( );

  glo_test ( );
  glo_sparse_test ( );

  gqn_test ( );
  gqn_sparse_test ( );

  gqu_test ( );
  gqu_sparse_test ( );

  kpn_test ( );
  kpn_sparse_test ( );

  kpu_test ( );
  kpu_sparse_test ( );

  nwspgr_size_test ( );
  nwspgr_time_test ( );
  nwspgr_test ( );

  order_report ( );

  symmetric_sparse_size_test ( );

  tensor_product_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_HW_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [xf, S, cnt] = LMFsolveOLD(varargin)
% Solve a Set of Overdetermined Nonlinear Equations in Least-Squares Sense.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A solution is obtained by a Fletcher version of the Levenberg-Maquardt
% algoritm for minimization of a sum of squares of equation residuals.
%
% [Xf, Ssq, CNT] = LMFsolveOLD(FUN,Xo,Options)
% FUN     is a function handle or a function M-file name that evaluates
%         m-vector of equation residuals,
% Xo      is n-vector of initial guesses of solution,
% Options is an optional set of Name/Value pairs of control parameters
%         of the algorithm. It may be also preset by calling:
%         Options = LMFsolveOLD('default'), or by a set of Name/Value pairs:
%         Options = LMFsolveOLD('Name',Value, ... ), or updating the Options
%                   set by calling
%         Options = LMFsolveOLD(Options,'Name',Value, ...).
%
%    Name   Values {default}         Description
% 'Display'     integer     Display iteration information
%                            {0}  no display
%                             k   display initial and every k-th iteration;
% 'Jacobian'    handle      Jacobian matrix function handle; {@finjac}
% 'FunTol'      {1e-7}      norm(FUN(x),1) stopping tolerance;
% 'XTol'        {1e-7}      norm(x-xold,1) stopping tolerance;
% 'MaxIter'     {100}       Maximum number of iterations;
% 'ScaleD'                  Scale control:
%               value        D = eye(m)*value;
%               vector       D = diag(vector);
%                {[]}        D(k,k) = JJ(k,k) for JJ(k,k)>0, or
%                                   = 1 otherwise,
%                                     where JJ = J.'*J
% Not defined fields of the Options structure are filled by default values.
%
% Output Arguments:
% Xf        final solution approximation
% Ssq       sum of squares of residuals
% Cnt       >0          count of iterations
%           -MaxIter,   did not converge in MaxIter iterations

% Example:
% The general Rosenbrock's function has the form
%    f(x) = 100(x(1)-x(2)^2)^2 + (1-x(1))^2
% Optimum solution gives f(x)=0 for x(1)=x(2)=1. Function f(x) can be
% expressed in the form
%    f(x) = f1(x)^2 =f2(x)^2,
% where f1(x) = 10(x(1)-x(2)^2) and f2(x) = 1-x(1).
% Values of the functions f1(x) and f2(x) can be used as residuals.
% LMFsolveOLD finds the solution of this problem in 5 iterations. The more
% complicated problem sounds:
% Find the least squares solution of the Rosenbrock valey inside a circle
% of the unit diameter centered at the origin. It is necessary to build
% third function, which is zero inside the circle and increasing outside it.
% This property has, say, the next function:
%    f3(x) = sqrt(x(1)^2 + x(2)^2) - r, where r is a radius of the circle.
% Its implementation using anonymous functions has the form
%    R  = @(x) sqrt(x'*x)-.5;    %   A distance from the radius r=0.5
%    ros= @(x) [10*(x(2)-x(1)^2); 1-x(1); (R(x)>0)*R(x)*1000];
%    [x,ssq,cnt]=LMFsolveOLD(ros,[-1.2,1],'Display',1,'MaxIter',50)
% Solution: x = [0.4556; 0.2059],  |x| = 0.5000
% sum of squares: ssq = 0.2966,
% number of iterations: cnt = 18.
%
% Note:
% Users with older MATLAB versions, which have no anonymous functions
% implemented, have to call LMFsolveOLD with named function for residuals.
% For above example it is
%
%   [x,ssq,cnt]=LMFsolveOLD('rosen',[-1.2,1]);
%
% where the function rosen.m is of the form
%
%   function r = rosen(x)
%%  Rosenbrock valey with a constraint
%   R = sqrt(x(1)^2+x(2)^2)-.5;
%%  Residuals:
%   r = [ 10*(x(2)-x(1)^2)  %   first part
%         1-x(1)            %   second part
%         (R>0)*R*1000.     %   penalty
%       ];

% Reference:
% Fletcher, R., (1971): A Modified Marquardt Subroutine for Nonlinear Least
% Squares. Rpt. AERE-R 6799, Harwell

% M. Balda,
% Institute of Thermomechanics,
% Academy of Sciences of The Czech Republic,
% balda AT cdm DOT cas DOT cz
% 2007-07-02
% 2007-10-08    formal changes, improved description
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   OPTIONS
    %%%%%%%
%               Default Options
if nargin==1 && strcmpi('default',varargin(1))
   xf.Display  = 0;         %   no print of iterations
   xf.Jacobian = @finjac;   %   finite difference Jacobian approximation
   xf.MaxIter  = 100;       %   maximum number of iterations allowed
   xf.ScaleD   = [];        %   automatic scaling by D = diag(diag(J'*J))
   xf.FunTol   = 1e-7;      %   tolerace for final function value
   xf.XTol     = 1e-4;      %   tolerance on difference of x-solutions
   return

%               Updating Options
elseif isstruct(varargin{1}) % Options=LMFsolveOLD(Options,'Name','Value',...)
    if ~isfield(varargin{1},'Jacobian')
        error('Options Structure not Correct for LMFsolveOLD.')
    end
    xf=varargin{1};          %   Options
    for i=2:2:nargin-1
        name=varargin{i};     %   option to be updated
        if ~ischar(name)
            error('Parameter Names Must be Strings.')
        end
        name=lower(name(isletter(name)));
        value=varargin{i+1};  %   value of the option
        if strncmp(name,'d',1), xf.Display  = value;
        elseif strncmp(name,'f',1), xf.FunTol   = value(1);
        elseif strncmp(name,'x',1), xf.XTol     = value(1);
        elseif strncmp(name,'j',1), xf.Jacobian = value;
        elseif strncmp(name,'m',1), xf.MaxIter  = value(1);
        elseif strncmp(name,'s',1), xf.ScaleD   = value;
        else   disp(['Unknown Parameter Name --> ' name])
        end
    end
    return

%               Pairs of Options
elseif ischar(varargin{1})  % check for Options=LMFsolveOLD('Name',Value,...)
   Pnames=char('display','funtol','xtol','jacobian','maxiter','scaled');
   if strncmpi(varargin{1},Pnames,length(varargin{1}))
      xf=LMFsolveOLD('default');  % get default values
      xf=LMFsolveOLD(xf,varargin{:});
      return
   end
end

%   LMFSOLVEOLD(FUN,Xo,Options)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%

FUN=varargin{1};            %   function handle
if ~(isvarname(FUN) || isa(FUN,'function_handle'))
   error('FUN Must be a Function Handle or M-file Name.')
end

xc=varargin{2};             %   Xo

if nargin>2                 %   OPTIONS
    if isstruct(varargin{3})
        options=varargin{3};
    else
        if ~exist('options','var')
            options = LMFsolveOLD('default');
        end
        for i=3:2:size(varargin,2)-1
            options=LMFsolveOLD(options, varargin{i},varargin{i+1});
        end
    end
else
    if ~exist('options','var')
        options = LMFsolveOLD('default');
    end
end

x=xc(:);
lx=length(x);

r=feval(FUN,x);             % Residuals at starting point
%~~~~~~~~~~~~~~
S=r'*r;
epsx=options.XTol(:);
epsf=options.FunTol(:);
if length(epsx)<lx, epsx=epsx*ones(lx,1); end
J=options.Jacobian(FUN,r,x,epsx);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
A=J.'*J;                    % System matrix
v=J.'*r;

D = options.ScaleD;
if isempty(D)
    D=diag(diag(A));        % automatic scaling
    for i=1:lx
        if D(i,i)==0, D(i,i)=1; end
    end
else
    if numel(D)>1
        D=diag(sqrt(abs(D(1:lx)))); % vector of individual scaling
    else
        D=sqrt(abs(D))*eye(lx);     % scalar of unique scaling
    end
end

Rlo=0.25; Rhi=0.75;
l=1;      lc=.75;     is=0;
cnt=0;
ipr=options.Display;
printit(-1);                %   Table header
d=options.XTol;             %   vector for the first cycle
maxit = options.MaxIter;    %   maximum permitted number of iterations

while cnt<maxit && ...      %   MAIN ITERATION CYCLE
    any(abs(d)>=epsx) && ...    %%%%%%%%%%%%%%%%%%%%
    any(abs(r)>=epsf)
    d=(A+l*D)\v;            % negative solution increment
    xd=x-d;
    rd=feval(FUN,xd);
%   ~~~~~~~~~~~~~~~~~
    Sd=rd.'*rd;
    dS=d.'*(2*v-A*d);       % predicted reduction
    R=(S-Sd)/dS;

    if R>Rhi
        l=l/2;
        if l<lc, l=0; end
    elseif R<Rlo
        nu=(Sd-S)/(d.'*v)+2;
        if nu<2
            nu=2;
        elseif nu>10
            nu=10;
        end
        if l==0
            lc=1/max(abs(diag(inv(A))));
            l=lc;
            nu=nu/2;
        end
        l=nu*l;
    end
    cnt=cnt+1;
    if ipr>0 && (rem(cnt,ipr)==0 || cnt==1)
        printit(cnt,[S,l,R,x(:).'])
        printit(0,[lc,d(:).'])
    end
    if Sd>S && is==0
        is=1;
        St=S; xt=x; rt=r; Jt=J; At=A; vt=v;
    end
    if Sd<S || is>0
        S=Sd; x=xd; r=rd;
        J=options.Jacobian(FUN,r,x,epsx);
%       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        A=J.'*J;   v=J.'*r;
    else
        S=St; x=xt; r=rt; J=Jt; A=At; v=vt;
    end
    if Sd<S, is=0; end
end
xf = x;                         %   finat solution
if cnt==maxit, cnt=-cnt; end    %   maxit reached
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   PRINTIT     Printing of intermediate results
%   %%%%%%%
function printit(cnt,mx)
%        ~~~~~~~~~~~~~~~
%  cnt	=  -1  print out the header
%           0  print out second row of results
%          >0  print out first row of results

if ipr>0
   if cnt<0                 %   table header
      disp('')
      disp(char('*'*ones(1,100)))
      disp(['  cnt   SUM(r^2)        l            R' blanks(21) 'x(i) ...'])
      disp([blanks(24) 'lc' blanks(32) 'dx(i) ...'])
      disp(char('*'*ones(1,100)))
      disp('')
   else                     %   iteration output
      if cnt>0 || rem(cnt,ipr)==0
          f='%12.4e ';
          form=[f f f f '\n' blanks(49)];
          if cnt>0
              fprintf(['%4.0f ' f f f ' x = '],[cnt,mx(1:3)])
              fprintf(form,mx(4:length(mx)))
          else
              fprintf([blanks(18) f blanks(13) 'dx = '], mx(1))
              fprintf(form,mx(2:length(mx)))
          end
          fprintf('\n')
      end
   end
end
end

%   FINJAC       numerical approximation to Jacobi matrix
%   %%%%%%
function J = finjac(FUN,r,x,epsx)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lx=length(x);
J=zeros(length(r),lx);
for k=1:lx
    dx=.25*epsx(k);
    xd=x;
    xd(k)=xd(k)+dx;
    rd=feval(FUN,xd);
%   ~~~~~~~~~~~~~~~~
    J(:,k)=((rd-r)/dx);
end
end
end


text:
function[z]=celladd(x,y)
%CELLADD  Addition acting on each element in a cell array.
%
%   Z=CELLADD(X,Y) where X and Y are both cell arrays of N arrays, with
%   corresponding elements in X and in Y having the same size, returns the
%   cell array Z containing their sums,
%
%       Z{1}=X{1}+Y{1}, Z{2}=X{2}+Y{2},..., Z{N}=X{N}+Y{N}.
%
%   One of X or Y may also be a scalar or a numeric array of the same
%   length as the other input argument.
%
%   Usage: z=celladd(x,y);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008--2019 J.M. Lilly --- type 'help jlab_license' for details

if ~iscell(x)&&iscell(y)
    if length(x)==1
        x=x.*ones(length(y),1);
    else
        x=x(:);
    end
elseif iscell(x)&&~iscell(y)
    if length(y)==1
        y=y.*ones(length(x),1);
    else
        y=y(:);
    end
end

if iscell(x)
    z=x;
else
    z=y;
end

for i=1:length(y)
     if iscell(x)&&iscell(y)
         z{i}=x{i}+y{i};
     elseif ~iscell(x)&&iscell(y)
         z{i}=x(i)+y{i};
     elseif iscell(x)&&~iscell(y)
         z{i}=x{i}+y(i);
     end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [A, b, Aeq, beq, lb, ub] = snd_l_constraints(k, net_load_data, varargin)

    u0_ref = cell2mat(varargin(2));
    f = 1.05;
    g = 0.95;
    p = 1.2;
    A   = []; %A*u<=b
    b   = [];
    Aeq = [1 1 1]; %A*u = b
    beq = net_load_data(k);

    if k ~= 12
        if u0_ref(1,k) * u0_ref(1,k+1) <= 0
            if u0_ref(1,k)>0
                ub = [u0_ref(1,k)*1.05];
                lb = [u0_ref(1,k+1)*1.05];
            else
                ub = [u0_ref(1,k+1)*1.05];
                lb = [u0_ref(1,k)*1.05];
            end
        elseif u0_ref(1,k)>0 && u0_ref(1,k+1)>0
            if u0_ref(1,k)>u0_ref(1,k+1)
                ub = [u0_ref(1,k)*1.05];
                lb = [u0_ref(1,k+1)*0.95];
            else
                ub = [u0_ref(1,k+1)*1.05];
                lb = [u0_ref(1,k)*0.95];
            end
        elseif u0_ref(1,k)<0 && u0_ref(1,k+1)<0
            if u0_ref(1,k)>u0_ref(1,k+1)
                ub = [u0_ref(1,k)*0.95];
                lb = [u0_ref(1,k+1)*1.05];
            else
                ub = [u0_ref(1,k+1)*0.95];
                lb = [u0_ref(1,k)*1.05];
            end
        end
    else
        if u0_ref(1,k) >=0
            ub = [u0_ref(1,k)*1.05];
            lb = [u0_ref(1,k)*0.95];
        else
            ub = [u0_ref(1,k)*0.95];
            lb = [u0_ref(1,k)*1.05];
        end

    end

    if u0_ref(2,k)>0 %
        ub  = [ub, u0_ref(2,k)*1.1,  10];
        lb  = [lb, 0, -10];
    else
        ub  = [ub, 0,  10];
        lb  = [lb, u0_ref(2,k)*1.1, -10];
    end

end

author: juchengquan


text:
function [expr morder] = sim_ex_Schelter_2009_3_1(varargin)
% Simulation:  Schelter 2009 Eq 3.1
%
% Description:
%
% 5-variate VAR[3] system of coupled oscillators.
% This system was first described in [1].
%
% The directed graph for this model is:
% x1 -> x3
% x1 -> x4
% x2 -> x1
% x2 -> x3
% x4 -> x5
% x5 -> x4
%
% The dependency structure of this model can
% be viewed by executing the following command:
%
% >>hlp_viewGraphicsResource('sim/Schelter_2009_3_1.jpg');
%
% Author Credits:
%
% Tim Mullen, 2011
%
% References and Code:
%
% [1] (Ex 3.1, Eq. 11-15) Schelter B, Timmer J, Eichler M (2009) Assessing the strength of directed influences among neural signals using renormalized partial directed coherence. Journal of neuroscience methods 179:121-30
%
% ------------------------------------------------------------------------

% specify the default system of equations
expr_def = { ...
    'x1(t) = 0.9*x1(t-1)  + 0.3*x2(t-2)  + e1(t)' ...
    'x2(t) = 1.3*x2(t-1)  + -0.8*x2(t-2) + e2(t)' ...
    'x3(t) = 0.3*x1(t-2)  + 0.6*x2(t-1)  + e3(t)' ...
    'x4(t) = -0.7*x4(t-3) + -0.7*x1(t-3) + 0.3*x5(t-3) + e4(t)' ...
    'x5(t) = 1*x5(t-1)    + -0.4*x5(t-2) + 0.3*x4(t-2) + e5(t)' ...
    };

% set up argument definitions
arg_define(varargin, ...
    arg({'expr','DynamicalEquations'},expr_def,[],'System of equations'), ...
    arg({'morder','ModelOrder'},3,[1 Inf],'Model order. This is mandatory'));

if isempty(morder)
    error('SIFT:sim_examples:badParam','ModelOrder must be specified');
end

repo: REST


text:
function L=comp_framelength_fusion(F,Ls);

% This is highly tricky: Get the minimal transform length for each
% subframe, and set the length as the lcm of that.
Lsmallest=1;
for ii=1:F.Nframes
    Lsmallest=lcm(Lsmallest,framelength(F.frames{ii},1));
end;
L=ceil(Ls/Lsmallest)*Lsmallest;

% Verify that we did not screw up the assumptions.
for ii=1:F.Nframes
    if L~=framelength(F.frames{ii},L)
        error(['%s: Cannot determine a frame length. Frame no. %i does ' ...
               'not support a length of L=%i.'],upper(mfilename),ii,L);
    end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function zeta_test ( )

%*****************************************************************************80
%
%% ZETA_TEST tests ZETA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 October 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZETA_TEST\n' );
  fprintf ( 1, '  ZETA evaluates the Riemann Zeta function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       N                 Zeta                  Zeta\n' );
  fprintf ( 1, '                        exact              computed\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, z1 ] = zeta_values ( n_data );

    if ( n_data == 0 )
      break
    end

    z2 = zeta ( n );

    fprintf ( 1, '  %6d  %20e  %20e\n', n, z1, z2 );

  end

  return
end


text:
function A = dh_calc(a,alpha,d,theta)
% DH_CALC
%
% Calculates the homogeneous transformation between two consecutive joint frames
%
% Input:  [a, alpha, d, theta] - classical DH parameters (not modified DH)
% Output: A - homogeneous transformation between two consecutive joint frames

v = [ a*cos(theta), a*sin(theta), d ];
Xx = cos(theta); Yx = -sin(theta) * cos(alpha); Zx =  sin(theta) * sin(alpha);
Xy = sin(theta); Yy =  cos(theta) * cos(alpha); Zy = -cos(theta) * sin(alpha);
Xz = 0.0;        Yz =  sin(alpha);              Zz =  cos(alpha);
A = [Xx, Yx, Zx, v(1);
     Xy, Yy, Zy, v(2);
     Xz, Yz, Zz, v(3);
     0,  0,  0,  1];
end

author: xuhuairuogu

repo: V-REP-Simulation-Projects

path: github-repos/MATLAB/xuhuairuogu-V-REP-Simulation-Projects/V-REP-Simulation-Projects-841b944af4ea3a8fb250578d36434515f577f411/admittance control(adapted from an example in the book--A Systematic Approach to Learning Robot Programming with ROS)/iiwa14_kinematics/dh_calc.m


text:
function [ i, j, S ] = pendulumDF( y )
%pendulumDF Jacobian of the pendulum equations.

    g = 1;
    m = 1;
    l = 1;

    it = [ 1  3         3          2  4         4         5      5      5      5];
    jt = [ 3  5         1          4  5         2         1      2      3      4 ];
    S = [1  -y(1,1)/m -y(5,1)/m  1  -y(2,1)/m -y(5,1)/m y(3,1) y(4,1) y(1,1) y(2,1) ];
    i = it;
    j = jt;

end


text:
function [E, N, U] = ch_LLA2ENU(lat, lon, h, lat0, lon0, h0)

% 经纬高 转 ENU
% lat0, lon0, h0: 起始点经纬高, 经纬度为rad， 高度为m
% lat, lon, h 终点经纬高, 经纬度为rad， 高度为m
% E, N ,U 系下增量，单位为m

%精确算法
% XYZ0 = ch_LLA2ECEF(lat0, lon0, h0);
% XYZ1 = ch_LLA2ECEF(lat, lon, h);
% dXYZ = XYZ1 - XYZ0;
%
%  [~, ~, C_ECEF2ENU, ~]= ch_earth(lat0, lon0, h0);
%  dENU = C_ECEF2ENU * dXYZ;
%  E = dENU(1);
%  N= dENU(2);
%  U = dENU(3);

 %近似算法
R_0 = 6378137; %WGS84 Equatorial radius in meters
clat = cos(lat0);
E = (lon - lon0) * clat * R_0;
N =  (lat - lat0) * R_0;
U = h - h0;
end

author: yandld


text:
classdef mlcnn
% Convolutional Multilayer Neural Network
%------------------------------------------------------------------------------
% Initialize, train, and test a multilayer convolutional neural network.
%
% Supports multiple activation functions including linear, sigmoid, tanh, and
% soft rectification (softplus).
%
% Also supports mean squared error (mse), binary (xent), and multi-class
% (mcxent) cross-entropy cost functions.
%
% Supports early stopping regularization.
%------------------------------------------------------------------------------
% DES
% stan_s_bury@berkeley.edu

properties
	class = 'mlcnn';
	nLayers;				% # OF UNIT LAYERS
	layers ;				% LAYER STRUCTS

	netOut = [];			% CURRENT OUTPUT OF THE NETWORK
	netError;				% CURRENT NETWORK OUTPUT ERROR

	costFun = 'mse';		% COST FUNCTION
	J = [];					% CURRENT COST FUNCTION VALUE

	nEpoch = 10;			% TOTAL NUMBER OF TRAINING EPOCHS
	epoch = 1;				% CURRENT EPOCH

	batchSize = 20;			% SIZE OF MINIBATCHES
	trainBatches = [];		% BATCH INDICES
	xValBatches = [];		% XVAL. DATA INDICES

	trainCost = [];			% TRAINING ERROR HISTORY (ALL DATA)
	xValCost = [];			% XVALIDATION EROR HISTORY
	bestxValCost = Inf		% TRACK BEST NETWORK ERROR

	nXVal = 0;				% PORTION OF DATA HELD OUT FOR XVALIDATION
	stopEarly = 5;			% # OF EARLY STOP EPOCHS (XVAL ERROR INCREASES)
	stopEarlyCnt = 0;		% EARLY STOPPING CRITERION COUNTER
	bestNet = [];			% STORAGE FOR BEST NETWORK

	denoise = 0;			% PROPORTION OF VISIBLE UNIT DROPOUT (SIGMOID ONLY)
	dropout = 'not used';	% (POTENTIAL) PROP. OF HID. UNIT DROPOUT (SIGMOID ONLY)

	wPenalty = 0;			% WEIGHT DECAY TERM
	beginWeightDecay=1;
	momentum = 'not used';	% (POTENTIAL) MOMENTUM

	saveEvery = realmax;	% SAVE PROGRESS EVERY # OF EPOCHS
	saveDir;
	visFun;					% VISUALIZATION FUNCTION HANDLE
	trainTime = Inf;		% TRAINING DURATION
	verbose = 500;			% DISPLAY THIS # OF WEIGHT UPDATES
	auxVars = [];			% AUXILIARY VARIABLES (VISUALIZATION, ETC)
	useGPU = 0;				% USE THE GPU, IF AVAILABLE
	gpuDevice = [];			% STORAGE FOR GPU DEVICE
	checkGradients = 0;		% NUMERICAL GRADIENT CHECKING
end

methods
	function self = mlcnn(arch)
	% net = mlnn(arch)
	%--------------------------------------------------------------------------
	%mlnn constructor method. Initilizes a mlnn object, <net> given a user-
	%provided architecture, <arch>.
	%--------------------------------------------------------------------------
		self = self.init(arch);
	end

	function print(self)
		properties(self)
		methods(self)
	end

	function self = train(self, data, targets)
	%net = train(data, targets)
	%--------------------------------------------------------------------------
	% Train a convolutional neural net using stochastic gradient descent.
	% INPUT:
	%     <data>:  - [#PixelsY x #PixelsX x #Channels x #Obs];
	%
	%  <targets>:  - [#Output x #Obs]
	%--------------------------------------------------------------------------

		% DISTRIBUTE VALUES TO THE GPU
		if self.useGPU
			self = gpuDistribute(self);
		end

		self = self.makeBatches(data);
		self.trainCost = zeros(self.nEpoch,1);
		self.xValCost = zeros(self.nEpoch,1);
		nBatches = numel(self.trainBatches);
		tic; cnt = 1;

		if self.checkGradients,	checkNNGradients(self,data,targets); end

		% MAIN
	    while 1
			if self.verbose, self.printProgress('epoch'); end
			batchCost = zeros(nBatches,1);
			wPenalty = self.wPenalty;

			if self.epoch >= self.beginWeightDecay
				self.wPenalty = wPenalty;
			else
				self.wPenalty = 0;
			end

			for iB = 1:nBatches
				% GET BATCH DATA
				batchIdx = self.trainBatches{iB};
				netInput = data(:,:,:,batchIdx);
				netTargets = (targets(:,batchIdx));

				% ADD NOISE TO INPUT (DENOISING CAE)
				if self.denoise > 0
				    netInput = netInput.*(rand(size(netInput))>self.denoise);
				end

				% BACKPROP MAIN
				self = self.fProp(netInput, netTargets);
				self = self.bProp;
				self = self.updateParams;

				% ASSESS BATCH COST
				batchCost(iB) = self.J;
				cnt = cnt + 1;
			end

	        % AVERAGE COST OVER ALL TRAINING POINTS
			self.trainCost(self.epoch) = mean(batchCost);

	        % CROSS-VALIDATION
			if ~isempty(self.xValBatches)
				self = self.crossValidate(data,targets);
				self = self.assessNet;
				if self.verbose, self.printProgress('xValCost');end
			end

			% EARLY STOPPING
			if (self.epoch >= self.nEpoch) || ...
				(self.stopEarlyCnt >= self.stopEarly),
				self.trainCost = self.trainCost(1:self.epoch);
				self.xValCost = self.xValCost(1:self.epoch);
				break;
			end

			% SAVE BEST NETWORK
			if ~mod(self.epoch,self.saveEvery) & ~isempty(self.saveDir)
				self.save;
			end

			% DISPLAY
			if self.verbose
				self.printProgress('trainCost');
				if ~mod(cnt,self.verbose);
					self.visLearning;
				end
			end
			self.epoch = self.epoch + 1;
	    end
		self.trainTime = toc;

		% RETURN VALUES FROM THE GPU
		if self.useGPU
			self = gpuGather(self);
			reset(self.gpuDevice);
			self.gpuDevice = [];
		end
	end

	function self = fProp(self,netInput, targets)
	%[net,out] = fProp(netInput, targets)
	%--------------------------------------------------------------------------
	%Forward propagation of input signals, <netInput>. Also updates state of
	%network cost, if provided with <targets>. Also returns network output
	%<out>, if requested.
	%--------------------------------------------------------------------------
		if notDefined('targets'), targets = []; end
		nObs = size(netInput, 1);
		self.layers{1}.fm = netInput;

		for lL = 2:self.nLayers
			switch self.layers{lL}.type
			case 'conv'
			% LOOP OVER LAYER FEATURE MAPS
			for jM = 1:self.layers{lL}.nFM

				[nInY,nInX,nInFM,nObs] = size(self.layers{lL-1}.fm);

				% INITIALIZE LAYER MAP -- [nY,nX,nM,nObs]
				featMap = zeros([self.layers{lL}.fmSize,1,nObs]);

				% POOL OVER INPUT FEATURE MAPS,
				% CALC LAYER PRE-ACTIVATION
				for iM = 1:self.layers{lL-1}.nFM
					featMap = featMap + convn(self.layers{lL-1}.fm(:,:,iM,:), ...
					self.layers{lL}.filter(:,:,iM,jM),'valid');
				end

				if any(isnan(self.layers{lL}.b(jM))), keyboard, end

				% ADD LAYER BIAS
				featMap = featMap + self.layers{lL}.b(jM);

				% COMPLETE FEATURE MAP
				self.layers{lL}.fm(:,:,jM,:) = self.calcAct(featMap,self.layers{lL}.actFun);

			end

			case 'subsample'
				stride = self.layers{lL}.stride;
				% DOWNSAMPLE THE FEATURE MAPS FROM LAYER (l-1)
				for jM = 1:self.layers{lL-1}.nFM
					layerIn = self.layers{lL-1}.fm(:,:,jM,:);
					self.layers{lL}.fm(:,:,jM,:) = self.DOWN(layerIn,stride);
				end

			case 'output'
				% UNPACK OUTPUT FEATURES & CALCULATE NETWORK OUTPUT
				self = self.calcOutput;
			case 'rect'
			case 'lcn'
			case 'pool'
			end
		end
		% COST FUNCTION & OUTPUT ERROR SIGNAL
		if ~isempty(targets)
			[self.J, self.netError] = self.cost(targets,self.costFun);
		end
		if nargout > 1
			out = self.layers{end}.act;
		end
	end

	function self = calcOutput(self);
	% net = calcOutput();
	%--------------------------------------------------------------------------
	% Calculate the network ouput given the current parameters
	%--------------------------------------------------------------------------

		[nY,nX,nM,nObs]= size(self.layers{end-1}.fm);

		% # OF ENTRIES IN EACH FEATURE MAP
		nMap = prod([nY,nX]);

		% INITIALIZE OUTPUT FEATURES
		self.layers{end}.features = zeros(nMap*nM,nObs);

		% UNPACK MAPS INTO A MATRIX FOR CALCULATING OUTPUT
		for jM = 1:self.layers{end-1}.nFM
			map = self.layers{end-1}.fm(:,:,jM,:);
			self.layers{end}.features((jM-1)*nMap+1:jM*nMap,:) = reshape(map,nMap,nObs);
		end

		% CALC NET OUTPUTS
		preAct = bsxfun(@plus,self.layers{end}.W* ...
						self.layers{end}.features, ...
		                self.layers{end}.b);

		self.netOut = self.calcAct(preAct,self.layers{end}.actFun);

		if any(isnan(self.netOut)), keyboard; end
	end

	function [J, dJ] = cost(self,targets,costFun)
	% [J, dJ] = cost(self,targets,costFun)
	%--------------------------------------------------------------------------
	% Calculate the cost function <J> and derivative thereof for a set of targets
	% and the current state of the network.
	%--------------------------------------------------------------------------

		netOut = self.netOut;

		[nTargets,nObs] = size(netOut);
		switch costFun
		case 'mse' % REGRESSION
			delta = targets - netOut;
			J = 0.5*sum(sum(delta.^2))/nObs;
			dJ = -delta;

		case 'xent' % BINARY CLASSIFICATION
			J = -sum(sum(targets.*log(netOut) + (1-targets).*log(1-netOut)))/nObs;
			dJ = (netOut - targets)./(netOut.*(1-netOut));

		case 'mcxent' % MULTI-CLASS CLASSIFICATION (UNDER DEVO)
			class = softMax(netOut);
			J = -sum(sum(targets.*log(class)))/nObs;
			dJ = sum(labels - targets);

		case {'class','classerr'}  % CLASSIFICATION ERROR (WINNER TAKE ALL)

			[~, class] = max(netOut,[],1);
			[~, t] = max(targets,[],1);
			J = sum((class ~= t))/nObs;
			dJ = 'no gradient';
		case {'correlation','cc'}
			J = corr2(netOut,targets);
			dJ = 'no gradient';
		end
		if any(isnan(self.J)), keyboard, end
	end

	function self = bProp(self)
	%net = bProp()
	%--------------------------------------------------------------------------
	%Perform gradient descent no the loss w.r.t. each of the model parameters
	%using the backpropagation algorithm. Returns updated network object, <net>
	%--------------------------------------------------------------------------

		% DERIVATIVE OF OUTPUT ACTIVATION FUNCTION
		dAct = self.calcActDeriv(self.netOut,self.layers{end}.actFun);

		% OUTPUT ERROR SIGNAL -- [#Out x #Obs]
		outES = self.netError.*dAct;

		% ERROR SIGNAL  -- [#Features x #Obs]
		es = self.layers{end}.W'*outES;

		% IN CASE LAST FEATURE LAYER IS CONV.
		if strcmp(self.layers{end-1}.type,'conv')
			dAct = self.calcActDeriv(self.layers{end}.features,self.layers{end-1}.actFun);
			es = es.*dAct;
		end

		% REPACK ERROR SIGNAL INTO 2-D FEATURE MAP REPRESENTATION
		[nY,nX,nM,nObs] = size(self.layers{end-1}.fm);
		self.layers{end-1}.es = zeros([nY,nX,nM,nObs]);

		nMap = prod([nY*nX]); % NUMBER OF ENTRIES PER 2D MAP

		for jM = 1:self.layers{end-1}.nFM
			self.layers{end-1}.es(:,:,jM,:) = ...
			reshape(es((jM-1)*nMap+1:jM*nMap,:),[nY,nX,1,nObs]);
		end

		% BACKPROPATE ERROR SIGNAL
		for lL = self.nLayers-2:-1:2
			switch self.layers{lL}.type
			case 'conv'
				stride = self.layers{lL + 1}.stride;
				mapSz = size(self.layers{lL}.fm);
				self.layers{lL}.es = zeros(mapSz);

				for jM = 1:self.layers{lL}.nFM
					switch self.layers{lL+1}.type
					case 'subsample'
						% UPSAMPLE ES FROM ABOVE SUBSAMPLE LAYER
						propES = self.UP(self.layers{lL+1}.es(:,:,jM,:), ...
						                 [stride(1), stride(2),1,1])/prod(stride);
					case 'rect'
					case 'lcn'
					end

                   % DERIVATIVE OF ACTIVATION FUNCTION
                   dAct = self.calcActDeriv(self.layers{lL}.fm(:,:,jM,:), ...
											self.layers{lL}.actFun);

                   % CALCULATE LAYER ERROR SIGNAL
					self.layers{lL}.es(:,:,jM,:) = propES.*dAct;
				end
			case 'rect'
			case 'lcn'
			case 'pool'

			case 'subsample'
				[nY,nX,nM,nObs] = size(self.layers{lL}.fm);
				self.layers{lL}.es = zeros([nY,nX,nM,nObs]);
				for jM = 1:self.layers{lL}.nFM
					% FORM FEATURE MAP ERROR SIGNAL
					propES = zeros(nY,nX,1,nObs);
					for kM = 1:self.layers{lL+1}.nFM
						rotFilt = self.ROT(self.layers{lL+1}.filter(:,:,jM,kM));
						es = self.layers{lL+1}.es(:,:,kM,:);
						propES = propES + convn(es,rotFilt,'full');
					end
					self.layers{lL}.es(:,:,jM,:) = propES;
					if any(isnan(self.layers{lL}.es(:))), keyboard; end
					if any(isinf(self.layers{lL}.es(:))), keyboard; end
				end
			end
		end

		% CALCULATE THE GRADIENTS
		for lL = 2:self.nLayers-1
			[nX,nY,nM,nObs] = size(self.layers{lL}.fm);
			switch self.layers{lL}.type
			case 'conv'
				for jM = 1:self.layers{lL}.nFM
					es = self.layers{lL}.es(:,:,jM,:);
					for iM = 1:self.layers{lL-1}.nFM
						input = self.FLIPDIMS(self.layers{lL-1}.fm(:,:,iM,:));
						dEdFilter = convn(input,es,'valid')/nObs;
						self.layers{lL}.dFilter(:,:,iM,jM) = dEdFilter;
						if isnan(any(dEdFilter)), keyboard; end
					end
					self.layers{lL}.db(jM) = sum(es(:))/nObs;
				end
			case 'rect'
			case 'lcn'
			case 'pool'
			end
		end

		% GRADIENTS FOR OUTPUT LAYER WEIGHTS AND BIASES
		self.layers{end}.dW = outES*self.layers{end}.features'/nObs;
		self.layers{end}.db = mean(outES,2);

	end

	function self = updateParams(self)
	%net = updateParams()
	%--------------------------------------------------------------------------
	%Update network parameters based on states of netowrk gradient, perform
	%regularization such as weight decay and weight rescaling
	%--------------------------------------------------------------------------
		wPenalty = 0;
		for lL = 2:self.nLayers-1
			switch self.layers{lL}.type
			% CURRENTLY, ONLY UPDATE FILTERS AND FM BIASES
			% PERHAPS, IN THE FUTURE, WE'LL BE FANCY, AND DO FANCY UPDATES
			case {'conv','output'}
				lRate = self.layers{lL}.lRate;
				for jM = 1:self.layers{lL}.nFM
					% UPDATE FEATURE BIASES
					self.layers{lL}.b(jM) = self.layers{lL}.b(jM) - ...
					                    lRate*self.layers{lL}.db(jM);

					% UPDATE FILTERS
					for iM = 1:self.layers{lL-1}.nFM
						if self.wPenalty > 0 % L2 REGULARIZATION
							wPenalty = self.layers{lL}.filter(:,:,iM,jM)*self.wPenalty;
						elseif self.wPenalty < 0 % L1-REGULARIZATION (SUB GRADIENTS)
							wPenalty = sign(self.layers{lL}.filter(:,:,iM,jM))*abs(self.wPenalty);
						end
						self.layers{lL}.filter(:,:,iM,jM) = ...
						self.layers{lL}.filter(:,:,iM,jM) - ...
						lRate*(self.layers{lL}.dFilter(:,:,iM,jM)+wPenalty);
					end
				end
			end
	    end
	end

	function out = calcAct(self,in,actFun)
	%out = calcAct(in,actFun)
	%--------------------------------------------------------------------------
	%Calculate the output activation <out> from an input <in> for activation
	%function <actFun>. Available activation functions include 'linear','exp',
	%'sigmoid', 'softmax', 'tanh', and 'softrect'.
	%--------------------------------------------------------------------------

		switch actFun
			case 'linear'
				out = self.stabilizeInput(in,1);

			case 'exp'
				in = self.stabilizeInput(in,1);
				out = exp(in);

			case 'sigmoid'
				in = self.stabilizeInput(in,1);
				out = 1./(1 + exp(-in));

			case 'softmax'
				in = self.stabilizeInput(in,1);
				maxIn = max(in, [], 2);
				tmp = exp(bsxfun(@minus,in,maxIn));
				out = bsxfun(@rdivide,tmp,sum(tmp,2));

			case 'tanh'
				in = self.stabilizeInput(in,1);
				out = tanh(in);

			case 'softrect'
				k = 8;
				in = self.stabilizeInput(in,k);
				out = 1/k.*log(1 + exp(k*in));
		end
	end

	function dAct = calcActDeriv(self,in,actFun)
	%dAct = calcActDeriv(in,actFun)
	%--------------------------------------------------------------------------
	%Calculate the output activation derivatives <dAct> from an input <in> for
	%activation function <actFun>. Available activation functions derivatives
	% include 'linear','exp', sigmoid','tanh', and 'softrect'.
	%--------------------------------------------------------------------------

		switch actFun
			case 'linear'
				dAct = ones(size(in));

			case 'exp';
				in = self.stabilizeInput(in,1);
				dAct = in;

			case 'sigmoid'
				in = self.stabilizeInput(in,1);
				dAct = in.*(1-in);

			case 'tanh'
				in = self.stabilizeInput(in,1);
				dAct = 1 - in.^2;

			case 'softrect'
				k = 8;
				in = self.stabilizeInput(in,k);
				dAct = 1./(1 + exp(-k*in));
		end
	end

	function out = DOWN(self,data,stride)
	%out = DOWN(data,stride)
	%--------------------------------------------------------------------------
	% Downsample 1st 2 dimensions of a tensor
	%--------------------------------------------------------------------------
		tmp = ones(stride(1),stride(2));
		tmp = tmp/prod(stride(:));
		out = convn(data,tmp,'valid');
		out = out(1:stride(1):end,1:stride(2):end,:,:,:);
	end

	function out = UP(self,data,scale);
	%out = UP(data,stride)
	%--------------------------------------------------------------------------
	% Upsample 1st 2 dimensions of a tensor
	%--------------------------------------------------------------------------
	% UPSAMPLE DIMESIONS OF A TENSOR
		dataSz = size(data);
		idx = cell(numel(dataSz),1);
		for iD = 1:numel(dataSz)
			tmp = zeros(dataSz(iD)*scale(iD),1);
			tmp(1:scale(iD):dataSz(iD)*scale(iD)) = 1;
			idx{iD} = cumsum(tmp);
		end
		out = data(idx{:});
	end

	function out = ROT(self,out)
	%out = ROT(self,in)
	%--------------------------------------------------------------------------
	% Rotate the 1st two dimensions of a tensor by one-half rotation
	%--------------------------------------------------------------------------
		out = out(end:-1:1,end:-1:1,:,:);
	end

	function out = FLIPDIMS(self,out)
	%out = FLIPDIMS(in)
	%--------------------------------------------------------------------------
	% Flip all dimensions of a tensor <in>
	%--------------------------------------------------------------------------
		for iD = 1:numel(size(out))
			out = flipdim(out,iD);
		end
	end

	function self = makeBatches(self,data);
	% net = makeBatches(data);
	%--------------------------------------------------------------------------
	% Create batches based on data. Observations are along the rows of data.
	% (ASSUME THAT data IS RANDOMIZED ACROSS SAMPLES (ROWS))
	%--------------------------------------------------------------------------

		nObs = size(data,4);

		if self.nXVal < 1
			nVal = round(self.nXVal*nObs);
		else
			nVal = self.nXVal;
		end

		xValIdx = nObs-nVal+1:nObs;
		nObs = nObs - nVal;
		nBatches = ceil(nObs/self.batchSize);
		idx = round(linspace(1,nObs+1,nBatches+1));

		for iB = 1:nBatches
			if iB == nBatches
				batchIdx{iB} = idx(iB):nObs;
			else
				batchIdx{iB} = idx(iB):idx(iB+1)-1;
			end
		end
		self.trainBatches = batchIdx;

		nxBatches = ceil(nVal/self.batchSize);
		if ~isempty(xValIdx)
			xValIdx = round(linspace(xValIdx(1),xValIdx(end)+1,nxBatches));

			for iB = 1:nxBatches-1
				if iB == nxBatches
					xBatchIdx{iB} = xValIdx(iB):nVal;
				else
					xBatchIdx{iB} = xValIdx(iB):xValIdx(iB+1)-1;
				end
			end
			self.xValBatches = xBatchIdx;
		end
	end

	function self = crossValidate(self,data,targets)
	%net = crossValidate(data,targets)
	%--------------------------------------------------------------------------
	%Run cross-validation on current model parameters.
	%--------------------------------------------------------------------------

		xValCost = 0;
		for iB = 1:numel(self.xValBatches)
			idx = self.xValBatches{iB};
			tmpCost = self.test(data(:,:,:,idx),targets(:,idx),self.costFun);
			xValCost = xValCost + mean(tmpCost);
		end
		% AVERAGE data-VALIDATION ERROR
		self.xValCost(self.epoch) = xValCost/iB;
	end

	% ASSES PREDICTIONS/ERRORS ON TEST DATA
	function [cost,pred] = test(self,data,targets,costFun)
		if notDefined('costFun') costFun=self.costFun; end

		for lL = 1:self.nLayers
			try
				self.layers{lL}.fm = [];
			catch
			end
		end
		self = self.fProp(data,targets);
		pred = self.netOut;
		cost = self.cost(targets,costFun);
	end

	function self = assessNet(self)
	%assessNet()
	%--------------------------------------------------------------------------
	%Utility function to assess the quality of current netork parameters and
	%store net, if necessary.
	%--------------------------------------------------------------------------

		if self.epoch > 1
			if self.xValCost(self.epoch) < self.bestxValCost
				self.bestNet = self.layers;
				self.bestxValCost = self.xValCost(self.epoch);
				self.stopEarlyCnt = 0;
			else
				self.stopEarlyCnt = self.stopEarlyCnt + 1;
			end
		else
			self.bestNet = self.layers; % STORE FIRST NET BY DEFAULT
		end
	end

	function printProgress(self,type)
	%printProgress(type)
	%--------------------------------------------------------------------------
	%Verbose utility function. <type> is the type of message to print.
	%--------------------------------------------------------------------------
		switch type
		case 'epoch'
			fprintf('Epoch: %i/%i',self.epoch,self.nEpoch);
		case 'trainCost'
			fprintf('\t%s: %2.3f\n',self.costFun,self.trainCost(self.epoch));
		case 'time'
			fprintf('\tTime: %g\n', toc);
		case 'xValCost'
			if ~self.stopEarlyCnt
				fprintf('\tCrossValidation Error:  %g (best net) \n',self.xValCost(self.epoch));
			else
				fprintf('\tCrossValidation Error:  %g\n',self.xValCost(self.epoch));
			end
		case 'gradCheck'
			netGrad = self.auxVars.netGrad;
			numGrad = self.auxVars.numGrad;
			gradFailed = self.auxVars.gradFailed;
			switch gradFailed
				case 1, gradStr = '(Failed)';
				otherwise, gradStr = '(Passed)';
			end
			fprintf('\tNetwork = %2.6f  \t Numerical = %2.6f  %s\n' ,netGrad,numGrad,gradStr);
		case 'save'
			fprintf('\nSaving...\n\n');
		end
	end

	function self = init(self,arch)
	%net = init(arch)
	%--------------------------------------------------------------------------
	%Utility function to used intitialize a neural network given an architecture
	%<arch> is a cell array of structs, one for each layer in the network. The
	%fields of each structure will depend on the type of layer. Supported lay-
	%ers include 'input','conv','subsample','output'.
	%
	%Note, the first and last layers should be of 'type' 'input' and 'output'
	%respectively.
	%
	% Returns a mlcnn object, <net>.
	%--------------------------------------------------------------------------
	% INTITIALIZE A NEURAL NETWORK GIVEN AN ARCHITECTURE
	% <arch> IS A CELLARRAY OF STRUCTS WITH LAYER-SPECIFIC PARAMETERS

		arch = self.ensureArchitecture(arch);
		self.nLayers = numel(arch);

		% INITIALIZE LAYERS FROM ARCHITECTURE
	    for lL = 1:numel(arch)
		    self.layers{lL}.type = arch{lL}.type;
		    switch arch{lL}.type
		    case 'input'
			    self.layers{lL}.dataSize = arch{lL}.dataSize;
			    self.layers{lL}.fmSize = arch{lL}.dataSize(1:2);
			    self.layers{lL}.nFM = arch{lL}.dataSize(3);

			case 'conv'
				if strcmp(arch{lL-1}.type,'input')
					nInY = arch{lL-1}.dataSize(1);
					nInX = arch{lL-1}.dataSize(2);
					nInFM = arch{lL-1}.dataSize(3);
				else
					nInX = self.layers{lL-1}.fmSize(2);
					nInY = self.layers{lL-1}.fmSize(1);
					nInFM = self.layers{lL-1}.nFM;

				end
				% INTERMEDIATE VARIABLES
				nFM = arch{lL}.nFM;
				filtSize = arch{lL}.filterSize;
				fmSize = [nInY,nInX] - filtSize + 1;
				fanIn = nInFM*prod(filtSize);
				fanOut = nFM*prod(filtSize);
				range = 2*sqrt(6/((fanIn + fanOut)));

				% INITIALIZE LAYER PARAMETERS
				self.layers{lL}.actFun = arch{lL}.actFun;
				self.layers{lL}.lRate = arch{lL}.lRate;
				self.layers{lL}.filterSize = filtSize;
				self.layers{lL}.fmSize = fmSize;
				self.layers{lL}.nFM = nFM;
				self.layers{lL}.fm = [];
				self.layers{lL}.filter = range*(rand([filtSize,nInFM,nFM])-.5);
				self.layers{lL}.dFilter = zeros(size(self.layers{lL}.filter));
				self.layers{lL}.b = zeros(nFM,1);
				self.layers{lL}.db = zeros(nFM,1);

			case 'subsample'
				% INTERMEDIATE VARIABLES
				stride = arch{lL}.stride;
				fmSize = floor(self.layers{lL-1}.fmSize./stride);
				nFM = self.layers{lL-1}.nFM;

				% INITIALIZE LAYER PARAMETERS
				self.layers{lL}.stride = stride;
				self.layers{lL}.fmSize = fmSize;
				self.layers{lL}.nFM = nFM;
				self.layers{lL}.fm = [];
				self.layers{lL}.b = zeros(nFM,1);
				self.layers{lL}.db = zeros(nFM,1);

			case 'output'

				nOut = arch{lL}.nOut;
				nFMIn = self.layers{lL-1}.nFM;
				nOutFeats = prod([self.layers{lL-1}.fmSize,nFMIn]);
				range = 2*sqrt(6/(nOutFeats + nOut));

				% ADJUST NETWORK OUTPUT LAYER PARAMETERS
				self.layers{lL}.nOut = nOut;
				self.layers{lL}.actFun = arch{lL}.actFun;
				self.layers{lL}.lRate = arch{lL}.lRate;
				self.layers{lL}.W = range*(rand(nOut,nOutFeats)-.5);
				self.layers{lL}.dW = zeros(size(self.layers{lL}.W));
				self.layers{lL}.b = zeros(nOut,1);
				self.layers{lL}.db = zeros(nOut,1);
		    end
		end
	end

	function arch = ensureArchitecture(self,arch)
	%arch = ensureArchitecture(arch)
	%--------------------------------------------------------------------------
	%Utility function to reprocess a supplied architecture, <arch>
	%--------------------------------------------------------------------------

		if ~iscell(arch), error('<arch> needs to be a cell array of layer params');end
		if ~strcmp(arch{1}.type,'input'), error('define an input layer'); end
		if ~strcmp(arch{end}.type,'output'), error('define an output layer'); end

		% ENSURE LAYER-SPECIFIC PARAMS
		for lL = 1:numel(arch)
			lParams = fields(arch{lL});
			switch arch{lL}.type
			case 'input'
				if ~any(strcmp(lParams,'dataSize')), error('must provide data size');end

			case 'conv'
				if ~any(strcmp(lParams,'filterSize')) || isempty(arch{lL}.filterSize)
					arch{lL}.filterSize = [5 5];
				elseif numel(arch{lL}.filterSize) == 1;
					arch{lL}.filterSize = repmat(arch{lL}.filterSize,[1,2]);
				end
				if ~any(strcmp(lParams,'nFM'));
					arch{lL}.nFM = 10;
				end
				if ~any(strcmp(lParams,'lRate'));
					arch{lL}.lRate = .5;
				end
				if ~any(strcmp(lParams,'actFun'));
					arch{lL}.actFun = 'sigmoid';
				end

			case 'subsample'
				if ~any(strcmp(lParams,'stride')) || isempty(arch{lL}.stride);
					arch{lL}.stride = [2 2];
				elseif numel(arch{lL}.stride) == 1;
					arch{lL}.stride = repmat(arch{lL}.stride,[1,2]);
				end
			case 'output'
				if ~any(strcmp(lParams,'nOut'))
					error('must provide number of outputs');
				end
				if ~any(strcmp(lParams,'actFun'));
					arch{lL}.actFun = 'sigmoid';
				end
				if ~any(strcmp(lParams,'lRate'));
					arch{lL}.lRate = .5;
				end
			case 'rect'
			case 'lcn'
			case 'pool'
			end
		end
	end

	function in = stabilizeInput(self,in,k);
	%in = stabilizeInput(in,k);
	%--------------------------------------------------------------------------
	%Utility function to ensure numerical stability. Clips values of <in>
	%such that exp(k*in) is within single numerical precision.
	%--------------------------------------------------------------------------
		cutoff = log(realmin('single'));
		in(in*k>-cutoff) = -cutoff/k;
		in(in*k<cutoff) = cutoff/k;
	end

	function visLearning(self)
	%visLearning()
	%--------------------------------------------------------------------------
	%Utility function to perform learning isualizations.
	%--------------------------------------------------------------------------
		try
			self.visFun(self);
		catch
			if ~isfield(self.auxVars,'printVisWarning')
				fprintf('\nWARNING: visualization failed.')
				self.auxVars.printVisWarning = true;
			end
		end
	end

	function save(self);
	%save();
	%--------------------------------------------------------------------------
	%Utility function to save current network.
	%--------------------------------------------------------------------------
		if self.verbose, self.printProgress('save'); end
		if ~isdir(self.saveDir)
			mkdir(self.saveDir);
		end
		fileName = fullfile(self.saveDir,sprintf('mlcnn.mat'));
		net = self.bestNet;
		save(fileName,'net');
	end
end % END METHODS
end % END CLASSDEF

author: dustinstansbury

repo: medal


text:
function [Y,I] = minnz(X)
  % MINNZ  find minimum nonzero entry in columns of X
  %
  % [Y,I] = minnz(X)
  %
  % Inputs:
  %   X  m by n sparse matrix
  % Outputs
  %   Y  n list of minimum non zero entries in each column of X
  %   I  n list of row indices to minimum non-zero entries in each column of X
  %
  % See also: min, max, maxnz
  %
  %

  % rebuild using 1/element
  [XI,XJ,XV] = find(X);
  X_i = sparse(XI,XJ,XV.^-1,size(X,1),size(X,2));
  [infX_i,infI_i] = max(X==inf);
  [maxX_i,maxI_i] = max(X_i);
  maxX_i(maxX_i==0) = inf;
  maxX_i(infX_i) = 0;
  maxI_i(infX_i) = infI_i(infX_i);
  [minX,minI] = min(X);
  minX(minX==0) = inf;
  Y = [maxX_i.^-1;minX];
  [Y,J] = min(Y);
  I = [maxI_i;minI];
  I = I(sub2ind(size(I),J,1:size(I,2)));

end

author: alecjacobson

repo: gptoolbox


text:
classdef CEC2010_F18 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{18};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 50;
            obj.upper    = zeros(1,obj.D) + 50;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum((Z(:,1:end-1)-Z(:,2:end)).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = mean((-Z.*sin(sqrt(abs(Z)))),2);
            PopCon(:,3) = abs(mean((-Z.*sin(sqrt(abs(Z)))),2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function dimensionMap = getMapOfDimensions(MRSIStruct)
    values = cell2mat(struct2cell(MRSIStruct.dims));
    dimsFieldNames = fieldnames(MRSIStruct.dims);
    zeroIndexs = values == 0;
    nonZeroFieldName = dimsFieldNames(~zeroIndexs);
    dimensionMap = containers.Map(values(~zeroIndexs), nonZeroFieldName);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A