text:
function [c,Ls]=dgt2(f,g1,p3,p4,p5,p6)
%DGT2  2-D Discrete Gabor transform
%   Usage: c=dgt2(f,g,a,M);
%          c=dgt2(f,g1,g2,[a1,a2],[M1,M2]);
%          c=dgt2(f,g1,g2,[a1,a2],[M1,M2],[L1,L2]);
%          [c,Ls]=dgt2(f,g1,g2,[a1,a2],[M1,M2]);
%          [c,Ls]=dgt2(f,g1,g2,[a1,a2],[M1,M2],[L1,L2]);
%
%   Input parameters:
%         f       : Input data, matrix.
%         g,g1,g2 : Window functions.
%         a,a1,a2 : Length of time shifts.
%         M,M1,M2 : Number of modulations.
%         L1,L2   : Length of transform to do
%
%   Output parameters:
%         c       : array of coefficients.
%         Ls      : Original size of input matrix.
%
%   `dgt2(f,g,a,M)` will calculate a separable two-dimensional discrete
%   Gabor transformation of the input signal *f* with respect to the window
%   *g* and parameters *a* and *M*.
%
%   For each dimension, the length of the transform will be the smallest
%   possible that is larger than the length of the signal along that dimension.
%   f will be appropriately zero-extended.
%
%   `dgt2(f,g,a,M,L)` computes a Gabor transform as above, but does
%   a transform of length *L* along each dimension. *f* will be cut or
%   zero-extended to length *L* before the transform is done.
%
%   `[c,Ls]=dgt2(f,g,a,M)` or `[c,Ls]=dgt2(f,g,a,M,L)` additionally returns
%   the length of the input signal *f*. This is handy for reconstruction::
%
%                [c,Ls]=dgt2(f,g,a,M);
%                fr=idgt2(c,gd,a,Ls);
%
%   will reconstruct the signal *f* no matter what the size of *f* is, provided
%   that *gd* is a dual window of *g*.
%
%   `dgt2(f,g1,g2,a,M)` makes it possible to use a different window along the
%   two dimensions.
%
%   The parameters *a*, *M*, *L* and *Ls* can also be vectors of length 2.
%   In this case the first element will be used for the first dimension
%   and the second element will be used for the second dimension.
%
%   The output *c* has *4* or *5* dimensions. The dimensions index the
%   following properties:
%
%   1. Number of translation along 1st dimension of input.
%
%   2. Number of channel along 1st dimension  of input
%
%   3. Number of translation along 2nd dimension of input.
%
%   4. Number of channel along 2nd dimension  of input
%
%   5. Plane number, corresponds to 3rd dimension of input.
%
%   See also:  dgt, idgt2, gabdual

complainif_argnonotinrange(nargin,4,6,mfilename);

L=[];

if prod(size(p3))>2
  % Two windows was specified.
  g2=p3;
  a=p4;
  M=p5;
  if nargin==6
    L=p6;
  end;
else
  g2=g1;
  a=p3;
  M=p4;
  if nargin==5
    L=p5;
  end;
end;

if isempty(L)
  L1=[];
  L2=[];
else
  L1=L(1);
  L2=L(2);
end;

% Expand 'a' and M if necessary to two elements
a=bsxfun(@times,a,[1 1]);
M=bsxfun(@times,M,[1 1]);

Ls=size(f);
Ls=Ls(1:2);

c=dgt(f,g1,a(1),M(1),L1);
c=dgt(c,g2,a(2),M(2),L2,'dim',3);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ ap, kpvt, info ] = sspfa ( ap, n )

%*****************************************************************************80
%
%% SSPFA factors a real symmetric matrix stored in packed form.
%
%  Discussion:
%
%    To solve A*X = B, follow SSPFA by SSPSL.
%
%    To compute inverse(A)*C, follow SSPFA by SSPSL.
%
%    To compute determinant(A), follow SSPFA by SSPDI.
%
%    To compute inertia(A), follow SSPFA by SSPDI.
%
%    To compute inverse(A), follow SSPFA by SSPDI.
%
%  Packed storage:
%
%    The following program segment will pack the upper triangle of a
%    symmetric matrix.
%
%      k = 0
%      do j = 1, n
%        do i = 1, j
%          k = k + 1
%          ap(k) = a(i,j)
%        end
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real AP(N*(N+1)/2), the packed form of a
%    symmetric matrix A.  The columns of the upper triangle are stored
%    sequentially in a one-dimensional array.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real AP(N*(N+1)/2), a block diagonal matrix and the
%    multipliers which were used to obtain it stored in
%    packed form.  The factorization can be written A = U*D*U' where U
%    is a product of permutation and unit upper triangular matrices, U'
%    is the transpose of U, and D is block diagonal with 1 by 1 and 2
%    by 2 blocks.
%
%    Output, integer KPVT(N), the pivot indices.
%
%    Output, integer INFO, error flag.
%    0, normal value.
%    K, if the K-th pivot block is singular.  This is not an error
%    condition for this subroutine, but it does indicate that SSPSL or
%    SSPDI may divide by zero if called.
%

%
%  ALPHA is used in choosing pivot block size.
%
  alpha = ( 1.0 + sqrt ( 17.0 ) ) / 8.0;

  info = 0;
%
%  Main loop on K, which goes from N to 1.
%
  k = n;
  ik = floor ( ( n * ( n - 1 ) ) / 2 );

  while ( 1 )
%
%  Leave the loop if K = 0 or K = 1.
%
    if ( k == 0 )
      break
    end

    if ( k == 1 )
      kpvt(1) = 1;
      if ( ap(1) == 0.0 )
        info = 1;
      end
      break
    end
%
%  This section of code determines the kind of elimination to be performed.
%  When it is completed, KSTEP will be set to the size of the pivot block,
%  and SWAP will be set to .true. if an interchange is required.
%
    km1 = k - 1;
    kk = ik + k;
    absakk = abs ( ap(kk) );
%
%  Determine the largest off-diagonal element in column K.
%
    imax = isamax ( k-1, ap(ik+1:ik+k-1), 1 );
    imk = ik + imax;
    colmax = abs ( ap(imk) );

    if ( alpha * colmax <= absakk )

      kstep = 1;
      swap = 0;
%
%  Determine the largest off-diagonal element in row IMAX.
%
    else

      rowmax = 0.0;
      imaxp1 = imax + 1;
      im = floor ( imax * ( imax - 1 ) ) / 2;
      imj = im + 2 * imax;

      for j = imaxp1 : k
        rowmax = max ( rowmax, abs ( ap(imj) ) );
        imj = imj + j;
      end

      if ( imax ~= 1 )
        jmax = isamax ( imax-1, ap(im+1:im+imax-1), 1 );
        jmim = jmax + im;
        rowmax = max ( rowmax, abs ( ap(jmim) ) );
      end

      imim = imax + im;

      if ( alpha * rowmax <= abs ( ap(imim) ) )
        kstep = 1;
        swap = 1;
      elseif ( alpha * colmax * ( colmax / rowmax ) <= absakk )
        kstep = 1;
        swap = 0;
      else
        kstep = 2;
        swap = ( imax ~= km1 );
      end

    end
%
%  Column K is zero.  Set INFO and iterate the loop.
%
    if ( max ( absakk, colmax ) == 0.0 )

      kpvt(k) = k;
      info = k;

    else

      if ( kstep ~= 2 )
%
%  1 x 1 pivot block.
%
        if ( swap )
%
%  Perform an interchange.
%
          temp(1:imax) = ap(im+1:im+imax);
          ap(im+1:im+imax) = ap(ik+1:ik+imax);
          ap(ik+1:ik+imax) = temp(1:imax);

          imj = ik + imax;

          for jj = imax : k
            j = k + imax - jj;
            jk = ik + j;
            t = ap(jk);
            ap(jk) = ap(imj);
            ap(imj) = t;
            imj = imj - ( j - 1 );
          end

        end
%
%  Perform the elimination.
%
        ij = ik - ( k - 1 );

        for jj = 1 : km1
          j = k - jj;
          jk = ik + j;
          mulk = -ap(jk) / ap(kk);
          t = mulk;
          ap(ij+1:ij+j) = saxpy ( j, t, ap(ik+1:ik+j), 1, ap(ij+1:ij+j), 1 );
          ijj = ij + j;
          ap(jk) = mulk;
          ij = ij - ( j - 1 );
        end
%
%  Set the pivot array.
%
        if ( swap ) then
          kpvt(k) = imax;
        else
          kpvt(k) = k;
        end

      else
%
%  2 x 2 pivot block.
%
        km1k = ik + k - 1;
        ikm1 = ik - ( k - 1 );
%
%  Perform an interchange.
%
        if ( swap )

          temp(1:imax) = ap(im+1:im+imax);
          ap(im+1:im+imax) = ap(ikm1+1:ikm1+imax);
          ap(ikm1+1:ikm1+imax) = temp(1:imax);

          imj = ikm1 + imax;

          for jj = imax : km1
            j = km1 + imax - jj;
            jkm1 = ikm1 + j;
            t = ap(jkm1);
            ap(jkm1) = ap(imj);
            ap(imj) = t;
            imj = imj - ( j - 1 );
          end

          t = ap(km1k);
          ap(km1k) = ap(imk);
          ap(imk) = t;

        end
%
%  Perform the elimination.
%
        if ( k-2 ~= 0 )

          ak = ap(kk) / ap(km1k);
          km1km1 = ikm1 + k - 1;
          akm1 = ap(km1km1) / ap(km1k);
          denom = 1.0 - ak * akm1;
          ij = ik - ( k - 1 ) - ( k - 2 );

          for jj = 1 : k-2

            j = km1 - jj;
            jk = ik + j;
            bk = ap(jk) / ap(km1k);
            jkm1 = ikm1 + j;
            bkm1 = ap(jkm1) / ap(km1k);
            mulk = ( akm1 * bk - bkm1 ) / denom;
            mulkm1 = ( ak * bkm1 - bk ) / denom;
            t = mulk;
            ap(ij+1:ij+j) = saxpy ( j, t, ap(ik+1:ik+j), 1, ap(ij+1:ij+j), 1 );
            t = mulkm1;
            ap(ij+1:ij+j) = saxpy ( j, t, ap(ikm1+1:ikm1+j), 1, ap(ij+1:ij+j), 1 );
            ap(jk) = mulk;
            ap(jkm1) = mulkm1;
            ijj = ij + j;
            ij = ij - ( j - 1 );
          end

        end
%
%  Set the pivot array.
%
        if ( swap )
          kpvt(k) = -imax;
        else
          kpvt(k) = 1 - k;
        end

        kpvt(k-1) = kpvt(k);

      end

    end

    ik = ik - ( k - 1 );
    if ( kstep == 2 )
      ik = ik - ( k - 2 );
    end

    k = k - kstep;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pts2 = grFaceToPolygon(varargin)
%GRFACETOPOLYGON Compute the polygon corresponding to a graph face
%
%   PTS2 = grFaceToPolygon(NODES, EDGES, FACES, INDF)
%   PTS2 = grFaceToPolygon(NODES, FACES, INDF)
%   Where NODES, EDGES, and FACES are internal data of graph, and INDF is
%   the index of the face to extract. The result is the (ordered) set of
%   points composing the face.
%
%
%   PTS2 = grFaceToPolygon(GRAPH, INDF)
%   use structure representation for graph. The structure GRAPH must
%   contain data for fields 'nodes' and 'faces'.
%
%   If several indices face indices are specified, result is a cell array
%   of polygons.
%
%   The number of columns of PTS2 is the same as for NODES.
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2005-11-30
% Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   27/07/2007: cleanup code

if length(varargin)==2
    % argument is a graph structure
    graph = varargin{1};
    nodes = graph.nodes;
    faces = graph.faces;
    indf  = varargin{2};

elseif length(varargin)==3
    % arguments are nodes, faces and indices
    nodes = varargin{1};
    faces = varargin{2};
    indf  = varargin{3};

elseif length(varargin)==4
    % arguments are nodes, edges, faces and indices, we forget edges
    nodes = varargin{1};
    faces = varargin{3};
    indf  = varargin{4};
end

if iscell(faces)
    % faces is a cell array
    if length(indf)==1
        face = faces{indf};
        pts2 = nodes(face, :);
    else
        pts2 = cell(length(indf), 1);
        for i=1:length(indf)
            face = faces{indf(i)};
            pts2{i} = nodes(face, :);
        end
    end
else
    % faces is an indices array: all faces have same number of vertices
    if length(indf)==1
        face = faces(indf, :);
        pts2 = nodes(face, :);
    else
        pts2 = cell(length(indf), 1);
        for i=1:length(indf)
            face = faces(indf(i), :);
            pts2{i} = nodes(face, :);
        end
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
classdef PTKVesselness < PTKPlugin
    % PTKVesselness. Plugin for detecting blood vessels
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKVesselness computes a mutiscale vesselness filter based on Frangi et
    %     al., 1998. "Multiscale Vessel Enhancement Filtering". The filter
    %     returns a value at each point which in some sense representes the
    %     probability of that point belonging to a blood vessel.
    %
    %     To reduce memory usage, the left and right lungs are filtered
    %     separately and each is further divided into subimages using the
    %     PTKImageDividerHessian function. This will compute the eigenvalues of
    %     the Hessian matrix for each subimage and use these to compute the
    %     vesselness using the PTKComputeVesselnessFromHessianeigenvalues
    %     function.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Vesselness'
        ToolTip = 'Shows the multiscale vesselness filter for detecting blood vessels'
        Category = 'Vessels'
        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)

        function results = RunPlugin(dataset, reporting)

            right_lung = dataset.GetResult('PTKGetRightLungROI');

            reporting.PushProgress;

            reporting.UpdateProgressStage(0, 2);
            vesselness_right = PTKVesselness.ComputeVesselness(right_lung, reporting, false);

            reporting.UpdateProgressStage(1, 2);
            left_lung = dataset.GetResult('PTKGetLeftLungROI');
            vesselness_left = PTKVesselness.ComputeVesselness(left_lung, reporting, true);

            reporting.PopProgress;

            left_and_right_lungs = dataset.GetResult('PTKLeftAndRightLungs');

            results = PTKCombineLeftAndRightImages(dataset.GetTemplateImage(PTKContext.LungROI), vesselness_left, vesselness_right, left_and_right_lungs);

            lung = dataset.GetResult('PTKLeftAndRightLungs');
            results.ChangeRawImage(results.RawImage.*single(lung.RawImage > 0));
            results.ImageType = PTKImageType.Scaled;
        end

        function results = GenerateImageFromResults(results, ~, ~)
            vesselness_raw = 3*uint8(results.RawImage > 5);
            results.ChangeRawImage(vesselness_raw);
            results.ImageType = PTKImageType.Colormap;
        end

    end

    methods (Static, Access = private)

        function vesselness = ComputeVesselness(image_data, reporting, is_left_lung)

            reporting.PushProgress;

            sigma_range = 0.5 : 0.5: 2;
            num_calculations = numel(sigma_range);
            vesselness = [];
            progress_index = 0;
            for sigma = sigma_range
                reporting.UpdateProgressStage(progress_index, num_calculations);
                progress_index = progress_index + 1;

                mask = [];
                vesselness_next = PTKImageDividerHessian(image_data.Copy, @PTKVesselness.ComputeVesselnessPartImage, mask, sigma, [], false, false, is_left_lung, reporting);
                vesselness_next.ChangeRawImage(100*vesselness_next.RawImage);
                if isempty(vesselness)
                    vesselness =  vesselness_next.Copy;
                else
                    vesselness.ChangeRawImage(max(vesselness.RawImage, vesselness_next.RawImage));
                end
            end

            reporting.PopProgress;

        end

        function vesselness_wrapper = ComputeVesselnessPartImage(hessian_eigs_wrapper, voxel_size)
            vesselness_wrapper = PTKComputeVesselnessFromHessianeigenvalues(hessian_eigs_wrapper, voxel_size);
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
filename='CantileverBeam_Triangle_Linear_Fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 50;
Vfrac_final = 0.3;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

monitoring_interval = 1;
optimalityInitial = 1e-3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function quadrule_test165 ( n, alpha )

%*****************************************************************************80
%
%% TEST165 compares GEN_LAGUERRE_EK_COMPUTE.
%
%  Discussion:
%
%    I used this test to generate tabular values of weights and
%    abscissas for generalized Gauss-Laguerre quadrature.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Input, real ALPHA, the parameter.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST165\n' );
  fprintf ( 1, '  GEN_LAGUERRE_EK_COMPUTE computes generalized Gauss-Laguerre rules;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The printed output of this test can be inserted into\n' );
  fprintf ( 1, '  a MATLAB program.\n' );

  [ x, w ] = gen_laguerre_ek_compute ( n, alpha );

  fprintf ( 1, '%\n' );
  fprintf ( 1, '%  Abscissas X and weights W for a generalized Gauss Laguerre rule\n' );
  fprintf ( 1, '%  of order = %d\n', n );
  fprintf ( 1, '%  with ALPHA = %f\n', alpha );
  fprintf ( 1, '%\n' );

  for i = 1 : n
    fprintf ( 1, '    x(%d) = %f;\n', i, x(i) );
  end
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '    w(%d) = %f;\n', i, w(i) );
  end

  return
end


text:
function mfprintf(fid, varargin)
% mfprintf(fid, varargin)
% As fprintf but allows fid to be a vector.
% For example:
% fid = fopen('foo.txt','a');
% mfprintf([fid 1], 'Hello %s\m','world')
% fclose(fid);
% prints 'Hello world' into file foo.txt and on the screen at the same
% time.

arrayfun(@(fid) fprintf(fid, varargin{:}), fid,'UniformOutput',0);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function submitWithConfiguration(conf)
  addpath('./lib/jsonlab');

  parts = parts(conf);

  fprintf('== Submitting solutions | %s...\n', conf.itemName);

  tokenFile = 'token.mat';
  if exist(tokenFile, 'file')
    load(tokenFile);
    [email token] = promptToken(email, token, tokenFile);
  else
    [email token] = promptToken('', '', tokenFile);
  end

  if isempty(token)
    fprintf('!! Submission Cancelled\n');
    return
  end

  try
    response = submitParts(conf, email, token, parts);
  catch
    e = lasterror();
    fprintf( ...
      '!! Submission failed: unexpected error: %s\n', ...
      e.message);
    fprintf('!! Please try again later.\n');
    return
  end

  if isfield(response, 'errorMessage')
    fprintf('!! Submission failed: %s\n', response.errorMessage);
  else
    showFeedback(parts, response);
    save(tokenFile, 'email', 'token');
  end
end

function [email token] = promptToken(email, existingToken, tokenFile)
  if (~isempty(email) && ~isempty(existingToken))
    prompt = sprintf( ...
      'Use token from last successful submission (%s)? (Y/n): ', ...
      email);
    reenter = input(prompt, 's');

    if (isempty(reenter) || reenter(1) == 'Y' || reenter(1) == 'y')
      token = existingToken;
      return;
    else
      delete(tokenFile);
    end
  end
  email = input('Login (email address): ', 's');
  token = input('Token: ', 's');
end

function isValid = isValidPartOptionIndex(partOptions, i)
  isValid = (~isempty(i)) && (1 <= i) && (i <= numel(partOptions));
end

function response = submitParts(conf, email, token, parts)
  body = makePostBody(conf, email, token, parts);
  submissionUrl = submissionUrl();
  params = {'jsonBody', body};
  responseBody = urlread(submissionUrl, 'post', params);
  response = loadjson(responseBody);
end

function body = makePostBody(conf, email, token, parts)
  bodyStruct.assignmentSlug = conf.assignmentSlug;
  bodyStruct.submitterEmail = email;
  bodyStruct.secret = token;
  bodyStruct.parts = makePartsStruct(conf, parts);

  opt.Compact = 1;
  body = savejson('', bodyStruct, opt);
end

function partsStruct = makePartsStruct(conf, parts)
  for part = parts
    partId = part{:}.id;
    fieldName = makeValidFieldName(partId);
    outputStruct.output = conf.output(partId);
    partsStruct.(fieldName) = outputStruct;
  end
end

function [parts] = parts(conf)
  parts = {};
  for partArray = conf.partArrays
    part.id = partArray{:}{1};
    part.sourceFiles = partArray{:}{2};
    part.name = partArray{:}{3};
    parts{end + 1} = part;
  end
end

function showFeedback(parts, response)
  fprintf('== \n');
  fprintf('== %43s | %9s | %-s\n', 'Part Name', 'Score', 'Feedback');
  fprintf('== %43s | %9s | %-s\n', '---------', '-----', '--------');
  for part = parts
    score = '';
    partFeedback = '';
    partFeedback = response.partFeedbacks.(makeValidFieldName(part{:}.id));
    partEvaluation = response.partEvaluations.(makeValidFieldName(part{:}.id));
    score = sprintf('%d / %3d', partEvaluation.score, partEvaluation.maxScore);
    fprintf('== %43s | %9s | %-s\n', part{:}.name, score, partFeedback);
  end
  evaluation = response.evaluation;
  totalScore = sprintf('%d / %d', evaluation.score, evaluation.maxScore);
  fprintf('==                                   --------------------------------\n');
  fprintf('== %43s | %9s | %-s\n', '', totalScore, '');
  fprintf('== \n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Service configuration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function submissionUrl = submissionUrl()
  submissionUrl = 'https://www-origin.coursera.org/api/onDemandProgrammingImmediateFormSubmissions.v1';
end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function[varargout]=trajchunk(varargin)
%TRAJCHUNK  Chunks Lagrangian trajectories based on the Coriolis period.
%
%   TRAJCHUNK is used to split float or drifter data into chunks such that
%   the length of each chunk is a fixed multiple of one over the mean
%   Coriolis frequency. This can be useful in spectral analysis.
%
%   [NUMO,LATO]=TRAJCHUNK(NUM,LAT,P), where NUM and LAT are date number and
%   latitude for Lagangian float or drifter data, re-organizes these into
%   chunks such that the average Coriolis frequency f_C in each chunk is at
%   least P times the Rayleigh frequency f_R for that chunk.
%
%   The Rayleigh frequency is f_R=2*pi/(DT*N), in units of radians per unit
%   time, where DT is the sample interval and N is the number of samples.
%   The Rayleigh frequency decreases as the chunk length increases.
%
%   The input fields NUM and LAT may either be numerical arrays, or
%   cell arrays of numerical arrays, e.g. NUM{1}=NUM1, NUM{2}=NUM2, etc.
%
%   The output variables NUM and LAT are cell arrays of numerical arrays,
%   with each cell being just long enough such that f_C > P * f_R.
%
%   Trajectories that are not long enough to satisfy this criterion are
%   discarded, as are short residual segments at the end of trajectories.
%
%   Each input trajectory is thus split into zero, one, or more than one
%   cells in the output variables.
%
%   TRAJCHUNK(...,P,LMIN) additionally specifies a mininum number of points
%   LMIN for each chunk.
%   __________________________________________________________________
%
%   Multiple input / output arguments
%
%   [NUMO,LATO,Y1,Y2,...,YM]=TRAJCHUNK(NUM,LAT,X1,X2,...XM,P) chunks the M
%   input arrays X1, X2,... XM in the same manner, and returns these as Y1,
%   Y2,... YM.  The input variables may either all be numerical arrays of
%   all the same size, or cell arrays of numerical arrays.
%
%   In the case of cell array input, some of the XM may be numerical arrays
%   of the same length as the cells NUM and LAT.  The corresponding output
%   variable will then also be a numerical array.  An example of such a
%   field is the identification number used in FLOATS.MAT and DRIFTERS.MAT.
%
%   TRAJCHUNK with no output arguments overwrites the original named output
%   variables.
%   __________________________________________________________________
%
%   Optional behaviors
%
%   By default, any data in short trajectories for which f_C < P * f_R are
%   discarded, as are data segments at the end of the trajectories.
%
%   TRAJCHUNK(...,'keep') keeps these instead.  Short trajectories are
%   returned in their own chunks, and leftover segments are appended to
%   the end of the preceding chunk.
%
%   TRAJCHUNK(...,'full') instead ensures that the output cells span the
%   full duration of the input fields.  This is done by appending a final
%   cell having f_C > P * f_R, like the others, but that ends at the final
%   data point, regardless of the degree of overlap with the previous cell.
%   Data from cells shorter than the specified length are discarded.
%   __________________________________________________________________
%
%   Overlap
%
%   TRAJCHUNK(...,'overlap',PCT) outputs chunks with a percentage PCT
%   overlap.  For example, TRAJCHUNK(...,'overlap',50) outputs chunks that
%   overlap by 50%.  The default behavior gives chunks with no overlap.
%   __________________________________________________________________
%
%   See also CELLCHUNK.
%
%   'trajchunk --t' runs a test.
%
%   Usage: [num,lat]=trajchunk(num,lat,P);
%          [num,lat,lon,cv]=trajchunk(num,lat,lon,cv,P);
%          [num,lat,lon,cv]=trajchunk(num,lat,lon,cv,P,lmin);
%          [num,lat,lon,cv]=trajchunk(num,lat,lon,cv,P,'overlap',50);
%          trajchunk(num,lat,lon,cv,P);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014--2019 J.M. Lilly --- type 'help jlab_license' for details

%   [...,II,KK]=TRAJCHUNK(...) in this case also outputs the indices of
%   the data locations within the input cells.  KK is not a cell array
%   like the other output arguments, but rather a row array of LENGTH(II).

%   [..,II]=TRAJCHUNK(...), with an extra final output argument,
%   outputs a cell array II of indices to the original time series.
%
%   As an example, LAT(II{1}) gives the latitudes of the data in the first
%   cell of the output, LATO{1}.

%   Equivalently, this means that the inertial period 2*pi/f_C is just less
%   than 1/M times the chunk duration DT*N, or 2*pi/f_C < (1/M) * (DT*N).
%   Thus M inertial oscillations fit into each chunk.
%   Not completely sure about the wording of this due to the definition of
%   'mean Coriolis frequency' etc.

if strcmpi(varargin{1}, '--t')
    trajchunk_test,return
end

factor=1;
opt='nokeep';     %What to do with leftover points

for i=1:2
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'kee')||strcmpi(varargin{end}(1:3),'nok')||strcmpi(varargin{end}(1:3),'ful')
            opt=varargin{end};
            varargin=varargin(1:end-1);
        end
    end
    if ischar(varargin{end-1})
        if strcmpi(varargin{end-1}(1:3),'ove')
            factor=1-varargin{end}/100;
            varargin=varargin(1:end-2);
        end
    end
end

if ~iscell(varargin{end-1})&&length(varargin{end-1})==1
    N=varargin{end-1};
    M=varargin{end};
    varargin=varargin(1:end-2);
else
    N=varargin{end};
    M=0;
    varargin=varargin(1:end-1);
end

num=varargin{1};
lat=varargin{2};
%Note: leave num and lat as first two entries also
na=length(varargin);

%na,N,M,size(num),size(lat),str,opt
if ~iscell(lat)
    %Create ii index
    varargin{na+1}=[1:length(lat)]';
    index=trajchunk_index(num,lat,N,M,opt,factor);
    n=0;
    if ~isempty(index)
        for k=1:length(index)
            n=n+1;
            for j=1:na+1
                varargout{j}{n,1}=varargin{j}(index{k});
            end
        end
    end
else
    %/**************
    %Put numerical array input into cell arrays
    bid=false(size(varargin));
    for i=3:length(bid) %Lat and lon are not allowed to be arrays
        if ~iscell(varargin{i})
            bid(i)=true;
            varargin{i}=celladd(varargin{i},cellmult(0,varargin{1}));
        end
    end
    %\**************

%     %Create ii and kk indices
%     for i=1:length(lat)
%         varargin{na+1}{i}=[1:length(lat{i})]';
%         varargin{na+2}{i}=i+0*lat{i};
%     end
    index=[];
    for i=1:length(lat)
        if length(lat)>1000
            if res((i-1)/1000)==0
                disp(['TRAJCHUNK working on cells ' int2str(i) ' to ' int2str(min(i+1000,length(lat))) ' of ' int2str(length(lat)) '.'])
            end
        end
        index{i,1}=trajchunk_index(num{i},lat{i},N,M,opt,factor);
    end
    n=0;
    for i=1:length(lat)
        if ~isempty(index{i})
            for k=1:length(index{i})
                n=n+1;
                for j=1:na
                    varargout{j}{n,1}=varargin{j}{i}(index{i}{k});
                end
            end
        end
    end

    %Return numerical array input back to numerical arrays
    for i=1:length(bid)
        if bid(i)
            varargout{i}=cellfirst(varargout{i});
        end
    end

%     %i,j,k,n
%     %Convert cell array kk into numeric array
%     temp=varargout{na+2};
%     varargout{na+2}=zeros(size(varargout{na+1}));
%     for i=1:length(varargout{na+1})
%         varargout{na+2}(i)=temp{i}(1);
%     end
end

eval(to_overwrite(na));

function[index]=trajchunk_index(num,lat,N,M,opt,factor)
if length(num)<2
    index=[];
else
    dt=num(2)-num(1);

    fo=abs(corfreq(lat))*24;
    meanfo=frac(cumsum(fo),[1:length(fo)]');
    fr=frac(2*pi,dt*[1:length(fo)]');
    %aresame(meanfo(end),vmean(fo,1))

    bdone=false;
    n=0;

    x=[1:length(fo)]';
    index=[];
    while ~bdone
        ii=max(find(N*fr<meanfo,1,'first'),M);
        n=n+1;
        if isempty(ii)||ii>=length(x)
            bdone=1;
        else
            index{n}=x(1:ii);
            %N*fr(ii)<meanfo(ii)
            if ii<length(x)
                x=x(floor(ii*factor)+1:end);
                fo=fo(floor(ii*factor)+1:end);
                %fr=fr(ii+1:end)-fr(ii+1)+frac(2*pi,dt);
                meanfo=frac(cumsum(fo),[1:length(fo)]');
                fr=frac(2*pi,dt*[1:length(fo)]');
                %aresame(meanfo(end),vmean(fo,1))
            end
        end
    end

    % if ~isempty(index)
    %     if strcmpi(opt(1:3),'kee')
    %         index{end}=index{end}(1):x(end);
    %     end
    % end
    if strcmpi(opt(1:3),'kee')
        if ~isempty(index)
            index{end}=index{end}(1):x(end);  %Append leftovers
        else
            index{1}=x;  %Keep short segments
        end
    elseif strcmpi(opt(1:3),'ful')
        indexlast=trajchunk_index(num,flipud(lat),N,M,'nokeep',1);
        indexlast=length(lat)-flipud(indexlast{1})+1;
        index{end+1}=indexlast;
        %maxmax(indexlast),length(lat)
    end
end
function[]=trajchunk_test

load ebasnfloats

use ebasnfloats
dt=1;
trajchunk(num,lat,lon,32);
meanfo=vmean(abs(corfreq(col2mat(cell2col(lat)))),1)'*24*dt;
fr=frac(2*pi,dt*cellength(lat));
reporttest('TRAJCHUNK no overlap',allall(meanfo>32*fr))
%length(num),length(cell2col(num))

use ebasnfloats
dt=1;
trajchunk(num,lat,lon,32,'overlap',50);
meanfo=vmean(abs(corfreq(col2mat(cell2col(lat)))),1)'*24*dt;
fr=frac(2*pi,dt*cellength(lat));
reporttest('TRAJCHUNK overlap',allall(meanfo>32*fr))

% load drifterheyerdahl
% use drifterheyerdahl
%
% [numc,latc]=trajchunk(num,lat,60,'overlap',50);
% [~,latc2,numc2]=trajchunk(num,flipud(lat),flipud(num),60,'overlap',50);
%
% min(cellmin(numc))-min(num)
% max(cellmax(numc))-max(num)
%
% numc{end+1}=flipud(numc2{1});
% latc{end+1}=flipud(latc2{1});
%
% for i=1:length(numc)
%     P(i)=sum(abs(corfreq(latc{i}))*24.*(numc{i}-numc{i}(1))/2/pi/24);
% end
%
% plot(cellmax(numc)-cellmin(numc))
% plot(60*2*pi./(cellmean(cellabs(corfreq(latc)))*24))

%length(num),length(cell2col(num))

% with minimum length
%
%     trajchunk(num,lat,lon,32,35);
%
% load drifters
% use drifters
% %trajchunk(num,lat,lon,32,'overlap');
% trajchunk(num,lat,lon,128);
%
% %meanfo=zeros(length(num),1);
% %fr=zeros(length(num),1);
%
% %dt=1/4;
%
% %Same as below but faster
% tic
% meanfo=vmean(abs(corfreq(col2mat(cell2col(lat)))),1)'*24*dt;
% fr=frac(2*pi,dt*cellength(lat));
% toc
%
% % tic
% % for i=1:length(num)
% %     meanfo(i)=vmean(abs(corfreq(lat{i})),1)*24*dt;
% %     fr(i)=frac(2*pi,dt*length(lat{i}));
% % end
% % toc
%
% figure,plot(meanfo,'b.'),hold on,plot(32*fr,'ro')
% figure,plot(2*pi./(meanfo/32),'b.'),hold on,plot(2*pi./fr,'ro')
% figure,plot(2*pi./(meanfo/32),2*pi./fr,'ro'),axis equal

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function out = gtdata2mhaconfig(gtdata)
  % This function writes the gain table data into the correct layout to be fed
  % to the mhaconfig file.
  %
  % Variables:
  % - gtdata
  %
  % Returns:
  % - out

  out = '[';
  for i=1:size(gtdata, 1)
    out = [out '[' sprintf('%.1f ',gtdata(i,:)) '];'];
  end
  out = [out ']'];
end

author: m-r-s

repo: hearingaid-prototype


text:
function path = funpath(~)
%FUNPATH Returns the path to the calling function.
% Usage:
%   path = funpath()
%   path = funpath(_) % returns the path to the parent directory
%
% See also: get_caller_name, which, dbstack

% Get the path to the calling function
path = get_caller_name('path', true);

% Return the path to the function's parent directory
if nargin > 0
    path = fileparts(path);
end
end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function owens_test03 ( )

%*****************************************************************************80
%
%% TEST03 demonstrates the use of ZNORM1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03:\n' );
  fprintf ( 1, '  ZNORM1 computes the normal CDF starting at 0.\n' );
  fprintf ( 1, '  Compare to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          X           P                         P                       DIFF\n' );
  fprintf ( 1, '                     (Tabulated)               (ZNORM1)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = normal_01_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx1 = fx1 - 0.5;

    fx2 = znorm1 ( x );

    fprintf ( 1, '  %12.8f  %24.16e  %24.16e  %10.4e\n', ...
    x, fx1, fx2, abs ( fx1 - fx2 ) );

  end

  return
end


text:
function [inference] = tapas_sem_hier_inference(inference, pars)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

if ~isfield(inference, 'estimate_method')
    inference.estimate_method = @tapas_mcmc_blocked_estimate;
end

if ~isfield(inference, 'initialize_states')
    inference.initialize_states = @tapas_sem_hier_init_states;
end

if ~isfield(inference, 'initialize_state')
    inference.initialize_state = @tapas_sem_hier_init_state;
end

if ~isfield(inference, 'sampling_methods')
    inference.sampling_methods = {
        @(d, m, i, s) tapas_mh_mc3_tempering_sample_node(d, m, i, s, 2), ...
        @(d, m, i, s) tapas_sampler_dlinear_gibbs_node(d, m, i, s, 3), ...
        @tapas_sampler_mc3_tempering ...
        ...  Use population mcmc step with generalized
        % Temperature for possibly Bayesian predictive distribution.
        };
end

if ~isfield(inference, 'metasampling_methods')
    inference.metasampling_methods = {@tapas_mcmc_meta_diagnostics, ...
        @tapas_mcmc_meta_adaptive, ...
%        @tapas_mcmc_meta_adaptive_ti, ...
        };
end

if ~isfield(inference, 'get_stored_state')
    inference.get_stored_state = @tapas_sem_hier_get_stored_state;
end

if ~isfield(inference, 'prepare_posterior')
    inference.prepare_posterior = @tapas_sem_hier_prepare_posterior;
end

if ~isfield(inference, 'mh_sampler')
    inference.mh_sampler = cell(4, 1);
end

if ~isfield(inference.mh_sampler{2}, 'propose_sample')
    inference.mh_sampler{2}.propose_sample = ...
        @tapas_mh_mc3_propose_gaussian_sample;
    inference.mh_sampler{2}.ar_rule = ...
        @tapas_mh_mc3g_arc;
end

if ~isfield(inference, 'niter')
    if isfield(pars, 'niter')
        inference.niter = pars.niter;
    else
        inference.niter = 5000;
    end
end

if ~isfield(inference, 'nburnin')
    if isfield(pars, 'nburnin')
        inference.nburnin = pars.nburnin;
    else
        inference.nburnin = 5000;
    end
end

if ~isfield(inference, 'mc3it')
    if isfield(pars, 'mc3it')
        inference.mc3it = pars.mc3it;
    else
        inference.mc3it = 10;
    end
end

if ~isfield(inference, 'thinning')
    if isfield(pars, 'thinning')
        inference.thinning = pars.thinning;
    else
        inference.thinning = 0;
    end
end

if ~isfield(inference, 'ndiag')
    if isfield(pars, 'ndiag')
        inference.ndiag = pars.ndiag;
    else
        inference.ndiag = 200;
    end
end

% Scale of the kernel
if ~isfield(inference, 'kernel_scale')
    inference.kernel_scale = 0.1 * 0.1;
end

end

author: translationalneuromodeling

repo: tapas


text:
function Edges = detectAreaBorders(images)

images(isnan(images))=0;
% Edges=double(edge(images, 'canny'));
images = images+1;
Edges = zeros(size(images));
for area=3:max(images(:))-1
    images_tmp = double(images==area);
    Edges=Edges | double(edge(images_tmp, 'canny'));
end
% figure;imagesc(Edges)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function removeLayer(obj, name, chainInputs)
%REMOVELAYER Remove a layer from the network
%   REMOVELAYER(OBJ, NAME) removes the layer NAME from the DagNN object
%   OBJ.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin < 3, chainInputs = true; end
f = find(strcmp(name, {obj.layers.name})) ;
if isempty(f), error('There is no layer ''%s''.', name), end
layer = obj.layers(f) ;

if chainInputs
  % chain input of l that has layer as input
  for l = 1:numel(obj.layers)
    for i = layer.outputIndexes
      sel = find(intersect(obj.layers(l).inputIndexes, i));
      if any(sel)
        obj.layers(l).inputs{sel} = layer.inputs{1};
      end
    end
  end
end
obj.layers(f) = [] ;
obj.rebuild() ;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%IM_BOX Fixed mapping determining a rectangular enclosing a blob (0/1 image)
%
%   B = IM_BOX(A)
%   B = A*IM_BOX
%
% If A is a 0/1 image then B is the same image with all empty (0) border
% columns and rows removed.
%
%   B = IM_BOX(A,N)
%
% If A is a 0/1 image then B is the same image, but having in each direction
% N empty (0) columns and rows around the object (1).
%
%   B = IM_BOX(A,[NX1 NX2 NY1 NY2])
%
% If A is a 0/1 image then B is the same image, but having NX1, NX2 empty
% columns (0) left, respectively right of the object (1) and NY1, NY2 empty
% rows (0) above, respectively below the object(1).
%
%   B = IM_BOX(A,N,ALF)
%
% Adds as many empty (0) columns or rows such that the aspect ratio of
% images (height/width) equals ALF. For ALF == 1, square images are returned.
% For ALF == 0, images are taken as they are and N rows and columns are
% added.
%
% EXAMPLE
% prdatafiles;            % make sure prdatafiles is in the path
% x = kimia_images;       % load kimia images
% x = x*im_box(0);        % remove all empty rows and columns
% x = x*im_box(0,1);      % add rows/columns to make images square
% x = x*im_resize([32 32]); % resample them
% x = x*im_box(1,0);      % add rows/columns and keep image square
% % now all images are 34x34 and no object touches the border.
% show(gendat(x,4))       % show a few
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function [b,J] = im_box(varargin)

	argin = shiftargin(varargin,'vector');
  argin = setdefaults(argin,[],[],[]);
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'Image bounding box');
  else
    [a,n,alf] = deal(argin{:});
    if isdataset(a)
      error('Command cannot be used for datasets as it may change image size')
    elseif isdatafile(a)
      isobjim(a);
      b = filtim(a,mfilename,{n,alf});
      %b = setfeatsize(b,getfeatsize(a));
    elseif isa(a,'double') || isa(a,'dip_image') % here we have a single image
      if isa(a,'dip_image'), a = double(a); end
      if isempty(n)
        jx = find(sum(a,1) ~= 0);
        jy = find(sum(a,2) ~= 0);
        J = [min(jy),max(jy),min(jx),max(jx)];
        b = a(min(jy):max(jy),min(jx):max(jx));
      else
        if (~isempty(alf) && alf == 0)
          c = a;
        else
          c = feval(mfilename,a);
        end
        [my,mx] = size(c);
        if length(n) == 1
        	n = [n n n n];
        elseif length(n) ~= 4
        	error('Second parameter should be scalar or vector of length 4')
        end
        b = zeros(my+n(3)+n(4),mx+n(1)+n(2));
        b(n(3)+1:n(3)+my,n(1)+1:n(1)+mx) = c;
      end
      if ~isempty(alf) && alf ~= 0
        [m,k] = size(b);
        r = round(m*alf) - k;
        if r == 0
          ;
        elseif r >= 1 % add r columns
          c = zeros(m,k+r);
          c(:,ceil(r/2):ceil(r/2)+k-1) = b;
          b = c;
        else % add rows
          r = round(k/alf) - m;
          c = zeros(m+r,k);
          c(ceil(r/2):ceil(r/2)+m-1,:) = b;
          b = c;
        end
      end
		end
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% t_meshFromFreesurfer
%
% Illustrates how to create mrVista compatible meshes from a freesurfer
% surface.
%
%
% Dependencies:
%   Remote Data Toolbox
%   Freesurfer paths
%
% This tutorial is part of a sequence. Run
%   t_initAnatomyFromFreesurfer
% prior to running this tutorial.
%
% Summary
%
% - Do actiom A
% - Do action B
% - Visualize
% - Clean up
%
% Tested 09/14/2016 - MATLAB r2015b, Mac OS 10.11.6
%
%  See t_initAnatomyFromFreesurfer
%
% NYU Winawer lab

%% Download ernie freesurfer directory
% Note that this is also done in the prior tutorial,
% t_initAnatomyFromFreesurfer. Doing it a second time won't hurt anyone,
% and it will make this tutorial work on its own, even without running the
% prior one in the sequence.

% Check whether freesurfer paths exist
fssubjectsdir = getenv('SUBJECTS_DIR');
if isempty(fssubjectsdir)
    error('Freesurfer paths not found. Cannot proceed.')
end

% Get ernie freesufer1 directory and install it in freesurfer subjects dir
%   If we find the directory, do not bother unzipping again
forceOverwrite = false;

% Do it
dFolder = mrtInstallSampleData('anatomy/freesurfer', 'ernie', ...
    fssubjectsdir, forceOverwrite);

fprintf('Freesurfer directory for ernie installed here:\n %s\n', dFolder)

%% Create Freesurfer meshes

% mrVista Project directory
erniePath      = fullfile(vistaRootPath, 'local', 'scratch', 'erniePRF');

if ~exist(erniePath, 'dir')
    % If we did not find the temporary directory created by prior
    % tutorials, then use the full directory, downloading if necessary
    erniePath = mrtInstallSampleData('functional', 'erniePRF');
end

% Install Ernie freesurfer directory and ernie vistasession
cd(erniePath);

% Create and save meshes
hemi = 'b';
surfaces = {'white' 'pial' 'sphere' 'inflated'};
[meshes, fnames] = meshImportFreesurferSurfaces('ernie', hemi, surfaces);

%% Visualize

% View the 4 left meshes: white, pial, inflated, sphere

for ii = 1:4
    meshVisualize(meshes(ii))
    % This call should open a new application, mrMeshMac or a related name
    % (depending on your OS), if it is not already open. Then it will open a
    % new window and show the mesh. You can rotate, zoom, and translate with
    % the appropriate mouse actions. We BELIEVE that the mouse buttons are as
    % follows:
    %   Rotation:  left mouse button plus movement
    %   Zoom:      right mouse button plus up/down movement
    %   Translate: left and right together plus movement

end

%% Alernative visualiztion in Matlab

for m = 1:4

    figure,

    % Faces (also called triangles) are defined by 3 points, each of
    % which is an index into the x, y, z vertices
    faces = meshes(m).triangles' + 1; % we need to 1-index rather than 0-index for Matlab

    % The vertices are the locations in mm spacing
    x     = meshes(m).vertices(1,:)';
    y     = meshes(m).vertices(2,:)';
    z     = meshes(m).vertices(3,:)';

    % The colormap will, by default, paint sulci dark and gyri light
    c     = meshes(m).colors(1,:)';

    % Render the triangle mesh
    tH = trimesh(faces, x,y,z);

    % Make it look nice
    set(tH, 'LineStyle', 'none', 'FaceColor', 'interp', 'FaceVertexCData',c)
    axis equal off; colormap gray; set(gca, 'CLim', [0 255])

    % Lighting to make it look glossy
    light('Position',100*[0 1 1],'Style','local')
    lighting gouraud

    % Which mesh are we plotting?
    title(meshes(m).name)

    % Rotate it
    set(gca, 'View', [-16.7000  -90.0000]);

end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Some examples using the tools in this distribution.
%%% Eero Simoncelli, 2/97.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Load an image, and downsample to a size appropriate for the machine speed.
oim = pgmRead('einstein.pgm');
tic; corrDn(oim,[1 1; 1 1]/4,'reflect1',[2 2]); time = toc;
imSubSample = min(max(floor(log2(time)/2+3),0),2);
im = blurDn(oim, imSubSample,'qmf9');
clear oim;

%%% ShowIm:
%% 3 types of automatic graylevel scaling, 2 types of automatic
%% sizing, with or without title and Range information.
help showIm
clf; showIm(im,'auto1','auto','Al')
clf; showIm('im','auto2')
clf; showIm(im,'auto3',2)

%%% Statistics:
mean2(im)
var2(im)
skew2(im)
kurt2(im)
entropy2(im)
imStats(im)

%%% Synthetic images.  First pick some parameters:
sz = 200;
dir = 2*pi*rand(1)
slope = 10*rand(1)-5
int = 10*rand(1)-5;
orig = round(1+(sz-1)*rand(2,1));
expt = 0.8+rand(1)
ampl = 1+5*rand(1)
ph = 2*pi*rand(1)
per = 20
twidth = 7

clf;
showIm(mkRamp(sz,dir,slope,int,orig));
showIm(mkImpulse(sz,orig,ampl));
showIm(mkR(sz,expt,orig));
showIm(mkAngle(sz,dir));
showIm(mkDisc(sz,sz/4,orig,twidth));
showIm(mkGaussian(sz,(sz/6)^2,orig,ampl));
showIm(mkZonePlate(sz,ampl,ph));
showIm(mkAngularSine(sz,3,ampl,ph,orig));
showIm(mkSine(sz,per,dir,ampl,ph,orig));
showIm(mkSquare(sz,per,dir,ampl,ph,orig,twidth));
showIm(mkFract(sz,expt));

%%% Point operations (lookup tables):
[Xtbl,Ytbl] = rcosFn(20, 25, [-1 1]);
plot(Xtbl,Ytbl);
showIm(pointOp(mkR(100,1,[70,30]), Ytbl, Xtbl(1), Xtbl(2)-Xtbl(1), 0));

%%% histogram Modification/matching:
[N,X] = histo(im, 150);
[mn, mx] = range2(im);
matched = histoMatch(rand(size(im)), N, X);
showIm(im + sqrt(-1)*matched);
[Nm,Xm] = histo(matched,150);
nextFig(2,1);
  subplot(1,2,1); plot(X,N); axis([mn mx 0 max(N)]);
  subplot(1,2,2);  plot(Xm,Nm); axis([mn mx 0 max(N)]);
nextFig(2,-1);

%%% Convolution routines:

%% Compare speed of convolution/downsampling routines:
noise = rand(400); filt = rand(10);
tic; res1 = corrDn(noise,filt(10:-1:1,10:-1:1),'reflect1',[2 2]); toc;
tic; ires = rconv2(noise,filt); res2 = ires(1:2:400,1:2:400); toc;
imStats(res1,res2)

%% Display image and extension of left and top boundaries:
fsz = [9 9];
fmid = ceil((fsz+1)/2);
imsz = [16 16];

% pick one:
im = eye(imsz);
im = mkRamp(imsz,pi/6);
im = mkSquare(imsz,6,pi/6);

% pick one:
edges='reflect1';
edges='reflect2';
edges='repeat';
edges='extend';
edges='zero';
edges='circular';
edges='dont-compute';

filt = mkImpulse(fsz,[1 1]);
showIm(corrDn(im,filt,edges));
line([0,0,imsz(2),imsz(2),0]+fmid(2)-0.5, ...
     [0,imsz(1),imsz(1),0,0]+fmid(1)-0.5);
title(sprintf('Edges = %s',edges));

%%% Multi-scale pyramids (see pyramids.m for more examples,
%%% and explanations):

%% A Laplacian pyramid:
[pyr,pind] = buildLpyr(im);
showLpyr(pyr,pind);

res = reconLpyr(pyr, pind); 		% full reconstruction
imStats(im,res);			% essentially perfect

res = reconLpyr(pyr, pind, [2 3]);  %reconstruct 2nd and 3rd levels only
showIm(res);

%% Wavelet/QMF pyramids:
filt = 'qmf9'; edges = 'reflect1';
filt = 'haar'; edges = 'qreflect2';
filt = 'qmf12'; edges = 'qreflect2';
filt = 'daub3'; edges = 'circular';

%[pyr,pind] = buildWpyr(im, 5-imSubSample, filt, edges);
[pyr,pind] = buildWpyr(im, 'auto', filt, edges);
showWpyr(pyr,pind,'auto2');

res = reconWpyr(pyr, pind, filt, edges);
clf;
showIm(im + i*res);
imStats(im,res);

res = reconWpyr(pyr, pind, filt, edges, 'all', [2]);  %vertical only
clf;
showIm(res);

%% Steerable pyramid:
%[pyr,pind] = buildSpyr(im,4-imSubSample,'sp3Filters');
[pyr,pind] = buildSpyr(im,'auto','sp3Filters');
showSpyr(pyr,pind);

%% Steerable pyramid, constructed in frequency domain:
%[pyr,pind] = buildSFpyr(im,5-imSubSample,4);  %5 orientation bands
[pyr,pind] = buildSFpyr(im);  %5 orientation bands
showSpyr(pyr,pind);
res = reconSFpyr(pyr,pind);
imStats(im,res);

repo: SelfExSR


text:
function S = timestr(D,precision)
%timestr          String representation of time.   HH:MM:SS.SSSS
%
% TIMESTR(D) converts D, a serial date number (as returned by DATENUM)
% into a time string with the format HH:MM:SS.SSSS
%
% TIMESTR(D,precision) uses precision values to the right of the decimal

% Copyright 2003 The MathWorks, Inc

if nargin==1
    precision=4;
end;
totalwidth = precision + 3;         %2 to left of decimal, plus decimal
precision = num2str(precision);
totalwidth = num2str(totalwidth);

D = D(:);

% Obtain components of date number
[y,mo,d,h,min,s] = datevecmx(D,1.1);  mo(mo==0) = 1;

% Generate formatted string
% sw = floor(s);      %Whole
% sf = floor((s-sw)*1000);          %Fraction
% M = [h';min';sw';sf'];       %sprintf works columnwise
% fmt = '%02d:%02d:%02d.%04d';
M = [h';min';s'];       %sprintf works columnwise

% Figure out how long to make seconds format.
% Since we are building a string array, every element must be the same
% length
sw = floor(s);      %Whole
sf = floor((s-sw)*1000);          %Fraction

fmt = ['%02d:%02d:%' totalwidth '.' precision 'f'];
S = sprintf(fmt,M(:,1));
for ii=2:length(D)
    t= sprintf(fmt,M(:,ii));
    S = [S;t];
end;

% My formatting is a bit messed up.  I can't figure out how to add
% zeros when necessary for seconds - I end up with blanks.
% Replace blanks with 0
blnks = find(double(S)==32);        % 32 - ASCII for 0
S(blnks) = '0';

author: Sable

repo: mcbench-benchmarks


text:
function [ft3] = l2ft3(l)
% Convert volume from liters to cubic feet.
% Chad Greene 2012
ft3 = l*0.035314666721;


text:
% Copyright (C) 2009 International Business Machines
% All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% $Id: hs071_c.c 699 2006-04-05 21:05:18Z andreasw $
%
% Author:  Andreas Waechter               IBM    2009-04-02
%
% This file is part of the Ipopt tutorial.  It is a version with
% mistakes for the matlab implemention of the coding exercise problem
% (in AMPL formulation):
%
% param n := 4;
%
% var x {1..n} <= 0, >= -1.5, := -0.5;
%
% minimize obj:
%   sum{i in 1..n} (x[i]-1)^2;
%
% subject to constr {i in 2..n-1}:
%   (x[i]^2+1.5*x[i]-i/n)*cos(x[i+1]) - x[i-1] = 0;
%
% The constant term "i/n" in the constraint is supposed to be input data

function [x, info] = TutorialMatlab

  % Size of the problem
  n = 5;

  % Problem data
  a = ((1:n)/n)';

  % Starting point
  x0 = -0.5*ones(n,1);

  % Lower and upper bounds for the variables
  options.lb = -1.5*ones(n,1);
  options.ub = zeros(n,1);

  % Constraint bounds
  options.cl = zeros(n-2,1);
  options.cu = zeros(n-2,1);

  % Set the Ipopt options
  % options.ipopt.mu_strategy = 'adaptive';
  options.ipopt.derivative_test = 'second-order';
  options.ipopt.max_iter = 10;

  % Set the callback functions
  funcs.objective         = @eval_f;
  funcs.constraints       = @eval_g;
  funcs.gradient          = @eval_grad_f;
  funcs.jacobian          = @eval_jac_g;
  funcs.jacobianstructure = @eval_jac_g_struct;
  funcs.hessian           = @eval_hess;
  funcs.hessianstructure  = @eval_hess_struct;

  [x info] = ipopt(x0, funcs, options);

  % End of main function
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Evaluate value of objective function
  function f = eval_f(x)

    tmp = (x - 1).^2;
    f = sum(tmp);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Evaluate gradient of objective function
  function df = eval_grad_f(x)

    df = (x - 1);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Evaluate value of constraint bodies
  function g = eval_g(x)

    g = (x(2:n-1).^2 + 1.5*x(1:n-2) - a(2:n-1)).*cos(x(3:n)) - x(1:n-2);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Return constraint Jacobian strcture
  function A = eval_jac_g_struct

    % tri-diagonal structure
    Diags = ones(n, 3);
    A = spdiags(Diags, 0:2, n-2, n);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Evaluate constraint Jacobian
  function A = eval_jac_g(x)

    Diags = -1*ones(n-2,3);
    Diags(:,2) = (2*x(2:n-1)+1.5).*cos(x(3:n));
    Diags(:,3) = -(x(2:n-1).^2 + 1.5*x(2:n-1) - a(2:n-1)).*sin(x(3:n));
    A = spdiags(Diags, 0:2, n-2, n);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Return Hessian of Lagrangian function structure
  function H = eval_hess_struct

    Diags = ones(n, 2);
    H = spdiags(Diags, -1:0, n, n);

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Evaluate Hessian of Lagrangian function
  function H = eval_hess(x, sigma, lambda)

    Diags = zeros(n,2);

    % part from the objective function
    Diags(:,2) = sigma * 2;

    % (x_i , x_i) part
    Diags(2:n-1,2) = Diags(2:n-1,2) + lambda.*cos(x(3:n));

    % (x_{i+1}, x_{i+1}) part
    Diags(3:n,2) = Diags(3:n,2) ...
      - lambda.*(x(2:n-1).^2 + 1.5*x(1:n-2) - a(2:n-1)).*cos(x(3:n));

    % (x_i, x_{i+1}) part
    Diags(2:n-1,1) = -lambda.*sin(x(3:n)).*(2*x(2:n-1)+1.5);

    H = spdiags(Diags, -1:0, n, n);

  end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function result = sumOfNeighbors(input,edges,edgeOffsets,numNeighbors)

numNodes = length(input);
result = zeros(size(input));
for n = 1:numNodes
   for e = [edgeOffsets(n):edgeOffsets(n)+numNeighbors(n)-1]
      result(n) = result(n) + input(edges(e));
   end
end

% ras 03/2007:
% let's make this agree with sumOfNeighbors.mex: it outputs a
% transverse:
result = result';

return


text:
classdef MFEA < ALGORITHM
% <single> <real/integer/label/binary/permutation> <large/none> <multitask>
% Multifactorial evolutionary algorithm
% rmp --- 0.3 --- Random mating probability

%------------------------------- Reference --------------------------------
% A. Gupta, Y. Ong, and L. Feng, Multifactorial evolution: towards
% evolutionary multitasking, IEEE Transactions on Evolutionary Computation,
% 2016, 20(3): 343-357.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            rmp = Algorithm.ParameterSet(0.3);

            %% Initialize population
            Population    = Problem.Initialization();
            SubPopulation = Divide(Population,length(Problem.SubD));

            %% Optimization
            while Algorithm.NotTerminated([SubPopulation{:}])
                SubParent    = SelectHalf(SubPopulation);
                SubOffspring = Child(Problem,SubParent,rmp);
                for i = 1 : length(Problem.SubD)
                    SubPopulation{i}= [SubPopulation{i},SubOffspring{i}];
                end
                SubPopulation = SelectPop(SubPopulation,Problem);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function t = spm_convmtx(v,n,OPT)
% as for convmtx but with boundary conditions
% FORMAT t = spm_convmtx(C,N,OPT)
%
% OPT  - 'circular' boundary conditions
%      - 'square'   top and tail convolution matrix
%
%--------------------------------------------------------------------------
%   CONVMTX(C,N) returns the convolution matrix for vector C.
%   If C is a column vector and X is a column vector of length N,
%   then CONVMTX(C,N)*X is the same as CONV(C,X).
%   If R is a row vector and X is a row vector of length N,
%   then X*CONVMTX(R,N) is the same as CONV(R,X).
%   See also CONV.%
%   With the circular option the convolution matrix is reduced to N X N
%__________________________________________________________________________
% Copyright (C) 1988-2004 The MathWorks, Inc.

% L. Shure and T. Krauss
% $Id: spm_convmtx.m 6122 2014-07-25 13:48:47Z karl $

if nargin < 3;
    OPT = 'none';
end

% create Toeplitz matrix
%--------------------------------------------------------------------------
[mv,nv] = size(v);
v    = full(v(:));                              % make v a column vector
c    = [v; zeros(n-1,1)];
r    = zeros(n,1);
m    = length(c);
x    = [r(n:-1:2) ; c(:)];                      % build vector of user data
cidx = (0:m-1)';
ridx = n:-1:1;
t    = cidx(:,ones(n,1)) + ridx(ones(m,1),:);   % Toeplitz subscripts
t(:) = x(t);                                    % actual data

% transpose if necessary
%--------------------------------------------------------------------------
if mv < nv
    t = t.';
end

% apply optional boundary conditions
%--------------------------------------------------------------------------
switch OPT

    case('circular')
        m      = fix((size(t,1) - n)/2);
        j      = (1:m) + m;
        t(j,:) = t(j,:) + t(j + n,:);
        j      = (1:m) + n;
        t(j,:) = t(j,:) + t(j - n,:);
        j      = 1:n;
        t      = t(j + m,:);

    case('square')
        m      = fix((size(t,1) - n)/2);
        j      = 1:n;
        t      = t(j + m,:);

    otherwise
end

author: spm


text:
%%******************************************************************
%% dwd: distance weighted discrimination.
%%
%% [w,beta,residp,residn,totalviolation,X,y,Z] = ...
%%             dwd(Ap,An,penalty);
%%
%% You're given two matrices, A+ and A-, whose columns
%% give points in R^n you want to separate "nicely" to allow
%% classification of future points (e.g., to see if a patient
%% has breast cancer or not). The "classical" solution puts a
%% hyperplane between the points to maximize the minimum distance
%% from a point to a hyperplane. Steve is interested in smaller
%% instances, and is uncomfortable with the criterion, which is
%% too subject to noise (like l_\infty fitting). He wants to use
%% a criterion that depends on all the points. Initially he wanted
%% to minimize the sum of the inverse squares of the distances, but
%% he had no strong reason for the inverse square, so I suggested
%% just the inverse distance. This turns into a nice socp!
%% If there are mp points on the positive side and mn on the negative,
%% the dimensions are: 2*(mp+mn)+1 constraints,
%% (n+1) + (2) + (mp+mn)*3 variables in SOCs, with the dimensions as
%% given (one large, the rest small), and mp+mn nonneg. variables.
%% The latter are artificial variables to account for points that
%% are on the wrong side of the hyperplane: they can be moved a
%% distance d at a cost penalty*d to put them on the right side.
%%
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [w,beta,residp,residn,totalviolation,X,y,Z] = ...
             dwd(Ap,An,penalty);

   [np,mp] = size(Ap);
   [nn,mn] = size(An);
   n = np;
   nv = 1 + n + 2 + 3*(mp + mn) + mp + mn;
   nc = 1 + 2*(mp + mn);
%%
   blk = cell(2,2);
   blk{1,1} = 'q';
   blk{1,2} = [n+1, 2, 3*ones(1, mp + mn)];
   blk{2,1} = 'l';
   blk{2,2} = mp + mn;
%%
   Avec = cell(2,1);
   A = sparse(nc,nv-mp-mn);
   A(1:mp,2:n+3) = [Ap',zeros(mp,1),ones(mp,1)];
   A(mp+1:mp+mn,2:n+3) = -[An',zeros(mn,1),ones(mn,1)];
   A(1:mp,n+4:3:n+3+3*mp) = - speye(mp,mp);
   A(1:mp,n+6:3:n+5+3*mp) = - speye(mp,mp);
   A(mp+1:mp+mn,3*mp+n+4:3:3*mp+n+3+3*mn) = - speye(mn,mn);
   A(mp+1:mp+mn,3*mp+n+6:3:3*mp+n+5+3*mn) = - speye(mn,mn);
   A(mp+mn+1,1) = 1;
   A(mp+mn+2:mp+mn+1+mp,n+5:3:n+4+3*mp) = speye(mp,mp);
   A(mp+mn+1+mp+1:mp+mn+1+mp+mn,3*mp+n+5:3:3*mp+n+4+3*mn) = speye(mn,mn);
%%
   Avec{1,1} = A;
   Avec{2,1} = [speye(mp+mn,mp+mn);zeros(1+mp+mn,mp+mn)];
   b = [zeros(mp+mn,1);ones(1+mp+mn,1)];
%%
   C = cell(2,1);
   c = zeros(nv-mp-mn,1);
   c(n+4:3:n+3+3*mp) = ones(mp,1);
   c(n+6:3:n+5+3*mp) = -ones(mp,1);
   c(3*mp+n+4:3:3*mp+n+3+3*mn) = ones(mn,1);
   c(3*mp+n+6:3:3*mp+n+5+3*mn) = -ones(mn,1);
   C{1,1} = c;
   C{2,1} = penalty*ones(mp+mn,1);
%%
   [obj,X,y,Z] = sqlp(blk,Avec,C,b);
   X1 = X{1};
   omega = X1(1);
   w = X1(2:n+1);
   alpha = X1(n+2);
   beta  = X1(n+3);
   residp = Ap'*w + beta*ones(mp,1);
   residn = An'*w + beta*ones(mn,1);
   X2 = X{2};
   totalviolation = sum(X2);
%%******************************************************************

author: yu-jiang

repo: radpbook

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


text:
function s = saveobj(obj)
%SAVEOBJ  Save a DagNN to a vanilla MATLAB structure
%   S = OBJ.SAVEOBJ() saves the DagNN OBJ to a vanilla MATLAB
%   structure S. This is particularly convenient to preserve future
%   compatibility and to ship networks that are pure structures,
%   instead of embedding dependencies to code.
%
%   The object can be reconstructe by `obj = DagNN.loadobj(s)`.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

device = obj.device ;
obj.move('cpu') ;
s.vars = struct(...
  'name', {obj.vars.name}, ...
  'precious', {obj.vars.precious}) ;
s.params = struct(...
  'name', {obj.params.name}, ...
  'value', {obj.params.value}, ...
  'learningRate', {obj.params.learningRate}, ...
  'weightDecay', {obj.params.weightDecay}) ;
s.layers = struct(...
  'name', {obj.layers.name}, ...
  'type', {[]}, ...
  'inputs', {obj.layers.inputs}, ...
  'outputs', {obj.layers.outputs}, ...
  'params', {obj.layers.params}, ...
  'block', {[]}) ;
s.meta = obj.meta ;

for l = 1:numel(obj.layers)
  block = obj.layers(l).block ;
  slayer = block.save() ;
  s.layers(l).type = class(block) ;
  s.layers(l).block = slayer ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function infrasoundEventOut = auto_measure_amplitudes(infrasoundEvent, wevent, titleStr, predicted_traveltime_seconds, figureOutDirectory)
numEvents = numel(infrasoundEvent);
close all
min_traveltime = min(predicted_traveltime_seconds);
relative_traveltimes = predicted_traveltime_seconds - min_traveltime;
infrasoundEventOut = [];
for eventNumber=1:numEvents
    thisEvent = infrasoundEvent(eventNumber);
    thisW=wevent{eventNumber};

    thisEvent = auto_measure_minmax3(thisW, thisEvent, relative_traveltimes)
    ah=get(gcf,'Children');
    title(ah(7),sprintf('Event %d',eventNumber));

    feval('print', '-dpng', sprintf('%s/%s_%03d',figureOutDirectory, titleStr, eventNumber) );
    close
    infrasoundEventOut = [infrasoundEventOut thisEvent];
end

author: geoscience-community-codes

repo: GISMO


text:
function [rstat]=AddObsErr_v000(finname, nrow, foutname, act_row, errdef, rstat)

%randomize
if (rstat(1)~=0)
    randn('state',rstat)
else
    rstat=randn('state');
end

finp=fopen(finname,'rb');
inp=fread(finp,[nrow,inf],'double');
fclose(finp);

for in=1:size(act_row,1)
    fout=fopen(foutname(1,:),'wb');
    act_in=act_row(in,:);
    nd=size(act_in,2);
    if (~isempty(errdef))
        sR=errdef(in).sR;
        for (in1=1:size(inp,2))
            out=inp(act_in,in1)+sR*randn(nd,1);
            fwrite(fout,[inp(1,in1);out],'double');
        end
    else
        out=inp([1 act_in],:);
        fwrite(fout,out,'double');
    end
    fclose(fout);
end

repo: ekfmonoslam


text:
function [rH, lH, isConnected, iStruct, iRightScout, iLeftScout] = tess_hemisplit(sSurf)
% HEMISPLIT: Find the right and left hemisphere indexes
%
% INPUT:
%     - sSurf: Brainstorm surface structure with fields: Vertices, VertConn, Atlas
%
% OUTPUT:
%     - rH          : right hemisphere indexes
%     - lH          : left hemisphere indexes
%     - isConnected : 1 if the two hemispheres are connected, else 0
%     - iStruct     : Index of the atlas "Structures"
%     - iLeftScout  : Index of the scout "lh" in the atlas "Structures"
%     - iRightScout : Index of the scout "rh" in the atlas "Structures"

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Guillaume Dumas, Francois Tadel, 2009-2013

% Initialize returned values
rH          = [];
lH          = [];
isConnected = 0;
iLeftScout  = [];
iRightScout = [];

%% ===== USE STRUCTURES ATLAS =====
% If there is a Structures atlas with left and right hemispheres: split in two spheres
iStruct = find(strcmpi({sSurf.Atlas.Name}, 'Structures'));
if ~isempty(iStruct)
    % === Method 1: Find the left and right hemispheres ===
    % Find left and right hemispheres
    iRightScout = [find(strcmpi({sSurf.Atlas(iStruct).Scouts.Label}, 'rh')), find(strcmpi({sSurf.Atlas(iStruct).Scouts.Label}, 'Cortex') & strcmpi({sSurf.Atlas(iStruct).Scouts.Region}, 'RU'))];
    iLeftScout  = [find(strcmpi({sSurf.Atlas(iStruct).Scouts.Label}, 'lh')), find(strcmpi({sSurf.Atlas(iStruct).Scouts.Label}, 'Cortex') & strcmpi({sSurf.Atlas(iStruct).Scouts.Region}, 'LU'))];
    % If both hemispheres are described here: get the indices
    if ~isempty(iRightScout) && ~isempty(iLeftScout)
        rH = sSurf.Atlas(iStruct).Scouts(iRightScout).Vertices(:)';
        lH = sSurf.Atlas(iStruct).Scouts(iLeftScout).Vertices(:)';
        return;
    end
    % === Method 2: Sum all the L/R scouts ===
    % Find all the L/R scouts
    scoutHemi = cellfun(@(c)c(1), {sSurf.Atlas(iStruct).Scouts.Region}, 'UniformOutput', 0);
    iRightScouts = find(strcmpi(scoutHemi, 'R'));
    iLeftScouts  = find(strcmpi(scoutHemi, 'L'));
    % If both hemispheres are described here: get the indices
    if ~isempty(iRightScouts) && ~isempty(iLeftScouts)
        rH = unique([sSurf.Atlas(iStruct).Scouts(iRightScouts).Vertices]);
        lH = unique([sSurf.Atlas(iStruct).Scouts(iLeftScouts).Vertices]);
        % Make sure these are row vectors
        rH = rH(:)';
        lH = lH(:)';
        return;
    end
end

%% ===== GROW FROM EACH SIDE =====
% Perform detection only in dimension 2
dim = 2;
% Consider that all the 30% extremes on the left-right axis belong only to one hemisphere
START_PERCENT = 0.3;
% Right
if ~isConnected
    % Get maximal y value
    [yMin, rH] = min(sSurf.Vertices(:,dim));
    % Get all the vertices that are at more than START_PERCENT of the maximum
    iNewR = find(sSurf.Vertices(:,dim) < START_PERCENT * yMin)';
    iNewR = setdiff(iNewR, rH);
    % If not displaying the whole brain
    if ~isempty(rH)
        % Grow region until getting all the hemisphere
        while ~isConnected && ~isempty(iNewR)
            rH = union(rH, iNewR);
            iNewR = tess_scout_swell(rH, sSurf.VertConn);
            % If it's including more than 80% of the surface: connected...
            if (length(rH) >= .8 * length(sSurf.Vertices))
                isConnected = 1;
            end
        end
    end
end
% Left
if ~isConnected
    % Get minimal y value
    [yMax, lH] = max(sSurf.Vertices(:,dim));
    % Get all the vertices that are at < than START_PERCENT of the minimum
    iNewL = find(sSurf.Vertices(:,dim) > START_PERCENT * yMax)';
    iNewL = setdiff(iNewL, lH);
    % If not displaying the whole brain
    if ~isempty(lH)
        % Grow region until getting all the hemisphere
        while ~isConnected && ~isempty(iNewL)
            lH = union(lH, iNewL);
            iNewL = tess_scout_swell(lH, sSurf.VertConn);
            % If it's including more than 80% of the surface: connected...
            if (length(lH) >= .8 * length(sSurf.Vertices))
                isConnected = 1;
            end
        end
    end
end

%% ===== SPLIT AT Y=0 =====
if isConnected
    lH = find(sSurf.Vertices(:,dim) < 0);
    rH = find(sSurf.Vertices(:,dim) >= 0);
end

author: brainstorm-tools


text:
function [kip] = kN2kip(kN)
% Convert force from kilonewtons to kip.
% Chad A. Greene 2012
kip = kN* 0.22480894387;


text:
classdef MultiSubspaceSignalGenerator < handle
    %MultiSubspaceSignalGenerator creates sparse signals
    % Usage:
    %   - Create generator with ambient dimension N and sparsity level K
    %   - Design the support for subspaces
    %   - Specify number of signals for each subspace
    %   - Generate signals for each subspace
    properties(SetAccess=private)
        % The ambient signal dimension
        N
        % The sparsity level
        K
        % The resultant vectors
        X
        % The Sparse support for each subspace
        Supports
        % The number of signals in each Subspace
        NumSignalsPerSubspace = []
        % Labels identify the subspace to which a signal belongs
        Labels = []
    end

    methods
        function self = MultiSubspaceSignalGenerator(N, K)
            self.N = N;
            self.K = K;
        end

        function createDisjointSupports(self, numSubspaces)
            n = self.N;
            k  = self.K;
            c = numSubspaces;
            if c*k > n
                error('Total cardinality of subspaces cannot exceed ambient dimensions');
            end
            % Select c*k indices randomly from 1:N
            q = randperm(n, c*k);
            % Arrange the support for each subspace.
            self.Supports = sort(reshape(q, k, c));
            % By default we will create one signal per subspace
            self.NumSignalsPerSubspace = ones(numSubspaces, 1);
        end

        function createOverlappingSupports(self, numSubspaces, overlap)
            % Overlap specifies number of indices which overlap between
            % two consecutive supports.
            n = self.N;
            k  = self.K;
            c = numSubspaces;
            l = c*(k-overlap)+overlap;
            if  l > n
                error('Total cardinality of subspaces cannot exceed ambient dimensions, n=%d, l=%d', n, l);
            end
            % Select l indices randomly from 1:N
            q = randperm(n, l);
            % Arrange the support for each subspace.
            qs = zeros(k,c);
            start = 0;
            for i=1:c
                qs(:, i) = q(start + (1:k));
                start = start + k - overlap;
            end
            self.Supports = sort(qs);
            % By default we will create one signal per subspace
            self.NumSignalsPerSubspace = ones(numSubspaces, 1);
        end

        function setNumSignalsPerSubspace(self, numSignalsPerSubspace)
            c = length(self.NumSignalsPerSubspace);
            data = numSignalsPerSubspace;
            if isscalar(data)
                self.NumSignalsPerSubspace = data *ones(c, 1);
                return;
            end
            % Verify that numSignalsPerSubspace has appropriate size
            if length(data) ~= c
                error('Number of subspaces mismatch');
            end
            % Make sure its a column vector
            if isrow(data)
                data = data';
            end
            % save it
            self.NumSignalsPerSubspace = data;
        end

        function result = uniform(self, a,b)
            % Generates sparse signals from uniform distribution
            if nargin < 3
                a = 0; b =1;
            elseif nargin < 4
                b = -a;
            end
            if b <= a
                error('b must be larger than a');
            end
            % sparsity level
            k = self.K;
            % number of signals per subspace
            signalCounts = self.NumSignalsPerSubspace;
            % total number of signals
            totalSignals = sum(signalCounts);
            % Ambient dimensions
            n = self.N;
            % supports
            qs = self.Supports;
            % Number of subspaces
            c = length(qs);
            % Create space for signals
            self.X = zeros(n, totalSignals);
            self.Labels = zeros(totalSignals, 1);
            % Iterate over subspaces
            start = 0;
            for i=1:c
                % support for this subspace
                q = qs(:, i);
                % number of signals
                s = signalCounts(i);
                % Generate signals
                xx = a + (b-a).*rand(k,s);
                % Assign into combined set
                self.X(q, start + (1:s)) = xx;
                self.Labels(start + (1:s)) = i;
                % update start for next set
                start = start + s;
            end
            % Store the results
            result = self.X;
        end

        function result = biUniform(self, a, b, zeroProbability)
            % Generates sparse vectors where each non-zero values
            % is picked up uniformly from the ranges
            % [-b, -a] and [a, b]
            if nargin < 3
                a = 1; b =2;
            elseif nargin < 4
                b = 2*a;
            end
            if a < 0 || b < 0
                error('a and b both must be +ve');
            end
            if b <= a
                error('b must be larger than a');
            end
            if ~exist('zeroProbability', 'var')
                zeroProbability = 0;
            end
            if zeroProbability > 1 || zeroProbability < 0
                error('zeroProbability must be in [0,1]');
            end
            % sparsity level
            k = self.K;
            % number of signals per subspace
            signalCounts = self.NumSignalsPerSubspace;
            % total number of signals
            totalSignals = sum(signalCounts);
            % Ambient dimensions
            n = self.N;
            % supports
            qs = self.Supports;
            % Number of subspaces
            c = size(qs,2);
            % Create space for signals
            self.X = zeros(n, totalSignals);
            self.Labels = zeros(totalSignals, 1);
            % Iterate over subspaces
            start = 0;
            for i=1:c
                % support for this subspace
                q = qs(:, i);
                % number of signals
                s = signalCounts(i);
                % Generate signals
                % unsigned result
                xx = a + (b-a).*rand(k,s);
                % sign
                sgn = sign(randn(k,s));
                if zeroProbability
                    % We need to introduce some zeros in between
                    sgn = sgn .* binornd(1, 1 - zeroProbability, k, s);
                end
                % Assign into combined set
                self.X(q, start + (1:s)) = sgn .* xx;
                self.Labels(start + (1:s)) = i;
                % update start for next set
                start = start + s;
            end
            % Store the results
            result = self.X;
        end
    end

end

author: indigits

repo: sparse-plex