text:
function [x,u,y] = spm_Markov_blanket(J,z,m,mj)
% FORMAT [x,u,y] = spm_Markov_blanket(J,z,m,mj)
% Markovian partition
% J  - Jacobian
% z  - {1 x N}  partition of states (indices)
% m  - number of internal states [default: 3]
%
% x  - {3 x n} particular partition of state indices
%     x{1,j} - active states of j-th partition
%     x{2,j} - sensory states of j-th partition
%     x{3,j} - internal states of j-th partition
%
% u  - location of partitions in scaling or embedding space
%
% y  - {3 x n} particular partition of partition indices
%     y{1,j} - active states of j-th partition
%     y{2,j} - sensory states of j-th partition
%     y{3,j} - internal states of j-th partition
%
% Partition or Grouping (coarse-scaling) operator
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_Markov_blanket.m 7655 2019-08-25 20:10:20Z karl $

% preliminaries
%--------------------------------------------------------------------------
GRAPHICS  = 1;                      % Graphics switch
nz        = length(z);              % number of partitions
if nargin < 3
    m = 3;                          % maximum size of internal states
end
if nargin < 4
   mj = ones(nz,1);                 % eligible internal states
end
if isempty(mj)
   mj = ones(nz,1);                 % eligible internal states
end

% Adjacency matrix (over z)
%--------------------------------------------------------------------------
for i = 1:nz
    for j = 1:nz
        Lij    = J(z{i},z{j});
        if any(any(Lij))
            L(i,j) = abs(norm(full(Lij)) > 1/128);
        else
            L(i,j) = 0;
        end
    end
end
L     = double(L);

% get Markov blanket
%--------------------------------------------------------------------------
B     = L + L' + L'*L;
B     = B - diag(diag(B));

% scaling space (defined by graph Laplacian)
%--------------------------------------------------------------------------
G     = L + L';
G     = G - diag(diag(G));
G     = G - diag(sum(G));
G     = expm(G);

% get principal dimensions of scaling space (u)
%--------------------------------------------------------------------------
if GRAPHICS
    [u,v] = eig(G,'nobalance');
    v     = abs(diag(v));
    for i = 1:nz
        [p,h] = hist(real(u(:,i)),16);
        dh    = h(2) - h(1) + exp(-16);
        p     = p(:)/sum(p)/dh;
        v(i)  = log(v(i)) - p'*log(p + exp(-16))*dh;
    end
    [v,j] = sort(real(v),'descend');
    u     = real(u(:,j));
end

