text:
%CORNERMINEIGENVAL  Calculates the minimal eigenvalue of gradient matrices for corner detection
%
%     dst = cv.cornerMinEigenVal(src)
%     dst = cv.cornerMinEigenVal(src, 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ Input single-channel 8-bit or floating-point image.
%
% ## Output
% * __dst__ Image to store the minimal eigenvalues. It has the same size as
%   `src` and the `single` type (single-channel).
%
% ## Options
% * __BlockSize__ Neighborhood size (see the details on
%   cv.cornerEigenValsAndVecs). default 5.
% * __KSize__ Aperture parameter for the cv.Sobel operator. default 3.
% * __BorderType__ Pixel extrapolation method. See cv.copyMakeBorder.
%   default 'Default'
%
% The function is similar to cv.cornerEigenValsAndVecs but it calculates and
% stores only the minimal eigenvalue of the covariance matrix of derivatives,
% that is, `min(lambda_1,lambda_2)` in terms of the formulae in the
% cv.cornerEigenValsAndVecs description.
%
% See also: cv.cornerEigenValsAndVecs, detectMinEigenFeatures, cornerPoints
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function x = is_Nd_distributed(x)
x = iscell(x) && length(x) == 2 && ischar(x{1}) && strcmp(x{1},'mvn') && iscell(x{2});

repo: REST


text:
function imgOut = RotateLLGUI(img)
%
%        imgOut = RotateLLGUI(img)
%
%        This tool rotates an environment map according to two picking
%        points.
%
%        Input:
%           -img: an environment map encoded as logitude-latitude
%
%        Output:
%           -imgOut: img rotated
%
%     Copyright (C) 2016 Francesco Banterle
%
%     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/>.
%

if(size(img, 2) > 2048)
   img = imresize(img, [2048, 1024], 'bilinear');
end

[r, c, col] = size(img);

figure(1);
imshow(img.^0.45);
hold on;

r_h = round(r / 2);
line([1, c], [r_h, r_h], 'Color', 'r', 'LineWidth', 4);

[x0, y0] = ginput(1);

plot(x0, y0, 'r+');

[x1, y1] = ginput(1);

plot(x1, y1, 'r+');
hold off;

%theta1 = ( (x1 - c/2) / c ) * pi;
%phi1   = ( (r/2 - y1) / r ) * pi * 0.5;

[theta1, phi1] = getThetaPhi(x1,y1,r,c);
vec1 = PolarVec3(theta1, phi1);

%thetar1 = ( (x0 - c/2) / c ) * pi;
%phir1   = ( (r/2 - y1) / r ) * pi * 0.5;
[thetar1, phir1] = getThetaPhi(x0,y1,r,c);
vecr1 = PolarVec3(thetar1, phir1);

M = getMatrixForVectorRotation(vec1, vecr1);

disp(M)
% disp(vec1);
% disp(vecr1);
% disp((M*vec1')');

D = LL2Direction(r, c);
D_rot = RotateMap(D, M');

[X1, Y1] = Direction2LL(D_rot, r, c);
[X, Y] = meshgrid(0:(c-1), 0:(r-1));
X1 = real(round(X1));
Y1 = real(round(Y1));
imgOut = zeros(size(img));
size(X)
size(X1)
for i=1:col
    imgOut(:,:,i) = interp2(X, Y, img(:,:,i), X1, Y1, 'spline');
end

hold on;
imshow(imgOut.^0.45);
plot([x0, x1], [y0, y1], 'r');
plot(x0, y0, 'go');
plot(x1, y1, 'go');
hold off;

end

function [theta, phi] = getThetaPhi(x,y,r,c)
    phi   =  pi * ((x / c) * 2 - 1) - pi * 0.5;
    theta =  pi * (y / r);
end

function M = getMatrixForVectorRotation(u, v)
    a = cross(u, v);
    len_a = norm(a);
    if(len_a > 0.0)
        a = a / len_a;

        alpha = acos(dot(u, v));
        c = cos(alpha);
        s = sin(alpha);

        ci = 1.0 - c;

        M = [a(1)^2*ci + c, a(1)*a(2)*ci - s * a(3), a(1)*a(3)*ci+s*a(2);...
             a(1)*a(2) * ci + s * a(3), a(2)^2 * ci + c, a(2)*a(3)*ci-s*a(1);...
             a(1)*a(3) * ci - s*a(2), a(1)*a(3)*ci+s*a(1), a(3)^2 * ci + c];
    else
        M = diag([1.0,1.0,1.0]);
    end
end

function D_rot = RotateMap(D, M)
    D_rot = zeros(size(D));
    D_rot(:,:,1) = D(:,:,1) * M(1,1) + D(:,:,2) * M(1,2) + D(:,:,3) * M(1,3);
    D_rot(:,:,2) = D(:,:,1) * M(2,1) + D(:,:,2) * M(2,2) + D(:,:,3) * M(2,3);
    D_rot(:,:,3) = D(:,:,1) * M(3,1) + D(:,:,2) * M(3,2) + D(:,:,3) * M(3,3);
end

author: banterle


text:
function z = maxdiff(x, y)
% Written by Mo Chen (sth4nth@gmail.com).
assert(all(size(x)==size(y)));
z = max(abs(x(:)-y(:)));

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


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

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

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

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

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

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

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

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

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

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

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

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

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

%% Post-process
if length(pDof)~=Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(p.*area)/sum(area);
    p = p - c;
end
w = R*u + utbd./vecMv;
psi = zeros(N,1);
% compute streamline function
if isfield(option,'stream') && option.stream
    [As,Ms] = assemblematrix(node,elem);
    [fixedNode,bdEdge,isBdNode] = findboundary(elem,bdFlag);
    freeNode = ~isBdNode;
    rhs = Ms*w;
    streamoption.freeDof = freeNode;
    psi(freeNode) = mg(As(freeNode,freeNode),rhs(freeNode),elemunSort,streamoption);
    % assume u\cdot n = 0.
end

%% Output information
info.assembleTime = assembleTime;

%% Output
soln = struct('u',u,'p',p,'w',w,'psi',psi);
eqn = struct('A',A0,'B',B0,'Me',Me,'Mv',Mv,'f',f0,'g',g0,...
             'edge',edge,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

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

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

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

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

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

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

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

repo: ifem


text:
function varargout = semilogy(varargin)
%SEMILOGY   Semi-log scale plot of a CHEBFUN.
%   SEMILOGY(...) is the same as PLOT(...), except a logarithmic (base 10) scale
%   is used for the Y-axis.
%
% See also PLOT, SEMILOGX, LOGLOG.

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

% Standard CHEBFUN/PLOT():
[h{1:3}] = plot(varargin{:});

% Find the CHEBFUNS:
isCheb = cellfun(@(v) isa(v, 'chebfun'), varargin);
% Choose a TOL based on vscale:
tol = max(cellfun(@(f) max(eps*vscale(f)), varargin(isCheb)));
% Loop over the different plot components:
for j = 1:numel(h)
    % Get the y data:
    yData = get(h{j}, 'yData');
    % Ensure it's a cell:
    if ( ~iscell(yData) )
        yData = {yData};
    end
    % Loop over each cell:
    for k = 1:numel(yData)
        if ( min(yData{k}(:)) > -tol )
            % If it's only a little negative, then take the absolute value:
            yData{k} = abs(yData{k});
        end
    end
    for k = 1:numel(h{j})
        % Put it back in h
        set(h{j}(k), 'yData', yData{k});
    end
end

% Set the YScale to be logarithmic:
set(gca, 'YScale', 'log');

% Output handle if requested:
if ( nargout > 0 )
    varargout = {h};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y=normminmax(x)
% bulld RBG image
xmin=min(x(:));
xmax=max(x(:));
y=(x-xmin)/(xmax -xmin);
return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% This file is the example of constrained formulation of SINDy-PI
% Coded By: K
% Last Update: 2019/07/09
%%
clc;close all;clear all;
%%
% Define the parameters for simulation
dt=0.01;T=3;

tspan=0:dt:T;

jx=0.6;Vmax=1.5;Km=0.3;

x0=1;

% Simulate
[t,x]=ode45(@(t,x)MMK_ODE(t,x,jx,Vmax,Km),tspan,x0);
dx=MMK_ODE(0,x',jx,Vmax,Km)';

% Now build library
Theta=[ones(size(x)) x x.^2 dx dx.*x dx.*x.^2];

[n,m]=size(Theta);
Diag_m=eye(m,m);
% Create the initial guess of the
C0=Diag_m(:)';

% Set the threshold
lambda=0.2;

% Determine how many iteration you need
N=5;

% Begin the optimization problem
C=C0;
tic
for iter=1:N
    % Run the optimization
    cvx_begin quiet
        variable xi(m,m)
            minimize(norm(Theta-Theta*xi));
        subject to
            diag(xi)==zeros(m,1);
            if iter>1
                xi(smallinds)==0;
            end
    cvx_end
    % Use thresholding
    fprintf('\n Iteration %d\n',iter)
    Xi=full(xi)
    smallinds = (abs(Xi)<lambda);
end
toc

author: dynamicslab

repo: SINDy-PI


text:
function k = kernCompute(kern, x, x2)

% KERNCOMPUTE Compute the kernel given the parameters and X.
% FORMAT
% DESC computes a kernel matrix for the given kernel type given an
% input data matrix.
% ARG kern : kernel structure to be computed.
% ARG X : input data matrix (rows are data points) to the kernel computation.
% RETURN K : computed elements of the kernel structure.
%
% FORMAT
% DESC computes a kernel matrix for the given kernel type given
% two input data matrices, one for the rows and one for the columns.
% ARG kern : kernel structure to be computed.
% ARG X : first input matrix to the kernel computation (forms the rows of the kernel).
% ARG X2 : second input matrix to the kernel computation (forms the columns of the kernel).
% RETURN K : computed elements of the kernel structure.
%
% SEEALSO : kernCreate, kernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

fhandle = str2func([kern.type 'KernCompute']);
if nargin < 3
  k = fhandle(kern, x);
else
  k = fhandle(kern, x, x2);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function u = poisson( f, varargin )
%POISSON   Fast Poisson solver for the rectangle.
%   POISSON(F) solves laplacian(U) = F on the domain of F with zero
%   Dirichlet boundary conditions. That is, U satisfies
%
%     U_{x,x} + U_{y,y} = F, on [a,b]x[c,d], with U = 0 on boundary
%
%   The equation is solved using an adaptively determined discretization
%   size.
%
%   POISSON(F, G) solves using Dirichlet boundary conditions given by G. G
%   can be a scalar, a function handle, or any chebfun2 object satisfying
%   the Dirichlet data.
%
%   POISSON(F, G, N) is the same as POISSON(F, G), but uses an N x N tensor
%   product discretization to solve the equation.
%
%   POISSON(F, G, M, N) is the same as POISSON(F, G, N), but with an M x N
%   tensor product discretization, where N is the number of coeffcieints in
%   the x-direction and M is the number in the y-direction.
%
%   POISSON(F, G, N, METHOD) or POISSON(F, G, M, N, METHOD) is the same as
%   POISSON(F, G, N) or POISSON(F, G, M, N), respectively, except the
%   underlying matrix equation is solved with METHOD. Available methods
%   are:
%
%     'adi'            - alternating direction implicit method
%     'fadi'           - factored alternating direction implicit method
%     'bartelsStewart' - Bartels-Stewart algorithm
%
%   If METHOD is not supplied, then this command selects one (based on the
%   discretization size and the rank of the righthand side).
%
% EXAMPLE:
%   f = chebfun2( @(x,y) 1 + 0*x, [-1 2 0 1]);
%   u = chebfun2.poisson(f);
%   plot(u)
%
% See also DISKFUN/POISSON, SPHEREFUN/POISSON.

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

% DEVELOPER'S NOTE:
%
% METHOD: Spectral method (in coefficient space). We use a C^{(3/2)} basis
% to discretize the equation, resulting in a discretization of the form
% AX + XA = F, where A is a symmetric tridiagonal matrix.
%
% LINEAR ALGEBRA: Matrix equations. The matrix equation is solved by the
% alternating direction implicit (ADI) method.
%
% SOLVE COMPLEXITY:  O(M*N*log(MAX(M,N))*log(1/eps)) with M*N = total
% degrees of freedom.
%
% AUTHORS: Dan Fortunato (dan.fortunato@gmail.com)
%          Alex Townsend (townsend@cornell.edu)
%
% The fast Poisson solver is based on:
%
% D. Fortunato and A. Townsend, Fast Poisson solvers for spectral methods,
% Submitted, 2017.

% Solve for u on the same domain as f, adjust diffmat to include scaling:
dom = f.domain;
scl_x = (2/(dom(2)-dom(1)))^2;
scl_y = (2/(dom(4)-dom(3)))^2;

% Not enough input arguments so we call chebop2 to adaptively select a
% discretization size:
if ( nargin == 1 )

    % Call is POISSON(F) so set G = 0 and use chebop2 to determine the
    % discretization size:
    g = 0;
    N = chebop2.setupLaplace( f.domain );
    N.bc = g;
    u = N \ f;
    return

elseif ( nargin == 2 )

    % Call is POISSON(F, G) so use chebop2 to determine the discretization
    % size:
    g = varargin{1};
    if ( ~isa(g, 'chebfun2') )
        g = chebfun2(g, f.domain);
    end
    N = chebop2.setupLaplace( f.domain );
    % Note: chebfun2/subsref (e.g. g(-1,:)) does not work here, so we use
    %       feval instead. Is this a bug?
    N.lbc = feval(g, f.domain(1), ':');
    N.rbc = feval(g, f.domain(2), ':');
    N.dbc = feval(g, ':', f.domain(3));
    N.ubc = feval(g, ':', f.domain(4));
    u = N \ f;
    return

end

% We are given a discretization size. It's solve time!
g = varargin{1};
m = varargin{2};
method = '';
if ( nargin == 3 )
    % Call must be POISSON(F, G, N) so employ an NxN discretization:
    n = m; % square discretization
elseif ( nargin == 4 )

    % The user typed one of the following:
    % POISSON(F, G, N, METHOD) or POISSON(F, G, M, N).

    if ( ischar( varargin{3} ) )
        % Call is POISSON(F, G, N, METHOD):
        method = varargin{3};
        n = m; % square discretization
    else
        % Call is POISSON(F, G, M, N):
        n = varargin{3};
    end

elseif ( nargin == 5 )

    % We are given a discretization size and a method. It's solve time!
    % Call must be POISSON(F, G, M, N, METHOD).
    n = varargin{3};
    method = varargin{4};

else
    error('CHEBFUN2:POISSON:NARGIN', ...
        'Too many input arguments to chebfun2.poisson().');
end

% Set the error tolerance for solve:
tol = chebfun2eps();

% Compute the Chebyshev coefficients of rhs:
[Cf, Df, Rf] = coeffs2(f, m, n);

% Solver only deals with zero homogeneous Dirichlet conditions. Therefore,
% if nonzero Dirichlet conditions are given, we solve lap(u) = f with u|bc = g
% as u = v + w, where v|bc = g, and lap(w) = f - lap(v), w|bc = 0:
if ( isa(g, 'double') || isa(g, 'chebfun2') || isa(g, 'function_handle') )

    % Make double or function handle into chebfun2:
    if ( isa(g, 'double') || isa(g, 'function_handle') )
        g = chebfun2(g, f.domain);
    end

    if ( g.domain == f.domain )
        % Adjust the rhs, if nonzero:
        lapg = lap(g);
        if ( ~iszero( lapg ) )
            [Cg, Dg, Rg] = coeffs2(lapg, m, n);
            Cf = [Cf Cg];
            Z = zeros(size(Df,1),size(Dg,1));
            Df = [Df Z ; Z' -Dg];
            Rf = [Rf Rg];
        end
    else
        error('CHEBFUN2:POISSON:BC', ...
            'Dirichlet data should be on the same domain as F.');
    end

else
    error('CHEBFUN2:POISSON', ...
        'Dirichlet data needs to be given as a scalar or function.')
end

% Convert rhs to C^{(3/2)} coefficients:
Cf = cheb2ultra( Cf );
Rf = cheb2ultra( Rf );

% Construct M, the multiplication matrix for (1-x^2) in the C^(3/2) basis
jj = (0:n-1)';
dsub = -1./(2*(jj+3/2)).*(jj+1).*(jj+2)*1/2./(1/2+jj+2);
dsup = -1./(2*(jj+3/2)).*(jj+1).*(jj+2)*1/2./(1/2+jj);
d = -dsub - dsup;
Mn = spdiags([dsub d dsup], [-2 0 2], n, n);
% Construct D^{-1}, which undoes the scaling from the Laplacian identity
invDn = spdiags(-1./(jj.*(jj+3)+2), 0, n, n);
Tn = scl_y * invDn * Mn;

jj = (0:m-1)';
dsub = -1./(2*(jj+3/2)).*(jj+1).*(jj+2)*1/2./(1/2+jj+2);
dsup = -1./(2*(jj+3/2)).*(jj+1).*(jj+2)*1/2./(1/2+jj);
d = -dsub - dsup;
Mm = spdiags([dsub d dsup], [-2 0 2], m, m);
invDm = spdiags(-1./(jj.*(jj+3)+2), 0, m, m);

% Construct T = D^{-1} * M:
Tm = scl_x * invDm * Mm;
Cf = invDm * Cf;
Rf = invDn * Rf;

switch lower(method)

    case 'bartelsstewart'

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%  Bartels-Stewart method %%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Solve TmX + XTn' = F using Bartels-Stewart, which requires O(n^3)
        % operations:

        X = chebop2.bartelsStewart(Tm, eye(n), eye(m), Tn, Cf*Df*Rf.', 0, 0);

        % Convert back to Chebyshev
        X = ultra1mx2cheb( ultra1mx2cheb( X ).' ).';
        u = chebfun2( X, f.domain, 'coeffs' );

    case {'adi', 'fadi', ''}

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%  Alternating Direction Implicit method %%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Solve TmX + XTn' = F using ADI, which requires
        % O(n^2log(n)log(1/eps)) operations:

        % An ADI method will be used (either given by the user, or selected
        % by us.)

        % Compute ADI shifts
        a =  -4/pi^2 * scl_y;
        b = -39*n^-4 * scl_y;
        c =  39*m^-4 * scl_x;
        d =   4/pi^2 * scl_x;
        [p, q] = chebop2.adiShifts(a, b, c, d, tol);

        if ( isempty(method) )
            % Let's go and pick a good method to use:
            % Test if we should use ADI or FADI:
            rho = size(Cf,2); % Rank of rhs
            adi_test = ( min(m,n) < rho*numel(p)/2 ); % Worth doing FADI?
            if ( adi_test )
                method = 'adi';
            else
                method = 'fadi';
            end
        end

        % Solve matrix equation:
        if ( strcmpi(method, 'adi') )
            % Run the ADI method:
            X = chebop2.adi(Tm, -Tn', Cf*Df*Rf.', p, q );

            % Convert back to Chebyshev
            X = ultra1mx2cheb( ultra1mx2cheb( X ).' ).';
            u = chebfun2( X, f.domain, 'coeffs' );

        else
            % Run the FADI method:
            [UX, DX, VX] = chebop2.fadi(Tm, -Tn, Cf*Df, Rf, p, q);

            % Convert back to Chebyshev:
            UX = ultra1mx2cheb(UX);
            VX = ultra1mx2cheb(VX);

            UX = chebfun(UX, dom(3:4), 'coeffs');
            VX = chebfun(VX, dom(1:2), 'coeffs');
            u = chebfun2();
            u.cols = UX;
            u.rows = VX;
            u.pivotValues = 1./diag(DX);
            u.domain = dom;
        end

    otherwise
        error('CHEBFUN2:POISSON:SOLVER', ...
            'Method supplied to chebfun2.poisson() is not recognized.');
end

% Add back in the boundary data:
u = u + g;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% CONVERSION CODES %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function X = cheb2ultra( X )
% CHEB2ULTRA   Convert vector of Chebyshev coefficients to C^(3/2).
%
%    CHEB2ULTRA(X) applies the conversion to each column of X if X is a
%    matrix.

% First convert the matrix of Chebyshev coefficients to a matrix of
% Legendre coefficients:
m = size( X, 1 );
if ( m <= 10000 ) % Determined experimentally
    S = cheb2leg_mat( m );
    X = S * X;
else
    X = cheb2leg( X );
end

% Now, convert the matrix of Legendre coefficients to a matrix of
% ultraspherical coefficients:
S = leg2ultra_mat( m );
X = S * X;

end

function X = ultra1mx2cheb( X )
% ULTRA1MX2CHEB    Convert vector of (1-x^2)C^(3/2) coefficients to
% Chebyshev.
%
%  ULTRA1MX2CHEB(X) applies the conversion each column of X if X is
%    matrix.

% First, convert the matrix of (1-x^2)C^(3/2)(x) coefficients
% to Legendre coefficients:
m = size( X, 1 );

S = ultra1mx2leg_mat( m );
X = S * X;

% Now, convert the matrix of Legendre coefficient to a matrix of Chebyshev
% coefficients:
if ( m <= 10000 ) % Determined experimentally
    S = leg2cheb_mat( m );
    X = S * X;
else
    X = leg2cheb( X );
end

end

function S = leg2ultra_mat( n )
% LEG2ULTRA_MAT Conversion matrix from Legendre coefficients to C^(3/2).
%
% Given coefficients in the Legendre basis the C^(3/2) coefficients
% can be computed via
%
%     c = rand(10, 1);    % Legendre coefficients
%     S = leg2ultra_mat( length(c) ); % conversion matrix
%     d = S * c;           % C^(3/2) coefficients

% Alex Townsend, 5th May 2016

lam = 1/2;
dg = lam./(lam + (2:n-1))';
v  = [1 ; lam./(lam+1) ; dg];
w  = [0 ; 0 ; -dg];
S  = spdiags( [v w], [0 2], n, n );

end

function S = ultra1mx2leg_mat( n )
% ULTRA1MX2LEG_MAT Conversion matrix for (1-x^2)C^(3/2) to Legendre.
%
% Given coefficients in the (1-x^2)C^(3/2) basis the Legendre coefficients
% can be computed via
%
%     c = rand(10, 1);     % (1-x^2)C^(3/2) coefficients
%     S = ultra1mx2leg_mat( length(c) ); % conversion matrix
%     c_leg = S * c;       % Legendre coefficients
%

% Alex Townsend, 5th May 2016

d = ones(n, 1);
S = spdiags(((1:n).*(2:(n+1))./2./(3/2:n+1/2))', 0, n, n);
S = spdiags( [d,-d], [0,-2], n, n ) * S;

end

function L = cheb2leg_mat( N )
% CHEB2LEG_MAT Construct the cheb2leg conversion matrix.

% This for-loop is a faster and more accurate way of doing:
% Lambda = @(z) exp(gammaln(z+1/2) - gammaln(z+1));
% vals = Lambda( (0:2*N-1)'/2 );
vals = zeros(2*N,1);
vals(1) = sqrt(pi);
vals(2) = 2/vals(1);
for i = 2:2:2*(N-1)
    vals(i+1) = vals(i-1)*(1-1/i);
    vals(i+2) = vals(i)*(1-1/(i+1));
end

L = zeros(N, N);
for j = 0:N-1
    for k = j+2:2:N-1
        L(j+1, k+1) = -k*(j+.5)*(vals((k-j-2)+1)./(k-j)).*(vals((k+j-1)+1)./(j+k+1));
    end
end
c = sqrt(pi)/2;
for j = 1:N-1
    L(j+1, j+1) = c./vals( 2*j+1 );
end
L(1,1) = 1;

end

function M = leg2cheb_mat( N )
% LEG2CHEB_MAT Construct the leg2cheb conversion matrix.

% This for-loop is a faster and more accurate way of doing:
% Lambda = @(z) exp(gammaln(z+1/2) - gammaln(z+1));
% vals = Lambda( (0:2*N-1)'/2 );
vals = zeros(2*N,1);
vals(1) = sqrt(pi);
vals(2) = 2/vals(1);
for i = 2:2:2*(N-1)
    vals(i+1) = vals(i-1)*(1-1/i);
    vals(i+2) = vals(i)*(1-1/(i+1));
end

M = zeros(N, N);
for j = 0:N-1
    for k = j:2:N-1
        M(j+1, k+1) = 2/pi*vals((k-j)+1).*vals((k+j)+1);
    end
end
M(1,:) = .5*M(1,:);

end


text:
classdef math_model_opf_accs_legacy < mp.math_model_opf_accs & mp.mm_shared_opf_legacy
%MP.MATH_MODEL_OPF_ACCS_LEGACY  MATPOWER mathematical model for AC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_ACCS_LEGACY ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        %% constructor
        function obj = math_model_opf_accs_legacy()
            obj@mp.math_model_opf_accs();
            if nargin > 0 && isstruct(mpc)
                obj.mpc = mpc;
            end

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if om.var is empty) in ADD_VAR(), DISPLAY() and
            %%              INIT_INDEXED_NAME(), after object construction,
            %%              but before object use.
        end

        function obj = add_named_set(obj, varargin)
            % call parent method (also checks for valid type for named set)
            add_named_set@mp.math_model_opf_accs(obj, varargin{:});
            obj.add_named_set_legacy(varargin{:});
        end

        function obj = def_set_types(obj)
            obj.def_set_types_legacy();
        end

        function obj = init_set_types(obj)
            init_set_types@mp.math_model_opf_accs(obj);
            obj.init_set_types_legacy();
        end

        function obj = build(obj, nm, dm, mpopt)
            obj.mpc = dm.source;
            build@mp.math_model_opf_accs(obj, nm, dm, mpopt);
            obj.build_legacy(nm, dm, mpopt);
        end

        function obj = add_vars(obj, nm, dm, mpopt)
            add_vars@mp.math_model_opf_accs(obj, nm, dm, mpopt);    %% call parent

            %% legacy user-defined variables
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_vars(nm, dm, mpopt);
            end
        end

        function add_system_costs(obj, nm, dm, mpopt)
            add_system_costs@mp.math_model_opf_accs(obj, nm, dm, mpopt);    %% call parent

            %% legacy user-defined costs
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_costs(nm, dm, 0);
            end
        end

        function obj = add_system_constraints(obj, nm, dm, mpopt)
            %% call parent
            add_system_constraints@mp.math_model_opf_accs(obj, nm, dm, mpopt);

            %% legacy user-defined constraints
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_constraints_ac(nm, dm, mpopt);
            end
        end

        function names = legacy_user_var_names(obj)
            names = {'Vr', 'Vi', 'Pg', 'Qg'};
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [tx, ty] = rectToPolygon(rect)
%RECTTOPOLYGON Convert a rectangle into a polygon (set of vertices)
%
%   POLY = rectToPolygon(RECT);
%   Converts rectangle given as [X0 Y0 W H] or [X0 Y0 W H THETA] into a
%   4-by-2 array double, containing coordinate of rectangle vertices.
%   X0 and Y0 are the coordinates of the "lower left" vertex (before
%   applying rotation), W and H are the width and the height of the
%   rectangle, and THETA is the rotation angle around the first vertex, in
%   degrees.
%
%   See also:
%   orientedBoxToPolygon, ellipseToPolygon, drawRect, drawPolygon
%
%

% ---------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% INRA - TPV URPOI - BIA IMASTE
% created the 06/04/2005.
%

% HISTORY

% extract rectangle parameters
theta = 0;
x0  = rect(1);
y0  = rect(2);
w   = rect(3);
h   = rect(4);
if length(rect) > 4
    theta = rect(5);
end

% precompute angular quantities
cot = cosd(theta);
sit = sind(theta);

% compute vertex coordinates
tx = zeros(4, 1);
ty = zeros(4, 1);
tx(1) = x0;
ty(1) = y0;
tx(2) = x0 + w * cot;
ty(2) = y0 + w * sit;
tx(3) = x0 + w * cot - h * sit;
ty(3) = y0 + w * sit + h * cot;
tx(4) = x0 - h * sit;
ty(4) = y0 + h * cot;

% format output
if nargout <= 1
    tx = [tx ty];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function s=sss(n_id_1,n_id_2,slot_num);

% s=sss(n_id_1,n_id_2,slot_num);
%
% Return the sss for slot slot_num for the specified n_id_1 and n_id_2.
%
% s is of length 62 and only includes the non-zero subcarriers. The calling
% function must insert the DC subcarrier and the 5 zero subcarriers that
% surround the sss.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU 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/>.

error(nargchk(3,3,nargin));
error(chk_param(n_id_1,'n_id_1','scalar','real','integer','>=',0,'<=',167));
error(chk_param(n_id_2,'n_id_2','scalar','real','integer','>=',0,'<=',2));
error(chk_param(slot_num,'slot_num','scalar','real','integer','>=',0,'<=',10));
if ((slot_num~=0)&&(slot_num~=10))
  error('slot_num must be either 0 or 10');
end

% Calculate m0 and m1 from n_id_1
qp=floor(n_id_1/30);
q=floor((n_id_1+qp*(qp+1)/2)/30);
mp=n_id_1+q*(q+1)/2;
m0=mod(mp,31);
m1=mod(m0+floor(mp/31)+1,31);

%s_td=[0 0 0 0 1];
%for t=1:26
%  s_td=[s_td mod(s_td(end-2)+s_td(end-4),2)];
%end
s_td=[0 0 0 0 1 0 0 1 0 1 1 0 0 1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 1 0 1];
s_td=1-2*s_td;

%c_td=[0 0 0 0 1];
%for t=1:26
%  c_td=[c_td mod(c_td(end-1)+c_td(end-4),2)];
%end
c_td=[0 0 0 0 1 0 1 0 1 1 1 0 1 1 0 0 0 1 1 1 1 1 0 0 1 1 0 1 0 0 1];
c_td=1-2*c_td;

%z_td=[0 0 0 0 1];
%for t=1:26
%  z_td=[z_td mod(z_td(end)+z_td(end-2)+z_td(end-3)+z_td(end-4),2)];
%end
z_td=[0 0 0 0 1 1 1 0 0 1 1 0 1 1 1 1 1 0 1 0 0 0 1 0 0 1 0 1 0 1 1];
z_td=1-2*z_td;

s0_m0=s_td(mod(m0:30+m0,31)+1);
s1_m1=s_td(mod(m1:30+m1,31)+1);

c0=c_td(mod(n_id_2:30+n_id_2,31)+1);
c1=c_td(mod(n_id_2+3:30+n_id_2+3,31)+1);

z1_m0=z_td(mod((0:30)+mod(m0,8),31)+1);
z1_m1=z_td(mod((0:30)+mod(m1,8),31)+1);

if (slot_num==0)
  s(2:2:62)=s1_m1.*c1.*z1_m0;
  s(1:2:62)=s0_m0.*c0;
elseif (slot_num==10)
  s(2:2:62)=s0_m0.*c1.*z1_m1;
  s(1:2:62)=s1_m1.*c0;
else
  error('Check code...');
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function RM = fmri_mrestriction(TestType,nH, nC, AC, CC, nHTest)
%
% RM = fmri_mrestriction(TestType, nH, nC, AC, CC, <nHTest>)
%
% Creates a restriction matrix
%
% TestType - String (case insensitive) indicating the test to perform:
%   t  - t-Test, RM will be a row vector
%   tm - t-Test, different p for each delay point (RM same as Fc)
%   Fm - F-Test, different p for each delay point
%   F0 - RM same as t-Test
%   Fd - F-test, different row for each delay, all conditions
%        for each delay on the same row
%   Fc - F-test, different row for each condition, all delays
%        for a condition on the same row
%   Fcd - F-test, different row for each condition for each delay
% nH - total number of elements in HDIR
% nC - total number of stimulus conditions (including fixation)
% AC - list of active conditions
% CC - list of control conditions
% nHTest - list of HDIR components to test (default: all).
%
% Eg: fmri_mrestriction('t',10, 6, [1 3], [2 5], [3 6:9])
% Generates an RM for testing conditions (1+3)-(2+5) using
% components 3,6,7,8, and 9 in the HDIR.  Conditions 0 and 4
% are not tested. There are 10 components in the HDIR, and 5
% stimulus conditions.
%
%
% $Id: fmri_mrestriction.m,v 1.1 2004/03/11 01:30:17 sayres Exp $

% Check for the correct number of arguments %
if(nargin ~= 5 & nargin ~= 6)
  msg = 'USAGE: RM = fmri_mrestriction(TestType, nH, nC, AC, CC, <nHTest>';
  qoe(msg);error(msg);
end

% Check the Test Type %
if( isempty( strmatch(upper(TestType),{'T ','TM','FM','F0','FD','FC','FCD','FDC'},'exact')))
  msg = sprintf('Unkown TestType %s',TestType);
  qoe(msg);error(msg);
end

%% check that the input parameters are correct %%
if(nH < 1)
  msg = 'nH must be greater than 0';
  qoe(msg);error(msg);
end
if(nC < 1)
  msg = 'Number of conditions must be greater than 0';
  qoe(msg);error(msg);
end
if(length(find(AC<0 | AC > nC)))
  msg = 'Invalid condition number in AC';
  qoe(msg);error(msg);
end
if(length(find(CC<0 | CC > nC)))
  msg = 'Invalid condition number in CC';
  qoe(msg);error(msg);
end

% Test that the same condition is not in both AC and CC %
C = [AC CC];
for n = 1:nC,
  if(length(find(C==n)) > 1)
    msg = 'Same condition is in both AC and CC';
    qoe(msg);error(msg);
  end
end

% Set defaults for nHTest, if needed, else test range %
if(nargin == 5) nHTest = [1:nH];
else
  if(length(find(nHTest<1 | nHTest > nH)))
    msg = 'Invalid nHTest Component';
    qoe(msg);error(msg);
  end
end

% Strip Condition 0 from AC (if there) %
iACnz = find(AC ~= 0);
if( ~ isempty(iACnz) )
  AC = AC(iACnz);
end
% Strip Condition 0 from CC (if there) %
iCCnz = find(CC ~= 0);
if( ~ isempty(iCCnz) )
  CC = CC(iCCnz);
end

%% -------- Generate the restriction matrix -------- %%
switch(upper(TestType))
  case {'T','F0'},
    RM = fmri_vrestriction(nH,nC,AC,CC,nHTest);
  case {'FD','TM','FM'},
    for n = 1:length(nHTest),
      RV = fmri_vrestriction(nH,nC,AC,CC,nHTest(n));
      RM(n,:) = reshape(RV', [1 prod(size(RV))]);
    end
  case {'FC'},
    for n = 1:length(iACnz),
      RV = fmri_vrestriction(nH,nC,AC(n),0,nHTest);
      RM(n,:) = reshape(RV', [1 prod(size(RV))]);
    end
    for n = 1:length(iCCnz),
      RV = fmri_vrestriction(nH,nC,0,CC(n),nHTest);
      RM(n+length(AC),:) = reshape(RV', [1 prod(size(RV))]);
    end
  case {'FCD','FDC'},
    RM = [];
    for n = 1:length(iACnz),
      RV = fmri_mrestriction('Fd',nH,nC,AC(n),0,nHTest);
      RM = [RM; RV;];
    end
    for n = 1:length(iCCnz),
      RV = fmri_mrestriction('Fd',nH,nC,0,CC(n),nHTest);
      RM = [RM; RV;];
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [phi,delta,alam,dipdir,ierr] = focal_nd2pl(wanx,wany,wanz,wdx,wdy,wdz)

    % compute strike, dip, rake and dip directions from Cartesian components of the outward normal and slip vectors
    %
    %     usage:
    %     call nd2pl(anx,any,anz,dx,dy,dz,strike,dip,rake,dipdir,ierr)
    %
    %     arguments:
    %     anx,any,anz    components of fault plane outward normal vector in the
    %                    Aki-Richards Cartesian coordinate system (INPUT)
    %     dx,dy,dz       components of slip vector in the Aki-Richards
    %                    Cartesian coordinate system (INPUT)
    %     strike         strike angle in degrees (OUTPUT)
    %     dip            dip angle in degrees (OUTPUT)
    %     rake           rake angle in degrees (OUTPUT)
    %     dipdir         dip direction angle in degrees (OUTPUT)
    %     ierr           error indicator (OUTPUT)
    %
    %     errors:
    %     1              input vectors not perpendicular among each other
    %
    %
    %       call fpsset
    amistr=-360.;
    amastr=360.;
    amidip=0.;
    amadip=90.;
    amirak=-360.;
    amarak=360.;
    amitre=-360.;
    amatre=360.;
    amiplu=0.;
    amaplu=90.;
    orttol=2.;
    ovrtol=0.001;
    tentol=0.0001;
    dtor=0.017453292519943296;
    c360=360.;
    c90=90.;
    c0=0.;
    c1=1.;
    c2=2.;
    c3=3.;

    ierr = 0;
    [ang] = focal_angle(wanx,wany,wanz,wdx,wdy,wdz);
    if (abs(ang - c90) > orttol)
        disp(['ND2PL: input vectors not perpendicular, angle=' num2str(ang)]);
        ierr = 1;
    end
    [anorm,anx,any,anz] = focal_norm(wanx,wany,wanz);
    [dnorm,dx,dy,dz] = focal_norm(wdx,wdy,wdz);
    if (anz > c0)
        [anx,any,anz] = focal_invert(anx,any,anz);
        [dx,dy,dz] = focal_invert(dx,dy,dz);
    end

    if(anz == -c1)
        wdelta = c0;
        wphi = c0;
        walam = atan2(-dy, dx);
    else
        wdelta = acos(-anz);
        wphi = atan2(-anx, any);
        walam = atan2(-dz/sin(wdelta),dx*cos(wphi)+dy*sin(wphi));
    end
    phi = wphi/dtor;
    delta = wdelta/dtor;
    alam = walam/dtor;
    phi = mod(phi+c360, c360);
    dipdir = phi+c90;
    dipdir = mod(dipdir+c360,c360);
end

author: CelsoReyes


text:
%kExpAny 'Eports khoros files into many different format '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ExpAny.pane file
%
% Parameters:
% InputFile: i 'Input File', required: 'First Input data object'
% OutputFile: o 'Output File', required: 'Resulting output data object'
% Toggle: color 'RGB-Color ', default: 0: 'generate a color-image'
% Toggle: nonormalize 'Do not normalize', default: 0: 'if selected, no normalization will be performed. Range is usually 0..255'
%
% Example: o = kExpAny(i, {'i','';'o','';'color',0;'nonormalize',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ExpAny - Eports khoros files into many different format
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kExpAny(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kExpAny(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'color', 0;'nonormalize', 0};
maxval={0,0,0,0};
minval={0,0,0,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'expany"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Laplacian Eigen maps and MDS
    close all; clear all;
	%[X, labels] = generate_data('swiss', 2000);

    % 1st group of points
    X1=[1 1
        1 0.5
        1 1.5
        0.5 1
        1.5 1
        ]
    [N1,r]=size(X1);
    labels([1:N1])=1;  % Labels of the class - 3 classes here

    % 2nd group of points
    X2=[-1 -1
        -1 -0.8
        -1 -1.2
        -0.8 -1
        -1.2 -1
        ]
     [N2,r]=size(X2);
       labels([N1+1:N1+N2])=2;

    % 3rd group of points (overlap with one and two)
%         X3=[-0.9 -1
%         -0.2 -0.8
%         -0.5 -0.5
%          0.5 0.8
%          0.2 0.0
%          0.1 -0.1
%          0.2 -0.2
% %         ]
        X3=X1;
     [N3,r]=size(X3);
   labels([N1+N2+1:N1+N2+N3])=3;

 figure; hold on;
 plot(X1(:,1),X1(:,2),'b*');
 plot(X2(:,1),X2(:,2),'r*');
 plot(X3(:,1),X3(:,2),'k*');
 legend('Class 1', 'Class 2', 'Class 3');
axis([-2 2 -2 2]);

X=[X1' X2' X3']';  %% Construct the entire matrix made of the 3 subgroups of points

[L,r]=size(X);
for i=1:L
 if labels(i)==1 c(i,:)=[0 0 1];
 elseif labels(i)==2; c(i,:)=[1 0 0];
 else c(i,:)=[0 0 0];
 end;
end;
%%
figure; hold on;

%% Perform MDS
   % [mappedX, mapping] = compute_mapping(X, 'MDS', 3);

   disp('Computing kernel matrix...');
   kernel       = 'gauss';
   param1       = 1; % sigma
   param2       = 0; % not used
   K            = gram(X, X, kernel, param1, param2);
   dims_keep    = 2;

   [mappedX, mapping] = cmdscale(K,dims_keep);

   subplot(2,2,1), scatter(mappedX(:,1), mappedX(:,2), 25, c, 'filled'); title(['Result of MDS']); drawnow

%% Perform PCA
    [mappedX, mapping] = compute_mapping(X, 'PCA', 2);
	subplot(2,2,2), scatter(mappedX(:,1), mappedX(:,2), 25, c, 'filled'); title(['Result of ' mapping.name]); drawnow

    %% Perform Laplacian
    [mappedX, mapping] = compute_mapping(X, 'Laplacian', 4);
	subplot(2,2,3), scatter(mappedX(:,1), mappedX(:,2), 25, c, 'filled'); title(['Result of ' mapping.name]); drawnow
    %% Isomap
    [mappedX, mapping] = compute_mapping(X, 'Isomap', 2);
	subplot(2,2,4), scatter(mappedX(:,1), mappedX(:,2), 25, c, 'filled'); title(['Result of ' mapping.name]); drawnow
return;

	figure, scatter3(X(:,1), X(:,2), X(:,3), 5, labels); title('Original dataset'), drawnow
	no_dims = round(intrinsic_dim(X, 'MLE'));
	disp(['MLE estimate of intrinsic dimensionality: ' num2str(no_dims)]);

    % This does not seem to be implemented in the toolbox, -> use matlab
    % dimensionality reduction: http://uk.mathworks.com/help/stats/dimensionality-reduction.html
    [mappedX, mapping] = compute_mapping(X, 'MDS', 2);
	figure, scatter(mappedX(:,1), mappedX(:,2), 5); title(['Result of ' mapping.name]); drawnow

    [mappedX, mapping] = compute_mapping(X, 'LaplacianEigenmaps', no_dims, 7);
	figure, scatter(mappedX(:,1), mappedX(:,2), 5, labels(mapping.conn_comp)); title('Result of Laplacian Eigenmaps'); drawnow

    %% Multi-dimensional scaling from matlab

    % Generate some points in 4-D space, but close to 3-D space, then reduce them to distances only.
    X = [normrnd(0,1,10,3) normrnd(0,.1,10,1)];
    D = pdist(X,'euclidean');

    [Y,e] = cmdscale(D);

    % Four, but fourth one small
    dim = sum(e > eps^(3/4));

    % Poor reconstruction
    maxerr2 = max(abs(pdist(X)-pdist(Y(:,1:2))))

    % Good reconstruction
    maxerr3 = max(abs(pdist(X)-pdist(Y(:,1:3))))

    % Exact reconstruction
    maxerr4 = max(abs(pdist(X)-pdist(Y)))

    %% Try matlab MDS with Swiss role data set

    [X, labels] = generate_data('swiss', 1000);

    % Use the euclidean distance metrix, rows of X should correspond to
    % observations
    D          = pdist(X,'euclidean');

    disp(['D: ' num2str(size(D,1)) ' x ' num2str(size(D,2))]);

    [Y,e]       = cmdscale(D);
  	figure, scatter(Y(:,1), Y(:,2), 5,labels); title(['Result of Matlab MDS']); drawnow

author: epfl-lasa


text:
function [dat] = read_neuralynx_cds(filename, hdr, begsample, endsample, chanindx)

% READ_NEURALYNX_CDS reads selected samples and channels from a combined Neuralynx dataset with separate subdirectories for the LFP, MUA and spike channels
%
% Use as
%    hdr = read_neuralynx_cds(parentdir)
%    dat = read_neuralynx_cds(parentdir, hdr, begsample, endsample, chanindx)
%
% This is not a formal Neuralynx file format, but at the F.C. Donders
% Centre we use it as a directory/file organization in conjunction
% with Neuralynx, SPIKESPLITTING and SPIKEDOWNSAMPLE.

% Copyright (C) 2006, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

needhdr = (nargin==1);
needdat = (nargin>=2);

dirlist = dir(filename);
haslfp   = any(filetype_check_extension({dirlist.name}, 'lfp'));
hasmua   = any(filetype_check_extension({dirlist.name}, 'mua'));
hasspike = any(filetype_check_extension({dirlist.name}, 'spike'));
%hasttl  = any(filetype_check_extension({dirlist.name}, 'ttl'));   % separate file with original Parallel_in
%hastsl  = any(filetype_check_extension({dirlist.name}, 'tsl'));   % separate file with original TimeStampLow
%hastsh  = any(filetype_check_extension({dirlist.name}, 'tsh'));   % separate file with original TimeStampHi

if needhdr
  % read the header from the LFP dataset
  if haslfp
    lfpdir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'lfp'))).name);
    lfphdr = read_header(lfpdir);
    lfpfil = lfphdr.filename;
    for i=1:lfphdr.nChans
      lfplab{i} = sprintf('lfp_%s', lfphdr.label{i});
    end
  else
    lfphdr = [];
    lfpfil = {};
    lfplab = {};
  end

  % read the header from the MUA dataset
  if hasmua
    muadir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'mua'))).name);
    muahdr = read_header(muadir);
    muafil = muahdr.filename;
    for i=1:muahdr.nChans
      mualab{i} = sprintf('mua_%s', muahdr.label{i});
    end
  else
    muahdr = [];
    muafil = {};
    mualab = {};
  end

  % read the header from the SPIKE dataset
  if hasspike
    spikedir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'spike'))).name);
    spikehdr = read_header(spikedir);
    spikefil = spikehdr.filename;
    for i=1:spikehdr.nChans
      spikelab{i} = sprintf('spike_%s', spikehdr.label{i});
    end
  else
    spikehdr = [];
    spikefil = {};
    spikelab = {};
  end

  if haslfp
    % assume that the LFP header describes the MUA and the spikes as well
    hdr = lfphdr;
  else hasmua
    % assume that the MUA header describes the spikes as well
    hdr = muahdr;
  end

  % concatenate the lfp/mua/spike channels, they have already been renamed
  hdr.label    = cat(1, lfplab(:), mualab(:), spikelab(:));
  hdr.filename = cat(1, lfpfil(:), muafil(:), spikefil(:));
  hdr.nChans = length(hdr.label);
  % remember the original header details
  hdr.orig    = {};
  hdr.orig{1} = lfphdr;
  hdr.orig{2} = muahdr;
  hdr.orig{3} = spikehdr;

  % return the header
  dat = hdr;
else

  % use the original header details
  lfphdr   = hdr.orig{1};
  muahdr   = hdr.orig{2};
  spikehdr = hdr.orig{3};
  % the spike header does not contain these, but requires them from the LFP or MUA
  spikehdr.FirstTimeStamp     = hdr.FirstTimeStamp;
  spikehdr.TimeStampPerSample = hdr.TimeStampPerSample;

  if ~isempty(lfphdr)
    Nlfp = lfphdr.nChans;
  else
    Nlfp = 0;
  end

  if ~isempty(muahdr)
    Nmua = muahdr.nChans;
  else
    Nmua = 0;
  end

  if ~isempty(spikehdr)
    Nspike = spikehdr.nChans;
  else
    Nspike = 0;
  end

  % determine which files should be read from each respective dataset
  selchan = zeros(hdr.nChans,1);
  selchan(chanindx) = 1;
  sellfp   = find(selchan(1:Nlfp));
  selmua   = find(selchan((Nlfp+1):(Nlfp+Nmua)));
  selspike = find(selchan((Nlfp+Nmua+1):(Nlfp+Nmua+Nspike)));

  if haslfp
    % the header contains the relevant file names
    lfpdat = read_neuralynx_ds([], lfphdr, begsample, endsample, sellfp);
  else
    lfpdat = [];
  end

  if hasmua
    % the header contains the relevant file names
    muadat = read_neuralynx_ds([], muahdr, begsample, endsample, selmua);
  else
    muadat = [];
  end

  if hasspike
    % the header contains the relevant file names
    spikedat = read_neuralynx_ds([], spikehdr, begsample, endsample, selspike);
  else
    spikedat = [];
  end

  % combine all data into a single array
  dat = cat(1, lfpdat, muadat, spikedat);

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% The COBRAToolbox: testloadBiGGModel.m
%
% Purpose:
%     - reads models from BiGG and the previously downloaded version in
%       the models directory and compares them.
%     - tests one model for equivalence wrt FBA results.
%
% Authors:
%     - Thomas Pfau
%

global CBTDIR
global CBT_MISSING_REQUIREMENTS_ERROR_ID

%Check the requirements (a LP solver is necessary)
solverPkgs = prepareTest('needsLP',true,'requireOneSolverOf',{'gurobi','ibm_cplex','glpk','mosek','quadMinos'},'needsWebAddress','http://bigg.ucsd.edu/api/v2/models');

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testLoadBiGGModel.m'));
cd(fileDir);

% initialize the test
cd([CBTDIR, filesep, 'test', filesep, 'models']);

% Models: FileName of local file, model ID and type for BiGG, along with
% FBA Min and Max value.
modelArr = {'iIT341.xml','iIT341','sbml','BiGGSBML',0,0.692812693473487;...
    'iJO1366.mat','iJO1366','mat','BiGG',NaN,NaN};

% set the tolerance
tol = 1e-6;

tested = false;

%loop through the models
for i = 1:size(modelArr,1)
    %reading the models takes quite a bit of time, so only do it once for
    %all solvers.
    % output a line before launching the test for model i
    fprintf('   Testing %s ...\n', modelArr{i,2});

    % load the model (actually supply the full filename of the path
    % where the model is found)
    model1 = getDistributedModel(modelArr{i,1});
    try
        model2 = loadBiGGModel(modelArr{i,2},modelArr{i,3});
        tested = true;
    catch ME
        if strcmp(ME.identifier,'MATLAB:webservices:Timeout')
            %Could not load the model, skip.
            continue;
        end
    end

    model3 = readCbModel(modelArr{i,2},'fileType',modelArr{i,4});

    %Check that the direct load is the same
    if 1
        printLevel=1;
        [isSame, nDiff, commonFields] = isSameCobraModel(model1, model2, printLevel);
        assert(isSame)
    else
        assert(isSameCobraModel(model1,model2));
    end

    %Check that the model loaded through readCbModel is the same.

    if 1
        printLevel=1;
        [isSame, nDiff, commonFields] = isSameCobraModel(model1, model3, printLevel);
        assert(isSame)
    else
        assert(isSameCobraModel(model1,model3));
    end

    if ~isnan(modelArr{i,5})
        for k = 1:length(solverPkgs.LP)
            fprintf(' -- Running testLoadBiGGModel using the solver interface: %s ... ', solverPkgs.LP{k});
            changeCobraSolver(solverPkgs.LP{k}, 'LP', 0);

            fprintf('   Testing loaded model ... \n');

            % solve the maximisation problem
            FBA = optimizeCbModel(model2, 'max');

            % test the maximisation solution
            assert(FBA.stat == 1);
            assert(abs(FBA.f - modelArr{i,6}) < tol);
            assert(norm(model2.S * FBA.x) < tol);

            % solve the minimisation problem
            FBA = optimizeCbModel(model2, 'min');

            % test the minimisation solution
            assert(FBA.stat == 1);
            assert(abs(FBA.f - modelArr{i,5}) < tol);
            assert(norm(model2.S * FBA.x) < tol);

            % print a line for success of loop i
            fprintf(' Done.\n');
        end
    end
end
%This should only ever happen, if the BiGG db could not be contacted.
if ~tested
    error(CBT_MISSING_REQUIREMENTS_ERROR_ID,'Could not connect to BiGG Database');
end
% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
clear all; close all;
bw=imread('text.png');
se=strel('line', 11, 90);
bw2=imdilate(bw, se);
figure;
subplot(121);  imshow(bw);
subplot(122);  imshow(bw2);

author: HuangCongQing


text:
function [Fh,Vh]=honeyCombMesh(minV,maxV,pointSpacing)

% function [Fh,Vh]=honeyCombMesh(minV,maxV,pointSpacing)
% -----------------------------------------------------------------------
% This function creates the faces (Fh) and vertices (Vh) for a hexagon
% (honey comb) mesh. The hexagons are created between the limits minV
% (containing desired minimum X and Y coordinates) and maxV (containing
% desired maximum X and Y coordinates). The size of the hexagons is set by
% the desired point spacing.
%
% -----------------------------------------------------------------------

%% CREATE TRIANGULATION

maxV(2)=maxV(2)+pointSpacing;
[F,V]=triMeshEquilateral(minV,maxV,pointSpacing);

%% GET DUAL FOR HONEY-COMB

[Vh,Fd]=patch_dual(V,F,0);
numVert=cellfun(@(x) size(x,2),Fd);
Fh=Fd{numVert==6};

[Fh,Vh]=patchCleanUnused(Fh,Vh);

%%
% _*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:
% -------------------------------------------------------------------------
% Author: Jai Juneja
% Date: 12/02/2013
% -------------------------------------------------------------------------
function World = ekfSLAM(handles, AxisDim, Lmks, Wpts, Obstacles)

    %%%%%%%%%%%%%%%%%%%%%%%%% INITIALISATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    cla; % Clear axes
    rob = Robot; % Create new robot object
    sen = Sensor;
    sen_rf = Sensor;
    sen_rf.noise = [0.1; 3 * pi/180];
    sen_rf.range = 10;
    sen.range = 30;
    sen.fov = 270 * pi/180;
    rob.q = [0.1;2*pi/180];

    [World, rob] = configuration(rob, sen_rf, Lmks, Wpts, AxisDim);

    Graphics = initGraphics(World, Obstacles, handles);

    % Determine which obstacles are moving
    numObstacles = length(Obstacles);
    numSegments = 0;
    movingObstacles = zeros(1, numObstacles);
    for i = 1:numObstacles
        if ~isequal([0; 0], Obstacles(i).velocity)
            movingObstacles(i) = 1;
        end
        numSegments = numSegments + size(Obstacles(i).vertices, 2) - 1;
    end
    movingObstacles = find(movingObstacles);

    %%%%%%%%%%%%%%%%%%%%%%%% TEMPORAL LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Initialise some loop variables
    pose_this_scan = [];
    this_scan_good = 0;
    big_turn = 0;
    r = World.r;    % Location of robot pose in mapspace

    for t = 1:World.tend

        World.t = t;

        %%%%%%%%%%%%%%%%%%%%%%%%% SIMULATE WORLD %%%%%%%%%%%%%%%%%%%%%%%%%%
        R_old = rob.R;
        if ~isempty(World.Wpts)
            rob.steer(World.Wpts);
        end
        for i = movingObstacles
            Obstacles(i).move(AxisDim);
        end
        % n = rob.q .* randn(2,1);    % Control noise in real world
        rob.R = rob.move(zeros(2,1), 'true');        % Move with noise
        for lid = 1:size(World.W,2)
            v = sen_rf.noise .* randn(2,1);
            World.y(:,lid) = scanPoint(rob.R, World.W(:,lid)) + v;
        end
        lmks_all = World.y;
        % Determine landmarks that are within the sensor range
        [World.y, lmks_visible] = sen_rf.constrainMeasurement(World.y);

        %%%%%%%%%%%%%%%%%%%%%%%%%%% START EKF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%% EKF PREDICTION STEP %%%%%%%%%%%%%%%%%%%%%%%

        % 1. POSE PREDICTION USING SCAN MATCHING %%%%%%%%%%%%%%%%%%%%%%%%%%
        % TODO: if scan_matching, then do this step. Scan match only if
        % ~isempty(Obstacles)
        enough_correlations = 0;

        % Current stored scan is now the last scan
        last_scan_good = this_scan_good;
        pose_last_scan = pose_this_scan;
        this_scan_good = 0; % Initialise this_scan_good for this iteration

        % Put the previous scan's data in scan_ref
        scan_ref = World.scan_data;
        % Initialise scan data vector at maximum possible size it can take
        World.scan_data = -ones(3, numSegments * (round(sen.fov/sen.angular_res) + 1));
        scan_ndx = 1;
        for i = 1:length(Obstacles)
            scan_data = Obstacles(i).getMeasured(sen, rob);
            % Build scan_data vector
            if ~isempty(scan_data)
                World.scan_data(:, scan_ndx:scan_ndx+size(scan_data, 2)-1) = scan_data;
                scan_ndx = scan_ndx + size(scan_data, 2);
            end
        end
        % Reduce scan to only those appended in above for loop
        World.scan_data = World.scan_data(:, World.scan_data(2, :) ~= -1);

        % Reduce scan so that each laser angle only has one associated
        % measurement
        World.scan_data = removeDuplicateLasers(World.scan_data);
        World.scan_data = World.scan_data(2:3, :); % Remove index row

        % First add Gaussian white noise to scan
        v = repmat(sen.noise, 1, size(World.scan_data, 2)) .* randn(2,size(World.scan_data, 2));
        World.scan_data = World.scan_data + v;

        if ~isempty(World.scan_data)
            obstaclesDetected = 1;

            % Convert to Cartesian co-ordinates for scan matching
            World.scan_data = getInvMeasurement(World.scan_data);

            % If the number of point scanned exceeds the specified
            % tolerance, it is suitable for scan matching
            if size(World.scan_data, 2) > World.scan_corr_tolerance
                this_scan_good = 1;
            end

            % Convert to global co-ordinates for graphical display
            % World.scan_global = transToGlobal(rob.R, World.scan_data);
        else
            obstaclesDetected = 0;
            World.scan_global = [];
        end
        % If both last and current scans are good for scan matching,
        % proceed to match scans:
        if last_scan_good && this_scan_good
            n = rob.q .* randn(2,1); % Noise in odometry measurement

            % Do ICP scan matching using odometry input as initial guess
            [R, T, correl, icp_var] = doICP(scan_ref, World.scan_data, 1, rob.u + n);
            % If there are enough points correlated between the two
            % scans:
            if length(correl) > World.scan_corr_tolerance
                enough_correlations = 1;
                % Determine estimated pose from scan match
                r_scan = zeros(3, 1);
                da = getPiToPi(asin(R(2,1)));
                r_scan(1:2) = transToGlobal(pose_last_scan, T);
                r_scan(3) = pose_last_scan(3) + da;

                % Transform scan data to global co-ordinates
                scan_data_corr = World.scan_data;
                scan_data_corr = transToGlobal(r_scan, scan_data_corr);
                % Grab data points that were successfully corellated
                scan_data_corr = scan_data_corr(:, correl(:, 2)');

                [r_odo, R_r, R_n] = rob.move(n);

                % Compute covariance matrix of scan match
                C = getICPCovariance(icp_var, scan_data_corr);
                J_u = [R zeros(2, 1); 0 0 1];
                cov_scan = J_u * C * J_u';
                cov_scan_norm = trace(cov_scan);
            end
        end

        % 2. POSE PREDICTION USING ODOMETRY MEASUREMENTS %%%%%%%%%%%%%%%%%%

        % If there isn't a useful scan, just use odometry data
        if ~enough_correlations
            n = rob.q .* randn(2,1);
            [r_odo, R_r, R_n] = rob.move(n);
            cov_scan = zeros(3, 3);
            dr_scan = zeros(3, 1);
        else
            dr_scan = r_scan - rob.r;
            dr_scan(3) = getPiToPi(dr_scan(3));
        end

        % 3. WEIGHTAGE OF ODOMETRY AND SCAN MATCH PREDICTIONS %%%%%%%%%%%%%

        cov_odo = R_n * World.Q * R_n';
        cov_odo_norm = trace(cov_odo);

        dr_odo = r_odo - rob.r;
        dr_odo(3) = getPiToPi(dr_odo(3));

        dr_true = rob.R - R_old;

        % Check if the robot is making a large turn
        if abs(dr_scan(2)) > pi/6 || abs(dr_odo(2)) > pi/6
            big_turn = 1;
        elseif abs(dr_odo(2)) < pi/18
            big_turn = 0;
        end

        % Determine weights:
        % If not enough correlations, or turning is large, use odometry
        if ~enough_correlations % || big_turn
            weight_odo = 1;
            weight_scan = 0;
        else
            weight_odo = cov_scan_norm / (cov_scan_norm + cov_odo_norm);
            weight_scan = 1 - weight_odo;
            weight_scan = 1;
            weight_odo = 0;
        end

        % Determine robot pose using weights:
        rob.r = weight_odo * dr_odo + weight_scan * dr_scan + rob.r;
        World.x(r) = rob.r;

        % Covariance update
        % Caution: this method of update is suboptimal for CPU processing.
        % Alternative is to have a single equation that updates both
        % the pose and landmark covariances in a single step, but this
        % equation is rather complicated.
        P_rr = World.P(r,r);
        World.P(r,:) = R_r * World.P(r,:);
        World.P(:,r) = World.P(r,:)';
        World.P(r,r) = R_r * P_rr * R_r' + ...
            weight_scan * cov_scan + ...
            weight_odo * cov_odo;

        %%%%%%%%%%%%%%%%%%%%%% EKF UPDATE STEP %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 1. CORRECT KNOWN VISIBLE LANDMARKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Find the visible landmarks. The below code circumvents the data
        % association problem:
        r_old = rob.r;

        lmks_visible_known = find(World.l(1,lmks_visible));
        lids = lmks_visible(lmks_visible_known);

        for lid = lids
            % Measurement prediction
            [e, E_r, E_l] = scanPoint(rob.r, World.x(World.l(:,lid)));

            E_rl = [E_r E_l];
            rl   = [r World.l(:,lid)'];
            E    = E_rl * World.P(rl,rl) * E_rl';

            % Actual Measurement
            yi = World.y(:,lid);

            % Innovation
            z = yi - e;
            z(2) = getPiToPi(z(2));
            Z = World.M + E;

            % Kalman gain
            K = World.P(:, rl) * E_rl' * Z^-1;

            % Update state and covariance
            World.x = World.x + K * z;
            World.P = World.P - K * Z * K'; % The complexity of this line is
                                            % very high (subtraction of two
                                            % large matrices)
            rob.r = World.x(r);
        end

        % 2. INITIALISE NEW LANDMARKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Check landmark availabiliy.
        % Again data association is avoided. A new landmark is initialised
        % in the storage space reserved specifically for that landmark:
        lmks_visible_unknown = find(World.l(1,lmks_visible)==0);
        lids = lmks_visible(lmks_visible_unknown);

        if ~isempty(lids)
            for lid = lids
                s = find(World.mapspace, 2);
                if ~isempty(s)
                    World.mapspace(s) = 0;
                    World.l(:,lid) = s';
                    % Measurement
                    yi = World.y(:,lid);

                    [World.x(World.l(:,lid)), L_r, L_y] = invScanPoint(rob.r, yi);
                    World.P(s,:) = L_r * World.P(r,:);
                    World.P(:,s) = World.P(s,:)'; % Cross variance of lmk with all other lmks
                    World.P(s,s) = L_r * World.P(r,r) * L_r' + L_y * World.M * L_y';
                end
            end
        end
        pose_this_scan = rob.r; % The corrected pose at which the scan was taken

        %%%%%%%%%%%%%%%%%%%%%%%%%% MAPPING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Determine pose correction during update step
        update = rob.r - r_old;

        if enough_correlations && ~isequal(rob.r, r_old)
            % Adjust latest correlated scan points to account for pose correction
            World.scan_data = transformPoints(World.scan_data, update);

            World = doMapping(World, sen, rob.r, World.scan_data, handles);
        elseif ~obstaclesDetected && isequal(mod(World.t, 5), 0)
            % If no obstacles have been detected, only compute map every 5
            % iterations
            World = doMapping(World, sen, rob.r, World.scan_data, handles);
        end

        if ~isempty(World.scan_data)
            World.scan_global = transToGlobal(rob.r, World.scan_data);
        end

        %%%%%%%%%%%%%%%%%%%%%%% HISTORICAL DATA COLECTION %%%%%%%%%%%%%%%%%
        World.R_hist(:, t) = rob.R(1:2);	% True position history
        World.r_hist(:, t) = rob.r(1:2);   % Estimated position history
        pose_error2 = (rob.R(1:2) - rob.r(1:2)).^2;
        World.error_hist(t) = sqrt(pose_error2(1) + pose_error2(2));
        World.Pr_hist(:, t) = [World.P(r(1),r(1)); World.P(r(2),r(2))];

        if enough_correlations
            scan_error = abs(dr_scan-dr_true); scan_error(3) = abs(getPiToPi(scan_error(3)));
            scan_error(1) = pdist([0 0; scan_error(1) scan_error(2)]); scan_error(2) = [];
        else
            scan_error = [0; 0];
        end

        odo_error = abs(dr_odo-dr_true); odo_error(3) = abs(getPiToPi(odo_error(3)));
        odo_error(1) = pdist([0 0; odo_error(1) odo_error(2)]); odo_error(2) = [];

        World.scan_error_hist(:, t) = scan_error;
        World.odo_error_hist(:, t) = odo_error;
        World.turning_hist(t) = rob.u(2);
        World.weight_scan_hist(t) = weight_scan;
        World.weight_odo_hist(t) = weight_odo;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%% GRAPHICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Graphics.lmks_all = lmks_all; Graphics.lmks_visible = lmks_visible;
        doGraphics(rob, World, Graphics, Obstacles(movingObstacles), AxisDim);

    end

    %%%%%%%%%%%%%%%%%%%%%%%% PLOT HISTORICAL DATA %%%%%%%%%%%%%%%%%%%%%%%%%
    % Plot trajectories;
    set(Graphics.R_hist, 'xdata', World.R_hist(1,:), 'ydata', World.R_hist(2, :))
    set(Graphics.r_hist, 'xdata', World.r_hist(1,:), 'ydata', World.r_hist(2, :))

    % Plot position uncertainties
    figure('color', 'white')
    subplot(2, 1, 1)
    plot(1:World.tend, World.error_hist)
    title('EKF SLAM: position error over time')
    xlabel('Time (s)'), ylabel('Position Error (m)')
    subplot(2, 1, 2)
    plot(1:World.tend, sqrt(World.Pr_hist(1, :)), 'r', ...
        1:World.tend, sqrt(World.Pr_hist(2, :)), 'b')
    title('EKF SLAM: position uncertainty over time')
    xlabel('Time (s)'), ylabel('Standard Deviation (m)')
    legend('St. Dev. in x', 'St. Dev. in y')

    figure('color', 'white')
    subplot(4, 1, 1)
    AX = plotyy(1:World.tend, World.scan_error_hist(1,:), ...
        1:World.tend, World.scan_error_hist(2,:), 'plot');
    set(get(AX(1),'Ylabel'),'String',['Position' char(10) 'Error (m)'])
    set(get(AX(2),'Ylabel'),'String',['Angular' char(10) 'Error (rad)'])
    title('Scan errors over time')
    xlabel('Time (s)')

    subplot(4, 1, 2)
    AX2 = plotyy(1:World.tend, World.odo_error_hist(1,:), ...
        1:World.tend, World.odo_error_hist(2,:), 'plot');
    set(get(AX2(1),'Ylabel'),'String',['Position' char(10) 'Error (m)'])
    set(get(AX2(2),'Ylabel'),'String',['Angular' char(10) 'Error (rad)'])
    title('Odometry error over time')

    subplot(4, 1, 3)
    plot(1:World.tend, World.weight_scan_hist, 'r', ...
        1:World.tend, World.weight_odo_hist, 'b')
    title('Scan and odometry weights over time')
    xlabel('Time (s)'), ylabel('Weight')
    legend('Scan', 'Odometry')
    axis tight

    subplot(4, 1, 4)
    plot(1:World.tend, World.turning_hist)
    title('Turning control over time')
    xlabel('Time (s)'), ylabel('Angle (rad)')
    axis tight
end

author: jaijuneja

repo: ekf-slam-matlab


text:
classdef ShFunc_Chomog_EnforceCh < ShFunc_Chomog

    properties (Access = public)
        ChTarget
    end

    properties (Access = private)
       C0
       weights
       difC
       difCp
    end

    properties (Access = protected)
        pNorm
    end

    methods (Access = public)

      function computeGradientValue(obj)
            obj.computeChDerivative();
            nElem  = obj.getnElem();
            nGaus  = obj.getnGaus();
            nStres = obj.getnStre();
            nComp  = obj.computeNcomp(nStres);
            grad = zeros(nElem,nGaus,nComp);
            dCh = permute(obj.dCh,[3,4,1,2]);
            for iStres = 1:nStres
                for jStres = 1:nStres
                    iv     = obj.vector2Voigt(iStres,jStres);
                    difCij = obj.difC(iStres,jStres);
                    C0ij   = obj.C0(iStres,jStres);
                    wij    = obj.weights(iv);
                    p      = obj.pNorm;
                    dChij   = squeeze(dCh(:,:,iStres,jStres));
                    dCostIv = (sum(obj.difCp))^(1/p-1)*wij*(difCij/C0ij)^(p-1)*(dChij/C0ij);
                    grad(:,:,iv) = dCostIv + grad(:,:,iv);
                end
            end
            obj.gradient = sum(grad,3);
        end

        function computeFunctionValue(obj)
            obj.difC  = obj.computeDifC();
            obj.difCp = obj.computeNormP(obj.difC);
            p = obj.pNorm;
            obj.value = sum(obj.difCp)^(1/p);
        end

        function v = getVariablesToPlot(obj)
            ChTargetP = obj.computeNormP(obj.ChTarget);
            p = obj.pNorm;
            normChTarget = sum(ChTargetP)^(1/p);
            v{1} = (obj.value*obj.value0)/normChTarget;
        end

    end

    methods (Access = protected)

        function difCp = computeNormP(obj,difC)
            nStres = obj.getnStre();
            nComp  = obj.computeNcomp(nStres);
            dCp = zeros(nComp,1);
            for iStres = 1:nStres
                for jStres = 1:nStres
                    iv = obj.vector2Voigt(iStres,jStres);
                    difCij = difC(iStres,jStres);
                    C0ij   = obj.C0(iStres,jStres);
                    wij    = obj.weights(iv);
                    p      = obj.pNorm;
                    costij = wij*(difCij/C0ij)^p;
                    dCp(iv) = costij;
                end
            end
            difCp = dCp;
        end

        function difC = computeDifC(obj)
            difC = (obj.Chomog - obj.ChTarget);
        end

        function computeChTarget(obj,cParams)
            obj.ChTarget = ChTargetFactory.create(cParams);
        end

        function computeC0(obj)
            C = ones(3,3);
            C(1,1) = obj.ChTarget(1,1);
%             C(2,2) = obj.ChTarget(2,2);
%             C(3,3) = obj.ChTarget(3,3);
             C(1,2) = 0.01;
             C(2,1) = 0.01;
            obj.C0 = C;
        end

        function computeWeights(obj)
            f = 1;%0000;
            obj.weights = [1,1,1,1,1,f]';
        end

    end

    methods (Access = private, Static)

        function n = computeNcomp(nStre)
           n = (nStre+1)*nStre/2;
        end

        function [iv] = vector2Voigt(iStre,jStre)
            T = [1 6 5; 6 2 4; 5 4 3];
            iv = T(iStre,jStre);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ rexp, sexp ] = poly_t6 ( )

%*****************************************************************************80
%
%% POLY_T6 returns the monomials associated with a 6 node triangle.
%
%  Reference Element T6:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  6  5
%    |  |   \
%    |  |    \
%    0  1--4--2
%    |
%    +--0--R--1-->
%
%  Formula:
%
%    Given coefficients A(I), the polynomial interpolant at (R,S) is
%
%      P(R,S) = sum ( 1 <= I <= N ) A(I) * R**REXP(I) * S**SEXP(I)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer REXP(6), SEXP(6), the powers of R and S associated
%    with each monomial.
%
  rexp(1:6) = [ 0, 0, 1, 0, 1, 2 ];
  sexp(1:6) = [ 0, 1, 0, 2, 1, 0 ];

  return
end

author: johannesgerer

repo: jburkardt-m