text:
function sequence_streak_display ( file_name, direction )

%*****************************************************************************80
%
%% SEQUENCE_STREAK_DISPLAY makes a "streak plot" of a sequence.
%
%  Discussion:
%
%    The sequence information is assumed to be available in
%    a text file, with one sequence entry per line.
%
%    One reason for making plots this way is to try to
%    see the "fill in" pattern of a sequence.  For instance,
%    the van der Corput sequence has the property that all
%    entries lie in [0,1], and that it first eliminates all
%    gaps of size 1/2, then of size 1/4, then 1/8 and so on.
%
%    Perhaps a second reason for looking at plots this way
%    is to notice patterns of convergence in a sequence.
%
%  Example:
%
%    The sequence 1, 3, 5, 7, 6, 2, 4 would plot as
%
%       *---
%        *--
%      *----
%          *
%     *-----
%         *-
%    *------
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of a text file containing
%    the sequence.
%
%    Input, character DIRECTION, defines the direction of the sequence.
%    'R', goes to the right;
%    'L', goes to the left;
%    'U', goes up;
%    'D', goes down.
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SEQUENCE_STREAK_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );

  y = load ( file_name );

  n = length ( y );
  x(1:n) = 1:n;

  if ( direction == 'R' | direction == 'r' )

    plot ( x, y, '*' );

    for i = 1 : n
      line ( [x(i), n], [y(i), y(i) ] )
    end

    xlabel ( 'Sequence Direction--->', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16  )

    ylabel ( 'Sequence Value', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16, 'Rotation', 90 );

    title ( 'Sequence "Streak" Plot', ...
      'FontName', 'Helvetica', 'FontWeight', ...
      'bold', 'FontSize', 16 );

  elseif ( direction == 'L' | direction == 'l' )

    plot ( x, y, '*' );

    for i = 1 : n
      line ( [1, x(i)], [y(i), y(i) ] )
    end

    xlabel ( 'Sequence Direction--->', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16  )

    ylabel ( 'Sequence Value', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16, 'Rotation', 90 );

    title ( 'Sequence "Streak" Plot', ...
      'FontName', 'Helvetica', 'FontWeight', ...
      'bold', 'FontSize', 16 );

  elseif ( direction == 'U' | direction == 'u' )

    plot ( y, x, '*' );

    for i = 1 : n
      line ( [y(i), y(i) ], [x(i), n] )
    end

    xlabel ( 'Sequence Value', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16  )

    ylabel ( 'Sequence Direction--->', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16, 'Rotation', 90 );

    title ( 'Sequence "Streak" Plot', ...
      'FontName', 'Helvetica', 'FontWeight', ...
      'bold', 'FontSize', 16 );

  elseif ( direction == 'D' | direction == 'd' )

    plot ( y, x, '*' );

    for i = 1 : n
      line ( [y(i), y(i) ], [1, x(i)] )
    end

    xlabel ( 'Sequence Value', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16  )

    ylabel ( 'Sequence Direction--->', ...
      'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16, 'Rotation', 90 );

    title ( 'Sequence Streak Plot', ...
      'FontName', 'Helvetica', 'FontWeight', ...
      'bold', 'FontSize', 16 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SEQUENCE_STREAK_DISPLAY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1,'\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Estimation_X] = notREKF_propagate(Estimation_X, OdometryFromThis2Next, Sigma_ODO )

% retrieve v and w
v = OdometryFromThis2Next(1:3);
w = OdometryFromThis2Next(4:6);

NumberOfLandmarks = size(Estimation_X.landmarks, 2);
Jrw = jaco_r(-w);
ExpMinusM = so3_exp(-w);

temp = repmat({Estimation_X.orientation}, NumberOfLandmarks+2,1 );
A = blkdiag(temp{:});
 A(4:6,1:3)=skew(Estimation_X.position)*Estimation_X.orientation;

if NumberOfLandmarks>0
   for i=1:NumberOfLandmarks
    A(6+3*i-2:6+3*i,1:3)=skew(Estimation_X.landmarks(1:3,i))*Estimation_X.orientation;
   end
end

B1=[-Jrw  zeros(3,3); -skew(v)*Jrw -eye(3)];
B=[B1; sparse(3*NumberOfLandmarks,6)];
B=sparse(B);

adA=A*B;

odoCov=diag([w.^2;v.^2])*Sigma_ODO^2;

% final update the covariance
Estimation_X.cov = Estimation_X.cov+ adA*odoCov*adA';

% update position and orientation
Estimation_X.position = Estimation_X.position+Estimation_X.orientation*v;
Estimation_X.orientation = Estimation_X.orientation*so3_exp(w);

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function [yOut, amp, phi]=ApplyHants(y,nb,nf,fet,dod,HiLo,low,high,delta)
if (max(size(size(y)))~=3)
    error('Input data must be three dimensional [time,lat,lon]')
end

[ni ny,nx]=size(y);

yOut= zeros(ni,ny,nx,'single');
amp = zeros(nf+1,ny,nx,'single');
phi = zeros(nf+1,ny,nx,'single');
ts=1:ni;

h= waitbar(0,'Total Progress:');
WBarOuterPosition=get(h,'OuterPosition');
WBarOuterPosition(2)=WBarOuterPosition(2)-WBarOuterPosition(4);
h2=waitbar(0,'Calculating, please wait ...');
set(h2,'OuterPosition',WBarOuterPosition);
for Sample=1:nx
    waitbar(Sample/nx,h);
    for Line=1:ny
        waitbar(Line/ny,h2,['Line:' num2str(Line) ', Sample:' num2str(Sample)]);
        data=y(:,Line,Sample);
        if (sum(isnan(data))~=ni)
            data(isnan(data))=low-1.0;
            [amp(:,Line,Sample),phi(:,Line,Sample),yOut(:,Line,Sample)] ...
                    = HANTS(ni,nb,nf,data,ts,HiLo,low,high,fet,dod,delta);
        end
    end
end
close(h);
close(h2);
end

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function performPCA(confs, objs)

numSbj = length(objs);
groupID = confs.GroupID;

% set parameters
fvecs = [];	% fvecs is accumulator of coefficients from each species.
keep = numSbj-1;	% number of eigenvectors to keep

% do processing of each file
for i = 1:numSbj
	% read in file.
	file = objs{i};
	load(file);
	[path, name, ext] = fileparts(file);
	clear('faces', 'vertices', 'sph_verts');

	[nrows ncols] = size(fvec);
	temp = fvec.';
	temp = reshape(temp, (nrows*3), 1);
	fvecs = [fvecs temp];
end
clear('fvec', 'temp', 'name');

fvecs = fvecs.';
[eigenvecs, scores, eigenvals] = princomp(fvecs);
% trace_eig = trace(eigenvals);
perc_variance_explained = eigenvals/sum(eigenvals);

% drop cells that are to not be kept
eigenvecs = eigenvecs(:,1:keep);
%scores = scores(:,1:keep);
eigenvals = eigenvals(1:keep,1);
perc_variance_explained = perc_variance_explained(1:keep,:);
cum_percent_explained = cumsum(perc_variance_explained);

new_name = [confs.OutDirectory '/' confs.OutputName];
if exist(new_name,'file')
    prompt = {'Enter new filename:'};
    dlg_title = 'New File Name';
    num_lines = 1;
    def = {new_name};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    new_name = answer{1};
end

save(new_name, 'fvecs', 'eigenvecs', 'eigenvals', 'perc_variance_explained', ...
    'cum_percent_explained', 'scores','groupID');

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function s = spmatrix(a)
%SPMATRIX Converts a two-way sparse tensor to sparse matrix.
%
%   SPMATRIX(X) converts a sparse tensor to a sparse matrix. The sparse
%   tensor must be two-dimensional.
%
%   See also SPTENSOR, SPTENSOR/RESHAPE, SPTENMAT
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ndims(a) ~= 2
    error('Sparse tensor must be two dimensional.');
end

if isempty(a.subs)
    s = sparse(a.size(1), a.size(2));
else
    s = sparse(a.subs(:,1), a.subs(:,2), a.vals, a.size(1), a.size(2));
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [eparams, errors] = calibrateEdgeClassifier(efeatures, adjlist, imsegs, eclassifier, ncv)
% [eparams, errors] = calibrateEdgeClassifier(efeatures, adjlist, imsegs,
% eclassifier, ncv)

nimages = numel(imsegs);
for k = 1:ncv
    testind = [(k-1)*nimages/ncv+1:k*nimages/ncv];
    trainind = setdiff([1:nimages], testind);
    [edata{k}, elab{k}] = formatData(efeatures(trainind), adjlist(trainind), {imsegs(trainind).labels});
    econf{k} = test_boosted_dt_mc(eclassifier, edata{k});
%    econf{k} = 1 ./ (1+exp(-econf{k}));
end

for k = 1:ncv
    disp(['iter: ' num2str(k)])
    traink = setdiff([1:ncv], k);
    eparams{k} = fminunc(@(x) objective(x, cat(1, econf{traink}), cat(1, elab{traink})), [-1 0], optimset('TolFun', 0.001));
end

for k = 1:ncv
    econf{k} = 1 ./ (1+exp(eparams{k}(1)*econf{k}+eparams{k}(2)));
end

elab = cat(1, elab{:});
econf = cat(1, econf{:});

eerror = mean((econf>0.5)~=elab);

econf2 = 1-abs(elab-econf);

ind1 = find(elab==0);
ind2 = find(elab==1);
px = [0.025:0.05:0.975];
f1 = ksdensity(econf(ind1), px, 'support', [0 1]);
f2 = ksdensity(econf(ind2), px, 'support', [0 1]);
fc = ksdensity(econf2, px, 'support', [0 1]);
%fc = fc;

errors.err = eerror;
errors.pneg = f1;
errors.ppos = f2;
errors.conf = fc;
errors.px = px;

medFS = 18;
bigFS = 20;

figure(1), hold on, plot(px, fc, 'y', 'LineWidth', 2);
%axis([0 1 0 1])
xlabel('Confidence in True Label', 'FontSize', medFS)
ylabel('Frequency', 'FontSize', medFS)
title('Same Label Confidence', 'FontSize', bigFS)
set(gca, 'FontSize', medFS)

figure(2), hold on, plot(px, f2 ./ (f1+f2), 'y', 'LineWidth', 2)
hold on, plot(px, px, '--k')
axis([0 1 0 1])
xlabel('Estimated Probability', 'FontSize', medFS)
ylabel('Empirical Probability', 'FontSize', medFS)
%title('Same Label Confidence', 'FontSize', bigFS)
set(gca, 'FontSize', medFS)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function err = objective(param, econf, elab)
econf = 1./(1+exp(param(1)*econf+param(2)));
px = [0.025:0.05:0.975];
f1 = ksdensity(econf(elab==0), px, 'support', [0 1])+1E-10;
f2 = ksdensity(econf(elab==1), px, 'support', [0 1])+1E-10;
f1 = f1 / sum(f1+f2);
f2 = f2 / sum(f1+f2);
err = sum((f1+f2).*(px - f2./(f1+f2)).^2);
disp(num2str([sum((f1+f2).*abs(px - f2./(f1+f2))) param]))

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

function [data, lab] = formatData(features, adjlist, labels)
% concatenate data and select ndata random datapoints

nimages = numel(features);

[tmp, nvars] = size(features{1});

% count edges
ne = 0;
for f = 1:nimages
    ne = ne + size(features{f}, 1);
end

data = zeros(ne, nvars);
lab = zeros(ne, 1);

% concatenate data
c = 0;
for f = 1:nimages
    cf = size(features{f}, 1);
    data(c+1:c+cf, :) = features{f};
    s1 = adjlist{f}(:, 1);
    s2 = adjlist{f}(:, 2);
    lab(c+1:c+cf) = (labels{f}(s1)==labels{f}(s2));
    c = c + cf;
end

author: zouchuhang

repo: LayoutNet


text:
function OUT2 = hillshade(DEM,varargin)

%HILLSHADE create hillshading from a digital elevation model (GRIDobj)
%
% Syntax
%
%     H = hillshade(DEM)
%     H = hillshade(DEM,'pn','pv',...)
%
% Description
%
%     Hillshading is a very powerful tool for relief depiction.
%     hillshade calculates a shaded relief for a digital elevation model
%     based on the angle between the surface and the incoming light beams.
%     If no output arguments are defined, the hillshade matrix will be
%     plotted with a gray colormap. The hillshading algorithm follows the
%     logarithmic approach to shaded relief representation of Katzil and
%     Doytsher (2003).
%
% Input
%
%     DEM       Digital elevation model (class: GRIDobj)
%
% Parameter name/value pairs
%
%     'azimuth'         azimuth angle, (default=315)
%     'altitude'        altitude angle, (default=60)
%     'exaggerate'      elevation exaggeration (default=1). Increase to
%                       pronounce elevation differences in flat terrain
%     'useblockproc'    true or {false}: use block processing
%                       (see function blockproc)
%     'useparallel'     true or {false}: use parallel computing toolbox
%     'blocksize'       blocksize for blockproc (default: 5000)
%
%
% Output
%
%     H         shaded relief (ranges between 0 and 1)
%
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     hillshade(DEM)
%
% References
%
%     Katzil, Y., Doytsher, Y. (2003): A logarithmic and sub-pixel approach
%     to shaded relief representation. Computers & Geosciences, 29,
%     1137-1142.
%
% See also: SURFNORM, IMAGESCHS
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 18. August, 2017

% Parse inputs
p = inputParser;
p.StructExpand  = true;
p.KeepUnmatched = false;
p.FunctionName = 'hillshade';
addParamValue(p,'azimuth',315,@(x) isscalar(x) && x>= 0 && x<=360);
addParamValue(p,'altitude',60,@(x) isscalar(x) && x>= 0 && x<=90);
addParamValue(p,'exaggerate',1,@(x) isscalar(x) && x>0);
addParamValue(p,'useparallel',true);
addParamValue(p,'blocksize',2000);
addParamValue(p,'useblockproc',true,@(x) isscalar(x));
parse(p,varargin{:});

OUT     = DEM;
OUT.Z   = [];

cs      = DEM.cellsize;
azimuth = p.Results.azimuth;
altitude = p.Results.altitude;
exaggerate = p.Results.exaggerate;

% Large matrix support. Break calculations in chunks using blockproc
if numel(DEM.Z)>(10001*10001) && p.Results.useblockproc;
    blksiz = bestblk(size(DEM.Z),p.Results.blocksize);
    padval = 'symmetric';
    Z      = DEM.Z;
    % The anonymous function must be defined as a variable: see bug 1157095
    fun   = @(x) hsfun(x,cs,azimuth,altitude,exaggerate);
    HS = blockproc(Z,blksiz,fun,...
                'BorderSize',[1 1],...
                'padmethod',padval,...
                'UseParallel',p.Results.useparallel);
    OUT.Z = HS;
else
    OUT.Z = hsfun(DEM.Z,cs,azimuth,altitude,exaggerate);
end

OUT.name = 'hillshade';
OUT.zunit = '';

if nargout == 0;
    OUT.Z = uint8(OUT.Z*255);
    imagesc(OUT);
    colormap(gray)
else
    OUT2 = OUT;
end

end
%% Subfunction
function H = hsfun(Z,cs,azimuth,altitude,exaggerate)

if isstruct(Z)
    Z = Z.data;
end

% correct azimuth so that angles go clockwise from top
azid = azimuth-90;

% use radians
altsource = altitude/180*pi;
azisource = azid/180*pi;

% calculate solar vector
[sx,sy,sz] = sph2cart(azisource,altsource,1);

% calculate surface normals
[Nx,Ny,Nz] = surfnorm(Z/cs*exaggerate);

% calculate cos(angle)
% H = [Nx(:) Ny(:) Nz(:)]*[sx;sy;sz];
% % reshape
% H = reshape(H,size(Nx));

H = Nx*sx + Ny*sy + Nz*sz;

% % usual GIS approach
% H = acos(H);
% % force H to range between 0 and 1
% H = H-min(H(:));
% H = H/max(H(:));

end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function vocabList = getVocabList()
%GETVOCABLIST reads the fixed vocabulary list in vocab.txt and returns a
%cell array of the words
%   vocabList = GETVOCABLIST() reads the fixed vocabulary list in vocab.txt
%   and returns a cell array of the words in vocabList.

%% Read the fixed vocabulary list
fid = fopen('vocab.txt');

% Store all dictionary words in cell array vocab{}
n = 1899;  % Total number of words in the dictionary

% For ease of implementation, we use a struct to map the strings => integers
% In practice, you'll want to use some form of hashmap
vocabList = cell(n, 1);
for i = 1:n
    % Word Index (can ignore since it will be = i)
    fscanf(fid, '%d', 1);
    % Actual Word
    vocabList{i} = fscanf(fid, '%s', 1);
end
fclose(fid);

end

author: Ayatans

repo: Machine-Learning-homework


text:
function []=v2structDemo2()
%  V2STRUCTDEMO2 is a usage example demo function for the function V2STRUCT. This demo
%  shows a suggestion for how to use v2struct in managing input to other functions with
%  improved usability.
%
%% Last update:
%    30.06.2011, Adi N.

fprintf(['Running usage demo function for vs2truct function:\n',...
         '\nThis demo shows a suggestion for how to use v2struct in managing input to',...
         ' other functions\nwith improved usability. Please read comments in the .m file for explanations\n',...
         'Press any key between breaks.\n']);
pause()

fprintf('\nOutput from four statements:\n');
pause()
% See the differences between calling foo and fooSmart using v2struct
% Note that we expect the same results in the following statements

% regular call to foo
foo(1,2,3,4);

% pack to structure and call with fooSmart
S = v2struct(1,2,3,4,{'fieldNames','num1','denum1','num2','denum2'});
fooSmart(S);

% regular call to fooSmartParamOrder
fooSmartParamOrder(1,2,3,4);

% call with structure to fooSmartParamOrder
fooSmartParamOrder(S);

% Consider the following functions

% foo is a function with regular construction using arguments passed in a defined order
% the result is that you have to call it with 4 arguments, and you must remember their order.
 function isFirstBiggerThanSecond = foo(num1,denum1,num2,denum2)
 isFirstBiggerThanSecond = (num1/denum1>num2/denum2);
 fprintf('\n using foo with a regular call:\n')
 if isFirstBiggerThanSecond
    fprintf('  %d/%d is bigger than %d/%d\n',num1,denum1,num2,denum2);
 else
    fprintf('  %d/%d is not bigger than %d/%d\n',num1,denum1,num2,denum2);
 end

% Now consider the function fooSmart which is modified just slightly to accept input
% packed in a structure. You get the same result but you can use a structure.
    function isFirstBiggerThansSecond = fooSmart(num1,denum1,num2,denum2)
 if isstruct(num1) % unpack if called with a structure
    v2struct(num1);
 end
 isFirstBiggerThansSecond = (num1/denum1>num2/denum2);
 fprintf('\n using fooSmart with a structure:\n')
 if isFirstBiggerThansSecond
    fprintf('  %d/%d is bigger than %d/%d\n',num1,denum1,num2,denum2);
 else
    fprintf('  %d/%d is not bigger than %d/%d\n',num1,denum1,num2,denum2);
 end

% Lastly the function fooSmartParamOrder represent a situation in which the order of the
% parameters is changed in relation to foo. Here calling regularly with the same variables
% as before leads to an erroneous result. However when using the structure as input the
% order is not important and you get the correct result.
     function isFirstBiggerThansSecond = fooSmartParamOrder(num1,num2,denum1,denum2)
 if isstruct(num1) % unpack if called with a structure
    v2struct(num1);
    fprintf('\n using fooSmartParamOrder with structure:\n')
 else
     fprintf('\n using fooSmartParamOrder with same variables:\n')
 end
 isFirstBiggerThansSecond = (num1/denum1>num2/denum2);
 if isFirstBiggerThansSecond
    fprintf('  %d/%d is bigger than %d/%d\n',num1,denum1,num2,denum2);
 else
    fprintf('  %d/%d is not bigger than %d/%d\n',num1,denum1,num2,denum2);
 end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function FPDF = GammaGenPDF (a)
% Return function handles to routines to calculate the area, mean,
% and second moment of a unit-variance, zero-mean, generalized
% Gamma probability density function with parameter a.
%                 b
%   Farea(a,b) = Int p(x) dx
%                 a
%                 b
%   Fmean(a,b) = Int x p(x) dx
%                 a
%                b
%   Fvar(a,b) = Int x^2 p(x) dx
%                a
% where p(x) = b/(2G(a) exp(-b|x|) (b|x|)^(a-1)
% with b=sqrt(a(a+1)) and where G(a) is the complete gamma function.

% global GGenPar
GGenPar = a;   % GGenPar available to nested functions

FPDF = {@GGenarea, @GGenmean, @GGenvar};

return

%----- ----- begin nested functions
function v = GGenarea (a, b)

% Evaluate the function so as to avoid taking differences
% between nearly equal quantities (e.g. when a<0 and b<0)
if (b >= 0)
  if (a >= 0)
    v = F0(b) - F0(a);                    % Both a and b positive
  else
    v = (F0(b) + 0.5) + (F0(-a) + 0.5);   % a negative, b positive
  end
else
  if (a < 0)
    v = F0(-a) - F0(-b);                  % Both a and b negative
  else
    v = (-0.5 - F0(a)) + (-0.5 - F0(-b)); % a positive, b negative
  end
end

return
end

% -----
% Evaluate the indefinite integral
% The definite integral between a and b is F0(b) - F0(a)

function v = F0(x)

% global GGenPar

v = -Gamma2aCCDF(x, GGenPar);

return
end

% ----- ------
function v = GGenmean (a, b)

% Since the integrand x*p(x) is odd, Gmean(A,B)=Gmean(abs(A),abs(B))
v = F1(abs(b)) - F1(abs(a));

return
end

% -----
% Evaluate the indefinite integral
% The definite integral between a and b is F1(b) - F1(a)

function v = F1 (x)

% global GGenPar

n = 1;
a1 = GGenPar;
a2 = a1 + n;
b1 = sqrt(a1*(a1+1));
b2 = sqrt(a2*(a2+1));
v = -gamma(a2)/(gamma(a1)*b1^n) * Gamma2aCCDF(b1*abs(x)/b2, a2);
% Simplifications for n=1,
%   G(a2)=(a1+1) G(a1)
% Then
%   G(a2)/(G(a1) b1) = sqrt((a1+1)/a1)
%              b1/b2 = sqrt(a1/(a1+1))

return
end

% ----- ------
function v = GGenvar (a, b)

% Evaluate the function so as to avoid taking differences
% between nearly equal quantities (e.g. when a<0 and b<0)
if (b >= 0)
  if (a >= 0)
    v = F2(b) - F2(a);                    % Both a and b positive
  else
    v = (F2(b) + 0.5) + (F2(-a) + 0.5);   % a negative, b positive
  end
else
  if (a < 0)
    v = F2(-a) - F2(-b);                  % Both a and b negative
  else
    v = (-0.5 - F2(a)) + (-0.5 - F2(-b)); % a positive, b negative
  end
end

return
end

% -----
% Evaluate the indefinite integral
% The definite integral between a and b is F2(b) - F2(a)

function v = F2 (x)

% global GGenPar

n = 2;
a1 = GGenPar;
a2 = a1 + n;
b1 = sqrt(a1*(a1+1));
b2 = sqrt(a2*(a2+1));
v = -gamma(a2)/(gamma(a1)*b1^n) * Gamma2aCCDF(b1*x/b2, a2);
% Simplifications for n=2,
%   G(a2)=(a1+2)(a1+1) G(a1)
% Then
%   G(a2)/(G(a1) b1^2) = (a1+2)/a1
%                b1/b2 = sqrt(a1(a1+1)/((a1+2)(a1+3))
return
end

% ---- end of the nested functions
end


text:
function res = clipPolyline(poly, box)
%CLIPPOLYLINE Clip an open polyline with a rectangular box.
%
%   POLY2 = clipPolyline(POLY, BOX);
%   POLY is N-by-2 array of vertex coordinates.
%   BOX has the form: [XMIN XMAX YMIN YMAX].
%   Returns the set of polylines created by the intersection of the
%   polyline POLY and the bounding box BOX. The result is a cell array with
%   as many cells as the number of curve clips.
%
%
%   Example
%     circle = [5 5 6];
%     poly = circleToPolygon(circle, 200);
%     box = [0 10 0 10];
%     res = clipPolyline(poly, box);
%     figure;
%     hold on; axis equal; axis([-2 12 -2 12]);
%     drawCircle(circle, 'b:')
%     drawBox(box, 'k')
%     drawPolyline(res, 'linewidth', 2, 'color', 'b')
%
%   See also
%     polygons2d, boxes2d, clipPolygon, clipEdge
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-05-14
% Copyright 2005-2022 INRA - Cepia Software Platform

% check case of polylines stored in cell array
if iscell(poly)
    res = cell(1, length(poly));
    for i = 1:length(poly)
        res{i} = clipPolyline(poly{i}, box);
    end
    return;
end

% check case of empty polylines
N = size(poly, 1);
if N == 0
    res = cell(0, 0);
    return
end

% create edges array of polyline
edges = [poly(1:N-1, :) poly(2:N, :)];

% clip edges
edges = clipEdge(edges, box);

% select non empty edges, and get their vertices
% find clipped edges within box
inds = sum(abs(edges), 2) > 1e-14;

% find list of adjacent edges within box
dinds = diff(inds);
inds0 = find(dinds == 1) + 1;
if inds(1) == 1
    inds0 = [1 inds0];
end
inds1 = find(dinds == -1);
if inds(end) == 1
    inds1 = [inds1 N-1];
end

nClips = length(inds0);
res = cell(1, nClips);
for iClip = 1:nClips
    range = inds0(iClip):inds1(iClip);
    poly2 = [edges(range, 1:2) ; edges(range(end), 3:4)];
    res{iClip} = poly2;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function determ = fibonacci2_determinant ( n )

%*****************************************************************************80
%
%% FIBONACCI2_DETERMINANT returns the determinant of the FIBONACCI2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  if ( n == 1 )
    determ = 0.0;
  else
    determ = -1.0;
  end

  return
end


text:
function vec = readBinaryDoubles(fileName, format)

% READBINARYDOUBLES Read information from a binary file in as doubles.
% FORMAT
% DESC reads in information from a binary file as a vector of doubles.
% ARG fileName : the name of the file.
% ARG format : the file format for 'fopen', defaults to 'ieee-le'.
% RETURN vec : vector of values in the file.
%
% COPYRIGHT : Neil D. Lawrence, 2009
%
% SEEALSO : fopen, fread, fclose

% NDLUTIL

  if nargin < 2
    format = 'ieee-le';
  end
  fid = fopen(fileName, 'r', format);
  vec = fread(fid, inf, 'double')';
  fclose(fid);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [y] = ge(x, val)

y = cell(size(x));
for k = 1:numel(x)
  y{k} = x{k}>=val;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function demo(IP)
%while(1)
    [plots,~,~]=ad9361_LTE(IP,'LTE1.4');
    close all;

    delete(plots{1});delete(plots{2});delete(plots{3});delete(plots{4});
    clear functions;
    [plots,~,~]=ad9361_LTE(IP,'LTE3');
    close all;

    delete(plots{1});delete(plots{2});delete(plots{3});delete(plots{4});
    clear functions;
    [plots,~,~]=ad9361_LTE(IP,'LTE5');
    close all;

    delete(plots{1});delete(plots{2});delete(plots{3});delete(plots{4});
    clear functions;
    [plots,~,~]=ad9361_LTE(IP,'LTE10');
    close all;

    delete(plots{1});delete(plots{2});delete(plots{3});delete(plots{4});
    clear functions;
%end

author: analogdevicesinc


text:
function S = slsum(A, d)
%SLSUM Compute the sum of values in subarrays
%
% $ Syntax $
%   - S = slsum(A)
%   - S = slsum(A, d)
%   - S = slsum(A, [d1 d2 ... dk])
%
% $ Arguments $
%   - A:        the input array
%   - d:        the dimensions along which the sum is performed
%   - S:        the resultant sum matrix
%
% $ Description $
%   - S = slsum(A) sums up the values in column vectors of A. It is
%     equivalent to S = sum(A)
%
%   - S = slsum(A, d) sums up the values along dimension d. It is
%     equivalent to S = sum(A, d)
%
%   - S = slsum(A, [d1 d2 ... dk]) sums up the values along dimension
%     d1, d2, ... dk.
%
% $ History $
%   - Created by Dahua Lin on Nov 18th, 2005
%

%% parse and verify input arguments
if nargin < 2 || isempty(d)
    d = 1;
end

%% compute
if isscalar(d)
    S = sum(A, d);
else
    k = length(d);
    S = A;
    for i = 1 : k
        S = sum(S, d(i));
    end
end

author: lmthang


function bool = versionCompare( vA, operator, vB )
%VERSIONCOMPARE Performs a version comparison operation
    switch operator
        case '<'
            bool = isVersionBelow(vA, vB);
        case '>'
            bool = isVersionBelow(vB, vA);
        case {'<=', '=<'}
            bool = ~isVersionBelow(vB, vA);
        case {'>=', '=>'}
            bool = ~isVersionBelow(vA, vB);
        case {'=', '=='}
            bool = ~isVersionBelow(vA, vB) && ~isVersionBelow(vB, vA);
        case {'~=', '!='}
            bool = isVersionBelow(vA, vB) || isVersionBelow(vB, vA);
        otherwise
            error('versionCompare:UnknownOperator',...
                  '"%s" is not a known comparison operator', operator);
    end
end


text:
function halton_dataset ( )

%*****************************************************************************80
%
%% HALTON_DATASET generates a Halton dataset and writes it to a file.
%
%  Discussion:
%
%    This program is meant to be used interactively.  It's also
%    possible to prepare a simple input file beforehand and use it
%    in batch mode.
%
%    The program requests input values from the user:
%
%    * NDIM, the spatial dimension,
%    * N, the number of points to generate,
%    * STEP, the index of the first subsequence element to be computed.
%    * SEED(1:NDIM), the Halton sequence index corresponding
%      to STEP = 0.
%    * LEAP(1:NDIM), the successive jumps in the Halton sequence.
%    * BASE(1:NDIM), the Halton bases (usually distinct primes).
%
%    The program generates the data, writes it to the file
%
%      halton_NDIM_N.txt
%
%    where "NDIM" and "N" are the numeric values specified by the user,
%    and then asks the user for more input.   To indicate that no further
%    computations are desired, it is enough to input a nonsensical
%    value, such as -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HALTON_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generate a Halton dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program is meant to be used interactively.\n' );
  fprintf ( 1, '  It is also possible to prepare a simple input \n' );
  fprintf ( 1, '  file beforehand and use it in batch mode.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * NDIM, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * STEP, the index of the first subsequence element.\n' );
  fprintf ( 1, '  * SEED(1:NDIM), the Halton sequence element\n' );
  fprintf ( 1, '    corresponding to STEP = 0\n' );
  fprintf ( 1, '  * LEAP(1:NDIM), the successive jumps in the\n' );
  fprintf ( 1, '    Halton sequence.\n' );
  fprintf ( 1, '  * BASE(1:M), the Halton bases,\n' );
  fprintf ( 1, '    usually distinct primes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program generates the data, writes it to the file\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    halton_NDIM_N.txt\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where "NDIM" and "N" are the numeric values specified\n' );
  fprintf ( 1, '  by the user, and then asks the user for more input.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To indicate that no further computations are \n' );
  fprintf ( 1, '  desired, it is enough to input a nonsensical value, \n' );
  fprintf ( 1, '  such as -1.\n' );

  while ( 1 )

    fprintf ( 1, '  *\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '*  Ready to generate a new dataset:\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '  *\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NDIM is the spatial dimension.\n' );
    fprintf ( 1, '  (Try ''2'' if you have no preference.)\n' );
    fprintf ( 1, '  Any value less than 1 terminates execution.\n' );
    ndim = input ( '  Enter NDIM:  ' );

    fprintf ( 1, '  User input NDIM = %d\n', ndim );

    if ( ~halham_ndim_check ( ndim ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HALTON_DATASET\n' );
      fprintf ( 1, '  The input value of NDIM = %d\n', ndim );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  N is the number of points.\n' );
    fprintf ( 1, '  (Try ''25'' if you have no preference.)\n' );
    fprintf ( 1, '  Any value less than 1 terminates execution.\n' );
    n = input ( '  Enter N:  ' );

    fprintf ( 1, '  User input N = %d\n', n );

    if ( ~halham_n_check ( n ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HALTON_DATASET\n' );
      fprintf ( 1, '  The input value of N = %d\n', n );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  STEP is the index of the first subsequence element.\n' );
    fprintf ( 1, '  (Try ''0'' or ''1'' if you have no preference.)\n' );
    fprintf ( 1, '  Any value less than 0 terminates execution.\n' );
    step = input ( '  Enter STEP:  ' );

    fprintf ( 1, '  User input STEP = %d\n', step );

    if ( ~halham_step_check ( step ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HALTON_DATASET\n' );
      fprintf ( 1, '  The input value of STEP = %d\n', step );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SEED(1:NDIM) is the starting element index\n' );
    fprintf ( 1, '  for each coordinate.\n' );
    fprintf ( 1, '  (Try ''[0,0,...,0]'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 0 terminates execution.)\n' );
    seed = [];
    seed(1:ndim) = input ( '  Enter SEED(1:NDIM):  ' );

    fprintf ( 1, '\n' );
    i4vec_transpose_print ( ndim, seed, 'SEED:' );

    if ( ~halham_seed_check ( ndim, seed ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HALTON_DATASET\n' );
      fprintf ( 1, '  The negative input value of at least one entry of\n' );
      fprintf ( 1, '  SEED is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  LEAP(1:NDIM) is the leaping multiplier\n' );
    fprintf ( 1, '  for each coordinate.\n' );
    fprintf ( 1, '  (Try ''[1,1,...,1]'' if you have no preference.)\n' );
    fprintf ( 1, '  (Another choice is a prime bigger than all the bases.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    leap = [];
    leap(1:ndim) = input ( '  Enter LEAP(1:NDIM):  ' );

    fprintf ( 1, '\n' );
    i4vec_transpose_print ( ndim, leap, 'LEAP:' );

    if ( ~halham_leap_check ( ndim, leap ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HALTON_DATASET\n' );
      fprintf ( 1, '  The input value of at least one entry of\n' );
      fprintf ( 1, '  LEAP is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  BASE(1:NDIM) is the base for each coordinate,\n' );
    fprintf ( 1, '  usually distinct primes.\n' );
    fprintf ( 1, '  (Try ''[2,3,5,7,11,13,...]'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 2 terminates execution.)\n' );
    base = [];
    base(1:ndim) = input ( '  Enter BASE:  ' );

    fprintf ( 1, '\n' );
    i4vec_transpose_print ( ndim, base, 'BASE:' );

    if ( ~halton_base_check ( ndim, base ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HALTON_DATASET\n' );
      fprintf ( 1, '  The input value of at least one entry of\n' );
      fprintf ( 1, '  BASE is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
    end

    r = i4_to_halton_sequence ( ndim, n, step, seed, leap, base );

    file_out_name = ...
      strcat ( 'halton_', num2str ( ndim ), '_', num2str ( n ), '.txt' );

    halham_write ( ndim, n, step, seed, leap, base, r, file_out_name );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  The data was written to the file "%s".\n', ...
      file_out_name );

  end

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function cvx_optval = sum_smallest( x, varargin )

%SUM_SMALLEST Sum of the smallest k elements of a vector.
%   For a real vector X and an integer k between 1 and length(X) inclusive,
%   y = SUM_SMALLEST(X,k) is the sum of the k smallest elements of X; e.g.,
%       temp = sort( x )
%       y = sum( temp( 1 : k ) )
%   If k=1, then SUM_SMALLEST(X,k) is equivalent to MIN(X); if k=length(X),
%   then SUM_SMALLEST(X,k) is equivalent to SUM(X).
%
%   Both X and k must be real, and k must be a scalar. But k is not, in
%   fact, constrained to be an integer between 1 and length(X); the
%   function is extended continuously and logically to all real k. For
%   example, if k <= 0, then SUM_SMALLEST(X,k)=0. If k > length(X), then
%   SUM_SMALLEST(X,k)=SUM(X). Non-integer values of k interpolate linearly
%   between their integral neighbors.
%
%   For matrices, SUM_SMALLEST(X,k) is a row vector containing the
%   application of SUM_SMALLEST to each column. For N-D arrays, the
%   SUM_SMALLEST operation is applied to the first non-singleton dimension
%   of X.
%
%   SUM_SMALLEST(X,k,DIM) performs the operation along dimension DIM of X.
%
%   Disciplined convex programming information:
%       SUM_SMALLEST(X,...) is concave and nondecreasing in X. Thus, when
%       used in CVX expressions, X must be concave (or affine). k and DIM
%       must both be constant.

error( nargchk( 2, 3, nargin ) );
cvx_optval = -sum_largest( -x, varargin{:} );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function draw_shape(x, y, color)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw face given x, y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% inverse y for better visual

bbox = getbbox([x y]);

radius = bbox(4) / 20;

if 1
    %hold on;

    n = size(x,1);

    for i = 1 : n

        r = x(i) - radius/2;
        c = y(i) - radius/2;

        rectangle('Position',[r,c,radius,radius],'Curvature',[1,1],...
            'FaceColor',color);

    end

    %plot(x, y, color);
    %plot(x, y, [color '.']);

    %hold off;
end

author: tntrung


text:
classdef SymBuilder < handle
%SYMBUILDER  Create a SYMBUILDER object for creating symbolic optimization problems
%
%   B = SymBuilder() creates a blank optimization problem
%
%   See also SymBuilder.AddObj SymBuilder.AddCon SymBuilder.Build
%
%   Copyright (C) 2012 Jonathan Currie (www.inverseproblem.co.nz)

    properties(SetAccess=private)
        vars        %Symbolic array of symbolic variables
        sobj        %Symbolic matrix of equations
        jac         %Symbolic jacobian of equations
        hess        %Structure of hessians of each equation
        hesslag     %Symbolic hessian of the lagrangian
        Opt         %OPTI object
        lastSolveOpts %Last Options passed when Solve() was called
    end

    properties(SetAccess=private)%,GetAccess=private)
        eqs         %Equations symbolic array
        cl          %Constraint Lower Bounds
        cu          %Constraint Upper Bounds
        constnt     %Cell array of constants
        exprsn      %Cell array of expressions
        bnds        %Cell array of bounds
        lb          %Double array of lower bounds
        ub          %Double array of upper bounds
        vartypes    %Cell array of variable types
        xtype       %Char array of variable types
        indobj      %Index of objectives in equations
        objlin      %Index of linear, quadratic and nonlinear objectives
        conlin      %Index of linear, quadratic and nonlinear constraints
        objbias     %Problem bias terms (for linear and quadratic problems)
        noObjs      %Counter of number objectives
        noCons      %Counter of number of constraints
        bldstat     %Build Status
        bldtime     %Build Time
        resgrp      %Cell array of result groups
        resexp      %Cell array of result expressions
        verbose     %Verbosity mode
    end

    methods
        function B = SymBuilder(verbose)
            if(nargin < 1), verbose = true; end
            if (~exist('syms.m','file')), error('The Symbolic Math Toolbox must be installed to use SymBuilder!'); end
            %Initialization
            B.noObjs = 0;
            B.noCons = 0;
            B.verbose = verbose;
            Unbuild(B);
        end

        function B = AddObj(B,str)
        %Add a General Objective (Linear or Nonlinear)

            %Check for string
            if(ischar(str))
                digits(16);
                wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
                s = sym(str);
                warning(wstate);
            elseif(isa(str,'sym'))
                s = str;
            else
                error('Unknown objective form');
            end

            %Concatenate with existing equations
            B.eqs = [B.eqs;s];
            B.cl = [B.cl;NaN];
            B.cu = [B.cu;NaN];
            if(isempty(B.indobj))
                B.indobj = length(B.cl);
            else
                B.indobj = [B.indobj;length(B.cl)];
            end
            %Add no objs
            B.noObjs = B.noObjs + 1;
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddCon(B,str_in,cl,cu)
        %Add a General Constraint (Linear or Nonlinear)

            %Check for string
            if(ischar(str_in))
                %Parse constraint string
                [str,l,u] = SymBuilder.parseConstraint(str_in);
                digits(16);
                wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
                s = sym(str);
                warning(wstate);
            elseif(isa(str_in,'sym'))
                if(nargin < 4 || ~isnumeric(cl) || ~isnumeric(cu) || length(cl)~=length(cu) || length(cl)~=length(str_in))
                    error('When supplying a symbolic expression to AddCon, you must supply as AddCon(sym_con,cl,cu) with cl, cu as numeric vectors of the same length.');
                end
                if(size(str_in,2) > 1)
                    s = str_in.';
                else
                    s = str_in;
                end
                l = cl; u = cu;
            else
                 error('Unknown constraint form');
            end
            %Concatenate with existing equations if not empty
            if(~isempty(s))
                B.eqs = [B.eqs;s];
                B.cl = [B.cl;l];
                B.cu = [B.cu;u];
                %Add no cons
                B.noCons = B.noCons + length(l);
                %Indicate Rebuild Required
                Unbuild(B);
            else
                optiwarn('symb:emptycon','The following constraint will be ignored by SymBuilder as both sides appear to be numeric.\n  %s\n',str_in);
            end
        end

        function B = AddLinCon(B,x,A,rl,ru)
        %Add a Linear Constraint

            if(nargin < 5 || ~isnumeric(rl) || ~isnumeric(ru) || length(rl)~=length(ru) || length(rl)~=size(A,1))
                error('When supplying linear constraints, you must supply as AddLinCon(x,A,rl,ru) with rl, ru as numeric vectors of the same length.');
            elseif(~isa(x,'sym') || length(x) ~= size(A,2))
                error('x must be a symbolic variable vector = size(A,2)');
            end
            %Evaluate Matrix & Store
            B.eqs = [B.eqs;A*x];
            B.cl = [B.cl;rl];
            B.cu = [B.cu;ru];
            %Add no cons
            B.noCons = B.noCons + length(rl);
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddConstraint(B,str)
        %Same as AddCon()
            B = AddCon(B,str);
        end

        function B = AddConstant(B,str,val,varargin)
        %Identify variable as a constant

            %Concatenate with existing constants
            B.constnt = [B.constnt; {str val}];
            %If user has passed multiple constants, process all
            if(nargin > 3 && ~isempty(varargin))
                if(mod(length(varargin),2))
                    error('You must supply a name and value for each constant');
                else
                    for i = 1:2:length(varargin)
                        B.constnt = [B.constnt; varargin(i:i+1)];
                    end
                end
            end
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddExpression(B,str,grp,ex)
        %Identify variable as an expression

            %Parse expression string
            [str,exp] = SymBuilder.parseExpression(str);
            %Concatenate with existing expressions
            B.exprsn = [B.exprsn; {str exp}];
            %If group passed, add to result group
            if(nargin > 3), B = AddResultExp(B,grp,exp,ex);
            elseif(nargin > 2), B = AddResultExp(B,grp,exp); end
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddBound(B,str)
        %Add variable bound

            %Parse expression string
            [var,llb,lub] = SymBuilder.parseBound(str);
            %Concatenate with existing bounds
            B.bnds = [B.bnds; {var llb lub}];
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddBounds(B,x,lb,ub)
        %Add numerical bounds

            bds = cell(length(x),3);
            for i = 1:length(x)
                bds{i,1} = char(x(i));
                bds{i,2} = sprintf('%1.15g',lb(i));
                bds{i,3} = sprintf('%1.15g',ub(i));
            end
            %Concatenate with existing bounds
            B.bnds = [B.bnds; bds];
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddInteger(B,str)
        %Add integer constraint

            %Parse expression string
            [var,xx] = SymBuilder.parseInteger(str);
            %Concatenate with existing constants
            B.vartypes = [B.vartypes; {var xx}];
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddIntegers(B,x,str)
        %Add integer constraints with character string

            xx = cell(length(x),2);
            for i = 1:length(x)
                xx{i,1} = char(x(i));
                xx{i,2} = str(i);
            end
            %Concatenate with existing constants
            B.vartypes = [B.vartypes; xx];
            %Indicate Rebuild Required
            Unbuild(B);
        end

        function B = AddResultGroup(B,name,str)
        %Add Result Group
            %Concatenate with existing groups
            B.resgrp = [B.resgrp; {name str}];
        end

        function B = AddResultExp(B,name,str,bin)
        %Add Result Expression
            %Parse expression string
            [group,name] = SymBuilder.parseResExp(name);
            %Concatenate with existing groups
            if(nargin > 3)
                B.resexp = [B.resexp; {group name str bin}];
            else
                B.resexp = [B.resexp; {group name str []}];
            end
        end

        function B = Draft(B)
        %Draft the object by solving for the system Jacobian
            t = tic;
            %Build Symbolic representation of the equations
            if(B.verbose), fprintf('\nGenerating Symbolic Representation of Equations...'); end
            buildSymRep(B);
            if(B.verbose), fprintf('Done\n'); end
            %Generate Equations Jacobian
            if(B.verbose), fprintf('Generating Symbolic Jacobian...'); end
            buildJac(B);
            if(B.verbose), fprintf('Done\n'); end
            %Determine Linear & Nonlinear Equations
            if(B.verbose), fprintf('Generating Equation Linearity...'); end
            detLinearity(B);
            if(B.verbose), fprintf('Done\n'); end
            %Save build time
            B.bldtime = toc(t);
            %Assign build Status
            B.bldstat = 'draft';
        end

        function B = Build(B)
        %Build the object by solving for the system Hessian
            t = tic;
            %Build Symbolic representation of the equations
            if(B.verbose), fprintf('\nGenerating Symbolic Representation of Equations...'); end
            buildSymRep(B);
            if(B.verbose), fprintf('Done\n'); end
            %Generate Equations Jacobian
            if(B.verbose), fprintf('Generating Symbolic Jacobian...'); end
            buildJac(B);
            if(B.verbose), fprintf('Done\n'); end
            %Generate Equations Hessians (also determines linearity)
            if(B.verbose), fprintf('Generating Symbolic Hessian...'); end
            buildHess(B);
            if(B.verbose), fprintf('Done\n'); end
            %Save build time
            B.bldtime = toc(t);
            %Assign build Status
            B.bldstat = 'built';
        end

        function B = Unbuild(B)
        %Unbuild object (free memory)

            B.jac = [];
            B.hess = [];
            B.hesslag = [];
            B.bldstat = 'unbuilt';
            B.bldtime = 0;
            B.Opt = [];

            clear symb_cb
        end

        function prob = GetLinProb(B)
        %Get Linear Problem Matrices [min f'*x s.t. rl <= A*x <= ru]

            %Check is built
            IsBuilt(B,1);
            prob = optiprob;
            %Get Linear Objective Index
            ind = B.objlin == 1;
            if(sum(ind) == 0)
                error('This object does not contain a linear objective');
            end
            prob.f = double(B.jac(ind,:))';
            %Get Linear Constraint Indices
            ind = B.conlin == 1;
            %Check we have some constraints
            if(B.noCons && sum(ind) > 0)
                if(isempty(B.objbias)), B.objbias = zeros(size(ind)); end %HACK
                %Get Linear Jacobian
                prob.A = sparse(double(B.jac(ind,:)));
                prob.rl = B.cl(ind) - B.objbias(ind);
                prob.ru = B.cu(ind) - B.objbias(ind);
            end
            %Bounds & Integer Vars
            prob.lb = B.lb;
            prob.ub = B.ub;
            prob.int = B.xtype;
            prob.objbias = B.objbias(B.indobj);
            %Problem Name
            if(any(B.xtype=='I'))
                prob.Name = 'SymBuilder MILP';
            elseif(any(B.xtype=='B'))
                prob.Name = 'SymBuilder BILP';
            else
                prob.Name = 'SymBuilder LP';
            end
        end

        function prob = GetQuadProb(B)
        %Get Quadratic Problem Matrices [min 0.5*x'*H'x + f'*x
        %s.t. rl <= A*x <= ru, qrl <= x'*Q*x + l'*x <= qru]

            %Check is built
            IsBuilt(B,1);
            prob = optiprob;
            %Get quadratic objective indices
            ind = B.objlin == 2;
            if(isempty(ind))
                error('This object does not contain an objective');
            elseif(sum(ind) > 1)
                error('This interface only supports a single quadratic objective');
            end
            %Get Hessian (DON'T NEED 0.5!)
            if(~any(ind)) %linear objective
                ind = B.objlin == 1;
                n = size(B.jac,2);
                prob.H = spalloc(n,n,0);
                prob.f = double(B.jac(ind,:));
            else
                prob.H = sparse(double(B.hess(ind).H));
                %Get Gradient, remove 2nd derivative parts
                f = B.jac(ind,:); v = symvar(f);
                prob.f = double(subs(f,v,{zeros(size(v))}).');
            end

            %Get Linear Constraint Indices
            ind = B.conlin == 1;
            %Check we have some constraints
            if(B.noCons && sum(ind) > 0)
                %Get Linear Jacobian
                prob.A = sparse(double(B.jac(ind,:)));
                prob.rl = B.cl(ind) - B.objbias(ind);
                prob.ru = B.cu(ind) - B.objbias(ind);
            end
            %Get Quadratic Constraint Indices
            ind = B.conlin == 2;
            %Check we have some constraints
            if(B.noCons && sum(ind) > 0)
                prob.qrl = B.cl(ind) - B.objbias(ind);
                prob.qru = B.cu(ind) - B.objbias(ind);
                ind = find(ind);
                %Get each quadratic constraint
                len = length(ind);
                prob.Q = cell(len,1);
                prob.l = cell(len,1);
                for i = 1:len
                    %Get Hessian (0.5* requried!)
                    prob.Q{i} = 0.5*sparse(double(B.hess(ind(i)).H));
                    %Get Gradient, remove 2nd derivative parts
                    f = B.jac(ind(i),:); v = symvar(f);
                    prob.l{i} = double(subs(f,v,{zeros(size(v))}).');
                end
            end
            %Bounds & Integer Vars
            prob.lb = B.lb;
            prob.ub = B.ub;
            prob.int = B.xtype;
            prob.objbias = B.objbias(B.indobj);
            %Problem Name
            if(any(B.xtype=='I'))
                prob.Name = 'SymBuilder MIQP';
            else
                prob.Name = 'SymBuilder QP';
            end
        end

        function prob = GetNLProb(B,opts)
        %Get Nonlinear Problem Callbacks, generating files as we go

            if(nargin < 2), opts = symbset; end
            opts.verbose = B.verbose;
            %Check is built
            IsBuilt(B,1);
            prob = optiprob;
            %Determine most efficient callback form IF auto selected
            if(strcmpi(opts.cbmode,'auto'))
                %Only relevant if we have any nonlinear terms & > 8 variables
                if((any(B.objlin > 2) || any(B.conlin > 2)) && length(B.vars) > 8)
                    %See if we have 2nd deriv info
                    if(~isempty(B.hess))
                        %Generate Hess Lag (will need it later anyway)
                        buildHessLag(B);
                        %Build Hessian Structure
                        lHstr = zeros(size(B.hesslag));
                        lHstr(logical(B.hesslag ~= 0)) = 1;
                        %If greater than 10% dense, use CppAD
                        if(nnz(lHstr)/numel(lHstr) > 0.1)
                            if(SymBuilder.CheckCppADCompile(false))
                                opts.cbmode = 'cppad';
                            else
                                opts.cbmode = 'mcode';
                            end
                        else
                            opts.cbmode = 'mcode';
                        end
                    %Only 1st deriv info
                    else
                        ind = B.conlin > 0;
                        ajac = B.jac(ind,:);
                        %Build Jacobian Structure
                        jacstr = zeros(size(ajac));
                        jacstr(logical(ajac ~= 0)) = 1;
                        %If greater than 10% dense, use CppAD
                        if(nnz(jacstr)/numel(jacstr) > 0.1)
                            if(SymBuilder.CheckCppADCompile(false))
                                opts.cbmode = 'cppad';
                            else
                                opts.cbmode = 'mcode';
                            end
                        else
                            opts.cbmode = 'mcode';
                        end
                    end
                else
                    opts.cbmode = 'mcode';
                end
            end
            %Determine whether we are generating C or M files
            if(strcmpi(opts.cbmode,'mcode'))
                fgen = 'M';
            elseif(strcmpi(opts.cbmode,'cppad'))
                fgen = 'A';
            else
                fgen = 'C';
            end
            %Create New Callback Function based on req'd callbacks
            opts.havCon = any(B.conlin > 0);
            if(fgen == 'M')
                B.buildMFun('new',[],[],opts);
            else
                B.buildCFun('new',[],[],opts);
            end
            %Get Objective (Nonlinear or not, returned as f(x))
            ind = B.objlin >= 1;
            if(sum(ind) == 0)
                error('This object does not contain an objective');
            end
            %Build Objective Callback
            if(fgen == 'M')
                prob.fun = B.buildMFun('obj',B.sobj(ind),B.vars,opts);
            else
                prob.fun = B.buildCFun('obj',B.sobj(ind),B.vars,opts);
            end
            %If requested, build gradient callback
            if(strcmpi(opts.use1stDerivs,'yes'))
                if(fgen == 'A')
                    prob.f = @(x) symb_ccb('grad',x); %no need to generate
                else
                    ajac = B.jac(ind,:);
                    if(fgen == 'M')
                        prob.f = B.buildMFun('grad',ajac,B.vars,opts);
                    else
                        prob.f = B.buildCFun('grad',ajac,B.vars,opts);
                    end
                end
            elseif(fgen=='C')
                B.buildCFun('grad'); %default empty
            end

            %If we have constraints, create callbacks
            if(opts.havCon)
                %Get ALL Constraint Indices
                ind = B.conlin > 0;
                %Build Con Callback + row constraints
                if(fgen == 'M')
                    prob.nlcon = B.buildMFun('con',B.sobj(ind),B.vars,opts);
                else
                    prob.nlcon = B.buildCFun('con',B.sobj(ind),B.vars,opts);
                end
                prob.cl = B.cl(ind);
                prob.cu = B.cu(ind);
                %If we want derivatives, create jacobian + structure
                if(strcmpi(opts.use1stDerivs,'yes'))
                    if(fgen == 'A')
                        prob.nljac = @(x) symb_ccb('jac',x); %no need to generate
                        prob.nljacstr = @() symb_ccb('jacstr');
                    else
                        ajac = B.jac(ind,:);
                        if(fgen == 'M')
                            prob.nljac = B.buildMFun('jac',ajac,B.vars,opts);
                        else
                            prob.nljac = B.buildCFun('jac',ajac,B.vars,opts);
                        end
                        %Build Jacobian Structure Callback
                        jacstr = zeros(size(ajac));
                        jacstr(logical(ajac ~= 0)) = 1;
                        prob.nljacstr = @() sparse(jacstr);
                    end
                elseif(fgen=='C')
                    B.buildCFun('jac'); %default empty
                end
            elseif(fgen=='C' || fgen=='A')
                B.buildCFun('con',[],[],opts); %default empty
                B.buildCFun('jac',[],[],opts); %default empty
            end
            %If we want 2nd derivatives, add them too
            if(strcmpi(opts.use2ndDerivs,'yes'))
                if(fgen == 'A')
                    prob.H = @(x,sigma,lambda) symb_ccb('hess',x,sigma,lambda); %no need to generate
                    prob.Hstr = @() symb_ccb('hstr');
                else
                    %Check if we have done a full build
                    if(isempty(B.hess))
                        if(B.verbose), optiwarn('No 2nd Derivatives Added - You may only obtain the Hessian of the Lagrangian if you called Build().\nDraft() only computes 1st Derivatives%s\n','.'); end
                        if(fgen=='C')
                            B.buildCFun('hess'); %default empty
                        end
                    else
                        %Build HessLag
                        buildHessLag(B);
                        %Build NL Hessian Callback
                        if(fgen=='M')
                            prob.H = B.buildMFun('hess',B.hesslag,B.vars,opts);
                        else
                            prob.H = B.buildCFun('hess',B.hesslag,B.vars,opts,B.noCons);
                        end
                        %Build Hessian Structure
                        lHstr = zeros(size(B.hesslag));
                        lHstr(logical(B.hesslag ~= 0)) = 1;
                        prob.Hstr = @() sparse(lHstr);
                    end
                end
            elseif(fgen=='C')
                B.buildCFun('hess'); %default empty
            end

            %Bounds & Integer Vars
            prob.lb = B.lb;
            prob.ub = B.ub;
            prob.int = B.xtype;
            %Problem Name
            if(any(B.xtype=='I'))
                prob.Name = 'SymBuilder MINLP';
            else
                prob.Name = 'SymBuilder NLP';
            end

            %Generate pretend x0 incase none given
            prob.x0 = (prob.ub - prob.lb)./2;
            prob.x0(isnan(prob.x0)) = 0;
            prob.x0(isinf(prob.x0)) = 0;

            %Compile if C Code
            if(fgen ~= 'M')
                if(B.verbose), fprintf('Compiling....'); end
                %If compiling with VS2015 but pre R2015a, need to manually add in UCRT location
                cc = mex.getCompilerConfigurations(); post = [];
                for i = 1:length(cc)
                    if(~isempty(strfind(cc(i).Name,'Microsoft Visual C++')) && str2double(cc(i).Version) >= 14 && verLessThan('matlab','8.5'))
                        post = opti_FindUCRT();
                        break;
                    end
                end
                if(fgen == 'A')
                    src = 'symb_ccb.cpp';
                    %find cppad source (expected in utilities folder)
                    str = which('asl');
                    if(isempty(str)), error('Cannot find asl to locate CppAD source!'); end
                    idx = strfind(str,filesep); str = str(1:idx(end));
                    str = [str 'Source'];
                    str = ['mex -largeArrayDims symb_ccb.cpp -I"' str '" ' post];
                    if(strcmpi(computer,'PCWIN'))
                        str = [str ' -DCPPAD_SIZE_T_SAME_UNSIGNED_INT=1'];
                    end
                	eval(str);
                else
                    src = 'symb_ccb.c';
                    str = ['mex -largeArrayDims symb_ccb.c ' post];
                    eval(str);
                end
                if(B.verbose), fprintf('Done\n'); end
                if(strcmpi(opts.srckeep,'no'))
                    delete(src);
                end
            end
        end

        %Get Build Status
        function tf = IsBuilt(B,doErr)
            if(nargin < 2), doErr = false; end
            if(strcmp(B.bldstat,'unbuilt'))
                tf = false;
                if(doErr)
                    error('Please build the object using Draft() or Build()');
                end
            else
                tf = true;
            end
        end

        function [x,fval,ef,info] = solve(B,x0,opts)
        %Lowercase version
            if(nargin < 3), opts = []; end
            if(nargin < 2), x0 = []; end
            [x,fval,ef,info] = Solve(B,x0,opts);
        end

        function [x,fval,ef,info] = Solve(B,x0,opts)
        %Solve Object using OPTI

            %Check is built
            if(~IsBuilt(B))
                error('Please build the object first using Draft() or Build()');
            end
            if(nargin < 3 || isempty(opts)), opts = symbset; end
            if(nargin < 2), x0 = []; end

            %Cusomtize settings based on solver
            switch(lower(opts.solver))
                %White box solvers
                case {'scip','baron'}
                     opts = symbset(opts,'use1stDerivs','no','use2ndDerivs','no','preallocate','no');
                %Derivative free solvers
                case {'nomad','pswarm','gmatlab'}
                    opts = symbset(opts,'use1stDerivs','no','use2ndDerivs','no');
                %No Hessian Support
                case {'filtersd','lbfgsb','nlopt'}
                    opts = symbset(opts,'use2ndDerivs','no');
            end
            if (isempty(B.Opt) || isempty(B.lastSolveOpts) || ~isequal(B.lastSolveOpts, opts))
                if(B.verbose), fprintf('\nGenerating OPTI Object....\n'); end
                B.Opt = GetOPTI(B,opts);
                B.lastSolveOpts = opts;
                if(B.verbose), fprintf('Done\n\n'); end
            else
                if(B.verbose), fprintf('\nRe-using OPTI Object\n'); end
            end
            %Solve
            [x,fval,ef,info] = solve(B.Opt,x0);
        end

        function B = plus(B,C)
        %Add objective, constraint, bound or integer declaration to SymBuilder Object
            if(~isa(B,'SymBuilder'))
                error('The LHS object must be a SymBuilder Object');
            end
            if(~ischar(C))
                error('The RHS object must be a string');
            end
            %Check for objective
            idx = strfind(C,'min');
            if(~isempty(idx))
                if(length(idx) > 1), error('min most only appear once in the objective'); end
                B.AddObj(C(idx(1)+4:end));
                return;
            end
            idx = strfind(C,'max');
            if(~isempty(idx))
                if(length(idx) > 1), error('max most only appear once in the objective'); end
                B.AddObj(['-(' C(idx(1)+4:end) ')']);
                return;
            end
            %Check for bound or integer (no mathematical ops)
            idx = regexp(C,{'*','+','-','/','^'});
            if(~isempty([idx{:}]))
                B.AddCon(C);
            else
                idx = regexp(C,{'<=','>='});
                if(~isempty([idx{:}]))
                    B.AddBound(C);
                else
                    B.AddInteger(C);
                end
            end

        end
    end

    methods(Access=private)

        function buildSymRep(B)
        %Build symbolic representation of the equations

            %Check if something new has been added
            if(strcmp(B.bldstat,'unbuilt'))
                %Build Symbolic Object
                symobj = B.eqs;
                %Substitute expressions
                if(~isempty(B.exprsn))
                    %Now subs into full equation system
                    wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
                    symobj = subs(symobj,B.exprsn(:,1),B.exprsn(:,2));
                    warning(wstate);
                    %Have to repeat until all nested expressions are sub'd
                    n = 10;  %max depth
                    no = size(B.exprsn,1);
                    se = sym(B.exprsn(:,1));
                    while n > 0
                        v = symvar(symobj); alldone = 1;
                        for i = 1:no
                            if(any(se(i) == v))
                                wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
                                symobj = subs(symobj,B.exprsn(i,1),B.exprsn(i,2));
                                warning(wstate);
                                alldone = 0;
                            else
                                if(i == no && alldone) %ensure we have checked them all
                                    n = -10;
                                    break;
                                end
                            end
                        end
                        n = n - 1;
                    end
                    if(n == 0)
                        error('Maximum expression recursion depth reached!');
                    else
%                         fprintf('%d exp recursions required\n',n+10);
                    end
                end
                %Convert constants to decimal form
%                 digits(16); %set vpa digits
%                 BV = cell(size(B.constnt,1),1);
%                 for i = 1:size(B.constnt,1)
%                     BV{i} = sym(B.constnt{i,2},'d');
%                 end
                %Substitute constants
                if(~isempty(B.constnt))
                    wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
                    symobj = subs(symobj,B.constnt(:,1),B.constnt(:,2));
                    warning(wstate);
                end
                %Save symbolic vector
                B.sobj = symobj;
                %Save Variables
                B.vars = symvar(B.sobj);

                %Build default bounds and xtype
                n = length(B.vars);
                B.lb = -1*Inf(n,1);
                B.ub = Inf(n,1);
                B.xtype = repmat('C',1,n);
                %Get names of variables and their indices
                names = SymBuilder.detVarNames(B.vars);
                if(isempty(names) && size(B.bnds,1) > 0)
                    error('Cannot process bounds without declared variables!');
                end

                %Process bound declarations
                for i = 1:size(B.bnds,1)
                    %Check if bounds are numerical
                    llb = str2double(B.bnds{i,2});
                    if(isnan(llb)) %not a number, check constants
                        try
                            ind = strcmp(B.constnt(:,1),B.bnds{i,2});
                            if(any(ind))
                                llb = B.constnt{ind,2};
                            else
                                error('Unknown lower bound: %s',B.bnds{i,2});
                            end
                        catch
                            error('Error processing lower bound for var %d',i);
                        end
                    end
                    lub = str2double(B.bnds{i,3});
                    if(isnan(lub)) %not a number, check constants
                        try
                            ind = strcmp(B.constnt(:,1),B.bnds{i,3});
                            if(any(ind))
                                lub = B.constnt{ind,2};
                            else
                                error('Unknown upper bound: %s',B.bnds{i,3});
                            end
                        catch
                            error('Error processing upper bound for var %d',i);
                        end
                    end
                    %Check if variable exists
                    ind = logical(B.vars == sym(B.bnds{i,1}));
                    if(any(ind))
                        if(~isinf(llb)), B.lb(ind) = llb; end
                        if(~isinf(lub)), B.ub(ind) = lub; end
                    else %could be applied to multiple variables
                        ind = strcmp(names(:,1),B.bnds{i,1});
                        if(isempty(ind) || all(ind  == 0))
                            error('Unknown bound: %s',B.bnds{i,1});
                        elseif(sum(ind) > 1)
                            error('Bound name is ambiguous and matches more than one variable name');
                        else
                            %Extract indices
                            start = names{ind,2}(1);
                            send = names{ind,2}(end);
                            if(~isinf(llb)), B.lb(start:send) = llb; end
                            if(~isinf(lub)), B.ub(start:send) = lub; end
                        end
                    end
                end

                %Process integer declarations
                for i = 1:size(B.vartypes,1)
                    %Check if variable exists
                    ind = logical(B.vars == sym(B.vartypes{i,1}));
                    if(any(ind))
                         B.xtype(ind) = B.vartypes{i,2};
                    else %could be applied to multiple variables
                        ind = strcmp(names(:,1),B.vartypes{i,1});
                        if(isempty(ind))
                            error('Unknown variable: %s',B.vartypes{i,1});
                        elseif(sum(ind) > 1)
                            error('Variable name is ambiguous and matches more than one variable name');
                        else
                            %Extract indices
                            start = names{ind,2}(1);
                            send = names{ind,2}(end);
                            B.xtype(start:send) = B.vartypes{i,2};
                        end
                    end
                end
            end
        end

        function buildJac(B)
        %Evaluate Jacobian of entire equation system

            if(isempty(B.sobj))
                error('You must build a symbolic representation of the system first!');
            end
            %Call SymToolbox jacobian function
            B.jac = jacobian(B.sobj);
        end

        function buildHess(B)
        %Evaluate Hessian(s) of entire equation system

            if(isempty(B.jac))
                error('You must evaluate the Jacobian of the system first!');
            end
            noeq = B.noObjs + B.noCons;
            B.hess = struct('H',cell(noeq,1),'lin',cell(noeq,1));
            %Build Index Vectors
            B.conlin = zeros(noeq,1);
            B.objlin = zeros(noeq,1);
            B.objbias = zeros(noeq,1);
            %Go through each equation
            for i = 1:noeq
                J = B.jac(i,:);
                %Determine if linear
                if(isempty(symvar(J)))
                    B.hess(i).lin = 1;
                    B.hess(i).H = 0;
                else
                    %Hessian is the Jacobian of the Jacobian (in this case)
                    H = jacobian(J,B.vars);
                    s = symvar(H);
                    %Determine if quadratic
                    if(isempty(s))
                        B.hess(i).lin = 2;
                        B.hess(i).H = H;
                    else
                        %Otherwise must be nonlinear, save entire symobj (for now)
                        B.hess(i).lin = 3;
                        B.hess(i).H = H;
                    end
                end
                %If linear or quadratic, save bias term
                if(B.hess(i).lin <= 2)
                    eq = B.sobj(i);
                    v = symvar(eq);
                    %Find bias term
                    B.objbias(i) = double(subs(eq,v,{zeros(size(v))}));
                end

                %Save linearity
                if(any(B.indobj == i))
                    B.objlin(i) = B.hess(i).lin;
                else
                    B.conlin(i) = B.hess(i).lin;
                end
            end
        end

        function buildHessLag(B)
        %Build Lagrangian Hessian assuming IPOPT (tril) format

            if(isempty(B.hess))
                error('You must evaluate the Hessian of the system first!');
            elseif(~isempty(B.hesslag))
                return; %already done
            end
            %Start with objective
            nvars = length(B.vars);
            ind = B.objlin >= 2;
            if(any(ind)) %check we have a nonlinear (or quadratic) objective
                if(sum(ind) > 1)
                    error('This interface only support single objective problems');
                end
                s = sym('sigma');
                objH = B.hess(ind).H;
                H = s*objH;
            else
                H = zeros(nvars);
            end
            %Now each nonlinear (or quadratic) constraint
            ind = B.conlin > 1;
            %Two sets of indexes, one for lambda, one for B.hess
            index_BH = find(ind);
            ind(B.indobj) = [];
            index_L = find(ind);
            wstate = warning('off','symbolic:sym:sym:DeprecateExpressions');
            for i = 1:length(index_BH)
                l = sym(sprintf('lambda(%d)',index_L(i)));
                H = H + l*B.hess(index_BH(i)).H;
            end
            warning(wstate);
            %Save resulting Hessian
            B.hesslag = tril(H);
        end

        function detLinearity(B)
        %Determine linearity of the supplied equations (without Hessian)

            if(isempty(B.jac))
                error('You must evaluate the Jacobian of the system first!');
            end
            %Build Index Vectors
            noeq = B.noCons + B.noObjs;
            B.conlin = zeros(noeq,1);
            B.objlin = zeros(noeq,1);
            B.objbias = zeros(noeq,1);
            %Determine Linearity (based on existence of symvars)
            for i = 1:noeq %must search all equations
                %Gather first derivative vars
                s = symvar(B.jac(i,:));
                %Objective
                if(any(B.indobj == i))
                    if(isempty(s))
                        B.objlin(i) = 1; %lin
                    else
                        B.objlin(i) = 3; %nonlin
                    end
                else
                    if(isempty(s))
                        B.conlin(i) = 1;
                    else
                        B.conlin(i) = 3;
                    end
                end
            end
            % Find objective bias terms
            for i = 1:noeq
                % If linear, save bias term
                if(B.objlin(i) == 1 || B.conlin(i) == 1)
                    eq = B.sobj(i);
                    v = symvar(eq);
                    %Find bias term
                    B.objbias(i) = double(subs(eq,v,{zeros(size(v))}));
                end
            end
        end
    end

    methods(Static)
        %Parse Constraint String
        [str,l,u] = parseConstraint(str);
        %Extract RHS of expression string, remove from eq, and return LHS
        [str,exp] = parseExpression(str);
        %Extract Group and Name from Result Expression Name
        [group,name] = parseResExp(name);
        %Parse bound string
        [var,lb,ub] = parseBound(str);
        %Parse integer string
        [var,xtype] = parseInteger(str);
        %Determine Variable names
        varn = detVarNames(svar);
        %Build MATLAB function file
        cb = buildMFun(mode,sobj,svar,opts);
        %Build C Code function file
        cb = buildCFun(mode,sobj,svar,opts,nocon);
        %Convert symbolic expression into matlab function
        fun = sym2fun(sobj,svar,var,skipSubs);

        %Check we have a compiler suitable for use with Cppad
        function ok = CheckCppADCompile(doerr)
            ok = true;
            c = mex.getCompilerConfigurations('c++');
            if(isempty(c) || (isempty(strfind(c.Name,'Microsoft Visual C++')) && isempty(strfind(c.Name,'Intel C++'))))
                str = sprintf(['CppAD compilation is only available with Microsoft Visual Studio or Intel as the C++ compiler.\n\nPlease either select it'...
                           ' as the MEX compiler (via mex -setup), or download from Microsoft:\n\n %s (2012, recommended)\n\nor\n\n %s (2013)'],...
                           'http://www.microsoft.com/en-nz/download/details.aspx?id=34673','http://www.visualstudio.com/en-us/products/visual-studio-express-vs.aspx');
                if(doerr)
                    error(str) %#ok<SPERR>
                else %warning, auto has selected it
                    ok = false;
                    optiwarn('symb:cppad',sprintf('SymBuilder has chosen CppAD has the most effective callback strategy, but %s',str));
                end
            end
        end
    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [V,D,permutation,failure] = dmpermblockeig(X,switchtosparse)

[permutation,aux1,aux2,blocks] = dmperm(X+speye(length(X)));
Xpermuted = X(permutation,permutation);

V = [];
D = [];
V = zeros(size(X,1),1);
top = 1;
left = 1;
anycholfail = 0;
failure = 0;

for i = 1:length(blocks)-1
    Xi = Xpermuted(blocks(i):blocks(i+1)-1,blocks(i):blocks(i+1)-1);
    [R,fail] = chol(Xi);
    anycholfail = anycholfail | fail;
    if fail && nnz(Xi)>0
        if length(Xi) >= switchtosparse
            [vi,di,eigfail] = eigs(Xi,5,'SA');
            if eigfail || isempty(di)
                res = 0;
                for j = 1:size(vi,2)
                    res(j) = norm(Xi*vi(:,j)-vi(:,j)*di(j,j));
                end
                % We only trust these
                notfailed = abs(res) <= 1e-12;
                vi = vi(:,notfailed);
                di = di(notfailed,notfailed);
                if length(vi) == 0
                    [vi,di,eigfail] = eigs(sparse(Xi),25,'SA');
                    if eigfail
                        res = 0;
                        for j = 1:size(vi,2)
                            res(j) = norm(Xi*vi(:,j)-vi(:,j)*di(j,j));
                        end
                        % We only trust these
                        notfailed = abs(res) <= 1e-12;
                        vi = vi(:,notfailed);
                        di = di(notfailed,notfailed);
                    end
                end
            end
        else
            [vi,di] = eig(full(Xi));
        end
        for j = 1:length(di)
            if di(j,j)<=0
                V(top:top+length(Xi)-1,left)=vi(:,j);
                left = left + 1;
                D = blkdiag(D,di(j,j));
            end
        end
    end
    top = top + length(Xi);
end

if (anycholfail && isempty(V)) || (anycholfail && all(diag(D)>0))
    % OK, we have a problem. The Cholesky factorization failed for some of
    % the matrices, but yet no eigenvalue decomposition revealed a negative
    % eigenvalue (due to convergence issues in the sparse eigs)
    failure = 1;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function aa=lpcdl2aa(dl)
%LPCDL2AA dct of log area to area coefficients AA=(DL)

% we do not correct for sinc errors

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcdl2aa.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p]=size(dl);
aa=[zeros(nf,1) exp(irdct(dl.').') ones(nf,1)];

author: decouples


text:
function m = max(T)

[e,v] = eig(T);
m = max(diag(v));

author: mtex-toolbox

repo: mtex

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


text:
clear, clf
fc= 9e8;          % 2e9;   % 2GHz
fs=5e4;            % 5[MHz]
speed_kmh=120;%3;     % 3[km/h]
Ts=1/fs;           % sampling frequency in sec
v_ms= speed_kmh/3.6;  % velocity[m/s]
wl_m= 3e8/fc;      % wavelength[m]
% channel parameters setting: SCM case 2
PDP_dB=[0. -1. -9. -10. -15. -20];
t_ns=[0 310 710 1090 1730 2510];
BS_theta_LOS_deg=0;
MS_theta_LOS_deg=0;
BS_AS_deg=2;                        % Laplacian PAS
BS_AoD_deg=50*ones(size(PDP_dB));
MS_AS_deg=35;                   % for Lapalcian PAS
DoT_deg=22.5;
MS_AoA_deg=67.5*ones(size(PDP_dB));
% generates phase of a subray
[BS_theta_deg,MS_theta_deg,BS_PHI_rad]=gen_phase(BS_theta_LOS_deg,BS_AS_deg,BS_AoD_deg,MS_theta_LOS_deg,MS_AS_deg,MS_AoA_deg);
PDP=dB2w(PDP_dB);
% generates coefficients
%for k=1:10000
%   t=(k-1)*Ts;
%   h(k,:)=ray_fading0(20,PDP,BS_PHI_rad,MS_theta_deg,v_ms,DoT_deg,wl_m,t);
%end
%plot([1:10000]*Ts,10*log10(abs(h(:,1))))
t=[0:9999]*Ts;
h= ray_fading(20,PDP,BS_PHI_rad,MS_theta_deg,v_ms,DoT_deg,wl_m,t);
plot(t,10*log10(abs(h(1,:))))
title(['Ray Channel Model, f_c=',num2str(fc),'Hz, T_s=',num2str(Ts),'s']);
xlabel('time[s]'), ylabel('Magnitude[dB]')

author: LyricYang


text:
% SquareRoot is the dagnn wrapper which performs square root normalization
% for the input features at each location

% Copyright (C) 2015 Tsung-Yu Lin, Aruni RoyChowdhury, Subhransu Maji.
% All rights reserved.
%
% This file is part of the BCNN. It uses MATCONVNET package and is made
% available under the terms of the BSD license (see the COPYING file).

classdef SquareRoot < dagnn.ElementWise
  properties
    param = 1e-8;
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnsqrt(inputs{1}, obj.param) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, param, derOutputs)
      derInputs{1} = vl_nnsqrt(inputs{1}, obj.param, derOutputs{1}) ;
      derParams = {} ;
    end

    function rfs = getReceptiveFields(obj)
      rfs.size = [1 1] ;
      rfs.stride = [1 1] ;
      rfs.offset = [1 1] ;
    end

    function obj = SquareRoot(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: DBCNN


text:
function img = discreteSquare(varargin)
%DISCRETESQUARE Discretize a planar square
%
%   IMG = discreteSquare(DIM, CENTER, SIDE)
%   DIM is the size of image, with the format [x0 dx x1;y0 dy y1]
%   CENTER is the center of the square
%   SIDE is the length of the side of the square.
%
%   IMG = discreteSquare(DIM, CENTER, SIDE, THETA)
%   Also specify spherical angle of the normal of a face of the square.
%   THETA is the angle with the horizontal, in degrees, counted counter-
%   clockwise in direct basis (and clockwise in image basis).
%
%   IMG = discreteSquare(DIM, SQUARE)
%   send parameters in a row vector, where SQUARE contains at least the
%   center coordinate, and possibly the other parameters.
%
%   IMG = discreteSquare(LX, LY, ...);
%   Specifes the pixels coordinates with the two row vectors LX and LY.
%
%   Example
%   img = discreteSquare([1 1 100;1 1 100], [50 50], 30);
%   img = discreteSquare([1 1 100;1 1 100], [50 50 30]);
%   img = discreteSquare([1 1 100;1 1 100], [50 50 30 30]);
%
%   See Also
%   imShapes, discreteDisc, discreteRectangle
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-05-16
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   12/10/2006: typo in the doc
%   04/01/2007: concatenate transforms before applying them
%   19/06/2007: update doc
%   04/03/2009: use meshgrid
%   29/04/2009: update transforms
%   29/05/2009: use more possibilities for specifying grid
%   22/01/2010: fix auto center with odd image size
%   2011-03-30 use degrees

% compute coordinate of image voxels
[lx, ly, varargin] = parseGridArgs(varargin{:});
[x, y]   = meshgrid(lx, ly);

% default parameters
center = [lx(ceil(end/2)) ly(ceil(end/2))];
side    = center;
theta   = 0;

% process input parameters
if length(varargin)==1
    % all paramater in the first argument
    var = varargin{1};
    center = var(:,1:2);
    if size(var, 2)>2
        side = var(:,3);
    end
    if size(var, 2)>3
        theta = var(:,4);
    end

elseif ~isempty(varargin)
    % first argument is center, optionnally followed by side length, then
    % by rotation angle
    center = varargin{1};
    if length(varargin)>1
        side = varargin{2};
    end
    if length(varargin)>2
        theta = varargin{3};
    end
end

% For squares, size of sides is the same in all directions
side    = [side(1) side(1)];

% transforms voxels according to square orientation
tra     = createTranslation(-center);
rot     = createRotation(-deg2rad(theta));
sca     = createScaling(1./side);
[x, y]  = transformPoint(x, y, sca*rot*tra);

% create image: simple threshold over 2 dimensions
img = abs(x) <= .5 & abs(y) <= .5;

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [geneKO] = calculateGeneKOMatrix(model, varargin)
% Build a rxn-gene matrix such that the i-th column indicates what
% reactions become inactive because of the i-th gene's knock-out.
%
% USAGE:
%
%    geneKO = calculateGeneKOMatrix(model, SeparateTranscript, printLevel)
%
% INPUT:
%    model:             The COBRA Model structure
%
% OPTIONAL INPUTS:
%    SeparateTranscript Character used to separate
%                       different transcripts of a gene. (default = '')
%                       Examples:
%                           - SeparateTranscript = ''
%                              - gene 10005.1    ==>    gene 10005.1
%                              - gene 10005.2    ==>    gene 10005.2
%                              - gene 10005.3    ==>    gene 10005.3
%                           - SeparateTranscript = '.'
%                              - gene 10005.1
%                              - gene 10005.2    ==>    gene 10005
%                              - gene 10005.3
%    printLevel:       Integer. 1 if the process is wanted to be shown
%                      on the screen, 0 otherwise. (default = 1)
%
% OUTPUT:
%    geneKO:            Struct which contains matrix with blocked
%                       reactions for each gene in the metabolic model,
%                       name of reactions and name of genes.
%
% .. Authors:
%       - Luis V. Valcarcel, Oct 2017, University of Navarra, CIMA & TECNUN School of Engineering.
%       - Luis V. Valcarcel, 03/11/2018, University of Navarra, CIMA & TECNUN School of Engineering.

p = inputParser; % check the input information
addRequired(p, 'model');
addOptional(p, 'SeparateTranscript', '', @ischar);
addOptional(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x));
parse(p, model, varargin{:});

% define gene set
genes = unique(strtok(model.genes, p.Results.SeparateTranscript));
ngenes = numel(genes);

% generate output matrix
ko_rxn_gene = zeros(numel(model.rxns),ngenes);

% affected reactions
if p.Results.printLevel > 0
    showprogress(0, 'Calculate Gene Knock-out matrix');
end
for gen = 1:ngenes
    if p.Results.printLevel > 0
        showprogress(gen/ngenes, 'Calculate Gene Knock-out matrix');
    end

    transcripts = model.genes(strcmp(strtok(model.genes,p.Results.SeparateTranscript),genes{gen})); % to support R2015b
    [~, hasEffect, constrRxnNames] = deleteModelGenes(model, transcripts);

    if hasEffect
        % search index of bloced reactions
        [~,idx] = ismember(constrRxnNames,model.rxns);
        ko_rxn_gene(idx,gen) = 1;
    end
end
if p.Results.printLevel > 0
    fprintf('\tGeneKOMatrix calculated\n');
end

% geneKO
geneKO.genes = genes;
geneKO.rxns = model.rxns;
geneKO.matrix = (ko_rxn_gene~=0);
end

author: opencobra

repo: cobratoolbox


text:
function [lags,use_ascend,E] = perform_extraction(points_ordered,options)

% PERFORM_EXTRACTION   Run lag extraction
%
%   Run lag extraction on a reordered raster plot (2D Image)
%
%   SYNTAX
%       [LAGS,USE_ASCEND,E] = PERFORM_EXTRACTION(POINTS_ORDERED,OPTIONS)
%

% $Id: perform_extraction.m 4 2009-08-15 21:10:35Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-08-15 17:10:35 -0400 (Sam, 15 aoû 2009) $
% $Revision: 4 $

if nargin<3
    options.null = 0;
end

if ~isfield(options, 'xwin')
    options.xwin = 1;
end
xwin = options.xwin;

if nargin == 0
    rmfield(options,'null')
    options = options
    return
end

%%%%%%%% END OPTIONS %%%%%%%%%

if xwin>0
    data_smooth = movav(points_ordered,[],2*xwin,1); % smoothing for better robustness
else
    data_smooth = points_ordered;
end

if options.disp_log
    disp('---- Running lag extraction');
end

[lags,use_ascend,dc,sc,E] = gc_lags(data_smooth,options);

if options.disp_log
    disp('-- Lag extraction done !');
end

if 0
    figure
    imagesc(data_smooth)
    colorbar
    xlabel('Time (ms)')
    ylabel('Trial')
    title('Points Ordered and Smoothed')
end


text:
% eeg_point2lat() - convert latency in data points to latency in ms relative
%                   to the time locking. Used in eeglab().
% Usage:
%       >> [newlat ] = eeg_point2lat( lat_array, [], srate);
%       >> [newlat ] = eeg_point2lat( lat_array, epoch_array,...
%                                 srate, timelimits, timeunit);
% Inputs:
%   lat_array   - latency array in data points assuming concatenated
%                 data epochs (see eeglab() event structure)
%   epoch_array - epoch number corresponding to each latency value
%   srate       - data sampling rate in Hz
%   timelimits  - [min max] timelimits in 'timeunit' units (see below)
%   timeunit    - time unit in second. Default is 1 = seconds.
%
% Outputs:
%   newlat      - converted latency values (in 'timeunit' units) for each epoch
%
% Example:
%   tmpevent = EEG.event;
%   eeg_point2lat( [ tmpevent.latency ], [], EEG.srate, [EEG.xmin EEG.xmax]);
%   % returns the latency of all events in second for a continuous
%   % dataset EEG
%
%   eeg_point2lat( [ tmpevent.latency ], [ tmpevent.epoch ],
%                 EEG.srate, [EEG.xmin EEG.xmax]*1000, 1E-3);
%   % returns the latency of all events in millisecond for a dataset
%   % containing data epochs.
%
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2 Mai 2002
%
% See also: eeg_lat2point(), eeglab(), pop_editieventvals(), pop_loaddat()

% Copyright (C) 2 Mai 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% 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 COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function newlat = eeg_point2lat( lat_array, epoch_array, srate, timewin, timeunit);

if nargin <3
    help eeg_point2lat;
    return;
end
if isempty( epoch_array )
    epoch_array = ones( size(lat_array) );
end
if nargin <4
    timewin = 0;
end
if nargin <5
	timeunit = 1;
end

if length(lat_array) ~= length(epoch_array)
	if length(epoch_array)~= 1
		disp('eeg_point2lat: latency and epoch arrays must have the same length'); return;
	else
		epoch_array = ones(1,length(lat_array))*epoch_array;
	end
end
if length(timewin) ~= 2
    disp('eeg_point2lat: timelimits array must have length 2'); return;
end
if iscell(epoch_array)
	epoch_array = [ epoch_array{:} ];
end
if iscell(lat_array)
	lat_array = [ lat_array{:} ];
end

timewin = timewin*timeunit;

if length(timewin) == 2
    pnts = (timewin(2)-timewin(1))*srate+1;
else
    pnts = 0;
end
newlat  = ((lat_array - (epoch_array-1)*pnts-1)/srate+timewin(1))/timeunit;
newlat = round(newlat*1E9)*1E-9;


text:
function h = r8mat_house_post ( n, a, row, col )

%*****************************************************************************80
%
%% R8MAT_HOUSE_POST computes a Householder post-multiplier matrix.
%
%  Discussion:
%
%    H(ROW,COL) has the property that the ROW-th column of
%    A*H(ROW,COL) is zero from entry COL+1 to the end.
%
%    In the most common case, where a QR factorization is being computed,
%    ROW = COL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrices.
%
%    Input, real A(N,N), the matrix whose Householder matrix
%    is to be computed.
%
%    Input, integer ROW, COL, specify the location of the
%    entry of the matrix A which is to be preserved.  The entries in
%    the same row, but higher column, will be zeroed out if
%    A is postmultiplied by H.
%
%    Output, real H(N,N), the Householder matrix.
%

%
%  Set up the vector V.
%
  a_row(1,1:col-1) = 0.0;
  a_row(1,col:n) = a(row,col:n);

  a_row = a_row';

  v = r8vec_house_column ( n, a_row, col );
%
%  Form the matrix H(V).
%
  h = r8mat_house_form ( n, v );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = cn_jac_monomial_integral ( n, alpha, beta, expon )

%*****************************************************************************80
%
%% CN_JAC_MONOMIAL_INTEGRAL: integral of a monomial with Jacobi weight over CN.
%
%  Discussion:
%
%    value = integral ( CN )
%      product ( 1 <= i <= n ) x(I)^expon(i) (1-x(i))^alpha (1+x(i))^beta dx(i)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, the exponent of (1-X) in the weight factor.
%
%    Input, real BETA, the exponent of (1+X) in the weight factor.
%
%    Input, integer EXPON(N), the exponents.
%
%    Output, real VALUE, the value of the integral.
%
  value = 1.0;
  for i = 1 : n
    value2 = c1_jac_monomial_integral ( alpha, beta, expon(i) );
    value = value * value2;
  end

  return
end