text:
% ZCA & l1-norm operation
function output = whitening_norm(features)
    features = im2double(features);
%     output = features;
    output = zca(features);
    output = l1_norm(output);
%     output = l2_norm(output);
%     output = nuclear(output);
end

% l1-norm
function output = l1_norm(features)
    [h,w,c] = size(features);
    output = zeros(h, w);
    unit = 5;
    % padding
    pad = (unit-1)/2;
    feature_temp = zeros(h+unit-1, w+unit-1, c);
    feature_temp(1+pad:h+pad, 1+pad:w+pad, :) = features;
    for i=1+pad:h+pad
        for j=1+pad:w+pad
            temp = feature_temp(i-pad:i+pad, j-pad:j+pad, :);
            norm = sum(sum(sum(abs(temp),3)))/(unit*unit);
            output(i-pad, j-pad) = norm;
        end
    end
end

% l2-norm
function output = l2_norm(features)
    [h,w,c] = size(features);
    output = zeros(h, w);
    unit = 5;
    % padding
    pad = (unit-1)/2;
    feature_temp = zeros(h+unit-1, w+unit-1, c);
    feature_temp(1+pad:h+pad, 1+pad:w+pad, :) = features;
    for i=1+pad:h+pad
        for j=1+pad:w+pad
            temp = feature_temp(i-pad:i+pad, j-pad:j+pad, :);
            temp = temp.*temp;
            norm = sum(sum(sqrt(sum(temp,3))))/(unit*unit);
            output(i-pad, j-pad) = norm;
        end
    end
end

% nuclear-norm
function output = nuclear(features)
    [h,w,c] = size(features);
    output = zeros(h, w);
    unit = 5;
    % padding
    pad = (unit-1)/2;
    feature_temp = zeros(h+unit-1, w+unit-1, c);
    feature_temp(1+pad:h+pad, 1+pad:w+pad, :) = features;
    for i=1+pad:h+pad
        for j=1+pad:w+pad
            temp = reshape(feature_temp(i-pad:i+pad, j-pad:j+pad, :), [unit*unit c]);
            [U, S, V] = svd(temp, 'econ');
            nu_norm = sum(diag(S));
            output(i-pad, j-pad) = nu_norm;
        end
    end

end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%  Figure 10.54      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_54a.m is a script to generate Fig. 10.54 (a)
% Kp=0.03333, and KI=0.01 corresponds to KpKs=2.0
clf;
[t,x,y,e]=sim('fig10_53and54a');

plot(t,e);
xlabel('Time (sec)');
ylabel('e');
title('Fig. 10.54(a)');
%grid
nicegrid

author: Sable

repo: mcbench-benchmarks


text:
function J = J_func_pTop(q, t, r, b, nv)
len = size(r, 1);
R = q2R(q);
J = 0;
for i = 1 : len
   rr = r(i, :).';
   bb = b(i, :).';
   J = J + 1 / len * (nv(i, :) * (R * rr + t - bb))^2;
end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function P = agm_pi(d)
% AGM_PI  Arithmetic-geometric mean for pi.
% Brent-Salamin algorithm.
% agm_pi(d) produces d decimal digits.
% See Cleve's Corner, "Computing Pi",
% http://www.mathworks.com/company/ ...
%    newsletters/news_notes/2011/

% Copyright 2011 MathWorks, Inc.

digits(d)
a = vpa(1,d);
b = 1/sqrt(vpa(2,d));
s = 1/vpa(4,d);
p = 1;
n = ceil(log2(d));
for k = 1:n
   c = (a+b)/2;
   b = sqrt(a*b);
   s = s - p*(c-a)^2;
   p = 2*p;
   a = c;
end
P = a^2/s;


text:
function [obs,nav]=readrnxobs(obs,nav,opt,fname)

global glc
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading obs file %s...',fname0);

fid=fopen(fname);
% read renix header
[headinfo,fid]=decode_rnxh(fid);

if headinfo.type~='O'
    error('\n The file %s is not renix observation file!!!\n',fname);
end

% decode file header
[headinfo,nav,obs,tobs,fid]=decode_obsh(headinfo,nav,obs,fid);

% decode body
obs=decode_obsb(headinfo,obs,tobs,opt,fid);
if obs.n==0,return;end

% sort obs
obs=sortobs(obs);

fprintf('over\n');

return

repo: GINav


text:
% MAIN  --  Quad-Rotor Simulation
%
% Simulates a quad-rotor using ode45, running a controller that will
% stabilize it to the origin.
%
%

% Dynamics paramters
p.g = 9.81; % (m/s^2) gravity
p.d = 0.3;  % (m) half-width of quad rotor
p.m = 0.2;  % (m) half-mass of the quad rotor

% Controller parameters:
p.wFast = 20;  % (rad/s) - char. freq. of orientation controller
p.wSlowX = 2;  % (rad/s) - char. freq. of horizontal controller
p.wSlowY = 5;  % (rad/s) - char. freq. of vertical controller
p.xi = 1.0;  % (1/1)  -  effective damping ratio in the controller
p.uMax = 5*(p.m*p.g);  % Maximum force available by each rotor

% Initial state and simulation duration
z0 = 2.0*randn(6,1);
tSpan = [0,5];

% Function handles for simulation
ctrlFun = @(z)(  controller(z, p)  );
dynFun = @(t,z)(  dynamics(z, ctrlFun(z), p)  );

% Run the simulation
soln = ode45(dynFun,tSpan,z0);

% Unpack the solution:
t = linspace(tSpan(1), tSpan(2), 150);
z = deval(soln,t);
[u, qRef] = ctrlFun(z);
x = z(1,:);
y = z(2,:);
q = z(3,:);
dx = z(4,:);
dy = z(5,:);
dq = z(6,:);
u1 = u(1,:);
u2 = u(2,:);

% Plot:
figure(1); clf;

subplot(2,2,1); hold on;
plot(tSpan,[0,0],'k--');
plot(t,x);
xlabel('t')
ylabel('x')

subplot(2,2,2); hold on;
plot(tSpan,[0,0],'k--');
plot(t,y);
xlabel('t')
ylabel('y')

subplot(2,2,3); hold on;
plot(t,qRef,'k--');
plot(t,q);
xlabel('t')
ylabel('q')

subplot(2,2,4); hold on;
plot(tSpan,p.uMax*[1,1],'k--');
plot(tSpan,-p.uMax*[1,1],'k--');
plot(t,u1);  plot(t,u2);
xlabel('t')
ylabel('u')
legend('u1','u2');

author: MatthewPeterKelly

repo: OptimTraj


text:
function setFEBioPath(FEBioPathSpec)

% function setFEBioPath(FEBioPathSpec)
% -----------------------------------------------------------------------
% This function sets the FEBio path defined by FEBioPathSpec. Setting the
% path means it is stored in the GIBBON configuration file:
% /GIBBON/config/FEBioPath.txt.
%
% -----------------------------------------------------------------------

%%
filePath=mfilename('fullpath');
toolboxPath=fileparts(fileparts(filePath));
configPath=fullfile(toolboxPath,'config');

fileName=fullfile(configPath,'FEBioPath.txt');

%Import text file containing paths

T{1}='#List paths for febio here, multiple can be given e.g. for use on several platforms. The first valid path found is always used.';

%Set paths in cell
switch class(FEBioPathSpec)
    case 'cell'
        for q=1:1:numel(FEBioPathSpec)
            pathNameNow=FEBioPathSpec{q};
            if exist(pathNameNow,'file')==2
                T{q+1}=FEBioPathSpec{q};
            else
                error('Path specified does not exist')
            end
        end
    otherwise
        if exist(FEBioPathSpec,'file')==2
            T{2}=FEBioPathSpec;
        else
            error('Path specified does not exist')
        end
end

cell2txtfile(fileName,T,0,0); %Write cell to config file

if contains(lower(FEBioPathSpec),'febio2')
        warning('FEBio2 detected. FEBio2 support is depricated. Please upgrade to FEBio3');
end

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function n=nodesize(t,j)
%NODESIZE Node size.
%   S=NODESIZE(T) returns an N-element vector S of the sizes of the
%   nodes in the tree T, where N is the number of nodes.  The size of a
%   node is defined as the number of observations from the data used to
%   create the tree that satisfy the conditions for the node.
%
%   S=NODESIZE(T,J) takes an array J of node numbers and returns the
%   sizes for the specified nodes.
%
%   See also CLASSREGTREE, CLASSREGTREE/NUMNODES.

%   Copyright 2006 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $  $Date: 2006/11/11 22:56:20 $

if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:nodesize:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

if nargin<2
    n = t.nodesize;
else
    n = t.nodesize(j,:);
end

author: Cloud-CV

repo: object-proposals


text:
classdef ChiSquareD
%%CHISQUARED Functions to handle the central or non-central chi-squared
%            distributions. These distributions arise when summing the
%            squares of independent normal random variables. The result is
%            central if the mean of the normal random variables was zero.
%Implemented methods are: mean, var, PDF, CDF, invCDF (only for the central
%                         chi squared distribution), rand, (only for the
%                         central chi squared distribution) entropy
%
%DEPENDENCIES: GammaD.m
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(nu,lambda)
%%MEAN Obtain the mean of chi squared probability distribution.
%
%INPUTS: nu The number of degrees of freedom of the chi-squared
%           distribution. Note that nu>0.
%    lambda The non-centrality parameter of the distribution. In the
%           central chi-squared distribution, this is zero.
%
%OUTPUTS: val The mean of the chi squared distribution under consideration.
%
%The noncentral chi-squared distribution is descirbed in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Noncentral Chi-Squared Distribution." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

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

    val=nu+lambda;
end

function val=var(nu,lambda)
%%VAR Obtain the variance of chi squared probability distribution for the
%     given number of degrees of freedom.
%
%INPUTS: nu The number of degrees of freedom of the chi-square
%           distribution. Note that nu>0.
%    lambda The non-centrality parameter of the distribution. In the
%           central chi-squared distribution, this is zero.
%
%OUTPUTS: val The variance of the chi squared distribution under
%             consideration.
%
%The noncentral chi-squared distribution is descirbed in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Noncentral Chi-Squared Distribution." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

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

    val=2*nu+4*lambda;
end

function val=PDF(x,nu,lambda)
%%PDF Evaluate the chi squared probability distribution function (PDF) at
%     the desired points.
%
%INPUTS: x The point(s) at which the chi-squared PDF is to be evaluated.
%          Note that x>=0 for nonzero probabilities.
%       nu The number of degrees of freedom of the chi-square distribution.
%          Note that nu>0.
%   lambda The non-centrality parameter of the distribution. In the central
%          chi-squared distribution, this is zero.
%
%OUTPUTS: val The value(s) of the chi squared PDF evaluated at x.
%
%The noncentral chi-squared distribution is descirbed in [1].
%
%When lambda is very small but nonzero, the PDF function could run into
%problems due to the logarithm of lambda. The logarithms were used instead
%of a simpler form to lessen the numerical instability. If an invalid value
%for the noncentral distribution is obtained, it is assumed that it is due
%to using an extremely small nonzero value of lambda, and the result is
%approximated using lambda=0.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Noncentral Chi-Squared Distribution." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(lambda))
        lambda=0;
    end

    if(lambda==0)
        val=(1/(2^(nu/2)*gamma(nu/2))).*x.^((nu-2)/2).*exp(-x/2);
    else
        val=0.5*exp(-(x+lambda)/2+(nu/4-0.5)*(log(x)-log(lambda))).*besseli(nu/2-1,sqrt(lambda*x));

        if(~isfinite(val))
            val=(1/(2^(nu/2)*gamma(nu/2))).*x.^((nu-2)/2).*exp(-x/2);
        end
    end

    val(x<0)=0;
end

function prob=CDF(x,nu,lambda)
%%CDF Evaluate the cumulative distribution function (CDF) of the central or
%     noncentral chi-squared distribution at desired points.
%
%INPUTS: x The point(s) at which the chi-squared CDF is to be evaluated.
%          Note that x>0 for nonzero values.
%       nu The number of degrees of freedom of the chi-square distribution.
%          Note that nu>0. When lambda is not zero, the CDF function is
%          numerically stabler when nu is an integer.
%   lambda The non-centrality parameter of the distribution. In the central
%          chi-squared distribution, this is zero.
%
%OUTPUTS: prob The value(s) of the CDF of the chi-squared distribution with
%              nu degrees of freedom and noncentrality parameter lambda
%              evaluated at x.
%
%The CDF of the central chi squared function with nu degrees of freedom
%evaluated at x is just a special case of the incomplete gamma function, as
%described at [1]. The incomplete gamma function is built into Matlab
%without the use of any toolboxes. When the noncentrality parameter is
%nonzero, the CDF is evaluated in terms of the equivalent noncentral gamma
%distribution. The noncentral gamma algorithm is documented further in
%GammaD.m.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Noncentral Chi-Squared Distribution." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%Updated to use the noncentral gamma distribution's function
%December 2014 David A. Karnick, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(lambda))
        lambda=0;
    end

    if(lambda==0)
        prob=gammainc(x/2,nu/2);
    else
        prob=GammaD.CDF(x,nu/2,2,lambda);
    end
    prob(x<=0)=0;
end

function val=invCDF(prob,nu)
%%INVCDF Evaluate the inverse of the cumulative distribution function (CDF)
%        of the central chi squared distribution. This is only implemented
%        for the central chi-squared distribution.
%
%INPUTS: prob The probability or probabilities (0<=prob<=1) at which the
%             argument of the CDF is desired.
%          nu The number of degrees of freedom of the chi-squared
%             distribution. Note that nu>0.
%
%OUTPUTS: val The argument(s) of the CDF that would give the probability or
%             probabilities in prob.
%
%The CDF of the chi squared function with nu degrees of freedom evaluated
%at x is just a special case of the incomplete gamma function, as described
%on Mathworld at [1]. The inverse incomplete gamma function is built into
%Matlab without the use of any toolboxes and thus is called here.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Noncentral Chi-Squared Distribution." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=2*gammaincinv(prob,nu/2);
end

function vals=rand(N,nu,lambda)
%%RAND Generate chi-squared distributed random variables with the given
%      parameters.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%       nu The number of degrees of freedom of the chi-squared
%          distribution. Note that nu>0. If lambda is not zero, nu >=1.
%   lambda The non-centrality parameter of the distribution. In the central
%          chi-squared distribution, this is zero. If this parameter is
%          omitted or an empty matrix is passed, then 0 is used.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated chi-squared random variables.
%
%The algorithm for the central chi squared distribution is an
%implementation of the inverse transform algorithm of Chapter 5.1 of [1].
%When the noncentral distribution is used, the random variables are
%generated by summing the squares of normally distributed random variables.
%
%A non-central chi squared distribution is the sum of a chi-squared
%distribution with nu-1 degrees of freedom plus the squared of a normal
%distribution with mean sqrt(lambda) and standard deviation 1. This is how
%the noncentral distribution is generated and why it is required that
%nu>=1.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, Ed. 4, Amsterdam: Elsevier, 2006.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(lambda))
        lambda=0;
    end

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

    if(lambda==0)
        U=rand(dims);
        vals=ChiSquareD.invCDF(U,nu);
    else
        U=rand(dims);
        vals=ChiSquareD.invCDF(U,nu-1)+(sqrt(lambda)+randn(N)).^2;
    end
end

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

    entropyVal=nu/2+log(2)+gammaln(nu/2)+(1-nu/2)*psi(nu/2);
end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function value = bst_prctile(vector, percentile)
% BST_PRCTILE: Returns the percentile value in vector
%
% USAGE: value = bst_prctile(vector, percentile)

% @=============================================================================
% 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: Martin Cousineau, 2020
%          Raymundo Cassani, 2022

% Try to use toolbox function
if exist('prctile','file')
    value = prctile(vector, percentile);
    return;
end

% Check inputs
if ~isvector(vector)
    error('Only vectors supported.');
end
if any(percentile < 0 | percentile > 100)
    error('Input percentile must be a real value between 0 and 100.');
end

% Custom implementation
vector = sort(vector);
rank   = percentile / 100 * length(vector);
lowerRank = floor(rank + 0.5);
upperRank = lowerRank + 1;
fraction  = rank - lowerRank;
lowerRank(lowerRank < 1) = 1;
upperRank(upperRank > length(vector)) = length(vector);
value = 0.5 * (vector(lowerRank) + vector(upperRank));

if fraction ~= 0
    value = value + fraction * (vector(upperRank) - vector(lowerRank));
end

author: brainstorm-tools


text:
% DEMO  --  Chebyshev Polynomial Multivariate Interpolation
% UPDATED  --  October 28, 2013
% Written by Matthew Kelly, Cornell University
%
% This script demonstrates how to approximate a multivariable function with
% a chebyshev polynomial interpolant. Here we use it to approximate the
% position and velocity of a simple harmonic oscillator.
%

%What order chebyshev polynomial to use?
order = 25;

%set up the time span for the system integration using ode45
tspan = [0,1];
nTime = 1001;

%Initial conditions for the simple harmonic oscillator
z0 = [0.25;0];   %[Position, Velocity]

%Physical parameters for the simple harmonic osciallator
P.m = 0.1;
P.c = 0.4;
P.k = 100;

%Function handle for the simple harmonic oscillator
userFunc = @(t,z)simpleHarmonicOscillator(z,P);

%Use ode45 to get the solution
Options = odeset(...
    'RelTol',1e-10,...   %ode45 complains if reltol is much smaller
    'AbsTol',1e-10);
sol = ode45(userFunc,tspan,z0,Options);

%Get the "exact" solution at each point in tspan
time = linspace(tspan(1),tspan(2),nTime);
z = deval(sol,time);
pos = z(1,:);
vel = z(2,:);

%Plot ode45 stuff
figure(404); clf
subplot(2,1,1); hold on
    plot(time,pos,'r-','LineWidth',2)
    title('Position')
    xlabel('Time')
    ylabel('Position')
subplot(2,1,2); hold on
    plot(time,vel,'r-','LineWidth',2)
    title('Velocity')
    xlabel('Time')
    ylabel('Velocity')

%% Chebyshev Stuff

IO.input = time;
IO.output = z;

%Check the grid points
[f,d,x] = chebyshevFit(IO, order);
subplot(2,1,1);
    plot(x,f(1,:),'ko','MarkerSize',7,'LineWidth',1)
subplot(2,1,2);
    plot(x,f(2,:),'ko','MarkerSize',7,'LineWidth',1)

%Check the interpolation
tic
y = chebyshevInterpolate(f,time,d);
toc
posCheck = y(1,:);
velCheck = y(2,:);

%Plots for the chebyshev interpolant
subplot(2,1,1);
    plot(time,posCheck,'k:','LineWidth',1)
    legend('ode45','chebPts','cheb-pos')
subplot(2,1,2);
    plot(time,velCheck,'k:','LineWidth',1)
legend('ode45','chebPts','cheb-vel')

repo: dscTutorials


text:
classdef CLFQP < Controller
    % This class defines a class of Control Lyapunov Functions (CLFs) that
    % use quadratic programming (QP)
    %
    % @todo implement the CLF-QP controller
    %
    % @author ayonga @date 2016-10-14
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    properties

    end

    methods

        function obj = CLFQP(name)
            % The controller class constructor function
            %
            % Parameters:
            % name: the controller name @type char

            % call superclass constructor
            obj = obj@Controller(name);
            warning('P matrix from the Ricatti equation not verified for relative degree one outputs and not efficiently computed');
%             error('This class has not been completely defined yet.');

        end

        function u = calcControl(obj, t, x, vfc, gfc, plant, params, logger)
            % Computes the control Lyapunov function based control laws
            % control law for virtual constraints
            %
            % Parameters:
            % t: the time instant @type double
            % x: the states @type colvec
            % vfc: the vector field f(x) @type colvec
            % gfc: the vector field g(x) @type colvec
            % plant: the continuous domain @type DynamicalSystem
            % params: the control parameters @type struct
            % logger: the data logger object @type SimLogger
            %
            % Return values:
            % u: the computed torque @type colvec

            nx = plant.numState;
            if strcmp(plant.Type,'SecondOrder')
                q = x(1:nx);
                dq = x(nx+1:end);
            else
                q = x;
                dq = []; % will not affect any computation
            end

            y = struct2array(plant.VirtualConstraints);
            ny = length(y);
            y_a = cell(ny,1);
            y_d = cell(ny,1);
            tau = cell(ny,1);

            % total dimension of the virtual constraints
            dim_y = sum([y.Dimension]);
            % total dimension of outputs (including relative degrees)
            dim_eta = sum([y.Dimension.*y.RelativeDegree]);
            % some constants required for CLF
            F_mat = zeros(dim_eta);
            G_mat = zeros(dim_eta,dim_y);
            I_mat = eye(dim_eta);
            % partial derivative of the highest order of derivative (y^n-1) w.r.t.
            % the state variable 'x'
            DLfy = zeros(dim_y,length(x));   % A = DLfy*gfc; Lf = DLfy*vfc;
            ddy = zeros(dim_y,1);
            mu = zeros(dim_y,1);    % The derivatives mu = k(1) y + k(2) y' + ... k(n) y^(n-1)
            eta = zeros(dim_eta,1);
            idx = 1; % indexing of outputs
            etaidx = 1; % indexing of eta
            for i=1:ny
                y_i = y(i);

                % returns y, y', y'', y^(n-1), Jy^n-1

                % calculate the actual outputs
                %         [y_a{i}{:}] = calcActual(y_i,q,dq);
                offset_param = y_i.OffsetParamName;
                if y_i.hasOffset
                    offset = params.(offset_param);
                    y_a{i} = calcActual(y_i, q, dq, offset);
                else
                    y_a{i} = calcActual(y_i, q, dq);
                end
                % extract the parameter values

                output_param = y_i.OutputParamName; % desired output parameters
                phase_param  = y_i.PhaseParamName;  % phase variable parameters

                if isfield(params,output_param)
                    a = params.(output_param);
                else
                    error('The parameter %s has not been specified in the ''params'' argument.\n', output_param);
                end

                if ~isempty(phase_param)
                    if isfield(params,phase_param)
                        p = params.(phase_param);
                    else
                        error('The parameter %s has not been specified in the ''params'' argument.\n', phase_param);
                    end
                else
                    p = [];
                end
                % calculate the desired outputs
                y_d{i} = calcDesired(y_i, t, q, dq, a, p);
                % calculate the phase variable
                tau{i} = calcPhaseVariable(y_i, t, q, dq, p);

                % control gain (k0,k1,...kN-1) for the feedback term
                if isfield(params, 'epsilon')
                    ep = params.epsilon;
                    K = ones(1, y_i.RelativeDegree);
                    for l= 1:y_i.RelativeDegree
                        K(l) = nchoosek(y_i.RelativeDegree,l-1)*ep^(y_i.RelativeDegree - l + 1);
                    end

                else
                    error('The control gain %s has not been specified in the ''params'' argument.\n', control_param);
                end

                % stack the partial derivatives of all outputs
                y_indices = idx:idx+y_i.Dimension-1;
                eta_indices = etaidx:etaidx+y_i.Dimension-1;

                G_mat(eta_indices+(y_i.RelativeDegree-1)*numel(eta_indices),y_indices) = eye(y_i.Dimension);

                if y_i.RelativeDegree > 1 % only modify the first degree outputs (0th derivative). not higher derivatives
                    I_mat(eta_indices,eta_indices) = ep*eye(numel(eta_indices));
                end

                if strcmp(y_i.PhaseType, 'TimeBased')
                    DLfy(y_indices,:) = y_a{i}{end};
                    ddy(y_indices) = y_d{i}{end};
                else
                    DLfy(y_indices,:) = y_a{i}{end} - y_d{i}{end};
                end
                for j=1:y_i.RelativeDegree
                    mu(y_indices) = mu(y_indices) + K(j)*(y_a{i}{j}-y_d{i}{j});
                    eta(eta_indices+(j-1)*numel(eta_indices)) = y_a{i}{j} - y_d{i}{j};

                    if j < y_i.RelativeDegree
                        F_mat(eta_indices+(j-1)*numel(eta_indices),eta_indices+j*numel(eta_indices)) = eye(numel(eta_indices));
                    end
                end
                % update the starting index for the next output
                idx = idx+y_i.Dimension;
                etaidx = etaidx+y_i.Dimension*y_i.RelativeDegree;
            end

            %% here is where the CLF based controller is designed
            % decoupling matrix
            A_mat  = DLfy*gfc;
            % feedforward term
            Lf_mat = DLfy*vfc;

            % Lyapunov function
            P = care(F_mat,G_mat,eye(dim_eta));
            Pep = I_mat' * P * I_mat;
            Veta = eta' * Pep * eta;

%             % CLF V for mu
%             LgVetamu = 2*eta'*Pep*G_mat;
%             LfVetamu = eta'*(F_mat'*Pep + Pep*F_mat)*eta + 0.3660*ep*Veta;

            % CLF V for u
            LgVetau = 2*eta'*Pep*G_mat*A_mat;
            LfVetau = eta'*(F_mat'*Pep+Pep*F_mat)*eta+0.3660*ep*Veta+2*eta'*Pep*G_mat*Lf_mat-2*eta'*Pep*G_mat*ddy;

            Hmat = A_mat' * A_mat;

            % feedforward controller
            if strcmp(y_i.PhaseType, 'TimeBased')
                bmat = Lf_mat'*A_mat - ddy'*A_mat;
%                 u_ff = - A_mat \ (Lf_mat - ddy);
            else
                bmat = Lf_mat'*A_mat;
%                 u_ff = - A_mat \ Lf_mat;
            end

            options = optimoptions('quadprog','Algorithm','interior-point-convex','display','off');

%             muqp = quadprog(eye(6),[],LgVetamu,-LfVetamu,[],[],[],[],[],options);
%             disp(max(mu - muqp))
%             mu = muqp;
%             u_fb = A_mat \ mu;
%             u = u_ff + u_fb;

            u = quadprog(Hmat,bmat,LgVetau,-LfVetau,[],[],[],[],[],options);
%             disp(max(u - uqp))
            % feedback controller

            if ~isempty(logger)
                calc = logger.calc;

                calc.mu = mu;

                for i=1:ny
                    y_i = y(i);
                    output_name = y_i.Name;

                    for j=1:y_i.RelativeDegree

                        if j > 1
                            ya_name = ['d' num2str(j-1) 'ya_' output_name];
                            yd_name = ['d' num2str(j-1) 'yd_' output_name];
                            tau_name = ['d' num2str(j-1) 'tau_' output_name];
                        else
                            ya_name = ['ya_' output_name];
                            yd_name = ['yd_' output_name];
                            tau_name = ['tau_' output_name];
                        end
                        calc.(ya_name) = y_a{i}{j};
                        calc.(yd_name) = y_d{i}{j};
                        calc.(tau_name) = tau{i}{j};

                    end

                end
%                 calc.u_ff = u_ff;
%                 calc.u_fb = u_fb;
                calc.u = u;

                logger.calc = calc;
            end

        end
    end

end

author: ayonga

repo: frost-dev

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


text:
% Compute the closeness centrality for every vertex: 1/sum(dist to all other nodes)
%
% INPUTs: adjacency matrix, nxn
% OUTPUTs: vector of closeness centralities, nx1
%
% Source: social networks literature (example: Wasserman, Faust, "Social Networks Analysis")
% Other routines used: simpleDijkstra.m
% GB: last updated, Sep 28, 2012

function C=closeness(adj)

C=zeros(length(adj),1);  % initialize closeness vector

for i=1:length(adj); C(i)=1/sum( simpleDijkstra(adj,i) ); end

author: aeolianine

repo: octave-networks-toolbox


text:
%by Davide Di Gloria
%
%Find displacement between 2 grayscale images by using cross-correlation
%
%[dr,dc]=findoff(unreg,ref)
%

function [yoffset,xoffset]=findoff(unreg , base)

c = normxcorr2(unreg,base);

[max_c, imax] = max(abs(c(:)));
[ypeak, xpeak] = ind2sub(size(c),imax(1));

corr_offset = round([(xpeak-(size(c,2)+1)/2) (ypeak-(size(c,1)+1)/2)]);

offset = corr_offset;
xoffset = offset(1);
yoffset = offset(2);


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

%sum(R .* (((X * Theta') .- Y) .^ 2)) 1*4 matrix
J = 1 / 2 * sum(sum(R .* (((X * Theta') .- Y) .^ 2)));
J = J + lambda / 2 * sum(sum(X .^ 2)) + lambda / 2 * sum(sum(Theta .^ 2));

% (((X * Theta') .- Y) * Theta)

X_grad = (R .* ((X * Theta') .- Y)) * Theta;
X_grad = X_grad + lambda .* X;
Theta_grad = (R .* ((X * Theta') .- Y))' * X;
Theta_grad = Theta_grad + lambda .* Theta;

% =============================================================

grad = [X_grad(:); Theta_grad(:)];

end

author: scruel

repo: Notes-ML-AndrewNg


text:
% Compares joint space vs. task space trajectories
% Copyright 2019 The MathWorks, Inc.

%% Setup
clc
createWaypointData;
figure, hold on
plot3(waypoints(1,:),waypoints(2,:),waypoints(3,:),'ko:','LineWidth',2);
title('Trajectory Waypoints');
xlabel('X [m]');
ylabel('Y [m]');
zlabel('Z [m]');
grid on
view([45 45]);
% Define IK solver
ik = inverseKinematics('RigidBodyTree',gen3);
ikWeights = [1 1 1 1 1 1];
% Use a small sample time for this example, so the difference between joint
% and task space is clear due to evaluation of IK in task space trajectory.
ts = 0.02;
trajTimes = 0:ts:waypointTimes(end);
% Initialize matrices for plots
qTask = zeros(numJoints,numel(trajTimes)); % Derived joint angles in task space trajectory
posJoint = zeros(3,numel(trajTimes)); % Derived end effector positions in joint space trajectory

%% Create and evaluate a task space trajectory
ikInitGuess = jointAnglesHome';
ikInitGuess(ikInitGuess > pi) = ikInitGuess(ikInitGuess > pi) - 2*pi;
ikInitGuess(ikInitGuess < -pi) = ikInitGuess(ikInitGuess < -pi) + 2*pi;

disp('Running task space trajectory generation and evaluation...')
tic;

% Trajectory generation
[posTask,velTask,accelTask] = trapveltraj(waypoints,numel(trajTimes), ...
    'AccelTime',repmat(waypointAccelTimes,[3 1]), ...
    'EndTime',repmat(diff(waypointTimes),[3 1]));

% Trajectory evaluation
for idx = 1:numel(trajTimes)
    % Solve IK
    tgtPose = trvec2tform(posTask(:,idx)');
    [config,info] = ik(eeName,tgtPose,ikWeights,ikInitGuess);
    ikInitGuess = config;
    qTask(:,idx) = config;
end

taskTime = toc;
disp(['Task space trajectory time : ' num2str(taskTime) ' s']);

%% Create and evaluate a joint space trajectory
ikInitGuess = jointAnglesHome';
ikInitGuess(ikInitGuess > pi) = ikInitGuess(ikInitGuess > pi) - 2*pi;
ikInitGuess(ikInitGuess < -pi) = ikInitGuess(ikInitGuess < -pi) + 2*pi;

disp('Running joint space trajectory generation and evaluation...')
tic;

% Solve IK for all waypoints
numWaypoints = size(waypoints,2);
numJoints = numel(gen3.homeConfiguration);
jointWaypoints = zeros(numJoints,numWaypoints);
for idx = 1:numWaypoints
    tgtPose = trvec2tform(waypoints(:,idx)');
    [config,info] = ik(eeName,tgtPose,ikWeights,ikInitGuess);
    cfgDiff = config - ikInitGuess;
    jointWaypoints(:,idx) = config';
end

% Trajectory Generation
[qJoint,qdJoint,qddJoint] = trapveltraj(jointWaypoints,numel(trajTimes), ...
    'AccelTime',repmat(waypointAccelTimes,[numJoints 1]), ...
    'EndTime',repmat(diff(waypointTimes),[numJoints 1]));

% Trajectory evaluation (only needed to find end effector position)
for idx = 1:numel(trajTimes)
    eeTform = getTransform(gen3,qJoint(:,idx)',eeName);
    posJoint(:,idx) = tform2trvec(eeTform)';
end
jointTime = toc;
disp(['Joint space trajectory time : ' num2str(jointTime) ' s']);

%% Create comparison plots
% Compare trajectories in Cartesian space
close all
figure, hold on
plot3(posTask(1,:),posTask(2,:),posTask(3,:),'b-');
plot3(posJoint(1,:),posJoint(2,:),posJoint(3,:),'r--');
plot3(waypoints(1,:),waypoints(2,:),waypoints(3,:),'ko','LineWidth',2);
title('Trajectory Comparison');
xlabel('X [m]');
ylabel('Y [m]');
zlabel('Z [m]');
legend('Task Space Trajectory','Joint Space Trajectory','Waypoints');
grid on
view([45 45]);

% Compare joint angles
% Plot each joint trajectory
for idx = 1:numJoints
    figure, hold on;
    plot(trajTimes,qTask(idx,:),'b-');
    plot(trajTimes,qJoint(idx,:),'r-');
    for wIdx = 1:numWaypoints
       xline(waypointTimes(wIdx),'k--');
    end
    title(['Joint ' num2str(idx) ' Trajectory']);
    xlabel('Time [s]');
    ylabel('Joint Angle [rad]');
    legend('Task Space Trajectory','Joint Space Trajectory');
end

author: mathworks-robotics

repo: trajectory-planning-robot-manipulators


text:
function im = VisualizeClassifier(jabfile)

%% params.
params = getParams;
npatches = params.npatches;
psize = params.psize;
nbins = params.nbins;
patchsz = params.patchsz;
scale = params.scale;

wd = params.wd;

%% load an example fly.

Q = load(jabfile,'-mat');

bdir = Q.x.expDirNames{1};
moviename = fullfile(bdir,'movie.ufmf');
trackfilename = fullfile(bdir,'trx.mat');
tracks = load(trackfilename);
tracks = tracks.trx;

[readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);
trackndx = 1; fly = 1;
fnum = tracks(fly).firstframe;
im1 = readfcn(fnum);

locy = round(tracks(fly).y(trackndx));
locx = round(tracks(fly).x(trackndx));
im1 = extractPatch(im1,...
  locy,locx,tracks(fly).theta(trackndx),patchsz);

im1 =[im1 im1];

%% Find Hog/HOF components.

ncl = numel(Q.x.classifierStuff);

for clnum = 1:ncl

  H{clnum} = zeros(npatches,npatches,nbins);
  F{clnum} = H{clnum};

  cl = Q.x.classifierStuff(clnum).params;
  fnames = Q.x.classifierStuff(clnum).featureNames;
  for ndx = 1:numel(cl)
    curf = fnames{cl(ndx).dim}{1};
    val = cl(ndx).alpha;
    prts = strsplit(curf,'_');
    loc = str2double(prts(end-2:end));
    if strcmp(curf(1:2),'hf')
      H{clnum}(loc(1),loc(2),loc(3))=H{clnum}(loc(1),loc(2),loc(3))+val;
    else
      F{clnum}(loc(1),loc(2),loc(3))=F{clnum}(loc(1),loc(2),loc(3))+val;
    end
  end
end
%% Draw them

% For HOG
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);
dt = mean(diff(bincenters));
binedges = [bincenters(1)-dt/2,(bincenters(1:end-1)+bincenters(2:end))/2,bincenters(end)+dt/2];

% For HOF
bincenters2 = bincenters*2;
dt = mean(diff(bincenters2));
binedges2 = [bincenters2(1)-dt/2,(bincenters2(1:end-1)+bincenters2(2:end))/2,bincenters2(end)+dt/2];

maxv = 3;
im = {};
for clnum = 1:ncl
  hfig = figure;
  clf;
  hax = axes('Position',[0,0,1,1]);
  set(hfig,'Units','pixels','Position',get(0,'ScreenSize'));

  im1curr = im1;

  him = imshow(imresize(im1curr,scale));
  axis image;
  truesize;
  colormap gray;
  hold on;
  axis off;

  colors = hsv(nbins);

  [nr,nc,~] = size(im1);
  nc = nc/2;

  h = [];
  for xi = 1:ceil(nc/psize),
    cx = (psize/2 + (xi-1)*psize)*scale + 1 ;
    if cx+psize/2 > nc*scale,
      break;
    end
    for yi = 1:ceil(nr/psize),
      cy = (psize/2 + (yi-1)*psize)*scale + 1;
      if cy+psize/2 > nr*scale,
        break;
      end

      for bini = 1:nbins,
        tmp = linspace(binedges2(bini),binedges2(bini+1),20);
        xcurr = cx + [0,psize/2*cos(tmp),0]*scale;
        ycurr = cy + [0,psize/2*sin(tmp),0]*scale;
        h(yi,xi,bini) = patch(xcurr,ycurr,colors(bini,:),'LineStyle','none','FaceAlpha',min(1,F{clnum}(yi,xi,bini)/maxv));
      end

    end
  end

  colors = hsv(nbins);
  colors = colors([ (end/2+1):end 1:end/2],:);
  hogpatch = [wd wd -wd -wd wd;-psize psize psize -psize -psize]/2;
  h = [];
  for xi = 1:ceil(nc/psize),
    cx = (psize/2 + 1 + (xi-1)*psize)*scale;
    if cx+psize/2 > nc*scale,
      break;
    end
    for yi = 1:ceil(nr/psize),
      cy = (psize/2 + 1 + (yi-1)*psize)*scale;
      if cy+psize/2 > nr*scale,
        break;
      end

      for bini = 1:nbins,
        tmp = bincenters(bini);
        curpatch = [cos(tmp) -sin(tmp); sin(tmp) cos(tmp)]*hogpatch;
        xcurr = nc*scale + cx + curpatch(1,:)*scale;
        ycurr = cy + curpatch(2,:)*scale;
        h(yi,xi,bini) = patch(xcurr,ycurr,colors(bini,:),'LineStyle','none','FaceAlpha',min(1,H{clnum}(yi,xi,bini)/maxv));
      end

    end
  end
  im{clnum} = getframe(hax);
end

author: kristinbranson

repo: JAABA


text:
function optsetup = spm_cfg_eeg_inv_optimize
%function optsetup = spm_cfg_eeg_inv_optimize
% configuration file to set up optimization routines for M/EEG source
% inversion
%_______________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes
% $Id: spm_cfg_eeg_inv_optimize.m 6499 2015-07-16 13:37:41Z gareth $

D = cfg_files;
D.tag = 'D';
D.name = 'M/EEG datasets';
D.filter = 'mat';
D.num = [1 Inf];
D.help = {'Select the M/EEG mat files.'};

val = cfg_entry;
val.tag = 'val';
val.name = 'Inversion index';
val.strtype = 'n';
val.help = {'Index of the cell in D.inv where the forward model can be found and the results will be stored.'};
val.val = {1};

priorname = cfg_entry;
priorname.tag = 'priorname';
priorname.name = 'Prefix for priors';
priorname.strtype = 's';
priorname.num = [1 Inf];
priorname.val = {'priorset1'};
priorname.help = {'Prefix for prior directory'};

REMLopt = cfg_entry;
REMLopt.tag = 'REMLopt';
REMLopt.name = 'REML parameters';
REMLopt.strtype = 'r';
REMLopt.num = [1 2];
REMLopt.val = {[-4,16]};
REMLopt.help = {'Select REML parameters'};

ARDopt = cfg_entry;
ARDopt.tag = 'ARDopt';
ARDopt.name = 'Threshold for ARD hyperparameter';
ARDopt.strtype = 'r';
ARDopt.num = [1 1];
ARDopt.val = {[128]};
ARDopt.help = {'ARD threshold for pruning hyperparameters (relative to max)'};

GSopt = cfg_entry;
GSopt.tag = 'GSopt';
GSopt.name = 'Set number of greed search iterations';
GSopt.strtype = 'i';
GSopt.num = [1 1];
GSopt.val = {[16]};
GSopt.help = {'Number of greedy search iterations'};

opttype = cfg_repeat;
opttype.tag = 'opttype';
opttype.name = 'Optimize';
opttype.help = {'Specify the optimization scheme'};
opttype.num  = [1 Inf];
opttype.values  = {REMLopt,ARDopt,GSopt};
opttype.val  = {REMLopt};

optsetup = cfg_exbranch;
optsetup.tag = 'optsetup';
optsetup.name = 'Inversion optimization';
optsetup.val = {D,val,priorname,opttype};
optsetup.help = {'Set optimization scheme for source reconstruction'};
optsetup.prog = @opt_priors;
optsetup.vout = @vout_opt_priors;

function  out = opt_priors(job)

D = spm_eeg_load(job.D{1});

%% get data specific terms sorted out
inverse=D.inv{job.val}.inverse;
AY=inverse.AY;

mesh=D.inv{job.val}.mesh.tess_mni;

[a1,b1,c1]=fileparts(D.fname);
priordir=[D.path filesep job.priorname '_' b1];

[priorfiles] = spm_select('FPListRec',priordir,'.*\.mat$');
%[postfiles] = spm_select('FPListRec',priordir,'^post.*\.mat$');

% USEPOST=1;
% if ~isempty(postfiles) && USEPOST,
%     fprintf('Using only posterior\n');
%     priorfiles=postfiles;
% end;
Nfiles=size(priorfiles,1);
fprintf('Found %d prior files\n',Nfiles);
if Nfiles==0,
    error('No prior file found in directory: %s', priordir);
end;

%% add in functional (from other modalities or experiment) hypotheses

Qe0=0;%% bounding ratio of noise to signal power
disp('NB NO min sensor noise level');  %% NO MIN SENSOR NOISE LEV

Qp_best=[];
Qe_best=[];
Fmax=-Inf;
F_aug=-Inf;
maxpriors=512;

for j=1:Nfiles, %% move through prior files
    %%% LOAD IN PRIOR FILE

    load(deblank(priorfiles(j,:)),'Qp','Qe','UL','F');
    fprintf('Optimizing priorfile %d of %d \n',j,Nfiles);
    %% ALSO CONSIDER AUGMENTED VERSION OF PRIORS IN FILE WITH BEST SO FAR
    Qp_aug=Qp; %% running with augmented Qp
    for k=1:length(Qp_best), %
        Qp_aug{length(Qp)+k}=Qp_best{k}; % augment with posterior from best so far
    end;
    Qe_aug=Qe; %% NB NOT AUGMENTING YET

    if length(Qp_aug)>maxpriors,
        Qp_aug=Qp_aug{1:maxpriors};
        warning('Limiting priors');
    end;
    %%%%%%%%%%%% NOW OPTIMIZE ORIGINAL AND AUGMENTED INDEPENDENTLY

    for k=1:length(job.opttype), %% move through optimization schemes
        %%% TERMS WHICH EVOLVE ARE Qe and Qp (M, Cq, Cp follow)
        priorcount(k)=length(Qp);

        [F,M,Cq,Cp,Qe,Qp] = spm_eeg_invert_EBoptimise(AY,UL,job.opttype(k),Qp,Qe,Qe0);

        if j>1, %% nothing to augment on 1st iteration

            [F_aug,M,Cq,Cp,Qe_aug,Qp_aug] = spm_eeg_invert_EBoptimise(AY,UL,job.opttype(k),Qp_aug,Qe_aug,Qe0);
        end;
    end;

    if F_aug>F,
        fprintf('Taking augmented set forward\n');
        Qp=Qp_aug;
        Qe=Qe_aug;
        F=F_aug;
    end;

    if F>Fmax, %% keep a record of these if they are best
        Qp_best=Qp;
        Qe_best=Qe;
        Fmax=F;
    end;

    Qp=Qp_best;
    Qe=Qe_best;
    save(deblank(priorfiles(j,:)),'Qp','Qe','UL','F');
    allF(j)=F;
end; % for j

%% Now get M, Cp, Cq etc based on Qp and Qe

[LCpL,Q,sumLCpL,QE,Cy,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp,{Qe});

inverse.F=Fmax;
inverse.M=M;
inverse.qC=Cq;
inverse.Cp=Cp; %% posterior source level
inverse.Qe=Qe; %% posterior sensor level
inverse.Qp=Qp;
inverse.Is=1:size(Cp,1); % %% temporary fix

%----------------------------------------------------------------------
% evaluate conditional expectation (of the sum over trials)
%----------------------------------------------------------------------
SSR   = 0;
SST   = 0;
J     = {};

UY=inverse.UY;

sourcevar=zeros(1,size(inverse.M,1));
for j = 1:numel(UY),

    % trial-type specific source reconstruction
    %------------------------------------------------------------------
    J{j} = inverse.M*UY{j};
    Jtime=J{j}*inverse.T'; %% J is Nvert* Nsamples
    sourcevar=sourcevar+var(Jtime'); %% ./inverse.qC';
    % sum of squares
    %------------------------------------------------------------------
    SSR  = SSR + sum(var((UY{j} - UL*J{j}))); %% changed variance calculation
    SST  = SST + sum(var( UY{j}));

end

inverse.J=J;

% accuracy; signal to noise (over sources)
%======================================================================
inverse.R2   = 100*(SST - SSR)/SST;
fprintf('Percent variance explained %.2f\n',full(inverse.R2));
D.inv{job.val}.inverse=inverse;

spm_eeg_invert_display(D);

rmind=find(allF<max(allF)-3);
fprintf('Removing %d poorest prior files\n',length(rmind));
for j=1:length(rmind),
    fprintf('Deleting %s\n',priorfiles(rmind(j),:));
    delete(deblank(priorfiles(rmind(j),:)));
end;

idnum=round(spm_data_id(sourcevar)*1000); %% get unique id for file
postfilename=[priordir filesep sprintf('post%d.mat',idnum)];
postgiftiname=[priordir filesep sprintf('post%d.gii',idnum)];
fprintf('Making posterior %s based on diagonal\n',postfilename);
Qp=[];
Qp{1}=sparse(diag(sourcevar));
Mmni=[];
Mmni.faces=D.inv{job.val}.mesh.tess_mni.face;
Mmni.vertices=D.inv{job.val}.mesh.tess_mni.vert;
Mmni.cdata=sourcevar';
Mmni=gifti(Mmni);

save(deblank(postfilename),'Qp','Qe','UL','Fmax','postgiftiname','-v7.3');
save(Mmni,postgiftiname);

figure;
spm_mip(sourcevar,mesh.vert',6);
title('Posterior');
colorbar;

D.save;
out.postname=postfilename;
out.postgiftiname=postgiftiname;
out.D = job.D;

function dep = vout_opt_priors(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'M/EEG dataset(s) after imaging source reconstruction';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function heat_oned_test01 ( )

%*****************************************************************************80
%
%% HEAT_ONED_TEST01 runs Jeff Borggaard's test case.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HEAT_ONED_TEST01:\n' );
  fprintf ( 1, '  Run Jeff Borggaard''s test case.\n' );
%
%  Set up the problem.
%
  mu = 0.1;
  h = 100.0;
  rho = 40.0;
  degree = 2;
  xa = 0.0;
  xb = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Viscosity parameter MU = %g\n', mu );
  fprintf ( 1, '  Heat transfer coefficient H = %g\n', h );
  fprintf ( 1, '  Mesh density RHO = %g\n', rho );
  fprintf ( 1, '  Requesting finite elements of degree %d\n', degree );
  fprintf ( 1, '  Spatial interval is [%g,%g].\n', xa, xb );

  [ A, B, M, n_nodes, x ] = heat_oned_setup ( mu, h, rho, degree, xa, xb );
%
%  Compute the solution over time.
%
  t_initial = 0.0;
  t_step    = 0.001;
  t_save    = 0.01;
  t_final   = 5.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial time =   %g\n', t_initial );
  fprintf ( 1, '  Final time =     %g\n', t_final );
  fprintf ( 1, '  Time step =      %g\n', t_step );
  fprintf ( 1, '  Data save step = %g\n', t_save );

  w_save = heat_oned ( A, B, M, n_nodes, x, t_initial, t_final, t_step, t_save );
%
%  Plot the results.
%
  t_plot = [ t_initial : t_save : t_final ];

  figure ( 1 )
  mesh ( t_plot, x, w_save )
  xlabel ( '<---Time--->' )
  ylabel ( '<---X--->' )
  zlabel ( '<--U(X,T)-->' )
  title ( 'U(X,T) as computed by HEAT\_ONED' )

  filename = 'heat_oned_test01.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved as "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%PRPROGRESS Report progress of some PRTools iterative routines
%
%  PRPROGRESS ON
%
% All progress of all routines will be written to the command window.
%
%  PRPROGRESS(FID)
%
% Progress reports will be written to the file with file descriptor FID.
%
%  PRPROGRESS(FID,FORMAT,...)
%
% Writes progress message to FID. If FID == [], the predefined destination
% (command window or file) is used.
%
%  PRPROGRESS OFF
%  PRPROGRESS(0)
%
% Progress reporting is turned off.
%
%  PRPROGRESS
%
% Toggles between PRPROGRESS ON and PRPROGRESS OFF
%
%  FID = PRPROGRESS
%
% Retrieves the status of PRPROGRESS
%
% Some routines (e.g. CLEVAL)  have a switch in the function call by which
% progress reporting for that routine only can be initiated.
%
% By default, PRPROGRESS is switched off. Interactive progress tracing can
% be best following by PRWAITBAR
%
% SEE ALSO (PRTools Guide)
% PRWAITBAR

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function fid = prprogress(par,varargin)

	persistent GLOBALPRPROGRESS
  mlock; % prevents GLOBALPRPROGRESS being cleared by CLEAR ALL

	if isempty(GLOBALPRPROGRESS)
		GLOBALPRPROGRESS = 0;
	end

	if nargin == 0

		if nargout == 0
			if GLOBALPRPROGRESS ~= 0,
				GLOBALPRPROGRESS = 0;
			else
				GLOBALPRPROGRESS = 1;
			end
		else
			fid = GLOBALPRPROGRESS;
		end

	elseif nargin == 1 & nargout == 0

		if isstr(par)
			switch par
			case {'on','ON'}
				GLOBALPRPROGRESS = 1;
			case {'off','OFF'}
				GLOBALPRPROGRESS = 0;
			otherwise
				error('Illegal input for PRPROGRESS')
			end
		else
			GLOBALPRPROGRESS = par;
		end

	elseif (GLOBALPRPROGRESS > 0) & (~isempty(par)) & (par>0)  %DXD, fid=0 means no printing...

		s = sprintf(varargin{:});
		n = fprintf(par,s);

		if nargout > 0
			fid = length(s);
		end

	elseif GLOBALPRPROGRESS > 0

		s = sprintf(varargin{:});
		n = fprintf(GLOBALPRPROGRESS,s);

		if nargout > 0
			fid = length(s);
		end

	else

		if nargout > 0
			fid = 0;
		end

	end

return

author: marianux

repo: ecg-kit

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


text:
function W = freqindex(W,varargin)

% WAVEFORM = FREQINDEX(WAVEFORM) calculates the frequency index of a waveform.
% The FI is stored as a new field called FFT_FI. This function must be run after
% WF_FFT.COMPUTE so that the necessary spectral fields exist in advance
%
% WAVEFORM = FREQINDEX(WAVEFORM,1) includes a plot of the waveform, and the
% frequency spectrum with the FI frequency bins shaded by their appropriate
% heights. If WAVEFORM is a vector, plots will only be included for the
% first five waveforms. This is a bit ad hoc, but it prevents accidental
% plotting of thousands of waveforms. To plot specific waveforms one can
% always use freqindex(w(n)) where n is the index of the nth waveform.
%
% See also wf_fft.compute wf_fft.plot

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

DOPLOT = 0;
if length(varargin) == 1
    DOPLOT = varargin{1};
end

% LOOP THROUGH WAVEFORMS (INEFFICIENT CODE!)
for i = 1:length(W)
    w = W(i);

    % GET FFT
    %W(i) = demean(W(i));
    %W(i) = wf_fft(W(i));
    F = get(W(i),'FFT_FREQ');
    A  = get(W(i),'FFT_AMP');
    if isempty(F) || isempty(A)
    	%W(i) = wf_fft.compute(W(i));
        error(['Trace ' num2str(i) ' does not have FFT_FREQ or FFT_AMP fields. Run WF_FFT.COMPUTE before WF_FFT.FREQINDEX' ]);
    end
    %A = A .* A;
    j = find(F >= 1 & F <= 2);
    Al = mean(A(j));
    j = find(F >= 10 & F <= 20);
    Ah = mean(A(j));
    fi = log10(Ah/Al);
    W(i) = addfield(W(i),'FFT_FI',fi);

    % PLOT IT
    if DOPLOT && i<=5

        figure('Color','w');
        set(gcf,'DefaultAxesFontSize',14);

        subplot(2,1,1);
        plot(w);
        title(datestr(get(w,'START'),31));

        subplot(2,1,2);
        hold on;
        fill( [1 1 2 2] , [0 Al Al 0] , [.7 .7 .7]);
        fill( [10 10 20 20] , [0 Ah Ah 0] , [.7 .7 .7]);
        plot( F , A , 'r-' );
        xlabel('Frequency (Hz)');
        ylabel('Power');
        set(gca,'XScale','Log');
        grid on; box on;
        xlim([0.1 50]);
        title(['Frequency index: ' num2str(fi,'%5.2f')]);

    end

end

author: geoscience-community-codes

repo: GISMO


text:
function [t0,t1,t2]=calculateInterpolationTimes(L,v,a)
%% About:
% calculates the times of the acceleration, stable motion, deceleration

%% Arreguments:
% L: length of the curve (mm).
% v: linear velocity of eef on the the curve (mm/sec).
% a: linear acceleration of eef on the the curve (mm/sec2).

%% Return values:
% t0; time for acceleration.
% t1: time for acceleration, constant velocity.
% t2: time for acceleration, constant velocity and for deceleration

% Copy right: Mohammad SAFEEA
% 16th-Oct-2017

taw=v/a;

if(L>a*taw*taw)
    t0=taw;
    deltat=(L-a*taw*taw)/v;
    t1=t0+deltat;
    t2=t1+taw;
else
    taw1=(L/a)^0.5;
    t0=taw1;
    t1=taw1;
    t2=2*taw1;
end

end

repo: KST-Kuka-Sunrise-Toolbox


text:
function D = gsp_distanz(X, Y, P)
%GSP_DISTANZ calculates the distances between all vectors in X and Y
%   Usage: D = gsp_distanz(X, Y);
%
%   Input parameters:
%       X   : matrix with col vectors
%       Y   : matrix with col vectors (default == X)
%       P   : distance matrix (default Identity)
%   Output parameters:
%       D   : distance matrix, not squared
%
%   This code computes the following
%
%   .. D = ( (X-Y)^T P (X-Y) )^(0.5)
%
%   .. math:: D = \sqrt{(X-Y)^T P (X-Y)}
%
%   for all vectors in X an Y!
%
%   This code is not optimized for memory, but for speed because it uses no
%   loops.
%

% Testing: test_gsp_distanz

% Handle Input parameters
if nargin<1, error('Not enought inputs parameters!'); end
if nargin<2,  Y=X; end

% Get the size
[rx, cx] = size(X);
[ry, cy] = size(Y);

% Verify the size
if rx~=ry, error('The sizes of x and y do not fit!'), end

if nargin < 3
    xx = sum(X.*X,1); % ||x||^2
    yy = sum(Y.*Y,1); % ||y||^2
    xy = X'*Y;        % <y,x>
    % \|x-y\|^2 = ||x||^2 +||y||^2 - 2 <y,x>
    D = abs(repmat(xx',[1 cy]) + repmat(yy,[cx 1]) - 2*xy);
else

    [rp,rp2] = size(P);
    if rx~=rp, error('The sizes of x and P do not fit!'), end
    if rp2~=rp, error('P must be square!'), end

    xx = sum(X .* (P* X),1 ); % x^T P x
    yy = sum(Y .* (P* Y),1 ); % y^T P y
    xy = X'*(P*Y);        % x^T P y
    yx = Y'*(P*X);        % x^T P y

    D = abs(repmat(xx',[1 cy]) + repmat(yy,[cx 1]) - xy-yx);
end

if sum(D(:)<0)
    warning('gsp_distanz: P is not semipositive or x is not real!')
end

% Take the square root
D = sqrt(D);

if nargin < 2   % if Y == X
    % The diagonal has to be zero!
    D(1:cx+1:end) = 0;
end

repo: gspbox


text:
function File = getFullPath(File, Style)
% getFullPath - Get absolute canonical path of a file or folder
% Absolute path names are safer than relative paths, when e.g. a GUI or TIMER
% callback changes the current directory. Only canonical paths without "." and
% ".." can be recognized uniquely.
% Long path names (>259 characters) require a magic initial key "\\?\" to be
% handled by Windows API functions, e.g. for Matlab's FOPEN, DIR and EXIST.
%
% FullName = getFullPath(Name, Style)
% INPUT:
%   Name:  String or cell string, absolute or relative name of a file or
%          folder. The path need not exist. Unicode strings, UNC paths and long
%          names are supported.
%   Style: Style of the output as string, optional, default: 'auto'.
%          'auto': Add '\\?\' or '\\?\UNC\' for long names on demand.
%          'lean': Magic string is not added.
%          'fat':  Magic string is added for short names also.
%          The Style is ignored when not running under Windows.
%
% OUTPUT:
%   FullName: Absolute canonical path name as string or cell string.
%          For empty strings the current directory is replied.
%          '\\?\' or '\\?\UNC' is added on demand.
%
% NOTE: The M- and the MEX-version create the same results, the faster MEX
%   function works under Windows only.
%   Some functions of the Windows-API still do not support long file names.
%   E.g. the Recycler and the Windows Explorer fail even with the magic '\\?\'
%   prefix. Some functions of Matlab accept 260 characters (value of MAX_PATH),
%   some at 259 already. Don't blame me.
%   The 'fat' style is useful e.g. when Matlab's DIR command is called for a
%   folder with les than 260 characters, but together with the file name this
%   limit is exceeded. Then "dir(getFullPath([folder, '\*.*], 'fat'))" helps.
%
% EXAMPLES:
%   cd(tempdir);                    % Assumed as 'C:\Temp' here
%   getFullPath('File.Ext')         % 'C:\Temp\File.Ext'
%   getFullPath('..\File.Ext')      % 'C:\File.Ext'
%   getFullPath('..\..\File.Ext')   % 'C:\File.Ext'
%   getFullPath('.\File.Ext')       % 'C:\Temp\File.Ext'
%   getFullPath('*.txt')            % 'C:\Temp\*.txt'
%   getFullPath('..')               % 'C:\'
%   getFullPath('..\..\..')         % 'C:\'
%   getFullPath('Folder\')          % 'C:\Temp\Folder\'
%   getFullPath('D:\A\..\B')        % 'D:\B'
%   getFullPath('\\Server\Folder\Sub\..\File.ext')
%                                   % '\\Server\Folder\File.ext'
%   getFullPath({'..', 'new'})      % {'C:\', 'C:\Temp\new'}
%   getFullPath('.', 'fat')         % '\\?\C:\Temp\File.Ext'
%
% COMPILE:
%   Automatic: InstallMex getFullPath.c uTest_getFullPath
%   Manual:    mex -O getFullPath.c
%   Download:  http://www.n-simon.de/mex
% Run the unit-test uTest_getFullPath after compiling.
%
% Tested: Matlab 6.5, 7.7, 7.8, 7.13, WinXP/32, Win7/64
%         Compiler: LCC2.4/3.8, BCC5.5, OWC1.8, MSVC2008/2010
% Assumed Compatibility: higher Matlab versions
% Author: Jan Simon, Heidelberg, (C) 2009-2013 matlab.THISYEAR(a)nMINUSsimon.de
%
% See also: CD, FULLFILE, FILEPARTS.

% $JRev: R-G V:032 Sum:7Xd/JS0+yfax Date:15-Jan-2013 01:06:12 $
% $License: BSD (use/copy/change/redistribute on own risk, mention the author) $
% $UnitTest: uTest_getFullPath $
% $File: Tools\GLFile\getFullPath.m $
% History:
% 001: 20-Apr-2010 22:28, Successor of Rel2AbsPath.
% 010: 27-Jul-2008 21:59, Consider leading separator in M-version also.
% 011: 24-Jan-2011 12:11, Cell strings, '~File' under linux.
%      Check of input types in the M-version.
% 015: 31-Mar-2011 10:48, BUGFIX: Accept [] as input as in the Mex version.
%      Thanks to Jiro Doke, who found this bug by running the test function for
%      the M-version.
% 020: 18-Oct-2011 00:57, BUGFIX: Linux version created bad results.
%      Thanks to Daniel.
% 024: 10-Dec-2011 14:00, Care for long names under Windows in M-version.
%      Improved the unittest function for Linux. Thanks to Paul Sexton.
% 025: 09-Aug-2012 14:00, In MEX: Paths starting with "\\" can be non-UNC.
%      The former version treated "\\?\C:\<longpath>\file" as UNC path and
%      replied "\\?\UNC\?\C:\<longpath>\file".
% 032: 12-Jan-2013 21:16, 'auto', 'lean' and 'fat' style.

% Initialize: ==================================================================
% Do the work: =================================================================

% #############################################
% ### USE THE MUCH FASTER MEX ON WINDOWS!!! ###
% #############################################

% Difference between M- and Mex-version:
% - Mex does not work under MacOS/Unix.
% - Mex calls Windows API function getFullPath.
% - Mex is much faster.

% Magix prefix for long Windows names:
if nargin < 2
    Style = 'auto';
end

% Handle cell strings: NOTE: It is faster to create a function
% @cell\getFullPath.m under Linux, but under Windows this would shadow the
% fast C-Mex.
if isa(File, 'cell')
    for iC = 1:numel(File)
        File{iC} = getFullPath(File{iC}, Style);
    end
    return;
end

% Check this once only:
isWIN    = strncmpi(computer, 'PC', 2);
MAX_PATH = 260;

% Warn once per session (disable this under Linux/MacOS):
persistent hasDataRead
if isempty(hasDataRead)
    % Test this once only - there is no relation to the existence of
    % DATAREAD!
    %if isWIN
    %   Show a warning, if the slower Matlab version is used - commented,
    %   because this is not a problem and it might be even useful when the
    %   MEX-folder is not inlcuded in the path yet.
    %   warning('JSimon:getFullPath:NoMex', ...
    %      ['getFullPath: Using slow Matlab-version instead of fast Mex.',
    %      ...
    %       char(10), 'Compile: InstallMex getFullPath.c']);
    %end

    % DATAREAD is deprecated in 2011b, but still available. In Matlab 6.5,
    % REGEXP does not know the 'split' command, therefore DATAREAD is
    % preferred:
    hasDataRead = ~isempty(which('dataread'));
end

if isempty(File)  % Accept empty matrix as input:
    if ischar(File) || isnumeric(File)
        File = cd;
        return;
    else
        error(['JSimon:', mfilename, ':BadTypeInput1'], ...
            ['*** ', mfilename, ': Input must be a string or cell string']);
    end
end

if ischar(File) == 0  % Non-empty inputs must be strings
    error(['JSimon:', mfilename, ':BadTypeInput1'], ...
        ['*** ', mfilename, ': Input must be a string or cell string']);
end

if isWIN  % Windows: --------------------------------------------------------
    FSep = '\';
    File = strrep(File, '/', FSep);

    % Remove the magic key on demand, it is appended finally again:
    if strncmp(File, '\\?\', 4)
        if strncmpi(File, '\\?\UNC\', 8)
            File = ['\', File(7:length(File))];  % Two leading backslashes!
        else
            File = File(5:length(File));
        end
    end

    isUNC   = strncmp(File, '\\', 2);
    FileLen = length(File);
    if isUNC == 0                        % File is not a UNC path
        % Leading file separator means relative to current drive or base
        % folder:
        ThePath = cd;
        if File(1) == FSep
            if strncmp(ThePath, '\\', 2)   % Current directory is a UNC path
                sepInd  = strfind(ThePath, '\');
                ThePath = ThePath(1:sepInd(4));
            else
                ThePath = ThePath(1:3);     % Drive letter only
            end
        end

        if FileLen < 2 || File(2) ~= ':'  % Does not start with drive letter
            if ThePath(length(ThePath)) ~= FSep
                if File(1) ~= FSep
                    File = [ThePath, FSep, File];
                else                        % File starts with separator:
                    File = [ThePath, File];
                end
            else                           % Current path ends with separator:
                if File(1) ~= FSep
                    File = [ThePath, File];
                else                        % File starts with separator:
                    ThePath(length(ThePath)) = [];
                    File = [ThePath, File];
                end
            end

        elseif FileLen == 2 && File(2) == ':'   % "C:" current directory on C!
            % "C:" is the current directory on the C-disk, even if the
            % current directory is on another disk! This was ignored in
            % Matlab 6.5, but modern versions considers this strange
            % behaviour.
            if strncmpi(ThePath, File, 2)
                File = ThePath;
            else
                try
                    File = cd(cd(File));
                catch    % No MException to support Matlab6.5...
                    if exist(File, 'dir')  % No idea what could cause an error then!
                        rethrow(lasterror);
                    else  % Reply "K:\" for not existing disk:
                        File = [File, FSep];
                    end
                end
            end
        end
    end

else         % Linux, MacOS: ---------------------------------------------------
    FSep = '/';
    File = strrep(File, '\', FSep);

    if strcmp(File, '~') || strncmp(File, '~/', 2)  % Home directory:
        HomeDir = getenv('HOME');
        if ~isempty(HomeDir)
            File(1) = [];
            File    = [HomeDir, File];
        end

    elseif strncmpi(File, FSep, 1) == 0
        % Append relative path to current folder:
        ThePath = cd;
        if ThePath(length(ThePath)) == FSep
            File = [ThePath, File];
        else
            File = [ThePath, FSep, File];
        end
    end
end

% Care for "\." and "\.." - no efficient algorithm, but the fast Mex is
% recommended at all!
if ~isempty(strfind(File, [FSep, '.']))
    if isWIN
        if strncmp(File, '\\', 2)  % UNC path
            index = strfind(File, '\');
            if length(index) < 4    % UNC path without separator after the folder:
                return;
            end
            Drive            = File(1:index(4));
            File(1:index(4)) = [];
        else
            Drive     = File(1:3);
            File(1:3) = [];
        end
    else  % Unix, MacOS:
        isUNC   = false;
        Drive   = FSep;
        File(1) = [];
    end

    hasTrailFSep = (File(length(File)) == FSep);
    if hasTrailFSep
        File(length(File)) = [];
    end

    if hasDataRead
        if isWIN  % Need "\\" as separator:
            C = dataread('string', File, '%s', 'delimiter', '\\');  %#ok<REMFF1>
        else
            C = dataread('string', File, '%s', 'delimiter', FSep);  %#ok<REMFF1>
        end
    else  % Use the slower REGEXP, when DATAREAD is not available anymore:
        C = regexp(File, FSep, 'split');
    end

    % Remove '\.\' directly without side effects:
    C(strcmp(C, '.')) = [];

    % Remove '\..' with the parent recursively:
    R = 1:length(C);
    for dd = reshape(find(strcmp(C, '..')), 1, [])
        index    = find(R == dd);
        R(index) = [];
        if index > 1
            R(index - 1) = [];
        end
    end

    if isempty(R)
        File = Drive;
        if isUNC && ~hasTrailFSep
            File(length(File)) = [];
        end

    elseif isWIN
        % If you have CStr2String, use the faster:
        %   File = CStr2String(C(R), FSep, hasTrailFSep);
        File = sprintf('%s\\', C{R});
        if hasTrailFSep
            File = [Drive, File];
        else
            File = [Drive, File(1:length(File) - 1)];
        end

    else  % Unix:
        File = [Drive, sprintf('%s/', C{R})];
        if ~hasTrailFSep
            File(length(File)) = [];
        end
    end
end

% "Very" long names under Windows:
if isWIN
    if ~ischar(Style)
        error(['JSimon:', mfilename, ':BadTypeInput2'], ...
            ['*** ', mfilename, ': Input must be a string or cell string']);
    end

    if (strncmpi(Style, 'a', 1) && length(File) >= MAX_PATH) || ...
            strncmpi(Style, 'f', 1)
        % Do not use [isUNC] here, because this concerns the input, which
        % can '.\File', while the current directory is an UNC path.
        if strncmp(File, '\\', 2)  % UNC path
            File = ['\\?\UNC', File(2:end)];
        else
            File = ['\\?\', File];
        end
    end
end

% return;
end

author: SysBioChalmers

repo: RAVEN


text:
function [kEst,thetaEst]=PoissonRateLikeGammaConjUpdate(xMeas,kEst,thetaEst)
%%POISSONRATELIKEGAMMACONJUPDATE When estimating the rate parameter of the
%               Poisson distribution, the conjugate prior distribution is a
%               central (lambda=0) gamma distribution. Given the prior
%               shape and rate parameters of the gamma distribution, this
%               function updates those parameters conditioned on the
%               measurements. The result is the parameters of the posterior
%               distribution, which is also a gamma distribution.
%
%INPUTS: xMeas An NX1 or 1XN set of N independent measurements of the
%              Poisson distribution.
%         kEst The shape parameter of the prior (central gamma)
%              distribution.
%     thetaEst The scale parameter of the prior distribution.
%
%OUTPUTS: kEst The shape parameter of the posterior (central gamma)
%              distribution.
%     thetaEst The scale parameter of the posterior distribution.
%
%A distribution that is conjugate prior to a particular likelihood function
%is such that after performing Bayes' rule to update using a measurement,
%one obtains the same type of distribution back again. This function just
%performs Bayes rule with the given measurements.
%
%The conjugate prior relation for the Poisson distribution is given in [1].
%
%REFERENCES:
%[1] D. Fink, "A compendium of conjugate priors," Montana State University,
%    Department of Biology, Environmental Statistics Group, Tech. Rep., May
%    1997. [Online].
%    Available: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.157.5540
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numMeas=length(xMeas);

kEst=kEst+sum(xMeas);
if(isfinite(thetaEst))
    thetaEst=thetaEst/(numMeas*thetaEst+1);
else%If an uninformativ prior is given (thetaEst=Inf), use the asymptotic
    %result.
    thetaEst=1/numMeas;
end
end

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


text:
function y = prtUtilLogisticSigmoid(x)
% y = prtUtilLogisticSigmoid(x)

y = 1./(1+exp(-x));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MATLAB Code for                                                  %
%                                                                   %
%  Multi-Objective Particle Swarm Optimization (MOPSO)              %
%  Version 1.0 - Feb. 2011                                          %
%                                                                   %
%  According to:                                                    %
%  Carlos A. Coello Coello et al.,                                  %
%  "Handling Multiple Objectives with Particle Swarm Optimization," %
%  IEEE Transactions on Evolutionary Computation, Vol. 8, No. 3,    %
%  pp. 256-279, June 2004.                                          %
%                                                                   %
%  Developed Using MATLAB R2009b (Version 7.9)                      %
%                                                                   %
%  Programmed By: S. Mostapha Kalami Heris                          %
%                                                                   %
%         e-Mail: sm.kalami@gmail.com                               %
%                 kalami@ee.kntu.ac.ir                              %
%                                                                   %
%       Homepage: http://www.kalami.ir                              %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function z=MyCost2(x)

    z1=sum(-10*exp(-0.2*sqrt(x(1:end-1).^2+x(2:end).^2)));

    z2=sum(abs(x).^0.8+5*(sin(x.^3)));

    z=[z1 z2];

end

author: sfvsfv

repo: Mathematical-modeling


text:
function xp = spm_dirichlet_exceedance(alpha,Nsamp)
% Compute exceedance probabilities for a Dirichlet distribution
% FORMAT xp = spm_dirichlet_exceedance(alpha,Nsamp)
%
% Input:
% alpha     - Dirichlet parameters
% Nsamp     - number of samples used to compute xp [default = 1e6]
%
% Output:
% xp        - exceedance probability
%__________________________________________________________________________
%
% This function computes exceedance probabilities, i.e. for any given model
% k1, the probability that it is more likely than any other model k2.
% More formally, for k1=1..Nk and for all k2~=k1, it returns p(x_k1>x_k2)
% given that p(x)=dirichlet(alpha).
%
% Refs:
% Stephan KE, Penny WD, Daunizeau J, Moran RJ, Friston KJ
% Bayesian Model Selection for Group Studies. NeuroImage (in press)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny & Klaas Enno Stephan
% $Id: spm_dirichlet_exceedance.m 3118 2009-05-12 17:37:32Z guillaume $

if nargin < 2
    Nsamp = 1e6;
end

Nk = length(alpha);

% Perform sampling in blocks
%--------------------------------------------------------------------------
blk = ceil(Nsamp*Nk*8 / 2^28);
blk = floor(Nsamp/blk * ones(1,blk));
blk(end) = Nsamp - sum(blk(1:end-1));

xp = zeros(1,Nk);
for i=1:length(blk)

    % Sample from univariate gamma densities then normalise
    % (see Dirichlet entry in Wikipedia or Ferguson (1973) Ann. Stat. 1,
    % 209-230)
    %----------------------------------------------------------------------
    r = zeros(blk(i),Nk);
    for k = 1:Nk
        r(:,k) = spm_gamrnd(alpha(k),1,blk(i),1);
    end
    sr = sum(r,2);
    for k = 1:Nk
        r(:,k) = r(:,k)./sr;
    end

    % Exceedance probabilities:
    % For any given model k1, compute the probability that it is more
    % likely than any other model k2~=k1
    %----------------------------------------------------------------------
    [y, j] = max(r,[],2);
    xp = xp + histc(j, 1:Nk)';

end
xp = xp / Nsamp;


text:
function varargout=lpcis2rf(varargin)
%
% For calling details please see v_lpcis2rf.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcis2rf.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcis2rf.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcis2rf(varargin{:});
else
    v_lpcis2rf(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function b = fastnnls(x,y,tol,b)

% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
%
% See also:
% 'unimodal' 'monreg' 'fastnnls'
%
%  FASTNNLS Fast non-negative least squares
%  The inputs are the matrix of predictor variables (x),
%  vector of predicted variable (y), and optional inputs
%  tolerance on the size of a regression coefficient that is
%  considered zero (tol), and initial guess for the regression
%  vector (b0). The output is the non-negatively constrained
%  least squares solution (b).
%
%  If tol is set to 0, the default tolerance will be used.
%
%  FASTNNLS is fastest when a good estimate of the regression
%  vector is input. This eliminates much of the computation
%  involved in determining which coefficients will be nonzero
%  in the final regression vector. This makes it very useful
%  in alternating least squares routines. Note that the input
%  b0 must be a feasible (i.e. nonnegative) solution.
%
%  The FASTNNLS algorithm is based on the one developed by
%  Bro and de Jong, J. Chemometrics, Vol. 11, No. 5, 393-401, 1997
%
%I/O: b = fastnnls(x,y,tol,b0);
%
%See also: MCR, PARAFAC

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% 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., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

[m,n] = size(x);
if (nargin < 3 | tol == 0)
  tol = max(size(x))*norm(x,1)*eps;
end
if nargin < 4
  b = zeros(n,1);
end

p = logical(zeros(1,n));
p(find(b>0)) = ones(size(find(b>0)));
r = ~p;

sp = x(:,p)\y;
b(find(p)) = sp;
while min(sp) < 0
  b(find(b<0)) = zeros(size(find(b<0)));
  p = logical(zeros(1,n));
  p(find(b>0)) = ones(size(find(b>0)));
  r = ~p;
  sp = x(:,p)\y;
  b(find(p)) = sp;
end

w = x'*(y-x*b);
[wmax,ind] = max(w);
flag = 0;
while (wmax > tol & any(r))
  p(ind) = 1;
  r(ind) = 0;
  sp = x(:,p)\y;
  while min(sp) < -tol
    tsp = zeros(n,1);
    tsp(find(p)) = sp;
    fb = find(b);
    rat = b(fb)./(eps+(b(fb)-tsp(fb)));
    alpha = min(rat(rat>tol));
    b = b + alpha*(tsp-b);
    p = b > tol;
    r = ~p;
    sp = x(:,p)\y;
  end
  b(find(p)) = sp;
  w = x'*(y-x*b);
  [wmax,ind] = max(w);
  if p(ind)
    wmax = 0;
  end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% Test whether it is necessary to rotate VB PCA
function nc2010_vbpca_mlexperiment(flatW, rotate, ncomps)

if nargin < 3
  ncomps = 100;
end

randn('state', 1);
rand('state', 1);

% Movie lens data with 100.000 ratings
U = load('/share/bayes/data/jaakko/movielens/u.data');
%Y = sparse(U(:,1), U(:,2), U(:,3));

% Divide into training (90%) and test sets (10%)
m = max(U(:,1));
n = max(U(:,2));
testset = (rand(rows(U),1) < 0.1);
Ytrain = full(sparse(U(~testset,1), U(~testset,2), U(~testset,3),m,n));
Ytest = full(sparse(U(testset,1), U(testset,2), U(testset,3),m,n));
Ytrain(Ytrain==0) = nan;
Ytest(Ytest==0) = nan;

%size(Ytrain)

if flatW
  disp('Using non-hierarchical flat prior for W');
  stringW = 'flatW';
else
  disp('Using hierarchical prior for W');
  stringW = 'hierW';
end

filename = sprintf(['/home/jluttine/matlab/neurocomputing2010/' ...
                    'nc2010_vbpca_mlexperiment_%s_ncomps=%d_rotate=%d'], ...
                   stringW, ncomps, rotate);

results = vbpcamv(Ytrain, ncomps, 'testset', Ytest, 'maxiters', 1e7, ...
                  'rotate', rotate, 'startupdatehyper', 1, ...
                  'startrotate', 1, 'autosavetime', 1, ...
                  'autosavefile', filename, 'fixw', flatW);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [u,Du,eqn,info] = PoissonQ1old(node,elem,pde,bdFlag, option)
%% POISSONQ1 Poisson equation: Q1 bilinear element.
%
%   u = PoissonQ1(node,elem,pde,bdFlag) produces the bilinear finite element
%   approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
%   The quad mesh is given by node and elem and the boundary edge is given by
%   bdFlag. See meshdoc, bddoc for details. The data is given by the
%   structure pde which contains function handles f, g_D, g_N, g_R, or d.
%   For general elliptic equations with convection and reaction
%   coefficients, see ellipticpde.
%
%   The function Poisson assembes the matrix equation AD*u = b and solves
%   it by the direct solver.
%
%   The diffusion coefficient d is a scalar function or a column array with
%   the same length as the elem array.
%
%   u = PoissonQ1(node,elem,pde,bdFlag,option) specifies the options.
%    - option.dquadorder: quadrature order for diffusion coefficients
%    - option.fquadorder: quadrature order for computing right hand side f
%
%   When only one type of boundary condition is imposed, the input argument
%   bdFlag can be skipped. The boundary condition is implicitly given in
%   the pde structure by specifying g_D or g_N only. See examples below.
%
%   [u,Du] = PoissonQ1(node,elem,pde,bdFlag) returns also the
%   non-modified stiffness matrix A, which is semi-definite. The kernel of
%   A consists of constant vectors. The matrix A can be used to evulate the
%   bilinear form A(u,v) = u'*A*v, especially the enery norm of a finite
%   element function u by sqrt(u'*A*u).
%
%   [u,Du,eqn] = PoissonQ1(node,elem,pde,bdFlag) returns also the equation
%   structure eqn, which includes:
%     - eqn.AD:  modified stiffness matrix AD;
%     - eqn.b:   the right hand side.
%   The solution u = AD\b. The output eqn can be used to test other solvers.
%
%
%   Example
%     clear all
%     node=  [0, 0; 1, 0; 1, 1; 0,1];
%     elem = [1,2,3,4];
%     for k = 1:4
%         [node,elem] = uniformrefinequad(node,elem);
%     end
%     % Homogenous Dirichlet boundary condition
%     pde.f = inline('ones(size(p,1),1)','p');
%     pde.g_D = 0;
%     u = PoissonQ1(node,elem,pde);
%     figure(1);
%     showsolution(node,elem,u);
%     % Non-homogenous Dirichlet boundary condition
%     pde.f = inline('-4*ones(size(p,1),1)','p');
%     pde.g_D = inline('p(:,1).^2 + p(:,2).^2','p');
%     u = PoissonQ1(node,elem,pde);
%     figure(2);
%     showsolution(node,elem,u);
%     % Homogenous Neumann boundary condition
%     clear pde
%     pde.f = inline('pi^2*cos(pi*p(:,1)).*cos(pi*p(:,2))','p');
%     u = PoissonQ1(node,elem,pde);
%     figure(3);
%     showsolution(node,elem,u);
%
% See also: Poisson3Q1
%
% Author: Huayi Wei < huayiwei1984@gmail.com>.
%
% Modified by Long Chen. Change the assembling process. Add mg solver and
% update boundary condition.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end
[N,Dim] = size(node);
[NT,NV] = size(elem);
Ndof = N;

tic;
%% Assemble stiffness matrix
% generate sparse pattern
ii = zeros(10*NT,1); jj = zeros(10*NT,1);
index = 0;
for i = 1:4
    for j = i:4
        ii(index+1:index+NT) = double(elem(:,i));
        jj(index+1:index+NT) = double(elem(:,j));
        index = index + NT;
    end
end
% quadrature points
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder')
    option.dquadorder = 2;        % default order is exact for quadratic function
end
[pts, w] = quadptsquad(option.dquadorder);
nQuad = size(pts,1);
% compute non-zeros
sA = zeros(10*NT,nQuad);
for p = 1:nQuad
    % Dphi at quadrature points
   [phi, Dphip, J] = quadbasis(node,elem,pts(p,:));
    index = 0;
    for i = 1:4
        for j = i:4
            Aij = 0;
            if isempty(pde.d) || isnumeric(pde.d)
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
            else
                pxy = zeros(NT, Dim);
                for ip = 1:Dim
                    xi = node(:,ip);
                    pxy(:,ip) = xi(elem)*phi;
                end
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
            end
            if ~isempty(pde.d) && isnumeric(pde.d) % d is piecewise constant
                Aij = pde.d.*Aij;
            end
            Aij = Aij.*J;
            sA(index+1:index+NT,p) = Aij;
            index = index + NT;
        end
    end
end
sA = sum(sA,2);
% assemble the matrix
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
A = A + AU + AU';
clear Aij ii jj Dphip

%% Assemble the right hand side
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end

if ~isempty(pde.f)
    [pts, weight] = quadptsquad(option.fquadorder);
    nQuad = size(pts,1);
    bt = zeros(NT,NV);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
        [phi, tempvar, J] = quadbasis(node,elem, pts(p,:));
        pxy = zeros(NT, Dim);
        for i = 1:Dim
            xi = node(:,i);
            pxy(:,i) = xi(elem)*phi; % ? questionable
        end
		fp = pde.f(pxy);
        bt = bt + (weight(p)*fp.*J)*phi';
    end
    b = accumarray(elem(:),bt(:),[N 1]);
end
clear pxy bt

%% Set up boundary conditions
if ~exist('bdFlag','var'), bdFlag = []; end
[AD,b,u,freeNode,isPureNeumann] = getbd(b);

%% Record assembling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeNode), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeNode) = AD(freeNode,freeNode)\b(freeNode);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
%         option.x0 = u;
        option.solver = 'CG';
        [u,info] = mg(AD,b,elem,option);
    case 'amg'
        option.solver = 'CG';
        [u(freeNode),info] = amg(AD(freeNode,freeNode),b(freeNode),option);
end
% post-process for pure Neumann problem
if isPureNeumann
    patchArea = accumarray(elem(:),[J;J;J;J]/4, [N 1]);
    uc = sum(u.*patchArea)/sum(J);
    u = u - uc;   % int u = 0
end

%% Output information
eqn = struct('A',AD,'b',b,'freeNode',freeNode);
info.assembleTime = assembleTime;

%% Compute Du
Du = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,b,u,freeNode,isPureNeumann] = getbd(b)
    %% Set up of boundary conditions.
    %
    % 1) Modify the matrix for Dirichlet boundary nodes, which are not degree
    % of freedom. Values at these nodes are evaluatation of pde.g_D. The
    % original stiffness matrix A is turn into the matrix AD by enforcing
    % AD(fixedNode,fixedNode)=I, AD(fixedNode,freeNode)=0, AD(freeNode,fixedNode)=0.
    %
    % 2) Modify the right hand side b. The Neumann boundary integral is added
    % to b. For Dirichlet boundary ndoes, b(fixedNode) is the evaluation of
    % pde.g_D.
    %
    % Special attentation should be given for the pure Neumann boundary
    % condition. To enforce the compatible condition, the vector b should have
    % mean value zero. To avoid a singular matrix, the 1st node is chosen as
    % fixedNode.
    %
    % The order of assigning Neumann and Dirichlet boundary condition is
    % important to get the right setting at the intersection nodes of Dirichlet
    % and Neumann boundary edges.
    %
    % Reference: Long Chen. Finite Element Methods and its Programming. Lecture
    % Notes.

    u = zeros(Ndof,1);
    %% Initial check
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    isRobin = (bdFlag(:) == 3);
    if any(isRobin)
        allEdge = [elem(:,[1,2]); elem(:,[2,3]); elem(:,[3,4]); elem(:,[4 1])];    % quad
        Robin = allEdge(isRobin,:);
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        ve = node(Robin(:,1),:) - node(Robin(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        mid = (node(Robin(:,1),:) + node(Robin(:,2),:))/2;
        % int g_R phi_iphi_j ds
        ii = [Robin(:,1),Robin(:,1),Robin(:,2),Robin(:,2)];
        jj = [Robin(:,1),Robin(:,2),Robin(:,1),Robin(:,2)];
        temp = pde.g_R(mid).*edgeLength;
        ss = [1/3*temp, 1/6*temp, 1/6*temp, 1/3*temp];
%         ss = [1/4*temp, 1/4*temp, 1/4*temp, 1/4*temp];
        A = A + sparse(ii,jj,ss,Ndof,Ndof);
    end

    % Find Dirichlet boundary nodes: fixedNode
    fixedNode = []; freeNode = [];
    if ~isempty(bdFlag) % find boundary edges and boundary nodes
        [fixedNode,bdEdge,isBdNode] = findboundary(elem,bdFlag);
        freeNode = find(~isBdNode);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given
        [fixedNode,bdEdge,isBdNode] = findboundary(elem);
        freeNode = find(~isBdNode);
    end
    isPureNeumann = false;
    if isempty(fixedNode) && isempty(Robin) % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedNode = 1;
        freeNode = 2:Ndof;    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedNode,fixedNode)=I, AD(fixedNode,freeNode)=0, AD(freeNode,fixedNode)=0.
    if ~isempty(fixedNode)
        bdidx = zeros(Ndof,1);
        bdidx(fixedNode) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    else
        AD = A;
    end

    %% Part 2: Find boundary edges and modify the right hand side b
    % Find boundary edges: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        Neumann = bdEdge;
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        [tempvar,Neumann] = findboundary(elem); %#ok<ASGLU>
    end

    % Neumann boundary condition
    if  isnumeric(pde.g_N) && all(pde.g_N == 0)
        pde.g_N = [];
    end
    if ~isempty(Neumann) && ~isempty(pde.g_N)
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 2;   % default order exact for linear gN
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        phigN = lambdagN;                 % linear bases
        nQuadgN = size(lambdagN,1);
        ge = zeros(size(Neumann,1),2);
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            for igN = 1:2
                ge(:,igN) = ge(:,igN) + weightgN(pp)*phigN(pp,igN)*gNp;
            end
        end
        ge = ge.*repmat(el,1,2);
        b = b + accumarray(Neumann(:), ge(:),[Ndof,1]);
    end
    % The case with non-empty Neumann edges but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann edges and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if isnumeric(pde.g_D) && all(pde.g_D == 0)   % zero g_D
        pde.g_D = [];
    end
    if ~isPureNeumann && ~isempty(fixedNode) && ~isempty(pde.g_D)
        if isnumeric(pde.g_D)  % pde.g_D could be a numerical array
            u(fixedNode) = pde.g_D(fixedNode);
        else % pde.g_D is a function handle
            u(fixedNode) = pde.g_D(node(fixedNode,:));
        end
        b = b - A*u;
    end
    if ~isPureNeumann % non-empty Dirichlet boundary condition
        b(fixedNode) = u(fixedNode);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b);   % compatilbe condition: sum(b) = 0
    end
    end % end of getbd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end % end of PoissonQ1

repo: ifem


text:
function [size1, size2] = size(f, varargin)
%SIZE   Size of a CLASSICFUN.
%   SIZE(F), where F is a CLASSICFUN, is the size of the ONEFUN of F.
%
% See also LENGTH.

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

% The size of a CLASSICFUN is the size of its onefun.
if ( nargout <= 1 )
    size1 = size(f.onefun, varargin{:});
else
    [size1, size2] = size(f.onefun, varargin{:});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [varargout]=ellipsefit(x,y)
%ELLIPSEFIT Stable Direct Least Squares Ellipse Fit to Data.
% [Xc,Yc,A,B,Phi,P]=ELLIPSEFIT(X,Y) finds the least squares ellipse that
% best fits the data in X and Y. X and Y must have at least 5 data points.
% Xc and Yc are the x- and y-axis center of the ellipse respectively.
% A and B are the major and minor axis of the ellipse respectively.
% Phi is the radian angle of the major axis with respect to the x-axis.
% P is a vector containing the general conic parameters of the ellipse.
% The conic representation of the ellipse is given by:
%
% P(1)*x^2 + P(2)*x*y + P(3)*y^2 + P(4)*x + P(5)*y + P(6) = 0
%
% S=ELLIPSEFIT(X,Y) returns the output data in a structure with field names
% equal to the variable names given above, e.g., S.Xc, S.Yc, S.A, S.B,
% S.Phi and S.P
%
% Reference: R. Halif and J. Flusser, "Numerically Stable Direct Least
% Squares FItting of Ellipses," Department of Software Engineering, Charles
% University, Czech Republic, 2000.

% Conversion from conic to conventional ellipse equation inspired by
% fit_ellipse.m on MATLAB Central

% D.C. Hanselman, University of Maine, Orono, ME 04469
% Mastering MATLAB 7
% 2005-02-28
% Rotation angle fixed 2005-08-09

%--------------------------------------------------------------------------
x=x(:); % convert data to column vectors
y=y(:);
if numel(x)~=numel(y) || numel(x)<5
   error('X and Y Must be the Same Length and Contain at Least 5 Values.')
end

D1=[x.*x x.*y y.*y]; % quadratic terms
D2=[x y ones(size(x))]; % linear terms
S1=D1'*D1;
S2=D1'*D2;

[Q2,R2]=qr(D2,0);
if condest(R2)>1.0e10
   warning('ellipsefit',...
      'Data is Poorly Conditioned and May Not Represent an Ellipse.')
end
T=-R2\(R2'\S2'); % -inv(S3) * S2'

M=S1+S2*T;
CinvM=[M(3,:)/2; -M(2,:); M(1,:)/2];
[V,na]=eig(CinvM);
c=4*V(1,:).*V(3,:) - V(2,:).^2;
A1=V(:,c>0);
P=[A1; T*A1];

% correct signs if needed
if ~isempty(P),
  P=sign(P(1))*P;

  Phi=atan(P(2)/(P(3)-P(1)))/2;
  c=cos(Phi);
  s=sin(Phi);

  % rotate the ellipse parallel to x-axis
  Pr=zeros(6,1);
  Pr(1)=P(1)*c*c - P(2)*c*s + P(3)*s*s;
  Pr(2)=2*(P(1)-P(3))*c*s + (c^2-s^2)*P(2);
  Pr(3)=P(1)*s*s + P(2)*s*c + P(3)*c*c;
  Pr(4)=P(4)*c - P(5)*s;
  Pr(5)=P(4)*s + P(5)*c;
  Pr(6)=P(6);

  % extract other data
  XcYc=[c s;-s c]*[-Pr(4)/(2*Pr(1));-Pr(5)/(2*Pr(3))];
  Xc=XcYc(1);
  Yc=XcYc(2);
  F=-Pr(6) + Pr(4)^2/(4*Pr(1)) + Pr(5)^2/(4*Pr(3));
  AB=sqrt(F./Pr(1:2:3));
  A=AB(1);
  B=AB(2);
  Phi=-Phi;
  if A<B % x-axis not major axis, so rotate it pi/2
    Phi=Phi-sign(Phi)*pi/2;
    A=AB(2);
    B=AB(1);
  end
  S.Xc=Xc;
  S.Yc=Yc;
  S.A=A;
  S.B=B;
  S.Phi=Phi;
  S.P=P;
else
  S.Xc=nan;
  S.Yc=nan;
  S.A=nan;
  S.B=nan;
  S.Phi=nan;
  S.P=nan;
end
if nargout==1
  varargout{1}=S;
else
  outcell=struct2cell(S);
  varargout=outcell(1:nargout);
end