text:
%IKINE Inverse kinematics of HAL objects
%
% Executes the inverse kinematics of a HAL (human arm-like) object.
% It returns the two sets of possible solutions for the shoulder and
% wrist in two of the four possible permutations - the other two may be
% created from the others.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% This file modifies file(s) from The Robotics Toolbox for MATLAB (RTB)
% by Peter Corke (www.petercorke.com), see file for statement
%
% Syntax:
%  (1) [q1, q2] = hal.ikine(Th, phi)
%  (2) [q1, q2] = hal.ikine(Th)
%
%  (2) is as per (1) but uses phi := (1)
%
% Outputs:
%  q1 : First family of solutions (mx7 matrix where m = size(Th,3))
%  q2 : Second family of solutions (mx7 matrix where m = size(Th,3))
%
% Inputs:
%  Th  : Hand frame(s), may be a 4x4xm array of m frames
%  phi : Swivel angle,
%        (0) 0     : Uses phi = 0
%        (1)       : Uses the z-axis of Tw for the z-axis of Ts, via
%                     projection
%        (2) 'h2g' : Is the hand-2-goal method, s.t. z-axis of Ts is
%                     the cross product of the swivel axis and g
%        (3) step  : step is a scalar >0 so that the swivel angle
%                      is added as an an extra dimension, with phi =
%                      median human range +/- step increments up to the
%                      human limits. The swivel angle increases with
%                      the increasing index of the extra dimension.
%        (4) PHI   : PHI uses an explicit value, the number of elements
%                     in PHI must be a multiple of the number of hand
%                     points/frames
%
% See also h2fsu HAL.gikine HAL.wikine

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE 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.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 1993-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license,
% Modified 16/6/2014 (HAL is a subclass of SerialLink)

function [q1, q2] = ikine(hal, Th, phi)

if nargin == 2
    [Tf, ~, Tu] = hal.h2fsu(Th);
else
    [Tf, ~, Tu] = hal.h2fsu(Th, phi);
end

[q1g, q2g] = hal.gikine(Tu);

qe = squeeze(pi/2 - acos(dot(Tu(1:3,2,:),Tf(1:3,1,:))));

[q1w, q2w] = wikine(Th, Tf);

q1 = [q1g, qe, q1w];
q2 = [q2g, qe, q2w];

end

author: Allopart

repo: rbpf-gmapping


text:
function cc_write ( prefix, ncc, n, icc, ccc, acc )

%*****************************************************************************80
%
%% CC_WRITE writes a sparse matrix in CC format to 3 files.
%
%  Discussion:
%
%    Three files will be created:
%    * prefix_icc.txt contains NCC ICC values;
%    * prefix_ccc.txt contains N+1 CCC values;
%    * prefix_acc.txt contains NCC ACC values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, a common prefix for the filenames.
%
%    Input, integer NCC, the number of CC elements.
%
%    Input, integer N, the number of columns in the matrix.
%
%    Input, integer ICC(NCC), the CC rows.
%
%    Input, integer CCC(N+1), the compressed CC columns.
%
%    Input, real ACC(NCC), the CC values.
%
  filename_icc = strcat ( prefix, '_icc.txt' );
  i4vec_write ( filename_icc, ncc, icc );

  filename_ccc = strcat ( prefix, '_ccc.txt' );
  i4vec_write ( filename_ccc, n + 1, ccc );

  filename_acc = strcat ( prefix, '_acc.txt' );
  r8vec_write ( filename_acc, ncc, acc );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestCircle
    %TestCircle

    methods (Static)
        function test_circle
            img = 255*ones(64,64,3,'uint8');
            out = cv.circle(img, [32,32], 20);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_options
            img = zeros([50 50 3],'uint8');
            out = cv.circle(img, [25 25], 15, ...
                'Color',[255 0 0], 'Thickness','Filled', 'LineType',8);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.circle();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Infrared small target detection utilizing the multiscale relative local contrast measure
clc;
clearvars;
close all;
for kk = 1
    fold = '.\data\';% 27 images
    try
        img = imread([fold, num2str(kk), '.jpg']);
    catch
        img = imread([fold, num2str(kk), '.bmp']);
    end
    img = img(:,:,1);
    out= RLCM(img);
    figure, imshow(img);
    figure, imshow(out, []);
    figure, surf(out);
end

author: daxjuanxiong

repo: infrared-small-target-detection


text:
%ControlCentor
% when you test the code,please correct the directory in next command
load('G:\K2\Sample.mat');

% Sample is a variable that saves our training database.
LGObj = ConstructLGObj( Sample); % construct an object

Order = [3 4 1 2 5 8 7 10 9 6]; % Order is the ordering of the input in K2 algorithm

u = 4; % u is the maximum edges of node in output graph.

[ DAG,K2Score ] = k2( LGObj,Order,u )
h = view(biograph( DAG ))

author: Sable

repo: mcbench-benchmarks


text:
function testFMinSearchNew
% Test function to show two things:
%
% - how the newly-modified fminsearchbnd works
% - how to use output functions and plot functions in fminsearch or fminsearchbnd (or other Matlab
% optimization routines) - heretofore documentation and examples for these have been sparse.
%
% note: use the fminsearchbnd that was modified 2007-Nov-29 by Ken Purchase, which handles plot and
% output functions properly.
%

    % set up the function that returns a sim structure (the simulation function)
    optFn = @(x, varargin) 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;

    x0 = [-1.2 1];
    mins = [-2 0]; %[2 -inf];
    maxs = [Inf Inf]; %[inf 3];

    % I'm going to pass in a useless extra parameter, just to show that fminsearchbnd and fminsearch
    % pass this back thru to your optimization function.  Your function could make use if this if
    % you watned.
    extraParams = 1;

    % Set up optimization options - you can leave any of these blank and fminsearch will use
    % defaults.
    searchOptions = struct(...
        'Display','none',...
        'MaxIter','200*numberOfVariables',...
        'MaxFunEvals','200*numberOfVariables',...
        'TolX',1e-6,...
        'TolFun',1e-6, ...
        'FunValCheck','off',...
        'OutputFcn', @firstOutputFunction,...
        'PlotFcns',@firstPlotFunction);
    % NOTE: you could add several output or plot functions by incluing a cell array of function
    % handles, such as {@firstOutputFunction, @secondOutputFunction}

    % Run the optimization:
    [outX,fval,exitflag,output] = fminsearchbnd(optFn, x0, mins, maxs, searchOptions, extraParams);

    % Finally, re-run the best case thru the simulation function and display result:
    outX
    finalValue = optFn(outX, extraParams)

    % Define output and print functions.  These functions are nexted WITHIN the overall routine so they
    % have access to variables in the above code if needed.
    %

    function stop = firstOutputFunction(xOutputfcn, optimValues, state, varargin)
        % create an output function for the fMinSearch
        %
        % inputs:
        % 1) xOutputfcn = the current x values
        % 2) optimValues - structure having:
        %         optimValues.iteration = iter;  % iteration number
        %         optimValues.funccount = numf;  % number of function eval's so far
        %         optimValues.fval = f;          % value of the function at current iter.
        %         optimValues.procedure = how;   % how is fminsearch current method (expand, contract, etc)
        % 3) State = 'iter','init' or 'done'     % where we are in the fminsearch algorithm.
        % 4) varargin is passed thru fminsearch to the user function and can be anything.
        %

        stop = false;

        % NOTE: this makes a bit of a messy display, but shows what you can do with an output
        % function.  You can get much of the same information using the fminsearch input option
        % 'Display', 'iter'
        disp(sprintf('Iteration: %d,  Evals: %d,  Current Min Value: %d', ...
            optimValues.iteration, optimValues.funccount, optimValues.fval));
        disp(['Best x so far: [' sprintf('%g ', xOutputfcn) ']']);

        % you could place plotting code here if you didn't want the automatic figure handling of the
        % plot functions.

        % you can also modify the value of 'stop' here to true if you want fminseach to terminate
        % based on any criteria you'd put here.
    end

    function stop = firstPlotFunction(xOutputfcn, optimValues, state, varargin)
        % create an print function for the fMinSearch
        %
        % NOTE: The plot functions do their own management of the plot and axes - if you want to
        % plot on your own figure or axes, just do the plotting in the output function, and leave
        % the plot function blank.
        %
        % One thing the plot function DOES have it that it installs STOP and PAUSE buttons on the
        % plot that allow you to interrupt the optimization to go in and see what's going on, and
        % then resume, or stop the iteration and still have it exit normally (and report output
        % values, etc).
        %
        % inputs:
        % 1) xOutputfcn = the current x values
        % 2) optimValues - structure having:
        %         optimValues.iteration = iter;  % iteration number
        %         optimValues.funccount = numf;  % number of function eval's so far
        %         optimValues.fval = f;          % value of the function at current iter.
        %         optimValues.procedure = how;   % how is fminsearch current method (expand, contract, etc)
        % 3) State = 'iter','init' or 'done'     % where we are in the fminsearch algorithm.
        % 4) varargin is passed thru fminsearch to the user function and can be anything.
        %

        stop = false;

        hold on;
        % this is fun - it simply plots the optimization variable (inverse figure of merit) as it
        % goes along, so you can see it improving, or stop the iterations if it stagnates.
        rectangle('Position', ...
            [(optimValues.iteration - 0.45) optimValues.fval, 0.9, 0.5*optimValues.fval]);
        set(gca, 'YScale', 'log');

        % when you run this, try pressing the 'stop' or 'pause' buttons on the plot.

        % you can add any code here that you desire.

    end

end % end of test code

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
  function distance = penalty_distance(offsets, sizes)
%|function distance = penalty_distance(offsets, sizes)
%|
%| convert scalar offsets to Euclidean distances to neighbors,
%| with the exception that '0' is mapped to 1 for identity case
%|
%| in
%|	offsets		[MM]
%|	sizes		[1 ndim]
%| out
%|	distance	[MM 1]
%|
%| Copyright 2006-12-6, Jeff Fessler, University of Michigan

if nargin == 1 && streq(offsets, 'test'), penalty_distance_test, return, end
if nargin < 2, ir_usage, end

displace = penalty_displace(offsets(:), sizes);
distance = sqrt(sum(displace.^2, 2));
distance(offsets == 0) = 1; % trick: for identity case

%
% penalty_distance_test()
%
function penalty_distance_test

nx = 100; ny = 80; % 2d
[ix iy] = ndgrid(-2:2, -2:2);
offsets = col(ix + iy * nx);
dis = penalty_distance(offsets, [nx ny]);
tru = sqrt(ix(:).^2 + iy(:).^2);
tru(ix == 0 & iy == 0) = 1; % trick: to match identity case
jf_equal(dis, tru)

nx = 10; ny = 8; nz = 7; % 3d
[ix iy iz] = ndgrid(-2:2, -2:2, -2:2);
offsets = col(ix + iy * nx + iz * nx * ny);
dis = penalty_distance(offsets, [nx ny nz]);
tru = sqrt(ix(:).^2 + iy(:).^2 + iz(:).^2);
tru(ix == 0 & iy == 0 & iz == 0) = 1; % trick: to match identity case
jf_equal(dis, tru)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ mean_llr ] = get_bpsk_llr_for_capacity( const_capacity )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

mean_llr = 0 * const_capacity;

load('bpsk_cap.mat');

for i = 1:size(const_capacity,1)
    for j = 1:size(const_capacity,2)
        for k = 1:size(const_capacity,3)
            for snr_index = 1 : length(snr_vec_db)
                if capacity(snr_index) >= const_capacity(i, j, k)
                    break;
                end
            end
            mean_llr(i, j, k) = 4 * 10^(snr_vec_db(snr_index)/10);
        end
    end
end
end

author: tavildar

repo: Polar

save_path: github-repos/MATLAB/tavildar-Polar


text:
function value = r8_epsilon ( )

%*****************************************************************************80
%
%% R8_EPSILON returns the R8 roundoff unit.
%
%  Discussion:
%
%    The roundoff unit is a number R which is a power of 2 with the
%    property that, to the precision of the computer's arithmetic,
%      1 < 1 + R
%    but
%      1 = ( 1 + R / 2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, the roundoff unit.
%
  value = eps;

  return
end


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

    properties
        targets(1,:) AbstractSensorTarget

        lvdData LvdData
    end

    methods
        function obj = LvdSensorTargetSet(lvdData)
            obj.lvdData = lvdData;
        end

        function addTarget(obj, target)
            obj.targets(end+1) = target;
        end

        function removeTarget(obj, target)
            obj.targets([obj.targets] == target) = [];
        end

        function [listBoxStr, targets] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.targets))
                listBoxStr{end+1} = obj.targets(i).getListboxStr(); %#ok<AGROW>
            end

            targets = obj.targets;
        end

        function indPoints = getTargetsForInds(obj, inds)
            indPoints = obj.targets(inds);
        end

        function inds = getIndsForTarget(obj, indTarget)
            inds = [];
            for(i=1:length(indTarget))
                for(j=1:length(obj.targets))
                    if(indTarget(i) == obj.targets(j))
                        inds(end+1) = j; %#ok<AGROW>
                        break;
                    end
                end
            end
        end

        function indSensor = getPointAtInd(obj, ind)
            if(ind > 0 && ind <= length(obj.targets))
                indSensor = obj.targets(ind);
            else
                indSensor = AbstractSensor.empty(1,0);
            end
        end

        function tf = usesGroundObj(obj, groundObj)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGroundObj(groundObj);
            end
        end

        function tf = usesGeometricPoint(obj, point)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGeometricPoint(point);
            end
        end

        function tf = usesGeometricVector(obj, vector)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGeometricVector(vector);
            end
        end

        function tf = usesGeometricCoordSys(obj, coordSys)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGeometricCoordSys(coordSys);
            end
        end

        function tf = usesGeometricRefFrame(obj, refFrame)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGeometricRefFrame(refFrame);
            end
        end

        function tf = usesGeometricAngle(obj, angle)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGeometricAngle(angle);
            end
        end

        function tf = usesGeometricPlane(obj, plane)
            tf = false;
            for(i=1:length(obj.targets))
                tf = tf || obj.targets(i).usesGeometricPlane(plane);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [h]=arrows(x,y,w,fac,color)
% function [h]=arrows(x,y,w,fac,color)
%    Draws arrows with their tails at each point corresponding
%	to identical indices in the matrices x,y.  The matrices u
%	and v are the components of the vector to be represented.
%	fac is the scaling factor
%
%  Geometry of arrowheads (choosing HEADA and HEADL):
%    If the arrow is defined by the points A B C B D where A is the base of
%  the arrow, B is the head, and C and D are the corners of the arrowhead, then
%  HEADA is the angle BAC (or BAD), and HEADL is the ratio of distances AC/AB.
%

% revision 3/20/97 to use nans for line breaks
% much more efficient and returns only a single handle

HEADA=10*pi/180; HEADL=.75;
z=x(:)+i*y(:);

if nargin < 5,color='red';end
if nargin < 4,help arrows,end
w=w(:)*fac;
r=w*HEADL; wr1=r*exp(+i*HEADA); wr2=r*exp(-i*HEADA);
wplot=ones(length(z),6);
wplot(:,1)=z;
wplot(:,[2,4])=(z+w)*ones(1,2);
wplot(:,3)=z+wr1; wplot(:,5)=z+wr2;
wplot(:,6)=z*nan;
wplot=wplot.';
wplot=wplot(:);
%z=eps*ones(size(wplot));
%h=line(real(wplot),imag(wplot),z,'color',color);
h=line(real(wplot),imag(wplot),'color',color);
set(h(1),'userdata',fac);

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function Amatrix=createDFTAmatrix(X,H);
% function Amatrix=createDFTAmatrix(X,H);
% Creates matrix of DFT-amplitudes for a signal X.
% In each column is a spectrum (1:N/2+1), row indices are time-frame indices.
% H is the analysis window. N is the frame length, 50% overlap is used.

% Turn into column vector
X=X(:);
H=H(:);
Nx=length(X);
N=length(H);
M=floor(2*Nx/N-1);
Amatrix=zeros(N/2+1,M);
for k=1:M
    index=(k-1)*N/2+1:(k+1)*N/2;
    F=abs(fft(X(index).*H'));
    Amatrix(:,k)=F(1:N/2+1);
end


text:
function plot_transform(~, eventObj)

hFig = gcf;
idxPC = eventObj.Source.Parent.Parent.UserData.idxPC; % index of PC(s)

% Get new value
prompt = {'UTM Zone:'};
dlgTitle = '';
noLines = 1;
defaultAnswer = {'33n'};
answer = inputdlg(prompt, dlgTitle, noLines, defaultAnswer);
if ~isempty(answer)
    UTMZone = answer{1};
else
    return;
end

% Transform
for i = 1:numel(idxPC)

    if i == 1, logLevelOrig = msg('O', 'GetLogLevel'); end
    msg('O', 'SetLogLevel', 'off');

    % mstruct
    mstruct       = defaultm('utm');
    mstruct.zone  = UTMZone;
    mstruct.geoid = referenceEllipsoid('GRS 80');
    mstruct       = defaultm(mstruct);

    % Transform!
    hFig.UserData.PC{idxPC(i)}.ecef2mapTrafo(mstruct);

    msg('O', 'SetLogLevel', logLevelOrig);

end

% Update plot!
for i = 1:numel(idxPC)
    plot_setPrmAndPlot(hFig, idxPC(i));
end

end

author: pglira


text:
function roomhyp = sample_roomtype1(n, lines, vp, imgwidth, imgheight)
% sample n room hyps
% 3 ways to create hypothesis
% recipe1: 1,left 2,bottom, junc2type,3
% recipe2: 1,left 3,bottom,left, junc2type,4
% recipe3: 2,bottom 3,bottom,left, junc2type,1

lc = [lines.lineclass]; % 0,1,2,3
lr = [lines.leftorright]; % 1 if right, -1 if left, 0 if neither
tb = [lines.above_horizon]; % 1 if above, -1 if below, 0 if neither

oneleft = find((lc==1) & (lr==-1));
twobottom = find((lc==2) & (tb==-1));
threebottomleft = find((lc==3) & (tb==-1) & (lr==-1));

linecode{1} = oneleft;
linecode{2} = twobottom;
linecode{3} = threebottomleft;

% recipe: [linecode1 linecode2 junc2type]
recipe(1,:) = [1 2 3];
recipe(2,:) = [1 3 4];
recipe(3,:) = [2 3 1];

%
% for i = 1:length(linecode),
%     num_linecode(i) = length(linecode{i});
% end
if (length(linecode{1})<1 || length(linecode{2})<1) && ...
   (length(linecode{1})<1 || length(linecode{3})<1) && ...
   (length(linecode{2})<1 || length(linecode{3})<1)
    roomhyp = [];
    return; % fail
end

count = 0;
num_continue_without_progress = 0;
while count < n && num_continue_without_progress < 30
    use_recipe_num = randsample(size(recipe,1), 1); % 1, 2, or 3
    cur_recipe = recipe(use_recipe_num,:);

    l1 = linecode{cur_recipe(1)};
    l2 = linecode{cur_recipe(2)};
    if length(l1)<1 || length(l2)<1
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    ls1 = l1(randsample(length(l1),1));
    ls2 = l2(randsample(length(l2),1));
    [cornerpt degen] = line_intersect(...
                        lines(ls1).point1, lines(ls1).point2, ...
                        lines(ls2).point1, lines(ls2).point2);
    if degen==1
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end
    type = getjunctype(cornerpt, lines(ls1), lines(ls2));
    if type ~= cur_recipe(3)
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    % check if inside img
    MARGIN = 5;
    if ~is_in_image(cornerpt, imgwidth, imgheight, MARGIN)
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    % found a valid way to generate a room hyp
    count = count+1;
    num_continue_without_progress = 0;
    roomhyp(count).corner(4).pt = []; % up to 4 corners
    roomhyp(count).corner(3).pt = cornerpt;
    roomhyp(count).type = 1;

%     global img;
%     disp_vanish(img, lines([ls1 ls2]), vp);
%     plot(cornerpt(1), cornerpt(2), 'x');
%     pause;
%     close;
end

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

author: zouchuhang

repo: LayoutNet


text:
%SETFILES Set files (directory names) for datafile
%
%   A = SETFILES(A,DIRNAMES)
%
% INPUT
%   A         - Datafile
%   DIRNAMES  - Character array with directory (or file) names
%
% OUTPUT
%   A         - Datafile
%
% DESCRIPTION
% Set the filenames for a datafile. This constructs a new datafile,
% except for the PREPROC and POSTPROC fields.
%
% SEE ALSO (PRTools Guide)
% DATAFILES

author: marianux

repo: ecg-kit

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


text:
function test_Script()

% GEOMETRIC PARAMETERS
pi = 4*atan(1);
L1 = 8;                              % length of computational domain (m)
N1 = 512;                            % number of Cartesian grid meshwidths at the finest level of the AMR grid
bell_length = 2;                     % bell length (m)
bell_circumference = pi;
npts_bell = ceil(2*(bell_length/L1)*N1);  % number of pos along the length of the bell
npts_circ = 1; %number of pos along the circumference (if in 3D)
npts = npts_bell*npts_circ;	    % total number pos
ds1 = bell_length/(npts_bell-1);   % mesh spacing(m) along length of bell

% Values from Alben, Peng, and Miller
betao = 0.5;
betam = 0.3;
to = 0.4;
Zs = L1/8;
xRef = L1/2;

t=0:0.025:5;

for i=1:length(t)

    %These are used to keep track of cycle number and time into the cycle
    pulse_time = t(i)-floor(t(i)); % determine time since beginning of first pulse

    % GIVE BELL STATES
    if (pulse_time<to)  %contract bell
        beta = betao+(betam-betao)*(pulse_time/to);
    else                % expand bell
        beta = betam+(betao-betam)*((pulse_time-to)/(1-to));
    end

    %---------------------
    %Xb and Yb are calculated here and will be used to determine new curvatures
    s1=0;
    zl = Zs;
    ro = xRef;

    % Pre-allocate memory for speed
    Xb_lam = zeros(npts,1);
    Yb_lam = zeros(npts,1);

    %top of bell
    Xb_lam(1)=ro;
    Yb_lam(1)=zl;

    %right side of bell
    for s1 = 2:ceil(npts_bell/2)
        thetaj = -1.55*(1-exp(-(s1)*ds1/beta));
        zl = zl + ds1*sin(thetaj);
        ro = ro + ds1*cos(thetaj);
        Xb_lam(s1)=ro;
        Yb_lam(s1)=zl;
    end

    zl = Zs;
    ro = xRef;

    %left side of bell
    for s1 = (ceil(npts_bell/2))+1:npts_bell
        s2=s1-(ceil(npts_bell/2));
        thetaj = -1.55*(1-exp(-(s2)*ds1/beta));
        zl = zl + ds1*sin(thetaj);
        ro = ro - ds1*cos(thetaj);
        Xb_lam(s1)=ro;
        Yb_lam(s1)=zl;
    end

    plot(Xb_lam,Yb_lam,'.'); hold on;
    axis([0 8 0 8]);
    pause(0.001);
    clf;

end

author: nickabattista


text:
function y = jed_to_year_hebrew ( jed )

%*****************************************************************************80
%
%% JED_TO_YEAR_HEBREW: the year in the Hebrew calendar when a JED occurred.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm H,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, page 331.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, the year in the Hebrew calendar that
%    included the JED.  If the input JED is less than the epoch of the Hebrew
%    calendar, then Y is always returned as -1.
%
  jed_epoch = epoch_to_jed_hebrew ( );

  if ( jed < jed_epoch )
    y = -1;
    return
  end
%
%  Using integer arithmetic in this computation may cause overflow.
%
%  Compute the number of months elapsed up to the date.
%
  m = 1 + floor ( ( 25920.0 * ( jed - jed_epoch + 2.5 ) ) / 765433.0 );
%
%  Estimate the number of years represented by these months.
%
  y = 19 * floor ( m / 235 ) + floor ( ( 19 * ( i4_modp ( m, 235 ) - 2 ) ) / 235 ) + 1;
  y = max ( y, 1 );
%
%  Determine the JED of the first day of that year.
%
  jed2 = new_year_to_jed_hebrew ( y );
%
%  We might have been off by 1 year.
%
  if ( jed < jed2 )
    y = y - 1;
  end

  return
end


text:
classdef ScalarProduct < handle

    properties (Access = public)
        epsilon
        Ksmooth
        Msmooth
    end

    properties (Access = private)
       nVariables
       mesh
       field
    end

    methods (Access = public)

        function obj = ScalarProduct(cParams)
            obj.init(cParams);
            obj.createMatrices(cParams);
        end

    end

    methods (Access = public)

        function sp = computeSP(obj,f,g)
            spM = obj.computeSP_M(f,g);
            spK = obj.computeSP_K(f,g);
            sp  = obj.epsilon^2*spK + spM;
        end

        function sp = computeSP_M(obj,f,g)
            sp = obj.computeProduct(obj.Msmooth,f,g);
        end

        function sp = computeSP_K(obj,f,g)
            sp = obj.computeProduct(obj.Ksmooth,f,g);
        end
    end

    methods (Access = private)

        function init(obj,cParams)
            obj.epsilon = cParams.epsilon;
            obj.nVariables = cParams.nVariables;
        end

        function createMatrices(obj,cParams)
            obj.mesh = cParams.mesh;
            M = obj.computeMassMatrix();
            K = obj.computeStiffnessMatrix();
            obj.Ksmooth = K;
            obj.Msmooth = M;
        end

        function n = computeProduct(obj,K,f,g)
            nx = length(f)/obj.nVariables;
            n = 0;
            for ivar = 1:obj.nVariables
                i0 = nx*(ivar-1) + 1;
                iF = nx*ivar;
                fs = f(i0:iF);
                gs = g(i0:iF);
                n = n + fs'*K*gs;
            end
        end

        function M = computeMassMatrix(obj)
            s.type  = 'MassMatrix';
            s.mesh  = obj.mesh;
            s.fun   = P1Function.create(obj.mesh, 1);
            s.quadratureOrder = 'QUADRATICMASS';
            LHS = LHSintegrator.create(s);
            M = LHS.compute();
        end

        function K = computeStiffnessMatrix(obj)
            s.type  = 'StiffnessMatrix';
            s.mesh  = obj.mesh;
            s.fun = P1Function.create(obj.mesh, 1);
            LHS = LHSintegrator.create(s);
            K = LHS.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function arrayIndices = sample2index(this, samplingPoints, iDims)
% Returns voxel samplingPoints corresponding to label (coordinate) samplingPoints
%
%   Y = MrDimInfo()
%   arrayIndices = sample2index(this, samplingPoints)
%
% Inverse operation to get_samplingPoints
%
% This is a method of class MrDimInfo.
%
% IN
%   samplingPoints     [nVoxels, nDims] of voxel samplingPoints
%                      (one per row) in coordinate system given by dimInfo
%   iDims              if specified, samplingPoints are assumed to be subset
%                      of all dimensions only
%
% OUT
%   arrayIndices        [nVoxels, nDims] of absolute
%                       voxel samplingPoints within array
%
% EXAMPLE
%   sample2index
%
%   See also MrDimInfo MrDimInfo.get_samplingPoints

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-01-23
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin < 3
    iDims = 1:this.nDims;
end

nDimsSubset = numel(iDims);

isRowVector = nDimsSubset == 1 && size(samplingPoints,1) == 1;

if isRowVector
   samplingPoints = samplingPoints(:); % allows row/column vectors as input for 1-dim trafo
end

nVoxels = size(samplingPoints,1);

arrayIndices = zeros(nVoxels,nDimsSubset);

for v = 1:nVoxels
    % find voxel index with closest (euclidean) sampling point in array
    for d = 1:nDimsSubset
        iDim = iDims(d);
        [~, arrayIndices(v,d)] = ...
            min(abs(this.samplingPoints{iDim} - samplingPoints(v,d)));
    end
end

if isRowVector % transform back to row vector for output
    arrayIndices = reshape(arrayIndices, 1, []);
end

author: translationalneuromodeling

repo: tapas


text:
%% SEMINAR: INTRODUCTION TO SIMULINK
% This seminar will show some of the basic and intermediate skills for Modeling,
% Designing and working with Simulink. This script will guide you through.
%%
%
%
Created by Roni Peer, 30.4.2012
Version 2.0

%

%% PART 1: BASICS
% This part will show how to use the basic abilities of Simulink to build, model,
% design and simulate a certain algorithm. Follow the steps bellow to do
% so.

%%% 1. Starting Simulink
% Let's open a new model and look what we can do with Simulink, and how
% it's built. We will build a small example to show the main features. The
% first model would show how to open Simulink, Drag blocks and connect
% them, configuring each block, running a simulation, and use of Stop
% Simulation.
simulink;

%%
% *To build SLWS_Step_01:*
%
%%
%
%
Step 01: Begin

%

%%
% # Open Simulink and show the Blocksets, search window, dragging of blocks.
% # Write the equation on the model: $F = ma ;  a = gravity; a = \dot{V}; V =
% \dot{X};$
% # Drag blocks: Constant, Integrator, Scope. Set $X_0 = 500; V_0 = 5; g =-9.81;$
% # Show how to zoom in/out/fit.
% # Find a way to stop the simulation when Altitude = 0;
% # Label the signals, and hide block names.
% # Change Simulation time to 100sec.
% # Simulate, and save the model.
%
% *To Open a complete step-1 model:*
%%
%
%
Step 01: Complete

%

%%% 2. Subsystems and libraries
% The next model would show some use of vector operations, subsystems,
% selector block, and finally to introduce Libraries.
% You can start with this model:
%%
%
%
Step 02: Begin

%

%%
% *To build SLWS_Step_02:*
%
% # Automaticaly align and resize the blocks.
% # Vectorize the gravity and signals: $g=$ [0 0 -9.81], $V_0=$ [8 -10 5],
% $X_0=$ [0 0 500]
% # Show the port data types.
% # Change solver to ODE23, Max. Step = 1sec.
% # Select the Z-component to stop simulation.
% # Create a subsystem for the "Stop Logic", and Change the Input signal to "Altitude"
% # Open Drag_Library and connect to the input with a "sum" block.
% # Define in MATLAB: $rho=1.2; Cd=0.333; A=1.5; m=150;$
% # Create Configurable Subsystem from the Drag Library
% In the end, you get this model:
%%
%
%
Step 02: Complete

%

%%% 3. Callbacks and BUSes
% The next model would show some use of BUS Signals, Masking Subsystems.
% You can start with this model:
%%
%
%
Step 03: Begin

%

%%
% *To build SLWS_Step_03:*
%
% # Create a callback which would load the parameters values ($Cd, A$...)
% # Create a mask for the Variable Density block. Show the constant one also.
% # Create a subsystem from the entire model - Name it Rocketeer.
% # Connect a BUS creator to velocity and Altitude, output them as "State".
% # Connect a constant block to acceleration. Value: $Thrust=$[0 0 1000]
% # Move the scope block to the top layer, and connect with a BUS Selector.
% # Run and verify.
% In the end, you get this model:
%%
%
%
Step 03: Complete

%

%%% 4. Closing the loop
% The next model would show how to continue to design the controller,
% and using a model reference.
% You can start with this model:
%%
%
%
Step 04: Begin

%

%%
% *To build SLWS_Step_04:*
%
% # Change "thrust" to "Desired Velocity", $V=$[0 0 -0.5]
% # Close the loop on velocity. Name the output "Error".
% # Connect a PID controller from Error to Plant. $PI =$ [50 5]
% # Connect scopes to "Thrust" and "Error".
% # Show block execution order. Discuss this concept. Treat Controller as
% Atomic.
% # Convert the controller to a model block, and replace.
% # Show difference between Controller and Plant - Hierarchy.
% # Run. View Error scope inside Controller. Why nothing? Why so slow?
% # Switch to Normal. Don't forget to inline parameters. What now?
% # Show the "Dependancy Viewer" and the "Generate Manifest"
% - You can talk more about controls in the link specified.
% In the end, you get this model:
%%
%
%
Step 04: Complete

%

%%% 5. "Slow Down!" - Control the rates and Embedded MATLAB
% The next model would show how to adjust the real-life settings of continuous
% and discrete systems. Also, how to incorprate Embedded MATLAB code into
% the model.
% You can start with this model:
%%
%
%
Step 05: Begin

%

%%
% *To build SLWS_Step_05:*
%
% # Set the PI Controller to Discrete, Sample Time  = 0.5sec.
% # Set the sample time display to Color.
% # Set Thrust outport and Velocity inport to inherit sample time.
% # Update, and look at the error.
% # Add rate transition blocks  before and after the model reference.
% # Change the controller's Solver to Fixed-Step Discrete.
% # Change the input to the controller to be a BUS. Define the BUS with the
% editor: StateBus - Position (3) and Velocity (3).
% # Save the BUS to a mat file, and save it. Use a callback to load it.
% # Include the Embedded MATLAB code in the Model.
% # Run.
%%
%
%
Open Embedded MATLAB Code

%
Step 05: Complete

%

%%% 6. Deploying the System
% Simulink allows you to generate different types of code.
% Open the following model, and generate these types of code:
% Embedded C-Code, HDL (VHDL or Verilog), and PLC Structured
% Text:
%%
%
%
Step 06: Deploy

%

%% PART 2: Multi-Domain Modeling
% This part will show how to use the more advanced abilities of Simulink to build, model,
% design and simulate your systems.
% Open the following models, and show how Simulink can aid in desiging
% multi-domain systems:

%%
%
%
Open Multi-Domain Model

%


text:
function x = one_null_right ( m, n )

%*****************************************************************************80
%
%% ONE_NULL_RIGHT returns a right null vector of the ONE matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of the matrix.
%
%    Output, real X(N,1), the null vector.
%
  if ( n == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ONE_NULL_RIGHT - Fatal error!\n' );
    fprintf ( 1, '  Matrix is nonsingular for N = 1.\n' );
    error ( 'ONE_NULL_RIGHT - Fatal error!' );
  end

  x = zeros ( n, 1 );

  x(1) =     1.0;
  x(n) =    -1.0;

  return
end


text:
% Q = vgg_line3d_Ppv(P)  Transforms a camera matrix to use it with Pluecker vector 3d line representation.
%
% P ... double(3,4), ordinary camera matrix
% Q ... double(6,3), transformed camera matrix (quadratic function of elements of P)
%   such that the projection of a 3d line by the camera is given by
%
%     l = Q*L
%
% where l is 1-by-3 image line vector, and L is 1-by-6 Pluecker vector of the 3d line.

% T.Werner

function Q = vgg_line3d_Ppv(P)

K = size(P,1)/3;
i1 = 1:3:3*K;
i2 = 2:3:3*K;
i3 = 3:3:3*K;

Q(:,i1) = vgg_line3d_pv_from_XY(P(i2,:)',P(i3,:)')';
Q(:,i2) = vgg_line3d_pv_from_XY(P(i3,:)',P(i1,:)')';
Q(:,i3) = vgg_line3d_pv_from_XY(P(i1,:)',P(i2,:)')';

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
classdef OptimizationMetricsPrinter_MMA < OptimizationMetricsPrinter

    methods (Access = protected)

        function printConvergenceVariables(obj,fid)
            kktnorm = obj.optimizer.historicalVariables.kktnorm;
            fprintf(fid,'Optimality tolerance: %f \n',kktnorm);
        end

    end

end


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

function Cstblood(mws)

material = invoke(mws,'material');
invoke(material,'Reset');
invoke(material,'Name','Blood'); %Blood
invoke(material,'FrqType','all');
invoke(material,'Type','Normal');
invoke(material,'SetMaterialUnit','GHz','mm');
invoke(material,'Epsilon','1');
invoke(material,'Mue','1');
invoke(material,'Sigma','0.0');
invoke(material,'TanD','0.0');
invoke(material,'TanDFreq','0.0');
invoke(material,'TanDGiven','False');
invoke(material,'TanDModel','ConstSigma');
invoke(material,'ConstTanDModelOrderEps','1');
invoke(material,'ReferenceCoordSystem','Global');
invoke(material,'CoordSystemType','Cartesian');
invoke(material,'SigmaM','0.0');
invoke(material,'TanDM','0.0');
invoke(material,'TanDMFreq','0.0');
invoke(material,'TanDMGiven','False');
invoke(material,'TanDMModel','ConstSigma');
invoke(material,'ConstTanDModelOrderMue','1');
invoke(material,'DispModelEps','None');
invoke(material,'DispModelMue','None');
invoke(material,'DispersiveFittingSchemeEps','Nth Order');
invoke(material,'MaximalOrderNthModelFitEps','3');
invoke(material,'ErrorLimitNthModelFitEps','0.01');
invoke(material,'DispersiveFittingSchemeMue','Nth Order');
invoke(material,'AddDispersionFittingValueEps','0.1', '76.8182454443825', '221.641134494173', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.2', '68.4738432550289', '115.059990632458', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.3', '65.6502510183341', '78.8589515742776', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.4', '64.1828321718796', '60.6489945539656', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.422222222222222', '63.9438618885212', '57.7730434712015', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.5', '63.257247554219', '49.7333267655039', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.6', '62.6027643145963', '42.5001241802731', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.7', '62.1030699559808', '37.3865932409505', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.744444444444444', '61.9142639813766', '35.5731389635495', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.8', '61.6997719763589', '33.6053481035262', '1.0');
invoke(material,'AddDispersionFittingValueEps','0.9', '61.3602871921475', '30.7161644656542', '1.0');
invoke(material,'AddDispersionFittingValueEps','1.06666666666667', '60.8862780821185', '27.206610582395', '1.0');
invoke(material,'AddDispersionFittingValueEps','1.38888888888889', '60.1510807057843', '23.0953297030954', '1.0');
invoke(material,'AddDispersionFittingValueEps','1.71111111111111', '59.5320763343363', '20.8376504484609', '1.0');
invoke(material,'AddDispersionFittingValueEps','2.03333333333333', '58.9648498262843', '19.545237054196', '1.0');
invoke(material,'AddDispersionFittingValueEps','2.35555555555556', '58.4208755792015', '18.8129113591261', '1.0');
invoke(material,'AddDispersionFittingValueEps','2.67777777777778', '57.8860468115685', '18.4294319222186', '1.0');
invoke(material,'AddDispersionFittingValueEps','3', '57.3529593274927', '18.4294319222186', '1.0');
invoke(material,'UseGeneralDispersionEps','True');
invoke(material,'UseGeneralDispersionMue','False');
invoke(material,'NLAnisotropy','False');
invoke(material,'NLAStackingFactor','1');
invoke(material,'NLADirectionX','1');
invoke(material,'NLADirectionY','0');
invoke(material,'NLADirectionZ','0');
invoke(material,'Rho','1060');
invoke(material,'ThermalType','Normal');
invoke(material,'ThermalConductivity','0.51');
invoke(material,'HeatCapacity','3.824');
invoke(material,'MetabolicRate','0');
invoke(material,'BloodFlow','1e+006');
invoke(material,'VoxelConvection','0');
invoke(material,'MechanicsType','Unused');
invoke(material,'Colour','1','0','0');
invoke(material,'Wireframe','False');
invoke(material,'Reflection','False');
invoke(material,'Allowoutline','False');
invoke(material,'Transparentoutline','False');
invoke(material,'Transparency','0');
invoke(material,'Create');
release(material);
end

repo: CST-MATLAB-API


text:
% acsobiro() - A. Chickocki's robust Second-Order Blind Identification (SOBI)
%              by joint diagonalization of the time-delayed covariance matrices.
%              NOTE: THIS CODE ASSUMES TEMPORALLY CORRELATED SIGNALS.
%              Thus, the estimated time-delayed covariance matrices
%              for at least some time delays must be nonsingular.
%
% Usage:  >>   [H] = acsobiro(X);
%         >> [H,S] = acsobiro(X,n,p);
% Inputs:
%         X - data matrix of dimension [m,N] where
%                    m is the number of sensors
%                    N is the number of samples
%         n - number of sources {Default: n=m}
%         p - number of correlation matrices to be diagonalized {default: 100}
%             For noisy data, use at least 100 time delays.
% Outputs:
%         H - matrix of dimension [m,n] an estimate of the *mixing* matrix
%         S - matrix of dimension [n,N] an estimate of the source activities
%             where  >> X [m,N] = H [m,n] * S [n,N]
%
% Authors: Implemented and improved by A. Cichocki on the basis of
%          the classical SOBI algorithm of Belouchrani and publications of:
%            A. Belouchrani et al., F. Cardoso et al.,
%            S. Choi, S. Cruces, S. Amari, and P. Georgiev
%          For references: see function body
%
% Note: Extended by Arnaud Delorme and Scott Makeig to process data epochs
%       (computes the correlation matrix respecting epoch boundaries).

% REFERENCES:
%  A. Belouchrani, K. Abed-Meraim, J.-F. Cardoso, and E. Moulines, ``Second-order
%  blind separation of temporally correlated sources,'' in Proc. Int. Conf. on
%  Digital Sig. Proc., (Cyprus), pp. 346--351, 1993.
%
%  A. Belouchrani, and A. Cichocki,
%  Robust whitening procedure in blind source separation context,
%  Electronics Letters, Vol. 36, No. 24, 2000, pp. 2050-2053.
%
%  A. Cichocki and S. Amari,
%  Adaptive Blind Signal and Image Processing, Wiley,  2003.

function [H,S,D]=acsobiro(X,n,p),

[m,N,ntrials]=size(X);
if nargin<1 | nargin > 3

  help acsobiro

elseif nargin==1,

 DEFAULT_LAGS = 100;
 n=m; % source detection (hum...)
 p=min(DEFAULT_LAGS,ceil(N/3)); % number of time delayed correlation matrices to be diagonalized
                                % Note: For noisy data, use at least p=100.
elseif nargin==2,

 p=min(DEFAULT_LAGS,ceil(N/3)); % number of correlation matrices to be diagonalized

end;

X(:,:)=X(:,:)-(mean(X(:,:)')'*ones(1,N*ntrials));        % Remove data means

for t = 1:ntrials
    if t == 1
        Rxx=(X(:,1:N-1,t)*X(:,2:N,t)')/(N-1)/ntrials; % Estimate the sample covariance matrix
                                          % for the time delay p=1, to reduce influence
                                          % of white noise.
    else
        Rxx=Rxx+(X(:,1:N-1,t)*X(:,2:N,t)')/(N-1)/ntrials; % Estimate the sample covariance matrix
                                          % for the time delay p=1, to reduce influence
                                          % of white noise.
    end;
end;

[Ux,Dx,Vx]=svd(Rxx);
        Dx=diag(Dx);

if n<m, % under assumption of additive white noise and when the number
         % of sources is known, or can be estimated a priori
  Dx=Dx-real((mean(Dx(n+1:m))));
  Q= diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
else    % under assumption of no additive noise and when the
        % number of sources is unknown
   n=max(find(Dx>1e-99)); % detect the number of sources
   fprintf('acsobiro(): Estimated number of sources is %d\n',n);
   Q= diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
end;
Xb = zeros(size(X));
Xb(:,:)=Q*X(:,:); % prewhitened data

% Estimate the time delayed covariance matrices:
 k=1;
 pn=p*n; % for convenience
 for u=1:m:pn,
   k=k+1;
   for t = 1:ntrials
       if t == 1
           Rxp=Xb(:,k:N,t)*Xb(:,1:N-k+1,t)'/(N-k+1)/ntrials;
       else
           Rxp=Rxp+Xb(:,k:N,t)*Xb(:,1:N-k+1,t)'/(N-k+1)/ntrials;
       end;
   end;
   M(:,u:u+m-1)=norm(Rxp,'fro')*Rxp;  % Frobenius norm =
 end;                                  % sqrt(sum(diag(Rxp'*Rxp)))

% Approximate joint diagonalization:
eps=1/sqrt(N)/100; encore=1; U=eye(n);
while encore, encore=0;
 for p=1:n-1,
  for q=p+1:n,
    % Givens rotations:
    g=[ M(p,p:n:pn)-M(q,q:n:pn)  ;
        M(p,q:n:pn)+M(q,p:n:pn)  ;
        i*(M(q,p:n:pn)-M(p,q:n:pn))];
   [Ucp,D] = eig(real(g*g')); [la,K]=sort(diag(D));
   angles=Ucp(:,K(3));angles=sign(angles(1))*angles;
   c=sqrt(0.5+angles(1)/2);
   sr=0.5*(angles(2)-j*angles(3))/c; sc=conj(sr);
   asr = abs(sr)>eps ;
   encore=encore | asr ;
   if asr , % Update the M and U matrices:
     colp=M(:,p:n:pn);
     colq=M(:,q:n:pn);
     M(:,p:n:pn)=c*colp+sr*colq;
     M(:,q:n:pn)=c*colq-sc*colp;
     rowp=M(p,:);
     rowq=M(q,:);
     M(p,:)=c*rowp+sc*rowq;
     M(q,:)=c*rowq-sr*rowp;
     temp=U(:,p);
     U(:,p)=c*U(:,p)+sr*U(:,q);
     U(:,q)=c*U(:,q)-sc*temp;
   end  %% if
  end  %% q loop
 end  %% p loop
end  %% while

% Estimate the mixing matrix H
H= pinv(Q)*U(1:n,1:n);

% Estimate the source activities S
if nargout>1
  S=[];
  W=U(1:n,1:n)'*Q;
  S= W*X(:,:);
end

author: PatternRecognition

repo: OpenBMI


text:
function [out] = prep_selectTrials(dat,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prep_selectTrials (Pre-processing procedure):
%
% Synopsis:
%     [out] = prep_selectTrials(DAT,<OPT>)
%
% Example :
%     out = prep_selectTrials(dat,{'Index',[20:35]});
%
% Arguments:
%     dat - Structure. epoched data
%         - Data which trials are to be selected
%     varargin - struct or property/value list of optional properties:
%          :  index - index of trials to be selected
%
% Returns:
%     out - Data structure which has selected channels from epoched data
%
% Description:
%     This function selects data of specified trials
%     from  epoched data.
%     continuous data should be [time * channels]
%     epoched data should be [time * channels * trials]
%
% See also 'https://github.com/PatternRecognition/OpenBMI'
%
% Min-ho Lee, 12-2017
% mh_lee@korea.ac.kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(varargin)
    warning('OpenBMI: Trials should be specified');
    out = dat;
    return
end
if iscell(varargin{:})
    opt=opt_cellToStruct(varargin{:});
elseif isstruct(varargin{:}) % already structure(x-validation)
    opt=varargin{:}
end

if ~isfield(dat, 'x')
    warning('OpenBMI: Data structure must have a field named ''x''');
    return
end
if ~isfield(dat, 't')
    warning('OpenBMI: Data structure must have a field named ''t''');
    return
end
if ~isfield(dat, 'y_dec') || ~isfield(dat, 'y_logic') || ~isfield(dat, 'y_class')
    warning('OpenBMI: Data structure must have a field named ''y_dec'',''y_logic'',''y_class''');
    return
end

idx = opt.Index;
nd = ndims(dat.x);
if nd == 3
    x = dat.x(:,idx,:);
elseif nd==2
    if size(dat.chan,2)==1
        x = dat.x(:,idx);
        warning('OpenBMI: just 1 channel data?')
    else
        x=dat.x;
    end
elseif nd ==1
    x = dat.x;
else
    warning('OpenBMI: Check for the data dimensionality')
    return
end
out = rmfield(dat,{'x','t','y_dec','y_logic','y_class'});
out.x = x;
out.t = dat.t(idx);
out.y_dec = dat.y_dec(idx);
out.y_logic = dat.y_logic(:,idx);
out.y_class = dat.y_class(idx);


text:
function [cl,cu] = matRad_getConstraintBounds(optiProb,cst)
% matRad IPOPT get constraint bounds function for direct aperture optimization
%
% call
%   [cl,cu] = matRad_daoGetConstBounds(optiProb,cst)
%
% input
%   optiProb:   option struct defining the type of optimization
%   cst:        matRad cst struct
%
% output
%   cl:         lower bounds on constraints
%   cu:         lower bounds on constraints
%
% References
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

apertureInfo = optiProb.apertureInfo;

% Initialize bounds
cl_dao = zeros(apertureInfo.totalNumOfLeafPairs,1);
cu_dao = inf*ones(apertureInfo.totalNumOfLeafPairs,1);

% get dosimetric bounds from cst (just like for conv opt) by call to
% superclass method
[cl_dos,cu_dos] = matRad_getConstraintBounds@matRad_OptimizationProblem(optiProb,cst);

% concatenate
cl = [cl_dao; cl_dos];
cu = [cu_dao; cu_dos];

repo: matRad


text:
function makefigure(y, ThreeDee, PlotOrClose)
% Copyright (c) 2009, The MathWorks, Inc.
%
% Creates / updates /closes a MATLAB figure
persistent figureHandle;

if ishandle(figureHandle)
    % If figure exists, use it
    figure(figureHandle);
    if (nargin > 2)&&(PlotOrClose < -eps)
        % close the figure if needed
        close(figureHandle);
        return;
    end;
else
    % If figure doesn't exist and plotting is required, create it.
    if (nargin > 2)&&(PlotOrClose > eps)
        figureHandle = figure;
    else
        return;
    end;
end;

% Clear the figure and plot in 2D or 3D
clf;
if ThreeDee
    surf(y);
else
    plot(y);
end;


text:
function [endpointsAll,bboxAll] = backprojectSpatialImg(expidx,firstidx,nImgs,bSave)

fprintf('backprojectSpatialImg()\n');

if (ischar(expidx))
    expidx = str2double(expidx);
end

if (ischar(firstidx))
    firstidx = str2double(firstidx);
end

if (ischar(nImgs))
    nImgs = str2double(nImgs);
end

if (nargin < 4)
    bSave = true;
elseif (ischar(bSave))
    bSave = str2double(bSave);
end

% load annolist with full size images
[annolistOrig, imgidxs] = getAnnolist(expidx);
for imgidx = 1:length(annolistOrig)
    annolistOrig(imgidx).annorect = [];
end

lastidx = firstidx + nImgs - 1;
if (lastidx > length(imgidxs))
    lastidx = length(imgidxs);
end

p = rcnn_exp_params(expidx);

if (isfield(p,'nms_thresh_det'))
    nms_thresh_det = p.nms_thresh_det;
else
    nms_thresh_det = 0.5;
end

fprintf('expidx: %d\n',expidx);
fprintf('firstidx: %d\n',firstidx);
fprintf('nImgs: %d\n',nImgs);
fprintf('nms_thresh: %1.1f\n',nms_thresh_det);

load(p.testGT);
if (~exist('annolist','var'))
    annolist = single_person_annolist;
end
annolist_crops = annolist;

[~,parts] = util_get_parts24();

saveTo = [fileparts(p.evalTest) '/pred_pairwise_backproj'];
if (~exist(saveTo,'dir'))
    mkdir(saveTo);
end

conf = rcnn_config('sub_dir', '/cachedir/test', 'exp_dir', [p.expDir '/' p.shortName]);
predPairwiseDir = [conf.cache_dir '/pred-pairwise'];

for imgidx = firstidx:lastidx
    endpointsAll(imgidx).imgname = annolistOrig(imgidxs(imgidx)).image.name;
    endpointsAll(imgidx).det = cell(length(p.spatidxs)+1,1);
end

sticksMap = [0 2; 2 4; 7 5; 9 7; nan nan; 22 23; 12 14; 14 16; 19 17; 21 19; 16 22; 17 22; 4 22; 5 22];

bVis = true;

ii = 1;
for i = firstidx:lastidx
    fprintf('.');

    [~,name] = fileparts(annolist_crops(ii).image.name);
    imgidx_crop = str2double(name(3:7));
    while (imgidx_crop ~= imgidxs(i))
        ii = ii + 1;
        [~,name] = fileparts(annolist_crops(ii).image.name);
        imgidx_crop = str2double(name(3:7));
    end

    % load bboxes
    while (imgidx_crop == imgidxs(i))

        fname = [predPairwiseDir '/imgidx_' padZeros(num2str(ii-1),5) '_nms'];
        load(fname,'det_nms_union_07_all','idxs_nms_union_07_all');
        det_spatial = det_nms_union_07_all;
        idxs_spatial = idxs_nms_union_07_all;

%         load(fname,'det_nms_inters_03_all','idxs_nms_inters_03_all');
%         det_spatial = det_nms_inters_03_all;
%         idxs_spatial = idxs_nms_inters_03_all;

%         load(fname,'det_nms_inters_05_all','idxs_nms_inters_05_all');
%         det_spatial = det_nms_inters_05_all;
%         idxs_spatial = idxs_nms_inters_05_all;

        % transformation matrix
        load([fileparts(p.testGT) '/T_' name(3:end)],'T');
        for j = 1:length(p.spatidxs)

            spatidxs = p.spatidxs{j};

            % backproject det
            det = det_spatial{j}(:,1:5);
            detNew1 = ([det(:,1:2) ones(size(det,1),1)]*T');
            detNew2 = ([det(:,3:4) ones(size(det,1),1)]*T');
            det(:,1:4) = [detNew1(:,1:2) detNew2(:,1:2)];

            idxStick = find(sum(abs(repmat(spatidxs,length(sticksMap),1)-sticksMap),2) == 0 | ...
                        sum(abs(repmat(spatidxs([2 1]),length(sticksMap),1)-sticksMap),2)==0);

            pidx1 = sticksMap(idxStick,1);
            pidx2 = sticksMap(idxStick,2);

            class_ids(1) = find(pidx1 == p.pidxs);
            class_ids(2) = find(pidx2 == p.pidxs);

            idx1 = find(idxs_spatial{j}(1,1:2) == class_ids(1));
            idx2 = find(idxs_spatial{j}(1,1:2) == class_ids(2));
            if (idx1 == 2 && idx2 == 1)
                det(:,1:4) = det(:,[3:4 1:2]);
            end

            % add to detections
            endpointsAll(i).det{idxStick} = [endpointsAll(i).det{idxStick}; det];
        end
%         x2y2 = (endpointsAll(i).det{11}(:,1:2) + endpointsAll(i).det{12}(:,1:2))/2;
%         x1y1 = (endpointsAll(i).det{13}(:,1:2) + endpointsAll(i).det{14}(:,1:2))/2;
%         endpointsAll(i).det{5} = [x1y1 x2y2 endpointsAll(i).det{11}(:,5)];
        ii = ii + 1;
        if (ii > length(annolist_crops))
            break;
        end
        [~,name] = fileparts(annolist_crops(ii).image.name);
        imgidx_crop = str2double(name(3:7));
    end

    if (bVis)
        clf; imagesc(imread(endpointsAll(i).imgname));axis equal;hold on;
        for j = 1:length(endpointsAll(i).det)
            det = endpointsAll(i).det{j};
            if (~isempty(det))
                [val,idx] = sort(det(:,end),'descend');
                plot(det(idx(1),[1 3]),det(idx(1),[2 4]),'r-','lineWidth',3);
%                 text(double(det(idx(1),1)) + 30,double(det(idx(1),2)),sprintf('%1.1f',det(idx(1),end)),'BackgroundColor','w','verticalalignment','top','horizontalalignment','left','fontSize',14);

                plot(det(idx(1),1),det(idx(1),2),'bo','MarkerSize',10,'MarkerFaceColor','b','MarkerEdgeColor','k');
                plot(det(idx(1),3),det(idx(1),4),'go','MarkerSize',10,'MarkerFaceColor','g','MarkerEdgeColor','k');
            end
        end
    end

    if (bSave)
        fname = [saveTo '/keypoints_' padZeros(num2str(imgidxs(i)),5)];
        keypoints = endpointsAll(i);
        save(fname,'keypoints');

        fname = [saveTo '/bboxes_' padZeros(num2str(imgidxs(i)),5)];
        bboxes = bboxAll(i);
        save(fname,'bboxes');
    end

    if (~mod(i, 100))
        fprintf(' %d/%d\n',i,length(imgidxs));
    end
end

fprintf(' done\n');

if (~exist('endpointsAll','var'))
    endpointsAll = [];
    bboxAll = [];
end
% path = fileparts(p.evalTest);
% fnameKeypoints = [path '/keypointsAll'];
% save(fnameKeypoints, 'keypointsAll');
%
% fnameBBox = [path '/bboxAll'];
% save(fnameBBox, 'bboxAll');
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function msm_to_mm_array_integer ( output_filename, a, symmetry )

%*****************************************************************************80
%
%% MSM_TO_MM_ARRAY_INTEGER writes a "matrix array integer" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
%    Input, string SYMMETRY, describes the symmetry type, as 'GENERAL',
%    'SYMMETRIC', 'SKEW-SYMMETRIC' or 'HERMITIAN'.
%
  if ( s_eqi ( symmetry, 'GENERAL' ) )
    msm_to_mm_array_integer_general ( output_filename, a );
  elseif ( s_eqi ( symmetry, 'SYMMETRIC' ) )
    msm_to_mm_array_integer_symmetric ( output_filename, a );
  elseif ( s_eqi ( symmetry, 'SKEW-SYMMETRIC' ) )
    msm_to_mm_array_integer_skewsymmetric ( output_filename, a );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_ARRAY_INTEGER - Fatal error!\n' );
    fprintf ( 1, '  Illegal symmetry option.\n' );
    error ( 'MSM_TO_MM_ARRAY_INTEGER - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function LU = getbounds(F,avoidequalitybounds)

if nargin == 1
    LU = getbounds(lmi(F));
else
    LU = getbounds(lmi(F),avoidequalitybounds);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% The COBRAToolbox: testLegalRxnFormula.m
%
% Purpose:
%     - testReactionEq tests the functionality of LegalRxnFormula in the
%     rBioNet extension
%
% Authors:
%     - Stefania Magnusdottir April 2017
%

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testReactionEq'));
cd(fileDir);

% load E. coli model
model = getDistributedModel('ecoli_core_model.mat');

% check validity of reaction formulas
formula = printRxnFormula(model, model.rxns(1));
result = LegalRxnFormula(formula, model.rxns(1));

% test output
assert(result)

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [ PD ] = PD_Orthogonal( X,a,SNR )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
[N,L] = size(X);
Rs = X*X'/L;
rou = SNR*abs(a'*Rs.'*a)^2;
Pfa = 1e-7;
ita = chi2inv(1-Pfa,2);
PD = 1-ncx2cdf(ita,2,rou);
end

author: yuanhao-cui

repo: Must-Reading-on-ISAC

path: github-repos/MATLAB/yuanhao-cui-Must-Reading-on-ISAC/Must-Reading-on-ISAC-34cd6615c52ebca121428a979e756c608b195040/Codes/Fan2018TSP/Codes for DFRC Waveform Design/Waveform Design With Given Radar Beampatterns/PD_Orthogonal.m


text:
function java_setcb(jObj, varargin)
%JAVA_SETCB: Associate Matlab callbacks to a Java object.
%
% USAGE:  java_setcb(jObj, CallbackName1, CallbackFunction1, CallbackName2, CallbackFunction2, ...)

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

% Old syntax, before Matlab2010b
% for i = 1:2:length(varargin)
%     set(jObj, varargin{i}, varargin{i+1});
% end

% New syntax: Matlab 2010b
hObj = handle(jObj, 'callbackProperties');
for i = 1:2:length(varargin)
    set(hObj, varargin{i}, varargin{i+1});
end

author: brainstorm-tools


text:
%CAMPOSE4	Camera pose estimation from 4 coplanar points.
%
%	T = CAMPOSE4(D, ci)
%	[T, UV0] = CAMPOSE4(D, ci)
%
% Input is a table of data points, D, with each row of the form [X Y u v]
% where (X, Y) are the world coordinate of the planar points with respect
% to some coordinate system whose origin lies within the plane,
% and (u, v) are the corresponding image  plane coordinate.
%
% ci is a structure of camera intrinsic parameters that must contain
% focal length (f), and pixel sizes (sx, sy).
%
% Output is a homogeneous transformation, T, of the camera's pose with
% respect to the coordinate frame of the planar points.
% The optional result, UV0, is the coordinate of the principal point (in
% distance units).
%
% from Ganapathy "Camera Location Determination Problem",
% Bell Labs Tech. Memo 11358-841102-20-TM, Nov 2 1984
%
% SEE ALSO: CAMPOSE4, CAMCALP, CAMERA, CAMCALT, INVCAMCAL

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

function [T, uv0] = campose4(uvXY, ci)

	t = camcald4(uvXY);
    k(1) = ci.f / ci.sx;
    k(2) = ci.f / ci.sy;
% K = [k1 k2] contains the X and Y-direction scale
% factors: k1 = kx * F, k2 = ky * F, where F is the focal length.
%

	lam1 = t(1,1)*t(3,2) - t(1,2)*t(3,1);
	lam2 = t(2,1)*t(3,2) - t(2,2)*t(3,1);
	if abs(lam1) < 1e-8,
		lam1 = 0;
	end
	if abs(lam2) < 1e-8,
		lam2 = 0;
	end

	if (lam1 == 0) & (lam2 == 0),
		r2 = (k(1)^2 + k(2)^2) / sqrt(t(1,1)^2 + t(1,2)^2 + t(2,1)^2 + t(2,2)^2);
		r = sqrt(r2);
		disp('Problem with solution: lam1 = lam2 = 0');
		i = 1;
		f = 0;
		c = 0;
		g = 0;
		h = 0;

		u0 = ci.u0;
		v0 = ci.v0;
	else
		r2 = (t(3,1)^2 + t(3,2)^2)/( (lam1/k(1))^2 + (lam2/k(2))^2 );
		r = sqrt(r2);

		i2 = 1 - r2*(t(3,1)^2 + t(3,2)^2);
		i = sqrt(i2);
		f = -lam1/k(1)*r2;
		c = lam2/k(2)*r2;

		u0 = (t(1,1)*t(3,1) + t(1,2)*t(3,2) + k(1)*c*i/r2) / (t(3,1)^2 + t(3,2)^2);
		v0 = (t(2,1)*t(3,1) + t(2,2)*t(3,2) + k(2)*f*i/r2) / (t(3,1)^2 + t(3,2)^2);

		g = t(3,1)*r;
		h = t(3,2)*r;
	end

	a = (t(1,1)*r - u0*g) / k(1);
	b = (t(1,2)*r - u0*h) / k(1);
	d = (t(2,1)*r - v0*g) / k(2);
	e = (t(2,2)*r - v0*h) / k(2);
	p = (t(1,3)*r - u0*r) / k(1);
	q = (t(2,3)*r - v0*r) / k(2);

	R = [a b c;d e f;g h i];
	P = [p q r];

	T = [R P'; 0 0 0 1];
	uv0 = [u0 v0];

author: petercorke

repo: machinevision-toolbox-matlab


text:
function perf=mse(e)
%
% calculate the mean squared error of the given errors
%
%  'perf = mse(E);'
%
% see also:
%    mae, linf, trimmedmse
%

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

perf = sum(sum(e.^2)) / numel(e);

author: peterhcharlton

repo: RRest