text:
%                           runBenchmarkSuccessRate.m
%
% This example benchmarks algorithms based on their ability to reconstruct
% a synthetic signal (random Gaussian) using synthetic measurements
% (random Gaussian).  The benchmark shows how the
% different methods behave as the number of measurements increases.  The
% y-axis plots the success rate (rate of exact signal recovery).  This is
% done by setting params.policy='successrate'.
%
% The algorithms currently used in this implementation are all instances of
% PhaseMax, but with different levels of accuracy in the initializer.  To
% control the level of initialization accuracy, we set the initializer to
% "angle", and specify the angle between the true signal, and the initial
% guess.
%
% This script does the following:
%
% 1. Set up parameters and create a list of algorithm structs.
%
% 2. Invoke the general benchmark function benchmarkPR. A graph of errors
% (under specified error metrics) of different algorithms will be shown.
%
% The benchmark program compares the performance of specified algorithms on
% 1D gaussian measurements that have different m/n ratio(i.e. the ratio
% between the number of measurements and the number of unknowns).
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

clc
clear
close all

%% 1.Set up parameters
% Choose x label (values shown on the x axis of the benchmark plot) and
% y label (values shown on the y-axis). The value on the x axis is the m/n
% ratio. The value on the y axis is 'reconerror', which is the relative
% 2-norm difference between the true and recovered signal.
xitem = 'm/n';
xvalues = 2:.5:12;
yitem = 'reconerror';

% Choose Dataset and set up dataSet '1DGaussian' specific parameters
dataSet = '1DGaussian';

% Set up general parameters
params.verbose = false;
params.numTrials = 20;         % run several random trials for each scenario, and report average results
params.n = 500;                 % num of unknown elements
params.isComplex = true;        % use complex matrices? or just stick to real?
params.policy = 'successrate';	% use the successrate
params.successConstant = 1e-4;

%  Each of these algorithm is an instance of PhaseMax.  However, they each
%  are initialized with starting points of different accuracies.  The
%  "angle" initializer grabs the "initAngle" entry from the options, and
%  produces an initializer that makes this angle with the true signal.
pmax25 = struct('algorithm','phasemax','initMethod','angle');
pmax25.tol=1e-6;
pmax25.initAngle=25/360*2*pi;

pmax36 = struct('algorithm','phasemax','initMethod','angle');
pmax36.tol=1e-6;
pmax36.initAngle=36/360*2*pi;

pmax45 = struct('algorithm','phasemax','initMethod','angle');
pmax45.tol=1e-6;
pmax45.initAngle=45/360*2*pi;

% Grab your pick of algorithms.
algorithms = {pmax25,pmax36,pmax45};

% Run benchmark
results = benchmarkSynthetic(xitem, xvalues, yitem, algorithms, dataSet, params);

author: tomgoldstein

repo: phasepack-matlab


text:
function [pyra, model_dp] = gdetect_pos_prepare_c(im, model, valid)
% Prepare a set of foreground examples in the same image for processing
% with gdetect_pos.m.
%   [pyra, model_dp] = gdetect_pos_prepare(im, model, boxes, fg_overlap)
%
% Return values
%   pyra          Feature pyramid for image im
%   model_dp      Model augmented with dynamic programming tables
%
% Arguments
%   im            Foreground image with one or more foreground examples
%   model         Object model

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% get feature pyramid
pyra = featpyramid(im, model);
pyra = project_pyramid(model, pyra);

pyra.valid_levels(:) = false;
pyra.valid_levels([valid.l valid.l-model.interval]) = true;

% compute dynamic programming tables (stored in model_dp)
model_dp = gdetect_dp(pyra, model);

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [rec,prec,ap] = VOCevallayout_pr(VOCopts,id,draw)

% load test set
[imgids,objids]=textread(sprintf(VOCopts.layout.imgsetpath,VOCopts.testset),'%s %d');

% hash image ids
hash=VOChash_init(imgids);

% load ground truth objects

tic;
n=0;
np=zeros(VOCopts.nparts,1);
for i=1:length(imgids)
    % display progress
    if toc>1
        fprintf('layout pr: load %d/%d\n',i,length(imgids));
        drawnow;
        tic;
    end

    % read annotation
    r=PASreadrecord(sprintf(VOCopts.annopath,imgids{i}));

    % extract object
    n=n+1;
    o=r.objects(objids(i));
    gt(n)=o;

    for j=1:numel(o.part)
        c=strmatch(o.part(j).class,VOCopts.parts,'exact');
        np(c)=np(c)+1;
    end
end

% load results

fprintf('layout pr: loading results\n');
xml=VOCreadxml(sprintf(VOCopts.layout.respath,id));

% test detections by decreasing confidence

[t,si]=sort(-str2double({xml.results.layout.confidence}));
nd=numel(si);

det=false(n,1);

ptp=[];
pfp=[];
pc=[];

for di=1:nd

    % display progress
    if toc>1
        fprintf('layout pr: compute: %d/%d\n',di,nd);
        drawnow;
        tic;
    end

    % match result to ground truth object
    d=xml.results.layout(si(di));
    ii=VOChash_lookup(hash,d.image);
    oi=ii(objids(ii)==str2num(d.object));

    if isempty(oi)
        warning('unrecognized layout: image %s, object %s',d.image,d.object);
        continue
    end

    if det(oi)
        warning('duplicate layout: image %s, object %s',d.image,d.object);
        continue
    end
    det(oi)=true;
    o=gt(oi);

    % assign parts to ground truth parts

    gtd=false(numel(o.part),1);
    da=zeros(numel(d.part),1);
    dc=zeros(numel(d.part),1);
    for i=1:numel(d.part)
        dc(i)=strmatch(d.part(i).class,VOCopts.parts,'exact');
        bb=str2double({d.part(i).bndbox.xmin d.part(i).bndbox.ymin ...
                       d.part(i).bndbox.xmax d.part(i).bndbox.ymax});

        ovmax=-inf;
        for j=1:numel(o.part)
            if strcmp(d.part(i).class,o.part(j).class)
                bbgt=o.part(j).bbox;
                bi=[max(bb(1),bbgt(1))
                    max(bb(2),bbgt(2))
                    min(bb(3),bbgt(3))
                    min(bb(4),bbgt(4))];
                iw=bi(3)-bi(1)+1;
                ih=bi(4)-bi(2)+1;
                if iw>0 & ih>0
                    % compute overlap as area of intersection / area of union
                    ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
                       (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
                       iw*ih;
                    ov=iw*ih/ua;
                    if ov>ovmax
                        ovmax=ov;
                        jmax=j;
                    end
                end
            end
        end
        if ovmax>=VOCopts.minoverlap && ~gtd(jmax)
            da(i)=jmax;
            gtd(jmax)=true;
        end
    end

    ptp=[ptp ; da~=0];
    pfp=[pfp ; da==0];
    pc=[pc ; dc];
end

% evaluate each part type

for i=1:VOCopts.nparts

    % compute precision/recall

    fpi=cumsum(pfp(pc==i));
    tpi=cumsum(ptp(pc==i));
    v=tpi+fpi>0;
    rec{i}=tpi(v)/np(i);
    prec{i}=tpi(v)./(fpi(v)+tpi(v));

    ap{i}=VOCap(rec{i},prec{i});

    if draw
        % plot precision/recall
        subplot(VOCopts.nparts,1,i);
        plot(rec{i},prec{i},'-');
        grid;
        xlabel 'recall'
        ylabel 'precision'
        title(sprintf('subset: %s, part: %s, AP = %.3f',VOCopts.testset,VOCopts.parts{i},ap{i}));
    end
end

author: quantombone

repo: exemplarsvm


text:
function [ mH ] = GenerateToeplitzConvMatrix( vH, numElementsData, convShape )
% ----------------------------------------------------------------------------------------------- %
%[ mH ] = GenerateToeplitzConvMatrix( vH, numElementsData, convShape )
% Generates matrix 'mH' such that mH * vX = conv(vX, vH, 'convShape').
% Input:
%   - vH                -   Input Vector.
%                           The coeeficients of the filter.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numElementsData   -   Number of Elements of the Data.
%                           Number of elements in the vector to convolve.
%                           In the equation above it is 'length(vX)'.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - convShape         -   Convolution Shape.
%                           Sets the convolution type. Using same
%                           definitions as in MATLAB 'conv()' function.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mH                -   Toeplitz Convolution Matrix.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  A
% Remarks:
%   1.  B
% Known Issues:
%   1.  C
% TODO:
%   1.  D
% Release Notes:
%   -   1.0.001     17/01/2020  Royi Avital
%       *   Direct use of 'toeplitz()' to create the
%           'CONVOLUTION_SHAPE_SAME' case.
%   -   1.0.000     11/01/2020  Royi Avital
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

numTaps = length(vH);

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        vR = [vH(1); zeros(numElementsData - 1, 1)];
        vC = [vH(:); zeros(numElementsData - 1, 1)];
    case(CONVOLUTION_SHAPE_SAME)
        sepIdx = ceil((numTaps + 1) / 2);
        vC = zeros(numElementsData, 1);
        vR = zeros(numElementsData, 1);
        vC(1:numTaps - sepIdx + 1) = vH(sepIdx:numTaps);
        vR(1:sepIdx) = vH(sepIdx:-1:1);
    case(CONVOLUTION_SHAPE_VALID)
        vR = [flip(vH(:), 1); zeros(numElementsData - numTaps, 1)];
        vC = [vH(end); zeros(numElementsData - numTaps, 1)];
end

mH = toeplitz(vC, vR);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [MB] = B2MB(B)
% Convert computery things from bytes to megabytes.
% Chad A. Greene 2012
MB = B*2^-20 ;

author: Sable

repo: mcbench-benchmarks


text:
% COPYRIGHT
% This file is part of ORCA: https://github.com/ayrna/orca
% Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
% Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
% Copyright:
%     This software is released under the The GNU General Public License v3.0 licence
%     available at http://www.gnu.org/licenses/gpl-3.0.html
% RUN SINGLE TESTS
% This script runs single tests using the ORCA API. The toy dataset is used
% with a 1 -hold-out experimental. Reference performance is compared. To
% add an script just place it in the test/singletests folder following
% other examples
%
clear;
addpath Algorithms
addpath Measures
addpath Utils

tests_dir = ['tests' filesep 'singletests'];
eval(['addpath ' tests_dir])

files = dir([tests_dir filesep '*.m']);

% Load the different partitions of the dataset
load ../exampledata/1-holdout/toy/matlab/train_toy.0
load ../exampledata/1-holdout/toy/matlab/test_toy.0

% "patterns" refers to the input variables and targets to the output one
train.patterns = train_toy(:,1:end-1);
train.targets = train_toy(:,end);
test.patterns = test_toy(:,1:end-1);
test.targets = test_toy(:,end);

for i=1:length(files)
    % Clear (almost) all the variables to avoid issues with some mex files in Octave +4.2.2
    if (exist ('OCTAVE_VERSION', 'builtin') > 0)
      clear -x tests_dir files i train test
    else
		  clearvars -except tests_dir files i train test
    end
    fprintf('==============================================\n');
    cmd = files(i).name;
    fprintf('Processing %s...\n', cmd);
    eval(cmd(1:end-2))
end

% error is called in individual tests, so if code reach this point all the
% tests have been run successfully
fprintf('\nAll tests ended successfully\n')

rmpath Algorithms
rmpath Measures
rmpath Utils
eval(['rmpath ' tests_dir])

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function [BSL,BKG] = BSLcalcShedsVis(Wsheds,nSheds,IMAGE)
%"BSLcalcShedsVis"
%   Strips away low values watersheds and calculates BSL
%   -- Returns BSL and Background estimates after each water shed is
%   removed and returns vectors of these values
%
% CRS 07/17/13
%
%Usage:
%   [BSL,BKG] = BSLcalcShedsVis(Wsheds,nSheds,IMAGE)
%       Wsheds  = Struture that holds:
%          Wsheds.PET    = PET VOI
%          Wsheds.Shed   = Watersheds ordered from low mean uptake to high
%          Wsheds.voxVol = voxel volume
%       nSheds = total number of sheds
%       Image  = Flag for display of smoothed BSL vector
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.%
%
%
%
%% BSL calc for all sheds
if (IMAGE == 1)
    figure,
end
bsl = zeros([nSheds 1]);
bkg = zeros([nSheds 1]);

PT = Wsheds.PET;
Shed = Wsheds.Shed;
voxVol = Wsheds.voxVol;
[pX pY pZ] = size(PT);

% build shed masks
ShedTmp = zeros(size(PT));
indNZ = find(Shed);
Shed(isnan(Shed)) = 0;
Shed(isinf(Shed)) = 0;
ShedTmp(indNZ) = Shed(indNZ);

SUVcut = 0.2;
ShedTmp(ShedTmp < SUVcut) = 0;
PT(ShedTmp <= 0) = 0;

% Shed Trimming
k = 0;
while ( sum(ShedTmp(:)) > 0 )
    voxVec = [];

    minShed = min(nonzeros(ShedTmp(:)));
    indS = find(ShedTmp == minShed);
    ShedTmp(indS) = 0;
    PT(indS) = 0;
    if (isempty(indS) == 1)
        continue
    end
    k = k + 1;
    voxVec = nonzeros(PT);

    if (sum(voxVec) <=0)
        bsl(k) = -1;
        bkg(k) = -1;
        BKG = bkg(1:k);
        BSL = bsl(1:k);
        if (IMAGE ~= 0)
            hold off
        end
        return
    end

    nVox = numel(voxVec);
    SUViqr = iqr(voxVec);
    IVHBinWidth = 2 * SUViqr * nVox^(-1/3);
    %%%
    if (numel(voxVec) < 5)
        BKG = bkg(1:k);
        BSL = bsl(1:k);
        if (IMAGE ~= 0)
            hold off
        end
        return
    end
    %%%
    [binHistV, volHistV]  = doseHist(voxVec, voxVol*ones(size(voxVec)), IVHBinWidth);
    nBins = numel(binHistV);

    [Hmax iBkg] = max(volHistV);
    Hbkg = binHistV(iBkg);
    if (iBkg > 1 && iBkg < numel(binHistV))
        Hbkg = ( volHistV(iBkg-1)*binHistV(iBkg-1) ...
            + volHistV(iBkg)*binHistV(iBkg) ...
            + volHistV(iBkg+1)*binHistV(iBkg+1) ) ...
            / ( volHistV(iBkg-1) + volHistV(iBkg) + volHistV(iBkg+1) );
    else
        if (iBkg < numel(binHistV))
            Hbkg = ( volHistV(iBkg)*binHistV(iBkg) ...
                + volHistV(iBkg+1)*binHistV(iBkg+1) ) ...
                / ( volHistV(iBkg) + volHistV(iBkg+1) );
        else
            Hbkg = ( volHistV(iBkg-1)*binHistV(iBkg-1) ...
                + volHistV(iBkg)*binHistV(iBkg) ) ...
                / ( volHistV(iBkg-1) + volHistV(iBkg) );
        end
    end
    SUVcut = Hbkg/3;

    voxVec(voxVec < SUVcut) = 0;
    voxVec = nonzeros(voxVec);
    ShedTmp(ShedTmp < SUVcut) = 0;
    PT(ShedTmp < SUVcut) = 0;

    nVox = numel(voxVec);
    SUViqr = iqr(voxVec);
    IVHBinWidth = 2 * SUViqr * nVox^(-1/3);

    [binHistV, volHistV]  = doseHist(voxVec, voxVol*ones(size(voxVec)), IVHBinWidth);
    nBins = numel(binHistV);
    if (IMAGE == 1)
        colorWidth = 2;
        colorPlot = mod(k,7);
        switch colorPlot
            case (1)
                plot(binHistV,volHistV,'-c','LineWidth',colorWidth)
            case (2)
                plot(binHistV,volHistV,'-g','LineWidth',colorWidth)
            case (3)
                plot(binHistV,volHistV,'-y','LineWidth',colorWidth)
            case (4)
                plot(binHistV,volHistV,'-r','LineWidth',colorWidth)
            case (5)
                plot(binHistV,volHistV,'-k','LineWidth',colorWidth)
            case (6)
                plot(binHistV,volHistV,'-m','LineWidth',colorWidth)
            otherwise
                plot(binHistV,volHistV,'-b','LineWidth',colorWidth)
        end
        drawnow, hold on
    end
    if (IMAGE == 2)
        plot3(binHistV,-k*ones(numel(binHistV)),volHistV), drawnow
        hold on
    end
    if (IMAGE == 3)
        colorWidth = 4;
        colorPlot = mod(k,7);
        switch colorPlot
            case (1)
                plot(binHistV,volHistV,'-c','LineWidth',colorWidth)
            case (2)
                plot(binHistV,volHistV,'-g','LineWidth',colorWidth)
            case (3)
                plot(binHistV,volHistV,'-y','LineWidth',colorWidth)
            case (4)
                plot(binHistV,volHistV,'-r','LineWidth',colorWidth)
            case (5)
                plot(binHistV,volHistV,'-k','LineWidth',colorWidth)
            case (6)
                plot(binHistV,volHistV,'-m','LineWidth',colorWidth)
            otherwise
                plot(binHistV,volHistV,'-b','LineWidth',colorWidth)
        end
        drawnow, hold on
    end
    %% Gaussian FITS
    fitOptions = optimset('Display','off');
    lBins = ceil(SUViqr/IVHBinWidth);

    plusBins  = iBkg + lBins;
    if (plusBins > numel(binHistV))
        plusBins = numel(binHistV);
    end
    minusBins = iBkg  - lBins;
    if (minusBins < 1)
        minusBins = 1;
    end
    if (numel(minusBins:plusBins) < 3)
        if (minusBins+2 <= numel(binHistV))
            plusBins = minusBins + 2;
        else
            minusBins = plusBins - 2;
        end
    end
    if (minusBins < 1)
        BKG = bkg(1:k);
        BSL = bsl(1:k);
        if (IMAGE ~= 0)
            hold off
        end
        return
    end

    sL2 = [];
    initVar = [ Hmax*sqrt(2*pi())*Hbkg Hbkg SUViqr ];
    lBound =  [ 0.01  0.01  0.01  ];
    uBound =  [ inf   inf   inf  ];

    sL2  = lsqnonlin(@(X) (...
        ( volHistV(minusBins:plusBins) - X(1)*normpdf(binHistV(minusBins:plusBins),X(2),X(3)) ) ...
        / sum(volHistV(minusBins:plusBins)) ).*sqrt(volHistV(minusBins:plusBins)), ...
        initVar,lBound,uBound,fitOptions);

    estGauss2 = sL2(1)*normpdf(binHistV(:),sL2(2),sL2(3));

    bkg(k) = sL2(2);

    %% Volume and TLG Estimation
    iEst2 = round( 2*sL2(3)/IVHBinWidth + sL2(2)/IVHBinWidth );
    if  (sL2(2) <= 0)
        iEst2 = round( 2*sL2(3)/IVHBinWidth + iBkg );
    end
    volGauss2 = estGauss2(iEst2:end);

    % Gauss volumes
    vecG2 = zeros(size(volGauss2));
    vecG2 = volHistV(iEst2:end)' - volGauss2;
    vecG2(vecG2 < 0) = 0;
    volG2 = sum(vecG2);

    bsl(k) = vecG2'*binHistV(iEst2:end)';
    if (bsl(k) <= 0 && k > nSheds/10)
        BKG = bkg(1:k);
        BSL = bsl(1:k);
        if (IMAGE ~= 0)
            hold off
        end
        return;
    end

end
BKG = bkg(1:k);
BSL = bsl(1:k);
if (IMAGE ~= 0)
    hold off
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function scalhand(H,scale)
    %
    %    scalhand(Handle,scalefactor)
    %
    %    Scale object Handle by amount scalefactor
    %
    %    See also ROTATE, TRNSLATE
    %
    %    Richard G. Cobb    3/96
    %
    origin=findcntr(H);
    %
    xold=get(H,'Xdata');
    xnew=scale*(xold-origin(1))+origin(1);
    yold=get(H,'Ydata');
    ynew=scale*(yold-origin(2))+origin(2);

    zold=get(H,'Zdata');
    znew=scale*(zold-origin(3))+origin(3);

    set(H,'Xdata',xnew,'Ydata',ynew,'Zdata',znew)

author: CelsoReyes


text:
%SE3 Lift SE(2) transform to SE(3)
%
% T3 = SE3(T2) returns a homogeneous transform (4x4) that represents
% the same X,Y translation and Z rotation as does T2 (3x3).
%
% See also SE2, SE2.SE3, TRANSL, ROTX.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math
function T = se3(x)

    if all(size(x) == [3 3])
        T = [x(1:2,1:2) [0 0]' x(1:2,3); 0 0 1 0; 0 0 0 1];
     end

author: petercorke

repo: spatialmath-matlab


text:
function L = listdijkstra(L,W,s,d)

index=size(W,1);
while index>0
    if W(2,d)==W(size(W,1),d)
        L=[L s];
        index=0;
    else
        index2=size(W,1);
        while index2>0
            if W(index2,d)<W(index2-1,d)
                L=[L W(index2,1)];
                L=listdijkstra(L,W,s,W(index2,1));
                index2=0;
            else
                index2=index2-1;
            end
            index=0;
        end
    end
end


text:
classdef LiftCoefficientModelEnum < matlab.mixin.SetGet
    %LiftCoefficientModelEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        KSPCylinder('KSP Cylinder',"Models the vehicle as a circular cylinder with constant radius and height whose long axis is the body X axis.  Uses the KSP 'lifting body' lift coefficient curves in physics.cfg.");
    end

    properties
        name(1,:) char
        desc(1,1) string
    end

    methods
        function obj = LiftCoefficientModelEnum(name, desc)
            obj.name = name;
            obj.desc = desc;
        end
    end

    methods(Static)
        function [listBoxStr, m] = getListBoxStr()
            m = enumeration('LiftCoefficientModelEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            m = m(I);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function tests = test_ft_preproc_polyremoval

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_polyremoval

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

flag = false;

result = [];
result{end+1} = ft_preproc_polyremoval(dat, 1, 1, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 2, 1, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 3, 1, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 4, 1, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 1, 101, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 2, 101, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 3, 101, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 4, 101, 1000, flag);
result{end+1} = ft_preproc_polyremoval(dat, 1, 1, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 2, 1, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 3, 1, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 4, 1, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 1, 101, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 2, 101, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 3, 101, 900, flag);
result{end+1} = ft_preproc_polyremoval(dat, 4, 101, 900, flag);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [d,x0,xv,idx,delay_offset] = driving_function_imp_localwfs_vss(x0,xs,src,conf)
%DRIVING_FUNCTION_IMP_LOCALWFS_VSS driving signal for local WFS using focused
%sources as virtual secondary sources
%
%   Usage: [d,x0,xv,idx,delay_offset] =
%               driving_function_imp_localwfs_vss(x0,xs,src,conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [nx7]
%       xs          - position of virtual source or direction of plane
%                     wave / m [1x3]
%       src         - source type of the virtual source
%                         'pw' - plane wave (xs is the direction of the
%                                plane wave in this case)
%                         'ps' - point source
%                         'ls' - line source
%                         'fs' - focused source
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       d             - driving signals [mxn]
%       x0            - position, direction, and weights of the real secondary
%                       sources / m [nx7]
%       xv            - position, direction, and weights of the virtual
%                       secondary sources / m [mx7]
%       idx           - index of the selected sources from the original x0
%                       matrix [mx1]
%       delay_offset  - additional added delay, so you can correct it
%
%   See also: plot_sound_field, sound_field_mono_wfs
%
%   References:
%       Spors and Ahrens (2010) - "Local Sound Field Synthesis by Virtual
%       Secondary Sources", 40th Conference of the Audio Engineering Society,
%       Paper 6-3, http://www.aes.org/e-lib/browse.cfm?elib=15561
%

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);
isargsecondarysource(x0);
isargxs(xs);
isargchar(src);
isargstruct(conf);

%% ===== Configuration ==================================================
virtualconf = conf;
virtualconf.usetapwin = conf.localwfs_vss.usetapwin;
virtualconf.tapwinlen = conf.localwfs_vss.tapwinlen;
virtualconf.secondary_sources.size = conf.localwfs_vss.size;
virtualconf.secondary_sources.center = conf.localwfs_vss.center;
virtualconf.secondary_sources.geometry = conf.localwfs_vss.geometry;
virtualconf.secondary_sources.number = conf.localwfs_vss.number;
virtualconf.wfs = conf.localwfs_vss.wfs;
virtualconf.nfchoa = conf.localwfs_vss.nfchoa;
virtualconf.driving_functions = conf.localwfs_vss.driving_functions;
method = conf.localwfs_vss.method;

%% ===== Computation ====================================================
if strcmp('fs',src)
    error(['%s: %s is not a supported method source type! Try to use a', ...
        ' point source, if the source is inside the secondary source array', ...
        ' but not inside the virtual secondary source array'], ...
        upper(mfilename),src);
end

% Determine driving functions of virtual array with different sfs methods
switch method
case 'wfs'
    % === Wave Field Synthesis ===
    % Create virtual source array
    xv = virtual_secondary_source_positions(x0,xs,src,conf);
    % Secondary_source_selection
    xv = secondary_source_selection(xv, xs, src);
    % Optional tapering
    xv = secondary_source_tapering(xv,virtualconf);
    % Driving functions for virtual source array
    [dv,~,~,delay_offset] = driving_function_imp_wfs(xv,xs,src,virtualconf);
case 'nfchoa'
    % === Near-Field-Compensated Higher Order Ambisonics ===
    % Create virtual source array
    xv = secondary_source_positions(virtualconf);
    % Driving functions for virtual source array
    [dv,~,delay_offset] = driving_function_imp_nfchoa(xv,xs,src,virtualconf);
otherwise
    error('%s: %s is not a supported method for localsfs!',upper(mfilename),...
      method);
end

% Select secondary sources
[x0,idx] = secondary_source_selection(x0,xv(:,1:6),'vss');
% Driving functions for real source array
[d,delay_vss] = driving_function_imp_wfs_vss(x0,xv,'fs',dv,conf);
% add delay
delay_offset = delay_offset + delay_vss;

author: sfstoolbox

repo: sfs-matlab


text:
%% Setting of the problem
global s
pde = fonedata; % f = 1;
option.theta = 0.3;
option.estType = 'star';
option.maxIt = 18;
option.maxN = 2e4;
option.solver = 'mg';
option.tol = 1e-6;
[node,elem] = squaremesh([-1,1,-1,1],0.5);
[node,elem] = delmesh(node,elem,'x>0 & y<0');
bdFlag = setboundary(node,elem,'Dirichlet');

%% s = 0.2
s = 0.2; %#ok<*NASGU>
err1 = afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.4
s = 0.4;
err2 = afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
err3 = afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
err4 = afemfracLap(node,elem,pde,bdFlag,option);

%% save data and plot the table
% save Lshapencf err1 err2 err3 err4
load Lshapencf
figure;
plot_error_table(err1.N,err1.energyError,err2.N,err2.energyError,...
                 err3.N,err3.energyError,err4.N,err4.energyError);
saveas(gcf, 'error_Lshapencf', 'pdf')
figure;
plot_error_table(err1.N,err1.eta,err2.N,err2.eta,err3.N,err3.eta,err4.N,err4.eta);
ylabel('Estimator','interpreter','latex', 'FontSize', 22)
saveas(gcf, 'eta_Lshapencf', 'pdf')

repo: ifem


text:
function rfParams = rmPlotGUI_getRFParams(model, modelName, viewType, coords, params)
% For the retinotopy model GUI: get PRF params for a single voxel.
%
%  rmParams = rmPlotGUI_getRFParams(model, modelName, viewType, coords, params)
%
% ras, 11/2008. Trying to make the large GUI code more modular and
% more manageable.
switch modelName,
    case {'2D pRF fit (x,y,sigma, positive only)',...
          '2D RF (x,y,sigma) fit (positive only)',...
          '1D pRF fit (x,sigma, positive only)'};
        rfParams = zeros(1,6);

        rfParams(1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(3) = rmCoordsGet(viewType, model, 'sigmamajor',coords);
        rfParams(5) = rmCoordsGet(viewType, model, 'sigmaminor',coords);
        rfParams(6) = rmCoordsGet(viewType, model, 'sigmatheta',coords);

    case {'2D pRF fit (x,y,sigma_major,sigma_minor)' ...
            'oval 2D pRF fit (x,y,sigma_major,sigma_minor,theta)'};
        rfParams = zeros(1,4);

        rfParams(1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(3) = rmCoordsGet(viewType, model, 'sigmamajor',coords);
        rfParams(5) = rmCoordsGet(viewType, model, 'sigmaminor',coords);
        rfParams(6) = rmCoordsGet(viewType, model, 'sigmatheta',coords);

    case 'unsigned 2D pRF fit (x,y,sigma)';
        rfParams = zeros(1,4);

        rfParams(1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(3) = rmCoordsGet(viewType, model, 'sigmamajor',coords);

    case {'Double 2D pRF fit (x,y,sigma,sigma2, center=positive)',...
          'Difference 2D pRF fit (x,y,sigma,sigma2, center=positive)',...
          'Difference 1D pRF fit (x,sigma, sigma2, center=positive)'},
        rfParams = zeros(2,4);

        % get RF parameters
        rfParams(1,1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(1,2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(1,3) = rmCoordsGet(viewType, model, 'sigmamajor',coords);
        rfParams(2,1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(2,2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(2,3) = rmCoordsGet(viewType, model, 'sigma2major',coords);

    case {'Two independent 2D pRF fit (2*(x,y,sigma, positive only))',...
          'Mirrored 2D pRF fit (2*(x,y,sigma, positive only))',...
          'Shifted 2D pRF fit (2*(x,y,sigma, positive only))'},
        rfParams = zeros(2,3);

        % get RF parameters
        rfParams(1,1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(1,2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(1,3) = rmCoordsGet(viewType, model, 's',coords);
        rfParams(2,1) = rmCoordsGet(viewType, model, 'x02', coords);
        rfParams(2,2) = rmCoordsGet(viewType, model, 'y02', coords);
        rfParams(2,3) = rmCoordsGet(viewType, model, 's2',coords);

    case {'Sequential 2D pRF fit (2*(x,y,sigma, positive only))'},
        rfParams = zeros(2,3);

        % get RF parameters
        rfParams(1,1) = rmCoordsGet(viewType, model, 'x0', coords);
        rfParams(1,2) = rmCoordsGet(viewType, model, 'y0', coords);
        rfParams(1,3) = rmCoordsGet(viewType, model, 's',coords);
        rfParams(2,1) = rmCoordsGet(viewType, model, 'x02', coords);
        rfParams(2,2) = rmCoordsGet(viewType, model, 'y02', coords);
        rfParams(2,3) = rmCoordsGet(viewType, model, 's2',coords);

    case {'fitprf' 'css' '2D nonlinear pRF fit (x,y,sigma,exponent, positive only)' ...
            '2D nonlinear pRF fit with boxcar (x,y,sigma,exponent, positive only)'}
        rfParams = zeros(1,6);

        rfParams(1) =  rmCoordsGet(viewType, model, 'x0', coords);        % x coordinate (in deg)
        rfParams(2) = rmCoordsGet(viewType, model, 'y0', coords);         % y coordinate (in deg)
        rfParams(3) =  rmCoordsGet(viewType, model, 'sigmamajor',coords); % sigma (in deg)
        rfParams(6) =  rmCoordsGet(viewType, model, 'sigmatheta',coords); % sigma theta (0 unless we have anisotropic Gaussians)
        rfParams(7) =  rmCoordsGet(viewType, model, 'exponent'  ,coords); % pRF exponent
        rfParams(8) =  rmCoordsGet(viewType, model, 'bcomp1',    coords); % gain ?
        rfParams(5) =  rfParams(3) / sqrt(rfParams(7));                   % sigma adjusted by exponent (not for calculations - just for diplay)

    otherwise,
        error('Unknown modelName: %s',modelName{M.modelNum});
end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs, non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

% Contraction Frequency
freq = 1.0;

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
springs_info(319:end,4) = abs( cos(freq*current_time*pi) );

%NOTE: not 2*pi*ft b/c of abs()

author: nickabattista


text:
function [ nf, xf, wf, vf ] = line_fekete_monomial ( m, a, b, n, x )

%*****************************************************************************80
%
%% LINE_FEKETE_MONOMIAL computes approximate Fekete points in an interval [A,B].
%
%  Discussion:
%
%    We use the uniform weight and the monomial basis:
%
%      P(j) = x^(j-1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alvise Sommariva, Marco Vianello,
%    Computing approximate Fekete points by QR factorizations of Vandermonde
%    matrices,
%    Computers and Mathematics with Applications,
%    Volume 57, 2009, pages 1324-1336.
%
%  Parameters:
%
%    Input, integer M, the number of basis polynomials.
%
%    Input, real A, B, the endpoints of the interval.
%
%    Input, integer N, the number of sample points.
%    M <= N.
%
%    Input, real X(N), the coordinates of the sample points.
%
%    Output, integer NF, the number of Fekete points.
%    If the computation is successful, NF = M.
%
%    Output, real XF(NF), the coordinates of the Fekete points.
%
%    Output, real WF(NF), the weights of the Fekete points.
%
%    Output, real VF(NF,M), the nonsingular Vandermonde submatrix.
%
  if ( n < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LINE_FEKETE_MONOMIAL - Fatal error!\n' );
    fprintf ( 1, '  N < M.\n' );
    error ( 'LINE_FEKETE_MONOMIAL - Fatal error!' );
  end
%
%  Destroy all row vectors!
%
  x = x(:);
%
%  Moments of monomials:
%
  mom = line_monomial_moments ( a, b, m );
%
%  Form the rectangular Vandermonde matrix V for the polynomial basis.
%
  v = zeros ( m, n );

  v(1,1:n) = 1.0;
  for i = 2 : m
    v(i,1:n) = v(i-1,1:n) .* x(1:n)';
  end
%
%  Solve the system for the weights W.
%
  w = v \ mom;
%
%  Locate nonzero W's.
%
  ind = find ( w ~= 0.0 );
%
%  Retain data associated with nonzero W's.
%
  nf = length ( ind );
  xf = x ( ind );
  wf = w ( ind );
  vf = v ( :, ind );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function candoit = canGetGradient(problem)
% Checks whether the gradient can be computed for a problem structure.
%
% function candoit = canGetGradient(problem)
%
% Returns true if the gradient of the cost function can be computed given
% the problem description, false otherwise.
%
% See also: canGetCost canGetDirectionalDerivative canGetHessian

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    candoit = isfield(problem, 'grad') || isfield(problem, 'costgrad') || ...
              canGetEuclideanGradient(problem);

end

author: andrewssobral

repo: lrslibrary


text:
function [W,result] = IHR(Ws,npos,gama,afa,yuzhi)
global eta;
global TrainSetX;
global TestSetX;
global TestSetY;
global kneighbor_testp;

global W0;
global s11;
global w1;
global index1;
W0 = Ws;
w1 = Ws;

tempsum1 = 0;
for i=1:size(TestSetX,2)
    tempKN = kneighbor_testp(i,:);
    tempTuples = TestSetX(:,tempKN);
    s1 = w1'*tempTuples;
    p1 = 1./(1 + exp(-s1));
    tempsum1 = tempsum1 + (sum(p1)/size(kneighbor_testp,2)-1/(1+exp(-w1'*TestSetX(:,i))))^2;
end

tempsum2 = 0;
for i=1:size(TestSetX,2)
    tempsum2 = tempsum2 + (1/(1+exp(-w1'*TestSetX(:,i))) - 0.5)^2;
end

s1 = w1'*TestSetX;
p1 = 1./(1 + exp(-s1));

f = ((w1)'*(w1))+afa*tempsum1/size(TestSetX,2)-gama*tempsum2/size(TestSetX,2)+eta*(sum(p1)-npos)^2/size(TestSetX,2); %07-12-31
fprintf('initial value...%g\n',f);
index1 = 0;

result(1,1) = getResult(p1,TestSetY);
ff(1,1) = f;
warning off;
while index1 < 100

    temp11 = zeros(size(TestSetX,1),1);
    s11 = zeros(size(TestSetX,1),1);
    for i=1:size(TestSetX,2)
        tempKN = kneighbor_testp(i,:);
        tempTuples = TestSetX(:,tempKN);
        s1 = w1'*tempTuples;
        p1 = 1./(1 + exp(-s1));
        tmp_p1 = sum(p1)/size(kneighbor_testp,2)-1/(1+exp(-w1'*TestSetX(:,i)));
        p2 = exp(-s1)./((1 + exp(-s1)).^2);
        for k = 1:size(kneighbor_testp,2)
            temp11 = temp11 + tempTuples(:,k)*p2(1,k);
        end
        temp11 = temp11/size(kneighbor_testp,2) - (exp(-w1'*TestSetX(:,i))/(1+exp(-w1'*TestSetX(:,i)))^2)*TestSetX(:,i);
        s11 = s11 + 2*tmp_p1*temp11;
    end
    s11 = afa*s11/size(TestSetX,2); %

    s12 = zeros(size(TestSetX,1),1);
    for i = 1:size(TestSetX,2)
        s12 = s12 + 2*gama*(1/(1+exp(-w1'*TestSetX(:,i))) - 0.5)*((exp(-w1'*TestSetX(:,i)))/(1+exp(-w1'*TestSetX(:,i)))^2)*TestSetX(:,i);
    end
    s12 = s12/size(TestSetX,2);

    s1 = w1'*TestSetX;
    p1 = 1./(1 + exp(-s1));
    tmps13 = 2*eta*(sum(p1)-npos)*(exp(-w1'*TestSetX)./(1+exp(-w1'*TestSetX)).^2);
    s13 = sum(scale_cols(TestSetX,tmps13),2);
    s13 = s13/size(TestSetX,2);

    if [s11-s12+s13+2*(w1)]'*[s11-s12+s13+2*(w1)] < yuzhi^2
        break;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if index1 == 0
        d0 = -s11+s12-s13-2*(w1);
    end
    if index1 > 0
        afa0 = ([s11-s12+s13+2*(w1)]'*[s11-s12+s13+2*(w1)-temp_s11-temp_s12])/([temp_s11+temp_s12]'*[temp_s11+temp_s12]);
        d1 = [-s11+s12-s13-2*(w1)] + afa0*d0;
        d0 = d1;
    end

    temp_s11 = s11;
    temp_s12 = -s12+s13+2*(w1);

    s11 = d0;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    beta = 0;

    options = optimset('LargeScale','on');
    [x,fval] = fminunc(@objfun,beta,options);

    beta = x;
    if beta == 0
        break;
    end
    w1 = w1 + beta*d0;

    tempsum1 = 0;
    for i=1:size(TestSetX,2)
        tempKN = kneighbor_testp(i,:);
        tempTuples = TestSetX(:,tempKN);
        s1 = w1'*tempTuples;
        p1 = 1./(1 + exp(-s1));
        tempsum1 = tempsum1 + (sum(p1)/size(kneighbor_testp,2)-1/(1+exp(-w1'*TestSetX(:,i))))^2;
    end

    tempsum2 = 0;
    for i=1:size(TestSetX,2)
        tempsum2 = tempsum2 + (1/(1+exp(-w1'*TestSetX(:,i))) - 0.5)^2;
    end

    s1 = w1'*TestSetX;
    p1 = 1./(1 + exp(-s1));

    f = ((w1)'*(w1))+afa*tempsum1/size(TestSetX,2)-gama*tempsum2/size(TestSetX,2)+eta*(sum(p1)-npos)^2/size(TestSetX,2); %07-12-31

    s1 = w1'*TestSetX;
    p1 = 1./(1 + exp(-s1));
    index1 = index1+1;
    result(1,index1+1) = getResult(p1,TestSetY);
    ff(1,index1+1) = f;

    fprintf('iterating...%g value : %g...the accuracy:%g\n',index1,f,getResult(p1,TestSetY));

end
W = w1;

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function varargout = imChainPixels(img, varargin)
%IMCHAINPIXELS Chain neighbor pixels in an image to form a contour
%
%   POINTS = imChainPixels(IMG)
%   [PX, PY] = imChainPixels(IMG)
%   returns a list of points, chaining points visible in the image.
%   IMG shoud be a binary image, containing a single 8-connected loop.
%
%   Example
%     img = imread('circles.png');
%     img = imFillHoles(img);
%     bnd = imSkeleton(imBoundary(img));
%     figure; subplot(1, 3, 1); imshow(img);
%     subplot(1, 3, 2); imshow(imOverlay(img, bnd));
%     pts = imChainPixels(bnd);
%     subplot(1, 3, 3); imshow(img); hold on;
%     drawPolygon(pts, 'g', 'linewidth', 2);
%
%   See Also
%     imFillHoles, imSkeleton, bwboundaries
%
%   -----
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 01/11/2003.
%

%   HISTORY
%   2004-04-06 add to graph lib, correct x-y ordering, and add doc.
%   2012-05-16 rename from findContour to imChainPixels

% set to logical
img = img ~= 0;

% find points in the image
[pty, ptx] = find(img);

% Initialize iteration: first set the initial point
points = zeros(length(ptx), 2);
x0 = ptx(1);
y0 = pty(1);
points(1, 1) = x0;
points(1, 2) = y0;

% Then set the second point, in the neighbourhood of the first point
vois = findNeighbors(img, [x0 y0]);
x = vois(1, 1);
y = vois(1, 2);
points(2, 1) = x;
points(2, 2) = y;

% For each point, find all neighbours (should be only 2)
% Compare with point previously processed, and choose the other one.
% then shift points references, and loop until all points are processed.
for p=3:length(points)
    vois = findNeighbors(img, [x y]);
    if vois(1,1)==x0 && vois(1,2)==y0
        x0 = x;
        y0 = y;
        x = vois(2, 1);
        y = vois(2, 2);
    else
        x0 = x;
        y0 = y;
        x = vois(1, 1);
        y = vois(1, 2);
    end
    points(p, 1) = x;
    points(p, 2) = y;
end

% format results to match output
if nargout==1
    varargout{1} = points;
elseif nargout==2
    varargout{1} = points(:, 1);
    varargout{2} = points(:, 2);
end

return

function neighList = findNeighbors(img, coord)

xp = coord(1);
yp = coord(2);
neighList = [];
nv = 0;

for x = xp-1:xp+1
    if img(yp-1, x)
        nv = nv+1;
        neighList(nv, 1) = x; %#ok<AGROW>
        neighList(nv, 2) = yp-1; %#ok<AGROW>
    end
    if img(yp+1, x)
        nv = nv+1;
        neighList(nv, 1) = x; %#ok<AGROW>
        neighList(nv, 2) = yp+1; %#ok<AGROW>
    end
end

if img(yp, xp-1)
    nv = nv+1;
    neighList(nv, 1) = xp-1;
    neighList(nv, 2) = yp;
end
if img(yp, xp+1)
    nv = nv+1;
    neighList(nv, 1) = xp+1;
    neighList(nv, 2) = yp;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function numLS=leapSecsOnDay(year,month,day)
%%LEAPSECSONDAY Returns the number of leap seconds in a day for dates from
%               1972 onward. Prior to 1972, non-integer leap second numbers
%               were used, meaning that the time of day matters. From 1972
%               onward, one could have a positive or negative leap second,
%               meaning that the last minute of the day in coordinated
%               universal time (UTC) might contain 59 or 61 seconds.
%
%INPUTS: year    An integer  year in the Gregorian calendar under UTC
%                time. year>= 1960 when UTC started.
%       month    An integer month in the Gregorian calendar under UTC
%                time. 1<=month<=12
%        day     An integer day in the Gregorian calendar under UTC
%                time. Days count from 1.
%
%OUTPUTS: numLS  The number of leap seconds on the given day. This can be
%                0, -1, or 1.
%
%This function calls the cumLeapSec function for the day in question and
%for one day later to see what the difference is. The function cumLeapSec
%relies on iaudat in the International Astronomical Union's (IAU)
%Standard's of Fundamental Astronomy library. Note that leap second
%information is not available for future dates as they can not be
%accurately predicted.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(year<1972)
        error('Dates before 1972 have a more complicated leap second methodology than is used here');
    end

%Determine the number of leap on the date.
    numSec=cumLeapSec(year,month,day);
    [Jul1,Jul2]=Cal2UTC(year,month,day,0,0,0);
    %Advance the UTC day count by one day.
    Jul1=Jul1+1;
    [year,month,day]=UTC2Cal(Jul1,Jul2,true);
    %Determine the number of leap seconds on day later.
    numSecNextDay=cumLeapSec(year,month,day);

    numLS=numSecNextDay-numSec;
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:
%  Figure 7.47      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% script to generate Fig. 7.47
%
%
clf;
dp=[1 10 16 0];
np=[10];
nc=94.5*conv([1 7.98],[1 2.52]);
dc=conv([1 8.56 59.5348],[1 10.6]);
num=conv(np,nc);
den=conv(dp,dc);
dcl=[0 0 0 0 num]+den; % closed-loop denominator
t=0:.1:5;
sys=tf(num,dcl);
y=step(sys,t);
plot(t,y)
grid;
xlabel('Time (sec)');
ylabel('y');
title('Fig.7.47 Step response for SRL compensation design: continuous');

[a,b,c,d]=tf2ss(nc,dc);
% Samplig period
ts=0.1;
% Convert controller to discrete
[ad,bd,cd,dd]=c2dm(a,b,c,d,ts,'zoh');
[ncd,dcd]=ss2tf(ad,bd,cd,dd)

% response comparison of continuous and digital control.
sim('fig7_46')
figure(1)
plot(ycd(:,1),ycd(:,2))
hold on
plot(ycd(:,1),ycd(:,3),'g--')
title('Figure 7.47(a) Output Responses of Digital and Continuous Controllers')
xlabel('Time (sec)');
ylabel('y');
text(1.1,0.95,'continuous controller');
text(1,1.25,'digital controller');
nicegrid;
hold off

figure(2)
plot(ucd(:,1),ucd(:,2))
hold on
plot(ucd(:,1),ucd(:,3),'g--')
title('Figure 7.47(b)  Control Responses of Digital and Continuous Controllers')
xlabel('Time (sec)');
ylabel('u');
grid on;
text(0.1,-1.4,'continuous controller')
text(0.7,0.5,'digital controller')
nicegrid
hold off


text:
function a = wathen ( nx, ny, n )

%*****************************************************************************80
%
%% WATHEN returns a finite element matrix.
%
%  Discussion:
%
%    The Wathen matrix is a finite element matrix which is sparse.
%
%    The entries of the matrix depend in part on a physical quantity
%    related to density.  That density is here assigned random values between
%    0 and 100.
%
%    The matrix order N is determined by the input quantities NX and NY,
%    which would usually be the number of elements in the X and Y directions.
%    The value of N is
%
%      N = 3*NX*NY + 2*NX + 2*NY + 1,
%
%    and sufficient storage in A must have been set aside to hold
%    the matrix.
%
%    A is the consistent mass matrix for a regular NX by NY grid
%    of 8 node serendipity elements.
%
%    Here is an illustration for NX = 3, NY = 2:
%
%     23-24-25-26-27-28-29
%      |     |     |     |
%     19    20    21    22
%      |     |     |     |
%     12-13-14-15-16-17-18
%      |     |     |     |
%      8     9    10    11
%      |     |     |     |
%      1--2--3--4--5--6--7
%
%    For this example, the total number of nodes is, as expected,
%
%      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29
%
%  Properties:
%
%    A is symmetric positive definite for any positive values of the
%    density RHO(NX,NY), which is here given the value 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, 1987, pages 449-457.
%
%  Parameters:
%
%    Input, integer NX, NY, values which determine the size of A.
%
%    Input, integer N, the order of N, as determined by NX and NY.
%
%    Output, real A(N,N), the matrix.
%
  em =  [ ...
     6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0; ...
    -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0; ...
     2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0; ...
    -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0; ...
     3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0; ...
    -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0; ...
     2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0; ...
    -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 ];

  a(1:n,1:n) = 0.0;

  for j = 1 : ny

    for i = 1 : nx
%
%  For the element (I,J), determine the indices of the 8 nodes.
%
      node(1) = 3 * j * nx + 2 * j + 2 * i + 1;
      node(2) = node(1) - 1;
      node(3) = node(1) - 2;

      node(4) = ( 3 * j - 1 ) * nx + 2 * j + i - 1;
      node(8) = node(4) + 1;

      node(5) = ( 3 * j - 3 ) * nx + 2 * j + 2 * i - 3;
      node(6) = node(5) + 1;
      node(7) = node(5) + 2;
%
%  The density RHO can also be set to a random positive value.
%
      for krow = 1 : 8
        for kcol = 1 : 8

          if ( node(krow) < 1 || n < node(krow) )
            fprintf ( 1, '\n' );
            fprintf ( 1, 'WATHEN - Fatal error!\n' );
            fprintf ( 1, '  Index NODE(KROW) out of bounds.\n' );
            fprintf ( 1, '  I = %d\n', i );
            fprintf ( 1, '  J = %d\n', j );
            fprintf ( 1, '  KROW = %d\n', krow );
            fprintf ( 1, '  NODE(KROW) = %d\n', node(krow) );
            error ( 'WATHEN - Fatal error!' );
          elseif ( node(kcol) < 1 || n < node(kcol) )
            fprintf ( 1, '\n' );
            fprintf ( 1, 'WATHEN - Fatal error!\n' );
            fprintf ( 1, '  Index NODE(KCOL) out of bounds.\n' );
            fprintf ( 1, '  I = %d\n', i );
            fprintf ( 1, '  J = %d\n', j );
            fprintf ( 1, '  KCOL = %d\n', kcol );
            fprintf ( 1, '  NODE(KCOL) = %d\n', node(kcol) );
            error ( 'WATHEN - Fatal error!' );
          end

          rho = 1.0;

          a(node(krow),node(kcol)) = a(node(krow),node(kcol)) ...
            + 20.0 * rho * em(krow,kcol) / 9.0;

        end
      end

    end
  end

  return
end


text:
function varargout = process_mne_maxwell( varargin )
% PROCESS_MNE_MAXWELL: MNE-Python call to mne.preprocessing.maxwell_filter: Maxwell filtering / SSS /tSSS
%
% USAGE:   sProcess = process_mne_maxwell('GetDescription')
%            sInput = process_mne_maxwell('Run', sProcess, sInput, method=[])

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2020

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'MNE-Python: maxwell_filter (SSS/tSSS)';
    sProcess.FileTag     = 'sss';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Pre-process';
    sProcess.Index       = 85;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw'};
    sProcess.OutputTypes = {'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.processDim  = [];    % Do not split matrix
    sProcess.Description = 'https://www.nmr.mgh.harvard.edu/mne/stable/generated/mne.preprocessing.maxwell_filter.html';
    % Definition of the options
    % Help
    sProcess.options.help.Comment = 'mne.preprocessing.maxwell_filter
For information about the parameters, click on "Online tutorial"

';
    sProcess.options.help.Type    = 'label';
    % int_order: Order of internal component of spherical expansion
    sProcess.options.int_order.Comment = 'int_order (default=8): ';
    sProcess.options.int_order.Type    = 'value';
    sProcess.options.int_order.Value   = {8,'',0};
    % ext_order: Order of external component of spherical expansion
    sProcess.options.ext_order.Comment = 'ext_order (default=3): ';
    sProcess.options.ext_order.Type    = 'value';
    sProcess.options.ext_order.Value   = {3,'',0};
    % origin: Origin of internal and external multipolar moment space in meters
    sProcess.options.origin.Comment = 'origin (auto or 3D point): ';
    sProcess.options.origin.Type    = 'text';
    sProcess.options.origin.Value   = 'auto';
    % coord_frame: Origin of internal and external multipolar moment space in meters
    sProcess.options.coord_frame.Comment = {'head', 'meg', 'coord_frame (default=head): '; ...
                                            'head', 'meg', ''};
    sProcess.options.coord_frame.Type    = 'radio_linelabel';
    sProcess.options.coord_frame.Value   = 'head';
    % destination: The destination location for the head.
    sProcess.options.destination.Comment = 'destination (empty or 3D-point):';
    sProcess.options.destination.Type    = 'value';
    sProcess.options.destination.Value   = {[], 'list', 3};
    % regularize
    sProcess.options.regularize.Comment = 'regularize (default=on)';
    sProcess.options.regularize.Type    = 'checkbox';
    sProcess.options.regularize.Value   = 1;
    % ignore_ref
    sProcess.options.ignore_ref.Comment = 'ignore_ref (default=off)';
    sProcess.options.ignore_ref.Type    = 'checkbox';
    sProcess.options.ignore_ref.Value   = 0;
    % st_duration: tSSS
    sProcess.options.st_duration.Comment = 'st_duration (0=disable tSSS, default=10): ';
    sProcess.options.st_duration.Type    = 'value';
    sProcess.options.st_duration.Value   = {0,'s',3};
    % st_correlation
    sProcess.options.st_correlation.Comment = 'st_correlation (default=0.98): ';
    sProcess.options.st_correlation.Type    = 'value';
    sProcess.options.st_correlation.Value   = {0.98,'',2};
    % st_fixed
    sProcess.options.st_fixed.Comment = 'st_fixed (default=on)';
    sProcess.options.st_fixed.Type    = 'checkbox';
    sProcess.options.st_fixed.Value   = 1;
    % st_only
    sProcess.options.st_only.Comment = 'st_only (default=off)';
    sProcess.options.st_only.Type    = 'checkbox';
    sProcess.options.st_only.Value   = 0;
    % mag_scale
    sProcess.options.mag_scale.Comment = 'mag_scale (default=100): ';
    sProcess.options.mag_scale.Type    = 'value';
    sProcess.options.mag_scale.Value   = {100,'',4};
    % skip_by_annotation
    sProcess.options.skip_by_annotation.Comment = 'skip_by_annotation: ';
    sProcess.options.skip_by_annotation.Type    = 'text';
    sProcess.options.skip_by_annotation.Value   = 'edge, bad_acq_skip';
    % fine_calibration_file (site specific)
    sProcess.options.calibration.Comment = 'fine-calibration file:';
    sProcess.options.calibration.Type    = 'filename';
    sProcess.options.calibration.Value = {...
        '', ...                               % Filename
        '', ...                               % FileFormat
        'open', ...                           % Dialog type: {open,save}
        'Import fine-calibration file...', ...     % Window title
        'ImportData', ...                     % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                         % Selection mode: {single,multiple}
        'files', ...                          % Selection mode: {files,dirs,files_and_dirs}
        {'.dat', 'fine-calibration file (*.dat)', 'calibration'}, ... % Specify file type
        'DataIn'};                          % DefaultFormats: {ChannelIn,DataIn,DipolesIn,EventsIn,MriIn,NoiseCovIn,ResultsIn,SspIn,SurfaceIn,TimefreqIn
    % cross_talk_file (site specific)
    sProcess.options.ctc.Comment = 'cross-talk file:';
    sProcess.options.ctc.Type    = 'filename';
    sProcess.options.ctc.Value   = {...
        '', ...                               % Filename
        '', ...                               % FileFormat
        'open', ...                           % Dialog type: {open,save}
        'Import cross-talk file...', ...     % Window title
        'ImportData', ...                     % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                         % Selection mode: {single,multiple}
        'files', ...                          % Selection mode: {files,dirs,files_and_dirs}
        {'.fif', 'cross-talk file (*.fif)', 'ctc'}, ... % Specify file type
        'DataIn'};                          % DefaultFormats: {ChannelIn,DataIn,DipolesIn,EventsIn,MriIn,NoiseCovIn,ResultsIn,SspIn,SurfaceIn,TimefreqIn
end

%% ===== RUN =====
function OutputFile = Run(sProcess, sInput) %#ok<DEFNU>
    OutputFile = process_mne_maxwell_py('Run', sProcess, sInput);
end

author: brainstorm-tools


text:
function varargout = asinh(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('increasing','callback');
        operator.derivative = @(x)((1 + x.^2).^-0.5);
        operator.inflection = [-inf 1 0 -1];
        operator.range = [-700 700];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [c, ceq] = ma_OutboundHyperVelVectConstraint(stateLog, eventID, lb, ub, bodyIDApply, celBodyData, maData, constType)
%ma_OutboundHyperVelVectConstraint Summary of this function goes here
%   Detailed explanation goes here
    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);

    bodyEventLog = eventLog(eventLog(:,8)==bodyIDApply,:);
    if(isempty(bodyEventLog))
        finalEntry = eventLog(end,:);
    else
        finalEntry = bodyEventLog(end,:);
    end

    bodyID = finalEntry(8);
    if(bodyID == bodyIDApply || bodyIDApply==-1)
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
        gmu = bodyInfo.gm;
        rVect = finalEntry(2:4)';
        vVect = finalEntry(5:7)';

        [sma, ecc, inc, raan, arg, tru] = getKeplerFromState(rVect,vVect,gmu);
        if(ecc > 1)
            [~, OUnitVector] = computeHyperSVectOVect(sma, ecc, inc, raan, arg, tru, gmu);

            switch constType
                case 'x'
                    vComp = OUnitVector(1);
                case 'y'
                    vComp = OUnitVector(2);
                case 'z'
                    vComp = OUnitVector(3);
                case 'mag'
                    vComp = sqrt(-gmu/sma); %from vis-viva equation; not actually a vector component but wanted to shoe-horn it in here
            end

            if(lb == ub)
                c = [0 0];
                ceq(1) = vComp - ub;
            else
                c(1) = lb - vComp;
                c(2) = vComp - ub;
                ceq = [0]; %#ok<NBRAK>
            end
        else
            c = [0 0];
            ceq = [0]; %#ok<NBRAK>
        end
    else
        c = [0 0];
        ceq = [0]; %#ok<NBRAK>
    end
end


text:
function [ beta ] =upadte_beta( H,W,K,grt)
    Hc=H*W;
    He=abs(Hc).^2;
    %% update beta
    beta=zeros(1,K);
    for k0=1:K
         tmp=Hc(k0,:);
         tmpe=He(k0,:);
         beta(k0)=sqrt(grt(k0))*tmp(k0)/(sum(tmpe)+1);
    end
end

author: guohuayan


text:
function findedgedof3(node,edge,range,varargin)
%% FINDEDGEDOF3 highlights some dofs associated to edges
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

hold on
N = size(node,1);
% set up range
if (nargin<=2) || (isempty(range))
    range = (1:size(edge,1))'+N;
end
if islogical(range)
    range = find(range);
end
if size(range,2)>size(range,1)
    range = range';
end
dofrange = range;
range = range - N;
midEdge = (node(edge(range,1),:)+node(edge(range,2),:))/2;
if length(range)<size(edge,1)
    h = line([node(edge(range,1),1); node(edge(range,2),1)],...
             [node(edge(range,1),2); node(edge(range,2),2)],...
             [node(edge(range,1),3); node(edge(range,2),3)],...
             'LineWidth',1.5,'Color','r');
else
    h = plot3(midEdge(:,1),midEdge(:,2),midEdge(:,3),'r.','MarkerSize', 22);
end
if nargin > 3
    if strcmp(varargin{1},'noindex') || strcmp(varargin{1},'index')
        startidx = 2;
    else
        startidx = 1;
    end
    set(h,varargin{startidx:end});
end
if (nargin <=3) || ~(strcmp(varargin{1},'noindex'))
    text(midEdge(:,1)+0.025,midEdge(:,2)-0.025,midEdge(:,3),int2str(dofrange), ...
         'FontSize',14,'FontWeight','bold','Color','k');
end
%% TODO: help and example


text:
function I = gray2rgb(I)

% converts graylevel image I to identical rgb image
% (useful for compatibility with certain algorithms
% which expect 3 color planes)
%
% Do nothing if image is already rgb
%

if size(I,3) == 1
  I(:,:,3) = I;
  I(:,:,2) = I(:,:,3);
end

author: Cloud-CV

repo: object-proposals


text:
function [RM,RN] = radius(lat)
% radius: calculates meridian and normal radii of curvature.
%
% INPUT
%   lat, Nx1 latitude (rad).
%
% OUTPUT
%   RM, Nx1 meridian radius of curvature (North-South)(m).
%   RN, Nx1 normal radius of curvature (East-West) (m).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%	Titterton, D.H. and Weston, J.L. (2004). Strapdown
% Inertial Navigation Technology (2nd Ed.). Institution
% of Engineering and Technology, USA. Eq. 2.6 and 2.7.
%
% 	Groves, P. (2008). Principles of GNSS, Inertial, and
% Multisensor Integrated Navigation Systems. Artech House, UK.
% Eq. 2.65 and 2.66.
%
%  	R. Gonzalez, J. Giribet, and H. Patiño. An approach to
% benchmarking of loosely coupled low-cost navigation systems,
% Mathematical and Computer Modelling of Dynamical Systems, vol. 21,
% issue 3, pp. 272-287, 2015. Eq. 11.
%
% Version: 006
% Date:    2021/03/19
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

a = 6378137.0;                  % WGS84 Equatorial radius in meters
e = 0.0818191908425;            % WGS84 eccentricity

e2 = e^2;
den = 1 - e2 .* (sin(lat)).^2;

% Meridian radius of curvature: radius of curvature for North-South motion.
RM = a .* (1-e2) ./ (den).^(3/2);

% Normal radius of curvature: radius of curvature for East-West motion.
% AKA transverse radius.
RN = a ./ sqrt(den);

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function [output] = generateRandomSample(model, n)

if (nargin < 1)
    error 'function [output] = generateRandomSample(model, n)';
end
if (nargin < 2)
    n = 5000;
end

m.A = model.S;
m.lb = model.lb;
m.ub = model.ub;
% sample until we have mixedfrac of .6 or less
m = gpSampler(m,10,[],0,0);
m.warmupPts = goodInitialPoint(model, n);
mf = 1;
while (mf > .52)
    [m,mf] = gpSampler(m,[],[],200,300);
    %mf
end
m.warmupPts = m.points;
m = rmfield(m, 'points');
mf = 1;
while (mf > .52)
    [m,mf] = gpSampler(m,[],[],200,300);
    %mf
end
m.warmupPts = m.points;
m = rmfield(m, 'points');
mf = 1;
while (mf > .52)
    [m,mf] = gpSampler(m,[],[],200,300);
    %mf
end

output.point = m.points;
output.mf = mf;

author: opencobra

repo: cobratoolbox


text:
function correct = check_doa_correctness(actual, est, tolerance)
%CHECK_DOA_CORRECTNESS Checks if estimated doa is reasonably close to the
%actual ones. Input must be in ascending order within range of
%(-pi/2, pi/2).
%More precisely, suppose the true DOAs are x_1, x_2, ..., x_k, the estimated
%DOAs must fall within the following regions:
% (-pi/2, x_1 + (x_2 - x_1) / 2),
% (x_1 + (x_2 - x_1) / 2, x_2 + (x_3 - x_2) / 2), ...
% (x_k - (x_k - x_{k-1}) / 2, pi/2)
%If tolerance is not inf, the estimated DOAs must also fall within the following
% regions:
% (x_1 - tolerance, x_1 + tolerance), ...
% (x_k - tolerance, x_k + tolerance)

if nargin == 2
    tolerance = inf;
end
n = length(actual);
if n ~= length(est)
    correct = false;
    return;
end
if n == 1
    correct = abs(actual(1) - est(1)) < tolerance;
elseif n > 1
    spaces = diff(actual);
    max_deviation = min(spaces/2, tolerance);
    correct = true;
    % first
    correct = correct && ...
        (est(1) >= max(-pi/2, actual(1)-tolerance) && ...
        est(1) < actual(1) + max_deviation(1));
    % last
    correct = correct && ...
        (est(end) <= min(pi/2, actual(end)+tolerance) && ...
        est(end) > actual(end) - max_deviation(end));
    % middle
    for ii = 2:(n-1)
        correct = correct && ...
        (est(ii) > actual(ii) - max_deviation(ii-1) && ...
        est(ii) < actual(ii) + max_deviation(ii));
    end
end

end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
% LOGNORMAL_LOGPDF
%
% Y = LOGNORMAL_LOGPDF(X)
% Y = LOGNORMAL_LOGPDF(X, MU, SIGMA)
%
% MU is the location parameter in (-INF, INF), default 0
% SIGMA is the scale parameter in (0, INF)

% Last modified 2010-11-12
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function y = lognormal_logpdf(x, varargin)

if nargin <= 3
  if nargin == 1
    mu = 0;
    sigma = 1;
  else
    mu = varargin{1};
    sigma = varargin{2};
  end
  logx = log(x);
  y = -0.5*log(2*pi) ...
      -logx ...
      -log(sigma) ...
      -0.5*((logx-mu)./sigma).^2;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function wb=writeSheet(wb,sheetName,sheetPosition,captions,units,raw,isIntegers)
% writeSheet
%   Writes a cell matrix to an Excel sheet into using the Java library Apache POI
%
%   workbook        Workbook object representing the Excel file
%   sheetName       name of the sheet
%   sheetPosition   0-based position of the sheet
%   captions        cell array of captions (opt)
%   units           WRITE INFO
%   raw             cell array with the data in the sheet
%   isIntegers      true if numeric values should be integers (opt, default
%                   true)
%
%   Usage: wb=writeSheet(wb,sheetName,sheetPosition,captions,units,raw)

if nargin<7
    isIntegers=true;
end

%Adds the required classes to the static Java path if not already added
addJavaPaths();

%Import required classes from Apache POI.
import org.apache.poi.ss.usermodel.*;
import org.apache.poi.ss.util.*;
import org.apache.poi.hssf.usermodel.*;
import org.apache.poi.xssf.usermodel.*;

%Set default style object and bold style object
defaultStyle = wb.createCellStyle();
if isIntegers
    idx = wb.getCreationHelper().createDataFormat().getFormat('0');
else
    idx = wb.getCreationHelper().createDataFormat().getFormat('##0.00');
end
defaultStyle.setDataFormat(idx);

boldFont=wb.createFont();
boldFont.setBoldweight(boldFont.BOLDWEIGHT_BOLD);
boldStyle=defaultStyle.clone();
boldStyle.setFont(boldFont);

s=wb.createSheet();
wb.setSheetName(sheetPosition, sheetName);

%Create the header row and units
counter=0;
if ~isempty(captions)
    r=s.createRow(counter);
    for i=0:numel(captions)-1
        c=r.createCell(i);
        c.setCellValue(captions{i+1});
        c.setCellStyle(boldStyle);
    end
    counter=counter+1;
end
if ~isempty(units)
    r=s.createRow(counter);
    for i=0:numel(units)-1
        c=r.createCell(i);
        content=units{i+1};
        if iscell(content) && numel(content)==2
            c.setCellValue(content{1});

            %This means that the cell should have a list of allowed values
            dvHelper=XSSFDataValidationHelper(s);

            cellRange=CellRangeAddressList(1,1,i,i); %Units are on the second line
            dvConstraint=dvHelper.createExplicitListConstraint(toValid(content{2}));
            validation=dvHelper.createValidation(dvConstraint, cellRange);
            validation.setShowErrorBox(true);
            s.addValidationData(validation);
        else
            c.setCellValue(content{1});
        end
        c.setCellStyle(boldStyle);
    end
    counter=counter+1;
end

%Loop through and fill in the values
for i=0:size(raw,1)-1
    r=s.createRow(i+counter);

    for j=0:size(raw,2)-1
        if ~isempty(raw{i+1,j+1})
            c=r.createCell(j);
            content=raw{i+1,j+1};
            if iscell(content) && numel(content)==2
                c.setCellValue(content{1});

                %This means that the cell should have a list of allowed
                %values
                dvHelper=XSSFDataValidationHelper(s);
                cellRange=CellRangeAddressList(i+1,i+1,j,j);
                dvConstraint=dvHelper.createExplicitListConstraint(toValid(content{2}));
                validation=dvHelper.createValidation(dvConstraint, cellRange);
                validation.setShowErrorBox(true);
                s.addValidationData(validation);
            else
                c.setCellValue(content);
            end
            c.setCellStyle(defaultStyle);
        end
    end
end

%Pre-determine column widths for each sheet
if strcmp(sheetName,'RXNS')
    widths=[786;2358;7860;15719;3406;7860;3406;3406;2358;3406;7860;7860;3668;7860;7860;4192];
elseif strcmp(sheetName,'METS')
    widths=[786;7860;7860;3668;7860;7860;7860;3406;3668;1834];
elseif strcmp(sheetName,'COMPS')
    widths=[786;3144;7860;3144;7860];
elseif strcmp(sheetName,'GENES')
    widths=[786;3144;7860;3144;3406];
elseif strcmp(sheetName,'MODEL')
    widths=[786;3144;7860;3668;3668;5240;5240;7860;7860;2620;7860];
end;

%Resize columns
for i=0:size(raw,2)-1
    s.setColumnWidth(i,widths(i+1));
end

%Add freeze panes
if counter>0
    s.createFreezePane(0,counter,0,counter);
end
end

function I=toValid(J)
I=cell(numel(J),1);
for i=1:numel(J)
    if ischar(J{i})
        I(i)=J(i);
    else
        if islogical(J{i})
            if J{i}==true
                I{i}='TRUE';
            else
                I{i}='FALSE';
            end
        else
            %Other, most likely numbers
            I{i}=num2str(J{i});
        end
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
classdef TestEstimator
    %TestEstimator

    properties (Constant)
        fields = {'aspect', 'focal', 'ppx', 'ppy', 'R', 't', 'K'};
    end

    methods (Static)
        function test_1
            img1 = imread(fullfile(mexopencv.root(),'test','tsukuba_l.png'));
            img2 = imread(fullfile(mexopencv.root(),'test','tsukuba_r.png'));

            finder = cv.FeaturesFinder('OrbFeaturesFinder');
            features = {finder.find(img1), finder.find(img2)};

            matcher = cv.FeaturesMatcher('BestOf2NearestMatcher');
            %matcher = cv.FeaturesMatcher('AffineBestOf2NearestMatcher');
            m = matcher.match_pairwise(features);

            obj = cv.Estimator('HomographyBasedEstimator');
            %obj = cv.Estimator('AffineBasedEstimator');
            typename = obj.typeid();
            validateattributes(typename, {'char'}, {'row', 'nonempty'});

            [cameras,success] = obj.estimate(features, m);
            validateattributes(success, {'logical'}, {'scalar'});
            if success
                validateattributes(cameras, {'struct'}, {'vector'});
                assert(all(ismember(TestEstimator.fields, fieldnames(cameras))));
                for i=1:numel(cameras)
                    validateattributes(cameras(i).aspect, {'numeric'}, {'scalar'});
                    validateattributes(cameras(i).focal, {'numeric'}, {'scalar'});
                    validateattributes(cameras(i).ppx, {'numeric'}, {'scalar'});
                    validateattributes(cameras(i).ppy, {'numeric'}, {'scalar'});
                    validateattributes(cameras(i).R, {'numeric'}, {'size',[3 3]});
                    validateattributes(cameras(i).t, {'numeric'}, {'vector','numel',3});
                    validateattributes(cameras(i).K, {'numeric'}, {'size',[3 3]});
                end
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [ p, q, pdf ] = normp ( z )

%*****************************************************************************80
%
%% NORMP computes the cumulative density of the standard normal distribution.
%
%  Discussion:
%
%    This is algorithm 5666 from Hart, et al.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by Alan Miller.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
%    Charles Mesztenyi, John Rice, Henry Thacher,
%    Christoph Witzgall,
%    Computer Approximations,
%    Wiley, 1968,
%    LC: QA297.C64.
%
%  Parameters:
%
%    Input, real Z, divides the real ( kind = 8 ) line into two
%    semi-infinite intervals, over each of which the standard normal
%    distribution is to be integrated.
%
%    Output, real P, Q, the integrals of the standard normal
%    distribution over the intervals ( - Infinity, Z] and
%    [Z, + Infinity ), respectively.
%
%    Output, real PDF, the value of the standard normal distribution
%    at Z.
%
  cutoff = 7.071;
  p0 = 220.2068679123761;
  p1 = 221.2135961699311;
  p2 = 112.0792914978709;
  p3 = 33.91286607838300;
  p4 = 6.373962203531650;
  p5 = 0.7003830644436881;
  p6 = 0.03526249659989109;
  q0 = 440.4137358247522;
  q1 = 793.8265125199484;
  q2 = 637.3336333788311;
  q3 = 296.5642487796737;
  q4 = 86.78073220294608;
  q5 = 16.06417757920695;
  q6 = 1.755667163182642;
  q7 = 0.08838834764831844;
  root2pi = 2.506628274631001;

  zabs = abs ( z );
%
%  37 < |Z|.
%
  if ( 37.0  < zabs )

    pdf = 0.0;
    p = 0.0;
%
%  |Z| <= 37.
%
  else

    expntl = exp ( - 0.5  * zabs * zabs );
    pdf = expntl / root2pi;
%
%  |Z| < CUTOFF = 10 / sqrt(2).
%
    if ( zabs < cutoff )

      p = expntl * (((((( ...
          p6   * zabs ...
        + p5 ) * zabs ...
        + p4 ) * zabs ...
        + p3 ) * zabs ...
        + p2 ) * zabs ...
        + p1 ) * zabs ...
        + p0 ) / ((((((( ...
          q7   * zabs ...
        + q6 ) * zabs ...
        + q5 ) * zabs ...
        + q4 ) * zabs ...
        + q3 ) * zabs ...
        + q2 ) * zabs ...
        + q1 ) * zabs ...
        + q0 );
%
%  CUTOFF <= |Z|.
%
    else

      p = pdf / ( ...
        zabs + 1.0  / ( ...
        zabs + 2.0  / ( ...
        zabs + 3.0  / ( ...
        zabs + 4.0  / ( ...
        zabs + 0.65  )))));

    end

  end

  if ( z < 0.0  )
    q = 1.0  - p;
  else
    q = p;
    p = 1.0  - q;
  end

  return
end


text:
function varargout = drawVertices(varargin)
%DRAWVERTICES Draw the vertices of a polygon or polyline.
%
%   drawVertices(POLY)
%   Draws the vertices of the given polygon, using pre-defined style.
%   Default is to draw vertices as squares, with the first vertex filled.
%
%   Example
%     poly = circleToPolygon([20 30 40], 16);
%     drawPolygon(poly);
%     hold on; axis equal;
%     drawVertices(poly);
%
%   See also
%   drawPoint, drawPolygon, drawPolyline

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-12-11, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% extract handle of axis to draw on
if isAxisHandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

var = varargin{1};

%% Manage cell arrays of polygons / polylines

% case of a set of polygons stored in a cell array
if iscell(var)
    N = length(var);
    h = zeros(N, 1);
    for i = 1:N
        state = ishold(gca);
        hold on;
        % check for empty polygons
        if ~isempty(var{i})
            h(i) = drawVertices(ax, var{i}, varargin{2:end});
        end
        if ~state
            hold off
        end
    end

    if nargout > 0
        varargout = {h};
    end

    return;
end

%% Parse coordinates and options

% Extract coordinates of vertices
if size(var, 2) > 1
    % first argument is a vertex array
    px = var(:, 1);
    py = var(:, 2);
    varargin(1) = [];

elseif length(varargin) >= 2 && isnumeric(varargin{1}) && isnumeric(varargin{2})
    % arguments 1 and 2 correspond to x and y coordinate respectively
    px = varargin{1};
    py = varargin{2};
    varargin(1:2) = [];

else
    % unknow input format
    error('Should specify either a N-by-2 array, or 2 N-by-1 vectors');
end

% concatenate input argument(s) with default options
defaults = {'MarkerSize', 6};
if length(varargin) == 1
    varargin = [varargin defaults];
else
    varargin = [{'ks'} defaults varargin];
end

%% Draw the vertices

% draw the vertices
h = plot(ax, px, py, varargin{:});

% draw the first vertex with a different style
plot(ax, px(1), py(1), 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 8);

% format output arg
if nargout > 0
    varargout = {h};
end

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function obj = prtUiSubplot(inputSpec,constructorFun)
% obj = prtUiSubplot(varargin)
%
% obj = prtUiSubplot({{subplot one inputs}, {subplot two inputs},...},{prtGuiManagerAxesType1, prtGuiMangerAxesType2,...});
% obj = prtUiSubplot([2,2,1:2],prtUiManagerAxesType1);
% obj = prtUiSubplot([2,2],prtUiManagerAxesType1);

obj = prtUiManagerMultiAxes;

if nargin < 2 || isempty(constructorFun)
    constructorFun = @prtUiManagerPlot;
end

if iscell(inputSpec)
    if ~iscell(constructorFun)
        constructorFun = repmat({constructorFun},length(inputSpec));
        assert(isa(constructorFun,'function_handle'),'constructorFun must be a function handle');
    else
        assert(length(constructorFun)==length(inputSpec),'length of specified prtGuiManager*s must match the number of specified subplots');
        assert(isa(constructorFun{1},'function_handle'),'constructorFun must be or contain a function handle');
    end

    nAxes = length(inputSpec);
    for iAxes = 1:nAxes
        cInputSpec = inputSpec{iAxes};

        if iscell(cInputSpec)
            assert(length(cInputSpec)==3,'prtUiSubplot inputs must be a cell');
            axesManagers(iAxes) = constructorFun{iAxes}('managedHandle',subplot(cInputSpec{1},cInputSpec{2},cInputSpec{3}, 'parent', obj.managedHandle)); %#ok<AGROW>
        else
            assert(prtUtilIsPositiveInteger(cInputSpec) & numel(cInputSpec)==3,'prtUiSubplot inputs must be a cell or a length(3) vector of positive integers');
            axesManagers(iAxes) = constructorFun{iAxes}('managedHandle',subplot(cInputSpec(1),cInputSpec(2),cInputSpec(3), 'parent', obj.managedHandle)); %#ok<AGROW>
        end
    end
else
    assert(prtUtilIsPositiveInteger(inputSpec) & ismember(numel(inputSpec),[2 3]),'prtUiSubplot inputs must be vector of positive integers');

    if ~iscell(constructorFun)
        assert(length(constructorFun)==1,'length of specified prtGuiManager*s must match the number of specified subplots');
    end

    switch numel(inputSpec)
        case 2
            nAxes = prod(inputSpec);
            if iscell(constructorFun)
                assert(numel(constructorFun)==nAxes,'length of specified prtGuiManager*s must match the number of specified subplots');
            end
            for iAxes = 1:nAxes
                if iscell(constructorFun)
                    cConstructorFun = constructorFun{iAxes};
                else
                    cConstructorFun = constructorFun;
                end

                axesManagers(iAxes) = cConstructorFun('managedHandle',subplot(inputSpec(1),inputSpec(2),iAxes, 'parent', obj.managedHandle)); %#ok<AGROW>
            end

        case 3
            axesManagers = constructorFun('managedHandle', subplot(inputSpec(1), inputSpec(2), inputSpec(3), 'parent', obj.managedHandle));
    end
end

obj.axesManagers = axesManagers;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
clear, clc
        % Definirea variabilelor simbolice
syms ud uq Rs Ls Rr Lr id iD iq iQ M w
syms diddt diqdt diDdt diQdt
        % Variabilele utilizate care nu corespund exact cu simbolurile utilizate
        % - w - omega
        % - diddt - d(id)/dt  si  diqdt - d(iq)/dt
        % - diDdt - d(iD)/dt  si  diQdt - d(iQ)/dt
        % Incarcarea celor patru ecuatii
eq1=-ud+Rs*id+Ls*diddt+M*diddt+M*diDdt;
eq2=-uq+Rs*iq+Ls*diqdt+M*diqdt+M*diQdt;
eq3=Rr*iD+Lr*diDdt+M*diddt+M*diDdt+w*Lr*iQ+w*M*(iq+iQ);
eq4=Rr*iQ+Lr*diQdt+M*diqdt+M*diQdt-w*Lr*iD-w*M*(id+iD);
        % Rezolvarea sistemului de ecuatii
        % avand ca variabile independente cele 4 derivate ale curentilor
        % Rezultatele obtinute sunt intocmai derivatele curentilor
SOL=solve(eq1,eq2,eq3,eq4,'diddt,diqdt,diDdt,diQdt');
        % Extragerea solutiei din structura returnata
diddt=SOL.diddt;
diqdt=SOL.diqdt;
diDdt=SOL.diDdt;
diQdt=SOL.diQdt;
        % Afisarea solutiei
disp('did/dt = '),pretty(diddt)
disp('diq/dt = '),pretty(diqdt)
disp('diD/dt = '),pretty(diDdt)
disp('diQ/dt = '),pretty(diQdt)
        % Verificarea solutiei prin inlocuirea rezultatelor in sistemul de ecuatii
Verif1=simple(-ud+Rs*id+Ls*diddt+M*diddt+M*diDdt);
Verif2=simple(-uq+Rs*iq+Ls*diqdt+M*diqdt+M*diQdt);
Verif3=simple(Rr*iD+Lr*diDdt+M*diddt+M*diDdt+w*Lr*iQ+w*M*(iq+iQ));
Verif4=simple(Rr*iQ+Lr*diQdt+M*diqdt+M*diQdt-w*Lr*iD-w*M*(id+iD));
        % Formarea unei vectori cu erorile obtinute
        % (daca toate elementele sale sunt nule, solutia este buna)
Eroare=[Verif1,Verif2,Verif3,Verif4]

author: Sable

repo: mcbench-benchmarks