text:
% GP_COV - This is a general documentation file for using covariance
%          functions.
%
% The covariance functions are used similarly. The covariance function is
% initialized by giving some function specific inputs, for instance,
%
%   COVFUNC = GP_COV_EXAMPLE(P1, P2, ...)
%
% where the number of inputs and their meaning depends on the covariance
% function. The input arguments can be, for instance, a (squared) distance
% matrix between the inputs, the dimensionality of the input space or
% another covariance function. Thus, be careful especially when changing
% covariance functions that you also change the inputs as needed! For
% instance, some covariance functions may use a distance matrix while
% some others a squared distance matrix.
%
% The returned covariance function COVFUNC can be used to obtain the
% covariance matrix by giving a parameter vector THETA:
%
%   K = COVFUNC(THETA)
%   [K, DK] = COVFUNC(THETA)
%
% where K is the covariance matrix and DK is a cell array of derivatives of
% K with respect to the parameters THETA. The length of THETA and the
% meaning of its elements depends on the covariance function. The parameters
% can be, for instance, characteristic length scale, magnitude or
% period. Again, be sure that you know which element of THETA corresponds to
% which parameter of the covariance function. Especially, if you modify your
% covariance function construction, also modify THETA appropriately!
%
% You can form complex covariance functions by using some functions to chain
% simple covariance functions. For instance, GP_COV_SCALE can be used to add
% a scaling parameter, GP_COV_SUM to sum two covariance functions and
% GP_COV_PRODUCT to multiply two covariance functions. Make sure how the
% additional parameters are added to the vector THETA and how the parameter
% vectors of multiple covariance functions are combined.
%
% For most covariance functions, you can fix some parameters:
%
%   COVFUNC = GP_COV_EXAMPLE(..., 'PARAMETER', VALUE)
%
% Then the parameter is not read from the vector THETA anymore. For
% instance, if you construct a covariance function as
%
%   COVFUNC = GP_COV_PERIODIC(SQRT(SQ_DIST(X1,X2)),'WAVELENGTH',10)
%
% then THETA should contain only one element, the smoothness parameter.
%
% You can obtain some dimensionality information by calling COVFUNC
% without any inputs:
%
%   [N_THETA, N1, N2] = COVFUNC()
%
%   N_THETA : Number of non-fixed parameters
%   N1      : Number of rows in the covariance matrix
%   N2      : Number of columns in the covariance matrix
%
% In order to obtain the covariance matrix of a covariance function with no
% parameters, you should call COVFUNC([]), not COVFUNC(). Note the
% difference!
%
% At least the following covariance functions are implemented:
%
%   GP_COV_SE       : Isotropic squared exponential
%   GP_COV_RQ       : Isotropic rational quadratic
%   GP_COV_PP       : Piecewise polynomial
%   GP_COV_DELTA    : Delta function, isotropic noise
%   GP_COV_SCALE    : Scales a covariance function
%   GP_COV_SUM      : Sums two covariance functions
%   GP_COV_PRODUCT  : Multiplies two covariance functions
%   GP_COV_JITTER   : Adds a small constant to the diagonal for numerical
%                     stability
%   GP_COV_TOEPLITZ : Creates a covariance matrix with Toeplitz structure
%   GP_COV_WRAP     : Makes a covariance function from a covariance matrix
%
% You can write your own covariance functions as long as they fulfill the
% specification described here. In some cases, it might be a good idea to
% write some wrapper functions for the covariance functions to make the
% usage more straightforward. For instance,
%
%   COVFUNC_SE = @(THETA,X1,X2) FEVAL(GP_COV_SE(SQ_DIST(X1,X2)),THETA)
%
% would be used as
%
%   [K, DK] = COVFUNC_SE(THETA, X1, X2)
%
% which some users may prefer. However, this kind of interface may reduce
% efficiency, as one may need to compute same things several times (squared
% distance matrix in the above example).
%
% Note: Despite the similar naming convention, GP_COV_PSEUDO is
% fundamentally different kind of covariance function and can not be used
% similarly to other covariance functions. See GP_COV_PSEUDO for more
% information.
%
% See also GP_LEARN, GP_LEARN_PSEUDO, GP_PREDICT, GP_PREDICT_PSEUDO.

% Last modified 2010-01-27
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function gp_cov()

help gp_cov

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
classdef vector3d < dynOption
%
% The class vector3d describes three dimensional vectors, given by
% their coordinates x, y, z and allows to calculate with them as
% comfortable as with real numbers.
%
% Syntax
%   v = vector3d(x,y,z)
%   v = vector3d(x,y,z,'antipodal')
%   v = vector3d.byPolar(theta,rho)
%
% Input
%  x,y,z - cart. coordinates
%
% Output
%  v - @vector3d
%
% Flags
%  antipodal - <VectorsAxes.html consider vector as an axis>
%
% Class Properties
%  x, y, z      - cart. coordinates
%  isNormalized - whether the vector is a direction
%  antipodal    - <VectorsAxes.html whether the vector is an axis>
%
% Dependent Class Properties
%  theta      - polar angle in radiant
%  rho        - azimuthal angle in radiant
%  resolution - mean distance between the points on the sphere
%  xyz        - cart. coordinates as matrix
%
% Derived Classes
%  @Miller - crystal directions
%  @S2Grid - sphercial grid
%
% See also
% VectorDefinition VectorsOperations VectorsAxes VectorsImport VectorsExport

  properties
    x = []; % x coordinate
    y = []; % y coordinate
    z = []; % z coordinate
    antipodal = false;
    isNormalized = false;
  end

  properties (Dependent = true)
    theta   % polar angle
    rho     % azimuth angle
    resolution % mean distance between the points on the sphere
    xyz
  end

  methods

    function v = vector3d(varargin)
      % constructor of the class vector3d

      if nargin >=3 && isnumeric(varargin{1})

        v.x = varargin{1};
        v.y = varargin{2};
        v.z = varargin{3};

      elseif nargin == 0
      elseif nargin <= 2
        if strcmp(class(varargin{1}),'vector3d') %#ok<STISA>

          v = varargin{1};

        elseif isa(varargin{1},'vector3d') % copy-constructor

          v.x = varargin{1}.x;
          v.y = varargin{1}.y;
          v.z = varargin{1}.z;
          v.antipodal = varargin{1}.antipodal;
          v.isNormalized = varargin{1}.isNormalized;
          v.opt = varargin{1}.opt;
          return

        elseif isa(varargin{1},'double')
          xyz = varargin{1};
          if all(size(xyz) == [1,3])
            xyz = xyz.';
          end
          v.x = xyz(1,:);
          v.y = xyz(2,:);
          v.z = xyz(3,:);
        else
          error('wrong type of argument');
        end
      elseif ischar(varargin{1})

        if strcmp(varargin{1},'polar')

          sy = sin(varargin{2});
          v.x = sy .* cos(varargin{3});
          v.y = sy .* sin(varargin{3});
          v.z = cos(varargin{2});

        else

          theta = get_option(varargin,{'theta','azimuth'});
          rho = get_option(varargin,{'rho','polar'});

          v.x = sin(theta).*cos(rho);
          v.y = sin(theta).*sin(rho);
          v.z = cos(theta);

        end

      end

      % ----------- check for equal size ------------------------
      if numel(v.x) ~= numel(v.y) || (numel(v.x) ~= numel(v.z))

        % find non singular size
        if numel(v.x) > 1
          s = size(v.x);
        elseif numel(v.y) > 1
          s = size(v.y);
        else
          s = size(v.z);
        end

        % try to correct
        if numel(v.x) == 1, v.x = repmat(v.x,s);end
        if numel(v.y) == 1, v.y = repmat(v.y,s);end
        if numel(v.z) == 1, v.z = repmat(v.z,s);end

        % check again
        if numel(v.x) ~= numel(v.y) || (numel(v.x) ~= numel(v.z))
          error('MTEX:Vector3d','Coordinates have different size.');
        end
      end

      % ------------------ options ------------------------------

      if nargin > 1

        % antipodal
        v.antipodal = check_option(varargin,'antipodal');

        % resolution
        if check_option(varargin,'resolution')
          v = v.setOption('resolution',get_option(varargin,'resolution'));
        end

        % normalize
       if check_option(varargin,'normalize'), v = normalize(v); end

      end
    end

    function n = numArgumentsFromSubscript(varargin)
      n = 0;
    end

    function rho = get.rho(v)
      try
        rho = v.opt.rho;
      catch
        rho = atan2(v.y,v.x);
      end
    end

    function theta = get.theta(v)
      try
        theta = v.opt.theta;
      catch
        theta = acos(v.z./v.norm);
      end
    end

    function xyz = get.xyz(v)

      xyz = [v.x(:),v.y(:),v.z(:)];

    end

    function res = get.resolution(v)

      if v.isOption('resolution')
        res = v.getOption('resolution');
      elseif length(v) <= 4
        res = 2*pi;
      elseif length(v) > 50000
        res = sqrt(40000 / length(v) / (1 + v.antipodal)) * degree;
      else
        try
          a = calcVoronoiArea(v);
          res = sqrt(median(a));
          assert(res>0);
        catch
            res = 2*pi;
        end
      end
    end

    function v = set.resolution(v,res)

      v = v.setOption('resolution',res);

    end

    function b = isnan(v)
      b = isnan(v.x) | isnan(v.y) | isnan(v.z);
    end

    function b = isinf(v)
      b = isinf(v.x) | isinf(v.y) | isinf(v.z);
    end

    function b = isfinite(v)
      b = ~(isinf(v) | isnan(v));
    end

    function v = real(v)
      v = vector3d(real(v.x),real(v.y),real(v.z));
    end

    function v = imag(v)
      v = vector3d(imag(v.x),imag(v.y),imag(v.z));
    end

  end

  methods (Static = true)

    v = nan(varargin)
    v = ones(varargin)
    v = zeros(varargin)
    v = rand(varargin)
    v = byPolar(polarAngle,azimuthAngle,varargin)
    [v,interface,options] = load(fname,varargin)

    function v = X(varargin)
      % the vector (1,0,0)
      %
      % Syntax
      %   x = vector3d.X % returns a single vector (1,0,0)
      %   x = vector3d.X(3,1) % returns 3 vectors (1,0,0)

      x = ones(varargin{:});
      v = vector3d(x,0,0);
    end

    function v = Y(varargin)
      x = ones(varargin{:});
      v = vector3d(0,x,0);
    end

    function v = Z(varargin)
      x = ones(varargin{:});
      v = vector3d(0,0,x);
    end

  end
end

author: mtex-toolbox

repo: mtex

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


text:
function  a = i4vec_indicator1 ( n )

%*****************************************************************************80
%
%% I4VEC_INDICATOR1 sets an I4VEC to the indicator vector (1,2,3,...).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Output, integer A(N), the vector.
%
  a = ( 1 : n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mlstripcommentsfile(infile, outfile)
%MLSTRIPCOMMENTSFILE Strip comments from a file with MATLAB code.
%
%   MLSTRIPCOMMENTSFILE(INFILE, OUTFILE) reads the m-file INFILE, line by line,
%   strips all MATLAB comments , and writes the result to the file OUTFILE.
%
%   Note that you should never specify the same value for INFILE and OUTFILE.
%   Since the output file is opened for writing before any data is read from
%   the input file, the file will be empty and the data will be lost.
%
%   See also MLSTRIPCOMMENTSFID, MLSTRIPCOMMENTSSTR.

%   Author:      Peter J. Acklam
%   Time-stamp:  2004-03-18 08:23:19 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % Check number of input arguments.
   error(nargchk(1, 2, nargin));

   % See if the input file exists.
   if ~exist(infile, 'file')
      error([infile ': no such file or directory.']);
   end

   % Open the input file for reading.
   ifid = fopen(infile, 'rt');
   if ifid < 0
      error([infile ': can''t open file for reading.']);
   end

   % Open the output file for writing.
   ofid = fopen(outfile, 'wt');
   if ofid < 0
      closefid(ifid, infile);
      error([outfile ': can''t open file for writing.']);
   end

   % Strip MATLAB comments.
   status = mlstripcommentsfid(ifid, ofid);

   % See if an error occurred while writing to file.
   if status < 0
      closefid(ifid, infile);
      closefid(ofid, outfile);
      error([outfile ': error while writing to file.']);
   end

   % Close input and output file handles.
   closefid(ifid, infile);
   closefid(ofid, outfile);

function closefid(fid, filename)
%CLOSEFID Close a file identifier.
%
%   CLOSEFID(IFID, FILENAME) closes the specified file identifier.

   status = fclose(fid);
   if status < 0
      error([filename ': can''t close file.']);
   end

author: Sable

repo: mcbench-benchmarks


text:
function variance = arcsin_variance ( a )

%*****************************************************************************80
%
%% ARCSIN_VARIANCE returns the variance of the Arcsin PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the CDF.
%    A must be positive.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = a * a / 2.0;

  return
end


text:
function [A,Q] = lti_disc(F,L,Q,dt)
% LTI_DISC Equivalent discrete-time solution of an LTI SDE
%
% Syntax:
%   [A,Q] = lti_disc(F,L,Qc,dt)
%
% In:
%   F  - NxN Feedback matrix
%   L  - NxL Noise effect matrix        (optional, default identity)
%   Qc - LxL Diagonal Spectral Density  (optional, default zeros)
%   dt - Time step                      (optional, default 1)
%
% Out:
%   A - Transition matrix
%   Q - Discrete process covariance matrix
%
% Description:
%   Equivalent discrete-time solution of an LTI SDE of form
%
%     dx/dt = F x + L w,
%
%   where w(t) is a white-noise process with spectral density Qc.
%   Results in the following model:
%
%     x[k] = A x[k-1] + q, q ~ N(0,Q).
%
%   Can be used for integrating the model exactly over time steps,
%   which are multiples of dt.
%
% Copyright:
%   2019 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

  % Check number of arguments and defaults
  if nargin < 1
    error('Too few arguments');
  end
  if nargin < 2
    L = [];
  end
  if nargin < 3
    Q = [];
  end
  if nargin < 4
    dt = [];
  end
  if isempty(L)
    L = eye(size(F,1));
  end
  if isempty(Q)
    Q = zeros(size(F,1),size(F,1));
  end
  if isempty(dt)
    dt = 1;
  end

  % Closed-form integration of the transition matrix
  A = expm(F*dt);

  % Closed-form integration of the covariance matrix
  n   = size(F,1);
  Phi = [F L*Q*L'; zeros(n,n) -F'];
  AB  = expm(Phi*dt)*[zeros(n,n);eye(n)];
  Q   = AB(1:n,:)*A'; % A' = inv(AB((n+1):(2*n),:));

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
function varargout = readMesh_off(fileName)
%READMESH_OFF Read mesh data stored in OFF format.
%
%   [VERTICES, FACES] = readMesh_off(FILENAME)
%   Read the data stored in file FILENAME and return the vertex and face
%   arrays as NV-by-3 array and NF-by-N array respectively, where NV is the
%   number of vertices and NF is the number of faces.
%
%   MESH = readMesh_off(FILENAME)
%   Read the data stored in file FILENAME and return the mesh into a struct
%   with fields 'vertices' and 'faces'.
%
%   Example
%     [v, f] = readMesh_off('mushroom.off');
%     figure; drawMesh(v, f, 'faceColor', [0 1 0], 'edgeColor', 'none')
%     view([5 80]); light; lighting gouraud
%
%   See also
%     meshes3d, readMesh, writeMesh_off, readMesh_obj, drawMesh
%

% ------
% 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

%% Read header

% open file
f = fopen(fileName, 'r');
if f == -1
    error('matGeom:readMesh_off:FileNotFound', ...
        ['Could not open input file: ' fileName]);
end

% check format
line = fgetl(f);   % -1 if eof
if ~strcmp(line(1:3), 'OFF')
    error('matGeom:readMesh_off:FileFormatError', ...
        'Not a valid OFF file');
end

% number of faces and vertices
line = fgetl(f);
vals = sscanf(line, '%d %d');
nVertices = vals(1);
nFaces = vals(2);

%% Read vertex data
[vertices, count] = fscanf(f, '%f ', [3 nVertices]);
if count ~= nVertices * 3
    error('matGeom:readMesh_off:FileFormatError', ...
        ['Could not read all the ' num2str(nVertices) ' vertices']);
end
vertices = vertices';

%% Read Face data
% First try to read faces as an homogeneous array. It if fails, start from
% face offset and parse each face individually. In the latter case, faces
% can have different number of vertices.

% keep position of face info within file
faceOffset = ftell(f);

% read first face to assess number of vertices per face
line = fgetl(f);
if line == -1
    error('matGeom:readMesh_off:FileFormatError', ...
        'Unexpected end of file');
end
tokens = split(line);
face1 = str2double(tokens(2:end))' + 1;
nv = length(face1);

try
    % attempt to read the remaining faces assuming they all have the same
    % number of vertices
    pattern = ['%d' repmat(' %d', 1, nv) '\n'];
    [faces, count] = fscanf(f, pattern, [(nv+1) (nFaces-1)]);
    if count ~= (nFaces-1) * (nv+1)
        error('matGeom:readMesh_off:FileFormatError', ...
            'Could not read all the %d faces', nFaces);
    end

    % transpose, remove first column, use 1-indexing, and concatenate with
    % first face
    faces = [face1 ; faces(2:end,:)'+1];

catch
    % if attempt failed, switch to slower face-by-face parsing
    disp('readMesh_off: Inhomogeneous number of vertices per face, switching to face-per-face parsing');

    fseek(f, faceOffset, 'bof');

    % allocate cell array
    faces = cell(1, nFaces);

    % iterate over faces
    for iFace = 1:nFaces
        % read next line
        line = fgetl(f);
        if line == -1
            error('matGeom:readMesh_off:FileFormatError', ...
                'Unexpected end of file');
        end

        % parse vertex indices for current face
        tokens = split(line);
        faces{iFace} = str2double(tokens(2:end))' + 1;
    end
end

%% Post-processing

% close the file
fclose(f);

% format output arguments
if nargout < 2
    mesh.vertices = vertices;
    mesh.faces = faces;
    varargout = {mesh};
else
    varargout = {vertices, faces};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = sin(varargin)
%SIN   Sine of a CHEBFUN2.

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

[varargout{1:nargout}] = sin@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function f = squish(m,num)

% function f = squish(m,num)
%
% <m> is a matrix
% <num> is the positive number of initial dimensions to squish together
%
% return <m> squished.
%
% example:
% isequal(squish([1 2; 3 4],2),[1 3 2 4]')

% get the size of m
msize = [size(m) ones(1,num-ndims(m))];  % add ones to end if necessary

% calculate the new dimensions
newdim = [prod(msize(1:num)) msize(num+1:end)];

% do the reshape
f = reshape(m,[newdim 1]);  % tack on a 1 to handle the special case of squishing everything together

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function [filt] = highpassfilter(dat,Fs,Fhp,N,type,dir)

% HIGHPASSFILTER removes low frequency components from EEG/MEG data
%
% Use as
%   [filt] = highpassfilter(dat, Fsample, Fhp, N, type, dir)
% where
%   dat        data matrix (Nchans X Ntime)
%   Fsample    sampling frequency in Hz
%   Fhp        filter frequency
%   N          optional filter order, default is 6 (but) or 25 (fir)
%   type       optional filter type, can be
%                'but' Butterworth IIR filter (default)
%                'fir' FIR filter using MATLAB fir1 function
%   dir        optional filter direction, can be
%                'onepass'         forward filter only
%                'onepass-reverse' reverse filter only, i.e. backward in time
%                'twopass'         zero-phase forward and reverse filter (default)
%
% Note that a one- or two-pass filter has consequences for the
% strength of the filter, i.e. a two-pass filter with the same filter
% order will attenuate the signal twice as strong.
%
% See also LOWPASSFILTER, BANDPASSFILTER

% Copyright (c) 2003, Robert Oostenveld
%
% 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$

% set the default filter order later
if nargin<4
    N = [];
end

% set the default filter type
if nargin<5
  type = 'but';
end

% set the default filter direction
if nargin<6
  dir = 'twopass';
end

% Nyquist frequency
Fn = Fs/2;

% compute filter coefficients
switch type
  case 'but'
    if isempty(N)
      N = 6;
    end
    [B, A] = butter(N, max(Fhp)/Fn, 'high');
  case 'fir'
    if isempty(N)
      N = 25;
    end
    B = fir1(N, max(Fhp)/Fn, 'high');
    A = 1;
end

% apply filter to the data
switch dir
  case 'onepass'
    filt = filter(B, A, dat')';
  case 'onepass-reverse'
    dat  = fliplr(dat);
    filt = filter(B, A, dat')';
    filt = fliplr(filt);
  case 'twopass'
    filt = filtfilt(B, A, dat')';
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function TC = rotate(TC,varargin)

%Rotate component orientations
% TC = ROTATE(TC,BEARING) Rotates the horizontal components of a threecomp
% object so that the first horizontal component (N, R or 1) points in the
% direction given by BEARING. If the rotation is to north/east or
% radial/transverse then the horizontral channel names will be renamed to
% ??N/??E or ??R/??T accordingly. If the rotation is to an arbitrary
% direction, the channels will be renamed to ??2/??1.
%
% TC = ROTATE(TC) is shorthand for ROTATE(TC,-BACKAZIMUTH). This common
% usage rotates horizontal traces to the radial and transverse directions,
% where radial is inline (and away) from the source. The transverse
% direction points +90 from the radial. For example,for a wave arriving
% from due south, there is no rotation and N+ becomes R+, while E+ becomes
% T+. This useage requires that the backAzimuth property be set.

% TODO: DISALLOW NO BACKAZIMUTH AND NO ARGUMENT

% CHECK ARGUMENTS AND GET FINAL BEARING
if length(varargin)>=1  % bearing provided by user
    bearing = varargin{1};
    if ~isa(bearing,'double')
       error('Threecomp:rotate:badArgumentType','Rotation angle must be a number');
    end
    if bearing<-360 || bearing>720
       error('Threecomp:rotate:badArgumentValue','Rotation angle must be between 0 and 360');
    end
else
    bearing = NaN * zeros(size(TC));
    for n = 1:numel(TC)
        if ~isempty(TC(n).backAzimuth)
            bearing(n) = TC(n).backAzimuth+180;
        end
        bearing = reshape(bearing,size(TC));
        bearing = mod(bearing,360);
    end
end

% COMPUTE ROTATION ANGLES
orientation = get(TC,'ORIENTATION');
for n = 1:length(TC)
   if isempty(orientation)
       error('ThreeComp:rotate:noOrientations','Component orientations must be provided');
   end
end
rotAngle = mod(bearing-orientation(:,3),360);

% APPLY ROTATION
for n = 1:length(TC)
    data = double(TC(n).traces(2:3))';
    rotMatrix = [cosd(rotAngle(n)) sind(rotAngle(n)) ; -sind(rotAngle(n)) cosd(rotAngle(n))];
    Dout = rotMatrix * data;
    TC(n).traces(2) = set(TC(n).traces(2),'DATA',Dout(1,:));
    TC(n).traces(3) = set(TC(n).traces(3),'DATA',Dout(2,:));
    newOrientation = TC(n).orientation;
    newOrientation([3 5]) = mod(newOrientation([3 5])+rotAngle(n),360);
    TC(n).orientation = newOrientation;
end

% ADJUST CHANNEL NAMES
for n = 1:length(TC)
    channels = get(TC(n),'CHANNEL');
    chan2 = channels{1}(2);
    chan2 = chan2{1};
    chan3 = channels{1}(3);
    chan3 = chan3{1};
    if abs(TC(n).orientation(3)-0)<0.5
        chan2(3) = 'N';
        chan3(3) = 'E';
    elseif abs(mod((TC(n).orientation(3)+180),360)-TC(n).backAzimuth)<0.5
        chan2(3) = 'R';
        chan3(3) = 'T';
    else
        chan2(3) = '2';
        chan3(3) = '1';
    end
    TC(n).traces(2) = set(TC(n).traces(2),'CHANNEL',chan2);
    TC(n).traces(3) = set(TC(n).traces(3),'CHANNEL',chan3);
end

author: geoscience-community-codes

repo: GISMO


text:
function [F,h,x] = sedumi2yalmip(At,b,c,K)

nvars = length(b);
x = sdpvar(nvars,1);

% No reason to try to do factor tracking here
x = flush(x);

if size(At,2)~=length(b)
    At = At';
end

F = ([]);
top = 1;

if isvalidfield(K,'f')
    X = c(top:top+K.f-1)-At(top:top+K.f-1,:)*x;
    F = F + (X(:) == 0);
    top = top + K.f;
end

if isvalidfield(K,'l')
    X = c(top:top+K.l-1)-At(top:top+K.l-1,:)*x;
    F = F + (X(:)>=0);
    top = top + K.l;
end

if isvalidfield(K,'q')
    for i = 1:length(K.q)
        X = c(top:top+K.q(i)-1)-At(top:top+K.q(i)-1,:)*x;
        F = F + (cone(X(2:end),X(1)));
        top = top + K.q(i);
    end
end

if isvalidfield(K,'r')
    for i = 1:length(K.r)
        X = c(top:top+K.r(i)-1)-At(top:top+K.r(i)-1,:)*x;
        F = F + (rcone(X(3:end),X(2),X(1)));
        top = top + K.r(i);
    end
end

if isvalidfield(K,'s')
    for i = 1:length(K.s)
        [ix,iy,iv] = find([c(top:top+K.s(i)^2-1) At(top:top+K.s(i)^2-1,:)]);
        off = (ix-1)/(K.s(i)+1);
        if all(off == round(off))
            X = c(top:top+K.s(i)^2-1)-At(top:top+K.s(i)^2-1,:)*x;
            if isa(X,'sdpvar')
                F = F + (diag(reshape(X,K.s(i),K.s(i))) >= 0);
            else
                X
                i
                'silly data!'
            end
            top = top + K.s(i)^2;
        else
            X = c(top:top+K.s(i)^2-1)-At(top:top+K.s(i)^2-1,:)*x;
            X = reshape(X,K.s(i),K.s(i));
            X = (X+X')/2;
            F = F + (X >= 0);
            top = top + K.s(i)^2;
        end
    end
end

h = -b'*x;

function ok = isvalidfield(K,fld)
ok = 0;
if isfield(K,fld)
    s = getfield(K,fld);
    if prod(size(s))>0
        if s(1)>0
            ok = 1;
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%WARPPERSPECTIVE  Applies a perspective transformation to an image
%
%     dst = cv.warpPerspective(src, M)
%     dst = cv.warpPerspective(src, M, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Input image.
% * __M__ 3x3 transformation matrix, floating-point.
%
% ## Output
% * __dst__ Output image that has the size `DSize` (with
%   `size(dst,3) == size(src,3)`) and the same type as `src`.
%
% ## Options
% * __DSize__ Size of the output image `[w,h]`. Default `[0,0]` means using
%   the same size as the input `[size(src,2) size(src,1)]`.
% * __Interpolation__ interpolation method, default 'Linear'. One of:
%   * __Nearest__ nearest neighbor interpolation
%   * __Linear__ bilinear interpolation
%   * __Cubic__ bicubic interpolation
%   * __Lanczos4__ Lanczos interpolation over 8x8 neighborhood
% * __WarpInverse__ Logical flag to apply inverse perspective transform,
%   meaning that `M` is the inverse transformation (`dst -> src`).
%   default false
% * __BorderType__ Pixel extrapolation method. When 'Transparent', it means
%   that the pixels in the destination image corresponding to the "outliers"
%   in the source image are not modified by the function. default 'Constant'
%   * __Constant__ `iiiiii|abcdefgh|iiiiiii` with some specified `i`
%   * __Replicate__ `aaaaaa|abcdefgh|hhhhhhh`
%   * __Reflect__ `fedcba|abcdefgh|hgfedcb`
%   * __Reflect101__ `gfedcb|abcdefgh|gfedcba`
%   * __Wrap__ `cdefgh|abcdefgh|abcdefg`
%   * __Transparent__ `uvwxyz|abcdefgh|ijklmno`
%   * __Default__ same as 'Reflect101'
% * __BorderValue__ Value used in case of a constant border. default 0
% * __Dst__ Optional initial image for the output. If not set, it is
%   automatically created by the function. Note that it must match the
%   expected size `DSize` and the type of `src`, otherwise it is ignored and
%   recreated by the function. This option is only useful when
%   `BorderType=Transparent`, in which case the transformed image is drawn
%   onto the existing `Dst` without extrapolating pixels. Not set by default.
%
% The function cv.warpPerspective transforms the source image using the
% specified matrix:
%
%     dst(x,y) = src((M_11*x + M_12*y + M_13) / (M_31*x + M_32*y + M_33),
%                    (M_21*x + M_22*y + M_23) / (M_31*x + M_32*y + M_33))
%
% when the `WarpInverse` option is true. Otherwise, the transformation is first
% inverted with cv.invert and then put in the formula above instead of `M`.
%
% See also: cv.warpAffine, cv.remap, cv.resize, cv.getRectSubPix,
%  cv.perspectiveTransform, imtransform, imwarp
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%  mcmc object makes diagnostic plots from mcmc draws
%
%
%    Reference page in Doc Center
%       doc mcmc
%
%

author: jmaih


text:
function [state, options,optchanged] = myfun(options,state,flag)
persistent history
persistent cost
optchanged = false;

switch flag
 case 'init'
        history(:,:,1) = state.Population;
        cost(:,1) = state.Score;
    case {'iter','interrupt'}
        ss = size(history,3);
        history(:,:,ss+1) = state.Population;
        cost(:,ss+1) = state.Score;
    case 'done'
        ss = size(history,3);
        history(:,:,ss+1) = state.Population;
        cost(:,ss+1) = state.Score;
        save history.mat history cost
end

author: dynamicslab


text:
% Generate time points for evaluating spline.
%
%  The convex-combination of the endpoitns with
%  five controls points are 80 percent the last cpt
%  and 20 percent the control point after that.
%
% Input
%   nland: number of landmarks
%   neval: number of evaluations
%
% Output
%   s: the time points used to evaluate spline
function [s,lb,ub] = bspline_gen_s(nland,neval)
    if nargin < 2
       neval = 200;
    end

    lb = 2;
    ub = nland+1;
    epts = [lb ub];

    len = epts(2)-epts(1);
    int = len/(neval-1);
    s = epts(1):int:epts(2);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [ y, m, d, f, seed ] = ymdf_uniform_common ( y1, m1, d1, f1, ...
  y2, m2, d2, f2, seed )

%*****************************************************************************80
%
%% YMDF_UNIFORM_COMMON: random Common YMDF date between two given dates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the first YMDF date.
%
%    Input, integer Y2, M2, D2, real F2,
%    the second YMDF date.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, integer Y, M, D, real F, the
%    random YMDF date.
%
  jed1 = ymdf_to_jed_common ( y1, m1, d1, f1 );
  jed2 = ymdf_to_jed_common ( y2, m2, d2, f2 );

  [ jed, seed ] = r8_uniform_ab ( jed1, jed2, seed );

  [ y, m, d, f ] = jed_to_ymdf_common ( jed );

  return
end


text:
% run all example sequences
%
% Stefan Kroboth <stefan.kroboth@uniklinik-freiburg.de>

% avoid cwd madness
od = cd(fileparts(which(mfilename)));
addpath(genpath('../matlab'));

% list of sequences to execute
fid
gre
epi_rs

% just in case the change of the cwd persists in the parent shell
cd(od);

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function F = le(X,Y)
% Internal class for constraint lists

% Try to evaluate
try
    if isa(X,'constraint')
        % (z > w) < y
        try
            Z = Y - X.List{end};
        catch
            Y = reshape(Y,[],1);
            Z = Y - X.List{end};
        end
        F = X;
        F.List{end+1} = '<=';
        F.List{end+1} = Y;
        F.Evaluated{end+1} = Z;
        F.ConstraintID(end+1) = yalmip('ConstraintID');
        F.strict(end+1) = 0;
    else
        % x < (w > y)
        try
            Z = Y.List{1} - X;
        catch
            X = reshape(X,[],1);
            Z = Y.List{1} - X;
        end
        F = Y;
        F.List = {X,'<=',F.List{:}};
        F.Evaluated = {Z,F.Evaluated{:}};
        F.ConstraintID = [yalmip('ConstraintID') F.ConstraintID];
        F.strict = [1 F.strict];
    end
catch
    error(lasterr);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%kvsamp 'Sample an Image Along an Arbitrary Line (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vsamp.pane file
%
% Parameters:
% InputFile: i 'Input image ', required: 'input image'
% Integer: x1 'x ', default: 0: 'starting x coordinate for line'
% Integer: x2 'x ', default: 511: 'ending x coordinate for line'
% Integer: y1 'y ', default: 255: 'starting y coordinate for line'
% Integer: y2 'y ', default: 255: 'ending y coordinate for line'
% Integer: n 'Number of samples ', default: 1: 'number of samples'
% OutputFile: o 'Output image', required: 'output image'
%
% Example: o = kvsamp(i, {'i','';'x1',0;'x2',511;'y1',255;'y2',255;'n',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vsamp - Sample an Image Along an Arbitrary Line  (K1)
%
%  DESCRIPTION
% .B NOTE:
% Program will NOT execute!
% vsamp functionality is replaced by the DATAMANIP::ksampline operator.
% The following command line arguments have changed.
%
%      -x1 (starting x coordinate for line), is now -wc1
%      -x2 (ending x coordinate for line), is now -wc2
%      -y1 (starting y coordinate for line), is now -hc1
%      -y2 (ending y coordinate for line), is now -hc1
%
% the values may be be slightly different because vsamp used bilinear
% interpolation whereas ksampline uses inverse-distance weighting
% interpolation (which is mush easier to extend to 5D and mask'd data).
%
% Given an image of arbitrary size and a line through that image, sample the
% image along the line.  The sampling interval is independent of
% image size and line orientation. The value at each sample point is obtained
% by bilinear interpolation.
%
% The first sample point is precisely at the starting point and the final sample
% is precisely at the ending point.
%
%
%
%  EXAMPLES
% vsamp -i diffract.xv -o line_data.xv -x1 0 -y1 0 -x2 0 -y2 0 -n 400
%
% Will sample the image diffract.xv to obtain a one dimensional data array
% in line_data.xv.  The line goes from the top left corner to the bottom
% right. 400 evenly spaced samples will be obtained via bilinear interpolation.
%
%  "SEE ALSO"
% DATAMANIP::ksampline
%
%  RESTRICTIONS
% .I vsamp
% will only work on images of type FLOAT.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kvsamp(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvsamp(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'x1', 0;'x2', 511;'y1', 255;'y2', 255;'n', 1;'o', '__output'};
maxval={0,1,1,1,1,2,0};
minval={0,1,1,1,1,2,0};
istoggle=[0,1,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [exp,changed,rules] = exp_replaceall(exp,rules)
% Apply substitution rules to some expression.
% [Out-Expression,Changed,Remaining-Rules] = exp_replaceall(In-Expression, Rule/Rule-List)
%
% Each part of the expression receives at most a single replacement, and each rule is used at most
% once rules may contain blanks and placeholders.
%
% In:
%   In-Expression : some pure expression in which substitution shall take place
%   Rule-List     : cell array of substitution rules
%
% Out:
%   Out-Expression     : In-Expression with substitutions performed
%
%   Changed            : whether a replacement has taken place
%   Remaining-Rules    : the list of rules that have not been applied
%
% Notes:
%   Impure expressions degrade into pure expressions, when substitutions are performed on them (i.e.
%   they lose their attached value). Named patterns inside the right-hand-sides of substitution
%   rules are treated as if just their symbols were given instead.
%
% See also:
%   exp_replacerepeated, exp_match
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-19

if ~exp_beginfun('symbolic') return; end

[exp,changed,rules] = utl_replaceall(exp,rules);

exp_endfun;

repo: REST


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Superposed Hawkes Processes: various source_num, D = 5
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
run setup
D = 5;
K = [2, 5, 10];
N = 5;
nTest = 10;
rate = 0;

Err = zeros(4, length(K), N, nTest);

para = cell(length(K), nTest);
Seqs = cell(K(end), nTest);
models = cell(4, length(K), N, nTest);

alg.LowRank = 0;
alg.Sparse = 0;
alg.GroupSparse = 0;
alg.outer = 10;
alg.rho = 0.1;
alg.inner = 10;
alg.thres = 1e-2;
alg.Tmax = [];
alg.w = 1;

algmu = alg;
algmu.Sparse = 1;
algmu.alphaS = 0.1;
algmu.w = 1;

options.N = 20; % the number of sequences
options.Nmax = 200; % the maximum number of events per sequence
options.Tmax = 100; % the maximum size of time window
options.tstep = 0.1;
options.dt = 0.1;
options.M = 250;
options.GenerationNum = 200;

for n = 1:nTest
    A = rand(D, D);
    A(A<0.4) = 0;
    A = A + eye(D);
    A = 0.5*A./max(abs(eig(A)));
    L = zeros(K(end),1);
    for k = 1:K(end)
        para{k, n}.kernel = 'exp';
        para{k, n}.landmark = 0;
        para{k, n}.w = 1;
        para{k, n}.mu = zeros(D,1);
        id = randperm(D);
        para{k, n}.mu(id(1)) = 0.1+rand;
        para{k, n}.A = reshape(A, [D, 1, D]);
        Seqs{k, n} = SimulationFast_Thinning_ExpHP(para{k, n}, options);
        L(k) = length(Seqs{k, n});
    end

    for k = 1:length(K)%:-1:1
        model.A = rand(D,1,D)./(D^2);
        model.mu = rand(D,1)./D;
        model.kernel = 'exp';
        model.w = 1;
        model.landmark = 0;

        for nn = 1:N

            num = min([nn*(options.N/N); L]);
            seq_sum = Seqs{1, n}(1:num);
            seq_super = Seqs{1,n}(1:num);
            for kk = 2:K(k)
                seq_sum = [seq_sum, Seqs{kk, n}(1:num)];
                seq_super = SuperPosition(seq_super, ...
                                             Seqs{kk, n}(1:num), rate);
            end

%             seq_super = [];
%             for m = 1:K(k)
%                 seq_super_tmp = Seqs{1, n}(1:num);
%                 for kk = 2:K(k)
%                     tmp = Seqs{kk, n}(1:num);
%                     index = randperm(num);
%                     seq_super_tmp = SuperPosition(seq_super_tmp, ...
%                                             Seqs{kk, n}(index), rate);
%                 end
%                 seq_super = [seq_super, seq_super_tmp];
%             end

            models{1, k, nn, n} = model;
            models{1, k, nn, n} = Learning_WLS_HP( ...
                Seqs{1, n}(1:num), models{1, k, nn, n}, alg, 'singlesource' );

            models{2, k, nn, n} = model;
            models{2, k, nn, n} = Learning_WLS_HP( ...
                seq_sum, models{2, k, nn, n}, alg, 'singlesource' );

            model2 = model;
            model2.mu = rand(D, num*K(k))./D;
            models{3, k, nn, n} = model2;
            models{3, k, nn, n} = Learning_WLS_HP( ...
                seq_sum, models{3, k, nn, n}, algmu, 'multisource' );

            models{4, k, nn, n} = model;
            alg2 = alg;
            alg2.w = 1/sqrt(K(k));
            models{4, k, nn, n} = Learning_WLS_HP( ...
                seq_super, models{4, k, nn, n}, alg2, 'singlesource' );

            Err(1, k, nn, n) = norm(models{1, k, nn, n}.A(:) - ...
                A(:))/norm(A(:));

            Err(2, k, nn, n) = norm(models{2, k, nn, n}.A(:) - ...
                A(:))/norm(A(:));

            Err(3, k, nn, n) = norm(models{3, k, nn, n}.A(:) - ...
                A(:))/norm(A(:));

            Err(4, k, nn, n) = norm(models{4, k, nn, n}.A(:) - ...
                A(:))/norm(A(:));

%             A1 = models{1,k,nn,n}.A;
%             A2 = models{2,k,nn,n}.A;
%             imagesc([A,reshape(A1,[D,D]),reshape(A2,[D,D])]);
%             axis tight;axis equal;colorbar

        end
    end
end

save('Result_superpose_error_syn_WLS_D5.mat', 'Err', 'models', ...
'para', 'options', 'K');

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function [cellIndex] = findUIDcellIndex(uidString,cellString)
% WRITTEN BY MARTIN: CREATE HEADER!

cellIndex = find(strcmp(uidString,cellString));
if isempty(cellIndex)
    cellIndex = numel(cellString) + 1; % Not present in cellString, so create a new position in the cell for the new UID
end

end

author: mvallieres

repo: radiomics


text:
%%*****************************************************************************
%% sqlpsummary: print summary
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************************

function  sqlpsummary(info,ttime,infeas_org,printlevel)

iter        = info.iter;
obj         = info.obj;
gap         = info.gap;
relgap      = info.relgap;
prim_infeas = info.pinfeas;
dual_infeas = info.dinfeas;
termcode    = info.termcode;
reldist     = info.reldist;
resid       = info.resid;
dimacs      = info.dimacs;
totaltime   = info.cputime;
%%
preproctime = ttime.preproc; pcholtime = ttime.pchol; dcholtime = ttime.dchol;
predtime  = ttime.pred; pstep_predtime = ttime.pred_pstep; dstep_predtime = ttime.pred_pstep;
corrtime  = ttime.corr; pstep_corrtime = ttime.corr_pstep; dstep_corrtime = ttime.corr_dstep;
misctime  = ttime.misc;
%%
if (printlevel >= 2)
    fprintf('\n------------------------------------------------');
    fprintf('-------------------\n');
    fprintf(' number of iterations   = %2.0f\n',iter);
end
if (termcode <= 0)
    if (printlevel >=2)
        fprintf(' primal objective value = %- 9.8e\n',obj(1));
        fprintf(' dual   objective value = %- 9.8e\n',obj(2));
        fprintf(' gap := trace(XZ)       = %3.2e\n',gap);
        fprintf(' relative gap           = %3.2e\n',relgap);
        fprintf(' actual relative gap    = %3.2e\n',-diff(obj)/(1+sum(abs(obj))));
        if norm(infeas_org)
            fprintf(' rel. primal infeas (scaled problem)   = %3.2e\n',prim_infeas);
            fprintf(' rel. dual     "        "       "      = %3.2e\n',dual_infeas);
            fprintf(' rel. primal infeas (unscaled problem) = %3.2e\n',infeas_org(1));
            fprintf(' rel. dual     "        "       "      = %3.2e\n',infeas_org(2));
        else
            fprintf(' rel. primal infeas     = %3.2e\n',prim_infeas);
            fprintf(' rel. dual   infeas     = %3.2e\n',dual_infeas);
        end
        fprintf(' norm(X), norm(y), norm(Z) = %3.1e, %3.1e, %3.1e\n',...
            info.normX,info.normy,info.normZ);
        fprintf(' norm(A), norm(b), norm(C) = %3.1e, %3.1e, %3.1e\n',...
            info.normA,info.normb,info.normC);
    end
elseif (termcode == 1)
    if (printlevel >=2)
        fprintf(' residual of primal infeasibility      \n')
        fprintf(' certificate (y,Z)      = %3.2e\n',resid);
        fprintf(' reldist to infeas.    <= %3.2e\n',reldist);
    end
elseif (termcode == 2)
    if (printlevel >=2)
        fprintf(' residual of dual infeasibility        \n')
        fprintf(' certificate X          = %3.2e\n',resid);
        fprintf(' reldist to infeas.    <= %3.2e\n',reldist);
    end
end
if (printlevel >=2)
    fprintf(' Total CPU time (secs)  = %3.2f  \n',totaltime);
    fprintf(' CPU time per iteration = %3.2f  \n',totaltime/iter);
    fprintf(' termination code       = %2.0f\n',termcode);
    fprintf(' DIMACS: %.1e  %.1e  %.1e  %.1e  %.1e  %.1e\n',dimacs);
    fprintf('------------------------------------------------');
    fprintf('-------------------\n');
    if (printlevel > 3)
        fprintf(' Percentage of CPU time spent in various parts \n');
        fprintf('------------------------------------------------');
        fprintf('-------------------\n');
        fprintf(' preproc Xchol Zchol pred pred_steplen  corr corr_steplen  misc\n')
        tt = [preproctime, pcholtime, dcholtime, predtime, pstep_predtime, dstep_predtime];
        tt = [tt, corrtime, pstep_corrtime, dstep_corrtime, misctime];
        tt = tt/sum(tt)*100;
        fprintf('   %3.1f   %3.1f   %3.1f   %3.1f  %3.1f  %3.1f ',...
            tt(1),tt(2),tt(3),tt(4),tt(5),tt(6));
        fprintf('    %3.1f  %3.1f  %3.1f     %3.1f\n',tt(7),tt(8),tt(9),tt(10));
        fprintf('------------------------------------------------');
        fprintf('-------------------\n');
    end
end
%%*****************************************************************************

author: yu-jiang

repo: radpbook

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


text:
clear

if(isunix)
    executable = '"../../build/bin/FeatureExtraction"';
else
    executable = '"../../x64/Release/FeatureExtraction.exe"';
end

output = './processed_features/';

in_file = '../../samples/2015-10-15-15-14.avi';

command = sprintf('%s -f "%s" -out_dir "%s" -gaze -verbose', executable, in_file, output);

if(isunix)
    unix(command);
else
    dos(command);
end

%% Demonstrating reading the output files
[~, out_filename,~] = fileparts(in_file);
out_filename = sprintf('%s/%s.csv', output, out_filename);

% First read in the column names
tab = readtable(out_filename);
column_names = tab.Properties.VariableNames;

all_params  = dlmread(out_filename, ',', 1, 0);

gaze_angle_ids = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'gaze_angle_'));

gaze = all_params(:,gaze_angle_ids);

plot(gaze(:,1), 'DisplayName', 'Left - right');
hold on;
plot(gaze(:,2), 'DisplayName', 'Up - down');
xlabel('Frame') % x-axis label
ylabel('Angle in radians') % y-axis label
legend('show');
hold off;

author: TadasBaltrusaitis

repo: OpenFace


text:
function output = structureStats(structNum, planC)
%"structureStats"
%   Returns output, a struct containing the xyz coordinates for the
%   Center of mass(COM) of the mask of the structure, and the min
%   and max extent of the structure.
%
%   If a transformation matrix exists for a structure's associated scan,
%   the coordinates are transformed using the matrix.
%
%JRA 3/12/04
%
%Usage:
%   function output = structureStats(structNum, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC;
end

indexS = planC{end};

%Get scan associated with structure.
scanSet = getStructureAssociatedScan(structNum, planC);

%Get mask of structure.
uniformStr = getUniformStr(structNum, planC);

%Get r,c,s of voxels inside uniformStr.
[r,c,s] = find3d(uniformStr);
nVoxInStruct = length(r);
clear uniformStr;

%Get scan's original x,y,z coordinates, and it's transformation matrix.
if isempty(planC{indexS.scan}(scanSet).uniformScanInfo)
    planC = setUniformizedData(planC);
end
[xV, yV, zV] = getUniformScanXYZVals(planC{indexS.scan}(scanSet));
transM = getTransM(planC{indexS.scan}(scanSet), planC);
voxVol = abs(xV(2)-xV(1)) * abs(yV(2)-yV(1)) * abs(zV(2)-zV(1));

%Get the x,y,z coords of points in the structure, and transform.
structXV = xV(c); clear c xV;
structYV = yV(r); clear r yV;
structZV = zV(s); clear s zV;
[xT, yT, zT] = applyTransM(transM, structXV, structYV, structZV);
clear structXV structYV structZV

%Find the bounding box around the structure.
minX = min(xT);
minY = min(yT);
minZ = min(zT);
maxX = max(xT);
maxY = max(yT);
maxZ = max(zT);

yCOM = mean(yT);
xCOM = mean(xT);
zCOM = mean(zT);

output.min = min([[minX, minY, minZ];[maxX, maxY, maxZ]]);
output.max = max([[minX, minY, minZ];[maxX, maxY, maxZ]]);
output.COM = [xCOM, yCOM, zCOM];
output.vol = voxVol * nVoxInStruct;
output.name = planC{indexS.structures}(structNum).structureName;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function  opt = FindOpt(model,Population,BU,BD)
% Find the minimum of the surrogate

%------------------------------- 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".
%--------------------------------------------------------------------------

    [~,I]  = sort(Population.objs,'ascend');
    init   = 1;
    preObj = rbf_predict(model,Population.decs,Population(I(init)).decs);
    while isnan(preObj)
        init   = init + 1;
        preObj = rbf_predict(model,Population.decs,Population(I(init)).decs);
    end
    opt = fmincon(@(x)rbf_predict(model,Population.decs,x),Population(I(init)).decs,[],[],[],[],BD,BU,[],optimoptions('fmincon','Display','off'));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%tstoolbox/mex/takens_estimator
%   Syntax:
%
%     * D = takens_estimator(pointset, query_indices, relative_range,
%       exclude)
%     * D = takens_estimator(atria, pointset, query_indices,
%       relative_range, exclude)
%
%   Input arguments:
%
%     * atria - output of nn_prepare for pointset (optional)
%       (cf. Section )
%     * pointset - a N by D double matrix containing the coordinates of
%       the point set, organized as N points of dimension D
%     * query_indices - query points are taken out of the pointset,
%       query_indices is a vector of length R which contains the indices
%       of the query points (indices may vary from 1 to N)
%     * relative_range - search radius, relative to attractor diameter (0
%       < relative_range < 1)
%     * exclude - in case the query points are taken out of the pointset,
%       exclude specifies a range of indices which are omitted from
%       search. For examples if the index of the query point is 124 and
%       exclude is set to 3, points with indices 121 to 127 are omitted
%       from search. Using exclude = 0 means: exclude self-matches
%
%   Output arguments:
%
%     * D - scalar value, estimation of correlation dimension
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% edge_response_2d.m
% examine variation in edge response for 2D CT
% Copyright 2011-02-23, Jeff Fessler, University of Michigan

if ~isvar('sino'), printm 'setup geometry, image, sinogram'
	down = 2;
	ig = image_geom('nx', 512, 'fov', 50, 'down', down);
	ig.mask = ig.circ > 0;
	sg = sino_geom('ge1', 'units', 'cm', 'down', down);

	% system object
	A = Gtomo2_dscmex(sg, ig);

	% image

	steps = [1 2 4 8] * 10;
        nstep = numel(steps);
	ell = [0 0 15 10 0 1000];
	leg = {};
	for ii=1:nstep
		ell = [ell; [-15+ii*6 0 2 2 0 steps(ii)]];
		leg{end+1} = sprintf('%d HU', steps(ii));
	end
	xtrue = ellipse_im(ig, ell, 'oversample', 2);

	im plc 2 2
	clim = [800 1200];
	im(1, xtrue, 'x', clim), cbar

%	sino = ellipse_sino(sg, ell, 'oversample', 2);
	sino = A * xtrue; % cheat
	im(2, sino, 'sino'), cbar
prompt
end

if ~isvar('fbp'), printm 'fbp 2d fan-beam reconstruction'
	tmp = fbp2(sg, ig);
	fbp = fbp2(sino, tmp, 'window', 'hanning,0.5');
	im(3, fbp, 'FBP', clim), cbar
prompt
end

if ~isvar('kappa'), printm 'kappa'
	wi = ones(size(sino)); % no statistical weighting
	kappa = ig.mask;
end

% use local psf to help select beta
if ~isvar('R'), printm 'R'
	f.l2b = 4; % maybe a bit too big, but ok for now
	f.delta = 10;
	R = Reg1(kappa, 'beta', 2^f.l2b, 'pot_arg', {'hyper3', f.delta});
	qpwls_psf(A, R, 1, ig.mask, Gdiag(wi), 'loop', 1);
prompt
end

if ~isvar('xpwls'), printm 'iterative reconstruction'
	Ab = Gblock(A, 41); % 41 subsets
	f.niter = 20;
	xpwls = pwls_sps_os(fbp(ig.mask), sino, wi, Ab, R, f.niter);
	xpwls = ig.embed(xpwls(:,end));
	im(4, xpwls, clim), cbar
prompt
end

im plc 3 3
pl = @(i) subplot(330+i);

im(1, 'notick', ig.x, ig.y, fbp, clim, 'FBP')
im(4, 'notick', ig.x, ig.y, xpwls, clim, 'PWLS')

tmp1 = fbp(:,end/2+1);
tmp2 = xpwls(:,end/2+1);

pl(2)
plot(ig.x, tmp1)
axis([-20 20 0 1200])
title 'Profile'
xtick([-1 0 1] * 15), ytick([0 1000])

pl(5)
plot(ig.x, tmp2)
axis([-20 20 0 1200])
title 'Profile'
xtick([-1 0 1] * 15), ytick([0 1000])

xx = outer_sum(ig.x, -ell(2:end,1));
pl(3)
plot(xx, (tmp1-1000) * (1 ./ steps), '.-')
ax = [-3 -1 -0.25 1.25];
axis(ax), legend(leg{:}, 'location', 'seo')
title 'Edge response'
xtick([-3 -2 -1]), ytick([0 1])

pl(6)
plot(xx, (tmp2-1000) * (1 ./ steps), '.-')
axis(ax), legend(leg{:}, 'location', 'seo')
title 'Edge response'
xtick([-3 -2 -1]), ytick([0 1])

% ir_savefig eps_c edge_response_2d

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%compute the xcoordinate of the head in the central coordinate system

function [data,units]=compute_xcentralhead_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
xcentralhead_mm=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    xcentralhead_mm{1,i}=trx(larva).xhead_mm-trx(larva).xcentral_mm;
end

units=parseunits('mm');
data=xcentralhead_mm;

author: kristinbranson

repo: JAABA


text:
function kernDisplay(kern, varargin)

% KERNDISPLAY Display the parameters of the kernel.
% FORMAT
% DESC displays the parameters of the kernel and the kernel type to
% the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO modelDisplay
%
% COPYRIGHT Neil D. Lawrence, 2006, 2005, 2004

% KERN

fhandle = str2func([kern.type 'KernDisplay']);
fhandle(kern, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function title = p26_title ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = r8_besy0 ( x )

%*****************************************************************************80
%
%% R8_BESY0 evaluates the Bessel function Y of order 0 of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the Bessel function Y of order 0 of X.
%
  persistent alnhaf
  persistent by0cs
  persistent nty0
  persistent twodpi
  persistent xsml

  alnhaf = -0.69314718055994530941723212145818;
  twodpi = 0.636619772367581343075535053490057;

  if ( isempty ( nty0 ) )
    by0cs = [ ...
      -0.1127783939286557321793980546028E-01, ...
      -0.1283452375604203460480884531838, ...
      -0.1043788479979424936581762276618, ...
      +0.2366274918396969540924159264613E-01, ...
      -0.2090391647700486239196223950342E-02, ...
      +0.1039754539390572520999246576381E-03, ...
      -0.3369747162423972096718775345037E-05, ...
      +0.7729384267670667158521367216371E-07, ...
      -0.1324976772664259591443476068964E-08, ...
      +0.1764823261540452792100389363158E-10, ...
      -0.1881055071580196200602823012069E-12, ...
      +0.1641865485366149502792237185749E-14, ...
      -0.1195659438604606085745991006720E-16, ...
      +0.7377296297440185842494112426666E-19, ...
      -0.3906843476710437330740906666666E-21, ...
      +0.1795503664436157949829120000000E-23, ...
      -0.7229627125448010478933333333333E-26, ...
      +0.2571727931635168597333333333333E-28, ...
      -0.8141268814163694933333333333333E-31 ]';

    nty0 = r8_inits ( by0cs, 19, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r8_mach ( 3 ) );

  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BESY0 - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R8_BESY0 - Fatal error!' )
  elseif ( x <= xsml )
    y = 0.0;
    value = twodpi * ( alnhaf + log ( x ) ) * r8_besj0 ( x ) ...
      + 0.375 + r8_csevl ( 0.125 * y - 1.0, by0cs, nty0 );
  elseif ( x <= 4.0 )
    y = x * x;
    value = twodpi * ( alnhaf + log ( x ) ) * r8_besj0 ( x ) ...
      + 0.375 + r8_csevl ( 0.125 * y - 1.0, by0cs, nty0 );
  else
    [ ampl, theta ] = r8_b0mp ( x );
    value = ampl * sin ( theta );
  end

  return
end


text:
function value = square_volume ( a, b )

%*****************************************************************************80
%
%% SQUARE_VOLUME: volume of a square in 2D.
%
%  Discussion:
%
%    The integration region is defined as:
%      A(1) <= X <= B(1)
%      A(2) <= Y <= B(2)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(2), B(2), the lower and upper limits.
%
%    Output, real VALUE, the volume.
%
  value = ( b(1) - a(1) ) * ( b(2) - a(2) );

  return
end


text:
function h = plotaxes(a, m, type, col, varargin)
% PLOTAXES  Plot axes of 3D coordinate system
%
% PLOTAXES(A, M, TYPE, COL)
%
%   A is a (3,3)-matrix with 3 column vectors. Each vector gives the
%   direction of one axis of a Coordinate system.
%
%   M is a 3-vector with the origin of coordinates. By default, M=[0,0,0].
%
%   TYPE is a cell array with the line types to plot (e.g. '--', ':'... see
%   "help plot" for more details). By default, type='-', i.e. a solid line.
%
%   COL is a cell array with the colours of each axis. By default,
%   col={'b', 'r', 'g'}. If only one string is provided, all axes are
%   displayed with the same colour.
%
% PLOTAXES(A, M, TYPE, COL, ...)
%
%   ... Any extra input arguments will be passed to the plot function.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010-2011 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(1, Inf, nargin, 'struct'));
error(nargoutchk(0, 1, nargout, 'struct'));

% defaults
if (nargin < 2 || isempty(m))
    m = [0 0 0];
end
if (nargin < 3 || isempty(type))
    type = {'-', '-', '-'};
end
if (nargin < 4 || isempty(col))
    col = {'b', 'r', 'g'};
end
if (~iscell(type) || length(type) == 1)
    if iscell(type)
        type = type{:};
    end
    type = {type, type, type};
end
if (~iscell(col) || length(col) == 1)
    if iscell(col)
        col = col{:};
    end
    col = {col, col, col};
end

% save hold state for later
PHOLD = ishold;

% plot axes
for I = 1:size(a, 2)
    h = plot3(m(1) + [0, a(1, I)], ...
        m(2) + [0, a(2, I)], m(3) + [0, a(3, I)], [ type{I} col{I} ], ...
        varargin{:});
    hold on
end

% recovert hold state
if PHOLD
    hold on
else
    hold off
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ S, f, Serr ]= mtspectrumc_unequal_length_trials( data, movingwin, params, sMarkers )

% This routine computes the multi-taper spectrum for a given set of unequal length segments. It is
% based on modifications to the Chronux routines. The segments are continuously structured in the
% data matrix, with the segment boundaries given by markers. Below,
% movingwin is used in a non-overlaping way to partition each segment into
% various windows. Th spectrum is evaluated for each window, and then the
% window spectrum estimates averaged. Further averaging is conducted by
% repeating the process for each segment.
%
% Inputs:
%
%   data = data( samples, channels )- here segments must be stacked
%   as explained in the email
%   movingwin = [window winstep] i.e length of moving
%              window and step size. Note that units here have
%              to be consistent with units of Fs. If Fs=1 (ie normalized)
%              then [window winstep]should be in samples, or else if Fs is
%              unnormalized then they should be in time (secs).
%   sMarkers = N x 2 array of segment start & stop marks. sMarkers(n, 1) = start
%           sample index; sMarkers(n,2) = stop sample index for the nth segment
%   params = see Chronux help on mtspecgramc
%
% Output:
%
%       S       frequency x channels
%       f       frequencies x 1
%       Serr    (error bars) only for err(1)>=1
%
%

iwAvg = 1; % 0=no weighted average, 1=weighted average
debug = 0; % will display intermediate calcs.

if nargin < 2; error('Unequal length trials:: Need data and window parameters'); end;
if nargin < 3; params=[]; end;
if isempty( sMarkers ), error( 'Unequal length trials:: Need Markers...' ); end
[ tapers, pad, Fs, fpass, err, trialave, params ] = getparams( params );
if nargout > 2 && err(1)==0;
%   Cannot compute error bars with err(1)=0. change params and run again.
    error('Unequal length trials:: When Serr is desired, err(1) has to be non-zero.');
end;

% Set moving window parameters to no-overlapping
if abs(movingwin(2) - movingwin(1)) >= 1e-6, disp( 'avgSpectrum:: Warming: Window parameters for averaging should be non-overlapping. Set movingwin(2) = movingwin(1).' ); end

wLength = round( Fs * movingwin(1) ); % number of samples in window
wStep = round( movingwin(2) * Fs ); % number of samples to step through

% Check whether window lengths satify segment length > NW/2
if ( wLength < 2*tapers(1) ), error( 'avgSpectrum:: movingwin(1) > 2*tapers(1)' ); end

% Left align segment markers for easier coding
sM = ones( size( sMarkers, 1 ), 2 );
sM( :, 2 ) = sMarkers( :, 2 ) - sMarkers( :, 1 ) + 1;

% min-max segments
Nmax = max( sM(:,2) ); Nmin = min( sM(:,2) );
if ( Nmin < 2*tapers(1) ), error( 'avgSpectrum:: Smallest segment length > 2*tapers(1). Change taper settings' ); end

% max time-sample length will be the window length.
nfft = 2^( nextpow2( wLength ) + pad );
[ f, findx ] = getfgrid( Fs, nfft, fpass);

% Precompute all the tapers
sTapers = tapers;
sTapers = dpsschk( sTapers, wLength, Fs ); % compute tapers for window length

nChannels = size( data, 2 );
nSegments = size( sMarkers, 1 );

if debug
    disp( ['Window Length = ' num2str(wLength)] );
    disp( ['Window Step = ' num2str(wStep)] );
    disp( ' ' );
end

s = zeros( length(f), nChannels );
serr = zeros( 2, length(f), nChannels );
S = zeros( length(f), nChannels );
Serr = zeros( 2, length(f), nChannels );
nWins = 0;
for sg = 1 : nSegments
    % Window lengths & steps fixed above
    % For the given segment, compute the positions & number of windows
    N = sM(sg,2);
    wStartPos = 1 : wStep : ( N - wLength + 1 );
    nWindows = length( wStartPos );
    if nWindows
        nWins = nWins + nWindows; % for averaging purposes

        w=zeros(nWindows,2);
        for n = 1 : nWindows
            w(n,:) = [ wStartPos(n), (wStartPos(n) + wLength - 1) ]; % nWindows x 2. just like segment end points
        end

        % Shift window limits back to original sample-stamps
        w(:, 1) = w(:,1) + (sMarkers( sg, 1 ) - 1);
        w(:, 2) = w(:,2) + (sMarkers( sg, 1 ) - 1);

        if debug
            disp( ['Segment Start/Stop = ' num2str( w(1,1) ) ' ' num2str( w(end,2) ) ] );
            disp( ['Min / Max Window Positions = ' num2str( min(w(:,1)) ) ' ' num2str( max(w(:,1)) ) ] );
            disp( ['Total Number of Windows = ' num2str(nWindows) ]);
            disp( ' ' );
        end

        % Pile up window segments similar to segment pileup
        wData = zeros( wLength, nChannels, nWindows ); %initialize to avoid fragmentation
        for n = 1:nWindows
            %wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ), 'constant' );
            wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ) );
        end

        % J1 = frequency x taper x nWindows
        % J2 = frequency x taper x nWindows x nChannels
        J2 = zeros( length(f), tapers(2), nWindows, nChannels ); J2 = complex( J2, J2 );
        for c = 1 : nChannels
            J1 = mtfftc( squeeze(wData( :, c, : )), sTapers, nfft, Fs ); % FFT for the tapered data
            J2( :, :, :, c ) = J1(findx,:,:);
        end
        % J2 = frequency x taper x nWindows x nChannels
        % Inner mean = Average over tapers => frequency x nWindows x nChannels
        % Outer mean = Average over windows => frequency x nChannels
        dim1 = [length(f), nWindows, nChannels];
        dim2 = [length(f), nChannels];
        % s = frequency x nChannels
        s = reshape( squeeze( mean( reshape( squeeze( mean( conj(J2).*J2, 2 ) ), dim1), 2 ) ), dim2 );

        % Now treat the various "windowed data" as "trials"
        % serr = 2 x frequency x channels. Output from specerr = 2 x frequency x 1
        for c = 1 : nChannels
            serr( :, :, c ) = specerr( squeeze( s(:, c ) ), squeeze( J2(:,:,:, c ) ), err, 1 );
        end

        if iwAvg
            % Segment Weighted error estimates.
            S = S + nWindows*s;
            Serr = Serr + nWindows*serr;
        else
            S = S + s;
            Serr = Serr + serr;
        end

    else
        if debug, disp(['avgSpectrum:: Zero windows for segment: ' num2str(sg) ]); end
    end
end

% Segment Weighted error estimates.
% Only over those that had non-zero windows
if nWins && iwAvg
    S=S/nWins; Serr=Serr/nWins;
end
if ~nWins
    if debug, disp(['avgCoherence:: No segment long enough with movingwin parameters found. Reduce movingwin.' ]); end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Note, in order to run this script, you will need
% the package sparsify v0.5 by Thomas Blumensath
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Signal space
N = 1024;
% Number of measurements
M = 100;
% Sparsity level
K = 8;

% solver instance
% solver = spx.pursuit.single.MatchingPursuit(Phi);
% Solve the sparse recovery problem
% result = solver.solve(y);
max_residual_norm = 1e-6;
options.max_residual_norm = max_residual_norm;
tstart = tic;
T = 1000;
t1 = 0;
t2 = 0;
for i=1:T
    if mod(i, 10) == 0
        fprintf('.');
    end
    % Construct the signal generator.
    gen  = spx.data.synthetic.SparseSignalGenerator(N, K);
    % Generate bi-uniform signals
    x = gen.biUniform(1, 2);
    % Sensing matrix
    Phi = spx.dict.simple.gaussian_mtx(M, N);
    % Measurement vectors
    y = Phi * x;
    tic;
    if 0
        result = spx.pursuit.single.mp(Phi, y,options);
    else
        [s, err_mse, iter_time] = greed_mp(y, Phi, N, ...
            'stopCrit', 'mse', 'stopTol', max_residual_norm^2, ...
            'maxIter', 4 * N);
    end
    % result.iterations
    t1 = t1 + toc;
    % norm(result.r) / norm(y)
    tic;
    result = spx.fast.mp(Phi, y,options);
    t2 = t2 + toc;
end
fprintf('\n');
fprintf('Time taken by MATLAB implementation: %0.3f\n', t1);
fprintf('Time taken by C++ implementation: %0.3f\n', t2);
fprintf('Gain: %0.2f\n', t1/t2);

author: indigits

repo: sparse-plex


text:
function [ y2, m2, d2, f2 ] = ymdf_prev_republican ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_PREV_REPUBLICAN returns the Republican YMDF date of the previous day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    yesterday's YMDF date.
%
  y2 = y1;
  m2 = m1;
  d2 = d1 - 1;
  f2 = f1;

  [ y2, m2, d2 ] = day_borrow_republican ( y2, m2, d2 );

  return
end


text:
% Fig. 8.10   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

numG=1;
denG=[1 0 0];

numD=[1 .2];
denD=[1 2];

num=conv(numG,numD);
den=conv(denG,denD);
poles=roots(den);
zeros=roots(num);

K1=0:.05:1.22;
K2=[1.25 1.28];  % K for break-in and break-away points
K3=1.5:5:100;
K=[K1 K2 K3];
Ko=.81;

r=rlocus(num,den,K);
ro=rlocus(num,den,Ko);

plot(r,'-'),
axis('square')
axis([-2.5 .5 -1.5 1.5])
hold on
plot(ro,'k*')
plot(-.2,0,'o')
plot(-2,0,'x')
plot(0,.01,'x')
plot(0,-.01,'x')
title('Fig. 8.10  s-plane locus vs. K')
xlabel('Re(s)')
ylabel('Im(s)')
text(-2.3,-.7,'*  K_c = 0.81')
nicegrid
hold off

author: Sable

repo: mcbench-benchmarks


text:
function row_new = catalan_row_next ( ido, n, row )

%*****************************************************************************80
%
%% CATALAN_ROW computes row N of Catalan's triangle.
%
%  Example:
%
%    I\J 0   1   2   3   4   5   6
%
%    0   1
%    1   1   1
%    2   1   2   2
%    3   1   3   5   5
%    4   1   4   9  14  14
%    5   1   5  14  28  42  42
%    6   1   6  20  48  90 132 132
%
%  Recursion:
%
%    C(0,0) = 1
%    C(I,0) = 1
%    C(I,J) = 0 for I < J
%    C(I,J) = C(I,J-1) + C(I-1,J)
%    C(I,I) is the I-th Catalan number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IDO, indicates whether this is a call for
%    the 'next' row of the triangle.
%    IDO = 0, this is a startup call.  Row N is desired, but
%    presumably this is a first call, or row N-1 was not computed
%    on the previous call.
%    IDO = 1, this is not the first call, and row N-1 was computed
%    on the previous call.  In this case, much work can be saved
%    by using the information from the previous values of IROW
%    to build the next values.
%
%    Input, integer N, the index of the row of the triangle desired.
%
%    Input, integer ROW(1:N), the previous row of coefficients, if
%    IDO = 1.
%
%    Output, integer ROW_NEW(1:N+1), the next row of coefficients.
%
  if ( n < 0 )
    row_new = [];
    return
  end

  if ( ido == 0 )

    row_new(1) = 1;
    row_new(2:n+1) = 0;

    for i = 1 : n

      row_new(1) = 1;

      for j = 1 : i-1
        row_new(j+1) = row_new(j+1) + row_new(j);
      end

      row_new(i+1) = row_new(i);

    end

  else

    row_new(1:n) = row(1:n);
    row_new(n+1) = 0;

    row_new(1) = 1;

    for j = 1 : n-1
      row_new(j+1) = row_new(j+1) + row_new(j);
    end

    if ( 1 <= n )
      row_new(n+1) = row_new(n);
    end

  end

  return
end


text:
function [in] = dm2in(dm)
% Convert length from decimeters to inches.
% I've still never used a decimeter.  Who uses this unit?
% Chad A. Greene 2012
in = dm/.254;


text:
classdef math_model_opf_acci < mp.math_model_opf_acc
%MP.MATH_MODEL_OPF_ACCI  MATPOWER mathematical model for AC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_ACCI ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

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

%     properties
%     end

    methods
        function tag = form_tag(obj)
            tag = 'acci';
        end

        function name = form_name(obj)
            name = 'AC-cartesian-current';
        end

        function add_node_balance_constraints(obj, nm, dm, mpopt)
            %% power balance constraints
            nn = nm.node.N;             %% number of nodes
            fcn_mis = @(x)obj.nodal_current_balance_fcn(x, nm);
            hess_mis = @(x, lam)obj.nodal_current_balance_hess(x, lam, nm);
            obj.add_nln_constraint({'rImis', 'iImis'}, [nn;nn], 1, fcn_mis, hess_mis);
        end

        function [lam_p, lam_q] = node_power_balance_prices(obj, nm)
            %% shadow prices on node power balance
            nne = obj.get_idx('nle');

            %% convert current balance shadow prices to equivalent lam_p and lam_q
            %% P + jQ = (Vr + jVi) * (M - jN)
            %% M = (Vr P + Vi Q) / (Vr^2 + Vi^2)
            %% N = (Vi P - Vr Q) / (Vr^2 + Vi^2)
            %% lam_p = df/dP = df/dM * dM/dP + df/dN + dN/dP
            %% lam_q = df/dQ = df/dM * dM/dQ + df/dN + dN/dQ
            V = nm.soln.v;
            lambda = obj.soln.lambda;
            VV = V ./ (V .* conj(V));   %% V / vm^2
            VVr = real(VV);
            VVi = imag(VV);
            lamM = lambda.eqnonlin(nne.i1.rImis:nne.iN.rImis);
            lamN = lambda.eqnonlin(nne.i1.iImis:nne.iN.iImis);
            lam_p = (VVr.*lamM + VVi.*lamN);
            lam_q = (VVi.*lamM - VVr.*lamN);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


% totalNum = 1867290;
%
% loopCompress(0.1, '/home/yh/mapModel/2018/08.03/weightVector.txt', '/home/yh/mapModel/2018/08.03/visMatrix/', 50, totalNum, 2040, '/home/yh/mapModel/2018/08.03/gurobi_compress_0.2/');
%
% ratio = salientNumCnt( '/home/yh/mapModel/2018/08.03/gurobi_compress_0.2/', totalNum ) / totalNum


text:
function build_trwsMex
% build_trwsMex builds package trwsMex
%
% Anton Osokin (firstname.lastname@gmail.com), 24.09.2014

mex trwsMex.cpp src/ordering.cpp src/MRFEnergy.cpp src/treeProbabilities.cpp src/minimize.cpp -output trwsMex -largeArrayDims

author: aosokin


text:
function json=savejson(rootname,obj,varargin)
%
% json=savejson(rootname,obj,filename)
%    or
% json=savejson(rootname,obj,opt)
% json=savejson(rootname,obj,'param1',value1,'param2',value2,...)
%
% convert a MATLAB object (cell, struct or array) into a JSON (JavaScript
% Object Notation) string
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%            created on 2011/09/09
%
% $Id: savejson.m 415 2013-10-07 16:38:31Z fangq $
%
% input:
%      rootname: name of the root-object, if set to '', will use variable name
%      obj: a MATLAB object (array, cell, cell array, struct, struct array)
%      filename: a string for the file name to save the output JSON data
%      opt: a struct for additional options, use [] if all use default
%        opt can have the following fields (first in [.|.] is the default)
%
%        opt.FileName [''|string]: a file name to save the output JSON data
%        opt.FloatFormat ['%.10g'|string]: format to show each numeric element
%                         of a 1D/2D array;
%        opt.ArrayIndent [1|0]: if 1, output explicit data array with
%                         precedent indentation; if 0, no indentation
%        opt.ArrayToStruct[0|1]: when set to 0, savejson outputs 1D/2D
%                         array in JSON array format; if sets to 1, an
%                         array will be shown as a struct with fields
%                         "_ArrayType_", "_ArraySize_" and "_ArrayData_"; for
%                         sparse arrays, the non-zero elements will be
%                         saved to _ArrayData_ field in triplet-format i.e.
%                         (ix,iy,val) and "_ArrayIsSparse_" will be added
%                         with a value of 1; for a complex array, the
%                         _ArrayData_ array will include two columns
%                         (4 for sparse) to record the real and imaginary
%                         parts, and also "_ArrayIsComplex_":1 is added.
%        opt.ParseLogical [0|1]: if this is set to 1, logical array elem
%                         will use true/false rather than 1/0.
%        opt.NoRowBracket [1|0]: if this is set to 1, arrays with a single
%                         numerical element will be shown without a square
%                         bracket, unless it is the root object; if 0, square
%                         brackets are forced for any numerical arrays.
%        opt.ForceRootName [0|1]: when set to 1 and rootname is empty, savejson
%                         will use the name of the passed obj variable as the
%                         root object name; if obj is an expression and
%                         does not have a name, 'root' will be used; if this
%                         is set to 0 and rootname is empty, the root level
%                         will be merged down to the lower level.
%        opt.Inf ['"$1_Inf_"'|string]: a customized regular expression pattern
%                         to represent +/-Inf. The matched pattern is '([-+]*)Inf'
%                         and $1 represents the sign. For those who want to use
%                         1e999 to represent Inf, they can set opt.Inf to '$11e999'
%        opt.NaN ['"_NaN_"'|string]: a customized regular expression pattern
%                         to represent NaN
%        opt.JSONP [''|string]: to generate a JSONP output (JSON with padding),
%                         for example, if opt.JSON='foo', the JSON data is
%                         wrapped inside a function call as 'foo(...);'
%        opt.UnpackHex [1|0]: conver the 0x[hex code] output by loadjson
%                         back to the string form
%        opt can be replaced by a list of ('param',value) pairs. The param
%        string is equivallent to a field in opt.
% output:
%      json: a string in the JSON format (see http://json.org)
%
% examples:
%      a=struct('node',[1  9  10; 2 1 1.2], 'elem',[9 1;1 2;2 3],...
%           'face',[9 01 2; 1 2 3; NaN,Inf,-Inf], 'author','FangQ');
%      savejson('mesh',a)
%      savejson('',a,'ArrayIndent',0,'FloatFormat','\t%.5g')
%
% license:
%     BSD license, see LICENSE_BSD.txt files for details
%
% -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==1)
   varname=inputname(1);
   obj=rootname;
   if(isempty(varname))
      varname='root';
   end
   rootname=varname;
else
   varname=inputname(2);
end
if(length(varargin)==1 && ischar(varargin{1}))
   opt=struct('FileName',varargin{1});
else
   opt=varargin2struct(varargin{:});
end
opt.IsOctave=exist('OCTAVE_VERSION');
rootisarray=0;
rootlevel=1;
forceroot=jsonopt('ForceRootName',0,opt);
if((isnumeric(obj) || islogical(obj) || ischar(obj) || isstruct(obj) || iscell(obj)) && isempty(rootname) && forceroot==0)
    rootisarray=1;
    rootlevel=0;
else
    if(isempty(rootname))
        rootname=varname;
    end
end
if((isstruct(obj) || iscell(obj))&& isempty(rootname) && forceroot)
    rootname='root';
end
json=obj2json(rootname,obj,rootlevel,opt);
if(rootisarray)
    json=sprintf('%s\n',json);
else
    json=sprintf('{\n%s\n}\n',json);
end

jsonp=jsonopt('JSONP','',opt);
if(~isempty(jsonp))
    json=sprintf('%s(%s);\n',jsonp,json);
end

% save to a file if FileName is set, suggested by Patrick Rapin
if(~isempty(jsonopt('FileName','',opt)))
    fid = fopen(opt.FileName, 'wt');
    fwrite(fid,json,'char');
    fclose(fid);
end

%%-------------------------------------------------------------------------
function txt=obj2json(name,item,level,varargin)

if(iscell(item))
    txt=cell2json(name,item,level,varargin{:});
elseif(isstruct(item))
    txt=struct2json(name,item,level,varargin{:});
elseif(ischar(item))
    txt=str2json(name,item,level,varargin{:});
else
    txt=mat2json(name,item,level,varargin{:});
end

%%-------------------------------------------------------------------------
function txt=cell2json(name,item,level,varargin)
txt='';
if(~iscell(item))
        error('input is not a cell');
end

dim=size(item);
len=numel(item); % let's handle 1D cell first
padding1=repmat(sprintf('\t'),1,level-1);
padding0=repmat(sprintf('\t'),1,level);
if(len>1)
    if(~isempty(name))
        txt=sprintf('%s"%s": [\n',padding0, checkname(name,varargin{:})); name='';
    else
        txt=sprintf('%s[\n',padding0);
    end
elseif(len==0)
    if(~isempty(name))
        txt=sprintf('%s"%s": []',padding0, checkname(name,varargin{:})); name='';
    else
        txt=sprintf('%s[]',padding0);
    end
end
for i=1:len
    txt=sprintf('%s%s%s',txt,padding1,obj2json(name,item{i},level+(len>1),varargin{:}));
    if(i<len) txt=sprintf('%s%s',txt,sprintf(',\n')); end
end
if(len>1) txt=sprintf('%s\n%s]',txt,padding0); end

%%-------------------------------------------------------------------------
function txt=struct2json(name,item,level,varargin)
txt='';
if(~isstruct(item))
	error('input is not a struct');
end
len=numel(item);
padding1=repmat(sprintf('\t'),1,level-1);
padding0=repmat(sprintf('\t'),1,level);
sep=',';

if(~isempty(name))
    if(len>1) txt=sprintf('%s"%s": [\n',padding0,checkname(name,varargin{:})); end
else
    if(len>1) txt=sprintf('%s[\n',padding0); end
end
for e=1:len
  names = fieldnames(item(e));
  if(~isempty(name) && len==1)
        txt=sprintf('%s%s"%s": {\n',txt,repmat(sprintf('\t'),1,level+(len>1)), checkname(name,varargin{:}));
  else
        txt=sprintf('%s%s{\n',txt,repmat(sprintf('\t'),1,level+(len>1)));
  end
  if(~isempty(names))
    for i=1:length(names)
	txt=sprintf('%s%s',txt,obj2json(names{i},getfield(item(e),...
             names{i}),level+1+(len>1),varargin{:}));
        if(i<length(names)) txt=sprintf('%s%s',txt,','); end
        txt=sprintf('%s%s',txt,sprintf('\n'));
    end
  end
  txt=sprintf('%s%s}',txt,repmat(sprintf('\t'),1,level+(len>1)));
  if(e==len) sep=''; end
  if(e<len) txt=sprintf('%s%s',txt,sprintf(',\n')); end
end
if(len>1) txt=sprintf('%s\n%s]',txt,padding0); end

%%-------------------------------------------------------------------------
function txt=str2json(name,item,level,varargin)
txt='';
if(~ischar(item))
        error('input is not a string');
end
item=reshape(item, max(size(item),[1 0]));
len=size(item,1);
sep=sprintf(',\n');

padding1=repmat(sprintf('\t'),1,level);
padding0=repmat(sprintf('\t'),1,level+1);

if(~isempty(name))
    if(len>1) txt=sprintf('%s"%s": [\n',padding1,checkname(name,varargin{:})); end
else
    if(len>1) txt=sprintf('%s[\n',padding1); end
end
isoct=jsonopt('IsOctave',0,varargin{:});
for e=1:len
    if(isoct)
        val=regexprep(item(e,:),'\\','\\');
        val=regexprep(val,'"','\"');
        val=regexprep(val,'^"','\"');
    else
        val=regexprep(item(e,:),'\\','\\\\');
        val=regexprep(val,'"','\\"');
        val=regexprep(val,'^"','\\"');
    end
    if(len==1)
        obj=['"' checkname(name,varargin{:}) '": ' '"',val,'"'];
	if(isempty(name)) obj=['"',val,'"']; end
        txt=sprintf('%s%s%s%s',txt,repmat(sprintf('\t'),1,level),obj);
    else
        txt=sprintf('%s%s%s%s',txt,repmat(sprintf('\t'),1,level+1),['"',val,'"']);
    end
    if(e==len) sep=''; end
    txt=sprintf('%s%s',txt,sep);
end
if(len>1) txt=sprintf('%s\n%s%s',txt,padding1,']'); end

%%-------------------------------------------------------------------------
function txt=mat2json(name,item,level,varargin)
if(~isnumeric(item) && ~islogical(item))
        error('input is not an array');
end

padding1=repmat(sprintf('\t'),1,level);
padding0=repmat(sprintf('\t'),1,level+1);

if(length(size(item))>2 || issparse(item) || ~isreal(item) || ...
   isempty(item) ||jsonopt('ArrayToStruct',0,varargin{:}))
    if(isempty(name))
    	txt=sprintf('%s{\n%s"_ArrayType_": "%s",\n%s"_ArraySize_": %s,\n',...
              padding1,padding0,class(item),padding0,regexprep(mat2str(size(item)),'\s+',',') );
    else
    	txt=sprintf('%s"%s": {\n%s"_ArrayType_": "%s",\n%s"_ArraySize_": %s,\n',...
              padding1,checkname(name,varargin{:}),padding0,class(item),padding0,regexprep(mat2str(size(item)),'\s+',',') );
    end
else
    if(isempty(name))
    	txt=sprintf('%s%s',padding1,matdata2json(item,level+1,varargin{:}));
    else
        if(numel(item)==1 && jsonopt('NoRowBracket',1,varargin{:})==1)
            numtxt=regexprep(regexprep(matdata2json(item,level+1,varargin{:}),'^\[',''),']','');
           	txt=sprintf('%s"%s": %s',padding1,checkname(name,varargin{:}),numtxt);
        else
    	    txt=sprintf('%s"%s": %s',padding1,checkname(name,varargin{:}),matdata2json(item,level+1,varargin{:}));
        end
    end
    return;
end
dataformat='%s%s%s%s%s';

if(issparse(item))
    [ix,iy]=find(item);
    data=full(item(find(item)));
    if(~isreal(item))
       data=[real(data(:)),imag(data(:))];
       if(size(item,1)==1)
           % Kludge to have data's 'transposedness' match item's.
           % (Necessary for complex row vector handling below.)
           data=data';
       end
       txt=sprintf(dataformat,txt,padding0,'"_ArrayIsComplex_": ','1', sprintf(',\n'));
    end
    txt=sprintf(dataformat,txt,padding0,'"_ArrayIsSparse_": ','1', sprintf(',\n'));
    if(size(item,1)==1)
        % Row vector, store only column indices.
        txt=sprintf(dataformat,txt,padding0,'"_ArrayData_": ',...
           matdata2json([iy(:),data'],level+2,varargin{:}), sprintf('\n'));
    elseif(size(item,2)==1)
        % Column vector, store only row indices.
        txt=sprintf(dataformat,txt,padding0,'"_ArrayData_": ',...
           matdata2json([ix,data],level+2,varargin{:}), sprintf('\n'));
    else
        % General case, store row and column indices.
        txt=sprintf(dataformat,txt,padding0,'"_ArrayData_": ',...
           matdata2json([ix,iy,data],level+2,varargin{:}), sprintf('\n'));
    end
else
    if(isreal(item))
        txt=sprintf(dataformat,txt,padding0,'"_ArrayData_": ',...
            matdata2json(item(:)',level+2,varargin{:}), sprintf('\n'));
    else
        txt=sprintf(dataformat,txt,padding0,'"_ArrayIsComplex_": ','1', sprintf(',\n'));
        txt=sprintf(dataformat,txt,padding0,'"_ArrayData_": ',...
            matdata2json([real(item(:)) imag(item(:))],level+2,varargin{:}), sprintf('\n'));
    end
end
txt=sprintf('%s%s%s',txt,padding1,'}');

%%-------------------------------------------------------------------------
function txt=matdata2json(mat,level,varargin)
if(size(mat,1)==1)
    pre='';
    post='';
    level=level-1;
else
    pre=sprintf('[\n');
    post=sprintf('\n%s]',repmat(sprintf('\t'),1,level-1));
end
if(isempty(mat))
    txt='null';
    return;
end
floatformat=jsonopt('FloatFormat','%.10g',varargin{:});
formatstr=['[' repmat([floatformat ','],1,size(mat,2)-1) [floatformat sprintf('],\n')]];

if(nargin>=2 && size(mat,1)>1 && jsonopt('ArrayIndent',1,varargin{:})==1)
    formatstr=[repmat(sprintf('\t'),1,level) formatstr];
end
txt=sprintf(formatstr,mat');
txt(end-1:end)=[];
if(islogical(mat) && jsonopt('ParseLogical',0,varargin{:})==1)
   txt=regexprep(txt,'1','true');
   txt=regexprep(txt,'0','false');
end
%txt=regexprep(mat2str(mat),'\s+',',');
%txt=regexprep(txt,';',sprintf('],\n['));
% if(nargin>=2 && size(mat,1)>1)
%     txt=regexprep(txt,'\[',[repmat(sprintf('\t'),1,level) '[']);
% end
txt=[pre txt post];
if(any(isinf(mat(:))))
    txt=regexprep(txt,'([-+]*)Inf',jsonopt('Inf','"$1_Inf_"',varargin{:}));
end
if(any(isnan(mat(:))))
    txt=regexprep(txt,'NaN',jsonopt('NaN','"_NaN_"',varargin{:}));
end

%%-------------------------------------------------------------------------
function newname=checkname(name,varargin)
isunpack=jsonopt('UnpackHex',1,varargin{:});
newname=name;
if(isempty(regexp(name,'0x([0-9a-fA-F]+)_','once')))
    return
end
if(isunpack)
    isoct=jsonopt('IsOctave',0,varargin{:});
    if(~isoct)
        newname=regexprep(name,'(^x|_){1}0x([0-9a-fA-F]+)_','${native2unicode(hex2dec($2))}');
    else
        pos=regexp(name,'(^x|_){1}0x([0-9a-fA-F]+)_','start');
        pend=regexp(name,'(^x|_){1}0x([0-9a-fA-F]+)_','end');
        if(isempty(pos)) return; end
        str0=name;
        pos0=[0 pend(:)' length(name)];
        newname='';
        for i=1:length(pos)
            newname=[newname str0(pos0(i)+1:pos(i)-1) char(hex2dec(str0(pos(i)+3:pend(i)-1)))];
        end
        if(pos(end)~=length(name))
            newname=[newname str0(pos0(end-1)+1:pos0(end))];
        end
    end
end

author: plotly


text:
filename='Cantilever_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.4;
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;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function K = dexpKernDiagCompute(kern, x)

% DEXPKERNDIAGCOMPUTE Compute diagonal of the double exponential kernel.
%
% FORMAT
% DESC computes the diagonal of the kernel matrix for the double
% exponential kernel given a column vector of inputs.
% ARG kern : the kernel structure for which the kernel matrix is computed.
% ARG x : input data in the form of a design matrix.
% RETURN K : a vector of the same size as x containing the diagonal of the
% kernel matrix computed at the given points.
%
% SEEALSO : dexpKernParamInit, kernDiagCompute, kernCreate, dexpKernCompute
%
% COPYRIGHT : David Luengo, 2009

% KERN

K = 0.5 * kern.variance * kern.decay * ones(size(x));