text:
function docinfo
% This is a placeholder matlab function to provide info about creating
% documentation in gerardus.
%
% HTML files created with the matlab publisher go in the /doc folder
% For them to appear as text/html mimetype on a server such as googlecode, the
% appropriate svn mimetypes need to be set. An example is shown below in this
% function's file.
%
% also see:
% http://manjeetdahiya.com/2010/09/29/serving-html-documentation-from-google-code-svn/
% https://code.google.com/p/flexlib/wiki/HowToBuild
% http://uk.mathworks.com/help/matlab/matlab_prog/display-custom-documentation.html
%
% Click here for Documentation with examples:

%%% The following is example lines from a configuration file as found in
%%% ~/.subversion.config on posix systems, to set appropriate mimetypes for
%%% common file-extensions.

% ### Section for configuring miscelleneous Subversion options.
% [miscellany]
% ### Set enable-auto-props to 'yes' to enable automatic properties
% ### for 'svn add' and 'svn import', it defaults to 'no'.
% ### Automatic properties are defined in the section 'auto-props'.
% enable-auto-props = yes
%
% ### Section for configuring automatic properties.
% [auto-props]
% ### The format of the entries is:
% ###   file-name-pattern = propname[=value][;propname[=value]...]
% ### The file-name-pattern can contain wildcards (such as '*' and
% ### '?').  All entries which match (case-insensitively) will be
% ### applied to the file.  Note that auto-props functionality
% ### must be enabled, which is typically done by setting the
% ### 'enable-auto-props' option.
%
% ### Mime-types for common file types (adapted from: https://code.google.com/p/flexlib/wiki/HowToBuild)
% ### Also see: http://manjeetdahiya.com/2010/09/29/serving-html-documentation-from-google-code-svn/
%
% # Scriptish formats
% *.bat        = svn:eol-style=native; svn:keywords="Rev Date Author"; svn-mine-type=text/plain
% *.bsh        = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/x-beanshell
% *.cgi        = svn:eol-style=native; svn:keywords="Rev Date Author"; svn-mine-type=text/plain
% *.cmd        = svn:eol-style=native; svn:keywords="Rev Date Author"; svn-mine-type=text/plain
% *.js         = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/javascript
% *.php        = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/x-php
% *.pl         = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/x-perl; svn:executable
% *.pm         = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/x-perl
% *.py         = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/x-python; svn:executable
% *.sh         = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/x-sh; svn:executable
%
% # Image formats
% *.bmp        = svn:mime-type=image/bmp
% *.gif        = svn:mime-type=image/gif
% *.ico        = svn:mime-type=image/ico
% *.jpeg       = svn:mime-type=image/jpeg
% *.jpg        = svn:mime-type=image/jpeg
% *.png        = svn:mime-type=image/png
% *.tif        = svn:mime-type=image/tiff
% *.tiff       = svn:mime-type=image/tiff
%
% # Data formats
% *.pdf        = svn:mime-type=application/pdf
% *.avi        = svn:mime-type=video/avi
% *.eps        = svn:mime-type=application/postscript
% *.gz         = svn:mime-type=application/gzip
% *.mov        = svn:mime-type=video/quicktime
% *.mp3        = svn:mime-type=audio/mpeg
% *.ps         = svn:mime-type=application/postscript
% *.psd        = svn:mime-type=application/photoshop
% *.rtf        = svn:mime-type=text/rtf
% *.swf        = svn:mime-type=application/x-shockwave-flash
% *.tgz        = svn:mime-type=application/gzip
% *.wav        = svn:mime-type=audio/wav
% *.zip        = svn:mime-type=application/zip
%
% # Microsoft-specific formats (see: http://filext.com/faq/office_mime_types.php)
% *.doc        = svn:mime-type=application/msword
% *.docx       = svn:mime-type=application/vnd.openxmlformats-officedocument.wordprocessingml.document
% *.xls        = svn:mime-type=application/vnd.ms-excel
% *.xlsx       = svn:mime-type=application/vnd.openxmlformats-officedocument.spreadsheetml.sheet
% *.ppt        = svn:mime-type=application/vnd.ms-powerpoint
% *.pptx       = svn:mime-type=application/vnd.openxmlformats-officedocument.presentationml.presentation
%
% # Matlab-specific binary formats (adapted from: http://www.mathworks.co.uk/help/matlab/matlab_prog/set-up-svn-source-control.html)
% *.mdl        = svn:mime-type=application/octet-stream
% *.mat        = svn:mime-type=application/octet-stream
% *.slx        = svn:mime-type= application/octet-stream
% *.mdlp       = svn:mime-type=application/octet-stream
% *.slxp       = svn:mime-type=application/octet-stream
% *.sldd       = svn:mime-type=application/octet-stream
% *.p          = svn:mime-type=application/octet-stream
% *.mexa64     = svn:mime-type=application/octet-stream
% *.mexw32     = svn:mime-type=application/octet-stream
% *.mexw64     = svn:mime-type=application/octet-stream
% *.mexmaci64  = svn:mime-type=application/octet-stream
%
% # Text formats
% .htaccess    = svn:mime-type=text/plain
% *.css        = svn:mime-type=text/css
% *.dtd        = svn:mime-type=text/xml
% *.html       = svn:mime-type=text/html
% *.ini        = svn:mime-type=text/plain
% *.sql        = svn:mime-type=text/x-sql
% *.txt        = svn:mime-type=text/plain
% *.xhtml      = svn:mime-type=text/xhtml+xml
% *.xml        = svn:mime-type=text/xml
% *.xsd        = svn:mime-type=text/xml
% *.xsl        = svn:mime-type=text/xml
% *.xslt       = svn:mime-type=text/xml
% *.xul        = svn:mime-type=text/xul
% *.yml        = svn:mime-type=text/plain
% CHANGES      = svn:mime-type=text/plain
% COPYING      = svn:mime-type=text/plain
% INSTALL      = svn:mime-type=text/plain
% Makefile*    = svn:mime-type=text/plain
% README       = svn:mime-type=text/plain
% TODO         = svn:mime-type=text/plain
%
% # Code formats
% *.c          = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain
% *.cpp        = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain
% *.h          = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain
% *.java       = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain
% *.as         = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain
% *.mxml       = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain
% *.m          = svn:eol-style=native; svn:keywords="Rev Date Author"; svn:mime-type=text/plain


text:
function [boxedX,Bx,By]=boxpdf(X)
% Forces the pdf of data to have a boxed distribution using a data adaptive lookup table.
%
% [boxedX,Bx,By]=boxpdf(X)
%
% boxedX=N(X) where N is an data adaptive monotonically increasing function.
% boxedX vary between zero and one.
%
% Bx,By describes the lookup table
%
% (c) Aslak Grinsted 2002-2004

%tested with: [prctile(x,boxpdf(x)*100)-x]

%home made fast unique:
[Bx,Js]=sort(X(:));
n=size(Bx,1);
d=(diff(Bx)~=0);

I=(1:n)';
I=I([d;logical(1)]);
%if (nargout>1)
    Bx=Bx(I);
%end

J=cumsum([1;d]);

n=length(X);
I=[0;I];

By=(.5/n)*(I(1:(end-1))+(I(2:end))); %percentile

%boxedX=By(J(Js([2:end 1])));
boxedX=interp1q(Bx,By,X);
%
% %old boxpdf:
%
% [sx i]=sort(X(:));
% [n, m] = size(sx);
% minx  = min(sx);
% maxx  = max(sx);
% range = maxx-minx;
%
%
% % Use the same Y vector if all columns have the same count
% eprob = [0.5./n:1./n:(n - 0.5)./n]';
%
%
% %group x's with same value
% I=[1; find(diff(sx)~=0)+1]; %unique indices
% Bx=sx(I); %unique values
% nI=diff([I;length(sx)+1])-1;
% By=zeros(length(Bx),1);
% for i=length(I):-1:1
%
%     By(i,1)=mean(eprob(I(i)+(0:nI(i))));
% end
% %plot(sx,y);
%
%
%
% boxedX=interp1q(Bx,By,X);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function plotGraph(vertices, edges, figureRef)

if nargin < 3
    figureRef = figure;
end

cla(figureRef);

figure(figureRef);
plot(vertices(:,1),vertices(:,2),'.');
hold on;

nv = size(vertices,1);
ne = size(edges,1);

v1 = vertices(edges(:,1),:);
v2 = vertices(edges(:,2),:);

vList = ones(ne*3,2);

vList(1:3:ne*3,:) = v1;
vList(2:3:ne*3,:) = v2;
vList(3:3:ne*3,:) = NaN;

plot(vList(:,1),vList(:,2));

% for i=1:size(edges,1)
% %     plot([vertices(edges(i,1),1) vertices(edges(i,2),1)], ...
% %          [vertices(edges(i,1),2) vertices(edges(i,2),2)]);
%
%     line([vertices(edges(i,1),1) vertices(edges(i,2),1)], ...
%          [vertices(edges(i,1),2) vertices(edges(i,2),2)]);
% end

axis equal;

author: Doraemonzzz


text:
%% OPTI Toolbox Mixed-Integer Linear Programming Examples
%
% This file contains a number of MILP problems and demonstrates how to
% solve them using the OPTI Toolbox. You should read and complete
% BasicUsage.m & LinearProgramming.m BEFORE running the below examples.
%
%   Copyright (C) 2014 Jonathan Currie (IPL)

% There is also a page on the Wiki which supplements this example:
web('https://www.inverseproblem.co.nz/OPTI/index.php/Probs/MILP');

%% Determing which Solver to Use
% OPTI Toolbox comes with a number of MILP solvers, thus to determine which
% ones are available on your system you can type:
clc
optiSolver('MILP')

%% Example 1
% This is a simple two decision variable MILP which will use for the next
% few examples
clc
f = -[6 5]';                %Objective Function (min f'x)
A = [1,4; 6,4; 2, -5];      %Linear Inequality Constraints (Ax <= b)
b = [16;28;6];
lb = [0;0];                 %Bounds on x (lb <= x <= ub)
ub = [10;10];
xtype = 'II';               %Integer Variables (I = integer, C = continuous, B = binary)

% Building an MILP problem is very similar to an LP, except just add the
% 'xtype' argument for integer variables
Opt = opti('f',f,'ineq',A,b,'bounds',lb,ub,'xtype',xtype)

%Call solve to solve the problem
[x,fval,exitflag,info] = solve(Opt)

%% Example 2 - Alternative Integer Setup
% You can also supply an vector of integer indices indicating the position of
% continuous and integer variables respectively (note no binary variables
% can be entered this way)

xtype = [1 2];

Opt = opti('f',f,'ineq',A,b,'bounds',lb,ub,'xtype',xtype);
solve(Opt);

%% Example 3 - Plotting the Solution
% Several problem types have a default plot command available. Note for
% MILP plots it will also plot the integer constraints.

plot(Opt)

%% Example 4 - Solving a slightly bigger MILP
%Build the opti problem:
clc

%Problem
f = -[1 2 3 1]';
A = [-1 1 1 10; 1 -3 1 0];
b = [20;30];
Aeq = [0 1 0 -3.5];
beq = 0;
lb = [0;0;0;2];
ub = [40;inf;inf;3];
xtype = 'CCCI';

Opt = opti('f',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'xtype',xtype)

[x,fval,exitflag,info] = solve(Opt)

%% Example 5 - Calling a solver directly
% All OPTI Linear and Quadratic solvers use standard MATLAB function
% prototypes (such as linprog or quadprog), so you can skip using the OPTI
% class all together. Note be sure to check the integer argument form for
% the particular solver by typing e.g. help opti_glpk
clc
[x,fval,exitflag,info] = opti_glpk(f,A,b,Aeq,beq,lb,ub,xtype)

%% Example 6 - Sparse MILPs
% As with LPs, all solvers are setup to directly solve sparse systems, which
% is the preferred format for most solvers:
clc
% A larger sparse MILP
load sparseMILP1;

opts = optiset('solver','glpk');    %Solve with GLPK
Opt = opti('f',f,'ineq',A,b,'eq',Aeq,beq,'xtype',find(xint),'options',opts)
[x,fval,exitflag,info] = solve(Opt);
fval
info

%% Problem 4
% MILP with Special Ordered Sets (SOS)
clc
%Problem
f = [-1 -1 -3 -2 -2]';
A = [-1 -1 1 1 0;
      1 0 1 -3 0];
b = [30;30];
lb = zeros(5,1);
ub = [40;1;inf;inf;1];

%SOS type 1
sos_type = '1';
sos_index = [1 2 3 4 5]';
sos_weight = [1 2 3 4 5]';

% Build the problem, specifying the three SOS fields. Note only some MILP
% solvers are setup to solve problems with SOS:
opts = optiset('solver','cbc');
Opt = opti('f',f,'ineq',A,b,'bounds',lb,ub,'sos',sos_type,sos_index,sos_weight,'options',opts)

[x,fval,exitflag,info] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef ElasticityMicroDefinedByLevelSet < CompositeResultsPrinter

    methods (Access = public)

        function obj = ElasticityMicroDefinedByLevelSet(d)
            obj.simulationStr = 'ElasticityMicroDefinedByLevelSet';
            obj.init(d);
        end

    end

    methods (Access = protected)

        function createPrinters(obj,d)
            obj.printers{1} = ResultsPrinter.create('DensityGauss',d);
            obj.printers{2} = ResultsPrinter.create('LevelSet',d);
            obj.printers{3} = ResultsPrinter.create('ElasticityMicro',d);
        end

        function storeFieldsToPrint(obj,d)
            obj.storeDensityGaussFields(d);
            obj.storeLevelSetFields(d);
            obj.storeMicroFields(d);
        end

    end

    methods (Access = protected)

        function storeDensityGaussFields(obj,d)
            di.fields = d.dens;
            obj.printers{1}.storeFieldsToPrint(di);
        end

        function storeLevelSetFields(obj,d)
            di.x = d.levelSet;
            obj.printers{2}.storeFieldsToPrint(di);
        end

        function storeMicroFields(obj,d)
            d.fields = d.microVar;
            obj.printers{3}.storeFieldsToPrint(d);
            obj.printers{3}.setStrVariablesMicroCase([])
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function SO3F = conv(SO3F1,SO3F2,varargin)
% convolution of an SO3FunHarmonic with a function or a kernel on SO(3)
%
% 1) SO3Fun * SO3Fun
% There are two SO3Funs $f: _{S_f^L\backslash}SO(3)_{/S_f^R} \to \mathbb{C}$
% where $S_f^L$ is the Left symmetry and $S_f^R$ is the Right symmetry and
% $g: _{S_g^L\backslash}SO(3)_{/S_g^R} \to \mathbb{C}$ given.
% Then the convolution $ f *_L g : _{S_f^L\backslash}SO(3)_{/S_g^R} \to
% \mathbb{C}$ is defined by
%
% $$ (f *_L g)(R) = \frac1{8\pi^2} \int_{SO(3)} f(q) \cdot g(q^{-1}\,R) \, dq $$
%
% and the convolution $ f *_R g : _{S_g^L\backslash}SO(3)_{/S_f^R} \to
% \mathbb{C}$ is defined by
%
% $$ (f *_R g)(R) = \frac1{8\pi^2} \int_{SO(3)} f(q) \cdot g(R\,q^{-1}) \, dq $$.
%
% with $vol(SO(3)) = \int_{SO(3)} 1 \, dR = 8\pi^2$.
% The convolution $*_L$ is used as default.
% The convolution of matrices of SO3Functions with matrices of SO3Functions
% works elementwise.
%
% 2) SO3Fun * S2Fun
% The convolution of an SO3Fun  $f: _{S_f^L\backslash}SO(3)_{/S_f^R} \to \mathbb{C}$
% with an S2Fun $h: \mathbb S^2_{/S_h} \to \mathbb{C}$ yields
% $f*h:\mathbb S^2_{/S_f^L} \to \mathbb{C}$ with
%
% $$ (f * h)(\xi) =  \frac1{8\pi^2} \int_{SO(3)} f(q) \cdot h(q^{-1}\,\xi) \, dq $$.
%
% 3) In particular we convolute an SO3Fun with an SO3Kernel similar to the
% first case. Therefore the Right and Left sided convolution are
% equivalent. The convolution of an SO3Fun with an S2Kernel works analogue
% to case 2.
%
%
% Syntax
%   SO3F = conv(SO3F1,SO3F2)
%   SO3F = conv(SO3F1,SO3F2,'Right')
%   SO3F = conv(SO3F1,psi)
%   sF2 = conv(SO3F1,sF1)
%   sF2 = conv(SO3F1,phi)
%
% Input
%  SO3F1, SO3F2 - @SO3Fun
%  psi          - convolution @SO3Kernel
%  sF1          - @S2Fun
%  phi          - convolution @S2Kernel
%
% Output
%  SO3F - @SO3FunHarmonic
%  sF2  - @S2FunHarmonic
%
% See also
% SO3Kernel/conv SO3FunHarmonic/conv SO3FunRBF/calcFourier S2FunHarmonic/conv S2Kernel/conv

% The convolution is defined like above. But in MTEX the convolution of two
% SO3Funs is mostly calculated by
%                    conv(inv(conj(SO3F1)),SO3F2).
%

if isnumeric(SO3F1)
  SO3F = conv(SO3F2,SO3F1,varargin{:});
  return
end
if isnumeric(SO3F2)
  v = 2*(1:length(SO3F2))'-1;
  SO3F2 = SO3Kernel(SO3F2.*v);
end

% ------------------- convolution with a S2Kernel -------------------
if isa(SO3F2,'S2Kernel')
  psi = SO3F2;
  L = min(SO3F1.bandwidth,psi.bandwidth);

  fhat = zeros((L+1)^2,1);
  for l = 0:L
    fhat(l^2+1:(l+1)^2) = 2*sqrt(pi)/(2*l+1) * ...
          SO3F1.fhat(deg2dim(l)+(2*l+1)*l+(1:2*l+1)')*psi.A(l+1);
  end

  warning(['There is no symmetry given for the S2Kernel function. But for convolution the ' ...
      'right symmetry of the SO3Fun has to be compatible with the symmetry of the S2Fun.'])
  SO3F = S2FunHarmonic(fhat);
  return

end

% ------------------- convolution with a S2Fun -------------------
if isa(SO3F2,'S2Fun')

  sF = S2FunHarmonic(SO3F2);
  L = min(SO3F1.bandwidth,sF.bandwidth);

  fhat = zeros((L+1)^2,1);
  for l = 0:L
    fhat(l^2+1:(l+1)^2) = reshape(SO3F1.fhat(deg2dim(l)+1:deg2dim(l+1)),2*l+1,2*l+1) * ...
          sF.fhat(l^2+1:(l+1)^2) ./ sqrt(2*l+1);
  end

  % we need that SO3F1.SRight == SO3F.Sym
  if numProper(SO3F1.SLeft) == 1
    SO3F = S2FunHarmonicSym(fhat,SO3F1.SRight);
  elseif isa(SO3F2,'S2FunHarmonicSym')
    ensureCompatibleSymmetries(SO3F1,SO3F2);
    SO3F = S2FunHarmonicSym(fhat,SO3F1.SRight);
  else
    warning(['There is no symmetry of the S2Fun given. But for convolution the ' ...
      'right symmetry of the SO3Fun has to be compatible with the unknown symmetry.'])
    SO3F = S2FunHarmonic(fhat);
  end

  return

end

% ------------------- convolution with a SO3Kernel -------------------
% ( Here *L is the same like *R )
if isa(SO3F2,'SO3Kernel')

  L = min(SO3F1.bandwidth,SO3F2.bandwidth);
  SO3F1.bandwidth = L;
  s = size(SO3F1); SO3F1 = SO3F1.subSet(':');

  % multiply Wigner-D coefficients of SO3F1
  % with the Chebyshev coefficients A of SO3F2
  A = SO3F2.A;
  for l = 0:L
    SO3F1.fhat(deg2dim(l)+1:deg2dim(l+1),:) = ...
      A(l+1)./(2*l+1) * SO3F1.fhat(deg2dim(l)+1:deg2dim(l+1),:);
  end

  SO3F = reshape(SO3F1,s);
  return

end

% ------------------- convolution of SO3Fun's -------------------
% i) right sided convolution
if check_option(varargin,'Right')
  SO3F = inv(conv(inv(SO3F1),inv(SO3F2)));
  return
end

% ii) left sided convolution (default)
ensureCompatibleSymmetries(SO3F1,SO3F2,'conv_Left');

L = min(SO3F1.bandwidth,SO3F2.bandwidth);

% compute Fourier coefficients
fhat1 = SO3F1.calcFourier('bandwidth',L);
fhat2 = SO3F2.calcFourier('bandwidth',L);

% construct SO3FunHarmonic by multiplying the Fourier coefficients
SO3F = SO3FunHarmonic(convSO3(fhat1,fhat2),SO3F2.SRight,SO3F1.SLeft);

end

author: mtex-toolbox

repo: mtex

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


text:
function output = callfiltersd(model)

options = model.options.filtersd;

model = yalmip2nonlinearsolver(model);

% Make sure the callback returns a dense Jacobian
model.dense = 1;

if ~model.derivative_available
    disp('Derivate-free call to filterSD not yet implemented')
    error('Derivate-free call to filterSD not yet implemented')
end
if model.options.savedebug
    save filtersddebug model
end
showprogress('Calling filterSD',model.options.showprogress);

