text:
function [points, pos, faceInds, lineInds] = intersectLineMesh3d(line, vertices, varargin)
%INTERSECTLINEMESH3D Intersection points of a 3D line with a mesh.
%
%   INTERS = intersectLineMesh3d(LINE, VERTICES, FACES)
%   Compute the intersection points between a 3D line and a 3D mesh defined
%   by vertices and faces. The mesh data is provided as a pair of arrays,
%   with VERTICES being a NV-by-3 array of vertex coordinates, and FACES
%   being a NF-by-3 or NF-by-4 array of face vertex indices.
%   The LINE data correspond to a 1-by-6 row vector concatenating the line
%   origin and direction. LINE can also be a NL-by-6 array representing a
%   collection of lines with various origins and directions.
%
%   INTERS = intersectLineMesh3d(LINE, MESH)
%   Provides the mesh data as a struct with the fields 'vertices' and
%   'faces'.
%
%   [INTERS, POS, IFACES] = intersectLineMesh3d(...)
%   Also returns the position of each intersection point on the input line,
%   and the index of the intersected faces.
%   If POS > 0, the point is also on the ray corresponding to the line.
%
%   [INTERS, POS, IFACES, ILINES] = intersectLineMesh3d(...)
%   Also returns the index of the line each intersection point belong to.
%
%   Example
%     [V, F] = createCube;
%     line = [.2 .3 .4 1 0 0];
%     pts = intersectLineMesh3d(line, V, F)
%     pts =
%         1.0000    0.3000    0.4000
%              0    0.3000    0.4000
%
%   See also
%   meshes3d, triangulateFaces, intersectLineTriangle3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2011-12-20, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% tolerance for detecting if a point is on line or within edge bounds
tol = 1e-12;

% parsing
if ~isempty(varargin)
    if isscalar(varargin{1})
        tol = varargin{1};
        [vertices, faces] = parseMeshData(vertices);
    else
        faces = varargin{1};
        varargin(1) = [];
        if ~isempty(varargin)
            tol = varargin{1};
        end
    end
else
    [vertices, faces] = parseMeshData(vertices);
end

% ensure the mesh has triangular faces
tri2Face = [];
if iscell(faces) || size(faces, 2) ~= 3
    [faces, tri2Face] = triangulateFaces(faces);
end

% find triangle edge vectors
t0  = vertices(faces(:,1), :);
u   = vertices(faces(:,2), :) - t0;
v   = vertices(faces(:,3), :) - t0;

% triangle normal
n   = crossProduct3d(u, v);

% direction vectors of lines and origins of lines
dv = permute(line(:,4:6), [3 2 1]);
d0 = permute(line(:,1:3), [3 2 1]);

% vector between triangle origin and line origin
w0 = d0 - t0;

a = -sum(n .* w0, 2); % negative dot product
b = sum(n .* dv, 2);  % dot product

valid = abs(b) > tol & vectorNorm3d(n) > tol;

% compute intersection point of line with supporting plane
% If pos < 0: point before ray
% IF pos > |dir|: point after edge
pos = a ./ b;

% coordinates of intersection point
points = d0 + (pos .* dv);

%% test if intersection point is inside triangle

% normalize direction vectors of triangle edges
uu  = dot(u, u, 2);
uv  = dot(u, v, 2);
vv  = dot(v, v, 2);

% coordinates of vector v in triangle basis
w   = points - t0;
wu  = sum(w .* u, 2);
wv  = sum(w .* v, 2);

% normalization constant
D = uv.^2 - uu .* vv;

% test first coordinate
s = (uv .* wv - vv .* wu) ./ D;
ind1 = s < -tol | s > (1.0 + tol);

% test second coordinate, and third triangle edge
t = (uv .* wu - uu .* wv) ./ D;
ind2 = t < -tol | (s + t) > (1.0 + tol);

% keep only interesting points
inds = ~ind1 & ~ind2 & valid;
[faceInds, lineInds] = find(permute(inds, [1 3 2]));

% Bit of an indexing trick to get points in appropriate order
points = points(sub2ind(size(points), ...
    faceInds+[0 0 0], faceInds*0+(1:3), lineInds+[0 0 0]) );

if nargout > 1
    pos = pos(inds);

    % convert to face indices of original mesh
    if ~isempty(tri2Face)
        faceInds = tri2Face(faceInds);
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function value = diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 )

%*****************************************************************************80
%
%% DIAEDG chooses a diagonal edge.
%
%  Discussion:
%
%    The routine determines whether 0--2 or 1--3 is the diagonal edge
%    that should be chosen, based on the circumcircle criterion, where
%    (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple
%    quadrilateral in counterclockwise order.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, real X0, Y0, X1, Y1, X2, Y2, X3, Y3, the
%    coordinates of the vertices of a quadrilateral, given in
%    counter clockwise order.
%
%    Output, integer VALUE, chooses a diagonal:
%    +1, if diagonal edge 02 is chosen;
%    -1, if diagonal edge 13 is chosen;
%     0, if the four vertices are cocircular.
%
  tol = 100.0 * r8_epsilon ( );

  dx10 = x1 - x0;
  dy10 = y1 - y0;
  dx12 = x1 - x2;
  dy12 = y1 - y2;
  dx30 = x3 - x0;
  dy30 = y3 - y0;
  dx32 = x3 - x2;
  dy32 = y3 - y2;

  tola = tol * max ( abs ( dx10 ), ...
               max ( abs ( dy10 ), ...
               max ( abs ( dx30 ), abs ( dy30 ) )));

  tolb = tol * max ( abs ( dx12 ), ...
               max ( abs ( dy12 ), ...
               max ( abs ( dx32 ), abs ( dy32 ) )));

  ca = dx10 * dx30 + dy10 * dy30;
  cb = dx12 * dx32 + dy12 * dy32;

  if ( tola < ca & tolb < cb )

    value = -1;

  elseif ( ca < -tola & cb < -tolb )

    value = 1;

  else

    tola = max ( tola, tolb );
    s = ( dx10 * dy30 - dx30 * dy10 ) * cb ...
      + ( dx32 * dy12 - dx12 * dy32 ) * ca;

    if ( tola < s )
      value = -1;
    elseif ( s < -tola )
      value = 1;
    else
      value = 0;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function CoreCompileMexFiles(mex_cache, output_directory, mex_files_to_compile, force_recompile, retry_instructions, reporting)
    % CoreCompileMexFiles. Checks if mex files are up to date and re-compiles if
    % necessary
    %
    %     CoreCompileMexFiles takes in a list of mex files to be compiled,
    %     with version numbers. The last compiled versions
    %     are cached using the supplied mex_cache. If the versions have changed, files
    %     are automatically recompiled. Recompilation also occurs if the output
    %     files are missing or the mex_cache file has been deleted.
    %
    %     CoreCompileMexFiles checks for the presence of compilers and warns the
    %     user if no compiler has been detected. CoreCompileMexFiles does not
    %     attempt to recompile a particular mex file again after a failed
    %     compile attempt. Recompilation can be achieved using the
    %     force_recompile flag or by deleting the mex_cache file.
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    cached_mex_file_info = mex_cache.MexInfoMap;
    framework_cache_was_missing = mex_cache.IsNewlyCreated;

    compiler = GetNameOfCppCompiler;
    if isempty(compiler)
        compiler = MexSetup(retry_instructions, reporting);
    end
    cuda_compiler = GetCudaCompiler;

    if ~isempty(compiler)
        CheckMexFiles(mex_files_to_compile, cached_mex_file_info, output_directory, framework_cache_was_missing, reporting);
        Compile(mex_files_to_compile, mex_cache, cached_mex_file_info, output_directory, compiler, cuda_compiler, force_recompile, retry_instructions, reporting);
    end

    mex_cache.UpdateCache(mex_files_to_compile, reporting);
end

function compiler = MexSetup(retry_instructions, reporting)
    cpp_compilers = mex.getCompilerConfigurations('C++', 'Installed');
    if isempty(cpp_compilers)
        compiler = [];
        reporting.ShowWarning('CoreCompileMexFiles:NoCompiler', ['I cannot compile mex files because no C++ compiler has been found on this system. Some parts of the Toolkit will not function. Install a C++ compiler.' retry_instructions], []);
    else
        disp('***********************************************************************');
        disp(' CoreMat - running MEX steup ');
        disp(' ');
        disp(' Follow the instructions below to select a C++ compiler. Normally the default options are good');
        disp(' ');
        disp('***********************************************************************');
        mex -setup;
        compiler = GetNameOfCppCompiler;
        if isempty(compiler)
            reporting.ShowWarning('CoreCompileMexFiles:NoCompiler', ['I cannot compile mex files because no C++ compiler has been selected. Run mex -setup to choose your C++ compiler.' retry_instructions], []);
        end
    end
end

function Compile(mex_files_to_compile, framework_cache, cached_mex_file_info, output_directory, compiler, cuda_compiler, force_recompile, retry_instructions, reporting)
    progress_message_showing_compile = false;
    for mex_file_s = mex_files_to_compile.values
        mex_file = mex_file_s{1};

        if (~mex_file.NeedsRecompile) && (~force_recompile)
            reporting.LogVerbose([mex_file.Name ' is up to date']);
            if ~strcmp(mex_file.StatusID, 'CoreCompileMexFiles:NoRecompileNeeded')
                reporting.Error('CoreCompileMexFiles:WrongStatus', 'Program error: mex status should be OK if a recompile is not required');
            end
        else
            if strcmp(mex_file.StatusID, 'CoreCompileMexFiles:VersionChanged')
                reporting.ShowMessage('CoreCompileMexFiles:VersionChanged', [mex_file.Name ' is out of date.']);
            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:TimestampChanged')
                reporting.ShowMessage('CoreCompileMexFiles:TimestampChanged', [mex_file.Name ' has been modified.']);
            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:CompiledFileRemoved')
                reporting.ShowMessage('CoreCompileMexFiles:CompiledFileRemoved', [mex_file.Name ': The compiled mex file was removed and must be recompiled.']);
            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:FileAdded')
                reporting.ShowMessage('CoreCompileMexFiles:FileAdded', ['A new mex file ' mex_file.Name ' has been found and requires compilation.']);
            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:NoCachedInfoForMex')
                reporting.Log([mex_file.Name ' requires compilation.']);
            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:CacheFileDeleted')
                reporting.ShowMessage('CoreCompileMexFiles:CacheFileDeleted', [mex_file.Name ' requires recompilation because it appears the cache file ' framework_cache.GetCacheFilename ' was deleted.']);
            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:CompiledFileNotFound')
                reporting.Log([mex_file.Name ' requires compilation.']);

            elseif strcmp(mex_file.StatusID, 'CoreCompileMexFiles:NoRecompileNeeded')
                if force_recompile
                    reporting.ShowMessage('CoreCompileMexFiles:ForcingRecompile', ['Forcing recompile of ' mex_file.Name '.']);
                else
                    reporting.Error('CoreCompileMexFiles:UnexpectedStatus', 'Program error: An unexpected status condiiton was found');
                end
            else
                reporting.Error('CoreCompileMexFiles:UnexpectedStatus', 'Program error: An unexpected status condiiton was found');
            end

            src_filename = [mex_file.Name '.' mex_file.Extension];
            src_fullfile = fullfile(mex_file.Path, src_filename);
            files = dir(src_fullfile);
            last_modified_datenum = files.date;

            if cached_mex_file_info.isKey(mex_file.Name)
                cached_mex_file = cached_mex_file_info(mex_file.Name);
                version_unchanged_since_last_compilation_attempt = isequal(cached_mex_file.LastAttemptedCompiledVersion, mex_file.CurrentVersion);
                compiler_unchanged_since_last_compilation_attempt = isequal(cached_mex_file.LastAttemptedCompiler, compiler);
                if isempty(cached_mex_file.LastAttemptedCompileDatenum)
                    datenum_unchanged_since_last_compilation_attempt = true;
                else
                    datenum_unchanged_since_last_compilation_attempt = isequal(cached_mex_file.LastAttemptedCompileDatenum, last_modified_datenum);
                end
                compiled_failed_last_time = cached_mex_file.LastCompileFailed;
                if isempty(compiled_failed_last_time)
                    compiled_failed_last_time = false;
                end

                if compiled_failed_last_time && version_unchanged_since_last_compilation_attempt && compiler_unchanged_since_last_compilation_attempt && datenum_unchanged_since_last_compilation_attempt;
                    try_compilation_again = false;
                else
                    try_compilation_again = true;
                end
            else
                try_compilation_again = true;
            end

            if ~(try_compilation_again || force_recompile)
                reporting.ShowWarning('CoreCompileMexFiles:NotRecompiling', ['The mex source file ' src_fullfile ' needs recompilation, but the previous attempt to compile failed so I am not going to try again. You need to force a recompilation.' retry_instructions], []);
            else
                if ~exist(src_fullfile, 'file')
                    reporting.ShowWarning('CoreCompileMexFiles:SourceFileNotFound', ['The mex source file ' src_fullfile ' was not found. If this file has been removed it should also be removed from the list of mex files'], []);
                else
                    if ~progress_message_showing_compile
                        progress_message_showing_compile = true;
                        reporting.ShowProgress('Compiling mex files');
                    end
                    use_cuda = false;
                    if isa(mex_file, 'CoreCudaInfo')
                        if isempty(cuda_compiler)
                            reporting.ShowWarning('CoreCompileMexFiles:NoCudaCompilerFound', ['The mex source file ' src_fullfile ' requires the CUDA SDK to be installed.'], []);
                        else
                            use_cuda = true;
                        end
                    end

                    if use_cuda
                        mex_result = CoreCudaCompile.Compile(cuda_compiler, mex_file, src_fullfile, output_directory);
                    else
                        mex_result = CoreMexCompile.Compile(compiler, mex_file, src_fullfile, output_directory);
                    end
                    mex_file.LastAttemptedCompiledVersion = mex_file.CurrentVersion;
                    mex_file.LastAttemptedCompiler = compiler;
                    mex_file.LastAttemptedCompileDatenum = last_modified_datenum;
                    if (mex_result == 0)
                        mex_file.LastCompileFailed = false;
                        mex_file.LastSuccessfulCompiledVersion = mex_file.CurrentVersion;
                        mex_file.LastSuccessfulCompileDatenum = last_modified_datenum;
                        mex_file.LastSuccessfulCompiler = compiler;
                        reporting.ShowMessage('CoreCompileMexFiles:MexCompilationSucceeded', [' - ' src_filename ' compiled successfully.']);
                    else
                        mex_file.LastCompileFailed = true;
                        reporting.ShowWarning('CoreCompileMexFiles:MexCompilationFailed', [' - ' src_fullfile ' failed to compile.'], []);
                    end
                end
            end
        end
    end
    if progress_message_showing_compile
        reporting.CompleteProgress;
    end
end

function CheckMexFiles(mex_files_to_compile, cached_mex_file_info, output_directory, framework_cache_was_missing, reporting)
    mex_extension = mexext;
    if ~isempty(mex_extension)
        mex_extension = ['.' mex_extension];
    end

    % Iterate through all the mex files to compile
    for mex_file_s = mex_files_to_compile.values
        mex_file = mex_file_s{1};

        mex_file.NeedsRecompile = false;

        % Check whether the currently compiled file is up to date
        if cached_mex_file_info.isKey(mex_file.Name)
            cached_mex_file = cached_mex_file_info(mex_file.Name);

            file_info = dir(fullfile(mex_file.Path, [mex_file.Name '.' mex_file.Extension]));
            if numel(file_info) > 0
                current_file_timestamp = file_info.date;
            else
                current_file_timestamp = [];
            end

            % Copy across the cached values
            mex_file.LastSuccessfulCompiledVersion = cached_mex_file.LastSuccessfulCompiledVersion;
            mex_file.LastSuccessfulCompiler = cached_mex_file.LastSuccessfulCompiler;
            mex_file.LastAttemptedCompiledVersion = cached_mex_file.LastAttemptedCompiledVersion;
            mex_file.LastAttemptedCompiler = cached_mex_file.LastAttemptedCompiler;
            mex_file.LastCompileFailed = cached_mex_file.LastCompileFailed;
            mex_file.LastAttemptedCompileDatenum = cached_mex_file.LastAttemptedCompileDatenum;
            mex_file.LastSuccessfulCompileDatenum = cached_mex_file.LastSuccessfulCompileDatenum;

            if ~isequal(cached_mex_file.LastSuccessfulCompiledVersion, mex_file.CurrentVersion)
                mex_file.StatusID = 'CoreCompileMexFiles:VersionChanged';
                mex_file.NeedsRecompile = true;

            elseif ~isequal(cached_mex_file.LastSuccessfulCompileDatenum, current_file_timestamp)
                mex_file.StatusID = 'CoreCompileMexFiles:TimestampChanged';
                mex_file.NeedsRecompile = true;

            else
                mex_file.StatusID = 'CoreCompileMexFiles:NoRecompileNeeded';
            end

            % Check whether compiled file exists. We only need to do this if the
            % cache was found, because otherwise we will be recompiling anyway
            if isa(mex_file, 'CoreCudaInfo')
                compiled_file_extension = '.ptx';
            else
                compiled_file_extension = mex_extension;
            end
            dst_fullfile = fullfile(output_directory, [mex_file.Name, compiled_file_extension]);
            if ~exist(dst_fullfile, 'file')
                mex_file.NeedsRecompile = true;

                % Only display a warning message here if there is a record of a
                % previous compilation
                if (~isempty(mex_file.LastSuccessfulCompiledVersion)) && (~mex_file.LastCompileFailed)
                    mex_file.StatusID = 'CoreCompileMexFiles:CompiledFileRemoved';
                else
                    mex_file.StatusID = 'CoreCompileMexFiles:CompiledFileNotFound';
                end
            end

        else
            % No cached file was found. No need to check for existance of
            % output file as we will try to recompile anyway
            if framework_cache_was_missing
                dst_fullfile = fullfile(output_directory, [mex_file.Name, mex_extension]);
                if exist(dst_fullfile, 'file')
                    mex_file.StatusID = 'CoreCompileMexFiles:CacheFileDeleted';
                else
                    mex_file.StatusID = 'CoreCompileMexFiles:NoCachedInfoForMex';
                end
            else
                mex_file.StatusID = 'CoreCompileMexFiles:FileAdded';
            end
            mex_file.NeedsRecompile = true;
        end

    end
end

function name = GetNameOfCppCompiler()
    cc = mex.getCompilerConfigurations('C++', 'Selected');
    if isempty(cc)
        name = [];
    else
        name = cc(1).Name;
    end
end

function cuda_compiler = GetCudaCompiler()
    cuda_compiler = CoreSystemUtilities.FindFirstExecutableOnPath('nvcc');
    if isempty(cuda_compiler)
        if ispc
            cuda_compiler = TryToFindCudaCompilerPc(fullfile(getenv('ProgramFiles'), 'NVIDIA GPU Computing Toolkit', 'CUDA'));
            if isempty(cuda_compiler)
                cuda_compiler = TryToFindCudaCompilerPc(fullfile(getenv('ProgramW6432'), 'NVIDIA GPU Computing Toolkit', 'CUDA'));
            end
            if isempty(cuda_compiler)
                cuda_compiler = TryToFindCudaCompilerPc(fullfile(getenv('ProgramFiles(x86)'), 'NVIDIA GPU Computing Toolkit', 'CUDA'));
            end
        else
            if 2 == exist('/usr/local/cuda/bin/nvcc', 'file')
                cuda_compiler = '/usr/local/cuda/bin/nvcc';
            else
                cuda_compiler = [];
            end
        end
    end
end

function compiler = TryToFindCudaCompilerPc(base_dir)
    compiler = [];
    directories = CoreDiskUtilities.GetListOfDirectories(base_dir);

    for dir_name = directories
        bin_dir = fullfile(base_dir, dir_name{1}, 'bin');
        if isdir(bin_dir)
            compiler = fullfile(bin_dir, 'NVCC.EXE');
            return;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [r_I_IS1_m,r_I_IS2_m,r_I_IS3_m,A_S1_I,A_S2_I,A_S3_I] = Transformation_Sensor(r_I_IK_m,A_K_I,r_K_KS1_m,r_K_KS2_m,r_K_KS3_m,phi0_K_KS1_rad,phi0_K_KS2_rad,phi0_K_KS3_rad)
%% Computation Sensor Position and Orientation
% Input parameters:
% phi   [rad]    Vehicle Rotation Angles [Phi Theta Psi]
% chi   [rad]    Road Cardan Angles [chi_x chi_y Psi]
% delta [rad]    Steering Angles [deltaFL deltaFR deltaRL deltaRR]
% X     [---]    Original Input Vectors in Vehicle Fixed Coordinate System [x1 x2x3 x4] [3x4]
% Output parameters:
% Y     [---]    Transformed Output Vectors [y1 y2 y3 y4] [3x4]

%% Definition of rotation matrices

% General Rotation Matrix between Coordinate Systems
A =@(phi) reshape([cos(phi(3))*cos(phi(2)),-cos(phi(1))*sin(phi(3))+cos(phi(3))*sin(phi(2))*sin(phi(1)),sin(phi(3))*sin(phi(1))+cos(phi(3))*cos(phi(1))*sin(phi(2)),cos(phi(2))*sin(phi(3)),cos(phi(3))*cos(phi(1))+sin(phi(3))*sin(phi(2))*sin(phi(1)),-cos(phi(3))*sin(phi(1))+cos(phi(1))*sin(phi(3))*sin(phi(2)),-sin(phi(2)),cos(phi(2))*sin(phi(1)),cos(phi(2))*cos(phi(1))],[3,3]);

% Rotation from Vehicle Fixed Reference Frame into Sensor Fixed Reference Frame of Sensor 1 (Kardan-Rotation)
A_S1_K = A(phi0_K_KS1_rad);
% Rotation from Vehicle Fixed Reference Frame into Sensor Fixed Reference Frame of Sensor 2 (Kardan-Rotation)
A_S2_K = A(phi0_K_KS2_rad);
% Rotation from Vehicle Fixed Reference Frame into Sensor Fixed Reference Frame of Sensor 3 (Kardan-Rotation)
A_S3_K = A(phi0_K_KS3_rad);

%% Output of transformed sensors

% Absolute Position of Sensors in Inertial Reference Frame
r_I_IS1_m = r_I_IK_m + transpose(A_K_I)*r_K_KS1_m;
r_I_IS2_m = r_I_IK_m + transpose(A_K_I)*r_K_KS2_m;
r_I_IS3_m = r_I_IK_m + transpose(A_K_I)*r_K_KS3_m;

% Absolute Orientation of Sensors in Inertial Reference Frame
A_S1_I = A_S1_K*A_K_I;
A_S2_I = A_S2_K*A_K_I;
A_S3_I = A_S3_K*A_K_I;

end

author: TUMFTM


text:
function SPM = spm_contrasts(SPM,Ic)
% Compute and store contrast parameters and inference SPM{.}
% FORMAT SPM = spm_contrasts(SPM,Ic)
%
% SPM  - SPM data structure
% Ic   - indices of xCon to compute
%
% This function fills in SPM.xCon and writes con_????, ess_???? and
% spm?_???? images.
%__________________________________________________________________________
% Copyright (C) 2002-2017 Wellcome Trust Centre for Neuroimaging

% Karl Friston, Will Penny & Guillaume Flandin
% $Id: spm_contrasts.m 7738 2019-12-02 12:45:37Z guillaume $

% Temporary copy of the SPM variable, to avoid saving it in SPM.mat unless
% it has changed (faster, read-only access)
%--------------------------------------------------------------------------
tmpSPM = SPM;

%-Change to results directory
%--------------------------------------------------------------------------
try, cd(SPM.swd); end

%-Get contrast definitions (if available)
%--------------------------------------------------------------------------
try
    xCon = SPM.xCon;
catch
    xCon = [];
end

%-Set all contrasts by default
%--------------------------------------------------------------------------
if nargin < 2
    Ic   = 1:length(xCon);
end
Ic(Ic == 0) = [];

%-Map parameter and hyperarameter files
%--------------------------------------------------------------------------
if ~isempty(xCon) && xCon(Ic(1)).STAT == 'P'

    %-Conditional estimators
    %----------------------------------------------------------------------
    Vbeta = SPM.VCbeta;
else

    %-OLS estimators and error variance estimate
    %----------------------------------------------------------------------
    Vbeta = SPM.Vbeta;
    VHp   = SPM.VResMS;
end

if spm_mesh_detect(Vbeta)
    file_ext = '.gii';
    g        = SPM.xY.VY(1).private;
    metadata = g.private.metadata;
    name     = {metadata.name};
    if any(ismember(name,'SurfaceID'))
        metadata = metadata(ismember(name,'SurfaceID'));
        metadata = {metadata.name, metadata.value};
    elseif isfield(g,'faces') && ~isempty(g.faces)
        metadata = {'SurfaceID', SPM.xY.VY(1).fname};
    else
        metadata = {};
    end
else
    file_ext = spm_file_ext;
    metadata = {};
end

%-Compute & store contrast parameters, contrast/ESS images, & SPM images
%==========================================================================
spm('Pointer','Watch')
XYZ  = SPM.xVol.XYZ;
iXYZ = cumprod([1,SPM.xVol.DIM(1:2)'])*XYZ - sum(cumprod(SPM.xVol.DIM(1:2)'));
for i = 1:length(Ic)

    %-Canonicalise contrast structure with required fields
    %----------------------------------------------------------------------
    ic  = Ic(i);
    if isempty(xCon(ic).eidf)
        X1o           = spm_FcUtil('X1o',xCon(ic),SPM.xX.xKXs);
        [trMV,trMVMV] = spm_SpUtil('trMV',X1o,SPM.xX.V);
        xCon(ic).eidf = trMV^2/trMVMV;
    end

    %-Write contrast/ESS images?
    %======================================================================
    if isempty(xCon(ic).Vcon)

        switch xCon(ic).STAT

            case {'T','P'}

                if strcmp(xCon(ic).STAT,'P') && strcmp(SPM.PPM.xCon(ic).PSTAT,'F')

                    % Bayes Factor for compound contrast
                    %------------------------------------------------------
                    disp('Bayes factor for compound contrast');
                    fprintf('\t%-32s: %30s',sprintf('LogBF image %2d',ic),...
                        '...computing');                                %-#

                    if isfield(SPM.PPM,'VB')
                        % First level Bayes
                        xCon = spm_vb_logbf(SPM,XYZ,xCon,ic);
                    else
                        % Second level Bayes
                        xCon = spm_bayes2_logbf(SPM,XYZ,xCon,ic);
                    end
                else
                    %-Implement contrast as linear combination of beta images
                    %------------------------------------------------------
                    fprintf('\t%-32s: %30s',sprintf('contrast image %2d',ic),...
                        '...computing');                                %-#

                    %-Prepare handle for contrast image
                    %------------------------------------------------------
                    xCon(ic).Vcon = struct(...
                        'fname',  [sprintf('con_%04d',ic) file_ext],...
                        'dim',    SPM.xVol.DIM',...
                        'dt',     [spm_type('float32'), spm_platform('bigend')],...
                        'mat',    SPM.xVol.M,...
                        'pinfo',  [1,0,0]',...
                        'descrip',sprintf('Contrast %d: %s',ic,xCon(ic).name),...
                        metadata{:});

                    xCon(ic).Vcon = spm_data_hdr_write(xCon(ic).Vcon);

                    %-Compute contrast
                    %------------------------------------------------------
                    Q      = find(abs(xCon(ic).c) > 0);
                    V      = Vbeta(Q);

                    cB     = zeros(1,size(XYZ,2));
                    for j=1:numel(V)
                        cB = cB + xCon(ic).c(Q(j)) * spm_data_read(V(j),'xyz',XYZ);
                    end

                    %-Write contrast image
                    %------------------------------------------------------
                    tmp = NaN(SPM.xVol.DIM');
                    tmp(iXYZ) = cB;
                    xCon(ic).Vcon = spm_data_write(xCon(ic).Vcon,tmp);

                    clear tmp cB
                    fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(...
                        '...written %s',spm_file(xCon(ic).Vcon.fname,'filename')))%-#

                end

            case 'F' %-Implement ESS as sum of squared weighted beta images
                %----------------------------------------------------------
                fprintf('\t%-32s: %30s',sprintf('ESS image %2d',ic),...
                    '...computing');                                    %-#

                %-Prepare handle for ESS image
                %----------------------------------------------------------
                xCon(ic).Vcon = struct(...
                    'fname',  [sprintf('ess_%04d',ic) file_ext],...
                    'dim',    SPM.xVol.DIM',...
                    'dt',     [spm_type('float32'), spm_platform('bigend')],...
                    'mat',    SPM.xVol.M,...
                    'pinfo',  [1,0,0]',...
                    'descrip',sprintf('ESS contrast %d: %s',ic,xCon(ic).name),...
                    metadata{:});

                xCon(ic).Vcon = spm_data_hdr_write(xCon(ic).Vcon);

                %-Compute ESS
                %----------------------------------------------------------
                % Residual (in parameter space) forming matrix
                h  = spm_FcUtil('Hsqr',xCon(ic),SPM.xX.xKXs);
                ss = zeros(numel(Vbeta),size(XYZ,2));
                for j=1:numel(Vbeta)
                    ss(j,:) = spm_data_read(Vbeta(j),'xyz',XYZ);
                end
                ss = sum((h*ss).^2,1);

                %-Write ESS image
                %----------------------------------------------------------
                tmp = NaN(SPM.xVol.DIM');
                tmp(iXYZ) = ss;
                xCon(ic).Vcon = spm_data_write(xCon(ic).Vcon,tmp);

                clear tmp ss
                fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(...
                        '...written %s',spm_file(xCon(ic).Vcon.fname,'filename')))%-#

            otherwise
                %----------------------------------------------------------
                error(['unknown STAT "',xCon(ic).STAT,'"'])

        end % (switch(xCon...)

    end % (if isempty(xCon(ic)...)

    %-Write inference SPM/PPM
    %======================================================================
    if isempty(xCon(ic).Vspm) || xCon(ic).STAT == 'P'
        % (always update PPM as size threshold, gamma, may have changed)

        fprintf('\t%-32s: %30s',sprintf('spm{%s} image %2d',xCon(ic).STAT,ic),...
            '...computing');                                            %-#

        switch(xCon(ic).STAT)

            case 'T'                                 %-Compute SPM{t} image
                %----------------------------------------------------------
                cB  = spm_data_read(xCon(ic).Vcon,'xyz',XYZ);
                l   = spm_data_read(VHp,'xyz',XYZ);    % get hyperparamters
                Vc  = xCon(ic).c'*SPM.xX.Bcov*xCon(ic).c;
                SE  = sqrt(l*Vc);                      % and standard error
                Z   = cB./SE;
                str = sprintf('[%.1f]',SPM.xX.erdf);

            case 'P'                                 %-Compute PPM{P} image
                %----------------------------------------------------------

                if all(strcmp({SPM.PPM.xCon(ic).PSTAT},'T'))
                    % Simple contrast - Gaussian distributed

                    c     = xCon(ic).c;
                    cB    = spm_data_read(xCon(ic).Vcon,'xyz',XYZ);
                    if isfield(SPM.PPM,'VB');
                        % If posterior sd image for that contrast does
                        % not already exist, then compute it
                        try
                            SPM.PPM.Vcon_sd(ic);
                        catch
                            SPM = spm_vb_contrasts(SPM,XYZ,xCon,ic);
                        end
                        % Read in posterior sd image for contrast
                        Vsd = spm_data_read(SPM.PPM.Vcon_sd(ic),'xyz',XYZ);
                        VcB = Vsd.^2;
                    else
                        VcB   = c'*SPM.PPM.Cby*c;
                        for j = 1:length(SPM.PPM.l)

                            % hyperparameter and Taylor approximation
                            %----------------------------------------------
                            l   = spm_data_read(SPM.VHp(j),'xyz',XYZ);
                            VcB = VcB + (c'*SPM.PPM.dC{j}*c)*(l - SPM.PPM.l(j));
                        end
                    end

                    % posterior probability cB > g
                    %------------------------------------------------------
                    Gamma          = xCon(ic).eidf;
                    Z              = 1 - spm_Ncdf(Gamma,cB,VcB);

                    % Convert probability to Log Odds Ratio
                    Z              = log( Z ./ (1 - Z+eps) );
                    str            = sprintf('[%.2f]',Gamma);
                    %xCon(ic).name = [xCon(ic).name ' ' str];
                else
                    % Compound contrast - Log Bayes Factor
                    fprintf('\t\t%-75s\n','Log Bayes Factor for compound contrast');
                    fprintf('\t%-32s: %29s\n',' ',' ');
                    Z = spm_data_read(xCon(ic).Vcon,'xyz',XYZ);

                    str = sprintf('[%1.2f]',xCon(ic).eidf);
                end

            case 'F'                                 %-Compute SPM{F} image
                %----------------------------------------------------------
                MVM = spm_data_read(xCon(ic).Vcon,'xyz',XYZ)/trMV;
                RVR = spm_data_read(VHp,'xyz',XYZ);
                Z   = MVM./RVR;
                str = sprintf('[%.1f,%.1f]',xCon(ic).eidf,SPM.xX.erdf);

            otherwise
                %----------------------------------------------------------
                error(['unknown STAT "',xCon(ic).STAT,'"']);

        end % (switch(xCon(ic)...)

        %-Write SPM - statistic image
        %------------------------------------------------------------------
        xCon(ic).Vspm = struct(...
            'fname',  [sprintf('spm%s_%04d',xCon(ic).STAT,ic) file_ext],...
            'dim',    SPM.xVol.DIM',...
            'dt',     [spm_type('float32'), spm_platform('bigend')],...
            'mat',    SPM.xVol.M,...
            'pinfo',  [1,0,0]',...
            'descrip',sprintf('SPM{%s_%s} - contrast %d: %s',...
                xCon(ic).STAT,str,ic,xCon(ic).name),...
            metadata{:});

        xCon(ic).Vspm = spm_data_hdr_write(xCon(ic).Vspm);

        tmp           = zeros(SPM.xVol.DIM');
        tmp(iXYZ)     = Z;
        xCon(ic).Vspm = spm_data_write(xCon(ic).Vspm,tmp);

        clear tmp Z
        cmd = sprintf(['[hReg,xSPM,SPM] = spm_results_ui(''Setup'',',...
            'struct(''swd'',''%s'',''Ic'',%d));',...
            'TabDat = spm_list(''List'',xSPM,hReg);'],pwd,ic);
        img = spm_file(spm_file(xCon(ic).Vspm.fname,'filename'),'link',cmd);
        n = 30; if length(img)>n, n = length(img)+n-13; end
        fprintf('%s%*s\n',repmat(sprintf('\b'),1,30),n,sprintf(...
            '...written %s',img)); %-#

    end % (if isempty(xCon(ic)...)

end % (for i = 1:length(Ic))
spm('Pointer','Arrow')

% place xCon back in SPM
%--------------------------------------------------------------------------
SPM.xCon = xCon;

% Check if SPM has changed. Save only if it has.
%--------------------------------------------------------------------------
if spm_check_version('matlab','8.0') >= 0, my_isequaln = @isequaln;
else my_isequaln = @isequalwithequalnans; end
if ~my_isequaln(tmpSPM,SPM)
    fprintf('\t%-32s: %30s','Saving SPM.mat','...writing');             %-#
    fmt = spm_get_defaults('mat.format');
    s = whos('SPM');
    if s.bytes > 2147483647, fmt = '-v7.3'; end
    save('SPM.mat', 'SPM', fmt);
    fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...SPM.mat saved')   %-#
end

author: spm


text:
function C = arr2cell(X, dim)
%ARR2CELL Splits an array into a cell across the specified dimension.
% Usage:
%   C = arr2cell(X) % default: last dimension
%   C = arr2cell(X, dim)
%
% See also: num2cell, stack2cell

if nargin < 2; dim = ndims(X); end

% Select every other dimension
dims = setdiff(1:ndims(X), dim);

% Convert to cell and squeeze
C = squeeze(num2cell(X,dims));

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function gui_check_vector(edithandle,resetstr)
% check whether the the given edit field contains a row vector; warn and reset if not
if ~exist('resetstr','var')
    resetstr = ''; end

str = get(edithandle,'String');
if ~isempty(str) && ~isvector(str2num(str))
    errordlg2('This field must be of the form [number, number, number, ...]');
    set(edithandle,'String',resetstr);
end

repo: REST


text:
function expInfo = exportInfo(guifile)
%EXPORTINFO    Extract useful info from a CHEBGUI object for exporting.
%
% Calling sequence
%
%   EXPINFO = EXPORTINFO(GUIFILE)
%
% where
%
%   GUIFILE:    A CHEBGUI object
%   EXPINFO:    A struct, containing fields with information for exporting to an
%               .m-file.

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

% Extract information from the GUI fields
dom = guifile.domain;
domNum = str2num(dom);
a = num2str(domNum(1));
b = num2str(domNum(end));
deInput = guifile.DE;
lbcInput = guifile.LBC;
rbcInput = guifile.RBC;
deRHSInput = 'u_t';
initInput = guifile.init;
tt = guifile.timedomain;

% Wrap all input strings in a cell (if they're not a cell already)
if ( isa(deInput, 'char') )
    deInput = cellstr(deInput);
end

if ( isa(lbcInput, 'char') )
    lbcInput = cellstr(lbcInput);
end

if ( isa(rbcInput, 'char') )
    rbcInput = cellstr(rbcInput);
end

if ( isa(deRHSInput, 'char') )
    deRHSInput = cellstr(deRHSInput);
end

if ( isa(initInput, 'char') )
    initInput = cellstr(initInput);
end

lbcRHSInput = cellstr(repmat('0', numel(lbcInput), 1));
rbcRHSInput = cellstr(repmat('0', numel(rbcInput), 1));
initRHSInput = cellstr(repmat('0', numel(initInput), 1));

% Obtain strings for setting up the problem:
[deString, allVarString, indVarNameDE, pdeVarName, pdeflag, ignored, allVarNames] ...
    = setupFields(guifile, deInput, 'DE');
if ( ~any(pdeflag) )
    error('CHEBFUN:CHEBGUIEXPORTERPDE:exportInfo:notPDE', ...
        ['Input does not appear to be a PDE, ', ...
         'or at least is not a supported type.']);
end

% Obtain the independent variable name appearing in the initial condition
if ( ~isempty(initInput{1}) )
    [initString, ignored, indVarNameInit] = ...
        setupFields(guifile, initInput, 'BC', allVarString);
else
    indVarNameInit = {''};
end

% Make sure we don't have a disrepency in indVarNames. Create a new
% variable indVarName which contains both independent variable, the first
% entry corresponds to space, the second to time.
if ( ~isempty(indVarNameInit{1}) && ~isempty(indVarNameDE{1}) )
    if ( strcmp(indVarNameDE{1}, indVarNameInit{1}) )
        indVarName{1} = indVarNameDE{1};
    else
        error('CHEBFUN:CHEBGUIEXPORTERPDE:exportInfo:solveGUIpde', ...
            'Independent variable names do not agree')
    end
elseif ( ~isempty(indVarNameInit{1}) && isempty(indVarNameDE{1}) )
    indVarName{1} = indVarNameInit{1};
elseif ( isempty(indVarNameInit{1}) && ~isempty(indVarNameDE{1}) )
    indVarName{1} = indVarNameDE{1};
else
    indVarName{1} = 'x'; % Default value
end

if ( ~isempty(indVarNameDE{2}) )
    indVarName{2} = indVarNameDE{2}; % This should never be empty for a PDE.
else
    indVarName{2} = 't'; % Default value
end

% Check that we're not using the same variable for space and time!
if ( strcmp(indVarName{1}, indVarName{2}) )
     error('CHEBFUN:CHEBGUIEXPORTERPDE:exportInfo:solveGUIpde', ...
        'The same variable appears to be used as space and time variable');
end

% Ensure that we have a PDE subscript:
if ( ~isempty(indVarNameDE{2}) )
    % Find what the name of the time variable specified is:
    indVarName{2} = indVarNameDE{2};
else
    error('CHEBFUN:CHEBGUIEXPORTERPDE:exportInfo:solveGUIpde', ...
        'No time variable specified, please do so via using subscript.');
end

% Create a string with the variables used in the problem
variableString = [indVarName{2}, ',', indVarName{1}, ','];

% The names of the variables used for space and time:
xName = indVarName{1};
tName = indVarName{2};

% Support for periodic boundary conditions
if ( (~isempty(lbcInput{1}) && strcmpi(lbcInput{1}, 'periodic')) ...
     || (~isempty(rbcInput{1}) && strcmpi(rbcInput{1}, 'periodic')) )
    lbcInput{1} = [];
    rbcInput{1} = [];
    periodic = true;
else
    periodic = false;
end

% Do we have a left BC specified?
if ( ~isempty(lbcInput{1}) )
    lbcString = setupFields(guifile, lbcInput, 'BC', allVarString);

    idx = strfind(lbcString, ')');
    if ( ~isempty(idx) )
        % Inject the t and x variables into the anonymous function:
        lbcString = [lbcString(1:2), tName, ',', lbcString(3:end)];

    end

else
    lbcString = [];
end

% Do we have a right BC specified?
if ( ~isempty(rbcInput{1}) )
    rbcString = setupFields(guifile, rbcInput, 'BC', allVarString);

    idx = strfind(rbcString, ')');
    if ( ~isempty(idx) )

        % Inject the t and x variables into the anonymous function:
        rbcString = [rbcString(1:2), tName,',', rbcString(3:end)];

    end

else
    rbcString = [];
end

% Glue the DESTRING together
deString = ['@(' variableString, deString(3:end)];

% Find what the dependent variables are:
s = allVarNames;
if ( (numel(deInput) == 1) && ~ischar(deInput) )

    % Get the strings of the dependent variable. Just use allVarNames.
    sol = s{1};
    sol0 = [sol '0'];

else

    % These can be changed
    sol0 = 'sol0';
    sol = 'sol';

end

%% Fill up the expInfo struct
expInfo.dom = dom;
expInfo.deInput = deInput;

% PDE specific information
expInfo.a = a;
expInfo.b = b;
expInfo.xName = xName;
expInfo.tName = tName;
expInfo.tt = tt;
expInfo.deRHSInput = deRHSInput;
expInfo.lbcRHSInput = lbcRHSInput;
expInfo.rbcRHSInput = rbcRHSInput;
expInfo.initRHSInput = initRHSInput;
expInfo.pdeflag = pdeflag;
expInfo.initInput = initInput;
expInfo.s = s;
expInfo.sol = sol;
expInfo.sol0 = sol0;

% And then some...
expInfo.deString = deString;
expInfo.pdeVarName = pdeVarName;
expInfo.initString = initString;
expInfo.lbcInput = lbcInput;
expInfo.rbcInput = rbcInput;
expInfo.allVarString = allVarString;
expInfo.allVarNames = allVarNames;
expInfo.indVarName = indVarName;
expInfo.periodic = periodic;
expInfo.lbcString = lbcString;
expInfo.rbcString = rbcString;
% Input related to options
expInfo.tol = guifile.tol;
expInfo.doplot = guifile.options.plotting;
expInfo.dohold = guifile.options.pdeholdplot;
expInfo.ylim1 = guifile.options.fixYaxisLower;
expInfo.ylim2 = guifile.options.fixYaxisUpper;
expInfo.fixN = guifile.options.fixN;
expInfo.pdeSolver = guifile.options.pdeSolver;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [f_hat]=gsp_gft(G,f)
%GSP_GFT  Graph Fourier transform
%   Usage:  f_hat=gsp_gft(G,f);
%
%   Input parameters:
%         G          : Graph or Fourier basis
%         f          : f (signal)
%   Output parameters:
%         f_hat      : Graph Fourier transform of *f*
%
%   'gsp_gft(G,f)' computes a graph Fourier transform of the signal $f$
%   with respect to the Fourier basis of the graph G: G.U. Alternatively,
%   one can provide directly the Fourier basis instead of the graph G.
%
%   .. f_hat = U' * f
%
%   .. math:: \hat{f}(\lambda_{\ell})=\langle f , u_{\ell} \rangle
%
%   To compute the Fourier basis of a graph G, you can use the function::
%
%           G = gsp_compute_fourier_basis(G);
%
%   Example:::
%
%       N = 30;
%       G = gsp_sensor(N);
%       G = gsp_compute_fourier_basis(G);
%       f = sin((1:N)'*2*pi/N);
%       fhat = gsp_gft(G,f);
%       gsp_plot_signal_spectral(G,fhat);
%
%   See also: gsp_igft, gsp_gwft, gsp_compute_fourier_basis
%

% Author : Nathanael Perraudin, David I Shuma
% Date:

s = size(f);

if ~isnumeric(G);
    if ~gsp_check_fourier(G)
       error(['GSP_GFT: You need first to compute the Fourier basis\n',...
           'You can do it with the function gsp_compute_fourier_basis']);
    end
    %U = G.U;
    f_hat=reshape(G.U'*reshape(f,G.N,[]),s);
else
    %U = G;
    f_hat=reshape(G'*reshape(f,G.N,[]),s);
end

repo: gspbox


text:
% DEMVOWELSFGPLVM1 Model the vowels data with a 2-D FGPLVM using RBF kernel and back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'vowels';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.back = 'mlp';
options.backOptions = mlpOptions;
options.numActive = 200;
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(dataSetName, experimentNo, 'vector')

errors = lvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function k = my_discrete_gaussian_curvature(V,F)
% MY_DISCRETE_GAUSSIAN_CURVATURE Compute discrete gaussian curvature
% using the angle defect
%
% k = my_discrete_gaussian_curvature(V,F)
%
% Inputs:
%   V  #V by 3 list of vertex positions
%   F  #F by 3 list of face indies
% Outputs:
%   k  #V by 1 list of discrete gaussian curvature values
%

k = 2*pi*ones(size(V,1),1);

%   for i=1:size(F,1)
%       for j=1:3
%           e1 = V(F(i,mod(j,3)+1),:) - V(F(i,j),:);
%           e2 = V(F(i,mod(j+1,3)+1),:) - V(F(i,j),:);
%           e1 = e1 / norm(e1);
%           e2 = e2 / norm(e2);
%           k(F(i,j)) = k(F(i,j)) - acos(e1*e2');
%       end
%   end

for i=1:size(F,1)
    e1 = V(F(i,2),:) - V(F(i,1),:);
    e2 = V(F(i,3),:) - V(F(i,1),:);
    e3 = V(F(i,3),:) - V(F(i,2),:);

    e1 = e1 / norm(e1);
    e2 = e2 / norm(e2);
    e3 = e3 / norm(e3);

    k(F(i,1)) = k(F(i,1)) - acos(e1*e2');
%
%     e1 = V(F(i,3),:) - V(F(i,2),:);
%     e2 = V(F(i,1),:) - V(F(i,2),:);
%     e1 = e1 / norm(e1);
%     e2 = e2 / norm(e2);
%     k(F(i,2)) = k(F(i,2)) - acos(e1*e2');
    k(F(i,2)) = k(F(i,2)) - acos(-e3*e1');

    e1 = V(F(i,1),:) - V(F(i,3),:);
    e2 = V(F(i,2),:) - V(F(i,3),:);
    e1 = e1 / norm(e1);
    e2 = e2 / norm(e2);
    k(F(i,3)) = k(F(i,3)) - acos(e1*e2');
end

% for t=1:size(F,1) %loops through each triagle
%     e1=V(F(t,2),:)-V(F(t,1),:); %finds each edge of triangle t
%     e2=V(F(t,3),:)-V(F(t,2),:);
%     e3=V(F(t,1),:)-V(F(t,3),:);
%
%     e1 = e1/norm(e1);
%     e2 = e2/norm(e2);
%     e3 = e3/norm(e3);
%
%     %finds theta between each pair of edges and subtracts from 2pi to vetex in common
%     k(F(t,1)) = k(F(t,1)) - acos(dot(e1,e3));
%     k(F(t,2)) = k(F(t,2)) - acos(dot(e1,e2));
%     k(F(t,3)) = k(F(t,3)) - acos(dot(e3,e2));
% end

end

author: odedstein

repo: sgi-introduction-course


text:
function [ans, Prefix] = PrefixStatus (Prefix)
%
%   [Status, Prefix] = PrefixStatus (Name)
%
%Purpose:
%   Determines if the name chosen is a valid for an AFNI Brick
%   If you pass a name containing prefix only and no view then
%   View is assumed to be +orig
%
%
%Input Parameters:
%   Name: a string containing an AFNI Prefix,
%          or Prefix and view moon+orig or
%           Prefix, view and extension moon+acpc.HEAD
%
%
%
%Output Parameters:
%   Status : -1 error
%          : 0 Read Only
%          : 1 Read Write
%   Prefix: The Prefix
%
%
%Key Terms:
%
%More Info :
%
%   The function checks for the existence of a brick with a similar
%   prefix
%
%
%     Author : Ziad Saad
%     Date : Wed Sep 19 10:30:05 EDT 2001
%     LBC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'CheckPrefix';

%Debug Flag
DBG = 1;

%initailize return variables
ans = 0;

%make sure you have no HEAD or BRIK as an extension
[Prefix, Ext] = RemoveExtension(Prefix, '.HEAD|.BRIK');

%make sure you have no view
[Prefix, ViewName] = RemoveExtension(Prefix, '+orig|+acpc|+tlrc');

%check for conflict
if (~isempty(Ext) & isempty(ViewName)),
	ans = -1;
	ErrEval(FuncName,'Err_Expected a correct view (+orig, +acpc or +tlrc) but found none.');
	return;
end

if (~isempty(ViewName)), View = ViewName; else View = '+orig'; end

%look for various files
if (exist(sprintf('%s%s.HEAD', Prefix, View)) == 2), return; end
if (exist(sprintf('%s%s.BRIK', Prefix, View)) == 2), return; end

ans = 1;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [TorF, vstr, rdate] = have_feature_osqp()
%HAVE_FEATURE_OSQP  Detect availability/version info for OSQP
%
%   Feature detection function implementing 'osqp' tag for HAVE_FEATURE
%   to detect availability/version of OSQP, (Operator Splitting QP solver)
%   (https://osqp.org)
%
%   See also HAVE_FEATURE, QPS_MASTER, OSQP.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = exist('osqp', 'file') == 2;
vstr = '';
rdate = '';
if TorF
    try
        o = osqp();
        vstr = o.version();
    catch
        TorF = 0;
        fprintf('OSQP Error!\n');
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function cvx_optval = quad_over_lin( x, y, dim )

%QUAD_OVER_LIN   Internal cvx version.

error( nargchk( 2, 3, nargin ) );
sx = size( x );
sy = size( y );
nx = length( sx );
ny = length( sy );
nd = max( nx, ny );
if nargin < 3 || isempty( dim ),
    dim = [ find( sx ~= 1 ), 1 ];
    dim = dim( 1 );
elseif ~isnumeric( dim ) || dim < 0 || dim ~= floor( dim ),
    error( 'Third argument must be a dimension.' );
elseif dim == 0,
    dim = nd + 1;
    nd = dim;
end

%
% Check sizes
%

sx = [ sx, ones( 1, nd - nx ) ];
sy = [ sy, ones( 1, nd - ny ) ];
sz = sx;
sz( dim ) = 1;
if sy( dim ) ~= 1 && any( sx ~= 1 ),
    error( 'Dimensions are not compatible.' );
end
if all( sz == sy ) || all( sy == 1 ),
    need_contraction = false;
elseif all( sz == 1 ),
    sz = sy;
    need_contraction = sx( dim ) ~= 1;
    sx = sz;
    dim = 0;
else
    error( 'Dimensions are not compatible.' );
end

%
% Check curvature
%

if ~cvx_isaffine( x ),
    error( 'The first argument must be affine.' );
elseif ~isreal( y ),
    error( 'The second argument must be real.' );
elseif ~cvx_isconcave( y ),
    error( 'The second argument must be concave.' );
end

%
% Construct the problem. Note that in order to support the case
% where y is convex, we actually have to introduce a new inequality
% y2 <= y, where y2 is a new variable. That is because set membership
% constraints require affine arguments. So this is an obscure case, then,
% where automatic curvature checking in epigraph functions breaks down.
%

if any( sz == 0 ),
    cvx_optval = cvx( zeros( sz ) );
elseif any( sx == 0 ),
    cvx_begin
        y == nonnegative( sz );
    cvx_end
else
    cvx_begin
        epigraph variable z( sz )
        y = cvx_accept_concave( y );
        if need_contraction,
            x = cvx_accept_convex( norms( x, 2, dim ) );
        end
        { x, y, z } == rotated_lorentz( sx, dim, ~isreal( x ) );
    cvx_end
end

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


text:
%  Figure 7.55      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% script for Fig. 7.55
%
clf;
dp=[1 1 0];
np=[1];
nc=conv([1 1],[8.32 0.8]);
dc=conv([1 4.08],[1 0.0196]);
num=conv(np,nc);
den=conv(dp,dc);
dcl=[0 0 num]+den; % closed-loop denominator
t=0:.1:5;
sys=tf(num,dcl);
y=step(sys,t);
plot(t,y)
grid;
xlabel('Time (sec)');
ylabel('y');
title('Fig.7.55 Step response for lag compensation design');
[a,b,c,d]=tf2ss(nc,dc);
% Samplig period
ts=0.1;
% Convert controller to discrete
[ad,bd,cd,dd]=c2dm(a,b,c,d,ts,'zoh');
[ncd,dcd]=ss2tf(ad,bd,cd,dd);
% response comparison of continuous and digital control.
sim('fig7_54')
figure(1)
plot(ycd(:,1),ycd(:,2))
hold on
plot(ycd(:,1),ycd(:,3),'--')
title('Figure 7.53(a) Output Responses of Digital and Continuous Controllers')
xlabel('Time (sec)');
ylabel('y');
gtext('continuous controller')
gtext('digital controller')
grid
hold off
pause;
figure(2)
plot(ucd(:,1),ucd(:,2))
hold on
plot(ucd(:,1),ucd(:,3),'--')
title('Figure 7.55(b)  Control Responses of Digital and Continuous Controllers')
xlabel('Time (sec)');
ylabel('u');
gtext('continuous controller')
gtext('digital controller')
grid
hold off

author: Sable

repo: mcbench-benchmarks


text:
function M = prepare_mesh_fittemplate(mesh1,mesh2)

% PREPARE_MESH_FITTEMPLATE computes an affine transformation matrix between 2 point clouds
%
% This function relies on cpd toolbox from  Myronenko, see https://sites.google.com/site/myronenko/research/cpd
%
% See also FT_PREPARE_MESH

% Copyright (C) 2019, Simon Homoelle
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% add toolbox cpd
ft_hastoolbox('cpd', 2);

%
opt.corresp = 0;
opt.method  = 'affine';
opt.max_it = 100;
opt.fgt = 0;
opt.tol = 10e-12;
opt.outliers = 0.0;
opt.outliers = 0;

figure
axis equal
axis vis3d

[transform,dum] = cpd_register(mesh1,mesh2, opt);

% figure
% plot3(mesh1(:,1), mesh1(:,2), mesh1(:,3), 'r.');
% hold on
% plot3(mesh2(:,1), mesh2(:,2), mesh2(:,3), 'bo');
% plot3(transform.Y(:,1), transform.Y(:,2), transform.Y(:,3), 'b.');
% axis equal
% axis vis3d

M = eye(4,4);
M(1:3,1:3) = transform.R;
M(1:3,4)   = transform.t;

end


text:
function rgb_contrast = contrast_gray ( rgb, s )

%*****************************************************************************80
%
%% CONTRAST_RGB enhances the contrast in an RGB image.
%
%  Discussion:
%
%    The current contrast at pixel (I,J) is the difference between
%    its value P(I,J) and the average AVERAGE(I,J) of its eight
%    immediate neighbors.
%
%    The enhanced image replaces each pixel value by:
%
%      p(i,j) <= s * p(i,j) + ( 1 - s ) * average(i,j)
%
%    where C is the contrast parameter.
%
%    Since this operation assumes a 1-pixel deep layer around each
%    processed pixel, the border pixels can't be handled in quite the
%    same way.  To avoid strange artifacts, the border pixels are kept
%    at their input values, although a more careful algorithm could
%    enhance their contrast as well.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 RGB(:,:,3), the original data.
%
%    Input, real S, the sharpness.  Normally, S is greater than 1.
%    S < 1, the image is blurred.
%    S = 1, the image is unchanged.
%    1 < S, the contrast is heightened.
%
%    Output, uint8 RGB_CONTRAST(:,:,3), the data after contrast enhancement.
%
  r = contrast_gray ( rgb(:,:,1), s );
  g = contrast_gray ( rgb(:,:,2), s );
  b = contrast_gray ( rgb(:,:,3), s );

  rgb_contrast = zeros ( size ( rgb ) );

  rgb_contrast = uint8 ( rgb_contrast );

  rgb_contrast(:,:,1) = r(:,:);
  rgb_contrast(:,:,2) = g(:,:);
  rgb_contrast(:,:,3) = b(:,:);

  return
end


text:
function [node,elem,bdFlag,HB] = uniformrefine3l(node,elem,bdFlag)
%% UNIFORMREFINE3L uniformly refine a 3-D triangulation.
%
% [node,elem] = uniformrefine3l(node,elem) divides each tetrahedra into
% eight small similar sub-tetrahedrons.
%
% [node,elem,bdFlag] = uniformrefine3l(node,elem,bdFlag) also update boundary
% conditions represented by bdFlag.
%
% [node,elem,~,HB] = uniformrefine3l(node,elem) outpus HB array which is
% useful for nodal interpolation.
%
% Difference with uniformrefine3: use the shortest diagonal.
%
% Example
%
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     elem = label3(node,elem);
%     figure(1); subplot(1,3,1);
%     set(gcf,'Units','normal'); set(gcf,'Position',[0.25,0.25,0.5,0.3]);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%     [node,elem] = uniformrefine3(node,elem);
%     figure(1); subplot(1,3,2);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%     bdFlag = setboundary3(node,elem,'Dirichlet');
%     [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
%     figure(1); subplot(1,3,3);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%
% See also uniformbisect, uniformrefine, bisect, bisect3
%
% Reference: S. Zhang. Successive subdivisions of tetrahedra and multigrid
% methods on t etrahedral meshes. Houston J. Math. 21, 541?556, 1995.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Construct data structure
[elem2dof,edge] = dof3P2(elem);
N = size(node,1); NT = size(elem,1); NE = size(edge,1);

%% Add new nodes
node(N+1:N+NE,:) = (node(edge(:,1),:)+node(edge(:,2),:))/2;
HB(:,[1 2 3]) = [(N+1:N+NE)', edge(:,1:2)];

%% Refine each tetrahedron into 8 tetrahedrons
t = 1:NT;
p(t,1:10) = elem2dof;
elem(8*NT,:) = [0 0 0 0]; % enlarge the elem array
elem(t,:)           = [p(t,1), p(t,5), p(t,6), p(t,7)];
elem(NT+1:2*NT,:)   = [p(t,5), p(t,2), p(t,8), p(t,9)];
elem(2*NT+1:3*NT,:) = [p(t,6), p(t,8), p(t,3), p(t,10)];
elem(3*NT+1:4*NT,:) = [p(t,7), p(t,9), p(t,10), p(t,4)];
% After cutting off four subtetrahedra at the corners, the remaining
% octahedron can be subdivided in three different ways by choosing one of
% three possible interior diagonals: [6 9], [7 8] or [5 10]. For better
% mesh quality, we choose the shortest interior diagonal.
d1 = (node(p(t,7),1) - node(p(t,8),1)).^2 + ...
     (node(p(t,7),2) - node(p(t,8),2)).^2 + ...
     (node(p(t,7),3) - node(p(t,8),3)).^2;
d2 = (node(p(t,6),1) - node(p(t,9),1)).^2 + ...
     (node(p(t,6),2) - node(p(t,9),2)).^2 + ...
     (node(p(t,6),3) - node(p(t,9),3)).^2;
d3 = (node(p(t,5),1) - node(p(t,10),1)).^2 + ...
     (node(p(t,5),2) - node(p(t,10),2)).^2 + ...
     (node(p(t,5),3) - node(p(t,10),3)).^2;
[tempvar,minindex] = min([d1 d2 d3],[],2); %#ok<ASGLU>
% diagonal 7-8 is shorter than diagonal 6-9: use 7-8
idx = (minindex == 1);
elem(4*NT+t(idx),:) = [p(t(idx),6), p(t(idx),8), p(t(idx),7), p(t(idx),5)];
elem(5*NT+t(idx),:) = [p(t(idx),9), p(t(idx),7), p(t(idx),8), p(t(idx),5)];
elem(6*NT+t(idx),:) = [p(t(idx),6), p(t(idx),8), p(t(idx),10), p(t(idx),7)];
elem(7*NT+t(idx),:) = [p(t(idx),9), p(t(idx),7), p(t(idx),10), p(t(idx),8)];
% use diagonal 6-9
idx = (minindex == 2);
elem(4*NT+t(idx),:) = [p(t(idx),5), p(t(idx),6), p(t(idx),7), p(t(idx),9)];
elem(5*NT+t(idx),:) = [p(t(idx),5), p(t(idx),6), p(t(idx),9), p(t(idx),8)];
elem(6*NT+t(idx),:) = [p(t(idx),6), p(t(idx),7), p(t(idx),9), p(t(idx),10)];
elem(7*NT+t(idx),:) = [p(t(idx),6), p(t(idx),8), p(t(idx),10), p(t(idx),9)];
% use diagonal 5-10
idx = (minindex == 3);
elem(4*NT+t(idx),:) = [p(t(idx),6), p(t(idx),8), p(t(idx),10), p(t(idx),5)];
elem(5*NT+t(idx),:) = [p(t(idx),9), p(t(idx),7), p(t(idx),10), p(t(idx),5)];
elem(6*NT+t(idx),:) = [p(t(idx),6), p(t(idx),10), p(t(idx),5), p(t(idx),7)];
elem(7*NT+t(idx),:) = [p(t(idx),9), p(t(idx),10), p(t(idx),5), p(t(idx),8)];
% elem = label3(node,elem);

%% Update boundary edges
if (nargin>=3) && (~isempty(bdFlag))
    bdFlag(8*NT,:) = [0 0 0 0]; % enlarge the bdFlag array
    bdFlag(NT+1:2*NT,[1 3 4])   = bdFlag(t,[1 3 4]);
    bdFlag(2*NT+1:3*NT,[1 2 4]) = bdFlag(t,[1 2 4]);
    bdFlag(3*NT+1:4*NT,[1 2 3]) = bdFlag(t,[1 2 3]);
    idx = (minindex == 1);    % diagonal 7-8
    bdFlag(4*NT+t(idx),3) = bdFlag(t(idx),4);
    bdFlag(5*NT+t(idx),3) = bdFlag(t(idx),3);
    bdFlag(6*NT+t(idx),2) = bdFlag(t(idx),2);
    bdFlag(7*NT+t(idx),2) = bdFlag(t(idx),1);
    idx = (minindex == 2);   % diagonal 5-10
    bdFlag(4*NT+t(idx),2) = bdFlag(t(idx),3);
    bdFlag(5*NT+t(idx),3) = bdFlag(t(idx),4);
    bdFlag(6*NT+t(idx),3) = bdFlag(t(idx),2);
    bdFlag(7*NT+t(idx),1) = bdFlag(t(idx),1);
    idx = (minindex == 3);
    bdFlag(4*NT+t(idx),3) = bdFlag(t(idx),4);
    bdFlag(5*NT+t(idx),3) = bdFlag(t(idx),3);
    bdFlag(6*NT+t(idx),3) = bdFlag(t(idx),2);
    bdFlag(7*NT+t(idx),3) = bdFlag(t(idx),1);
    % change t in the last
    bdFlag(t,1) = 0;
else
    bdFlag = [];
end

repo: ifem


text:
function errorbar(x, y, l,u)
%ERRBAR Plot error bars.
%	ERRBAR(X,Y,L,U) plots errorbars on an existing graph of vector X
%	vs. vector Y with error bars specified by the vectors L and U.  The vectors
%	X,Y, L and U must be the same length.  If X,Y, L and U are matrices then
%	each column produces a separate line.  The error bars are each drawn a
%	distance of U(i) above and L(i) below the points in (X,Y) so that each bar
%	is L(i) + U(i) long.
%	ERRBAR(X,Y,L) puts error bars on a graph of X vs Y that are symmetric
%	about Y and are 2*L(i) long.
%
%	ERRBAR(Y,L) puts error bars [Y-L Y+L] on a graph of Y.
%
%	For example,
%
%	   x = 1:10;
%	   y = sin(x);
%	   e = std(y)*ones(size(x));
%	   plot(x,y);
%	   errbar(x,y,e)
%
%	draws symmetric error bars of unit standard deviation.

%	L. Shure 5-17-88, 10-1-91 B.A. Jones 4-5-93
%	Copyright (c) 1984-94 by The MathWorks, Inc.

hold on
if min(size(x))==1,
  npt = length(x);
  x = x(:);
  y = y(:);
  if nargin == 3,
			l = l(:);
		elseif nargin == 4 | nargin == 5
			l = l(:);
			u = u(:);
		end
else
  [npt,n] = size(x);
end

if nargin == 3,
	u = l;
end

if nargin == 2
  l = y;
		u = y;
  y = x;
  [m,n] = size(y);
  x(:) = [1:npt]'*ones(1,n);;
end
if isstr(x) | isstr(y) | isstr(l) | isstr(u)
	error('Arguments must be numeric.')
end

if any(size(x)~=size(y)) | any(size(x)~=size(l)) |  any(size(x)~=size(u)),
  error('The sizes of X, Y, L and U must be the same.');
end

tee = (max(x(:))-min(x(:)))/100;  % make tee .02 x-distance for error bars
xl = x - tee;
xr = x + tee;
n = size(y,2);

% Plot graph and bars
cax = newplot;
next = lower(get(cax,'NextPlot'));
% build up nan-separated vector for bars
xb = [];
yb = [];
nnan = nan*ones(1,n);
for i = 1:npt
    ytop = y(i,:) + u(i,:);
    ybot = y(i,:) - l(i,:);
    xb = [xb; x(i,:); x(i,:) ; nnan; xl(i,:);xr(i,:);nnan ;xl(i,:);xr(i,:);nnan];
    yb = [yb; ytop;ybot;nnan;ytop;ytop;nnan;ybot;ybot;nnan];
end

plot(xb,yb)
hold off

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function res=teslaLeastLogistic(A, B, y, b, lambda1, lambda2, opts)
%
%%
% Function teslaLeastLogistic:
%      Logistic Loss for ...... (to be added)
%
%% Problem
%
%   to be added
%% Input parameters:
%
%  A-         Matrix of size m x n
%               (Logistic)
%  B-         Matrix of size mm x n
%               (Least Squares)
%  y -        Response vector (of size m x 1)
%               (for A, discrete, 1, -1)
%  b -        Response vector (of size mm x 1)
%               (real value)
%  lambda1 -
%  lambda2 -
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x k
%  c-         The obtained intercept if size 1 x k
%  funVal-    Function value during iterations
%
%% Related papers
%
% [1]  Jun Liu, Lei Yuan, and Jieping Ye, An Efficient Algorithm for
%      a Class of Fused Lasso Problems, KDD, 2010.
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last February 5, 2010.
%
% Related functions:
%  ....
%%

%% Verify and initialize the parameters
%%
if (nargin <7)
    error('\n Inputs: A, B, y, b, lambda1,  lambda2, and opts.ind should be specified!\n');
end

[m,n]=size(A);
mm=size(B,1);

if (length(y) ~=m || length(b)~=mm)
    error('\n Check the length of y and b!\n');
end

if (lambda1<0 || lambda2 <0)
    error('\n lambda1 should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%%

% Initialize ind
if ~isfield(opts,'ind1') || ~isfield(opts,'ind2')
    error('\n In teslaLeastLogistic, .ind1 and ind2 should be specified');
else
    ind1=opts.ind1;
    k=length(ind1)-1;

    ind2=opts.ind2;

    if ind1(k+1)~=m || length(ind2)-1~=k || ind2(k+1)~=mm
        error('\n Check opts.ind1, opts.ind2');
    end
end

% %% Normalization
%
% % Please refer to sll_opts for the definitions of mu, nu and nFlag
% %
% % If .nFlag =1, the input matrix A is normalized to
% %                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
% %
% % If .nFlag =2, the input matrix A is normalized to
% %                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
% %
% % Such normalization is done implicitly
% %     This implicit normalization is suggested for the sparse matrix
% %                                    but not for the dense matrix
% %
%
% if (opts.nFlag~=0)
%     if (isfield(opts,'mu'))
%         mu=opts.mu;
%         if(size(mu,2)~=n)
%             error('\n Check the input .mu');
%         end
%     else
%         mu=mean(A,1);
%     end
%
%     if (opts.nFlag==1)
%         if (isfield(opts,'nu'))
%             nu=opts.nu;
%             if(size(nu,1)~=n)
%                 error('\n Check the input .nu!');
%             end
%         else
%             nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
%         end
%     else % .nFlag=2
%         if (isfield(opts,'nu'))
%             nu=opts.nu;
%             if(size(nu,1)~=m)
%                 error('\n Check the input .nu!');
%             end
%         else
%             nu=(sum(A.^2,2)/n).^(0.5);
%         end
%     end
%
%     ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
%     % If some values in nu is typically small, it might be that,
%     % the entries in a given row or column in A are all close to zero.
%     % For numerical stability, we set the corresponding value to 1.
% end
%
% if (~issparse(A)) && (opts.nFlag~=0)
%     fprintf('\n -----------------------------------------------------');
%     fprintf('\n The data is not sparse or not stored in sparse format');
%     fprintf('\n The code still works.');
%     fprintf('\n But we suggest you to normalize the data directly,');
%     fprintf('\n for achieving better efficiency.');
%     fprintf('\n -----------------------------------------------------');
% end
%
% %% Starting point initialization
%
% p_flag=(y==1);                  % the indices of the postive samples
%
% for i=1:k
%     ind_i=(ind(i)+1):ind(i+1);  % indices for the i-th group
%
%     m1(1,i)=sum(p_flag(ind_i));      % the total number of the positive samples
%     m2(1,i)=length(ind_i)-m1(1,i);   % the total number of the negative samples
% end
%
% % process the regularization parameter
% if (opts.rFlag==1) % lambda1 here is the scaling factor lying in [0,1]
%     if (lambda1<0 || lambda1>1)
%         error('\n opts.rFlag=1, and lambda1 should be in [0,1]');
%     end
%
%     % we compute ATb for computing lambda_max, when the input lambda1 is a ratio
%
%     p_flag=(y==1);                  % the indices of the postive samples
%
%     for i=1:k
%         ind_i=(ind(i)+1):ind(i+1);  % indices for the i-th group
%
%         b(ind_i,1)=p_flag(ind_i)*m1(i)/(m1(i)+m2(i))-...
%             (~p_flag(ind_i))*m2(i)/(m1(i)+m2(i));
%     end
%
%     ATb=zeros(n, k);
%     % compute AT b
%     for i=1:k
%         ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
%
%         if (opts.nFlag==0)
%             tt =A(ind_i,:)'*b(ind_i,1);
%         elseif (opts.nFlag==1)
%             tt= A(ind_i,:)'*b(ind_i,1) - sum(b(ind_i,1)) * mu';
%             tt=tt./nu(ind_i,1);
%         else
%             invNu=b(ind_i,1)./nu(ind_i,1);
%             tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
%         end
%
%         ATb(:,i)= tt;
%     end
%
%     q_bar=Inf;
%
%     lambda_max=0;
%     for i=1:n
%         lambda_max=max(lambda_max,...
%             norm(  ATb(i,:), q_bar) );
%     end
%
%     lambda_max=lambda_max / m;
%     lambda1=lambda1*lambda_max;
%
%     if isfield(opts,'rFlag2')
%         if (opts.rFlag2==1)
%             lambda2=lambda1*lambda_max;
%         end
%     end
% end
%
% % initialize a starting point
% if opts.init==2
%     x=zeros(n,k); c=zeros(1,k);
% else
%     if isfield(opts,'x0')
%         x=opts.x0;
%         if ( size(x,1)~=n && size(x,2)~=k )
%             error('\n Check the input .x0');
%         end
%     else
%         x=zeros(n,k);
%     end
%
%     if isfield(opts,'c0')
%         c=opts.c0;
%
%         if ( length(c)~=k )
%             error('\n Check the input .c0');
%         end
%     else
%         c=log(m1./m2);
%     end
% end

x=zeros(n,k); c=zeros(1,k);

Ax=zeros(m,1);      % m x 1
Bx=zeros(mm,1);   % mm x 1
% compute Ax: Ax_i= A_i * x_i
for i=1:k
    ind_i=(ind1(i)+1):ind1(i+1);     % indices for the i-th group
    m_i=ind1(i+1)-ind1(i);           % number of samples in the i-th group
    Ax(ind_i,1)=A(ind_i,:)* x(:,i);

    ind_i=(ind2(i)+1):ind2(i+1);     % indices for the i-th group
    m_i=ind2(i+1)-ind2(i);           % number of samples in the i-th group
    Bx(ind_i,1)=B(ind_i,:)* x(:,i);
end

%% The main program
% The Armijo Goldstein line search schemes + accelearted gradient descent

z0=zeros(n-1,k);

bFlag=0; % this flag tests whether the gradient step only changes a little

L=1/m; % the intial guess of the Lipschitz continuous gradient

% assign xp with x, and Axp with Ax
xp=x; Axp=Ax; xxp=zeros(n,k);
cp=c;         ccp=zeros(1,k);

alphap=0; alpha=1;

for iterStep=1:opts.maxIter
    % --------------------------- step 1 ---------------------------
    % compute search point s based on xp and x (with beta)
    beta=(alphap-1)/alpha;    s=x + beta* xxp;   sc=c + beta* ccp;

    % --------------------------- step 2 ---------------------------
    % line search for L and compute the new approximate solution x

    % compute the gradient (g) at s
    As=Ax + beta* (Ax-Axp);

    % aa= - diag(y) * (A * s + sc)
    vec_sc=zeros(m,1);
    for i=1:k
        ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
        vec_sc(ind_i,1)=sc(i);
    end
    aa=- y.*(As+ vec_sc);

    % fun_s is the logistic loss at the search point
    bb=max(aa,0);
    fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

    % compute prob=[p_1;p_2;...;p_m]
    prob=1./( 1+ exp(aa) );

    % b= - diag(y) * (1 - prob)
    b= -y.*(1-prob) / m;

    gc=zeros(1,k);                     % the gradient of c
    for i=1:k
        ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
        gc(1,k)=sum(b(ind_i));
    end

    % compute g= AT b, the gradient of x
    for i=1:k
        ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

        if (opts.nFlag==0)
            tt =A(ind_i,:)'*b(ind_i,1);
        elseif (opts.nFlag==1)
            tt= A(ind_i,:)'*b(ind_i,1) - sum(b(ind_i,1)) * mu';
            tt=tt./nu(ind_i,1);
        else
            invNu=b(ind_i,1)./nu(ind_i,1);
            tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
        end

        g(:,i)= tt;
    end

    % copy x and Ax to xp and Axp
    xp=x;    Axp=Ax;
    cp=c;

    while (1)
        % let s walk in a step in the antigradient of s to get v
        % and then do the Lq/L1-norm regularized projection
        v=s-g/L; c= sc- gc/L;

        [x, z, gap]=tesla_proj(v, z0, ...
            lambda1/L, lambda2/L, n, k,...
            1000, 1e-8, 1, 6);
        z0=z;

        v=x-s;  % the difference between the new approximate solution x
        % and the search point s

        % compute Ax: Ax_i= A_i * x_i
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
            m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

            if (opts.nFlag==0)
                Ax(ind_i,1)=A(ind_i,:)* x(:,i);
            elseif (opts.nFlag==1)
                invNu=x(:,i)./nu; mu_invNu=mu * invNu;
                Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
            else
                Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
                Ax(ind_i,1)=Ax./nu(ind_i,1);
            end
        end

        % aa= - diag(y) * (A * x + c)
        vec_sc=zeros(m,1);
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
            vec_sc(ind_i,1)=c(i);
        end
        aa=- y.*(Ax+ vec_sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

        r_sum=(norm(v,'fro')^2 + norm(c-sc,2)^2) / 2;
        l_sum=fun_x - fun_s - sum(sum(v.* g)) - (c-sc)* gc';

        if (r_sum <=1e-20)
            bFlag=1; % this shows that, the gradient step makes little improvement
            break;
        end

        % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
        %                           + L/2 * (<v,v> + <c-sc,c-sc> )
        if(l_sum <= r_sum * L)
            break;
        else
            L=max(2*L, l_sum/r_sum);
            % fprintf('\n L=%5.6f',L);
        end
    end

    % --------------------------- step 3 ---------------------------
    % update alpha and alphap, and check whether converge
    alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

    ValueL(iterStep)=L;

    xxp=x-xp;   ccp=c-cp;

    funVal(iterStep)=fun_x;

    funVal(iterStep)=funVal(iterStep)+ lambda1 * sum(abs(x(:))) + lambda2 * sum( sum( abs(x(:,1:(k-1)) - x(:,2:k) ) ) );

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    switch(opts.tFlag)
        case 0
            if iterStep>=2
                if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iterStep>=2
                if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                        opts.tol* funVal(iterStep-1))
                    break;
                end
            end
        case 2
            if ( funVal(iterStep)<= opts.tol)
                break;
            end
        case 3
            norm_xxp=sqrt(xxp'*xxp);
            if ( norm_xxp <=opts.tol)
                break;
            end
        case 4
            norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
            if ( norm_xxp <=opts.tol * max(norm_xp,1))
                break;
            end
        case 5
            if iterStep>=opts.maxIter
                break;
            end
    end
    res.gap(:,iterStep)=gap;
end
res.x=x;
res.c=c;
res.funVal=funVal;
res.ValueL=ValueL;

author: thomaskuestner


text:
function circle_segment_test14 ( )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_TEST14 demonstrates R_JACOBI.
%
%  Discussion:
%
%    R_JACOBI returns recursion coefficients ALPHA and BETA for rules
%    using a Jacobi type weight w(x) = (1-x)^A * (1+x)^B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference
%
%    Walter Gautschi,
%    Orthogonal Polynomials: Computation and Approximation,
%    Oxford, 2004,
%    ISBN: 0-19-850672-4,
%    LC: QA404.5 G3555.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_SEGMENT_TEST14\n' );
  fprintf ( 1, '  R_JACOBI computes recursion coefficients ALPHA and BETA\n' );
  fprintf ( 1, '  for a Jacobi weight w(x)=(1-x)^A * (1+x)^B.\n' );
%
%  Legendre rule.
%
  n = 10;
  a = 0.0;
  b = 0.0;
  ab = r_jacobi ( n, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Legendre\n' );
  fprintf ( 1, '  A = %g, B = %g\n', a, b );
  fprintf ( 1, '  Alpha   Beta\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', ab(i,1), ab(i,2) );
  end
%
%  Chebyshev Type 1 rule.
%
  n = 10;
  a = -0.5;
  b = -0.5;
  ab = r_jacobi ( n, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Chebyshev Type 1\n' );
  fprintf ( 1, '  A = %g, B = %g\n', a, b );
  fprintf ( 1, '  Alpha   Beta\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', ab(i,1), ab(i,2) );
  end
%
%  Chebyshev Type 2 rule.
%
  n = 10;
  a = 0.5;
  b = 0.5;
  ab = r_jacobi ( n, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Chebyshev Type 2\n' );
  fprintf ( 1, '  A = %g, B = %g\n', a, b );
  fprintf ( 1, '  Alpha   Beta\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', ab(i,1), ab(i,2) );
  end
%
%  General Jacobi.
%
  n = 10;
  a = 0.5;
  b = 1.5;
  ab = r_jacobi ( n, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  General Jacobi weight\n' );
  fprintf ( 1, '  A = %g, B = %g\n', a, b );
  fprintf ( 1, '  Alpha   Beta\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', ab(i,1), ab(i,2) );
  end

  return
end


text:
function browse_mdl(bpath)
% browse_mdl(bpath), Browse models
%
% A quick and dirty simulink model file browser
% Copies pictures of the simulink models files specified by bpath into a new word document with the model name
% If bpath is not specified it lists all the simulink files in the current directory.

if nargin < 1
    files=dir('*.mdl');
else
    [pathstr,name] = fileparts(bpath);
    if exist(bpath,'dir')
        name = [name '\*'];
    end
    ext = '.mdl';
    files=dir(fullfile(pathstr,[name ext]));
end

if ~isempty(files)
    whandle = actxserver('word.application');
    whandle.Visible = 1;
    Add(whandle.Documents);
    sel=whandle.Selection;
    set(sel.ParagraphFormat,'Alignment','wdAlignParagraphLeft');
    set(sel.PageSetup,'LeftMargin',28.35,'RightMargin',28.35,'BottomMargin',28.35*1.43);
    for index=1:size(files,1)
        title=files(index).name;
        eval(title(1:end-4));
        print('-dmeta',['-s', title(1:end-4)]);
        bdclose;
        close('all');
        TypeText(sel,title);
        set(sel.ParagraphFormat,'KeepWithNext',-1)
        TypeParagraph(sel);
        Paste(sel);
        set(sel.ParagraphFormat,'KeepWithNext',0)
        TypeParagraph(sel);
    end
    release(whandle);
end

return


text:
function phase = ph(params,freq)
% This function computes the phase of an analog allpass network with 2 legs and N sections.
% The phase is evalutated at the set of frequencies given in the vector: freq.
% function phase = ph(params,freq)
%
% Inputs:
% params (Nx2)	Matrix of parameters (1/RC for 2 legs and N sections)
% freq	(1xM)	frequency vector (in Hz)
%
% Output:
% phase	(2xM)	output phase vector in deg.

[N,L] = size(params);
[dummy,M] = size(freq);

temp = zeros(N,2*M);		% set array size
%temp(:) = 180 - (360/pi)*atan(2*pi*(1./(10.^params(:)))*freq);	% compute phase matrix
temp(:) = (360/pi)*atan(2*pi*(1./(10.^params(:)))*freq);	% compute phase matrix (ignor the fact that each stage inverts signal)
phase = zeros(2,M);			% set array size
phase(:) = sum(temp,1);		% sum up the phase over N cascaded sections
return


text:
function []=plot_sb_baselines(ix)
%PLOT_SB_BASELINES plot the small baselines in small_baselines.list
%Optional an input argument ix can be specified containing a vector of the
%small baseline interferograms to keep which will be plotted in the baseline plot.
%
%   Andy Hooper, June 2007
%
%   ======================================================================
%   09/2010 AH: Add option to plot in MERGED directory
%   09/2010 AH: For SMALL_BASELINES/MERGED don't plot dropped ifgs
%   12/2012 DB: Added meaning of ix to the syntax of the code
%   04/2013 DB: Command variable
%   03/2014 DB: Suppress command line output
%   01/2017 DB: Relax the name for SB such it could have a prefix or suffix
%   01/2017 DB: add fontsize
%   ======================================================================

% 05/2019 MG: Modified to plot colorize network based on interferogram standard
% deviation

if nargin <1
   ix=[];
end

fontsize = 16;

currdir=pwd;
dirs=strread(currdir,'%s','delimiter','/');

if ~isempty(findstr(dirs{end},'SMALL_BASELINES'))

    if strcmp(dirs{end},'SMALL_BASELINES')
        [a,b] = system(['\ls -d [1,2]* | sed ''' 's/_/ /''' ' > small_baselines.list']);
    else
        cd ../SMALL_BASELINES
        [a,b] = system(['\ls -d [1,2]* | sed ''' 's/_/ /''' ' > ' currdir filesep 'small_baselines.list']);
        cd(currdir)
    end
    load ../psver
    psname=['../ps',num2str(psver)];
    small_baseline_flag='y';

elseif strcmp(dirs{end},'MERGED')
    cd ../SMALL_BASELINES
    [a,b] = system(['\ls -d [1,2]* | sed ''' 's/_/ /''' ' > ../MERGED/small_baselines.list']);
    cd ../MERGED
    load ../psver
    psname=['../ps',num2str(psver)];
    small_baseline_flag='y';
else
    load psver
    psname=['ps',num2str(psver)];
    small_baseline_flag='n';
end

sb=load('small_baselines.list');
n_ifg=size(sb,1);
if small_baseline_flag=='y' & isempty(ix) & exist('./parms.mat','file')
    drop_ifg_index=getparm('drop_ifg_index');
    if ~isempty(drop_ifg_index)
       ix=setdiff([1:n_ifg],drop_ifg_index);
    end
end

if ~isempty(ix)
    sb=sb(ix,:);
else
    ix=1:size(sb,1);
end

ps=load(psname);

% loading the std information
    ifgstdname=['ifgstd',num2str(2)];
    if exist([ifgstdname,'.mat'],'file')
       stdin=load(ifgstdname);
       ifg_std=stdin.ifg_std;
       clear stdin
    else
    ifg_std=zeros(ps.n_ifg,1);
    end

    if ~isempty(ix)
    ifg_std = ifg_std(ix,:);
    else
    ix=1:size(ifg_std,1);
    end

n_ifg=size(sb,1);
[yyyymmdd,I,J]=unique(sb);
ifg_ix=reshape(J,n_ifg,2);
x=ifg_ix(:,1);
y=ifg_ix(:,2);

day=str2num(datestr(ps.day,'yyyymmdd'));
[B,I]=intersect(day,yyyymmdd);

x=I(x);
y=I(y);

%create colormap
%colors
C1 = [1 0 0]; %red
C2 =[1 1 1];  %white
C3 = [0 0 1];  %blue

invert_flag = 0;

%invert colormap
if invert_flag ==1
  CM1=[linspace(C1(1),C2(1),100)',linspace(C1(2),C2(2),100)', linspace(C1(3),C2(3),100)'];
  CM2=[linspace(C2(1),C3(1),100)',linspace(C2(2),C3(2),100)', linspace(C2(3),C3(3),100)'];
  CM = [CM1;CM2];
elseif invert_flag ==0
  CM1=[linspace(C2(1),C1(1),100)',linspace(C2(2),C1(2),100)', linspace(C2(3),C1(3),100)'];
  CM2=[linspace(C3(1),C2(1),100)',linspace(C3(2),C2(2),100)', linspace(C3(3),C2(3),100)'];
  CM = [CM2;CM1];
else
end

%cm = colormap('parula'); % take your pick (doc colormap)
%cm = CM;
%cm = interp1(linspace(min(ifg_std),max(ifg_std),length(cm)),cm,ifg_std); % map color to y values
%cm = uint8(cm'*255); % need a 4xN uint8 array
%cm(4,:) = 255; % last column is transparency

h_fig=figure('name','Baseline plot','position',[ 440   170   821   628]);
for i=1:length(x)
    hold on
    x1 = [ps.day(x(i)),ps.day(y(i))]';
    y1 = [ps.bperp(x(i)),ps.bperp(y(i))]';
    xx = [x1 x1];
    yy = [y1 y1];
    zz = zeros(size(xx));
    cc = [ifg_std(i) ifg_std(i)];
    l=surf(xx,yy,zz,[cc;cc],'EdgeColor','interp','FaceColor','none');
    %l=plot([ps.day(x(i)),ps.day(y(i))],[ps.bperp(x(i)),ps.bperp(y(i))],'-');
    text((ps.day(x(i))+ps.day(y(i)))/2,(ps.bperp(x(i))+ps.bperp(y(i)))/2,num2str(ix(i)),'fontsize',fontsize-7);
    set(l,'linewidth',2)
    %set(l.Edge,'ColorBinding','interpolated','ColorData',cm)
    %drawnow i
end

if length(ix) ~= 1
colormap(CM);
cb = colorbar;
ylabel(cb,'Ifg std [deg]');
else
end
shading flat
view(2)
p=plot(ps.day,ps.bperp,'ko');
set(p,'markersize',6,'linewidth',2)
hold off
box on
set(gca,'fontsize',fontsize)
datetick('x',12)
ylabel('B_{perp}','fontsize',fontsize)
try
    print(h_fig,'-dpng','baseline_plot.png')
catch
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% Copyright (C) 1993-2013, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for MATLAB (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% We will show how to acquire images from files, movies, the web and attached
% video cameras.

% We start by opening a movie file, a standard MPG format movie file
camera = Movie('traffic_sequence.mpg');
about camera
% and we see that the result is a Movie object

% The displayed value of this object
camera
% shows the frame dimension, frame rate, number of frames etc.

% To grab a frame from the movie
im = camera.grab();
% and then display it
idisp(im)

% We can display subsequent frames from the movie
idisp( camera.grab() );
idisp( camera.grab() );
idisp( camera.grab() );
idisp( camera.grab() );

% Finally we destroy the camera object which deallocates resources
clear camera

% We can grab an image from an internet webcam if we know its URL
camera = AxisWebCamera('http://wc2.dartmouth.edu')
about camera
% and we see that the result is an AxisWebCamera  object

% We can display a frame from that image source
idisp( camera.grab() );
clear camera

% We can grab an image from Google Maps
ev = EarthView();
% which does require you to have a Google key (available for free) and given
% in the environment GOOGLE_KEY
idisp( ev.grab(-27.475722,153.0285, 17) )
% Alternatively we can a geographic lookup service to find a place by name
% rather than latitude/longitude
idisp( ev.grab('brisbane', 14) )

% Finally, if you have a camera in, or attached to, your computer you can grab
% a live image frame.
% First create a VideoCamera object
camera = VideoCamera();
idisp( camera.grab() )
% and close the camera
clear camera

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [ y, j ] = j_carry_hebrew ( y, j )

%*****************************************************************************80
%
%% J_CARRY_HEBREW carries year-days to years in a Hebrew date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, a YJ date.
%
  while ( 1 )

    days = year_length_hebrew ( y );

    if ( j < days )
      break
    end

    j = j - days;
    y = y + 1;

  end

  return
end


text:
%
% Please cite the following article when using this source code:
%

%
% Author:
% - Mehrtash Harandi (mehrtash.harandi at gmail dot com)
%
% This file is provided without any warranty of
% fitness for any purpose. You can redistribute
% this file and/or modify it under the terms of
% the GNU General Public License (GPL) as published
% by the Free Software Foundation, either version 3
% of the License or (at your option) any later version.

clear;
clc;

Beta = 1e1;
SPARSE_SOLVER = 'SPAMS'; %options are CVX and SPAMS

load toy_data; %The toy data is indeed a split of Brodatz-16c dataset.

TestSet;
TrainSet;

%Generate Kernel
galleryKernel = exp(-Beta*Stein_Divergence(TrainSet.X,TrainSet.X));
probeKernel = exp(-Beta*Stein_Divergence(TestSet.X,TrainSet.X))';

%Normalize dictionary
KD = galleryKernel./repmat(sqrt(sum(galleryKernel.^2)),[size(galleryKernel,1) 1]);

%
KX = probeKernel./repmat(sqrt(sum(probeKernel.^2)),[size(probeKernel,1) 1]);
L1 = size(probeKernel,2);
if (strcmpi(SPARSE_SOLVER,'CVX'))
    SR_Lambda = 5e-1;
    A = KD;
    Number_Visual_Codes = size(A,2);
    scX = zeros(Number_Visual_Codes,L1);
    for tmpC1=1:L1
        if (mod(tmpC1,10) == 1)
            fprintf('Working on sample %d/%d\n',tmpC1,L1);
        end
        cvx_begin quiet
        variable tmp_Sparse_X(Number_Visual_Codes);
        minimize( -2*tmp_Sparse_X'*KX(:,tmpC1)+tmp_Sparse_X'*KD*tmp_Sparse_X+SR_Lambda*norm(tmp_Sparse_X,1));
        % if you are interested in positive sparse codes, uncomment the following two lines
        %             subject to
        %                 tmp_Sparse_X >= 0;
        cvx_end
        scX(:,tmpC1) = tmp_Sparse_X;
    end %endfor
elseif (strcmpi(SPARSE_SOLVER,'SPAMS'))
    SR_Lambda = 1e-2;
    [KD_U,KD_D,~] = svd(KD);
    A = diag(sqrt(diag(KD_D)))*KD_U';
    D_Inv = KD_U*diag(1./sqrt(diag(KD_D)));
    KX = D_Inv'*KX;
    scX = full(mexLasso(KX,A,struct('mode',2,'lambda',SR_Lambda,'lambda2',0)));
else
    error('Sparse solver is not defined. This version accepts CVX or SPAMS');
end

%% Classification
qLabel = Classify_SRC(A,TrainSet.y,scX,KX);

CRR = sum((qLabel - TestSet.y == 0))/L1;

fprintf('Clasiffication Ratio: %f\n', CRR);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% TEST_cartPoleDynamics.m
%
% A simple script to run a forward simulation of the cart-pole dynamics, to
% make sure that everything looks reasonable.
clc; clear;

p.m1 = 1;
p.m2 = 1;
p.g = 9.81;
p.l = 1;

z0 = [...
    0;              %Initial horizontal position
    (pi/180)*170;   %Initial angle of the pendulum (from stable config)
    0;              %Initial velocity of cart
    0];             %initial rotation rate of pendulum

%Simple PD controller
u = @(z)( -[ 10*z(1) + 2*z(3) ; 10*z(2) + 2*z(4) ]  );

tSpan = [0;2];   %Time span for simulation

dynFun = @(t,z)( cartPoleDynamics(z,u(z),p) );   %Function handle to give ode45
sol = ode45(dynFun,tSpan,z0);  %Run simulation

%Extract uniformly spaced samples via high-order interpolation
t = linspace(tSpan(1),tSpan(2),250);
z = deval(sol,t);

% Make some simple plots:
figure(1); clf;

subplot(2,2,1);
plot(t,z(1,:));
xlabel('t')
ylabel('x')

subplot(2,2,2);
plot(t,z(2,:));
xlabel('t')
ylabel('q')

subplot(2,2,3);
plot(t,z(3,:));
xlabel('t')
ylabel('dx')

subplot(2,2,4);
plot(t,z(4,:));
xlabel('t')
ylabel('dq')

author: MatthewPeterKelly

repo: dscTutorials


text:
function islands=bwislands(img)
%
% islands=bwislands(img)
%
% return the indices of non-zero elements in a 2D or 3D image
% grouped by connected regions in a cell array
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%	 img: a 2D or 3D array
% output:
%	 islands: a cell array, each cell records the indices
%		  of the non-zero elements in img for a connected
%		  region (or an island)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

img=logical(1-img);
idx=find(1-img(:));
islands={};

count=1;
while(length(idx))
        [I,J]=ind2sub(size(img),idx(1));
	imgnew=imfill(img,[I,J]);
	islands{count}=find(imgnew~=img);
	count=count+1;
	img=imgnew;
	idx=find(1-img(:));
end


text:
%This Matlab script can be used to reproduce Figure 4.11 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Number of BS antennas
M = 100;

%Define the pilot reuse factor
f = 2;

%Select the number of setups with random UE locations
nbrOfSetups = 10;

%Select the number of channel realizations per setup
nbrOfRealizations = 1000;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Noise figure at the BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the exact Gaussian local scattering model
accuracy = 1;

%Range of angular standard deviation in the local scattering model (in degrees)
ASDdegRange = [1e-4 5:5:50];

%Prepare to save simulation results
sumSE_MR_theorem = zeros(length(ASDdegRange),nbrOfSetups);
sumSE_MR_UatF_corollary = zeros(length(ASDdegRange),nbrOfSetups);
sumSE_MR_UatF_alt1 = zeros(length(ASDdegRange),nbrOfSetups);
sumSE_MR_UatF_alt2 = zeros(length(ASDdegRange),nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup and many ASDs
    [R,channelGaindB] = functionExampleSetup(L,K,M,accuracy,ASDdegRange);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Go through all ASDs
    for spr = 1:length(ASDdegRange)

        %Output simulation progress
        disp([num2str(spr) ' ASDs out of ' num2str(length(ASDdegRange))]);

        %Generate channel realizations with estimates and estimation error
        %correlation matrices
        [Hhat,C,tau_p,~,H] = functionChannelEstimates(R(:,:,:,:,:,spr),channelGainOverNoise,nbrOfRealizations,M,K,L,p,f);

        %Compute SE using Monte Carlo simulations
        [signal,interference,prelogFactor,SE_MR] = functionComputeSE_UL_MR_MC(Hhat,H,C,tau_c,tau_p,nbrOfRealizations,M,K,L,p);

        %Save results
        sumSE_MR_theorem(spr,n) = mean(sum(SE_MR,1));

        sumSE_MR_UatF_corollary(spr,n) = prelogFactor*mean(sum(log2(1+signal(:,:,1)./interference(:,:,1)),1));

        sumSE_MR_UatF_alt1(spr,n) = prelogFactor*mean(sum(log2(1+signal(:,:,2)./interference(:,:,2)),1));

        sumSE_MR_UatF_alt2(spr,n) = prelogFactor*mean(sum(log2(1+signal(:,:,3)./interference(:,:,3)),1));

        %Delete large matrices
        clear Hhat H C;

    end

    %Delete large matrices
    clear R;

end

%Compute the average sum SE over setups
sumSE_MR_theorem_opt = mean(sumSE_MR_theorem,2);
sumSE_MR_UatF_corollary_opt = mean(sumSE_MR_UatF_corollary,2);
sumSE_MR_UatF_alt1_opt = mean(sumSE_MR_UatF_alt1,2);
sumSE_MR_UatF_alt2_opt = mean(sumSE_MR_UatF_alt2,2);

%% Plot the simulation results
figure(1);
hold on; box on;

plot(ASDdegRange,sumSE_MR_theorem_opt,'rs--','LineWidth',1);
plot(ASDdegRange,sumSE_MR_UatF_alt2_opt,'b*-','LineWidth',1);
plot(ASDdegRange,sumSE_MR_UatF_alt1_opt,'kd:','LineWidth',1);
plot(ASDdegRange,sumSE_MR_UatF_corollary_opt,'kp--','LineWidth',1);

xlabel('ASD $$(\sigma_{\varphi})$$ [degree]','Interpreter','latex');
ylabel('Average sum SE [bit/s/Hz/cell]');
legend('Original bound (Th. 4.1)','UatF: h/||h||^2','UatF: h/||h||','UatF: h (Corr. 4.5)','Location','NorthEast');
ylim([0 50]);

author: emilbjornson

repo: massivemimobook


text:
function face = compute_delaunay(vertex)

% compute_delaunay - compute faces of a Delaunay triangulation
%
%   face = compute_delaunay(vertex);
%
%   Copyright (c) 2008 Gabriel Peyre

if size(vertex,1)>size(vertex,2)
    vertex = vertex';
end

face = delaunay(vertex(1,:), vertex(2,:))';

author: gpeyre

repo: matlab-toolboxes