text:
function varargout = adjustIntensityRangeGui(varargin)
% ADJUSTINTENSITYRANGEGUI M-file for adjustIntensityRangeGui.fig
%      ADJUSTINTENSITYRANGEGUI, by itself, creates a new ADJUSTINTENSITYRANGEGUI or raises the existing
%      singleton*.
%
%      H = ADJUSTINTENSITYRANGEGUI returns the handle to a new ADJUSTINTENSITYRANGEGUI or the handle to
%      the existing singleton*.
%
%      ADJUSTINTENSITYRANGEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ADJUSTINTENSITYRANGEGUI.M with the given input arguments.
%
%      ADJUSTINTENSITYRANGEGUI('Property','Value',...) creates a new ADJUSTINTENSITYRANGEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before adjustIntensityRangeGui_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to adjustIntensityRangeGui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help adjustIntensityRangeGui

% Last Modified by GUIDE v2.5 11-Aug-2004 09:30:47

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @adjustIntensityRangeGui_OpeningFcn, ...
                   'gui_OutputFcn',  @adjustIntensityRangeGui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
%--------------------------------------------------------
% --- Executes just before adjustIntensityRangeGui is made visible.
%--------------------------------------------------------
function adjustIntensityRangeGui_OpeningFcn(hObject, eventdata, handles, varargin)

% Choose default command line output for adjustIntensityRangeGui
handles.output = hObject;

hMainGui = getappdata(0,'hMainGui');
intensityParam = getappdata(hMainGui,'intensityParam');
set(handles.lowerEdit,'String',num2str(intensityParam.lowerLimit));
set(handles.upperEdit,'String',num2str(intensityParam.upperLimit));
handles.intensityParam = intensityParam;

% Set the slider's value and step size
handax=findobj(hMainGui,'tag','axes1');
cimg = getimage(handax);
cmin = intensityParam.minIntensity;
cmax = intensityParam.maxIntensity;
rng = cmax-cmin;

if cmax == 0
    cmin = 0; cmax = 1;
    step_size = 1;
elseif cmax <= 1
    step_size = double(rng/256);
else
    step_size = 1/double(rng-1);
end
set(handles.lowerSlider,'Min',cmin,'Max',cmax,...
    'Value', intensityParam.lowerLimit,'SliderStep',[step_size step_size]);
set(handles.upperSlider,'Min',cmin,'Max',cmax,...
    'Value',intensityParam.upperLimit,'SliderStep',[step_size step_size]);
set(handles.upperEdit,'String',num2str(intensityParam.upperLimit));
set(handles.lowerEdit,'String',num2str(intensityParam.lowerLimit));

guidata(hObject, handles);

%--------------------------------------------------------
% --- Outputs from this function are returned to the command line.
%--------------------------------------------------------
function varargout = adjustIntensityRangeGui_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

%--------------------------------------------------------
% --- Executes during object creation, after setting all properties.
%--------------------------------------------------------
function lowerSlider_CreateFcn(hObject, eventdata, handles)
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
%--------------------------------------------------------
function lowerSlider_Callback(hObject, eventdata, handles)
lwLimit= round(get(hObject,'Value'));
upLimit = round(get(handles.upperSlider,'Value'));
if lwLimit<upLimit
    set(handles.lowerEdit,'String',num2str(lwLimit));
    hMainGui = getappdata(0,'hMainGui');
    handax=findobj(hMainGui,'tag','axes1');
    ll = get(handax,'Clim');
    set(handax,'Clim',[lwLimit ll(2)]);
    handles.intensityParam.lowerLimit = lwLimit;
else
    errordlg('Lower limit cannot be larger than the upper');
end
guidata(hObject, handles);

%--------------------------------------------------------
% --- Executes during object creation, after setting all properties.
function upperSlider_CreateFcn(hObject, eventdata, handles)
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%--------------------------------------------------------
% --- Executes on slider movement.
%--------------------------------------------------------
function upperSlider_Callback(hObject, eventdata, handles)
upLimit= round(get(hObject,'Value'));
lwLimit = round(get(handles.lowerSlider,'Value'));
if lwLimit < upLimit
    set(handles.upperEdit,'String',num2str(upLimit));
    hMainGui = getappdata(0,'hMainGui');
    handax=findobj(hMainGui,'tag','axes1');
    ll = get(handax,'Clim');
    set(handax,'Clim',[ll(1) upLimit]);
    handles.intensityParam.upperLimit = upLimit;
else
    errordlg('Lower limit cannot be larger than upper limit');
end
guidata(hObject, handles);
%--------------------------------------------------------
% --- Executes on button press in pushbuttonApply.
%--------------------------------------------------------
function pushbuttonApply_Callback(hObject, eventdata, handles)
hMainGui = getappdata(0,'hMainGui');
setappdata(hMainGui,'intensityParam',handles.intensityParam);
% getappdata(hMainGui,'intensityParam')
if ishandle(handles.figure1),
    close(handles.figure1);
end
%--------------------------------------------------------
% --- Executes during object creation, after setting all properties.
%--------------------------------------------------------
function lowerEdit_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%--------------------------------------------------------
%--------------------------------------------------------
function lowerEdit_Callback(hObject, eventdata, handles)
strg = get(hObject,'String');
lwLimit =str2num(strg);
upLimit = get(handles.upperSlider,'Value');
maxLimit = get(handles.lowerSlider,'Max');
if lwLimit <= maxLimit
    if lwLimit < upLimit
        set(handles.lowerSlider,'Value',lwLimit);
        hMainGui = getappdata(0,'hMainGui');
        handax=findobj(hMainGui,'tag','axes1');
        ll = get(handax,'Clim');
        set(handax,'CLim',[lwLimit ll(2)]);
        handles.intensityParam.lowerLimit = lwLimit;
    else
        errordlg('Lower limit cannot be larger than upper limit');
    end
else
    errordlg('Larger than maximum limit');
end
guidata(hObject, handles);
%--------------------------------------------------------
% --- Executes during object creation, after setting all properties.
%--------------------------------------------------------
function upperEdit_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%--------------------------------------------------------
%--------------------------------------------------------
function upperEdit_Callback(hObject, eventdata, handles)
strg = get(hObject,'String');
upLimit =str2num(strg);
lwLimit = get(handles.lowerSlider,'Value');
maxLimit = get(handles.lowerSlider,'Max');
if upLimit <=maxLimit
    if upLimit > lwLimit
        set(handles.upperSlider,'Value',upLimit);
        hMainGui = getappdata(0,'hMainGui');
        handax=findobj(hMainGui,'tag','axes1');
        ll = get(handax,'CLim');
        set(handax,'CLim',[ll(1) upLimit]);
        handles.intensityParam.upperLimit = upLimit;
    else
        errordlg('Lower limit cannot be larger than the upper');
    end
else
    errordlg('Larger than maximum limit');
end
guidata(hObject, handles);

%--------------------------------------------------------
% --- Executes on button press in pushbuttonCancel.
function pushbuttonCancel_Callback(hObject, eventdata, handles)
close(handles.figure1);


text:
function test_bug2338

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY ft_prepare_bemmodel ft_prepare_headmodel ft_prepare_leadfield ft_compute_leadfield ft_headmodel_openmeeg

% 4 Layers
r = [85 88 92 100];
c = [1 1/20 1/80 1];

order = [1 2 3 4];

% Description of the spherical mesh
[pnt, tri] = mesh_sphere(42);

% Create a set of electrodes on the outer surface
elec.elecpos = max(r) * pnt;
elec.label = {};
nelec = size(elec.elecpos,1);
for ii=1:nelec
  elec.label{ii} = sprintf('vertex%03d', ii);
end

% Create one triangulated mesh for each boundary, the first boundary is inside
mesh = [];
for ii=1:length(r)
  mesh.bnd(ii).pnt = pnt * r(ii);
  mesh.bnd(ii).tri = tri;
end

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

% Compute the BEM model using the old code
cfg = [];
cfg.method = 'openmeeg';

cfg.conductivity = c(order);
mesh1 = mesh;
mesh1.bnd = mesh1.bnd(order);
vol1 = ft_prepare_bemmodel(cfg, mesh1);

% flip the geometry and conductivity around
order = [4 3 2 1];

cfg.conductivity = c(order);
mesh2 = mesh;
mesh2.bnd = mesh2.bnd(order);
vol2 = ft_prepare_bemmodel(cfg, mesh2);

cfg = [];
cfg.sourcemodel.pos = [0 0 70];
cfg.elec = elec;
cfg.headmodel = vol1;
lf1 = ft_prepare_leadfield(cfg);
cfg.headmodel = vol2;
lf2 = ft_prepare_leadfield(cfg);

assert(isalmostequal(lf1.leadfield{1}, lf2.leadfield{1}, 'reltol', 1e-6));

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

% Compute the BEM model using the new code
cfg = [];
cfg.method = 'openmeeg';

cfg.conductivity = c(order);
mesh1 = mesh;
mesh1.bnd = mesh1.bnd(order);
vol1 = ft_prepare_headmodel(cfg, mesh1);

% flip the geometry and conductivity around
order = [4 3 2 1];

cfg.conductivity = c(order);
mesh2 = mesh;
mesh2.bnd = mesh2.bnd(order);
vol2 = ft_prepare_headmodel(cfg, mesh2);

cfg = [];
cfg.sourcemodel.pos = [0 0 70];
cfg.elec = elec;
cfg.headmodel = vol1;
lf1 = ft_prepare_leadfield(cfg);
cfg.headmodel = vol2;
lf2 = ft_prepare_leadfield(cfg);

assert(isalmostequal(lf1.leadfield{1}, lf2.leadfield{1}, 'reltol', 1e-6));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% SIOWriter is a class implementing a specific FlatfileImageWriter for SICD
% files using an SIO "carrier".  The SICD metadata is written in a user
% data field called "SICDMETA" in the standard SIO file format.  The caller
% can also turn off the user data, since not all SIO readers handle this.
%
% Written by: Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////
classdef SIOWriter < FlatfileImageWriter

    methods
        % Constructor
        function obj = SIOWriter(filename, sicdmeta, varargin)
            % Parse inputs
            % These inputs are really parsed later in the superclass
            % constructor.  Here just need to 1) extract
            % include_sicd_metadata, which is specific to SIOWriter class
            % and 2) make sure that no additional arguments are passed--
            % like header_skip, which could break things, if passed on to
            % the superclass constructor.
            p = inputParser;
            p.addParamValue('data_type',0);
            p.addParamValue('is_complex',0);
            p.addParamValue('include_sicd_metadata',true,@(x) isscalar(x)&&islogical(x));
            p.FunctionName = 'SIOWriter';
            p.parse(varargin{:});

            % We'll use a fixed definition SIO header with (possibly) one user-data segment.
            % The header will the look like this:
            %    Core SIO header:
            %       Magic key   (4 byte uint, fixed value of 'FF027FFD')
            %       Rows        (4 byte uint)
            %       Columns     (4 byte uint)
            %       Data type   (4 byte uint)
            %       Data size   (4 byte uint, # bytes per element)
            %    Optional "user data":
            %       Num pairs   (4 byte uint, # pairs of user data, fixed at 1)
            %       Name bytes  (4 byte uint, # bytes in name of user element,
            %                    fixed at 8 for name "SICDMETA")
            %       Name        (8 bytes containing "SICDMETA")
            %       Value bytes (4 byte uint, value is length of XML string)
            %       Value       (XML string holding SICD metadata)
            header_skip = 5*4; % SIO header with no user data is 5 uint32 words
            if p.Results.include_sicd_metadata
                XML_meta_string = sicdstruct2xml(sicdmeta, 'inc_newline', false, ...
                    'inc_padding', false, 'file_type', 'SICD');
                header_skip = header_skip + 4 + 4 + 8 + 4 + length(XML_meta_string); % Add user data length
            end

            % Construct base object
            obj = obj@FlatfileImageWriter(filename, sicdmeta, 'header_skip', header_skip, varargin{:});

            % Open file and write SIO header
            obj.FID = fopen(filename,'w', 'b'); % We always write big-endian
            % Magic key
            if p.Results.include_sicd_metadata
                magic_key = hex2dec('FF027FFD'); % Indicates big endian, with user-data
            else
                magic_key = hex2dec('FF017FFE'); % Indicates big endian, with no user-data
            end
            fwrite(obj.FID, magic_key, 'uint32');
            % Rows and columns
            fwrite(obj.FID, obj.image_size(2), 'uint32');
            fwrite(obj.FID, obj.image_size(1), 'uint32');
            % Data type and size
            [ element_type, element_length ] = matlabtype2sio( obj.data_type, obj.is_complex );
            fwrite(obj.FID, element_type, 'uint32');
            fwrite(obj.FID, element_length, 'uint32');
            % User data
            if p.Results.include_sicd_metadata
                fwrite(obj.FID, 1, 'uint32');                       % Num pairs of user data
                fwrite(obj.FID, 8, 'uint32');                       % Pair 1 name length
                fwrite(obj.FID, 'SICDMETA', 'char');                % Pair 1 name
                fwrite(obj.FID, length(XML_meta_string), 'uint32'); % Pair 1 value length
                fwrite(obj.FID, XML_meta_string, 'char');           % Pair 1 value
            end
        end

        % Destructor
        function delete(obj)
            % If data has not been written to the end of the pre-defined
            % image data size, fill out rest of file space with zeros.
            image_data_size = (obj.is_complex + 1) * ...
                prod(double(obj.image_size)) * ...
                obj.data_size;
            obj.fseek(obj.FID,obj.header_skip+image_data_size,'bof');

            fclose(obj.FID);
        end
    end

end

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

author: ngageoint


text:
%{
Authors:
Jonatas Lopes de Paiva
Claudio Fabiano Motta Toledo
Helio Pedrini

%}

function f = crossover(p1, p2)

r = randi(3);
[linhas, colunas] = size(p1);
if (linhas == 0 || colunas == 0)
    fprinf('Error\n');
end

switch r

    case 1
        idx = randi(linhas);
        tmp = vertcat(p1(1:idx,:), p2(idx+1:linhas,:));
    case 2
        idx = randi(colunas);
        tmp = horzcat(p1(:,1:idx), p2(:,idx+1:colunas));
    case 3
        tmp = zeros(linhas, colunas);
        for i = 1: linhas
            for j = 1:colunas
                if(rand < 0.5)
                    tmp(i,j) = p1(i,j);
                else
                    tmp(i,j) = p2(i,j);
                end
            end
        end

end

f = tmp;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Feature sign search
% code by Wang Jinjun @ NEC Research Lab America
% reference
% Efficient sparse coding algorithms
%   Honglak Lee Alexis Battle Rajat Raina Andrew Y. Ng
%       Computer Science Department
%       Stanford University
%       Stanford, CA 94305

function [x]=feature_sign(B,y,lambda,init_x)

nbases=size(B,2);

OptTol = 1e-5;

if nargin < 4,
    x=zeros(nbases, 1);
else
    x = init_x;
end;

theta=sign(x);          %sign flag
a=(x~=0);               %active set

optc=0;

By=B'*y;
B_h=B(:,a);
x_h=x(a);
Bx_h=B_h*x_h;
all_d=2*(B'*Bx_h-By);
[ma mi]=max(abs(all_d).*(~a));

while optc==0,

    optc=1;

    if all_d(mi)>lambda+1e-10,
        theta(mi)=-1;
        a(mi)=1;
        b=B(:,mi);
        x(mi)=(lambda-all_d(mi))/(b'*b*2);
    elseif all_d(mi)<-lambda-1e-10,
        theta(mi)=1;
        a(mi)=1;
        b=B(:,mi);
        x(mi)=(-lambda-all_d(mi))/(b'*b*2);
    else
        if sum(a)==0,
            lambda=ma-2*1e-10;
            optc=0;
            b=B(:,mi);
            x(mi)=By(mi)/(b'*b);
            break;
        end
    end

    opts=0;
    B_h=B(:,a);
    x_h=x(a);
    theta_h=theta(a);

    while opts==0,
        opts=1;

        if size(B_h,2)<=length(y),
            BB=B_h'*B_h;
            x_new=BB\(B_h'*y-lambda*theta_h/2);
            o_new=L1_cost(y,B_h,x_new,lambda);

            %cost based on changing sign
            s=find(sign(x_new)~=theta_h);
            x_min=x_new;
            o_min=o_new;
            for j=1:length(s),
                zd=s(j);
                x_s=x_h-x_h(zd)*(x_new-x_h)/(x_new(zd)-x_h(zd));
                x_s(zd)=0;  %make sure it's zero
                o_s=L1_cost(y,B_h,x_s,lambda);
                if o_s<o_min,
                    x_min=x_s;
                    o_min=o_s;
                end
            end
        else
            d=x_h-B_h'*((B_h*B_h')\(B_h*x_h));
            q=x_h./(d+eps);
            x_min=x_h;
            o_min=L1_cost(y,B_h,x_h,lambda);
            for j=1:length(q),
                zd=q(j);
                x_s=x_h-zd*d;
                x_s(j)=0;       %make sure it's zero
                o_s=L1_cost(y,B_h,x_s,lambda);
                if o_s<o_min,
                   x_min=x_s;
                   o_min=o_s;
                end
            end
        end

        x(a)=x_min;

        a=(x~=0);
        theta=sign(x);

        B_h=B(:,a);
        x_h=x(a);
        theta_h=theta(a);
        Bx_h=B_h*x_h;

        active_d=2*(B_h'*(Bx_h-y))+lambda*theta_h;

        if ~isempty(find(abs(active_d)>OptTol)),
            opts=0;
        end
    end

    all_d=2*(B'*Bx_h-By);

    [ma mi]=max(abs(all_d).*(~a));
    if ma>lambda+OptTol,
        optc=0;
    end
end

return;

function cost=L1_cost(y,B,x,lambda)
    tmp = y-B*x;
    cost = tmp'*tmp+lambda*norm(x,1);
return

author: thomas-koehler

repo: SupER


text:
function [Bt]=GetCentrifugalMatrix_2(T,Pcii,Icii,mcii,dq)
%% About the function: this is a function that is used to calculate
% Centrifugal matrix of the serially linked manipulator,  the return value of
% the function is (nxn) Centrifugal matrix from which when multiplied by
% the joints' velocity vector the torques of the centrifugal forces
% is going to be calculated
% The input parameters are as the following:
% T is (4x4xn) transformation matrix of the serially linked robot, each
% (4x4) matrix represents the transform for each link in the base frame.
% Pcii is 3Xn matrix while each column represents the local coordinates
% of the center of mass of each link.
% Icii is (3x3xn) matrix, each 3x3 matrix of which represnets the
% associated link inertial tensor represented in its local inertial frame
% mcii is (1xn) vector, each element of which specifies a mass of one of
% the links

% Copyright Mohammad SAFEEA 2nd,April,2018

n=max(size(mcii));
%% Initialization of Ai and Bi.
Bi=zeros(3,n,n);
Di=zeros(3,n,n);
Kj=zeros(3,n);
half_Kj=zeros(3,n);
%% Calculate === some auxuliary variables
Pcii_A=zeros(3,n);
mcii_Pcii_A=zeros(3,n);
Pcii_A(:,1)=T(1:3,1:3,1)*Pcii(:,1);
mcii_Pcii_A(:,1)=mcii(1)*Pcii_A(:,1);
Kj(:,1)=T(1:3,3,1);
half_Kj(:,1)=0.5*Kj(:,1);
for i=2:n
        Pcii_A(:,i)=T(1:3,1:3,i)*Pcii(:,i);
        mcii_Pcii_A(:,i)=mcii(i)*Pcii_A(:,i);
        Kj(:,i)=T(1:3,3,i);
        half_Kj(:,i)=0.5*Kj(:,i);
end
%% calculating the links model, Mci and ddPci
for i=1:n
    %% calculating the Mci term
    Pci=Pcii_A(:,i)+T(1:3,4,i);
    L=T(1:3,1:3,i)*(trace(Icii(:,:,i))*eye(3)-2*Icii(:,:,i))*T(1:3,1:3,i)';
    for j=1:i
        %Calculating inertial moment due to normal acceleration due to
        %effect of each frame j.
        Bi(:,j,i)=cross(L*half_Kj(:,j),Kj(:,j));
        %Calculating acceleration of center of mass of link i due to
        %injection of each frame j
        Pcij=Pci-T(1:3,4,j);
        Di(:,j,i)=Kj(:,j)*(Kj(:,j)'*Pcij)-Pcij;
    end
end
Bt=zeros(n,n);
Fac_D=zeros(3,n);
Mac_B=zeros(3,n);
Pjp1_j=zeros(3,1);
%% calculating Mac for all of the links, then calculating two by filling At
%% and Bt

start=n-1;
j=n;
%% recursive rprocedure on moments and forces
        for k=1:n %% iterate through the matrix
            Mac_B(:,k)=Bi(:,k,j)+cross(mcii_Pcii_A(:,j),Di(:,k,j));
        end
        %% on forces
        Fac_D=mcii(j)*Di(:,:,j);
        Bt(j,:)=T(1:3,3,j)'*Mac_B;

for j=start:-1:1 %% iterate through the joints
    Pjp1_j=T(1:3,4,j+1)-T(1:3,4,j);
    %% recursive rprocedure on moments and forces
        for k=1:j %% iterate through the matrix
            Mac_B(:,k)=Mac_B(:,k)+Bi(:,k,j)+cross(Pjp1_j,Fac_D(:,k))+cross(mcii_Pcii_A(:,j),Di(:,k,j));
        end
        for k=j+1:n
            Mac_B(:,k)=Mac_B(:,k)+cross(Pjp1_j,Fac_D(:,k));
        end
    %% on forces
    Fac_D(:,1:j)=Fac_D(:,1:j)+mcii(j)*Di(:,1:j,j);
    Bt(j,:)=T(1:3,3,j)'*Mac_B;

end
%% close the function
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function r8_cosd_test ( )

%*****************************************************************************80
%
%% R8_COSD_TEST tests R8_COSD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_COSD_TEST\n' );
  fprintf ( 1, '  R8_COSD computes the cosine of an angle\n' );
  fprintf ( 1, '  given in degrees.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ANGLE    R8_COSD(ANGLE)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 15 : 360
    angle = i;
    fprintf ( 1, '  %8.2f  %14.6g\n', angle, r8_cosd ( angle ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% LRR | EALM | Exact ALM (Lin et al. 2009)
% process_video('LRR', 'EALM', 'dataset/demo.avi', 'output/demo_LRR-EALM.avi');

alg_path_aux = fullfile(lrs_conf.lrr_path,'ALM');
addpath(genpath(alg_path_aux));

A = mean(M,2);
lambda = 0.01;
[Z,E] = solve_lrr(M,A,lambda,0,0,1);
% M_hat = A*Z + E;
L = A*Z;
S = E;

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function area = sphere_cap_area_3d ( r, h )

%*****************************************************************************80
%
%% SPHERE_CAP_AREA_3D computes the surface area of a spherical cap in 3D.
%
%  Discussion:
%
%    Draw any radius of the sphere and note the point P where the radius
%    intersects the sphere.  Consider the point on the radius line which is
%    H units from P.  Draw the circle that lies in the plane perpendicular to
%    the radius, and which intersects the sphere.  The circle divides the sphere
%    into two pieces, and the corresponding disk divides the solid sphere into
%    two pieces.  The spherical cap is the part of the solid sphere that
%    includes the point P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real H, the "height" of the spherical cap.
%
%    Output, real AREA, the area of the spherical cap.
%
  if ( h <= 0.0 )
    area = 0.0;
  elseif ( 2.0 * r <= h )
    area = 4.0 * pi * r * r;
  else
    area = 2.0 * pi * r * h;
  end

  return
end


text:
%Voice Based Biometric System
%By Ambavi K. Patel.

function disply(nmatch)
%indicate result
mx=max(max(nmatch));

   if  mx<1
             h1 = msgbox('INCORRECT','your password is','error');
       soundsc(wavread('denied.wav'),50000);

    uiwait(h1,1);
close(h1);

   else
      h1 = msgbox('CORRECT','your password is');
      soundsc(wavread('allow.wav'),50000);

     uiwait(h1,1);
close(h1);
   end;

author: Sable

repo: mcbench-benchmarks


text:
% this function compute the spatial covariance matrix from complex Fourier
% transform of array signals.
% Inputs:
%   X: D x N x T matrix of Fourier transform coefficients. D is the number
%   of frequency bins, N is the number of microphone channels, and T is the
%   number of frames
%   context_size: number of context frames used to compute the spatial
%   covariance at each time-frequency bin.
% Outputs:
%   R: N x N x D x T matrix of spatial covariance matrixes
%
% Author: Xiong Xiao, Nanyang Technological University, Singapore
% Last Modified: 22 Feb 2016.
%
function R = ComplexSpectrum2SpatialCov(X, context_size, shift, useGPU)
if nargin<4
    useGPU = 1;
end

half_ctx = (context_size-1)/2;
[D, N, T] = size(X);    % D is number of freq bin, N is number of channel

if context_size == 0    % get global spatial covariance matrix
    X2 = permute(X, [2 1 3]);
    XX = outProdND(X2);
    R = mean(XX,4);
else    % get windowed spatial covariance
    if 0
        nBlock = length(1:shift:T);
        R = zeros(N,N,D,nBlock);
        for d = 1:D
            for t = 1:shift:T
                frame_selector = t-half_ctx:t+half_ctx;
                frame_selector = min(T,max(1,frame_selector));
                X2 = squeeze(X(d,:,frame_selector));
                R(:,:,d,(t-1)/shift+1) = X2*X2' / context_size;
            end
        end
    else
        X2 = permute(X, [2 1 3]);
        XX = outProdND(X2);
%         X2cell = num2cell(X2, [1]);
%         XXcell = cellfun(@(x) (reshape(x*x', N^2,1)), X2cell, 'UniformOutput', 0);
%         XX = cell2mat(XXcell);
        XX2 = reshape(XX, N^2*D, T);
        if 1
            if useGPU
                XX2 = gpuArray(XX2);
            end
            idx = [ones(1,half_ctx) 1:T ones(1,half_ctx)*T];
            SCM = conv2(XX2(:,idx), ones(1,context_size, class(gather(X)))/context_size, 'valid');
%             SCM = SCM(:,half_ctx+1:end-half_ctx);
        else
            fake_layer.a = XX2;
            XX3 = F_splice(fake_layer, context_size);
            XX4 = reshape(XX3,  N^2*D, context_size, T);
            SCM = mean(XX4,2);
        end
        SCM2 = reshape(SCM, N^2, D, T);
        SCM3 = reshape(SCM2, N, N, D, T);
        R = SCM3(:,:,:,1:shift:end);
    end
end

author: singaxiong

repo: SignalGraph


text:
function [out, idxtrain,eigvector, eigvalue] = kpca(spectraldata, No_Train, dimension, kernel, parameter)

%
%  Kernel principal component analysis, KPCA
%
% [out, eigvector, eigvalue] = kpca(spectraldata, No_Train, dimension, kernel, parameter)
%
% Input
%
% spectraldata       - input hyperspectral data with 3-D, ncols by nrows by nbands
% No_Train           - randomly selected samples for training (<= 5000),
%                       depending on the memory of your PC
% dimension          - the dimension of extracted kernelized PCs (with largest eigen values)
% kernel             - kernel function
% parameter          - parameters for kernel fuction
%                        - kernel = 'linear';% linear kernel
%                        - kernel = 'Gaussian'; parameter=1; %Gausian kernel
%                        - kernel = 'poly'; parameter=[1,3];% third order polynomial
%
%
% Output
%
% out                - the extracted kernelized PCs
% eigvector          - the eigenvectors divided by square root of corresponding eigenvalues
% eigvalue           - the first dimension largest eigenvalues

if nargin <  3, error('not enough input'); end
if nargin <  4
    if strncmp(kernel,'Gaussian',1)
        parameter=1;
    elseif strncmp(kernel,'poly',1)
        parameter=[1, 3];
    end
end

[nrows,ncols,nbands] = size(spectraldata);
X0 = reshape(spectraldata,nrows*ncols,nbands);
clear spectraldata

%% sub-sample for training, select No_Train samples randomly
rand('state',4711007);% initialization of rand
if No_Train>nrows*ncols, No_Train = nrows*ncols; end
idxtrain = randsample(nrows*ncols, No_Train);
X = double(X0(idxtrain,:));
ntrain = size(X,1);

Xtest = X0;
ntest = size(Xtest,1);
clear X0;

%% kernelized training data and centering the kernel matrix
[K scale sums] = kernelize_training(kernel, X, parameter);
meanK = mean(K(:));
meanrowsK = mean(K);
K = centering(K);

%% select the first dimension eigvectors
dimout = ntrain;
if dimout>dimension, dimout=dimension; end
[eigvector,eigvalue,flagk] = eigs(K, dimout, 'LM');

if flagk~=0, warning('*** Convergence problems in eigs ***'); end
eigvalue = diag(abs(eigvalue))';
eigvector = sqrt(ntrain-1)*eigvector*diag(1./sqrt(eigvalue));
clear K

out = NaN(ntest,dimout);
%% kernelized the test samples, center kernel and calculate kernel PCs
for rr=1:nrows
     idx = (rr-1)*ncols+1;
     idx = idx:(idx+ncols-1);
     Xk = kernelize_test(kernel, X, Xtest(idx,:), scale, sums);
     Xk = Xk - repmat(meanrowsK,ncols,1)' - repmat(mean(Xk),ntrain,1) + meanK;
     out(idx,:) = Xk'*eigvector;
end

out = reshape(out,nrows,ncols,dimout);

author: BehnoodRasti


text:
%Save_Figure_Script

%Variables that must be defined in the workspace:
% P
% iterationNum
%

%Now, save the figure if desired:
if P.disp.saveIntermediatePlots
   BaseName = ['MS_Fig_Iter_' num2str(iterationNum)];
   FigureName = [BaseName '.pdf'];
   DirName = 'img/';  %Which subdirectory are the images in?

   % %Force Figure Location (Good for extracting figures -> LaTeX)
   set(gcf,'position',800*[0,0,2,1])

   %This is a magical function that I got on file exchange
   %This should not be so hard to do in Matlab...
   save2pdf(FigureName,gcf,400)

   %Now, create the latex code if desired:
   if P.disp.createTex
       FileName = P.disp.TexFileName;
       fid = fopen(FileName,'a');  %Open file for appending
        fprintf(fid,'\n');
        fprintf(fid,'\\begin{figure} \n');
        fprintf(fid,'    \\centering \n');
        fprintf(fid,['    \\includegraphics[width = \\columnwidth]{' DirName FigureName '}    \n']);
        fprintf(fid,'    \\caption{}     \n');
        fprintf(fid,['    \\label{fig: ' BaseName '}     \n']);
        fprintf(fid,'\\end{figure} \n');
        fclose(fid);
   end
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [out, bestL] = Peubseparate(k, d, PeCh, e)
%Excess-distortion probability for the transmission of a Gaussian source
%over a Gaussian channel as achieved by a separated scheme - optimized with respect to the number of
%messages allowable between the channel encoder and the channel decoder (2^L)

%PeCh - array of channel error probabilities
%k - source blocklength
%d - target distortion
%e - excess distortion probability

%outputs:
%out - the excess distortion probability
%bestL - the optimal length of messages exchanged between the encoder and
%the decoder

%
%   Created in 2012 by Victoria Kostina (vkostina@caltech.edu)
%

Kbig = 20;
tol = 1e-3;

out = Inf;
bestL = NaN;
for L = find(PeCh <= e)
    cur = PeubSource(L, k, d) + PeCh(L);
    if out > cur
        out = cur;
        bestL = L;
    end
end
%fprintf('ASeparate: k = %i, Peub = %f\n', k, out);

    function out = PeubSource(L, k, d)
        %achievability via sphere covering

        r0 = sqrt( max( 0, 1 - d) );
        a = r0 - sqrt(d);
        b = r0 + sqrt(d);

        out = chi2cdf((a^2 + tol)*k, k) + 1 - chi2cdf((b^2 - tol)*k, k)...  %nontypical source realization
            + quad(@Pes, a^2 + tol, b^2- tol);                              %source error

        function out = logPdball(x)
            %x - distance squared from the origin (x*k central chi square k)
            cosa = (r0^2 + x - d)./(2*r0.*sqrt(x));
            sina = (1 - cosa.^2).^(1/2);
            out = loggamma(k/2+1, 'LB') - loggamma((k-1)/2 +1, 'UB')- .5*log(pi)-log(k) + (k-1)*log(sina);
        end

        function out = Pes(x)
            %source error conditioned on the distance squared from the origin (x*k central chi square k)
            if k <  Kbig
                out = (1 - exp(logPdball(x))).^exp(L);
            else
                out = exp( - exp( logPdball(x) + L));
            end
            out = out.*chi2pdf(k*x, k).*k;
        end
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
% demo_example_lmgp_data - demonstration script file for LMGP model identification.
%
%% Description
% Demo to present how to obtain and compose data for the identification
% of GP model with incorporated local models (LMGP model).
%
% There are two options how to obtain local model data:
% a. analytically (set fIdentifyLm to 0)
% b. with identification (set fIdentifyLm to 1)
%
% It can be used only with the Gaussian covariance function and
% with the white noise model (sum of covSEard and covNoise).
%
% See Also
% example.m, demo_example_LMGP_training.m, demo_example_LMGP_simulation.m,
% example_LM_ident.m, example_derivative.m
%
%%
% Written by K. Azman, 2007
%
% Changelog:
%
% 16.2.2015, Martin Stepancic:
%		 	-included the gathering of both equilibrium and non-equilibrium
%			 points in the same matlab script;
%			-changed .mat filenames and their structure contents to:
%			 {train|valid}_data.{x|u|y}
%
clear;
global flag_LM_data_ident
% LM generation data
flag_LM_data_ident = 0; % 0=calculated, 1=identified using IV

me = mfilename;                                            % what is my filename
mydir = which(me); mydir = mydir(1:end-2-numel(me));        % where am I located
addpath([mydir 'system']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                 Points at Equilibrium                              %%%%

% local model ident
ord = 1;    % order of identified LMs
eq.U = [-1.5 -1.1 -0.7 -0.3 0 0.3 0.7 1.1 1.5]';
dU = 0.1;  % depending on diference between the elements of eq.U, valid only if LMs identified!

noiseStd = 0.1;
% lm calculated/identified
if (flag_LM_data_ident)
    % *** local model's parameters identified and not calculated ***

    for ii = 1:length(eq.U)
        [lm{ii},y0] = demo_example_LM_ident(eq.U(ii), dU, noiseStd, 2000+ii);

        eq.Y(ii,1) = y0;
        eq.dfdy(ii,1) = -lm{ii}.a(2);
        eq.dfdu(ii,1) = -lm{ii}.b(2);
    end
    eq.lm = lm;

    filename = 'example_data_lmgp_eq_ident';

else
    % *** local model's parameters calculated ***
    [eq.dfdy, eq.dfdu, eq.Y] = demo_example_derivative(eq.U);

    filename = 'example_data_lmgp_eq_anal';

end

figure('Name', 'EQ Curve');
plot(eq.U,eq.Y);

save(filename,'eq');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                 Non-eqilibrium points                              %%%%
umax = 1.5;
umin = -umax;
N_training_points = 200;
N_validation_points = 200;
noise_std = 0.05;

utrain = sig_prs_minmax(N_training_points,2,umin,umax); % create input singal

% simulate example system, to obtain delayed outputs and delayed outputs
[ytrain, xtrain, utrain] = demo_example(utrain);
% add noise to training data
[ytrain, ytrain_no_noise] = add_noise_to_vector(ytrain, noise_std);
figure('Name', 'Training Input Signal');
t=1:length(utrain);
plot(t, utrain, t, xtrain);
legend('u', 'y');
grid on;

% validation data --> less excited signal
uvalid = sig_prs_minmax(N_validation_points,5,umin,umax);
[yvalid, xvalid, uvalid] = demo_example(uvalid);

figure('Name', 'Validation Input Signal');
plot(uvalid);

% training data
train_data.u = utrain;
train_data.x = xtrain;
train_data.y = ytrain;

% validation data
valid_data.u = uvalid;
valid_data.x = xvalid;
valid_data.y = yvalid;

save example_data_lmgp_oeq valid_data train_data umax umin noise_std
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

save(filename,'eq')
rmpath([mydir 'system'])
return

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function value = i4_factorial2 ( n )

%*****************************************************************************80
%
%% I4_FACTORIAL2 computes the double factorial N!!
%
%  Discussion:
%
%    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
%                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the double factorial function.
%    If N is less than 1, the value is returned as 1.
%
%    Output, integer VALUE, the value of N!!.
%
  if ( n < 1 )
    value = 1;
    return
  end

  value = 1;

  while ( 1 < n )
    value = value * n;
    n = n - 2;
  end

  return
end


text:
function[result]=Forward(docbatch,parameter,isTraining)
    [result]=Forward_Source_Word(docbatch,parameter,isTraining);
    % forward calculation at word level for each sentence
    [result]=Forward_Source_Sen(result,docbatch,parameter,isTraining);
    % forward calculation at sentence level for each document
    if isTraining==1
        [result]=Forward_Target(result,docbatch,parameter,isTraining);
    end
end

function[result]=Forward_Source_Word(docbatch,parameter,isTraining)
    % forward calculation at word level for each sentence
    sourceBatch=docbatch.source_smallBatch;
    result.source_sen_vector=[];
    for i=1:length(sourceBatch)
        batch=sourceBatch{i};
        T=batch.max_length;
        h_t_source_word=cell(parameter.layer_num,T);
        % store h_t
        result.c_t_source_word{i}=cell(parameter.layer_num,T);
        % store c_t
        result.lstms_source_word{i} = cell(parameter.layer_num,T);
        % store gate values for lstm units
        N=size(batch.Word,1);
        zeroState=zeroMatrix([parameter.hidden,N]);
        for ll=1:parameter.layer_num
            for tt=1:T
                h_t_source_word{ll,tt}=zeroState;
                result.c_t_source_word{i}{ll,tt}=zeroState;
            end
        end
        for t=1:T
            for ll=1:parameter.layer_num
                W=parameter.Word_S{ll};
                % W for word level composition in source
                if t==1
                    h_t_1=zeroState;
                    c_t_1=zeroState;
                else
                    c_t_1=result.c_t_source_word{i}{ll, t-1};
                    h_t_1=h_t_source_word{ll, t-1};
                end
                if ll==1
                    x_t=parameter.vect(:,batch.Word(:,t));
                else
                    x_t=h_t_source_word{ll-1,t};
                end
                x_t(:,batch.Delete{t})=0;
                h_t_1(:,batch.Delete{t})=0;
                c_t_1(:,batch.Delete{t})=0;
                % set postion that do not have words to zero
                [result.lstms_source_word{i}{ll, t},h_t_source_word{ll, t},result.c_t_source_word{i}{ll, t}]=lstmUnit(W,parameter,x_t,[],h_t_1,c_t_1,ll,t,isTraining);
                % compute lstm unit
                if t==T && ll==parameter.layer_num
                    result.source_sen_vector=[result.source_sen_vector,h_t_source_word{ll,t}];
                    % store vector embeddings for each source sentence
                    clear h_t_source_word;
                end
            end
        end
    end
    clear x_t;
    clear h_t_1;
    clear c_t_1;
end

function[result]=Forward_Source_Sen(result,docbatch,parameter,isTraining)
    % forward calculation at sentence level for each document
    T=docbatch.max_source_sen;
    h_t_source_sen=cell(parameter.layer_num,T);
    result.c_t_source_sen=cell(parameter.layer_num,T);
    result.lstms_source_sen=cell(parameter.layer_num,T);
    result.source_sen=cell(parameter.layer_num,1);
    result.source_each_sen=gpuArray();

    N=size(docbatch.source_sen_matrix,1);
    zeroState=zeroMatrix([parameter.hidden,N]);
    for ll=1:parameter.layer_num
        for tt=1:T
            h_t_source_sen{ll,tt}=zeroState;
            result.c_t_source_sen{ll,tt}=zeroState;
        end
    end
    for t=1:T
        for ll=1:parameter.layer_num
            W=parameter.Sen_S{ll};
            % sentence-level compostion
            if t==1
                h_t_1=zeroState;
                c_t_1 =zeroState;
            else
                c_t_1 =result.c_t_source_sen{ll, t-1};
                h_t_1 =h_t_source_sen{ll, t-1};
            end
            if ll==1
                x_t=result.source_sen_vector(:,docbatch.source_sen_matrix(:,t));
                % compute sentence-level embedding
            else
                x_t=h_t_source_sen{ll-1,t};
            end
            x_t(:,docbatch.source_delete{t})=0;
            h_t_1(:,docbatch.source_delete{t})=0;
            c_t_1(:,docbatch.source_delete{t})=0;
            % set values for deleted postion to 0
            [result.lstms_source_sen{ll, t},h_t_source_sen{ll, t},result.c_t_source_sen{ll, t}]=lstmUnit(W,parameter,x_t,[],h_t_1,c_t_1,ll,t,isTraining);
            % lstm unit calculation at sentence level
            if t==T
                result.source_sen{ll,1}=h_t_source_sen{ll, t};
                % store vector embeddings for documents
            end
            if ll==parameter.layer_num
                result.source_each_sen=[result.source_each_sen,h_t_source_sen{parameter.layer_num,t}];
            end
        end
    end
    clear result.source_sen_vector;
    clear h_t_source_sen;
    clear x_t;
    clear h_t_1;
    clear c_t_1;
end

function[result]=Forward_Target(result,docbatch,parameter,isTraining)
    % forward for target documents
    T=docbatch.max_target_sen;
    result.h_t_target_sen=cell(parameter.layer_num,T);
    % sentence level embeddings
    result.c_t_target_sen=cell(parameter.layer_num,T);
    result.lstms_target_sen=cell(parameter.layer_num,T);
    N=size(docbatch.target_sen_matrix,1);
    zeroState=zeroMatrix([parameter.hidden,N]);

    for ll=1:parameter.layer_num
        for tt=1:T
            result.h_t_target_sen{ll,tt}=zeroState;
            result.c_t_target_sen{ll,tt}=zeroState;
        end
    end
    result.Target_sen={};
    for sen_tt=1:T
        for ll=1:parameter.layer_num
            W=parameter.Sen_T{ll};
            % sentnece compositions for target
            if sen_tt==1
                % if sentence index is 1, h_t_1 and c_t_1 are from outputs from source sentences
                h_t_1=result.source_sen{ll,1};
                dim=size(result.c_t_source_sen);
                c_t_1=result.c_t_source_sen{ll,dim(2)};
            else
                % otherwise, c_t_1 are lstm outputs from last time step
                c_t_1 =result.c_t_target_sen{ll, sen_tt-1};
                h_t_1 =result.h_t_target_sen{ll, sen_tt-1};
            end
            if ll==1
                Word_List=docbatch.target_word{sen_tt}.Word;
                Word_Delete=docbatch.target_word{sen_tt}.Delete;
                if sen_tt==1
                    M1=result.source_sen(:,1);
                    dim=size(result.c_t_source_sen);
                    M2=result.c_t_source_sen(:,dim(2));
                else
                    M1=result.h_t_target_sen(:,sen_tt-1);
                    M2=result.c_t_target_sen(:,sen_tt-1);
                end
                result=Forward_Target_Word(result,M1,M2,Word_List,Word_Delete,docbatch,parameter,isTraining,sen_tt);
                x_t=result.Target_sen{sen_tt}.h_t_target_word{parameter.layer_num,size(Word_List,2)};
                result=Attention(docbatch,result,M1{parameter.layer_num,1},parameter,sen_tt);
                % compute attention values
                m_t=result.sum_vector{sen_tt};
            else
                x_t=result.h_t_target_sen{ll-1,sen_tt};
                m_t=[];
            end
            x_t(:,docbatch.target_delete{sen_tt})=0;
            h_t_1(:,docbatch.target_delete{sen_tt})=0;
            c_t_1(:,docbatch.target_delete{sen_tt})=0;
            % set deleted postions to 0 value
            [result.lstms_target_sen{ll,sen_tt},result.h_t_target_sen{ll,sen_tt},result.c_t_target_sen{ll,sen_tt}]=lstmUnit(W,parameter,x_t,m_t,h_t_1,c_t_1,ll,sen_tt,isTraining);
            % lstm unit calculation
        end
    end
    clear target_sen_vector;
    clear x_t;
    clear h_t_1;
    clear c_t_1;
end

function[result]=Forward_Target_Word(result,h_t_sen,c_t_sen,Word_List,Word_Delete,docbatch,parameter,isTraining,sen_tt)
    % obtain sentence level embeddings for target sentences
    N=size(Word_List,1);
    T=size(Word_List,2);
    target_sen.h_t_target_word=cell(parameter.layer_num,T);
    target_sen.c_t_target_word=cell(parameter.layer_num,T);
    target_sen.lstms=cell(parameter.layer_num,T);
    zeroState=zeroMatrix([parameter.hidden,N]);

    for ll=1:parameter.layer_num
        for tt=1:T
            target_sen.h_t_target_word{ll,tt}=zeroState;
            target_sen.c_t_target_word{ll,tt}=zeroState;
        end
    end
    for t=1:T
        for ll=1:parameter.layer_num
            W=parameter.Word_T{ll};
            if t==1
                h_t_1=h_t_sen{ll,1};
                c_t_1=c_t_sen{ll,1};
            else
                c_t_1 =target_sen.c_t_target_word{ll, t-1};
                h_t_1 =target_sen.h_t_target_word{ll, t-1};
            end
            if ll==1
                x_t=parameter.vect(:,Word_List(:,t));
                % if ll=1, x_t are correspondent word embeddings
            else
                x_t=target_sen.h_t_target_word{ll-1,t};
                % else x_t are outputs from previous layer
            end
            x_t(:,Word_Delete{t})=0;
            h_t_1(:,Word_Delete{t})=0;
            c_t_1(:,Word_Delete{t})=0;
            % set deleted positions to 0
            [target_sen.lstms{ll, t},target_sen.h_t_target_word{ll, t},target_sen.c_t_target_word{ll, t}]=lstmUnit(W,parameter,x_t,[],h_t_1,c_t_1,ll,t,isTraining);
            % lstm unit calculations
            if t~=1
                target_sen.h_t_target_word{ll, t}(:,Word_Delete{t})=target_sen.h_t_target_word{ll,t-1}(:,Word_Delete{t});
                 % sentence representations stay the same for deleted positions
            end
        end
    end
    result.Target_sen{sen_tt}=target_sen;
    clear h_t_1;
    clear c_t_1;
    clear x_t;
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
filename='Bridge_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% setdiff_bc - setdiff backward compatible with Matlab versions prior to 2013a

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

function [C,IA] = setdiff_bc(A,B,varargin)

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end
end

if v >= 7.14
    [C,IA] = setdiff(A,B,varargin{:},'legacy');
    if errorFlag
        [C2,IA2] = setdiff(A,B,varargin{:});
        if (~isequal(C, C2) || ~isequal(IA, IA2))
            warning('backward compatibility issue with call to setdiff function');
        end
    end
else
    [C,IA] = setdiff(A,B,varargin{:});
end


text:
function [f,J,Q] = spm_fx_gen(x,u,P,M)
% generic state equations for a neural mass models
% FORMAT [f,J,D] = spm_fx_gen(x,u,P,M)
% FORMAT [f,J]   = spm_fx_gen(x,u,P,M)
% FORMAT [f]     = spm_fx_gen(x,u,P,M)
% x  - neuronal states
% u  - exogenous input
% P  - model parameters
% M  - model structure
%
% This routine compiles equations of motion for multiple nodes or neural
% masses in the cell array of hidden states. To include a new sort of node,
% it is necessary to update the following routines:
%
% spm_dcm_neural_priors: to specify the intrinsic parameters of a new model
% spm_dcm_x_neural:      to specify its initial states
% spm_L_priors:          to specify which hidden states generate signal
% spm_fx_gen (below):    to specify how different models interconnect
%
% This routine deal separately with the coupling between nodes (that depend
% upon extrinsic connectivity, sigmoid activation functions and delays -
% and coupling within nodes (that calls on the model specific equations of
% motion).
%
% In generic schemes one can mix and match different types of sources;
% furthermore, they can have different condition-specific modulation of
% intrinsic connectivity parameters and different, source-specific-
% contribution to the lead field (or electrode gain). Source-specific
% models are specified by a structure array model, For the i-th source:
%
% model(i).source  = 'ECD','CMC',...  % source model
% model(i).B       = [i j k ,...]     % free parameters that have B effects
% model(i).J       = [i j k ,...]     % cardinal states contributing to L
% model(i).K       = [i j k ,...]     % other states contributing to L
% ...
%__________________________________________________________________________
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_gen.m 7409 2018-08-27 11:39:00Z bernadette $

% model-specific parameters
%==========================================================================

% model or node-specific state equations of motions
%--------------------------------------------------------------------------
fx{1} = @spm_fx_erp;                     % ERP model
fx{2} = @spm_fx_cmc;                     % CMC model
fx{3} = @spm_fx_bgt;                     % basal ganglia circuit
fx{4} = @spm_fx_mmc;                     % motor micro circuit

% indices of extrinsically coupled hidden states
%--------------------------------------------------------------------------
efferent(1,:) = [9 9 9 9];               % sources of ERP connections
afferent(1,:) = [4 8 5 8];               % targets of ERP connections

efferent(2,:) = [3 3 7 7];               % sources of CMC connections
afferent(2,:) = [2 8 4 6];               % targets of CMC connections

efferent(3,:) = [9 9 9 9];               % sources of BGT connections (thalamus)
afferent(3,:) = [2 6 2 6];               % targets of BGT connections (striatum & STN)

efferent(4,:) = [3 3 7 7];               % sources of MMC connections
afferent(4,:) = [2 4 8 0];               % targets of MMC connections

% scaling of afferent extrinsic connectivity (Hz)
%--------------------------------------------------------------------------
E(1,:) = [64 2 32 8]*1000;               % ERP connections
E(2,:) = [64 4 64 8]*1000;               % CMC connections
E(3,:) = [1.8 1.2 1.8 1.2]*100000;       % BGC connections
E(4,:) = [.9 .9 .11 0]*100000;           % MMC connections

if isfield(M,'ERP_E'); E(1,:)= M.ERP_E; end
if isfield(M,'CMC_E'); E(2,:)= M.CMC_E; end
if isfield(M,'BGT_E'); E(3,:)= M.BGT_E; end
if isfield(M,'MMC_E'); E(4,:)= M.MMC_E; end

% intrinsic delays ms (scaling)
%--------------------------------------------------------------------------
D(1) = 2;                                % ERP connections
D(2) = 1;                                % CMC connections
D(3) = 4;                                % BGT connections
D(4) = 1;                                % MMC connections

% get model specific operators
%==========================================================================
if isvector(x)
    x = spm_unvec(x,M.x);                               % neuronal states
end

% get the neural mass models {'ERP','CMC'}
%--------------------------------------------------------------------------
n     = numel(x);
model = M.dipfit.model;
for i = 1:n
    if     strcmp(model(i).source,'ERP')
        nmm(i) = 1;
    elseif strcmp(model(i).source,'CMC')
        nmm(i) = 2;
    elseif strcmp(model(i).source,'BGT')
        nmm(i) = 3;
    elseif strcmp(model(i).source,'MMC')
        nmm(i) = 4;
    end
end

% synaptic activation function
%--------------------------------------------------------------------------
R     = 2/3;                      % gain of sigmoid activation function
B     = 0;                        % bias or background (sigmoid)
R     = R.*exp(P.S);              % posterior gain
S     = @(x,R,B)1./(1 + exp(-R*x(:) + B)) - 1/(1 + exp(B));
dSdx  = @(x,R,B)(R*exp(B - R*x(:)))./(exp(B - R*x(:)) + 1).^2;
for i = 1:n
    Sx{i} = S(x{i},R,B);          % sigmod firing rate function
end

% Extrinsic connections
%--------------------------------------------------------------------------
for i = 1:numel(P.A)
    A{i} = exp(P.A{i});
end

% detect and reduce the strength of reciprocal (lateral) connections
%--------------------------------------------------------------------------
TOL   = exp(-8);
for i = 1:numel(A)
    L    = (A{i} > TOL) & (A{i}' > TOL);
    A{i} = A{i}./(1 + 4*L);
end

% and scale of extrinsic connectivity (Hz)
%--------------------------------------------------------------------------
for i = 1:n
    for k = 1:numel(P.A)
        A{k}(i,:) = E(nmm(i),k)*A{k}(i,:);
    end
end

% assemble flow
%==========================================================================
N     = M;
for i = 1:n

    % intrinsic flow
    %----------------------------------------------------------------------
    N.x  = M.x{i};
    if isfield(M,'u')
        ui = u(i,:);
    else
        ui = u;
    end
    Q    = P.int{i};
    Q.C  = P.C(i,:);
    f{i} = fx{nmm(i)}(x{i},ui,Q,N);

    % extrinsic flow
    %----------------------------------------------------------------------
    for j = 1:n
        for k = 1:numel(P.A)
            if A{k}(i,j) > TOL
                ik       = afferent(nmm(i),k);
                jk       = efferent(nmm(j),k);
                f{i}(ik) = f{i}(ik) + A{k}(i,j)*Sx{j}(jk);
            end
        end
    end
end

% concatenate flow
%--------------------------------------------------------------------------
f  = spm_vec(f);

if nargout < 2; return, end

% Jacobian
%==========================================================================
for i = 1:n
    for j = 1:n

        J{i,j} = zeros(numel(x{i}),numel(x{j}));
        if i == j

            % intrinsic flow
            %--------------------------------------------------------------
            N.x    = M.x{i};
            if isfield(M,'u')
                ui = u(i,:);
            else
                ui = u;
            end
            Q      = P.int{i};
            Q.C    = P.C(i,:);
            J{i,i} = spm_diff(fx{nmm(i)},x{i},ui,Q,N,1);

        else

            % extrinsic flow
            %--------------------------------------------------------------
            for k = 1:numel(P.A)
                if A{k}(i,j) > TOL
                    ik  = afferent(nmm(i),k);
                    jk  = efferent(nmm(j),k);
                    dS  = dSdx(x{j}(jk),R,B);
                    J{i,j}(ik,jk) = J{i,j}(ik,jk) + A{k}(i,j)*dS;
                end
            end
        end
    end
end

J  = spm_cat(J);

if nargout < 3; return, end

% delays
%==========================================================================
% Delay differential equations can be integrated efficiently (but
% approximately) by absorbing the delay operator into the Jacobian
%
%    dx(t)/dt     = f(x(t - d))
%                 = Q(d)f(x(t))
%
%    J(d)         = Q(d)df/dx
%--------------------------------------------------------------------------
% Implement: dx(t)/dt = f(x(t - d)) = inv(1 + D.*dfdx)*f(x(t))
%                     = Q*f = Q*J*x(t)
%--------------------------------------------------------------------------

% source specific (intrinsic) delays
%--------------------------------------------------------------------------
for i = 1:n
    P.D(i,i) = P.D(i,i) + log(D(nmm(i)));
end

% N-th order Taylor approximation to delay
%--------------------------------------------------------------------------
N     = 256;
for i = 1:8
    Q = spm_dcm_delay(P,M,J,N);
    if isfinite(norm(Q,'inf')) && (norm(Q*J,'inf') < norm(J,'inf')*4)
        break
    else
        N = N/2;
    end
end

return

author: spm


text:
clc; clear all; close all;
I = checkerboard(50,3,3);
h = fspecial('gaussian',[5 5],2);
harris(I,0.05,0.01,h);

author: decouples

path: github-repos/MATLAB/decouples-Matlab_deep_learning/Matlab_deep_learning-1b823b82686080e32b03e1f1a4648896bd6e3c44/第 17 章 基于 Harris 的角点特征检测/main.m


text:
function [] = checkNumericalGradient()
% This code can be used to check your numerical gradient implementation
% in computeNumericalGradient.m
% It analytically evaluates the gradient of a very simple function called
% simpleQuadraticFunction (see below) and compares the result with your numerical
% solution. Your numerical gradient implementation is incorrect if
% your numerical solution deviates too much from the analytical solution.

% Evaluate the function and gradient at x = [4; 10]; (Here, x is a 2d vector.)
x = [4; 10];
[value, grad] = simpleQuadraticFunction(x);

% Use your code to numerically compute the gradient of simpleQuadraticFunction at x.
% (The notation "@simpleQuadraticFunction" denotes a pointer to a function.)
numgrad = computeNumericalGradient(@simpleQuadraticFunction, x);

% Visually examine the two gradient computations.  The two columns
% you get should be very similar.
disp([numgrad grad]);
fprintf('The above two columns you get should be very similar.\n(Left-Your Numerical Gradient, Right-Analytical Gradient)\n\n');

% Evaluate the norm of the difference between two solutions.
% If you have a correct implementation, and assuming you used EPSILON = 0.0001
% in computeNumericalGradient.m, then diff below should be 2.1452e-12
diff = norm(numgrad-grad)/norm(numgrad+grad);
disp(diff);
fprintf('Norm of the difference between numerical and analytical gradient (should be < 1e-9)\n\n');
end

function [value,grad] = simpleQuadraticFunction(x)
% this function accepts a 2D vector as input.
% Its outputs are:
%   value: h(x1, x2) = x1^2 + 3*x1*x2
%   grad: A 2x1 vector that gives the partial derivatives of h with respect to x1 and x2
% Note that when we pass @simpleQuadraticFunction(x) to computeNumericalGradients, we're assuming
% that computeNumericalGradients will use only the first returned value of this function.

value = x(1)^2 + 3*x(1)*x(2);

grad = zeros(2, 1);
% 对X1 求偏导
grad(1)  = 2*x(1) + 3*x(2);
% 对X2 求偏导
grad(2)  = 3*x(1);

end

repo: MachineLearning


text:
function pass = test_transpose( )
% Test transpose and ctranspose

tol = 10*chebfunpref().cheb2Prefs.chebfun2eps;

% Test empty diskfunv.
f = diskfunv;
pass(1) = isempty(f');
pass(2) = isempty(f.');

% Test function
f = diskfun(@(x,y) cos((x+.1).*y));
% diskfunv
u = grad(f);
[m,n,p] = size(u');
pass(3) = all( ([isinf(m) isinf(n) p==2]) );
[m,n,p] = size(u.');
pass(4) = all( ([isinf(m) isinf(n) p==2]) );

% Check transpose and ctranspose give the same results for real-valued
% diskfunv objects.
w = u'-u.';
rng(10); th0 = rand; r0 = rand;
pass(5) = ( norm(w(th0,r0, 'polar')) < tol );

% Check transpose of transpose is the original u.
v = u.';
pass(6) = ( norm(u-v.') < tol );

% Check u'*u is a diskfun.
pass(7) = isa(u'*u,'diskfun');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% $Id: realign_synth_pca.m 2 2009-06-16 19:24:10Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-06-16 15:24:10 -0400 (Mar, 16 jui 2009) $
% $Revision: 2 $

%% build synthetic data

% synthetic data
randn('seed',0);
rand('seed',0);

clear EEG

EEG.pnts = 500;
EEG.srate = 1000; % 0.5 sec. recording
EEG.trials = 100;
EEG.nbchan = 1;

EEG.times = linspace(0,EEG.pnts / EEG.srate,EEG.pnts); % tf*100

noiselevel = 5;
[EEG,ref,latencies,snr] = generate_synth_eeg_data(EEG,noiselevel);

X = squeeze(EEG.data)';
X = X - ones(size(X,1),1)*mean(X,1);
norms = sqrt(sum(X.^2,2));

X = diag(1 ./ norms) * X;
[PC, SCORE, LATENT]=princomp(X);

[tmp,order] = sort(SCORE(:,1));

data = squeeze(EEG.data)';

%% view results
figure(400)
plot(SCORE(:,1),SCORE(:,2),'+');
xlabel('X');
ylabel('Y');
% title('Points in 2D with PCA');
% savefig('pca_embedding')

figure(401)
imagesc(data,[min(data(:)),max(data(:))])
colorbar
xlabel('Time (ms)')
ylabel('Trial')
title('Data')
set(gca,'XTick',[])

figure(402)
imagesc(data(order,:),[min(data(:)),max(data(:))])
colorbar
xlabel('Time (ms)')
ylabel('Trial')
% title('Data Ordered with PCA')
% set(gcf,'WindowButtonDownFcn','callbackRaster(EEG)'); % set the callback
set(gca,'XTick',[])
% savefig('pca_realigned')

author: fieldtrip


text:
function y = nansum(x,dim)
% nansum - Sum ignoring NaNs.
%
% Synopsis:
%   y = nansum(x)
%   y = nansum(x,dim)
%
% Arguments:
%   x: Matrix or vector
%   dim: Dimension along which sum operates. Default: First non-singleton
%       dimension.
%
% Returns:
%   y: Sum along the chosen dimension, treating NaNs as missing values.
%
% Description:
%   For vectors, nansum(X) is the sum of the non-NaN elements in
%   X. For matrices, nansum(X) is a row vector containing the sum
%   of the non-NaN elements in each column of X. For N-D arrays,
%   nansum(X) operates along the first non-singleton dimension.
%
%   nansum(X,dim) sums along the dimension dim.
%
% Examples:
%   nansum([1 2 NaN]) returns 3.
%   nansum([1 2 NaN], 1) returns [1 2 NaN].
%   nansum([1 2 NaN], 2) returns 3.
%
% See also: nanmean,nanstd
%

% Author(s): Anton Schwaighofer, Feb 2005

if nargin<2,
  % Operate along the first non-singleton dimension
  dim = min(find(size(x)~=1));
  if isempty(dim),
    dim = 1;
  end
end
% Replace NaNs with zeros.
nans = isnan(x);
x(nans) = 0;

% Protect against an entire column of NaNs
y = sum(x, dim);
allNaNs = all(nans, dim);
y(allNaNs) = NaN;

author: bbci


text:
function vw = replaceROIQuadrilateralCoords(vw,corners)
%
% vw = replaceROIQuadrilateralCoords(vw,corners)
%
% Author:  BW,AAB
% Purpose:
%   Replace the coordinates in the currently selected roi with a
%   quadrilateral region defined by the corners.

curSlice =viewGet(vw, 'Current Slice');

% is this a flat view with potential rotations/flips?
% if so, we need to rotate the corners to reflect how they would be
% displayed given the current settings
if strcmp(vw.viewType, 'Flat')
	% we need to do some reorienting to get the corners into the standard
	% view coordinate system, then rotate them:
	cornerCoords = [fliplr(corners) repmat(curSlice, [size(corners, 1) 1])]';
	cornerCoords = (rotateCoords(vw, cornerCoords, 0));
	corners = cornerCoords([2 1],:)';
end

% Create a binary image with 1s within this quadrilateral
polyIm = roipoly(vw.ui.image,corners(:,1),corners(:,2));

% % markPoly returns an image with 1's marking the polygon
% dims=size(vw.ui.image);
% polyIm = markPoly(dims);

% Compute image coordinates
polyImIndices = find(polyIm);
polyImCoords = indices2Coords(polyImIndices,size(polyIm));

% Add curSlice as 3rd row to get volume coordinates
polyCoords = [polyImCoords; curSlice*ones(1,size(polyImCoords,2))];

% Do an (inverse) rotation if necessary
if (strcmp(vw.viewType,'Flat'))
    polyCoords = (rotateCoords(vw,polyCoords,1));
end

% Convert coords to canonical frame of reference
polyCoords = curOri2CanOri(vw,polyCoords);

vw.ROIs(vw.selectedROI).coords = polyCoords;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

function [X, residual, cost] = amen( L, F, X, opts )

% set default opts
if ~exist( 'opts', 'var');        opts = struct();     end
if ~isfield( opts, 'nSweeps');    opts.nSweeps = 4;    end
if ~isfield( opts, 'maxrank');    opts.maxrank = 20;   end
if ~isfield( opts, 'maxrankRes'); opts.maxrankRes = 4; end
if ~isfield( opts, 'tolRes');     opts.tolRes = 1e-13;  end
if ~isfield( opts, 'tol');        opts.tol = 1e-13;     end

d = X.order;
n = X.size;

normF = norm(F);
cost = cost_function( L, X, F );
residual = norm( apply(L, X) - F ) / normF;

for sweep = 1:opts.nSweeps
    X = orthogonalize(X, 1);
    for mu = 1:d-1
        disp( ['Current core: ', num2str(mu)] )

        % STEP 1: Solve mu-th core opimization
        L_mu = contract( L, X, mu );
        F_mu = contract( X, F, mu );

        U_mu = L_mu \ F_mu(:);
        X.U{mu} = reshape( U_mu, size(X.U{mu}) );

        % STEP 2: Calculate current residual and cost function
        res =  F - apply(L, X);
        residual = [residual; norm( res ) / normF];
        cost = [cost; cost_function( L, X, F )];
        disp(['Rel. residual: ' num2str(residual(end)) ', Current rank: ' num2str(X.rank) ]);

        % STEP 3: Augment mu-th and (mu+1)-th core with (truncated) residual
        R = contract( X, res, [mu, mu+1] );
        R_combined = unfold(R{1},'left') * unfold(R{2},'right');

        [uu,ss,~] = svd( R_combined, 'econ');
        s = find( diag(ss) > opts.tolRes*norm(diag(ss)), 1, 'last' );
		if opts.maxrankRes ~= 0
			s = min( s, opts.maxrankRes );
		end
        R{1} = reshape( uu(:,1:s)*ss(1:s,1:s), [X.rank(mu), n(mu), s]);
        %R{2} = reshape( vv(:,1:s)', [s, n(mu+1), X.rank(mu+2)]);

        left = cat(3, X.U{mu}, R{1});
        %right = cat(1, X.U{mu+1}, R{2});

        % STEP 4: Move orthogonality to (mu+1)-th core while performing rank truncation
        [U,S,~] = svd( unfold(left,'left'), 'econ' );
        t = find( diag(S) > opts.tol*norm(diag(S)), 1, 'last' );
        t = min( t, opts.maxrank );
        X.U{mu} = reshape( U(:,1:t), [X.rank(mu), n(mu), t] );
        %X.U{mu+1} = tensorprod_ttemps( right, S(1:t,1:t)*V(:,1:t)', 1);
        X.U{mu+1} = rand( t, n(mu+1), X.rank(mu+2));

    end
    for mu = d:-1:2
        disp( ['Current core: ', num2str(mu)] )

        % STEP 1: Solve mu-th core opimization
        L_mu = contract( L, X, mu );
        F_mu = contract( X, F, mu );

        U_mu = L_mu \ F_mu(:);
        X.U{mu} = reshape( U_mu, size(X.U{mu}) );

        % STEP 2: Calculate current residual and cost function
        res =  F - apply(L, X);
        residual = [residual; norm( res ) / normF];
        disp(['Rel. residual: ' num2str(residual(end)) ', Current rank: ' num2str(X.rank) ]);
        cost = [cost; cost_function( L, X, F )];

        % STEP 3: Augment mu-th and (mu+1)-th core with (truncated) residual
        R = contract( X, res, [mu-1, mu] );
        R_combined = unfold(R{1},'left') * unfold(R{2},'right');

        [~,ss,vv] = svd( R_combined, 'econ');
        s = find( diag(ss) > opts.tolRes*norm(diag(ss)), 1, 'last' );
		if opts.maxrankRes ~= 0
			s = min( s, opts.maxrankRes );
		end
        R{2} = reshape( ss(1:s,1:s)*vv(:,1:s)', [s, n(mu), X.rank(mu+1)]);

        right = cat(1, X.U{mu}, R{2});

        % STEP 4: Move orthogonality to (mu+1)-th core while performing rank truncation
        [~,S,V] = svd( unfold(right,'right'), 'econ' );
        t = find( diag(S) > opts.tol*norm(diag(S)), 1, 'last' );
        t = min( t, opts.maxrank );
        X.U{mu} = reshape( V(:,1:t)', [t, n(mu), X.rank(mu+1)] );
        %X.U{mu+1} = tensorprod_ttemps( right, S(1:t,1:t)*V(:,1:t)', 1);
        X.U{mu-1} = rand( X.rank(mu-1), n(mu-1), t);

        %residuum = [residuum; norm( apply(L, X) - F ) / normF];

        %L_mu = contract( L, X, mu );
        %F_mu = contract( X, F, mu );

        %U_mu = L_mu \ F_mu(:);
        %X.U{mu} = reshape( U_mu, size(X.U{mu}) );
        %X = orth_at( X, mu, 'right', true );
        %residuum = [residuum; norm( apply(L, X) - F ) / normF];
        %cost = [cost; cost_function( L, X, F )];
        %disp(['Rel. residual: ' num2str(residuum(end)) ', Current rank: ' num2str(X.rank) ]);
    end

end

end

function res = cost_function( L, X, F )
res = 0.5*innerprod( X, apply(L, X) ) - innerprod( X, F );
end

function res = euclid_grad( L, X, F )
res = apply(L, X) - F;
end

author: NicolasBoumal

repo: manopt


text:
function numOsc = populateOtherSCCombo(handles, hRefSpacecraftCombo)
    maData = getappdata(handles.ma_MainGUI,'ma_data');

    oSCStr = cell(0,1);
    if(isfield(maData.spacecraft,'otherSC'))
        otherSC = maData.spacecraft.otherSC;
        if(isempty(otherSC))
            set(hRefSpacecraftCombo, 'Enable','off');
            set(hRefSpacecraftCombo,'String', 'None Available');
        else
            set(hRefSpacecraftCombo, 'Enable','on');

            for(i=1:length(otherSC)) %#ok<*NO4LP>
                oSCStr{end+1} = otherSC{i}.name; %#ok<AGROW>
            end
            set(hRefSpacecraftCombo,'String',oSCStr);
        end
    else
        set(hRefSpacecraftCombo, 'Enable','off');
        set(hRefSpacecraftCombo,'String', 'None Available');
    end

    numOsc = length(oSCStr);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [GAmp,GTime]=GzAreaTrapezoid(p) % unit m-1
% Trapezoid with prescribed area & time span & maximum slew rate
% note: don't support ramp sampling

global VCtl
global VObj

tStart=p.tStart; % start time
tEnd=p.tEnd; % end time
Area=abs(p.Area);   % trapezoid area m-1
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

Area  =Area/(VObj.Gyro/(2*pi));
GzAmp =Area /(tEnd - tStart);
tRamp=max(VCtl.MinUpdRate,GzAmp/VCtl.MaxSlewRate);   % ramp time

[GAmp,GTime]=StdTrap(tStart-tRamp, ...
                     tEnd+tRamp,   ...
                     tStart,               ...
                     tEnd,                 ...
                     GzAmp*sign(p.Area),2,2,2);

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http ://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

function mx = maxn(V)
% Max of V treated as a vector, equal to max(V(:))
mx = max(V(:));

author: mahmoudnafifi


text:
clc;
clearvars;
close all;
%rng default;

N = 1000;
A = spx.dict.simple.gaussian_mtx(N, N);
A = A' * A;
x = randn(N, 1);
b  = A * x;

[x, res, iter ] = cgsolve(A, b, 1e-2, N*5, 0);
spx.io.print.vector(x, 3);

options.tolerance = 1e-2;
options.max_iterations = N*5;
fprintf('\nTrying fast solver: \n\n')
result = spx.fast.cg(A, b, options);
x = result.x;
spx.io.print.vector(x, 3);

author: indigits

repo: sparse-plex


text:
filt_opt.J = 3;
filters1d = morlet_filter_bank_1d(32, filt_opt);
dual_filters1d = dual_filter_bank(filters1d);

rs = RandStream.create('mt19937ar','Seed',1234);
x = rs.randn(32,1);

wavelet_opt.oversampling = 2;
[x_phi,x_psi,meta_phi,meta_psi] = wavelet_1d(x(1:32),filters1d,wavelet_opt);
x1 = inverse_wavelet_1d(32, x_phi, x_psi, meta_phi, meta_psi, dual_filters1d);
assert(norm(x-x1)/norm(x)<0.02)

wavelet_opt.oversampling = 100;
[x_phi,x_psi,meta_phi,meta_psi] = wavelet_1d(x(1:32),filters1d,wavelet_opt);
x1 = inverse_wavelet_1d(32, x_phi, x_psi, meta_phi, meta_psi, dual_filters1d);
assert(norm(x-x1)/norm(x)<1e-14)

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% A demo script that demonstrates how to process a single video file using
% OpenFace and extract and visualize all of the features

clear

% The location executable will depend on the OS
if(isunix)
    executable = '"../../build/bin/FeatureExtraction"';
else
    executable = '"../../x64/Release/FeatureExtraction.exe"';
end

% Input file
in_dir = '../../samples/image_sequence';

% Where to store the output
output_dir = './processed_features/';

% This will take directory after -fdir and output all the features to directory
% after -out_dir
command = sprintf('%s -fdir "%s" -out_dir "%s" -verbose', executable, in_dir, output_dir);

if(isunix)
    unix(command);
else
    dos(command);
end

%% Demonstrating reading the output files

% Most of the features will be in the csv file in the output directory with
% the same name as the input file
[~,name,~] = fileparts(in_dir);
output_csv = sprintf('%s/%s.csv', output_dir, name);

% First read in the column names, to know which columns to read for
% particular features
tab = readtable(output_csv);
column_names = tab.Properties.VariableNames;

% Read all of the data
all_params  = dlmread(output_csv, ',', 1, 0);

% This indicates which frames were succesfully tracked

% Find which column contains success of tracking data and timestamp data
valid_ind = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'success'));
frame_ind = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'frame'));

% Extract tracking success data and only read those frame
valid_frames = logical(all_params(:,valid_ind));

% Get the timestamp data
frame_nums = all_params(valid_frames, frame_ind);

%% Finding which header line starts with p_ (basically model params)
shape_inds = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'p_'));

% Output rigid (first 6) and non-rigid shape parameters
shape_params  = all_params(valid_frames, shape_inds);

figure
plot(frame_nums, shape_params);
title('Shape parameters');
xlabel('Frame');

%% Demonstrate 2D landmarks
landmark_inds_x = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'x_'));
landmark_inds_y = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'y_'));

xs = all_params(valid_frames, landmark_inds_x);
ys = all_params(valid_frames, landmark_inds_y);

eye_landmark_inds_x = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'eye_lmk_x_'));
eye_landmark_inds_y = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'eye_lmk_y_'));

eye_xs = all_params(valid_frames, eye_landmark_inds_x);
eye_ys = all_params(valid_frames, eye_landmark_inds_y);

figure

for j = 1:size(xs,1)
    plot(xs(j,:), -ys(j,:), '.');
    hold on;
    plot(eye_xs(j,:), -eye_ys(j,:), '.r');
    hold off;

    xlim([min(xs(1,:)) * 0.5, max(xs(2,:))*1.4]);
    ylim([min(-ys(1,:)) * 1.4, max(-ys(2,:))*0.5]);
    xlabel('x (px)');
    ylabel('y (px)');
    drawnow
end

%% Demonstrate 3D landmarks
landmark_inds_x = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'X_'));
landmark_inds_y = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'Y_'));
landmark_inds_z = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'Z_'));

xs = all_params(valid_frames, landmark_inds_x);
ys = all_params(valid_frames, landmark_inds_y);
zs = all_params(valid_frames, landmark_inds_z);

eye_landmark_inds_x = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'eye_lmk_X_'));
eye_landmark_inds_y = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'eye_lmk_Y_'));
eye_landmark_inds_z = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'eye_lmk_Z_'));

eye_xs = all_params(valid_frames, eye_landmark_inds_x);
eye_ys = all_params(valid_frames, eye_landmark_inds_y);
eye_zs = all_params(valid_frames, eye_landmark_inds_z);

figure
for j = 1:size(xs,1)
    plot3(xs(j,:), ys(j,:), zs(j,:), '.');axis equal;
    hold on;
    plot3(eye_xs(j,:), eye_ys(j,:), eye_zs(j,:), '.r');
    hold off;
    xlabel('X (mm)');
    ylabel('Y (mm)');
    zlabel('Z (mm)');
    drawnow
end

%% Demonstrate AUs
au_reg_inds = cellfun(@(x) ~isempty(x) && x==5, strfind(column_names, '_r'));

aus = all_params(valid_frames, au_reg_inds);
figure
plot(frame_nums, aus);
title('Facial Action Units (intensity)');
xlabel('Time (s)');
ylabel('Intensity');
ylim([0,6]);

au_class_inds = cellfun(@(x) ~isempty(x) && x==5, strfind(column_names, '_c'));

aus = all_params(valid_frames, au_class_inds);
figure
plot(frame_nums, aus);
title('Facial Action Units (presense)');
xlabel('Time (s)');
ylim([0,2]);
%% Demo pose
pose_inds = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'pose_'));

pose = all_params(valid_frames, pose_inds);
figure
plot(frame_nums, pose);
title('Pose (rotation and translation)');
xlabel('Frame number');

%% Demo gaze
gaze_inds = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'gaze_angle'));

% Read gaze (x,y,z) for one eye and (x,y,z) for another
gaze = all_params(valid_frames, gaze_inds);

plot(frame_nums, gaze(:,1), 'DisplayName', 'Left - right');
hold on;
plot(frame_nums, gaze(:,2), 'DisplayName', 'Up - down');
xlabel('Frame number') % x-axis label
ylabel('Angle radians') % y-axis label
legend('show');
hold off;

%% Output HOG files
output_hog_file = sprintf('%s/%s.hog', output_dir, name);
[hog_data, valid_inds] = Read_HOG_file(output_hog_file);

%% Output aligned images
output_aligned_dir = sprintf('%s/%s_aligned/', output_dir, name);
img_files = dir([output_aligned_dir, '/*.bmp']);
imgs = cell(numel(img_files, 1));
for i=1:numel(img_files)
   imgs{i} = imread([ output_aligned_dir, '/', img_files(i).name]);
   imshow(imgs{i})
   drawnow
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function plane2 = normalizePlane(plane1)
%NORMALIZEPLANE Normalize parametric representation of a plane
%
%   PLANE2 = normalizePlane(PLANE1);
%   Transforms the plane PLANE1 in the following format:
%   [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where:
%   - (X0, Y0, Z0) is a point belonging to the plane
%   - (DX1, DY1, DZ1) is a first direction vector
%   - (DX2, DY2, DZ2) is a second direction vector
%   into another plane, with the same format, but with:
%   - (x0 y0 z0) is the closest point of plane to the origin
%   - (DX1 DY1 DZ1) has norm equal to 1
%   - (DX2 DY2 DZ2) has norm equal to 1 and is orthogonal to (DX1 DY1 DZ1)
%
%   See also:
%   planes3d, createPlane
%
%   ---------
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY
%   21/08/2009 compute origin after computation of vectors (more precise)
%       and add support for several planes.

% compute first direction vector
d1  = normalizeVector3d(plane1(:,4:6));

% compute second direction vector
n   = normalizeVector3d(planeNormal(plane1));
d2  = -normalizeVector3d(crossProduct3d(d1, n));

% compute origin point of the plane
origins = repmat([0 0 0], [size(plane1, 1) 1]);
p0 = projPointOnPlane(origins, [plane1(:,1:3) d1 d2]);

% create the resulting plane
plane2 = [p0 d1 d2];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [pdf,cdf] = VBA_binomial(x,n,p)
% PDF and CDF of binomial distribution
% function [pdf,cdf] = VBA_binomial(x,n,p)
% Note: CDF is P(X<x), i.e. it excludes x, whereas PDF is P(X=x)!
% IN:
%   - x: number of successes
%   - n: total number of draws
%   - p: probability of a success
% OUT:
%   - pdf: P(X=x), where X is a binomial random variable
%   - cdf: P(X<x), where X is a binomial random variable

fullpdf = zeros(max(x),1);
for i=0:max(x)
    fullpdf(i+1) = nchoosek(n,i).*(p.^i).*((1-p).^(n-i));
end
pdf = fullpdf(x+1);
fullcdf = cumsum([0;fullpdf]);
cdf = fullcdf(x+1);

author: MBB-team

repo: VBA-toolbox


text:
% brainstorm2eeglab - convert Brainstorm structures to EEGLAB dataset
%
% EEG = brainstorm2eeglab(bst, epochType);
%
% Inputs:
%    bst - Brainstorm string for a folder containing data epoch files or Brainstorm
%          data epoch structure
%    epochType - [string] epoch type. Only import these epochs which
%                Comment field contains this string
%
% Limitations:
%   - Only import data epochs
%   - Cannot import linked file
%
% Output:
%    EEG     - EEGLAB structure
%
% Author: Arnaud Delorme, UCSD

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

function EEG = brainstorm2eeglab(bst, epochType)

if nargin < 1
    help brainstorm2eeglab;
    return;
end
if nargin < 2
    epochType = '';
end

if ischar(bst)
    oriBst = bst;
    bst = dir(fullfile(bst, '*.mat'));
    if isempty(bst)
        error('Files not found')
    end
else
    oriBst = '';
    % Need brainstrorm
    try
        prot = bst_get('ProtocolInfo');
    catch
        brainstorm nogui
        try
            prot = bst_get('ProtocolInfo');
        catch
            error('Brainstorm not found. Brainstorm must be installed and running.')
        end
    end
    brainstorm_path = bst_get('BrainstormDbDir');
    protocol = prot.Comment;

    for iEpoch = 1:length(bst)
        bst(iEpoch).name = bst(iEpoch).FileName;
        bst(iEpoch).folder = fullfile(brainstorm_path, protocol, 'data');
    end
end

disp('Reading data epochs...');
allEpochs = {};
for iEpoch = 1:length(bst)
    if ~isequal(bst(iEpoch).name, 'brainstormstudy.mat')
        epochStruct = load('-mat', fullfile(bst(iEpoch).folder, bst(iEpoch).name));
        % epochStruct.Comment(1) contain the type of event
        % events not imported (in epochStruct.event)
        if ~isempty(epochType)
            if contains(epochStruct.Comment, epochType)
                allEpochs{end+1} = epochStruct.F;
            end
        else
            allEpochs{end+1} = epochStruct.F;
        end
    end
end
data = reshape([allEpochs{:}], size(allEpochs{1},1), size(allEpochs{1},2), length(allEpochs));

% create EEG structure
EEG = eeg_emptyset;
EEG.srate = mean(1./diff(epochStruct.Time));
EEG.data  = data; % microvolts
EEG.xmin  = epochStruct.Time(1);

% channel locations
if isempty(oriBst)
    chanFile = fullfile(bst(1).folder, bst(1).ChannelFile);
else
    chanFile = fullfile(oriBst, '..', '@default_study', 'channel.mat');
end
if exist(chanFile, 'file')
    chans = load('-mat', chanFile);
    for iChan = 1:length(chans.Channel)
        EEG.chanlocs(iChan).labels = chans.Channel(iChan).Name;
        if ~isempty(chans.Channel(iChan).Loc)
            EEG.chanlocs(iChan).X      = chans.Channel(iChan).Loc(1);
            EEG.chanlocs(iChan).Y      = chans.Channel(iChan).Loc(2);
            EEG.chanlocs(iChan).Z      = chans.Channel(iChan).Loc(3);
        end
        EEG.chanlocs(iChan).type   = chans.Channel(iChan).Type;
    end
    EEG.chanlocs = convertlocs(EEG.chanlocs, 'cart2all');
end
EEG = eeg_checkset(EEG);

% rescale data to microvolts
if ~isempty(EEG.chanlocs)
   chanTypes = strmatch('EEG', { EEG.chanlocs.type });
   if ~isempty(chanTypes)
       EEG.data(chanTypes,:,:) = EEG.data(chanTypes,:,:)*1e6; % assumes volt (this code should not be reached as channel types should always be defined)
   end
else
    EEG.data = EEG.data*1e6; % assumes volt (this code should not be reached as channel types should always be defined)
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [ n_data, n, c ] = bell_values ( n_data )

%*****************************************************************************80
%
%% BELL_VALUES returns some values of the Bell numbers.
%
%  Discussion:
%
%    The Bell number B(N) is the number of restricted growth functions on N.
%
%    Note that the Stirling numbers of the second kind, S^m_n, count the
%    number of partitions of N objects into M classes, and so it is
%    true that
%
%      B(N) = S^1_N + S^2_N + ... + S^N_N.
%
%    The Bell numbers were named for Eric Temple Bell.
%
%    In Mathematica, the function can be evaluated by
%
%      Sum[StirlingS2[n,m],{m,1,n}]
%
%    The Bell number B(N) is defined as the number of partitions (of
%    any size) of a set of N distinguishable objects.
%
%    A partition of a set is a division of the objects of the set into
%    subsets.
%
%  Example:
%
%    There are 15 partitions of a set of 4 objects:
%
%      (1234),
%      (123) (4),
%      (124) (3),
%      (12) (34),
%      (12) (3) (4),
%      (134) (2),
%      (13) (24),
%      (13) (2) (4),
%      (14) (23),
%      (1) (234),
%      (1) (23) (4),
%      (14) (2) (3),
%      (1) (24) (3),
%      (1) (2) (34),
%      (1) (2) (3) (4).
%
%    and so B(4) = 15.
%
%  First values:
%
%     N         B(N)
%     0           1
%     1           1
%     2           2
%     3           5
%     4          15
%     5          52
%     6         203
%     7         877
%     8        4140
%     9       21147
%    10      115975
%
%  Recursion:
%
%    B(I) = sum ( 1 <= J <=I ) Binomial ( I-1, J-1 ) * B(I-J)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the Bell number.
%
%    Output, integer C, the value of the Bell number.
%
  n_max = 11;

  c_vec = [ ...
    1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 ];

  n_vec = [ ...
     0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data);
    c = c_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m