text:
%computes the novelty measure used by Hainsworth
%> called by ::ComputeNoveltyFunction
%>
%> @param X: spectrogram (dimension FFTLength X Observations)
%> @param f_s: sample rate of audio data (unused)
%>
%> @retval d_hai novelty measure
% ======================================================================
function [d_hai] = NoveltyHainsworth (X, f_s)

    epsilon     = 1e-5;

    % difference spectrum
    X           = [X(:,1), X];
    X(X<=0)     = epsilon;

    % flux
    d_hai       = sum(log2(X(:,2:end)./X(:,1:end-1)))/size(X,1);
end

author: alexanderlerch

repo: ACA-Code


text:
%KNNM Trainable K-Nearest Neighbour density estimation
%
%   W = KNNM(A,KNN)
%   W = A*KNNM([],KNN)
%   W = A*KNNM(KNN)
%
%   D = B*W
%
% INPUT
% 	A     Dataset used for training
%   B     Dataset used for evaluation
%   KNN   Number of nearest neighbours
%
% OUTPUT
%   W     Density estimate
%
% DESCRIPTION
% A density estimator is constructed based on the k-Nearest Neighbour rule
% using the objects in A. In case A is labeled, density estimates are
% performed classwise and combined by the class priors. The default KNN is
% the square root of the size of the class. The data is scaled by variance
% normalisation determined by the training set.
%
% The mapping W may be applied to a new dataset B using DENSITY = B*W.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, KNNC, PARZENM, GAUSSM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function w = knnm(varargin)

  mapname = 'KNN density estimation';
  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],[]);

  if mapping_task(argin,'definition')

    w = define_mapping(argin,'untrained',mapname);

	elseif mapping_task(argin,'training')			% Train a mapping.

		[a,knn] = deal(argin{:});
		if isa(a,'prdataset')
			labname = getname(a);
		else
			labname = '';
		end
		if (~isdataset(a) && ~isdatafile(a)) || getsize(a,3) ~= 1
			w = mclassm(a,prmapping(mfilename,knn),'weight');
			w = setlabels(w,labname);
			w = setname(w,mapname);
			return
		end
		islabtype(a,'crisp');
		isvaldfile(a,1);
		a = testdatasize(a,'objects');
		a = remclass(a);
		[m,k] = size(a);
    knn = setdefaults({knn},round(sqrt(m)));
		w = prmapping(mfilename,'trained',{a,knn},labname,k,1);
		w = setname(w,mapname);

	else															% Execute a trained mapping V.

    [a,v] = deal(argin{:});
    [b,knn] = getdata(v);
		[m,k] = size(b);
		u = scalem(b,'variance');
		a = a*u;
		b = b*u;
		d = sqrt(distm(+a,+b));					% Calculate squared distances
		s = sort(d,2);						    	% and find nearest neighbours.
		ss = s(:,knn);									% Find furthest neighbour
		ss(ss==0) = realmin.^(1/k);	    % Avoid zero distances
		f = knn./(m*nsphere(k)*(ss.^k));% Normalize by the volume of sphere
																		% defined by the furthest neighbour.
		f = f*prod(u.data.rot);         % Correct for scaling

		w = setdat(a,f,v);							% proper normalisation has still to be done.

	end

return

function v = nsphere(k)

	% Volume k-dimensional sphere
	v = (2*(pi^(k/2))) / (k*gamma(k/2));

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function lp = create_lp_solve_model(A,b,f,xint,LB,UB,e,options);

[m,n] = size(A);
lp = lp_solve('make_lp', m, n);

lp_solve('set_mat', lp, A);
lp_solve('set_rh_vec', lp, b);
lp_solve('set_obj_fn', lp, f);
lp_solve('set_maxim', lp); % default is solving minimum lp.
for i = 1:length(e)
    if e(i) < 0
        con_type = 1;
    elseif e(i) == 0
        con_type = 3;
    else
        con_type = 2;
    end
    lp_solve('set_constr_type', lp, i, con_type);
end
for i = 1:length(LB)
    %if ~isinf(LB(i))
        lp_solve('set_lowbo', lp, i, LB(i));
    %end
end
for i = 1:length(UB)
   if ~isinf(UB(i))
        lp_solve('set_upbo', lp, i, UB(i));
   end
end
for i = 1:length(xint)
    lp_solve('set_int', lp, xint(i), 1);
end

if options.lpsolve.scalemode~=0
    lp_solve('set_scaling', lp, scalemode);
end

% for i = 1:length(sos)
%     lp_solve('add_SOS', lp, ['dummy' num2str(i)], 1, i, sos{i}, 1:length(sos{i}));
% end

switch options.verbose
    case 0
        lp_solve('set_verbose', lp, 0);%options.verbose)
    case 1
        lp_solve('set_verbose', lp, 4);%options.verbose)
    case 2
        lp_solve('set_verbose', lp, 5);%options.verbose)
    otherwise
        lp_solve('set_verbose', lp, 6);%options.verbose)
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function score = Min_value(Population,~)
% <min> <single> <real/integer/label/binary/permutation> <large/none> <constrained/none> <expensive/none> <sparse/none> <dynamic/none>
% The minimum objective value

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

	score = Population.best.objs;
    if isempty(score)
        score = nan;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [h0, h1, g0, g1] = daubf(K,str);
% [h0, h1, g0, g1] = daubf(K);
% h0 - low-pass analysis
% h1 - high-pass analysis
% g0 - low-pass analysis
% g1 - high-pass analysis
%
% K zeros at z=-1
% Use daubf(K,'mid') for mid-phase type

% Ivan Selesnick
% selesi@nyu.edu
% NYU - School of Engineering

[h,s,g] = maxflatI(K,K-1);

r = roots(g);
r = r(abs(r) < 1);
q = real(poly(r));

if nargin > 1
	if strcmp(str,'mid')
		q = sfact_mid(g);
	end
end

q = q/sum(q);                 % normalize
h0 = q;                       % set  h0 = q;
for k = 1:K                   % make h0 = q * [(z^(-1)+1)/2]^K
   h0 = conv(h0,[1 1]/2);
end
h0 = sqrt(2)*h0;              % normalize so that sum(h0) = sqrt(2)

h1 = h0(end:-1:1);
h1(2:2:end) = -h1(2:2:end);

g0 = h0(end:-1:1);
g1 = h1(end:-1:1);

author: thomaskuestner


text:
function success = write_bspline_coeff_file(bspFileName,bsplineParamsS)
% function success = write_bspline_coeff_file(bspFileName,bsplineParamsS)
%
% APA, 07/17/2012

bsp_img_origin          = bsplineParamsS.bsp_img_origin;
bsp_img_spacing         = bsplineParamsS.bsp_img_spacing;
bsp_img_dim             = bsplineParamsS.bsp_img_dim;
bsp_roi_offset          = bsplineParamsS.bsp_roi_offset;
bsp_roi_dim             = bsplineParamsS.bsp_roi_dim;
bsp_vox_per_rgn         = bsplineParamsS.bsp_vox_per_rgn;
bsp_direction_cosines   = bsplineParamsS.bsp_direction_cosines;
bsp_coefficients        = bsplineParamsS.bsp_coefficients;

fid = fopen(bspFileName,'wb');

fprintf(fid,'MGH_GPUIT_BSP <experimental>\n');
fprintf(fid,['img_origin = ', sprintf('%0.20g ',bsp_img_origin),'\n']);
fprintf(fid,['img_spacing = ', sprintf('%0.20g ',bsp_img_spacing),'\n']);
fprintf(fid,['img_dim = ', sprintf('%0.20g ',bsp_img_dim),'\n']);
fprintf(fid,['roi_offset = ', sprintf('%0.20g ',bsp_roi_offset),'\n']);
fprintf(fid,['roi_dim = ', sprintf('%0.20g ',bsp_roi_dim),'\n']);
fprintf(fid,['vox_per_rgn = ', sprintf('%0.20g ',bsp_vox_per_rgn),'\n']);
fprintf(fid,['direction_cosines = ', sprintf('%0.20g ',bsp_direction_cosines),'\n']);
fprintf(fid,sprintf('%0.20g\n',bsp_coefficients));

%cell2file(fileC,bspFileName);
fclose(fid);

success = 1;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function b = subsref(a,S)
% B = A.FIELDNAME; % SUBSREF - subscripted indexing
%
% See also GET

switch length(S)
case 1
    switch S.type
    case '()'
        b = get(a,S.subs{:});
    case '.'
        b = get(a,S.subs);
    end;
case 2
    if strcmp(S(1).type,'()')&strcmpi(S(2).type,'.')    % a(i).field
        b = get(a,S(1).subs{:},S(2).subs);
    elseif strcmp(S(1).type,'.')&strcmpi(S(2).type,'()')% a.field(i)
        b = get(a,S(1).subs,S(2).subs{:});
    elseif strcmp(S(1).type,'.')&strcmpi(S(2).type,'.') % a.field1.field2
        b = get(a,[S(1).subs,'.',S(2).subs]);
    else
        error('Sorry - output structure not supported');
    end;
otherwise
    error('Sorry - output structure not supported');
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function tf = ischarlike(v)
    % returns true if argument is either a char array or a string scalar
    tf = ischar(v) || (isstring(v) && isscalar(v));

author: CelsoReyes


text:
function line_fekete_rule_test ( )

%*****************************************************************************80
%
%% LINE_FEKETE_RULE_TEST tests the LINE_FEKETE_RULE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_FEKETE_RULE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LINE_FEKETE_RULE library.\n' );

  for m = [ 5, 11, 21 ]
    line_fekete_rule_test01 ( m );
  end

  for m = [ 5, 11, 21 ]
    line_fekete_rule_test02 ( m );
  end

  m = 21;
  line_fekete_rule_bos_levenberg_test ( m );

  for m = [ 5, 11, 21 ]
    line_fekete_rule_test03 ( m );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_FEKETE_RULE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Z = minus(A,B)
%minus           Overloaded

% Standard case A.cx + A.logdetP - (B.cx + B.logdetP)

if isa(A,'sdpvar') | isa(A,'double')
    if prod(size(A))>1
        error('Only scalar terms can be added to a logdet term');
    end
    Acx = A;
    Alog = {};
    Again = [];
else
    Z = A;
    Acx  = A.cx;
    Alog = A.P;
    Again = A.gain;
end

if isa(B,'sdpvar') | isa(B,'double')
    if prod(size(B))>1
        error('Only scalar terms can be added to a logdet term');
    end
    Bcx = B;
    Blog = {};
    Bgain = [];
else
    Z = B;
    Bcx  = B.cx;
    Blog = B.P;
    Bgain = B.gain;
end

if isempty(Acx)
    Acx = 0;
end
if isempty(Bcx)
    Bcx = 0;
end

if isequal(Acx-Bcx,0)
    Z.cx = [];
else
    Z.cx = Acx - Bcx;
end

Z.P = {Alog{:},Blog{:}};
Z.gain = [Again -Bgain];


text:
function [association] = ft_nonlinearassociation(cfg, data)

% NONLINEARASSOCIATION calculate the association coefficient as a
% function of delay.
%
% In order to estimate the amount of association between all possible
% pairs of MEG sensors the nonlinear association analysis is used.
% It was first developed for EEG data analysis by Pijn and co-workers
% (Lopes da Silva, et al. 1989; Pijn, et al. 1990). The basic principle
% is similar to that of coherence and (cross) correlation, with the
% exception that this nonlinear association method can be applied
% independent of the type of relationship (linear or nonlinear) in
% the data.
%
% The method is based on the idea that if two signals x and y are
% correlated, a nonlinear regression curve can be calculated that
% represents their relationship. In practice, that regression curve
% is estimated by creating a scatterplot of y versus x, dividing the
% data in segments and describing each segment with a linear regression
% curve. The estimated correlation ratio h2, which gives the reduction
% in variance of y as a result of predicting its values according to
% the regression curve, can be calculated as follows:
%
% h^2 = (sum(Yi - mean(Y))^2 - sum(Yi - f(Xi))^2) / sum(Yi - mean(Y))^2
%
% With the sum going over N samples and f(Xi) the estimated value of
% Yi according to the regression line. The h2 coefficient has values
% between 0 (y is completely independent of x) and 1 (y is completely
% determined by x). In the case of a linear relationship between x
% and y, h2 is equal to the well known Pearson correlation coefficient
% (r2). As is the case with cross-correlation, it is possible to
% estimate h2 as a function of time shift () between the signals. The
% h2 is then iteratively calculated for different values of , by
% shifting the signals in comparison to each other, and the value for
% which the maximal h2 is reached can be used as an estimate of the
% time lag between both signals. In deciding what epoch length to use
% in the association analysis, a trade-off has to be made between
% successfully determining the correct delay and h2-value (for which
% large epoch lengths are necessary) and a small enough time-resolution
% (for which small epoch lengths are necessary).
%
% Use as
%   [association] = ft_nonlinearassociation(cfg, data)
%
% The input data should be organised in a structure as obtained from
% the PREPROCESSING function.
%
% The configuration should contain
%   cfg.channel    = Nx1 cell-array with selection of channels (default = 'all'), see CHANNELSELECTION for details
%   cfg.keeptrials = 'yes' or 'no', process the individual trials or the concatenated data (default = 'no')
%   cfg.trials     = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.fsample    = 1200
%   cfg.maxdelay   = 32/cfg.fsample
%   cfg.delaystep  = 2/cfg.fsample
%   cfg.nr_bins    = 7
%   cfg.offset     = 0
%   cfg.order      = 'Hxy'
%   cfg.timwin     = 0.2
%   cfg.toi        = []
%
% References
% - Lopes da Silva F, Pijn JP, Boeijinga P. (1989): Interdependence of
% EEG signals: linear vs. nonlinear associations and the significance
% of time delays and phase shifts. Brain Topogr 2(1-2):9-18.
% - Pijn JP, Vijn PC, Lopes da Silva FH, Van Ende Boas W, Blanes W.
% (1990): Localization of epileptogenic foci using a new signal
% analytical approach. Neurophysiol Clin 20(1):1-11.

% Copyright (C) 2007, Inge Westmijse

ft_defaults

% set the defaults
if ~isfield(cfg, 'trials'),        cfg.trials  = 'all';             end
if ~isfield(cfg, 'keeptrials'),    cfg.keeptrials  = 'no';          end
if ~isfield(cfg, 'channel'),       cfg.channel = 'all';             end
if ~isfield(cfg, 'toi'),           cfg.toi = [];                    end
if ~isfield(cfg, 'timwin'),        cfg.timwin = 0.2;                end
if ~isfield(cfg, 'fsample'),       cfg.fsample = 1200;              end
if ~isfield(cfg, 'delaystep'),     cfg.delaystep = 2/cfg.fsample;   end
if ~isfield(cfg, 'maxdelay'),      cfg.maxdelay = 32/cfg.fsample;   end
if ~isfield(cfg, 'offset'),        cfg.offset = 0;                  end
if ~isfield(cfg, 'nr_bins'),       cfg.nr_bins = 7;                 end
if ~isfield(cfg, 'order'),         cfg.order = 'Hxy';               end
if ~isfield(cfg, 'feedback'),      cfg.feedback = 'textbar';        end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do some bookkeeping on the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% select trials of interest
if ~strcmp(cfg.trials, 'all')
  if islogical(cfg.trials),  cfg.trials=find(cfg.trials);  end
  fprintf('selecting %d trials\n', length(cfg.trials));
  data.trial  = data.trial(cfg.trials);
  data.time   = data.time(cfg.trials);
end
Ntrials  = numel(data.trial);

% select channels of interest
cfg.channel = channelselection(cfg.channel, data.label);
chansel = match_str(data.label, cfg.channel);
fprintf('selecting %d channels\n', length(chansel));
for trial=1:Ntrials
  data.trial{trial} = data.trial{trial}(chansel,:);
end
data.label = data.label(chansel);
Nchans     = length(chansel);

% determine the size of each trial, they can be variable length
Nsamples = zeros(1,Ntrials);
for trial=1:Ntrials
  Nsamples(trial) = size(data.trial{trial},2);
end

if strcmp(cfg.keeptrials, 'no')
  % concatenate all the data into a 2D matrix
  fprintf('concatenating data');
  dat = zeros(Nchans, sum(Nsamples));
  for trial=1:Ntrials
    fprintf('.');
    begsample = sum(Nsamples(1:(trial-1))) + 1;
    endsample = sum(Nsamples(1:trial));
    dat(:,begsample:endsample) = data.trial{trial};
  end
  fprintf('\n');
  fprintf('concatenated data matrix size %dx%d\n', size(dat,1), size(dat,2));
  time = [1/cfg.fsample : size(dat,1)/cfg.fsample : size(dat,1)];
  data.trial = {dat};
  data.time  = {time};
  Ntrials    = 1;
else
  % replace the time axis, since shifted time-axes over trials are not supported
  for i = 1:Ntrials
    data.time{i} = [1/cfg.fsample : 1/cfg.fsample : size(data.trial{i},2) / cfg.fsample];
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare all data selection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% item 1: channel combinations
numchannelcmb = Nchans*(Nchans-1)/2;
channelcmb = zeros(numchannelcmb, 2);
k = 1;

if (strcmp(cfg.order, 'Hxy'))
  for i=1:Nchans
    for j=(i+1):Nchans
      channelcmb(k, :) = [i j];
      k = k+1;
    end
  end
elseif (strcmp(cfg.order, 'Hyx'))
  for i = Nchans:-1:1
    for j = Nchans : -1:(i+1)
      channelcmb(k,:) = [i j];
      k = k+1;
    end
  end
end
numchannelcmb = size(channelcmb,1);

% item 2: time windows
% this requires
%    cfg.toi = [t1 t2 t3 t4]    center time of each window
%    cfg.timwin = 0.2           in seconds
% TODO implement "time, offset, fsample"

for j = 1:Ntrials
  time = data.time{j};
  numtimwin = size(cfg.toi,2); % initial guess, not all might fit in this trial
  timwin = zeros(numtimwin,2);
  sel = zeros(numtimwin,1);
  for i=1:numtimwin
    timwin(i,1) = cfg.toi(i) - cfg.timwin/2;  % begin of time window
    timwin(i,2) = cfg.toi(i) + cfg.timwin/2;  % end of time window
    sel(i) = timwin(i,1)>=time(1) && timwin(i,2)<=time(end);  % does it fit in?
  end
  timwin        = timwin(find(sel==1),:);
  toi_mat{j}    = cfg.toi(find(sel == 1));                        % update the configuration
  timwin_mat{j} = round((timwin - cfg.offset) * cfg.fsample);     % convert to samples
end

timwin = timwin_mat;
cfg.toi = toi_mat;

% item 3: delays within each timewindow
% this requires
%   cfg.delaystep
%   cfg.maxdelay
delay = -cfg.maxdelay:cfg.delaystep:cfg.maxdelay;
numdelay = length(delay);
% convert to samples
delay = round(delay * cfg.fsample);

if strcmp(cfg.keeptrials, 'yes')
  for trllop=1:Ntrials
    association.trial(trllop) = Do_Association_Calculation( trllop , Ntrials , cfg , timwin , numchannelcmb , numdelay , data , channelcmb , delay );
  end % for each trial
else
  trllop = 1;
  association = Do_Association_Calculation( trllop , Ntrials , cfg , timwin , numchannelcmb , numdelay , data , channelcmb , delay );
end

% add the version details of this function call to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id   = '$Id$';
% remember the configuration details of the input data
try, cfg.previous = data.cfg; end
% remember the exact configuration details in the output
association.cfg = cfg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that does the  computation for each trial
function [association] = Do_Association_Calculation( trllop , Ntrials , cfg , timwin , numchannelcmb , numdelay , data , channelcmb , delay )

fprintf('\nprocessing trial %d from %d\n', trllop, Ntrials);

association = [];
numtimwin = size(timwin{trllop},1);
h2 = zeros(numchannelcmb, numtimwin, numdelay);

l = 0;
maxl = numchannelcmb*numtimwin*numdelay;
progress('init', cfg.feedback, 'Computing nonlinear association for each delay');

for i=1:numchannelcmb
  dat1_chan = data.trial{trllop}(channelcmb(i,1),:);
  dat2_chan = data.trial{trllop}(channelcmb(i,2),:);
  progress(l/maxl,'Computing nonlinear association %d from %d', l, maxl);

  for j=1:numtimwin
    dat1_timwin = dat1_chan(timwin{trllop}(j,1):timwin{trllop}(j,2));
    dat2_timwin = dat2_chan(timwin{trllop}(j,1):timwin{trllop}(j,2));

    for k=1:numdelay
      if delay(k)==0
        % select the complete (unshifted) window
        dat1_delay = dat1_timwin;
        dat2_delay = dat2_timwin;
      elseif delay(k)<0
        % channel 1 is shifted w.r.t. channel 2
        dat1_delay = dat1_timwin((abs(delay(k))+1):end);
        dat2_delay = dat2_timwin(1:(end-abs(delay(k))));
      elseif delay(k)>0
        % channel 2 is shifted w.r.t. channel 1
        dat1_delay = dat1_timwin(1:(end-delay(k)));
        dat2_delay = dat2_timwin((delay(k)+1):end);
      end

      % remove the mean of each snippet of data
      dat1_delay = dat1_delay - mean(dat1_delay);
      dat2_delay = dat2_delay - mean(dat2_delay);

      % do the computation

      [sorted_data1,id] = sort(dat1_delay);
      sorted_data2     = dat2_delay(id);

      data_is = sorted_data1;
      data_js = sorted_data2;

      % divided data_i in bins
      [H,mp_i] = hist(data_is,cfg.nr_bins);

      % Divide data_i and data_j in bins
      bp = 1;
      gem_j = zeros (cfg.nr_bins,1);
      for s = 1:cfg.nr_bins
        gem_j(s) = mean(data_js(bp:bp+H(s)-1));
        bp = bp + H(s);
      end

      % Calculation of line segment and variance per bin
      p=1;
      sp = 1;

      clear unex_var tot_var;

      for u = 1:cfg.nr_bins
        data_is_bin = data_is(sp:sp+H(u)-1)';
        data_js_bin = data_js(sp:sp+H(u)-1)';

        if u == 1
          fx1 = gem_j(u) + (gem_j(u+1) - gem_j(u))/(mp_i(u+1) - mp_i(u))*(data_is_bin(1:round(length(data_is_bin)/2)) - mp_i(u));
        else
          fx1 = gem_j(u-1) + (gem_j(u) - gem_j(u-1))/(mp_i(u) - mp_i(u-1))*(data_is_bin(1:round(length(data_is_bin)/2)) - mp_i(u-1));
        end
        if u == cfg.nr_bins
          fx2 = gem_j(u-1) + (gem_j(u) - gem_j(u-1))/(mp_i(u) - mp_i(u-1))*(data_is_bin(round(length(data_is_bin)/2)+1:end) - mp_i(u-1));
        else
          fx2 = gem_j(u) + (gem_j(u+1) - gem_j(u))/(mp_i(u+1) - mp_i(u))*(data_is_bin(round(length(data_is_bin)/2)+1:end) - mp_i(u));
        end

        % calculation unexplained variance
        ftot = [fx1; fx2];
        unex_var(p:p+length(data_js_bin)-1,:) = (data_js_bin - ftot).^2;

        % calculation total variance
        tot_var(p:p+length(data_js_bin)-1,:) = (data_js_bin - mean(data_js)).^2;

        p = p+length(data_is_bin);
        sp = sp + H(u);
      end

      % Calculation of association coefficient h2
      h2(i,j,k) = ((sum(tot_var) - sum(unex_var)) / sum(tot_var))*100;
      l=l+1;
    end
  end
end

progress('close');

% collect the results
if (size(h2, 1) ~= numchannelcmb)
  error('number of channel combinations is not right');
end
if (size(h2, 2) ~= numtimwin)
  error('number of timewindows does not match input');
end
if (size(h2, 3) ~= numdelay)
  error('number of delays does not match input');
end
h2 = reshape(h2, numchannelcmb*numtimwin, numdelay);
[m, indx] = max(h2, [], 2);
h2 = reshape(m, numchannelcmb, numtimwin);
d = (delay(indx)./cfg.fsample)*1000; % convert to miliseconds (to compare with original method)
d = reshape(d, numchannelcmb, numtimwin);

% FIXME this does not work: one is not supposed to rely on data.cfg.trl, use data.time please!
% association.time = cfg.toi{trllop} + data.cfg.trl(trllop,1);

association.h2 = h2;
association.delay = d;

