text:
classdef linear_svm
% This file defines l2-regularized SVM problem
%
% Inputs:
%       x_train     train data matrix of x of size dxn.
%       y_train     train data vector of y of size 1xn.
%       x_test      test data matrix of x of size dxn.
%       y_test      test data vector of y of size 1xn.
%       lambda      l2-regularized parameter.
% Output:
%       Problem     problem instance.
%
%
% The problem of interest is defined as
%
%           min f(w) = 1/n * sum_i^n f_i(w),
%           where
%           f_i(w) = 1/2 * (max(0.0, 1 - y_i .* (w'*x_i) )^2 + lambda/2 * w^2.
%
% "w" is the model parameter of size d vector.
%
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Feb. 17, 2016
% Modified by H.Kasai on Mar. 23, 2018

    properties
        name;
        dim;
        samples;
        lambda;
        classes;
        hessain_w_independent;
        d;
        n_train;
        n_test;
        x_train;
        y_train;
        x_test;
        y_test;
        x_norm;
        x;
    end

    methods
        function obj = linear_svm(x_train, y_train, x_test, y_test, varargin)

            obj.x_train = x_train;
            obj.y_train = y_train;
            obj.x_test = x_test;
            obj.y_test = y_test;

            if nargin < 5
                obj.lambda = 0.1;
            else
                obj.lambda = varargin{1};
            end

            obj.d = size(x_train, 1);
            obj.n_train = length(y_train);
            obj.n_test = length(y_test);

            obj.name = 'linear_svm';
            obj.dim = obj.d;
            obj.samples = obj.n_train;
            obj.classes = 2;
            obj.hessain_w_independent = false;
            obj.x_norm = sum(obj.x_train.^2,1);
            obj.x = obj.x_train;
        end

        function f = cost(obj, w)

            f_sum = 0.5 * sum(max(0.0, 1 - obj.y_train' .*(w'*obj.x_train)').^2);
            f = f_sum/obj.n_train + obj.lambda/2 * w(:)'*w(:);

        end

        function g = grad(obj, w, indices)

            alpha = w' * obj.x_train(:,indices);
            flag = obj.y_train(indices) .* alpha;
            flag(flag<1.0) = 1;
            flag(flag>1.0) = 0;

            coeff = flag' .* (1 - obj.y_train(indices)' .* alpha') .* obj.y_train(indices)';
            coeff = coeff';
            coeff = repmat(coeff,[obj.d 1]);
            g = obj.lambda * w - sum(coeff .* obj.x_train(:,indices),2)/length(indices);

        end

        function g = full_grad(obj, w)

            g = obj.grad(w, 1:obj.n_train);

        end

        function h = hess(obj, w, indices)

            alpha = w' * obj.x_train(:,indices);
            flag = obj.y_train(indices) .* alpha;
            flag_indices = flag<1.0;

            x_part_new = obj.x_train(:,indices);
            x_part_new = x_part_new(:,flag_indices);

            h = obj.lambda * eye(obj.d) + x_part_new * x_part_new'/length(indices);
        end

        function hv = hess_vec(obj, w, v, indices)

            alpha = w' * obj.x_train(:,indices);
            flag = obj.y_train(indices) .* alpha;
            flag_indices = flag<1.0;

            x_part_new = obj.x_train(:,indices);
            x_part_new = x_part_new(:,flag_indices);

            hv = obj.lambda * v + (x_part_new * (x_part_new'*v) )/length(indices);
        end

        function p = prediction(obj, w)

            p = w' * obj.x_test;

            class1_idx = p>0;
            class2_idx = p<=0;
            p(class1_idx) = 1;
            p(class2_idx) = -1;

        end

        function a = accuracy(obj, y_pred)

            a = sum(y_pred == obj.y_test) / obj.n_test;

        end

        function w_opt = calc_solution(obj, maxiter, method)

            if nargin < 3
                method = 'lbfgs';
            end

            options.max_iter = maxiter;
            options.verbose = true;
            options.tol_optgap = 1.0e-24;
            options.tol_gnorm = 1.0e-16;
            options.step_alg = 'backtracking';

            if strcmp(method, 'sd')
                [w_opt,~] = sd(obj, options);
            elseif strcmp(method, 'cg')
                [w_opt,~] = ncg(obj, options);
            elseif strcmp(method, 'newton')
                options.sub_mode = 'INEXACT';
                options.step_alg = 'non-backtracking';
                [w_opt,~] = newton(obj, options);
            else
                options.step_alg = 'backtracking';
                options.mem_size = 5;
                [w_opt,~] = lbfgs(obj, options);
            end
        end
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
%MDL_FANUC10L  Create kinematic model of Fanuc AM120iB/10L robot
%
% MDL_FANUC10L is a script that creates the workspace variable R which
% describes the kinematic characteristics of a Fanuc AM120iB/10L robot
% using standard DH conventions.
%
% Also defines the workspace vector:
%   q0   mastering position.
%
% Notes::
% - SI units of metres are used.
%
% Author::
%  Wynand Swart,
%  Mega Robots CC, P/O Box 8412, Pretoria, 0001, South Africa,
%  wynand.swart@gmail.com
%
% See also mdl_irb140, mdl_m16, mdl_motomanHP6, mdl_puma560, SerialLink.

% MODEL: Fanuc, AM120iB/10L, 6DOF, standard_DH

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

%Cell: 073-1555-430
%30 Sep 2007
%Fanuc AM120iB/10L robot

%            theta    d      a    alpha
clear L
L(1) = Link([0        0     0.15  -pi/2   0  ]);
L(2) = Link([0        0     0.77   pi     0  ]);
L(3) = Link([0        0     0.1   -pi/2   0  ]);
L(4) = Link([0       -0.96  0      pi/2   0  ]);
L(5) = Link([0        0     0     -pi/2   0  ]);
L(6) = Link([0       -0.1   0      0      0  ]);
%##########################################################
%Pose 0; At MASTERING position;
%##########################################################
q0 =[0   -pi/2   0   0   0   0];
R=SerialLink(L, 'name', 'Fanuc AM120iB/10L');
%##########################################################

author: petercorke

repo: robotics-toolbox-matlab


text:
suite = matlab.unittest.TestSuite.fromFolder(pwd);
suite.run;

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function varargout=ndstest(TOL)
%Performs numerous tests of ndSparse math operations,
%
%  ndstest(TOL)
%
%TOL is a tolerance value on the percent error. Execution will pause in debug
%mode for inspection if any one of the tests exhibits an error greater than
%TOL.

if nargin<1
 TOL=inf; %default tolerance value on discrepancies
end

CHECKTYPES=false;

%%function for measuring error

 err=@(x,y) DiscrepancyMeasure(x,y,TOL,CHECKTYPES);

Pf=srand(3,3,2,4);
Qf=srand(size(Pf));
Af=srand(3,2)*1i;
Bf=srand(3,2);
Cf=srand(3,2);
Sf=srand(3);
Vrf=1:(numel(Pf)-3);
Vcf=Vrf.';
three=single(3);

%%Representations of the above as ndSparse objects
P=ndSparse(Pf(:),size(Pf));   % =Pf
Q=ndSparse(Qf(:),size(Qf));   % =Qf
A=ndSparse(Af);   % =Af
B=ndSparse(Bf);   % =Bf
C=ndSparse(Cf);   % =Cf
S=ndSparse(Sf);   % =Sf
Vc=ndSparse(Vcf);   % =Vcf
Vr=ndSparse(Vrf);   % =Vrf

%%%%%%%%%%%%%%%%%%%%%%%%%%%TESTS%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Test of full()
Error(1)=     err( full(P), Pf  );
Error(end+1)= err( full(Q), Qf  );
Error(end+1)= err( full(S), Sf );

%%Test of sparse2d, sparse

  %Error(end+1)= ~isequal( sparse(A), A ); %obsolete test of sparse(ndSparse)

Error(end+1)= ~isequal( sparse2d(A), reshape(sparse(Af),[],size(Af,ndims(Af)) ) );
Error(end+1)= ~isequal( sparse(A), sparse(Af) );

   %%No sense in proceding if the tests so far didn't pass - the error
   %%calculations rely on the functionality of full@ndSparse()
    if max(Error)>TOL,
        Error,
        error 'Something wrong with sparse() and full() methods';
    end

%Test of logical, double

Error(end+1)=err( class(logical(P)), class(logical(Pf)));
Error(end+1)=err( class(double(P)), class(double(Pf)));

%%Test isnumeric, islogical, isempty, issparse, isfloat, isreal,isinf,
%      isnan, isfinite, isequal, isequalwithequalnans
Error(end+1)= err( isnumeric(P),1);
Error(end+1)= err( ~isnumeric(P>0),1);
Error(end+1)= err( isfloat(P),1);
Error(end+1)= err( ~isfloat(P>0),1);
Error(end+1)= err( ~islogical(P),1);
Error(end+1)= err( islogical(P>0),1);
Error(end+1)= err( ~isempty(P),1);
Error(end+1)= err( isempty(ndSparse([])),1);
Error(end+1)= err( issparse(P),1);
Error(end+1)= err( issparse(A), 1 );
Error(end+1)= err( ~isreal(A), 1 );
Error(end+1)= err( isreal(B), 1 );

     Z=P; Z(1)=nan; Z(2)=inf;
     Zf=full(Z);

Error(end+1)= err( isnan(Z), isnan(Zf) );
Error(end+1)= err( isinf(Z), isinf(Zf) );
Error(end+1)= err( isfinite(Z), isfinite(Zf) );
Error(end+1)= err( ~isequal(Z,Zf), 1 );
Error(end+1)= err( ~isequal(Zf,Z), 1 );
Error(end+1)= err( ~isequal(Z,Zf,Z), 1 );
Error(end+1)= err( ~isequal(Zf,Z,Z), 1 );
Error(end+1)= err( isequalwithequalnans(Z,Zf), 1 );
Error(end+1)= err( isequalwithequalnans(Zf,Z), 1 );
Error(end+1)= err( isequalwithequalnans(Z,Zf,Z), 1 );
Error(end+1)= err( isequalwithequalnans(Zf,Z,Zf), 1 );

%Test of real, imag, conj, abs,sqrt

     Z=ndSparse(Bf+Af); %complex result
     Zf=full(Z);

Error(end+1)= err( real(Z), real(Zf) );
Error(end+1)= err( imag(Z), imag(Zf) );
Error(end+1)= err( conj(Z), conj(Zf) );
Error(end+1)= err( abs(Z), abs(Zf) );
Error(end+1)= err( sqrt(Z), sqrt(Zf) );

%%Test of size

Error(end+1)=err( size(P), size(Pf) );
Error(end+1)=err( size(P,2), size(Pf,2) );

       [mm,nn]  =size(Q);
       [mmm,nnn]=size(Qf);

Error(end+1)=err( [mm,nn] , [mmm,nnn] );

%%Test of reshape

Error(end+1)= err( reshape(P,size(P,1),[]) , reshape(Pf,size(Pf,1),[]));

%%Test of permute, ipermute

       ord=randperm(ndims(Pf));
       Z=permute(P,ord);
       Zf=permute(Pf,ord);

Error(end+1)= err( Z , Zf );
Error(end+1)= err( ipermute(Z,ord) ,  Pf );

%%Test of transpose, ctranspose
Error(end+1)= err( A.' ,  Af.'   );
Error(end+1)= err( A'  ,  Af'   );

%%Test of uplus, uminus

Error(end+1)= err( +P  , +Pf );
Error(end+1)= err( -Q  , -Qf );

%%Test of plus, minus

Error(end+1)= err( P+Q  , Pf+Qf );
Error(end+1)= err( P-Q ,  Pf-Qf );

%%Test of inv

Error(end+1)= err( inv(S)  , inv(Sf) );

%%Test of find

            [II,JJ,KK]=find(P);
            [IIf,JJf,KKf]=find(Pf);

Error(end+1)= err( II  , IIf );
Error(end+1)= err( JJ  , JJf );
Error(end+1)= err( KK  , KKf );

%%Test of mtimes

Error(end+1)= err( A*three , Af*three  );  %scalar with ndSparse
Error(end+1)= err( three*P , three*Pf  );

Error(end+1)= err( S*A , Sf*Af  ); % 2 ndSparses

   x=Af(end,:).';
   y=Af(:,end).';

Error(end+1)= err( A*[x,x] , Af*[x,x]  ); %pre-mult with columnized data
Error(end+1)= err( [y;y]*A , [y;y]*Af  ); %post-mult with columnized data

Error(end+1)= err( S*Af , Sf*Af  ); %mixed op

%%Test of mldivide
Error(end+1)= err( three\P , three\Pf  );  %scalar with ndSparse
Error(end+1)= err( three\A , three\Af  );
Error(end+1)= err( B\A , Bf\Af  ); % 2 ndSparses
Error(end+1)= err( B\Af , Bf\Af  );

%%Test of mrdivide

    Bt=B.'; Btsp=Bf.';

Error(end+1)= err( Bt/three , Btsp/three  );%scalar with ndSparse

Error(end+1)= err( P/three , Pf/three  );%scalar with ndSparse
Error(end+1)= err( A.'/Bt , Af.'/Btsp ); % 2 ndSparses
Error(end+1)= err( Af.'/Bt , Af.'/Btsp );

%%Test of times

Error(end+1)= err( P.*three , Pf.*three  );  %scalar with ndSparse
Error(end+1)= err( three.*P , three.*Pf  );

Error(end+1)= err( P.*Q , Pf.*Qf  ); % 2 ndSparses
Error(end+1)= err( Pf.*Q , Pf.*Qf  );

%%Test of rdivide

Error(end+1)= err( Q./three , Qf./three  );  %scalar with ndSparse
Error(end+1)= err( three./Q , three./Qf  );

Error(end+1)= err( P./Q , Pf./Qf  ); %2 ndSparses
Error(end+1)= err( Q./P , Qf./Pf  ); %2 ndSparses
Error(end+1)= err( Qf./P , Qf./Pf  ); %mixed

%%Test of ldivide

Error(end+1)= err( Q.\three , Qf.\three  );  %scalar with ndSparse
Error(end+1)= err( three.\Q , three.\Qf  );

Error(end+1)= err( Q.\P , Qf.\Pf  ); %2 ndSparses
Error(end+1)= err( P.\Q , Pf.\Qf  ); %2 ndSparses
Error(end+1)= err( Pf.\Q , Pf.\Qf  );

%%Test of power, mpower
Error(end+1)= err( P.^three ,  Pf.^three   );
Error(end+1)= err( three.^P ,  three.^Pf   );
Error(end+1)= err( Q.^P ,  Qf.^Pf   );
Error(end+1)= err( Qf.^P ,  Qf.^Pf   );

%%Test of mpower
Error(end+1)= err( S^three  ,  Sf^three   );

%%Test of relops

Error(end+1)= err( P>P ,  Pf>Pf   );
Error(end+1)= err( P>Pf ,  Pf>Pf   );
Error(end+1)= err( P>three/6  , Pf>three/6   );

Error(end+1)= err( P>=P ,  Pf>=Pf   );
Error(end+1)= err( P>=Pf ,  Pf>=Pf   );
Error(end+1)= err( P>=three/6  , Pf>=three/6   );

Error(end+1)= err( P<P ,  Pf<Pf   );
Error(end+1)= err( P<Pf ,  Pf<Pf   );
Error(end+1)= err( P<three/6  , Pf<three/6   );

Error(end+1)= err( P<=P ,  Pf<=Pf   );
Error(end+1)= err( P<=Pf ,  Pf<=Pf   );
Error(end+1)= err( P<=three/6  , Pf<=three/6   );

Error(end+1)= err( P==P ,  Pf==Pf   );
Error(end+1)= err( P==Pf ,  Pf==Pf   );
Error(end+1)= err( P==three/6  , Pf==three/6   );

Error(end+1)= err( P~=P ,  Pf~=Pf   );
Error(end+1)= err( P~=Pf ,  Pf~=Pf   );
Error(end+1)= err( P~=three/6  , Pf~=three/6   );

%%Test of logical ops

Error(end+1)= err( P&P ,  Pf&Pf   );
Error(end+1)= err( P&Pf ,  Pf&Pf   );
Error(end+1)= err( P|three/6  , Pf|three/6   );

Error(end+1)= err( ~P ,   ~Pf   );
Error(end+1)= err( ~Q  ,  ~Qf  );

%%Test of sum

Error(end+1)= err( sum(P,1)  , sum(Pf,1) );
Error(end+1)= err( sum(Q,2)  , sum(Qf,2) );
Error(end+1)= err( sum(P)    , sum(Pf) );
Error(end+1)= err( sum(P,'native')    , sum(Pf,'native') );
Error(end+1)= err( sum(P,3,'double')    , sum(Pf,3,'double') );
Error(end+1)= err( sum(P,4)    , sum(Pf,4) );
Error(end+1)= err( sum(P,5)    , sum(Pf,5) );

%%Test of cat, horzcat, vertcat

Error(end+1)= err([P,P] ,[Pf,Pf]);
Error(end+1)= err([Pf,P] ,[Pf,Pf]);
Error(end+1)= err([P,Pf] ,[Pf,Pf]);

Error(end+1)= err([P;P] ,[Pf;Pf]);
Error(end+1)= err([Pf;P] ,[Pf;Pf]);
Error(end+1)= err([P;Pf] ,[Pf;Pf]);

Error(end+1)= err([P,P] ,[Pf,Pf]);
Error(end+1)= err([Pf,P] ,[Pf,Pf]);
Error(end+1)= err([P,Pf] ,[Pf,Pf]);

Error(end+1)= err(cat(3,P,P,P) ,cat(3,Pf,Pf,Pf));
Error(end+1)= err(cat(3,Pf,P) ,cat(3,Pf,Pf));
Error(end+1)= err(cat(3,P,Pf) ,cat(3,Pf,Pf));

Error(end+1)= err(cat(4,P,Pf) ,cat(4,Pf,Pf));
Error(end+1)= err(cat(5,P,Pf) ,cat(5,Pf,Pf));

%%Test of spfun

 f=@(x) cos(x).^2;

Error(end+1)= err( spfun(f,P) ,    reshape(full(spfun(f,Pf(:))) ,size(Pf)) );

%%Test of subsindex

   idx=1:3;
   Zf=rand(3);

Error(end+1)= err( Zf(ndSparse(idx)) , Zf(idx) );

%%Test of subsref

Error(end+1)= err( P(P<.5) ,  Pf(Pf<.5)   ); %logical indexing
Error(end+1)= err( Pf(P<.5) ,  Pf(Pf<.5)   );
Error(end+1)= err( P(Pf<.5) ,  Pf(Pf<.5)   );

Error(end+1)= err( P(:) ,  Pf(:)   );  %linear indexing
Error(end+1)= err( P(1) ,  Pf(1)   );
Error(end+1)= err( P(1:4) ,  Pf(1:4)   );
Error(end+1)= err( P((1:4).') ,  Pf((1:4).')   );
Error(end+1)= err( Vc(1:3) ,  Vcf(1:3)   );
Error(end+1)= err( Vr(1:3) ,  Vrf(1:3)   );

Error(end+1)= err( P(Vc(1:3)) ,  Pf(Vcf(1:3))   ); %indexing vectors test different shaping rules
Error(end+1)= err( Pf(Vc(1:3)) ,  Pf(Vcf(1:3))   );
Error(end+1)= err( P(Vcf(1:3)) ,  Pf(Vcf(1:3))   );
Error(end+1)= err( Pf(Vcf(1:3)) ,  Pf(Vcf(1:3))   );

Error(end+1)= err( P(2,1,2,2) ,  Pf(2,1,2,2)   ); %subscript indexing
Error(end+1)= err( P(2,:,2,2) ,  Pf(2,:,2,2)   );
Error(end+1)= err( P(:,1,:,2) ,  Pf(:,1,:,2)   );
Error(end+1)= err( P(:,2,2,2) ,  Pf(:,2,2,2)   );
Error(end+1)= err( P(1,2,2,:) ,  Pf(1,2,2,:)   );
Error(end+1)= err( P(:,:,:,2) ,  Pf(:,:,:,2)   );

Error(end+1)= err( P(2,1,2) ,  Pf(2,1,2)   ); %truncated subscript indexing
Error(end+1)= err( P(2,:,2) ,  Pf(2,:,2)   );
Error(end+1)= err( P(:,1,:) ,  Pf(:,1,:)   );
Error(end+1)= err( P(:,2,2) ,  Pf(:,2,2)   );

            lidx=logical([1 0]);

Error(end+1)= err( P(lidx,lidx,1:2) ,  Pf(lidx,lidx,1:2)   ); %combine all types of indexing
Error(end+1)= err( P(lidx,lidx,:) ,  Pf(lidx,lidx,:)   );
Error(end+1)= err( P(lidx,:,lidx) ,  Pf(lidx,:,lidx)   );

%%Test of subsasgn

        Z=P;
        Zf=Pf;

        Z(Z<.5)=6;   Zf(Zf<.5)=6;

Error(end+1)= err(  Z, Zf); %logical indexing

        Z(1)=7;   Zf(1)=7;

Error(end+1)= err(  Z, Zf); %scalar indexing

        Z(Zf<.5)=13; Zf(Zf<.5)=13;

Error(end+1)= err( Z ,  Zf   );

           Z(:)=999; Zf(:)=999;

Error(end+1)= err(  Z, Zf    );

           Z(1:4)=88;  Zf(1:4)=88;

Error(end+1)= err( Z ,  Zf   );

          Z((1:4).')=77;  Zf((1:4).')=77;

Error(end+1)= err(Z,Zf);

          d=rand;
          Z(2,1,2,2)=d;   Zf(2,1,2,2)=d;

Error(end+1)= err(Z,Zf); %subscript indexing

           d=rand;
          Z(2,:,2,2)=d;  Zf(2,:,2,2)=d;

Error(end+1)= err(Z,Zf);

          d=rand;
          Z(:,2,2,2)=d;  Zf(:,2,2,2)=d;

Error(end+1)= err(Z,Zf);

          d=rand;
          Z(:,:,:,2)=d;  Zf(:,:,:,2)=d;

Error(end+1)= err(Z,Zf);

           mm=size(Z,1)+3;
           nn=size(Z,2)+3;

          Z(mm,nn,:,:)=5;  Zf(mm,nn,:,:)=5; %matrix expansion test

Error(end+1)= err(Z,Zf);

          Z(:,1:2:end,:,:)=[];  Zf(:,1:2:end,:,:)=[]; %null assignment test

Error(end+1)= err(Z,Zf);

           Z=P; Zf=Pf;
           Z(:,mod(1:2:end,2)==1,:,:)=[];  %null assignment test logical indexing
           Zf(:,mod(1:2:end,2)==1,:,:)=[];

Error(end+1)= err(Z,Zf);

      Z=P; Zf=Pf;
      Z(:,:,:,:,1)=0;  Zf(:,:,:,:,1)=0; %matrix expansion with n-1 colons
      Z(:,:,:,:,1)=[];  Zf(:,:,:,:,1)=[]; %null assignment with n-1 colons

Error(end+1)= err(Z,Zf);

           Z=P; Zf=Pf;

           d=rand(size(Zf(:,:,1)));

           Z(:,:,1)=d;  Zf(:,:,1)=d; %nonscalar assignment

Error(end+1)= err( Z , Zf   );

            idx={':'};

            Z=P; Zf=Pf;
            d=rand(size(Zf(idx{:})));

            Z(idx{:})=d;  Zf(idx{:})=d;%nonscalar assignment

Error(end+1)= err(Z,Zf);

            lidx=logical([1 0]);

            idx={lidx,lidx,1:2};

            Z=P; Zf=Pf;
            d=rand(size(Zf(idx{:})));

            Z(idx{:})=d;  Zf(idx{:})=d;%combine all types of indexing

Error(end+1)= err( Z , Zf   );

%%Test of all(), any(), mean(), max/min

      Z=P>.5; Zf=Pf>.5;

      Z(:,1,1,1)=0; Zf(:,1,1,1)=0;

      Args={ {},{1},{2},{3},{4},{5} };
      for jj=1:length(Args)

          args=Args{jj};

Error(end+1)= err( all(Z,args{:})  , all(Zf,args{:}) );

Error(end+1)= err( any(Z,args{:})  , any(Zf,args{:}) );

Error(end+1)= err( mean(Z,args{:})  , mean(Zf,args{:}) );

Error(end+1)= err( all(P,args{:})  , all(Pf,args{:}) );

Error(end+1)= err( any(P,args{:})  , any(Pf,args{:}) );

Error(end+1)= err( mean(P,args{:})  , mean(Pf,args{:}) );
      end

        Args={ {1},{2},{3},{4},{5} };

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=max(P,[],args{:});
          [Zf,idxf]=max(Pf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( max(P,Q) , max(Pf,Qf) );
Error(end+1)= err( max(P,Qf) , max(Pf,Qf) );
Error(end+1)= err( max(Pf,Q) , max(Pf,Qf) );

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=min(P,[],args{:});
          [Zf,idxf]=min(Pf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( min(P,Q) , min(Pf,Qf) );
Error(end+1)= err( min(P,Qf) , min(Pf,Qf) );
Error(end+1)= err( min(Pf,Q) , min(Pf,Qf) );

       %%%%%%2D cases%%%%%

      Z=B>.5; Zf=Bf>.5;

      Z(:,1,1,1)=0; Zf(:,1,1,1)=0;

      Args={ {},{1},{2},{3},{4},{5} };
      for jj=1:length(Args)

          args=Args{jj};

Error(end+1)= err( all(Z,args{:})  , all(Zf,args{:}) );

Error(end+1)= err( any(Z,args{:})  , any(Zf,args{:}) );

Error(end+1)= err( mean(Z,args{:})  , mean(Zf,args{:}) );

      end

        Args={ {1},{2},{3},{4},{5} };

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=max(B,[],args{:});
          [Zf,idxf]=max(Bf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( max(B,C) , max(Bf,Cf) );
Error(end+1)= err( max(B,Cf) , max(Bf,Cf) );
Error(end+1)= err( max(Bf,C) , max(Bf,Cf) );

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=min(B,[],args{:});
          [Zf,idxf]=min(Bf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( min(B,C) , min(Bf,Cf) );
Error(end+1)= err( min(B,Cf) , min(Bf,Cf) );
Error(end+1)= err( min(Bf,C) , min(Bf,Cf) );

%Test of numel
Error(end+1)= err( numel(P) , numel(Pf) );

%Test of repmat, circshift

       N=ndims(Pf)+1;
       z=ones(1,N);

       for kk=0:N

           idx=nchoosek(1:N,kk);

           for jj=1:size(idx,1)

             arg2=z;
             arg2(idx(jj,:))=2;

Error(end+1)= err( repmat(P,arg2) , repmat(Pf,arg2) );
Error(end+1)= err( repmat(P,arg2-1) , repmat(Pf,arg2-1) );
Error(end+1)= err( circshift(P,arg2) , circshift(Pf,arg2) );
Error(end+1)= err( circshift(P,arg2-1) , circshift(Pf,arg2-1) );
Error(end+1)= err( circshift(A,arg2) , circshift(Af,arg2) );
Error(end+1)= err( circshift(A,arg2-1) , circshift(Af,arg2-1) );

           end
       end

              Z=P*0; Zf=full(Z);

Error(end+1)= err( repmat(Z,arg2) , repmat(Zf,arg2) );
Error(end+1)= err( circshift(Z,arg2) , circshift(Zf,arg2) );
Error(end+1)= err( repmat(Z,arg2-1) , repmat(Zf,arg2-1) );
Error(end+1)= err( circshift(Z,arg2-1) , circshift(Zf,arg2-1) );

%Test of squeeze, shiftdim

     N=ndims(Pf);
     z=ones(1,2*N); z(2:2:end)=size(P);
     Z=reshape(P,[1 1 z]); %add some singleton dimensions
     Zf=full(Z);

Error(end+1)= err( squeeze(Z) , squeeze(Zf) );

    [Z,n]=shiftdim(Z); [Zf,nf]=shiftdim(Zf);

Error(end+1)= err( Z , Zf);
Error(end+1)= err( n , nf);

Error(end+1)= err( shiftdim(Z,3) , shiftdim(Zf,3));
Error(end+1)= err( shiftdim(Z,-3) , shiftdim(Zf,-3));

%Test of bsxfun

     funcs={@plus,@minus,@times,@rdivide,@ldivide,@power,...
            @max,@min,@rem,@mod,@atan2,@hypot,@(a,b) a.^2-b};

     lfuncs={@eq,@ne,@lt,@le,@gt,@ge,@and,@or,@xor};

     HH=P;
     for qq=1:2

     %test non-logical funcs
     H=HH; Hf=full(H);
     Z=mean(H,1); Z=mean(Z,3);   Zf=full(Z);

     for kk=1:2
        for ii=1:length(funcs)

          fun=funcs{ii};

Error(end+1)= err( bsxfun(fun,Z,H) ,     bsxfun(fun,Zf,Hf)  );
Error(end+1)= err( bsxfun(fun,Z,Hf) ,     bsxfun(fun,Zf,Hf)  );
Error(end+1)= err( bsxfun(fun,Zf,H) ,     bsxfun(fun,Zf,Hf)  );
Error(end+1)= err( bsxfun(fun,three,H) ,     bsxfun(fun,three,Hf)  );
Error(end+1)= err( bsxfun(fun,Zf,three) ,     bsxfun(fun,Zf,three)  );

       end

      Z=mean(H,2); Zf=full(Z);

      end

      %test logical lfuncs
      H=(HH<0.5); Hf=full(H);
      Z=mean(H,1); Z=mean(Z,3);
      Z=(Z>0.5);  Zf=full(Z);

      onebit=(rand>=0.5);

      for kk=1:2
       for ii=1:length(lfuncs)

          fun=lfuncs{ii};

Error(end+1)= err( bsxfun(fun,Z,H) ,     bsxfun(fun,Zf,Hf)  );
Error(end+1)= err( bsxfun(fun,Z,Hf) ,     bsxfun(fun,Zf,Hf)  );
Error(end+1)= err( bsxfun(fun,Zf,H) ,     bsxfun(fun,Zf,Hf)  );
Error(end+1)= err( bsxfun(fun,onebit,H) ,     bsxfun(fun,onebit,Hf)  );
Error(end+1)= err( bsxfun(fun,Zf,onebit) ,     bsxfun(fun,Zf,onebit)  );

       end

      Z=mean(H,2);
      Z=(Z>0.5); Zf=full(Z);

      end

      HH=repmat(P,[1,1,1,0]);

     end

%Test of ndSparse.build, nzmax, nonzeros

      nzm=10;
      Z=ndSparse.build([1 1 1; 2 2 2; 3 3 3; 4 4 4], 1:4,[4 4 6],nzm);
      Zf=zeros(4,4,6); for ii=1:4, Zf(ii,ii,ii)= ii; end

Error(end+1)= err( Z , Zf  );
Error(end+1)= err( nzmax(Z) , nzm );
Error(end+1)= err( nonzeros(Z) , nonzeros(Zf) );

      nzm=10;
      Z=ndSparse.build([1 1 1; 2 2 2; 3 3 3; 4 4 4], 5,[4 4 6],nzm);
      Zf=zeros(4,4,6); for ii=1:4, Zf(ii,ii,ii)= 5; end

Error(end+1)= err( Z , Zf  );
Error(end+1)= err( nzmax(Z) , nzm );
Error(end+1)= err( nonzeros(Z) , nonzeros(Zf) );

       Z=ndSparse.build([4 4 6]);
       Zf=zeros(4,4,6);

Error(end+1)= err( Z , Zf  );
Error(end+1)= err( ndSparse.build([],[],[4 4 6]) , Zf  );
Error(end+1)= err( nzmax(Z) , 1 );
Error(end+1)= err( nonzeros(Z) , nonzeros(Zf) );

%Test of ndSparse.spalloc

      nzm=10;
      Z=ndSparse.spalloc([3,5,3],nzm);
      Zf=full(Z);

Error(end+1)= err( nzmax(Z) , nzm );
Error(end+1)= err( Z , Zf );

%Test of ndSparse.accumarray

         args={[1 1 1; 2 2 2; 3 3 3;2 2 2], [5,4,3,2]};
         Z=ndSparse.accumarray(args{:});
         Zf=accumarray(args{:});

Error(end+1)= err( Z , Zf );

         args={[1 1 1; 2 2 2; 3 3 3;2 2 2], [5,4,3,2],[4 4 4]};
         Z=ndSparse.accumarray(args{:});
         Zf=accumarray(args{:});

Error(end+1)= err( Z , Zf );

         args={[1;2;3;2], [5,4,3,2]};
         Z=ndSparse.accumarray(args{:});
         Zf=accumarray(args{:});

Error(end+1)= err( Z , Zf );

         args={[1;2;3;2], [5,4,3,2],[4,1]};
         Z=ndSparse.accumarray(args{:});
         Zf=accumarray(args{:});

Error(end+1)= err( Z , Zf );

         args={[1;2;3;2], [5,4,3,2],[4,1],@prod};
         Z=ndSparse.accumarray(args{:});
         Zf=accumarray(args{:});

Error(end+1)= err( Z , Zf );

%Test of spones

          Z=ndSparse.build([1,1,1], 5,[3,5,3]);
          Zf=full(Z); Zf(~~Zf)=1;

Error(end+1)= err( spones(Z) , Zf );

%Test of length

Error(end+1)= err( length(P) , length(Pf) );

%Test of triu,triul

Error(end+1)= err( triu(A) , triu(Af) );
Error(end+1)= err( triu(A,1) , triu(Af,1) );
Error(end+1)= err( triu(A,-1) , triu(Af,-1) );

Error(end+1)= err( tril(A) , tril(Af) );
Error(end+1)= err( tril(A,1) , tril(Af,1) );
Error(end+1)= err( tril(A,-1) , tril(Af,-1) );

Error(end+1)= err( triu(S) , triu(Sf) );
Error(end+1)= err( triu(S,1) , triu(Sf,1) );
Error(end+1)= err( triu(S,-1) , triu(Sf,-1) );

Error(end+1)= err( tril(S) , tril(Sf) );
Error(end+1)= err( tril(S,1) , tril(Sf,1) );
Error(end+1)= err( tril(S,-1) , tril(Sf,-1) );

%Test of flipdim

Error(end+1)= err( flipdim(P,1) , flipdim(Pf,1) );
Error(end+1)= err( flipdim(P,2) , flipdim(Pf,2) );
Error(end+1)= err( flipdim(P,3) , flipdim(Pf,3) );

%Test of  fliplr

   Z=P(:,:,:,1); Zf=full(Z);   for ii=1:size(Zf,3), Zf(:,:,ii)=fliplr( Zf(:,:,ii) ); end

Error(end+1)= err( fliplr(Z) ,  Zf );

%Test of flipud

   Z=P(:,:,:,1);  Zf=full(Z);    for ii=1:size(Zf,3), Zf(:,:,ii)=flipud( Zf(:,:,ii) ); end

Error(end+1)= err( flipud(Z) ,  Zf );

%Test of rot90

  Z=P(:,:,:,1);  Zf=full(Z);    for ii=1:size(Zf,3), Zf(:,:,ii)=rot90( Zf(:,:,ii) ); end

Error(end+1)= err( rot90(Z) ,  Zf );

             for kk=-5:5

   Z=P(:,:,:,1);  Zf=full(Z);    for ii=1:size(Zf,3), Zf(:,:,ii)=rot90( Zf(:,:,ii) ,kk); end

Error(end+1)= err( rot90(Z,kk) , Zf);

             end

%Test of convn

            args={{},{'full'},{'same'},{'valid'}};
            Zs=sparse([0;0;1;2]); Zf=full(Zs);

            for ii=1:length(args)

Error(end+1)= err( convn(P,P,args{ii}{:}) , convn(Pf,Pf,args{ii}{:}) );
Error(end+1)= err( convn(P,A,args{ii}{:}) , convn(Pf,Af,args{ii}{:}) );

Error(end+1)= err( convn(Pf,P,args{ii}{:}) , convn(Pf,Pf,args{ii}{:}) );
Error(end+1)= err( convn(Pf,A,args{ii}{:}) , convn(Pf,Af,args{ii}{:}) );

Error(end+1)= err( convn(P,Pf,args{ii}{:}) , convn(Pf,Pf,args{ii}{:}) );
Error(end+1)= err( convn(P,Af,args{ii}{:}) , convn(Pf,Af,args{ii}{:}) );

Error(end+1)= err( convn(P,Zs,args{ii}{:}) , convn(Pf,Zf,args{ii}{:}) );
Error(end+1)= err( convn(Zs,P,args{ii}{:}) , convn(Zf,Pf,args{ii}{:}) );

Error(end+1)= err( convn(P,Zf,args{ii}{:}) , convn(Pf,Zf,args{ii}{:}) );
Error(end+1)= err( convn(Zf,P,args{ii}{:}) , convn(Zf,Pf,args{ii}{:}) );

            end

            Hf={rand(2,1),rand(1,3)};
            H=Hf;
            H{1}=ndSparse(H{1});

            Zf=diag([0 0 0 0 0 1 1]);
            Zs=sparse(Zf);

           for ii=1:length(args)

Error(end+1)= err( convn(H{:},A,args{ii}{:}) , conv2(Hf{:},Af,args{ii}{:}) );

Error(end+1)= err( convn(H{:},Zs,args{ii}{:}) , conv2(Hf{:},Zf,args{ii}{:}) );

Error(end+1)= err( convn(H{:},Zf,args{ii}{:}) , conv2(Hf{:},Zf,args{ii}{:}) );

           end

%%Test of allml(), anyml(), meanml(), maxml/minml, summl, catml

      Z=P>.5; Zf=Pf>.5;

      Z(:,1,1,1)=0; Zf(:,1,1,1)=0;

      Args={ {},{1},{2},{3},{4},{5} };
      for jj=1:length(Args)

          args=Args{jj};

Error(end+1)= err( allml(Z,args{:})  , all(Zf,args{:}) );

Error(end+1)= err( anyml(Z,args{:})  , any(Zf,args{:}) );

Error(end+1)= err( meanml(Z,args{:})  , mean(Zf,args{:}) );

Error(end+1)= err( allml(P,args{:})  , all(Pf,args{:}) );

Error(end+1)= err( anyml(P,args{:})  , any(Pf,args{:}) );

Error(end+1)= err( meanml(P,args{:})  , mean(Pf,args{:}) );
      end

        Args={ {1},{2},{3},{4},{5} };

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=maxml(P,[],args{:});
          [Zf,idxf]=max(Pf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( maxml(P,Q) , max(Pf,Qf) );
Error(end+1)= err( maxml(P,Qf) , max(Pf,Qf) );
Error(end+1)= err( maxml(Pf,Q) , max(Pf,Qf) );

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=minml(P,[],args{:});
          [Zf,idxf]=min(Pf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( minml(P,Q) , min(Pf,Qf) );
Error(end+1)= err( minml(P,Qf) , min(Pf,Qf) );
Error(end+1)= err( minml(Pf,Q) , min(Pf,Qf) );

       %%%%%%2D cases%%%%%

      Z=B>.5; Zf=Bf>.5;

      Z(:,1,1,1)=0; Zf(:,1,1,1)=0;

      Args={ {},{1},{2},{3},{4},{5} };
      for jj=1:length(Args)

          args=Args{jj};

Error(end+1)= err( allml(Z,args{:})  , all(Zf,args{:}) );

Error(end+1)= err( anyml(Z,args{:})  , any(Zf,args{:}) );

Error(end+1)= err( meanml(Z,args{:})  , mean(Zf,args{:}) );

      end

        Args={ {1},{2},{3},{4},{5} };

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=maxml(B,[],args{:});
          [Zf,idxf]=max(Bf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( maxml(B,C) , max(Bf,Cf) );
Error(end+1)= err( maxml(B,Cf) , max(Bf,Cf) );
Error(end+1)= err( maxml(Bf,C) , max(Bf,Cf) );

      for jj=1:length(Args)

          args=Args{jj};

          [Z,idx]=minml(B,[],args{:});
          [Zf,idxf]=min(Bf,[],args{:});

Error(end+1)= err( Z ,Zf );
Error(end+1)= err( idx ,idxf );

      end

Error(end+1)= err( minml(B,C) , min(Bf,Cf) );
Error(end+1)= err( minml(B,Cf) , min(Bf,Cf) );
Error(end+1)= err( minml(Bf,C) , min(Bf,Cf) );

Error(end+1)= err( summl(P,1)  , sum(Pf,1) );
Error(end+1)= err( summl(Q,2)  , sum(Qf,2) );
Error(end+1)= err( summl(P)    , sum(Pf) );
Error(end+1)= err( summl(P,'native')    , sum(Pf,'native') );
Error(end+1)= err( summl(P,3,'double')    , sum(Pf,3,'double') );
Error(end+1)= err( summl(P,4)    , sum(Pf,4) );
Error(end+1)= err( summl(P,5)    , sum(Pf,5) );

Error(end+1)= err(catml(3,P,P,P) ,cat(3,Pf,Pf,Pf));
Error(end+1)= err(catml(3,Pf,P) ,cat(3,Pf,Pf));
Error(end+1)= err(catml(3,P,Pf) ,cat(3,Pf,Pf));

Error(end+1)= err(catml(4,P,Pf) ,cat(4,Pf,Pf));
Error(end+1)= err(catml(5,P,Pf) ,cat(5,Pf,Pf));

       N=ndims(Pf)+1;
       z=ones(1,N);

       for kk=0:N

           idx=nchoosek(1:N,kk);

           for jj=1:size(idx,1)

             arg2=z;
             arg2(idx(jj,:))=2;

Error(end+1)= err( circshiftml(P,arg2) , circshift(Pf,arg2) );
Error(end+1)= err( circshiftml(P,arg2-1) , circshift(Pf,arg2-1) );
Error(end+1)= err( circshiftml(A,arg2) , circshift(Af,arg2) );
Error(end+1)= err( circshiftml(A,arg2-1) , circshift(Af,arg2-1) );

           end
       end

              Z=P*0; Zf=full(Z);

Error(end+1)= err( circshiftml(Z,arg2) , circshift(Zf,arg2) );
Error(end+1)= err( circshiftml(Z,arg2-1) , circshift(Zf,arg2-1) );

%%%%%%%%%%%%%%%%%%%%%%%END OF TESTS%%%%%%%%%%%%%%%%%%%%%%%%%%

MAX_ERROR=max(Error);

if any(~isfinite(Error)), warning 'There were improper Error values'; keyboard; end

disp(['Maximum observed error was   ' num2str(MAX_ERROR) ' percent.'])

if nargout, varargout{1}=Error; end

function errval=DiscrepancyMeasure(X,Y,TOL,CHECKTYPES)

  if Discrepancy(0,Y)
    errval=Discrepancy(X,Y)/Discrepancy(0,Y)*100; %normalize
  else
    errval=Discrepancy(X,Y);
  end

  isndSparse=@(c) ~isempty(strfind(class(c),'ndSparse'));

  if CHECKTYPES
    if ~isndSparse(X),
       warning(['X is not ndSparse class, but rather class ' class(X)])
    end
  end

  if errval>TOL || ~isfinite(errval),
      disp ' '; disp 'Discrepancy detected'
      errval,
      x=full(X); y=full(Y);
      keyboard;
  end

function errval=Discrepancy(X,Y)
%Primary error measurement function

  fin=@(a) reshape( a(isfinite(a)),[],1);
  nonfin=@(a) reshape(  a(~isfinite(a))  ,[],1);

  x=full(X); y=full(Y);

  if ( isequal(x,y));
     errval=0; return
  elseif xor(isempty(x),isempty(y))
      errval=1; return
  end

  errval= norm( fin(x-y) , inf)+...
       ~isequalwithequalnans(nonfin(x),nonfin(y))*...
       ~isempty([nonfin(x);nonfin(y)])+ ~isequal(size(X),size(Y));

function out=srand(varargin)

   out=rand(varargin{:})-.5;

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% l1_tv_restore1
% illustrate image restoration using a l1-type data fit term
% (for noise robustness) and a (anisotropic) TV-type regularizer
% based on l1_regress_example.m
% Copyright 2010-05-25, Jeff Fessler, University of Michigan

if ~isvar('xtrue')
	f.dir = [path_find_dir('mri') '/../data/mri/'];
	f.xtrue = [f.dir 'brainweb_t1.jpg'];
	xtrue = single(imread(f.xtrue)');
	xtrue = xtrue(2:end-1,2:end-1);
	im(xtrue), cbar
end

	snr = @(a,b) 20 * log10(norm(a(:))/norm(b(:)));

if ~isvar('A')
%	psf = fspecial('Gaussian', [7 7], 5)
	psf = gaussian_kernel(5 * sqrt(log(256)), 3);
	psf = psf * psf'; psf = psf / sum(psf(:));
	mask = true(size(xtrue));
	A = Gblur(mask, 'psf', psf);

	yb = A * xtrue;
	im(yb), cbar
	f.snr_db = 25;
	sig = 10^(-f.snr_db/20) * norm(yb(:)) / sqrt(numel(yb));

	rng(0)
	yi = yb + sig * randn(size(yb));
	printm('data snr = %g dB', snr(yb, yi-yb))

	im(yi), cbar
end

if ~isvar('xh')
	xh = l1_tv_restore1_fun(yi, A, 'l2b', -2);
	xh = embed(xh, mask);
	im(xh), cbar
end

	printm('raw snr = %g dB', snr(xtrue, yi - xtrue))
	printm('new snr = %g dB', snr(xtrue, xh - xtrue))

if 0
	beta = 1e-3;
	gamma = 1e-6;
	mu = 1/25; % for 25% random-valued noise
%	sqrt((y - A f)^2 + gamma) + mu * sqrt((Cx f)^2 + (Cy *f)^2 + beta)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function y=bvpme(t,x,flag,U)
A=[1-2*cos(2*t) 0 1+2*sin(2*t);
    0 2 0;
    -1+2*sin(2*t) 0 1+2*cos(2*t);];
q=exp(t)*[-1+2*cos(2*t)-2*sin(2*t);
    -1;
    1-2*cos(2*t)-2*sin(2*t);];
y = A*x+q;
end

author: Sable

repo: mcbench-benchmarks


text:
function vs = sphere_triangle_vertices_to_centroid ( r, v1, v2, v3 )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_VERTICES_TO_CENTROID gets a spherical triangle centroid in 3D.
%
%  Discussion:
%
%    A sphere centered at 0 in 3D satisfies the equation:
%
%      X*X + Y*Y + Z*Z = R*R
%
%    A spherical triangle is specified by three points on the sphere.
%
%    The (true) centroid of a spherical triangle is the point
%
%      VT = (XT,YT,ZT) = Integral ( X, Y, Z ) dArea / Integral 1 dArea
%
%    Note that the true centroid does NOT, in general, lie on the sphere.
%
%    The "flat" centroid VF is the centroid of the planar triangle defined by
%    the vertices of the spherical triangle.
%
%    The "spherical" centroid VS of a spherical triangle is computed by
%    the intersection of the geodesic bisectors of the triangle angles.
%    The spherical centroid lies on the sphere.
%
%    VF, VT and VS lie on a line through the center of the sphere.  We can
%    easily calculate VF by averaging the vertices, and from this determine
%    VS by normalizing.
%
%    Of course, we still will not have actually computed VT, which lies
%    somewhere between VF and VS!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real V1(3), V2(3), V3(3), the vertices of the triangle.
%
%    Output, real VS(3), the coordinates of the "spherical
%    centroid" of the spherical triangle.
%
  dim_num = 3;

  vs(1:dim_num) = ( v1(1:dim_num) + v2(1:dim_num) + v3(1:dim_num) ) / 3.0;

  norm = sqrt ( sum ( vs(1:dim_num).^2 ) );

  vs(1:dim_num) = r * vs(1:dim_num) / norm;

  return
end

author: johannesgerer

repo: jburkardt-m


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

%*****************************************************************************80
%
%% LD0302 computes the 302 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(302,1);
  y = zeros(302,1);
  z = zeros(302,1);
  w = zeros(302,1);
  a = 0.0;
  b = 0.0;
  v = 0.8545911725128148E-03;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.3599119285025571E-02;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.3515640345570105;
  v = 0.3449788424305883E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6566329410219612;
  v = 0.3604822601419882E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4729054132581005;
  v = 0.3576729661743367E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.9618308522614784E-01;
  v = 0.2352101413689164E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2219645236294178;
  v = 0.3108953122413675E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.7011766416089545;
  v = 0.3650045807677255E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2644152887060663;
  v = 0.2982344963171804E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.5718955891878961;
  v = 0.3600820932216460E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.2510034751770465;
  b = 0.8000727494073952;
  v = 0.3571540554273387E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1233548532583327;
  b = 0.4127724083168531;
  v = 0.3392312205006170E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end


text:
function slwritetext(T, filename)
%SLWRITETEXT Writes a text to file
%
% $ Syntax $
%   - slwritetext(T, filename)
%
% $ Arguments $
%   - T             the cell array representing the text
%   - filename      the path of the file to write to
%
% $ Description $
%   - slwritetext(T, filename) Writes in a text file stored in a cell array
%   to a text file specified by filename.
%
% $ History $
%   Created by Dahua Lin, on 2005-06-02
%

if ~iscell(T)
    error('T should be a text cell array');
end
fid = fopen(filename, 'w');
if (fid <= 0)
    error(['Fail to open ', filename]);
end

n = length(T);
for i = 1 : n
    fprintf(fid, '%s\r\n', T{i});
end

fclose(fid);

author: lmthang


text:
function fixed = p15_fixed_points ( m, fixed_num )

%*****************************************************************************80
%
%% P15_FIXED_POINTS returns the fixed points in problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer FIXED_NUM, the number of fixed points.
%
%    Output, real FIXED(M,FIXED_NUM), the fixed points.
%
  fixed = [    ...
   -8.0, -1.0; ...
    2.0, -1.0; ...
    2.0,  0.0; ...
    8.0,  0.0; ...
    8.0,  1.0; ...
   -2.0,  1.0; ...
   -2.0,  0.0; ...
   -8.0,  0.0 ]';

  return
end


text:
function[]=makefigs_maternoise
%MAKEFIGS_MATERNOISE  Makes a sample figure for MATERNOISE.

N=1000;
alpha=[0.6 1 1.5 2 3 4];
h=[.01 .02 .05 .2 1];
[alpha,h]=meshgrid(alpha,h);
h=h.*alpha;

rng(1);  %set seed
z=maternoise(1,N,10,alpha,h);
z=z./vrep(std(z,1,1),size(z,1),1);  %Set to unit std
y=detrend(cumsum(z),'constant');
y=y./vrep(std(y,1,1),size(y,1),1);  %Set to unit std

%This is just to make an offset
[xo,yo]=meshgrid(1:6,1:5);
zo=xo(:)'+1i*yo(:)';
zo=vrep(zo,length(z),1);

figure,
plot(y+zo*3),axis equal,noxlabels,noylabels
xlabel('Increasing $\alpha \rightarrow$'),ylabel('Increasing $\lambda \rightarrow$')
set(gca,'xticklabel',[]),set(gca,'xticklabel',[])
title('Example of Matern Processes')

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function p = get_params(gait)

p.predHorizon = 6;
p.simTimeStep = 1/200;
p.Tmpc = 4/100;             % MPC prediction step time
p.gait = gait;
p.Umax = 50;
p.decayRate = 1;
p.freq = 30;
p.Rground = eye(3);
p.Qf = diag([1e5 2e5 3e5 5e2 1e3 150 1e3 1e4 800 40 40 10]);

% ---- gait ----
if gait == 1                % 1 - bound
    p.Tst = 0.1;
    p.Tsw = 0.18;
    p.predHorizon = 7;
    p.simTimeStep = 1/100;
    p.Tmpc = 2/100;
    p.decayRate = 1;
    p.R = diag(repmat([0.1 0.1 0.1]',[4,1]));
    p.Q = diag([5e4 2e4 1e6 4e3 5e2 5e2 1e4 5e4 1e3 1e2 5e2 1e2]);
    p.Qf = diag([2e5 5e4 5e6 8e3 5e2 5e2 1e4 5e4 5e3 1e2 1e2 1e2]);
elseif gait == 2            % 2 - pacing
    p.Tst = 0.12;
    p.Tsw = 0.12;
    p.R = diag(repmat([0.1 0.2 0.1]',[4,1]));
    p.Q = diag([5e3 5e3 9e4 5e2 5e2 5e2 7e3 7e3 7e3 5e1 5e1 5e1]);
elseif gait == 3            % 3 - gallop
    p.Tst = 0.08;
    p.Tsw = 0.2;
    p.R = diag(repmat([0.1 0.2 0.1]',[4,1]));
    p.Q = diag([3e3 3e3 4e6 5e2 1e3 150 1e4 1e4 800 1e2 5e1 5e1]);
elseif gait == 4            % 4 - trot run
    p.Tst = 0.12;
    p.Tsw = 0.2;
    p.Tmpc = 3/100;
    p.predHorizon = 6;
    p.decayRate = 1;
    p.R = diag(repmat([0.1 0.18 0.08]',[4,1]));
    p.Q = diag([1e5 1e5 1e5 1e3 1e3 1e3 2e3 1e4 800 100 40 10]);
    p.Qf = diag([1e5 1.5e5 2e4 1.5e3 1e3 100 2e3 2e3 800 100 60 10]);
elseif gait == 5            % 4 - crawl
    p.Tst = 0.3;
    p.Tsw = 0.1;
    p.R = diag(repmat([0.1 0.2 0.1]',[4,1]));
    p.Q = diag([5e5 5e5 9e5 5 5 5 3e3 3e3 3e3 3 3 3]);
else                        % 0 - trot
    p.predHorizon = 6;
    p.simTimeStep = 1/100;
    p.Tmpc = 8/100;
    p.Tst = 0.3;
    p.Tsw = 0.15;
    p.R = diag(repmat([0.1 0.2 0.1]',[4,1]));
    p.Q = diag([1e5 2e5 3e5 5e2 1e3 1e3 1e3 1e4 800 40 40 10]);
    p.Qf = p.Q;
end

%% Physical Parameters
p.mass = 5.5;
p.J = diag([0.026,0.112,0.075]);
p.g = 9.81;
p.mu = 1;       % friction coefficient
p.z0 = 0.2;     % nominal COM height
p.pf34 = [[0.15;0.094;0],[0.15;-0.094;0],[-0.15;0.094;0],[-0.15;-0.094;0]];

p.L = 0.301;    % body length
p.W = 0.088;    % body width
p.d = 0.05;     % ABAD offset
p.h = 0.05;     % body height
p.l1 = 0.14;    % link1 length
p.l2 = 0.14;    % link2 length

%% Swing phase
p.Kp_sw = 300;  % Kp for swing phase

%% color
p.body_color    = [42 80 183]/255;
p.leg_color     = [7 179 128]/255;
p.ground_color  = [195 232 243]/255;

author: YanranDing

repo: RF-MPC

save_path: github-repos/MATLAB/YanranDing-RF-MPC


text:
function T67=T67(h)

T67=-1.96e-3*h+373.2;


text:
function y = R3RSR(x,batch)
% R3RSR will R3R a vector x follow by a repeated splitting
%    (repeated SR) of the results from R3R.
%
% Usage: y = R3RSR(x)
%
%   where  x is the input vector with length at least 3
%          y is the return vector of the same length as x
%

% See John Wilder Tukey, "EXPLORATORY DATA ANALYSIS".
% Addison-Wesley Publishing Co. 1977. Chpt 7A, 7D

if min(size(x)) ~= 1, error('Input must be a vector'); end
x = x(:); xlen = length(x); cont = 1; y = x;
if nargin<2, batch = 0; end
if xlen < 3,
   if ~batch, warning('Input vector must be at least 3 elements long'); end
   cont = 0;
end
xtemp1 = R3R(x,batch);
xtemp2 = SR(xtemp1,batch);
while any(xtemp2-xtemp1),
   xtemp1 = xtemp2;
   xtemp2 = SR(xtemp1,batch);
end

return


text:
function [posterior] = tapas_sem_hier_prepare_posterior(data, model, ...
    inference, states)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

T = model.graph{1}.htheta.T;

posterior = struct('data', data, 'model', model, 'inference', inference, ...
    'samples_theta', [], 'fe', [], 'llh', []);

np = numel(states);

theta = cell(floor(np/inference.thinning), 1);
nc = 1;
for i = 1:inference.thinning:np
    theta{nc} = states{i}.graph{2}(:, end);
    nc = nc + 1;
end

% Contains all the terms related to the likelihood
cllh = cell(2, 1);

[ns, nc] = size(states{1}.llh{1});

% Log likelihood under the posteriors
llh = zeros(ns, nc, floor(np/inference.thinning));
nc = 1;

for i = 1:inference.thinning:np
    llh(:, :, nc) = states{i}.llh{1};
    nc = nc + 1;
end

cllh{1} = llh;

posterior.llh = cllh;
if size(T, 2) > 1
    fe = trapz(T(1, :), mean(squeeze(sum(cllh{1}, 1)), 2));
else
    fe = nan;
end

posterior.fe = fe;

% Compute the WAIC / sum of the variance of the log likelihood
% (gradient of the FE) Take only the last chain (second dimension)
[waic, accuracy] = tapas_waic(squeeze(llh(:, end, :)));

posterior.waic = waic;
posterior.accuracy = accuracy;

posterior.T = T;

theta = horzcat(theta{:});
jm = model.graph{1}.htheta.model.jm;
p0 = model.graph{1}.htheta.model.p0;

for i = 1:numel(theta)
    theta{i} = p0 + jm * theta{i};
end

posterior.ps_theta = theta;

% Add the summaries automatically
posterior.summary = tapas_sem_posterior_summary(posterior);

end

author: translationalneuromodeling

repo: tapas


text:
%XML_PARSEANY  parses XML strings and convert content into Matlab structure
%
% SYNTAX
%           V = xml_parseany( xmlstr )
%
% INPUT
%   xmlstr   xml string, possibly from file with function xml_load.m
%
% OUTPUT
%   v        Matlab variable or cell array.
%
% RELATED
%   xml_format, xml_formatany, xml_parse, xml_load, xml_save, (xmlread, xmlwrite)

% Copyright (C) 2002-2005, University of Southampton
% Author: Dr Marc Molinari <m.molinari@soton.ac.uk>
% $Revision: 1.1 $ $Date: 2005/04/15 17:12:14 $ $Tag$

author: Cloud-CV

repo: object-proposals


text:
%%% Computes the extrinsic parameters for all the active calibration images

check_active_images;

N_points_views = zeros(1,n_ima);

for kk = 1:n_ima,

    if exist(['x_' num2str(kk)]),

        eval(['x_kk = x_' num2str(kk) ';']);
        eval(['X_kk = X_' num2str(kk) ';']);

        if (isnan(x_kk(1,1))),
            if active_images(kk),
                fprintf(1,'Warning: Cannot calibrate with image %d. Need to extract grid corners first.\n',kk)
                fprintf(1,'         Set active_images(%d)=1; and run Extract grid corners.\n',kk)
            end;
        end;
        if active_images(kk),
            N_points_views(kk) = size(x_kk,2);
            [omckk,Tckk] = compute_extrinsic_init(x_kk,X_kk,fc,cc,kc,alpha_c);
            [omckk,Tckk,Rckk,JJ_kk] = compute_extrinsic_refine(omckk,Tckk,x_kk,X_kk,fc,cc,kc,alpha_c,20,thresh_cond);
            if check_cond,
                if (cond(JJ_kk)> thresh_cond),
                    active_images(kk) = 0;
                    omckk = NaN*ones(3,1);
                    Tckk = NaN*ones(3,1);
                    fprintf(1,'\nWarning: View #%d ill-conditioned. This image is now set inactive.\n',kk)
                    desactivated_images = [desactivated_images kk];
                end;
            end;
            if isnan(omckk(1,1)),
                %fprintf(1,'\nWarning: Desactivating image %d. Re-activate it later by typing:\nactive_images(%d)=1;\nand re-run optimization\n',[kk kk])
                active_images(kk) = 0;
            end;
        else
            omckk = NaN*ones(3,1);
            Tckk = NaN*ones(3,1);
        end;

    else

        omckk = NaN*ones(3,1);
        Tckk = NaN*ones(3,1);

        if active_images(kk),
            fprintf(1,'Warning: Cannot calibrate with image %d. Need to extract grid corners first.\n',kk)
            fprintf(1,'         Set active_images(%d)=1; and run Extract grid corners.\n',kk)
        end;

        active_images(kk) = 0;

    end;

    eval(['omc_' num2str(kk) ' = omckk;']);
    eval(['Tc_' num2str(kk) ' = Tckk;']);

end;

check_active_images;

repo: ekfmonoslam


text:
function opengz(gz)
% Open imtool3D by double click on a .nii.gz file in Matlab file browser
if strcmp(gz(max(end-6,1):end),'.nii.gz')
    imtool3D_nii_3planes(gz);
else
    gunzip(gz);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function value = p21_exact ( dim_num )

%*****************************************************************************80
%
%% P21_EXACT returns the exact integral for problem 21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  c = 0.0;
  c = p21_r8 ( 'G', 'C', c );

  e = [];
  e = p21_i4vec ( 'G', 'E', dim_num, e );

  if ( any ( mod ( e(1:dim_num), 2 ) == 1 ) )
    value = 0.0;
    return
  end

  value = 2.0 * c;
  for i = 1 : dim_num
    arg = ( e(i) + 1 ) / 2.0;
    value = exact * gamma ( arg );
  end

  arg = ( sum ( e(1:dim_num) ) + dim_num ) / 2.0;
  value = value / gamma ( arg );

  return
end


text:
function solution = createSolution(dim,N)

    solution.lambda = zeros(dim.lambda,   N);
    solution.mu     = zeros(dim.mu,       N);
    solution.u      = zeros(dim.u,        N);
    solution.x      = zeros(dim.x,        N);
    solution.z      = zeros(dim.z,        N);
    solution.LAMBDA = zeros(dim.x, dim.x, N);

end

author: deng-haoyang

repo: ParNMPC


text:
function [y,deriv] = sums_of_squares(w,m)
% This is a MV2DF. See MV2DF_API_DEFINITION.readme.
% Does:
% (i) Reshapes w to m-by-n.
% (ii) Computes sum of squares of each of n columns.
% (iii) Transposes to output n-vector.

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)sums_of_squares(w,m);
    return;
end

if isa(w,'function_handle')
    outer = sums_of_squares([],m);
    y = compose_mv(outer,w,[]);
    return;
end

M = reshape(w,m,[]);
y = sum(M.^2,1);
y = y(:);

deriv = @(g2) deriv_this(g2,M);

function [g,hess,linear] = deriv_this(g2,M)

g = 2*bsxfun(@times,M,g2.');
g = g(:);
linear = false;
hess = @(d) hess_this(d,g2,M);

function [h,Jv] = hess_this(d,g2,M)
h = deriv_this(g2,reshape(d,size(M)));
if nargout>1
    Jv = 2*sum(reshape(d,size(M)).*M,1);
    Jv = Jv(:);
end

function test_this()
f = sums_of_squares([],10);
test_MV2DF(f,randn(10*4,1));

author: nesl


text:
function [ shape2D, global_params, local_params, final_lhood, landmark_lhoods, view_used, early_term ] = Fitting_from_bb_multi_hyp( Image, DepthImage, bounding_box, PDM, patches, clmParams, views, early_term_params)

    num_points = numel(PDM.M)/3;

    shapes = zeros(num_points, 2, size(views,1));
    ls = zeros(size(views,1),1);
    lmark_lhoods = zeros(num_points,size(views,1));
    views_used = zeros(size(views,1),1);
    globals = cell(size(views,1),1);
    locals = cell(size(views,1),1);

    % If early termination parameters are defined use clever multi-hypothesis setup
    % otherwise evaluate all of them in a simplistic manner
    if(nargin > 7)

        patches_small = patches(1);
        patches_large = patches(2:end);

        clmParams_small = clmParams;
        clmParams_small.numPatchIters = 1;
        clmParams_large = clmParams;
        clmParams_large.regFactor = clmParams_large.regFactor(2:end);
        clmParams_large.sigmaMeanShift = clmParams_large.sigmaMeanShift(2:end);
        clmParams_large.tikhonov_factor = clmParams_large.tikhonov_factor(2:end);
        clmParams_large.numPatchIters = 3;
        clmParams_large.window_size = clmParams_large.window_size(2:end,:);

        early_term = 0;
        % Find the best orientation
        for v = 1:size(views,1)
            [shapes(:,:,v),globals{v},locals{v},ls(v),lmark_lhoods(:,v),views_used(v)]...
                = Fitting_from_bb(Image, DepthImage, bounding_box, PDM, patches_small, clmParams_small, 'orientation', views(v,:));

            ls(v) = ls(v) * early_term_params.weights_scale(views_used(v)) + early_term_params.weights_add(views_used(v));

            if(ls(v) > early_term_params.cutoffs(views_used(v)))
                early_term = v;
                [shape,global_params,local_params,final_lhood,landmark_lhoods,view_used]...
                    = Fitting_from_bb(Image, DepthImage, bounding_box, PDM, patches_large, clmParams_large, 'gparam', globals{v}, 'lparam', locals{v});
                break;
            end
        end

        if(early_term == 0)

            % Discard the really unreliable starts
            [~, ids] = sort(ls);
            to_refine = ids(end-3:end);

            % If there are any starts left try refining them in a
            % multiscale approach
            if(~isempty(to_refine))
                shapes = zeros(num_points, 2, numel(to_refine));
                ls = zeros(numel(to_refine),1);
                lmark_lhoods = zeros(num_points,numel(to_refine));
                views_used = zeros(numel(to_refine),1);
                globals_n = cell(size(to_refine,1),1);
                locals_n = cell(size(to_refine,1),1);
                ind = 1;
                for v = to_refine'
                    [shapes(:,:,ind), globals_n{ind}, locals_n{ind},...
                        ls(ind),lmark_lhoods(:,ind),views_used(ind)] ...
                        = Fitting_from_bb(Image, DepthImage, bounding_box, PDM, patches_large, clmParams_large, 'gparam', globals{v}, 'lparam', locals{v});
                    ind = ind+1;
                end
            end

            % TODO aggregation?
            [final_lhood, v_ind] = max(ls);
            landmark_lhoods = lmark_lhoods(:,v_ind);
            global_params = globals_n{v_ind};
            local_params = locals_n{v_ind};
            shape = shapes(:,:,v_ind);
            view_used = views_used(v_ind);
        end
        shape2D = shape;
    else

        % Find the best orientation
        for v = 1:size(views,1)
            [shapes(:,:,v),globals{v},locals{v},ls(v),lmark_lhoods(:,v),views_used(v)] = Fitting_from_bb(Image, [], bounding_box, PDM, patches, clmParams, 'orientation', views(v,:));
        end

        [final_lhood, v_ind] = max(ls);
        landmark_lhoods = lmark_lhoods(:,v_ind);

        global_params = globals{v_ind};
        local_params = locals{v_ind};

        shape2D = shapes(:,:,v_ind);
        view_used = views_used(v_ind);

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
classdef TestInvertAffineTransform
    %TestInvertAffineTransform

    methods (Static)
        function test_1
            src = [0 1 0; -1 0 1];
            ref = [0 -1 1; 1 0 0];
            dst = cv.invertAffineTransform(src);
            validateattributes(dst, {'numeric'}, {'size',[2 3]});
            assert(all(abs(dst(:) - ref(:)) < 1e-5));
        end

        function test_error_argnum
            try
                cv.invertAffineTransform();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [P, resort] = sort_image_filenames(P)
% :Usage:
% ::
%
%     [P, indices] = sort_image_filenames(P)
%
% Not all image name listing functions return imgs in the correct numbered
% order!
%
% This function resorts a string matrix of image file names by image number,
% in ascending order
% At most, filename can have one number in it, or error is returned
%
% P can be a string matrix of a cell of string matrices
%
% ..
%    Tor Wager, April 2005
% ..

% if iscell(P) % if cell, call this function recursively
%     for i = 1:length(P)
%         fprintf(1,'%3.0f',i);
%         P{i} = sort_image_filenames(P{i});
%     end
%     fprintf(1,'\n')
%     return
% end

if size(P,1) == 1, return, end

nums = [];
for j = 1:size(P,1)

    if iscell(P)
        [dd,ff,ee]=fileparts(P{j, 1});
    else
        % str matrix
        [dd,ff,ee]=fileparts(P(j,:));
    end

    % find the last real number in the name
    tmp = nums_from_text(ff);
    tmp(tmp ~= real(tmp)) = [];
    nums(j) = tmp(end);
end

[n2,resort] = sort(nums);

if any(n2-nums)
    fprintf(1,'Resorting.')

    P = P(resort,:);
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function h = assignColumns(f, colIdx, g)
%ASSIGNCOLUMNS   Extract columns (or rows) of an array-valued CHEBTECH.
%   G = ASSIGNCOLUMNS(F, COLIDX, G) assigns the columns specified by the row
%   vector COLIDX from the FUN F so that F(:, COLIDX) = G. COLIDX need not be
%   increasing in order or unique, but must contain only integers in the range
%   [1, M] (where F has M columns) and satisfy LENGTH(COLIDX) = SIZE(G, 2) or
%   ISEMPTY(G);
%
% See also EXTRACTCOLUMNS, MAT2CELL.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% g is empty - remove columns:
if ( isempty(g) )
    h = f;
    h.coeffs(:,colIdx) = [];
    return
end

% Prolong so that f and g have the same length:
if ( length(f) > length(g) )
    g = prolong(g, length(f));
elseif ( length(g) > length(f) )
    f = prolong(f, length(g));
end

% Assign the columns of h.coeffs:
h = f;
h.coeffs(:, colIdx) = g.coeffs;

% Update ishappy:
h.ishappy = f.ishappy && g.ishappy;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function square_hex_grid_test03 ( )

%*****************************************************************************80
%
%% SQUARE_HEX_GRID_TEST03 tests HEX_GRID_01_N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 17;

  nodes_per_layer_test = [ ...
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 21, 41, 81, 101, 1001, 10001 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_HEX_GRID_TEST03\n' );
  fprintf ( 1, '  For a hexagonal grid of points in the unit square,\n' );
  fprintf ( 1, '  given NODES_PER_LAYER, the number of grid points\n' );
  fprintf ( 1, '  along the first layer,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  HEX_GRID_01_N computes N, the total number of grid\n' );
  fprintf ( 1, '  points in the unit square.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    NODES   LAYERS	     N\n' );
  fprintf ( 1, '      PER\n' );
  fprintf ( 1, '    LAYER\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    nodes_per_layer = nodes_per_layer_test ( test );

    layers = hex_grid_01_layers ( nodes_per_layer );

    n = hex_grid_01_n ( nodes_per_layer );

    fprintf ( 1, '  %6d  %6d  %12d\n', nodes_per_layer, layers, n );

  end

  return
end


text:
function f = flops_solve(a,b,c)
% FLOPS_SOLVE    Flops for matrix left division.
% FLOPS_SOLVE(a,b) returns the number of flops for a\b.
% FLOPS_SOLVE(n,m,c) returns the number of flops for ones(n,m)\ones(m,c).

if nargin == 2
  f = flops_solve(rows(a),cols(a),cols(b));
  return;
end
if a == b
  if a == 1
    % scalar division
    f = c*flops_div;
    return;
  end
  f = flops_chol(a) + 2*flops_solve_tri(a,b,c);
elseif a > b
  % this comes from Ax=b, x = (A'*A)\(A'*b)
  f = flops_mul(b,a,b) + flops_mul(b,a,c) + flops_solve(b,b,c);
else
  % this comes from Ax=b, x = A'*(A*A')\b
  f = flops_mul(a,b,a) + flops_mul(b,a,c) + flops_solve(a,a,c);
end


text:
function test_suite = test_normalizeVector3d
%TESTNORMALIZEVECTOR3D  One-line description here, please.
%
%   output = testNormalizeVector3d(input)
%
%   Example
%   testNormalizeVector3d
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-11-16,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testOrthoVectors(testCase) %#ok<*DEFNU>

vec1    = [3 0 0];
exp1    = [1 0 0];
norm1   = normalizeVector3d(vec1);
testCase.assertEqual(exp1, norm1);

vec2    = [0 3 0];
exp2    = [0 1 0];
norm2   = normalizeVector3d(vec2);
testCase.assertEqual(exp2, norm2);

vec3    = [0 0 3];
exp3    = [0 0 1];
norm3   = normalizeVector3d(vec3);
testCase.assertEqual(exp3, norm3);

function testDiagoVector(testCase)

vec = [6 8 0];
norm = normalizeVector3d(vec);
exp = [3/5 4/5 0];
testCase.assertEqual(exp, norm);

vec = [0 6 8];
norm = normalizeVector3d(vec);
exp = [0 3/5 4/5];
testCase.assertEqual(exp, norm);

function testArray(testCase)

vecs = [2 0 0;0 3 0;0 0 4];
exp = [1 0 0;0 1 0;0 0 1];
norm = normalizeVector3d(vecs);
testCase.assertEqual(exp, norm);

vecs = [2 0 0;0 3 0;0 0 4;3 0 4];
exp = [1 0 0;0 1 0;0 0 1;3/5 0 4/5];
norm = normalizeVector3d(vecs);
testCase.assertEqual(exp, norm);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function P = fdkPreFilter(P, D, Filter)
% This function applies the weighting and filtering steps necessary for FDK
% reconstruction of a conebeam CT sinogram
%
% P: conebeam CT sinogram
% D: rotation radius
% Filter: filter to apply ('Shepp-Logan', 'Cosine', 'Hamming', 'Hann' or
% 'None')
%
% Author: Rene Willemink (Signals and Systems Group, University of Twente)
% Date: 2009-03-11

if nargin<3
    Filter = 'None';
end

% Pre-calculate some things
[U V] = ndgrid(1:size(P,1), 1:size(P,2));
U = U-(size(P,1)-1)/2;
V = V-(size(P,2)-1)/2;

% Calculate filter to use in backprojection (in frequency domain)
% Frequency range of filter
w = (0:(size(P,1)-1))' * 2*pi/size(P,1);
H = w;

% Filter with hamming window
switch Filter
  case 'Shepp-Logan'
    % Sinc window
    wind = [0; sin(w(2:end)/2)./(w(2:end)/2)];
  case 'Cosine'
    wind = cos(w/2);
  case 'Hamming'
    alpha = 0.54;
    wind = alpha + (1-alpha)*cos(w);
  case 'Hann'
    wind = (1 + cos(w))/2;
  case 'None'
    wind = ones(size(w));
end
H = H.*wind;

% Step 1 Weighting
P = P .* repmat(D./sqrt(D^2 + U.^2 + V.^2), [1 1 size(P,3)]);

% Step 2 Filtering
P = ifft(fft(P) .* repmat(H, [1 size(P,2) size(P,3)]), 'symmetric');

% Step 3 Normalizing for the backprojection:
P = double(0.5*1/size(P,3)*P);


text:
function JTotal=calcSpherInvJacob(z,systemType)
%%SPHERINVJACOB Calculate the Jacobian for a 3D Cartesian position with
%          respect to monostatic spherical range, azimuth, and elevation
%          components. This produces derivatives of (x,y,z) with respect to
%          (range,Az,El). The function calcSpherJacob produces derivatives
%          of (range,Az,El) with respect to (x,y,z) in the more general
%          bistatic case.
%
%INPUTS: z The 3XN position vectors in the global spherical coordinate
%          system, each with [range;Az;El] components.
% systemType An optional parameter specifying the axis from which the
%          angles are measured in radians. Possible values are
%          0 (The default if omitted) Azimuth is measured
%            counterclockwise from the x-axis in the x-y plane. Elevation
%            is measured up from the x-y plane (towards the z-axis). This
%            is consistent with common spherical coordinate systems for
%            specifying longitude (azimuth) and geocentric latitude
%            (elevation).
%          1 Azimuth is measured counterclockwise from the z-axis in the
%            z-x plane. Elevation is measured up from the z-x plane
%            (towards the y-axis). This is consistent with some spherical
%            coordinate systems that use the z axis as the boresight
%            direction of the radar.
%          2 This is the same as 0 except instead of being given
%            elevation, one desires the angle away from the z-axis, which
%            is (pi/2-elevation).
%          3 This is the same as 0 except azimuth is measured clockwise
%            from the y-axis in the x-y plane instead of counterclockwise
%            from the x-axis. This coordinate system often arises when
%            given "bearings" in a local East-North-Up coordinate system,
%            where the bearing directions are measured East of North.
%
%OUTPUTS: JTotal A 3X3XN Jacobian matrices where the rows in each matrix
%           are [x;y;z] in that order and the columns take the derivative
%           of the row component with respect to [r,azimuth,elevation] in
%           that order.
%
%This function evaluates analytic expressions for the Jacobian matrix that
%were derived by differentiating standard equations for the spherical
%coordinate systems.
%
%EXAMPLE:
%Here, we verify that the inverse of calcSpherJacob is equal to
%calcSpherInvJacob within reasonable finite precision limits.
% z=[1e3;0.3;0.25];
% systemType=2;
% J=calcSpherInvJacob(z,systemType)
% zC=spher2Cart(z,systemType);
% J1Inv=calcSpherJacob(zC,systemType);
% J1=inv(J1Inv)
%The difference between J and J1 values is small (and the values themselves
%are reasonable in magnitude), so we see that the results of this function
%are consistent with inverting the Jacobian from calcSpherJacob.
%
%June 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

if(nargin<2||isempty(systemType))
   systemType=0;
end

N=size(z,2);

JTotal=zeros(3,3,N);

for curPoint=1:N
    r=z(1,curPoint);
    Az=z(2,curPoint);
    El=z(3,curPoint);

    sinAz=sin(Az);
    cosAz=cos(Az);
    sinEl=sin(El);
    cosEl=cos(El);

    J=zeros(3,3);
    switch(systemType)
        case 0
            %dx/dr
            J(1,1)=cosAz*cosEl;
            %dy/dr
            J(2,1)=cosEl*sinAz;
            %dz/dr
            J(3,1)=sinEl;
            %dx/dAz
            J(1,2)=-r*cosEl*sinAz;
            %dy/dAz
            J(2,2)=r*cosAz*cosEl;
            %dz/dAz
            J(3,2)=0;
            %dx/dEl
            J(1,3)=-r*cosAz*sinEl;
            %dy/dEl
            J(2,3)=-r*sinAz*sinEl;
            %dz/dEl
            J(3,3)=r*cosEl;
        case 1
            %dx/dr
            J(1,1)=cosEl*sinAz;
            %dy/dr
            J(2,1)=sinEl;
            %dz/dr
            J(3,1)=cosAz*cosEl;
            %dx/dAz
            J(1,2)=r*cosAz*cosEl;
            %dy/dAz
            J(2,2)=0;
            %dz/dAz
            J(3,2)=-r*cosEl*sinAz;
            %dx/dEl
            J(1,3)=-r*sinAz*sinEl;
            %dy/dEl
            J(2,3)=r*cosEl;
            %dz/dEl
            J(3,3)=-r*cosAz*sinEl;
        case 2
            %dx/dr
            J(1,1)=cosAz*sinEl;
            %dy/dr
            J(2,1)=sinAz*sinEl;
            %dz/dr
            J(3,1)=cosEl;
            %dx/dAz
            J(1,2)=-r*sinAz*sinEl;
            %dy/dAz
            J(2,2)=r*cosAz*sinEl;
            %dz/dAz
            J(3,2)=0;
            %dx/dEl
            J(1,3)=r*cosAz*cosEl;
            %dy/dEl
            J(2,3)=r*cosEl*sinAz;
            %dz/dEl
            J(3,3)=-r*sinEl;
        case 3
            %dx/dr
            J(1,1)=sinAz*cosEl;
            %dy/dr
            J(2,1)=cosAz*cosEl;
            %dz/dr
            J(3,1)=sinEl;
            %dx/dAz
            J(1,2)=r*cosEl*cosAz;
            %dy/dAz
            J(2,2)=-r*sinAz*cosEl;
            %dz/dAz
            J(3,2)=0;
            %dx/dEl
            J(1,3)=-r*sinAz*sinEl;
            %dy/dEl
            J(2,3)=-r*cosAz*sinEl;
            %dz/dEl
            J(3,3)=r*cosEl;
        otherwise
            error('Invalid system type specified.')
    end
    JTotal(:,:,curPoint)=J;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Vi]=surface_intersect(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,interpMethod)

% function [Xi1,Yi1,Zi1]=surface_intersect(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,interpMethod)
% ------------------------------------------------------------------------
% Determines the intersection points of 3 surfaces. X1, Y1, Z1 need to be a
% monotonic plaid surface.
%
% %%EXAMPLE
%
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 30/08/2012
% ------------------------------------------------------------------------

%% Compute intersection curves of surface set 1 and 2

[Xi_12,Yi_12,Zi_12]=surfacePairIntersect(X1,Y1,Z1,X2,Y2,Z2,interpMethod);

%% Compute intersection curves of surface set 1 and 3

[Xi_13,Yi_13,Zi_13]=surfacePairIntersect(X1,Y1,Z1,X3,Y3,Z3,interpMethod);

%% Computer intersection points for the two curve sets 12 and 13

Xi=[]; Yi=[]; Zi=[]; %Unknown number of intersection points
for i=1:numel(Xi_12)
    for j=1:numel(Xi_13);
        Xi_12=cell2mat(Xi_12(i,1));
        Yi_12=cell2mat(Yi_12(i,1));
        Zi_12=cell2mat(Zi_12(i,1));
        Xc2_i=cell2mat(Xi_13(j,1));
        Yc2_i=cell2mat(Yi_13(j,1));
        Zc2_i=cell2mat(Zi_13(j,1));

        %Finding intersection sites in XY projection of curves
        [xi,yi] = polyxpoly(Xi_12,Yi_12,Xc2_i,Yc2_i);
        Xi=[Xi; xi]; Yi=[Yi; yi];
        if ~isempty(xi)
            [Xi_12, m, n] = unique(Xi_12);
            Zi_12=Zi_12(m);
            zi=interp1(Xi_12,Zi_12,xi,'cubic');
            Zi=[Zi; zi];
        end

    end
end

Vi=[Xi Yi Zi];

end

%% OLD VERSION
%
%
% %% CHECK INPUT
%
% % Check for plaid data.
% x1 = X1(1,:); y1 = Y1(:,1);
% if (size(X1,2)>1 && ~isequal(repmat(x1,size(X1,1),1),X1)) || ...
%         (size(Y1,1)>1 && ~isequal(repmat(y1,1,size(Y1,2)),Y1)),
%     error('MATLAB:interp2:meshgrid',...
%         ['X1 and Y1 must be matrices produced by MESHGRID. Use' ...
%         ' GRIDDATA instead \nof INTERP2 for scattered data.']);
% end
%
% %% STEP 1:
% % L=~isnan(X1) & ~isnan(Y1) & ~isnan(Z1);
% % ZI_12 = griddata(X1(L),Y1(L),Z1(L),X2,Y2,'nearest');
%
% ZI_12 = interp2(X1,Y1,Z1,X2,Y2,'linear');
%
% C=(ZI_12-Z2);
% set(figure,'Visible','off');
% [D,h] = contour(X2,Y2,C,[0 0],'Visible','off','LineColor','none');
%
% Xc1={};Yc1={};
% A1=get(h,'children');
% for i=1:1:length(A1)
%     B=get(A1(i));
%     Xc1_i=B.XData(~isnan(B.XData)); Xc1{i,1}=Xc1_i;
%     Yc1_i=B.YData(~isnan(B.YData)); Yc1{i,1}=Yc1_i;
%     Zc1_i=interp2(X1,Y1,Z1,Xc1_i,Yc1_i,'linear'); Zc1{i,1}=Zc1_i;
% %     plot3(Xc1_i(:), Yc1_i(:),Zc1_i(:),'k-','LineWidth',2);
% end
%
% close gcf;
% %% STEP 2:
%
% % ZI_13 = griddata(X1(L),Y1(L),Z1(L),X3,Y3,'nearest');
% ZI_13 = interp2(X1,Y1,Z1,X3,Y3,'linear');
%
% C=(ZI_13-Z3);
% set(figure,'Visible','off');
% [D,h] = contour(X3,Y3,C,[0 0],'Visible','off','LineColor','none');
%
% Xc2={};Yc2={};
% A2=get(h,'children');
% for i=1:1:length(A2)
%     B=get(A2(i));
%     Xc2_i=B.XData(~isnan(B.XData)); Xc2{i,1}=Xc2_i;
%     Yc2_i=B.YData(~isnan(B.YData)); Yc2{i,1}=Yc2_i;
%     Zc2_i=interp2(X1,Y1,Z1,Xc2_i,Yc2_i,'linear'); Zc2{i,1}=Zc2_i;
% %     plot3(Xc2_i(:), Yc2_i(:),Zc2_i(:),'k-','LineWidth',2);
% end
%
% close gcf;
%
% %% STEP 3:
%
% Xi1=[]; Yi1=[]; Zi1=[];
% Xi2=[]; Zi2=[];
% Yi3=[]; Zi3=[];
% for i=1:length(A1)
%     for j=1:length(A2);
%         Xc1_i=cell2mat(Xc1(i,1));
%         Yc1_i=cell2mat(Yc1(i,1));
%         Zc1_i=cell2mat(Zc1(i,1));
%         Xc2_i=cell2mat(Xc2(j,1));
%         Yc2_i=cell2mat(Yc2(j,1));
%         Zc2_i=cell2mat(Zc2(j,1));
%
%         %Finding intersection sites in XY projection of curves
%         [xi,yi] = polyxpoly(Xc1_i,Yc1_i,Xc2_i,Yc2_i);
%         Xi1=[Xi1; xi]; Yi1=[Yi1; yi];
%         if ~isempty(xi)
%             [Xc1_i, m, n] = unique(Xc1_i);
%             Zc1_i=Zc1_i(m);
%             zi=interp1(Xc1_i,Zc1_i,xi,'cubic');
%             Zi1=[Zi1; zi];
%         end
%     end
% end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function yp = p00_fun ( test, neqn, t, y )

%*****************************************************************************80
%
%% P00_FUN evaluates the right hand side of any test problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the test problem index.
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, the value of the independent variable.
%
%    Input, real Y(NEQN,1), the value of the dependent variables.
%
%    Output, real YP(NEQN,1), the value of the derivative of the
%    dependent variables, as specified by the ODE.
%
  yp = zeros ( neqn, 1 );

  if ( test == 1 )
    yp = p01_fun ( neqn, t, y );
  elseif ( test == 2 )
    yp = p02_fun ( neqn, t, y );
  elseif ( test == 3 )
    yp = p03_fun ( neqn, t, y );
  elseif ( test == 4 )
    yp = p04_fun ( neqn, t, y );
  elseif ( test == 5 )
    yp = p05_fun ( neqn, t, y );
  elseif ( test == 6 )
    yp = p06_fun ( neqn, t, y );
  elseif ( test == 7 )
    yp = p07_fun ( neqn, t, y );
  elseif ( test == 8 )
    yp = p08_fun ( neqn, t, y );
  elseif ( test == 9 )
    yp = p09_fun ( neqn, t, y );
  elseif ( test == 10 )
    yp = p10_fun ( neqn, t, y );
  elseif ( test == 11 )
    yp = p11_fun ( neqn, t, y );
  elseif ( test == 12 )
    yp = p12_fun ( neqn, t, y );
  elseif ( test == 13 )
    yp = p13_fun ( neqn, t, y );
  elseif ( test == 14 )
    yp = p14_fun ( neqn, t, y );
  elseif ( test == 15 )
    yp = p15_fun ( neqn, t, y );
  elseif ( test == 16 )
    yp = p16_fun ( neqn, t, y );
  elseif ( test == 17 )
    yp = p17_fun ( neqn, t, y );
  elseif ( test == 18 )
    yp = p18_fun ( neqn, t, y );
  elseif ( test == 19 )
    yp = p19_fun ( neqn, t, y );
  elseif ( test == 20 )
    yp = p20_fun ( neqn, t, y );
  elseif ( test == 21 )
    yp = p21_fun ( neqn, t, y );
  elseif ( test == 22 )
    yp = p22_fun ( neqn, t, y );
  elseif ( test == 23 )
    yp = p23_fun ( neqn, t, y );
  elseif ( test == 24 )
    yp = p24_fun ( neqn, t, y );
  elseif ( test == 25 )
    yp = p25_fun ( neqn, t, y );
  elseif ( test == 26 )
    yp = p26_fun ( neqn, t, y );
  elseif ( test == 27 )
    yp = p27_fun ( neqn, t, y );
  elseif ( test == 28 )
    yp = p28_fun ( neqn, t, y );
  elseif ( test == 29 )
    yp = p29_fun ( neqn, t, y );
  elseif ( test == 30 )
    yp = p30_fun ( neqn, t, y );
  elseif ( test == 31 )
    yp = p31_fun ( neqn, t, y );
  elseif ( test == 32 )
    yp = p32_fun ( neqn, t, y );
  elseif ( test == 33 )
    yp = p33_fun ( neqn, t, y );
  elseif ( test == 34 )
    yp = p34_fun ( neqn, t, y );
  elseif ( test == 35 )
    yp = p35_fun ( neqn, t, y );
  elseif ( test == 36 )
    yp = p36_fun ( neqn, t, y );
  elseif ( test == 37 )
    yp = p37_fun ( neqn, t, y );
  elseif ( test == 38 )
    yp = p38_fun ( neqn, t, y );
  elseif ( test == 39 )
    yp = p39_fun ( neqn, t, y );
  elseif ( test == 40 )
    yp = p40_fun ( neqn, t, y );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_FUN - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem index TEST = %d\n', test );
    error ( 'P00_FUN - Fatal error!' );
  end

  return
end


text:
function a = gt( x, y )

%Disciplined convex programming information for GT (>):
%   The right-hand side of a less-than constraint must be convex. The
%   left-hand side must be concave. Of course, real constant and affine
%   expressions are both convex and concave and can be used on either
%   side as well.
%
%Disciplined geometric programming information for GT (>):
%   The right-hand side of a less-than constraint must be log-convex---
%   including positive constants, monomials, posynomials, generalized
%   posynomials, and products thereof. The left-hand side must be log-
%   concave---including positive constants, monomials, reciprocals of
%   log-convex expressions, and products thereof.
%
%Note that CVX does not distinguish between strict greater-than (>) and
%greater-than-or-equal (>=) constraints; they are treated identically.
%Feasible interior-point solvers tend to return points which satisfy
%strict inequality, but not all solvers do.

warning( 'CVX:StrictInequalities', cvx_error( 'The use of strict inequalities in CVX is strongly discouraged, because solvers treat them as non-strict inequalities. Please consider using ">=" instead.', [66,75], false, '' ) );
b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '>' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ r, z, c, s ] = zchex ( r, ldr, p, k, l, z, ldz, nz, job )

%*****************************************************************************80
%
%% ZCHEX updates a Cholesky factorization.
%
%  Discussion:
%
%    ZCHEX updates a Cholesky factorization
%
%      A = hermitian(R) * R
%
%    of a positive definite matrix A of order P under diagonal
%    permutations of the form
%
%      E' * A * E
%
%    where E is a permutation matrix.  Specifically, given
%    an upper triangular matrix R and a permutation matrix
%    E (which is specified by K, L, and JOB), ZCHEX determines
%    a unitary matrix U such that
%
%      U * R * E = RR,
%
%    where RR is upper triangular.  At the user's option, the
%    transformation U will be multiplied into the array Z.
%
%    If A = hermitian(X)*X, so that R is the triangular part of the
%    QR factorization of X, then RR is the triangular part of the
%    QR factorization of X * E, that is, X with its columns permuted.
%
%    For a less terse description of what ZCHEX does and how
%    it may be applied, see the LINPACK guide.
%
%    The matrix Q is determined as the product U(L-K)*...*U(1)
%    of plane rotations of the form
%
%      (    C(I)       S(I) )
%      (                    ) ,
%      ( -conj(S(i))  C(I) )
%
%    where C(I) is real, the rows these rotations operate on
%    are described below.
%
%    There are two types of permutations, which are determined
%    by the value of job.
%
%    JOB = 1, right circular shift:
%    The columns are rearranged in the following order.
%
%      1, ..., K-1, L, K, K+1, ..., L-1, L+1, ..., P.
%
%    U is the product of L-K rotations U(I), where U(I)
%    acts in the (L-I,L-I+1)-plane.
%
%    JOB = 2, left circular shift:
%    The columns are rearranged in the following order
%
%      1, ..., K-1, K+1, K+2, ..., L, L, L+1, ..., P.
%
%    U is the product of L-K rotations U(I), where U(I)
%    acts in the (K+I-1,K+I)-plane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%
%  Parameters:
%
%    Input, complex R(LDR,P); the upper triangular factor that is to be updated.
%    Elements below the diagonal are not referenced.
%
%    Input, integer LDR, the leading dimension of R, which is at least P.
%
%    Input, integer P, the order of the matrix.
%
%    Input, integer K, the first column to be permuted.
%
%    Input, integer L, the last column to be permuted.
%    L must be strictly greater than K.
%
%    Input, complex Z(LDZ,NZ), an array of NZ P-vectors into
%    which the transformation U is multiplied.
%
%    Input, integer LDZ, the leading dimension of Z, which must
%    be at least P.
%
%    Input, integer NZ, the number of columns of the matrix Z.
%
%    Input, integer JOB, determines the type of permutation.
%    1, right circular shift.
%    2, left circular shift.
%
%    Output, complex R(LDR,P), the updated factor.
%
%    Output, complex Z(LDZ,NZ), the updated matrix.  Z is not referenced
%    if NZ = 0.
%
%    Output, real C(P), the cosines of the transforming rotations.
%
%    Output, complex S(P), the sines of the transforming rotations.
%
  if ( job == 1 )
%
%  Right circular shift.
%
%  Reorder the columns.
%
    for i = 1 : l
      ii = l - i + 1;
      s(i) = r(ii,l);
    end

    for jj = k : l - 1
      j = l - 1 - jj + k;
      r(1:j,j+1) = r(1:j,j);
      r(j+1,j+1) = 0.0;
    end

    for i = 1 : k - 1
      ii = l - i + 1;
      r(i,k) = s(ii);
    end
%
%  Calculate the rotations.
%
    t = s(1);
    for i = 1 : l - k
      [ c(i), s(i), s(i+1) ] = zrotg ( s(i+1), t );
      t = s(i+1);
    end

    r(k,k) = t;
    for j = k+1 : p
      il = max ( 1, l - j + 1 );
      for ii = il : l - k
        i = l - ii;
        t = c(ii) * r(i,j) + s(ii) * r(i+1,j);
        r(i+1,j) = c(ii) * r(i+1,j) - conj ( s(ii) ) * r(i,j);
        r(i,j) = t;
      end
    end
%
%  If required, apply the transformations to Z.
%
    for j = 1 : nz
      for ii = 1 : l - k
        i = l - ii;
        t = c(ii) * z(i,j) + s(ii) * z(i+1,j);
        z(i+1,j) = c(ii) * z(i+1,j) - conj ( s(ii) ) * z(i,j);
        z(i,j) = t;
      end
    end

  else
%
%  Left circular shift.
%
%  Reorder the columns.
%
    for i = 1 : k
      ii = l - k + i;
      s(ii) = r(i,k);
    end

    for j = k : l - 1
      r(1:j,j) = r(1:j,j+1);
      jj = j - k + 1;
      s(jj) = r(j+1,j+1);
    end

    for i = 1 : k
      ii = l - k + i;
      r(i,l) = s(ii);
    end

    r(k+1:l,l) = 0.0;
%
%  Reduction loop.
%
    for j = k : p
%
%  Apply the rotations.
%
      if ( j ~= k )
        iu = min ( j - 1, l - 1 );
        for i = k : iu
          ii = i - k + 1;
          t = c(ii) * r(i,j) + s(ii) * r(i+1,j);
          r(i+1,j) = c(ii) * r(i+1,j) - conj ( s(ii) ) * r(i,j);
          r(i,j) = t;
        end
      end

      if ( j < l )
        jj = j - k + 1;
        t = s(jj);
        [ c(jj), s(jj), r(j,j) ] = zrotg ( r(j,j), t );
      end

    end
%
%  Apply the rotations to Z.
%
    for j = 1 : nz
      for i = k : l - 1
        ii = i - k + 1;
        t = c(ii) * z(i,j) + s(ii) * z(i+1,j);
        z(i+1,j) = c(ii) * z(i+1,j) - conj ( s(ii) ) * z(i,j);
        z(i,j) = t;
      end
    end

  end

  return
end

author: johannesgerer


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Specify project path
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
root_path = get_project_path();
export_path = fullfile(root_path,'export','sim');
if ~exist(export_path,'dir')
    mkdir(export_path);
end
addpath(export_path);

%%
robot = acrobot();
robot.configureDynamics();

[frame, fric_coef, geom] = get_contact_frame(robot);
robot = addContact(robot,frame,fric_coef,geom);

%% relative degree two outputs:
ya_2 = robot.States.x(4:end);
t = SymVariable('t');
p = SymVariable('p',[2,1]);
tau = (t-p(2))/(p(1)-p(2));
y2 = VirtualConstraint(robot,ya_2,'joints','DesiredType','Bezier','PolyDegree',5,...
    'RelativeDegree',2,'PhaseType','TimeBased',...
    'PhaseVariable',tau,'PhaseParams',p,'Holonomic',true);

robot = addVirtualConstraint(robot,y2);

%%
robot.compile(export_path);
%%
params.ajoints = [ones(1,6)*0.1;zeros(2,6)];
params.pjoints = [0,1];
params.kjoints = [100,20]; %[kp,kd]

t0 = 0;
tf = 10;
eventnames = [];
sim_opts = [];
logger = SimLogger(robot);
x0 = [zeros(6,1);zeros(6,1)];
io_control  = IOFeedback('IO');
tic
robot.simulate(t0, x0, tf, io_control, params, logger, eventnames, sim_opts);
toc

%%
[conGUI] = load_animator(robot, logger.flow.t, logger.flow.states.x)

%%
t = logger.flow.t;
force = logger.flow.inputs.ConstraintWrench.ffoot;
plot_zmp(t,force);

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev