text:
function status = besa_save2Avr(file_path, file_name, data_matrix, ...
    time_samples, channel_labels, data_scale_factor, time_scale_factor)
% BESA_SAVE2AVR writes a data matrix into ASCII-vectorized file format
%
% Parameters:
%     [file_path]
%         Full path to the folder where the file should be saved.
%
%     [file_name]
%         The name of the file where the output should be written.
%
%     [data_matrix]
%         A 2D data matrix to be saved. It should have the size
%         [NumberOfChannels x NumberOfTimeSamples].
%
%     [time_samples]
%         An array with the time samples corresponding to the data. The
%         size of the array should be [1 x NumberOfTimeSamples].
%
%     [channel_labels]
%         An array containing the channel labels for the current data.
%
%     [data_scale_factor]
%         A constant used for converting the data values from Tesla into fT
%         (10e15) or from Volt into uV (10e6). If the data is already in fT
%         or uV the just set it to 1.
%
%     [time_scale_factor]
%         A constant used for converting the time samples from seconds into
%         milliseconds (1000). If the time samples are already in
%         milliseconds just set the constant to 1.
%
% Return:
%     [status]
%         The status of the writing process. Still unused.
%

% Copyright (C) 2013, BESA GmbH
%
% File name: besa_save2Avr.m
%
% This file is part of MATLAB2BESA.
%
%    MATLAB2BESA 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.
%
%    MATLAB2BESA 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 MATLAB2BESA. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Todor Jordanov
% Created: 2012-04-17

n_tbins = size(time_samples, 2);

% Check if the matrix is [NumberOfChannels x NumberOfTimeSamples]
% or [NumberOfTimeSamples x NumberOfChannels] instead and if required
% transpose the matrix
tmp_size = size(data_matrix, 1);
% Case [NumberOfTimeSamples x NumberOfChannels]
if(n_tbins == tmp_size)

    % Transpose
    data_matrix = data_matrix';

end

n_channels = size(data_matrix, 1);

% Scale the time.
start_time = time_scale_factor * time_samples(1);
time_step = time_scale_factor * (time_samples(2) - time_samples(1));

% Scale the data.
data_matrix = data_matrix .* data_scale_factor;

% Change to the directory where the data should be saved.
cd(file_path)

% Open file for writing.
fid = fopen(file_name, 'w');

% Write the first line of the header.
fprintf(fid, ...
    'Npts= %d   TSB= %f DI= %f SB= 1.000 SC= 200.0 Nchan= %d\n', ...
    n_tbins, start_time, time_step, n_channels);

if(iscell(channel_labels))

    tmp = channel_labels{1};
    for iCh = 2:n_channels

        tmp = [tmp ' ' channel_labels{iCh}];

    end

    channel_labels = tmp;

end
% Write the second line of the header.
fprintf(fid, '%s\n', channel_labels);

% Write the data matrix to the file.
for i=1:n_channels

    for j=1:n_tbins

        fprintf(fid, '%f ', data_matrix(i, j));

    end

    fprintf(fid, '\n');

end

fclose(fid);

status = 1;


