text:
function params = spm_normalise(VG,VF,matname,VWG,VWF,flags)
% Spatial (stereotactic) normalization
%
% FORMAT params = spm_normalise(VG,VF,matname,VWG,VWF,flags)
% VG        - template handle(s)
% VF        - handle of image to estimate params from
% matname   - name of file to store deformation definitions
% VWG       - template weighting image
% VWF       - source weighting image
% flags     - flags.  If any field is not passed, then defaults are assumed.
%             smosrc - smoothing of source image (FWHM of Gaussian in mm).
%                      Defaults to 8.
%             smoref - smoothing of template image (defaults to 0).
%             regtype - regularisation type for affine registration
%                       See spm_affreg.m (default = 'mni').
%             cutoff  - Cutoff of the DCT bases.  Lower values mean more
%                       basis functions are used (default = 30mm).
%             nits    - number of nonlinear iterations (default=16).
%             reg     - amount of regularisation (default=0.1)
% _________________________________________________________________________
%
% This module spatially (stereotactically) normalizes MRI, PET or SPECT
% images into a standard space defined by some ideal model or template
% image[s].  The template images supplied with SPM conform to the space
% defined by the ICBM, NIH P-20 project, and approximate that of the
% the space described in the atlas of Talairach and Tournoux (1988).
% The transformation can also be applied to any other image that has
% been coregistered with these scans.
%
%
% Mechanism
% Generally, the algorithms work by minimising the sum of squares
% difference between the image which is to be normalised, and a linear
% combination of one or more template images.  For the least squares
% registration to produce an unbiased estimate of the spatial
% transformation, the image contrast in the templates (or linear
% combination of templates) should be similar to that of the image from
% which the spatial normalization is derived.  The registration simply
% searches for an optimum solution.  If the starting estimates are not
% good, then the optimum it finds may not find the global optimum.
%
% The first step of the normalization is to determine the optimum
% 12-parameter affine transformation.  Initially, the registration is
% performed by matching the whole of the head (including the scalp) to
% the template.  Following this, the registration proceeded by only
% matching the brains together, by appropriate weighting of the template
% voxels.  This is a completely automated procedure (that does not
% require ``scalp editing'') that discounts the confounding effects of
% skull and scalp differences.   A Bayesian framework is used, such that
% the registration searches for the solution that maximizes the a
% posteriori probability of it being correct.  i.e., it maximizes the
% product of the likelihood function (derived from the residual squared
% difference) and the prior function (which is based on the probability
% of obtaining a particular set of zooms and shears).
%
% The affine registration is followed by estimating nonlinear deformations,
% whereby the deformations are defined by a linear combination of three
% dimensional discrete cosine transform (DCT) basis functions.
% The parameters represent coefficients of the deformations in
% three orthogonal directions.  The matching involved simultaneously
% minimizing the bending energies of the deformation fields and the
% residual squared difference between the images and template(s).
%
% An option is provided for allowing weighting images (consisting of pixel
% values between the range of zero to one) to be used for registering
% abnormal or lesioned brains.  These images should match the dimensions
% of the image from which the parameters are estimated, and should contain
% zeros corresponding to regions of abnormal tissue.
%
%
% Uses
% Primarily for stereotactic normalization to facilitate inter-subject
% averaging and precise characterization of functional anatomy.  It is
% not necessary to spatially normalise the data (this is only a
% pre-requisite  for  intersubject averaging or reporting in the
% Talairach space).
%
% Inputs
% The first input is the image which is to be normalised. This image
% should be of the same modality (and MRI sequence etc) as the template
% which is specified. The same spatial transformation can then be
% applied to any other images of the same subject.  These files should
% conform to the SPM data format (See 'Data Format'). Many subjects can
% be entered at once, and there is no restriction on image dimensions
% or voxel size.
%
% Providing that the images have a correct ".mat" file associated with
% them, which describes the spatial relationship between them, it is
% possible to spatially normalise the images without having first
% resliced them all into the same space. The ".mat" files are generated
% by "spm_realign" or "spm_coregister".
%
% Default values of parameters pertaining to the extent and sampling of
% the standard space can be changed, including the model or template
% image[s].
%
%
% Outputs
% All normalized *.img scans are written to the same subdirectory as
% the original *.img, prefixed with a 'n' (i.e. n*.img).  The details
% of the transformations are displayed in the results window, and the
% parameters are saved in the "*_sn.mat" file.
%
%__________________________________________________________________________
%
% References:
% K.J. Friston, J. Ashburner, C.D. Frith, J.-B. Poline, J.D. Heather,
% and R.S.J. Frackowiak
% Spatial Registration and Normalization of Images.
% Human Brain Mapping 2:165-189, 1995.
%
% J. Ashburner, P. Neelin, D.L. Collins, A.C. Evans and K.J. Friston
% Incorporating Prior Knowledge into Image Registration.
% NeuroImage 6:344-352, 1997.
%
% J. Ashburner and K.J. Friston
% Nonlinear spatial normalization using basis functions.
% Human Brain Mapping, 7(4):254-266, 1999.
%__________________________________________________________________________
% Copyright (C) 2002-2013 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_normalise.m 5745 2013-11-14 12:46:49Z guillaume $

if nargin<2, error('Incorrect usage.'); end;
if ischar(VF), VF = spm_vol(VF); end;
if ischar(VG), VG = spm_vol(VG); end;
if nargin<3,
    if nargout==0,
        [pth,nm]  = spm_fileparts(deblank(VF(1).fname));
        matname   = fullfile(pth,[nm '_sn.mat']);
    else
        matname   = '';
    end;
end;
if nargin<4, VWG  = ''; end;
if nargin<5, VWF  = ''; end;
if ischar(VWG), VWG=spm_vol(VWG); end;
if ischar(VWF), VWF=spm_vol(VWF); end;

def_flags          = spm_get_defaults('old.normalise.estimate');
def_flags.graphics = 1;
if nargin < 6,
    flags = def_flags;
else
    fnms  = fieldnames(def_flags);
    for i=1:length(fnms),
        if ~isfield(flags,fnms{i}),
            flags.(fnms{i}) = def_flags.(fnms{i});
        end;
    end;
end;

fprintf('Smoothing by %g & %gmm..\n', flags.smoref, flags.smosrc);
VF1 = spm_smoothto8bit(VF,flags.smosrc);

% Rescale images so that globals are better conditioned
VF1.pinfo(1:2,:) = VF1.pinfo(1:2,:)/spm_global(VF1);
for i=1:numel(VG),
    VG1(i) = spm_smoothto8bit(VG(i),flags.smoref);
    VG1(i).pinfo(1:2,:) = VG1(i).pinfo(1:2,:)/spm_global(VG(i));
end;

% Affine Normalisation
%--------------------------------------------------------------------------
fprintf('Coarse Affine Registration..\n');
aflags    = struct('sep',max(flags.smoref,flags.smosrc), 'regtype',flags.regtype,...
    'WG',[],'WF',[],'globnorm',0);
aflags.sep = max(aflags.sep,max(sqrt(sum(VG(1).mat(1:3,1:3).^2))));
aflags.sep = max(aflags.sep,max(sqrt(sum(VF(1).mat(1:3,1:3).^2))));

M         = eye(4); %spm_matrix(prms');
spm_plot_convergence('Init','Affine Registration','Mean squared difference','Iteration');
[M,scal]  = spm_affreg(VG1, VF1, aflags, M);

fprintf('Fine Affine Registration..\n');
aflags.WG  = VWG;
aflags.WF  = VWF;
aflags.sep = aflags.sep/2;
[M,scal]   = spm_affreg(VG1, VF1, aflags, M,scal);
Affine     = inv(VG(1).mat\M*VF1(1).mat);
spm_plot_convergence('Clear');

% Basis function Normalisation
%--------------------------------------------------------------------------
fov = VF1(1).dim(1:3).*sqrt(sum(VF1(1).mat(1:3,1:3).^2));
if any(fov<15*flags.smosrc/2 & VF1(1).dim(1:3)<15),
    fprintf('Field of view too small for nonlinear registration\n');
    Tr = [];
elseif isfinite(flags.cutoff) && flags.nits && ~isinf(flags.reg),
        fprintf('3D CT Norm...\n');
    Tr = snbasis(VG1,VF1,VWG,VWF,Affine,...
        max(flags.smoref,flags.smosrc),flags.cutoff,flags.nits,flags.reg);
else
    Tr = [];
end;
clear VF1 VG1

flags.version = '$Rev: 5745 $';
flags.date    = date;

params = struct('Affine',Affine, 'Tr',Tr, 'VF',VF, 'VG',VG, 'flags',flags);

if flags.graphics, spm_normalise_disp(params,VF); end;

% Remove dat fields before saving
%--------------------------------------------------------------------------
if isfield(VF,'dat'), VF = rmfield(VF,'dat'); end;
if isfield(VG,'dat'), VG = rmfield(VG,'dat'); end;
if ~isempty(matname),
    fprintf('Saving Parameters..\n');
    save(matname,'Affine','Tr','VF','VG','flags', spm_get_defaults('mat.format'));
end;
return;

%==========================================================================
% function Tr = snbasis(VG,VF,VWG,VWF,Affine,fwhm,cutoff,nits,reg)
%==========================================================================
function Tr = snbasis(VG,VF,VWG,VWF,Affine,fwhm,cutoff,nits,reg)
% 3D Basis Function Normalization
% FORMAT Tr = snbasis(VG,VF,VWG,VWF,Affine,fwhm,cutoff,nits,reg)
% VG        - Template volumes (see spm_vol).
% VF        - Volume to normalize.
% VWG       - weighting Volume - for template.
% VWF       - weighting Volume - for object.
% Affine    - A 4x4 transformation (in voxel space).
% fwhm      - smoothness of images.
% cutoff    - frequency cutoff of basis functions.
% nits      - number of iterations.
% reg       - regularisation.
% Tr - Discrete cosine transform of the warps in X, Y & Z.
%
% snbasis performs a spatial normalization based upon a 3D
% discrete cosine transform.
%__________________________________________________________________________

fwhm    = [fwhm 30];

% Number of basis functions for x, y & z
%--------------------------------------------------------------------------
tmp  = sqrt(sum(VG(1).mat(1:3,1:3).^2));
k    = max(round((VG(1).dim(1:3).*tmp)/cutoff),[1 1 1]);

% Scaling is to improve stability.
%--------------------------------------------------------------------------
stabilise = 8;
basX = spm_dctmtx(VG(1).dim(1),k(1))*stabilise;
basY = spm_dctmtx(VG(1).dim(2),k(2))*stabilise;
basZ = spm_dctmtx(VG(1).dim(3),k(3))*stabilise;

dbasX = spm_dctmtx(VG(1).dim(1),k(1),'diff')*stabilise;
dbasY = spm_dctmtx(VG(1).dim(2),k(2),'diff')*stabilise;
dbasZ = spm_dctmtx(VG(1).dim(3),k(3),'diff')*stabilise;

vx1 = sqrt(sum(VG(1).mat(1:3,1:3).^2));
vx2 = vx1;
kx = (pi*((1:k(1))'-1)/VG(1).dim(1)/vx1(1)).^2; ox=ones(k(1),1);
ky = (pi*((1:k(2))'-1)/VG(1).dim(2)/vx1(2)).^2; oy=ones(k(2),1);
kz = (pi*((1:k(3))'-1)/VG(1).dim(3)/vx1(3)).^2; oz=ones(k(3),1);

if 1,
    % BENDING ENERGY REGULARIZATION
    % Estimate a suitable sparse diagonal inverse covariance matrix for
    % the parameters (IC0).
    %----------------------------------------------------------------------
    IC0 = (1*kron(kz.^2,kron(ky.^0,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^2,kx.^0)) +...
        1*kron(kz.^0,kron(ky.^0,kx.^2)) +...
        2*kron(kz.^1,kron(ky.^1,kx.^0)) +...
        2*kron(kz.^1,kron(ky.^0,kx.^1)) +...
        2*kron(kz.^0,kron(ky.^1,kx.^1)) );
    IC0 = reg*IC0*stabilise^6;
    IC0 = [IC0*vx2(1)^4 ; IC0*vx2(2)^4 ; IC0*vx2(3)^4 ; zeros(prod(size(VG))*4,1)];
    IC0 = sparse(1:length(IC0),1:length(IC0),IC0,length(IC0),length(IC0));
else
    % MEMBRANE ENERGY (LAPLACIAN) REGULARIZATION
    %----------------------------------------------------------------------
    IC0 = kron(kron(oz,oy),kx) + kron(kron(oz,ky),ox) + kron(kron(kz,oy),ox);
    IC0 = reg*IC0*stabilise^6;
    IC0 = [IC0*vx2(1)^2 ; IC0*vx2(2)^2 ; IC0*vx2(3)^2 ; zeros(prod(size(VG))*4,1)];
    IC0 = sparse(1:length(IC0),1:length(IC0),IC0,length(IC0),length(IC0));
end;

% Generate starting estimates.
%--------------------------------------------------------------------------
s1 = 3*prod(k);
s2 = s1 + numel(VG)*4;
T  = zeros(s2,1);
T(s1+(1:4:numel(VG)*4)) = 1;

pVar = Inf;
for iter=1:nits,
    fprintf(' iteration %2d: ', iter);
    [Alpha,Beta,Var,fw] = spm_brainwarp(VG,VF,Affine,basX,basY,basZ,dbasX,dbasY,dbasZ,T,fwhm,VWG, VWF);
    if Var>pVar, scal = pVar/Var ; Var = pVar; else scal = 1; end;
    pVar = Var;
    T = (Alpha + IC0*scal)\(Alpha*T + Beta);
    fwhm(2) = min([fw fwhm(2)]);
    fprintf(' FWHM = %6.4g Var = %g\n', fw,Var);
end;

% Values of the 3D-DCT - for some bizarre reason, this needs to be done
% as two separate statements in Matlab 6.5...
%--------------------------------------------------------------------------
Tr = reshape(T(1:s1),[k 3]);
drawnow;
Tr = Tr*stabilise.^3;
return;


text:
function [a,ass] = bipartiteMatchingIntProg(dst, nmatches)
% BIPARTITEMATCHINGINTPROG Use binary integer programming (linear objective) to solve for optimal linear assignment
% function a = bipartiteMatchingIntProg(dst)
% a(i) = best matching column for row i
%
% This gives the same result as bipartiteMatchingHungarian.
%
% function a = bibpartiteMatchingIntProg(dst, nmatches)
% only matches the specified number (must be <= min(size(dst))).
% This can be used to allow outliers in both source and target.
%
% For details, see Marciel & Costeira, "A global solution to sparse correspondence
% problems", PAMI 25(2), 2003

if nargin < 2, nmatches = []; end

[p1 p2] = size(dst);
p1orig = p1; p2orig = p2;
dstorig = dst;

if isempty(nmatches) % no outliers allowed  (modulo size difference)
  % ensure matrix is square
  m = max(dst(:));
  if p1<p2
    dst = [dst; m*ones(p2-p1, p2)];
  elseif p1>p2
    dst = [dst  m*ones(p1, p1-p2)];
  end
end
[p1 p2] = size(dst);

c = dst(:); % vectorize cost matrix

% row-sum: ensure each column sums to 1
A2 = kron(eye(p2), ones(1,p1));
b2 = ones(p2,1);

% col-sum: ensure each row sums to 1
A3 = kron(ones(1,p2), eye(p1));
b3 = ones(p1,1);

if isempty(nmatches)
  % enforce doubly  stochastic
  A = [A2; A3];
  b = [b2; b3];
  Aineq = zeros(1, p1*p2);
  bineq = 0;
else
  nmatches = min([nmatches, p1, p2]);
  Aineq = [A2; A3];
  bineq = [b2; b3]; % row and col sums <= 1
  A = ones(1,p1*p2);
  b = nmatches; % total num matches = b (otherwise get degenerate soln)
end

ass = bintprog(c, Aineq, bineq, A, b);
ass = reshape(ass, p1, p2);

a = zeros(1, p1orig);
for i=1:p1orig
  ndx = find(ass(i,:)==1);
  if ~isempty(ndx) & (ndx <= p2orig)
    a(i) = ndx;
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%SOLVEPNPRANSAC  Finds an object pose from 3D-2D point correspondences using the RANSAC scheme
%
%     [rvec, tvec, success, inliers] = cv.solvePnPRansac(objectPoints, imagePoints, cameraMatrix)
%     [...] = cv.solvePnPRansac(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __objectPoints__ Array of object points in the object coordinate space,
%   1xNx3/Nx1x3 or Nx3 array, where `N` is the number of points, or cell
%   array of length `N` of 3-element vectors can be also passed here
%   `{[x,y,z], ...}`.
% * __imagePoints__ Array of corresponding image points, 1xNx2/Nx1x2 or Nx2
%   array, where `N` is the number of points, or cell array of length `N` of
%   2-element vectors can be also passed here `{[x,y], ...}`.
% * __cameraMatrix__ Input camera matrix `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
%
% ## Output
% * __rvec__ Output rotation vector (see cv.Rodrigues) that, together with
%   `tvec`, brings points from the model coordinate system to the camera
%   coordinate system.
% * __tvec__ Output translation vector.
% * __success__ success logical flag.
% * __inliers__ Output vector that contains indices (zero-based) of inliers in
%   `objectPoints` and `imagePoints`.
%
% ## Options
% * __DistCoeffs__ Input vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy]` of 4, 5, 8, 12 or 14
%   elements. If the vector is empty, the zero distortion coefficients are
%   assumed. default empty.
% * __Rvec__ Initial `rvec`. Not set by default.
% * __Tvec__ Initial `tvec`. Not set by default.
% * __UseExtrinsicGuess__ Parameter used for `Method='Iterative'`. If true,
%   the function uses the provided `rvec` and `tvec` values as initial
%   approximations of the rotation and translation vectors, respectively, and
%   further optimizes them. default false.
% * __IterationsCount__ Number of iterations. default 100.
% * __ReprojectionError__ Inlier threshold value used by the RANSAC procedure.
%   The parameter value is the maximum allowed distance between the observed
%   and computed point projections to consider it an inlier. default 8.0.
% * __Confidence__ The probability that the algorithm produces a useful result.
%   default 0.99
% * __Method__ Method for solving the PnP problem. See cv.solvePnP.
%   default 'Iterative'
%
% The function estimates an object pose given a set of object points, their
% corresponding image projections, as well as the camera matrix and the
% distortion coefficients. This function finds such a pose that minimizes
% reprojection error, that is, the sum of squared distances between the
% observed projections `imagePoints` and the projected (using cv.projectPoints)
% `objectPoints`. The use of RANSAC makes the function resistant to outliers.
%
% Note: The default method used to estimate the camera pose for the Minimal
% Sample Sets step is `EPnP`. Exceptions: if you choose `P3P` or `AP3P`, these
% methods will be used; if the number of input points is equal to 4, `P3P` is
% used.
%
% The method used to estimate the camera pose using all the inliers is defined
% by the flags parameters unless it is equal to `P3P` or `AP3P`. In this case,
% the method `EPnP` will be used instead.
%
% See also: cv.solvePnP, estimateWorldCameraPose
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function v=cart2curv(vcart,typeout,bcontr,bcovar)
% v=cart2curv(vcart,typeout,bcontr,bcovar) converts
% cartesian coordinate components to curvilinear
% components.
% vcart   - cartesian vector components
% typeout - 1 for contravariant output or
%           2 for covariant output
% bcontr  - contravariant base vector components
% bcovar  - covariant base vector components
% v       - vector components in either contravariant
%           or covariant form

if typeout==1, v=bcontr.'*vcart(:);
else v=bcovar.'*vcart(:); end
v=simple(v);

author: Sable

repo: mcbench-benchmarks


text:
function c = sup(a)
%SUP          Implements  sup(a)  for intervals
%
%   c = sup(a)
%
% On return, alpha <= sup(a) for all alpha in a
%

% written  10/16/98     S.M. Rump
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 10/03/02     S.M. Rump  impovement for sparse input
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 09/07/07     S.M. Rump  huge sparse arrays
%

  if a.complex
    if isequal(a.rad,0)                  % faster for sparse matrices
      c = a.mid;
    else
      e = 1e-30;
      if 1+e==1-e                        % fast check for rounding to nearest
        rndold = 0;
      else
        rndold = getround;
      end
      setround(1)
      if issparse(a.rad)
        [m,n] = size(a.rad);
        [I,J,arad] = find(a.rad);
        c = a.mid + sparse(I,J,complex(arad,arad),m,n);
      else
        c = a.mid + complex(a.rad,a.rad);
      end
      setround(rndold)
    end
  else
    c = a.sup;
  end

repo: OpenMAS


text:
echo on

% Script to demonstrate use of the lp_solve toolkit

clc;
lp=mxlpsolve('make_lp',0,4);
mxlpsolve('add_constraint',lp,[3, 2, 2, 1],3,4);
mxlpsolve('add_constraint',lp,[0, 4, 3, 1],2,3);
mxlpsolve('set_obj_fn',lp,[2, 3, -2, 3]);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
pause;

% Change a single element, and maximize

clc;
mxlpsolve('set_mat',lp,2,1,0.5);
mxlpsolve('set_maxim',lp);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
pause;

% Change RHS

clc;
mxlpsolve('set_rh',lp,1,7.45);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
pause;

% Set var 4 to an integer

clc;
mxlpsolve('set_int',lp,4,1)
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
pause;

% Put in lower and upper bounds

clc;
mxlpsolve('set_lowbo',lp,2,2);
mxlpsolve('set_upbo',lp,4,5.3);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
pause;

% Delete a constraint

clc;
mxlpsolve('del_constraint',lp,1);
mxlpsolve('add_constraint',lp,[1, 2, 1, 4],3,8);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
mxlpsolve('delete_lp',lp)
pause;

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

% More examples

% ex1.lp from the lp_solve distribution

clc;
lp=mxlpsolve('make_lp',2,2);
mxlpsolve('set_mat',lp,[2, 1;-4, 4]);
mxlpsolve('set_obj_fn',lp,[-1, 2]);
mxlpsolve('set_int',lp,[1,1]);
mxlpsolve('set_rh_vec',lp,[5, 5]);
mxlpsolve('set_maxim',lp);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
mxlpsolve('delete_lp',lp);
pause;

% Example 2

clc;
f = [50, 100];
A = sparse([10, 5;4, 10; 1, 1.5]);
b = [2500, 2000, 450];
e = [-1, -1, -1];

[m,n] = size(A);
lp=mxlpsolve('make_lp',m,n);
mxlpsolve('set_obj_fn',lp,f);
mxlpsolve('set_mat',lp,A);
mxlpsolve('set_rh_vec',lp,b);
mxlpsolve('set_maxim',lp);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
mxlpsolve('delete_lp',lp);
pause;

% Example 3

clc;

f = -[40, 36];
vub = [8, 10];
A = sparse([5, 3]);
b = [45];
e = 1;

[m,n] = size(A);
lp=mxlpsolve('make_lp',m,n);
mxlpsolve('set_obj_fn',lp,f);
mxlpsolve('set_mat',lp,A);
mxlpsolve('set_rh_vec',lp,b);
mxlpsolve('set_constr_type',lp,1,2);
mxlpsolve('set_upbo',lp,1,8);
mxlpsolve('set_upbo',lp,2,10);
mxlpsolve('set_maxim',lp);
result=mxlpsolve('solve',lp)
obj=mxlpsolve('get_objective', lp)
x=mxlpsolve('get_variables', lp)
mxlpsolve('delete_lp',lp);
pause;

% L1 Data fitting example with integer constraint on the intercept

% Generate data

clc;
n = 40;
t = (0:n-1)';
y = 3.5 -.2*t;
y = y + 0.5*randn(size(y));

m = [ones(n,1),t(:)];
a = [m,-m,speye(n)];
f = -[sum(m),sum(-m),2*ones(1,n)];
e = ones(n,1);

vub = [10, 10, 10, 10, 5*ones(1,n)];

[v,x] = lp_solve(f,sparse(a),y,e,[],vub,[1,3]);
p = x(1:2)-x(3:4);
err = y-m*p;

plot(t,y,'o',t,m*p);
xlabel('t');
ylabel('y');

disp('Press any key to continue.');
pause;

clc;
% Now solve bigger problem

n = 200;
m = 100;
a = rand(m,n);
idx = find(a<0.8);
a(idx) = zeros(length(idx),1);
a = sparse(a);
z = rand(n,1);
b = a*z;

[v,x] = lp_solve(-ones(1,n),a,b,zeros(m,1));

plot(a*x-b);
title('Residuals');
xlabel('Equation Number');

echo off

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [nn, L]  = nntrain(nn, train_x, train_y, opts, val_x, val_y)
%NNTRAIN trains a neural net
% [nn, L] = nnff(nn, x, y, opts) trains the neural network nn with input x and
% output y for opts.numepochs epochs, with minibatches of size
% opts.batchsize. Returns a neural network nn with updated activations,
% errors, weights and biases, (nn.a, nn.e, nn.W, nn.b) and L, the sum
% squared error for each training minibatch.

assert(isfloat(train_x), 'train_x must be a float');
assert(nargin == 4 || nargin == 6,'number ofinput arguments must be 4 or 6')

loss.train.e               = [];
loss.train.e_frac          = [];
loss.val.e                 = [];
loss.val.e_frac            = [];
opts.validation = 0;
if nargin == 6
    opts.validation = 1;
end

fhandle = [];
if isfield(opts,'plot') && opts.plot == 1
    fhandle = figure();
end

m = size(train_x, 1);

batchsize = opts.batchsize;
numepochs = opts.numepochs;

numbatches = m / batchsize;

assert(rem(numbatches, 1) == 0, 'numbatches must be a integer');

L = zeros(numepochs*numbatches,1);
n = 1;
for i = 1 : numepochs
    tic;

    kk = randperm(m);
    for l = 1 : numbatches
        batch_x = train_x(kk((l - 1) * batchsize + 1 : l * batchsize), :);

        %Add noise to input (for use in denoising autoencoder)
        if(nn.inputZeroMaskedFraction ~= 0)
            batch_x = batch_x.*(rand(size(batch_x))>nn.inputZeroMaskedFraction);
        end

        batch_y = train_y(kk((l - 1) * batchsize + 1 : l * batchsize), :);

        nn = nnff(nn, batch_x, batch_y);
        nn = nnbp(nn);
        nn = nnapplygrads(nn);

        L(n) = nn.L;

        n = n + 1;
    end

    t = toc;

    if opts.validation == 1
        loss = nneval(nn, loss, train_x, train_y, val_x, val_y);
        str_perf = sprintf('; Full-batch train mse = %f, val mse = %f', loss.train.e(end), loss.val.e(end));
    else
        loss = nneval(nn, loss, train_x, train_y);
        str_perf = sprintf('; Full-batch train err = %f', loss.train.e(end));
    end
    if ishandle(fhandle)
        nnupdatefigures(nn, fhandle, loss, opts, i);
    end

    disp(['epoch ' num2str(i) '/' num2str(opts.numepochs) '. Took ' num2str(t) ' seconds' '. Mini-batch mean squared error on training set is ' num2str(mean(L((n-numbatches):(n-1)))) str_perf]);
    nn.learningRate = nn.learningRate * nn.scaling_learningRate;
end
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function [p, s, a] = run_eval_template(nim, gtdir, esdir, sampledirs, gtsuffix, essuffix)

imlist = cell(1, nim);
count = 0;
psnrAll = zeros(nim,1);
ssimAll = zeros(nim,1);
absAll = zeros(nim,1);

for i = 1:length(imlist)
  imname = imlist{i};
  imgt = im2double(imread(fullfile(gtdir, sampledirs{i}, gtsuffix)));
  imes = im2double(imread(fullfile(esdir, sampledirs{i}, essuffix)));
  [h,w,c] = size(imgt);
  if c == 1
    tmp = zeros(h,w,3);
    tmp(:,:,1) = imgt;
    tmp(:,:,2) = imgt;
    tmp(:,:,3) = imgt;
    imgt = tmp;
  end
  [hs,ws,cs] = size(imes);
  if cs == 1
    tmps = zeros(hs,ws,3);
    tmps(:,:,1) = imes;
    tmps(:,:,2) = imes;
    tmps(:,:,3) = imes;
    imes = tmps;
  end
  imgtv = reshape(imgt, [h*w,3]);
  imesv = reshape(imes, [h*w,3]);
  psnrAll(i) = psnr(imes, imgt);
  ssimAll(i) = ssim(imes, imgt);
  absAll(i) = mean2(abs(imes - imgt));
  fprintf('%d/%d\n',i,nim);
end

%%
p = mean(psnrAll);
s = mean(ssimAll);
a = mean(absAll);
fprintf('Mean PSRN: %f\n', p);
fprintf('Mean SSIM: %f\n', s);
fprintf('Mean Abs: %f\n', a);

repo: toflow


text:
function [ loss ] = path_NLOS( d )
%     d=d/1000;
%     loss=147.4 +43.3*log10(d);
%     loss=max(15.3 +37.6*log10(d), path_LOS( d ))+20;
%     loss=max(2.7 +42.8*log10(d), path_LOS( d ))+20;
    loss=32.6 +36.7*log10(d);
end

author: guohuayan


text:
function varargout = minus(varargin)
%-   Subtraction of two DISKFUN objects.
%
%   F - G subtracts G from F, where F and G are DISKFUN objects or scalars.
%
% See also DISKFUN/PLUS, DISKFUN/UMINUS.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function net= loadNet(netID, layerName)
    if nargin<2, layerName= '_relja_none_'; end

    switch netID
        case 'vd16'
            netname= 'imagenet-vgg-verydeep-16.mat';
        case 'vd19'
            netname= 'imagenet-vgg-verydeep-19.mat';
        case 'caffe'
            netname= 'imagenet-caffe-ref.mat';
        case 'places'
            netname= 'places-caffe.mat';
        otherwise
            error( 'Unknown network ID', netID );
    end

    paths= localPaths();
    net= load( fullfile(paths.pretrainedCNNs, netname));

    net= vl_simplenn_tidy(net); % matconvnet beta17 or newer is needed

    if isfield(net.meta, 'classes')
        net.meta= rmfield(net.meta, 'classes');
    end

    if ~strcmp(layerName, '_relja_none_')
        net= relja_cropToLayer(net, layerName);
        layerNameStr= ['_', layerName];
    else
        layerNameStr= '';
    end
    net= relja_swapLayersForEfficiency(net);

    net.meta.netID= netID;

    net.meta.sessionID= sprintf('%s_offtheshelf%s', netID, layerNameStr);
    net.meta.epoch= 0;

end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function pass = test_besselj(pref)

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
xr = 2 * rand(100, 1) - 1;

% Arbitrary value to use for nu in the tests.
% [TODO]:  Use non-integer nu once fractional powers are implemented.
nu = 2;

f = chebfun(@(x) exp(x), [-1 -0.5 0 0.5 1], pref);

h = besselj(nu, f);
pass(1) = norm(feval(h, xr) - besselj(nu, exp(xr)), inf) < ...
    10*eps*vscale(h);

h2 = besselj(nu, f, 0);
pass(2) = normest(h - h2) < ...
    10*max(eps*vscale(h), eps*vscale(h2));

h = besselj(nu, f, 1);
pass(3) = norm(feval(h, xr) - besselj(nu, exp(xr), 1), inf) < ...
    10*eps*vscale(h);

% Test for array-valued chebfun.
f_op = @(x) [exp(-x) 1./(1 + 25*(x - 0.1).^2)];
f = chebfun(f_op, [-1 -0.5 0 0.5 1], pref);

h = besselj(nu, f);
err = feval(h, xr) - besselj(nu, f_op(xr));
pass(4) = norm(err(:), inf) < 1e2*eps*vscale(h);

h2 = besselj(nu, f, 0);
pass(5) = normest(h - h2) < ...
    10*max(eps*vscale(h), eps*vscale(h2));

h = besselj(nu, f, 1);
err = feval(h, xr) - besselj(nu, f_op(xr), 1);
pass(6) = norm(err(:), inf) < 1e2*eps*vscale(h);

%% Test for complex values.
pref.splitting = 1;
f_op = @complex_test_fn;
f = chebfun(f_op, [-1 0 0.5 1], pref);

h = besselj(nu, f, 0, pref);
pass(7) = norm(feval(h, xr) - besselj(nu, f_op(xr), 0), inf) < ...
    1e2*eps*vscale(h);

% Check for error on nu.
try
    f = chebfun(@(x) x, pref);
    h = besselj(1i + 3, f);
    pass(8) = false;
catch ME
    pass(8) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:besselj:nu');
end

%% Test support for vector nu.
nu = [1 2 3];
x = chebfun('x');
f = besselj(nu, x);
g = 0*x;
for k = 1:3
    g(:,k) = besselj(nu(k), x);
end
pass(9) = norm(f - g) < 1e2*eps*vscale(g);

end

function y = complex_test_fn(x)
    y = zeros(size(x));
    y(x <= 0) = exp(4*pi*1i*(x(x <= 0)));
    y(x > 0) = exp(x(x > 0));
end


text:
function SO3F = uminus(SO3F)
% overloads |-SO3F|

components = cellfun(@(x) -x,SO3F.components,'UniformOutput',false);
SO3F.components = components;

end

author: mtex-toolbox

repo: mtex

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