cu = [ repmat(0,length(model.bnonlinineq),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(0,length(model.b),1);
    repmat(0,length(model.beq),1)];

cl = [ repmat(-inf,length(model.bnonlinineq),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(-inf,length(model.b),1);
    repmat(0,length(model.beq),1)];

if isempty(cu)
    Flow = [];
    Fupp = [];
end

% These are needed to avoid recomputation due to ipopts double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G = [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

funcs.objective = @(x)ipopt_callback_f(x,model);
funcs.gradient = @(x)ipopt_callback_df(x,model);
if ~isempty(cu)
    funcs.constraints = @(x)ipopt_callback_g(x,model);
    funcs.jacobian  = @(x)ipopt_callback_dg(x,model);
else
    funcs.constraints = [];
    funcs.jacobian = [];
end

lb  = model.lb(:)';
ub = model.ub(:)';

if ~model.options.usex0
    model.x0 = (lb+ub)/2;
    model.x0(isinf(ub)) = lb(isinf(ub))+1;
    model.x0(isinf(lb)) = ub(isinf(lb))-1;
    model.x0(isinf(model.x0)) = 0;
end

options.display = model.options.verbose;
solvertime = tic;
[xout,fval,exitflag,stats,lambda] = filtersd(funcs.objective, funcs.gradient, model.x0, lb, ub, funcs.constraints, funcs.jacobian, cl, cu, options);
solvertime = toc(solvertime);

% Duals currently not supported
lambda = [];

x = RecoverNonlinearSolverSolution(model,xout);

switch exitflag
    case {0}
        problem = 0;
    case {1}
        problem = 2;
    case {2,3}
        problem = 1;
    case {5}
        problem = 3;
    case {4,9}
        problem = 4;
    case 105
        problem = 16;
    otherwise
        problem = -1;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.exitflag = exitflag;
    solveroutput.stats = stats;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(x,D_struc,[],problem,'filterSD',solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [P,f,alpha] = fastlomb(x,t,varargin)
% FASTLOMB caculates the Lomb normalized periodogram (aka Lomb-Scargle, Gauss-Vanicek or Least-Squares spectrum) of a vector x with coordinates in t.
%   The coordinates need not be equally spaced. In fact if they are, it is
%   probably preferable to use PWELCH or SPECTRUM. For more details on the
%   Lomb normalized periodogram, see the excellent section 13.8 in [1], pp.
%   569-577.
%
%   This code is a transcription of the Fortran subroutine fasper in [1]
%   (pp.575-577), so it is a really fast (albeit not really exact)
%   implementation of the Lomb periodogram. Also Matlab's characteristics
%   have been taken into account in order to make it even faster for
%   Matlab. For an exact calculation of the Lomb periodogram use LOMB,
%   which is however about 100 times slower.
%
%   If the 'fig' flag is on, the periodogram is created, along with some
%   default statistical significance levels: .001, .005, .01, .05, .1, .5.
%   If the user wants more significance levels, they can give them as input
%   to the function. Those will be red.
%
% SYNTAX
%   [P,f,alpha] = fastlomb(x,t,fig,hifac,ofac,a);
%   [P,f,alpha] = fastlomb(x,t,fig,hifac,ofac);
%   [P,f,alpha] = fastlomb(x,t,fig,hifac);
%   [P,f,alpha] = fastlomb(x,t,fig);
%   [P,f,alpha] = fastlomb(x,t);
%
% INPUTS
%   x:     the vector whose spectrum is wanted.
%   t:     coordinates of x (should have the same length).
%   fig:   if 0 (default), no figure is created.
%   hifac: the maximum frequency returned is
%               (hifac) x (average Nyquist frequency)
%          See "THE hifac PARAMETER" in the NOTES section below for more
%          details on the hifac parameter. Default is 1 (i.e. max frequency
%          is the Nyquist frequency).
%   ofac:  oversampling factor. Typically it should be 4 or larger. Default
%          is 4. See "INTERPRETATION AND SELECTION OF THE ofac PARAMETER"
%          in the NOTES section below for more details on the choice of
%          ofac parameter.
%   a:     additional significance levels to be drawn on the figure.
%
% OUTPUTS
%   P:     the Lomb normalized periodogram
%   f:     respective frequencies
%   alpha: statistical significance for each value of P
%
% NOTES
% A. INTERPRETATION AND SELECTION OF THE ofac PARAMETER [1]
%    The lowest independent frequency f to be examined is the inverse of
%    the span of the input data,
%               1/(tmax-tmin)=1/T.
%    This is the frequency such that the data can include one complete
%    cycle. In an FFT method, higher independent frequencies would be
%    integer multiples of 1/T . Because we are interested in the
%    statistical significance of ANY peak that may occur, however, we
%    should over-sample more finely than at interval 1/T, so that sample
%    points lie close to the top of ANY peak. This oversampling parameter
%    is the ofac. A value ofac >~4 might be typical in use.
%
% B. THE hifac PARAMETER [1]
%    Let fhi be the highest frequency of interest. One way to choose fhi is
%    to compare it with the Nyquist frequency, fc, which we would obtain, if
%    the N data points were evenly spaced over the same span T, that is
%               fc = N/(2T).
%    The input parameter hifac, is defined as fhi/fc. In other words, hifac
%    shows how higher (or lower) that the fc we want to go.
%
% REFERENCES
%   [1] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
%       "Numerical recipes in Fortran 77: the art of scientific computing",
%       2nd ed., vol. 1, Cambridge University Press, NY, USA, 2001.
%
% C. Saragiotis, Nov 2008

%% Inputs check and initialization
if nargin < 2, error('%s: there must be at least 2 inputs.',mfilename); end
MACC = 10;   % Number of interpolation points per 1/4 cycle of the highest frequency

[x,t,hifac,ofac,a_usr,f,fig] = init(x,t,varargin{:});
nt = length(t);

mx = mean(x);
x  = x-mx;
vx = var(x);
if vx==0, error('x has zero variance'); end

nf = length(f);
nfreq  = 2^nextpow2(ofac*hifac*nt*MACC);
% ndim = 2*nfreq;

tmin = t(1);
tmax = t(end);
T = tmax-tmin;

%% Extirpolation
fac = 2*nfreq/(T*ofac);
[wk1,wk2] = extirpolate(t-tmin,x,fac,nfreq,MACC);

%% Take the FFT's
W = fft(wk1);
reft1 = real(W(2:nf+1)); % only positive frequencies, without f=0
imft1 = imag(W(2:nf+1));

W = fft(wk2);
reft2 = real(W(2:nf+1));
imft2 = imag(W(2:nf+1));

%% Main
hypo  = sqrt(reft2.^2+imft2.^2);
hc2wt = 0.5*reft2./hypo;
hs2wt = 0.5*imft2./hypo;
cwt   = sqrt(0.5+hc2wt);
swt   = (sign(hs2wt)+(hs2wt==0)).*(sqrt(0.5-hc2wt));
den   = 0.5*nt + hc2wt.*reft2 + hs2wt.*imft2;
cterm = (cwt.*reft1 + swt.*imft1).^2./den;
sterm = (cwt.*imft1 - swt.*reft1).^2./(nt-den);
P  = (cterm+sterm)/(2*vx);
P  = P(:);

%% Significance
M = 2*nf/ofac;
alpha = 1 - (1-exp(-P)).^M;              % statistical significance
alpha(alpha<0.1) = M*exp(-P(alpha<0.1)); % (to avoid round-off errors)

%% Figure
if fig
    figure
    styles = {':','-.','--'};

    a = [0.001 0.005 0.01 0.05 0.1 0.5];
    La = length(a);
    z = -log(1-(1-a).^(1/M));
%     widths = [.5*ones(1,floor(La/3)) 1*ones(1,La-2*floor(La/3)) 2*ones(1,floor(La/3))];
    hold on;
    for i=1:La
        line([f(1),0.87*f(end)],[z(i),z(i)],'Color','k','LineStyle',styles{ceil(i*3/La)});%,'LineWidth',widths(i));
        text(0.9*f(end),z(i),strcat('\alpha = ',num2str(a(i))),'fontsize',8); %,'fontname','symbol');
%         lgd{i}=strcat('\alpha=',num2str(a(i)));
    end
    if ~isempty(a_usr)
        [tmp,ind] = intersect(a_usr,a);
        a_usr(ind)=[];
        La_usr = length(a_usr);
        z_usr  = -log(1-(1-a_usr).^(1/M));
        for i = 1:La_usr
            line([f(1),0.87*f(end)],[z_usr(i),z_usr(i)],'Color','r','LineStyle',styles{ceil(i*3/La_usr)});%,'LineWidth',widths(i));
            text(0.9*f(end),z_usr(i),strcat('\alpha = ',num2str(a_usr(i))),'fontsize',8); %,'fontname','symbol');
    %         lgd{La+i}=strcat('\alpha=',num2str(a_usr(i)));
        end
    z = [z z_usr];
    end
%     legend(lgd);
    plot(f,P,'k');
    title('Lomb-Scargle normalized periodogram')
    xlabel('f (Hz)'); ylabel('P(f)')
    xlim([0 f(end)]); ylim([0,1.2*max([z'; P])]);
    hold off;
end

end

%% #### Local functions

%% init (initialize)
function [x,t,hifac,ofac,a,f,fig] = init(x,t,varargin)
    if nargin < 6, a = [];    % set default value for a
    else           a = sort(varargin{4});
                   a = a(:)';
    end
    if nargin < 5, ofac = 4;  % set default value for ofac
    else           ofac = varargin{3};
    end
    if nargin < 4, hifac = 1; % set default value for hifac
    else           hifac = varargin{2};
    end
    if nargin < 3, fig = 0;   % set default value for hifac
    else           fig = varargin{1};
    end

    if isempty(ofac),  ofac  = 4; end
    if isempty(hifac), hifac = 1; end
    if isempty(fig),   fig   = 0; end

    if ~isvector(x) ||~isvector(t),
        error('%s: inputs x and t must be vectors',mfilename);
    else
        x = x(:); t = t(:);
        nt = length(t);
        if length(x)~=nt
            error('%s: Inputs x and t must have the same length',mfilename);
        end
    end

    [t,ind] = unique(t);    % discard double entries and sort t
    x = x(ind);
    if length(x)~=nt, disp(sprintf('WARNING %s: Double entries have been eliminated',mfilename)); end

    T = t(end) - t(1);
    nf = round(0.5*ofac*hifac*nt);
    f = (1:nf)'/(T*ofac);
end

%% extirpolation
function [wk1,wk2] = extirpolate(t,x,fac,nfreq,MACC)
    nt = length(x);
    wk1 = zeros(2*nfreq,1);
    wk2 = zeros(2*nfreq,1);

    for j = 1:nt
        ck  = 1 + mod(fix(t(j)*fac),2*nfreq);
        ckk = 1 + mod(2*(ck-1),     2*nfreq);
        wk1 = spread(x(j),wk1,ck, MACC);
        wk2 = spread(1   ,wk2,ckk,MACC);
    end
end

%% spread
function yy = spread(y,yy,x,m)

    nfac = [1 1 2 6 24 120 720 5040 40320 362880];
    n = length(yy);

    if x == round(x)  % Original Fortran 77 code: ix=x
                      %                           if(x.eq.float(ix)) then...
        yy(x) = yy(x) + y;
    else
        % Original Fortran 77 code: i1 = min(  max( int(x-0.5*m+1.0),1 ), n-m+1  )
        i1 = min([  max([ floor(x-0.5*m+1),1 ]),  n-m+1  ]);
        i2 = i1+m-1;
        nden = nfac(m);
        fac = x-i1;

        fac = fac*prod(x - (i1+1:i2));
        yy(i2) = yy(i2) + y*fac/(nden*(x-i2));

        for j = i2-1:-1:i1
            nden = (nden/(j+1-i1))*(j-i2);
            yy(j) = yy(j) + y*fac/(nden*(x-j));
        end

    end
end

author: Sable

repo: mcbench-benchmarks


text:
filename='Gripping_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function u = nvecs(X,n,r,opts)
%NVECS Compute the leading mode-n vectors for a ktensor.
%
%   U = NVECS(X,n,r) computes the r leading eigenvalues of Xn*Xn'
%   (where Xn is the mode-n matricization of X), which provides
%   information about the mode-n fibers. In two-dimensions, the r
%   leading mode-1 vectors are the same as the r left singular vectors
%   and the r leading mode-2 vectors are the same as the r right
%   singular vectors.
%
%   U = NVECS(X,n,r,OPTS) specifies options:
%   OPTS.eigsopts: options passed to the EIGS routine [struct('disp',0)]
%   OPTS.flipsign: make each column's largest element positive [true]
%
%   See also KTENSOR, TENMAT, EIGS.
%
%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 ~exist('opts','var')
    opts = struct;
end

if isfield(opts,'eigsopts')
    eigsopts = opts.eigsopts;
else
    eigsopts.disp = 0;
end

% Compute Xn * Xn' excluding the nth factor
M = X.lambda * X.lambda';
for i = 1:ndims(X)
    if i == n, continue, end;
    M = M .* (X.u{i}' * X.u{i});
end

% Compute Xn * Xn'
Y = X.u{n} * M * X.u{n}';

[u,d] = eigs(Y, r, 'LM', eigsopts);

if isfield(opts,'flipsign')
    flipsign = opts.flipsign;
else
    flipsign = true;
end

if flipsign
    % Make the largest magnitude element be positive
    [val,loc] = max(abs(u));
    for i = 1:r
        if u(loc(i),i) < 0
            u(:,i) = u(:,i) * -1;
        end
    end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%%******************************************************************
%% mybicgstab
%%
%% [xx,resnrm,flag] = mybicgstab(A,b,M1,tol,maxit)
%%
%% iterate on  bb - (M1)*AA*x
%%
%% r = b-A*xtrue;
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************
  function [xx,resnrm,flag] = mybicgstab(A,b,M1,tol,maxit,printlevel)

  N = length(b);
  if (nargin < 6); printlevel = 1; end
  if (nargin < 5) | isempty(maxit); maxit = max(30,length(A.mat22)); end;
  if (nargin < 4) | isempty(tol); tol = 1e-10; end;
  tolb = min(1e-4,tol*norm(b));
  flag = 1;

  x = zeros(N,1);
  if (norm(x))
     if isstruct(A); r = b-matvec(A,x); else; r = b-mexMatvec(A,x); end;
  else
     r =b;
  end
  err = norm(r); resnrm(1) = err;  minresnrm = err; xx = x;
  %%if (err < 1e-3*tolb); return; end

  omega  = 1.0;
  r_tld = r;
%%
%%
%%
  breakyes = 0;
  smtol = 1e-40;
  for iter = 1:maxit,

     rho   = (r_tld'*r);
     if (abs(rho) < smtol)
        flag = 2;
        if (printlevel); fprintf('*'); end;
        breakyes = 1;
        break;
     end
     if (iter > 1)
        beta  = (rho/rho_1)* (alp/omega);
        p = r + beta*(p - omega*v);
     else
        p = r;
     end
     p_hat = precond(A,M1,p);
     if isstruct(A); v = matvec(A,p_hat); else; v = mexMatvec(A,p_hat); end;
     alp = rho / (r_tld'*v);
     s = r - alp*v;
     %%
     s_hat = precond(A,M1,s);
     if isstruct(A); t = matvec(A,s_hat); else; t = mexMatvec(A,s_hat); end;
     omega = (t'*s) / (t'*t);
     x = x + alp*p_hat + omega*s_hat;
     r = s - omega*t;
     rho_1 = rho;
     %%
     %% check convergence
     %%
     err = norm(r); resnrm(iter+1) = err;
     if (err < minresnrm);
        xx = x; minresnrm = err;
     end
     if (err < tolb)
        break;
     end
     if (err > 10*minresnrm)
        if (printlevel); fprintf('^'); end
        breakyes = 2;
        break;
     end
     if (abs(omega) < smtol)
        flag = 2;
        if (printlevel); fprintf('*'); end
        breakyes = 1;
        break;
     end
  end
  if (~breakyes) & (printlevel >=3); fprintf(' '); end
%%
%%*************************************************************************
%%*************************************************************************
%% precond:
%%*************************************************************************

   function Mx = precond(A,L,x)

   m = L.matdim; m2 = length(x)-m;
   Mx = zeros(length(x),1);

   for iter = 1
      if norm(Mx); r = x - matvec(A,Mx); else; r = x; end
      if (m2 > 0)
         r1 = full(r(1:m));
      else
         r1 = full(r);
      end
      if (m2 > 0)
         r2 = r(m+[1:m2]);
         w = linsysolvefun(L,r1);
         z = mexMatvec(A.mat12,w,1) - r2;
         z = L.Mu \ (L.Ml \ (L.Mp*z));
         r1 = r1 - mexMatvec(A.mat12,z);
      end
      d = linsysolvefun(L,r1);
      if (m2 > 0)
         d = [d; z];
      end
      Mx = Mx + d;
   end
%%*************************************************************************
%%*************************************************************************
%% matvec: matrix-vector multiply.
%% matrix = [A.mat11, A.mat12; A.mat12', A.mat22]
%%*************************************************************************

   function Ax = matvec(A,x);

   m = length(A.mat11); m2 = length(x)-m;
   if issparse(x); x = full(x); end
   if (m2 > 0)
      x1 = x(1:m);
   else
      x1 = x;
   end
   Ax = mexMatvec(A.mat11,x1);
   if (m2 > 0)
      x2 = x(m+[1:m2]);
      Ax = Ax + mexMatvec(A.mat12,x2);
      Ax2 = mexMatvec(A.mat12,x1,1) + mexMatvec(A.mat22,x2);
      Ax = [full(Ax); full(Ax2)];
   end
%%*************************************************************************

author: yu-jiang

repo: radpbook

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


text:
function imsegs2 = transferLabels(imsegs1, imsegs2, l2v, l2h)
%
% imsegs2 = transferLabels(imsegs1, imsegs2)
%
% Transfers labeling from imsegs1 to imsegs2.  imsegs1.segimage and ground
% truth labels are used to create a label map, which is then used to label
% the superpixels in imsegs2.  l2v and l2h give conversions from full label
% index to vertical and horz sublabel indices, respectively.
% (default: l2v = [1 2 2 2 2 2 3], l2h = [0 1 2 3 4 5 0])

if ~exist('l2v')
    l2v = [1 2 2 2 2 2 3];
end
if ~exist('l2h')
    l2h = [0 1 2 3 4 5 0];
end

for f = 1:numel(imsegs1)

    disp(num2str(f))

    stats = regionprops(imsegs1(f).segimage, 'PixelIdxList');
    idx = {stats.PixelIdxList};
    lmap = zeros(size(imsegs1(f).segimage));
    for s = 1:imsegs1(f).nseg
        lmap(idx{s}) = imsegs1(f).labels(s);
    end

    % set imsegs2
    nseg = imsegs2(f).nseg;
    imsegs2(f).vlabels = cell(nseg, 1);
    imsegs2(f).hlabels = cell(nseg, 1);
    imsegs2(f).labels = zeros(nseg, 1);
    imsegs2(f).vert_labels = zeros(nseg, 1);
    imsegs2(f).horz_labels = zeros(nseg, 1);
    imsegs2(f).label_names = imsegs1(f).label_names;
    imsegs2(f).vert_names = imsegs1(f).vert_names;
    imsegs2(f).horz_names = imsegs1(f).horz_names;

    stats = regionprops(imsegs2(f).segimage, 'PixelIdxList');
    idx = {stats.PixelIdxList};

    for s = 1:nseg
        mcl = mode(lmap(idx{s}));

        imsegs2(f).labels(s) = mcl;
        if mcl==0
            imsegs2(f).vert_labels(s) = 0;
            imsegs2(f).horz_labels(s) = 0;
        else
            imsegs2(f).vert_labels(s) = l2v(mcl);
            imsegs2(f).horz_labels(s) = l2h(mcl);
        end

        if imsegs2(f).vert_labels(s)==0
            imsegs2(f).vlabels{s} = '---';
        else
            imsegs2(f).vlabels{s} = imsegs2(f).vert_names{imsegs2(f).vert_labels(s)};
        end
        if imsegs2(f).horz_labels(s)==0
            imsegs2(f).hlabels{s} = '---';
        else
            imsegs2(f).hlabels{s} = imsegs2(f).horz_names{imsegs2(f).horz_labels(s)};
        end
    end
end

author: zouchuhang

repo: LayoutNet


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

temp = (X*Theta').*R;
J = sum( sum( (temp - Y.*R).^2) )/2.0 + (lambda/2) * ( sum(sum( X.^2 )) + sum(sum( Theta.^2 )) );
% J = sum( sum( (temp - Y.*R).^2) )/2.0 + (lambda/2) * ( sum(sum( X.^2 )) + sum(sum( Theta.^2 )) ) ;

X_grad = (temp - Y.*R) * Theta + lambda * X;
Theta_grad = (temp - Y.*R)' * X + lambda * Theta;

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

grad = [X_grad(:); Theta_grad(:)];

end

author: Ayatans

repo: Machine-Learning-homework


text:
function [ n_data, n, x, fx ] = h_polynomial_values ( n_data )

%*****************************************************************************80
%
%% H_POLYNOMIAL_VALUES: tabulated values of H(i,x).
%
%  Discussion:
%
%    H(i,x) is the physicist's Hermite polynomial of degree I.
%
%    In Mathematica, the function can be evaluated by:
%
%      HermiteH[n,x]
%
%  Differential equation:
%
%    Y'' - 2 X Y' + 2 N Y = 0
%
%  First terms:
%
%      1
%      2 X
%      4 X^2     -  2
%      8 X^3     - 12 X
%     16 X^4     - 48 X^2     + 12
%     32 X^5    - 160 X^3    + 120 X
%     64 X^6    - 480 X^4    + 720 X^2    - 120
%    128 X^7   - 1344 X^5   + 3360 X^3   - 1680 X
%    256 X^8   - 3584 X^6  + 13440 X^4  - 13440 X^2   + 1680
%    512 X^9   - 9216 X^7  + 48384 X^5  - 80640 X^3  + 30240 X
%   1024 X^10 - 23040 X^8 + 161280 X^6 - 403200 X^4 + 302400 X^2 - 30240
%
%  Recursion:
%
%    H(0,X) = 1,
%    H(1,X) = 2*X,
%    H(N,X) = 2*X * H(N-1,X) - 2*(N-1) * H(N-2,X)
%
%  Norm:
%
%    Integral ( -oo < X < +oo ) exp ( - X^2 ) * H(N,X)^2 dX
%    = sqrt ( PI ) * 2^N * N!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the polynomial.
%
%    Output, real X, the point where the polynomial is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 18;

  fx_vec = [ ...
      0.1000000000000000E+01, ...
      0.1000000000000000E+02, ...
      0.9800000000000000E+02, ...
      0.9400000000000000E+03, ...
      0.8812000000000000E+04, ...
      0.8060000000000000E+05, ...
      0.7178800000000000E+06, ...
      0.6211600000000000E+07, ...
      0.5206568000000000E+08, ...
      0.4212712000000000E+09, ...
      0.3275529760000000E+10, ...
      0.2432987360000000E+11, ...
      0.1712370812800000E+12, ...
      0.0000000000000000E+00, ...
      0.4100000000000000E+02, ...
     -0.8000000000000000E+01, ...
      0.3816000000000000E+04, ...
      0.3041200000000000E+07 ];

  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10, 11, ...
    12,  5,  5, ...
     5,  5,  5 ];

  x_vec = [ ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     5.0E+00, ...
     0.0E+00, ...
     0.5E+00, ...
     1.0E+00, ...
     3.0E+00, ...
     1.0E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [lineHandles textHandles] = hline(y, lineType, label, textPosition, axesHandle)

% Draws horizontal lines (specified by a vector of y)
%
% Except the first argument 'y', everything else are optional.
% To preserve input argument order, enter [] for unused input arguments
%
% lineType:     same as the ones used in plot(...)
% label:        simple strings applies to all lines
%               cell strings applies to each line
% textPosition: [x, y] bewteen 0 and 1 relative to the 'lower' end
%
% Author:      Hoi Wong (wonghoi.ee@gmail.com)
% Date:        02/14/2008
%
% Acknowledgement:  It was based on hline() written by Brandon Kuczenski.

    if( ~isvector(y) )
        error('y must be a vector');
    else
        y=y(:);     % Make sure it's column vector
    end

    if( ~exist('label', 'var') )        label = [];         end
    if( ~exist('lineType', 'var') )     lineType = [];      end
    if( ~exist('axesHandle', 'var') )   axesHandle = gca;   end
    if( ~exist('textPosition', 'var') ) textPosition = [];  end

    if( isempty(axesHandle) )           axesHandle = gca;   end

    if( isempty(textPosition) )
        textPositionX = 0.02;
        textPositionY = 0.02;
    elseif( isscalar(textPosition) )
        textPositionX = textPosition;
        textPositionY = 0.02;
    elseif( length(textPosition)~=2 )
        error('Invalid textPosition');
    else
        textPositionX = textPosition(1);
        textPositionY = textPosition(2);
    end
    clear textPosition;     % Avoid typos for similarly named variables

    holdState = ishold(axesHandle);
    hold(axesHandle, 'on');

    xLimits=get(axesHandle,'xlim');             % Row vector
    Xlimits=repmat(xLimits', 1, length(y));

    % Example: for horizontal lines
    % X = [2 2        Y = [3 4
    %      5 5];           3 4];

    if( isempty(lineType) )
        lineHandles = plot(axesHandle, Xlimits, [y';y']);
    else
        lineHandles = plot(axesHandle, Xlimits, [y';y'], lineType);
    end

    if( ~isempty(label) )
        yLimits = get(axesHandle,'ylim');
        yLowerLimit = yLimits(2);
        yUpperLimit = yLimits(1);
        yRange      = yUpperLimit - yLowerLimit;
        yPosition   = y-textPositionY*yRange;

        xUpperLimit = xLimits(2);
        xLowerLimit = xLimits(1);
        xRange      = xUpperLimit - xLowerLimit;
        xPosition   = xLowerLimit+textPositionX*xRange;
        Xposition   = repmat(xPosition, length(y), 1);

        % textHandle is a vector correspond to the line
        textHandles = text(Xposition, yPosition, label, 'Parent', axesHandle);

        for k=1:length(y)
            % Set the text colors to be identical to line colors
            set( textHandles(k), 'color', get(lineHandles(k), 'color') );
        end
    end

    if( holdState==false )
        hold(axesHandle, 'off');
    end
    % this last part is so that it doesn't show up on legends
    set(lineHandles,'tag','hline','handlevisibility','off')

repo: REST


text:
function [Dr] = Dmatrix1D(N,r,V)

% function [Dr] = Dmatrix1D(N,r,V)
% Purpose : Initialize the (r) differentiation matrices on the interval,
%	        evaluated at (r) at order N

Vr = GradVandermonde1D(N, r);
Dr = Vr/V;
return

repo: Aero-matlab


text:
function [ pos ] = scoremap_to_coord( p, locations, crop, locref )

if nargin < 4
    locref = [];
end

scale_factor = p.scale_factor;
stride = p.stride;
half_stride = stride/2;
locref_scale = p.locref_scale;

crd = (locations-1)*stride;
if p.res_net
    crd = crd + half_stride;
end
if ~isempty(locref)
    crd = crd + squeeze(locref)*locref_scale;
end
pos = bsxfun(@plus, crd/scale_factor, double(crop(1:2)-1));

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [omega,u,eqn,info] = fourCurl3P2(node,elem,bdFlag,pde,option)

%% fourCurl3P2: the fourth order curl problem in 3D
%
% [w,u,eqn,info] = fourCurl3P2(node,elem,bdFlag,pde,option)
% uses the first P2 elements to approximate the velocity u and
% the stream function w.
%
% We solve the following equations:
%  -w + curl curl u = 0
%  curl curl w + u  = f
% u\times n = (curl u )\times n  = 0  on \partial \Omega
%
% please check fourCurl3doc for details.
%
% Lin Zhong, June 2013.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Set up optional input arguments
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end

%% Sort elem to asend ordering
elemunSort = elem;
[elem,bdFlag] = sortelem3(elemunSort,bdFlag);

%% Construct Data Structure
[Dlambda,volume] = gradbasis3(node,elem);
volume = abs(volume);
%---------------------------------------------%
% locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
% locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];
locBasesIdx = [1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % phi
               1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % psi
               3 2 4; 3 1 4; 2 1 4; 2 1 3; ...
               4 2 3; 4 1 3; 4 1 2; 3 1 2]; % chi

%---------------------------------------------%
[elem2edge,edge] = dof3edge(elem);
[elem2face,face] = dof3face(elem);
%---------------------------------------------%
NE = size(edge,1);
NF = size(face,1);
NT = size(elem,1);
Ndof = 2*(NE+NF);
%---------------------------------------------%
elem2dof = [elem2edge elem2edge+NE elem2face+2*NE elem2face+2*NE+NF];
face2edge = zeros(NF,3,'int32');
face2edge(elem2face(:,1),:) = elem2edge(:,[4 5 6]);
face2edge(elem2face(:,2),:) = elem2edge(:,[2 3 6]);
face2edge(elem2face(:,3),:) = elem2edge(:,[1 3 5]);
face2edge(elem2face(:,4),:) = elem2edge(:,[1 2 4]);

tic;

%% Assemble Matrix
DiDjcross = zeros(NT,3,4,4);
for i = 1:4
    for j = i+1:4
        DiDjcross(:,:,i,j) = mycross(Dlambda(:,:,i),Dlambda(:,:,j),2);
        DiDjcross(:,:,j,i) = -DiDjcross(:,:,i,j);
    end
end
DiDj = zeros(NT,4,4);
for i = 1:4
    for j = i:4
        DiDj(:,i,j) = dot(Dlambda(:,:,i),Dlambda(:,:,j),2);
        DiDj(:,j,i) = DiDj(:,i,j);
    end
end

ii = zeros(210*NT,1); jj = zeros(210*NT,1);
index = 0;
for i = 1:20
    for j = i:20
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); i3 = locBasesIdx(i,3);
        j1 = locBasesIdx(j,1); j2 = locBasesIdx(j,2); j3 = locBasesIdx(j,3);
        Aij = zeros(NT,1);  Mij = zeros(NT,1);
        %% curl-curl matrix
        if (i<=6) && (j<=6)
            Aij = 4*dot(DiDjcross(:,:,i1,i2),DiDjcross(:,:,j1,j2),2);
        end
        if (i<=6) && (j>12)
            Aij = dot(DiDjcross(:,:,i1,i2),...
            DiDjcross(:,:,j1,j3)-DiDjcross(:,:,j1,j2)+2*DiDjcross(:,:,j2,j3),2)/2;
        end
        if (i>12) && (j>12)
        % curl chi_i =  Dlambda_{i1}mycross phi_{i2,i3} + lambda_{i1}curl phi_{i2,i3}
        % curl chi_j =  Dlambda_{j1}mycross phi_{j2,j3} + lambda_{j1}curl phi_{j2,j3}
        % (Dlambda_{i1}mycross phi_{i2,i3}) dot (Dlambda_{j1}mycross phi_{j2,j3})
            temp11 = ((1+(i2==j2))*dot(DiDjcross(:,:,i1,i3),DiDjcross(:,:,j1,j3),2) ...
                    - (1+(i2==j3))*dot(DiDjcross(:,:,i1,i3),DiDjcross(:,:,j1,j2),2) ...
                    - (1+(i3==j2))*dot(DiDjcross(:,:,i1,i2),DiDjcross(:,:,j1,j3),2) ...
                    + (1+(i3==j3))*dot(DiDjcross(:,:,i1,i2),DiDjcross(:,:,j1,j2),2))/20;
            % lambda_{i1}curl phi_{i2,i3} dot lambda_{j1}curl phi_{j2,j3}
            temp22 = (1+(i1==j1))/5*dot(DiDjcross(:,:,i2,i3),DiDjcross(:,:,j2,j3),2);
            % Dlambda_{i1}mycross phi_{i2,i3} dot lambda_{j1}curl phi_{j2,j3}
            temp12 = dot( (1+(j1==i2))*DiDjcross(:,:,i1,i3) ...
                        - (1+(j1==i3))*DiDjcross(:,:,i1,i2), ...
                                       DiDjcross(:,:,j2,j3),2)/10;
            % Dlambda_{j1}mycross phi_{j2,j3} dot lambda_{i1}curl phi_{i2,i3}
            temp21 = dot( (1+(i1==j2))*DiDjcross(:,:,j1,j3) ...
                        - (1+(i1==j3))*DiDjcross(:,:,j1,j2), ...
                                       DiDjcross(:,:,i2,i3),2)/10;
            Aij = temp11 + temp22 + temp12 + temp21;
        end
        if (i<=6) && (j<=6)
            % block 1: (phi_i,phi_j)
            Mij = 1/20*((1+(i1==j1))*DiDj(:,i2,j2) ...
                      - (1+(i1==j2))*DiDj(:,i2,j1) ...
                      - (1+(i2==j1))*DiDj(:,i1,j2) ...
                      + (1+(i2==j2))*DiDj(:,i1,j1));
        end
        if (i<=6) && (7<=j) && (j<=12)
            % block 2: (psi_j,phi_i)
            Mij = 1/20*( (1+(j1==i1))*DiDj(:,j2,i2) ...
                       - (1+(j1==i2))*DiDj(:,j2,i1) ...
                       + (1+(j2==i1))*DiDj(:,j1,i2) ...
                       - (1+(j2==i2))*DiDj(:,j1,i1));

        end
        if (7<=i) && (i<=12) && (7<=j) && (j<=12)
            % block 3: (psi_j,psi_i)
            Mij = 1/20*((1+(i1==j1))*DiDj(:,i2,j2) ...
                      + (1+(i1==j2))*DiDj(:,i2,j1) ...
                      + (1+(i2==j1))*DiDj(:,i1,j2) ...
                      + (1+(i2==j2))*DiDj(:,i1,j1));
        end
        if (i<=6) && (j>12)
            % block 4: (chi_j,phi_i)
            Mij = intlambda([j1,i1,j2],3)*DiDj(:,i2,j3) ...
                 -intlambda([j1,i1,j3],3)*DiDj(:,i2,j2) ...
                 -intlambda([j1,i2,j2],3)*DiDj(:,i1,j3) ...
                 +intlambda([j1,i2,j3],3)*DiDj(:,i1,j2);
        end
        if (7<=i) && (i<=12) && (j>12)
            % block 5: (chi_j,psi_i)
            Mij = intlambda([j1,i1,j2],3)*DiDj(:,i2,j3) ...
                 -intlambda([j1,i1,j3],3)*DiDj(:,i2,j2) ...
                 +intlambda([j1,i2,j2],3)*DiDj(:,i1,j3) ...
                 -intlambda([j1,i2,j3],3)*DiDj(:,i1,j2);
        end
        if (i>12) && (j>12)
            % block 6: (chi_j,chi_i)
            Mij = intlambda([i1,j1,i2,j2],3)*DiDj(:,i3,j3) ...
                 -intlambda([i1,j1,i2,j3],3)*DiDj(:,i3,j2) ...
                 -intlambda([i1,j1,i3,j2],3)*DiDj(:,i2,j3) ...
                 +intlambda([i1,j1,i3,j3],3)*DiDj(:,i2,j2);
        end
        Aij = Aij.*volume;
        Mij = Mij.*volume;
        sA(index+1:index+NT) = Aij;
        sM(index+1:index+NT) = Mij;
        index = index + NT;
    end
end
clear curlBasis_i curlBasis_j basis_i basis_j
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
A = A + AU + AU';
M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),Ndof,Ndof);
MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),Ndof,Ndof);
M = M + MU + MU';
clear AU MU

% Whole matrix
bigA = [-M A; A M];

%% Righthand side
f = zeros(Ndof,1);
if ~isfield(pde,'f') || (isfield(pde,'f') && isreal(pde.f) && all(pde.f==0))
    pde.f = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 4;   % default order is 4
end
if isfield(pde,'f') && ~isempty(pde.f)
    [lambda,w] = quadpts3(option.fquadorder); % quadrature order is 4
    nQuad = size(lambda,1);
    bt = zeros(NT,20);
    for p = 1:nQuad
        % quadrature points in the x-y-z coordinate
        pxyz = lambda(p,1)*node(elem(:,1),:) ...
             + lambda(p,2)*node(elem(:,2),:) ...
             + lambda(p,3)*node(elem(:,3),:) ...
             + lambda(p,4)*node(elem(:,4),:);
        fp = pde.f(pxyz);
        for k = 1:20
            k1 = locBasesIdx(k,1);
            k2 = locBasesIdx(k,2);
            k3 = locBasesIdx(k,3);
            % evaluate basis at quadrature points
            if k<=6
            % phi_k = lambda_{k1}Dlambda_{k2} - lambda_{k2}Dlambda_{k1};
                basis_k = (lambda(p,k1)*Dlambda(:,:,k2) ...
                          -lambda(p,k2)*Dlambda(:,:,k1));
            elseif k<=12
            % psi_k = lambda_{k1}Dlambda_{k2} + lambda_{k2}Dlambda_{k1};
                basis_k = (lambda(p,k1)*Dlambda(:,:,k2) ...
                          +lambda(p,k2)*Dlambda(:,:,k1));
            else
            % chi_k = lambda_{k1}phi_{k2,k3};
                basis_k = lambda(p,k1)*(lambda(p,k2)*Dlambda(:,:,k3) ...
                                       -lambda(p,k3)*Dlambda(:,:,k2));
            end
            rhs = dot(basis_k,fp,2);
            bt(:,k) = bt(:,k) + w(p)*rhs;
        end
    end
    bt = bt.*repmat(volume,1,20);
    f = accumarray(elem2dof(:),bt(:),[Ndof 1]);
end
clear pxy fp bt rhs basis_k

bigf = [zeros(Ndof,1); f];
assembleTime = toc;

%% Boundary condition
% Find Dirichlet boundary dof: fixedDof
if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N)
    % Dirichlet boundary condition only
    bdFlag = setboundary3(node,elem,'Dirichlet');
end
isBdDof = false(Ndof,1);

if ~isempty(bdFlag)
    %% Dirichlet boundary condition on edge dofs
    % Find boundary faces, edges and nodes
    isBdFace = false(NF,1);
    isBdFace(elem2face(bdFlag(:,1) == 1,1)) = true;
    isBdFace(elem2face(bdFlag(:,2) == 1,2)) = true;
    isBdFace(elem2face(bdFlag(:,3) == 1,3)) = true;
    isBdFace(elem2face(bdFlag(:,4) == 1,4)) = true;
%     boundaryFace = face(isBdFace,:);
    bdFace2edge = face2edge(isBdFace,:);
    isBdEdge = false(NE,1);
    isBdEdge(bdFace2edge(:)) = true;
    edgeBdDof = [find(isBdEdge); NE + find(isBdEdge)];
%     bdEdge = edge(isBdEdge,:);
%     isBdNode(bdEdge) = true;
%     bdNode = find(isBdNode);
    faceBdDof = 2*NE + [find(isBdFace); NF+find(isBdFace)];
%     edgeIdxMap = zeros(NE,1);
%     edgeIdxMap(isBdEdge) = 1:size(bdEdge,1);
%     bdFace2edge = edgeIdxMap(bdFace2edge);
    isBdDof(edgeBdDof) = true;
    isBdDof(faceBdDof) = true;
end

freeDof = ~isBdDof;

if any(freeDof)
    idx = [true(Ndof,1); freeDof];
    bigA_bd= bigA(idx,idx);
    bigf_bd = bigf(idx);
end

%% Solve
t = cputime;
bigu = bigA_bd\bigf_bd;
omega = bigu(1:Ndof);
u = zeros(Ndof,1);
u(freeDof) = bigu(Ndof+1:end);
info.solverTime = cputime -t;
display(info.solverTime);

%% Output information
eqn = struct('elem2edge',elem2edge,'elem2face',elem2face,'face2edge',face2edge,'A',A,'edge',edge,'face',face,'M',M,'f',f);
info.assembleTime = assembleTime;
end

repo: ifem


text:
function F = feval(d, loc, varargin)
%FEVAL     Evaluation functional.
%   This function is deprecated. Use FUNCTIONALBLOCK.FEVAL.

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

F = linop( functionalBlock.feval(loc, double(d), varargin{:}) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef LRN < dagnn.ElementWise
  properties
    param = [5 1 0.0001/5 0.75]
  end

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

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

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

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [R,eff] = randmio_und_signed(W, ITER)
%RANDMIO_UND     Random graph with preserved signed degree distribution
%
%   R       = randmio_und_signed(W,ITER);
%   [R,eff] = randmio_und_signed(W,ITER);
%
%   This function randomizes an undirected network with positively and
%   negatively signed connections, while preserving the positively and
%   negatively signed degree distribution. The function does not preserve
%   the strength distribution in weighted networks.
%
%   Input:      W,      undirected (binary/weighted) connection matrix
%               ITER,   rewiring parameter
%                       (each edge is rewired approximately ITER times)
%
%   Output:     R,      randomized network
%               eff,    number of actual rewirings carried out
%
%   Reference:  Maslov and Sneppen (2002) Science 296:910
%
%
%   2011-2015
%   Dani Bassett, UCSB
%   Olaf Sporns,  Indiana U
%   Mika Rubinov, U Cambridge

%   Modification History:
%   Mar 2011: Original (Dani Bassett, based on randmio_und.m)
%   Mar 2012: Limit number of rewiring attempts,
%             count number of successful rewirings (Olaf Sporns)
%   Dec 2015: Rewritten the core of the rewiring algorithm to allow
%             unbiased exploration of all network configurations. The new
%             algorithm allows positive-positive/negative-negative
%             rewirings, in addition to the previous positive-positive/0-0
%             and negative-negative/0-0 rewirings (Mika Rubinov).

if nargin('randperm')==1
    warning('This function requires a recent (>2011) version of MATLAB.')
end

R     = double(W);              % sign function requires double input
n     = size(R,1);
ITER  = ITER*n*(n-1)/2;

% maximal number of rewiring attempts per 'iter'
maxAttempts = round(n/2);
% actual number of successful rewirings
eff = 0;

for iter=1:ITER
    att=0;
    while (att<=maxAttempts)    %while not rewired
        %select four distinct vertices
        nodes = randperm(n,4);
        a = nodes(1);
        b = nodes(2);
        c = nodes(3);
        d = nodes(4);

        r0_ab = R(a,b);
        r0_cd = R(c,d);
        r0_ad = R(a,d);
        r0_cb = R(c,b);

        %rewiring condition
        if      (sign(r0_ab)==sign(r0_cd)) && ...
                (sign(r0_ad)==sign(r0_cb)) && ...
                (sign(r0_ab)~=sign(r0_ad))

            R(a,d)=r0_ab; R(a,b)=r0_ad;
            R(d,a)=r0_ab; R(b,a)=r0_ad;
            R(c,b)=r0_cd; R(c,d)=r0_cb;
            R(b,c)=r0_cd; R(d,c)=r0_cb;

            eff = eff+1;
            break;
        end %rewiring condition
        att=att+1;
    end %while not rewired
end %iterations

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function imData=bigread2(path_to_file,sframe,num2read)
%reads tiff files in Matlab bigger than 4GB, allows reading from sframe to sframe+num2read-1 frames of the tiff - in other words, you can read page 200-300 without rading in from page 1.
%based on a partial solution posted on Matlab Central (http://www.mathworks.com/matlabcentral/answers/108021-matlab-only-opens-first-frame-of-multi-page-tiff-stack)
%Darcy Peterka 2014, v1.0
%Darcy Peterka 2014, v1.1
%Darcy Peterka 2016, v1.2(bugs to dp2403@columbia.edu)
%Eftychios Pnevmatikakis 2016, v1.3 (added hdf5 support)
%Program checks for bit depth, whether int or float, and byte order.  Assumes uncompressed, non-negative (i.e. unsigned) data.
%
% Usage:  my_data=bigread('path_to_data_file, start frame, num to read);
% "my_data" will be your [M,N,frames] array.
%Will do mild error checking on the inputs - last two inputs are optional -
%if nargin == 2, then assumes second is number of frames to read, and
%starts at frame 1

[~,~,ext] = fileparts(path_to_file);

if strcmpi(ext,'.tiff') || strcmpi(ext,'.tif');

    %get image info
    info = imfinfo(path_to_file);

    % if ~isfield(info,'ImageDescription')
    %    blah=size(info);
    %    numFrames= blah(1);
    % else
    %     he=info.ImageDescription;
    %     numFramesStr = regexp(he, 'images=(\d*)', 'tokens');
    %     numFrames = str2double(numFramesStr{1}{1});
    % end
    blah=size(info);
    numFrames= blah(1);

    num_tot_frames=numFrames;

    %should add more error checking for args... very ugly code below.  works
    %for me after midnight though...
    if nargin<2
        sframe = 1;
    end
    if nargin<3
        num2read=numFrames-sframe+1;
    end
    if sframe<=0
        sframe=1;
    end
    if num2read<1
        num2read=1;
    end
    if sframe>num_tot_frames
        sframe=num_tot_frames;
        num2read=1;
        display('starting frame has to be less than number of total frames...');
    end
    if (num2read+sframe<= num_tot_frames+1)
        lastframe=num2read+sframe-1;
    else
        num2read=numFrames-sframe+1;
        lastframe=num_tot_frames;
        display('Hmmm...just reading from starting frame until the end');
    end

    bd=info.BitDepth;
    he=info.ByteOrder;
    bo=strcmp(he,'big-endian');
    if (bd==64)
        form='double';
    elseif(bd==32)
        form='single'
    elseif (bd==16)
        form='uint16';
    elseif (bd==8)
        form='uint8';
    end

    % Use low-level File I/O to read the file
    fp = fopen(path_to_file , 'rb');
    % The StripOffsets field provides the offset to the first strip. Based on
    % the INFO for this file, each image consists of 1 strip.

    he=info.StripOffsets;
    %finds the offset of each strip in the movie.  Image does not have to have
    %uniform strips, but needs uniform bytes per strip/row.
    idss=max(size(info(1).StripOffsets));
    ofds=zeros(numFrames);
    for i=1:numFrames
        ofds(i)=info(i).StripOffsets(1);
        %ofds(i)
    end
    sframemsg = ['Reading from frame ',num2str(sframe),' to frame ',num2str(num2read+sframe-1),' of ',num2str(num_tot_frames), ' total frames'];
    disp(sframemsg)
    pause(.2)
    %go to start of first strip
    fseek(fp, ofds(1), 'bof');
    %framenum=numFrames;
    framenum=num2read;
    imData=cell(1,framenum);

    he_w=info.Width;
    he_h=info.Height;
    % mul is set to > 1 for debugging only
    mul=1;
    if strcmpi(form,'uint16') || strcmpi(form,'uint8')
        if(bo)
            for cnt = sframe:lastframe
                %cnt;
                fseek(fp,ofds(cnt),'bof');
                tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-be')';
                imData{cnt-sframe+1}=cast(tmp1,form);
            end
        else
            for cnt=sframe:lastframe
                % cnt;
                fseek(fp,ofds(cnt),'bof');
                tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-le')';
                imData{cnt-sframe+1}=cast(tmp1,form);
            end
        end
    elseif strcmpi(form,'single')
        if(bo)
            for cnt = sframe:lastframe
                %cnt;
                fseek(fp,ofds(cnt),'bof');
                tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-be')';
                imData{cnt-sframe+1}=cast(tmp1,'single');
            end
        else
            for cnt = sframe:lastframe
                %cnt;
                fseek(fp,ofds(cnt),'bof');
                tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-le')';
                imData{cnt-sframe+1}=cast(tmp1,'single');
            end
        end
        elseif strcmpi(form,'double')
            if(bo)
                for cnt = sframe:lastframe
                    %cnt;
                    fseek(fp,ofds(cnt),'bof');
                    tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-be.l64')';
                    imData{cnt-sframe+1}=cast(tmp1,'single');
                end
            else
                for cnt = sframe:lastframe
                    %cnt;
                    fseek(fp,ofds(cnt),'bof');
                    tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-le.l64')';
                    imData{cnt-sframe+1}=cast(tmp1,'single');
                end
            end
    end
            %ieee-le.l64

        imData=cell2mat(imData);
        imData=reshape(imData,[he_h*mul,he_w,framenum]);
        fclose(fp);
        display('Finished reading images')
elseif strcmpi(ext,'.hdf5') || strcmpi(ext,'.h5');
    info = hdf5info(path_to_file);
    dims = info.GroupHierarchy.Datasets.Dims;
    if nargin < 2
        sframe = 1;
    end
    if nargin < 3
        num2read = dims(end)-sframe+1;
    end
    num2read = min(num2read,dims(end)-sframe+1);
    imData = h5read(path_to_file,'/mov',[ones(1,length(dims)-1),sframe],[dims(1:end-1),num2read]);
else
    error('Unknown file extension. Only .tiff and .hdf5 files are currently supported');
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [dist,PP0] = pointTriangleDistance(TRI,P)
% POINTTRIANGLEDISTANCE  Calculate distance between a point and a triangle
% in 3D
%
% SYNTAX
%   dist = pointTriangleDistance(TRI,P)
%   [dist,PP0] = pointTriangleDistance(TRI,P)
%
% DESCRIPTION
%   Calculate the distance of a given point P from a triangle TRI.
%   Point P is a row vector of the form 1x3. The triangle is a matrix
%   formed by three rows of points TRI = [P1;P2;P3] each of size 1x3.
%   dist = pointTriangleDistance(TRI,P) returns the distance of the point P
%   to the triangle TRI.
%   [dist,PP0] = pointTriangleDistance(TRI,P) additionally returns the
%   closest point PP0 to P on the triangle TRI.

% Author: Gwendolyn Fischer
% Release: 1.0
% Release date: 09/02/02
% Release: 1.1 Fixed Bug because of normalization
% Release: 1.2 Fixed Bug because of typo in region 5 20101013
% Release: 1.3 Fixed Bug because of typo in region 2 20101014

% Minor modifications by Ramon Casero <rcasero@gmail.com> for the Gerardus
% Project.
% Version: 0.1.0

% Possible extention could be a version tailored not to return the distance
% and additionally the closest point, but instead return only the closest
% point. Could lead to a small speed gain.

% Example:
% %% The Problem
% P0 = [0.5 -0.3 0.5];
%
% P1 = [0 -1 0];
% P2 = [1  0 0];
% P3 = [0  0 0];
%
% vertices = [P1; P2; P3];
% faces = [1 2 3];
%
% %% The Engine
% [dist,PP0] = pointTriangleDistance([P1;P2;P3],P0);
%
% %% Visualization
% [x,y,z] = sphere(20);
% x = dist*x+P0(1);
% y = dist*y+P0(2);
% z = dist*z+P0(3);
%
% figure
% hold all
% patch('Vertices',vertices,'Faces',faces,'FaceColor','r','FaceAlpha',0.8);
% plot3(P0(1),P0(2),P0(3),'b*');
% plot3(PP0(1),PP0(2),PP0(3),'*g')
% surf(x,y,z,'FaceColor','b','FaceAlpha',0.3)
% view(3)

% The algorithm is based on
% "David Eberly, 'Distance Between Point and Triangle in 3D',
% Geometric Tools, LLC, (1999)"
% http:\\www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
%
%        ^t
%  \     |
%   \reg2|
%    \   |
%     \  |
%      \ |
%       \|
%        *P2
%        |\
%        | \
%  reg3  |  \ reg1
%        |   \
%        |reg0\
%        |     \
%        |      \ P1
% -------*-------*------->s
%        |P0      \
%  reg4  | reg5    \ reg6

%% Do some error checking
if nargin<2
  error('pointTriangleDistance: too few arguments see help.');
end
P = P(:)';
if size(P,2)~=3
  error('pointTriangleDistance: P needs to be of length 3.');
end

if size(TRI)~=[3 3]
  error('pointTriangleDistance: TRI needs to be of size 3x3.');
end

% ToDo: check for colinearity and/or too small triangles.

% rewrite triangle in normal form
B = TRI(1,:);
E0 = TRI(2,:)-B;
%E0 = E0/sqrt(sum(E0.^2)); %normalize vector
E1 = TRI(3,:)-B;
%E1 = E1/sqrt(sum(E1.^2)); %normalize vector

D = B - P;
a = dot(E0,E0);
b = dot(E0,E1);
c = dot(E1,E1);
d = dot(E0,D);
e = dot(E1,D);
f = dot(D,D);

det = a*c - b*b; % do we have to use abs here?
s   = b*e - c*d;
t   = b*d - a*e;

% Terible tree of conditionals to determine in which region of the diagram
% shown above the projection of the point into the triangle-plane lies.
if (s+t) <= det
  if s < 0
    if t < 0
      %region4
      if (d < 0)
        t = 0;
        if (-d >= a)
          s = 1;
          sqrDistance = a + 2*d + f;
        else
          s = -d/a;
          sqrDistance = d*s + f;
        end
      else
        s = 0;
        if (e >= 0)
          t = 0;
          sqrDistance = f;
        else
          if (-e >= c)
            t = 1;
            sqrDistance = c + 2*e + f;
          else
            t = -e/c;
            sqrDistance = e*t + f;
          end
        end
      end %of region 4
    else
      % region 3
      s = 0;
      if e >= 0
        t = 0;
        sqrDistance = f;
      else
        if -e >= c
          t = 1;
          sqrDistance = c + 2*e +f;
        else
          t = -e/c;
          sqrDistance = e*t + f;
        end
      end
    end %of region 3
  else
    if t < 0
      % region 5
      t = 0;
      if d >= 0
        s = 0;
        sqrDistance = f;
      else
        if -d >= a
          s = 1;
          sqrDistance = a + 2*d + f;% GF 20101013 fixed typo d*s ->2*d
        else
          s = -d/a;
          sqrDistance = d*s + f;
        end
      end
    else
      % region 0
      invDet = 1/det;
      s = s*invDet;
      t = t*invDet;
      sqrDistance = s*(a*s + b*t + 2*d) ...
                  + t*(b*s + c*t + 2*e) + f;
    end
  end
else
  if s < 0
    % region 2
    tmp0 = b + d;
    tmp1 = c + e;
    if tmp1 > tmp0 % minimum on edge s+t=1
      numer = tmp1 - tmp0;
      denom = a - 2*b + c;
      if numer >= denom
        s = 1;
        t = 0;
        sqrDistance = a + 2*d + f; % GF 20101014 fixed typo 2*b -> 2*d
      else
        s = numer/denom;
        t = 1-s;
        sqrDistance = s*(a*s + b*t + 2*d) ...
                    + t*(b*s + c*t + 2*e) + f;
      end
    else          % minimum on edge s=0
      s = 0;
      if tmp1 <= 0
        t = 1;
        sqrDistance = c + 2*e + f;
      else
        if e >= 0
          t = 0;
          sqrDistance = f;
        else
          t = -e/c;
          sqrDistance = e*t + f;
        end
      end
    end %of region 2
  else
    if t < 0
      %region6
      tmp0 = b + e;
      tmp1 = a + d;
      if (tmp1 > tmp0)
        numer = tmp1 - tmp0;
        denom = a-2*b+c;
        if (numer >= denom)
          t = 1;
          s = 0;
          sqrDistance = c + 2*e + f;
        else
          t = numer/denom;
          s = 1 - t;
          sqrDistance = s*(a*s + b*t + 2*d) ...
                      + t*(b*s + c*t + 2*e) + f;
        end
      else
        t = 0;
        if (tmp1 <= 0)
            s = 1;
            sqrDistance = a + 2*d + f;
        else
          if (d >= 0)
              s = 0;
              sqrDistance = f;
          else
              s = -d/a;
              sqrDistance = d*s + f;
          end
        end
      end
      %end region 6
    else
      % region 1
      numer = c + e - b - d;
      if numer <= 0
        s = 0;
        t = 1;
        sqrDistance = c + 2*e + f;
      else
        denom = a - 2*b + c;
        if numer >= denom
          s = 1;
          t = 0;
          sqrDistance = a + 2*d + f;
        else
          s = numer/denom;
          t = 1-s;
          sqrDistance = s*(a*s + b*t + 2*d) ...
                      + t*(b*s + c*t + 2*e) + f;
        end
      end %of region 1
    end
  end
end

% account for numerical round-off error
if (sqrDistance < 0)
  sqrDistance = 0;
end

dist = sqrt(sqrDistance);

if nargout>1
  PP0 = B + s*E0 + t*E1;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function SO3VF = minus(SO3VF1,SO3VF2)
% overloads |SO3VF1 - SO3VF2|
%
% Syntax
%   SO3VF = SO3VF1 - SO3VF2
%   SO3VF = a - SO3VF1
%   SO3VF = SO3VF1 - a
%
% Input
%  SO3VF1, SO3VF2 - @SO3VectorField
%  a - double
%
% Output
%  SO3VF - @SO3VectorField
%

SO3VF = plus(SO3VF1,-SO3VF2);

end

repo: mtex


text:
function [Gal] = G2Gal(G)
% Convert acceleration from average acceration due to Earth's gravity to
% galileos
% Chad A. Greene 2012
Gal = G*9.80665e+2;
end


text:
function stroud_test34 ( )

%*****************************************************************************80
%
%% TEST34 tests SPHERE_UNIT_AREA_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST34\n' );
  fprintf ( 1, '  In N dimensions:\n' );
  fprintf ( 1, '  SPHERE_UNIT_AREA_ND computes the area of \n' );
  fprintf ( 1, '  the unit sphere;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N    Area\n' );
  fprintf ( 1, '\n' );

  for n = 2 : 10
    fprintf ( 1, '  %2d  %12f\n', n, sphere_unit_area_nd ( n ) );
  end

  return
end


text:
function stackOut = WardAlignment(stack, bStackOut, folder_name, format, target_exposure)
%
%
%       stackOut = WardAlignment(stack, bStackOut, folder_name, format)
%
%       This function shifts pixels on the right with wrapping of the moved
%       pixels. This can be used as rotation on the Y-axis for environment
%       map encoded as longituted-latitude encoding.
%
%       Input:
%           -stack: a stack (4D) containing all images.
%           -bStackOut: if it is true (1), it outputs an aligned stack in
%           stackOut. Otherwise, stackOut is set to empty, [].
%           -folder_name: the folder name where the stack is stored. This flag
%           is valid if stack is empty, [].
%           -format: the file format of the stack. This flag is valid if
%           stack is empty, [].
%           -target_exposure: The index of the target exposure for aligning
%           images. If stack is empty, [], it contains the name of the file
%           for the alignment. If not provided the stack will be analyzed.
%
%       Output:
%           -stackOut: the aligned stack as output
%
%     Copyright (C) 2012-15  Francesco Banterle
%
%     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 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('folder_name', 'var'))
    folder_name = '';
end

if(~exist('format', 'var'))
    format = '';
end

lst = [];

bStack = ~isempty(stack);

if(~bStack)
    lst = dir([folder_name, '/*.', format]);
    n = length(lst);
else
    n = size(stack, 4);
    stack = normalizeFromAnything(stack);
end

if(n < 2)
    return;
end

if(~exist('target_exposure', 'var'))
    target_exposure = GalloReferenceImage(stack, folder_name, format);
else
    if(~bStack)
        target_exposure = findNameInList(lst, target_exposure);
    end
end

if(bStack)
    img = stack(:,:,:,target_exposure);
else
    img = ldrimread([folder_name, '/', lst(target_exposure).name]);
end

[r, c, col] = size(img);

alignment = zeros(n, 2);

stackOut = [];

if(bStackOut)
    stackOut = zeros(r, c, col, n);
    stackOut(:,:,:,target_exposure) = img;
end

for i=1:n
    shift_ret = [0, 0];

    if(i ~= target_exposure)
        disp(['Aligning image ', num2str(i), ' to image ', num2str(target_exposure)]);

        if(~bStack)
            img_work = ldrimread([folder_name, '/', lst(i).name]);
        else
            img_work = stack(:,:,:,i);
        end

        img_work_shifted = WardImageAlignment(img, img_work);

        if(bStackOut)
            stackOut(:,:,:,i) = img_work_shifted;
        end

        if(~bStack)
            name = strrep(lst(i).name, ['.', format], ['_shifted.', format]);

            if(strcmp(lst(i).name, name) == 1)
                name = [name, '_shifted.', format];
            end

            imwrite(img_work_shifted,[folder_name, '/', name]);
        end

        clear('img_work_shifted');
        clear('img_work');
    end

    alignment(i,:) = shift_ret;
end

end

author: banterle


text:
function bbci_log_close(data)
%BBCI_LOG_CLOSE - Close log file(s)
%
%Synopsis:
%  bbci_log_close(DATA)
%
%Arguments:
% DATA - DATA strcut as in bbci_apply, bbci_calibrate. The subfield
%        DATA.log.fid holds the file identifiers of the log files.

% 12-2011 Benjamin Blankertz

if isfield(data, 'log'),
  fid_to_close= data.log.fid(find(data.log.fid>2));
  for fid= fid_to_close,
    fclose(fid);
  end
end

author: bbci


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% 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 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function nw = resample_rotation_rate(w, old_time, new_time)

dT = old_time(2:end) - old_time(1:end-1);
theta = cumsum( w(2:end) .* dT );

new_theta = interp1(old_time(2:end), theta(:,1), new_time);

%{
figure(1);
plot(old_time(2:end), theta, 'rx-'); hold on;
plot(new_time, new_theta, 'bo-'); hold off;
%}

nw = (new_theta(2:end) - new_theta(1:end-1)) ./ (new_time(2:end) - new_time(1:end-1));
nw(isnan(nw)) = 0;
nw = [0; nw];

%{
figure(2);
plot(old_time, w(:,1), 'rx-'); hold on;
plot(new_time, nw, 'bo-'); hold off;
%}

author: alex-golts

repo: Video-Stabilization


text:
function op = prox_l1_deadzone( q, epsilon )

%PROX_L1_DEADZONE    L1 norm with deadzone [-eps,eps]
%    OP = PROX_L1_DEADZONE( q, eps ) implements the nonsmooth function
%        OP(X) = q*sum_i max(0,x-eps)+max(0,-x-eps) = q*sum_i
%        max(0,|x|-eps)
%    Q is optional; if omitted, Q=1 is assumed. But if Q is supplied,
%    then it must be a positive real scalar (or must be same size as X).
%

% New March 2 2016

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) ||  any( q(:) < 0 ) || all(q(:)==0) %|| numel( q ) ~= 1
    if q==0
        op = prox_0;
        warning('TFOCS:zeroQ','q=0 so returning the proximal operator for the zero function');
        return;
    else
        error( 'Argument must be positive.' );
    end
end

% This is Matlab and Octave compatible code
op = tfocs_prox( @(x)f(q, epsilon,x), @(x,t)prox_f(q, epsilon,x,t) , 'vector' );
end

% These are now subroutines, that are NOT in the same scope
function v = f(qq,epsilon, x)
    v = sum( qq(:).*max(0,abs(x(:))-epsilon), 1 );
end

function x = prox_f(qq,epsilon, x,t)
% p = x if |x|<eps,
% p = sign(x)*eps if eps < |x| < tq
% p = sign(x)*( |x| - (eps+tq) ) if |x| > tq
% i.e., if |x|>eps,  p = sign(x).*max( |x|-tq, eps )
    tq = t .* qq;
    x  = sign(x).*( abs(x).*( abs(x)<epsilon ) + max(abs(x)-tq, epsilon).*( abs(x)>= epsilon) );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [ pl, pld ] = p12_legendre_val ( t, dtdx, npolys )

%*****************************************************************************80
%
%% P12_LEGENDRE_VAL evaluates the Legendre polynomials and derivatives.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T, the argument of the Legendre polynomials, in
%    the normalized interval [-1,1].
%
%    Input, real DTDX, the value of the quantity dTdX at the point
%    X.  In the most common case, this is simply the relationship
%    between the width of the normalized T interval (2), and the
%    width of the X interval to which the Legendre polynomial
%    arguments have been mapped.  DTDX is needed so that the
%    computed values PLD can be converted from dPL/dT to dPL/dX.
%
%    Input, integer NPOLYS, the number of Legendre polynomials to
%    evaluate.  If NPOLYS is 1, then only the constant polynomial
%    is evaluated,  NPOLYS = 2 means the constant and linear, and so on.
%
%    Output, real PL(NPOLYS), PLD(NPOLYS), the values of PL(X)
%    and dPL(X)/dX at the point X which has normalized coordinate T.
%
  if ( 1 <= npolys )
    pl(1) = 1.0;
    pld(1) = 0.0;
  end

  if ( 2 <= npolys )
    pl(2) = t;
    pld(2) = 1.0;
  end

  a = 0.0;
  for i = 3 : npolys

    a = a + 1.0;

    pl(i) = ( ( 2.0 * a + 1.0 ) * t * pl(i-1) - a * pl(i-2) ) / ( a + 1.0 );

    pld(i) = ( ( 2.0 * a + 1.0 ) * ( t * pld(i-1) + pl(i-1) ) ...
           - a * pld(i-2) ) / ( a + 1.0 );

  end

  pld(1:npolys) = dtdx * pld(1:npolys);

  return
end


text:
 % Copyright 2001, Brown University, Providence, Rhode Island.
 %
 % All Rights Reserved
 %
 % Permission to use this software for noncommercial research and
 % educational purposes is hereby granted without fee.
 % Redistribution, sale, or incorporation of this software into a
 % commercial product is prohibited.
 %
 % BROWN UNIVERSITY DISCLAIMS ANY AND ALL WARRANTIES WITH REGARD TO
 % THIS SOFTWARE,INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
 % AND FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL BROWN
 % UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 % DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE,
 % DATA OR PROFITS.

function [dmodedr, dmodeds] = am282jacobideriv2dab(a,b,ID,JD)

  dims = size(a);

  j2deriv = zeros(dims);

  fa  = am282jacobi1d(a,ID,0,0);
  dfa = am282jacobideriv1d(a,ID,0,0);

  gb  = am282jacobi1d(b,JD,2*ID+1,0);
  dgb = am282jacobideriv1d(b,JD,2*ID+1,0);

  % r-derivative
  % d/dr = da/dr d/da + db/dr d/db = (2/(1-s)) d/da = (2/(1-b)) d/da
  dmodedr = dfa.*gb;
  if(ID>0)
   dmodedr = dmodedr.*((0.5*(1-b)).^(ID-1));
  end

  % s-derivative
  % d/ds = da/ds d/da + db/ds d/db = (2(r+1)/(1-s)^2) d/da + d/db
  %                                = ((1+a)/2)/((1-b)/2) d/da + d/db
  dmodeds = dfa.*(gb.*(0.5*(1+a)));
  if(ID>0)
   dmodeds = dmodeds.*((0.5*(1-b)).^(ID-1));
  end

  tmp = dgb.*((0.5*(1-b)).^ID);
  if(ID>0)
   tmp = tmp-0.5*ID*gb.*((0.5*(1-b)).^(ID-1));
  end

  dmodeds = dmodeds+fa.*tmp;

   % normalise
  dmodedr = dmodedr*sqrt((ID+0.5)*(ID+JD+1));
  dmodeds = dmodeds*sqrt((ID+0.5)*(ID+JD+1));


text:
%kklocxform 'Rotate, Translate, and Scale Location Data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros klocxform.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'the input object to be processed'
% Double: xrot ' X ', default: 0: 'rotate around x axis'
% Double: yrot ' Y ', default: 0: 'rotate around y axis'
% Double: zrot ' Z ', default: 0: 'rotate around z axis'
% Toggle: fx 'Flip X', default: 0: 'flip x coordinates'
% Double: xtrans ' X ', default: 0: 'translate in x direction'
% Toggle: fy 'Flip Y', default: 0: 'flip y coordinates'
% Double: ytrans ' Y ', default: 0: 'translate in y direction'
% Toggle: fz 'Flip Z', default: 0: 'flip z coordinates'
% Double: ztrans ' Z ', default: 0: 'translate in z direction'
% OutputFile: o 'Output', required: 'the output object'
% Double: scale ' XYZ ', default: 1: 'orthogonal scale factor'
% Double: xscale ' X ', default: 1: 'scale x coordinates'
% Double: yscale ' Y ', default: 1: 'scale y coordinates'
% Double: zscale ' Z ', default: 1: 'scale z coordinates'
%
% Example: o = kklocxform(i, {'i','';'xrot',0;'yrot',0;'zrot',0;'fx',0;'xtrans',0;'fy',0;'ytrans',0;'fz',0;'ztrans',0;'o','';'scale',1;'xscale',1;'yscale',1;'zscale',1})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% klocxform - Rotate, Translate, and Scale Location Data
%
%  DESCRIPTION
% This routine performs a homogeneous coordinate
% transformation on the location data in the source
% object, and saves the transformed location data into
% the destination object.
%
% The source object must have location data, and the
% destination object must not have location data for
% this routine to function properly.  This routine will
% create location data of the proper grid type in the
% destination.
%
% The coordinate transformation can be specfied by nine
% parameters which specify an x-y-z rotation, an x-y-z
% translation, and an x-y-z scaling.  The rotations are
% all performed with respect to the (0,0,0) origin.  A
% scaling matrix, a translation matrix, and three
% rotation matrices are constructed from these
% specified parameters and composited into a single
% transformation matrix in the following order :
%
%    [1] Scaling
%    [2] X Rotation
%    [3] Y Rotation
%    [4] Z Rotation
%    [5] Translation
%
% A value of 1.0 should be used for any scaling
% parameters where no scaling is desired.  A value of
% 0.0 should be used for any rotation parameters where
% no rotation is desired or for any translation
% parameters where no translation is desired.
%
% The location points are interpreted to be in
% three-space for the transformations.  Explicit
% location data located in dimensions higher than three
% will be lost during the transformation.  The
% transformation will be performed using the datatype
% KFLOAT.
%
% In general, the grid type of the source object will
% be preserved in the destination object.  For example,
% uniform location data, when transformed and stored in
% the destination object, will still be uniform
% location data.  The exception to this rule is
% rectilinear data, which must be converted to
% curvilinear data in the event of a rotation.  In this
% case the destination location data will be stored as
% explicit curvilinear data points.
%
% Reference : Foley & VanDam "Computer Graphics :
% Principles and Practice", Second Edition, 1990,
% Chapter 5.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kklocxform(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kklocxform(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'xrot', 0;'yrot', 0;'zrot', 0;'fx', 0;'xtrans', 0;'fy', 0;'ytrans', 0;'fz', 0;'ztrans', 0;'o', '__output';'scale', 1;'xscale', 1;'yscale', 1;'zscale', 1};
maxval={0,360,360,360,0,0,0,0,0,0,0,2,2,2,2};
minval={0,-360,-360,-360,0,0,0,0,0,0,0,2,2,2,2};
istoggle=[0,1,1,1,1,1,1,1,1,1,0,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','Double','Double','Double','Toggle','Double','Toggle','Double','Toggle','Double','OutputFile','Double','Double','Double','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'klocxform"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function qwgw_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests QWGW for the Gegenbauer weight.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2014
%
%  Author:
%
%    John Burkardt
%

%
%  Set the quadrature interval and number of points.
%
  a = -1.0;
  b = +1.0;
  n = 5;
%
%  Set the weight function parameter.
%
  alpha = 0.25;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03:\n' );
  fprintf ( 1, '  Compute points and weights for Gauss quadrature\n' );
  fprintf ( 1, '  with the Gegenbauer weight w(x) = (1-x^2)^alpha.\n' );
  fprintf ( 1, '  Order N = %d\n', n );
  fprintf ( 1, '  ALPHA = %g\n', alpha );
  fprintf ( 1, '  Interval = [%g,%g]\n', a, b );
%
%  Set the recursion coefficients.
%
  aj = zeros ( n, 1 );
  bj = zeros ( n, 1 );

  aj(1:n) = 0.0;

  for j = 1 : n - 1
    jr = j;
    bj(j) = ( jr * ( 2.0 * alpha + jr ) ) ...
      / ( 4.0 * ( alpha + jr )^2 - 1.0 );
  end
  bj(n) = 0.0;

  bj(1:n) = sqrt ( bj(1:n) );

  mu0 = gamma ( alpha + 1.0 ) * gamma ( 0.5 ) / gamma ( alpha + 1.5 );
%
%  Compute the points and weights.
%
  [ x, w ] = sgqf ( n, aj, bj, mu0 );

  r8vec_print ( n, x, '  Abscissas:' );
  r8vec_print ( n, w, '  Weights:' );

  return
end


text:
function bnet = mk_car_bnet()
% MK_CAR_BNET Make the car trouble-shooter bayes net.
%
% This network is from p13 of "Troubleshooting under uncertainty", Heckerman, Breese and
% Rommelse, Microsoft Research Tech Report 1994.

BatteryAge = 1;
Battery = 2;
Starter = 3;
Lights = 4;
TurnsOver = 5;
FuelPump = 6;
FuelLine = 7;
FuelSubsys =8;
Fuel = 9;
Spark = 10;
Starts = 11;
Gauge = 12;

n = 12;
dag = zeros(n);
dag(1,2) = 1;
dag(2,[4 5])=1;
dag(3,5) = 1;
dag(6,8) = 1;
dag(7,8) = 1;
dag(8,11) = 1;
dag(9,12) = 1;
dag(10,11) = 1;

arity = 2;
ns = arity*ones(1,n);
bnet = mk_bnet(dag, ns);
for i=1:n
    bnet.CPD{i} = tabular_CPD(bnet, i);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [laglength, AIC, logL] = VARlag(ENDO,maxlag,const,EXOG,lag_ex)
% =======================================================================
% Determine VAR lag length wit Akaike criterion
% =======================================================================
% [laglength, AIC, logL] = VARlag(ENDO,maxlag,const,EXOG,lag_ex)
% -----------------------------------------------------------------------
% INPUT
%   - ENDO: an (nobs x nvar) matrix of endogenous variables.
%	- maxlag: the maximum lag length over which Akaike information criterion
%       is computed
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - const: 0 no const; 1 const ; 2 const&trend; 3 const&trend^2; (dflt = 1)
%	- EXOG: optional matrix of variables (nobs x nvar_ex)
%   - nlag_ex: number of lags for exogeonus variables (dflt = 0)
% -----------------------------------------------------------------------
% OUTPUT
%	- laglength: preferred lag lenghth according to Akaike information criterion
%   - AIC: vector [maxlag x 1] of Akaike information criterion
%   - logL: vector [maxlag x 1] of loglikelihood
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% Note:
% The determinant of the residual covariance is computed without adjusting
% for the degrees of freedom. The log likelihood value is computed assuming
% a multivariate normal (Gaussian) distribution

%% Check inputs
%=========================================================
[nobs, ~] = size(ENDO);

% Check if ther are constant, trend, both, or none
if ~exist('const','var')
    const = 1;
end

% Check if there are exogenous variables
if exist('EXOG','var')
    [nobs2, num_ex] = size(EXOG);
    % Check that ENDO and EXOG are conformable
    if (nobs2 ~= nobs)
        error('var: nobs in EXOG-matrix not the same as y-matrix');
    end
    clear nobs2
else
    num_ex = 0;
end

% Check if there is lag order of EXOG, otherwise set it to 0
if ~exist('lag_ex','var')
    lag_ex = 0;
end

% number of exogenous variables per equation
nvar_ex = num_ex*(lag_ex+1);

%% Compute log likelihood and Akaike criterion
%=========================================================
logL = zeros(maxlag,1);
AIC  = zeros(maxlag,1);
for i=1:maxlag
    X = ENDO(maxlag+1-i:end,:);
    aux = VARmodel(X,i,const);
    if nvar_ex>0
        Y = EXOG(maxlag+1-i:end,:);
        aux = VARmodel(X,i,const,Y,lag_ex);
    end
    NOBS = aux.nobs;
    NEQS = aux.nvar;
    NTOTCOEFF = aux.ntotcoeff;
    RES = aux.residuals;
    SIGMA = (1/(NOBS)).*(RES)'*(RES); %use non dof adjusted version of sigma as in Eviews
    logL(i) = -(NOBS/2)* (NEQS*(1+log(2*pi)) + log(det(SIGMA)));
    AIC(i)  = -2*(logL(i)/NOBS) + 2*(NEQS*NTOTCOEFF/NOBS);
end

laglength = find(AIC==min(AIC));

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox