text:
% statistic_image: An object that allows for storage and manipulation of
% images containing t-values, p-values, and thresholded results
% - multiple images can be stored in a single object
% - this is a subclass of image_vector and inherits all its methods
% - easy thresholding without changing underlying image (threshold method)
% - easy visualization (orthviews, surface, and montage methods)
%
% Usage:
% -------------------------------------------------------------------------
%  object constructor: this function creates an object and populates it
%  with data.  Entering image name(s) loads in the data.
%
%  obj = statistic_image(varargin)
%  obj = statistic_image(nifti/img filename, or allowed fields followed by values)
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% Author and copyright information:
% -------------------------------------------------------------------------
%     Copyright (C) 2010 Tor Wager
%
%     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/>.
%
% Inputs:
% -------------------------------------------------------------------------
% :properties: any valid attribute/property of statistic_image, followed by
% data values to store
%   - see properties(statistic_image)
%
% :image names: a character array with one or more filenames for
% Analyze/NIFTI images (.img/nii).
%
%
% Outputs:
% -------------------------------------------------------------------------
% obj           A statistic_image object
%
% Examples:
% -------------------------------------------------------------------------
% As with any object class in this toolbox, you can create an object by
% specifying names of fields paired with values. You can also enter
% filenames when you call statistic_image and create an image by loading a
% file. Finally, some methods performed on other objects (e.g., predict,
% regress, and ttest for fmri_data objects) will return statistic_image
% objects.
% For example, the line of code below creates a map of random t-values:
% t = statistic_image('dat', rand(1000, 1), 'type', 't', 'dfe', 29);
%
% If you want to use the full functionality of the object type without
% errors, however, a .volInfo field will need to be attached with
% information about the image space.  The easiest way to create this is to
% load an existing image, which will automatically read in the space along
% with other info:
%
% Read in a p-value image as a p-value object type:
% p = which('nonnoc_v11_4_137subjmap_weighted_pvalue.nii');
% pstat = statistic_image(p, 'type', 'p');
% pstat = threshold(pstat, .05, 'fdr');
%
% Start with the name of a statistic image we're interested in, with arbitrary values
% name = 'salientmap.nii';
%
% % Load it into a statistic_image object
% img = statistic_image('image_names', name);
%
% % Plot a histogram of its values
% figure; histogram(img)
%
% % Threshold the image values; save values < 0 or > 5
% img = threshold(img, [0 5], 'raw-outside');
%
% % Plot a montage and spm-orthviews of the thresholded image
% montage(img);
% orthviews(img);
%
% % Remove the threshold; include all values between -Inf and Inf
% img = threshold(img, [-Inf Inf], 'raw-between');
%
% % Show the orthviews plot of the unthresholded image
% orthviews(img);
%
% You can force an image to be of a particular type, e.g., 'T', in which
% case p-values will automatically be added. This is useful for
% thresholding.
% e.g., for an SPM T-map:
% wh = 1;
% load SPM
% img = sprintf('spmT_00%02d.img', wh(1))
% t = statistic_image('image_names', img, 'type', 'T', 'dfe', SPM.xX.erdf);
%
% If 'type' is 'robreg', then assumes user is in a robust regression
% directory and will load in the robust_beta_000X.img and load in
% robust_p_000X.img as the p values.  X is assumed = 1, but can be
% overloaded by passing in a dat_descrip field. i.e.
% deltadon=statistic_image('dat_descrip', 4, 'type', 'robreg')
%
%
% See also:
% image_vector.m
% fmri_data.m
% region.m
% fmridisplay.m
%
classdef statistic_image < image_vector

    properties

        type        % String with image type: 'generic', 't', 'p', 'robreg'
        p           % Matrix of p-values for images, [voxels x images]
        p_type      % String with source info for p-values
        ste         % Matrix of standard error values for images, [voxels x images]
        threshold   % Latest statistical threshold applied
        thr_type    % Information about statistical threshold applied
        sig         % Logical matrix of which voxels are significant, [voxels x images]
        N           % Sample size
        dfe         % Error degrees of freedom for test
        image_labels ={''}; % Semantic labels for each image, cell array
    end

    methods

        function obj = statistic_image(varargin)

            % Enter fieldname', value pairs in any order to create class
            % instance

            % ---------------------------------
            % Create empty image_vector object, and return if no additional
            % arguments
            % ---------------------------------
            obj.type = 'generic';
            %obj.X = []; alread defined

            obj.p = [];
            obj.p_type = [];
            obj.ste = [];
            obj.threshold = [];
            obj.thr_type = [];
            obj.sig = [];
            obj.dfe = [];
            obj.N = [];
            obj.image_labels = {};
            doverbose = false;

            wh = strcmp(varargin, 'verbose');
            if any(wh), doverbose = true; varargin(wh) = []; end

            % The code below can be generic to any class definition
            % It parses 'fieldname', value pairs of inputs
            % and returns a warning if unexpected strings are found.

            if nargin == 0
                return
            end

            % all valid fieldnames
            valid_names = fieldnames(obj);

            for i = 1:length(varargin)
                if ischar(varargin{i})

                    % Look for a field (attribute) with the input name
                    wh = strmatch(varargin{i}, valid_names, 'exact');

                    % behaviors for valid fields
                    if ~isempty(wh)

                        obj.(varargin{i}) = varargin{i + 1};

                        % eliminate strings to prevent warnings on char
                        % inputs
                        if ischar(varargin{i + 1})
                            varargin{i + 1} = [];
                        end

                        % special methods for specific fields
                        switch varargin{i}

                        end

                    else
                        % Assume is char array of image to load
                        % could also use obj = read_from_file(obj);

                        if doverbose
                            fprintf('Reading from file: %s\n', varargin{i});
                            tmp = fmri_data(varargin{i});
                        else
                            tmp = fmri_data(varargin{i}, 'noverbose');
                        end

                        obj.dat = tmp.dat;
                        obj.volInfo = tmp.volInfo;
                        obj.image_names = tmp.image_names;
                        obj.fullpath = tmp.fullpath;
                        obj.files_exist = tmp.files_exist;
                        obj.history = tmp.history;
                        %warning('inputargs:BadInput', sprintf('Unknown field: %s', varargin{i}));

                    end
                end % string input
            end % process inputs

            % load data, if we don't have data yet
            % But we do have valid image names.
            % ----------------------------------------

            % if loading in a robust regression image, overload image_names
            %    Yoni originally added this part (2014-07-24), and Wani
            %    added if-statement to fix a bug related to this part (2014-08-26).
            if strcmp(obj.type, 'robreg')
                if isempty(obj.dat_descrip), obj.dat_descrip=1; end %by default, assume user wants rob img 1 if none entered
                obj.image_names = sprintf('rob_tmap_%04d.img',obj.dat_descrip);
            end

            if doverbose
            obj = check_image_filenames(obj);
            else
                 obj = check_image_filenames(obj, 'noverbose');
            end

            if isempty(obj.dat) && any(obj.files_exist)
                obj = read_from_file(obj);

            elseif isempty(obj.dat)
                disp('Warning: .dat is empty and files cannot be found.  No image data in object.');
            end

            % 6/22/13 Tor Added to enforce consistency in objects across usage cases
            if isempty(obj.sig) || (numel(obj.sig) == 1 && ~obj.sig)
                obj.sig = true(size(obj.dat));
            end

            % type - handle specific forced types

            switch obj.type
                case {'T', 't'}
                    obj.type = 'T';

                    if isempty(obj.dfe)
                        error('If forcing object type = ''T'', enter dfe in call to object constructor.')
                    end

                    obj.p = 2 * (1 - tcdf(abs(obj.dat), obj.dfe));
                    obj.p_type = '2-tailed P-value from input dfe';

                case 'p'
                    obj.type = 'p';
                    obj.p = obj.dat;

                case 'r'
                    obj.type = 'r';

                case 'robreg'
                    obj.type = 'robust regression:  .dat is beta values, .p is robust p values';

                    pimg = image_vector('image_names', sprintf('rob_p_%04d.img',obj.dat_descrip));
                    % bimg = image_vector('image_names', sprintf('rob_beta_%04d.img',obj.dat_descrip));
                    timg = image_vector('image_names', sprintf('rob_tmap_%04d.img',obj.dat_descrip));

                    obj.p = pimg.dat;
                    obj.p_type = 'robust';
                    obj.dat = timg.dat;%bimg.dat;
                    obj.dat_descrip = ['robust regression for covariate number ' num2str(obj.dat_descrip)];

                case 'generic'
                    % do nothing

                otherwise
                    error('User forced unknown statistic image type.')
            end

        end % class constructor

    end  % methods

end  % classdef

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%%  PosControl_Sim
clear
path('./icon/',path);
Init;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
%Attitude PID parameters
Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
%Position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
PosControl_Sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
% OP_GRADU_GRADV: assemble the stiffness matrix A = [a(i,j)], a(i,j) = (epsilon grad u_j, grad v_i).
%
%   mat = op_gradu_gradv (spu, spv, msh, epsilon);
%   [rows, cols, values] = op_gradu_gradv (spu, spv, msh, epsilon);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_scalar/sp_evaluate_col)
%   spv:   structure representing the space of test functions (see sp_scalar/sp_evaluate_col)
%   msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%   epsilon: diffusion coefficient
%
% OUTPUT:
%
%   mat:    assembled stiffness matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    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/>.

function varargout = op_gradu_gradv (spu, spv, msh, coeff)

  gradu = reshape (spu.shape_function_gradients, spu.ncomp, [], ...
		   msh.nqn, spu.nsh_max, msh.nel);
  gradv = reshape (spv.shape_function_gradients, spv.ncomp, [], ...
		   msh.nqn, spv.nsh_max, msh.nel);

  ndir = size (gradu, 2);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      gradu_iel = reshape (gradu(:,:,:,:,iel), spu.ncomp*ndir, msh.nqn, 1, spu.nsh_max);
      gradv_iel = reshape (gradv(:,:,:,:,iel), spv.ncomp*ndir, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      jacdet_gradu = bsxfun (@times, jacdet_iel, gradu_iel);
      tmp1 = sum (bsxfun (@times, jacdet_gradu, gradv_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_gradu_gradv: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_gradu_gradv: wrong number of output arguments')
  end

end

%% COPY OF THE FIRST VERSION OF THE FUNCTION (MORE UNDERSTANDABLE)
%
% function mat = op_gradu_gradv (spu, spv, msh, coeff)
%
%   mat = spalloc (spv.ndof, spu.ndof, 1);
%
%   gradu = reshape (spu.shape_function_gradients, spu.ncomp, [], msh.nqn, spu.nsh_max, msh.nel);
%   gradv = reshape (spv.shape_function_gradients, spv.ncomp, [], msh.nqn, spv.nsh_max, msh.nel);
%
%   ndir = size (gradu, 2);
%
%   for iel = 1:msh.nel
%     if (all (msh.jacdet(:,iel)))
%       mat_loc = zeros (spv.nsh(iel), spu.nsh(iel));
%       for idof = 1:spv.nsh(iel)
%         ishg = reshape(gradv(:,:,:,idof,iel),spv.ncomp * ndir, []);
%         for jdof = 1:spu.nsh(iel)
%           jshg = reshape(gradu(:,:,:,jdof,iel),spu.ncomp * ndir, []);
% % The cycle on the quadrature points is vectorized
%           %for inode = 1:msh.nqn
%           mat_loc(idof, jdof) = mat_loc(idof, jdof) + ...
%              sum (msh.jacdet(:,iel) .* msh.quad_weights(:, iel) .* ...
%                   sum (ishg .* jshg, 1).' .* coeff(:,iel));
%           %end
%         end
%       end
%       mat(spv.connectivity(:, iel), spu.connectivity(:, iel)) = ...
%         mat(spv.connectivity(:, iel), spu.connectivity(:, iel)) + mat_loc;
%     else
%       warning ('geopdes:jacdet_zero_at_quad_node', 'op_gradu_gradv: singular map in element number %d', iel)
%     end
%   end
%
% end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [h,g,a,info] = wfilt_hden(K)
%WFILT_HDEN  Higher DENsity dwt filters (tight frame, frame)
%   Usage: [h,g,a] = wfilt_hden(K);
%
%   `[h,g,a]=wfilt_hden(K)` with $K \in {1,2,3,4}$ returns Higher DENsity
%   dwt filters (tight frame, frame) from the reference. The filterbanks
%   have 3 channels and unusual non-uniform subsamplig factors `[2,2,1]`.
%
%   Examples:
%   ---------
%   :::
%     wfiltinfo('hden3');
%
%   :::
%     wfiltinfo('ana:hden4');
%
%   References: selesnick2006higher
%

% AUTHOR: Zdenek Prusa

a= [2;2;1];
info.istight = 1;
switch(K)
case 1
% from the paper Example 1.
garr = [
    0                   0                    0
    0.353553390593274   0.353553390593274    0.5
    0.707106781186548   0                    -0.5
    0.353553390593274   -0.353553390593274   0
];
d = [-2,-2,-2];
case 2
% from the paper Example 2.
garr = [
    0                 0                0
    0.189604909379    0.025752563665   0.010167956157
    0.631450512121    0.075463998066   0.046750380120
    0.655505518357   -0.064333341412  -0.009172584871
    0.099615139800   -0.327704691428  -0.354664087684
    -0.163756210215   0.228185687127   0.499004628714
    -0.023958870736   0.252240693362  -0.192086292435
    0.025752563665   -0.189604909379   0
];

d = [-3,-5,-5];
case 3
% from the paper Example 3.
garr = [
   0                 0                0
   0.022033327573    0.048477254777   0.031294135831
   0.015381522616    0.019991451948   0.013248398005
  -0.088169084245   -0.304530024033  -0.311552292833
   0.051120949834    0.165478923930   0.497594326648
   0.574161374258    0.308884916012  -0.235117092484
   0.717567366340   -0.214155508410  -0.020594576659
   0.247558418377   -0.074865474330   0.015375249485
  -0.076963057605    0.028685132531   0.009751852004
  -0.048477254777    0.022033327573   0
];
 d = [-6,-4,-4];
case 4
    info.istight = 0;
    % from the paper Example 5. Is not a tight frame!
    harr = [
       0          0           0
       0          0           0
       0.027222   0.044889    0
       0.011217   0.005671    0.027671
      -0.112709  -0.286349    0.007159
       0.096078   0.235789   -0.277671
       0.685299   0.235789    0.485682
       0.685299   -0.286349  -0.277671
       0.096078   0.005671    0.007159
      -0.112709   0.044889    0.027671
       0.011217   0           0
       0.027222   0           0
    ];
    d = [-5,-5,-5];

    harr = flipud(harr);
    h=mat2cell(harr,size(harr,1),ones(1,size(harr,2)));
    h=cellfun(@(gEl,dEl) struct('h',gEl,'offset',dEl),h,num2cell(d),...
              'UniformOutput',0);

        garr = [
            0          0          0
           -0.039237   0.023794   0.011029
           -0.073518   0.037784   0.019204
            0.181733  -0.070538  -0.020024
            0.638129  -0.253037  -0.269204
            0.638129   0.261996   0.517991
            0.181733   0.261996  -0.269204
           -0.073518  -0.253037  -0.020024
           -0.039237  -0.070538   0.019204
            0          0.037784   0.011029
            0          0.023794   0
            0          0          0
        ];

        g=mat2cell(garr,size(garr,1),ones(1,size(garr,2)));
        g=cellfun(@(gEl,dEl) struct('h',gEl,'offset',dEl),g,num2cell(d),...
                  'UniformOutput',0);

    return;

otherwise
        error('%s: No such Higher Density Wavelet Transform Filters..',upper(mfilename));
end

%garr = flipud(harr);
g=mat2cell(garr,size(garr,1),ones(1,size(garr,2)));
g = cellfun(@(gEl,dEl) struct('h',gEl,'offset',dEl),g,num2cell(d),'UniformOutput',0);
h = g;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
close all;
clear variables
pde = HodgeLaplacianEdata1;
% pde = HodgeLaplacianFdata1;
[node,elem] = squaremesh([0,1,0,1],1/8);
% bdFlag = setboundary(node,elem,'Neumann');
% Pure Neumann boundary condition doesn't work.
% bdFlag = setboundary(node,elem,'Dirichlet');
bdFlag = setboundary(node,elem,'Dirichlet','x==0','Neumann','~(x==0)');
err = zeros(4,1); N = zeros(4,1);
option.solver = 'direct';
for i = 1:4
    [sigma,u] = HodgeLaplacianE(node,elem,bdFlag,pde,option);
%     [sigma,u,AD] = HodgeLaplacianF(node,elem,pde,bdFlag,option);
    err(i) = getL2error(node,elem,pde.sigma,sigma);
    N(i) = size(u,1);
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
end
showrate(N,err,2);

repo: ifem


text:
% housekeeping
close all
clear
clc
%% load estimated models

load tut13_bayestimation

load tut01_data

endog=models.ve.endogenous;

%% posterior sampling of parameters
params=struct();
params.ve=models.ve.estim_.sampler(1000);
params.ve_lr=models.ve_lr.estim_.sampler(1000);

%% forecast over all sampled parameters
clc
myfkst=struct();
% date_start=[];
date_start='2003Q1';
myfkst.ve=forecast(models.ve,db,date_start,params.ve);

myfkst.ve_lr=forecast(models.ve_lr,db,date_start,params.ve_lr);

%% set environment

ci=[30,50,68,90];

%% Forecast plots

modelnames=fieldnames(models);

for jj=1:numel(modelnames)

    modname=modelnames{jj};

    figure('name',['model (',modname,') Forecasts of Norwegian Data']);

    for ii=1:numel(endog)

        subplot(3,2,ii)

        d=myfkst.(modname).(endog{ii});

        out=fanchart(d,ci);

        plot_fanchart(out,[244, 122, 66]/244)%'c'

        %     hold on
        %
        %     plot('1995:2004',db.(endog{ii}))

        title(tex.(endog{ii}))

        axis tight

    end

    xrotate(45)

end

author: jmaih


text:
% Computes the step response of a linear system

function X = simple_step(A,B,DT,N)
n  = size(A,1);
Ad = expm( full( A * DT ) );
Bd = ( Ad - eye(n) ) * B;
Bd = A \ Bd;
X  = zeros(n,N);
for k = 2 : N,
    X(:,k) = Ad*X(:,k-1)+Bd;
end

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [fx,fy] = forward_gradient(f,fx,fy, nx, ny)

for  i = 1: ny-2

          for j = 0: nx-2
                p  = i * nx + j;
                p1 = p + 1;
                p2 = p + nx;
if p <= length(f)
                fx(round(p)) = f(round(p1)) - f(round(p));
                fy(round(p)) = f(round(p2)) - f(round(p));
end
          end
end

    for ( j = 1: nx-2)
         p = round((ny-1) * nx + j);
if(p < length(f))
        fx(p) = f(p+1) - f(p);
        fy(p) = 0;
end

    for  i = 1: ny-1
      p = i * nx-1;
if((p <= length(f)) && (p > 0))
    if(p <= length(fx) && p <= length(fy))

        fx(round(p)) = 0;
        if( p+nx <= length(f))
        fy(round(p)) = (f(round(p+nx)) - f(round(p)));
        end
    end
    end
    end
    fx(round(ny * nx - 1)) = 0;
    fy(round(ny * nx - 1)) = 0;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function blend_test12 ( )

%*****************************************************************************80
%
%% BLEND_TEST12 tests BLEND_113.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
  n1 = 3;
  n2 = 5;
  n3 = 4;

  d(1:n1,1:n2,1:n3) = 0.0;

  for i = 1 : n1
    for j = 1 : n2
      for k = 1 : n3
        if ( ...
          ( ( i == 1 | i == n1 ) & ( j == 1 | j == n2 ) ) | ...
          ( ( i == 1 | i == n1 ) & ( k == 1 | k == n3 ) ) | ...
          ( ( j == 1 | j == n2 ) & ( k == 1 | k == n3 ) ) )
          d(i,j,k) = ( i + j + k );
        end
      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLEND_TEST12\n' );
  fprintf ( 1, '  BLEND_113 blends edge values into a table.\n' );
  fprintf ( 1, '\n' );

  r8block_print ( n1, n2, n3, d, '  Initial data array' );

  for i = 1 : n1

    r = ( i - 1 ) / ( n1 - 1 );

    for j = 1 : n2

      s = ( j - 1 ) / ( n2 - 1 );

      for k = 1 : n3

        t = ( k - 1 ) / ( n3 - 1 );

        if ( ...
          ( ( i == 1 | i == n1 ) & ( j == 1 | j == n2 ) ) | ...
          ( ( i == 1 | i == n1 ) & ( k == 1 | k == n3 ) ) | ...
          ( ( j == 1 | j == n2 ) & ( k == 1 | k == n3 ) ) )
          continue
        end

        d(i,j,k) = blend_113 ( r, s, t, ...
          d(1,1,1), d(1,1,n3), d(1,n2,1), d(1,n2,n3), ...
          d(n1,1,1), d(n1,1,n3), d(n1,n2,1), d(n1,n2,n3), ...
          d(i,1,1), d(i,1,n3), d(i,n2,1), d(i,n2,n3), ...
          d(1,j,1), d(1,j,n3), d(n1,j,1), d(n1,j,n3), ...
          d(1,1,k), d(1,n2,k), d(n1,1,k), d(n1,n2,k) );

      end

    end

  end

  r8block_print ( n1, n2, n3, d, '  Interpolated data array' );

  return
end

author: johannesgerer

repo: jburkardt-m


%--------------------------------------------------------------------------
% Programed by:
%   - Lucas Rath (lucasrm25@gmail.com)
%   -
%   -
%--------------------------------------------------------------------------

classdef MotionModelGP_SingleTrack_nominal < MotionModelGP
%--------------------------------------------------------------------------
%   xk+1 = fd(xk,uk) + Bd * ( d(zk) + w ),
%
%       where: zk = [Bz_x*xk ; Bz_u*uk],
%              d ~ N(mean_d(zk),var_d(zk))
%              w ~ N(0,sigmaw)
%
%
%   x = [I_x              (x position in global coordinates),
%        I_y              (y position in global coordinates),
%        psi              (yaw angle),
%        V_vx             (longitudinal velocity in vehicle coordinates)
%        V_vy             (lateral velocity in vehicle coordinates)
%        psi_dot          (yaw rate),
%        track_dist       (distance traveled in the track centerline)
%        ]
%
%   u = [delta          (steering angle),
%        T              (wheel torque gain),  -1=max.braking, 1=max acc.
%        track_vel      (velocity in the track centerline)
%       ]
%
%--------------------------------------------------------------------------

    properties(Constant)
        M    = 500 *1.5      % vehicle mass
        I_z  = 600 *1.5      % vehicle moment of inertia (yaw axis)
        g    = 9.81     % gravitation
        l_f  = 0.9      % distance of the front wheel to the center of mass
        l_r  = 1.5      % distance of the rear wheel to the center of mass

        deltamax    = deg2rad(30)   % maximum steering amplitude
        % deltadotmax = deg2rad(20)   % maximum steering velocity amplitude

        maxbrakeWForce = 8000 % = 2*g*M;  % allow ~ 2g brake
        maxmotorWForce = 4000 % = 1*g*M;  % allow ~ 1g acc

        % Pacejka lateral dynamics parameters
        c_f = 14000 * 2.5 % = 1*g*M/deltamax  % front coornering stiffness (C*delta=Fy~M*a)
        c_r = 14000 * 2.5 % = 2*g*M/deltamax  % rear coornering stiffness
    end

    properties(Constant)
        % keep in mind the dimensions:  xk+1 = fd(xk,uk) + Bd*(d(z)+w)),
        % where z = [Bz_x*x;Bz_u*u]
        Bz_x = [zeros(3), eye(3), zeros(3,1)]
        Bz_u = [1 0 0;
                0 1 0]
        Bd = [zeros(3);
              eye(3);
              zeros(1,3)]
        n  = 7   % number of outputs x(t)
        m  = 3   % number of inputs u(t)
        nz = 5   % dimension of z(t)
        nd = 3   % output dimension of d(z)
    end

    methods
        function obj = MotionModelGP_SingleTrack_nominal(d,sigmaw)
        %------------------------------------------------------------------
        %   object constructor. Create model and report model stability
        %   analysis
        %------------------------------------------------------------------
            % call superclass constructor
            obj = obj@MotionModelGP(d,sigmaw);
            % report single track dynamics analysis
            % obj.analyseSingleTrack();

            % add folder CODEGEN to path. Here there will be some functions
            % generated with the method generate_grad_functions()
            addpath(fullfile(pwd,'CODEGEN'))
        end

        function xdot = f (obj, x, u)
        %------------------------------------------------------------------
        %   Continuous time dynamics of the single track (including
        %   disturbance):
        %------------------------------------------------------------------

            %--------------------------------------------------------------
            % Model parameters
            %--------------------------------------------------------------
            g = obj.g;
            % switch to symbolic variables if input is symbolic
            % (to be used with automatic gradient code generation)
            if isa(x,'sym')
                syms M I_z l_f l_r deltamax maxbrakeWForce maxmotorWForce c_f c_r real
            else
                M = obj.M;
                I_z  = obj.I_z;
                l_f  = obj.l_f;
                l_r  = obj.l_r;
                deltamax = obj.deltamax;
                maxbrakeWForce = obj.maxbrakeWForce;
                maxmotorWForce = obj.maxmotorWForce;
                c_f = obj.c_f;
                c_r = obj.c_r;
            end

            %--------------------------------------------------------------
            % State Vector
            %--------------------------------------------------------------
            I_x         = x(1);
            I_y         = x(2);
            psi         = x(3);
            V_vx        = x(4);
            V_vy        = x(5);
            psi_dot     = x(6);
            track_dist  = x(7);
            beta = atan2(V_vy,V_vx);

            %--------------------------------------------------------------
            % Inputs
            %--------------------------------------------------------------
            delta     = u(1);   % steering angle
            T         = u(2);   % wheel torque gain,  -1=max.braking, 1=max acc.
            track_vel = u(3);   % track centerline velocity

            %--------------------------------------------------------------
            % Saturate inputs and
            %--------------------------------------------------------------
            % saturate steering angle
            if isa(x,'sym') % switch to differentiable function
                delta = obj.sclip(delta, -deltamax, deltamax);
            else
                delta = obj.clip(delta, -deltamax, deltamax); % (NOT DIFFERENTIABLE)
            end

            % saturate pedal input
            if isa(x,'sym') % switch to differentiable function
                T = obj.sclip( T, -1, 1);
            else
                T = obj.clip( T, -1, 1);    % % (NOT DIFFERENTIABLE)
            end

            %--------------------------------------------------------------
            % Wheel slip angles (slip ration not being used for now)
            %--------------------------------------------------------------
            a_r = atan2(V_vy-l_r*psi_dot,V_vx);
            a_f = atan2(V_vy+l_f*psi_dot,V_vx) - delta;

            %--------------------------------------------------------------
            % Tyre forces
            %--------------------------------------------------------------
            % desired total wheel torque to be applied
            if isa(x,'sym') % switch to differentiable function
                totalWForce = T*(  (obj.gez(T)).*maxmotorWForce ...
                              +(obj.lez(T)).*maxbrakeWForce.*obj.ssign(V_vx));
            else
                totalWForce = T * ( (T>0)*maxmotorWForce+(T<0)*maxbrakeWForce*sign(V_vx) ); % % (NOT DIFFERENTIABLE)
            end
            % longitudinal wheel torque distribution
            zeta = 0.5;

            % Wheel forces in wheel coordinates (z-axis points down, x-axis front)
            % This means positive W_Fy_f turns vehicle to the right
            W_Fx_r = zeta * totalWForce;
            W_Fx_f = (1-zeta) * totalWForce;

            W_Fy_r = c_r * a_r;
            W_Fy_f = c_f * a_f;
            % W_Fy_r = D_r*sin(C_r*atan(B_r*a_r-E_r*(B_r*a_r -atan(B_r*a_r)))); % rear lateral force
            % W_Fy_f = D_f*sin(C_f*atan(B_f*a_f-E_f*(B_f*a_f -atan(B_f*a_f)))); % front lateral force

            % Wheel forces in vehicle coordinates (z-axis points up, x-axis front)
            V_Fx_r = W_Fx_r;
            V_Fx_f = W_Fx_f;
            V_Fy_r = - W_Fy_r;
            V_Fy_f = - W_Fy_f;

            %--------------------------------------------------------------
            % Calculate state space time derivatives
            %--------------------------------------------------------------
            % vector field (right-hand side of differential equation)
            I_x_dot = V_vx*cos(psi) - V_vy*sin(psi); % longitudinal velocity
            I_y_dot = V_vx*sin(psi) + V_vy*cos(psi); % lateral velocity
            V_vx_dot = 1/M * (V_Fx_r + V_Fx_f*cos(delta) - V_Fy_f*sin(delta) + V_vy*psi_dot);
            V_vy_dot = 1/M * (V_Fy_r + V_Fx_f*sin(delta) + V_Fy_f*cos(delta) - V_vy*psi_dot);
            psi_dot_dot = 1/I_z * (V_Fy_f*l_f*cos(delta) + V_Fx_f*l_f*sin(delta) - V_Fy_r*l_r);
            track_dist_dot = track_vel; % Traveled distance in the track centerline

            %--------------------------------------------------------------
            % write outputs
            %--------------------------------------------------------------
            xdot  = [I_x_dot; I_y_dot; psi_dot; V_vx_dot; V_vy_dot; psi_dot_dot; track_dist_dot];

            %--------------------------------------------------------------
            % check model integrity
            %--------------------------------------------------------------
            if ~isa(x,'sym')    % sym input used for code generation
                if any(isnan(xdot)) || any(isinf(xdot)) || any(imag(xdot)~=0)
                    error('Single Track Model evaluated to Inf of NaN... CHECK MODEL!!!')
                end
            end
        end

        function gradx = gradx_f(obj, x, u)
        %------------------------------------------------------------------
        %   Continuous time dynamics.
        %   out:
        %       gradx: <n,n> gradient of xdot w.r.t. x
        %------------------------------------------------------------------
            % gradx = zeros(obj.n);
            params = [obj.M obj.I_z obj.l_f obj.l_r obj.deltamax ...
                      obj.maxbrakeWForce obj.maxmotorWForce obj.c_f obj.c_r]';
            gradx = singletrack_gradx_f(x,u,params);
        end

        function gradu = gradu_f(obj, x, u)
        %------------------------------------------------------------------
        %   Continuous time dynamics.
        %   out:
        %       gradu: <m,n> gradient of xdot w.r.t. u
        %------------------------------------------------------------------
            % gradu = zeros(obj.m,obj.n);
            params = [obj.M obj.I_z obj.l_f obj.l_r obj.deltamax ...
                      obj.maxbrakeWForce obj.maxmotorWForce obj.c_f obj.c_r]';
            gradu = singletrack_gradu_f(x,u,params);
        end

        function generate_grad_functions(obj)
        %------------------------------------------------------------------
        %   Generate external files for the evaluation of the gradient of
        %   the continuous time dynamics. (Make use of symbolic toolbox)
        %   Please ensure that your dynamics only contains smooth
        %   diferentiable functions.
        %
        %   To generate files, simply call:
        %       nomModel = MotionModelGP_SingleTrack_nominal(@(z)deal(zeros(3,1),zeros(3)), zeros(3));
        %       nomModel.generate_grad_functions()
        %------------------------------------------------------------------
            syms I_x I_y vpsi V_vx V_vy psi_dot track_dist real
            x = [I_x I_y vpsi V_vx V_vy psi_dot track_dist]';
            syms delta T track_vel real
            u = [delta T track_vel]';
            syms M I_z l_f l_r deltamax maxbrakeWForce maxmotorWForce c_f c_r real
            params = [M I_z l_f l_r deltamax maxbrakeWForce maxmotorWForce c_f c_r]';

            xdot = obj.f(x,u);
            gradx = jacobian(xdot,x)';
            gradu = jacobian(xdot,u)';
            % gradx = simplify(expand(gradx));  % does not work. eqs are too complex
            % gradu = simplify(expand(gradu));  % does not work. eqs are too complex

            % Save model
            folder = fullfile(pwd,'CODEGEN');
            if ~exist(folder,'dir')
                mkdir(folder);
            end
            addpath(folder)

            matlabFunction(gradx,'Vars',{x,u,params},'File',fullfile('CODEGEN','singletrack_gradx_f'),'Optimize',true);
            matlabFunction(gradu,'Vars',{x,u,params},'File',fullfile('CODEGEN','singletrack_gradu_f'),'Optimize',true);

            disp('FINISHED! functions CODEGEN_singletrack_gradx_f and CODEGEN_singletrack_gradu_f generated!!')
        end

        function analyseSingleTrack(obj)
        %------------------------------------------------------------------
        %   report single track model dynamics analysis
        %------------------------------------------------------------------
            fprintf('Single track model created!!! Summary:\n');
            % vehicle size
            l = obj.l_f + obj.l_r;
            % Eigenlenkgradient = Yaw gradient (d_delta/d_ay)
            EG = obj.M/l*(obj.l_r/obj.c_f - obj.l_f/obj.c_r);
            fprintf(2,'\tYaw gradient EG=(d_delta/d_ay) ~ %.1f [deg/g]\n',rad2deg(EG)*obj.g);
            if rad2deg(EG)*obj.g < 10
                fprintf(2,'\tBE CAREFULL... EG is too low (< 10 [deg/g])... increase l_r,c_r or decrease l_f,c_f\n');
            end
            if EG < 0
                v_cr = sqrt( obj.c_f*obj.c_r*l^2 / (obj.M*(obj.c_f*obj.l_f - obj.c_r*obj.l_r)));
                fprintf(2,'\tVehicle is intrinsically OVERSTEER (c_r*l_r < c_f*l_f)\n');
                fprintf('\tCritical velocity v_cr ~ %.1f [m/s]\n', v_cr);
                fprintf('\tMax ss. Yaw (at delta_max) at v= 5[m/s] ~ %.1f [deg/s]\n',rad2deg(obj.deltamax)* 5/(l+EG* 5^2));
                fprintf('\tMax ss. Yaw (at delta_max) at v=10[m/s] ~ %.1f [deg/s]\n',rad2deg(obj.deltamax)*10/(l+EG*10^2));
                fprintf('\tMax ss. Yaw (at delta_max) at v=15[m/s] ~ %.1f [deg/s]\n',rad2deg(obj.deltamax)*15/(l+EG*15^2));
            else
                v_ch = sqrt(l/EG);
                fprintf(2,'\tVehicle is intrinsically UNDERSTEER (c_r*l_r > c_f*l_f)\n');
                fprintf('\tCharacteristic velocity v_ch ~ %.1f [m/s]\n',v_ch);
                fprintf('\tMax ss. Yaw rate gain (at v_ch)          (psidot/delta)ss,max ~ %.1f [1/s]\n',v_ch/(l+EG*v_ch^2));
                fprintf('\tMax ss. Yaw rate (at v_ch and delta_max) (psidot)ss,max~ %.1f [deg/s]\n',rad2deg(obj.deltamax)*v_ch/(l+EG*v_ch^2));
            end
            fprintf('\tax_{max} ~ %.1f [g]\n',obj.maxmotorWForce/obj.M/obj.g);
            fprintf('\tax_{min} ~ -%.1f [g]\n',obj.maxbrakeWForce/obj.M/obj.g);
        end
    end

    %----------------------------------------------------------------------
    % Smooth alternatives for nonmooth functions
    %----------------------------------------------------------------------
    methods
        function x = sclip(obj,x,lb,ub)
            % Smooth (differentiable) clip (saturation) function
            x = x.*obj.gez(x-lb).*obj.lez(x-ub) + ub*obj.gez(x-ub) + lb*obj.lez(x-lb);
        end
    end
    methods(Static)
        function x = clip(x,lb,ub)
            % standard nonsmooth clip (saturation) function
            x = min(max(x,lb),ub);
        end
        function x = srec(x,lb,ub)
            % Smooth rectangular function
            alpha = 50; % the larger the sharper the rectangular function
            x = 0.5*(tanh((x-lb)*alpha)-tanh((x-ub)*alpha));
        end
        function x = gez(x)
            % Smooth >=0 boolean function
            alpha = 50; % the larger the sharper the clip function
            x = (1+exp(-alpha*x)).^-1;
        end
        function x = lez(x)
            % Smooth <=0 boolean function
            alpha = 50; % the larger the sharper the clip function
            x = 1-(1+exp(-alpha*x)).^-1;
        end
        function x = ssign(x)
            % Smooth sign(x) boolean function
            alpha = 100; % the larger the sharper the clip function
            x = tanh(alpha*x);
        end
    end
end


text:
function view = makeMultiLevelButtons(view);
%
% view = makeHemisphereButtons(view);
%
% For flat multi-level view:
% Installs radio buttons that set whether to
% display flat data from tSeries averaged across
% gray levels, or whether to display each gray level
% separately
%
% djh, 1/97 (makeHemisphereButtons)
% ras, 8/04 (makeMultiLevelButtons)

% Multi button Callback:
%   selectButton(view.ui.sliceButtons,1);
%   view=refreshScreen(view);
% The hemisphere buttons (off of which these buttons are hacked)
% pass a '2' in the recompute field of refreshView.
% This lets refreshView deal with the image rotation slider bar.
% I'm not sure if this is relevant to the gray levels, but I'll
% leave it in for now. -ras 08/04
callbackStr = ...
    ['selectButton(',view.name,'.ui.levelButtons,1); ',...
     'switch2AvgLevels(',view.name,'); ',...
	view.name,'=refreshScreen(',view.name,',2);'];

view.ui.levelButtons(1) = ...
    uicontrol('Style','radiobutton',...
    'String','Average',...
    'Units','normalized',...
    'BackgroundColor',get(gcf,'Color'),...
    'Position',[0 .65 .1 .05],...
    'Callback', callbackStr);

% Separate button Callback:
%   selectButton(view.ui.sliceButtons,1);
%   view=refreshScreen(view);
callbackStr = ...
    ['selectButton(',view.name,'.ui.levelButtons,2); ',...
     'switch2SeparateLevels(',view.name,'); ',...
	view.name,'=refreshScreen(',view.name,',2);'];
view.ui.levelButtons(2) = ...
    uicontrol('Style','radiobutton',...
    'String','Separate',...
    'BackgroundColor',get(gcf,'Color'),...
    'Units','normalized',...
    'Position',[0 .6 .1 .05],...
    'Callback', callbackStr);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ar = psd2ar(h,p)

%PSD2AR Power Spectral Density Function to AR parameters.
%
%  psd2ar(h,p) calculates an AR(p) model from equidistant samples of a power
%  spectrum h.
%
%  See also: ARMA2PSD, MODERR, SIMUARMA.

%S. de Waele, November 2001.

hs = [h fliplr(h(2:end))];
corip = ipifft(hs);
corip = corip(1:1+p)/corip(1);
ar = cov2arset(corip);

%------------------------------------------------------------------
function rho = ipifft(h)

%IPIFFT InterPolated Inverse Fast Fourier Transform
%
%  ipfft(h) considers the vector h as samples of a spectrum
%  between 0 and 2*PI. The spectrum is interpolated using
%  Nearest Neighbor interpolation.
%  If the samples in h are all positive, the resulting
%  continuous spectrum is also positive.
%
%  See also: IFFT, FFT.

%S. de Waele, October 2001.

n = length(h);
d_om = 2*pi/n; %Omega interval
r = 1:n-1;
cf = 1/d_om*2*imag(exp(i*d_om/2.*r))./r;
cf = [1 cf];
rho = real(cf.*ifft(h));

author: Sable

repo: mcbench-benchmarks


text:
function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport] = getINITModel2(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, removeGenes, taskFile, useScoresForTasks, printReport, taskStructure, params, paramsFT)
% getINITModel2
%   Note: There is a newer improved version of this function in RAVEN Toolbox
%   called ftINIT that can be used together with the helper functions in this repo.
%
%   Generates a model using the INIT algorithm, based on proteomics and/or
%   transcriptomics and/or metabolomics and/or metabolic tasks.
%   This is the newer version, which has updated handling of gene rules to
%   differentiate between isozymes and enzyme complexes.
%
%   refModel            a model structure. The model should be in the
%                       closed form (no exchange reactions open). Import
%                       using import(filename,false). If the model is not
%                       loaded using importModel, it might be that there
%                       is no "unconstrained" field. In that case,
%                       manually add the field like:
%                       model.unconstrained = false(numel(model.mets),1);
%   tissue              tissue to score for. Should exist in either
%                       hpaData.tissues or arrayData.tissues
%   celltype            cell type to score for. Should exist in either
%                       hpaData.celltypes or arrayData.celltypes for this
%                       tissue (opt, default is to use the max values
%                       among all the cell types for the tissue.
%   hpaData             HPA data structure from parseHPA (opt if arrayData
%                       is supplied, default [])
%   arrayData           gene expression data structure (opt if hpaData is
%                       supplied, default [])
%       genes           cell array with the unique gene names
%       tissues         cell array with the tissue names. The list may not
%                       be unique, as there can be multiple cell types per
%                       tissue
%       celltypes       cell array with the cell type names for each tissue
%       levels          GENESxTISSUES array with the expression level for
%                       each gene in each tissue/celltype. NaN should be
%                       used when no measurement was performed
%       threshold       a single value or a vector of gene expression
%                       thresholds, above which genes are considered to be
%                       "expressed". (opt, by default, the mean expression
%                       levels of each gene across all tissues in arrayData
%                       will be used as the threshold values)
%       singleCells     binary value selecting whether to use the
%                       single-cell algorithm to identify expressed genes.
%                       If used, specify cell subpopulations in CELLTYPES
%                       (opt, default [])
%       plotResults     true if single cell probability distributions
%                       should be plotted (opt, default = false)
%   metabolomicsData    cell array with metabolite names that the model
%                       should produce (opt, default [])
%   removeGenes         if true, low-abundance genes will be removed from
%                       grRules, unless they are the only gene associated
%                       with a reaction, or a subunit of an enzyme complex
%                       (see "removeLowScoreGenes" function for details).
%                       If false, grRules will not be modified; however,
%                       genes that were associated only with removed
%                       reactions will not be present in the final model.
%                       (opt, default true).
%   taskFile            a task list in Excel format. See parseTaskList for
%                       details (opt, default [])
%   useScoresForTasks   true if the calculated reaction scored should be
%                       used as weights when fitting to tasks (opt, default
%                       true)
%   printReport         true if a report should be printed to the screen
%                       (opt, default true)
%   taskStructure       task structure generated from parseTaskList. Can be
%                       used as an alternative way to define tasks when Excel
%                       sheets are not suitable. Overrides taskFile (opt,
%                       default [])
%   params              parameter structure as used by getMILPParams. This
%                       is for the INIT algorithm. For the the MILP
%                       problems solved to fit tasks, see paramsFT (opt,
%                       default [])
%   paramsFT            parameter structure as used by getMILPParams. This
%                       is for the fitTasks step. For the INIT algorithm,
%                       see params (opt, default [])
%
%   model                   the resulting model structure
%   metProduction           array that indicates which of the
%                           metabolites in metabolomicsData that could be
%                           produced. Note that this is before the
%                           gap-filling process to enable defined tasks. To
%                           see which metabolites that can be produced in
%                           the final model, use canProduce.
%                           -2: metabolite name not found in model
%                           -1: metabolite found, but it could not be produced
%                           1: metabolite could be produced
%   essentialRxnsForTasks   cell array of the reactions which were
%                           essential to perform the tasks
%   addedRxnsForTasks       cell array of the reactions which were added in
%                           order to perform the tasks
%   deletedDeadEndRxns      cell array of reactions deleted because they
%                           could not carry flux (INIT requires a
%                           functional input model)
%   deletedRxnsInINIT       cell array of the reactions which were deleted by
%                           the INIT algorithm
%   taskReport              structure with the results for each task
%   	id                  cell array with the id of the task
%       description         cell array with the description of the task
%       ok                  boolean array with true if the task was successful
%       essential           cell array with cell arrays of essential
%                           reactions for the task
%       gapfill             cell array of cell arrays of reactions included
%                           in the gap-filling for the task
%
%   This is the main function for automatic reconstruction of models based
%   on the INIT algorithm (PLoS Comput Biol. 2012;8(5):e1002518). Not all
%   settings are possible using this function, and you may want to call the
%   functions scoreComplexModel, runINIT and fitTasks individually instead.
%
%   NOTE: Boundary metabolites should normally not be removed from the model
%   when using this approach, since checkTasks/fitTasks rely on putting specific
%   constraints for each task. The INIT algorithm will remove boundary metabolites
%   if any are present. Use the addBoundaryMets function to add boundary
%   metabolites to a model.
%
%   Usage: [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks,...
%               deletedDeadEndRxns, deletedRxnsInINIT, taskReport] = ...
%               getINITModel2(refModel, tissue, celltype, hpaData, arrayData,...
%               metabolomicsData, removeGenes, taskFile, useScoresForTasks, ...
%               printReport, taskStructure, params, paramsFT);
%

if nargin < 5
    arrayData = [];
end
if nargin < 6
    metabolomicsData = [];
end
if nargin < 7 || isempty(removeGenes)
    removeGenes = true;
end
if nargin < 8
    taskFile = [];
end
if nargin < 9 || isempty(useScoresForTasks)
    useScoresForTasks = true;
end
if nargin < 10 || isempty(printReport)
    printReport = true;
end
if nargin < 11
    taskStructure = [];
end
if nargin < 12
    params = [];
end
if nargin < 13
    paramsFT = [];
end

%Check that the model is in the closed form
if ~isfield(refModel,'unconstrained')
    EM = 'Boundary metabolites should normally be present in the model when using getINITModel2. Use addBoundaryMets(model) to add boundary metabolites to a model';
    dispEM(EM);
end

%Create the task structure if not supplied
if any(taskFile) && isempty(taskStructure)
    taskStructure = parseTaskList(taskFile);
end

% sc-tINIT to identify confidence levels of gene expression
if isfield(arrayData,'singleCells')
    if arrayData.singleCells == 1
        % Check to ensure cell type is defined
        if ~isfield(arrayData,'celltypes')
            dispEM('arrayData must contain cell type information if sc-tINIT is to be used','false');
        end
        if ~ismember(upper(celltype),upper(arrayData.celltypes))
            dispEM('The cell type name does not match');
        end

        % Analyze only cell type of interest
        J = strcmpi(arrayData.celltypes,celltype);

        % Analyze only genes included in the reference model
        I = ismember(arrayData.genes,refModel.genes);

        % Convert expression data to population fractions
        binary_levels = arrayData.levels(I,J) ~= 0; % Classify each gene as detected (1) or not (0) in each cell
        cell_count_levels = sum(binary_levels,2); % Number of cells expressing each transcript
        cell_frac_levels = cell_count_levels/size(binary_levels,2); % Number of cells expressing each transcript

        % Bin cell_frac_counts manually
        x = 0:.01:1;
        for i = 1:length(x)
            cell_frac_count(i) = sum(round(cell_frac_levels,2) == x(i));
        end

        % Fit four beta distributions
        cell_frac_count(cell_frac_count == 0) = NaN; % Remove zeros from optimization
        cell_frac_count(1) = NaN; % Remove non-expressed genes from optimization
        x_lim = 1; % Somewhat arbitrary parameter to fit left tail of distr.
        myfun = @(par) nansum((cell_frac_count(1:find(x >= x_lim,1)) - ...
            abs(par(1))*betapdf(x(1:find(x >= x_lim,1)),abs(par(2)),abs(par(3))) - ...
            abs(par(4))*betapdf(x(1:find(x >= x_lim,1)),abs(par(5)),abs(par(6))) - ...
            abs(par(7))*betapdf(x(1:find(x >= x_lim,1)),abs(par(8)),abs(par(9))) - ...
            abs(par(10))*betapdf(x(1:find(x >= x_lim,1)),abs(par(11)),abs(par(12)))).^2);

        par0 = [4,2,100,7,2,30,7,5,20,5,15,20];
        opts = optimset('Display','off');
        [par,f_val] = fminsearch(myfun,par0,opts);
        par = abs(par);

        % Plot results
        if (isfield(arrayData,'plotResults'))
            if arrayData.plotResults == true
                figure(); hold on; plot(x,cell_frac_count,'ko','MarkerSize',5);
                plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))),'b-','LineWidth',1)
                plot(x,abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))),'b-','LineWidth',1)
                plot(x,abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))),'b-','LineWidth',1)
                plot(x,abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'b-','LineWidth',1)
                plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))) + ...
                    abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))) + ...
                    abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))) + ...
                    abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'-','Color',[.5 .5 .5],'LineWidth',2)
                xlabel('Expression Probability');ylabel('# of genes');set(gca,'FontSize',14,'LineWidth',1.25);
                title('Expression prediction','FontSize',18,'FontWeight','bold')
            end
        end

        % Score genes based on population expression (p = .05)
        exprs_cutoff_1 = find(cdf('beta',x,par(2),par(3)) >.95,1)-1; % Find index of no confidence genes
        exprs_cutoff_2 = find(cdf('beta',x,par(5),par(6)) >.95,1)-1; % Find index of low confidence genes
        exprs_cutoff_3 = find(cdf('beta',x,par(8),par(9)) >.95,1)-1; % Find index of low confidence genes
        exprs_cutoffs = sort([exprs_cutoff_1,exprs_cutoff_2,exprs_cutoff_3]);
        gene_scores = cell_frac_levels*0;
        gene_scores(cell_frac_levels <= x(exprs_cutoffs(1))) = 4; % Not detected
        gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(1))).*(cell_frac_levels < x(exprs_cutoffs(2))))) = 3; % Low detection
        gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(2))).*(cell_frac_levels < x(exprs_cutoffs(3))))) = 2; % Medium detection
        gene_scores(cell_frac_levels > x(exprs_cutoffs(3))) = 1; % High detection

        % Replace hpaData with singleCellData
        if printReport == true
            dispEM('Single cell data is not currently compatible with HPA data. \n         Replacing hpaData with single cell-based scoring.',false);
        end
        hpaData.genes = arrayData.genes;
        hpaData.tissues = arrayData.tissues;
        hpaData.celltypes = arrayData.celltypes;
        hpaData.levels = [{'High'},{'Medium'},{'Low'},{'None'}];
        hpaData.gene2Level = zeros(length(arrayData.genes),length(arrayData.celltypes));
        for i = 1:length(find(J))
            find_var = find(J,i);
            hpaData.gene2Level(I,find_var(end)) = gene_scores;
        end

        % Remove arrayData from the analysis (Might be a bad idea)
        clear arrayData
        arrayData = [];
    end
end

if printReport == true
    if any(celltype)
        fprintf(['***Generating model for: ' tissue ' - ' celltype '\n']);
    else
        fprintf(['***Generating model for: ' tissue '\n']);
    end
    if ~isempty(hpaData)
        fprintf('-Using HPA data\n');
    end
    if ~isempty(arrayData)
        fprintf('-Using array data\n');
    end
    if ~isempty(metabolomicsData)
        fprintf('-Using metabolomics data\n');
    end
    if ~isempty(taskFile) || ~isempty(taskStructure)
        fprintf('-Using metabolic tasks\n');
    end
    fprintf('\n');

    printScores(refModel,'Reference model statistics',hpaData,arrayData,tissue,celltype);
end

%Remove dead-end reactions to speed up the optimization and to
%differentiate between reactions removed by INIT and those that are
%dead-end
if isfield(arrayData,'deletedDeadEndRxns')                  % temporary shortcut, should be revised.
    deletedDeadEndRxns = arrayData.deletedDeadEndRxns;      % temporary shortcut, should be revised.
elseif isfield(hpaData,'deletedDeadEndRxns')                % temporary shortcut, should be revised.
    deletedDeadEndRxns = hpaData.deletedDeadEndRxns;        % temporary shortcut, should be revised.
else                                                        % temporary shortcut, should be revised.
    [~, deletedDeadEndRxns]=simplifyModel(refModel,true,false,true,true,true);
end                                                         % temporary shortcut, should be revised.
cModel = removeReactions(refModel,deletedDeadEndRxns,false,true);

%Store the connected model like this to keep track of stuff
if printReport == true
    printScores(cModel,'Pruned model statistics',hpaData,arrayData,tissue,celltype);
end

%If tasks have been defined, then go through them and get essential
%reactions
if ~isempty(taskStructure)

    if isfield(arrayData,'taskReport') && isfield(arrayData,'essentialRxnMat')  % temporary shortcut, should be revised.
        taskReport = arrayData.taskReport;                                      % temporary shortcut, should be revised.
        essentialRxnMat = arrayData.essentialRxnMat;                            % temporary shortcut, should be revised.
    elseif isfield(hpaData,'taskReport') && isfield(hpaData,'essentialRxnMat')  % temporary shortcut, should be revised.
        taskReport = hpaData.taskReport;                                        % temporary shortcut, should be revised.
        essentialRxnMat = hpaData.essentialRxnMat;                              % temporary shortcut, should be revised.
    else                                                                        % temporary shortcut, should be revised.
        [taskReport, essentialRxnMat]=checkTasks(cModel,[],printReport,true,true,taskStructure);
    end                                                                         % temporary shortcut, should be revised.

    essentialRxnsForTasks = cModel.rxns(any(essentialRxnMat,2));

    %Find metabolites present in taskStruct. We want to avoid removing
    %these metabolites from the final model (even though they may no longer
    %participate in any reacitons) so that the final model is still able to
    %complete all of the tasks without any errors.
    taskMets = union(vertcat(taskStructure.inputs),vertcat(taskStructure.outputs));
    taskMets = union(taskMets, parseRxnEqu(vertcat(taskStructure.equations)));
    modelMets = strcat(cModel.metNames,'[',cModel.comps(cModel.metComps),']');
    [inModel,metInd] = ismember(taskMets,modelMets);
    essentialMetsForTasks = cModel.mets(metInd(inModel));

    %Remove tasks that cannot be performed
    taskStructure(taskReport.ok == false) = [];
    if printReport == true
        printScores(removeReactions(cModel,setdiff(cModel.rxns,essentialRxnsForTasks),true,true),'Reactions essential for tasks',hpaData,arrayData,tissue,celltype);
    end
else
    essentialRxnsForTasks = {};
    essentialMetsForTasks = {};
end

% Score the connected model
rxnScores = scoreComplexModel(cModel,hpaData,arrayData,tissue,celltype);

%Run the INIT algorithm. The exchange reactions that are used in the final
%model will be open, which doesn't fit with the last step. Therefore
%delete reactions from the original model instead of taking the output. The
%default implementation does not constrain reversible reactions to only
%carry flux in one direction. Runs without the constraints on reversibility
%and with all output allowed. This is to reduce the complexity of the
%problem.
[~, deletedRxnsInINIT, metProduction] = runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
initModel = removeReactions(cModel,deletedRxnsInINIT,false,true);

% remove metabolites separately to avoid removing those needed for tasks
unusedMets = initModel.mets(all(initModel.S == 0,2));
initModel = removeMets(initModel, setdiff(unusedMets,essentialMetsForTasks));

if printReport == true
    printScores(initModel,'INIT model statistics',hpaData,arrayData,tissue,celltype);
    printScores(removeReactions(cModel,setdiff(cModel.rxns,deletedRxnsInINIT),true,true),'Reactions deleted by INIT',hpaData,arrayData,tissue,celltype);
end

%The full model has exchange reactions in it. fitTasks calls on fillGaps,
%which automatically removes exchange (boundary) metabolites (because it
%assumes that the reactions are constrained when appropriate). In this case
%the uptakes/outputs are retrieved from the task sheet instead. To prevent
%exchange reactions being used to fill gaps, they are delete from the
%reference model here.
initModel.id = 'INITModel';

%If gaps in the model should be filled using a task list
if ~isempty(taskStructure)
    %Remove exchange reactions and reactions already included in the INIT
    %model
    refModelNoExc = removeReactions(refModel,union(initModel.rxns,getExchangeRxns(refModel)),true,true);

    %At this stage the model is fully connected and most of the genes with
    %good scores should have been included. The final gap-filling should
    %take the scores of the genes into account, so that "rather bad"
    %reactions are preferred to "very bad" reactions. However, reactions
    %with positive scores will be included even if they are not connected
    %in the current formulation. Therefore, such reactions will have to be
    %assigned a small negative score instead.
    if useScoresForTasks == true
        refRxnScores = scoreComplexModel(refModelNoExc,hpaData,arrayData,tissue,celltype);
        [outModel,addedRxnMat] = fitTasks(initModel,refModelNoExc,[],true,min(refRxnScores,-0.1),taskStructure,paramsFT);
    else
        [outModel,addedRxnMat] = fitTasks(initModel,refModelNoExc,[],true,[],taskStructure,paramsFT);
    end
    if printReport == true
        printScores(outModel,'Functional model statistics',hpaData,arrayData,tissue,celltype);
        printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,arrayData,tissue,celltype);
    end

    addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2));
else
    outModel = initModel;
    addedRxnMat = [];
    addedRxnsForTasks = {};
end

% The model can now perform all the tasks defined in the task list.
model = outModel;

% If requested, attempt to remove negative-score genes from the model,
% depending on their role (isozyme or complex subunit) in each grRule.
% See the "removeLowScoreGenes" function more more details, and to adjust
% any default parameters therein.
if ( removeGenes )
    [~, geneScores] = scoreComplexModel(model,hpaData,arrayData,tissue,celltype);
    model = removeLowScoreGenes(model,geneScores);
end

% At this stage the model will contain some exchange reactions but probably
% not all (and maybe zero). This can be inconvenient, so all exchange
% reactions from the reference model are added, except for those which
% involve metabolites that are not in the model.

% First delete any included exchange reactions in order to prevent the order
% from changing
model = removeReactions(model,getExchangeRxns(model));

% Create a model with only the exchange reactions in refModel
excModel = removeReactions(refModel,setdiff(refModel.rxns,getExchangeRxns(refModel)),true,true);

% Find the metabolites there which are not boundary metabolites and which do
% not exist in the output model
I = ~ismember(excModel.mets,model.mets) & excModel.unconstrained == 0;

% Then find those reactions and delete them
[~, J] = find(excModel.S(I,:));
excModel = removeReactions(excModel,J,true,true);

% Merge with the output model
model = mergeModels({model;excModel},'metNames');
model.id = 'INITModel';
model.description = ['Automatically generated model for ' tissue];
if any(celltype)
    model.description = [model.description ' - ' celltype];
end

if printReport == true
    printScores(model,'Final model statistics',hpaData,arrayData,tissue,celltype);
end

% Add information about essential reactions and reactions included for
% gap-filling and return a taskReport
if ~isempty(taskStructure)
    I = find(taskReport.ok); %Ignore failed tasks
    for i = 1:numel(I)
        taskReport.essential{I(i),1} = cModel.rxns(essentialRxnMat(:,I(i)));
        taskReport.gapfill{I(i),1} = refModelNoExc.rxns(addedRxnMat(:,i));
    end
else
    taskReport = [];
end

%Fix grRules and reconstruct rxnGeneMat
[grRules,rxnGeneMat] = standardizeGrRules(model,true);
model.grRules = grRules;
model.rxnGeneMat = rxnGeneMat;
end

%This is for printing a summary of a model
function [rxnS, geneS] = printScores(model,name,hpaData,arrayData,tissue,celltype)
[a, b] = scoreComplexModel(model,hpaData,arrayData,tissue,celltype);
rxnS = mean(a);
geneS = mean(b,'omitnan');
fprintf([name ':\n']);
fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
end

author: SysBioChalmers

repo: Human-GEM


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Locally Weighted Regression 2D Example  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%          1) Generate 2D Regression Datasets                %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generate target function to learn from and visualise
clear all; close all; clc;
f       = @(x,y)sin(x).*cos(y);
r       = @(a,b,N,M)a + (b-a).*rand(N,M);
N       = 2000;
X       = r(-5,5,N,2);
y       = f(X(:,1),X(:,2));

% Plot data
options             = [];
options.points_size = 15;
options.cmap        = y;
options.plot_labels = {'$x_1$','$x_2$','y'};
options.title       = 'Example 2D Non-linear Training Data';
if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data([X,y],options); hold on;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 2) "Train" LWR Model                       %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Train LWR
options             = [];
options.dim         = 2;
options.bUseKDT     = true;
options.D           = ones(2,1) .* ((0.2).^2);
options.K           = 100;
lwr                 = LWR(options);

% Plot Train Data and train LWR
lwr.train(X,y);

f = @lwr.f;

% Plot Estimated Function
options           = [];
options.title     = 'Estimated y=f(x) from Locally Weighted Regression';
options.surf_type = 'surf';
if exist('h2','var') && isvalid(h2), delete(h2);end
h2 = ml_plot_value_func(X,f,[1 2],options);hold on

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           3) Test LWR Model on Query points                %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Predict values on a smaller grid x_1:[-2,2] x_2:[-2,2]

[x1_test,x2_test] = meshgrid(linspace(-2,2,50),linspace(-2,2,50));
y_test     = lwr.f([x1_test(:),x2_test(:)]);
M = length(y_test);

% Plotting Options for Regressive Function
options           = [];
options.title     = 'Estimated y=f(x) from Locally Weighted Regression';
options.surf_type = 'surf';
if exist('h3','var') && isvalid(h3), delete(h3);end
h3 = ml_plot_value_func([x1_test(:),x2_test(:)],f,[1 2],options);hold on

% Plot Training Data on Top
options = [];
options.plot_figure = true;
options.points_size = 12;
options.labels = zeros(M,1);
options.plot_labels = {'$x_1$','$x_2$','y'};
ml_plot_data([x1_test(:),x2_test(:),y_test(:)],options);

author: epfl-lasa


text:
function [ fea, out ] = ex_convdiff2( varargin )
%EX_CONVDIFF2 1D Time dependent convection and diffusion equation example.
%
%   [ FEA, OUT ] = EX_CONVDIFF2( VARARGIN ) 1D time dependent convection and diffusion equation on
%   a line with exact solution exp(-k^2*nu*t)*sin(k*(x-a*t)) and periodic boundary conditions.
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       k           scalar {2*pi}          Simulation parameter
%       a           scalar {1}             Convection velocity
%       nu          scalar {0.1}           Diffusion coefficient
%       hmax        scalar {1/25}          Max grid cell size
%       dt          scalar {0.01}          Time step size
%       ischeme     scalar {2}             Time stepping scheme
%       sfun        string {sflag1}        Shape function
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'k',        2*pi; ...
  'a',        1; ...
  'nu',       0.1; ...
  'hmax',     1/25; ...
  'dt'        0.01; ...
  'ischeme'   2; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'tol',      1e-1; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

refsol = ['exp(-',num2str(opt.k^2*opt.nu),'*t)*sin(',num2str(opt.k),'*(x-',num2str(opt.a),'*t))'];

% Grid generation.
fea.grid = linegrid( 1/opt.hmax, 0, 1 );

% Problem definition.
fea.sdim  = { 'x' };
fea = addphys( fea, @convectiondiffusion );
fea.phys.cd.sfun = { opt.sfun };
fea.phys.cd.eqn.coef{2,4} = { opt.nu };
fea.phys.cd.eqn.coef{3,4} = { opt.a  };
fea = parsephys(fea);

% Parse and solve problem.
fea = parseprob( fea );

x = fea.grid.p';
n = length(x);
if( strcmp( opt.sfun,'sflag2' ) )
  x = [ x; (x(2:end)+x(1:end-1))/2 ];
end
t  = 0;
u0 = eval( refsol );

% Assembly.
[M,A,f] = assembleprob( fea, 'f_m', 1, 'imass', 1, 'f_a', 1, 'f_f', 1, 'f_sparse', 1 );
M = spdiags( full(sum(M')'), 0, size(M,1), size(M,1) );
fea.sol.u = u0;
dt = opt.dt;
it = 0;
tlist = 0;
tmax = 1;
if( opt.ischeme==2 )       % Crank-Nicolson.
  C = [ M + dt/2*A ];
elseif( opt.ischeme==1 )   % Backward Euler.
  C = [ M + dt*A ];
end
v0 = zeros(size(C,1),1);
v0(1) =  1;
v0(n) = -1;
C = [C v0; v0' 0];

% Solver loop.
while 1
  t  = t  + dt;
  it = it + 1;

  u_r = eval( refsol );
  if( opt.ischeme==2 )       % Crank-Nicolson.
    b = [ [ M - dt/2*A ]*u0 + dt*f ];
  elseif( opt.ischeme==1 )   % Backward Euler.
    b = [ M*u0 + dt*f ];
  end

  u1 = C\[b;0];
  u1(end) = [];

  if( t>=tmax )
    break
  end

  err = norm( u1 - u_r )/norm( u_r );
  errnm(it) = err;
  if( ~isempty(fid) )
    fprintf( fid, 'Time = %f, error norm = %d\n', t, err );
  end

  fea.sol.u = [ fea.sol.u u1 ];
  tlist = [ tlist t ];
  u0 = u1;
end

% Postprocessing.
if( opt.iplot>0 )
  figure;
  if( opt.iplot>1 )
    i_sol_list = 1:numel(tlist);
  else
    i_sol_list = numel(tlist);
  end
  [~,ix] = sort( x );
  for i_sol=i_sol_list
    t = tlist(i_sol);
    clf
    postplot( fea, 'surfexpr', 'c', 'solnum', i_sol );
    hold on
    u_r = eval( refsol );
    plot( sort(x), u_r(ix), 'r--' );
    title( ['Solution at time ',num2str(t)])
    xlabel( 'x' )
    drawnow
  end
end

% Error checking.
out.err  = errnm;
out.pass = all( errnm<opt.tol );

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
i_fish=8;
% also load fish 8 in GUI
% change alpha_max to 1 in 'DrawCellsOnAnatProj.mat'

CR = getappdata(hfig,'CellResp');

%% plot basic stats of time series
TS = CR(:,1:3000);
TS_max = max(TS,[],2);
TS_avr = mean(TS,2);
figure;
hist(TS_max,0:0.1:5)
figure
hist(TS_avr,0:0.01:1)

%% color-code based on average activity
cIX_plot = 1:size(CR,1);
TS_avr = mean(TS,2);

TSmax = 0.3;
TSmin = min(TS_avr);
TS_avr(TS_avr>TSmax) = TSmin;

score = (TS_avr-TSmin)/(TSmax-TSmin);
gIX_plot = ceil(score*64);
gIX_plot(gIX_plot==0)=1;
wIX = score.^1.8; % set transparency weight, and enhance contrast

numK = max(gIX_plot);
clrmap = parula(numK);

setappdata(hfig,'isWeighAlpha','1');
setappdata(hfig,'wIX',wIX);
figure
DrawCellsOnAnatProj(hfig,0,1,cIX_plot,gIX_plot,clrmap);

author: xiuyechen

repo: FishExplorer


text:
function [flag, rect] = Judge_Direction(bw)
[L, num] = bwlabel(bw);
stats = regionprops(bw, 'Area', 'BoundingBox');
Area = cat(1, stats.Area);
[Area, ind] = sort(Area, 'descend');
if num == 1
    rect = stats.BoundingBox;
else
    rect1 = stats(ind(1)).BoundingBox;
    rect2 = stats(ind(2)).BoundingBox;
    s1 = [rect1(1); rect2(1)];
    s2 = [rect1(2); rect2(2)];
    s = [min(s1) min(s2) rect1(3)+rect2(3) rect1(4)+rect2(4)];
    rect = s;
end
rate = rect(3)/rect(4);
if rate > 1
    flag = 1;
else
    flag = 2;
end

author: decouples


text:
function compute_label_map(varargin)
% COMPUTE_LABEL_MAP Determines the ILVSRC mapping used in SENets
%   COMPUTE_LABEL_MAP aims to invert the ILSVRC mapping from the
%   order used to train the public SENEt models to the standard
%   ASCII-ordering convention
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  opts.numClasses = 1000 ;
  opts.dataDir = fullfile(vl_rootnn, 'data/datasets/ILSVRC2012') ;
  opts.miscDir = fullfile(vl_rootnn, 'contrib/mcnSENets/misc') ;
  opts = vl_argparse(opts, varargin) ;

  imdbPath = fullfile(vl_rootnn, 'data', 'imagenet12', 'imdb.mat');

  % set paths to source and storage destination for label map
  valFile = fullfile(opts.miscDir, 'ILSVRC2017_val.txt') ;
  mapFile = fullfile(opts.miscDir, 'label_map.txt') ;

  if exist(imdbPath, 'file')
    imdb = load(imdbPath) ;
    imdb.imageDir = fullfile(opts.dataDir, 'images');
  else
    imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
    mkdir(opts.expDir) ;
    save(imdbPath, '-struct', 'imdb') ;
  end

  % read examples
  fid = fopen(valFile, 'r') ; data = textscan(fid, '%s %d') ; fclose(fid) ;
  gtImNames = data{1} ; gtLabels = data{2} + 1 ; % fix offset

  % trim stems for comparison
  valIdx = imdb.images.set == 2 ; valNames = imdb.images.name(valIdx) ;
  imNames = cellfun(@getImName, valNames, 'Uni', 0) ;
  valLabels = imdb.images.label(valIdx) ;

  % create map - favour safety over speed here
  labelMap = zeros(1, opts.numClasses) ;
  for ii = 1:numel(imNames)
    imName = imNames{ii} ;
    [member, idx] = ismember(imName, gtImNames) ;
    assert(member, 'unrecognised val image') ;
    srcLabel = valLabels(ii) ; destLabel = gtLabels(idx) ;
    if labelMap(srcLabel) ~= 0
      assert(labelMap(srcLabel) == destLabel, 'inconsistent label map') ;
    else
      labelMap(srcLabel) = destLabel ;
      fprintf('(%d/%d) remapping %d\n', ii, numel(imNames), srcLabel) ;
    end
  end

  % check for completeness
  assert(numel(unique(labelMap)) == numel(labelMap), 'incomplete mapping') ;

  % store label map to disk
  fid = fopen(mapFile, 'w') ; fprintf(fid, '%d\n', labelMap) ; fclose(fid) ;

% -------------------------------
function imName = getImName(path)
% -------------------------------
  [~,stem] = fileparts(path) ;
  imName = strcat(stem, '.JPEG') ;

repo: MFT


text:
%Test kinematics

%Test robot at this position
px=0.0;
py=0.0;
pz=0.6;
%pz=1;

T=eye(4);
T(1,4)=px;
T(2,4)=py;
T(3,4)=pz;

%call inverse kinematics in position
q=inversekinematic(robot, T);
robot.robot1.graphical.draw_transparent=1;
robot.robot2.graphical.draw_transparent=1;
robot.robot3.graphical.draw_transparent=1;

%draw the robot!
drawrobot3d(robot, q(:,8))

%define position as computed by the inverse kinematic function. Only active
%joints
th1=q(1,8);
th2=q(4,8);
th3=q(7,8);

%Joint speeds
th1d=0;
th2d=0;
th3d=0;

%Joint accelerations
th1dd=0;
th2dd=0;
th3dd=0;

% This inverse dynamic model uses a Lagrangian formulation
% thus, we need to use a different approach to define the kinematic state
% of the arm, we first solve the directkinematic problem in position

% Next, we solve it for speed and acceleration
% Two different approaches are here valid. Consider first that we want the
% mechanism to be placed at q1, q2 and q3 known, and we want each joint to
% accelerate at a rate qdd1, qdd2, qdd3. The Lagrangian formulation is
% based on a redundant set of coordinates, thus, we need to specify the
% position px, py, pz of the end effector. Also, the acceleration of the
% end effector pddx pddy pddz must be specified
%[T, pxyzd, pxyzdd]=directkinematic_Maryland3DOF(robot, q(:,8), qd, qdd)
[T, pxyzd, pxyzdd]=directkinematic_Maryland3DOF(robot, [th1 0 0 th2 0 0 th3 0 0], [th1d 0 0 th2d 0 0 th3d 0 0], [th1dd 0 0 th2dd 0 0 th3dd 0 0]);

%We then use pxdd, pydd, pzdd, th1, th2, th3, th1dd,th2dd,th3dd
%tau=inversedynamic_Maryland3DOF(robot, [px, py, pz, th1 th2 th3], [0 0 0 0 0 0 0 0 0], [0 0 0 0 0 0 0 0 0], [0 0 0])
% You can also specify a set of external forces applied at the end
% effector. These forces are specified in the base reference system.
% Typically, the force exerted by a load should be specified as -m*g being
% g=9.81 m/s^2. For example fext = [0 0 -10*9.81] for a 10kg load placed on
% the end effector's centre.
%
% gc=9.81 if the robot is placed upwards
% gc=-9.81 if the robot is hanging.
%
% Consider a load of 10kg placed on the end effector. Write:
%   fext = [0 0 -10*9.81] if the robot is placed upwards, or
%   fext = [0 0 10*9.81] if the robot is hanging.

%Example 1: robot placed upwards
gc = 9.81;
fext = [0 0 -10*9.81];
tau=inversedynamic_Maryland3DOF(robot, [px, py, pz, th1 th2 th3], [pxyzdd(1) pxyzdd(2) pxyzdd(3) th1dd th2dd th3dd], gc, fext)

%Example 2: robot hanging
gc = -9.81;
fext = [0 0 10*9.81];
tau=inversedynamic_Maryland3DOF(robot, [px, py, pz, th1 th2 th3], [pxyzdd(1) pxyzdd(2) pxyzdd(3) th1dd th2dd th3dd], gc, fext)

repo: ARTE


text:
function [varargout, peval] = separcomp(dpix, peval, winit_pix, hinit, verbose)
% separcomp(dpix, p_script, savethis, winit, hinit, verbose)
% separate components
% V ~ WH
% V -> N_pix x N_t
% W -> N_pix x N_comp - xth pixel of the ith components
% H -> N_copm x N_t - contribution of the i-th component in the time t
tstart = tic;
if ~exist('verbose','var')
    verbose = 1;
end

if ~isfield(peval, 'w_fixvec') peval.w_fixvec=[]; end
if ~isfield(peval, 'h_fixvec') peval.h_fixvec=[]; end

% Fixing background components:
if ~isfield(peval, 'fixbackground')
    peval.fixbackground = 'wh';
end

% Fixing w component of the background
if sum(peval.fixbackground == 'w')>0
    if ~sum(ismember(peval.ncomp+1, peval.w_fixvec))
        peval.w_fixvec = [peval.w_fixvec, peval.ncomp + 1];
    end
end

% Fixing h component of the background
if sum(peval.fixbackground == 'h')>0
    if ~sum(ismember(peval.ncomp+1, peval.h_fixvec))
        peval.h_fixvec = [peval.h_fixvec, peval.ncomp + 1];
    end
end

[peval.nx, peval.ny, peval.nt] = size(dpix);
peval.numpix = peval.nx*peval.ny;
peval.meandata = mean(dpix(:));

dvec = reshape(dpix,peval.numpix, peval.nt);

% Background subtraction with empirical values:
if ~isfield(peval, 'bg')
    [out_nobg, peval]=backgroundestimation(dpix, peval);
end

% Estimation of the number of componenets
if ~isfield(peval, 'ncomp')
    peval.ncomp=estimate_ncpomp(dvec);
end

% Initialization of w and h
[winit, hinit] = initwh(winit_pix, hinit, peval, verbose);

% This is to fix zeros in the data as the divergence is not defined for that...
dvec=max(dvec, eps);

eval (['[w,h,peval, dtrace, htrace]=' peval.method '(dvec,winit,hinit,peval,verbose);']);
peval.ddiv_end = ddivergence(dvec, w*h);
peval.gradddiv_end = dtrace(end-1)-dtrace(end);
%better with feval....
% varargout = struct('w',w,'h',h, 'dtrace', dtrace, 'htrace', htrace);
varargout = struct('w',w,'h',h);
peval.elapsedtime_separcomp = toc(tstart);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [pb1, pb2, theta] = getOrientedOcclusionProbs(bndinfo)
% [pb1, pb2, theta] = getOrientedOcclusionProbs(bndinfo)

if isstruct(bndinfo)
    bndinfo = bndinfo{1};
end

N = numel(bndinfo);
pocc = zeros([bndinfo{1}.imsize N], 'single');
pocc_left = zeros([bndinfo{1}.imsize N], 'single');

for k = 1:N

    tmpim = zeros(bndinfo{1}.imsize);
    tmpim1 = zeros(bndinfo{1}.imsize);

    if isfield(bndinfo{k}, 'result')
        pbk = 1-bndinfo{k}.result.edgeProb(:, 1);
        pbk1 = bndinfo{k}.result.edgeProb(:, 2) ./ pbk;
    elseif isfield(bndinfo{k}, 'pbnd')
        ne = bndinfo{k}.ne;
        pbk = 1-bndinfo{k}.pbnd(1:ne, 1);
        pbk1 = bndinfo{k}.pbnd(1:ne, 2) ./ pbk;
    else
        error('pb info not found')
    end

    for e = 1:bndinfo{k}.ne
        ind = bndinfo{k}.edges.indices{e};
        tmpim(ind(2:end-1)) = pbk(e); %tmpim(ind(2:end-1)) + pbk(e);
        tmpim(ind([1 end])) = max(tmpim(ind([1 end])), pbk(e)); % avoid overcounting junctions
        tmpim1(ind(1:end-1)) = pbk1(e);
        %tmpim1(ind([1 end])) = tmpim1(ind([1 end])) + pbk1(e); % take average for junctions
    end

    pocc(:, :, k) = tmpim;
    pocc_left(:, :, k) = tmpim1;
end

theta = bndinfo{1}.edges.thetaDirected;

% compute average pb for left and right side
ind = getBoundaryCenterIndices(bndinfo{1});
pb1 = mean(pocc.*pocc_left, 3);
pb1 = pb1(ind);
pb2 = mean(pocc.*(1-pocc_left), 3);
pb2 = pb2(ind);

author: Cloud-CV

repo: object-proposals


text:
function [ TreeMatric,Cost ] =  DirectedMinimalSpanningTree( OriginalCostMatric,Root )

    MaxTree =  DirectedMaximumSpanningTree( -OriginalCostMatric,Root );
    TreeMatric = -MaxTree;
    Cost = sum(sum(TreeMatric));

end


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

%*****************************************************************************80
%
%% WEDGE01_SAMPLE samples points uniformly from the unit wedge in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 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(3,N), the points.
%
  x = zeros ( 3, n );

  for j = 1 : n

    [ e, seed ] = r8vec_uniform_01 ( 4, seed );

    e(1:3) = - log ( e(1:3) );

    e_sum = sum ( e(1:3) );

    x(1,j) = e(1,1) / e_sum;
    x(2,j) = e(2,1) / e_sum;
    x(3,j) = 2.0 * e(4) - 1.0;

  end

  return
end


text:
function d =  generate(a)

x1=zeros(a.m,1);
y1=zeros(a.m,1);

x2=zeros(a.m,1);
y2=zeros(a.m,1);

step=1/100;
t=step:step:a.n*2*pi;

x1= t.*cos(t) ;
y1= t.*sin(t);
x2= (1+t).*cos(t) ;
y2= (1+t).*sin(t) ;

xfinal1=[];
yfinal1=[];
xfinal2=[];
yfinal2=[];

nn=a.noise;
 for j=1:a.m
     i=ceil(length(t)*rand);
     xfinal1=[xfinal1,x1(i)+nn*rand];
     yfinal1=[yfinal1,y1(i)+nn*rand];
     xfinal2=[xfinal2,x2(i)+nn*rand];
     yfinal2=[yfinal2,y2(i)+nn*rand];
 end

 d=data([get_name(a) ' '] ,[xfinal1' yfinal1';xfinal2' yfinal2'],[ ones(a.m,1);-1*ones(a.m,1)]);

author: canlab

repo: CanlabCore


text:
% The COBRAToolbox: testChangeObjective.m
%
% Purpose:
%     - Tests the changeObjective function
% Author:
%     - Original file: Stefania Magnusdottir

% save the current path
currentDir = pwd;

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

% define the test model
toyModel = struct;
toyModel.rxns = {'Rxn1'; 'Rxn2'; 'Rxn3'};
toyModel.c = [0 ; 0 ; 0];

% test default coefficient (1)
modelNew = changeObjective(toyModel, 'Rxn1');
assert(modelNew.c(1) == 1 && modelNew.c(2) == 0 && modelNew.c(3) == 0)

% test that error is thrown if objective Rxn4 not in model
try
    modelNew = changeObjective(toyModel, 'Rxn4');
catch ME
    assert(length(ME.message) > 0)
end

% test multiple rxns and 1 rxn not in model, default coefficient (1)
modelNew = changeObjective(toyModel,{'Rxn1'; 'Rxn2'; 'Rxn4'});
assert(modelNew.c(1) == 1 && modelNew.c(2) == 1 && modelNew.c(3) == 0)

% test multiple rxns, different coefficients
modelNew = changeObjective(toyModel,{'Rxn1'; 'Rxn2'},[0.3 ; 0.7]);
assert(modelNew.c(1) == 0.3 && modelNew.c(2) == 0.7 && modelNew.c(3) == 0)

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function varargout = setcodebitGUI(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @setcodebitGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @setcodebitGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

%--------------------------------------------------------------------------
function setcodebitGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
set(handles.radiobutton_byreset, 'Value', 1)
set(handles.radiobutton_bybit, 'Value', 0)
set(handles.togglebutton_lowerbyte, 'Enable', 'on')
set(handles.togglebutton_upperbyte, 'Enable', 'on')
set(handles.edit_set2zero, 'Enable', 'off')
set(handles.edit_set2one, 'Enable', 'off')

version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   Set code bits GUI'])

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

% help
helpbutton

uiwait(handles.gui_chassis);

%--------------------------------------------------------------------------
function varargout = setcodebitGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
delete(handles.gui_chassis);
pause(0.1)

%--------------------------------------------------------------------------
function pushbutton_help_Callback(hObject, eventdata, handles)
doc pop_setcodebit

%--------------------------------------------------------------------------
function pushbutton_cancel_Callback(hObject, eventdata, handles)
handles.output = [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);
%--------------------------------------------------------------------------
function pushbutton_ok_Callback(hObject, eventdata, handles)

if get(handles.radiobutton_byreset, 'Value')

        if get(handles.togglebutton_lowerbyte, 'Value') && ~get(handles.togglebutton_upperbyte, 'Value')
                output = {'lower'};
        elseif ~get(handles.togglebutton_lowerbyte, 'Value') && get(handles.togglebutton_upperbyte, 'Value')
                output = {'upper'};
        elseif get(handles.togglebutton_lowerbyte, 'Value') && get(handles.togglebutton_upperbyte, 'Value')
                output = {'all'};
        else
                msgboxText =  'You must press one button at least!';
                title = 'ERPLAB: setcodebitGUI few inputs';
                errorfound(msgboxText, title);
                return
        end
else
        bit2zero = str2num(get(handles.edit_set2zero, 'String'));
        bit2one  = str2num(get(handles.edit_set2one, 'String'));

        if isempty(bit2zero) && isempty(bit2one)
                msgboxText =  'You must specify one valid bit (1 to 16) at least!';
                title = 'ERPLAB: setcodebitGUI few inputs';
                errorfound(msgboxText, title);
                return
        end
        if ~isempty(bit2zero)
                if max(bit2zero)>16 || min(bit2zero)<1
                        msgboxText =  'For bits you must specify any integer between 1 to 16!';
                        title = 'ERPLAB: setcodebitGUI few inputs';
                        errorfound(msgboxText, title);
                        return
                end
        end
        if ~isempty(bit2one)
                if max(bit2one)>16 || min(bit2one)<1
                        msgboxText =  'For bits you must specify any integer between 1 to 16!';
                        title = 'ERPLAB: setcodebitGUI few inputs';
                        errorfound(msgboxText, title);
                        return
                end
        end
        output = {'bit' bit2zero bit2one};
end
handles.output = output;
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function edit_set2zero_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_set2zero_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function edit_set2one_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_set2one_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function radiobutton_byreset_Callback(hObject, eventdata, handles)
if get(hObject, 'Value')
        set(handles.radiobutton_bybit, 'Value', 0)
        set(handles.togglebutton_lowerbyte, 'Enable', 'on')
        set(handles.togglebutton_upperbyte, 'Enable', 'on')
        set(handles.edit_set2zero, 'Enable', 'off')
        set(handles.edit_set2one, 'Enable', 'off')
else
        set(handles.radiobutton_byreset, 'Value', 1)
end

%--------------------------------------------------------------------------
function radiobutton_bybit_Callback(hObject, eventdata, handles)
if get(hObject, 'Value')
        set(handles.radiobutton_byreset, 'Value', 0)
        set(handles.togglebutton_lowerbyte, 'Enable', 'off')
        set(handles.togglebutton_upperbyte, 'Enable', 'off')
        set(handles.edit_set2zero, 'Enable', 'on')
        set(handles.edit_set2one, 'Enable', 'on')
else
        set(handles.radiobutton_bybit, 'Value', 1)
end

%--------------------------------------------------------------------------
function togglebutton_lowerbyte_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function togglebutton_upperbyte_Callback(hObject, eventdata, handles)

% -----------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)

if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')
        %The GUI is still in UIWAIT, us UIRESUME
        handles.output = '';
        %Update handles structure
        guidata(hObject, handles);
        uiresume(handles.gui_chassis);
else
        % The GUI is no longer waiting, just close it
        delete(handles.gui_chassis);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
 function [xs, ni] = tml_convex(x, G, yi, bi, ri, niter, pixmax, type)
%function [xs, ni] = tml_convex(x, G, yi, bi, ri, niter, pixmax, type)
% T-ML-Convex algorithms for transmission Poisson problem
% model: Y_i ~ Poisson(b_i exp(-[G x]_i) + r_i)
% in
%	x	[np,1]	initial guess
%	G	[nd,np]	transpose of system matrix
%	yi	[nd,1]	transmission sinogram
%	bi	[nd,1]	blank scan factors
%	ri	[nd,1]	background (randoms, scatter, crosstalk, etc)
%		bi,ri:	optional (can use empty matrices)
%		yi,bi,ri must have identical dimensions
%	niter		# of iterations
%	pixmax	upper constraint for pixel values
%		can be scalar (e.g. 'inf') or an array the size of x
%	type	'ps'	paraboloidal surrogate
%		'nr2'	newton raphson (2nd choice of zij's)
%		'log1'	-log method with amax
%		'log2'	-log with better choice fo zij's
% out
%	x [np,niter]	updated image vectors each iteration
%
% Copyright 2001-10-17, Jeff Fessler, The University of Michigan

if nargin < 3, ir_usage, end

if ~isvar('bi') || isempty(bi)
	bi = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end
if ~isvar('niter') || isempty(niter)
	niter = 2;
end
if ~isvar('pixmax') || isempty(pixmax)
	pixmax = inf;
end
pixmin = 0;
if ~isvar('type') || isempty(type)
	type = 'ps';
end

trl_check(yi, bi, ri);

% column sums
gi = sum(G')';	% [nd,1] g_i = sum_j g_ij

if streq(type, 'log1')
	step = 1 / max(gi);
	printf('log1 step = %g', step)

elseif streq(type, 'log2')
	gj = sum(G)';				% [np,1] g_j = sum_i g_ij
	tmp = G * (1 ./ gj);			% [nd,1] sum_j g_ij / gj
	fj = max(tmp) * gj;
	tmp = G * (1 ./ fj);			% [nd,1] sum_j g_ij / gj
	tmp = full(tmp);
	printf('check inequality = %g, %g', min(tmp), max(tmp))
	step = 1 ./ fj;
	printf('log2 steps = %g, %g', min(step), max(step))

elseif streq(type, 'nr2')
	warning 'newton can be non-monotone'

elseif ~streq(type, 'ps')
	error 'unknown convex type'
end

% initialize
xs = zeros(length(x), niter);
x = max(x,pixmin);
x = min(x,pixmax);
xs(:,1) = x;

%
% loop over iterations
%
for it=2:niter
	li = reshape(G * x(:), size(yi));	% l=G*x "line integrals"
	bel = bi .* exp(-li);
	yb = bel + ri;			% predicted measurement means

	if streq(type, 'log1') || streq(type, 'log2')
		top = G' * bel(:);
		Ni = bel ./ yb .* yi;
		bottom = G' * Ni(:);

		x = x + step .* log(top ./ bottom);	% the update!

	else
		dothi = (1 - yi ./ yb) .* bel;

		if streq(type, 'ps')
			Ni = yi .* bel ./ yb;			% \bar{N}_i
			ci = trl_curvature(Ni, bi, 0, li, 'oc');	% trick!
		elseif streq(type, 'nr2')
			ci = bel;
		else
			error bug
		end

                denom = G' * (gi .* ci(:));

		x = x + (G' * dothi(:)) ./ denom;	% the update!
	end

	x = max(x,pixmin);			% lower bound
	x = min(x,pixmax);			% upper bound

	xs(:,it) = x;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function b = st_mv ( m, n, nst, ist, jst, ast, x )

%*****************************************************************************80
%
%% ST_MV multiplies an R8SP matrix by an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of
%    the matrix.
%
%    Input, integer NST, the number of nonzero elements in
%    the matrix.
%
%    Input, integer IST(NST), JST(NST), the row and
%    column indices of the nonzero elements.
%
%    Input, real AST(NST), the nonzero elements of the matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(M), the product vector A*X.
%
  b = zeros ( m, 1 );

  for k = 1 : nst
    i = ist(k);
    j = jst(k);
    b(i) = b(i) + ast(k) * x(j);
  end

  return
end


text:
function z_callBack_virtualTeachPendant_J3_minus()
%% Soft realtime joints control Joint X button callback
% Copyright: Mohammad SAFEEA, 27th-July-2018
global virtual_smartPad_w_joints_command;
disp('Moving J3 in negative direction');
virtual_smartPad_w_joints_command(3)=-1;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function value = halham_n_check ( n )

%*****************************************************************************80
%
%% HALHAM_N_CHECK checks N for a Halton or Hammersley sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points in the subsequence.
%
%    Output, logical VALUE, is true if N is legal.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALHAM_N_CHECK - Fatal error!\n' );
    fprintf ( 1, '  Input value of N < 1!\n' );
    fprintf ( 1, '  N = %d\n', n );
    value = 0;
  else
    value = 1;
  end

  return
end