text:
classdef VNT4 < PROBLEM
% <multi> <real> <constrained>
% Benchmark MOP proposed by Viennet

%------------------------------- Reference --------------------------------
% R. Viennet, C. Fonteix, and I. Marc, Multicriteria optimization using a
% genetic algorithm for determining a Pareto set, International Journal of
% Systems Science, 1996, 27(2): 255-260.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M        = 3;
            obj.D        = 2;
            obj.lower    = [-4,-4];
            obj.upper    = [4,4];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj(:,1) = (PopDec(:,1)-2).^2/2 + (PopDec(:,2)+1).^2/13 + 3;
            PopObj(:,2) = (PopDec(:,1)+PopDec(:,2)-3).^2/175 + (2*PopDec(:,2)-PopDec(:,1)).^2/17 - 13;
            PopObj(:,3) = (3*PopDec(:,1)-2*PopDec(:,2)+4).^2/8 + (PopDec(:,1)-PopDec(:,2)+1).^2/27 + 15;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            PopCon(:,1) = PopDec(:,2) + 4*PopDec(:,1) - 4;
            PopCon(:,2) = -1 - PopDec(:,1);
            PopCon(:,3) = PopDec(:,1) - 2 - PopDec(:,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            N = ceil(sqrt(N));
            x = linspace(-1,1.2,N);
            X = [];
            for i = 1 : N
                X = [X;repmat(x(i),N,1),linspace(x(i)-2,-4*x(i)+4,N)'];
            end
            R = obj.CalObj(X);
            R = R(NDSort(R,1)==1,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            x = linspace(-1,1.2,40);
            X = [];
            for i = 1 : 40
                X = [X;repmat(x(i),40,1),linspace(x(i)-2,-4*x(i)+4,40)'];
            end
            R = obj.CalObj(X);
            R(NDSort(R,1)>1,:) = nan;
            R = {reshape(R(:,1),40,40),reshape(R(:,2),40,40),reshape(R(:,3),40,40)};
        end
    end
end

author: BIMK

repo: PlatEMO


text:
function q = qinv(p)

% QINV  quaternion inverse
%
%   Q = QINV(P) returns a quaternion Q which is the quaternion inverse of
%     quaternion P.
%     - P is a quaternion. It is a 4-vector or a 4*N array (column i
%        represents quaternion i) where N is the number of quaternions.
%     - Q is the quaternion inverse of quaternion P. It is a 4*N array.

sp = size(p);
if sp == [1 4]
    p = p.';
    sp = size(p);
end

% wrong format
if sp(1) ~= 4
    error('DualQuaternion:Qinv:wrongsize',...
        '%d rows in the P array. It should be 4.',sp(1));
end

normp = qnorm(p);
q = qconj(p)./repmat(normp.^2,4,1);

author: Sable

repo: mcbench-benchmarks


text:
function spiral_test01 ( )

%*****************************************************************************80
%
%% SPIRAL_TEST01 generates a field and estimates its range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPIRAL_TEST01:\n' );
  fprintf ( 1, '  Sample a spiral velocity field and estimate\n' );
  fprintf ( 1, '  the range of the solution values.\n' );

  n = 1000;

  x_lo = 0.0;
  x_hi = 1.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, x_lo, x_hi, seed );
  y_lo = 0.0;
  y_hi = 1.0;
  [ y, seed ] = r8vec_uniform_ab ( n, y_lo, y_hi, seed );

  c = 1.00;

  [ u, v ] = uv_spiral ( n, x, y, c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '           Minimum       Maximum\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  U:  %14.6g  %14.6g\n', min ( u ), max ( u ) );
  fprintf ( 1, '  V:  %14.6g  %14.6g\n', min ( v ), max ( v ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = SC_MMA(y,doOverlap,scaleRange,qRange)
% SC_MMA   Physionet implementation of multiscale multifractal analysis
%
% Scale-dependent estimates of multifractal scaling in a time series.

% ------------------------------------------------------------------------------
% Modified by Ben Fulcher for use in hctsa, 2015-05-12.
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Copyright (C) 2014 Jan Gieraltowski
% ------------------------------------------------------------------------------
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE.  See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place - Suite 330, Boston, MA 02111-1307, USA.
%
% Author: Jan Gieraltowski
% Warsaw University of Technology
% Faculty of Physics
% gieraltowski@if.pw.edu.pl
% http://gieraltowski.fizyka.pw.edu.pl/
%
% Method was first proposed in:
% J. Gieraltowski, J. J. Zebrowski, and R. Baranowski,
% Multiscale multifractal analysis of heart rate variability recordings
% with a large number of occurrences of arrhythmia,
% Phys. Rev. E 85, 021915 (2012).
% http://dx.doi.org/10.1103/PhysRevE.85.021915
%
% Please cite the above publication when referencing this material,
% and also include the standard citation for PhysioNet:
% Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG,
% Mietus JE, Moody GB, Peng C-K, Stanley HE.
% PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource
% for Complex Physiologic Signals.
% Circulation 101(23):e215-e220
% [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215]; 2000 (June 13).
% ------------------------------------------------------------------------------

% Generate a plot:
doPlot = false;

% Time-series length:
N = length(y);

% ------------------------------------------------------------------------------
% Check Inputs:
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(doOverlap)
    doOverlap = false;
    % 0 - time series is partitioned into non overlapping windows of analysis,
    % 1- time series is partitioned into overlapping windows of analysis, step between consecutive windows is = 1 (much longer calculations)
end

if nargin < 3 || isempty(scaleRange)
    scaleRange = [10,round(N/40)];
    % minimal s scale used, when calculating Fq(s) functions family (default 10)
    % maximal s scale used, when calculating Fq(s) functions family, has to be multiple of 5 (default 600; in general should be near to N/50, where N is a time series length)
end
minScale = scaleRange(1);
maxScale = scaleRange(2);
if (maxScale/5) < minScale
    warning('Time-series (N=%u) too short for multiscale multifractal analysis',N);
    out = NaN;
    return
elseif rem(maxScale,5)~=0
    maxScale = round(maxScale/5)*5;
    fprintf(1,'adjusted maxScale to %u\n',maxScale);
end

if nargin < 4 || isempty(qRange)
    qRange = [-5,5];
    % minimal/maximal multifractal parameter q used (default -5)
end
qMin = qRange(1);
qMax = qRange(2);

qList = qMin:0.1:qMax;
qList(qList == 0) = 0.0001;

% if nargin < 5
%     precisionMode = 0;
%     % 0 (default) better looking plot, smaller files, faster calculations;
%     % set 1 for enhanced precision (smaller q and s steps)
% end

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

signal = y;

prof = cumsum(signal);
slength = size(prof);
fqs = [];

numIncrements = 20;
sListFull = unique(round(linspace(minScale,maxScale,numIncrements)));

timer = tic;
for s = sListFull

    if doOverlap
        vec = [0:s-1];
        ind = [1:slength-s+1]';
        coordinates = bsxfun(@plus, vec, ind);
    else
        ind2 = [1:size(prof,1)];
        coordinates = reshape(ind2(1:(size(prof,1)-mod(size(prof,1),s))),s,(size(prof,1)-mod(size(prof,1),s))/s)';
    end

    segments = prof(coordinates);
    xbase = [1:1:s];
    f2nis = [];

    for ni = 1:size(segments,1)
        seg = segments(ni,:);
        fit = polyfit(xbase,seg,2);
        variance = mean((seg - polyval(fit,xbase)).^2);
        f2nis(end+1) = variance;
    end

    for q = qList
        fqs = [fqs; q s (mean(f2nis.^(q/2)))^(1/q)];
    end
end
% fprintf(1,'Detrended fluctuations computed in %s\n',BF_TheTime(toc(timer)));

fqsll = [fqs(:,1) fqs(:,2) log(fqs(:,2)) log(fqs(:,3))];

% ------------------------------------------------------------------------------
% Now compute hurst exponents as the gradients of F(q) curves
% ------------------------------------------------------------------------------
hqs = [];

% if precisionMode == 0
%     % Take 11 points through the space
%     sspacing = ((maxScale/5)-minScale)/10;
%     sList = minScale:sspacing:(maxScale/5);
%     % Coarser sampling of q space
%     qList = qMin:1:qMax; qList(qList == 0) = 0.0001;
% else
    % sspacing = 1;
% sList = minScale:sspacing:(maxScale/5);

if sum(sListFull<=maxScale/5)>=10
    sList = sListFull(sListFull<=maxScale/5);
else
    % sample higher in the scale dimension:
    sspacing = ((maxScale/5)-minScale)/10;
    sList = minScale:sspacing:(maxScale/5);
end

% Coarser sampling of q space
qList = qMin:0.5:qMax; qList(qList == 0) = 0.0001;
% end

% sList = minScale:sspacing:(maxScale/5);

hqs = zeros(length(qList),length(sList));
for si = 1:length(sList)
    sit = sList(si);
    for qi = 1:length(qList)
        qit = qList(qi);

        fitTemp = fqsll(fqsll(:,1) == qit & fqsll(:,2) >= sit & fqsll(:,2) <= 5*sit,:);
        hTemp = polyfit(fitTemp(:,3),fitTemp(:,4),1);

        hqs(qi,si) = hTemp(1);
        % hqs = [hqs; qit 3*sit hTemp(1)];
    end
end

sListScaled = sList*3; % Not completely on top of the algorithm, but for some reason
                       % this was recorded as a multiple of 3 in the original algorithm

% hqsPlotData = reshape(hqs(:,3),size(qList,2),length(minScale:sspacing:(maxScale/5)));

% ------------------------------------------------------------------------------
% Plotting
% ------------------------------------------------------------------------------

if doPlot
    if max(max(hqs)) < 1.5
        hLim = 1.5;
    elseif max(max(hqs)) < 2.5
        hLim = 2.5;
    else
        hLim = ceil((max(max(hqs))*10))/10;
    end

    f = figure('color','w'); box('on');
    hqsplot = surf(sListScaled,qList,hqs);
    colormap(jet);
    colorbar;
    caxis([0,hLim]);
    set(gca,'YDir','reverse');
    view(-62,50);
    axis([sListScaled(1),sListScaled(end),qMin,qMax,0,hLim]);
    xlabel('scale')
    ylabel('q')
    zlabel('h')
end

% ------------------------------------------------------------------------------
% Output statistics:
% ------------------------------------------------------------------------------

% hqsPlotData (scale,q)

% Global properties:
allExponents = hqs(:);
out.meanHurstExponent = mean(allExponents);
out.stdHurstExponent = std(allExponents);
out.minHurstExponent = min(allExponents);
out.maxHurstExponent = max(allExponents);

% Changes with scale:
out.scaleHurstStd = std(mean(hqs,1));
out.scaleHurstTrend = GiveMeGradient(sListScaled,mean(hqs,1));

% Changes with q:
out.qHurstStd = std(mean(hqs,2));
out.qHurstTrend = GiveMeGradient(qList,mean(hqs,2));

% max/min points are where in scale/q space?
[qi,si] = find(hqs==max(hqs(:)),1);
out.maxHurstQ = qList(qi);
out.maxHurstScale = sListScaled(si);
[qi,si] = find(hqs==min(hqs(:)),1);
out.minHurstQ = qList(qi);
out.minHurstScale = sListScaled(si);

% phase transitions
% there is some peak or trough somewhere, so the standard deviation across
% scales or q is inconsistent
stdS = std(hqs,[],1);
stdQ = std(hqs,[],2);
out.stdStdHurstQ = std(stdQ); % will be large if variance changes alot with Q
out.stdStdHurstScale = std(stdS); % will be large if variance changes alot with scale

% saveas(hqsplot, [dirname filesep 'MMA_results' filesep 'MMA_' n1 '.jpg']);

% ------------------------------------------------------------------------------
function m = GiveMeGradient(xData,yData)
    if size(xData,1) ~= size(yData,1);
        yData = yData';
    end
    p = polyfit(xData,yData,1);
    m = p(1);
end
% ------------------------------------------------------------------------------

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function K = multiKernCompute(kern, varargin)
% MULTIKERNCOMPUTE Compute the MULTI kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the multiple output block
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the multiple output block
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : multiKernParamInit, kernCompute, kernCreate, multiKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Pei Gao, 2007

% KERN

if iscell(varargin{1})  % If x is a cell array of input points
  if length(varargin{1}) ~= kern.numBlocks
    error('Time information is not matched among blocks!');
  end
  for i = 1:kern.numBlocks
    dim1(i) = size(varargin{1}{i}, 1);
    if length(varargin)>1  % If x2 is provided it must also be a cell array
      if length(varargin{1}) ~= length(varargin{2})
        error('Time information is not matched within the block!');
      end
      dim2(i) = size(varargin{2}{i}, 1);
    else
      dim2(i) = dim1(i); % If x2 is provided is not provided just use x1
    end

  end

  K = zeros(sum(dim1), sum(dim2));

  for i = 1:kern.numBlocks
    % Compute the i:th diagonal block
    startOne = sum(dim1(1:(i-1))) + 1;
    endOne = sum(dim1(1:i));
    startThree = sum(dim2(1:(i-1))) + 1;
    endThree = sum(dim2(1:i));

    % Compute the block if it has a nonzero amount of rows and columns
    if ((endOne >= startOne) && (endThree >= startThree)),
      if length(varargin)<2
        K(startOne:endOne, startThree:endThree) = multiKernComputeBlock(kern, ...
                                                      varargin{1}{i}, i, i);
      else
        K(startOne:endOne, startThree:endThree) = multiKernComputeBlock(kern, ...
                                       varargin{1}{i}, varargin{2}{i}, i, i);
      end
    end;

    % Compute the off-diagonal blocks in the same row/column as the i:th diagonal block
    for j = 1:i-1
      if ~isempty(kern.block{i}.cross{j})
        startTwo = sum(dim2(1:(j-1))) + 1;
        endTwo =  sum(dim2(1:j));

        if length(varargin)<2
          % If a separate set of column-inputs has not been given,
          % row and column inputs are the same

	  %size(K)
	  %startOne
	  %endOne
	  %startTwo
	  %endTwo
	  %size(varargin{1})
	  %i
	  %j
	  %K(startOne:endOne, startTwo:endTwo)
	  %varargin{1}{i}
	  %varargin{1}{j}
          %kern
          %tempK = multiKernComputeBlock(kern, ...
          %                      varargin{1}{i}, varargin{1}{j}, i, j)

	  % Compute the block if it has a nonzero amount of rows and columns
	  if ((endOne >= startOne) && (endTwo >= startTwo)),
            K(startOne:endOne, startTwo:endTwo) = multiKernComputeBlock(kern, ...
                                  varargin{1}{i}, varargin{1}{j}, i, j);
  	    % Since the row inputs are the same as the column inputs, the
            % corresponding block on the other side of the diagonal is just
            % the transpose of the one just computed.
            K(startTwo:endTwo, startOne:endOne) = K(startOne:endOne, ...
                                                  startTwo:endTwo)';
          end;
        else
          % If a separate set of column-inputs has been given, use that

	  % Compute the block if it has a nonzero amount of rows and columns
	  if ((endOne >= startOne) && (endTwo >= startTwo)),
            K(startOne:endOne, startTwo:endTwo) = multiKernComputeBlock(kern, ...
                                  varargin{1}{i}, varargin{2}{j}, i, j);
          end;

	  % Since the row inputs are not the same as the column inputs, the
          % corresponding block on the other side of the diagonal is not just
          % a transpose, so we compute it separately here.
          startFour = sum(dim1(1:(j-1))) + 1;
          endFour =  sum(dim1(1:j));
	  % Compute the block if it has a nonzero amount of rows and columns
	  if ((endFour >= startFour) && (endThree >= startThree)),
            K(startFour:endFour, startThree:endThree) = ...
                multiKernComputeBlock(kern, varargin{2}{i}, varargin{1}{j}, j, i)';
          end;
        end
      end
    end
  end

else
  dim1 = size(varargin{1}, 1);

  if length(varargin)>1
    dim2 = size(varargin{2}, 1);
  else
    dim2 = dim1;
  end

  K = zeros(kern.numBlocks*dim1, kern.numBlocks*dim2);

  for i = 1:kern.numBlocks
    startOne = (i-1)*dim1 + 1;
    endOne = i*dim1;
    startThree = (i-1)*dim2 + 1;
    endThree = i*dim2;
    K(startOne:endOne, startThree:endThree) = multiKernComputeBlock(kern, ...
                                                      varargin{:}, i, i);
    for j = 1:i-1
      if ~isempty(kern.block{i}.cross{j})
        startTwo = (j-1)*dim2 + 1;
        endTwo = j*dim2;
        K(startOne:endOne, startTwo:endTwo) = multiKernComputeBlock(kern, ...
                                                          varargin{:}, i, j);
        if length(varargin)<2
          K(startTwo:endTwo, startOne:endOne) = K(startOne:endOne, ...
                                                startTwo:endTwo)';
        else
          startFour = (j-1)*dim1 + 1;
          endFour = j*dim1;
          K(startFour:endFour, startThree:endThree) = ...
              multiKernComputeBlock(kern, varargin{end:-1:1}, j, i)';
        end
      end
    end
  end

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function st = type( x, usegeo )

df  = x.dof_;
if nargin > 1 && usegeo,
    geo = any( df < 0 );
else
    geo = false;
end
df  = abs( df );
s   = x.size_;
len = prod( s );

if len == 1,
    isstruct = 0;
    nzs = len;
    st = 'scalar';
else
    isstruct = ~isempty( df ) & ( sum( df ) < ( 2 - isreal( x.basis_ ) ) * len );
    if ~isstruct,
        nzs = nnz( any( x.basis_, 1 ) );
    end
    nd = length( s );
    st = sprintf( '%dx', s );
    st = st( 1 : end - 1 );
    if nd > 2,
        st = [ st, ' array' ];
    elseif any( s == 1 ),
        st = [ st, ' vector' ];
    else
        st = [ st, ' matrix' ];
    end
end

if isstruct,
    st = sprintf( '%s, %d d.o.f.', st, df );
elseif nzs < len,
    if nzs > 1,
        st = sprintf( '%s, %d nonzeros', st, nzs );
    elseif nzs == 1,
        st = sprintf( '%s, 1 nonzero', st );
    end
end
if geo,
    st = [ st, ', geometric' ];
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ATilde,P,R,C]=equilibrateMat(A,algorithm,maxIter,RelTol,AbsTol)
%%EQUILIBRATEMAT Equilibrate a general (not necessarily symmetric) real or
%                complex matrix. This finds a permutation matrix and two
%                diagonal matrices such that the equilibrated matrix,
%                ATilde=P*R*A*C has the infinity-norm of each row and
%                column equal to 1.
%
%INPUTS: A An nXn real or complex matrix. Unless algorithm=1, in which case
%          A must be real and all positive.
% algorithm An optional parameter selecting the algorithm to use for
%          equilibration. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) Use the
%            algorithm involving bipartite matching, as described in
%            Section 2.1.2 of [1], which is one of the algorithmic variants
%            given in [2].
%          1 This is the heuristic ALG3 in [1]. The algorithm should work
%            with all-positive, real matrices. It does not work with
%            complex matrices. Also, the algorithm usually does not work
%            with matrices containing negative elements. This algorithm
%            always produces P=eye(n,n).
%  maxIter If algorithm=1, this is the maximum number of iterations to
%          perform. The default if omitted or an empty matrix is passed is
%          500. This is not used if algorithm=0.
% RelTol, AbsTol The absolute and relative convergence criteria for
%         algorithm 1. The defaults if omitted or empty matrices are passed
%         are 1e-10 and 1e-13.
%
%OUTPUTS: ATilde The nXn equlibrated matrix.
%              P The nNn permutation matrix.
%              R The left-multiplying diagonal matrix.
%              C The right-multiplying diagonal matrix.
%
%EXAMPLE:
%Here, we create a symmetric matrix and equilibrate it. The equilibration
%makes the largest magnitude values in each row and column equal 1.
%However, it also changes the condition number (usually for the better), as
%we shall show:
% %A random complex matrix.
% A=10*rand(10,10)+10*1j*rand(10,10);
% ATilde=equilibrateMat(A);
%
% %The original maximum absolute value elements in the columns and rows.
% max(ATilde,[],1)
% max(ATilde,[],2)
%
% %After equilibration, the maximum mangitude elements in the columns and
% % %rows are all 1.
% max(ATilde,[],1)
% max(ATilde,[],2)
%
% %The condition number often but not always improves.
% cond(A,Inf)
% cond(ATilde,Inf)
%
%REFERENCES:
%[1] P. Liu, "An exploration of matrix equilibration," University of
%    British Columbia. CPSC 517, Tech. Rep., 30 Sep. 2015. [Online].
%    Available: http://dx.doi.org/10.14288/1.0103601
%[2] I. S. Duff and J. Koster, "On algorithms for permuting large entries
%    to the diagonal of a sparse matrix," SIAM Journal on Matrix Analysis
%    and Applications, vol. 22, no. 4, pp. 973-996, 2001.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(AbsTol))
    AbsTol=1e-13;
end

if(nargin<4||isempty(RelTol))
    RelTol=1e-10;
end

if(nargin<3||isempty(maxIter))
    maxIter=500;
end

if(nargin<2||isempty(algorithm))
    algorithm=0;
end

switch(algorithm)
    case 0%The algorithm of [2].
        n=size(A,1);
        ATrans=-log(A);

        %In the current version of the Tracker Component Library, the compiled
        %version of assign2D does not support complex numbers in the cost matrix.
        %Thus, if the mex file is in the path, we have to remove it from the path
        %and use the m file. We can add the compiled version back to the path after
        %calling assign2D.
        [filepath,~,ext]=fileparts(which('assign2D'));

        if(~strcmp(ext,'m'))
            %If it isn't the m file that will be used, we have to remove the
            %filepath from the current path.
            rmpath(filepath);
        else
            filepath=[];
        end

        %Note that u and v are reversed; the definition in [1] is different from
        %that used in assign2D. ATrans will be complex if A had any negative terms
        %or if A were complex to begin with. However, the complex values are
        %correctly handled in the m function implementation of assign2D.
        [col4row,~,~,v,u]=assign2D(ATrans,false);

        %If a compiled implementation of assign2D had to be removed from the
        %current path.
        if(~isempty(filepath))
            addpath(filepath)
        end

        if(isempty(col4row))
            %The equilibration problem is infeasible.
            P=[];
            R=[];
            C=[];
            return
        end

        %We have to obtain the dual variables for the original problem, not the
        %transformed one that is sovled in assign2D.
        ADelta=min(ATrans(:));
        %If A is all positive, the problem wasn't transformed.
        if(ADelta>0)
            ADelta=0;
        end
        u=u+ADelta;

        %Construct the permutation matrix.
        P=zeros(n,n);
        for k=1:n
            P(k,col4row(k))=1;
        end

        R=diag(exp(u));
        C=diag(exp(v));

        %R and C should be real if A is real. This deals with any possible finite
        %precision issues.
        if(isreal(A))
            R=real(R);
            C=real(C);
        end

        ATilde=P*R*A*C;
    case 1%ALG3 in [1].
        if(~isreal(A))
            error('This algorithm only works with real matrices.')
        end

        if(any(A(:)<0))
            warning('This algorithm often does not satisfy the desired norm criterion when A contains negative elements.')
        end

        n=size(A,1);

        R=eye(n,n);
        C=eye(n,n);
        ATilde=A;

        D=eye(n,n);
        for curIter=1:maxIter
            ATildeOld=ATilde;
            for i=1:n
                alphaR=1/sqrt(norm(ATilde(i,:),Inf));
                alphaC=1/sqrt(norm(ATilde(:,i),Inf));

                Dr=D;
                Dr(i,i)=alphaR;
                Dc=D;
                Dc(i,i)=alphaC;

                R=R*Dr;
                C=C*Dc;
                ATilde=Dr*ATilde*Dc;
            end

            diffMag=abs(ATilde-ATildeOld);
            if(all((diffMag(:)<=RelTol*abs(ATilde(:)))|(diffMag(:)<=AbsTol)))
               return;
            end
        end

        P=eye(n,n);

    otherwise
        error('Unknown algorithm specified.')
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = ChristoffelTensor(S,varargin)
% Christoffel tensor of an elasticity tensor for a given direction
%
% Formula: E_jk = C_ijkl n_j n_l
%
% Input
%  S - elatic compliance @tensor
%  x - list of @vector3d
%
% Output
%  E - Christoffel @tensor
%
% See also
% tensor/directionalMagnitude tensor/rotate

% take formula using stiffness
[varargout{1:nargout}] = ChristoffelTensor(inv(S),varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function eeg_write_ascii(file,data,format)

% eeg_write_ascii - Write matrix of EEG data into ascii data format
%
% Useage: eeg_write_ascii(file,data,format)
%
% where:    file is the path + filename, eg "c:\data.dat"
%           data is a data matrix (floating point)
%           format is an fprintf format string ('%12.6f' default)
%
% comment:  Uses the fprintf method, which is quicker than the
%           dlmwrite() function. It remains to be seen whether it
%           is quicker or provides more format options than the
%           built-in SAVE function.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no express or implied warranties
% History:  07/2001, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

    % transpose data for output, see note below
    data = data';

    [r,c] = size(data);         % size of matrix: r rows and c columns

    % create a format string long enough to provide a format for every
    % row of the data - put a line return at the end of this string
    formats = '';
    for i = 1:r
        formats = strcat(formats, ' ', format);
    end
    formats = strcat(formats,'\n');

    % Output data
    fid = fopen(file,'w');  fprintf(fid,formats,data);  fclose(fid);

    fprintf('Completed ascii write to:\t%s\n', file);
    return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   NOTE on TRANSPOSE
%   this routine transposes in the input data matrix to
%   accommodate the matlab fprintf routine, which outputs in
%   column major format to effectively transpose the input.  For
%   example, 'help fprintf' provides this example:
%
%   x = 0:.1:1; y = [x; exp(x)];
%   fprintf('%6.2f  %12.8f\n',y);
%
%   where y is a 2row, 11 column matrix and the output is:
%
%   0.00    1.00000000
%   0.10    1.10517092
%        ...
%   1.00    2.71828183

author: PatternRecognition

repo: OpenBMI


text:
% RecoverByGraInitial.m
% RecoverByGra.m
% -------------------------------------------------------------------
% This function is just the joint of RecoverByGra and initial f0
% Date:    17/03/2013
% Last modified: 16/04/2015
% -------------------------------------------------------------------

function [imgRec, rms, gObj] = RecByGraInitial(img1, img2, ww1, ww2, dxdy, iter, res, alpha, iniMode)

    % ------------ Check parameter --------------
%     narginchk(9, 9);
%     if size(img1, 3) ~=1 || size(img2, 3) ~=1,
%         error('The image should be gray');
%     end
    % -------------------------------------------
    switch lower(iniMode),
        case 'avg',
            f0 = (img1+img2)/2;
        case 'weight',
            f0 = ww1.*img1+ww2.*img2;
        otherwise
            error('There only two mode');
    end

    [imgRec, rms, gObj] = RecoverByGra(f0, dxdy, iter, res, alpha);

end

%%
function [imgRec, rms, gObj] = RecoverByGra(imgOri, Obj, iter, res, alpha)

    lp = [0 1 0;1 -4 1;0 1 0];
    if isreal(Obj),
%         disp('The original is the Laplace');
        gObj = Obj;
    else
%         disp('The original is the Gradient');
        gObj = LaplaceZ(Obj);
    end

    rms = [];
%     f0 = Boundary(imgOri, dxObj, dyObj);
    f0 = imgOri;
    for ii = 1:iter,
%         f0 = Boundary(f0, dxObj, dyObj);
        deltaF = imfilter(int16(f0), lp, 'replicate', 'corr');
        deltaF = double(deltaF);
        delta = alpha * (deltaF-gObj);

        f1 = f0 + delta;

        f1 = max(f1, 0);
        f1 = min(f1, 255);

        f0 = f1;
    end

    disp(['The ' num2str(ii) ' iteration is complete.']);
    imgRec = f1;
end

%%
function lap = LaplaceZ(IMGORGRA)

    if isreal(IMGORGRA),
%         disp('The input should be IMAGE')
        lh=[0,  1, 0;...
            1, -4, 1;...
            0,  1, 0];
        lap = imfilter(IMGORGRA, lh, 'replicate', 'corr');
    else
%         disp('The input should be GRADIENT');
        fx=[0, 0, 0; 0, -1, 1; 0, 0, 0];
        fy=[0, 0, 0;  0, -1, 0; 0, 1, 0];

        dx = real(IMGORGRA);
        dy = imag(IMGORGRA);

        ddx = imfilter(dx, fx, 0, 'corr');
        ddy = imfilter(dy, fy, 0, 'corr');

        lap=ddx + ddy;
    end
end

author: thfylsty


text:
function [y]=fconv(x, h)
% Fast Parallelised Convolution
%   [y] = FCONV(x, h) convolves x and h in the frequency domain
%         to +-1.
%
%      x = input vector
%      h = input vector
%
%      See also CONV
%

Ly=size(x,1)+size(h,1)-1;  %
Ly2=pow2(nextpow2(Ly));    % Find smallest power of 2 that is > Ly

if isa(x, 'gpuArray')
Ly  = gpuArray(Ly);
Ly2 = gpuArray(Ly2);
end

X=fft(x, Ly2);             % Fast Fourier transform
H=fft(h, Ly2);	           % Fast Fourier transform
Y=X.*H;        	           %
y=real(ifft(Y, Ly2));      % Inverse fast Fourier transform
y=y(1:1:Ly,:);             % Take just the first N elements

author: jtkim-kaist

repo: Speech-enhancement


text:
function y = cvx_blkdiag( varargin )

% A reimplemenetation of MATLAB's blkdiag function, which is broken in
% certain versions of MATLAB.

if nargin == 1,
    y = varargin{1};
else
    isYsparse = false;
    for k = 1 : nargin,
        x = varargin{k};
        [ p2(k+1), m2(k+1) ] = size(x); %#ok
        if issparse(x), isYsparse = true; end
    end
    p1 = cumsum(p2);
    m1 = cumsum(m2);
    if isYsparse
        y = sparse( varargin{1} );
        for k = 2 : nargin,
            y = [y sparse(p1(k),m2(k+1)); sparse(p2(k+1),m1(k)) varargin{k}]; %#ok
        end
    else
        y = zeros(p1(end),m1(end));
        for k=1:nargin
            y(p1(k)+1:p1(k+1),m1(k)+1:m1(k+1)) = varargin{k};
        end
    end
end

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


text:
function convertNetwork(obj)
% convertNetwork(obj)
%
% Converts a network from test to train, inserting lossed, dropout and
% adopting the classification layer.
%
% Copyright by Holger Caesar, 2015

fprintf('Converting test network to train network (dropout, loss, etc.)...\n');

% Add dropout layers after relu6 and relu7
fprintf('Adding dropout to network...\n');
dropout6Layer = dagnn.DropOut();
dropout7Layer = dagnn.DropOut();
insertLayer(obj.net, 'relu6', 'fc7', 'dropout6', dropout6Layer);
insertLayer(obj.net, 'relu7', 'fc8', 'dropout7', dropout7Layer);

% Rename variable prob to x%d+ style to make insertLayer work
% (Required for beta18 or later matconvnet default networks)
predVarIdx = obj.net.getVarIndex('prediction');
if ~isnan(predVarIdx)
    freeVarName = getFreeVariable(obj.net);
    obj.net.renameVar('prediction', freeVarName);
end

% Replace softmax with correct loss for training (default: softmax)
switch obj.nnOpts.lossFnObjective
    case 'softmaxlog'
        softmaxlossBlock = dagnn.LossWeighted('loss', 'softmaxlog');
        replaceLayer(obj.net, 'prob', 'softmaxloss', softmaxlossBlock, 'label');
        objectiveName = obj.net.layers(obj.net.getLayerIndex('softmaxloss')).outputs;
        obj.net.renameVar(objectiveName, 'objective');
    case 'hinge'
        hingeLossBlock = dagnn.Loss('loss', 'hinge');
        replaceLayer(obj.net, 'prob', 'hingeloss', hingeLossBlock, 'label');
        objectiveName = obj.net.layers(obj.net.getLayerIndex('hingeloss')).outputs;
        obj.net.renameVar(objectiveName, 'objective');
    otherwise
        error('Unknown loss specified');
end

% Adapt number of classes in softmaxloss layer from 1000 to numClasses
fc8Idx = obj.net.getLayerIndex('fc8');
obj.net.layers(fc8Idx).block.size(4) = obj.imdb.numClasses;
newParams = obj.net.layers(fc8Idx).block.initParams();
obj.net.params(obj.net.layers(fc8Idx).paramIndexes(1)).value = newParams{1} / std(newParams{1}(:)) * 0.01; % Girshick initialization
obj.net.params(obj.net.layers(fc8Idx).paramIndexes(2)).value = newParams{2}';

% Rename input
if ~isnan(obj.net.getVarIndex('x0'))
    % Input variable x0 is renamed to input
    obj.net.renameVar('x0', 'input');
else
    % Input variable already has the correct name
    assert(~isnan(obj.net.getVarIndex('input')));
end

% Modify for Fast Rcnn (ROI pooling, bbox regression etc.)
if obj.nnOpts.fastRcnn
    obj.convertNetworkToFastRcnn();
end

author: nightrome

repo: matconvnet-calvin


text:
function rhs=ch_rom_deim_rhs(tspan, a,dummy,P_NL,P_Psi,L)
N=P_Psi*a;
rhs=L*a + i*P_NL*((abs(N).^2).*N);

author: dynamicslab


text:
function [cpulse, verbose] = tapas_physio_get_cardiac_pulses_manual_template(...
    c, t, pulse_detect_options, verbose)
% Detects R-peaks via matched-filter smoothing & peak detection using a
% manually defined QRS-wave (or R-peak environment)
%
%   [cpulse, verbose] = tapas_physio_get_cardiac_pulses_manual_template(...
%    c, t, thresh_min, dt120, verbose)
%
% IN
%   c               [nSamples, 1] raw pulse oximeter samples
%   t               [nSamples, 1] time vector corresponding to samples (un seconds)
%   pulse_detect_options
%                   physio.thresh.cardiac.initial_cpulse_select-substructure
%                   with elements
%                   .method 'manual' or 'load'/'load_template'
%                           'manual' - select template manually
%                           'load'/'load_template' - load from template
%                   .min     threshold for correlation with QRS-wave to find cardiac pulses
%                   .file   variable saving an example cardiac QRS-wave to correlate with ECG time series
%   verbose         Substructure of Physio, holding verbose.level and
%                   verbose.fig_handles with plotted figure handles
%                   debugging plots for thresholding are only provided, if verbose.level >=2
%
% OUT
%
% EXAMPLE
%   tapas_physio_get_cardiac_pulses_manual_template
%
%   See also tapa_physio_new

% Author: Lars Kasper
% Created: 2012-02-20
% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the physIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

if nargin < 5
    verbose.level = 0;
    verbose.fig_handles = [];
end

% manual peak selection, if no file selected and loading is
% specified

hasKrpeakLogfile = exist(pulse_detect_options.file,'file') || ...
    exist([pulse_detect_options.file '.mat'],'file');

% if no file exists, also do manual peak-find
doSelectTemplateManually = any(strcmpi(...
    pulse_detect_options.method, ...
    {'manual', 'manual_template'})) || ~hasKrpeakLogfile;

if doSelectTemplateManually
    pulse_detect_options.kRpeak = [];
    hasECGMin = isfield(pulse_detect_options, 'min') && ~isempty(pulse_detect_options.min);
    if ~hasECGMin
        pulse_detect_options.min = 0.5;
    end
else
    fprintf('Loading %s\n', pulse_detect_options.file);
    ECGfile = load(pulse_detect_options.file);
    pulse_detect_options.min = ECGfile.ECG_min;
    pulse_detect_options.kRpeak = ECGfile.kRpeak;
end

inp_events = [];
ECG_min = pulse_detect_options.min;
kRpeak = pulse_detect_options.kRpeak;
if doSelectTemplateManually
    while ECG_min
        [cpulse, ECG_min_new, kRpeak] = tapas_physio_find_ecg_r_peaks(t,c, ECG_min, [], inp_events);
        fprintf('Press 0, then return, if right ECG peaks were found\n');
        ECG_min = input('otherwise type next numerical choice for ECG_min and continue the selection: ');
    end
else
    [cpulse, ECG_min_new, kRpeak] = tapas_physio_find_ecg_r_peaks(t,c, ECG_min, kRpeak, inp_events);
end
ECG_min = ECG_min_new;
cpulse = t(cpulse);

% save manually found peak parameters to file
if doSelectTemplateManually
    save(pulse_detect_options.file, 'ECG_min', 'kRpeak');
end

author: translationalneuromodeling

repo: tapas


text:
function [] = computeRerankScores(modelFiles, outputFile,varargin)
% Test a trained LSTM model by generating translations.
% Arguments:
%   modelFiles: single or multiple models to decode. Multiple models are
%     separated by commas.
%   batchSize: number of sentences decoded simultaneously. We only ensure
%     accuracy of batchSize = 1 for now.
%   outputFile: output translation file.
%   varargin: other optional arguments.
%
% Thang Luong @ 2016, <lmthang@stanford.edu>

  addpath(genpath(sprintf('%s/..', pwd)));

  %% Argument Parser
  p = inputParser;
  % required
  addRequired(p,'modelFiles',@ischar);
  addRequired(p,'outputFile',@ischar);

  % optional
  addOptional(p,'gpuDevice', 0, @isnumeric); % choose the gpuDevice to use: 0 -- no GPU
  addOptional(p,'align', 0, @isnumeric); % 1 -- output aignment from attention model
  addOptional(p,'assert', 0, @isnumeric); % 1 -- assert
  addOptional(p,'debug', 0, @isnumeric); % 1 -- debug
  addOptional(p,'testPrefix', '', @ischar); % to specify a different file for decoding
  addOptional(p,'continueId', 0, @isnumeric); % > 0: start decoding from this continueId (base 1) sent and append the results
  addOptional(p,'computePpl', 0, @isnumeric); % 1 -- compute ppl of each participating model over the test set

  % useful for rescoring if we have many sentence pairs with the same source
  addOptional(p,'reuseEncoder', 0, @isnumeric);
  % print decoding scores, require stackSize = 1
  addOptional(p,'printScore', 1, @isnumeric);

  p.KeepUnmatched = true;
  parse(p,modelFiles,outputFile,varargin{:})
  decodeParams = p.Results;

  % GPU settings
  decodeParams.isGPU = 0;
  if decodeParams.gpuDevice
    n = gpuDeviceCount;
    if n>0 % GPU exists
      fprintf(2, '# %d GPUs exist. So, we will use GPUs.\n', n);
      decodeParams.isGPU = 1;
      gpuDevice(decodeParams.gpuDevice)
      decodeParams.dataType = 'single';
    else
      decodeParams.dataType = 'double';
    end
  else
    decodeParams.dataType = 'double';
  end
  printParams(2, decodeParams);

  %% load multiple models
  tokens = strsplit(decodeParams.modelFiles, ',');
  numModels = length(tokens);
  models = cell(numModels, 1);
  for mm=1:numModels
    [models{mm}] = loadDecodeModel(tokens{mm}, decodeParams);

    % make sure all models have the same vocab, and the number of layers
    if mm>1
      for ii=1:models{mm}.params.srcVocabSize
        assert(strcmp(models{mm}.params.srcVocab{ii}, models{1}.params.srcVocab{ii}), '! model %d, mismatch src word %d: %s vs. %s\n', mm, ii, models{mm}.params.srcVocab{ii}, models{1}.params.srcVocab{ii});
      end
      for ii=1:models{mm}.params.tgtVocabSize
        assert(strcmp(models{mm}.params.tgtVocab{ii}, models{1}.params.tgtVocab{ii}), '! model %d, mismatch tgt word %d: %s vs. %s\n', mm, ii, models{mm}.params.tgtVocab{ii}, models{1}.params.tgtVocab{ii});
      end
      models{mm}.params = rmfield(models{mm}.params, {'srcVocab', 'tgtVocab'});
    end
  end
  params = models{1}.params;

  if params.continueId > 0 % appending
    fileOpt = 'a';
  else
    fileOpt = 'w';
  end

  params.fid = fopen(params.outputFile, fileOpt);
  params.logId = fopen([outputFile '.log'], fileOpt);
  % align
  if params.align
    params.alignId = fopen([params.outputFile '.align'], fileOpt);
  end
  % print score
  if params.printScore
    params.scoreFid = fopen([params.outputFile '.score'], fileOpt);
  end

  printParams(2, params);

  % default params
  params.hasTgt = 1;
  batchSize = 1;
  params.prefixDecoder = 0;

  % load test data
  [srcSents, tgtSents, numSents]  = loadBiData(params, params.testPrefix, params.srcVocab, params.tgtVocab, -1, params.hasTgt);

  % compute perplexities
  if decodeParams.computePpl
    [testData] = prepareData(srcSents, tgtSents, 1, params);
    for mm=1:numModels
      evalValidTestSimple(models{mm}, testData, models{mm}.params);
    end
  end

  %%%%%%%%%%%%
  %% decode %%
  %%%%%%%%%%%%
  numBatches = floor((numSents-1)/batchSize) + 1;

  fprintf(2, '# Decoding %d sents, %s\n', numSents, datestr(now));
  fprintf(params.logId, '# Decoding %d sents, %s\n', numSents, datestr(now));
  startTime = clock;

  assert(length(models) == 1);
  totalScore = 0;
  totalPredict = 0;
  for batchId = 1 : numBatches
    % prepare batch data
    startId = (batchId-1)*batchSize+1;
    endId = batchId*batchSize;
    if endId > numSents
      endId = numSents;
    end

    % continue training
    if params.continueId > startId
      continue;
    end

    % prepare data
    [decodeData] = prepareData(srcSents(startId:endId), tgtSents(startId:endId), 1, params);
    decodeData.startId = startId;

    testCosts = lstmCostGrad(models{1}, decodeData, params, 1);
    totalScore = totalScore + testCosts.word;
    totalPredict = totalPredict + length(testCosts.indLosses);

    if params.printScore
      printSent(2, decodeData.srcInput, params.srcVocab, ['  src ' startId ': ']);
      printSent(2, decodeData.tgtOutput, params.tgtVocab, ['  tgt ' startId ': ']);
      fprintf(2, 'score %g, ind scores [%s]\n', -testCosts.word, num2str(-testCosts.indLosses'));
      fprintf(params.scoreFid, '%g\n', params.scoreFid);
    end

    if mod(batchId, 100) == 0
      fprintf(2, '# Batch %d\n', batchId);
    end
    % printDecodeResults(decodeData, {tgtSents(startId)}, -testCosts.word, [], params, 1, []);
  end

  endTime = clock;
  timeElapsed = etime(endTime, startTime);
  fprintf(2, '# Complete reranking %d sents, num predict %d, ppl %g, time %.0fs, %s\n', numSents, totalPredict, ...
    exp(totalScore/totalPredict), timeElapsed, datestr(now));
  fprintf(params.logId, '# Complete reranking %d sents, num predict %d, ppl %g, time %.0fs, %s\n', numSents, totalPredict, ...
    exp(totalScore/totalPredict), timeElapsed, datestr(now));

  fclose(params.fid);
  fclose(params.logId);

  if params.printScore
    fclose(params.scoreFid);
  end
end

function [perp] = evalValidTestSimple(model, testData, params)
  fprintf(2, '  evaluating ...');
%   perp = 0;
  [testCosts] = evalCost(model, testData, params); % run on the test data
  testCounts = initCosts();
  testCounts = updateCounts(testCounts, testData);
  testCosts = scaleCosts(testCosts, testCounts);
  perp = exp(testCosts.word);
  fprintf(2, '  perp=%f\n', perp);
end

%     % convert local paths to absolute paths
%     fieldNames = fields(models{mm}.params);
%     for ii=1:length(fieldNames)
%       field = fieldNames{ii};
%       if ischar(models{mm}.params.(field))
%         if strfind(models{mm}.params.(field), '~lmthang/') ==1
%           models{mm}.params.(field) = strrep(models{mm}.params.(field), '~lmthang/', '/afs/ir/users/l/m/lmthang/');
%         end
%         if strfind(models{mm}.params.(field), '~lmthang/') ==1
%           models{mm}.params.(field) = strrep(models{mm}.params.(field), '~lmthang/', '/afs/cs.stanford.edu/u/lmthang/');
%         end
%         if strfind(models{mm}.params.(field), '~lmthang/') ==1
%           models{mm}.params.(field) = strrep(models{mm}.params.(field), '~lmthang/', '/home/lmthang/');
%         end
%       end
%     end

author: lmthang


text:
function [p] = elec_write_emse(p)

% elec_write_emse - Write an EMSE probe file (*.elp)
%
% Usage: [p] = elec_write_emse(p)
%
% This script outputs x,y,z values to an EMSE
% probe (*.elp) file.
%
% EMSE *.elp files are in meters.  Also, when using
% EMSE *.elp files in the eeg_toolbox, it is necessary
% to swap X and Y.
%
% See also: emse_read_elp, elec_open
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no express or implied warranties
% History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[path,name,ext] = fileparts(strcat(p.elec.path,filesep,p.elec.file));
ext = '.elp';
file = fullfile(path,[name ext]);

p.elec.file = [name ext];

[fid,msg] = fopen(file,'w','ieee-le');
if ~isempty(msg), error(msg); end

tic

fprintf('\nELEC_WRITE_EMSE...\n');
fprintf('...writing .elp data to:\n\t%s\n',file);

write_elp(fid,p);

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [elp] = write_elp(fid,p)

    Nelec = size(p.elec.data.x,1);

    % Probe files contain position information for electrode locations
    % and/or gradiometer locations. The file consists of a prolog, a
    % header, and a list of one or more sensor fields.

    % Any line beginning with '//' is a comment line, which is ignored

    % Write the prolog
    fprintf(fid,'%d %d\n%d\n',3,2,1);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Write the header
    % The header consists of one optional entry and 2 entries in
    % mandatory sequence and one optional entry:
    % Name [optional] > %N %s replace %s with name string (8 or fewer characters)
    % Type Code       > %x    replace %x with 1 (all electric), 2 (all magnetic) or 4 (mixed).
    % #Channels       > %d    number of points per channel per epoch [???? DLW]

    fprintf(fid,'//TypeCode nsensors\n');
    fprintf(fid,'%d %d\n',1,Nelec + 1);     % Nelectrodes plus REF

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Fiducial points may be included optionally. They are required
    % for MRI registration. If they are included, they must be in
    % the obligatory order : nasion, left preauricular point,
    % right preauricular point. Table A-2 defines the format for
    % representing fiduciary points.

    fprintf(fid,'//Fiducials (Nasion, Left Preauricular, Right Preauricular)\n');
    fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','%F',...
        p.elec.data.nasion(1),...
        p.elec.data.nasion(2),...
        p.elec.data.nasion(3));
    fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','%F',...
        p.elec.data.lpa(1),...
        p.elec.data.lpa(2),...
        p.elec.data.lpa(3));
    fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','%F',...
        p.elec.data.rpa(1),...
        p.elec.data.rpa(2),...
        p.elec.data.rpa(3));

    % Each electrode is represented by an electric sensor,
    % and consists of 5 fields, of which 1 (the name) is
    % optional.
    % Name              Format      Description
    % Type Code         %S          %x replace %x with 400 (electrode) or 1c00 if reference channel
    % Name [optional]   %N          %s replace %s with name string (8 or fewer characters)
    % Position          %g %g %g    electrode location with respect to head frame (Cartesian, meters)
    % Orientation       %g %g %g    not used, replace with 0 0 1
    %
    % Sensor state (which appears in the 'type code' field) may
    % be obtained by logically OR-ing suitable combinations from
    % the table below. Note that not all combinations are physically valid.
    %
    % type/state        type code
    % magnetic          200
    % electric          400
    % off               800
    % reference         1000        actually '1c00' [DLW]
    % optical           4000
    % trigger           8000
    % other             10000
    % named point       20000
    %
    % Other types (such as named points, trigger, and optical) should
    % be represented in the same pattern as electrodes, with the type
    % code set to identify the type. Even those types (e.g. trigger)
    % which do not have a true location, should have a nominal
    % location, (e.g. 0 0 0).

    % Note below that eeg_toolbox x,y,z are equivalent to
    % emse y,-x,z

    typecode = 400;

    for n = 1:Nelec,
        fprintf(fid,'\n//ecSensor typecode/state---------------------\n');
        fprintf(fid,'%s\t%d\n','%S',typecode);
        fprintf(fid,'//ecSensor name:\n');
        fprintf(fid,'%s\t%s\n','%N',char(p.elec.data.label(n)));
        fprintf(fid,'//sphere origin\n');
        fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','%O',...
            p.elec.data.centroid(2),...
            p.elec.data.centroid(1) * -1,...
            p.elec.data.centroid(3));
        fprintf(fid,'//ecSensor location (origin)\n');
        fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','  ',...
            p.elec.data.y(n),...
            p.elec.data.x(n) * -1,...
            p.elec.data.z(n));
    end

    % Output the REF
    typecode = '1c00'; % This is not consistent with the
                       % table above, but it works!
    fprintf(fid,'\n//ecSensor typecode/state---------------------\n');
    fprintf(fid,'%s\t%s\n','%S',typecode);
    fprintf(fid,'//ecSensor name:\n');
    fprintf(fid,'%s\t%s\n','%N','Ref');
    fprintf(fid,'//sphere origin\n');
    fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','%O',...
        p.elec.data.centroid(2),...
        p.elec.data.centroid(1) * -1,...
        p.elec.data.centroid(3));
    fprintf(fid,'//ecSensor location (origin)\n');
    fprintf(fid,'%s\t%12.8f\t%12.8f\t%12.8f\n','  ',...
        p.elec.data.ref(2),...
        p.elec.data.ref(1) * -1,...
        p.elec.data.ref(3));

    fclose(fid);

return


text:
% Laser to Camera calibration parameters (II optim stages)
% 22-Jan-2017 22:26:16
%
% Transformation matrix specifies laser coordinate frame
% in the reference frame of the camera
%
%-- Translation vector (t)
t = [ -0.193510 ; 0.031471 ; -0.094037 ]
%-- Rotation matrix (R)
R = ...
[ -0.523357  0.852009  -0.013370 ;...
  0.117409  0.056561  -0.991472 ;...
  -0.843986  -0.520463  -0.129635 ]

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
function [ a, b ] = p01_lim ( dim_num )

%*****************************************************************************80
%
%% P01_LIM returns the integration limits for problem 01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
function [tt] = tt_qtoepl(x, d)

% returns the multilevel Toeplitz matrix tt generated by the multi-dimensional input vector x
% in the QTT format
%
% The size of the input vector is 2^(d(1) + 1) x ... x 2^(d(D) + 1),
% The output matrix is square of order 2^d(1) x ... x 2^d(D).
% The (i_1,...i_D)-th component of x is not used if at least one of i_1,...,i_D is equal to 1
% (e. g., in the one-dimensional case the first component of x is not used)
%
% April 20, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For details please see the preprint
% http://www.mis.mpg.de/publications/preprints/2011/prepr2011-36.html
% Vladimir A. Kazeev, Boris N. Khoromskij and Eugene E. Tyrtyshnikov
% Multilevel Toeplitz matrices generated by QTT tensor-structured vectors and convolution with logarithmic complexity
% January 12, 2012
% Vladimir Kazeev,
% Seminar for Applied Mathematics, ETH Zurich
% vladimir.kazeev@sam.math.ethz.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D = numel(d);
tt = tt_qshiftstack(d);

sz=[];
for K=1:D
	sz = [sz; 4*ones(d(K),1), 2*ones(d(K),1); 1, 2];
end

tt = tt_qreshape(tt, 3, sz);
tt = tt_matrix(tt) * x;

tt1 = cell(sum(d),1);
ind = 0;
ind1 = 0;
for K = 1:D
	for k = 1:d(K)-1
		tt1{ind1 + k} = tt{ind + k};
	end

	p = size(tt{ind + d(K)}, 1);
	q = size(tt{ind + d(K)}, 3);
	r = size(tt{ind + d(K) + 1}, 3);
	tt1{ind1 + d(K)} = reshape(reshape(permute(tt{ind+d(K)},[2, 1, 3]),[4*p,q])*...
        reshape(permute(tt{ind+d(K)+1},[2, 1, 3]),[q,r]),[4,p,r]);
    tt1{ind1 + d(K)} = permute(tt1{ind1 + d(K)}, [2, 1, 3]);
	ind=ind+d(K)+1;
	ind1=ind1+d(K);
end
tt = tt1;
tt1 = tt_tensor;
tt = cell2core(tt1, tt);
tt = tt_matrix(tt);
return
end

author: oseledets

repo: TT-Toolbox


text:
function [dimLabels, resolutions, nSamples, units, firstSamplingPoint] = ...
    read_nifti(~, fileName)
% Reads nifti files and extracts properties for dimInfo
%
%   Y = MrDimInfo()
%   [dimLabels, resolutions, nSamples, units] = Y.read_nifti(fileName)
%   used in Y.load(fileName)
%
% This is a method of class MrDimInfo.
%
% IN
%
% OUT
%
% EXAMPLE
%   read_nifti
%
%   See also MrDimInfo

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2017-11-02
% Copyright (C) 2017 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

%% read header info
V = spm_vol(fileName);

%% nSamples
% for nifti files, nSamples should always has 3 entries, even if just a
% single line or slice is in the data (the nifti header will have
% information on the remaining singleton dimensions)
nSamples = [1, 1, 1];

% fill with info from file
tempSamples = V(1).private.dat.dim;
nSamples(1:numel(tempSamples)) = tempSamples;
tempDimInfo = MrDimInfo('nSamples', nSamples);
nDims = numel(nSamples);

%% dimLabels
dimLabels = tempDimInfo.dimLabels;

%% units
units = tempDimInfo.units;

%% resolutions
N = floor(abs(log10(eps('double'))));
P = round(tapas_uniqc_spm_imatrix(V(1).mat),N);
resolution_mm  = P(7:9);
% some nifti formats supply timing information (for files with more than 3
% dimension)
if isfield(V(1), 'private') && nDims > 3
    if isstruct(V(1).private.timing)
        TR_s = V(1).private.timing.tspace;
        tStart = 0;
    else
        TR_s = 1;
        units{4} = 'samples';
        tStart = 1;
    end
end
% need nifti to reference first sampling point as offcenter
resolutions = ones(1, nDims);
resolutions(1, 1:3) = resolution_mm;
if nDims > 3
    resolutions(1,4) = TR_s;
end

%% firstSamplingPoint
% voxel position by voxel center
% first sampling points for x,y,z are -FOV/2+res/2 such that [0 0 0] is in
% the centre of the matrix
% time starts at 0 seconds or 1 (1st) sample
FOV = nSamples(1:3) .* resolutions(1:3);
firstSamplingPoint = zeros(1, nDims);
firstSamplingPoint(1:3) = -FOV/2 + resolutions(1:3)/2;
if nDims > 3
    firstSamplingPoint(4) = tStart;
end
end


text:
function [s] = newid()
%returns a unique identifier string
%generated by AFNI's 3dnewid
   [ee,s] = system('3dnewid -fun');
   s = s(1:length(s)-1);
   return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function X = ncvar(X)
X = ncvar(struct(X));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function s = addrice(s)
%ADDRICE Add the Rician distribution.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.10 $  $Date: 2004/02/01 22:10:34 $

j = length(s) + 1;
s(j).name = 'Rician';
s(j).code = 'rician';
s(j).pnames = {'s' 'sigma'};
s(j).pdescription = {'noncentrality' 'scale'};
s(j).prequired = [false false];
s(j).fitfunc = @ricefit;
s(j).likefunc = @ricelike;
s(j).cdffunc = @ricecdf;
s(j).pdffunc = @ricepdf;
s(j).invfunc = @riceinv;
s(j).statfunc = @ricestat;
s(j).loginvfunc = [];
s(j).logcdffunc = [];
s(j).hasconfbounds = false;
s(j).censoring = true;
s(j).paramvec = true;
s(j).support = [0 Inf];
s(j).closedbound = [false false];
s(j).iscontinuous = true;
s(j).islocscale = false;
s(j).uselogpp = false;

% ==== Rician distribution functions ====

% these distribution functions do not yet handle arrays of parameters

function y = ricepdf(x,s,sigma)
%RICEPDF Rician probability density function (pdf).
s(s < 0) = NaN;
sigma(sigma <= 0) = NaN;

x(x<0) = 0;
sigsq = sigma.^2;
expon = (x.^2 + s.^2)./(2.*sigsq);
y = (x./sigsq) .* exp(-expon) .* besseli(0, x.*s./sigsq);
y(expon > (log(realmax(class(x)))-1)) = 0; % fix up 0*Inf

function p = ricecdf(x,s,sigma)
%RICECDF Rician cumulative distribution function (cdf).
s(s < 0) = NaN;
sigma(sigma <= 0) = NaN;

x(x<0) = 0;
p = ncx2cdf((x./sigma).^2, 2, (s./sigma).^2);

function x = riceinv(p,s,sigma)
%RICEINV Inverse of the Rician cumulative distribution function (cdf).
s(s < 0) = NaN;
sigma(sigma <= 0) = NaN;

x = sigma .* sqrt(ncx2inv(p, 2, (s./sigma).^2));

function r = ricernd(s,sigma,varargin)
%RICERND Random arrays from the Rician distribution.
s(s < 0) = NaN;
sigma(sigma <= 0) = NaN;

[err, sizeOut] = statsizechk(2,s,sigma,varargin{:});
if err > 0
    error('stats:ricernd:InconsistentSizes','Size information is inconsistent.');
end

r = sigma .* sqrt(ncx2rnd(2, (s./sigma).^2, sizeOut));

function [m,v] = ricestat(s,sigma)
%RICESTAT Mean and variance for the Rician distribution.
s(s < 0) = NaN;
sigma(sigma <= 0) = NaN;

t = .5 .* (s./sigma).^2;
m = sigma.*sqrt(.5.*pi).*exp(-.5.*t) .* ((1+t).*besseli(0,.5.*t) + t.*besseli(1,.5.*t));
v = 2.*sigma.^2 + s.^2 - m.^2;

function [nlogL,acov] = ricelike(params,data,cens,freq)
%RICELIKE Negative log-likelihood for the Rician distribution.
if nargin < 4 || isempty(freq), freq = ones(size(data)); end
if nargin < 3 || isempty(cens), cens = zeros(size(data)); end

nlogL = rice_nloglf(params, data, cens, freq);
if nargout > 1
    acov = mlecov(params, data, 'nloglf',@rice_nloglf, 'cens',cens, 'freq',freq);
end

% ==== Rician fitting functions ====

function [phat,pci] = ricefit(x,alpha,cens,freq,opts)
%NAKAFIT Parameter estimates and confidence intervals for Rician data.

if nargin < 2 || isempty(alpha), alpha = .05; end
if nargin < 3 || isempty(cens), cens = zeros(size(x)); end
if nargin < 4 || isempty(freq), freq = ones(size(x)); end
if nargin < 5, opts = []; end

if any(x <= 0)
    error('stats:ricefit:BadData','The data in X must be positive');
end

% Moment estimators of the uncensored data as starting point
% E[x.^2] = s.^2 + 2.*sigma.^2
% E[x.^4] = s.^4 + 8.*s.^2.*sigma.^2 + 8.*sigma.^4
xsqunc = x(cens == 0).^2;
meanxsq = mean(xsqunc); meanx4th = mean(xsqunc);
if meanxsq.^2 < meanx4th && meanx4th < 2.*meanxsq.^2
    s4th = 2.*meanxsq.^2 - meanx4th;
    ssq = sqrt(s4th);
    sigsq = .5.*(meanxsq - ssq);
    start = [sqrt(ssq) sqrt(sigsq)];
else
    start = cast([1 1],class(x));
end

% The default options include turning fminsearch's display off.  This
% function gives its own warning/error messages, and the caller can turn
% display on to get the text output from fminsearch if desired.
options = statset(statset('ricefit'), opts);
tolBnd = options.TolBnd;
options = optimset(options);

% Maximize the log-likelihood with respect to s and sigma.
[phat,nll,err,output] = ...
    fminsearch(@rice_nloglf, start, options, x, cens, freq, tolBnd);
if (err == 0)
    % fminsearch may print its own output text; in any case give something
    % more statistical here, controllable via warning IDs.
    if output.funcCount >= options.MaxFunEvals
        wmsg = 'Maximum likelihood estimation did not converge.  Function evaluation limit exceeded.';
    else
        wmsg = 'Maximum likelihood estimation did not converge.  Iteration limit exceeded.';
    end
    warning('stats:ricefit:IterOrEvalLimit',wmsg);
elseif (err < 0)
    error('stats:ricefit:NoSolution',...
          'Unable to reach a maximum likelihood solution.');
end

% Compute CIs using a normal approximation for phat.
if nargout > 1
    acov = mlecov(phat, x, 'nloglf',@rice_nloglf, 'cens',cens, 'freq',freq);
    probs = [alpha/2; 1-alpha/2];
    se = sqrt(diag(acov))';
    pci = norminv([probs probs], [phat; phat], [se; se]);
end

function nll = rice_nloglf(parms, x, cens, freq, tolBnd)
%RICE_NLOGLF Objective function for Rician maximum likelihood.
s = parms(1);
sigma = parms(2);
sigsq = sigma.^2;

% Restrict sigma to the open interval (0, Inf).
if nargin > 4
    if s < tolBnd || sigma < tolBnd
        nll = Inf;
        return
    end
end

bess0 = besseli(0, x.*s./sigsq);
rsq = (x.^2 + s.^2)./(2.*sigsq);
L = -rsq + log(bess0) + log(x./sigsq);
ncen = sum(freq.*cens);
if ncen > 0
    cen = (cens == 1);
    xcen = x(cen);
    L(cen) = log(marcumq(s./sigma,xcen./sigma));
end
nll = -sum(freq .* L);

% Don't have derivatives of the Marcum's Q, so can't compute an analytic
% gradient with censoring.
%
% if nargout > 1
%     dlogbess0 = besseli(1, x.*s./sigsq) ./ bess0;
%     dL1 = (-s + dlogbess0.*x) ./ sigsq;
%     dL2 = (rsq - 1 - dlogbess0.*x.*s./sigsq) ./ sigma;
%     if ncen > 0
% %         dL1(cen) = ;
% %         dL2(cen) = ;
%     end
%     ngrad = -[sum(freq .* dL1) sum(freq .* dL2)];
% end

function Q = marcumq(a,b)
% Q = MARCUMQ(A,B) returns Marcum's "Q" function.

if isa(a,'single') || isa(b,'single')
   Q = repmat(single(NaN), size(a));
else
   Q = repmat(NaN, size(a));
end
Q(a~=Inf & b==0) = 1;
Q(a~=Inf & b==Inf) = 0;
Q(a==Inf & b~=Inf) = 1;
z = (isnan(Q) & a==0 & b~=Inf);
if (any(z))
   Q(z) = exp((-b(z).^2)./2);
end

z = isnan(Q) & ~isnan(a) & ~isnan(b);
if (any(z(:)))
%    aa = (a(z).^2)./2;
   aa = (a.^2)./2;
   bb = (b(z).^2)./2;

   d = exp(-aa);
   h = d;
   f = bb.*exp(-bb);
   k = 1;
   delta = f .* h;
   sum = delta;
   j = (delta > sum.*eps(class(delta)));
   while any(j)
      d = aa.*d./k;
      h = h + d;
      f = bb.*f./(k+1);
      delta = f .* h;
      sum(j) = sum(j) + delta(j);
      j = (delta > sum.*eps(class(delta)));
      k = k + 1;
   end
   Q(z) = 1 - sum;
end

author: zouchuhang

repo: LayoutNet


text:
% This is used to set the evaluation metrics that to be computed
%
% Author:Xingchen Zhang, Ping Ye, Gang Xiao
% Contact: xingchen.zhang@imperial.ac.uk

function metrics=configMetrics

 metricVIFB={struct('name','Avg_gradient'),...
     struct('name','Cross_entropy'),...
     struct('name','Edge_intensity'),...
     struct('name','Entropy'),...
     struct('name','Mutinf'),...
     struct('name','Psnr'),...
     struct('name','Qabf'),...
     struct('name','Qcb'),...
     struct('name','Qcv'),...
     struct('name','Rmse'),...
     struct('name','Spatial_frequency'),...
     struct('name','Ssim'),...
     struct('name','Variance'),...
     };

 metrics = [metricVIFB];

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [alpha,b,eta,zeta,kappa] = solve_hkl_square_dual_fast(eta,data,Y,lambda,varargin);

% solve the milasso for a given regularization parameter
n = size(data.X,1);
p =length(data.affinity);
pX = size(data.X,2);
affs = data.affinity;
weights = data.weights;

% optional parameters
mingap = 1e-3;
display = 0;
gapeta = 1e-3;

% READ OPTIONAL PARAMETERS
args = varargin;
nargs = length(args);
for i=1:2:nargs
    switch args{i},
        case 'mingap',        mingap = args{i+1};
        case 'gapeta',        gapeta = args{i+1};
        case 'display',        display = args{i+1};
    end
end

if isempty(eta)
    eta = 1/p ./ weights.^2;
end
x = eta .* weights.^2;
x = max(x,0);
x = x / sum(x);

if display==1
display = Inf;
end
beta_ARMIJO = .25 ;
sigma_ARMIJO = .1 ;
kmax_ARMIJO = 20;
kmax = 200;
tol_dx = 1e-8;
k = 1;
alpha=1;

% % optimization parameters for the dual problem
% optparam.tol = 1e-16;
% optparam.kmax = 200;
% optparam.tol_df = 1e-20;
% optparam.display= 0;
fx0 = Inf;

% solve the regular kernel learning problem
mY = mean(Y);

%nu =  1e-6; % added ridge on kernels
while k<=kmax;
    x = max(x,0);
    x = x / sum(x);

    % compute value of function
    eta = ( x*(1-gapeta) + gapeta/p )./ weights.^2 ;
    zeta = zeros(p,1);
    for i=1:p
        zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
    end
    zeta = 1./zeta;

    % solve the regular kernel learning problem
    K = center_gram_matrix(double(devectorize_single( single(data.kernels * zeta ))));
    alphaK = (K + n * lambda * eye(n) ) \ ( Y - mean(Y) );
    fx =  lambda/2 * ( Y - mean(Y))'* alphaK;

    gradient_zeta = zeros(p,1);
    for i=1:p
        gradient_zeta(i) = - lambda/2 *  vectorize_quad_single(data.kernels(:,i),alphaK );
    end

    gradient_eta = zeros(p,1);
    for i=1:p;
        gradient_eta(i) = sum( gradient_zeta(affs{i}) .* zeta(affs{i}).^2 ) / eta(i)^2;
    end
    gradient = (1-gapeta) * gradient_eta ./ weights.^2;

    % check optimality condition
    optcond = max(x .* ( gradient - min(gradient)));
    if max(x .* ( gradient - min(gradient))) < mingap, break; end

    fx0=fx;

    % projected gradient
    xbar = x -   gradient;
    d = - gradient;
    x0=x;
    ialpha = 1;
    while ialpha<kmax_ARMIJO
        xalpha = simplex_projection(x + d * alpha);
        % fxalpha = feval(fun,xalpha,varargin{:});
        eta = ( xalpha*(1-gapeta) + gapeta/p )./ weights.^2 ;
        zeta = zeros(p,1);
        for i=1:p
            zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
        end
        zeta = 1./zeta;

    K = center_gram_matrix(double(devectorize_single( single(data.kernels * zeta ))));
    alphaK = (K + n * lambda * eye(n) ) \ ( Y - mean(Y) );
    fxalpha =  lambda/2 * ( Y - mean(Y))'* alphaK;

        if fx - fxalpha >= - sigma_ARMIJO * gradient' * ( xalpha - x );

            % start to go up
            while ( fx - fxalpha >= - sigma_ARMIJO * gradient' * ( xalpha - x ) ) & ialpha<kmax_ARMIJO;
                alpha = alpha / beta_ARMIJO;
                xalpha = simplex_projection(x + d * alpha);
                %fxalpha = feval(fun,xalpha,varargin{:});
                eta = ( xalpha*(1-gapeta) + gapeta/p )./ weights.^2 ;
                zeta = zeros(p,1);
                for i=1:p
                    zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
                end
                zeta = 1./zeta;

    K = center_gram_matrix(double(devectorize_single( single(data.kernels * zeta ))));
    alphaK = (K + n * lambda * eye(n) ) \ ( Y - mean(Y) );
    fxalpha =  lambda/2 * ( Y - mean(Y))'* alphaK;

                % fprintf('+');
                ialpha = ialpha+1;
            end
            alpha = alpha * beta_ARMIJO;
            xalpha = simplex_projection(x + d * alpha);
            x = xalpha;
            break;
        else
            alpha = alpha * beta_ARMIJO;

        end
        ialpha = ialpha+1;
    end

    if ialpha==kmax_ARMIJO, alpha=1; end

    if ~isinf(display)
        if display==1,
            fprintf('k=%d - f=%f - armijo=%d - dx=%e - gap=%f\n',k,fx,ialpha,norm(x-x0),optcond);
        end
        if display>1,
            if mod(k,display)==1
                fprintf('k=%d - f=%f - armijo=%d - dx=%e - gap=%f\n',k,fx,ialpha,norm(x-x0),optcond);

            end
        end
    end
    if norm(x-x0)<tol_dx,  break; end

    if sum(x) > 1+1e-1, keyboard; end

    k=k+1;
end

% if k==kmax+1
%     keyboard;
% end

eta = ( x*(1-gapeta) + gapeta/p )./ weights.^2 ;
zeta = zeros(p,1);
for i=1:p
    zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
end
zeta = 1./zeta;

    K = center_gram_matrix(double(devectorize_single( single(data.kernels * zeta ))));
    alpha = (K + n * lambda * eye(n) ) \ ( Y - mean(Y) );

    pred = K * alpha;
    b=mean(Y-pred );

kappa = zeros(p,p);
for i=1:p
    kappa(i,affs{i}) = zeta(affs{i})./eta(i);
end

 if isinf(display)
            fprintf('k=%d - f=%f -  gap=%f\n',k,fx,optcond);
        end

repo: REST


text:
% Test file for @deltafun/times.m.

function pass = test_times(pref)

if (nargin < 1)
    pref = chebfunpref();
end
%%
% Get the tolerance:
dTol = pref.deltaPrefs.deltaTol;

d = deltafun(bndfun(@sin), struct('deltaMag', 1, 'deltaLoc', 0));
pass(1) = isempty(deltafun() .* deltafun());
pass(2) = isempty(deltafun() .* d) && isempty(d .* deltafun());

f = fun.constructor(@(x) exp(-x));
g = bndfun(@sin);
df1 = deltafun(f, []);
df2 = deltafun(g, struct('deltaMag', [0; 0; 0; 0; 1], 'deltaLoc', 0));
s = df1.*df2;
pass(3) = norm(s.deltaMag - [1, 4, 6, 4, 1].', inf) < dTol;

f = fun.constructor(@(x) exp(x));
g = bndfun(@sin);
df1 = deltafun(f, []);
df2 = deltafun(g, struct('deltaMag', [0; 0; 0; 1], 'deltaLoc', 0));
s = df1.*df2;
pass(4) = norm(s.deltaMag - [-1, 3, -3, 1].', inf) < dTol;

a = -4;
b = 4;

f1 = fun.constructor(@(x) exp(sin(x)), struct('domain', [a, b]));
d1 = .5*[0 1 0;
         0 1 0;
         1 0 1; ];
l1 = sort(.9*(a + (b-a)/2*rand(1,3)));

f2 = fun.constructor(@(x) exp(cos(x)), struct('domain', [a, b]));

df1 = deltafun(f1, struct('deltaMag', d1, 'deltaLoc', l1));
df2 = deltafun(f2, []);

s = df1 .* df2;
pass(5) = iszero(s.funPart - f1.*f2);
pass(6) = norm(s.deltaLoc - sort(l1), inf) == 0;

c1 = [ feval(diff(f2,2), l1(1));
       -2*feval(diff(f2,1), l1(1));
       feval(diff(f2,0), l1(1));    ];

c2 = [ feval(diff(f2,0)-diff(f2,1), l1(2));
       feval(diff(f2,0), l1(2));
       0;                       ];

c3 = [feval(diff(f2,2), l1(3));
      -2*feval(diff(f2,1), l1(3));
      feval(diff(f2,0), l1(3))];

deltas1 = .5*[c1, c2, c3];

err = s.deltaMag - deltas1;
pass(7) = norm(err(:), inf) < dTol;

d2 = .8*[1 0 1;
         0 1 0;
         0 1 0; ];
l2 =  .0010751808168034 + sort(.9*(b-a)/2*rand(1,3));
df2 = deltafun(f2, struct('deltaMag', d2, 'deltaLoc',l2));
%%
s = df1 .* df2;
pass(8) = iszero(s.funPart - f1.*f2);
pass(9) = norm(s.deltaLoc - sort(union(l1,l2)), inf) == 0;
%%
c1 = [ feval(diff(f1,0), l2(1));
       0;
       0; ];
c2 = [ feval(diff(f1, 2)-diff(f1, 1), l2(2));
       feval(diff(f1, 0)-2*diff(f1, 1), l2(2));
       feval(diff(f1, 0), l2(2));  ];

c3 = [feval(diff(f1,0), l2(3));
      0;
      0; ];

deltas2 = .8*[c1, c2, c3];

err = [deltas1, deltas2] - s.deltaMag;
pass(10) = norm(err(:), inf) < dTol;
%%
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [results, success, raw] = opf_execute(om, mpopt)
%OPF_EXECUTE  Executes the OPF specified by an OPF model object.
%   [RESULTS, SUCCESS, RAW] = OPF_EXECUTE(OM, MPOPT)
%
%   RESULTS are returned with internal indexing, all equipment
%   in-service, etc.
%
%   See also OPF, OPF_SETUP.

%   MATPOWER
%   Copyright (c) 2009-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%%-----  setup  -----
%% options
dc  = strcmp(upper(mpopt.model), 'DC');
alg = upper(mpopt.opf.ac.solver);
switch alg
    case {'PDIPM', 'TRALM', 'MINOPF', 'SDPOPF'}
        legacy_solver = 1;
    otherwise
        legacy_solver = 0;
end
sdp = strcmp(alg, 'SDPOPF');
vcart = ~dc && mpopt.opf.v_cartesian;

%% get indexing
[vv, ll, nne, nni] = om.get_idx();

if mpopt.verbose > 0
    v = mpver('all');
    fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
end

if dc
    %%-----  run DC OPF solver  -----
    if mpopt.verbose > 0
        fprintf(' -- DC Optimal Power Flow\n');
    end
    [results, success, raw] = dcopf_solver(om, mpopt);
else
    %%-----  run AC OPF solver  -----
    if mpopt.verbose > 0
        fprintf(' -- AC Optimal Power Flow\n  AC OPF formulation: ');
        if sdp
            fprintf('SDP relaxation\n');
        else
            if vcart
                v = 'cartesian';
            else
                v = 'polar';
            end
            if mpopt.opf.current_balance
                v2 = 'current';
            else
                v2 = 'power';
            end
            fprintf('%s voltages, %s balance eqns\n', v, v2);
        end
    end

    %% ZIP loads?
    if legacy_solver && ( ...
            (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
              any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ...
            (~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
              any(mpopt.exp.sys_wide_zip_loads.qw(2:3))) )
        warning('opf_execute: ''%s'' solver does not support ZIP load model. Converting to constant power loads.', alg)
        mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ...
                          struct('pw', [], 'qw', []));
    end

    %% run specific AC OPF solver
    if legacy_solver
        switch alg
            case 'PDIPM'
                if mpopt.pdipm.step_control
                    if ~have_feature('scpdipmopf')
                        error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires SCPDIPMOPF (see http://www.pserc.cornell.edu/tspopf/)', alg);
                    end
                else
                    if ~have_feature('pdipmopf')
                        error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires PDIPMOPF (see http://www.pserc.cornell.edu/tspopf/)', alg);
                    end
                end
                [results, success, raw] = tspopf_solver(om, mpopt);
            case 'TRALM'
                if ~have_feature('tralmopf')
                    error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires TRALM (see http://www.pserc.cornell.edu/tspopf/)', alg);
                end
                [results, success, raw] = tspopf_solver(om, mpopt);
            case 'MINOPF'
                if ~have_feature('minopf')
                    error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires MINOPF (see http://www.pserc.cornell.edu/minopf/)', alg);
                end
                [results, success, raw] = mopf_solver(om, mpopt);
            case 'SDPOPF'
                if ~have_feature('yalmip')
                    error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires YALMIP (see https://yalmip.github.io)', alg);
                end
                [results, success, raw] = sdpopf_solver(om, mpopt);
            otherwise
                error('opf_execute: MPOPT.opf.ac.solver = ''%s'' is not a valid AC OPF solver selection', alg);
        end
    else    %% not legacy solver
        switch alg
            case 'IPOPT'
                if ~have_feature('ipopt')
                    error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires IPOPT (see https://github.com/coin-or/Ipopt)', alg);
                end
            case 'FMINCON'
                if ~have_feature('fmincon')
                    error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires FMINCON (Optimization Toolbox 2.x or later)', alg);
                end
            case 'KNITRO'
                if ~have_feature('knitro')
                    error('opf_execute: MPOPT.opf.ac.solver = ''%s'' requires Artelys Knitro (see https://www.artelys.com/solvers/knitro/)', alg);
                end
        end
        [results, success, raw] = nlpopf_solver(om, mpopt);
    end     %% if legacy_solver
end %% if dc
if ~isfield(raw, 'output') || ~isfield(raw.output, 'alg') || isempty(raw.output.alg)
    raw.output.alg = alg;
end

if success
  if ~dc && ~sdp
    %% copy bus voltages back to gen matrix
    results.gen(:, VG) = results.bus(results.gen(:, GEN_BUS), VM);

    %% cartesian voltage magnitude multipliers
    if vcart
        results.bus(:, MU_VMIN) = results.bus(:, MU_VMIN) .* results.bus(:, VM) * 2;
        results.bus(:, MU_VMAX) = results.bus(:, MU_VMAX) .* results.bus(:, VM) * 2;
    end

    %% gen PQ capability curve multipliers
    if ll.N.PQh > 0 || ll.N.PQl > 0
      mu_PQh = results.mu.lin.l(ll.i1.PQh:ll.iN.PQh) - results.mu.lin.u(ll.i1.PQh:ll.iN.PQh);
      mu_PQl = results.mu.lin.l(ll.i1.PQl:ll.iN.PQl) - results.mu.lin.u(ll.i1.PQl:ll.iN.PQl);
      Apqdata = om.get_userdata('Apqdata');
      results.gen = update_mupq(results.baseMVA, results.gen, mu_PQh, mu_PQl, Apqdata);
    end

    %% compute g, dg, f, df, d2f if requested by opf.return_raw_der = 1
    if mpopt.opf.return_raw_der
      %% move from results to raw if using v4.0 of MINOPF or TSPOPF
      if isfield(results, 'dg')
        raw.dg = results.dg;
        raw.g = results.g;
      end
      %% compute g, dg, unless already done by post-v4.0 MINOPF or TSPOPF
      if ~isfield(raw, 'dg')
        [g, geq, dg, dgeq] = nlp_consfcn(om, results.x);
        raw.g = [ geq; g];
        raw.dg = [ dgeq'; dg'];   %% true Jacobian organization
      end
      %% compute df, d2f
      [f, df, d2f] = nlp_costfcn(om, results.x);
      raw.df = df;
      raw.d2f = d2f;
    end
  end

  %% delete g and dg fields from results if using v4.0 of MINOPF or TSPOPF
  if isfield(results, 'dg')
    rmfield(results, 'dg');
    rmfield(results, 'g');
  end

  %% angle limit constraint multipliers
  iang = om.get_userdata('iang');
  if ~sdp && length(iang)
    if vcart
      iang = om.get_userdata('iang');
      results.branch(iang, MU_ANGMIN) = results.mu.nli(nni.i1.angL:nni.iN.angL) * pi/180;
      results.branch(iang, MU_ANGMAX) = results.mu.nli(nni.i1.angU:nni.iN.angU) * pi/180;
    else
      iang = om.get_userdata('iang');
      results.branch(iang, MU_ANGMIN) = results.mu.lin.l(ll.i1.ang:ll.iN.ang) * pi/180;
      results.branch(iang, MU_ANGMAX) = results.mu.lin.u(ll.i1.ang:ll.iN.ang) * pi/180;
    end
  end
else
  %% assign empty g, dg, f, df, d2f if requested by opf.return_raw_der = 1
  if ~dc && mpopt.opf.return_raw_der
    raw.dg = [];
    raw.g = [];
    raw.df = [];
    raw.d2f = [];
  end
end

if ~sdp
  %% assign values and limit shadow prices for variables
  om_var_order = om.get('var', 'order');
  for k = 1:length(om_var_order)
    name = om_var_order(k).name;
    if om.getN('var', name)
      idx = vv.i1.(name):vv.iN.(name);
      results.var.val.(name) = results.x(idx);
      results.var.mu.l.(name) = results.mu.var.l(idx);
      results.var.mu.u.(name) = results.mu.var.u(idx);
    end
  end

  %% assign shadow prices for linear constraints
  om_lin_order = om.get('lin', 'order');
  for k = 1:length(om_lin_order)
    name = om_lin_order(k).name;
    if om.getN('lin', name)
      idx = ll.i1.(name):ll.iN.(name);
      results.lin.mu.l.(name) = results.mu.lin.l(idx);
      results.lin.mu.u.(name) = results.mu.lin.u(idx);
    end
  end

  %% assign shadow prices for nonlinear constraints
  if ~dc
    om_nle_order = om.get('nle', 'order');
    for k = 1:length(om_nle_order)
      name = om_nle_order(k).name;
      if om.getN('nle', name)
        results.nle.lambda.(name) = results.mu.nle(nne.i1.(name):nne.iN.(name));
      end
    end

    om_nli_order = om.get('nli', 'order');
    for k = 1:length(om_nli_order)
      name = om_nli_order(k).name;
      if om.getN('nli', name)
        results.nli.mu.(name) = results.mu.nli(nni.i1.(name):nni.iN.(name));
      end
    end
  end

  %% assign values for components of quadratic cost
  om_qdc_order = om.get('qdc', 'order');
  for k = 1:length(om_qdc_order)
    name = om_qdc_order(k).name;
    if om.getN('qdc', name)
      results.qdc.(name) = om.eval_quad_cost(results.x, name);
    end
  end

  %% assign values for components of general nonlinear cost
  om_nlc_order = om.get('nlc', 'order');
  for k = 1:length(om_nlc_order)
    name = om_nlc_order(k).name;
    if om.getN('nlc', name)
      results.nlc.(name) = om.eval_nln_cost(results.x, name);
    end
  end

  %% assign values for components of legacy user cost
  om_cost_order = om.get('cost', 'order');
  for k = 1:length(om_cost_order)
    name = om_cost_order(k).name;
    if om.getN('cost', name)
      results.cost.(name) = om.eval_legacy_cost(results.x, name);
    end
  end

  %% if single-block PWL costs were converted to POLY, insert dummy y into x
  %% Note: The "y" portion of x will be nonsense, but everything should at
  %%       least be in the expected locations.
  pwl1 = om.get_userdata('pwl1');
  if ~isempty(pwl1) && ~strcmp(alg, 'TRALM') && ~(strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control)
    %% get indexing
    vv = om.get_idx();
    if dc
      nx = vv.iN.Pg;
    else
      nx = vv.iN.Qg;
    end
    y = zeros(length(pwl1), 1);
    raw.xr = [ raw.xr(1:nx); y; raw.xr(nx+1:end)];
    results.x = [ results.x(1:nx); y; results.x(nx+1:end)];
  end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% SB2_LIKELIHOODS  Convenience function to encapsulate likelihood types
%
% LIKELIHOOD = SB2_LIKELIHOODS(TYPE)
%
% OUTPUT ARGUMENTS:
%
%	LIKELIHOOD	A structure containing an enumeration of possible
%				likelihood types, and specification of the type in use.
%
% INPUT ARGUMENTS:
%
%	TYPE		A string. One of:
%					'Gaussian'
%					'Bernoulli'
%					'Poisson'
%
% EXAMPLE:
%
%	LIKELIHOOD	= SB2_Likelihoods('Gaussian');
%
% NOTES:
%
% This function is "for internal use only" and exists purely for convenience
% of coding and to facilitate potential future expansion.
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the Free
% Software Foundation; either version 2 of the License, or (at your option)
% any later version.
%
% SPARSEBAYES 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 SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function LIKELIHOOD = SB2_Likelihoods(likelihood_)

% For convenience in code, and later expandability,
% we "enum" the available likelihood models
%
% Continuous
LIKELIHOOD.Gaussian		= 1;
% Discrete
LIKELIHOOD.Bernoulli	= 2;
LIKELIHOOD.Poisson		= 3;

% Feel free to add your own ... and work out the maths :-)

% Determine the likelihood to be used
%
LIKELIHOOD.InUse		= 0;
%
switch upper(likelihood_(1:5))
 case 'GAUSS',
  LIKELIHOOD.InUse	= LIKELIHOOD.Gaussian;
  %
 case 'BERNO',
  LIKELIHOOD.InUse	= LIKELIHOOD.Bernoulli;
  %
 case 'POISS',
  LIKELIHOOD.InUse	= LIKELIHOOD.Poisson;
  %
 otherwise
  error('Unknown likelihood ''%s'' specified\n', likelihood_)
end


text:
function bigpot = divide_by_pot(bigpot, smallpot)
% DIVIDE_BY_POT bigpot /= smallpot for cpot
% bigpot = divide_by_pot(bigpot, smallpot)
%
% smallpot's domain must be a subset of bigpot's domain.
if bigpot.g ==-Inf && smallpot.g == -Inf %bug fix by Sacha Gunaratne 2017-10-23
    bigpot.g = 0;
else
    bigpot.g = bigpot.g - smallpot.g;
end
if sum(smallpot.sizes) > 0
  mask = find_equiv_posns(smallpot.domain, bigpot.domain);
  u = block(mask, bigpot.sizes);
  bigpot.h(u) = bigpot.h(u) - smallpot.h;
  bigpot.K(u, u) = bigpot.K(u, u) - smallpot.K;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%
% FILE NAME:    curateLipidPools.m
%
% PURPOSE: This script revises the formulation of reactions involving lipid
%          and cholesterol-ester pool formation and breakdown in HumanGEM
%          to achieve a more generic and mass-balanced model. In the
%          process, the following changes were made:
%             - 6 new pool reactions
%             - 2 new transport reaction
%             - 16 new pool metabolites (3 unique, ignoring compartment)
%             - 23 pool metabolites were updated to the new metabolites
%
%          Full details on model changes can be found in the generated
%          model change log files (lipidPool_modelChanges_mets.tsv and
%          lipidPool_modelChanges_rxns.tsv)
%
%          This script also removes 3 model fields that are now deprecated:
%             1. proteins
%             2. prRules
%             3. rxnProtMat
%
%          The protein-related fields can simply be regenerated using the
%          translateGrRules function.
%

%% Load Human-GEM and annotation files

% load Human-GEM model
load('HumanGEM.mat');
ihuman_orig = ihuman;  % keep copy of original version

% load metabolite and reaction annotation data
metAssoc = jsondecode(fileread('humanGEMMetAssoc.JSON'));
rxnAssoc = jsondecode(fileread('humanGEMRxnAssoc.JSON'));
changeNotes = {};

% verify that HumanGEM and annotation structures are aligned
if ~isequal(ihuman.mets, metAssoc.mets) || ~isequal(ihuman.rxns, rxnAssoc.rxns)
    error('HumanGEM is not synced with the metAssoc and/or rxnAssoc structure!');
end

%% Remove three deprecated model fields

ihuman = rmfield(ihuman, {'proteins';'prRules';'rxnProtMat'});

%% Add new metabolites

% first need to rename an existing "fatty acid pool" metabolite to avoid confusion
ihuman.metNames(ismember(ihuman.metNames,'fatty acid pool')) = {'fatty acid biomass pool'};

% load new metabolite information from file
fid = fopen('../../ComplementaryData/modelCuration/lipidPools/newLipidPoolMets.tsv');
metData = textscan(fid,'%s%s%s%s%d%s%s%s%s%s%s','Delimiter','\t','Headerlines',1);
fclose(fid);

% verify that none of the metabolite IDs (ignoring compartment) exist in the current model
if any(startsWith(ihuman.mets, regexprep(metData{1},'.$','')))
    error('One or more metabolite IDs to be added already exist in the model.');
end

% add metabolites to the model
metsToAdd = {};
metsToAdd.mets = metData{1};
metsToAdd.metNames = metData{2};
metsToAdd.compartments = metData{3};
metsToAdd.metFormulas = metData{4};
metsToAdd.metCharges = metData{5};
ihuman = addMets(ihuman, metsToAdd);

% add metabolites to the metAssoc structure
numOrigMets = numel(metAssoc.mets);
numNewMets = numel(metData{1});
newMetInd = (numOrigMets+1:numOrigMets+numNewMets)';
f = fieldnames(metAssoc);
for i = 1:numel(f)
    % initialize structure with empty entries
    metAssoc.(f{i})(newMetInd) = {''};
end
metAssoc.mets(newMetInd) = metData{1};
metAssoc.metsNoComp(newMetInd) = regexprep(metData{1},'.$','');
metAssoc.metKEGGID(newMetInd) = metData{6};
metAssoc.metBiGGID(newMetInd) = metData{7};
metAssoc.metChEBIID(newMetInd) = metData{8};
metAssoc.metMetaNetXID(newMetInd) = metData{9};
metAssoc.metLipidMapsID(newMetInd) = metData{10};

changeNotes = [changeNotes; [metsToAdd.mets, metData{11}]];

%% Add new pool reactions

% load new reaction information from file
fid = fopen('../../ComplementaryData/modelCuration/lipidPools/newLipidPoolRxns.tsv');
rxnData = textscan(fid,'%s%s%s%s','Delimiter','\t','Headerlines',1);
fclose(fid);

% verify that none of the reactions exist in the current model
if any(ismember(rxnData{1},ihuman.rxns))
    error('One or more reactions to be added already exist in the model.');
end

% add reactions to the model
rxnsToAdd = {};
rxnsToAdd.rxns = rxnData{1};
rxnsToAdd.rxnNames = rxnData{2};
rxnsToAdd.subSystems = cellfun(@(s) {{s}},rxnData{3});
rxnsToAdd.equations = rxnData{4};
ihuman = addRxns(ihuman, rxnsToAdd, 3);
ihuman.priorCombiningGrRules(end+1:end+numel(rxnsToAdd.rxns)) = {''};

% add reactions to the rxnAssoc structure
numOrigRxns = numel(rxnAssoc.rxns);
numNewRxns = numel(rxnData{1});
newRxnInd = (numOrigRxns+1:numOrigRxns+numNewRxns)';
f = fieldnames(rxnAssoc);
for i = 1:numel(f)
    % initialize structure with empty entries
    rxnAssoc.(f{i})(newRxnInd) = {''};
end
rxnAssoc.rxns(newRxnInd) = rxnData{1};

%% Update cholesterol-ester pool metabolites
% The "cholesterol-ester pool" metabolites will be updated to the new
% "cholesterol-ester plasma pool" metabolites ONLY in the old cholesterol-
% ester pool formation/degradation reactions (HMR_3537 and HMR_3622). All
% other instances of the "cholesterol-ester pool" metabolite will remain
% unchanged.

% HMR_3537 (lysosome)
oldMetInd = getIndexes(ihuman, 'cholesterol-ester pool[l]', 'metscomps');
newMetInd = getIndexes(ihuman, 'cholesterol-ester plasma pool[l]', 'metscomps');
rxnInd = getIndexes(ihuman,'HMR_3537','rxns');
ihuman.S(oldMetInd, rxnInd) = 0;
ihuman.S(newMetInd, rxnInd) = -1;

% HMR_3622 (endoplasmic reticulum)
oldMetInd = getIndexes(ihuman, 'cholesterol-ester pool[r]', 'metscomps');
newMetInd = getIndexes(ihuman, 'cholesterol-ester plasma pool[r]', 'metscomps');
rxnInd = getIndexes(ihuman,'HMR_3622','rxns');
ihuman.S(oldMetInd, rxnInd) = 0;
ihuman.S(newMetInd, rxnInd) = -1;

%% Update fatty acid, 1-acylglycerol-3P, and acyl-CoA pool metabolites
% The old version of these pool metabolites will be updated to their new
% versions in ALL reactions EXCEPT those that involve their formation/
% degradation. The formation/degradation reactions involving the old pool
% metabolites will remain unchanged.

% load list of existing pool reactions related to lipid pools
fid = fopen('../../ComplementaryData/modelCuration/lipidPools/oldLipidPoolRxns.tsv');
filedata = textscan(fid,'%s','Headerlines',1);
fclose(fid);
oldPoolRxns = filedata{1};
oldPoolRxnInd = getIndexes(ihuman, oldPoolRxns, 'rxns');

% load array specifying the conversion from old lipid pool mets to new mets
fid = fopen('../../ComplementaryData/modelCuration/lipidPools/lipidPoolMetReplacements.tsv');
filedata = textscan(fid,'%s%s','Delimiter','\t','Headerlines',1);
fclose(fid);
poolMets = [filedata{1}, filedata{2}];

% replace old pool metabolites with new pool metabolites
for i = 1:size(poolMets,1)
    metComps = unique(ihuman.metComps(ismember(ihuman.metNames, poolMets(i,1))));
    for j = 1:numel(metComps)
        oldMetInd = find(ismember(ihuman.metNames,poolMets(i,1)) & ihuman.metComps == metComps(j));
        newMetInd = find(ismember(ihuman.metNames,poolMets(i,2)) & ihuman.metComps == metComps(j));
        rxnInd = setdiff(find(ihuman.S(oldMetInd,:) ~= 0), oldPoolRxnInd);
        if ~isempty(rxnInd)
            stoichCoeffs = ihuman.S(oldMetInd, rxnInd);
            ihuman.S(oldMetInd, rxnInd) = 0;
            ihuman.S(newMetInd, rxnInd) = stoichCoeffs;
        end
    end
end

% Rename the old pool metabolites, specifying that their composition is
% specific to liver tissue (see Mardinoglu et al. Nat Commun 2014, PMID:24419221)
for i = 1:size(poolMets,1)
    metInd = getIndexes(ihuman, poolMets(i,1), 'metnames');
    ihuman.metNames(metInd) = regexprep(ihuman.metNames(metInd), 'pool$', 'pool (liver tissue)');
end

% Remove old pool metabolites that are no longer used
ihuman = removeReactions(ihuman,[],true);

% Remove 23 unused pool (12 unique) mets from the metAssoc structure
removedMets = setdiff(metAssoc.mets, ihuman.mets);
removedMetsIdx = find(ismember(metAssoc.mets, removedMets));
f = fieldnames(metAssoc);
for i = 1:numel(f)
    metAssoc.(f{i})(removedMetsIdx) = [];
end

%% Identify and document changes to the model

% verify that HumanGEM and annotation structures are aligned
if ~isequal(ihuman.mets, metAssoc.mets) || ~isequal(ihuman.rxns, rxnAssoc.rxns)
    error('HumanGEM is not synced with the metAssoc and/or rxnAssoc structure!');
end

% export reaction and metabolite annotation structures to JSON
jsonStr = jsonencode(rxnAssoc);
fid = fopen('../../ComplementaryData/annotation/humanGEMRxnAssoc.JSON', 'w');
fwrite(fid, prettyJson(jsonStr));
fclose(fid);
jsonStr = jsonencode(metAssoc);
fid = fopen('../../ComplementaryData/annotation/humanGEMMetAssoc.JSON', 'w');
fwrite(fid, prettyJson(jsonStr));
fclose(fid);

% determine and document model changes
modelChanges = docModelChanges(ihuman_orig,ihuman,changeNotes);
writeModelChanges(modelChanges,'../../ComplementaryData/modelCuration/lipidPools/lipidPool_modelChanges.tsv');

% export HumanGEM
exportHumanGEM(ihuman,'HumanGEM','../../',{'mat','yml'},false,false);

% clear unneeded variables
clearvars -except ihuman_orig ihuman modelChanges

author: SysBioChalmers

repo: Human-GEM