text:
function [ output ] = fista( par )
% fista for stably principal component pursuit problem
%{
par.M = M;
par.iter = 1000;
%lambda = 1/sqrt(max(size(M))); % default lambda
par.lambda_1 = 0.1;
par.lambda_2 = 0.001;
out = fista(par);
subplot(1,3,1),imagesc(M);
subplot(1,3,2),imagesc(out.L);
subplot(1,3,3),imagesc(out.S);
show_results(M,out.L,out.S,[],size(M,2),m,n);
%}
%%  Cun Mu and John Wright, Mar '14

M = par.M; % data matrix
[m,n] = size(M); d = min(m,n);
Omega = ones(m,n); % if not specified, full observation
obj_target = -inf;

if isfield(par,'Omega') Omega = par.Omega;  end
if isfield(par,'iter') iter_no = par.iter; end
if isfield(par,'obj') obj_target = par.obj; end

lambda_1 = par.lambda_1;
lambda_2 = par.lambda_2;
%display = par.display;
lips = 2;

L = zeros(m,n); S = zeros(m,n);
sv = round(d/10);
iter = 1; hist = 0;

L_hat = zeros(m,n); S_hat = zeros(m,n); t = 1;
L_pre = zeros(m,n); S_pre = zeros(m,n);

while true

  temp_L = L_hat - Omega.*(L_hat+S_hat-M)/lips;
  temp_S = S_hat - Omega.*(L_hat+S_hat-M)/lips;

  S = max(temp_S - lambda_2/lips, 0);
  S = S+min(temp_S + lambda_2/lips, 0);

  if choosvd(n, sv) == 1
    [U s V] = lansvd(temp_L, sv, 'L');
    %[U s V] = svd(temp_L, 'econ');
  else
    [U s V] = svd(temp_L, 'econ');
  end
  diagS = diag(s);
  svp = size(find(diagS > lambda_1/lips),1);
  if svp < sv
    flag = 1;
    sv = min(svp + 1, d);
  else
    flag = 0;
    sv = min(svp + round(0.05*d), d);
  end

  newSingularValues = diagS(1:svp) - lambda_1/lips;
  newNuclearNorm = sum(newSingularValues);

  L = U(:, 1:svp) * diag(diagS(1:svp) - lambda_1/lips) * V(:, 1:svp)';

  obj_value = 0.5*norm(Omega.*(L+S-M),'fro')^2 + lambda_1*newNuclearNorm ...
    +lambda_2*norm(vec(S),1);

  fprintf('this is the %d th iteration; sv = %d; obj. value = %d \n',...
    iter, sv, obj_value);

  hist(iter) = obj_value;

  if flag
    if obj_value < obj_target
      break;
    end
  end

  t_pre = t;
  t = (1+4*t^2)^0.5/2 + 0.5;
  alpha = (t_pre-1)/t;

  L_hat = L + alpha*(L-L_pre);
  S_hat = S + alpha*(S-S_pre);

  iter = iter +1;
  L_pre = L;  S_pre = S;

  if(iter > iter_no)
    break;
  end
end

output.L = L;
output.S = S;
output.hist = hist;
end

author: andrewssobral

repo: lrslibrary


text:
function [ cfg ] = default_parameters_dat()
%DEFAULT_PARAMETERS_DAT Default parametrization
  cfg = struct('show_figures', false);

  % Image scaling
  cfg.img_scale_target_diagonal = 75; % Length of object hypothesis diagonal (Used to downscale image)

  % Search and surrounding regions
  cfg.search_win_padding = 2; % Search win = hypothesis + X * max(w,h of hypothesis)
  cfg.surr_win_factor = 1.9;  % Surrounding win = X * surr_win_factor

  % Appearance model
  cfg.color_space = 'rgb';                % 'rgb', 'hsv', or 'lab'
  cfg.num_bins = 16;                      % Number of bins per channel
  cfg.bin_mapping = getBinMapping(cfg.num_bins); % Maps pixel values from [0, 255] to the corresponding bins
  cfg.prob_lut_update_rate = .05;         % Update rate for LUT
  cfg.distractor_aware = true;            % Toggle distractor-awareness
  cfg.adapt_thresh_prob_bins = 0:0.05:1;  % Bins for adaptive threshold.

  % Motion estimation
  cfg.motion_estimation_history_size = 5;  % Motion based on past X frames (Set to 0 to disable motion estimation)

  % NMS-based localization
  cfg.nms_scale = 1;                  % NMS yields rects scaled by X - no scaling typically achieves best localization (Lower scales yield more distractors - so choose wisely)
  cfg.nms_overlap = .9;               % Overlap between candidate rectangles for NMS
  cfg.nms_score_factor = .5;          % Report all rectangles with score >= X * best_score
  cfg.nms_include_center_vote = true; % Prefer hypothesis with highly confident center regions

  % camera correction part
  cfg.transform = '';
end

author: flyers

repo: drone-tracking


text:
% instantiate the library
disp('Loading the library...');
lib = lsl_loadlib();

% resolve a stream...
disp('Resolving an EEG stream...');
result = {};
while isempty(result)
    result = lsl_resolve_byprop(lib,'type','EEG'); end

% create a new inlet
disp('Opening an inlet...');
inlet = lsl_inlet(result{1});

disp('Now receiving chunked data...');
while true
    % get chunk from the inlet
    [chunk,stamps] = inlet.pull_chunk();
    for s=1:length(stamps)
        % and display it
        fprintf('%.2f\t',chunk(:,s));
        fprintf('%.5f\n',stamps(s));
    end
    pause(0.05);
end

author: PatternRecognition

repo: OpenBMI


text:
function final_score = aggregation_attention(score, attention, T)

if size(score, 1)
    score = cat(1, score, score);
    attention = cat(1, attention, attention);
end

if size(score, 3) > 1
score = double(squeeze(mean(score, 2)));
attention = double(squeeze(mean(attention, 2)));
end

attention = attention/T;
softmax_attention = exp(attention);
softmax_attention = bsxfun(@rdivide, softmax_attention, sum(softmax_attention));

final_score = bsxfun(@times, score, softmax_attention);
final_score = sum(final_score, 1);

end

author: wanglimin

repo: UntrimmedNet


text:
function [condMu,condCov] = prtConditionalMvnMuCov(x,indices,globalMu,globalCov)
%[condMu,condCov] = conditionalMuCov(x,indices,globalMu,globalCov)

indices2 = indices;
indices1 = setdiff(1:length(globalMu),indices);

mu2 = globalMu(indices2);
mu1 = globalMu(indices1);

cov22 = globalCov(indices2,indices2);
cov12 = globalCov(indices1,indices2);
cov21 = globalCov(indices2,indices1);
cov11 = globalCov(indices1,indices1);

condMu = mu1 + (cov12*cov22^-1*(x - mu2)')';
condCov = cov11 - cov12*cov22^-1*cov21;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%  LU 分解的紧凑方式
function [L,U,x] = my_lu(A,b)

[m,n]=size(A);
if (m~=n)
    fprintf('\n Error: A 不是方阵!\n'); return;
end

L = eye(n);
U = zeros(n);

for k = 1:n
    U(k,k) = A(k,k)-sum(L(k,1:k-1)'.*U(1:k-1,k));
    if (U(k,k) == 0)
        fprintf('\n Error: 第 %d 个主元＝0 !\n',k); return;
    end
    for j = k+1:n
        U(k,j) = A(k,j) - sum(L(k,1:k-1)'.*U(1:k-1,j));
        L(j,k) = (A(j,k) - sum( L(j,1:k-1)'.*U(1:k-1,k) ))/U(k,k);
    end
end

% 求解Ly=b 和 Ux=y
y(1) = b(1);
for i=2:n
    summary = 0;
    for t=1:i-1
        summary = summary + L(i,t) * y(t);
    end
    y(i)=b(i) - summary;
end
x(n) = y(n)/U(n,n);
for i = n-1 : -1 : 1
    summary = 0;
    for t=i+1:n
        summary = summary + U(i,t) * x(t);
    end
    x(i)=(y(i) - summary)/U(i,i);
end

% END

repo: NumericalAnalysis


text:
function phi = basis_brick20 ( n, p )

%*****************************************************************************80
%
%% BASIS_BRICK20: BRICK20 basis functions at natural coordinates.
%
%  Discussion:
%
%        8----19---7
%       /|        /|
%     20 |      18 |        t   s
%     /  16     /  15       |  /
%    5----17---6   |        | /
%    |   |     |   |        |/
%    |   4--11-|---3        0---------r
%   13  /     14  /
%    | 12      | 10
%    |/        |/
%    1----9----2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real P(3,N), natural coordinates of evaluation points.
%
%    Output, real PHI(20,N), the basis function values.
%
  phi(1,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( - p(1,1:n) - p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(2,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( + p(1,1:n) - p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(3,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( + p(1,1:n) + p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(4,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( - p(1,1:n) + p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(5,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( - p(1,1:n) - p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;
  phi(6,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( + p(1,1:n) - p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;
  phi(7,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( + p(1,1:n) + p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;
  phi(8,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( - p(1,1:n) + p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;

  phi(9,1:n) =  ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(10,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(11,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 + p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(12,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(13,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(14,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(15,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(16,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(17,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;
  phi(18,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;
  phi(19,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 + p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;
  phi(20,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function linplus_test58 ( )

%*****************************************************************************80
%
%% TEST58 tests R8SS_MXV, R8SS_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 9;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST58\n' );
  fprintf ( 1, '  For a symmetric skyline storage matrix,\n' );
  fprintf ( 1, '  R8SS_MXV computes A*x,\n' );
  fprintf ( 1, '  R8SS_PRINT prints it.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ na, diag, a, seed ] = r8ss_random ( n, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nonzero entries stored is %d\n', na );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Diagonal storage indices:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '%6d  %6d\n', i, diag(i) );
  end
%
%  Replace the random entries by marker values.
%
  ij = 0;
  for j = 1 : n

    if ( j == 1 )
      ilo = 1;
    else
      ilo = diag(j-1) - diag(j) + j + 1;
    end

    for i = ilo : j
      ij = ij + 1;
      a(ij) = 10 * i + j;
    end

  end

  r8ss_print ( n, na, diag, a, '  The R8SS matrix:' );
%
%  Copy the matrix into a general matrix.
%
  a2 = r8ss_to_r8ge ( n, na, diag, a );
%
%  Set the vector X.
%
  x = r8vec_indicator ( n );
%
%  Compute the product.
%
  b = r8ss_mxv ( n, na, diag, a, x );
%
%  Compute the product using the general matrix.
%
  b2 = r8ge_mxv ( n, n, a2, x );
%
%  Compare the results.
%
  r8vec2_print_some ( n, b, b2, 10, '  R8SS_MXV verse R8GE_MXV' );

  return
end


text:
classdef CEC2013_F12 < PROBLEM
% <single> <real> <large>
% Shifted Rosenbrock's function

%------------------------------- Reference --------------------------------
% X. Li, K. Tang, M. N. Omidvar, Z. Yang, and K. Qin, Benchmark functions
% for the CEC'2013 special session and competition on large-scale global
% optimization, RMIT University, Australia, 2013.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        Xopt;	% Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2013.mat'),'Data');
            obj.Xopt = Data{12};
            obj.M    = 1;
            obj.D    = 1000;
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = Rosenbrock(PopDec-repmat(obj.Xopt,size(PopDec,1),1));
        end
    end
end

function F = Rosenbrock(X)
    F = sum(100*(X(:,1:end-1).^2-X(:,2:end)).^2+(X(:,1:end-1)-1).^2,2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function xy = padua_points ( l )

%*****************************************************************************80
%
%% PADUA_POINTS returns the Padua points of level L.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marco Caliari, Stefano de Marchi, Marco Vianello,
%    Bivariate interpolation on the square at new nodal sets,
%    Applied Mathematics and Computation,
%    Volume 165, Number 2, 2005, pages 261-274.
%
%  Parameters:
%
%    Input, integer L, the level of the set.
%    0 <= L
%
%    Output, real XY(2,N), the Padua points of level L.
%
  n = ( ( l + 1 ) * ( l + 2 ) ) / 2;

  xy = zeros ( 2, n );

  if ( l == 0 )
    xy(1,1) = 0.0;
    xy(2,1) = 0.0;
    return
  end

  k = 0;

  for i = 0 : l

    j_hi = floor ( l / 2 ) + 1;
    if ( mod ( l, 2 ) == 1 && mod ( i, 2 ) == 1 )
      j_hi = j_hi + 1;
    end;

   for j = 1 : j_hi

     k = k + 1;

      if ( i * 2 == l )
        xy(1,k) = 0.0;
      else
        angle1 = i * pi / l;
        xy(1,k) = cos ( angle1 );
      end

      if ( mod ( i, 2 ) == 0 )
        if ( 2 * ( 2 * j - 1 ) == l + 1 )
          xy(2,k) = 0.0;
        else
          angle2 = ( 2 * j - 1 ) * pi / ( l + 1 );
          xy(2,k) = cos ( angle2 );
        end
      else
        if ( 2 * ( 2 * j - 2 ) == l + 1 )
          xy(2,k) = 0.0;
        else
          angle2 = ( 2 * j - 2 ) * pi / ( l + 1 );
          xy(2,k) = cos ( angle2 );
        end
      end

    end

  end

  return
end


text:
function [beta_gibbs, sigma_gibbs, theta_gibbs, ss_record,indH,beta_theta_gibbs]=TVEmagibbs(data_endo,It,Bu,beta0,omega0,psi0,lambda0,Y,X,n,T,k1,q1,p,regimeperiods,names,TVEH)

% function [beta_gibbs psi_gibbs sigma_gibbs ss_record delta_gibbs]=magibbs(data_endo,data_exo,It,Bu,beta0,omega0,psi0,lambda0,Y,X,Z,n,m,T,k1,k3,q1,q2,q3,p)
% performs the Gibbs algortihm 3.5.1 for a MABVAR model, and returns draws from posterior distribution
% inputs:  - matrix 'data_endo': the matrix storing the endogenous time series data used to estimate the model
%          - matrix 'data_exo': the matrix storing the exogenous time series data used to estimate the model
%          - integer 'It': the total number of iterations run by the Gibbs sampler
%          - integer 'Bu': the number of initial iterations discared as burn-in sample
%          - vector 'beta0': the vector containing the mean of the prior distribution for beta, defined in (3.5.17)
%          - matrix 'omega0': the variance-covariance matrix for the prior distribution of beta, defined in (3.5.17)
%          - vector 'psi0': the vector containing the mean of the prior distribution for psi, defined in (3.5.19)
%          - matrix 'lambda0': the variance-covariance matrix for the prior distribution of psi, defined in (3.5.19)
%          - matrix 'Y': the matrix of endogenous variables, defined in (3.5.10)
%          - matrix 'X': the matrix of endogenous regressors, defined in (3.5.10)
%          - matrix 'Z': the matrix of exogenous regressors, defined in (3.5.10)
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'm': the number of exogenous variables in the model
%          - integer 'T': the sample size, i.e. the number of time periods used to estimate the model
%          - integer 'k1': the number of coefficients related to the endogenous variables for each equation in the model
%          - integer 'k3': the number of coefficients related to the exogenous variables for each equation, in the reformulated model (3.5.5)
%          - integer 'q1': the total number of VAR coefficients related to the endogenous variables
%          - integer 'q2': the total number of VAR coefficients related to the exogenous variables
%          - integer 'q3': the total number of VAR coefficients related to the exogenous variables, in the reformulated model (3.5.5)
%          - integer 'p': the number of lags in the model
% outputs: - matrix 'beta_gibbs': the matrix recording the post-burn draws of beta
%          - matrix 'psi_gibbs': the matrix recording the post-burn draws of psi
%          - matrix 'sigma_gibbs': the matrix recording the post-burn draws of sigma
%          - matrix 'delta_gibbs': the matrix recording the post-burn draws of delta
%          - cell 'ss_record': the cell storing the post-burn steady-state values

% this function implements algorithm...

% preliminary tasks

% take the alternative exogenous regime
if isempty(regimeperiods)
    data_exo=[ones(size(data_endo,1),1)];
else
data_exo=[ones(size(data_endo,1),1) zeros(size(data_endo,1),1)];
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),2)=1;
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),1)=0;
end

indH=zeros(T+p,1);
for iR=1:size(data_exo,2)
    indH(data_exo(:,iR)==1)=iR;
end

% create the cell that will store the records for the steady-state
ss_record=cell(n,1); %change

% invert omega0
invomega0=diag(1./diag(omega0));

% invert lambda0
invlambda0=lambda0\eye(length(lambda0));

% step 2: set initial values
% set initial values for B, beta and sigma; as no OLS estimates are available, simply set the value as zeros for B and beta, and identity for sigma
B=zeros(k1,n);

% define the initial value for the inverse of sigma: beacause sigma is identity, this is also identity
invsigma=eye(n);

% preallocate space for the matrix with the equilibrium values
eq=zeros(T+p,n);

q2=length(psi0);

hbar = bear.parfor_progressbar(It,'Progress of the Gibbs sampler');  %create the progress bar

% start iterations
for ii=1:It

hbar.iterate(1);   % update progress by one iteration

% step 3: at iteration ii, draw psi from N, conditional on beta and sigma
% obtain the lambdabar matrix
Ybar=Y-X*B;
Ypsi=bear.vec(Ybar');

t=1;
Fsimple=TVEH(:,:,indH(t+p),t+p);
for k=2:p+1
    Fsimple=Fsimple-B((k-2)*n+1:(k-1)*n,:)'*TVEH(:,:,indH(t+p-(k-1)),t+p-(k-1));
end

for t=2:T
    Ftemp=TVEH(:,:,indH(t+p),t+p);
    for k=2:p+1
        Ftemp=Ftemp-B((k-2)*n+1:(k-1)*n,:)'*TVEH(:,:,indH(t+p-(k-1)),t+p-(k-1));
    end
    Fsimple=cat(1,Fsimple,Ftemp);
end

invOmega=kron(eye(T),invsigma);

CT=(invlambda0+Fsimple'*invOmega*Fsimple)\eye(length(lambda0));
mT=CT*(Fsimple'*invOmega*Ypsi++invlambda0*psi0);

% draw from N(psibar,lambdabar);
theta=mT+chol(bear.nspd(CT),'lower')*randn(q2,1);

% recover equilibrium values from psi
for it=1:T+p
    eq(it,:)=(squeeze(TVEH(:,:,indH(it),it))*theta)'; % compute the equilibrium values given theta
end

% step 4: now that psi/F has been drawn, it is possible to generate Yhat, Xhat and yhat
temp2=data_endo-eq;
temp3=bear.lagx(temp2,p);
Yhat=temp3(:,1:n);
yhat=Yhat(:);
Xhat=temp3(:,n+1:end);

% step 5: next, at iteration ii, draw sigma from IW, conditional on most recent draw for psi and beta
% obtain first Stilde
Stilde=(Yhat-Xhat*B)'*(Yhat-Xhat*B);
% next draw from IW(Stilde,T)
sigma=bear.iwdraw(Stilde,T);
% invert sigma
C=bear.trns(chol(bear.nspd(sigma),'Lower'));
invC=C\speye(n);
invsigma=invC*invC';

% step 6: finally, at iteration ii, draw beta from a N, conditional on most recent draw for psi and sigma
% first obtain the omegabar matrix
invomegabar=invomega0+kron(invsigma,Xhat'*Xhat);
C=bear.trns(chol(bear.nspd(invomegabar),'Lower'));
invC=C\speye(q1);
omegabar=invC*invC';
% following, obtain betabar
betabar=omegabar*(invomega0*beta0+kron(invsigma,Xhat')*yhat);
% draw from N(betabar,omegabar);
beta=betabar+chol(bear.nspd(omegabar),'lower')*randn(q1,1);
% reshape to obtain B
B=reshape(beta,k1,n);
% update U, using the draw obtained for B

   % record the values if the number of burn-in iterations is exceeded
   if ii>Bu
   % record the value of beta
   beta_gibbs(:,ii-Bu)=beta;
   % record the value of sigma (in vectorized form)
   sigma_gibbs(:,ii-Bu)=sigma(:);

   % compute and record the steady-state
   ss=eq;
   % trim p initial conditions to be consistent with the sample
   ss=ss(p+1:end,:);
      % then record in the corresponding cell
      for jj=1:n
      ss_record{jj,1}(ii-Bu,:)=ss(:,jj)';
      end
   % finally, record theta
   theta_gibbs(:,ii-Bu)=theta;

   beta_theta_gibbs(:,ii-Bu)=[beta;theta];
   % if current iteration is still a burn iteration, do not record the result
   else
   end

% step 7: go for next iteration
end

close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
function value = r8_lngam ( x )

%*****************************************************************************80
%
%% R8_LNGAM evaluates the log of the absolute value of gamma of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the logarithm of the absolute value of
%    the gamma function of X.
%
  persistent dxrel
  persistent sq2pil
  persistent sqpi2l
  persistent xmax

  sq2pil = 0.91893853320467274178032973640562;
  sqpi2l = +0.225791352644727432363097614947441;

  if ( isempty ( xmax ) )
    xmax = r8_mach ( 2 ) / log ( r8_mach ( 2 ) );
    dxrel = sqrt ( r8_mach ( 4 ) );
  end

  y = abs ( x );

  if ( y <= 10.0 )
    value = log ( abs ( r8_gamma ( x ) ) );
    return
  end

  if ( xmax < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_LNGAM - Fatal error!\n' );
    fprintf ( 1, '  Result overflows, |X| too big.\n' );
    error ( 'R8_LNGAM - Fatal error!' )
  end

  if ( 0.0 < x )
    value = sq2pil + ( x - 0.5 ) * log ( x ) - x + r8_lgmc ( y );
    return
  end

  sinpiy = abs ( sin ( pi * y ) )

  if ( sinpiy == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_LNGAM - Fatal error!\n' );
    fprintf ( 1, '  X is a negative integer.\n' );
    error ( 'R8_LNGAM - Fatal error!' )
  end

  value = sqpi2l + ( x - 0.5 ) * log ( y ) - x - log ( sinpiy ) ...
    - r8_lgmc ( y );

  if ( abs ( ( x - r8_aint ( x - 0.5 ) ) * r8_lngam / x ) < dxrel )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_LNGAM - Warning!\n' );
    fprintf ( 1, '  Result is half precision because\n' );
    fprintf ( 1, '  X is too near a negative integer.\n' );
  end

  return
end


text:
function [varargout]=textView(varargin)

%%

defaultOptionStruct.BackgroundColor=[44 48 55]/255;
defaultOptionStruct.ForegroundColor=[168 177 190]/255;
defaultOptionStruct.HorizontalAlignment='Left';
defaultOptionStruct.FontName='Monospaced';
defaultOptionStruct.FontWeight='Bold';
defaultOptionStruct.FontSize=12;

switch nargin
    case 1
        T=varargin{1};
        optionStruct=[];
    case 2
        T=varargin{1};
        optionStruct=varargin{2};
end

fileName='';
if ~iscell(T)
    if exist(T,'file')==2
        fileName=T;
        T=txtfile2cell(T);
    end
end

%Fix option structure, complete and remove empty values
[optionStruct]=structComplete(optionStruct,defaultOptionStruct,1);

%%

t=(1:1:numel(T));
n=numel(sprintf('%d',max(t)));
t=sprintf(['%-',num2str(n),'.0d    \n'],t);
t=t(1:end-1);
t=strsplit(t,'\n')';
T=strcat(t,T);

%Open figure
figStruct.Name=fileName;
figStruct.Color=optionStruct.BackgroundColor;
figStruct.MenuBar='none';
figStruct.vcw=0;
figStruct.efw=0;

hf = cFigure(figStruct);

hPan = uipanel(hf, 'Title',fileName, ...
    'Units','pixels','Position',[10 10 hf.Position(3)-20 hf.Position(4)-20],...
    'BorderType','none','BackgroundColor',defaultOptionStruct.BackgroundColor,...
    'ForegroundColor',defaultOptionStruct.ForegroundColor);
hPan.Units='Normalized';
hEdit = uicontrol(hPan, 'Style','edit', 'FontSize',12, ...
    'Min',0, 'Max',2, 'HorizontalAlignment','left', ...
    'Units','normalized', 'Position',[0 0 1 1], ...
    'String',T);

optionSet=fieldnames(optionStruct);
for q=1:1:numel(optionSet)
   hEdit.(optionSet{q})=optionStruct.(optionSet{q});
end

% enable horizontal scrolling
try
jEdit = findjobj(hEdit);
jEditbox = jEdit.getViewport().getComponent(0);
jEditbox.setWrapping(false);                % turn off word-wrapping
jEditbox.setEditable(false);                % non-editable
set(jEdit,'HorizontalScrollBarPolicy',30);  % HORIZONTAL_SCROLLBAR_AS_NEEDED

% maintain horizontal scrollbar policy which reverts back on component resize
hjEdit = handle(jEdit,'CallbackProperties');
set(hjEdit, 'ComponentResizedCallback',...
    'set(gcbo,''HorizontalScrollBarPolicy'',30)')
catch
end
drawnow;

%%
if nargout>0
    varargout{1}=hf;
end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function script = ma_updateOptimScript(x, script, variables)
    if(any(isnan(x)))
        error(sprintf('Recieved NaN from optimizer algorithm when attempting to update script with new variable values.  Generally this means the algorithm has failed in some way.\n\nPlease try a different optimizer algorithm via the menu:\n\nScript -> Execution Settings -> Optimizer Algorithm.'));
    end

    i = 1;
    for(varC = variables) %#ok<*NO4LP>
        var = varC{1};
        event = script{var.eventNum};

        xSubset = x(i:i+length(var.x0)-1);
        event = ma_updateOptimEvent(xSubset, event);

        script{var.eventNum} = event;

        i = i + length(var.x0);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [FO,ntrials] = getFractionalOccupancy (Gamma,T,options,dim,alignment)
%
% Computes the fractional occupancy
% - across trials if dim==1 (i.e. FO is time by no.states)
%   This is indicating how the fractional occupancy evolves as a function
%   of time across trials.
% - across time if dim==2 (i.e. FO is no.trials/subjects by no.states)
%   This is indicating how much of each state each subject has.
% (default, dim=2)
%
% The parameter 'options' must be the same than the one supplied to the
% hmmmar function for training. The structure hmm.train can also be
% supplied
%
% If dim==1, how is the across-time average defined when trials have not
% the same length? The trials can be aligned in two different ways: such
% that they all start at the same time point, or such that they all finish
% at the same time point. The parameter 'alignment' defines this, and can
% take values 'start' or 'finish'.
%
% If dim==1, also, ntrials will return a (time by 1) vector telling how many
% trials were used to compute the average at that time point. This is
% useful when trials have different lengths.
%
% Note: this can be applied to the data as well, if you want to look at the
% evoked response.
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2016)

if nargin<3, options = struct(); options.Fs = 1; options.downsample = 1; end
if ~isfield(options,'Fs'), options.Fs = 1; end
if ~isfield(options,'downsample'), options.downsample = options.Fs; end
if nargin<4, dim = 2; end
if nargin<5 && dim == 1, alignment = 'start'; end

is_vpath = (size(Gamma,2)==1 && all(rem(Gamma,1)==0));
if iscell(T)
    if size(T,1) == 1, T = T'; end
    for i = 1:length(T)
        if size(T{i},1)==1, T{i} = T{i}'; end
    end
    Nsubj = length(T);
    trials2subjects = zeros(length(cell2mat(T)),1); ii = 1;
    for i = 1:length(T)
        Ntrials = length(T{i});
        trials2subjects(ii:ii+Ntrials-1) = i;
        ii = ii + Ntrials;
    end
    T = cell2mat(T);
elseif dim == 2
    Nsubj = length(T);
    trials2subjects = 1:Nsubj;
end
N = length(T);

r = 1;
if isfield(options,'downsample') && options.downsample>0
    r = (options.downsample/options.Fs);
end

if ~isfield(options,'order') && ~isfield(options,'embeddedlags')
   options.order = (sum(T) - size(Gamma,1)) / length(T);
end

if isfield(options,'tuda') && options.tuda | ...
        isfield(options, 'order') && options.order == 0
    T = ceil(r * T);
elseif isfield(options,'order') && options.order > 0
    T = ceil(r * T);
    T = T - options.order;
elseif isfield(options,'embeddedlags') && length(options.embeddedlags) > 1
    d1 = -min(0,options.embeddedlags(1));
    d2 = max(0,options.embeddedlags(end));
    T = T - (d1+d2);
    T = ceil(r * T);
end

if is_vpath % viterbi path
    vpath = Gamma;
    if options.dropstates==1
        K = length(unique(vpath));
    else
        K = options.K;
    end
    Gamma = zeros(length(vpath),K);
    for k = 1:K
       Gamma(vpath==k,k) = 1;
    end
else
    K = size(Gamma,2);
    %Gamma = Gamma > (2/3);
end

if dim == 2
    FO = zeros(Nsubj,K);
    for j = 1:N
        t0 = sum(T(1:j-1));
        ind = (1:T(j)) + t0;
        jj = trials2subjects(j);
        if length(ind) > 1
            FO(jj,:) = FO(jj,:) + sum(Gamma(ind,:));
        else
            FO(jj,:) = FO(jj,:) + Gamma(ind,:);
        end
    end
    FO = FO ./ repmat(sum(FO,2),1,K);
    ntrials = [];
else
    if all(T==T(1))
        Gamma = reshape(Gamma,[T(1),N,K]);
        FO = squeeze(mean(Gamma,2));
        ntrials = N * ones(T(1),1);
    else
        maxT = max(T);
        FO = zeros(maxT,K);
        ntrials = zeros(maxT,1);
        for j = 1:N
            t0 = sum(T(1:j-1));
            ind_gamma = (1:T(j)) + t0;
            if strcmpi(alignment,'start')
                ind_FO = 1:length(ind_gamma);
            else
                ind_FO = (1:length(ind_gamma)) + (maxT-length(ind_gamma));
            end
            FO(ind_FO,:) = FO(ind_FO,:) + Gamma(ind_gamma,:);
            ntrials(ind_FO) = ntrials(ind_FO) + 1;
        end
        FO = FO ./ repmat(ntrials,1,K);
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function C = exp(X,p)
%EXP          Long exponential function
%
%  C = exp(X,p)
%
%for long number X. Input parameter p is optional; if specified,
%approximate accuracy of C approximately p decimals, otherwise that of X.
%Input X must be less than beta.
%

% written  12/30/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_LONG_BETA = getappdata(0,'INTLAB_LONG_BETA');
  INTLAB_LONG_LOGBETA = getappdata(0,'INTLAB_LONG_LOGBETA');

  % convert decimal precision in beta-digits
  if nargin==2
    longprecision(p);
    p = ceil( p/log10(INTLAB_LONG_BETA) );
  else
    p = size(X.mantissa,2) + 1;
  end

  if any( X.exponent>1 )
    error('exponent too big for long exponential')
  end

  index = ( X.exponent>0 );
  if any( index )
    large = 1;
    E = zeros(size(X.sign));
    while any(index)
      E(index) = E(index) + 1;
      X(index) = X(index)/2;
      index = ( X.exponent>0 );
    end
  else
    large = 0;
  end

  C = 1 + X;
  T = X;
  i = 1;
  while 1
    i = i+1;
    T = T * X / i;
    if all( T.exponent<-p )
      if isequal( longinit('ErrorTerm',0) , 'WithErrorTerm' )
        C.error = errorupdate( 1 , C.error , 0 , 1 , 1 , -p );
      end
      break
    end
    C = C + T;
  end

  if large
    index = ( E~=0 );
    while any(index)
      C(index) = C(index)*C(index);
      E(index) = E(index) - 1;
      index = ( E~=0 );
    end
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function test_convert_event

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY artifact2boolvec artifact2event artifact2trl boolvec2artifact boolvec2event boolvec2trl event2artifact event2boolvec event2trl trl2artifact trl2boolvec trl2event

% The functionality of converting between representations used to be covered by
% fieldtrip/private/convert_event but has been split over multiple functions to
% provide a better overview of the different conversions.

[ftver, ftpath] = ft_version;
cd(fullfile(ftpath, 'private'));

if false
  % this section helps to get all files opened in the editor
  edit	artifact2boolvec.m
  edit	artifact2event.m
  edit	artifact2trl.m
  edit	boolvec2artifact.m
  edit	boolvec2event.m
  edit	boolvec2trl.m
  edit	event2artifact.m
  edit	event2boolvec.m
  edit	event2trl.m
  edit	trl2artifact.m
  edit	trl2boolvec.m
  edit	trl2event.m
end

% there are 4 representations, hence there are 4x3=12 conversions to be tested

%%

artifact = [
  1 9
  11 19
  21 29
  31 39
  41 49
  ];

if true
  artifact = array2table(artifact);
  artifact.Properties.VariableNames = {'begsample', 'endsample'};
  artifact.type = {'a', 'b', 'c', 'd', 'e'}';
  artifact.value = {1, '2', 3, 4, 5}';
end

boolvec   = artifact2boolvec(artifact, 'endsample', 100)
event     = artifact2event(artifact)
trl       = artifact2trl(artifact)

% convert back
artifact1  = boolvec2artifact(boolvec)
artifact2  = event2artifact(event)
artifact3  = trl2artifact(trl)

%%

boolvec = zeros(1,100);
boolvec(1:10)  = 1;
boolvec(21)    = 2;

artifact  = boolvec2artifact(boolvec)
event     = boolvec2event(boolvec)
trl       = boolvec2trl(boolvec)

% convert back
boolvec1 = artifact2boolvec(artifact, 'endsample', 100)
boolvec2 = event2boolvec(event, 'endsample', 100)
boolvec3 = trl2boolvec(trl, 'endsample', 100)

%%

event = [];

event(1).type     = 'trigger';
event(1).value    = 1;
event(1).sample   = 1;
event(1).duration = [];
event(1).offset   = 0;

event(2).type     = 'trigger';
event(2).value    = 2;
event(2).sample   = 10;
event(2).duration = 2;
event(2).offset   = [];

artifact  = event2artifact(event)
boolvec   = event2boolvec(event, 'endsample', 100)
trl       = event2trl(event)

% convert back
event1 = artifact2event(artifact)
event2 = boolvec2event(boolvec)
event3 = trl2event(trl)

%%

trl = [
  1 9 0
  11 19 0
  21 29 0
  31 39 0
  41 49 0
  ];

if false
  trl = array2table(trl);
  trl.Properties.VariableNames = {'begsample', 'endsample', 'offset'};
  trl.type = {'a', 'b', 'c', 'd', 'e'}';
  trl.value = {1, '2', 3, 4, 5}';
end

artifact  = trl2artifact(trl)
boolvec   = trl2boolvec(trl, 'endsample', 100)
event     = trl2event(trl)

% convert back
trl1 = artifact2trl(artifact)
trl2 = boolvec2trl(boolvec)
trl3 = event2trl(event)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Support function for tarjan.m
% "Performs a single depth-first search of the graph, finding all
% successors from the node vi, and reporting all strongly connected
% components of that subgraph."
% See: http://en.wikipedia.org/wiki/Tarjan's_strongly_connected_components_algorithm
%
% INPUTs: start node, vi;
%         graph structure (list), L
%         tarjan.m variables to update: S, ind, v, GSCC
% OUTPUTs: updated tarjan.m variables: S, ind, v, GSCC
%
% Note: Contains recursion.
% Other routines used: strongConnComp.m
% GB: last updated, Sep 22 2012

function [GSCC,S,ind,v]=strongConnComp(vi,S,ind,v,L,GSCC)

v(vi).index = ind;                % Set the depth index for vi
v(vi).lowlink = ind;
ind = ind + 1;

S = [vi S];                         % Push vi on the stack

for ll=1:length(L{vi})
  vj = L{vi}(ll);                   % Consider successors of vi

  if isempty(v(vj).index)           % Was successor vj visited?

    [GSCC,S,ind,v]=strongConnComp(vj,S,ind,v,L,GSCC);   % Recursion
    v(vi).lowlink = min([v(vi).lowlink, v(vj).lowlink]);

  elseif not(isempty(find(S==vj)))            % Is vj on the stack?
    v(vi).lowlink = min([v(vi).lowlink, v(vj).index]);

  end
end

if v(vi).lowlink == v(vi).index     % Is v the root of an SCC?

  SCC = [vi];
  while 1
    vj = S(1); S = S(2:length(S));
    SCC = [SCC vj];

    if vj==vi; SCC = unique(SCC); break; end
  end

  GSCC{length(GSCC)+1} = SCC;

end

author: aeolianine

repo: octave-networks-toolbox


text:
function F = asec(F, varargin)
%ASEC   Inverse secant of a CHEBFUN.
%   ASEC(F) computes the inverse secant of the CHEBFUN F.
%
%   ASEC(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also SEC, ASECD.

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

% Call the compose method:
F = compose(F, @asec, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ t_start, y_start ] = p34_start ( neqn )

%*****************************************************************************80
%
%% P34_START returns the starting point for problem p34.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start = [ 2.0; 2.0 ];

  return
end


text:
function [sigma,u,eqn,info] = elasticitymfemP1P0(node,elem,pde,bdFlag,option)

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

time = cputime;  % record assembling time

d = 2; % two dimensions
N = size(node,1);
NT = size(elem,1);
Nsigma = 3*N;
Nu = 2*NT;
Ndof = Nsigma + Nu;

%% Assemble matrix
[Dphi,area] = gradbasis(node,elem);

%% Mass matrix of linear (P1) element
M = sparse(N,N);
for i = 1:3
    for j = i:3
        ii = double(elem(:,i));
        jj = double(elem(:,j));
        if (j==i)
            M = M + sparse(ii,jj,area/6,N,N);
        else
            M = M + sparse([ii;jj],[jj;ii],[area/12; area/12],N,N);
        end
    end
end

%% Compliance tensor
lambda = pde.lambda;
mu = pde.mu;
C1 = 1/(2*mu);
C2 = lambda/(2*mu + d*lambda);
% E = mu*(3*lambda+2*mu)/(lambda+mu);
% nu = lambda/(2*(lambda+mu));
A = sparse(C1*(eye(3,3) - C2*([1 1 0]'*[1 1 0])));

%% Matrix for (Asigma,tau)
Am = kron(A,M);

%% Div operator
elem2dofsigma = zeros(NT,3,3);
elem2dofsigma(:,:,1) = double(elem);
for k = 2:3
    elem2dofsigma(:,:,k) = elem2dofsigma(:,:,k-1)+N;
end
elemIdx = (1:NT)';
elem2dofu = [elemIdx elemIdx+NT];
Dx = squeeze(Dphi(:,1,:)).*repmat(area,1,3);
Dy = squeeze(Dphi(:,2,:)).*repmat(area,1,3);
clear Dphi
% u1: dx sigma(1) + dy sigma(3)
% u2: dx sigma(3) + dy sigma(2)
B = sparse(repmat(elem2dofu(:,1),1,3),elem2dofsigma(:,:,1),Dx,Nu,Nsigma) ...
  + sparse(repmat(elem2dofu(:,1),1,3),elem2dofsigma(:,:,3),Dy,Nu,Nsigma) ...
  + sparse(repmat(elem2dofu(:,2),1,3),elem2dofsigma(:,:,3),Dx,Nu,Nsigma) ...
  + sparse(repmat(elem2dofu(:,2),1,3),elem2dofsigma(:,:,2),Dy,Nu,Nsigma);

%% Stabilization
T = auxstructure(elem);
edge2elem = T.edge2elem;
elem2edge = T.elem2edge;
[normal,edgeLength,unitNormal] = edgenormal(node,T.edge);
clear T;
harea = edgeLength.^2;
% elementwise part: [u_in_k']:[u_jn_k']=n_k(i)n_k(j);
% Use Dphi as a scaled outwards normal vector to compute - int_F h[u n'][v n']
C = sparse(Nu,Nu);
for i = 1:2
    for j = i:2
        ii = elem2dofu(:,i);
        jj = elem2dofu(:,j);
        Cij = zeros(NT,1);
        for k = 1:3 % sum of all four faces
            fk = elem2edge(:,k);
            Cij = Cij + harea(fk).*((i==j) + unitNormal(fk,i).*unitNormal(fk,j));
        end
        if (j==i)
            C = C + sparse(ii,jj,-Cij,Nu,Nu);
        else
            C = C + sparse([ii;jj],[jj;ii],repmat(-Cij,1,2),Nu,Nu);
        end
    end
end
% cross face part
for i = 1:2
    for j = 1:2
        ii = elem2dofu(edge2elem(:,1),i);
        jj = elem2dofu(edge2elem(:,2),j);
        Cij = harea.*((i==j) + unitNormal(:,i).*unitNormal(:,j));
        C = C + sparse([ii;jj],[jj;ii],repmat(-Cij,1,2),Nu,Nu);
    end
end

%% Assemble the right hand side
F = zeros(Ndof,1);
fu = zeros(NT,2);
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order
end
if ~isempty(pde.f)
	[lambda,weight] = quadpts3(option.fquadorder);
	nQuad = size(lambda,1);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
		fp = pde.f(pxy);
		fu = fu - weight(p)*fp;
    end
    fu = fu.*repmat(area,1,2);
end
clear fp
F((Nsigma+1):Ndof,1) = fu(:);

%% Boundary Conditions
if ~exist('bdFlag','var'), bdFlag = []; end
eqn = struct('Am',Am,'B',B,'C',C,'f',F(1:Nsigma),'g',F(Nsigma+1:end));
assembleTime = cputime - time;
%% Solver
bigA = [Am B'; B C];
bigu = bigA\F;
sigma = bigu(1:Nsigma);
u = bigu(Nsigma+1:end);

%% Output information
info.assembleTime = assembleTime;

repo: ifem


text:
function [post nlZ dnlZ] = infGrid(hyp, mean, cov, lik, x, y, opt)

% Inference for a GP with Gaussian likelihood and covGrid covariance.
% The (Kronecker) covariance matrix used is given by:
%   K = kron( kron(...,K{2}), K{1} ) = K_p x .. x K_2 x K_1.
%
% Compute a parametrization of the posterior, the negative log marginal
% likelihood and its derivatives w.r.t. the hyperparameters.
% The result is exact for complete grids, otherwise results are approximate.
% See also "help infMethods".
%
% The function takes a specified covariance function (see covFunctions.m) and
% likelihood function (see likFunctions.m), and is designed to be used with
% gp.m and in conjunction with covFITC and likGauss.
%
% Copyright (c) by Hannes Nickisch and Andrew Wilson, 2014-12-03.
%
% See also INFMETHODS.M, COVGRID.M.

if iscell(lik), likstr = lik{1}; else likstr = lik; end
if ~ischar(likstr), likstr = func2str(likstr); end
if ~strcmp(likstr,'likGauss')               % NOTE: no explicit call to likGauss
  error('Inference with infGrid only possible with Gaussian likelihood.');
end
cov1 = cov{1}; if isa(cov1, 'function_handle'), cov1 = func2str(cov1); end
if ~strcmp(cov1,'covGrid'); error('Only covGrid supported.'), end    % check cov

xg = cov{3}; p = numel(xg);                                    % underlying grid
N = 1; D = 0; for i=1:p, N = N*size(xg{i},1); D = D+size(xg{i},2); end    % dims
[K,M,xe] = feval(cov{:}, hyp.cov, x);     % evaluate covariance mat constituents
xe = M*xe; n = size(xe,1);
m = feval(mean{:}, hyp.mean, xe);                         % evaluate mean vector

if nargin<=6, opt = []; end                        % make opt variable available
if isfield(opt,'cg_tol'),   cgtol = opt.cg_tol;       % stop conjugate gradients
else cgtol = 1e-5; end
if isfield(opt,'cg_maxit'), cgmit = opt.cg_maxit;      % number of cg iterations
else cgmit = min(n,20); end
for j=1:numel(K)
  if iscell(K{j}) && strcmp(K{j}{1},'toep'), K{j} = toeplitz(K{j}{2}); end
end

sn2 = exp(2*hyp.lik);                               % noise variance of likGauss
V = cell(size(K)); E = cell(size(K));                 % eigenvalue decomposition
for j=1:numel(K)
  if iscell(K{j}) && strcmp(K{j}{1},'toep')
       [V{j},E{j}] = eigr(toeplitz(K{j}{2}));
  else [V{j},E{j}] = eigr(K{j});
  end
end
e = 1; for i=1:numel(E), e = kron(diag(E{i}),e); end       % eigenvalue diagonal
if numel(m)==N
  s = 1./(e+sn2); ord = 1:N;                  % V*diag(s)*V' = inv(K+sn2*eye(N))
else
  [eord,ord] = sort(e,'descend');              % sort long vector of eigenvalues
  s = 1./((n/N)*eord(1:n) + sn2);               % approx using top n eigenvalues
end

if numel(m)==N               % decide for Kronecker magic or conjugate gradients
  L = @(k) kronmvm(V,repmat(-s,1,size(k,2)).*kronmvm(V,k,1));     % mvm callback
else                                                % go for conjugate gradients
  mvm = @(t) mvmK(t,K,sn2,M);                             % single parameter mvm
  L = @(k) -solveMVM(k,mvm,cgtol,cgmit);
end
alpha = -L(y-m);

post.alpha = alpha;                            % return the posterior parameters
post.sW = ones(n,1)/sqrt(sn2);                         % sqrt of noise precision
post.L = L;                           % function to compute inv(K+sn2*eye(N))*Ks

if nargout>1                               % do we want the marginal likelihood?
  lda = -sum(log(s));                                               % exact kron
  % exact: lda = 2*sum(log(diag(chol(M*kronmvm(K,eye(N))*M'+sn2*eye(n)))));
  nlZ = (y-m)'*alpha/2 + n*log(2*pi)/2 + lda/2;
  if nargout>2                                         % do we want derivatives?
    dnlZ = hyp;     % allocate space for derivatives, define Q=inv(M*K*M'+sn2*I)
    Mtal = M'*alpha;                          % blow up alpha vector from n to N
    for i = 1:numel(hyp.cov)
      dK = feval(cov{:}, hyp.cov, x, [], i);
      for j=1:numel(dK)                                        % expand Toeplitz
        if iscell(dK{j}) && strcmp(dK{j}{1},'toep')
          dK{j} = toeplitz(dK{j}{2});
        end
      end
      P = cell(size(dK));               % auxiliary Kronecker matrix P = V'*dK*V
      for j = 1:numel(dK), P{j} = sum((V{j}'*dK{j}).*V{j}',2); end
      p = 1; for j=1:numel(P), p = kron(P{j},p); end               % p = diag(P)
      p = (n/N)*p(ord(1:n));        % approximate if incomplete grid observation
      dnlZ.cov(i) = (p'*s - Mtal'*kronmvm(dK,Mtal))/2;
    end
    dnlZ.lik = sn2*(sum(s) - alpha'*alpha);                  % sum(s) = trace(Q)
    for i = 1:numel(hyp.mean)
      dnlZ.mean(i) = -feval(mean{:}, hyp.mean, xe, i)'*alpha;
    end
  end
end

% Solve x=A*b with symmetric A(n,n), b(n,m), x(n,m) using conjugate gradients.
% The method is along the lines of PCG but suited for matrix inputs b.
function [x,flag,relres,iter,r] = conjgrad(A,b,tol,maxit)
if nargin<3, tol = 1e-10; end
if nargin<4, maxit = min(size(b,1),20); end
x0 = zeros(size(b)); x = x0;
if isnumeric(A), r = b-A*x; else r = b-A(x); end, r2 = sum(r.*r,1); r2new = r2;
nb = sqrt(sum(b.*b,1)); flag = 0; iter = 1;
relres = sqrt(r2)./nb; todo = relres>=tol; if ~any(todo), flag = 1; return, end
on = ones(size(b,1),1); r = r(:,todo); d = r;
for iter = 2:maxit
  if isnumeric(A), z = A*d; else z = A(d); end
  a = r2(todo)./sum(d.*z,1);
  a = on*a;
  x(:,todo) = x(:,todo) + a.*d;
  r = r - a.*z;
  r2new(todo) = sum(r.*r,1);
  relres = sqrt(r2new)./nb; cnv = relres(todo)<tol; todo = relres>=tol;
  d = d(:,~cnv); r = r(:,~cnv);                           % get rid of converged
  if ~any(todo), flag = 1; return, end
  b = r2new./r2;                                               % Fletcher-Reeves
  d = r + (on*b(todo)).*d;
  r2 = r2new;
end

% solve q = mvm(p) via conjugate gradients
function q = solveMVM(p,mvm,varargin)
  [q,flag,relres,iter] = conjgrad(mvm,p,varargin{:});                 % like pcg
  if ~flag,error('Not converged after %d iterations, r=%1.2e\n',iter,relres),end

% mvm so that q = M*K*M'*p + sn2*p using the Kronecker representation
function q = mvmK(p,K,sn2,M)
  q = M*kronmvm(K,M'*p) + sn2*p;

% Perform a matrix vector multiplication b = A*x with a matrix A being a
% Kronecker product given by A = kron( kron(...,As{2}), As{1} ).
function b = kronmvm(As,x,transp)
if nargin>2 && ~isempty(transp) && transp   % transposition by transposing parts
  for i=1:numel(As), As{i} = As{i}'; end
end
m = zeros(numel(As),1); n = zeros(numel(As),1);                  % extract sizes
for i=1:numel(n)
  if iscell(As{i}) && strcmp(As{i}{1},'toep')
    m(i) = size(As{i}{2},1); n(i) = size(As{i}{2},1);
  else [m(i),n(i)] = size(As{i});
  end
end
d = size(x,2);
b = x;
for i=1:numel(n)
  a = reshape(b,[prod(m(1:i-1)), n(i), prod(n(i+1:end))*d]);    % prepare  input
  tmp = reshape(permute(a,[1,3,2]),[],n(i))*As{i}';
  b = permute(reshape(tmp,[size(a,1),size(a,3),m(i)]),[1,3,2]);
end
b = reshape(b,prod(m),d);                        % bring result in correct shape

% Real eigenvalues and eigenvectors up to the rank of a real symmetric matrix.
% Decompose A into V*D*V' with orthonormal matrix V and diagonal matrix D.
% Diagonal entries of D obave the rank r of the matrix A as returned by
% the call rank(A,tol) are zero.
function [V,D] = eigr(A,tol)
[V,D] = eig((A+A')/2); n = size(A,1);    % decomposition of strictly symmetric A
d = max(real(diag(D)),0); [d,ord] = sort(d,'descend');        % tidy up and sort
if nargin<2, tol = size(A,1)*eps(max(d)); end, r = sum(d>tol);     % get rank(A)
d(r+1:n) = 0; D = diag(d);                % set junk eigenvalues to strict zeros
V(:,1:r) = real(V(:,ord(1:r))); V(:,r+1:n) = null(V(:,1:r)'); % ortho completion

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function pf = calcPoleFigure(pf,odf,varargin)
% simulate pole figure
%
% Syntax
%   pf = calcPoleFigure(pf,odf)
%
% Input
%  pf  - meassured @PoleFigure
%  odf - @SO3Fun
%
% Output
%  pf - recomputed @PoleFigure
%
% See also
% SO3Fun/calcPoleFigure

pf = calcPoleFigure(odf,pf.allH,pf.allR,'superposition',pf.c,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function handleFileName(obj)

%--IF EMPTY FILENAME, CHECK FOR PLOT TITLES--%
if isempty(obj.PlotOptions.FileName)
    for t = 1:obj.State.Figure.NumTexts
        if obj.State.Text(t).Title
            str = get(obj.State.Text(t).Handle,'String');
            interp = get(obj.State.Text(t).Handle,'Interpreter');
            obj.PlotOptions.FileName = parseString(str,interp);
        end
    end
end

%--IF FILENAME IS STILL EMPTY SET TO UNTITLED--%
if isempty(obj.PlotOptions.FileName)
    obj.PlotOptions.FileName = 'untitled';
end

end

author: goGPS-Project


text:
function g = std(f, varargin)
%STD   Standard deviation of a CHEBFUN3 along one variable.
%   G = STD(F) returns the standard deviation of F in the x-variable
%   (default). If F is defined on the cuboid [a, b] x [c, d] x [e, g] then
%
%                      b
%                     /
%     G.^2 = 1/(b-a)  | (F(x,y,z) - mean(F, 1))^2 dx
%                     /
%                     a
%
%   The output G is a CHEBFUN2 object over the rectangle [c, d] x [e, g].
%
%   G = STD(F, FLAG, DIM) takes the standard deviation along the x, y, or
%   z directions if DIM = 1, 2, or 3, respectively.  FLAG is ignored and
%   kept in this function so the syntax agrees with the Matlab STD command.
%
% See also CHEBFUN/STD, CHEBFUN2/STD, CHEBFUN2/STD2, CHEBFUN3/STD2 and
% CHEBFUN3/STD3.

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

% Empty check:
if ( isempty(f) )
    g = chebfun2();
    return
end

dom = f.domain;
if ( nargin < 3 )
    dim = 1;   % default to std over the x-variable.
elseif ( nargin == 3 )
    dim = varargin{ 2 };
else
    error( 'CHEBFUN:CHEBFUN3:std:nargin', 'Too many input arguments.' );
end

%   std(X) = sqrt( mean( (X - mean(X) )^2 ) ).
if ( dim == 1 )          % x-variable.
    mx = chebfun3(@(x,y,z) feval(mean(f, 1), y, z), dom);
    g = sqrt(1/(diff(dom(1:2))) * sum((f - mx).^2, 1)) ;

elseif ( dim == 2 )      % y-variable.
    my = chebfun3(@(x,y,z) feval(mean(f, 2), x, z), dom);
    g = sqrt(1/(diff(dom(3:4))) * sum((f - my).^2, 2));

elseif ( dim == 3 )      % z-variable.
    mz = chebfun3(@(x,y,z) feval(mean(f, 3), x, y), dom);
    g = sqrt(1/(diff(dom(5:6))) * sum((f - mz).^2, 3));

else
    error('CHEBFUN:CHEBFUN3:std:dim', ...
        'Third argument should have value 1, 2 or 3.');
end

end


text:
function [total, uxref, yref] = flowvolume(qmethod,xy,xst,xref,fig)
%flowvolume   plots/explores flow solution on vertical cross-section
%   [total, uxref, yref] = flowvolume(qmethod,xy,xst,xref,fig);
%   input
%          qmethod    mixed method
%          xy         velocity nodal coordinate vector
%          xst        flow solution vector
%          xref       x-location of grid line
%          fig        figure number
%
%   IFISS function: DJS; 7 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nvtx=length(xy(:,1));
ux=xst(1:nvtx); uy=xst(nvtx+1:2*nvtx);
kk=find(xy(:,1)==xref);
if isempty(kk), error('location xref is not a grid-line!'), end
uxref=ux(kk)'; uyref=uy(kk)'; yref=xy(kk,2);
figure(fig)
plot(yref,uxref,'-k'), axis('square'), title('x-section of flow')
%%
%% compute volume of flow using appropriate quadrature
nny=length(yref); hy=yref(2)-yref(1);
if qmethod >1,
%% Simpson's rule
ww=ones(nny,1); ww(2:2:nny-1)=4; ww(3:2:nny-1)=2;
total=(hy/3)*sum(ww.*uxref');
else
%% Trapezium rule
ww=ones(nny,1); ww(2:2:nny-1)=2; ww(3:2:nny-1)=2;
total=(hy/2)*sum(ww.*uxref');
end
fprintf('\nvolume of flow is %10.6e \n',total)
return


text:
function call_price=american_call_bjerkesun_stensland(S, K, r, b, sigma, T)

%--------------------------------------------------------------------------
%
% DESCRIPTION:
%
% Approximation of American call due to Bjerksund and Stensland (1993)
%
%
% Reference:
%
% Petter Bjerksund and Gunnar Stensland,
% "Closed form approximations of american options",
% Scandinavian Journal of Management, 20(5):761-764, 1993.
%
%--------------------------------------------------------------------------
%
% INPUTS:
%
%  S:       spot price
%  K:       exercice price
%  r:       interest rate
%  b:       dividend yield
%  sigma:   volatility
%  T:       time to maturity
%
% This function uses phi1.m
%
%--------------------------------------------------------------------------
%
% OUTPUT:
%
% call_price: price of a call option
%
%--------------------------------------------------------------------------
%
% Author:  Paolo Z., February 2012
%
%--------------------------------------------------------------------------

sigma_sqr=sigma^2;
B0=max(K,(r/(r-b)*K));
beta = (0.5 - b/sigma_sqr) + sqrt( ((b/sigma_sqr-0.5)^2) + 2.0 * r/sigma_sqr);
Binf = beta/(beta-1.0)*K;
hT= - (b*T + 2.0*sigma*sqrt(T))*((K*K)/(Binf-B0));
XT = B0+(Binf-B0)*(1.0-exp(hT));
alpha = (XT-K)*(XT^-beta);
C=alpha*(S^beta) ...
    -alpha*phi1(S,T,beta,XT,XT,r,b,sigma)...
    +phi1(S,T,1,XT,XT,r,b,sigma)...
    -phi1(S,T,1,K,XT,r,b,sigma)...
    -K*phi1(S,T,0,XT,XT,r,b,sigma)...
    +K*phi1(S,T,0,K,XT,r,b,sigma);

c=european_call_contpay(S,K,r,b,sigma,T);

call_price = max(c,C);

author: Sable

repo: mcbench-benchmarks


text:
function [A] = testIJmatlabhypre (inputmatlab, matrix_or_vector, filename, varargin)
%TESTIJMATLABHYPRE Test IJ Martix Conversion between MATLAB and HYPRE Formats.
%
%    This program was tested on Matlab 7.1.
%
%   Description:
%     [A] = testIJmatlabhypre (inputmatlab, matrix_or_vector, filename,
%     varargin) calculates the accuracy  of MATLAB to HYPRE and HYPRE to
%     MATLAB conversion of input matrix or vector.
%
%     To calculate the accuracy of conversion [A], this function converts
%     matrix (or vector) bidirectionally between MATLAB and HYPRE formats
%     and calculates the difference between the input matrix and the result
%     of bidirectional conversion of the input matrix.
%
%     The input format can be either MATLAB or HYPRE. If the input format
%     is MATLAB then the intermediate format is HYPRE. If the input format
%     is HYPRE then the intermediate format is MATLAB.
%
%     The function performs two conversions: forward conversion from input
%     format to intermediate format and backward conversion from
%     intermediate format to input format. The result of such bidirectional
%     conversion is a matrix in the same format as the input matrix.
%
%     [A] = testIJmatlabhypre (inputmatlab, matrix_or_vector, filename,
%     numprocs) calculates the accuracy of MATLAB-HYPRE-MATLAB conversion
%     performed on number of processors given by numprocs.
%
%     [A] = testIJmatlabhypre (inputmatlab, matrix_or_vector, filename, ,
%     imax) calculates the accuracy of HYPRE-MATLAB-HYPRE matrix conversion,
%     where imax is the size of a sparse square HYPRE matrix.
%
%   Inputs:
%     Inputmatlab = 'matlab' if input file is in MATLAB format
%                    'hypre' if input file is in HYPRE format
%     matrix_or_vector = 'matrix' if input file is a matrix
%                      = 'vector' if input file is a vector
%     filename = full file name of the input file
%
%   Example:
%     A = testIJmatlabhypre( 'matlab', 'matrix', 'mymatrix' )
%     will produce the difference [A] between input MATLAB matrix
%     and the result of MATLAB-HYPRE-MATLAB conversion of the input
%     matrix.
%
%     A = testIJmatlabhypre( 'hypre', 'vector', 'mymatrix' )
%     will produce the difference [A] between input HYPRE vector
%     and the result of HYPRE-MATLAB-HYPRE conversion of the input
%     vector.
%
%     A = testIJmatlabhypre( 'matlab', 'matrix', 'mymatrix', 2 )
%     will produce the difference [A] between input MATLAB matrix
%     and the result of MATLAB-HYPRE-MATLAB conversion of the input
%     matrix ran on two processors.
%
%     A = testIJmatlabhypre( 'hypre', 'matrix', 'mymatrix' )
%     will produce the difference [A] between input sparse square HYPRE
%     matrix of size given by imax and the result of HYPRE-MATLAB-HYPRE
%     conversion of the input matrix.
%
%   See also matlabIJ2hypre, matlab2hypreParVectors.m,
%      hypreIJ2matlab.m, hypreParVectors2matlab.m
%
%   Author: Diana Zakaryan, Dept of Mathematics, University of Colorado,
%      Denver, 15-Mar-2005.
%

if nargin > 3
    numprocs = varargin{1}
end

if nargin > 4
    imax = varargin{2}
end

if nargin > 5
    matrix_filename = varargin{3};
    vector_filename = varargin{4};
    A=hypreIJ2matlab(matrix_filename,imax);
    B=hypreParVectors2matlab(vector_filename,1);
    if size(A,2)~=size(B,1)
        error('Dimentions must agree')
    else
    size(A*B)
    end
    clear ('A', 'B');
end

% find out if the check is for matlab format to hypre format
if size(inputmatlab,2) == 6
    % find out if the check is for matrix
    if size(matrix_or_vector,2) == 6
        % code for matlab to hypre in matrix
        % todo add code to check if .mat file or ascii file
        load (filename);
        % converts a matlab formatted matrix to
        % hypre formatted one

        matlab2hypreIJ (matlab_smatrix, numprocs, 'temp_m_hypre');
        imax = size(matlab_smatrix, 1);

        % converts a hypre formatted matrix to
        % matlab formatted one
        matlab_out_m = hypreIJ2matlab ('temp_m_hypre', imax);

        % compare input and output files;
        A1=matlab_smatrix - matlab_out_m
        fprintf('the min of the difference is: %\n', num2str(min(A1)));
        fprintf('the max of the difference is: %\n', num2str(max(A1)));
        clear('matlab_out_m');
        delete ('../temp_m_hypre.*');

    % find out if the check is for vector
    elseif size(matrix_or_vector,2) == 7
        %code for matlab to hypre in vector
        load (filename);

        % converts a matlab formatted vector to
        % hypre formatted one
        matlab2hypreParVectors (matlab_vector,numprocs,'temp_v_hypre');
        % converts a hypre formatted vector to
        % matlab formatted one
        matlab_out_v = hypreParVectors2matlab('temp_v_hypre',1);

        % compare input and output files;
        A2=matlab_vector - matlab_out_v;
        fprintf('the min of the difference is: %\n', num2str(min(A2)));
        fprintf('the max of the difference is: %\n', num2str(max(A2)));
        clear('matlab_out_v');
        delete ('../temp_v_hypre.*');
    else
        %return error message
    end

% find out if the check is for hypre format to matlab format
elseif size(inputmatlab,2) == 5
    % find out if the check is for matrix
    if size(matrix_or_vector,2) == 6
        % code for hypre to matlab in matrix
        % converts a hypre formatted matrix to
        % matlab formatted one
        A=hypreIJ2matlab('matrix',imax);
        % converts a matlab formatted matrix to
        % hypre formatted one
        matlab2hypreIJ(A,numprocs,'matrix_test') %note that for this particular comparison
        % numprocs should be initialized as 2.
        for i=0:numprocs
            s=int2str(i);
            hypre_data_m1=dlmread(strcat('matrix','.0000',s));
            hypre_data_m2=dlmread(strcat('matrix_test','.0000',s));
            % compare input and output files;
            B1=hypre_data_m1-hypre_data_m2
        end

    % find out if the check is for vector
    elseif size(matrix_or_vector,2) == 7
        % converts a hypre formatted vector to
        % matlab formatted one
        A=hypreParVectors2matlab('vectors',1);
        % converts a matlab formatted vector to
        % hypre formatted one
        matlab2hypreParVectors(A,numprocs,'vector_test');
        for i=0:(numprocs-1)
            s=int2str(i);
            hypre_data_v1 = dlmread( strcat('vectors','.0.',s));
            hypre_data_v2 = dlmread( strcat('vector_test','.0.',s));
            % compare input and output files;
            B2=hypre_data_v1-hypre_data_v2
        end
        for i=0:(numprocs-1)
            s=int2str(i);
            hypre_data_v1 = dlmread( strcat('vectors','.0','.INFO.',s));
            hypre_data_v2 = dlmread( strcat('vector_test','.0','.INFO.',s));
            % compare 'INFO' input and output files;
            B3=hypre_data_v1-hypre_data_v2
        end
    else
        %return error message
    end
end


text:
function [B,twom] = modularitydir_f(A,gamma)
% MODULARITYDIR_F returns monolayer Leicht-Newman modularity matrix for directed network given by adjacency matrix A, function handle version
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
%
%   Input: A:  NxN adjacency matrices of a directed network
%          gamma: resolution parameter
%
%   Output: B: function handle where B(i) returns the ith column of
%          [N]x[N] modularity matrix of the monolayer network
%           with adjacency matrix A
%           twom: normalisation constant
%
%   Example of usage: [B,twom]=modularitydir_f(A,gamma);
%          [S,Q]= genlouvain(B);
%          Q=Q/twom;
%   Notes:
%     The matrix A is assumed to be square. This assumption is not checked
%     here.
%
%     For smaller systems, it is potentially more efficient (and easier) to
%     directly use the sparse quality/modularity matrix B in MODULARITY. For
%     large systems with undirected networks, use MODULARITY_F.
%
%     This code serves as a template and can be modified for situations
%     with other wrinkles (e.g., different null models).
%
%     By using this code, the user implicitly acknowledges that the authors
%     accept no liability associated with that use.  (What are you doing
%     with it anyway that might cause there to be a potential liability?!?)
%
%   References
%     Elizabeth A. Leicht and Mark E. J. Newman. "Community structure in
%     Directed Networks", Physical Review Letters 100, 118703 (2008).

if nargin<2||isempty(gamma)
	gamma=1;
end

k=sum(A,2);
d=sum(A,1);
twom=sum(k);
A=(A+A')/2;

B=@(i) full(A(:,i)-gamma/2*(k*d(i)+d'*k(i))/twom);

end

author: GenLouvain

repo: GenLouvain


text:
function out = cryptosubdecode(in,key)

numbers = double(lower(in))-96;
decodednumbers = zeros(1,length(numbers));
for i = 1:length(numbers)
    if numbers(i) >= 1 & numbers(i) <= 26
        decodednumbers(i) = find(key==numbers(i));
    else
        decodednumbers(i) = numbers(i);
    end
end

out = char(decodednumbers+96);


text:
function x = zero_null_right ( m, n )

%*****************************************************************************80
%
%% ZERO_NULL_RIGHT returns a right null vector of the ZERO matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of A.
%
%    Output, real X(N,1), the null vector.
%
  x = ones ( n, 1 );

  return
end


text:
function t = isprtchr(c)
%ISPRTCHR True for printable characters.
%
%   For a string C, ISPRTCHR(C) is 1 for printable characters and 0
%   otherwise.
%
%   See also: ISALNUM, ISALPHA, ISASCII, ISDIGIT, ISLOWER, ISUPPER,
%   ISXDIGIT.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:50:57 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 1, nargin));

   t = ischar(c) & (   ( (  32 <= c ) & ( c <= 126 ) ) ...
                     | ( ( 159 <= c ) & ( c <= 255 ) ) );

   % whether this is correct depends on the coding system used
%   if strncmp(computer, 'PC', 2)
%      t = t | ( ( 130 <= c ) & ( c <= 140 ) ) ...
%            | ( ( 145 <= c ) & ( c <= 156 ) );
%   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function self = gui_test_result(progress_bar, text_runs, error_listbox, max_number_of_tests)
%gui_test_result contructor.
%  The constructor creates an object of the class gui_test_result with the
%  following parameters:
%    progress_bar: handle of the progess bar.
%    text_runs: handle of the text area, which shows the number of test ran,
%               errors and failures.
%    error_listbox: handle of the listbox, which contains the list of
%                   all errors and failures.
%    max_number_of_tests: the number of tests.
%
%  Class Info
%  ==========
%  The class gui_test_result is inherited from test_result and shows
%  the test results and the progress at the mlUnit gui. Normally an
%  instance of the class is only created by gui_test_runner, but not by the
%  user.
%
%  Example
%  =======
%  A gui_test_result is created in gui_test_runner/gui:
%         result = gui_test_result(handles.gui_progress_bar, ...
%             handles.gui_text_runs, ...
%             handles.gui_error_list, ...
%             0);
%
%  See also TEST_RESULT, GUI_TEST_RESULT, GUI_TEST_RUNNER/GUI.

self.progress_bar = progress_bar;
self.text_runs = text_runs;
self.error_listbox = error_listbox;
self.max_number_of_tests = max_number_of_tests;
result = test_result;
self = class(self, 'gui_test_result', result);
reset_progress_bar(self);
reset_texts(self);

repo: REST


text:
classdef nnspnorm < nntest
  methods (Test)
   function basic(test)
      h = 13 ;
      w = 17 ;
      d = 4 ;
      n = 5 ;
      param = [3, 3, 0.1, 0.75] ;
      x = test.randn(h,w,d,n) ;
      y = vl_nnspnorm(x, param) ;
      dzdy = test.rand(h, w, d, n) ;
      dzdx = vl_nnspnorm(x, param, dzdy) ;
      test.der(@(x) vl_nnspnorm(x,param), x, dzdy, dzdx, test.range * 1e-3) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function W = spm_vb_adjacency(edges,weights,N)
% (Weighted) adjacency (or weight) matrix of a graph
% FORMAT W = spm_vb_adjacency(edges,weights,N)
%
% edges    [Nedges x 2] list of neighboring voxel indices
% weights  [Nedges x 1] list of edge weights (unity of not specified)
% N        number of nodes (cardinality of node set)
%
% W        [N x N] matrix of (weighted) edges
% Wij      edge weight between nodes i and j if they are neighbors, otherwise 0
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% Lee Harrison
% $Id: spm_vb_adjacency.m 6079 2014-06-30 18:25:37Z spm $

% Number of edges
Ne = size(edges,1);

% Uniform weights if not specified
if nargin < 2
    weights = ones(Ne,1);
end

% Number of nodes (if N is not specified)
if nargin < 3
    N = max(edges(:));
end

% (Weighted) adjacency matrix
W  = sparse([edges(:,1);edges(:,2)],[edges(:,2);edges(:,1)],...
            [weights;weights],N,N);

author: spm


text:
function g = sphf2cartf(f, lam, th, coord)
%SPHF2CARTF    Wrapper for evaluating a function defined in Cartesian
%
%  G = SPHF2CARTF(F, LAM, TH) evaluates the function handle F = F(x,y,z)
%  at x = cos(lam).*sin(th), y = sin(lam).*sin(th), z = cos(th). This is
%  the co-latitude spherical coordinate system.
%
%  G = SPHF2CARTF(F, LAM, TH, 0) same as SPHF2CARTF(F, LAM, TH).
%
%  G = SPHF2CARTF(F, LAM, TH, 1) evaluates F = F(x,y,z) at
%  at x = cos(lam).*cos(th), y = sin(lam).*cos(th), z = sin(th). This is
%  the latitude spherical coordinate system.

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

% For Developers, recall that:
%   coord - Type of spherical coordinate system:
%          coord = 0 (co-latitude) --> -pi <= lam < pi, 0 <= th <= pi
%          coord = 1 (latitude)    --> -pi <= lam < pi, -pi/2 <= th <= pi/2
% (Default above is co-latitude.)

% TODO: Create separate wrapper functions for latitude/co-latitude so that
% this if can be removed (thus improving performance.

if ( nargin == 3 )
    coord = 0;  % Default is to use co-latitude.
end

if ( coord == 0 )
    % Latitude: 0 <= theta < pi
    x = cos(lam).*sin(th);
    y = sin(lam).*sin(th);
    z = cos(th);
elseif ( coord == 1 )
    % Latitude: -pi/2 <= theta < pi/2
    x = cos(lam).*cos(th);
    y = sin(lam).*cos(th);
    z = sin(th);
else
    error('SPHEREFUN:sphf2cartf:CoordSysUnknown', ['Unknown coordinate '...
        'system for the sphere.']);
end

g = feval(f, x, y, z);

end


text:
%fsqfindpn(lambda,mu,c,m,n)
%   This function finds the probability there are n machines down
%   (in the system) for a machine repair problem (finite source queue).

function out = fsqfindpn(lambda,mu,c,m,n)

p0 = fsqfindp0(lambda,mu,c,m);

if n == 0
    pn = p0;
end

if n >= 1
    if n < c
        pn = nchoosek(m,n)*(lambda/mu)^n*p0;
    end
end

if n >= c
    if n <= m
        pn = nchoosek(m,n)*(lambda/mu)^n*p0*factorial(n)/(factorial(c)*c^(n-c));
    end
end

if n > m
    pn = 0;
end

out = pn;


text:
function varargout = geotiffread_lcog(filename, varargin)
%GEOTIFFREAD Read GeoTIFF file
%
%       GEOTIFFREAD is not recommended. Use READGEORASTER instead.
%
%   [A, R] = GEOTIFFREAD(FILENAME) reads a georeferenced grayscale, RGB, or
%   multispectral image or data grid from the GeoTIFF file specified by
%   FILENAME into A and constructs a spatial referencing object, R.
%
%   [X, CMAP, R] = GEOTIFFREAD(FILENAME) reads an indexed image into X and
%   the associated colormap into CMAP, and constructs a spatial referencing
%   object, R. Colormap values in the image file are rescaled into the
%   range [0,1].
%
%   FILENAME is a string scalar or character vector that specifies the name
%   of the GeoTIFF file. FILENAME can include the folder name. Otherwise,
%   the file must be in the current folder or in a folder on the MATLAB
%   path. If the named file includes the extension '.TIF' or '.TIFF'
%   (either upper or lower case), you can omit the extension from FILENAME.
%
%   A is a two-dimensional array if the file contains a grayscale image or
%   data grid. A is an M-by-N-by-P array if the file contains a color
%   image, multispectral image, hyperspectral image, or data grid. The
%   class of A depends on the storage class of the pixel data in the file
%   which is related to the BitsPerSample property as returned by the
%   IMFINFO function.
%
%   R is a geographic raster reference object if the image or data grid is
%   referenced to a geographic coordinate system, or a map raster reference
%   object if it is referenced to a projected coordinate system.
%
%   [A, REFMAT, BBOX] = GEOTIFFREAD(FILENAME) reads a georeferenced
%   grayscale, RGB, or multispectral image or data grid into A; the
%   corresponding referencing matrix into REFMAT; and the bounding box into
%   BBOX.
%
%   [X, CMAP, REFMAT, BBOX] = GEOTIFFREAD(FILENAME) reads an indexed image
%   into X, the associated colormap into CMAP, the referencing matrix into
%   REFMAT, and the bounding box into BBOX.  The referencing  matrix must
%   be unambiguously defined by the GeoTIFF file, otherwise  it and the
%   bounding box are returned empty.
%
%   [___] = GEOTIFFREAD(URL) reads the GeoTIFF image from a URL. The
%   URL must include the protocol type (e.g., "https://").
%
%   [___] = GEOTIFFREAD(___, IDX) reads one image from a multi-image
%   GeoTIFF file (or URL). IDX is an integer value that specifies the order
%   that the image appears in the file. For example, if IDX is 3,
%   GEOTIFFREAD reads the third image in the file. If you omit this
%   argument, GEOTIFFREAD reads the first image.
%
%   Note
%   ----
%   GEOTIFFREAD imports pixel data using the TIFF-reading capabilities of
%   the MATLAB function IMREAD and likewise shares any limitations of
%   IMREAD.  Consult the IMREAD documentation for specific information on
%   TIFF image support.
%
%   Example
%   -------
%   % Read and display the Boston GeoTIFF image.
%   % Includes material (c) GeoEye, all rights reserved.
%   [boston, R] = geotiffread('boston.tif');
%   figure
%   mapshow(boston, R)
%   axis image off
%
%   See also GEOSHOW, GEOTIFFINFO, GEOTIFFWRITE, IMREAD, MAPSHOW, READGEORASTER

% Copyright 1996-2020 The MathWorks, Inc.

% Verify the input and output argument count.
narginchk(1,2);
nargoutchk(0,4);

% Parse the inputs.
filename = convertStringsToChars(filename);
[filename, url, idx] = parseInputs(filename, varargin);

% Read the info fields from the filename.
info = geotiffinfo_lcog(filename);

% Read the image from the filename.
[A, cmap] = imread(filename,idx);

% If the RefMatrix is empty, try to obtain the spatial information from a
% corresponding worldfile.
if isempty(info.RefMatrix)
    worldfilename = getworldfilename(filename);
    info = getSpatialInfoFromWorldfile(worldfilename, info, size(A));
end

% Delete temporary file from Internet download.
if (url)
    deleteDownload(filename);
end

% Assign output arguments.
varargout = assignOutputArguments(A, cmap, info, nargout);

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

function [filename, url, idx] = parseInputs(filename, inputs)
% Parse the inputs from the cell array, INPUTS.

% Verify the filename and obtain the full pathname.
extensions = {'tif', 'tiff'};
[filename, url] = internal.map.checkfilename(filename, extensions, mfilename, 1, true);

% Check and set the image index number.
if ~isempty(inputs)
  idx = inputs{1};
  attributes = {'real' 'scalar' 'positive'};
  validateattributes(idx, {'numeric'}, attributes, mfilename, 'IDX', 2);
else
  idx = 1;
end

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

function info = getSpatialInfoFromWorldfile(worldFileName, info, rasterSize)
% Obtain the referencing matrix from a world file, if it exits. If so,
% compute the bounding box and construct a spatial referencing object from
% the referencing matrix and update the fields of the INFO structure.
% WORLDFILENAME is a string denoting the name of the world file.

if exist(worldFileName,'file')
    % Check that worldFileName is a file and that it can be opened.
    worldFileName = internal.map.checkfilename(worldFileName, {}, mfilename, 1, false);

    % Obtain the referencing matrix from the world file.
    W = load(worldFileName);
    map.internal.assert(isnumeric(W) && numel(W) == 6, 'map:fileio:expectedSixNumbers');
    W = reshape(W,[2 3]);
    refmat = map.internal.referencingMatrix(W);

    % Calculate the spatial referencing object and bounding box from the
    % referencing matrix if it is not empty.
    if ~isempty(refmat)
        if strcmp(info.ModelType, 'ModelTypeGeographic')
            R = refmatToGeoRasterReference(refmat, rasterSize);
        elseif strcmp(info.ModelType, 'ModelTypeProjected')
            R = refmatToMapRasterReference(refmat, rasterSize);
        else
            R = [];
        end
        h = rasterSize(1);
        w = rasterSize(2);
        outline = [(0.5 + [0 0; 0 w; h w; h 0]), ones(4,1)] * refmat;
        info.BoundingBox = [min(outline); max(outline)];
        info.RefMatrix = refmat;
        info.SpatialRef = R;
    end
end

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

function outputs = assignOutputArguments(A, cmap, info, numOutputs)
% Assign the output arguments based on the number of arguments requested.

outputs{1} = A;
switch numOutputs
    case 2
        if strcmp(info.ColorType, 'indexed')
            % [X, CMAP] = GEOTIFFREAD(...)
            outputs{2} = cmap;
        else
            % [A, R] = GEOTIFFREAD(...)
            outputs{2} = info.SpatialRef;
        end

    case 3
        if strcmp(info.ColorType, 'indexed')
            % [X, CMAP, R] = GEOTIFFREAD(...)
            outputs{2} = cmap;
            outputs{3} = info.SpatialRef;
        else
            % [A, REFMAT, BBOX] = GEOTIFFREAD(...)
            outputs{2} = info.RefMatrix;
            outputs{3} = info.BoundingBox;
        end

    case 4
        % [X, CMAP, REFMAT, BBOX] = GEOTIFFREAD(...)
        outputs{2} = cmap;
        outputs{3} = info.RefMatrix;
        outputs{4} = info.BoundingBox;
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
%FIG_DESCI_CASSI Demonstrate decompress snapshot compressive imaging
%(DeSCI) results of simulated coded aperture snapshot spectral imaging
%(CASSI) `toy` dataset.
% Note: Run this demonstration script after completing the reconstruction
%       process TEST_DESCI_CASSI and obtaining the saved .mat file.
%   See also TEST_DESCI_CASSI, GAPDENOISE_CACTI, GAPDENOISE.
clear; clc;
% close all
addpath('../'); % repository root
addpath(genpath('../packages/')); % packages
% [1] run the corresponding test file for recovery or load the saved
% results
% test_desci_cassi % run the corresponding test file for recovery
load('../results/savedmat/desci_cassi_toy.mat'); % load the saved results

methods   = {'gaptv','desci'}; % all methods for comparison
methnames = {'GAP-TV','DeSCI'}; % corresponding names of the methods
nim = nframe*nmask; % number of images

%% [2.1] demonstrate the reconstructed spectum
colors = {'Brown', 'Orange', 'Blue', 'Red'}; % colors of the four birds
crops  =  [150 280 24 24;  % brown [x y width height] as insertShape
           340 210 24 24;  % orange
           349 396 24 24;  % blue
           349 439 24 24]; % red
locos = {'northwest','northwest','northwest','northwest'};

markers   = {'bx--','r*-'}; % markers

linewidth = 2;
markersize = 4;
textfontsize = 14;
legendfontsize = 12;

% Change default axes fonts.
set(0,'DefaultAxesFontName','Arial');
set(0,'DefaultAxesFontSize',textfontsize);
set(0,'DefaultAxesFontWeight','normal');
% Change default text fonts.
set(0,'DefaultAxesFontName','Arial');
set(0,'DefaultTextFontSize',textfontsize);
set(0,'DefaultTextFontWeight','normal');

f = figure('position',[25 50 700 900]); % PSNR plot
h = tight_subplot(3,2,[0.07 0.08],[.07 .035],[.1 .02]);
% subplot(3,2,1); % show crops on original RGB image
axes(h(1));
imshow(insertShape(orig_rgb,'Rectangle',crops,'Color','cyan','Linewidth',5));
title('Original');
% subplot(3,2,2); % coded measurement
axes(h(2));
imshow(meas,[]);
title('Coded frame');
for ic = 1:length(colors)
    color = colors{ic};
    crop = crops(ic,:);
    sp_truth = squeeze(sum(sum(orig(crop(2):crop(2)+crop(4)-1,crop(1):crop(1)+crop(3)-1,:)/MAXB,2),1));
    % figure;
    % subplot(3,2,ic+2);
    axes(h(ic+2));
    plot(wavelength,sp_truth/max(sp_truth),'k-','linewidth',linewidth,'markersize',markersize,'DisplayName','Ground truth'); hold on;

    for imeth = 1:length(methods)
        meth = methods{imeth};
        eval(sprintf('sp_%s = squeeze(sum(sum(v%s(crop(2):crop(2)+crop(4)-1,crop(1):crop(1)+crop(3)-1,:),2),1));',meth,meth));
        eval(sprintf('corrall(ic,imeth)=corr(sp_%s,sp_truth);',meth)); % corr is scale invariant
        eval(sprintf('plot(wavelength,rescale(sp_%s,min(sp_truth)/max(sp_truth),1),markers{imeth},''linewidth'',linewidth,''markersize'',markersize,''DisplayName'',[methnames{imeth} '', corr: %s'']);',meth,num2str(corrall(ic,imeth),'%.4f')));
    end
    title(color);
    if ic==3 || ic==4
        xlabel('Wavelength (nm)');
    end
    if ic==1 || ic==3
        ylabel('Intensity (a.u.)');
    end
    ylim([0 1]);
    hlg = legend('show');
    set(hlg,'FontSize',legendfontsize);
    legend('location',locos{ic});
    legend('boxoff');
end

%% [2.2] show exemplar frames (wavelength)
bands = [1 10 20 31];

nb = length(bands);
nmeth = length(methods);
f = figure('Position',[50 50 700 950]);
h = tight_subplot(nb,nmeth+1,[.002 .002],[.03 .03],[.02 .02]);
set(f, 'Color', 'white');
for ib = 1:nb
    band = bands(ib);
    lambda = wavelength(band);
    wlmap = (gray*kron(ones(3,1),spectrumRGB(lambda))); % colormap with the corresponding RGB wavelength
    wlmap = wlmap/max(wlmap(:));

    axes(h(1+(ib-1)*(nmeth+1)));
    imshow(uint8(orig(:,:,band)),'colormap',wlmap);
    if ib==1
        title('Ground truth');
    end
    for imeth = 1:nmeth
        meth = methods{imeth};
        axes(h(imeth+1+(ib-1)*(nmeth+1)));
        eval(sprintf('imshow(v%s(:,:,band),''colormap'',wlmap);',meth));
        if ib==1
            title(sprintf('%s',methnames{imeth}));
        end
    end
end

repo: DeSCI


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function:         PlotIt(FVr_temp,iter,S_struct)
% Author:           Rainer Storn
% Description:      PlotIt can be used for special purpose plots
%                   used in deopt.m.
% Parameters:       FVr_temp     (I)    Paramter vector
%                   iter         (I)    counter for optimization iterations
%                   S_Struct     (I)    Contains a variety of parameters.
%                                       For details see Rundeopt.m
% Return value:     -
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotIt(FVr_temp,iter,S_struct)

figure(1);
plot(S_struct.FM_pop(:,1),S_struct.FM_pop(:,2),'ro');
xlabel('x');
ylabel('y');
xlim([-10;10]);
ylim([-10;10]);
grid on;

ff = vertcat(S_struct.S_val.FVr_oa);
f1 = ff(:,1);
f2 = ff(:,2);
figure(2);
plot(f1,f2,'bo');
xlabel('f_1');
ylabel('f_2');
xlim([-100;100]);
ylim([-100;100]);
grid on;

drawnow;
pause(0.5); %wait for one second to allow convenient viewing
return

author: Ewenwan

repo: Mathematics

save_path: github-repos/MATLAB/Ewenwan-Mathematics


text:
function [Z,A,B,rss] = arch2(X, na, mode, sil)

%   [Z,A,B,rss] = arch2(X, na, mode, silent )
%
%   archetypal analysis of column orientated data set <X>
%
%   input arguments :
%
%   - each column of data is one 'observation', e.g. the sample values of
%     all channels in a multichannel measurement at one point in time
%
%   - na : number of generated archetypes
%
%   - mode can be one of the following : 'normalized' (default),
%     'mean', 'raw'
%     - in mode 'normalized' each column of data is centered by removing its mean
%       and then normalized by dividing through its standard deviation before
%       the covariance matrix is calculated
%     - in mode 'mean' only the mean of every column of data is removed
%     - in mode 'raw' no preprocessing is applied to data
%     - in mode 'scale' <X> is divided by max(abs(X))
%
%   - silent is an optional flag which supresses output of text and plot on the matlab
%     screen. Returned values (see below) are in no way affected
%
%
%   output arguments :
%
%   - Z : each column of Z is an archetype
%
%   - A : the columns of A are the coefficients of the archetypes to
%         the constrains ||X-Z*A|| -> min
%
%   - B : the columns of B create the X-mixtures (the archtypes)
%         Z=X*B
%
%   - rss : the residual sum of squares for each iteration
%
%   Christian Merkwirth & Joerg Wichard
%   Februar 1998

global silent

narginchk(2,3);

if nargin < 3
    mode = 'normalized';
end

if nargin < 4
    silent = 0;
else
    silent = 1;
end

[m,n] = size(X);

printline('archetypal analysis')
printline(['on data set of size ' num2str(n) 'x' num2str(m)]);

if (na < 1)
  printline('number of archetypes must be greater than zero');
end

if strncmp(mode, 'r',1)
  mode = 'raw';
  printline('no data preprocessing');
elseif strncmp(mode, 'm',1)
  mode = 'mean';
  printline('removing mean from data set');
  mn = mean(X,2);
  X = X - repmat(mn, 1, n);
elseif strncmp(mode, 's',1)
  mode = 'scale';
  printline('scaling data set');
  xm = max(max(abs(X)));
  X = (1/xm)*X;
else
  mode = 'normalized';
  printline('removing mean and normalizing data');
  mn = mean(X,2);
  X = X - repmat(mn, 1, n);
  dv = std(X,0,2);
  X = X ./ repmat(dv, 1, n);
end

gew1=20*m;              % Gewichtung der Convexit�tsbedingung
gew2=5*m;               % Gewichtung der Convexit�tsbedingung
tol=0.01;              %%Toleranz f�r die Abbruchbedingung
numb = 20;              %% Max. Anzahl der Iterationen

%%  x zuf�llig ausw�len
B=eye(n);
rp=randperm(na);
B=B(:,rp);
Z=X*B;

%% Hier beginnt die Alternierende Optimierung
rss(1,:)=[ 0 sum(sum(X .* X))]
count=0;

for c=1:numb;

  for l=1:na;

    %% Maximum an Z anf�gen
    MX=max(max(Z));
    Z(m+1,:)=gew1*MX;
    X(m+1,:)=gew1*MX;

    %% A suchen bei konstantem Z
    for i=1:n;
      A(:,i)= lsqnonneg(Z,X(:,i));
    end;

    %% Gewichtung entfernen
    X=X(1:m,:);
    Z=Z(1:m,:);

    %% Z suchen bei konstantem A
    for i=1:n;
      V(:,i)=A(l,i)*(X(:,i) - Z*A(:,i) + A(l,i)*Z(:,l));
    end;

    %% Singul�re Archetypen durch max ||X-Z*A|| ersetzen
    if ( (sum(A(l,:)) .* sum(A(l,:))) ==0)
      VT=sum((X-Z*A).*(X-Z*A));
      [VTC,VTI]=max(VT);
      B(:,l)=0;
      B(VTI,l)=1;

    else
      VS=(sum(A(l,:).*A(l,:)))*sum(V,2);
      MV=max(max(X));
      X(m+1,:)=gew2*MV;
      VS(m+1)=gew2*MV;

      B(:,l)=nnls(X,VS, 4 * max(size(X)) * norm(X,1) * eps);

      %% Gewichtung entfernen
      X=X(1:m,:);
      VS=VS(1:m);
    end;
  end;

  Z=X*B;

  %% norm RSS berechnen
  rss((c+1),:)=[ c sum(sum((X-Z*A).*(X-Z*A))) ];

  if (abs(rss(c+1,2)-rss(c,2)) < tol*rss(c,2) )
    break;
  end;

  count=count+1;
end;

%% Ergebnis auf Konsistenz �berpr�fen
[C,I]=max(rss);
if( (I < count) & (rss(I) < rss(count)) )
  printline('Minimum not reached')
end

if( count == numb )
  printline('Maximum Number of Iterations')
end

if silent
  figure(1)
  subplot(4,1,1)
  plot(rss(:,1), rss(:,2))
  title('Residual Sum of Squares')

  subplot(4,1,2)
  plot(X)
  title('Input Data');

  subplot(4,1,3)
  plot((Z),'k')
  title('Archetypes');

  subplot(4,1,4)
  plot((Z),'k')
  hold on;
  plot(X)
  hold off;
  title('Archetypes & Input Data');
end

function printline(string)
global silent
if silent~=1
        disp(string)
end


text:
function [varargout]=patchEdgeAngles(F,V)

% function [A,Ad]=patchEdgeAngles(F,V)
% -----------------------------------------------------------------------
% Computes the edge angles (A) for the patch data specified by the faces
% (F) and vertices (V) arrays.
%
%
% See also:
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2015/08/18
%------------------------------------------------------------------------

%%

numNodes=size(F,2);

A=zeros(size(F));
for q=1:1:numNodes

    q1=q;
    q0=q-1;
    if q0<1
        q0=numNodes;
    end
    q2=q+1;
    if q2>numNodes
        q2=1;
    end

    P0=V(F(:,q0),:)-V(F(:,q1),:);
    P0=vecnormalize(P0);
    P2=V(F(:,q2),:)-V(F(:,q1),:);
    P2=vecnormalize(P2);
    A(:,q)=acos(dot(P0,P2,2));

end

varargout{1}=A;

if nargout>1
    varargout{2}=180*(A./pi);
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

repo: GIBBON


text:
function outcl = merge_clusters(c2m, subcl)
% Function for synchronizing the field list of cluster structures
% and merging them
%
% :Usage:
% ::
%
%    subclusters = merge_clusters(clusters_to_match, subclusters_to_change)
%

    if isempty(c2m)
        disp('Merge clusters: Clusters to match is empty.');
        outcl = subcl;
        return
    end

    if isempty(subcl)
        outcl = c2m;
        return
    end

    if ~iscell(subcl)
        s{1} = subcl;
        subcl = s;
    end
    outcl = c2m;

    for i = 1:length(subcl)
        for j = 1:length(subcl{i})
            % make all field names match.
            newcl{i}(j) = c2m(1);
            N = fieldnames(newcl{i}(j));
            for NN = 1:length(N)
                eval(['newcl{i}(j).' N{NN} ' = [];'])
            end

            %N = fieldnames(clusters(i));
            for NN = 1:length(N)
                str = ['newcl{i}(j).' N{NN} ' = subcl{i}(j).' N{NN} ';'];
                try
                    eval(str)
                catch
                    disp(['Making ' N{NN}])
                    str = ['newcl{i}(j).' N{NN} ' = [];'];
                    eval(str)
                end
            end
        end

        outcl = ([outcl, newcl{i}]);
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Template Matlab script to create an BIDS compatible _bold.json file
% This example lists all required and optional fields.
% When adding additional metadata please use CamelCase
% Use version of DICOM ontology terms whenever possible.
%
% DHermes, 2017

% Writing json files relies on the JSONio library
% https://github.com/bids-standard/bids-matlab
%
% Make sure it is in the matab/octave path
try
    bids.bids_matlab_version;
catch
    warning('%s\n%s\n%s\n%s', ...
            'Writing the JSON file seems to have failed.', ...
            'Make sure that the following library is in the matlab/octave path:', ...
            'https://github.com/bids-standard/bids-matlab');
end

%%

clear;

this_dir = fileparts(mfilename('fullpath'));
root_dir = fullfile(this_dir, '..', filesep, '..');

project = 'templates';

sub_label = '01';
ses_label = '01';
task_label = 'FullExample';
run_label = '01';

name_spec.modality = 'func';
name_spec.suffix = 'bold';
name_spec.ext = '.json';
name_spec.entities = struct('sub', sub_label, ...
                            'ses', ses_label, ...
                            'task', task_label, ...
                            'run', run_label);

% using the 'use_schema', true
% ensures that the entities will be in the correct order
bids_file = bids.File(name_spec, 'use_schema', true);

% Contrust the fullpath version of the filename
json_name = fullfile(root_dir, project, bids_file.bids_path, bids_file.filename);

%% Adding metadata
% to get the definition of each metadata,
% you can use the bids.Schema class from bids matlab

% For example
schema = bids.Schema;
def = schema.get_definition('TaskName');
fprintf(def.description);

%% Required fields
json.TaskName = 'FullExample';
json.RepetitionTime = [];

%% Required fields if using a fieldmap
json.PhaseEncodingDirection = '';
json.EffectiveEchoSpacing = [];
json.EchoTime = [];

%% Recommended fields:

%% Scanner Hardware metadata fields
json.Manufacturer = '';
json.ManufacturersModelName = '';
json.MagneticFieldStrength = [];
json.DeviceSerialNumber = '';
json.StationName = ' ';
json.SoftwareVersions = ' ';
json.HardcopyDeviceSoftwareVersion = '';
json.ReceiveCoilName = ' ';
json.ReceiveCoilActiveElements = '';
json.GradientSetType = '';
json.MRTransmitCoilSequence = '';
json.MatrixCoilMode = '';
json.CoilCombinationMethod = '';

%% Sequence Specifics metadata fields
json.PulseSequenceType = '';
json.ScanningSequence = '';
json.SequenceVariant = '';
json.ScanOptions = '';
json.SequenceName = '';
json.PulseSequenceDetails = '';
json.NonlinearGradientCorrection = '';

%% In-Plane Spatial Encoding metadata fields
json.NumberShots = [];
json.ParallelReductionFactorInPlane = '';
json.ParallelAcquisitionTechnique = '';
json.PartialFourier = '';
json.PartialFourierDirection = '';
json.WaterFatShift = '';
json.EchoTrainLength = '';

%% Timing Parameters metadata fields
json.InversionTime = [];
json.SliceEncodingDirection = '';
json.DwellTime = [];
json.TotalReadoutTime = [];
json.DelayAfterTrigger = [];
json.NumberOfVolumesDiscardedByScanner = [];
json.NumberOfVolumesDiscardedByUser = [];

%% RF & Contrast metadata field
json.FlipAngle = [];

%% Slice Acceleration metadata field
json.MultibandAccelerationFactor = [];

%% Task metadata field
json.Instructions = '';
json.TaskDescription = '';
json.CogAtlasID = 'https://www.cognitiveatlas.org/FIXME';
json.CogPOID = 'http://www.cogpo.org/ontologies/CogPOver1.owl#FIXME';

%% Institution information metadata fields
json.InstitutionName = '';
json.InstitutionAddress = '';
json.InstitutionalDepartmentName = '';

%% Write JSON
% Make sure the directory exists
bids.util.mkdir(fileparts(json_name));
bids.util.jsonencode(json_name, json);

author: bids-standard

repo: bids-starter-kit