text:
function [x, c, funVal, ValueL]=LogisticC(A, y, z, opts)
%
%%
% Function LogisticC
%      Logistic Loss with the L1 ball constraint
%
%% Problem
%
%  min  f(x,c) = \sum_i - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
%  s.alpha. \|x\|_1 <=z
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  weight_i denotes the weight for the i-th sample
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        Radius of the L1 ball
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x 1
%  c-         The obtained intercept (scalar)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
%% Related papers
%
% [1]  Jun Liu, Jianhui Chen, and Jieping Ye, Large-Scale Sparse Logistic
%      Regression, KDD 2009
%
% [2]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
%% Related functions
%
%  sll_opts, initFactor, pathSolutionLeast,
%  LogisticR, nnLogisticR, nnLogisticC,
%  eplb
%
%%

%% Verify and initialize the parameters
%%
if (nargin <3)
    error('\n Inputs: A, y and z should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z<=0)
    error('\n z should be positive!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% The parameter 'weight' contains the weight for each training sample.
% See the definition of the problem above.
% The summation of the weights for all the samples equals to 1.
if (isfield(opts,'sWeight'))
    sWeight=opts.sWeight;

    if ( length(sWeight)~=2 || sWeight(1) <=0 || sWeight(2) <= 0)
        error('\n Check opts.sWeight, which contains two positive values');
    end

    % we process the weight, so that the summation of the weights for all
    % the samples is 1.

    p_flag=(y==1);                  % the indices of the postive samples
    m1=sum(p_flag) * sWeight(1);    % the total weight for the positive samples
    m2=sum(~p_flag) * sWeight(2);   % the total weight for the positive samples

    weight(p_flag,1)=sWeight(1)/(m1+m2);
    weight(~p_flag,1)=sWeight(2)/(m1+m2);
else
    weight=ones(m,1)/m;             % if not specified, we apply equal weight
end

%% Starting point initialization

% L2 norm regularization
if isfield(opts,'rsL2')
    rsL2=opts.rsL2;
    if (rsL2<0)
        error('\n opts.rsL2 should be nonnegative!');
    end
else
    rsL2=0;
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1); c=0;
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,1);
    end

    if isfield(opts,'c0')
        c=opts.c0;
    else
        c=0;
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);    Ax=Ax./nu;
end

%% The main program

%% Nemirovski's line search

if (opts.lFlag==0)

    lambda0=0;
    % a guess of the root

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m + rsL2; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;

    %% The Armijo Goldstein line search schemes + accelearted gradient descent

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute As=A*s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
            rsL2/2 * s'*s;

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y.* weight) * (1 - prob)
        b= -weighty.*(1-prob);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g=A'*b - sum(b) * mu';  g=g./nu;
        else
            invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
        end

        g=g+ rsL2 * s;  % add the squared L2 norm regularization term

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L;  c= sc- gc/L;

            % projection
            [x, lambda, zf_step]=eplb(v, n, z, lambda0);
            lambda0=lambda;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ c);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
                rsL2/2 * x'*x;

            r_sum= (v'*v + (c-sc)^2) / 2;
            l_sum=fun_x - fun_s - v'* g - (c-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;    ccp=c-cp;
        funVal(iterStep)=fun_x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% adaptive line search

if (opts.lFlag==1)

    lambda0=0;
    % a guess of the root

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m + rsL2; % the intial guess of the Lipschitz continuous gradient

    gamma=1;
    % we shall set the value of gamma = L,
    % and L is appropriate for the starting point

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;

    for iterStep=1:opts.maxIter
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp;  sc=c + beta* ccp;
                As=Ax + beta* (Ax-Axp);
            else
                alpha= (-1+ sqrt(5)) / 2;  beta=0;
                s=x; sc= c;
                As=Ax;
            end

            % aa= - diag(y) * (A * s + sc)
            aa=- y.*(As+ sc);

            % fun_s is the logistic loss at the search point
            bb=max(aa,0);
            fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
                rsL2/2 * s'*s;

            % compute prob=[p_1;p_2;...;p_m]
            prob=1./( 1+ exp(aa) );

            % b= - diag(y.* weight) * (1 - prob)
            b= -weighty.*(1-prob);

            gc=sum(b); % the gradient of c

            % compute g= AT b, the gradient of x
            if (opts.nFlag==0)
                g=A'*b;
            elseif (opts.nFlag==1)
                g=A'*b - sum(b) * mu';  g=g./nu;
            else
                invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
            end

            g=g+ rsL2 * s;  % add the squared L2 norm regularization term

            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L;  cnew= sc- gc/L;

            % projection
            [xnew, lambda, zf_step]=eplb(v, n, z, lambda0);
            lambda0=lambda;

            v=xnew-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A xnew
            if (opts.nFlag==0)
                Axnew=A* xnew;
            elseif (opts.nFlag==1)
                invNu=xnew./nu; mu_invNu=mu * invNu;
                Axnew=A*invNu -repmat(mu_invNu, m, 1);
            else
                Axnew=A*xnew-repmat(mu*xnew, m, 1);     Axnew=Axnew./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Axnew+ cnew);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
                rsL2/2 * xnew'*xnew;

            r_sum= (v'*v + (cnew-sc)^2) / 2;
            l_sum=fun_x - fun_s - v'* g - (cnew-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;
        % store values for L

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;     x=xnew; xxp=x-xp;
        Axp=Ax;   Ax=Axnew;
        cp=c;     c=cnew; ccp=c-cp;
        % update x and Ax with xnew and Axnew

        funVal(iterStep)=fun_x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

author: thomaskuestner


text:
function Ref = RefSelect(Population,k)
% Reference solutions selection by RSEA strategy

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    k      = min(k,length(Population));
    PopObj = Population.objs;
	[FrontNO,MaxFNO] = NDSort(PopObj,k);
    Next = find(FrontNO<=MaxFNO);
    Pmin = min(PopObj,[],1) + 1e-6;
    Pmax = max(PopObj,[],1);
    if Pmax > Pmin
        PopObj = (PopObj-repmat(Pmin,size(PopObj,1),1))./repmat(Pmax-Pmin,size(PopObj,1),1);
    end

    %% Environmental selection
    Choose = LastSelection(PopObj(Next,:),ismember(Next,find(FrontNO<MaxFNO)),ceil(sqrt(k)),k);
    Ref    = Population(Next(Choose));
end

function Choose = LastSelection(PopObj,Choose,div,k)
% Select part of the solutions based on the radar grid

    %% Identify the extreme solutions
	[~,Extreme] = min(sqrt(sum(PopObj.^2,2)).*sqrt(1-(1-pdist2(PopObj,ones(1,size(PopObj,2)),'cosine')).^2),[],1); %Calculate the extreme points based on PBI
    Choose      = Choose | ismember(1:size(PopObj,1),Extreme);

    %% Calculate the convergence of each solution
	Con = sum(PopObj.^1,2).^1;
    Con = Con./max(Con);

    %% Calculate the radar grid of each solution
    [Site,RLoc] = RadarGrid(PopObj,div);
    RDis        = pdist2(RLoc,RLoc);
    RDis(logical(eye(length(RDis)))) = inf;
    CrowdG      = zeros(1,max(Site));
    temp        = tabulate(Site(Choose));
    CrowdG(temp(:,1)) = temp(:,2);

    %% Select k solutions
    while sum(Choose) < k
        % Delete outline solutions
        remainS  = find(~Choose);
        remainG  = unique(Site(remainS));
        bestG    = CrowdG(remainG) == min(CrowdG(remainG));
        current  = remainS(ismember(Site(remainS),remainG(bestG)));
        fitness  = 0.1.*size(PopObj,2).*Con(current) - min(RDis(current,Choose),[],2); % - 0.1.* min(Dis(current,Choose),[],2);
        [~,best] = min(fitness);
        Choose(current(best))       = true;
        CrowdG(Site(current(best))) = CrowdG(Site(current(best))) + 1;
    end
end

function [Site,RLoc] = RadarGrid(P,div)

	[N,M] = size(P);

    %% Calculate the radar coordinate of each solution
    theta     = 0 : 2*pi/M : 2*pi/M*(M-1);
    RLoc(:,1) = sum(P.*repmat(cos(theta),N,1),2)./sum(P,2);
    RLoc(:,2) = sum(P.*repmat(sin(theta),N,1),2)./sum(P,2);
    RLoc      = (RLoc+1)/2;
    YL        = min(RLoc,[],1);                             % Lower bounary of the transferred points
    YU        = max(RLoc,[],1);                             % Upper bounary of the transferred points
    NRLoc     = (RLoc-repmat(YL,N,1))./repmat(YU-YL,N,1);	% Normalized points

    %% Identify the index of grid of each solution
    GLoc            = floor(NRLoc.*div);
    GLoc(GLoc>=div) = div - 1;
    UniqueGLoc      = sortrows(unique(GLoc,'rows'));
    [~,Site]        = ismember(GLoc,UniqueGLoc,'rows');
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef TestANN_MLP
    %TestANN_MLP

    properties (Constant)
        X = [randn(10,3)+1; randn(10,3)-1];
        Y = [ones(10,1); -ones(10,1)];
        YReg = [ones(10,1); -ones(10,1)] + randn(20,1)*0.5;
    end

    methods (Static)
        function test_classification1
            model = cv.ANN_MLP();
            model.LayerSizes = [3,1];  % 1 output node
            model.TermCriteria.maxCount = 100;
            model.TermCriteria.epsilon = 1e-5;
            model.TrainMethod = 'RProp';
            model.BackpropWeightScale = 0.05;
            model.BackpropMomentumScale = 0.05;
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(TestANN_MLP.X, TestANN_MLP.Y);
            assert(model.isTrained());
            for i=1:numel(model.LayerSizes)
                w = model.getWeights(i-1);
                validateattributes(w, {'numeric'}, {});  %TODO: sizes
            end
            assert(isequal(model.getVarCount(), size(TestANN_MLP.X,2)));

            Yhat = model.predict(TestANN_MLP.X);
            validateattributes(Yhat, {'numeric'}, ...
                {'vector', 'numel',numel(TestANN_MLP.Y)});
            Yhat = sign(Yhat);
            acc = nnz(Yhat == TestANN_MLP.Y) / numel(Yhat);
        end

        function test_classification2
            Y = double([TestANN_MLP.Y==-1, TestANN_MLP.Y==+1]);
            model = cv.ANN_MLP();
            model.LayerSizes = [3,2];  % 2 output nodes (1-of-N encoding)
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(TestANN_MLP.X, Y);
            Yhat = model.predict(TestANN_MLP.X);
            validateattributes(Yhat, {'numeric'}, {'size',size(Y)});
            [~,pred] = max(Yhat, [], 2);
            pred = pred*2 - 3;    % [1 2] -> [-1 1]
            acc = nnz(pred == TestANN_MLP.Y) / numel(pred);
        end

        function test_classification3
            % we load data from Neural Network toolbox
            if mexopencv.isOctave() || ~mexopencv.require('nnet')
                error('mexopencv:testskip', 'toolbox');
            end

            load simpleclass_dataset
            X = simpleclassInputs';   % 1000x2
            Y = simpleclassTargets';  % 1000x4 (1-of-N encoded)
            [~,labels] = max(Y, [], 2);
            trainIdx = 1:500;
            testIdx = 501:1000;

            model = cv.ANN_MLP();
            model.LayerSizes = [size(X,2) 10 size(Y,2)];
            model.TrainMethod = 'Backprop';
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(X(trainIdx,:), Y(trainIdx,:));
            Yhat = model.predict(X(testIdx,:));
            validateattributes(Yhat, {'numeric'}, ...
                {'size',[numel(testIdx) size(Y,2)]});

            [~,pred] = max(Yhat, [], 2);
            acc = nnz(labels(testIdx,:) == pred) / numel(testIdx);
        end

        function test_regression1
            model = cv.ANN_MLP();
            model.LayerSizes = [3,1];
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(TestANN_MLP.X, TestANN_MLP.YReg);
            Yhat = model.predict(TestANN_MLP.X);
            validateattributes(Yhat, {'numeric'}, ...
                {'vector', 'real' 'numel',numel(TestANN_MLP.YReg)});
            err = norm(Yhat - TestANN_MLP.YReg);
        end

        function test_regression2
            % we load data from Neural Network toolbox
            if mexopencv.isOctave() || ~mexopencv.require('nnet')
                error('mexopencv:testskip', 'toolbox');
            end

            load simplefit_dataset
            X = simplefitInputs';   % 94x1
            Y = simplefitTargets';  % 94x1
            trainIdx = 1:3:numel(X);
            testIdx = setdiff(1:numel(X), trainIdx);

            model = cv.ANN_MLP();
            model.LayerSizes = [1 5 1];
            model.TrainMethod = 'Backprop';
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(X(trainIdx,:), Y(trainIdx));
            Yhat = model.predict(X(testIdx,:));
            validateattributes(Yhat, {'numeric'}, ...
                {'vector', 'numel',numel(testIdx)});
            err = norm(Yhat - Y(testIdx));
        end

        function test_data_options
            model = cv.ANN_MLP();
            model.LayerSizes = [3,1];

            N = size(TestANN_MLP.X, 1);
            model.train(TestANN_MLP.X, TestANN_MLP.Y, 'Data',{...
                'Layout','Row', 'VarType','NNNN', ...
                'VarIdx',[], 'SampleIdx',[], 'SampleWeights',ones(N,1), ...
                'TrainTestSplitRatio',1/3, 'TrainTestSplitShuffle',true});
        end

        function test_storage
            fname = tempname();
            model = cv.ANN_MLP();
            model.LayerSizes = [3,1];
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(TestANN_MLP.X, TestANN_MLP.Y);

            model.save([fname '.xml']);
            cleanObj = onCleanup(@() delete([fname '.xml']));
            model1 = cv.ANN_MLP();
            model1.load([fname '.xml']);
            %isequal(model, model1)

            model.save([fname '.yaml']);
            cleanObj = onCleanup(@() delete([fname '.yaml']));
            model2 = cv.ANN_MLP();
            model2.load([fname '.yaml']);
            %isequal(model, model2)

            model1.clear();
            model2.clear();
        end

        function test_serialization
            model = cv.ANN_MLP();
            model.LayerSizes = [3,1];
            model.setActivationFunction('Sigmoid', 'Param1',1, 'Param2',1);
            model.train(TestANN_MLP.X, TestANN_MLP.Y);
            strXML = model.save('.xml');
            strYML = model.save('.yml');

            model2 = cv.ANN_MLP();
            model2.load(strXML, 'FromString',true);
            Yhat = model2.predict(TestANN_MLP.X);

            model3 = cv.ANN_MLP();
            model3.load(strYML, 'FromString',true);
            Yhat = model3.predict(TestANN_MLP.X);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [D, tree, found] = addWordnetFolders(D);
%
% Adds to the object names, synonyms from wordnet.
% Uses is a simple algorithm for sense desambiguation.
%
% Dw = addWordnet(D);
% LMdbshowobjects(LMquery(Dw, 'object.name', 'animal', 'word'), HOMEIMAGES)
%
% You need to install wordnet: http://wordnet.princeton.edu/

% addwordnet3
% 0) get word frequencies in getWordnetSense
% 1) get all words an compute co-ocurrences between words
% 2) for each word select sense common with co-ocurrent words (exclude from
% tree the current word!)

tree.parent = []; % node = 0 is the root node
tree.nodewords = {};
tree.scores = [];

found = [];
clear folders;
for i = 1:length(D)
    folders{i} = D(i).annotation.folder;
end
folders = unique(folders);
for n = 1:length(folders)
    disp(sprintf('Folder %d (out of %d)', n, length(folders)))

    [Dq, q]    = LMquery(D, 'folder', folders{n});
    if length(q)>3
        [Dq, tree, f,des] = addWordnet(Dq, tree);
        D(q) = Dq;
        found(q) = f;

        disp(sprintf('   there are %d descriptions, and we found senses for %d.', length(des), sum(found)))
        Nnodes = length(tree.parent);
        leaves = setdiff(1:Nnodes, tree.parent); % these are all the nodes without parents
        disp(sprintf('   there are %d leaves, and there are %d objects with senses.', length(leaves), sum(tree.scores(leaves))))
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
% std_maketrialinfo() - create trial information structure using the
%                       .epoch structure of EEGLAB datasets
%
% Usage:
%   >> STUDY = std_maketrialinfo(STUDY, ALLEEG);
%
% Inputs:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of the EEG datasets included in the STUDY structure
%
% Inputs:
%   STUDY      - EEGLAB STUDY set updated. The fields which is created or
%                updated is STUDY.datasetinfo.trialinfo
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, April 2010

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

function STUDY = std_maketrialinfo(STUDY, ALLEEG);

%% test if .epoch field exist in ALLEEG structure
epochfield = cellfun(@isempty, { ALLEEG.epoch });
if any(epochfield)
    fprintf('Warning: some datasets are continuous and trial information cannot be created\n');
    return;
end;

%% Make trial info
for index = 1:length(ALLEEG)
    eventlat = abs(eeg_point2lat( [ ALLEEG(index).event.latency ], [ ALLEEG(index).event.epoch ], ALLEEG(index).srate, [ALLEEG(index).xmin ALLEEG(index).xmax]));
    events   = ALLEEG(index).event;
    ff = fieldnames(events);
    ff = setdiff(ff, { 'latency' 'urevent' 'epoch' });
    trialinfo = [];

    % process time locking event fields
    % ---------------------------------
    indtle    = find(eventlat < 0.02);
    epochs    = [ events(indtle).epoch ];
    extractepoch = true;
    if length(epochs) ~= ALLEEG(index).trials

        % special case where there are not the same number of time-locking
        % event as there are data epochs
        if length(unique(epochs)) ~= ALLEEG(index).trials
            extractepoch = false;
            disp('std_maketrialinfo: not the same number of time-locking events as trials, trial info ignored');
        else
            % pick one event per epoch
            [tmp tmpind] = unique(epochs(end:-1:1)); % reversing the array ensures the first event gets picked
            tmpind = length(epochs)+1-tmpind;
            indtle = indtle(tmpind);
            if length(indtle) ~= ALLEEG(index).trials
                extractepoch = false;
                disp('std_maketrialinfo: not the same number of time-locking events as trials, trial info ignored');
            end;
        end;
    end;
    if extractepoch
        commands = {};
        for f = 1:length(ff)
            eval( [ 'eventvals = {events(indtle).' ff{f} '};' ]);
            %if isnumeric(eventvals{1})
            %    eventvals = cellfun(@num2str, eventvals, 'uniformoutput', false);
            %end;
            commands = { commands{:} ff{f} eventvals };
        end;
        trialinfo = struct(commands{:});
        STUDY.datasetinfo(index).trialinfo = trialinfo;
    end;

%    % same as above but 10 times slower
%     for e = 1:length(ALLEEG(index).event)
%         if eventlat(e) < 0.0005 % time locking event only
%             epoch = events(e).epoch;
%             for f = 1:length(ff)
%                 fieldval  = getfield(events, {e}, ff{f});
%                 trialinfo = setfield(trialinfo, {epoch}, ff{f}, fieldval);
%             end;
%         end;
%     end;
end;

author: PatternRecognition

repo: OpenBMI


%  Script file: compare.m
%
%  Purpose:
%    To compare binary and formatted I/O operations.
%    This program generates an array of 10,000 random
%    values and writes it to disk both as a binary and
%    as a formatted file.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/22/07    S. J. Chapman        Original code
%
% Define variables:
%   count     -- Number of values read / written
%   fid       -- File id
%   in_array  -- Input array
%   msg       -- Open error message
%   out_array -- Output array
%   status    -- Operation status
%   time      -- Elapsed time in seconds

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate the data array.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
out_array = randn(1,10000);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First, time the binary output operation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reset timer
tic;

% Loop for 10 times
for ii = 1:10

   % Open the binary output file for writing.
   [fid,msg] = fopen('unformatted.dat','w');

   % Write the data
   count = fwrite(fid,out_array,'float64');

   % Close the file
   status = fclose(fid);

end

% Get the average time
time = toc / 10;
fprintf ('Write time for unformatted file = %6.3f\n',time);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, time the formatted output operation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reset timer
tic;

% Loop for 10 times
for ii = 1:10

   % Open the formatted output file for writing.
   [fid,msg] = fopen('formatted.dat','wt');

   % Write the data
   count = fprintf(fid,'%23.15e\n',out_array);

   % Close the file
   status = fclose(fid);
end

% Get the average time
time = toc / 10;
fprintf ('Write time for formatted file =   %6.3f\n',time);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time the binary input operation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reset timer
tic;

% Loop for 10 times
for ii = 1:10

   % Open the binary file for reading.
   [fid,msg] = fopen('unformatted.dat','r');

   % Read the data
   [in_array, count] = fread(fid,Inf,'float64');

   % Close the file
   status = fclose(fid);

end

% Get the average time
time = toc / 10;
fprintf ('Read time for unformatted file =  %6.3f\n',time);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time the formatted input operation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reset timer
tic;

% Loop for 10 times
for ii = 1:10

   % Open the formatted file for reading.
   [fid,msg] = fopen('formatted.dat','rt');

   % Read the data
   [in_array, count] = fscanf(fid,'%f',Inf);

   % Close the file
   status = fclose(fid);

end

% Get the average time
time = toc / 10;
fprintf ('Read time for formatted file =    %6.3f\n',time);


text:
function varargout = panel_anova(varargin)
% PANEL_ANOVA: Creation and management of list of files to apply some batch proccess.
%
% USAGE:  bstPanelNew = panel_anova('CreatePanel')

% @=============================================================================
% 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, 2009

eval(macro_method);
end

%% ===== CREATE PANEL =====
function bstPanelNew = CreatePanel()
    panelName = 'ANOVA';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    % Create list of nodes (tree-table)
    nodelist = panel_nodelist('CreatePanel', 'Anova', 'Files', 'table');
    % Create the BstPanel object that is returned by the function
    % => constructor BstPanel(jHandle, panelName, sControls)
    bstPanelNew = BstPanel(panelName, ...
                           nodelist.jPanel, ...
                           struct('nodelist', nodelist));
end

%% =========================================================================
%  ===== PROCESSING FUNCTIONS ==============================================
%  =========================================================================

author: brainstorm-tools


text:
 function y = ifft_sym(varargin)
%function y = ifft_sym(varargin)
%|
%| matlab 7.0 introduced a 'symmetric' option to ifft to handle
%| spectra that are (circularly) hermitian symmetric (real signal).
%| this glue routine is to provide backward compatibility for matlab 6.5.
%| Caution: v7 ifft with 'symmetric' just uses the first half of the spectrum
%| along whichever dimension is requested.  Here, for pre v7, I just take
%| the real part.  The difference is neglible in the cases where this
%| routine is expected to be used, where the spectrum should be exactly
%| symmetric but has slight asymmetry due to numerical precision.
%| If the spectrum is severely asymmetric, then "real(ifft())" and
%| ifft(..., 'symmetric') will differ substantially.  (But one should
%| not call this routine in such cases.)

if ~nargin, ir_usage, end
if nargin == 1 && streq(varargin{1}, 'test'), ifft_sym_test, return, end

if ir_is_octave || is_pre_v7
	y = ifft(varargin{:});
	if isa(varargin{1}, 'double')
		tol = 1e-11;
	else
		tol = 1e-6;
	end

	y = reale(y, tol, 'prompt');
else
	y = ifft(varargin{:}, 'symmetric');
end

function y = ifft_sym_test
del = 10^5*eps;
format compact
x1 = [4 2+0i*del 8 2-1i*del]
y1 = ifft(x1)
y2 = ifft_sym(x1)
x2 = fft(y2)
y1 - y2

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function n=random_n;
% random unit vector n

while 1
    % random from cube 3*3*3:
    rfc=3*rand(3,1)-1.5;
    rfcl=sqrt(rfc'*rfc); % length
    if (0.5<=rfcl)&&(rfcl<=1)
        break
    end
end

n=rfc/rfcl;

author: Sable

repo: mcbench-benchmarks


text:
function OutputFiles = bst_project_grid(ResultsFiles, iSubjectDest, isInteractive)
% BST_PROJECT_GRID: Project volume-based source files the template (files must have been computed from a template grid).
%
% USAGE:  OutputFiles = bst_project_grid(ResultsFiles, iSubjectDest=[group], isInteractive=1)
%
% INPUT:
%    - ResultsFiles  : Relative path to sources file to project
%    - iSubjectDest  : Index of destination subject (by default: "Group analysis")
%    - isInteractive : If 1, displays questions and dialog boxes
%                      If 0, consider it is running from a process (no user interactions)

% @=============================================================================
% 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, 2016

% ===== INITIALIZATION ======
% Parse inputs
if (nargin < 3) || isempty(isInteractive)
    isInteractive = 1;
end
if (nargin < 2) || isempty(iSubjectDest)
    [sSubjectDest, iSubjectDest] = bst_get('NormalizedSubject');
else
    sSubjectDest = bst_get('Subject', iSubjectDest);
end
OutputFiles = {};
errMsg = [];
isFirstWarning = 1;

% ===== GET DESTINATION SUBJECT =====
% Check for errors
if isempty(sSubjectDest) || isempty(sSubjectDest.Anatomy) || isempty(sSubjectDest.iAnatomy)
    errMsg = 'Cannot get destination subject.';
    ResultsFiles = [];
end
% Get destination MRI
MriFileDest = sSubjectDest.Anatomy(sSubjectDest.iAnatomy).FileName;
sMriDest = in_mri_bst(MriFileDest);
if isempty(sMriDest) || ~isfield(sMriDest, 'NCS') || ((~isfield(sMriDest.NCS, 'R') || isempty(sMriDest.NCS.R)) && (~isfield(sMriDest.NCS, 'y') || isempty(sMriDest.NCS.y)))
    errMsg = 'You should compute the MNI normalization for the destination MRI first.';
    ResultsFiles = [];
end

% ===== PROJECT FILES =====
iUpdatedStudies = [];
% Progress bar
bst_progress('start', 'Project sources', 'Projecting files...', 0, length(ResultsFiles));
% Loop on input source files
for iRes = 1:length(ResultsFiles)
    % === GET SOURCE SUBJECT ===
    % Find file in database
    [sStudySrc, iStudySrc, iFile, DataType] = bst_get('AnyFile', ResultsFiles{iRes});
    if isempty(sStudySrc)
        errMsg = ['File is not resgistered in database: ' ResultsFiles{iRes}];
        break;
    end
    % Get corresponding subject
    [sSubjectSrc, iSubjectSrc] = bst_get('Subject', sStudySrc.BrainStormSubject);
    if isempty(sSubjectSrc) || isempty(sSubjectSrc.Anatomy) || isempty(sSubjectSrc.iAnatomy)
        errMsg = 'Cannot get source subject.';
        break;
    elseif (iSubjectSrc == iSubjectDest)
        errMsg = 'Source and destination anatomies are the same: nothing ot project.';
        break;
    end
    % Get destination MRI
    MriFileSrc = sSubjectSrc.Anatomy(sSubjectSrc.iAnatomy).FileName;
    sMriSrc = in_mri_bst(MriFileSrc);
    if isempty(sMriSrc) || ~isfield(sMriSrc, 'NCS') || ((~isfield(sMriSrc.NCS, 'R') || isempty(sMriSrc.NCS.R)) && (~isfield(sMriSrc.NCS, 'y') || isempty(sMriSrc.NCS.y)))
        errMsg = 'You should compute the MNI normalization for the source MRI first.';
        ResultsFiles = [];
        break;
    end

    % === LOAD SOURCE FILE ===
    switch (DataType)
        case 'timefreq'
            % Load file
            ResultsMat = in_bst_timefreq(ResultsFiles{iRes}, 1);
            ResFile = ResultsFiles{iRes};
            % Check the data type: timefreq must be source/surface based, and no kernel-based file
            if ~strcmpi(ResultsMat.DataType, 'results')
                errMsg = 'Only source results can be projected.';
                break;
            end
            % Get source model type
            if isfield(ResultsMat, 'HeadModelType') && ~isempty(ResultsMat.HeadModelType)
                % Nothing to do
            elseif ~isempty(ResultsMat.DataFile)
                % Load related source file
                [AssociateMat,AssociateFile] = in_bst_results(ResultsMat.DataFile, 0, 'HeadModelType', 'HeadModelFile');
                ResultsMat.HeadModelType = AssociateMat.HeadModelType;
                ResultsMat.HeadModelFile = AssociateMat.HeadModelFile;
            else
                % By default: volume
                ResultsMat.HeadModelType = 'volume';
            end
            % Remove link to parent file
            ResultsMat.DataFile = [];
            ResultsMat.SurfaceFile = [];
        case {'results', 'link'}
            % Load file
            [ResultsMat, ResFile] = in_bst_results(ResultsFiles{iRes}, 1);
            % If it depends on the data file: try to use the comment of the data file
            if ~isempty(ResultsMat.DataFile) && strcmpi(file_gettype(ResultsMat.DataFile), 'data')
                % Find parent file in database
                [sStudy, iStudy, iData] = bst_get('DataFile', ResultsMat.DataFile);
                % If file was found: use its comment
                if ~isempty(sStudy) && ~isempty(sStudy.Data(iData).Comment)
                    ResultsMat.Comment = sStudy.Data(iData).Comment;
                end
            end
            % Unconstrained sources: Make a flat map
            if (ResultsMat.nComponents ~= 1)
                ResultsMat = process_source_flat('Compute', ResultsMat);
            end
            % Remove link to parent file
            ResultsMat.DataFile = [];
            ResultsMat.ChannelFlag = [];
            ResultsMat.GoodChannel = [];
            ResultsMat.SurfaceFile = [];
            ResultsMat.SourceDecompSa = [];
            ResultsMat.SourceDecompVa = [];
    end

    % Check source model type
    if strcmpi(ResultsMat.HeadModelType, 'surface')
        errMsg = 'To process surface-based source maps, use function bst_project_sources.m instead.';
        break;
    elseif strcmpi(ResultsMat.HeadModelType, 'mixed')
        errMsg = 'Mixed head models are not supported yet.';
        break;
    elseif ~isfield(ResultsMat, 'GridLoc') || isempty(ResultsMat.GridLoc)
        errMsg = ['No field "GridLoc" in the file to project.' 10 'To process surface-based source maps, use function bst_project_sources.m instead.'];
        break;
    end
    % Remove link to parent file
    if isfield(ResultsMat, 'DataFile') && ~isempty(ResultsMat.DataFile)
        ResultsMat.DataFile = [];
    end

    % === GET EXISTING GROUP COORDINATES ===
    GridLoc = [];
    SurfaceFile = [];
    isDestGroup = strcmpi(sSubjectDest.Name, bst_get('NormalizedSubjectName'));
    % Load forward model (only when projecting TO group analysis subject)
    if isDestGroup && isfield(ResultsMat, 'HeadModelFile') && ~isempty(ResultsMat.HeadModelFile)
        HeadModelMat = in_bst_headmodel(ResultsMat.HeadModelFile, 0, 'GridOptions');
        % If this file was computed from a template grid
        if ~isempty(HeadModelMat.GridOptions) && isfield(HeadModelMat.GridOptions, 'Method') && strcmpi(HeadModelMat.GridOptions.Method, 'group') && isfield(HeadModelMat.GridOptions, 'FileName') && ~isempty(HeadModelMat.GridOptions.FileName)
            % Get template source grid
            TemplateFile = HeadModelMat.GridOptions.FileName;
            % If this template file still exists
            if ~isempty(bst_get('HeadModelFile', TemplateFile))
                % Load template source grid
                TemplateMat = in_bst_headmodel(TemplateFile, 0, 'GridLoc', 'SurfaceFile');
                % If the file exists and is not empty
                if ~isempty(TemplateMat) && ~isempty(TemplateMat.GridLoc)
                    GridLoc     = TemplateMat.GridLoc;
                    SurfaceFile = TemplateMat.SurfaceFile;
                end
            else
                disp(['BST> Warning: Template grid file does not exist anymore: ' TemplateFile]);
            end
        end
    end

    % === CONVERT COORDINATES ===
    % If there is no group grid: project with the MNI transformation
    if isempty(GridLoc)
        % Convert grid coordinates: Source SCS => MNI
        GridLoc = cs_convert(sMriSrc, 'scs', 'mni', ResultsMat.GridLoc);
        if isempty(ResultsMat.GridLoc)
            errMsg = 'Cannot convert grid locations to MNI coordinates.';
            break;
        end
        % Convert grid coordinates: MNI => Destination SCS
        GridLoc = cs_convert(sMriDest, 'mni', 'scs', GridLoc);
        % This is an error only when projecting TO the group analysis subject
        if isDestGroup
            % Display warning
            disp('BST> Warning: This source file was not computed based on a template grid.');
            if isInteractive && isFirstWarning
                isFirstWarning = 0;
                java_dialog('warning', ['Warning: This source file was not computed based on a template grid.' 10 ...
                    'You will not be able to average sources from multiple subjects.' 10 10 ...
                    'In order to get a common grid of points for all the subjects:' 10 ...
                    ' - Right-click on "Group analysis" > Generate source grid' 10 ...
                    ' - For each subject, compute the forward model with option "Use template grid"'], 'Project sources');
            end
        end
    end
    % If the surface is not defined (or not on the correct subject): get the default cortex in destination folder
    if (isempty(SurfaceFile) || ~strcmpi(bst_fileparts(SurfaceFile), bst_fileparts(sSubjectDest.FileName))) && ~isempty(sSubjectDest.Surface) && ~isempty(sSubjectDest.iCortex)
        SurfaceFile = sSubjectDest.Surface(sSubjectDest.iCortex).FileName;
    end
    % Update grid in the file structure
    ResultsMat.GridLoc     = GridLoc;
    ResultsMat.SurfaceFile = SurfaceFile;

    % === OUTPUT STUDY ===
    % New condition name
    NewCondition = strrep(sStudySrc.Condition{1}, '@raw', '');
    % Get condition
    [sStudyDest, iStudyDest] = bst_get('StudyWithCondition', [sSubjectDest.Name '/' NewCondition]);
    % Create condition if doesnt exist
    if isempty(iStudyDest)
        iStudyDest = db_add_condition(sSubjectDest.Name, NewCondition, 0);
        if isempty(iStudyDest)
            error(['Cannot create condition: "' sSubjectDest.Name  '/' NewCondition '".']);
        end
        sStudyDest = bst_get('Study', iStudyDest);
    end

    % === SAVE NEW RESULTS ===
    % Use source filename as a base
    [fPath, fBase] = bst_fileparts(ResFile);
    % Remove KERNEL tag for saving full source files
    fBase = strrep(fBase, '_KERNEL', '');
    % Build full filename
    newResultsFileFull = bst_fullfile(bst_fileparts(file_fullpath(sStudyDest.FileName)), [fBase, '_', file_standardize(sSubjectSrc.Name), '.mat']);
    newResultsFileFull = file_unique(newResultsFileFull);
    newResultsFile = file_short(newResultsFileFull);
    % Update comment (add source filename)
    ResultsMat.Comment = [sSubjectSrc.Name '/' ResultsMat.Comment];
    % History: project source
    ResultsMat = bst_history('add', ResultsMat, 'project', ['Project source grid: ' sSubjectSrc.Name ' => ' sSubjectDest.Name]);
    % Save new results file
    bst_save(newResultsFileFull, ResultsMat, 'v6');

    % === ADD FILE IN DATABASE ===
    % Create Results/Timefreq structure for database
    switch (DataType)
        case 'timefreq'
            sNewResults = db_template('Timefreq');
            sNewResults.FileName = newResultsFile;
            sNewResults.Comment  = ResultsMat.Comment;
            sNewResults.DataFile = ResultsMat.DataFile;
            sNewResults.DataType = ResultsMat.DataType;
            % Add new result
            sStudyDest.Timefreq(end + 1) = sNewResults;
        case {'results', 'link'}
            sNewResults = db_template('Results');
            sNewResults.FileName = newResultsFile;
            sNewResults.Comment  = ResultsMat.Comment;
            sNewResults.DataFile = ResultsMat.DataFile;
            sNewResults.isLink   = 0;
            sNewResults.HeadModelType = ResultsMat.HeadModelType;
            % Add new result
            sStudyDest.Result(end + 1) = sNewResults;
    end
    % Update study in database
    bst_set('Study', iStudyDest, sStudyDest);
    iUpdatedStudies = [iUpdatedStudies, iStudyDest];
    % Add to list of returned files
    OutputFiles{end+1} = newResultsFile;

    % Increment progress bar
    bst_progress('inc', 1);
end
% Error handling
if ~isempty(errMsg)
    if isInteractive
        bst_error(errMsg, 'Project sources', 0);
        bst_progress('stop');
    else
        bst_report('Error', 'process_project_sources', ResultsFiles, errMsg);
    end
    return;
end

% ===== UDPATE DISPLAY =====
if isInteractive
    bst_progress('stop');
end
if isempty(OutputFiles)
    return;
end
% Update tree display
% panel_protocols('UpdateTree');
% Update node
panel_protocols('UpdateNode', 'Study', unique(iUpdatedStudies));
% Select first output study
panel_protocols('SelectStudyNode', iUpdatedStudies(1));
% Select first output file
panel_protocols('SelectNode', [], OutputFiles{1});
% Save database
db_save();


text:
clc
clear all
bb=8; % block size
K=256; % number of atoms in the dictionary
img =imread('Test_Fig2_Missing.png');

[N,M,dim]=size(img);
img = double(img);
%Compute mask and extracting its patches
Mask = double(~(img(:,:,1)==0));
blkMask=im2col(Mask,[bb,bb],'sliding');  % distinct  sliding
img_yuv = rgb2ycbcr(uint8(img));
img_inpaint_yuv = zeros(size(img_yuv));
% Interpolation CbCr Componet
img_inpaint_yuv(:,:,2) = Interpolation(double(img_yuv(:,:,2)),~Mask);
img_inpaint_yuv(:,:,3) = Interpolation(double(img_yuv(:,:,3)),~Mask);
IMin0 = double(img_yuv(:,:,1));

load Dict
load Coeff

% Creating the output image
imag_Y=ImageRecover(IMin0,Dict,Coeff);
imag_Y=max(min(imag_Y,255),0);

img_inpaint_yuv(:,:,1) = imag_Y;
img_inpaint_rgb = ycbcr2rgb(uint8(img_inpaint_yuv));
imshow(uint8(img_inpaint_rgb))
imwrite(img_inpaint_rgb,strcat('KSVD_Result_','iter_25','.png'),'png')

author: chongyangtao

repo: Color-Image-Inpainting


text:
function peak_out=pss_sss_foe(peak,capbuf,fc,sampling_carrier_twist,tdd_flag)

% Use (only) the PSS and SSS to calculate the frequency offset.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU Affero 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 Affero General Public License for more details.
%
% You should have received a copy of the GNU Affero General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% fc*k_factor is the receiver's actual RX center frequency.
if sampling_carrier_twist==1
    k_factor=(fc-peak.freq)/fc;
    k_factor_vec = (fc-peak.freq+(3:-2:-3).*1e3)./fc;
% else
%     k_factor=1;
else
    k_factor = peak.k_factor;
    k_factor_vec = k_factor.*ones(1, 4);
end

if tdd_flag == 1
    [peak.freq, k_factor_idx] = refine_fo(capbuf, peak.cp_type, peak.n_id_2, peak.freq, 1920000, peak.frame_start, k_factor_vec);
    k_factor = k_factor_vec(k_factor_idx);
end

%%%%%%
% Calculate the frequency offset
%%%%%%
% Find the location of the first SSS where we can take a DFT.
if (strcmpi(peak.cp_type,'normal'))
    if tdd_flag==1
        pss_sss_dist=round((3*(128+9)+1)*k_factor); % TDD
        first_sss_dft_location=peak.frame_start+(1920-128)*k_factor; % TDD
    else
        pss_sss_dist=round((128+9)*k_factor); % FDD
        first_sss_dft_location=peak.frame_start+(960-128-9-128)*k_factor; % FDD
    end
elseif (strcmpi(peak.cp_type,'extended'))
    if tdd_flag==1
        pss_sss_dist=round(3*(128+32)*k_factor); %TDD
        first_sss_dft_location=peak.frame_start+(1920-128)*k_factor; %TDD
    else
        pss_sss_dist=round((128+32)*k_factor); %FDD
        first_sss_dft_location=peak.frame_start+(960-128-32-128)*k_factor; %FDD
    end
else
  error('Check code...');
end
first_sss_dft_location=wrap(first_sss_dft_location,0.5,9600*2+0.5);
if (first_sss_dft_location-9600*k_factor>0.5)
%if (first_sss_dft_location>=9601)
  first_sss_dft_location=first_sss_dft_location-9600*k_factor;
  sn=10;
else
  sn=0;
end
sss_dft_loc_set=first_sss_dft_location:9600*k_factor:length(capbuf)-127-pss_sss_dist-100;
n_sss=length(sss_dft_loc_set);

h_raw_fo_pss=NaN(n_sss,62);
sss_raw_fo=NaN(n_sss,62);
%M1=0;
%M2=0;
sn=(1-(sn/10))*10;
M=0;
for k=1:n_sss
  sn=(1-(sn/10))*10;
  sss_dft_location=round(sss_dft_loc_set(k));

  % Find where to take the pss dft
  pss_dft_location=sss_dft_location+pss_sss_dist;
  %if (pss_dft_location+127>length(capbuf))
  %  break;
  %end

  % Find the PSS and use it to estimate the channel.
  dft_in_pss=fshift(capbuf(pss_dft_location:pss_dft_location+127),-peak.freq,fs_lte/16);
  % TOC
  %disp('toc diabled');
  dft_in_pss=[dft_in_pss(3:end) dft_in_pss(1:2)];
  dft_out_pss=dft(dft_in_pss);
  h_raw_fo_pss(k,:)=[dft_out_pss(end-30:end) dft_out_pss(2:32)];
  h_raw_fo_pss(k,:)=h_raw_fo_pss(k,:).*conj(pss(peak.n_id_2));

  % Smoothening... Basic...
  for t=1:62
    %arm_length=min([6 t-1 62-t]);
    lt=max([1 t-6]);
    rt=min([62 t+6]);
    % Growing matrix...
    %h_sm(k,t)=mean(h_raw_fo_pss(k,t-arm_length:t+arm_length));
    h_sm(k,t)=mean(h_raw_fo_pss(k,lt:rt));
  end
  %h_sm(k,:)=1;

  % Estimate noise power.
  % Growing matrix...
  pss_np(k)=sigpower(h_sm(k,:)-h_raw_fo_pss(k,:));

  %M2=M2+sum(conj(h_sm(k,:)).*h_raw_fo_pss(k,:))/pss_np(k);

  % Calculate the SSS in the frequency domain
  dft_in_sss=fshift(capbuf(sss_dft_location:sss_dft_location+127),-peak.freq,fs_lte/16)*exp(j*pi*-peak.freq/(fs_lte/16/2)*-pss_sss_dist);
  % TOC
  %disp('toc diabled');
  dft_in_sss=[dft_in_sss(3:end) dft_in_sss(1:2)];
  dft_out_sss=dft(dft_in_sss);
  sss_raw_fo(k,:)=[dft_out_sss(end-30:end) dft_out_sss(2:32)];
  sss_raw_fo(k,:)=sss_raw_fo(k,:).*conj(sss(peak.n_id_1,peak.n_id_2,sn));

  %M1=M1+sum(conj(h_sm(k,:)).*sss_raw_fo(k,:))/pss_np(k);
  M=M+sum(conj(sss_raw_fo(k,:)).*h_raw_fo_pss(k,:).*absx2(h_sm(k,:))./(2*absx2(h_sm(k,:))*pss_np(k)+pss_np(k)^2));
end
%f_off_est=peak.freq+angle(conj(M1)*M2)/(2*pi)/(1/(fs_lte/16)*pss_sss_dist);
f_off_est=peak.freq+angle(M)/(2*pi)/(1/(fs_lte/16)*pss_sss_dist);

peak_out=peak;
peak_out.freq_fine=f_off_est;

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function cvx_optval = sigma_max( x )

%SIGMA_MAX   Internal cvx version.

narginchk(1,1);
if ndims( x ) > 2, %#ok
    error( 'lambda_max is not defined for N-D arrays.' );
elseif ~cvx_isaffine( x ),
    error( 'Input must be affine.' );
end

%
% Construct problem
%

[ m, n ] = size( x );
cvx_optval = lambda_max( [ zeros( m, m ), x ; x', zeros( n, n ) ] );

% 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 determ = neumann_determinant ( row_num, col_num )

%*****************************************************************************80
%
%% NEUMANN_DETERMINANT returns the determinant of the NEUMANN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ROW_NUM, COL_NUM, the number of rows and columns in the grid.
%
%    Output, real DETERM, the determinant.
%
  determ = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [AB]=vectorTensorProductArray(a,b)

%%

%Determine multiplication order
if size(a,2)==3 %type B*A
   B=a;
   A=b;
   multType='post';
else %type A*B
    B=b;
    A=a;
    multType='pre';
end

%Perform product
switch multType
    case 'pre' %type A*B
        %  A1_1*B1 + A1_2*B2 + A1_3*B3
        %  A2_1*B1 + A2_2*B2 + A2_3*B3
        %  A3_1*B1 + A3_2*B2 + A3_3*B3
        AB=[A(:,1).*B(:,1)+A(:,4).*B(:,2)+A(:,7).*B(:,3) ...
            A(:,2).*B(:,1)+A(:,5).*B(:,2)+A(:,8).*B(:,3) ...
            A(:,3).*B(:,1)+A(:,6).*B(:,2)+A(:,9).*B(:,3)];
    case 'post' %type B*A
        %  A1_1*B1 + A1_2*B2 + A1_3*B3
        %  A2_1*B1 + A2_2*B2 + A2_3*B3
        %  A3_1*B1 + A3_2*B2 + A3_3*B3
        AB=[A(:,1).*B(:,1)+A(:,2).*B(:,2)+A(:,3).*B(:,3) ...
            A(:,4).*B(:,1)+A(:,5).*B(:,2)+A(:,6).*B(:,3) ...
            A(:,7).*B(:,1)+A(:,8).*B(:,2)+A(:,9).*B(:,3)];
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 scores = score_system(W,system,K,for_linear_backend)

% Compute scores, using W-matrix as trained by train_mclr
% Inputs:
%
%   W: system parameters
%
%   system: function handle that maps parameters to scores
%           (input data is wrapped in handle)
%
%   K: number of classes
%
%   for_linear_backend: if true, will remove linear dependency from scores
%                       by zero-ing last score and removing the zeros.
%
%  Outputs:
%
%    scores: K-by-N,       if for_linear_backend = false.
%            (K-1)-by-N,   if for_linear_backend = true, in which case
%                          the K-th log-likelihood is 0 for every trial.

if nargin==0
    fprintf('do test of "train_system" to also test this function\n');
    return;
end

if nargin<4
    for_linear_backend = false;
end

scores = system(W);
scores = reshape(scores,K,[]);

if for_linear_backend
    scores = bsxfun(@minus,scores,scores(end,:));
    scores(end,:) = [];
end

author: nesl


text:
function animacion

global p1 p2 p3 p4 p5 p6 p7 p1aux p2aux

global TD_gripper
global robot

p1=[[1.2254, 1.1375, 1.3],[0.0000, -0.0027, -1.0000, -0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p1aux=[[1.2254, 1.1375, 1.3],[0.0000, -0.0027, -1.0000, -0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p2=[[1.1254, 1.975, 1.2],[0.0000, -0.0027, -1.0000, -0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p2aux=[[1.1254, 1.975, 1.3],[0.0000, -0.0027, -1.0000, -0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p3=[[2.2500, -0.150, 0.080],[0.0000, 1.0000, 0.0000, 0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p4=[[1.5000, -0.150, 0.080],[0.0000, 1.0000, 0.0000, 0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p5=[[1, -1.5, 1.3],[0.0000, -0.6946, 0.7194, 0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p6=[[0.05, -1.5, 1.3],[0.0000, -0.6946, 0.7194, 0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
p7=[[0.05, -1.5, 1.2],[0.0000, -0.6946, 0.7194, 0.0000], [0, 0, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

robot = load_robot('KUKA','KR120_R3200_PA');
robot.T0 = [1 0 0 0;0 1 0 0;0 0 1 1;0 0 0 1];
robot.graphical.draw_axes=0;
robot.equipment{1} = load_robot('equipment/tables','table_big');
robot.equipment{1}.graphical.color = [75 75 75]./255;
robot.equipment{1}.graphical.draw_axes=0;
robot.equipment{1}.T0 = [1 0 0 -0.5;0 1 0 -0.5; 0 0 1 0;0 0 0 1];
robot.tool = load_robot('equipment/end_tools','pallet_gripper');
robot.tool.graphical.draw_axes=0;
robot.equipment{2} = load_robot('equipment','conveyor_belt_2');
robot.equipment{2}.T0 = [0 -1 0 0;1 0 0 2.4; 0 0 1 0; 0 0 0 1];
robot.equipment{2}.graphical.draw_axes=0;
robot.equipment{2}.graphical.color = [1 80 30]./255;
robot.equipment{3} = load_robot('equipment','conveyor_belt_2');
robot.equipment{3}.graphical.color = [1 80 30]./255;
robot.equipment{3}.graphical.draw_axes=0;
robot.equipment{3}.T0 = [1 0 0 0; 0 1 0 -2.4; 0 0 1 0; 0 0 0 1];
robot.equipment{4} = load_robot('equipment/tables','palletizing');
robot.equipment{4}.T0 = [-1 0 0 2.5; 0 -1 0 0; 0 0 1 0; 0 0 0 1];
robot.equipment{4}.graphical.draw_axes=0;
robot.piece{1} = load_robot('equipment/objects','loaded_pallet');
robot.piece{1}.T0 = [0 -1 0 1;1 0 0 2.2;0 0 1 1.23;0 0 0 1];
robot.piece{1}.graphical.draw_axes=0;
robot.piece{1}.graphical.color = [90 70 10]./255;
robot.piece{1}.graphical.draw_transparent=0;
drawrobot3d(robot);

MoveJ(p1,'vmax','z100',TD_gripper, 'wobj0');
MoveJ(p2,'vmax','fine',TD_gripper, 'wobj0');
simulation_grip_piece;
WaitTime(2);
MoveJ(p2aux,'vmax','fine',TD_gripper, 'wobj0');
MoveJ(p1aux,'vmax','z100',TD_gripper, 'wobj0');
MoveJ(p3,'vmax','fine',TD_gripper, 'wobj0');
WaitTime(2);
simulation_release_piece;
MoveJ(p4,'vmax','fine',TD_gripper, 'wobj0');
WaitTime(20);
MoveJ(p3,'vmax','fine',TD_gripper, 'wobj0');
simulation_grip_piece;
WaitTime(2);
MoveJ(p5,'vmax','z100',TD_gripper, 'wobj0');
MoveJ(p6,'vmax','fine',TD_gripper, 'wobj0');
MoveJ(p7,'vmax','fine',TD_gripper, 'wobj0');
simulation_release_piece;
WaitTime(2);
MoveJ(p5,'vmax','z100',TD_gripper, 'wobj0');
end

repo: ARTE


text:
function [Data, para] = LoadWavMask_Simu_CHiME4(para, step, dataset)
nCh = 6;
wavlist_noisy = unique(my_cat('local/wavlist_tr05_simu'));
wavlist_clean = findFiles(para.local.wavroot_clean, 'wav');

for i=1:length(wavlist_clean)   % build an index of clean files so we can find them quickly by utterance ID
    [~,curr_uttID] = fileparts(wavlist_clean{i});
    words = ExtractWordsFromString_v2(curr_uttID, '_');
    clean_struct.(['U_' words{2}]) = wavlist_clean{i};
end

wavlist_noisy = wavlist_noisy(1:step:end);  % choose a portion of the data
if strcmpi(dataset,'train')     % split to training and cross validation data
    wavlist_noisy(50:50:end) = [];
else
    wavlist_noisy = wavlist_noisy(50:50:end);
end

% Load data file list
for i=1:length(wavlist_noisy)
    PrintProgress(i, length(wavlist_noisy), 500);
    % provide the file name of the noisy speech,
    if para.topology.nChMask>1
        ch_idx = randperm(nCh);
        ch_idx(ch_idx==2) = [];     % we don't use channel 2
        ch_idx = sort(ch_idx(1:para.topology.nChMask));
        clear wav wavfile
        for j=1:length(ch_idx)
            wavfile{j} = [para.local.wavroot_noisy '/' wavlist_noisy{i} 'CH' num2str(ch_idx(j)) '.wav'];
            if j==1 || ~para.local.useFileName
                [wav(:,j),fs] = audioread(wavfile{j});
            end
        end
        wav_noisy = wav(:,1);
    else
        wavfile{1} = [para.local.wavroot_noisy '/' wavlist_noisy{i}];
        [wav, fs] = audioread(wavfile{1});
        wav_noisy = wav;
    end

    [~,curr_uttID] = fileparts(wavlist_noisy{i});
    words = ExtractWordsFromString_v2(curr_uttID, '_');
    clean_uttID = words{2};
    wavfile_clean = clean_struct.(['U_' clean_uttID]);
    wav_clean = audioread(wavfile_clean);

    % Compute the mask based on local SNR estimate

    % The clean signal and noisy signal have different gains, so we need to
    % first roughly normalize the absolute power of clean speech in wav_noisy and
    % wav_clean
    noise_power0 = mean(wav_noisy(1:fs*0.2).^2);
    noisy_power0 = mean(wav_noisy.^2);
    clean_power0 = noisy_power0 - noise_power0;
    clean_power_true = mean(wav_clean.^2);
    scale = sqrt(max(0,clean_power0) / clean_power_true);
    scale = max(scale, 1);
    wav_clean = wav_clean*scale;

    [~,spec_clean] = wav2abs(wav_clean,fs);
    [~,spec_noisy] = wav2abs(wav_noisy,fs);
    nFr = min(size(spec_clean,2), size(spec_noisy,2));
    spec_clean = spec_clean(1:257,1:nFr);
    spec_noisy = spec_noisy(1:257,1:nFr);
    spec_noise = spec_noisy-spec_clean;
    power_clean = abs(spec_clean).^2;
    power_noisy = abs(spec_noisy).^2;
    if 0
        power_noise = power_noisy - power_clean;       % use the noisy-clean as noise estimate. This is not stable as sometimes the noisy and clean are not in the same scale
    else
        power_noise = repmat(mean(abs(spec_noisy(:,1:20)).^2,2), 1, nFr);   % use the first 20 frames as noise estimate
    end
    SNR = 10*log10(power_clean ./ power_noise);
    mask{i} = logical(SNR>5);   % 0dB threshold gives us too many speech TF bins. So 5dB is used as threshold. The threshold may not be critical, as we are only using the mask to initialize the mask subnet.
    if mod(i,100)==0
        subplot(5,1,1:2); imagesc(log([abs(spec_clean).^2; abs(spec_noisy).^2; power_noise]));
        subplot(5,1,3); imagesc(SNR);
        subplot(5,1,4); imagesc(mask{i});
        subplot(5,1,5); plot(wav_noisy); hold on; plot(wav_clean,'r'); hold off;
        pause(.01);
    end

    % make sure that the features generated from the wave have the same
    % length as the mask
    nSampleRequired = DecideWavLen4XFrames(nFr, para.topology.frame_len, para.topology.frame_shift);
    if para.local.useFileName
        for j=1:length(wavfile)
            wavfile{j} = sprintf('%s 0 %2.3f', wavfile{j}, nSampleRequired/fs);
        end
        wavInt{i} = wavfile;
    else
        wavInt{i} = StoreWavInt16(wav(1:nSampleRequired,:))';      % note that the input waveform to the network should be a row vector
    end
end

Data(1).data = wavInt;
Data(2).data = mask;
para.IO.context = [1 1];    % context size (for splicing) of data streams
para.IO.sparse = [0 0];     % is the streams stored in sparse format?
para.IO.DataSyncSet{1} = [];    % do we need to synchronize any two streams?
para.IO.frame_rate = [16000 100];   % frame rate of the streams, i.e. how many frames per second.
para.IO.isTensor = [1 1];   % whether we intend to use the input streams as tensor or just matrix.

if para.local.useFileName
    para.IO.inputFeature = [0 1];
    para.IO.fileReader(1).name = 'wavfile';     % define the configurations used for reading the data from files
    para.IO.fileReader(1).multiArrayFiles = 1;
    para.IO.fileReader(1).array = 1;
    para.IO.fileReader(1).fs = 16000;
    para.IO.fileReader(1).precision = 'int16';  % store waveforms in int16 rather than floating points to save space.
    para.IO.fileReader(2).name = '';
else
    para.IO.inputFeature = [1 1];
end
end

author: singaxiong

repo: SignalGraph


text:
if( ishandle(UCP_struct.fig_hdl) )
    %classifier arguments
    aux_str = get(UCP_struct.recording_name, 'String');

    %Check if reprocessing is needed
    if( ~strcmpi(recording_name, aux_str) )
        bRecordingChanged = true;
    end
    recording_name = aux_str;

    recording_format = cKnownFormats{get(UCP_struct.recording_fmt, 'Value')};
    tmp_path = get(UCP_struct.tmp_path, 'String');
    op_mode = cKnownModesOfOperation{get(UCP_struct.op_mode, 'Value')};

    %Labeling options
    aux_str = cellstr(get(UCP_struct.labelings, 'String'));
    aux_val = get(UCP_struct.labelings, 'value');

    %Check if reprocessing is needed
    if( ~strcmpi(class_labeling, aux_str{aux_val}) )
        bLabelingChanged = true;
    end
    class_labeling = aux_str{aux_val};

    aux_str = cellstr(get(UCP_struct.class_filter, 'string'));
    aux_val = get(UCP_struct.class_filter, 'value');

    class_filter = char(aux_str(aux_val));

    %clustering options
    CantClusters = get(UCP_struct.CantClusters, 'Value');
    iter_times = get(UCP_struct.iter_times, 'Value');
    cluster_presence = get(UCP_struct.cluster_presence, 'Value');

    %User interface options
    CantSamplesClose = str2double(get(UCP_struct.CantSamplesClose, 'String'));
    CantSamplesFar = str2double(get(UCP_struct.CantSamplesFar, 'String'));
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [A,B,C,D] = vibsBlockOperator(omega,U,sigma,u,cL,cT,rhoS,c0,rho0,f)
%+========================================================================+
%|                                                                        |
%|                 OPENVIBS - LIBRARY FOR VIBRO-ACOUSTIC                  |
%|           openVibs is part of the GYPSILAB toolbox for Matlab          |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal, Marc Bakry (c) 2017-2019.                 |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             marc.bakry@polytechnique.edu                               |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : vibsBlockOperator.m                           |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Marc Bakry                  |
%|  ( # )  |   CREATION   : 14.03.2019                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+

% Constants
mu     = rhoS*cT^2;
lambda = rhoS*(cL^2 - 2*cT^2);
w      = 2*pi*f;
k      = w/c0;

% Dimension
n = size(omega.msh.elt,2)-1;

% Green kernel function
if (n == 2)
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[H0(kr)]',k);
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]1',k);
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]2',k);
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]3',k);
    G0          = '[log(r)]';
    gradyG0     = 'grady[log(r)]';
    cteGxy      = 1i/4;
    cteG0       = -1/(2*pi);

elseif (n == 3)
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[exp(ikr)/r]',k);
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',k);
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',k);
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',k);
    G0          = '[1/r]';
    gradyG0     = 'grady[1/r]';
    cteGxy      =  1/(4*pi);
    cteG0       =  1/(4*pi);

else
    error('vibsNeumannBW.m : unavailable case.')
end

% Coupling coeff for Brackage-Werner simulation
beta = 1i*k;

%%%%%%%%%%%%%%%%%%%%%%%%%%% ELASTO-ELASTO (A11) %%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialization
A = cell(n,n);

% Static part
GG = integral(omega,grad(U),grad(U));
for i = 1:n
    for j = 1:n
        % Operator div(U):div(U)
        DD = integral(omega,grad(U,i),grad(U,j));

        % Operator e(U):e(U)
        EE = integral(omega,grad(U,j),grad(U,i));
        if (i==j)
            EE = EE + GG;
        end

        % Summation
        A{i,j} = lambda.*DD + mu.*EE;
    end
end

% Dynamic part
if (f ~= 0)
    Id = integral(omega,U,U);
    for i = 1:n
       A{i,i} = A{i,i} - (rhoS*w^2) .* Id;
    end
end

% Final form (sparse)
A = cell2mat(A);

%%%%%%%%%%%%%%%%%%%%%%%%%% ELASTO-ACOUSTIC (A12) %%%%%%%%%%%%%%%%%%%%%%%%%%

% Coupling to FEM
B = cell(n,1);
for i = 1:n
    % Collocation mass operator
    Id = integral(sigma,ntimes(U,i),u);

    % Collocation boundary operator
    S = cteGxy .* integral(sigma,sigma,ntimes(U,i),Gxy,u) + ...
        cteG0  .* regularize(sigma,sigma,ntimes(U,i),G0,u);

    % Collocation boundary operator
    D = cteGxy .* integral(sigma,sigma,ntimes(U,i),gradyGxy,ntimes(u)) + ...
        cteG0  .* regularize(sigma,sigma,ntimes(U,i),gradyG0,ntimes(u));

    % Final operator Brackage-Werner : [1i*k*beta*S - (Id/2 + D)]
    B{i} = beta.*S - (0.5*Id + D);
end

% Final form (full matrix)
B = cell2mat(B);

%%%%%%%%%%%%%%%%%%%%%%%%%% ELASTO-ACOUSTIC (A21) %%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialization
C = cell(1,n);

% Coupling FEM
for i = 1:n
    C{i} = (rho0*w^2) .* integral(sigma,u,ntimes(U,i));
end

% Final format (sparse)
C = cell2mat(C);

%%%%%%%%%%%%%%%%%%%%%%%%%% ACOUSTIC-ACOUSTIC (A22) %%%%%%%%%%%%%%%%%%%%%%%%

% Finite element mass matrix
Id = integral(sigma,u,u);

% Finite element boundary operator
H  = cteGxy .* (k^2 * integral(sigma,sigma,ntimes(u),Gxy,ntimes(u)) ...
    - integral(sigma,sigma,nxgrad(u),Gxy,nxgrad(u)));
Hr = cteG0  .* (k^2 * regularize(sigma,sigma,ntimes(u),G0,ntimes(u)) ...
    - regularize(sigma,sigma,nxgrad(u),G0,nxgrad(u)));

% Finite element boundary operator
D  = cteGxy .* integral(sigma,sigma,u,gradyGxy,ntimes(u));
Dr = cteG0  .* regularize(sigma,sigma,u,gradyG0,ntimes(u));

% Final operator Brackage-Werner : - [1i*k*beta*(-Id/2 + Dt) - H]
D = - (beta.*(-0.5*Id + (D+Dr).') - (H+Hr));
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function J = J_pnp_loss(image_pt, world_pt, K, R, t)
numPoints = size(image_pt, 1);
J = 0;
for i = 1 : numPoints
    world_point = [world_pt(i, :), 1]; % homogeneous coordinates
    image_point = image_pt(i, :);

    cameraMatrix = [R; t.'] * K;
    projectedPoint = world_point * cameraMatrix;

    if(~isnumeric(R))
        image_point = image_point .* projectedPoint(3);
        projectedPoint = projectedPoint(1 : 2);
    else
        projectedPoint = projectedPoint(1 : 2) / projectedPoint(3);
    end
    d = image_point - projectedPoint;
    J = J + d * d';
end
J = J ./ numPoints;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function qout = vgg_quat_mul(r,q)
% VGG_QUAT_MUL  Quaternion multiplication
%               pq = vgg_quat_mul(p,q)

% awf@robots.ox.ac.uk, 30/01/05

qout = ...
     [(r(1)*q(1) - r(2)*q(2) - r(3)*q(3) - r(4)*q(4))
      (r(1)*q(2) + r(2)*q(1) + r(3)*q(4) - r(4)*q(3))
      (q(3)*r(1) - q(4)*r(2) + q(1)*r(3) + q(2)*r(4))
      (q(4)*r(1) + q(3)*r(2) - q(2)*r(3) + q(1)*r(4)) ];

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function pos = cartPolePosition(z,p)
% pos = cartPolePosition(z,p)
%
% This function computes the position of the cart and the pole, given the
% state of the system
%
% INPUTS:
%   z = [4, n] = [x;q;dx;dq] = state of the system
%   p = parameter struct
%       .g = gravity
%       .m1 = cart mass
%       .m2 = pole mass
%       .l = pendulum length
% OUTPUTS:
%   pos = [4, n] = [x1;y1;x2;y2]; = position of [cart; pole]
%

%%%% unpack the state
x = z(1,:);   %Cart position (Not used in dynamics)
q = z(2,:);   % pendulum (pole) angle, measure from gravity vector

%%%% Unpack the physical parameters
l = p.l;  %Pendulum length

%%%% Position of the cart:
x1 = x;
y1 = zeros(size(x));

%%%% Position of the pole:
x2 = x1 + l*sin(q);
y2 = y1 - l*cos(q);

%%%% Pack up position vector:
pos = [x1;y1;x2;y2];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% This demo shows the working steps for artefact reduction using MARA.
% It requires
%    - the function <fastica.m>
%    - the data fiels which contain the classifier,
%       <fv_training_MARA.mat>  and  <inv_matrix_icbm152.mat>
%       http://www.user.tu-berlin.de/irene.winkler/artifacts/MARAtrdata.zip

eeg_file= fullfile(BTB.DataDir, 'demoMat', 'VPiac_10_10_13', ...
                   'calibration_CenterSpellerMVEP_VPiac');

% Load data
try
  [cnt, mrk, mnt] = file_loadMatlab(file);
catch
  error('You need to run ''demo_convert_ERPSpeller'' first');
end

%High-pass filter the data (ICA otherwise finds drifts)
high_pass_fs = 1; %at 1 Hz in this example
[b, a] = butter(2, high_pass_fs/(cnt.fs/2), 'high');
cnt = proc_filtfilt(cnt, b, a);

%Do FastICA
[icasig, A_ica, W_ica] = fastica(cnt.x'); % 'maxNumIterations', 50);
cnt_ica = cnt; cnt_ica.x = cnt.x *W_ica';

%Classify components (artefact vs. neuro) using MARA
disp('MARA Artefact Classification ...')
[goodcomp, info] = proc_MARA(cnt_ica,mnt.clab,A_ica);
fprintf('MARA identified %d artifactual components. \n', ... ...
        length(A_ica(1,:))- length(goodcomp));
fprintf('Kept %d components \n',length(goodcomp));

%Plot each classified components to check if ICA separation and MARA
%classification was successful
plot_BSScomponents(cnt_ica, mnt, W_ica, A_ica, 'goodcomp', goodcomp, 'out', info.out_p);

%Remove activity to get (hopefully) cleaner EEG signals: Reconstruct EEG
%only with the good components in goodcomp
cnt_clean = cnt;
cnt_clean.x =  cnt.x * W_ica(goodcomp, :)' * A_ica(:, goodcomp)';

author: bbci


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef simulationClass<handle
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The  ``simulationClass`` creates a ``simu`` object saved to the MATLAB
    % workspace. The ``simulationClass`` includes properties and methods used
    % to define WEC-Sim's simulation parameters and settings.
    %
    %.. autoattribute:: objects.simulationClass.simulationClass
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (SetAccess = 'public', GetAccess = 'public') % WEC-Sim input
        adjMassFactor       = 2                                            % (`integer`) Weighting function for adjusting added mass term in the translational direction. Default = ``2``
        b2b                 = 0                                            % (`integer`) Flag for body2body interactions, Options: 0 (off), 1 (on). Default = ``0``
        cicDt               = []                                           % (`float`) Time step to calculate Convolution Integral. Default = ``dt``
        cicEndTime          = 60                                           % (`float`) Convolution integral time. Default = ``60`` s
        domainSize          = 200                                          % (`float`) Size of free surface and seabed. This variable is only used for visualization. Default = ``200`` m
        explorer            = 'on'                                         % (`string`) SimMechanics Explorer 'on' or 'off'. Default = ``'on'``
        dt                  = 0.1                                          % (`float`) Simulation time step. Default = ``0.1`` s
        dtOut               = []                                           % (`float`) Output sampling time. Default = ``dt``
        endTime             = []                                           % (`float`) Simulation end time. Default = ``[]``
        gravity             = 9.81                                         % (`float`) Acceleration due to gravity. Default = ``9.81`` m/s
        mcrMatFile          = []                                           % (`string`) mat file that contain a list of the multiple conditions runs with given conditions. Default = ``[]``
        mcrExcelFile        = []                                           % (`string`) File name from which to load wave statistics data. Default = ``[]``
        mode                = 'normal'                                     % (`string`) Simulation execution mode, 'normal', 'accelerator', 'rapid-accelerator'. Default = ``'normal'``
        morisonDt           = []                                           % (`float`) Sample time to calculate Morison Element forces. Default = ``dt``
        nonlinearDt         = []                                           % (`float`) Sample time to calculate nonlinear forces. Default = ``dt``
        paraview            = struct(...                                   % (`structure`) Defines the Paraview visualization.
            'option',           0,...                                      %
            'startTime',        0, ...                                     %
            'endTime',          100, ...                                   %
            'dt',               0.1, ...                                   %
            'path',             'vtk')                                     % (`structure`) Defines the Paraview visualization. ``option`` (`integer`) Flag for paraview visualization, and writing vtp files, Options: 0 (off) , 1 (on). Default = ``0``. ``startTime`` (`float`) Start time for the vtk file of Paraview. Default = ``0``. ``endTime`` (`float`) End time for the vtk file of Paraview. Default = ``100``.  ``dt`` (`float`) Timestep for Paraview. Default = ``0.1``. ``path`` (`string`) Path of the folder for Paraview vtk files. Default = ``'vtk'``.
        pressure            = 0                                            % (`integer`) Flag to save pressure distribution, Options: 0 (off), 1 (on). Default = ``0``
        rampTime            = 100                                          % (`float`) Ramp time for wave forcing. Default = ``100`` s
        rateTransition      = 'on'                                         % (`string`) Flag for automatically handling rate transition for data transfer, Opyions: 'on', 'off'. Default = ``'on'``
        reloadH5Data        = 0                                            % (`integer`) Flag to re-load hydro data from h5 file between runs, Options: 0 (off), 1 (on). Default = ``0``
        rho                 = 1000                                         % (`float`) Density of water. Default = ``1000`` kg/m^3
        saveStructure       = 0                                            % (`integer`) Flag to save results as a MATLAB structure, Options: 0 (off), 1 (on). Default = ``0``
        saveText            = 0                                            % (`integer`) Flag to save results as ASCII files, Options: 0 (off), 1 (on). Default = ``0``
        saveWorkspace       = 1                                            % (`integer`) FLag to save .mat file for each run, Options: 0 (off), 1 (on). Default = ``1``
        simMechanicsFile    = 'NOT DEFINED'                                % (`string`) Simulink/SimMechanics model file. Default = ``'NOT DEFINED'``
        solver              = 'ode4'                                       % (`string`) PDE solver used by the Simulink/SimMechanics simulation. Any continuous solver in Simulink possible. Recommended to use 'ode4, 'ode45' for WEC-Sim. Default = ``'ode4'``
        stateSpace          = 0                                            % (`integer`) Flag for convolution integral or state-space calculation, Options: 0 (convolution integral), 1 (state-space). Default = ``0``
        startTime           = 0                                            % (`float`) Simulation start time. Default = ``0`` s
        zeroCross           = 'DisableAll'                                 % (`string`) Disable zero cross control. Default = ``'DisableAll'``
        outputDir           = 'output'                                     % (`string`) Data output directory name. Default = ``'output'``
    end

    properties (SetAccess = 'public', GetAccess = 'public') % internal WEC-Sim
        caseDir             = []                                           % (`string`) WEC-Sim case directory. Default = dependent
        numCables           = 0                                            % (`integer`) Number of cables in the wec model. Default = ``0``
        numConstraints      = 0                                            % (`integer`) Number of contraints in the wec model. Default = ``0``
        numDragBodies       = 0                                            % (`integer`) Number of drag bodies that comprise the WEC device (excluding hydrodynamic bodies). Default = ``0``
        numHydroBodies      = 0                                            % (`integer`) Number of hydrodynamic bodies that comprise the WEC device. Default = ``0``
        numMoorings         = 0                                            % (`integer`) Number of moorings in the wec model. Default = ``0``
        numPtos             = 0                                            % (`integer`) Number of power take-off elements in the model. Default = ``0``
        numPtoSim           = 0                                            % (`integer`) Number of PTO-Sim elements in the model. Default = ``0``
        time                = 0                                            % (`float`) Simulation time [s]. Default = ``0`` s
    end

    properties (SetAccess = 'private', GetAccess = 'public') % internal WEC-Sim
        caseFile            = []                                           % (`string`) .mat file with all simulation information. Default = dependent
        cicLength           = []                                           % (`integer`) Number of timesteps in the convolution integral length. Default = dependent
        cicTime             = []                                           % (`float vector`) Convolution integral time series. Default = dependent
        date                = datetime                                     % (`string`) Simulation date and time
        gitCommit           = []                                           % (`string`) GitHub commit
        maxIt               = []                                           % (`integer`) Total number of simulation time steps. Approximate for variable step solvers. Default = dependent
        wsVersion           = '5.0.1'                                        % (`string`) WEC-Sim version
    end

    methods
        function obj = simulationClass()
            % This method initializes the ``simulationClass``.
            fprintf('WEC-Sim: An open-source code for simulating wave energy converters\n')
            fprintf('Version: %s\n\n',obj.wsVersion)
            fprintf('Initializing the Simulation Class...\n')
            obj.caseDir = pwd;
            fprintf('\tCase Dir: %s \n',obj.caseDir)
            obj.outputDir = ['.' filesep obj.outputDir];
        end

        function checkInputs(obj)
            % This method checks WEC-Sim user inputs and generates error messages if parameters are not properly defined.

            % Check that simMechanics file exists
            obj.simMechanicsFile = [obj.caseDir filesep obj.simMechanicsFile];
            if exist(obj.simMechanicsFile,'file') ~= 4
                error('The simMechanics file, %s, does not exist in the case directory',obj.simMechanicsFile)
            end
            % Check 'simu.paraview' fields
            if length(fieldnames(obj.paraview)) ~=5
                error(['Unrecognized method, property, or field for class "simulationClass", ' ...
                    '"simulationClass.paraview" structure must only include fields: "option", "startTime", "endTime", "dt", "outputDir"']);
            end
            % Check that visualization is off when using accelerator modes
            if (strcmp(obj.mode,'accelerator') || strcmp(obj.mode,'rapid-accelerator')) ...
                    && strcmp(obj.explorer,'on')
                warning('Mechanics explorer not allowed in accelerator or rapid-accelerator modes. Turning mechanics explorer off.');
                obj.explorer = 'off';
            end

             % Checks user input to ensure that ``simu.endTime`` is specified
            if isempty(obj.endTime) || obj.endTime < obj.startTime
                error('simu.endTime, the simulation end time must be specified in the wecSimInputFile and must be greater than simu.startTime.')
            end
            % Checks user input to ensure that output timestep is larger
            % than or equal to simulation timestep
            if obj.dtOut < obj.dt
                warning('Simu.dtOut, the output timestep, should be greater than or equal to simu.dt, the simulation timestep. The smaller output timestep will be overwritten with the value of the simulation timestep to avoid errors.')
            end
        end

        function setup(obj)
            % Sets simulation properties based on values specified in input file
            obj.time = obj.startTime:obj.dt:obj.endTime;
            obj.maxIt = floor((obj.endTime - obj.startTime) / obj.dt);
            % Set dtOut if it was not specificed in input file
            if isempty(obj.dtOut) || obj.dtOut < obj.dt
                obj.dtOut = obj.dt;
            end
            % Set nonlinearDt if it was not specificed in input file
            if isempty(obj.nonlinearDt) || obj.nonlinearDt < obj.dt
                obj.nonlinearDt = obj.dt;
            end
            % Set cicDt if it was not specificed in input file
            if isempty(obj.cicDt) || obj.cicDt < obj.dt
                obj.cicDt = obj.dt;
            end
            % Set morisonDt if it was not specificed in input file
            if isempty(obj.morisonDt) || obj.morisonDt < obj.dt
                obj.morisonDt = obj.dt;
            end
            % Set paraview.dt if it was not specified in input file
            if isempty(obj.paraview.dt) || obj.paraview.dt < obj.dtOut
                obj.paraview.dt = obj.dtOut;
            end
            obj.cicTime = 0:obj.cicDt:obj.cicEndTime;
            obj.cicLength = length(obj.cicTime);
            obj.caseFile = [obj.simMechanicsFile(length(obj.caseDir)+1:end-4) '_matlabWorkspace.mat'];
            % Remove existing output folder
            if exist(obj.outputDir,'dir') ~= 0
                try
                    rmdir(obj.outputDir,'s')
                catch
                    warning('The output directory could not be removed. Please close any files in the output directory and try running WEC-Sim again')
                end
            end
            mkdir(obj.outputDir)
            obj.getGitCommit;
        end

        function listInfo(obj,waveTypeNum)
            % Lists simulation info
            fprintf('\nWEC-Sim Simulation Settings:\n');
            %fprintf('\tTime Marching Solver                 = Fourth-Order Runge-Kutta Formula \n')
            fprintf('\tStart Time                     (sec) = %G\n',obj.startTime)
            fprintf('\tEnd Time                       (sec) = %G\n',obj.endTime)
            fprintf('\tTime Step Size                 (sec) = %G\n',obj.dt)
            fprintf('\tRamp Function Time             (sec) = %G\n',obj.rampTime)
            if waveTypeNum > 10
                fprintf('\tConvolution Integral Interval  (sec) = %G\n',obj.cicEndTime)
            end
            fprintf('\t Number of Time Steps     = %u \n',obj.maxIt)
        end

        function obj = loadSimMechModel(obj,fName)
            % This method loads the simulink model and sets parameters
            %
            % Parameters
            % ------------
            %   fname : string
            %       the name of the SimMechanics ``.slx`` file
            %

            load_system(fName);
            [~,modelName,~] = fileparts(fName);
            set_param(modelName,'Solver',obj.solver,...
                'StopTime',num2str(obj.endTime),...
                'SimulationMode',obj.mode,...
                'StartTime',num2str(obj.startTime),...
                'FixedStep',num2str(obj.dt),...
                'MaxStep',num2str(obj.dt),...
                'AutoInsertRateTranBlk',obj.rateTransition,...
                'ZeroCrossControl',obj.zeroCross,...
                'SimCompilerOptimization','on',...
                'ReturnWorkspaceOutputs','off',...
                'SimMechanicsOpenEditorOnUpdate',obj.explorer);
        end

        function rhoDensitySetup(obj,rho,gravity)
            % Assigns density and gravity values
            %
            % Parameters
            % ------------
            %   rho : float
            %       density of the fluid medium (kg/m^3)
            %   gravity : float
            %       gravitational acceleration constant (m/s^2)
            %
            obj.rho = rho;
            obj.gravity   = gravity;
        end

        function getGitCommit(obj)
            % Determines GitHub commit tag
            try
                ws_exe = which('wecSim');
                ws_dir = fileparts(ws_exe);
                git_ver_file = [ws_dir '/../.git/refs/heads/master'];
                obj.gitCommit = textread(git_ver_file,'%s');
            catch
                obj.gitCommit = 'No git commit tag available';
            end
        end
    end
end

author: WEC-Sim

repo: WEC-Sim

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


text:
function backDiffUx = backdiffx(u)
% backward diff in x direction in equation (2.9a)
% \Delta^x_{+}u_{ij} = u_{i+1,j} - u_{ij}

rows = size(u,1);
% u([2:rows rows], :) is u_{i+1,j} in (2.9a)
backDiffUx = u([2:rows rows], :) - u;

author: YimianDai

path: github-repos/MATLAB/YimianDai-Image-Processing-Codes-for-Easier-Understanding/Image-Processing-Codes-for-Easier-Understanding-874302799e48852624bc3760b58b46bd9360f238/src/(Physica D 1992) Nonlinear Total Variation based noise removal algorithms/version_1/support/backdiffx.m


text:
function obj = setName(obj, name)
    % Specifies the name of the NLP function
    %
    % Parameters:
    % name: the name character @type char

    assert(ischar(name), 'NlpFunction:invalidNameStr', ...
        'The name must be a string.');

    % validate name string
    assert(isempty(regexp(name, '\W', 'once')),...
        'NlpFunction:invalidNameStr', ...
        'Invalid name string, can NOT contain special characters.');

    obj.Name = name;

end

author: ayonga

repo: frost-dev

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


function G = rand_inverse_wishart(m, v, H_inv_upper_chol)

% function G = rand_inverse_wishart(m, v, H_inv_upper_chol)
% rand_inverse_wishart  Pseudo random matrices drawn from an
% inverse Wishart distribution
% G = rand_inverse_wishart(m, v, H_inv_upper_chol)
% Returns an m-by-m matrix drawn from an inverse-Wishart distribution.
%
% INPUTS:
%     m:          dimension of G and H_inv_upper_chol.
%     v:          degrees of freedom, greater or equal than m.
%     H_inv_chol: upper cholesky decomposition of the inverse of the
%                 matrix parameter.
%                 The upper cholesky of the inverse is requested here
%                 in order to avoid to recompute it at every random draw.
%                 H_inv_upper_chol = chol(inv(H))
% OUTPUTS:
%     G:          G ~ IW(m, v, H) where H = inv(H_inv_upper_chol'*H_inv_upper_chol)
%                 or, equivalently, using the correspondence between Wishart and
%                 inverse-Wishart: inv(G) ~ W(m, v, S) where
%                 S = H_inv_upper_chol'*H_inv_upper_chol = inv(H)
%
% SPECIAL REQUIREMENT
%     none
%

% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
% Dynare 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.
%
% Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.

X = randn(v, m) * H_inv_upper_chol;

% At this point, X'*X is Wishart distributed
% G = inv(X'*X);

% Rather compute inv(X'*X) using the SVD
[U,S,V] = svd(X, 0);
SSi = 1 ./ (diag(S) .^ 2);
G = (V .* repmat(SSi', m, 1)) * V';