text:
function varargout = read_meas_prot(filename, varargin)
%READ_MEAS_PROT  read pulse sequence protocol from VB13A- and VB15A-style "meas.dat"
%
% YAPS = read_meas_prot(filename)

% jonathan polimeni <jonp@padkeemao.nmr.mgh.harvard.edu>, 2007/jan/03
% $Id: read_meas_prot.m,v 1.1 2008/05/21 19:01:01 greve Exp $
%**************************************************************************%

  VERSION = '$Revision: 1.1 $';
  if ( nargin == 0 ), help(mfilename); return; end;

  header = '';
  if ( nargin >= 2 ),
    if ( ~isstruct(varargin{1}) ),
      header = varargin{1};
    end;
  end;

  %------------------------------------------------------------------------%
  % error checking

  if ( ~exist(filename, 'file') ),
    error('file [%s] does not exist', filename);
  end;

  VERBOSE = 0;

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

  if ( isempty(header) ),
    [fp, errstr] = fopen(filename, 'r', 'l');
    if ( fp == -1 ),
      error(errstr);
    end;

    % determine size (in bytes) of ascii header files stored in the 'meas.dat'
    % format (i.e., "Config_.evp", "Dicom_.evp", etc) to skip over them all.
    % [note: this works for VB11A also---the first integer is 32, which is
    % the number of bytes to skip at the beginning!]
    data_start = fread(fp, 1, 'uint32');

    header_files = fread(fp, 1, 'uint32');

    % read header into one string for parsing
    header = fscanf(fp, '%c', data_start-24);

    fclose(fp);

  end;

  %------------------------------------------------------------------------%
  % ICE parameters from Config_.evp

  param_list = {
      ... %%% k-space sample counts  (no "RawSeg" for some reason)
      'RawCol','RawLin','RawCha','RawSet','RawEco','RawPhs','RawRep', ...
      'RawPar','RawSlc','RawIda','RawIdb','RawIdc','RawIdd','RawIde','RawAve', ...
      ... %%%
      ... %%% iPAT / parallel imaging acquisition parameters
      'NAFLin','NAFPar','NFirstLin','NFirstPar', ...
      'NRefLin','NRefPar','NFirstRefLin','NFirstRefPar', ...
      ... %%%
      ... %%% effective k-space dimensions
      'NColMeas','NLinMeas','NChaMeas','NSetMeas','NEcoMeas','NPhsMeas','NRepMeas','NSegMeas', ...
      'NParMeas','NSlcMeas','NIdaMeas','NIdbMeas','NIdcMeas','NIddMeas','NIdeMeas','NAveMeas', ...
      ... %%%
      ... %%% final image dimensions
      'NImageLins','NImageCols','NImagePar'};

  Config_evp = cell2struct(cell(length(param_list),1), param_list, 1);
  Config_evp_values = [];

  for ind = 1:length(param_list),
    param = param_list{ind};

    %%%%%%% MATCH TYPE #0: integer dimension size, e.g., <ParamLong."RawCol">  { 256  }
    match = regexp(header, ['<Param\w+\."(?<param>' param ')">\s*{\s*(?<value>[-]*\d*)\s*}'], 'names');

    % check if no match is found
    if ( isempty(match) ),
      if ( VERBOSE ),
        warning('SIEMENS:IO:versioning', 'missing protocol parameter "%s"---check sequence', param);
      end;
      continue;
    end;

    % consider only first match
    match = match(1);

    % empty means value for this parameter dimension = 0
    if ( isempty(match.value) ),
      match.value = '0';
    end;

    % save out struct and numerical array
    param_values(ind) = str2num(match.value);

    %%%-------------------------------%
    %%% SPECIAL CASES: "NFirstLin", "NFirstPar", "NFirstRefLin", "NFirstRefPar" (0-based indices needing conversion to 1-based)
    if ( regexp(param, '^NFirst') ),
      param_values(ind) = param_values(ind) + 1;
    end;

    Config_evp = setfield(Config_evp, param_list{ind}, param_values(ind));

  end;  %% FOR

  %  %%%-------------------------------%
  %  %%% SPECIAL CASES: "NFirstRefLin", "NFirstRefPar"
  %  if ( Config_evp.NFirstRefLin == 0 ),
  %    Config_evp.NRefLin = 0;
  %  end;
  %  if ( Config_evp.NFirstRefPar == 0 ),
  %    Config_evp.NRefPar = 0;
  %  end;

  %========================================================================%
  %------------------------------------------------------------------------%
  %========================================================================%
  % IDEA parameters from YAPS, a.k.a., ASCCONV (contained in meas.asc)

  %  TODO: determine phase encoding or readout direction, and polarities,
  %  to convert reconstructed matrix into proper viewing position

  %------------------------------------------------------------------------%
  % anastasia's list

  %  p.SequenceFileName      = 'tSequenceFileName';
  %  p.ProtocolName          = 'tProtocolName';
  %  p.FlipAngleDegree       = 'adFlipAngleDegree';
  %  p.TxFrequency           = 'sTXSPEC.lFrequency';
  %  p.BaseResolution        = 'sKSpace.lBaseResolution';
  %  p.PhaseEncodingLines    = 'sKSpace.lPhaseEncodingLines';
  %  p.FourierRows           = 'iNoOfFourierLines';
  %  p.FourierColumns        = 'iNoOfFourierColumns';
  %  p.NumChannels           = 'iMaxNoOfRxChannels';
  %  p.PhaseResolution       = 'sKSpace.dPhaseResolution';
  %  p.Partitions            = 'sKSpace.lPartitions';
  %  p.TR                    = 'alTR\[[0-9]*\]';
  %  p.TI                    = 'alTI\[[0-9]*\]';
  %  p.TE                    = 'alTE\[[0-9]*\]';
  %  p.SliceArraySize        = 'sSliceArray.lSize';
  %  p.SlicePositionSag      = 'sSliceArray.asSlice\[[0-9]*\].sPosition.dSag';
  %  p.SlicePositionCor      = 'sSliceArray.asSlice\[[0-9]*\].sPosition.dCor';
  %  p.SlicePositionTra      = 'sSliceArray.asSlice\[[0-9]*\].sPosition.dTra';
  %  p.SliceNormalSag        = 'sSliceArray.asSlice\[[0-9]*\].sNormal.dSag';
  %  p.SliceNormalCor        = 'sSliceArray.asSlice\[[0-9]*\].sNormal.dCor';
  %  p.SliceNormalTra        = 'sSliceArray.asSlice\[[0-9]*\].sNormal.dTra';
  %  p.Thickness             = 'sSliceArray.asSlice\[[0-9]*\].dThickness';
  %  p.PhaseFOV              = 'sSliceArray.asSlice\[[0-9]*\].dPhaseFOV';
  %  p.ReadoutFOV            = 'sSliceArray.asSlice\[[0-9]*\].dReadoutFOV';
  %  p.InPlaneRot            = 'sSliceArray.asSlice\[[0-9]*\].dInPlaneRot';
  %  p.DiffDirections        = 'sDiffusion.lDiffDirections';
  %  p.bValue                = 'sDiffusion.alBValue\[1\]';
  %  p.b0Repetitions         = 'sWiPMemBlock.alFree\[8\]';
  %  p.Repetitions           = 'lRepetitions';

  % dwelltime = The base resolution and the effective dwelltime are specified without oversampling
  % (true number of samples = base_resolution * oversampling,
  % real dwelltime = effective dwelltime/oversampling).

  % duration of the ADC (columns * dwelltime)

  % YAPS.lRampTime, YAPS.lFlatTime, YAPS.lADCDuration, YAPS.lRampMode

  [YAPS, param_list] = read_meas_prot__struct;

  param_values = zeros(length(param_list),1);

  % scan through header for each of the protocol parameter values
  for ind = 1:length(param_list),
    param = param_list{ind};

    %%%-------------------------------%
    %%% SPECIAL CASE:

    if ( strcmp(param, 'sCoilElementID_tElement') ),

      match = regexp(header, ['asCoilSelectMeas\S*(?<param>lRxChannelConnected)\s*=\s*(?<value>\d*\.?x?\d*)'], 'names');
      channel_val = str2double({match.value});
      [channel_unique, channel_index] = unique(channel_val);

      match = regexp(header, ['asCoilSelectMeas\S*(?<param>tElement)\s*=\s*"(?<string>\S*)"'], 'names');

      YAPS = setfield(YAPS, 'sCoilElementID_tElement', {match(channel_index).string});
      continue;

    end;

    %%%-------------------------------%
    %%% SPECIAL CASE: slice array params, "dPhaseFOV", "dReadoutFOV"

    if ( strcmp(param, 'dPhaseFOV') || strcmp(param, 'dReadoutFOV') ),
      param = strcat('asSlice\S*', param);
    end;

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

    % exploit MATLAB regexp machinery to pull out parameter/value pairs

    %%%%%%% MATCH TYPE #1: string parameter, e.g.,  tSequenceFileName                        = "%SiemensSeq%\gre"
    match = regexp(header, ['\W(?<param>' param ')\s*=\s*"(?<string>\S*)"'], 'names');
    if ( ~isempty(match) ),
      % consider only first match
      match = match(1);

      % set string value
      YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), match.string);
      continue;
    end;

    %%%%%%% MATCH TYPE #2: numerical parameter but single element, e.g., iMaxNoOfRxChannels                       = 9
    match = regexp(header, ['\W(?<param>' param ')\s*=\s*(?<value>(-*)\d*\.?x?\w*)'], 'names', 'once');

    if ( isempty(match) ),

      %%%%%%% MATCH TYPE #3: numerical parameter array, possibly negative, e.g., sRXSPEC.alDwellTime[0]                   = 6500
      match = regexp(header, ['\W(?<param>' param ')\[(?<index>\d*)\]\s*=\s*(?<value>(-*)\d*\.?\d*)'], 'names');

      if ( ~isempty(match) ),

        %%%-------------------------------%
        %%% SPECIAL CASE: "alRegridMode"

        % the "alRegridMode" parameter consists of an array of values indicating
        % which basic waveform components comprise the readout waveform. the
        % presence of a fifth component (i.e., "alRegridMode[4]") indicates a
        % sinusoidal waveform. the parameter value itself is irrelevant here.
        if ( strcmp(param, 'alRegridMode') ),
          sinusoid_gradient_waveform_flag = 0;
          for im = 1:length(match),
            if ( match(im).index == 4 ),
              % sinusoidal waveform found!
              % save out struct and numerical array
              sinusoid_gradient_waveform_flag = 1;
            end;
          end;
          if ( sinusoid_gradient_waveform_flag == 1 ),
            param_values(ind) = 4;  % 'REGRID_SINUSOIDAL'
          else,
            param_values(ind) = 2;  % 'REGRID_TRAPEZOIDAL'
          end;
          YAPS = setfield(YAPS, param_list{ind}, param_values(ind));
          continue;
        end;

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

        indices = str2num([match.index].');
        values = [];

        for ii = 1:length(match),
          % empty means value for this parameter = 0
          if ( isempty(match(ii).value) ),
            match(ii).value = '0';
          end;
          values(str2num(match(ii).index)+1) = str2num(match(ii).value);
        end;

        YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), values);
        continue;

      end;

    end;

    % check if still no match is found
    if ( isempty(match) ),
      if ( VERBOSE ),
        warning('SIEMENS:IO:versioning', 'missing protocol parameter "%s"---check sequence', param);
      end;

      if ( regexp(param, '^fl.*') )
	YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), 0.0);
      end;

      continue;
    end;

    %========================================================================%
    %========================================================================%
    % consider only first match
    match = match(1);

    % empty means value for this parameter = 0
    if ( isempty(match.value) ),
      match.value = '0';
    end;

    % convert hex entries to decimal
    if ( regexp(match.value, '^0x\d*') ),
      match.value = num2str(hex2dec( regexprep(match.value, '^0x', '') ));
    end;

    %%%-------------------------------%
    %%% ENUMERATED TYPE: "ucPhasePartialFourier" and "ucSlicePartialFourier"
    if ( strcmp(param, 'ucPhasePartialFourier') || strcmp(param, 'ucSlicePartialFourier') ),

      % "enum PartialFourierFactor" snarfed from <n4/pkg/MrServers/MrProtSrv/MrProt/SeqDefines.h>
      switch str2num(match.value),
       case hex2dec('01'),  % PF_HALF = 0x01
        match.value = '4/8';
       case hex2dec('02'),  % PF_5_8  = 0x02
        match.value = '5/8';
       case hex2dec('04'),  % PF_6_8  = 0x04
        match.value = '6/8';
       case hex2dec('08'),  % PF_7_8  = 0x08
        match.value = '7/8';
       case hex2dec('10'),  % PF_OFF  = 0x10
        match.value = '8/8';
       otherwise,
        warning(sprintf('unrecognized mode for "%s":  [ %d ]', param, match.value));
      end;

    end;

    %%%-------------------------------%
    %%% ENUMERATED TYPE: "ucMultiSliceMode"

    if ( strcmp(param, 'ucMultiSliceMode') ),
      % "enum MultiSliceMode" snarfed from <n4/pkg/MrServers/MrProtSrv/MrProt/SeqDefines.h>
      switch str2num(match.value),
       case hex2dec('01'),  % MSM_SEQUENTIAL  = 0x01
        match.value = 'MSM_SEQUENTIAL';
       case hex2dec('02'),  % MSM_INTERLEAVED = 0x02
        match.value = 'MSM_INTERLEAVED';
       case hex2dec('04'),  % MSM_SINGLESHOT  = 0x04
        match.value = 'MSM_SINGLESHOT';
       otherwise,
        warning(sprintf('unrecognized mode for "%s":  [ %d ]', param, match.value));
      end;

      YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), match.value);
      continue;
    end;

    %%%-------------------------------%
    %%% ENUMERATED TYPE: "ucPATMode"
    if ( strcmp(param, 'ucPATMode') ),

      % "enum PATSelMode" snarfed from <n4/pkg/MrServers/MrProtSrv/MrProt/SeqDefines.h>
      switch str2num(match.value),
       case hex2dec('01'),
        match.value = 'PAT_MODE_NONE';
       case hex2dec('02'),
        match.value = 'PAT_MODE_GRAPPA';
       case hex2dec('04'),
        match.value = 'PAT_MODE_SENSE';
       case hex2dec('08'),
        match.value = 'PAT_MODE_2D';
       otherwise,
        warning(sprintf('unrecognized mode for "%s":  [ %d ]', param, match.value));
      end;

      YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), match.value);
      continue;
    end;

    %%%-------------------------------%
    %%% ENUMERATED TYPE: "ucRefScanMode"
    if ( strcmp(param, 'ucRefScanMode') ),

      % "enum PATRefScanMode" snarfed from <n4/pkg/MrServers/MrProtSrv/MrProt/SeqDefines.h>
      switch str2num(match.value),
       case hex2dec('01'),
        match.value = 'PAT_REF_SCAN_UNDEFINED';      % e.g. if no PAT is selected
       case hex2dec('02'),
        match.value = 'PAT_REF_SCAN_INPLACE';        % sequence supplies inplace reference lines
       case hex2dec('04'),
        match.value = 'PAT_REF_SCAN_EXTRA';          % sequence supplies extra reference lines
       case hex2dec('08'),
        match.value = 'PAT_REF_SCAN_PRESCAN';        % sequence does not supply reference lines, the data must have been acquired with a previous measurement
       case hex2dec('10'),
        match.value = 'PAT_REF_SCAN_INTRINSIC_AVE';  % The sequence contains intrinsic ref.lines due to sharing e.g. in the averages dimension
       case hex2dec('20'),
        match.value = 'PAT_REF_SCAN_INTRINSIC_REP';  % The sequence contains intrinsic ref.lines due to sharing e.g. in the repetition or phases dimension (i.e., TSENSE)
       case hex2dec('40'),
        match.value = 'PAT_REF_SCAN_INTRINSIC_PHS';  % The sequence contains intrinsic ref.lines due to sharing e.g. in the repetition or phases dimension (i.e., TSENSE)
       case hex2dec('80'),
        match.value = 'PAT_REF_SCAN_INPLACE_LET';    % A single (L)ong (E)cho (T)rain acquires reference lines and imaging lines
       otherwise,
        warning(sprintf('unrecognized mode for "%s":  [ %d ]', param, match.value));
      end;

      YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), match.value);
      continue;
    end;

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

    % save out struct and numerical array
    param_values(ind) = str2num(match.value);
    YAPS = setfield(YAPS, regexprep(param_list{ind}, '\.', '_'), param_values(ind));

  end;  %% FOR

  %------------------------------------------------------------------------%
  %------------------------------------------------------------------------%
  %%% post-processing

  if ( ~strncmp(YAPS.ulVersion, '0x', 2) ),
    YAPS.ulVersion = lower(['0x' dec2hex(YAPS.ulVersion)]);
  end;

  switch YAPS.ulVersion,
   case '0xbee332',
    YAPS.ulVersion = [YAPS.ulVersion, ':  "N4_VA25A_LATEST_20040724"'];
   case '0x1421cf5',
    YAPS.ulVersion = [YAPS.ulVersion, ':  "N4_VB11D_LATEST_20040724"'];
   case '0x1452a3b',
    YAPS.ulVersion = [YAPS.ulVersion, ':  "N4_VB13A_LATEST_20060607"'];
   case '0x1483779',
    YAPS.ulVersion = [YAPS.ulVersion, ':  "N4_VB15A_LATEST_20070519"'];
   otherwise,
    disp(sprintf('==> [%s]: version string "%s" not on record!', mfilename, ...
                 YAPS.ulVersion));
  end;

  %%%-------------------------------%
  %%% SPECIAL CASE: "alTE"

  YAPS.alTE = YAPS.alTE(1:YAPS.lContrasts);

  %%%-------------------------------%
  %%% SPECIAL CASE: regrid parameters (only a pain in VB15)

  if ( isempty(YAPS.alRegridMode) ),

    param_list = {'alRegridRampupTime', 'alRegridRampdownTime', ...
                  'alRegridFlattopTime', 'alRegridDelaySamplesTime', ...
                  'alRegridMode', 'aflRegridADCDuration'};

    for ind = 1:length(param_list),
      param = param_list{ind};

      match = regexp(header, ['<Param\w+\."(?<param>' param ')">\s*{\s*(?<value>[\d ]*)\s*}'], 'names');

      % SPECIAL CASE (within a special case): pesky 'Precision' string for "aflRegridADCDuration"
      if ( isempty(match) ),
        match = regexp(header, ['<Param\w+\."(?<param>' param ')">\s*{\s*<Precision>[ \d]*\s*(?<value>[\d \.]*)\s*}'], 'names');
      end;

      match.value = str2num(match.value);
      YAPS = setfield(YAPS, param_list{ind}, match.value(1));

    end;

  end;

  %%%-------------------------------%
  %%% ENUMERATED TYPE: "alRegridMode"

  switch YAPS.alRegridMode,
    % "enum RegriddingMode" snarfed from <n4/pkg/MrServers/MrProtSrv/MrProt/SeqDefines.h>

   case hex2dec('01'),
    YAPS.ucRegridMode = 'REGRID_NONE';
   case hex2dec('02'),
    YAPS.ucRegridMode = 'REGRID_TRAPEZOIDAL';
   case hex2dec('04'),
    YAPS.ucRegridMode = 'REGRID_SINUSOIDAL';
   otherwise,
    warning(sprintf('unrecognized mode for "%s":  [ %d ]', 'alRegridMode', YAPS.alRegridMode));
  end;

  %%%-------------------------------%
  %%% SPECIAL CASE: flBandwidthPerPixelPhaseEncode and iEffectiveEpiEchoSpacing (only present in VB15)

  param_list = {'flBandwidthPerPixelPhaseEncode', 'iEffectiveEpiEchoSpacing'};

  for ind = 1:length(param_list),
    param = param_list{ind};

    match = regexp(header, ['<Param\w+\."(?<param>' param ')">\s*{\s*(?<value>[\d.?\d]*)\s*}'], 'names');

    if ( isempty(match) ),
      continue;
    end;

    match.value = str2num(match.value);

    if ( isempty(match.value) ),
      match.value = 0;
    end;

    YAPS = setfield(YAPS, param_list{ind}, match.value(1));

  end;

  %%%-------------------------------%
  %%% SPECIAL CASE: sSliceArray

  param_list = {'sPosition.dSag', 'sPosition.dCor', 'sPosition.dTra', ...
                  'sNormal.dSag',   'sNormal.dCor',   'sNormal.dTra', ...
                'dThickness', 'dPhaseFOV', 'dReadoutFOV', 'InPlaneRot'};

  % empty sub-struct for slice array
  asSlice = cell2struct(cell(length(param_list),1), regexprep(param_list, '\.', '_'), 1);

  for slc = 1:YAPS.sSliceArray_lSize,

    % assign new empty sub-struct for each slice
    YAPS.sSliceArray(slc) = asSlice;

    for ind = 1:length(param_list),

    param = param_list{ind};
    param = regexprep(param, '\.', '\\.');

    match = regexp(header, ['sSliceArray\.asSlice\[' num2str(slc-1) '\]\.(?<param>' param ')\s*=\s*(?<value>(-*)\d*\.?x?\w*)'], 'names');

    % check if no match is found (e.g., only one of three possible normals will be represented, so empty match is common)
    if ( isempty(match) ),
      match(1).value = '0';
    end;

    % consider only first match
    match = match(1);

    YAPS.sSliceArray(slc) = setfield(YAPS.sSliceArray(slc), regexprep(param_list{ind}, '\.', '_'), str2num(match.value));
    end;
  end;

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

  if ( nargout > 0 ),
    varargout{1} = YAPS;
    varargout{2} = Config_evp;
    varargout{3} = header;
  end;

  return;

  %************************************************************************%
  %%% $Source: /space/repo/1/dev/dev/matlab/read_meas_prot.m,v $
  %%% Local Variables:
  %%% mode: Matlab
  %%% fill-column: 76
  %%% comment-column: 0
  %%% End:

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function x=irdct(y,n,a,b)
%IRDCT    Inverse discrete cosine transform of real data X=(Y,N)
% Data is truncated/padded to length N.
%
% This routine is equivalent to multiplying by the matrix
%
%    irdct(eye(n)) = cos((0.5:n)'*(0:n-1)*pi/n)*diag([sqrt(0.5)*a/b repmat(a,1,n-1)])/n
%
%
% Default values of the scaling factors are A=sqrt(2N) and B=1 which
% results in an orthogonal matrix. Other common values are A=1 or N and/or B=1 or sqrt(2).
% If b~=1 then the columns are no longer orthogonal.
% See also RCDT

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: irdct.m 713 2011-10-16 14:45:43Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fl=size(y,1)==1;
if fl y=y(:); end
[m,k]=size(y);
if nargin<4 b=1;
    if nargin<3 a=sqrt(2*m);
        if nargin<2 n=m;
        end
    end
end
if n>m y=[y; zeros(n-m,k)];
elseif n<m y(n+1:m,:)=[];
end

x=zeros(n,k);
w=ones(1,k);
m=fix((n+1)/2);
p=n-m;
z=0.5*exp((0.5i*pi/n)*(1:p)).';
u=(y(2:p+1,:)-1i*y(n:-1:m+1,:)).*z(:,w)*a;
y=[y(1,:)*sqrt(0.5)*a/b; u(1:m-1,:)];
if m==p
    z=-0.5i*exp((2i*pi/n)*(0:m-1)).';
    y=(z(:,w)+0.5).*(conj(flipud(u))-y)+y;
    z=ifft(y,[],1);
    u=real(z);
    y=imag(z);
    q=m/2;
    h=rem(m,2)/2;
    x(1:4:n,:)=u(1:q+h,:);
    x(2:4:n,:)=y(m:-1:q+1-h,:);
    x(3:4:n,:)=y(1:q-h,:);
    x(4:4:n,:)=u(m:-1:q+1+h,:);
else
    z=real(ifft([y; conj(flipud(u))]));
    x(1:2:n,:)=z(1:m,:);
    x(2:2:n,:)=z(n:-1:m+1,:);
end

if fl x=x.'; end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
% UIGETFILE2 - same as uigetfile but remember folder location.
%
% Usage: >> uigetfile2(...)
%
% Inputs: Same as uigetfile
%
% Author: Arnaud Delorme & Hilit Serby, Scott Makeig, SCCN, UCSD, 2004
%         Thanks to input from Bas Kortmann
%
% Copyright (C) Arnaud Delorme & Hilit Serby, Scott Makeig, SCCN, UCSD, 2004

% 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.

function varargout = uigetfile2(varargin);
mlock
persistent tmp_fld

if nargin < 1
    help uigetfile2;
    return;
end;

% remember old folder
olddir = pwd;
try
    eeglab_options;
    if option_rememberfolder
        if isempty(tmp_fld)
            tmp_fld = getpref('eeglab','lastdir',olddir);
        end
        cd(tmp_fld);
    end
end

%% Show the open dialog and save the latest directory to the file
% ---------------------------------------------------------------
[varargout{1} varargout{2}] = uigetfile(varargin{:});
try
    if option_rememberfolder
        if ischar(varargout{1})
            tmp_fld = varargout{2};
            setpref('eeglab','lastdir',tmp_fld);
        end
    end
end
cd(olddir)

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [x,w,x_variables,w_variables,aux_variables,F,failure] = robust_classify_variables(F,h,ops,w);

failure = 0;
% Variables before expanding nonlinear operators
initial_variables = unique([depends(F) depends(h)]);

% Uncertain variables
w_variables = getvariables(w);

% Decision variables
x_variables = 1:yalmip('nvars');
x_variables = setdiff(x_variables,w_variables);

% Any fancy modelling goin on here? If so, we need to expand the model, and
% associate the new lifted variables either to the uncertainty set or to
% the set of general variables. There are some classes of variables
% 1. Uncertain variables
% 2. Original decision variables
% 4. Lifted variables used for uncertain variables in uncertainty description
% 5. Lifted variables used for uncertain variables in uncertain model
% 6. Lifted variables used for mixed  variables in uncertain model
% FIX : This code is currently a miserable hack and should be cleaned up
extended = yalmip('extvariables');
aux_variables = [];
if ~isempty(extended) & any(ismember(initial_variables,extended))

    % Remove auxilliary variables (defined by YALMIP)
    aux_variables = intersect(extended,x_variables);
    % These variables are original decision variables
    x_variables = setdiff(x_variables,aux_variables);

    % Original constraint index
    id = getlmiid(F(find(~lifted(F))));

    % This is a tweak to allow epxansion of bilinear terms in robust problems,
    % is expression such as abs(x*w) < 1 for all -1 < w < 1
    ops.expandbilinear = 1;

    % Expand the model to model norms etc
    [F,failure,cause] = expandmodel(F,h,ops,w);
    if failure % Convexity propgation failed
        interfacedata = [];
        recoverdata = [];
        solver = '';
        diagnostic.solvertime = 0;
        diagnostic.problem = 14;
        diagnostic.info = yalmiperror(14,cause);
        x = [];
        w = [];
        x_variables=[];
        w_variables=[];
        F = [];
        return
    end

    % Auxillary variables (introduced by YALMIP). These variables are lifting
    % variables that are used to model advanced constructs suck as norms etc in
    % the nonlinear operator framework
    new_aux_variables = setdiff(unique([depends(F) depends(h)]),initial_variables);
    aux_variables = setdiff(union(aux_variables,new_aux_variables),w_variables);

    % We know try to find the auxillary variables that are used to model
    % functions that only depends on w
    F_lifted = F(find(~ismember(getlmiid(F),id)));
    aux_variables_w = aux_variables;
    for i = 1:length(F_lifted)
        variables = depends(F_lifted(i));
        with_x = any(ismember(variables,x_variables));
        with_w = any(ismember(variables,w_variables));
        if with_x
            aux_variables_w = setdiff(aux_variables_w,intersect(variables,aux_variables));
            x_variables = union(x_variables,intersect(variables,aux_variables));
        end
    end

    aux_variables_x = setdiff(aux_variables,aux_variables_w);
    x_variables = union(aux_variables_x, x_variables);
    w_variables = union(aux_variables_w, w_variables);
end

x_variables = intersect([depends(F) depends(h)],x_variables);
w_variables = intersect([depends(F) depends(h)],w_variables);

x = recover(x_variables);
x = recover(setdiff(depends(x),w_variables));
x_variables = getvariables(x);

% % Could be some of these two cases
% % F = (abs(x) + w < 1) + (norm(w,1) < 0.3);
% % F = (x + abs(w) < 1) + (-2 < w < 2) ;
%
% w_variables = setdiff(w_variables,heh);
% aux_variables = union(aux_variables,heh);

w = recover(w_variables);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%GEOM_LOFT_FACES Loft faces to solid.
%
%   [ SOUT, NEW_TAG ] = GEOM_LOFT_FACES( SIN, TAGS, FACES ) Loft
%   selected FACES from geometry objects with TAGS to a solid geometry
%   object.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function h = conv(f, g)
%CONV   Convolution of DELTAFUN objects.
%   H = CONV(F, G) produces the convolution of DELTAFUN objects F and G:
%                     -
%                    /
%           H(x) =   |    F(t) G(x-t) dt,  x in [a + c, b + d]
%                    /
%   The output H is a cell array of DELTAFUNs or CLASSICFUNS. If the result of
%   the convolution of F and G is zero, an empty cell is returned. The cell
%   array returned is used by higher level convolutions.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Return empty for an empty input:
if ( isempty(f) || isempty(g) )
    h = deltafun();
    return
end

if ( ~isa(f, 'deltafun') )
    deltaMagF = [];
    deltaLocF = [];
    funF = f;
else
    f = simplifyDeltas(f);
    if ( isa(f, 'deltafun') )
        deltaMagF = f.deltaMag;
        deltaLocF = f.deltaLoc;
        funF = f.funPart;
    else
        deltaMagF = [];
        deltaLocF = [];
        funF = f;
    end
end

if ( ~isa(g, 'deltafun') )
    deltaMagG = [];
    deltaLocG = [];
    funG = g;
else
    g = simplifyDeltas(g);
    if ( isa(g, 'deltafun') )
        deltaMagG = g.deltaMag;
        deltaLocG = g.deltaLoc;
        funG = g.funPart;
    else
        deltaMagG = [];
        deltaLocG = [];
        funG = g;
    end
end

% Extract the domains of f and g:
domF = funF.domain;
domG = funG.domain;
a = domF(1);
b = domF(end);
c = domG(1);
d = domG(end);

% Get the threshold for deltafunctions:
pref = chebfunpref();
deltaTol = pref.deltaPrefs.deltaTol;

% Compute the convolution of funParts and append it to the output cell:
h = conv(funF, funG);

% Remove zero funs for simplicity:
zeroIndices = cellfun( @(hk) iszero(hk), h);
h(zeroIndices) = [];

%% Get all the deltafunction contributions
% (f + df) * (g + dg) = f*g + dg * (f + df) + df * (g + dg) - df * dg
%                     = f*g + dg * (f + df/2) + df * (g + dg/2)

% Contributions due to deltafunctions in F:
% df * (g + dg/2):
h = convDeltas(deltaMagF, deltaLocF, g, c, d, deltaTol, h);

% Contributions due to delta functions in G:
% dg * (f + df/2);
h = convDeltas(deltaMagG, deltaLocG, f, a, b, deltaTol, h);

% Check for emptiness:
if ( isempty(h) )
    h = {fun.constructor(0, struct('domain', [a, b]))};
    return
end

%% Make sure that h is a cell array with non-overlapping funs.

% First extract the breakpints from funs forming h:
doms = cellfun(@(hk) hk.domain, h, 'uniformoutput', 0);
breakPoints = sort(unique(horzcat(doms{:})));

% Coalesce breakpoints that are extremely close together.
if ( any(isinf(breakPoints)) )
    tol = 2*eps;  % hscale of functions on unbounded domains is 1.
else
    tol = 2*eps*norm(breakPoints, Inf);
end

% TODO:  Take the average of points which are close together instead of just
% picking one more or less arbitrarily?
closeBreaks = abs(diff(breakPoints)) < tol;
breakPoints(closeBreaks) = [];

% Restrict each FUN to lie within the new breakpoints.
nFuns = length(h);
H = {};
deltaTol = pref.deltaPrefs.proximityTol;
for k = 1:nFuns
    hk = h{k};
    dom = hk.domain;

    % The domain endpoints of each computed FUN must lie _exactly_ within the
    % set of breakpoints or RESTRICT will fail.
    idx1 = find(abs(breakPoints - dom(1)) < tol);
    idx2 = find(abs(breakPoints - dom(2)) < tol);
    hk = changeMap(hk, breakPoints([idx1, idx2]));

    % If we just changed the map of a DELTAFUN, the locations of deltas at the
    % domain endpoints also need to be updated.
    if ( isa(hk, 'deltafun') )
        deltaLoc = hk.deltaLoc;
        domk = hk.domain;
        if( ~isempty(deltaLoc) )
            if ( abs(deltaLoc(1) - domk(1) ) < deltaTol )
                deltaLoc(1) = domk(1);
            end
            if ( abs(deltaLoc(end) - domk(2)) < deltaTol )
                deltaLoc(end) = domk(2);
            end
            hk.deltaLoc = deltaLoc;
        end
    end

    % Do the restriction.
    hk = restrict(hk, breakPoints(idx1:idx2));
    if ( ~isa(hk, 'cell') )
        hk = {hk};
    end
    H = [H, hk]; %#ok<AGROW>
end

% If H has the same number of funs, nothing further needs to be done:
if ( length(H) == length(breakPoints) - 1 )
    return
end

% H has more funs as a result of restrict. Re-initialize h with 0 BNDFUNS:
nFuns = length(breakPoints) - 1;
h = {};
for k = 1:nFuns
    data.domain = breakPoints(k:k+1);
    h = [h, {bndfun(0, data)}]; %#ok<AGROW>
end

% Loop through H and add each fun to the relevant domain piece in h:
for k = 1:length(H)
    hk = H{k};
    dom = hk.domain;
    if ( abs(dom(2) - dom(1)) > tol )
        idx = find(breakPoints == dom(1));
        h{idx} = h{idx} + hk; %#ok<AGROW>
    end
end

end

function h = convDeltas(deltaMagF, deltaLocF, g, c, d, deltaTol, h)
%   H = CONVDELTAS() convolves the function G with deltafunctions described by
%   deltaMagF and deltaLocF. G originally lives on [c, d] and the output is
%   returned in H. This function is called twice so deltafunctions in g are
%   scaled by half to get the correct delta-delta contributions.

%%
% Contributions due to deltafunctions in F:
% df * (g + dg/2):
[m, n] = size(deltaMagF);
% Loop through the delta function matrix:
for i = 1:m
    for j = 1:n
        if ( abs(deltaMagF(i, j)) > deltaTol )
            % Take appropriate derivative, scale and shift the function:
            hij = deltaMagF(i, j) * changeMap(diff(g,i-1), deltaLocF(j) + [c d]);
            % The half below is to make sure that delta-delta interaction
            % is not counted twice:
            if ( isa(hij, 'deltafun') )
                hij.deltaMag = hij.deltaMag/2;
            end

            % If the result is non-zero, append it:
            if ( ~iszero(hij) )
                h = [h, {hij}]; %#ok<AGROW>
            end
        end
    end
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function k = lfmvKernDiagCompute(lfmKern, t)

% LFMVKERNDIAGCOMPUTE Compute diagonal of LFMV kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for velocity- velocity in
% the switching dynamical latent force kernel.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG t : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

% Get length scale out.
sigma2 = 2/lfmKern.inverseWidth;
sigma = sqrt(sigma2);

% Parameters of the kernel
alpha(1) = lfmKern.damper./(2*lfmKern.mass);
omega(1) = sqrt(lfmKern.spring./lfmKern.mass - alpha(1)*alpha(1));

% Precomputations to increase the speed
preExp1 = zeros(length(t),2);
gamma1_p = alpha(1) + j*omega(1);
gamma1_m = alpha(1) - j*omega(1);
preGamma(1) = gamma1_p + gamma1_p;
preGamma(2) = gamma1_p + gamma1_m;
preGamma(3) = gamma1_m + gamma1_p;
preGamma(4) = gamma1_m + gamma1_m;
preConst = 1./preGamma;
preFactors(1) = preConst(2) - preConst(1);
preFactors(2) = preConst(3) - preConst(4);
preFactors(3) = preConst(3) - preConst(1);
preFactors(4) = preConst(2) - preConst(4);
preExp1(:,1) = gamma1_p.*exp(-gamma1_p*t);
preExp1(:,2) = gamma1_m.*exp(-gamma1_m*t);
% Actual computation of the kernel
sk =lfmDiagComputeH3VV(gamma1_p, gamma1_m, sigma2, t, preFactors([1 2]), 1) + ...
    lfmDiagComputeH3VV(gamma1_p, gamma1_m, sigma2, t, preFactors([3 4]), 0) + ...
    lfmDiagComputeH4VV(gamma1_p, gamma1_m, sigma2, t, preGamma([1 2 4 3]), preExp1 ) + ...
    lfmDiagComputeH4VV(gamma1_p, gamma1_m, sigma2, t, preGamma([1 3 4 2]), preExp1 );

if lfmKern.isNormalised
    k0 = lfmKern.sensitivity^2/(8*sqrt(2)*lfmKern.mass^2*omega^2);
else
    k0 = sqrt(pi)*sigma*lfmKern.sensitivity^2/(8*lfmKern.mass^2*omega^2);
end
k = k0*sk;

function h = lfmDiagComputeH3VV(gamma1_p, gamma1_m, sigma2, t, preFactor, mode)

h = preFactor(1)*lfmvvComputeUpsilonDiagVector(gamma1_p, sigma2, t, mode) ...
    + preFactor(2)*lfmvvComputeUpsilonDiagVector(gamma1_m,sigma2, t, mode);

function h = lfmDiagComputeH4VV(gamma1_p, gamma1_m, sigma2, t, ...
    preFactor, preExp)

h =   lfmvpComputeUpsilonVector(gamma1_p,sigma2, t, 0).*( preExp(:,2)/preFactor(2) - preExp(:,1)/preFactor(1)) ...
    + lfmvpComputeUpsilonVector(gamma1_m,sigma2, t, 0).*( preExp(:,1)/preFactor(4) - preExp(:,2)/preFactor(3));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [engine, loglik] = dbn_update_bel1(engine, evidence)
% DBN_UPDATE_BEL1 Update  the initial belief state (bk)
% engine = dbn_update_bel1(engine, evidence)
%
% evidence{i} has the evidence on node i for slice 1

bnet = bnet_from_engine(engine);
ss = length(bnet.intra);
CPDpot = cell(1,ss);
t = 1;
for n=1:ss
  fam = family(bnet.dag, n);
  e = bnet.equiv_class(n, 1);
  CPDpot{n} = convert_to_pot(bnet.CPD{e}, engine.pot_type, fam(:), evidence);
end

onodes = find(~isemptycell(evidence));

[clpot, loglik] = enter_soft_evidence(engine.sub_engine1, engine.clq_ass_to_node1, CPDpot, onodes, engine.pot_type);

C  = length(engine.clusters);
newbel = cell(1,C);
for c=1:C
  k = engine.clq_ass_to_cluster1(c);
  newbel{c} = marginalize_pot(clpot{k}, engine.clusters{c});
end

engine.bel = newbel;
engine.bel_clpot = clpot;
engine.slice1 = 1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% ATANDEMO  This program creates the ArcTangent example in Chapter 2.
%
% C. T. Kelley, November 13, 2002.
%
x=10;
tol=[1.d-6,1.d-6];
%
% Solve the problem with Newton's method.
%
params=[40, 1, 0,0];
[an, errsn, ierrn]=nsold(x, 'fatan', tol, params);
%
% Set the default Jacobian updating strategy.
%
params=[40, 1000, .5, 0];
[ac, errsd, ierrd]=nsold(x, 'fatan', tol, params);
inewt=length(errsn(:,1)); cnewt=0:inewt-1;
idefault=length(errsd(:,1)); cdefault=0:idefault-1;
%
% Plot the residual histories.
%
semilogy(cnewt,errsn(:,1),'-',cdefault,errsd(:,1),'--')
legend('Newton','default')

repo: Aero-matlab


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Many Optimizing Liaisons (MOL) optimization method is a
% simplification of PSO originally by Eberhart et al. (1, 2).
% MOL does not have any attraction to the particle's own best
% known position. It is similar to the "Social Only" PSO
% suggested by Kennedy (3) and was studied more thoroguhly
% by Pedersen et al. (4) who found it to sometimes outperform
% PSO and have behavioural parameters that were easier to tune.
% This is a parallelized version.
% Literature references:
% (1) J. Kennedy and R. Eberhart. Particle swarm optimization.
%     In Proceedings of IEEE International Conference on Neural
%     Networks, volume IV, pages 1942-1948, Perth, Australia, 1995
% (2) Y. Shi and R.C. Eberhart. A modified particle swarm optimizer.
%     In Proceedings of the IEEE International Conference on
%     Evolutionary Computation, pages 69-73, Anchorage, AK, USA, 1998.
% (3) J. Kennedy. The particle swarm: social adaptation of knowledge,
%     In: Proceedings of the IEEE International Conference on
%     Evolutionary Computation, Indianapolis, USA, 1997.
% (4) M.E.H. Pedersen and A.J. Chipperfield. Simplifying particle swarm
%     optimization. Applied Soft Computing, volume 10, pages 618-628, 2010.
% Parameters:
%     problem; name or handle of optimization problem, e.g. @myproblem.
%     data; data-struct, see e.g. the file myproblemdata.m
%     parameters; behavioural parameters for optimizer,
%                 see file molparameters.m
% Returns:
%     bestX; best found position in the search-space.
%     bestFitness; fitness of bestX.
%     evaluations; number of fitness evaluations performed.
function [bestX, bestFitness, evaluations] = molparallel(problem, data, parameters)

    % Copy data contents to local variables for convenience.
    n = data.Dim;
    acceptableFitness = data.AcceptableFitness;
    maxEvaluations = data.MaxEvaluations;
    lowerBound = data.LowerBound;
    upperBound = data.UpperBound;

    % Behavioural parameters for this optimizer.
    s = parameters(1);        % Swarm-size
    omega = parameters(2);    % Inertia weight.
    phiG = parameters(3);     % Swarm's best weight.

    % Initialize the velocity boundaries.
    range = upperBound-lowerBound;
    lowerVelocity = -range;
    upperVelocity = range;

    % Initialize swarm.
    x = initpopulation(s, n, data.LowerInit, data.UpperInit);    % Particle positions.
    v = initpopulation(s, n, lowerVelocity, upperVelocity);      % Velocities.

    % Compute fitness of initial particle positions. (Parallel)
    fitness = zeros(1, s); % Preallocate array for efficiency.
    parfor i=1:s
        fitness(i) = feval(problem, x(i,:), data);
    end

    % Determine best particle.
    [bestFitness, bestIndex] = min(fitness);
    bestX = x(bestIndex,:);

    % Perform optimization iterations until acceptable fitness
    % is achieved or the maximum number of fitness evaluations
    % has been performed.
    evaluations = s; % Fitness evaluations above count as iterations.
    while (evaluations < maxEvaluations) && (bestFitness > acceptableFitness)

        % Update particle velocities and positions. (Non-parallel)
	for i=1:s
            % Pick random weight.
            rG = rand(1, 1);

            % Update velocity for i'th particle.
            v(i,:) = omega * v(i,:) + rG * phiG * (bestX - x(i,:));

            % Bound velocity.
            v(i,:) = bound(v(i,:), lowerVelocity, upperVelocity);

            % Update position for i'th particle.
            x(i,:) = x(i,:) + v(i,:);

            % Bound position to search-space.
            x(i,:) = bound(x(i,:), lowerBound, upperBound);
        end

        % Compute fitness. (Parallel)
        % Only this fitness evaluation is parallelized
        % which makes synchronization easier.
        parfor i=1:s
            fitness(i) = feval(problem, x(i,:), data);
        end

        % Update swarm's best-known position. (Non-parallel)
        [newBestFitness, newBestIndex] = min(fitness);
        if newBestFitness < bestFitness
            % Update swarm's best-known fitness.
            bestFitness = newBestFitness;

            % Update swarm's best-known position.
            % This must be copied because the particles
            % will continue to move in the search-space.
            bestX = x(newBestIndex,:);
        end

        % Increment counter.
        evaluations = evaluations + s;
    end
end

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

author: Sable

repo: mcbench-benchmarks


text:
function [IND_mask]=maskfind(M,IND,MASK_I,MASK_J,MASK_K)

% function [IND_mask]=maskfind(M,IND,I,J,K)
% ------------------------------------------------------------------------
% This function finds the indices of the elements found inside the mask
% defined by MASK_I, MASK_J and MASK_K. The matrix IND_mask is numel(IND) x
% numel(MASK_I) in size. N.B. zeros are used when indices are not found.
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 07/04/2010
% ------------------------------------------------------------------------

%%

%Get image size
siz=size(M);
if size(M,3)==1
    siz(3)=1;
end

%Force row vector input
MASK_I=MASK_I(:)';
MASK_J=MASK_J(:)';
MASK_K=MASK_K(:)';

%Get subscript indices for IND
[Im,Jm,Km]=ind2sub(siz,IND(:));

%Create adjacency subscript indices
Im=Im(:,ones(1,numel(MASK_I)))+MASK_I(ones(numel(Im),1),:);
Jm=Jm(:,ones(1,numel(MASK_J)))+MASK_J(ones(numel(Jm),1),:);
Km=Km(:,ones(1,numel(MASK_K)))+MASK_K(ones(numel(Km),1),:);

%Get logic for valid indices.
L_valid= ~((Im<1 | Im>siz(1)) | (Jm<1 | Jm>siz(2)) | (Km<1 | Km>siz(3)));

%Get linear indiced of adjacenty indices
[ind_m]=sub2ind(siz,Im(L_valid),Jm(L_valid),Km(L_valid));

IND_mask=zeros(size(Im));
IND_mask(L_valid)=ind_m;

%%
% _*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:
%% amfStruct2xml
% Below is a demonstration of the features of the |amfStruct2xml| function

%%
clear; close all; clc;

%% Syntax
% |[domNode]=amfStruct2xml(amf,fileName,optionStruct);|

%% Description
% This function provides the basis for coding AMF files.

%% Example:

%%
% Create example mesh

[F,V]=geoSphere(2,50); %Faces and vertices
c=[1 0 0]; %RGB Color

volume1_Name='sphere';
material1_Name='sphere material';

%%
cFigure;
title('Model to export to AMF')
gpatch(F,V,c);
axisGeom; camlight headlight;
drawnow;

%%
% Define file name
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');
fileName=fullfile(savePath,'tempModel.amf');

%%
% Create AMF structure

% Define AMF units and version
amf.ATTR.unit='mm'; %Set units
amf.ATTR.version='1.1'; %Set AMF version

% Add optional metadata here
amf.metadata{1}.ATTR.type='name';
amf.metadata{1}.VAL='Sphere mesh';

amf.metadata{2}.ATTR.type='author';
amf.metadata{2}.VAL='GIBBON';

% object definition

%-> vertices
amf.object{1}.ATTR.id=1;
for q=1:1:size(V,1) %Loop over all vertices
    amf.object{1}.mesh.vertices.vertex{q}.coordinates.x=V(q,1);
    amf.object{1}.mesh.vertices.vertex{q}.coordinates.y=V(q,2);
    amf.object{1}.mesh.vertices.vertex{q}.coordinates.z=V(q,3);
end

%-> volume -> triangle
amf.object{1}.mesh.volume{1}.ATTR.materialid=1;
amf.object{1}.mesh.volume{1}.metadata{1}.ATTR.type='name';
amf.object{1}.mesh.volume{1}.metadata{1}.VAL=volume1_Name;

for q=1:1:size(F,1) %Loop over all triangles (note zero indexing hence -1)
    amf.object{1}.mesh.volume{1}.triangle{q}.v1=F(q,1)-1;
    amf.object{1}.mesh.volume{1}.triangle{q}.v2=F(q,2)-1;
    amf.object{1}.mesh.volume{1}.triangle{q}.v3=F(q,3)-1;
end

% material
amf.material{1}.ATTR.id=1;
amf.material{1}.metadata.ATTR.type='name';
amf.material{1}.metadata.VAL=material1_Name;
amf.material{1}.color.r=c(1);
amf.material{1}.color.g=c(2);
amf.material{1}.color.b=c(3);

%% Export xml
amfStruct2xml(amf,fileName);

%% zip file
[pathName,fileNameClean,c]=fileparts(fileName);
zipName=fullfile(pathName,[fileNameClean,'.zip']);
zip(zipName,fileName);
movefile(zipName,fileName);


text:
% WAVELET_FACTORY_1D Create wavelet cascade
%
% Usage
%    [Wop, filters] = WAVELET_FACTORY_1D(N)
%
%    [Wop, filters] = WAVELET_FACTORY_1D(N, filt_opt)
%
%    [Wop, filters] = WAVELET_FACTORY_1D(N, filt_opt, scat_opt)
%
% Input
%    N (int): The size of the signals to be transformed.
%    filt_opt (struct): The filter options, same as FILTER_BANK.
%    scat_opt (struct): The scattering options, same as WAVELET_LAYER_1D.
%
% Output
%    Wop (cell): A cell array of wavelet layer transforms needed for the
%       scattering transform.
%    filters (cell): A cell array of the filters used in defining the
%       wavelets.
%
% Description
%    In order to calculate the scattering coefficients of a signal, a set of
%    filter banks need to be defined first. Given these filter banks, wavelet
%    transforms can be defined on scattering layers. To obtain the former,
%    WAVELET_FACTORY_1D calls the FILTER_BANK function using the parameters
%    provided in filt_opt, the result of which is returned as filters.
%    Each filter bank is then used to create a layer operator, using the
%    function WAVELET_LAYER_1D. This function takes a layer, a filter bank
%    and a set of parameters as input. WAVELET_FACTORY_1D fixes the filter
%    bank as to the element of filters corresponding to the layer order, while
%    scat_opt is used as the parameters. The result is a cell array Wop,
%    of layer operators which take one layer as input, and return two layers
%    as output, A and V corresponding to the "average" and "variation" of the
%    input layer.
%
% See also
%    WAVELET_LAYER_1D, WAVELET_1D, FILTER_BANK

function [Wop, filters] = wavelet_factory_1d(N, filt_opt, scat_opt)
	if nargin < 2
		filters = filter_bank(N);
    else
        filters = filter_bank(N, filt_opt);
    end

	if nargin < 3
		scat_opt = struct();
    end
    scat_opt = fill_struct(scat_opt, 'M', 2); % M is the scattering order

    Wop = cell(1,scat_opt.M);
	for m = 0:scat_opt.M
		filt_ind = min(numel(filters), 1+m);
		Wop{1+m} = @(X)(wavelet_layer_1d(X, filters{filt_ind}, scat_opt));
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [Q2,Bg2,delta_X2,Pk2]=ESKF(Q1,Bg1,delta_X1,ImuData,t,Vm,Pk1)

% ESKF( Error State Kalman Filter)
% derivation <Quaternion Kinematics for the Error-state Kalman Filter>
% author Zhang xin

if isempty(Pk1)
    Para_delta_theta0=0.00001;
    Para_delta_Bg0=0.00001;
    Pk1=diag([[1,1,1]*Para_delta_theta0,[1,1,1]*Para_delta_Bg0]);
end

norm_a=norm(ImuData(1,2:4));
%norm_g=norm(ImuData(1,5:7)-Bg1);
wx=(ImuData(1,5)-Bg1(1))*t;
wy=(ImuData(1,6)-Bg1(2))*t;
wz=(ImuData(1,7)-Bg1(3))*t;

Qp=[ 1    , -wx/2 , -wy/2 , -wz/2  ;...
     wx/2 ,   1   ,  wz/2 , -wy/2  ;...
     wy/2 , -wz/2 ,   1   ,  wx/2  ;...
     wz/2 ,  wy/2 , -wx/2 ,   1   ]*Q1';
Qp=Qp/norm(Qp);

delta_theta=norm([wx,wy,wz]);

if delta_theta==0
    R_u_delta_theta=eye(3);
else
    u= [wx,wy,wz]'/delta_theta;
    R_u_delta_theta=eye(3)-Skew_symmetric(u)*sin(delta_theta)+Skew_symmetric(u)*Skew_symmetric(u)*(1-cos(delta_theta));
end

F_delta_X=[R_u_delta_theta, -eye(3)*t;...
           zeros(3,3),       eye(3) ];
Para_delta_theta=0.000001;
Para_delta_Bg=0.000001;
Q_delta_X=diag([[1,1,1]*Para_delta_theta,[1,1,1]*Para_delta_Bg]);

delta_Xp = F_delta_X*delta_X1';

P_k=F_delta_X*Pk1*F_delta_X'+Q_delta_X;

if abs(norm_a-9.8)<2 %&& norm_g< 2

   % R1=[2*Qp(1)^2+2*Qp(2)^2-1;...
   %     2*Qp(2)*Qp(3)-2*Qp(1)*Qp(4);...
   %     2*Qp(1)*Qp(3)+2*Qp(2)*Qp(4)];

    R2=[2*Qp(2)*Qp(3)+2*Qp(1)*Qp(4);...
        2*Qp(1)^2+2*Qp(3)^2-1;...
        2*Qp(3)*Qp(4)-2*Qp(1)*Qp(2)];

    R3=[2*Qp(2)*Qp(4)-2*Qp(1)*Qp(3);...
        2*Qp(3)*Qp(4)+2*Qp(1)*Qp(2);...
        2*Qp(1)^2+2*Qp(4)^2-1 ];

   % J1=2*[ Qp(1), Qp(2),-Qp(3),-Qp(4);...
   %       -Qp(4), Qp(3), Qp(2), Qp(1);...
   %        Qp(3), Qp(4), Qp(1), Qp(2)];

    J2=2*[ Qp(4), Qp(3), Qp(2) , Qp(1);...
           Qp(1),-Qp(2), Qp(3) ,-Qp(4);...
          -Qp(2),-Qp(1), Qp(4) , Qp(3)] ;

    J3=2*[-Qp(3), Qp(4),-Qp(1) , Qp(2);...
           Qp(2), Qp(1), Qp(4) , Qp(3);...
           Qp(1),-Qp(2),-Qp(3) , Qp(4)] ;

    h_acc= R3;
    h_mag = Vm(2)*R2+Vm(3)*R3;

    H_acc = [J3,zeros(3,3)];
    H_mag = [Vm(2)*J2+Vm(3)*J3,zeros(3,3)];

    Q_delta_theta= 1/2*[ -Qp(2), -Qp(3), -Qp(4) ;...
                          Qp(1), -Qp(4),  Qp(3) ;...
                          Qp(4),  Qp(1), -Qp(2) ;...
                         -Qp(3),  Qp(2),  Qp(1) ];

    X_delta_x=[Q_delta_theta, zeros(4,3);...
                zeros(3,3),      eye(3)];

    Hk= [H_acc;H_mag]* X_delta_x;

    para_Rk_acc=0.1;
    para_Rk_mag=0.2;
    Rk=diag([[1,1,1]*para_Rk_acc,[1,1,1]*para_Rk_mag]);

    Kk=P_k*Hk'*inv(Hk*P_k*Hk'+Rk);

    delta_X_hat=Kk*([ImuData(1,2:4)'/norm_a;ImuData(1,8:10)'/norm(ImuData(1,8:10))]-[h_acc;h_mag]);

    Q2=quaternProd(Qp, [1;delta_X_hat(1:3)/2])';
    Q2=Q2/norm(Q2);
    Bg2=Bg1+delta_X_hat(4:6)';

    Pk2_=(eye(6)-Kk*Hk)*P_k;

    delta_X2=delta_Xp'+delta_X_hat';

    G=[eye(3)-Skew_symmetric(-delta_X2(1:3)/2),zeros(3,3);...
        zeros(3,3)   ,                 eye(3)];

    Pk2=G*Pk2_*G';

else
    Q2=Qp';
    Bg2=Bg1;
    delta_X2=delta_Xp';
    Pk2=P_k;

end
    delta_X2(1:3)=[0,0,0];

    if Q2(1)<0
        Q2=-Q2;
    end

end

function S=Skew_symmetric(u)
%Skew_Operator

S=[  0  , -u(3) ,  u(2) ;...
    u(3),   0   , -u(1) ;...
   -u(2),  u(1),    0 ];

end

function ab = quaternProd(a, b)
    ab(1,1) = a(1)*b(1)-a(2)*b(2)-a(3)*b(3)-a(4)*b(4);
    ab(2,1) = a(1)*b(2)+a(2)*b(1)+a(3)*b(4)-a(4)*b(3);
    ab(3,1) = a(1)*b(3)-a(2)*b(4)+a(3)*b(1)+a(4)*b(2);
    ab(4,1) = a(1)*b(4)+a(2)*b(3)-a(3)*b(2)+a(4)*b(1);
    if ab(1)<0
        ab=-ab;
    end
end

author: shenshikexmu

repo: IMUCalibration-Gesture


text:
%% RRT for Planar Redundant Manipulator
%
%

% Do not change anything in rrt.m,
% rrt_star.m and rrt_star_fn.m
%

%

%% Getting started
% Create *main_redundant_rrt_star.m* file. You can create m-file with any other name, that is
% perfectly okay, this will not affect to the solution of a path/motion
% planning problem. All the sources could be found in examples/ directory
% of the distribution.
%

%% Step 1: Choosing the map
% Firstly, we should define what map to use. It is done by defining map
% structure the sample code is goes after.
%
%   map = struct('name', 'bench_redundant_3.mat', 'start_point', [0 0], 'goal_point', [35 35]);

%%
%
% * *name field* defines the file name in maps/ directory
% * *start_point* field defines where the initial point of the problem
% * *goal_point* field define where is the goal point on the map
%

%% Step 2: Setting maximum number of iterations
%
%   max_iter = 20e3;
%
% * *max_iter* variable defines how much iteration should be done to solve
% the path planning problem.

%% Step 3: Do we have to benchmark
%
%   is_benchmark = false;
%
% * *is_benchmark* enables benchmarking. For more details please read the
% sources of rrt.m, rrt_star.m and rrt_star_fn.m

%% Setting random seed
%
%   rand_seed = 40;
%
% * *rand_seed* variable is a random seed for random number generator, it
% is used for sampling nodes. It is usually used for benchmarking. We set
% the same map, however the random seed is different for the same map. You
% can use *now* if you don't care about random seed.
%
%   rand_seed = now;
%

%% Choosing the class (model) we want
%
%   variant = 'FNRedundantManipulator';
%
% * *variant* defines from what class we should instantiate the object. In
% other words it defines what model we choose for application of RRT.
%
% *FNSimple2D* is a name of a class which contains all necessary methods
% and fields in order to represent simple 2D Mobile Robot model.

%% RRT*
%
%   rrt_star(map, max_iter, is_benchmark, rand_seed, variant);
%
% Line above runs RRT with given parameters. In addition, *rrt* function
% returns the class object with a certain solution.
%

%% Sources of *main_redundant_rrt_star.m*
% Press
% <matlab:edit('examples/main_redundant_rrt_star.m') here>
% to play with example code.
%
%   % RRT for Redundant manipulator example.
%   % by Olzhas Adiyatov
%   % 08/28/2013
%
%   map = struct('name', 'bench_redundant_3.mat', 'start_point', [0 0], 'goal_point', [35 35]);
%   max_iter = 20e3;
%   is_benchmark = false;
%   rand_seed = 40;
%   variant = 'FNRedundantManipulator';
%   result = rrt_star(map, max_iter, is_benchmark, rand_seed, variant);
%

author: olzhas


text:
classdef nnscale < nntest
  methods (Test)

    function basic(test)
      batchSize = 10 ;
      x1 = test.randn([5 5 3 batchSize]) ;
      x2 = test.randn([1 1 3 1]) ; % match along non-singletons
      b = [] ;
      y = vl_nnscale(x1, x2, b) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      [dzdx1, dzdx2, dzdb] = vl_nnscale(x1, x2, b, dzdy) ;
      test.der(@(x1) vl_nnscale(x1, x2, b), x1, dzdy, dzdx1, 1e-3*test.range) ;
      test.der(@(x2) vl_nnscale(x1, x2, b), x2, dzdy, dzdx2, 1e-3*test.range) ;
      assert(isempty(dzdb), 'bias derivative should be empty') ;
    end

    function basicBias(test)
      batchSize = 10 ;
      x1 = test.randn([5 5 3 batchSize]) ;
      x2 = test.randn([1 1 3 1]) ; % match along non-singletons
      b = test.randn([1 1 3 1]) ; % match along non-singletons
      y = vl_nnscale(x1, x2, b) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      [dzdx1, dzdx2, dzdb] = vl_nnscale(x1, x2, b, dzdy) ;
      test.der(@(x1) vl_nnscale(x1, x2, b), x1, dzdy, dzdx1, 1e-3*test.range) ;
      test.der(@(x2) vl_nnscale(x1, x2, b), x2, dzdy, dzdx2, 1e-3*test.range) ;
      test.der(@(b) vl_nnscale(x1, x2, b), b, dzdy, dzdb, 1e-3*test.range) ;
    end

  end
end

repo: MFT


text:
function Loglike = Loglike_Basis( Seqs, model, alg )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Learning Hawkes processes via maximum likelihood estimation
% Different regularizers (low-rank, sparse, group sparse) of parameters and
% their combinations are considered, which are solved via ADMM.
%
% Reference:
% Xu, Hongteng, Mehrdad Farajtabar, and Hongyuan Zha.
% "Learning Granger Causality for Hawkes Processes."
% International Conference on Machine Learning (ICML). 2016.
%
% Provider:
% Hongteng Xu @ Georgia Tech
% June. 10, 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial
Aest = model.A;
muest = model.mu;

%D = size(Aest, 1);

tic;

Loglike = 0; % negative log-likelihood

% E-step: evaluate the responsibility using the current parameters
for c = 1:length(Seqs)
    Time = Seqs(c).Time;
    Event = Seqs(c).Mark;
    Tstart = Seqs(c).Start;

    if isempty(alg.Tmax)
        Tstop = Seqs(c).Stop;
    else
        Tstop = alg.Tmax;
        indt = Time < alg.Tmax;
        Time = Time(indt);
        Event = Event(indt);
    end

    %Amu = Amu + Tstop - Tstart;

    dT = Tstop - Time;
    GK = Kernel_Integration(dT, model);

    Nc = length(Time);

    for i = 1:Nc

        ui = Event(i);

        ti = Time(i);

        lambdai = muest(ui);
        %pii = muest(ui);
        %pij = [];

        if i>1

            tj = Time(1:i-1);
            uj = Event(1:i-1);

            dt = ti - tj;
            gij = Kernel(dt, model);
            auiuj = Aest(uj, :, ui);
            pij = auiuj .* gij;
            lambdai = lambdai + sum(pij(:));
        end

        Loglike = Loglike - log(lambdai);

    end

    Loglike = Loglike + (Tstop-Tstart).*sum(muest);
    Loglike = Loglike + sum( sum( GK.*sum(Aest(Event,:,:),3) ) );

end

Loglike = -Loglike;

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function i4vec_indexed_heap_d_max_test ( )

%*****************************************************************************80
%
%% I4VEC_INDEXED_HEAP_D_MAX_TEST tests I4VEC_INDEXED_HEAP_D_MAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
  m = 20;
  n_max = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_INDEXED_HEAP_D_MAX_TEST\n' );
  fprintf ( 1, '  For an indexed I4VEC,\n' );
  fprintf ( 1, '  I4VEC_INDEXED_HEAP_D_MAX reports the maximum value.\n' );
%
%  Set the data array.  To keep things easy, we will use the indicator vector.
%
  a = i4vec_indicator1 ( m );
%
%  The index array will initially be a random subset of the numbers 1 to M,
%  in random order.
%
  n = 5;
  indx(1:11,1) = [ 9, 2, 8, 14, 5, 7, 15, 1, 19, 20, 3 ]';

  i4vec_print ( m, a, '  The data vector:' );
  i4vec_print ( n, indx, '  The index vector:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(INDX):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %4d\n', i, a(indx(i)) );
  end
%
%  Create a descending heap from the indexed array.
%
  indx = i4vec_indexed_heap_d ( n, a, indx );

  i4vec_print ( n, indx, '  The index vector after heaping:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(INDX) after heaping:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %4d\n', i, a(indx(i)) );
  end
%
%  Insert five entries, and monitor the maximum.
%
  for i = 1 : 5

    indx_insert = indx(n+1);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Inserting value %d\n', a(indx_insert) );

    [ n, indx ] = i4vec_indexed_heap_d_insert ( n, a, indx, indx_insert );

    indx_max = i4vec_indexed_heap_d_max ( n, a, indx );

    fprintf ( 1, '  Current maximum is %d\n', a(indx_max) );

  end
  i4vec_print ( m, a, '  The data vector after insertions:' )
  i4vec_print ( n, indx, '  The index vector after insertions:' )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(INDX) after insertions:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %4d\n', i, a(indx(i)) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym eval (@var{f})
%% Symbolic expression to double, taking values from workspace.
%%
%% For expressions without symbols, @code{eval} does the same thing
%% as @code{double}:
%% @example
%% @group
%% f = 2*sin(sym(3))
%%   @result{} f = (sym) 2⋅sin(3)
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% eval(f)
%%   @result{} ans = 0.2822
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% double(f)
%%   @result{} ans = 0.2822
%% @end group
%% @end example
%%
%% For an expression containing symbols, @code{eval} looks in the
%% workspace for variables whose names match the symbols.  It then
%% evaluates the expression using the values from those variables.
%% For example:
%% @example
%% @group
%% syms x y
%% f = x*sin(y)
%%   @result{} f = (sym) x⋅sin(y)
%% @end group
%%
%% @group
%% x = 2.1
%%   @result{} x = 2.1000
%% y = 2.9
%%   @result{} y = 2.9000
%% @end group
%%
%% @group
%% f
%%   @result{} f = (sym) x⋅sin(y)
%%
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% eval(f)
%%   @result{} ans = 0.5024
%% @end group
%% @end example
%%
%% @seealso{@@sym/subs}
%% @end defmethod

function g = eval(f)

  if (nargin ~= 1)
    print_usage ();
  end

  %% take values of x from the workspace
  in = findsymbols (f);
  out = {};
  i = 1;
  while (i <= length (in))
    xstr = char (in{i});
    try
      xval = evalin ('caller', xstr);
      foundit = true;
    catch
      foundit = false;
    end
    if (foundit)
      out{i} = xval;
      i = i + 1;
    else
      in(i) = [];  % erase that input
    end
  end

  try
    %% Fails if the workspace doesn't have values for all symbols.
    % Could also fail for fcns with broken "roundtrip"
    fh = function_handle(f, 'vars', in);
    g = fh(out{:});
    return
  catch
    % no-op
  end

  %% Instead, try substituting and converting to double.
  g = subs (f, in, out);
  try
    g = double (g);
  catch
    % just g then
  end

end

%!error eval (sym(1), 2)

%!assert (isnumeric (eval (sym(3))))
%!assert (isnumeric (eval (sin (sym(3)))))

%!test
%! syms x y
%! f = 2*x*y;
%! x = 3;
%! y = 4;
%! g = eval (f);
%! assert (isequal (g, 24))

%!test
%! syms x y
%! f = 2*x*y;
%! clear y
%! x = 3;
%! g = eval (f);
%! assert (isequal (g, 6*sym('y')))

%!test
%! % do not convert inputs to sym, for SMT compat
%! nearpi = pi + 1e-14;  % sym could make this pi
%! x = sym('x');
%! f = 2*x;
%! x = nearpi;
%! d = eval (f);
%! assert (abs (d - 2*pi) > 1e-15)

repo: octsympy


text:
function [s, t, tobj] = cs_must_compile (srcdir, f, suffix, obj, hfile, force)
%CS_MUST_COMPILE return 1 if source code f must be compiled, 0 otherwise
%   Used by cs_make, and MATLAB/Test/cs_test_make.m.
%
%   Example:
%        none, not meant for end users.
%   See also: CS_MAKE.

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

dobj = dir ([f obj]) ;
if (force | isempty (dobj))                                                 %#ok
    s = 1 ;
    t = Inf ;
    tobj = -1 ;
    return
end
dsrc = dir ([srcdir f suffix '.c']) ;
dh = dir (hfile) ;
t = max (datenum (dsrc.date), datenum (dh.date)) ;
tobj = datenum (dobj.date) ;
s = (tobj < t) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [model,output] = normalizeExponentialCone(model)

% N.B: YALMIP Definition % x2*exp(x1/x2)  <= x3

% Temporarily extract the data in a more standard name min c'*y, b + Ay >= 0
output.problem = 0;
if ~isempty(model.evalMap)

    % Temporarily extract the data in a more standard name min c'*y, b + Ay >= 0
    data.A = model.F_struc(:,2:end);
    data.b = full(model.F_struc(:,1));
    data.c = full(model.c);
    cones = model.K;
    % Clear away the SDP cone to easier add new rows. Will be put back in
    % place later
    data.A(end-sum(model.K.s.^2)+1:end,:)=[];
    data.b(end-sum(model.K.s.^2)+1:end) = [];
    sdpData = model.F_struc(end-sum(model.K.s.^2)+1:end,:);

    % First check that we have only exponentials
    convexFunctions = [];
    concaveFunctions = [];
    vectorentropy = [];
    vectorkullback = [];
    vectorlogsumexp = [];
    for i = 1:length(model.evalMap)
        switch model.evalMap{i}.fcn
            case {'exp','pexp'}
                convexFunctions = [convexFunctions model.evalMap{i}.computes];
            case {'log','plog','slog'}
                concaveFunctions = [concaveFunctions model.evalMap{i}.computes];
            case 'entropy'
                concaveFunctions = [concaveFunctions model.evalMap{i}.computes];
                if length(model.evalMap{i}.variableIndex) > 1
                    vectorentropy = [vectorentropy i];
                end
            case 'kullbackleibler'
                convexFunctions = [convexFunctions model.evalMap{i}.computes];
                if length(model.evalMap{i}.variableIndex) > 2
                    vectorkullback = [vectorkullback i];
                end
            case 'logsumexp'
                convexFunctions = [convexFunctions model.evalMap{i}.computes];
                vectorlogsumexp = [vectorlogsumexp i];
            otherwise
                % Standard interface, return solver not applicable
                % This should not be able to happen as we check this
                % earlier
                output = createoutput([],[],[],-4,model.solver.tag,[],[],0);
                return
        end
    end
    % Check that all exp/log enter in a convex fashion
    if model.K.f > 0
        if nnz(data.A(1:model.K.f,convexFunctions))>0 || nnz(data.A(1:model.K.f,concaveFunctions))>0
            output = createoutput([],[],[],19,model.solver.tag,[],[],0);
            return
        end
    end
    % Check sign in objective on exp/log terms
    if any(data.c(convexFunctions) < 0) || any(data.c(concaveFunctions) > 0)
        output = createoutput([],[],[],19,model.solver.tag,[],[],0);
        return
    end
    % Check sign in elementwise inequalities on exp/log terms
    if model.K.l > 0
        if nnz(data.A(1+model.K.f:model.K.f+model.K.l,convexFunctions)>0) || nnz(data.A(1+model.K.f:model.K.f+model.K.l,concaveFunctions)<0)
            output = createoutput([],[],[],19,model.solver.tag,[],[],0);
            return
        end
    end
    % Check so there are no exp/log terms in quadratic or SDP cone
    if sum(model.K.q) + sum(model.K.s) > 0
        if nnz(data.A(1+model.K.f+model.K.l:end,convexFunctions))>0 || nnz(data.A(1+model.K.f+model.K.l:end,concaveFunctions))>0
            output = createoutput([],[],[],19,model.solver.tag,[],[],0);
            return
        end
    end
    % Scalarize entropy operator
    if ~isempty(vectorentropy)
        for i = vectorentropy(:)'
            involved = model.evalMap{i}.variableIndex;
            computes = model.evalMap{i}.computes;
            n = length(data.c);
            m = length(involved);
            data.A = [sparse(ones(1,m+1),[computes n+(1:m)],[-1 ones(1,m)],1,n+m);
                      data.A spalloc(size(data.A,1),m,0)];
            data.b = [0;data.b];
            data.c = [data.c;zeros(m,1)];
            cones.f = cones.f + 1;
            % The original strucuture now just relates to the first new term
            model.evalMap{i}.computes = n+1;
            model.evalMap{i}.variableIndex = involved(1);
            % and then we add some
            for j = 2:length(involved)
                model.evalMap{end+1} = model.evalMap{i};
                model.evalMap{end}.variableIndex = involved(j);
                model.evalMap{end}.computes = n+j;
            end
        end
    end

    % Scalarize kullbackleibler operator
    if ~isempty(vectorkullback)
        for i = vectorkullback(:)'
            involved = model.evalMap{i}.variableIndex;
            computes = model.evalMap{i}.computes;
            n = length(data.c);
            m = length(involved)/2;
            data.A = [sparse(ones(1,m+1),[computes n+(1:m)],[-1 ones(1,m)],1,n+m);
                      data.A spalloc(size(data.A,1),m,0)];
            data.b = [0;data.b];
            data.c = [data.c;zeros(m,1)];
            cones.f = cones.f + 1;
            % The original strucuture now just relates to the first new term
            model.evalMap{i}.computes = n+1;
            model.evalMap{i}.variableIndex = involved([1 m+1]);
            % and then we add some
            for j = 2:length(involved)/2
                model.evalMap{end+1} = model.evalMap{i};
                model.evalMap{end}.variableIndex = involved([j j+m]);
                model.evalMap{end}.computes = n+j;
            end
        end
    end

    % Scalarize logsumexp operator
    % write log(expx1+expx2..)<= z as exp(x1-z)+exp(x2-z)+...<=1
    if ~isempty(vectorlogsumexp)
        for i = vectorlogsumexp(:)'
            involved = model.evalMap{i}.variableIndex;
            computes = model.evalMap{i}.computes;
            % Orignal #variables
            n = length(data.c);
            % Number of terms in logsumexp
            m = length(involved);

            % exp(x1-z)+exp(x2-z)+...<=1
            data.A = [data.A(1:cones.f,:) spalloc(cones.f,m,0);
                      spalloc(1,n,0) -ones(1,m);
                      data.A(cones.f+1:end,:) spalloc(cones.l+sum(cones.q)+sum(cones.s),m,0)];
            data.b = [data.b(1:cones.f);
                      1;
                      data.b(cones.f+1:end)];
            data.c = [data.c;zeros(m,1)];
            cones.l = cones.l + 1;

            % The original strucuture now just relates to the first new term
            model.evalMap{i}.computes = n+1;
            model.evalMap{i}.variableIndex = [involved(1) computes];
            model.evalMap{i}.fcn = 'expdiff';
            % and then we add some new
            for j = 2:length(involved)
                model.evalMap{end+1} = model.evalMap{i};
                model.evalMap{end}.variableIndex = [involved(j) computes];
                model.evalMap{end}.computes = n+j;
            end
        end
    end

    % Describe all exponential cones
    m = length(model.evalMap);
    cones.e = model.K.e + m;
    for i = 1:m
        switch model.evalMap{i}.fcn
            case 'exp'
                % 1*exp(xv/1) <= xc
                % y*exp(x/y)  <= z.  y new variable, xv the original variable, and xc the "computed"
                x = model.evalMap{i}.variableIndex;
                z = model.evalMap{i}.computes;
                data.A = [data.A;
                          -sparse([1;3],[x z],-1,3,size(data.A,2))];
                data.b = [data.b;[0;1;0]];
            case 'pexp'
                % xv(1)*exp(xv(2)/xv(1)) <= xc
                x = model.evalMap{i}.variableIndex(2);
                y = model.evalMap{i}.variableIndex(1);
                z = model.evalMap{i}.computes;
                data.A = [data.A;
                          -sparse([1;2;3],[x y z],[-1 -1 -1],3,size(data.A,2))];
                data.b = [data.b;[0;0;0]];
            case 'log'
                % log(xv) >= xc i.e. xv >= exp(xc/1)*1
                z = model.evalMap{i}.variableIndex;
                x = model.evalMap{i}.computes;
                data.A = [data.A;
                          -sparse([1;3],[x z],-1,3,size(data.A,2))];
                data.b = [data.b;[0;1;0]];
            case 'plog'
                % xv(1)log(xv(2)/xv(1))>=xc i.e.
                % -xc >= -xv(1)log(xv(2)/xv(1))
                z = model.evalMap{i}.computes;
                y = model.evalMap{i}.variableIndex(1);
                x = model.evalMap{i}.variableIndex(2);
                data.A = [data.A;
                          -sparse([1;2;3],[x y z],[1 -1 1],3,size(data.A,2))];
                data.b = [data.b;[0;0;0]];
            case 'slog'
                % log(1+xv) >= xc i.e. (1+xv) >= exp(xc/1)*1
                z = model.evalMap{i}.variableIndex;
                x = model.evalMap{i}.computes;
                data.A = [data.A;
                          -sparse([1;3],[x z],-1,3,size(data.A,2))];
                data.b = [data.b;0;1;1];
            case 'entropy'
                % -xv*log(xv)>=xc i.e. 1 >= exp(xc/xv)*xv
                x = model.evalMap{i}.computes;
                y = model.evalMap{i}.variableIndex;
                data.A = [data.A;
                          -sparse([1;2],[x y],[-1 -1],3,size(data.A,2))];
                data.b = [data.b;[0;0;1]];
            case 'kullbackleibler'
                % -xv1*log(xv1/xv2)>=-xc i.e. xv2 >= exp(-xc/xv1)*xv1
                x = model.evalMap{i}.computes;
                y = model.evalMap{i}.variableIndex(1);
                z = model.evalMap{i}.variableIndex(2);
                data.A = [data.A;
                          -sparse([1;2;3],[x y z],[1 -1 -1],3,size(data.A,2))];
                data.b = [data.b;[0;0;0]];
            case 'expdiff'
                % exp(xv(1)-xv(2)) <= xc
                x = model.evalMap{i}.variableIndex;
                z = model.evalMap{i}.computes;
                data.A = [data.A;
                          -sparse([1;1;3],[x(1) x(2) z],[-1 1 -1],3,size(data.A,2))];
                data.b = [data.b;[0;1;0]];
            otherwise
        end
    end
    model.F_struc = [data.b data.A];
    model.c = data.c;
    model.K = cones;
    % Put back the SDP cone in place
    if model.K.s(1)>0
        if size(sdpData,2) < size(model.F_struc,2)
            sdpDAta(end,size(model.F_struc,2)) = 0;
        end
        model.F_struc = [model.F_struc;sdpData];
    end
end
n = size(model.F_struc,2)-1;
if n > length(model.lb)
    missing = n - length(model.lb);
    model.lb = [model.lb;-inf(missing,1)];
end
if size(model.F_struc,2)-1 > length(model.ub)
    missing = n - length(model.ub);
    model.ub = [model.ub;inf(missing,1)];
end


text:
function jed = ymdf_to_jed_gregorian ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_GREGORIAN converts a Gregorian YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, real F, the YMDF date.
%
%    Output, real JED, the corresponding JED.
%

%
%  Check the date.
%
  [ y, m, d, ierror ] = ymd_check_gregorian ( y, m, d );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Account for the missing year 0 by moving negative years up one.
%
  y2 = y_common_to_astronomical ( y );
%
%  Convert the calendar date to a computational date.
%
%  This frap-dapping formula fails for 1 March 1900!
%
  if ( m == 3 && d == 1 && ~year_is_leap_gregorian ( y2 ) )
    m = m - 1;
    d = 29;
  end

  y_prime = y2 + 4716 - floor ( ( 14 - m ) / 12 );
  m_prime = mod ( m + 9, 12 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  j1 = floor ( ( 1461 * y_prime ) / 4 );

  j2 = floor ( ( 153 * m_prime + 2 ) / 5 );

  g = floor ( 3 * floor ( ( y_prime + 184 ) / 100 ) / 4 ) - 38;

  jed = j1 + j2 + d_prime - 1401 - g - 0.5;
  jed = jed + f;

  return
end


text:
function ax_out = subplot_tight(numRows,numCols,currIdx,spacing)
% tight version of subplot, by xiuye
% FUNCTION: subplot_tight(numRow,numCol,ax_curr,spacing)
%
% INPUT:
% - numRows: number of rows
% - numCols: number of columns
% - currIdx: the index of the current selected suplot, counting
%     left to right and top to bottom, like
%     1 | 2 | 3
%     4 | 5 | 6
%     This parameter can also take a vector of several indices. The user
%     would typically specify contiguous indices, and the corresponding
%     subplots will be merged into a single axis (subplot).
% - spacing: optional vector to specify spacing parameters. Can be of length
%     1, 2, 3 or 4.
%     The order of the parameters is [frameWidth,axBorderWidth,bottomframe
%
% OUTPUT:
% - ax_out: axis handle of the current subplot axis, if output requested
%
% USAGE:
% ax_out = subplot_tight(m,n,p), as in the built-in function SUBPLOT, breaks
%     the Figure window into an m-by-n matrix of small axes, selects the
%     p-th axes for the current plot, and returns the axes handle.  The axes
%     are counted along the top row of the Figure window, then the second
%     row, etc.
%     ax_out, returns the current axis if specified [i.e. if no semicolon
%     is used, the axis info is not printed automatically].
%
% subplot_tight(m,n,P), where P is a vector, specifies an axes position
%     that covers all the subplot positions listed in P.
%
% subplot_tight(m,n,p,frameWidth), where frameWidth is a number between 0
%     and 1 (normalized to size of the full figure)that specifies all 4
%     margins of the figure.
%
% subplot_tight(m,n,p,[frameWidth,axBorderWidth]), where axBorderWidth is
%     a number between 0 and 1 (normalized to size of the full figure) that
%     specifies the margins between adjacent subplots (in all 4 directions).
%
% subplot_tight(m,n,p,[frameWidth,axBorderWidth,bottomFrameWidth]), where
%     bottomFrameWidth is a number between 0 and 1 that overwrites the
%     frameWidth value for the bottom margin of the figure.
%
% subplot_tight(m,n,p,[fWidth,aWidth,bottomFrameWidth,topFrameWidth]), where
%     topFrameWidth is a number between 0 and 1 that overwrites the
%     frameWidth value for the top margin of the figure.

%% set params
% set default (before user override)
frameWidth = 0.05;
bottomFrameWidth = frameWidth;
topFrameWidth = frameWidth;

axBorderWidth = 0.02;

% user override:
if exist('spacing','var')
    frameWidth = spacing(1);
    bottomFrameWidth = frameWidth;
    topFrameWidth = frameWidth;

    if length(spacing)>=2,
        axBorderWidth = spacing(2);
    end
    if length(spacing)>=3
        bottomFrameWidth = spacing(3);
    end
    if length(spacing)>=4
        topFrameWidth = spacing(4);
    end
end

leftFrameWidth = frameWidth;
rightFrameWidth = frameWidth;

widthCol = (1-leftFrameWidth-rightFrameWidth)/numCols;
heightRow = (1-topFrameWidth-bottomFrameWidth)/numRows;
netwidthCol = (1-leftFrameWidth-rightFrameWidth-(numCols-1)*axBorderWidth)/numCols;
netheightRow = (1-topFrameWidth-bottomFrameWidth-(numRows-1)*axBorderWidth)/numRows;

[i_col, i_row] = ind2sub([numCols,numRows],currIdx);

if length(currIdx)==1 % default
    % left, bottom, width, height
    pos1 = leftFrameWidth+(i_col-1)*widthCol;
    pos2 = 1-topFrameWidth-i_row*heightRow;
    pos3 = netwidthCol;
    pos4 = netheightRow;

else % merge "grids": assign the combined space to this axis
    left = leftFrameWidth+(i_col-1)*widthCol;
    bottom = 1-topFrameWidth-i_row*heightRow;
    right = left + netwidthCol;
    top = bottom + netheightRow;

    for i_sub = 1:length(currIdx)
        left(i_sub) = leftFrameWidth+(i_col(i_sub)-1)*widthCol;
        bottom(i_sub) = 1-topFrameWidth-i_row(i_sub)*heightRow;
        right(i_sub) = left(i_sub) + netwidthCol;
        top(i_sub) = bottom(i_sub) + netheightRow;
    end
    pos1 = min(left);
    pos2 = min(bottom);
    pos3 = max(right)-pos1;
    pos4 = max(top)-pos2;
end

%% make axis
ax = axes('Position',[pos1,pos2,pos3,pos4]);
% % subplot('Position',[pos1,pos2,pos3,pos4]);
% % small difference from using axes: subplot deletes any underlying subplots

% return identifier, if requested:
if(nargout > 0)
    ax_out = ax;
end

author: xiuyechen

repo: FishExplorer


text:
function [ f ] = getEEF_Moment( t )
%% This function is used to get the endeffector moment of the KUKA iiwa 7 R 800.

%% Syntax:
% [ M ] = getEEF_Moment( t )

%% About:
% This function is used to get the endeffector moment of the KUKA iiwa 7 R 800.
% The moment is measured at the media flange of the robot.

%% Arreguments:
% t: is the TCP/IP connection
% M: is 1x3 cell array. Represinting the  X,Y and Z components of the moment acting on the
% media flange of the robot. The components of the vector of moment are
% described in the robot base frame. The moment is measured in
% (Newton.meter)

% Copyright, Mohammad SAFEEA, 3rd of May 2017

theCommand='Eef_moment';
fprintf(t, theCommand);
message=fgets(t);
%disp(message)
[P1,N]=getDoubleFromString(message);

for i=1:3
    f{i}=P1{i};
end

end

function [jPos,j]=getDoubleFromString(message)
n=max(max(size(message)));
j=0;
numString=[];
for i=1:n
    if message(i)=='_'
        j=j+1;
        jPos{j}=str2num(numString);
        numString=[];
    else
        numString=[numString,message(i)];
    end
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function B = Bspline(t,k,u,v,ForceSup)
% :Usage:
% ::
%
%     Bspline(t,k,u[,v])
%
% Create B-Spline basis of order k, with knots u, evaluated at t.
% If control verticies v are specified then then B is the spline
% function instead of the basis.
%
% u must be at least length(k)+1
%
% ..
%    $Id: Bspline.m,v 1.3 1998/11/09 00:55:06 nicholst Exp $
% ..

global iB

if (nargin<3)
  help Bspline
  return
end

if (k+1>length(u))
  error('u must be at least length k+1')
end

% Silent flag for forcing the support to be defined between u(k)
% and u(end-k+1)
if (nargin<5)
  ForceSup = 1;
end

if ((nargin>3) & (length(v)+k ~= length(u)) & (ForceSup))
    error(sprintf('%d knots requires %d control verticies', ...
	length(u), length(u)-k))
end

% columnize
t=t(:);
u=u(:);

nBasis = length(u)-k;
B = zeros(length(t),nBasis);
iB = zeros(1,nBasis);

for i=1:nBasis
   B(:,i) = recu(t,i,k,u);
   iB(i) = (u(i+k)-u(i))/k;
end

% zero outside of valid range, if there's enough
if (length(u)>=2*k)
  if (ForceSup)
    B(t<u(k) | u(end-k+1)<t,:) = 0;
  end
else
  warning('Insufficient knots to be a proper spline basis')
end

if (nargin>3)
  B = B*v(:);
end

return

function B = recu(t,i,k,u)

if (k==1)

  if (u(i)==u(i+1))
    B = zeros(size(t));
  else
    B = (u(i)<=t) & (t<u(i+1));
  end

else

  B = w(t,i,k,u).*recu(t,i,k-1,u)+(1-w(t,i+1,k,u)).*recu(t,i+1,k-1,u);

end

function wt = w(t,i,k,u)

if (u(i)~=u(i+k-1))
  wt = (t-u(i))./(u(i+k-1)-u(i));
else
  wt = 0;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function code=sparse_code_wrapper(test)
    tic;
	n=length(test.X); dim=size(test.X{1},1);
    if iscell(test.B)
        code=zeros(length(test.B),n);
    else
        code = zeros(size(test.B,2),n);
    end
	if ~isfield(test,'gamma')
		test.gamma=0;
	end
	for i=1:n
		if isfield(test,'X')
			out = sdivLasso(test.B, test.X{i}, test.lambda, test.numiter, test.algo, test.gamma);
			code(:,i) = out.x;
		end
	end
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [params, names] = rbfperiodic2KernExtractParam(kern)

% RBFPERIODIC2KERNEXTRACTPARAM Extract parameters from the RBFPERIODIC2 kernel structure.
% FORMAT
% DESC extracts parameters from the RBF periodic covariance with variying period
% kernel structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the RBF periodic covariance with variying period
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO rbfperiodic2KernParamInit, rbfperiodic2KernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2009
%
%
% MODIFICATIONS : Andreas C. Damianou, 2011
%
% MODIFICATIONS : Michalis K. Titsias, 2011

% KERN

params = [kern.inverseWidth kern.variance kern.factor];
if nargout > 1
  names={'inverse width', 'variance', 'factor'};
end


text:
function CM = Tpm_CM(P,X)
% TPM_CM  Chenoweth-Martin two-phase multiplier
% TPM_CM(P,X) Returns the Chenoweth-Martin two-phase
% multipier for a steam-water system
%  Called function: h2o_rhof(P), h2o_rhog(P)
%  Required Inputs are: P - pressure (kPa), X - quality (fraction)
% ---------------------------------------------------------------
% The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

format long g;  % set the format of the calculations

if X<=0
   CM=1;
else
   RF=h2o_rhof(P);
   RG=h2o_rhog(P);
% Liquid Volume Fraction
   LVF=1-(X*RF/(X*RF+(1-X)*RG));
%The two-phase multiplier is
   CM=LVF^(-0.8642);
end

return


text:
function out1 = Sindy_ODE_RHS(in1,in2)
%SINDY_ODE_RHS
%    OUT1 = SINDY_ODE_RHS(IN1,IN2)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    31-Jul-2019 21:06:56

dz1 = in2(:,1);
dz2 = in2(:,2);
dz4 = in2(:,4);
z1 = in1(:,1);
z2 = in1(:,2);
t2 = z2.*-1.0;
t3 = t2+z1;
out1 = dz1.*-3.235001066692433e-2+dz2.*6.00492565527997e-3+sin(z1).*3.796608727034097e1-dz2.^2.*sin(t3).*3.229288227581633e-1-dz4.*cos(t3).*3.229288227581664e-1;

author: dynamicslab

repo: SINDy-PI


text:
function [prod,deriv] = mm_special(w,extractA,extractB)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
%   [vec(A);vec(B)] --> vec(A*B)
%
% where
%   A is extractA(w)
%   B is extractB(w)

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

if isempty(w)
    prod = @(w)mm_special(w,extractA,extractB);
    return;
end

if isa(w,'function_handle')
    outer = mm_special([],extractA,extractB);
    prod = compose_mv(outer,w,[]);
    return;
end

w = w(:);
A = extractA(w);
[m,k] = size(A);

B = extractB(w);
[k2,n] = size(B);
assert(k==k2,'inner matrix dimensions must agree');

M = A*B;
prod = M(:);

deriv = @(g2) deriv_this(g2);

function [g,hess,linear] = deriv_this(g2)
g = vJ_this(g2,A,B);
linear = false;
hess = @(w) hess_this(g2,w);
end

function [h,Jv] = hess_this(g2,w)
h = vJ_this(g2,extractA(w),extractB(w));
if nargout>=2
    Jv = Jv_this(w);
end
end

function prod = Jv_this(w)
Aw = extractA(w);
Bw = extractB(w);
M = Aw*B + A*Bw;
prod = M(:);
end

function w = vJ_this(prod,A,B)
M = reshape(prod,m,n);
Bp = A.'*M;
Ap = M*B.';
w = [Ap(:);Bp(:)];
end

end

function A = extractA_this(w,m,k)
A = w(1:m*k);
A = reshape(A,m,k);
end

function B = extractB_this(w,m,k,n)
B = w(m*k+(1:k*n));
B = reshape(B,k,n);
end

function test_this()

m = 4;
k = 5;
n = 6;

A = randn(m,k);
B = randn(k,n);

w = [A(:);B(:)];

extractA = @(w) extractA_this(w,m,k);
extractB = @(w) extractB_this(w,m,k,n);

f = mm_special([],extractA,extractB);
test_MV2DF(f,w);
end

author: nesl


text:
function view=segmentationMenu(view)
% view = segmentationMenu(view)
%
% Set up the callbacks for the segmentation/unfold menu
%
% djh, 2/14/2001
segmentationMenu = uimenu('Label','Segmentation','separator','on');

% Install new segmentation
%    installSegmentation;
callBackstr='installSegmentation;';
uimenu(segmentationMenu,'Label','(Re-)install segmentation','Separator','off',...
    'CallBack',callBackstr);

if strcmp(view.viewType,'Gray') || strcmp(view.viewType,'Volume')
    % Wrapper for mrFlatMesh
    %    flattenFromROI(view, ROI)
    callBackstr=['flattenFromROI(',view.name,');'];
    uimenu(segmentationMenu,'Label','Flatten (start from ROI)','Separator','on',...
        'CallBack',callBackstr);
end

% New flat callback:
%  newFlat;
callBackstr='newFlat;';
uimenu(segmentationMenu,'Label','Install New Unfold','Separator','off',...
    'CallBack',callBackstr);

% If flat view, pass arg to specify which flat to re-install
% Otherwise, user will be prompted to pick one of the flat subdirectories
if strcmp(view.viewType,'Flat')
    % Install new unfold
    %    installUnfold(view);
    callBackstr=['installUnfold(',view.name,'.subdir);'];
    uimenu(segmentationMenu,'Label','Reinstall unfold','Separator','off',...
        'CallBack',callBackstr);
else
    callBackstr='installUnfold;';
    uimenu(segmentationMenu,'Label','Reinstall unfold','Separator','off',...
        'CallBack',callBackstr);
end

if strcmp(view.viewType,'Volume') || strcmp(view.viewType,'Gray') || strcmp(view.viewType,'Flat')
    % Message box with segmentation information
    %    segInfo = segmentInfo(view)
    callBackstr=['segmentInfo(',view.name,');'];
    uimenu(segmentationMenu,'Label','Segmentation info','Separator','on',...
        'CallBack',callBackstr);
end

if strcmp(view.viewType,'Flat')
    % Xform Verify Gray-Flat correspondence callback:
    %    verifyGrayFlat(view);
    callBackstr=['checkCoordsNodes(',view.name,');'];
    uimenu(segmentationMenu,'Label','Verify Gray-Flat match','Separator','off',...
        'CallBack',callBackstr);
end

return

author: vistalab

repo: vistasoft


text:
function cFiles = defaultCFileList;
% Returns a list of all .c source files found in the VISTASOFT repository.
%
% cFiles = defaultCFileList;
%
% This list was found using grep (around October 2008). cFiles is a cell
% array of strings. Each entry is relative to the default VISTASOFT root
% path.
%
% Some mex files need to concatenate multiple source files. You can
% concatenate the source file names separated by a ';'. E.g.:
% '/some/dir/first.c;some/dir/second.c' will compile first.c as the mex
% file 'first.mexext' and link in the code from second.c.
%
% Some mex files need options for compilations (e.g. linking against system
% provided libraries). To pass in an option use a string that starts with
% '-' (e.g. '-lz').  This will then be skipped in the filepath expansion
% and just based to the mex command line.
%
% ras 10/2008.
% 2009.07.01 RFD: cleaned up a few details. Can now compile things that
% require multiple source files (e.g., niftiRead)

% This will ultimately be the output:
cFiles = {};

% We search our path for these files. Files grouped together, need to be
% compiled in that order: e.g. 'mex matToQuat.c nifti1_io.c znzlib.c' (with
% their respective paths, of course).

cFileNames = {...
    {'-lz','matToQuat.c', 'nifti1_io.c', 'znzlib.c'}...
    {'-lz','readFileNifti.c', 'nifti1_io.c', 'znzlib.c'}...
    {'-lz','writeFileNifti.c', 'nifti1_io.c', 'znzlib.c'}...
    {'mrManDist.c'}...
    {'assignToNearest.c'}...
    {'spm_hist2_weighted.c'}...
    {'spm_hist2_weighted_MI.c'}...
    {'coranal.c'} ...
    {'sumOfNeighbors.c'} ...
    {'updateTinC.c'} ...
    {'varOfNeighbors.c'} ...
    {'dtiFitTensor.c'} ...
    {'dtiJointHist.c'} ...
    {'dtiSplitTensor.c'} ...
    {'dtiTensorInterp_Pajevic.c'} ...
    {'ndfun.c'} ...
    {'magicwand1.c'} ...
    {'initInvDiagsInC.c'} ...
    {'trilin.c'}...
    {'dtiLoadConTrackPaths.c'}...
    {'mrAnatFastInterp3.c'}...
    {'myCinterp3.c'}...
    {'findNearest.c'}...
    {'pnet.c'}...
    {'mrManDist.c'}...
    {'assignToNearest.c'}...
    {'edgesCross2d.c'}...
    {'mCheck.c'}...
    {'mrCheckNode2dIntersect.c'}...
    {'evaluate.c'}...
    {'regHistogram.c'}...
    {'SagittalRot.c'}...
    {'ObliqueRot.c'}...
    {'SagittalRot.c'}...
    {'ObliqueRot.c'}...
    {'SagittalRot.c'}...
    {'regInplanes.c'}...
    {'corrDn.c'}...
    {'histo.c'}...
    {'corrDn.c'}...
    {'histo.c'}...
    {'wrap.c'}...
    {'convolve.c'}...
    {'upConv.c'}...
    {'range2.c'}...
    {'edges-orig.c'}...
    {'edges.c'}...
    {'pointOp.c'}...
    {'wrap.c'}...
    {'convolve.c'}...
    {'upConv.c'}...
    {'range2.c'}...
    {'edges-orig.c'}...
    {'edges.c'}...
    {'pointOp.c'}...
    {'corrDn.c'}...
    {'histo.c'}...
    {'wrap.c'}...
    {'convolve.c'}...
    {'upConv.c'}...
    {'range2.c'}...
    {'edges-orig.c'}...
    {'edges.c'}...
    {'pointOp.c'}...
    {'varOfNeighbors.c'}};

for ii=1:length(cFileNames)
    this_files = [];
    for jj=1:length(cFileNames{ii})

        %JMA - This is a bit kludgy and not robust but quick.
        %If the first character of the string is a "-" then treat it as a
        %option string and don't try to expand the full path of a filename
        if strncmp(cFileNames{ii}{jj},'-',1)
            stringToAdd =     cFileNames{ii}{jj};
        else
            stringToAdd = which(cFileNames{ii}{jj});
        end

        this_files = [this_files stringToAdd ' '];
    end
    cFiles{ii} = this_files;
end

return


text:
function vl_bench_bnorm(gpu)
  if nargin < 1
    gpu = false ;
  end

  T = 100 ;
  x = randn(64,64,32,32,'single') ;
  g = randn(32,1,'single') ;
  b = randn(32,1,'single') ;

  if gpu
    x = gpuArray(x) ;
    g = gpuArray(g) ;
    b = gpuArray(b) ;
  end

  tic
  for t=1:T
    y = vl_nnbnorm(x,g,b) ;
  end
  if gpu, wait(gpuDevice) ; end
  fprintf('new: %f\n',toc);

  tic
  for t=1:T
    y_ = vl_nnbnorm_old(x,g,b) ;
  end
  if gpu, wait(gpuDevice) ; end
  fprintf('old: %f\n',toc);

  dzdy = randn(size(y),'single') ;
  if gpu
    dzdy = gpuArray(dzdy) ;
  end

  tic
  for t=1:T
    [a,b,c] = vl_nnbnorm(x,g,b,dzdy) ;
  end
  if gpu, wait(gpuDevice) ; end
  fprintf('new deriv: %f\n',toc);

  tic
  for t=1:T
    [a_,b_,c_] = vl_nnbnorm_old(x,g,b,dzdy) ;
  end
  if gpu, wait(gpuDevice) ; end
  fprintf('old deriv: %f\n',toc);

  vl_testsim(y,y_);
  vl_testsim(a,a_);
  vl_testsim(b,b_);
  vl_testsim(c,c_);
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function model = foo(inputs, outputs, varargin)

model.inputs = inputs;
model.outputs = outputs;


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

m = size(X,1);

for i = 1:m
    for j = 1:K
        idx(i,j) = (X(i,:)-centroids(j,:))*(X(i,:)-centroids(j,:))';
    end
end
[dummy idx] = min(idx,[],2);

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

end

repo: MachineLearningEx


text:
function gainS = computeGainS(n,w)
    gainS = w*eye(n);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function cMap=spectral(varargin)

%%

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

%%

cMap=[0.368627450980392,0.309803921568627,0.635294117647059;...
      0.196078431372549,0.533333333333333,0.741176470588235;...
      0.400000000000000,0.760784313725490,0.647058823529412;...
      0.670588235294118,0.866666666666667,0.643137254901961;...
      0.901960784313726,0.960784313725490,0.596078431372549;...
      1,1,0.749019607843137;...
      0.996078431372549,0.878431372549020,0.545098039215686;...
      0.992156862745098,0.682352941176471,0.380392156862745;...
      0.956862745098039,0.427450980392157,0.262745098039216;...
      0.835294117647059,0.243137254901961,0.309803921568627;...
      0.619607843137255,0.00392156862745098,0.258823529411765];

if n~=size(cMap,1)
    cMap=resampleColormap(cMap,n);
end

%%
% _*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/>.

repo: GIBBON