% recursive (particular) partition into internal, sensory and active states
%--------------------------------------------------------------------------
nn    = zeros(nz,1);
for i = 1:nz

    % internal states (defined by graph Laplacian)
    %----------------------------------------------------------------------
    jj = ~(B*nn) & ~nn & mj;
    if any(jj)

        % find densely coupled internal states (using the graph Laplacian)
        %------------------------------------------------------------------
        [g,j] = max(diag(G).*jj);
        if m > 1
            g      = G(:,j);
            g(j)   = 0;
            g(~jj) = 0;
            [g,k]  = sort(g,'descend');
            try
                j = [j; k(1:m - 1)];
            end
        end

        jj    = sparse(j,1,1,size(L,1),1) & jj;         % internal states
        bb    = B*jj & ~jj & ~nn;                       % Markov blanket
        ee    =  ~bb & ~jj & ~nn;                       % external states
        b     = find(bb);
        e     = find(ee);
        s     = b(find( any(L(b,e),2)));
        a     = b(find(~any(L(b,e),2)));

        % indices of individual states in the i-th particle
        %------------------------------------------------------------------
        x{1,i} = spm_cat(z(a));
        x{2,i} = spm_cat(z(s));
        x{3,i} = spm_cat(z(j));

        % states accounted for (nn)
        %------------------------------------------------------------------
        nn   = nn | bb | jj;

    else

        % no internal states - find active states (not influenced by e)
        %------------------------------------------------------------------
        j = ~any(L(~nn,nn),2);
        if any(j)

            % sensory states connected with active states
            %--------------------------------------------------------------
            a  = find(~nn);
            a  = a(find(j,1));
            aa = sparse(a,1,1,size(L,1),1);
            ss = (L*aa | L'*aa) & ~aa & ~nn;
            a  = find(aa);
            s  = find(ss);
            j  = [];

            % indices of individual states in the i-th particle
            %--------------------------------------------------------------
            x{1,i} = spm_cat(z(a));
            x{2,i} = spm_cat(z(s));
            x{3,i} = [];

            % states accounted for (nn)
            %--------------------------------------------------------------
            nn   = nn | aa | ss;

        elseif any(~nn)

            % sensory states connected with sensory states
            %--------------------------------------------------------------
            s  = find(~nn);
            ss = sparse(s(1),1,1,nz,1);
            ss = ss | B*ss & ~nn;
            s  = find(ss);
            a  = [];
            j  = [];

            % indices of individual states in the i-th particle
            %--------------------------------------------------------------
            x{1,i} = [];
            x{2,i} = spm_cat(z(s));
            x{3,i} = [];

            % states accounted for (nn)
            %--------------------------------------------------------------
            nn   = nn | ss;
        end
    end

    % indices of partitions (i.e., n-states) in the i-th particle
    %----------------------------------------------------------------------
    y{1,i} = a;
    y{2,i} = s;
    y{3,i} = j;

    % plot
    %----------------------------------------------------------------------
    if all(nn) && numel(u) > 1

        % remove isolated (internal) states
        %--------------------------------------------------------------
        j     = [];
        for n = 1:size(x,2)
            if any(x{1,n}) || any(x{2,n})
                j = [j,n];
            end
        end
        x  = x(:,j);
        y  = y(:,j);

        if GRAPHICS,clf

            % colours for different particles
            %--------------------------------------------------------------
            nx            = size(x,2);
            [col,bol,msz] = spm_MB_col(nx);

            % plot partitions in embedding space (which particle)
            %--------------------------------------------------------------
            subplot(3,2,3)
            for k = 1:nx
                plot(u(y{1,k},1),u(y{1,k},2),'.','color',bol{k},'MarkerSize',msz), hold on
                plot(u(y{2,k},1),u(y{2,k},2),'.','color',bol{k},'MarkerSize',msz), hold on
                plot(u(y{3,k},1),u(y{3,k},2),'.','color',col{k},'MarkerSize',msz), hold on
            end
            axis square
            title(sprintf('Particles [%i n-states]',nz),'Fontsize',16)

            % plot particles in embedding space (which sort of state)
            %--------------------------------------------------------------
            subplot(3,2,4)
            for k = 1:nx
                plot(u(y{1,k},1),u(y{1,k},2),'.r','MarkerSize',msz), hold on
                plot(u(y{2,k},1),u(y{2,k},2),'.m','MarkerSize',msz), hold on
                plot(u(y{3,k},1),u(y{3,k},2),'.b','MarkerSize',msz), hold on
            end
            axis square
            title(sprintf('Markov partition [%i particles]',nx),'Fontsize',16)

            % plot particles in three embedding dimensions
            %--------------------------------------------------------------
            subplot(3,2,2)
            try
                for k = 1:nx
                    plot3(u(y{1,k},1),u(y{1,k},2),u(y{1,k},3),'.r','MarkerSize',msz), hold on
                    plot3(u(y{2,k},1),u(y{2,k},2),u(y{2,k},3),'.m','MarkerSize',msz), hold on
                    plot3(u(y{3,k},1),u(y{3,k},2),u(y{3,k},3),'.b','MarkerSize',msz), hold on
                end
            catch
                for k = 1:nx
                    plot(u(y{1,k},1),u(y{1,k},2),'.r','MarkerSize',msz), hold on
                    plot(u(y{2,k},1),u(y{2,k},2),'.m','MarkerSize',msz), hold on
                    plot(u(y{3,k},1),u(y{3,k},2),'.b','MarkerSize',msz), hold on
                end
            end

            axis square
            title('Embedding space','Fontsize',16)
            rotate3d(gca,'on')

            % Jacobian (ordered by partition and type)
            %--------------------------------------------------------------
            j = spm_vec(x');
            k = spm_vec(x );
            subplot(3,2,5),imagesc(-log(abs(J(k,k)) + exp(-4))),axis square
            subplot(3,2,6),imagesc(-log(abs(J(j,j)) + exp(-4))),axis square

            % Colors
            %--------------------------------------------------------------
            nj   = spm_length(x);
            msz  = fix(16 + 128/nj);
            j    = 1:nj;
            k    = spm_unvec(j,x')';
            j    = spm_unvec(j,x);
            subplot(3,2,5),hold on
            for q = 1:nx
                plot(j{1,q},ones(size(x{1,q})),'.','color',bol{q},   'MarkerSize',msz)
                plot(j{2,q},ones(size(x{2,q})),'.','color',bol{q},   'MarkerSize',msz)
                plot(j{3,q},ones(size(x{3,q})),'.','color',col{q},   'MarkerSize',msz)
                plot(j{1,q},zeros(size(x{1,q})) + nj,'.','color','r','MarkerSize',msz)
                plot(j{2,q},zeros(size(x{2,q})) + nj,'.','color','m','MarkerSize',msz)
                plot(j{3,q},zeros(size(x{3,q})) + nj,'.','color','b','MarkerSize',msz)
            end
            title(sprintf('Jacobian (by %i particles)',nx),'Fontsize',16)

            subplot(3,2,6),hold on
            for q = 1:nx
                plot(k{1,q},ones(size(x{1,q})),'.','color',bol{q},   'MarkerSize',msz)
                plot(k{2,q},ones(size(x{2,q})),'.','color',bol{q},   'MarkerSize',msz)
                plot(k{3,q},ones(size(x{3,q})),'.','color',col{q},   'MarkerSize',msz)
                plot(k{1,q},zeros(size(x{1,q})) + nj,'.','color','r','MarkerSize',msz)
                plot(k{2,q},zeros(size(x{2,q})) + nj,'.','color','m','MarkerSize',msz)
                plot(k{3,q},zeros(size(x{3,q})) + nj,'.','color','b','MarkerSize',msz)
            end
            title('Jacobian (by type)','Fontsize',16)

        end
        break
    end

end

return

author: spm


text:
%% IMAGE_BATCH_LOCAL uses the BATCH command to run the IMAGE code locally.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'IMAGE_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run IMAGE_SCRIPT locally\n' );
%
%  Set the number of workers in the pool.
%
  n = 3;
%
%  BATCH sends the script for execution.
%
  job = batch ( 'image_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'image_fun', 'balloons.tif' }, ...
    'matlabpool', n );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( job );
%
%  LOAD makes the script's workspace avaiable.
%
%  x = output #1.
%  y = output #2.
%  z = output #3.
%
  load ( job );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( job );
%
%  Display the original, noisy, and filtered versions.
%
  figure;

  subplot ( 1, 3, 1 );
  imshow ( x );
  title ( 'Original image' );

  subplot ( 1, 3, 2 );
  imshow( y );
  title( 'Noisy Image' );

  subplot ( 1, 3, 3 );
  imshow ( z );
  title ( 'Median Filtered Image' );

author: johannesgerer

repo: jburkardt-m


text:
function [M,P,K,MU,S,LH] = ghkf_update(M,P,Y,h,R,h_param,p)
% GHKF_UPDATE - Gauss-Hermite Kalman filter update step
%
% Syntax:
%   [M,P,K,MU,S,LH] = GHKF_UPDATE(M,P,Y,h,R,param,p)
%
% In:
%   M  - Mean state estimate after prediction step
%   P  - State covariance after prediction step
%   Y  - Measurement vector.
%   h  - Measurement model function as a matrix H defining
%        linear function h(x) = H*x, inline function,
%        function handle or name of function in
%        form h(x,param)
%   R  - Measurement covariance
%   h_param - Parameters of h
%   p  - Degree of approximation (number of quadrature points)
%
% Out:
%   M  - Updated state mean
%   P  - Updated state covariance
%   K  - Computed Kalman gain
%   MU - Predictive mean of Y
%   S  - Predictive covariance Y
%   LH - Predictive probability (likelihood) of measurement.
%
% Description:
%   Perform additive form Gauss-Hermite Kalman filter (GHKF)
%   measurement update step. Assumes additive measurement
%   noise.
%
%   Function h(.) should be such that it can be given a
%   DxN matrix of N sigma Dx1 points and it returns
%   the corresponding measurements for each sigma
%   point. This function should also make sure that
%   the returned sigma points are compatible such that
%   there are no 2pi jumps in angles etc.
%
% Example:
%   h = inline('atan2(x(2,:)-s(2),x(1,:)-s(1))','x','s');
%   [M2,P2] = ghkf_update(M1,P1,Y,h,R,S);
%
% See also:
%   GHKF_PREDICT, GHRTS_SMOOTH, GH_TRANSFORM

% History:
%   Jun 18, 2009 - Initial version
%   May 24, 2010 - Fixed parameter input and description (asolin)
%   Aug 5,  2010 - Renamed from 'gh_update' to 'ghkf_update' (asolin)

% Copyright (C) 2009 Hartikainen, Särkkä, Solin
%
% $Id: gh_update.m,v 1.2 2009/07/01 06:34:41 ssarkka Exp $
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
%%

  %
  % Check that all arguments are there
  %
  if nargin < 5
     error('Too few arguments');
  end
  if nargin < 6
     h_param = [];
  end
  if nargin < 7
     p = [];
  end
  if isempty(p)
     p = 3;
  end

  %
  % Do the transform and make the update
  %
  tr_param = {p};
  [MU,S,C,X] = gh_transform(M,P,h,h_param,tr_param);
  S = S + R;
  K = C / S;
  M = M + K * (Y - MU);
  P = P - K * S * K';

  if nargout > 5
    LH = gauss_pdf(Y,MU,S);
  end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
classdef DTLZ9 < PROBLEM
% <multi/many> <real> <large/none> <constrained> <expensive/none>
% Benchmark MOP proposed by Deb, Thiele, Laumanns, and Zitzler

%------------------------------- Reference --------------------------------
% K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, Scalable test problems
% for evolutionary multiobjective optimization, Evolutionary multiobjective
% Optimization. Theoretical Advances and Applications, 2005, 105-145.
%------------------------------- 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)
            if isempty(obj.M); obj.M = 2; end
            if isempty(obj.D); obj.D = 10*obj.M; end
            obj.D        = ceil(obj.D/obj.M)*obj.M;
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            PopDec = varargin{1};
            PopDec = max(min(PopDec,repmat(obj.upper,size(PopDec,1),1)),repmat(obj.lower,size(PopDec,1),1));
            X      = PopDec;
            PopDec = PopDec.^0.1;
            PopObj = zeros(size(PopDec,1),obj.M);
            for m = 1 : obj.M
                PopObj(:,m) = sum(PopDec(:,(m-1)*obj.D/obj.M+1:m*obj.D/obj.M),2);
            end
            PopCon = 1 - repmat(PopObj(:,obj.M).^2,1,obj.M-1) - PopObj(:,1:obj.M-1).^2;
            Population = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            Temp = (0:1/(N-1):1)';
            R    = [repmat(cos(0.5.*pi.*Temp),1,obj.M-1),sin(0.5.*pi.*Temp)];
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M < 4
                R = obj.GetOptimum(100);
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverse kinematics for the 3dof planar robot
% T: homogeneous matrix
% robot: structure with arm parameters
% returns: all possible solutions or q = [q1 q2] that place the end effectors at the
% position specified by T. Two possible solutions q1 and q2 q3 and q4 are returned,
% generally called elbow up and elbow down, combined with
%   Author: Arturo Gil Aparicio arturo.gil@umh.es
%   Date: 08/03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_3dofplanar(robot, T)

fprintf('\nComputing inverse kinematics for the %s robot', robot.name);

%Initialize q
% q = [q1 q2 q3 q4], atends for two possible solutions
q=zeros(3,4);

a = eval(robot.DH.a);

%Link lengths
L1=abs(a(1));
L2=abs(a(2));
L3=abs(a(2));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Q=T(1:3,4);

%find angle Phi
x3 = T(1:3,1);

cphi = x3'*[1 0 0]';
sphi = x3'*[0 1 0]';
phi = atan2(sphi, cphi);

%Find point P from P and vector (nx, ny, nz)
P = Q - a(3)*T(1:3,1);

%Distance of the point to the origin.
R= sqrt(P(1)^2+P(2)^2);

if R > (L1+L2)
   disp('\ninversekinematic_3dofplanar: unfeasible solution. The point cannot be reached');
end

%compute geometric solution
beta = atan2(P(2),P(1));
gamma = real(acos((L1^2+R^2-L2^2)/(2*R*L1)));
delta = real(acos((L1^2+L2^2-R^2)/(2*L1*L2)));

%arrange possible combinations for q(1) and q(2)
%elbow down     elbow up solutions
q =[beta+gamma beta-gamma;
    delta-pi   pi-delta];

%in this case, phi = q(1) + q(2) + q(3) and
%q(3) can be computed as q(3) = phi - q(1) - q(2)
%corresponding to each of the previous solutions, a unique q(3) can be
%computed for each case
for i=1:2 %iterate through columns
    q(3,i) = phi - q(1,i) - q(2,i);
end

repo: ARTE


text:
function [M] = covar(x1,x2,k,n)
%
% this function calculates the lag k covariance matrice
% for three time series. Each time series is a 1xn vector
%
if k == 0
   %
   % lag 0 covariance matrix
   %
   p1=find(x1~=0);
   p2=find(x2~=0);
   %p3=find(x3~=0);
   p12=find(x1~=0 & x2~=0);
   %p13=find(x1~=0 & x3~=0);
   %p23=find(x2~=0 & x3~=0);
 %  m1=mean(x1(p1));
 %  m2=mean(x2(p2));
 %  m3=mean(x3(p3));
   m1=0;    % x1, x2 and x3 are residuals characterized by N(0,1)
   m2=0;
   %m3=0;
 %  M(1,1)=std(x1(p1));
 %  M(2,2)=std(x2(p2));
 %  M(3,3)=std(x3(p3));
   M(1,1)=1;   % x1, x2 and x3 are residuals characterized by N(0,1)
   M(2,2)=1;
   %M(3,3)=1;
   M(1,2)=sum((x1(p12)-m1).*(x2(p12)-m2))/(length(p12)-1);
   M(2,1)=M(1,2);
   %M(1,3)=sum((x1(p13)-m1).*(x3(p13)-m3))/(length(p13)-1);
   %M(3,1)=M(1,3);
   %M(2,3)=sum((x2(p23)-m2).*(x3(p23)-m3))/(length(p23)-1);
   %M(3,2)=M(2,3);
else
   %
   % lag k covariance matrix
   %
   x1t=x1(1:n-k);
   x1tk=x1(k+1:n);
   x2t=x2(1:n-k);
   x2tk=x2(k+1:n);
   %x3t=x3(1:n-k);
   %x3tk=x3(k+1:n);
   p1t=find(x1t~=0);
   p2t=find(x2t~=0);
   %p3t=find(x3t~=0);
   p1tk=find(x1tk~=0);
   p2tk=find(x2tk~=0);
   %p3tk=find(x3tk~=0);
 %  m1t=mean(x1t(p1t));
 %  m1tk=mean(x1tk(p1tk));
 %  m2t=mean(x2t(p2t));
 %  m2tk=mean(x2tk(p2tk));
 %  m3t=mean(x3t(p3t));
 %  m3tk=mean(x3tk(p3tk));
 m1t=0;   % x1, x2 and x3 are residuals characterized by N(0,1)
 m1tk=0;
 m2t=0;
 m2tk=0;
 %m3t=0;
 %m3tk=0;
   %
   p11=find(x1t~=0 & x1tk~=0);
   p22=find(x2t~=0 & x2tk~=0);
   %p33=find(x3t~=0 & x3tk~=0);
   p12=find(x1t~=0 & x2tk~=0);
   %p13=find(x1t~=0 & x3tk~=0);
   p21=find(x2t~=0 & x1tk~=0);
   %p23=find(x2t~=0 & x3tk~=0);
   %p31=find(x3t~=0 & x1tk~=0);
   %p32=find(x3t~=0 & x2tk~=0);
   %
   M(1,1)=sum((x1t(p11)-m1t).*(x1tk(p11)-m1tk))/(length(p11)-k-1);
   M(2,2)=sum((x2t(p22)-m2t).*(x2tk(p22)-m2tk))/(length(p22)-k-1);
   %M(3,3)=sum((x3t(p33)-m3t).*(x3tk(p33)-m3tk))/(length(p33)-k-1);
   M(1,2)=sum((x1t(p12)-m1t).*(x2tk(p12)-m2tk))/(length(p12)-k-1);
   M(2,1)=sum((x2t(p21)-m2t).*(x1tk(p21)-m1tk))/(length(p21)-k-1);
   %M(1,3)=sum((x1t(p13)-m1t).*(x3tk(p13)-m3tk))/(length(p13)-k-1);
   %M(3,1)=sum((x3t(p31)-m3t).*(x1tk(p31)-m1tk))/(length(p31)-k-1);
   %M(2,3)=sum((x2t(p23)-m2t).*(x3tk(p23)-m3tk))/(length(p23)-k-1);
   %M(3,2)=sum((x3t(p32)-m3t).*(x2tk(p32)-m2tk))/(length(p32)-k-1);
end

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 the National Renewable Energy Laboratory and Sandia Corporation
%
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plots solver-to-solver comparison for the newly run simulation for the
% regular waves.

%% Load Data
cd regular
load('regular_org.mat')       % Load WEC-Sim Run Data
cd .. ; cd regularCIC
load('regularCIC_org.mat')    % Load WEC-Sim Run Data
cd .. ; cd regularSS
load('regularSS_org.mat')     % Load WEC-Sim Run Data
cd ..

%% Plot Heave Comparisons
h=figure('units','normalized','outerposition',[0 0 1 1]);
% First Row = All simulations, all times
% Second Row = All simulations, t=350-400s
% Third Row = Original and New WEC-Sim runs only, t=350-400s

%First Column: Body 1 Heave
startTime = 100;
subplot(2,3,1)
m = plot(regular.B1.WEC_Sim_new.time,regular.B1.WEC_Sim_new.heave,'k:',...
	regularCIC.B1.WEC_Sim_new.time,regularCIC.B1.WEC_Sim_new.heave,'r--',...
    regularSS.B1.WEC_Sim_new.time,regularSS.B1.WEC_Sim_new.heave,'b-.');
a=get(m(1),'LineWidth')+1;
set(m(1),'LineWidth',a)
xlabel('time(s)')
ylabel('Heave(m)')
title('Float 1DOF 1200PTO')
xlim([0 150])
ylim([-2 2])

subplot(2,3,4)
m = plot(   regular.B1.WEC_Sim_new.time(find(    regular.B1.WEC_Sim_new.time==startTime):end),...
          regular.B1.WEC_Sim_new.heave(find(   regular.B1.WEC_Sim_new.time==startTime):end),'k:',...
       regularCIC.B1.WEC_Sim_new.time(find( regularCIC.B1.WEC_Sim_new.time==startTime):end),...
       regularCIC.B1.WEC_Sim_new.heave(find(regularCIC.B1.WEC_Sim_new.time==startTime):end),'r--',...
        regularSS.B1.WEC_Sim_new.time(find(  regularSS.B1.WEC_Sim_new.time==startTime):end),...
        regularSS.B1.WEC_Sim_new.heave(find( regularSS.B1.WEC_Sim_new.time==startTime):end),'b-.');
a=get(m(1),'LineWidth')+1;
set(m(1),'LineWidth',a)
xlabel('time(s)')
ylabel('Heave(m)')
xlim([100 150])
ylim([-2 2])

%Second Column: Body 2 Heave
subplot(2,3,2)
m = plot(   regular.B2.WEC_Sim_new.time,   regular.B2.WEC_Sim_new.heave,'k:',...
     regularCIC.B2.WEC_Sim_new.time,regularCIC.B2.WEC_Sim_new.heave,'r--',...
      regularSS.B2.WEC_Sim_new.time, regularSS.B2.WEC_Sim_new.heave,'b-.');
a=get(m(1),'LineWidth')+1;
set(m(1),'LineWidth',a)
xlabel('time(s)')
ylabel('Heave(m)')
title('Spar/Plate 1DOF 1200PTO')
xlim([0 150])
ylim([-0.2 0.2])

subplot(2,3,5)
m=plot(   regular.B2.WEC_Sim_new.time(find(    regular.B2.WEC_Sim_new.time==startTime):end),...
          regular.B2.WEC_Sim_new.heave(find(   regular.B2.WEC_Sim_new.time==startTime):end),'k:',...
       regularCIC.B2.WEC_Sim_new.time(find( regularCIC.B2.WEC_Sim_new.time==startTime):end),...
       regularCIC.B2.WEC_Sim_new.heave(find(regularCIC.B2.WEC_Sim_new.time==startTime):end),'r--',...
        regularSS.B2.WEC_Sim_new.time(find(  regularSS.B2.WEC_Sim_new.time==startTime):end),...
        regularSS.B2.WEC_Sim_new.heave(find( regularSS.B2.WEC_Sim_new.time==startTime):end),'b-.');
a=get(m(1),'LineWidth')+1;
set(m(1),'LineWidth',a)
xlabel('time(s)')
ylabel('Heave(m)')
xlim([100 150])
ylim([-0.2 0.2])

%Third Column: Relative Heave
subplot(2,3,3)
m= plot(   regular.Rel.WEC_Sim_new.time,   regular.Rel.WEC_Sim_new.heave,'k:',...
     regularCIC.Rel.WEC_Sim_new.time,regularCIC.Rel.WEC_Sim_new.heave,'r--',...
      regularSS.Rel.WEC_Sim_new.time, regularSS.Rel.WEC_Sim_new.heave,'b-.');
a=get(m(1),'LineWidth')+1;
set(m(1),'LineWidth',a)
xlabel('time(s)')
ylabel('Heave(m)')
title('Relative Motion 1DOF 1200PTO')
xlim([0 150])
ylim([-2 2])
l=legend('"regular"','"regularCIC"','"regularCIC" with State Space');
set(l,'Position',[0.765 0.84 0.07 0.07],'Units','normalized',...
    'FontSize',12);

subplot(2,3,6)
m=plot(   regular.Rel.WEC_Sim_new.time(find(    regular.Rel.WEC_Sim_new.time==startTime):end),...
          regular.Rel.WEC_Sim_new.heave(find(   regular.Rel.WEC_Sim_new.time==startTime):end),'k:',...
       regularCIC.Rel.WEC_Sim_new.time(find( regularCIC.Rel.WEC_Sim_new.time==startTime):end),...
       regularCIC.Rel.WEC_Sim_new.heave(find(regularCIC.Rel.WEC_Sim_new.time==startTime):end),'r--',...
        regularSS.Rel.WEC_Sim_new.time(find(  regularSS.Rel.WEC_Sim_new.time==startTime):end),...
        regularSS.Rel.WEC_Sim_new.heave(find( regularSS.Rel.WEC_Sim_new.time==startTime):end),'b-.');
a=get(m(1),'LineWidth')+1;
set(m(1),'LineWidth',a)
xlabel('time(s)')
ylabel('Heave(m)')
xlim([100 150])
ylim([-2 2])

clear a h l m B1_H_Tolerance B2_H_Tolerance Rel_H_Tolerance str1 str2

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
classdef WeibullD
%%WEIBULLD Functions to handle the Weibull distribution.
%Implemented methods are: mean, var, PDF, CDF, invCDF, rand, entropy
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(lambda,k)
%%MEAN Obtain the mean of the Weibull distribution for given scale and
%      shape parameters.
%
%INPUTS: lambda The scale parameter of the distribution. lambda>0.
%             k The shape parameter of the distribution. k>0.
%
%OUTPUTS:val The mean of the Weibull distribution.
%
%The mean of the Weibull distribution is given on the inside of the front
%cover of [1].
%
%EXAMPLE:
%We verify the computed mean by comparing it to the sample mean.
% lambda=2.2;
% k=3.1;
% meanVal=WeibullD.mean(lambda,k)
% numSamp=1e6;
% sampMeanVal=mean(WeibullD.rand([numSamp,1],lambda,k))
%One will find both values are about 1.967
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

   val=lambda*gamma(1+1/k);
end

function val=var(lambda,k)
%%VAR Obtain the variance of the Weibull distribution for given scale and
%     shape parameters/
%
%INPUTS: lambda The scale parameter of the distribution. lambda>0.
%             k The shape parameter of the distribution. k>0.
%
%OUTPUTS:val The variance of the Weibull distribution.
%
%The variance of the Weibull distribution is given on the inside of the front
%cover of [1].
%
%EXAMPLE:
%We verify the computed variance by comparing it to the sample variance.
% lambda=2.2;
% k=3.1;
% varVal=WeibullD.var(lambda,k)
% numSamp=1e6;
% sampVarVal=var(WeibullD.rand([numSamp,1],lambda,k))
%One will find both values are about 0.4821
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

   val=lambda^2*(gamma(1+2/k)-(gamma(1+1/k))^2);
end

function val=PDF(x,lambda,k)
%%PDF Evaluate the probability density function (PDF) of the Weibull
%     distribution at one or more points.
%
%INPUTS: x The point(s) at which the Weibull PDF is to be evaluated. x>=0
%          for nonzero PDF values.
%   lambda The scale parameter of the distribution. lambda>0.
%        k The shape parameter of the distribution. k>0.
%
%OUTPUTS: val The value(s) of the Weibull PDF.
%
%The PDF of the Weibull distribution is given on the inside of the front
%cover of [1].
%
%EXAMPLE:
%Here, we validate the PDF by generating random samples and comparing the
%PDF plot with a histogram of the random samples.
% lambda=2.2;
% k=3.1;
% numSamples=1e5;
% figure(1)
% clf
% histogram(WeibullD.rand([numSamples,1],lambda,k),'Normalization','pdf')
% hold on
% numPoints=1000;
% x=linspace(0,6,numPoints);
% vals=WeibullD.PDF(x,lambda,k);
% plot(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=(k/lambda)*(x/lambda).^(k-1).*exp(-(x/lambda).^k);
    val(x<0)=0;
end

function val=CDF(x,lambda,k)
%%CDF Evaluate the cumulative distribution function (CDF) of the Weibull
%     distribution at one or more desired points.
%
%INPUTS: x The point(s) at which the Weibull CDF is to be evaluated.
%   lambda The scale parameter of the distribution. lambda>0.
%        k The shape parameter of the distribution. k>0.
%
%OUTPUTS: val The value(s) of the Weibull CDF.
%
%The CDF is the integral from 0 to x of the PDF. The Weibull PDF is
%fairely simple and the integral is not hard. The Weibull PDF is given on
%the inside of the front cover of [1].
%
%EXAMPLE:
%We validate the CDF value by comparing it to a value computed from random
%samples.
% lambda=2.2;
% k=3.1;
% x=2;
% numSamples=1e5;
% prob=WeibullD.CDF(x,lambda,k)
% probSamp=mean(WeibullD.rand([numSamples,1],lambda,k)<=x)
%One will see that both values are about 0.5249.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=1-exp(-(x/lambda).^k);
    val(x<0)=0;
end

function val=invCDF(prob,lambda,k)
%%INVCDF Evaluate the inverse of the cumulative distribution function (CDF)
%        of the Weibull distribution.
%
%INPUTS: prob The probability or probabilities (0<=prob<=1) at which the
%             argument of the CDF is desired.
%      lambda The scale parameter of the distribution. lambda>0.
%           k The shape parameter of the distribution. k>0.
%
%OUTPUTS: val The argument(s) of the CDF that would give the probability or
%             probabilities in prob.
%
%The CDF of the Weibull distribution can be easily algebraicly inverted,
%which is what is done here. The Weibull CDF is the integral from 0 to x of
%the Weibull PDF, which is given on the inside of the front cover of [1].
%
%EXAMPLE:
%Here, we validate the inverse CDF by showing it to be the inverse of the
%CDF.
% lambda=2.2;
% k=3.1;
% x=2;
% xBack=WeibullD.invCDF(WeibullD.CDF(x,lambda,k),lambda,k)
%One will see that xBack is the same as x.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

val=lambda*(-log(1-prob)).^(1/k);

end

function vals=rand(N,lambda,k,lowerBound)
%%RAND Generate Weibull distributed random variables with the given
%      parameters. This function can also be used to generate random
%      variable from a Weibull distribution that is clipped so that only
%      values above a certain threshold are generated. That can be useful
%      when performing simulations and generating false alarms that are all
%      above a given threshold.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%   lambda The scale parameter of the distribution. lambda>0.
%        k The shape parameter of the distribution. k>0.
% lowerBound An optional parameter specifying a lower bound below which
%          none of the randomly generated values should go. If this is
%          omitted or an empty matrix is passed, then lowerBound=0 and a
%          non-clipped Weibull distribution is used.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated Weibull random variables.
%
%This is an implementation of the inverse transform algorithm of Chapter
%5.1 of [1]. To clip values to only being above lowerBound, the CDF value
%of the lower bound is determined and the uniform random variables as used
%in the inverse transform algorithm are all only generated above the
%threshold.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, Ed. 4, Amsterdam: Elsevier, 2006.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N, N];
    else
        dims=N;
    end

    if(nargin<4||isempty(lowerBound)||lowerBound<=0)
        %If we are evaluating a standard Weibull-distributed random
        %variable rather than one that has been clipped to a particular
        %region.
        U=rand(dims);
        vals=WeibullD.invCDF(U,lambda,k);
    else
        %Only generate Weibull random variables that are >=lowerBound.
        p=WeibullD.CDF(lowerBound,lambda,k);
        %Generate a uniform random variable between p and 1.
        U=p+(1-p)*rand(dims);
        vals=WeibullD.invCDF(U,lambda,k);
    end
end

function entropyVal=entropy(lambda,k)
%%ENTROPY Obtain the differential entropy of the Weibull distribution
%         given in nats. The differential entropy of a continuous
%         distribution is entropy=-int_x p(x)*log(p(x)) dx where the
%         integral is over all values of x. Units of nats mean that the
%         natural logarithm is used in the definition. Unlike the Shannon
%         entropy for discrete variables, the differential entropy of
%         continuous variables can be both positive and negative.
%
%INPUTS: lambda The scale parameter of the distribution. lambda>0.
%        k The shape parameter of the distribution. k>0.
%
%OUTPUTS: entropyVal The value of the differential entropy in nats.
%
%Differential entropy is defined in Chapter 8 of [1].
%
%REFERENCES:
%[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    gammaVal=-psi(1);
    entropyVal=gammaVal*(1-1/k)+log(lambda/k)+1;
end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function determ = lotkin_determinant ( n )

%*****************************************************************************80
%
%% LOTKIN_DETERMINANT returns the determinant of the LOTKIN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  delta = 1.0;

  for i = 2 : n
    delta = - r8_choose ( 2 * i - 2, i - 2 ) * r8_choose ( 2 * i - 2, i - 1 ) ...
      * ( 2 * i - 1 ) * delta;
  end

  determ = 1.0 / delta;

  return
end


text:
% display a matrix in a format of %0.mf
function outStrMat = mat2strf(inNumMat,dispFmt)
[nrow,ncol] = size(inNumMat);
outStrMat = '[';
for i=1:nrow
    for j=1:ncol
        if j~=ncol
            outStrMat = [outStrMat,sprintf(dispFmt,inNumMat(i,j)),','];
        else
            if nrow~=1
                outStrMat = [outStrMat,sprintf(dispFmt,inNumMat(i,j)),';'];
            else
                outStrMat = [outStrMat,sprintf(dispFmt,inNumMat(i,j))];
            end
        end
    end
end
outStrMat = [outStrMat,']'];
end

repo: paper-simulation


text:
classdef ProjectorToP1discont < handle
    % OBJECTIU: encara no fa falta resoldre el sistema matricial. Simplement
    % fer switch depenent de la funcio d'entrada (P0 o P1) i trobar "a lo
    % cutre" la P1discont. Ja mes endavant fem el cas general. Ja hi ha
    % fragments de codi escampats que resolen P0toP1discont i P1toP1discont.

    properties (Access = public)

    end

    properties (Access = private)
        mesh
%         connec
%         quadOrder
    end

    properties (Access = private)
        quadrature
%         field
%         meshD
    end

    methods (Access = public)

        function obj = ProjectorToP1discont(cParams)
            obj.init(cParams);
%             obj.createDiscontinuousMesh();
            obj.createQuadrature();
%             obj.createField();
        end

%         function xFun = project(obj, x)
%             LHS = obj.computeLHS();
%             RHS = obj.computeRHS(x);
%         end

        function xFun = project(obj,cParams) % THIS WILL BE DELETED ONCE RHS IS OK
            origin = cParams.origin;
            x      = cParams.x;
            connec = obj.mesh.connec;
            switch origin
                case {'P0'}
                    dim      = 1;
                    ndim  = size(x.fValues, 1);
                    nnodeElem = size(connec,2);
                    fEl = squeeze(x.fValues(dim,:,:));
                    fRepeated = zeros(ndim, size(fEl,1), nnodeElem);
                    for idim = 1:ndim
                        fEl = squeeze(x.fValues(idim,:,:));
                        for iNode = 1:nnodeElem
                            fRepeated(idim, :,iNode) = fEl;
                        end
                    end
                    fD = permute(fRepeated, [1 3 2]);
                    fD = fD(1,:,:);
                    s.fValues = fD;
                    s.mesh    = obj.mesh;
                    xFun      = P1DiscontinuousFunction(s);
                case {'P1'}
                    f = x.fValues;
                    nNode  = size(connec,2);
                    nDime  = size(f,2);
                    nodes = reshape(connec',1,[]);
                    fe = f(nodes,:)';
                    fNodeElem = reshape(fe,nDime,nNode,[]);

                    s.fValues = fNodeElem;
                    s.mesh    = obj.mesh;
                    xFun = P1DiscontinuousFunction(s);
            end
        end

    end

    methods (Access = private)

        function init(obj, cParams)
            obj.mesh   = cParams.mesh;
        end

%         function createDiscontinuousMesh(obj)
%             obj.meshD = obj.mesh.createDiscontinuousMesh();
%         end

        function q = createQuadrature(obj)
            quadOrder = 'LINEAR';
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature(quadOrder);
            obj.quadrature = q;
        end

%         function createField(obj)
%             s.mesh               = obj.meshD;
%             s.ndimf              = 1;
%             s.interpolationOrder = 'LINEAR';
%             s.quadratureOrder    = 'QUADRATIC';
%             obj.field = Field(s);
%         end
%
%         function LHS = computeLHS(obj)
%             s.type  = 'MassMatrix';
%             s.mesh  = obj.meshD;
%             s.field = obj.field;
%             lhs = LHSintegrator.create(s);
%             LHS = lhs.compute();
%         end
%
%         function RHS = computeRHS(obj, fun)
%             fDisc = fun.computeDiscontinuousField();
%             fVals = fDisc.fValues;
%             dV = obj.mesh.computeDvolume(obj.quadrature);
%             xV = obj.quadrature.posgp;
%             nGaus  = obj.quadrature.ngaus;
%             nF     = size(fVals,1);
%             nElem  = size(obj.mesh.connec,1);
%             rhs = zeros(nElem,nF);
%
%             % Separate in two loops
%             for igaus = 1:nGaus
%                 xGaus = xV(:,igaus);
%                 fGaus = fDisc.evaluate(xGaus);
%                 dVg(:,1) = dV(igaus,:);
%                 for iF = 1:nF
%                     fGausF = squeeze(fGaus(iF,:,:));
%                     Ni = 1;
%                     int = Ni*fGausF.*dVg;
%                     rhs(:,iF) = rhs(:,iF) + int;
%                 end
%             end
%         end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function F = spm_ncFcdf(x,df,d)
% Cumulative Distribution Function (CDF) of non-central F-distribution
% FORMAT f = spm_ncFcdf(x,df,d)
% x  - F-variate (F has range [0,Inf) )
% df - degrees of freedom, df = [v,w] with v>0 and w>0
% d  - non-centrality parameter
% F  - CDF of non-central F-distribution with [v,w] d.f. at points x
%
% Reference:
% https://en.wikipedia.org/wiki/Noncentral_F-distribution
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_ncFcdf.m 7354 2018-06-22 10:44:22Z guillaume $

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
%-Unpack degrees of freedom v & w from single df parameter (df)
if numel(df) == 2
    v = df(1);
    w = df(2);
elseif size(df,2) == 2
    v = df(:,1);
    w = df(:,2);
else
    error('Cannot unpack degrees of freedom.');
end

%-Check argument sizes
ad = [ndims(x);ndims(v);ndims(w);ndims(d)];
rd = max(ad);
as = [[size(x),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))];...
      [size(w),ones(1,rd-ad(3))];...
      [size(d),ones(1,rd-ad(4))];];
rs = max(as);
xa = prod(as,2)>1;
if all(xa) && any(any(diff(as(xa,:)),1))
    error('Non-scalar arguments must match in size.');
end

%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
F = zeros(rs);

%-Only defined for v>0, w>0 and d>=0. Return NaN if undefined.
md = true(size(F)) & v>0 & w>0 & d>=0;
if any(~md(:))
    F(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Compute where defined
if ~any(md(:)), return, end
Q = md;
if xa(1), Qx=Q; else Qx=1; end
if xa(2), Qv=Q; else Qv=1; end
if xa(3), Qw=Q; else Qw=1; end
if xa(4), Qd=Q; else Qd=1; end

a = exp(-d(Qd)/2);
x = v(Qv) .* x(Qx) ./ (w(Qw) + v(Qv) .* x(Qx));
F(Q) = a .* betainc(x(Qx), v(Qv)/2, w(Qw)/2);
for i=1:1024
    a = a .* d(Qd) / (2 * i);
    % could use recursion with:
    % Ix(a+1,b)=Ix(a,b)-G(a+b)/(G(a+1)*G(b))*x^a*(1-x)^b and G(a+1)=a*G(a)
    e = a .* betainc(x(Qx), v(Qv)/2+i, w(Qw)/2);
    if max(e) < 1e-12, break; end
    F = F + e;
end


text:
function out=testDependencyStructure()
out=0;

which('initCobraToolbox.m')
global CBTDIR

%excluded folders and subfolders
exclusionFolderList={'external'};

exclusion_list

function Build_Dependency_Map(root_dir, exclusion_list);

author: opencobra

repo: cobratoolbox


text:
% example of usage of package bruteForceBinaryPairwiseMex
%
% Anton Osokin,  03.04.2015

numNodes = 4;
numEdges = 5;

% [Dp(1), Dp(2)] - unary terms
unaryPotentials=[
    0,16;
    0,13;
    20,0;
    4,0
];

% [p, q, Vpq(1, 1), Vpq(1, 2), Vpq(2, 1), Vpq(2, 2)] - pairwise terms
pairwisePotentials=[
    1,2,0,10,4,0;
    1,3,0,12,-1,0;
    3,2,0,9,-1,0;
    2,4,0,14,0,0;
    3,4,0,0,7,0
    ];

[energy, labels] = bruteForceBinaryPairwiseMex(unaryPotentials, pairwisePotentials);

% % correct answer:
% energy = 22;
% labels = [1; 1; 2; 1];

if ~isequal(energy, 22)
    warning('Wrong value of the energy!')
end
if ~isequal(labels, [1; 1; 2; 1])
    warning('Wrong values of the labels!')
end

author: aosokin


text:
function H=prh(H)
%function H=prh(H)
%
% sets the current figure to print on full page in LANDSCAPE mode.
%
% By Michael Wu  --  waftingpetal@yahoo.com (Oct 2001)
%
% ====================

if ~exist('H','var');
  H = gcf;
end;

papMarg=0.1;
set(H,'PaperOrientation','landscape','PaperPosition',[0+papMarg, 0+papMarg, 11-2*papMarg, 8.5-2*papMarg]);


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured_SYM';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volume'};
optimizer = 'MMA';
incrementFactor = 1;
filterType = 'P1';
designVariable = 'Density';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';
showBC = true;

nsteps = 1;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
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;


text:
function [y,dzdw] = vl_nnscalenorm(x, w, varargin)
%VL_NNSCALENORM Feature normalization with scaling.
%   Y = VL_NNSCALENORM(X, W) normalizes input features across
%   across channels to have an L2 norm of 1, and scales each
%   channel by a learnable weight according to the following
%   formula:
%
%      Y(i,j,c,n) = X(i,j,c,n) * W(1,1,c) / sum_{k=1}^C (X(i,j,k,n)^2)
%
%   where X is an H x W x C x N 4D tensor input, Y is the
%   H x W x C x N tensor output and W is a 1 x 1 x Q tensor of
%   weights whose number of channels `Q` depends on the choice of
%   `channelShared` option value described below.
%
%   [DX, DW] = VL_NNSCALENORM(X, W, DY) computes the derivatives of
%   the operator projected onto P. DX and DW have the same
%   dimensions as X and W. The derivatives with respect to the weights
%   W are given by
%
%      dzDW(1,1,c,1) = X(i,j,k,n) / sum_{k=1}^K (X(i,j,k,n)^2)
%
%   NOTES: This layer was introduced in the paper: "ParseNet: Looking Wider
%   To See Better". It is useful when combining activations from different
%   layers of the network that might possess different scales.
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  [~, dzdy] = vl_argparsepos(struct(), varargin) ;

  sz = [1 1 1 1] ;
  sz(1:numel(size(x))) = size(x) ;
  multipliers = repmat(w, [sz(1:2) 1 size(x,4)]) ;

  %Set Local Response Normalization parameters
  kappa = 0 ; alpha = 1 ; beta = 0.5 ; N = 2 * size(x,3) ;
  normParams = [N kappa alpha beta] ;
  normalizedFeats = vl_nnnormalize(x, normParams) ;

  if isempty(dzdy)
    y = multipliers .* vl_nnnormalize(x, [N kappa alpha beta]) ;
  else
    dzdw_ = dzdy{1} .* normalizedFeats ;
    dzdw = sum(sum(sum(dzdw_,1), 2), 4) ;
    dzdx = vl_nnnormalize(x, normParams, multipliers .* dzdy{1}) ;
    y = dzdx ;
  end

repo: MFT


text:
function [u,p,w,edge,eqn,info] = StokesRT0(node,elem,bdFlag,pde,option)
%% STOKESRT0 Stokes equations: the lowest order Raviart-Thomas element in 2D.
%
%  [u,p,w,rotuI,Me,Mv,eqn,info] = StokesRT0(node,elem,bdFlag,pde,option)
%  uses the lowest order of Raviart-Thomas element (RT0) to approximate
%  velocity u and piecewise constant element (P0)  to approximate pressure
%  p, repectively.
%
%  We solve the following equation:
%       - grad div u + curl rot u + grad p  = f   in \Omega
%                                  - div u  = 0   in \Omega
%                                        u  = g   on \Gamma
%
% ifem StokesRT0doc
%
%  See also StokesBDM1B
%
% Based on a version by Ming Wang. Revised by Lin Zhong. Discussed with Jie
% Zhou and Long Chen. Further clean up by Long Chen.

if ~exist('option','var'), option = []; end

%% Data structure
elemunSort = elem;
[elem,bdFlag] = sortelem(elemunSort,bdFlag);
[elem2edge,edge] = dofedge(elem);
[Clambda,area,elemSign] = curlbasis(node,elem);
N = size(node,1); NT = size(elem,1); NE = size(edge,1);
Nu = NE; Np = NT; Ndof = Nu + Np;

%% Assemble matrix
tic; % record assemble time
% Mv: Lumped mass matrix for vertex: P1 element
vecMv = accumarray([elem(:,1);elem(:,2);elem(:,3)],[area;area;area]/3,[N,1]);
% invMv = spdiags(1./vecMv,0,N,N);
Mv = spdiags(vecMv,0,N,N);

% Me: Mass matrix for RT0 element
Me = getmassmatvec(elem2edge,area,Clambda,'RT0');

%invMt: the inverse of Mass matrix for P0 element
invMt = spdiags(1./area,0,NT,NT);

% B: negative divergence operator
B = -icdmat(double(elem2edge),elemSign*[1 -1 1]);

% C: curl operator
C = icdmat(double(edge),[-1 1]);

% R: weak rot operator
R = spdiags(1./vecMv,0,N,N)*C'*Me;

% Vector Laplacian
A = B'*invMt*B + R'*Mv*R;

%% Assemble right hand side
locEdge = [2,3; 1 3; 1,2]; % ascend ordering
fu = zeros(Nu,1);% the right hand side of u
g = zeros(Np,1); % the right hand side of p
if ~isfield(pde,'f') || (isfield(pde,'f') && isreal(pde.f) && all(pde.f==0))
    pde.f = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order is 3
end
if isfield(pde,'f') && ~isempty(pde.f)
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
    bt = zeros(NT,3);
    for p = 1:nQuad
        % quadrature points in the x-y-z coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        fp = pde.f(pxy);
        for k = 1:3
            i = locEdge(k,1); j = locEdge(k,2);
            % phi_k = lambda_iClambda_j - lambda_jClambda_i;
            phi_k = lambda(p,i)*Clambda(:,:,j)-lambda(p,j)*Clambda(:,:,i);
            rhs = dot(phi_k,fp,2);
            bt(:,k) = bt(:,k) + w(p)*rhs;
        end
    end
    bt = bt.*repmat(area,1,3);
    fu = accumarray(elem2edge(:),bt(:),[Nu 1]);
end
clear pxy fp bt rhs phi_k psi_k

%% Boundary condition and graddiv part of A
[u,p,ufreeDof,pDof,utbd] = getbdStokesRT0;
assembleTime = toc;

%% Solve the system of linear equations
% set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 1e5  % Direct solver for small size systems
        solver = 'direct';
    else             % Multigrid-type  solver for large size systems
        solver = 'mg';
    end
else
    solver = option.solver;
end
% solve the system
% get submatrices of ufreeDof
A0 = A(ufreeDof,ufreeDof);
B0 = B(:,ufreeDof);
f0 = fu(ufreeDof);
g0 = g;
if strcmp(solver,'direct') && ~isempty(ufreeDof)
    tic;
    bigA = [A0, B0'; ...
            B0, sparse(Np,Np)];
    bigF = [f0; g0];
    bigu = [u; p];
    bigFreeDof = [ufreeDof; Nu+pDof];
    bigu(bigFreeDof) = bigA(1:end-1,1:end-1)\bigF(1:end-1);
    u = bigu(1:Nu);
    p = bigu(Nu+1:end);
    info.solverTime = toc;
elseif strcmp(solver,'mg')
%   option.solver = 'vcycle';
    option.solver  = 'WCYCLE';
    [u(ufreeDof),p,info] = mgstokesRT0(A0,B0,f0,g0,u,p,node,elemunSort,ufreeDof,option);
end

%% Post-process
if length(pDof)~=Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(p.*area)/sum(area);
    p = p - c;
end
% w = invMv*(C'*Me*u + utbd);
w = R*u + utbd./vecMv;

%% Output information
eqn = struct('A',A0,'B',B0,'f',f0,'g',g0,'ufreeDof',ufreeDof,'pDof',pDof,'Me',Me,'Mv',Mv);
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdStokesRT0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [u,p,ufreeDof,pDof,utbd] = getbdStokesRT0
        %% Boundary condition of Stokes equation: RT0-P0 elements

        % Initial set up
        utbd = zeros(N,1); % the line integral of u on the boundary (u.t, tau)|_\partial \Omega
        u = zeros(Nu,1);
        p = zeros(Np,1);
        ufreeDof = (1:Nu)';
        pDof = (1:Np-1)';

        if ~exist('bdFlag','var'), bdFlag = []; end
        if ~isfield(pde,'g_D'), pde.g_D = []; end
        if ~isfield(pde,'g_N'), pde.g_N = []; end
        if ~isfield(pde,'g_R'), pde.g_R = []; end
        if isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
            bdFlag = [];
        end

        % Find Dirichlet boundary dof: fixedDof and pDof
        isFixedDof = false(Nu,1);
        if ~isempty(bdFlag)
            isDirichlet(elem2edge(bdFlag(:)==1)) = true;
            Dirichlet = edge(isDirichlet,:);
            isFixedDof(isDirichlet) = true;
            fixedDof = find(isFixedDof);
            ufreeDof = find(~isFixedDof);
        end

        % Set up edge sign
        % edgeSign records the inconsistency of asecond orientation and
        % induced orientation for each boundary edges
        edgeSign = ones(NE,1);
        idx = (bdFlag(:,1) ~= 0 ) & (elemSign == -1) ; % the first edge is on boundary
        edgeSign(elem2edge(idx,1)) = -1;
        idx = (bdFlag(:,2) ~= 0 ) & (elemSign ==  1) ; % the second edge is on boundary
        edgeSign(elem2edge(idx,2)) = -1;
        idx = (bdFlag(:,3) ~= 0 ) & (elemSign == -1) ; % the third edge is on boundary
        edgeSign(elem2edge(idx,3)) = -1;

        % Compute the boundary integral
        if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
            % else no bddof or g_D = 0 (no modification needed)
            % 1. Normal component of u is imposed strongly
            if (isnumeric(pde.g_D) && length(pde.g_D) == NE)
                u(fixedDof) = pde.g_D(fixedDof);
            else
                u(fixedDof) = faceinterpolate(pde.g_D,node,edge(fixedDof,:),'RT0');
            end
            % 2. Tangential component of u is imposed weakly
            [lambdagD,wgD] = quadpts1(3);
            nQuadgD = size(lambdagD,1);
            % quadratic bases 1--3--2
            bdphi = lambdagD;
            ve = node(Dirichlet(:,2),:) - node(Dirichlet(:,1),:);
            ge = zeros(size(Dirichlet,1),2);
            int_left = zeros(size(Dirichlet,1),2);
            int_right = zeros(size(Dirichlet,1),2);
            for pp = 1:nQuadgD
                ppxy = lambdagD(pp,1)*node(Dirichlet(:,1),:) ...
                     + lambdagD(pp,2)*node(Dirichlet(:,2),:);
                gDp = pde.g_D(ppxy);
                int_left = int_left + wgD(pp)*gDp*bdphi(pp,1);
                int_right = int_right + wgD(pp)*gDp*bdphi(pp,2);
            end
            ge(:,1) = dot(int_left,ve,2).*edgeSign(fixedDof);
            ge(:,2) = dot(int_right,ve,2).*edgeSign(fixedDof);
            utbd = accumarray(Dirichlet(:), [ge(:,1); ge(:,2)],[N,1]);
        end

        % Modify the right hand side
        fu = fu - A*u - Me*(C*(utbd./vecMv));
         g = g - B*u;
         g = g - mean(g);
    end
end

repo: Aero-matlab


text:
close all;
clear all;
clc;

data_file_path = 'bin/ra_mmv_phase_transition_snr_40db_s_8.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'ra_mmv_snr_40_db_s_8';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'Rank Aware, SNR=40dB s=8';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function [ x, y, z, w ] = ld0194 ( )

%*****************************************************************************80
%
%% LD0194 computes the 194 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(194,1);
  y = zeros(194,1);
  z = zeros(194,1);
  w = zeros(194,1);
  a = 0.0;
  b = 0.0;
  v = 0.1782340447244611E-02;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.5716905949977102E-02;
  [ n, x, y, z, w ] = gen_oh ( 2, n, a, b, v, x, y, z, w );
  v = 0.5573383178848738E-02;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.6712973442695226;
  v = 0.5608704082587997E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2892465627575439;
  v = 0.5158237711805383E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4446933178717437;
  v = 0.5518771467273614E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1299335447650067;
  v = 0.4106777028169394E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3457702197611283;
  v = 0.5051846064614808E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.1590417105383530;
  b = 0.8360360154824589;
  v = 0.5530248916233094E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end


text:
clear
clc
PATH = rand(3,3);
tau_vec = rand(size(PATH,1)-1,1)*2+1;
t = rand(10,1)*sum(tau_vec);
%t = 0;
K = length(tau_vec);
% Mapping matrix from p to equality constraint, b
A = A_const_appA(tau_vec);
A_yaw = augmentA_yaw(tau_vec);
% Permutation matrix to rearrange b
C = permut_mat(K);
C_yaw = permutMat_yaw(tau_vec);
% Cost matrix
Q = augmentQ(tau_vec,4);
Q_yaw = augmentQ(tau_vec,2);
% C*A^(-T)*Q*A^(-1)*C^(T) = R
R = comp_R_appA(A,C,Q);
R_yaw = comp_R_appA(A_yaw,C_yaw,Q_yaw);
% Known derivatives including position.
bF = const_bf(PATH);
bFyaw = bF_yaw(0, tau_vec);
% Unknown b terms are optimized.
b_sorted = comp_b_sorted(R, bF);
b_srtdYaw = b_srtd_yaw(R_yaw, bFyaw);
% C*A*P = b_sorted
% P = [p0 p1 ... p8 p9| ... ] where P(t) = p9t^9 + p8t^8 + ... + p0
P = (C*A)\b_sorted;
states_matlab = [];
for i=1:length(t)
    des_s = desired_state(tau_vec, t(i), PATH, P);
    ParticleState = [des_s.pos des_s.vel des_s.acc des_s.jerk]';
    Snap = des_s.snap';
    states_matlab = [states_matlab; ParticleState;Snap];
end

dlmwrite('t_matlab.csv',t, 'precision', '%2.16f')
dlmwrite('path_matlab.csv',PATH, 'precision', '%2.16f')
dlmwrite('tau_vec_matlab.csv',tau_vec, 'precision', '%2.16f')
dlmwrite('states_matlab.csv',states_matlab, 'precision', '%2.16f')

author: yorgoon

repo: minimum-snap-geometric-control


text: any(cellfun(@exist,{'arraydeal', 'chopdeal', 'cummax'})==3)

repo: REST


text:
classdef MatNatConfiguration < handle
    % MatNatConfiguration An object for providing details of the host, username, password etc for an XNAT server
    %
    %     Licence
    %     -------
    %     Part of MatNat. https://github.com/tomdoel/matnat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (Access = private)
        applicationDirectory
        serverName
        baseUrl
        userName
        password
    end

    methods
        function baseUrl = getBaseUrl(obj)
            baseUrl = obj.baseUrl;
        end

        function serverName = getServerName(obj)
            serverName = obj.serverName;
        end

        function userName = getUserName(obj)
            userName = obj.userName;
        end

        function password = getPassword(obj)
            password = obj.password;
        end

        function applicationDirectory = getApplicationDirectory(obj)
            applicationDirectory = obj.applicationDirectory;
        end

        function setBaseUrl(obj, baseUrl)
            obj.baseUrl = baseUrl;
        end

        function setServerName(obj, serverName)
            obj.serverName = serverName;
        end

        function setUserName(obj, userName)
            obj.userName = userName;
        end

        function setPassword(obj, password)
            obj.password = password;
        end

        function setApplicationDirectory(obj, applicationDirectory)
            obj.applicationDirectory = applicationDirectory;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% shell to start SCNlab BrainTools

global EXPT
global cl

% for timeseries viewer
global VOL
global f
global f2

hewma_gui;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [C_tri]=transitivity_bu(A)
%TRANSITIVITY_BU    Transitivity
%
%   T = transitivity_bu(A);
%
%   Transitivity is the ratio of 'triangles to triplets' in the network.
%   (A classical version of the clustering coefficient).
%
%   Input:      A       binary undirected connection matrix
%
%   Output:     T       transitivity scalar
%
%   Reference: e.g. Humphries et al. (2008) Plos ONE 3: e0002051
%
%
%   Alexandros Goulas, Maastricht University, 2010

    C_tri = trace(A^3) / (sum(sum(A^2)) - trace(A^2));

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function iI = iminvert(I)
%
% IMINVERT inverts the grayscaled input image in the range [0 255]

iI = 255 - I;


text:
function [label] = mne_read_label_file(filename)
%
% [label] = mne_read_label_file(filename)
%
% Reads a label file. The returned structure has the following fields
%
%     comment        comment from the first line of the label file
%     vertices       vertex indices (0 based, column 1)
%     pos            locations in meters (columns 2 - 4 divided by 1000)
%     values         values at the vertices (column 5)
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%

%
% This is based on the FreeSurfer read_label routine
% SUBJECTS_DIR environment variable is not consulted for the standard location
%

me='MNE:mne_read_label_file';
if(nargin ~= 1)
   error(me,'usage: mne_read_label_file(filename)');
end

[fid,message] = fopen(filename,'r');
if (fid < 0)
   error(me,'Cannot open file %s (%s)', filename,message);
end

comment = fgets(fid) ;
line = fgets(fid) ;
nv = sscanf(line, '%d') ;
data = fscanf(fid, '%d %f %f %f %f\n') ;
data = reshape(data, 5, nv);

for k = 2:length(comment)
   if comment(k) ~= ' '
      break;
   end
end
if comment(length(comment)) == 10
   comment = comment(1:end-1);
end
label.comment  = comment(k:end);
label.vertices = int32(data(1,:));
label.pos      = 1e-3*data(2:4,:)';
label.values   = data(5,:);
fclose(fid);


text:
function [f1,colors,axh] = tor_fig(varargin)
% function [f1,colors,axh] = tor_fig(varargin)
%
% Creates a figure, hold on, fontsize 18
%
% Var args are number of rows, then columns for subplot
fs=18;
if length(varargin)>2;,varargin{3}=fs;,end

f1 = figure('Color','w'); set(gca,'FontSize',18),hold on
colors = {'ro' 'go' 'bo' 'yo' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};

if length(varargin) > 0
    i = varargin{1};
    j = varargin{2};
    np = i * j;
    for k = 1:np
        axh(k) = subplot(i,j,k);
        set(gca,'FontSize',fs),hold on
    end
    axes(axh(1));
end

return


text:
% Fig. 2.15   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
close all;

g=9.81;     % m/sec^2
L=1;        % m
m=1;        % Kg
r2d=57.295; % radians to degrees

num = 1/(m*L^2);
den = [1 0 g/L];
t=0:.02:10;

y = step(num,den,t);  % output in radians
plot(t,r2d*y),grid
xlabel('Time (sec)')
ylabel('Pendulum angle \theta (deg)')
title('Fig. 2.15')
nicegrid


text:
function [ eddies ] = bottom_up_single(ssh_data, lat, lon, areamap, cyc, varargin)
%BOTTOM_UP_SINGLE Finds eddies using the Bottom Up method
%   Will return an array of struct's that contain the eddy data.
%   ssh_data: A 2D array of double's that contain the sea surface heights (latsxlons)
%   lat: A 1D array of double's that gives the latitude for a given index (dimension should match
%         that of ssh_data)
%   lon: A 1D array of double's that gives the longitude for a given index (dimension should match
%         that of ssh_data)
%   areamap: A 2D array that refer to the area of each pixel in SSH data (should have same size as ssh), or 1D array
%   that refer to area of each pixel for a specific lat in a regular grid (pixeld have same area for the same
%   latitude)
%   cyc: Pass 1 to output anticyclonic eddies or -1 to output cyclonic eddies
%   Optional parameters:
%   'minimumArea': minimum number of pixels for an eddy, used for validating eddies
%   'thresholdStep': the minimum step for thresholding, the unit is SSH's unit
%   'isPadding': whether or not to pad SSH data, should be true when scanning SSH data of the whole map. Set to false if
%   only partial SSH data is used.
%   'sshUnits': The units the SSH data is in. bottom_up_single is built to work natively on centimeter SSH data.
%   Valid parameters are 'meters' and 'centimeters'. If the paramater passed in is 'meters', the SSH data will
%   be multiplied by 100. No changes will be made if the paramater passed in is 'centimeters'.
%   The default value of 'sshUnits' is centimeters.
    p = inputParser;
    defaultMinPixelSize = 9;
    defaultThresholdStep = 0.05;
    defaultSSHUnits = 'centimeters';
    defaultPaddingFlag = true;
    addRequired(p, 'ssh_data');
    addRequired(p, 'lat');
    addRequired(p, 'lon');
    addRequired(p, 'areamap');
    addRequired(p, 'cyc');
    addParameter(p, 'minimumArea', defaultMinPixelSize, @isnumeric);
    addParameter(p, 'thresholdStep', defaultThresholdStep, @isnumeric);
    addParameter(p, 'isPadding', defaultPaddingFlag);
    addParameter(p, 'sshUnits', defaultSSHUnits);
    parse(p, ssh_data, lat, lon, areamap, cyc, varargin{:});
    minimumArea = p.Results.minimumArea;
    thresholdStep = p.Results.thresholdStep;
    isPadding = p.Results.isPadding;
    SSH_Units = p.Results.sshUnits;

    if strcmp(SSH_Units, 'meters')
        ssh_data = ssh_data * 100;
    elseif strcmp(SSH_Units, 'centimeters')
        max_val = max(ssh_data(:));
        min_val = max(ssh_data(:));
        if max_val < 1.5 && min_val > -1.5
            ssh_data = ssh_data * 100;
        elseif max_val < 150 && min_val > -150

        else
            disp('Could not figure out what units the SSH data provided is in. Running the scan assuming units of centimeters');
            disp('To specify SSH unit data, include an additional parameter of sshUnits, followed by the unit type, E.G. meters');
        end
    end

    %Check if the grid is regular (differences between lats and lons are equal)
    lat_diffs = lat(2:end) - lat(1:end-1);
    lat_diffs2 = lat_diffs(2:end) - lat_diffs(1:end-1);
    lon_diffs = lon(2:end) - lon(1:end-1);
    lon_diffs(lon_diffs <= -180) = lon_diffs(lon_diffs <= -180) + 360;
    lon_diffs(lon_diffs >= 180) = lon_diffs(lon_diffs >= 180) - 360;
    lon_diffs = abs(lon_diffs);
    lon_diffs2 = lon_diffs(2:end) - lon_diffs(1:end-1);
    if all(lat_diffs2 == 0) && all(lon_diffs2 == 0)
        % Regular grid, create a georasterref object to get eddy's centroid
        geo_raster_lat_limit = [lat(1) lat(end)];
        if lon(1) > lon(end)
            geo_raster_lon_limit = [lon(1) (360 + lon(end))];
        else
            geo_raster_lon_limit = [lon(1) lon(end)];
        end

        R = georasterref('LatLim', geo_raster_lat_limit, 'LonLim', geo_raster_lon_limit, 'RasterSize', ...
         size(ssh_data), 'ColumnsStartFrom', 'south', 'RowsStartFrom', 'west');
    else
        % Use normal indexing to get eddy's centroid
        R = [];
    end

    extrema = get_extrema(ssh_data, cyc);

    if isPadding
        origExtrema = extrema;
        extrema = [zeros(size(extrema, 1), 200), extrema, zeros(size(extrema, 1), 200)];
        sshExtended = [ssh_data(:, end-199:end), ssh_data(:, :), ssh_data(:, 1:200)];
        [extrema_lat_indexes, extrema_lon_indexes] = ind2sub(size(extrema), find(extrema == 1));

        extrema(:, 1:200) = origExtrema(:, end-199:end);
        extrema(:, end-199:end) = origExtrema(:, 1:200);
    else
        [extrema_lat_indexes, extrema_lon_indexes] = ind2sub(size(extrema), find(extrema == 1));
        sshExtended = ssh_data;
    end

    eddies = new_eddy();
    eddies(length(extrema_lat_indexes)).Date = NaN;
    cyc_sshExtended = sshExtended * cyc;
    parfor i = 1:length(extrema_lat_indexes)
        curr_lat_index = extrema_lat_indexes(i); curr_lon_index = extrema_lon_indexes(i);
        e = thresholdBU(cyc, curr_lat_index-5, curr_lat_index+5, curr_lon_index-5, curr_lon_index+5, ...
            sshExtended, extrema, curr_lat_index, curr_lon_index, sshExtended(curr_lat_index, curr_lon_index), ...
            thresholdStep, NaN, ...
            zeros(size(sshExtended)), lat, lon, R, areamap, minimumArea, isPadding, cyc_sshExtended);
        if ~isempty(e)
            eddies(i) = e;
        end
    end
    mask = cellfun('isempty', {eddies.Lat});
    eddies = eddies(~mask);
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [Dc,rho,dD,drho,d2psi] = MIspline(Tc,Rc,omega,m,varargin);
%
% Mutual Information based distance measure  (using a matlab implementation for the joint
% density, see MIspline for a C version)  using a general Gauss-Newton type framework.
%
% The MI distance measure is implemented using a Parzen-window estimator
% for the joint density rho. Given rho and drho, MI is computed as
% Dc    = psi(rho) = rho' * log(rho + tol) + ...
% dpsi  = log(rho + tol) + rho./(rho + tol) + ...
% d2psi = (rho + 2*tol)./(rho + tol)^2 + ...
% where the "..." indicates two more summands of similar type
%
% Example: run('MIcc')
%   setup2DhandData;
%   xc = getCellCenteredGrid(omega,m);
%   Tc = linearInter(dataT,omega,xc);
%   Rc = linearInter(dataR,omega,xc);
%   Dc = MIspline(Tc,Rc,omega,m);
%
% Input:
%  Tc, Rc            template and reference
%  omega, m     represents domain and discretization
%  varargin        optional parameters, e.g. doDerivative=0
%
% Output:
%  Dc           MI(Tc,Rc)
%  rho          joint density estimator
%  dD           dpsi*drrho
%  drho         derivative of joint density estimator
%  d2psi        log2(length(Tc))*sqrt(nT*nR);
%
% see also distances/contents
%==============================================================================

function [Dc,rho,dD,drho,d2psi] = MIspline(Tc,Rc,omega,m,varargin);

if nargin == 0,
  help(mfilename);
  setup2DhandData
  xc = getCellCenteredGrid(omega,m);
  Tc = linearInter(dataT,omega,xc);
  Rc = linearInter(dataR,omega,xc);
  D0 = feval(mfilename,Tc,Rc,omega,m);
  fprintf('%s: distance = %s\n',mfilename,num2str(D0));
  Dc = 'endOfMinimalExample';
  return;
end;

% initialize the output variables
Dc  = []; dD = []; rho  = []; drho = []; d2psi = [];

doDerivative = (nargout > 3);  % boolean for derivative computation

%setup default parameter for joint density estimator
tol  = 1e-6;  % tolerance for logarithm in entropy
minT = 0;     % (assumed) smallest value in Tc
maxT = 256;   % (assumed) largest value in Tc
nT   = 32;    % number of "bins" for Tc
minR = 0;     % (assumed) smallest value in Rc
maxR = 256;   % (assumed) largest value in Rc
nR   = 30;    % number of "bins" for Rc

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

widthT = (maxT-minT)/nT;
widthR = (maxR-minR)/nR;

% to ensure nothing is missed at the boundary, two artifical bins are added
minT  = minT-2*widthT;  maxT  = maxT+2*widthT;   nT = nT+4;
minR  = minR-2*widthR;  maxR  = maxR+2*widthR;   nR = nR+4;

% compute the Parzen-window estimator for the density and its derivative
% rho  is (nT*nR) x 1
% drho is (nT*nR) x length(Tc), large but sparse
[rho,drho] = rhoSpline(Tc,Rc,minT,maxT,nT,minR,maxR,nR);

% reformat rho for efficient computation of marginal densities
% rhoT and rhoR and undo formating
% note, summation can be described via matrix muliplication:
% rhoT = ST*rho, rhoR = SR*rho, where
% ST   = kron( eye(nT,nT), ones(1,nR) ) and
% SR   = kron( ones(1,nT), eye(nR,nR) );
% this is used for the computation of the derivative

rho  = reshape(rho,nT,nR);
rhoT = sum(rho,2);
rhoR = sum(rho,1)';
rho  = rho(:);

% compute MI
Dc = rhoT'*log(rhoT+tol) + rhoR'*log(rhoR+tol) - rho'*log(rho+tol);

if ~doDerivative, return; end;

% build the matrices ST, SR
ST    = sparse(kron(ones(1,nR),speye(nT,nT)));
SR    = sparse(kron(speye(nR,nR),ones(1,nT)));

dpsi  = (log(rhoT+tol)+rhoT./(rhoT+tol))'*ST ...
  +     (log(rhoR+tol)+rhoR./(rhoR+tol))'*SR ...
  -     (log(rho +tol)+rho ./(rho +tol))';

