text:
% CHANPROJ - make a detailed plot of data returned from PLOTPROJ
%              for given channel. Returns the data plotted.
% Usage:
%  >> [chandata] = chanproj(projdata,chan);
%  >> [chandata] = chanproj(projdata,chan,ncomps,framelist,limits,title,colors);
%
% Inputs:
%   projdata    = data returned from PLOTPROJ
%   chan        = single channel to plot
%   ncomps      = number of component projections in projdata
%
% Optional:
%   framelist   = data frames to plot per epoch Ex: [1:128] (0|def -> all)
%   limits      = [xmin xmax ymin ymax]  (x's in msec)
%                                        (0, or y's 0 -> data limits)
%   title       = fairly short single-quoted 'plot title'  (0|def -> chan)
%   colors      = file of color codes, 3 chars per line  (NB: '.' = space)
%                                        (0|def -> white is original data)
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1996

% Copyright (C) 1996 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 11-30-96 Scott Makeig  CNL / Salk Institute, La Jolla as plotprojchan.m
% 03-19-97 changed VAR to diag(COV) -sm
% 03-26-97 fix (== -> =) -sm
% 04-03-97 allow framelist to be a col vector, allow 32 traces, fix pvaf, made
%           ncomps mandatory  -sm
% 04-04-97 shortened name to CHANPROJ -sm
% 05-20-97 added read of icadefs.m -sm
% 11-05-97 disallowed white traces unless default axis color is white -sm & ch
% 11-13-97 rm'ed errcode variable -sm
% 12-08-97 added LineWidth 2 to data trace, changed whole plot color to BACKCOLOR -sm
% 01-25-02 reformated help & license -ad

function [chandata] = chanproj(projdata,chan,ncomps,framelist,limits,titl,colorfile)

icadefs; % read default MAXPLOTDATACHANS

if nargin < 7,
    colorfile =0;
end
if nargin < 6,
    titl = 0;
end
if nargin < 5,
    limits = 0;
end
if nargin < 4,
    framelist = 0;
end

if nargin < 3,
    fprintf('chanproj(): requires at least three arguments.\n\n');
    help chanproj
    return
end

[chans,framestot] = size(projdata);
frames = fix(framestot/(ncomps+1));

if ncomps < 1 || frames*(ncomps+1) ~= framestot,
    fprintf(...
  'chanproj(): data length (%d) not a multiple of ncomps (%d).\n',...
                                   framestot,ncomps);
    return
end

if chan> chans,
    fprintf(...
 'chanproj(): specified channel (%d) cannot be > than nchans (%d).\n',...
                            chan,chans);
    return
else
    chandata = projdata(chan,:);
    epochs = fix(length(chandata)/frames);
end
if epochs > MAXPLOTDATACHANS
  fprintf(...
 'chanproj(): maximum number of traces to plot is %d\n',...
              MAXPLOTDATACHANS);
   return
end
if framestot~=epochs*frames,
   fprintf('chanproj(): projdata is wrong length.\n');  % exit
   return
end

%
%%%%%%%%%%%%%%%% Find or read plot limits %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if limits==0,
  xmin=0;xmax=0;ymin=0;ymax=0;
else
  if length(limits)~=4,
      fprintf( ...
'chanproj():^G limits should be 0 or an array [xmin xmax ymin ymax].\n');
      return
  end
  xmin = limits(1);
  xmax = limits(2);
  ymin = limits(3);
  ymax = limits(4);
end

if xmin == 0 && xmax == 0,
  x = [0:frames-1];
  xmin = 0;
  xmax = frames-1;
else
  dx = (xmax-xmin)/(frames-1);
  x=xmin*ones(1,frames)+dx*(0:frames-1);          % construct x-values
end

if ymax == 0 && ymin == 0,
  ymax=max(max(projdata(chan,:)));
  ymin=min(min(projdata(chan,:)));
end
%
%%%%% Reshape the projdata for plotting and returning to user %%%%%%%%%
%
chandata=reshape(chandata,frames,epochs);
chandata=chandata';

if framelist ==0,
  framelist = [1:frames];  % default
end
if size(framelist,2)==1,
     if size(framelist,1)>1,
         framelist = framelist';
     else
         fprintf(...
 'chanproj(): framelist should be a vector of frames to plot per epoch.\n');
         return
     end
end

if framelist(1)<1 || framelist(length(framelist))>frames,
    fprintf('chanproj(): framelist values must be between 1-%d.\n',frames);
    return
else
    chandata=chandata(:,framelist);
end
%
%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if colorfile ~=0,
    if ~ischar(colorfile)
       fprintf('chanproj(); color file name must be a string.\n');
       return
    end
    cid = fopen(colorfile,'r');
    if cid <3,
        fprintf('chanproj(): cannot open file %s.\n',colorfile);
        return
    end
    colors = fscanf(cid,'%s',[3 MAXPLOTDATACHANS]);
    colors = colors';
    [r c] = size(colors);
    for i=1:r
        for j=1:c
            if colors(i,j)=='.',
                colors(i,j)=' ';
            end
        end
    end
else    % default color order - no yellow
    colors =['w  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  '];
end
%
% Make vector of x-values
%
x=[xmin:(xmax-xmin)/(frames-1):xmax+0.00001];
xmin=x(framelist(1));
xmax=x(framelist(length(framelist)));
x = x(framelist(1):framelist(length(framelist)));
%
%%%%%%%%% Compute percentage of variance accounted for %%%%%%%%%%%%%%
%
if epochs>1,
    sumdata = zeros(1,length(framelist));
    for e=2:epochs
        sumdata= sumdata + chandata(e,:);  % sum the component projections
    end
    sigvar = diag(cov(chandata(1,:)'));
    difvar = diag(cov(((chandata(1,:)-sumdata)')));
                                    % percent variance accounted for
    pvaf   = round(100.0*(1.0-difvar/sigvar));
end
%
%%%%%%%%%%%%%%%%%%%%%% Plot traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
fprintf('chanproj(): Drawing trace ');
for e=1:epochs,
    fprintf ('%d ',e);
    set(gcf,'Color',BACKCOLOR); % set the background color to grey
    set(gca,'Color','none');    % set the axis color = figure color
    if e==1
      plot(x,chandata(e,:),colors(e),'LineWidth',2);  % plot it!
    else
      plot(x,chandata(e,:),colors(e),'LineWidth',1);  % plot it!
    end
    hold on;
end
fprintf('\n');
%
%%%%%%%%%%%%%%%%%%%% Fix axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
axis([xmin xmax ymin-0.1*(ymax-ymin) ymax+0.1*(ymax-ymin)]);
%
%%%%%%%%%%%%%%%%%%% Add title and labels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if titl==0,
    titl = ['channel ' int2str(chan)];
end
titl = [ titl '   (p.v.a.f. ' int2str(pvaf) '%)' ];
title(titl);
xlabel('Time (msec)');
ylabel('Potential (uV)');

author: sccn

repo: eeglab


text:
function ekg_filt(up)
%EKG_FILT extracts respiratory signals using various filtering techniques
% from the ECG signal as specified in PC's literature review.
%	            ekg_filt(option, up)
%
%	Inputs:
%		option          the option which has led to this function being used
%       up              universal parameters structure
%
%	Outputs:
%       ...
%

fprintf('\n--- Extracting Respiratory Signals using Filtering Techniques ');
log_int_respSig = 0;             % Has value 1 unless this is a final respiratory signal

for subj = up.paramSet.subj_list

    %% Cycle through each method
    for filt_no = 1 : length(up.al.options.ekg_filt)

        %% Skip if this processing has been done previously
        eval(['save_name = ''ekg', up.paths.filenames.filt '_' up.al.options.ekg_filt{filt_no} ''';']);
        iden_resp_sig_file_ending
        savepath = [up.paths.data_save_folder, num2str(subj), ending];
        exist_log = check_exists(savepath, save_name);
        if exist_log
            continue
        end

        %% Load relevant data
        if ~exist('data', 'var')
            load([up.paths.data_load_folder, up.paths.data_load_filename]);
        end
        % Extract EKG data
        rel_data = data(subj).ekg;
        rel_data.fs = data(subj).ekg.fs;

        %% Filter the raw signal using this method
        respWave = feval(up.al.options.ekg_filt{filt_no}, rel_data, up);

        %% Band-pass filter
        filtered_data = bpf_signal_to_remove_non_resp_freqs(respWave, respWave.fs, up);
        eval([save_name ' = filtered_data;']);

        %% Save processed data
        save_or_append_data
    end

end

end

author: peterhcharlton

repo: RRest


text:
function value = bvec_enum ( n )

%*****************************************************************************80
%
%% BVEC_ENUM enumerates the binary vectors of length N.
%
%  Discussion:
%
%    A BVEC is a vector of binary digits representing an integer.
%
%    BVEC(1) is 0 for positive values and 1 for negative values, which
%    are stored in 2's complement form.
%
%    For positive values, BVEC(N) contains the units digit, BVEC(N-1)
%    the coefficient of 2, BVEC(N-2) the coefficient of 4 and so on,
%    so that printing the digits in order gives the binary form of the number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the length of the vectors.
%
%    Output, integer BVEC_ENUM, the number of binary vectors.
%
  value = 2 ^ n;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gauss_seidel_poisson_1d ( k )

%*****************************************************************************80
%
%% GAUSS_SEIDEL_POISSON_1D uses Gauss-Seidel iteration for 1D Poisson equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Commandline input, integer K, the grid index.
%    K specifies the number of nodes, by the formula NK = 2^K + 1.
%    K must be at least 0.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSS_SEIDEL_POISSON_1D:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use Gauss-Seidel iteration for the 1D Poisson equation.\n' );

  if ( nargin < 1 )
    k = input ( '  Enter K >= 0, the grid index, so that NK = 2^K + 1:  ' );
  elseif ( ischar ( k ) )
    k = str2num ( k );
  end
%
%  Set boundaries.
%
  a = 0.0;
  b = 1.0;
%
%  Set boundary conditions.
%
  ua = 0.0;
  ub = 0.0;
%
%  Get NK.
%
  nk = 2^k + 1;
%
%  Set XK.
%
  xk = ( linspace ( a, b, nk ) )';
%
%  Get HK.
%
  hk = ( b - a ) / ( nk - 1 );
%
%  Set FK.
%
  fk = force ( xk );
  fk(1) = ua;
  fk(nk) = ub;
%
%  Set the -1/2/-1 entries of Ak.
%
%  In order that the operator Ak approximate the Poisson operator,
%  and in order that we can compare linear systems for successive grids,
%  we should NOT multiply through by hk^2.
%
  Uk = sparse ( 2:nk-1, 3:nk,   -1.0, nk, nk ) / hk^2;
  Dk = sparse ( 2:nk-1, 2:nk-1,  2.0, nk, nk ) / hk^2;
  Dk(1,1) = 1.0;
  Dk(nk,nk) = 1.0;
  Lk = sparse ( 2:nk-1, 1:nk-2, -1.0, nk, nk ) / hk^2;
  Ak = Uk + Dk + Lk;
%
%  Just because we can, ask MATLAB to get the exact solution of the linear system
%  directly.
%
  udk = Ak \ fk;
%
%  Sample the solution to the continuous problem.
%
  uek = exact ( xk );
%
%  Use Gauss-Seidel iteration to solve the linear system to the given tolerance.
%
  ujk = zeros ( nk, 1 );
  it_max = 10000000;
  res_max = 0.000001;

  [ ujk, it_num, res_norm ] = gauss_seidel ( nk, Lk, Dk, Uk, fk, ujk, it_max, ...
    res_max );
%
%  Examine errors:
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using grid index K = %d\n', k );
  fprintf ( 1, '  System size NK was %d\n', nk );
  fprintf ( 1, '  Tolerance for linear residual %g\n', res_max );
  fprintf ( 1, '  RMS of linear residual %g\n', res_norm );
  fprintf ( 1, '  Number of GS iterations taken was %d\n', it_num );
  fprintf ( 1, '  RMS GS error in solution of linear system = %g\n', ...
    norm ( udk - ujk ) / sqrt ( nk ) );
  fprintf ( 1, '  RMS discretization error in Poisson solution = %g\n', ...
    norm ( uek - ujk ) / sqrt ( nk ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I        X        U_Exact    U_Direct    U_GS\n' );
  fprintf ( 1, '\n' );
  for i = 1 : nk
    fprintf ( 1, '  %4d  %10.4f  %10.4g  %10.4g  %10.4g\n', ...
      i, xk(i), uek(i), udk(i), ujk(i) );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSS_SEIDEL_POISSON_1D:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function uex = exact ( x )

%*****************************************************************************80
%
%% UEX evaluates the solution of the continuous problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the evaluation points.
%
%    Output, real UEX(*), the solution of the continuous problem
%    at the evaluation points.
%
  uex = x .* ( x - 1 ) .* exp ( x );

  return
end
function f = force ( x )

%*****************************************************************************80
%
%% FORCE evaluates the forcing term.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the evaluation points.
%
%    Output, real F(*), the forcing term at the evaluation points.
%
  f = - x .* ( x + 3 ) .* exp ( x );

  return
end
function [ u, it_num, res_norm ] = gauss_seidel ( n, Lk, Dk, Uk, f, u, ...
  it_max, res_max )

%*****************************************************************************80
%
%% GAUSS_SEIDEL carries out the Gauss-Seidel iteration.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real sparse Lk(*,*), Dk(*,*), Uk(*,*), the strict lower triangle,
%    the diagonal, and the strict upper triangle of A.
%
%    Input, real F(N), the right hand side.
%
%    Input, real U(N), the estimated solution.
%
%    Input, integer IT_MAX, the maximum number of iterations.
%
%    Input, real RES_MAX, a convergence tolerance on the residual.
%
%    Output, real U(N), the improved estimate of the solution.
%
%    Output, integer IT_NUM, the number of iterations taken.
%
%    Output, real RES_NORM, the RMS norm of the residual.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Step    Residual      Change\n' );
  fprintf ( 1, '\n' );

  it_num = it_max;

  for it = 1 : it_max

    u_old = u;
    u = ( Lk + Dk ) \ ( f - Uk * u_old );
    r = ( Lk + Dk + Uk ) * u - f;
    res_norm = norm ( r ) / sqrt ( n );

    fprintf ( 1, '  %6d  %10.4g  %10.4g\n', ...
      it, res_norm, norm ( u - u_old ) / sqrt ( n ) );

    if ( res_norm <= res_max )
      it_num = it;
      break;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function [theta] = normalEqn(X, y)
%NORMALEQN Computes the closed-form solution to linear regression
%   NORMALEQN(X,y) computes the closed-form solution to linear
%   regression using the normal equations.

theta = zeros(size(X, 2), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the code to compute the closed form solution
%               to linear regression and put the result in theta.
%

% ---------------------- Sample Solution ----------------------

theta = pinv(X' * X) * X' * y;

% -------------------------------------------------------------

% ============================================================

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
%%
%   Author: Yash Bansod
%
% GitHub: <https://github.com/YashBansod>
%
% This is the main program.

%% Clear the environment and the command line
clear;
clc;
close all;

%% Add the directory containing relevant functions to the path variables
addpath('./INV-functions/')

%% Define the input parameters and simulate

% Set the length of the links of the manipulator robot.
L1 = 5;
L2 = 5;
L3 = 5;

% Set the initial orientation of the robot.
theta1 = 10;
theta2 = 0;
theta3 = 15;

% Define the radius of the circle the end effector should follow
radius = 10;
r_sq = radius ^ 2;

hold on;
% Code for drawing a circle
for i = -radius: radius/10: radius
    expX = i;
    expY = sqrt(r_sq - expX^2);
    cla;
    images.roi.Circle(gca,'Center',[0 0],'Radius',radius, 'Facealpha', 0.05);
    % You can modify the Kp value defined in PLANAR_INV_KIN_3DOF to modify
    % the inverse jacobian controller behavior.
    [expPoint, Joint, Theta] = PLANAR_INV_KIN_3DOF(L1, L2, L3, expX, ...
                                            expY, theta1, theta2, theta3);
    scatter(Joint(end,1), Joint(end,2));
    theta1 = Theta(1, 1);
    theta2 = Theta(2, 1);
    theta3 = Theta(3, 1);
end

for i = radius: -radius/10: -radius
    expX = i;
    expY = -sqrt(r_sq - expX^2);
    cla;
    images.roi.Circle(gca,'Center',[0 0],'Radius',radius, 'Facealpha', 0.05);
    % You can modify the Kp value defined in PLANAR_INV_KIN_3DOF to modify
    % the inverse jacobian controller behavior.
    [expPoint, Joint, Theta] = PLANAR_INV_KIN_3DOF(L1, L2, L3, expX, ...
                                            expY, theta1, theta2, theta3);
    theta1 = Theta(1, 1);
    theta2 = Theta(2, 1);
    theta3 = Theta(3, 1);
end

author: YashBansod

repo: Robotics-Planning-Dynamics-and-Control


text:
% FT_PREAMBLE_INIT is a helper script that is used at the start of all FieldTrip main
% functions. It checks whether the user specified at lease one input arguments (i.e.
% the cfg) or shows the help of the calling function. It merges the global defaults
% with the cfg. It checks whether the output file already exists and whether it is OK
% to overwrite it. It tracks the function call.
%
% Use as
%   ft_preamble init
%
% See also FT_PREAMBLE, FT_POSTAMBLE, FT_TRACKUSAGE

% Copyright (C) 2011-2016, Robert Oostenveld, DCCN
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% disabled for now, see further down
global ft_default

% Which high-level FieldTrip function is lowest on the stack, and which is
% highest? Lowest could e.g. be ft_selectdata called by another FieldTrip
% function, highest should always be the top-level FieldTrip function that
% was explicitly called by the user. These are used further down below,
% e.g. in the reproducescript functionality.
% (Note that we prepend random strings to local variables in pre- and
% postambles to avoid contaminating the caller's namespace.)
FjmoT6aA_stack = dbstack('-completenames');
% FjmoT6aA_stack(1) is this script
% FjmoT6aA_stack(2) is the calling ft_postamble function
% FjmoT6aA_stack(3) is the lowest FieldTrip function on the stack, the one currently being evaluated
% the highest FieldTrip function on the stack is given by the highest n
% such that FjmoT6aA_stack(n).path begins with the FieldTrip path
FjmoT6aA_lowest_ft = FjmoT6aA_stack(3).name;
FjmoT6aA_highest_ft = FjmoT6aA_stack(3).name;
[FjmoT6aA_ft_ver, FjmoT6aA_ft_path] = ft_version;
for FjmoT6aA_k = 3:numel(FjmoT6aA_stack)
  if startsWith(FjmoT6aA_stack(FjmoT6aA_k).file, FjmoT6aA_ft_path)
    FjmoT6aA_highest_ft = FjmoT6aA_stack(FjmoT6aA_k).name;
  else
    % we are operating under the assumption here that a private FieldTrip
    % function (i.e. one whose name does not start with "ft_") never calls
    % another high-level FieldTrip function.
    break;
  end
end
% check whether we are currently in a top-level user-called FT function
% (note that checking of names only probably is not sufficient, as recursion
% could occur)
FjmoT6aA_current_ft_toplevel = strcmp(FjmoT6aA_lowest_ft, FjmoT6aA_highest_ft) && ...
  (FjmoT6aA_k == 4 || numel(FjmoT6aA_stack) == 3);

if ft_nargin==0
  help(FjmoT6aA_lowest_ft);
  % throw the error as if it happened in the original function
  msg.message     = 'This function requires one or multiple input arguments, please refer to the documentation above';
  msg.identifier  = '';
  msg.stack       = FjmoT6aA_stack;
  error(msg);
end % if nargin

% check if there are fieldnames in the cfg that suggest as if the user
% erroneously inputted a data argument
checkdatafields = isfield(cfg, {'cfg' 'label' 'dimord' 'trialinfo' 'avg' 'powspctrm'});
if any(checkdatafields)
  help(FjmoT6aA_lowest_ft);
  % throw the error as if it happened in the original function
  msg.message     = 'It seems as if the first input argument is a FieldTrip data structure, while a cfg is expected';
  msg.identifier  = '';
  msg.stack       = FjmoT6aA_stack;
  error(msg);
end

% convert automatically from cell-array to structure
if iscell(cfg)
  cfg = ft_keyval2cfg(cfg);
end

% check that it is an struct or empty numeric array
assert(isstruct(cfg) || (isnumeric(cfg) && isempty(cfg)), 'The configuration must be a structure or empty');

% this script requires some options that can be user-specified, but otherwise are obtained from ft_default
% merge the default options into the configuration, except the preamble field which is used for passing arguments
cfg = mergestruct(cfg, ft_default);

% determine whether function execution should be aborted or continued
if isfield(cfg, 'outputfile') && ~isempty(cfg.outputfile)
  assert(any(strcmp(fieldnames(cfg), 'outputfilepresent')), 'cfg.outputfilepresent is a required option, please see FT_DEFAULTS');
  % check whether the output file already exists
  if ischar(cfg.outputfile)
    chiL7fee_outputfile = {cfg.outputfile};
  else
    chiL7fee_outputfile = cfg.outputfile;
  end
  for i=1:numel(chiL7fee_outputfile)
    [p, f, x] = fileparts(chiL7fee_outputfile{i});
    if isempty(p)
      % the relative path was speciield
      chiL7fee_outputfile{i} = fullfile(pwd, chiL7fee_outputfile{i});
    else
      % the absolute path was specified
      chiL7fee_outputfile{i} = chiL7fee_outputfile{i};
    end
    if ~exist(chiL7fee_outputfile{i}, 'file')
      ft_abort = false;
    else
      % the output file exists, determine how to deal with it
      switch cfg.outputfilepresent
        case 'keep'
          if ft_nargout>0
            % continue executing the parent function
            ft_warning('output file %s is already present, but you also requested an output argument: continuing function execution', chiL7fee_outputfile{i});
            ft_abort = false;
          else
            % stop executing the parent function
            ft_warning('output file %s is already present: aborting function execution', chiL7fee_outputfile{i});
            ft_abort = true;
          end
        case 'overwrite'
          ft_warning('output file %s is already present: it will be overwritten', chiL7fee_outputfile{i});
          ft_abort = false;
        case 'error'
          ft_error('output file %s is already present', chiL7fee_outputfile{i});
        otherwise
          ft_error('invalid option for cfg.outputfilepresent');
      end % case
    end
  end % for each of the output files
else
  % there is no reason to abort execution
  ft_abort = false;
end % if chiL7fee_outputfile{i}

if isfield(cfg, 'reproducescript') && ~isempty(cfg.reproducescript)
  % the reproducescript code should only be executed in a top-level FT function
  if ~FjmoT6aA_current_ft_toplevel
    % we are in a FT function that was called by another FT function
    cfg = rmfield(cfg, 'reproducescript');
  else
    % we are in a top-level FT function
    if ~isfolder(cfg.reproducescript)
      mkdir(cfg.reproducescript);
    end
    % user-specified cfg.inputfile or cfg.outputfile are not compatible
    % with cfg.reproducescript functionality, throw error to make them
    % mutually exclusive
    if (isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)) || ...
       (isfield(cfg, 'outputfile') && ~isempty(cfg.outputfile))
      ft_error('cfg.reproducescript cannot be used together with a user-specified cfg.inputfile or cfg.outputfile');
    end
    % this variable is used in loadvar, savevar and savefig
    Fief7bee_reproducescript = cfg.reproducescript;
    cfg = rmfield(cfg, 'reproducescript');
    % pause one second to ensure that subsequent file names (which contain the time stamp) are unique
    pause(1);
  end
end

if false
  % this is currently generating too much data and therefore disabled
  if isfield(cfg, 'trackusage') && ~(isequal(cfg.trackusage, false) || isequal(cfg.trackusage, 'no') || isequal(cfg.trackusage, 'off'))
    % track the usage of the calling function
    stack = dbstack('-completenames');
    % stack(1) is this script
    % stack(2) is the calling ft_postamble function
    % stack(3) is the main FieldTrip function that we are interested in
    ft_trackusage('function call', 'function', stack(3).name);
  end % if trackusage
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef TestGetStructuringElement
    %TestGetStructuringElement

    methods (Static)
        function test_1
            elem = cv.getStructuringElement();
        end

        function test_2
            shapes = {'Rect', 'Cross', 'Ellipse'};
            for i=1:numel(shapes)
                elem = cv.getStructuringElement(...
                    'Shape',shapes{i}, 'KSize',[3 3]);
                validateattributes(elem, {'numeric'}, {'size',[3 3]});
            end
        end

        function test_error_unrecognized_option
            try
                cv.getStructuringElement('foo','bar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% MAIN.m  --  Lesson 1 -- Simple Plot
%
% This script performs a simulation of the cart as it passively moves from
% some initial state.
%
% The purpose of this lesson is to learn how use ode45 to generate data for
% a simulation, and then plot the state of the system to see what is going
% on.
%

clc; clear;

%%%% Initial State
z0 = [
    0.0;   %horizontal position
    (pi/180)*80;  %pendulum angle (wrt gravity)
    0.3;   %horizontal velocity
    0.5];  %pendulum angular rate

%%%% Physical Parameters
p.m1 = 1.0;  % (kg) Cart mass
p.m2 = 0.3;  % (kg) pole mass
p.g = 9.81;  % (m/s^2) gravity
p.l = 0.5;   % (m) pendulum (pole) length

%%%% Time vector
t = linspace(0,2,250);  %Simulation time stamps

%%%% Function Handle
dynFun = @(t,z)( cartPoleDynamics(z, p) );

%%%% Simulate the system!
options = odeset(...
    'RelTol',1e-8, ...
    'AbsTol',1e-8);
[~, z] = ode45(dynFun, t, z0, options);   %  <-- This is the key line!
z = z';

%%%% Unpack the state:
x = z(1,:);
q = z(2,:);
dx = z(3,:);
dq = z(4,:);

%%%% Plots:
figure(1); clf;

subplot(2,2,1);
plot(t,x)
ylabel('x')
title('Position')

subplot(2,2,2);
plot(t,q)
ylabel('q')
title('Angle')

subplot(2,2,3);
plot(t,dx)
ylabel('dx')
title('Velocity')

subplot(2,2,4);
plot(t,dq)
ylabel('dq')
title('Angle Rate')

author: MatthewPeterKelly

repo: dscTutorials


text:
function Y=diff(varargin)
%DIFF (overloaded)

X = varargin{1};
n = X.dim(1);
m = X.dim(2);

Y = X;
x_lmi_variables = X.lmi_variables;
lmi_variables = [];

% Fast diff of large matrices. Convert to optmized case
switch nargin
    case 1
        if (min([n m]) > 1)
            Y = diff(X',1,2)';
            return
        end

    case 2
        if varargin{2} == 1 & (min([n m]) > 1)
            Y = diff(X',1,2)';
            return
        end
    case 3
        if varargin{2}==1 & varargin{3} == 2 & (min([n m]) > 1)
            shift = [-speye(m-1) spalloc(m-1,1,0)] + [spalloc(m-1,1,0) speye(m-1)];
            %            shift = shift(1:m-1,1:m);
            shift = kron(shift,speye(n));
            Y.basis = shift*X.basis;
            Y.dim(1) = n;
            Y.dim(2) = m - 1;
            % Reset info about conic terms
            Y.conicinfo = [0 0];
            Y = clean(Y);
            return
        elseif varargin{2}==1 & varargin{3} == 1  & (min([n m]) > 1)
            Y = diff(X',1,2)';
            return
        end

    otherwise
        error('To many input arguments.')
end

% Slow but safe case. Used for higher order diff
try
    j = 1;
    diffX = diff(reshape(X.basis(:,1),n,m),varargin{2:end});
    Y.basis=diffX(:);
    for i = 1:length(x_lmi_variables)
        diffX = diff(reshape(X.basis(:,i+1),n,m),varargin{2:end});
        if (norm(diffX,inf)>0)
            Y.basis(:,j+1) = diffX(:);
            lmi_variables = [lmi_variables x_lmi_variables(i)];
            j = j+1;
        end
    end
    if j~=1
        Y.lmi_variables = lmi_variables;
    else
        Y = full(reshape(Y.basis(:,1),size(diffX,1),size(diffX,2)));
        return
    end
catch
    error(lasterr)
end
Y.dim(1) = size(diffX,1);
Y.dim(2) = size(diffX,2);
% Reset info about conic terms
Y.conicinfo = [0 0];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [out] = read_tmsi_poly5(filename, varargin)

% READ_TMSI_POLY5
%
% Use as
%   hdr = read_tmri_poly5(filename)
%   dat = read_tmsi_poly5(filename, hdr, begblock, endblock)
%
% This implementation is as closely as possible based on the original "tms_read",
% which contains the comments
%
%   Changed on 08-10-2014 by TL, TMSi
%   - Now supports loading a file from different directory than the script file
%   - Feedback on the validity of arguments and whether a file could be found or not.
%   - Dialogue is opened when no argument was given.
%
%   Changed on 18-10-2022 by JMS, DCCN
%   - Massive speed up: no intermediate double->single->double conversion ,and
%   - Don't store metadata that is not broadcasted to outside function in a struct array, and
%   - Allow for a selection of channels to be read, reducing memory footprint, and calibration step
%
% 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.

needhdr = (nargin==1);
needdat = (nargin>1);

if needdat
  signal = varargin{1};
  begblock = varargin{2};
  endblock = varargin{3};
end

if numel(varargin)==4
  chanindx = varargin{4};
else
  chanindx = [];
end

fid = fopen_or_error(filename, 'rb');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if needhdr
  pos = 31;
  fseek(fid,pos,'bof');
  version = fread(fid,1,'int16');
  if version == 203
    frewind(fid);
    signal.header.FID                   = fread(fid,31,'uchar');
    signal.header.VersionNumber         = fread(fid, 1,'int16');
  else % version 204
    frewind(fid);
    signal.header.FID                   = fread(fid,32,'uchar');
    signal.header.VersionNumber         = fread(fid, 1,'int16');
  end
  signal.header.MeasurementName       = fread(fid,81,'uchar');
  signal.header.FS                    = fread(fid, 1,'int16');
  signal.header.StorageRate           = fread(fid, 1,'int16');
  signal.header.StorageType           = fread(fid, 1,'uchar');
  signal.header.NumberOfSignals       = fread(fid, 1,'int16');
  signal.header.NumberOfSamplePeriods = fread(fid, 1,'int32');
  signal.header.EMPTYBYTES            = fread(fid, 4,'uchar');
  signal.header.StartMeasurement      = fread(fid,14,'uchar');
  signal.header.NumberSampleBlocks    = fread(fid, 1,'int32');
  signal.header.SamplePeriodsPerBlock = fread(fid, 1,'uint16');
  signal.header.SizeSignalDataBlock   = fread(fid, 1,'uint16');
  signal.header.DeltaCompressionFlag  = fread(fid, 1,'int16');
  signal.header.TrailingZeros         = fread(fid,64,'uchar');

  %conversion to char of text values
  signal.header.FID               = char(signal.header.FID);
  signal.header.MeasurementName   = char(signal.header.MeasurementName(2:signal.header.MeasurementName(1)+1));
  %    signal.header.StartMeasurement  = typecast(int8(signal.header.StartMeasurement),'int16'); % commented by RS 1feb11

  signal.fs = signal.header.FS;

  %check for right fileversion
  % if signal.header.VersionNumber ~= 203
  %    disp('Wrong file version! Imput file must be a Poly5/TMS version 2.03 file!');
  %    return;
  % end

  % Signal description
  for g=1:signal.header.NumberOfSignals
    signal.description(g).SignalName        = fread(fid,41,'uchar');
    signal.description(g).Reserved          = fread(fid, 4,'uchar');
    signal.description(g).UnitName          = fread(fid,11,'uchar');
    signal.description(g).UnitLow           = fread(fid, 1,'float32');
    signal.description(g).UnitHigh          = fread(fid, 1,'float32');
    signal.description(g).ADCLow            = fread(fid, 1,'float32');
    signal.description(g).ADCHigh           = fread(fid, 1,'float32');
    signal.description(g).IndexSignalList   = fread(fid, 1,'int16');
    signal.description(g).CacheOffset       = fread(fid, 1,'int16');
    signal.description(g).Reserved2         = fread(fid,60,'uchar');

    % conversion of char values (to right format)
    signal.description(g).SignalName = char(signal.description(g).SignalName(2:signal.description(g).SignalName(1)+1));
    signal.description(g).UnitName   = char(signal.description(g).UnitName(2:signal.description(g).UnitName(1)+1));
  end

  % return the header details
  out = signal;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%read data blocks
NB = signal.header.NumberSampleBlocks;
SD = signal.header.SizeSignalDataBlock;
NS = signal.header.NumberOfSignals;
NS_32bit = NS/2;

if isempty(chanindx)
  chanindx = 1:NS_32bit;
end

if needdat
  nblock = endblock + 1 - begblock;
  nsmp = SD/(NS_32bit*4); % number of samples per data block
  on  = 1;
  off = nsmp;
  dat = nan(numel(chanindx), nsmp*nblock);
  PI  = nan(1,NB);
  BT  = nan(7,NB);
  for g=begblock:endblock

    % jump to right position in file
    if signal.header.VersionNumber == 203
      pos = 217 + NS*136 + (g-1) *(86+SD);
    else
      pos = 218 + NS*136 + (g-1) *(86+SD);
    end
    fseek(fid,pos,'bof');

    PI(g) = fread(fid,1,'int32'); % if no struct array is created, this saves a lot of time, data is not broadcasted anyway
    fread(fid,4,'uchar'); %reserved for extension of previous field to 8 bytes
    BT(:,g) = fread(fid,14/2,'int16'); %dostime
    fread(fid,64,'uchar'); %reserved
    data = fread(fid,SD/4,'float32'); % no reason to convert to single, and back to double downstairs

    % Convert data to 32bit values.
    % In case also 16bit values have to be measured, these values are typecasted below:
    % data = fread(fid,SD/2,'int16'); %read data
    % data = int16(data);
    % data = typecast(data,'int32');
    % signal.block(g).DATA = data;
    % signal.data{g} = reshape(data,NS_32bit,SD/(NS_32bit*4));
    tmp = reshape(data, [NS_32bit, nsmp]);
    dat(:, on:off) = tmp(chanindx, :);
    on  = off+1;
    off = off+nsmp;
  end % for

  % Converting data to a usable format
  ADCLow   = [signal.description(1:2:end).ADCLow];   ADCLow   = ADCLow(chanindx);
  ADCHigh  = [signal.description(1:2:end).ADCHigh];  ADCHigh  = ADCHigh(chanindx);
  UnitLow  = [signal.description(1:2:end).UnitLow];  UnitLow  = UnitLow(chanindx);
  UnitHigh = [signal.description(1:2:end).UnitHigh]; UnitHigh = UnitHigh(chanindx);
  for g = 1:size(dat,1)  % represent data in [uV]
    dat(g,:) = (dat(g,:) - ADCLow(g))./(ADCHigh(g) - ADCLow(g))  .* (UnitHigh(g) - UnitLow(g)) + UnitLow(g) ; %correction for uV
  end

  % return the data
  out = dat;
end

fclose(fid);


text:
function [SF,SFI,SFJ] = remove_duplicate_simplices(F)
  % REMOVE_DUPLICATE_SIMPLICES Remove duplicate simplices (regardless of order)
  %
  % Inputs:
  %   F  #F by ss list of simplex indices
  % Outputs:
  %   SF  #SF by ss list of unique simplices (maintaining original order of
  %     first simplex found)
  %   SFI  #SF list so that SF = F(SFI,:)
  %   SFJ  #SFI list so that F = SF(SFJ,:)
  [~,SFI,SFJ] = unique(sort(F,2),'rows','stable');
  SF = F(SFI,:);
end

author: alecjacobson

repo: gptoolbox


text:
%rate_c = converse_simo(nn,P,error,rx,K) computes the converse bound on the
%maximal channel coding rate over a quasi-static single-input
%multiple-output (SIMO) Rician fading channel with full channel state information.
%
%The function takes the following inputs:
%
%nn: blocklength, scalar/vector;
%P: input power (in linear scale), scalar;
%error: target block error probability, scalar, should be greater than 10^(-6);
%rx: number of receive antennas, scalar;
%K: rician K-factor, scalar; the default value is 0
function rate_c = converse_simo(nn,P,error,rx,K)

if (nargin < 5) || isempty(K)
	K = 0;
end

loop = min(1000/error, 10^7);    % number of samples of channel gain for Monte Carlo

%Generate samples of channel gain for Monte Carlo simulation
%channel gain is G = |H_1|^2 + ... +|H_{rx}|^2
if rx ==1
    G= abs(sqrt(K/(K+1))+sqrt(1/(K+1)/2)*(randn(rx,loop)+1i*randn(rx,loop))).^2 ;
else
    G = sum(abs(sqrt(K/(K+1))+sqrt(1/(K+1)/2)*(randn(rx,loop)+1i*randn(rx,loop))).^2 );
end

rate_c=[];

loop21=1000; %number of samples in numerical integration, set larger value for better precision
loop22=1000;

%samples of channel gain for numerical computation
%Todo: optimize the way of sampling G from ncx2 (noncentral chi square)
g_mid = ncx2inv(error*2, 2*rx,2*rx*K)/(2*K+2);
g_max = ncx2inv(1-10^(-5),2*rx,2*rx*K )/(2*K+2);

G1 = (g_mid/loop21):(g_mid/loop21):g_mid;
G2 = (g_mid+ (g_max -g_mid)/loop22):(g_max -g_mid)/loop22:g_max;
G_num=[G1,G2];

%compute the pdf of channel gain
pdf_G = ncx2pdf(G_num*(2*K+2), 2*rx, 2*rx*K)*(2*K+2);

%pdf_G = (K+1) * exp(sqrt(4*K*2*(K+1)*G)-(K+1).*G-2*K).* sqrt(2*(K+1).*G/4/K).*besseli(1,sqrt(4*K*2*(K+1)*G),1);

precision = 10000; %number of samples for numerical integration; set larger value to increase precision

for n = nn

    %compute gamma_n through Monte Carlo
    S_n = n.*log(1+P*G) + n - ncx2rnd(2*n, 2*n/P./G).*P.*G./(1+P*G)/2; %information density
    S_n_sort=sort(S_n);
    gamma_n = S_n_sort(error*loop); % Pr[S_n <= gamma_n]=error
    %%%Todo: Change to numerical computations

    % Calculate Q[log dP/dQ >=gamma_n]; Given G, log dP/dQ is ncx2 distributed
    y_a=zeros(1,length(G_num)); % conditional probability P(Ln\geq n\gamma_n) given g for different g values

    for ii=1:length(G_num)
        gg= G_num(ii);
        k_n = 2/P/gg*( n*log(1+P*gg) + n-gamma_n );
        noncentrality = 2* n*(1+P*gg)/P/gg;
        %use ncx2cdf to compute cdf of noncentral chi-square
        cond_cdf = ncx2cdf(k_n , 2*n, noncentrality);

        if cond_cdf ==0 && k_n>0
            %Numerically compute the cdf of ncx2; this is done by convolving chi^2 pdf and |Gaussian(noncentrality,1)|^2

            step = sqrt(k_n)/precision;

            t = (sqrt(noncentrality)-sqrt(k_n)) :step:(sqrt(noncentrality) + sqrt(k_n));

            y_a(ii) = sum( exp(-1*(t).^2/2 + log( gammainc((k_n -(sqrt(noncentrality) - t).^2 )/2,(2*n-1)/2) ) ) ) * step./sqrt(2*pi);
        else
            y_a(ii)=cond_cdf;
        end
    end

    beta = sum( y_a(1:length(G1)) .* pdf_G(1:length(G1)))*g_mid/loop21 + sum(y_a(length(G1)+1: end) .*pdf_G(length(G1)+1:end)) *(g_max- g_mid)/loop22; % average over G.

    rate=-1*log2(beta)/(n-1);
    rate_c=[rate_c,rate];
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
% FT is the enumeration class for the typel of fields.
classdef FT < Enumerated
	enumeration
		e('E')
		h('H')
	end

	methods (Static)
		function elems = elems(ind)
			elems = [FT.e, FT.h];
			if nargin > 0  % ind
				elems = elems(ind);
			end
		end

		function count = count()
			count = length(FT.elems);
		end
	end

	methods
		function f = alter(this)
			% This function simply returned FT.h if this == FT.e and FT.e if
			% this == FT.h originally, but it is modified to handle arguments
			% given as arrays.
			chkarg(istypeof(this, 'FT'), '"this" should have FT as elements.');
			if isempty(this)  % alter([]) makes no sense, but alter(FT.empty()) does
				f = [];
			else
				f = this.elems(FT.count + 1 - int(this));
			end
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function v = polynomialCurveDerivative(t, varargin)
%POLYNOMIALCURVEDERIVATIVE Compute derivative vector of a polynomial curve
%
%   VECT = polynomialCurveLength(T, XCOEF, YCOEF);
%   XCOEF and YCOEF are row vectors of coefficients, in the form:
%       [a0 a1 a2 ... an]
%   VECT is a 1x2 array containing direction of derivative of polynomial
%   curve, computed for position T. If T is a vector, VECT has as many rows
%   as the length of T.
%
%   VECT = polynomialCurveLength(T, COEFS);
%   COEFS is either a 2xN matrix (one row for the coefficients of each
%   coordinate), or a cell array.
%
%   Example
%   polynomialCurveDerivative
%
%   See also
%   polynomialCurves2d, polynomialCurveNormal, polynomialCurvePoint,
%   polynomialCurveCurvature
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-02-23
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Extract input parameters

% polynomial coefficients for each coordinate
var = varargin{1};
if iscell(var)
    xCoef = var{1};
    yCoef = var{2};
elseif size(var, 1)==1
    xCoef = varargin{1};
    yCoef = varargin{2};
else
    xCoef = var(1,:);
    yCoef = var(2,:);
end

%% compute derivative

% compute derivative of the polynomial
dx = polynomialDerivate(xCoef);
dy = polynomialDerivate(yCoef);

% convert to polyval convention
dx = dx(end:-1:1);
dy = dy(end:-1:1);

% numerical integration of the Jacobian of parametrized curve
v = [polyval(dx, t) polyval(dy, t)];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function varargout = rosenbruck(X)
% Extended Rosenbruck's Banana-function, for N-dimensional input
%
%   ROSENBRUCK([x1, x2, .., xn]) returns the value of the Rosenbruck
%   function at the specified points. All [xi] may be vectors. The search
%   domain is
%
%               -100 < x_i < 100
%
%   The global minimum is
%
%               f(x1, x2, ..., xn) = f(1, 1, ..., 1) = 0

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

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

    % otherwise, output function value
    else

        % keep all values within the domain
        X(X < -100) = inf;  X(X > 100) = inf;

        % split input vector X into X1, X2
        % NOTE: proper orientation can not be determined automatically
        % the sum is taken by default over the rows:
        X1 = X(1:2:end-1, :);        X2 = X(2:2:end, :);

        % output rowsum
        varargout{1} = sum(  100*(X2 - X1.^2).^2 + (1 - X1).^2, 1);
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function [splitPsi, logQ] = getSplitConfig_SeqAllocation( Psi, data_struct, hyperparams, model, anchorObjIDs, featIDs, TargetPsi )

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

if isfield( Psi, 'theta' );
    Psi = rmfield( Psi, 'theta');
    Psi = rmfield( Psi, 'TS');
end
nObj = size(Psi.F,1);
logQ = struct( 'F', 0, 'z', 0 );

obsModel = model.obsModel;

ii = anchorObjIDs(1);
jj = anchorObjIDs(2);

kold = featIDs(1);

otherFeatIDs = setdiff( 1:size(Psi.F,2), featIDs );

otherObjIDs = setdiff( find( Psi.F( :, kold ) > 0   )', anchorObjIDs );
allObjIDs = [anchorObjIDs otherObjIDs];

kA = size(Psi.F,2) + 1;
kB = size(Psi.F,2) + 2;
propFeatIDs = [kA kB];

propF = Psi.F;
propF(:, [kold] ) = 0;
propF(ii, [kA kB] ) = [1 0];
propF(jj, [kA kB] ) = [0 1];

% First, draw starting value for theta from sequence ii and jj
propStateSeq = Psi.stateSeq;
propStateSeq( ii ).z(  propStateSeq( ii ).z == kold )  = kA;
propStateSeq( jj ).z(  propStateSeq( jj ).z == kold )  = kB;

% Create deterministic trans params Eta
propTS = getEta_PriorMean( propF, [], hyperparams, model, allObjIDs, propFeatIDs );

% Create deterministic Theta for ALL features: others + [kA kB]
[Xstats] = getXSuffStats( propF, propStateSeq, data_struct, model, 1:nObj, [otherFeatIDs propFeatIDs] );
myTheta = setThetaToPosteriorMean( [], Xstats, obsModel, [otherFeatIDs propFeatIDs] );

if exist( 'TargetPsi', 'var' ) && ~isempty( TargetPsi )
    TargetPsi.externalFeatIDs = propFeatIDs;
end

% Keep *only* the params for the feature kA, kB being split
%   this is all we update
AXstats  = getXSuffStats( propF, propStateSeq, data_struct, model, [ii jj], [propFeatIDs] );

% Initialize some key sufficient stats
featCounts = sum( propF( anchorObjIDs, propFeatIDs ), 1 );
nObj = length( anchorObjIDs );
counter=0;

for aa = [ otherObjIDs( randperm(length(otherObjIDs)) ) ii jj]
    counter=counter+1;

    if aa == ii || aa == jj
        % Make sure to update F suff. stats for anchors
        %  since we've already "seen" them before, unlike other sequences
        featCounts = featCounts - propF(aa, propFeatIDs);
        nObj = nObj - 1;

        Xdelta = getXSuffStats( propF, propStateSeq, data_struct, model, aa, propFeatIDs );
        AXstats = decXStats( AXstats, Xdelta, model.obsModel );
    end

    % Calc Soft Evidence for the current sequence, including both new feats
    curFeatIDs = find( propF(aa,:));
    allPropFeatIDs = union( curFeatIDs, propFeatIDs );
    LL = calcLogCondPrObsGivenTheta( data_struct(aa), myTheta, obsModel, allPropFeatIDs, size(propF,2) );

    % Propose new feature assignments for the cur. seq.
    myFPr = featCounts ./ ( nObj + hyperparams.c );
    [newFaa, logQ_Faa] = sampleFfromRestrictedSet( aa, anchorObjIDs, LL, propFeatIDs, propTS.obj( aa ), myFPr, TargetPsi );
    propF(aa, propFeatIDs) = newFaa;

    % Keep track of probability of making this specific proposal
    logQ.F = logQ.F + logQ_Faa;

    % In case we chose not to include one of the proposed new split feats,
    %   we need to remove that feature's transition weights from propTS
    %   and also remove that feature's entries in soft ev. matrix LL
    availFeatIDs = propTS.obj(aa).availFeatIDs;
    kiloc = find( availFeatIDs == propFeatIDs(1) );
    kjloc = find( availFeatIDs == propFeatIDs(2) );
    switch ( newFaa * [10;1] )
        case 01
            keepIDs = [1:kiloc-1 kiloc+1:length(availFeatIDs)];
        case 10
            keepIDs = [1:kjloc-1 kjloc+1:length(availFeatIDs)];
        case 11
            keepIDs = 1:length(availFeatIDs);
    end
    curEta.availFeatIDs = propTS.obj(aa).availFeatIDs( keepIDs );
    curEta.pi_z = propTS.obj(aa).pi_z( keepIDs, keepIDs );
    activeFeatIDs = union( curFeatIDs, propFeatIDs(newFaa==1) );
    activeLL = LL( activeFeatIDs, : );

    % Propose new State Seq for current sequence
    %[propStateSeq, logQaa] = sampleZ_RGS(data_struct, propTS, propF, myTheta, obsModel, propStateSeq, aa, TargetPsi);
    [propStateSeq(aa).z, logPrQ_z] = sampleStateSeq_PrecompSoftEv( aa, curEta, activeLL, TargetPsi);
    logQ.z = logQ.z + logPrQ_z;

    % Update suff. stats for F sampling
    featCounts = featCounts + propF(aa, propFeatIDs);
    nObj = nObj + 1;

    % Update Emission Params for new proposed features
    %   using the newly allocated data from the current sequence
    Xdelta = getXSuffStats( propF, propStateSeq, data_struct, model, aa, propFeatIDs );
    AXstats = incXStats( AXstats, Xdelta, model.obsModel );
    myTheta = setThetaToPosteriorMean( myTheta, AXstats, obsModel, [propFeatIDs] );

end

% Record final launch state
splitPsi.F = propF;
splitPsi.stateSeq = propStateSeq;
splitPsi.activeFeatIDs = propFeatIDs;
splitPsi.anchorObjIDs  = anchorObjIDs;
end

author: michaelchughes

repo: NPBayesHMM


text:
function[str_all, Irec] = icassp12(example, Img, Nloops, SHOW_IMGS);

if nargin < 4
  SHOW_IMGS = false;
  if nargin < 3
    Nloops = 1;
    if nargin < 2
      Img = lena;
      if nargin < 1;
	example = 'gauss_sp_denoise';
      end
    end
  end
end

%  SHOW_IMGS = true;

%  example = 'icassp12'
%  example = 'gauss_sp_denoise';
%  example = 'gauss_random_denoise';

%  Img = 'Clena';
%  Img = 'peppers';
%  Img = 'lena';
%  Img = 'barb';
%  Img = 'cman';
%  Img = 'Cboats';

%  Img = 'all';

%  Nloops = 1; % This value will be change to 10 if example == icassp12

%%%%%%%%%%%%%%%%%%%%%%%%
%        Setup         %
%%%%%%%%%%%%%%%%%%%%%%%%

  % Check for SSIM
  SSIM_CODE = check4ssim();

  % Setup simulation
  [Nloops, Img] = setupSim(example, Nloops, Img);

  % open string (for reports)
  str_all = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Input images        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( strcmp(lower(Img), 'all') )

  all = 1;
  if( strcmp(example, 'icassp12') )
    Img = {'lena', 'goldhl', 'cman512', 'barb', 'Clena', 'Cbarb', 'Cboats', 'Cgoldhl' };
  else
    Img = {'lena', 'peppers', 'goldhl.', 'cman512', 'Clena', 'Cpeppers', 'Cmandrill'};
  end

  NImgs = length(Img);

else

  all = 0;
  NImgs = 1;

end  % _END_ IF(strcmp(lower(Img))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          loop over Img
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for t = 1:NImgs,

  if( all == 1 )
    tImg = Img{t};
  else
    tImg = lower(Img);
  end

  [Ig Color] = selectImage(tImg);

  if Color == 1
    SSIM_CODE = 0;
  end

  disp(' ');
  disp(sprintf('%s Color %d', tImg, Color));
  disp(' ');

  [Nrows Ncols Ndims] = size(Ig);

  % Regularization parameter of the L2 / L1
  [lambdaS_l2 Lambda_O Lambda_update] = setRegParameter(example, Ndims);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  loops for average results    %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for loops=1:Nloops

    %%%%%%%%%%%%%%%%%
    %  add noise    %
    %%%%%%%%%%%%%%%%%

    % Noise level
    sigmaGauss = [5/255, 10/255, 15/255 ];
    NGauss = length(sigmaGauss);

    spnoise = [0.3 0.5 0.7];
    Nspnoise = length(spnoise);

    randomnoise = [0.1 0.2 0.3];
    NRnoise = length(randomnoise);

    % noisy images
    clear InCell;
    [InCell noiseName vNoise nNoise] = addNoise(Ig, example, sigmaGauss, NGauss, spnoise, Nspnoise, randomnoise, NRnoise);
    nNoisyIm = length(InCell);

    % loop over noisy images
    for k=1:length(InCell)

      In = InCell{k};

%^^^^^% ^^^^^^^^^^^^^^^^^^^^^
%     %   Clock is ticking
%     % ---------------------
      tstart = tic;

      % First estimate the pixels that are corrupted
      switch lower(example)

        case{'gauss_sp_denoise'}
          S0 = adaptMedian_colfilt(In);		% estimate S&P corrupted pixels

        case{'gauss_random_denoise'}	% estimate random value (impulse noise) corrupted pixels
          S0 = estimateRValImpulse(In);

      end	% _END_ SWITCH

      % Set lambda (regularization) matrix
      lambda_l1l2 = zeros(Nrows, Ncols, 2*Ndims);

      for d=1:Ndims,
          lambda_l1l2(:,:,d)  =  Lambda_O*(S0(:,:,d) ~= 0);
          lambda_l1l2(:,:,Ndims+d) = lambdaS_l2(k)*(S0(:,:,d)==0);
      end

      % =================================================
      % ---------- Setup for IRN_L1L2  ------------------

      pars_irn_adapt = irntvInputPars('l1tv_l2tv');

      pars_irn_adapt.adapt_epsR   = 1;
      pars_irn_adapt.epsR_cutoff  = 0.01;
      pars_irn_adapt.adapt_epsF   = 1;
      pars_irn_adapt.epsF_cutoff  = 0.05;

      pars_irn_adapt.pcgtol_ini = 1e-4;

      switch lower(example)
        case{'gauss_sp_denoise'}
          pars_irn_adapt.loops      = 1;
          IRNloops = 7;
        case{'gauss_random_denoise'}
          pars_irn_adapt.loops      = 2;
          IRNloops = 4;
      end

      pars_irn_adapt.U0      = In;   % necessary the for adapt case

      % =================================================

      % 1st call
      Ig_L1L2 = irntv(In, {}, lambda_l1l2, pars_irn_adapt);

      % loop

      for m = 2:IRNloops

        for d=1:Ndims
          lambda_l1l2(:,:,d)  =  Lambda_update*lambda_l1l2(:,:,d);
        end

        pars_irn_adapt.U0      = Ig_L1L2;	% reset initial sol. to previous sol.

        Ig_L1L2 = irntv(In, {}, lambda_l1l2, pars_irn_adapt);

      end % _END_ FOR(m)

      % denoised image
      Irec{k} = Ig_L1L2;

%^^^^^% ^^^^^^^^^^^^^^^^^^^^^
%     %    Time performance
      l1l2_time(k, loops) = toc(tstart);
%     % ---------------------

      if(SHOW_IMGS)
        figure; imagesc( Normalize(Irec{k}) ); axis image; axis off; colormap gray;
      end

      % Reconstruction quality assessment
      [psnr_l1l2(k, loops), snr_l1l2(k, loops), ssim_l1l2(k, loops)] = qualityAssessment(Ig, Irec{k}, SSIM_CODE);

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      %           Generate reports

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      if(k==1)

        if(Color == 0) str = sprintf('Img: grayscale %s (denoising)\n', lower(tImg));
        else str = sprintf('Img: color %s (denoising)\n', lower(tImg));
        end

        str_all = [ str_all str ];

        str = sprintf(' sigma/%s \t SNR \t SSIM \t PSNR    Time\n', noiseName);
        str_all = [ str_all str ];

      end

        str = sprintf(' %1.0f/255 -- %1.1f \t %1.2f \t %1.2f \t %1.2f    %1.2f\n', ...
                   255*sigmaGauss( rem(k,NGauss) + NGauss*(rem(k,NGauss)==0) ), vNoise( floor((k-1)/nNoise)+1 ), ...
                   snr_l1l2(k, loops), ssim_l1l2(k, loops), psnr_l1l2(k, loops), l1l2_time(k, loops) );

%        if(k<=nNoisyIm)
%          str = sprintf(' %1.0f/255 -- %1.1f \t %1.2f \t %1.2f \t %1.2f    %1.2f\n', ...
%                     255*sigmaGauss( rem(k,NGauss) + NGauss*(rem(k,NGauss)==0) ), vNoise( floor((k-1)/nNoise)+1 ), ...
%                     snr_l1l2(k, loops), ssim_l1l2(k, loops), psnr_l1l2(k, loops), l1l2_time(k, loops) );
%        else
%          str = sprintf(' %1.0f/255 -- %1.1f \t %1.2f \t %1.2f \t %1.2f    %1.2f\n', ...
%                     255*sigmaGauss( rem(k,NGauss) + NGauss*(rem(k,NGauss)==0) ), vNoise( floor((k-1-nNoisyIm)/nNoise)+1 ), ...
%                     snr_l1l2(k, loops), ssim_l1l2(k, loops), psnr_l1l2(k, loops), l1l2_time(k, loops) );
%        end
      str_all = [ str_all str ];

      if(k==length(InCell))
        str = sprintf(' \n\n');
      end

    end % _END_ FOR(length(InCell))

  end % FOR(loops)

  disp(str_all)

end % FOR(t)  -- test Image

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[SSIM_CODE] = check4ssim()

  SSIM_CODE = exist('ssim_index');
  if( SSIM_CODE == 0 )
    disp('NOTE:');
    disp('  The function ssim_index (code for [SSIM]) is not in your path...');
    disp('  You may download it from:');
    disp('  http://www.ece.uwaterloo.ca/~z70wang/research/ssim/');
    disp('  ');
    disp('  Disabling SSIM reports');
    disp('  ');
    disp('  [SSIM]  Z. Wang, A. Bovik, H. Sheikh and E. Simoncelli, ');
    disp('         "Image quality assessment: From error visibility to');
    disp('         structural similarity " ');
    disp('         IEEE Transactions on Image Processing, 2004, 13:4(600-612).');

  end

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[Nloops, Img] = setupSim(example, Nloops, Img)

  tmpN = Nloops;
  tmpI = Img;

  if( strcmp(example, 'icassp12') )

    disp('  Running the "icassp11" example...');
    disp('  This script reproduces the Exprimental Results of ');
    disp('  "Mixed Gaussian-impulse Noise Image Restoration via');
    disp('    Total Variation"');
    disp(' ')
    disp('  This script takes sometime (ten-trial average)');
    disp(' ')
    tmpN = 10;

    if( ~strcmp(tmpI, 'all') )

      disp('NOTE:');
      disp('  If you run the "icassp12" example, you should set Img = "all" ...');
      disp('  setting Img = "all" ');
      tmpI = 'all';
    end

  end

  Nloops = tmpN;
  Img    = tmpI;

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[Ig Color] = selectImage(tImg)

  switch lower(tImg)

    case{'lena'}
      Ig = double( imread('gray_imgs/lena_gray_512.png') ) / 255;
      Color = 0;

    case{'peppers'}
      Ig = double( imread('gray_imgs/peppers_gray.png') ) / 255;
      Color = 0;

    case{'cman'}
      Ig = double( imread('gray_imgs/cameraman_256x256.tiff') ) / 255;
      Color = 0;

    case{'cman512'}
      Ig = double( imread('gray_imgs/cameraman_512x512.tiff') ) / 255;
      Color = 0;

    case{'goldhl'}
      Ig = double( imread('gray_imgs/goldhill_gray.png') ) / 255;
      Color = 0;

    case{'barb'}
      Ig = double( imread('gray_imgs/barbara_gray512.png') ) / 255;
      Color = 0;

    case{'clena'}
      Ig = double( imread('color_imgs/lena_color_512.png') ) / 255;
      SSIM_CODE = 0;
      Color = 1;

    case{'cbarb'}
      Ig = double( imread('color_imgs/barbara_color.png') ) / 255;
      SSIM_CODE = 0;
      Color = 1;

    case{'cboats'}
      Ig = double( imread('color_imgs/boats_color.png') ) / 255;
      SSIM_CODE = 0;
      Color = 1;

    case{'cgoldhl'}
      Ig = double( imread('color_imgs/goldhill_color.png') ) / 255;
      SSIM_CODE = 0;
      Color = 1;

    case{'none'}
      disp(' ');
      disp('Select a test image (see code for an example)...');
      disp(' ');
      disp('Exiting ICASSP12 test code...');
      return;

    otherwise
      error('Not a valid image\n');

  end  % _END_ SWITCH(Img)

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[lambdaS_l2, Lambda_O, Lambda_update] = setRegParameter(example, Ndims)

    switch lower(example)

    case{'gauss_sp_denoise'}
      if(Ndims == 1)
        lambdaS_l2 = [ 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025 ];
      else
        lambdaS_l2 = 1.5*[ 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025 ];
      end

      Lambda_O = 5;
      Lambda_update = 1.15;

    case{'gauss_random_denoise'}
      if(Ndims == 1)
        lambdaS_l2 = [ 1.25*0.0075, 0.0185, 0.025, ...
                 1.25*0.0075, 0.0185, 0.025, ...
                 1.25*0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025 ];
%                   0.0075, 0.0185, 0.025 ]/8;
      else
        lambdaS_l2 = 1.5*[ 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025, ...
                 0.0075, 0.0185, 0.025 ]/8;
      end

      Lambda_O = 5;
      Lambda_update = 1.15;

    end

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[InCell, noiseName, vNoise, nNoise] = addNoise(Ig, example, sigmaGauss, NGauss, spnoise, Nspnoise, randomnoise, NRnoise)

    % Add Gaussian noise
    Ig_05L2 = imnoise(Ig,'gaussian', 0, (sigmaGauss(1)^2) );
    Ig_10L2 = imnoise(Ig,'gaussian', 0, (sigmaGauss(2)^2) );
    Ig_15L2 = imnoise(Ig,'gaussian', 0, (sigmaGauss(3)^2) );

    % Add impulse noise (S&P or random value)
    switch lower(example)

      case{'gauss_sp_denoise'}

	noiseName = 'S&P';
	vNoise    = spnoise;
	nNoise    = Nspnoise;

	for k=1:Nspnoise,

	  InCell{NGauss*(k-1)+1} = imnoise(Ig_05L2, 'salt & pepper', spnoise(k));
	  InCell{NGauss*(k-1)+2} = imnoise(Ig_10L2, 'salt & pepper', spnoise(k));
	  InCell{NGauss*(k-1)+3} = imnoise(Ig_15L2, 'salt & pepper', spnoise(k));
	end

      case{'gauss_random_denoise'}

	noiseName = 'RandVal';
	vNoise    = randomnoise;
	nNoise    = NRnoise;

	[Nrows Ncols Ndims] = size(Ig);
	dummy = 0.5*ones(Nrows, Ncols, Ndims);

	for k=1:NRnoise,

	  mask = (imnoise(dummy, 'salt & pepper', randomnoise(k)) ~= 0.5);
	  InCell{NGauss*(k-1)+1} = Ig_05L2.*(1-mask) + mask.*rand(Nrows, Ncols, Ndims);

	  mask = (imnoise(dummy, 'salt & pepper', randomnoise(k)) ~= 0.5);
	  InCell{NGauss*(k-1)+2} = Ig_10L2.*(1-mask) + mask.*rand(Nrows, Ncols, Ndims);

	  mask = (imnoise(dummy, 'salt & pepper', randomnoise(k)) ~= 0.5);
	  InCell{NGauss*(k-1)+3} = Ig_15L2.*(1-mask) + mask.*rand(Nrows, Ncols, Ndims);
	end

    end % _END_ SWITCH

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[S0] = estimateRValImpulse(In)

          thres_diffmed= 0.2;
          window_side= 2;

          road_size= 7;
          thres_ranked= 0.05;

          [Nrows Ncols Ndims] = size(In);

          if(Ndims == 1)

            mask_diff_median = diff_median_blk( In, window_side, thres_diffmed);
            mask_ranked_mean = ranked_mean( In, window_side, road_size, thres_ranked);
            S0 = mask_ranked_mean & mask_diff_median;

          else

            for d=1:Ndims,
              mask_diff_median = diff_median_blk( In(:,:,d), window_side, thres_diffmed);
              mask_ranked_mean = ranked_mean( In(:,:,d), window_side, road_size, thres_ranked);
              S0(:,:,d) = mask_ranked_mean & mask_diff_median;
            end

          end

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[psnr_l1l2, snr_l1l2, ssim_l1l2] = qualityAssessment(Ig, Ig_L1L2, SSIM_CODE)

      psnr_l1l2 = psnr(Ig, Ig_L1L2);
      snr_l1l2  = snr(Ig, Ig_L1L2);
      if(SSIM_CODE)
        ssim_l1l2 = ssim_index(255*Normalize(Ig), 255*Normalize(Ig_L1L2));
      else
        ssim_l1l2 = NaN;
      end

end

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

function[y] = Normalize(x)

  y = (x - min(x(:)))/(max(x(:)) - min(x(:)));

end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [k,v]=v_findpeaks(x,m,w)
%V_FINDPEAKS finds peaks with optional quadratic interpolation [K,V]=(X,M,W)
%
%  Inputs:  X        is the input signal (does not work with UInt datatype)
%           M        is mode:
%                       'q' performs quadratic interpolation
%                       'v' finds valleys instead of peaks
%           W        is the width tolerance; a peak will be eliminated if there is
%                    a higher peak within +-W samples
%
% Outputs:  K        are the peak locations in X (fractional if M='q')
%           V        are the peak amplitudes: if M='q' the amplitudes will be interpolated
%                    whereas if M~='q' then V=X(K).

% Outputs are column vectors regardless of whether X is row or column.
% If there is a plateau rather than a sharp peak, the routine will place the
% peak in the centre of the plateau. When the W input argument is specified,
% the routine will eliminate the lower of any pair of peaks whose separation
% is <=W; if the peaks have exactly the same height, the second one will be eliminated.
% All peak locations satisfy 1<K<length(X).
%
% If no output arguments are specified, the results will be plotted.
%

%	   Copyright (C) Mike Brookes 2005
%      Version: $Id: v_findpeaks.m 3601 2013-10-11 15:27:30Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2
    m=' ';
end
nx=length(x);
if any(m=='v')
    x=-x(:);        % invert x if searching for valleys
else
    x=x(:);        % force to be a column vector
end
dx=diff(x);
r=find(dx>0);
f=find(dx<0);

if ~isempty(r) && ~isempty(f)    % we must have at least one rise and one fall
    dr=r;
    dr(2:end)=diff(r);
    rc=ones(nx,1);
    rc(r+1)=1-dr;
    rc(1)=0;
    rs=cumsum(rc); % = time since the last rise

    df=f;
    df(2:end)=diff(f);
    fc=ones(nx,1);
    fc(f+1)=1-df;
    fc(1)=0;
    fs=cumsum(fc); % = time since the last fall

    rp=repmat(-1,nx,1);
    rp([1; r+1])=[dr-1; nx-r(end)-1];
    rq=cumsum(rp);  % = time to the next rise

    fp=repmat(-1,nx,1);
    fp([1; f+1])=[df-1; nx-f(end)-1];
    fq=cumsum(fp); % = time to the next fall

    k=find((rs<fs) & (fq<rq) & (floor((fq-rs)/2)==0));   % the final term centres peaks within a plateau
    v=x(k);

    if any(m=='q')         % do quadratic interpolation
        b=0.5*(x(k+1)-x(k-1));
        a=x(k)-b-x(k-1);
        j=(a>0);            % j=0 on a plateau
        v(j)=x(k(j))+0.25*b(j).^2./a(j);
        k(j)=k(j)+0.5*b(j)./a(j);
        k(~j)=k(~j)+(fq(k(~j))-rs(k(~j)))/2;    % add 0.5 to k if plateau has an even width
    end

    % now purge nearby peaks

    if nargin>2
        j=find(diff(k)<=w);
        while any(j)
            j=j+(v(j)>=v(j+1));
            k(j)=[];
            v(j)=[];
            j=find(diff(k)<=w);
        end
    end
else
    k=[];
    v=[];
end
if any(m=='v')
    v=-v;    % invert peaks if searching for valleys
end
if ~nargout
    if any(m=='v')
        x=-x;    % re-invert x if searching for valleys
        ch='v';
    else
        ch='^';
    end
    plot(1:nx,x,'-',k,v,ch);
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%% isrounded
% Below is a demonstration of the features of the |isrounded| function

%%
clear; close all; clc;

%% Syntax
% |L=isrounded(A);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function pass = test_constructor(pref)
% Test the Chebfun3v constructor when performing simple arithmetic
% operations.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 100*pref.cheb3Prefs.chebfun3eps;

doms = [-1 1 -1 1 -1 1;
       -2 3 0 1 2 4];

for j = 1:size(doms, 1)
    dom = doms(j, :);
    % Check the constructor works with lots of different syntax:
    f1 = @(x,y,z) cos(x.*y.*z);
    f2 = @(x,y,z) sin(x+z) + cos(y);
    f3 = @(x,y,z) exp(x) - exp(-y) + z;
    g1 = chebfun3(f1, dom);
    g2 = chebfun3(f2, dom);
    g3 = chebfun3(f3, dom);

    % Two components:
    H1 = chebfun3v(f1, f2, dom);
    H2 = chebfun3v(f1, g2, dom);
    H3 = chebfun3v(g1, f2, dom);
    H4 = chebfun3v(g1, g2, dom);

    pass(j, 1) = norm(H1 - H2) < tol;
    pass(j, 2) = norm(H1 - H3) < tol;
    pass(j, 3) = norm(H1 - H4) < tol;

    % Three components:
    H1 = chebfun3v(f1, f2, f3, dom);
    H2 = chebfun3v(f1, f2, g3, dom);
    H3 = chebfun3v(f1, g2, f3, dom);
    H4 = chebfun3v(g1, f2, f3, dom);
    H5 = chebfun3v(g1, g2, f3, dom);
    H6 = chebfun3v(g1, f2, g3, dom);
    H7 = chebfun3v(f1, g2, g3, dom);
    H8 = chebfun3v(g1, g2, g3, dom);

    pass(j, 4) = norm(H1 - H2) < tol;
    pass(j, 5) = norm(H1 - H3) < tol;
    pass(j, 6) = norm(H1 - H4) < tol;
    pass(j, 7) = norm(H1 - H5) < tol;
    pass(j, 8) = norm(H1 - H6) < tol;
    pass(j, 9) = norm(H1 - H7) < tol;
    pass(j, 10) = norm(H1 - H8) < tol;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%...!
%/*************************************************************************
% Function Name : wrl2mat
% Author: Alireza Bossaghzadeh

    % This Function Extract Pointcloud from a wrl file and change it to a 2D
    % mesh matrix
    % In the process of changing to mesh, it will interpolated the hole parts
    % (Parts without any value)

%/*************************************************************************
    % Example: [ZI pointcloud]=wrl2mat('mywrl.wrl');
%     Input:         fname 	 wrl File Name
%     Outputs:       ZI          	 2D mesh File
%                    pointcloud 	 3xN  pointcloud data
    % Suggestion:
    % to have a better output you should use crop3dface.m that is available in my
    % files in mathworks file exchange.

% In the case of any problem you can call me by
% Email:Alibossagh@yahoo.co.uk
% This code is a modified version of the read_vrml.m written by "G. Akroyd"

% Version:    1.00       Published: 2008 June 07

function [ZI pointcloud]=wrl2mat(fname)

disp(['Reading The File  ', fname])
        vrfile=fopen(fname);
        p=1;counter=0;pointcloud=[0 0 0];
 while counter~=-1
    data=fgets(vrfile);
    fpoint=findstr(data,'point');% 2 checkers to find out the begining
    f2point=findstr(data,'[');   %of the x,y,z " point [ "
    while ~isempty(fpoint) & ~isempty(f2point)
      data=fgets(vrfile);
        if isempty(findstr(data,']'));
          t=sscanf(data,'%f %f %f');
          pointcloud(p,:)=t';
          p=p+1;
        else fpoint=[];counter=-1;
        end
    end
 end
 disp(['Convert to 2D Image'])
X=pointcloud(:,1);Y=pointcloud(:,2);Z=pointcloud(:,3);%Load positions to X,Y,Z
X=X-min(X);Y=Y-min(Y);%change X,Y to start from 0
x=[min(X):1:max(X)];y=[min(Y):1:max(Y)];%create image dimentions
[XI,YI] = meshgrid(x,y);%create matrix of positions
YI=max(YI(:))-YI;
Z=Z-min(Z(:));%Change Z data to start from 0
ZI=griddata(X,Y,Z,XI,YI,'linear'); % calculate Z values from noninteger positions
    figure;imshow(ZI,[]);title('Extracted Image')


text:
%% Sachs Model
%
%%

author: mtex-toolbox

repo: mtex

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


text:
function Mi = get_block_col(M, C, col_range)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: 1/27/2016 2:39:32 AM
	% Last modified	: 1/27/2016 2:39:36 AM
	% Description	: Get column blocks of a big matrix M, blocks indexed by C
	% 	INPUT
	%		M        : the big matrix. M = [M1 , M2 , ... , MC]
	%		C        : block indexes
	%		col_range: a vector store the last index of each block. col_range(1) = 0.
	%					i-th block is indexed by col_range(i)+1: col_range(i+1).
	% 	OUTPUT
	%		Mi: output block matrix
	%
	%% ================== end File info ==========================

	id_sel = [];
	for i = 1: numel(C)
		c = C(i);
		id_sel = [id_sel, col_range(c) + 1: col_range(c+1)];
	end
	Mi = M(:, id_sel, :);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [ key, seed, matrix, ierror ] = rcont2 ( nrow, ncol, nrowt, ncolt, ...
  key, seed )

%*****************************************************************************80
%
%% RCONT2 constructs a random two-way contingency table with given sums.
%
%  Discussion:
%
%    It is possible to specify row and column sum vectors which
%    correspond to no table at all.  As far as I can see, this routine does
%    not detect such a case.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2009
%
%  Author:
%
%    Original FORTRAN77 version by WM Patefield.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    WM Patefield,
%    Algorithm AS 159:
%    An Efficient Method of Generating RXC Tables with
%    Given Row and Column Totals,
%    Applied Statistics,
%    Volume 30, Number 1, 1981, pages 91-97.
%
%  Parameters:
%
%    Input, integer NROW, NCOL, the number of rows and columns
%    in the table.  NROW and NCOL must each be at least 2.
%
%    Input, integer NROWT(NROW), NCOLT(NCOL), the row and column
%    sums.  Each entry must be positive.
%
%    Input, logical KEY, a flag that indicates whether data has
%    been initialized for this problem.  Set KEY = .FALSE. before the first
%    call.
%
%    Input, integer SEED, a seed for the random number
%    generator.
%
%    Output, logical KEY, a flag that indicates whether data has
%    been initialized for this problem.  Set KEY = .FALSE. before the first
%    call.
%
%    Output, integer SEED, a seed for the random number
%    generator.
%
%    Output, integer MATRIX(NROW,NCOL), the matrix.
%
%    Output, integer IERROR, an error flag, which is returned
%    as 0 if no error occurred.
%

  persistent fact;
  persistent ntotal;

  ierror = 0;
%
%  On user's signal, set up the factorial table.
%
  if ( ~key )

    key = 1;

    if ( nrow <= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'RCONT - Fatal error!\n' );
      fprintf ( 1, '  Input number of rows is less than 2.\n' );
      ierror = 1;
      return
    end

    if ( ncol <= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'RCONT - Fatal error!\n' );
      fprintf ( 1, '  The number of columns is less than 2.\n' );
      ierror = 2;
      return
    end

    for i = 1 : nrow
      if ( nrowt(i) <= 0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'RCONT - Fatal error!\n' );
        fprintf ( 1, '  An entry in the row sum vector is not positive.\n' );
        ierror = 3;
        return
      end
    end

    for j = 1 : ncol
      if ( ncolt(j) <= 0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'RCONT - Fatal error!\n' );
        fprintf ( 1, '  An entry in the column sum vector is not positive.\n' );
        ierror = 4;
        return
      end
    end

    if ( sum ( ncolt(1:ncol) ) ~= sum ( nrowt(1:nrow) ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'RCONT - Fatal error!\n' );
      fprintf ( 1, '  The row and column sum vectors do not have the same sum.\n' );
      ierror = 6;
      return
    end

    ntotal = sum ( ncolt(1:ncol) );

    fact = zeros(ntotal,1);
%
%  Calculate log-factorials.
%
    x = 0.0;
    fact(1) = 0.0;
    for i = 1 : ntotal
      x = x + log ( i );
      fact(i+1) = x;
    end

  end
%
%  Construct a random matrix.
%
  jwork(1:ncol-1) = ncolt(1:ncol-1);

  jc = ntotal;

  for l = 1 : nrow - 1

    nrowtl = nrowt(l);
    ia = nrowtl;
    ic = jc;
    jc = jc - nrowtl;

    for m = 1 : ncol - 1

      id = jwork(m);
      ie = ic;
      ic = ic - id;
      ib = ie - ia;
      ii = ib - id;
%
%  Test for zero entries in matrix.
%
      if ( ie == 0 )
        ia = 0;
        matrix(l,m:ncol) = 0;
        break
      end
%
%  Generate a pseudo-random number.
%
      [ r, seed ] = r8_uniform_01 ( seed );
%
%  Compute the conditional expected value of MATRIX(L,M).
%
      done1 = 0;

      while ( 1 );

        nlm = floor ( ia * id / ie + 0.5 );
        iap = ia + 1;
        idp = id + 1;
        igp = idp - nlm;
        ihp = iap - nlm;
        nlmp = nlm + 1;
        iip = ii + nlmp;
        x = exp ( fact(iap) + fact(ib+1) + fact(ic+1) + fact(idp) - ...
          fact(ie+1) - fact(nlmp) - fact(igp) - fact(ihp) - fact(iip) );

        if ( r <= x )
          break;
        end

        sumprb = x;
        y = x;
        nll = nlm;
        lsp = 0;
        lsm = 0;
%
%  Increment entry in row L, column M.
%
        while ( ~lsp )

          j = ( id - nlm ) * ( ia - nlm );

          if ( j == 0 )

            lsp = 1;

          else

            nlm = nlm + 1;
            x = x * j / ( nlm * ( ii + nlm ) );
            sumprb = sumprb + x;

            if ( r <= sumprb )
              done1 = 1;
              break
            end

          end

          done2 = 0;

          while ( ~lsm )
%
%  Decrement the entry in row L, column M.
%
            j = nll * ( ii + nll );

            if ( j == 0 )
              lsm = 1;
              break
            end

            nll = nll - 1;
            y = y * j / ( ( id - nll ) * ( ia - nll ) );
            sumprb = sumprb + y;

            if ( r <= sumprb )
              nlm = nll;
              done2 = 1;
              break
            end

            if ( ~lsp )
              break
            end

          end

          if ( done2 )
            break
          end

        end

        if ( done1 )
          break
        end

        if ( done2 )
          break
        end

        [ r, seed ] = r8_uniform_01 ( seed );
        r = sumprb * r;

      end

      matrix(l,m) = nlm;
      ia = ia - nlm;
      jwork(m) = jwork(m) - nlm;

    end

    matrix(l,ncol) = ia;

  end
%
%  Compute the last row.
%
  matrix(nrow,1:ncol-1) = jwork(1:ncol-1);
  matrix(nrow,ncol) = ib - matrix(nrow,ncol-1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Example_6b_GBAY_w_floodplain
% Continue on from Example_6_GBAY.m by building on
% a floodplain onto the mesh.

% Most mesh sizing functions can be enforced in a
% topographic elevation (-1*depth) range such as:
%
% sizing_parameter = [size_value1, z_min1, z_max1;
%                     size_value2, z_min2, z_max2;
%                    ];
% Note: multiple elevation ranges must be delineated with semi-colons ;

clearvars; clc;

addpath('..')
addpath(genpath('../utilities/'))
addpath(genpath('../datasets/'))
addpath(genpath('../m_map/'))

%% STEP 1: set mesh extents and set parameters for mesh.
min_el    = 60;  		     % Minimum mesh resolution in meters.
max_el    = [1e3,-inf,0 ;    % Underwater maximum mesh resolution in meters.
             500,0,+inf];    % Overland maximum mesh resolution in meters.
grade     = [0.25,-inf,0;    % Underwater gradation rate
             0.05,0,+inf] ;  % Overland gradation rate
angleOfReslope = 60 ;       % Control width of channel by changing angle of reslope.
ch = 0.1 ;                  % Scale resolution propotional to depth nearby thalweg.
fs = 3 ;                    % Place 3 vertices per width of shoreline feature.
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
coastline = 'us_medium_shoreline_polygon';
demfile   = 'galveston_13_mhw_2007.nc';

gdatuw = geodata('shp',coastline,...
                 'dem',demfile,...
                 'h0',min_el);

load ECGC_Thalwegs.mat % Load the Channel thalweg data
%% STEP 3: create an edge function class
fh = edgefx('geodata',gdatuw,...
            'fs',fs,...
            'ch',ch,...
            'AngOfRe',angleOfReslope,...% control the width
            'Channels',pts2,...
            'g',grade,...
            'max_el',max_el);
%% STEP 4: Pass your edgefx class object along with some meshing options and
% build the mesh...
mshopts = meshgen('ef',fh,'bou',gdatuw,'plot_on',1,'proj','lambert');
% now build the mesh with your options and the edge function.
mshopts = mshopts.build;
%% STEP 5: Get fixed constraints and update gdat with overland meshing domain.
muw = mshopts.grd ;
muw = makens(muw,'auto',gdatuw) ; % apply so that extractFixedConstraints only grabs the shoreline constraints.

[pfix,egfix] = extractFixedConstraints(muw) ;

% 10-m contour extracted from the Coastal Relief Model.
coastline = 'us_coastalreliefmodel_10mLMSL';
demfile   = 'galveston_13_mhw_2007.nc';
landuse   = 'galveston_2016_CCAP.nc';

gdat = geodata('shp',coastline,...
               'dem',demfile,...
               'h0',min_el);

%gdat.inpoly_flip =  mod(1,gdat.inpoly_flip) ; % if the meshing domain is inverted, you can always flip it .

% Here we pass our constraints to the mesh generator (pfix and egfix).
mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'proj','lambert',...
                  'pfix',pfix,'egfix',egfix);

% now build the mesh with your options and the edge function.
mshopts = mshopts.build;

m = mshopts.grd ;

% plot resolution on the mesh
plot(m,'type','resomesh','colormap',[10 0 1e3])

% interpolate bathy using special constraining technique
% for overland and underwater
m = interpFP(m,gdat,muw,gdatuw);
plot(m,'type','bmesh') % plot the bathy on the mesh

% computing mannings based on CCAP landcover data using
% cell-averaged interpolation with stencil 4*grid_size for stability
m = Calc_Mannings_Landcover(m,landuse,'ccap','N',4);
plot(m,'type','mann') % plot the mannings on the mesh

author: CHLNDDEV


text:
function [y,pm0]=anaqpsk(N,Ncomp,f0);
%ANAQPSK Quaternary Phase Shift Keying (QPSK) signal.
% 	[Y,PM]=ANAQPSK(N,NCOMP,F0) returns a complex phase modulated signal
% 	of normalized frequency F0, whose phase changes every NCOMP point according
%	to a discrete uniform law, between the values (0, pi/2, pi, 3*pi/2).
% 	Such signal is only 'quasi'-analytic.
%
%	N     : number of points
%	NCOMP : number of points of each component (default: N/5)
% 	F0    : normalized frequency.              (default: 0.25)
% 	Y     : signal
% 	PM0   : initial phase of each component	   (optional).
%
%	Example :
%	 [signal,pm0]=anaqpsk(512,64,0.05); clf; figure(gcf);
%  	 subplot(211); plot(real(signal)); subplot(212); plot(pm0);
%
%	See also ANAFSK, ANABPSK, ANAASK.

%	O. Lemoine - October 1995
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('The number of parameters must be at least 1.');
elseif (nargin == 1),
 Ncomp=round(N/5); f0=0.25;
elseif (nargin == 2),
 f0=0.25;
end;

if (N <= 0),
 error('The signal length N must be strictly positive' );
elseif (f0<0)|(f0>0.5),
 error('f0 must be between 0 and 0.5');
end;

MatlabVersion=version; MatlabVersion=str2num(MatlabVersion(1));
if (MatlabVersion==4), rand('uniform'); end;

m=ceil(N/Ncomp);

% jumps=round(3*rand(m,1));
% This is a modification proposed by Alpesh Patel, McMaster University
jumps=floor(4*rand(m,1)); jumps(jumps==4)=3;

pm0=pi*kron(jumps,ones(Ncomp,1))/2; pm0=pm0(1:N,1);
tm=(1:N)'-1;
pm=(2.0*pi*f0*tm+pm0);

y = exp(j*pm);

author: HeLiangHIT


text:
% Test for @chebfun/comet3.m

function pass = test_comet3(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

% Real scalar functions.
fsr1 = chebfun(@sin, [-1 0 1], pref);
fsr2 = chebfun(@cos, [-1 0.5 1], pref);
fsr3 = chebfun(@exp, [-1 -0.5 1], pref);

% Real array-valued functions.
far1 = chebfun(@(x) [sin(x) cos(x) exp(x)], [-1 0 1], pref);

% Real quasimatrices.
fqr1 = cheb2quasi(far1);

% Obviously, we can't check if the plots are correct without human
% intervention, so all these tests are meant to do is make sure that we don't
% crash.

hfig = figure('Visible', 'off');

% Check COMET and COMET3.
pass(1) = doesNotCrash(@() comet3(fsr1, fsr2, fsr3));
pass(2) = doesNotCrash(@() comet3(far1));
pass(3) = doesNotCrash(@() comet3(fqr1));

close(hfig);

end

function pass = doesNotCrash(fn)

try
    fn();
    pass = true;
catch ME;
    pass = false;
end

end


text:
function particles = prediction_step(particles, u, noise)
% Updates the particles by drawing from the motion model
% Use u.r1, u.t, and u.r2 to access the rotation and translation values
% which have to be pertubated with Gaussian noise.
% The position of the i-th particle is given by the 3D vector
% particles(i).pose which represents (x, y, theta).

% noise parameters
% Assume Gaussian noise in each of the three parameters of the motion model.
% These three parameters may be used as standard deviations for sampling.
r1Noise = noise(1);
transNoise = noise(2);
r2Noise = noise(3);

numParticles = length(particles);

for i = 1:numParticles

  % append the old position to the history of the particle
  particles(i).history{end+1} = particles(i).pose;

  % sample a new pose for the particle
  r1 = normrnd(u.r1, r1Noise);
  r2 = normrnd(u.r2, r2Noise);
  trans = normrnd(u.t, transNoise);
  particles(i).pose(1) = particles(i).pose(1) + trans*cos(particles(i).pose(3) + r1);
  particles(i).pose(2) = particles(i).pose(2) + trans*sin(particles(i).pose(3) + r1);
  particles(i).pose(3) = normalize_angle(particles(i).pose(3) + r1 + r2);
end

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function [] = exportInputsMustOrder2ToGAMS(model, minFluxesW, maxFluxesW, constrOpt, excludedRxns, mustSetFirstOrder, inputFolder)
% This function exports all the inputs needed to run the GAMS functions to
% find second order Must Sets (`MustUU`, `MustLL`, `MustUL`). The inputs will
% be stored in inputFolder. Some inputs will be exported using GDXMRW and
% others will be exported as simple .txt files.
%
% USAGE:
%
%    exportInputsMustOrder2ToGAMS(model, minFluxesW, maxFluxesW, constrOpt, excludedRxns, mustSetFirstOrder, inputFolder)
%
% INPUTS:
%    model:                   (structure) a metabolic model with at least
%                             the following fields:
%
%                               * .rxns - Reaction IDs in the model
%                               * .mets - Metabolite IDs in the model
%                               * .S -    Stoichiometric matrix (sparse)
%                               * .b -    RHS of `Sv = b` (usually zeros)
%                               * .c -    Objective coefficients
%                               * .lb -   Lower bounds for fluxes
%                               * .ub -   Upper bounds for fluxes
%    minFluxesW:              (double array of size `n_rxns x 1`) minimum
%                             fluxes for each reaction in the model for
%                             wild-type strain. This can be obtained by
%                             running the function `FVAOptForce` e.g.:
%                             `minFluxesW = [-90; -56];`
%    maxFluxesW:              (double array of size `n_rxns x 1`) maximum
%                             fluxes for each reaction in the model for
%                             wild-type strain. This can be obtained by
%                             running the function `FVAOptForce` e.g.:
%                             `maxFluxesW = [90; 56];`
%    constrOpt:               (structure) Structure containing additional
%                             contraints. Include here only reactions
%                             whose flux is fixed, i.e., reactions whose
%                             lower and upper bounds have the same value.
%                             Do not include here reactions whose lower
%                             and upper bounds have different values. Such
%                             contraints should be defined in the lower
%                             and upper bounds of the model. The structure
%                             has the following fields:
%
%                               * .rxnList - Reaction list (cell array)
%                               * .values -  Values for constrained reactions
%                                 (double array). E.g.: `struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');`
%    excludedRxns:            (cell array) Reactions to be excluded to the
%                             `MustUL` set. This could be used to avoid
%                             finding transporters or exchange reactions
%                             in the set. Default: empty.
%    mustSetFirstOrder:       (cell array) Reactions that belong to `MustU`
%                             and `MustL` (first order sets). Default: empty.
%    inputFolder:             (string) Folder where inputs for GAMS
%                             function will be stored
%
% OUTPUTS:
%    Reactions.txt:           (file) File containing the
%                             identifiers for reactions
%    Metabolites.txt:         (file) File containing the
%                             identifiers for metabolites
%    Constrains.txt:          (file) File containing the
%                             identifiers for constrained reactions
%    Excluded.txt:            (file) File containing the
%                             identifiers for excluded reactions. These
%                             reactions will not be considered in when
%                             running `findMustXX.gms` (`XX = UU` or `LL` or `UL`
%                             depending on the case)
%    MustSetFirstOrder.txt:   (file) File containing the
%                             identifiers for reactions in `MustL` and
%                             `MustU`
%    MtoG.gdx:                (file) File containing the
%                             parameters which will be read by GAMS
%                             (lower bounds, upper bounds, stoichiometrix
%                             matrix `S`, minimum and maximun fluxes for
%                             each reaction in the previous step of FVA,
%                             and the values for contrained reactions)
%
% .. Author: - Sebastian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl

if nargin < 8  % inputs handling
    error('Optforce: All inputs must be specified when using exportInputsMustOrder2ToGAMS');
end

%Create a temporaty folder for inputs
if ~exist(inputFolder, 'dir')
   mkdir(inputFolder);
end
current = pwd;
cd(inputFolder);

%Export Sets
exportSetToGAMS(model.rxns, 'Reactions.txt')
exportSetToGAMS(model.mets, 'Metabolites.txt')
exportSetToGAMS(constrOpt.rxnList, 'Constraints.txt')
exportSetToGAMS(excludedRxns, 'Excluded.txt')
exportSetToGAMS(mustSetFirstOrder, 'MustSetFirstOrder.txt')

%Export parameters
s.name = 's';
s.val = full(model.S);
s.type = 'parameter';
s.form = 'full';
s.uels = {model.mets',model.rxns'};

basemin.name = 'basemin';
basemin.val = minFluxesW;
basemin.type = 'parameter';
basemin.form = 'full';
basemin.uels = model.rxns';

basemax.name = 'basemax';
basemax.val = maxFluxesW;
basemax.type = 'parameter';
basemax.form = 'full';
basemax.uels = model.rxns';

lb.name = 'lb';
lb.val = model.lb';
lb.type = 'parameter';
lb.form = 'full';
lb.uels = model.rxns';

ub.name = 'ub';
ub.val = model.ub';
ub.type = 'parameter';
ub.form = 'full';
ub.uels = model.rxns';

b.name = 'b';
b.val = constrOpt.values;
b.type = 'parameter';
b.form = 'full';
b.uels = {'EX_gluc','R75','EX_suc'};

%Using GDXMRW to export inputs
fileName = ['MtoG' setType];
wgdx(fileName, s, basemin, basemax, lb, ub, b)

cd(current);

end

author: opencobra

repo: cobratoolbox