text:
function [im, gx, gy, gz, midx] = scimat_intersect_plane(scimat, m, v, interp)
% SCIMAT_INTERSECT_PLANE  Intersection of a plane with an image volume.
%
% [IM, GX, GY, GZ, MIDX] = scimat_intersect_plane(SCIMAT, M, V)
%
%   SCIMAT is a struct with the 3D volume that we want to intersect. We use
%   SCIMAT structs widely in Gerardus, because that way we have the image
%   data and metainformation (e.g. voxel size) together in the same
%   variable. For details on SCIMAT structs, see "help scimat".
%
%   IM is an image that displays the intersection of the plane with the
%   image volume. Voxels that fall outside the image volume are returned as
%   NaN.
%
%   GX, GY, GZ are matrices of the same size as IM, and contain the
%   Cartesian coordinates of the voxels in IM (i.e. the coordinates of the
%   points at which the SCIMAT volume was sampled to obtain IM). You can
%   visualize the resulting plane using
%
%     >> surf(gx, gy, gz, im, 'EdgeColor', 'none')
%
%   Note: When INTERP='nn' (see below), coordinates in the plane are
%   rounded to the nearest image voxel centre. Hence, if you plot the plane
%   as above, it will in general look slightly "jagged" instead of flat.
%
%   The plane is uniquely defined in 3D space using a point and a vector:
%
%     * M is a 3-vector with the coordinates of a point contained in the
%       plane. It is assumed that M is within the image boundaries
%
%     * V is a 3-vector that represents a normalized vector orthogonal to
%       the plane
%
%   MIDX is a 2-vector with the row and colum of M in the intersecting
%   plane. That is, the X, Y, Z coordinates of M are
%
%     [gx(midx(1), midx(2)), ...
%      gy(midx(1), midx(2)), ...
%      gz(midx(1), midx(2))]
%
% ... = scimat_intersect_plane(..., INTERP)
%
%   INTERP is a string with the interpolation method when sampling the
%   volume with the plane:
%
%     'nn' (default): nearest neighbour. Good for binary segmentation masks
%     'linear': linear interpolation. Good for grayscale images

% Authors: Ramon Casero <rcasero@gmail.com>,
% Pablo Lamata <pablo.lamata@dpag.ox.ac.uk>
% Copyright © 2010-2015 University of Oxford
% Version: 0.4.4
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(3, 4);
nargoutchk(0, 5);

% defaults
if (nargin < 4 || isempty(interp))
    interp = 'nn';
end
% remove dummy dimension of data if necessary
scimat = scimat_squeeze(scimat, true);

% the longest distance within the image volume is the length of the largest
% diagonal. We compute it in voxel units
lmax = ceil(sqrt(sum(([scimat.axis.size] - 1).^2)));

% create horizontal grid for the plane sampling points. We are going to have one
% sampling point per voxel. The grid is centered on 0
[gc, gr] = meshgrid(-lmax:lmax, -lmax:lmax);
gs = 0 * gc;

% compute a rotation matrix to map the Z-axis to v. Note that v is given in
% x, y, z coordinates, but we want to work with r, c, s indices, i.e. y, x,
% z coordinates
v = v([2 1 3]);
rotmat = vec2rotmat(v(:));

% rotate the grid sampling points accordingly
grcs = rotmat * [gr(:) gc(:) gs(:)]';

% compute index coordinates of real world coordinates of the plane point
idxm = scimat_world2index(m, scimat);

% center rotated grid on the plane point m
grcs(1, :) = grcs(1, :) + idxm(1);
grcs(2, :) = grcs(2, :) + idxm(2);
grcs(3, :) = grcs(3, :) + idxm(3);

% sampling points that are outside the image domain
idxout = (grcs(1, :) < 1) | (grcs(1, :) > scimat.axis(1).size) ...
    | (grcs(2, :) < 1) | (grcs(2, :) > scimat.axis(2).size) ...
    | (grcs(3, :) < 1) | (grcs(3, :) > scimat.axis(3).size);

% sampling points that are inside
idxin = ~idxout;
im = nan(size(grcs, 2), 1);
switch interp
    case 'nn' % nearest neighbour
        % round coordinates, so that we are sampling at voxel centers and don't
        % need to interpolate
        grcs = round(grcs);
        % rcs indices => linear indices
        idx = sub2ind(size(scimat.data), grcs(1, idxin), grcs(2, idxin), ...
        grcs(3, idxin));
        % sample image volume with the rotated and translated plane
        im(idxin) = scimat.data(idx);
    case 'linear' % linear interpolation
        xi = grcs(1, idxin);
        yi = grcs(2, idxin);
        zi = grcs(3, idxin);
        im(idxin) = interp3(scimat.data, yi, xi, zi, interp);
    otherwise
        error('Interpolation method not implemented')
end

% compute real world coordinates for the sampling points
gxyz = scimat_index2world(grcs', scimat)';

% reshape the sampled points to get again a grid distribution
im = reshape(im, size(gc));
gx = reshape(gxyz(1, :), size(gr));
gy = reshape(gxyz(2, :), size(gc));
gz = reshape(gxyz(3, :), size(gs));

% keep track of where the rotation point is using a boolean vector for the
% row position, and another one for the column position. The rotation point
% is at the center of the grid
v0row = false(size(gr, 1), 1);
v0row((size(gr, 1) + 1) / 2) = true;
v0col = false(1, size(gc, 2));
v0col((size(gc, 2) + 1) / 2) = true;

% find columns where all elements are NaNs
idxout = all(isnan(im), 1);

% remove those columns
im = im(:, ~idxout);
gx = gx(:, ~idxout);
gy = gy(:, ~idxout);
gz = gz(:, ~idxout);
v0col = v0col(~idxout);

% find rows where all elements are NaNs
idxout = all(isnan(im), 2);

% remove those columns
im = im(~idxout, :);
gx = gx(~idxout, :);
gy = gy(~idxout, :);
gz = gz(~idxout, :);
v0row = v0row(~idxout);

% find row, column coordinates of the rotation point in the intersection
% plane
midx = [find(v0row) find(v0col)];

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% Machine Learning Online Class - Exercise 1: Linear Regression

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%
% x refers to the population size in 10,000s
% y refers to the profit in $10,000s
%

%% Initialization
clear ; close all; clc

%% ==================== Part 1: Basic Function ====================
% Complete warmUpExercise.m
fprintf('Running warmUpExercise ... \n');
fprintf('5x5 Identity Matrix: \n');
warmUpExercise()

%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ======================= Part 2: Plotting =======================
fprintf('Plotting Data ...\n')
data = load('ex1data1.txt');
X = data(:, 1); y = data(:, 2);
m = length(y); % number of training examples

% Plot Data
% Note: You have to complete the code in plotData.m
plotData(X, y);

%fprintf('Program paused. Press enter to continue.\n'); pause;

%% =================== Part 3: Gradient descent ===================
fprintf('Running Gradient Descent ...\n')

X = [ones(m, 1), data(:,1)]; % Add a column of ones to x
theta = zeros(2, 1); % initialize fitting parameters

% Some gradient descent settings
iterations = 1500;
alpha = 0.01;

% compute and display initial cost
computeCost(X, y, theta)

% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent: ');
fprintf('%f %f \n', theta(1), theta(2));

% Plot the linear fit
hold on; % keep previous plot visible
plot(X(:,2), X*theta, '-')
legend('Training data', 'Linear regression')
hold off % don't overlay any more plots on this figure

% Predict values for population sizes of 35,000 and 70,000
predict1 = [1, 3.5] *theta;
fprintf('For population = 35,000, we predict a profit of %f\n',...
    predict1*10000);
predict2 = [1, 7] * theta;
fprintf('For population = 70,000, we predict a profit of %f\n',...
    predict2*10000);

%fprintf('Program paused. Press enter to continue.\n'); pause;

%% ============= Part 4: Visualizing J(theta_0, theta_1) =============
fprintf('Visualizing J(theta_0, theta_1) ...\n')

% Grid over which we will calculate J
theta0_vals = linspace(-10, 10, 100);
theta1_vals = linspace(-1, 4, 100);

% initialize J_vals to a matrix of 0's
J_vals = zeros(length(theta0_vals), length(theta1_vals));

% Fill out J_vals
for i = 1:length(theta0_vals)
    for j = 1:length(theta1_vals)
	  t = [theta0_vals(i); theta1_vals(j)];
	  J_vals(i,j) = computeCost(X, y, t);
    end
end

% Because of the way meshgrids work in the surf command, we need to
% transpose J_vals before calling surf, or else the axes will be flipped
J_vals = J_vals';
% Surface plot
figure;
surf(theta0_vals, theta1_vals, J_vals)
xlabel('\theta_0'); ylabel('\theta_1');

% Contour plot
figure;
% Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100
contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))
xlabel('\theta_0'); ylabel('\theta_1');
hold on;
plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);

author: worldveil

repo: coursera-ml


text:
classdef  VideoMaker_TopOpt_density3DZmirror < VideoMaker_TopOpt
  properties

  end

  methods (Access = public)

        function Make_video_design_variable(obj,output_video_name)
            post = Postprocess_TopOpt_density();
            field2print = post.density_name;
            componentfield = post.density_name_component;
            obj.Make_video_density_field(field2print,componentfield,output_video_name)
        end

        function Make_video_design_variable_reg(obj,output_video_name)
            post = Postprocess_TopOpt_density();
            field2print = post.density_name_reg;
            componentfield = post.density_name_component_reg;
            obj.Make_video_density_field(field2print,componentfield,output_video_name)
        end

  end

  methods (Access = private)
        function Make_video_density_field(obj,field2print,componentfield,output_video_name)
            file_tcl_name = 'tcl_gid.tcl';
            file_list = obj.create_file_list(obj.iterations_to_print,obj.file_name,obj.files_folder);
            [output_video_name] = obj.replace_special_character(output_video_name);
            [output_photo]=obj.replace_special_character(strcat(obj.files_folder,'\',obj.file_name,'.png'));
            [file_list] = obj.replace_special_character(file_list);

            file_tcl_name_with_path = fullfile(obj.files_folder,file_tcl_name);
            file_path_in = fullfile(pwd,'FEM','PostProcess','Make_Video_density3DZmirror.tcl');
            filepath = obj.replace_special_character(file_path_in);

            fid = fopen(file_tcl_name_with_path,'w+');

            fprintf(fid,'GiD_Process PostProcess \n');
            fprintf(fid,'%s\n',['set arg1 "',file_list,'"']);
            fprintf(fid,'%s\n',['set arg2 "',output_video_name,'"']);
            fprintf(fid,'%s\n',['set arg3 "',field2print,'"']);
            fprintf(fid,'%s\n',['set arg4 "',componentfield,'"']);
            fprintf(fid,'%s\n',['set arg5 "',output_photo,'"']);

            fprintf(fid,'%s\n',['source "',filepath,'"']);
            fprintf(fid,'%s\n',['Make_Video_density $arg1 $arg2 $arg3 $arg4 $arg5']);
            fprintf(fid,'%s\n',['GiD_Process Mescape Quit']);
            fclose(fid);
            obj.execute_tcl_files(obj.gidPath,file_tcl_name_with_path)
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function h=hor_slice(xvect,yvect,data,const)
%FUNCTION H=VERT_SLICE(HVECT,VVECT,DATA,CONST)
% plots the 2D data of, say a seismic survey, GPR,
% resistivity survey in 3D as vertical slices
% For multiple vertical slice note that this only work
% for surveys that are mutually perpendicular.
% You have to provide:
% - HVECT: a vector with the 1st horizontal data (say offset in GPR)
% - VVECT: a vector with values for the 2nd horizontal axis.
% - DATA: the data to be displayed (make sure that the dimensions
%   match with the two vectors
% - CONST: a constant defining the depth of the slice (intersection point
%   with the z-axis.
%
% The function returns the handle of the surface plot.
%
% CAREFUL: Plotting too many slices kills your computer when dong graphical
% adjustments (rotation rescaling, ...) - you will have to be very patient
% ... or you need a faster machine.
%
% Ulrich Theune, 2003
%
%

if ~ishold
    hold on
end

[m,n]=size(data);

z=const*ones(m,n);
[x,y]=meshgrid(xvect,yvect);
h=surf(x,y,z,data);
plot3([min(xvect) min(xvect) max(xvect) max(xvect) min(xvect)],...
    [min(yvect) max(yvect) max(yvect) min(yvect) min(yvect)],...
    [const const const const const],'k','linewidth',0.1);
shading interp

clear x y z
view(3)

author: Sable

repo: mcbench-benchmarks


text:
function index = findPoint(coord, points, varargin)
%FINDPOINT Find index of a point in an set from its coordinates.
%
%   IND = findPoint(POINT, ARRAY)
%   Returns the index of point whose coordinates match the 1-by-2 row array
%   POINT in the N-by-2 array ARRAY. If the point is not found, returns 0.
%   If several points are found, keep only the first one.
%
%   If POINT is a M-by-2 array, the result is a M-by-1 array, containing
%   the index in the array of each point given by COORD, or 0 if the point
%   is not found.
%
%   IND = findPoint(POINT, ARRAY, TOL)
%   use specified tolerance, to find point within a distance of TOL.
%   Default tolerance is zero.
%
%   See also
%    points2d, minDistancePoints, distancePoints, findClosestPoint

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-07-17
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

% number of points
np = size(coord, 1);

% allocate memory for result
index = zeros(np, 1);

% specify the tolerance
tol = 0;
if ~isempty(varargin)
    tol = varargin{1};
end

if tol == 0
    for i = 1:np
        % indices of matches
        ind = find(points(:,1) == coord(i,1) & points(:,2) == coord(i,2));

        % format current result
        if isempty(ind)
            index(i) = 0;
        else
            index(i) = ind(1);
        end
    end
else
    for i = 1:np
        % indices of matches
        ind = find(sqrt(sum(bsxfun(@minus, points, coord) .^ 2, 2)) <= tol);

        % format current result
        if isempty(ind)
            index(i) = 0;
        else
            index(i) = ind(1);
        end
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function XYZtransform = transformPointCloud(XYZ,Rt)
    XYZtransform = Rt(1:3,1:3) * XYZ + repmat(Rt(1:3,4),1,size(XYZ,2));
end

author: jianxiongxiao

repo: ProfXkit


text:
function c = mrdivide(a,b)
% MRDIVIDE implements a/b, where either a or b is an adiff object.
% It is mapped to a./b where possible (i.e. when b has length 1), and
% yields an error otherwise.

if prod(size(b))==1
   c = times(a,1./b);
else
   error('Cannot use / for adiff objects when denominator is nonscalar');
end


text:
function experiments = stack_vot2019()

set_global_variable('bundle', 'http://data.votchallenge.net/vot2019/main/description.json');

baseline.name = 'baseline';
baseline.converter = [];
baseline.type = 'supervised';
baseline.tags = {'camera_motion', 'illum_change', 'occlusion', 'size_change', 'motion_change', 'empty'};

baseline.parameters.repetitions = 15;
baseline.parameters.burnin = 10;
baseline.parameters.skip_initialize = 5;
baseline.parameters.failure_overlap = 0;

baseline.analysis = {'expected_overlap', 'ar', 'speed'};

unsupervised.name = 'unsupervised';
unsupervised.converter = [];
unsupervised.type = 'unsupervised';
unsupervised.tags = {'camera_motion', 'illum_change', 'occlusion', 'size_change', 'motion_change', 'empty'};

unsupervised.parameters.repetitions = 1;
unsupervised.parameters.burnin = 0;

unsupervised.analysis = {'overlap', 'speed'};

realtime.name = 'realtime';
realtime.converter = [];
realtime.type = 'realtime';
realtime.tags = {'camera_motion', 'illum_change', 'occlusion', 'size_change', 'motion_change', 'empty'};

realtime.parameters.repetitions = 1;
realtime.parameters.default_fps = 20;
realtime.parameters.grace = 3;
realtime.parameters.burnin = 10;
realtime.parameters.override_fps = true;
realtime.parameters.skip_initialize = 5;
realtime.parameters.realtime_type = 'real';

realtime.analysis = {'expected_overlap'};

experiments = {baseline, unsupervised, realtime};

end

author: votchallenge

repo: toolkit-legacy


text:
%--- help for ts/pdecomp ---
%
%  Parametric decomposition into trend, seasonal and irregular components
%
%  ::
%
%    out=pdecomp(y)
%    out=pdecomp(y,doLog)
%    out=pdecomp(y,doLog,dorder)
%
%  Args:
%
%     y (ts): time series to decompose
%     doLog (true | {false}): if log, do a multiplicative decomposition
%       otherwise the decomposition is additive
%     dorder (integer | {2}): detrending order
%
%  Returns:
%     :
%
%     - **out** [struct] :
%
%       - **trend** [ts] : estimated trend
%       - **sc**    [ts] : estimated seasonal component
%       - **sa**    [ts] : seasonally adjusted data
%       - **ic**    [ts] : estimated irregular component
%
%  Note:
%
%     If there are many variables and the variables are named, the first level
%     of the structure will be the names of the different variables.
%
%  See also:
%     - :func:`npdecomp <ts.npdecomp>`
%
%

author: jmaih


text:
function x3 = mltply ( xx, x2, npl )

%*****************************************************************************80
%
%% MLTPLY multiplies two Chebyshev series.
%
%  Discussion:
%
%    This routine multiplies two given Chebyshev series, XX and X2,
%    to produce an output Chebyshev series, X3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Roger Broucke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Broucke,
%    Algorithm 446:
%    Ten Subroutines for the Manipulation of Chebyshev Series,
%    Communications of the ACM,
%    October 1973, Volume 16, Number 4, pages 254-256.
%
%  Parameters:
%
%    Input, real XX(NPL), the first Chebyshev series.
%
%    Input, real X2(NPL), the second Chebyshev series.
%
%    Input, integer NPL, the number of terms in the
%    Chebyshev series.
%
%    Output, real X3(NPL), the Chebyshev series of the
%    product.
%
  x3(1:npl) = 0.0;

  for k = 1 : npl
    ex = 0.0;
    mm = npl - k + 1;
    for m = 1 : mm
      l = m + k - 1;
      ex = ex + xx(m) * x2(l) + xx(l) * x2(m);
    end
    x3(k) = 0.5 * ex;
  end

  x3(1) = x3(1) - 0.5 * xx(1) * x2(1);

  for k = 3 : npl
    ex = 0.0;
    mm = k - 1;
    for m = 2 : mm
      l = k - m + 1;
      ex = ex + xx(m) * x2(l);
    end
    x3(k) = 0.5 * ex + x3(k);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function r8_min_test ( )

%*****************************************************************************80
%
%% R8_MIN_TEST tests R8_MIN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_MIN_TEST\n' );
  fprintf ( 1, '  R8_MIN returns the minimum of two R8''s.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '       A       B      C=R8_MIN(A,B)\n' );
  fprintf ( 1, '\n' );

  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  for i = 1 : 10
    [ a, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
    [ b, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
    c = r8_min ( a, b );
    fprintf ( 1, '  %8.4f  %8.4f  %8.4f\n', a, b, c );
  end

  return
end


text:
function cl = hlp_bm_updateLights(BMopts,hlights)
% update lighting for a brainmovie
% returns handles to camlight objects

if nargin<2
    hlights = []; end

if isempty(hlights) || ~ishandle(hlights(1))
    lighting(BMopts.vars.hbrainax,BMopts.facelighting);

    % delete existing lights
    hlights = findobj(BMopts.vars.hbrainax,'type','light');
    delete(hlights)
    hlights = [];
end

if ~isempty(BMopts.theme)
    cl=feval(BMopts.theme.lightingfcn,BMopts.vars.hbrainax,BMopts.theme.name,hlights);
    set(cl,'tag','brain_camlight','parent',BMopts.vars.hbrainax);
else
    cl(1) = camlight(BMopts.vars.hbrainax,'left');
    cl(2) = camlight(BMopts.vars.hbrainax,'right');
    set(cl,'tag','brain_camlight','parent',BMopts.vars.hbrainax);
end

repo: REST


text:
%Figure 1.1       Feedback Control of Dynamic Systems, 6e
%                      Franklin, Powell, Emami

% script to generate figure 1.1
clf
sim('fig1_01bsim',24)
plot(roomtemp(:,1),roomtemp(:,2))
axis([0 24 -1 21]);
title('Figure 1.1(b) Plot of room temperature and thermostat setting')
hold on
plot(heatin(:,1),heatin(:,2))
xlabel('Time (hours)');
ylabel('Temperature (degrees F)');
gtext('Room Temperature')
gtext(' Furnace on')
gtext(' Furnace off')
nicegrid;


text:
function test_suite = test_midwt
  disp("midwt")
  test_midwt_1D
  test_midwt_2D

function test_midwt_1D
       x = makesig('LinChirp',8);
       h = daubcqf(4,'min');
       L = 2;
       [y,L] = mdwt(x,h,L);
       [x_new,L] = midwt(y,h,L);
assertVectorsAlmostEqual(x, x_new,'relative',0.0001);

function test_midwt_2D
       load ../lena512;
       x = lena512;
       h = daubcqf(6);
       [y,L] = mdwt(x,h);
       [x_new,L] = midwt(y,h);
assertEqual(L,9);
assertVectorsAlmostEqual(x, x_new,'relative',0.0001);

author: ricedsp


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

    properties
        srcTank LaunchVehicleTank
        tgtTank LaunchVehicleTank
        initFlowRate(1,1) double = 0;
    end

    properties(Dependent)
        lvdData
    end

    methods
        function obj = TankToTankConnection(srcTank, tgtTank)
            if(nargin > 0)
                obj.srcTank = srcTank;
                obj.tgtTank = tgtTank;
            end
        end

        function lvdData = get.lvdData(obj)
            if(not(isempty(obj.srcTank)))
                lvdData = obj.srcTank.stage.launchVehicle.lvdData;
            else
                lvdData = obj.tgtTank.stage.launchVehicle.lvdData;
            end
        end

        function nameStr = getName(obj)
            if(isempty(obj.srcTank))
                srcNameStr = '<External>';
            else
                srcNameStr = obj.srcTank.name;
            end

            if(isempty(obj.tgtTank))
                tgtNameStr = '<External>';
            else
                tgtNameStr = obj.tgtTank.name;
            end

            nameStr = sprintf('%s to %s (%.3f mT/s)', srcNameStr, tgtNameStr, obj.initFlowRate);
        end

        function tf = isInUse(obj)
            tf = obj.lvdData.usesTankToTankConn(obj);
        end
    end

    methods(Static)
        function tankMassDots = getTankMassFlowRatesFromTankToTankConnections(tankStates, tankStatesMasses, t2tConnStates)
            activeT2TConnStates = t2tConnStates([t2tConnStates.active]);

            tankMassDots = zeros(size(tankStates));
            tankMassDots = tankMassDots(:);

            if(isempty(activeT2TConnStates))
                return;
            end

            tanks = [tankStates.tank];

            for(i=1:length(activeT2TConnStates))
                t2tConnState = activeT2TConnStates(i);

                srcTank = t2tConnState.conn.srcTank;
                if(not(isempty(srcTank)))
                    srcBool = tanks == srcTank;
                    if(any(srcBool))
                        srcTankState = tankStates(srcBool);
                        srcStageState = srcTankState.stageState.active;
                        srcIsTank = true;
                    else
                        srcIsTank = true;
                        srcStageState = false;
                    end
                else
                    srcStageState = true;
                    srcIsTank = false;
                end

                tgtTank = t2tConnState.conn.tgtTank;
                if(not(isempty(tgtTank)))
                    tgtBool = tanks == tgtTank;
                    if(any(tgtBool))
                        tgtTankState = tankStates(tgtBool);
                        tgtStageState = tgtTankState.stageState.active;
                        tgtIsTank = true;
                    else
                        tgtStageState = false;
                        tgtIsTank = true;
                    end
                else
                    tgtStageState = true;
                    tgtIsTank = false;
                end

                if(srcStageState && tgtStageState)
                    flowRate = t2tConnState.flowRate;

                    if(srcIsTank && tgtIsTank)
                        if(tankStatesMasses(srcBool) > 0)
                            tankMassDots(srcBool) = tankMassDots(srcBool) - flowRate;
                            tankMassDots(tgtBool) = tankMassDots(tgtBool) + flowRate;
                        end

                    elseif(srcIsTank && tgtIsTank == false)
                        if(tankStatesMasses(srcBool) > 0)
                            tankMassDots(srcBool) = tankMassDots(srcBool) - flowRate;
                        end

                    elseif(srcIsTank == false && tgtIsTank)
                        tankMassDots(tgtBool) = tankMassDots(tgtBool) + flowRate;
                    end
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function C4 = C4f(epsi, C4x)
%C4F  Evaluate C_4
%
%   C4 = C4F(EPSI, C4X) evaluates C_{4,l} in the expansion for the area
%   (Eq. (65) expressed in terms of n and epsi) using the coefficient
%   vector C4X.  K2 is a K x 1 array.  C4X is a 1 x 15 array.  C4 is a K x
%   5 array.

  nC4 = 6;
  nC4x = size(C4x, 2);
  j = nC4x;
  C4 = zeros(length(epsi), nC4);
  for k = nC4 : -1 : 1
    t = C4(:, k);
    for i = nC4 - k : -1 : 0
      t = epsi .* t + C4x(j);
      j = j - 1;
    end
    C4(:, k) = t;
  end
  mult = ones(length(epsi), 1);
  for k = 2 : nC4
    mult = mult .* epsi;
    C4(:, k) = C4(:, k) .* mult;
  end
end


text:
function [ctcell,ct_xmesh,ct_ymesh,ct_zmesh] = dicomrt_loadct(filename,sorting)
% dicomrt_loadct(filename,sorting)
%
% Read dicom CT for a case using MATLAB native function dicomread. Light Version.
%
% filename contains a list of CT slices to import
%
% CT are stored in a single 3D matrix: case_study
% x-y-z coordinates of the center of the ct-pixel are stored in ct_xmesh, ct_ymesh and ct_zmesh.
%
% NOTE:  CT numbers range from -1000 for air to 1000 for bone with that for water set to 0.
% CT numbers normalized in this manner are called Hounfield numbers or units (HU):
%
% HU = ((mu_tissue-mu_water)/mu_water)*1000
%
% Following DICOM RT standard HU = m*(SV)+b where m is RescaleSlope, b RescaleIntercept
% and SV are pixel (stored) values.
%
% Often CT scale is shifted so that HU(water)=1000 (=CToffset) instead of 0.
%
% NOTE: as opposed to dicomrt_loadct ct_xmesh, ct_ymesh and ct_zmesh are vectors and not
% matrices. This allow to run this functions also in "low" memory pcs.
%
% See also dicomrt_loaddose, dicomrt_getPatientPosition
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Sort CT slices first
%
% load ct slices info

% Check number of argument and set-up some parameters and variables
error(nargchk(1,2,nargin))

if exist('sorting')==0
    sorting=1; % perform sorting check
end

% Retrieve the number of images
nlines=dicomrt_nASCIIlines(filename);

if nlines>1
    disp('Loading CT slice information ...');
    [filelist,xlocation,ylocation,zlocation,modifUID,user_request] = dicomrt_loadctlist(filename);
    if user_request==1;
        return
    end
    disp('CT slice information loaded.');

    if sorting==1
        % sort ct slices
        disp('Sorting CT slices ...');
        [modifZ,user_request]=dicomrt_sortct(filelist,xlocation,ylocation,zlocation,filename);
        if user_request==1;
            return
        end
        disp('CT slices sorted.');
    else
        modifZ=0;
        modifUID=0;
    end

    if modifZ~=0 | modifUID~=0 % a modification to the filelist was made by dicomrt_sortct
        filename_sorted=[filename,'.sort.txt'];
    else % no modification have been made
        filename_sorted=filename;
    end
else
    disp('Loading one image. Zmesh will have no sense. Slice position can be loaded using dicominfo.');
    filename_sorted=filename;
end

% Now get CT images and create 3D Volume

% Set parameters
CToffset=1000;
nct=0;

% Define cell array to store info
case_study_info=cell(1);

%loop until the end-of-file is reached and build 3D CT matrix
disp('Loading ct volume ...');

% Progress bar
h = waitbar(0,['Loading progress:']);
set(h,'Name','dicomrt_loadct: loading CT images');

fid=fopen(filename_sorted);

while (feof(fid)~=1);
    clear info_temp temp
    nct=nct+1; % counting
    nctcheck=nct; % check for eof
    ct_file_location{1,nct}=fgetl(fid);
    if isnumeric(ct_file_location{1,nct}), nct=nct-1; break, end %end of line reached

    temp=dicomread(deblank(ct_file_location{1,nct}));

    dictFlg = checkDictUse;
    if dictFlg
        if isdeployed
            info_temp = dicominfo(deblank(ct_file_location{1,nct}),...
                'dictionary', fullfile(getCERRPath,'bin','ES - IPT4.1CompatibleDictionary.mat'));
        else
            info_temp = dicominfo(deblank(ct_file_location{1,nct}),...
                'dictionary', 'ES - IPT4.1CompatibleDictionary.mat');
        end
    else
        info_temp=dicominfo(deblank(ct_file_location{1,nct}));
    end

    zmesh(nct)=info_temp.ImagePositionPatient(3);

    if isfield(info_temp,'RescaleSlope')~=0 | isfield(info_temp,'RescaleIntercept')~=0
        temp=double(temp)*info_temp.RescaleSlope+info_temp.RescaleIntercept+CToffset;
    else
        warning('dicomrt_loadct: no DICOM Rescale data were found. Assuming RescaleSlope = 1, RescaleIntercept = 0 and CToffset = 1000');
        temp=double(temp);
    end

    case_study_info{nct}=info_temp;
    case_study(:,:,nct)=uint16(temp);
    waitbar(nct/nlines,h);
end

% Make zmesh a column vector
zmesh=dicomrt_makevertical(zmesh);

fclose(fid);

if nctcheck~=nct;
    warning('dicomrt_loadct: End of file was prematurely reached. Check the expected dimensions of your data. It may be OK to continue');
end

[PatientPositionCODE]=dicomrt_getPatientPosition(info_temp);

% PatientOrientation     ImageOrientationPatient
%
% HFS                    (1,0,0) (0,1,0)
% FFS                    (1,0,0) (0,1,0)
% HFP                    (-1,0,0) (0,-1,0)
% FFP                    (-1,0,0) (0,-1,0)
%
%
if PatientPositionCODE == 1 | PatientPositionCODE == 2
    min_x=info_temp.ImagePositionPatient(1);
    pixel_spacing_x=info_temp.PixelSpacing(1);

    min_y=info_temp.ImagePositionPatient(2);
    pixel_spacing_y=info_temp.PixelSpacing(2);

    [xmesh] = dicomrt_create1dmesh(min_x,pixel_spacing_x,size(temp,2),0);
    [ymesh] = dicomrt_create1dmesh(min_y,pixel_spacing_y,size(temp,1),0);

else
    max_x=info_temp.ImagePositionPatient(1);
    pixel_spacing_x=info_temp.PixelSpacing(1);

    max_y=info_temp.ImagePositionPatient(2);
    pixel_spacing_y=info_temp.PixelSpacing(2);

    [xmesh] = dicomrt_create1dmesh(max_x,pixel_spacing_x,size(temp,1),1);
    [ymesh] = dicomrt_create1dmesh(max_y,pixel_spacing_y,size(temp,2),1);

end

ct_zmesh=dicomrt_mmdigit(zmesh*0.1,7);
ct_ymesh=dicomrt_mmdigit(ymesh*0.1,7);
ct_xmesh=dicomrt_mmdigit(xmesh*0.1,7);

disp('Loading complete ...');

% 3D CT matrix and mesh matrix imported
%
% Some CT images info have private tags or fragmented info.
% If so, go into manual mode.
%

% Check support for current Patient Position
if PatientPositionCODE == 1 % supported Patient Position: HFS
    disp('dicomrt_loadct: Patient Position is Head First Supine (HFS).');
elseif PatientPositionCODE == 2 % supported Patient Position: FFS
    disp('dicomrt_loadct: Patient Position is Feet First Supine (FFS).');
elseif PatientPositionCODE == 3 % unsupported Patient Position: HFP
    disp('dicomrt_loadct: Patient Position is Head First Prone (HFP).');
elseif PatientPositionCODE == 4 % unsupported Patient Position: FFP
    disp('dicomrt_loadct: Patient Position is Feet First Prone (FFP).');
else
    warning('Unable to determine Patient Position');
    PatientPosition=input('dicomrt_loadct: Please specify Patient Position: HFS(default),FFS,HFP,FFP: ','s');
    if isempty(PatientPosition)==1
        PatientPosition='HFS';
    end

    if strcmpi(PatientPosition, 'HFS')
        PatientPositionCODE = 1;
    elseif strcmpi(PatientPosition, 'FFS')
        PatientPositionCODE = 2;
    elseif strcmpi(PatientPosition, 'HFP')
        PatientPositionCODE = 3;
    elseif strcmpi(PatientPosition, 'FFP')
        PatientPositionCODE = 4;
    end

end

% Create cell array
ctcell=cell(3,1);
ctcell{1,1}=case_study_info;
ctcell{2,1}=case_study;
ctcell{3,1}=[];

% Close progress bar
close(h);
pack

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function test_rb_inverse_dynamics()

path_to_urdf = 'ur10e.urdf';
ur10 = parse_urdf(path_to_urdf);

rbt = importrobot('ur10e.urdf');
rbt.DataFormat = 'column';
rbt.Gravity = [0 0 -9.81];

no_iter = 100;
for i = 1:no_iter
    q = -2*pi + 4*pi*rand(6,1);
    q_d = zeros(6,1);
    q_2d = zeros(6,1);

    Ylgr = standard_regressor_UR10E(q,q_d,q_2d);

    tau_matlab = inverseDynamics(rbt,q,q_d,q_2d);
    tau_reg = Ylgr*reshape(ur10.pi,[60,1]);
    tau_manip = M_mtrx_fcn(q, ur10.pi(:))*q_2d + ...
                C_mtrx_fcn(q, q_d, ur10.pi(:))*q_d + ...
                G_vctr_fcn(q, ur10.pi(:));

%   verifying if regressor is computed correctly
    assert(norm(tau_matlab - tau_reg) < 1e-8);
    assert(norm(tau_matlab - tau_manip) < 1e-8);
end

fprintf("Rigid Body Inverse Dynamics Test - OK!\n");

author: shamilmamedov


text:
function [u,s,U_r,U_s,U_pk,U_pd,U_l]  = projAhmPntIntoPinHoleOnRob(Rf, Sf, Spk, Spd, l)

% PROJAHMPNTINTOPINHOLEONROB Project Ahm pnt into pinhole on robot.
%    [U,S] = PROJAHMPNTINTOPINHOLEONROB(RF, SF, SPK, SPD, L) projects 3D
%    anchored homogeneous points into a pin-hole camera mounted on a robot,
%    providing also the non-measurable depth. The input parameters are:
%       RF : robot frame
%       SF : pin-hole sensor frame in robot
%       SPK: pin-hole intrinsic parameters [u0 v0 au av]'
%       SPD: radial distortion parameters [K2 K4 K6 ...]'
%       L  : 3D anchored homog. point [x y z vx vy vz rho]'
%    The output parameters are:
%       U  : 2D pixel [u v]'
%       S  : non-measurable depth
%
%    The function accepts an ahm points matrix L = [L1 ... Ln] as input.
%    In this case, it returns a pixels matrix U = [U1 ... Un] and a depths
%    row-vector S = [S1 ... Sn].
%
%    [U,S,U_R,U_S,U_K,U_D,U_L] = ... gives also the jacobians of the
%    observation U wrt all input parameters. Note that this only works for
%    single points.
%
%    See also PINHOLE, TOFRAME, PROJEUCPNTINTOPINHOLEONROB.

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

if nargout <= 2  % No Jacobians requested

    p     = ahm2euc(l);
    [u,s] = projEucPntIntoPinHoleOnRob(Rf, Sf, Spk, Spd, p);

else            % Jacobians requested

    if size(l,2) == 1

        % function calls
        [p,P_l]                     = ahm2euc(l);
        [u,s,U_r,U_s,U_pk,U_pd,U_p] = projEucPntIntoPinHoleOnRob(Rf, Sf, Spk, Spd, p);

        % chain rule
        U_l = U_p*P_l;

    else
        error('??? Jacobians not available for multiple AHM points.')

    end

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [nyquist_number,Bezugsfrequenz,Signal_amp] = Fast_Harmonics(x,fs,duration)
xx = 0;
if xx == 1
    % Sampling
    fs = 1000;     % Sampling rate [Hz]
    Ts = 1/fs;     % Sampling period [s]
    fNy = fs / 2;  % Nyquist frequency [Hz]
    duration = 10; % Duration [s]
    t = 0 : Ts : duration-Ts; % Time vector
    noSamples = length(t);    % Number of samples
    % Original signal
    x = 220.*sin(2 .* pi .* 50 .* t);
    % Harmonics
    x1 = 100.*sin(2 .* pi .* 100 .* t);
    x2 = 100.*sin(2 .* pi .* 200 .* t);
    x3 = 100.*sin(2 .* pi .* 300 .* t);
    % Contaminated signal
    xn = x + x1 + x2 + x3;
    % Frequency analysis
    f = 0 : fs/noSamples : fs - fs/noSamples; % Frequency vector
    % FFT
    x_fft = abs(fft(x));
    xn_fft = abs(fft(xn));
    % Plot
    figure(1);
    subplot(2,2,1);
    plot(t, x);
    subplot(2,2,2);
    plot(t, xn);
    subplot(2,2,3);
    plot(f,x_fft);
    xlim([0 fNy]);
    subplot(2,2,4);
    plot(f,xn_fft);
    xlim([0 fNy]);
else
%     x is the sampled data
%     m = length(x) is the window length (number of samples)
%     fs is the samples/unit time
%     dt = 1/fs is Time increment per sample
%     t = (0:m-1)/fs Time range for data
%     y = fft(x,n) is the Discrete Fourier Transform (DFT)
%     abs(y) is the Amplitude of the DFT
%     (abs(y).^2)/n is power of the DFT
%     fs/n Frequency increment
%     f = (0:n-1)*(fs/n) is frequency range
%     fs/2 is Nyquist frequency
%     n = pow2(nextpowe2(m)) is the transfrom length
    Ts = 1/fs;     % Sampling period [s]
    fNy = fs / 2;  % Nyquist frequency [Hz]
    t = 0 : Ts : duration-Ts; % Time vector
    noSamples = length(t);    % Number of samples
    % Original signal
    % Frequency analysis
    f = 0 : fs/noSamples : fs - fs/noSamples; % Frequency vector
    % FFT
    x_fft = abs(fft(x));
    xx = x';
    X = fft(xx(2:end));
    if rem(length(X),2)==0
        nyquist_number = length(X)/2;
    else
        nyquist_number = (length(X)+1)/2;
    end
    Bezugsfrequenz = 1/(t(end)-t(1));
    Signal_amp        = zeros(1,nyquist_number+1);
    Signal_phi        = zeros(1,nyquist_number+1);
    Signal_amp(1)     = X(1)/length(X);									                % Gleichanteil des Signals [unit]
    Signal_amp(2:end-1) = 2*abs(X(2:nyquist_number))/length(X);
    Signal_amp(end) = abs(X(nyquist_number+1))/length(X);
    %Signal_amp(2:end) = [2*abs(X(2:nyquist_number)) abs(X(nyquist_number+1))]/length(X);   % Amplituden der Harmonischen des Signals [unit]
    Signal_phi(1)     = 0;
    Signal_phi(2:end) = angle(X(2:nyquist_number+1));		                                % Phasenwinkel der Harmonischen des Signals [rad]

    %xn_fft = abs(fft(xn));
    % Plot
end


text:
function [Fmax, iMax] = bst_max(F, dim)
% BST_MAX: Get the maximum in magnitude, but with the correct sign.
%
% USAGE:  [Fmax, iMax] = bst_max(F, dim=1)

% @=============================================================================
% 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, 2014

% Parse inputs
if (nargin < 2) || isempty(dim)
    dim = [];
end

% Get the sign of each maximum
if isempty(dim)
    [Fmax, iMax] = max(abs(F(:)));
    Fmax = sign(F(iMax)) .* Fmax;
elseif (dim <= 5)
    % Permute with first dimension
    if (dim > 1)
        permdim = [dim, 2:ndims(F)];
        permdim(dim) = 1;
        F = permute(F, permdim);
    end
    % Reshape to ensure that the matrix has 2 dimensions
    nd = ndims(F);
    if (nd > 2)
        oldSize = [size(F,1), size(F,2), size(F,3), size(F,4), size(F,5)];
        F = reshape(F, size(F,1), []);
    end
    % Get maximum absolute values
    [Fmax, iMax] = max(abs(F), [], 1);
    % Build indices of the values to read
    iF = sub2ind(size(F), iMax, 1:size(F,2));
    Fmax = sign(F(iF)) .* Fmax;
    % Restore initial shape
    if (nd > 2)
        Fmax = reshape(Fmax, 1, oldSize(2), oldSize(3), oldSize(4), oldSize(5));
        iMax = reshape(iMax, 1, oldSize(2), oldSize(3), oldSize(4), oldSize(5));
    end
    % Restore initial dimensions
    if (dim > 1)
        Fmax = permute(Fmax, permdim);
        iMax = permute(iMax, permdim);
    end
else
    error('Not supported yet.');
end

author: brainstorm-tools


text:
function test_bug2232

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY gifti ft_read_headshape ft_read_atlas

ft_hastoolbox('gifti', 1);

datadir = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2232');
d       = dir(datadir);
subdir  = {d.name}';
sel     = cellfun(@numel, subdir)>2;
subdir  = subdir(sel);

% check whether we can read *.surf.gii
for k = 1:numel(subdir)
  d = dir(fullfile(datadir,subdir{k},'*surf*'));
  if ~isempty(d)
    for m = 1:numel(d)
      fprintf('reading %s\n',d(m).name);
      bnd = ft_read_headshape(fullfile(datadir, subdir{k}, d(m).name));
      assert(isfield(bnd, 'pos') && isfield(bnd, 'tri'));

      hemi = strfind(d(m).name, '.L.');
      if ~isempty(hemi)
        hemi = '.L.';
      else
        hemi = '.R.';
      end
      dlabel = dir(fullfile(datadir,subdir{k},['*',hemi,'*label*']));
      if ~isempty(dlabel)
        for q = 1:numel(dlabel)
          fprintf('reading parcellation based on %s\n', dlabel(q).name);
          bnd = ft_read_atlas({fullfile(datadir,subdir{k},dlabel(q).name) fullfile(datadir,subdir{k},d(m).name)});
        end
      end
    end
  end

  % check whether there are subdirs
  d   = dir(fullfile(datadir,subdir{k}));
  sel = cellfun(@numel, {d.name}')>2;
  d   = d(sel);

  isdir     = find([d.isdir]);
  d         = d(isdir);
  if ~isempty(d)
    subsubdir = {d.name}';
    for p = 1:numel(subsubdir)
      d = dir(fullfile(datadir,subdir{k},subsubdir{p},'*surf*'));
      if ~isempty(d)
        for m = 1:numel(d)
          fprintf('reading %s\n',d(m).name);
          bnd = ft_read_headshape(fullfile(datadir, subdir{k}, subsubdir{p}, d(m).name));
          assert(isfield(bnd, 'pos') && isfield(bnd, 'tri'));

          hemi = strfind(d(m).name, '.L.');
          if ~isempty(hemi)
            hemi = '.L.';
          else
            hemi = '.R.';
          end
          dlabel = dir(fullfile(datadir,subdir{k},['*',hemi,'*label*']));
          if ~isempty(dlabel)
            for q = 1:numel(dlabel)
              fprintf('reading parcellation based on %s\n', dlabel(q).name);
              bnd = ft_read_atlas({fullfile(datadir,subdir{k},subsubdir{p},dlabel(q).name) fullfile(datadir,subdir{k},subsubdir{p},d(m).name)});
            end
          end

        end
      end
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%DNADamage
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 4/19/2011
classdef DNADamage
    methods (Static = true)
        function run(sim, fileName)
            import edu.stanford.covert.cell.sim.analysis.DNADamage;
            import edu.stanford.covert.cell.sim.util.PrintUtil;

            %excel file
            [content, colLabels, indentation] = DNADamage.calcExpectedRates(sim);
            if nargin == 1
                PrintUtil.printToStdIO(content, colLabels, struct('indentation', indentation));
            else
                PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'State', struct('indentation', indentation));
            end
        end

        function [content, colLabels, indentation] = calcExpectedRates(sim)
            d = sim.process('DNADamage');

            reactionVulnerableMotifs = d.reactionVulnerableMotifs;
            tfs = ismember(d.reactionVulnerableMotifTypes, {'damagedBases','damagedSugarPhosphates','intrastrandCrossLinks'});
            reactionVulnerableMotifs(tfs) = ...
                d.metabolite.wholeCellModelIDs(cell2mat(d.reactionVulnerableMotifs(tfs)));
            tfs = ismember(d.reactionVulnerableMotifTypes, {'abasicSites','gapSites','hollidayJunctions', 'strandBreaks'});
            reactionVulnerableMotifs(tfs) = cell(sum(tfs), 1);

            rates = d.calcExpectedReactionRates();

            reactionRadiation = cell(size(d.reactionRadiation));
            reactionRadiationLevels = zeros(size(d.reactionRadiation));
            reactionRadiation(d.reactionRadiation~=0) = ...
                d.substrateWholeCellModelIDs(d.reactionRadiation(d.reactionRadiation~=0));
            reactionRadiationLevels(d.reactionRadiation~=0) = ...
                d.substrates(d.reactionRadiation(d.reactionRadiation~=0));

            colLabels = {'ID', 'Name', ...
                'Damage Type', 'Damage Motif', 'No. Vulnerable Site', ...
                'Specific Rate (eg 1/s/nt/gy)', 'Rate (1/s/cell)', 'Rate (1/cell cycle)', ...
                'Radiation ID', 'Radiation Value'};

            content = [...
                d.reactionWholeCellModelIDs d.reactionNames ...
                d.reactionVulnerableMotifTypes reactionVulnerableMotifs num2cell(d.calcNumberVulnerableSites()) ...
                num2cell(d.reactionBounds(:, 2)) num2cell(rates) num2cell(rates * ((9.0 + 2) * 3600)) ...
                reactionRadiation num2cell(reactionRadiationLevels) ...
                ];

            indentation = zeros(size(content, 1), 1);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%Test script to run SVP for matrix completion using uniformly sampled random matrices
%Written by: Prateek Jain (pjain@cs.utexas.edu) and Raghu Meka (meka@cs.utexas.edu)
%Last updated: October, 2009

m = 5000;
n = 2000;
k = 10; %rank of the optimal matrix
p = 0.2; %fraction of known entries
params.tol=1e-3;
params.vtol=1e-4;
params.mxitr=100;
params.verbosity=1;

M = randn(m,k)*randn(k,n); %Generate the underlying matrix
random_permutation=randperm(m*n);
Omega=random_permutation(1:p*m*n)'; %Generate samples with uniform density p
Omega=sort(Omega);

fprintf('Running SVP with m=%d, n=%d, p=%f...\n',m,n,p);
t=cputime;
[U,S,V,num_iter] = svp(Omega, M(Omega), m, n, k, params); %Run SVP
t = cputime-t;
fprintf('SVP finished with RMSE=%f in time t=%f\n',norm(M-U*diag(S)*V','fro')/sqrt(length(Omega)),t);

author: andrewssobral

repo: lrslibrary


text:
function logPr = calcMargLogPrData_ARMatrixNormalInvWishart( Xstats, PP )

[D DR K] = size( Xstats.XY );
LOG_PI =  1.144729885849400;

logPr = -Inf( 1, K );
for jj = 1:K
    XX = Xstats.XX(:,:,jj);
    XY = Xstats.XY(:,:,jj);
    YY = Xstats.YY(:,:,jj);
    CC = PP.invAScaleMatrix;

    degFreeN = PP.degFree + Xstats.nObs(jj);

    ScaleMatrixN = PP.ScaleMatrix + XX - ( XY/(YY + CC) )*XY';
    %ScaleMatrixN = PP.ScaleMatrix + XX - XY*( (YY+CC)\(XY') );

    logPr(jj) = logMvGamma( 0.5*degFreeN, D ) - logMvGamma( 0.5*PP.degFree, D ) ...
                  + 0.5*PP.degFree*log( det( PP.ScaleMatrix ) )   ...
                  - 0.5*degFreeN  *log( det(   ScaleMatrixN ) ) ...
                  + 0.5*D*log( det( CC ) ) ...
                  - 0.5*D*log( det( CC + YY ) );
end

N = sum(Xstats.nObs);
logPr = sum( logPr ) + -0.5*N*D*LOG_PI;

author: michaelchughes

repo: NPBayesHMM


text:
function output = base64decode(input)
%BASE64DECODE Decode Base64 string to a byte array.
%
%    output = base64decode(input)
%
% The function takes a Base64 string INPUT and returns a uint8 array
% OUTPUT. JAVA must be running to use this function. The result is always
% given as a 1-by-N array, and doesn't retrieve the original dimensions.
%
% See also base64encode

error(nargchk(1, 1, nargin));
error(javachk('jvm'));
if ischar(input), input = uint8(input); end

output = typecast(org.apache.commons.codec.binary.Base64.decodeBase64(input), 'uint8')';

end


text:
% demo_example_lmgp_training - demonstration script file for LMGP model identification - training part.

%% Description
% Demo to present the training, i.e. identification, of LMGP model.
% Note:
% It can be used only with the Gaussian covariance function and
% with the white noise model (sum of covSEard and covNoise).

% See Also
% demo_example_LMGP_data.m, demo_example_LMGP_simulation.m, trainLMGP.m,
% gpSD00.m

%
% Changelog:
%
% 16.2.2015, Martin Stepancic:
%		 	-added compatibility with gpml>=3.1
%

clear;
global flag_LM_data_ident
close all;

%%%%% off-equlibrium data
mat_data=load('example_data_lmgp_oeq.mat');
ind_oeq = 5:2:length(mat_data.train_data.u);

utrain = mat_data.train_data.u(ind_oeq);
xtrain = mat_data.train_data.x(ind_oeq);
ytrain = mat_data.train_data.y(ind_oeq);

%%%%% identified local-model data
if(flag_LM_data_ident == 1)
    load example_data_lmgp_eq_ident
else
     load example_data_lmgp_eq_anal
end

Yeq = eq.Y;
Ueq = eq.U;
dfdy = eq.dfdy;
dfdu = eq.dfdu;

%***************************
% training data construction
% functional
input = [Yeq Ueq; xtrain utrain];
target = [Yeq; ytrain];
targetvar = NaN*ones(size(target));
% derivative

inputDer = [Yeq Ueq];
targetDer = [dfdy dfdu];

% LM variance

for ii=1:size(inputDer,1)
    if(flag_LM_data_ident == 1)
%         eq.lm{ii}.CovarianceMatrix
        derivevar(ii,:) = reshape(eq.lm{ii}.CovarianceMatrix,1,4);
    else
        derivevar(ii,:) = reshape(0.01*eye(2),1,4);
    end
end

% covariance function: SE + white noise
covfunc = {'covSum',{'covSEard','covNoise'}};

%************************************************************************
% training

% set the mean, covariance, inference method, and likelihood function. These
% are all unused and dummy, since the gpSD00 function has its own implementation.
inffunc= 'infExact';
meanfunc='meanZero';
covfunc= 'covSEard';
likfunc= 'likGauss';

% initialize hyperparameters:
hyp0.cov=ones(size(input,2)+1,1);
hyp0.lik=-1;

lag = 1;

hyp=trainLMGP([],inffunc,meanfunc,covfunc,likfunc,input,target,inputDer, targetDer, derivevar);

% validation on ident data
for ii=1:size(input,1)
    test = input(ii,:);
    [mug(ii) s2(ii)]=gpSD00(hyp,inffunc,meanfunc,covfunc,likfunc,input,target,targetvar, inputDer, targetDer, derivevar, test);
end

figure('Name', 'Validation on Identification Data');
plot(Ueq,Yeq,'r',utrain(:,1),ytrain,'or',[input(:,2)],mug,'k*',...
    [input(:,2)],mug+2*s2,'k.',[input(:,2)],mug-2*s2,'k.');

disp('Hyperparameters: ');
disp(num2str(exp(unwrap(hyp))));
disp(' ');
disp(['Number of local models: ', num2str(length(targetDer))]);
disp(['Number of points out of equilibrium: ', num2str(length(utrain))]);
disp(['Number of training vectors: ', num2str(length(target)+numel(targetDer))]);
disp(' ');

    save example_lmgp_trained hyp input target targetvar ...
        inputDer targetDer derivevar covfunc meanfunc inffunc likfunc lag;
return

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function value = r8_normal_01_sample ( )

%*****************************************************************************80
%
%% R8_NORMAL_01_SAMPLE returns a unit pseudonormal R8.
%
%  Discussion:
%
%    The standard normal probability distribution function (PDF) has
%    mean 0 and standard deviation 1.
%
%    The Box-Muller method is used, which is efficient, but
%    generates two values at a time.
%
%    Typically, we would use one value and save the other for the next call.
%    However, the fact that this function has saved memory makes it difficult
%    to correctly handle cases where we want to re-initialize the code,
%    or to run in parallel.  Therefore, we will instead use the first value
%    and DISCARD the second.
%
%    EFFICIENCY must defer to SIMPLICITY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a sample of the standard normal PDF.
%
  r1 = r8_uniform_01_sample ( );
  r2 = r8_uniform_01_sample ( );

  x = sqrt ( - 2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 );

  value = x;

  return
end


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, spv, vel, gnum, spp, press, gnump] = ...
              mp_solve_stokes_3d (problem_data, method_data)

warning ('geopdes:obsolete','Function MP_SOLVE_STOKES_3D is obsolete. Using MP_SOLVE_STOKES instead')

function [geometry, msh, spv, vel, gnum, spp, press, gnump] = ...
              mp_solve_stokes (problem_data, method_data)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%MDL_S4ABB2p8  Create kinematic model of ABB S4 2.8robot
%
% MDL_S4ABB2p8 is a script that creates the workspace variable s4 which
% describes the kinematic characteristics of an ABB S4 2.8 robot using
% standard DH conventions.
%
% Also defines the workspace vector:
%   q0   mastering position.
%
% Author::
%  Wynand Swart,
%  Mega Robots CC, P/O Box 8412, Pretoria, 0001, South Africa,
%  wynand.swart@gmail.com
%
% See also mdl_fanuc10l, mdl_m16, mdl_motormanHP6, mdl_irb140, mdl_puma560, SerialLink.

% MODEL: ABB, S4_2.8, S4 2.8m reach version, 6DOF, standard_DH

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.

%Cell: 073-1555-430
%30 Sep 2007
%S4 ABB 2.8 robot

%            theta    d      a    alpha
clear L
L(1) = Link([ 0      0.9    0.188  -pi/2   0]);
L(2) = Link([ 0      0      0.95    0      0]);
L(3) = Link([ 0      0      0.225  -pi/2   0]);
L(4) = Link([ 0      1.705  0       pi/2   0]);
L(5) = Link([ 0      0      0      -pi/2   0]);
L(6) = Link([ 0      0.2    0      -pi/2   0]);
%##########################################################
%Pose 0; At SYNCHRONISATION position
%##########################################################
q0 = [0     -pi/2         0       0      0     -pi/2];
s4=SerialLink(L, 'name', 'S4 ABB 2.8');
%##########################################################

author: petercorke

repo: robotics-toolbox-matlab


text:
% this class computes sentence-to-vec model, given word-to-vec model.
classdef Model < handle
    properties(Constant)
        MODE_AVG      = 0   % encode by average
        MODE_FISHER   = 1   % encode by fisher vectors
        MODE_BOW      = 2   % bag of words

        % improved fv options

        POWER_OFF           = 0   % no power transform
        POWER_SQRT_ALL      = 1   % sqrt all coordinates
        POWER_SQRT_GMM_ONLY = 2   % sqrt only the gmm coordinates
		POWER_EXPERIMENTAL  = 3

        NORMALIZE_OFF             = 0 % no normalization
        NORMALIZE_L2_ALL          = 1 % L2 normalization
        NORMALIZE_L1_L2_SELECTIVE = 2 % L2 for gmm coordinates, L1 for lmm
    end

    properties

        % improved fv options
        power_option;
        norm_option;

        % Word2vec model
        w2v;

        % GMM / LMM model
        mm_type
        mm_num_of_centers

        % will be used for GMM / LMM
        mm_means
        mm_covariances
        mm_priors

        % will be used for HGLMM
        hglmm_m
        hglmm_s
        hglmm_mu
        hglmm_sigma
        hglmm_b
        hglmm_prior

        % will be used for GG
        gg_p
        gg_mu
        gg_sigma

        % options paramter for gg_fv function
        gg_options

        % multiplication factor (workaround for numeric issues with LMM / HGLMM)
        mult_fctr;

        % transition matrix, for PCA / ICA
        trans_mat;

    end

    methods

        % constructor
        % mm_type: one of: 'gmm', 'lmm', 'hcg', 'hcl', 'hglmm'
        % lowercase: a flag indicating that our dictionary is in lowercase
        % (this is the case with GloVe)
        function this = Model(words_file_name, vectors_file_name, freq_file_name, mm_file_name, mm_type, trans_matrix, lowercase, power_option, norm_option)

            if ~exist('lowercase', 'var')
                lowercase = false;
            end

            if ~exist('power_option', 'var')
                power_option = this.POWER_SQRT_ALL;
            end

            if ~exist('norm_option', 'var')
                norm_option = this.NORMALIZE_L2_ALL;
            end

            this.gg_options = [];

            this.power_option = power_option;
            this.norm_option = norm_option;

            this.trans_mat = trans_matrix;

			% load word2vec
            this.w2v = Str2vec(words_file_name, vectors_file_name, freq_file_name, lowercase);

            if ~isempty(mm_file_name)

                this.mm_type = mm_type;

                lg(1, 'MM type: %s\n', this.mm_type);

                % read GMM/LMM file
                lg(1, 'reading %s file %s\n', this.mm_type, mm_file_name);
                load(mm_file_name);
                lg(1, 'done\n');

                % if mult_factor was saved in the mm file (it is relevnt
                % for LMM / HGLMM)
                if exist('mult_factor', 'var')
                    this.mult_fctr = mult_factor;
                else
                    switch this.mm_type
                        case {'lmm', 'hglmm'}
                            % default is 10 (backward compatibility)
                            this.mult_fctr = 10;
                        otherwise
                            this.mult_fctr = 1;
                    end
                end

                switch this.mm_type
                    case 'hglmm'
                        this.hglmm_m = m_out;
                        this.hglmm_s = s_out;
                        this.hglmm_mu = mu_out;
                        this.hglmm_sigma = sigma_out;
                        this.hglmm_b = b_out;
                        this.hglmm_prior = prior_out;

                        this.mm_num_of_centers = length(this.hglmm_prior);

                    case 'gg'
                        this.gg_p = gg_p;
                        this.gg_mu = gg_mu;
                        this.gg_sigma = gg_sigma;

                        this.mm_num_of_centers = size(gg_mu,2);

                    % gmm or lmm
                    otherwise
                        this.mm_means = n_means;
                        this.mm_covariances = n_covariances;
                        this.mm_priors = n_priors;

                        this.mm_num_of_centers = length(this.mm_priors);
                end
            else
                lg(1, 'no MM file was provided\n');
            end
        end

        function set_gg_options(this, gg_options)
            this.gg_options = gg_options;
        end

        function str = mode2str(this, mode)
            switch mode
                case this.MODE_AVG
                    str = 'AVG';
                case this.MODE_FISHER
                    str = 'FISHER';
                case this.MODE_BOW
                    str = 'BOW';
                otherwise
                    str = 'UNKNOWN';
            end
        end

        function d = encoding_size(this, mode)
            dim = size(this.trans_mat, 1);

            % trans_mat is empty
            if dim == 0
                dim = this.w2v.dimension;
            end

            switch mode
                case this.MODE_AVG
                    d = dim;
                case this.MODE_FISHER
                    switch this.mm_type
                        case 'gg'
                            d = 2 * this.mm_num_of_centers * dim;
                        otherwise
                            d = 2 * this.mm_num_of_centers * dim;
                    end
                case this.MODE_BOW
                    d = this.w2v.num_strings;
                otherwise
                    d = 0;
            end
        end

        % encodes a sentence using fisher vectors
        function encoding = encode_sen_fv(this, sentence, num_infreq_words, sum_pair_mode)

			if ~exist('sum_pair_mode', 'var')
				sum_pair_mode = false;
            end

            X = this.w2v.sen_to_word_vecs(sentence, num_infreq_words);
            if isempty(X)
                encoding = [];
                return;
            end

            % encode the sentence using fisher vectors

            if sum_pair_mode
                Z = zeros(size(X,1), size(X,2) - 1);
                for i = 1:size(Z,2)
                    Z(:,i) = X(:,i) + X(:,i+1);
                end

                X = [X, Z];
            end

            if ~isempty(this.trans_mat)
                X = this.trans_mat * X;

                % in pca and ica, we normalized the vectors after
                % transformation
                X = normc(X);
            end

            switch this.mm_type
                case {'gmm', 'hcg'}

					switch this.power_option
                        case {this.POWER_OFF, this.POWER_EXPERIMENTAL}
                            error('need to add support for POWER_OFF, POWER_EXPERIMENTAL');
                    end

                    if this.norm_option == this.NORMALIZE_OFF
                        error('need to add support for NORMALIZE_OFF');
                    end

                    encoding = vl_fisher(X, this.mm_means, this.mm_covariances, this.mm_priors, 'SquareRoot', 'Normalized');

                case 'gg'

					switch this.power_option
                        case {this.POWER_OFF, this.POWER_EXPERIMENTAL}
                            error('need to add support for POWER_OFF, POWER_EXPERIMENTAL');
                    end

                    if this.norm_option == this.NORMALIZE_OFF
                        error('need to add support for NORMALIZE_OFF');
                    end

                    encoding = gg_fv(X, this.gg_p, this.gg_mu, this.gg_sigma, this.gg_options);

                case {'lmm', 'hcl'}
                    encoding = FVMex(double(this.mult_fctr * X'),double(this.mm_means), double(this.mm_covariances), double(this.mm_priors),1.0);

                    switch this.power_option
                        case this.POWER_SQRT_ALL
                            encoding = sign(encoding) .* sqrt(abs(encoding));
                        case this.POWER_EXPERIMENTAL
                            %encoding = sign(encoding) .* (abs(encoding).^experimental_power);
                            %encoding = sign(encoding) .* log(abs(encoding) + 1);
                            encoding = asinh(encoding);
                    end

                    switch this.norm_option
                        case this.NORMALIZE_L2_ALL
                            encoding = encoding / norm(encoding, 2);
                        case this.NORMALIZE_L1_L2_SELECTIVE
                            encoding = encoding / norm(encoding, 1);
                    end

                case 'hglmm'
                    numOfCpus = 8;

                    encoding = HybridFV(double(this.mult_fctr * X'), this.hglmm_m, this.hglmm_s, this.hglmm_mu, this.hglmm_sigma, this.hglmm_b, this.hglmm_prior, numOfCpus);

                    I_lmm = this.hglmm_b';
                    I_lmm = I_lmm(:);
                    I_lmm = [I_lmm; I_lmm];
                    I_gmm = ~I_lmm;

                    switch this.power_option
                        case this.POWER_OFF
                            I_power = false(length(encoding), 1);
                        case this.POWER_SQRT_ALL
                            I_power = true(length(encoding), 1);
                        case this.POWER_SQRT_GMM_ONLY
                            I_power = I_gmm;
                    end

                    switch this.norm_option
                        case this.NORMALIZE_OFF
                            I_norm2 = false(length(encoding), 1);
                            I_norm1 = false(length(encoding), 1);
                        case this.NORMALIZE_L2_ALL
                            I_norm2 = true(length(encoding), 1);
                            I_norm1 = false(length(encoding), 1);
                        case this.NORMALIZE_L1_L2_SELECTIVE
                            I_norm2 = I_gmm;
                            I_norm1 = I_lmm;
                    end

                    if this.power_option == this.POWER_EXPERIMENTAL
                        %encoding = sign(encoding) .* (abs(encoding).^experimental_power);
                        %encoding = sign(encoding) .* log(abs(encoding) + 1);
                        encoding = asinh(encoding);
                    else
                        encoding(I_power) = sign(encoding(I_power)) .* sqrt(abs(encoding(I_power)));
                    end

                    encoding(I_norm2) = encoding(I_norm2) / norm(encoding(I_norm2),2);
                    encoding(I_norm1) = encoding(I_norm1) / norm(encoding(I_norm1),1);

                otherwise
                    error('illegal mm_type');
            end
        end

        % encodes a sentence using average
        function encoding = encode_sen_avg(this, sentence, num_infreq_words)

            X = this.w2v.sen_to_word_vecs(sentence, num_infreq_words);
            if isempty(X)
                encoding = [];
                return;
            end

            if ~isempty(this.trans_mat)
                X = this.trans_mat * X;

                % in pca and ica, we normalized the vectors after
                % transformation
                X = normc(X);
            end

            % encode the sentence using average
            encoding = sum(X, 2) / size(X, 2);
            % normalize
            encoding = norma(encoding);
        end

        % bag of words (normalized)
        function bow = encode_sen_norm_bow(this, sentence, num_infreq_words)
            if ~isempty(num_infreq_words)
                error('encode_sen_norm_bow does not support num_infreq_words yet');
            end

            valid_words_in_sen = this.w2v.all_valid_words(sentence, false);
            num_valid_words_in_sen = length(valid_words_in_sen);

            bow = zeros(this.w2v.num_strings, 1);

            for j = 1:num_valid_words_in_sen
                word = valid_words_in_sen{j};
                word_idx = this.w2v.str2idx(word);
                bow(word_idx) = bow(word_idx) + 1;
            end

            % normalize
            bow = norma(bow);
        end

        % 'sentence' may be either a string, or a vector with words indices
        % (according to their order in the word2vec) of the valid words of
        % a sentence. for MODE_BOW, the words indices option is currenlty
        % not supported.
        function encoding = encode_sen(this, mode, sentence, num_infreq_words)

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

            switch mode
                case this.MODE_FISHER
                    encoding = this.encode_sen_fv(sentence, num_infreq_words);
                case this.MODE_AVG
                    encoding = this.encode_sen_avg(sentence, num_infreq_words);
                case this.MODE_BOW
                    if ~ischar(sentence)
                        error('for MODE_BOW, the words indices option is currenlty not supported');
                    end
                    encoding = this.encode_sen_norm_bow(sentence, num_infreq_words);
                otherwise
                    error('illegal mode');
            end
        end

        % words_I: a vector of indices of valid words (return by
        % all_valid_words with paramter word_indices==true).
        % this function encodes each word as an entire sentence and retures
        % the matrix of all encodings
        function words_encoding = encode_words_as_sentences(this, mode, words_I)
            num_words = length(words_I);

            words_encoding = zeros(this.encoding_size(mode), num_words);

            for i = 1:num_words
                words_encoding(:,i) = this.encode_sen(mode, words_I(i));
            end
        end

        % returns a matrix containing the encodings of the given sentences.
        % I_bad: indicators for sentences which could not be encoded since
        % none of their words is in the vocabulary. in such cases, the
        % vector of the sentence will be all zeros
        function [encoding, I_bad] = encode_sentences(this, mode, sentences)
            % if only one sentence - make it a list
            if ~iscell(sentences)
                sentences = {sentences};
            end

            num_sentences = length(sentences);

            encoding = zeros(this.encoding_size(mode), num_sentences);
            I_bad = false(1, num_sentences);

            lg(1, 'encoding %d sentences...\n', num_sentences);

            for i = 1:num_sentences
                v = this.encode_sen(mode, sentences{i});

                if isempty(v)
                    I_bad(i) = true;
                    lg(0, '\n** Sentence number %d could not be encoded since none of its words is in the vocabulary:\n', i);
                    lg(0, '%s\n', sentences{i});
                else
                    encoding(:,i) = v;
                end

                printmark(i, 100, num_sentences)
            end

            lg(1, 'done\n');
        end

    end

    methods(Static)

        function str = power2str(power_option)
            switch power_option
                case Model.POWER_OFF
                    str = 'POWER_OFF';
                case Model.POWER_SQRT_ALL
                    str = 'POWER_SQRT_ALL';
                case Model.POWER_SQRT_GMM_ONLY
                    str = 'POWER_SQRT_GMM_ONLY';
                case Model.POWER_EXPERIMENTAL
                    str = 'POWER_EXPERIMENTAL';
                otherwise
                    str = 'UNKNOWN';
            end
        end

        function str = normalize2str(normalize_option)
            switch normalize_option
                case Model.NORMALIZE_OFF
                    str = 'NORMALIZE_OFF';
                case Model.NORMALIZE_L2_ALL
                    str = 'NORMALIZE_L2_ALL';
                case Model.NORMALIZE_L1_L2_SELECTIVE
                    str = 'NORMALIZE_L1_L2_SELECTIVE';
                otherwise
                    str = 'UNKNOWN';
            end
        end

    end

end

author: layumi

repo: Image-Text-Embedding


text:
 function y = mtimes(ob, x)
%function y = mtimes(ob, x)
% y = G * x	or x = G' * y

%
% full forward projection
%
if ~ob.is_transpose

	% if needed, convert 2d array to concise column
	if ob.is_masked & size(x,1) ~= ob.dims(2);
		error 'not done'
%		if size(x,1) ~= sum(sum(ob.mask));
%			error 'size mismatch'
%		end
%		x = embed(x, ob.mask);
%		x = x(:);
	end

	if ob.apower ~= 1
		y = (ob.G).^(ob.apower) * x;
	else
		y = ob.G * x;
	end

%
% full backprojection
%
else
	if ob.apower ~= 1
		y = (ob.G .^ ob.apower)' * x;
	else
		y = (x' * ob.G)';
	end
%	if ob.is_masked
%		y = y(ob.mask,:);
%	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [breadth, labels] = imMeanBreadthEstimate(img, varargin)
% Estimate mean breadth of a binary structure.
%
%   BREADTH = imMeanBreadthEstimate(IMG)
%   The aim of this function is to be called by the "imMeanBreadthDensity"
%   function, for providing an estimate of mean breadth density within a
%   representative volume of interest.
%
%   BREADTH = imMeanBreadthEstimate(IMG, RESOL)
%   Also specifies the resolution of input image.
%
%   BREADTH = imMeanBreadthEstimate(IMG, RESOL, NDIRS)
%   Also specifies the number of directions used for computation. Can be
%   either 3 or 13 (the default).
%
%   Example
%
%   See also
%     imMeanBreadthEstimate, imSurfaceAreaEstimate
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2015-04-20,    using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015 INRAE - Cepia Software Platform.

%% Basic error checking

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    % identify the labels in image
    labels = imFindLabels(img);

    % allocate result array
    nLabels = length(labels);
    breadth = zeros(nLabels, 1);

    % compute bounding box of each label
    boxes = imBoundingBox(img);

    % Compute mean breadth of each label considered as binary image
    % The computation is performed on a subset of the image for reducing
    % memory footprint.
    for i = 1:nLabels
        label = labels(i);

        % convert bounding box to image extent, in x, y and z directions
        box = boxes(i,:);
        i0 = ceil(box([3 1 5]));
        i1 = floor(box([4 2 6]));

        % crop image of current label
        bin = img(i0(1):i1(1), i0(2):i1(2), i0(3):i1(3)) == label;
        breadth(i) = imMeanBreadthEstimate(bin, varargin{:});
    end

    return;
end

%% Process input arguments

% in case of binary image, compute only one label...
labels = 1;

% default number of directions
nDirs = 13;

% default image resolution
delta = [1 1 1];

% Process user input arguments
while ~isempty(varargin)
    var = varargin{1};
    if ~isnumeric(var)
        error('option should be numeric');
    end

    % option is either connectivity or resolution
    if isscalar(var)
        nDirs = var;
    else
        delta = var;
    end
    varargin(1) = [];
end

%% Compute mean breadth estimate using Look-up-Table

histo = imBinaryConfigHisto(img);

lut = imMeanBreadthLut(delta, nDirs);

breadth = sum(histo .* lut);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [mfg, mbg] = im_modes(im)
% IM_MODES  Estimate typical values (modes) of foreground/background
% intensities in a greyscale image.
%
% [MFG, MBG] = im_modes(IM)
%
%   IM is an N-array of grayscale intensity values.
%
%   MFG, MBG are scalars with the typical intensity values of foreground
%   (darker) and foreground (lighter) voxels. "Typical" means the mode of
%   the corresponding distributions.
%
%   The modes are estimated by computing histograms with up to 500 bins for
%   the intensity values, and smoothing the distribution until only 1 or 2
%   peaks are visible. In case of 1 peak, it is assume that the image
%   contains only background voxels, and MFG=NaN.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

DEBUG = 0;

% check arguments
narginchk(1, 1);
nargoutchk(0, 2);

% ignore intensity values = 0. Those are considered to be masked out
im2 = im(im > 0);

% if the input image is empty, or too small, we assume that we don't have
% enough information to estimate the background and foreground
if (nnz(im2) < 100)
    mbg = nan;
    mfg = nan;
    return
end

% number of bins to use. We assume that we need at least 10 samples per
% bin, but we don't need more than 500 bins in total. Histograms with more
% bins are slower to process and smooth
nbin = min(500, ceil(nnz(im2)/10));

% compute histogram of the intensity values
[fhist, xhist] = hist(im2, nbin);
fhist = fhist / sum(fhist);

% initial estimation of peaks in the histogram
[pks, loc] = findpeaks(fhist, 'minpeakheight', 0.5e-3);

% number of peaks found
npks = length(pks);

% smooth the histogram until we find just 1 or 2 peaks
tol = 0.5e-10;
while ~((npks == 1) ...
        || ((npks == 2) && (abs(diff(loc)) >= 50)))

    % increase the smoothing parameter
    tol = tol * 2;

    % smooth the histogram
    [~, fhist2] = spaps(xhist, fhist, tol);

    % find the peaks
    [pks, loc] = findpeaks(fhist2, 'minpeakheight', 0.5e-3);

    % number of peaks found
    npks = length(pks);

end

% deal with the number of peaks
if (npks == 2)

    % we have background and tissue. The background is lighter. Sort the
    % peaks in darker to lighter order
    loc = sort(loc, 'ascend');

    % extract typical intensities for background and foreground
    mfg = xhist(loc(1));
    mbg = xhist(loc(2));

    % DEBUG
    if (DEBUG)
        subplot(2, 1, 1)
        hold off
        plot(xhist, fhist, 'b')
        hold on
        plot(xhist, fhist2, 'r')
        plot(mfg*[1 1], [0 max(fhist)], 'g')
        plot(mbg*[1 1], [0 max(fhist)], 'k')

        subplot(2, 1, 2)
        hold off
        imagesc(im(:, :, round((size(im, 3)+1)/2)))
    end

elseif (npks == 1)

    % we assume there's only background (although note that this could be a
    % case with only tissue)
    mfg = nan;
    mbg = xhist(loc);

    % DEBUG
    if (DEBUG)
        hold off
        plot(xhist, fhist, 'b')
        hold on
        plot(xhist, fhist2, 'r')
        plot(mbg*[1 1], [0 max(fhist)], 'r')
    end

else

    error('Assertion fail: The histogram has no peaks')

end

end

author: vigente

repo: gerardus


text:
function tc_plotSubplots(tc, plotMeanFlag, parent);
% tc_plotSubplots(tc, <plotMeanFlag=1>, <parent=tc.ui.plot panel>);
%
% Plots each condition in a separate subplot of
% the tc UI figure. If plotMeanFlag is set to
% 1 (default), plots mean TC with error bars;
% otherwise, plots all trials.
%
% 10/04 ras.
% 03/05 changed colors so they reflect trial #
% kgs 031605 changed colors to reflect trial number
% ras 03/16/05 amended this: will color diff't trials if 'legend'
% is selected
if nargin<1,    tc = get(gcf,'UserData');     end
if nargin<2,    plotMeanFlag = 1;             end
if nargin<3,    parent = tc.ui.plot;          end
if parent==gcf | parent==get(gcf, 'CurrentAxes')
    % make a uipanel to fit on the target
    parent = uipanel('Parent', parent, ...
                     'Units', 'normalized', ...
                     'BackgroundColor', get(gcf, 'Color'), ...
                     'Position', [0 0 1 1])
end

% params
sameAxes = 1; % if 1, will set all subplots to
              % have the same axis bounds
columnMajor = 1;   % plot conditions along columns

% init axes
axes('Parent', parent);
cla
% otherAxes = findobj('Parent', parent, 'Type', 'axes');
% delete(otherAxes);

conds = find(tc_selectedConds(tc));
nConds = length(conds);
nrows = ceil(sqrt(nConds));
ncols = ceil(nConds/nrows);

maxTrials = max(squeeze(sum(~isnan(tc.allTcs(1,:,conds)))));
% altColors = (tc.params.legend==1 & plotMeanFlag==0 & maxTrials < 10);
try
    altColors = tc.params.markEachTrial;
catch
    altColors = 0 % back compatibility
end

%% get peri-stimulus time sample points
t = tcGet(tc, 't');

%% main loop: plot each condition
for c = 1:nConds
    if columnMajor==1
        col = ceil(c/nrows);
        row = mod(c-1,nrows) + 1;
        pos = (row-1)*ncols + col;
    else
        pos = c;
    end

    subplotHandles(c) = subplot(nrows,ncols,pos);

    i = conds(c);
    cond = tc.trials.condNums(i);
    col = tc.trials.condColors{i};
    name = tc.trials.condNames{i};

    if plotMeanFlag==1
        htmp = errorbar(t, tc.meanTcs(:,i), tc.sems(:,i));
        set(htmp, 'Color', col, 'LineWidth', 2.5);
    else
        htmp = plot(t, tc.allTcs(:,:,i));
        if altColors==0
            set(htmp, 'Color', col, 'LineWidth', 1);
        else
            setLineColors( jet(maxTrials) );
        end
    end

    % indicate conditions with titles
    title(name);

    if tc.params.grid==1
        grid on
    end

    if isfield(tc.params, 'axisBounds') & ~isempty(tc.params.axisBounds)
        axis(tc.params.axisBounds);
    else
        axis tight
    end

    % if normalizing axis bounds, record them for this subplot
    if sameAxes==1
        allAX(c,:) = axis;
    end
end

% normalize axis bounds if selected
if sameAxes==1
    maxAX = [allAX(1,1) allAX(1,2) min(allAX(:,3)) max(allAX(:,4))];
    for c = subplotHandles
        subplot(c);
        axis(maxAX);
        axis off
    end
    if columnMajor==1
        corner = nrows;
    else
        corner = (nrows-1)*ncols + 1;
    end
    subplot(subplotHandles(corner)); axis on; set(gca, 'Box', 'off');
end

% label only appropriate subplots
for c = 1:nConds
    subplot(subplotHandles(c));
    if columnMajor==1
        if ceil(c/nrows)==1
			ylabel('% Signal');
        end
        if mod(c-1, nrows)+1==nrows
        	xlabel('Trial time, secs');
        end
    else
        if mod(c, ncols)==1
			ylabel('% Signal');
        end
        if ceil(c/nrows)==nrows
        	xlabel('Trial time, secs');
        end
    end
end

% indicate the peak and baseline periods, if selected
if tc.params.showPkBsl==1
    for c = 1:nConds
        subplot(subplotHandles(c));
        hold on
        AX = axis;
        plot(tc.bslPeriod, repmat(AX(3), size(tc.bslPeriod)), ...
            'k', 'LineWidth', 3);
        plot(tc.peakPeriod, repmat(AX(4), size(tc.peakPeriod)), ...
            'r', 'LineWidth', 4);
    end
end

if altColors==1
% 	% add legend that indicates trial number to last subject plot
	for j=1:tc_numTrials(tc)
        legendtext{j} = num2str(j);
    end

	hleg = legend(legendtext, -1);
    legend(hleg, 'Boxoff', 'Location', 'BestOutside')

%     axes(hleg)
%     xlabel('Trial #')
%     legendPanel(legendtext);
end

% add the handles to the user data
tc.ui.subplotHandles = subplotHandles;
set(gcf, 'UserData', tc);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function W = upgrading_merge(W, miu)
N = size(W, 2);%N is not a constant if merge is needed
if N <= miu
    return
else
%         sum_before_merge = 0;
%         for i = 1 : N/2
%           sum_before_merge = sum_before_merge + norm(W(:, i) - W(end:-1:1, N - i + 1));
%         end
%         sum_before_merge
%
%         if sum_before_merge > 0.01
%             sum(W(1, :))
%              sum(W(2, :))
%              W
%         end
%
%         if sum_before_merge ~= 0
%             sum_before_merge
%             W
%             sum(W(1, :))
%             sum(W(2, :))
%             LLR = log(W(1,:)) - log(W(2,:))
%             [~, orderd] = sort(LLR, 'descend')
%         end

    epsilon = 1e-3;
    while(N > miu)

        W_first_half = W(:, 1 : N/2);
        LR = W_first_half(1, :)./W_first_half(2, :);
        numerical_warning = 0;

        for i = 1 : N/2 - 1
            ratio = LR(i)/LR(i + 1);
            if ratio < 1 + epsilon
                numerical_warning = 1;
                break;
            end
        end

        if numerical_warning == 1

            min_deltaI = realmax;
            min_index = 0;
            for i = 1 : N/2 - 1
                a2 = W(1, i);
                b2 = W(2, i);
                a1 = W(1,i + 1);
                b1 = W(2, i + 1);
                deltaI = delta_capacity_lemma9(a1, a2, b1, b2);
                if deltaI < min_deltaI %find minimum delta I
                    min_deltaI = deltaI;
                    min_index = i;
                end
            end

            if min_index == 0
                for k = 1 : N/2
                    if sum(W(:, k)) < 1e-20
                        min_index = k;
                        W(:, min_index) = [];
                        W(:, N - min_index) = [];
                        N = size(W, 2);
                        break;
                    end
                end
                continue;
            end

            a2 = W(1, min_index);
            b2 = W(2, min_index);
            a1 = W(1, min_index + 1);
            b1 = W(2, min_index + 1);

            if a2/b2 < inf
                lambda2 = a2/b2;
                alpha2 = lambda2 * (a1 + b1)/(lambda2 + 1);
                beta2 = (a1 + b1)/(lambda2 + 1);
            else
                alpha2 = a1 + b1;
                beta2 = 0;
            end

            W(1, min_index) =  a2 + alpha2;
            W(2, min_index) =  b2 + beta2;
            W(1, N - min_index + 1) =  b2 + beta2;
            W(2, N - min_index + 1) =  a2 + alpha2;

            W(:, min_index + 1) = [];
            W(:, N - min_index - 1) = [];
            N = size(W, 2);

        else

            min_deltaI = realmax;
            min_index = 0;

            for i = 1 : N/2 - 2
                a3 = W(1, i);
                b3 = W(2, i);
                a2 = W(1,i + 1);
                b2 = W(2, i + 1);
                a1 = W(1,i + 2);
                b1 = W(2, i + 2);

                deltaI = delta_capacity_lemma11(a1, a2, a3, b1, b2, b3);

                if deltaI < min_deltaI %find minimum delta I
                    min_deltaI = deltaI;
                    min_index = i;
                end
            end

            if min_index == 0
                for k = 1 : N/2
                    if sum(W(:, k)) < 1e-20
                        min_index = k;
                        W(:, min_index) = [];
                        W(:, N - min_index) = [];
                        N = size(W, 2);
                        break;
                    end
                end
                continue;
            end

            a3 = W(1, min_index);
            b3 = W(2, min_index);
            a2 = W(1, min_index + 1);
            b2 = W(2, min_index + 1);
            a1 = W(1, min_index + 2);
            b1 = W(2, min_index + 2);

            lambda1 = a1/b1;
            if a3/b3 < inf
                lambda3 = a3/b3;
                alpha1 = lambda1 * (lambda3 * b2 - a2) / (lambda3 - lambda1);
                beta1 = (lambda3 * b2 - a2) / (lambda3 - lambda1);
                alpha3 = lambda3 * (a2 - lambda1 * b2) / (lambda3 - lambda1);
                beta3 = (a2 - lambda1 * b2) / (lambda3 - lambda1);
            else
                alpha1 = lambda1 * b2;
                beta1 = b2;
                alpha3 = a2 - lambda1 * b2;
                beta3 = 0;
            end

            W(1, min_index) =  a3 + alpha3;
            W(2, min_index) =  b3 + beta3;

            W(1, min_index + 1) =  a1 + alpha1;
            W(2, min_index + 1) =  b1 + beta1;

            W(1, N - min_index + 1) =  b3 + beta3;
            W(2, N - min_index + 1) =  a3 + alpha3;

            W(1, N - min_index) =  b1 + beta1;
            W(2, N - min_index) =  a1 + alpha1;

            W(:, min_index + 2) = [];
            W(:, N - min_index - 2) = [];
            N = size(W, 2);

%             if N == miu
%                 W
%             end
        end
    end

%         sum_after_merge = 0;
%         for i = 1 : N/2
%             sum_after_merge = sum_after_merge + norm(W(:, i) - W(end:-1:1, N - i + 1));
%         end
%         sum_after_merge
% W
end
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [x,y,xy] = q2grid(x,y,xy,mv,bound);
%q2grid   biquadratic element grid generator
%   [x,y,xy] = q2grid(x,y,xy,mv,bound);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          xy         vertex coordinate vector
%          mv         Q2 macroelement mapping matrix
%          bound      boundary vertex vector
%
% postpocesses Q2 element partitioning information to
% give standard approximation in the case of stretched grids
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
xx=xy(:,1); yy=xy(:,2); nvtx=length(xx);
%
%% recompute mid-side points in the case of stretched grids
% y-direction
yv=yy; ny=length(y);
for k=2:2:ny;
yold=y(k); ynew=0.5*(y(k+1)+y(k-1));
l=find(yy==yold); yv(l)=ynew; y(k)=ynew;
end
% x-direction
xv=xx; nx=length(x);
for k=2:2:nx;
xold=x(k); xnew=0.5*(x(k+1)+x(k-1));
l=find(xx==xold); xv(l)=xnew; x(k)=xnew;
end
xy=[xv,yv];
%
% plotting of the grid
mel=length(mv(:,1));
%if mel <= 256,
adj=sparse(nvtx,nvtx);
for i=1:mel
adj(mv(i,1),mv(i,2))=1;
adj(mv(i,2),mv(i,3))=1;
adj(mv(i,3),mv(i,4))=1;
adj(mv(i,4),mv(i,1))=1;
end
figure(10)
gplot(adj,xy,'b')
axis('square')
hold on
plot(xy(:,1),xy(:,2),'ro')
xybd=xy(bound,:);
plot(xybd(:,1),xybd(:,2),'ko')
hold off
title('Q2 finite element subdivision')
drawnow
%end
return

author: johannesgerer

repo: jburkardt-m


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

clc;
close all;
%%This example shows a trajectory planning example for the
%%antropomorphic arm

%%Link length
a2=1;
a3=1;

%%End effector initial point
%% by inverse kinematics we have the joint starting variables
[theta1s,theta2s,theta3s]=inverseAntro(1,0.5,0.3,a2,a3);

%%End effector final point
%% by inverse kinematics we have the joint ends variables
[theta1e,theta2e,theta3e]=inverseAntro(0.1,0.5,0.3,a2,a3);

%%3 joint variables: 3 trajectories
%% We want the movement be completed in 6 seconds
tstart=0;
tend=6;
time=[tstart tend];

%
% %%Check the forward kinematics and the multiple inverse solutions
figure(1);
T01=DHmatrix(theta1s(1),0,0,45);
T12=DHmatrix(theta2s(1),0,a2,0);
T23=DHmatrix(theta3s(1),0,a3,0);
Tuh1=T01*T12*T23;
figure(1);
hold on;
plotT(Tuh1);
T01=DHmatrix(theta1s(2),0,0,45);
T12=DHmatrix(theta2s(2),0,a2,0);
T23=DHmatrix(theta3s(2),0,a3,0);
Tuh2=T01*T12*T23;
plotT(Tuh2);

%%Once we have the kinematics and paths we can plot the movements!
%%Use a spline interpolation

pp1 = spline(time,[0 theta1s(1) theta1e(1) 0]);
pp2 = spline(time,[0 theta2s(1) theta2e(1) 0]);
pp3 = spline(time,[0 theta3s(1) theta3e(1) 0]);
time=linspace(tstart,tend);
figure(2);
subplot(3,1,1);
title('Position Theta1');
plot(time,fnval(pp1,time),'b');
xlabel('Time');
ylabel('Theta1');
subplot(3,1,2);
title('Position Theta2');
plot(time,fnval(pp2,time),'b');
xlabel('Time');
ylabel('Theta2');
subplot(3,1,3);
title('Position Theta3');
plot(time,fnval(pp3,time),'b');
xlabel('Time');
ylabel('Theta3');

figure(3);
for k=1:1:length(time)
clf;
theta1=fnval(pp1,time(k));
theta2=fnval(pp2,time(k));
theta3=fnval(pp3,time(k));
T01=DHmatrix(theta1,0,0,45);
T12=DHmatrix(theta2,0,a2,0);
T23=DHmatrix(theta3,0,a3,0);
Tuh1=T01*T12*T23;
plotT(T01);
plotT2(T01,T01*T12);
plotT2(T01*T12,T01*T12*T23);
pause(0.1);
title('Arm trajectory');
end

author: Sable

repo: mcbench-benchmarks


text:
function[pixAcc, meanClassPixAcc] = evaluatePixAccHierarchy_preload(imageList, probs, ~, overlapListCell, superPixelLabelHistosCell, varargin)
% [pixAcc, meanClassPixAcc] = evaluatePixAccHierarchy_preload(imageList, probs, ~, varargin)
%
% Same as evaluatePixAccPaint, but much faster by taking into account the SS hierarchy.
% To use this hierarchy, run reconstructSelSearchHierarchyFromFz().
%
% Copyright by Holger Caesar, 2015

% Parse input
p = inputParser;
addParameter(p, 'printStatus', true);
parse(p, varargin{:});

printStatus = p.Results.printStatus;

% Init
labelCount = size(probs{1}, 2);
assert(labelCount > 1);
pixCorrectHisto = zeros(labelCount, 1);
pixTotalHisto = zeros(labelCount, 1);

imageCount = numel(imageList);
for imageIdx = 1 : imageCount,
    if printStatus,
        printProgress('Evaluating pixel accuracy for image', imageIdx, imageCount);
    end;

    % Skip images without ground-truth
    if isempty(probs{imageIdx}),
        continue;
    end;

    % Precompute maximum over labels
    [maxProbs, maxInds] = max(probs{imageIdx}, [], 2);

    % Compute maximum over regions (that contain a superpixel) and count pixels
    [pixCorrectHisto, pixTotalHisto] = evaluatePixAccHierarchy_loop(maxProbs, maxInds, full(overlapListCell{imageIdx}), superPixelLabelHistosCell{imageIdx}, pixCorrectHisto, pixTotalHisto);
end;

% Compute overall accuracies
pixAcc = sum(pixCorrectHisto) / sum(pixTotalHisto);
classPixAcc = pixCorrectHisto ./ pixTotalHisto;
meanClassPixAcc = nanmean(classPixAcc);

author: nightrome

repo: matconvnet-calvin