% apply chain rule for MI = psi(rho(Tc))
dD    = dpsi * drho;
d2psi = log2(length(Tc))*sqrt(nT*nR);
%==============================================================================


text:
function [ y1, m1, d1, h1, n1, s1, y2, m2, d2, h2, n2, s2 ] = ymdhms_swap ( ...
  y1, m1, d1, h1, n1, s1, y2, m2, d2, h2, n2, s2 )

%*****************************************************************************80
%
%% YMDHMS_SWAP swaps the data defining two YMDHMS dates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y1, M1, D1, H1, M1, S1, the first date.
%
%    Input/output, integer Y2, M2, D2, H2, M2, S2, the second date.
%
  tt = y1;
  y1 = y2;
  y2 = tt;

  tt = m1;
  m1 = m2;
  m2 = tt;

  tt = d1;
  d1 = d2;
  d2 = tt;

  tt = h1;
  h1 = h2;
  h2 = tt;

  tt = n1;
  n1 = n2;
  n2 = tt;

  tt = s1;
  s1 = s2;
  s2 = tt;

  return
end

repo: jburkardt-m


text:
function [x, infos] = recursive_nmu(V, rank, in_options)
% Recursive non-negative matrix underapproximation (Recursive-NMU).
%
% The problem of interest is defined as
%
%       min || V - WH ||_F^2,
%       where
%       {V, W, H} > 0 and WH <= V.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
%           Cnorm       Choice of the norm 1 or 2, default = 2.
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       N. Gillis and F. Glineur,
%       "Using Underapproximations for Sparse Nonnegative Matrix Factorization,"
%       Pattern Recognition 43 (4), pp. 1676-1687,
%       2010.
%
%       N. Gillis and R.J. Plemmons,
%       "Dimensionality Reduction, Classification, and Spectral Mixture Analysis
%       using Nonnegative Underapproximationm,"
%       Optical Engineering 50, 027001,
%       2011.
%
%
% This file is part of NMFLibrary.
%
% This file has been ported from
%       recursiveNMU.m at https://gitlab.com/ngillis/nmfbook/-/tree/master/algorithms
%       by Nicolas Gillis (nicolas.gillis@umons.ac.be)
%
% Ported by T.Fukunaga and H.Kasai on June 24, 2022 for NMFLibrary
%
% Change log:
%
%

    % set dimensions and samples
    [m, n] = size(V);
    % set local options
    local_options = [];
    local_options.Cnorm = 2;
    local_options.inner_max_epoch = 200;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;
    M = V;

    % initialize
    method_name = 'Recursive-NMU';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % select disp_freq
    disp_freq = set_disp_frequency(options);

    % initialize for this algorithm
    % (here)

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: k = 00, Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();
    prev_time = start_time;

    % main loop
    for k = 1 : rank

        % initialize epoch
        epoch = 0;

        % initialize (x,y) with an optimal rank-one NMF of M
        [w, s, h] = svds(M, 1);
        ws = abs(w) * sqrt(s);
        hs = abs(h) * sqrt(s);
        W(:, k) = ws;
        H(k, :) = hs';

        % initialize Lagrangian variable lambda
        R = M - ws * hs';
        lambda = max(zeros(size(R)), -R);

        % inner loop
        while (optgap > options.tol_optgap) && (epoch < options.inner_max_epoch)

            % update ws and hs
            A = M - lambda;

            if options.Cnorm == 1
                % l_1 norm minimization
                ws = max(0, (wmedian(A, hs)));
                hs = max(0, (wmedian(A', ws)));

             elseif options.Cnorm == 2
                % l_2 norm minimization
                ws = max(0, A * hs);
                ws = ws / (max(ws) + 1e-16);
                hs = max(0, (A' * ws) / (ws' * ws));
            end

            % update lambda
            if sum(ws) ~= 0 && sum(hs) ~= 0
                R = M - ws * hs';
                W(:, k) = ws;
                H(k, :) = hs';
                lambda = max(0, lambda - R / ((epoch + 1) + 1));
            else
                lambda = lambda / 2;
                ws = W(:, k);
                hs = H(k, :)';
            end

            % measure elapsed time
            elapsed_time = toc(start_time);

            % measure gradient calc count
            grad_calc_count = grad_calc_count + m*n;

            % update epoch
            epoch = epoch + 1;

            % store info
            % total iteration is computed as (k - 1) * options.inner_max_epoch + epoch
            W_rec = W(:, 1:k);
            H_rec = H(1:k, :);
            [infos, f_val, optgap] = store_nmf_info(V, W_rec, H_rec, [], options, infos, (k - 1) * options.inner_max_epoch + epoch, grad_calc_count, elapsed_time);

            % display infos
            if options.verbose > 2
                if ~mod(epoch, disp_freq)
                    fprintf('%s: k = %02d, Epoch = %04d, cost = %.16e, optgap = %.4e, time = %e\n', method_name, k, (k - 1) * options.inner_max_epoch + epoch, f_val, optgap, elapsed_time - prev_time);
                end
            end

        end

        M = max(0, M - ws * hs');

        % store info
        % total iteration is computed as (k - 1) * options.inner_max_epoch + epoch
        W_rec = W(:, 1:k);
        H_rec = H(1:k, :);
        [infos, f_val, optgap] = store_nmf_info(V, W_rec, H_rec, [], options, infos, (k - 1) * options.inner_max_epoch + epoch, grad_calc_count, elapsed_time);

        % display infos
        if options.verbose > 1
            if ~mod(epoch, disp_freq)
                fprintf('%s: k = %02d, Epoch = %04d, cost = %.16e, optgap = %.4e, time = %e\n', method_name, k, (k - 1) * options.inner_max_epoch + epoch, f_val, optgap, elapsed_time - prev_time);
            end
        end

        prev_time = elapsed_time;

    end

    if options.verbose > 0
        if optgap < options.tol_optgap
            fprintf('# Recursive-NMU: Optimality gap tolerance reached: f_val = %.4e < f_opt = %.4e (%.4e)\n', f_val, options.f_opt, options.tol_optgap);
        elseif (k - 1) * options.inner_max_epoch + epoch == rank * options.inner_max_epoch
            fprintf('# Recursive-NMU: Max epoch reached (%g).\n', rank * options.inner_max_epoch);
        end
    end

    x.W = W;
    x.H = H;

end

% WMEDIAN computes an optimal solution of
%
% min_x  || A - xy^T ||_1, y >= 0
%
% where A has dimension (m x n), x (m) and y (n),
% in O(mn log(n)) operations. Should be done in O(mn)...

function x = wmedian(A,y)

    % Reduce the problem for positive entries of y
    indi = y > 1e-16;
    A = A(:, indi);
    y = y(indi);
    [m, n] = size(A);
    A = A ./ repmat(y', m, 1);
    y = y / sum(y);

    % Sort rows of A, m*O(n log(n)) operations
    [As, Inds] = sort(A, 2);

    % Construct matrix of ordered weigths
    Y = y(Inds);

    % Extract the median
    actind = 1 : m;
    i = 1;
    sumY = zeros(m, 1);
    x = zeros(m, 1);
    while ~isempty(actind) % O(n) steps... * O(m) operations
        % sum the weitghs
        sumY(actind, :) = sumY(actind, :) + Y(actind, i);
        % check which weitgh >= 0
        supind = (sumY(actind, :) >= 0.5);
        % update corresponding x
        x(actind(supind)) = As(actind(supind), i);
        % only look reminding x to update
        actind = actind(~supind);
        i = i + 1;
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% Demo: OTVCA Hyperspectral Feature Extraction
%
load Indian_site
FN=16; % Number of features to extract
[FE]=OTVCA_V3(R,FN);
%% Extrected Features
figure(1)
subplot(1,3,1),imagesc(FE(:,:,1));colormap(gray);axis image;axis off;title('Feature 1');
subplot(1,3,2),imagesc(FE(:,:,2));colormap(gray);axis image;axis off;title('Feature 2');
subplot(1,3,3),imagesc(FE(:,:,3));colormap(gray);axis image;axis off;title('Feature 3');

author: BehnoodRasti


text:
function [V,X,Z,W] = spm_DEM_int(M,z,w,c)
% Integrates/evaluates a hierarchical model given innovations z{i} and w{i}
% FORMAT [V,X,Z,W] = spm_DEM_int(M,z,w,c);
%
% M{i}    - model structure
% z{i}    - innovations (causes)
% w{i}    - innovations (states)
% c{i}    - exogenous causes
%
% V{i}    - causal states (V{1} = y = response)
% X{i}    - hidden states
% Z{i}    - fluctuations (causes)
% W{i}    - fluctuations (states)
%
% The system is evaluated at the prior expectation of the parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_int.m 6132 2014-08-06 19:59:46Z karl $

% set model indices and missing fields
%--------------------------------------------------------------------------
M      = spm_DEM_M_set(M);

% innovations
%--------------------------------------------------------------------------
z = spm_cat(z(:)) + spm_cat(c(:));
w = spm_cat(w(:));

% number of states and parameters
%--------------------------------------------------------------------------
nt   = size(z,2);                        % number of time steps
nl   = size(M,2);                        % number of levels
nv   = sum(spm_vec(M.l));                % number of v (casual states)
nx   = sum(spm_vec(M.n));                % number of x (hidden states)

% order parameters (n= 1 for static models)
%==========================================================================
dt   = M(1).E.dt;                        % time step
n    = M(1).E.n + 1;                     % order of embedding
nD   = M(1).E.nD;                        % number of iterations per sample
td   = dt/nD;                            % integration time for D-Step

% initialize cell arrays for derivatives z{i} = (d/dt)^i[z], ...
%--------------------------------------------------------------------------
u.v      = cell(n,1);
u.x      = cell(n,1);
u.z      = cell(n,1);
u.w      = cell(n,1);
[u.v{:}] = deal(sparse(nv,1));
[u.x{:}] = deal(sparse(nx,1));
[u.z{:}] = deal(sparse(nv,1));
[u.w{:}] = deal(sparse(nx,1));

% hyperparameters
%--------------------------------------------------------------------------
ph.h   = {M.hE};
ph.g   = {M.gE};

% initialize with starting conditions
%--------------------------------------------------------------------------
vi     = {M.v};
xi     = {M.x};
u.v{1} = spm_vec(vi);
u.x{1} = spm_vec(xi);

% derivatives for Jacobian of D-step
%--------------------------------------------------------------------------
Dx    = kron(spm_speye(n,n,1),spm_speye(nx,nx,0));
Dv    = kron(spm_speye(n,n,1),spm_speye(nv,nv,0));
D     = spm_cat(spm_diag({Dv,Dx,Dv,Dx}));
dfdw  = kron(eye(n,n),eye(nx,nx));

% initialize conditional estimators of states to be saved (V and X)
%--------------------------------------------------------------------------
mnx   = 0;
mnv   = 0;
for i = 1:nl

    V{i} = sparse(M(i).l,nt);
    X{i} = sparse(M(i).n,nt);
    Z{i} = sparse(M(i).l,nt);
    W{i} = sparse(M(i).n,nt);

    % check for state-dependent precision
    %----------------------------------------------------------------------
    mnx = mnx | length(M(i).pg);
    mnv = mnv | length(M(i).ph);

end

% defaults for state-dependent precision
%--------------------------------------------------------------------------
Sz = 1;
Sw = 1;

% iterate over sequence (t) and within for static models
%==========================================================================
for t  = 1:nt
    for iD = 1:nD

        % Get generalised motion of random fluctuations
        %==================================================================

        % sampling time
        %------------------------------------------------------------------
        ts   = (t + (iD - 1)/nD)*dt;

        % evaluate state-dependent precision
        %------------------------------------------------------------------
        if mnx || mnv

            vi{nl} = vi{nl} + c{nl}(:,t);
            pu.x   = {spm_vec(xi(1:end - 1))};
            pu.v   = {spm_vec(vi(1 + 1:end))};
            p      = spm_LAP_eval(M,pu,ph);

            if mnv, Sz = sparse(diag(exp(-p.h/2))); end
            if mnx, Sw = sparse(diag(exp(-p.g/2))); end

        end

        % derivatives of innovations (and exogenous input)
        %------------------------------------------------------------------
        u.z  = spm_DEM_embed(Sz*z,n,ts,dt);
        u.w  = spm_DEM_embed(Sw*w,n,ts,dt);

        % Evaluate and update states
        %==================================================================

        % evaluate functions
        %------------------------------------------------------------------
        [u,dg,df] = spm_DEM_diff(M,u);

        % tensor products for Jacobian
        %------------------------------------------------------------------
        dgdv = kron(spm_speye(n,n,1),dg.dv);
        dgdx = kron(spm_speye(n,n,1),dg.dx);
        dfdv = kron(spm_speye(n,n,0),df.dv);
        dfdx = kron(spm_speye(n,n,0),df.dx);

        % Save realization
        %==================================================================
        vi   = spm_unvec(u.v{1},{M.v});
        xi   = spm_unvec(u.x{1},{M.x});
        zi   = spm_unvec(u.z{1},{M.v});
        wi   = spm_unvec(u.w{1},{M.x});
        if iD == 1
            for i = 1:nl
                if M(i).l, V{i}(:,t) = spm_vec(vi{i}); end
                if M(i).n, X{i}(:,t) = spm_vec(xi{i}); end
                if M(i).l, Z{i}(:,t) = spm_vec(zi{i}); end
                if M(i).n, W{i}(:,t) = spm_vec(wi{i}); end
            end
        end

        % break for static models
        %------------------------------------------------------------------
        if nt == 1, break, end

        % Jacobian for update
        %------------------------------------------------------------------
        J      = spm_cat({dgdv dgdx Dv  []  ;
                          dfdv dfdx []  dfdw;
                          []   []   Dv  []  ;
                          []   []   []  Dx});

        % update states u = {x,v,z,w}
        %------------------------------------------------------------------
        du     = spm_dx(J,D*spm_vec(u),td);

        % and unpack
        %------------------------------------------------------------------
        u      = spm_unvec(spm_vec(u) + du,u);

    end % iterations over iD

end % iterations over t

author: spm


text:
function [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% Applies Image Convolution in the Frequency Domain.
% Input:
%   - mI                -   Input Image.
%                           Structure: Matrix (numRows x numCols).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - mH                -   Filtering Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - convShape         -   Convolution Shape.
%                           Sets the convolution shape.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mI                -   Output Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'conv2()' - https://www.mathworks.com/help/matlab/ref/conv2.html.
% Remarks:
%   1.  For the transformation the (0, 0) of the image and the kernel is
%       the top left corner of teh matrix (i = 1, j = 1). Hence it can be
%       thought that a shifting is needed (To center the kernel).
%   1.  The function supports single channel image only.
% TODO:
%   1.
%   Release Notes:
%   -   1.0.000     29/04/2021  Royi Avital     RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONV_SHAPE_FULL     = 1;
CONV_SHAPE_SAME     = 2;
CONV_SHAPE_VALID    = 3;

numRows     = size(mI, 1);
numCols     = size(mI, 2);

numRowsKernel = size(mH, 1);
numColsKernel = size(mH, 2);

switch(convShape)
    case(CONV_SHAPE_FULL)
        numRowsFft  = numRows + numRowsKernel - 1;
        numColsFft  = numCols + numColsKernel - 1;
        firstRowIdx = 1;
        firstColIdx = 1;
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
    case(CONV_SHAPE_SAME)
        numRowsFft  = numRows + numRowsKernel;
        numColsFft  = numCols + numColsKernel;
        firstRowIdx = ceil((numRowsKernel + 1) / 2);
        firstColIdx = ceil((numColsKernel + 1) / 2);
        lastRowIdx  = firstRowIdx + numRows - 1;
        lastColdIdx = firstColIdx + numCols - 1;
    case(CONV_SHAPE_VALID)
        numRowsFft = numRows;
        numColsFft = numCols;
        firstRowIdx = numRowsKernel;
        firstColIdx = numColsKernel;
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
end

mO = ifft2(fft2(mI, numRowsFft, numColsFft) .* fft2(mH, numRowsFft, numColsFft), 'symmetric');
mO = mO(firstRowIdx:lastRowIdx, firstColIdx:lastColdIdx);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [x,state] = struct_poly(z,task,t)
%STRUCT_POLY Matrix with columns as polynomials.
%   [x,state] = struct_poly(z,[],t) computes a matrix x in which the
%   jth column is equal to the polynomial
%
%      polyval(z(j,:),s)
%
%   evaluated at the points s, defined as
%
%      (t-0.5*(min(t)+max(t)))/(0.5*(max(t)-min(t))).
%
%   The degree of the polynomial is equal to size(z,2)-1. The structure
%   state stores information which is reused in computing the right and
%   left Jacobian-vector products.
%
%   struct_poly(z,task,t) computes the right or left Jacobian-vector
%   product of this transformation, depending on the structure task. Use
%   the structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_rational, struct_rbf.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
if nargin < 3
    error('struct_poly:t','Please supply evaluation points.');
end

if isempty(task)
    [x,state] = struct_rational({z,[]},task,t);
elseif ~isempty(task.r)
    task.r = {task.r,[]};
    [x,state] = struct_rational({z,[]},task,t);
elseif ~isempty(task.l)
    [x,state] = struct_rational({z,[]},task,t); x = x{1};
end

end

author: andrewssobral

repo: lrslibrary