return % from Do_Association_Calculation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that computes the mean over all values without checking the dimensions
% this function is approximately 2x faster on 1 dimensional data than the matlab version
function y = mean(x)
y = sum(x(:))./numel(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that bins the elements of Y into equally spaced containers
% this function is approximately 2x faster on 1 dimensional data than the matlab version
function [nn, x] = hist(y, x)
miny = min(y);
maxy = max(y);
binwidth = (maxy - miny) ./ x;
xx = miny:binwidth:maxy;
x  = xx(1:(end-1)) + binwidth/2;
% Shift bins so the interval is ( ] instead of [ ).
xx = full(real(xx)); y = full(real(y)); % For compatibility
bins = xx + eps(xx);
nn = histc(y',[-inf bins],1);
% Combine first bin with 2nd bin and last bin with next to last bin
nn(2,:) = nn(2,:)+nn(1,:);
nn(end-1,:) = nn(end-1,:)+nn(end,:);
nn = nn(2:end-1,:);


text:
%IMPEXP_GMSH Import/export grid in GMSH format.
%
%   [ GRID ] = IMPEXP_GMSH( FILE_NAME, MODE, DATA, IIND, FID_LOG, BIX )
%   Import and export of Gmsh (.msh) Ascii mesh and grid file format.
%
%   FILE_NAME is a string specifying the (root) file name to process.
%
%   MODE can either be a string indicating IMPORT (where boundary
%   numbering are reconstructed from Gmsh physical or elemental tags if
%   available), IMPORT_BDR (boundary reconstruction with gridbdr), or
%   EXPORT.
%
%   With the IMPORT mode option the function returns a GRID struct. While
%   for exporting, the input DATA can be given as either a valid fea
%   finite element data struct, or a FEATool grid struct.
%
%   IIND is an optional array for physical domains (subdomains in 2D and
%   boundaries in 3D) to extract and import (default all).
%
%   FID_LOG is an optional log file handle for message output
%   (negative for Gui output or empty for no output).
%
%   BIX is an optional index vector to correct boundary numbering.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function graph = prtGraphDataGenKarate
% graph = prtGraphDataGenKarate
%   Read karate.gml
%
% From the README:
%
% The file karate.gml contains the network of friendships between the 34
% members of a karate club at a US university, as described by Wayne
% Zachary in 1977.  If you use these data in your work, please cite W. W.
% Zachary, An information flow model for conflict and fission in small
% groups, Journal of Anthropological Research 33, 452-473 (1977).

baseDir = prtGraphDataDir;
gmlFile = 'karate.gml';
file = fullfile(baseDir,gmlFile);

[sparseGraph,names] = prtUtilReadGml(file);
graph = prtDataTypeGraph(sparseGraph,names);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function convertNetworkToFastRcnn(obj, varargin)
% convertNetworkToFastRcnn(obj, varargin)
%
% Modify network for Fast R-CNN's ROI pooling.
%
% Copyright by Holger Caesar, 2015
% Updated by Jasper Uijlings:
%  - Extra flexibility and possible bounding box regression
%  - Added instanceWeights to loss layer

% Initial settings
p = inputParser;
addParameter(p, 'lastConvPoolName', 'pool5');
addParameter(p, 'finalFCName', 'fc8');
parse(p, varargin{:});

lastConvPoolName = p.Results.lastConvPoolName;
finalFCName = p.Results.finalFCName;

fprintf('Converting object classification network to Fast R-CNN network (ROI pooling, box regression, etc.)...\n');

%%% Add weights to loss layer. Note that this field remains empty when not
% given as input. So the loss layers should ignore empty weights.
softmaxInputs = obj.net.layers(obj.net.getLayerIndex('softmaxloss')).inputs;
if ~ismember('instanceWeights', softmaxInputs)
    softmaxInputs{end+1} = 'instanceWeights';
    obj.net.setLayerInputs('softmaxloss', softmaxInputs);
end

%%% Replace pooling layer of last convolution layer with roiPooling
lastConvPoolIdx = obj.net.getLayerIndex(lastConvPoolName);
assert(~isnan(lastConvPoolIdx));
roiPoolName = ['roi', lastConvPoolName];
firstFCIdx = obj.net.layers(lastConvPoolIdx).outputIndexes;
assert(length(firstFCIdx) == 1);
roiPoolSize = obj.net.layers(firstFCIdx).block.size(1:2);
roiPoolBlock = dagnn.RoiPooling('poolSize', roiPoolSize);
replaceLayer(obj.net, lastConvPoolName, roiPoolName, roiPoolBlock, {'oriImSize', 'boxes'}, {'roiPoolMask'});

%%% Add bounding box regression layer
if obj.nnOpts.bboxRegress
    finalFCLayerIdx = obj.net.getLayerIndex(finalFCName);
    inputVars = obj.net.layers(finalFCLayerIdx).inputs;
    finalFCLayerSize = size(obj.net.params(obj.net.layers(finalFCLayerIdx).paramIndexes(1)).value);
    regressLayerSize = finalFCLayerSize .* [1 1 1 4]; % Four times bigger than classification layer
    regressName = [finalFCName 'regress'];
    obj.net.addLayer(regressName, dagnn.Conv('size', regressLayerSize), inputVars, {'regressionScore'}, {'regressf', 'regressb'});
    regressIdx = obj.net.getLayerIndex(regressName);
    newParams = obj.net.layers(regressIdx).block.initParams();
    obj.net.params(obj.net.layers(regressIdx).paramIndexes(1)).value = newParams{1} / std(newParams{1}(:)) * 0.001; % Girshick initialization with std of 0.001
    obj.net.params(obj.net.layers(regressIdx).paramIndexes(2)).value = newParams{2};

    obj.net.addLayer('regressLoss', dagnn.LossRegress('loss', 'Smooth', 'smoothMaxDiff', 1), ...
        {'regressionScore', 'regressionTargets', 'instanceWeights'}, 'regressObjective');
end

%%% Set correct learning rates and biases (Girshick style)
if obj.nnOpts.fastRcnnParams
    % Biases have learning rate of 2 and no weight decay
    for lI = 1 : length(obj.net.layers)
        if isa(obj.net.layers(lI).block, 'dagnn.Conv')
            biasI = obj.net.layers(lI).paramIndexes(2);
            obj.net.params(biasI).learningRate = 2;
            obj.net.params(biasI).weightDecay = 0;
        end
    end

    % First convolutional layer should not learn
    % Note that this is different from Fast R-CNN, but gives good results
    % Also note that there is no speedup as Matconvnet still computes the
    % gradients, but does not update the parameters.
    conv1I = obj.net.getLayerIndex('conv1'); % AlexNet-style networks
    if isnan(conv1I)
        conv1I = obj.net.getLayerIndex('conv1_1'); % VGG-16 style networks
    end
    obj.net.params(obj.net.layers(conv1I).paramIndexes(1)).learningRate = 0;
    obj.net.params(obj.net.layers(conv1I).paramIndexes(2)).learningRate = 0;
end

author: nightrome

repo: matconvnet-calvin


text:
% StackExchange Signal Processing Q688
% https://dsp.stackexchange.com/questions/688
% What Is the Algorithm Behind Photoshop's �Black and White� Adjustment Layer?
% References:
%   1.  A
% Remarks:
%   1.  HSL Picker - http://hslpicker.com/.
%   2.  Colorizer Color Picker - http://colorizer.org/.
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     22/12/2018
%   *   First release.

%% General Parameters

subStreamNumberDefault = 2123;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;
generateImages  = OFF;

%% Simulation Parameters

vPhotoshopValues = randi([0, 100], [6, 1]);
% vPhotoshopValues = [56; 61; 50; 78; 76; 54];
vPhotoshopValues = [30; 98; 51; 73; 5; 53]; %<! This is used in Photoshop

% 11 Colors - [H, 1, 0.5] Where H goes in 30% Steps -> [0, 30, 60, ..., 330]
% Each color is a row.
mBaseColors = [1, 0, 0; 1, 0.5, 0; 1, 1, 0; 0.5, 1, 0; 0, 1, 0; 0, 1, 0.5; 0, 1, 1; 0, 0.5, 1; 0, 0, 1; 0.5, 0, 1; 1, 0, 1; 1, 0, 0.5];

mBaseColors = rand(21, 3);

cellSize    = 50; %<! Pixels

%% Generate Data

numColors = size(mBaseColors, 1);

numRows = cellSize;
numCols = cellSize * numColors;

mI = zeros(numRows, numCols, 3);

for ii = 1:numColors
    vCurrColor = mBaseColors(ii, :);
    firstColIdx = ((ii - 1) * cellSize) + 1;
    lastColIdx  = ii * cellSize;
    mI(:, firstColIdx:lastColIdx, :) = repmat(reshape(vCurrColor, [1, 1, 3]), [cellSize, cellSize 1]);
end

% vCoeffValues = [0.5; 0; 0; 0; 0; 0];
vCoeffValues = (vPhotoshopValues - 50) ./ 50;
mO = ApplyBlackWhiteFilter(mI, vCoeffValues);

% figure;
% imshow(mI);
% figure;
% imshow(mO);
%
% figure;
% imshow(im2uint8(mO));

if(generateImages == ON)
    imwrite(im2uint8(mI), 'ReferenceImage2.png');
end

%% Analysis vs. Photoshop

mORef   = im2single(imread('PhotoshopImage2.png'));
mE      = mORef - mO;

maxAbsDev = max(abs(mE(:)));

%% Display Results

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = subplot(4, 1, 1);
hImgObj = imshow(mI);
set(get(hAxes, 'Title'), 'String', {['Reference Image']}, ...
    'FontSize', fontSizeTitle);

hAxes   = subplot(4, 1, 2);
hImgObj = imshow(mO);
set(get(hAxes, 'Title'), 'String', {['MATLAB Output'], ['Black & White Adjustment Layer Input Values - ', num2str(vPhotoshopValues.')]}, ...
    'FontSize', fontSizeTitle);

hAxes   = subplot(4, 1, 3);
hImgObj = imshow(mORef);
set(get(hAxes, 'Title'), 'String', {['Photoshop Output'], ['Black & White Adjustment Layer Input Values - ', num2str(vPhotoshopValues.')]}, ...
    'FontSize', fontSizeTitle);

hAxes   = subplot(4, 1, 4);
hImgObj = imshow(mE);
set(get(hAxes, 'Title'), 'String', {['Error Image (Absolute Deviation)'], ['Max Absolute Deviation - ', num2str(maxAbsDev)]}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function newRxnScores = groupRxnScores(model, origRxnScores, origRxnIds, groupIds, origRxnsToZero)
% groupRxnScores
% This function sums up the reaction scores for all reactions that were merged
% into one by the linear merge.
%
% model          The model with linearly merged rxns.
% origRxnScores  The rxnScores from the model before the linear merge.
% origRxnIds     The rxn ids of the model before the linear merge.
% groupIds       The groupIds vector output from linearMerge.
%                There is one integer for each rxn in origRxnIds. 0 means
%                the reaction was merged. A non-zero integer means that the
%                reaction was merged with all other rxns having the same integer.
% origRxnsToZero A logical vector saying which of the original rxns that should not
%                be part of the problem. The way this is solved is that all
%                such reactions have a rxnScore of 0. If any original rxnScore
%                value should be zero (which is very unlikely) it is changed to 0.01.
%                If the sum of the rxnScores for a merged rxn becomes zero
%                while some of them are nonzero, the new value will also be 0.01,
%                to distinguish the rxn from rxns with only rxns to zero.
%                There are two reasons why we don't want zeros in the reaction
%                scores unless these reactions should be ignored:
%                  1) we want to be able to separate those
%                  2) it is difficult to handle a zero value in the MILP -
%                     the on/off of such a reaction can be random, so better
%                     to fix it in one direction.

newRxnScores = zeros(length(model.rxns),1);
[~,ia,ib] = intersect(model.rxns,origRxnIds);
grpIdsMerged = nan(length(model.rxns),1);
grpIdsMerged(ia) = groupIds(ib);
%check if any of the original scores are 0, in that case change them to 0.01 (unlikely)
origRxnScores(origRxnScores == 0) = 0.01;
%Then set the rxn scores for rxns to zero to 0
origRxnScores(origRxnsToZero) = 0;

%fill in original scores
newRxnScores(ia) = origRxnScores(ib);

for i = 1:length(model.rxns)
    %for reactions that are not merged with anything, just keep score as it is
    if grpIdsMerged(i) ~= 0
        %find all original rxns in the group
        sel = groupIds == grpIdsMerged(i);
        newRxnScores(i) = sum(origRxnScores(sel));
        if (newRxnScores(i) == 0 && any(origRxnScores(sel) ~= 0))
           %special unfortunate case, where the reactions happen to sum to 0 while some of them are nonzero
           %set to 0.01 in this case (proabably pretty unusual)
           newRxnScores(i) = 0.01;
        end
    end
end

end

author: SysBioChalmers

repo: RAVEN


text:
function cS = olivine(cs)
% olivine crystal shape

if nargin == 0
  cs = crystalSymmetry('mmm',[4.762 10.225 5.994],'mineral','olivine');
end

N = Miller({1,1,0},{1,2,0},{0,1,0},{1,0,1},{0,2,1},{0,0,1},cs);

cS = crystalShape(N,2.0,[6.0,1.0,6.0]);

end

author: mtex-toolbox

repo: mtex

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


text:
function M = globMatrixIFE3DFace(fun1,fun2,fem1,fem2,femIF,d1,j1,d2,j2)

%% USAGE: generate global matrix on Interface Faces
%             e.g.  \int_{Omega^I} [coef1 Du*n] {coef2 v} dx
%
% INPUTS:
% fun1 --- coefficient function in front of test function
% fun2 --- coefficient function in front of trial function
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
% femIF --- quadrature info on interface faces
% d1 --- derivative info for test function
% d2 --- derivative info for trial function
%            d = 0: function value
%            d = 1: normal gradient value Du*n
% j1 --- jump info for test function
% j2 --- jump info for trial function
%            j = 0: average e.g. {u} = 1/2*(uL+uR)
%            j = 1: average e.g. [u] = uL-uR
%            Note: on Dirichlet boundary: [u] = {u} = u
% OUTPUTS:
% [I J X] --- triplets of the sparse matrix.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initialization
dof1 = fem1.ldof; dof2 = fem2.ldof;
if strcmp(fem1.type,'P1')||strcmp(fem1.type,'DGP1')||strcmp(fem1.type,'CR')
    feEvalBas1 = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem1.type,'DGP2')
    feEvalBas1 = @evalP2Bas3D;
end

if strcmp(fem2.type,'P1')||strcmp(fem2.type,'DGP1')||strcmp(fem2.type,'CR')
    feEvalBas2 = @evalP1Bas3D;
elseif strcmp(fem2.type,'P2')||strcmp(fem2.type,'DGP2')
    feEvalBas2 = @evalP2Bas3D;
end

%% 1. Matrix on Internal Interface Faces
A = femIF.area; nm = femIF.normal;
gw = femIF.gw; gx = femIF.gx; gy = femIF.gy; gz = femIF.gz;
nt = size(femIF.tL,1); % NOT # of interface elements, but quadrature element
ntot = 4*nt*dof1*dof2;
I = zeros(ntot,1); J = zeros(ntot,1); X = zeros(ntot,1);

coef1 = feval(fun1,gx,gy,gz);
coef2 = feval(fun2,gx,gy,gz);

IbasL = cell(dof1,1); IbasR = cell(dof1,1);
if d1 == 0
    for i = 1:dof1
        IbasL{i} = feEvalBas1(femIF.basL(:,:,i), gx, gy, gz, [0,0,0]);
        IbasR{i} = feEvalBas1(femIF.basR(:,:,i), gx, gy, gz, [0,0,0]);
    end
elseif d1 == 1
    for i = 1:dof1
        IbasLx = feEvalBas1(femIF.basL(:,:,i), gx, gy, gz, [1,0,0]);
        IbasLy = feEvalBas1(femIF.basL(:,:,i), gx, gy, gz, [0,1,0]);
        IbasLz = feEvalBas1(femIF.basL(:,:,i), gx, gy, gz, [0,0,1]);

        IbasRx = feEvalBas1(femIF.basR(:,:,i), gx, gy, gz, [1,0,0]);
        IbasRy = feEvalBas1(femIF.basR(:,:,i), gx, gy, gz, [0,1,0]);
        IbasRz = feEvalBas1(femIF.basR(:,:,i), gx, gy, gz, [0,0,1]);

        IbasL{i} = IbasLx.*nm(:,1) + IbasLy.*nm(:,2) + IbasLz.*nm(:,3);
        IbasR{i} = IbasRx.*nm(:,1) + IbasRy.*nm(:,2) + IbasRz.*nm(:,3);
    end
end

JbasL = cell(dof2,1);  JbasR = cell(dof2,1);
if d2 == 0
    for j = 1:dof2
        JbasL{j} = feEvalBas2(femIF.basL(:,:,j), gx, gy, gz, [0,0,0]);
        JbasR{j} = feEvalBas2(femIF.basR(:,:,j), gx, gy, gz, [0,0,0]);
    end
elseif d2 == 1
    for j = 1:dof2
        JbasLx = feEvalBas2(femIF.basL(:,:,j), gx, gy, gz, [1,0,0]);
        JbasLy = feEvalBas2(femIF.basL(:,:,j), gx, gy, gz, [0,1,0]);
        JbasLz = feEvalBas2(femIF.basL(:,:,j), gx, gy, gz, [0,0,1]);
        JbasRx = feEvalBas2(femIF.basR(:,:,j), gx, gy, gz, [1,0,0]);
        JbasRy = feEvalBas2(femIF.basR(:,:,j), gx, gy, gz, [0,1,0]);
        JbasRz = feEvalBas2(femIF.basR(:,:,j), gx, gy, gz, [0,0,1]);
        JbasL{j} = JbasLx.*nm(:,1) + JbasLy.*nm(:,2) + JbasLz.*nm(:,3);
        JbasR{j} = JbasRx.*nm(:,1) + JbasRy.*nm(:,2) + JbasRz.*nm(:,3);
    end
end

% Jump and Average
if j1 == 1
    for i = 1:dof1
        IbasL{i} = IbasL{i}; IbasR{i} = -1*IbasR{i};
    end
elseif j1 == 0
    for i = 1:dof1
        IbasL{i} = 1/2*IbasL{i}; IbasR{i} = 1/2*IbasR{i};
    end
end
if j2 == 1
    for j = 1:dof2
        JbasL{j} = JbasL{j};  JbasR{j} = -1*JbasR{j};
    end
elseif j2 == 0
    for j = 1:dof2
        JbasL{j} = 1/2*JbasL{j}; JbasR{j} = 1/2*JbasR{j};
    end
end

ind = 0;
for i = 1:dof1
    for j = 1:dof2
        I(ind+1:ind+4*nt) = [femIF.tL(:,i);femIF.tL(:,i);femIF.tR(:,i);femIF.tR(:,i)];
        J(ind+1:ind+4*nt) = [femIF.tL(:,j);femIF.tR(:,j);femIF.tL(:,j);femIF.tR(:,j)];
        X(ind+1:ind+4*nt) = [...
            A.*sum((((coef1.*IbasL{i}).*(coef2.*JbasL{j})).*gw'),2); ...
            A.*sum((((coef1.*IbasL{i}).*(coef2.*JbasR{j})).*gw'),2); ...
            A.*sum((((coef1.*IbasR{i}).*(coef2.*JbasL{j})).*gw'),2); ...
            A.*sum((((coef1.*IbasR{i}).*(coef2.*JbasR{j})).*gw'),2)];
        ind = ind + 4*nt;
    end
end
M = sparse(I,J,X,size(fem1.p,1),size(fem2.p,1));

%% 2. Matrix on Boundary Interface Faces
if size(femIF.tB,1) ~= 0
    nmB = femIF.normalB;
    AB = femIF.areaB; gxB = femIF.gxB; gyB = femIF.gyB; gzB = femIF.gzB;
    ntB = size(femIF.tB,1); % not number of interface element, but quadrature element
    ntotB = ntB*dof1*dof2;
    IB = zeros(ntotB,1); JB = zeros(ntotB,1); XB = zeros(ntotB,1);

    coef1B = feval(fun1,gxB,gyB,gzB);
    coef2B = feval(fun2,gxB,gyB,gzB);
    IbasB = cell(dof1,1);
    if d1 == 0
        for i = 1:dof1
            IbasB{i} = feEvalBas1(femIF.basB(:,:,i), gxB, gyB, gzB, [0,0,0]);
        end
    elseif d1 == 1
        for i = 1:dof1
            IbasBx = feEvalBas1(femIF.basB(:,:,i), gxB, gyB, gzB, [1,0,0]);
            IbasBy = feEvalBas1(femIF.basB(:,:,i), gxB, gyB, gzB, [0,1,0]);
            IbasBz = feEvalBas1(femIF.basB(:,:,i), gxB, gyB, gzB, [0,0,1]);
            IbasB{i} = IbasBx.*nmB(:,1) + IbasBy.*nmB(:,2) + IbasBz.*nmB(:,3);
        end
    end

    JbasB = cell(dof2,1);
    if d2 == 0
        for j = 1:dof2
            JbasB{j} = feEvalBas2(femIF.basB(:,:,j), gxB, gyB, gzB, [0,0,0]);
        end
    elseif d2 == 1
        for j = 1:dof2
            JbasBx = feEvalBas2(femIF.basB(:,:,j), gxB, gyB, gzB, [1,0,0]);
            JbasBy = feEvalBas2(femIF.basB(:,:,j), gxB, gyB, gzB, [0,1,0]);
            JbasBz = feEvalBas2(femIF.basB(:,:,j), gxB, gyB, gzB, [0,0,1]);
            JbasB{j} = JbasBx.*nmB(:,1) + JbasBy.*nmB(:,2) + JbasBz.*nmB(:,3);
        end
    end

    ind = 0;
    for i = 1:dof1
        for j = 1:dof2
            IB(ind+1:ind+ntB) = femIF.tB(:,i);
            JB(ind+1:ind+ntB) = femIF.tB(:,j);
            XB(ind+1:ind+ntB) = AB.*sum((((coef1B.*IbasB{i}).*(coef2B.*JbasB{j})).*gw'),2);
            ind = ind + ntB;
        end
    end
    MB = sparse(IB,JB,XB,size(fem1.p,1),size(fem2.p,1));
    M = M + MB;
end

repo: ifem


text:
function [g,h] = compute_spline_filter(s,N)

% compute_spline_filter - compute the spline filter of a wavelete transform
%
% function [g,h] = compute_spline_filter(s,N);
%
%    'N' is the length of the filter (must be even).
%    's' is the order of the spline.
%    'g' is the low pass filter.
%    'h' is the high pass filter.
%

x = (0:N-1)*2*pi/(N-1);
gg = ; % compute here the fourier transform
g = ifft(gg);
g = [g(N/2+1:end),g(1:N/2)];
h = (-1).^(0:N-1) .* g;

author: gpeyre

repo: matlab-toolboxes


text:
% Local Regression and Likelihood, Figure 4.4.
% Author: Catherine Loader
%
% AIC plot for a local poisson regression.

load mine;
a = (0.4:0.05:1)';
figure('Name','fig4_4: AIC plot for local poisson regression' );
aicplot(a,extrp,frac,'family','poisson','deg',1);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [Population,Fitness,Next] = EnvironmentalSelection(Population,N,isOrigin)
% The environmental selection of SPEA2

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Kangjia Qiao

    %% Calculate the fitness of each solution
    if isOrigin==1
        Fitness = CalFitness(Population.objs,Population.cons);
    else
        Fitness = CalFitness(Population.objs);
    end

    %% Environmental selection
    Next = Fitness < 1;
    if sum(Next) < N
        [~,Rank] = sort(Fitness);
        Next(Rank(1:N)) = true;
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
    end

    % Population for next generation
    Population = Population(Next);
    Fitness    = Fitness(Next);
    % Sort the population
    [Fitness,rank] = sort(Fitness);
    Population = Population(rank);
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end


text:
function fp = hlp_fingerprint(data,hashed)
% Make a fingerprint (hash) of the given data structure.
% Fingerprint = hlp_fingerprint(Data)
%
% This includes all contents; however, large arrays (such as EEG.data) are only spot-checked. For
% thorough checking, use hlp_cryptohash.
%
% In:
%   Data   : some data structure
%
%   Hashed : whether to return a hash of the fingerprint instead of the full fingerprint string
%            (default: true)
%
% Out:
%   Fingerprint : an integer that identifies the data (or string if Hashed=false)
%
% Notes:
%   The fingerprint is not unique and identifies the data set only with a certain (albeit high)
%   probability.
%
%   On MATLAB versions prior to 2008b, hlp_fingerprint cannot be used concurrently from timers,
%   and also may alter the random generator's state if cancelled via Ctrl+C.
%
% Examples:
%   % calculate the hash of a large data structure
%   hash = hlp_fingerprint(data);
%
% See also:
%   hlp_cryptohash
%
%                                   Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                   2010-04-02

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

if nargin < 2
    hashed = true; end

warning off MATLAB:structOnObject

if hlp_matlab_version >= 707
    fp = fingerprint(data,RandStream('swb2712','Seed',5183),hashed);
else
    try
        % save & override random state
        randstate = rand('state'); %#ok<*RAND>
        rand('state',5183);
        % make the fingerprint
        fp = fingerprint(data,0,hashed);
        % restore random state
        rand('state',randstate);
    catch e
        % restore random state in case of an error...
        rand('state',randstate);
        rethrow(e);
    end
end

% make a fingerprint of the given data structure
function fp = fingerprint(data,rs,hashed)
% convert data into a string representation
data = summarize(data,rs);
% make sure that it does not contain 0's
data(data==0) = 'x';
if hashed
    % obtain a hash code via Java (MATLAB does not support proper integer arithmetic...)
    str = java.lang.String(data);
    fp = str.hashCode()+2^31;
else
    fp = data;
end

% get a recursive string summary of arbitrary data
function x = summarize(x,rs)
if isnumeric(x)
    % numeric array
    if ~isreal(x)
        x = [real(x) imag(x)]; end
    if issparse(x)
        try
            x = [find(x) nonzeros(x)];
        catch
            x = [find(x)' nonzeros(x)];
        end
    end
    if numel(x) <= 4096
        % small matrices are hashed completely
        try
            x = ['n' typecast([size(x) x(:)'],'uint8')];
        catch
            if hlp_matlab_version <= 702
                x = ['n' typecast([size(x) double(x(:))'],'uint8')]; end
        end
    else
        % large matrices are spot-checked
        ne = numel(x);
        count = floor(256 + (ne-256)/1000);
        if hlp_matlab_version < 707
            indices = 1+floor((ne-1)*rand(1,count));
        else
            indices = 1+floor((ne-1)*rand(rs,1,count));
        end
        if size(x,2) == 1
            % x is a column vector: reindexed expression needs to be transposed
            x = ['n' typecast([size(x) x(indices)'],'uint8')];
        else
            % x is a matrix or row vector: shape follows that of indices
            x = ['n' typecast([size(x) x(indices)],'uint8')];
        end
    end
elseif iscell(x)
    % cell array
    sizeprod = cellfun('prodofsize',x(:));
    if all(sizeprod <= 1) && any(sizeprod)
        % all scalar elements (some empty, but not all)
        if all(cellfun('isclass',x(:),'double')) || all(cellfun('isclass',x(:),'single'))
            % standard floating-point scalars
            if cellfun('isreal',x)
                % all real
                x = ['cdr' typecast([size(x) x{:}],'uint8')];
            else
                % some complex
                x = ['cdc' typecast([size(x) real([x{:}]) imag([x{:}])],'uint8')];
            end
        elseif cellfun('isclass',x(:),'logical')
            % all logical
            x = ['cl' typecast(uint32(size(x)),'uint8') uint8([x{:}])];
        elseif cellfun('isclass',x(:),'char')
            % all single chars
            x = ['ccs' typecast(uint32(size(x)),'uint8') x{:}];
        else
            % generic types (structs, cells, integers, handles, ...)
            tmp = cellfun(@summarize,x,repmat({rs},size(x)),'UniformOutput',false);
            x = ['cg' typecast(uint32(size(x)),'uint8') tmp{:}];
        end
    elseif isempty(x)
        % empty cell array
        x = ['ce' typecast(uint32(size(x)),'uint8')];
    else
        % some non-scalar elements
        dims = cellfun('ndims',x(:));
        size1 = cellfun('size',x(:),1);
        size2 = cellfun('size',x(:),2);
        if all((size1+size2 == 0) & (dims == 2))
            % all empty and nondegenerate elements
            if all(cellfun('isclass',x(:),'double'))
                % []'s
                x = ['ced' typecast(uint32(size(x)),'uint8')];
            elseif all(cellfun('isclass',x(:),'cell'))
                % {}'s
                x = ['cec' typecast(uint32(size(x)),'uint8')];
            elseif all(cellfun('isclass',x(:),'struct'))
                % struct()'s
                x = ['ces' typecast(uint32(size(x)),'uint8')];
            elseif length(unique(cellfun(@class,x(:),'UniformOutput',false))) == 1
                % same class
                x = ['cex' class(x{1}) typecast(uint32(size(x)),'uint8')];
            else
                % arbitrary class...
                tmp = cellfun(@summarize,x,repmat({rs},size(x)),'UniformOutput',false);
                x = ['cg' typecast(uint32(size(x)),'uint8') tmp{:}];
            end
        elseif all((cellfun('isclass',x(:),'char') & size1 <= 1) | (sizeprod==0 & cellfun('isclass',x(:),'double')))
            % all horizontal strings or proper empty strings, possibly some []'s
            x = ['cch' [x{:}] typecast(uint32(size2'),'uint8')];
        else
            % arbitrary sizes...
            if all(cellfun('isclass',x(:),'double')) || all(cellfun('isclass',x(:),'single'))
                % all standard floating-point types...
                tmp = cellfun(@vectorize,x,'UniformOutput',false);
                % treat as a big vector...
                x = ['cn' typecast(uint32(size(x)),'uint8')  summarize([tmp{:}],rs)];
            else
                tmp = cellfun(@summarize,x,repmat({rs},size(x)),'UniformOutput',false);
                x = ['cg' typecast(uint32(size(x)),'uint8') tmp{:}];
            end
        end
    end
elseif ischar(x)
    % char array
    x = ['c' x(:)'];
elseif isstruct(x)
    % struct
    fn = fieldnames(x)';
    if numel(x) > length(fn)
        % summarize over struct fields to expose homogeneity
        x = cellfun(@(f)summarize({x.(f)},rs),fn,'UniformOutput',false);
        x = ['s' [fn{:}] ':' [x{:}]];
    else
        % summarize over struct elements
        x = ['s' [fn{:}] ':' summarize(struct2cell(x),rs)];
    end
elseif islogical(x)
    % logical array
    x = ['l' typecast(uint32(size(x)),'uint8') uint8(x(:)')];
elseif isa(x,'function_handle')
    if strncmp(char(x),'@(',2)
        f = functions(x);
        x = ['f' f.function summarize(f.workspace{1},rs)];
    else
        x = ['f ' char(x)];
    end
elseif isobject(x)
    x = ['o' class(x) ':' summarize(struct(x),rs)];
else
    try
        x = ['u' class(x) ':' summarize(struct(x),rs)];
    catch
        warning('BCILAB:hlp_fingerprint:unsupported_type','Unsupported type: %s',class(x));
        error; %#ok<LTARG>
    end
end

function x = vectorize(x)
x = x(:)';

repo: REST


text:
function [data] = simpleFill2D(data, flags, fun)
% SYNTAX:
%    [data_filled] = simpleFill1D(data, flags)
%
% DESCRIPTION:
%    fill flagged data with a simple interpolation using MATLAB
%    interp1 'pchip', 'extrap'
%
% NOTE: data can be a matrix, the operation is executed column by column
%

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

if nargin < 3
        % Correlation function
        fun = @(dist) exp(-(dist).^2);
end

    x_g = 1 : size(data,2);
    y_g = 1 : size(data,1);

    [x_mesh, y_mesh] = meshgrid(x_g, y_g);

    x_p = x_mesh(:);
    x_o = x_mesh(~flags);

    y_p = y_mesh(:);
    y_o = y_mesh(~flags);

    data_o = data(~flags);

    for x = 1 : size(data,2)
        id_p = (x_p == x);
        if sum(id_p) > 0
            tmp_x_p = x_p(id_p);
            tmp_y_p = y_p(id_p);
            d = sqrt((repmat(tmp_x_p, 1, size(x_o,1)) - repmat(x_o', size(tmp_x_p, 1), 1)).^2 + ...
                (repmat(tmp_y_p, 1, size(y_o,1)) - repmat(y_o', size(tmp_y_p, 1), 1)).^2);
            %d = d / sqrt(sum(size(data).^2));
            data(id_p) = ((fun(d) * data_o) ./ sum(fun(d),2));
        end
    end
end

author: goGPS-Project


text:
function lambda = fourier_eigenvalues ( n )

%*****************************************************************************80
%
%% FOURIER_EIGENVALUES returns the eigenvalues of the FOURIER matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, complex LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  lambda(1,1) = 1.0;

  lambda(2:4:n,1) = - 1.0;
  lambda(3:4:n,1) =   i;
  lambda(4:4:n,1) =   1.0;
  lambda(5:4:n,1) = - i;

  return
end


text:
% sga
%
% This script implements the Simple Genetic Algorithm described
% in the examples section of the GA Toolbox manual.
%

% Author:     Andrew Chipperfield
% History:    23-Mar-94     file created

NIND = 40;           % Number of individuals per subpopulations
MAXGEN = 300;        % maximum Number of generations
GGAP = .9;           % Generation gap, how many new individuals are created
NVAR = 20;           % Generation gap, how many new individuals are created
PRECI = 20;          % Precision of binary representation

% Build field descriptor
   FieldD = [rep([PRECI],[1, NVAR]); rep([-512;512],[1, NVAR]);...
              rep([1; 0; 1 ;1], [1, NVAR])];

% Initialise population
   Chrom = crtbp(NIND, NVAR*PRECI);

% Reset counters
   Best = NaN*ones(MAXGEN,1);	% best in current population
   gen = 0;			% generational counter

% Evaluate initial population
   ObjV = objfun1(bs2rv(Chrom,FieldD));

% Track best individual and display convergence
   Best(gen+1) = min(ObjV);
   plot(log10(Best),'ro');xlabel('generation'); ylabel('log10(f(x))');
   text(0.5,0.95,['Best = ', num2str(Best(gen+1))],'Units','normalized');
   drawnow;

% Generational loop
   while gen < MAXGEN,

    % Assign fitness-value to entire population
       FitnV = ranking(ObjV);

    % Select individuals for breeding
       SelCh = select('sus', Chrom, FitnV, GGAP);

    % Recombine selected individuals (crossover)
       SelCh = recombin('xovsp',SelCh,0.7);

    % Perform mutation on offspring
       SelCh = mut(SelCh);

    % Evaluate offspring, call objective function
       ObjVSel = objfun1(bs2rv(SelCh,FieldD));

    % Reinsert offspring into current population
       [Chrom ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);

    % Increment generational counter
       gen = gen+1;

    % Update display and record current best individual
       Best(gen+1) = min(ObjV);
       plot(log10(Best),'ro'); xlabel('generation'); ylabel('log10(f(x))');
       text(0.5,0.95,['Best = ', num2str(Best(gen+1))],'Units','normalized');
       drawnow;
   end
% End of GA

author: sfvsfv

repo: Mathematical-modeling


text:
function [xVals, yVals, zVals] = getScanXYZVals(scanStruct, slice)
%"getScanXYZVals"
%   Returns the x, y, and z values of the cols, rows, and slices of the
%   passed scanStruct.  If specified, slice determines which element of
%   scanInfo is used to determine the x and y vals.  If slice is invalid or
%   not specified, 1 is used.
%
%   REMINDER: These x,y,z values are the coordinates of the MIDDLE of the
%   voxels of the scan.  They are not the coordinates of the dividers
%   between the voxels.
%
%   By JRA 12/26/03
%
%   scanStruct      : ONE planC{indexS.scan} struct.
%   slice           : slice to base x,y vals on
%
% xVals yVals zVals : x,y,z Values for scan.
%
% Usage:
%   function [xVals, yVals, zVals] = getScanXYZVals(scanStruct, slice)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if exist('slice','var')
    scanInfo = scanStruct.scanInfo(slice);
else
    scanInfo = scanStruct.scanInfo(1);
end

sizeDim1 = scanInfo.sizeOfDimension1-1;
sizeDim2 = scanInfo.sizeOfDimension2-1;

MATLABVer = version;
if MATLABVer(1) ~= '6'
    sizeDim1 = double(sizeDim1);
    sizeDim2 = double(sizeDim2);
end

xVals = scanInfo.xOffset - (sizeDim2*scanInfo.grid2Units)/2 : scanInfo.grid2Units : scanInfo.xOffset + (sizeDim2*scanInfo.grid2Units)/2;
yVals = fliplr(scanInfo.yOffset - (sizeDim1*scanInfo.grid1Units)/2 : scanInfo.grid1Units : scanInfo.yOffset + (sizeDim1*scanInfo.grid1Units)/2);
sI = scanStruct.scanInfo;
zVals = [sI.zValue];


text:
function Es = PQspreadCB (E, Ver)
% Spread an excitation vector (pitch pattern) - FFT model
% Both E and Es are powers

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:32:58 $

persistent Bs Version

if (~ strcmp (Ver, Version))
    Version = Ver;
    Nc = length (E);
    Bs = PQ_SpreadCB (ones(1,Nc), ones(1,Nc), Version);
end

Es = PQ_SpreadCB (E, Bs, Version);

%-------------------------
function Es = PQ_SpreadCB (E, Bs, Ver);

persistent Nc dz fc aL aUC Version

% Power law for addition of spreading
e = 0.4;

if (~ strcmp (Ver, Version))
    Version  = Ver;
    [Nc, fc, fl, fu, dz] = PQCB (Version);
end

% Allocate storage
aUCEe = zeros (1, Nc);
Ene = zeros (1, Nc);
Es = zeros (1, Nc);

% Calculate energy dependent terms
aL = 10^(-2.7 * dz);
for (m = 0:Nc-1)
    aUC = 10^((-2.4 - 23 / fc(m+1)) * dz);
    aUCE = aUC * E(m+1)^(0.2 * dz);
    gIL = (1 - aL^(m+1)) / (1 - aL);
    gIU = (1 - aUCE^(Nc-m)) / (1 - aUCE);
    En = E(m+1) / (gIL + gIU - 1);
    aUCEe(m+1) = aUCE^e;
    Ene(m+1) = En^e;
end

% Lower spreading
Es(Nc-1+1) = Ene(Nc-1+1);
aLe = aL^e;
for (m = Nc-2:-1:0)
    Es(m+1) = aLe * Es(m+1+1) + Ene(m+1);
end

% Upper spreading i > m
for (m = 0:Nc-2)
    r = Ene(m+1);
    a = aUCEe(m+1);
    for (i = m+1:Nc-1)
       r = r * a;
       Es(i+1) = Es(i+1) + r;
    end
end

for (i = 0:Nc-1)
    Es(i+1) = (Es(i+1))^(1/e) / Bs(i+1);
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function t = toeplitz( c, r )

%   Disciplined convex/geometric programming information for TOEPLITZ:
%      TOEPLITZ imposes no convexity restrictions on its arguments.
%      Instead of using the TOEPLITZ function, however, consider
%      creating a matrix variable using the 'toeplitz' keyword; e.g.
%          variable X(5,5) toeplitz;

%
% Check arguments
%

if nargin < 2,
    c    = vec( c );
    m    = length( c );
    p    = m;
    x    = [ cvx_subsref( c, p : -1 : 1 ) ; conj( cvx_subsref( c, 2 : p ) ) ];
else
    temp = cvx_subsref( r, 1 ) - cvx_subsref( c, 1 );
    if ~cvx_isconstant( temp ) || cvx_constant( temp ) ~= 0,
        warning('MATLAB:toeplitz:DiagonalConflict',['First element of ' ...
               'input column does not match first element of input row. ' ...
               '\n         Column wins diagonal conflict.'])
    end
    r = vec( r );
    c = vec( c );
    p = length( r );
    m = length( c );
    x = [ cvx_subsref( r, p : -1 : 2 ) ; c ];
end

%
% Construct matrix
%

cidx = ( 0 : m - 1 )';
ridx = p : -1 : 1;
t    = cidx( :, ones( p, 1 ) ) + ridx( ones( m, 1 ) , : );
t    = reshape( cvx_subsref( x, t ), size( t ) );

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% DTFT properties

% time difference

syms w
n=0:4;
x=[1,1,1,1,0];
x1=[0,1,1,1,1];
L=x-x1;
Left=sum(L.*exp(-j*w*n));
subplot(211);
ezplot(abs(Left));
title('| DTFT of x[n]-x[n-1] |')
w1=-2*pi:.1:2*pi;
Left=subs(Left,w,w1);
subplot(212);
plot(w1,angle(Left));
title('Phase of DTFT of x[n]-x[n-1]')
xlim([-2*pi  2*pi])

figure
X= sum(x.*exp(-j*w*n));
Right=(1-exp(-j*w))*X;
subplot(211);
ezplot(abs(Right));
title('| (1-exp(jw))*X(w) |')
subplot(212);
Right=subs(Right,w,w1);
plot(w1,angle(Right));
title('\angle (1-exp(jw))*X(w)')
xlim([-2*pi 2*pi])

author: Sable

repo: mcbench-benchmarks


text:
%COMPUTEINITIALPOTENTIALS Sets up the cliques in the clique tree that is
%passed in as a parameter.

%   P = COMPUTEINITIALPOTENTIALS(C) Takes the clique tree C which is a
%   struct with three fields:
%   - nodes: represents the cliques in the tree.
%   - edges: represents the adjacency matrix of the tree.
%   - factorList: represents the list of factors that were used to build
%   the tree.
%
%   It returns a compact form of a clique tree P that we will use through
%   the rest of the assigment. P is struct with two fields:
%   - cliqueList: represents an array of cliques with appropriate factors
%   from factorList assigned to each clique. Where the .val of each clique
%   is initialized to the initial potential of that clique.
%   - edges: represents the adjacency matrix of the tree.

% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function P = ComputeInitialPotentials(C)

% number of cliques
N = length(C.nodes);

% initialize cluster potentials
P.cliqueList = repmat(struct('var', [], 'card', [], 'val', []), N, 1);

% Find the values and cards for each clique, compute the initial potential
% for each clique by find the joint distribution over the factors in it,
% and check that the clique tree satisfies family preservation
assignments = zeros(length(C.factorList), 1);
for i = 1:length(C.factorList)
    % Assign each factor to a clique
    for j = 1:N
        % Iterate through nodes until a node with all of the variables in
        % the clique has been found
        if isempty(setdiff(C.factorList(i).var, C.nodes{j}))
            % A clique with the variables has been indentified
            assignments(i) = j;
            break
        end
    end
end
% Determine the cardinaliteis for each node in the clique that will be
% created
cardinalities = zeros(1, length(unique([C.factorList(:).var])));
for i = 1:length(cardinalities)
    % Iterate through variables
    for j = 1:length(C.factorList)
        % Iterate through factors
        if ~isempty(find(C.factorList(j).var == i))
            % If the current variable is in the factor, find its
            % cardinality
            cardinalities(i) = C.factorList(j).card(find(C.factorList(j).var == i));
            break;
        end
    end
end
for i = 1:N
    % Iterate through nodes and put each node into the clique list
    P.cliqueList(i).var = C.nodes{i};
    P.cliqueList(i).card = cardinalities(C.nodes{i});
    P.cliqueList(i).val = ones(1, prod(P.cliqueList(i).card));
end
for i = 1:length(assignments)
    % Iterate through assignments and multiply the appropriate factors
    % in to the clique
    P.cliqueList(assignments(i)) = FactorProduct(P.cliqueList(assignments(i)), C.factorList(i));
end
P.edges = C.edges;
for i = 1:length(C.factorList)
    % Check that every factor has a clique
    factorVars = C.factorList(i).var;
    allVarsInFactor = 0;
    for j = 1:N
        % Find the clique with the variables from the factor
        varNotFound = 0;
        for k = 1:length(factorVars)
            % Find the clique with all of the factorVariables
            if isempty(find(ismember(P.cliqueList(j).var, factorVars(k))))
                % The variable is not in the factor
                varNotFound = 1;
                break
            end
        end
        if varNotFound == 0
            % All of the variables are in the current factor
            allVarsInFactor = 1;
            break
        end
    end
    if allVarsInFactor == 0
        % Family intersection property has been violated
        disp('Family intersection property violated!');
    end
end

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function ARDL = ARDLmodel(ENDO,nlag,const,EXOG,nlag_ex)
% =======================================================================
% Estimate ARDL models with OLS
% =======================================================================
% ARDL = ARDLmodel(ENDO,nlag,const,EXOG,nlag_ex)
% -----------------------------------------------------------------------
% INPUT
%	- ENDO: an (nobs x 1) vector of endogenous
%	- nlag: lag length
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%	- const: 0 no constant; 1 constant; 2 constant and trend; 3 constant,
%       trend, and trend^2 [dflt = 0]
%	- EXOG: optional vector of exogenous variable (nobs x 1)
%	- nlag_ex: number of lags for exogeonus variable [dflt = 0]
% -----------------------------------------------------------------------
% OUTPUT
%   - VAR: structure including VAR estimation results
%   - VARopt: structure including VAR options (see VARoption)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Check inputs
%===============================================
[nobs, ~] = size(ENDO);
ARDL.ENDO = ENDO;
ARDL.nlag = nlag;

% Check if ther are constant, trend, both, or none
if ~exist('const','var')
    const = 1;
end
ARDL.const = const;

% Check if there is exogenous variable
if exist('EXOG','var')
    [nobs_ex, nvar_ex] = size(EXOG);
    % Check that ENDO and EXOG are conformable
    if (nobs_ex ~= nobs)
        error('var: nobs in EXOG-matrix not the same as y-matrix');
    end
    clear nobs_ex
    % Check if there is lag order of EXOG, otherwise set it to 0
    if ~exist('nlag_ex','var')
        nlag_ex = 0;
    end
    ARDL.EXOG = EXOG;
else
    nvar_ex = 0;
    nlag_ex = 0;
    ARDL.EXOG = [];
end

%% Save some parameters and create data matrices
%===============================================
    nobse        = nobs - max(nlag,nlag_ex);
    ARDL.nobs    = nobse;
    ARDL.nlag    = nlag;
    ARDL.nlag_ex = nlag_ex;
    ncoeff       = nlag;
    ARDL.ncoeff  = ncoeff;
    ncoeff_ex    = nvar_ex + nvar_ex*nlag_ex;
    nvar         = ncoeff + ncoeff_ex + const;
    ARDL.nvar    = nvar;
    ARDL.nvar_ex = nvar_ex;
    ARDL.const   = const;

% Create independent vector and lagged dependent matrix
[Y, X] = VARmakexy(ENDO,nlag,const);

% Create (lagged) exogenous matrix
if exist('EXOG','var')
    X_EX  = VARmakelags(EXOG,nlag_ex);
    if nlag == nlag_ex
        X = [X X_EX];
    elseif nlag > nlag_ex
        diff = nlag - nlag_ex;
        X_EX = X_EX(diff+1:end,:);
        X = [X X_EX];
    elseif nlag < nlag_ex
        diff = nlag_ex - nlag;
        Y = Y(diff+1:end,:);
        X = [X(diff+1:end,:) X_EX];
    end
end

ARDL.meth = 'ols';
ARDL.y = Y;
ARDL.x = X;

% xpxi = (X'X)^(-1)
if nobse < 10000
  [~, r] = qr(X,0);
  xpxi = (r'*r)\eye(nvar);
else
  xpxi = (X'*X)\eye(nvar);
end;

% OLS estimator
beta = xpxi*(X'*Y);
ARDL.beta = beta;

% Predicted values & residuals
ARDL.yhat = X*ARDL.beta;
ARDL.resid = Y - ARDL.yhat;

% Covariance matrix of residuals
sigu = ARDL.resid'*ARDL.resid;
ARDL.sige = sigu/(nobse-nvar);

% Covariance matrix of beta
sigbeta = ARDL.sige*xpxi;
ARDL.sigbeta = sigbeta;

% Std errors of beta, t-stats, and intervals
tmp = (ARDL.sige)*(diag(xpxi));
bstd = sqrt(tmp);
ARDL.bstd = bstd;
tcrit=-tdis_inv(.025,nobse);
ARDL.bint=[ARDL.beta-tcrit.*bstd, ARDL.beta+tcrit.*bstd];
ARDL.tstat = ARDL.beta./(sqrt(tmp));
ARDL.tprob = tdis_prb(ARDL.tstat,nobs);

% R2
ym = Y - mean(Y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
ARDL.rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(nobse-nvar);
rsqr2 = rsqr2/(nobse-1.0);
if rsqr2 ~= 0
    ARDL.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
else
    ARDL.rbar = ARDL.rsqr;
end;

% Durbin-Watson
ediff = ARDL.resid(2:nobse) - ARDL.resid(1:nobse-1);
ARDL.dw = (ediff'*ediff)/sigu; % durbin-watson
ARDL.const = const;

% F-test
if const>0
    fx = X(:,1);
    fxpxi = (fx'*fx)\eye(1);
    fbeta = fxpxi*(fx'*Y);
    fyhat = fx*fbeta;
    fresid = Y - fyhat;
    fsigu = fresid'*fresid;
    fym = Y - mean(Y);
    frsqr1 = fsigu;
    frsqr2 = fym'*fym;
    frsqr = 1.0 - frsqr1/frsqr2; % r-squared
    ARDL.F = ((frsqr-ARDL.rsqr)/(1-nvar)) / ((1-ARDL.rsqr)/(nobse-nvar));
end

% % Long-run coefficients
% q = ncoeff_ex;
% p = ncoeff;
%
% sumendo = sum(beta(const+1:const+p)); % sum of lagged endo
% sumexog = sum(beta(const+p+1:end)); % sum of cont and lagged exog
%
% theta = sumexog/(1-sumendo);
% ARDL.theta = theta;
%
% aux1(1:p,1) = sumexog/((1-sumendo)^2);
% aux2(1:q,1) = 1/(1-sumendo);
% dtheta = [aux1; aux2];
% sigbeta_noconst = sigbeta(const+1:nvar,const+1:nvar);
% ARDL.sigtheta = dtheta'*sigbeta_noconst*dtheta;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function val=numInvolutions(n)
%%NUMINVOLUTIONS Determine the number of involutions of length n. An
%                involution is a permution of items 1:n such that the
%                permutation is its own inverse. The number of involutions
%                also equals the number of possible Young tableau having n
%                cells.
%
%INPUTS: n A positive integer >=0.
%
%OUTPUTS: val The number of involutions of length n. For n=0, this is
%             defined as 1.
%
%This function uses the recurrence relation given in Equation 40 in Section
%5.1.4 of [1].
%
%REFERENCES:
%[1] D. Knuth, The Art of Computer Programming: Sorting and Searching, 2nd
%    ed. Boston: Addison-Wesley, 2018, vol. 3.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(n<0)
    error('n must be >=0.')
end

if(n<=1)
    val=1;
    return;
end

an1=1;
an=1;
for k=2:n
    an2=an1;
    an1=an;
    an=an1+(k-1)*an2;
end
val=an;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [us] = min2us(min)
% Convert time from minutes to microseconds.
% Chad Greene 2012
us = min*60000000;


text:
function targets = sledgeset2adjlist(n, edges, sch)
%SLEDGESET2ADJLIST Converts edge set to adjacency list
%
% $ Syntax $
%   - targets = sledgeset2adjlist(n, edges, sch)
%
% $ Arguments $
%   - n:        The number of (source) nodes
%   - edges:    The set of edges (nedges x 2 or nedges x 3)
%   - sch:      The id of scheme of conversion to take
%               - 0:    no value -> no value
%               - 1:    no value -> has value
%               - 2:    has value -> no value
%               - 3:    has value -> has value
%   - targets:  The cell array of targets in adj list
%
% $ Remarks $
%   - an internal function for graph representation conversion.
%     no checking of input arguments would be performed.
%
% $ History $
%   - Created by Dahua Lin, on Sep 9, 2006
%

%% prepare storage

targets = cell(n, 1);

%% sort edges

I = edges(:, 1);
J = edges(:, 2);

[I, si] = sort(I);
J = J(si);

if sch == 3
    V = edges(si, 3);
elseif sch == 1
    nedges = length(I);
    V = ones(nedges, 1);
end

%% group targets

nums = slcount(I);
ni = length(nums);
[spos, epos] = slnums2bounds(nums);

if sch == 1 || sch == 3
    for i = 1 : ni
        s = spos(i); e = epos(i);
        if s <= e
            targets{I(s)} = [J(s:e), V(s:e)];
        end
    end
else
    for i = 1 : ni
        s = spos(i); e = epos(i);
        if s <= e
            targets{I(s)} = J(s:e);
        end
    end
end

author: lmthang


text:
%ISCALESPACE  Scale-space image sequence
%
% [G,L,S] = ISCALESPACE(IM, N, SIGMA) is a scale space image sequence of
% length N derived from IM (HxW).  The standard deviation of the smoothing
% Gaussian is SIGMA.  At each scale step the variance of the Gaussian increases
% by SIGMA^2.  The first step in the sequence is the original image.
%
% G (HxWxN) is the scale sequence, L (HxWxN) is the absolute value of the
% Laplacian of Gaussian (LoG) of the scale sequence, corresponding to each
% step of the sequence, and S (Nx1) is the vector of scales.
%
% [G,L,S] = ISCALESPACE(IM, N) as above but SIGMA=1.
%
% Examples::
%  Create a scale-space image sequence
%         im = iread('lena.png', 'double', 'grey');
%         [G,L,s] = iscalespace(im, 50, 2);
%  Then find scale-space maxima, an array of ScalePointFeature objects.
%         f = iscalemax(L, s);
%  Look at the scalespace volume
%         slice(L, [], [], 5:10:50); shading interp
%
% Notes::
% - The Laplacian is approximated by the the difference of adjacent Gaussians.
%
% See also ISCALEMAX, ISMOOTH, ILAPLACE, KLOG.

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

function [G, L, S] = scalespace(im, n, sigma)

    if nargin < 3
        sigma = 1;
    end

    g(:,:,1) = im;
    scale = 0.5;
    scales = scale;
    for i=1:n-1,
        im = ismooth(im, sigma);
        scale = sqrt(scale^2 + sigma^2);
        scales = [scales; scale];
        g(:,:,i+1) = im;
        lap(:,:,i) = scale^2 * ( g(:,:,i+1) - g(:,:,i) );
    end

    % return results as requested
    if nargout > 2,
        S = scales;
    end
    if nargout > 0,
        G = g;
    end
    if nargout > 1,
        L = lap;
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function cl = atlas_spmanatomy_get_cluster(region_name)
% cl = atlas_spmanatomy_get_cluster(region_name)
%
% Returns a REGION object, cl, so that a named area from the SPM anatomy
% toolbox can be visualized.
%
% Requires SPM ANATOMY TOOLBOX V1.7 by Eickhoff (see code to update this.)
%
% For a list of valid region names, do not enter region_name:
% atlas_spmanatomy_get_cluster
%
% Example:
% ba44_cl = atlas_spmanatomy_get_cluster('Area 44');
%
% Visualize:
% cluster_orthviews(ba44_cl, {[0 0 1]}, 'add');
%
% Tor wager, Jan 2011

cl = region;

basename = 'AllAreas_v17_MPM.mat';
filename = which(basename);

if ~exist(filename, 'file')
    fprintf('Looking for atlas filename %s\nCannot find it.\n', basename)
    error('Exiting.')
end

load(filename)

imagename = which('AllAreas_v17.img');

if ~exist(imagename, 'file')
    fprintf('Looking for %s\nCannot find it.\n', imagename)
    error('Exiting.')
end

volInfo = iimg_read_img(imagename, 2);

names = char(MAP.name);

if nargin < 1 || isempty(region_name)
   disp('Run this function and enter a region name:')
   disp(names);
   return

end

%% Get info about chosen region

wh = strmatch(region_name, names, 'exact');

if isempty(wh)
    disp('Warning! No match.')
    return
end

info = MAP(wh);

% Transfer info
N = fieldnames(cl);
M = fieldnames(info);
for i = 1:length(M)
    iname = M{i};
    if any(strmatch(iname, N))
        cl.(iname) = info.(iname);
    end
end

% Add other info

cl.M = volInfo.mat;
cl.dim = volInfo.dim;
cl.voxSize = diag(volInfo.mat(1:3, 1:3));
cl.center = mean(cl.XYZ', 1);
cl.mm_center = mean(cl.XYZmm', 1);
cl.source_images = imagename;
cl.numVox = size(cl.XYZ, 2);
cl.shorttitle = sprintf('%s', region_name);
cl.title = sprintf('%s from %s', region_name, filename);
cl.descrip1 = sprintf('%s, entry %3.0f in %s, mapped using %s', region_name, wh, filename, imagename);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function timestamps = edgeDetect(signal, threshold, type)

% edgeDetect
%
% Detects rising or falling edges of a signal whenever the signal crosses
% the threshold.
%
%
% Use timestamps = edgeDetect(signal, threshold, type)
%
% Variables fname and threshold are optional.
%
%   signal:       Name of the file to be opened. If the fname is omitted
%                 the user will be prompted to select a file.
%
%   threshold:    The value used for threshold crossings. This function
%                 will detect when the signal crosses this value.
%                 DEFAULT: It will automatically choose the threshold at
%                 80% of the maximum value in the signal.
%
%
%   type:         Contains the points in the signal where the threshold
%                 crossing occurs.
%                 DEFAULT: It will find the rising edges.
%
%   Example 1:
%   timestamps = edgeDetect(signal, 1000);
%
%   In the example above, the points where signal crosses value 1000
%   (rising edge) is detected and returned in variable timestamps.
%
%   timestamps = edgeDetect(signal, 1000, 'falling');
%
%   In the example above, the points where signal crosses value 1000
%   (falling edge) is detected and returned in variable timestamps.
%
%   Kian Torab
%   kian@blackrockmicro.com
%   Blackrock Microsystems
%   Version 1.1.0.0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0:
%   - Initial release.
%
% 1.1.0.0:
%   - Added automatic edge detection.
%   - Updated help.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Validations

% Validating input arguments.
if nargin < 1 || nargin > 3
    disp('Invalid number of input arguments. Use ''help edgeDetect'' for more information.');
    return;
end

if ~exist('threshold', 'var')
    threshold = 0.8 * max(signal);
    disp(['Threshold was not provided. It wss automatically calculated and set at ' num2str(threshold) '.']);
end

% Validating variable 'type'
if ~exist('type', 'var')
    type = 'rising';
end

% Validating type and determining the threshold crossing points
if strcmpi(type, 'rising')
    timestamps = signal>threshold;
elseif strcmpi(type, 'falling')
    timestamps = signal<threshold;
else
    disp('Type does not exist. Please type ''help edgeDetect'' to see all available types.');
    timestamps = 0;
    return;
end

% Finding all the points where the signal crosses the threshold
timestamps = diff(timestamps);
timestamps = find(timestamps==1);

author: fieldtrip

repo: fieldtrip


text:
function params = IEEE80211a_settings(...
    symbolPeriod, ...
    OFDMSymPerFrame, ...
    OFDMTrainPerFrame, ...
    vtbd);

% Note: A sample period of 0.08e-6 corresponds to approx. 12 Msym/s,
% which, in turn, corresponds to 54 Mb/s for 3/4-coded 64-QAM

% IEEE 802.11a - fundamental sizes
p.NSD = 48;  % number of data symbols in OFDM symbol (variable name as in 802.11a standard)
p.NST = 52;  % number of data symbols and pilots in OFDM symbol (variable name as in 802.11a standard)
p.NFFT = 64; % number of points on FFT
p.NcyclicPrefix = 16;

p.NFFT2 = p.NFFT + p.NcyclicPrefix;

% TX indices
p.TXFFTShiftIndices = [p.NST/2+1:p.NFFT 1:p.NST/2];
p.TXCyclicPrefixIndices = [p.NFFT-[p.NcyclicPrefix-1:-1:0] 1:p.NFFT];

% RX indices
p.RXCyclicPrefixIndices = [p.NcyclicPrefix+1:p.NFFT2];
p.RXSelectFFTIndices = [p.NFFT-[p.NST/2-1:-1:0] 1:p.NST/2+1];

% OFDM symbols
p.OFDMSymPerFrame = OFDMSymPerFrame;
p.OFDMTrainPerFrame = OFDMTrainPerFrame;
p.OFDMTotSymPerFrame = OFDMSymPerFrame + OFDMTrainPerFrame;

% Constellation symbols
p.numTxSymbols = p.NSD * OFDMSymPerFrame;
p.numTrainingSymbols = p.NSD * OFDMTrainPerFrame;

% Training sequence
p.long_training_seq = ...
    [1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 0 ...
     1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1].';

% Modulator/demodulator banks
p.numModulators = 8;
p.txBitsPerSymbol = [1 1 2 2 4 4 6 6];
p.txBitsPerBlock = p.numTxSymbols * p.txBitsPerSymbol;
p.modOrder = 2.^p.txBitsPerSymbol;
p.codeRate = [1/2 3/4 1/2 3/4 1/2 3/4 2/3 3/4];
p.bitsPerBlock = p.txBitsPerBlock .* p.codeRate;
p.bitsPerSymbol = p.txBitsPerSymbol .* p.codeRate;
p.maxBitsPerBlock = max(p.bitsPerBlock);

% Frame size for variable rate source
p.nSource = min( gcd( min(p.bitsPerBlock), p.bitsPerBlock ) );

% Source blocks per TX frame
p.nS = p.bitsPerBlock/p.nSource;

% Timing-related parameters
p.symbolPeriod = symbolPeriod;
p.blockPeriod = p.numTxSymbols * symbolPeriod;
p.bitPeriod = symbolPeriod ./ p.bitsPerSymbol;
p.minBitPeriod =  min(p.bitPeriod);
p.chanSamplePeriod = p.blockPeriod/(p.OFDMTotSymPerFrame * p.NFFT2);

% Viterbi trace back depth and link delay
vtbd_set = vtbd(ones(1, p.numModulators));
p.vtbd_set = vtbd_set;
p.link_delay = vtbd_set;

params = p;


text:
function [relaxRxnBool, solutionRelax] = minCardinalityConservationRelaxationVector(S, param, printLevel)
% DC programming for solving the cardinality optimization problem
%
% .. math::
%
%    min  ~& \lambda ||x||_0 \\
%    s.t. ~& x + S^T z = 0 \\
%         ~& -\infty \leq x \leq \infty, \\
%         ~& 1 \leq z \leq 1 / \epsilon
%
% USAGE:
%
%    [relaxRxnBool, solutionRelax] = minCardinalityConservationRelaxationVector(S, param, printLevel)
%
% INPUT:
%    S:                `m` x `n` stoichiometric matrix
%
% OPTIONAL INPUTS:
%    param:           structure with:
%
%                        * param.epsilon - (getCobraSolverParams('LP', 'feasTol')*100) 1/epsilon is the largest flux expected
%                        * param.eta - (`feasTol` * 100), cutoff for mass leak/siphon
%                        * param.nonRelaxBool - (false(n, 1)), `n` x 1 boolean vector for reactions not to relax
%    printLevel:       verbose level
%
% OUTPUTS:
%    relaxRxnBool:     `n` x 1 boolean vector where true correspond to relaxation
%    solutionRelax:    structure with:
%
%                        * solutionRelax.stat - solution status
%                        * solutionRelax.x - `n` x 1 vector where nonzeros>eta correspond to relaxations
%                        * solutionRelax.z - `m` x 1 vector where positives correspond to molecular mass

[mlt,nlt]=size(S');

if ~exist('param','var') || isempty(param)
    param.epsilon=getCobraSolverParams('LP', 'feasTol')*100;
    feasTol = getCobraSolverParams('LP', 'feasTol');
    param.eta=feasTol*100;
    param.nonRelaxBool=false(mlt,1);
    param.checkConsistency=0;
else
    if ~isfield(param,'epsilon')
        param.epsilon=getCobraSolverParams('LP', 'feasTol')*100;
    end
    if ~isfield(param,'eta')
        feasTol = getCobraSolverParams('LP', 'feasTol');
        param.eta=feasTol*100;
    end
    if ~isfield(param,'nonRelaxBool')
        param.nonRelaxBool=false(mlt,1);
    end
    if ~isfield(param,'checkConsistency')
        param.checkConsistency=1;
    end
end

if ~exist('printLevel','var')
    printLevel =0;
end

if param.checkConsistency
    % Check the stoichiometric consistency of the network without relaxation by
    % solving the following linear problem
    %       min sum(m_i)
    %           s.t     S'*m = 0
    %                   m >= 1
    % where l  is is a  mx1 vector of the molecular mass of m molecular species
    cardProblem.A=S';
    cardProblem.b=zeros(size(cardProblem.A,1),1);
    cardProblem.lb=ones(size(cardProblem.A,2),1);
    cardProblem.ub=inf*ones(size(cardProblem.A,2),1);
    cardProblem.c=1*ones(size(cardProblem.A,2),1);
    cardProblem.osense=1;
    cardProblem.csense(1:size(cardProblem.A,1),1)='E';

    solutionRelax = solveCobraLP(cardProblem,'printLevel',printLevel);
    if solutionRelax.stat==1
        solutionRelax.z = solutionRelax.full;
        solutionRelax.x = S'*solutionRelax.z;
    end
    done=1;
else
    done=0;
end

%relaxation problem
if done==0
    cardProblem.p=mlt;
    cardProblem.q=0;
    cardProblem.r=nlt;
    cardProblem.c=zeros(nlt+mlt,1);
    if 1
        cardProblem.A=[speye(mlt,mlt),S'];
    else
        cardProblem.A=[sparse(mlt,mlt),S'];
    end
    cardProblem.b=zeros(mlt,1);
    cardProblem.lb=[-inf*ones(mlt,1);ones(nlt,1)];
    %cardProblem.lb=[zeros(mlt,1);epsilon*ones(nlt,1)];
    cardProblem.ub=[inf*ones(nlt,1);(1/param.epsilon)*ones(mlt,1)];
    %omits flux from this reaction - perhaps not a good way to do it.
    if any(param.nonRelaxBool)
        %prevent relaxation of specified reactions
        cardProblem.lb([param.nonRelaxBool;false(mlt,1)])=0;
        cardProblem.ub([param.nonRelaxBool;false(mlt,1)])=0;
    end
    cardProblem.csense(1:mlt,1)='E';
    cardProblem.lambda0=1;
    cardProblem.lambda1=getCobraSolverParams('LP', 'feasTol')*100;% sensitive to this value, 1e-4 works for Recon3Model.
    cardProblem.delta=0;
    solutionRelax = optimizeCardinality(cardProblem,param);
    %  problem                  Structure containing the following fields describing the problem
    %       p                   size of vector x
    %       q                   size of vector y
    %       r                   size of vector z
    %       c                   (p+q+r) x 1 linear objective function vector
    %       lambda              trade-off parameter of ||x||_0
    %       delta               trade-off parameter of ||y||_0
    %       A                   s x (p+q+r) LHS matrix
    %       b                   s x 1 RHS vector
    %       csense              s x 1 Constraint senses, a string containting the constraint sense for
    %                           each row in A ('E', equality, 'G' greater than, 'L' less than).
    %       lb                  (p+q+r) x 1 Lower bound vector
    %       ub                  (p+q+r) x 1 Upper bound vector
    %
    % OPTIONAL INPUTS
    % param                    parameters structure
    %       nbMaxIteration      stopping criteria - number maximal of iteration (Defaut value = 1000)
    %       epsilon             stopping criteria - (Defaut value = 10e-6)
    %       theta               parameter of the approximation (Defaut value = 2)
    %
    % OUTPUT
    % solution                  Structure containing the following fields
    %       x                   p x 1 solution vector
    %       y                   q x 1 solution vector
    %       z                   r x 1 solution vector
    %       stat                status
    %                           1 =  Solution found
    %                           2 =  Unbounded
    %                           0 =  Infeasible
    %                           -1=  Invalid input
end

%check optimality
if printLevel>2
    fprintf('%g%s\n',norm(solutionRelax.x + S'*solutionRelax.z),' = ||x + S''*z||')
    fprintf('%g%s\n',min(solutionRelax.z),' = min(z_i)')
    fprintf('%g%s\n',max(solutionRelax.z),' = min(z_i)')
    fprintf('%g%s\n',min(solutionRelax.x),' = min(x_i)')
    fprintf('%g%s\n',max(solutionRelax.x),' = max(x_i)')
end

if solutionRelax.stat==1
    %conserved if relaxation is below epsilon
    relaxRxnBool=abs(solutionRelax.x)>=param.eta;
    if printLevel>1
        fprintf('%g%s\n',norm(S(:,~relaxRxnBool)'*solutionRelax.z),' = ||N''*z|| (should be zero)')
    end
    if printLevel>1
        fprintf('%s\n',[int2str(nnz(relaxRxnBool)) '/' int2str(length(relaxRxnBool)) ' reactions relaxed.'])
    end
else
    disp(solutionRelax)
    error('solve for minimum cardinality of conservation relaxation vector failed')
end

author: opencobra

repo: cobratoolbox


text:
function z = Parabola(x)
    z = sum(x.^2);
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function f = times(f, g)
%.*   CHEBFUN multiplication.
%   F.*G multiplies F and G, where F and G may be CHEBFUN objects or scalars.
%   If F and/or G is array-valued, the dimensions must match.
%
% See also MTIMES.

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

% Deal with the special cases:
if ( ~isa(f, 'chebfun') )      % ??? * CHEBFUN

    % Ensure CHEBFUN is the first input:
    if ( ~g(1).isTransposed )
        f = times(g, f);
    else
        f = times(g.', f.').';
    end

elseif ( isempty(g) )          % CHEBFUN * []

    f = [];
    return

elseif ( isempty(f) )          % empty CHEBFUN * CHEBFUN

    % Nothing to do here. (Return an empty CHEBFUN as output.)
    return

elseif ( isnumeric(g) )        % CHEBFUN * double

    if ( numel(f) == 1 )
        % Array-valued case:

        % Loop over the funs:
        for k = 1:numel(f.funs)
            f.funs{k} = times(f.funs{k}, g);
        end

        % Multiply the pointValues:
        if ( numel(g) > 1 )
            f.pointValues = bsxfun(@times, f.pointValues, g);
        else
            f.pointValues = f.pointValues .* g;
        end

    else
        % Quasimatrix case:

        numCols = numel(f);
        % Promote g if required:
        if ( ~isscalar(g) )
            error('CHEBFUN:CHEBFUN:times:dim', 'Matrix dimensions must agree.');
        end
        % Loop over the columns:
        for k = 1:numCols
            f(k) = f(k).*g;
        end

    end

elseif ( ~isa(g, 'chebfun') )  % CHEBFUN * ???

    error('CHEBFUN:CHEBFUN:times:unknown', ...
          ['Undefined function ''times'' for input arguments of type ' ...
           '%s and %s.'], class(f), class(g));

else                           % CHEBFUN .* CHEBFUN

    % Check to see if one of the CHEBFUNs is transposed:
    if ( xor(f(1).isTransposed, g(1).isTransposed) )
        error('CHEBFUN:CHEBFUN:times:matdim', ...
            'Matrix dimensions must agree. (One input is transposed).');
    end

    dimCheck(f, g);

    if ( numel(f) == 1 && numel(g) == 1 )
        % CHEBFUN case:

        % If one of the two CHEBFUNs uses a PERIODICTECH reprensetation,
        % cast it to a NONPERIODICTECH.
        if ( ~isPeriodicTech(f.funs{1}) && isPeriodicTech(g.funs{1}) )
            g = chebfun(g, g.domain, 'tech', get(f.funs{1}, 'tech'));
        elseif ( isPeriodicTech(f.funs{1}) && ~isPeriodicTech(g.funs{1}) )
            f = chebfun(f, f.domain, 'tech', get(g.funs{1}, 'tech'));
        end

        % Overlap:
        [f, g] = overlap(f, g);

        % Loop over the FUNs:
        for k = 1:numel(f.funs)
            f.funs{k} = times(f.funs{k}, g.funs{k});
        end
        f.pointValues = f.pointValues .* g.pointValues;

    else
        % QUASIMATRIX case:

        % Convert to cell for simplicity
        f = cheb2cell(f);
        g = cheb2cell(g);

        % Loop over the columns:
        if ( numel(f) == 1 )
            for k = numel(g):-1:1
                h(k) = f{1} .* g{k};
            end
        elseif ( numel(g) == 1 )
            for k = numel(f):-1:1
                h(k) = f{k} .* g{1};
            end
        else % numel(f) = numel(g)
            for k = numel(f):-1:1
                h(k) = f{k} .* g{k};
            end
        end
        f = h;
    end
end

% Set small breakpoint values to zero:
f = thresholdBreakpointValues(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%*********************************************************************
%% sortA: sort columns of At{p} in ascending order according to the
%%        number of nonzero elements.
%%
%% [At,C,b,X0,Z0,permA,permZ] = sortA(blk,At,C,b,X0,Z0);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*********************************************************************

   function [At,C,Cnew,X0,Z0,permA,invpermA,permZ] = sortA(blk,At,C,Cnew,b,X0,Z0);

   global spdensity smallblkdim
%%
   numblk = size(blk,1);
   m = length(b);
   nnzA = zeros(numblk,m);
   permA = kron(ones(numblk,1),[1:m]);
   invpermA = kron(ones(numblk,1),[1:m]);
   permZ = cell(size(blk,1),1);
%%
   for p=1:size(blk,1)
      pblk = blk(p,:);
      n = sum(pblk{2});
      numblk = length(pblk{2});
      if strcmp(pblk{1},'s') & (max(pblk{2}) > smallblkdim)
         n2 = sum(pblk{2}.*pblk{2});  n22 = sum(pblk{2}.*(pblk{2}+1))/2;
         m1 = size(At{p,1},2);
         if (length(pblk{2}) == 1)
            tmp = abs(C{p}) + abs(Z0{p});
            if  (~isempty(At{p,1}))
                tmp = tmp + smat(blk(p,:),abs(At{p,1})*ones(m1,1),1);
            end
            if (nnz(tmp) < spdensity*n22);
               per = symamd(tmp);
               invper = zeros(n,1); invper(per) = [1:n];
               permZ{p} = invper;
               if (~isempty(At{p,1}))
                  isspAt = issparse(At{p,1});
                  for k = 1:m1
                     Ak = smat(pblk,At{p,1}(:,k),1);
                     At{p,1}(:,k) = svec(pblk,Ak(per,per),isspAt);
                  end
               end
               C{p}  = C{p}(per,per);
               Z0{p} = Z0{p}(per,per);
               X0{p} = X0{p}(per,per);
               Cnew{p} = Cnew{p}(per,per);
            else
               per = [];
            end
            if (length(pblk) > 2) & (~isempty(per))
               m2 = length(pblk{3});
               P = spconvert([(1:n)', per', ones(n,1)]);
               At{p,2} = P*At{p,2};
            end
         end
         if ~isempty(At{p,1}) & (mexnnz(At{p,1}) < m*n22/2)
            for k = 1:m1
               Ak = At{p,1}(:,k);
               nnzA(p,k) = length(find(abs(Ak) > eps));
            end
            [dummy,permAp] = sort(nnzA(p,1:m1));
            At{p,1}  = At{p,1}(:,permAp);
            permA(p,1:m1) = permAp;
            invpermA(p,permAp) = [1:m1];
         end
      elseif strcmp(pblk{1},'q') | strcmp(pblk{1},'l') | strcmp(pblk{1},'u');
         if ~issparse(At{p,1});
            At{p,1} = sparse(At{p,1});
         end
      end
   end
%%*********************************************************************


text:
function v = times(varargin)
% ensure that result is not S2Grid anymore

v = vector3d(times@vector3d(varargin{:}));

author: mtex-toolbox

repo: mtex