text:
%   CLASS Settings_Interface
% =========================================================================
%
% DESCRIPTION
%   Abstract class with basic properties and methods that a setting class
%   must have
%
% COMMENTS
% Settings have been build with multiple inheritance
% A standard abstract interface have been created: Setting_Interface
% it add to each subclass the object log
% force the subclasses to implement three basic methods:
%  - import:       when a setting object of identical class, or that inherits from the same class is passed to this function, the relative parameters are copied in the calling object
%                  when an ini file is the input of the function, the object is updated with the settings contained into the file
%  - toString:     display the content of the object, in a human readable way, a goGPS user can "ask" for the value of the settings on screen
%  - export:       create a cell array of strings containing the settings in plain text ini format. The variable it's the raw data format of Ini_Manager
%

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef Settings_Interface < Exportable
    properties (Abstract)
    end

    methods  (Abstract)
        import(obj, settings);
        % This function import the settings of the current class from another setting object having the same properties, or an ini file

        str = toString(obj, str)
        % Display the content of the class in a human readable format

        str_cell = export(obj, str_cell)
        % Conversion to string ini format of the minimal information needed to reconstruct the obj
    end

    methods (Access = 'public')
        function ini = save(this, file_path)
            % Save to a file (in INI fomat) the content of the Settings object
            % SYNTAX: <ini> = this.save(file_path);
            % return optionally the ini manager object used by the save function
            ini = Ini_Manager(file_path, this.export());
        end

        function [home_reset] = importIniFile(this, file_path)
            % Import from an INI file the content of the Settings object
            %
            % OUTPUT
            %   home_reset is empty if the home have not been changed
            %              otherwise it contains the old path
            % SYNTAX
            %   home_reset = this.importIniFile(file_path);
            ini = Ini_Manager(file_path);
            fnp = File_Name_Processor;
            prj_home = fnp.checkPath(ini.getData('prj_home'));
            home_reset = prj_home;
            if ~isempty(prj_home) && ~exist(prj_home, 'dir')
                [path_str, ~, ~] = fileparts(fnp.getFullDirPath(file_path));
                prj_home = fnp.getFullDirPath(strcat(path_str, [filesep '..' filesep]), pwd);
                if exist(prj_home, 'dir')
                    ini.setData('prj_home', prj_home);
                else
                    home_reset = '';
                end
            else
                home_reset = '';
            end
            this.import(ini);
        end

        function updateSettingsFile(this, file_path)
            this.importIniFile(file_path);
            this.save(file_path);
        end

    end

    % =========================================================================
    %%  SETTERS IO
    % =========================================================================
    methods
        function value = setProperty(this, prop_name, value)
            this.(prop_name) = value;
        end
    end

    % =========================================================================
    %%  GETTERS IO
    % =========================================================================
    methods
        function value = getProperty(this, prop_name)
            props = fieldnames(this);
            if any(strcmp(props,prop_name))
                value = this.(prop_name);
            else
                value = [];
            end
        end
    end

    % =========================================================================
    %  TEST PARAMETERS VALIDITY
    % =========================================================================
    methods (Access = 'protected')

        function checkLogicalField(this, field_name)
            % Check if a logical field of the object is a valid logical number
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            % This superclass function must be copied in each child (to
            % have read / write permission on the parameters
            %
            % SYNTAX
            %   this.checkLogicalField(string_field_name);
            this.(field_name) = this.checkLogical(field_name, this.(field_name), this.(upper(field_name)));
        end

        function is_existing = checkStringField(this, field_name, empty_is_valid, check_existence)
            % Check if a string field of the object is a valid string
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            % This superclass function must be copied in each child (to
            % have read / write permission on the parameters
            %
            % SYNTAX
            %   this.checkStringField(string_field_name, <empty_is_valid == false>, <check_existence == false>);
            switch nargin
                case 2, [this.(field_name), is_existing] = this.checkString(field_name, this.(field_name), this.(upper(field_name)));
                case 3, [this.(field_name), is_existing] = this.checkString(field_name, this.(field_name), this.(upper(field_name)), empty_is_valid);
                case 4, [this.(field_name), is_existing] = this.checkString(field_name, this.(field_name), this.(upper(field_name)), empty_is_valid, check_existence);
                otherwise, error('Settings checkStringField called with the wrong number of parameters');
            end
        end

        function checkNumericField(this, field_name, limits, valid_val)
            % Check if a numeric field of the object is valid
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            % This superclass function must be copied in each child (to
            % have read / write permission on the parameters
            %
            % SYNTAX
            %   this.checkNumericField(string_field_name, <limits>, <valid_values>);
            switch nargin
                case 2, this.(field_name) = this.checkNumber(field_name, this.(field_name), this.(upper(field_name)));
                case 3, this.(field_name) = this.checkNumber(field_name, this.(field_name), this.(upper(field_name)), limits);
                case 4, this.(field_name) = this.checkNumber(field_name, this.(field_name), this.(upper(field_name)), limits, valid_val);
                otherwise, error('Settings checkNumericField called with the wrong number of parameters');
            end
        end

        function checked_val = checkLogical(this, field_name, field_val, default_val)
            % Check if a logical is a valid logical
            % This superclass function must be called in each child
            %
            % SYNTAX
            %   checked_val = this.checkLogicalField(string_field_name);
            checked_val = default_val;
            if (~isempty(field_val)) && (~isnan(field_val))
                checked_val = logical(field_val);
            else
                Core.getLogger.addWarning(sprintf('The settings field %s is not valid => using default %d', field_name, checked_val));
            end
        end

        function [checked_val, is_existing] = checkCellString(this, field_name, field_val, default_val, empty_is_valid, check_existence)
            % Check if a string is a valid string
            % check_existence 0: do not check
            %                 1: check and do nothing
            %                 2: check and try to correct
            % SYNTAX
            %   checked_val = this.checkString(string_field_name, <empty_is_valid == false>, <check_existence == false>);
            if (nargin < 5)
                empty_is_valid = false;
            end
            if (nargin < 6)
                check_existence = false;
            end

            [checked_val, is_existing] = this.checkString(field_name, field_val, default_val, empty_is_valid, check_existence);
            if ~iscell(checked_val)
                checked_val = {checked_val};
            end
        end

        function [checked_val, is_existing] = checkString(this, field_name, field_val, default_val, empty_is_valid, check_existence)
            % Check if a string is a valid string
            % check_existence 0: do not check
            %                 1: check and do nothing
            %                 2: check and try to correct
            %
            % SYNTAX
            %   [checked_val, is_existing] = this.checkString(string_field_name, <empty_is_valid == false>, <check_existence == 0>);
            if (nargin < 5)
                empty_is_valid = false;
            end
            if (nargin < 6)
                check_existence = 0;
            end

            is_existing = false;
            checked_val = default_val;
            default_val = field_val;
            if ~isempty(field_val) && iscell(field_val) % A cell of strings must contain at least one string
                field_val = field_val{1};
            end
            %if (ischar(field_val) || (empty_is_valid && isempty(field_val))) && ((~isempty(field_val)) || empty_is_valid) && ((exist(field_val,'file') || exist(field_val,'dir')) || ~check_existence)
            if ischar(field_val)  && (empty_is_valid || ~isempty(field_val)) && (~check_existence || exist(field_val, 'file'))
                is_existing = true;
                checked_val = default_val;
            else
                if check_existence
                    if exist(default_val, 'file')
                        checked_val = default_val;
                    else
                        checked_val = field_val;
                    end
                end

                if isempty(field_val) || ~any(strcmp(checked_val, field_val))
                    log = Core.getLogger();
                    if check_existence > 1 || ~empty_is_valid
                        if iscell(checked_val)
                            log.addWarning(sprintf('The value "%s" of the settings field %s is not valid => using default: "%s"', iif(isempty(field_val), 'empty', field_val), field_name, Ini_Manager.strCell2Str(checked_val)));
                        else
                            log.addWarning(sprintf('The value "%s" of the settings field %s is not valid => using default: "%s"', iif(isempty(field_val), 'empty', field_val), field_name, checked_val));
                        end
                    else
                        if empty_is_valid
                            checked_val = field_val;
                        else
                            log.addWarning(sprintf('The value "%s" of the settings field %s is not valid!!!', iif(isempty(field_val), 'empty', field_val), field_name));
                        end
                    end
                end
            end
        end

        function checked_val = checkNumber(this, field_name, field_val, default_val, limits, valid_val)
            % Check if a number is valid
            %
            % SYNTAX
            %   checked_val = this.checkNumericField(string_variable_name, value, default_value, <limits>, <valid_values>);
            checked_val = default_val;
            log = Core.getLogger();
            if isnumeric(field_val) && (~isempty(field_val)) && (any(~isnan(field_val)))
                checked_val = field_val;
                % if I have limits to check => check for out of bound
                for i = 1 : numel(checked_val)
                    if (nargin >= 5) && (numel(limits) == 2) && ...
                            ((checked_val(i) > limits(2)) || (checked_val(i) < limits(1)))
                        checked_val(i) = max(limits(1), min(limits(2), checked_val(i)));
                        log.addWarning(sprintf('The value %g of the settings field %s is not within the valid limits (%g .. %g) => updating it to %g', field_val(i), field_name, limits(1),limits(2), checked_val(i)));
                    end
                    % if I have a set of values => check for set intersection
                    if (nargin >= 6) && (~isempty(valid_val)) && (~ismember(checked_val(i), valid_val))
                        checked_val(i) = default_val;
                        log.addWarning(sprintf('The value %g for the settings field %s is not valid => using default %g. It should be one of: %s', field_val(i), field_name, checked_val(i), sprintf('%g ', valid_val)));
                    end
                end
            else
                log.addWarning(sprintf('The settings field %s is not valid => using default [%s ]', field_name, sprintf(' %g', checked_val)));
            end
        end
    end

    % =========================================================================
    %  TEST
    % =========================================================================

    methods (Access = 'protected')
        function testInterfaceRoutines(this)
            % test the class (Interface Routines)
            % SINTAX: this.testInterfaceRoutines();

            log = Core.getLogger;
            try
                vl = log.getVerbosityLev();
                log.setVerbosityLev(1e3);
                test = this;
                raw_data = test.export();

                ini = Ini_Manager('test__.ini', raw_data);
                ini.showData();

                test.import(ini);
                test_copy = repmat(test,2,1); test_copy = test_copy(2);
                test.import(test_copy);
                clear test_copy;
                fprintf('\n');
                disp(test.toString());
                %delete('test__.ini');
            catch ex
                log.addError(['Test failed: ' ex.message]);
                Core_Utils.printEx(ex);
            end
            log.setVerbosityLev(vl);
        end
    end
end

author: goGPS-Project


text:
function [obj] = update(obj)
    % Updates the NLP problems before load it to NLP solver

    % update variable indices
    index_offset = 0;

    num_vars = numel(obj.VariableArray);

    for i = 1:num_vars
        dim = obj.VariableArray(i).Dimension;

        % set the index
        setIndices(obj.VariableArray(i),index_offset + cumsum(ones(dim,1)));

        % increments (updates) offset
        index_offset = index_offset + dim;
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function z =fonction(input)

x = input(1);
y = input(2);

z =  exp(-(x-0.5)^2)+exp(-(y-0.7)^2)+(1/25)*exp(-(x+0.2)^2)+(1/25)*exp(-(x-0.2)^2)+(1/25)*exp(-(y+1)^2)+(1/25)*exp(-(y-1)^2);

drawnow;

author: Sable

repo: mcbench-benchmarks


text:
function [kWh] = Nm2kWh(Nm)
% Convert energy or work from newtons-meters to kilowatt-hours.
% Chad A. Greene 2012
kWh = Nm*2.7777777778e-7;


text:
function value = scnrm2 ( n, x, incx )

%*****************************************************************************80
%
%% SCNRM2 returns the euclidean norm of a complex vector.
%
%  Discussion:
%
%    SCNRM2 := sqrt ( sum ( conjg ( x(1:n) ) * x(1:n) ) )
%            = sqrt ( dot_product ( x(1:n), x(1:n) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for FORTRAN usage,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, pages 308-323, 1979.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, complex X(*), the vector.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Output, real VALUE, the norm of the vector.
%
  value = sqrt ( sum ( ( real ( x(1:incx:1+(n-1)*incx) ) ).^2 ) ...
               + sum ( ( imag ( x(1:incx:1+(n-1)*incx) ) ).^2 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


%% planRRT
%%  - basic coverage algorithm
%%
%% Last Modified - 11/15/2010 - R. Beard

function path=planCoverRRTDubins(wpp_start, R, map)

    % desired down position is down position of start node
    pd = wpp_start(3);

    % specify start node from wpp_start
    start_node = [wpp_start(1), wpp_start(2), pd, 0, 0, 0];
    % format is [N, E, D, chi, cost, parent]

    % return map
    returnMapSize = 30;  % this is a critical parameter!
    return_map = 50*ones(returnMapSize,returnMapSize)+ rand(returnMapSize,returnMapSize);
    plotReturnMap(return_map), %pause

    % construct search path by doing N search cycles
    SEARCH_CYCLES = 100;  % number of search cycles

    % look ahead tree parameters
    L = 2*R+1;  % segment Length
    vartheta = pi/4; % search angle
    depth = 5; % depth of look ahead tree

    % initialize path and tree
    path = start_node;
    for i=1:SEARCH_CYCLES,
        tree = extendTree(path(end,:),L,vartheta,depth,map,return_map,pd);
        next_path = findMaxReturnPath(tree);
        path = [path; next_path(1,:)];
        % update the return map
        return_map = updateReturnMap(next_path(1,:),return_map,map);
        plotReturnMap(return_map), %pause
        % set the end of the path as the root of the tree
    end

    % remove path segments where there is no turn
    path_=path;
    path = path(1,:);
    for i=2:size(path_,1),
        if path_(i,4)~=path_(i-1,4);
            path = [path; path_(i,:)];
        end
    end

%     % specify that these are straight-line paths.
%     for i=1:size(path,1),
%         path(i,4)=-9999;
%     end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% extendTree
%%   extend tree by randomly selecting point and growing tree toward that
%%   point
function tree = extendTree(start_node,L,vartheta,depth,map,return_map,pd)

  tree_ = [start_node, 0];  % the last variable marks node as expanded

  % extend tree initially
  for d = 1:depth,
      newnodes = [];
      for j=1:size(tree_,1),
          if tree_(j,7)~=1, % process unexpanded nodes
              for i=1:3,
                  if     i==1, theta = tree_(j,4)-vartheta;
                  elseif i==2, theta = tree_(j,4);
                  elseif i==3, theta = tree_(j,4)+vartheta;
                  end
                  newnode_ = [...
                      tree_(j,1)+L*cos(theta),...
                      tree_(j,2)+L*sin(theta),...
                      tree_(j,3),...
                      theta,...
                      0,...
                      j,...
                      0,...
                      ];
                  if collision(tree_(j,:), newnode_, map)==0,
                      newnode_(5) = tree_(j,5)...
                          +findReturn(newnode_(1),newnode_(2),return_map,map);
                      newnodes = [newnodes; newnode_];
                  end
              end
              tree_(j,7)=1; % mark as expanded
          end
      end
      tree_ = [tree_; newnodes];
  end
  tree = tree_(:,1:6);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% collision
%%   check to see if a node is in collsion with obstacles
function collision_flag = collision(start_node, end_node, map)
    collision_flag = 0;
    sigma = 0:.1:1; % define where collisions will be checked
    for i=1:length(sigma),
        X = (1-sigma(i))*start_node(1) + sigma(i)*end_node(1);
        Y = (1-sigma(i))*start_node(2) + sigma(i)*end_node(2);
        Z = (1-sigma(i))*start_node(3) + sigma(i)*end_node(3);
        if Z >= downAtNE(map, X, Y),
            collision_flag = 1;
        end
        % check to see if outside of world
        if (X>map.width) | (X<0) | (Y>map.width) | (Y<0),
            collision_flag = 1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% downAtNE
%%   find the map down coordinate at a specified (n,e) location
function down = downAtNE(map, n, e)

      [d_n,idx_n] = min(abs(n - map.buildings_n));
      [d_e,idx_e] = min(abs(e - map.buildings_e));

      if (d_n<=map.BuildingWidth) & (d_e<=map.BuildingWidth),
          down = -map.heights(idx_e,idx_n);
      else
          down = 0;
      end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% findReturn
%%   compute the return value at a particular location
function return_value = findReturn(pn,pe,return_map,map);

  [pn_max,pe_max] = size(return_map);
  fn = pn_max*pn/map.width;
  fn = min(pn_max,round(fn));
  fn = max(1,fn);
  fe = pe_max*pe/map.width;
  fe = min(pe_max,round(fe));
  fe = max(1,fe);
  return_value = return_map(fn,fe);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% findMaxReturnPath
%%   find the maximum return path in the tree
function path = findMaxReturnPath(tree)

     % find node with max return
    [tmp,idx] = max(tree(:,5));

    % construct path with maximum return
    path = tree(idx,:);
    parent_node = tree(idx,6);
    %while parent_node>1,
    while parent_node>1,
        path = [tree(parent_node,:); path];
        parent_node = tree(parent_node,6);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% updateReturnMap
%%   update the return map to indicate where MAV has been
function new_return_map = updateReturnMap(path,return_map,map)

  new_return_map = return_map;
  for i=1:size(path,1),
    pn = path(i,1);
    pe = path(i,2);
    [pn_max,pe_max] = size(return_map);
    fn = pn_max*pn/map.width;
    fn = min(pn_max,round(fn));
    fn = max(1,fn);
    fe = pe_max*pe/map.width;
    fe = min(pe_max,round(fe));
    fe = max(1,fe);

    new_return_map(fn,fe) = return_map(fn,fe) - 50;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plotReturnMap
%%   plot the return map
function plotReturnMap(return_map)

  figure(2), clf
  mesh(return_map)
end


text:
function ua = diagonal_pointer_cr ( n, nz_num, ia, ja )

%*****************************************************************************80
%
%% DIAGONAL_POINTER_CR finds diagonal entries in a sparse compressed row matrix.
%
%  Discussion:
%
%    The matrix A is assumed to be stored in compressed row format.  Only
%    the nonzero entries of A are stored.  The vector JA stores the
%    column index of the nonzero value.  The nonzero values are sorted
%    by row, and the compressed row vector IA then has the property that
%    the entries in A and JA that correspond to row I occur in indices
%    IA[I] through IA[I+1]-1.
%
%    The array UA can be used to locate the diagonal elements of the matrix.
%
%    It is assumed that every row of the matrix includes a diagonal element,
%    and that the elements of each row have been ascending sorted.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2008
%
%  Author:
%
%    Original C version by Lili Ju
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the system.
%
%    Input, integer NZ_NUM, the number of nonzeros.
%
%    Input, integer IA(N+1), JA(NZ_NUM), the row and column indices
%    of the matrix values.  The row vector has been compressed.  On output,
%    the order of the entries of JA may have changed because of the sorting.
%
%    Output, integer UA(N), the index of the diagonal element of each row.
%
  ua(1:n) = -1;

  for i = 1 : n
    for j = ia(i) : ia(i+1) - 1
      if ( ja(j) == i )
        ua(i) = j;
      end
    end
  end

  return
end


text:
function [km] = A2km(A)
% Convert length from angstroms to kilometers.
% Chad A. Greene 2012
km = A*1e-13;


text:
function logPr = calcLogJointPr_RGS_BPHMM( Psi, data_struct, hyperparams, model, objIDs, Xstats )
% Calculate logarithm of the joint probability
%   of the state Psi output by RESTRICTED scan
% OUTPUT
%   logPr := struct with fields
%              .F
%              .z
%              .obs

logPr.eta = 0;
logPr.z = 0;
logPr.theta = 0;

featIDs = unique(Psi.activeFeatIDs);

alpha0 = hyperparams.alpha;
kappa0 = hyperparams.kappa;

logPr.F = calcLogPrFeatureMatrix( Psi.F, hyperparams.gamma, hyperparams.c );

Zstats = getZSuffStats( Psi.F, Psi.stateSeq, model, objIDs );
switch model.HMMmodel.transType % ---------------------------------------------- log p( z | F, alph, kappa)
    case 'byObject'
        for ii = objIDs
            Pz = Psi.TS.obj(ii).pi_z;
            Pi = bsxfun( @rdivide, Pz, sum(Pz,2) );
            logPr.eta = logPr.eta + calcLogPrGamma( Pz, [], [], alpha0, kappa0 );
            logPr.z = logPr.z + sum(sum( Zstats.obj(ii).Nz .* log(Pi) ) );
        end
    otherwise
        error( 'To Do.' );
end

if ~exist( 'Xstats', 'var' )
    Xstats = getXSuffStats( Psi.F, Psi.stateSeq, data_struct, model );
end
logPr.obs = 0;
switch model.obsModel.type
    case 'Multinomial'
        lambda = model.obsModel.params.lambda';
        logp = Psi.theta.logp;
        %logPr.theta = calcLogPrDirichlet( Psi.theta.logp( featIDs,: ), lambda, 1 );
        %logPr.obs = sum( sum( Xstats.Nkv .* logp ) );
        %Xall = vertcat(  data_struct(:).obsHist );
        %Zall = horzcat( Psi.stateSeq(:).z  );
        logPr.obs = sum( sum( Xstats.Nkv .* logp ) );
        logPr.theta = 0;
        for kk = 1:size( Psi.F, 2 )
           if sum( Psi.F(:,kk) ) > 0
              logPr.theta = logPr.theta + calcLogPrDirichlet( logp(kk,:), lambda, 1 );
           end
        end
    case 'Gaussian'
        PP = model.obsModel.params;
        Mu = Psi.theta.Mu;
        invSigma = Psi.theta.invSigma;
        logPr.theta = 0;
        logPr.obs = 0;
        Xall = vertcat(  data_struct(:).obs );
        Zall = horzcat( Psi.stateSeq(:).z  );
        for kk = 1:size( Psi.F, 2 )
           Xkk = Xall( Zall==kk, : );
           logPr.obs = logPr.obs + sum( calcLogPrGaussian( Xkk, Mu(kk,:), invSigma(:,:,kk) ) );
           if sum( Psi.F(:,kk) ) > 0
              logPr.theta = logPr.theta + calcLogPrNormalInvWishart( Mu(kk,:), invSigma(:,:,kk), PP );
           end
        end
    case 'AR-Gaussian'
        PP = model.obsModel.params;
        A = Psi.theta.A;
        invSigma = Psi.theta.invSigma;
        logPr.theta = 0;
        logPr.obs = 0;
        Xall = vertcat(  data_struct(:).obs );
        XallR = vertcat(  data_struct(:).obsPrevR );
        Zall = horzcat( Psi.stateSeq(:).z  );
        for kk = 1:size( Psi.F, 2)
            keepINDS = Zall==kk;
            Xkk = Xall(   keepINDS, : );
            XkkR = XallR( keepINDS, : );
            logPr.obs = logPr.obs + sum( calcLogPrGaussian( Xkk, XkkR*A(:,:,kk)', invSigma(:,:,kk)  ) );

            if sum( Psi.F(:,kk) ) > 0
                logPr.theta = logPr.theta + calcLogPrMatrixNormalInvWishart( A(:,:, kk), invSigma(:,:,kk), PP );
            end
        end
    otherwise
        error( 'to do.' );
end

logPr.all = logPr.obs + logPr.z + logPr.F + logPr.eta + logPr.theta;

end

author: michaelchughes

repo: NPBayesHMM


text:
function [ a, ipvt, rcond ] = zgeco ( a, lda, n )

%*****************************************************************************80
%
%% ZGECO factors a complex matrix and estimates its condition.
%
%  Discussion:
%
%    If RCOND is not needed, ZGEFA is slightly faster.
%
%    To solve A*X = B, follow ZGECO by ZGESL.
%
%    To compute inverse(A)*C, follow ZGECO by ZGESL.
%
%    To compute determinant(A), follow ZGECO by ZGEDI.
%
%    To compute inverse(A), follow ZGECO by ZGEDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complexA(LDA,N), the matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N), an upper triangular matrix and the multipliers
%    which were used to obtain it.  The factorization can be written A = L*U
%    where L is a product of permutation and unit lower triangular matrices
%    and U is upper triangular.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, real RCOND, an estimate of the reciprocal condition of A.
%    For the system A*X = B, relative perturbations in A and B of size
%    EPSILON may cause relative perturbations in X of size (EPSILON/RCOND).
%    If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate
%    underflows.
%

%
%  Compute the 1-norm of A.
%
  anorm = 0.0;
  for j = 1 : n
    anorm = max ( anorm, dzasum ( n, a(1:n,j), 1 ) );
  end
%
%  Factor.
%
  [ a, ipvt, info ] = zgefa ( a, lda, n );
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and hermitian(A)*Y = E.
%
%  Hermitian(A) is the conjugate transpose of A.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where hermitian(U)*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve hermitian(U)*W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( zabs1 ( z(k) ) ~= 0.0 )
      ek = zsign1 ( ek, -z(k) );
    end

    if ( zabs1 ( a(k,k) ) < zabs1 ( ek - z(k) ) )
      s = zabs1 ( a(k,k) ) / zabs1 ( ek - z(k) );
      z(1:n) = z(1:n) * s;
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = zabs1 ( wk );
    sm = zabs1 ( wkm );

    if ( zabs1 ( a(k,k) ) ~= 0.0 )
      wk = wk / conj ( a(k,k) );
      wkm = wkm / conj ( a(k,k) );
    else
      wk = 1.0;
      wkm =1.0;
    end

    for j = k+1 : n
      sm = sm + zabs1 ( z(j) + wkm * conj ( a(k,j) ) );
      z(j) = z(j) + wk * conj ( a(k,j) );
      s = s + zabs1 ( z(j) );
    end

    if ( s < sm )
      t = wkm - wk;
      wk = wkm;
      z(k+1:n) = z(k+1:n) + t * conj ( a(k,k+1:n) );
    end

    z(k) = wk;

  end

  s = 1.0 / dzasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
%
%  Solve hermitian(L) * Y = W.
%
  for k = n : -1 : 1

    if ( k < n )
      z(k) = z(k) + z(k+1:n) * conj ( a(k+1:n,k) );
    end

    if ( 1.0 < zabs1 ( z(k) ) )
      s = 1.0 / zabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
    end

    l = ipvt(k);

    t    = z(l);
    z(l) = z(k);
    z(k) = t;

  end

  s = 1.0 / dzasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = 1.0;
%
%  Solve L * V = Y.
%
  for k = 1 : n

    l = ipvt(k);

    t    = z(l);
    z(l) = z(k);
    z(k) = t;

    if ( k < n )
      z(k+1:n) = z(k+1:n) + t * transpose ( a(k+1:n,k) );
    end

    if ( 1.0 < zabs1 ( z(k) ) )
      s = 1.0 / zabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
      ynorm = s * ynorm;
    end

  end

  s = 1.0 / dzasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;
%
%  Solve U * Z = V.
%
  for k = n : -1 : 1

    if ( zabs1 ( a(k,k) ) < zabs1 ( z(k) ) )
      s = zabs1 ( a(k,k) ) / zabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
      ynorm = s * ynorm;
    end

    if ( zabs1 ( a(k,k) ) ~= 0.0 )
      z(k) = z(k) / a(k,k);
    else
      z(k) = 1.0;
    end

    t = -z(k);
    z(1:k-1) = z(1:k-1) + t * transpose ( a(1:k-1,k) );

  end
%
%  Make ZNORM = 1.
%
  s = 1.0 / dzasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
function test04 ( dim_num, level_max )

%*****************************************************************************80
%
%% TEST04 sums the weights and compares them to 2^DIM_NUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the level.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04:\n' );
  fprintf ( 1, '  Compute the weights of a Gauss-Legendre sparse grid .\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  As a simple test, sum these weights.\n' );
  fprintf ( 1, '  They should sum to exactly 2^DIM_NUM.\n' );

  level_min = max ( 0, level_max + 1 - dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MIN = %d\n', level_min );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
%
%  Determine the number of points.
%
  point_num = sparse_grid_gl_size ( dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of unique points in the grid = %d\n', point_num );
%
%  Compute the weights and points.
%
  [ grid_weight, grid_point ] = sparse_grid_gl ( dim_num, level_max, point_num );
%
%  Sum the weights.
%
  weight_sum = sum ( grid_weight(1:point_num) );

  weight_sum_exact = 2.0^dim_num;

  weight_sum_error = abs ( weight_sum - weight_sum_exact );

  fprintf ( 1, '\n' ) ;
  fprintf ( 1, '    Weight sum     Exact sum    Difference\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %12e  %12e  %12e\n', ...
    weight_sum, weight_sum_exact, weight_sum_error );

  return
end


text:
function CartStates=stateRuv2Cart(x)
%%STATERUV2CART Convert a target state in 3D space in monostatic
%               r-u-v coordinates with first and possibly second time
%               derivatives of the position components into Cartesian
%               coordinates. The state has the format
%               [r;u;v;rDot;uDot;vDot;rDDot;uDDot;vDDot],
%               where two Ds indicate a second derivative with respect to
%               time. In Cartesian coordinates, the converted state has the
%               form [x;y;z;xDot;yDot;zDot;xDDot;yDDot;zDDot]. The target
%               is assumed to be in front of the radar.
%
%INPUTS: x The 6XN or 9XN set of r-u-v state vectors consisting of
%          position, velocity and possibly acceleration.  The range is a
%          one-way range.
%
%OUTPUTS: CartStates The 6XN or 9XN set of target states given in Cartesian
%                    coordinates.
%
%The derivation is given in [1]. The function aCVRuv is an implementation
%of a related linear dynamic model.
%
%EXAMPLE:
%Here we note that the results are consistent with the inverse function:
%stateCart2Ruv.
% x=[100;-60;200;-3;12;160;108;-116;2];
% xRet=stateRuv2Cart(stateCart2Ruv(x));
% max(abs(x(:)-xRet(:)))
%One will see that the error is less than 1e-13, indicating good agreement.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic linear Cartesian dynamic models in local
%    coordinates," Naval Research Laboratory, Washington, DC, Tech. Rep.
%    NRL/MR/5344-19-9882, 24 Aug. 2019.
%
%August 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(x,1);
N=size(x,2);
CartStates=zeros(numDim,N);

r=x(1,:);
u=x(2,:);
v=x(3,:);
rDot=x(4,:);
uDot=x(5,:);
vDot=x(6,:);

w2=max(0,1-u.^2-v.^2);
w=sqrt(w2);
diffV2=1-v.^2;
diffV=sqrt(1-v.^2);
denom2=w2.*diffV2;
denom=sqrt(denom2);

u1=[u;v;w];
u2=[w./diffV;zeros(1,N);-u./diffV];
u3=[-u.*v./diffV;diffV;-v.*(w./diffV)];

%Position
CartStates(1:3,:)=bsxfun(@times,r,u1);

c1=(uDot.*diffV2+u.*v.*vDot)./denom;
c2=vDot./diffV;

%Velocity
CartStates(4:6,:)=bsxfun(@times,rDot,u1)+bsxfun(@times,r.*c1,u2)+bsxfun(@times,r.*c2,u3);

if(numDim>6)
    %Acceleration
    rDDot=x(7,:);
    uDDot=x(8,:);
    vDDot=x(9,:);

    c3=-(((w+u.^2./w).*(uDot-uDot.*v.^2+u.*v.*vDot))./diffV.^3);
    c4=v.*(-u.^2.*(1./w)-w).*(-uDot.*diffV2-u.*v.*vDot)./diffV2.^2;
    c5=-vDot./diffV;
    c6=(w./diffV).*(-v.*uDot.*diffV2-u.*vDot)./diffV.^3+u.*(diffV./w).*(-u.*v.*uDot.*diffV2-u.^2.*vDot+vDot.*diffV2.^2)./diffV.^5;

    c1Dot=((uDot.*diffV2+u.*v.*vDot).*(v.*vDot.*(2-u.^2-2*v.^2)+u.*uDot.*diffV2))./denom.^3 ...
          +(uDDot.*diffV2-v.*uDot.*vDot+u.*(v.*vDDot+vDot.^2))./denom;
    c2Dot=(vDDot.*diffV2+v.*vDot.^2)./diffV.^3;

    a1=rDDot+r.*(c1.*c3+c2.*c5);
    a2=2*rDot.*c1+r.*(c1Dot+c2.*c6);
    a3=2*rDot.*c2+r.*(c2Dot+c1.*c4);

    CartStates(7:9,:)=bsxfun(@times,a1,u1)+bsxfun(@times,a2,u2)+bsxfun(@times,a3,u3);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%IM2COL Convert an image to pixel per row format
%
% OUT = IM2COL(IM) is a matrix (NxP) where each row represents a single
% of the image IM (HxWxP).  The pixels are in image column order (ie. column
% 1, column 2 etc) and there are N=WxH rows.
%
% OUT = IM2COL(IM, MASK) as above but only includes pixels if:
% - the corresponding element of MASK (HxW) is non-zero
% - the corresponding element of MASK (N) is non-zero where N=HxW
% - the pixel index is included in the vector MASK
%
% See also COL2IM.

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

function c = im2col(im, mask)

    if ndims(im) == 3
        c = reshape(shiftdim(im, 2), 3, [])';
    else
        c = reshape(im, 1, [])';
    end

    if nargin > 1 && ~isempty(mask)
        d = size(im);
        if ndims(mask) == 2 && all(d(1:2) == size(mask))
            k = find(mask);
        elseif isvector(mask)
            k = mask;
        else
            error('MVTB:im2col:badarg', 'mask must be same size as image or a vector');
        end

        c = c(k,:);
    end
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (kalman)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% maximize - if 1, does max-product (same as sum-product for Gaussians!), else sum-product [0]
% filter -   if 1, do filtering, else smoothing [0]
%
% e.g., engine = enter_evidence(engine, ev, 'maximize', 1)

maximize = 0;
filter = 0;

% parse optional params
args = varargin;
nargs = length(args);
if nargs > 0
  for i=1:2:nargs
    switch args{i},
     case 'maximize', maximize = args{i+1};
     case 'filter', filter = args{i+1};
     otherwise,
      error(['invalid argument name ' args{i}]);
    end
  end
end

assert(~maximize);

bnet = bnet_from_engine(engine);
n = length(bnet.intra);
onodes = bnet.observed;
hnodes = mysetdiff(1:n, onodes);
T = size(evidence, 2);
ns = bnet.node_sizes;
O = sum(ns(onodes));
data = reshape(cat(1, evidence{onodes,:}), [O T]);

A = engine.trans_mat;
C = engine.obs_mat;
Q = engine.trans_cov;
R = engine.obs_cov;
init_x = engine.init_state;
init_V = engine.init_cov;

if filter
  [x, V, VV, loglik] = kalman_filter(data, A, C, Q, R, init_x, init_V);
else
  [x, V, VV, loglik] = kalman_smoother(data, A, C, Q, R, init_x, init_V);
end

% Wrap the posterior inside a potential, so it can be marginalized easily
engine.one_slice_marginal = cell(1,T);
engine.two_slice_marginal = cell(1,T);
ns(onodes) = 0;
ns(onodes+n) = 0;
ss = length(bnet.intra);
for t=1:T
  dom = (1:n);
  engine.one_slice_marginal{t} = mpot(dom+(t-1)*ss, ns(dom), 1, x(:,t), V(:,:,t));
end
% for t=1:T-1
%   dom = (1:(2*n));
%   mu = [x(:,t); x(:,t)];
%   Sigma = [V(:,:,t) VV(:,:,t+1)';
% 	   VV(:,:,t+1) V(:,:,t+1)];
%   engine.two_slice_marginal{t} = mpot(dom+(t-1)*ss, ns(dom), 1, mu, Sigma);
% end
for t=2:T
  %dom = (1:(2*n));
  current_slice = hnodes;
  next_slice = hnodes + ss;
  dom = [current_slice next_slice];
  mu = [x(:,t-1); x(:,t)];
  Sigma = [V(:,:,t-1) VV(:,:,t)';
	   VV(:,:,t) V(:,:,t)];
  engine.two_slice_marginal{t-1} = mpot(dom+(t-2)*ss, ns(dom), 1, mu, Sigma);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function A = warmUpExercise()
% function [y1,...,yN] = myfun(x1,...,xN)
% above function 'myfun' takes argument (x1,...,xN) and returns y1,...,yN
% Return the 5x5 identity matrix in octave

A = eye(5);
end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function phasediag
% Phase diagram visualisation
%    using MATLAB expm
%
%   $Ekkehard Holzbecher  $Date: 2006/04/15 $
%--------------------------------------------------------------------------
T = 10;                % maximum time
C = [-1 1; 1 -3];      % matrix
f = [1; 0];            % input vector
cc = 1;                % initial concentrations (absolute value of the vector)
N = 60;                % discretization of time
M = 16;                % no. of trajectories

%----------------------execution & output----------------------------------
equilibrium = -(inv(C)*f);
t = linspace (0,T,N);
for angle = linspace (0,pi+pi,M)
    c0 = equilibrium + cc*[sin(angle); cos(angle)]; c = c0;
    for i = 2:N
        E = expm(C*t(i));
        c = [c E*c0-(eye(size(C,1))-E)*inv(C)*f];
    end
    plot (c(1,:)',c(2,:)'); hold on;
end

plot (equilibrium(1),equilibrium(2),'s');
xlabel ('variable 1'); ylabel ('variable 2')
title ('phase diagram')


text:
% Textual representation of any MatLab value.

% Copyright 2009 Levente Hunyadi
function s = var2str(value)

if islogical(value) || isnumeric(value)
    s = num2str(value);
elseif ischar(value) && isvector(value)
    s = reshape(value, 1, numel(value));
elseif isjava(value)
    s = char(value);  % calls java.lang.Object.toString()
else
    try
        s = char(value);
    catch %#ok<CTCH>
        s = '[no preview available]';
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [ vG ] = CalcFunGrad( vX, hObjFun, difMode, epsVal )
% ----------------------------------------------------------------------------------------------- %
% [ vG ] = CalcFunGrad( vX, hObjFun, difMode, epsVal )
%   Calculating the Gradient Vector of a function using Finite Differences
%   method.
% Input:
%   - vX            -   Input Vector.
%                       The point the gradient is caclaulated at.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - hObjFun       -   Objective Function.
%                       Function handler which evaluates the Objective
%                       Function at a given point - hObjFun(vX).
%                       Structure: Function Handler.
%                       Type: Function Handler.
%                       Range: NA.
%   - difMode       -   Difference Mode.
%                       Sets the mode of operation of the Finite
%                       Differences Method - Forward, Backward, Central or
%                       Complex.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, 3, 4}.
% Output:
%   - vG            -   Gradient Vector.
%                       The numerical approximation of the gradient of the
%                       Objective Function at the input point 'vX'.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  https://en.wikipedia.org/wiki/Finite_difference_coefficient.
%   2.  Nick Higham Talk (https://www.youtube.com/watch?v=Q9OLOqEhc64) At 02:15.
%   3.  Complex Step Differentiation - https://blogs.mathworks.com/cleve/2013/10/14/complex-step-differentiation/.
%   4.  Complex Step Derivative - https://timvieira.github.io/blog/post/2014/08/07/complex-step-derivative/.
%   5.  The Complex Step Derivative Approximation - https://dl.acm.org/citation.cfm?id=838251.
% Remarks:
%   1.  If the Complex Mode is selected the function must return complex
%       values in order to work. For instance, if the input function is
%       'norm(vX)' use 'sum(vX .^ 2)' and if the input function is
%       'sum(abs(vX))' use 'sum(sqrt(vX .^ 2))'.
% TODO:
%   1.  U.
% Release Notes:
%   -   1.1.000     22/08/2018  Royi Avital
%       *   Added Complex Mode.
%   -   1.0.001     13/06/2018  Royi Avital
%       *   Small fixes for better readability.
%   -   1.0.000     08/07/2017  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

DIFF_MODE_FORWARD   = 1;
DIFF_MODE_BACKWARD  = 2;
DIFF_MODE_CENTRAL   = 3;
DIFF_MODE_COMPLEX   = 4;

numElements = size(vX, 1);

objFunRef   = hObjFun(vX);
vG          = zeros([numElements, 1]);
vPertVal    = zeros([numElements, 1]);

switch(difMode)
    case(DIFF_MODE_FORWARD)
        hCalcGradFun = @(vPertVal) (hObjFun(vX + vPertVal) - objFunRef) / epsVal;
    case(DIFF_MODE_BACKWARD)
        hCalcGradFun = @(vPertVal) (objFunRef - hObjFun(vX - vPertVal)) / epsVal;
    case(DIFF_MODE_CENTRAL)
        hCalcGradFun = @(vPertVal) (hObjFun(vX + vPertVal) - hObjFun(vX - vPertVal)) / (2 * epsVal);
    case(DIFF_MODE_COMPLEX)
        hCalcGradFun = @(vPertVal) imag(hObjFun(vX + (1i * vPertVal))) / epsVal;
end

for ii = 1:numElements
    vPertVal(ii)    = epsVal;
    vG(ii)          = hCalcGradFun(vPertVal);
    vPertVal(ii)    = 0;
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [ x, seed ] = triangle01_sample ( n, seed )

%*****************************************************************************80
%
%% TRIANGLE01_SAMPLE samples the interior of the unit triangle in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Reuven Rubinstein,
%    Monte Carlo Optimization, Simulation, and Sensitivity
%    of Queueing Networks,
%    Krieger, 1992,
%    ISBN: 0894647644,
%    LC: QA298.R79.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real X(2,N), the points.
%
  m = 2;

  for j = 1 : n

    [ e, seed ] = r8vec_uniform_01 ( m + 1, seed );

    e(1:m+1) = - log ( e(1:m+1) );

    x(1:m,j) = e(1:m) / sum ( e(1:m+1) );

  end

  return
end


text:
function triangulation_test026 ( )

%*****************************************************************************80
%
%% TEST026 tests DIAEDG.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2009
%
%  Author:
%
%    John Burkardt
%
  node_num = 4;
  triangle_num = 2;
  triangle_order = 3;
  seed = 123456789;
  test_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST026\n' );
  fprintf ( 1, '  DIAEDG determines whether two triangles\n' );
  fprintf ( 1, '  with a common edge need to "swap" diagonals.\n' );
  fprintf ( 1, '  If swapping is indicated, then ALPHA_MIN should increase.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Swap   ALPHA_MIN   ALPHA_MIN\n' );
  fprintf ( 1, '         Unswapped   Swapped\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
%
%  Generate a random quadrilateral (1,2,3,4).
%
    [ node_xy, seed ] = quad_convex_random ( seed );
%
%  Does it need swapping?
%
    value = diaedg ( ...
      node_xy(1,1), node_xy(2,1), ...
      node_xy(1,2), node_xy(2,2), ...
      node_xy(1,3), node_xy(2,3), ...
      node_xy(1,4), node_xy(2,4) );

    if ( value == 1 )
      swap = 0;
    else
      swap = 1;
    end
%
%  Compute ALPHA_MIN unswapped.
%
    triangle_node(1:3,1) = [ 1, 2, 3 ]';
    triangle_node(1:3,2) = [ 1, 3, 4 ]';

    [ alpha_min_unswapped, alpha_ave, alpha_area ] = alpha_measure ( ...
      node_num, node_xy, triangle_order, triangle_num, triangle_node );
%
%  Compute ALPHA_MIN swapped.
%
    triangle_node(1:3,1) = [ 1, 2, 4 ]';
    triangle_node(1:3,2) = [ 2, 3, 4 ]';

    [ alpha_min_swapped, alpha_ave, alpha_area ] = alpha_measure ( ...
      node_num, node_xy, triangle_order, triangle_num, triangle_node );

    if ( 0 )
      r8mat_transpose_print ( 2, node_num, node_xy, '  Quadrilateral' );
    end

    fprintf ( 1, '     %1d  %10f  %10f\n', ...
      swap, alpha_min_unswapped, alpha_min_swapped );

  end

  return
end


text:
function mObject3 = directSTRAIGHTmorphing(mObject1,mObject2,mRate,mixMethod);
%   Morphing based on direct mixing of STRAIGHT parameters
%   (Without time alignment)
%   mObject3 = directSTRAIGHTmorphing(mObject1,mObject2,mRate);

%   Designed and coded by Hideki Kawahara
%   27/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara

if mObject1.samplingFrequency ~= mObject2.samplingFrequency
    mObject3 = [];
    return
end;
if mObject1.frameUpdateInterval ~= mObject2.frameUpdateInterval
    mObject3 = [];
    return
end;

nw1 = length(mObject1.F0);
nw2 = length(mObject2.F0);
[nr1,nc1] = size(mObject1.spectrogram);
[nr2,nc2] = size(mObject2.spectrogram);
nr3 = max(nr1,nr2);
nc3 = max(max(nc1,nc2),max(nw1,nw2));
nsg = zeros(nr3,nc3);
nSgram = zeros(nr3,nc3);
ap = zeros(nr3,nc3);
f0 = zeros(nc3);
nVoice = zeros(nc3);

switch mixMethod
    case 'linear'
        nsg(1:nr1,1:nc1) = (1-mRate)*mObject1.spectrogram;
        nSgram(1:nr1,1:nc1) = nSgram(1:nr1,1:nc1)+(1-mRate);
        nsg(1:nr2,1:nc2) = mRate*mObject2.spectrogram+nsg(1:nr2,1:nc2);
        nSgram(1:nr2,1:nc2) = nSgram(1:nr2,1:nc2)+mRate;
        nsg = nsg./nSgram;
    case 'log'
        nsg(1:nr1,1:nc1) = (1-mRate)*log(mObject1.spectrogram);
        nSgram(1:nr1,1:nc1) = nSgram(1:nr1,1:nc1)+(1-mRate);
        nsg(1:nr2,1:nc2) = mRate*log(mObject2.spectrogram)+nsg(1:nr2,1:nc2);
        nSgram(1:nr2,1:nc2) = nSgram(1:nr2,1:nc2)+mRate;
        nsg = exp(nsg./nSgram);
end;
ap(1:nr1,1:nc1) = (1-mRate)*mObject1.aperiodicityIndex;
ap(1:nr2,1:nc2) = mRate*mObject2.aperiodicityIndex+ap(1:nr2,1:nc2);

f0(mObject1.F0>0) = (1-mRate)*log(mObject1.F0(mObject1.F0>0));
nVoice(mObject1.F0>0) = nVoice(mObject1.F0>0)+(1-mRate);
f0(mObject2.F0>0) = mRate*log(mObject2.F0(mObject2.F0>0))+f0(mObject2.F0>0);
nVoice(mObject2.F0>0) = nVoice(mObject2.F0>0)+mRate;
f0(nVoice>0) = exp(f0(nVoice>0)./nVoice(nVoice>0));

mObject3 = createMobject;
mObject3 = updateFieldOfMobject(mObject3,'spectrogram',nsg);
mObject3 = updateFieldOfMobject(mObject3,'aperiodicityIndex',ap);
mObject3 = updateFieldOfMobject(mObject3,'F0',f0);

author: HidekiKawahara


text:
%%% Launch algorithms
% results = run_algorithm(method_id, algorithm_id, A, params)
%
% method_id - string
% algorithm_id - string
% A - matrix/tensor
% params - struct (optional)
%
% params.Idx - observed indexes (vector)
% params.Omega - observed entries (binary matrix/tensor)
%
% results - struct
% results.L - low-rank component
% results.S - sparse component
% results.O - hard thresholding
% results.cputime - CPU time;
%
function results = run_algorithm(method_id, algorithm_id, data, params)
  if(nargin < 3)
    error('method, algorithm and data must be defined');
  end
  if(nargin < 4)
    params = [];
  end
  lrs_load_conf;

  method_id = strtrim(method_id);
  algorithm_id = strtrim(algorithm_id);

  method_name = get_method_name_by_id(method_id);
  algorithm_name = get_algorithm_name_by_id(algorithm_id);
  displog(['Running ' method_name ' with ' algorithm_name]);

  switch method_id
		case 'RPCA'
			method_path = lrs_conf.rpca_path;
		case 'ST'
			method_path = lrs_conf.st_path;
		case 'MC'
			method_path = lrs_conf.mc_path;
		case 'LRR'
			method_path = lrs_conf.lrr_path;
		case 'TTD'
			method_path = lrs_conf.ttd_path;
		case 'NMF'
			method_path = lrs_conf.nmf_path;
		case 'NTF'
			method_path = lrs_conf.ntf_path;
		case 'TD'
			method_path = lrs_conf.td_path;
		otherwise
			error('Undefined method!');
	end

  alg_path = fullfile(method_path,algorithm_id);
  addpath(genpath(alg_path));

  params.A = data; M = data; T = data;
  L = zeros(size(data)); % low-rank component
  S = zeros(size(data)); % sparse component
  results.cputime = 0;

  if(~isfield(params,'rows') && ~isfield(params,'cols'))
    params.rows = size(data,1);
    params.cols = size(data,2);
  end

  % For matrix/tensor completion
  if(~isfield(params,'Idx') && ~isfield(params,'Omega'))
    [params.Idx, params.Omega] = subsampling(data, 0.5);
  end
  if(isfield(params,'Idx') && ~isfield(params,'Omega')) % Build "params.Omega" from "params.Idx"
    params.Omega = zeros(size(data));
    params.Omega(params.Idx) = 1;
  end
  if(isfield(params,'Omega')) % Build "params.Idx" from "params.Omega"
    % params.Omega = ones(size(A));
    % params.Omega = randi([0 1],size(A));
    params.Idx = find(params.Omega);
  end
  Idx = params.Idx; % MIdx = M(Idx);
  Omega = params.Omega; % MOmega = M.*Omega;

  timerVal = tic;
  % warning('off','all');
  run_alg;
  % warning('on','all');
  cputime = toc(timerVal);

  rmpath(genpath(alg_path));

  results.L = L; % low-rank component
  results.S = S; % sparse component
  results.O = hard_threshold(S); % apply hard thresholding
  results.Omega = params.Omega; % for matrix/tensor completion
  results.Idx = params.Idx; % for matrix/tensor completion
  results.cputime = cputime;
	displog('Decomposition finished');
end

author: andrewssobral

repo: lrslibrary


text:
function val=StirlingNumber2(n,k)
%%STIRLINGNUMBER2 Compute the Stirling number of the second kind {n,k},
%                 which is also called a Stirling set number. This is the
%                 number of ways to partition a set of n objects into k
%                 non-empty subsets. The number is computed using a
%                 recursion that will not overflow unless the actual number
%                 overflows. When n is very large, the function can be
%                 slow.
%
%INPUTS: n The integer total number of items in the set n>=0.
%        k The desired integer number of non-empty subsets of the set.
%          k>=0.
%
%OUTPUTS: val The number of ways of partitioning n items into k subsets. If
%             an overflow occurs, this will be infinite.
%
%The implementation uses the recurrence relation from [1].
%
%The more commonly seen expression in terms of a sum of binomials is not
%desirable, because overflows will occur even when the final number does
%not overflow.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Stirling Number of the Second Kind." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/StirlingNumberoftheSecondKind.html
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(n<k)
    val=0;
    return;
end

curN=zeros(k+1,1);
prevN=zeros(k+1,1);

prevN(0+1)=1;%{0,0}=1
for nCur=1:n
    for kCur=1:min(nCur,k)
        curN(kCur+1)=prevN(kCur-1+1)+kCur*prevN(kCur+1);
    end
    prevN=curN;
end
val=prevN(k+1);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
% test for image warping
M = load_image('barb');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% extract eyes %%%
clf;
imagesc([0,1],[0,1],M);
title('Click on the each eye');
colormap gray(256);
axis image; axis off;
[x,y,b] = ginput(2);

% position of the eye
u = [x(1) y(1)];
v = [x(2) y(2)];

u1 = [0.3,0.3];
v1 = [0.7,0.3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% computations %%%
M1 = perform_image_similitude(M,u,u1,v,v1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% display %%%
clf;
subplot(1,2,1);
hold on;
imagesc([0,1],[0,1],M);
plot([u(1) v(1)], [u(2) v(2)], '*');
hold off;
title('Original');
colormap gray(256);
axis image; axis off; axis ij

subplot(1,2,2);
hold on;
imagesc([0,1],[0,1],M1);
plot([u1(1) v1(1)], [u1(2) v1(2)], '*');
hold on;
title('Warped');
colormap gray(256);
axis image; axis off; axis ij

author: gpeyre

repo: matlab-toolboxes


text:
function planC = flipAlongX(scanNum,planC)
% function planC = flipAlongX(scanNum,planC)
%
% Flip scan, structures and doses along Y direction
%
% USAGE:
% global planC
% scanNum = 1;
% planC = flipAlongY(scanNum)
%
% APA, 05/17/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS
if ~exist('planC','var')
    global planC
end
indexS = planC{end};

[xVals, yVals, zVals] = getScanXYZVals(planC{indexS.scan}(scanNum));

xMin = min(xVals);
xMax = max(xVals);

[assocScansV, relStructNumV] = getStructureAssociatedScan(1:length(planC{indexS.structures}),planC);

strV = find(assocScansV == scanNum);

%% Flip structures
for strNum = 1:length(strV)
    for slcNum = 1:length(planC{indexS.structures}(strV(strNum)).contour)
        for segNum = 1:length(planC{indexS.structures}(strV(strNum)).contour(slcNum).segments)
            pointsM = planC{indexS.structures}(strV(strNum)).contour(slcNum).segments(segNum).points;
            if ~isempty(pointsM)
                pointsM(:,1) = xMin + xMax - pointsM(:,1);
                planC{indexS.structures}(strV(strNum)).contour(slcNum).segments(segNum).points = pointsM;
            end
        end
    end
end

%% Flip Scan
planC{indexS.scan}(scanNum).scanArray = flipdim(planC{indexS.scan}(scanNum).scanArray,2);

%% Flip Dose
for doseNum = 1:length(planC{indexS.dose})
    assocScanNum = getAssociatedScan(planC{indexS.dose}(doseNum).assocScanUID,planC);
    if scanNum == assocScanNum
        planC{indexS.dose}(doseNum).doseArray = flipdim(planC{indexS.dose}(doseNum).doseArray,2);
        xV = getDoseXYZVals(planC{indexS.dose}(doseNum));
        planC{indexS.dose}(doseNum).coord1OFFirstPoint = ...
           2*planC{indexS.scan}(scanNum).scanInfo(1).xOffset - xV(end);
    end
end

%ReRaster and ReUniformize
planC = reRasterAndUniformize(planC,scanNum);

if isfield(stateS,'handle') && isfield(stateS.handle,'CERRSliceViewer') && ...
        ishandle(stateS.handle.CERRSliceViewer) && ...
        isfield(stateS,'planLoaded') && stateS.planLoaded
    CERRRefresh
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function mat = matrix(rot)
% quaternion to direction cosine matrix conversion
% converts direction cosine matrix to quaternion
%
% Syntax
%   mat = matrix(q)
%
% Input
%
%  q - @quaternion
%
% Output
%
%  mat - vector of matrixes
%
% See also
% mat2quat Euler axis2quat hr2quat

mat = matrix@quaternion(rot);

mat = repmat(reshape(1-2*rot.i,1,1,[]),3,3) .* mat;

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function MRIVertices = ctf_head2mri(HeadVertices,mri)

% ctf_head2mri - convert a CTF head into a voxel coordinate
%
% MRIVertices = ctf_mri2head(HeadVertices,mri)
%
% This function converts CTF head coordinates into the CTF MRI voxel
% coordinates.  The input 'HeadVertices' are head space locations; they are
% Nx3 coordinates in the CTF head space (see below). The 'mri' input
% is a struct returned by ctf_read_mri, which contains the fiducials
% and a transformation matrix between head and MRI voxel coordinates.
%
% CTF MRI volumes are 256^3 voxels, 1mm^3 each. The volume index has
% an origin at the left, anterior, superior voxel, such that:
%
% Sag increases from left to right (+X Right)
% Cor increases from anterior to posterior (+Y Posterior)
% Axi increases from superior to inferior (+Z Inferior)
%
% The head space coordinates are defined in relation to the MRI fiducial
% locations (nasion, left preauricular and right preauricular).  The origin
% lies half way between the left and right preauricular points and the
% coordinate axes are given as:
%
% +X is through the nasion,
% +Y is left
% +Z is superior
%
% CTF head coordinate axes are othogonalized by first taking the cross
% product of the +X and +Y vectors (from the origin through the nasion and
% left preauricular, respectively) to find the +Z vector.  Then, the Y axis
% is orthogonalized to the X/Z plane using the cross product and right hand
% rule.  Hence, +Y can be slightly offset from the left preauricular.
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $

% Copyright (C) 2004  Darren L. Weber
%
% 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 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% History:  04/2004, Darren.Weber_at_radiology.ucsf.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('in development!'); return

ver = '$Revision: 1.1 $';
fprintf('CTF_HEAD2MRI [v %s]\n',ver(11:15)); tic;

%--------------------------------------------------------------------
% Extract the coordinate transform matrices from the mri struct.  These are
% designed to be left multiplied into the vertices (I don't like this
% because it means the column arrays need to be transposed below)
%T.ctfHEAD2MRI = mri.hdr.transformMatrixHead2MRI;
%T.ctfMRI2HEAD = mri.hdr.transformMatrixMRI2Head;

% OK, so I got fed up with trying to figure out how to use these
% bloody matrices and decided to calculate them myself, so
% that the entries would work in a way that I can understand!
[trans,rot] = calc_tranfer_matrix(mri);

%--------------------------------------------------------------------
fprintf('...converting from cm to mm\n');
MRIVertices = HeadVertices * 10;

%--------------------------------------------------------------------
% Pad out the HeadVertices to Nx4 matrix

Nvertices = size(MRIVertices,1);

right_column = ones( Nvertices, 1 );

MRIVertices = [ MRIVertices right_column ];

%--------------------------------------------------------------------
% Convert CTF head space coordinates into voxel indices

MRIVertices = (-1 * MRIVertices * trans) * rot;

MRIVertices = MRIVertices(:,1:3);

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

%-------------------------------------------------------
function [trans,rot] = calc_tranfer_matrix(mri)

% these are the translations!
%ctf.mri.hdr.headOrigin_sagittal: 130
%ctf.mri.hdr.headOrigin_coronal: 123
%ctf.mri.hdr.headOrigin_axial: 152

% This is how they are calculated from the fiducials:

nas(1) = mri.hdr.HeadModel_Info.Nasion_Sag;
nas(2) = mri.hdr.HeadModel_Info.Nasion_Cor;
nas(3) = mri.hdr.HeadModel_Info.Nasion_Axi;
lpa(1) = mri.hdr.HeadModel_Info.LeftEar_Sag;
lpa(2) = mri.hdr.HeadModel_Info.LeftEar_Cor;
lpa(3) = mri.hdr.HeadModel_Info.LeftEar_Axi;
rpa(1) = mri.hdr.HeadModel_Info.RightEar_Sag;
rpa(2) = mri.hdr.HeadModel_Info.RightEar_Cor;
rpa(3) = mri.hdr.HeadModel_Info.RightEar_Axi;

headOffset_Sag = ( rpa(1) - lpa(1) ) / 2;
headOrigin_Sag = lpa(1) + headOffset_Sag;

headOffset_Cor = ( rpa(2) - lpa(2) ) / 2;
headOrigin_Cor = lpa(2) + headOffset_Cor;

headOffset_Axi = ( rpa(3) - lpa(3) ) / 2;
headOrigin_Axi = lpa(3) + headOffset_Axi;

headOrigin = [headOrigin_Sag headOrigin_Cor headOrigin_Axi];

%-------------from here to...

% calculate voxel space vectors, in slices
voxNASvector = nas - headOrigin;
voxLPAvector = lpa - headOrigin;
voxRPAvector = rpa - headOrigin;

% calculate head space vectors, in mm
headNASvector = -1 * [ voxNASvector(2), voxNASvector(1), voxNASvector(3) ];
headLPAvector = -1 * [ voxLPAvector(2), voxLPAvector(1), voxLPAvector(3) ];
headRPAvector = -1 * [ voxRPAvector(2), voxRPAvector(1), voxRPAvector(3) ];

%-------------here; is encapsulated in trans:

trans = eye(3);
trans(1,1) = 0;
trans(2,1) = 1;
trans(1,2) = 1;
trans(2,2) = 0;
trans(4,1:3) = -1 * [ headOrigin(2) headOrigin(1) headOrigin(3) ];

%headVector = -1 * ctfVox * trans;

% At this point, the voxel indices are now in the head space, in
% mm; however, the head space at this point is not orthogonal

%--------------now the head space is orthogonalized

headNASunit = headNASvector / norm(headNASvector);
headLPAunit = headLPAvector / norm(headLPAvector);
headRPAunit = headRPAvector / norm(headRPAvector);

headVERTEXunit = cross( headNASunit, headLPAunit );

% revise the LPA unit vector, so it is orthogonal to Nasion
headLPAunit = cross( headVERTEXunit, headNASunit );

% CHECK, these dot products = 0
%dot( headNASunit, headLPAunit )
%dot( headNASunit, headVERTEXunit )
%dot( headLPAunit, headVERTEXunit )

% Note that the LPA/RPA moves!  This has the effect of
% rotating the coordinates in the XY plane.

rot = [ headNASunit; headLPAunit; headVERTEXunit ];

return

author: PatternRecognition

repo: OpenBMI


text:
% visualize - Visualize
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
% Original code written by Tapani Raiko
%
% 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 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function visualize(W, horiz, do_sort, borders)
if (nargin < 2)
    horiz = 0;
end

if (nargin < 3)
    do_sort = 0;
end

if (do_sort == 1)
    Wnorms = sum(W.^2,1);
    [B,IX] = sort(Wnorms,'descend');
    W = W(:,IX);
end

% how many pixels for borders?
if nargin < 4
    borders = 1;
end
%fprintf('Visualizing patches\n');
[ndim,nunits]=size(W);
npix = floor(sqrt(ndim)+0.999);
npix2 = floor(sqrt(nunits)+0.999);
minW=min(W(:));
maxW=max(W(:));
bigpic = -(minW+maxW)/2*ones(((npix+borders)*npix2+borders));
if (nunits/npix2<=npix2-1),
    bigpic = bigpic(:,1:(npix+borders)*(npix2-1)+borders);
end;
for i=1:nunits;
    if (horiz)
    bigpic(mod(i-1,npix2)*(npix+borders)+borders+1:mod(i-1,npix2)*(npix+borders)+borders+npix,...
        floor((i-1)/npix2)*(npix+borders)+borders+1:floor((i-1)/npix2)*(npix+borders)+borders+npix)...
        = reshape(W(:,i),npix,npix)';
    else
    bigpic(mod(i-1,npix2)*(npix+borders)+borders+1:mod(i-1,npix2)*(npix+borders)+borders+npix,...
        floor((i-1)/npix2)*(npix+borders)+borders+1:floor((i-1)/npix2)*(npix+borders)+borders+npix)...
        = reshape(W(:,i),npix,npix);
    end
end;
imagesc(bigpic);
colormap(gray);
axis off;
axis equal;
%fprintf('done.\n');

author: kyunghyuncho

repo: deepmat


text:
function res = imMergeCells(img, c1, c2)
%IMMERGECELLS merge labeled cell of an image
%
%   Deprecated: use 'imMergeLabels' instead.
%
%   Usage:
%   LBL2 = imMergeCells(LBL, C1, C2);
%   LBL is a label image, C1 and C2 are labels of the 2 cells to merge.
%   The indices of two regions are set to the same value, and the boundary
%   between the 2 regions is also merged to the new region, preserving
%   topology of adjacent regions.
%
%
%   -----
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 19/07/2004.
%

warning('imael:deprecated', ...
    'imMergeCells is deprecated, use imMergeLabels instead');

ind = imBoundaryIndices(img, c1, c2);
res = img;
res(img==c2)=c1;
res(ind)=c1;

return;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function varargout = levi13(X)
% Levi function, #13
%
%   LEVI13([x1, x2]) returns the value of the value of the 13th Levi
%   function at the specified points. [x1] and [x2] may be vectors.
%   The search domain is
%
%               -10 < x_i < 10
%
%   The global minimum is
%
%               f(x1, x2) = f(1, 1) = 0.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;  % # dims
        varargout{2} = [-10, -10]; % LB
        varargout{3} = [+10, +10]; % UB
        varargout{4} = [1,1]; % solution
        varargout{5} = 0; % function value at solution

    % otherwise, output function value
    else

        % keep values in teh search domain
        X(X < -10) = inf;       X(X > 10) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = sin(3*pi*x1).^2 + (x1-1).^2.*(1 + sin(3*pi*x2).^2) + ...
            (x2-1).^2.*(1 + sin(2*pi*x2).^2);

    end

end

author: Sable


text:
cur_dir = pwd;
cd(fileparts(mfilename('fullpath')));

fprintf('Downloading SPP_net_release1_VOC2007_selective_search_data...\n');
urlwrite('https://onedrive.live.com/download?resid=D7AF52BADBA8A4BC!112&authkey=!ALMzBRGL1B9Eqz8&ithint=file%2czip', ...
    'SPP_net_release1_VOC2007_selective_search_data.zip');

fprintf('Unzipping...\n');
unzip('SPP_net_release1_VOC2007_selective_search_data.zip', '.');

fprintf('Done.\n');
system('del SPP_net_release1_VOC2007_selective_search_data.zip');

cd(cur_dir);

author: ShaoqingRen


text:
% Chapter 4 - Electromagnetic Waves and Optical Resonators.
% Program_4a - Iteration of the Ikeda Map.
% Copyright Birkhauser 2013. Stephen Lynch.

% Chaotic attractor for the Ikeda map (Figure 4.11(b)).
clear
echo off
A=10;
B=0.15;
N=10000;
E=zeros(1,N);x=zeros(1,N);y=zeros(1,N);
E(1)=A;x(1)=A;y(1)=0;
for n=1:N
    E(n+1)=A+B*E(n)*exp(1i*abs(E(n))^2);
    x(n+1)=real(E(n+1));
    y(n+1)=imag(E(n+1));
end
axis([8 12 -2 2])
axis equal
plot(x(50:N),y(50:N),'.','MarkerSize',1);
fsize=15;
set(gca,'XTick',8:1:12,'FontSize',fsize)
set(gca,'YTick',-2:1:2,'FontSize',fsize)
xlabel('Real E','FontSize',fsize)
ylabel('Imag E','FontSize',fsize)

% End of Program_4a.


text:
function decision = cutBlocks(structIn,cnst,currImName)
    blockIn = structIn.data;

    currSize = size(blockIn);
    if (currSize(1) == currSize(2)) % discard non-square blocks

    % keep only non-background blocks, according to Coudray et al., Nat Med 2018
    % criterion
    bwIm = mean(blockIn,3); % black and white image
    if sum(sum(bwIm>=220))<(0.5*numel(bwIm))

    % resize block to 512x512
    blockIn = imresize(blockIn,[cnst.finalBlockSize,cnst.finalBlockSize]);

    if cnst.verbose
    imshow(blockIn);
    drawnow
    end

    % clean up file name
    currDot = strfind(currImName,'.');
    currImName = currImName(1:(currDot(1)-1));

    % identify failed (black) images
    if sum(sum(bwIm<=1))<(numel(bwIm))
    % save image
        imwrite(blockIn,[cnst.outDir,currImName,'-blk-',randseq(12,'alphabet','amino'),'.png']);
        decision = 4; % code for success
    else
        warning('moving block to failures');
        imwrite(blockIn,[cnst.failDir,currImName,'-blk-',randseq(12,'alphabet','amino'),'.png']);
        decision = 1; % code for fail
    end
    else
        if cnst.verbose
        disp('discarding background block');
        end
        decision = 2; % code for background
    end
    else
        if cnst.verbose
        disp('discarding edge');
        end
        decision = 3; % code for edge
    end

end

author: jnkather

repo: MSIfromHE

save_path: github-repos/MATLAB/jnkather-MSIfromHE


text:
function [model, lbls] = modelReadFromFile(fileName, varargin)

% MODELREADFROMFILE Read model from a file FID produced by the C++ implementation.
% FORMAT
% DESC loads in from a file a model produced by C++ code.
% C++ GP implementation.
% ARG fileName : the file ID from where the data is loaded.
% RETURN model : the model loaded in from the file.
%
% COPYRIGHT : Neil D. Lawrence, 2008
%
% SEEALSO : modelReadFromFile

% MLTOOLS

FID = fopen(fileName);
if FID==-1
  error(['Cannot find file ' fileName])
end
model = modelReadFromFID(FID, varargin{:});
fclose(FID);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function  setRelPosition( hObject,units,hObject_pos )
%this function invokes set(hObject,'Position',hObject_pos) while with 'units' provided

old_units=get(hObject,'units');
set(hObject,'units',units);
set(hObject,'Position',hObject_pos);
set(hObject,'units',old_units);

end

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_snr_20db_s_8.mat';
N = 256;
S = 8;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 20;

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function [Dps_eff, Dns_eff] = solidPhaseDiffusionCoefficients(T,param)
% solidPhaseDiffusionCoefficients evaluates diffusion coefficients of the solid phase [m^2 /s].
% The user may modify the script to meet specific requirements.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

if(param.TemperatureEnabled>=1)
    Dps_eff     = param.Dps*exp(-param.EaDps/param.R*(1./T(param.Nal+1:param.Nal+param.Np)-1/param.Tref));
else
    Dps_eff     = param.Dps*ones(param.Np,1);
end

if(param.TemperatureEnabled>=1)
    Dns_eff     = param.Dns*exp(-param.EaDns/param.R*(1./T(param.Nal+param.Np+param.Ns+1:param.Nal+param.Np+param.Ns+param.Nn)-1/param.Tref));
else
    Dns_eff     = param.Dns*ones(param.Nn,1);
end

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function display( x )
long = ~isequal(get(0,'FormatSpacing'),'compact');
if long, disp( ' ' ); end
disp(x);
if long, disp( ' ' ); end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [F,sA] = spm_multinomial_log_evidence(qA,pA,rA)
% Bayesian model reduction for multinomial distibutions
% FORMAT [F,sA] = spm_multinomial_log_evidence(qA,pA,rA)
%
% qA  - parameter of posterior of full model
% pA  - parameter of prior of full model
% rA  - parameter of prior of reduced model
%
%
% F   - (negative) free energy or log evidence of reduced model
% sA  - parameter of reduced posterior
%
% This routine computes the negative log evidence of a reduced model of a
% mutinomial distribution. This also applies for Bernoulli, Binomial, and
% Categorical distributions.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Thomas Parr
% $Id: spm_multinomial_log_evidence.m 7679 2019-10-24 15:54:07Z spm $

% reduced posteriors
%--------------------------------------------------------------------------
sA = spm_softmax(qA + rA - pA);

% change in free energy or log model evidence
%--------------------------------------------------------------------------
k = find(rA);
F  = log(rA(k(1))) - log(pA(k(1))) + log(qA(k(1))) - log(sA(k(1)));

author: spm


text:
function A = cs_frand (n,nel,s)                                             %#ok
%CS_FRAND generate a random finite-element matrix
% A = cs_frand (n,nel,s) creates an n-by-n sparse matrix consisting of nel
% finite elements, each of which are of size s-by-s with random symmetric
% nonzero pattern, plus the identity matrix.
%
% Example
%   A = cs_frand (100, 100, 3) ;
% See also cs_demo.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_frand mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function main()
    addpath('../include/vlfeat-0.9.20/toolbox/');
    vl_setup();
    addpath(genpath('../include/eccv14text'));
    %data_infos
    data_infos.img_path = '../data/msra_torch/im/';
    data_infos.map_path = '../data/msra_torch/multiscale/';
    data_infos.res_path = '../data/msra_torch/proposal_res/';

    %global variable
    global globalVar

    %gen_proposals param
    global param
    param.workPath = '/Users/zhangzheng/Documents/FCN_FULL/genProposal';
    param.reuse_mser = false;
    param.debug = false;
    param.minRegionProb = 0.2;
    param.minRegionCompCoveredArea = 0.7;
    param.maxRegionCompArea = 1;
    param.secondaryMinRegionCompArea = 0.5;
    param.orient_param.orientationInterval = 2;
    param.orient_param.minOrientation = -90;
    param.orient_param.maxOrientation = 90;
    param.minCompHeightSimilarity = 0.7;
    param.maxCompOrientationDiff = 3;
    param.minIoUDiff = 0.85;
    param.maxDistRatio = 2;
    %minmaxRFs = [137, 68, 5];
    minmaxRFs = [137, 32, 5];

    %% Only used for normal_mser3
    param.mser_info.delta = 1;
    param.mser_info.minArea = 0.002;
    param.mser_info.maxArea = 1;
    param.mser_info.minDiversity = 0.8;
    param.mser_info.maxVariation = 0.15;

    mkdir(data_infos.res_path);

    extension = '.jpg';
    imgData = dir([data_infos.img_path,'*.jpg']);% original image.
    if(length(imgData) == 0)
       imgData = dir([data_infos.img_path,'*.JPG']);% original image.
       extension = '.JPG';
    end
    nImg = length(imgData);
    for ii= 1:nImg
        disp(ii);
        [~, name, ~] = fileparts(imgData(ii).name);
%         if(~strcmp(name, 'img_11'))
%             continue;
%         end

        globalVar.imgName = name;
        img_path = [data_infos.img_path, imgData(ii).name];%the original image.

        proposalsSavePath = [data_infos.res_path, name, '.txt'];
        nMap = 3;
        proposals = cell(nMap, 1);
        for jj = 1 : nMap
            map_path = [data_infos.map_path, name, '_', num2str(jj), extension];%the res image from last phase.
            img = imread(img_path);
            map = imread(map_path);
            map = double(map) / 255;

            [map_h,map_w,~]=size(map);
            resizeRatio = size(img, 1) / map_h;
            img = imresize(img, [map_h, map_w], 'bilinear');

            proposals_tmp = gen_proposals(img, map, resizeRatio, minmaxRFs(jj));
            if(isempty(proposals_tmp) == false)
                proposals_tmp(:, 1 : 8) = proposals_tmp(:, 1 : 8) * resizeRatio;
            end
            proposals{jj} = proposals_tmp;

            if(false)
                imshow(img);
                hold on;
                for nProposal = 1 : size(proposals_tmp, 1)
                    x_arr = proposals_tmp(nProposal, 1 : 2 : 8);
                    y_arr = proposals_tmp(nProposal, 2 : 2 : 8);
                    plot([x_arr, x_arr(1)], [y_arr, y_arr(1)], 'color', rand(3,1));
                end
                hold off;
                saveas(gcf, [data_infos.res_path, name, '_', num2str(jj), '.jpg'], 'jpg');
            end
        end
        proposals = cell2mat(proposals);

        if(true)
            img = imread(img_path);
            imshow(img);
            hold on;
            for nProposal = 1 : size(proposals, 1)
                x_arr = proposals(nProposal, 1 : 2 : 8);
                y_arr = proposals(nProposal, 2 : 2 : 8);
                plot([x_arr, x_arr(1)], [y_arr, y_arr(1)], 'color', rand(3,1));
            end
            hold off;
            saveas(gcf, [data_infos.res_path, name, '_all.jpg'], 'jpg');
        end

        %% to adapter old code
        proposalsToSave = zeros(10, size(proposals, 1));
        if(size(proposalsToSave, 2) ~= 0)
            proposalsToSave(1 : 8, :) = round(proposals(:, 1 : 8))';
            proposalsToSave(10, :) = proposals(:, 9);
        end
        fid = fopen(proposalsSavePath, 'w');
        fprintf(fid, '%d %d %d %d %d %d %d %d %d %f\n', proposalsToSave);
        fclose(fid);
    end
end

author: stupidZZ


text:
function legendre_exactness ( n, x, w, p_max )

%*****************************************************************************80
%
%% LEGENDRE_EXACTNESS investigates exactness of Legendre quadrature.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points in the rule.
%
%    Input, real X(N), the quadrature points.
%
%    Input, real W(N), the quadrature weights.
%
%    Input, integer P_MAX, the maximum exponent.
%    0 <= P_MAX.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature rule for Legendre integral.\n' );
  fprintf ( 1, '  Rule of order N = %d\n', n );
  fprintf ( 1, '  Degree          Relative Error\n' );
  fprintf ( 1, '\n' );

  for p = 0 : p_max

    s = legendre_integral ( p );

    v(1:n,1) = x(1:n) .^ p;

    q = w' * v;

    if ( s == 0.0 )
      e = abs ( q - s );
    else
      e = abs ( ( q - s ) / s );
    end

    fprintf ( 1, '  %6d  %24.16f\n', p, e );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all; close all;
I=imread('onion.png');
I=rgb2gray(I);
I=im2double(I);
LEN=25;
THETA=20;
PSF=fspecial('motion', LEN, THETA);
J=imfilter(I, PSF, 'conv', 'circular');
NSR=0;
K=deconvwnr(J, PSF, NSR);
figure;
subplot(131);  imshow(I);
subplot(132);  imshow(J);
subplot(133);  imshow(K);

author: HuangCongQing


text:
function fluxlimiter(N, method)
% N = number of grid points

% grid spacing on interval (-1,1)
dx = 2./N;

% location of cell centers
x = linspace(-1+0.5*dx, 1-0.5*dx, N);

% advection speed
u = 1;

% cfl condition + safety margin
dt = 0.25*dx;

% final time
FinalTime = 2;
Nsteps = ceil(FinalTime/dt);
dt = FinalTime/Nsteps;
lam = u*dt/dx;

% initial condition
q = fluxlimiterexact(x);

% time step
for n=1:Nsteps

  qp1  = [q(2:N),q(1)];    % q_{m+1}
  qm1 = [q(N),q(1:N-1)]; % q_{m-1}
  qm2 = [qm1(N), qm1(1:N-1)]; % q_{m-2}

  % compute limited cell average jumps based on method choice
  % theta_{i-1/2} = (q_{i-1}-q_{i-2})/(q_{i}-q_{i-1})
  tol = 1e-10;
  dq1 = q-qm1;
  dq2 = qm1-qm2;
  % have to test in case dq1 is small because of finite precision effects

  % are both dq1, and dq2 small (smooth)
  ids = find(abs(dq1)+abs(dq2)<tol); % both jumps small
  thetaL(ids) = 1;

  % is dq1 small and dq2 relatively large
  ids = find(abs(dq1)<tol & abs(dq2)>tol );
  thetaL(ids) = 100; % choose theta large

  % hopefully the remainder
  ids = find(abs(dq1)>= tol);
  thetaL(ids) = dq2(ids)./dq1(ids);

  if(strcmp(method,'minmod'))
      phiL = minmod(1, thetaL);
  elseif(strcmp(method,'MC'))
      phiL = max(0, min(min(0.5*(1+thetaL),2), 2*thetaL));
  elseif(strcmp(method,'vanleer'))
      phiL = (thetaL+abs(thetaL))./(1+abs(thetaL));
  elseif(strcmp(method,'superbee'))
      phiL = max(0, max(min(1,2*thetaL), min(2,thetaL)));
  elseif(strcmp(method,'Fromm'))
      phiL = 0.5*(1+thetaL);
  elseif(strcmp(method,'BeamWarming'))
      phiL = thetaL; % upwind
  elseif(strcmp(method,'LaxWendroff'))
      phiL = ones(size(q)); % downwind
  elseif(strcmp(method, 'upwind'))
      phiL = zeros(size(q));
  elseif(strcmp(method, 'custom'))

      phiL = customphi(thetaL);
  else
      disp('WARNING: wrong choice of flux limiter');
  end
  phiR = [phiL(2:N), phiL(1)];

  % update cell averages
  q = q - lam*(q-qm1) - 0.5*lam*(1-lam)*(phiR.*(qp1-q) -phiL.*(q-qm1));

  if(mod(n,40)==0)
    plot(x,q, '.');
    hold on; plot(x-0.5*dx, phiL, 'g-'); hold off;pause(0.02);
  end
end

plot(x, q, 'b.');
% have gone one exact period
hold on; plot(x, fluxlimiterexact(x), 'r-'); hold off;
hold on; plot(x-0.5*dx, phiL, 'g-'); hold off;
xlabel('x');
legend(sprintf('%s N=%d T=%f', method, N, FinalTime), ...
           'Exact solution', sprintf('flux limiter function: phi_{%s}', method));
axis([-1 1 -1 3])

repo: Aero-matlab


text:
function [output] = compareBinsOfFluxes(xglc, model, sammin, sammax, metabolites)
% Takes the overall sammin and sammax samples, bins them into
% separate bin sizes and compares them, then compares the
% results to the largest bin size.
% calls `[totalz, zscore, mdv1, mdv2] = compareTwoSamp(xglc, model, samp1, samp2, measuredMetabolites)`
% sammin and sammax each contain bins of fluxes in `x.samps(r,1).points`
%
% USAGE:
%
%    [output] = compareBinsOfFluxes(xglc, model, sammin, sammax, metabolites)
%
% INPUTS:
%    xglc:           sugar distribution
%    model:          model structure
%    sammin:         samples containing bins of fluxes
%    sammax:         samples containing bins of fluxes
%
% OPTIONAL INPUT:
%    metabolites:    list of metabolites
%
% OUTPUT:
%    output:         result of comparison
%
% .. Author: - Wing Choi 3/7/08

output = 0;

if (nargin < 4)
    disp '[output] = compareBinsOfFluxes(xglc,model,samplo,samphi,metabolites)'
    return;
end

if (nargin < 5)
    metabolites = [];
end

if (isempty(xglc))
    % random glucose
    %xglc = rand(64,1);
    %xglc = xglc/sum(xglc);
    %xglc = idv2cdv(6)*xglc;

    glucose = rand(8,1);
    glucose = glucose/sum(glucose);
    %glc = idv2cdv(6)*glc;

    % glc 1-6 = carbon 1-6
    % glc 7 = carbon 1+2 (really 5 and 6)
    % glc 8 = unlabeled
    % glc 9 = fully labeled
    glc = zeros(64,9);
    glc(1+1,1) = 1;
    glc(2+1,2) = 1;
    glc(4+1,3) = 1;
    glc(8+1,4) = 1;
    glc(16+1,5) = 1;
    glc(32+1,6) = 1;
    glc(32+16+1,7) = 1;
    glc(0+1,8) = 1;
    glc(63+1,9) = 1;

    xGlc = zeros(64,1);
    for i = 1:8
        xGlc = xGlc + glucose(i)*glc(:,i);
    end

    xglc = idv2cdv(6)*xGlc;

end

nbins = size(sammin.samps,1);
npoints = size(sammin.samps(1,1).points,2);

disp (sprintf('found %d samples in input',npoints));
disp (sprintf('numbins  : %d',nbins));

for bin = 1:nbins
    samp1.points = sammin.samps(bin,1).points;
    samp2.points = sammax.samps(bin,1).points;
    [totalz,zscore,mdv1,mdv2] = compareTwoSamp(xglc,model,samp1,samp2,metabolites);
    output.totalz(bin,1) = totalz;
end

return
end

author: opencobra

repo: cobratoolbox


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018           |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtHmxFemBemEFIEhalf.m                        |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Francois Alouges            |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Solve fem/bem with PEC scatering problem      |
%|  `---'  |                Volumic ans surfacic EFIE formulation         |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

%%% PREPARATION
disp('~~~~~~~~~~~~~ PREPARATION ~~~~~~~~~~~~~')

% Spherical mesh
mesh = msh('sphere_1e3.msh');

% Normalization interior sphere
mesh.vtx = mesh.vtx/0.8;

% Interior boundary
bound = mesh.bnd;
Iint  = (sum(bound.ctr.*bound.nrm,2) <= 0);
int   = swap(bound.sub(Iint));

% Final volumn
ctr = mesh.ctr;
vol = mesh.sub(ctr(:,3)>=0);

% Boundaries
tmp  = vol.bnd;
tmp1 = setdiff(int,tmp);
dvol = setdiff(tmp,int);
ext  = union(tmp1,dvol);
int  = intersect(int,tmp);

figure
plot(ext)
hold on
plotNrm(ext,'w')
plot(dvol,'b')
plot(int,'r')
axis equal
alpha(0.5)

% Domaine volumique
omega = dom(vol,4);

% Domaines surfacic
sigma = dom(dvol,3);
gamma = dom(ext,3);

% Frequency adjusted to maximum edge size
stp = mesh.stp;
k   = 1/stp(2);
c   = 299792458;
f   = (k*c)/(2*pi);
disp(['Frequency : ',num2str(f/1e6),' MHz']);

% Accuracy
tol = 1e-3;
disp(['Accuracy : ',num2str(tol)]);

% Incident direction and field
X0 = [0 0 -1];
E  = [0 1  0]; % Polarization (+x for Theta-Theta and +y for Phi-Phi)
H  = cross(X0,E);

% Incident Plane wave (electromagnetic field)
PWE{1} = @(X) exp(1i*k*X*X0') * E(1);
PWE{2} = @(X) exp(1i*k*X*X0') * E(2);
PWE{3} = @(X) exp(1i*k*X*X0') * E(3);

% Incident wave representation
figure
plot(ext,real(PWE{2}(ext.vtx)))
title('Incident wave')
xlabel('X');   ylabel('Y');   zlabel('Z');
axis equal
view(0,10)

%%% FORMUMATIONS
disp('~~~~~~~~~~~~~ FORMULATIONS ~~~~~~~~~~~~~')

% Green kernel
Gxy    = @(X,Y) femGreenKernel(X,Y,'[exp(ikr)/r]',k);
Hxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',k) ;
Hxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',k) ;
Hxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',k) ;

% Surfacic finite elements for magnetic current
Jh = fem(ext,'RWG');
N1 = size(Jh.unk,1);

% Volumique finite elements for electric field
Eh = fem(vol,'NED');
Eh = dirichlet(Eh,int);
N2 = size(Eh.unk,1);

% Restriction matrix for regularization
[~,I1,I2] = intersect(Jh.unk,Eh.unk,'rows');
P         = sparse(I1,I2,1,size(Jh.unk,1),size(Eh.unk,1));

% Sparse operators
tic
rotErotE = integral(omega,curl(Eh),curl(Eh));
EE       = integral(omega,Eh,Eh);
JE       = integral(sigma,Jh,Eh);
toc

% Full operators
tic
JTJ = 1/(4*pi) .* integral(gamma, gamma, Jh, Gxy, Jh,tol) ...
    - 1/(4*pi*k^2) * integral(gamma, gamma, div(Jh), Gxy, div(Jh),tol) ;
JTJ = JTJ + 1/(4*pi) * regularize(gamma, gamma, Jh, '[1/r]', Jh) ...
      - 1/(4*pi*k^2) * regularize(gamma, gamma, div(Jh), '[1/r]', div(Jh));
toc

tic
JKExn = - 1/(4*pi) * integral(gamma, sigma, Jh, Hxy, nx(Eh),tol);
JKExn = JKExn - 1/(4*pi) * regularize(gamma, gamma, Jh, 'grady[1/r]', Jh) * P;
toc

% LHS = [A B ; C D]
A = 1i*k *  JTJ;
size(A)

B = JKExn - 0.5 * JE;
size(B)

C = JE.';
size(C)

D = 1/(1i*k) * (rotErotE - k^2*EE);
size(D)

% Right hand side
Y   = - integral(gamma,Jh,PWE);
RHS = [Y;zeros(size(D,1),1)];

%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% Factorization LU H-Matrix
tic
[La,Ua] = lu(A);
toc

figure
subplot(1,2,1)
spy(La)
subplot(1,2,2)
spy(Ua)

% Shurr complement resolution
tic
Sm1V = @(V) Ua\(La\V);
SV   = @(V) A*V - B*(D\(C*V));
J    = gmres(SV,Y,[],tol,100,Sm1V);
E    = - D\(C*J);
toc

%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
Ninf  = 1e3;
theta = 2*pi/1e3 .* (1:Ninf)';
nu    = [sin(theta),zeros(size(theta)),cos(theta)];

% Green kernel function
xdoty = @(X,Y) X(:,1).*Y(:,1) + X(:,2).*Y(:,2) + X(:,3).*Y(:,3);
Ginf  = @(X,Y) exp(-1i*k*xdoty(X,Y));

% Magnetic current radiation
Tinf = integral(nu,gamma,Ginf,Jh);
Jinf = 1i*k/(4*pi)*cross(nu, cross([Tinf{1}*J, Tinf{2}*J, Tinf{3}*J], nu));

% Electric field radiation
Kinf = integral(nu,sigma,Ginf,nx(Eh));
Einf = 1i*k/(4*pi) * cross(nu, [-Kinf{1}*E, -Kinf{2}*E, -Kinf{3}*E] );

% Total electric field radiation
sol = Jinf - Einf;

% Radiation infinie de reference, convention e^(+ikr)/r
nMax = 100; refInf = zeros(Ninf,1);
if E(1) == 1
    for jj = 1:Ninf
        refInf(jj,:) = sphereMaxwell(1, -f, theta(jj), 0.0, nMax);
    end
else
    for jj = 1:Ninf
        [~,refInf(jj)] = sphereMaxwell(1, -f, theta(jj), pi/2, nMax);
    end
end
refInf = refInf ./ sqrt(4*pi);

% Radiations infinies en X
if E(1) == 1
    sol = sin(theta)'.*sol(:,3) - cos(theta)'.*sol(:,1);
else
    sol = sol(:,2);
end

% Erreur
eL2   = norm(refInf-sol,2)/norm(refInf,2)
eLINF = norm(refInf-sol,'inf')/norm(refInf,'inf')

% Representation graphique
figure
subplot(1,2,1)
plot(theta,20*log10(abs(sol)),'b',theta,20*log10(abs(refInf)),'r--')

subplot(1,2,2)
plot(theta,real(sol),'--b', theta,imag(sol),'--r', theta, real(refInf),':b', theta,imag(refInf),':r');
drawnow

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan