text:
function [LS,LD,I] = isolines(V,F,S,iso)
  % ISOLINES compute a list of isolines for a scalar field S defined on the
  % mesh (V,F)
  %
  % [LS,LD,I] = isolines(V,F,S,iso)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle indices
  %   S  #V list of scalar values defined on V
  %   iso #iso list of iso values
  % Outputs:
  %   LS  #L by dim list of isoline start positions
  %   LD  #L by dim list of isoline destination positions
  %   I  #L list of indices into iso revealing corresponding iso value
  %
  % alternative: tracing isolines should result in smaller plot data (every
  % vertex only appears once
  %

  % make sure iso is a ROW vector
  iso = iso(:)';
  % number of isolines
  niso = numel(iso);

  % number of domain positions
  n = size(V,1);
  % number of dimensions
  dim = size(V,2);

  % number of faces
  m = size(F,1);

  % Rename for convenience
  S1 = S(F(:,1),:);
  S2 = S(F(:,2),:);
  S3 = S(F(:,3),:);

  % t12(i,j) reveals parameter t where isovalue j falls on the line from
  % corner 1 to corner 2 of face i
  t12 = bsxfun(@rdivide,bsxfun(@minus,iso,S1),S2-S1);
  t23 = bsxfun(@rdivide,bsxfun(@minus,iso,S2),S3-S2);
  t31 = bsxfun(@rdivide,bsxfun(@minus,iso,S3),S1-S3);

  % replace values outside [0,1] with NaNs
  t12( (t12<-eps)|(t12>(1+eps)) ) = NaN;
  t23( (t23<-eps)|(t23>(1+eps)) ) = NaN;
  t31( (t31<-eps)|(t31>(1+eps)) ) = NaN;

  % masks for line "parallel" to 12
  l12 = ~isnan(t23) & ~isnan(t31);
  l23 = ~isnan(t31) & ~isnan(t12);
  l31 = ~isnan(t12) & ~isnan(t23);

  % find non-zeros (lines) from t23 to t31
  [F12,I12,~] = find(l12);
  [F23,I23,~] = find(l23);
  [F31,I31,~] = find(l31);
  % indices directly into t23 and t31 corresponding to F12 and I12
  %ti12 = sub2ind(size(l12),F12,I12);
  %ti23 = sub2ind(size(l23),F23,I23);
  %ti31 = sub2ind(size(l31),F31,I31);
  % faster sub2ind
  ti12 = F12+(I12-1)*size(l12,1);
  ti23 = F23+(I23-1)*size(l23,1);
  ti31 = F31+(I31-1)*size(l31,1);

  % compute actual position values
  LS = [ ...
    ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t23(ti12), V(F(F12,2),:)) + ...
    bsxfun(@times,  t23(ti12), V(F(F12,3),:)) ...
    ; ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t31(ti23), V(F(F23,3),:)) + ...
    bsxfun(@times,  t31(ti23), V(F(F23,1),:)) ...
    ;... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t12(ti31), V(F(F31,1),:)) + ...
    bsxfun(@times,  t12(ti31), V(F(F31,2),:)) ...
    ];
  LD = [...
    ... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t31(ti12), V(F(F12,3),:)) + ...
    bsxfun(@times,  t31(ti12), V(F(F12,1),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t12(ti23), V(F(F23,1),:)) + ...
    bsxfun(@times,  t12(ti23), V(F(F23,2),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t23(ti31), V(F(F31,2),:)) + ...
    bsxfun(@times,  t23(ti31), V(F(F31,3),:)) ...
    ];

  % I is just concatenation of each I12
  I = [I12;I23;I31];

end

author: gpeyre


text:
function MS = createMesh1D(varargin)
% MeshStructure = createMesh1D(Nx, Width)
% MeshStructure = createMesh1D(facelocationX)
% creates a uniform 1D mesh:
% Nx is the number of cells in x (horizontal) direction
% Width is the domain length in x direction
%
% SYNOPSIS:
%   MeshStructure = createMesh1D(Nx, Width)
%
% PARAMETERS:
%   Nx: number of cells in the domain
%   Width: domain length
%
% RETURNS:
%   MeshStructure.
%                 dimensions=1 (1D problem)
%                 numbering: shows the indexes of cellsn from left to right
%                 dx: cell size (=Width/Nx)
%                 cellcenters.x: location of each cell in the x direction
%                 facecenters.x: location of interface between cells in the
%                 x direction
%                 numberofcells: Nx
%
%
% EXAMPLE:
%   L = 1.0; % length of the domain
%   Nx = 10; % number of cells in the domain
%   m = createMesh1D(Nx, L);
%   plot(m.cellcenters.x, ones(size(m.cellcenters.x)), 'o', ...
%        m.facecenters.x, ones(size(m.facecenters.x)), '-+');
%   legend('cell centers', 'face centers');
%
% SEE ALSO:
%     createMesh2D, createMesh3D, createMeshCylindrical1D, ...
%     createMeshCylindrical2D, createCellVariable, createFaceVariable

% Written by Ali A. Eftekhari
% See the license file

if nargin==2
  % uniform 1D mesh
  Nx=varargin{1};
  Width=varargin{2};
  % cell size is dx
  dx = Width/Nx;
  CellSize.x= dx*ones(Nx+2,1);
  CellSize.y= [0.0];
  CellSize.z= [0.0];
  CellLocation.x= [1:Nx]'*dx-dx/2;
  CellLocation.y= [0.0];
  CellLocation.z= [0.0];
  FaceLocation.x= [0:Nx]'*dx;
  FaceLocation.y= [0.0];
  FaceLocation.z= [0.0];
elseif nargin==1
  % nonuniform 1D mesh
  facelocationX=varargin{1};
  n=size(facelocationX);
  if n(1)==1
      facelocationX=facelocationX';
  end
  Nx = length(facelocationX)-1;
  CellSize.x= [facelocationX(2)-facelocationX(1); ...
    facelocationX(2:end)-facelocationX(1:end-1); ...
    facelocationX(end)-facelocationX(end-1)];
  CellSize.y= [0.0];
  CellSize.z= [0.0];
  CellLocation.x= 0.5*(facelocationX(2:end)+facelocationX(1:end-1));
  CellLocation.y= [0.0];
  CellLocation.z= [0.0];
  FaceLocation.x= facelocationX;
  FaceLocation.y= [0.0];
  FaceLocation.z= [0.0];
end

MS=MeshStructure(1, ...
  [Nx,1], ...
  CellSize, ...
  CellLocation, ...
  FaceLocation, ...
  [1], ...
  [1]);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [groups, orphans, R, C] = connectedComponents(model, type, figures, files)
% Assuming two reactions are connected if they share metabolites, calculates the connected components
% in the stoichiometric matrix, that is, the sets of reactions that share a set of metabolites
%
% USAGE:
%
%    [groups, orphans, R, C] = connectedComponents(model, type, figures)
%
% INPUT:
%    model:
%
% OPTIONAL INPUTS:
%    type:       {('allComponents'), 'largestComponent'}
%    figures:    Will generate plots of the grouping algorithm as it creates block diagonal
%                groups in from top left to bottom right in `W`.
%    files:      Indicator, whether several files containing indicator
%                matrices are generated.
% OUTPUTS:
%    groups:     a structure array (the number of distinct groups is length(groups)) with fields:
%
%                  * `groups(i).num_els` - number of reactions in group `i`.
%                  * `groups(i).block` - sub-block identity of group `i`.
%                  * `groups(i).elements` - reactions of W that are in group `i`.
%                  * `groups(i).degrees` - degrees of connection for each reaction in group `i`.
%    orphans:    elements of W that were not in any group, becasue they did not meet the constraints.
%    R:          reaction adjacency
%    C:          compound adjacency
%
% All components require:
% Connected Component Analysis on an Undirected Graph by Tristan Ursell
% http://www.mathworks.com/matlabcentral/fileexchange/35442-connected-component-analysis-on-an-undirected-graph.
%
% Largest component requires:
% gaimc : Graph Algorithms In Matlab Code by David Gleich
% http://www.mathworks.com/matlabcentral/fileexchange/24134-gaimc-graph-algorithms-in-matlab-code.
%
%  .. Author:
%        - Ronan Fleming, 2012
%        - Thomas Pfau May 2017, Speedup and addition of files indicator

if ~exist('type','var')
    type='allComponents';
end
if ~exist('figures','var')
    figures=0;
end
if ~exist('files','var')
    files=0;
end

model=findSExRxnInd(model);

%stoichiometric matrix
S=model.S;
%dont include exchange reactions
S(:,~model.SIntRxnBool)=0;

[m,n]=size(S);

%binary form
B=sparse(m,n);
B(S~=0)=1;

%Compound adjacency
C1=B*B';

%number of reactions a species participates in
nReactionsSpeciesParticipatesIn=diag(C1,0);

%take out connections by cofactors
[nReactionsSpeciesParticipatesInSorted,IX] = sort(nReactionsSpeciesParticipatesIn,'descend');
%model.mets(IX(1:80))

% %omit reactions connected by cofactors
% omitMet=false(m,1);
% for i=1:m
%     metAbbr=model.mets{i};
%     if strcmp(metAbbr(1:2),'h[')
%         omitMet(i)=1;
%         continue;
%     end
%     if strcmp(metAbbr(1:3),'k[')
%         omitMet(i)=1;
%         continue;
%     end
%     if length(metAbbr)>3
%         if strcmp(metAbbr(1:3),'pi[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:3),'cl[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:3),'o2[')
%             omitMet(i)=1;
%             continue;
%         end
%     end
%
%     if length(metAbbr)>4
%         if strcmp(metAbbr(1:4),'na1[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'h2o[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'co2[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'atp[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'adp[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'utp[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'gtp[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'gdp[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'amp[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'nad[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'fad[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'coa[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'ppi[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'nh4[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'ACP[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'thf[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:4),'crn[')
%             omitMet(i)=1;
%             continue;
%         end
%     end
%
%     if length(metAbbr)>5
%         if strcmp(metAbbr(1:5),'nadh[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:5),'fadh[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:5),'nadp[')
%             omitMet(i)=1;
%             continue;
%         end
%     end
%
%     if length(metAbbr)>6
%         if strcmp(metAbbr(1:6),'nadph[')
%             omitMet(i)=1;
%             continue;
%         end
%         if strcmp(metAbbr(1:6),'accoa[')
%             omitMet(i)=1;
%             continue;
%         end
%     end
%
% end
%
% %omit these metabolites
% B(omitMet,:)=0;

%Reaction adjacency
R1=B'*B;

%number of species in a reaction
nMolecularSpeciesInReaction=diag(R1,0);

R=R1;
%Set the diagonal to 0
R(logical(eye(size(R1,1)))) = 0;

if files
    R2=triu(R);
    fid=fopen('reactionAdjacencyOtherThanCofactors.txt','w');
    for j=1:n
        fprintf(fid,'%s\t',model.rxns{j});
    end
    fprintf(fid,'\n');
    for j=1:n
        fprintf(fid,'%s\t',model.rxns{j});
        for k=1:n
            fprintf(fid,'%u\t',full(R2(j,k)));
        end
        fprintf(fid,'\n');
    end
    fclose(fid);
end

C=C1;
%Set the diagonal to 0
C(logical(eye(size(C1,1)))) = 0;

if strcmp(type,'largestComponent')
    if ~exist('largest_component')
        error('Install gamic and add it to your path. (http://www.mathworks.com/matlabcentral/fileexchange/24134-gaimc-graph-algorithms-in-matlab-code)')
    end
    [Acc,p] = largest_component(R);
    degree=sum(Acc);
    groups(1).num_els=nnz(degree);
    groups(1).block='largest';
    groups(1).elements=find(p);
    groups(1).degrees=degree;
    orphans=[];
else
    if figures==1
        [groups,orphans]=graph_analysis(R,'plot',1);
    else
        [groups,orphans]=graph_analysis(R);
    end

end

if files
    fid=fopen('reactionsNotConnectedByAnything.txt','w');
    bool=model.SIntRxnBool;
    bool(groups.elements) = 0;
    for j=1:n
        if bool(j)
            fprintf(fid,'%s\n',model.rxns{j});
        end
    end
    fclose(fid);
end

author: opencobra

repo: cobratoolbox


text:
function [C,nc] = graph_coloring(A,nc)
  % GRAPH_COLORING Color a graph given an adjacency matrix. This heuristic is
  % designed to efficiently find a conservative coloring (i.e., using too many
  % colors). It does *not* attempt to find the minimal coloring. It is best used
  % when the provided nc is a few more than the optimal coloring. For example,
  % planar graphs can be colored with 4 colors, (slow) heuristics for finding the
  % optimal coloring will often find 6 colors, but this function will very
  % quickly find a 7-coloring (nc=7). On such a graph, if you passed nc=4, this
  % function will likely waste time attemping to find a 4-coloring, 5-coloring,
  % 6-coloring, and then very quickly find a 7-coloring. It is designed to give
  % up searching for parsimonious colorings and increase the number of colors
  % (throwing a warning).
  %
  % [C,nc] = graph_coloring(A,nc)
  %
  % Inputs:
  %   A  #V by #V adjacency matrix
  %   nc  desired number of colors (for planar/triangle meshes, 7 is a fast
  %     choice; for tetrahedral meshes, 13 appears to often be fast)
  % Outputs:
  %   C  #V list of color ids into (1:nc)
  %   nc  effective number of colors >= input nc
  %
  % Examples:
  %   C = graph_coloring(facet_adjacency_matrix(F),9);
  %   tsurf(F,V,'CData',C);
  %   colormap(cbrewer('Set1',max(C)));

  n = size(A,1);
  [FI,FJ] = find(triu(A,1));
  % always marking the higher valence vertex as bad definitely makes things
  % worse 50x.
  %
  % if nc is much greater than optimal, always marking the lower valence vertex
  % doesn't much effect. But for nc close to optimal, this seems to make a very
  % big difference (300x)
  val = sum(A,2);
  swap = val(FI)<val(FJ);
  EI = FI.*swap + FJ.*~swap;
  EJ = FJ.*swap + FI.*~swap;

  onc = nc;
  % worst case just increase number of colors
  max_outer = 10;
  for outer = 1:max_outer
    % Might as well retry a few times for this nc
    for retry = 1:10
      C = ceil(nc*rand(n,1));
      for iter = 1:2*ceil(sqrt(size(A,1)))
        bad = unique(EI(C(EI)==C(EJ)));
        if numel(bad) == 0
          break;
        end
        C(bad) = nan;
        C(bad) = ceil(nc*rand(numel(bad),1));
      end
      if numel(bad) == 0
        break;
      end
    end
    if numel(bad) == 0
      break;
    end
    nc = nc+1;
    if nc>max_outer
      error('Failed to converge.\n');
    end
  end
  if nc ~= onc
    warning('Had to increasing number of colors')
  end

  %V2C = sparse(1:n,C,1,n,nc);
  % find all edges with the same color
  %A*V2C;
end

author: alecjacobson

repo: gptoolbox


text:
function exact = p14_exact ( )

%*****************************************************************************80
%
%% P14_EXACT returns the estimated integral for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 1.0634618101722400407;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% update_vel_pos(p,n);
function update_vel_pos()
global vel El El_lim_up El_lim_low
global n w k_pb_i R iteration k_pb_f  k_gb_i k_gb_f k_pb k_gb c_gb c_pb delta_t;
global NumVer p pbest_loc gbest_loc
R = n /iteration;
k_pb = k_pb_i + R*(k_pb_f - k_pb_i);
k_gb = k_gb_i + R*(k_gb_f - k_gb_i);
c_pb = rand*k_pb;
c_gb = rand*k_gb;

      vel(1,1:NumVer,p) = w*vel(1,1:NumVer,p) + c_pb*(pbest_loc(1,1:NumVer,p) - El(1,1:NumVer,p)) + c_gb*(gbest_loc(1,1:NumVer) - El(1,1:NumVer,p));

%Velocity Upper Bound Velocity Lower Bound
     for chgdir=1:NumVer
      if((vel(1,chgdir,p) > (El_lim_up(1,chgdir) - El(1,chgdir,p))))
          vel(1,chgdir,p) = -0.5*vel(1,chgdir,p);%(El_lim_up(1,1:NumVer) - El(1,1:NumVer,p));%             x(i)=x(i)+2*v(i);
      end
      if((vel(1,chgdir,p) < (El_lim_low(1,chgdir) - El(1,chgdir,p))))
          vel(1,chgdir,p) = -0.5*vel(1,chgdir,p);%(El_lim_up(1,1:NumVer) - El(1,1:NumVer,p));
      end
     end

 El(1,1:NumVer,p)= El(1,1:NumVer,p)+delta_t*vel(1,1:NumVer,p);
end

author: Sable

repo: mcbench-benchmarks


text:
function ai = r8cvv_rget ( mn, a, m, roff, i )

%*****************************************************************************80
%
%% R8CVV_RGET gets row I from an R8CVV.
%
%  Discussion:
%
%    An R8CVV is a "vector of vectors" of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer MN, the size of the cell array.
%
%    Input, real A(MN), the cell array.
%
%    Input, integer M, the number of rows in the array.
%
%    Input, integer ROFF(M+1), the row offsets.
%
%    Input, integer I, the row.
%    1 <= I <= M.
%
%    Output, real AI(NR(I)), the value of A(I,*).
%
  k1 = roff(i) + 1;
  k2 = roff(i+1);
  nv = k2 + 1 - k1;
  ai(1:nv) = a(k1:k2);

  return
end


text:
function [policy, EU] = upot_to_opt_policy(pot)
% UPOT_TO_OPT_POLICY Compute an optimal deterministic policy given a utility potential
% [policy, EU] = upot_to_opt_policy(pot)
%
% policy(a,b, ..., z) = P(do z | a, b, ..), which will be a delta function
% EU is the contraction of this potential, i.e., P .* U

sz = pot.sizes; % mysize(pot.p);
if isempty(sz)
  EU = pot.u;
  policy = [];
  return;
end

parent_size = prod(sz(1:end-1));
self_size = sz(end);
C = pot.p .* pot.u; % contraction
C = reshape(C, parent_size, self_size);
policy = zeros(parent_size, self_size);
for i=1:parent_size
  act = argmax(C(i,:));
  policy(i, act) = 1;
end
policy = myreshape(policy, sz);
EU = sum(C(:));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [modes, weights,centerId] = calcComponents(SO3F,varargin)
% heuristic to find modal orientations
%
% Syntax
%   [modes, volume] = calcComponents(SO3F)
%   [modes, volume, centerId] = calcComponents(SO3F,'seed',ori)
%
% Input
%  SO3F - @SO3Fun
%  ori - initial list of @orientation
%
% Output
%  modes     - modal @orientation
%  volume    - volume of the component
%  centerId  - list of ids to which each initial ori converged to
%
% Options
%  resolution - search-grid resolution
%  angle      - maximum component width used for volume computation
%  exact      - do not dismiss very small modes at the end
%
% See also
% SO3Fun/max

% extract options
maxIter = get_option(varargin,'maxIter',100);
res = get_option(varargin,'resolution',0.05*degree);
omega = 1.5.^(-7:1:4) * degree;
omega(omega<res) = [];
omega = [0,omega];
tol = get_option(varargin,'tolerance',0.5*degree);
maxAngle = get_option(varargin,{'radius','angle'},inf);

% initial seed
if check_option(varargin,'seed')
  seed = reshape(get_option(varargin,'seed'),[],1);
  weights = get_option(varargin,'weights',SO3F.eval(seed));
elseif isa(SO3F,'SO3FunRBF')
  seed = SO3F.center;
  weights = SO3F.weights;
else
  seed = equispacedSO3Grid(SO3F.CS,SO3F.SS,'resolution',2.5*degree);
  weights = ones(length(seed),1) ./ length(seed);
end
id = weights>0;
seed = reshape(seed(id),[],1);
weights = weights(id);
weights = weights ./ sum(weights);

centerId = 1:length(seed);
modes = seed;

% join orientations if possible
%[modes,~,id2] = unique(modes,'tolerance',tol);
%centerId = id2(centerId);
%weights = accumarray(id2,weights);

finished = false(size(modes));

G = SO3F.grad;

for k = 1:maxIter
  progress(k,maxIter,' finding ODF components: ');

  % gradient
  g = normalize(G.eval(modes(~finished)),1);

  % prepare for linesearch
  line_ori = exp(repmat(modes(~finished),1,length(omega)),g * omega);

  % evaluate along lines
  line_v = SO3F.eval(line_ori);

  % take the maximum
  [v_max(~finished),id] = max(line_v,[],2);

  % update orientions
  modes(~finished) = line_ori(sub2ind(size(line_ori),(1:length(g)).',id));

  %nnz(id>1)
  if all(id == 1), break; end

  % join orientations if possible
  [~,~,id2] = unique(modes,'tolerance',tol);

  modes = modes(maxVote(id2,v_max(~finished)));

  centerId = id2(centerId);
  finished = accumarray(id2,finished,[],@any);
  if maxAngle == inf, weights = accumarray(id2,weights); end
  length(modes);

end

% accumulate only weights that are sufficiently close to the centers
if maxAngle < inf
  inRadius = angle(seed,modes(centerId))<maxAngle;
  weights = accumarray(centerId(inRadius), weights(inRadius));
end

% sort components according to volume
[weights,id] = sort(weights,'descend');
modes = modes(id);
iid(id) = 1:length(id);
centerId = iid(centerId);

if ~check_option(varargin,'exact')
  id = weights > min([0.01, numSym(SO3F.CS) * maxAngle^3 ./ 8*pi^2,0.5 * max(weights)]);
  weights = weights(id);
  modes = modes(id);
  ids = 1:length(id);
  centerId(ismember(centerId,ids(~id))) = 0;
end

% weights = [2 1 5 3 4] -> [1 2 3 4 5]
% id -> [2 1 5 3 4]
% centerId = 3 -> 5
end

function test
% testing code

cs = crystalSymmetry('432');
cs2 = specimenSymmetry;
center = orientation.rand(5,cs,cs2);
odf = unimodalODF(center,'halfwidth',5*degree)
ori = discreteSample(odf,2000);
odf2 = calcDensity(ori,'noFourier','exact','halfwidth',2.5*degree)

cs2 = crystalSymmetry('432')
center = orientation.rand(5,cs,cs2);
odf = unimodalODF(center,'halfwidth',2.5*degree)
ori = discreteSample(odf,1000);
odf2 = calcDensity(ori,'noFourier','exact','halfwidth',3*degree)

[modes,vol,cId] = odf2.calcComponents;

for i = 1:length(modes)
  plot(ori(cId==i),'axisAngle')
  hold on
  plot(modes(i),'MarkerFaceColor','k','MarkerSize',10)
end
hold off

min(angle_outer(modes,center) ./ degree)

end

author: mtex-toolbox

repo: mtex

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


text:
function raw_y = adpcm_decoder(adpcm_y)

% This m-file is based on the app note: AN643, Adaptive differential pulse
% code modulation using PICmicro microcontrollers, Microchip Technology
% Inc. The app note is avaialbe from www.microchip.com
% Example:  Y = wavread('test.wav');
%           y = adpcm_encoder(Y);
%           YY = adpcm_decode(y);
IndexTable = [-1, -1, -1, -1, 2, 4, 6, 8, -1, -1, -1, -1, 2, 4, 6, 8];

StepSizeTable = [7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 19, 21, 23, 25, 28, 31, 34, 37, 41, 45, 50, 55, 60, 66, 73, 80, 88, 97, 107, 118, 130, 143, 157, 173, 190, 209, 230, 253, 279, 307, 337, 371, 408, 449, 494, 544, 598, 658, 724, 796, 876, 963, 1060, 1166, 1282, 1411, 1552, 1707, 1878, 2066, 2272, 2499, 2749, 3024, 3327, 3660, 4026, 4428, 4871, 5358, 5894, 6484, 7132, 7845, 8630, 9493, 10442, 11487, 12635, 13899, 15289, 16818, 18500, 20350, 22385, 24623, 27086, 29794, 32767];

prevsample = 0;
previndex = 1;

Ns = length(adpcm_y);
n = 1;

while (n <= Ns)
    predsample = prevsample;
    index = previndex;
    step = StepSizeTable(index);
    code = adpcm_y(n);

    diffq = bitshift(step, -3);
    if (bitand(code, 4))
        diffq = diffq + step;
    end
    if (bitand(code, 2))
        diffq = diffq + bitshift(step, -1);
    end
    if (bitand(code, 1))
        diffq = diffq + bitshift(step, -2);
    end

    if (bitand(code, 8))
        predsample = predsample - diffq;
    else
        predsample = predsample + diffq;
    end

    if (predsample > 32767)
        predsample = 32767;
    elseif (predsample < -32768)
        predsample = -32768;
    end

    index = index + IndexTable(code+1);

    if (index < 1)
        index = 1;
    end
    if (index > 89)
        index = 89;
    end

    prevsample = predsample;
    previndex = index;

    raw_y(n) = predsample / 32767;
    n = n + 1;
end


text:
function to_be_implemented(mfile)
%TO_BE_IMPLEMENTED throws an error indicating missing code
%
%   Usage: to_be_implemented([mfilename])
%
%   Input parameters:
%       mfilename   - string containing the name of the calling m file.
%                   NOTE: this variable is already available in all m files as
%                   mfilename!

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ====================================
nargmin = 0;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Main ============================================================
if exist('mfile','var')
    isargchar(mfile)
    error('%s: This functionality has to be implemented yet!',upper(mfile));
else
    error('This functionality has to be implemented yet!');
end

author: sfstoolbox

repo: sfs-matlab


text:
function [km] = ly2km(ly)
% Convert length from light years to kilometers.
% Chad A. Greene 2012
km = ly*9460730472581;


text:
function y = cvx_isconcave( x, full )
persistent remap
if isempty( remap ),
	remap = cvx_remap('concave');
end
y = remap( cvx_classify( x ) );
if nargin < 2 || ~full,
	y = all( y );
else
	y = reshape( y, x.size_ );
end

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

author: yu-jiang

repo: radpbook

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


text:
function exact = p02_exact ( )

%*****************************************************************************80
%
%% P02_EXACT returns the exact integral for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = sqrt ( pi );

  return
end


text:
function [centered_data,centroid] = CenterRowData(input)
% each row is a data point

centroid = mean(input);

centered_data = input-(ones(size(input, 1), 1)*centroid);

author: zzstefan

repo: BrainNetClass


text:
function [out] = random(msg,BSID,DIUC,Frame);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                       %
%%     Name: random.m                                                    %
%%                                                                       %
%%     Description: This function is realisez uniquely by the way        %
%%      defined in the standard. It is one more stage in the encoding.   %
%%                                                                       %
%%                                                                       %
%%     Parameters:                                                       %
%%       msg  --> Sequence of bits to encode with the algorithm BSID,    %
%%       DIUC, Frame--> Identifiers of the connection that is carried    %
%%       out.(Base Station, Downlink and the number of frames)           %
%%                                                                       %
%%                                                                       %
%%     Result: It gives back the chain of bits encoded accordin to       %
%%      the need (following DIUC, BSID and Frame). The encoding and      %
%%      decoding are an identical process.                               %
%%                                                                       %
%%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

seed = zeros (1,15);
seed=[de2bi(BSID,4,'left-msb') 1 1 de2bi(DIUC,4,'left-msb') 1 de2bi(Frame,4,'left-msb')];

out=zeros(1,length(msg));

for i=1:length(msg)
    next = xor(seed(14),seed(15));
    out(i) = xor( msg(i),next);
    seed = [next seed(1,1:14)];
end

return;


text:
function value = gud ( x )

%*****************************************************************************80
%
%% GUD evaluates the Gudermannian function.
%
%  Discussion:
%
%    The Gudermannian function relates the hyperbolic and trigonometric
%    functions.  For any argument X, there is a corresponding value
%    G so that
%
%      sinh(x) = tan(g).
%
%    The value G is called the Gudermannian of X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the Gudermannian.
%
%    Output, real VALUE, the value of the Gudermannian.
%
  value = 2.0 * atan ( tanh ( 0.5 * x ) );

  return
end


text:
function ZTemp = mysympower(vartable,Z);
% MYSYMPOWER --- Fast symbolic power.
%
% ZSYM = mysympower(VARTABLE,ZNUM)
%

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%

% 03/01/02 - SP

ZTemp = '[';
for j = 1:size(Z,1)
    idx = find(Z(j,:));
    charexpr = [];
    for i = idx
        if Z(j,i)>1
            charexpr = [charexpr,char(vartable(i)),'^',int2str(Z(j,i)),'*'];
        else
            charexpr = [charexpr,char(vartable(i)),'*'];
        end;
    end;
    if isempty(idx)
        charexpr = '1*';
    end;
    charexpr = [charexpr(1:end-1),';'];
    ZTemp = [ZTemp,charexpr];
end;
ZTemp = [ZTemp(1:end-1),']'];
ZTemp = sym(ZTemp);


text:
function bratubvp
%BRATUBVP  Exercise for Example 1 of the BVP tutorial.
%   The BVP  y'' + exp(y) = 0, y(0) = 0 = y(1) is a standard example
%   of a problem with two solutions.  It is easy enough to solve, but
%   some experimentation with the guess may be necessary to get both.

% Copyright 2002, The MathWorks, Inc.

options = bvpset('stats','on');
solinit = bvpinit(linspace(0,1,5),[0.1 0]);
sol1 = bvp4c(@bratuode,@bratubc,solinit,options);

fprintf('\n');

% Change the initial guess to converge to a different solution.
solinit = bvpinit(linspace(0,1,5),[3 0]);
sol2 = bvp4c(@bratuode,@bratubc,solinit,options);

figure
plot(sol1.x,sol1.y(1,:),sol2.x,sol2.y(1,:))
title('Bratu''s equation has two solutions when \lambda = 1.')
xlabel('x')
ylabel('y')

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

function dydx = bratuode(x,y)
%BRATUODE  ODE function for the exercise of Example 1 of the BVP tutorial.
dydx = [  y(2)
         -exp(y(1))];

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

function res = bratubc(ya,yb)
%BRATUBC  Boundary conditions for the exercise of Example 1 of the BVP tutorial.
res = [ya(1)
       yb(1)];


text:
function [type,center,theta,foci,ecc,normform] = conic_properties(A,B,C,D,E,F)

[type,x0,y0,theta,A2,B2,C2,D2,E2,F2] = conic_standard_form(A,B,C,D,E,F);

% change from translation followed by rotation to rotation followed by
% translation
% x <- (x + u)*R
% = x*R + u*R
R = [cos(theta),sin(theta);-sin(theta),cos(theta)];
center = [x0,y0]*R;
normform = [A2,B2,C2,D2,E2,F2];

if strcmpi(type,'circle')
  foci = [0,0];
  ecc = 0;
elseif strcmpi(type,'ellipse') || strcmpi(type,'hyperbola'),
  a2 = 1 / min(abs(A2),abs(C2));
  b2 = 1 / max(abs(A2),abs(C2));
  c2 = sqrt(a2 - b2);
  if A2 < C2,
    foci = [-c2,0;c2,0];
  else
    foci = [0,-c2;0,c2];
  end
  ecc = c2 / a2;
elseif strcmpi(type,'parabola'),
  if abs(A2) > abs(C2),
    foci = [0,-E2/4];
  else
    foci = [-D2/4,0];
  end
  ecc = 1;
else % other
  foci = zeros(0,2);
  ecc = 0;
end

nfoci = size(foci,2);
if nfoci > 0,
  % translate
  foci(:,1) = foci(:,1) + x0;
  foci(:,2) = foci(:,2) + y0;
  % rotate
  foci = foci*R;
end

author: kristinbranson

repo: JAABA


text:
% EASY3	  Read RINEX navigation file, version 3.03, and reformat
%             into a Matlab Eph matrix. Open a RINEX observation file,
%             version 3.03, analyse the header and identify observation
%             types. The call fgetl finds the information for the epoch time.
%             Next we read observations line  by line. Finally recpo_ls
%             estimates the  (stand alone) receiver position.

% Kai Borre 31-10-2001
% Copyright (c) by Kai Borre
% $Revision: 1.0 $  $Date: 2001/10/31  $
% Total revision 2.0, February 13, 2016

% RINEX version 3.03

% Read RINEX ephemerides file and convert to internal Matlab format
rinexe('log_24h.15n','eph.dat');
% rinexe('D:\Work\Laboratory\Trana\rec\r302_long\Bas_log_2016_05_07_12.00.00.16N','eph.dat');
Eph = get_eph('eph.dat');

% Open the observation file
ofile2 = 'log_r.15o'; % log_24h.15o
% ofile2 = 'D:\Work\Laboratory\Trana\rec\r302_long\Bas_log_2016_05_07_12.00.00.16o';
fid2 = fopen(ofile2,'rt');

% The selection of observation type is set
ss = 'C1W'%;
linjer = 0;

while 1			   % Gobbling the header
    linjer = linjer +1;
    line = fgetl(fid2);
    answer = strfind(line,'END OF HEADER');
    if  ~isempty(answer), break; end;
    if (line == -1), eof = 1; break; end;
    answer = strfind(line,'ANT # / TYPE');
    if ~isempty(answer)
        delta = textscan(fid2,'%6.4f','Delimiter','\n');
        delt = delta{1};
        delx = delt(1,1);
        dely = delt(2,1);
        delz = delt(3,1);
    end

    answer = strfind(line,'SYS / # / OBS TYPES');
    if ~isempty(answer)
        tline1 = strsplit(line);
        line = fgetl(fid2);
        tline2 = strsplit(line);
        tt = horzcat(tline1,tline2);
        i = strcmp(tt,ss); % if tt equals ss, i =1, else 0
        ii = find(i == 1) ;
        ii = ii-2;
        if ii > 15,  ii = ii-7; end;
        % the cell array of strings tline1 originally contains two strings
        % which describe the system and number of observation types.
        % Both tline1 and tline2 terminates with six additional strings.
        % An extra string appears at  the start of tline2; it originates
        % from concatenation of  the two lines. The indexing does not
        % change, even if you empty the cells. They remain empty
        % cells and keep a place
        obs_col = ii%;
    end;
    answer = strfind(line,'INTERVAL');
    if ~isempty(answer)
        interval = strtok(line);
        int = str2double(interval);
    end;
end % end reading header

% the string arrays for the tline1 and tline2 contain an integer after the
% carrier phase observation. We account for this by the following
% correctional table
if         strcmp(ss(2:3), '1C'), obs_col = obs_col +1;
elseif   strcmp(ss(2:3), '1W'), obs_col = obs_col+2;
elseif  strcmp(ss(2:3), '2X'), obs_col = obs_col +3;
else    strcmp(ss(2:3), '2W'), obs_col = obs_col +4;
end

Pos = [];
epoch = 0;
dt = [];
Tline = [];

while ~feof(fid2)
    epoch = epoch +1;
    %time = 0;
    sats = [];
    sats0 = [];

    % We read the first line in every  epoch and get sow and
    % number of SVs.
    [time,post] = textscan(fid2,'%s %d8','Delimiter','\n');
    tid = time{1}{1};
    year = str2double(tid(3:6));
    month = str2double(tid(8:9));
    day = str2double(tid(11:12));
    hour = str2double(tid(14:15));
    minute = str2double(tid(17:18));
    second = str2double(tid(20:29));
    static = str2double(tid(31:32));
    NoSvs = str2double(tid(34:36));
    dte  = str2double(tid(38:56));
    dt = [dt dte];
    h = hour+minute/60+second/3600;
    jd = julday(year, month, day, h);
    [~, sec_of_week] = gps_time(jd);
    time = sec_of_week; % sow

    Obs = zeros(NoSvs,length(obs_col));
    for i = 1:NoSvs
        obs = textscan(fid2,'%s %d8','Delimiter','\n');
        obsy = obs{1}{1};
        obs = strsplit(obsy);
        sat = obs{1};
        sats(i,:) = str2double(sat(2:3));
        Obs(i,1) = str2double(obs(obs_col));
    end

    % Next we test if all observed sats have an ephemeris.
    % sats contains the SVs as read in the observation lines.
    % The intersect command delivers Sats in sorted order!
    % Therefore we must be careful in the follwing manipulations
    Sats = intersect(sats,Eph(1,:));

    %The command ismember does not change the sequence of entries in sats
    lia = ismember(sats,Sats);
    % A 0 (zero) in lia indicates that a SV has been deleted. We delete the
    % corresponding row in the observations
    sats(lia==0) = [];
    Obs(lia==0) = [];
    % All book-keeping has prepared the data so that we can nake the call
    % for a  final position computation
    pos = recpo_ls(Obs,sats,time,Eph);
    Pos = [Pos pos];
end % while
fclose(fid2);

me = mean(Pos,2);
fprintf('\n\nMean Position as Computed From %d Epochs:', epoch)
fprintf('\n\nX: %12.3f  Y: %12.3f  Z: %12.3f\n\n', me(1,1), me(2,1), me(3,1))

figure(1);
plot(1:epoch,(Pos(1,:)-Pos(1,1)*ones(1,epoch))','-',...
    1:epoch,(Pos(2,:)-Pos(2,1)*ones(1,epoch))','-.',...
    1:epoch,(Pos(3,:)-Pos(3,1)*ones(1,epoch))','--','linewidth',.25)
title('Positions over time','fontsize',16)
legend('X','Y','Z')
xlabel('Epochs [1 s interval]','fontsize',16)
ylabel('Changes in {\itX}, {\itY}, {\itZ} since the first epoch [m]','fontsize',16)
set(gca,'fontsize',16)
legend
print -dpdf easy31

figure(2);
plot(1:epoch,dt*10^6')
ylabel('Receiver clock off-set in {\mu}s','fontsize',16),
xlabel('Epochs [1 s interval]','fontsize',16)
print -dpdf  easy32

% transformation from (X, Y, Z) to geographical coordinates
[phi, lambda, h] = cart2geo(me(1,1), me(2,1), me(3,1), 5)
%%%%%%%%%%%%%%%%%%%%% end easy3.m %%%%%%%%%%%%%%%

author: yandld


text:
function sphere_llt_grid_test ( )

%*****************************************************************************80
%
%% SPHERE_LLT_GRID_TEST tests the SPHERE_LLT_GRID library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( )
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_LLT_GRID_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SPHERE_LLT_GRID library.\n' );

  sphere_llt_grid_point_count_test ( );
  sphere_llt_grid_points_test ( );
  sphere_llt_grid_line_count_test ( );
  sphere_llt_grid_lines_test ( );
  sphere_llt_grid_display_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_LLT_GRID_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function c = agc(c,varargin)

% This function applies automatic gain control (AGC) to each trace. This
% process, commonly used in seismic reflection processing, applies a
% variable scale to each trace such that the amplitude along the entire
% trace is roughly uniform. By minimizing amplitude variations along the
% trace, well-correlated but low-amplitude phases become more visible. A
% time window may be specified to control how tightly this scaling is
% applied. Use a longer window for lower frequecy signals. The function
% returns a correlation object.
%
% EXAMPLES:
%  c=agc(c)    				apply agc using the default time window (0.5 s)
%  c=agc(c,0.8)            	apply agc using a window of 0.8 s

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if (length(varargin) >= 2)
    error('Too many inputs');
end

if ~isa(c,'correlation')
    error('First input must be a correlation object');
end

if length(varargin)==1
	agcwin = varargin{1};
else
	agcwin = 0.5;
end;
agcsamp = round( agcwin * get(c.W(1),'Fs') );

% LOOP THROUGH TRACES APPLYING GAIN
for tracenum = 1:length(c.W)
   w = get(c.W(tracenum),'DATA');
   scale=zeros( length(w)-2*agcsamp , 1 );
   for index=-1*agcsamp:agcsamp
      scale=scale + abs( w(agcsamp+index+1:agcsamp+index+length(scale)) );
   end;
   scale = scale/mean(abs(scale));
   scale = [ones(agcsamp,1)*scale(1) ; scale ; ones(agcsamp,1)*scale(end)];
   w = w./scale;
   c.W(tracenum) = set(c.W(tracenum),'DATA',w);
end;

author: geoscience-community-codes

repo: GISMO


text:
% @authors:     Fuxin Li
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

function s_feat = compute_all_relevant_features( sp_seg, edge_graph )
% somehow, all the rest are computed on a halfed image, while the last one
% is on full image...
    s_feat = zeros(max(sp_seg(:)), 6);
%    sp_seg_rszd = imresize(sp_seg,0.5,'nearest');
    [cent, bb, area, perim, secondorder] = superpix_regionprops(uint16(sp_seg));
    bb = double(bb);
    bb(:,1:2) = bb(:,1:2) - 0.5;
    bb(:,3:4) = bb(:,3:4) + 0.5;
    % Left and top both have -0.5 as per MATLAB convention
    s_feat(:,1:2) = double(bb(:,1:2));
    % Eccentricity is a bit tricky...
    % Now the secondorder information is non-central, to get centralized
    % moments we need to use (x - centx)^2 = x^2 - 2 x*centx + centx^2
    % but \sum x /N = centx, therefore we get (x-centx)^2 = secondorder_x - 2
    % centx * centx * N + N * centx * centx = secondorder_x - N * centx *
    % centx
    %1/12 is the normalized second central moment of a pixel with unit length.
    uxx = secondorder(:,1) ./ single(area) - cent(:,1).^2 + 1/12;
    uyy = secondorder(:,2) ./ single(area) - cent(:,2).^2 + 1/12;
    % for uxy, we compute (x - centx) (y-centy)  = xy - centx y - centy x -
    % centx centy, using the same idea we can get that is xy - N centx
    % centy
    uxy = secondorder(:,3) ./ single(area) - cent(:,1) .* cent(:,2);
    common = sqrt((uxx - uyy).^2 + 4*uxy.^2);
    MajorAxisLength =  2*sqrt(2)*sqrt(uxx + uyy + common);
    MinorAxisLength = 2*sqrt(2)*sqrt(uxx + uyy - common);
    Eccentricity = 2*sqrt((MajorAxisLength/2).^2 - (MinorAxisLength/2).^2) ./ MajorAxisLength;
    s_feat(:,3) = Eccentricity;
    % Extent, area / bounding box area
    s_feat(:,4) = single(area) ./ single((bb(:,4) - bb(:,2)) .* (bb(:,3) - bb(:,1)));
    % Perimeter
    s_feat(:,5) = perim;
    s_feat(:,6) = sum(edge_graph,2);
    % Left here commented for verification of correctness
%     for i=1:max(sp_seg(:))
% %          s = regionprops(sp_seg==i,'BoundingBox','Eccentricity','EquivDiameter','Extent','Perimeter');
% %          s = s(1);
% %          s_feat(i,1:5) = [s.BoundingBox(1) s.BoundingBox(2) s.Eccentricity s.Extent s.Perimeter];
% %         % Inter-contour energy, change this just to superpixel boundary
%         % energy
%         all_bw = imdilate(sp_seg==i, ones(5,5)) & imdilate(~(sp_seg==i), ones(5,5));
%         s_feat(i,6) = sum(pb_thin(all_bw));
%     end
end

author: Cloud-CV

repo: object-proposals


text:
function [fitness, class, Simil] = calcfitness(data, ideals, y, w)
[nc, v_dim] = size(ideals);
d_dim = size(data,1);

if nargin==3
   w=ones(1,v_dim);
end
[class, Simil] = classifier(data(:,1:v_dim), ideals, y, w);
fitness = length(find(class-data(:, v_dim +1) == 0))/d_dim;


text:
function y = ipermute( x, order )

%   Disciplined convex/geometric programming information for IPERMUTE:
%      IPERMUTE(A,ORDER) imposes no convexity restrictions on A. ORDER
%      must be constant.

%
% Determine the permutation
%

s = x.size_;
ndxs = 1 : prod( s );
ndx2 = ipermute( reshape( ndxs, s ), order );

%
% Permute the data
%

b = x.basis_;
try
    b = b( :, ndx2 );
catch
    ndxs( ndx2( : ).' ) = ndxs;
    [ r, c, v ] = find( b );
    b = sparse( r, ndxs( c ), v, size( b, 1 ), size( b, 2 ) );
    clear r c v
end
y = cvx( size( ndx2 ), b );

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


text:
function pot = mpot(members, sizes, logp, mu, Sigma)
% MPOT Make a moment Gaussian potential.
% pot = mpot(members, sizes, logp, mu, Sigma)
%
% All params default to 0 if omitted.

n = sum(sizes);
if nargin < 3, logp = 0; end
if nargin < 4, mu = zeros(n,1); end
if nargin < 5, Sigma = zeros(n,n); end

pot.domain = members;
pot.sizes = sizes;
pot.logp = logp;
pot.mu = mu;
pot.Sigma = Sigma;zeros(n,n);
pot = class(pot, 'mpot');


text:
function pascal_car_grammar(dotrainval, testyear)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% Set configuration override
global VOC_CONFIG_OVERRIDE;
if isempty(VOC_CONFIG_OVERRIDE)
  VOC_CONFIG_OVERRIDE = @voc_config_car_grammar;
end

cls = 'car';
conf = voc_config();
cachedir = conf.paths.model_dir;
testset = conf.eval.test_set;

if nargin < 1
  dotrainval = false;
end

if nargin < 2
  % which year to test on -- a string, e.g., '2007'.
  testyear = conf.pascal.year;
end

timestamp = datestr(datevec(now()), 'dd.mmm.yyyy:HH.MM.SS');

% set the note to the training time if none is given
if nargin < 3
  note = timestamp;
end

% record a log of the training and test procedure
diary(conf.training.log([cls '-' timestamp]));

th = tic;
model = pascal_train_car_grammar(note);
toc(th);
% Free feature vector cache memory
fv_cache('free');

% lower threshold to get high recall
model.thresh = min(conf.eval.max_thresh, model.thresh);
model.interval = conf.eval.interval;

ds = pascal_test(model, testset, testyear, testyear);
ap1 = pascal_eval(cls, ds, testset, testyear, testyear);
%[ap1, ap2] = bboxpred_rescore(cls, testset, testyear);

fprintf('AP = %.4f (without bounding box prediction)\n', ap1)
%fprintf('AP = %.4f (with bounding box prediction)\n', ap2)

% Clear the override
clearvars -global VOC_CONFIG_OVERRIDE;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function u = sphs(Pf,K,t)
% Schroeder-phased harmonic sequence
% Pf = 1; % Fundamental period
% K = 10;
% u = zeros(size(time_test));
u = zeros(size(t));
for i = 1:K
    theta = 2*pi/K * sum([1:i]);
    u = u + sqrt(2/(K))*cos(2*pi * i*t/Pf + theta) ;
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function y = norms( varargin )

%NORMS   Computation of multiple vector norms.
%   NORMS( X ) provides a means to compute the norms of multiple vectors
%   packed into a matrix or N-D array. This is useful for performing
%   max-of-norms or sum-of-norms calculations.
%
%   All of the vector norms, including the false "-inf" norm, supported
%   by NORM() have been implemented in the NORMS() command.
%     NORMS(X,P)           = sum(abs(X).^P).^(1/P)
%     NORMS(X)             = NORMS(X,2).
%     NORMS(X,inf)         = max(abs(X)).
%     NORMS(X,-inf)        = min(abs(X)).
%   If X is a vector, these computations are completely identical to
%   their NORM equivalents. If X is a matrix, a row vector is returned
%   of the norms of each column of X. If X is an N-D matrix, the norms
%   are computed along the first non-singleton dimension.
%
%   NORMS( X, [], DIM ) or NORMS( X, 2, DIM ) computes Euclidean norms
%   along the dimension DIM. NORMS( X, P, DIM ) computes its norms
%   along the dimension DIM.
%
%   Disciplined convex programming information:
%       NORMS is convex, except when P<1, so an error will result if these
%       non-convex "norms" are used within CVX expressions. NORMS is
%       nonmonotonic, so its input must be affine.

persistent P
if isempty( P ),
    P.map = cvx_remap( { 'constant' ; 'l_convex' ; ...
        { 'p_convex', 'n_concave', 'affine' } } );
    P.funcs = { @norms_1, @norms_1, @norms_2 };
    P.zero = 0;
    P.reduce = true;
    P.reverse = false;
    P.constant = 1;
    P.fname = 'norms';
    P.dimarg = 3;
end
[ sx, x, p, dim ] = cvx_get_dimension( varargin, 3 );
if nargin < 2 || isempty(p),
    p = 2;
elseif ~( isnumeric(p) && numel(p)==1 && isreal(p) && p >= 1 ),
    cvx_throw( 'Second argument must be a scalar between 1 and +Inf, inclusive.' );
end
if sx(dim) == 0,
    sx(dim) = 1;
    y = zeros( sx, 1 );
    if isa( x, 'cvx' ), y = cvx( y ); end
elseif sx(dim) == 1,
    y = abs( x );
elseif p == 1,
    y = sum( abs( x ), dim );
elseif p == Inf,
    y = max( abs( x ), [], dim );
else
    y = cvx_reduce_op( P, x, p, dim );
end

function y = norms_1( x, p )
y = sum( abs( x ) .^ p, 1 ) .^ ( 1 / p );

function y = norms_2( x, p ) %#ok
[nx,nv] = size(x);
if p == 2,
    cvx_begin
        epigraph variable y( 1, nv ) nonnegative_
        { linearize(x), y } == lorentz( [ nx, nv ], 1, ~isreal( x ) ); %#ok
    cvx_end
else
    cvx_begin
        variable z( nx, nv )
        epigraph variable y( 1, nv ) nonnegative_
        if isreal(x), cmode = 'abs'; else cmode = 'cabs'; end
        { cat( 3, z, repmat(y,[nx,1]) ), linearize(x) } ...
            == geo_mean_cone( [nx,nv,2], 3, [1/p,1-1/p], cmode ); %#ok
        sum( z ) == y; %#ok
    cvx_end
end

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


text:
function this = unlink(this)
% Unlinks the object from the data file
% FORMAT this = unlink(this)
% _________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: unlink.m 5025 2012-10-31 14:44:13Z vladimir $

this.data = [];
this      = check(this);

author: spm


text:
function metabModel = extractMetModel(model, metabNames, nLayers, allCompFlag, nRxnsMetThr)
% Creates a subnetwork model around one or more metabolites
%
% USAGE:
%
%    metabModel = extractMetModel(model, metabNames, nLayers, allCompFlag, nRxnsMetThr)
%
% INPUTS:
%    model:          COBRA model structure
%    metabNames:     Metabolites to build subnetwork model around
%    nLayers:        Number of layers
%    allCompFlag:    Use all metabolites regardless of compartment
%
% OPTIONAL INPUT:
%    nRxnsMetThr:    Ignore metabolites which appear in more than `nRxnMetThr` (Default = 100)
%
% OUTPUT:
%    metabModel:     COBRA model around one or more metabolites
%
% .. Author: - Markus Herrgard 3/1/06
%              Uri David Akavia 26-Mar-2017

if (nargin < 5)
    nRxnsMetThr = 100;
end

% Filter out high degree metabolites
nRxnsMet = full(sum(model.S~=0,2));
baseMetNames = parseMetNames(model.mets);
for i = 1:length(baseMetNames)
    nRxnsMetComp(i) = sum(nRxnsMet(strcmp(baseMetNames,baseMetNames{i})));
end
nRxnsMetComp = nRxnsMetComp';
selLowDegMet = nRxnsMetComp < nRxnsMetThr;
model.S(~selLowDegMet,:) = 0;
%model.mets = model.mets(selLowDegMet); Commented out, because other fields
%aren't shrunk, which can lead to mismatch later

if (~iscell(metabNames))
    tmpMetName = metabNames;
    clear metabNames;
    metabNames{1} = tmpMetName;
end

if (allCompFlag)
    allMetNames = parseMetNames(model.mets);
    metabNames = parseMetNames(metabNames);
else
    allMetNames = model.mets;
end

metabNames = intersect(metabNames, allMetNames(selLowDegMet));
selMets = find(ismember(allMetNames,metabNames));

metS = model.S(selMets,:);
[nMet,tmp] = size(metS);
if (nMet > 1)
    selRxns = any(full(metS) ~= 0)';
elseif (~nMet)
	selRxns = false(size(model.rxns));
else
    selRxns = (full(metS) ~= 0)';
end
for i = 1:nLayers+1
    metS = model.S(selMets,:);
    [nMet,tmp] = size(metS);
    if (nMet > 1)
        selRxns = any(full(metS) ~= 0)';
    elseif (~nMet)
        selRxns = false(size(model.rxns));
	else
        selRxns = (full(metS) ~= 0)';
    end

    if (isfield(model,'c'))
        selRxns = selRxns & ~ (model.c == 1);
    end
    nextLayerMets = find(any(model.S(:,selRxns) ~= 0,2));
    selMets = union(selMets,nextLayerMets);
end

rxnNames = model.rxns(selRxns);
metNames = model.mets(selMets);

metabModel = extractSubNetwork(model,rxnNames,metNames);


text:
function[S,tau,tc] = staogram(data_spk,data_lfp,smp,plt,Tc,Tinc,Tw,w,D)
%
% staogram : calculates a moving window spike triggered ave %
%   Usage:[S,tau,tc] = staogram(data_spk,data_lfp,smp,plt,Tc,Tinc,Tw,w,D)
%
%                 ******** INPUT *********
% Note that all times have to be consistent. If data_spk
% is in seconds, so must be sig and t. If data_spk is in
% samples, so must sig and t. The default is seconds.
%
% data_spk    - strucuture array of spike times data
% data_lfp    - array of lfp data(samples x trials)
% smp         - lfp times of samples
%
% Optional...
%
% Parameter
%
% plt     'y'|'n'
%
% 'y' standard staogram
% 'n' no plot
%
% Tc = start and end times (centres)           whole trial
% Tinc = time increment between windows             0.1
% Tw = time window width                          0.3
% w = smoothing width in seconds
% D = plot sta out to on axis [D(1) D(2)] s
%
%                ******** OUTPUT ********
%  S spike triggered average
%  tau - lag
%  tc  - bin centers

% setup defaults...
if nargin < 3;error('Require spike, lfp and lfptimes ');end
[data_spk]=padNaN(data_spk); % create a zero padded data matrix from input structural array
data_spk=data_spk'; % transposes data to get it in a form compatible with Murray's routine
if nargin < 4; plt = 'y';end
if nargin < 6; Tinc = 0.1; end
if nargin < 7; Tw = 0.5;end
if nargin < 8; w = 0.01;end
if nargin < 9; D = 0.15*[-1 1]; end
if nargin < 5;
    Tc(1) = min(data_spk(:,1)) + Tw/2;
    Tc(2) = max(max(data_spk)) - Tw/2;
end

if isempty(plt); plt = 'y';end
if isempty(Tinc); Tinc = 0.1; end
if isempty(Tw); Tw = 0.5;end
if isempty(w); w = 0.01;end
if isempty(D); D = 0.15*[-1 1]; end
if isempty(Tc);
    Tc(1) = min(data_spk(:,1)) + Tw/2;
    Tc(2) = max(max(data_spk)) - Tw/2;
end

%  round to nearest tinc...

t = smp;
Tc(1) = ceil(Tc(1)/Tinc)*Tinc;
Tc(2) = floor(Tc(2)/Tinc)*Tinc;
tc = Tc(1):Tinc:Tc(2);
for tt=1:length(tc)
  T = [tc(tt)-Tw/2 tc(tt)+Tw/2];
  if tt == 1
    [SS,tau] = sta(data_spk,data_lfp,t,'y',w,T,D,0);
    S = zeros(length(tc),length(SS));
  else
    [SS,tau] = sta(data_spk,data_lfp,t,'y',w,T,D,0);
  end
  S(tt,:) = SS';
  S(tt,:) = SS';
end

if ~strcmp(plt,'n')
  imagesc(tc,tau,squeeze(S)')
  set(gca,'ydir','normal')
  xlabel('time (s)')
  ylabel('frequency (Hz)')
  colorbar;
%  axes(h)
%  line(get(h,'xlim'),conf_C*[1 1],'color','k','linewidth',5)
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function varargout = main(varargin)
% MAIN M-file for main.fig
%      MAIN, by itself, creates a new MAIN or raises the existing
%      singleton*.
%
%      H = MAIN returns the handle to a new MAIN or the handle to
%      the existing singleton*.
%
%      MAIN('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MAIN.M with the given input arguments.
%
%      MAIN('Property','Value',...) creates a new MAIN or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before main_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to main_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help main

% Last Modified by GUIDE v2.5 06-Mar-2013 13:59:04

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

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

% --- Executes just before main is made visible.
function main_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to main (see VARARGIN)

% Choose default command line output for main
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes main wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = main_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in Browse_Image.
function Browse_Image_Callback(hObject, eventdata, handles)
% hObject    handle to Browse_Image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.output = hObject;
[filename,pathname]=uigetfile('*.jpg');
if ~filename
    errordlg('Select an Image File.');
     return;
end
a=imread(strcat(pathname,filename));
handles.my_data=a;
guidata(hObject, handles);

% --- Executes on button press in Press_for_haar.
function Press_for_haar_Callback(hObject, eventdata, handles)
% hObject    handle to Press_for_haar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
a=handles.my_data;
[v,w]=haaar(a);
handles.my_out=v;
handles.w=w;
guidata(hObject,handles);
axes(handles.image_out);
axis off;
axes(handles.image_out);
imshow(handles.my_out);

% --- Executes on button press in press_for_LL_pixels.
function press_for_LL_pixels_Callback(hObject, eventdata, handles)
% hObject    handle to press_for_LL_pixels (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.my_out=handles.w;
axes(handles.image_out);
axis off;
axes(handles.image_out);
imshow(handles.my_out);

% --- Executes on button press in Close.
function Close_Callback(hObject, eventdata, handles)
% hObject    handle to Close (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close all;
clc;
clear all;


text:
function spherical_harmonic_test01 ( )

%*****************************************************************************80
%
%% SPHERICAL_HARMONIC_TEST01 tests LEGENDRE_ASSOCIATED_NORMALIZED,  LEGENDRE_ASSOCIATED_NORMALIZED_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERICAL_HARMONIC_TEST01:\n' );
  fprintf ( 1, '  LEGENDRE_ASSOCIATED_NORMALIZED evaluates the\n' );
  fprintf ( 1, '  associated Legrendre functions, using a\n' );
  fprintf ( 1, '  normalization appropriate for spherical harmonics.\n' );
  fprintf ( 1, '  LEGENDRE_ASSOCIATED_NORMALIZED_VALUES returns some exact values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      N       M    X     Exact F     PNM(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, m, x, fx ] = legendre_associated_normalized_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = legendre_associated_normalized ( n, m, x );

    fprintf ( 1, '  %6d  %6d  %6f  %12f  %12f\n', n, m, x, fx, fx2(n+1) );

  end

  return
end


text:
function g = h2g(h, epsilon, WordToReport)

% G = h2g(H,EPSILON)
%
% Hybrid-H to Hybrid-G transformation
% H and G are matrices of size [2,2,F]
% where F is the number of frequencies
% (the number of ports is always 2)
%
% EPSILON is a limit used in finding correspondent Hybrid-H matrices in the
% vicinity of singularities; by default 1e-12, should be enough for most
% realistic problems; could be increased for a gain in speed
%
% written by tudor dima, tudima@zahoo.com, change the z into y
if nargin < 3, WordToReport = 'hybrid-G'; end;
if nargin < 2, epsilon = []; end;
if isempty(epsilon), epsilon = 1e-12; end;
% 30.09.2011    - clear bug in D calculation

d = h(1,1,:).*h(2,2,:) - h(1,2,:).*h(2,1,:);
[n,i] = min(abs(d));
flagApproxG = false;
while n <= epsilon
    % 'fix' only the i-th D
    flagApproxG = true;
    p1 = 1+round(rand); p2 = 1+round(rand);
    h(p1,p2,i) = h(p1,p2,i)+(rand-0.5)*epsilon;
    % was:
    %h(1+csi,1+~csi,i) = h(1+csi,1+~csi,i)+(rand-0.5)*epsilon;
    d = h(1,1,i).*h(2,2,i) - h(1,2,i).*h(2,1,i);
    [n,i] = min(abs(d));
end

if flagApproxG
    fprintf(1,'%s%s%s\n%s\n', 'caution: correspondent ', ...
        WordToReport, ' matrix non-existent', ...
        'an approximation is produced');
end

g(1,1,:) = h(2,2,:)./d;
g(1,2,:) = -h(1,2,:)./d;
g(2,1,:) = -h(2,1,:)./d;
g(2,2,:) = h(1,1,:)./d;


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym besseljn (@var{alpha}, @var{x})
%% Symbolic Spherical Bessel function of the first kind.
%%
%% Example:
%% @example
%% @group
%% syms n x
%% A = besseljn(n, x)
%%   @result{} A = (sym) jn(n, x)
%% diff(A)
%%   @result{} ans = (sym)
%%
%%                      (n + 1)⋅jn(n, x)
%%       jn(n - 1, x) - ────────────────
%%                           x
%% @end group
%% @end example
%%
%% @seealso{@@sym/besselyn, @@sym/besselj}
%% @end defmethod

function Y = besseljn(n, x)
  if (nargin ~= 2)
    print_usage ();
  end
  Y = elementwise_op ('jn', sym(n), sym(x));
end

%!test
%! % roundtrip
%! syms x
%! A = double(besseljn(sym(2), sym(9)));
%! q = besseljn(sym(2), x);
%! h = function_handle(q);
%! B = h(9);
%! assert (abs (A - B) <= eps)

%!error jn(sym('x'))

repo: octsympy


text:
function R=q2dcm(q)
% Q2DCM(Q) converts quaternions into direction cosine matrices.
%
%     The resultant DCM(s) will perform the same transformations as the
%     quaternion(s) in Q, i.e.:
%
%       R*v = qvxform(q, v)
%
%     where R is the DCM, V is a vector, and Q is the quaternion.  Note that
%     for purposes of quaternion-vector multiplication, a vector is treated
%     as a quaterion with a scalar element of zero.
%
%     If the input, Q, is a vector of quaternions, the output, R, will be
%     3x3xN where input quaternion Q(k,:) corresponds to output DCM
%     R(:,:,k).
%
%     Note that the input Q will be processed by QNORM to ensure normality.
%
% See also DCM2Q, QNORM.

% Release: $Name: quaternions-1_2_2 $
% $Revision: 1.13 $
% $Date: 2002/01/21 06:46:20 $

% Copyright (C) 2000-02, Jay A. St. Pierre.  All rights reserved.

if nargin~=1
  error('q2dcm() requires one input argument');
else
  qtype=isq(q);
  if ( qtype == 0 )
    error(['Invalid input: must be a quaternion or a vector of' ...
          ' quaternions'])
  end
end

% Make sure input is a column of quaternions
if( qtype==1 )
  q=q.';
end

% Make sure quaternion is normalized to prevent skewed DCM
q=qnorm(q);

% Build quaternion element products
q1q1=q(:,1).*q(:,1);
q1q2=q(:,1).*q(:,2);
q1q3=q(:,1).*q(:,3);
q1q4=q(:,1).*q(:,4);

q2q2=q(:,2).*q(:,2);
q2q3=q(:,2).*q(:,3);
q2q4=q(:,2).*q(:,4);

q3q3=q(:,3).*q(:,3);
q3q4=q(:,3).*q(:,4);

q4q4=q(:,4).*q(:,4);

% Build DCM
R(1,1,:) =  q1q1 - q2q2 - q3q3 + q4q4;
R(1,2,:) = 2*(q1q2 + q3q4);
R(1,3,:) = 2*(q1q3 - q2q4);

R(2,1,:) = 2*(q1q2 - q3q4);
R(2,2,:) = -q1q1 + q2q2 - q3q3 + q4q4;
R(2,3,:) = 2*(q2q3 + q1q4);

R(3,1,:) = 2*(q1q3 + q2q4);
R(3,2,:) = 2*(q2q3 - q1q4);
R(3,3,:) = -q1q1 - q2q2 + q3q3 + q4q4;

author: christianwengert


text:
%TEST_TIMING_DGT_FAC  Test timing factorization DGTs
%
%   This script test the timing SPREADADJs by comparing the results to
%   spreadadj in the main toolbox. Therefore, the correctness of
%   spreadadj must be verified first.

routinemax=6;

spfraction=.1;

test_failed=0;

disp('--- Used subroutines ---');

for rtype=1:2

  if rtype==1
    rname='REAL ';
  else
    rname='CMPLX';
  end;

  for sptype=1:2

    if sptype==1
      spname='FULL  ';
    else
      spname='SPARSE';
    end;

    for L=12:13

      if rtype==1
        if sptype==1
          coef=rand(L,L);
        else
          coef=sprand(L,L,spfraction);
        end;
      else
        if sptype==1
          coef=crand(L,L);
        else
          coef=spcrand(L,L,spfraction);
        end;
      end;

      cadj=spreadadj(coef);

      for rout=1:routinemax
        cadj2=feval(['ref_spreadadj_',num2str(rout)],coef);

        rdiff=cadj-cadj2;

        res=norm(rdiff(:));

        fail='';
        if res>10e-10
          fail='FAILED';
          test_failed=test_failed+1;
        end;

        s=sprintf('ADJ %s %s %i L:%3i %0.5g %s',rname,spname,rout,L,res,fail);
        disp(s)
      end;

    end;

  end;

end;

test_failed

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function TS_Compute(doParallel,ts_id_range,op_id_range,computeWhat,customFile,howVocal)
% TS_Compute    Computes missing elements of TS_DataMat
%
%---EXAMPLE USAGE:
% TS_Compute;
%
%---INPUTS:
% doParallel: if true, attempts to use the Parallel Computing Toolbox to run
%               computations in parallel over multiple cores.
% ts_id_range: a custom range of time series IDs to compute (default: [] -- compute all)
% op_id_range: a custom range of operation IDs to compute (default: [] -- compute all)
% computeWhat: whether to compute just missing values ('missing', default), or
% 				ALSO retry results that previously threw an error ('error'), or
% 				ALSO retry any result that previously did not return a good value ('bad')
% customFile: reads in and writes to a custom output file
% howVocal:   {'fast','minimal','full'}: how much output to show about the calculation of operations.
%
%---OUTPUTS:
% Writes output to customFile (HCTSA.mat by default)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% --------------------------------------------------------------------------
%% Check inputs and set defaults
% --------------------------------------------------------------------------

% Use Matlab's Parallel Computing toolbox?
if nargin < 1
	doParallel = false;
end

% Custom range of ts_ids to compute
if nargin < 2
    ts_id_range = []; % compute all ts_ids in the file by default
end

% Custom range of op_ids to compute
if nargin < 3
    op_id_range = []; % compute all op_ids in the file by default
end

if nargin < 4 || isempty(computeWhat)
	computeWhat = 'missing';
end
if ~ismember(computeWhat,{'missing','error','bad'})
	error('Unknown setting ''%s''',computeWhat);
end

% Custom HCTSA.mat file:
if nargin < 5 || isempty(customFile)
    customFile = 'raw';
end

% Be vocal?
if nargin < 6
    howVocal = 'full'; % Write back full information on all calculations to screen
    % prints every piece of code evaluated (nice for error checking)
end

% --------------------------------------------------------------------------
%% Load information from local files
% --------------------------------------------------------------------------
[TS_DataMat,TimeSeries,Operations,customFile] = TS_LoadData(customFile);
fileVarsStruct = whos('-file',customFile);
fileVars = {fileVarsStruct.name};
if ~all(ismember({'TimeSeries','Operations','MasterOperations','TS_DataMat'},fileVars))
	error('\nCannot compute on %s: there are missing variables',customFile);
end
MasterOperations = TS_GetFromData(customFile,'MasterOperations');
if ismember('TS_CalcTime',fileVars)
	TS_CalcTime = TS_GetFromData(customFile,'TS_CalcTime');
end
if ismember('TS_Quality',fileVars)
	TS_Quality = TS_GetFromData(customFile,'TS_Quality');
end

%-------------------------------------------------------------------------------
% Tell the user about settings
%-------------------------------------------------------------------------------
if doParallel
    fprintf(1,['Computation will be performed across multiple cores' ...
            ' using Matlab''s Parallel Computing Toolbox.\n']);
else % use single-threaded for loops
    fprintf(1,'Computations will be performed serially without parallelization.\n');
end

% ------------------------------------------------------------------------------
% Get indices if computing a subset
% ------------------------------------------------------------------------------
allIDs = TimeSeries.ID;
if isempty(ts_id_range)
    ts_id_range = allIDs;
    tsIndex = 1:height(TimeSeries);
else
    ts_id_range = intersect(ts_id_range,allIDs);
    tsIndex = find(ismember(allIDs,ts_id_range));
end
allIDs = Operations.ID;
if isempty(op_id_range)
    op_id_range = allIDs;
    opCompute = ones(1,height(Operations));
else
    op_id_range = intersect(op_id_range,allIDs);
    opCompute = ismember(allIDs,op_id_range)';
end

% Definitions
numTimeSeries = length(ts_id_range); % Number of time series
numOps = length(op_id_range); % Number of operations

% Check that some computable range exists
if numTimeSeries==0 || numOps==0
    fprintf(1,'%u time series and %u operations match the ids provided. Exiting.\n',...
						numTimeSeries,numOps);
    return
end

if (nargin < 6) && (numOps < 30)
    howVocal = 'fast';
end

%-------------------------------------------------------------------------------
fprintf(1,'[%s]: Extracting %u features from each of %u time series.\n',...
                    datestr(now),numOps,numTimeSeries);

% The times vector stores the time taken for each time series to have its
% operations calculated (for determining time remaining)
times = zeros(numTimeSeries,1);

% Initialize TS_CalcTime and TS_Quality if they don't yet exist
if ~exist('TS_CalcTime','var')
    TS_CalcTime = zeros(size(TS_DataMat));
end
if ~exist('TS_Quality','var')
    TS_Quality = zeros(size(TS_DataMat));
end

% --------------------------------------------------------------------------
%% Computation
% --------------------------------------------------------------------------
if strcmp(howVocal,'fast')
    % Just show a progress bar over time series (not across features for a given time series)
    % makes sense for, say, catch22, where computations are super fast
    BF_ProgressBar('new')
end

numCalc_all = zeros(numTimeSeries,1);

for i = 1:numTimeSeries

    tsInd = tsIndex(i);
	bigTimer = tic;

    % ----
    % Which operations need calculating for this time series?:
    % ----
	qualityLabels = TS_Quality(tsInd,:); % The calculation states of any existing results for the current time series, a line of TS_Quality
	   	  % NaN indicates a value never before calculated, 1 indicates fatal error before (try again now)

	% Determine which operations are awaiting calculation for this time series:
	switch computeWhat
	case 'missing'
		% try to compute missing values (i.e, never previously computed for this time series)
	    toCalc = (opCompute & isnan(qualityLabels));
	case 'error'
		% compute missing or previously threw an error
		toCalc = (opCompute & (isnan(qualityLabels) | qualityLabels == 1));
	case 'bad'
		% compute missing, or anything that wasn't previously a good value
		toCalc = (opCompute & (isnan(qualityLabels) | qualityLabels > 0)); % Operations awaiting calculation
	end
    numCalc = sum(toCalc); % Number of operations to evaluate
	numCalc_all(i) = numCalc; % keep a record of how many were calculated at each iteration

    % -----
    % Check that all operations have a master ID attached:
    % -----
    if sum(toCalc) < numCalc
        % Error in the database structure; some operations are missing MasterID assignment
        error('?? Database structure error: some operations have not been assigned a valid master operation');
    end

    switch howVocal
    case 'minimal'
        fprintf(1,'- - - - - - - Time series %u / %u: %s - - - - -\n',i,numTimeSeries,TimeSeries.Name{tsInd});
    case 'full'
        fprintf(1,'\n\n=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=\n');
        fprintf(1,'; ; ; : : : : ; ; ;    %s     ; ; ; : : : ; ; ;\n',datestr(now));
        fprintf(1,'- - - - - - - - - - - Time series %u / %u - - - - - - - - - - -\n',i,numTimeSeries);
    end

    if numCalc > 0 % some to calculate
		try
	        [featureVector,calcTimes,calcQuality] = TS_CalculateFeatureVector(TimeSeries(tsInd,:),...
								doParallel,Operations(toCalc,:),MasterOperations,true,howVocal);
		catch ME
			% Skip to the next time series; the entries for this time series in TS_DataMat etc. will remain NaNs
			warning('Calculation for time series %u / %u failed:\n%s',i,numTimeSeries,ME.message)
			continue
		end

        % ------------------------------------------------------------------------------
		%% Store the calculated information back to local matrices
        % ------------------------------------------------------------------------------
        TS_DataMat(tsInd,toCalc) = featureVector; % store outputs in TS_DataMat
		TS_CalcTime(tsInd,toCalc) = calcTimes; % store calculation times in TS_CalcTime
		TS_Quality(tsInd,toCalc) = calcQuality; % store quality labels in TS_Quality
		% NB: the calculation time assigned for individual operations is the total calculation
		% time taken to evaluate the master code.
    else
    	fprintf(1,'Nothing calculated! All %u operations already complete!!\n',numOps);
	end

    % The time taken to calculate (or not, if numCalc = 0) all operations for this time series:
    times(i) = toc(bigTimer);

    % Update progress to the user
    switch howVocal
    case 'fast'
        BF_ProgressBar(i/numTimeSeries)
    case 'minimal'

    case 'full'
        if i < numTimeSeries
            fprintf(1,'- - - - - - - -  %u time series remaining - - - - - - - -\n',numTimeSeries-i);
            fprintf(1,'- - - - - - - -  %s remaining - - - - - - - - -\n', ...
                                                BF_TheTime(((numTimeSeries-i)*mean(times(1:i))),1));
        else % The final time series
            fprintf(1,'- - - - - - - - - - All %u time series calculated! - - - - - - - - - -\n', ...
                                                        numTimeSeries);
        end
        fprintf(1,'********************************************************************\n');
    end

end

% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
%% Finished calculating!!
% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
if strcmp(howVocal,'fast')
    BF_ProgressBar('close')
    fprintf(1,'Calculations complete in a total of %s.\n',BF_TheTime(sum(times),1));
else
    fprintf(1,'!! !! !! !! !! !! Calculation completed !! !! !! !! !!\n');
    fprintf(1,'[%s]: Calculations complete in a total of %s.\n',...
                        datestr(now),BF_TheTime(sum(times),1));
end

% Save back to local files (if results were computed):
if any(numCalc_all > 0)
	fprintf(1,'Saving all results to %s...',customFile);
	saveTimer = tic;
	save(customFile,'TS_DataMat','TS_CalcTime','TS_Quality','-append')
	fprintf(1,' Saved in %s.\n',BF_TheTime(toc(saveTimer)));
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function f = smiley(varargin)

if nargin == 0 || ~isa(varargin{1},'vector3d')
  if check_option(varargin,'exact')
    f = S2FunHandle(@(v) S2Fun.smiley(v));
  else
    f = S2FunHarmonic.quadrature(@(v) S2Fun.smiley(v),varargin{:});
  end
  return;
end

v = varargin{1};
v = v(:)';

% Radial test function: quadratic spline
f_r = @(z,h) (z>h).*(z-h).^2./(1-h).^2;
if check_option(varargin,'even')
  f_r = @(z,h) f_r(z,h)+f_r(-z,h); % f has to be even
end

x_0 = [pi/2,0; 0.6,-0.6; 0.6,0.6; -0.5,-1; -0.5,-0.5; -0.5,0; -0.5,0.5; -0.5,1];
h_0 = [0.7; 0.96; 0.96; 0.93; 0.93; 0.93; 0.93; 0.93];
c_0 = [0.5; -0.5; -0.5; 0.25; 0.25; 0.25; 0.25; 0.25];

centers = vector3d.byPolar(x_0(:, 1), x_0(:, 2));
if strcmpi(getMTEXpref('xAxisDirection'),'east')
  centers = rotate(centers,90*degree);
end

% TODO: upper line can be replaced by lower line with Matlab 2017
fh = @(v) (sum(repmat(c_0,1,length(v)) .* f_r(dot(repmat(v,length(centers),1), repmat(centers,1,length(v))), repmat(h_0,1,length(v))), 1))';
%fh = @(v) (sum(c_0.*f_r(dot(v, centers), h_0), 1))';

f = fh(v);

end


text:
%script: copy correct dll-files

 if(exist('./pco_uint32.m','file'))
  delete('./pco_uint32.m');
 end

 if(exist('./pco_uint32err.m','file'))
  delete('./pco_uint32err.m');
 end

 if verLessThan('matlab','8.0.0')
  copyfile('./ver_7/pco_uint32.m','./pco_uint32.m');
  copyfile('./ver_7/pco_uint32err.m','./pco_uint32err.m');
 else
  copyfile('./ver_8/pco_uint32.m','./pco_uint32.m');
  copyfile('./ver_8/pco_uint32err.m','./pco_uint32err.m');
 end

 if(exist('../runtime/include','dir'))
  copyfile('../runtime/include/*.h','./');
 end

 pco_camera_create_deffile();

 if(strcmp(computer('arch'),'win32'))

  if(exist('./sc2_cam.dll','file'))
   delete('./sc2_cam.dll');
  end
  if(exist('../runtime/bin/sc2_cam.dll','file'))
   copyfile('../runtime/bin/sc2_cam.dll','./');
  end

  if(exist('./sc2_cl_me4.dll','file'))
   delete('./sc2_cl_me4.dll');
  end
  if(exist('../runtime/bin/sc2_cl_me4.dll','file'))
   copyfile('../runtime/bin/sc2_cl_me4.dll','./');
  end

  if(exist('./sc2_cl_nat.dll','file'))
   delete('./sc2_cl_nat.dll');
  end
  if(exist('../runtime/bin/sc2_cl_nat.dll','file'))
   copyfile('../runtime/bin/sc2_cl_nat.dll','./');
  end

  if(exist('./sc2_cl_mtx.dll','file'))
   delete('./sc2_cl_mtx.dll');
  end
  if(exist('../runtime/bin/sc2_cl_mtx.dll','file'))
   copyfile('../runtime/bin/sc2_cl_mtx.dll','./');
  end

  if(exist('./sc2_clhs.dll','file'))
   delete('./sc2_clhs.dll');
  end
  if(exist('../runtime/bin/sc2_clhs.dll','file'))
   copyfile('../runtime/bin/sc2_clhs.dll','./');
  end

  if(exist('./PCO_Recorder.dll','file'))
   delete('./PCO_Recorder.dll');
  end
  if(exist('../runtime/bin/PCO_Recorder.dll','file'))
   copyfile('../runtime/bin/PCO_Recorder.dll','./');
  end

  if(exist('./PCO_File.dll','file'))
   delete('./PCO_File.dll');
  end
  if(exist('../runtime/bin/PCO_File.dll','file'))
   copyfile('../runtime/bin/PCO_File.dll','./');
  end

 elseif(strcmp(computer('arch'),'win64'))

  if(exist('./sc2_cam.dll','file'))
   delete('./sc2_cam.dll');
  end
  if(exist('../runtime/bin64/sc2_cam.dll','file'))
   copyfile('../runtime/bin64/sc2_cam.dll','./');
  end

  if(exist('./sc2_cl_me4.dll','file'))
   delete('./sc2_cl_me4.dll');
  end
  if(exist('../runtime/bin64/sc2_cl_me4.dll','file'))
   copyfile('../runtime/bin64/sc2_cl_me4.dll','./');
  end

  if(exist('./sc2_cl_nat.dll','file'))
   delete('./sc2_cl_nat.dll');
  end
  if(exist('../runtime/bin64/sc2_cl_nat.dll','file'))
   copyfile('../runtime/bin64/sc2_cl_nat.dll','./');
  end

  if(exist('./sc2_cl_mtx.dll','file'))
   delete('./sc2_cl_mtx.dll');
  end
  if(exist('../runtime/bin64/sc2_cl_mtx.dll','file'))
   copyfile('../runtime/bin64/sc2_cl_mtx.dll','./');
  end

  if(exist('./sc2_clhs.dll','file'))
   delete('./sc2_clhs.dll');
  end
  if(exist('../runtime/bin64/sc2_clhs.dll','file'))
   copyfile('../runtime/bin64/sc2_clhs.dll','./');
  end

  if(exist('./PCO_Recorder.dll','file'))
   delete('./PCO_Recorder.dll');
  end
  if(exist('../runtime/bin64/PCO_Recorder.dll','file'))
   copyfile('../runtime/bin64/PCO_Recorder.dll','./');
  end

  if(exist('./PCO_File.dll','file'))
   delete('./PCO_File.dll');
  end
  if(exist('../runtime/bin64/PCO_File.dll','file'))
   copyfile('../runtime/bin64/PCO_File.dll','./');
  end

 else
  error('This platform is not supported');
 end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
% OtTablePost Post changes to the database table.
%
% Usage:
%    OtTablePost(htable)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-11 Created (Johan Meijdam, VORtech)

function iret = OtTablePost(htable)
   iret = Ot('table.post', htable);
end

author: Sable


text:
function [Ynorm, Ymean] = normalizeRatings(Y, R)
%NORMALIZERATINGS Preprocess data by subtracting mean rating for every
%movie (every row)
%   [Ynorm, Ymean] = NORMALIZERATINGS(Y, R) normalized Y so that each movie
%   has a rating of 0 on average, and returns the mean rating in Ymean.
%

[m, n] = size(Y);
Ymean = zeros(m, 1);
Ynorm = zeros(size(Y));
for i = 1:m
    idx = find(R(i, :) == 1);
    Ymean(i) = mean(Y(i, idx));
    Ynorm(i, idx) = Y(i, idx) - Ymean(i);
end

end

author: Ayatans

repo: Machine-Learning-homework


text:
function test_checkcode

% WALLTIME 00:20:00
% MEM 2gb
% DEPENDENCY

% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3309

% this is a subset of the LINT warnings, I removed the ones about unused variables
%
%     ''findstr' is not recommended. Use 'strfind' instead.'
%     ''flipdim' is not recommended. Use 'flip' instead.'
%     ''isempty(strfind(str1, str2))' is not recommended. Use '~contains(str1, str2)' instead.'
%     ''isequalwithequalnans' is not recommended. Use 'isequaln' instead.'
%     ''resampled_tsdif' produces a value that might be unused.'
%     ''skip' produces a value that might be unused.'
%     ''spectrum.welch' is not recommended. Use 'pwelch' instead.'
%     ''strfind(str1, str2)' is not recommended. Use 'contains(str1, str2)' instead.'
%     ''strread' is not recommended. Use 'textscan' instead.'
%     ''textread' is not recommended. Use 'textscan' instead.'
%     ''||' produces a value that might be unused.'
%     ''~cellfun('isempty', strfind(str1, str2))' is not recommended. Use 'contains(str1, str2)' instead.'
%     ''~cellfun(@isempty, strfind(str1, str2))' is not recommended. Use 'contains(str1, str2)' instead.'
%     ''~isempty(strfind(str1, str2))' is not recommended. Use 'contains(str1, str2)' instead.'
%     'Assign the onCleanup output argument to a variable. Do not use the tilde operator (~) in place of a variable.'
%     'Best practice is to separate output variables with commas.'
%     'Calling AXES(h) in a loop can be slow. Consider moving the call to AXES outside the loop.'
%     'Consider replacing GRIDDATA with SCATTEREDINTERPOLANT for better performance.'
%     'Consider using ISA instead of comparing the class name.'
%     'Consider using newline, semicolon, or comma before this statement for readability.'
%     'DISP(SPRINTF(...)) can usually be replaced by FPRINTF(...\n).'
%     'ERROR takes SPRINTF-like arguments directly.'
%     'EXIST with two input arguments is generally faster and clearer than with one input argument.'
%     'Extra comma is unnecessary in CASE statement before newline.'
%     'Extra comma is unnecessary in CATCH statement before newline.'
%     'Extra comma is unnecessary in ELSE statement before newline.'
%     'Extra comma is unnecessary in ELSEIF statement before newline.'
%     'Extra comma is unnecessary in END statement before newline.'
%     'Extra comma is unnecessary in FOR statement before newline.'
%     'Extra comma is unnecessary in IF statement before newline.'
%     'Extra comma is unnecessary in OTHERWISE statement before newline.'
%     'Extra comma is unnecessary in SWITCH statement before newline.'
%     'Extra comma is unnecessary in WHILE statement before newline.'
%     'Extra comma is unnecessary.'
%     'Extra semicolon is unnecessary in CASE statement before newline.'
%     'Extra semicolon is unnecessary in ELSE statement before newline.'
%     'Extra semicolon is unnecessary in ELSEIF statement before newline.'
%     'Extra semicolon is unnecessary in END statement before newline.'
%     'Extra semicolon is unnecessary in FOR statement before newline.'
%     'Extra semicolon is unnecessary in FUNCTION statement before newline.'
%     'Extra semicolon is unnecessary in IF statement before newline.'
%     'Extra semicolon is unnecessary in WHILE statement before newline.'
%     'Extra semicolon is unnecessary.'
%     'FOR might not be aligned with its matching END (line 169).'
%     'FOR might not be aligned with its matching END (line 31).'
%     'FOR might not be aligned with its matching END (line 437).'
%     'FOR might not be aligned with its matching END (line 81).'
%     'FREAD(FID,...,'*char') is more efficient than CHAR(FREAD(...)).'
%     'For improved robustness, consider replacing i and j by 1i.'
%     'For readability and performance, consider using 'newline' instead of 'sprintf('\n')'.'
%     'Function ERROR might be called with too few arguments.'
%     'IF might not be aligned with its matching END (line 172).'
%     'IF might not be aligned with its matching END (line 30).'
%     'IF might not be aligned with its matching END (line 37).'
%     'IF might not be aligned with its matching END (line 54).'
%     'INV(A)*b can be slower and less accurate than A\b. Consider using A\b for INV(A)*b or b/A for b*INV(A).'
%     'Instead of using transpose ('), consider using a different DIMENSION input argument to MAX.'
%     'Instead of using transpose ('), consider using a different DIMENSION input argument to MEAN.'
%     'Instead of using transpose ('), consider using a different DIMENSION input argument to MIN.'
%     'Instead of using transpose ('), consider using a different DIMENSION input argument to SUM.'
%     'LASTERR and LASTERROR are not recommended. Use an identifier on the CATCH block instead.'
%     'Loop index 'i' is changed inside of a FOR loop.'
%     'NARGCHK is not recommended. Use NARGINCHK without ERROR instead.'
%     'NUMEL(x) is usually faster than PROD(SIZE(x)).'
%     'Operator '+' is seldom used in a logical context.'
%     'Overloading DISPLAY is not recommended.'
%     'Parentheses are not needed in a FOR statement.'
%     'Parse error at FOR: usage might be invalid MATLAB syntax.'
%     'Possible inappropriate use of - operator. Use = if assignment is intended.'
%     'RAND or RANDN with the 'seed', 'state', or 'twister' inputs is not recommended. Use RNG instead.'
%     'STRMATCH is not recommended. Use STRCMP instead.'
%     'STRMATCH is not recommended. Use STRNCMP or VALIDATESTRING instead.'
%     'TRY might not be aligned with its matching END (line 299).'
%     'TRY statement should have a CATCH statement to check for unexpected errors.'
%     'Terminate statement with semicolon to suppress output (in functions).'
%     'Terminate statement with semicolon to suppress output (within a script).'
%     'The comparison will likely fail due to case mismatch.'
%     'The format might not agree with the argument count.'
%     'The warning with tag 'MATLAB:divideByZero' has been removed from MATLAB, so this statement has no effect.'
%     'This sparse indexing expression is likely to be slow.'
%     'This statement (and possibly following ones) cannot be reached.'
%     'This use of MAT2CELL should probably be replaced by a simpler, faster call to NUM2CELL.'
%     'To improve performance, replace ISEMPTY(FIND(X)) with ISEMPTY(FIND( X, 1 )).'
%     'Unexpected use of ' in a scalar context.'
%     'Unexpected use of '.*' in a scalar context.'
%     'Unexpected use of '.^' in a scalar context.'
%     'Unexpected use of '[' in a scalar context.'
%     'Use FIND with the 'first' or 'last' option.'
%     'Use STRCMP instead of == or ~= to compare character vectors.'
%     'Use STRCMPI(str1,str2) instead of using UPPER/LOWER in a call to STRCMP.'
%     'Use STRTRIM(str) instead of nesting FLIPLR and DEBLANK calls.'
%     'Use dynamic fieldnames with structures instead of GETFIELD.'
%     'Use dynamic fieldnames with structures instead of SETFIELD.'
%     'Use of 'nargout' in a script will be removed in a future release.'
%     'Use of brackets [] is unnecessary. Use parentheses to group, if needed.'
%     'Use one call to TEXTSCAN instead of calling STRTOK in a loop.'
%     'Using ISEMPTY is usually faster than comparing LENGTH to 0.'
%     'Using REGEXP(str, pattern, 'ONCE') is faster in this case.'
%     'Variable 'ElecPres' might be set by a nonlogical operator.'
%     'Variable 'ElecPres' might be set by a nonscalar operator.'
%     'Variable 'H', apparently a structure, is changed but the value seems to be unused.'
%     'Variable 'Npoints' might be set by a nonscalar operator.'
%     'Variable 'Nyq' might be set by a nonscalar operator.'
%     'Variable 'binwidth' might be set by a nonscalar operator.'
%     'Variable 'cfg', apparently a structure, is changed but the value seems to be unused.'
%     'Variable 'cond' might be set by a nonlogical operator.'
%     'Variable 'cond' might be set by a nonscalar operator.'
%     'Variable 'hdr', apparently a structure, is changed but the value seems to be unused.'
%     'Variable 'info', apparently a structure, is changed but the value seems to be unused.'
%     'Variable 'locktrllen' might be set by a nonscalar operator.'
%     'Variable 'nchans' might be set by a nonscalar operator.'
%     'Variable 'ncmb' might be set by a nonscalar operator.'
%     'Variable 'npad' might be set by a nonscalar operator.'
%     'Variable 'origrank' might be set by a nonscalar operator.'
%     'Variable 'surface', apparently a structure, is changed but the value seems to be unused.'
%     'WARNING takes SPRINTF-like arguments directly.'
%     'When checking if a variable is a matrix consider using ISMATRIX.'
%     '{ A{:} B } can often be replaced by [ A {B}], which can be much faster.'

%%

[v, p] = ft_version;
p = {p};
f = {};

while ~isempty(p)
  fprintf('looking in directory %s\n', p{1});
  nf = dir(p{1});
  p(1) = []; % remove this one
  for i=1:numel(nf)
    if nf(i).isdir && ~isequal(nf(i).name(1), '.') && ~isequal(nf(i).name, 'external') && ~isequal(nf(i).name, 'test')
      p{end+1} = fullfile(nf(i).folder, nf(i).name);
    else
      if ~nf(i).isdir && nf(i).name(end-1)=='.' && nf(i).name(end)=='m'
        f{end+1} = fullfile(nf(i).folder, nf(i).name);
      end
    end
  end
end

filelist = f;
clear f p

fprintf('found %d *.m files\n', numel(filelist));

%%

invalid = {
  'Use || instead of | as the OR operator in (scalar) conditional statements.'
  'Use && instead of & as the AND operator in (scalar) conditional statements.'
  'is known to MATLAB by its file name'
  'Use ISCELL instead of comparing the class to ''cell''.'
  'Use ISLOGICAL instead of comparing the class to ''logical''.'
  'Use TRUE or FALSE instead of LOGICAL(1) or LOGICAL(0).'
  'Best practice is to separate output variables with commas.'
  };

m = {};
for i=1:numel(filelist)
  f = filelist{i};
  s = checkcode(f);
  m = union(m, {s.message});
  for j=1:numel(s)
    for k=1:numel(invalid)
      if contains(s(j).message, invalid{k})
        % pretty display
        fprintf('================================================================================\n');
        fprintf('ERROR: %s\n', f); % invalid{k});
        fprintf('================================================================================\n');

        checkcode(f)
        % error('please fix line %d in %s', s(j).line, f);
      end
    end
  end
end

for i=1:numel(filelist)
  f = filelist{i};
  assert(checktilde(f), sprintf('Do not use tilde to ignore outputs in %s', f));
end

%%

% this shows the full list of warning messages
if false
  disp(m);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION, see https://github.com/fieldtrip/fieldtrip/issues/831
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bool = checktilde(filename)
[p, f, x] = fileparts(filename);
if ismember(f, {'matlab2markdown', 'markdown2matlab'})
  % these are known false positives, they should not trigger an error
  bool = true;
  return
end
fid = fopen(filename, 'r');
str = fread(fid, [1, inf], 'char=>char');
bool = true; % this is the return value when all is fine
bool = bool & isempty(regexp(str, '\[~ ', 'once'));
bool = bool & isempty(regexp(str, '\[~,', 'once'));
bool = bool & isempty(regexp(str, '~\]', 'once'));
bool = bool & isempty(regexp(str, ', *~ *,', 'once'));
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip