text:
classdef ConicalSensor < AbstractSensor
    %ConicalSensor Summary of this class goes here
    %   Detailed explanation goes here

    properties
        name(1,:) char = 'Untitled Sensor';

        %Initial sensor properties (can be changed in sensor state later)
        angle(1,1) double = deg2rad(10) %rad
        range(1,1) double = 1000;       %km
        origin AbstractGeometricPoint
        steeringModel AbstractSensorSteeringModel
        initActiveTf(1,1) logical = true;

        %LVD Data
        lvdData LvdData

        %drawing properties
        color(1,1) ColorSpecEnum = ColorSpecEnum.Green;
        alpha(1,1) double = 0.3;
        showMeshEdges(1,1) logical = false;
    end

    properties(Constant)
        typeEnum = SensorEnum.ConicalSensor;
    end

    methods
        function obj = ConicalSensor(name, angle, range, origin, steeringModel, lvdData)
            arguments
                name(1,:) char
                angle(1,1) double {mustBeGreaterThanOrEqual(angle,0)}
                range(1,1) double {mustBeGreaterThan(range, 0)}
                origin(1,1) AbstractGeometricPoint
                steeringModel(1,1) AbstractSensorSteeringModel
                lvdData(1,1) LvdData
            end

            if(angle > pi/2)
                angle = pi/2;
            end

            obj.name = name;
            obj.angle = angle;
            obj.range = range;
            obj.origin = origin;
            obj.steeringModel = steeringModel;
            obj.lvdData = lvdData;
        end

        function [V,F] = getSensorMesh(obj, sensorState, scElem, dcm, inFrame)
            arguments
                obj(1,1) ConicalSensor
                sensorState(1,1) ConicalSensorState
                scElem(1,1) CartesianElementSet
                dcm(3,3) double
                inFrame(1,1) AbstractReferenceFrame
            end

            active = sensorState.getSensorActiveState();
            if(active)
                sensorSteering = sensorState.getSensorSteeringMode();
                time = scElem.time;
                sensorRange = sensorState.getSensorMaxRange();
                sensorAngle = sensorState.getSensorAngle();

%                 theta = linspace(pi/2, pi/2 + sensorAngle, 10);
%                 xPts = sensorRange*cos(theta);
%                 yPts = sensorRange*sin(theta);
%
%                 xPts = [0, 0, xPts];
%                 yPts = [0, sensorRange/2, yPts];
%
%                 v = normVector([xPts(end); yPts(end)]);
%                 xPts = [xPts 0.5*v(1)*sensorRange];
%                 yPts = [yPts 0.5*v(2)*sensorRange];
%
%                 xPts = [xPts, 0];
%                 yPts = [yPts, 0];
%
%                 PV = [xPts(:), yPts(:)];
%
%                 [X,Y,Z] = revolutionSurface(PV, 15, [0,0, 0, 1]);
%                 fvc = surf2patch(X,Y,Z);
%                 fvc.faces = triangulateFaces(fvc.faces);

                rVectSensorOrigin = obj.getOriginInFrame(time, scElem, inFrame);

                MM = sensorSteering.getSensorDcmToInertial(time, scElem, dcm, inFrame);
                r = rotm2axang(MM*roty(90));
                M = makehgtform('translate',rVectSensorOrigin(:)', 'axisrotate',r(1:3),r(4));

%                 fvc = transformPoint3d(fvc, M);
%                 [V, F] = meshVertexClustering(fvc, 1);

                S = [0,0,0, sensorRange];
                sphere = sphereMesh(S, 'nTheta', ceil(5*2*pi/sensorAngle), 'nPhi', 16);
                sPtsRaw = unique(sphere.vertices,'rows');
                sPtsAngs = dang(repmat([0;0;1],1,size(sPtsRaw,1)), sPtsRaw');
                sPts = sPtsRaw(sPtsAngs <= sensorAngle+1E-10, :);

                sPts2 = transformPoint3d(sPts(:,1), sPts(:,2), sPts(:,3), M);

                V = vertcat(rVectSensorOrigin(:)', sPts2);
                F = convhull(V);

                [V, F] = meshVertexClustering(V,F, 0.001);
            else
                V = [];
                F = [];
            end
        end

        function boreDir = getSensorBoresightDirection(~, sensorState, time, scElem, dcm, inFrame)
            boreDir = sensorState.getSensorSteeringMode().getBoresightVector(time, scElem, dcm, inFrame);
        end

        function sensorDcm = getSensorDcmToInertial(obj, sensorState, scElem, dcm, inFrame)
            time = scElem.time;
            sensorDcm = sensorState.getSensorSteeringMode().getSensorDcmToInertial(time, scElem, dcm, inFrame);
        end

        function rVectOrigin = getOriginInFrame(obj, time, scElem, inFrame)
            newCartElem = obj.origin.getPositionAtTime(time, scElem, inFrame);
            rVectOrigin = [newCartElem.rVect];
        end

        function tf = isVehDependent(obj, sensorState)
            tf = obj.origin.isVehDependent() || sensorState.getSensorSteeringMode().isVehDependent();
        end

        function listboxStr = getListboxStr(obj)
            listboxStr = obj.name;
        end

        function color = getMeshColor(obj)
            color = obj.color;
        end

        function alpha = getMeshAlpha(obj)
            alpha = obj.alpha;
        end

        function tf = getDisplayMeshEdges(obj)
            tf = obj.showMeshEdges;
        end

        function name = getName(obj)
            name = obj.name;
        end

        function tf = isInUse(obj, lvdData)
            tf = lvdData.usesSensor(obj);
        end

        function useTf = openEditDialog(obj)
            output = AppDesignerGUIOutput({false});
            lvd_EditConicalSensorGUI_App(obj, obj.lvdData, output);
            useTf = output.output{1};
        end

        function state = getInitialState(obj)
            state = ConicalSensorState(obj, obj.initActiveTf, obj.steeringModel, obj.angle, obj.range);
        end

        function tf = usesGeometricPoint(obj, point)
            tf = obj.origin == point;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function lambda = ill3_eigenvalues ( )

%*****************************************************************************80
%
%% ILL3_EIGENVALUES returns the eigenvalues of the ILL3 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Output, real LAMBDA(3,1), the eigenvalues.
%
  lambda = [ 3.0; 2.0; 1.0 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym log2 (@var{x})
%% Symbolic log base 2 function.
%%
%% Examples:
%% @example
%% @group
%% log2(sym(256))
%%   @result{} ans = (sym) 8
%%
%% syms x
%% log2(x)
%%   @result{} ans = (sym)
%%       log(x)
%%       ──────
%%       log(2)
%% @end group
%% @end example
%%
%% @seealso{@@sym/log, @@sym/log10}
%% @end defmethod

function z = log2(x)

  z = elementwise_op ('lambda x: sp.log(x, 2)', x);

end

%!assert (isequal (log2 (sym (1024)), sym (10)))

%!assert (isequal (log2 (sym ([2 16; 32 1])), sym ([1 4; 5 0])))

%!test
%! % round-trip
%! syms x
%! f = log2 (x);
%! h = function_handle (f);
%! A = h (1.1);
%! B = log2 (1.1);
%! assert (A, B, -5*eps)

repo: octsympy


text:
classdef AbstractHomogenizedTensorPrinter < CompositeResultsPrinter

    properties (Access = protected)
        nstre
    end

    methods (Access = protected)

        function storeFieldsToPrint(obj,d)
            obj.storeMicroProblemsFields(d);
          %  obj.storeRegularizedDensity(d);
        end

        function createPrinters(obj,d)
            obj.printers        = obj.createMicroProblemsPrinters(d);
          %  obj.printers{end+1} = obj.createRegularizedDensityPrinter(d);
        end

        function computeNstre(obj,ndim)
            if ndim == 2
                obj.nstre = 3;
            elseif ndim == 3
                obj.nstre = 6;
            end
        end

    end

    methods (Access = private)

        function p = createMicroProblemsPrinters(obj,d)
            for istre = 1:obj.nstre
                p{istre} = ResultsPrinter.create('ElasticityMicro',d);
                obj.printers{istre} = p;
            end
        end

        function p = createRegularizedDensityPrinter(obj,d)
            p = ResultsPrinter.create('DensityGauss',d);
        end

        function storeRegularizedDensity(obj,d)
            d.fields = d.regDensity;
            obj.printers{obj.nstre+1}.storeFieldsToPrint(d);
        end

    end

    methods (Access = protected, Abstract)
       storeMicroProblemsFields(obj)
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function af = naca5gen(iaf)
%
% "naca5gen" Generates the NACA 5 digit airfoil coordinates with desired no.
% of panels (line elements) on it.
%      Author : Divahar Jayaraman (j.divahar@yahoo.com)
%
% INPUTS-------------------------------------------------------------------
%       iaf.designation = NACA 5 digit designation (eg. '23012') - STRING !
%                 iaf.n = no of panels (line elements) PER SIDE (upper/lower)
% iaf.HalfCosineSpacing = 1 for "half cosine x-spacing"
%                       = 0 to give "uniform x-spacing"
%          iaf.wantFile = 1 for creating airfoil data file (eg. 'naca2412.dat')
%                       = 0 to suppress writing into a file
%       iaf.datFilePath = Path where the data  file has to be created
%                         (eg. 'af_data_folder/naca5digitAF/')
%                         use only forward slash '/' (Just for OS portability)
%
% OUTPUTS------------------------------------------------------------------
% Data:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%       af.x = x cordinate (nx1 array)
%       af.z = z cordinate (nx1 array)
%      af.xU = x cordinate of upper surface (nx1 array)
%      af.zU = z cordinate of upper surface (nx1 array)
%      af.xL = x cordinate of lower surface (nx1 array)
%      af.zL = z cordinate of lower surface (nx1 array)
%      af.xC = x cordinate of camber line (nx1 array)
%      af.zC = z cordinate of camber line (nx1 array)
%    af.name = Name of the airfoil
%  af.header = Airfoil name ; No of panels ; Type of spacing
%              (eg. 'NACA23012 : [50 panels,Uniform x-spacing]')
%
%
% File:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
% First line : Header eg. 'NACA23012 : [50 panels,Half cosine x-spacing]'
% Subsequent lines : (2*iaf.n+1) rows of x and z values
%
% Typical Inputs:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
% iaf.designation='23012';
% iaf.n=56;
% iaf.HalfCosineSpacing=1;
% iaf.wantFile=1;
% iaf.datFilePath='./'; % Current folder
% iaf.is_finiteTE=0;

% % [[Calculating key parameters-----------------------------------------]]
cld=str2num(iaf.designation(1))*(3/2)/10;
p=0.5*str2num(iaf.designation(2:3))/100;
t=str2num(iaf.designation(4:5))/100;

a0= 0.2969;
a1=-0.1260;
a2=-0.3516;
a3= 0.2843;

if iaf.is_finiteTE ==1
    a4=-0.1015; % For finite thick TE
else
    a4=-0.1036;  % For zero thick TE
end

% % [[Giving x-spacing---------------------------------------------------]]
if iaf.HalfCosineSpacing==1
    beta=linspace(0,pi,iaf.n+1)';
    x=(0.5*(1-cos(beta))); % Half cosine based spacing
    iaf.header=['NACA' iaf.designation ' : [' num2str(2*iaf.n) 'panels,Half cosine x-spacing]'];
else
    x=linspace(0,1,iaf.n+1)';
    iaf.header=['NACA' iaf.designation ' : [' num2str(2*iaf.n) 'panels,Uniform x-spacing]'];
end

yt=(t/0.2)*(a0*sqrt(x)+a1*x+a2*x.^2+a3*x.^3+a4*x.^4);

P=[  0.05     0.1     0.15    0.2     0.25  ];
M=[  0.0580   0.1260  0.2025  0.2900  0.3910];
K=[361.4     51.64   15.957   6.643   3.230 ];

m=spline(P,M,p);
k1=spline(M,K,m);

xc1=x(find(x<=p));
xc2=x(find(x>p));
xc=[xc1 ; xc2];

if p==0
    xu=x;
    yu=yt;

    xl=x;
    yl=-yt;

    zc=zeros(size(xc));
else
    yc1=(1/6)*k1*( xc1.^3-3*m*xc1.^2+m^2*(3-m)*xc1 );
    yc2=(1/6)*k1*m^3*(1-xc2);
    zc=(cld/0.3)*[yc1 ; yc2];

    dyc1_dx=(1/6)*k1*( 3*xc1.^2-6*m*xc1+m^2*(3-m) );
    dyc2_dx=repmat((1/6)*k1*m^3,size(xc2));
    dyc_dx=[dyc1_dx ; dyc2_dx];
    theta=atan(dyc_dx);

    xu=x-yt.*sin(theta);
    yu=zc+yt.*cos(theta);

    xl=x+yt.*sin(theta);
    yl=zc-yt.*cos(theta);
end
af.name=['NACA ' iaf.designation];

af.x=[flipud(xu) ; xl(2:end)];
af.z=[flipud(yu) ; yl(2:end)];

indx1=1:min( find(af.x==min(af.x)) );  % Upper surface indices
indx2=min( find(af.x==min(af.x)) ):length(af.x); % Lower surface indices
af.xU=af.x(indx1); % Upper Surface x
af.zU=af.z(indx1); % Upper Surface z
af.xL=af.x(indx2); % Lower Surface x
af.zL=af.z(indx2); % Lower Surface z

af.xC=xc;
af.zC=zc;

lecirFactor=0.8;
af.rLE=0.5*(a0*t/0.2)^2;

le_offs=0.5/100;
dyc_dx_le=(1/6)*k1*( 3*le_offs.^2-6*m*le_offs+m^2*(3-m) );
theta_le=atan(dyc_dx_le);
af.xLEcenter=af.rLE*cos(theta_le);
af.yLEcenter=af.rLE*sin(theta_le);

% % [[Writing iaf data into file------------------------------------------]]
if iaf.wantFile==1
    F1=iaf.header;
    F2=num2str([af.x af.z]);
    F=strvcat(F1,F2);
    fileName=[iaf.datFilePath 'naca' iaf.designation '.dat'];
    dlmwrite(fileName,F,'delimiter','')
end

author: Sable

repo: mcbench-benchmarks


text:
function [out] = idv2idv(n)
% Outputs a transformation matrix for changing from forward to reverse order.
%
% order 1 (Jennie's)
%
% 000, 001, 010, 011, 100, 101, 110, 111
%
% order 2 (mine)
%
% 000, 100, 010, 110, 001, 101, 011, 111
%
% USAGE:
%
%    [out] = idv2idv(n)
%
% INPUT:
%    n:      matrix, size of matrix `(2^n x 2^n)`
%
% OUTPUT:
%    out:    transforamtion matrix

warning('are you sure you want to call this function?');

if n <= 0
    out = 1;
    return;
end

out = sparse(2^n,2^n);
for i = 0:(2^n-1)
    t = dec2bin(i,n);
    t2 = t(end:-1:1);
    i2 = bin2dec(t2);
    out(i+1,i2+1) = 1;
end

author: opencobra

repo: cobratoolbox


text:
function [b,ndx,pos] = uniquesafe(a,flag)
% Ripped unique
% For safety, we use 7.0 (6.1 sort differently...)

nIn = nargin;

if nIn < 1
  error('MATLAB:UNIQUE:NotEnoughInputs', 'Not enough input arguments.');
elseif nIn > 2
  error('MATLAB:UNIQUE:TooManyInputs', 'Too many input arguments.');
end

if nIn == 1
  flag = [];
end

rows = size(a,1);
cols = size(a,2);

rowvec = (rows == 1) & (cols > 1);

numelA = prod(size(a));
nOut = nargout;

if isempty(flag)

  % Handle empty: no elements.

  if (numelA == 0)
    % Predefine b to be of the correct type.
    b = a([]);
    if max(size(a)) > 0
      b = reshape(b,0,1);
      ndx = zeros(0,1);
      pos = zeros(0,1);
    else
      ndx = [];
      pos = [];
    end
    return

  elseif (numelA == 1)
    % Scalar A: return the existing value of A.
    b = a; ndx = 1; pos = 1;
    return

    % General handling.
  else

    % Convert to columns
    a = a(:);

    % Convert to double array for purposes of faster sorting.
    % Additionally, UNIQUE calls DIFF, which requires double input.

    whichclass = class(a);
    isdouble = strcmp(whichclass,'double');

    if ~isdouble
      a = double(a);
    end

    % Sort if unsorted.  Only check this for long lists.

    checksortcut = 1000;

    if numelA <= checksortcut | ~(issorted(a))
      if nOut <= 1
        b = sort(a);
      else
        [b,ndx] = sort(a);
      end
    else
      b = a;
      if nOut > 1
        ndx = (1:numelA)';  % If presorted, indices are 1,2,3,...
      end
    end

    % d indicates the location of non-matching entries.

    db = diff(b);

    % Since DIFF returns NaN in both Inf and NaN cases,
    % use slower method of detection if NaN's detected in DIFF(b).
    % After sort, Infs or NaNs will be at ends of list only.

    if (isnan(db(1)) | isnan(db(numelA-1)))
      d = b(1:numelA-1) ~= b(2:numelA);
    else
      d = db ~= 0;
    end

    d(numelA,1) = 1;    % Final element is always member of unique list.

    b = b(d);         % Create unique list by indexing into sorted list.

    if nOut == 3
      pos = cumsum([1;full(d)]);  % Lists position, starting at 1.
      pos(numelA+1) = [];         % Remove extra element introduced by d.
      pos(ndx) = pos;             % Re-reference POS to indexing of SORT.
    end

    % Create indices if needed.
    if nOut > 1
      ndx = ndx(d);
    end

    % Re-convert to correct output data type using FEVAL.
    if ~isdouble
      b = feval(whichclass,b);
    end
  end

  % If row vector, return as row vector.
  if rowvec
    b = b.';
    if nOut > 1
      ndx = ndx.';
      if nOut > 2
        pos = pos.';
      end
    end
  end

else    % 'rows' case
  if ~strcmpi(flag,'rows')
    error('MATLAB:UNIQUE:UnknownFlag', 'Unknown flag.');
  end

  % Handle empty: no rows.

  if (rows == 0)
    % Predefine b to be of the correct type.
    b = a([]);
    ndx = [];
    pos = [];
    b = reshape(b,0,cols);
    if cols > 0
      ndx = reshape(ndx,0,1);
    end
    return

    % Handle scalar: one row.

  elseif (rows == 1)
    b = a; ndx = 1; pos = 1;
    return
  end

  % General handling.
  % Conversion to double not done: SORTROWS is slower for doubles
  % than other types.

  if nOut > 1
    [b,ndx] = sortrows(a);
  else
    b = sortrows(a);
  end

  % d indicates the location of non-matching entries.

  d = b(1:rows-1,:)~=b(2:rows,:);

  % d = 1 if differences between rows.  d = 0 if the rows are equal.

  d = any(d,2);
  d(rows,1) = 1;      % Final row is always member of unique list.

  b = b(d,:);         % Create unique list by indexing into sorted list.

  % Create position mapping vector using CUMSUM.

  if nOut == 3
        pos = cumsum([1;full(d)]);  % Lists position, starting at 1.
        pos(rows+1) = [];           % Remove extra element introduced by d.
        pos(ndx) = pos;             % Re-reference POS to indexing of SORT.
  end

  % Create indices if needed.
  if nOut > 1
    ndx = ndx(d);
  end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function posinpu()
    %[lat1 lon1 lat2,  lon2] =posinpu
    % Position input of two coordinates to build a crossection
    %
    % Alexander Allmann
    % updated: 20.10.2004
    % j.woessner@sed.ethz.ch

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    %...

    hInpuCoord = figure_w_normalized_uicontrolunits( ...
        'Name','Cross-section Input Coordinates',...
        'tag','InpuCoordX','Position',[260 504 300 300],...
        'NumberTitle','off', ...
        'backingstore','on',...
        'Visible','on');
    set(gca,'Box','off','visible','off')

    txt1 = text(...
        'Position',[0. 0.65 0 ],...
        'FontSize',12,...
        'String','Point 1: ');
    txt2 = text(...
        'Position',[0. 0.45 0 ],...
        'FontSize',12,...
        'String','Point 2: ');

    txt3 = text(...
        'Position',[0.35 0.8 0 ],...
        'FontSize',12,...
        'String','Longitude:');
    txt4 = text(...
        'Position',[0.75 0.8 0 ],...
        'FontSize',12,...
        'String','Latitude:');

    inp1_field=uicontrol('Style','edit',...
        'Position',[.70 .60 .25 .10],...
        'Units','normalized','String',num2str(lat1,5),...
        'callback',@callbackfun_001);

    inp2_field=uicontrol('Style','edit',...
        'Position',[.70 .40 .25 .10],...
        'Units','normalized','String',num2str(lat2,5),...
        'callback',@callbackfun_002);

    inp3_field=uicontrol('Style','edit',...
        'Position',[.40 .60 .25 .10],...
        'Units','normalized','String',num2str(lon1,6),...
        'callback',@callbackfun_003);

    inp4_field=uicontrol('Style','edit',...
        'Position',[.40 .40 .25 .10],...
        'Units','normalized','String',num2str(lon2,6),...
        'callback',@callbackfun_004);
    go_button=uicontrol('Style','Pushbutton',...
        'Position',[.40 .05 .15 .15],...
        'Units','normalized',...
        'callback',@callbackfun_005,...
        'String','Go');

    cancel_button=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .15],...
        'Units','normalized',...
        'callback',@callbackfun_006,...
        'String','Cancel');

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lat1=str2double(inp1_field.String);
        inp1_field.String=num2str(lat1);
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lat2=str2double(inp2_field.String);
        inp2_field.String=num2str(lat2);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lon1=str2double(inp3_field.String);
        inp3_field.String=num2str(lon1);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lon2=str2double(inp4_field.String);
        inp4_field.String=num2str(lon2);
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        [xsecx xsecy,  inde] =mysect(tmp1,tmp2,ZG.primeCatalog.Depth,ZG.xsec_defaults.WidthKm,0,lat1,lon1,lat2,lon2);
        nlammap2;
        close(hInpuCoord);
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        figure(mapl);
    end

end

author: CelsoReyes


text:
function adjust_bvecs(ecclog,bvecsfile,newbvecsfile)
% This code is written by Liang Zhan (zhan.liang@ucla.edu)
% This is to rotate gradient table using eddy_correct output
% Input variables are:
%  1) ecclog --- this is one output from eddy_correct
%  2) bvecsfile --- this is bvec generated from dicom
%  3) newbvecsfile --- this is output name for adjusted bvecs

fid=fopen(ecclog);
mat=[];
while ~feof(fid)
    % skip first three lines
    for i=1:3
      fgetl(fid);
    end
    % read four lines
    for i=1:4
     x=str2num(fgetl(fid))
    mat=[mat
         x];
    end
    % skip one line
    fgetl(fid);
end
fclose(fid);

% read bvecs file
bvecs = load(bvecsfile);
if(size(bvecs,2)==3 && size(bvecs,1)>3)
    bvecs = bvecs';
end

% rotate bvecs
rotbvecs = zeros(size(bvecs));
for i = 1:size(bvecs,2)
    %M = mat((i-1)*4+1:i*4,:);
    %M = mat(((i-1)*3+1):i*4,:);
    %M = M(1:3,1:3);
    M = mat(1:3,1:3);
    % extract rotation matrix
    [u,s,v] = svd(M*M');
    R = inv(u*sqrt(s)*v')*M;

    rotbvecs(:,i) = R*bvecs(:,i);
end

save(newbvecsfile,'rotbvecs','-ascii');
end

author: HennyJie

repo: BrainGB

save_path: github-repos/MATLAB/HennyJie-BrainGB


text:
% SCRIPT TEST FOR THE UR10 ROBOT KINEMATICS

% 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/>.
close all
robot = load_robot('practicals', 'UR10');
adjust_view(robot)

% DESCOMENTE LA L�NEA SIGUIENTE PARA PROBAR UN PUNTO SINGULAR
q0 = [0.1 0.1 0.1 0.1 0.1 0.1]';
q = [0.1 -pi/2 pi/2 pi/4 pi/4 pi/4]';

T = directkinematic(robot, q);

fprintf('\nSimple test: try to reach T')
% Llame a la cinemática inversa
qinv = inversekinematic(robot, T, q0)

T_reach = directkinematic(robot, qinv)
'diff T-Treach'
T-T_reach

repo: ARTE


text:
function make_animation ( frame_total )

%*****************************************************************************80
%
%% MAKE_ANIMATION creates an AVI animation file from a sequence of frame files.
%
%  Discussion:
%
%    Before running this program close all figures and don't interfere
%    with the generated figures until the recording process is done.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAKE_ANIMATION\n' );
  fprintf ( 1, '  Create an AVI animation from a sequence of frame files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here we assume a sequence of frame files was created by\n' )
  fprintf ( 1, '  the MAKE_FRAMES program already.\n' );
%
%  Make sure we have the number of frames.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    frame_total = input ( 'Enter the number of frames to view:  ' );
  end
%
%  Set the total number of frames that we will generate.
%
  numframes = frame_total + 1;
%
%  Control how fast final movie will play; total time
%  of movie in seconds will be about (numframes)/(num_frames_per_second)
%
  fps = 10;
  seconds = numframes / fps;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  This movie will contain %d frames.\n', numframes );
  fprintf ( 1, '  The number of frames per second will be %f\n', fps );
  fprintf ( 1, '  so the movie should take %f seconds to play.\n', seconds );
%
%  The AVIFILE function creates a new Audio Video Interleaved (AVI) file.
%  We specify its name and FPS.
%
  aviobj = avifile ( 'peak_movie.avi', 'fps', fps );
%
%  Retrieve each file.
%
  for frame_id = 0 : frame_total
    filename = sprintf ( 'frame%04d.jpg', frame_id );
    H = imshow ( filename );
    set ( H, 'EraseMode', 'xor' );
    set ( gca, 'nextplot', 'replacechildren' );
    shading interp;
    axis tight square off;
    frame = getframe ( gca );
    aviobj = addframe ( aviobj, frame );
  end
%
%  Tell MATLAB we have completed the movie.
%
  aviobj = close ( aviobj );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAKE_ANIMATION:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function [state] = tapas_mcmc_meta_adaptive_ti(data, model, inference, ...
    state, states, si)
%%
%
% Input
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

if state.nsample > inference.nburnin || mod(si, inference.ndiag)
    return;
end

% Number of subjects and chains
[ns, nc] = size(states{6}.llh{1});

% No real adapation is possible
if nc <= 4
    return;
end

llh = zeros(nc, inference.ndiag);
g = [];
for i = 1:inference.ndiag
    llh(:, i) = sum(states{i}.llh{1}, 1);
    if any(isnan(llh(:, i))) || any(abs(llh(:, i)) == inf);
        g = [g i];
    end
end

% Remove nans if any
llh(:, g) = [];
t = state.T{1}(1, :)';

% Maybe something is not working with the first derivatives.
% Let's try to find a workable schedule
for i = 1:nc - 4
    try
        [~, ~, nt] = tapas_genpath(t, llh, i, 'hermite');
        % It worked
        break
    catch
        % Ok keep going
    end
end

% It is a mess and we could not estimate anything
if i == nc - 4
    return
else
    nt = [t(1:i - 1); nt];
end

% Update the temperature schedule
state.T{1} = repmat(nt', ns, 1);
hold on; plot(mean(llh, 2));
end

author: translationalneuromodeling

repo: tapas


text:
% plots a basic tire model given via the vector PacParam
function plotTireModel(PacParam, alpha_data, F_data)
  figure;
  grid on; hold on;
  if(nargin > 2)
    % if tire model data is available, plot it
    scatter(alpha_data, F_data, 10);
    alphaMax = max(alpha_data);
    alphaMin = min(alpha_data);
  else
    alphaMax = 0.5;
    alphaMin = -0.5;
  end
  alpha_sample = alphaMin:0.002:alphaMax;
  FyF = PacParam(3).*sin(PacParam(2).*atan(PacParam(1).*alpha_sample - PacParam(4).*(PacParam(1).*alpha_sample - atan(PacParam(1).*alpha_sample))));
  scatter(alpha_sample, FyF, 'LineWidth', 2);
  xlabel('Side slip angle in rad');
  ylabel('Tire force in N');
  ylim([-6000, 6000]);
  hold on;
  plot(alpha_sample, PacParam(3).*sin(PacParam(2).*atan(PacParam(1).*alpha_sample)));
  legend('Detailed Model', 'Simple Model for Control');
end

author: TUMFTM


text:
function [u,x]=initialvalues(uflux,a,b,init_f,N);
%
%   Approximates the function init_f on the interval [a b] by a piecewise
%   constant function taking values in the set {uflux}.
%   The output is the location of the discontinuities x and u such that
%   the discontinuity between u(i-1) and u(i) is located at x(i).
%
if nargin<5,
  N=1/(uflux(2)-uflux(1));
end;
x=linspace(a,b,N);
xf=0.5*(x(2:N)+x(1:N-1));
u1=feval(init_f,xf);
A=ones(size(u1))'*uflux;
B=ones(size(uflux))'*u1;
[m ind]=min(abs(A'-B));
h=uflux(ind);
x=x(2:N-1);
d=diff(h);
ind=d~=0;
[i s]=find(ind);
u=[h(s) h(N-1)];
x=x(s);

repo: Aero-matlab


text:
function cnn_update_net_group(net_config, update_fn, update_opts)

assert(~isempty(net_config));

group_num=length(net_config.ref.group_infos);

for g_idx=1:group_num

    one_group_info=net_config.ref.group_infos{g_idx};

    one_group_info=update_fn(one_group_info, update_opts);

    net_config.ref.group_infos{g_idx}=one_group_info;

end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
% Test file for SPINCHEME/STARTMULTISTEP:

function pass = test_startMultistep()

tol = 1e-10;

%% 1D:

% Solve 1D KDV equation with a one-step method (ETDRK4) and with a multistep
% method (PECEC736) which needs to do 5 initialization time-steps:
S = spinop('KDV');
dt = 1e-7;
N = 256;
S.tspan = [0 10*dt];

% Solve with PECEC736:
umulti = spin(S, N, dt, 'plot', 'off', 'scheme', 'pecec736');

% Solve with ETDRK4:
u = spin(S, N, dt, 'plot', 'off', 'scheme', 'etdrk4');

% Compare:
pass(1) = norm(u - umulti, inf)/norm(u, inf) < tol;

%% 2D:

% Solve 2D GS equation with a one-step method (ETDRK4) and with a multistep
% method (PECEC736) which needs to do 5 initialization time-steps:
S = spinop2('GS');
dt = 1e-2;
N = 64;
S.tspan = [0 10*dt];

% Solve with ETDRK4:
u = spin2(S, N, dt, 'plot', 'off', 'scheme', 'etdrk4');

% Solve with PECEC736:
umulti = spin2(S, N, dt, 'plot', 'off', 'scheme', 'pecec736');

% Compare:
pass(2) = norm(u-umulti, inf)/norm(u, inf) < tol;

%% SPHERE:

% Solve the AC equation with a one-step method (LIRK4) and with a multistep
% method (IMEXBDF4) which needs to do 3 initialization time-steps:
S = spinopsphere('AC');
dt = 5e-4;
N = 128;
S.tspan = [0 10*dt];

% Solve with LIRK4:
u = spinsphere(S, N, dt, 'plot', 'off', 'scheme', 'lirk4');

% Solve with IMEXBDF4:
umulti = spinsphere(S, N, dt, 'plot', 'off', 'scheme', 'imexbdf4');

% Compare:
pass(3) = norm(u-umulti, inf)/norm(u, inf) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function exact = p05_exact ( )

%*****************************************************************************80
%
%% P05_EXACT returns the estimated integral for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 1.5822329637296729331;

  return
end


text:
function [M]=unpadeven(M_pad,L_even_dim)

num_dims=ndims(M_pad);
siz=size(M_pad)-L_even_dim;
switch num_dims
    case 2
        M=M_pad(1:1:siz(1),1:1:siz(2));
    case 3
        M=M_pad(1:1:siz(1),1:1:siz(2),1:1:siz(3));
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function pass = test_constructor( )
% Test the spherefun constructor

% Get tolerance:
tol = 2e3*chebfunpref().cheb2Prefs.chebfun2eps;

f = @(x,y,z) x.^2 + y.^2 + z.^2;
f = redefine_function_handle(f);
g = spherefun(f);
pass(1) = ( SampleError(f, g) < tol );

f = @(x,y,z) exp(-cos(pi*(x+y+z)));
f = redefine_function_handle(f);
g = spherefun(f);
pass(2) = ( SampleError(f, g) < tol );

f = @(x,y,z) 1-exp(x);
g = spherefun(f);
f = redefine_function_handle(f);
pass(3) = ( SampleError(f, g) < tol );

f = @(x,y,z) exp(y);
g = spherefun(f);
f = redefine_function_handle(f);
pass(4) = ( SampleError(f, g) < tol );

f = @(x,y,z) exp(z);
g = spherefun(f);
f = redefine_function_handle(f);
pass(5) = ( SampleError(f, g) < tol );

f = @(x,y,z) cos(x.*y);
g = spherefun(f);
f = redefine_function_handle(f);
pass(6) = ( SampleError(f, g) < tol );

f = @(x,y,z) sin(x.*y.*z);
g = spherefun(f);
f = redefine_function_handle(f);
pass(7) = ( SampleError(f, g) < tol );

f = @(x,y,z) sin(x+ y.*z);
f = redefine_function_handle(f);
g = spherefun(f);
pass(8) = ( SampleError(f, g) < tol );

f = @(x,y,z) sin(x+ y.*z) + 1;
f = redefine_function_handle(f);
g = spherefun(f);
pass(9) = ( SampleError(f, g) < tol );

f = @(x,y,z) 0*x;
g = spherefun(f);
pass(10) = ( norm(g, inf) == 0 );

% Test the vectorize flag is working:
f = spherefun(@(x,y,z) cos(z));
g = spherefun(@(x,y,z) cos(z), 'vectorize');
pass(11) = ( norm(f - g) < tol );

f = spherefun(@(x,y,z) x.*y.*z);
g = spherefun(@(x,y,z) x*y*z, 'vectorize');
pass(12) = ( norm(f - g) < tol );

f = spherefun(@(x,y,z) 1);
g = spherefun(@(x,y,z) 1, 'vectorize');
pass(13) = ( norm(f - g) < tol );

% Test construction from samples
f = spherefun(@(x,y,z) 1 + x.*sin(x.*y));
[m,n] = length(f);
F = sample(f, m + mod(m, 2), n);
g = spherefun(F);
pass(14) = ( norm(f - g) < tol );

f = spherefun(@(x,y,z) 1 + 0*x);
F = ones(2, 2);
g = spherefun(F);
pass(15) = ( norm(f - g) < tol );

F = ones(1, 2);
try
    g = spherefun(F);
    pass(16) = false;
catch ME
    pass(16) = strcmp(ME.identifier,'CHEBFUN:SPHEREFUN:constructor:poleSamples');
end

% Test construction from coefficients.
f = spherefun(@(x,y,z) exp(-10*((x-1/sqrt(2)).^2 + (z-1/sqrt(2)).^2 + y.^2)));
C = coeffs2(f);
g = spherefun(C,'coeffs');
pass(17) = ( norm(f - g) < tol );

% Test fixed rank construction
ff = @(x,y,z) exp(-10*((x-1/sqrt(2)).^2 + (z-1/sqrt(2)).^2 + y.^2));
f = spherefun(ff,5);
pass(18) = ( rank(f) == 5 );
f = spherefun(ff,6);
pass(19) = ( rank(f) == 6 );
f = spherefun(ff);
g = spherefun(f,7);
pass(20) = ( rank(g) == 7 );

% Test zero rank construction gives zero.
g = spherefun(f,0);
pass(21) = ( rank(g) == 0 );
pass(22) = norm(g) < tol;

try
    f = spherefun(ff,-1);
    pass(23) = false;
catch ME
    pass(23) = strcmp(ME.identifier,'CHEBFUN:SPHEREFUN:constructor:parseInputs:domain3');
end

% Check the 'eps' flag works.
ff = @(x,y,z) exp(-10*((x-1/sqrt(2)).^2 + (z-1/sqrt(2)).^2 + y.^2));
f = spherefun(ff);
g = spherefun(ff,'eps', 1e-5);
pass(24) = rank(g) < rank(f);
[mf,nf] = length(f);
[mg,ng] = length(g);
pass(25) = ( (mg < mf) && (ng < nf) );

% Construction from a single value
f = spherefun(1);
pass(26) = norm(f-1,inf) == 0;

% Construction from a string with (x,y,z) variables
f = spherefun(@(x,y,z) exp(-10*((x-1/sqrt(2)).^2 + (z-1/sqrt(2)).^2 + y.^2)));
g = spherefun('exp(-10*((x-1/sqrt(2)).^2 + (z-1/sqrt(2)).^2 + y.^2))');
pass(27) = norm(f-g,inf) == 0;

% Construction from a string with (l,t) variables
f = spherefun(@(l,t) cos(l).*sin(t) );
g = spherefun('cos(l).*sin(t)');
pass(28) = norm(f-g,inf) == 0;

try
    f = spherefun('x.*y.*z.*w');
    pass(29) = false;
catch ME
    pass(29) = strcmp(ME.identifier,'CHEBFUN:SPHEREFUN:constructor:str2op:depvars');
end

% Construction from zeros matrix should maintain a zeros coefficient
% matrix.
f = spherefun(zeros(5,4));
[n,m] = length(f);
pass(30) = (m == 5) && (n == 4);

end

function f = redefine_function_handle(f)
% nargin(f) = 2, then we are already on the sphere, if nargin(f) = 3,
% then do change of variables:

if ( nargin(f) == 3 )
    % Wrap f so it can be evaluated in spherical coordinates
    f = @(lam, th) spherefun.sphf2cartf(f, lam, th, 0);
end

end

function sample_error = SampleError(h, g)
m = 6;
n = m;
[x, y] = getPoints(m, n);
[L2, T2] = meshgrid(x, y);
F = h(L2, T2);
approx = fevalm(g, x, y);
sample_error = norm(F(:) - approx(:), inf);
end

function [x, y] = getPoints(m, n)

x = trigpts(2*n, [-pi pi]);
y = linspace(0, pi, m).';

end


text:
function a = full(a)
%FULL         Convert Taylor to full
%

% written  05/21/09     S.M. Rump
%

  a.t = full(a.t);

repo: OpenMAS


text:
%% clear all workspace variables
clear ;

%% this is where all the algorithm parameters are
addpath(genpath('./input_data'));
% this is where all the problems are defined
addpath(genpath('./problemdef'));
% this is where all the legacy rng related stuffs are.
% THIS IS NOT VECTORIZED, SO DO NOT USE, SLOW !!!
addpath(genpath('./rand'));

%% global variables that may be used here
global popsize ;
global nreal ;
global nbin ;
global nbits ;
global nobj ;
global ncon ;
global ngen ;

%% load algorithm parameters
load_input_data('input_data/zdt4.in');
pprint('\nInput data successfully entered, now performing initialization\n\n');

%% for debugging puproses
% global min_realvar ;
% global max_realvar ;
% popsize = 24 ;
% nreal = 3 ;
% ngen = 400 ;
% min_realvar = min_realvar(1:nreal);
% max_realvar = max_realvar(1:nreal);
% global min_binvar ;
% global max_binvar ;
% nbin = 2;
% nbits = [3;3];
% min_binvar = [0;0];
% max_binvar = [5;5];

if(nreal > 0)
    obj_col = nreal + 1 : nreal + nobj ;
elseif(nbin > 0)
    obj_col = sum(nbits) + 1 : sum(nbits) + nobj ;
end

%% this is the objective function that we are going to optimize
obj_func = @zdt4 ;

%% allocate memory for pops
if(nreal > 0)
    child_pop = zeros(popsize, nreal + nobj + ncon + 3);
    mixed_pop = zeros(2 * popsize, nreal + nobj + ncon + 3);
elseif(nbin > 0)
    child_pop = zeros(popsize, sum(nbits) + nobj + ncon + 3);
    mixed_pop = zeros(2 * popsize, sum(nbits) + nobj + ncon + 3);
end

%% you need to warm-up the cache if you need to do profiling
% for k = 1:50000
%     tic(); elapsed = toc();
% end

%% start nsga2
tic;
% initialize population
parent_pop = initialize_pop(90);
pprint('Initialization done, now performing first generation\n\n');
parent_pop = evaluate_pop(parent_pop, obj_func);
parent_pop = assign_rank_and_crowding_distance(parent_pop);

do_plot = false ;
% plot the pareto front
if(do_plot); show_plot(1, parent_pop, false, [1 2 3]); end;

do_save = false ;
% save the current pop
if(do_save); save_pop(1, parent_pop, false); end;

for i = 2:ngen
    fprintf('gen = %d\n', i)
    child_pop = selection(parent_pop, child_pop);
    child_pop = mutation_pop(child_pop);
    child_pop(:,obj_col) = 0;
    child_pop = evaluate_pop(child_pop, obj_func);
    mixed_pop = merge_pop(parent_pop, child_pop);
    parent_pop = fill_nondominated_sort(mixed_pop);

    % plot the current pareto front
    % if(do_plot)
    %   show_plot(i, parent_pop, false, [1 2 3], [], ...
    %       [0.0, 1.0], [0.0, 1.0]);
    % end;
    if(do_plot); show_plot(i, parent_pop, false, [1 2 3]); end;

    % save the current pop
    if(do_save); save_pop(i, parent_pop, false); end;
end
toc;
fprintf('Generations finished, now reporting solutions\n');

if(do_save)
    % save the final pop
    save_pop(i, parent_pop, false, 'final');
    % save the best pop
    save_pop(i, parent_pop, false, 'best');
end
fprintf('Routine successfully exited\n');

author: chudur-budur


text:
function [row,header] = genTablePCP(pcp,name)

assert(length(pcp)==11)
header = sprintf(' Torso & Upper & Lower & Upper & Fore- & Head  & Total %s\n       & Leg   & Leg&  & Arm   & arm   &       &       %s\n','\\','\\');

row = sprintf('%s& %1.1f  & %1.1f  & %1.1f  & %1.1f  & %1.1f  & %1.1f & %1.1f %s\n',name,pcp(10),(pcp(2)+pcp(3))/2,(pcp(1)+pcp(4))/2,(pcp(6)+pcp(7))/2,(pcp(5)+pcp(8))/2,pcp(9),pcp(end),'\\');

fprintf('\n%s %s',blanks(length(name)),header);
fprintf('%s\n',row);
end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
% Proces MSE results into MSE curves: average out and/or limit to a
% subset of temporal indices.
%
% Input:
% - results: a configresults cell produced by kafbox_profiler.
% - inds: array of indices for which to calculate output
%
% Output:
% - MSE_avg_setups: cell containing averaged out MSE curves, structure
% corresponds to "results" cell structure.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

function MSE_avg_setups = kafbox_profiler_msecurves(results,inds)

num_setups = length(results);
MSE_avg_setups = cell(num_setups,1);

for setup_ind = 1:num_setups
    setup_results = results{setup_ind};

    num_configs = length(setup_results);
    MSE_avg_configs = cell(num_configs,1);

    for config_ind = 1:num_configs
        config_results = setup_results{config_ind};

        if isfield(config_results{1},'NMSE')
            N = length(config_results{1}.NMSE); % temporary
        else
            N = length(config_results{1}.MSE);
        end
        MSE = zeros(N,1);

        num_sim = length(config_results);
        for sim_ind = 1:num_sim
            simresults = config_results{sim_ind};

            if isfield(simresults,'NMSE') % temporary
                simresults.MSE = simresults.NMSE;
            end
            if nargin<2
                inds = 1:length(simresults.MSE);
            end

            mm = min(length(simresults.MSE(inds)),N);
            MSE = MSE(1:mm) + simresults.MSE(inds(1:mm))/num_sim;
        end
        MSE_avg_configs{config_ind} = MSE;
    end
    MSE_avg_setups{setup_ind} = MSE_avg_configs;
end

repo: kafbox


text:
function my_check_valid_numeric(feat_data)

try

    tmp_check=any(isnan(feat_data)) | any(isinf(feat_data));
    while numel(tmp_check)>1
        tmp_check=any(tmp_check);
    end

    assert(~gather(tmp_check));

catch err_info
    disp('my_check_valid_numeric failed');
    dbstack;
    keyboard;
end

end


text:
                    function ekin
% ***************************************************************
%                   P r o g r a m   EKIN
% ***************************************************************
%
% PURPOSE:
%   Derive differential equation of motion of a mechanical
%   system with one degree of freedom by means of the theo-
%   rem of Kinetic Energy     dEk/dt = N(t,q,qt)
%   Resolve the equation numerically and plots the graphics
%   of the coordinate, velocity and phase plane.
%   If possible, the program could solve the problem analytically.
%
% INPUT DATA:
%   Ek   - expression of the kinetic energy Ek = Ek(q,qt);
%   N    - power of the forces and moments N = N(t,q,qt);
%   q0   - initial value of the coordinate;
%   qt0  - initial value of the velocity;
%   Tend - upper bound of the integration;
%   eps  - precision of the calculations;
%   np   - number of parameters .
%   P{1}, P{2}, ..., P{np} - names of the parameters (array of cells);
%
%  NOTES:
%   1. The coordinate is designed by the symbol 'q' and velocity by 'qt';
%   2. The physical names of the parameters are assigned to the
%      cells of the array P like this: P{1}='m', P{2}='c',...;
%   3. For analytical solution the values of Tend, eps, np and P are not
%      needed.
%   4. Initial values q0, qt0 must to be entered as strings, even though
%      they represent numbers!
%   5. All the data can be input from file or in interactive mode.
%
%  EXAMPLE of DATA FILE:
%   % Data for problem ...
%     Ek   = '1/2*a*qt^2';
%     N    = '-(k*qt + 9.81*sin(q))*qt';
%     q0   = 'q0';  ( or q0  = '0.12';)
%     qt0  = 'qt0'; ( or qt0 = '7.5';)
%     Tend = 20;
%     eps  = 1.e-8;
%     np   = 2;
%     P{1} = 'a';
%     P{2} = 'k';

% ---------------------------------------------------------
%                DATA INPUT OF THE PROBLEM
% =========================================================

 clear
 disp(' ');
 disp(' How will you input the data ?    ');
 disp('     1. From a data file;         ');
 disp('     2. In interactive mode.      ');
 ans = input(' Number of Your choice : '  );
 flag = 0;
 if ans == 1
    while 1
       disp(' ');
       indat = input(' Input the name of the data file :', 's');
       if exist([cd,'\',indat]) % Search only in current directory
          eval(indat);
          flag = 1; break  % Successful
        else               % Unsuccessful
          disp(' ');
          disp([' File ',indat,' not exist!'])
          disp(' You have to:')
          disp(' 1. Enter another DATA file name, or')
          disp(' 2. Input the DATA interactively !')
          ans2 = input(' Your choice, please: ');
          if ans2 == 2, break , end
       end
    end
 end
if flag == 0
    %Input of the data in-line mode
    Ek   = input(' Expression of the kinetic energy Ek : ','s');
    N    = input(' Power of the forces and moments N : ','s'  );
    q0   = input(' Initial value of the coordinate q0 : '     );
    qt0  = input(' Initial value of the velocity qt0 : '      );
    Tend = input(' Upper bound of the integration Tend : '    );
    eps  = input(' Precision of the calculations eps : '      );
    np   = input(' Number of parameters  np : '               );
    % Asigning names of the parameters
    if np > 0
       disp(' ');
       disp(' Enter names of the parameters:')
       for i = 1:np
           ii = num2str(i);
           P{i} = input([' Name of parameter P',ii,': '],'s');
       end
    end
 end

% ---------------------------------------------------------
%             Differential Equation of Motion
% =========================================================

           syms q qt qtt Dq D2q
           dEkdt = diff(Ek,q)*qt + diff(Ek,qt)*qtt;
           deq = (dEkdt - N)/qt;
           deq = subs(deq,{qt,qtt},{Dq,D2q});
           deq = simplify(deq)

% ---------------------------------------------------------
%                  ANALYTICAL SOLUTION
% =========================================================

 disp(' ');
 ans = input(' Would you like analytical solution? (Y/N): ','s');
 if ans =='Y' | ans == 'y'
    if ~isstr(q0) , q0  = num2str(q0) ; end % Repairing user
    if ~isstr(qt0), qt0 = num2str(qt0); end % input errors!
    inicond = ['q(0)=',q0,',Dq(0)=',qt0];
    q = dsolve(char(deq), inicond, 't');
    if ~isempty(q)
        disp(' ');
        disp('   Low of Motion   ');
        disp(' ***************** ');
        disp(' ');
        disp('q = '); pretty(q)
        disp(' ');
        fname = input(' Name of file to write solution: ','s');
        save(fname, 'q');
    end
    disp(' ');
    ans = input(' Would you like numerical solution? (Y/N): ','s');
    if ans == 'N' | ans == 'n', return, end
    q = 'q'; % Clear analitical solution from q !
 end

% ---------------------------------------------------------
%                  NUMERICAL SOLUTION
% =========================================================

% Input the name of the file-function
disp(' ');
fname = input(' Name of the File-function to be generated: ','s');
flag1 = 'Y';
if exist([cd,'\',fname]) % Search only in curent directory!
    disp(' ');
    disp([' A file-function with name ',fname,' already exist !'])
    flag1 = input(' Overwrite it ? (Y/N): ', 's');
end

% ---------------------------------------------------------
%              GENERATING THE FILE-FUNCTION
% ---------------------------------------------------------

if ( flag1 == 'Y' | flag1 == 'y' )
   qtt = solve(deq,'D2q');
   qtt = subs(qtt,{q,Dq},{'y(1)','y(2)'});
% Opening the file to write file-function
   [Fid,mes] = fopen([fname,'.m'],'wt');
% Generating the string with physical parameters: m, c ...
   strpar = '';
   for j = 1:np
      strpar = [strpar,',',P{j}];
   end
   disp(' ');
   titl = input(' Denomination of the Problem: ','s');
%       Writing the headline of the File-function
   fprintf(Fid,['function yt = ',fname,'(t,y',strpar,')\n']);
   fprintf(Fid,['%% ',titl]);
%       Writing the first derivatives
   fprintf(Fid,'\n%% The first derivatives\n');
   fprintf(Fid, '  yt(1) = y(2); \n');
   fprintf(Fid,['  yt(2) = ',char(qtt),'; \n']);
   fprintf(Fid,'  yt = yt''; \n');
   fprintf(Fid,['%% *** End of File-function ',fname,' ***']);
   fclose(Fid);
   edit(fname)
end

% ---------------------------------------------------------
%        INTEGRATION AND VISUALIZATION OF THE REZULTS
% ---------------------------------------------------------

flag2 = 0;
% Initial entering values of the parameters and generating
% the string with parameters 'P{1}, P{2}, ..., P{np}' to be
% passed to the File-function as actual arguments
if np > 0
   PP = P; % Saving the physical names of the parameters in PP
   parameters = ' ';
   disp(' ');
   disp(' Input the numerical values of the parameters: ')
   for i = 1:np
       i = num2str(i);
       eval(['P{',i,'}=input([''   '',P{',i,'},'' = '']);']);
       parameters = [parameters,',P{',i,'}'];
   end
 else
   parameters = [];
end
% Check-up type of q0 and qt0 and correct
% it if needed
if ischar(q0)
    q0 = str2num(q0);
    if isempty(q0), q0 = input(' q0 = '); end
end
if ischar(qt0)
    qt0 = str2num(qt0);
    if isempty(qt0), qt0 = input(' qt0 = '); end
end
while 1
    if flag2 == 1
        disp(' ');
        eps  = input(' Precision of the computations eps: ');
        Tend = input(' Upper bound of the integration Tend: ');
        q0   = input(' Initial coordinate q0: ');
        qt0  = input(' Initial velocity qt0: ');
        if np > 0
          P = PP; % Restoring the names of the parameters !
          disp(' ');
          disp(' Input the numerical values of the parameters: ')
          for i = 1:np
              i = num2str(i);
              eval(['P{',i,'}=input([''   '',P{',i,'},'' = '']);']);
          end
        end
    end
    y0 = [q0 qt0]; % initial conditions
    options = odeset('AbsTol',eps,'RelTol',100*eps);
    % Choosing of the Solver
    disp('                                        ');
    disp('      Choose the proper Solver:         ');
    disp('  -------------------------------       ');
    disp(' A. Non stiff differential equations    ');
    disp('   1. ode45   - middle precision;       ');
    disp('   2. ode23   - low precision;          ');
    disp('   3. ode113  - from low to upper.      ');
    disp('                                        ');
    disp(' B. Stiff differential equations        ');
    disp('   1. ode15s  - from low to upper;      ');
    disp('   2. ode23s  - low precision;          ');
    disp('   3. ode23t  - middle precision;       ');
    disp('   4. ode23tb - low precision.          ');
    disp('                                        ');
    solver = input(' The name of the Solver: ','s');

    % Integration of the Differential Equations

    eval(['[t,y] = feval(solver,eval([''@'',fname]),',...
                  '[0 Tend],y0,options',parameters,');']);
    % Plotting graphs
    tmin  = min(t);
    tmax  = max(t);
    y1min = min(y(:,1));
    y1max = max(y(:,1));
    y2min = min(y(:,2));
    y2max = max(y(:,2));
    dy1   = y1max - y1min;
    dy2   = y2max - y2min;
    xmin  = y1min - 0.1*dy1;
    xmax  = y1max + 0.1*dy1;
    ymin  = y2min - 0.1*dy2;
    ymax  = y2max + 0.1*dy2;
    % Coordinate q = q(t)
    figure % 1
    comet(t,y(:,1))
    plot(t,y(:,1),[tmin tmax],[0 0],'k'), grid on
    axis([tmin, tmax, xmin, xmax]);
    set(gca,'FontName','Arial Cyr','FontSize',12);
    title('Low of motion {\itq} = {\itq}({\itt})')
    xlabel('{\itt}'); ylabel('{\itq}'); pause
    % Velocity qt = qt(t)
    figure % 2
    comet(t,y(:,2))
    plot(t,y(:,2),[tmin tmax],[0 0],'k'), grid on
    axis([tmin, tmax, ymin, ymax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Velocity {\it qt} = {\it qt}({\itt})')
    xlabel('{\itt}'); ylabel('{\it qt}'); pause
    % Coordinate and Velocity
    figure % 3
    subplot(2,1,1), plot(t,y(:,1),[tmin tmax],[0 0],'k')
    grid on, axis([tmin, tmax, xmin, xmax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Low of motion {\itq} = {\itq}({\itt})')
    subplot(2,1,2), plot(t,y(:,2),[tmin tmax],[0 0],'k')
    grid on, axis([tmin, tmax, ymin, ymax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Velocity {\it qt} = {\it qt}({\itt})')
    pause
    % Phase Plane
    figure % 4
    subplot(1,1,1)
    comet(y(:,1),y(:,2))
    plot(y(:,1),y(:,2), [xmin,xmax],[0 0],'k',...
                  [0 0],[ymin,ymax],'k'), grid on
    axis([xmin, xmax, ymin, ymax]);
    set(gca,'FontName','Arial Cyr','FontSize',12);
    title(' Phase Plane {\it qt} = {\it qt}({\itq})')
    xlabel('{\itq}'), ylabel('{\it qt}'), pause
    flag2 = 1;
    close all
    disp(' ');
    ans = input(' Would you like to continue? (Y/N): ','s');
    if ans == 'n' | ans == 'N', break, end
end

%  **************** End of Program EKIN ******************


text:
function value = r8_factorial ( n )

%*****************************************************************************80
%
%% R8_FACTORIAL returns N!.
%
%  Definition:
%
%    N! = Product ( 1 <= I <= N ) I
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the function.
%    0 <= N.
%
%    Output, real VALUE, the factorial of N.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_FACTORIAL - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'R8_FACTORIAL - Fatal error!' );
  end

  value = 1.0;

  for i = 2 : n
    value = value * i;
  end

  return
end


%% REW_lawn
% %{
imfolder = 'images\REW_lawn';
im_n = 2;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'lawn_01.jpg'];
imfile{2} = [imfolder '\' 'lawn_02.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 1, 'persp', 0, imfolder );
%}


text:
% RGB normal distributions fit to colors around each pixel

function [meanImage, covarMat] = localRGBnormalDistributions(image, windowRadius, epsilon)

    if ~exist('windowRadius', 'var') || isempty(windowRadius)
        windowRadius = 1;
    end
    if ~exist('epsilon', 'var') || isempty(epsilon)
        epsilon = 1e-8;
    end

    [h, w, ~] = size(image);
    N = h * w;
    windowSize = 2 * windowRadius + 1;

    meanImage = imboxfilt(image, windowSize);
    covarMat = zeros(3, 3, N);

    for r = 1 : 3
        for c = r : 3
            temp = imboxfilt(image(:, :, r).*image(:, :, c), windowSize) - meanImage(:,:,r) .*  meanImage(:,:,c);
            covarMat(r, c, :) = temp(:);
        end
    end

    for i = 1 : 3
        covarMat(i, i, :) = covarMat(i, i, :) + epsilon;
    end

    for r = 2 : 3
        for c = 1 : r - 1
            covarMat(r, c, :) = covarMat(c, r, :);
        end
    end

end

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
function [train_sample1, train_sample2, label1, label2] = gen_train_sample_xqda(ID_cam_train, train_feat)

label = ID_cam_train(:, 1);
cam = ID_cam_train(:, 2);   % camera label of all training data

uni_label = unique(label);
ind1 = [];
ind2 = [];
label1 = [];
label2 = [];
for n = 1:length(uni_label)
    pos = find(label == uni_label(n));
    perm = randperm(length(pos));
    tmp1 = pos(perm(1:floor(length(pos)/2)));
    tmp2 = pos(perm(floor(length(pos)/2)+1:floor(length(pos)/2)*2));
    cam1 = cam(tmp1);
    cam2 = cam(tmp2);
    pos2 = find(cam1~=cam2);
    tmp1 = tmp1(pos2);
    tmp2 = tmp2(pos2);
    ind1 = [ind1; tmp1];
    ind2 = [ind2; tmp2];
    label1 = [label1; repmat(uni_label(n), [length(tmp1), 1])];
    label2 = [label2; repmat(uni_label(n), [length(tmp2), 1])];
end
train_sample1 = train_feat(:, ind1)';
train_sample2 = train_feat(:, ind2)';

repo: PRW-baseline


text:
function an=euler2axan(phi,theta,psi,M)
% From Euler angles (phi,theta,psi) to axis angle (x y z alpha)/(gamma,delta):

% v=[(cos(al-gm)+1)*(1-cos(bt));
%    sin(al-gm)*(1-cos(bt));
%    sin(bt)*(sin(al)+sin(gm))]

cost=cos(theta);
sint=cos(theta);
omcost=(1-cost);
if sint==0
    x=0;
    y=0;
    z=1;
else
    x = (cos(phi-psi)+1)*omcost;
    y = sin(phi-psi)*omcost;
    z = sint*(sin(phi)+sin(psi));
end

if (x==0)&&(y==0)&&(z==0)
    x=1;
    y=0;
    z=0;
else
    vn=sqrt(x^2+y^2+z^2);
    x=x/vn;
    y=y/vn;
    z=z/vn;
end

v=[x y z]';

%[THETA,PHI,R] = cart2sph(X,Y,Z)

[gamma,delta,tmp] = cart2sph(x,y,z);

cppp=cos(phi+psi);
Tr=cppp+cppp*cost+cost;
alpha=atan2(-(M(1,2)*v(3)+M(2,3)*v(1)+M(3,1)*v(2)-M(1,3)*v(2)-M(3,2)*v(1)-M(2,1)*v(3)),Tr-1);
% if (-(M(1,2)*v(3)+M(2,3)*v(1)+M(3,1)*v(2)-M(1,3)*v(2)-M(3,2)*v(1)-M(2,1)*v(3)))>=0
%     alpha=acos((Tr-1)/2);
% else
%     alpha=-acos((Tr-1)/2);
% end

an={{gamma,delta,alpha},v};


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Set:
%   Set the specified variable to one.
%   Example: robot = Set(robot, do1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 Set(variable)

global configuration robot

%robot.variable = 1;

%do1 is assigned to the
if variable=='do1'
    T07=directkinematic(robot, robot.q)*robot.tool.TCP;
    robot.tool.Trel=inv(T07)*robot.piece.T0;
end


text:
%% dataNorm
% Below is a demonstration of the features of the |dataNorm| function

%%
clear; close all; clc;

%% Syntax
% |Mn=dataNorm(M,normOpt);|

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


text:
function [ C ] = compute_pid_controller( Td, P )

% COMPUTE_PID_CONTROLLER

[Va_idx, gamma_idx, R_idx]  = Td2idx(Td);
xyu_T = table_T(:,Va_idx, gamma_idx, R_idx);
param_trim_chap5;
[a,TF] = compute_tf_model(xyu_T, P);
Va_T = xyu_T(13);  % Extract trimmed state

% Maximum performance of actuators
P.de_max = deg2rad(20);
P.de_teeth = 60;

P.da_max = deg2rad(20);
P.da_teeth = 60;

P.dr_max = deg2rad(25);
P.dr_teeth = 60;

P.dt_max = 1;
P.dt_freq = 50; % Hz

P.phi_max   =   deg2rad(60);
P.theta_max =   deg2rad(45);

% Altitude state machine parameters
P.altitude_take_off_zone = 40;
P.h_theta_hold_zone      = 10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tuning PID loops
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lateral dynamics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Roll_aileron-hold
ephi_max = deg2rad(15);
P.kp_phi = P.da_max / ephi_max;
ksi_phi = KSI_OPT;
wn_phi = sqrt(P.kp_phi * a.phi2);
P.kd_phi = (2*ksi_phi*wn_phi - a.phi1) / a.phi2;

P.tau_phi = 0.05; % [seconds] <=> 1 autopilot loop

% Course_roll-hold
W_chi_phi = 10; % Bandwith separation factor
                % Between 5 and 10. Safety corresponds to 10.
wn_chi = wn_phi / W_chi_phi; %[rad/sec]
ksi_chi = KSI_OPT;
P.kp_chi = 2*ksi_chi*wn_chi*Va_T/P.gravity;
P.ki_chi = wn_chi^2*Va_T/P.gravity;

% Sideslip_rudder-hold
% TODO

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Longitudinal dynamics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pitch_elevator-hold
etheta_max  = deg2rad(10);
P.kp_theta  = P.de_max / etheta_max * sign(a.theta3);
wn_theta = sqrt(a.theta2 + P.kp_theta*a.theta3);
ksi_theta = KSI_OPT;
if isfinite(a.theta3)
    P.kd_theta  = (2*ksi_theta*wn_theta - a.theta1) / a.theta3;
else
    P.kd_theta  = 0;
    disp("kd_theta is null!");
end
P.tau_theta = 0.1; % [seconds] <=> 10 autopilot loops

% Altitude_pitch-hold
% Scale factor between theta_c and theta
K_theta_DC = P.kp_theta*a.theta3/(a.theta2+P.kp_theta*a.theta3);
W_h_theta = 10; % Between 5 and 15. Safety corresponds to 15.
wn_h = wn_theta / W_h_theta; % [rad/sec]
ksi_h = 1;
P.kp_h  = 2*ksi_h*wn_h/(K_theta_DC*Va_T);
P.ki_h  = wn_h^2 / (K_theta_DC*Va_T);

% Airspeed_thrust-hold
wn_va_dt = 10; % [rad/sec] Need to be tuned.
ksi_va_dt = 1;
P.kp_va_dt = (2*ksi_va_dt*wn_va_dt - a.va1)/a.va2;
P.ki_va_dt = wn_va_dt^2/a.va2;

% Airspeed_pitch-hold
W_va_theta = 7; % [rad/sec] 10 by default, can be lower.
wn_va_theta = wn_theta / W_va_theta;
ksi_va_theta = KSI_OPT;
P.kp_va_theta = (a.va1 - 2*ksi_va_theta*wn_va_theta)/(K_theta_DC*P.gravity);
P.ki_va_theta = -wn_va_theta^2/(K_theta_DC*P.gravity);

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function i4row_sort2_d_test ( )

%*****************************************************************************80
%
%% I4ROW_SORT2_D_TEST tests I4ROW_SORT2_D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 6;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4ROW_SORT2_D_TEST\n' );
  fprintf ( 1, '  For a rectangular integer matrix:\n' );
  fprintf ( 1, '  I4ROW_SORT2_D sorts the elements of the rows.\n' );

  seed = 123456789;

  for i = 1 : m
    for j = 1 : n
      a(i,j) = 10 * i + j;
    end
  end

  i4mat_print ( m, n, a, '  The original matrix:' );

  [ a, seed ] = i4mat_perm2_uniform ( m, n, a, seed );

  i4mat_print ( m, n, a, '  The matrix, permuted by I4MAT_PERM2_UNIFORM:' );

  a = i4row_sort_d ( m, n, a );

  i4mat_print ( m, n, a, '  The row-sorted matrix:' );

  a = i4row_sort2_d ( m, n, a );

  i4mat_print ( m, n, a, '  The element-sorted row-sorted matrix:' );

  return
end

repo: jburkardt-m


text:
seed = 0;
rand('state', seed);
randn('state', seed);

nrows = 3;
ncols = 3;
npixels = nrows*ncols;

% we number pixels in transposed raster scan order (top to bottom, left to right)

% hidden var
HV = reshape(1:npixels, nrows, ncols);
% observed var
OV = reshape(1:npixels, nrows, ncols) + length(HV(:));

% observed factor
OF = reshape(1:npixels, nrows, ncols);
% vertical edge factor VEF(i,j) is the factor for edge HV(i,j) - HV(i+1,j)
VEF = reshape((1:(nrows-1)*ncols), nrows-1, ncols) + length(OF(:));
% horizontal edge factor HEF(i,j) is the factor for edge HV(i,j) - HV(i,j+1)
HEF = reshape((1:nrows*(ncols-1)), nrows, ncols-1) + length(OF(:)) + length(VEF(:));

nvars = length(HV(:))+length(OV(:));
assert(nvars == 2*npixels);
nfac = length(OF(:)) + length(VEF(:)) + length(HEF(:));

K = 2; % number of discrete values for the hidden vars
%O = 1; % each observed pixel is a scalar
O = 2; % each observed pixel is binary

factors = cell(1,3);

% hidden states generate observed 0 or 1 plus noise
%factors{2} = cond_gauss1_kernel(K, O, 'mean', [0 1], 'cov', [0.1 0.1]);
pnoise = 0.2;
factors{1} = tabular_kernel([K O], [1-pnoise pnoise; pnoise 1-pnoise]);
ofactor = 1;

% encourage compatibility between neighboring vertical pixels
factors{2} = tabular_kernel([K K], [0.8 0.2; 0.2 0.8]);
vedge_factor = 2;

%% no constraint between neighboring horizontal pixels
%factors{3} = tabular_kernel([K K], [0.5 0.5; 0.5 0.5]);

factors{3} = tabular_kernel([K K], [0.8 0.2; 0.2 0.8]);
hedge_factor = 3;

factor_ndx = zeros(1, 3);
G = zeros(nvars, nfac);
ns = [K*ones(1,length(HV(:))) O*ones(1,length(OV(:)))];

N = length(ns);
%cnodes = OV(:);
cnodes = [];
dnodes = 1:N;

for i=1:nrows
  for j=1:ncols
    G([HV(i,j), OV(i,j)], OF(i,j)) = 1;
    factor_ndx(OF(i,j)) = ofactor;

    if i < nrows
      G(HV(i:i+1,j), VEF(i,j)) = 1;
      factor_ndx(VEF(i,j)) = vedge_factor;
    end

    if j < ncols
      G(HV(i,j:j+1), HEF(i,j)) = 1;
      factor_ndx(HEF(i,j)) = hedge_factor;
    end

  end
end

fg = mk_fgraph(G, ns, factors, 'discrete', dnodes, 'equiv_class', factor_ndx);

if 1
  % make image with vertical stripes
  I = zeros(nrows, ncols);
  for j=1:2:ncols
    I(:,j) = 1;
  end
else
  % make image with square in middle
  I = zeros(nrows, ncols);
  I(3:6,3:6) = 1;
end

% corrupt image
O = mod(I + (rand(nrows,ncols)> (1-pnoise)), 2);

maximize = 1;
engine = belprop_fg_inf_engine(fg, 'maximize', maximize, 'max_iter', npixels*5);

evidence = cell(1, nvars);
onodes = OV(:);
evidence(onodes) = num2cell(O+1); % values must be in range {1,2}

engine = enter_evidence(engine, evidence);

for i=1:nrows
  for j=1:ncols
    m = marginal_nodes(engine, HV(i,j));
    Ihat(i,j) = argmax(m.T)-1;
  end
end

Ihat

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function field=meg_forward(dip_par,forwpar)
% calculates the magnetic field of n dipoles
% in a realistic volume conductor
% usage: field=meg_forward(dip_par,forwpar)
%
% input:
% dip_par nx6 matrix where each row contains location (first 3 numbers)
%         and moment (second 3 numbers) of a dipole
% forwpar structure containing all information relevant for this
%         calculation; forwpar is calculated with meg_ini
%         You have here an option to include linear transformations in
%         the forward model by specifying forpwar.lintrafo=A
%         where A is an NxM matrix. Then field -> A field
%         You can use that, e.g., if you can write the forward model
%         with M magnetometer-channels plus a matrix multiplication
%         transforming this to a (eventually higher order) gradiometer.
%
% output:
% field  mxn matrix where the i.th column is the field in m channels
%        of the i.th dipole
%
% note:  No assumptions about units are made (i.e. no scaling factors)
%
% Copyright (C) 2003, Guido Nolte
%
% 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$

device_sens=forwpar.device_sens;

field_sens_sphere=getfield_sphere(dip_par,forwpar.device_sens,forwpar.center);
field=field_sens_sphere;clear field_sens_sphere;
if isfield(forwpar,'device_ref')
  field=field-getfield_sphere(dip_par,forwpar.device_ref,forwpar.center);
end
if isfield(forwpar,'device_weights')
  field=field+forwpar.weights*getfield_sphere(dip_par,forwpar.device_weights,forwpar.center);
end

if forwpar.order>0
  coeff=forwpar.coeff_sens;
  if isfield(forwpar,'device_ref')
    coeff=coeff-forwpar.coeff_ref;
  end
  if isfield(forwpar,'device_weights')
    coeff=coeff+forwpar.coeff_weights*forwpar.weights';
  end
  field=field+getfield_corr(dip_par,coeff,forwpar.center,forwpar.order);
end

if isfield(forwpar,'lintrafo');
  field=forwpar.lintrafo*field;
end

return % main function

function field=getfield_sphere(source,device,center)
[ndip,ndum]=size(source);
[nchan,ndum]=size(device);
x1=source(:,1:3)-repmat(center',ndip,1);
n1=source(:,4:6);
x2=device(:,1:3)-repmat(center',nchan,1);
n2=device(:,4:6);
%spherical
bt=leadsphere_all(x1',x2',n2');
n1rep=reshape(repmat(n1',1,nchan),3,ndip,nchan);
b=dotproduct(n1rep,bt);
field=b';
return

function field=getfield_corr(source,coeffs,center,order)
[ndip,ndum]=size(source);
x1=source(:,1:3)-repmat(center',ndip,1);
n1=source(:,4:6);
%correction
if order>0
  scale=10;
  [bas,gradbas]=legs(x1,n1,order,scale);
  nbasis=(order+1)^2-1;
  coeffs=coeffs(1:nbasis,:);
  field=-(gradbas*coeffs)';
end
return

function out=crossproduct(x,y)
% usage: out=testprog(x,y)
% testprog calculates the cross-product of vector x and y
[n,m,k]=size(x);
out=zeros(3,m,k);
out(1,:,:)=x(2,:,:).*y(3,:,:)-x(3,:,:).*y(2,:,:);
out(2,:,:)=x(3,:,:).*y(1,:,:)-x(1,:,:).*y(3,:,:);
out(3,:,:)=x(1,:,:).*y(2,:,:)-x(2,:,:).*y(1,:,:);
return

function out=dotproduct(x,y)
% usage: out=dotproduct(x,y)
% testprog calculates the dotproduct of vector x and y
[n,m,k]=size(x);
outb=x(1,:,:).*y(1,:,:)+x(2,:,:).*y(2,:,:)+x(3,:,:).*y(3,:,:);
out=reshape(outb,m,k);
return

function result=norms(x)
[n,m,k]=size(x);
resultb=sqrt(x(1,:,:).^2+x(2,:,:).^2+x(3,:,:).^2);
result=reshape(resultb,m,k);
return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function xyd_xyn=xydistort_by_xynormalized(xyn, cam)
k1=cam.kc(1);
k2=cam.kc(2);
p1=cam.kc(3);
p2=cam.kc(4);
k3=cam.kc(5);
xn=xyn(1);
yn=xyn(2);
xyd_xyn=zeros(2);
% xyd to xyn
xyd_xyn(1,1)=k2*(xn^2 + yn^2)^2 + k3*(xn^2 + yn^2)^3 + 6*p2*xn + 2*p1*yn + xn*(2*k1*xn ...
+ 4*k2*xn*(xn^2 + yn^2) + 6*k3*xn*(xn^2 + yn^2)^2) + k1*(xn^2 + yn^2) + 1;
% xd to yn
xyd_xyn(1,2)=2*p1*xn + 2*p2*yn + xn*(2*k1*yn + 4*k2*yn*(xn^2 + yn^2) + 6*k3*yn*(xn^2 + yn^2)^2);
% yd to (xn, yn)
xyd_xyn(2,1)=2*p1*xn + 2*p2*yn + yn*(2*k1*xn + 4*k2*xn*(xn^2 + yn^2) + 6*k3*xn*(xn^2 + yn^2)^2);
xyd_xyn(2,2)=k2*(xn^2 + yn^2)^2 + k3*(xn^2 + yn^2)^3 + 2*p2*xn + 6*p1*yn + yn*(2*k1*yn ...
+ 4*k2*yn*(xn^2 + yn^2) + 6*k3*yn*(xn^2 + yn^2)^2) + k1*(xn^2 + yn^2) + 1;

repo: ekfmonoslam


text:
function errSS = dtiComputeCurveDiff(param,origCurve,targetCurve)
% Compute error sum of squares between 2 curves given the warping params
%
%    errSS = dtiComputeCurveDiff(param,origCurve,targetCurve)
%
%    param       - 4x1 vector [xScale xShift yScale yShift]
%    origCurve   - n1x2 matrix [origX origY]
%    targetCurve - n2x2 matrix [targetX targetY]
%
%    errSS - sum of squares of differences between the 2 curves
%
% History:
%    2007/01/17 shc wrote it.
%

if ieNotDefined('param'),       error('Require warping parameters!');  end
if ieNotDefined('origCurve'),   error('Require original curve data!'); end
if ieNotDefined('targetCurve'), error('Require target curve data!');   end

x1 = origCurve(:,1);
x2 = targetCurve(:,1);

startPt = max([min(x1) dtiUnwarpStep(min(x2),param(1:2))]);
endPt   = min([max(x1) dtiUnwarpStep(max(x2),param(1:2))]);

xs = linspace(startPt,endPt,5000);

y1Hat = (spline(origCurve(:,1),origCurve(:,2),xs) + param(4)) * param(3);
y2Hat = spline(targetCurve(:,1),targetCurve(:,2),dtiWarpStep(xs,param));

errSS = sum((y1Hat - y2Hat) .^ 2) / (endPt - startPt);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [output,cost,psave,timing] = solvelower(p,options,lowersolver,xmin,upper,timing)

psave = p;
removeThese = find(p.InequalityConstraintState==inf);
p.F_struc(p.K.f + removeThese,:) = [];
p.K.l = p.K.l - length(removeThese);

removeThese = find(p.EqualityConstraintState==inf);
p.F_struc(removeThese,:) = [];
p.K.f = p.K.f - length(removeThese);

if p.options.bmibnb.cut.bilinear
    p_cut = addBilinearVariableCuts(p);
end
if p.options.bmibnb.cut.evalvariable
    p_cut = addEvalVariableCuts(p_cut);
    psave.evalMap = p_cut.evalMap;
end
if p.options.bmibnb.cut.monomial
    p_cut = addMonomialCuts(p_cut);
end
if p.options.bmibnb.cut.multipliedequality
    p_cut = addMultipliedEqualityCuts(p_cut);
end
if p.options.bmibnb.cut.convexity
    p_cut = addConvexityCuts(p_cut);
end
if p.options.bmibnb.cut.complementarity
    p_cut = addComplementarityCuts(p_cut);
end
% **************************************
% SOLVE NODE PROBLEM
% **************************************
if any(p_cut.ub+1e-8<p_cut.lb)
    output.problem=1;
    cost = inf;
else
    % We are solving relaxed problem (penbmi might be local solver)
    p_cut.monomtable = eye(length(p_cut.c));

    if p.solver.lowersolver.objective.quadratic.convex
        % Setup quadratic
        [p_cut.Q,p_cut.c] = compileQuadratic(p.c,p);

        if nonconvexQuadratic(p_cut.Q);
            p_cut.Q = p.Q;
            p_cut.c = p.c;
        end
    end

    fixed = p_cut.lb >= p_cut.ub;
    if nnz(fixed) == length(p.c)
        % All variables are fixed to a bound
        output.Primal = p.lb;
        res = constraint_residuals(p,output.Primal);
        eq_ok = all(res(1:p.K.f)>=-p.options.bmibnb.eqtol);
        iq_ok = all(res(1+p.K.f:end)>=p.options.bmibnb.pdtol);
        feasible = eq_ok & iq_ok;
        if feasible
            output.problem = 0;
        else
            output.problem = 1;
        end
        cost = output.Primal'*p.Q*output.Primal + p.c'*output.Primal + p.f;
    else

        if nnz(fixed)==0

            if ~isempty(p_cut.bilinears) & 0
                top = size(p_cut.F_struc,1);
                if length(p_cut.K.s)==1 & p_cut.K.s(1)==0
                    p_cut.K.s = [];
                end
                usedterms = zeros(size(p_cut.bilinears,1),1);
                for i = 1:size(p_cut.bilinears,1)
                    if ~usedterms(i)
                        windex = p_cut.bilinears(i,1);
                        xindex = p_cut.bilinears(i,2);
                        yindex = p_cut.bilinears(i,3);
                        if xindex ~=yindex
                            % OK, we have a bilinear term
                            xsquaredindex = find(p_cut.bilinears(:,2)==xindex & p_cut.bilinears(:,3)==xindex);
                            ysquaredindex = find(p_cut.bilinears(:,2)==yindex & p_cut.bilinears(:,3)==yindex);
                            if ~isempty(xsquaredindex) & ~isempty(ysquaredindex)
                                usedterms(i) = 1;
                                usedterms(xsquaredindex) = 1;
                                usedterms(ysquaredindex) = 1;
                                xsquaredindex =  p_cut.bilinears(xsquaredindex,1);
                                ysquaredindex =  p_cut.bilinears(ysquaredindex,1);
                                if 0
                                    Z = zeros(9,size(p_cut.F_struc,2));
                                    Z(1,xsquaredindex+1) = 1;
                                    Z(2,windex+1) = 1;
                                    Z(4,windex+1) = 1;
                                    Z(5,ysquaredindex+1) = 1;
                                    Z(3,xindex+1) = 1;
                                    Z(7,xindex+1) = 1;
                                    Z(6,yindex+1) = 1;
                                    Z(8,yindex+1) = 1;
                                    Z(9,1)=1;
                                else
                                    xL = p.lb(xindex);
                                    yL = p.lb(yindex);

                                    Z = zeros(9,size(p_cut.F_struc,2));
                                    Z(1,xsquaredindex+1) = 1;
                                    Z(2,windex+1) = 1;
                                    Z(4,windex+1) = 1;
                                    Z(5,ysquaredindex+1) = 1;
                                    Z(3,xindex+1) = 1;
                                    Z(7,xindex+1) = 1;
                                    Z(6,yindex+1) = 1;
                                    Z(8,yindex+1) = 1;
                                    Z(9,1)=1;
                                    Z(3,1) = -xL;
                                    Z(7,1) = -xL;
                                    Z(6,1) = -yL;
                                    Z(8,1) = -yL;

                                    Z(1,xindex+1) = -2*xL;
                                    Z(5,yindex+1) = -2*yL;

                                    Z(1,1) = xL^2;
                                    Z(5,1) = yL^2;

                                    Z(4,xindex+1) = -yL;
                                    Z(4,yindex+1) = -xL;
                                    Z(4,1) = xL*yL;

                                    Z(2,xindex+1) = -yL;
                                    Z(2,yindex+1) = -xL;
                                    Z(2,1) = xL*yL;

                                end
                                p_cut.F_struc = [p_cut.F_struc;Z];
                                p_cut.K.s = [p_cut.K.s 3];
                            end
                        end
                    end
                end
            end

            p_cut.linearindicies = 1:length(p.c);
            p_cut.nonlinearindicies = [];
            p_cut.variabletype = zeros(1,length(p.c));
            p_cut.deppattern = eye(length(p.c));
            p_cut.linears = 1:length(p.c);
            p_cut.bilinears = [];
            p_cut.nonlinears = [];
            p_cut.monomials = [];
            p_cut.evaluation_scheme = [];

            tstart = tic;
            output = feval(lowersolver,removenonlinearity(p_cut));
            psave.counter.lowersolved = psave.counter.lowersolved + 1;
            timing.lowersolve = timing.lowersolve + toc(tstart);
            cost = output.Primal'*p_cut.Q*output.Primal + p_cut.c'*output.Primal + p.f;
            % Minor clean-up
            pp=p;
            output.Primal(output.Primal<p.lb) = p.lb(output.Primal<p.lb);
            output.Primal(output.Primal>p.ub) = p.ub(output.Primal>p.ub);
            x=output.Primal;
            return
        else
            pp = p_cut;
            removethese = fixed;
            if ~isempty(p_cut.F_struc)
                p_cut.F_struc(:,1)=p_cut.F_struc(:,1)+p_cut.F_struc(:,1+find(fixed))*p_cut.lb(fixed);
                p_cut.F_struc(:,1+find(fixed))=[];

                rf = find(~any(p_cut.F_struc,2));
                rf = rf(rf<=(p_cut.K.f + p_cut.K.l));
                p_cut.F_struc(rf,:) = [];
                p_cut.K.l = p_cut.K.l - nnz(rf>p_cut.K.f);
                p_cut.K.f = p_cut.K.f - nnz(rf<=p_cut.K.f);
            end
            p_cut.c(removethese)=[];
            if nnz(p_cut.Q)>0
                p_cut.c = p_cut.c + 2*p_cut.Q(find(~removethese),find(removethese))*p_cut.lb(removethese);
                p_cut.Q(:,find(removethese))=[];
                p_cut.Q(find(removethese),:)=[];
            else
                p_cut.Q = spalloc(length(p_cut.c),length(p_cut.c),0);
            end

            if ~isempty(p_cut.binary_variables)
                new_bin = [];
                new_var = find(~fixed);
                for i = 1:length(p_cut.binary_variables)
                    temp = find(p_cut.binary_variables(i) == new_var);
                    new_bin =  [new_bin temp(:)'];
                end
                p_cut.binary_variables = new_bin;
            end
            if ~isempty(p_cut.integer_variables)
                new_bin = [];
                new_var = find(~fixed);
                for i = 1:length(p_cut.integer_variables)
                    temp = find(p_cut.integer_variables(i) == new_var);
                    new_bin =  [new_bin temp(:)'];
                end
                p_cut.integer_variables = new_bin;
            end

            p_cut.lb(removethese)=[];
            p_cut.ub(removethese)=[];
            p_cut.x0(removethese)=[];
            p_cut.monomtable(:,find(removethese))=[];
            p_cut.monomtable(find(removethese),:)=[];

            % The model can become absolutely trivial in some case
            % For instance in ex9_2_2 everything is presolved
            if nnz(p_cut.c)==0 & nnz(p_cut.Q)==0 & size(p_cut.F_struc,1)==0
                % No objective and no constraints
                if all(p_cut.lb <= p_cut.ub)
                    output.Primal = (p_cut.lb + p_cut.ub)/2;
                    output.Primal(isinf(p_cut.lb) & isinf(p_cut.ub))=0;
                    output.problem = 0;
                else
                    output.Primal = zeros(length(p_cut.lb),1);
                    output.problem = 1;
                end
            else
                try
                    tstart = tic;
                    output = feval(lowersolver,removenonlinearity(p_cut));
                    psave.counter.lowersolved = psave.counter.lowersolved + 1;
                    timing.lowersolve = timing.lowersolve + toc(tstart);
                catch
                    1
                end
            end
            x=full(p.c*0);
            x(removethese)=p.lb(removethese);
            x(~removethese)=output.Primal;
            output.Primal = x;
            cost = output.Primal'*pp.Q*output.Primal + pp.c'*output.Primal + p.f;
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function geometry_test0185 ( )

%*****************************************************************************80
%
%% TEST0185 tests CIRCLE_PPPR2IMP_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  p_hi =  10.0;
  p_lo = -10.0;
  test_num = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0185\n' );
  fprintf ( 1, '  CIRCLE_PPPR2IMP_3D is given 3D points P1, P2, P3,\n' );
  fprintf ( 1, '  and a radius R,\n' );
  fprintf ( 1, '  and determines the centers C of two circles\n' );
  fprintf ( 1, '  of the given radius, passing through P1 and P2\n' );
  fprintf ( 1, '  and lying in the plane of P1, P2 and P3.\n' );

  seed = 123456789;

  for test = 1 : test_num

    [ p1, seed ] = r8vec_uniform ( dim_num, p_lo, p_hi, seed );
    [ p2, seed ] = r8vec_uniform ( dim_num, p_lo, p_hi, seed );
    [ p3, seed ] = r8vec_uniform ( dim_num, p_lo, p_hi, seed );

    r_lo = sqrt ( sum ( ( p1(1:dim_num) - p2(1:dim_num) ).^2 ) );
    r_hi = r_lo + 5.0;
    [ r, seed ] = r8_uniform ( r_lo, r_hi, seed );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Radius R = %f\n', r );

    fprintf ( 1, '  Point #1: %f  %f  %f\n', p1(1:dim_num) );
    fprintf ( 1, '  Point #2: %f  %f  %f\n', p2(1:dim_num) );
    fprintf ( 1, '  Point #3: %f  %f  %f\n', p3(1:dim_num) );

    [ pc, normal ] = circle_pppr2imp_3d ( p1, p2, p3, r );

    fprintf ( 1, '  Center #1: %f  %f  %f\n', pc(1:dim_num,1) );
    fprintf ( 1, '  Center #2: %f  %f  %f\n', pc(1:dim_num,2) );
%
%  Check that the points are the right distance from the center.
%
    d11 = sqrt ( sum ( ( p1(1:dim_num) - pc(1:dim_num,1)' ).^2 ) );
    d21 = sqrt ( sum ( ( p2(1:dim_num) - pc(1:dim_num,1)' ).^2 ) );
    d12 = sqrt ( sum ( ( p1(1:dim_num) - pc(1:dim_num,2)' ).^2 ) );
    d22 = sqrt ( sum ( ( p2(1:dim_num) - pc(1:dim_num,2)' ).^2 ) );

    fprintf ( 1, '  %f  %f  %f  %f\n', d11, d21, d12, d22 );
%
%  Check that the radial vector to the point is perpendicular to NORMAL.
%
    d11 = normal(1:dim_num) * ( p1(1:dim_num) - pc(1:dim_num,1)' )';
    d21 = normal(1:dim_num) * ( p2(1:dim_num) - pc(1:dim_num,1)' )';
    d12 = normal(1:dim_num) * ( p1(1:dim_num) - pc(1:dim_num,2)' )';
    d22 = normal(1:dim_num) * ( p2(1:dim_num) - pc(1:dim_num,2)' )';

    fprintf ( 1, '  %f  %f  %f  %f\n', d11, d21, d12, d22 );

  end

  return
end

author: johannesgerer


text:
function varargout = svd( f )
%SVD    Singular value decomposition of a SEPARABLEAPPROX.
%   SVD(F) returns the singular values of F. The number of singular values
%   returned is equal to the rank of the SEPARABLEAPPROX.
%
%   S = SVD(F) returns the singular values of F. S is a vector of singular
%   values in decreasing order.
%
%   [U, S, V] = SVD(F) returns the SVD of F. U and V are quasi-matrices of
%   orthogonal CHEBFUN objects and S is a diagonal matrix with the singular
%   values on the diagonal.
%
%   The length and rank of a SEPARABLEAPPROX are slightly different quantities.
%   LENGTH(F) is the number of pivots used by the constructor, and
%   RANK(F) is the number of significant singular values of F. The relation
%   RANK(F) <= LENGTH(F) should always hold.

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

% Empty check:
if ( isempty( f ) )
    varargout = { chebfun, [], chebfun };
    return
end

% if ( iszero( f ) )
%     varargout = {0};
%     return
% end

% Get the low rank representation for f.
[cols, D, rows] = cdr(f);
d = diag(D);

% Extract information:
dom = f.domain;
width = diff( dom( 1:2 ) );
height = diff( dom( 3:4 ) );

% If the function is the zero function then special care is required.
if ( norm( d ) == 0 )
    if ( nargout > 1 )
        f = 1 + 0*f;
        U = 1/sqrt( width )*simplify(f.cols, [], 'globaltol');
        V = 1/sqrt( height )*simplify(f.rows, [], 'globaltol');
        varargout = { U, 0, V };
    else
        varargout = { 0 };
    end

else

    % If the function is non-zero then do the standard stuff.
    %
    % Algorithm:
    %   f = C D R'                 (cdr decomposition)
    %   C = Q_C R_C                (qr decomposition)
    %   R = Q_R R_R                (qr decomposition)
    %   f = Q_C (R_C D R_R') Q_R'
    %   R_C D R_R' = U S V'        (svd)

    [Qleft, Rleft] = qr( cols );
    [Qright, Rright] = qr( rows );
    [U, S, V] = svd( Rleft * D * Rright.' );
    U = Qleft * U;
    V = Qright * V;

    % Output just like the svd of a matrix.
    if ( nargout > 1 )
        varargout = { U, S, V };
    else
        varargout = { diag( S ) };
    end

end

end


text:
function [seedRow, seedCol] = FindSeedFromUnlabeled( regionMatrix )
%FINDSEEDFROMUNLABELED Summary of this function goes here
%   Detailed explanation goes here

[seedRow, seedCol] = find(regionMatrix == 0, 1);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [lo, hi] = afb(x, af)

% Analysis filter bank
%
% USAGE:
%    [lo, hi] = afb(x, af)
% INPUT:
%    x - N-point vector, where
%            1) N is even
%            2) N >= length(af)
%    af - analysis filters
%    af(:, 1) - lowpass filter (even length)
%    af(:, 2) - highpass filter (even length)
% OUTPUT:
%    lo - Low frequecy output
%    hi - High frequency output
% EXAMPLE:
%    [af, sf] = farras;
%    x = rand(1,64);
%    [lo, hi] = afb(x, af);
%    y = sfb(lo, hi, sf);
%    err = x - y;
%    max(abs(err))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

N = length(x);
L = length(af)/2;
x = cshift(x,-L);

% lowpass filter
lo = upfirdn(x, af(:,1), 1, 2);
lo(1:L) = lo(N/2+[1:L]) + lo(1:L);
lo = lo(1:N/2);

% highpass filter
hi = upfirdn(x, af(:,2), 1, 2);
hi(1:L) = hi(N/2+[1:L]) + hi(1:L);
hi = hi(1:N/2);

author: thomaskuestner


text:
%% tetVolMeanEst
% Below is a demonstration of the features of the |tetVolMeanEst| function
%
%%
clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceAlpha1=0.5;
faceAlpha2=1;
edgeColor=0.25*ones(1,3);
edgeWidth=1.5;
patchColor=[1 0.5 0];

%% Estimating tetrahedral volume based on mean edge length and regular face assumption

% Get a regular tentrahedron
[V,F]=platonic_solid(1,1);

% Calculate true volume
VE=tetVol([1 2 3 4],V)

%Estimated volume for regular tets based on mean edge lengths
[VE_est]=tetVolMeanEst(F,V)

%%
% Plotting model
hf=cFigure;
title('A regular tetrahedron','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

hp=patch('Faces',F,'Vertices',V);
set(hp,'FaceColor',patchColor,'FaceAlpha',faceAlpha1,'lineWidth',edgeWidth,'edgeColor',edgeColor);
camlight headlight;
set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal;  grid on;

%%
% The two metrics coincide in this case but for irregular meshes they may
% divergerge.

%% Using |tetVolMeanEst| to set desired mesh volume for tetgen meshing
% For tetrahedral meshing schemes surface geometry is usually provided. For
% instance triangulated surface data. If the desired element volume can be
% specified then in this case |tetVolMeanEst| can be used to estimate the
% desired element volume given the input surface mesh (provided the surface
% mesh is not remeshed). This is highlighted in the following example.

%%
% Building a geodesic dome surface model
[F,V,~]=geoSphere(2,1);

%%
% Plotting model
hf=cFigure;
title('A triangulated surface model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

hp=patch('Faces',F,'Vertices',V);
set(hp,'FaceColor',patchColor,'FaceAlpha',faceAlpha1,'lineWidth',edgeWidth,'edgeColor',edgeColor);
camlight headlight;
set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal;  grid on;

%%
% The triangles are quite regular and can be used to estimate desired tetrahedral element volume
[regionA]=tetVolMeanEst(F,V); %Volume for regular tets

%%
% Defining input structure

inputStruct.stringOpt='-pq1.2AaYQ';
inputStruct.Faces=F;
inputStruct.Nodes=V;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=ones(size(F,1),1); %Face boundary markers
inputStruct.regionPoints=[0 0 0]; %region points
inputStruct.regionA=regionA;
inputStruct.minRegionMarker=2; %Minimum region marker

%%
% Mesh model using tetrahedral elements using tetGen
[meshOutput]=runTetGen(inputStruct); %Run tetGen

%%
% Access model element and patch data
F=meshOutput.faces;
V=meshOutput.nodes;
C=meshOutput.faceMaterialID;
E=meshOutput.elements;

%%
% PLOTTING MODEL

%Selecting half of the model to see interior
Y=V(:,2); YE=mean(Y(E),2);
L=YE>mean(Y);
[Fs,Cs]=element2patch(E(L,:),C(L));

hf1=cFigure;
subplot(1,2,1);
title('Solid tetrahedral mesh model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;
hps=patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',C,'lineWidth',edgeWidth,'edgeColor',edgeColor);
view(3); axis tight;  axis equal;  grid on;
colormap(autumn);
camlight headlight;
set(gca,'FontSize',fontSize);
subplot(1,2,2);
title('Cut view of Solid tetrahedral mesh model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;
hps=patch('Faces',Fs,'Vertices',V,'FaceColor','flat','CData',Cs,'lineWidth',edgeWidth,'edgeColor',edgeColor);
view(3); axis tight;  axis equal;  grid on;
colormap(autumn);
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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


text:
% Create some data.  First pick out the test dimensions, and the locations
% of some files.
D = 4;
PythonCmd = 'python2.6';
LSHCmd = 'lsh.py';
tmpFile = '/tmp/lshtest.out';
tmpFile2 = '/tmp/lshtest.out2';
subPlots = 1;

%%
% Now run the python command to create the data.
cmd = sprintf('%s %s -d %d -create', PythonCmd, LSHCmd, D);
fprintf('Running the command: %s\n', cmd);
system(cmd);
%%
% Now run the python command to measure the distance data
cmd = sprintf('%s %s -d %d -histogram', PythonCmd, LSHCmd, D);
fprintf('Running the command: %s\n', cmd);
system(cmd);
%%
% Load in the distance data and calculate the distance histograms.
testData = load(sprintf('testData%03d.distances', D));
nBins = 40;
[dnnHist, dnnBins] = hist(testData(:,1), nBins);
[danyHist, danyBins] = hist(testData(:,2), nBins);
if subPlots
    subplot(2,2,1);
end
plot(dnnBins, dnnHist, danyBins, danyHist);
legend('Nearest Neighbor', 'Any Neighbor');
xlabel('Distance')
ylabel('Frequency of Occurance');
title(sprintf('Distance Histogram for %d-D data', D));
%%
% Now calculate the optimum LSH parameters.
N=100000;
deltaTarget = 0.5;
r = 0;
uHash = 1;
uCheck = 1;
results = CalculateMPLSHParameters(N, ...
    dnnHist, dnnBins, danyHist, danyBins, deltaTarget, r, uHash, uCheck);
%%
% Now let's run the W test.
w = results.exactW;
cmd = sprintf('%s %s -d %d -w %g -wTest > %s', PythonCmd, LSHCmd, ...
    D, w, tmpFile);
fprintf('Running the command: %s\n', cmd);
system(cmd);

%%
% And load the W-test results.
cmd = sprintf('egrep -v "[a-zA-Z]" %s | tail -20 > %s', tmpFile, tmpFile2);
[s,res] = system(cmd);
wTest = load(tmpFile2);

%%
% And now create the w-test plot.
if subPlots
    subplot(2,2,2);
end

semilogx(results.wList/results.dScale, results.binNnProb, ...
    wTest(:,1), wTest(:,2), 'bx', ...
    results.wList/results.dScale, results.binAnyProb, ...
    wTest(:,1), wTest(:,3), 'gx', ...
    [results.exactW, results.exactW], [0 1], 'r--');

xlabel('W');
ylabel('Probability of Collision');
legend('p_{nn}', 'p_{nn} experimental', ...
    'p_{any}', 'p_{any} experimental', ...
    'Optimum W', ...
    'Location', 'SouthEast');
title(sprintf('wTest for %d-D data', D));

%%
% Now let's run the K test
cmd = sprintf('%s %s -d %d -w %g -kTest > %s', PythonCmd, LSHCmd, ...
    D, w, tmpFile);
fprintf('Running the command: %s\n', cmd);
system(cmd);

%%
% And grab the K-test results.
cmd = sprintf('grep " 10 " %s > %s', tmpFile, tmpFile2);
system(cmd);
kTest = load(tmpFile2);

%%
% Now plot the k-test results.
if subPlots
    subplot(2,2,3);
end

semilogy(kTest(:,2), kTest(:,4), 'bx', ...
    kTest(:,2), kTest(1,4).^kTest(:,2), 'g-', ...
    kTest(:,2), ...
        results.binNnProb(results.exactBin).^kTest(:,2), 'r--');
xlabel('Number of Projections (K)');
ylabel('Probability');
legend('Experimental Data', 'Extrapolated Prediction', ...
    'Theoretical Prediction', 'Location', 'SouthWest');
title(sprintf('kTest for %d-D data', D));

%%
% Now let's run the L test.
cmd = sprintf('%s %s -d %d -w %g -lTest > %s', PythonCmd, LSHCmd, ...
    D, w, tmpFile);
fprintf('Running the command: %s\n', cmd);
system(cmd);

%%
% And grab the L-test results
cmd = sprintf('grep " 10 " %s > %s', tmpFile, tmpFile2);
system(cmd);
lTest = load(tmpFile2);

%%
% Now plot the results.
if subPlots
    subplot(2,2,4);
end

baseNN = results.binNnProb(results.exactBin).^lTest(1,2);
baseAny = results.binAnyProb(results.exactBin).^lTest(1,2);
semilogy(lTest(:,3), lTest(:,4), 'rx',  ...
    lTest(:,3), lTest(:,5), 'kx', ...
    lTest(:,3), (1-(1-baseNN).^lTest(:,3)), ...
    lTest(:,3), (1-(1-baseAny).^lTest(:,3)), ...
    lTest(:,3), (1-(1-kTest(8,4)).^lTest(:,3)), ...
    lTest(:,3), (1-(1-kTest(8,5)).^lTest(:,3)));

legend('p_{nn} Experimental', 'p_{any} Experimental', ...
    'p_{nn} Theory', 'p_{any} Theory',...
    'p_{nn} k-Prediction', 'p_{any} k-Prediction', ...
    'Location', 'SouthEast');
xlabel('Number of Tables (L)');
ylabel('Probability');
title(sprintf('kTest for %d-D data', D));

%%
save TestLSHCode wTest lTest kTest results D subplots ...
    dnnHist dnnBins danyHist danyBins

%%
% set(gcf,'Position', [100 100 900 800])
pictureFile = sprintf('TestLSHCode-%s.eps', date);
print('-depsc', pictureFile);

author: YahooArchive

repo: Optimal-LSH