text:
function reducedModel = removeGenes(model,genesToRemove,removeUnusedMets,removeBlockedRxns,standardizeRules)
% removeGenes
%   Deletes a set of genes from a model
%
%   model                   a model structure
%   genesToRemove           either a cell array of gene IDs, a logical vector
%                           with the same number of elements as genes in the model,
%                           or a vector of indexes to remove
%   removeUnusedMets        remove metabolites that are no longer in use (opt, default
%                           false)
%   removeBlockedRxns       remove reactions that get blocked after deleting the genes
%                           (opt, default false)
%   standardizeRules        format gene rules to be compliant with standard format
%                           (opt, default true)
%
%   reducedModel            an updated model structure
%
%   Usage: reducedModel = removeGenes(model,genesToRemove,removeUnusedMets,removeBlockedRxns)

if nargin<3
    removeUnusedMets = false;
end
if nargin<4
    removeBlockedRxns = false;
end
if nargin<5
    standardizeRules = true;
end
%Format grRules and rxnGeneMatrix:
if standardizeRules
    [grRules,rxnGeneMat,toCheck] = standardizeGrRules(model,true);
    model.grRules    = grRules;
    model.rxnGeneMat = rxnGeneMat;
else
    toCheck    = [];
    rxnGeneMat = model.rxnGeneMat;
end
reducedModel = model;
%Only remove genes that are actually in the model
if ~(islogical(genesToRemove) || isnumeric(genesToRemove))
    genesToRemove=convertCharArray(genesToRemove);
    genesToRemove=genesToRemove(ismember(genesToRemove,model.genes));
end
if ~isempty(genesToRemove)
    indexesToRemove = getIndexes(model,genesToRemove,'genes');
    if ~isempty(indexesToRemove)
        %Make 0 corresponding columns from rxnGeneMat:
        reducedModel.rxnGeneMat(:,indexesToRemove) = 0;
        genes = model.genes(indexesToRemove);
        %Check if conflicting grRules contain any of the genes to remove
        if ~isempty(toCheck)
            for i=1:numel(toCheck)
                index   = toCheck(i);
                gIdxs   = find(rxnGeneMat(index,:));
                g2check = model.genes(gIdxs);
                if any(ismember(g2check,genes))
                    warning(['Conflicting grRule #' num2str(index) ' (' model.rxns{index} ') contains at least one of the genes to be removed, this grRule will be bypassed in order to avoid logical errors'])
                end
            end
        end
        canCarryFlux = true(size(model.rxns));
        %Loop through genes and adapt rxns:
        for i = 1:length(genes)
            %Find all reactions for this gene and loop through them:
            geneRxns = find(rxnGeneMat(:,indexesToRemove(i)));
            if ~isempty(geneRxns)
                for j = 1:numel(geneRxns)
                    index  = geneRxns(j);
                    grRule = reducedModel.grRules{index};
                    ruleGenes = reducedModel.genes(logical(rxnGeneMat(index,:)));
                    if ~ismember(index,toCheck) && canCarryFlux(index) && ~isempty(grRule)
                        %Check if rxn can carry flux without this gene:
                        canCarryFlux(index) = canRxnCarryFlux(ruleGenes,grRule,genes{i});
                        %Adapt gene rule & gene matrix:
                        grRule = removeGeneFromRule(grRule,genes{i});
                        reducedModel.grRules{index} = grRule;
                    end
                end
            end
        end
        %Delete or block the reactions that cannot carry flux:
        if removeBlockedRxns
            rxnsToRemove = reducedModel.rxns(~canCarryFlux);
            reducedModel = removeReactions(reducedModel,rxnsToRemove,removeUnusedMets,true);
        else
            reducedModel = removeReactions(reducedModel,[],removeUnusedMets,true);
            reducedModel.lb(~canCarryFlux) = 0;
            reducedModel.ub(~canCarryFlux) = 0;
        end
    end
end
%Format grRules and rxnGeneMatrix after all modifications
if standardizeRules
    [grRules,rxnGeneMat] = standardizeGrRules(reducedModel,true);
    reducedModel.grRules = grRules;
    reducedModel.rxnGeneMat = rxnGeneMat;
end
end

function canIt = canRxnCarryFlux(ruleGenes,geneRule,geneToRemove)
%This function converts a gene rule to a logical statement, and then
%asseses if the rule is true (i.e. rxn can still carry flux) or not (cannot
%carry flux).
geneRule = [' ', geneRule, ' '];
for i = 1:length(ruleGenes)
    if strcmp(ruleGenes{i},geneToRemove)
        geneRule = strrep(geneRule,[' ' ruleGenes{i} ' '],' false ');
        geneRule = strrep(geneRule,['(' ruleGenes{i} ' '],'(false ');
        geneRule = strrep(geneRule,[' ' ruleGenes{i} ')'],' false)');
    else
        geneRule = strrep(geneRule,[' ' ruleGenes{i} ' '],' true ');
        geneRule = strrep(geneRule,['(' ruleGenes{i} ' '],'(true ');
        geneRule = strrep(geneRule,[' ' ruleGenes{i} ')'],' true)');
    end
end
geneRule = strtrim(geneRule);
geneRule = strrep(geneRule,'and','&&');
geneRule = strrep(geneRule,'or','||');
canIt    = eval(geneRule);
end

function geneRule = removeGeneFromRule(geneRule,geneToRemove)
%This function receives a standard gene rule and it returns it without the
%chosen gene.
geneSets = strsplit(geneRule,' or ');
hasGene  = ~cellfun(@isempty,strfind(geneSets,geneToRemove));
geneSets = geneSets(~hasGene);
geneRule = strjoin(geneSets,' or ');
end

author: SysBioChalmers

repo: RAVEN


text:
function [dWx,dWy] = vl_dwaffine(x,y)
% VL_DWAFFINE  Derivative of an affine warp
%   [DWX,DWY]=VL_DWAFFINE(X,Y) returns the derivative of the 2-D affine
%   warp [WX; WY] = [A T] [X; Y] with respect to the parameters A,T
%   computed at points X,Y.
%
%   See also: VL_WAFFINE(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% dW = [ kron(x',I) I ]
%    |
%    = [ x1  0  x2  0 1 0 ]
%      [  0 x1   0 x2 0 1 ]

z = zeros(length(x(:)),1) ;
o =  ones(length(x(:)),1) ;

dWx = [ x(:) z      y(:) z      o z ] ;
dWy = [ z    x(:)   z    y(:)   z o ] ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% SP_H1_ERROR: Evaluate the error in H^1 norm.
%
%   [errh1, errl2, errh1s] = sp_h1_error (space, msh, u, uex, graduex)
%
% INPUT:
%
%    space:   object defining the space of discrete functions (see sp_vector)
%    msh:     object defining the domain partition and the quadrature rule (see msh_cartesian)
%    u:       vector of dof weights
%    uex:     function handle to evaluate the exact solution
%    graduex: function handle to evaluate the gradient of the exact solution
%
% OUTPUT:
%
%     errh1:  error in H^1 norm
%     errl2:  error in L^2 norm
%     errh1s: error in H^1 seminorm
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2011, 2015 Rafael Vazquez
%
% 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [errh1, errl2, errh1s] = sp_h1_error (space, msh, u, uex, graduex)

  if (numel(u) ~= space.ndof)
    error ('Wrong size of the vector of degrees of freedom')
  end

  errl2 = 0; errh1s = 0;
  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    sp_col  = sp_evaluate_col (space, msh_col, 'value', true, 'gradient', true);

    [~, err_l2, err_h1s] = sp_h1_error (sp_col, msh_col, u, uex, graduex);

    errh1s = errh1s + err_h1s.^2;
    errl2 = errl2 + err_l2.^2;
  end

  errh1 = sqrt (errl2 + errh1s);
  errl2 = sqrt (errl2);
  errh1s = sqrt (errh1s);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function yPred = objFcn(p , tObs, drug)
%Copyright (c) 2011, The MathWorks, Inc.

L0  = p(1)   ; % Drug-independent parameter 1
L1  = p(2)   ; % Drug-independent parameter 2
k1  = p(3)   ; % Drug-independent parameter 3

k2_A  = p(4) ; % Drug dependent parameter (drug A)
k2_B  = p(5) ; % Drug dependent parameter (drug B)

% Simulate model for drug A
yPred_A = evalTumorWeight(tObs(drug == 'A'), [L1, L0, k1, k2_A]) ;

% Simulate model for drug B
yPred_B = evalTumorWeight(tObs(drug == 'B'), [L1, L0, k1, k2_B]) ;

% Combine prediction
yPred = [yPred_A; yPred_B] ;

end

author: Sable

repo: mcbench-benchmarks


text:
function [] = isotopomerViewer(mdv1, mdv2, names)
% Takes in an "experiment" and views the isotopomer as distributions
% between `mdv1` and `mdv2`. No output.
%
% USAGE:
%
%    isotopomerViewer(mdv1, mdv2, names)
%
% INPUTS:
%    mdv1, mdv2:    structures
%    names:         names in the plot

bins = round(sqrt(size(mdv1,2)));

for i = 1:size(mdv1,1)
    subplot(4,4,mod(i-1,16)+1);
    [x1,x2]=hist(mdv1(i,:), bins);
    plot(x2,x1);
    if max(x2) < .02
        plot(x2,x1, 'k');
    end
    hold on;
    [x1,x2]=hist(mdv2(i,:), bins);
    plot(x2,x1, 'g');
    if max(x2) < .02
        plot(x2,x1, 'k');
    end
    title(names{i});
    hold off;
    if mod(i,16) == 0
        pause;
    end
end

author: opencobra

repo: cobratoolbox


text:
function X = real(X)
% REAL (overloaded)

X.basis = real(X.basis);
X.conicinfo = [0 0];
X = clean(X);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [u,p] = StokesLSCDGS(u,p,f,g,A,B,auxMat,elem,smootherOpt,Ai,Si,SSi,Res,Pro)
%% STOKESLSCDGSTRI Matrix-vecotr form DGS relaxation for Stokes eqns
%
%************************** Algorithm Description *************************
%
%  In matrix form, we need to solve the equations
%
%         Lx =  |A B' | |u| = |f|
%               |B -D | |p| = |0|
%  Distributive matrix will be given by
%
%         M = |I            B'|
%             |0 -inv(BB')BAB'|,
%
%  Therefore, the transformed matrix T will be
%
%         T = L*M = |A  PAB'               |    with P = I-B'inv(BB')B
%                   |B  BB'+D*inv(BB')BAB' |
%
%         tildeT = |Su   0|    with Su and Sp smoother for A and BB'.
%                  |B   Sp|
%
%  The matrix form DGS update can be written as
%
%    |uk+1|   |uk|                 |ru|     with ru = f-Auk -B'pk,
%    |    | = |  | + M*inv(tildeT)*|  |
%    |pk+1|   |pk|                 |rp|     and  rp = 0-Buk
%
%********************** End-of- Algorithm Description *********************
%
%  Created by Ming Wang (with discussion with Long Chen) at Jan, 2012.

%% Set up of smoothers
smoothingstep = smootherOpt.smoothingstep;
smootherSp = upper(smootherOpt.smootherSp);
smootherbarSp = upper(smootherOpt.smootherbarSp);
if isfield(smootherOpt,'smootherbarSpPara')
    smootherbarSpPara = smootherOpt.smootherbarSpPara;
else
    smootherbarSpPara = 1;
end
if strcmp(smootherbarSp,'VCYCLE') || strcmp(smootherSp,'VCYCLE')
    % mg Parameters for barSp
    optionmg.solvermaxit = 1;
%     optionmg.tol = 0.1;
    optionmg.solver = 'VCYCLE';
    optionmg.smoothingstep = 2;
    optionmg.printlevel=0;
    optionmg.setupflag = 0;
end

%% initialize smoother
Bt = auxMat.Bt;
BBt = auxMat.BBt;
BABt = auxMat.BABt;
Su = auxMat.Su;
Sp = auxMat.Sp;
Spt = auxMat.Spt;
DSp = auxMat.DSp;

%% DGS relaxation step
for k = 1: smoothingstep
    % Step 1: relax Momentum eqns
    u = u + Su\(f-Bt*p-A*u);
    % Step 2: relax transformed Continuity eqns
    rp = g - B*u;
    switch(smootherSp)
        case 'SGS'
            dq = Spt\(DSp.*(Sp\rp)); % symmetric Gauss-Seidel iteration
        case 'GS'
            dq = Sp\rp; % Gauss-Seidel iteration
        case 'VCYCLE'
            if exist('elem','var')
                dq = mg(BBt,rp,elem,optionmg,Ai,Si,SSi,Res,Pro); % mg vcycle
            else
                dq = amg(BBt,rp,optionmg); % mg vcycle
            end
    end
    % Step 3: transform the correction back to the original variables
    %   Step 3.1: velocity
    u = u + Bt*dq;
    %   Step 3.2: pressure
    dq = BABt*dq;
%     dq = dq - mean(dq);
    switch(smootherbarSp)
        case 'GS'
            dp = Spt\dq;
        case 'SGS'
            dp = Sp\(DSp.*(Spt\dq));
%             dp = dp + Sp\(DSp.*(Spt\(dq - BBt*dp)));
        case 'VCYCLE'
%            dp = BBt\dq;
            if exist('elem','var')
%                 dq =dq-mean(dq);
                dp = mg(BBt,dq,elem,optionmg,Ai,Si,SSi,Res,Pro); % mg vcycle;
            else
                dp = amg(BBt,dq,optionmg);
            end
    end
    p = p - smootherbarSpPara*dp;

end

repo: Aero-matlab


text:
function [xyZ,sigma,gammaVal]=ellips2PolarStereo(latLon,lambda0,k0,xyPole,a,f)
%%ELLIPS2POLARSTEREO Convert a point on the reference ellipsoid in
%            ellipsoidal (geodetic) coordinates gives as a latitude and a
%            longitude into a point on the polar stereographic projection.
%            Stereographic projections tend to be used near the poles. The
%            Uniform Polar Stereographic (UPS) coordinate system, which is
%            a specific realization of the polar stereographic coordinate
%            system, is typically only used at latitudes >=84 degrees and
%            those < -80 degrees.
%
%INPUTS: latLon A 2XN set of N [latitude;longitude] pairs given in radians
%               to be converted into polar stereographic coordinates.
%       lambda0 The location in raidns of the central meridian. If this
%               parameter is omitted or an empty matrix is passed, the
%               default of 0 is used, which is the value used in the UPS
%               system.
%            k0 The unitless central scale factor. If this parameter is
%               omitted or an empty matrix is passed, then the default of
%               0.994 is used, which is the value used in the UPS.
%        xyPole A 2X1 vector holding the false Easting and False
%               Northing of the origin. This is typically given in meters.
%               If this parameter is omitted or an empty matrix is passed,
%               then the default of [2000000;2000000] is used, which is the
%               value in meters used for the UPS system.
%             a The semi-major axis of the reference ellipsoid. If this
%               argument is omitted, the value in
%               Constants.WGS84SemiMajorAxis is used.
%             f The flattening factor of the reference ellipsoid. If this
%               argument is omitted, the value in Constants.WGS84Flattening
%               is used.
%
%OUTPUTS: xyZ A 3XN set of points in polar stereographic coordinates
%             consisting of [Easting;Northing;Z], where Z=sign(latitude),
%             so 1 for the northern hemisphere and -1 for the southern
%             hemisphere and 0 on the equator. Easting and Northing will be
%             in meters if a and xyPole are in meters.
%       sigma The NX1 set of point scales. This indicates how the map
%             projections enlarges or reduces small distances in this type
%             of map projection.
%    gammaVal The NX1 set of convergence of meridians in radians. This
%             gives the angles of intersection between the meridian and
%             vertical lines on the map projection that are given by
%             x=constant value. This value
%
%If lambda0, k0, xyPole, a and f are all omitted, then the conversion is
%for the UPS on the WGS-84 reference ellipsoid.
%
%Chapter 21 of [1] discusses the polar stereographic projection. The
%implementation is taken from Sections 8 and 9 of [1]. Section 10 of 2
%provides the values for the UPS.
%
%REFERENCES:
%[1] J. P. Snyder, "Map projections- a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%[2] Office of Geomatics, "National geospatial-intelligence agency
%    standardization document: Implementation practice: The universal
%    grids and the transverse mercator and polar stereographic map
%    projections," National Geospatial-Intelligence Agency, Tech. Rep.
%    NGA.SIG.0012 2.0.0 UTMUPS, 25 Mar. 2014. [Online]. Available:
%    http://earth-info.nga.mil/GandG/publications/NGA_SIG_0012_2_0_0_UTMUPS/NGA.SIG.0012_2.0.0_UTMUPS.pdf
%
%July 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(lambda0))
   lambda0=0;%Central meridian location for UPS.
end

if(nargin<3||isempty(k0))
   k0=0.994;%Central scale factor for UPS
end

if(nargin<4||isempty(xyPole))
   xyPole=[2000000;2000000];%Easting and Northing of the pole for UPS.
end

if(nargin<5||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

if(nargin<6||isempty(f))
    f=Constants.WGS84Flattening;
end

%The first numerical eccentricity of the ellipsoid.
e=sqrt(2*f-f^2);

N=size(latLon,2);

xyZ=zeros(3,N);
sigma=zeros(N,1);
gammaVal=zeros(N,1);
for curPoint=1:N
    phi=latLon(1,curPoint);
    lambda=latLon(2,curPoint);
    lambda=lambda-lambda0;

    Z=sign(phi);
    if(Z==-1)
        phi=-phi;
    end

    sinPhi=sin(phi);
    sinLambda=sin(lambda);
    cosLambda=cos(lambda);

    [sinChi,cosChi]=ellipsLat2SinCosConformLat(phi,f);

    k90=sqrt(1-e^2)*exp(e*atanh(e));
    denom=k90*(1+sinChi);
    x=2*k0*a*sinLambda*cosChi/denom;
    y=-2*k0*a*cosLambda*cosChi/denom;
    sigma(curPoint)=k0*2*sqrt(1-e^2*sinPhi^2)*exp(e*atanh(e*sinPhi))/(k90*(1+sinPhi));
    gammaVal(curPoint)=lambda;

    if(Z==-1)
        y=-y;
        gammaVal(curPoint)=-gammaVal(curPoint);
    end

    xyZ(:,curPoint)=[x;y;Z]+[xyPole;0];
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 Brockett_Wen_bfgs_test()

    n = 500;
    p = 4;
    B = randn(n, n); B = B + B';
    D = sparse(diag(p : -1 : 1));

    M = stiefelfactory(n, p);

    problem.M = M;
    problem.cost  = @cost;
    problem.egrad = @egrad;

    % Cost function
    function f = cost(X)
        f = trace(X'*(B*X*D));
    end

    % Euclidean gradient of the cost function
    function g = egrad(X)
        g = 2*B*X*D;
    end

    problem.precon = preconBFGS(problem);
    problem.linesearch = @(x, xdot, storedb, key) 1;
    options.beta_type ='steep';
    options.tolgradnorm = 1e-5;
    conjugategradient(problem, [], options);

    pause;

    % For comparison, run RTR-FD
    warning('off', 'manopt:getHessian:approx');
    problem = rmfield(problem, 'precon');
    trustregions(problem);

end

author: NicolasBoumal

repo: manopt


text:
%--- help for sym/hessian ---
%
% HESSIAN Hessian matrix.
%    HESSIAN(f,x) computes the Hessian of the scalar f with respect
%    to the vector x. The (i,j)-th entry of the resulting matrix is
%    (d^2f/(dx(i)dx(j)). Note that scalar x is allowed although this
%    is just DIFF(f,x,2).
%
%    Example:
%        syms x y z; hessian(x*y + 2*z*x, [x, y, z])
%        returns  [0, 1, 2; 1, 0, 0; 2, 0, 0]
%
%    See also SYM/CURL, SYM/DIFF, SYM/GRADIENT, SYM/JACOBIAN, SYM/DIVERGENCE,
%    SYM/POTENTIAL, CURL, DIVERGENCE, HESSIAN, JACOBIAN, LAPLACIAN,
%    VECTORPOTENTIAL, SUBS.
%
%    Other functions named hessian
%
%       generic/hessian
%

author: jmaih


text:
function paracellularAbsorption = calculateParacellularAbsorption(DrugLogP, DrugMolecularWeight)
% calculate the paracellular absorption rate using the heuristics presented
% in Eq. 2 in:
% Peters, S. A. (2008). Evaluation of a generic physiologically based pharmacokinetic model for lineshape analysis. Clinical pharmacokinetics, 47(4), 261-75.
% Copyright 2012 The MathWorks, Inc.

paracellularAbsorption = 0;
if DrugLogP > 0.7,
    paracellularAbsorption = 0;
elseif DrugLogP < 0.1 && (DrugMolecularWeight > 200 && DrugMolecularWeight < 360),
    paracellularAbsorption = 0.1;
elseif DrugMolecularWeight < 200,
    paracellularAbsorption = -0.0045*DrugMolecularWeight+1;
end

%paracellularAbsorption = paracellularAbsorption*5;


text:
function quad = fejer2_integrate_fast ( f, n )

%*****************************************************************************80
%
%% FEJER2_INTEGRATE_FAST approximates an integral using a Fejer type 2 rule.
%
%  Discussion:
%
%    The function is integrated over the interval [-1,1].
%
%    The abscissas of the rule of order N are
%
%      x(1:n) = cos ( ( 0:n-1 ) / n )
%
%    The computation should be very efficient in MATLAB.
%
%  Modified:
%
%    08 October 2006
%
%  Author:
%
%    Joerg Waldvogel
%
%  Reference:
%
%    Charles Clenshaw, Alan Curtis,
%    A Method for Numerical Integration on an Automatic Computer,
%    Numerische Mathematik,
%    Volume 2, Number 1, December 1960, pages 197-205.
%
%    Lloyd Trefethen,
%    Is Gauss Quadrature Better than Clenshaw-Curtis?,
%    SIAM Review,
%    Volume 50, Number 1, 2008, pages 67-87.
%
%    Joerg Waldvogel,
%    Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules
%    BIT Numerical Mathematics
%    Volume 43, Number 1, pages 1-18, 2003.
%
%  Parameters:
%
%    Input, function F, an expression, or the name of a function
%    to integrate.
%
%    Input, integer N, the order of the rule to use.
%    N must be at least 1.
%
%    Output, real QUAD, the approximate integral of the function F
%    over the interval [-1,1], using the quadrature rule.
%
  if ( n == 1 )

    x = 0.0;
    w = 2.0;
    quad = w * feval ( f, x );

  else

    N = [ 1 : 2 : n-1 ]';
    L = length ( N );
    m = n - L;
    v0 = [ 2./N./(N-2); 1/N(end); zeros(m,1) ];
    v2 = -v0(1:end-1) - v0(end:-1:2);
    w = ifft ( v2 );

    x = cos ( pi * ( 0:n-1)' / n );
    fx = feval ( f, x ) ;

    quad = w' * fx;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function weekday_test ( )

%*****************************************************************************80
%
%% WEEKDAY_TEST tests the WEEKDAY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WEEKDAY_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the WEEKDAY library.\n' );

  weekday_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WEEKDAY_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function weekday_test01 ( )

%*****************************************************************************80
%
%% WEEKDAY_TEST01 tests YMD_TO_WEEKDAY_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WEEKDAY_TEST01\n' );
  fprintf ( 1, '  For dates in the Common calendar:\n' );
  fprintf ( 1, '  YMD_TO_WEEKDAY_COMMON returns the day of the week.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMD                   Weekday    Weekday\n' );
  fprintf ( 1, '                        Tabulated  Computed\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, y, m, d, w1 ] = weekday_values ( n_data );

    if ( n_data == 0 )
      break;
    end

    s3 = ymd_to_s_common ( y, m, d );
    w2 = ymd_to_weekday_common ( y, m, d );
    s1 = weekday_to_name_common ( w1 );
    s2 = weekday_to_name_common ( w2 );

    fprintf ( 1, '  %20s  %9s  %9s\n', s3, s1, s2 );

  end

  return
end


text:
function AInterior=ellipsInsideConcEllips(A1,A2,gammaVal)
%%ELLIPSINSIDECONCELLIPS Given two concentric ellipsoids (they have the
%               same center), find the largest ellipsoid that fits inside
%               both of them. It will also have the same center.
%
%INPUTS: A1,A2 Two numDimXnumDim symmetric, positive definite matrices that
%           specify the size and shape of the ellipse or ellipsoid, where
%           a point zp is on the ellipse/ellipsoid if
%           (zp-z)'*A*(zp-z)=gammaVal.
%           %Both ellipsoids share the same center point, z.
%  gammaVal The threshold for declaring a point to be in the ellipsoid. If
%           this parameter is omitted or an empty matrix is passed, the
%           default value of 1 is used.
%
%OUTPUTS: AInterior A numDimXnumDim symmetric, positive definite matrix
%            that specifies the size and shape of the largest ellipse that
%            is within both of the ellipses.
%
%This function implements the algorithm that is given in Section 17 of [1].
%
%EXAMPLE:
% A1=[30,0;
%     0,1];
% A2=[1,0;
%     0,21];
% M=[0.413074198133900,  0.910697373904216;
%   -0.910697373904216,   0.413074198133900];%A rotation matrix
% A2=M*A2*M';
% AInterior=ellipsInsideConcEllips(A1,A2);
% figure(1)
% clf
% hold on
% drawEllipse([0;0],A1,1,'b','linewidth',4)
% drawEllipse([0;0],A2,1,'b','linewidth',4)
% drawEllipse([0;0],AInterior,1,'r','linewidth',2)
%
%REFERENCES:
%[1] S. B. Pope, "Algorithms for ellipsoids," Cornell University, Tech.
%    Rep. FDA-08-01, Feb. 2008. [Online].
%    Available: https://tcg.mae.cornell.edu/pubs/Pope_FDA_08.pdf
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(gammaVal))
   gammaVal=1;
end

A1=A1/gammaVal;
A2=A2/gammaVal;

L1=chol(A1,'lower');
L2=chol(A2,'lower');

L2Prime=L1\L2;

[U,Sigma,~]=svd(L2Prime);
SigmaTilde=diag(max(diag(Sigma).^2,1));
APrime=U*SigmaTilde*U';
LPrime=chol(APrime,'lower');

L=L1*LPrime;
AInterior=L*L';

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 [OutImg,OutImgIdx] = PCA_output(InImg, InImgIdx, PatchSize, NumFilters, V)

ImgZ = length(InImg);
mag = (PatchSize-1)/2;
OutImg = cell(NumFilters*ImgZ,1);
cnt = 0;
for i = 1:ImgZ
    [ImgX, ImgY, NumChls] = size(InImg{i});
    img = zeros(ImgX+PatchSize-1,ImgY+PatchSize-1, NumChls);
    img((mag+1):end-mag,(mag+1):end-mag,:) = InImg{i};
    im = im2col_mean_removal(img,[PatchSize PatchSize]); % collect all the patches of the ith image in a matrix, and perform patch mean removal
    for j = 1:NumFilters
        cnt = cnt + 1;
        OutImg{cnt} = reshape(V(:,j)'*im,ImgX,ImgY);  % convolution output
    end
    InImg{i} = [];
end
OutImgIdx = kron(InImgIdx,ones(NumFilters,1));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function matrix = genMatVIFE3DmgHcurl(A,mesh,fem,femI)
%% Generate global matrices and load vector of PPIFEM for 3D ellipic eq
%     -div(A grad u)  = f,    x\in \Omega
%      where A is a piecewise constant on Omega^+ and Omega^-.
% INPUTS:
% pde --- given data function from equation, e.g.
%         pde.A --- diffusion coefficient
%         pde.f --- right hand side function
%         pde.gD --- Dirichlet boundary value function
%         pde.one --- constant function 1.
% mesh --- mesh structure.
% fem --- global degree of freedom of FEM
% femI --- quadrature info on interface cells
% femIF --- quadrature info on interface faces, required in PPIFE.
% PPtype --- partial penalty type possible value
%            'N' : Nonsymmetric PPIFE (e = 1)
%            'S' : Symmetric PPIFE (e = -1)
%            'I' : Incomplete PPIFE (e = 0)
% OUTPUTS:
% matrix.S --- stiffness matrix (w/o boundary condition)
% matrix.E --- consistence matrix (w/o boundary condition)
% matrix.P --- penalty matrix (w/o boundary condition)
% matrix.A --- final FEM matrix (after boundary condition)
% matrix.rhsF --- load vector (w/o boundary condition)
% matrix.rhsE --- consistence vector (=0 if no interface face on boundary)
% matrix.rhsP --- penalty vector (=0 if no interface face on boundary)
% matrix.f --- final RHS matrix (after boundary condition)

% Last Modified: 08/07/2020 by Xu Zhang

%% 1. Stiffness Matrix
S = globMatrixVIFE3DStiff(A,mesh,femI,fem,fem);

Ndof = size(femI.p,1);
bdidx = zeros(Ndof,1);
isBdEdge = true(Ndof,1);
isBdEdge(femI.mapper) = false;
bdidx(isBdEdge) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*S*T + Tbd;

%% 5. Outputs
matrix = struct('A',A);

repo: ifem


text:
function y = testall( func, x )
y = do_test( func, x.value_ );

function y = do_test( func, x )
switch class( x ),
    case 'struct',
        y = do_test( func, struct2cell( x ) );
    case 'cell',
        y = all( cellfun( func, x ) );
    otherwise,
        y = feval( func, x );
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function Y = clfDecTreeFwd( clf, X )
% Apply the decision tree to data X.
%
% USAGE
%  Y = clfDecTreeFwd( clf, X )
%
% INPUTS
%  clf     - trained model
%  X       - nxp data array
%
% OUTPUTS
%  Y       - nx1 vector of labels predicted according to the clf
%
% EXAMPLE
%
% See also CLFDECTREE
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2008 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

if(~strcmp(clf.type,'dectree')); error( ['incor. type: ' clf.type] ); end;
if( size(X,2)~= clf.p ); error( 'Incorrect data dimension' ); end;
T = clf.T;

[Y,d,cnames] = treeval( T, X );
Y = str2double( cnames ); % convert Y back to an int format

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function DEM_HB_and_LE
%--------------------------------------------------------------------------
% This routine is a numerical examination of the relationship between
% entropy, mutual information and the exponential divergence of
% trajectories as the Rayleigh parameter of a Lorenz attractoris increased
% - through a pitchfork bifurcation and subsequent (subcritical) Hopf
% bifurcation. The (stochastic) Lorentz system is integrated for different
% values of the Rayleigh parameter. The nonequilibrium steady-state density
% is then estimated by embedding into a discrete state space; while the
% bifurcations are characterised in terms of the maximal Lyapunov exponent.
% The key thing to observe is the decrease in entropy of blanket states
% prior to the Hopf bifurcation and implicit exponential divergence of
% trajectories. This is scored by the maximal Lyapunov exponent crossing
% zero. Here, the form of the Lorenz attractor defines the three states as
% active, sensory and hidden. Note that there are no internal states in
% this example and blanket states become the particular states (i.e., the
% states of a particle).
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_HB_and_LE.m 7502 2018-12-02 12:28:03Z karl $

% generative model
%==========================================================================                       % switch for demo
spm_figure('GetWin','DEM'); clf

% flow and Jacobian functions
%--------------------------------------------------------------------------
f  = @(x,v,P,G) v(:) + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x/64;
Df = @(x,v,P,G) [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]/64 + ...
                [0 0 0; -x(3) 0 0; 0 x(1) 0]/64;

% model with initial states
%--------------------------------------------------------------------------
G.x  = [1; 1; 24];
G.f  = f;

% set up
%--------------------------------------------------------------------------
b     = 32;                         % number of bins for density estimation
T     = 2^20;                       % length of trajectory
W     = 32;                         % precision of intrinsic fluctuations
P     = exp((-16:32)*log(32)/32);   % Rayleigh parameter range
for k = 1:length(P)

    % integrated timeseries
    %----------------------------------------------------------------------
    U.u   = randn(T + 256,3)/W;
    Pk    = [10; -8/3; P(k)];
    t     = spm_int_L(Pk,G,U);

    % remove intial transients
    %----------------------------------------------------------------------
    t     = t(256:end,:);
    t     = t(1:T,:);
    tt{k} = t;

    % sample density
    %----------------------------------------------------------------------
    for i = 1:3
        t(:,i) = t(:,i) - min(t(:,i));
        t(:,i) = (b - 2)*t(:,i)/max(t(:,i));
    end
    t     = t + 1;
    p     = zeros(b,b,b) + 1;
    S     = 0;
    for i = 1:T

        % accumulate in bins
        %------------------------------------------------------------------
        j = round(t(i,:));
        p(j(3),j(2),j(1)) = p(j(3),j(2),j(1)) + 1;

        % Lyapunov exponents
        %------------------------------------------------------------------
        dfdx = Df(tt{k}(i,:)',zeros(1,3)',Pk);
        S    = S + sort(real(eig(dfdx,'nobalance')),'descend');

    end
    p     = p/sum(p(:));
    pp{k} = p;

    % mutual informations of sample density
    %----------------------------------------------------------------------
    [I,Ii,Ie] = spm_self_entropy(p);
    MI(1,k)   = I;
    MI(2,k)   = Ii;
    MI(3,k)   = Ie;
    LE(:,k)   = S/T;

    % plot mutual informations
    %----------------------------------------------------------------------
    subplot(3,2,3)
    semilogx(P(1:k),MI(1,:),'b',P(1:k),MI(2,:),'b-.',P(1:k),MI(3,:),'b:')
    axis square xy
    title('Expected surpise','Fontsize',16)
    xlabel('Control parameter'),ylabel('Entropies (nats)'),drawnow

    % plot maximal Lyapunov exponent
    %----------------------------------------------------------------------
    subplot(3,2,4)
    semilogx(P(1:k),LE(1,:),'b')
    axis square xy
    title('Lyapunov exponent','Fontsize',16)
    xlabel('Control parameter'),ylabel('Principal exponent'),drawnow

end

% lines and thresholds
%--------------------------------------------------------------------------
j     = find(LE(1,:) > 0,1);
j     = P(j);
subplot(3,2,3)
hold on, plot([j j],[ 0 6],':'), hold off
hold on, plot([1 1],[ 0 6],':'), hold off
legend({'H(B)','H(B|E)','I(B,E)'})
subplot(3,2,4)
hold on, plot([j j],[-1 1]*0.03,':'), hold off
hold on, plot([1 1],[-1 1]*0.03,':'), hold off
hold on, plot(P,zeros(size(P)),'--'), hold off

% illustrate trajectories and ergodic density
%--------------------------------------------------------------------------
j     = [8 38 length(P)];
for i = 1:length(j)

    % exemplar trajectory (plot)
    %----------------------------------------------------------------------
    subplot(3,3,i)
    t    = tt{j(i)};
    plot(t(1:1024,2),t(1:1024,3),'k')
    axis square xy
    title('Trajectory','Fontsize',16)
    axis([-30 30 -5 60])

    % image format
    %----------------------------------------------------------------------
    subplot(3,3,6 + i)
    p = pp{j(i)};
    imagesc(1-squeeze(sum(p,2))'),axis xy square
    title('Marginal density','Fontsize',16)
    ylabel('State'), xlabel('State')

%     subplot(3,3,6 + i)
%     [W,S,X,beta] = spm_power_law(t');
%     plot(W,S,'b.',W,X*beta,'b','LineWidth',1)
%     title(sprintf('alpha = %-2.2f',beta(2)),'FontSize',16)
%     ylabel('Log power'), xlabel('Log frequency')
%     axis square, axis xy, spm_axis tight

end

return

function [HB,HBH,IBH] = spm_self_entropy(pHxB)
% FORMAT [HB,HBH,IBH] = spm_self_entropy(pHxB)
% Entropies
% HS  = H(B)              % self entropy
% HBH = H(B|H)            % conditional entropy
% IBH = I(B,H)            % mutual information
%
% This subroutine assumes that the first dimension of the joint density
% corresponds to a hidden or external state and the rest are particular
% or blanket states

% evaluate joint density and posterior
%--------------------------------------------------------------------------
pB    = sum(pHxB,1);
pH    = sum(sum(sum(pHxB,2),3),4);

% inline functions
%--------------------------------------------------------------------------
ln    = @(p)log(spm_vec(p) + 1e-16);
H     = @(p)-spm_vec(p)'*ln(p);

% relative entropies
%--------------------------------------------------------------------------
HB    = H(pB);
IBH   = H(pH) + H(pB) - H(pHxB);
HBH   = HB - IBH;

return

function D = spm_Kaplan_Yorke(LE)
% FORMAT Kaplan Yorke estimate of dimensional complexity
%--------------------------------------------------------------------------
L     = real(LE);
for i = 1:size(L,2)
    l = sort(L(:,i),'descend');
    j = find(cumsum(l) > 0,1);
    if isempty(j), j = 0; end
    D(i) = j + sum(l(1:j))/abs(l(j + 1));
end

function [W,S,X,beta] = spm_power_law(x)
% FORMAT spm_power_law(x)

% illustrate power law scaling
%--------------------------------------------------------------------------
N     = floor(log2(size(x,2)));
s     = abs(fft(x(1,:)')).^2;
w     = (1:2^12)';
W     = w;
S     = s(w + 1);

S     = decimate(log(S),N - 4);
W     = log(decimate(W,N - 4));
X     = [ones(size(W)),W];

% plot part of trajectory
%--------------------------------------------------------------------------
[~,i] = max(abs(diff(spm_conv(x(1,:),2^(N - 8)))));
nn    = 2^10;
i     = (-nn:nn) + i;
i     = i(i > 0 & i < size(x,2));

% estimate exponent (alpha)
%--------------------------------------------------------------------------
[~,~,beta] = spm_ancova(X,[],S,[0;1]);

author: spm


text:
function [sol,info] = solvebilevel(varargin)
%SOLVEBILEVEL Simple global bilevel solver
%
%   min        Objective(Pouter)(x,y)
%   subject to Constraints(Pouter)(x,y)>0
%              y = arg min Objective(Pinner)(x,y)
%              subject to Constraints(Pinner)(x,y)>0
%
%   [DIAGNOSTIC,INFO] = SOLVEBILEVEL(Pouter, Pinner, y, options)
%
%   diagnostic : Struct with standard YALMIP diagnostics
%   info       : Bilevel solver specific information
%
%   Input
%      Pouter   : Outer problem
%      Pinner   : Inner problem
%      y        : Inner variables
%      options  : solver options from SDPSETTINGS.
%
%   The behaviour of the bilevel solver can be controlled
%   using the field 'bilevel' in SDPSETTINGS
%
%      bilevel.outersolver : Solver for outer problems with inner KKT removed
%      bilevel.innersolver : Solver for inner problem
%      bilevel.rootcut     : Number of cuts (based on complementary
%                            constraints) added in root (experimental)
%      bilevel.relgaptol   : Termination tolerance
%      bilevel.compslacktol: Tolerance for accepting complementary slackness
%      bilevel.feastol     : Tolerance for feasibility in outer problem
%
%
%   See also SDPVAR, SDPSETTINGS, SOLVESDP

if isa(varargin{2},'optimizer')
    Pout = varargin{1};
    s = struct(varargin{2});
    if nargin < 3
        options = sdpsettings;
    else
        options = varargin{3};
    end
    [sol,info] = solvebilevel(Pout.Constraints,Pout.Objective,s.F,s.h,s.output.expression,options)
else
    Pout = varargin{1};
    Pinn = varargin{2};
    y = varargin{3};

    if nargin < 4
        options = sdpsettings;
    else
        options = varargin{4};
    end

    [sol,info] = solvebilevel(Pout.Constraints,Pout.Objective,Pinn.Constraints,Pinn.Objective,y,options)
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   q = solve_spherical_wrist_flexion_n2(robot, q, T, wrist)
%   Solves the inverse kinematic problem for a spherical wrist of the robot
%   Epson Flexion N2
%   robot: robot structure.
%   q: vector containing the values of the joints 1, 2 and 3.
%   T: orientation of the last reference system.
%   wrist: select -1 or 1 for two possible solutions (wrist up, wrist down)
%
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

function q = solve_spherical_wrist_flexion_n2(robot, q, T, wrist, method)

switch method

     %algebraic solution
    case 'algebraic'
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);

        Q=inv(T23)*inv(T12)*inv(T01)*T;

        %detect the degenerate case when q(5)=0, this leads to zeros
        % in Q13, Q23, Q31 and Q32 and Q33=1
        thresh=1e-12;
        %detect if q(5)==0
        % this happens when cos(q5) in the matrix Q is close to 1
        if abs(Q(3,3)-1)>thresh
            %normal solution
            if wrist==1 %wrist up
                q(4)=atan2(-Q(2,3),-Q(1,3))+pi;
                q(6)=atan2(-Q(3,2),Q(3,1));
                %q(5)=atan2(-Q(3,2)/sin(q(6)),Q(3,3));
            else %wrist down
                q(4)=atan2(-Q(2,3),-Q(1,3))-pi+pi;
                q(6)=atan2(-Q(3,2),Q(3,1))+pi;
                %q(5)=atan2(-Q(3,2)/sin(q(6)),Q(3,3));
            end
            if abs(cos(q(6)+q(4)))>thresh
                cq5=(Q(1,1)+Q(2,2))/cos(q(4)+q(6))-1;
            end
            if abs(sin(q(6)+q(4)))>thresh
                cq5=(-Q(1,2)+Q(2,1))/sin(q(4)+q(6))-1;
            end
            if abs(sin(q(6)))>thresh
                sq5=-Q(3,2)/sin(q(6));
            end
            if abs(cos(q(6)))>thresh
                sq5=Q(3,1)/cos(q(6));
            end
            q(5)=atan2(sq5,cq5)+pi;

        else %degenerate solution, in this case, q4 cannot be determined,
             % so q(4)=0 is assigned
            if wrist==1 %wrist up
                q(4)=0;
                q(5)=0;
                q(6)=atan2(-Q(1,2)+Q(2,1),Q(1,1)+Q(2,2))-pi/2;
            else %wrist down
                q(4)=-pi;
                q(5)=0;
                q(6)=atan2(-Q(1,2)+Q(2,1),Q(1,1)+Q(2,2))+pi-pi/2;
            end

        end

       %geometric solution
    case 'geometric'
        % T is the noa matrix defining the position/orientation of the end
        % effector's reference system
        vx6=T(1:3,1);
        vz5=T(1:3,3); % The vector a z6=T(1:3,3) is coincident with z5

        % Obtain the position and orientation of the system 3
        % using the already computed joints q1, q2 and q3
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);
        T03=T01*T12*T23;

        vx3=T03(1:3,1);
        vy3=T03(1:3,2);
        vz3=T03(1:3,3);

        % find z4 normal to the plane formed by z3 and a
        vz4=cross(vz3, vz5);	% end effector's vector a: T(1:3,3)

        % in case of degenerate solution,
        % when vz3 and vz6 are parallel--> then z4=0 0 0, choose q(4)=0 as solution
        if norm(vz4) <= 0.000001
            if wrist == 1 %wrist up
                q(4)=0;
            else
                q(4)=-pi; %wrist down
            end
        else
            %this is the normal and most frequent solution
            cosq4=wrist*dot(-vy3,vz4);
            sinq4=wrist*dot(vx3,vz4);
            q(4)=atan2(sinq4, cosq4)+pi;
        end
        %propagate the value of q(4) to compute the system 4
        T34=dh(robot, q, 4);
        T04=T03*T34;
        vx4=T04(1:3,1);
        vy4=T04(1:3,2);

        % solve for q5
        cosq5=dot(vy4,vz5);
        sinq5=dot(-vx4,vz5);
        q(5)=atan2(sinq5, cosq5)+pi;

        %propagate now q(5) to compute T05
        T45=dh(robot, q, 5);
        T05=T04*T45;
        vx5=T05(1:3,1);
        vy5=T05(1:3,2);

        % solve for q6
        cosq6=dot(vx6,vx5);
        sinq6=dot(vx6,vy5);
        q(6)=atan2(sinq6, cosq6)-pi/2;

    otherwise
        disp('no method specified in solve_spherical_wrist');
end

repo: ARTE


text:
function b = r8gd_to_r8ge ( n, ndiag, offset, a )

%*****************************************************************************80
%
%% R8GD_TO_R8GE copies a R8GD matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R8GD storage format is suitable for matrices whose only nonzero entries
%    occur along a few diagonals, but for which these diagonals are not all
%    close enough to the main diagonal for band storage to be efficient.
%
%    In that case, we assign the main diagonal the offset value 0.
%    Each successive superdiagonal gets an offset value 1 higher, until
%    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
%    Similarly, the subdiagonals are assigned offsets of -1 through -(N-1).
%
%    Now, assuming that only a few of these diagonals contain nonzeros,
%    then for the I-th diagonal to be saved, we stored its offset in
%    OFFSET(I), and its entries in column I of the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer NDIAG, the number of diagonals of the matrix
%    that are stored in the array.
%    NDIAG must be at least 1, and no more than 2 * N - 1.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal storage.
%
%    Input, real A(N,NDIAG), the R8GD matrix.
%
%    Output, real B(N,N), the R8GE matrix.
%
  b(1:n,1:n) = 0.0E+00;

  for i = 1 : n
    for diag = 1 : ndiag
      j = i + offset(diag);
      if ( 1 <= j & j <= n )
        b(i,j) = a(i,diag);
      end
    end
  end

  return
end


text:
function calpak_test25 ( )

%*****************************************************************************80
%
%% CALPAK_TEST25 tests JED_TO_YMDF_REPUBLICAN and YMDF_TO_JED_REPUBLICAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST25\n' );
  fprintf ( 1, '  For the Republican calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_REPUBLICAN: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_REPUBLICAN: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF                JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_republican ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_republican ( jed1 );

      s2 = ymdf_to_s_republican ( y2, m2, d2, f2 );

      jed3 = ymdf_to_jed_republican ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
% This m-file produces the scales used in realized_range.  It can be run with more simulations to
% increase the precision of the estimated scales used.  If you do this you will have to copy and
% paste the results into the function since it uses a lookup table and interpolation.

% Clear
clear all
% Number of simulations
BB = 1000000;

% Number of prices to use in each interval.  All integer factors of 23400.  The asymptotic value is
% 4*log(2)

m1 = factor(23400);
m1 = [1 m1];
ms = m1';
for i=2:length(m1);
    temp = nchoosek(1:length(m1),i);
    ms = [ms;    unique(prod(m1(temp),2))]; %#ok<AGROW>
end
ms = unique(ms);
maxM = max(ms);

% Turn the number of prices per window in to a step size
gap = maxM./ms;
ms = ms+1;
% Initialize a place to hold results
M = length(ms);
MC = zeros(BB,M);
% BB simultions
tic
for j = 1:BB
    % Resample a single BM 1, 2, ..., 23400 times.
    x = [0;cumsum(randn(maxM,1)/sqrt(maxM))];
    for i=1:M
        x2 = x(1:gap(i):(maxM+1));
        MC(j,i)  = (max(x2)-min(x2))^2;
    end
    % Display the count and the expected time remaining.
    if mod(j,10000)==0
        t=toc;
        str = [num2str(j) ' iterations complete.  The elapsed time is ' num2str(t) '. The expected time remaining is ' num2str(BB*t/j) ];
        disp(str)
    end
end

% Some smoothing using a concave regression
y = mean(MC);
y(1) = 1;
orig_mean_MC=y;
n = length(y);
x = eye(n);
A1 = [eye(n-1) zeros(n-1,1)] + [zeros(n-1,1) -eye(n-1)];
A2 = [eye(n-2) zeros(n-2,2)] + [zeros(n-2,1) -2*eye(n-2) zeros(n-2,1)] + [zeros(n-2,2) eye(n-2)];
b = zeros(2*n-3,1);
concave_mean_MC=lsqlin(x,y,[A1;A2],b);
% Save the results.
%save realized_range_simulation_results orig_mean_MC concave_mean_MC ms

author: bashtage

repo: mfe-toolbox


text:
function [varargout] = buffer(varargin)

% BUFFER manages and accesses the realtime data acquisition buffer
% This function is implented as a mex file.
%
% Use as
%   retval = buffer(cmd, detail, host, port)
%
% To read data from a buffer server over the network
%   hdr = buffer('get_hdr', [],     host, port)
%   dat = buffer('get_dat', datsel, host, port)
%   evt = buffer('get_evt', evtsel, host, port)
%
% The selection for data and events should be zero-offset and contain
%   datsel = [begsample endsample]
%   evtsel = [begevent  endevent ]
%
% To write data to a buffer server over the network
%   buffer('put_hdr', hdr, host, port)
%   buffer('put_dat', dat, host, port)
%   buffer('put_evt', evt, host, port)
%
% To implement a local buffer server and have other clients
% connect to it at a specified network port
%   buffer('tcpserver', 'init', [], port)
%   buffer('tcpserver', 'exit', [], port)

% Copyright (C) 2008-2010, 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$

error('Could not locate mex file');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Mean-Shift Video Tracking
% by Sylvain Bernhardt
% July 2008
%% Description
% This is a simple example of how to use
% the Mean-Shift video tracking algorithm
% implemented in 'MeanShift_Algorithm.m'.
% It imports the video 'Ball.avi' from
% the 'Videos' folder and tracks a selected
% feature in it.
% The resulting video sequence is played after
% tracking, but is also exported as a AVI file
% 'Movie_out.avi' in the 'Videos' folder.

clear all
close all

%% Import movie
% and time it with tic/toc
tic
[Length,height,width,Movie]=...
    Import_mov('Videos/Ball.avi');
toc

%% Play the movie
% % Put the figure in the center of the screen,
% % without axes and menu bar.
% scrsz = get(0,'ScreenSize');
% figure('Position',[scrsz(3)/2-width/2 scrsz(4)/2-height/2 width height],...
%     'MenuBar','none');
% axis off
% % Image position inside the figure
% set(gca,'Units','pixels','Position',[1 1 width height])
% % Play the movie
% movie(Movie);

%% Variables
index_start = 1;
% Similarity Threshold
f_thresh = 0.16;
% Number max of iterations to converge
max_it = 5;
% Parzen window parameters
kernel_type = 'Gaussian';
radius = 1;

%% Target Selection in Reference Frame
[T,x0,y0,H,W] = Select_patch(Movie(index_start).cdata,0);
pause(0.2);

%% Run the Mean-Shift algorithm
% Calculation of the Parzen Kernel window
[k,gx,gy] = Parzen_window(H,W,radius,kernel_type,0);
% Conversion from RGB to Indexed colours
% to compute the colour probability functions (PDFs)
[I,map] = rgb2ind(Movie(index_start).cdata,65536);
Lmap = length(map)+1;
T = rgb2ind(T,map);
% Estimation of the target PDF
q = Density_estim(T,Lmap,k,H,W,0);
% Flag for target loss
loss = 0;
% Similarity evolution along tracking
f = zeros(1,(Length-1)*max_it);
% Sum of iterations along tracking and index of f
f_indx = 1;
% Draw the selected target in the first frame
Movie(index_start).cdata = Draw_target(x0,y0,W,H,...
    Movie(index_start).cdata,2);
%%%% TRACKING
WaitBar = waitbar(0,'Tracking in progress, be patient...');
% From 1st frame to last one
for t=1:Length-1
    % Next frame
    I2 = rgb2ind(Movie(t+1).cdata,map);
    % Apply the Mean-Shift algorithm to move (x,y)
    % to the target location in the next frame.
    [x,y,loss,f,f_indx] = MeanShift_Tracking(q,I2,Lmap,...
        height,width,f_thresh,max_it,x0,y0,H,W,k,gx,...
        gy,f,f_indx,loss);
    % Check for target loss. If true, end the tracking
    if loss == 1
        break;
    else
        % Drawing the target location in the next frame
        Movie(t+1).cdata = Draw_target(x,y,W,H,Movie(t+1).cdata,2);
        % Next frame becomes current frame
        y0 = y;
        x0 = x;
        % Updating the waitbar
        waitbar(t/(Length-1));
    end
end
close(WaitBar);
%%%% End of TRACKING

%% Export/Show the processed movie
% Export the video sequence as an AVI file in the Videos folder
% WaitBar = waitbar(0,'Exporting the output AVI file, be patient...');
% movie2avi(Movie,'Videos\Movie_out','Quality',50);
% waitbar(1);
% close(WaitBar);

% Put a figure in the center of the screen,
% without menu bar and axes.
scrsz = get(0,'ScreenSize');
figure(1)
set(1,'Name','Movie Player','Position',...
    [scrsz(3)/2-width/2 scrsz(4)/2-height/2 width height],...
    'MenuBar','none');
axis off
% Image position inside the figure
set(gca,'Units','pixels','Position',[1 1 width height])
% Play the movie
movie(Movie);

%% End of File
%=============%


text:
function[j]=mat2vect(i)
j = [];
j =[i(1,:).';i(2,:).';i(3,:).';i(4,:).';i(5,:).';i(6,:).';i(7,:).'; ...
    i(8,:).';i(9,:).';i(10,:).';i(11,:).';i(12,:).';i(13,:).' ...
    ;i(14,:).';i(15,:).';i(16,:).';i(17,:).';i(18,:).' ...
    ;i(19,:).';i(20,:).';i(21,:).';i(22,:).';i(23,:).';i(24,:).';i(25,:).'...
    ;i(26,:).';i(27,:).';i(28,:).';i(29,:).';i(30,:).';i(31,:).';i(32,:).'...
    ;i(33,:).';i(34,:).';i(35,:).'];
end


text:
function demoDrawTubularMesh
%DEMODRAWTUBULARMESH Draw a tubular mesh using matGeom library
%
%   Example
%     demoDrawTubularMesh
%
%   See also
%     drawSoccerBall, demoRevolutionSurface, demoEquivalentEllipsoid
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2012-04-25,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

%% Initialisation

% extract vertices, edges, and faces of soccerball polyhedron
[vertices, edges, faces] = createSoccerBall;

% prepare figure
figure(1); clf; hold on;
axis equal;
view(3);

% draw the polyhedron as basis
drawPolyhedron(vertices, faces);

%% Draw each edge as a green cylinder

figure; hold on;
axis equal;
view(3);

width = .05;
radius = .1;

for i = 1:size(edges, 1)
    p1 = vertices(edges(i, 1), :);
    p2 = vertices(edges(i, 2), :);

    drawCylinder([p1 p2 width], 'FaceColor', 'g');
end

%% Draw each vertex as a blue ball

for i = 1:size(vertices, 1)
    drawSphere([vertices(i, :) radius], 'FaceColor', 'b');
end

%% More setup on display

set(gcf, 'renderer', 'opengl')
light;
view([10 30]);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%% maxwell_run
% Run MaxwellFDFD.

%%% Syntax
%  [E_cell, H_cell] = maxwell_run(OSC, DOM, OBJ, SRC, [inspect_only])
%  [E_cell, H_cell, obj_array, src_array] = maxwell_run(...)
%  [E_cell, H_cell, obj_array, src_array, extra] = maxwell_run(...)

%%% Description
% |maxwell_run(OSC, DOM, OBJ, SRC)| constructs a simulation domain from given
% objects and sources, and launches a finite-difference frequency-domain (FDFD)
% solver to solve Maxwell's equations in the simulation domain.
%
% Each of |OSC|, |DOM|, |OBJ|, and |SRC| represents a group of parameters.
% Each group supports several flexible expressions.  See the following sections
% about the input parameter groups for more details.
%
% |maxwell_run| can take an optional argument |inspect_only|, which is a logical
% argument (i.e., |true| or |false|).  If |inspect_only = true|, |maxwell_run|
% runs without launching the solver.  This is useful to inspect input arguments
% before starting expensive computation.
%
% |[E_cell, H_cell] = maxwell_run(...)| returns the E- and H-field solutions of
% Maxwell's equations.  |E_cell| and |H_cell| are length-3 row cell arrays whose
% elements are the x-, y-, and z-components of the respective field solutions.
% Each x-, y-, z-component of the fields are provided as instances of
% <Scalar3d.html Scalar3d>.
%
% |[E_cell, H_cell, obj_array, src_array] = maxwell_run(...)| returns arrays of
% instances of <EMObject.html |EMObject|> and <Source.html |Source|>.  The |EMObject|
% and |Source| elements represent the objects and sources placed in the
% simulation domain, so they can be used to visualize the simulation domain.
%
% |[E_cell, H_cell, obj_array, src_array, extra] = maxwell_run(...)| returns
% |extra|, which is a structure containing the fields |grid3d|, |J|, |M|, |eps|,
% and |mu|: |grid3d| is an instance of <Grid3d.html Grid3d> that contains the
% grid information; |J| and |M| are the electric and magnetic current source
% distributions that have the same structure as |E_cell| and |H_cell|; |eps| and
% |mu| are instances of <Scalar3d.html Scalar3d> that contain the electric
% permittivity and magnetic permeability defined at the centers of the grid
% cells.

%%% Input Parameter Group - OSC
% OSC group specifies oscillation parameters.  The group begins with |'OSC'| and
% ends with one of the followings:
%
% * |L0, wvlen|
% * |osc|
%
% |L0|: unit of length in meters. For example, |L0 = 1e-9| means 1 nm.
% All lengths in other parameters are described in this unit.
%
% |wvlen|: vacuum wavelength in the unit of |L0|.  In other words, the actual
% wavelength is |wvlen * L0|.
%
% |osc|: instance of <Oscillation |Oscillation|>, which contains the information
% about the length unit and vacuum wavelength.

%%% Input Parameter Group - DOM
% DOM group specifies domain parameters.  The group begins with |'DOM'| and ends
% with one of the followings:
%
% * |material, box, dl, bc, Lpml, [deg_pml, [R_pml]], [withuniformgrid]|
% * |domain, dl, bc, Lpml, [deg_pml, [R_pml]], [withuniformgrid]|
%
% |material|: background material filling the simulation domain.  See
% <maxwell_run.html#7 Material Description> for various ways to describe a
% material.
%
% |box|: range of the simulation doamin in the form of |[xmin xmax; ymin ymax;
% zmin zmax]|.
%
% |dl|: default grid cell size in the form of |[dx dy dz]|.  Can be abbreviated
% to a scalar |dl| if |dx == dy == dz|.
%
% |bc|: boundary conditions in the form of |[bc_xn bc_yn bc_zn]|, whose each
% element specifies the boundary condition at the negative end in one of the x-,
% y-, z-axes (e.g., |bc_xn| for the negative end in the x-axis.).  Each element
% of |bc| is an instance of <BC.html |BC|>.  The boundary conditions at the
% positive ends are automatically determined by those at the negative ends:
% |bc_wp == BC.p| if |bc_wn == BC.p|, and |bc_wp == BC.m| otherwise.  Can be
% further abbreviated to a single BC instance if |bc_xn == bc_yn == bc_zn|.
%
% |Lpml|: thicknesses of PML in the form of |[Lpml_xn Lpml_xp; Lpml_yn Lpml_yp;
% Lpml_zn Lpml_zp]|, whose each element specifies the thickness at either
% negative or positive end in one of the x-, y-, z-axes (e.g., |Lpml_xn| for the
% negative end in the x-axis).  Can be abbreviated to |[Lpml_x Lpml_y Lpml_z]|
% if |Lpml_wn == Lpml_wp| for |w = x, y, z|.  Can be further abbreviated to a
% single scalar thickness if |Lpml_x == Lpml_y == Lpml_z|.  All thicknesses are
% in the unit of |L0|.
%
% |domain|: instance of <Domain.html |Domain|>, which contains the information
% about |material| and |box|.
%
% |deg_pml|: optional argument to specify the degrees of the polynomial gradings
% of PML loss parameters in the form of |[deg_xn deg_xp; deg_yn deg_yp; deg_zn
% deg_zp]|.  Can be abbreviated to |[deg_x deg_y deg_z]| if |deg_wn == deg_wp|
% for |w = x, y, z|.  Can be further abbreviated to a single scalar if |deg_x ==
% deg_y == deg_z|.  The default value is |deg_pml = 4|, meaning that the
% polynomials of degree 4 are used to smoothly increase the PML loss parameters
% for all boundaries with PML.
%
% |R_pml|: optional argument to specify the target reflectances in the form of
% |[R_xn R_xp; R_yn R_yp; R_zn R_zp]|.  To specify |R_pml|, |deg_pml| should be
% specified.  Can be abbreviated to |[R_x R_y R_z]| if |R_wn == R_wp| for |w =
% x, y, z|.  Can be further abbreviated to a single scalar if |R_x == R_y ==
% R_z|.  The default value is |R_pml = exp(-16)|, meaning that the reflectance
% from PML is aimed to be as small as |exp(-16) ~= 1e-7|.
%
% |withuniformgrid|: optional argument to select a grid generation algorithm.
% If |withuniformgrid = true|, a grid is generated uniformly; otherwise a grid
% is generated dynamically to best fit the objects in the simulation domain.
% The default value is |withuniformgrid = false|, i.e., dynamic grid generation
% is preferred to uniform grid generation.

%%% Input Parameter Group - OBJ and SOBJ
% OBJ group specifies obect parameters.  The group begins with |'OBJ'| and ends
% with one of the followings:
%
% * |material_1, shapes_1, ..., material_N, shapes_N|
% * |obj_1, ..., obj_N|
% * |eps_node_array|
%
% |material_i|: material filling |shapes_i|.  See <maxwell_run.html#7 Material
% Description> for various ways to describe a material.
%
% |shapes_i|: shapes made of |material_i|.  It can be a single shape, an array
% of shapes (i.e., |[shape_a, shape_b, ...]|), or can be spanned into several
% shape arguments (i.e., |shape_1, shape_2, ...|).  Each shape is an instance of
% <Shape.html |Shape|>.
%
% |obj_i|: instance or array of instances of <EMObject.html |EMObject|>.  Each
% |EMObject| is composed of a material and shape.
%
% |eps_node_array| is a 3D array of permittivities defined at the centers of
% grid cells.  The size of the array should be the same as the number of grid
% cells in a generated grid.  Because the number of grid cells generated by
% dynamic grid generation is hard to predict, uniform grid generation is
% recommended (see DOM group) when |eps_node_array| is used.
%
% There is a similar, optional parameter group SOBJ that begins with |'SOBJ'|.
% SOBJ stands for _scattering objects_, and it is used to define scatterers for
% total-field/scattered-field (TF/SF) simulation.  When a TF/SF source is used
% as a source, the objects defined in SOBJ group are treated as scatterers
% whereas the objects defined in OBJ group are treated as background objects.
% If a TF/SF source is not used, SOBJ group works the same as OBJ group.

%%% Input Parameter Group - SRC
% SRC group specifies current source parameters. Depending on whether the
% sources are electric or magnetic, the group begins with |'SRCJ'| or |'SRCM'|
% and ends with
%
% * |src_1, ..., src_M|
%
% |src_i|: instance or array of instances of <Source.html |Source|>.

%%% Material Description
% Each material is described by one of the followings:
%
% * |{name, color, permittivity}|
% * |{material_datapath, color}|
% * |material|
%
% |name|: string describing the name of the material (e.g., 'vacuum', 'silver',
% 'Ag').
%
% |color|: color to be used in visualizing objects made of the material.  It can
% be either a standard color specifier character used in MATLAB (e.g., |'w'| for
% white, |'k'| for black, etc.) or |[r, g, b]| where |r|, |g|, |b| are RGB color
% scales ranging from |0.0| to |1.0| (e.g., |[0.5, 0.5, 0.5]| for gray).  Use
% |color = 'none'| to not draw the material.
%
% |permittivity|: electric permittivity value in the unit of the vacuum
% permittivity.  It is a complex number.
%
% |material_datapath|: path to the file of the permittivity data of the
% material.  It is a path relative to |material| directory.  For example,
% |'Palik/Si'| refers to silicon's permittivity data stored in |Si.m| file in
% |material/Palik| directory.  The permittivity for the frequency of interest is
% automatically retrieved from the data file.

%%% Example
%   gray = [0.5 0.5 0.5];  % [r g b]
%   inspect_only = true;
%   [E, H, obj_array, extra] = maxwell_run(...
%       'OSC', 1e-9, 1550, ...
%       'DOM', {['Palik/SiO2'], 'none'}, [-700, 700; -600, 600; -200, 1700], 20, BC.p, 200, ...
%       'OBJ', ...
%           {['Palik/SiO2'], 'none'}, Box([-50, 50; -50, 50; -200, 1700], [2, 2, 20]), ...  % OBJ1
%           {['CRC/Ag'], gray}, [Box([-700, -25; -25, 25; -200, 1700], 20), Box([25, 700; -25, 25; -200, 1700], 20)], ...  % OBJ2
%       'SRCJ', PointSrc(Axis.x, [0, 0, 200]), ...
%       inspect_only);

function [E_cell, H_cell, obj_array, src_array, extra] = maxwell_run(varargin)
	DEFAULT_METHOD = 'direct';  % 'direct', 'gpu', 'aws', 'inputfile'

	% Set solver options.
	iarg = nargin; arg = varargin{iarg};
	inspect_only = false;
	if istypesizeof(arg, 'logical')
		inspect_only = arg;
		iarg = iarg - 1; arg = varargin{iarg};
	end

	is_solveropts = false;
	if istypesizeof(arg, 'struct')
		solveropts = arg;
		is_solveropts = true;
		iarg = iarg - 1;  % arg = varargin{iarg};
	end

	if ~is_solveropts || ~isfield(solveropts, 'showstruct')
		solveropts.showstruct = true;
	end

	if ~is_solveropts || ~isfield(solveropts, 'method')
		solveropts.method = DEFAULT_METHOD;
	end

	if is_solveropts && isequal(solveropts.method, 'inputfile')
		chkarg(isfield(solveropts, 'filenamebase'), '"solveropts" should have "filenamebase" field.');
	end

	if ~is_solveropts || ~isfield(solveropts, 'maxit')
% 		solveropts.maxit = intmax;
		solveropts.maxit = 1e6;
	else
		chkarg(istypesizeof(solveropts.maxit, 'real') && solveropts.maxit > 0, ...
			'solveropts.maxit should be positive.');
	end

	if ~is_solveropts || ~isfield(solveropts, 'tol')
		solveropts.tol = 1e-6;
	else
		chkarg(istypesizeof(solveropts.tol, 'real') && solveropts.tol > 0, ...
			'solveropts.tol should be positive.');
	end

	if ~is_solveropts || ~isfield(solveropts, 'eqtype')
		solveropts.eqtype = EquationType(FT.e, GT.prim);
	else
		chkarg(istypesizeof(solveropts.eqtype, 'EquationType'), ...
			'solveropts.eqtype should be instance of EquationType.');
	end

	if ~is_solveropts || ~isfield(solveropts, 'pml')
		solveropts.pml = PML.sc;
	else
		chkarg(istypesizeof(solveropts.pml, 'PML'), ...
			'solveropts.pml should be instance of PML.');
	end

	if ~is_solveropts || ~isfield(solveropts, 'returnAandb')
		solveropts.returnAandb = false;
	else
		chkarg(istypesizeof(solveropts.returnAandb, 'logical'), ...
			'solveropts.returnAandb should be logical.');
	end

	if ~is_solveropts || ~isfield(solveropts, 'returnDiv')
		solveropts.returnDiv = false;
	else
		chkarg(istypesizeof(solveropts.returnDiv, 'logical'), ...
			'solveropts.returnDiv should be logical.');
	end

	chkarg(iarg > 0, 'first argument is not correct.');

	if inspect_only
		fprintf('%s begins (inspection only).\n', mfilename);
	else
		fprintf('%s begins.\n', mfilename);
	end
	ge = solveropts.eqtype.ge;
	fprintf('E-field grid type: %s\n', char(ge));
	pm = ProgMark();

	% Build the system.
	% Make sure to pass the first consecutive elements of varargin to
	% build_system() for correct error reports.
	[osc, grid3d, s_factor, eps, mu, J, M, obj_array, src_array, mat_array, eps_node, mu_node, isiso] = ...
		build_system(ge, solveropts.pml, varargin{1:iarg}, pm);

	if ~is_solveropts || ~isfield(solveropts, 'F0')
		solveropts.F0 = 'zero';  % 'rand' is the other choice
	else
		chkarg(isequal(solveropts.F0, 'zero') || isequal(solveropts.F0, 'rand') ...
			|| istypesizeof(solveropts.F0, 'complexcell', [1 Axis.count], grid3d.N), ...
			'solveropts.F0 should be length-%d cell array whose each element is %d-by-%d-by-%d array with complex numbers.', ...
			Axis.count, grid3d.Ncell{:});
	end

	if inspect_only  % inspect objects and sources
		if solveropts.showstruct
			figure;
			set(gcf, 'units','normalized','position',[0.5 0 0.5 0.5]);
			withpml = true;
			visobjsrc(grid3d, obj_array, src_array, withpml);
			drawnow
			pm.mark('domain visualization');
		end

		% Visualize modes.
		is_modalsrc = false;
		for src = src_array
			if istypesizeof(src, 'ModalSrc')
				is_modalsrc = true;
				modalsrc = src;
				opts.withabs = true;

				E2d = modalsrc.E2d;
				cmax = max(abs([E2d{Axis.x}.array(:); E2d{Axis.y}.array(:); E2d{Axis.z}.array(:)]));
				opts.cmax = cmax;
				for w = Axis.elems
					figure;
					set(gcf, 'units','normalized','position',[(int(w)-1)/3 1/2 1/3 1/3]);
					vis2d(E2d{w}, obj_array, opts);
					drawnow;
				end

				H2d = modalsrc.H2d;
				cmax = max(abs([H2d{Axis.x}.array(:); H2d{Axis.y}.array(:); H2d{Axis.z}.array(:)]));
				opts.cmax = cmax;
				for w = Axis.elems
					figure;
					set(gcf, 'units','normalized','position',[(int(w)-1)/3 0 1/3 1/3]);
					vis2d(H2d{w}, obj_array, opts);
					drawnow;
				end
			end
		end

		if is_modalsrc
			pm.mark('distributed source visualization');
		end
		fprintf('%s finishes (inspection only).\n\n', mfilename);
		E = {};
		H = {};
	else  % inspect_only == false
		if isequal(solveropts.method, 'inputfile')
			% Define eps_node_array at vertices of the E-field edges.
			chkarg(isiso, 'anisotropic materials are not supported for solveropts.method == ''inputfile''.');

			if ge == GT.prim
				eps_node = eps_node{Axis.x};  % ad hoc solution

				% Old implementation (without anisotropy support) starts
				% here.
				eps_node_array = eps_node.data_expanded();  % (Nx+2) x (Ny+2) x (Nz+2)
				eps_node_array = eps_node_array(1:end-1, 1:end-1, 1:end-1);  % (Nx+1) x (Ny+1) x (Nz+1)

				% The below line is an ad hoc solution for the error.  It
				% is just to pass an array with correct size to
				% write_input(), but eps_node_array is not used in
				% write_input() currently.
				eps_node_array = eps_node_array(1:end-1, 1:end-1, 1:end-1);  % Nx x Ny x Nz
%
% 				% (Anisotropy support? begins)
% % 				eps_node_cell = eps_cell;
% % 				for w = Axis.elems
% % 					ind_g = {':',':',':'};
% % 					if grid3d.bc(w) == BC.p
% % 						ind_g{w} = grid3d.N(w);
% % 					else
% % 						ind_g{w} = 1;
% % 					end
% % 					eps_node_cell{w} = cat(int(w), eps_node_cell{w}(ind_g{:}), eps_node_cell{w});
% % 				end
% % 				eps_node_cell{Axis.x} = 2./(1./eps_node_cell{Axis.x}(1:end-1, :, :) + 1./eps_node_cell{Axis.x}(2:end, :, :));
% % 				eps_node_cell{Axis.y} = 2./(1./eps_node_cell{Axis.y}(:, 1:end-1, :) + 1./eps_node_cell{Axis.y}(:, 2:end, :));
% % 				eps_node_cell{Axis.z} = 2./(1./eps_node_cell{Axis.z}(:, :, 1:end-1) + 1./eps_node_cell{Axis.z}(:, :, 2:end));
% %
% % 				eps_node_array = (eps_node_cell{Axis.x} + eps_node_cell{Axis.y} + eps_node_cell{Axis.z})./3;
% %
% % 				eps_node_array = (eps_node_array(1:end-1,1:end-1,1:end-1) + eps_node_array(2:end,1:end-1,1:end-1) ...
% % 							+ eps_node_array(1:end-1,2:end,1:end-1) + eps_node_array(1:end-1,1:end-1,2:end) ...
% % 							+ eps_node_array(1:end-1,2:end,2:end) + eps_node_array(2:end,1:end-1,2:end) ...
% % 							+ eps_node_array(2:end,2:end,1:end-1) + eps_node_array(2:end,2:end,2:end))./8;
% 				% (Anisotropy support? ends)
%
% 				eps_node_array = 8./(1./eps_node_array(1:end-1,1:end-1,1:end-1) + 1./eps_node_array(2:end,1:end-1,1:end-1) ...
% 							+ 1./eps_node_array(1:end-1,2:end,1:end-1) + 1./eps_node_array(1:end-1,1:end-1,2:end) ...
% 							+ 1./eps_node_array(1:end-1,2:end,2:end) + 1./eps_node_array(2:end,1:end-1,2:end) ...
% 							+ 1./eps_node_array(2:end,2:end,1:end-1) + 1./eps_node_array(2:end,2:end,2:end));
			else
				eps_node_array = eps_node.data_original();  % Nx x Ny x Nz
			end

			write_input(solveropts.filenamebase, solveropts.eqtype, osc, grid3d, s_factor, ...
				eps_node_array, eps, mu, J, M, solveropts.F0, solveropts.tol, solveropts.maxit);

			pm.mark('input file creation');
			fprintf('%s finishes. (input file created)\n\n', mfilename);
			E = {};
			H = {};
		else  % solveropts.method ~= 'inputfile'
			if isequal(solveropts.method, 'direct')
				[E, H] = solve_eq_direct(solveropts.eqtype, solveropts.pml, osc.in_omega0(), eps, mu, s_factor, J, M, grid3d);
			elseif isequal(solveropts.method, 'iterative')
				[E, H, relres, iter, resvec] = solve_eq_iterative(solveropts.maxit, solveropts.tol, solveropts.F0, solveropts.eqtype, solveropts.pml, osc.in_omega0(), eps, mu, s_factor, J, M, grid3d);
			else  % for solvers using E-field based equation
				%% Apply spatial inversion.
% 				d_prim = grid3d.dl(:, GT.prim);
% 				d_dual = grid3d.dl(:, GT.dual);
% 				s_prim = s_factor(:, GT.prim);
% 				s_dual = s_factor(:, GT.dual);
				d_prim = flip_vec(grid3d.dl(:, GT.dual));  % GT.dual, not GT.prim
				d_dual = flip_vec(grid3d.dl(:, GT.prim));  % GT.prim, not GT.dual
				s_prim = flip_vec(s_factor(:, GT.dual));  % GT.dual, not GT.prim
				s_dual = flip_vec(s_factor(:, GT.prim));  % GT.prim, not GT.dual
				mu = flip_vec(mu);
				eps = flip_vec(eps);
				J = neg_vec(flip_vec(J));  % pseudovector

				if isequal(solveropts.method, 'gpu')
					ds_prim = mult_vec(d_prim, s_prim);
					ds_dual = mult_vec(d_dual, s_dual);
					figure;
					F0 = {zeros(grid3d.N), zeros(grid3d.N), zeros(grid3d.N)};
					[E, H] = fds(osc.in_omega0(), ...
									ds_prim, ds_dual, ...
									mu, eps, ...
									F0, J, ...
									solveropts.maxit, solveropts.tol, 'plot');
					%   norm(A2 * ((1./e) .* (A1 * y)) - omega^2 * m .* y - A2 * (b ./ (-i*omega*e))) / norm(b) % Error for H-field wave equation.
				elseif isequal(solveropts.method, 'aws')
					ds_prim = mult_vec(d_prim, s_prim);
					ds_dual = mult_vec(d_dual, s_dual);
					F0 = {zeros(grid3d.N), zeros(grid3d.N), zeros(grid3d.N)};
% 					F0 = {rand(grid3d.N), rand(grid3d.N), rand(grid3d.N)};
% 					F0 = {rand(1)*ones(grid3d.N), rand(1)*ones(grid3d.N), rand(1)*ones(grid3d.N)};
					callback = maxwell(osc.in_omega0(), ...
									ds_prim, ds_dual, ...
									mu, eps, ...
									F0, J, ...
									solveropts.maxit, solveropts.tol);
					while ~callback(); end
					[~, E, H] = callback();
				end

				E = neg_vec(flip_vec(E));  % pseudovector
				J = neg_vec(flip_vec(J));  % pseudovector
				H = flip_vec(H);
			end

			pm.mark('solution calculation');
			if isequal(solveropts.method, 'iterative')
				%% Report the iterative solution process.
				semilogy(1:length(resvec), resvec);
				fprintf('\t# of iteration steps = %d\n', iter);
				fprintf('\trelative residual error = %e\n', relres);
			end
			fprintf('solve for: %s-field\n', char(solveropts.eqtype.f));
			fprintf('%s finishes.\n\n', mfilename);
		end
	end

	if solveropts.returnAandb
%		[A, b, g_from_f] = create_eq(solveropts.eqtype, solveropts.pml, osc.in_omega0(), eps, mu, s_factor, J, M, grid3d);
		eq = MatrixEquation(solveropts.eqtype, solveropts.pml, osc.in_omega0(), eps, mu, s_factor, J, M, grid3d);

		[extra.A, extra.b] = eq.matrix_op();
		[~, ~, extra.GfromF] = eq.matrixfree_op();
	end

	if solveropts.returnDiv
		[Dive, Divm] = create_divs(ge, s_factor, grid3d);
		extra.Dive = Dive;
		extra.Divm = Divm;
	end

	% Construct Scalar3d objects.
	E_cell = cell(1, Axis.count);
	H_cell = cell(1, Axis.count);
	J_cell = cell(1, Axis.count);
	M_cell = cell(1, Axis.count);
	for w = Axis.elems
		if ~isempty(E)
			E_cell{w} = array2scalar(E{w}, grid3d, ge, FT.e, w, osc, PhysQ.E);
		end

		if ~isempty(H)
			H_cell{w} = array2scalar(H{w}, grid3d, ge, FT.h, w, osc, PhysQ.H);
		end

		J_cell{w} = array2scalar(J{w}, grid3d, ge, FT.e, w, osc, PhysQ.J);
		M_cell{w} = array2scalar(M{w}, grid3d, ge, FT.h, w, osc, PhysQ.M);
	end

	extra.grid3d = grid3d;
	extra.J = J_cell;
	extra.M = M_cell;
	extra.eps = eps_node;
	extra.mu = mu_node;
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
        % Rezolvarea ecuatiei diferentiale
y=dsolve('D2y=12*x^2-10','y(0)=4','Dy(0)=0','x')
        % Vizualizarea pe ecran a solutiei
pretty(y)
        % Reprezentarea grafica a solutiei pe intervalul dat
ezplot(y,[0,2.5])


text:
% test for subdivision curves

options.h = [1 4 6 4 1];

rep = 'results/subdivision-curve/';
if not(exist(rep))
    mkdir(rep);
end

%%%%%%% subdivision function %%%%%%%%
name = 'rand';
name = 'dirac';
name = 'curve';
name = 'square';

switch name
    case 'rand'
        n0 = 9;
        f0 = rescale(rand(1,n0), .05,.95);
    case 'dirac'
        f0 = [0 0 0 1 0 0];
    case 'curve'
        f0 = []; b = 1;
        while b==1
            clf;
            if size(f0,2)>1
                plot(f0(1,:), f0(2,:), '.-');
            end
            axis([0 1 0 1]); box on;
            [x,y,b] = ginput(1);
            f0(:,end+1) = [x;y];
        end
    case 'square'
        f0 = [0 0 1 1; 0 1 1 0];
        f0 = rescale(f0,.05,.95);
end
n0 = size(f0,2);
x0 = linspace(0,1,n0+1);

Jmax = 5; ms = 20; lw = 1.5;
for j=1:Jmax
    f = perform_curve_subdivision(f0, j, options);
    x = linspace(0,1,size(f,2)+1);
    clf;
    hold on;
    if size(f0,1)>1
        h = plot([f(1,:) f(1,1)], [f(2,:) f(2,1)], 'k.-');
    else
        h = plot(x, [f f(1)], 'k.-');
    end
    set(h, 'MarkerSize', ms);
    set(h, 'LineWidth', lw);
    if size(f0,1)>1
        h = plot([f0(1,:) f0(1,1)],[f0(2,:) f0(2,1)], 'r.--');
    else
        h = plot(x0,[f0 f0(1)], 'r.--');
    end
    set(h, 'LineWidth', lw);
    hold off;
    axis([0 1 0 1]); box on;
    saveas(gcf, [rep 'subdivision-func-' name '-' num2str(j) '.png'], 'png');
end

author: gpeyre

repo: matlab-toolboxes


text:
function out=icassoGet(sR,field,index)
%function out=icassoGet(sR,field,[index])
%
%PURPOSE
%
%Auxiliary function for obtaining various information from the
%Icasso result data structure. Using icassoGet, you can return
%information related to original data, e.g.: data matrix,
%(de)whitening matrix, total number of estimates, number of
%estimates on each round. You can also return specified estimated
%independent components (sources), and rows of demixing matrices
%from. (However, it is easier to use function icassoResult to return
%the final estimation results from the complete Icasso procedure.)
%
%EXAMPLES OF BASIC USAGE
%
%   M=icassoGet(sR,'m');
%
%returns total number of estimates.
%
%   nic=icassoGet(sR,'numOfIC')
%
%returns number of estimated IC components on each round in an Nx1
%vector (The number may differ in deflatory estimation mode due to
%convergence problems.)
%
%   r=icassoGet(sR,'round',[25 17 126]);
%
%workspace variable r contains now a 3x2 matrix where the first
%column shows the number of estimation cycle. The second column
%is the order of appearance of the estimate within the cycle:
%e.g. 2  5 : estimate 25 was 5th estimate in cycle 2
%     1 17 :          17     1st                   1
%     7  6 :         126     6th                   7
%
%   W=icassoGet(sR,'demixingmatrix',[25 17 126]);
%
%returns the demixing matrix rows for estimates 25, 17, and 126.
%
%   s=icassoGet(sR,'source');
%
%returns _all_ M estimated independent components (sources), i.e.,
%estimates from all resampling cycles.
%
%INPUT
%
%[An argument in brackets is optional. If it isn't  given or it's
% an empty matrix/string, the function will use a default value.]
%
% sR    (struct) Icasso result data structure (from icassoEst).
% field (string) (case insensitive) determines what is returned:
%    'data'           the original data centered
%    'N'              the number of estimation cycles
%    'M'              the total number of estimates available;
%                     equal to sum(icassoGet(sR,'numofic'));
%    'numofic'        number of ICs extracted on each cycle (Nx1
%                     vector)
%    'odim'           original data dimension
%    'rdim'           reduced data dimension (after PCA)
%    'round'          output is an Mx2 matrix that identifies from
%                     which estimation  round each estimate
%                     originates:
%                     out(i,1) is the estimation round for
%                     estimate i,
%                     out(i,2) is the ordinal number of the
%                     estimate within the round.
%    'source'         estimated source signals
%    'dewhitemat'     dewhitening matrix for the original data
%    'whitemat'       dewhitening matrix for the original data
%    'demixingmatrix' demixing matrix rows ('W' works as well)
%
% [index] (vector of integers) specifies which estimates to return
%          Applies only for 'demixingmatrix', and 'source'
%          Default: leaving this argument out corresponds to
%          selecting all estimates, i.e., giving [1:M] as index
%          vector.
%
%OUTPUT
%
% out   (type and meaning vary according to input argument 'field')
%
%SEE ALSO
% icassoResult

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% ver 1.2 100105

M=size(sR.index,1);

if nargin<3,
   index=1:M;
end

switch lower(field)
 case 'm'
   out=size(sR.index,1);
 case 'n'
  out=length(sR.A);
 case 'odim'
  out=size(sR.whiteningMatrix,2);
 case 'rdim'
  out=size(sR.whiteningMatrix,1);
 case 'numofic'
  for i=1:length(sR.A),
    out(i,1)=size(sR.A{i},2);
  end
 case {'w','demixingmatrix'}
  out=getDeMixing(sR,index);
 case 'data'
  out=sR.signal;
 case {'dewhitemat'}
  out=sR.dewhiteningMatrix;
 case {'whitemat'}
  out=sR.whiteningMatrix;
 case 'source'
  out=getSource(sR,index);
 case 'round'
  out=sR.index(index,:)
 otherwise
   error('Unknown operation.');
end

%function A=getMixing(sR,index);
%
%function A=getMixing(sR,index);
%

% reindex
%index2=sR.index(index,:);
%N=size(index2,1); A=[];
%for i=1:N,
%   A(:,i)=sR.A{index2(i,1)}(:,index2(i,2));
%end

function W=getDeMixing(sR,index);
%
%function W=getDeMixing(sR,index);
%

% reindex
index2=sR.index(index,:);
N=size(index2,1); W=[];
for i=1:N,
   W(i,:)=sR.W{index2(i,1)}(index2(i,2),:);
end

function S=getSource(sR,index)
%function S=getSource(sR,index)
%

X=sR.signal;
W=getDeMixing(sR,index);
S=W*X;

%% Old stuff
%function dWh=getDeWhitening(sR,index);
%
%function dWh=getDeWhitening(sR,index);
%
%
%index=sR.index(index,:);
%Nindex=size(index,1); W=[];
%for i=1:Nindex,
%   dWh(:,i)=sR.dewhiteningMatrix{1}(:,index(i,2));
%end

%function W=getWhitening(sR,index);
%
%function W=getWhitening(sR,index);
%
%
%index=sR.index(index,:);
%Nindex=size(index,1); W=[];
%for i=1:Nindex,
%   W(:,i)=sR.whiteningMatrix{index(i,1)}(:,index(i,2));
%end


text:
%CentralCamera.flowfield Optical flow
%
% C.flowfield(V) displays the optical flow pattern for a sparse grid
% of points when the camera has a spatial velocity V (6x1).
%
% See also QUIVER.

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

function flowfield(cam, vel)
    a = 50:100:1000;
    [U,V] = meshgrid(a, a);
    du=[]; dv=[];
    for i=1:numcols(U)
        for j=1:numrows(U)
            pdot = cam.visjac_p( [U(i,j); V(i,j)], 2 ) * vel(:);
            du(i,j) = pdot(1); dv(i,j) = pdot(2);
        end
    end

    quiver(U, V, du, dv, 0.4)
    axis equal
        axis([1 cam.npix(1) 1 cam.npix(2)]);

    set(gca, 'Ydir', 'reverse');
    xlabel('u (pixels)');
    ylabel('v (pixels)');
    grid

author: petercorke

repo: machinevision-toolbox-matlab


text:
% AssignmentToIndex Convert assignment to index.
%
%   I = AssignmentToIndex(A, D) converts an assignment, A, over variables
%   with cardinality D to an index into the .val vector for a factor.
%   If A is a matrix then the function converts each row of A to an index.
%
%   See also IndexToAssignment.m and SampleFactors.m
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function I = AssignmentToIndex(A, D)

D = D(:)'; % ensure that D is a row vector
if (any(size(A) == 1)),
    I = cumprod([1, D(1:end - 1)]) * (A(:) - 1) + 1;
else
    I = sum(bsxfun(@times, A - 1, cumprod([1, D(1:end - 1)])), 2) + 1;
end;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function test_bug2443

% WALLTIME 00:20:00
% MEM 4gb
% DEPENDENCY ft_multiplotER

% get some data
filename = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_neuromag306.mat');
load(filename);

data_pre  = data;
data_post = data;

cfg = [];
cfg.method = 'mtmfft';
cfg.output = 'powandcsd'; % the CSD is needed for source reconstruction
cfg.taper  = 'dpss';
cfg.foi    = 5:5:100;
cfg.tapsmofrq = 10;       % we apply plenty of frequency smoothing
freq_pre  = ft_freqanalysis(cfg, data_pre);
freq_post = ft_freqanalysis(cfg, data_post);

cfg = [];
cfg.layout      = 'neuromag306planar.lay';
cfg.parameter   = 'powspctrm';
ft_multiplotER(cfg, freq_pre, freq_post);


text:
function [V,E,F,H] = readPOLY_pyramid(filename)
  % READPOLY_PYRAMID Read .poly file of pyramid/svr output
  %
  % Inputs:
  %   filename  path to .poly file
  % Outputs:
  %   V  #V list of vertices
  %   E  #E by 2+boundary_markers list of segment indices, indexing V, and
  %     optional boundary markers
  %   F  #F struct containing polygon information arrays
  %     .facets  a #facets list of facets,  each facet is a polygon
  %       **NOTE: facets index E *not* V, contrary to typical (V,F) meshes and
  %       contrary to writePOLY_tetgen prototype
  %     .boundary_markers a #facets list of boundary_markers
  %     OR
  %   F  #F by constant-degree+boundary_markers  list of facets
  %   H  #H by dim list of holes (not supported)
  %

  fp = fopen(filename,'r');
  line = eat_comments(fp,'#');
  [Vhead,count] = sscanf(line,'%d %d %d %d');
  assert(count==4);
  [V,count] = fscanf(fp,'%g',[(1+sum(Vhead(2:end))) Vhead(1)]);
  V = V(2:Vhead(2)+1,:)';
  assert(count == (Vhead(1)*(1+sum(Vhead(2:end)))));
  line = eat_comments(fp,'#');
  [Ehead,count] = sscanf(line,'%d %d');
  assert(count == 2);
  [E,count] = fscanf(fp,'%d',[3+Ehead(2) Ehead(1)]);
  % only segments are supported
  BME = E(4:3+Ehead(2),:)';
  E = E(2:3,:)';
  offset = 1-min(min(E(:,1)),1);
  if offset
    warning('Offseting indices by %d',offset);
    E = E+offset;
  end
  assert(count == prod([3+Ehead(2) Ehead(1)]));
  line = eat_comments(fp,'#');
  [Fhead,count] = sscanf(line,'%d %d');
  % Not sure if boundary markers are allowed here
  if count == 1
    Fhead(2) = 0;
  end
  % preallocate
  %F = repmat(struct('facets',[],'boundary_markers',[]),Fhead(1),1);
  F.facets = cell(Fhead(1),1);
  F.boundary_markers = cell(Fhead(1),1);
  min_F = inf;
  for f = 1:Fhead(1)
    line = eat_comments(fp,'#');
    [EF,count] = sscanf(line,'%d');
    min_F = min(min_F,EF(1));
    assert(count >= 2);
    assert(count == 2+EF(2));
    F.facets{f} = EF(3:end);
    % not supported
    F.boundary_markers{f} = [];
  end

  fs = cell2mat(cellfun(@size,F.facets,'UniformOutput',false));
  % try to collapse F
  if all(fs(:,1) == 1) && all(fs(:,2) == fs(1,2))
      F = cell2mat(F.facets);
  end

  offset = 1-min(min_F,1);
  if offset
    warning('Offseting indices by %d',offset);
    F.facets = cellfun(@plus,F.facets, ...
      mat2cell(offset*ones(numel(F.facets),1),ones(numel(F.facets),1),1), ...
      'UniformOutput',false);
  end

  line = eat_comments(fp,'#');
  [Hhead,count] = sscanf(line,'%d %d');
  assert(Hhead(1) == 0);
  line = eat_comments(fp,'#');
  if ~isempty(line)
    [Rhead,count] = sscanf(line,'%d %d');
    assert(Rhead(1) == 0);
  end

  fclose(fp);

end

author: alecjacobson

repo: gptoolbox


text:
function [L,L_u,L_k] = pix2PluckerRay(k,u)

% PIX2PLUCKERRAY  Plucker ray from optical center through pixel.
%   PIX2PLUCKERRAY(K,U) is a Plucker line passing over pixel U and the
%   optical center, in a pin-Hole camera of intrinsic vector K.
%
%   [R, R_k, R_u] = PIX2PLUCKERRAY(...) returns the Jacobians wrt K and U.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1

    n = [0;0;0];
    v = invPinHole(u,1,k);

    L = [n;v];

else

    n = [0;0;0];
    [v,V_u,~,V_k] = invPinHole(u,1,k);

    L = [n;v];
    L_u = [zeros(3,2);V_u];
    L_k = [zeros(3,4);V_k];

end

return

%%
syms u1 u2 u0 v0 au av real
k = [u0;v0;au;av];
u = [u1;u2];

[L,L_u,L_k] = pix2PluckerRay(k,u);

simplify(L_u - jacobian(L,u))
simplify(L_k - jacobian(L,k))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb