text:
function BrainNet=SR(BOLD,lambda)
% FC Network Construction using Sparse Representation (SR)
%
% Inpute:
% BOLD. A cell array with size of N x 1, each cell is a matrix of BOLD signals (#time points x #ROIs) from one of N subject
%       Each subject may have different #time points but should have the same #ROIs
% lambda. parameter for sparsity
%
% Output:
% BrainNet. FC network (#ROIs x #ROIs x #Subjects)
%
% Requires SLEP toolbox: http://www.yelab.net/software/SLEP/
%
% By Yu Zhang, zhangyu0112@gmail.com
% IDEA lab https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, UNC Chapel Hill

%% Initialize parameters for SLEP toolbox
opts=[];
opts.init=2;    % starting point: starting from a zero point here
opts.tFlag=0;   % termination criterion
opts.nFlag=0;   % normalization option: 0-without normalization
opts.rFlag=1;   % regularization % the input parameter 'rho' is a ratio in (0, 1)
opts.rsL2=0;    % the squared two norm term in min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * ||x||_1
opts.mFlag=0;   % treating it as compositive function
opts.lFlag=0;   % Nemirovski's line search

%% Construct FC networks using sparse representation
[nTime,nROI]=size(BOLD{1});
nSubj=length(BOLD);
BrainNet=zeros(nROI,nROI,nSubj,'single');
for i=1:nSubj
    temp=BOLD{i};
    temp=temp-repmat(mean(temp,2),1,nROI);     % centralization
    tempNet=zeros(nROI,nROI);
    for j=1:nROI
        ndic=setdiff(1:nROI,j);
        y=temp(:,j);
        A=temp(:,ndic);
        x=LeastR(A,y,lambda,opts);
        tempNet(ndic,j)=x;
    end
    tempNet=(tempNet+tempNet')/2;              % form symmetric adjancency matrix (main diag are zeros)
    tempNet=tempNet-diag(diag(tempNet));       % ignore self connections
    BrainNet(:,:,i)=tempNet;
end

author: zzstefan

repo: BrainNetClass


text:
clear all
addpath('../ParNMPC/')
%% Formulate an OCP using Class OptimalControlProblem

% Create an OptimalControlProblem object
OCP = OptimalControlProblem(2,... % dim of inputs
                            6,... % dim of states
                            7,... % dim of parameters
                            48);  % N: num of discritization grids

% Give names to x, u, p
[X,Theta1,Theta2,dX,dTheta1,dTheta2] = ...
    OCP.setStateName({'X','Theta1','Theta2','dX','dTheta1','dTheta2'});
[F,slack] = ...
    OCP.setInputName({'F','slack'});

% Set the prediction horizon T
OCP.setT(1.5);

% Set the dynamic function f
g  = 9.81;
m0 = 1.0;
m1 = 0.8;
m2 = 0.5;
L1 = 0.3;
L2 = 0.45;
d1 = m0+m1+m2;
d2 = (0.5*m1+m2)*L1;
d3 = 0.5*m2*L2;
d4 = (1/3*m1+m2)*L1*L1;
d5 = 0.5*m2*L1*L2;
d6 = 1/3*m2*L2*L2;
f1 = (0.5*m1+m2)*L1*g;
f2 = 0.5*m2*L2*g;
Ds = [d1              d2*cos(X)        d3*cos(Theta1);...
      d2*cos(X)       d4               d5*cos(X-Theta1);...
      d3*cos(Theta1)  d5*cos(X-Theta1) d6];
Cs = [0 -d2*sin(X)*dX      -d3*sin(Theta1)*dTheta1;...
      0  0                    d5*sin(X-Theta1)*dTheta1;...
      0 -d5*sin(X-Theta1)*dX  0];
Gs = [0;...
     -f1*sin(X);...
     -f2*sin(Theta1)];
Hs = [1 0 0].';
f = [ dX;...
      dTheta1;
      dTheta2;
      Hs*F-Gs-Cs*[dX dTheta1 dTheta2].'];
M = blkdiag(eye(3),Ds);
OCP.setf(f);
OCP.setM(M);
OCP.setDiscretizationMethod('Euler');

% Set the cost function L
Q = diag(OCP.p(1:6));
R = diag(OCP.p(7));
xRef  = [0;0;0;0;0;0];
uRef  = 0;
L     =   0.5*(OCP.x-xRef).'*Q*(OCP.x-xRef)...
        + 0.5*(OCP.u(1)-uRef).'*R*(OCP.u(1)-uRef)...
        +1e6*slack^2;
OCP.setL(L);

% Set the linear constraints G(u,x,p)>=0
G =[ OCP.u(1) + 10;...
    -OCP.u(1) + 10;...
     slack;...
     X + slack + 0.09;...
    -X + slack + 0.09];
OCP.setG(G);

% Generate necessary files
OCP.codeGen();
%% Configrate the solver using Class NMPCSolver

% Create a NMPCSolver object
nmpcSolver = NMPCSolver(OCP);

% Configurate the Hessian approximation method
nmpcSolver.setHessianApproximation('GaussNewton');

% Generate necessary files
nmpcSolver.codeGen();
%% Solve the very first OCP for a given initial state and given parameters using Class OCPSolver

% Set the initial state
x0   = [0;pi;pi;0;0;0];

% Set the parameter p
dim      = OCP.dim;
N        = OCP.N;
QDiagVal = [10;10;10;1;1;1];
RDiagVal = 0.1;
p        = zeros(dim.p,N);
p(1:7,:) = repmat([QDiagVal;RDiagVal],1,N);
p(1:7,end) = [100;100;100;10;10;10;0.1]; % terminal penlty

% init
solutionInitGuess.lambda = [randn(dim.lambda,1),zeros(dim.lambda,1)];
solutionInitGuess.mu     = [randn(dim.mu,1),randn(dim.mu,1)];
solutionInitGuess.u      = [uRef;1];
solutionInitGuess.x      = [x0,xRef];
solutionInitGuess.z      = ones(dim.z,N);
solutionInitGuess.LAMBDA = zeros(dim.lambda,dim.lambda,N);
solution = NMPC_SolveOffline(x0,p,solutionInitGuess,0.01,1000);

plot(solution.x([2 3],:).');
figure(2);
plot(solution.u(1,:).');
figure(3);
plot(solution.x(1,:).');

% Save to file
save GEN_initData.mat  dim x0 p N
global ParNMPCGlobalVariable
ParNMPCGlobalVariable.solutionInitGuess = solution;
%% Define the controlled plant using Class DynamicSystem

% M(u,x,p) \dot(x) = f(u,x,p)
% Create a DynamicSystem object
plant = DynamicSystem(1,6,0);

% Give names to x, u
[X,Theta1,Theta2,dX,dTheta1,dTheta2] = ...
    plant.setStateName({'X','Theta1','Theta2','dX','dTheta1','dTheta2'});
[F] = ...
    plant.setInputName({'F'});

% Set the dynamic function f
plant.setf(f); % same model
plant.setM(M); % same model

% Generate necessary files
plant.codeGen();

author: deng-haoyang

repo: ParNMPC


text:
function nnz = mm_nnz_set ( rep, symm, nrow, ncol )

%*****************************************************************************80
%
%% MM_NNZ_SET sets the value of NNZ for the ARRAY representation.
%
%  Discussion:
%
%    If the representation is not "ARRAY", then NNZ is returned as 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ( len = 10 ) REP, the Matrix Market 'representation'
%    indicator.  Possible values include:
%    'coordinate'   (for sparse data)
%    'array'        (for dense data)
%    'elemental'    (to be added)
%
%    Input, character ( len = 19 ) SYMM, the Matrix Market symmetry.
%    Possible values include:
%    'symmetric'
%    'hermitian'
%    'skew-symmetric'
%    'general'
%
%    Input, integer NROW, the number of rows in the matrix.
%
%    Input, integer NCOL, the number of columns in the matrix.
%
%    Output, integer NNZ, the number of nonzero entries required to store
%    the matrix.
%
  nnz = 0;

  if ( s_eqi ( rep, 'coordinate' ) )

  elseif ( s_eqi ( rep, 'array' ) )

    if ( s_eqi ( symm, 'general' ) )
      nnz = nrow * ncol;
    elseif ( s_eqi ( symm, 'symmetric' ) || s_eqi ( symm, 'hermitian' ) )
      nnz = floor ( ( nrow * ncol - nrow ) / 2 ) + nrow;
    elseif ( s_eqi ( symm, 'skew-symmetric' ) )
      nnz = floor ( ( nrow * ncol - nrow ) / 2 );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = optimiMinimize(objectiveGradient, params, options, varargin);

% OPTIMIMINIMIZE Wrapper for Carl Rasmussen's minimize function.
% FORMAT
% DESC is a nother wrapper for Carl Rasmussen's minimize function,
% but this time using functions which return the gradient and the
% function value together (such as gpObjectiveGradient).
% ARG objectiveGradient : function that returns the value (as a
% ARG params : the initial parameters to be optimised.
% scalar) and the gradients (as a row vector) of the function to be
% optimised.
% ARG options : a NETLAB style options vector.
% ARG P1, P2, P3 ... : further arguments to be passed to
% OBJECTIVEGRADIENT.
% RETURN params : row vector of optimised parameters.
%
% SEEALSO : minimize, cgcarl
%
% COPYRIGHT : Neil D. Lawrence, 2006

% OPTIMI

x = minimize(params', objectiveGradient, options(14), ...
             varargin{:})';

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Control of a water tank temperature
% Case of Proportional-Integral Control
% Author's Data: Housam BINOUS
% Department of Chemical Engineering
% National Institute of Applied Sciences and Technology
% Tunis, TUNISIA
% Email: binoushousam@yahoo.com

% function PI_control

function xdot=PI_control(t,x)

V = 100;
F = 10;
rho = 1;
cp = 4.19;
Q0 = 2500;
TRset = 80;
T0 = 20;
TauSens = 2;
Kp = 300;
TauQ = 1;
TauI=9;

QC= Q0 + Kp*(TRset - x(3))+Kp/TauI*x(4);

xdot(1) = 1/(V*rho*cp)*(F*rho*cp*(T0-x(1))+x(2));
xdot(2) = (QC-x(2))/TauQ;
xdot(3) = 1/TauSens*(x(1)-x(3));
xdot(4) = TRset - x(3);

xdot=xdot';

author: Sable

repo: mcbench-benchmarks


text:
function prob = sedumi2opti(prob)
%SEDUMI2OPTI  Converts a SEDUMI Problem (At, b, C, K) to OPTI Format
% prob = sedumi2opti(prob)

%   Copyright (C) 2013 Jonathan Currie (IPL)

if(~isfield(prob,'sdcone') || isempty(prob.sdcone) || ~isstruct(prob.sdcone))
    error('This function expects a prob.sdcone to contain a structure');
end
sdcone = prob.sdcone;
if((~isfield(sdcone,'A') && ~isfield(sdcone,'At')) || ~isfield(sdcone,'b') || (~isfield(sdcone,'C') && ~isfield(sdcone,'c')) || ~isfield(sdcone,'K'))
    error('This function expects prob.sdcone to contain the fields At (or A), b, C and K');
end
if(~isstruct(sdcone.K))
    error('SeDuMi prob.sdcone.K must be a structure!');
end
if(~isfield(sdcone.K,'s') && ~isfield(sdcone.K,'l'))
    error('SeDuMi prob.sdcone.K must contain the fields l or s (or both)');
end
if((isfield(sdcone.K,'scomplex') && ~isempty(sdcone.K.scomplex)) || (isfield(sdcone.K,'ycomplex') && ~isempty(sdcone.K.ycomplex)) || (isfield(sdcone.K,'xcomplex') && ~isempty(sdcone.K.xcomplex)))
    error('Complex variables are not supported in this interface - use SeDuMi directly to solve!');
end
if(isfield(sdcone.K,'f') && ~isempty(sdcone.K.f) && any(sdcone.K.f > 0))
    error('SeDuMi free variables are not supported in this interface - use SeDuMi directly to solve!');
end
if(isfield(sdcone.K,'q') && ~isempty(sdcone.K.q))
    error('SeDuMi SOCP constraints are not supported in this interface - use SeDuMi directly to solve!');
end
if(isfield(sdcone.K,'r') && ~isempty(sdcone.K.r))
    error('SeDuMi rotated cone constraints are not supported in this interface - use SeDuMi directly to solve!');
end

%Read SeDuMi problem, convert to OPTI standard primal form
% MINIMIZE   dot(f,x) SUCH THAT  SUM A_i*x_i - C >= 0.
top = 1; dims = 0;
haveSDP = false; haveLP = false;
K = sdcone.K;
if(isfield(sdcone,'C'))
    c = sdcone.C;
else
    c = sdcone.c;
end
if(isfield(sdcone,'A'))
    A = sdcone.A;
else
    A = sdcone.At;
end
b = full(sdcone.b);

%Check c, b is a vector
if(size(c,2) > 1 && size(c,1) > 1)
    error('SeDuMi c field must be a column vector!');
end
if(size(b,2) > 1 && size(b,1) > 1)
    error('SeDuMi b field must be a column vector!');
end

%Flip A based on b dimensions
if(size(A,2) ~= length(sdcone.b))
    if(size(A,1) == length(sdcone.b))
        A = A';
    else
        error('SeDuMi A (At) appears to  have incorrect dimensions.\nExpected one dimension of A to be the same as length as %s.','b');
    end
end
if(size(c,2) > size(c,1))
    c = c';
end
%Check length of C against A
if(length(c) ~= size(A,1))
    error('SeDuMi c does not have the same number of elements as rows in A (At)');
end
%Check dims entered in K
if(isfield(K,'s') && ~isempty(K.s) && K.s(1) > 0)
    haveSDP = true;
    dims = dims + sum(K.s.^2);
end
if(isfield(K,'l') && ~isempty(K.l) && K.l(1) > 0)
    if(length(K.l) > 1)
        error('The SeDuMi K.l field must contain a scalar, the number of linear constraints');
    end
    haveLP = true;
    dims = dims + K.l;
end
%Check dims against rows
if(dims ~= size(A,1))
    error('SeDuMi A (At) does not have correct dimensions as per the specifications in K. Expected %d rows.',dims);
end

%Extract Objective
prob.f = full(-b);

%Extract Linear Constraints
if(haveLP)
    n = K.l;
    if(isfield(prob,'A') && ~isempty(prob.A)) %concatenate
        prob.A = [prob.A; A(top:top+n-1,:)];
        prob.b = full([prob.b; prob.c(top:top+n-1,:)]);
    else
        prob.A = A(top:top+n-1,:);
        prob.b = full(c(top:top+n-1,:));
    end
    top = top+n;
end
%Extract Semidefinite Constraints
if (haveSDP)
    prob.sdcone = cell(length(K.s),1);
    for i = 1:length(K.s)
        n = K.s(i);
        prob.sdcone{i} = [-c(top:top+n^2-1,:) -A(top:top+n^2-1,:)];
        top = top+n^2;
    end
else
    prob.sdcone = [];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [I,B1,B2,B3] = in_element(V,F,P,varargin)
  % IN_ELEMENT test for each p in P whether it lies inside each f in F defined
  % over V.
  %
  % I = in_element(V,F,P)
  % [I,B1,B2,B3] = in_element(V,F,P,'ParameterName',ParameterValue,...)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by dim+1 list of element indices
  %   P  #P by dim list of query positions
  %   Optional:
  %     'Method' followed by one of the following {'knn'}:
  %       'brute-force' no acceleration O(#P * #F)
  %       'edge-walk' walk along edges ~O(#P * sqrt(#F)) Starting with a random
  %         barycenter step along the edges that intersect with the ray toward
  %         the query point. If a boundary is reached then search along all
  %         boundary edges and jump to farthest hit. **dim=2 only**
  %       'knn' use knnsearch to find closest element barycenters
  %       'spatial-hash' spatial hash on regular grid ~O(#P * sqrt(#F)) **dim=2
  %         only**
  %     'First'  only keep first match {false}
  %     'Quiet' suppress warnings {false}
  %     'Epsilon' epsilon used for determining inclusion {eps}
  % Outputs:
  %   I  #P by #F matrix of bools
  %   B1  #P by #F list of barycentric coordinates
  %   B2  #P by #F list of barycentric coordinates
  %   B3  #P by #F list of barycentric coordinates
  %
  % Example:
  %   P = bsxfun(@plus,min(V),bsxfun(@times,rand(100,2),max(V)-min(V)));
  %   [I,B1,B2,B3] = in_element(V,F,P);
  %   % Only keep first
  %   [mI,J] = max(I,[],2);
  %   I = sparse(1:size(I,1),J,mI,size(I,1),size(I,2));
  %   % Mask barycentric coordinates
  %   B1 = B1.*I;
  %   B2 = B2.*I;
  %   B3 = B3.*I;
  %   Q = B1*V(F(:,1),:) + B2*V(F(:,2),:) + B3*V(F(:,3),:);
  %   Q = Q(any(I,2),:);
  %   tsurf(F,V);
  %   hold on;
  %   plot(Q(:,1),Q(:,2),'or','LineWidth',6);
  %   plot(P(:,1),P(:,2),'ob','LineWidth',2);
  %   hold off;
  %

  function [I] = in_element_brute_force(V,F,P)
    dim = size(V,2);
    assert(dim+1 == size(F,2));

    % number of elements
    m = size(F,1);
    % number of query points
    np = size(P,1);

    switch dim
    case 3
      T = F;
      % tet face ares
      vol = abs(volume(V,T));
      allF = [ ...
        T(:,2) T(:,4) T(:,3); ...
        T(:,1) T(:,3) T(:,4); ...
        T(:,1) T(:,4) T(:,2); ...
        T(:,1) T(:,2) T(:,3); ...
        ];
      % List of tets for each face f of each tet t for each point p
      TP = cat(2, ...
        repmat(allF,[1 1 np]), ...
        permute(repmat(size(V,1)+(1:np),m*4,1),[1 3 2]));
      TP = reshape(permute(TP,[1 3 2]),m*4*np,dim+1);
      Pvol = abs(volume([V;P],TP));
      % Pvol(t,f,p) --> volume of tet t, face f with point p
      Pvol = reshape(Pvol,[size(T,1) dim+1 np]);
      % sumvol(p,t) --> sum of volumes of tets made with faces of t and p
      sumvol = permute(sum(Pvol,2),[3 1 2]);
      I = sparse(abs(bsxfun(@minus,sumvol,vol')) < sqrt(epsilon));
    case 2
      % triangle side lengths
      l = [ ...
        sqrt(sum((V(F(:,2),:)-V(F(:,3),:)).^2,2)) ...
        sqrt(sum((V(F(:,3),:)-V(F(:,1),:)).^2,2)) ...
        sqrt(sum((V(F(:,1),:)-V(F(:,2),:)).^2,2)) ...
        ];

      B = zeros([np m dim+1]);
      for ii = 1:(dim+1)
        jj = mod(ii+1,dim+1)+1;
        kk = mod(ii,dim+1)+1;
        ljj = pdist2(P,V(F(:,jj),:));
        lkk = pdist2(P,V(F(:,kk),:));

        % semiperimeters
        s = bsxfun(@plus,l(:,ii)',ljj + lkk)*0.5;
        % Heron's formula for area
        B(:,:,ii) = 2*sqrt(s.*(bsxfun(@minus,s,l(:,ii)').*(s-ljj).*(s-lkk)));
      end
      % sum of barycentric coordinates
      sumA = sum(B,3);
      % area of element
      dblA = doublearea(V,F);
      %% check whether sum is more than true are
      %I = ~bsxfun(@gt,sumA,dblA');
      I = sparse((bsxfun(@minus,sumA,dblA')) < sqrt(epsilon));
    case 1
      % To-do: this is actually being computed in a dense way
      I = sparse( ...
        ((P <= V(F(:,1))') & (P >= V(F(:,2))')) | ...
        ((P <= V(F(:,2))') & (P >= V(F(:,1))')));
    end
    %B1 = sparse(B(:,:,1));
    %B2 = sparse(B(:,:,2));
    %B3 = sparse(B(:,:,3));
  end

  function I = in_element_hash_helper(V,F,P)
    assert(size(F,2) == 3, ...
      'F must contain triangles for Method=''spatial-hash''');
    num_bins = ceil(sqrt(size(F,1)));
    bin_x = ceil(sqrt(num_bins));
    bin_y = ceil(num_bins/bin_x);
    num_bins = bin_x*bin_y;
    % spatial hash
    function VH = hash(V,MN,MX,bin_x,bin_y)
      [~,X] = histc(V(:,1),linspace(MN(1),MX(1),bin_x));
      [~,Y] = histc(V(:,2),linspace(MN(2),MX(2),bin_y));
      VH = sub2ind([bin_x bin_y],X,Y);
    end
    %% http://stackoverflow.com/a/5929567/148668
    %primes = [ 40960001, 59969537 45212177];
    %hash = @(X) mod( ...
    %  bitxor( int32(V(:,1)*primes(1)), int32(V(:,2)*primes(2))), ...
    %  num_bins);
    MN = min([V;P]);
    MX = max([V;P]);
    VH = hash(V,MN,MX,bin_x,bin_y);
    PH = hash(P,MN,MX,bin_x,bin_y);
    % This is wrong for triangles that span more hash cells than their vertices:
    % any time a triangle lands on the corner....
    FH = sparse(repmat(1:size(F,1),1,size(F,2))',VH(F(:)),1,size(F,1),num_bins)~=0;
    [~,FHI,FHJ] = find(FH);
    [FHJX,FHJY] = ind2sub([bin_x,bin_y],FHJ);
    % This assumes that #P >> #F
    I = sparse(size(P,1),size(F,1));
    for h = 1:num_bins
      %Vh = V(VH==h,:);
      IFh = FH(:,h);
      if any(IFh)
        IPh = PH==h;
        Ph = P(IPh,:);
        if any(IPh)
          Fh = F(IFh,:);
          Ih = in_element_brute_force(V,Fh,Ph);
          I(IPh,IFh) = Ih;
        end
      end
    end
  end

  % default values
  method = [];
  first = false;
  quiet = false;
  epsilon = eps;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'Method','First','Quiet','Epsilon'}, ...
    {'method','first','quiet','epsilon'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this
      % workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  if isempty(method)
    switch size(V,2)
    case 1
      method = 'brute-force';
    otherwise
      method = 'knn';
    end
  end

  switch method
  case 'brute-force'
    I = in_element_brute_force(V,F,P);
  case 'spatial-hash'
    % Try 45?? spatial grid, too (reduce corner cases)
    switch size(V,2)
    case 2
      R = [cos(pi/4) -sin(pi/4);sin(pi/4) cos(pi/4)];
    case 3
      R = [cos(pi/4) -sin(pi/4) 0 ;sin(pi/4) cos(pi/4) 0; 0 0 1];
    end
    I = in_element_hash_helper(V,F,P) | in_element_hash_helper(V*R,F,P*R);

    % Find any obviously incorrect values: not inside but winding number says
    % inside. Could still missing something if mesh overlaps itself.
    NI = ~any(I,2);
    switch size(F,2)
    case 3
      O = outline(F);
    case 4
      O = boundary_faces(F);
    end
    WI = abs(winding_number(V,O,P(NI,:))/(2*pi))>0.5;
    WI = sparse(find(NI),1,WI,size(P,1),1)~=0;
    % redo any that currently are not in any but winding number says are inside
    RI = NI & WI;
    I(RI,:) = in_element_brute_force(V,F,P(RI,:));

  case 'edge-walk'

    assert(size(F,2) == 3,'F must contain triangles for Method=''edge-walk''');
    % List of all "half"-edges: 3*#F by 2
    allE = [F(:,[2 3]); F(:,[3 1]); F(:,[1 2])];
    % Sort each row
    sortallE = sort(allE,2);
    % IC(i) tells us where to find sortallE(i,:) in uE:
    % so that sortallE(i,:) = uE(IC(i),:)
    [uE,~,IC] = unique(sortallE,'rows');
    % uE2F(e,f) = i means face f's ith edge is unique edge e
    uE2F = sparse(IC(:),repmat(1:size(F,1),1,3)',reshape(repmat(1:3,size(F,1),1),[],1));
    % uE2F(e,f) = 1 means face f is adjacent to unique edge e
    uE2F1 = sparse(IC(:),repmat(1:size(F,1),1,3)',1);
    % Face-face Adjacency matrix
    A = uE2F1'*uE2F;
    % A(f,g) = i means face f's ith edge is shared with g
    A = A-diag(diag(A));

    I = sparse(size(P,1),size(F,1));
    B1 = sparse(size(I,1),size(I,2));
    B2 = sparse(size(I,1),size(I,2));
    B3 = sparse(size(I,1),size(I,2));

    [~,is_b] = on_boundary(F);
    EF = repmat(1:size(F,1),1,3)';
    EFI = reshape(repmat(1:3,size(F,1),1),[],1);
    O = allE(is_b(:),:);
    OF = EF(is_b(:));
    OFI = EFI(is_b(:));

    % initial closest vertex
    BC = barycenter(V,F);
    % Centroid
    f_init = snap_points(mean(V),BC);
    %% Random initial guess
    %f_init = ceil(rand(1,1)*size(F,1));
    for p = 1:size(P,1)
      % current closest face, barycenter
      f = f_init;
      q = BC(f_init,:);
      % incoming edge
      e_in = [];
      while true
        % current point
        % edges to test
        E = ... %mod(bsxfun(@plus,setdiff([1;2;3],e_in),-1+(1:2)),3)+1;
          [2 3;3 1;1 2];
        out = ...
          lineSegmentIntersect([q P(p,:)],[V(F(f,E(:,1)),:) V(F(f,E(:,2)),:)]);
        out.intAdjacencyMatrix(e_in) = false;
        e_out = find(out.intAdjacencyMatrix);
        if isempty(e_out)
          % no hits so we're in the element
          I(p,f) = 1;
          B = barycentric_coordinates( ...
            P(p,:),V(F(f,1),:),V(F(f,2),:),V(F(f,3),:));
          B1(p,f) = B(1); B2(p,f) = B(2); B3(p,f) = B(3);
        else
          f_prev = f;
          if numel(e_out) > 1
            % Take farthests
            [sd,si] = sort(out.intNormalizedDistance1To2(e_out),'descend');
            e_out = e_out(si(1));
            % TODO: recurse on all in f
          end
          f = find(A(:,f_prev)==e_out);
          % Todo if
          if isempty(f)
            % no neighbors so we're shooting outside.
            out = ...
              lineSegmentIntersect([q P(p,:)],[V(O(:,1),:) V(O(:,2),:)]);
            hits = find(out.intAdjacencyMatrix);
            [sd,si] = sort(out.intNormalizedDistance1To2(hits),'descend');
            f = OF(hits(si(1)));
            e_in = OFI(hits(si(1)));
            if f==f_prev
              error('Point is outside');
            end
          else
            if numel(f) > 1
              f = f(1);
              if ~quiet
                warning('Ignoring non-manifold edge: might miss multiply inside.');
              end
              % TODO: recurse on all in f
            end
            e_in = A(f_prev,f);
          end
        end
        %tsurf(F,V);
        %hold on;
        %tsurf(F(f,:),V,'CData',-1);
        %tsurf(F(f_prev,:),V,'CData',1);
        %plot_edges([q;P(p,:)],[1 2],'r');
        %plot_edges(V,[F(f_prev,E(:,1));F(f_prev,E(:,2))]','b');
        %plot(P(p,1),P(p,2),'*');
        %plot_edges(V,[F(f_prev,E(e_out,1));F(f_prev,E(e_out,2))]', ...
        %  'y','LineWidth',2);
        %hold off;
        %drawnow;
        %input('');
      end
    end

    return
  case 'knn'
    % assumes samples are all inside exactly one element
    BC = barycenter(V,F);
    I = sparse(size(P,1),size(F,1));
    B1 = sparse(size(P,1),size(F,1));
    B2 = sparse(size(P,1),size(F,1));
    B3 = sparse(size(P,1),size(F,1));
    B4 = sparse(size(P,1),size(F,1));
    % indices of points we haven't found yet
    IP = 1:size(P,1);

    prev_k = 0;
    k = 2;
    while true
      K = knnsearch(BC,P(IP,:),'K',k);
      K = K(:,prev_k+1:end);
      for ki = 1:size(K,2)
        switch size(F,2)
        case 3
          B = abs(barycentric_coordinates( ...
            P(IP,:), ...
            V(F(K(:,ki),1),:), ...
            V(F(K(:,ki),2),:), ...
            V(F(K(:,ki),3),:)));
        case 4
          B = abs(barycentric_coordinates( ...
            P(IP,:), ...
            V(F(K(:,ki),1),:), ...
            V(F(K(:,ki),2),:), ...
            V(F(K(:,ki),3),:), ...
            V(F(K(:,ki),4),:)));
        end
        found = abs(sum(B,2)-1)<sqrt(epsilon);
        I(sub2ind(size(I),IP,K(:,ki)')) = found;
        % DOESN'T REALLY PAY OFF TO COMPUTE THESE HERE
        %B1(sub2ind(size(I),IP,K(:,ki)')) = found.*B(:,1);
        %B2(sub2ind(size(I),IP,K(:,ki)')) = found.*B(:,2);
        %B3(sub2ind(size(I),IP,K(:,ki)')) = found.*B(:,3);
        %if size(F,2) == 4
        %  B4(sub2ind(size(I),IP,K(:,ki)')) = found.*B(:,4);
        %end
        % Peel off found
        IP = IP(~found);
        if isempty(IP)
          break;
        end
        K = K(~found,:);
      end

      %for p = 1:numel(IP)
      %  Kp = K(p,:);
      %  U = V(F(Kp,:),:);
      %  FKp = reshape((1:size(U,1))',size(U,1)/size(F,2),size(F,2));
      %  I(IP(p),Kp) = in_element_brute_force(U,FKp,P(IP(p),:));
      %end
      %IP = find(~any(I,2));

      if isempty(IP)
        break;
      end
      prev_k = k;
      k = min(prev_k*2,size(BC,1));
      if k == size(BC,1)
        if ~quiet
          warning('Some points not found');
        end
        break;
      end
    end

    %return;

  end

  if first
    % Only keep first
    [mI,J] = max(I,[],2);
    I = sparse(1:size(I,1),J,mI,size(I,1),size(I,2));
  end

  % Compute barycentric coordinates
  if nargout > 1
    B1 = sparse(size(I,1),size(I,2));
    B2 = sparse(size(I,1),size(I,2));
    B3 = sparse(size(I,1),size(I,2));
    B4 = sparse(size(I,1),size(I,2));
    work_I = I;
    while true
      % Peel off a layer
      [mIf,Jf] = max(work_I,[],2);
      If = find(mIf);
      Jf = Jf(If);
      mIf = mIf(If);
      if isempty(If)
        break
      end
      work_I = work_I - sparse(If,Jf,mIf,size(work_I,1),size(work_I,2));;
      Pf = P(If,:);
      Ff = F(Jf,:);
      switch size(Ff,2)
      case 4
        Bf = barycentric_coordinates(Pf, ...
          V(Ff(:,1),:),V(Ff(:,2),:),V(Ff(:,3),:),V(Ff(:,4),:));
      case 3
        Bf = barycentric_coordinates(Pf,V(Ff(:,1),:),V(Ff(:,2),:),V(Ff(:,3),:));
      case 2
        Bf = barycentric_coordinates(Pf,V(Ff(:,1),:),V(Ff(:,2),:));
      end
      B1(sub2ind(size(B1),If,Jf)) = Bf(:,1);
      B2(sub2ind(size(B2),If,Jf)) = Bf(:,2);
      if size(Bf,2) >= 3
        B3(sub2ind(size(B3),If,Jf)) = Bf(:,3);
      end
      if size(Bf,2) >= 4
        B4(sub2ind(size(B3),If,Jf)) = Bf(:,4);
      end
    end
  end

end


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

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

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

for c = 1:m
	one_example = X(c,:);
	[max_value, max_index] = max(all_theta * one_example');
	p(c,1) = max_index;

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

end

author: sfvsfv

repo: Mathematical-modeling


text:
function sessionGUI_close;
% Close a Session GUI, clearing the relevant global variables along the
% way.
%
% sessionGUI_close;
%
% ras, 07/06
mrGlobals2;

% close any subsidiary viewers
for v = GUI.viewers
    mrViewClose(v);
end

% close the main figure
delete(GUI.fig);

INPLANE{1} = [];
VOLUME{1} = [];
if isempty(cellfind(INPLANE)), INPLANE = {}; end
if isempty(cellfind(VOLUME)), VOLUME = {}; end

% Check if there are no more open views: if not, clean the workspace
if isempty(cellfind(INPLANE)) & isempty(cellfind(VOLUME))
    mrvCleanWorkspace;
end

clear global GUI;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function y = entr( x )

%ENTR   Scalar entropy.
%   ENTR(X) returns an array of the same size as X with the unnormalized
%   entropy function applied to each element:
%                { -X.*LOG(X) if X > 0,
%      ENTR(X) = { 0          if X == 0,
%                { -Inf       otherwise.
%   If X is a vector representing a discrete probability distribution, then
%   SUM(ENTR(X)) returns its entropy.
%
%   Disciplined convex programming information:
%       ENTR(X) is concave and nonmonotonic in X. Thus when used in CVX
%       expressions, X must be real and affine. Its use will effectively
%       constrain X to be nonnegative: there is no need to add an
%       additional X >= 0 to your model in order to enforce this.

narginchk(1,1);
cvx_expert_check( 'entr', x );
y = -rel_entr( x, 1 );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% calc_frap - makes the actual analysis of the FRAP data
%   [k,D1k,D1,D2,gamma2,gamma0,v_conv,a_conv]=...
%   calc_frap(file_in,dx,Ipre,Idark,istart,t,trc,Rp_max,nr) where:
%
%   k =  spatial frequencies [um^-1]
%   D1k = D1(k2), the diffusion coefficient determined for each k
%       separately [um^2s^-1]
%   D1 and D2 = effective values of the diffusion coefficients D1 and
%       D2 [um^2s^-1]
%   gamma2 = the intensity fraction of component 2
%   gamma0 = the intensity fraction of immobile molecules
%   v_conv = the absolute velocity of the tracked center of mass [um/s]
%   a_conv = the angle the tracked center of mass moves in, with 0 degrees
%       corresponding to moving horizontally from the left to the right
%       [degrees]
%
%   file_in = the filename + directory of the image stack
%   dx = the pixel size (square pixels assumed)
%   Ipre = the pre-bleach intensity
%   Idark = the dark count intensity
%   istart = the first post-bleach frame
%   t = the times for each frame
%   trc = parameter that if equal to 'y' leads to tracking of the center of
%       mass for each frame (otherwise the center of mass is determined from
%       the first post-bleach frame)
%   Rp_max = the maximum radial value used in the analysis [pixels]
%   nr = a vector for all frames in the image stack. If nr==1 then this
%       image will be part of the analysis.

function [k,D1k,D1,D2,gamma2,gamma0,v_conv,a_conv]=...
    calc_frap(file_in,dx,Ipre,Idark,istart,t,trc,Rp_max,nr)

% closes the figures
close(figure(1))
close(figure(2))

% Meshpoints
x=(1:1:size(Ipre,2));
y=(size(Ipre,1):-1:1);
[Xx,Yy]=meshgrid(x,y);

% Calculates the centre of mass from the first frame after bleaching
I=double(imread(file_in,istart))-Idark;
I2=I./Ipre;

fig=figure(1);
set(fig,'Color','w');
I2_min=max(0.5,min(I2(:)));
imshow(I2,[I2_min,min([1+(1-I2_min)*0.1,max(I2(:))])]);
ax1 = gca;
set(ax1,'FontSize',8)
axis square
h=title({'Define a polygon around the bleached spot (close with the right mouse button).';...
    'Exit by double clicking the inner part of the polygon.'},...
    'Color','k','BackgroundColor','w');
hp=get(h);
pos=hp.Position;
pos(2)=pos(2)+50;
set(h,'Position',pos);
BW=roipoly; % the region in which the bleached spot is situated

% Calculates the centre of mass (x_cm,y_cm)
x_cm=sum(sum(Xx.*(1-I2).*BW))/sum(sum((1-I2).*BW));
y_cm=sum(sum(Yy.*(1-I2).*BW))/sum(sum((1-I2).*BW));
Rr=sqrt((Xx-x_cm).^2+(Yy-y_cm).^2);
wid=max(Rr(BW==1)); % max r-value to use to calculate the center of mass

figure(1)
hold on
plot3([x_cm-10,x_cm+10],[max(Yy(:))-y_cm+1,max(Yy(:))-y_cm+1],max(I2(:))*ones(1,2),...
    [x_cm,x_cm],[max(Yy(:))-y_cm+1-10,max(Yy(:))-y_cm+1+10],max(I2(:))*ones(1,2),'LineWidth',2,'Color',[1 1 1])
hold off
pause(0.2)

% Defines the options for nonlinear fitting
options=optimset('Display','off');

Nrp=512;    % the length of r

% Determines the radial distance r and the radius of the field of view R
[test,r,R,Rp_max]=I_radial(Nrp,x_cm,y_cm,Xx,Yy,Xx,Rp_max);
r=r*dx;
R=R*dx;

% Initiates variables
wf=zeros(sum(nr==1),1);
x_cm2=zeros(sum(nr==1),1);
y_cm2=zeros(sum(nr==1),1);
Ixx=zeros(size(x_cm2));
Iyy=zeros(size(x_cm2));
Ixy=zeros(size(x_cm2));
Af=zeros(sum(nr==1),1);
Isum=zeros(sum(nr==1),1);
Ir=zeros(sum(nr==1),Nrp);
g=zeros(sum(nr==1),Nrp);
i=0;

t=t(nr==1)';    % updates the time vector

for n_t=1:length(nr)
    % n_t==0 if a frame is omitted from the analysis
    if nr(n_t)==1
        i=i+1;

        I=double(imread(file_in,istart+n_t-1))-Idark;   % Reads the image
        I2=I./Ipre; % The relative intensity

        if i==1
            x_cm2(i)=x_cm;
            y_cm2(i)=y_cm;

            I20=I2;
        else
            x_cm2(i)=x_cm2(i-1);
            y_cm2(i)=y_cm2(i-1);
        end

        if i>2
            x_cm2(i)=(x_cm2(i-1)-x_cm2(i-2))/(t(i-1)-t(i-2))*(t(i)-t(i-1))+x_cm2(i-1);
            y_cm2(i)=(y_cm2(i-1)-y_cm2(i-2))/(t(i-1)-t(i-2))*(t(i)-t(i-1))+y_cm2(i-1);
        end

%         % Take away the immobile fraction
%         if trc=='y'
%             gamma_g=0.020;   % immobile fraction
%             I2=I2+(I2_inf-I20)*gamma_g;
%         end

        % Determines I2(inf) using values at the edge of the image
        nx=round(length(x)*0.1);
        ny=round(length(y)*0.1);
        I2s=sum(sum(I2(1:ny,:)))+sum(sum(I2(end:end-ny,:)))+...
            sum(sum(I2(:,1:nx)))+sum(sum(I2(:,end:end-nx)));
        I2n=sum(sum(~isnan(I2(1:ny,:))))+sum(sum(~isnan(I2(end:end-ny,:))))+...
            sum(sum(~isnan(I2(:,1:nx))))+sum(sum(~isnan(I2(:,end:end-nx))));
        I2_inf=I2s/I2n;
        beta=I2_inf;

        % If trc=='y' then the program tracks (and moves) the center of mass
        % for each frame
        Iin=(I2_inf-I2);
        if trc=='y'
            for j=1:3
                Rr2=(Xx-x_cm2(i)).^2+(Yy-y_cm2(i)).^2;
                BW=Rr2<(wid^2);
                cm_w=Iin(BW==1);
                Xx_w=Xx(BW==1);
                Yy_w=Yy(BW==1);

                % Calculates the centre of mass (x_cm,y_cm)
                x_cm2(i)=sum(Xx_w.*cm_w)/sum(cm_w);
                y_cm2(i)=sum(Yy_w.*cm_w)/sum(cm_w);
            end
        else
            Rr2=(Xx-x_cm2(i)).^2+(Yy-y_cm2(i)).^2;
            BW=Rr2<(wid^2);
            cm_w=Iin(BW==1);
            Xx_w=Xx(BW==1);
            Yy_w=Yy(BW==1);
            x_cm2(i)=x_cm;
            y_cm2(i)=y_cm;
        end

        % Calculates the moments (Ixx and Iyy) and product (Ixy) of inertia for the bleached spot
        Ixx(i)=sum((Yy_w-y_cm2(i)).^2.*cm_w);
        Iyy(i)=sum((Xx_w-x_cm2(i)).^2.*cm_w);
        Ixy(i)=sum((Xx_w-x_cm2(i)).*(Yy_w-y_cm2(i)).*cm_w);

        % Performs the radial averaging of I2 around the centre of mass
        Ir(i,:)=I_radial(Nrp,x_cm2(i),y_cm2(i),Xx,Yy,I2,Rp_max);

        % Determines Ir(R,0) using values where r>0.9R
        if i==1
            Ir0=mean(Ir(1,r>(max(r)*0.9)));
        end

        % Determines the integral of Ir at r<R
        Isum(i)=2*pi*Ir(i,1)*r(1)*r(1);
        for j=2:length(r)
            Isum(i)=Isum(i)+2*pi*(Ir(i,j)*r(j)+Ir(i,j-1)*r(j-1))/2*(r(j)-r(j-1));
        end

        % Determines starting parameters used to fit the tail of Ir
        y=Ir(i,:)/Ir0;
        x=r;
        pos=1;
        y0=y(pos);
        while (1-y0)>(1-y(1))*exp(-1) && pos<length(y)
            pos=pos+1;
            y0=y(pos);
        end
        w_g0=sqrt(x(pos)^2-x(1)^2);

        y=Ir(i,r>max(r)*1/2)/Ir0;
        x=r(r>max(r)*1/2);
        pos=1;
        y0=y(pos);
        while (1-y0)>(1-y(1))*exp(-1) && pos<length(y)
            pos=pos+1;
            y0=y(pos);
        end
        w_g=max(sqrt(x(pos)^2-x(1)^2),w_g0);

        % Makes a curve fit to the tail of Ir (for r>R/2) for each frame
        p0=[min(1-Ir(i,1)/Ir0,abs(1-y(1))*exp(x(1)^2/w_g^2)),w_g];
        p0_2=p0;
        if i>1
            p0_2=p;
        end

        [p1,sse1]=lsqcurvefit(@(p,r) fkn_gauss(p,r,Isum(i),Isum(1),Ir0),p0,x,y*Ir0,...
            [0,0],[Inf,Inf],options);

        [p2,sse2]=lsqcurvefit(@(p,r) fkn_gauss(p,r,Isum(i),Isum(1),Ir0),p0_2,x,y*Ir0,...
            [0,0],[Inf,Inf],options);

        if sse1<sse2
            p=p1;
        else
            p=p2;
        end

        wf(i)=p(2); % w(t)
        Af(i)=p(1); % A(t)
        [y,beta]=fkn_gauss(p,r,Isum(i),Isum(1),Ir0);    % y = Ir|r>R

        % Compensates for temporal variations using beta and for
        % a net influx of molecules with y
        g(i,:)=(y-Ir(i,:))/beta;

        if mod(i,10)==0 || i==1 || i==sum(nr==1)
            figure(2)
            ax1 = gca;
            set(ax1,'FontSize',12)
            plot(r,Ir(i,:)/beta,'b.',r(r>(max(r)*1/2)),y(r>(max(r)*1/2))/beta,'r-','LineWidth',2)
            xlabel('r [\mum]','FontSize',12)
            ylabel('Intensity [a.u.]','FontSize',12)
            title(['Frame ',num2str(i),' out of ',num2str(sum(nr==1))],'FontSize',12)
            xlim([min(r),max(r)])
            pause(0.1)
        end

    end
end

% Calculates the Hankel transform and performs the curve fits yielding D1,
% D2, gamma2 and gamma0
[k,D1k,D1,D2,gamma2,gamma0,v_conv,a_conv]=...
    hankel_diff(r,t,R,g,wf,Af,Nrp,trc,x_cm2,y_cm2,dx,Ixx,Iyy,Ixy);

end


text:
function [Q,w,S] = quadrature_points(V,F,k)
  % QUADRATURE_POINTS Sample quadrature points on a triangle or tetrahedron.
  %
  % Q = quadrature_points(V,F)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by simplex-size list of indices into V
  %   k  number of quadrature points {1,3,6}
  % Outputs:
  %   Q  #F by dim by k list of quadrature points
  %   w  #k list of weights
  %   S  #F*k by #V sparse matrix so that Q ~= S*V (~= means after
  %     reordering)
  %
  % Example:
  %   Q = quadrature_points(V,F);
  %   QQ = reshape(permute(Q,[1 3 2]),[],3);
  %   QQQ = permute(reshape(QQ,permute(size(Q),[1 3 2])),[1 3 2]);
  %   max(abs(Q-QQQ))
  %
  % See also: barycenter
  %

  % Compute element quadrature points
  % "A SET OF SYMMETRIC QUADRATURE RULES ON TRIANGLES AND TETRAHEDRA", Table 2.2
  dim = size(V,2);
  ss = size(F,2);
  switch ss
  case 3
    switch k
    case 1
      lambda = [1/3 1/3 1/3];
    case 3
      %a = 0.5;
      a = 0.15;
      lambda = unique(perms([a a 1-2*a]),'rows');
    case 6
      a = 0.05;
      b = 0.1;
      lambda = unique(perms([a b 1-a-b]),'rows');
    otherwise
    error(sprintf('Unsupported number of quadratures %d',k));
    end
    % uniform weights
    w = ones(size(lambda,1),1)/size(lambda,1);
    QQ = cat(4,V(F(:,1),:),V(F(:,2),:),V(F(:,3),:));
  otherwise
    error(sprintf('Unsupported simplex-size %d',ss));
  end

  % #F by #lambda matrix of indices of all of the new points
  I = ((1:size(F,1))') + (0:size(lambda,1)-1)*size(F,1);
  S = sparse( ...
      repmat(I,size(F,2),1), ...
      repmat(F(:),1,size(lambda,1)), ...
      reshape(permute(repmat(lambda,[1 1 size(F,1)]),[3 2 1]),numel(F),size(lambda,1)), ...
      size(F,1)*size(lambda,1), ...
      size(V,1));
  % Q2 = S*V;

  lambda = permute(lambda,[3 4 1 2]);
  Q = sum(bsxfun(@times,QQ,lambda),4);

  % The Q ordering is silly. S*V ~= Q after gnarly reordering
  assert(max(max(max(abs(Q-permute(reshape(S*V,[size(Q,1) size(Q,3) size(Q,2)]),[1 3 2])))))<1e-12)

end

author: alecjacobson

repo: gptoolbox


text:
function n = p05_n ( )

%*****************************************************************************80
%
%% P05_N returns the number of variables N for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N, the number of variables.
%
  n = 10;

  return
end


text:
function length = polyloop_length_nd ( dim_num, nk, pk )

%*****************************************************************************80
%
%% POLYLOOP_LENGTH_ND computes the length of a polyloop in ND.
%
%  Discussion:
%
%    A polyloop of order NK is the geometric structure consisting of
%    the NK line segments that lie between successive elements of a list
%    of NK points, with the last point joined to the first.
%
%    Warning: I just made up the word "polyloop".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NK, the number of points defining the polyloop.
%
%    Input, real PK(DIM_NUM,NK), the points defining the polyloop.
%
%    Output, real LENGTH, the length of the polyloop.
%
  length = 0.0;

  for i = 2 : nk + 1

    if ( i <= nk )
      j = i;
    else
      j = 1;
    end if

    length = length + sqrt ( sum ( ( pk(1:dim_num,j) - pk(1:dim_num,i-1) ).^2 ) );

  end

  return
end


text:
function [data_mean,data_diff,md,sd] = bland_altman(data1,data2)
% Function to generate Bland Altman plots. Barry Greene, September 2008
% Bland, J.M., Altman, D.G. 'Statistical methods for assessing agreement ...
% between two methods of clinical measurement'(1986) Lancet, 1 (8476), pp. 307-310.
% Inputs: data1: Data from first instrument
%         data2: Data from second instument
% Produces Bland Altman plot with mean difference and mean difference +/-
% 2*SD difference lines.

[m,n] = size(data1);
if(n>m)
    data1 = data1';
end

if(size(data1)~=size(data2))
    error('Data matrices must be the same size')
end

data_mean = mean([data1,data2],2);  % Mean of values from each instrument
data_diff = data1 - data2;              % Difference between data from each instrument
md = mean(data_diff);               % Mean of difference between instruments
sd = std(data_diff);                % Std dev of difference between instruments

% figure;
plot(data_mean,data_diff,'ok','MarkerSize',3,'LineWidth',1);   % Bland Altman plot
hold on; x = [min(data_mean),max(data_mean)];
y = md*ones(1,2); plot(x,y,'-k'); % Mean difference line
y = md+2*sd*ones(1,2);
plot(x,y,'--k'); % Mean plus 2*SD line
% text(x(2),y(2),'+2 SD');
y = md-2*sd*ones(1,2);
plot(x,y,'--k'); % Mean minus 2*SD line

% text(x(2),y(2),'-2 SD');
% grid on
% title('Bland Altman plot','FontSize',9)
% xlabel('Mean of two measures','FontSize',8)
% ylabel('Difference between two measures','FontSize',8)

author: zhouyuanzxcv

repo: Hyperspectral


text:
function varargout = sv_gridMc(varargin)
% SV_GRIDMc Application M-file for sv_gridMv.fig
%    FIG = SV_GRID launch sv_grid GUI.
%    SV_GRIDMc('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.0 22-Aug-2002 10:30:09

if nargin == 1  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Use system color scheme for figure:
	set(fig,'Color',get(groot,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it.
	handles = guihandles(fig);
	guidata(fig, handles);

if varargin{1} == 1
    set(handles.lblSpacingDepth, 'String', 'Latitude spacing [deg]:');
    set(handles.txtSpacingDepth, 'String', '0.1');
    set(handles.lblSpacingHorizontal, 'String', 'Longitude spacing [deg]:');
    set(handles.txtSpacingHorizontal, 'String', '0.1');
else
    set(handles.lblSpacingDepth, 'String', 'Depth spacing [km]:');
    set(handles.txtSpacingDepth, 'String', '1');
    set(handles.lblSpacingHorizontal, 'String', 'Horizontal spacing [km]:');
    set(handles.txtSpacingHorizontal, 'String', '1');
end
if get(handles.radNumber, 'Value') == 1
    set(handles.txtNumber, 'Enable', 'on');
    set(handles.radRadius, 'Value', [0]);
    set(handles.txtRadius, 'Enable', 'off');
else
    set(handles.txtNumber, 'Enable', 'off');
    set(handles.radRadius, 'Value', [1]);
    set(handles.txtRadius, 'Enable', 'on');
end

% if get(handles.radSplitTime, 'Value') == 1
%     set(handles.txtSplitTime, 'Enable', 'on');
%     set(handles.txtTimePeriod, 'Enable', 'on');
% else
%     set(handles.radSplitTime, 'Value', [1]);
%     set(handles.txtSplitTime, 'Enable', 'on');
%     set(handles.txtTimePeriod, 'Enable', 'on');
% end
	% Wait for callbacks to run and window to be dismissed:
	uiwait(fig);

	if nargout > 0
		varargout{1} = fig;
	end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end

%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and
%| sets objects' callback properties to call them through the FEVAL
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

% --------------------------------------------------------------------
% Radiobutton: Number of events
function varargout = radNumber_Callback(h, eventdata, handles, varargin)

if get(handles.radNumber, 'Value') == 1
  set(handles.txtNumber, 'Enable', 'on');
  set(handles.txtMaxRadius, 'Enable', 'on');
  set(handles.radRadius, 'Value', [0]);
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radRectangle, 'Value', [0]);
  set(handles.txtSizeRectHorizontal, 'Enable', 'off');
  set(handles.txtSizeRectDepth, 'Enable', 'off');
else
  set(handles.txtNumber, 'Enable', 'off');
  set(handles.txtMaxRadius, 'Enable', 'off');
  set(handles.radRadius, 'Value', [1]);
  set(handles.txtRadius, 'Enable', 'on');
  set(handles.radRectangle, 'Value', [0]);
  set(handles.txtSizeRectHorizontal, 'Enable', 'off');
  set(handles.txtSizeRectDepth, 'Enable', 'off');
end

% --------------------------------------------------------------------
% Radiobutton: Constant radius
function varargout = radRadius_Callback(h, eventdata, handles, varargin)

if get(handles.radRadius, 'Value') == 1
  set(handles.txtRadius, 'Enable', 'on');
  set(handles.radNumber, 'Value', [0]);
  set(handles.txtMaxRadius, 'Enable', 'off');
  set(handles.txtNumber, 'Enable', 'off');
  set(handles.radRectangle, 'Value', [0]);
  set(handles.txtSizeRectHorizontal, 'Enable', 'off');
  set(handles.txtSizeRectDepth, 'Enable', 'off');
else
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.txtMaxRadius, 'Enable', 'on');
  set(handles.radNumber, 'Value', [1]);
  set(handles.txtNumber, 'Enable', 'on');
  set(handles.radRectangle, 'Value', [0]);
  set(handles.txtSizeRectHorizontal, 'Enable', 'off');
  set(handles.txtSizeRectDepth, 'Enable', 'off');
end

% --------------------------------------------------------------------
% Radiobutton: Rectangle size
function varargout = radRectangle_Callback(h, eventdata, handles, varargin)

if get(handles.radRectangle, 'Value') == 1
  set(handles.txtSizeRectHorizontal, 'Enable', 'on');
  set(handles.txtSizeRectDepth, 'Enable', 'on');
  set(handles.radRadius, 'Value', [0]);
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radNumber, 'Value', [0]);
  set(handles.txtNumber, 'Enable', 'off');
  set(handles.txtMaxRadius, 'Enable', 'off');
else
  set(handles.radRadius, 'Value', [0]);
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radNumber, 'Value', [1]);
  set(handles.txtNumber, 'Enable', 'on');
  set(handles.txtMaxRadius, 'Enable', 'on');
end

% --------------------------------------------------------------------
% Radiobutton: Split Time
function varargout = radSplitTime_Callback(h, eventdata, handles, varargin)

if get(handles.radSplitTime, 'Value') == 1
   set(handles.txtSplitTime, 'Enable', 'on');
   set(handles.txtTimePeriod, 'Enable', 'on');
else
  set(handles.txtTimePeriod, 'Enable', 'on');
  set(handles.txtSplitTime, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkTimePeriod_Callback(h, eventdata, handles, varargin)

if get(handles.chkTimeperiod, 'Value') == 1
  set(handles.txtTimeperiod, 'Enable', 'on');
else
  set(handles.txtTimeperiod, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkSaveParameter_Callback(h, eventdata, handles, varargin)

if get(handles.chkSaveParameter, 'Value') == 1
  set(handles.btnSaveParameter, 'Enable', 'on');
else
  set(handles.btnSaveParameter, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = btnSaveParameter_Callback(h, eventdata, handles, varargin)

[newfile, newpath] = uiputfile('*.mat','Choose parameter file');
% Cancel pressed?
if isequal(newfile, 0)  ||  isequal(newpath, 0)
  return;
end
% Everything ok?
newfile = [newpath newfile];
if length(newfile) > 1
  set(handles.lblSaveParameter, 'String', newfile);
end

% --------------------------------------------------------------------
function varargout = btnOK_Callback(h, eventdata, handles, varargin)

handles.answer = 1;
guidata(h, handles);
uiresume(handles.figGrid);

% --------------------------------------------------------------------
function varargout = btnCancel_Callback(h, eventdata, handles, varargin)

handles.answer = 0;
guidata(h, handles);
uiresume(handles.figGrid);

author: CelsoReyes


text:
% go to every subject in the BA44-45 list
% load each fiber group in the BA44-45\dissections dir
% transpose  fg.fibers
% save as .._fixed.mat

baseDir = 'U:\data\reading_longitude\dti_adults\';
subDirs = {'ab050307','as050307','aw040809','bw040806','da050311','gm050308',...
        'jl040902','ka040923','mbs040503', 'me050126', 'mz040828',...
        'pp050208', 'rd040630','sn040831','sp050303'};
fiberDir = 'fibers\BA44-45\dissections\';

for ii = 1:length(subDirs)
    workDir = fullfile(baseDir,subDirs{ii},fiberDir);
    cd(workDir);
    d = dir('?BA*.mat');
    fileNames = {d.name};
    for jj = 1:length(fileNames)
        curFileName = fileNames{jj};
        load(curFileName);
        fg.fibers = fg.fibers';
        newFileName = [curFileName(1:end-4) '_fixed.mat'];
        save(newFileName,'fg','versionNum','coordinateSpace');
        clear fg versionNum coordinateSpace
    end
end


text:
function [ a, b ] = p00_lim ( prob )

%*****************************************************************************80
%
%% P00_LIM returns the integration limits for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the desired test problem.
%
%    Output, real A, B, the limits of integration.
%
  if ( prob == 1 )
    [ a, b ] = p01_lim ( );
  elseif ( prob == 2 )
    [ a, b ] = p02_lim ( );
  elseif ( prob == 3 )
    [ a, b ] = p03_lim ( );
  elseif ( prob == 4 )
    [ a, b ] = p04_lim ( );
  elseif ( prob == 5 )
    [ a, b ] = p05_lim ( );
  elseif ( prob == 6 )
    [ a, b ] = p06_lim ( );
  elseif ( prob == 7 )
    [ a, b ] = p07_lim ( );
  elseif ( prob == 8 )
    [ a, b ] = p08_lim ( );
  elseif ( prob == 9 )
    [ a, b ] = p09_lim ( );
  elseif ( prob == 10 )
    [ a, b ] = p10_lim ( );
  elseif ( prob == 11 )
    [ a, b ] = p11_lim ( );
  elseif ( prob == 12 )
    [ a, b ] = p12_lim ( );
  elseif ( prob == 13 )
    [ a, b ] = p13_lim ( );
  elseif ( prob == 14 )
    [ a, b ] = p14_lim ( );
  elseif ( prob == 15 )
    [ a, b ] = p15_lim ( );
  elseif ( prob == 16 )
    [ a, b ] = p16_lim ( );
  elseif ( prob == 17 )
    [ a, b ] = p17_lim ( );
  elseif ( prob == 18 )
    [ a, b ] = p18_lim ( );
  elseif ( prob == 19 )
    [ a, b ] = p19_lim ( );
  elseif ( prob == 20 )
    [ a, b ] = p20_lim ( );
  elseif ( prob == 21 )
    [ a, b ] = p21_lim ( );
  elseif ( prob == 22 )
    [ a, b ] = p22_lim ( );
  elseif ( prob == 23 )
    [ a, b ] = p23_lim ( );
  elseif ( prob == 24 )
    [ a, b ] = p24_lim ( );
  elseif ( prob == 25 )
    [ a, b ] = p25_lim ( );
  elseif ( prob == 26 )
    [ a, b ] = p26_lim ( );
  elseif ( prob == 27 )
    [ a, b ] = p27_lim ( );
  elseif ( prob == 28 )
    [ a, b ] = p28_lim ( );
  elseif ( prob == 29 )
    [ a, b ] = p29_lim ( );
  elseif ( prob == 30 )
    [ a, b ] = p30_lim ( );
  elseif ( prob == 31 )
    [ a, b ] = p31_lim ( );
  elseif ( prob == 32 )
    [ a, b ] = p32_lim ( );
  elseif ( prob == 33 )
    [ a, b ] = p33_lim ( );
  elseif ( prob == 34 )
    [ a, b ] = p34_lim ( );
  elseif ( prob == 35 )
    [ a, b ] = p35_lim ( );
  elseif ( prob == 36 )
    [ a, b ] = p36_lim ( );
  elseif ( prob == 37 )
    [ a, b ] = p37_lim ( );
  elseif ( prob == 38 )
    [ a, b ] = p38_lim ( );
  elseif ( prob == 39 )
    [ a, b ] = p39_lim ( );
  elseif ( prob == 40 )
    [ a, b ] = p40_lim ( );
  elseif ( prob == 41 )
    [ a, b ] = p41_lim ( );
  elseif ( prob == 42 )
    [ a, b ] = p42_lim ( );
  elseif ( prob == 43 )
    [ a, b ] = p43_lim ( );
  elseif ( prob == 44 )
    [ a, b ] = p44_lim ( );
  elseif ( prob == 45 )
    [ a, b ] = p45_lim ( );
  elseif ( prob == 46 )
    [ a, b ] = p46_lim ( );
  elseif ( prob == 47 )
    [ a, b ] = p47_lim ( );
  elseif ( prob == 48 )
    [ a, b ] = p48_lim ( );
  elseif ( prob == 49 )
    [ a, b ] = p49_lim ( );
  elseif ( prob == 50 )
    [ a, b ] = p50_lim ( );
  elseif ( prob == 51 )
    [ a, b ] = p51_lim ( );
  elseif ( prob == 52 )
    [ a, b ] = p52_lim ( );
  elseif ( prob == 53 )
    [ a, b ] = p53_lim ( );
  elseif ( prob == 54 )
    [ a, b ] = p54_lim ( );
  elseif ( prob == 55 )
    [ a, b ] = p55_lim ( );
  elseif ( prob == 56 )
    [ a, b ] = p56_lim ( );
  elseif ( prob == 57 )
    [ a, b ] = p57_lim ( );
  elseif ( prob == 58 )
    [ a, b ] = p58_lim ( );
  elseif ( prob == 59 )
    [ a, b ] = p59_lim ( );
  elseif ( prob == 60 )
    [ a, b ] = p60_lim ( );
  elseif ( prob == 61 )
    [ a, b ] = p61_lim ( );
  elseif ( prob == 62 )
    [ a, b ] = p62_lim ( );
  elseif ( prob == 63 )
    [ a, b ] = p63_lim ( );
  elseif ( prob == 64 )
    [ a, b ] = p64_lim ( );
  elseif ( prob == 65 )
    [ a, b ] = p65_lim ( );
  elseif ( prob == 66 )
    [ a, b ] = p66_lim ( );
  elseif ( prob == 67 )
    [ a, b ] = p67_lim ( );
  elseif ( prob == 68 )
    [ a, b ] = p68_lim ( );
  elseif ( prob == 69 )
    [ a, b ] = p69_lim ( );
  elseif ( prob == 70 )
    [ a, b ] = p70_lim ( );
  elseif ( prob == 71 )
    [ a, b ] = p71_lim ( );
  elseif ( prob == 72 )
    [ a, b ] = p72_lim ( );
  elseif ( prob == 73 )
    [ a, b ] = p73_lim ( );
  elseif ( prob == 74 )
    [ a, b ] = p74_lim ( );
  elseif ( prob == 75 )
    [ a, b ] = p75_lim ( );
  elseif ( prob == 76 )
    [ a, b ] = p76_lim ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_LIM - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_LIM - Fatal error!' );
  end

  return
end


text:
function status = allEq( X, Y, THR )
% allEq() : rapidly assert X and Y are (almost exactly) the same
% Check if matrices X, Y are equal to within given tolerance
%   defined on the percent difference scale
%INPUT:
%   X, Y: N x D matrices
%   THR: numeric tolerance for equality (default 1e-8)
%OUTPUT
%  status : boolean
%     true if EVERY entry has perc diff <= THR
%     false o.w.

if ~exist( 'THR', 'var' )
    THR = 1e-8;
end

percDiff = abs( X - Y )./abs(X);
if max(percDiff) > THR
    status = false;
else
    status = true;
end

author: michaelchughes

repo: NPBayesHMM


text:
function [radps] = GHz2radps(GHz)
% Convert frequency from gigahertz to radians per second.
% Chad A. Greene 2012
radps = GHz*(2*pi)*1e+9;


text:
%% Read shape

[X,T] = readOff('../data/meshes/sphere_102.off');
M = getMeshData(X,T,10); % compute 10 LB eigenfunctions for fun

%% Set up Gaussian blur function

blurTime = .005; % if this gets too small, distances get noisy
blurSteps = 50;

% blur = @(x) blurOnMesh(x,M,blurTime,blurSteps)+1e-300; % faster than pre-factored?
% blurTranspose = @(x) blurOnMesh(x,M,blurTime,blurSteps,1)+1e-300;

% structure = prefactorMeshBlur(M,blurTime,blurSteps,5000);
% blur = @(x) prefactoredBlur(x,structure,0)+1e-300;
% blurTranspose = @(x) prefactoredBlur(x,structure,1)+1e-300;

h = blurTime/blurSteps;
nv = M.numVertices;

% Sphere is small enough that we can just write down the heat kernel explicitly
blurInverse = spdiags(M.areaWeights,0,nv,nv) - h*M.cotLaplacian;
mtx = full(blurInverse) \ diag(M.areaWeights);
mtx = mtx^blurSteps;

mtx(mtx<1e-50) = 1e-50;

blur = @(x) mtx*x;
blurTranspose = @(x) mtx'*x;

%% Make boundary

[~,fixedVerts(1)] = min(X(:,1));
[~,fixedVerts(2)] = min(X(:,2));
[~,fixedVerts(3)] = min(X(:,3));
[~,fixedVerts(4)] = max(X(:,1));
[~,fixedVerts(5)] = max(X(:,2));
[~,fixedVerts(6)] = max(X(:,3));

nFixed = length(fixedVerts);

fixedDistributions = zeros(M.numVertices,nFixed);

close all;
for i=1:nFixed
    fixedDistributions(fixedVerts(i),i) = 1/M.areaWeights(fixedVerts(i));
end

fixedDistributions = fixedDistributions+1e-10;
% fixedDistributions = blur(fixedDistributions);
showDescriptor(M,fixedDistributions(:,1));

%% Compute average entropy of data

averageEntropy = -mean(sum(bsxfun(@times,M.areaWeights,(fixedDistributions.*log(fixedDistributions))),1));

%% Solve barycenter problem using this machinery

edges = [2 1; 3 1; 4 1];
edgeWeights = [1 1 1 1 1 1];

targetEntropy = averageEntropy + 1;

result = convolutionalPropagation(edges, edgeWeights, [2 3 4], ...
    fixedDistributions(:,1:3),M.areaWeights,blur,blurTranspose,targetEntropy);

barycenter = convolutionalBarycenter(fixedDistributions(:,1:3),[1 1 1],...
    M.areaWeights,blur,blurTranspose,targetEntropy);

close all

f = subplot(1,3,1);
showDescriptor(M,sum(fixedDistributions(:,1:3),2),[],[],[],f);
title('Boundary distributions');
colorbar off;

f = subplot(1,3,2);
showDescriptor(M,result(:,1),[],[],[],f);
title('Barycenter computed using graph algorithm');
colorbar off;

f = subplot(1,3,3);
showDescriptor(M,barycenter,[],[],[],f);
title('Barycenter computed using specialized barycenter algorithm')
colorbar off;

%% Solve displacement problem using this machinery

k = 9;

edges = [(1:(k-1))' (2:k)'];
edges = [edges ; edges(:,2) edges(:,1)];
edgeWeights = ones(1,size(edges,1));

targetEntropy = averageEntropy + 1.5;

result = convolutionalPropagation(edges, edgeWeights, [1 k], ...
    fixedDistributions(:,[1 6]),M.areaWeights,blur,blurTranspose,targetEntropy);

figure

for i=1:k
    f = subplot(1,k,i);
    showDescriptor(M,result(:,i),[],[],[],f);
    colorbar off;
end

%% Solve soft mapping problem

edges = [T(:,1) T(:,2); T(:,2) T(:,3); T(:,3) T(:,1)];
edges = sort(edges,2);
edges = unique(edges,'rows');
edges = [edges; edges(:,2) edges(:,1)];

edgeWeights = ones(size(edges,1),1);

targetEntropy = averageEntropy+1.5;

result = convolutionalPropagation(edges, edgeWeights, fixedVerts, ...
                    fixedDistributions,M.areaWeights,blur,blurTranspose,targetEntropy);

%% Show result

testVerts = randperm(M.numVertices);
nTests = 5;
testVerts = testVerts(1:nTests);

figure;
for i=1:length(testVerts)
    f = subplot(2,length(testVerts),i);
    ind = zeros(M.numVertices,1);
    ind(testVerts(i)) = 1;
    showDescriptor(M,ind,[],[],[],f);
    colorbar off;
    colormap hot;

    f = subplot(2,length(testVerts),i+length(testVerts));
    showDescriptor(M,result(:,testVerts(i)),[],[],[],f);
    colorbar off;
    colormap hot;
end

author: gpeyre


text:
function imagescg(im)
% imagescg(im)
% Uses colormap('gray') for the imagecs
imagesc(im)
colormap('gray')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function updateScatterhistogram(obj, plotIndex)

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

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    [xsource, ysource, xoverlay, yoverlay] = findSourceAxis(obj,axIndex);

    [~, ~, groupName] = getTraceData(plotData);

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

    %-SET MAIN SCATTER PLOT-%

    %-set plotly layout-%
    updateMainScatterAxis(obj, plotIndex);
    updateTitle(obj, plotIndex);
    updateLegend(obj, plotIndex, groupName);
    obj.layout.barmode = 'overlay';
    obj.layout.bargap = 0.05;

    %-set plotly data-%
    updateMainScatter(obj, plotIndex);

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

    %-SET MARGINAL PLOTS-%

    %-set plotly layout-%
    updateXMarginalAxis(obj, plotIndex);
    updateYMarginalAxis(obj, plotIndex);

    %-set plotly data-%
    switch plotData.HistogramDisplayStyle
        case 'stairs'
            updateMarginalHistogram(obj, plotIndex, 'X');
            updateMarginalHistogram(obj, plotIndex, 'Y');
        case 'smooth'
            updateMarginalSmooth(obj, plotIndex, 'X');
            updateMarginalSmooth(obj, plotIndex, 'Y');
    end

    %-------------------------------------------------------------------------%
end

function updateMainScatter(obj, plotIndex)

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

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    [xSource, ySource] = findSourceAxis(obj,axIndex);

    %-get trace data-%
    traceIndex = plotIndex;
    [xData, yData, groupName] = getTraceData(plotData);

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

    %-SET EACH ALL TRACES-%

    for t = 1:length(xData)

        %---------------------------------------------------------------------%
        %-get current trace index-%
        if t > 1
            obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
            traceIndex = obj.PlotOptions.nPlots;
        end

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

        %-set current trace-%
        obj.data{traceIndex}.type = 'scatter';
        obj.data{traceIndex}.mode = 'markers';
        obj.data{traceIndex}.xaxis = sprintf('x%d', xSource);
        obj.data{traceIndex}.yaxis = sprintf('y%d', ySource);
        obj.data{traceIndex}.visible = strcmp(plotData.Visible,'on');

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

        %-set current trace data-%
        obj.data{traceIndex}.x = xData{t};
        obj.data{traceIndex}.y = yData{t};

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

        %-scatter marker-%
        childmarker = extractScatterhistogramMarker(plotData, t);
        obj.data{traceIndex}.marker = childmarker;

        %-legend-%

        if ~isempty(groupName)
            try
              obj.data{traceIndex}.name = char(groupName(t));
            catch
              obj.data{traceIndex}.name = char(string(groupName(t)));
            end

            obj.data{traceIndex}.legendgroup = obj.data{traceIndex}.name;
            obj.data{traceIndex}.showlegend = true;
        end

        %---------------------------------------------------------------------%
    end

    %-------------------------------------------------------------------------%
end

function updateMainScatterAxis(obj, plotIndex)

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

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    [xSource, ySource] = findSourceAxis(obj,axIndex);

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

    %-set x axis-%
    xaxis = getMainScatterAxis(plotData, 'X');
    xaxis.anchor = sprintf('y%d', xSource);
    obj.layout = setfield(obj.layout, sprintf('xaxis%d',xSource), xaxis);

    %-set y axis-%
    yaxis = getMainScatterAxis(plotData, 'Y');
    yaxis.anchor = sprintf('x%d', ySource);
    obj.layout = setfield(obj.layout, sprintf('yaxis%d',xSource), yaxis);

    %-------------------------------------------------------------------------%
end

function ax = getMainScatterAxis(plotData, axName)

    axisPos = plotData.Position;
    axisColor = 'rgba(0,0,0, 0.9)';
    axisLim = eval(sprintf('plotData.%sLimits', axName));
    axisPlot = eval(sprintf('plotData.%sData', axName));
    axisLabel = eval(sprintf('plotData.%sLabel', axName));

    switch axName
        case 'X'
            axisDomain = min([axisPos(1) sum(axisPos([1 3]))], 1);
        case 'Y'
            axisDomain = min([axisPos(2) sum(axisPos([2 4]))], 1);
    end

    %-general-%
    ax.domain = axisDomain;
    ax.linecolor = axisColor;
    ax.zeroline = true;
    ax.showgrid = false;
    ax.mirror = 'ticks';

    %-ticks-%
    ax.showticklabels = true;
    ax.ticks = 'inside';
    ax.tickfont.size = 1.2*plotData.FontSize;
    ax.tickcolor = axisColor;
    ax.tickfont.family = matlab2plotlyfont(plotData.FontName);

    %-label-%
    ax.title.text = axisLabel;
    if ~isempty(axisLabel), axisLabel = parseString(axisLabel); end
    ax.title.font.size = 1.2*plotData.FontSize;
    ax.title.font.color = 1.2*axisColor;
    ax.title.font.family = matlab2plotlyfont(plotData.FontName);

    %-range and ticklabels-%
    if ~iscategorical(axisPlot)
        ax.range = axisLim;
        ax.nticks = 10;
    else
        [axCateg, ~, axisPlot] = unique(axisPlot);
        ax.range = [min(axisPlot)-0.5, max(axisPlot)+0.5];
        ax.tickvals = 1:max(axisPlot);
        for n=1:length(axCateg), ax.ticktext{n} = char(axCateg(n)); end
    end
end

function updateMarginalHistogram(obj, plotIndex, axName)

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

    %-INITIALIZATIONS-%
    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    if strcmp(axName, 'X'), xySource = 1; else, xySource = 2; end
    xySource = obj.State.Figure.NumAxes + xySource;

    %-get trace data-%
    [xData, yData, groupName] = getTraceData(plotData);

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

    %-SET ALL TRACES-%

    for t = 1:length(xData)

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

        %-get current trace index-%
        obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
        traceIndex = obj.PlotOptions.nPlots;

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

        %-set current trace-%
        obj.data{traceIndex}.type = 'histogram';
        obj.data{traceIndex}.xaxis = sprintf('x%d', xySource);
        obj.data{traceIndex}.yaxis = sprintf('y%d', xySource);

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

        %-set current plot data-%
        obj.data{traceIndex}.x = xData{t};
        obj.data{traceIndex}.y = yData{t};

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

        %-set other trace properties-%
        traceColor = getStringColor(plotData.Color(t,:), 0.7);

        obj.data{traceIndex}.marker.color = traceColor;
        obj.data{traceIndex}.histnorm = 'probability';
        obj.data{traceIndex}.histfunc = 'count';
        obj.data{traceIndex}.showlegend = false;

        switch axName
            case 'X'
                obj.data{traceIndex}.orientation = 'v';
                try obj.data{traceIndex}.nbinsx = plotData.NumBins(1,t); end
            case 'Y'
                obj.data{traceIndex}.orientation = 'h';
                try obj.data{traceIndex}.nbinsy = plotData.NumBins(2,t); end
        end

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

        %-link legend-%
        if ~isempty(groupName)
            try
                obj.data{traceIndex}.name = char(groupName(t));
            catch
                obj.data{traceIndex}.name = char(string(groupName(t)));
            end

            obj.data{traceIndex}.legendgroup = obj.data{traceIndex}.name;
        end

        %---------------------------------------------------------------------%
    end

    %-------------------------------------------------------------------------%
end

function updateMarginalSmooth(obj, plotIndex, axName)

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

    %-INITIALIZATIONS-%
    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    if strcmp(axName, 'X'), xySource = 1; else, xySource = 2; end
    xySource = obj.State.Figure.NumAxes + xySource;

    %-get trace data-%
    [xData, yData, groupName] = getTraceData(plotData);
    axisLim = getAxisLim(plotData, axName);
    evalPoints = linspace(axisLim(1), axisLim(2), 500);

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

    %-SET ALL TRACES-%

    for t = 1:length(xData)

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

        %-get current trace index-%
        obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
        traceIndex = obj.PlotOptions.nPlots;

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

        %-set current trace-%
        obj.data{traceIndex}.type = 'scatter';
        obj.data{traceIndex}.mode = 'lines';
        obj.data{traceIndex}.xaxis = sprintf('x%d', xySource);
        obj.data{traceIndex}.yaxis = sprintf('y%d', xySource);

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

        %-get current trace data-%
        if strcmp(axName, 'X')
            [ySmooth, xSmooth] = ksdensity(xData{t}, evalPoints);
        elseif strcmp(axName, 'Y')
            [xSmooth, ySmooth] = ksdensity(yData{t}, evalPoints);
        end

        %-set current plot data-%
        obj.data{traceIndex}.x = xSmooth;
        obj.data{traceIndex}.y = ySmooth;

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

        %-set other trace properties-%
        traceColor = getStringColor(plotData.Color(t,:), 0.7);
        lineStyle = plotData.LineStyle(t);

        obj.data{traceIndex}.line.color = traceColor;
        obj.data{traceIndex}.line.width = 2*plotData.LineWidth(t);
        obj.data{traceIndex}.line.dash = getLineDash(lineStyle);
        obj.data{traceIndex}.showlegend = false;

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

        %-link legend-%
        if ~isempty(groupName)
            try
                obj.data{traceIndex}.name = char(groupName(t));
            catch
                obj.data{traceIndex}.name = char(string(groupName(t)));
            end

            obj.data{traceIndex}.legendgroup = obj.data{traceIndex}.name;
        end

        %---------------------------------------------------------------------%
    end

    %-------------------------------------------------------------------------%
end

function updateXMarginalAxis(obj, plotIndex)
    %-------------------------------------------------------------------------%

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    % [xSource, ySource] = findSourceAxis(obj,axIndex);

    xySource = obj.State.Figure.NumAxes + 1;

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

    %-set x axis-%
    xaxis = getXMarginalAxis(plotData, 'X');
    xaxis.anchor = sprintf('y%d', xySource);
    obj.layout = setfield(obj.layout, sprintf('xaxis%d',xySource), xaxis);

    %-set y axis-%
    yaxis = getXMarginalAxis(plotData, 'Y');
    yaxis.anchor = sprintf('x%d', xySource);
    obj.layout = setfield(obj.layout, sprintf('yaxis%d',xySource), yaxis);

    %-------------------------------------------------------------------------%
end

function ax = getXMarginalAxis(plotData, axName)

    switch axName
        case 'X'
            ax.showline = true;
            ax.linecolor = 'black';
            ax.range = getAxisLim(plotData, 'X');
        case 'Y'
            ax.showline = false;
    end

    ax.domain = getXMarginalDomain(plotData, axName);
    ax.showgrid = false;
    ax.showticklabels = false;
    ax.zeroline = false;
end

function axisDomain = getXMarginalDomain(plotData, axName)
    axisPos = plotData.Position;
    plotLocation = plotData.ScatterPlotLocation;
    isTitle = ~isempty(plotData.Title);

    switch axName
        case 'X'
            axisDomain = min([axisPos(1) sum(axisPos([1,3]))], 1);
        case 'Y'
            if ~isempty(strfind(plotLocation, 'South'))
                yo = axisPos(2) + axisPos(4) + 0.01;
                if isTitle, h=0.9-yo; else, h = 0.96 - yo; end
            elseif ~isempty(strfind(plotLocation, 'North'))
                yo = 0.02; h = axisPos(2)*0.7-yo;
            end

            axisDomain = min([yo yo+h], 1);
    end
end

function updateYMarginalAxis(obj, plotIndex)
    %-------------------------------------------------------------------------%

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    % [xSource, ySource] = findSourceAxis(obj,axIndex);

    xySource = obj.State.Figure.NumAxes + 2;

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

    %-set x axis-%
    xaxis = getYMarginalAxis(plotData, 'X');
    xaxis.anchor = sprintf('y%d', xySource);
    obj.layout = setfield(obj.layout, sprintf('xaxis%d',xySource), xaxis);

    %-set y axis-%
    yaxis = getYMarginalAxis(plotData, 'Y');
    yaxis.anchor = sprintf('x%d', xySource);
    obj.layout = setfield(obj.layout, sprintf('yaxis%d',xySource), yaxis);

    %-------------------------------------------------------------------------%
end

function ax = getYMarginalAxis(plotData, axName)

    switch axName
        case 'X'
            ax.showline = false;
        case 'Y'
            ax.showline = true;
            ax.linecolor = 'black';
            ax.range = getAxisLim(plotData, 'Y');
    end

    ax.domain = getYMarginalDomain(plotData, axName);
    ax.showgrid = false;
    ax.showticklabels = false;
    ax.zeroline = false;
end

function axisDomain = getYMarginalDomain(plotData, axName)
    axisPos = plotData.Position;
    plotLocation = plotData.ScatterPlotLocation;

    switch axName
        case 'X'
            if ~isempty(strfind(plotLocation, 'West'))
                xo = axisPos(1) + axisPos(3) + 0.01;
                w = 0.96-xo;
            elseif ~isempty(strfind(plotLocation, 'East'))
                xo = 0.02; w = axisPos(1)*0.7-xo;
            end

            axisDomain = min([xo xo+w], 1);

        case 'Y'

            axisDomain = min([axisPos(2) sum(axisPos([2,4]))], 1);

    end
end

function lineDash = getLineDash(lineStyle)
    switch lineStyle
        case '-'
            lineDash = 'solid';
        case '--'
            lineDash = 'dash';
        case ':'
            lineDash = 'dot';
        case '-.'
            lineDash = 'dashdot';
    end
end

function axisLim = getAxisLim(plotData, axName)
    axisLim = eval(sprintf('plotData.%sLimits', axName));
    axisPlot = eval(sprintf('plotData.%sData', axName));

    if iscategorical(axisPlot)
        axisPlot = eval(sprintf('plotData.%sData', axName));
        [~, ~, axisPlot] = unique(axisPlot);
        axisLim = [min(axisPlot)-0.5, max(axisPlot)+0.5];
    end
end

function [xData, yData, groupName] = getTraceData(plotData)
    %-parcing data-%
    xPlot = plotData.XData;
    yPlot = plotData.YData;

    if iscategorical(xPlot), [~, ~, xPlot] = unique(xPlot); end
    if iscategorical(yPlot), [~, ~, yPlot] = unique(yPlot); end

    xData = {}; yData = {};
    groupData = plotData.GroupData;
    isByGroups = ~isempty(groupData);
    groupName = {};

    if isByGroups
        if iscellstr(groupData), groupData = string(groupData); end
        groupName = unique(groupData,'stable');

        for g = 1:length(groupName)
            groudInd = groupData == groupName(g);
            xData{g} = xPlot(groudInd);
            yData{g} = yPlot(groudInd);
        end

        if isnumeric(groupName), groupName=num2str(groupName); end;
    else
        xData{1} = xPlot;
        yData{1} = yPlot;
    end
end

function updateTitle(obj, plotIndex)
    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    [xSource, ~] = findSourceAxis(obj,axIndex);
    isTitle = ~isempty(plotData.Title);

    obj.layout.annotations{1}.text = plotData.Title;
    obj.layout.annotations{1}.showarrow = false;

    if isTitle
        titleText = sprintf('%s', parseString(plotData.Title));
        titleFamily = matlab2plotlyfont(plotData.FontName);
        xaxis = eval(sprintf('obj.layout.xaxis%d', xSource));

        obj.layout.annotations{1}.text = titleText;
        obj.layout.annotations{1}.x = mean(xaxis.domain);
        obj.layout.annotations{1}.y = 0.96;
        obj.layout.annotations{1}.xref = 'paper';
        obj.layout.annotations{1}.yref = 'paper';
        obj.layout.annotations{1}.yanchor = 'top';
        obj.layout.annotations{1}.xanchor = 'middle';

        obj.layout.annotations{1}.font.color = 'black';
        obj.layout.annotations{1}.font.family = titleFamily;
        obj.layout.annotations{1}.font.size = 1.5*plotData.FontSize;
    end
end

function updateLegend(obj, plotIndex, groupName)
    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    isTitle = isempty(plotData.Title);

    if ~isempty(groupName)
        fontFamily = matlab2plotlyfont(plotData.FontName);
        legTitle = plotData.LegendTitle;
        plotLocation = plotData.ScatterPlotLocation;

        obj.layout.showlegend = true;
        obj.layout.legend.xref = 'paper';
        obj.layout.legend.valign = 'middle';
        obj.layout.legend.borderwidth = 1;
        obj.layout.legend.bordercolor = 'rgba(0,0,0,0.2)';
        obj.layout.legend.font.size = 1.0*plotData.FontSize;
        obj.layout.legend.font.family = fontFamily;

        if ~isempty(legTitle) > 0
            legTitle = sprintf('%s', parseString(legTitle));

            obj.layout.legend.title.text = legTitle;
            obj.layout.legend.title.side = 'top';
            obj.layout.legend.title.font.size = 1.2*plotData.FontSize;
            obj.layout.legend.title.font.color = 'black';
            obj.layout.legend.title.font.family = fontFamily;
        end

        if ~isempty(strfind(plotLocation, 'SouthWest'))
            obj.layout.legend.x = 0.96;
            obj.layout.legend.y = 0.96;
            obj.layout.legend.xanchor = 'right';
            obj.layout.legend.yanchor = 'top';
        elseif ~isempty(strfind(plotLocation, 'NorthEast'))
            obj.layout.legend.x = 0.02;
            obj.layout.legend.y = 0.02;
            obj.layout.legend.xanchor = 'left';
            obj.layout.legend.yanchor = 'bottom';
        end
    end
end

author: plotly


text:
function obj = coef_Jacob3_qt_pTop_func_new(in1)
coef_J3 = in1(:,3);
coef_J6 = in1(:,6);
coef_J8 = in1(:,8);
coef_J9 = in1(:,9);
coef_J16 = in1(:,16);
coef_J18 = in1(:,18);
coef_J19 = in1(:,19);
coef_J25 = in1(:,25);
coef_J26 = in1(:,26);
coef_J27 = in1(:,27);
coef_J28 = in1(:,28);
coef_J29 = in1(:,29);
coef_J30 = in1(:,30);
coef_J31 = in1(:,31);
coef_J38 = in1(:,38);
coef_J40 = in1(:,40);
coef_J41 = in1(:,41);
coef_J47 = in1(:,47);
coef_J48 = in1(:,48);
coef_J49 = in1(:,49);
coef_J50 = in1(:,50);
coef_J51 = in1(:,51);
coef_J52 = in1(:,52);
coef_J53 = in1(:,53);
coef_J59 = in1(:,59);
coef_J60 = in1(:,60);
coef_J61 = in1(:,61);
coef_J62 = in1(:,62);
coef_J63 = in1(:,63);
coef_J64 = in1(:,64);
coef_J65 = in1(:,65);
coef_J66 = in1(:,66);
coef_J67 = in1(:,67);
coef_J68 = in1(:,68);
coef_J69 = in1(:,69);
coef_J70 = in1(:,70);
obj = [coef_J3,coef_J6,coef_J8.*2.0,coef_J9,coef_J16,coef_J18.*2.0,coef_J19,coef_J25.*3.0,coef_J26.*2.0,coef_J27,coef_J28,coef_J29,coef_J30,coef_J31,coef_J38,coef_J40.*2.0,coef_J41,coef_J47.*3.0,coef_J48.*2.0,coef_J49,coef_J50,coef_J51,coef_J52,coef_J53,coef_J59.*4.0,coef_J60.*3.0,coef_J61.*2.0,coef_J62.*2.0,coef_J63.*2.0,coef_J64.*2.0,coef_J65.*2.0,coef_J66,coef_J67,coef_J68,coef_J69,coef_J70];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function model = forestTrain(X, Y, opts)
    % X is NxD, where rows are data points
    % for convenience, for now we assume X is 0 mean unit variance. If it
    % isn't, preprocess your data with
    %
    % X= bsxfun(@rdivide, bsxfun(@minus, X, mean(X)), var(X) + 1e-10);
    %
    % If this condition isn't satisfied, some weak learners won't work out
    % of the box in current implementation.
    %
    % Y is discrete Nx1 vector of labels
    % model can be plugged into forestTest()
    %
    % decent default opts are:
    % opts.depth= 9;
    % opts.numTrees= 100; %(but more is _ALWAYS_ better, monotonically, no exceptions)
    % opts.numSplits= 5;
    % opts.classifierID= 2
    % opts.classifierCommitFirst= true;
    %
    % which means use depth 9 trees, train 100 of them, use 5 random
    % splits when training each weak learner. The last option controls
    % whether each node commits to a weak learner at random and then trains
    % it best it can, or whether each node tries all weak learners and
    % takes the one that does best. Empirically, it appears leaving this as
    % true (default) gives slightly better looking results.
    %

    numTrees= 100;
    verbose= false;

    if nargin < 3, opts= struct; end
    if isfield(opts, 'numTrees'), numTrees= opts.numTrees; end
    if isfield(opts, 'verbose'), verbose= opts.verbose; end

    treeModels= cell(1, numTrees);
    for i=1:numTrees

        treeModels{i} = treeTrain(X, Y, opts);

        % print info if verbose
        if verbose
            p10= floor(numTrees/10);
            if mod(i, p10)==0 || i==1 || i==numTrees
                fprintf('Training tree %d/%d...\n', i, numTrees);
            end
        end
    end

    model.treeModels = treeModels;
end

author: karpathy

repo: Random-Forest-Matlab


text:
function y = beta_lpdf(x,a,b)
%BETA_LPDF    Beta log-probability density function (lpdf).
%
%   Y = BETA_LPDF(X,A,B) Returns the log of the Beta pdf with
%   parameters A and B, at the values in X.
%
%   The size of Y is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default value for A and B is 1.

% Copyright (c) 2005 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3,
  a = 1;
end

if nargin < 2;
  b = 1;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y= (a-1).*log(x) +(b-1).*log(1-x) -betaln(a,b);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% largely drawn from run_readlsl of BCILAB
function lslin(handles)
% find stream
result = lsl_resolve_bypred(handles.lsllib, ['name=''' handles.streamName '''']);

% open inlet
inlet = lsl_inlet(result{1});
% info = inlet.info();

% create online stream data structure in base workspace (using appropriate meta-data)
onl_newstream(parseStreamName(handles.streamName), 'srate', result{1}.nominal_srate(), ...
    'chanlocs', {handles.chanlocs.labels}, 'buffer_len', 10);

% state variables for recursive least squares jitter correction
P = 1e10*eye(2);        % precision matrix (inverse covariance matrix of predictors)
w = [0 0]';             % linear regression coefficients [offset,slope]
lam = 2^(-1/(128 * 30)); % forget factor in RLS calculation
n = 0;                  % number of samples observed so far
numeric_offset = [];    % time-stamp offset to keep numerics healthy; will be initialized with first measured time stamp

% start reading
onl_read_background(parseStreamName(handles.streamName), @read_data, 20);

    % reads from inlet
    function results = read_data()
        [chunk, stamps] = inlet.pull_chunk();
        data_clock = inlet.time_correction([], 'median', 30);
        stamps = stamps + data_clock;
        stamps = update_regression(stamps);
        chunk = double(chunk);

        % this is the source of grief
        % taking only the last timestamp allows REST to run but the
        % timestamps are garbage
        % giving all the timestamps soft errors in onl_append. onl_append
        % seems to expect only one value so perhaps this idea of many
        % timestamps if not correct.
        try %#ok<TRYNC>
            stamps = stamps(end);
        end
        results = {chunk, stamps};
    end

    % perform RLS block update of regression coefficients
    % this is a regression from sample index onto timestamp of the sample
    function y = update_regression(y)
        if ~isempty(y)
            % sanitize numerics (all done relative to the first observed time stamp)
            if isempty(numeric_offset)
                numeric_offset = y(1); end
            y = y - numeric_offset;
            % define predictor matrix (bias, sample index)
            X = [ones(1,length(y)); n + (1:length(y))];
            n = n + length(y);
            % apply updates...
            for t=1:length(y)
                u = X(:,t);
                d = y(t);
                pi = u'*P;
                gam = lam + pi*u;
                k = pi'/gam;
                al = d - w'*u;
                w = w + k*al;
                Pp = k*pi;
                P = (1/lam)*(P-Pp);
            end
            % predict y
            y = w'*X + numeric_offset;
        end
    end

    % parse streamname
    function streamnames = parseStreamName(streamnames)
        if ~isvarname(streamnames)
            streamnames = streamnames(~ismember(streamnames,['-' ' ']));
        end
    end

end

repo: REST


text:
function Cct=clg2cct(Clg,lat,lon)
% CLG2CCT  Convert local geodetic covariance matrix to CT.
%   Vectorized. Requires rotlg2ct. See also CCT2CLG.
% Version: 2011-04-05
% Useage:  Cct=clg2cct(Clg,lat,lon)
% Input:   Clg - LG covariance matrix
%          lat - vector of station latitudes (rad)
%          lon - vector of station longitudes (rad)
% Output:  Cct - CT covariance matrix

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

n=length(lat);
if (n*3 ~= max(size(Clg)) )
  error('Size of lat,lon does not match size of Clg');
end

for i=1:n
  sinlat=sin(lat(i));
  coslat=cos(lat(i));
  sinlon=sin(lon(i));
  coslon=cos(lon(i));
  Ji=rotlg2ct(lat(i),lon(i));
  indi=(i-1)*3+[1:3];
  for j=1:n
    sinlat=sin(lat(j));
    coslat=cos(lat(j));
    sinlon=sin(lon(j));
    coslon=cos(lon(j));
    Jj=rotlg2ct(lat(j),lon(j));
    indj=(j-1)*3+[1:3];
    Cct(indi,indj)=Ji*Clg(indi,indj)*Jj';
  end
end


text:
% blob_unity1.m
% find blob parameters that make the best "partition of unity"
% jeff fessler

% 1d
% J m alpha	alpha/J	(max/min-1)%
% 2 2 2.502	1.251	0.0192		1 min
% 3 2 2.472	0.824	0.0021		2 min
% 3 2 7.464	2.488	0.0587776%	"
% 4 2 2.461	0.61525 0.000315658%	3 min
% 4 2 8.687	2.17175	0.0102284%	"
% 4 2 11.17	2.7925	0.00177428%	"

% 2d
% J m alpha	alpha/J	(max/min-1)%
% 4 2 7.888	1.972	0.0411645
% 4 2 10.83	2.7075	0.0180456%

% 3d
% J m alpha	alpha/J	(max/min-1)%
% 4 2 10.4	2.6	0.0290578	and mae and rmse < 0.01%.
% for CT this is a fraction of a HU!

if ~isvar('J')
	J = 4; % a = J/2

	alf_list = linspace(0.5*J, 3.0*J, 26+0*101);
	alf_list = linspace(7, 12, 5*10+1); % 2d
	alf_list = linspace(10.5, 11, 5*10+1); % 2d
	alf_list = linspace(2, 12, 10001); % 1d
	alf_list = linspace(7.0, 8.0, 1001); % 1d fine, J=3
alf_list = linspace(0.1, 20, 1001); % 1d coarse
	alf_list = linspace(8.0, 9.0, 1001); % 1d fine, J=4
%	alf_list = linspace(10.3, 10.5, 3); % 3d
%	alf_list = 10.4; % 3d
%	m_list = linspace(1.9, 2.1, 21);
	m_list = 2.0;
	%alf_list = 2.34 * J;
	%m_list = 2;
	[aa mm] = ndgrid(alf_list, m_list);

	x1 = linspace(0, J/2, 101)';
	x2 = linspace(0, J/2, 101)';
	x3 = linspace(0, J/2, 101)';
	x2 = 0; % 1d
	x3 = 0; % 2d
	[xx1 xx2 xx3] = ndgrid(x1, x2, x3);

	% get all j that affect x in [0,J/2]
	j1max = ceil(J-1);
	j1min = floor(-J/2+1);
	ndim = 1 + (length(x2) > 1);
	nj1 = j1max-j1min+1;

	if ndim >= 3
		nj3 = nj1;
		j3min = j1min;
		j3max = j1max;
	else
		nj3 = 1;
		j3min = 0;
		j3max = 0;
	end

	if ndim >= 2
		nj2 = nj1;
		j2min = j1min;
		j2max = j1max;
	else
		nj2 = 1;
		j2min = 0;
		j2max = 0;
	end
end

%
% precompute r samples for each j offset
%
if ~isvar('rr'), disp 'do rr'
	rr = zeros([length(x1) length(x2) length(x3) nj1 nj2 nj3]);
	for j3=j3min:j3max
		for j2=j2min:j2max
			for j1=j1min:j1max
				rr(:,:,:,j1-j1min+1,j2-j2min+1,j3-j3min+1) = ...
				sqrt((xx1-j1).^2 + (xx2-j2).^2 + (xx3-j3).^2);
			end
		end
	end
	%im(x1, x2, rr)
end

if ~isvar('bad'), disp 'do bad'
	bad = zeros(size(mm));
	bb = zeros(length(x1), length(x2), length(x3), nj1*nj2*nj3);

	for ii=1:numel(aa)
		ticker(mfilename, ii, numel(aa))
		kb_a = aa(ii);
		kb_m = mm(ii);

		bb = kaiser_bessel(rr, J, kb_a, kb_m);
		bsum = sum(sum(sum(bb, 3), 4), 5);
		bad(ii) = max(bsum(:)) / min(bsum(:));
	end
	bad = abs(bad-1);
end

if 1
	im clf, im(121, alf_list/J, m_list, bad), cbar
	axis normal
	hold on
	ibest = imin(bad, 2);
	plot(alf_list(ibest(1))/J, m_list(ibest(2)), '*')
	hold off
	xlabel '\alpha/J', ylabel 'm'
end

if ~isvar('bmean')
	kb_a = alf_list(ibest(1));
	kb_m = m_list(ibest(2));

	bb = kaiser_bessel(rr, J, kb_a, kb_m);
	bsum = sum(sum(sum(bb, 3), 4), 5);
	bmean = mean(bsum(:));
end

printf('J m alpha alpha/J ratio-1')
printf('%g %g %g %g %g%%', ...
	J, kb_m, kb_a, kb_a / J, (max(bsum(:)) / min(bsum(:))-1)*100)
printf('mea=%g%%', mean(abs(bsum(:)-bmean))/bmean * 100)
printf('rmse=%g%%', sqrt(mean(abs(bsum(:)-bmean).^2))/bmean * 100)

if 1
	subplot(122)
	if ndim == 1
		plot(x1, reshape(bb, length(x1), []), '--', x1, bsum, '-')
		clf, semilogy(alf_list, bad, '.-'), xlabel '\alpha'
	else
		im(x1, x2, bsum/bmean, 'bsum/bmean'), cbar
	end
end

if 0
	clf
	for j1=1:nj1
		for j2=1:nj2
			subplot(nj1, nj2, j1 + (j2-1)*nj1)
			if ndim == 1
				plot(x1, bb(:,1,j1,j2))
			else
				im(x1, x2, bb(:,:,j1,j2))
			end
		end
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function SO3F = plus(SO3F1,SO3F2)
% overloads |SO3F1 + SO3F2|
%
% Syntax
%   SO3F = SO3F1 + SO3F2
%   SO3F = a + SO3F1
%   SO3F = SO3F1 + a
%
% Input
%  SO3F1, SO3F2 - @SO3Fun
%  a - double
%
% Output
%  SO3F - @SO3Fun
%

% uniform component + SO3Fun
if isa(SO3F1,'SO3FunRBF') && isempty(SO3F1.center) && isa(SO3F2,'SO3Fun')
  SO3F = SO3F1.c0 + SO3F2;
  return
end
if isa(SO3F2,'SO3FunRBF') && isempty(SO3F2.center) && isa(SO3F1,'SO3Fun')
  SO3F = SO3F1 + SO3F2.c0;
  return
end

ensureCompatibleSymmetries(SO3F1,SO3F2);
SO3F = SO3FunComposition(SO3F1, SO3F2);

end

author: mtex-toolbox

repo: mtex

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


text:
function [fullname,ftime]=findobsrfile(filepath,sitename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global gls
fullname='';
ftime.ts=gls.gtime; ftime.te=gls.gtime; ftime.week=0; ftime.sow=0;
ftime.year=0; ftime.doy=0;

% search the obsr file
fileExt = '*.*O';
files = dir(fullfile(filepath,fileExt));
if isempty(files)
    fileExt = '*.*o';
    files = dir(fullfile(filepath,fileExt));
end
nf=size(files,1);
if nf==0,return,end
for i=1:nf
    if any(sitename)
        if ~isempty(strfind(files(i).name,sitename))&&...
                isempty(strfind(files(i).name,'base'))
            fullname=[filepath,files(i).name];
            break;
        end
    else
        if isempty(strfind(files(i).name,'base'))
            fullname=[filepath,files(i).name];
            break;
        end
    end
end

% get the observe time from obsr file
if ~isempty(fullname)
    ftime=getobstime(fullname,ftime);
end

return

repo: GINav


text:
%Copyright 2008-2009  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
%Copyright 2008-2009  David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
%
%THE BSD LICENSE
%
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions
%are met:
%
%1. Redistributions of source code must retain the above copyright
%   notice, this list of conditions and the following disclaimer.
%2. Redistributions in binary form must reproduce the above copyright
%   notice, this list of conditions and the following disclaimer in the
%   documentation and/or other materials provided with the distribution.
%
%THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
%INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
%NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
%THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function [indices, dists] = flann_search(data, testset, n, search_params)
%NN_SEARCH  Fast approximate nearest neighbors search
%
% Performs a fast approximate nearest neighbor search using an
% index constructed using flann_build_index or directly a
% dataset.

% Marius Muja, January 2008

    algos = struct( 'linear', 0, 'kdtree', 1, 'kmeans', 2, 'composite', 3, 'saved', 254, 'autotuned', 255 );
    center_algos = struct('random', 0, 'gonzales', 1, 'kmeanspp', 2 );
    log_levels = struct('none', 0, 'fatal', 1, 'error', 2, 'warning', 3, 'info', 4);
    function value = id2value(map, id)
        fields = fieldnames(map);
        for i = 1:length(fields),
            val = cell2mat(fields(i));
            if map.(val) == id
                value = val;
                break;
            end
        end
    end
    function id = value2id(map,value)
        id = map.(value);
    end

    default_params = struct('algorithm', 'kdtree' ,'checks', 32, 'eps', 0.0, 'sorted', 1, 'max_neighbors', -1, 'cores', 1, 'trees', 4, 'branching', 32, 'iterations', 5, 'centers_init', 'random', 'cb_index', 0.4, 'target_precision', 0.9,'build_weight', 0.01, 'memory_weight', 0, 'sample_fraction', 0.1, 'log_level', 'warning', 'random_seed', 0);

    if ~isstruct(search_params)
        error('The "search_params" argument must be a structure');
    end

    params = default_params;
    fn = fieldnames(search_params);
    for i = [1:length(fn)],
        name = cell2mat(fn(i));
        params.(name) = search_params.(name);
    end
    if ~isnumeric(params.algorithm),
        params.algorithm = value2id(algos,params.algorithm);
    end
    if ~isnumeric(params.centers_init),
        params.centers_init = value2id(center_algos,params.centers_init);
    end
    if ~isnumeric(params.log_level),
        params.log_level = value2id(log_levels,params.log_level);
    end

    if (size(data,1)==1 && size(data,2)==1)
        % we already have an index
        [indices,dists] = nearest_neighbors('index_find_nn', data, testset, n, params);
    else
        % create the index and search
        [indices,dists] = nearest_neighbors('find_nn', data, testset, n, params);
    end
end

repo: gspbox


text:
vm.f=@(A,B).5*norm(A(:)+B(:)-Z(:))^2;
vm.g=@(A,B)A+B-Z;

author: andrewssobral

repo: lrslibrary


text:
function varargout = domainCheck(varargin)
%DOMAINCHECK   True if the domains of two CHEBFUN2 objects are the same.
%   DOMAINCHECK(F, G) returns TRUE if the domains of the two
%   CHEBFUN2 objects F and G coincide up to a tolerance depending on their
%   horizontal scales or if both F and G are empty CHEBFUN objects.
%
% See also CHEBFUN/DOMAINCHECK.

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

[varargout{1:nargout}] = domainCheck@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun