text:
classdef ElasticMaterial < Material

    properties (Access = public)
        C
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function P=mg_ns_prolong_step(nc,x,y,outbnd)
%mg_ns_prolong_step  GMG prolongation operator for step domain (Navier-Stokes)
%   P=mg_ns_prolong_step(nc,x,y,outbnd)
%   input
%          nc      grid parameter
%          x       x coordinate vector for coarse grid
%          y       y coordinate vector for coarse grid
%          outbnd  location of outflow boundary
%   output
%          P       prolongation operator
%
%   IFISS function: HCE; 5 May 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

nelf=2^nc;nelc=nelf/2;nelcx=outbnd*nelc;
nnf=(nelc+1)*((outbnd+1)*nelc+1)+nelc*(outbnd*nelc+1);
nnc=(nelc/2+1)*(((outbnd+1)/2)*nelc+1)+(nelc/2)*(outbnd*nelc/2+1);

% split coordinates
x1=x(1:nelc+1); x2=x(nelc+1:end);
y1=y(nelc+1:end); y2=y;

% block size
nnf1=nelc*(nelc+1);
nnc1=nelc/2*(nelc/2+1);

start_col=(nelc/2+1)*(nelc/2+1);

% Prolongation matrix has the form
% P = [P_1, 0, P_2]
%     [   0 ,  P_3]

% nonzero components of upper block, for upper left of step and coupling
[P_1,P_2]=mg_ellblock(nelc,nelc,x1,y1);
% P(1:nnf1,1:nnc1)=P_1;                % upper block
% P(1:nnf1,nnc-start_col+1:nnc)=P_2;   % coupling

% nonzero component of lower block
P_3=mg_prolong(nelcx,nelf,x2,y2);
% P(nnf1+1:end,nnc1+1:end)=P_3;

P = [[P_1,sparse(nnf1,nnc-start_col-nnc1),P_2];[sparse(nnf-nnf1,nnc1),P_3]];

author: johannesgerer

repo: jburkardt-m


text:
function computeAllModelChoice_HN(pathSet,pathModels,training,param,maxOrder,nBoot,imbalance,seed)
% -------------------------------------------------------------------------
% function computeAllModelChoice_HN(pathSet,pathModels,training,param,maxOrder,nBoot,imbalance,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection for a given feature set type
% and for all experiments with different degrees of freedom. See ref. [1,2]
% for more details.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathSet: Full path to the Feature Set folder of the corresponding experiment.
%             --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/FSET'
% 2. pathModels: Full path to the Models folder of the corresponding experiment.
%                --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/MODELS'
% 3. training: Structure defining all parameters for the given experiment
%              to perform. See masterScript_HN.m for more details.
% 4. param: Cell of two strings, defining 1) The feature set type/name; and
%           2) the outcome to model
% 5. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 6. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 7. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 8. batchNum: (optional input). If present, integer that specifies the
%              batch number for parallelization purposes.
%              --> Ex: 6
% -------------------------------------------------------------------------
% OUTPUTS: Mutivariable models are saved in a folder named 'MODELS' in the
% corresponding folder of "pathExperiment"
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
warning off

cd(pathSet)
fSetName = param{1};
outcomeName = param{2};
outcome = training.(outcomeName).outcome;
fSet = load(['FSET_',fSetName,'_',outcomeName]); fSet = struct2cell(fSet); fSet = fSet{1};

tic
fprintf(['\n--> SELECTING FEATURES (MODEL ORDERS OF 1 to % u) FOR "',fSetName,'" FEATURE SET AND "',outcomeName,'" OUTCOME ... '],maxOrder)
[models] = featureSelection_HN(fSet.Data,outcome,maxOrder,nBoot,fSet.Info,imbalance,seed);
cd(pathModels), save(['MODELS_',fSetName,'_',outcomeName],'models')
fprintf('DONE!\n'), toc

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [u,p] = StokesIUzawa(u,p,f,g,A,B,auxMat,elem,smootherOpt)
%% STOKESIUzawa Inexact Uzawa smoother for the Stokes eqns.
%
%  In matrix form, we solve the Stokes equations
%
%               |A B'| |u| = |f|
%               |B 0 | |p| = |g|
%
%  Define
%            M = |I -S_u^-1 B'|   T = |S_u    0|
%                |     I      |       |B   -S_p|
%  The matrix form DGS update can be written as
%
%    |uk+1|   |uk|            |ru|     with ru = f-Auk -B'pk,
%    |    | = |  | + M*inv(T)*|  |
%    |pk+1|   |pk|            |rp|     and  rp = 0-Buk
%
%  Created by Ming Wang (with discussion of Long Chen) at Jan, 2012.
%  Revised at Sept 2012.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

%% Parameters
itStep = smootherOpt.itStep;
Su = auxMat.DA; % DA = 2*diag(A)
Bt = auxMat.Bt;
BinvDABt = auxMat.BinvDABt; % invDA = 1/(2*diag(A));
% mg options
mgoption.printlevel = 0;
mgoption.maxIt = 1;
mgoption.N0=50;
mgoption.solver = 'Vcycle';
mgoption.mu = 2;

%% DGS relaxation step
for k = 1: itStep
    % Step 1: relax Momentum eqns
    if isempty(u) && isempty(p)
        ru = f;
        u = zeros(size(f));
        p = zeros(size(g));
    else
        ru = f-Bt*p-A*u;
    end
    u = u + ru./Su;
    rp = g-B*u;
    rp = rp - mean(rp);
    % Step 2: relax transformed Continuity eqns
    dp = mg(BinvDABt,-rp,elem,mgoption);
    % Step 3: transform the correction back to the original variables.
    p = p + dp;
    u = u - (Bt*dp)./Su;
end
end

repo: ifem


text:
%
% Perform an iteration of MCMC
% for re-sampling all of the type-level
% variables, excluding relation types
%
function mcmc_iter_type(MH,M,lib)

    % shape types
    for sid=1:M.ns
        for bid=1:M.S{sid}.nsub
            MH.mh_shape_type(sid,bid,M,lib);
        end
    end

    % scale types
    for sid=1:M.ns
        for bid=1:M.S{sid}.nsub
            MH.mh_scale_type(sid,bid,M,lib);
        end
    end

    % specific relation parameters
    for sid=1:M.ns
       if strcmp(M.S{sid}.R.type,'unihist')
          MH.mh_gobal_position(sid,M,lib);
       elseif strcmp(M.S{sid}.R.type,'mid')
          MH.mh_eval_spot_type(sid,M,lib);
       end
    end

    % sub-stroke ids
    for sid=1:M.ns
        for bid=1:M.S{sid}.nsub
             MH.gibbs_substroke_id(sid,bid,M,lib);
        end
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
% vgg_line3d_from_lP_nonlin  Non-linear estimation of (possibly constrained) 3D line segment from image line segments.
%
% SYNOPSIS
% L = vgg_line3d_from_lP_nonlin(s,P [,imsize] [,L0] [,X] [,nonlin_opt])
%
% s ... cell(K) of double(3,3), inv. covariance matrices of the K image line segments:-
%   - If the segments are estimated from edges, it is s(:,k) = x*x',
%     where x (3-by-N) are homog. coordinates of the edgels with last components 1.
%   - If only end points are available, s(:,k) = d*x*y' where x, y (column 2-vectors)
%     are the segment's end points and d its length.
%
% P ... K-cell with 3-by-4 camera matrices
%
% imsize ... size (2,K), image size(s) for preconditiong.
%   Omit if s and P are already preconditioned.
%
% L0 ... double(4,2), initial scene line (optional). Homogeneous points L0(:,i) span
%   the line. If omitted, linear estimation is done first.
%
% X ... constraint on L :-
%   - if X is omitted:     no constraint on L
%   - if X is double(4,1): L goes thru point X
%   - if X is double(4,2): L goes thru 3D line spanned by X(:,i)
%
% nonlin_opt ... options for Levenberg-Marquardt. It is comma-separated list
%   of pairs ['option',value]. Possible options are :-
% opt ... options structure with possible fields :-
%   - verbose ... 1 or 0 (default: 0)
%   - niter_term ... maximum number of iterations
%   - rmsstep_term ... terminating step of rms of residuals
%   - lambda_term ... terminating value of lambda (default: 1e10)
%   - lambda_init ... initial value of lambda
% E.g., vgg_line3d_from_lP_nonlin(...,'lambda_init',1e-9,'niter_term',5).
%
% L ... double(4,2), estimated 3D line. Points L(:,i) span the line.
%
% Note: use [] if you want to omit a parameter and use a later one, e.g.
%   vgg_line3d_from_lP_nonlin(s,P,imsize,[],[],'verbose',1,'lam_init',1e-9)
%
% ALGORITHM
% - Minimization is done by Levenberg-Marquardt.
% - 3D line L is parameterized by image lines in the first two images.
%   The positions of these image lines are possibly constrained by X.

% T.Werner, Feb 2002

function L = vgg_line3d_from_lP_nonlin(s,P,imsize,L,X,varargin)

if nargin < 3, imsize = []; end
if nargin < 4, L = []; end
if nargin < 5, X = []; end

if isempty(L)
  L = vgg_line3d_from_lP_lin(s,P,imsize);
end
K = length(P); % number of images
if K<2
  error('Cannot reconstruct 3D line from 1 image');
end
if isempty(X) & K==2 % no need for non-linear minimization
  return
end

% Prepare square root of covariance matrices; now s{k}(:,n) has meaning of 3 homogeneous image points
for k = 1:K
  [us,ss,vs] = svd(s{k},0);
  s{k} = us*sqrt(ss);
end

% Preconditioning
if ~isempty(imsize)
  for k = 1:K
    H = vgg_conditioner_from_image(imsize(:,k));
    P{k} = H*P{k};
    s{k} = H*s{k};
    scale(k) = H(1,1); % save the scales for evaluating objective function
  end
else
  scale = ones(1,K);
end

switch size(X,2)
 case 0
  % Scene line L is unconstrained, having thus 4 DOF.
  % L is parameterized by two image lines in images 1 and 2, each having 2 DOF, as follows:
  %   l1 = l0(1,:) + p(1:2)'*ldelta{1}
  %   l2 = l0(2,:) + p(3:4)'*ldelta{2}
  % where row 4-vector p represents 4 DOF of L.
  for k = 1:2
    l0(k,:) = normx(vgg_wedge(P{k}*L)')';
    ldelta{k} = null(l0(k,:))';
  end

  % optimization
  p = levmarq(@F, {vertcat(P{:}),s,scale,l0,ldelta},...
              @normsolve,...
              [0;0;0;0],...
              varargin{:});
  l = l12_from_p(p,l0,ldelta);

 case 1
  % Scene line L is constrained to intersect the scene point X, having thus 2 DOF.
  % L is parameterized by two image lines in images 1 and 2, each having 1 DOF, as follows:
  %   l1 = l0(1,:) + p(1)*ldelta{1}
  %   l2 = l0(2,:) + p(2)*ldelta{2}
  % where 2-vector p represents 2 DOF of L.
  for k = 1:2
    l0(k,:) = normx(vgg_wedge(P{k}*L)')';
    x = P{k}*X;

    % Since L might not intersect X, move l0(k,:) 'as little as possible' to intersect x.
    l0(k,:) = l0(k,:) - (l0(k,:)*x)/(x'*x).*x';

    Q = null(x')';
    ldelta{k} = null(l0(k,:)*pinv(Q))'*Q;
  end

  % optimization
  p = levmarq(@F, {vertcat(P{:}),s,scale,l0,ldelta},...
              @normsolve,...
              [0;0],...
              varargin{:});
  l = l12_from_p(p,l0,ldelta);

 case 2
  % Scene line L is constrained to intersect the scene line given by 2 points X.
  % This constraint is given by
  %   l1*G*l2' = 0
  % where G is 3x3 rank 2 matrix (analogical in fact to fundamental matrix) and
  %   G = P{1}*M*P{2}'
  % where M is Pluecker matrix of line given by X, M = X(:,1)*X(:,2)'-X(:,2)*X(:,1).
  %
  % This constraint can be written as
  %   p(1:2)*D*p(3:4)' + p(1:2)*d{2}' + d{1}*p(3:4)' + c = 0
  % where D, d, c are given below and 4-vector p are 4 parameters of L.
  %
  % L is parameterized by two image lines in images 1 and 2, each having 2 DOF, as follows:
  %   l1 = l0(1,:) + p(1:2)'*ldelta{1}
  %   l2 = l0(2,:) + p(3:4)'*ldelta{2}
  % where p(1:3) are chosen freely and p(4) is computed from the above formula as
  %   p(4) = -(p(1:2)'*(D(:,1)*p(3)+d{1})+p(3)*d{2}(1)+c)/(p(1:2)'*D(:,2)+d{2}(2)).
  Lpm = X(:,1)*X(:,2)' - X(:,2)*X(:,1)';
  G = P{1}*Lpm*P{2}';
  for k = 1:2
    l0(k,:) = normx(vgg_wedge(P{k}*L)')';
  end

  % As L might not intersect line X, move l0(2,:) 'as little as possible' to enforce l0(1,:)*G*l0(2,:)'==0.
  x = (l0(1,:)*G)';
  l0(2,:) = l0(2,:) - (l0(2,:)*x)/(x'*x).*x';

  for k = 1:2
    ldelta{k} = null(l0(k,:))';
  end

  D = ldelta{1}*G*ldelta{2}';
  d{1} = ldelta{1}*G*l0(2,:)';
  d{2} = ldelta{2}*G'*l0(1,:)';
  c = l0(1,:)*G*l0(2,:)';

  % optimization
  p = levmarq(@F, {vertcat(P{:}),s,scale,l0,ldelta,D,d,c},...
              @normsolve,...
              [0;0;0],...
              varargin{:});
  l = l12_from_p(p,l0,ldelta,D,d,c);

end
if all(~isnan(l(:)))
  L = null([l(1,:)*P{1}; l(2,:)*P{2}]);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% objective function

% Objective function of Levenberg-Marquardt
function [y,w,J] = F(p,P,s,scale,varargin)
K = length(s);

% l := lines in images 1, 2 from p
l = l12_from_p(p,varargin{:});

% l := reprojection of lines in images 1, 2 to all images
if all(abs(p) < inf)
  [dummy,dummy,L]= svd([l(1,:)*P(1:3,:); l(2,:)*P(4:6,:)],0);
else
  L = inf*ones(4);
end
l = norml(vgg_wedge(reshape(P*L(:,3),[3 K]),reshape(P*L(:,4),[3 K])));

% compute residual function
y = [];
for k = 1:K
  y = [y l(k,:)*s{k}];
end
y = y';

w = [1;1;1] * (1./scale);
w = w(:);

% else, compute also jacobian
if nargout < 2
  return
end
dif = 1e-6;
J = zeros(length(y),length(p));
for i = 1:length(p)
  pdif = p;
  pdif(i) = pdif(i) + dif;
  J(:,i) = (F(pdif,P,s,scale,varargin{:}) - y)/dif;
end
return

% The following function computes lines in the first two images
% from parameters p. Explanation see above.
function l = l12_from_p(p,l0,ldelta,D,d,c)
switch length(p)
 case 4 % unconstrained
  l = [l0(1,:) + p(1:2)'*ldelta{1}
       l0(2,:) + p(3:4)'*ldelta{2}];
 case 2 % going thru X
  l = [l0(1,:) + p(1)*ldelta{1}
       l0(2,:) + p(2)*ldelta{2}];
 case 3 % going thru L
  p(4) = -(p(1:2)'*(D(:,1)*p(3)+d{1})+p(3)*d{2}(1)+c)/(p(1:2)'*D(:,2)+d{2}(2));
  l = [l0(1,:) + p(1:2)'*ldelta{1}
       l0(2,:) + p(3:4)'*ldelta{2}];
end
return

function dp = normsolve(J,Y,w,lambda)
OLDWARN = warning('off');
dp = -( J'*diag(w)*J + lambda*eye(size(J,2)) ) \ ( J'*(Y.*w) );
warning(OLDWARN);
return

function x = normx(x)
if ~isempty(x)
  x = x./(ones(size(x,1),1)*sqrt(sum(x.*x)));
end

function l = norml(l)
% l = norml(l)  Multiplies hyperplane l by scalar so that for each n, norm(l(1:end-1,n))==1.
l = l./(sqrt(sum(l(:,1:end-1).^2,2))*ones(1,size(l,2)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a = levmarq(@RES,PARAMS,@NORMSOLVE,a [,opt])  Non-linear least-squares by Levenberg-Marquardt.
%
% Minimizes f(a)'*W*f(a) over a.
%
% @RES ... residual function f called like [e,w,J] = RES(a,PARAMS{:}), where
%    - a ... double(M,1), parameter vector
%    - e ... double(N,1), residual vector
%    - J ... double(N,M), derivative of e wrt a
%    - w ... double(N,1), weights of e; covariance matrix of e is diag(1/e.^2).
%            Use 1 instead of ones(N,1).
% For efficiency, RES should not compute the jacobian if called with two output parameters only.
%
% @NORMSOLVE ... function solving normal equations, called like da = NORMSOLVE(J,e,W,lambda).
% a ... initial parameter vector
% opt ... options structure, see code

function [a,w] = levmarq(RES,PARAMS,NORMSOLVE,a0,varargin)

% options
[opt,rem_opt] = vgg_argparse( { 'niter_term',     +inf,...
                                'drmsrel_term',    0,...
                                'loglambda_term',  6,...
                                'loglambda_init', -4,...
                                'verbose',         0 },...
                              varargin );
if ~isempty(rem_opt), if ~isempty(fieldnames(rem_opt))
  error(['Unknown option(s) ' fieldnames(rem_opt)]);
end, end

% Initial statistics
if opt.verbose
  [e0,w0] = feval(RES,a0,PARAMS{:});
  ssd0 = sum( (e0.*w0).^2 );
  fprintf( '                         [rms=%14.12g] [maxabs=%14.12g]\n',...
          sqrt(ssd0/length(e0)),...
          max(abs(e0.*w0)) );
end

loglambda = opt.loglambda_init;
niter = 0;
while 1

  % Compute actual residual and jacobian
  [e0,w0,J] = feval(RES,a0,PARAMS{:});

  % Update a as a := a0 + da, by finding
  % optimal lambda and solving normal equations for da.
  nfail = 1;
  while (loglambda < opt.loglambda_term)
    niter = niter + 1;

    a = a0 + feval(NORMSOLVE,J,e0,w0,10^loglambda);
    [e,w] = feval(RES,a,PARAMS{:});

    if sum((e.*w).^2) < sum((e0.*w0).^2) % success
      a0 = a;
      loglambda = loglambda - 1;
      break
    end

    if opt.verbose
      fprintf('%4i.%.2i: [loglambda=%3i] [REJECTED]\n',niter,nfail,loglambda);
    end

    loglambda = loglambda + 1;
    nfail = nfail + 1;
  end

  % Print statistic after successful iteration
  ssd0 = sum( (e0.*w0).^2 );
  ssd = sum( (e.*w).^2 );
  if opt.verbose
    fprintf( '%4i   : [loglambda=%3i] [rms=%14.12g] [maxabs=%14.12g] [drmsrel=%4g%%]\n',...
             niter,...
             round(loglambda),...
             sqrt(ssd/length(e)),...
             max(abs(e.*w)),...
             100*(1-sqrt(ssd/ssd0)) );
  end

  % Termination criteria
  test(1) = loglambda <  opt.loglambda_term;
  test(2) = ssd0-ssd  >= opt.drmsrel_term^2*ssd0;
  test(3) = niter     <  opt.niter_term;
  if any(test==0)
    break
  end
end

if opt.verbose
  onoff = {'YES','no'};
  fprintf( ' Levenberg-Marquardt finished succesfully.\n Reason for termination:\n' );
  fprintf( '   lambda  = %s\n', onoff{test(1)+1} );
  fprintf( '   drmsrel = %s\n', onoff{test(2)+1} );
  fprintf( '   niter   = %s\n', onoff{test(3)+1} );
end

return

function print_statistics(niter,loglambda,e0,w0,e,w,opt)
if opt.verbose
  ssd0 = sum( (e0.*w0).^2 );
  ssd = sum( (e.*w).^2 );
  fprintf( '%4i   : [loglambda=%3i] [rms=%14.12g] [maxabs=%14.12g] [drmsrel=%11.5g]\n',...
           niter,...
           round(loglambda),...
           sqrt(ssd/length(e)),...
           max(abs(e.*w)),...
           sqrt(1-ssd/ssd0) );
end
return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function monotonicity = DeriveMonotonicityFromStationary(properties,xL,xU)

monotonicity = 'none';
if isequal(properties.convexity,'convex')
    if xU <= properties.stationary(1)
        monotonicity = 'decreasing';
    elseif xL >= properties.stationary(1)
        monotonicity = 'increasing';
    end
elseif isequal(properties.convexity,'concave')
    if xU <= properties.stationary(1)
        monotonicity = 'increasing';
    elseif xL >= properties.stationary(1)
        monotonicity = 'decreasing';
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function obj = subobjective(weight, ind, idealpoint, method)
%SUBOBJECTIVE function evaluate a point's objective with a given method of
%decomposition.

%   Two method are implemented by far is Weighted-Sum and Tchebesheff.
%   weight: is the decomposition weight.(column wise vector).
%   ind: is the individual point(column wise vector).
%   idealpoint: the idealpoint for Tchebesheff decomposition.
%   method: is the decomposition method, the default is 'te' when is
%   omitted.
%
%   weight and ind can also be matrix. in which have two scenairos:
%   When weight is a matrix, then it's treated as a column wise set of
%   weights. in that case, if ind is a size 1 column vector, then the
%   subobjective is computed with every weight and the ind; if ind is also
%   a matrix of the same size as weight, then the subobjective is computed
%   in a column-to-column, with each column of weight computed against the
%   corresponding column of ind.
%   A row vector of subobjective is return in both case.

    if (nargin==2)
        obj = ws(weight, ind);
    elseif (nargin==3)
        obj = te(weight, ind, idealpoint);
    else
        if strcmp(method, 'ws')
            obj=ws(weight, ind);
        elseif strcmp(method, 'te')
            obj=te(weight, ind, idealpoint);
        else
            obj= te(weight, ind, idealpoint);
        end
    end
end

function obj = ws(weight, ind)
    if size(ind, 2) == 1
       obj = (weight'*ind)';
    else
       obj = sum(weight.*ind);
    end
end

function obj = te(weight, ind, idealpoint)
    s = size(weight, 2);
    indsize = size(ind,2);

    weight((weight == 0))=0.00001;

    if indsize==s
        part2 = abs(ind-idealpoint(:,ones(1, indsize)));
        obj = max(weight.*part2);
    elseif indsize ==1
        part2 = abs(ind-idealpoint);
        obj = max(weight.*part2(:,ones(1, s)));
    else
        error('individual size must be same as weight size, or equals 1');
    end
end

author: sfvsfv

repo: Mathematical-modeling


text:
% Generate some shark data
animGT=generateToyAnimation( 4.2,'noise', 0.0, 'fixT', true );
dimSSpace = 2;
anim = computeNRSFM( 2, animGT.W);
[ err ] = computeNRSFMError( animGT, anim );
fprintf( 'Error with Xiao-Kanade: %e \n', sum(sum(err(1:2,:))) );
fprintf( 'Reprojection error: %e \n', computeSFMError( 'reproj', 0, 'anim', anim ) );

% Force some constraints
anim.SBasis(:,:,1)=-anim.SBasis(:,:,1);

anim.l(1,:)=1; anim.t(:)=0;
anim=anim.setFirstPRtToId();
anim = bundleAdjustment( anim, 'nItr', 100 );
anim=anim.setFirstPRtToId();

[ err ] = computeNRSFMError( animGT, anim );
fprintf( 'Error with Xiao-Kanade: %e \n', sum(sum(err(1:2,:))) );
fprintf( 'Reprojection error: %e \n', computeSFMError( 'reproj', 0, 'anim', anim ) );

% Set back to the original camera
for t=1:anim.nFrame
  anim.R(:,:,t)=animGT.R(:,:,t)'*anim.R(:,:,t);
end

animGT=generateToyAnimation( 4,'noise', 0.0, 'fixT', true );

anim.W=animGT.W;

[ err ] = computeNRSFMError( animGT, anim );
fprintf( 'Error with Xiao-Kanade: %e \n', sum(sum(err(1:2,:))) );
fprintf( 'Reprojection error: %e \n', computeSFMError( 'reproj', 0, 'anim', anim ) );

% Normalize the whole thingy
%  anim.SBasis(:,:,i)/mean(anim.l(1,:)); anim.l(1,:)=1;
for i=2:3
  anim.l(i,:)=norm(anim.SBasis(:,:,i),'fro')*anim.l(i,:);
  anim.SBasis(:,:,i)=anim.SBasis(:,:,i)/norm(anim.SBasis(:,:,i),'fro');
  anim.SBasis(:,:,1)=anim.SBasis(:,:,1)+mean(anim.l(i,:))*...
    anim.SBasis(:,:,i);
  anim.l(i,:)=anim.l(i,:)-mean(anim.l(i,:));
end

save('jawSource.mat','anim');

author: vrabaud


text:
function T_old = ReInitializeSolution(T_new,new_Temperature,x_intervals,y_intervals)

% For temperatures at points in the domain which are not at the boundary
for x_index = 2:1:x_intervals-1
    for y_index = 2:1:y_intervals-1

        T_old(x_index,y_index) = T_new(x_index,y_index);

    end
end

% for points at the boundary
% Left BC
for x_index = 1
    for y_index = 1:1:y_intervals
        T_old(x_index,y_index) =  new_Temperature;
    end
end

% Right BC
for x_index = x_intervals
    for y_index = 1:1:y_intervals
        T_old(x_index,y_index) =  new_Temperature;
    end
end

% Bottom BC
for y_index = 1
    for x_index = 1:1:x_intervals
        T_old(x_index,y_index) =  new_Temperature;
    end
end

% Top BC
for y_index = y_intervals
    for x_index = 1:1:x_intervals
        T_old(x_index,y_index) =  new_Temperature;
    end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function test2 (nmat)
%TEST2 test for BTF
% Requires CSparse and UFget
% Example:
%   test2
% See also btf, maxtrans, strongcomp, dmperm, UFget,
%   test1, test2, test3, test4, test5.

% Copyright 2007, Timothy A. Davis, University of Florida

index = UFget ;
f = find (index.nrows == index.ncols) ;

% too much time:
skip = [1514 1297 1876 1301] ;
f = setdiff (f, skip) ;

[ignore i] = sort (index.nnz (f)) ;
f = f (i) ;

if (nargin < 1)
    nmat = 1000 ;
end
nmat = min (nmat, length (f)) ;
f = f (1:nmat) ;

T0 = zeros (nmat,1) ;
T1 = zeros (nmat,1) ;
Anz = zeros (nmat,1) ;
figure (1) ;
clf
MN = zeros (nmat, 2) ;
Nzdiag = zeros (nmat,1) ;

% warmup
p = maxtrans (sparse (1)) ;             %#ok
p = btf (sparse (1)) ;                  %#ok
p = cs_dmperm (sparse (1)) ;            %#ok
a = cs_transpose (sparse (1)) ;         %#ok

h = waitbar (0, 'BTF test 2 of 6') ;

try
    for k = 1:nmat

        Prob = UFget (f (k), index) ;
        A = Prob.A ;

        waitbar (k/nmat, h) ;

        Nzdiag (k) = nnz (diag (A)) ;

        [m n] = size (A) ;
        Anz (k) = nnz (A) ;
        MN (k,:) = [m n] ;

        tic
        [p,q,r] = btf (A) ;
        t0 = toc ;
        s0 = sum (q > 0) ;
        T0 (k) = max (1e-9, t0) ;

        tic
        [p2,q2,r2] = cs_dmperm (A) ;
        t1 = toc ;
        s1 = sum (dmperm (A) > 0) ;
        T1 (k) = max (1e-9, t1) ;

        fprintf ('%4d btf %10.6f cs_dmperm %10.6f', f(k), t0, t1) ;
        if (t1 ~= 0)
            fprintf (' rel: %8.4f', t0 / t1) ;
        end
        fprintf ('\n') ;

        if (s0 ~= s1)
            error ('!') ;
        end

        C = A (p, abs (q)) ;
        subplot (1,2,1) ;
        cspy (C) ;
        z = find (q < 0) ;
        zd = nnz (diag (C (z,z))) ;
        if (zd > 0)
            error ('?') ;
        end

        minnz = Anz (1) ;
        maxnz = nnz (A) ;

        subplot (1,2,2) ;
        loglog (Anz (1:k), T0 (1:k) ./ T1 (1:k), ...
            'o', [minnz maxnz], [1 1], 'r-') ;
        drawnow

        clear C A Prob
    end

catch
    % out-of-memory is OK, other errors are not
    disp (lasterr) ;
    if (isempty (strfind (lasterr, 'Out of memory')))
        error (lasterr) ;                                                   %#ok
    else
        fprintf ('test terminated early, but otherwise OK\n') ;
    end
end

close (h) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function h = p02_h ( n, x )

%*****************************************************************************80
%
%% P02_H evaluates the Hessian for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  for i = 1 : 13

    c = - i / 10.0;

    fi =   x(3) * exp ( c * x(1) )  - x(4) * exp ( c * x(2) ) ...
          + x(6) * exp ( c * x(5) ) -        exp ( c ) ...
          + 5.0 * exp ( 10.0 * c ) -  3.0 * exp ( 4.0 * c );

    dfdx1 =     c     * x(3) * exp ( c * x(1) );
    d2fdx11 =   c * c * x(3) * exp ( c * x(1) );
    d2fdx13 =   c     *        exp ( c * x(1) );
    dfdx2 =   - c     * x(4) * exp ( c * x(2) );
    d2fdx22 = - c * c * x(4) * exp ( c * x(2) );
    d2fdx24 = - c     *        exp ( c * x(2) );
    dfdx3 =                    exp ( c * x(1) );
    d2fdx31 =   c     *        exp ( c * x(1) );
    dfdx4 =   -                exp ( c * x(2) );
    d2fdx42 = - c     *        exp ( c * x(2) );
    dfdx5 =     c     * x(6) * exp ( c * x(5) );
    d2fdx55 =   c * c * x(6) * exp ( c * x(5) );
    d2fdx56 =   c     *        exp ( c * x(5) );
    dfdx6 =                    exp ( c * x(5) );
    d2fdx65 =   c     *        exp ( c * x(5) );

    h(1,1) = h(1,1) + 2.0 * dfdx1 * dfdx1 + 2.0 * fi * d2fdx11;
    h(1,2) = h(1,2) + 2.0 * dfdx2 * dfdx1;
    h(1,3) = h(1,3) + 2.0 * dfdx3 * dfdx1 + 2.0 * fi * d2fdx13;
    h(1,4) = h(1,4) + 2.0 * dfdx4 * dfdx1;
    h(1,5) = h(1,5) + 2.0 * dfdx5 * dfdx1;
    h(1,6) = h(1,6) + 2.0 * dfdx6 * dfdx1;

    h(2,1) = h(2,1) + 2.0 * dfdx1 * dfdx2;
    h(2,2) = h(2,2) + 2.0 * dfdx2 * dfdx2 + 2.0 * fi * d2fdx22;
    h(2,3) = h(2,3) + 2.0 * dfdx3 * dfdx2;
    h(2,4) = h(2,4) + 2.0 * dfdx4 * dfdx2 + 2.0 * fi * d2fdx24;
    h(2,5) = h(2,5) + 2.0 * dfdx5 * dfdx2;
    h(2,6) = h(2,6) + 2.0 * dfdx6 * dfdx2;

    h(3,1) = h(3,1) + 2.0 * dfdx1 * dfdx3 + 2.0 * fi * d2fdx31;
    h(3,2) = h(3,2) + 2.0 * dfdx2 * dfdx3;
    h(3,3) = h(3,3) + 2.0 * dfdx3 * dfdx3;
    h(3,4) = h(3,4) + 2.0 * dfdx4 * dfdx3;
    h(3,5) = h(3,5) + 2.0 * dfdx5 * dfdx3;
    h(3,6) = h(3,6) + 2.0 * dfdx6 * dfdx3;

    h(4,1) = h(4,1) + 2.0 * dfdx1 * dfdx4;
    h(4,2) = h(4,2) + 2.0 * dfdx2 * dfdx4 + 2.0 * fi * d2fdx42;
    h(4,3) = h(4,3) + 2.0 * dfdx3 * dfdx4;
    h(4,4) = h(4,4) + 2.0 * dfdx4 * dfdx4;
    h(4,5) = h(4,5) + 2.0 * dfdx5 * dfdx4;
    h(4,6) = h(4,6) + 2.0 * dfdx6 * dfdx4;

    h(5,1) = h(5,1) + 2.0 * dfdx1 * dfdx5;
    h(5,2) = h(5,2) + 2.0 * dfdx2 * dfdx5;
    h(5,3) = h(5,3) + 2.0 * dfdx3 * dfdx5;
    h(5,4) = h(5,4) + 2.0 * dfdx4 * dfdx5;
    h(5,5) = h(5,5) + 2.0 * dfdx5 * dfdx5 + 2.0 * fi * d2fdx55;
    h(5,6) = h(5,6) + 2.0 * dfdx6 * dfdx5 + 2.0 * fi * d2fdx56;

    h(6,1) = h(6,1) + 2.0 * dfdx1 * dfdx6;
    h(6,2) = h(6,2) + 2.0 * dfdx2 * dfdx6;
    h(6,3) = h(6,3) + 2.0 * dfdx3 * dfdx6;
    h(6,4) = h(6,4) + 2.0 * dfdx4 * dfdx6;
    h(6,5) = h(6,5) + 2.0 * dfdx5 * dfdx6 + 2.0 * fi * d2fdx65;
    h(6,6) = h(6,6) + 2.0 * dfdx6 * dfdx6;

  end

  return
end


text:
function  [Truecol,x,y]=m_shadedrelief(x,y,Z,varargin)
% M_SHADEDRELIEF Shaded relief topography in an image
%  M_SHADEDRELIEF(X,Y,Z) presents a shaded relief topography, as would
%  be seen if a 3D model was artifically lit from the side. Slopes
%  facing the light are lightened, and slopes facing away from the
%  light are darkened. X and Y are horizontal and vertical coordinate
%  VECTORS, and these should be in the same units as the height Z
%  (e.g., all in meters), otherwise the slope angle calculations will be
%  in error.
%
%  Usage notes:
%
%  (1) M_SHADEDRELIEF is a replacement for a low-level call to IMAGE
%  displaying a true-colour image so it MUST be preceded by COLORMAP and
%  CAXIS calls.
%
%  (2) M_SHADEDRELIEF probably is most useful as a backdrop to maps with
%  a rectangular outline box - either a cylindrical projection, or some
%  other projection with M_PROJ(...'rectbox','on').
%
%  (3) Finally, the simplest way of not running into problems:
%     - if your elevation data is in LAT/LON coords (i.e. in a matrix where
%       each row has points with the same latitude, and each column has points
%       with the same longitude), use
%                  M_PROJ('equidistant cylindrical',...)
%     - if your elevation data is in UTM coords (meters E/N), i.e. in a matrix
%       where each row has the same UTM northing and the each column has the
%       same UTM easting, use
%                  M_PROJ('utm',....)
%
%  M_SHADEDRELIEF(...,'parameter',value) lets you set various properties.
%  These are:
%       'coords' : Coordinates of X/Y/Z:
%                     'geog' for lat/lon, Z meters,  (default)
%                     'map'  for X/Y map coordinates, Z meters
%                     'Z' if X/Y/Z are all in same units (e.g., meters)
%       'lightangle' : true direction (degrees) of light source (default
%                      -45, i.e. from the north-west)
%       'gradient': Shading effects increase with slope angle
%                   until slopes reach this value (in degrees), and are
%                   held constant for higher slopes (default 10). Reduce
%                   for smoother surfaces.
%       'clipval' : Fractional change in shading for slopes>='gradient'.
%                  0 means no change, 1 means saturation to white or black
%                  if slope is facing directly towards or away from light
%                  source, (default 0.9).
%       'nancol'  : RGB colour of NaN values (default [1 1 1]);
%       'lakecol' : RGB colour of lakes (flat sections) (default NaN)
%                   If set to NaN lakes are ignored.
%
%   IM=M_SHADEDRELIEF(...) returns a handle to the image.
%
%   [SR,X,Y]=m_SHADEDRELIEF(...) does not create an image but only returns
%   the  true-color matrix SR of the shaded relief, as well as the X/Y
%   vectors needed to display it.
%
%   Example:
%           load topo
%           subplot(2,1,1);  % Example without it
%           imagesc(topolonlim,topolatlim,topo);
%           caxis([-5000 5000]);
%           colormap([m_colmap('water',64);m_colmap('gland',64)]);
%           set(gca,'ydir','normal');
%
%           subplot(2,1,2);  % Example with it
%           caxis([-5000 5000]);
%           colormap([m_colmap('water',64);m_colmap('gland',64)]);
%           m_shadedrelief(topolonlim,topolatlim,topo,'gradient',5e2,'coord','Z');
%           axis tight
%

% Rich Pawlowicz (rich@eoas.ubc.ca) Dec/2017
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
% Changes:
% Jan/2018 - changed outputs for flexibility
%            and added 'map' coordinate handling
% Mar/2019 - added alphamapping for out-of-map areas, started using colormap
%            local to AXES not to FIGURE.
% Apr/2019 - some parts relied on the ones-expansion; went back to meshgrid
%            for compatibility with older matlab versions (thanks P. Grahn)

global MAP_PROJECTION MAP_VAR_LIST

lighthead=-45;
gradfac=10;
clipval=.9;
nancol=[1 1 1];
lakecol=NaN; %[.7 .9 1];
geocoords='geog';
scfac=6400000;  % Used for geo coordinates if on sphere radius 1

while ~isempty(varargin)
    switch lower(varargin{1}(1:3))
        case 'coo'
            switch lower(varargin{2}(1))
                case 'g'
                    geocoords='geog';
                case 'm'
                    geocoords='map';
                case {'z','u'}
                    geocoords='z';
                otherwise
                    error('Unknown coordinate specification');
            end
        case 'lig'
            lighthead=varargin{2};
        case 'gra'
            gradfac=varargin{2};
        case 'cli'
            clipval=varargin{2};
        case 'nan'
            nancol=varargin{2};
        case 'lak'
            lakecol=varargin{2};
        otherwise
	        error(['m_shadedrelief: Unknown option: ' varargin{1}]);
     end
     varargin(1:2)=[];
end

% All kinds of issues dealing with coords:
% First, we need VECTOR x/y as an input to gradient function.

if isvector(x) && size(x,2)==1
     x=x';
elseif ~isvector(x)       % Can't be a matrix
    error('Input X must be a VECTOR');
end
if isvector(y) && size(y,1)==1
    y=y';
elseif ~isvector(y)
    error('Input Y must be a VECTOR');
end

%  Now handle the case if we just give a starting and an ending point
if length(x)==2
    x=linspace(0,1,size(Z,2))*diff(x)+x(1);
end

if length(y)==2
    y=linspace(0,1,size(Z,1))'*diff(y)+y(1);
end

if strcmp(geocoords,'geog')                           %  If its Lat/Long points

    % Have to have initialized a map first
    if isempty(MAP_PROJECTION)
      disp('No Map Projection initialized - call M_PROJ first!');
      return;
    end

    % Convert to X/Y
    [X,Y]=m_ll2xy(x(1,:),repmat(mean(y(:,1)),1,size(x,2)));
    [X2,Y2]=m_ll2xy(repmat(mean(x(1,:)),size(y,1),1),y(:,1));
    x=X;
    y=Y2;

end

% Note - 'image' spaces points evenly, so we should just check that they
% are even otherwise the image won't line up with coastlines...

if max(abs( x - linspace(x(1),x(end),length(x)) ) )/abs(x(end)-x(1)) >.005
    warning(['********** Image will be distorted in X direction!! use M_IMAGE to re-map? *************']);
end

if max(abs( y - linspace(y(1),y(end),length(y))' ) )/abs(y(end)-y(1)) >.005
    warning(['********** Image will be distorted in Y direction!! use M_IMAGE to re-map? *************']);
end

% Convert colours to uint8s
if all(nancol<=1)
    nancol=uint8(nancol*255);
end

% Convert colours to uint8s
if all(lakecol<=1)
    lakecol=uint8(lakecol*255);
end

% Get caxis

clims=caxis;
if all(clims==[0 1])   % Not set
    clims=[min(Z(:)) max(Z(:))];
end

% Get colormap for the current axes
cc=colormap(gca);
cc2=round(cc*255);  % we need these in 0-255 range to get Truecolor
lcc=size(cc,1);

%inan=isnan(Z);

% Get slopes
% If we are using a normal ellipsoid we need to rescale
% x/y to get true slope angles
if  (strcmp(geocoords,'map') || strcmp(geocoords,'geog')) && strcmp(MAP_VAR_LIST.ellipsoid,'normal')
   scfac=6370997;
   [Fx,Fy]=gradient( Z, x*scfac, y*scfac);
else
   [Fx,Fy]=gradient( Z, x, y);
end

% Find NaN
[inan,jnan]=find(isnan(Z) | isnan(Fx) | isnan(Fy) );

% Probable lakes
[islake,jlake]=find(Fx==0 & Fy==0);

% Convert z levels into a colormap index.

% Some iteration to discover the exact formula that matlab uses for mapping to
% color indices (from 1 to lcc)
%idx=min( floor( min(max( (Z-clims(1))/(clims(2)-clims(1)),0) ,1 )*lcc )+1,lcc);
idx=max(min( floor(   (Z-clims(1))/(clims(2)-clims(1))*lcc  )+1  ,lcc),1);

% The slope angle relative to the light direction in degrees.
Fnw=atand(imag(-(Fx+i*Fy)*exp(i*lighthead*pi/180)));

%Put an upper and lower limit on the angles
%%Fnw=min(clipval,max(-clipval,Fnw/gradfac));
Fnw=clipval*tanh(Fnw/gradfac);

% Now get the colormap for each pixel and scale the RGB value 'c'.
% If the correction is -0.1 then scale  c*(1 - |-0.1|)
% If the correction is  0.1 then scale  c*(1 - |+0.1|) + 0.1*255
%depending on the slope.

%Truecol=uint8(max(0,min(255,   reshape([cc2(idx,:)],[size(idx) 3]).*repmat(1-abs(Fnw),1,1,3)+repmat(255*Fnw.*(Fnw>0),1,1,3) ) ));
Truecol=uint8( reshape([cc2(idx,:)],[size(idx) 3]).*repmat(1-abs(Fnw),1,1,3)+repmat(255*Fnw.*(Fnw>0),1,1,3) ) ;

%Truecol=uint8(max(0,min(255,   reshape([cc2(idx,:)],[size(idx) 3]).*repmat(1+Fnw/gradfac,1,1,3) ) ));

% Colour Lakes
if any(islake) && isfinite(lakecol(1))
    Truecol(sub2ind(size(Truecol),islake,jlake,  ones(size(islake))))=lakecol(1);
    Truecol(sub2ind(size(Truecol),islake,jlake,1+ones(size(islake))))=lakecol(2);
    Truecol(sub2ind(size(Truecol),islake,jlake,2+ones(size(islake))))=lakecol(3);
end

% Colour the NaNs
if any(inan)
    Truecol(sub2ind(size(Truecol),inan,jnan,  ones(size(inan))))=nancol(1);
    Truecol(sub2ind(size(Truecol),inan,jnan,1+ones(size(inan))))=nancol(2);
    Truecol(sub2ind(size(Truecol),inan,jnan,2+ones(size(inan))))=nancol(3);
end

if strcmp(geocoords,'map')   % Have to "make invisible" the points outside the map limits.

    [xm,ym]=meshgrid(x,y);
    [HLG,HLT]=m_xy2ll(xm,ym);

    % Find pixels outside the limits of the actual map (if the boundary
    % isn't a rectangle)
    if strcmp(MAP_VAR_LIST.rectbox,'off')
        [I,J]=find(HLT<MAP_VAR_LIST.lats(1) | HLT>MAP_VAR_LIST.lats(2) | HLG<MAP_VAR_LIST.longs(1) | HLG>MAP_VAR_LIST.longs(2));
    elseif strcmp(MAP_VAR_LIST.rectbox,'circle')
        R=(xm.^2 +ym.^2);
        [I,J]=find(R>MAP_VAR_LIST.rhomax.^2);
    else
        I=[];J=[];
    end

    backcolor=uint8(get(gcf,'color')*255);

   if any(I)                          % if some pixels are outside the map area
       for k=1:3
           IJ=sub2ind(size(Truecol),I,J,repmat(k,size(I)));
           Truecol(IJ)=backcolor(k);           % Set them to the background colour
       end
   end
else
    I=[];J=[];
end

if nargout<=1
  x_ok = ~isnan(x);
  y_ok = ~isnan(y);
  if any(I)  % make pixels outside the map area transparent, if needed.
     alphadata=  ones(size(Truecol,1),size(Truecol,2),'logical');
     IJ=sub2ind(size(alphadata),I,J);
     alphadata(IJ)=0;
     Truecol=image('xdata',x(x_ok),'ydata',y(y_ok),'cdata', Truecol(y_ok, x_ok,:),'alphadata',alphadata,'tag','m_shadedrelief');
  else
     Truecol=image('xdata',x(x_ok),'ydata',y(y_ok),'cdata', Truecol(y_ok, x_ok,:),'tag','m_shadedrelief');
  end
end

author: goGPS-Project


text:
function fields = TmodelFields()
% Returns the field names of the `Tmodel` structure. This is
% to help alleviate control between different scripts.
% Called by `cleanTmodel`, `readCbTmodel`, `verifyModel`, `organizeModelCool`, `TmodelStats`, `orderModelFieldsBorg`.
%
% USAGE:
%
%    fields = TmodelFields()
%
% OUTPUTS:
%    Fields:       A column cell array that contains the following in order:
%
%                    * rxnFields
%                    * rNumField
%                    * metFields
%                    * mNumFields
%                    * allFields
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

fields = cell(4,1) ; % Declare fields

% Reaction related fields that are cell arrays.
rxnFields = {'rxns' 'rxnID' 'rxnNames' 'subSystems' 'rxnECNumbers' ...
             'rxnKEGGID' 'rxnSEEDID' 'rxnEquations' ...
             'rxnReferences' 'rxnNotes' 'grRules'}' ;
fields{1,1} = rxnFields ;

% Reaction related fields that are numeric arrays.
rNumFields = {'lb' 'ub' 'c'}' ;
fields{2,1} = rNumFields ;

% Metabolite related fields that are cell arrays
metFields = {'mets' 'metID' 'metNames' 'metFormulas' ...
             'metKEGGID' 'metSEEDID' ...
             'metChEBIID' 'metPubChemID' 'metInChIString'}' ;
fields{3,1} = metFields ;

% Metabolite related fields that are numeric arrays.
mNumFields = {'metCharge'} ;
fields{4,1} = mNumFields ;

% All field names in correct order (27 fields).
allFields = {'rxns' 'mets' 'S' 'lb' 'ub' 'c' ...
             'rxnID' 'rxnNames' 'subSystems' 'rxnEquations' ...
             'rxnECNumbers' 'rxnKEGGID' 'rxnSEEDID' ...
             'rxnReferences' 'rxnNotes' 'grRules' ...
             'metID' 'metNames' 'metFormulas' 'metCharge' ...
             'metKEGGID' 'metSEEDID' 'metChEBIID' 'metPubChemID' ...
             'metInChIString' ...
             'genes' 'description'}' ;
fields{5,1} = allFields ;

author: opencobra

repo: cobratoolbox


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 flann_set_distance_type(type, order)
%FLANN_LOAD_INDEX  Loads an index from disk
%
% Marius Muja, March 2009

    distances = struct('euclidean', 1, 'manhattan', 2, 'minkowski', 3, 'max_dist', 4, 'hik', 5, 'hellinger', 6, 'chi_square', 7, 'cs', 7, 'kullback_leibler', 8, 'kl', 8);
    function id = value2id(map,value)
        id = map.(value);
    end

    if ~isnumeric(type),
        type = value2id(distances,type);
    end
    if type~=3
        order = 0;
    end
    nearest_neighbors('set_distance_type', type, order);
end

repo: gspbox


text:
% EUCL: Calculates the euclidean distances among a set of points, or between a
%       reference point and a set of points, or among all possible pairs of two
%       sets of points, in P dimensions.  Returns a single distance for two points.
%
%     Syntax: dists = eucl(crds1,crds2)
%
%        crds1 = [N1 x P] matrix of point coordinates.  If N=1, it is taken to
%                   be the reference point.
%        crds2 = [N2 x P] matrix of point coordinates.  If N=1, it is taken to
%                   be the reference point.
%        -----------------------------------------------------------------------
%        dists = [N1 x N1] symmetric matrix of pairwise distances (if only crds1
%                   is specified);
%                [N1 x 1]  col vector of euclidean distances (if crds1 & ref
%                   are specified);
%                [1 x N2]  row vector of euclidean distances (if ref & crds2
%                   are specified);
%                [N1 x N2] rectangular matrix of pairwise distances (if crds1
%                   & crds2 are specified);
%                [1 x 1]   scalar (if crds1 is a [2 x P] matrix or ref1 & ref2
%                   are specified);
%

% RE Strauss, 5/4/94
%   10/28/95 - output row (rather than column) vector for the (reference
%               point)-(set of points) case; still outputs column vector for the
%               (set of points)-(reference point) case.
%   10/30/95 - for double for-loops, put one matrix-col access in outer loop
%               to increase speed.
%   10/12/96 - vectorize inner loop to increase speed.
%    6/12/98 - allow for P=1.
%   11/11/03 - initialize dists to NaN for error return.

function dists = eucl(crds1,crds2)
  if (~nargin) help eucl; return; end;
  dists = NaN;

  if (nargin < 2)                     % If only crds1 provided,
    [N,P] = size(crds1);
    if (N<2)
      error('  EUCL: need at least two points');
    end;

    crds1 = crds1';                     % Transpose crds
    dists = zeros(N,N);                 % Calculate pairwise distances

    for i = 1:N-1
      c1 = crds1(:,i) * ones(1,N-i);
      if (P>1)
        d = sqrt(sum((c1-crds1(:,(i+1:N))).^2));
      else
        d = abs(c1-crds1(:,(i+1:N)));
      end;
      dists(i,(i+1:N)) = d;
      dists((i+1:N),i) = d';
    end;
    if (N==2)                            % Single distance for two points
      dists = dists(1,2);
    end;

  else                                % If crds1 & crds2 provided,
    [N1,P1] = size(crds1);
    [N2,P2] = size(crds2);
    if (P1~=P2)
      error('  EUCL: sets of coordinates must be of same dimension');
    else
      P = P1;
    end;

    crds1 = crds1';                     % Transpose crds
    crds2 = crds2';

    if (N1>1 & N2>1)                    % If two matrices provided,
      dists = zeros(N1,N2);               % Calc all pairwise distances between them
      for i = 1:N1
        c1 = crds1(:,i) * ones(1,N2);
        if (P>1)
          d = sqrt(sum((c1-crds2).^2));
        else
          d = abs(c1-crds2);
        end;
        dists(i,:) = d;
      end;
    end;

    if (N1==1 & N2==1)                  % If two vectors provided,
      dists = sqrt(sum((crds1-crds2).^2));  % Calc scalar
    end;

    if (N1>1 & N2==1)                   % If matrix & reference point provided,
      crds1 = crds1 - (ones(N1,1)*crds2')'; % Center points on reference point
      if (P>1)                              % Calc euclidean distances in P-space
         dists = sqrt(sum(crds1.^2))';
      else
         dists = abs(crds1)';
      end;
    end;                                    % Return column vector

    if (N1==1 & N2>1)                   % If reference point & matrix provided,
      crds2 = crds2 - (ones(N2,1)*crds1')'; % Center points on reference point
      if (P>1)                              % Calc euclidean distances in P-space
        dists = sqrt(sum(crds2.^2));
      else
        dists = abs(crds2);
      end;
    end;                                    % Return row vector
  end;

  return;

author: PatternRecognition

repo: OpenBMI


text:
function c = times(a,b)
% TIMES implements a.*b, where either a or b is an adiff object.

switch [class(a),class(b)]

case 'adiffdouble'
   c = adiff(b.*a.x, rowmult(b, a.dx), a.root);

case 'doubleadiff'
   c = adiff(a.*b.x, rowmult(a, b.dx), b.root);

case 'adiffadiff'
   checkroot(a,b);
   c = adiff( a.x.*b.x, rowmult(b.x, a.dx)+rowmult(a.x, b.dx), a.root);

otherwise
   error(['Can''t multiply (.*) ',class(a),' and ',class(b)]);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function o = niftistruc
% Create a data structure describing NIFTI headers
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

persistent org;
if ~isempty(org),
    o = org;
    return;
end;
t = struct('conv',{ @char , @uint8 , @int16 , @int32 , @single },...
           'prec',{'uint8', 'uint8', 'int16', 'int32', 'single'},...
           'size',{       1,      1,      2,      4,       4 });
c = t(1);
b = t(2);
s = t(3);
i = t(4);
f = t(5);

table = {...
    i, 1,'sizeof_hdr',348
    c,10,'data_type',[]
    c,18,'db_name',[]
    i, 1,'extents',[]
    s, 1,'session_error',[]
    c, 1,'regular','r'
    b, 1,'dim_info',[]
    s, 8,'dim',[3 0 0 0  1 1 1 1 1]
    f, 1,'intent_p1',0
    f, 1,'intent_p2',0
    f, 1,'intent_p3',0
    s, 1,'intent_code',0
    s, 1,'datatype',2
    s, 1,'bitpix',8
    s, 1,'slice_start',[]
    f, 8,'pixdim',[0 1 1 1]
    f, 1,'vox_offset',0
    f, 1,'scl_slope',1
    f, 1,'scl_inter',0
    s, 1,'slice_end',[]
    b, 1,'slice_code',[]
    b, 1,'xyzt_units',10
    f, 1,'cal_max',[]
    f, 1,'cal_min',[]
    f, 1,'slice_duration',[]
    f, 1,'toffset',[]
    i, 1,'glmax',[]
    i, 1,'glmin',[]
    c,80,'descrip','NIFTI-1 Image'
    c,24,'aux_file',''
    s, 1,'qform_code',0
    s, 1,'sform_code',0
    f, 1,'quatern_b',0
    f, 1,'quatern_c',0
    f, 1,'quatern_d',0
    f, 1,'qoffset_x',0
    f, 1,'qoffset_y',0
    f, 1,'qoffset_z',0
    f, 4,'srow_x',[1 0 0 0]
    f, 4,'srow_y',[0 1 0 0]
    f, 4,'srow_z',[0 0 1 0]
    c,16,'intent_name',''
    c, 4,'magic','ni1'};
org = struct('label',table(:,3),'dtype',table(:,1),'len',table(:,2),...
    'offset',0,'def',table(:,4));
os  = 0;
for j=1:length(org)
    os  = os + org(j).dtype.size*ceil(os/org(j).dtype.size);
    fun = org(j).dtype.conv;
    def = [org(j).def zeros(1,org(j).len-length(org(j).def))];
    org(j).def    = feval(fun,def);
    org(j).offset = os;
    os  = os + org(j).len*org(j).dtype.size;
end;
o = org;
return;


text:
function [ area, radout, side ] = polygon_inrad_data_2d ( n, radin )

%*****************************************************************************80
%
%% POLYGON_INRAD_DATA_2D determines polygonal data from its inner radius in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sides of the polygon.
%    N must be at least 3.
%
%    Input, real RADIN, the inner radius of the polygon, that is,
%    the radius of the largest circle that can be inscribed within
%    the polygon.
%
%    Output, real AREA, the area of the regular polygon.
%
%    Output, real RADOUT, the outer radius of the polygon, that is,
%    the radius of the smallest circle that can be described about
%    the polygon.
%
%    Output, real SIDE, the length of one side of the polygon.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_INRAD_DATA_2D - Fatal error!\n' );
    fprintf ( 1, '  Input value of N must be at least 3\n' );
    fprintf ( 1, '  but your input value was N = %d\n', n );
    error ( 'POLYGON_INRAD_DATA_2D - Fatal error!' );
  end

  angle = pi / n;
  area = n * radin * radin * tan ( angle );
  side = 2.0 * radin * tan ( angle );
  radout = 0.5 * side / sin ( angle );

  return
end


text:
function [pos, tri] = octahedron(varargin)

% OCTAHEDRON
%
% Use as
%   [pos tri] = octahedron;
%
% See also TETRAHEDRON ICOSAHEDRON

pos = [
  0  0 +1
  +1  0  0
  0 +1  0
  -1  0  0
  0 -1  0
  0  0 -1
  ];

tri = [
  1 2 3
  1 3 4
  1 4 5
  1 5 2
  6 3 2
  6 4 3
  6 5 4
  6 2 5
  ];

if nargin>0
  n = varargin{1};
  % perform an n-fold refinement
  for i=1:n
    [pos, tri] = refine(pos, tri, 'banks');
  end
  % scale all vertices to the unit sphere
  pos = pos ./ repmat(sqrt(sum(pos.^2,2)), 1,3);
end


text:
function [aVal,aJacob,aHess,papt]=aPoly(x,numDim)
%%APOLY The drift function for a linear continuous-time motion model of a
%       given order in a specified number of Cartesian dimensions. The
%       order of the linear filter, that is the number of moments of
%       position, does not need to be explicitly specified.
%
%INPUTS: x The xDimXN state vector of N targets in the order of
%          [position;velocity;acceleration;etc] for however many
%          derivatives of position there are.
%   numDim The number of dimensions of the simulation problem. If the
%          numDim parameter is omitted, then numDim=3 (3D motion) is
%          assumed. The dimensionality of the state must be an integer
%          multiple of numDim.
%
%OUTPUTS: aVal The xDimXN time-derivative of the N state vectors under the
%              linear motion model.
%       aJacob The xDimXxDim Jacobian of aVal (it is the same for all x and
%              is not repeated N times). This is such that aJacob(:,k) is
%              the partial derivative of aVal with respect to the kth
%              element of x.
%        aHess The xDimXxDimXxDim hypermatrix such that aHess(:,k1,k2) is
%              the second partial derivative of aVal with respect to
%              elements k1 and k2 of x (all zero in this instance). It is
%              the same for all x.
%         papt The xDimX1 partial derivative of aVal with respect to time
%              (all zero in this instance). It is the same for all x.
%
%The drift function corresponds to the state transition given in
%discrete-time by the function FPolyKal.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    numDim=3;
end

numTar=size(x,2);
xDim=size(x,1);

aVal=[x((numDim+1):(xDim),:);zeros(numDim,numTar)];

if(nargout>1)
    aJacob=[zeros(xDim,numDim),[eye(xDim-numDim,xDim-numDim);zeros(numDim,xDim-numDim)],zeros(xDim,numDim)];

    if(nargout>2)
        aHess=zeros(xDim,xDim,xDim);

        if(nargout>3)
            papt=zeros(xDim,1);
        end
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function inspect_bug1995

% DEPENDENCY ft_voumerealign ft_sourceplot

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/mri/dicom/19112010_JHORSCHIG.MR.FCDC_SEQUENCES_STANDARD_SEQUENCES.0002.0176.2010.11.19.12.08.01.265625.73007255.IMA');
mri = ft_read_mri(filename);

cfg =[];
cfg.method = 'interactive';
ft_volumerealign(cfg, mri);

cfg =[];
cfg.method = 'ortho';
% cfg.interactive = 'yes';
ft_sourceplot(cfg, mri);


text:
function [se_within, stats] = barplot_get_within_ste(dat,varargin)
% :Usage:
% ::
%
%     [se_within, stats] = barplot_get_within_ste(dat)
%
% Compute within-subjects standard errors
%
% :Note: The old version of this used average s.e.'s for contrasts
% of interest, which depend on contrast scaling and are less appropriate
% for use as error bars.
%
% The new version as of 12/10 uses the Loftus & Masson 1994 method and has
% been checked against the data in that paper.
%
% See help below for the L & M data and additional
% code for getting the ANOVA decomposition.
%
% Useful for barplots of conditions when calculating within-subject SEs
%
% :Examples:
% ::
%
%    %% data from Loftus & Masson, Table 2
%
%    mtx = [1 10 13 13 12.00
%    2 6 8 8 7.33
%    3 11 14 14 13.00
%    4 22 23 25 23.33
%    5 16 18 20 18.00
%    6 15 17 17 16.33
%    7 1 1 4 2.00
%    8 12 15 17 14.67
%    9 9 12 12 11.00
%    10 8 9 12 9.67];
%
%    dat = mtx(:, 2:4);  % data from Loftus & Masson, Table 2
%
%    [se_within, stats] = barplot_get_within_ste(dat)
%    fprintf('Within ste: %3.2f, 95%% CI: mean +/- %3.2f\n', se_within, stats.ci);
%
% Extra stuff from ANOVA table
%
% Mean square for condition: Variance of condition means * sample
% size...average squared variance accounted for by condition means
% ::
%
%    [n, k] = size(dat);
%    MS_cond = var( mean(dat) - mean(dat(:)) ) * n;
%
%    MS_subject = var( mean(dat') - mean(dat(:)) ) * k;
%
%    datv = dat(:);
%    MS_total = scale(datv, 1)' * scale(datv, 1);

% notes: 6/3/2018 - tor added wh_omit line to omit all-NaN cases and
% row-wise replacement to make error bars robust to presence of NaNs in
% data.  Previously returned NaN/zero error bars in presence of NaNs

if nargin > 1 && ~isempty(varargin{1})

    covs = varargin{1};

    if size(covs,1) ~= size(dat,1), error('Covs and dat must have same no. of observations.'); end

    disp('Warning: covariates are no longer removed for within-subject SE bars.');

end

% remove cases with all NaN values
% otherwise, these cases could influence calculations and interfere with
% scaling below.
wh_omit = all(isnan(dat), 2); dat(wh_omit, :) = [];

% Replace remaining NaNs with row means
m = nanmean(dat')';
m = repmat(m, 1, size(dat, 2));
wh = isnan(dat);
dat(wh) = m(wh);

% number of subjects and conditions
[n, k] = size(dat);

df_s = n - 1;  % degrees of freedom for subjects
df_c = k - 1;  % degrees of freedom for conditions
df_sxc = df_s * df_c;  % degrees of freedom for subj x condition

% within-subject, within-condition errors
errs = scale(scale(dat, 1)', 1)';

% Sums of squared errors
SS_sxc = errs(:)' * errs(:);

% Mean squared error, accounting for reduced df due to estimating marginal means
MS_sxc  = SS_sxc / df_sxc;

se_within = sqrt(MS_sxc / n);

ci = se_within * tinv(1 - (.05 / 2), df_sxc);

stats.se_within = se_within;
stats.ci = ci;
stats.ci_descrip = '95% confidence interval';
stats.SS_sxc = SS_sxc;
stats.MS_sxc = MS_sxc;

% Extra stuff from ANOVA table
% ---------------------------------------------------------------
% Mean square for condition: Variance of condition means * sample
% size...average squared variance accounted for by condition means
stats.MS_cond = var( mean(dat) - mean(dat(:)) ) * n;

stats.MS_subject = var( mean(dat') - mean(dat(:)) ) * k;

datv = dat(:);
stats.MS_total = scale(datv, 1)' * scale(datv, 1);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [x,label] = spm_samp_mix (mix, N)
% Sample from a Gaussian Mixture PDF
% FORMAT [x,label] = spm_samp_mix (mix, N)
%
% mix   Data structure for mixture model (see spm_mix for info)
% N     Number of samples
%
% x     [N x d] matrix of samples
% label [N x 1] vector of sample labels
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_samp_mix.m 1143 2008-02-07 19:33:33Z spm $

priors = rand(1, N);
cum_prior = 0;
total_samples = 0;
label = zeros(N, 1);
for j=1:mix.m,
    num_samples = sum(priors >= cum_prior & ...
        priors < cum_prior + mix.state(j).prior);
    m=mix.state(j).m;
    C=squeeze(mix.state(j).C);
    x(total_samples+1:total_samples+num_samples, :) = ...
        spm_samp_gauss(m, C, num_samples);
    cum_prior = cum_prior + mix.state(j).prior;
    total_samples = total_samples + num_samples;
    label(total_samples+1:total_samples+num_samples) = j;
end

author: spm


text:
clear all; close all; clear classes; clc;

%% Set flags.
inspect_only = false;

%% Create shapes.
a = 420;  % lattice constant
t = 1;  % slab thickness

% permittivity, r/a, and omega*a/(2*pi*c) are taken from p.234 of John D.
% Joannopoulos? et al., "Photonic Crystals: Molding the Flow of Light," 2nd
% edition so that the frequency lies in a band gap.
eps_diel = 11.4;
r = 0.25*a;  % hole radius
wvlen = a/0.3;  % omega*a/(2*pi*c) = k*a/(2*pi) = a/lambda = 0.3

ad = 25;  % divider for a
td = 10;  % divider for t
dd = 10;  % divider for d = 2*r

mx = 25.5;  % half integer puts domain boundary between cylinders and makes PML works better
my = 7.5;
slab_yn = Box([-mx*a mx*a; -my*a -0.5*a; 0 t], [a/ad, a/ad, t]);
slab_yp = Box([-mx*a mx*a; 0.5*a my*a; 0 t], [a/ad, a/ad, t]);

rod = CircularCylinder(Axis.z, t, [0 0 t/2], r, [2*r/dd, 2*r/dd, t]);

%% Solve the system.
gray = [0.5 0.5 0.5];  % [r g b]
src_loc = 6*a;
[E, H, obj_array, src_array, J] = maxwell_run(...
	'OSC', 1e-9, wvlen, ...
	'DOM', {'vacuum', 'white', 1.0}, [-mx*a mx*a; -my*a my*a; 0 t], [a/ad a/ad t], BC.p, [7*a 2*a 0], 2, 1e-4, ...
	'OBJ', ...
		{'dielectric', gray, eps_diel}, periodize_shape(rod, {[a 0 0], [0 a 0], [0 0 t]}, slab_yn), ...
		{'dielectric', gray, eps_diel}, periodize_shape(rod, {[a 0 0], [0 a 0], [0 0 t]}, slab_yp), ...
	'SRCJ', PointSrc(Axis.z, [src_loc, 0, 0.5]), ...
	inspect_only);

%% Visualize the solution.
if ~inspect_only
	figure;
	clear opts
% 	opts.withgrid = true;
	opts.withobjsrc = false;
	opts.withabs = true;
	opts.withpml = false;
	opts.phase = pi/2;
	figure(1)
	vis2d(E{Axis.z}, Axis.z, 0.5, obj_array, src_array, opts)
	%%
	figure(2)
	vis2d(H{Axis.y}, Axis.z, 0.5, obj_array, src_array, opts)

	%%
	flux_loc = 3*a;
	power_right = powerflux_patch(E, H, Axis.x, src_loc + flux_loc);
	power_left = -powerflux_patch(E, H, Axis.x, src_loc - flux_loc);
	fprintf('power:\n');
	fprintf('right = %s\n', num2str(power_right));
	fprintf('left = %s\n', num2str(power_left));
	fprintf('error = %s%%\n',num2str((power_left-power_right)/power_right*100));

	%%
	Sx = poynting(Axis.x, E{Axis.y}, E{Axis.z}, H{Axis.y}, H{Axis.z}, Axis.y, 0);
	[array, l] = Sx.data_original;
	figure(3)
	plot(l{2}, abs(array))
	mx*a - 10*a

% 	%%
% 	xs = E{Axis.z}.grid3d.l{Axis.x,GT.prim};
% 	ez = NaN(size(xs));
% 	i = 0;
% 	for x = xs
% 		i = i + 1;
% 		ez(i) = E{Axis.z}.value([x 0 t/2]);
% 	end
% 	%%
% 	figure(4);
% 	plot(xs, imag(ez));
%
% 	%%
% 	xs = E{Axis.z}.grid3d.l{Axis.x,GT.prim};
% 	hy = NaN(size(xs));
% 	i = 0;
% 	for x = xs
% 		i = i + 1;
% 		hy(i) = H{Axis.y}.value([x 0 t/2]);
% 	end
% 	%%
% 	figure(5);
% 	plot(xs, real(hy));
%
% 	%%
% % 	sx = real(imag(ez) .* real(hy));
% % 	sx = imag(ez .* conj(hy));
% 	sx = real(ez .* conj(hy));
% 	figure(6);
% 	plot(xs, abs(sx));
%
% 	% Need to verify if the same problem occurs for precisely selected frequency
% 	% for defect waveguide.
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function best_weights = GetBestMixingWeightForHarmonics(standard_deviation_list)
% Best mixing of instantaneous frequencies with common mean
% and built-in constraint (safeguard version designed for harmonics information mixing)
%
% best_weights = GetBestMixingWeightForHarmonics(sdListOrg)
%
% Input argument
%   standard_deviation_list    : standard deviation of each variable
%
% Output argument
%   best_weights   : best set of weights for mixing random variables
%
% Note
%   This mixing consists of heuristics. Re-design is needed.

% Copyright 2016 Google Inc. All Rights Reserved.
% Author: hidekik@google.com (Hideki Kawahara)

original_list_length = length(standard_deviation_list);
safe_list = (1:original_list_length)';
safe_list = safe_list(standard_deviation_list > 0 & ...
  standard_deviation_list < 100 * min(abs(standard_deviation_list)));
sdLsafe_SD_list = standard_deviation_list(safe_list);
n_of_safe_SD = length(sdLsafe_SD_list);
if n_of_safe_SD > 0
  H = ones(n_of_safe_SD-1, n_of_safe_SD-1) * sdLsafe_SD_list(n_of_safe_SD) ^ 2;
  H = H + diag(sdLsafe_SD_list(1:n_of_safe_SD - 1) .^ 2);
  v = ones(n_of_safe_SD - 1, 1) * sdLsafe_SD_list(n_of_safe_SD) ^ 2;
  a = H \ v; % revised 23/May/2016 from inv(H) * v
  if sum(a) > 1
    a = a / sum(a);
  end;
  w = [a; 1 - sum(a)];
  best_weights = zeros(original_list_length, 1);
  best_weights(safe_list) = w;
else
  best_weights = zeros(original_list_length, 1);
  best_weights(1) = 1;
end;
end

author: google


text:
[datapoints, numpoints] = px4_read_binary_file('e4_ekf_A.bin');
roll_ekf = datapoints(1, :);
pitch_ekf = datapoints(2, :);
roll_cf = datapoints(3, :);
pitch_cf = datapoints(4, :);
roll_px4 = datapoints(5, :);
pitch_px4 = datapoints(6, :);
t = datapoints(7, :);
rad2deg = 180/pi;
figure(1)
plot(t, roll_ekf*rad2deg, 'g', t, roll_cf*rad2deg, 'b', t, roll_px4*rad2deg, 'r')
legend('ekf', 'cf', 'px4');
xlabel('time/s');
ylabel('roll/deg')
title('roll-Filter effect comparison')

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function k = rbfard2KernDiagCompute(kern, x)

% RBFARD2KERNDIAGCOMPUTE Compute diagonal of RBFARD2 kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the automatic relevance determination radial basis function kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : rbfard2KernParamInit, kernDiagCompute, kernCreate, rbfard2KernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Michalis K. Titsias, 2009

% KERN

k = repmat(kern.variance, size(x, 1), 1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [img, mmPerPix, vSize, fname, format] = loadVolume(pth)
%Read one of a variety of types of volume anatomies
%
% [img, mmPerPix, vSize, fname, format] = loadVolume([file or directory])
%
% An expanded version of loadVolumeDataAnyFormat (with a slightly
% less-unwieldy name). A generic function for loading MRI volume data.
% Hopefully this will be expanded as more file types are being worked
% with.
%
% The first argument can be the name/path of a single file (for formats
% like Analyze .img, mrVista .dat, etc containing and MRI volume, the first
% file in a directory, for formats like DICOM or Genesis I-files in which
% the volume is spread across several files (e.g. Ifiles/I.001), or else
% the directory containing these files (e.g.'Ifiles). If omitted, a dialog
% pops up.
%
% Supported formats are:
%   * mrVista .dat (e.g. vAntomy.dat)
%   * ANALYZE .img or .hdr
%   * Gensis I-file (I.###)
%   * DICOM   .dcm
%   * Freesurfer/FS-FAST .bshort
%   * Freesurfer/FS-FAST .bfloat
%   * Lucas P*.mag fuctional files
%
% Examples:
%    [img, mmPerPix, vSize, fname, format] = loadVolume;
%
% 02/14/05 ras.

img = []; mmPerPix = []; vSize = []; fname = []; format = [];

% if pth is empty, have the user select the proper file.
if ieNotDefined('pth')
    pth = mrvSelectFile('r');
    if isempty(pth), return;
    elseif ~exist(pth,'file'), error(sprintf('%s not found.\n',pth));
    end
else
    % If the user sent in pth, then it might be a directory containing I-files
    if exist(pth,'dir')
        ifileCheck = dir(fullfile(pth,'I*'));
        if ~isempty(ifileCheck)
            pth = fullfile(pth,ifileCheck(1).name);
        end
    end
end

[p,fname,ext] = fileparts(pth);

switch lower(ext)
    case '.dat',
        % vanat .dat
        format = 'vanat';
        [img, mmPerPix, vSize] = readVolAnat(pth);
    case {'.hdr','.img'},
        % ANALYZE
        format = 'analyze';
        [img, mmPerPix, hdr] = analyzeRead(pth);
		% we'll go ahead and reorient this to I|P|R orientation from R|A|S
		img = mrAnatRotateAnalyze(img);
		mmPerPix = mmPerPix([3 2 1]);
        vSize = size(img);
	case {'.nii', '.gz'},
		% NIFTI (or compressed NIFTI)
		format = 'nifti';
		mr = mrLoad(pth);
		% we'll go ahead and reorient this to I|P|R orientation from R|A|S
		img = mrAnatRotateAnalyze(mr.data);
		mmPerPix = mr.voxelSize([3 2 1]);
		vSize = size(img);
		fname = mr.path;
    case '.bshort',
        % FS-FAST .bshort
        format = 'bshort';
        img = loopOverAnatFiles(p,'*.bshort','readBshort');
    case '.bfloat',
        % FS-FAST .bfloat
        format = 'bshort';
        img = loopOverAnatFiles(p,'*.bshort','readBfloat');
    case '.mag',
        % P*.mag functional file
        format = 'pmag';
        img = readMagFile(pth);
    case '.dcm',
        % DICOM
        format = 'dicom';
        img = loopOverAnatFiles(p,'*.dcm','ReadMRImage');
    otherwise,
        % check if it's a GENESIS I-file: extension should be numeric
        if ~isempty(str2num(ext(2:end)))
            % numeric extension -- try to load I-files
            format = 'ifile';
            img = loopOverAnatFiles(p,'I*','ReadMRImage');
        else
            error('Unknown file format!');
        end
end

return
% /-------------------------------------------------------------------/ %

% /-------------------------------------------------------------------/ %
function img = loopOverAnatFiles(parent,pattern,func);
% img = loopOverAnatFiles(parent,pattern,func);
% Within the specified parent directory, find all files
% of the specified pattern and load them with the specified
% function.
img = [];
if ~exist(parent,'dir')
    warning(sprintf('Directory %s not found.',parent));
    return
end
flag4D = 0;
callingDir = pwd;
cd(parent);
w = dir(pattern);
fnames = {w.name};
cd(callingDir);
fprintf('Looping across files');
for i = 1:length(fnames)
    filepath = fullfile(parent,fnames{i});
    subvol = eval(sprintf('%s(''%s'');',func,filepath));
    if size(subvol,3) <= 1 % 2D image
        img(:,:,i) = subvol;
    else
        % try to build up a 4D matrix --
        % assume 3rd dimension is time
        img(:,:,:,i) = subvol;
        flag4D = 1; % will permute later
    end
    fprintf('.')
end
fprintf('done.\n');
if flag4D==1
    img = permute(img,[1 2 4 3]);
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [newLine, isOrthogonal] = projLineOnPlane(line, plane)
%PROJLINEONPLANE Return the orthogonal projection of a line on a plane
%
%   NEWLINE = PROJLINEONPLANE(LINE, PLANE) Returns the orthogonal
%   projection of LINE or multiple lines on the PLANE.
%
%   [..., ISORTHOGONAL] = PROJLINEONPLANE(LINE, PLANE) Also returns if the
%   LINE is orthogonal to the PLANE.
%
%   Example
%     plane = [.1 .2 .3 .4 .5 .6 .7 .8 .9];
%     lines = [0 .3 0 1 0 0;0  .5 .5 0 0 1;...
%         .4 .1 .5 1 0 2;.2 .7 .1 0 1 0;...
%         plane(1:3) planeNormal(plane)];
%     [newLines, isOrthogonal] = projLineOnPlane(lines, plane);
%     figure('color','w'); axis equal; view(3)
%     drawLine3d(lines,'b')
%     drawPlane3d(plane)
%     drawLine3d(newLines(~isOrthogonal,:), 'r')
%
%   See also:
%   planes3d, lines3d, intersectLinePlane, projPointOnPlane
%
% ---------
% Author: oqilipo
% Created: 2017-08-06
% Copyright 2017

p1 = projPointOnPlane(line(:,1:3), plane);
p2 = projPointOnPlane(line(:,1:3)+line(:,4:6), plane);

newLine=createLine3d(p1, p2);
isOrthogonal = ismembertol(p1,p2,'ByRows',true);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function rnrb = nrbreverse(nrb)
%
% Function Name:
%
%   nrbreverse - Reverse the evaluation direction of a NURBS curve or surface.
%
% Calling Sequence:
%
%   rnrb = nrbreverse(nrb);
%
% Parameters:
%
%   nrb		: NURBS data structure, see nrbmak.
%
%   rnrb		: Reversed NURBS.
%
% Description:
%
%   Utility function to reverse the evaluation direction of a NURBS
%   curve or surface.

%  D.M. Spink
%  Copyright (c) 2000.

if nargin ~= 1
  error('Incorrect number of input arguments');
end

if iscell(nrb.knots)

  % reverse a NURBS surface
  coefs = nrb.coefs(:,:,end:-1:1);
  rnrb = nrbmak(coefs(:,end:-1:1,:), {1.0-fliplr(nrb.knots{1}),...
                1.0-fliplr(nrb.knots{2})});

else

  % reverse a NURBS curve
  rnrb = nrbmak(fliplr(nrb.coefs), 1.0-fliplr(nrb.knots));

end


text:
function cls = spm_mars_newSeg(res,tc,bf,df,mode,tcm)
% FORMAT cls = spm_mars_newSeg(res,tc,bf,df,mode,tcm)
%
% Write out the results from toolbox 'New Segment'. This function is only for
% the purpose of making toolbox 'MARS' backwards compatible with 'New Segment'.
%
% Adapted from spm_preproc_write8 by John Ashburner
% $Id: spm_preproc_write8.m 4337 2011-05-31 16:59:44Z john $
% Copyright (C) 2008 Wellcome Department of Imaging Neuroscience
%
% Modified by Yu (Andy) Huang
% $Id: spm_mars_newSeg.m 2015-07-27 andy$
% Neural Engineering Lab, Dept. of Biomedical Engineering, City College of New York
% yhuang16@citymail.cuny.edu

% Read essentials from tpm (it will be cleared later)
tpm = res.tpm;
if ~isstruct(tpm) || ~isfield(tpm, 'bg1'),
    tpm = spm_load_priors8(tpm);
end
d1        = size(tpm.dat{1});
d1        = d1(1:3);
M1        = tpm.M;
[bb1 vx1] = spm_get_bbox(tpm.V(1), 'old');

if isfield(res,'mg'),
    lkp = res.lkp;
    Kb  = max(lkp);
else
    Kb  = size(res.intensity(1).lik,2);
end

N   = numel(res.image);
if nargin<2, tc = true(Kb,4); end % native, import, warped, warped-mod
if nargin<3, bf = true(N,2);  end % field, corrected
if nargin<4, df = true(1,2);  end % inverse, forward
% if nargin<5, mrf= 2;          end % MRF parameter  % andy 2013-05-03

[pth,nam]=fileparts(res.image(1).fname);
ind  = res.image(1).n;
d    = res.image(1).dim(1:3);

[x1,x2,o] = ndgrid(1:d(1),1:d(2),1);
x3  = 1:d(3);

chan(N) = struct('B1',[],'B2',[],'B3',[],'T',[],'Nc',[],'Nf',[],'ind',[]);
for n=1:N,
    d3         = [size(res.Tbias{n}) 1];
    chan(n).B3 = spm_dctmtx(d(3),d3(3),x3);
    chan(n).B2 = spm_dctmtx(d(2),d3(2),x2(1,:)');
    chan(n).B1 = spm_dctmtx(d(1),d3(1),x1(:,1));
    chan(n).T  = res.Tbias{n};

    [pth1,nam1,ext1] = fileparts(res.image(n).fname);
    chan(n).ind      = res.image(n).n;

    if bf(n,2),
        chan(n).Nc      = nifti;
        chan(n).Nc.dat  = file_array(fullfile(pth1,['m', nam1, '.nii']),...
                                     res.image(n).dim(1:3),...
                                     [spm_type('float32') spm_platform('bigend')],...
                                     0,1,0);
        chan(n).Nc.mat  = res.image(n).mat;
        chan(n).Nc.mat0 = res.image(n).mat;
        chan(n).Nc.descrip = 'Bias corrected';
        create(chan(n).Nc);
    end

    if bf(n,1),
        chan(n).Nf      = nifti;
        chan(n).Nf.dat  = file_array(fullfile(pth1,['BiasField_', nam1, '.nii']),...
                                     res.image(n).dim(1:3),...
                                     [spm_type('float32') spm_platform('bigend')],...
                                     0,1,0);
        chan(n).Nf.mat  = res.image(n).mat;
        chan(n).Nf.mat0 = res.image(n).mat;
        chan(n).Nf.descrip = 'Estimated Bias Field';
        create(chan(n).Nf);
    end
end

do_cls   = any(tc(:)) || nargout>1;
tiss(Kb) = struct('Nt',[]);
for k1=1:Kb,
    if tc(k1,4) || any(tc(:,3)) || tc(k1,2) || nargout>=1,
        do_cls  = true;
    end
    if tc(k1,1),
        tiss(k1).Nt      = nifti;
        tiss(k1).Nt.dat  = file_array(fullfile(pth,['c', num2str(k1), nam, '.nii']),...
                                      res.image(n).dim(1:3),...
                                      [spm_type('uint8') spm_platform('bigend')],...
                                      0,1/255,0);
        tiss(k1).Nt.mat  = res.image(n).mat;
        tiss(k1).Nt.mat0 = res.image(n).mat;
        tiss(k1).Nt.descrip = ['Tissue class ' num2str(k1)];
        create(tiss(k1).Nt);
        do_cls = true;
    end;
end

prm     = [3 3 3 0 0 0];
Coef    = cell(1,3);
Coef{1} = spm_bsplinc(res.Twarp(:,:,:,1),prm);
Coef{2} = spm_bsplinc(res.Twarp(:,:,:,2),prm);
Coef{3} = spm_bsplinc(res.Twarp(:,:,:,3),prm);

do_defs = any(df);
do_defs = do_cls | do_defs;
if do_defs,
    if df(1),
        [pth,nam,ext1]=fileparts(res.image(1).fname);
        Ndef      = nifti;
        Ndef.dat  = file_array(fullfile(pth,['iy_', nam1, '.nii']),...
                               [res.image(1).dim(1:3),1,3],...
                               [spm_type('float32') spm_platform('bigend')],...
                               0,1,0);
        Ndef.mat  = res.image(1).mat;
        Ndef.mat0 = res.image(1).mat;
        Ndef.descrip = 'Inverse Deformation';
        create(Ndef);
    end
    if df(2) || any(any(tc(:,[2,3,4]))) || nargout>=1,
        y = zeros([res.image(1).dim(1:3),3],'single');
    end
end

spm_progress_bar('init',length(x3),['Working on ' nam],'Planes completed');
M = M1\res.Affine*res.image(1).mat;

if do_cls
    Q = zeros([d(1:3),Kb],'single');
end

for z=1:length(x3),

    % Bias corrected image
    cr = cell(1,N);
    for n=1:N,
        f          = spm_sample_vol(res.image(n),x1,x2,o*x3(z),0);
        bf         = exp(transf(chan(n).B1,chan(n).B2,chan(n).B3(z,:),chan(n).T));
        cr{n}      = bf.*f;
        if ~isempty(chan(n).Nc),
            % Write a plane of bias corrected data
            chan(n).Nc.dat(:,:,z,chan(n).ind(1),chan(n).ind(2)) = cr{n};
        end;
        if ~isempty(chan(n).Nf),
            % Write a plane of bias field
            chan(n).Nf.dat(:,:,z,chan(n).ind(1),chan(n).ind(2)) = bf;
        end;
    end

    if do_defs,
        [t1,t2,t3] = defs(Coef,z,res.MT,prm,x1,x2,x3,M);
        if exist('Ndef','var'),
            tmp = M1(1,1)*t1 + M1(1,2)*t2 + M1(1,3)*t3 + M1(1,4);
            Ndef.dat(:,:,z,1,1) = tmp;
            tmp = M1(2,1)*t1 + M1(2,2)*t2 + M1(2,3)*t3 + M1(2,4);
            Ndef.dat(:,:,z,1,2) = tmp;
            tmp = M1(3,1)*t1 + M1(3,2)*t2 + M1(3,3)*t3 + M1(3,4);
            Ndef.dat(:,:,z,1,3) = tmp;
        end

        if exist('y','var'),
            y(:,:,z,1) = t1;
            y(:,:,z,2) = t2;
            y(:,:,z,3) = t3;
        end

        if do_cls,
            msk = (f==0) | ~isfinite(f);

            if isfield(res,'mg'),
                q   = zeros([d(1:2) Kb]);
                q1  = likelihoods(cr,[],res.mg,res.mn,res.vr);
                q1  = reshape(q1,[d(1:2),numel(res.mg)]);
                b   = spm_sample_priors8(tpm,t1,t2,t3);
                for k1=1:Kb,
                    q(:,:,k1) = sum(q1(:,:,lkp==k1),3).*b{k1};
                end
            else
                q   = spm_sample_priors8(tpm,t1,t2,t3);
                q   = cat(3,q{:});
                for n=1:N,
                    tmp = round(cr{n}*res.intensity(n).interscal(2) + res.intensity(n).interscal(1));
                    tmp = min(max(tmp,1),size(res.intensity(n).lik,1));
                    for k1=1:Kb,
                        likelihood = res.intensity(n).lik(:,k1);
                        q(:,:,k1)  = q(:,:,k1).*likelihood(tmp);
                    end
                end
            end
            Q(:,:,z,:) = reshape(q,[d(1:2),1,Kb]); % final probabilities not normalized

        end
    end
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');

cls   = cell(1,Kb);
if do_cls
    P = zeros([d(1:3),Kb],'uint8');
    for z=1:length(x3),
        sq = sum(Q(:,:,z,:),4) + eps^2;
        for k1=1:Kb
            P(:,:,z,k1) = uint8(round(255 * Q(:,:,z,k1)./sq)); % nomalized probabilities
        end
    end
%     if mrf~=0, nmrf_its = 10; else nmrf_its = 1; end % andy 2013-05-03

%     if mrf==0, % Done this way as spm_mrf is not compiled for all platforms yet
    if mode == 0  % andy 2013-05-03
        sQ = (sum(Q,4)+eps)/255;
        for k1=1:size(Q,4)
            P(:,:,:,k1) = uint8(round(Q(:,:,:,k1)./sQ));
        end
        clear sQ
    else
        nmrf_its = 10; % andy 2013-05-03
        spm_progress_bar('init',nmrf_its,['MRF: Working on ' nam],'Iterations completed');
        vx2 = single(sum(res.image(1).mat(1:3,1:3).^2));
        if iscellstr(tcm) % if user specifies TCM, load it  % andy 2013-05-03
            load(tcm{1},'C');
            if ndims(C)~=2
                warning('Warn:convert',...
                    'Local or regional TCM detected. For MRF-based clean-up provided by SPM8, user-provided TCM must be a global TCM!\nSPM8 default config is used instead...\n');
                G   = ones([Kb,1],'single')*2;
            else
            G = single(C);  % andy 2013-04-22
            end
        else
            G   = ones([Kb,1],'single')*2; % andy 2013-05-03
            %save PQG P Q G tiss Kb x3 ind
        end

        for iter=1:nmrf_its,
            spm_mrf(P,Q,G,vx2);
            spm_progress_bar('set',iter);
        end
    end

    clear Q

    for k1=1:Kb,
        if ~isempty(tiss(k1).Nt),
            for z=1:length(x3),
                tmp = double(P(:,:,z,k1))/255;
                tiss(k1).Nt.dat(:,:,z,ind(1),ind(2)) = tmp;
            end
        end
    end
    spm_progress_bar('clear');

    for k1=1:Kb,
        if tc(k1,4) || any(tc(:,3)) || tc(k1,2) || nargout>=1,
            cls{k1} = P(:,:,:,k1);
        end
    end
    clear P
end

clear tpm
M0 = res.image(1).mat;

if any(tc(:,2)),

    bb = nan(2,3);
    vx = 1.5;
    % Sort out bounding box etc
    bb(~isfinite(bb)) = bb1(~isfinite(bb));
    if ~isfinite(vx), vx = abs(prod(vx1))^(1/3); end;
    bb(1,:) = vx*round(bb(1,:)/vx);
    bb(2,:) = vx*round(bb(2,:)/vx);

    % Figure out the mapping from the volumes to create to the original
    mm = [[
        bb(1,1) bb(1,2) bb(1,3)
        bb(2,1) bb(1,2) bb(1,3)
        bb(1,1) bb(2,2) bb(1,3)
        bb(2,1) bb(2,2) bb(1,3)
        bb(1,1) bb(1,2) bb(2,3)
        bb(2,1) bb(1,2) bb(2,3)
        bb(1,1) bb(2,2) bb(2,3)
        bb(2,1) bb(2,2) bb(2,3)]'; ones(1,8)];

    vx2  = M1\mm;
    odim = abs(round((bb(2,1:3)-bb(1,1:3))/vx))+1;
    vx3  = [[
        1       1       1
        odim(1) 1       1
        1       odim(2) 1
        odim(1) odim(2) 1
        1       1       odim(3)
        odim(1) 1       odim(3)
        1       odim(2) odim(3)
        odim(1) odim(2) odim(3)]'; ones(1,8)];

    x      = affind(rgrid(d),M0);
    y1     = affind(y,M1);
    ind    = find(tc(:,2));
    [M,R]  = spm_get_closest_affine(x,y1,single(cls{ind(1)})/255);
    clear x y1

    M      = M0\inv(R)*M1*vx2/vx3;
    mat0   =         R\M1*vx2/vx3;
    mat    = mm/vx3;

    fwhm = max(vx./sqrt(sum(res.image(1).mat(1:3,1:3).^2))-1,0.01);
    for k1=1:size(tc,1),
        if tc(k1,2),
            tmp1     = decimate(single(cls{k1}),fwhm);
            [pth,nam,ext1]=fileparts(res.image(1).fname);
            VT      = struct('fname',fullfile(pth,['rc', num2str(k1), nam, '.nii']),...
                'dim',  odim,...
                'dt',   [spm_type('float32') spm_platform('bigend')],...
                'pinfo',[1.0 0]',...
                'mat',mat);
            VT = spm_create_vol(VT);

            Ni             = nifti(VT.fname);
            Ni.mat0        = mat0;
            Ni.mat_intent  = 'Aligned';
            Ni.mat0_intent = 'Aligned';
            create(Ni);

            for i=1:odim(3),
                tmp = spm_slice_vol(tmp1,M*spm_matrix([0 0 i]),odim(1:2),[1,NaN])/255;
                VT  = spm_write_plane(VT,tmp,i);
            end
            clear tmp1
        end
    end
end

if any(tc(:,3)),
    C = zeros([d1,Kb],'single');
end

if any(tc(:,3)) || any(tc(:,4)) || nargout>=1,
    spm_progress_bar('init',Kb,'Warped Tissue Classes','Classes completed');
    for k1 = 1:Kb,
        if ~isempty(cls{k1}),
            c = single(cls{k1})/255;
            if any(tc(:,3)),
                [c,w]  = dartel3('push',c,y,d1(1:3));
                vx          = sqrt(sum(M1(1:3,1:3).^2));
                C(:,:,:,k1) = optimNn(w,c,[1  vx  1e-4 1e-6 0  3 2]);
                clear w
            else
                c      = dartel3('push',c,y,d1(1:3));
            end
            if nargout>=1,
                cls{k1} = c;
            end
            if tc(k1,4),
                N      = nifti;
                N.dat  = file_array(fullfile(pth,['mwc', num2str(k1), nam, '.nii']),...
                                    d1,...
                                    [spm_type('float32') spm_platform('bigend')],...
                                    0,1,0);
                N.mat  = M1;
                N.mat0 = M1;
                N.descrip = ['Jac. sc. warped tissue class ' num2str(k1)];
                create(N);
                N.dat(:,:,:) = c*abs(det(M0(1:3,1:3))/det(M1(1:3,1:3)));
            end
            spm_progress_bar('set',k1);
        end
    end
    spm_progress_bar('Clear');
end

if any(tc(:,3)),
    spm_progress_bar('init',Kb,'Writing Warped Tis Cls','Classes completed');
    C = max(C,eps);
    s = sum(C,4);
    for k1=1:Kb,
        if tc(k1,3),
            N      = nifti;
            N.dat  = file_array(fullfile(pth,['wc', num2str(k1), nam, '.nii']),...
                                d1,'uint8',0,1/255,0);
            N.mat  = M1;
            N.mat0 = M1;
            N.descrip = ['Warped tissue class ' num2str(k1)];
            create(N);
            N.dat(:,:,:) = C(:,:,:,k1)./s;
        end
        spm_progress_bar('set',k1);
    end
    spm_progress_bar('Clear');
    clear C s
end

if df(2),
    y         = spm_invert_def(y,M1,d1,M0,[1 0]);
    N         = nifti;
    N.dat     = file_array(fullfile(pth,['y_', nam1, '.nii']),...
                           [d1,1,3],'float32',0,1,0);
    N.mat     = M1;
    N.mat0    = M1;
    N.descrip = 'Deformation';
    create(N);
    N.dat(:,:,:,:,:) = reshape(y,[d1,1,3]);
end

return;
%=======================================================================

%=======================================================================
function [x1,y1,z1] = defs(sol,z,MT,prm,x0,y0,z0,M)
iMT = inv(MT);
x1  = x0*iMT(1,1)+iMT(1,4);
y1  = y0*iMT(2,2)+iMT(2,4);
z1  = (z0(z)*iMT(3,3)+iMT(3,4))*ones(size(x1));
x1a = x0    + spm_bsplins(sol{1},x1,y1,z1,prm);
y1a = y0    + spm_bsplins(sol{2},x1,y1,z1,prm);
z1a = z0(z) + spm_bsplins(sol{3},x1,y1,z1,prm);
x1  = M(1,1)*x1a + M(1,2)*y1a + M(1,3)*z1a + M(1,4);
y1  = M(2,1)*x1a + M(2,2)*y1a + M(2,3)*z1a + M(2,4);
z1  = M(3,1)*x1a + M(3,2)*y1a + M(3,3)*z1a + M(3,4);
return;
%=======================================================================

%=======================================================================
function t = transf(B1,B2,B3,T)
if ~isempty(T)
    d2 = [size(T) 1];
    t1 = reshape(reshape(T, d2(1)*d2(2),d2(3))*B3', d2(1), d2(2));
    t  = B1*t1*B2';
else
    t = zeros(size(B1,1),size(B2,1),size(B3,1));
end;
return;
%=======================================================================

%=======================================================================
function p = likelihoods(f,bf,mg,mn,vr)
K  = numel(mg);
N  = numel(f);
M  = numel(f{1});
cr = zeros(M,N);
for n=1:N,
    if isempty(bf),
        cr(:,n) = double(f{n}(:));
    else
        cr(:,n) = double(f{n}(:).*bf{n}(:));
    end
end
p  = ones(numel(f{1}),K);
for k=1:K,
    amp    = mg(k)/sqrt((2*pi)^N * det(vr(:,:,k)));
    d      = cr - repmat(mn(:,k)',M,1);
    p(:,k) = amp * exp(-0.5* sum(d.*(d/vr(:,:,k)),2));
end
%=======================================================================

%=======================================================================
function dat = decimate(dat,fwhm)
% Convolve the volume in memory (fwhm in voxels).
lim = ceil(2*fwhm);
x  = -lim(1):lim(1); x = spm_smoothkern(fwhm(1),x); x  = x/sum(x);
y  = -lim(2):lim(2); y = spm_smoothkern(fwhm(2),y); y  = y/sum(y);
z  = -lim(3):lim(3); z = spm_smoothkern(fwhm(3),z); z  = z/sum(z);
i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;
spm_conv_vol(dat,dat,x,y,z,-[i j k]);
return;
%=======================================================================

%=======================================================================
function y1 = affind(y0,M)
y1 = zeros(size(y0),'single');
for d=1:3,
    y1(:,:,:,d) = y0(:,:,:,1)*M(d,1) + y0(:,:,:,2)*M(d,2) + y0(:,:,:,3)*M(d,3) + M(d,4);
end
%=======================================================================

%=======================================================================
function x = rgrid(d)
x = zeros([d(1:3) 3],'single');
[x1,x2] = ndgrid(single(1:d(1)),single(1:d(2)));
for i=1:d(3),
    x(:,:,i,1) = x1;
    x(:,:,i,2) = x2;
    x(:,:,i,3) = single(i);
end
%=======================================================================


text:
function checkFactors(min_number, max_number)
%CHECKFACTORS   Return the maximum prime factor for a range of numbers.
%
% DESCRIPTION:
%       checkFactors loops through the given range of numbers and finds the
%       numbers with the smallest maximum prime factors. This allows
%       suitable grid sizes to be selected to maximise the speed of the FFT
%       (this is fastest for FFT lengths with small prime factors).
%
% USAGE:
%       checkFactors(min_number, max_number)
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 20th April 2011
%       last update - 8th July 2013
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

% extract factors
facs = zeros(1, max_number - min_number);
fac_max = facs;
for index = min_number:max_number;
    facs(index - min_number + 1) = length(factor(index));
    fac_max(index - min_number + 1) = max(factor(index));
end

% plot factors
figure;
subplot(2, 1, 1), bar(min_number:max_number, facs);
set(gca, 'XLim', [(min_number -0.5) (max_number + 0.5)]);
title('number of factors');
subplot(2, 1, 2), bar(min_number:max_number, fac_max);
set(gca, 'XLim', [(min_number -0.5) (max_number + 0.5)]);
title('largest factor');

% compute best factors in range
disp('Numbers with a maximum prime factor of 2');
ind = min_number + find(fac_max == 2) - 1;
disp(num2str(ind));
disp('Numbers with a maximum prime factor of 3');
ind = min_number + find(fac_max == 3) - 1;
disp(num2str(ind));
disp('Numbers with a maximum prime factor of 5');
ind = min_number + find(fac_max == 5) - 1;
disp(num2str(ind));
disp('Numbers with a maximum prime factor of 7');
ind = min_number + find(fac_max == 7) - 1;
disp(num2str(ind));
disp('Numbers to avoid (prime numbers)');
nums = min_number:max_number;
disp(num2str(nums(fac_max == nums)));

repo: Aero-matlab


text:
function varargout = setdiff(varargin)
% Adds support for 'stable' input flag with one output
if iscell(varargin{1}) && ~iscell(varargin{2}), varargin{2} = {varargin{2}}; end
if iscell(varargin{2}) && ~iscell(varargin{1}), varargin{1} = {varargin{1}}; end
if nargin>=3 && strcmp(varargin{end},'stable')
    if strcmp(varargin{3},'rows')
        y = [varargin{1}; varargin{2}];
        [a, ~, x] = unique(y,'rows');
        x1 = x(1:size(varargin{1},1));
        x2 = x(size(varargin{1},1)+1:end);
    else
        y = [varargin{1}(:); varargin{2}(:)];
        [a, ~, x] = unique(y);
        x1 = x(1:numel(varargin{1}));
        x2 = x(numel(varargin{1})+1:end);
    end
    ind = ~ismember(x1, x2);
    x = x1(ind);
    if ~isempty(x)
        x = x(~any(triu(bsxfun(@eq, x, x.'),1)));
    end
    if strcmp(varargin{3},'rows')
        y = a(x,:);
    else
        y = a(x);
    end
    if ~strcmp(varargin{3},'rows') && isrow(varargin{1}) % row output if first input is row
        y = reshape(y,1,[]);
    end
    varargout{1} = y;
else
    if strcmp(varargin{end},'sorted'), varargin(end) = []; end
    varargout = cell(1,max(nargout,1)); % if called without outputs, nargout will be zero. In that case we return one output
    [varargout{:}] = builtin('setdiff', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
% Runs the large tests for the whole cell simulation.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 8/13/2010
function runLargeTests()
%% initialize
warning('off','WholeCell:warning');
setWarnings();
setPath();
setPreferences();

%% run tests
import edu.stanford.covert.test.XMLTestRunDisplay;
import edu.stanford.covert.test.runtests;
monitor = XMLTestRunDisplay('Large Tests','Whole cell simulation large tests', 'output/runLargeTests/tests.xml');
runtests(monitor, {
    'edu.stanford.covert.cell.sim.Simulation_Large_Test'
    });

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [nch, npoints, freq, d] = ddt_v(filename)
% ddt_v(filename) Read data from a .ddt file returning samples in mV
%
% [nch, npoints, freq, d] = ddt_v(filename)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%
% OUTPUT:
%   nch - number of channels
%   npoints - number of data points for each channel
%   freq - A/D frequency
%   d - [nch npoints] data array (in mV)

if nargin < 1
    error 'Expected 1 input argument';
end
if (isempty(filename))
   [fname, pathname] = uigetfile('*.ddt', 'Select a Plexon .ddt file');
   if isequal(fname,0)
     error 'No file was selected'
   end
   filename = fullfile(pathname, fname);
end

[nch, npoints, freq, d] = mexPlex(20, filename);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function bits_per_pixel = fmf_get_bits_per_pixel( data_format )

if strcmp(data_format,'MONO8'),
  bits_per_pixel = 8;
elseif strcmp(data_format,'YUV422'),
  bits_per_pixel = 16;
else,
  error( 'Only MONO8 supported for now' );
end

author: kristinbranson

repo: JAABA


text:
function Q = diffusionRWR(A, maxiter, restartProb)
	n = size(A, 1);

	% Add self-edge to isolated nodes
	A = A + diag(sum(A) == 0);

	% Normalize the adjacency matrix
	renorm = @(M) bsxfun(@rdivide, M, sum(M));
	P = renorm(A);

	% Personalized PageRank
	restart = eye(n);
	Q = eye(n);
	for i = 1 : maxiter
		Q_new = (1 - restartProb) * P * Q + restartProb * restart;
		delta = norm(Q - Q_new, 'fro');
		%fprintf('Iter %d. Frobenius norm: %f\n', i, delta);
		Q = Q_new;
		if delta < 1e-6
			%fprintf('Converged.\n');
			break;
        end
    end
end

author: luoyunan

repo: DTINet

save_path: github-repos/MATLAB/luoyunan-DTINet


text:
function x = c8vec_sort_a2 ( n, x )

%*****************************************************************************80
%
%% C8VEC_SORT_A2 ascending sorts a complex array by L2 norm.
%
%  Discussion:
%
%    The L2 norm of A+Bi is sqrt ( A**2 + B**2 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, length of input array.
%
%    Input, complex X(N), an unsorted array.
%
%    Output, complex X(N), a sorted copy of the array.
%
  indx = 0;
  isgn = 0;

  while ( 1 )

    [ indx, i, j ] = sort_heap_external ( n, indx, isgn );

    if ( 0 < indx )

      t = x(i);
      x(i) = x(j);
      x(j) = t;

    elseif ( indx < 0 )

      if ( c8_le_l2 ( x(i), x(j) ) )
        isgn = -1;
      else
        isgn = +1;
      end

    elseif ( indx == 0 )

      break;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m