text:
function funReset1(object_handle, event, start_button1, pause_button1, save_button1)

    global mode
    mode = 'Reset1';

    set(start_button1, 'BackGroundColor', 'g', 'Enable', 'on');
    set(object_handle, 'BackGroundColor', [0.4,0.4,0.4], 'Enable', 'off');
    set(pause_button1, 'BackGroundColor', [0.4,0.4,0.4], 'Enable', 'off');

    set(save_button1, 'BackGroundColor', [0.4,0.4,0.4], 'Enable', 'off');

end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
function out = ctranspose3x3(in)

% compute ctranspose of multiple 3x3 matrices, input is 3x3xN

out = conj(in);
out(1,2,:,:) = conj(in(2,1,:,:));
out(2,1,:,:) = conj(in(1,2,:,:));
out(1,3,:,:) = conj(in(3,1,:,:));
out(3,1,:,:) = conj(in(1,3,:,:));
out(2,3,:,:) = conj(in(3,2,:,:));
out(3,2,:,:) = conj(in(2,3,:,:));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s = i4_to_s_zero ( intval, s_len )

%*****************************************************************************80
%
%% I4_TO_S_ZERO converts an I4 to a string, with zero padding.
%
%  Discussion:
%
%    An I4 is an integer.
%
%  Example:
%
%    Assume that S is 6 characters long:
%
%    INTVAL  S
%
%         1  000001
%        -1  -00001
%         0  000000
%      1952  001952
%    123456  123456
%   1234567  ******  <-- Not enough room%
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer INTVAL, an integer to be converted.
%
%    Input, integer S_LEN, the length of the string to be used.
%
%    Output, string S, the representation of the integer.
%    The integer will be right justified, and zero padded.
%    If there is not enough space, the string will be filled with stars.
%
  s = [];

  if ( s_len <= 0 )
    return
  end

  ilo = 1;
%
%  Make a copy of the integer.
%
  ival = intval;
%
%  Handle the negative sign.
%
  if ( ival < 0 )

    if ( s_len <= 1 )
      s(1) = '*';
      s = char ( s );
      return
    end

    ival = - ival;
    s(1) = '-';
    ilo = 2;

  end
%
%  Working from right to left, strip off the digits of the integer
%  and place them into S(ILO:S_LEN).
%
  ipos = s_len;

  while ( ival ~= 0 | ipos == s_len )

    idig = mod ( ival, 10 );
    ival = floor ( ival / 10 );

    if ( ipos < ilo )
      s(1:s_len) = '*';
      s = char ( s );
      return
    end

    c = digit_to_ch ( idig );

    s(ipos) = c;
    ipos = ipos - 1;

  end
%
%  Fill the empties with zeroes.
%
  s(ilo:ipos) = '0';
  s = char ( s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ transform ] = genTransMat( front, up )
%GENTRANSMAT

axes = cross(front,up);
axes = axes/norm(axes,2);
transform = [front(:)';up(:)';axes(:)'];

end

repo: GRAINS


text:
%----UPDATE AXIS DATA/LAYOUT----%

function obj = updateAxisMultipleYAxes(obj,axIndex,yaxIndex)

    %-STANDARDIZE UNITS-%
    axisUnits = get(obj.State.Axis(axIndex).Handle,'Units');
    set(obj.State.Axis(axIndex).Handle,'Units','normalized')

    try
        fontUnits = get(obj.State.Axis(axIndex).Handle,'FontUnits');
        set(obj.State.Axis(axIndex).Handle,'FontUnits','points')
    catch
        % TODO
    end

    %-AXIS DATA STRUCTURE-%
    axisData = get(obj.State.Axis(axIndex).Handle);

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

    %-xaxis-%
    xaxis = extractAxisData(obj,axisData, 'X');

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

    %-yaxis-%
    [yaxis, yAxisLim] = extractAxisDataMultipleYAxes(obj, axisData, yaxIndex);

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

    %-getting and setting postion data-%

    xo = axisData.Position(1);
    yo = axisData.Position(2);
    w = axisData.Position(3);
    h = axisData.Position(4);

    if obj.PlotOptions.AxisEqual
        wh = min(axisData.Position(3:4));
        w = wh;
        h = wh;
    end

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

    %-xaxis domain-%
    xaxis.domain = min([xo xo + w],1);
    scene.domain.x = min([xo xo + w],1);

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

    %-yaxis domain-%
    yaxis.domain = min([yo yo + h],1);
    scene.domain.y = min([yo yo + h],1);

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

    [xsource, ysource, xoverlay, yoverlay] = findSourceAxis(obj, axIndex, yaxIndex);

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

    %-xaxis anchor-%
    xaxis.anchor = ['y' num2str(ysource)];

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

    %-yaxis anchor-%
    yaxis.anchor = ['x' num2str(xsource)];

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

    %-xaxis overlaying-%
    if xoverlay
        xaxis.overlaying = ['x' num2str(xoverlay)];
    end

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

    %-yaxis overlaying-%
    if yoverlay
        yaxis.overlaying = ['y' num2str(yoverlay)];
    end

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

    % update the layout field (do not overwrite source)
    if xsource == axIndex
        obj.layout = setfield(obj.layout,['xaxis' num2str(xsource)],xaxis);
    end

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

    % update the layout field (do not overwrite source)
    obj.layout = setfield(obj.layout,['yaxis' num2str(ysource)],yaxis);

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

    %-REVERT UNITS-%
    set(obj.State.Axis(axIndex).Handle,'Units',axisUnits);

    try
        set(obj.State.Axis(axIndex).Handle,'FontUnits',fontUnits);
    catch
        % TODO
    end

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

    %-do y-axes visibles-%
    obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
    plotIndex = obj.PlotOptions.nPlots;

    obj.data{plotIndex}.type = 'scatter';
    obj.data{plotIndex}.xaxis = ['x' num2str(xsource)];
    obj.data{plotIndex}.yaxis = ['y' num2str(ysource)];

    %-------------------------------------------------------------------------%
end

author: plotly


text:
function C = ojw_bsxfun(func, A, B)
%OJW_BSXFUN  Provides backwards compatibility for Matlab's bsxfun
%
%   C = ojw_bsxfun(func, A, B)
%
% Each dimension of A and B must be equal to each other, or equal to one.
% Whenever a dimension of one of A or B is singleton (equal to 1), the
% array is virtually replicated along that dimension to match the other
% array (or diminished if the corresponding dimension of the other array is
% 0). The two resized arrays are then passed, as vectors, to func(), which
% must be able to accept as input either two column vectors of the same
% size, or one column vector and one scalar, and return as output a column
% vector of the same size as the input(s).
%
% IN:
%    func - Function handle.
%    A - Input array.
%    B - Input array of same type as A.
%
% OUT:
%    C - max(size(A),size(B)).*(size(A)>0 & size(B)>0) sized output array.
%
% See also BSXFUN.

% $Id: ojw_bsxfun.m,v 1.1 2007/12/07 11:27:52 ojw Exp $

if exist('bsxfun', 'builtin')
    % Do the obvious - Matlab versions below 7.3 don't have it though
    C = bsxfun(func, A, B);
else
    % Backwardly compatible method
    if isscalar(A)
        C = reshape(func(A, B(:)), size(B));
        return
    elseif isscalar(B)
        C = reshape(func(A(:), B), size(A));
        return
    end
    sA = size(A);
    sB = size(B);
    % Enlarge the smaller of the two sizes
    s = numel(sA) - numel(sB);
    if s < 0
        sA = [sA ones(1, -s)];
    elseif s > 0
        sB = [sB ones(1, s)];
    end
    % Calculate the output array size
    sMax = max(sA, sB) .* (sA > 0 & sB > 0);
    % Make sure arrays have same size of dimension or size 1
    a = sA ~= sMax;
    b = sB ~= sMax;
    if any(a & sA ~= 1) || any(b & sB ~= 1)
        error('A and B must have equal or singleton dimensions');
    end
    if ~all(sMax)
        % Some entries are zero, so array is empty
        C = zeros(sMax);
        return
    end
    % Resize the arrays to have the same size, sMax
    if any(a)
        A = repmat(A, sMax ./ sA);
    end
    if any(b)
        B = repmat(B, sMax ./ sB);
    end
    % Apply the function
    C = reshape(func(A(:), B(:)), sMax);
end

repo: SelfExSR


text:
function [roitSeries, subCoords] = getTseriesOneROI(vw, ROIcoords, scanNum, getRawTseries, removeRedundant)
%
%   [roitSeries, subCoords] = getTseriesOneROI(vw,ROIcoords,scanNum, getRawTseriesFlag(=0 default), removeRedundantFlag(=1 default) )
%
% For a single ROI, we loop through slices and acquire all of the
% tSeries(voxel by voxel). Return these time courses, removing redundant
% tSeries (which occur because several anatomical coords may reference the
% same functional voxel). Also returns the locations of the non-redundant
% voxels in subCoords.
% The removing redundant can be turned off by removeRedundantFlag = 0;
%
% ras, 04/04: remove redundant time courses from upsampling the functional.
% ras, 05/04: also, return the sub-coords, with the redundant voxels
% removed. (They'll be in the coordinates of the inplane anatomy, but
% there will only be one voxel contained within a given functional voxel.)
% ras, 5/2004 -- now it returns the tSeries in the order that corresponds
% to the voxels specified in ROIcoords (before, it would return in them
% in a semi-shuffled order if you were in the inplane view -- b/c of
% the loop through slices).
% jl, 5/2005 - make "remove redundant from upsampling" optional

% Make sure the input is in the right format
if notDefined('ROIcoords'),         ROIcoords       = viewGet(vw, 'ROI coords'); end
if notDefined('scanNum'),           scanNum         = viewGet(vw, 'Current Scan'); end
if notDefined('getRawTseries'),     getRawTseries   = 0; end
if notDefined('removeRedundant'),   removeRedundant = 1; end

ROIcoords = ROIcoords2cellArray(vw,ROIcoords);

% b/c of upsampling b/w the functionals and anats,
% there tend to be redundant coords specifying the
% same functional voxel -- remove these:
switch vw.viewType
    case 'Inplane',
        if removeRedundant
            % b/c of upsampling b/w the functionals and anats, there tend
            % to be redundant coords specifying the same functional voxel
            % -- remove these:
            subCoords = roiSubCoords(vw,ROIcoords{1});
        end
        ROIcoords = {subCoords};
    case 'Flat',
        % the tSeries will only cover points
        % on the flat map mapped by the ROI,
        % ignoring ones on blank/masked nodes

        % jw: comment out
        %         ind = sub2ind(size(vw.indices),ROIcoords{1}(1,:),...
        %                         ROIcoords{1}(2,:),ROIcoords{1}(3,:));
        %         ind = find(vw.indices(ind)>0);
        %         subCoords = ROIcoords{1}(:,ind);

        % jw: comment in
        subCoords = ROIcoords;
    otherwise,
        subCoords = ROIcoords;
end

roitSeries = [];

if ~getRawTseries
    detrend = detrendFlag(vw,scanNum);
    smoothFrames = detrendFrames(vw,scanNum);
end

% Find the slice indices for this collection of ROIs
sliceInds = getSlicesROI(vw,ROIcoords);

for iSlice = 1:length(sliceInds)
    slice = sliceInds(iSlice);

    % Load tSeries & divide by mean, but don't detrend yet.
    % Otherwise, detrending the entire tSeries is much slower. DJH
    % this crashes on the volume - out of memory - kgs 06/05
    % rewrote this section
    if(getRawTseries)
        vw = percentTSeries(vw,scanNum,slice,0,0,0,1);
     else
        vw = percentTSeries(vw,scanNum,slice,0);
     end

    % Extract time-series from the current slice
    %  The third argument of getTSeriesROI is 'preserveCoords', which is
    %  the logical opposite of removeRedundant. That is, if we remove
    %  redundant voxels, we will not preserve the coordinates (we will get
    %  back fewer time series than coordinates), and vice versa.
    preserveCoords = ~removeRedundant;
    [subtSeries subIndices] = getTSeriesROI(vw,ROIcoords{1}, preserveCoords);

    if removeRedundant & strcmpi(viewGet(vw, 'view type'), 'gray')
        coords = viewGet(vw, 'coords');
        subCoords{1} = coords(:, subIndices);
    end

    if ~isempty(subtSeries) & ~getRawTseries %#ok<*AND2>
        % Norm and Detrend (faster and doesn't crash after extracting subtSeries for a small subset of the voxels)
        subtSeries = detrendTSeries(subtSeries,detrend,smoothFrames);
    end

    switch vw.viewType
        case 'Inplane',
            % assign to the columns in roitSeries that correspond to the selected
            % voxels:
            voxelsInSlice = ROIcoords{1}(3,:) == vw.tSeriesSlice;
            roitSeries(:,voxelsInSlice) = subtSeries;
        case 'Gray',
            roitSeries = subtSeries;
        case 'Volume',
            roitSeries = subtSeries;
        case 'Flat',
            roitSeries = [roitSeries subtSeries];
    end
end

% make a cell (to be in line with multi-ROI analyses)
roitSeries = {roitSeries};
clear vw.tSeries; % free memory

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% SERVER Write a message over the specified port
%
% Usage - server(message, output_port, number_of_retries)
function server(port, number_of_retries)

    import java.net.ServerSocket
    import java.io.*

    if (nargin < 2)
        number_of_retries = 2; % set to -1 for infinite
    end
    retry             = 0;

    server_socket  = [];
    output_socket  = [];

    port

    while true

        retry = retry + 1;

        try
            if ((number_of_retries > 0) && (retry > number_of_retries))
                fprintf(1, 'Too many retries\n');
                break;
            end

            fprintf(1, ['Try %d waiting for client to connect to this ' ...
                        'host on port : %d\n'], retry, port);

            % wait for 1 second for client to connect server socket
            server_socket = ServerSocket(port);
            server_socket.setSoTimeout(10000);

            fprintf(1, 'Waiting for client to connect\n');

            output_socket = server_socket.accept;

            fprintf(1, 'Client connected\n');

            input_stream   = output_socket.getInputStream;
            %d_output_stream = DataOutputStream(output_stream);

            for i=1:10
                % read data from the socket - wait a short time first
            pause(0.5);
            bytes_available = input_stream.available;
            fprintf(1, 'Reading %d bytes\n', bytes_available);
            end

            % clean up
            server_socket.close;
            output_socket.close;
            break;

        catch me
            me
            if ~isempty(server_socket)
                server_socket.close
            end

            if ~isempty(output_socket)
                output_socket.close
            end

            % pause before retrying
            pause(1);
   end
end
end

author: Allopart

repo: rbpf-gmapping


text:
function varargout = drawSquareMesh(nodes, edges, faces, varargin) %#ok<INUSL>
%DRAWSQUAREMESH Draw a 3D square mesh given as a graph.
%
%   drawSquareMesh(NODES, EDGES, FACES)
%   Draw the mesh defined by NODES, EDGES and FACES. FACES must be a N-by-4
%   array of vertex indices.
%
%   See also
%   boundaryGraph, drawGraph

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-06-28
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% input size check up
if size(faces, 2) ~= 4
    error('Requires a face array with 4 columns');
end

% number of faces
Nf = size(faces, 1);

% allocate memory for vertex coordinates
px = zeros(4, Nf);
py = zeros(4, Nf);
pz = zeros(4, Nf);

% initialize vertex coordinates of each face
for f = 1:Nf
    face = faces(f, 1:4);
    px(1:4, f) = nodes(face, 1);
    py(1:4, f) = nodes(face, 2);
    pz(1:4, f) = nodes(face, 3);
end

p = patch(px, py, pz, 'r');

if nargout > 0
    varargout = {p};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [operator_pt]=pt(operator, partition, dimensions)
%returns the partial transpose of the operator 'operator' with respect to
%the particles indicated by a one in the binary vector 'partition'.
%the optional parameter 'dimensions' includes the dimensions of the systems
%in array form. per default, it is assumed that each system is a qubit.

%otherwise, number of systems is given by the length of 'dimensions'
n=size(partition);
n=n(2);

if (nargin == 2)
    %if 'dimensions' are not given, assume qubits, i.e. always dimension 2
    dimensions = 2*ones(1,n);
else
    %throw an error if 'partion' and 'dimensions' have different length
    if (max(size(partition) ~= size(dimensions))==1)
        error('partition array and dimensions have different length');
    end
end

%throw an error if 'operator' is not a square matrix
opdims=size(operator);
if (opdims(1) ~= opdims(2))
    error('first argument is no square matrix');
end

%throw an error if 'operator' is not hermitian (within a certain precision)
if (max(max(abs(ctranspose(operator)-operator))) > 1e-12)
    error('first argument is not hermitian');
end

%throw an error if any value in 'dimensions' is smaller than two
if (min(dimensions)<2)
    error('dimensions must be larger than 1');
end

if (max(dimensions)==min(dimensions))
    %if all particles have the same dimension, we can simply use another
    %numerical system, e.g. the binary system. this speeds things up when
    %compared to the case in which the systems have different dimensions

    %obtain the dimensionality of each system (equals the first system's
    %dimensionsality, since all are the same)
    dim = dimensions(1);

    %throw an error if operator dimensions do not match the length of
    %'partition'
    if (opdims(1) ~= dim^n)
        error('operator dimensions do not match the partition array');
    end

    %******************* this is the main part *******************
    %start with the zero matrix to build up the partially transposed
    %operator
    oppt = zeros(dim^n,dim^n);

    %define the identity on the space of 'operator'
    id=eye(size(operator));

    for rowind=1:opdims(1) %loop through rows ...
        for colind=(rowind+1):opdims(2) %... and columns of the operator
                               %due to hermiticity and invariance of trace
                               %under partial transpose, only loop through
                               %upper right half
            col=dec2base(colind-1,dim,n); %determine current row and ...
            row=dec2base(rowind-1,dim,n); %column index in the base given by
                                      %'dimensions'

            %determine new row and column index by transposing the systems
            %indicated by 'partition'

            %for the new column index, take the ith digit from col if
            %partition(i) is zero. if partition(i) is one, take the ith
            %digit of row ...
            newcol=transpose((1-partition(:)).*str2num(col(:))+partition(:).*str2num(row(:)));
            %... and vice versa for the new row index
            newrow=transpose((1-partition(:)).*str2num(row(:))+partition(:).*str2num(col(:)));

            %note that newrow and newcol are row vectors. therefore,
            %convert them into strings and drop the white spaces
            newcol=strrep(num2str(newcol),' ','');
            newrow=strrep(num2str(newrow),' ','');

            %convert row and column index back into the decimal system
            newcolind=base2dec(newcol,dim)+1;
            newrowind=base2dec(newrow,dim)+1;

            %add the matrix element on its new place to oppt
            oppt=operator(rowind,colind)*(id(:,newrowind)*id(newcolind,:))+oppt;
        end
    end

    %add elements due to hermiticity
    oppt=oppt+ctranspose(oppt);
    %add diagonal elements (which did not change through the partial
    %transposition)
    oppt = oppt + diag(diag(operator));

else
    %if the system has different dimensions, the program is a bit more complex
    %and slower

    %throw an error if operator dimensions do not match the length of
    %'partition'
    if (opdims(1) ~= prod(dimensions))
        error('operator dimensions do not match the partition array');
    end

    %******************* this is the other main part *******************

    %start with the zero matrix to build up the partially transposed
    %operator
    oppt = zeros(size(operator));

    %define the identity on the space of 'operator'
    id=eye(size(operator));
    %define the n x n - identity
    idnxn=eye(n,n);

    %to loop through all matrix elements of 'operator', we need to create
    %the indices strings of the basis vectors in the usual notation. now,
    %however, the different digits run from zero to the corresponding
    %system's dimension (minus one), which differs from system to system.
    indices=zeros(n,1); %first index has only zeros

    for k=1:opdims(1) %loop through whole matrix 'operator'
        inddims=size(indices);
        last = indices(:,inddims(2)); %get the last index string

        for l=n:-1:1 %loop through digits of last index string
            if (last(l) < dimensions(l)-1) %the first digit from the right
                                           %which is still smaller than the
                                           %dimension (minus one) ...
                newvec=last+idnxn(:,l); %... must be increased by one ...
                newvec=newvec.*(vertcat(ones(l,1),zeros(n-l,1))); % ... and
                                       %all digits to the right set to zero
                indices=horzcat(indices,newvec); %append new index string
                break;
            end
        end
    end

    for rowind=1:opdims(1) %loop through rows ...
        for colind=(rowind+1):opdims(2) %... and columns of the operator
                               %due to hermiticity and invariance of trace
                               %under partial transpose, only loop through
                               %upper right half
            col=indices(:,colind); %write current row and ...
            row=indices(:,rowind); %column index as index string.
                                   %this time, row and col are column
                                   %vectors

            %determine new row and column index by transposing the systems
            %indicated by 'partition'

            %for the new column index, take the ith digit from col if
            %partition(i) is zero. if partition(i) is one, take the ith
            %digit of row ...
            newcol=transpose((1-partition(:)).*col(:)+partition(:).*row(:));
            %... and vice versa for the new row index
            newrow=transpose((1-partition(:)).*row(:)+partition(:).*col(:));

            %since newrow and newcol are row vectors denoting an index
            %string, we need to convert them back to a decimal number that
            %denotes the element's new position
            newcolind=find(ismember(transpose(indices),newcol,'rows'));
            newrowind=find(ismember(transpose(indices),newrow,'rows'));

            %add the matrix element on its new place to oppt
            oppt=operator(rowind,colind)*(id(:,newrowind)*id(newcolind,:))+oppt;
        end
    end

    %build the partial transpose of 'operator'
    %elements that used to be in the upper right half
    %oppt=full(sparse(colarray,rowarray,valarray,dim^n,dim^n,((dim^n)^2)/2-(dim^n)/2));
    %add elements due to hermiticity
    oppt=oppt+ctranspose(oppt);
    %add diagonal elements (which did not change through the partial
    %transposition)
    oppt = oppt + diag(diag(operator));
end

operator_pt=oppt;

author: Sable

repo: mcbench-benchmarks


text:
% Test file for fractional calculus.

function pass = test_fracInt(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

tol = 100*pref.chebfuneps;
dom = [0, 1];

% Test values:
xx = linspace(dom(1)+.1, dom(2)-.1, 10)';

%% Polynomials
x = chebfun('x', dom);
q = sqrt(2)/2;
k = 1;
for n = [1, 4]

    U = diff(x.^n, q);
    tru = gamma(n+1)./gamma(n+1-q)*chebfun(@(x) x.^(n-q), dom, 'exps', [n-q, 0]);

    err(k) = norm(tru(xx) - U(xx), inf);
    tol(k) = 10*eps*vscale(U)*hscale(U);
    k = k + 1;

end

%% Exponential
u = chebfun('exp(x)', dom);
trueC = chebfun('erf(sqrt(x)).*exp(x) + 1./sqrt(pi*x)', dom, 'exps', [-.5 0]);
trueRL = chebfun('erf(sqrt(x)).*exp(x)', dom, 'exps', [.5, 0]);

% RL
U = diff(u, .5, 'Caputo');
err(3) = norm(trueRL(xx) - U(xx), inf);
tol(3) = 10*eps*vscale(U)*hscale(U);

% Caputo
U = diff(u, .5, 'RL');
err(4) = norm(trueC(xx) - U(xx), inf);
tol(4) = 1e2*eps*vscale(U)*hscale(U);

%% Integrate twice:
xx = linspace(-sqrt(2)*pi+.1, pi-.1, 10).';
f = chebfun(@sin, [-sqrt(2)*pi, pi]);
F = cumsum(f);
G = fracInt(fracInt(f, .3), .7);
err(5) = norm(feval(F, xx) - feval(G, xx), inf);
tol(5) = 10*eps*vscale(G)*hscale(G);

%% Differentiate twice:
xx = linspace(-sqrt(2)*pi+.1, pi-.1, 10)';
f = chebfun(@sin, [-sqrt(2)*pi, pi]);
F = diff(f);
G = fracDiff(fracDiff(f, .3), .7);
err(6) = norm(feval(F, xx) - feval(G, xx), inf);
tol(6) = 1e2*eps*vscale(G)*hscale(G);

%% Test quasimatrix support
x = chebfun('x', [0 1]);
V = cheb2quasi(vander(x, 5));
F = diff(V, .5);
err(7) = 0;
for k = 1:numel(V)
    err(7) = err(7) + norm(diff(V(:,k), .5) - F(:,k));
end
tol(7) = 10*eps;

%%
pass = err < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nodes, edges] = createTestGraph03(varargin)
%CREATETESTGRAPH03  One-line description here, please.
%
%   output = createTestGraph03(input)
%
%   Example
%   [nodes edges] = createTestGraph03
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

nodes = [ ...
     30 10; ...
     10 30; ...
     80 30; ...
     30 50; ...
     60 50; ...
    100 50; ...
     10 70; ...
     80 70; ...
     30 90; ...
     60 90; ...
    100 90; ...
    130 90];

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


text:
function wathen_test11 ( )

%*****************************************************************************80
%
%% WATHEN_TEST11 assemble, factor and solve using WATHEN_ST + CG_ST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WATHEN_TEST11\n' );
  fprintf ( 1, '  Assemble, factor and solve a Wathen system\n' );
  fprintf ( 1, '  defined by WATHEN_ST and CG_ST.\n' );
  fprintf ( 1, '\n' );

  nx = 1;
  ny = 1;
  fprintf ( 1, '  Elements in X direction NX = %d\n', nx );
  fprintf ( 1, '  Elements in Y direction NY = %d\n', ny );
  fprintf ( 1, '  Number of elements = %d\n', nx * ny );
%
%  Compute the number of unknowns.
%
  n = wathen_order ( nx, ny );
  fprintf ( 1, '  Number of nodes N = %d\n', n );
%
%  Compute the matrix size.
%
  nz_num = wathen_st_size ( nx, ny );
  fprintf ( 1, '  Number of nonzeros = %d\n', nz_num );
%
%  Set up a random solution X1.
%
  seed = 123456789;
  [ x1, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Compute the matrix.
%
  seed = 123456789;
  [ row, col, a, seed ] = wathen_st ( nx, ny, nz_num, seed );
%
%  Compute the corresponding right hand side B.
%
  b = mv_st ( n, n, nz_num, row, col, a, x1 );
%
%  Solve the linear system.
%
  x2 = ones ( n, 1 );
  x2 = cg_st ( n, nz_num, row, col, a, b, x2 );
%
%  Compute the maximum solution error.
%
  e = max ( abs ( x1 - x2 ) );
  fprintf ( 1, '  Maximum solution error is %g\n', e );

  return
end


text:
%%
% MetaboRePort:

%% correct CheBI ID's in rBionet metabolite_structure
if 0
    load('C:\Users\0123322S\Desktop\X\met_strc_rBioNet_new_30_06_2022.mat')
    F = fieldnames(metabolite_structure_rBioNet);
    % reset chebiid
    for i = 1 : length(F)
        if strcmp(metabolite_structure_rBioNet.(F{i}).cheBIId,'1')
            metabolite_structure_rBioNet.(F{i}).cheBIId = NaN;
        end
    end
    % search hmdb for all chebi ids missing
    metabolite_structure_rBioNet2 = metabolite_structure_rBioNet;
    for i = 1 : length(F)
        if isnan(metabolite_structure_rBioNet2.(F{i}).cheBIId)
            [metabolite_structure_rBioNet2,IDsAddedHMDB,IDsMismatchHMDB] = parseHmdbWebPage(metabolite_structure_rBioNet2,i,i);
            i
            metabolite_structure_rBioNet2.(F{i}).cheBIId
        end
    end
    metabolite_structure_rBioNet = metabolite_structure_rBioNet2;
    save('C:\Users\0123322S\Desktop\X\met_strc_rBioNet_new_18_10_22.mat','metabolite_structure_rBioNet')
end

%% creating the files from scratch
load('C:\Users\0123322S\Desktop\X\met_strc_rBioNet_new_18_10_22.mat')
%folder = 'C:\Users\0123322S\Downloads\AGORA2_mat_files\AGORA2_mat_files\AGORA2_mat_files_1\';
folder = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\';
%folderUpdated = 'C:\Users\0123322S\Downloads\AGORA2_mat_files\AGORA2_mat_files\AGORA2_mat_files_1\';
folderUpdated = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\AGORA2_annotated\';
%reportDir = 'C:\Users\0123322S\Documents\GitHub\MetaboRePort\reports';
reportDir = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\AGORA2_annotated\metaboRePorts\';

dInfo=dir(folder);
modelList={dInfo.name};
modelList=modelList';
modelList(~contains(modelList(:,1),'.mat'),:)=[];

dInfoUpdated=dir(folderUpdated);
modelListUpdated={dInfoUpdated.name};
modelListUpdated=modelListUpdated';

dInfoReport=dir(reportDir);
modelListReport={dInfoReport.name};
modelListReport=modelListReport';

if 0
    tic;
    cnt = 1
    for i = s : e%length(modelList)
        i
        if isempty(find(ismember(modelListUpdated,modelList{i})))
            load(strcat(folder, modelList{i}));
            %[modelProp1,ScoresOverall1] = generateMemoteLikeScore(model);
            [modelUpdated] = populateModelMetStr(model, metabolite_structure_rBioNet,1);
            [modelUpdated] = annotateSBOTerms(modelUpdated);
            modelUpdated = rmfield(modelUpdated,'metInChIString'); % wrongly in microbe models

            [modelUpdated] = populateModelwithRxnIDs(modelUpdated);
            try
                [modelProp2,ScoresOverall2] = generateMetaboScore(modelUpdated);

                modelProperties.(regexprep(modelList{i},'.mat','')).ScoresOverall = ScoresOverall2;
                modelProperties.(regexprep(modelList{i},'.mat','')).modelUpdated = modelUpdated;
                modelProperties.(regexprep(modelList{i},'.mat','')).modelProp2 = modelProp2;
                ScoresOverall{i,1} = regexprep(modelList{i},'.mat','');
                ScoresOverall{i,2} = num2str(ScoresOverall2);
                % if mod(i,10)
                save(strcat(reportDir,'AGORA2_MetaboRePorts_',modelList{i}),'modelProperties','ScoresOverall');
                clear modelProperties
                %  end
                %% save updated mat file
            catch
                modelList{i}
                missing{cnt,1} = modelList{i};
                cnt = cnt + 1;
            end
            model =modelUpdated;
            save(strcat(folderUpdated,modelList{i}),'model');
            %     %%generate sbml file
            %     %remove description from model structure as this causes issues
            %
            %     modelUpdated = rmfield(modelUpdated,'description');
            %     fileName = regexprep(modelList{i},'.mat','');
            %     if isempty(find(ismember(modelListUpdated,strcat(fileName,'.xml'))))
            %         outmodel = writeCbModel(modelUpdated, 'format','sbml', 'fileName', strcat(folderUpdated,'',fileName));
            %     end
        end
    end
    toc;
    save('AGORA2_MetaboRePorts_ScoresOverall.mat','ScoresOverall');

end

%% generate missing metaboReports for existing updated annotated models
if 0
    tic;
    cnt = 1
    for i = s : e%length(modelList)
        i
        modelListReport = regexprep(modelListReport,'AGORA2_MetaboRePorts_','');
        if isempty(find(ismember(modelListReport,modelList{i})))
            load(strcat(folderUpdated, modelList{i}));
            try
                [modelProp2,ScoresOverall2] = generateMetaboScore(model);

                modelProperties.(regexprep(modelList{i},'.mat','')).ScoresOverall = ScoresOverall2;
                modelProperties.(regexprep(modelList{i},'.mat','')).modelUpdated = model;
                modelProperties.(regexprep(modelList{i},'.mat','')).modelProp2 = modelProp2;
                ScoresOverall{i,1} = regexprep(modelList{i},'.mat','');
                ScoresOverall{i,2} = num2str(ScoresOverall2);
                % if mod(i,10)
                save(strcat(reportDir,'AGORA2_MetaboRePorts_',modelList{i}),'modelProperties','ScoresOverall');
                clear modelProperties
                %  end
                %% save updated mat file
            catch
                modelList{i}
                missing{cnt,1} = modelList{i};
                cnt = cnt + 1;
            end

            %     %%generate sbml file
            %     %remove description from model structure as this causes issues
            %
            %     modelUpdated = rmfield(modelUpdated,'description');
            %     fileName = regexprep(modelList{i},'.mat','');
            %     if isempty(find(ismember(modelListUpdated,strcat(fileName,'.xml'))))
            %         outmodel = writeCbModel(modelUpdated, 'format','sbml', 'fileName', strcat(folderUpdated,'',fileName));
            %     end
        end
    end
    toc;
end

if 0
    % get all the report m files
    dInfoReport=dir(reportDir);
    modelListReport={dInfoReport.name};
    modelListReport=modelListReport';

    cnt = 1;
    for i = s :e% length(modelListReport)
        report = regexprep(modelListReport{i},'\.mat','\.html');
        report = regexprep(report,'AGORA2_MetaboRePorts_','modelreport_');
        i
        if ~isempty(strfind(modelListReport{i},'.mat')) && isempty(find(ismember(modelListReport,report)))
            load(strcat(reportDir,filesep,modelListReport{i}))
            orgName{1} = regexprep(modelListReport{i},'\.html','');
            orgName{1} = regexprep(orgName{1},'AGORA2_MetaboRePorts_','');
            evalc('generateMetaboReport(modelProperties,reportDir,orgName)');
            % get scores
            report
            F = fieldnames(modelProperties);
            ScoresAll(cnt,1) = modelProperties.(F{1}).modelProp2.Scores.Consistency;
            ScoresAll(cnt,2) = modelProperties.(F{1}).modelProp2.Scores.AnnotationMetabolites;
            ScoresAll(cnt,3) = modelProperties.(F{1}).modelProp2.Scores.AnnotationReactions;
            ScoresAll(cnt,4) = modelProperties.(F{1}).modelProp2.Scores.AnnotationGenes;
            ScoresAll(cnt,5) = modelProperties.(F{1}).modelProp2.Scores.AnnotationSBO;
            ScoresAll(cnt,6) = modelProperties.(F{1}).modelProp2.Scores.Overall;
            cnt = cnt +1;
        end
    end
end
%% updating only Chebii id's in the files
load('C:\Users\0123322S\Desktop\X\met_strc_rBioNet_new_18_10_22.mat')

folder = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\';
%folderUpdated = 'C:\Users\0123322S\Downloads\AGORA2_mat_files\AGORA2_mat_files\AGORA2_mat_files_1\';
folderUpdated = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\AGORA2_annotated\';
folderUpdatedNew = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\AGORA2_annotatedNew\';
folderUpdatedNewSBML = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\AGORA2_annotatedNewSBML\';

%reportDir = 'C:\Users\0123322S\Documents\GitHub\MetaboRePort\reports';
reportDir = 'C:\Users\0123322S\Dropbox\MY PAPERS\SUBMITTED\Submitted\AGORA2\NatureBiotech\resubmission\AGORA2_recreated\AGORA2_annotated\metaboRePorts\';

dInfo=dir(folder);
modelList={dInfo.name};
modelList=modelList';
modelList(~contains(modelList(:,1),'.mat'),:)=[];

dInfoUpdated=dir(folderUpdated);
modelListUpdated={dInfoUpdated.name};
modelListUpdated=modelListUpdated';

dInfoUpdatedSBML=dir(folderUpdatedNewSBML);
modelListUpdatedSBML={dInfoUpdatedSBML.name};
modelListUpdatedSBML=modelListUpdatedSBML';

dInfoReport=dir(reportDir);
modelListReport={dInfoReport.name};
modelListReport=modelListReport';

if 1
    tic;
    Fori = fieldnames(metabolite_structure_rBioNet);
    F = regexprep(Fori,'VMH_','');
    for i =  s: e%length(modelListUpdated)
        i

        dInfoUpdatedSBML=dir(folderUpdatedNewSBML);
        modelListUpdatedSBML={dInfoUpdatedSBML.name};
        modelListUpdatedSBML=modelListUpdatedSBML';
        modelName = regexprep(modelListUpdated{i},'.mat','');

        if ~isempty(strfind(strcat(folderUpdated, modelListUpdated{i}),'.mat')) && isempty(strmatch(strcat(modelName,'.xml'),modelListUpdatedSBML))
            i
            % load annotated model
            load(strcat(folderUpdated, modelListUpdated{i}));
            for k = 1 : length(model.mets)
                % replace metChEBIID
                met = split(model.mets{k},'[');
                if strcmp('sertrna(sec)',met{1}) % replace the abbr of this metabolite
                    met{1} = 'sertrna';
                end
                match = strmatch(met{1},F,'exact');
                if length(find(isnan(metabolite_structure_rBioNet.(Fori{match}).cheBIId)))==0 && ~isempty(metabolite_structure_rBioNet.(Fori{match}).cheBIId)
                    model.metChEBIID{k} = strcat('CHEBI:',metabolite_structure_rBioNet.(Fori{match}).cheBIId);
                end
            end
            % annotate model properties correctly

            model.modelName = regexprep(modelName,'_',' ');
            model.modelID = modelName;
            model.modelAnnotation = {model.modelID;...
                'version 2.0';...
                strcat('Created: ',model.description.date)
                'Molecular Systems Physiology Group, University of Galway, Ireland. https://thielelab.eu';...
                strcat('This is a metabolism reconstruction of:', model.modelName);...
                'The draft reconstruction was created using Kbase/ModelSeed and refined with the DEMETER pipeline from the COBRA toolbox v3';...
                'When using or adapting this reconstruction, please cite:';...
                'Heinken, A., Hertel, J., Acharya, G., Ravcheev, D.A., Nyga, M., Okpala, O.E., Hogan, M., Magnusdottir, S., Martinelli, F., Nap, B., Preciat, G., Edirisinghe, J.N., Henry, C.S., Fleming, R.M.T., Thiele, I., "AGORA2: Knowledge-driven genome-scale reconstruction of 7,302 human microbes for personalised medicine", Nat Biotech, accepted (2022), https://doi.org/10.1101/2020.11.09.375451.';...
                'Obtained from https://vmh.life/';...
                };
            model = rmfield(model,'description');
            model.description = modelName;

            save(strcat(folderUpdatedNew, modelListUpdated{i}),'model');

            %%generate sbml file
            %remove description from model structure as this causes issues

            if isempty(find(ismember(modelListUpdatedSBML,strcat(modelName,'.xml'))))
                outmodel = writeCbModel(model, 'format','sbml', 'fileName', strcat(folderUpdatedNewSBML,'',modelName));
            end
            clear outmodel model
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function varargout = xmkdir(varargin)
% VL_XMKDIR  Create a directory recursively.
%   VL_XMKDIR(PATH) creates all directory specified by PATH if they
%   do not exist (existing directories are skipped).
%
%   The function accepts the following options:
%
%   Pretend:: [false]
%     Set to true to avoid creating any directory but print which
%     directories would be created (implies 'Verbose',true).
%
%   Verbose:: [false]
%     Set to true to print the operations performed.
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_xmkdir(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
function disp(c)

   % CORRELATION/DISPLAY Command window display of a correlation object
   % See help correlation for fields

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   disp(' ');
   disp([inputname(1),' = '])
   disp(' ');
   %
   showline('WAVEFORMS', c.traces, 'vector');
   showline('TRIG', c.trig, 'vector');
   showline('CORR', c.corrmatrix, 'square matrix');
   showline('LAG', c.lags, 'square matrix');
   showline('STAT', c.stat, 'matrix');
   showline('LINK', c.link, 'matrix');
   showline('CLUST', c.clust, 'vector');
end

function showline(fname, value, desc)
   fprintf('%11s: %s %s\n',fname, getsizestr(value), desc);
end

function s = getsizestr(val)
   s = num2str(size(val),'%dx');
   s(s == ' ') = [];
   s = s(1:end-1);
end

author: geoscience-community-codes

repo: GISMO


text:
function E =extA(A, c)
%------------------------------------------------------------------------------
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 6, 1999.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
E=extU(extD(extR(extL(A, c), c), c), c);
%------------------------------------------------------------------------------


text:
function six_j_values_test ( )

%*****************************************************************************80
%
%% SIX_J_VALUES_TEST demonstrates the use of SIX_J_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIX_J_VALUES_TEST:\n' );
  fprintf ( 1, '  SIX_J_VALUES returns values of \n' );
  fprintf ( 1, '  the Wigner 6J coefficient.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '      J1      J2      J3      J4      J5      J6        SIX_J\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, j1, j2, j3, j4, j5, j6, fx ] = six_j_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %6.2f  %6.2f  %6.2f  %6.2f  %6.2f  %6.2f  %24.16f\n', ...
      j1, j2, j3, j4, j5, j6, fx );

  end

  return
end


text:
function idx = FindRoute(to)
% return the list of joint number connecting ROOT to 'to'
global uLINK

i = uLINK(to).mother;
%if i == 0
%    idx = [];    % search failed
%    fprintf('FindRoute: search failed\n');
%else
if i == 1
    idx = [to];
else
    idx = [FindRoute(i) to];
%     idx2 = FindRoute(i);
%     if isempty(idx2)
%         idx = [];
%     else
%         idx = [idx2 to];
%     end
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function  A = umADEIGSwhole(Ncell,p)

Dhat = zeros(p+1);
for n=0:p
    for m=n+1:2:p
        Dhat(n+1,m+1) =  (2*n+1);
    end
end

% M = int_{-1}^{1} L_n(x)*L_m(x) dx (note dx=1 assumed)
M = diag(2./(2*(0:p)+1));

% D_nm = int_{-1}^{1}  L_n(x)*L_m(x)*dx
D = M*Dhat;

% F_nm = L_n(-1)*L_m(-1)
F = (-1).^[0:p]'*(-1).^[0:p];

% G_nm = -L_n(-1)*L_m(1)
G = -(-1).^[0:p]'*ones(1,p+1);

% build entire matrix
Ntotal = Ncell*(p+1);
A    = zeros(Ntotal);

for cell = 1:Ncell
    id1 = ((cell-1)*(p+1)+1):cell*(p+1);
    A(id1,id1) = (1/2)*M\(D+F);

    if(cell~=1)
        id2 = id1-(p+1);
    else
        id2 = id1+(Ncell-1)*(p+1);
    end

    A(id1, id2) = (1/2)*M\G;
end

eA = eig(A);

plot(-real(eA), -imag(eA), 'bo');

repo: Aero-matlab


text:
classdef DepMat
    % DepMat A class used to update git repositories
    %
    %
    %
    %     Licence
    %     -------
    %     Part of DepMat. https://github.com/tomdoel/depmat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        RepoList
        RootSourceDir
        RepoUpdaterList
        RepoDirList
        RepoNameList
    end

    methods
        function obj = DepMat(repoList, rootSourceDir)
            obj.RepoList = repoList;
            obj.RootSourceDir = rootSourceDir;
            DepMat.fixCurlPath;

            obj.RepoDirList = cell(1, numel(obj.RepoList));
            obj.RepoNameList = cell(1, numel(obj.RepoList));
            obj.RepoUpdaterList = DepMatRepositoryUpdater.empty;

            for repoIndex = 1 : numel(obj.RepoList)
                repo = obj.RepoList(repoIndex);
                repoCombinedName = repo.FolderName;
                repoSourceDir = fullfile(obj.RootSourceDir, repoCombinedName);
                repo = DepMatRepositoryUpdater(repoSourceDir, repo);
                obj.RepoUpdaterList(repoIndex) = repo;
                obj.RepoDirList{repoIndex} = repoSourceDir;
                obj.RepoNameList{repoIndex} = repoCombinedName;
            end
        end

        function statusList = getAllStatus(obj)
            statusList = DepMatStatus.empty;
            for repoIndex = 1 : numel(obj.RepoList)
                repo = obj.RepoUpdaterList(repoIndex);
                statusList(repoIndex) = repo.getStatus;
            end
        end

        function success = updateAll(obj)
            success = true;
            for repoIndex = 1 : numel(obj.RepoList)
                repo = obj.RepoUpdaterList(repoIndex);
                success = success && repo.updateRepo;
            end
        end

        function anyChanged = cloneOrUpdateAll(obj)
            anyChanged = false;

            for repoIndex = 1 : numel(obj.RepoList)
                repo = obj.RepoUpdaterList(repoIndex);
                [~, changed] = repo.cloneOrUpdate;
                anyChanged = anyChanged || changed;
            end

        end
    end

    methods (Static)
        function [success, output] = execute(command)

            % Fixes issue #5. Apparent hang on launch if GitHub has not
            % been added to authorized_hosts. Ensure the `Are you sure you
            % want to continue connecting (yes/no)?` message is reported to
            % the user
            if ~ispc()
                command = [command ' > /dev/null'];
            end
            [return_value, output] = system(command, '-echo');

            success = return_value == 0;
            if ~success
                if strfind(output, 'Protocol https not supported or disabled in libcurl')
                    disp('! You need to modify the the DYLD_LIBRARY_PATH environment variable to point to a newer version of libcurl. The version installed with Matlab does not support using https with git.');
                end
            end
        end

        function installed = isGitInstalled
            if ispc
                command = 'where git';
            else
                command = 'which git';
            end

            installed = DepMat.execute(command);
        end

        function fixCurlPath
            % Matlab's curl configuration doesn't include https so git will not work.
            % We need to add the system curl configuration directory earlier in the
            % path so that it picks up this one instead of Matlab's

            try
                if ismac
                    pathName = 'DYLD_LIBRARY_PATH';
                    binDir = '/usr/lib';
                elseif isunix
                    pathName = 'LD_LIBRARY_PATH';
                    binDir = '/usr/lib';
                else
                    pathName = [];
                    binDir = [];
                end

                if ~isempty(pathName)
                    currentLibPath = getenv(pathName);
                    if (7 == exist(binDir, 'dir')) && ~strcmp(currentLibPath(1:length(binDir) + 1), [binDir ':'])
                        setenv(pathName, [binDir ':' currentLibPath]);
                    end
                end
            catch exception
                disp(['DepMat:fixCurlPath error: ' exception.message]);
            end
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function [E_Img, Par]   =  MCWNNM_ADMM1_NL_Denoising( N_Img, O_Img, Par )
E_Img           = N_Img;   % Estimated Image
[h, w, ch]  = size(E_Img);
Par.h = h;
Par.w = w;
Par.ch = ch;
Par = SearchNeighborIndex( Par );
% noisy image to patch
NoiPat =	Image2Patch( N_Img, Par );
Par.TolN = size(NoiPat, 2);
Sigma = ones(Par.ch, length(Par.SelfIndex));
for iter = 1 : Par.Iter
    Par.iter = iter;
    % iterative regularization
    E_Img =	E_Img + Par.delta * (N_Img - E_Img);
    % image to patch
    CurPat =	Image2Patch( E_Img, Par );
    % estimate local noise variance
    for c = 1:Par.ch
        if(iter == 1)
            TempSigma_arrCh = Par.lambda * Par.nSig0(c) * Sigma(c, :);
        else
            TempSigma_arrCh = Par.lambda * Sigma(c, :);
        end
        Sigma_arrCh((c-1)*Par.ps2+1:c*Par.ps2, :) = repmat(TempSigma_arrCh, [Par.ps2, 1]);
    end

    if (mod(iter-1, Par.Innerloop) == 0)
        Par.nlsp = Par.nlsp - 10;  % Lower Noise level, less NL patches
        NL_mat  =  Block_Matching_Real(CurPat, Par);% Caculate Non-local similar patches for each
    end
    [Y_hat, W_hat, Sigma]  =  MCWNNM_ADMM1_NL_Estimation( NL_mat, Sigma_arrCh, NoiPat, Par );   % Estimate all the patches
    E_Img = PGs2Image(Y_hat, W_hat, Par);
    %%
    PSNR  = csnr( O_Img, E_Img, 0, 0 );
    SSIM      =  cal_ssim( O_Img, E_Img, 0, 0 );
    fprintf( 'Iter = %2.3f, PSNR = %2.2f, SSIM = %2.2f \n', iter, PSNR, SSIM );
    Par.PSNR(iter, Par.image)  =   PSNR;
    Par.SSIM(iter, Par.image)      =  SSIM;
end
return;

author: csjunxu


text:
classdef TestScharr
    %TestScharr

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','blox.jpg');
    end

    methods (Static)
        function test_1
            img = imread(TestScharr.im);
            result = cv.Scharr(img);
            validateattributes(result, {class(img)}, {'size',size(img)});
        end

        function test_2
            img = imread(TestScharr.im);
            result = cv.Scharr(img, ...
                'XOrder',0, 'YOrder',1, 'Scale',1, 'Delta',0, ...
                'BorderType','Default', 'DDepth','double');
            validateattributes(result, {'double'}, {'size',size(img)});
        end

        function test_3
            img = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
            result = cv.Scharr(img);
            validateattributes(result, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.Scharr();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = psabr32(t,T,f,F,alpha,A,beta,nu,rho)
% Wu expansion of SABR density (2nd order)
tau = (T-t)/T;
if beta == 1    % account for log normal setting
    u = (log(f) - log(F))./(alpha * sqrt(T));
else
    u = (f^(1-beta)-F.^(1-beta))/(alpha*(1-beta)*sqrt(T));
end
v = log(alpha./A)./(nu*sqrt(T));
a11 = -beta * F^(beta-1).*A./nu.*(u-rho*v);
a10 = u^2*v-rho*u*v.^2;

a = 2*rho+beta*F^(beta-1)*A./nu;
b = 2+beta*(1-beta)*F^(2*(beta-1))*A.^2./nu^2;
c = beta*F^(beta-1)*A./nu;

a23 = rho^4*(20-6*b)-12*rho^3*a+rho^2*(3*a.^2-28+12*b) ...
    +12*a*rho+8-3*a.^2-6*b;
a22 = u^2*(3*a.^2-12*a*rho+6*b*(-1+rho^2)+2*rho^2+10)...
    - 2*u*v.*(rho^3*(2+3*b)+rho^2*(-9*a+3*c) ...
    +rho*(10+3*a.^2-3*b)-(3*a+3*c)) + v.^2.*((2+3*(a-2*rho).^2)*rho^2 ...
    + 6*c*rho*(-1+rho^2)-2);
a21 = u^4+v.^4*rho^2+u^3*v.*(8*rho-6*a)+u*v.^3.*(8*rho^3-6*a*rho^2) ...
    +u^2*v.^2.*(-14*rho^2+12*a*rho-4);
a20 = 3*u^4*v.^2-6*u^3*v.^3*rho+3*u^2*v.^4*rho^2;
y = 1./(nu*T*F^beta*A.^2).*(1+nu*sqrt(T)./(2*(-1+rho^2)).*(a11+a10/tau) ...
    +(nu^2*T)./(24*(1-rho^2)^2).*(a23*tau+a22+a21/tau+a20/tau^2))...
    ./(2*pi*tau*sqrt(1-rho^2)).*exp(-(u^2-2*rho*u*v+v.^2) ...
    ./(2*tau*(1-rho^2)));
end


text:
function rootPath=mrvTestRootPath()
%Determine path to root of the mrVista unit testing directory.
%
%        rootPath = mrvTestRootPath;
%
% This file and similarly named ones are used to identify the root
% directory of important repositories or files.
%
% This one is used as part of the testing of vistasoft
%
% This function MUST reside in the directory at the base of the vistatest
% directory structure.
%
% mrVista Team

rootPath = which('mrvTestRootPath');

rootPath = fileparts(rootPath);

return


text:
% Test code for the stroke utility procedures

% Fake character
stk1 = [1 2; 1 2; 1 2; 2 2];
stk2 = [1 2];
stk3 = [2 2];
stk4 = [0 0; 20 20];
stk1 = stk1 + randn(size(stk1));
stk2 = stk2 + randn(size(stk2));
stk3 = stk3 + randn(size(stk3));
stk4 = stk4 + randn(size(stk4));

vstk = [stk1; stk2; stk3; stk4];
S = {stk1; stk2; stk3; stk4};
S2 = {2*stk1; 2*stk2; 2*stk3; 2*stk4}; % rescale by 2
S_reverse = {stk1(end:-1:1,:); stk2; stk3; stk4(end:-1:1,:)};
com = mean(vstk);
rg = range(vstk,1);

% Check center of mass
com2 = com_char(S);
if aeq(com,com2)
   fprintf(1,'Center of mass test passed\n');
else
   error('COM failed');
end

% Check range of character
myrange = range_char(S);
if isequal(myrange,rg)
    fprintf(1,'Range test passed\n');
else
    error('Range test failed');
end

% Check removing strokes
minlen = 3;
mindist = 10;
RR = remove_short_stk(S,minlen,mindist);
if isequal(RR,S([1 4]))
   fprintf(1,'Remove stroke test passed\n');
else
   error('Remove stroke test failed');
end

% Check rescaling
R = apply_each_stroke(S,@(x)rescale_stk(x,[2 2]));
if isequal(R,S2)
   fprintf(1,'Apply each stroke test passed\n');
   fprintf(1,'Rescale test passed\n');
else
   error('Rescale test failed');
end

% Check reverse
RV = apply_each_stroke(S,@(x)reverse_stk(x));
if isequal(RV,S_reverse)
   fprintf(1,'Reverse test passed\n');
else
   error('Reverse test failed');
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
%% Housekeeping
close all
clc
%% load RISE before proceeding if you have not done so already

%% read the models and their calibrations
sw0=rise('usmodel','solve_linear',true);

%% solve the model
sw=solve(sw0,'solver','mn');

%% print results
sw.print_solution()

%% print solution for a subset of variables only
sw.print_solution({'a','b','c','cf','dc','dinve','dw','dy'})

%% compute regime-specific impulse responses
myirfs0=irf(sw,'irf_periods',20);

%% compute generalized impulse responses
myirfs1=irf(sw,'irf_periods',20,'irf_type','girf');

%% plot the impulse responses
close all
var_list={'dc','dinve','dw','dy','lab','pinf','r'};
figure('name','Impulse responses to a wage markup shock');
for ii=1:numel(var_list)
    subplot(3,3,ii)
	v=var_list{ii};
	aggregate=[myirfs0.ew.(v),myirfs1.ew.(v)];
    plot(aggregate,'linewidth',2);
    title(var_list{ii})
    if ii==1
        legend(aggregate.varnames)
    end
    axis tight
end

author: jmaih


text:
function invT = TransInv(T)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes a transformation matrix T.
% Returns its inverse.
% Uses the structure of transformation matrices to avoid taking a matrix
% inverse, for efficiency.
% Example Input:
%
% clear; clc;
% T = [[1, 0, 0, 0]; [0, 0, -1, 0]; [0, 1, 0, 3]; [0, 0, 0, 1]];
% invT = TransInv(T)
%
% Ouput:
% invT =
%     1     0     0     0
%     0     0     1    -3
%     0    -1     0     0
%     0     0     0     1

[R, p] = TransToRp(T);
invT = [R', -R' * p; 0, 0, 0, 1];
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function result = mse(A, B)
% MSE: mean-squared-error between two matrices/images.
%
% result=mse(A,B)
%

result = mean(mean(abs(A-B).^2));


text:
function DataSet = prtDataGenSwissRoll
% prtDataGenSwissRoll  Generates data from the Swiss Roll data set.
%
%   DataSet = prtDataGenSwissRoll generates a prtDataSetRegress from the
%   swiss roll data set. This data is drawn from a 2-D manifold embedded in
%   3-dimensions. For more information on ths data set see the following:
%
%   http://isomap.stanford.edu/code/Readme
%
%   Example:
%
%   ds = prtDataGenSwissRoll;
%
%   See also: prtDataSetClass, prtDataGenBiModal, prtDataGenIris,
%   prtDataGenManual, prtDataGenMary, prtDataGenNoisySinc,
%   prtDataGenOldFaithful,prtDataGenProtate, prtDataGenSprial,
%   prtDataGenSpiral3Regress, prtDataGenUnimodal, prtDataGenSwissRoll,
%   prtDataGenUnimodal, prtDataGenXor

% The 2-D manifold locations are in ds.getTargets, and the 3-dimensional
% embedding is in ds.getObservations.

swissRollFile = fullfile(prtRoot,']beta','dataGen','swissRoll','swiss_roll_data.mat');
swiss = load(swissRollFile);
X = swiss.X_data';
Y = swiss.Y_data';

DataSet = prtDataSetRegress(X,Y,'name','Standard Swiss Roll Data');

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
 function yn = jf_has_hct2(pn)
%function yn = jf_has_hct2(pn)
% determine if this matlab session has access to the binary "hct2" executable,
% for internal UM use only.

yn = 0;

if ~isunix, return, end

try
	% unix returns nonzero on failure
	if unix('which hct2 >& /dev/null')
		return
	end

	yn = 1;
	return

catch
	printm('the last error was "%s"', lasterr)
	return
end

error 'bug'

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function c = tapas_unitsq_sgm_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the unit square sigmoid observation model for binary responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The unit square sigmoid (ussgm) is the function
%
% f(x) = x^zeta/(x^zeta + (1-x)^zeta) = ustapas_sgm(x; zeta),
%
% where x is in the unit interval, and zeta > 0 is a parameter that determines the shape of the
% sigmoid. Since both its argument and value are always in the unit interval, its graph is
% restricted to the unit square, hence the name unit square sigmoid.
%
% In the application here, the ussgm is the probability of observing a decision y=1 (rather than
% the only alternative y=0) given the current probability mu1hat (or value) of input u=1:
%
% p(y=1|mu1hat) = ustapas_sgm(mu1hat; zeta)
%
% The parameter zeta regulates the steepness of the sigmoid such that it forms the diagonal of
% the unit square for zeta=1 and approaches a step function at 0.5 as zeta approaches infinity.
% Values of 0 < zeta < 1 lead to sigmoids with reverse concavity than usual, but they still
% represent valid observation models.
%
% Zeta can be interpreted as inverse decision noise. To have a shrinkage prior on this, choose a
% high value. It is estimated log-space since it has a natural lower bound at zero.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF 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/>.

% Config structure
c = struct;

% Is the decision based on predictions or posteriors? Comment as appropriate.
c.predorpost = 1; % Predictions
%c.predorpost = 2; % Posteriors

% Model name
c.model = 'tapas_unitsq_sgm';

% Sufficient statistics of Gaussian parameter priors

% Zeta
c.logzemu = log(48);
c.logzesa = 1;

% Gather prior settings in vectors
c.priormus = [
    c.logzemu,...
         ];

c.priorsas = [
    c.logzesa,...
         ];

% Model filehandle
c.obs_fun = @tapas_unitsq_sgm;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_unitsq_sgm_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
% Euclidean distance between polyhedra in 2D
% Section 8.2.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/09/05
% (a figure is generated)
%
% Given two polyhedra C = {x | A1*x <= b1} and D = {x | A2*x <= b2}, the
% distance between them is the optimal value of the problem:
%           minimize    || x - y ||_2
%               s.t.    A1*x <= b1
%                       A2*y <= b2
% Note: here x is in R^2

% Input data
randn('seed',0);
n = 2;
m = 2*n;
A1 = randn(m,n);
b1 = randn(m,1);
A2 = randn(m,n);
b2 = randn(m,1);

fprintf(1,'Computing the distance between the 2 polyhedra...');
% Solution via CVX
cvx_begin
    variables x(n) y(n)
    minimize (norm(x - y))
    norm(x,1) <= 2; %#ok
    norm(y-[4;3],inf) <=1; %#ok
cvx_end

fprintf(1,'Done! \n');

% Displaying results
disp('------------------------------------------------------------------');
disp('The distance between the 2 polyhedra C and D is: ' );
disp(['dist(C,D) = ' num2str(cvx_optval)]);
disp('The optimal points are: ')
disp('x = '); disp(x);
disp('y = '); disp(y);

%Plotting
figure;
fill([-2; 0; 2; 0],[0;2;0;-2],'b', [3;5;5;3],[2;2;4;4],'r')
axis([-3 6 -3 6])
axis square
hold on;
plot(x(1),x(2),'k.')
plot(y(1),y(2),'k.')
plot([x(1) y(1)],[x(2) y(2)])
title('Euclidean distance between 2 polyhedron in R^2');
xlabel('x_1');
ylabel('x_2');

author: yu-jiang

repo: radpbook

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


text:
function vw = switch2Vol(vw)
%
% Switches Volume view to volume mode:
% - changes vw.subdir to 'Volume'
% - sets vw.ui.grayVolButtons=2
% - computes/loads vw.coords
% - empties co, amp, ph, map
%
% djh, 7/98
% ress, 6/03: added UI controls for 3D buttons

if isempty(vw)
  myErrorDlg('Volume window must be open to switch it to volume mode.');
end
if ~strcmp(vw.viewType,'Volume') && ~strcmp(vw.viewType,'Gray')
  myErrorDlg('Only volume window can be switched to volume mode.');
end

% Set viewType
vw = viewSet(vw, 'viewType', 'Volume');
vw = viewSet(vw, 'subdir', 'Volume');

% Compute/load VOLUME.coords
vw = getVolCoords(vw);

% Empty data matrices
vw = viewSet(vw, 'co',  []);
vw = viewSet(vw, 'amp',  []);
vw = viewSet(vw, 'ph',  []);
vw = viewSet(vw, 'map',  []);

% Only do this for VOLUME, not for hiddenVolume
if strfind(viewGet(vw, 'name') ,'VOLUME')
  % Set vw.ui.grayVolButtons
  selectButton(vw.ui.grayVolButtons,2);
  % Hide the 3D buttons:
%   HideButtons(vw.ui.gray3dButtons);
%   vw = Disable3D(vw);
  % Load user preferences
  vw = loadPrefs(vw);
end

author: vistalab

repo: vistasoft


text:
% Minimize sidelobe level of a uniform linear array via spectral factorization
% "FIR Filter Design via Spectral Factorization and Convex Optimization" example
% by S.-P. Wu, S. Boyd, and L. Vandenberghe
% (figures are generated)
%
% Designs a uniform linear antenna array using spectral factorization method where:
% - it minimizes sidelobe level outside the beamwidth of the pattern
% - it has a constraint on the maximum ripple around unit gain in the beamwidth
%
%   minimize   max |y(theta)|                   for theta in the stop-beamwidth
%       s.t.   1/delta <= |y(theta)| <= delta   for theta in the pass-beamwidth
%
% We first replace the look-angle variable theta with the "frequency"
% variable omega, defined by omega = -2*pi*d/lambda*cos(theta).
% This transforms the antenna pattern y(theta) into a standard discrete
% Fourier transform of array weights w. Then we apply another change of
% variables: we replace w with its auto-correlation coefficients r.
%
% Now the problem can be solved via spectral factorization approach:
%
%   minimize   max R(omega)                        for omega in the stopband
%       s.t.   (1/delta)^2 <= R(omega) <= delta^2  for omega in the passband
%              R(omega) >= 0                       for all omega
%
% where R(omega) is the squared magnitude of the y(theta) array response
% (and the Fourier transform of the autocorrelation coefficients r).
% Variables are coefficients r. delta is the allowed passband ripple.
% This is a convex problem (can be formulated as an LP after sampling).
%
% Written for CVX by Almir Mutapcic 02/02/06

%********************************************************************
% problem specs: a uniform line array with inter-element spacing d
%                antenna element locations are at d*[0:n-1]
%                (the array pattern will be symmetric around origin)
%********************************************************************
n = 20;               % number of antenna elements
lambda = 1;           % wavelength
d = 0.45*lambda;      % inter-element spacing

% passband direction from 30 to 60 degrees (30 degrees bandwidth)
% transition band is 15 degrees on both sides of the passband
theta_pass = 40;
theta_stop = 50;

% passband max allowed ripple
ripple = 0.1; % in dB (+/- around the unit gain)

%********************************************************************
% construct optimization data
%********************************************************************
% number of frequency samples
m = 30*n;

% convert passband and stopband angles into omega frequencies
omega_zero = -2*pi*d/lambda;
omega_pass = -2*pi*d/lambda*cos(theta_pass*pi/180);
omega_stop = -2*pi*d/lambda*cos(theta_stop*pi/180);
omega_pi   = +2*pi*d/lambda;

% build matrix A that relates R(omega) and r, ie, R = A*r
omega = linspace(-pi,pi,m)';
A = exp( -j*omega(:)*[1-n:n-1] );

% passband constraint matrix
Ap = A(omega >= omega_zero & omega <= omega_pass,:);

% stopband constraint matrix
As = A(omega >= omega_stop & omega <= omega_pi,:);

%********************************************************************
% formulate and solve the magnitude design problem
%********************************************************************
cvx_begin
  variable r(2*n-1,1) complex
  % minimize stopband attenuation
  minimize( max( real( As*r ) ) )
  subject to
    % passband ripple constraints
    (10^(-ripple/20))^2 <= real( Ap*r ) <= (10^(+ripple/20))^2;
    % nonnegative-real constraint for all frequencies
    % a bit redundant: the passband frequencies are already constrained
    real( A*r ) >= 0;
    % auto-correlation symmetry constraints
    imag(r(n)) == 0;
    r(n-1:-1:1) == conj(r(n+1:end));
cvx_end

% check if problem was successfully solved
if ~strfind(cvx_status,'Solved')
    return
end

% find antenna weights by computing the spectral factorization
w = spectral_fact(r);

% divided by 2 since this is in PSD domain
min_sidelobe_level = 10*log10( cvx_optval );
fprintf(1,'The minimum sidelobe level is %3.2f dB.\n\n',...
          min_sidelobe_level);

%********************************************************************
% plots
%********************************************************************
% build matrix G that relates y(theta) and w, ie, y = G*w
theta = [-180:180]';
G = kron( cos(pi*theta/180), [0:n-1] );
G = exp(2*pi*i*d/lambda*G);
y = G*w;

% plot array pattern
figure(1), clf
ymin = -40; ymax = 5;
plot([-180:180], 20*log10(abs(y)), ...
     [theta_stop theta_stop],[ymin ymax],'r--',...
     [-theta_pass -theta_pass],[ymin ymax],'r--',...
     [-theta_stop -theta_stop],[ymin ymax],'r--',...
     [theta_pass theta_pass],[ymin ymax],'r--');
xlabel('look angle'), ylabel('mag y(theta) in dB');
axis([-180 180 ymin ymax]);

% polar plot
figure(2), clf
zerodB = 50;
dBY = 20*log10(abs(y)) + zerodB;
plot(dBY.*cos(pi*theta/180), dBY.*sin(pi*theta/180), '-');
axis([-zerodB zerodB -zerodB zerodB]), axis('off'), axis('square')
hold on
plot(zerodB*cos(pi*theta/180),zerodB*sin(pi*theta/180),'k:') % 0 dB
plot( (min_sidelobe_level + zerodB)*cos(pi*theta/180), ...
      (min_sidelobe_level + zerodB)*sin(pi*theta/180),'k:')  % min level
text(-zerodB,0,'0 dB')
text(-(min_sidelobe_level + zerodB),0,sprintf('%0.1f dB',min_sidelobe_level));
plot([0 60*cos(theta_pass*pi/180)], [0 60*sin(theta_pass*pi/180)], 'k:')
plot([0 60*cos(-theta_pass*pi/180)],[0 60*sin(-theta_pass*pi/180)],'k:')
plot([0 60*cos(theta_stop*pi/180)], [0 60*sin(theta_stop*pi/180)], 'k:')
plot([0 60*cos(-theta_stop*pi/180)],[0 60*sin(-theta_stop*pi/180)],'k:')
hold off

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ l, u ] = l1nn_lu ( n, h )

%*****************************************************************************80
%
%% L1NN_LU computes the LU factors of the 1D NN Laplacian.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real L(N,N), U(N,n), the LU factors.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1NN_LU - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1NN_LU - Fatal error!' );
  end

  l = zeros ( n, n );

  for i = 1 : n
    l(i,i) = 1.0;
  end

  for i = 2 : n
    l(i,i-1) = - 1.0;
  end

  u = zeros ( n, n );

  for i = 1 : n - 1
    u(i,i) = 1.0;
  end
  u(n,n) = 0.0;

  for i = 1 : n - 1
    u(i,i+1) = - 1.0;
  end

  u(1:n,1:n) = u(1:n,1:n) / h / h;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function title = p21_title ( )

%*****************************************************************************80
%
%% P21_TITLE returns the title of problem P21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 21, Enright and Pryce #E1';

  return
end


text:
function produceVASARIplot_LGG(pathResults,nameOutcomes,maxOrder,pathFig)
% -------------------------------------------------------------------------
% function produceVASARIplot_LGG(pathResults,nameOutcomes,maxOrder,pathFig)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function produce the Figure 4 of the LGG study.
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathResults: Full path to the directory containing all experiments.
%                     --> Ex: '/myProject/WORKSPACE/VASARI/RESULTS'
% 2. nameOutcomes: Cell of strings specifying the outcome names to analyze.
%                  --> Ex: {'nonIDH1','IDHcodel','progression','lowGrade'}
% 3. maxOrder: Integer specifying the maximal multivariable model order.
%              --> Ex: 10
% 4. pathFig: (optional).  Full path to where figure is saved without
%             displaying it. Put '' for displaying the figure and not
%             saving it to 'pathFig' (default).
%             --> Ex: ''
% -------------------------------------------------------------------------
% OUTPUTS: Final prediction models saved in a folder named
%          '/myProject/WORKSPACE/VASARI/FINAL_MODELS/nameOutcome/fSetName'
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathResults)

if nargin < 4
    pathFig = '';
end

signs = {'-r',':b','--g','-.y'};
nOutcomes = numel(nameOutcomes);

if isempty(pathFig)
    figure
else
    h = figure('visible','off');
end
for o = 1:nOutcomes
    nameOutcome = nameOutcomes{o};
    plotAllPredictionResults_VASARI_LGG(pathResults,nameOutcome,{'VASARI'},'AUC632',maxOrder,signs{o})
    hold on
end
legend(nameOutcomes)
title(['VASARI'],'FontSize',30,'FontWeight','bold')

if ~isempty(pathFig)
    cd(pathFig)
    saveas(h,['VASARI_PLOT'],'fig')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function self = despike(self, spiketype, maxRatio)
    % rsam.despike Despike a rsam object by comparing ratios of
    % concurrent samples. Checks for spikes lasting 1 or 2 samples.
    %   rsamdespiked = rsamobject.despike(spiketype, maxRatio)

    %   Example 1: Remove spikes which are at least 10 times
    %   adjacent samples. Store these in s.spikes.
    %        s = s.despike('spikes', 10)
    %
    %   Example 2: Remove spikes which are at least 3 times
    %   adjacent samples. Store these in s.events.
    %        s = s.despike('events', 3)
    %
    %   Inputs:
    %       maxRatio - Maximum ratio that defines "normal" data
    %                       compared to surrounding samples
    %   Outputs:
    %       s = rsam object with spikes removed

    % NOTE: THIS METHOD REQUIRES EXTRA PROPERTIES FOR rsam CLASS
        %spikes = []; % a vector of rsam objects that describe large spikes
        % in the data. Populated after running 'despike' method. These are
        % removed simultaneously from the data vector.
        %transientEvents = []; % a vector of rsam objects that describe
        % transient events in the data that might correspond to vt, rf, lp
        % etc. Populated after running 'despike' method with the
        % 'transientEvents' argument. These are not removed from the data
        % vector, but are instead returned in the continuousData vector.

    % find spikes lasting 1 sample only
    y= self.data;
    spikeNumber = 0;
    for i=2:length(self.data)-1
        if self.data(i)>maxRatio*self.data(i-1)
            if self.data(i)>maxRatio*self.data(i+1)
                %sample i is an outlier
                y(i) = mean([self.data(i-1) self.data(i+1)]);
                spikeNumber = spikeNumber + 1;
                %spikes(spikeNumber) = spike(self.dnum(i), self.data(i) - y(i), y(i), '');
                spikes(spikeNumber) = rsam('dnum', self.dnum(i), 'data', self.data(i) - y(i),  'sta', self.sta, 'chan', self.chan, 'seismogram_type', self.seismogram_type, 'units', self.units, 'measure', self.measure);
                disp(sprintf('%s: sample %d, time %s, before %f, this %f, after %f. Replacing with %f',upper(spiketype),i, datestr(self.dnum(i)), self.data(i-1), self.data(i), self.data(i+1), y(i)));
            end
        end
    end

    % find spikes lasting 2 samples
    for i=2:length(self.data)-2
        if self.data(i)>maxRatio*self.data(i-1) & self.data(i+1)>maxRatio*self.data(i-1)
            if self.data(i)>maxRatio*self.data(i+2) & self.data(i+1)>maxRatio*self.data(i+2)
                %samples i & i+1 are outliers
                y(i:i+1) = mean([self.data(i-1) self.data(i+2)]);
                spikeNumber = spikeNumber + 1;
                % spikes(spikeNumber) = spike( ...
                %     self.dnum(i:i+1), ...
                %     self.data(i:i+1) - y(i:i+1), ...
                %     y(i:i+1), '' );
                spikes(spikeNumber) = rsam('dnum', self.dnum(i:i+1), 'data', self.data(i:i+1) - y(i:i+1),  'sta', self.sta, 'chan', self.chan, 'seismogram_type', self.seismogram_type, 'units', self.units, 'measure', self.measure);
                disp(sprintf('%s: sample %d, time %s, before %f, these %f %f, after %f. Replacing with %f',upper(spiketype), i, datestr(self.dnum(i)), self.data(i-1), self.data(i), self.data(i+1), self.data(i+2), y(i)));
            end
        end
    end

    % find spikes lasting 3 samples - could be a short as 62
    % seconds
    if exist('spikes', 'var')
        if ~strcmp(spiketype, 'spikes') % only makes sense for events - an actual telemetry spike will be 1 or 2 samples long only (a few seconds)
            for i=2:length(self.data)-3
                if self.data(i)>maxRatio*self.data(i-1) & self.data(i+1)>maxRatio*self.data(i-1) & self.data(i+2)>maxRatio*self.data(i-1)
                    if self.data(i)>maxRatio*self.data(i+3) & self.data(i+1)>maxRatio*self.data(i+3) & self.data(i+2)>maxRatio*self.data(i+3)
                        %samples i & i+1 are outliers
                        y(i:i+2) = mean([self.data(i-1) self.data(i+3)]);
                        spikeNumber = spikeNumber + 1;
                        % spikes(spikeNumber) = spike( ...
                        %     self.dnum(i:i+1), ...
                        %     self.data(i:i+1) - y(i:i+1), ...
                        %     y(i:i+1), '' );
                        spikes(spikeNumber) = rsam('dnum', self.dnum(i:i+2), 'data', self.data(i:i+2) - y(i:i+2), 'sta', self.sta, 'chan', self.chan, 'seismogram_type', self.seismogram_type, 'units', self.units, 'measure', self.measure);
                        disp(sprintf('%s: sample %d, time %s, before %f, these %f %f %f, after %f. Replacing with %f',upper(spiketype), i, datestr(self.dnum(i)), self.data(i-1), self.data(i), self.data(i+1), self.data(i+2), self.data(i+3), y(i)));
                    end
                end
            end
        end
    end

    self.data = y;
    if exist('spikes', 'var')
        if strcmp(spiketype, 'spikes')
            self.spikes = spikes;
        else
            self.transientEvents = spikes;
        end
    end
end


text:
function r = csc(a)
%CSC          Taylor cosecans  csc(a)
%
%Thanks to George Corliss for providing the Taylor expansion
%

% written  06/03/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                   % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = getappdata(0,'INTLAB_TAYLOR_ORDER');

  ct = a.t;
  st = a.t;
  r = a;
  N = size(a.t,2);
  st(1,:) = sin(a.t(1,:));
  ct(1,:) = cos(a.t(1,:));
  r.t(1,:) = 1 ./ sin(a.t(1,:));
  st1 = -st(1,:);
  for j=2:K+1
    at_ = a.t(2:j,:);           % some 3 % faster
    st(j,:) = sum( repmat((1:j-1)',1,N).*ct(j-1:-1:1,:).*at_ , 1 ) ./ (j-1);
    if j~=K+1
      ct(j,:) = - sum( repmat((1:j-1)',1,N).*st(j-1:-1:1,:).*at_ , 1 ) ./ (j-1);
    end
    r.t(j,:) = sum( r.t(1:j-1,:).*st(j:-1:2,:) , 1 ) ./ st1;
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [bcx,bcy] = specific_flow(xbd,ybd)
%regcavity_flow   Reference problem 5.3 default inflow condition
%   [bcx,bcy] = specific_flow(xbd,ybd);
%   input
%          xbd          x coordinate vector
%          ybd          y coordinate vector
%
%   specifies regularized cavity flow boundary condition
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bcx=0*xbd; bcy=0*xbd;
k=find(ybd==1 & xbd>-1 & xbd<1); bcx(k)=(1-xbd(k).*xbd(k)).*(1+xbd(k).*xbd(k));
return


text:
function [ tokens ] = strtokenize( str, delim )
%STRTOKENIZE Summary of this function goes here
%   Detailed explanation goes here

    tokens = cell(0);
    remain = str;

    while true
        [token remain] = strtok(remain, delim);
        if (isempty(token))
            break;
        end;
        tokens(end+1) = {token};
    end;

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

n = length(theta);

% X.shape = [m, n], y.shape = [m, 1], theta.shape = [n, 1]
% y_pred = sigmoid(X*theta), y_pred.shape = [m, 1]
y_pred = sigmoid(X*theta);
J = -1.0 / m * sum( y.*log(y_pred) + (1-y).*log(1-y_pred) );
J = J + lambda / (2.0*m) * sum( theta(2:n, 1).^2 );  % ignore theta(1)
grad = 1.0 / m .* (X' * (y_pred - y));
grad(2:n, 1) = grad(2:n, 1) + lambda / m * theta(2:n, 1);   % ignore theta(1)

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

end

author: imLogM


text:
function Results = wrapper_FMINCON(P,X,Xlow,Xupp)

%Function Handles:
objfunc = @(x)Objective_Function(x,P);
nonlcon = @(x)NonLinCon(x,P);

%No linear constraints:
A=[]; Aeq=[]; B=[]; Beq=[];

%Store things in a Problem struct:
  Problem.objective = objfunc;
  Problem.x0 = convertVecStruct(X,P.structDat);
  Problem.Aineq = A;
  Problem.bineq = B;
  Problem.Aeq = Aeq;
  Problem.beq = Beq;
  Problem.lb = convertVecStruct(Xlow,P.structDat);
  Problem.ub = convertVecStruct(Xupp,P.structDat);
  Problem.nonlcon = nonlcon;
  Problem.options = P.options;
  Problem.solver = 'fmincon';

%Run fmincon:
tic
  [xSoln, Fval, ExitFlag, Output] = fmincon(Problem);
cpuTime = toc;

%Get the solution struct:
Xsoln = convertVecStruct(xSoln,P.structDat);

%Store the solver output:
name = 'fmincon';
input = makeStruct(name, Problem);
output = makeStruct(xSoln, Fval, ExitFlag, Output, cpuTime);
Results = makeStruct(P,input,output,Xsoln);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [out] = GLCM_Features4(glcmin,pairs)
%
%
% This is an update of GLCM_Features2 (vectorized) without ismember()
%
% GLCM_Features2 helps to calculate the features from the different GLCMs
% that are input to the function. The GLCMs are stored in a i x j x n
% matrix, where n is the number of GLCMs calculated usually due to the
% different orientation and displacements used in the algorithm. Usually
% the values i and j are equal to 'NumLevels' parameter of the GLCM
% computing function graycomatrix(). Note that matlab quantization values
% belong to the set {1,..., NumLevels} and not from {0,...,(NumLevels-1)}
% as provided in some references
% http://www.mathworks.com/access/helpdesk/help/toolbox/images/graycomatrix
% .html
%
% This vectorized version of GLCM_FEatures1.m reduces the 19 'for' loops
% used in the earlier code to 5 'for' loops
% http://blogs.mathworks.com/loren/2006/07/12/what-are-you-really-measuring
% /
% Using tic toc and cputime as in above discussion
%
% Although there is a function graycoprops() in Matlab Image Processing
% Toolbox that computes four parameters Contrast, Correlation, Energy,
% and Homogeneity. The paper by Haralick suggests a few more parameters
% that are also computed here. The code is not fully vectorized and hence
% is not an efficient implementation but it is easy to add new features
% based on the GLCM using this code. Takes care of 3 dimensional glcms
% (multiple glcms in a single 3D array)
%
% If you find that the values obtained are different from what you expect
% or if you think there is a different formula that needs to be used
% from the ones used in this code please let me know.
% A few questions which I have are listed in the link
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/239608
%
%
%
% Features computed
% Autocorrelation: [2]                      (out.autoc)
% Contrast: matlab/[1,2]                    (out.contr)
% Correlation: matlab                       (out.corrm)
% Correlation: [1,2]                        (out.corrp)
% Cluster Prominence: [2]                   (out.cprom)
% Cluster Shade: [2]                        (out.cshad)
% Dissimilarity: [2]                        (out.dissi)
% Energy: matlab / [1,2]                    (out.energ)
% Entropy: [2]                              (out.entro)
% Homogeneity: matlab                       (out.homom)
% Homogeneity: [2]                          (out.homop)
% Maximum probability: [2]                  (out.maxpr)
% Sum of squares: Variance [1]              (out.sosvh)
% Sum average [1]                           (out.savgh)
% Sum variance [1]                          (out.svarh)
% Sum entropy [1]                           (out.senth)
% Difference variance [1]                   (out.dvarh)
% Difference entropy [1]                    (out.denth)
% Information measure of correlation1 [1]   (out.inf1h)
% Informaiton measure of correlation2 [1]   (out.inf2h)
% Inverse difference (INV) is homom [3]     (out.homom)
% Inverse difference normalized (INN) [3]   (out.indnc)
% Inverse difference moment normalized [3]  (out.idmnc)
%
% The maximal correlation coefficient was not calculated due to
% computational instability
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
%
% Formulae from MATLAB site (some look different from
% the paper by Haralick but are equivalent and give same results)
% Example formulae:
% Contrast = sum_i(sum_j(  (i-j)^2 * p(i,j) ) ) (same in matlab/paper)
% Correlation = sum_i( sum_j( (i - u_i)(j - u_j)p(i,j)/(s_i.s_j) ) ) (m)
% Correlation = sum_i( sum_j( ((ij)p(i,j) - u_x.u_y) / (s_x.s_y) ) ) (p[2])
% Energy = sum_i( sum_j( p(i,j)^2 ) )           (same in matlab/paper)
% Homogeneity = sum_i( sum_j( p(i,j) / (1 + |i-j|) ) ) (as in matlab)
% Homogeneity = sum_i( sum_j( p(i,j) / (1 + (i-j)^2) ) ) (as in paper)
%
% Where:
% u_i = u_x = sum_i( sum_j( i.p(i,j) ) ) (in paper [2])
% u_j = u_y = sum_i( sum_j( j.p(i,j) ) ) (in paper [2])
% s_i = s_x = sum_i( sum_j( (i - u_x)^2.p(i,j) ) ) (in paper [2])
% s_j = s_y = sum_i( sum_j( (j - u_y)^2.p(i,j) ) ) (in paper [2])
%
%
% Normalize the glcm:
% Compute the sum of all the values in each glcm in the array and divide
% each element by it sum
%
% Haralick uses 'Symmetric' = true in computing the glcm
% There is no Symmetric flag in the Matlab version I use hence
% I add the diagonally opposite pairs to obtain the Haralick glcm
% Here it is assumed that the diagonally opposite orientations are paired
% one after the other in the matrix
% If the above assumption is true with respect to the input glcm then
% setting the flag 'pairs' to 1 will compute the final glcms that would result
% by setting 'Symmetric' to true. If your glcm is computed using the
% Matlab version with 'Symmetric' flag you can set the flag 'pairs' to 0
%
% References:
% 1. R. M. Haralick, K. Shanmugam, and I. Dinstein, Textural Features of
% Image Classification, IEEE Transactions on Systems, Man and Cybernetics,
% vol. SMC-3, no. 6, Nov. 1973
% 2. L. Soh and C. Tsatsoulis, Texture Analysis of SAR Sea Ice Imagery
% Using Gray Level Co-Occurrence Matrices, IEEE Transactions on Geoscience
% and Remote Sensing, vol. 37, no. 2, March 1999.
% 3. D A. Clausi, An analysis of co-occurrence texture statistics as a
% function of grey level quantization, Can. J. Remote Sensing, vol. 28, no.
% 1, pp. 45-62, 2002
% 4. http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
%
%
% Example:
%
% Usage is similar to graycoprops() but needs extra parameter 'pairs' apart
% from the GLCM as input
% I = imread('circuit.tif');
% GLCM2 = graycomatrix(I,'Offset',[2 0;0 2]);
% stats = GLCM_features4(GLCM2,0)
% The output is a structure containing all the parameters for the different
% GLCMs
%
% [Avinash Uppuluri: avinash_uv@yahoo.com: Last modified: 04/05/2010]

% If 'pairs' not entered: set pairs to 0
if ((nargin > 2) || (nargin == 0))
   error('Too many or too few input arguments. Enter GLCM and pairs.');
elseif ( (nargin == 2) )
    if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1))
       error('The GLCM should be a 2-D or 3-D matrix.');
    elseif ( size(glcmin,1) ~= size(glcmin,2) )
        error('Each GLCM should be square with NumLevels rows and NumLevels cols');
    end
elseif (nargin == 1) % only GLCM is entered
    pairs = 0; % default is numbers and input 1 for percentage
    if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1))
       error('The GLCM should be a 2-D or 3-D matrix.');
    elseif ( size(glcmin,1) ~= size(glcmin,2) )
       error('Each GLCM should be square with NumLevels rows and NumLevels cols');
    end
end

format long e
if (pairs == 1)
    newn = 1;
    for nglcm = 1:2:size(glcmin,3)
        glcm(:,:,newn)  = glcmin(:,:,nglcm) + glcmin(:,:,nglcm+1);
        newn = newn + 1;
    end
elseif (pairs == 0)
    glcm = glcmin;
end

size_glcm_1 = size(glcm,1);
size_glcm_2 = size(glcm,2);
size_glcm_3 = size(glcm,3);

% checked
out.autoc = zeros(1,size_glcm_3); % Autocorrelation: [2]
out.contr = zeros(1,size_glcm_3); % Contrast: matlab/[1,2]
out.corrm = zeros(1,size_glcm_3); % Correlation: matlab
out.corrp = zeros(1,size_glcm_3); % Correlation: [1,2]
out.cprom = zeros(1,size_glcm_3); % Cluster Prominence: [2]
out.cshad = zeros(1,size_glcm_3); % Cluster Shade: [2]
out.dissi = zeros(1,size_glcm_3); % Dissimilarity: [2]
out.energ = zeros(1,size_glcm_3); % Energy: matlab / [1,2]
out.entro = zeros(1,size_glcm_3); % Entropy: [2]
out.homom = zeros(1,size_glcm_3); % Homogeneity: matlab
out.homop = zeros(1,size_glcm_3); % Homogeneity: [2]
out.maxpr = zeros(1,size_glcm_3); % Maximum probability: [2]

out.sosvh = zeros(1,size_glcm_3); % Sum of sqaures: Variance [1]
out.savgh = zeros(1,size_glcm_3); % Sum average [1]
out.svarh = zeros(1,size_glcm_3); % Sum variance [1]
out.senth = zeros(1,size_glcm_3); % Sum entropy [1]
out.dvarh = zeros(1,size_glcm_3); % Difference variance [4]
%out.dvarh2 = zeros(1,size_glcm_3); % Difference variance [1]
out.denth = zeros(1,size_glcm_3); % Difference entropy [1]
out.inf1h = zeros(1,size_glcm_3); % Information measure of correlation1 [1]
out.inf2h = zeros(1,size_glcm_3); % Informaiton measure of correlation2 [1]
%out.mxcch = zeros(1,size_glcm_3);% maximal correlation coefficient [1]
%out.invdc = zeros(1,size_glcm_3);% Inverse difference (INV) is homom [3]
out.indnc = zeros(1,size_glcm_3); % Inverse difference normalized (INN) [3]
out.idmnc = zeros(1,size_glcm_3); % Inverse difference moment normalized [3]

glcm_sum  = zeros(size_glcm_3,1);
glcm_mean = zeros(size_glcm_3,1);
glcm_var  = zeros(size_glcm_3,1);

% http://www.fp.ucalgary.ca/mhallbey/glcm_mean.htm confuses the range of
% i and j used in calculating the means and standard deviations.
% As of now I am not sure if the range of i and j should be [1:Ng] or
% [0:Ng-1]. I am working on obtaining the values of mean and std that get
% the values of correlation that are provided by matlab.
u_x = zeros(size_glcm_3,1);
u_y = zeros(size_glcm_3,1);
s_x = zeros(size_glcm_3,1);
s_y = zeros(size_glcm_3,1);

% checked p_x p_y p_xplusy p_xminusy
p_x = zeros(size_glcm_1,size_glcm_3); % Ng x #glcms[1]
p_y = zeros(size_glcm_2,size_glcm_3); % Ng x #glcms[1]
p_xplusy = zeros((size_glcm_1*2 - 1),size_glcm_3); %[1]
p_xminusy = zeros((size_glcm_1),size_glcm_3); %[1]
% checked hxy hxy1 hxy2 hx hy
hxy  = zeros(size_glcm_3,1);
hxy1 = zeros(size_glcm_3,1);
hx   = zeros(size_glcm_3,1);
hy   = zeros(size_glcm_3,1);
hxy2 = zeros(size_glcm_3,1);

corm = zeros(size_glcm_3,1);
corp = zeros(size_glcm_3,1);

for k = 1:size_glcm_3

    glcm_sum(k) = sum(sum(glcm(:,:,k)));
    glcm(:,:,k) = glcm(:,:,k)./glcm_sum(k); % Normalize each glcm
    glcm_mean(k) = mean2(glcm(:,:,k)); % compute mean after norm
    glcm_var(k)  = (std2(glcm(:,:,k)))^2;

    for i = 1:size_glcm_1

        for j = 1:size_glcm_2
            p_x(i,k) = p_x(i,k) + glcm(i,j,k);
            p_y(i,k) = p_y(i,k) + glcm(j,i,k); % taking i for j and j for i
            %if (ismember((i + j),[2:2*size_glcm_1]))
                p_xplusy((i+j)-1,k) = p_xplusy((i+j)-1,k) + glcm(i,j,k);
            %end
            %if (ismember(abs(i-j),[0:(size_glcm_1-1)]))
                p_xminusy((abs(i-j))+1,k) = p_xminusy((abs(i-j))+1,k) +...
                    glcm(i,j,k);
            %end
        end
    end

end

% marginal probabilities are now available [1]
% p_xminusy has +1 in index for matlab (no 0 index)
% computing sum average, sum variance and sum entropy:

%Q    = zeros(size(glcm));

i_matrix  = repmat([1:size_glcm_1]',1,size_glcm_2);
j_matrix  = repmat([1:size_glcm_2],size_glcm_1,1);
% i_index = [ 1 1 1 1 1 .... 2 2 2 2 2 ... ]
i_index   = j_matrix(:);
% j_index = [ 1 2 3 4 5 .... 1 2 3 4 5 ... ]
j_index   = i_matrix(:);
xplusy_index = [1:(2*(size_glcm_1)-1)]';
xminusy_index = [0:(size_glcm_1-1)]';
mul_contr = abs(i_matrix - j_matrix).^2;
mul_dissi = abs(i_matrix - j_matrix);
%div_homop = ( 1 + mul_contr); % used from the above two formulae
%div_homom = ( 1 + mul_dissi);

for k = 1:size_glcm_3 % number glcms

    out.contr(k) = sum(sum(mul_contr.*glcm(:,:,k)));
    out.dissi(k) = sum(sum(mul_dissi.*glcm(:,:,k)));
    out.energ(k) = sum(sum(glcm(:,:,k).^2));
    out.entro(k) = - sum(sum((glcm(:,:,k).*log(glcm(:,:,k) + eps))));
    out.homom(k) = sum(sum((glcm(:,:,k)./( 1 + mul_dissi))));
    out.homop(k) = sum(sum((glcm(:,:,k)./( 1 + mul_contr))));
    % [1] explains sum of squares variance with a mean value;
    % the exact definition for mean has not been provided in
    % the reference: I use the mean of the entire normalized glcm
    out.sosvh(k) = sum(sum(glcm(:,:,k).*((i_matrix - glcm_mean(k)).^2)));
    out.indnc(k) = sum(sum(glcm(:,:,k)./( 1 + (mul_dissi./size_glcm_1) )));
    out.idmnc(k) = sum(sum(glcm(:,:,k)./( 1 + (mul_contr./(size_glcm_1^2)))));
    out.maxpr(k) = max(max(glcm(:,:,k)));

    u_x(k)       = sum(sum(i_matrix.*glcm(:,:,k)));
    u_y(k)       = sum(sum(j_matrix.*glcm(:,:,k)));
    % using http://www.fp.ucalgary.ca/mhallbey/glcm_variance.htm for s_x
    % s_y : This solves the difference in value of correlation and might be
    % the right value of standard deviations required
    % According to this website there is a typo in [2] which provides
    % values of variance instead of the standard deviation hence a square
    % root is required as done below:
    s_x(k)  = (sum(sum( ((i_matrix - u_x(k)).^2).*glcm(:,:,k) )))^0.5;
    s_y(k)  = (sum(sum( ((j_matrix - u_y(k)).^2).*glcm(:,:,k) )))^0.5;

   corp(k) = sum(sum((i_matrix.*j_matrix.*glcm(:,:,k))));
   corm(k) = sum(sum(((i_matrix - u_x(k)).*(j_matrix - u_y(k)).*glcm(:,:,k))));

   out.autoc(k) = corp(k);
   out.corrp(k) = (corp(k) - u_x(k)*u_y(k))/(s_x(k)*s_y(k));
   out.corrm(k) = corm(k) / (s_x(k)*s_y(k));

   out.cprom(k) = sum(sum(((i_matrix + j_matrix - u_x(k) - u_y(k)).^4).*...
                glcm(:,:,k)));
   out.cshad(k) = sum(sum(((i_matrix + j_matrix - u_x(k) - u_y(k)).^3).*...
                glcm(:,:,k)));

   out.savgh(k) = sum((xplusy_index + 1).*p_xplusy(:,k));
   % the summation for savgh is for i from 2 to 2*Ng hence (i+1)
   out.senth(k) =  - sum(p_xplusy(:,k).*...
            log(p_xplusy(:,k) + eps));

    % compute sum variance with the help of sum entropy
    out.svarh(k) = sum((((xplusy_index + 1) - out.senth(k)).^2).*...
        p_xplusy(:,k));
        % the summation for savgh is for i from 2 to 2*Ng hence (i+1)

    % compute difference variance, difference entropy,
    % out.dvarh2(k) = var(p_xminusy(:,k));
    % but using the formula in
    % http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
    % we have for dvarh
    out.denth(k) = - sum((p_xminusy(:,k)).*...
        log(p_xminusy(:,k) + eps));
    out.dvarh(k) = sum((xminusy_index.^2).*p_xminusy(:,k));

    % compute information measure of correlation(1,2) [1]
    hxy(k) = out.entro(k);
    glcmk  = glcm(:,:,k)';
    glcmkv = glcmk(:);

    hxy1(k) =  - sum(glcmkv.*log(p_x(i_index,k).*p_y(j_index,k) + eps));
    hxy2(k) =  - sum(p_x(i_index,k).*p_y(j_index,k).*...
        log(p_x(i_index,k).*p_y(j_index,k) + eps));
     hx(k) = - sum(p_x(:,k).*log(p_x(:,k) + eps));
     hy(k) = - sum(p_y(:,k).*log(p_y(:,k) + eps));

    out.inf1h(k) = ( hxy(k) - hxy1(k) ) / ( max([hx(k),hy(k)]) );
    out.inf2h(k) = ( 1 - exp( -2*( hxy2(k) - hxy(k) ) ) )^0.5;

    %     eig_Q(k,:)   = eig(Q(:,:,k));
    %     sort_eig(k,:)= sort(eig_Q(k,:),'descend');
    %     out.mxcch(k) = sort_eig(k,2)^0.5;
    % The maximal correlation coefficient was not calculated due to
    % computational instability
    % http://murphylab.web.cmu.edu/publications/boland/boland_node26.html

end

%       GLCM Features (Soh, 1999; Haralick, 1973; Clausi 2002)
%           f1. Uniformity / Energy / Angular Second Moment (done)
%           f2. Entropy (done)
%           f3. Dissimilarity (done)
%           f4. Contrast / Inertia (done)
%           f5. Inverse difference
%           f6. correlation
%           f7. Homogeneity / Inverse difference moment
%           f8. Autocorrelation
%           f9. Cluster Shade
%          f10. Cluster Prominence
%          f11. Maximum probability
%          f12. Sum of Squares
%          f13. Sum Average
%          f14. Sum Variance
%          f15. Sum Entropy
%          f16. Difference variance
%          f17. Difference entropy
%          f18. Information measures of correlation (1)
%          f19. Information measures of correlation (2)
%          f20. Maximal correlation coefficient
%          f21. Inverse difference normalized (INN)
%          f22. Inverse difference moment normalized (IDN)

author: Sable

repo: mcbench-benchmarks