text:
% CircularMat test cases
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/5/2010
classdef CircularMat_Test < TestCase
    %constructor
    methods
        function this = CircularMat_Test(name)
            this = this@TestCase(name);
        end
    end

    %tests
    methods
        function testConstructor(this)
            import edu.stanford.covert.util.CircularMat;

            %0 arguments
            cmat = CircularMat();

            %1 argument
            mat = zeros(4,3,2);
            mat(4, 3, 1)=1;
            mat(3, 2, 2)=2;
            mat(1, 1, 2)=3;
            mat(3, 1, 1)=4;

            cmat = CircularMat(mat);

            %2 arguments
            mat = zeros(4,3,2);
            mat(4, 3, 1)=1;
            mat(3, 2, 2)=2;
            mat(1, 1, 2)=3;
            mat(3, 1, 1)=4;

            cmat = CircularMat(mat,1);
        end

        function testSubscriptReferenceAssignment(this)
            import edu.stanford.covert.util.CircularMat;

            %reference
            mat = zeros(100,2);
            mat(2,2)=1;
            cmat = CircularMat(mat,1);
            cmat2 = CircularMat(mat,2);
            assertEqual(CircularMat([0; 0; 0; 1]), cmat(99:102,2));
            assertEqual(CircularMat([0;0;0;1]), cmat(99:102,2));

            assertEqual(CircularMat(zeros(1,2)), cmat(1,:));
            assertEqual(CircularMat(zeros(1,2),2), cmat2(1,:));
            assertEqual(CircularMat(zeros(100,1),1), cmat(:,1));
            assertEqual(CircularMat([0;1;zeros(98,1)],1), cmat(:,2));
            assertEqual(CircularMat(zeros(100,1)), cmat(1:100,1));
            assertEqual(CircularMat([0;1;zeros(98,1)]), cmat(1:100,2));

            %assignment
            cmat = CircularMat(zeros(100,2),1);
            cmat(102,2)=1;
            assertEqual(CircularMat([0 0; 0 1; zeros(98,2)],1), cmat);

            cmat = CircularMat(zeros(100,2),1);
            cmat(100,3)=1;
            mat = zeros(100,2);
            mat(100,3)=1;
            assertEqual(CircularMat(mat,1), cmat);
        end

        function testEnd(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            assertEqual(CircularMat(mat(2:end,:,:,:)), cmat(2:end,:,:,:));
            assertEqual(CircularMat(mat(:,:,end,:)), cmat(:,:,end,:));
            assertEqual(CircularMat(mat(:,:,:,end-1)), cmat(:,:,:,end-1));
            assertEqual(CircularMat(mat(1:end)), cmat(1:end));
        end

        function testSize(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            %isempty
            assertEqual(false, isempty(cmat));
            assertEqual(true, isempty(CircularMat));

            %length
            assertEqual(4, length(cmat));

            %size
            assertEqual([4 1 3 2], size(cmat));
            assertEqual(4, size(cmat,1));
            assertEqual(1, size(cmat,5));

            %numel
            assertEqual(4*1*3*2, numel(cmat));

            %ndims
            assertEqual(4, ndims(cmat));

            %permute
            assertEqual(CircularMat(permute(mat,[2 4 1 3])), permute(cmat, [2 4 1 3]));
            assertEqual(CircularMat(permute(mat,[2 5 4 1 3])), permute(cmat, [2 5 4 1 3]));

            %transpose
            mat2 = permute(sum(mat,4),[1 3 2]);
            assertEqual(CircularMat(mat2'), transpose(CircularMat(mat2)));
            assertEqual(CircularMat(mat2'), CircularMat(mat2)');

            %ctranspose
            mat2 = permute(sum(mat,4),[1 3 2])+sqrt(-1);
            assertEqual(CircularMat(ctranspose(mat2)), ctranspose(CircularMat(mat2)));

            %reshape
            assertEqual(CircularMat(reshape(mat,[3 2 4])), reshape(CircularMat(mat),[3 2 4]));
            assertEqual(CircularMat(reshape(zeros(0,2),[0 3])), reshape(CircularMat(zeros(0,2)),[0 3]));

            %squeeze
            assertEqual(CircularMat(squeeze(mat)), squeeze(cmat));
        end

        function testFind(this)
            import edu.stanford.covert.util.CircularMat;

            %example 1
            mat = zeros(4,3);
            mat(4, 3)=1;
            mat(3, 2)=2;
            mat(1, 1)=3;
            mat(3, 1)=4;
            cmat = CircularMat(mat);

            [i1,j1,k1]=find(mat);
            [i2,j2,k2]=find(cmat);

            assertEqual(i1, i2);
            assertEqual(j1, j2);
            assertEqual(k1, k2);

            %example 2
            mat = zeros(4,3,'int32');
            mat(4, 3)=1;
            mat(3, 2)=2;
            mat(1, 1)=3;
            mat(3, 1)=4;
            cmat = CircularMat(mat);

            [i1,j1,k1]=find(mat);
            [i2,j2,k2]=find(cmat);

            assertEqual(i1, i2);
            assertEqual(j1, j2);
            assertEqual(k1, k2);
        end

        %concatenation
        function testConcatenation(this)
            import edu.stanford.covert.util.CircularMat;

            %Vertcat_nullIdentity
            assertEqual(...
                CircularMat,...
                vertcat(CircularMat));

            %Vertcat_identity
            t = CircularMat([1 2; 3 3], 1);
            assertEqual(t, vertcat(t));

            %Vertcat_threeTensors
            u = CircularMat([1 2; 3 3], 1);
            v = CircularMat([6 3; 4 2; 1 4], 1);

            assertEqual(CircularMat([...
                1 2;
                3 3;
                6 3;
                4 2;
                1 4;
                1 2;
                3 3]),...
                [u; v; u]);

            %Vertcat_nullFirst
            t = CircularMat([1 2; 3 3], 1);
            assertEqual(...
                CircularMat([1 2; 3 3]),...
                vertcat(CircularMat(), t));

            %Vertcat_nullLast
            t = CircularMat([1 2; 3 3], 1);
            assertEqual(...
                CircularMat([1 2; 3 3]),...
                vertcat(t, CircularMat()));

            %Vertcat_zeros
            assertEqual(...
                CircularMat(zeros(7,2)),...
                vertcat(...
                CircularMat(zeros(3,2)), CircularMat(zeros(4,2))));

            %Horzcat_nullIdentity
            assertEqual(...
                CircularMat,...
                horzcat(CircularMat));

            %Horzcat_identity
            t = CircularMat([1 2; 3 3], 1);
            assertEqual(t, horzcat(t));

            %Horzcat_nullFirst
            t = CircularMat([1 2; 3 3], 1);
            assertEqual(...
                CircularMat([zeros(2,3) [1 2; 3 3]]),...
                horzcat(CircularMat(zeros(2,3)), t));

            %Horzcat_nullLast
            t = CircularMat([1 2; 3 3], 1);
            assertEqual(...
                CircularMat([[1 2; 3 3] zeros(2,3)],1),...
                horzcat(t, CircularMat(zeros(2,3))));

            %Horzcat_zeros
            assertEqual(...
                CircularMat(zeros([2 7])),...
                horzcat(...
                CircularMat(zeros([2 3])), CircularMat(zeros([2 4]))));

            %padarray
            x = reshape(1:15,5,3);
            cmat = CircularMat(x);

            assertEqual(CircularMat(padarray(x,3)),padarray(cmat,3));
            assertEqual(CircularMat(padarray(x,[0 2])),padarray(cmat,[0 2]));
            assertEqual(CircularMat(padarray(x,[1 0 2])),padarray(cmat,[1 0 2]));
            assertEqual(CircularMat(padarray(x,[1 0 0])),padarray(cmat,[1 0 0]));
            assertEqual(CircularMat(padarray(x,[0 1 1],1)),padarray(cmat,[0 1 1],1));
            assertEqual(CircularMat(padarray(x,[0 1 1],2,'pre')),padarray(cmat,[0 1 1],2,'pre'));
            assertEqual(CircularMat(padarray(x,[0 1 1],2,'post')),padarray(cmat,[0 1 1],2,'post'));
            assertEqual(CircularMat(padarray(x,[0 1 1],2,'both')),padarray(cmat,[0 1 1],2,'both'));

            assertEqual(CircularMat(padarray(x,[2 2],'circular','pre')),padarray(cmat,[2 2],'circular','pre'));
            assertEqual(CircularMat(padarray(x,[2 2],'replicate','pre')),padarray(cmat,[2 2],'replicate','pre'));
            assertEqual(CircularMat(padarray(x,[2 2],'symmetric','pre')),padarray(cmat,[2 2],'symmetric','pre'));
            assertEqual(CircularMat(padarray(x,[2 2],'circular','post')),padarray(cmat,[2 2],'circular','post'));
            assertEqual(CircularMat(padarray(x,[2 2],'replicate','post')),padarray(cmat,[2 2],'replicate','post'));
            assertEqual(CircularMat(padarray(x,[2 2],'symmetric','post')),padarray(cmat,[2 2],'symmetric','post'));
            assertEqual(CircularMat(padarray(x,[2 2],'circular','both')),padarray(cmat,[2 2],'circular','both'));
            assertEqual(CircularMat(padarray(x,[2 2],'replicate','both')),padarray(cmat,[2 2],'replicate','both'));
            assertEqual(CircularMat(padarray(x,[2 2],'symmetric','both')),padarray(cmat,[2 2],'symmetric','both'));
        end

        function testAndOrXorNot(test)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            %mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            %not
            assertEqual(CircularMat(~mat), ~cmat);
            assertEqual(CircularMat(~mat2), ~cmat2);

            %or
            assertEqual(CircularMat(mat | 0), cmat | 0);
            assertEqual(CircularMat(mat | 2), cmat | 2);
            assertEqual(CircularMat(mat | Inf), cmat | Inf);
            assertEqual(CircularMat(mat | -Inf), cmat | -Inf);
            assertEqual(CircularMat(mat | mat2), cmat | cmat2);
            assertEqual(CircularMat(mat | mat2), cmat | mat2);

            assertEqual(CircularMat(0 | mat), 0 | cmat);
            assertEqual(CircularMat(2 | mat), 2 | cmat);
            assertEqual(CircularMat(Inf | mat), Inf | cmat);
            assertEqual(CircularMat(-Inf | mat), -Inf | cmat);
            assertEqual(CircularMat(mat2 | mat), cmat2 | cmat);
            assertEqual(CircularMat(mat2 | mat), mat2 | cmat);

            %and
            assertEqual(CircularMat(mat & 0), cmat & 0);
            assertEqual(CircularMat(mat & 2), cmat & 2);
            assertEqual(CircularMat(mat & Inf), cmat & Inf);
            assertEqual(CircularMat(mat & -Inf), cmat & -Inf);
            assertEqual(CircularMat(mat & mat2), cmat & cmat2);
            assertEqual(CircularMat(mat & mat2), cmat & mat2);

            assertEqual(CircularMat(0 & mat), 0 & cmat);
            assertEqual(CircularMat(2 & mat), 2 & cmat);
            assertEqual(CircularMat(Inf & mat), Inf & cmat);
            assertEqual(CircularMat(-Inf & mat), -Inf & cmat);
            assertEqual(CircularMat(mat2 & mat), cmat2 & cmat);
            assertEqual(CircularMat(mat2 & mat), mat2 & cmat);

            %xor
            assertEqual(CircularMat(xor(mat,0)), xor(cmat,0));
            assertEqual(CircularMat(xor(mat,2)), xor(cmat,2));
            assertEqual(CircularMat(xor(mat,Inf)), xor(cmat,Inf));
            assertEqual(CircularMat(xor(mat,-Inf)), xor(cmat,-Inf));
            assertEqual(CircularMat(xor(mat,mat2)), xor(cmat,cmat2));
            assertEqual(CircularMat(xor(mat,mat2)), xor(cmat,mat2));

            assertEqual(CircularMat(xor(0,mat)), xor(0,cmat));
            assertEqual(CircularMat(xor(2,mat)), xor(2,cmat));
            assertEqual(CircularMat(xor(Inf,mat)), xor(Inf,cmat));
            assertEqual(CircularMat(xor(-Inf,mat)), xor(-Inf,cmat));
            assertEqual(CircularMat(xor(mat2,mat)), xor(cmat2,cmat));
            assertEqual(CircularMat(xor(mat2,mat)), xor(mat2,cmat));
        end

        function testPlusMinus(test)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            %unary
            assertEqual(CircularMat(+mat), +cmat);
            assertEqual(CircularMat(-mat), -cmat);

            %addition
            assertEqual(CircularMat(mat+0), cmat+0);
            assertEqual(CircularMat(mat+2), cmat+2);
            assertEqual(CircularMat(mat+NaN), cmat+NaN);
            assertEqual(CircularMat(mat+Inf), cmat+Inf);
            assertEqual(CircularMat(mat+-Inf), cmat+-Inf);
            assertEqual(CircularMat(mat+mat2), cmat+cmat2);
            assertEqual(CircularMat(mat+mat2), cmat+mat2);

            assertEqual(CircularMat(0+mat), 0+cmat);
            assertEqual(CircularMat(2+mat), 2+cmat);
            assertEqual(CircularMat(NaN+mat), NaN+cmat);
            assertEqual(CircularMat(Inf+mat), Inf+cmat);
            assertEqual(CircularMat(-Inf+mat), -Inf+cmat);
            assertEqual(CircularMat(mat2+mat), cmat2+cmat);
            assertEqual(CircularMat(mat2+mat), mat2+cmat);

            %subtraction
            assertEqual(CircularMat(mat-0), cmat-0);
            assertEqual(CircularMat(mat-2), cmat-2);
            assertEqual(CircularMat(mat-NaN), cmat-NaN);
            assertEqual(CircularMat(mat-Inf), cmat-Inf);
            assertEqual(CircularMat(mat--Inf), cmat--Inf);
            assertEqual(CircularMat(mat-mat2), cmat-cmat2);
            assertEqual(CircularMat(mat-mat2), cmat-mat2);

            assertEqual(CircularMat(0-mat), 0-cmat);
            assertEqual(CircularMat(2-mat), 2-cmat);
            assertEqual(CircularMat(NaN-mat), NaN-cmat);
            assertEqual(CircularMat(Inf-mat), Inf-cmat);
            assertEqual(CircularMat(-Inf-mat), -Inf-cmat);
            assertEqual(CircularMat(mat2-mat), cmat2-cmat);
            assertEqual(CircularMat(mat2-mat), mat2-cmat);
        end

        function testMultiplyDivideExponentiate(test)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            mat3 = [
                2 Inf 0;
                0 2 -3;
                0 0 NaN];
            mat4 = [
                4 0 -7;
                -Inf 3 0
                NaN 0 0];
            cmat3 = CircularMat(mat3);
            cmat4 = CircularMat(mat4);

            %element-wise multiplication
            assertEqual(CircularMat(mat.*0), cmat.*0);
            assertEqual(CircularMat(mat.*2), cmat.*2);
            assertEqual(CircularMat(mat.*NaN), cmat.*NaN);
            assertEqual(CircularMat(mat.*Inf), cmat.*Inf);
            assertEqual(CircularMat(mat.*-Inf), cmat.*-Inf);
            assertEqual(CircularMat(mat.*mat2), cmat.*cmat2);
            assertEqual(CircularMat(mat.*mat2), cmat.*mat2);

            assertEqual(CircularMat(0.*mat), 0.*cmat);
            assertEqual(CircularMat(2.*mat), 2.*cmat);
            assertEqual(CircularMat(NaN.*mat), NaN.*cmat);
            assertEqual(CircularMat(Inf.*mat), Inf.*cmat);
            assertEqual(CircularMat(-Inf.*mat), -Inf.*cmat);
            assertEqual(CircularMat(mat2.*mat), cmat2.*cmat);
            assertEqual(CircularMat(mat2.*mat), mat2.*cmat);

            %matrix multiplication
            assertEqual(CircularMat(mat*0), cmat*0);
            assertEqual(CircularMat(mat*2), cmat*2);
            assertEqual(CircularMat(mat*NaN), cmat*NaN);
            assertEqual(CircularMat(mat*Inf), cmat*Inf);
            assertEqual(CircularMat(mat*-Inf), cmat*-Inf);
            assertEqual(CircularMat(mat3*mat4), cmat3*cmat4);
            assertEqual(CircularMat(mat3*mat4), cmat3*mat4);

            assertEqual(CircularMat(0*mat), 0*cmat);
            assertEqual(CircularMat(2*mat), 2*cmat);
            assertEqual(CircularMat(NaN*mat), NaN*cmat);
            assertEqual(CircularMat(Inf*mat), Inf*cmat);
            assertEqual(CircularMat(-Inf*mat), -Inf*cmat);
            assertEqual(CircularMat(mat4*mat3), cmat4*cmat3);
            assertEqual(CircularMat(mat4*mat3), mat4*cmat3);

            %element-wise right division
            assertEqual(CircularMat(mat./0), cmat./0);
            assertEqual(CircularMat(mat./2), cmat./2);
            assertEqual(CircularMat(mat./NaN), cmat./NaN);
            assertEqual(CircularMat(mat./Inf), cmat./Inf);
            assertEqual(CircularMat(mat./-Inf), cmat./-Inf);
            assertEqual(CircularMat(mat./mat2), cmat./cmat2);
            assertEqual(CircularMat(mat./mat2), cmat./mat2);

            assertEqual(CircularMat(0./mat), 0./cmat);
            assertEqual(CircularMat(2./mat), 2./cmat);
            assertEqual(CircularMat(NaN./mat), NaN./cmat);
            assertEqual(CircularMat(Inf./mat), Inf./cmat);
            assertEqual(CircularMat(-Inf./mat), -Inf./cmat);
            assertEqual(CircularMat(mat2./mat), cmat2./cmat);
            assertEqual(CircularMat(mat2./mat), mat2./cmat);

            %matrix right division
            % assertEqual(CircularMat(mat/0), cmat/0);
            % assertEqual(CircularMat(mat/2), cmat/2);
            % assertEqual(CircularMat(mat/NaN), cmat/NaN);
            % assertEqual(CircularMat(mat/Inf), cmat/Inf);
            % assertEqual(CircularMat(mat/-Inf), cmat/-Inf);
            % assertEqual(CircularMat(mat/mat2), cmat/cmat2);
            % assertEqual(CircularMat(mat/mat2), cmat/mat2);

            % assertEqual(CircularMat(0/mat), 0/cmat);
            % assertEqual(CircularMat(2/mat), 2/cmat);
            % assertEqual(CircularMat(NaN/mat), NaN/cmat);
            % assertEqual(CircularMat(Inf/mat), Inf/cmat);
            % assertEqual(CircularMat(-Inf/mat), -Inf/cmat);
            % assertEqual(CircularMat(mat2/mat), cmat2/cmat);
            % assertEqual(CircularMat(mat2/mat), mat2/cmat);

            %element-wise left division
            assertEqual(CircularMat(mat.\0), cmat.\0);
            assertEqual(CircularMat(mat.\2), cmat.\2);
            assertEqual(CircularMat(mat.\NaN), cmat.\NaN);
            assertEqual(CircularMat(mat.\Inf), cmat.\Inf);
            assertEqual(CircularMat(mat.\-Inf), cmat.\-Inf);
            assertEqual(CircularMat(mat.\mat2), cmat.\cmat2);
            assertEqual(CircularMat(mat.\mat2), cmat.\mat2);

            assertEqual(CircularMat(0.\mat), 0.\cmat);
            assertEqual(CircularMat(2.\mat), 2.\cmat);
            assertEqual(CircularMat(NaN.\mat), NaN.\cmat);
            assertEqual(CircularMat(Inf.\mat), Inf.\cmat);
            assertEqual(CircularMat(-Inf.\mat), -Inf.\cmat);
            assertEqual(CircularMat(mat2.\mat), cmat2.\cmat);
            assertEqual(CircularMat(mat2.\mat), mat2.\cmat);

            %matrix left division
            % assertEqual(CircularMat(mat\0), cmat\0);
            % assertEqual(CircularMat(mat\2), cmat\2);
            % assertEqual(CircularMat(mat\NaN), cmat\NaN);
            % assertEqual(CircularMat(mat\Inf), cmat\Inf);
            % assertEqual(CircularMat(mat\-Inf), cmat\-Inf);
            % assertEqual(CircularMat(mat\mat2), cmat\cmat2);
            % assertEqual(CircularMat(mat\mat2), cmat\mat2);

            % assertEqual(CircularMat(0\mat), 0\cmat);
            % assertEqual(CircularMat(2\mat), 2\cmat);
            % assertEqual(CircularMat(NaN\mat), NaN\cmat);
            % assertEqual(CircularMat(Inf\mat), Inf\cmat);
            % assertEqual(CircularMat(-Inf\mat), -Inf\cmat);
            % assertEqual(CircularMat(mat2\mat), cmat2\cmat);
            % assertEqual(CircularMat(mat2\mat), mat2\cmat);

            %element-wise exponentiation
            assertEqual(CircularMat(mat.^0), cmat.^0);
            assertEqual(CircularMat(mat.^2), cmat.^2);
            assertEqual(CircularMat(mat.^NaN), cmat.^NaN);
            assertEqual(CircularMat(mat.^Inf), cmat.^Inf);
            assertEqual(CircularMat(mat.^-Inf), cmat.^-Inf);
            assertEqual(CircularMat(mat.^mat2), cmat.^cmat2);
            assertEqual(CircularMat(mat.^mat2), cmat.^mat2);

            assertEqual(CircularMat(0.^mat), 0.^cmat);
            assertEqual(CircularMat(2.^mat), 2.^cmat);
            assertEqual(CircularMat(NaN.^mat), NaN.^cmat);
            assertEqual(CircularMat(Inf.^mat), Inf.^cmat);
            assertEqual(CircularMat(-Inf.^mat), -Inf.^cmat);
            assertEqual(CircularMat(mat2.^mat), cmat2.^cmat);
            assertEqual(CircularMat(mat2.^mat), mat2.^cmat);

            %matrix exponentiation
            mat = reshape(1:9,3,3);
            mat2 = reshape(10:18,3,3);
            cmat = CircularMat(mat);
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat^0), cmat^0);
            assertEqual(CircularMat(mat^2), cmat^2);
            assertEqual(CircularMat(mat^NaN), cmat^NaN);
            assertEqual(CircularMat(mat^Inf), cmat^Inf);
            assertEqual(CircularMat(mat^-Inf), cmat^-Inf);
            % assertEqual(CircularMat(mat^mat2), cmat^cmat2);
            % assertEqual(CircularMat(mat^mat2), cmat^mat2);

            assertEqual(CircularMat(0^mat), 0^cmat);
            assertEqual(CircularMat(2^mat), 2^cmat);
            assertEqual(CircularMat(NaN^mat), NaN^cmat);
            assertEqual(CircularMat(Inf^mat), Inf^cmat);
            assertEqual(CircularMat(-Inf^mat), -Inf^cmat);
            % assertEqual(CircularMat(mat2^mat), cmat2^cmat);
            % assertEqual(CircularMat(mat2^mat), mat2^cmat);
        end

        function testEquality(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat==0), cmat==0);
            assertEqual(CircularMat(mat==2), cmat==2);
            assertEqual(CircularMat(mat==NaN), cmat==NaN);
            assertEqual(CircularMat(mat==Inf), cmat==Inf);
            assertEqual(CircularMat(mat==-Inf), cmat==-Inf);
            assertEqual(CircularMat(mat==mat2), cmat==cmat2);
            assertEqual(CircularMat(mat==mat2), cmat==mat2);

            assertEqual(CircularMat(0==mat), 0==cmat);
            assertEqual(CircularMat(2==mat), 2==cmat);
            assertEqual(CircularMat(NaN==mat), NaN==cmat);
            assertEqual(CircularMat(Inf==mat), Inf==cmat);
            assertEqual(CircularMat(-Inf==mat), -Inf==cmat);
            assertEqual(CircularMat(mat2==mat), cmat2==cmat);
            assertEqual(CircularMat(mat2==mat), mat2==cmat);
        end

        function testInequality(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat~=0), cmat~=0);
            assertEqual(CircularMat(mat~=2), cmat~=2);
            assertEqual(CircularMat(mat~=NaN), cmat~=NaN);
            assertEqual(CircularMat(mat~=Inf), cmat~=Inf);
            assertEqual(CircularMat(mat~=-Inf), cmat~=-Inf);
            assertEqual(CircularMat(mat~=mat2), cmat~=cmat2);
            assertEqual(CircularMat(mat~=mat2), cmat~=mat2);

            assertEqual(CircularMat(0~=mat), 0~=cmat);
            assertEqual(CircularMat(2~=mat), 2~=cmat);
            assertEqual(CircularMat(NaN~=mat), NaN~=cmat);
            assertEqual(CircularMat(Inf~=mat), Inf~=cmat);
            assertEqual(CircularMat(-Inf~=mat), -Inf~=cmat);
            assertEqual(CircularMat(mat2~=mat), cmat2~=cmat);
            assertEqual(CircularMat(mat2~=mat), mat2~=cmat);
        end

        function testGreaterThan(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat>0), cmat>0);
            assertEqual(CircularMat(mat>2), cmat>2);
            assertEqual(CircularMat(mat>NaN), cmat>NaN);
            assertEqual(CircularMat(mat>Inf), cmat>Inf);
            assertEqual(CircularMat(mat>-Inf), cmat>-Inf);
            assertEqual(CircularMat(mat>mat2), cmat>cmat2);
            assertEqual(CircularMat(mat>mat2), cmat>mat2);

            assertEqual(CircularMat(0>mat), 0>cmat);
            assertEqual(CircularMat(2>mat), 2>cmat);
            assertEqual(CircularMat(NaN>mat), NaN>cmat);
            assertEqual(CircularMat(Inf>mat), Inf>cmat);
            assertEqual(CircularMat(-Inf>mat), -Inf>cmat);
            assertEqual(CircularMat(mat2>mat), cmat2>cmat);
            assertEqual(CircularMat(mat2>mat), mat2>cmat);
        end

        function testGreaterThanOrEqualTo(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat>=0), cmat>=0);
            assertEqual(CircularMat(mat>=2), cmat>=2);
            assertEqual(CircularMat(mat>=NaN), cmat>=NaN);
            assertEqual(CircularMat(mat>=Inf), cmat>=Inf);
            assertEqual(CircularMat(mat>=-Inf), cmat>=-Inf);
            assertEqual(CircularMat(mat>=mat2), cmat>=cmat2);
            assertEqual(CircularMat(mat>=mat2), cmat>=mat2);

            assertEqual(CircularMat(0>=mat), 0>=cmat);
            assertEqual(CircularMat(2>=mat), 2>=cmat);
            assertEqual(CircularMat(NaN>=mat), NaN>=cmat);
            assertEqual(CircularMat(Inf>=mat), Inf>=cmat);
            assertEqual(CircularMat(-Inf>=mat), -Inf>=cmat);
            assertEqual(CircularMat(mat2>=mat), cmat2>=cmat);
            assertEqual(CircularMat(mat2>=mat), mat2>=cmat);
        end

        function testLessThan(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat<0), cmat<0);
            assertEqual(CircularMat(mat<2), cmat<2);
            assertEqual(CircularMat(mat<NaN), cmat<NaN);
            assertEqual(CircularMat(mat<Inf), cmat<Inf);
            assertEqual(CircularMat(mat<-Inf), cmat<-Inf);
            assertEqual(CircularMat(mat<mat2), cmat<cmat2);
            assertEqual(CircularMat(mat<mat2), cmat<mat2);

            assertEqual(CircularMat(0<mat), 0<cmat);
            assertEqual(CircularMat(2<mat), 2<cmat);
            assertEqual(CircularMat(NaN<mat), NaN<cmat);
            assertEqual(CircularMat(Inf<mat), Inf<cmat);
            assertEqual(CircularMat(-Inf<mat), -Inf<cmat);
            assertEqual(CircularMat(mat2<mat), cmat2<cmat);
            assertEqual(CircularMat(mat2<mat), mat2<cmat);
        end

        function testLessThanOrEqualTo(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=5;
            mat2(3, 1, 2, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(3, 1, 3, 1)=4;
            cmat2 = CircularMat(mat2);

            assertEqual(CircularMat(mat<=0), cmat<=0);
            assertEqual(CircularMat(mat<=2), cmat<=2);
            assertEqual(CircularMat(mat<=NaN), cmat<=NaN);
            assertEqual(CircularMat(mat<=Inf), cmat<=Inf);
            assertEqual(CircularMat(mat<=-Inf), cmat<=-Inf);
            assertEqual(CircularMat(mat<=mat2), cmat<=cmat2);
            assertEqual(CircularMat(mat<=mat2), cmat<=mat2);

            assertEqual(CircularMat(0<=mat), 0<=cmat);
            assertEqual(CircularMat(2<=mat), 2<=cmat);
            assertEqual(CircularMat(NaN<=mat), NaN<=cmat);
            assertEqual(CircularMat(Inf<=mat), Inf<=cmat);
            assertEqual(CircularMat(-Inf<=mat), -Inf<=cmat);
            assertEqual(CircularMat(mat2<=mat), cmat2<=cmat);
            assertEqual(CircularMat(mat2<=mat), mat2<=cmat);
        end

        function testScalarFunctions(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            %abs
            assertEqual(CircularMat(abs(mat)), abs(cmat));

            %sign
            assertEqual(CircularMat(sign(mat)), sign(cmat));

            %sqrt
            assertEqual(CircularMat(sqrt(mat)), sqrt(cmat));

            %realsqrt
            assertEqual(CircularMat(realsqrt(abs(mat))), realsqrt(abs(cmat)));

            %exp
            assertEqual(CircularMat(exp(mat)), exp(cmat));

            %expm1
            assertEqual(CircularMat(expm1(mat)), expm1(cmat));

            %log
            assertEqual(CircularMat(log(mat)), log(cmat));

            %log1p
            assertEqual(CircularMat(log1p(mat)), log1p(cmat));

            %log2
            assertEqual(CircularMat(log2(mat)), log2(cmat));

            %log10
            assertEqual(CircularMat(log10(mat)), log10(cmat));

            %reallog
            assertEqual(CircularMat(reallog(abs(mat))), reallog(abs(cmat)));
        end

        function testComplexNumbers(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3+sqrt(-1);
            mat(3, 1, 1, 1)=4-sqrt(-1);
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            %abs
            assertEqual(CircularMat(abs(mat)), abs(cmat));

            %sign
            assertEqual(CircularMat(sign(mat)), sign(cmat));

            %real
            assertEqual(CircularMat(real(mat)), real(cmat));

            %imag
            assertEqual(CircularMat(imag(mat)), imag(cmat));

            %angle
            assertEqual(CircularMat(angle(mat)), angle(cmat));

            %conj
            assertEqual(CircularMat(conj(mat)), conj(cmat));

            %isreal
            assertEqual(isreal(mat), isreal(cmat));
        end

        function testSetOperators(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2=ones(4,1);
            cmat2=CircularMat(mat2);

            %unique
            assertEqual(CircularMat(unique(mat)), unique(cmat));
            assertEqual(CircularMat(unique(mat2)), unique(cmat2));
            assertEqual(CircularMat(unique(zeros(0,1))), unique(CircularMat(zeros(0,1))));
            assertEqual(CircularMat(unique(zeros(1,0))), unique(CircularMat(zeros(1,0))));
        end

        function testAnyAll(this)
            import edu.stanford.covert.util.CircularMat;

            fl_x=zeros(3,2);
            fl_x(1,3)=2;
            fl_x(2,:)=-1;

            st_x=CircularMat(fl_x);

            %all
            assertEqual(CircularMat(all(fl_x)), all(st_x));
            assertEqual(CircularMat(all(fl_x,1)), all(st_x,1));
            assertEqual(CircularMat(all(fl_x,2)), all(st_x,2));
            assertEqual(CircularMat(all(fl_x,3)), all(st_x,3));

            %any
            assertEqual(CircularMat(any(fl_x)), any(st_x));
            assertEqual(CircularMat(any(fl_x,1)), any(st_x,1));
            assertEqual(CircularMat(any(fl_x,2)), any(st_x,2));
            assertEqual(CircularMat(any(fl_x,3)), any(st_x,3));
        end

        function testReduction(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = zeros(4,1,3,2);
            mat2(4, 1, 3, 1)=-1;
            mat2(2, 1, 3, 2)=2;
            mat2(1, 1, 1, 2)=3;
            mat2(1, 1, 2, 1)=4;
            mat2(3, 1, 2, 2)=Inf;
            mat2(4, 1, 2, 1)=NaN;
            mat2(3, 1, 2, 1)=-Inf;
            cmat2=CircularMat(mat2);

            %min
            assertEqual(CircularMat(min(mat,0)), min(cmat,0));
            assertEqual(CircularMat(min(mat,2)), min(cmat,2));
            assertEqual(CircularMat(min(mat,-2)), min(cmat,-2));
            assertEqual(CircularMat(min(mat,NaN)), min(cmat,NaN));
            assertEqual(CircularMat(min(mat,Inf)), min(cmat,Inf));
            assertEqual(CircularMat(min(mat,-Inf)), min(cmat,-Inf));

            assertEqual(CircularMat(min(0,mat)), min(0,cmat));
            assertEqual(CircularMat(min(2,mat)), min(2,cmat));
            assertEqual(CircularMat(min(-2,mat)), min(-2,cmat));
            assertEqual(CircularMat(min(NaN,mat)), min(NaN,cmat));
            assertEqual(CircularMat(min(Inf,mat)), min(Inf,cmat));
            assertEqual(CircularMat(min(-Inf,mat)), min(-Inf, cmat));

            assertEqual(CircularMat(min(mat,mat2)), min(cmat,cmat2));
            assertEqual(CircularMat(min(mat,mat2)), min(cmat,mat2));
            assertEqual(CircularMat(min(mat2,mat)), min(cmat2,cmat));
            assertEqual(CircularMat(min(mat,mat2)), min(mat2,cmat));

            assertEqual(CircularMat(min(mat)), min(cmat));
            assertEqual(CircularMat(min(mat,[],1)), min(cmat,[],1));
            assertEqual(CircularMat(min(mat,[],2)), min(cmat,[],2));
            assertEqual(CircularMat(min(mat,[],3)), min(cmat,[],3));
            assertEqual(CircularMat(min(mat,[],4)), min(cmat,[],4));
            assertEqual(CircularMat(min(mat,[],5)), min(cmat,[],5));

            %max
            assertEqual(CircularMat(max(mat,0)), max(cmat,0));
            assertEqual(CircularMat(max(mat,2)), max(cmat,2));
            assertEqual(CircularMat(max(mat,-2)), max(cmat,-2));
            assertEqual(CircularMat(max(mat,NaN)), max(cmat,NaN));
            assertEqual(CircularMat(max(mat,Inf)), max(cmat,Inf));
            assertEqual(CircularMat(max(mat,-Inf)), max(cmat,-Inf));

            assertEqual(CircularMat(max(0,mat)), max(0,cmat));
            assertEqual(CircularMat(max(2,mat)), max(2,cmat));
            assertEqual(CircularMat(max(-2,mat)), max(-2,cmat));
            assertEqual(CircularMat(max(NaN,mat)), max(NaN,cmat));
            assertEqual(CircularMat(max(Inf,mat)), max(Inf,cmat));
            assertEqual(CircularMat(max(-Inf,mat)), max(-Inf, cmat));

            assertEqual(CircularMat(max(mat,mat2)), max(cmat,cmat2));
            assertEqual(CircularMat(max(mat,mat2)), max(cmat,mat2));
            assertEqual(CircularMat(max(mat2,mat)), max(cmat2,cmat));
            assertEqual(CircularMat(max(mat,mat2)), max(mat2,cmat));

            assertEqual(CircularMat(max(mat)), max(cmat));
            assertEqual(CircularMat(max(mat,[],1)), max(cmat,[],1));
            assertEqual(CircularMat(max(mat,[],2)), max(cmat,[],2));
            assertEqual(CircularMat(max(mat,[],3)), max(cmat,[],3));
            assertEqual(CircularMat(max(mat,[],4)), max(cmat,[],4));
            assertEqual(CircularMat(max(mat,[],5)), max(cmat,[],5));

            %range
            assertEqual(CircularMat(range(mat)), range(cmat));
            assertEqual(CircularMat(range(mat,1)), range(cmat,1));
            assertEqual(CircularMat(range(mat,2)), range(cmat,2));
            assertEqual(CircularMat(range(mat,3)), range(cmat,3));
            assertEqual(CircularMat(range(mat,4)), range(cmat,4));
            assertEqual(CircularMat(range(mat,5)), range(cmat,5));

            %sum
            assertEqual(CircularMat(sum(mat)), sum(cmat));
            assertEqual(CircularMat(sum(mat,1)), sum(cmat,1));
            assertEqual(CircularMat(sum(mat,2)), sum(cmat,2));
            assertEqual(CircularMat(sum(mat,3)), sum(cmat,3));
            assertEqual(CircularMat(sum(mat,4)), sum(cmat,4));
            assertEqual(CircularMat(sum(mat,5)), sum(cmat,5));

            %mean
            assertEqual(CircularMat(mean(mat)), mean(cmat));
            assertEqual(CircularMat(mean(mat,1)), mean(cmat,1));
            assertEqual(CircularMat(mean(mat,2)), mean(cmat,2));
            assertEqual(CircularMat(mean(mat,3)), mean(cmat,3));
            assertEqual(CircularMat(mean(mat,4)), mean(cmat,4));
            assertEqual(CircularMat(mean(mat,5)), mean(cmat,5));

            %var
            assertElementsAlmostEqual(var(mat), double(var(cmat)));
            assertElementsAlmostEqual(var(mat,1), double(var(cmat,1)));

            % assertElementsAlmostEqual(var(mat,0,1), double(var(cmat,0,1)));
            % assertElementsAlmostEqual(var(mat,0,2), double(var(cmat,0,2)));
            % assertElementsAlmostEqual(var(mat,0,3), double(var(cmat,0,3)));
            % assertElementsAlmostEqual(var(mat,0,4), double(var(cmat,0,4)));
            % assertElementsAlmostEqual(var(mat,0,5), double(var(cmat,0,5)));

            assertElementsAlmostEqual(var(mat,1,1), double(var(cmat,1,1)));
            assertElementsAlmostEqual(var(mat,1,2), double(var(cmat,1,2)));
            assertElementsAlmostEqual(var(mat,1,3), double(var(cmat,1,3)));
            assertElementsAlmostEqual(var(mat,1,4), double(var(cmat,1,4)));
            assertElementsAlmostEqual(var(mat,1,5), double(var(cmat,1,5)));

            %std
            assertEqual(CircularMat(std(mat)), std(cmat));
            assertEqual(CircularMat(std(mat,1)), std(cmat,1));

            % assertEqual(CircularMat(std(mat,0,1)), std(cmat,0,1));
            % assertEqual(CircularMat(std(mat,0,2)), std(cmat,0,2));
            % assertEqual(CircularMat(std(mat,0,3)), std(cmat,0,3));
            % assertEqual(CircularMat(std(mat,0,4)), std(cmat,0,4));
            % assertEqual(CircularMat(std(mat,0,5)), std(cmat,0,5));

            assertEqual(CircularMat(std(mat,1,1)), std(cmat,1,1));
            assertEqual(CircularMat(std(mat,1,2)), std(cmat,1,2));
            assertEqual(CircularMat(std(mat,1,3)), std(cmat,1,3));
            assertEqual(CircularMat(std(mat,1,4)), std(cmat,1,4));
            assertEqual(CircularMat(std(mat,1,5)), std(cmat,1,5));

            %norm
            mat = reshape(1:20,4,5);
            cmat = CircularMat(mat);
            assertEqual(norm(mat,1), norm(cmat,1));
            assertEqual(norm(mat,2), norm(cmat,2));
            assertEqual(norm(mat,Inf), norm(cmat,Inf));
            assertEqual(norm(reshape(mat,[],1),-Inf), norm(reshape(cmat,[],1),-Inf));
            assertEqual(norm(mat,'fro'), norm(cmat,'fro'));
        end

        function testIss(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            assertEqual(CircularMat(isnan(mat)), isnan(cmat));
            assertEqual(CircularMat(isinf(mat)), isinf(cmat));
            assertEqual(CircularMat(isfinite(mat)), isfinite(cmat));
        end

        %type casting
        function testTypeCasting(this)
            import edu.stanford.covert.util.CircularMat;

            mat = zeros(4,1,3,2);
            mat(4, 1, 3, 1)=-1;
            mat(3, 1, 2, 2)=2;
            mat(1, 1, 1, 2)=3;
            mat(3, 1, 1, 1)=4;
            mat(3, 1, 2, 1)=Inf;
            mat(1, 1, 2, 1)=NaN;
            mat(2, 1, 2, 1)=-Inf;
            cmat = CircularMat(mat);

            mat2 = sparse(permute(sum(mat,4),[1 3 2]));
            cmat2 = CircularMat(mat2);

            %sparse
            assertEqual(mat2, sparse(cmat2));

            %cast
            assertEqual(cast(mat,'double'), cast(cmat,'double'));
            assertEqual(cast(mat,'int32'), cast(cmat,'int32'));

            warning('off','MATLAB:nonIntegerTruncatedInConversionToChar')
            assertEqual(char(mat), char(cmat));
            warning('on','MATLAB:nonIntegerTruncatedInConversionToChar')

            exception=[];
            try
                assertEqual(logical(mat), logical(cmat));
            catch exception
            end
            if isempty(exception)
                throw(MException('CircularMat:typeConversion','NaN''s cannot be converted to logicals'));
            end
            assertEqual(int8(mat), int8(cmat));
            assertEqual(int16(mat), int16(cmat));
            assertEqual(int32(mat), int32(cmat));
            assertEqual(int64(mat), int64(cmat));
            assertEqual(uint8(mat), uint8(cmat));
            assertEqual(uint16(mat), uint16(cmat));
            assertEqual(uint32(mat), uint32(cmat));
            assertEqual(uint64(mat), uint64(cmat));
            assertEqual(single(mat), single(cmat));
            assertEqual(double(mat), double(cmat));
        end
    end
end


text:
function [type,beats,beatsuct] = beat_class(ecg,qrs,T_LENGTH)
% this function is used to contruct a template ecg based on the location of
% the R-peaks. A series of peaks that match with each other are stacked to
% build a template. This template can then be used for ecg morphological
% analysis or further processing. Note that the qrs location inputed must
% be as precise as possible. This approach for building the template ECG
% seems to be the best of the alternatives that were tested and leaves the
% freedom of having more than one mode (i.e. multiple ECG template can be
% built if there are different cycle morphology such as PVC)  but it is not
% particularly fast.
% The procedure for building the template is:
% 1. create average wrapped template
% 2. identify different modes present
%
% inputs
%   ecg:            the ecg channel(s)
%   qrs:            qrs location [number of samples]
%   T_LENGTH        template length
%
% --
% ECG classification from single-lead segments using Deep Convolutional Neural
% Networks and Feature-Based Approaches - December 2017
%
% Released under the GNU General Public License
%
% Copyright (C) 2017  Fernando Andreotti, Oliver Carr
% University of Oxford, Insitute of Biomedical Engineering, CIBIM Lab - Oxford 2017
% fernando.andreotti@eng.ox.ac.uk
%
%
% For more information visit: https://github.com/fernandoandreotti/cinc-challenge2017
%
% Referencing this work
%
% Andreotti, F., Carr, O., Pimentel, M.A.F., Mahdi, A., & De Vos, M. (2017).
% Comparing Feature Based Classifiers and Convolutional Neural Networks to Detect
% Arrhythmia from Short Segments of ECG. In Computing in Cardiology. Rennes (France).
%
% Last updated : December 2017
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% == manage inputs
if nargin<2; error('ecg_template_build: wrong number of input arguments \n'); end;

if size(ecg,1)>size(ecg,2), ecg = ecg';end

% == constants
% qrs(1) = []; qrs(end) = [];                                % remove extremity peaks
extremities = (qrs <= 2*T_LENGTH | qrs >= length(ecg)-2*T_LENGTH);        % test if there are peaks on the border that may lead to error
qrs = round(qrs(~extremities)); % remove extremity peaks

%% Phase wrap to get all within 2*pi
NB_BINS = 250; % number of bins for wrapping
[M,RES] = stackbeats(ecg,qrs,T_LENGTH,NB_BINS);

% beatsuct = median(M')';
% dem = round((length(beatsuct)-2*T_LENGTH-1)/2);
% beats{1} =  beatsuct(dem+1:end-dem);
% type = ones(length(qrs),1);

%% Clustering
Nclus1 = 5;
Nclus2 = 2;
M = bsxfun(@times,M,hamming(size(M,1)));  % do it twice, it's the
M = bsxfun(@times,M,hamming(size(M,1)));  % magic sauce
M = bsxfun(@times,M,hamming(size(M,1)));  % magic sauce
M = M';
% [~, clusters, ~] = fkmeans(M',Nclus);
rng(1);
[~,clusters,~] = kmeans(M,Nclus1); % may need to provide centroid..
[~,clusters,~] = kmeans(clusters,Nclus2); % may need to provide centroid..
clusters = clusters';

% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% TODO: allow more than two clusters!!!!!!!!!!!!!!!!!!!!
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

corel = corrcoef(clusters(50:200,1),clusters(50:200,2));

%% In case beats are correlated to each other
if corel(1,2) > 0.7
    type = ones(size(qrs));
    beats{1} = resample(mean([clusters(:,1), clusters(:,2)],2),RES,NB_BINS);
    dem = (length(beats{1})-2*T_LENGTH-1)/2;
    beatsuct = beats;
    beats = cellfun(@(x) x(dem+1:end-dem),beats,'UniformOutput',0);
    return;
end

%% In case one cluster has low amplitude

idx = var(clusters)<0.1*var(ecg); % remove clusters with variance less than 10% of signal
if any(idx)
   type = ones(size(qrs));
   beats{1} = resample(clusters(:,~idx),RES,NB_BINS);
   dem = round((length(beats{1})-2*T_LENGTH-1)/2);
   beatsuct = beats;
   beats = cellfun(@(x) x(dem+1:end-dem),beats,'UniformOutput',0);
   return;
end

%% Classifying beats

beats{1} = resample(clusters(:,1),RES,NB_BINS);
beats{2} = resample(clusters(:,2),RES,NB_BINS);
beatsuct = beats;
dem = (length(beats{1})-2*T_LENGTH-1)/2;
beats = cellfun(@(x) x(dem+1:end-dem),beats,'UniformOutput',0);

for i = 1:length(qrs)
    cor = corrcoef(ecg(qrs(i)-T_LENGTH:qrs(i)+T_LENGTH),beats{1});
    corel(1,i) = cor(1,2);
    cor = corrcoef(ecg(qrs(i)-T_LENGTH:qrs(i)+T_LENGTH),beats{2});
    corel(2,i) = cor(1,2);
end
[~,type]=max(corel);

if sum(type==2)>sum(type==1)
    type(type==2) = 3;
    type(type==1) = 2;
    type(type==3) = 2;
end

if all(type==2)
    type = ones(size(type));
end
type = type';

end

author: fernandoandreotti


text:
function [mu, varsigma] = modelPosteriorMeanVar(model, X);

% MODELPOSTERIORMEANVAR Mean and variances of the posterior at points given by X.
% FORMAT
% DESC returns the posterior mean and variance for a given set of
% points.
% ARG model : the model for which the posterior will be computed.
% ARG x : the input positions for which the posterior will be
% computed.
% RETURN mu : the mean of the posterior distribution.
% RETURN sigma : the variances of the posterior distributions.
%
% SEEALSO : modelCreate
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS

  fhandle = str2func([model.type 'PosteriorMeanVar']);
  if str2num(version('-release'))>13
    [mu, varsigma] = fhandle(model, X);
  else
    [mu, varsigma] = feval(fhandle, model, X);
  end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% problem 6 - Computation of linear convolution

% a)
x1=[1 2 3 4];
x2=[ 5 4 3 2 1];
y=conv(x1,x2)

% b)
N1=length(x1);
N2=length(x2);
M=N1+N2-1;
x1(M)=0;
x2(M)=0;
y=circonv(x1,x2)

% c)
X1=fft(x1,M);
X2=fft(x2,M);
y=ifft(X1.*X2)

author: Sable

repo: mcbench-benchmarks


text:
close all
clear

%------------------------------------------------------------------------------------------
nb = 1000;%Total number of initial states
nmpcd=zeros(1,nb);%Store computing time for one time-step for NMPC-DCBF
impcd=zeros(1,nb);%Store computing time for one time-step for iMPC-DCBF
nmpcinf=zeros(1,1);%Store infeasible rate for one time-step for NMPC-DCBF
impcinf=zeros(1,1);%Store infeasible rate for one time-step for iMPC-DCBF
i=8;%i is number of horozon
ts = 1;%Total time steps
load(gamma1_0p6gamma2_0p6\InitialStateData');%Load the generated 1000 initial states
[a, b, c, d]=compare(i,ts,nb,InitialMat);%Involves NMPC-DCBF and iMPC-DCBF
nmpcd(1,:)=a;
impcd(1,:)=b;
nmpcinf(1,1)=c;
impcinf(1,1)=d;
nmpcplot1=nmpcd(1,:);
nmpcplot1(nmpcplot1==-1)=[];%Eliminate infeasible one time-step trajectories
impcplot1=impcd(1,:);
impcplot1(impcplot1==-1)=[];%Eliminate infeasible one time-step trajectories

save('timecom8','nmpcplot1','impcplot1');
save('feasibility8','nmpcinf','impcinf');

distnmpc = fitdist(nmpcplot1','Normal');
distimpc = fitdist(impcplot1','Normal');
mu1=distnmpc.mu;%Get mean of sample of computing time for NMPC-DCBF
mu2=distimpc.mu;%Get mean of sample of computing time for iMPC-DCBF
sigma1=distnmpc.sigma;%Get variance of sample of computing time for NMPC-DCBF
sigma2=distimpc.sigma;%Get variance of sample of computing time for iMPC-DCBF

%Initialize atmosphere parameters
function [tnmpc, timpc, ratiotnmpc, ratiotimpc]=compare(N11,ttsim,samplen,InitialMat)
tnmpc=[];
timpc=[];
for i=1:samplen
N1=N11;
tsim=ttsim;
xini1=InitialMat(1,i);
yini1=InitialMat(2,i);
thetaini1=InitialMat(3,i);
vini1=InitialMat(4,i);
x01=[xini1;yini1;thetaini1;vini1];
x02=[xini1 yini1 thetaini1 vini1];
t1 = nmpcdcbf(x01, N1, tsim);
t2 = impcdcbf(x02, N1, tsim);
tindex=[N11 i];
disp(tindex);
tnmpc=[tnmpc t1];%Computing time for NMPC-DCBF
timpc=[timpc t2];%Computing time for iMPC-DCBF
end
nnmpc1 = length(tnmpc);
nimpc1 = length(timpc);
tnmpcs = tnmpc;
timpcs = timpc;
tnmpcs(tnmpcs==-1)=[];
timpcs(timpcs==-1)=[];
nnmpc2 = length(tnmpcs);
nimpc2 = length(timpcs);
ratiotnmpc = (nnmpc1-nnmpc2)/nnmpc1;%Infeasible rate for NMPC-DCBF
ratiotimpc = (nimpc1-nimpc2)/nimpc1;%Infeasible rate for iMPC-DCBF
end

function t1 = nmpcdcbf(x00, N1, ttsim)
%% General Flags
run_nmpc_dcbf_one = true;
run_nmpc_dcbf_multiple = true;

%% Setup and Parameters
x0 = x00;
dt = 0.1;
time_total = ttsim *dt;
P = zeros(4,4);%Weight matrix P
P(1,1) = 10; P(2,2) = 10; P(3,3) = 10;P(4,4) = 10;
Q = zeros(4,4);%Weight matrix Q
Q(1,1) = 10; Q(2,2) = 10; Q(3,3) = 10;Q(4,4) = 10;
R = zeros(4,4);%Weight matrix R
R(1,1) = 1; R(2,2) = 1;R(3,3) = 1000;R(4,4) = 1000;
%Variables range as below
xmin = [-10; -10; -10;-10];
xmax = [10; 10; 10;10];
umin = [-7; -5;-inf;-inf];
umax = [7; 5;inf;inf];

%% Discrete-time unicycle model
system.dt = dt;
system.A = [1 0 0 0;
    0 1 0 0;
    0 0 1 0;
    0 0 0 1];
system.B = [x0(4,1)*cos(x0(3,1))*dt;
    x0(4,1)*sin(x0(3,1))*dt;
    0;
    0];
system.C =[0 0 0 0;
    0 0 0 0;
    1*dt 0 0 0;
    0 1*dt 0 0];
system.xl = xmin;
system.xu = xmax;
system.ul = umin;
system.uu = umax;

%% NMPC-DCBF parameters
params_nmpc_dcbf.Q = Q;
params_nmpc_dcbf.R = R;
params_nmpc_dcbf.P = P;
params_nmpc_dcbf.gamma1 = 0.6;
params_nmpc_dcbf.gamma2 = 0.6;

%% Obstacle
obs.pos1 = [0; 0];
obs.r1 = 1;

%% Simulate NMPC-DCBF with other N
params_nmpc_dcbf.N = N1;
if run_nmpc_dcbf_one
    fprintf('Run NMPC-DCBF-8\n');
    controller_nmpc_dcbf_multiple8 = NMPCDCBF2(x0, system, params_nmpc_dcbf);%Highest order mcbf=2
    controller_nmpc_dcbf_multiple8.obs = obs;
    controller_nmpc_dcbf_multiple8.sim(time_total);
end
t1=controller_nmpc_dcbf_multiple8.tt;
end

function t2=impcdcbf(x00, N1, ttsim)
t2=0;
xo = 0;
yo = 0;
r = 1;
N = N1;%Number of horizon
K1 = 1000;%Maximum iteration times, jmax
K = 1000;%Maximum iteration times, jmax
dt = 0.1;
gamma1 = 6;
gamma2 = 6;
nsim = ttsim;

tic;

% Constraints

umin = [-7; -5; -inf; -inf;];
umax = [7; 5; Inf; Inf];
xmin = [-10; -10; -10; -10];
xmax = [ 10;  10;  10;  10];

% Objective function
Q = diag([10 10 10 10]);
QN = Q;
R = 1 * eye(4);
R(3,3) = 1000;
R(4,4) = 1000;
% Initial and reference states
x0 = x00;
xr = [3; 0.01; 0; 0];
ur = [0; 0; 1; 1];

% Dynamic system initialization
BB = [0 0 0 0;0 0 0 0;dt 0 0 0;0 dt 0 0];

% Convex MPC frame
[nx, nu] = size(BB);
% Initialize states set (x00,x01,,,x0N）
x_0 = [];
x0 = x0';
u_0 = zeros(nu, N);
u0 = zeros(nu, 1);% we may need to change this for a better warm-up start
impc = zeros(nx, nsim + 1);
impc(:, 1) = x0;
x0new = x0;
for i=1 : (N+1)
    x_0 = [x_0 x0new];
    x0new = [dt*x0new(4)*cos(x0new(3))+x0new(1);dt*x0new(4)*sin(x0new(3))+x0new(2);dt*u0(1)+x0new(3);dt*u0(2)+x0new(4)];
end

abc = tangentline(x_0, r);
Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);%First order CBF constraints
lcbf = -abc(3, :);
lcbf(1) = [];
lcbf = lcbf';
ucbf = inf * ones(N, 1);
A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);%Second order CBF constraints
lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
ucbf2 = inf * ones(N-1, 1);

% Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
% - quadratic objective
P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
% - linear objective
q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
% - linear dynamics
Ax = getax(x_0, dt, N, nx);
Bu = kron([sparse(1, N); speye(N)], BB);
Cx = getcx(x_0, u_0, dt, N, nx);
Aeq = [Ax, Bu];
leq = [-x0; -Cx];
ueq = leq;
% - input and state constraints
Aineq = speye((N+1)*nx + N*nu);
lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

% - OSQP constraints
A = [Aeq; Aineq; Acbf; A2cbf];
l = [leq; lineq; lcbf; lcbf2];
u = [ueq; uineq; ucbf; ucbf2];

% Create an OSQP object
prob = osqp;
% Setup workspace
prob.setup(P, q, A, l, u, 'warm_start', true);
% Solve
res = prob.solve();

% Check solver status
if ~strcmp(res.info.status, 'solved')
%    error('OSQP did not solve the problem!')
        t2=-1;
        return
end
ctrlx = res.x(1:(N+1)*nx);
ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);

for j = 1 : (K1-1)%Iterations for the first time-step
    storagex = ctrlx;%store the state variables and inputs in order to compare them with the next iteration
    storageu = ctrlu;
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
end % Move for the first step
ctrl = ctrlu(1:nu);
x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
ctrlu = rewrite(ctrlu, nu, N);
t2=t2+toc;
return

x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
u_0 = transu(ctrlu, nu, N);
impc(:, 2) = x0;
storagex = ctrlx;
storageu = ctrlu;

for i = 1 : (nsim-1)%Iterations for the left time steps

   for j = 1 : K
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    x0 = ctrlx(1:nx);
    testx0 = [testx0 x0];
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
    storagex = ctrlx;
    storageu = ctrlu;
   end
   ctrl = ctrlu(1:nu);
   x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
   ctrlu = rewrite(ctrlu, nu, N);
   x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
   u_0 = transu(ctrlu, nu, N);
   impc(:, i+2) = x0;
end
end

% Linerize the CBF constraints (get a, b, c for lines)
function abc = tangentline(xy, r)% x and y from initialize states set, abc are coeeficients for linear equation a*x+b*y+c=0
[xx, ~] = size(xy);%xx=2,yy=N+1
xy(xx,:) = []; % this part should be changed for other case
xy((xx-1),:) = []; % this part should be changed other case
[xx, yy] = size(xy);%xx=2,yy=N+1
xyjiao = zeros(xx, yy);%intersection points
for i = 1 : xx
    for j = 1 : yy
        xyjiao(i, j) = r * xy(i, j) * (1 / (xy(:, j)' * xy(:, j)))^(0.5);%calculate coordinates of intersection points
    end
end
cc = -r^2 * ones(1, yy);
abc = [xyjiao; cc];
end

% Get CBF constraints matrix
function Acbf = cbfmat(abc, nx, nu, dt, gamma, x0)
[~, yy] = size(abc);
Acbfx = zeros((yy-1), yy*nx);
Acbfu = zeros((yy-1), (yy-1)*nu);
for i = 1 : (yy-1)
    Acbfx(i, (i*nx)+1) = abc(1, (i+1));
    Acbfx(i, (i*nx)+2) = abc(2, (i+1));
end
for i = 1 : (yy-1)
    Acbfu(i, ((i-1)*nu+3)) = - (1 - dt * gamma)^(i) * (abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
Acbf = [Acbfx Acbfu];
end

% Transfer vector x into matrix
function res = trans(x0, vector, nxx, nyy)%nxx=nx,nyy=N+1
 res = zeros(nxx, nyy);
 res(:,1) = x0;
 for i = 1 : (nyy -1)
     res(:,i+1)= vector(((i)*nxx+1):(i+1)*nxx);
 end
end

% Transfer vector u into matrix
function resu = transu(vector, nxx, nyy)%nxx=nu,nyy=N
 resu = zeros(nxx, nyy);
 for i = 1 : (nyy)
     resu(:,i)= vector(((i-1)*nxx+1):(i)*nxx);
 end
end

% Rewrite u vector
function reu = rewrite(vector, nu, N)
append = vector((N-1)*nu+1:N*nu);
vector(1:nu) = [];
reu = [vector;append];
end

% Get new x_0
function x_0 = newinit(x0, ctrlu, nx, nu, N, dt)
 x_0 = zeros(nx, N+1);
 x_0(:, 1) = x0;
 for i=1 : N
    u0 = ctrlu((i-1)*nu+1:i*nu);
    x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*u0(1)+x0(3);dt*u0(2)+x0(4)];
    x_0(:, i + 1) = x0;
 end
end

% Get AA matrix
function AA = getaa(x0, dt)
 AA = [1 0 -x0(4)*sin(x0(3))*dt cos(x0(3))*dt;0 1 x0(4)*cos(x0(3))*dt sin(x0(3))*dt;0 0 1 0;0 0 0 1];
end
% Get CC matrix
function CC = getcc(x0, x1, u0, dt)
 CC = [x0(4)*sin(x0(3))*x0(3)*dt-x0(4)*cos(x0(3))*dt+x1(1)-x0(1);-x0(4)*cos(x0(3))*x0(3)*dt-x0(4)*sin(x0(3))*dt+x1(2)-x0(2);-u0(1)*dt+x1(3)-x0(3);-u0(2)*dt+x1(4)-x0(4)];
end
% Get Ax matrix
function Ax = getax(x_0, dt, N, nx)
x0 = x_0(:,1);
AA = getaa(x0, dt);
Ax = kron(speye(N+1), -speye(nx)) + kron(sparse(diag(ones(N, 1), -1)), AA);
for i = 1 : (N-1)
   x0 = x_0(:,i+1);
   AA = getaa(x0, dt);
   Ax(nx*(i+1)+1:nx*(i+1)+nx,nx*i+1:nx*i+nx) = AA;
end
end
% Get Cx matrix
function Cx = getcx(x_0, u_0, dt, N, nx)
Cx = zeros(N*nx, 1);
for i = 1 : N
    u0 = u_0(:,i);
    x0 = x_0(:,i);
    x1 = x_0(:,i+1);
    CC = getcc(x0, x1, u0, dt);
    Cx((i-1)*nx+1:(i-1)*nx+nx) = CC;
end
end
% Get A2cbf
function A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0)
[~, yy] = size(abc);
Acbfx2 = zeros((yy-2), yy*nx);
Acbfx22 = zeros((yy-2), yy*nx);
Acbfu2 = zeros((yy-2), (yy-1)*nu);
for i = 1 : (yy-2)
    Acbfx2(i, (i*nx)+1) = (gamma1-1/dt)*abc(1, (i+1));
    Acbfx2(i, (i*nx)+2) = (gamma1-1/dt)*abc(2, (i+1));
    Acbfx2(i, ((i+1)*nx)+1) = (1/dt)*abc(1, (i+2));
    Acbfx2(i, ((i+1)*nx)+2) = (1/dt)*abc(2, (i+2));
end
for i = 1 : (yy-2)
    Acbfx22(i, (1*nx)+1) = -(1-dt*gamma2)^(i)/dt*abc(1, (1+1));
    Acbfx22(i, (1*nx)+2) = -(1-dt*gamma2)^(i)/dt*abc(2, (1+1));
end
Acbfx2 = Acbfx2 + Acbfx22;
for i = 1 : (yy-2)
    Acbfu2(i, ((i-1)*nu+4)) = - (1 - dt * gamma2)^(i)*(gamma1-1/dt)*(abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
A2cbf = [Acbfx2 Acbfu2];
end
% Get lcbf2
function lcbf2 = getlcbf2(abc, dt, gamma1, gamma2)
[~, yy] = size(abc);
lcbf2 = zeros((yy-2),1);
for i = 1 : (yy-2)
    lcbf2(i, 1) = -abc(3, (i+2))/dt-(gamma1-1/dt)*abc(3, (i+1))+(1 - dt * gamma2)^(i)*abc(3, (1+1))/dt;
end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function bnet = bayes_update_params(bnet, cases, clamped)
% BAYES_UPDATE_PARAMS Bayesian parameter updating given completely observed data
% bnet = bayes_update_params(bnet, cases, clamped)
%
% If there is a missing data, you must use EM.
% cases(i,m) is the value assigned to node i in case m (this can also be a cell array).
% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros).
% Clamped nodes are not updated.
% If there is a single case, clamped is a list of the clamped nodes, not a bit vector.

%if iscell(cases), usecell = 1; else usecell = 0; end

n = length(bnet.dag);
ncases = size(cases, 2);
if n ~= size(cases, 1)
  error('data must be of size nnodes * ncases');
end

if ncases == 1 % clamped is a list of nodes
  if nargin < 3, clamped = []; end
  clamp_set = clamped;
  clamped = zeros(n,1);
  clamped(clamp_set) = 1;
else % each row of clamped is a bit vector
  if nargin < 3, clamped = zeros(n,ncases); end
end

for i=1:n
  e = bnet.equiv_class(i);
  if adjustable_CPD(bnet.CPD{e})
    u = find(clamped(i,:)==0);
    ps = parents(bnet.dag, i);
    bnet.CPD{e} = bayes_update_params(bnet.CPD{e}, cases(i,u), cases(ps,u));
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function value = v_hofstadter ( n )

%*****************************************************************************80
%
%% V_HOFSTADTER computes the Hofstadter V sequence.
%
%  Discussion:
%
%    V(N) = 0                          if N = 0
%         = 1                          if 1 <= N <= 4
%         = V(N-V(N-1)) + V(N-V(N-4)), otherwise.
%
%    V(N) is defined for all nonnegative integers.
%
%  Table:
%
%     N  V(N)
%    --  ----
%
%     0     0
%     1     1
%     2     1
%     3     1
%     4     1
%     5     2
%     6     3
%     7     4
%     8     5
%     9     5
%    10     6
%    11     6
%    12     7
%    13     8
%    14     8
%    15     9
%    16     9
%    17    10
%    18    11
%    19    11
%    20    11
%    21    12
%    22    12
%    23    13
%    24    14
%    25    14
%    26    15
%    27    15
%    28    16
%    29    17
%    30    17
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the function.
%
%    Output, integer VALUE, the value of the function.
%
  if ( n <= 0 )
    value = 0;
  elseif ( n <= 4 )
    value = 1;
  else
    value = v_hofstadter ( n - v_hofstadter(n-1) ) ...
          + v_hofstadter ( n - v_hofstadter(n-4) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Demo for DCM for fMRI (without Balloon observer)
% This demo inverts the DCM for fMRI model, without the ballon model, which
% is replaced by a nonlinear sigmoid observation function.

% close all
clear variables

%% First get neural level dynamics
n_t = 1.5e2;
TR = 2e0;
microDT = 1e-1;
nreg = 3;
nu = 2;
u       = zeros(2,n_t);         % inputs
u(1,2:15) = 1;
u(1,30:31) = 1;
u(1,50:51) = 1;
u(2,30:31) = 1;
u(2,60:61) = 1;
alpha   = Inf;
sigma   = Inf;
% DCM specification
A = [0 1 1
     1 0 1
     0 1 0];
B{1} = zeros(nreg,nreg);
B{2} = [0 0 0
        1 0 0
        0 0 0];
C = [1 0
     0 0
     0 0];
D{1} = [0 0 0
        0 0 0
        0 1 0];
D{2} = zeros(nreg,nreg);
D{3} = zeros(nreg,nreg);
% Prepare precalculated matrices for DCM inversion
[options.inF] = prepare_dcm(A,B,C,D);
inG.G0 = 50;
inG.k = 1;
% A matrix
t_A = exp([ -1.5
            -1.5
            -0.5
            -2.5
            -1.5 ]);
% self-inhibition gain
t_Aself = -1;
% B matrices
t_B{1} = [];
t_B{2} = exp([ -0.5 ]);
% C matrix
t_C = exp([ +0.1 ]);
% D matrices
t_D{1} = 0*exp([ -2 ]);
t_D{2} = [];
t_D{3} = [];
% fill in theta
theta(options.inF.indA) = t_A;
for i=1:nu
    theta(options.inF.indB{i}) = t_B{i};
end
theta(options.inF.indC) = t_C;
for i=1:nreg
    theta(options.inF.indD{i}) = t_D{i};
end
theta(options.inF.indself) = t_Aself;
theta = theta(:);
phi     = [];
% Build options structure for temporal integration of SDE
options.inG     = inG;
options.decim = max([1,floor(TR./microDT)]);
options.inF.deltat = TR./options.decim;
dim.n_theta         = length(theta);
dim.n_phi           = 0;
dim.n               = nreg;
dim.p               = nreg;
dim.n_t             = n_t;
f_fname = @f_dcm4fmri;
g_fname = @g_Id;
% get micro-time dynamics
out0.options = options;
out0.options.dim = dim;
out0.options.OnLine = 0;
out0.options.microU = 0;
out0.suffStat.dx = [];
out0.options.f_fname = f_fname;
out0.options.g_fname = g_fname;
[xx] = VBA_microTime(...
    struct('muTheta',theta,'muPhi',phi,'muX0',[0;0;0]),u,out0);

%% then get full dynamics (including hemodynamics)
[thetaHRF,phiHRF]   = get_HRFparams(TR,options.inF.deltat);
f_fname2             = @f_DCMwHRF;
g_fname2             = @g_HRF3;
[options]           = prepare_fullDCM(A,B,C,D,TR,options.inF.deltat);
options.inG.TE      = 0.04;
% simulated evolution parameters: hemodynamic level
t_E0 = thetaHRF(1)*ones(nreg,1);       % HbO2 extraction fraction gain
t_tau0 = thetaHRF(2)*ones(nreg,1);     % mean blood transit time gain
t_kaf = thetaHRF(3)*ones(nreg,1);      % vasodilatory signal feedback regulation
t_kas = thetaHRF(4)*ones(nreg,1);      % vasodilatory signal decay gain
t_alpha = thetaHRF(6)*ones(nreg,1);    % vessel stifness gain
% simulated observation parameters
p_E0 = phiHRF(1)*ones(nreg,1);       % HbO2 extraction fraction gain
p_epsilon = phiHRF(2)*ones(nreg,1);  % ratio of intra- and extravascular signal
% fill in theta
theta(options.inF.ind1) = t_E0;
theta(options.inF.ind2) = t_tau0;
theta(options.inF.ind3) = t_kaf;
theta(options.inF.ind4) = t_kas;
theta(options.inF.ind5) = t_alpha;
phi(options.inG.ind1) = p_E0;
phi(options.inG.ind2) = p_epsilon;
phi = phi(:);
dim.n_theta         = options.inF.ind5(end);
dim.n_phi           = 2*nreg;
dim.n               = 5*nreg;
dim.p               = nreg;
dim.n_t             = n_t;
X0 = kron(ones(nreg,1),[0;0;0;0;0]);
% get micro-time dynamics
out.options = options;
out.options.dim = dim;
out.options.OnLine = 0;
out.options.microU = 0;
out.suffStat.dx = [];
out.options.f_fname = f_fname2;
out.options.g_fname = g_fname2;
[xx0,y] = VBA_microTime(...
    struct('muTheta',theta,'muPhi',phi,'muX0',X0),u,out);

%% Finally, convolve neural dynamics with canonical HRF
f_fname3     = @f_HRF2;               % Ballon model evolution function
g_fname3     = @g_HRF3;               % Balloon model observation function
options.inF.fullDCM = 0;
options.inG.fullDCM = 0;
options.inG.TE      = 0.04;
options.microU      = 1;
dim.n_theta         = 6;
dim.n_phi           = 2;
dim.n               = 4;
dim.p               = 1;
dim.n_t             = n_t;
dim.u               = 1;
% get micro-time dynamics
out1.options = options;
out1.options.dim = dim;
out1.dim = dim;
out1.options.OnLine = 0;
out1.options.microU = 1;
out1.suffStat.dx = [];
out1.options.f_fname = f_fname3;
out1.options.g_fname = g_fname3;
yhrf = zeros(size(xx));
xx1 = zeros(size(xx0));
xx1(options.inF.n5,:) = xx;
for i=1:nreg
    u = [xx(i,:)];
    [xxhrf,yhrf(i,:)] = VBA_microTime(...
        struct('muTheta',thetaHRF,'muPhi',phiHRF,'muX0',[0;0;0;0]),u,out1);
    xx1(options.inF.n1(i),:) = xxhrf(1,:);
    xx1(options.inF.n2(i),:) = xxhrf(2,:);
    xx1(options.inF.n3(i),:) = xxhrf(3,:);
    xx1(options.inF.n4(i),:) = xxhrf(4,:);
end

author: MBB-team

repo: VBA-toolbox


text:
function f = flops_det(n)
% FLOPS_DET     Flops for matrix determinant.
% FLOPS_DET(n) returns the number of flops for det(eye(n)).

if n == 1
  f = 1;
else
  % this is from logdet
  f = flops_chol(n) + n;
end

author: Cloud-CV

repo: object-proposals


text:
classdef nnnormalizelp < nntest
  properties (TestParameter)
    h = {1 2 3 4}
    w = {1 2 3 4}
    d = {2 3 4}
    p = {2 4}
  end

  methods (Test)
    function basicl2(test, h,w,d)
      x = test.randn(h,w,d,3,'single') ;
      y = vl_nnnormalizelp(x) ;
      dzdy = test.rand(size(y),'single')-0.5 ;
      dzdx = vl_nnnormalizelp(x,dzdy) ;
      test.der(@(x) vl_nnnormalizelp(x), x, dzdy, dzdx, test.range * 1e-3, 0.3) ;
    end

    function lp(test, p)
      x = test.randn(2,3,5,3,'single') / test.range ;
      y = vl_nnnormalizelp(x, [], 'p', p) ;
      dzdy = test.rand(size(y),'single')-0.5 ;
      dzdx = vl_nnnormalizelp(x,dzdy, 'p', p) ;
      test.der(@(x) vl_nnnormalizelp(x,[],'p',p), x, dzdy, dzdx, 1e-4, 0.3) ;
    end

  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
% FHN evolution function:
%
%                   [fx] = f_CaBBI_FHN(Xt,Theta,Isyn,inF)
%
% this function contains the equations and parameters of the FHN neuron
% model, and the over-time integration scheme

function [fx] = f_CaBBI_FHN(Xt,Theta,Isyn,inF)

dt = inF.dt;
Tau_Ca0 = inF.Tau_Ca0;

% Hidden variables
V = Xt(1);                                                                 % voltage
W = Xt(2);                                                                 % recovery variable
Ca = Xt(3);                                                                % [Ca2+] kinetics

% evolution parameters
Tau_Ca  = Tau_Ca0* exp(Theta(1));                                          % decay time-constant of [Ca2+] kinetics
Ca_base = Theta(2);                                                        % basal value of [Ca2+]
scale   = inF.k_Ca0* exp(Theta(3));                                        % scale parameter of [Ca2+] kinetis

% parameters of the evolution mode
g_Ca  = 5;                                                                 % maximal conductance of the calcium channel
E_Ca  = 120;                                                               % reversal potential of the calcium channel

% compute the current of HVA calcium channel
V_scaled = (30*V) - 40;
s_inf = 1./( 1 + exp(- (V_scaled+45)/5 ) );
I_Ca  = g_Ca .* s_inf.* (V_scaled-E_Ca);                                   % HVA calcium current

% build the right hand side (RHS) of the evolution's differential equations
xdot = [ ( V - (power(V,3)/3) - W  + Isyn )
         (0.08*(V+0.7-(0.8*W)))
         ( -scale*I_Ca - ( (Ca-Ca_base)/Tau_Ca) )
         ];

fx = Xt + dt.*xdot;


text:
function L = image_laplacian(varargin)
  % IMAGE_LAPLACIAN Compute the image laplacian for a given image in the manner
  % of "Colorization using Optimization" by [Levin et al. 2004]. This is a sort
  % of amalgamation of the literal description in the paper but the knowledge
  % that a Laplacian (rather than a bi-Laplacian) is being used.
  %
  % L = IMAGE_LAPLACIAN(im)
  % L = IMAGE_LAPLACIAN(im,'ParameterName',ParameterValue)
  %
  % Inputs:
  %   im  h by w by (3|1) image (expects double)
  %   Optional:
  %     'Omega' followed by an omega value
  % Outputs:
  %   L w*h by w*h sparse laplacian matrix
  %
  % See also: cotmatrix, levin, lischinski
  %

  im = varargin{1};

  % width and height
  h = size(im,1);
  w = size(im,2);
  %   image size
  %  w     h    c  harmonic biharmonic
  %  185   138  1  0.05     0.1
  %  185   138  3           0.1
  %  93    69   3           0.1
  omega = 0.1;
  % number of vertices/pixels
  n = h*w;

  if ~strcmp(class(im),'double')
    warning('Casting input to double: `im = im2double(im)`');
    im = im2double(im);
  end

  ii = 2;
  while(ii <= nargin)
    switch varargin{ii}
    case 'Omega'
      ii = ii + 1;
      assert(ii<=nargin)
      omega = varargin{ii};
    otherwise
      error('Unsupported parameter');
    end
    ii = ii+1;
  end

  % runs across height fastest
  I = (1:n)';
  I = reshape(I,h,w);
  Ileft = I(1:h,1:(w-1));
  Iright = I(1:h,2:w);
  Ibottom = I(1:(h-1),1:w);
  Itop= I(2:h,1:w);
  % determine neighborhood via FD stencil
  E = [ ...
    Ileft(:) Iright(:);...
    Iright(:) Ileft(:);...
    Ibottom(:) Itop(:);...
    Itop(:) Ibottom(:);...
    ];
  assert(size(E,1) == ((h-1)*w + (w-1)*h)*2);
  r = E(:,1);
  s = E(:,2);
  imlong = reshape(im,size(im,1)*size(im,2),size(im,3));

  wrs = exp(-sum((imlong(r,:)-imlong(s,:)).^2,2)./(2*omega.^2));

  L = sparse(r,s,wrs,n,n);
  % We've made L(r,s) = wrs, now we need L(r,r) = ∑L(r,s) and L(r,s) = -wrs
  L = diag(sum(L,2))-L;
end

author: alecjacobson

repo: gptoolbox


text:
function [param, names] = gpsimMapExtractParam(model)

% GPSIMMAPEXTRACTPARAM Extract the parameters of a GPSIMMAP model.
% FORMAT
% DESC extracts the model parameters from a structure containing
% the information about a Gaussian process for single input motif
% modelling.
% ARG model : the model structure containing the information about
% the model.
% RETURN params : a vector of parameters from the model.
%
% SEEALSO : gpsimMapCreate, gpsimMapExpandParam, modelExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFIED : Pei Gao, 2008

% SHEFFIELDML

if nargout>1
  [param, names] = kernExtractParam(model.kern);
else
  param = kernExtractParam(model.kern);
end

fhandle = str2func([model.Transform 'Transform']);

for i = 1:model.numGenes
  if isfield(model,'bTransform') && isempty(model.bTransform);
    param = [param model.B(i)];
  else
    param = [param fhandle(model.B(i), 'xtoa')];
  end

  param = [param fhandle(model.S(i), 'xtoa') fhandle(model.D(i), ...
                                                    'xtoa')];

  if  isfield(model, 'includeRepression') && model.includeRepression
    if isfield(model,'alphaTransform') && isempty(model.alphaTransform);
      param = [param model.alpha(i)];
    else
      param = [param fhandle(model.alpha(i), 'xtoa')];
    end
  end

  if model.ngParam > 0
    param = [param (fhandle(model.gParam(:,i), 'xtoa'))'];
  end

  if nargout >1
    names{end+1} = ['Basal' num2str(i)];
    names{end+1} = ['Sensitivity' num2str(i)];
    names{end+1} = ['Decay' num2str(i)];

    if isfield(model, 'isGroupNonlinearity') && model.isGroupNonlinearity
      if strcmp(model.nonLinearity{i}, 'repression')
        names{end+1} = ['alpha' num2str(i)];
      end
    else
      names{end+1} = ['alpha' num2str(i)];
    end

    if model.ngParam > 0
      ngParamk = model.ngParam/model.numGenes;
      names{(end+1):(end+ngParamk)} = ['Parameters for g' num2str(i)];
    end
  end
end

if isfield(model,'includeNoise') && model.includeNoise
  param = [param fhandle(sqrt(model.noiseVar), 'xtoa')];
  if nargout > 1
    for i=1:model.numGenes
      names = {names{:} ['noiseSd' num2str(i)]};
    end
  end
end

if isfield(model, 'fix')
  for i = 1:length(model.fix)
    param(model.fix(i).index) = model.fix(i).value;
  end
end

param = real(param);

% Check if there is a mean function.
if isfield(model, 'meanFunction') & ~isempty(model.meanFunction)
  if nargout>1
    [meanFuncParams, meanFuncNames] = modelExtractParam(model.meanFunction);
  else
    meanFuncParams = modelExtractParam(model.meanFunction);
  end
else
  meanFuncParams =[];
  meanFuncNames = {};
end

param = [param meanFuncParams];
if nargout > 1
  names = {names{:} meanFuncNames{:}};
end


text:
function test_failed=test_firwin
%TEST_FIRWIN  Test the firwin windows
%
%  This test script verifies the properties listed in the help of firwin

allwins = getfield(arg_firwin,'flags','wintype');

test_failed=0;

disp(' ===============  TEST_FIRWIN ================');

for L=[18,19,20,21]

  for cent=0:1
    if cent==0
      centtype='wp';
    else
      centtype='hp';
    end;
    for ii=1:length(allwins);
      winname=allwins{ii};

      [g,info]=firwin(winname,L,centtype);

      res = 1-isevenfunction(fir2long(g,2*L),centtype);
      [test_failed,fail]=ltfatdiditfail(res,test_failed);

      s=sprintf(['SYMM %10s %s L: %i %0.5g %s'],winname,centtype,L,res,fail);
      disp(s);

      if cent==0
        res=1-g(1);

        [test_failed,fail]=ltfatdiditfail(res,test_failed);

        s=sprintf(['PEAK %10s %s L: %i %0.5g %s'],winname,centtype,L,res,fail);
        disp(s);

      end;

      if mod(L,2)==0
        if info.ispu
          gpu=g+fftshift(g);
          res=norm(gpu-gpu(1)*ones(L,1));

          [test_failed,fail]=ltfatdiditfail(res,test_failed);

          s=sprintf(['PU   %10s %s L: %i %0.5g %s'],winname,centtype,L,res,fail);
          disp(s);

        end;

        if info.issqpu
          gpu=g.^2+fftshift(g.^2);
          res=norm(gpu-gpu(1)*ones(L,1));

          [test_failed,fail]=ltfatdiditfail(res,test_failed);

          s=sprintf(['SQPU %10s %s L: %i %0.5g %s'],winname,centtype,L,res,fail);
          disp(s);

        end;
      end;
    end;

  end;

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ExportFile, sFileOut] = export_data(DataFile, ChannelMat, ExportFile, FileFormat)
% EXPORT_DATA: Exports a recordings file to one of the supported file formats.
%
% USAGE:  [ExportFile, sFileOut] = export_data( DataFile,        [], ExportFile=[ask], FileFormat=[detect] )
%         [ExportFile, sFileOut] = export_data( DataMat, ChannelMat, ExportFile=[ask], FileFormat=[detect] ) : Save a data structure
%                                  export_data( DataFiles{}, ... )                      : Batch process over multiple files
%
% INPUT:
%     - DataFile   : Brainstorm data file name to be exported
%     - DataMat    : Brainstorm data structure to be exported
%     - ExportFile : Full path to target file (extension will determine the format)
%                    If not specified: asked to the user

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

% ===== PARSE INPUTS =====
if (nargin < 4) || isempty(FileFormat)
    FileFormat = [];
end
if (nargin < 3) || isempty(ExportFile)
    ExportFile = [];
end
if (nargin < 2) || isempty(ChannelMat)
    ChannelMat = [];
end
% Initialize returned variables
sFileOut = [];

% ===== RECURSIVE CALLS =====
% CALL: export_data( DataFiles, ... )
if iscell(DataFile)
    % Single file
    if (length(DataFile) == 1)
        DataFile = DataFile{1};
    % Multiple files
    else
        AllOutputs = cell(1,length(DataFile));
        % Call function once to get the output path
        [AllOutputs{1}, sFileOut1] = export_data(DataFile{1}, ChannelMat, ExportFile);
        if isempty(AllOutputs{1})
            ExportFile = [];
            return;
        end
        % Get output path
        [outPath, outBase, outExt] = bst_fileparts(AllOutputs{1});
        % Loop on the other files
        for i = 2:length(DataFile)
            % Build output file name
            [dataPath, dataBase, dataExt] = bst_fileparts(DataFile{i});
            newFile = bst_fullfile(outPath, [dataBase, outExt]);
            newFile = strrep(newFile, '_data', '');
            newFile = strrep(newFile, 'data_', '');
            newFile = strrep(newFile, '0raw_', '');
            AllOutputs{i} = newFile;
            % Export file
            export_data(DataFile{i}, ChannelMat, AllOutputs{i}, sFileOut1.format);
        end
        ExportFile = AllOutputs;
        return;
    end
end

% ===== LOAD DATA =====
% CALL: export_data( DataFile, [], ExportFile )
if ischar(DataFile)
    isRawIn = ~isempty(strfind(DataFile, '_0raw'));
    % Load initial file
    DataMat = in_bst_data(DataFile);
    % Get raw file structure
    if isRawIn
        sFileIn = DataMat.F;
        if (length(sFileIn.epochs) > 1)
            error('Cannot export epoched files.');
        end
    else
        sFileIn = in_fopen(DataFile, 'BST-DATA');
    end
    % Load channel file
    ChannelFile = bst_get('ChannelFileForStudy', DataFile);
    ChannelMat = in_bst_channel(ChannelFile);
% CALL: export_data( DataMat, ChannelMat, ExportFile )
elseif isstruct(DataFile)
    DataMat = DataFile;
    DataFile = [];
    isRawIn = 0;
    % Create a fake sFile structure
    sFileIn = in_fopen(DataMat, 'BST-DATA');
else
    error('Invalid call.');
end

% ===== SELECT OUTPUT FILE =====
if isempty(ExportFile)
    % === Build a default filename ===
    % Get default directories and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    % Get default extension
    switch (DefaultFormats.DataOut)
        case {'', 'BST-BIN'},   DefaultExt = '.bst';
        case 'FT-TIMELOCK',     DefaultExt = '.mat';
        case 'SPM-DAT',         DefaultExt = '.mat';
        case 'EEG-BRAINAMP',    DefaultExt = '.eeg';
        case 'EEG-CARTOOL-EPH', DefaultExt = '.eph';
        case 'EEG-EGI-RAW',     DefaultExt = '.raw';
        case 'EEG-EDF',         DefaultExt = '.edf';
        case 'NIRS-SNIRF',      DefaultExt = '.snirf';
        case 'ASCII-CSV',       DefaultExt = '.csv';
        case 'ASCII-CSV-HDR',   DefaultExt = '.csv';
        case 'ASCII-CSV-HDR-TR',DefaultExt = '.csv';
        case 'ASCII-TSV',       DefaultExt = '.tsv';
        case 'ASCII-TSV-HDR',   DefaultExt = '.tsv';
        case 'ASCII-TSV-HDR-TR',DefaultExt = '.tsv';
        case 'ASCII-SPC',       DefaultExt = '.txt';
        case 'ASCII-SPC-HDR',   DefaultExt = '.txt';
        case 'EXCEL',           DefaultExt = '.xlsx';
        case 'EXCEL-TR',        DefaultExt = '.xlsx';
        case 'BST',             DefaultExt = '_timeseries.mat';
        otherwise,              DefaultExt = '_timeseries.mat';
    end
    % Build default output filename
    if ~isempty(DataFile)
        [BstPath, BstBase, BstExt] = bst_fileparts(DataFile);
    else
        BstBase = file_standardize(DataMat.Comment);
    end
    DefaultExportFile = bst_fullfile(LastUsedDirs.ExportData, [BstBase, DefaultExt]);
    DefaultExportFile = strrep(DefaultExportFile, '_data', '');
    DefaultExportFile = strrep(DefaultExportFile, 'data_', '');
    DefaultExportFile = strrep(DefaultExportFile, '0raw_', '');

    % === Ask user filename ===
    % RAW file or imported
    if isRawIn
        FileFilters = bst_get('FileFilters', 'rawout');
    else
        FileFilters = bst_get('FileFilters', 'dataout');
    end
    % Put file
    [ExportFile, FileFormat, FileFilter] = java_getfile( 'save', ...
        'Export MEG/EEG recordings...', ... % Window title
        DefaultExportFile, ...              % Default directory
        'single', 'files', ...              % Selection mode
        FileFilters, ...
        DefaultFormats.DataOut);
    % If no file was selected: exit
    if isempty(ExportFile)
        return;
    end
    % Save new default export path
    LastUsedDirs.ExportData = bst_fileparts(ExportFile);
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default export format
    DefaultFormats.DataOut = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);

% Guess file format based on its extension
elseif isempty(FileFormat)
    [BstPath, BstBase, BstExt] = bst_fileparts(ExportFile);
    switch lower(BstExt)
        case '.bst',   FileFormat = 'BST-BIN';
        case '.eeg',   FileFormat = 'EEG-BRAINAMP';
        case '.eph',   FileFormat = 'EEG-CARTOOL-EPH';
        case '.raw',   FileFormat = 'EEG-EGI-RAW';
        case '.edf',   FileFormat = 'EEG-EDF';
        case '.snirf', FileFormat = 'NIRS-SNIRF';
        case '.txt',   FileFormat = 'ASCII-SPC';
        case '.csv',   FileFormat = 'ASCII-CSV-HDR';
        case '.tsv',   FileFormat = 'ASCII-TSV-HDR';
        case '.xlsx',  FileFormat = 'EXCEL';
        case '.mat',   FileFormat = 'BST';
        otherwise,     error('Unsupported file extension.');
    end
end
% Show progress bar
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'Export EEG/MEG recordings', 'Exporting file...');
end
% Option for input raw file
if isRawIn
    ImportOptions = db_template('ImportOptions');
    ImportOptions.ImportMode      = 'Time';
    ImportOptions.DisplayMessages = 0;
    ImportOptions.UseCtfComp      = 0;
    ImportOptions.UseSsp          = 0;
    ImportOptions.RemoveBaseline  = 'no';
end

% ===== REMOVE ANNOTATION CHANNELS =====
% Detect annotation channels (exclude BST-BIN output, because we want to keep everything)
if ~strcmpi(FileFormat, 'BST-BIN')
    iAnnot = channel_find(ChannelMat.Channel, {'EDF', 'BDF', 'KDF'});
else
    iAnnot = [];
end
% List of input and output channels
iChannelsIn = setdiff(1:length(ChannelMat.Channel), iAnnot);
iChannelsOut = 1:length(iChannelsIn);
% Selected channel mats
ChannelMatIn = ChannelMat;
ChannelMatOut = ChannelMat;
% Remove unwanted channels from channel file
if ~isempty(iAnnot)
    ChannelMatOut.Channel = ChannelMatOut.Channel(iChannelsIn);
    for iProj = 1:length(ChannelMatOut.Projector)
        if isequal(ChannelMatOut.Projector(iProj).SingVal, 'REF')
            ChannelMatOut.Projector(iProj).Components = ChannelMatOut.Projector(iProj).Components(iChannelsIn, iChannelsIn);
        else
            ChannelMatOut.Projector(iProj).Components = ChannelMatOut.Projector(iProj).Components(iChannelsIn, :);
        end
    end
end
% Remove unwanted channels from input data file
if ~isRawIn
    if isfield(DataMat, 'F') && ~isempty(DataMat.F)
        DataMat.F = DataMat.F(iChannelsIn,:);
    end
    if isfield(DataMat, 'ChannelFlag') && ~isempty(DataMat.ChannelFlag)
        DataMat.ChannelFlag = DataMat.ChannelFlag(iChannelsIn);
    end
end

% ===== CREATE OUTPUT RAW FILE =====
% Output data as raw file (continuous writers routines)
isRawOut = ismember(FileFormat, {'BST-BIN', 'EEG-EGI-RAW', 'SPM-DAT', 'EEG-EDF', 'EEG-BRAINAMP'});
% Open output file
if isRawOut
    [sFileOut, errMsg] = out_fopen(ExportFile, FileFormat, sFileIn, ChannelMatOut, iChannelsIn);
    % Error management
    if isempty(sFileOut) && ~isempty(errMsg)
        error(errMsg);
    elseif ~isempty(errMsg)
        disp(['BST> Warning: ' errMsg]);
    end
end

% ===== RAW IN / RAW OUT =====
if isRawIn && isRawOut
    % Check that in!=out
    if file_compare(sFileIn.filename, sFileOut.filename)
        error('Input and output files are the same.');
    end
    % Get default epoch size
    EpochSize = bst_process('GetDefaultEpochSize', sFileOut);
    % Process by sample blocks
    nSamples = round((sFileOut.prop.times(2) - sFileOut.prop.times(1)) * sFileOut.prop.sfreq) + 1;
    nBlocks = ceil(nSamples / EpochSize);
    % Show progress bar
    if ~isProgress
        bst_progress('start', 'Export EEG/MEG recordings', 'Exporting file...', 0, nBlocks);
    end
    % Copy files by block
    for iBlock = 1:nBlocks
        % Get sample indices
        SamplesBounds = sFileOut.prop.times(1) * sFileOut.prop.sfreq + [(iBlock-1) * EpochSize, min(iBlock*EpochSize-1, nSamples-1)];
        % Read from input file
        F = in_fread(sFileIn, ChannelMatIn, 1, SamplesBounds, iChannelsIn, ImportOptions);
        % Save to output file
        sFileOut = out_fwrite(sFileOut, ChannelMatOut, 1, SamplesBounds, iChannelsOut, F);
        % Increase progress bar
        if ~isProgress
            bst_progress('inc', 1);
        end
    end

% ===== SAVE FULL FILES =====
else
    % Load full file
    if isRawIn
        F = in_fread(sFileIn, ChannelMatIn, 1, [], iChannelsIn, ImportOptions);
    else
        if isfield(DataMat, 'F') && ~isempty(DataMat.F)
            F = DataMat.F;
        elseif isfield(DataMat, 'ImageGridAmp') && ~isempty(DataMat.ImageGridAmp)
            F = DataMat.ImageGridAmp;
        else
            error('No relevant data to save found in structure.');
        end
    end

    % Save full file
    if isRawOut
        out_fwrite(sFileOut, ChannelMatOut, 1, [], iChannelsOut, F);
    else
        % Switch between file formats
        switch FileFormat
            case 'BST'
                DataMat.F = F;
                bst_save(ExportFile, DataMat, 'v6');
            case 'FT-TIMELOCK'
                DataMat.F = F;
                ftData = out_fieldtrip_data(DataMat, ChannelMatOut, [], 1);
                bst_save(ExportFile, ftData, 'v6');
            case 'EEG-CARTOOL-EPH'
                % Get sampling rate
                samplingFreq = round(1/(DataMat.Time(2) - DataMat.Time(1)));
                % Write header : nb_electrodes, nb_time, sampling_freq
                dlmwrite(ExportFile, [size(F,1), size(F,2), samplingFreq], 'newline', 'unix', 'precision', '%d', 'delimiter', ' ');
                % Write data
                dlmwrite(ExportFile, F' * 1000, 'newline', 'unix', 'precision', '%0.7f', 'delimiter', '\t', '-append');
            case 'NIRS-SNIRF'
                if isRawIn
                    DataMat.Events = DataMat.F.events;
                end
                DataMat.F = F;
                out_data_snirf(ExportFile, DataMat, ChannelMatOut);
            case {'ASCII-SPC', 'ASCII-CSV', 'ASCII-TSV', 'ASCII-SPC-HDR', 'ASCII-CSV-HDR', 'ASCII-TSV-HDR', 'ASCII-CSV-HDR-TR', 'ASCII-TSV-HDR-TR', 'EXCEL', 'EXCEL-TR'}
                out_matrix_ascii(ExportFile, F, FileFormat, {ChannelMatOut.Channel.Name}, DataMat.Time, []);
            otherwise
                error('Unsupported format.');
        end
    end
end

% Hide progress bar
if ~isProgress
    bst_progress('stop');
end

end

author: brainstorm-tools


text:
function a = ijfact1 ( n )

%*****************************************************************************80
%
%% IJFACT1 returns the IJFACT1 matrix.
%
%  Formula:
%
%    A(I,J) = (I+J)!
%
%  Example:
%
%    N = 4
%
%     2   6   24   120
%     6  24  120   720
%    24 120  720  5040
%   120 720 5040 40320
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is a Hankel matrix: constant along anti-diagonals.
%
%    A is integral: int ( A ) = A.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    MJC Gover,
%    The explicit inverse of factorial Hankel matrices,
%    Department of Mathematics, University of Bradford, 1993.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  fact = 1;

  for k = 2 : 2 * n

    fact = fact * k;
    ilo = max ( 1, k - n );
    ihi = min ( n, k - 1 );

    for i = ilo : ihi
      a(i,k-i) = fact;
    end

  end

  return
end


text:
function [ x, w ] = moment_method ( n, moment )

%*****************************************************************************80
%
%% MOMENT_METHOD computes a quadrature rule by the method of moments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, John Welsch,
%    Calculation of Gaussian Quadrature Rules,
%    Mathematics of Computation,
%    Volume 23, Number 106, April 1969, pages 221-230.
%
%  Parameters:
%
%    Input, integer N, the order of the quadrature rule.
%
%    Input, real MOMENT(2*N+1), moments 0 through 2*N.
%
%    Output, real X(N,1), W(N,1), the points and weights of the quadrature rule.
%
  debug = 0;

  if ( debug )
    r8vec_print ( 2 * n + 1, moment, '  Moments:' );
  end
%
%  Define the N+1 by N+1 Hankel matrix H(I,J) = moment(I+J).
%
  h = zeros ( n + 1, n + 1 );
  for i = 1 : n + 1
    for j = 1 : n + 1
      h(i,j) = moment(i+j-1);
    end
  end

  if ( debug )
    r8mat_print ( n + 1, n + 1, h, '  Hankel matrix:' )
  end
%
%  Compute R, the upper triangular Cholesky factor of H.
%
  r = chol ( h );

  if ( debug )
    r8mat_print ( n + 1, n + 1, r, '  Cholesky factor:' )
  end
%
%  Compute ALPHA and BETA from R, using Golub and Welsch's formula.
%
  alpha = zeros ( n, 1 );
  alpha(1) = r(1,2) / r(1,1);
  for i = 2 : n
    alpha(i) = r(i,i+1) / r(i,i) - r(i-1,i) / r(i-1,i-1);
  end

  beta = zeros ( n - 1, 1 );
  for i = 1 : n - 1
    beta(i) = r(i+1,i+1) / r(i,i);
  end
%
%  Compute the points and weights from the moments.
%
  jacobi = diag ( alpha, 0 ) + diag ( beta, -1 ) + diag ( beta, +1 );

  if ( debug )
    r8mat_print ( n, n, jacobi, '  Jacobi matrix: ' );
  end
%
%  Get the eigendecomposition of the Jacobi matrix.
%
  [ evec, eval ] = eig ( jacobi );

  x = diag ( eval );
  w = moment(1) * evec(1,:).^2;
%
%  Destroy all row vectors.
%
  x = x(:);
  w = w(:);

  return
end