text:
function [success,msg,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  pxpermm] = ...
  ReadArenaParameters_Ctrax(varargin)

success = false;
msg = ''; %#ok<NASGU>

[inmoviefile,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  pxpermm,intrxfile,annfile,inperframedir] = myparse(varargin,...
  'inmoviefile','',...
  'arenatype','None','arenacenterx',0,'arenacentery',0,...
  'arenaradius',123,'arenawidth',123,'arenaheight',123,...
  'pxpermm',1,'intrxfile','','annfile','','inperframedir','');

% try to read from the trx file first
if isempty(intrxfile),
  msg = 'Input trxfile not yet set';
  return;
end

try
  [trx,~,success1] = load_tracks(intrxfile,inmoviefile,...
    'dosave',false,'annname',annfile,'verbose',false);
  if ~success1,
    msg = sprintf('Could not load tracks from trxfile %s',intrxfile);
    return;
  end
catch ME,
  msg = sprintf('Could not load from trxfile %s: %s',intrxfile,getReport(ME));
  return;
end

readarenafrom = '';
readpxpermmfrom = '';
readoff = false;
if isfield(trx,'pxpermm'),
  newpxpermm = nanmean([trx.pxpermm]);
  if ~isnan(newpxpermm) && newpxpermm > 0,
    pxpermm = newpxpermm;
    readpxpermmfrom = 'intrxfile';
  end
  if all(isfield(trx,{'x','x_mm','y','y_mm'})),
    offx = nanmean(([trx.x] - [trx.x_mm])*pxpermm);
    offy = nanmean(([trx.y] - [trx.y_mm])*pxpermm);
    readoff = true;
  end
end

if isfield(trx,'arena'),

  % in px
  if all(isfield(trx(1).arena,{'x','y','r'})),
    arenacenterx = trx(1).arena.x;
    arenacentery = trx(1).arena.y;
    arenaradius = trx(1).arena.r;
    arenatype = 'Circle';
    readarenafrom = 'intrxfile';
  elseif ~isempty(readpxpermmfrom) && readoff && ...
      all(isfield(trx(1).arena,{'arena_radius_mm','arena_center_mm_x','arena_center_mm_y'})),
    arenacenterx = trx(1).arena_center_mm_x * pxpermm + offx;
    arenacentery = trx(1).arena_center_mm_y * pxpermm + offy;
    arenaradius = trx(1).arena_radius_mm * pxpermm;
    readarenafrom = 'intrxfile';
  elseif all(isfield(arena,{'tl_px','tr_px','bl_px','br_px'})),
    tl_px = trx(1).arena.tl;
    tr_px = trx(1).arena.tr;
    bl_px = trx(1).arena.bl;
    br_px = trx(1).arena.br;
    arenawidth = mean([abs(tr_px(1)-tl_px(1)),abs(br_px(1)-bl_px(1))]);
    arenaheight = mean([abs(tr_px(2)-br_px(2)),abs(tl_px(2)-bl_px(2))]);
    arenacenterx = mean([tl_px(1),tr_px(1),bl_px(1),br_px(1)]);
    arenacentery = mean([tl_px(2),tr_px(2),bl_px(2),br_px(2)]);
    arenatype = 'Rectangle';
    readarenafrom = 'intrxfile';
  elseif ~isempty(readpxpermmfrom) && readoff && all(isfield(arena,{'tl','tr','bl','br'})),
    tl_px = trx(1).arena.tl(:)'*pxpermm + [offx,offy];
    tr_px = trx(1).arena.tr(:)'*pxpermm + [offx,offy];
    bl_px = trx(1).arena.bl(:)'*pxpermm + [offx,offy];
    br_px = trx(1).arena.br(:)'*pxpermm + [offx,offy];
    arenawidth = mean([abs(tr_px(1)-tl_px(1)),abs(br_px(1)-bl_px(1))]);
    arenaheight = mean([abs(tr_px(2)-br_px(2)),abs(tl_px(2)-bl_px(2))]);
    arenacenterx = mean([tl_px(1),tr_px(1),bl_px(1),br_px(1)]);
    arenacentery = mean([tl_px(2),tr_px(2),bl_px(2),br_px(2)]);
    arenatype = 'Rectangle';
    readarenafrom = 'intrxfile';
  end

end

if isempty(readarenafrom) && ~isempty(annfile),
  % try to read from annfile
  [newarenacenterx,newarenacentery,newarenaradius,dosetcirculararena] = ...
    read_ann('arena_center_x','arena_center_y','arena_radius','do_set_circular_arena');
  if dosetcirculararena && newarenaradius > 0 && ~any(isnan([newarenaradius,newarenacenterx,newarenacentery])),
    arenacenterx = newarenacenterx;
    arenacentery = newarenacentery;
    arenaradius = newarenaradius;
    arenatype = 'Circle';
    readarenafrom = 'annfile';
  end
end

success = ~isempty(readarenafrom);
if ~isempty(readarenafrom);
  msg = sprintf('Read arena from %s, ',readarenafrom);
else
  msg = sprintf('Did not read arena, ');
end
if ~isempty(readpxpermmfrom);
  msg = [msg, sprintf('read pxpermm from %s.',readpxpermmfrom)];
else
  msg = [msg, sprintf('Did not read pxpermm.')];
end


text:
function r = taylor(a,str)
%TAYLOR       Taylor class constructor
%
%  r = taylor(a)
%
%An explicit call of the constructor is only necessary to initialize
%  a constant to be of type taylor. Otherwise, any operation
%  with a dependent variable produces a result of type taylor.
%
%For more details try
%
%  help taylorinit
%
%and demotaylor.
%

%taylor.size is size of input
%taylor.t stored as column vector of length INTLAB_TAYLOR_ORDER; in case of
%vector or matrix input, columns of taylor.t are the Taylor coefficients
%

% written  05/21/09     S.M. Rump
% modified 02/28/10     S.M. Rump  rounding
% modified 08/26/12     S.M. Rump  global variables removed
%

  superiorto('intval');

  if nargin==0
    r.size = [];
    r.t = [];
    r = class(r,'taylor');
    return
  end

  INTLAB_TAYLOR_ORDER = getappdata(0,'INTLAB_TAYLOR_ORDER');

  if INTLAB_TAYLOR_ORDER==0
    error('no dependent variables initialized for use of Taylor')
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin==1

    if isa(a,'taylor')
      r = a;
    else
      r.size = size(a);
      len = prod(r.size);
      r.t = [ a(:).' ; zeros(INTLAB_TAYLOR_ORDER,len) ];
      r = class(r,'taylor');
    end

  elseif nargin==2

    if ischar(str)

      if isequal(str,'taylorinit')         % call by taylorinit

        r.size = size(a.init);
        len = prod(r.size);
        r.t = [ a.init(:).' ; ones(1,len) ; zeros(INTLAB_TAYLOR_ORDER-1,len) ];
        r = class(r,'taylor');

      elseif isequal(str,'taylor')         % call by @intval\taylor

        r.size = size(a.init);
        len = prod(r.size);
        r.t = intval([ a.init(:).' ; zeros(INTLAB_TAYLOR_ORDER,len) ]);
        r = class(r,'taylor');

      elseif isequal(str,'taylorintval')   % call by @intval\intval

        r = a;
        r.t = intval(r.t);

      elseif isequal(str,'random')          % generates .t randomly, only for test purposes

        if isa(a,'struct')                  % input interval
          a = a.init;
        end
        r.size = size(a);
        len = prod(r.size);
        r.t = [ a(:).' ; randn(INTLAB_TAYLOR_ORDER,len) ];
        r = class(r,'taylor');

      else
        error('invalid call of taylor constructor')
      end

    end

  else

    error('invalid call of constructor taylor')

  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function	[out,dout]=rodrigues(in)

% RODRIGUES	Transform rotation matrix into rotation vector and viceversa.
%
%		Sintax:  [OUT]=RODRIGUES(IN)
% 		If IN is a 3x3 rotation matrix then OUT is the
%		corresponding 3x1 rotation vector
% 		if IN is a rotation 3-vector then OUT is the
%		corresponding 3x3 rotation matrix
%

%%
%%		Copyright (c) March 1993 -- Pietro Perona
%%		California Institute of Technology
%%

%% ALL CHECKED BY JEAN-YVES BOUGUET, October 1995.
%% FOR ALL JACOBIAN MATRICES !!! LOOK AT THE TEST AT THE END !!

%% BUG when norm(om)=pi fixed -- April 6th, 1997;
%% Jean-Yves Bouguet

%% Add projection of the 3x3 matrix onto the set of special ortogonal matrices SO(3) by SVD -- February 7th, 2003;
%% Jean-Yves Bouguet

% BUG FOR THE CASE norm(om)=pi fixed by Mike Burl on Feb 27, 2007

[m,n] = size(in);
%bigeps = 10e+4*eps;
bigeps = 10e+20*eps;

if ((m==1) & (n==3)) | ((m==3) & (n==1)) %% it is a rotation vector
    theta = norm(in);
    if theta < eps
        R = eye(3);

        %if nargout > 1,

        dRdin = [0 0 0;
            0 0 1;
            0 -1 0;
            0 0 -1;
            0 0 0;
            1 0 0;
            0 1 0;
            -1 0 0;
            0 0 0];

        %end;

    else
        if n==length(in)  in=in'; end; 	%% make it a column vec. if necess.

        %m3 = [in,theta]

        dm3din = [eye(3);in'/theta];

        omega = in/theta;

        %m2 = [omega;theta]

        dm2dm3 = [eye(3)/theta -in/theta^2; zeros(1,3) 1];

        alpha = cos(theta);
        beta = sin(theta);
        gamma = 1-cos(theta);
        omegav=[[0 -omega(3) omega(2)];[omega(3) 0 -omega(1)];[-omega(2) omega(1) 0 ]];
        A = omega*omega';

        %m1 = [alpha;beta;gamma;omegav;A];

        dm1dm2 = zeros(21,4);
        dm1dm2(1,4) = -sin(theta);
        dm1dm2(2,4) = cos(theta);
        dm1dm2(3,4) = sin(theta);
        dm1dm2(4:12,1:3) = [0 0 0 0 0 1 0 -1 0;
            0 0 -1 0 0 0 1 0 0;
            0 1 0 -1 0 0 0 0 0]';

        w1 = omega(1);
        w2 = omega(2);
        w3 = omega(3);

        dm1dm2(13:21,1) = [2*w1;w2;w3;w2;0;0;w3;0;0];
        dm1dm2(13: 21,2) = [0;w1;0;w1;2*w2;w3;0;w3;0];
        dm1dm2(13:21,3) = [0;0;w1;0;0;w2;w1;w2;2*w3];

        R = eye(3)*alpha + omegav*beta + A*gamma;

        dRdm1 = zeros(9,21);

        dRdm1([1 5 9],1) = ones(3,1);
        dRdm1(:,2) = omegav(:);
        dRdm1(:,4:12) = beta*eye(9);
        dRdm1(:,3) = A(:);
        dRdm1(:,13:21) = gamma*eye(9);

        dRdin = dRdm1 * dm1dm2 * dm2dm3 * dm3din;

    end;
    out = R;
    dout = dRdin;

    %% it is prob. a rot matr.
elseif ((m==n) & (m==3) & (norm(in' * in - eye(3)) < bigeps)...
        & (abs(det(in)-1) < bigeps))
    R = in;

    % project the rotation matrix to SO(3);
    [U,S,V] = svd(R);
    R = U*V';

    tr = (trace(R)-1)/2;
    dtrdR = [1 0 0 0 1 0 0 0 1]/2;
    theta = real(acos(tr));

    if sin(theta) >= 1e-4,

        dthetadtr = -1/sqrt(1-tr^2);

        dthetadR = dthetadtr * dtrdR;
        % var1 = [vth;theta];
        vth = 1/(2*sin(theta));
        dvthdtheta = -vth*cos(theta)/sin(theta);
        dvar1dtheta = [dvthdtheta;1];

        dvar1dR =  dvar1dtheta * dthetadR;

        om1 = [R(3,2)-R(2,3), R(1,3)-R(3,1), R(2,1)-R(1,2)]';

        dom1dR = [0 0 0 0 0 1 0 -1 0;
            0 0 -1 0 0 0 1 0 0;
            0 1 0 -1 0 0 0 0 0];

        % var = [om1;vth;theta];
        dvardR = [dom1dR;dvar1dR];

        % var2 = [om;theta];
        om = vth*om1;
        domdvar = [vth*eye(3) om1 zeros(3,1)];
        dthetadvar = [0 0 0 0 1];
        dvar2dvar = [domdvar;dthetadvar];

        out = om*theta;
        domegadvar2 = [theta*eye(3) om];

        dout = domegadvar2 * dvar2dvar * dvardR;

    else
        if tr > 0; 			% case norm(om)=0;

            out = [0 0 0]';

            dout = [0 0 0 0 0 1/2 0 -1/2 0;
                0 0 -1/2 0 0 0 1/2 0 0;
                0 1/2 0 -1/2 0 0 0 0 0];
        else

            % case norm(om)=pi;
            if(0)

                %% fixed April 6th by Bouguet -- not working in all cases!
                out = theta * (sqrt((diag(R)+1)/2).*[1;2*(R(1,2:3)>=0)'-1]);
                %keyboard;

            else

                % Solution by Mike Burl on Feb 27, 2007
                % This is a better way to determine the signs of the
                % entries of the rotation vector using a hash table on all
                % the combinations of signs of a pairs of products (in the
                % rotation matrix)

                % Define hashvec and Smat
                hashvec = [0; -1; -3; -9; 9; 3; 1; 13; 5; -7; -11];
                Smat = [1,1,1; 1,0,-1; 0,1,-1; 1,-1,0; 1,1,0; 0,1,1; 1,0,1; 1,1,1; 1,1,-1;
                    1,-1,-1; 1,-1,1];

                M = (R+eye(3,3))/2;
                uabs = sqrt(M(1,1));
                vabs = sqrt(M(2,2));
                wabs = sqrt(M(3,3));

                mvec = [M(1,2), M(2,3), M(1,3)];
                syn  = ((mvec > 1e-4) - (mvec < -1e-4)); % robust sign() function
                hash = syn * [9; 3; 1];
                idx = find(hash == hashvec);
                svec = Smat(idx,:)';

                if ( isempty(svec) )
                    out = theta * (sqrt((diag(R)+1)/2).*[1;2*(R(1,2:3)>=0)'-1]);
                else
                    out = theta * [uabs; vabs; wabs] .* svec;
                end

            end;

            if nargout > 1,
                fprintf(1,'WARNING!!!! Jacobian domdR undefined!!!\n');
                dout = NaN*ones(3,9);
            end;
        end;
    end;

else
    error('Neither a rotation matrix nor a rotation vector were provided');
end;

return;

%% test of the Jacobians:

%%%% TEST OF dRdom:
om = randn(3,1);
dom = randn(3,1)/1000000;

[R1,dR1] = rodrigues(om);
R2 = rodrigues(om+dom);

R2a = R1 + reshape(dR1 * dom,3,3);

gain = norm(R2 - R1)/norm(R2 - R2a)

%%% TEST OF dOmdR:
om = randn(3,1);
R = rodrigues(om);
dom = randn(3,1)/10000;
dR = rodrigues(om+dom) - R;

[omc,domdR] = rodrigues(R);
[om2] = rodrigues(R+dR);

om_app = omc + domdR*dR(:);

gain = norm(om2 - omc)/norm(om2 - om_app)

%%% OTHER BUG: (FIXED NOW!!!)

omu = randn(3,1);
omu = omu/norm(omu)
om = pi*omu;
[R,dR]= rodrigues(om);
[om2] = rodrigues(R);
[om om2]

%%% NORMAL OPERATION

om = randn(3,1);
[R,dR]= rodrigues(om);
[om2] = rodrigues(R);
[om om2]

return

% Test: norm(om) = pi

u = randn(3,1);
u = u / sqrt(sum(u.^2));
om = pi*u;
R = rodrigues(om);

R2 = rodrigues(rodrigues(R));

norm(R - R2)

author: zouchuhang

repo: LayoutNet


text:
function u_est = CG_MultiScale_LAP3D(I1, I2, FilterSizes, PreFilt, MedFilt, uin)
% The function implements a multi-scale framework for the 3D version of the
% LAP optical flow algorithm. Instead of downsampling the images, the
% framework changes the size of the all-pass filters used in the LAP
% algorithm. The filter basis used in the LAP algorithm spans the
% derivatives of a Gaussian filter.
% Note that this implementation is for greyscale images only.
%
% Input parameters:
%       I1 and I2   -> Input images of size M by N (Greyscale)
%       FilterSizes -> Controls the size of the filters (vector with
%                      values for R
%                      i.e. filter size is [2R + 1] by [2R + 1]
%       PreFilt     -> Optional parameter to decide whether to highpass
%                       filter the images. (Default = 1 (Yes), 0 = (No))
%       MedFilt     -> Optional parameter to decide whether to median
%                       filter the flow at fine scales. (Default = 1 (Yes), 0 = (No))
%
% Outputs:
%       u_est       -> Estimate of the optical flow
%
%
% Date: 08/07/2015, Author: Chris Gilliam
%
% modifications for MRI application:
% - different input parameter: Filter size instead of upper parameter R
% - initialization in single precision (except filter kernel)
% 26.08.2015, Verena Neumann/Thomas Kuestner

if nargin <= 3,
    PreFilt = 1;
    MedFilt = 1;
elseif nargin <= 4,
    MedFilt = 1;
end

% Obtain the dimensions of the images:
[M, N, P] = size(I1);

% Initialise filter functions:
funs = Filter_Functions;

% Local Gaussian filter (used in high pass filtering)
h = exp(-(-2:2).^2/2);
h = h./sum(h(:));

% Initial optical flow estimate
if(nargin < 6)
    u_holder = repmat({zeros(M,N,P,'single')}, 1, 3);
else
    u_holder = uin;
end
% u_holder = repmat({zeros(M,N,P)}, 1, 3);

% define half support of filters (i.e. R)
amp_array = FilterSizes;    % 2.^(Level_Max:-1:Level_Min);

% Initialise local counter:
num_level = 0;

% Local I1 variable:
if PreFilt == 1,
    im1 = I1 - funs.Filter_General(I1, h); % High pass filtering using gaussian filter h
%       im1 = funs.Filter_Laplacian(im1);  % high-pass using generalised Laplacian filter
else
    im1 = I1;
end

% Start estimating the optical flow
for l = 1:length(amp_array),
    num_level = num_level + 1;
    disp(['Level ', int2str(num_level), '/', int2str(size(amp_array,2))]);

    % Define filter parameter R at each iteration
    amp_size = amp_array(l);

    % Load Filter Basis:
    Basis_Set = loadbasis(3,amp_size);

    if l == 1,
        % No warping for first iteration:
        if PreFilt == 1,
            I2_shift = I2 - funs.Filter_General(I2, h); % High pass filtering using gaussian filter h
%             I2_shift = funs.Filter_Laplacian(I2);       % high-pass using generalised Laplacian filter
        else
            I2_shift = I2;
        end
    else
        % Warp I1 closer to I2 using current optical flow estimate
        if amp_size > 2,
            I2_shift = imshift_3D(I2,{-u_holder{1}, -u_holder{2}, -u_holder{3}}, 'shiftedlinear');
        else
            I2_shift = imshift_3D(I2,{-u_holder{1}, -u_holder{2}, -u_holder{3}});
        end
        if PreFilt == 1,
            I2_shift = I2_shift - funs.Filter_General(I2_shift, h); % High pass filtering using gaussian filter h
%           I2_shift = funs.Filter_Laplacian(I2_shift); % high-pass using generalised Laplacian filter
        end
    end

    % Using basis functions estimate optical flow on a local scale:
    [uest_Orig, ~] = optiflowFilter3D(im1, I2_shift, Basis_Set);

    % Clean optical flow
    % Stage 1: Remove nan's in the optical flow using inpainting:
    if sum(isnan(uest_Orig{1}(:))) >= M*N*P,
        error('All NaN. Suggest reduce Level_Num');
    end
    [uest_Clean, ~] = cleanOF3D(uest_Orig);

    % Stage 2: Remove flow elements that corresponding to large warping
    % errors.
    R = round(2*amp_size);
    k1 = -R:R;
    uest_Clean = cellfun(@funs.Filter_Gauss, uest_Clean, {R,R,R}, {k1,k1,k1}, 'uni',false);

    % Rescale optical flow and add to estimate
    u_holder = cellfun(@plus,u_holder,uest_Clean,'uniformoutput',false);

    % Refinement of OF at highest level based on errors in shifted
    %           images
    if amp_size <= 2 && MedFilt == 1,
        u_holder = Cleaning_Procedure(u_holder);
    end

end

u_est = u_holder;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Embedded functions %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u_out = Cleaning_Procedure(u_in)
% function cleans the estimate of the optical flow. The cleaning process
% comprises a robust smoothing stage using two Median filters (of different
% sizes)
%

% Define size of median filters
B1 = 11;
B2 = 3;

% Two part median filtering:
% fine scale
u_out = cellfun(@medfilt3, u_in, {B2,B2,B2}, {'symmetric','symmetric','symmetric'},'uni',false);
% coarse scale
u_out = cellfun(@medfilt3, u_out, {B1,B1,B1}, {'symmetric','symmetric','symmetric'},'uni',false);
end

author: thomaskuestner


text:
% errordlg2() - Makes a popup dialog box with the specified message and (optional)
%               title.
%
% Usage:
%   errordlg2(Prompt, Title);
%
% Example:
%   errordlg2('Explanation of error','title of error');
%
% Input:
%   Prompt  -   A text string explaning why the user is seeing this error message.
%   Title   _   A text string that appears in the title bar of the error message.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 12 August 2002
%
% See also: inputdlg2(), questdlg2()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function errordlg2(Prompt, Title);

if exist('beep') == 5
	beep;
else
	disp(char(7));
end;
if nargin <2
	Title = 'Error';
end;
if ~ismatlab, error(Prompt); end;
questdlg2(Prompt, Title, 'OK', 'OK');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [U,Ud,data] = linear_elasticity(V,F,b,bc,varargin)
  % LINEAR_ELASTICITY Compute the deformation of a 2D solid object according to
  % a linear model of elasticity, assuming a linear isotropic material.
  %
  % U = linear_elasticity(V,F,b,bc)
  % [U,Ud,K] = linear_elasticity(V,F,b,bc,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   V  #V by d list of vertex positions
  %   F  #F by d+1 list of element indices into V
  %   b  #b list of indices into V of fixed vertices
  %   bc #bc by d list of fixed vertex positions
  %   Optional:
  %     'Lambda'  followed by first Lamé parameter {1.7423333}, scalar
  %       (homogeneous) or #F by 1 list of per-element values
  %     'Mu'  followed by shear modulus {0.0115}, scalar (homogeneous) or #F by
  %       1 list of per-element values
  %     'Young'  followed by Young's modulus, scalar (homogeneous) or #F by 1
  %       list of per-element values
  %     'Nu'  followed by Poisson's ratio, scalar (homogeneous) or #F by 1 list
  %       of per-element values
  %     'U0'  followed by #V by d list of previous displacements
  %     'Ud0'  followed by #V by d list of previous velocities: (U0 - Um1)/dt
  %     'BodyForces'  followed by #V by d list of body forces
  %     'TimeStep' followed by time step {0.1}
  %     'Data'  see output {[]}
  % Outputs:
  %   U  #V by d list of vertex displacements
  %   Ud  #V by d list of vertex velocities
  %   data  precomputation data
  %     data.A  #F*(d*(d+1)/2) by #F*(d*(d+1)/2) diagonal element area matrix
  %     data.K  #V*d by #V*d sparse stiffness matrix
  %     data.M  #V*d by #V*d sparse mass matrix
  %     data.strain  #F*(d*(d+1)/2) by #V*d sparse strain matrix
  %     data.dt  timestep
  %     data.C  #F**(d*(d+1)/2) by #F**(d*(d+1)/2) sparse constituitive model matrix
  %     data.mqwf  precomputation for implicit solve (from min_quad_with_fixed)
  %     data.solve  function handle for conducting implicit step
  %
  % Example:
  %   % Fit to half the unit square
  %   V = V/(2*max(max(V)-min(V)));
  %   % Initialize as stretched object
  %   U = 1.5*V-V;
  %   Ud = zeros(size(V));
  %   t = tsurf(F,V+U);
  %   axis equal;
  %   axis manual;
  %   while true
  %     [U,Ud] = linear_elasticity(V,F,[],[],'U0',U,'Ud0',Ud);
  %     t.Vertices = V+U;
  %     drawnow;
  %   end
  %

  data = [];
  % Time step
  dt = 0.1;
  % Silicone rubber: http://www.azom.com/properties.aspx?ArticleID=920
  mu = 0.0115;
  % Bulk modulus
  K = 1.75;
  lambda = K-2/3*mu;
  young = [];
  nu = [];
  U0 = zeros(size(V));
  Ud0 = zeros(size(V));
  fext = zeros(size(V));
  %% Parameters so that off-diagonals _should_ be zero
  %lambda = 1;
  %mu = -2*lambda;

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Lambda','Mu','Nu','Young','U0','Ud0','BodyForces','TimeStep','Data'}, ...
    {'lambda','mu','nu','young','U0','Ud0',      'fext',      'dt','data'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  first_solve = false;
  if isempty(data)
    first_solve = true;
    tic;
    data.dt = dt;
    assert( ...
      (~isempty(lambda) && ~isempty(mu))||(~isempty(young) && ~isempty(nu)), ...
      'Must define either lambda and mu or young and nu');
    if (~isempty(young) && ~isempty(nu))
      lambda = young.*nu./((1+nu).*(1-2.*nu));
      mu = .5.*young./(1+nu);
    end
    % young = mu*(3*lambda+2*mu)/(lambda+mu);
    % nu = lambda/(2*(lambda+mu));

    [data.K,data.C,data.strain,data.A,data.M] = ...
      linear_elasticity_stiffness(V,F,'Lambda',lambda,'Mu',mu);

    % ∇⋅σ + F = ü
    % Ku₂ + MF = M(u₂-2u₁+u₀)/dt²
    % dt²Ku₂ + dt²MF = M(u₂-2u₁+u₀)
    % (dt²K - M)u₂  = -dt²MF + M(-2u₁+u₀)
    % -(dt²K - M)u₂  = dt²MF - M(-2u₁+u₀)
    % (M-dt²K)u₂  = dt²MF + M(2u₁-u₀)
    % (M-dt²K)u₂  = M*(dt²F + 2u₁-u₀)
    % ud₀ = (u₁-u₀)/dt
    % dt*ud₀ = u₁-u₀
    % (M-dt²K)u₂  = M*(dt²F + u₁ + u₁-u₀)
    % (M-dt²K)u₂  = M*(dt²F + u₁ + dt*ud₀)
    A = data.M+data.dt^2*data.K;
    % ud = (u - u0)/dt
    % udd = ((u - u0)-(u0-um1)/dt²
    % udd = (u - 2u0 +um1)/dt²
    % ud*dt = (u - u0)
    % ud*dt - u = -u0
    % u - ud*dt = u0
    % udd*dt² = u - 2u0 +um1
    % udd*dt² - u + 2u0 = um1
    % 2u0-um1
    % 2(u - ud*dt)-(udd*dt² - u + 2u0)
    % 2u - 2ud*dt-udd*dt² + u - 2u0
    % 3u - 2ud*dt-udd*dt² - 2(u - ud*dt)
    % 3u - 2ud*dt-udd*dt² - 2u + ud*dt
    % u-ud*dt-udd*dt²

  end
  B = data.M*(data.dt^2*fext(:) + U0(:) + data.dt*Ud0(:));
  if first_solve
    % Fix each coordinate
    bb = reshape(bsxfun(@plus,reshape(b,[],1),(0:size(V,2)-1)*size(V,1)),1,[]);
    % solve once to set data.mqwf
    [U,data.mqwf] = min_quad_with_fixed(A,-2*B,bb,bc(:));
    data.solve  = @(B,bc) min_quad_with_fixed(A,-2*B,bb,bc(:),[],[],data.mqwf);
  else
    U = data.solve(B,bc(:));
  end
  U = reshape(U,size(V));
  Ud = (U-U0)/data.dt;

end

author: alecjacobson

repo: gptoolbox


text:
function v = nan2zero(v)

ind = isnan(v.x);
v.x(ind) = 0;
v.y(ind) = 0;
v.z(ind) = 0;

end

author: mtex-toolbox

repo: mtex

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


text:
function title = p06_title ( )

%*****************************************************************************80
%
%% P06_TITLE returns a title for problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = '#6: The superellipse with superelliptical hole.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hout=suptitle(str, varargin)
%SUPTITLE Puts a title above all subplots.
%	SUPTITLE('text') adds text to the top of the figure
%	above all subplots (a "super title"). Use this function
%	after all subplot commands.

% Drea Thomas 6/15/95 drea@mathworks.com
% John Cristion 12/13/00 modified
% Mark Histed 03/13/04 histed@mit.edu: fix disappearing legend on last plot
%
% $Id: suptitle.m,v 1.2 2004/03/13 22:17:47 histed Exp $

% Warning: If the figure or axis units are non-default, this
% will break.

if exist('varargin', 'var')
    for ii = 1:2:length(varargin)
       switch lower(varargin{ii})
           case {'fs', 'fontsize'}
               fs = varargin{ii+1};
       end
    end
end

% Parameters used to position the supertitle.

% Amount of the figure window devoted to subplots
plotregion = .92;

% Y position of title in normalized coordinates
titleypos  = .95;

% Fontsize for supertitle
if ~exist('fs', 'var') fs = get(gcf,'defaultaxesfontsize'); end

% Fudge factor to adjust y spacing between subplots
fudge=1;

haold = gca;
figunits = get(gcf,'units');

% Get the (approximate) difference between full height (plot + title
% + xlabel) and bounding rectangle.

	if (~strcmp(figunits,'pixels')),
		set(gcf,'units','pixels');
		pos = get(gcf,'position');
		set(gcf,'units',figunits);
	else,
		pos = get(gcf,'position');
	end
	ff = (fs-4)*1.27*5/pos(4)*fudge;

        % The 5 here reflects about 3 characters of height below
        % an axis and 2 above. 1.27 is pixels per point.

% Determine the bounding rectange for all the plots

% h = findobj('Type','axes');

% findobj is a 4.2 thing.. if you don't have 4.2 comment out
% the next line and uncomment the following block.

h = findobj(gcf,'Type','axes');  % Change suggested by Stacy J. Hills

% If you don't have 4.2, use this code instead
%ch = get(gcf,'children');
%h=[];
%for i=1:length(ch),
%  if strcmp(get(ch(i),'type'),'axes'),
%    h=[h,ch(i)];
%  end
%end

max_y=0;
min_y=1;

oldtitle =0;
for i=1:length(h),
	if (~strcmp(get(h(i),'Tag'),'suptitle')),
		pos=get(h(i),'pos');
		if (pos(2) < min_y), min_y=pos(2)-ff/5*3;end;
		if (pos(4)+pos(2) > max_y), max_y=pos(4)+pos(2)+ff/5*2;end;
	else,
		oldtitle = h(i);
	end
end

if max_y > plotregion,
	scale = (plotregion-min_y)/(max_y-min_y);
	for i=1:length(h),
		pos = get(h(i),'position');
		pos(2) = (pos(2)-min_y)*scale+min_y;
		pos(4) = pos(4)*scale-(1-scale)*ff/5*3;
		set(h(i),'position',pos);
	end
end

np = get(gcf,'nextplot');
set(gcf,'nextplot','add');
if (oldtitle),
	delete(oldtitle);
end
ha=axes('pos',[0 1 1 1],'visible','off','Tag','suptitle');
ht=text(.5,titleypos-1,str);set(ht,'horizontalalignment','center','fontsize',fs);
set(gcf,'nextplot',np);
axes(haold);

% fix legend if one exists
legH = legend;
if ~isempty(legH)
    axes(legH);
end

if nargout,
	hout=ht;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function net = resnet52_new_hope()

%----------------------img cnn----------------------
netStruct = load('./data/imagenet-resnet-50-dag.mat') ;
net = dagnn.DagNN.loadobj(netStruct) ;
net.removeLayer('fc1000');
net.removeLayer('prob');

for i = 1:numel(net.params)
    if(mod(i,2)==0)
        net.params(i).learningRate= 0; %0.02;
    else net.params(i).learningRate= 0; %0.001;
    end
    net.params(i).weightDecay=0; %1;
end

net.params(1).learningRate = 0; %1e-5;

fc1Block = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc1',fc1Block,{'pool5'},{'fc1'},{'fc1f'});
net.addLayer('fc1bn',dagnn.BatchNorm(),{'fc1'},{'fc1bn'},...
    {'fc1bn_w','fc1bn_b','fc1bn_m'});
net.addLayer('fc1x',dagnn.ReLU(),{'fc1bn'},{'fc1bnx'});
fc1Block = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc1_1',fc1Block,{'fc1bnx'},{'fc1_1'},{'fc1_1f'});
net.addLayer('fc1_1bn',dagnn.BatchNorm(),{'fc1_1'},{'fc1_1bn'},...
    {'fc1_1bn_w','fc1_1bn_b','fc1_1bn_m'});
net.addLayer('fc1_1x',dagnn.ReLU(),{'fc1_1bn'},{'fc1_1bnx'});
net.addLayer('dropout',dagnn.DropOut('rate',0.5),{'fc1_1bnx'},{'fc1_1bnxd'});

%----------------------char cnn----------------------
% input is 1*32*20074*16
fc2Block = dagnn.Conv('size',[1 1 29972 300],'hasBias',true,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2',fc2Block,{'data2'},{'fc2'},{'fc2f','fc2b'});
net.addLayer('fc2bn',dagnn.BatchNorm(),{'fc2'},{'fc2bn'},...
    {'fc2bn_w','fc2bn_b','fc2bn_m'});
%net.addLayer('fc2x',dagnn.ReLU(),{'fc2bn'},{'fc2bnx'});
%net.addLayer('dropout_diction',dagnn.DropOut('rate',0.5),{'fc2bn'},{'fc2bnd'});
% 32*256
convBlock = dagnn.Conv('size',[1 1 300 128],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1_1',convBlock,{'fc2bn'},{'fc2_1_1'},{'fc2_1_1f'});
net.addLayer('fc2_1_1bn',dagnn.BatchNorm(),{'fc2_1_1'},{'fc2_1_1bn'},...
    {'fc2_1_1bn_w','fc2_1_1bn_b','fc2_1_1bn_m'});
net.addLayer('fc2_1_1x',dagnn.ReLU(),{'fc2_1_1bn'},{'fc2_1_1bnx'});

convBlock = dagnn.Conv('size',[1 2 128 128],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
net.addLayer('fc2_1_2',convBlock,{'fc2_1_1bnx'},{'fc2_1_2'},{'fc2_1_2f'});
net.addLayer('fc2_1_2bn',dagnn.BatchNorm(),{'fc2_1_2'},{'fc2_1_2bn'},...
    {'fc2_1_2bn_w','fc2_1_2bn_b','fc2_1_2bn_m'});
net.addLayer('fc2_1_2x',dagnn.ReLU(),{'fc2_1_2bn'},{'fc2_1_2bnx'});

convBlock = dagnn.Conv('size',[1 1 128 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1_3',convBlock,{'fc2_1_2bnx'},{'fc2_1_3'},{'fc2_1_3f'});
net.addLayer('fc2_1_3bn',dagnn.BatchNorm(),{'fc2_1_3'},{'fc2_1_3bn'},...
    {'fc2_1_3bn_w','fc2_1_3bn_b','fc2_1_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 300 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1b',convBlock,{'fc2bn'},{'fc2_1b'},{'fc2_1bf'});
net.addLayer('fc2_1bbn',dagnn.BatchNorm(),{'fc2_1b'},{'fc2_1bbn'},...
    {'fc2_1bbn_w','fc2_1bbn_b','fc2_1bbn_m'});

net.addLayer('fc2_1sum',dagnn.Sum(),{'fc2_1_3bn','fc2_1bbn'},...
    {'fc2_1sum'});
net.addLayer('fc2_1x',dagnn.ReLU(),{'fc2_1sum'},{'fc2_1sumx'});

% 32*256
for i  = 2:3
    convBlock = dagnn.Conv('size',[1 1 256 64],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc2_%d_1',i),convBlock,{sprintf('fc2_%dsumx',i-1)},{sprintf('fc2_%d_1',i)}, ...
        {sprintf('fc2_%d_1f',i)});
    net.addLayer(sprintf('fc2_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc2_%d_1',i)},{sprintf('fc2_%d_1bn',i)},...
        {sprintf('fc2_%d_1bn_w',i),sprintf('fc2_%d_1bn_b',i),sprintf('fc2_%d_1bn_m',i)});
    net.addLayer(sprintf('fc2_%d_1x',i),dagnn.ReLU(),{sprintf('fc2_%d_1bn',i)},{sprintf('fc2_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 64 64],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc2_%d_2',i),convBlock,{sprintf('fc2_%d_1bnx',i)},{sprintf('fc2_%d_2',i)}, ...
        {sprintf('fc2_%d_2f',i)});
    net.addLayer(sprintf('fc2_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc2_%d_2',i)},{sprintf('fc2_%d_2bn',i)},...
        {sprintf('fc2_%d_2bn_w',i),sprintf('fc2_%d_2bn_b',i),sprintf('fc2_%d_2bn_m',i)});
    net.addLayer(sprintf('fc2_%d_2x',i),dagnn.ReLU(),{sprintf('fc2_%d_2bn',i)},{sprintf('fc2_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 64 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc2_%d_3',i),convBlock,{sprintf('fc2_%d_2bnx',i)},{sprintf('fc2_%d_3',i)}, ...
        {sprintf('fc2_%d_3f',i)});
    net.addLayer(sprintf('fc2_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc2_%d_3',i)},{sprintf('fc2_%d_3bn',i)},...
        {sprintf('fc2_%d_3bn_w',i),sprintf('fc2_%d_3bn_b',i),sprintf('fc2_%d_3bn_m',i)});

    net.addLayer(sprintf('fc2_%dsum',i),dagnn.Sum(),{sprintf('fc2_%dsumx',i-1),sprintf('fc2_%d_3bn',i)},...
        {sprintf('fc2_%dsum',i)});
    net.addLayer(sprintf('fc2_%dx',i),dagnn.ReLU(),{sprintf('fc2_%dsum',i)},{sprintf('fc2_%dsumx',i)});
end

%32*256
convBlock = dagnn.Conv('size',[1 1 256 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_4a_1',convBlock,{'fc2_3sumx'},{'fc2_4a_1'},{'fc2_4a_1f'});
net.addLayer('fc2_4a_1bn',dagnn.BatchNorm(),{'fc2_4a_1'},{'fc2_4a_1bn'},...
    {'fc2_4a_1bn_w','fc2_4a_1bn_b','fc2_4a_1bn_m'});
net.addLayer('fc2_4a_1x',dagnn.ReLU(),{'fc2_4a_1bn'},{'fc2_4a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 512 512],'hasBias',false,'stride',[2,2],'pad',[0,0,1,0]);
net.addLayer('fc2_4a_2',convBlock,{'fc2_4a_1bnx'},{'fc2_4a_2'},{'fc2_4a_2f'});
net.addLayer('fc2_4a_2bn',dagnn.BatchNorm(),{'fc2_4a_2'},{'fc2_4a_2bn'},...
    {'fc2_4a_2bn_w','fc2_4a_2bn_b','fc2_4a_2bn_m'});
net.addLayer('fc2_4a_2x',dagnn.ReLU(),{'fc2_4a_2bn'},{'fc2_4a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 512 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_4a_3',convBlock,{'fc2_4a_2bnx'},{'fc2_4a_3'},{'fc2_4a_3f'});
net.addLayer('fc2_4a_3bn',dagnn.BatchNorm(),{'fc2_4a_3'},{'fc2_4a_3bn'},...
    {'fc2_4a_3bn_w','fc2_4a_3bn_b','fc2_4a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 256 512],'hasBias',false,'stride',[2,2],'pad',[0,0,0,0]);
net.addLayer('fc2_4b',convBlock,{'fc2_3sumx'},{'fc2_4b'},{'fc2_4bf'});
net.addLayer('fc2_4bbn',dagnn.BatchNorm(),{'fc2_4b'},{'fc2_4bbn'},...
    {'fc2_4bbn_w','fc2_4bbn_b','fc2_4bbn_m'});

%16*512
net.addLayer('fc2_4sum',dagnn.Sum(),{'fc2_4a_3bn','fc2_4bbn'},...
    {'fc2_4sum'});
net.addLayer('fc2_4x',dagnn.ReLU(),{'fc2_4sum'},{'fc3_1sumx'});

for i  = 2:4
    convBlock = dagnn.Conv('size',[1 1 512 128],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc3_%d_1',i),convBlock,{sprintf('fc3_%dsumx',i-1)},{sprintf('fc3_%d_1',i)}, ...
        {sprintf('fc3_%d_1f',i)});
    net.addLayer(sprintf('fc3_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc3_%d_1',i)},{sprintf('fc3_%d_1bn',i)},...
        {sprintf('fc3_%d_1bn_w',i),sprintf('fc3_%d_1bn_b',i),sprintf('fc3_%d_1bn_m',i)});
    net.addLayer(sprintf('fc3_%d_1x',i),dagnn.ReLU(),{sprintf('fc3_%d_1bn',i)},{sprintf('fc3_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 128 128],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc3_%d_2',i),convBlock,{sprintf('fc3_%d_1bnx',i)},{sprintf('fc3_%d_2',i)}, ...
        {sprintf('fc3_%d_2f',i)});
    net.addLayer(sprintf('fc3_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc3_%d_2',i)},{sprintf('fc3_%d_2bn',i)},...
        {sprintf('fc3_%d_2bn_w',i),sprintf('fc3_%d_2bn_b',i),sprintf('fc3_%d_2bn_m',i)});
    net.addLayer(sprintf('fc3_%d_2x',i),dagnn.ReLU(),{sprintf('fc3_%d_2bn',i)},{sprintf('fc3_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 128 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc3_%d_3',i),convBlock,{sprintf('fc3_%d_2bnx',i)},{sprintf('fc3_%d_3',i)}, ...
        {sprintf('fc3_%d_3f',i)});
    net.addLayer(sprintf('fc3_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc3_%d_3',i)},{sprintf('fc3_%d_3bn',i)},...
        {sprintf('fc3_%d_3bn_w',i),sprintf('fc3_%d_3bn_b',i),sprintf('fc3_%d_3bn_m',i)});

    net.addLayer(sprintf('fc3_%dsum',i),dagnn.Sum(),{sprintf('fc3_%dsumx',i-1),sprintf('fc3_%d_3bn',i)},...
        {sprintf('fc3_%dsum',i)});
    net.addLayer(sprintf('fc3_%dx',i),dagnn.ReLU(),{sprintf('fc3_%dsum',i)},{sprintf('fc3_%dsumx',i)});
end

convBlock = dagnn.Conv('size',[1 1 512 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc3_5a_1',convBlock,{'fc3_4sumx'},{'fc3_5a_1'},{'fc3_5a_1f'});
net.addLayer('fc3_5a_1bn',dagnn.BatchNorm(),{'fc3_5a_1'},{'fc3_5a_1bn'},...
    {'fc3_5a_1bn_w','fc3_5a_1bn_b','fc3_5a_1bn_m'});
net.addLayer('fc3_5a_1x',dagnn.ReLU(),{'fc3_5a_1bn'},{'fc3_5a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 1024 1024],'hasBias',false,'stride',[2,2],'pad',[0,0,1,0]);
net.addLayer('fc3_5a_2',convBlock,{'fc3_5a_1bnx'},{'fc3_5a_2'},{'fc3_5a_2f'});
net.addLayer('fc3_5a_2bn',dagnn.BatchNorm(),{'fc3_5a_2'},{'fc3_5a_2bn'},...
    {'fc3_5a_2bn_w','fc3_5a_2bn_b','fc3_5a_2bn_m'});
net.addLayer('fc3_5a_2x',dagnn.ReLU(),{'fc3_5a_2bn'},{'fc3_5a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 1024 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc3_5a_3',convBlock,{'fc3_5a_2bnx'},{'fc3_5a_3'},{'fc3_5a_3f'});
net.addLayer('fc3_5a_3bn',dagnn.BatchNorm(),{'fc3_5a_3'},{'fc3_5a_3bn'},...
    {'fc3_5a_3bn_w','fc3_5a_3bn_b','fc3_5a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 512 1024],'hasBias',false,'stride',[2,2],'pad',[0,0,0,0]);
net.addLayer('fc3_5b',convBlock,{'fc3_4sumx'},{'fc3_5b'},{'fc3_5bf'});
net.addLayer('fc3_5bbn',dagnn.BatchNorm(),{'fc3_5b'},{'fc3_5bbn'},...
    {'fc3_5bbn_w','fc3_5bbn_b','fc3_5bbn_m'});

%8*1024
net.addLayer('fc3_5sum',dagnn.Sum(),{'fc3_5a_3bn','fc3_5bbn'},...
    {'fc3_5sum'});
net.addLayer('fc3_5x',dagnn.ReLU(),{'fc3_5sum'},{'fc4_1sumx'});

for i  = 2:6
    convBlock = dagnn.Conv('size',[1 1 1024 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc4_%d_1',i),convBlock,{sprintf('fc4_%dsumx',i-1)},{sprintf('fc4_%d_1',i)}, ...
        {sprintf('fc4_%d_1f',i)});
    net.addLayer(sprintf('fc4_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc4_%d_1',i)},{sprintf('fc4_%d_1bn',i)},...
        {sprintf('fc4_%d_1bn_w',i),sprintf('fc4_%d_1bn_b',i),sprintf('fc4_%d_1bn_m',i)});
    net.addLayer(sprintf('fc4_%d_1x',i),dagnn.ReLU(),{sprintf('fc4_%d_1bn',i)},{sprintf('fc4_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 256 256],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc4_%d_2',i),convBlock,{sprintf('fc4_%d_1bnx',i)},{sprintf('fc4_%d_2',i)}, ...
        {sprintf('fc4_%d_2f',i)});
    net.addLayer(sprintf('fc4_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc4_%d_2',i)},{sprintf('fc4_%d_2bn',i)},...
        {sprintf('fc4_%d_2bn_w',i),sprintf('fc4_%d_2bn_b',i),sprintf('fc4_%d_2bn_m',i)});
    net.addLayer(sprintf('fc4_%d_2x',i),dagnn.ReLU(),{sprintf('fc4_%d_2bn',i)},{sprintf('fc4_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 256 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc4_%d_3',i),convBlock,{sprintf('fc4_%d_2bnx',i)},{sprintf('fc4_%d_3',i)}, ...
        {sprintf('fc4_%d_3f',i)});
    net.addLayer(sprintf('fc4_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc4_%d_3',i)},{sprintf('fc4_%d_3bn',i)},...
        {sprintf('fc4_%d_3bn_w',i),sprintf('fc4_%d_3bn_b',i),sprintf('fc4_%d_3bn_m',i)});

    net.addLayer(sprintf('fc4_%dsum',i),dagnn.Sum(),{sprintf('fc4_%dsumx',i-1),sprintf('fc4_%d_3bn',i)},...
        {sprintf('fc4_%dsum',i)});
    net.addLayer(sprintf('fc4_%dx',i),dagnn.ReLU(),{sprintf('fc4_%dsum',i)},{sprintf('fc4_%dsumx',i)});
end

convBlock = dagnn.Conv('size',[1 1 1024 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7a_1',convBlock,{'fc4_6sumx'},{'fc4_7a_1'},{'fc4_7a_1f'});
net.addLayer('fc4_7a_1bn',dagnn.BatchNorm(),{'fc4_7a_1'},{'fc4_7a_1bn'},...
    {'fc4_7a_1bn_w','fc4_7a_1bn_b','fc4_7a_1bn_m'});
net.addLayer('fc4_7a_1x',dagnn.ReLU(),{'fc4_7a_1bn'},{'fc4_7a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
net.addLayer('fc4_7a_2',convBlock,{'fc4_7a_1bnx'},{'fc4_7a_2'},{'fc4_7a_2f'});
net.addLayer('fc4_7a_2bn',dagnn.BatchNorm(),{'fc4_7a_2'},{'fc4_7a_2bn'},...
    {'fc4_7a_2bn_w','fc4_7a_2bn_b','fc4_7a_2bn_m'});
net.addLayer('fc4_7a_2x',dagnn.ReLU(),{'fc4_7a_2bn'},{'fc4_7a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7a_3',convBlock,{'fc4_7a_2bnx'},{'fc4_7a_3'},{'fc4_7a_3f'});
net.addLayer('fc4_7a_3bn',dagnn.BatchNorm(),{'fc4_7a_3'},{'fc4_7a_3bn'},...
    {'fc4_7a_3bn_w','fc4_7a_3bn_b','fc4_7a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 1024 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7b',convBlock,{'fc4_6sumx'},{'fc4_7b'},{'fc4_7bf'});
net.addLayer('fc4_7bbn',dagnn.BatchNorm(),{'fc4_7b'},{'fc4_7bbn'},...
    {'fc4_7bbn_w','fc4_7bbn_b','fc4_7bbn_m'});

%8*2048
net.addLayer('fc4_7sum',dagnn.Sum(),{'fc4_7a_3bn','fc4_7bbn'},...
    {'fc4_7sum'});
net.addLayer('fc4_7x',dagnn.ReLU(),{'fc4_7sum'},{'fc5_1sumx'});

poolBlock = dagnn.Pooling('poolSize',[1 8]);
net.addLayer('fc5_1',poolBlock,{'fc5_1sumx'},{'fc5_1bnx'});
%net.addLayer('fc5_1bn',dagnn.BatchNorm(),{'fc5_1'},{'fc5_1bn'},...
 %   {'fc5_1bn_w','fc5_1bn_b','fc5_1bn_m'});
%net.addLayer('fc5_1x',dagnn.ReLU(),{'fc5_1bn'},{'fc5_1bnx'});

fc5_2Block = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc5_2',fc5_2Block,{'fc5_1bnx'},{'fc5_2'},{'fc5_2f'});
net.addLayer('fc5_2bn',dagnn.BatchNorm(),{'fc5_2'},{'fc5_2bn'},...
    {'fc5_2bn_w','fc5_2bn_b','fc5_2bn_m'});
net.addLayer('fc5_2x',dagnn.ReLU(),{'fc5_2bn'},{'fc5_2bnx'});
net.addLayer('dropout2',dagnn.DropOut('rate',0.5),{'fc5_2bnx'},{'fc5_2bnxd'});

%----------------------add share layer----------------------
%1

fc_imgBlock = dagnn.Conv('size',[1 1 2048 113287],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc_img',fc_imgBlock,{'fc1_1bnxd'},{'prediction_img'},{'fcsharef'});
net.addLayer('softmaxloss_img',dagnn.Loss('loss','softmaxlog'),{'prediction_img','label_img'},'objective_img');
net.addLayer('top1err_img', dagnn.Loss('loss', 'classerror'), ...
    {'prediction_img','label_img'}, 'top1err_img') ;
net.addLayer('top5err_img', dagnn.Loss('loss', 'topkerror', ...
    'opts', {'topK',5}), ...
    {'prediction_img','label_img'}, 'top5err_img') ;
%2
fc_txtBlock = dagnn.Conv('size',[1 1 2048 113287],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc_txt',fc_txtBlock,{'fc5_2bnxd'},{'prediction_txt'},{'fcsharef'});
net.addLayer('softmaxloss_txt',dagnn.Loss('loss','softmaxlog'),{'prediction_txt','label_txt'},'objective_txt');
net.addLayer('top1err_txt', dagnn.Loss('loss', 'classerror'), ...
    {'prediction_txt','label_txt'}, 'top1err_txt') ;
net.addLayer('top5err_txt', dagnn.Loss('loss', 'topkerror', ...
    'opts', {'topK',5}), ...
    {'prediction_txt','label_txt'}, 'top5err_txt') ;
%}

%}
net.initParams();
%----------NOTICE--------------
first = net.getParamIndex('fc2f');

net.params(first).learningRate = 1e-3;  %w
net.params(first+1).learningRate = 1e-3;  %b

load('./dataset/MSCOCO-prepare/COCO_dictionary.mat');
net.params(first).value = reshape(single(subset.features'),1,1,29972,300);

end

author: layumi

repo: Image-Text-Embedding


text:
% This example shows how to find the modes of a birefringent
% (uniaxial) waveguide.  This waveguide is the same as the one
% considered in 'uniaxial_channel.m', except that the c-axis
% is now rotated by an angle of -pi/4 relative to the x axis.

n1 = 1.55;
n2x = 2.156;
n2y = 2.232;
n2z = 2.232;
theta = -pi/4;

e2xx = n2x^2*cos(theta)^2 + n2y^2*sin(theta)^2;
e2yy = n2y^2*cos(theta)^2 + n2x^2*sin(theta)^2;
e2xy = cos(theta)*sin(theta)*(n2x^2-n2y^2);
e2yx = e2xy;

Rx = 0.30;
Ry = 0.20;
side = 0.2;

dx = 0.0025;        % grid size (x)
dy = dx;            % grid size (y)

lambda = 1.00;      % wavelength
nmodes = 2;         % number of modes to compute

[x,y,xc,yc,nx,ny,epsxx,edges] = ...
    waveguidemeshfull([n1,sqrt(e2xx),n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epsxy,edges] = ...
    waveguidemeshfull([0,sqrt(e2xy),0],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epsyx,edges] = ...
    waveguidemeshfull([0,sqrt(e2yx),0],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epsyy,edges] = ...
    waveguidemeshfull([n1,sqrt(e2yy),n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epszz,edges] = ...
    waveguidemeshfull([n1,n2z,n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);

% Now we stretch out the mesh at the boundaries:
[x,y,xc,yc,dx,dy] = stretchmesh(x,y,[80,80,80,80],[4,4,4,4]);

[Hx,Hy,neff] = wgmodes (lambda, n2y, nmodes, dx, dy, epsxx, epsxy, epsyx, epsyy, epszz, '0000');

fprintf(1,'neff = %7.5f\n',neff);

figure(1);

for ii = 1:nmodes,
  subplot(nmodes,2,2*(ii-1)+1);
  contourmode(x,y,Hx(:,:,ii));
  title(sprintf('Hx (mode %d)',ii));
  for v = edges, line(v{:}); end
  subplot(nmodes,2,2*(ii-1)+2);
  contourmode(x,y,Hy(:,:,ii));
  title(sprintf('Hy (mode %d)',ii));
  for v = edges, line(v{:}); end
end

author: Sable

repo: mcbench-benchmarks


text:
function Cd = som_cldist(D,clinds1,clinds2,cldist,q,mask)

% SOM_CLDIST Distances between two clusters.
%
%   Cd = som_cldist(Md,c1,c2,'single')
%   Cd = som_cldist(Md,c1,c2,'average')
%   Cd = som_cldist(Md,c1,c2,'complete')
%   Cd = som_cldist(Md,c1,c2,'neighf',H)
%   Cd = som_cldist(Md,c1,[],...)
%   Cd = som_cldist(D,c1,c2,'centroid',q,mask)
%   Cd = som_cldist(D,c1,c2,'ward',q,mask)
%   Cd = som_cldist(D,c1,[],...)
%
%  Input and output arguments ([]'s are optional):
%   D        (matrix) size dlen x dim, the data set
%            (struct) map or data struct
%   Md       (matrix) size dlen x dlen, mutual distance matrix, see SOM_MDIST
%   c1       (cell array) size n1 x 1, indices of clusters from which
%                     the distances should be calculated, each cell
%                     contains indices of vectors that belong to that
%                     cluster (indices are between 1...dlen)
%   c2       (cell array) size n2 x 1, same as c1 but have the clusters
%                     to which the distances should be calculated
%            (empty)  c1 is used in place of c2
%   [q]      (scalar) distance norm, default = 2
%   [mask]   (vector) size dim x 1, the weighting mask, a vector of ones
%                     by default
%   H        (matrix) size dlen x dlen, neighborhood function values
%
%   Cd       (matrix) size n1 x n2, distances between the clusters
%
% See also SOM_MDIST.

% Copyright (c) 2000 by Juha Vesanto
% Contributed to SOM Toolbox on XXX by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 250800

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

[dlen dim] = size(D);
if nargin<5, q = 2; end
if nargin<6, mask = ones(dim,1); end
if ~iscell(clinds1), clinds1 = {clinds1}; end
if ~isempty(clinds2) && ~iscell(clinds2), clinds2 = {clinds2}; end

n1 = length(clinds1);
n2 = length(clinds2);
if n2>0, Cd = zeros(n1,n2); else Cd = zeros(n1); end
if n1==0, return; end

switch cldist,

 % centroid distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 'centroid',

  C1 = zeros(n1,dim); for i=1:n1, C1(i,:) = mean(D(clinds1{i},:),1); end
  C2 = zeros(n2,dim); for i=1:n2, C2(i,:) = mean(D(clinds2{i},:),1); end
  if n2==0,
    for i=1:n1-1,
      for j=i+1:n1,
	diff = C1(i,:)-C1(j,:);
	switch q,
	 case 1,    Cd(i,j)=abs(diff)*mask;
	 case 2,    Cd(i,j)=sqrt((diff.^2)*mask);
	 case Inf,  Cd(i,j)=max(diag(mask)*abs(diff),[],2);
	 otherwise, Cd(i,j)=((abs(diff).^q)*mask).^(1/q);
	end
      end
      Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
    end
  else
    for i=1:n1,
      for j=1:n2,
	diff = C1(i,:)-C2(j,:);
	switch q,
	 case 1,    Cd(i,j)=abs(diff)*mask;
	 case 2,    Cd(i,j)=sqrt((diff.^2)*mask);
	 case Inf,  Cd(i,j)=max(diag(mask)*abs(diff),[],2);
	 otherwise, Cd(i,j)=((abs(diff).^q)*mask).^(1/q);
	end
      end
    end
  end

 % ward distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 'ward',

  C1 = zeros(n1,dim); nn1 = zeros(n1,dim);
  for i=1:n1, C1(i,:) = mean(D(clinds1{i},:),1); nn1(i) = length(clinds1{i}); end
  C2 = zeros(n2,dim); nn2 = zeros(n2,dim);
  for i=1:n2, C2(i,:) = mean(D(clinds2{i},:),1); nn2(i) = length(clinds2{i}); end
  if n2==0,
    for i=1:n1-1,
      for j=i+1:n1,
	diff = C1(i,:) - C1(j,:);
	f = 2*nn1(i)*nn1(j) / (nn1(i)+nn1(j));
	switch q,
	 case 1,    Cd(i,j)=f*abs(diff)*mask;
	 case 2,    Cd(i,j)=f*sqrt((diff.^2)*mask);
	 case Inf,  Cd(i,j)=f*max(diag(mask)*abs(diff),[],2);
	 otherwise, Cd(i,j)=f*((abs(diff).^q)*mask).^(1/q);
	end
      end
      Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
    end
  else
    for i=1:n1,
      for j=1:n2,
	diff = C1(i,:) - C2(j,:);
	f = 2*nn1(i)*nn2(j) / (nn1(i)+nn2(j));
	switch q,
	 case 1,    Cd(i,j)=f*abs(diff)*mask;
	 case 2,    Cd(i,j)=f*sqrt((diff.^2)*mask);
	 case Inf,  Cd(i,j)=f*max(diag(mask)*abs(diff),[],2);
	 otherwise, Cd(i,j)=f*((abs(diff).^q)*mask).^(1/q);
	end
      end
    end
  end

 % single linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 'single',

  if n2==0,
    for i=1:n1-1,
      for j=i+1:n1,
	vd = D(clinds1{i},clinds1{j});
	fi = isfinite(vd(:));
	if any(fi), Cd(i,j) = min(vd(fi)); else Cd(i,j) = Inf; end
      end
      Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
    end
  else
    for i=1:n1,
      for j=1:n2,
	vd = D(clinds1{i},clinds2{j});
	fi = isfinite(vd(:));
	if any(fi), Cd(i,j) = min(vd(fi)); else Cd(i,j) = Inf; end
      end
    end
  end

 % average linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 'average',

  if n2==0,
    for i=1:n1-1,
      for j=i+1:n1,
	vd = D(clinds1{i},clinds1{j});
	fi = isfinite(vd(:));
	if any(fi), Cd(i,j) = mean(vd(fi)); else Cd(i,j) = Inf; end
      end
      Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
    end
  else
    for i=1:n1,
      for j=1:n2,
	vd = D(clinds1{i},clinds2{j});
	fi = isfinite(vd(:));
	if any(fi), Cd(i,j) = mean(vd(fi)); else Cd(i,j) = Inf; end
      end
    end
  end

 % complete linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 'complete',

   if n2==0,
    for i=1:n1-1,
      for j=i+1:n1,
	vd = D(clinds1{i},clinds1{j});
	fi = isfinite(vd(:));
	if any(fi), Cd(i,j) = max(vd(fi)); else Cd(i,j) = Inf; end
      end
      Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
    end
  else
    for i=1:n1,
      for j=1:n2,
	vd = D(clinds1{i},clinds2{j});
	fi = isfinite(vd(:));
	if any(fi), Cd(i,j) = max(vd(fi)); else Cd(i,j) = Inf; end
      end
    end
  end

 % neighborhood function linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 'neighf',

  if n2==0,
    for i=1:n1-1,
      for j=i+1:n1,
	vd = D(clinds1{i},clinds1{j});
	fi = isfinite(vd(:));
	if any(fi),
	  hd = q(clinds1{i},clinds1{j});
	  hd = hd(fi);
	  Cd(i,j) = sum(hd.*vd(fi))/sum(hd);
	else Cd(i,j) = Inf;
	end
      end
      Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
    end
  else
    for i=1:n1,
      for j=1:n2,
	vd = D(clinds1{i},clinds2{j});
	fi = isfinite(vd(:));
	if any(fi),
	  hd = q(clinds1{i},clinds2{j});
	  hd = hd(fi);
	  Cd(i,j) = sum(hd.*vd(fi))/sum(hd);
	else Cd(i,j) = Inf;
	end
      end
    end
  end

 otherwise, error(['Unknown cluster distance metric: ' cldist]);
end

return;

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

author: ilarinieminen

repo: SOM-Toolbox


text:
classdef BT5 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP with bias feature

%------------------------------- Reference --------------------------------
% H. Li, Q. Zhang, and J. Deng, Biased multiobjective optimization and
% decomposition algorithm, IEEE Transactions on Cybernetics, 2017, 47(1):
% 52-66.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            I1    = 2 : 2 : D;
            I2    = 3 : 2 : D;
            Y     = X - sin(repmat(1:D,N,1)*pi/2/D);
            PopObj(:,1) = X(:,1)                                     + sum(Y(:,I1).^2+(1-exp(-Y(:,I1).^2/1e-10))/5,2);
            PopObj(:,2) = (1-X(:,1)).*(1-X(:,1).*sin(8.5*pi*X(:,1))) + sum(Y(:,I2).^2+(1-exp(-Y(:,I2).^2/1e-10))/5,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = (1-R(:,1)).*(1-R(:,1).*sin(8.5*pi*R(:,1)));
            R      = R(NDSort(R,1)==1,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,500)';
            R(:,2) = (1-R(:,1)).*(1-R(:,1).*sin(8.5*pi*R(:,1)));
            R(NDSort(R,1)~=1,:) = nan;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function model = ivmOptimise(model, options);

% IVMOPTIMISE Optimise the IVM.
% FORMAT
% DESC optimises an IVM by iterating between selecting points and
% optimising kernel and noise parameters.
% ARG model : the model to be optimised.
% ARG options : options structure as returned by ivmOptions.
%
% SEEALSO : ivmOptimiseIvm, ivmOptimiseKernel, ivmOptimiseNoise
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% IVM

for i = 1:options.extIters
  if options.kernIters
    % Update the kernel if required.
    model = ivmOptimiseIvm(model, options.display);
    model = ivmOptimiseKernel(model, options.display, options.kernIters);
  end
  if options.noiseIters
    % Update the noise model if required.
    model = ivmOptimiseIvm(model, options.display);
    model = ivmOptimiseNoise(model, options.display, options.noiseIters);
  end
  if options.display
    ivmDisplay(model);
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [output] = minus(input1, input2)

if ~iscell(input1)
  input1 = repmat({input1}, size(input2));
end
if ~iscell(input2)
  input2 = repmat({input2}, size(input1));
end
output = cellfun(@minus,input1,input2,'uniformoutput',false);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ y, m, d, f, ierror ] = ymdf_check_republican ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_CHECK_REPUBLICAN checks a Republican YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, real F, the
%    YMDF date, which may be corrected if necessary and possible.
%
%    Output, integer IERROR, is 0 if the date is legal.
%
  ierror = 0;

  [ y, m, d, ierror ] = ymd_check_republican ( y, m, d );

  if ( ierror ~= 0 )
    return
  end

  [ y, m, d, f ] = frac_borrow_republican ( y, m, d, f );

  [ y, m, d, f ] = frac_carry_republican ( y, m, d, f );

  return
end


text:
function c8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% C8MAT_PRINT_SOME prints some of a C8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the matrix.
%
%    Input, complex A(M,N), the matrix.
%
%    Input, integer ILO, JLO, IHI, JHI, the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  incx = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
%
%  Print the columns of the matrix, in strips of INCX.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );
    for j = j2lo : j2hi
      j2 = j + 1 - j2lo;
      fprintf ( 1, '%10d          ', j );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi
%
%  Print out (up to) INCX entries in row I, that lie in the current strip.
%
      for j2 = 1 : inc

        j = j2lo - 1 + j2;

        if ( a(i,j) == 0.0 )
          fprintf ( 1, '       0.0          ' );
        elseif ( imag ( a(i,j) ) == 0.0 )
          fprintf ( 1, ' %9f          ', real ( a(i,j) ) );
        else
          fprintf ( 1, ' %9f %9f', real ( a(i,j) ), imag ( a(i,j) ) );
       end

      end

      fprintf ( 1, '\n' );

    end

  end

  fprintf ( 1, '\n' );

  return
end


text:
function  [a,resnorm]=Optimize_my_LM2(Loss_fun,a0,data,TolX,TolFun,MaxIter)

% author  Zhang Xin

tao=1e-10;
xk=a0;

v=2;

Jacobi=Get_Jacobi(Loss_fun,xk,data);
Ek=Loss_fun(xk,data);

g=Jacobi'*Ek;

found=logical(norm(g)<=TolFun);
mou=tao*max(diag(Jacobi'*Jacobi));

k=0;
while (~found && k<MaxIter+1)

    %delta_x=-(Jacobi'*Jacobi+mou*eye(size(a0,2)))\Jacobi'*Ek;

    delta_x=-[Jacobi;sqrt(mou)*eye(size(a0,2))]\[Ek;zeros(size(a0,2),1)];

    if (norm(delta_x)<=TolX*(norm(xk)+TolX))
        found=true;

    else
        xk_new=xk+delta_x';
        Ek=Loss_fun(xk,data);
        Ek_new=Loss_fun(xk_new,data);
        L0=delta_x'*Jacobi'*Ek;
        L_delta=delta_x'*Jacobi'*Jacobi*delta_x;
        rho=(Ek'*Ek-Ek_new'*Ek_new)/(-L0-L_delta);

        if rho>0

            fprintf('Iterations: %d, Residual: %d, Step: %d \n',k,Ek'*Ek,norm(delta_x));
            k=k+1;

            xk=xk_new;
            Jacobi=Get_Jacobi(Loss_fun,xk,data);
            Ek=Loss_fun(xk,data);
            g=Jacobi'*Ek;
            found=(norm(g)<=TolFun);
            mou=mou*max([1/3,1-(2*rho-1)^3]);
            v=2;
        else
            mou=mou*v;
            v=2*v;

        end

    end

end

xk=xk+delta_x';
Ek=Loss_fun(xk,data);
fprintf('Iterations: %d, Residual: %d, Step: %d \n',k,Ek'*Ek,norm(delta_x));

a=xk;

resnorm=Ek'*Ek;

end

function Jacobi=Get_Jacobi(Loss_fun,xk,data)

scale=1e-4;

%Ek=Loss_fun(xk,data);

for i=1:length(xk)
    x_temp1=xk;
    x_temp2=xk;
    if abs(x_temp1(i))>scale

        delta=x_temp1(i)*scale;
    else
        delta=scale;
    end
    x_temp1(i)=x_temp1(i)+delta;
    x_temp2(i)=x_temp2(i)-delta;

    E_temp1=Loss_fun(x_temp1,data);

    E_temp2=Loss_fun(x_temp2,data);

    Jacobi(:,i)=(E_temp1-E_temp2)/delta/2;

end

end

author: shenshikexmu

repo: IMUCalibration-Gesture


text:
% PURPOSE: sets ERPLAB Foreground Color
%
% FORMAT
%
% Fcolorerplab(c)
%
% INPUT
%
% c      - color. e.g. [0 1 0] -> green
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2011

function Fcolorerplab(c)
if nargin<1
        c = uisetcolor([0.83 0.82 0.79],'ERPLAB Foreground Color') ;
end
erpworkingmemory('ColorF', c);

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function output=user7sym()

%**************************************************************************
 %*************************************************************************
                %//input signal form transmitter side for symbol based\\
 clear all;
close all;
clc;
% us=input('enter the number of users');
% nb1=input('number of bits for user1');
% m1=input('enter the message sequence1');
% nb2=input('number of bits for user2');
% m2=input('enter the message sequence2');
% nb3=input('number of bits for user3');
% m3=input('enter the message sequence3');
% nb4=input('number of bits for user4');
% m4=input('enter the message sequence4');
% nb5=input('number of bits for user5');
% m5=input('enter the message sequence5');
us=3;
nb1=20;
nb2=20;
nb3=20;
nb4=20;
nb5=20;
nb6=20;
nb7=20;
m1=[1 1 1 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 0];
m2=[1 1 0 1 1 1 1 0 1 1 1 0 1 1 0 1 1 1 1 0 ];
m3=[1 0 1 1 0 1 1 0 1 1 1 0 1 0 1 1 0 1 1 0];
m4=[1 1 1 1 0 1 1 1 0 1 1 0 1 1 1 1 1 1 0 1 ];
m5=[1 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 0 1 1 1 ];
m6=[1 0 1 1 0 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 ];
m7=[1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1];

% nb1=20;
% nb2=20;
% nb3=20;
% nb4=20;
% nb5=20;
% nb6=20;
% nb7=20;
% m1=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
% m2=[1 0 1 0 1  0 1 1 0 1 1 1 1 0 1 0 1 0 1 0];
% m3=[1 0 1 0 1 1 0 0 1 1 1 0 1 0 1 0 1 0 1 1 ];
% m4=[1 1 0 1 0  1 0 1 0 1 1 1 1 0 1 1 0 1 0 1];
% m5=[1 0 1 1 0  0 1 1 0 1 1 0 1 1 1 0 0 1 0 0];
% m6=[1 0 1 0 1  1 0 1 0 1 1 1 0 0 1 0 1 0 1 1];
% m7=[1 1 0 1 1 0 1 1 1 0 1 0 1 0 1 1 0 1 1 0];

                    %//message1 into polar form\\

for i=1:nb1
    if m1(i)==1
        mb1(i)=m1(i);
    else
        mb1(i)=-1;
    end
end
  % disp(mb1);

                %//message2 into polar form\\

for i=1:nb2
    if m2(i)==1
        mb2(i)=m2(i);
    else
        mb2(i)=-1;
    end
end
  %disp(mb2);

                %//message3 into polar form\\

for i=1:nb3
    if m3(i)==1
        mb3(i)=m3(i);
    else
        mb3(i)=-1;
    end
end
  %disp(mb3);

       %//message4 into polar form\\

for i=1:nb4
    if m4(i)==1
        mb4(i)=m4(i);
    else
        mb4(i)=-1;
    end
end
  %disp(mb4);

    %//message5 into polar form\\

for i=1:nb5
    if m5(i)==1
        mb5(i)=m5(i);
    else
        mb5(i)=-1;
    end
end
  %disp(mb5);

       %//message6 into polar form\\

for i=1:nb6
    if m6(i)==1
        mb6(i)=m6(i);
    else
        mb6(i)=-1;
    end
end
  %disp(mb6);

    %//message7 into polar form\\

for i=1:nb7
    if m7(i)==1
        mb7(i)=m7(i);
    else
        mb7(i)=-1;
    end
end
  %disp(mb7);

                  %//generation of maximal length sequence1\\

  f1=1;f2=1;f3=0;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f2,f5);
    f5=f4;
    op1(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op1);

                    %//msequence1 into polar form\\

    for i=1:31

         if op1(i)==1
             pb1(i)=1;
         else
         pb1(i)=-1;
        end
    end
    %disp(op1);
    %disp(pb1);

                    %//generation of maximal length sequence2\\

  f1=1;f2=0;f3=0;f4=1;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f3,f5);
    f5=f4;
    op2(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op2);

                    %//msequence2 into polar form\\

    for i=1:31
        if op2(i)==1
             pb2(i)=1;
         else
         pb2(i)=-1;
        end
    end
    %disp(op2);
    %disp(pb2);

                     %//generation of maximal length sequence3\\

  f1=1;f2=1;f3=1;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op3(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op3);

                    %//msequence3 into polar form\\

    for i=1:31

         if op3(i)==1
             pb3(i)=1;
         else
         pb3(i)=-1;
        end
    end
    %disp(op3);
    %disp(pb3);

                     %//generation of maximal length sequence4\\

  f1=0;f2=1;f3=1;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op4(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op4);

                    %//msequence4 into polar form\\

    for i=1:31

         if op4(i)==1
             pb4(i)=1;
         else
         pb4(i)=-1;
        end
    end
    %disp(op4);
    %disp(pb4);

                     %//generation of maximal length sequence5\\

  f1=1;f2=1;f3=0;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op5(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op5);

                    %//msequence5 into polar form\\

    for i=1:31

         if op5(i)==1
             pb5(i)=1;
         else
         pb5(i)=-1;
        end
    end
    %disp(op5);
    %disp(pb5);

                     %//generation of maximal length sequence6\\

  f1=0;f2=0;f3=1;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op6(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op6);

                            %//msequence6 into polar form\\

    for i=1:31

         if op6(i)==1
             pb6(i)=1;
         else
         pb6(i)=-1;
        end
    end
    %disp(op6);
    %disp(pb6);

                 %//generation of maximal length sequence7\\

  f1=1;f2=1;f3=0;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op7(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op7);

                      %//msequence7 into polar form\\

    for i=1:31

         if op7(i)==1
             pb7(i)=1;
         else
         pb7(i)=-1;
        end
    end
    %disp(op7);
    %disp(pb7);

                    %//spreading signals from transmitter side\\

                         %//first transmit bit\\

k=1;
for i=1:nb1
    for j=1:31
        tb1(k)=mb1(i)*pb1(j);
        k=k+1;
    end
end
%disp(tb1);

                     %//second transmit bit\\

k=1;
for i=1:nb2
    for j=1:31
        tb2(k)=mb2(i)*pb2(j);
        k=k+1;
    end
end
%disp(tb2);

                    %//third transmit bit\\

k=1;
for i=1:nb3
    for j=1:31
        tb3(k)=mb3(i)*pb3(j);
        k=k+1;
    end
end
%disp(tb3);

                %//fourth transmit bit\\

k=1;
for i=1:nb4
    for j=1:31
        tb4(k)=mb4(i)*pb4(j);
        k=k+1;
    end
end
%disp(tb4);

             %//fifth transmit bit\\

k=1;
for i=1:nb5
    for j=1:31
        tb5(k)=mb5(i)*pb5(j);
        k=k+1;
    end
end
%disp(tb5);

     %//sixth transmit bit\\

k=1;
for i=1:nb6
    for j=1:31
        tb6(k)=mb6(i)*pb6(j);
        k=k+1;
    end
end
%disp(tb6);

     %//seventh transmit bit\\

k=1;
for i=1:nb7
    for j=1:31
        tb7(k)=mb7(i)*pb7(j);
        k=k+1;
    end
end
%disp(tb7);

                 %//Addition of seven  signals transmitted in the channel\\

n=1;
for i=1:k-1
    tb(n)=tb1(n)+tb2(n)+tb3(n)+tb4(n)+tb5(n)+tb6(n)+tb7(n);
    n=n+1;
end
%disp(tb);
chn=awgn(tb,10);

                         %//%receiver side\\
                    %//reconstruction of user2 signal\\

s2=0;
t2=1;
for i=0:31:k-32
    for j=1:31
        ob2(t2)=chn(i+j)*pb2(j)+s2;
        s2=ob2(t2);
        if ob2(t2)>0
            ob2(t2)=1;
        else ob2(t2)=0;
        end
    end
    t2=t2+1;
    s2=0;
end
%disp('second message is')
%disp(ob2);
n1=awgn(ob2,10);
 figure(1);
grid on;
subplot(2,2,1);
plot(ob2);
title('Received signal with out noise for SB');
xlabel('Time');
ylabel('Amplitude');
%g=awgn(ob2,10);
%disp(g);
subplot(2,2,2);
plot(chn);
title('Received signal with noise for SB');
xlabel('Time');
ylabel('Amplitude');

                    %//reconstruction of received signal from noise\\

for i=1:nb2
     if chn(i)>0
            g2(i)=1;
        else g2(i)=0;
        end
    end

% disp(g2);

                 % //calculating optimum weight using minimum variance\\

 ob=randint(1,31);
q=length(10);
%q1=input('enter the sequence');
s=randsrc(1,q);
c1=complex(ob,s);
%disp(c);
d=conj(c1);
d1=d*d';
%disp(d1);
m=mean(d1);
lamda=5;
teta=45;
k=((2*pi)/lamda);
K=4;
for i=0:K-1
    v=(exp(j*k*i*d*sin(teta)))';
end
%disp(v);
v1=conj(v)';
teta=30;
k=((2*pi)/lamda);
K=5;
for i=0:K-1
    eta=exp(sqrt(-1)*k*i*d*sin(teta))';
end
%u1=[2+2j 3-1j 4+3j 5+2j]'
%m=randsrc(1,124);
%a=randsrc(1,124);
l=length(10);
%l1=input('enter the sequence');
u1=complex(ob2,l);
u=eta*u1;
u2=conj(u)';
u3=u*u2;
Ru=mean(u3);

%Ru1=imresize(Ru,[124 124]);
%a1=inv(Ru1);
%a1=Ru1';
a1=Ru';
%disp(a1);
%calculating beta value
%g=1;
dem1=v*v1*a1;
beta1=dem1.^-1;
 %optimum weight
%Wopt=beta*a1*v';
temp1=beta1*v';
Wopt1=temp1*a1;
%disp(Wopt1);
 subplot(2,2,3);
plot(Wopt1);
title('Received signal with optimum weight for SB');
xlabel('Time');
ylabel('weight');

                    %//beamformer output\\

o1=conj(Wopt1)*ob2;
%disp('Beamformer output for SB')
%disp(o1);
 subplot(2,2,4);
plot(o1);
title('Beamformer output  for SB');
xlabel('Time');
ylabel('Beam former output');

                    %//SINR CALCULATION FOR SB CONFIGURATION\\

%s=length(m2);
%q1=input('enter the sequence');
s=randsrc(1,q);
c=complex(ob2,s);

hop1=conj(pb2)';
hob1=conj(c)';
nr3=hop1*hob1';
sum2=nr3'*Wopt1;
%disp('sum2');
%disp(sum2);
cm2=sum2'*sum2;
signal1=mean(cm2);
%disp('cm2');
%disp(cm2);
hop2=conj(pb2)';
hx2=conj(n1)';
nr5=hop2*hx2';
 sum3=nr5'*Wopt1;
 %disp('sum');
 %disp(sum);
  cm3=sum3'*sum3;
  noise1=mean(cm3);
 %disp('cm3');
 %disp(cm3);
 SSINR=signal1/noise1;
 disp('SSINR');
 disp(SSINR);

                %//GRAPH FOR DOA VERSUS SINR\\

    DOA=-80:10:80;
    figure(4);
    plot(DOA,SSINR,'B-*');
    title('simulation for DOA versus SSINR');
    xlabel('DOA in degree');
    ylabel('SINR in db');

                %// BER CALCULATION FOR SB CONFIGURATION\\

       k1=biterr(m2,g2);
       output=k1;
       disp('biterror rate for Sb configuration');
       disp(k1);

                 %//GRAPH FOR NUMBER USERS VERSUS BER\\

     figure(5);
     plot(us,k1,'R-*');

     title('simulation for No. of users versus BER SYMBOL BASED ');
     xlabel('NUMBER OF USERS');
     ylabel('BER');


text:
function [v1,v2,v3,v4,v5,v6,v7]=lsqnonlinSym(funhandle,varargin)
% LSQNONLINSYM   A wrap of lsqnonlin to use symbolic toolbox to obtain Jacobian.
% Limitation: The cost function and nonlinear constraint function cannot
% have conditional operation, such as if, max, min etc.
% Its usage is the same as the original lsqnonlin.
%
% Example: see lsqnonlinSym_example
%
% By Yi Cao, Cranfield University, 25/09/2007
%
x0=varargin{1};
n=numel(x0);

x=[];
for k=1:n
    eval(['syms x' num2str(k)]);
    eval(['x=[x;x' num2str(k) '];']);
end

nv=numel(varargin);
if nv>4
    F=feval(funhandle,x,varargin(5:end));
else
    F=feval(funhandle,x);
end
F=simplify(F);
Fx=simplify(jacobian(F,x));

if nv>=4
    opt = varargin{4};
    opt = optimset(opt,'jacobian','on');
else
    opt = optimset('jacobian','on');
end

switch numel(varargin)
    case 1
        [v1,v2,v3,v4,v5,v6,v7]=lsqnonlin(@symFun,x0,[],[],opt,F,Fx);
    case 2
        LB=varargin{2};
        [v1,v2,v3,v4,v5,v6,v7]=lsqnonlin(@symFun,x0,LB,[],opt,F,Fx);
    otherwise
        LB=varargin{2};
        UB=varargin{3};
        [v1,v2,v3,v4,v5,v6,v7]=lsqnonlin(@symFun,x0,LB,UB,opt,F,Fx);
end

function [E,J]=symFun(x0,F,Fx)
n=numel(x0);
for k=1:n
    eval(['x' num2str(k) '=x0(k);']);
end
E = eval(F);
if nargout>1
    J  = eval(Fx);
end


text:
%% Using Cook's Distance to Detect Outliers
% Copyright (c) 2011, The MathWorks, Inc.

% Create a vector of X values
clear all
clc
hold off

X = 1:100;
X = X';

% Create a noise vector
noise = randn(100,1);

% Create a second noise value where sigma is much larger
noise2 = 10*randn(100,1);

% Substitute noise2 for noise1 at obs# (11, 31, 51, 71, 91)
% Many of these points will have an undue influence on the model

noise(11:20:91) = noise2(11:20:91);

% Specify Y = F(X)
Y = 3*X + 2 + noise;

% Cook's Distance for a given data point measures the extent to
% which a regression model would change if this data point
% were excluded from the regression. Cook's Distance is
% sometimes used to suggest whether a given data point might be an outlier.

% Use regstats to calculate Cook's Distance
stats = regstats(Y,X,'linear');

% if Cook's Distance > n/4 is a typical treshold that is used to suggest
% the presence of an outlier
potential_outlier = stats.cookd > 4/length(X);

% Display the index of potential outliers and graph the results
X(potential_outlier)
scatter(X,Y, 'b.')
hold on
scatter(X(potential_outlier),Y(potential_outlier), 'r.')


text:
function v = basis_serene ( xq, yq, xw, ys, xe, yn, xx, yy )

%*****************************************************************************80
%
%% BASIS_SERENE evaluates the serendipity basis functions.
%
%  Discussion:
%
%    This procedure assumes that a serendipity element has been defined,
%    whose sides are parallel to coordinate axes.
%
%    The local element numbering is
%
%      YN  3--2--1
%       |  |     |
%       |  4     8
%       |  |     |
%      YS  5--6--7
%       |
%       +--XW---XE--
%
%    We note that each basis function can be written as the product of
%    three linear terms, which never result in an x^2y^2 term.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XQ, YQ, the evaluation point.
%
%    Input, real XW, YS, the coordinates of the lower left corner.
%
%    Input, real XE, YN, the coordinates of the upper right corner.
%
%    Input, real XX(8), YY(8), the coordinates of the 8 nodes.
%
%    Output, real V(8), the value of the basis functions at (XQ,YQ).
%
  v = zeros(8,1);

  v(1) = not1 ( xq, xw, xx(1) ) ...
       * not1 ( yq, ys, yy(1) ) ...
       * not2 ( xq, yq, xx(8), yy(8), xx(2), yy(2), xx(1), yy(1) );

  v(2) = not1 ( xq, xw, xx(2) ) ...
       * not1 ( xq, xe, xx(2) ) ...
       * not1 ( yq, ys, yy(2) );

  v(3) = not1 ( xq, xe, xx(3) ) ...
       * not1 ( yq, ys, yy(3) ) ...
       * not2 ( xq, yq, xx(2), yy(2), xx(4), yy(4), xx(3), yy(3) );

  v(4) = not1 ( xq, xe, xx(4) ) ...
       * not1 ( yq, yn, yy(4) ) ...
       * not1 ( yq, ys, yy(4) );

  v(5) = not1 ( xq, xe, xx(5) ) ...
       * not1 ( yq, yn, yy(5) ) ...
       * not2 ( xq, yq, xx(4), yy(4), xx(6), yy(6), xx(5), yy(5) );

  v(6) = not1 ( xq, xe, xx(6) ) ...
       * not1 ( xq, xw, xx(6) ) ...
       * not1 ( yq, yn, yy(6) );

  v(7) = not1 ( xq, xw, xx(7) ) ...
       * not1 ( yq, yn, yy(7) ) ...
       * not2 ( xq, yq, xx(6), yy(6), xx(8), yy(8), xx(7), yy(7) );

  v(8) = not1 ( yq, ys, yy(8) ) ...
       * not1 ( yq, yn, yy(8) ) ...
       * not1 ( xq, xw, xx(8) );

  return
end


text:
function LBP_Im = LBP(Input_Im, R)

% %=======================================================================
% %=======================================================================
% This function computes the LBP transformation of the input image
% Input_Im.
% The parameters of the LBP operator are (P = 8, R),
% where P - the number of sampling points in the region with the radius R.
% Radius R is the input parameter of the function.
% Possible values for R = 1, 2, 3, etc.
% If input image is COLOR, then the grayscale transformation is performed.
% %=======================================================================
% %=======================================================================

if size(Input_Im, 3) == 3
    Input_Im = rgb2gray(Input_Im);
end;
L = 2*R + 1; %% The size of the LBP label
C = round(L/2);
Input_Im = uint8(Input_Im);
row_max = size(Input_Im,1)-L+1;
col_max = size(Input_Im,2)-L+1;
LBP_Im = zeros(row_max, col_max);
for i = 1:row_max
    for j = 1:col_max
        A = Input_Im(i:i+L-1, j:j+L-1);
        A = A+1-A(C,C);
        A(A>0) = 1;
        LBP_Im(i,j) = A(C,L) + A(L,L)*2 + A(L,C)*4 + A(L,1)*8 + A(C,1)*16 + A(1,1)*32 + A(1,C)*64 + A(1,L)*128;
    end;
end;


text:
% Section 8.5.3: Analytic center of a set of linear inequalities
% Boyd & Vandenberghe "Convex Optimization"
% Joëlle Skaf - 04/29/08
%
% The analytic center of a set of linear inequalities and equalities:
%           a_i^Tx <= b_i   i=1,...,m,
%           Fx = g,
% is the solution of the unconstrained minimization problem
%           minimize    -sum_{i=1}^m log(b_i-a_i^Tx).

% Input data
randn('state', 0);
rand('state', 0);
n = 10;
m = 50;
p = 5;
tmp = randn(n,1);
A = randn(m,n);
b = A*tmp + 10*rand(m,1);
F = randn(p,n);
g = F*tmp;

% Analytic center
cvx_begin
    variable x(n)
    minimize -sum(log(b-A*x))
    F*x == g; %#ok
cvx_end

disp(['The analytic center of the set of linear inequalities and ' ...
      'equalities is: ']);
disp(x);

author: yu-jiang

repo: radpbook

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


text:
function sR = extractSphericalRegion(varargin)

% TODO: maybe consider option 'antipodal'
sR = getClass(varargin,'sphericalRegion',[],'last');
sR = sphericalRegion(sR,varargin{:});


text:
function batch_efficiency(dwcard)
% Start in directory above individual model/results directories
%
% :Usage:
% ::
%
%     function batch_efficiency(dwcard)
%
% :Inputs:
%
%   **dwcard:**
%        is a wildcard for directories to probe, e.g., 'subject*'
%
% ..
%    Tor Wager, 10/8/02
% ..

d = dir(dwcard);

mypwd = pwd;
eval(['diary ' mypwd filesep 'model_efficiency_stats.out'])

for i = 1:length(d)

	try
		spm_efficiency([d(i).name filesep 'SPM.mat']);
	catch
		disp(['Problem with subject ' d(i).name])
	end

end

diary off

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%*************************************************************************
%% iterrefine: Iterative refinement.
%% This step is crucial to ensure that computed solution
%% is sufficiently accuraete.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

   function [x,resnrm,solve_ok] = iterrefine(A,b,L,x0);

     tol = 1e-6;
     maxit = 10;
     bnorm = max(1,norm(b));
     tolb = tol*bnorm;
     x = x0;
     solve_ok = 1;
     resnrm(1) = bnorm;
%%
     for iter = 1:maxit
        x0 = x;
        if isstruct(A); r = b-matvec(A,x); else; r=b-A*x; end;
        err = norm(r);
        resnrm(iter+1) = err;

        if (err < tolb); break; end;
        if (iter > 1) & (resnrm(iter+1)/resnrm(iter) > 0.9)
           x = x0; solve_ok = 0; break;
        end
	d = linsysolvefun(L,r);
        x = x + d;
     end
%%*************************************************************************
%% matvec: matrix-vector multiply.
%% matrix = [A.mat11 A.mat12; A.mat12' A.mat22]
%%*************************************************************************

   function Ax = matvec(A,x);

   m = length(A.mat11); m2 = length(x)-m;
   x1 = x(1:m);
   Ax = A.mat11*x1;
   if (m2 > 0)
      x2 = x(m+[1:m2]);
      Ax = Ax + A.mat12*x2;
      Ax2 = (x1'*A.mat12)' + A.mat22*x2;
      Ax = [Ax; Ax2];
   end
   return;
%%*************************************************************************


text:
function [Q, Yh] = hadsst3_experiment_gpfa_run()

% GP priors for loadings, PCA/isotropic prior for states.

% Number of components
D = 80;

%
% Process data
%

data = hadsst3_load_data();

% Form the data matrix
Y = data.observations;
[M,N] = size(Y);
Obs = ~isnan(Y);

%
% GP model for spatial W
%

ind = 0;

covfunc_w = cell(D,1);
theta_w = cell(D,1);
is_pseudos_w = false(D,1);

% Pseudo inputs (uniformly with respect to area size)
pseudo_w = points_on_sphere(18); % uniform points by number of latitudes
% Remove pseudo inputs that are on land (the nearest grid point is land)
ind_pseudo_w = mohsst5_points_to_grid_index(pseudo_w);
pseudo_w(:,mohsst5_is_land_index(ind_pseudo_w)) = [];
% $$$ % This code shows the pseudo inputs on the map
% $$$ figure
% $$$ map_projection('global-ellipse');
% $$$ map_plot(pseudo_w,'r+');
% $$$ map_coast()
% $$$ map_grid()
% $$$ return

% Transform inputs to 3-D Euclidean coordinates
in_w = data.coordinates;
in_w = geographic_to_euclidean(in_w);
pseudo_w = geographic_to_euclidean(pseudo_w);

% Squared distance matrices for the covariance functions
D2_ww = sq_dist(in_w);
D2_pp = sq_dist(pseudo_w);
D2_pw = sq_dist(pseudo_w, in_w);
d2_w = diag(D2_ww);

ind = 1:D;
fprintf('%d slow components for W (using %d pseudo inputs)\n', length(ind), ...
        size(pseudo_w,2));

% Covariance function (scaled squared exponential) with pseudo inputs
covfunc = @(D2) gp_cov_se(D2);
covfunc_w(ind) = {gp_cov_pseudo(...
    gp_cov_scale(gp_cov_jitter(covfunc(D2_pp), 1e-3)), ...
    gp_cov_scale(covfunc(D2_pw)), ...
    gp_cov_scale(covfunc(d2_w)))};

% Hyperparameters for the covariance functions
theta_w(ind) = columns_to_cells(...
    [linspace(1,0.1,length(ind));       % magnitudes
     linspace(5000,1000,length(ind))]); % lengthscales
is_pseudos_w(ind) = true;

%% Short scale components: piecewise polynomial in 3-D

ind = (ind(end)+1):D;
fprintf('%d fast components for W\n', length(ind));

% Use block-Toeplitz structure for the covariance function
[lat,lon0] = meshgrid(data.latitude,...
                      data.longitude(1));
in_w0 = geographic_to_euclidean([lon0(:)';lat(:)']);
d_ww = sqrt(sq_dist(in_w0, ...
                    geographic_to_euclidean(data.coordinates)));
covfunc = gp_cov_toeplitz_block(gp_cov_pp(d_ww,3));
% Add scaling and jitter
covfunc_w(ind) = {gp_cov_scale(gp_cov_jitter(covfunc))};
% Hyperparameters
theta_w(ind) = columns_to_cells(...
    [linspace(0.5,0.1,length(ind));     % magnitude
     linspace(3000,2000,length(ind))]); % lengthscale
is_pseudos_w(ind) = false;

%
% GPFA inference
%

% PCA module for X
%X_module = factor_module_iid(D, N);
covfunc_x = cell(D,1);
theta_x = cell(D,1);
covfunc_x(:) = {gp_cov_delta(N)};
theta_x(:) = {[]};
X_module = factor_module_gp_factorized(N, covfunc_x, theta_x, ...
                                       'update_hyperparameters', [], ...
                                       'init', zeros(D,N));

% Component-wise factorization for W
W_module = factor_module_gp_factorized(M, covfunc_w, theta_w, ...
                                       'update_hyperparameters', [5 10:10:100 100:25:2000], ...
                                       'maxiter_hyperparameters', 5, ...
                                       'is_pseudo', is_pseudos_w);

% Isotropic noise (precisions weighted proportionally to grid size)
weights = mohsst5_weights();
weights = repmat(weights, [1, N]);
noise_module = noise_module_product(...
    noise_module_isotropic(M, N, ...
                           'prior', struct('a_tau', 1e-3, ...
                                           'b_tau', 1e-3), ...
                           'init', struct('a_tau', 1, ...
                                          'b_tau', 0.01)), ...
    noise_module_fixed(M, N, weights));

% Filename for saving the results
folder = sprintf('/share/bayes/jluttine/results/hadsst3/gpfa');
mkdir(folder);
filename = sprintf('%s/results_hadsst3_gpfa_D=%d_%s', ...
                   folder, ...
                   D, ...
                   datestr(now,30));

% Run GPFA
Q = vbfa(D, Y, W_module, X_module, noise_module, ...
         'maxiter', 200, ...
         'rotate', 1:50, ...
         'autosavefile', filename, ...
         'autosave', [1 5:10:2000]);

% Reconstruct
Yh = Q.W'*Q.X;

% Save the results
save(filename, '-struct', 'Q');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [data,units] = compute_dwing_angle_imbalance(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  imbalancer = trx(fly).wing_angler+trx(fly).wing_anglel;
  imbalancel = -imbalancer;
  dimbalancer = diff(imbalancer);
  dimbalancel = diff(imbalancel);
  data{i} = dimbalancel;
  idx = imbalancer(1:end-1) > imbalancel(1:end-1);
  data{i}(idx) = dimbalancer(idx);

  data{i} = data{i} ./ trx(fly).dt;

end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

X = [ones(m, 1) X];
t1 = sigmoid(X * Theta1');
t1 = [ones(m, 1) t1];

t2 = sigmoid( t1 * Theta2');

[~, p] = max(t2, [], 2);

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function [B,twom]=multicatdir_f(A,gamma,omega)
%MULTICATDIR_F  returns multilayer Leicht-Newman modularity matrix for categorical directed layers, function handle version
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
%   Input: A: Cell array of NxN adjacency matrices for each layer of a
%          categorical directed multilayer network
%          gamma: intralayer resolution parameter
%          omega: interlayer coupling strength
%
%   Output: B: function handle where B(i) returns the ith column of
%          [NxT]x[NxT] flattened modularity tensor for the
%           multilayer network with uniform ordinal coupling (T is
%           the number of layers of the network)
%           twom: normalisation constant
%
%   Example of usage: [B,twom]=multicatdir_f(A,gamma,omega);
%          [S,Q]= genlouvain(B); % see iterated_genlouvain.m and
%          postprocess_categorical_multilayer.m for how to improve output
%          multilayer partition
%          Q=Q/twom;
%          S=reshape(S,N,T);
%
%   [B,twom] = MULTICATDIR_F(A,GAMMA, OMEGA) with A a cell array of square
%   matrices of equal size each representing an directed network "layer"
%   computes the Leicht-Newman multilayer modularity matrix B using the
%   quality function described in Mucha et al. 2010, with intralayer
%   resolution parameter GAMMA, and with interlayer coupling OMEGA
%   connecting all-to-all categorical layers. Once the mulilayer modularity
%   matrix is computed, optimization can be performed by the generalized
%   Louvain code GENLOUVAIN or ITERATED_GENLOUVAIN. The output B can be used
%   with other heuristics, provided the same mapping is used to go from the
%   multilayer tensor to the multilayer flattened matrix. That is, the
%   node-layer tuple (i,s) is mapped to i + (s-1)*N. [Note that we can
%   define a mapping between a multilayer partition S_m stored as an N by T
%   matrix and the corresponding flattened partition S stored as an NT by 1
%   vector. In particular S_m = reshape(S,N,T) and S = S_m(:).]
%
%   See also
%       genlouvain heuristics:      GENLOUVAIN, ITERATED_GENLOUVAIN
%       multilayer wrappers:        MULTICAT, MULTICAT_F, MULTIORD
%       other heuristics:           SPECTRAL23
%       Kernighan-Lin improvement:  KLNB
%
%   Notes:
%     The matrices in the cell array A are assumed to be square,
%     and of equal size.  These assumptions are not checked here.
%
%     For smaller systems, it is potentially more efficient (and easier) to
%     directly use the sparse quality/modularity matrix B in MULTICAT. For
%     large systems with undirected layer networks, use MULTICAT_F.
%
%     This code serves as a template and can be modified for situations
%     with other wrinkles (e.g., different intralayer null models,
%     different numbers of nodes from layer-to-layer, or systems which are
%     both multiplex and longitudinal).  That is, this code is only a
%     starting point; it is by no means exhaustive.
%
%     By using this code, the user implicitly acknowledges that the authors
%     accept no liability associated with that use.  (What are you doing
%     with it anyway that might cause there to be a potential liability?!?)
%
%   References:
%     Blondel, Vincent D., Jean-Loup Guillaume, Renaud Lambiotte, and
%     Etienne Lefebvre, "Fast unfolding of communities in large networks,"
%     Journal of Statistical Mechanics: Theory and Experiment, P10008
%     (2008).
%
%     Fortunato, Santo, "Community detection in graphs," Physics Reports
%     486, 75-174 (2010).
%
%     Good, Benjamin H., Yves-Alexandre de Montjoye, and Aaron Clauset,
%     "Performance of modularity maximization in practical contexts,"
%     Physical Review E 81, 046106 (2010).
%
%     Mucha, Peter J., Thomas Richardson, Kevin Macon, Mason A. Porter, and
%     Jukka-Pekka Onnela. "Community Structure in Time-Dependent,
%     Multiscale, and Multiplex Networks," Science 328, 876-878 (2010).
%
%     Elizabeth A. Leicht and Mark E. J. Newman. "Community structure in
%     Directed Networks", Physical Review Letters 100, 118703 (2008).
%
%     Porter, M. A., J. P. Onnela, and P. J. Mucha, "Communities in
%     networks," Notices of the American Mathematical Society 56, 1082-1097
%     & 1164-1166 (2009).
%
%   Acknowledgments:
%     Thank you to Dani Bassett, Jesse Blocher, Bruce Rogers, and Simi Wang
%     for their collaborative help which led to significant cleaning up
%     of earlier versions of our multilayer community detection codes.

if nargin<2||isempty(gamma)
    gamma=1;
end

if nargin<3||isempty(omega)
    omega=1;
end

N=length(A{1});
T=length(A);

if length(gamma)==1
    gamma=repmat(gamma,T,1);
end

m=zeros(T,1);
for i=1:T
    m(i)=sum(A{i}(:));
end
A=blkdiag(A{:});
kout=sum(A,1);
koutmat=sparse(1:(N*T),kron(1:T,ones(1,N)),kout);
kin=sum(A,2);
kinmat=sparse(1:(N*T),kron(1:T,ones(1,N)),kin);
A=(A+A')./2;
all2all = N*[(-T+1):-1,1:(T-1)];
A = A + omega*spdiags(ones(N*T,2*T-2),all2all,N*T,N*T);

B=@(i) A(:,i)-gamma(ceil(i./(N+eps))).*(kout(i).*kinmat(:,ceil(i./(N+eps)))+kin(i).*koutmat(:,ceil(i./(N+eps))))./(2*m(ceil(i./(N+eps))));

twom=sum(m)+omega*2*N*(T-1);
end

author: GenLouvain

repo: GenLouvain