text:
classdef L2Function < handle

    properties (Constant, Access = public)
        fType = 'L2';
    end

    properties (Access = protected)
        mesh
    end

    properties (Access = private)

    end

    properties (Access = private)

    end

    methods (Access = public)

        function obj = L2Function(cParams)

        end

        function fun = project(obj,target)
            s.mesh          = obj.mesh;
            s.projectorType = target;
            proj = Projector.create(s);
            fun = proj.project(obj);
        end

    end

    methods (Access = private)

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%kkgt 'IF Input 1 is Greater Than (Input 2 or Constant), Output = TRUE'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgt.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
%
% Example: o = kkgt({i1, i2}, {'i1','';'o','';'i2',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgt - IF Input 1 is Greater Than (Input 2 or Constant), Output = TRUE
%
%  DESCRIPTION
% The GREATER THAN comparison operator compares each point
% in \fBInput 1" either with the corresponding point in \fBInput 2\fP
% or with the \fBConstant" value, which ever is specified by the user. If the
% value from \fBInput 1" is greater than the value from \fBInput 2\fP
% or the \fBConstant", the corresponding point in \fBOutput\fP will be assigned
% the value specified by \fBTRUE Value".  If the value from \fBInput 1\fP
% is not greater than the value against which it is being compared, the
% corresponding point in \fBOutput" will be the value specified by \fBFALSE
% Value".
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
%
%  "Map Data" 5
% The comparison operations have not been written to be fully polymorphic yet.
% They does not
% check for map data, and will therefore always operate on the value data,
% even if a map exists.  This will most likely corrupt indexing into the map.
% In the case of a single input, it is recommended to use the "Copy to
% Value" (kcptoval) segment operator to temporarily move the map data into
% the value segment, run the operator on the data, then move it back to the
% map with the "Copy from Value" (kcpfromval) operator.
% When operating with two input objects, where at least one of the objects
% contains map data, data should be mapped prior to running the comparison
% operator.  The "Map Data" operator (kmapdata) can be used to perform
% the mapping operation.
%
%  "Validity Mask - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Validity Mask - Two Input Objects" 5
% Masking has not been implemented yet for the comparison operators.  Therefore,
% only the mask from the first input object will be transferred to the output.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The value used to pad the data when the input files are not the same size is
% zero.
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% The comparison operations have not been extended to understand location and
% time data.  Therefore, only location and time data present in the first input
% object will be transferred to the output.
%
% Executing the Greater Than operator runs the program "kcompare" with
% the -gt flag.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kcompare
%
%  RESTRICTIONS
% This routine is still under development, and has not been modified to fully
% support the polymorphic data model.  See paragraphs above for discussions
% concerning the polymorphic data segments.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkgt(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkgt(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input'};
maxval={0,0,1};
minval={0,0,1};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kcompare"  -gt'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [index, tracker] = find_tracker(trackers, identifier)
% find_tracker Find a tracker by its identifier
%
% Find a tracker by its identifier in a cell array of tracker structures. Returns index of tracker and its structure.
% If a tracker is not found the function returns an empty matrix.
%
% Input:
% - trackers: Cell array of tracker structures.
% - identifiers: A string containing tracker identifier.
%
% Output:
% - index: Index of the tracker in the cell array or empty matrix if not found.
% - tracker: The tracker structure.

index = find(cellfun(@(t) strcmp(t.identifier, identifier), trackers, 'UniformOutput', true), 1);

if isempty(index)
    tracker = [];
else
    tracker = trackers{index};
end;

author: votchallenge

repo: toolkit-legacy


text:
function P = xyz2polar(X)
% XYZ2POLAR Transformation of cartesian coordinates to polar coordinates.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% Transformation of cartesian coordinates [x, y, z] to polar coordinates
% [distance, horizontal angle, vertical angle].
% ------------------------------------------------------------------------------
% INPUT
% X
%   Matrix of cartesian coordinates x, y and z. Each point in a row.
% ------------------------------------------------------------------------------
% OUTPUT
% P
%   Matrix of polar coordinates distace, horizontal angle and vertical angle.
%   Each point in a row.
%   Horizontal angle is defined mathematically positive starting from x-axis.
%   Vertical   angle is defined as angle between point and zenith.
% ------------------------------------------------------------------------------
% EXAMPLES
% Transformation of randomly generated points into polar coordinates and back.
%   X1   = [rand(10,1)*10 rand(10,1)*5 rand(10,1)];
%   P    = xyz2polar(X1);
%   X2   = polar2xyz(P);
%   Diff = X1-X2;
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired('X', @(x) isreal(x) && size(x,2) == 3);
p.parse(X);
p = p.Results;
% Clear required inputs to avoid confusion
clear X

% Distance ---------------------------------------------------------------------

d = sqrt(p.X(:,1).^2 + p.X(:,2).^2 + p.X(:,3).^2);

% Vertical angle ---------------------------------------------------------------

va = acosg(p.X(:,3) ./ d); % always positive [0,200]

% Horizontal angle -------------------------------------------------------------

% Compute all 4 solutions for the horizontal angle
hSol(:,1) = acosg( p.X(:,1) ./ (d.*sing(va)) ); % sol1: always positive [   0,200]
hSol(:,2) = 400 - hSol(:,1);                    % sol2: always positive [ 200,400]
hSol(:,3) = asing( p.X(:,2) ./ (d.*sing(va)) ); % sol3: can be negative [-100,100]
hSol(:,4) = 200 - hSol(:,3);                    % sol4: always positive [ 100,300]

% Eliminate complex part (occure if argument of acosg/asing is not in range
% [-1, 1]); should be the case just for few points
hSol = real(hSol);

% If some values of sol3 are negativ, transform them into positive range
idxNeg = hSol(:,3) < 0;
hSol(idxNeg,3) = hSol(idxNeg,3) + 400; % now positive

% Compute the differences between the solutions
% Pairs for comparisons: 1,3 1,4 2,3 2,4
diff = [hSol(:,1)-hSol(:,3), hSol(:,1)-hSol(:,4), hSol(:,2)-hSol(:,3), hSol(:,2)-hSol(:,4)];
% Find for each point the pair with the lowest difference
[~, idx] = min(abs(diff), [], 2);

% If pair 1 (1,3) take column 1 of hSol (i.e. sol1)
% If pair 2 (1,4) take column 1 of hSol (i.e. sol1)
% If pair 3 (2,3) take column 2 of hSol (i.e. sol2)
% If pair 4 (2,4) take column 2 of hSol (i.e. sol2)
idx(idx == 2) = 1;
idx(idx == 3) = 2;
idx(idx == 4) = 2;

% Convert subscripts to linear indices to access the elements
idxLin = sub2ind(size(hSol), [1:size(hSol,1)]', idx);

% Get values
ha = hSol(idxLin);

% Output -----------------------------------------------------------------------

P = [d ha va];

end

author: pglira


text:
function [xs,ys] = plot_gmm_1D(Priors,Mu,Sigmas)
%PLOT_GMM Summary of this function goes here
%   Detailed explanation goes here

D = size(Mu,1);
K = size(Mu,2);

[Y,I] = max(Mu);
Max = Y + 4.0*sqrt(Sigmas(I));
[Y,I] = min(Mu);
Min = Y - 4.0*sqrt(Sigmas(I));
nbSamples = 401;
xs=linspace(Min,Max, nbSamples);
ys = zeros(1,nbSamples);

for k=1:K
    ys=ys+(Priors(k).*gaussPDF(xs,Mu(k),Sigmas(k)))';
end

end

author: epfl-lasa


text:
function idx = cellregexp(cell_in,regexpstring)

    idx = logical(1 - cellfun(@isempty, regexp(cell_in, regexpstring)));
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%E2H Euclidean to homogeneous
%
% H = E2H(E) is the homogeneous version (K+1xN) of the Euclidean
% points E (KxN) where each column represents one point in R^K.
%
% Reference::
% - Robotics, Vision & Control: Second Edition, P. Corke, Springer 2016; p604.
%
% See also H2E.

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

function h = e2h(e)
    h = [e; ones(1,numcols(e))];

author: petercorke

repo: spatialmath-matlab


text:
function [] = fig2u3d(ax, fname, imgtype, addaxes, varargin)
%FIG2U3D   Convert figure to U3D file.
%   FIG2U3D saves the current axes as a U3D file for inclusion as an
%   interactive 3-dimensional figure within a PDF. Either LaTeX or Adobe
%   Acrobat can be used to embed the U3D file in the PDF.
%
%   A VWS file is also created, which contains the current camera view of
%   the axes saved. This file can be used to set the figure's default view
%   in the PDF to be the same with the open figure window in MATLAB.
%
%   The media9 LaTeX package can import U3D files with their associated VWS
%   files in a PDF document. It can be found here:
%       http://www.ctan.org/tex-archive/macros/latex/contrib/media9
%
%   For PDF readers which do not render 3D figures, it is possible to
%   include an alternative 2D image as a substitute to the 3D object.
%   For conveniency, the script saves a 2D image together with U3D file.
%   File type and other options for exporting this 2D image can be
%   specified as additional arguments.
%
%   Either export_fig (if available) or the standard print function are
%   used for 2D export. Arguments defining 2D export options need conform
%   to the export function used.
%
%   The IDTF2U3D executables are needed from:
%       http://sourceforge.net/projects/u3d/
%   To obtain them download:
%       http://www.mathworks.com/matlabcentral/fileexchange/25383-matlab-mesh-to-pdf-with-3d-interactive-object
%   and place the "bin" directory in the "idtf2u3d" directory of the
%   fig2u3d distribution.
%
% usage
%   FIG2U3D
%   FIG2U3D(ax, fname, imgtype, addaxes, varargin)
%
% optional input
%   ax = axes object handle to export (default = gca)
%   fname = file name string, e.g. 'myfigure' (default = 'surface')
%           Note: Filename extensions are appended automatically.
%   imgtype = save axes also as an image, to be used as a substitute in
%             PDF readers which cannot render the interactive 3D figure.
%           = 'filetype' | 'none'
%             where:
%                   filetype = any file format supported by the export_fig
%                              or print function (whichever is used)
%                   'none' = in this case a 2D image is not saved
%   addaxes = show axes in u3d file (default = 0)
%           = 0 | 1 (do not show/ show, respectively)
%             (depending on which one is available and used).
%   varargin = additional options to be passed to print or export_fig
%              functions, depending on which one is used for saving the
%              2D image for PDF viewers without 3D support. These must
%              conform to the inputs accepted by the export function used.
%
% output
%   This M-function does not return any data.
%   It saves a U3D file containing the axes object with handle ax.
%
% examples
%   peaks
%   FIG2U3D(gca, 'peaks', '-dpdf') % if using print
%   FIG2U3D(gca, 'peaks', '-dpng', 0, '-r300') % if using print
%
%   FIG2U3D(gca, 'peaks', '-pdf') % if using export_fig
%   FIG2U3D(gca, 'peaks', '-png') % if using export_fig
%
% remark
%   When the axes are unequal, the exporter maintains the aspect ratio by
%   changing the exported coordinates to achieve the same view as in the
%   MATLAB plot window. To avoid this, edit the "fix_daspect = 1;" line in
%   fig2u3d to "fix_daspect = 0;".
%
% dependency
%   IDTF2U3D executables, download:
%       http://www.mathworks.com/matlabcentral/fileexchange/25383-matlab-mesh-to-pdf-with-3d-interactive-object
%   and place the "bin" directory within the "idtf2u3d" directory of the
%   fig2u3d distribution.
%
% optional dependency
%   export_fig, for saving an accompanying 2D image to substitute the 3D
%   inetractive figure in PDF readers which do not render it.
%   This can be downloaded from the MATLAB Central File Exchange:
%       http://www.mathworks.com/matlabcentral/fileexchange/23629
%   Otherwise, the print function is used for the 2D image.
%
% reference
%   IDTF (Intermediate Data Text File) Format Description, Version 100,
%   Intel Corporation, 2005, available at:
%       http://u3d.svn.sourceforge.net/viewvc/u3d/releases/Gold12Update/Docs/IntermediateFormat/IDTF%20Format%20Description.pdf
%
% See also EXAMPLE_FIG2U3D, FIG2PDF3D, FIG2IDTF, IDTF2U3D, VIEW2VWS, PRINT,
%          EXPORT_FIG.
%
% File:      fig2u3d.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.10 - 2012.07.26
% Language:  MATLAB R2012a
% Purpose:   convert figure to U3D file format (lines, surfaces, points)
% Copyright: Ioannis Filippidis, 2012-

% depends
%   u3d_pre_surface, u3d_pre_line, u3d_pre_quivergroup, u3d_pre_contourgroup
%   fig2idtf, idtf2u3d, view2vws
%   axes_extremal_xyz, takehold, restorehold, quivermd, check_file_extension

% optional depends
%   export_fig

% known limitations
%   surface texturemapping not exported
%   filled contourplots not exported

% todo
%   use viewmtx

%% input
if nargin < 1
    ax = gca;
else
    % not a handle ?
    if ~ishandle(ax)
        warning('fig2u3d:handle',...
                'AX argument is not a handle. Exporting current axes.')
        ax = gca;
    end

    % not axes handle ?
    handle_type = get(ax, 'Type');
    if strcmp(handle_type, 'figure')
        warning('fig2u3d:handle',...
               ['AX argument is a figure handle.',...
                'Exporting the first axes child of this figure.'] )
        ax = get(ax, 'Children');
    end
end

% no axes object ?
if isempty(ax)
    warning('fig2u3d:emptyaxes',...
            'AX axes handle is empty. Noting to export.')
    return
end

% filename provided ?
if nargin < 2
    fname = 'surface';
elseif ~ischar(fname) || isempty(fname)
    fname = 'surface';
end

% save image substitute ?
saveimg = 1;
if nargin < 3
    imgtype = '-png';
elseif strcmp(imgtype, 'none') || isempty(imgtype)
    saveimg = 0;
elseif ~ischar(imgtype)
    imgtype = '-png';
else
    imgtype = '-pdf';
end

% show axes ?
if nargin < 4
    addaxes = 0;
elseif isempty(addaxes)
    addaxes = 0;
elseif (addaxes ~= 1) && (addaxes ~= 0)
    error('fig2u3d:addaxes', 'Argument "addaxes" can be 0 or 1.')
end

plot_axes(ax, addaxes)
delete_idtf = 1; % delete intermediate IDTF file

% nothing in the plot ?
obj = get(ax, 'Children');
if isempty(obj)
    warning('axes:empty', 'Nothing in the plot. U3D file not exported.')
    return
end

% are the axes equal ?
dar = daspect(ax);
dar = dar /dar(1);
if any(dar ~= 1)
    msg = 'Axis aspect ratio unequal. Exported U3D will look different.';
    warning('axes:aspect', msg)
end

fix_daspect = 1; % when axes are unequal, change aspect ratio of exported
                 % graphics, so that the exported figure has the same view

%% convert graphics objects to meshes, line_sets and point_sets
[surf_vertices, surf_faces, surf_facevertexcdata, surf_renderers] = u3d_pre_surface(ax);
[patch_vertices, patch_faces, patch_facevertexcdata, patch_renderers] = u3d_pre_patch(ax);
[line_vertices, line_edges, line_colors,...
                line_points, line_point_colors] = u3d_pre_line(ax);
[quiver_vertices, quiver_edges, quiver_colors] = u3d_pre_quivergroup(ax);
[contour_vertices, contour_edges, contour_colors] = u3d_pre_contourgroup(ax);

%% group meshes, line_sets and point_sets
% aggregate meshes
mesh_vertices = [surf_vertices, patch_vertices];
mesh_faces = [surf_faces, patch_faces];
mesh_colors = [surf_facevertexcdata, patch_facevertexcdata];

% aggregate lines
line_vertices = [line_vertices, quiver_vertices, contour_vertices];
line_edges = [line_edges, quiver_edges, contour_edges];
line_colors = [line_colors, quiver_colors, contour_colors];

% aggregate pointsets
pointset_points = line_points;
pointset_colors = line_point_colors;

%% fix axis aspect ratio
if fix_daspect == 1
    dar = daspect(ax);

    if ~isequal(dar, ones(1, 3) )
        msg = ['Unequal axes data aspect ratio, to be fixed by ',...
               'scaling the axes. To change this behavior, edit ',...
               '"fix_daspect = 1;" in fig2u3d to "fix_daspect = 0;"'];
        warning('fig2u3d:daspect', msg)
    end

    mesh_vertices = axis_rescale(mesh_vertices, dar);
    line_vertices = axis_rescale(line_vertices, dar);
    pointset_points = axis_rescale(pointset_points, dar);
end

%% export
fig2idtf(fname,...
          mesh_vertices, mesh_faces, mesh_colors,...
          line_vertices, line_edges, line_colors,...
          pointset_points, pointset_colors)

idtf2u3d(fname)
rm_idtf(fname, delete_idtf)

part_renderers = [surf_renderers, patch_renderers];
view2vws(ax, fname, part_renderers, fix_daspect)

save_png_substitute(ax, fname, saveimg, imgtype, varargin{:} )

function rm_idtf(fname, yesno)
if yesno == 0
    return
end

% fname extensions ok ?
fname = check_file_extension(fname, '.idtf');
fname  = fullfile(cd, fname);

system(['rm ', fname] )

function plot_axes(ax, addaxes)
if addaxes == 0
    return
end

[xyz_minmax, dim] = axes_extremal_xyz(ax);
xyz_minmax = 1.3 *xyz_minmax;

x0 = [0, 0, 0].';
vx = [xyz_minmax(2), 0, 0].';
vy = [0, xyz_minmax(4), 0].';

if dim == 3
    vz = [0, 0, xyz_minmax(6) ].';
end

held = takehold(ax);
    quivermd(ax, x0, vx, 'k')
    quivermd(ax, x0, vy, 'k')
    quivermd(ax, x0, vz, 'k')
restorehold(ax, held)

function [] = save_png_substitute(ax, fname, saveimg, imgtype, varargin)
% save image substitute ?
if saveimg == 0
    return
end

% export_fig available ?
fighandle = get(ax, 'Parent');
if exist('export_fig', 'file') == 2
    col = get(fighandle, 'Color'); % current fig background color

    % set white background for png ?
    if strcmp(imgtype, '-png')
        set(fighandle, 'Color', 'w')
    end

    export_fig(ax, imgtype, varargin{:}, fname)

    set(fighandle, 'Color', col) % restore color
else
    if strfind(imgtype, 'pdf')
        imgtype = '-dpdf';
    elseif strfind(imgtype, 'png')
        imgtype = '-dpng';
    end
    print(fighandle, imgtype, varargin{:}, fname)
end

function [a] = axis_rescale(a, dar)
if isempty(a)
   return
end

a = cellfun(@(x) cell_divider(x, dar), a, 'UniformOutput', false);

function [x] = cell_divider(x, dar)
x = [x(1, :)/dar(1);
     x(2, :)/dar(2);
     x(3, :)/dar(3) ];

author: Sable

repo: mcbench-benchmarks


text:
function [c,f,s]=rovnice(x,t,u,DuDx)
global a
% evaluates the quantities defining the
%     differential equation. The input arguments are scalars X and T and
%     vectors U and DUDX that approximate the solution and its partial
%     derivative with respect to x, respectively. PDEFUN returns column
%     vectors: C (containing the diagonal of the matrix c(x,t,u,Dx/Du)),
%     F, and S (representing the flux and source term, respectively).
c=1;
f=a*DuDx;
s=0;


text:
function producetSNEPlots(propertiesFolder,infoFilePath,reconVersion,customFeatures)
% This function plots reaction presence and uptake and secretion potential
% by taxon with t-SNE.
%
% USAGE
%   producetSNEPlots(propertiesFolder,infoFilePath,reconVersion)
%
% INPUTS
% propertiesFolder      Folder where the reaction presences and uptake and
%                       secretion potential to be analyzed are stored
%                       (default: current folder)
% infoFilePath          Path to spreadsheet with taxonomical information of
%                       the refined strains
% reconVersion          Name assigned to the reconstruction resource
% OPTIONAL INPUT
% customFeatures        Features other than taxonomy to cluster microbes
%                       by. Need to be a table header in the file with
%                       information on reconstructions.
%
%   - AUTHOR
%   Almut Heinken, 06/2020

% euclidean should work for most
distance='euclidean';
alg='barneshut';

currentDir=pwd;
cd(propertiesFolder)
mkdir('tSNE_Plots')

Summary=struct;

tol=0.0000001;

% read the file with organism information
if contains(infoFilePath,'.xlsx')
    [~,infoFile,~]=xlsread(infoFilePath);
else
    gettab=tdfread(infoFilePath);
    getcols=fieldnames(gettab);
    infoFile={};
    for j=1:length(getcols)
        infoFile{1,j}=getcols{j};
        if isnumeric(gettab.(getcols{j}))
            infoFile(2:size(gettab.(getcols{j}),1)+1,j)=cellstr(num2str(gettab.(getcols{j})));
        else
            infoFile(2:size(gettab.(getcols{j}),1)+1,j)=cellstr(gettab.(getcols{j}));
        end
    end
end

% define files to analyze
analyzedFiles={
    'Reaction presence' ['ReactionMetabolitePresence' filesep 'ReactionPresence_' reconVersion]
    'Metabolite presence' ['ReactionMetabolitePresence' filesep 'MetabolitePresence_' reconVersion]
    'Uptake and secretion potential' ['ComputedFluxes' filesep 'UptakeSecretion_' reconVersion]
    'Internal metabolite production' ['ComputedFluxes' filesep 'InternalProduction_' reconVersion]
    };

for k=1:size(analyzedFiles,1)
    if isfile([propertiesFolder filesep analyzedFiles{k,2} '.txt'])
        DataToAnalyze = readInputTableForPipeline([propertiesFolder filesep analyzedFiles{k,2} '.txt']);
        DataToAnalyze=DataToAnalyze';

        [C,I]=setdiff(DataToAnalyze(1,:),infoFile(:,1),'stable');
        DataToAnalyze(:,I(2:end))=[];

        % can only be performed if there are enough strains with taxonomical information
        if size(DataToAnalyze,2) >= 10

            if contains(version,'(R202') % for Matlab R2020a and newer
                rp=cell2mat(DataToAnalyze(2:end,2:end));
            else
                rp=str2double(DataToAnalyze(2:end,2:end));
            end
            orgs=DataToAnalyze(1,2:end)';

            taxonlevels={
                'Phylum'
                'Class'
                'Order'
                'Family'
                'Genus'
                'Species'
                };

            for i=1:length(taxonlevels)
                % plot on different taxon levels
                taxa={};
                taxcol=find(strcmp(infoFile(1,:),taxonlevels{i}));
                for j=2:size(DataToAnalyze,2)
                    if ~any(strcmp(infoFile(:,1),DataToAnalyze{1,j}))
                        taxa{j-1,1}='N/A';
                    else
                        taxa{j-1,1}=infoFile{find(strcmp(infoFile(:,1),DataToAnalyze{1,j})),taxcol};
                    end
                end

                data=rp';
                red_orgs=orgs;

                % remove entries that are all zeros
                toDel=sum(data,1)<tol;
                data(:,toDel)=[];

                % remove entries that are NaNs or empty
                findnans=any(isnan(data));
                data(:,findnans==1)=[];

                % remove unclassified organisms
                data(find(strcmp(taxa,'N/A')),:)=[];
                red_orgs(strcmp(taxa,'N/A'),:)=[];
                taxa(find(strcmp(taxa,'N/A')),:)=[];
                data(find(strcmp(taxa,'NA')),:)=[];
                red_orgs(strcmp(taxa,'NA'),:)=[];
                taxa(find(strcmp(taxa,'NA')),:)=[];

                % remove empty cells
                data(find(strcmp(taxa,'')),:)=[];
                red_orgs(strcmp(taxa,''),:)=[];
                taxa(find(strcmp(taxa,'')),:)=[];

                % remove unclassified organisms
                data(find(strncmp(taxa,'unclassified',length('unclassified'))),:)=[];
                red_orgs(find(strncmp(taxa,'unclassified',length('unclassified'))),:)=[];
                taxa(find(strncmp(taxa,'unclassified',length('unclassified'))),:)=[];

                if i==6
                    % remove unclassified species
                    toDel=[];
                    cnt=1;
                    for j=1:size(data,1)
                        if strcmp(taxa{j,1}(length(taxa{j,1})-2:length(taxa{j,1})),' sp')
                            toDel(cnt)=j;
                            cnt=cnt+1;
                        end
                    end
                    data(toDel,:)=[];
                    red_orgs(toDel,:)=[];
                    taxa(toDel,:)=[];
                end

                % remove taxa with too few members
                [uniqueXX, ~, J]=unique(taxa) ;
                occ = histc(J, 1:numel(uniqueXX));

                if length(uniqueXX) >15
                    % sort by number of entries and remove the ones with the least
                    % entries
                    [B,I]=sort(occ,'descend');
                    uniqueXX=uniqueXX(I);

                    if sum(B==1) > length(B)-15
                        % remove all that are just one entry
                        uniqueXX(B==1)=[];
                    else
                        % remove all but 20 highest
                        uniqueXX=uniqueXX(1:20);
                    end

                    [C,IA]=setdiff(taxa,uniqueXX);
                    data(find(ismember(taxa,C)),:)=[];
                    red_orgs(ismember(taxa,C),:)=[];
                    taxa(find(ismember(taxa,C)),:)=[];
                end

                % sort alphabetically
                [B,I]=sortrows(taxa,'ascend');
                data = data(I,:);
                red_orgs = red_orgs(I,:);
                taxa = taxa(I,:);

                if size(data,1)>10

                    % adjust perplicity to number of variables
                    if size(data,1) > 150
                        perpl=50;
                    elseif size(data,1) >= 50
                        perpl=30;
                    elseif size(data,1) >= 20
                        perpl=10;
                    else
                        perpl=5;
                    end

                    Y = tsne(data,'Distance',distance,'Algorithm',alg,'Perplexity',perpl,'NumDimensions',2);
                    %                 Y = tsne(data,'Distance',distance,'Algorithm',alg,'Perplexity',perpl,'NumDimensions',3);
                    Summary.([strrep(analyzedFiles{k,1},' ','_') '_' taxonlevels{i}])(:,1)=red_orgs;
                    Summary.([strrep(analyzedFiles{k,1},' ','_') '_' taxonlevels{i}])(:,2)=taxa;
                    Summary.([strrep(analyzedFiles{k,1},' ','_') '_' taxonlevels{i}])(:,3:size(Y,2)+2)=cellstr(string(Y));

                    if size(data,1) == size(Y,1) && size(Y,2) > 1
                        f=figure;
                        cols=hsv(length(unique(taxa)));
                        % define markers to better distinguish groups
                        cmarkers='';
                        for j=1:7:length(unique(taxa))
                            cmarkers=[cmarkers '+o*xsdp'];
                        end
                        cmarkers=cmarkers(1:length(unique(taxa)));
                        %                     gscatter3(Y(:,1),Y(:,2),Y(:,3),taxa,cols,{},30);
                        h=gscatter(Y(:,1),Y(:,2),taxa,cols,cmarkers);
                        set(h,'MarkerSize',4)
                        hold on
                        title(analyzedFiles{k,1})

                        h=legend('Location','northeastoutside');
                        if length(uniqueXX) < 12
                            set(h,'FontSize',12)
                        elseif length(uniqueXX) < 20
                            set(h,'FontSize',11)
                        else
                            set(h,'FontSize',8)
                        end
                        set(h, 'Interpreter', 'none')
                        grid off
                        f.Renderer='painters';
                        print(['tSNE_Plots' filesep taxonlevels{i} '_' strrep(analyzedFiles{k,1},' ','_') '_' reconVersion],'-dpng','-r300')
                    else
                        warning('Not enough strains with available organism information. Cannot cluster based on taxonomy.')
                    end
                end
            end

            % if the data should be clustered by any custom features from the info file
            if nargin > 3
                for i=1:length(customFeatures)
                    % plot on different taxon levels
                    feats={};
                    cuscol=find(strcmp(infoFile(1,:),customFeatures{i}));
                    if ~isempty(cuscol)
                        for j=2:size(DataToAnalyze,2)
                            if ~any(strcmp(infoFile(:,1),DataToAnalyze{1,j}))
                                feats{j-1,1}='N/A';
                            else
                                feats{j-1,1}=infoFile{find(strcmp(infoFile(:,1),DataToAnalyze{1,j})),cuscol};
                            end
                        end

                        data=rp';
                        red_orgs=orgs;

                        % remove organisms with no data
                        data(find(strcmp(taxa,'N/A')),:)=[];
                        red_orgs(strcmp(taxa,'N/A'),:)=[];
                        taxa(find(strcmp(taxa,'N/A')),:)=[];

                        if size(data,1) >= 10

                            % remove features with too few members
                            [uniqueXX, ~, J]=unique(feats) ;
                            occ = histc(J, 1:numel(uniqueXX));
                            toofew=uniqueXX(occ<sum(occ)/2000);
                            data(find(ismember(feats,toofew)),:)=[];
                            red_orgs(ismember(feats,toofew),:)=[];
                            feats(find(ismember(feats,toofew)),:)=[];

                            Y = tsne(data,'Distance',distance,'Algorithm',alg,'Perplexity',perpl,'NumDimensions',3);
                            Summary.([strrep(analyzedFiles{k,1},' ','_') '_' strrep(customFeatures{i},' ','_')])(:,1)=red_orgs;
                            Summary.([strrep(analyzedFiles{k,1},' ','_') '_' strrep(customFeatures{i},' ','_')])(:,2)=feats;
                            Summary.([strrep(analyzedFiles{k,1},' ','_') '_' strrep(customFeatures{i},' ','_')])(:,3:4)=cellstr(string(Y));

                            f=figure;
                            hold on
                            gscatter3(Y(:,1),Y(:,2),Y(:,3),feats);
                            title(analyzedFiles{k,1})
                            h=legend('Location','northeastoutside');
                            if length(uniqueXX) < 12
                                set(h,'FontSize',11)
                            elseif length(uniqueXX) < 20
                                set(h,'FontSize',9)
                            else
                                set(h,'FontSize',6)
                            end
                            grid off
                            set(h, 'Interpreter', 'none')
                            f.Renderer='painters';
                            print(['tSNE_Plots' filesep customFeatures{i} '_' strrep(analyzedFiles{k,1},' ','_') '_' reconVersion],'-dpng','-r300')
                        else
                            warning('Not enough strains with available organism information. Cannot cluster based on features.')
                        end
                    end
                end
            end
        end
    end
end

% export computed distance underlying the plots
datasets=fieldnames(Summary);
if length(datasets)>0
    mkdir(['tSNE_Plots' filesep 'ComputedDistances'])
    for i=1:length(datasets)
        data = vertcat({'Data point','Stratifying feature','x-axis value','y-axis value'},Summary.(datasets{i}));
        cell2csv(['tSNE_Plots' filesep 'ComputedDistances' filesep datasets{i} '.csv'],data)
    end
end

cd(currentDir)

end

author: opencobra

repo: cobratoolbox


text:
function [fusion,params] = quality_fuser_v1(w,scores,train_vecs,test_vecs,train_ndx,test_ndx,ddim)
%
% Inputs:
%
%    scores: the primary detection scores, for training
%            D-by-T matrix of T scores for D input systems
%
%    train_vecs: K1-by-M matrix, one column-vector for each of M training
%                segments
%
%    test_vecs: K2-by-N matrix, one column-vector for each of N training
%               segemnts
%
%    train_ndx: 1-by-T index where train_ndx(t) is the index into train_vecs
%               for trial t.
%
%    test_ndx: 1-by-T index where test_ndx(t) is the index into test_vecs
%              for trial t.
%    ddim: dimension of subspace for quality distandce calculation,
%          where ddim <= min(K1,K2)

if nargin==0
    test_this();
    return;
end

% Check data dimensions
[K1,M] = size(train_vecs);
[K2,N] = size(test_vecs);
assert(ddim<min(K1,K2));

[D,T] = size(scores);
assert(T == length(train_ndx));
assert(T == length(test_ndx));
assert(max(train_ndx)<=M);
assert(max(test_ndx)<=N);

% Create building blocks

[linfusion,params1] = linear_fuser([],scores);

[train_quality,params2] = sigmoid_logdistance(params1.tail,train_vecs,ddim);
train_distributor = duplicator_fh(train_ndx,size(train_vecs,2));
train_quality = train_distributor(train_quality);

[test_quality,params3] = sigmoid_logdistance(params2.tail,test_vecs,ddim);
test_distributor = duplicator_fh(test_ndx,size(test_vecs,2));
test_quality = test_distributor(test_quality);

params.get_w0 = @(ssat) [params1.get_w0(); params2.get_w0(ssat); params3.get_w0(ssat)];
params.tail = params3.tail;

% Assemble building blocks

%  combine train and test quality
quality = dottimes_of_functions([],train_quality,test_quality);

%  modulate linear fusion with quality
fusion = dottimes_of_functions([],quality,linfusion);

if ~isempty(w)
    fusion = fusion(w);
end

end

function test_this()

D = 2;
N = 5;
T = 3;
Q = 4;
ndx = ceil(T.*rand(1,N));
scores = randn(D,N);
train = randn(Q,T);
test = randn(Q,T);

ddim = 2;
ssat = 0.999;
[fusion,params] = quality_fuser_v1([],scores,train,test,ndx,ndx,ddim);

w0 = params.get_w0(ssat);
test_MV2DF(fusion,w0);

quality_fuser_v1(w0,scores,train,test,ndx,ndx,ddim),

end

author: nesl


text:
function [ x, seed ] = sphere_unit_sample_3d_2 ( seed )

%*****************************************************************************80
%
%% SPHERE_UNIT_SAMPLE_3D_2 is a BAD method for sampling the unit sphere in 3D.
%
%  Discussion:
%
%    The unit sphere in 3D satisfies:
%
%      X * X + Y * Y + Z * Z = 1
%
%    Points on the unit sphere have coordinates ( PHI, THETA ) where
%    PHI varies from 0 to PI, and THETA from 0 to 2 PI, so that:
%
%    x = cos ( theta ) * sin ( phi )
%    y = sin ( theta ) * sin ( phi )
%    z =                 cos ( phi )
%
%    This routine implements a sampling of the sphere that simply
%    picks PHI and THETA uniformly at random from their ranges.
%    This is a uniform sampling on the cylinder, but it is NOT
%    a uniform sampling on the sphere.  I implement it here just
%    so I can run some tests against the code in SPHERE_UNIT_SAMPLE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(3), the sample point.
%
%    Output, integer SEED, a seed for the random number generator.
%
  dim_num = 3;

  [ phi, seed ] = r8_uniform_01 ( seed );
  phi = pi * phi;

  [ theta, seed ] = r8_uniform_01 ( seed );
  theta = 2.0 * pi * theta;

  x(1) = cos ( theta ) * sin ( phi );
  x(2) = sin ( theta ) * sin ( phi );
  x(3) = cos ( phi );

  return
end

author: johannesgerer

repo: jburkardt-m


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

    properties
        groundObj LaunchVehicleGroundObject
        celBodyData CelestialBodyData

        timesArr(1,:) cell = {};
        xInterps(1,:) cell = {};
        yInterps(1,:) cell = {};
        zInterps(1,:) cell = {};
%         xInterpsSc(1,:) cell = {};
%         yInterpsSc(1,:) cell = {};
%         zInterpsSc(1,:) cell = {};

        vehPosVelData LaunchVehicleViewPosVelInterp

        viewFrames(1,:) AbstractReferenceFrame
        showGrdObjLoS(1,:) logical

        markerPlot(1,:) cell = {}
        losMarkerPlot(1,:) cell = {};
    end

    methods
        function obj = LaunchVehicleViewProfileGroundObjData(groundObj, vehPosVelData, celBodyData)
            obj.groundObj = groundObj;
            obj.vehPosVelData = vehPosVelData;
            obj.celBodyData = celBodyData;
        end

        function addData(obj, times, rVectsGrdObj, viewInFrame, showGrdObjLoS)
            obj.timesArr(end+1) = {times};

            obj.xInterps{end+1} = griddedInterpolant(times, rVectsGrdObj(1,:), 'spline', 'linear');
            obj.yInterps{end+1} = griddedInterpolant(times, rVectsGrdObj(2,:), 'spline', 'linear');
            obj.zInterps{end+1} = griddedInterpolant(times, rVectsGrdObj(3,:), 'spline', 'linear');

%             obj.xInterpsSc{end+1} = griddedInterpolant(times, rVectsSc(1,:), 'spline', 'linear');
%             obj.yInterpsSc{end+1} = griddedInterpolant(times, rVectsSc(2,:), 'spline', 'linear');
%             obj.zInterpsSc{end+1} = griddedInterpolant(times, rVectsSc(3,:), 'spline', 'linear');

            obj.viewFrames(end+1) = viewInFrame;
            obj.showGrdObjLoS(end+1) = showGrdObjLoS;

            obj.markerPlot{end+1} = [];
            obj.losMarkerPlot{end+1} = [];
        end

        function plotBodyMarkerAtTime(obj, time, hAx)
            [rVect, ~] = obj.vehPosVelData.getPositionVelocityAtTime(time);
            rVectInd = 1;

            for(i=1:length(obj.timesArr))
                times = obj.timesArr{i};

                if(time >= min(floor(times)) && time <= max(ceil(times)))
                    xInterp = obj.xInterps{i};
                    xGrd = xInterp(time);

                    yInterp = obj.yInterps{i};
                    yGrd = yInterp(time);

                    zInterp = obj.zInterps{i};
                    zGrd = zInterp(time);

                    if(not(isempty(obj.markerPlot{i})) && isvalid(obj.markerPlot{i}) && isa(obj.markerPlot{i}, 'matlab.graphics.primitive.Transform'))
                        M = makehgtform('translate',[xGrd,yGrd,zGrd]);
                        set(obj.markerPlot{i},'Matrix',M);
                        obj.markerPlot{i}.Visible = 'on';

                    else
                        hold(hAx,'on');
                        obj.markerPlot{i} = hgtransform('Parent', hAx);
                        hBM = plot3(hAx, 0,0,0, 'MarkerEdgeColor','k', 'Marker',obj.groundObj.markerShape.shape, 'MarkerFaceColor',obj.groundObj.markerColor.color);
                        set(hBM,'Parent',obj.markerPlot{i});

                        M = makehgtform('translate',[xGrd,yGrd,zGrd]);
                        set(obj.markerPlot{i},'Matrix',M);
                        obj.markerPlot{i}.Visible = 'on';

                        hold(hAx,'off');
                    end

                    if(obj.showGrdObjLoS(i))
%                         xInterp = obj.xInterpsSc{i};
%                         xSc = xInterp(time);
%
%                         yInterp = obj.yInterpsSc{i};
%                         ySc = yInterp(time);
%
%                         zInterp = obj.zInterpsSc{i};
%                         zSc = zInterp(time);

                        xSc = rVect(1, rVectInd);
                        ySc = rVect(2, rVectInd);
                        zSc = rVect(3, rVectInd);
                        rVectInd = rVectInd + 1;

                        if(not(isempty(obj.losMarkerPlot{i})) && isvalid(obj.losMarkerPlot{i})) % && isa(obj.markerPlot{i}, 'matlab.graphics.primitive.Transform')
                            obj.losMarkerPlot{i}.XData = [xGrd xSc];
                            obj.losMarkerPlot{i}.YData = [yGrd ySc];
                            obj.losMarkerPlot{i}.ZData = [zGrd zSc];
                        else
                            hold(hAx,'on');
                            hLoS = plot3(hAx, [xGrd xSc], [yGrd ySc], [zGrd zSc], 'g-');
                            obj.losMarkerPlot{i} = hLoS;

                            hold(hAx,'off');
                        end

                        scFrame = obj.viewFrames(i);
                        bodyInfo = scFrame.getOriginBody;
                        bodyInfoInertialFrame = bodyInfo.getBodyCenteredInertialFrame();

                        scElemSet = CartesianElementSet(time, [xSc; ySc; zSc], [0;0;0], scFrame);
                        scElemSet = scElemSet.convertToFrame(bodyInfoInertialFrame);
                        rVectScConverted = scElemSet.rVect;
                        rVectScConverted = rVectScConverted(:)';

                        stateLogEntry = [time, rVectScConverted];

                        targetBodyInfo = obj.groundObj.centralBodyInfo;

                        allBodyInfo = obj.celBodyData.getAllBodyInfo();
                        hasLoSAll = true;
                        for(j=1:length(allBodyInfo))
                            eclipseBodyInfo = allBodyInfo(j);

                            hasLoS = LoS2Target(stateLogEntry, bodyInfo, eclipseBodyInfo, targetBodyInfo, obj.celBodyData, obj.groundObj);

                            if(hasLoS == 0)
                                hasLoSAll = false;
                                break;
                            end
                        end

                        if(hasLoSAll)
                            obj.losMarkerPlot{i}.Visible = 'on';
                        else
                            obj.losMarkerPlot{i}.Visible = 'off';
                        end
                    else
                        if(not(isempty(obj.losMarkerPlot{i})) && isvalid(obj.losMarkerPlot{i})) %&& isa(obj.losMarkerPlot{i}, 'matlab.graphics.primitive.Transform')
                            obj.losMarkerPlot{i}.Visible = 'off';
                        end
                    end
                else
                    if(not(isempty(obj.markerPlot{i})) && isvalid(obj.markerPlot{i}) && isa(obj.markerPlot{i}, 'matlab.graphics.primitive.Transform'))
                        obj.markerPlot{i}.Visible = 'off';
                    end

                    if(not(isempty(obj.losMarkerPlot{i})) && isvalid(obj.losMarkerPlot{i})) %&& isa(obj.losMarkerPlot{i}, 'matlab.graphics.primitive.Transform')
                        obj.losMarkerPlot{i}.Visible = 'off';
                    end
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function map = mapcolor2(A, breaks)
%MAPCOLOR2 Create a blue-gray-red colormap.
%   MAPCOLOR2(A, BREAKS), when A is a matrix and BREAKS a vector, returns a
%   colormap that can be used as a parameter in COLORMAP function. BREAKS
%   has to contain six break values for a 7-class colormap. The break
%   values define the points at which the scheme changes color. The first
%   three colors are blue, the middle one gray and the last three ones
%   are red.
%
%   Example: A = ones(100, 100);
%            for i=1:10:100
%               A(i:end, i:end) = A(i,i) + 5;
%            end
%            A = A + rand(100, 100);
%            map = mapcolor2(A, [10 15 20 25 30 35]);
%            pcolor(A), shading flat
%            colormap(map), colorbar
%

% Copyright (c) 2006 Markus Siivola

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% there must be 6 break values for seven colors
if length(breaks) ~= 6
   error('Break point vector must have 6 elements.');
end

% sort break values in ascending order if necessary
if ~issorted(breaks)
    breaks = sort(breaks);
end

% check the value range consumed by the data
rng = [min(A(:)) max(A(:))];

if any(breaks > rng(2)) | any(breaks < rng(1))
    error('Break point out of value range');
end

if ispc
    n = 256;
else
    n = 1000;
end

% a color scheme from blue through gray to red
colors = flipud([ 140 0   0;
                  194 80  68;
                  204 143 151;
                  191 191 191;
                  128 142 207;
                  70  84  158;
                  7   39  115 ] / 255);

% create a colormap with 256 colors
map = []; n = 256; ibeg = 1;
for i = 1:6
    iend = round(n * (breaks(i) - rng(1)) / (rng(2) - rng(1)));
    map  = [map; repmat(colors(i,:), [iend-ibeg+1 1])];
    ibeg = iend+1;
end
map = [map; repmat(colors(7,:), [n - size(map, 1) 1])];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function rect = convert_rect( reccell )
    rect = zeros(0, 4, 'int16');
    for i = 1 : length(reccell)
        rect(i, :) = round([reccell{i}.x, reccell{i}.y, reccell{i}.w, reccell{i}.h]);
    end

end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function rtk=udbias_rtkins(rtk,obsr,obsb,nav,ind)

global glc;
tt=rtk.tt;  nf=rtk.NF; MAXSAT=glc.MAXSAT; NFREQ=glc.NFREQ;
sat=ind.sat; ns=ind.ns; ir=ind.ir; ib=ind.ib;

% detect cycle slip
for i=1:ns

    for f=1:rtk.opt.nf
        tmp=rtk.sat(sat(i)).slip(f);
        rtk.sat(sat(i)).slip(f)=bitand(tmp,252);%0xFC=252
    end

    % detect cycle slip using loss of lock indicator (LLI)
    rtk = detslip_LLI(rtk,obsr,ir(i),1);
    rtk = detslip_LLI(rtk,obsb,ib(i),2);

    % detect cycle slip using geometry-free method
    rtk = detslip_gf_L1L2(rtk,obsr(ir(i)),obsb(ib(i)),nav);
    rtk = detslip_gf_L1L5(rtk,obsr(ir(i)),obsb(ib(i)),nav);

    % detect cycle slip using doppler integration method
    %rtk = detslip_dop(rtk,obsr(i),1);
    %rtk = detslip_dop(rtk,obsb(i),2);

end

% update ambiguity
for f=1:nf

    % reset ambiguity if instantaneous AR mode or expire obs outage counter
    for i=1:MAXSAT
        rtk.sat(i).outc(f)=rtk.sat(i).outc(f)+1;
        if rtk.sat(i).outc(f)>rtk.opt.maxout,reset=1;else,reset=0;end

        if rtk.opt.modear==glc.ARMODE_INST && rtk.x(rtk.ib+i)~=0
            rtk=initx(rtk,0,0,rtk.ib+(f-1)*MAXSAT+i);
        elseif reset==1 && rtk.x(rtk.ib+i)~=0
            rtk=initx(rtk,0,0,rtk.ib+(f-1)*MAXSAT+i);
            rtk.sat(i).outc(f)=0;
        end

        if rtk.opt.modear~=glc.ARMODE_INST && reset==1
            rtk.sat(i).lock(f)=-rtk.opt.minlock;
        end
    end

    % reset ambiguity if detecting cycle slip
    for i=1:ns
        j=rtk.ib+(f-1)*MAXSAT+sat(i);
        rtk.P(j,j)=rtk.P(j,j)+rtk.opt.prn(1)^2*abs(tt);
        slip=rtk.sat(sat(i)).slip(f);
        if rtk.opt.ionoopt==glc.IONOOPT_IFLC
            slip=bitor(slip,rtk.sat(sat(i)).slip(f));
        end
        if rtk.opt.modear==glc.ARMODE_INST||~bitand(slip,1),continue;end
        rtk.x(j)=0;
        rtk.sat(sat(i)).lock(f)=-rtk.opt.minlock;
    end

    bias=zeros(ns,1);

    % estimate approximate ambiguity by code adn phase comparison method
    j=0; offset=0;
    for i=1:ns
        if rtk.opt.ionoopt~=glc.IONOOPT_IFLC
            cp=sdobs(obsr(ir(i)),obsb(ib(i)),f);
            pr=sdobs(obsr(ir(i)),obsb(ib(i)),NFREQ+f);
            lami=nav.lam(sat(i),f);
            if cp==0||pr==0||lami<=0,continue;end
            bias(i)=cp-pr/lami;
        else
            cp1=sdobs(obsr(ir(i)),obsb(ib(i)),1);
            cp2=sdobs(obsr(ir(i)),obsb(ib(i)),2);
            pr1=sdobs(obsr(ir(i)),obsb(ib(i)),NFREQ+1);
            pr2=sdobs(obsr(ir(i)),obsb(ib(i)),NFREQ+2);
            lam1=nav.lam(sat(i),1);
            lam2=nav.lam(sat(i),2);
            if cp1==0||cp2==0||pr1==0||pr2==0||lam1<=0||lam2<=0,continue;end

            C1= lam2^2/(lam2^2-lam1^2);
            C2=-lam1^2/(lam2^2-lam1^2);
            bias(i)=(C1*lam1*cp1+C2*lam2*cp2)-(C1*pr1+C2*pr2);
        end

        if rtk.x(rtk.ib+(f-1)*MAXSAT+sat(i))~=0
            offset=offset+(bias(i)-rtk.x(rtk.ib+(f-1)*MAXSAT+sat(i)));
            j=j+1;
        end
    end

    % correct phase-bias offset to enssure phase-code coherency
    if j>0
        for i=1:MAXSAT
            idx=rtk.ib+(f-1)*MAXSAT+i;
            if rtk.x(idx)~=0
                rtk.x(idx)=rtk.x(idx)+offset/j;
            end
        end
    end

    % set initial ambiguity
    for i=1:ns
        idx=rtk.ib+(f-1)*MAXSAT+sat(i);
        if bias(i)==0 || rtk.x(idx)~=0,continue;end
        rtk=initx(rtk,bias(i),rtk.opt.std(1)^2,idx);
    end

end

return

repo: GINav


text:
function model = mlpExpandParam(model, params)

% MLPEXPANDPARAM Update mlp model with new vector of parameters.
% FORMAT
% DESC takes a vector of MLP weights and places them in their
% respective positions in the MLP model. For single hidden layer
% neural networks the function is a wrapper for the mlpunpak command.
% ARG model : the model in which the weights are to be placed.
% ARG params : a vector of the weights to be placed in the model.
% RETURN model : the model with the weights distributed in the
% correct places.
%
% SEEALSO : mlpunpak, mlpCreate, mlpExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2007

% MLTOOLS

if length(model.hiddenDim) == 1
  model = mlpunpak(model, params);
else
  startVal = 1;
  endVal = model.inputDim*model.hiddenDim(1);
  model.w{1} = reshape(params(startVal:endVal, model.inputDim, ...
                            model.hiddenDim(1)));
  startVal = endVal + 1;
  endVal = endVal + model.hiddenDim(1);
  model.b{1} = params(startVal:endVal);
  for i = 2:length(model.hiddenDim)
    startVal = endVal + 1;
    endVal = endVal + model.hiddenDim(i-1)*model.hiddenDim(i);
    model.w{i} = reshape(params(startVal:endVal), model.hiddenDim(i-1), ...
                                model.hiddenDim(i));
    startVal = endVal + 1;
    endVal = endVal + model.hiddenDim(i);
    model.b{i} = params(startVal:endVal);
  end
  i = length(model.hiddenDim);
  startVal = endVal + 1;
  endVal = endVal + model.hiddenDim(i)*model.outputDim;
  model.w{i+1} = resphape(params(startVal:endVal), model.hiddenDim(i), ...
                          model.outputDim);
  startVal = endVal + 1;
  endVal = endVal + model.outputDim;
  model.b{i+1} = params(startVal:endVal);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function sensor_data = kspaceFirstOrder_extractSensorData(dim, sensor_data, file_index, sensor_mask_index, record, p, ux_sgx, uy_sgy, uz_sgz)
% DESCRIPTION:
%       Function to extract the required sensor data from the acoustic
%       variables at each time step. This is defined as a function rather
%       than a script to avoid the computational overhead of using scripts
%       to access variables local to another function. For k-Wave < V1.1,
%       this code was included directly in the simulation functions.
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 9th July 2013
%       last update - 6th September 2013
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% =========================================================================
% GRID STAGGERING
% =========================================================================

% shift the components of the velocity field onto the non-staggered
% grid if required for output
if (...
        record.u_non_staggered || ...
        record.I || ...
        record.I_avg)

    switch dim
        case 1
            ux_shifted = real(ifft(bsxfun(@times, record.x_shift_neg, fft(ux_sgx, [], 1)), [], 1));
        case 2
            ux_shifted = real(ifft(bsxfun(@times, record.x_shift_neg, fft(ux_sgx, [], 1)), [], 1));
            uy_shifted = real(ifft(bsxfun(@times, record.y_shift_neg, fft(uy_sgy, [], 2)), [], 2));
        case 3
            ux_shifted = real(ifft(bsxfun(@times, record.x_shift_neg, fft(ux_sgx, [], 1)), [], 1));
            uy_shifted = real(ifft(bsxfun(@times, record.y_shift_neg, fft(uy_sgy, [], 2)), [], 2));
            uz_shifted = real(ifft(bsxfun(@times, record.z_shift_neg, fft(uz_sgz, [], 3)), [], 3));
    end
end

% =========================================================================
% BINARY SENSOR MASK
% =========================================================================

if record.binary_sensor_mask

    % store the time history of the acoustic pressure
    if record.p || record.I || record.I_avg
        if ~record.compute_directivity
            sensor_data.p(:, file_index) = p(sensor_mask_index);
        end
    end

    % store the maximum acoustic pressure
    if record.p_max
        if file_index == 1
            sensor_data.p_max = p(sensor_mask_index);
        else
            sensor_data.p_max = max(sensor_data.p_max, p(sensor_mask_index));
        end
    end

    % store the minimum acoustic pressure
    if record.p_min
        if file_index == 1
            sensor_data.p_min = p(sensor_mask_index);
        else
            sensor_data.p_min = min(sensor_data.p_min, p(sensor_mask_index));
        end
    end

    % store the rms acoustic pressure
    if record.p_rms
        sensor_data.p_rms = sqrt((sensor_data.p_rms.^2*(file_index - 1) + p(sensor_mask_index).^2)./file_index);
    end

    % store the time history of the particle velocity on the staggered grid
    if record.u
        switch dim
            case 1
                sensor_data.ux(:, file_index) = ux_sgx(sensor_mask_index);
            case 2
                sensor_data.ux(:, file_index) = ux_sgx(sensor_mask_index);
                sensor_data.uy(:, file_index) = uy_sgy(sensor_mask_index);
            case 3
                sensor_data.ux(:, file_index) = ux_sgx(sensor_mask_index);
                sensor_data.uy(:, file_index) = uy_sgy(sensor_mask_index);
                sensor_data.uz(:, file_index) = uz_sgz(sensor_mask_index);
        end
    end

    % store the time history of the particle velocity
    if record.u_non_staggered || record.I || record.I_avg
        switch dim
            case 1
                sensor_data.ux_non_staggered(:, file_index) = ux_shifted(sensor_mask_index);
            case 2
                sensor_data.ux_non_staggered(:, file_index) = ux_shifted(sensor_mask_index);
                sensor_data.uy_non_staggered(:, file_index) = uy_shifted(sensor_mask_index);
            case 3
                sensor_data.ux_non_staggered(:, file_index) = ux_shifted(sensor_mask_index);
                sensor_data.uy_non_staggered(:, file_index) = uy_shifted(sensor_mask_index);
                sensor_data.uz_non_staggered(:, file_index) = uz_shifted(sensor_mask_index);
        end
    end

    % store the maximum particle velocity
    if record.u_max
        if file_index == 1
            switch dim
                case 1
                    sensor_data.ux_max = ux_sgx(sensor_mask_index);
                case 2
                    sensor_data.ux_max = ux_sgx(sensor_mask_index);
                    sensor_data.uy_max = uy_sgy(sensor_mask_index);
                case 3
                    sensor_data.ux_max = ux_sgx(sensor_mask_index);
                    sensor_data.uy_max = uy_sgy(sensor_mask_index);
                    sensor_data.uz_max = uz_sgz(sensor_mask_index);
            end
        else
            switch dim
                case 1
                    sensor_data.ux_max = max(sensor_data.ux_max, ux_sgx(sensor_mask_index));
                case 2
                    sensor_data.ux_max = max(sensor_data.ux_max, ux_sgx(sensor_mask_index));
                    sensor_data.uy_max = max(sensor_data.uy_max, uy_sgy(sensor_mask_index));
                case 3
                    sensor_data.ux_max = max(sensor_data.ux_max, ux_sgx(sensor_mask_index));
                    sensor_data.uy_max = max(sensor_data.uy_max, uy_sgy(sensor_mask_index));
                    sensor_data.uz_max = max(sensor_data.uz_max, uz_sgz(sensor_mask_index));
            end
        end
    end

    % store the minimum particle velocity
    if record.u_min
        if file_index == 1
            switch dim
                case 1
                    sensor_data.ux_min = ux_sgx(sensor_mask_index);
                case 2
                    sensor_data.ux_min = ux_sgx(sensor_mask_index);
                    sensor_data.uy_min = uy_sgy(sensor_mask_index);
                case 3
                    sensor_data.ux_min = ux_sgx(sensor_mask_index);
                    sensor_data.uy_min = uy_sgy(sensor_mask_index);
                    sensor_data.uz_min = uz_sgz(sensor_mask_index);
            end
        else
            switch dim
                case 1
                    sensor_data.ux_min = min(sensor_data.ux_min, ux_sgx(sensor_mask_index));
                case 2
                    sensor_data.ux_min = min(sensor_data.ux_min, ux_sgx(sensor_mask_index));
                    sensor_data.uy_min = min(sensor_data.uy_min, uy_sgy(sensor_mask_index));
                case 3
                    sensor_data.ux_min = min(sensor_data.ux_min, ux_sgx(sensor_mask_index));
                    sensor_data.uy_min = min(sensor_data.uy_min, uy_sgy(sensor_mask_index));
                    sensor_data.uz_min = min(sensor_data.uz_min, uz_sgz(sensor_mask_index));
            end
        end
    end

    % store the rms particle velocity
    if record.u_rms
        switch dim
            case 1
                sensor_data.ux_rms(:) = sqrt((sensor_data.ux_rms(:).^2*(file_index - 1) + ux_sgx(sensor_mask_index).^2)./file_index);
            case 2
                sensor_data.ux_rms(:) = sqrt((sensor_data.ux_rms(:).^2*(file_index - 1) + ux_sgx(sensor_mask_index).^2)./file_index);
                sensor_data.uy_rms(:) = sqrt((sensor_data.uy_rms(:).^2*(file_index - 1) + uy_sgy(sensor_mask_index).^2)./file_index);
            case 3
                sensor_data.ux_rms(:) = sqrt((sensor_data.ux_rms(:).^2*(file_index - 1) + ux_sgx(sensor_mask_index).^2)./file_index);
                sensor_data.uy_rms(:) = sqrt((sensor_data.uy_rms(:).^2*(file_index - 1) + uy_sgy(sensor_mask_index).^2)./file_index);
                sensor_data.uz_rms(:) = sqrt((sensor_data.uz_rms(:).^2*(file_index - 1) + uz_sgz(sensor_mask_index).^2)./file_index);
        end
    end

% =========================================================================
% CARTESIAN SENSOR MASK
% =========================================================================

% extract data from specified Cartesian coordinates using interpolation
% (record.record.tri and record.record.bc are the Delaunay triangulation
% and Barycentric coordinates returned by gridDataFast3D)
else

    % store the time history of the acoustic pressure
    if record.p || record.I || record.I_avg
        if dim == 1
            sensor_data.p(:, file_index) = interp1(record.grid_x, p, record.sensor_x);
        else
            sensor_data.p(:, file_index) = sum(p(record.tri) .* record.bc, 2);
        end
    end

    % store the maximum acoustic pressure
    if record.p_max
        if dim == 1
            if file_index == 1
                sensor_data.p_max = interp1(record.grid_x, p, record.sensor_x);
            else
                sensor_data.p_max = max(sensor_data.p_max, interp1(record.grid_x, p, record.sensor_x));
            end
        else
            if file_index == 1
                sensor_data.p_max = sum(p(record.tri) .* record.bc, 2);
            else
                sensor_data.p_max = max(sensor_data.p_max, sum(p(record.tri) .* record.bc, 2));
            end
        end
    end

    % store the minimum acoustic pressure
    if record.p_min
        if dim == 1
            if file_index == 1
                sensor_data.p_min = interp1(record.grid_x, p, record.sensor_x);
            else
                sensor_data.p_min = min(sensor_data.p_min, interp1(record.grid_x, p, record.sensor_x));
            end
        else
            if file_index == 1
                sensor_data.p_min = sum(p(record.tri) .* record.bc, 2);
            else
                sensor_data.p_min = min(sensor_data.p_min, sum(p(record.tri) .* record.bc, 2));
            end
        end
    end

    % store the rms acoustic pressure
    if record.p_rms
        if dim == 1
            sensor_data.p_rms = sqrt((sensor_data.p_rms.^2*(file_index - 1) + (interp1(record.grid_x, p, record.sensor_x)).^2)./file_index);
        else
            sensor_data.p_rms(:) = sqrt((sensor_data.p_rms(:).^2*(file_index - 1) + (sum(p(record.tri) .* record.bc, 2)).^2)./file_index);
        end
    end

    % store the time history of the particle velocity on the staggered grid
    if record.u
        switch dim
            case 1
                sensor_data.ux(:, file_index) = interp1(record.grid_x, ux_sgx, record.sensor_x);
            case 2
                sensor_data.ux(:, file_index) = sum(ux_sgx(record.tri) .* record.bc, 2);
                sensor_data.uy(:, file_index) = sum(uy_sgy(record.tri) .* record.bc, 2);
            case 3
                sensor_data.ux(:, file_index) = sum(ux_sgx(record.tri) .* record.bc, 2);
                sensor_data.uy(:, file_index) = sum(uy_sgy(record.tri) .* record.bc, 2);
                sensor_data.uz(:, file_index) = sum(uz_sgz(record.tri) .* record.bc, 2);
        end
    end

    % store the time history of the particle velocity
    if record.u_non_staggered || record.I || record.I_avg
        switch dim
            case 1
                sensor_data.ux_non_staggered(:, file_index) = interp1(record.grid_x, ux_shifted, record.sensor_x);
            case 2
                sensor_data.ux_non_staggered(:, file_index) = sum(ux_shifted(record.tri) .* record.bc, 2);
                sensor_data.uy_non_staggered(:, file_index) = sum(uy_shifted(record.tri) .* record.bc, 2);
            case 3
                sensor_data.ux_non_staggered(:, file_index) = sum(ux_shifted(record.tri) .* record.bc, 2);
                sensor_data.uy_non_staggered(:, file_index) = sum(uy_shifted(record.tri) .* record.bc, 2);
                sensor_data.uz_non_staggered(:, file_index) = sum(uz_shifted(record.tri) .* record.bc, 2);
        end
    end

    % store the maximum particle velocity
    if record.u_max
        if file_index == 1
            switch dim
                case 1
                    sensor_data.ux_max = interp1(record.grid_x, ux_sgx, record.sensor_x);
                case 2
                    sensor_data.ux_max = sum(ux_sgx(record.tri) .* record.bc, 2);
                    sensor_data.uy_max = sum(uy_sgy(record.tri) .* record.bc, 2);
                case 3
                    sensor_data.ux_max = sum(ux_sgx(record.tri) .* record.bc, 2);
                    sensor_data.uy_max = sum(uy_sgy(record.tri) .* record.bc, 2);
                    sensor_data.uz_max = sum(uz_sgz(record.tri) .* record.bc, 2);
            end
        else
            switch dim
                case 1
                    sensor_data.ux_max = max(sensor_data.ux_max, interp1(record.grid_x, ux_sgx, record.sensor_x));
                case 2
                    sensor_data.ux_max = max(sensor_data.ux_max, sum(ux_sgx(record.tri) .* record.bc, 2));
                    sensor_data.uy_max = max(sensor_data.uy_max, sum(uy_sgy(record.tri) .* record.bc, 2));
                case 3
                    sensor_data.ux_max = max(sensor_data.ux_max, sum(ux_sgx(record.tri) .* record.bc, 2));
                    sensor_data.uy_max = max(sensor_data.uy_max, sum(uy_sgy(record.tri) .* record.bc, 2));
                    sensor_data.uz_max = max(sensor_data.uz_max, sum(uz_sgz(record.tri) .* record.bc, 2));
            end
        end
    end

    % store the minimum particle velocity
    if record.u_min
        if file_index == 1
            switch dim
                case 1
                    sensor_data.ux_min = interp1(record.grid_x, ux_sgx, record.sensor_x);
                case 2
                    sensor_data.ux_min = sum(ux_sgx(record.tri) .* record.bc, 2);
                    sensor_data.uy_min = sum(uy_sgy(record.tri) .* record.bc, 2);
                case 3
                    sensor_data.ux_min = sum(ux_sgx(record.tri) .* record.bc, 2);
                    sensor_data.uy_min = sum(uy_sgy(record.tri) .* record.bc, 2);
                    sensor_data.uz_min = sum(uz_sgz(record.tri) .* record.bc, 2);
            end
        else
            switch dim
                case 1
                    sensor_data.ux_min = min(sensor_data.ux_min, interp1(record.grid_x, ux_sgx, record.sensor_x));
                case 2
                    sensor_data.ux_min = min(sensor_data.ux_min, sum(ux_sgx(record.tri) .* record.bc, 2));
                    sensor_data.uy_min = min(sensor_data.uy_min, sum(uy_sgy(record.tri) .* record.bc, 2));
                case 3
                    sensor_data.ux_min = min(sensor_data.ux_min, sum(ux_sgx(record.tri) .* record.bc, 2));
                    sensor_data.uy_min = min(sensor_data.uy_min, sum(uy_sgy(record.tri) .* record.bc, 2));
                    sensor_data.uz_min = min(sensor_data.uz_min, sum(uz_sgz(record.tri) .* record.bc, 2));
            end
        end
    end

    % store the rms particle velocity
    if record.u_rms
        switch dim
            case 1
                sensor_data.ux_rms = sqrt((sensor_data.ux_rms.^2*(file_index - 1) + (interp1(record.grid_x, ux_sgx, record.sensor_x)).^2)./file_index);
            case 2
                sensor_data.ux_rms(:) = sqrt((sensor_data.ux_rms(:).^2*(file_index - 1) + (sum(ux_sgx(record.tri) .* record.bc, 2)).^2)./file_index);
                sensor_data.uy_rms(:) = sqrt((sensor_data.uy_rms(:).^2*(file_index - 1) + (sum(uy_sgy(record.tri) .* record.bc, 2)).^2)./file_index);
            case 3
                sensor_data.ux_rms(:) = sqrt((sensor_data.ux_rms(:).^2*(file_index - 1) + (sum(ux_sgx(record.tri) .* record.bc, 2)).^2)./file_index);
                sensor_data.uy_rms(:) = sqrt((sensor_data.uy_rms(:).^2*(file_index - 1) + (sum(uy_sgy(record.tri) .* record.bc, 2)).^2)./file_index);
                sensor_data.uz_rms(:) = sqrt((sensor_data.uz_rms(:).^2*(file_index - 1) + (sum(uz_sgz(record.tri) .* record.bc, 2)).^2)./file_index);
        end
    end
end

% =========================================================================
% RECORDED VARIABLES OVER ENTIRE GRID
% =========================================================================

% store the maximum acoustic pressure over all the grid elements
if record.p_max_all
    switch dim
        case 1
            if file_index == 1
                sensor_data.p_max_all = p(record.x1_inside:record.x2_inside);
            else
                sensor_data.p_max_all = max(sensor_data.p_max_all, p(record.x1_inside:record.x2_inside));
            end
        case 2
            if file_index == 1
                sensor_data.p_max_all = p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside);
            else
                sensor_data.p_max_all = max(sensor_data.p_max_all, p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside));
            end
        case 3
            if file_index == 1
                sensor_data.p_max_all = p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
            else
                sensor_data.p_max_all = max(sensor_data.p_max_all, p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
            end
    end
end

% store the minimum acoustic pressure over all the grid elements
if record.p_min_all
    switch dim
        case 1
            if file_index == 1
                sensor_data.p_min_all = p(record.x1_inside:record.x2_inside);
            else
                sensor_data.p_min_all = min(sensor_data.p_min_all, p(record.x1_inside:record.x2_inside));
            end
        case 2
            if file_index == 1
                sensor_data.p_min_all = p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside);
            else
                sensor_data.p_min_all = min(sensor_data.p_min_all, p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside));
            end
        case 3
            if file_index == 1
                sensor_data.p_min_all = p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
            else
                sensor_data.p_min_all = min(sensor_data.p_min_all, p(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
            end
    end
end

% store the maximum particle velocity over all the grid elements
if record.u_max_all
    switch dim
        case 1
            if file_index == 1
                sensor_data.ux_max_all = ux_sgx(record.x1_inside:record.x2_inside);
            else
                sensor_data.ux_max_all = max(sensor_data.ux_max_all, ux_sgx(record.x1_inside:record.x2_inside));
            end
        case 2
            if file_index == 1
                sensor_data.ux_max_all = ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside);
                sensor_data.uy_max_all = uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside);
            else
                sensor_data.ux_max_all = max(sensor_data.ux_max_all, ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside));
                sensor_data.uy_max_all = max(sensor_data.uy_max_all, uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside));
            end
        case 3
            if file_index == 1
                sensor_data.ux_max_all = ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
                sensor_data.uy_max_all = uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
                sensor_data.uz_max_all = uz_sgz(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
            else
                sensor_data.ux_max_all = max(sensor_data.ux_max_all, ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
                sensor_data.uy_max_all = max(sensor_data.uy_max_all, uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
                sensor_data.uz_max_all = max(sensor_data.uz_max_all, uz_sgz(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
            end
    end
end

% store the minimum particle velocity over all the grid elements
if record.u_min_all
    switch dim
        case 1
            if file_index == 1
                sensor_data.ux_min_all = ux_sgx(record.x1_inside:record.x2_inside);
            else
                sensor_data.ux_min_all = min(sensor_data.ux_min_all, ux_sgx(record.x1_inside:record.x2_inside));
            end
        case 2
            if file_index == 1
                sensor_data.ux_min_all = ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside);
                sensor_data.uy_min_all = uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside);
            else
                sensor_data.ux_min_all = min(sensor_data.ux_min_all, ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside));
                sensor_data.uy_min_all = min(sensor_data.uy_min_all, uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside));
            end
        case 3
            if file_index == 1
                sensor_data.ux_min_all = ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
                sensor_data.uy_min_all = uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
                sensor_data.uz_min_all = uz_sgz(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside);
            else
                sensor_data.ux_min_all = min(sensor_data.ux_min_all, ux_sgx(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
                sensor_data.uy_min_all = min(sensor_data.uy_min_all, uy_sgy(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
                sensor_data.uz_min_all = min(sensor_data.uz_min_all, uz_sgz(record.x1_inside:record.x2_inside, record.y1_inside:record.y2_inside, record.z1_inside:record.z2_inside));
            end
    end
end

repo: Aero-matlab


text:
% Copyright Mohammad SAFEEA, 17th-Aug-2017
function [ y ] = getTheKSTDirectory( x )
%% This function is used to get the directory of the KST matlab client library
n=max(size(x));
KSTfolderName='Matlab_client';
for i=1:n
    if(x(end)=='\')
        break;
    else
        x=x(1:end-1);
    end
end
y=[x,KSTfolderName];
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [logww] = free_energy(vv,beta_k1,beta_k0,visbiases_base,visbiases,hidbias,vishid)

 [numcases,numdims] = size(vv);
 temp_base = vv*visbiases_base';
 temp_bias = vv*visbiases';
 tempvWh = hidbias + (vv*vishid);

 aa = beta_k1;
 p_star_k1  = (1-aa)*temp_base  + aa*temp_bias + sum(log(1+exp(aa*tempvWh)),2);

 aa = beta_k0;
 p_star_k0  =  (1-aa)*temp_base + aa*temp_bias + sum(log(1+exp(aa*tempvWh)),2);

 logww = (p_star_k1 - p_star_k0);

repo: DeepFi