text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the CHANNEL_CHANNEL-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Channel_Channel()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  128;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  128;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
ds= min(Lx/(2*Nx),Ly/(2*Ny));  % Lagrangian spacing
L = 0.9*Lx;                    % Length of Channel
w = 0.2*Ly;                    % Width of Channel
struct_name = 'channel'; % Name for .vertex, .spring, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,L,w,Lx,Ly);

% Plot Geometry to test
plot(xLag(1:end/2),yLag(1:end/2),'r-'); hold on;
plot(xLag(end/2+1:end),yLag(end/2+1:end),'r-'); hold on;

plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis square;

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
%k_Spring = 1e7;
%print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name);

% Prints .beam file!
%k_Beam = 0.5; C = 0.0;
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1e7;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,L,w,Lx,Ly)

% The immsersed structure is a channel %
x = (Lx-L)/2:ds:(L+(Lx-L)/2);  %xPts
yBot = (Ly-w)/2;               %yVal for bottom of Channel
yTop = Ly - (Ly-w)/2;          %yVal for top of Channel

xLag = [x x];
yLag = [yBot*ones(1,length(x)) yTop*ones(1,length(x))];

author: nickabattista


text:
function [fevd_record,fevd_estimates]=panel2fevd(struct_irf_record,gamma_record,It,Bu,IRFperiods,n,FEVDband)

% because fevd is obtained from IRFs, and because IRFs are common to all units, fevd can be run only once
[fevd_record]=bear.fevd(struct_irf_record,gamma_record,It,Bu,n,IRFperiods,FEVDband);
% compute posterior estimates
[fevd_estimates]=bear.fevdestimates(fevd_record,n,IRFperiods,FEVDband);

author: european-central-bank

repo: BEAR-toolbox


text:
function op = prox_Sl1( lambda )
%PROX_SL1    Sorted/Ordered L1 norm.
%    OP = PROX_L1( lambda ) implements the nonsmooth function
%        OP(X) = sum(lambda.*sort(abs(X),'descend'))
%    where lambda is strictly positive and sorted in decreasing order,
%    in which case this function is a norm (and hence convex).
%    If lambda is a scalar, it will be expanded to all elements, but in this
%    case it is equivalent to the (scaled) usual l1 norm.
%
%    Notice: this function uses mex files. Some pre-compiled binaries
%       for common systems are included; if these do not work for you,
%       then please install yourself. In the mexFiles/ subdirectory,
%       run the file "makeMex.m"
% Reference:
%  http://www-stat.stanford.edu/~candes/OrderedL1/
% "Statistical Estimation and Testing via the Ordered L1 Norm"
% by M. Bogdan, E. van den Berg, W. Su, and E. J. Candès
% 2013
%
% See also prox_l1.m, makeMex.m

if nargin == 0,
	lambda = 1;
elseif ~isnumeric( lambda ) || ~isreal( lambda ) ||  any( lambda(:) <= 0 )
	error( 'Argument lambda must have all positive entries.' );
end
if ~issorted(flipud(lambda(:)))
    error( 'Argument lambda must be sorted in decreasing order.');
end
if numel(lambda)==1
    warning('TFOCS:prox_SL1','When lambda is a scalar, we recommend prox_l1.m instead pf prox_SL1.m');
end

% The mex file is in the child directory mexFiles/
% Check for its existence. First, add the right paths
addpath(fullfile(tfocs_where,'mexFiles'));
if 3 ~= exist('proxAdaptiveL1Mex','file')
    makeMex;
    % check that it worked
    if 3 ~= exist('proxAdaptiveL1Mex','file')
        disp('Compilation of mex files for prox_SL1.m failed; please report this error');
    end
end

f       = @(x) sum( lambda(:) .* sort(abs(x(:)), 'descend') );
prox_f  = @(x,t) proxOrderedL1(x,t.*lambda);

op      = tfocs_prox( f, prox_f , 'vector' ); % Allow vector stepsizes

end

% -- subroutines --
function x = proxOrderedL1(y,lambda)
    % Normalization
    lambda = lambda(:);
    y      = y(:);
    sgn    = sign(y);
    [y,idx] = sort(abs(y),'descend');

    % Simplify the problem
    k = find(y > lambda,1,'last');

    % Compute solution and re-normalize
    n = numel(y);
    x = zeros(n,1);

    if (~isempty(k))
        v1 = y(1:k);
        if numel(lambda) > 1
            v2 = lambda(1:k);
        else
            v2 = lambda*ones(k,1); % if lambda is a scalar, implicity make it lambda*ones(size(y))
        end
        v = proxAdaptiveL1Mex(v1,v2);
        x(idx(1:k)) = v;
    end

    % Restore signs
    x = sgn .* x;
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% Separating ellipsoids in 2D
% Joelle Skaf - 11/06/05
% (a figure is generated)
%
% Finds a separating hyperplane between 2 ellipsoids {x| ||Ax+b||^2<=1} and
% {y | ||Cy + d||^2 <=1} by solving the following problem and using its
% dual variables:
%               minimize    ||w||
%                   s.t.    ||Ax + b||^2 <= 1       : lambda
%                           ||Cy + d||^2 <= 1       : mu
%                           x - y == w              : z
% the vector z will define a separating hyperplane because z'*(x-y)>0

% input data
n = 2;
A = eye(n);
b = zeros(n,1);
C = [2 1; -.5 1];
d = [-3; -3];

% solving for the minimum distance between the 2 ellipsoids and finding
% the dual variables
cvx_begin
    variables x(n) y(n) w(n)
    dual variables lam muu z
    minimize ( norm(w,2) )
    subject to
    lam:    square_pos( norm (A*x + b) ) <= 1;
    muu:    square_pos( norm (C*y + d) ) <= 1;
    z:      x - y == w;
cvx_end

t = (x + y)/2;
p=z;
p(1) = z(2); p(2) = -z(1);
c = linspace(-2,2,100);
q = repmat(t,1,length(c)) +p*c;

% figure
nopts = 1000;
angles = linspace(0,2*pi,nopts);
[u,v] = meshgrid([-2:0.01:4]);
z1 = (A(1,1)*u + A(1,2)*v + b(1)).^2 + (A(2,1)*u + A(2,2)*v + b(2)).^2;
z2 = (C(1,1)*u + C(1,2)*v + d(1)).^2 + (C(2,1)*u + C(2,2)*v + d(2)).^2;
contour(u,v,z1,[1 1]);
hold on;
contour(u,v,z2,[1 1]);
axis square
plot(x(1),x(2),'r+');
plot(y(1),y(2),'b+');
line([x(1) y(1)],[x(2) y(2)]);
plot(q(1,:),q(2,:),'k');

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function fem1d_spectral_numeric_test ( )

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_SPECTRAL_NUMERIC_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FEM1D_SPECTRAL_NUMERIC library.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Run FEM1D_SPECTRAL_NUMERIC with an increasing number of basis functions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N        L2-error        H1-error\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 2 : 11
    [ l2, h1 ] = fem1d_spectral_numeric ( n );
    fprintf ( 1, '  %2d  %14.6g  %14.6g\n', n, l2, h1 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_SPECTRAL_NUMERIC_TEST:\n' );
  fprintf ( 1, '  Normal end of execution\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f=burger(u,deriv)
if nargin<2,
	deriv=0;
end;
if deriv==0,
	f=0.5*u.^2;
else
	f=u;
end;

repo: Aero-matlab


text:
function printPostSolver(fid, expInfo)
%PRINTPOSTSOLVER   Print commands after solution has been found
%
% Calling sequence:
%   PRINTPOSTSOLVER(FID, EXPINFO)
% where
%   FID:        ID of a file-writing stream.
%   EXPINFO:    Struct containing information for printing the problem.

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

% Extract info from the EXPINFO struct:
allVarNames = expInfo.allVarNames;
allVarString = expInfo.allVarString;
numVars = expInfo.numVars;
indVarNameSpace = expInfo.indVarNameSpace;

% Print commands that will create a plot of the solution obtained:
fprintf(fid, '\n%%%% Plot the solution.\n');
fprintf(fid, 'figure\n');
if ( numVars == 1 ) % Scalar case
    fprintf(fid, 'plot(%s, ''LineWidth'', 2)\n', allVarString);
else
    fprintf(fid, 'plot([%s], ''LineWidth'', 2)\n', allVarString);
end

% Title and xlabel
fprintf(fid, 'title(''Final solution''), xlabel(''%s'')', indVarNameSpace);

% Deal with ylabel (scalar problem) or legend (systems):
if ( numVars == 1 )
    % Scalar problem:
    fprintf(fid, ', ylabel(''%s'')', allVarNames{:});
else
    % Coupled system. Create a legend.
    leg = '';
    for k = 1:numel(allVarNames)-1
        leg = [leg '''' allVarNames{k} '''' ','];
    end
    leg = [leg '''' allVarNames{k+1} ''''];
    fprintf(fid, ', legend(%s)\n', leg);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
classdef Tz < ETS3

    methods
        function obj = Tz(varargin)
            obj@ETS3('Tz', varargin{:});

        end

        % function s = char(ets)
        %     s = ets.what;
        % end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [R, scale]=arqr(v, p, mcor)
%ARQR	QR factorization for least squares estimation of AR model.
%
%  [R, SCALE]=ARQR(v,p,mcor) computes the QR factorization needed in
%  the least squares estimation of parameters of an AR(p) model. If
%  the input flag mcor equals one, a vector of intercept terms is
%  being fitted. If mcor equals zero, the process v is assumed to have
%  mean zero. The output argument R is the upper triangular matrix
%  appearing in the QR factorization of the AR model, and SCALE is a
%  vector of scaling factors used to regularize the QR factorization.
%
%  ARQR is called by ARFIT.
%
%  See also ARFIT.

%  Modified 29-Dec-99
%  Author: Tapio Schneider
%          tapio@gps.caltech.edu

  % n: number of time steps; m: dimension of state vectors
  [n,m] = size(v);

  ne    = n-p;                  % number of block equations of size m
  np    = m*p+mcor;             % number of parameter vectors of size m

  % If the intercept vector w is to be fitted, least squares (LS)
  % estimation proceeds by solving the normal equations for the linear
  % regression model
  %
  %                  v(k,:)' = Aaug*u(k,:)' + noise(C)
  %
  % with Aaug=[w A] and `predictors'
  %
  %              u(k,:) = [1 v(k-1,:) ...  v(k-p,:)].
  %
  % If the process mean is taken to be zero, the augmented coefficient
  % matrix is Aaug=A, and the regression model
  %
  %                u(k,:) = [v(k-1,:) ...  v(k-p,:)]
  %
  % is fitted.
  % The number np is the dimension of the `predictors' u(k).

  % Assemble the data matrix K (of which a QR factorization will be computed)
  K = zeros(ne,np+m);                 % initialize K
  if (mcor == 1)
    % first column of K consists of ones for estimation of intercept vector w
    K(:,1) = ones(ne,1);
  end

  % Assemble `predictors' u in K
  for j=1:p
    K(:, mcor+m*(j-1)+1:mcor+m*j) = [v(p-j+1:n-j, :)];
  end
  % Add `observations' v (left hand side of regression model) to K
  K(:,np+1:np+m) = [v(p+1:n, :)];

  % Compute regularized QR factorization of K: The regularization
  % parameter delta is chosen according to Higham's (1996) Theorem
  % 10.7 on the stability of a Cholesky factorization. Replace the
  % regularization parameter delta below by a parameter that depends
  % on the observational error if the observational error dominates
  % the rounding error (cf. Neumaier, A. and T. Schneider, 2001:
  % "Estimation of parameters and eigenmodes of multivariate
  % autoregressive models", ACM Trans. Math. Softw., 27, 27--57.).
  q     = np + m;             % number of columns of K
  delta = (q^2 + q + 1)*eps;  % Higham's choice for a Cholesky factorization
  scale = sqrt(delta)*sqrt(sum(K.^2));
  R     = triu(qr([K; diag(scale)]));

author: Sable

repo: mcbench-benchmarks


text:
classdef PTKFissurenessHessianFactor < PTKPlugin
    % PTKFissureApproximation. Plugin to detect fissures using analysis of the
    %     Hessian matrix
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     This is an intermediate stage towards lobar segmentation.
    %
    %     PTKFissurenessHessianFactor computes the components of the fissureness
    %     generated using analysis of eigenvalues of the Hessian matrix.
    %
    %     For more information, see
    %     [Doel et al., Pulmonary lobe segmentation from CT images using
    %     fissureness, airways, vessels and multilevel B-splines, 2012]
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Fissureness <BR>(Hessian part)'
        ToolTip = 'The part of the fissureness filter which uses Hessian-based analysis'
        Category = 'Fissures'

        AllowResultsToBeCached = false
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
        Version = 2

        MemoryCachePolicy = 'Temporary'
        DiskCachePolicy = 'Off'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            reporting.UpdateProgressValue(0);
            left_and_right_lungs = dataset.GetResult('PTKLeftAndRightLungs');

            right_lung = dataset.GetResult('PTKGetRightLungROI');

            fissureness_right = PTKFissurenessHessianFactor.ComputeFissureness(right_lung, left_and_right_lungs, reporting, false);

            reporting.UpdateProgressValue(50);
            left_lung = dataset.GetResult('PTKGetLeftLungROI');
            fissureness_left = PTKFissurenessHessianFactor.ComputeFissureness(left_lung, left_and_right_lungs, reporting, true);

            reporting.UpdateProgressValue(100);
            results = PTKCombineLeftAndRightImages(dataset.GetTemplateImage(PTKContext.LungROI), fissureness_left, fissureness_right, left_and_right_lungs);

            results.ImageType = PTKImageType.Scaled;
        end
    end

    methods (Static, Access = private)

        function lung = DuplicateImageInMask(lung, mask_raw)
            mask_raw = mask_raw > 0;
            if any(mask_raw(:) > 0)
                [~, labelmatrix] = bwdist(mask_raw);
                lung(~mask_raw(:)) = lung(labelmatrix(~mask_raw(:)));
            else
                lung(:) = 0;
            end
        end

        function fissureness = ComputeFissureness(image_data, left_and_right_lungs, reporting, is_left_lung)

            left_and_right_lungs = left_and_right_lungs.Copy;
            left_and_right_lungs.ResizeToMatch(image_data);
            image_data.ChangeRawImage(PTKFissurenessHessianFactor.DuplicateImageInMask(image_data.RawImage, left_and_right_lungs.RawImage));

            mask = [];
            fissureness = PTKImageDividerHessian(image_data, @PTKFissurenessHessianFactor.ComputeFissurenessPartImage, mask, 1.5, [], false, false, is_left_lung, reporting);
        end

        function fissureness_wrapper = ComputeFissurenessPartImage(hessian_eigs_wrapper, voxel_size)
            fissureness_wrapper = PTKComputeFissurenessFromHessianeigenvalues(hessian_eigs_wrapper, voxel_size);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [hPa] = mmHg2hPa(mmHg)
% Convert pressure from millimeters of mercury at 0 degrees C to
% hectopascals
% Chad Greene 2012
hPa = mmHg*1.33322;


text:
function sae = saesetup(size)
    for u = 2 : numel(size)
        sae.ae{u-1} = nnsetup([size(u-1) size(u) size(u-1)]);
    end
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function [ fx ] = f_BSLinGame(x,theta,u,in)
% evolution function for BSL learner engaging in dyadic games
% [ fx ] = f_BSLinGame(x,theta,u,in)
% BSL is simply tracking the log-odds of P(o_t=1|o_{t-1}), where o is the
% opponent's binary action. This variable is updated according to a
% Laplace-Kalman filter, yielding 2 sufficient statistics (m and V) per
% combination of past outcome. BSL can learn sequences of arbitrary depth
% (K). For example, if K=1, then BSL tracks 2 probabilities, namely:
% P(o_t=1|o_{t-1}=1) and P(o_t=1|o_{t-1}=0). More generally, BSL tracks 2^K
% probabilities. In this scheme, the only evolution param (theta) is BSL's
% prior volatity about the log-odds.
% Note: unsampled sequences will eventually be "forgotten", since the
% prediction step in the Laplace-Kalman update will dilute any previously
% sampled evidence.
% IN:
%   - x: sufficient statistics of log-odds of P(o=1):
%       x(1:2^K)= E[log-odds]
%       x((2^K)+1:2^(K+1))= log V[log-odds]
%   - theta: BSL's prior volatity
%   - u: sequence of past outcomes:
%       u(1)= opponent's last move
%       u(2)= learner's last move
%       u(3:K+2) = sequence of K past opponent's moves
%   - in: depth of sequence learning
% OUT:
%   - fx: updated sufficient statistics of log-odds of P(o=1)

if VBA_isWeird (u) % e.g., 1st trial
    fx = x;
    return
end

u(2) = []; % remove agent's previous move
[fx] = f_BSL(x,theta,u,in);

author: MBB-team

repo: VBA-toolbox


text:
function net = gtminit(net, options, data, samp_type, varargin)
%GTMINIT Initialise the weights and latent sample in a GTM.
%
%	Description
%	NET = GTMINIT(NET, OPTIONS, DATA, SAMPTYPE) takes a GTM NET and
%	generates a sample of latent data points and sets the centres (and
%	widths if appropriate) of NET.RBFNET.
%
%	If the SAMPTYPE is 'REGULAR', then regular grids of latent data
%	points and RBF centres are created.  The dimension of the latent data
%	space must be 1 or 2.  For one-dimensional latent space, the
%	LSAMPSIZE parameter gives the number of latent points and the
%	RBFSAMPSIZE parameter gives the number of RBF centres.  For a two-
%	dimensional latent space, these parameters must be vectors of length
%	2 with the number of points in each of the x and y directions to
%	create a rectangular grid.  The widths of the RBF basis functions are
%	set by a call to RBFSETFW passing OPTIONS(7) as the scaling
%	parameter.
%
%	If the SAMPTYPE is 'UNIFORM' or 'GAUSSIAN' then the latent data is
%	found by sampling from a uniform or Gaussian distribution
%	correspondingly.  The RBF basis function parameters are set by a call
%	to RBFSETBF with the DATA parameter as dataset and the OPTIONS
%	vector.
%
%	Finally, the output layer weights of the RBF are initialised by
%	mapping the mean of the latent variable to the mean of the target
%	variable, and the L-dimensional latent variale variance to the
%	variance of the targets along the first L principal components.
%
%	See also
%	GTMINIT2, GTM, GTMEM, PCA, RBFSETBF, RBFSETFW
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

% Check type of sample
stypes = {'regular', 'uniform', 'gaussian'};
if (strcmp(samp_type, stypes)) == 0
  error('Undefined sample type.')
end

if net.dim_latent > size(data, 2)
  error('Latent space dimension must not be greater than data dimension')
end
nlatent = net.gmmnet.ncentres;
nhidden = net.rbfnet.nhidden;

% Create latent data sample and set RBF centres

switch samp_type
case 'regular'
   if nargin ~= 6
      error('Regular type must specify latent and RBF shapes');
   end
   l_samp_size = varargin{1};
   rbf_samp_size = varargin{2};
   if round(l_samp_size) ~= l_samp_size
      error('Latent sample specification must contain integers')
   end
   % Check existence and size of rbf specification
   if any(size(rbf_samp_size) ~= [1 net.dim_latent]) | ...
         prod(rbf_samp_size) ~= nhidden
      error('Incorrect specification of RBF centres')
   end
   % Check dimension and type of latent data specification
   if any(size(l_samp_size) ~= [1 net.dim_latent]) | ...
         prod(l_samp_size) ~= nlatent
      error('Incorrect dimension of latent sample spec.')
   end
   if net.dim_latent == 1
      net.X = [-1:2/(l_samp_size-1):1]';
      net.rbfnet.c = [-1:2/(rbf_samp_size-1):1]';
      net.rbfnet = rbfsetfw(net.rbfnet, options(7));
   elseif net.dim_latent == 2
      net.X = gtm_rctg(l_samp_size);
      net.rbfnet.c = gtm_rctg(rbf_samp_size);
      net.rbfnet = rbfsetfw(net.rbfnet, options(7));
   else
      error('For regular sample, input dimension must be 1 or 2.')
   end

case {'uniform', 'gaussian'}
   if strcmp(samp_type, 'uniform')
      net.X = 2 * (rand(nlatent, net.dim_latent) - 0.5);
   else
      % Sample from N(0, 0.25) distribution to ensure most latent
      % data is inside square
      net.X = randn(nlatent, net.dim_latent)/2;
   end
   net.rbfnet = rbfsetbf(net.rbfnet, options, net.X);
otherwise
   % Shouldn't get here
   error('Invalid sample type');

end

% Latent data sample and basis function parameters chosen.
% Now set output weights
[PCcoeff, PCvec] = pca(data);

% Scale PCs by eigenvalues
A = PCvec(:, 1:net.dim_latent)*diag(sqrt(PCcoeff(1:net.dim_latent)));

[temp, Phi] = rbffwd(net.rbfnet, net.X);
% Normalise X to ensure 1:1 mapping of variances and calculate weights
% as solution of Phi*W = normX*A'
normX = (net.X - ones(size(net.X))*diag(mean(net.X)))*diag(1./std(net.X));
net.rbfnet.w2 = Phi \ (normX*A');
% Bias is mean of target data
net.rbfnet.b2 = mean(data);

% Must also set initial value of variance
% Find average distance between nearest centres
% Ensure that distance of centre to itself is excluded by setting diagonal
% entries to realmax
net.gmmnet.centres = rbffwd(net.rbfnet, net.X);
d = dist2(net.gmmnet.centres, net.gmmnet.centres) + ...
  diag(ones(net.gmmnet.ncentres, 1)*realmax);
sigma = mean(min(d))/2;

% Now set covariance to minimum of this and next largest eigenvalue
if net.dim_latent < size(data, 2)
  sigma = min(sigma, PCcoeff(net.dim_latent+1));
end
net.gmmnet.covars = sigma*ones(1, net.gmmnet.ncentres);

% Sub-function to create the sample data in 2d
function sample = gtm_rctg(samp_size)

xDim = samp_size(1);
yDim = samp_size(2);
% Produce a grid with the right number of rows and columns
[X, Y] = meshgrid([0:1:(xDim-1)], [(yDim-1):-1:0]);

% Change grid representation
sample = [X(:), Y(:)];

% Shift grid to correct position and scale it
maxXY= max(sample);
sample(:,1) = 2*(sample(:,1) - maxXY(1)/2)./maxXY(1);
sample(:,2) = 2*(sample(:,2) - maxXY(2)/2)./maxXY(2);
return;

author: ilarinieminen

repo: SOM-Toolbox


text: %GT Datafile overload

author: marianux

repo: ecg-kit

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


text:
function transfo = createBasisTransform(source, target)
%CREATEBASISTRANSFORM Compute matrix for transforming a basis into another basis
%
%   TRANSFO = createBasisTransform(SOURCE, TARGET)
%   Both SOURCE and TARGET represent basis, in the following form:
%   [x0 y0  ex1 ey1  ex2 ey2]
%   [y0 y0] is the origin of the basis, [ex1 ey1] is the first direction
%   vector, and [ex2 ey2] is the second direction vector.
%
%   The result TRANSFO is a 3-by-3 matrix such that a point expressed with
%   coordinates of the first basis will be represented by new coordinates
%   P2 = transformPoint(P1, TRANSFO) in the target basis.
%
%   TRANSFO = createBasisTransform(TARGET)
%   Assumes the source is the standard (Oij) basis, with origin at (0,0),
%   first direction vector equal to (1,0) and second direction  vector
%   equal to (0,1).
%
%
%   Example
%     % define source and target bases
%     src = [ 0 0   1  0    0  1];
%     tgt = [20 0  .5 .5  -.5 .5];
%     trans = createBasisTransform(src, tgt);
%     % create a polygon in source basis
%     poly = [10 10;30 10; 30 20; 20 20;20 40; 10 40];
%     figure;
%     subplot(121); drawPolygon(poly, 'b'); axis equal; axis([-10 50 -10 50]);
%     hold on; drawLine([0 0 1 0], 'k'); drawLine([0 0 0 1], 'k');
%     drawLine([20 0 1 1], 'r'); drawLine([20 0 -1 1], 'r');
%     t = -1:5; plot(t*5+20, t*5, 'r.'); plot(-t*5+20, t*5, 'r.');
%     % transform the polygon in target basis
%     poly2 = transformPoint(poly, trans);
%     subplot(122); drawPolygon(poly2, 'b'); axis equal; axis([-10 50 -10 50]);
%     hold on; drawLine([0 0 1 0], 'r'); drawLine([0 0 0 1], 'r');
%     t = -1:5; plot(t*10, zeros(size(t)), 'r.'); plot(zeros(size(t)), t*10, 'r.');
%
%   See also
%   transforms2d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-12-03,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% init basis transform to identity
t1 = eye(3);
t2 = eye(3);

if nargin == 2
    % from source to reference basis
    t1(1:2, 1) = source(3:4);
    t1(1:2, 2) = source(5:6);
    t1(1:2, 3) = source(1:2);
else
    % if only one input, use first input as target basis, and leave the
    % first matrix to identity
    target = source;
end

% from reference to target basis
t2(1:2, 1) = target(3:4);
t2(1:2, 2) = target(5:6);
t2(1:2, 3) = target(1:2);

% compute transform matrix
transfo = zeros(3, 3);
maxSz = 1;
for i = 1:maxSz
    % coordinate of three reference points in source basis
    po = t1(1:2, 3, i)';
    px = po + t1(1:2, 1, i)';
    py = po + t1(1:2, 2, i)';

    % express coordinates of reference points in the new basis
    t2i = inv(t2(:,:,i));
    pot = transformPoint(po, t2i);
    pxt = transformPoint(px, t2i);
    pyt = transformPoint(py, t2i);

    % compute direction vectors in new basis
    vx = pxt - pot;
    vy = pyt - pot;

    % concatenate result in a 3-by-3 affine transform matrix
    transfo(:,:,i) = [vx' vy' pot' ; 0 0 1];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function ind = rhoInside(rho,minRho,maxRho)

minr = mod(minRho+1e-6,2*pi)-3e-6; % in [-2e-6,2*pi-3e-6]
maxr = mod(maxRho-1e-6,2*pi)+3e-6; % in [2e-6,2*pi+3e-6]

if minr < 0
  rho = mod(rho+1e-6,2*pi);
else
  rho = mod(rho-1e-6,2*pi);
end

if minr < maxr
  ind = rho > minr & rho < maxr;
else
  ind = rho > minr | rho < maxr;
end

author: mtex-toolbox

repo: mtex

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


text:
function [F] = Filter1D(N,Nc,s)

% function [F] = Filter1D(N,Nc,s)
% Purpose : Initialize 1D filter matrix of size N.
%           Order of exponential filter is (even) s with cutoff at Nc;

Globals1D;
filterdiag = ones(N+1,1);
alpha = -log(eps);

% Initialize filter function
for i=Nc:N
    filterdiag(i+1) = exp(-alpha*((i-Nc)/(N-Nc))^s);
end;

F = V*diag(filterdiag)*invV;
return;


text:
function value = p08_autonomous ( )

%*****************************************************************************80
%
%% P08_AUTONOMOUS reports whether problem P08 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %

    h = X * theta;
    errors = h - y;
    delta = X' * errors;
    theta = theta - (alpha / m) * delta;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCost(X, y, theta);

end

end

author: zsiciarz

repo: ml-coursera


text:
function res = stackSlice(img, dir, slice)
%STACKSLICE Extract a planar slice from a 3D image
%
%   SLICE = stackSlice(IMG, DIR, INDEX)
%   IMG is either a 3D or a 4D (3D+color) image.
%   DIR is 1, 2 or 3, or 'x', 'y', 'z'
%   INDEX is the slice index, between 1 and the number of voxels in the DIR
%   direction (using physical indexing).
%
%   If the input image has size NY-by-NX-by-NZ (gray scale image), or
%   NY-by-NX-by-3-by-NZ (color image), then the resulting slice has size:
%   - NZ-by-NY(-by-3) if DIR is 1 or DIR is 'x'
%   - NX-by-NZ(-by-3) if DIR is 2 or DIR is 'y'
%   - NY-by-NX(-by-3) if DIR is 3 or DIR is 'z'
%
%   The functions returns a planar image with the same type as the
%   original.
%
%
%   Example
%   % Display 3 slices of a MRI head
%     img = analyze75read(analyze75info('brainMRI.hdr'));
%     figure(1); clf; hold on;
%     for i=1:3
%         subplot(3, 1, i);
%         imshow(stackSlice(img, 'x', 20+20*i)');
%         set(gca, 'ydir', 'normal')
%     end
%
%   See also
%     imStacks, imMiddleSlice

%   ------
%   author: David Legland, david.legland@inra.fr
%   INRA - Cepia Software Platform
%   Created: 2007-08-14,    using Matlab 7.4.0.287 (R2007a)
%   http://www.pfl-cepia.inra.fr/index.php?page=slicer

%   HISTORY
%   2010-12-02 code cleanup, fix dimension order of output slices
%   2011-04-26 use xyz convention for slice index

% image size
dim = size(img);

% convert to an index between 1 and 3, in xyz order
dir = parseAxisIndex(dir);

% check dimension of input array
nd = length(dim);
if nd < 3 || nd > 4
    error('Input array should have dimension 3 or 4');
end

% Extract either a grayscale or a color slice
if nd == 3
    % gray-scale image
    switch dir
        case 1
            % X-slice, result is a YZ slice with size NZ-by-NY
            res = squeeze(img(:, slice, :))';
        case 2
            % Y-slice, result is a ZX slice with size NX-by-NZ
            res = squeeze(img(slice, :, :));
        case 3
            % Z-slice, result is a XY slice with size NY-by-NX
            res = img(:, :, slice);
        otherwise
            error('Direction should be comprised between 1 and 3');
    end

else
    % color images: keep channel dim as third dimension, and put slice
    % dimension as last dimension.
    switch dir
        case 1
            % X-slice, result is a YZ color slice with size NZ-by-NY-by-3
            res = squeeze(permute(img(:, slice, :, :), [4 1 3 2]));
        case 2
            % Y-slice, result is a ZX color slice with size NX-by-NZ-by-3
            res = squeeze(permute(img(slice, :, :,:), [2 4 3 1]));
        case 3
            % Z-slice, result is a XY color slice with size NY-by-NX-by-3
            res = img(:, :, :, slice);
        otherwise
            error('Direction should be comprised between 1 and 3');
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function dG = filterDoog( dims, sigmas, nderivs, show )
% n-dim difference of offset Gaussian DooG filter (Gaussian derivative).
%
% Creates a nd derivative of Gaussian kernel.  For all but d==1 use
% primarily for visualization purposes -- for filtering better to use the
% indvidiual seperable kernels for efficiency purposes.
%
% USAGE
%  dG = filterDoog( dims, sigmas, nderivs, [show] )
%
% INPUTS
%  dims    - nd element vector of dimensions of final Gaussian
%  sigmas  - sigmas for nd Gaussian
%  nderivs - order of derivative along each of the nd dimensions
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%  dG      - The derivative of Gaussian filter
%
% EXAMPLE
%  dG1 = filterDoog( 43, 2, 3, 1 ); %1D
%  dG2 = filterDoog( [41 41], [3 3], [1,1], 2 ); %2D
%  dG3 = filterDoog( [101 101 101], [4,4,10], [1,1,0], 3 ); %3D
%
% See also FILTERGAUSS, NORMPDF2, FILTERDOG2D, FILTERGABOR2D
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

nd = length( dims );
if( nargin<4 || isempty(show) ); show=0; end;
if( length(sigmas)~=nd ); error('invalid sigmas'); end
if( length(nderivs)~=nd ); error('invalid nderivs'); end

% get initial Gaussian
dG = filterGauss( dims, [], sigmas.^2, 0 );

% compute derivatives along each axis
for d=1:nd
  if( d==1 );
    dOp=.5*[-1 0 1]';
  else
    dOp = .5*permute( [-1 0 1]', d:-1:1 );
  end
  if( nd==1 || nd==2 )
    for i=1:nderivs(d); dG = conv2( dG, dOp, 'same' ); end
  else
    for i=1:nderivs(d); dG = convn( dG, dOp, 'same' ); end
  end
end

% normalize (don't need to adjust mean since DOOG always have 0 mean)
dG=dG/norm(dG(:),1);

% display
if( show && nd<=3 )
  filterVisualize( dG, show );
  title( ['sigs=[' num2str(sigmas) '], derivs=[' num2str( nderivs ) ']']);
end

author: zouchuhang

repo: LayoutNet


text:
function mtrExportFibers(fg, fileName, xformToAcPc, xform, fiberToStrFlag, version)
% Exports fiber group to .pdb format file
%
%  mtrExportFibers(fg, fileName, [xformToAcPc], [xform=eye()],[fiberToStrFlag=true],[version=3])
%
% Writes a mrDiffusion fiber group structure to a pdb file.  You should
% always write to a version 3 file (the default) unless you know what you
% are doing and really need a version 2 file.
%
% If your fibers are already in AcPc space then you probably want to store
% them that way in the file as it makes everything else easier.  In that
% case don't set xformToAcPc and xform.  They will default to identity and
% leave your fibers in the AcPc space.
%
% fg:             Fiber group to save into PDB file
% fileName:       Filename for PDB file
% xformToAcPc:    Transform current fiber coordinates into AcPc space
% xform:          Transforms current fiber coordinates into the space we would like
%                 to store on file
% fiberToStrFlag: set to zero if want to use legacy code for saving pdb
%                 file
% version:        PDB 3 (default) or 2 (not desired)
%
% Examples: (See Debug configuration for testing)
%   fgName = fullfile(mrvDataRootPath,'quench','fibers','LeftArcuate_ver3.pdb');
%   fg = mtrImportFibers(fgName);
%   v = 3; mtrExportFibers(fg, 'deleteMe.pdb', [], [], [], v);
%
% Stanford VISTA Team

if ~exist('fiberToStrFlag', 'var')|| isempty(fiberToStrFlag) || fiberToStrFlag~=0
    fiberToStrFlag=true; %By defaultconvert the Matlab fiber data structure into one big contiguous array.
end

if ~exist('fileName','var') || isempty(fileName)
    fileName = mrvSelectFile('w', 'pdb', 'Write pdb fibers');
end

% Make sure extension is pdb
[p,n] = fileparts(fileName);
fileName = [fullfile(p,n),'.pdb'];

if ~exist('version', 'var')|| isempty(version), version = 3; end

if version ~= 2 && version ~= 3
	error('PDB version should be 2 or 3');
end

if(~exist('xformToAcPc','var') || isempty(xformToAcPc))
    xformToAcPc = eye(4);
end

% comment and explain, please
if(exist('xform','var') && ~isempty(xform) && (isstruct(xform) || ~all(all(xform==eye(4)))))
    fg = dtiXformFiberCoords(fg, xform);
    % xformToAcPc = xformToAcPc*inv(xform);
    xformToAcPc = xformToAcPc/xform;
end

% Create blank database for finding some of the sizes to write out
pdb = mtrPathwayDatabase();
% PDB structure will have stat header info, but fg will retain the actual
% statistics per path to record

%fg.params = {};
% This routine (see below) is empty.  If we want to check something, then
% let's figure out what and then check it. For now, I am commenting it out.
% [pdb,fg] = enforceRightStats(pdb,fg);
numstats = length( pdb.pathway_statistic_headers );
numalgs = 0;

% Set offset and open the file.
% Tragically the getAlgorithms... function is a single number that doesn't
% depend on the input argument.
%
offset = 4*4 + 16*8 + numstats * getStatisticsHeaderSize(pdb) + ...
    numalgs * getAlgorithmsHeaderSize(pdb); % 4 uint, 16 doubles
fid = fopen(fileName, 'wb');

% Writing header
fwrite(fid,offset,'uint');  % An integer that points to something ...
xformToAcPc = xformToAcPc';
fwrite(fid,xformToAcPc(:),'double');  % The xform
fwrite(fid,numstats,'uint');          % The number of statistics saved with the fibers

% Write statistics to the PDB header
saveStatisticsHeader(pdb,fid);

fwrite(fid,numalgs,'uint');  % Well, this is weird.  I guess numalgs is always zero.
% Write algorithm header
% saveAlgorithmHeader(pdb,fid);

fwrite(fid,version,'uint');

% Writing path info
numpaths = length(fg.fibers);
fwrite(fid,numpaths,'uint');
fileOffsets = zeros(1,numpaths);
for ff = 1:numpaths
    % Record current file position to write later
    fileOffsets(ff) = ftell(fid);

    path_offset = 3 * 4  + numstats*8; % 4 int, numstats double
    fwrite(fid,path_offset,'int');
    nn = size( fg.fibers{ff}, 2 );
    fwrite(fid,nn,'int');
    fwrite(fid,0,'int'); % algo type HACK
    fwrite(fid,1,'int'); % seedpointindex HACK ???

    % Write stats per path
    for as = 1:numstats
        fwrite(fid,fg.params{as}.stat(ff),'double');
    end

    % Writing path nodes
    %pos = fg.fibers{ff}.coords;
    pos = fg.fibers{ff};
    %pos = mrAnatXformCoords(xformFromAcpc, pos')';
    % Change from 1 based to 0 based positions in mm space
    %pos = pos - repmat([2,2,2]',1,size(pos,2));
    fwrite(fid,pos,'double');

    % Writing stats values per position
    for as = 1:numstats
        % Just write garbage as fiber groups can't handle per point stats
        if( pdb.pathway_statistic_headers(as).is_computed_per_point )
            error('Cannot export fiber group with per point statistics!');
            %fwrite(fid,this.pathways(p).point_stat_array(as,:),'double');
        end
    end
end

% Database Footer
fwrite(fid,fileOffsets,'uint64');
fclose(fid);

% Convert to string depending on the version.
if fiberToStrFlag
	if version == 2, str=fiberToStr(fg); end
	if version == 3, str=fiberToStr3(fg);end
    myf = fopen(fileName,'wb');fwrite(myf,str,'uint8');fclose(myf);%fclose(fid);return;
end
return;

%-----
function [pdb,fg] = enforceRightStats(pdb,fg)
%if ~isempty(fg.params) && ~ieNotDefined('fg.params{1}.name') && strcmpi(fg.params{1}.name,'weight')
% The above check is wrong because ieNotDefined cant work this way. This
% wrong check ends up destroying the params. Nor do I understand why this
% function always ends up overwriting the fibers.
return;

if ~isempty(fg.params) && isfield(fg.params{1},'name') && ~isempty(fg.params{1}.name)
    tempFGParams = fg.params(1);
    fg.params = {};
    fg.params{1}.name = 'Length';
    for ff=1:length(fg.fibers)
        fiber = fg.fibers{ff};
        if(size(fiber,2)<size(fiber,1)); fiber=fiber'; end
        mmPerStep = sqrt(sum((fiber(:,1)-fiber(:,2)).^2));
        fg.params{1}.stat(ff) = mmPerStep*(size(fiber,2)-1);
    end
    s=length(pdb.pathway_statistic_headers)+1;
    pdb = addStatisticHeader(pdb,fg.params{1}.name,'NA', 1, 0, 1, 8877+s);
    fg.params(2) = tempFGParams;
    s=length(pdb.pathway_statistic_headers)+1;
    pdb = addStatisticHeader(pdb,fg.params{2}.name,'NA', 1, 0, 1, 8877+s);
    for ii=3:3
        pdb = addBlankStat(pdb,'Blank');
        fg.params{ii}.name = 'Blank';
        fg.params{ii}.stat = zeros(1,length(fg.fibers));
    end
else
    fg.params = {};
    fg.params{1}.name = 'Length';
    for ff=1:length(fg.fibers)
        fiber = fg.fibers{ff};
        if(size(fiber,2)<size(fiber,1)); fiber=fiber'; end
        mmPerStep = sqrt(sum((fiber(:,1)-fiber(:,2)).^2));
        fg.params{1}.stat(ff) = mmPerStep*(size(fiber,2)-1);
    end
    for ii=2:3
        pdb = addBlankStat(pdb,'Blank');
        fg.params{ii}.name = 'Blank';
        fg.params{ii}.stat = zeros(1,length(fg.fibers));
        %     if numstats>=ii
        %         pdb = addBlankStat(pdb,fg.params{ii}.name);
        %     else
        %         pdb = addBlankStat(pdb,'Blank');
        %         fg.params{ii}.stat = zeros(1,length(fg.fibers));
        %     end
    end
end
return;

%-----
function pdb = addBlankStat(pdb,agg_name)
s=length(pdb.pathway_statistic_headers)+1;
pdb = addStatisticHeader(pdb,agg_name,'NA', 1, 0, 1, 8877+s);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [inF, n_r] = extend_dcm(inF,hA,hB,hC,hD,dim,sources)

%%- get dimensions
n=dim.n;
n_u=dim.n_u;
n_r = get_dims(hA,hB,hC,hD) ;

sourcesDim = cellfun(@numel,{sources(2:end).out});
for i=1:sum(sourcesDim)
    respSource(i) = sum(i>cumsum(sourcesDim))+1 ;
end

%%
[inFtemp] = prepare_dcm(hA,hB,hC,hD,respSource);

%% = save extended structure
%- raw matrices
inF.hA = inFtemp.A ;
inF.hB = inFtemp.B ;
inF.hC = inFtemp.C ;
inF.hD = inFtemp.D ;

%- parameter indices
inF.indhA = inFtemp.indA + inF.indself ;

for i=1:n_u
    inF.indhB{i} = inFtemp.indB{i} + inF.indself ;
end
inF.indhC = inFtemp.indC + inF.indself ;
for i=1:n
    inF.indhD{i} = inFtemp.indD{i} + inF.indself ;
end
inFtemp.indself = inFtemp.indself ;
inF.indhself = inFtemp.indself + inF.indself ;

%- indicators
inF.dhA = inFtemp.dA;
inF.dhB = inFtemp.dB;
inF.dhC = inFtemp.dC;
inF.dhD = inFtemp.dD;

%- save dimensions
% dim.n_r=n_r;

end

%% ========== subfunctions =============

%- compute number of predicted responses
function nr=get_dims(hA,hB,hC,hD)
nr=0;
nr=max(nr,size(hA,1));
for i=1:length(hB)
    nr=max(nr,size(hB{i},1));
end
nr=max(nr,size(hC,1));
for i=1:length(hD)
    nr=max(nr,size(hD{i},1));
end
end


text:
function morph_profile=GetMP(img,SEs)
% Extract Morphological Profiles of the input image with the given SEs
%2016-10-20, jlfeng
[nr,nc]=size(img);
numSE=length(SEs);
morph_profile=zeros(nr,nc,numSE*8);
%%
% Opening and Closing
mp1=zeros(nr, nc,numSE);mp2=zeros(nr, nc,numSE);
for kk=1:numSE
    mp1(:,:,kk)=imclose(img,SEs{kk});
    mp2(:,:,kk)=imopen(img,SEs{numSE+1-kk});
end
dmp=diff(cat(3,mp1,img,mp2),[],3);
morph_profile(:,:,1:numSE*4)=cat(3,mp1,mp2,dmp);
% OpenRec and CloseRec
for kk=1:numSE
    marker=imerode(img,SEs{numSE+1-kk});
    mp1(:,:,kk)=imreconstruct(marker,img);
    marker=imdilate(img,SEs{kk});
    mp2(:,:,kk)=imreconstruct(imcomplement(marker),imcomplement(img));
end
dmp=diff(cat(3,mp1,img,mp2),[],3);
morph_profile(:,:,(numSE*4+1):numSE*8)=cat(3, mp1,mp2,dmp);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function outputImage = cumsum(this, applicationDimension)
% Computes cumulative sum along specified dimension, uses Matlab cumsum function
%
%
%   Y = MrImage()
%   Y.cumsum(applicationDimension)
%
% This is a method of class MrImage.
%
% IN
%   applicationDimension    image dimension along which operation is
%                           performed (e.g. 4 = time, 3 = slices)
%                           default: The last dimension with more than one
%                           value is chosen
%                           (i.e. 3 for 3D image, 4 for 4D image)
%
% OUT
%   outputImage             cumsum of all images along application dimension
%
% EXAMPLE
%   cumsum
%
%   See also MrImage MrImage.perform_unary_operation

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

if nargin < 2
    applicationDimension = this.dimInfo.nDims;
else
    applicationDimension = this.dimInfo.convert_application_dimensions(...
        applicationDimension);
end

outputImage = this.perform_unary_operation(@(x) cumsum(x, applicationDimension));

author: translationalneuromodeling

repo: tapas


text:
function edges=meshedge(elem,varargin)
%
% edges=meshedge(elem,opt)
%
% return all edges in a surface or volumetric mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2011/02/26
%
% input:
%    elem:  element table of a mesh (support N-d space element)
%    opt: optional input, giving the additional options. If opt
%         is a struct, it can have the following field:
%       opt.nodeorder: if 1, assuming the elem node indices is in CCW
%                      orientation; 0 use nchoosek() output to order edges
%         you can replace opt by a series of ('param', value) pairs.
%
% output:
%    edge:  edge list; each row is an edge, specified by the starting and
%           ending node indices, the total edge number is
%           size(elem,1) x nchoosek(size(elem,2),2). All edges are ordered
%           by looping through each element first.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

dim=size(elem);
edgeid=nchoosek(1:dim(2),2);
len=size(edgeid,1);
edges=zeros(dim(1)*len,2);
for i=0:len-1
    edges((i*dim(1)+1):((i+1)*dim(1)),:)=[elem(:,edgeid(i+1,1)) elem(:,edgeid(i+1,2))];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [solver,diagnostic] = setupBNB(solver,ProblemClass,options,solvers,socp_are_really_qc,F,h,logdetStruct,parametric,evaluation_based,F_vars,exponential_cone,allsolvers)

diagnostic = [];

temp_options = options;
temp_options.solver = options.bnb.solver;
tempProblemClass = ProblemClass;
tempProblemClass.constraint.binary  = 0;
tempProblemClass.constraint.integer = 0;
tempProblemClass.constraint.semicont = 0;
localsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc,allsolvers);
if isempty(localsolver) | strcmpi(localsolver.tag,'bnb')
    if isempty(temp_options.bnb.solver)
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
    else
        % User has specified a lower-bound solver, but we failed to use
        % this. Could be that it doesn't exist, or that it is not
        % applicable
        for i = 1:length(solvers)
            if strcmpi(solvers(i).tag,temp_options.bnb.solver)
                % Solver exist, hence it is not applicable
                diagnostic.solvertime = 0;
                diagnostic.info = yalmiperror(-4,temp_options.bnb.solver);
                diagnostic.problem = -4;
                return
            end
        end
        % Solver was not found in list of available solvers, hance user
        % has specified a solver which doesn't exist
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-3,temp_options.bnb.solver);
        diagnostic.problem = -3;
        return
    end
    return
elseif strcmpi(localsolver.tag,'bmibnb')
    [localsolver,diagnostics] = setupBMIBNB(localsolver,tempProblemClass,options,solvers,socp_are_really_qc,F,h,logdetStruct,parametric,evaluation_based,F_vars,[],allsolvers);
end
solver.lower = localsolver;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function map = mrvColorMaps(name, nC, varargin)
%
% map = mrvColorMaps(name, [nColors=256]);
%
% Color maps used in mrVista.
%
% list = mrvColorMaps without input arguments
% returns a list of the currently-defined
% color maps.
%
% map = mrvColorMaps([name or num]) returns a 3x256
% color map corresponding to the name, or its index
% number in the list of all colormaps (useful for
% popup menus etc).
%
% map = mrvColorMaps('user') prompts for a user-defined
% colormap. This can be specified in one of two ways:
%   (1) Enter the # of colors to use. A prompt pops up
%       with a GUI to edit manually each of the colors.
%       You can also invoke the matlab CMAPEDITOR on this
%       for further options.
%
%   (2) Save a matlab file which contains the variable
%       'cmap'. At the dialog, leave the # of colors empty.
%       A second dialog will appear prompting for this file.
%
% map = mrvColorMaps([name or num], nC) returns a
% 3xnColors preset color map.
%
% For the case 'redgreenblue', you can specify an optional 3rd argument,
% gamma, which specifies the sharpness of the color transitions. (gamma > 1
% is sharper, gamma==1 is linear, gamma < 1 is more gradual).
%
% mrvColorMaps('demo') shows a set of all current colormaps in a figure.
%
% ras 07/05.
% ras 10/05: added support for external color maps.
% ras 05/08: added 'short HSV' color map, support for MATLAB expression
if nargin==0
	% output a list of names
	map = {'hot' 'cool' 'jet' 'hsv' 'hsvshort' 'gray' 'autumn' 'winter'...
		'red' 'green' 'blue' 'yellow' 'cyan' 'magenta' ...
		'red binary' 'green binary' 'blue binary' ...
		'redgreenblue' 'redyellowblue' 'greenbluered' ...
		'redgreenblueyellow' 'redmagentablue' ...
		'blueredyellow' 'greenred' ...
		'revjet' 'bluered' 'coolhot' 'user'};
	return
end

if ~exist('nC','var') || isempty(nC), nC = 256; end

if isnumeric(name)
	% convert from index # into name
	list =  mrvColorMaps;
	name = list{name};
end

map = zeros(nC,3);
switch lower(name)
	case 'hot', map = hot(nC);
	case 'cool', map = cool(nC);
	case 'jet', map = jet(nC);
	case 'hsv', map = hsv(nC);
	case 'hsvshort', map = circshift( hsv(round(nC*1.4)), 10 ); map = map(1:nC,:);
	case 'gray', map = gray(nC);
	case 'autumn', map = autumn(nC);
	case 'winter', map = winter(nC);
	case 'red binary', map(:,1) = 1;
	case 'green binary', map(:,2) = 1;
	case 'blue binary', map(:,3) = 1;
	case 'red', map(:,1) = linspace(0, 1, nC);
	case 'green', map(:,2) = linspace(0, 1, nC);
	case 'blue',  map(:,3) = linspace(0, 1, nC);
	case 'yellow', map = mrvColorMaps('red',nC) + mrvColorMaps('green',nC);
	case 'cyan', map = mrvColorMaps('blue',nC) + mrvColorMaps('green',nC);
	case 'magenta', map = mrvColorMaps('red',nC) + mrvColorMaps('blue',nC);

	case {'redgreen', 'greenred'},
		map(:,1) = linspace(0.2,1,nC);
		map(:,2) = linspace(1,0.2,nC);

	case {'redgreenblue', 'rgb'},
		% do some gamma adjust so the color trabsitions are more prominent
		% gamma coefficient (<1 = sharp, 1=linear, >1 = gradual)
		if ~isempty(varargin)
			gamma = varargin{1};
		else
			gamma = 2.8;
		end

		% make ramps for each color
		r = mrvColorMaps('red', nC) .^ gamma;
		g = [mrvColorMaps('green', ceil(nC/2)) .^ gamma] .* 0.7;
		g = [g; flipud(g)]; g = g(1:nC,:);
		b = flipud( mrvColorMaps('blue', nC) ) .^ gamma;
		map = r + g + b;
		map = map(1:nC,:);

	case {'redyellowblue', 'ryb'}
		% complement to rgb color map, for representing the ipsilateral
		% visual field.

		% do some gamma adjust so the color trabsitions are more prominent
		% gamma coefficient (<1 = sharp, 1=linear, >1 = gradual)
		if ~isempty(varargin)
			gamma = varargin{1};
		else
			gamma = 2.8;
		end

		% make ramps for each color
		r = mrvColorMaps('red', nC) .^ gamma;
		y = [mrvColorMaps('yellow', ceil(nC/2)) .^ gamma] .* 0.7;
		y = [y; flipud(y)]; y = y(1:nC,:);
		b = flipud( mrvColorMaps('blue', nC) ) .^ gamma;
		map = r + y + b;
		map = map(1:nC,:);

	case {'pinwheel' 'rgb_alt' 'redgreenblueyellow' 'rgby'}
		% an alternate version of the RGB color map above: still trying to
		% keep (red, green, blue) mapping to (down, horizontal, up) in
		% representing polar angle, but with an accommodation for
		% representing activations at ipsilateral polar angles as well.
		n = ceil(nC/4); m = nC - 2*n;
		map = [ mrvColorMaps('blue', n); ...
			flipud(mrvColorMaps('hsvshort', m)); ...
			flipud(mrvColorMaps('red', n)) ];
% 		n = ceil(nC/2); m = nC - n;
% 		map = [ mrvColorMaps('rgb', n); ...
% 				flipud(mrvColorMaps('ryb', m)) ];

	case {'greenbluered', 'gbr'},
		% do some gamma adjust so the color trabsitions are more prominent
		if ~isempty(varargin)
			gamma = varargin{1};
		else
			gamma = 2.8;  % gamma coefficient (<1 = sharp, 1=linear, >1 = gradual)
		end

		g = flipud( mrvColorMaps('green', nC) ) .^ gamma;
		b = mrvColorMaps('blue', ceil(nC/2)) .^ gamma;
		b = [b; flipud(b)]; b = b(1:nC,:);
		r = mrvColorMaps('red', nC) .^ gamma;    % make ramps for each color

		I = floor(nC/3):ceil(2*nC/3);

		map = r + g + b;

		map = map(1:nC,:);

	case {'redmagentablue', 'rmb', 'redblue', 'rb'},
		map(:,1) = linspace(1, 0, nC);
		map(:,2) = zeros(nC, 1);
		map(:,3) = linspace(0, 1, nC);

	case {'blueredyellow' 'bry' 'rby' 'redblueyellow'}
		r = linspace(0,1,nC/2)';
		z = zeros(nC/2,1);
		o = ones(nC/2,1);
		map = [[r;o] [z;r] flipud([z;r])];

	case {'revjet' 'reversejet'}
		map = 1 - jet(nC);

	case {'bluered' 'corr'}
		% blue-red colormap which works well for correlation coefficients
		blue = flipud( mrvColorMaps('blue', ceil(nC/2)) );
		red = mrvColorMaps('red', nC - ceil(nC/2));
		map = [blue; red];

	case {'coolhot' 'blueredmagenta' 'corr2'}
		% modified version of the above correlation colormap:
		% saturates to magenta at both ends
		nQ = ceil(nC / 4);  % one quarter of the cmap

		mag1 = [ zeros(1, nQ)' linspace(1, 0, nQ)' ones(1, nQ)'];
		blue = flipud( mrvColorMaps('blue', nQ) );
		red = mrvColorMaps('red', nC - 3*nQ);
		mag2 = [ones(1, nQ)' linspace(0, 1, nQ)' zeros(1, nQ)'];

		map = [mag1; blue; red; mag2];

	case 'user', map = mrvUserColormap;

	case 'demo', % show a demo of the available color maps
		map = mrvColorMaps;
		figure('Name', 'Color Maps in mrvColorMaps');
		for i = 1:27
			subplot(7, 4, i);
			image(ind2rgb(1:256, mrvColorMaps(i)));
			title(map{i});
			axis off;
		end
		return

	case 'expression' % let user input a MATLAB expression
		msg = ['Enter a MATLAB expression for the color map. ' ...
			'The color map should be a matrix of 3 x ncolors, ' ...
			'with values from 0 - 255. ' ...
			'You can directly enter the expression, or assign ' ...
			'it to a variable named "map". ' ...
			'E.g.:    [gray(128); hot(128)]; ' ...
			'or:      tmp = hsv(256);  map = tmp(1:100,:); '];
		def = sprintf('gray(%i); ', nC);
		resp = inputdlg({msg}, mfilename, 4, {def});
		if ~isempty(strfind(resp{1}, '=')) && ...
			~isempty(strfind(resp{1}, 'map'))
			% presume that the expression is assigning a value to the
			% 'map' variable:
			eval(resp{1})
		else
			% the map assignment is implied
			map = eval(resp{1});
		end

	otherwise
		warning('Invalid name specified. Returning empty color map');
end

% bounds check
map(map>1) = 1; map(map<0) = 0;

return
% /-------------------------------------------------------------------/ %

% /-------------------------------------------------------------------/ %
function map = mrvUserColormap;
% Use the colormap editor to get a colormap
%
% map = mrvUserColormap;
%
% ras 10/2005.
q='Enter # of colors in new colormap (leave empty to load from file)';
resp = inputdlg(q, 'mrvColormaps', 1, {'256'});
if isempty(resp{1})
	% load from file
	p = 'Select a .mat file with the colormap as a ''cmap'' variable.';
	[f pth] = uigetfile('*.mat',p);
	try
		load(fullfile(pth,f),'cmap');
		map = cmap;
	catch
		msg='Can''t find file w/ cmap variable, returning empty cmap.';
		warning(msg); map=[]; return
	end

else
	nC = str2num(resp{1});
	cbar = cbarEdit( cbarDefault(hot(nC)) );
	map = cbar.cmap;

end

return

% old cmaps
%     case 'selectivity'   % for AdaptNSelect experiment (4 colors), gum
%         colors = {[1 0 0] [1 1 0] [0 1 1] [0 0 1]};
%         colPerCond = nC / length(colors);
%         map = [];
%         for i = 1:length(colors)
%             for j = 1:colPerCond
%                 col = colors{i};
%                 w = 0.3 + 0.5*(j-1) / colPerCond; % weight of color
%                 map(end+1,:) = (w*col+ ((colPerCond-j+1) / (colPerCond))*[.7 .7 .7]);
%             end
%         end
%
%         % shift so that, when this is mapped to an integer range,
%         % integer values (like 1.0) will map to the max of the previous
%         % range, not the first of the next range -- also add black:
%         map = [0 0 0; map(1:nC-1,:)];


text:
function [f,h, whichSlice] = sagittalSlices(tv, ccfCoords, f)

if nargin<3
    f = figure;
    p = get(f, 'Position');
    set(f, 'Position',   [p(1) p(2)      1597         700]);
    set(f, 'Color', 'k');
    slicesExist = false;
else
    slicesExist = true;
end

slicePoints = 0.25:0.25:3.0; %mm
bregma = allenCCFbregma(); %voxels
sliceCCF = -slicePoints*1000+bregma(3)*10;

cc = ccfCoords;

r = randn(size(cc,1),2)*45;

[~,whichSlice] = min(abs(cc(:,3)-sliceCCF),[],2);

h = [];

apCoords = (0:size(tv,1)-1)*10;
dvCoords = (0:size(tv,2)-1)*10;

inclAP = apCoords>1000 & apCoords<11000;
inclDV = dvCoords<6000;

for s = 1:numel(sliceCCF)
    subtightplot(3,4,s);
    if ~slicesExist
        imagesc(apCoords(inclAP), dvCoords(inclDV), ...
            tv(inclAP,inclDV,round(sliceCCF(s)/10))');
        colormap gray
        axis image
        caxis([0 400]);
        axis off
    end

%     title(sprintf('ML = %.2f', slicePoints(s)));

    hold on;

    inclN = whichSlice==s & cc(:,3)>0;

    h(s) = scatter(cc(inclN,1)+r(inclN,1),cc(inclN,2)+r(inclN,2),4,'r','filled');
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function [montage, cfg] = ft_nirs_prepare_ODtransformation(cfg, data)

% FT_NIRS_PREPARE_ODTRANSFORMATION returns the transformation matrix from
% optical densities (OD) to chromophore concentrations such as (de-)
% oxygenated hemoglobin.
%
% Use as
%   [montage] = ft_prepare_ODtransformation(cfg, data)
%
% It is necessary to input the data on which you want to perform the
% inverse computations, since that data generally contain the optode
% information and information about the channels that should be included in
% the transformation. The data structure can be either obtained
% from FT_PREPROCESSING, FT_FREQANALYSIS, FT_TIMELOCKANALYSIS or
% FT_COMPONENTANALYSIS. If the data is empty, all channels will be included
% in transformation.
%
% The configuration should contain
%   cfg.channel            = Nx1 cell-array with selection of channels
%                            (default = 'nirs'), see FT_CHANNELSELECTION for
%                            more details
%
% Optional configuration settings are
%   cfg.age                = scalar, age of the subject (necessary to
%                            automatically select the appropriate DPF, or
%   cfg.dpf                = scalar, differential path length factor
%   cfg.dpffile            = string, location to a lookup table for the
%                            relation between participant age and DPF
%
% Note that the DPF might be different across channels, and is usually
% contained in the optode structure contained in the data.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_NIRS_TRANSFORM_ODS, FT_APPLY_MONTAGE

% options to be implemented:
%
% The NIRS positions can be present in the data or can be specified as
% cfg.opto      = structure with optode positions, see FT_DATATYPE_SENS
% cfg.siunits   = ft_getopt(cfg, 'siunits', 'no');   % yes/no, ensure that SI units are used consistently
% cfg.logarithm = string, can be 'natural' or 'base10' (default =

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2015-2016 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com
%
% Main programmer: Jörn M. Horschig
% $Id$

revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function

% the ft_preamble function works by calling a number of scripts from
% fieldtrip/utility/private that are able to modify the local workspace

ft_defaults                   % this ensures that the path is correct and that the ft_defaults global variable is available
ft_preamble init              % this will reset ft_warning and show the function help if nargin==0 and return an error
ft_preamble debug             % this allows for displaying or saving the function name and input arguments upon an error
ft_preamble loadvar    data   % this reads the input data in case the user specified the cfg.inputfile option
ft_preamble provenance data   % this records the time and memory usage at the beginning of the function

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% check if the input data is valid for this function
data = ft_checkdata(data);

% set the defaults
cfg.channel = ft_getopt(cfg, 'channel', 'nirs');
cfg.age     = ft_getopt(cfg, 'age', []);
cfg.dpf     = ft_getopt(cfg, 'dpf', []);

% get the optode definition
opto = ft_fetch_sens(cfg, data);

% select the appropriate channels
if isfield(data, 'topolabel')
  % the data reflects a componentanalysis, where the topographic and the
  % timecourse labels are different
  cfg.channel = ft_channelselection(cfg.channel, data.topolabel);
elseif isfield(data, 'label')
  % In the subsequent code, the matching channels in the sensor array and
  % in the configuration will be selected. To ensure that these channels
  % are also present in the data, update the configuration to match the data.
  cfg.channel = ft_channelselection(cfg.channel, data.label);
else
  % update the selected channels based on the optode definition
  cfg.channel = ft_channelselection(cfg.channel, opto.label);
end
cfg.channel = ft_channelselection('nirs', cfg.channel); % only NIRS valid

% check current data type by looking what's inside of []
typeidx = regexp(cfg.channel, sprintf('%s%s', regexptranslate('wildcard', '[*]'), '$'));

% remove channels that do not obey to labelling guidelines
illegalidx = cellfun(@isempty, typeidx);
cfg.channel(illegalidx) = [];

if isempty(cfg.channel)
  ft_error('no valid NIRS channels found')
end

% channel indices wrt optode structure
chanidx     = match_str(opto.label, cfg.channel);

% check on DPF values
if ~isempty(cfg.age) && ~isempty(cfg.dpf)
  error('unsure if age or dpf in configuration should take precedence');
elseif ~isempty(cfg.age)
  error('the use of cfg.age is not implemented, yet');
elseif ~isempty(cfg.dpf)
  % this is where they should be
  dpfs = repmat(cfg.dpf, size(cfg.channel));
elseif isfield(opto, 'DPF')
  % this is for backward compatibility
  dpfs = opto.DPF(chanidx);
elseif isfield(data, 'opto') && isfield(data.opto, 'DPF')
  % this is for backward compatibility
  dpfs = data.opto.DPF(chanidx);
elseif isfield(data, 'hdr') && isfield(data.hdr, 'opto') && isfield(data.hdr.opto, 'DPF')
  % this is for backward compatibility
  dpfs = data.hdr.opto.DPF(chanidx);
end

% which chromophores are desired
chromophoreIdx = [3 2];
chromophoreName = {'O2Hb' 'HHb'};

%% transform to concentrations or to OD
% read in the coefficient table
% FIXME put this into an own function to read this out
fid = fopen(fullfile(fileparts(mfilename('fullpath')), 'private', 'Cope_ext_coeff_table.txt'));
coefs = cell2mat(textscan(fid, '%f %f %f %f %f'));

% extract all optode combinations that are relevant here
tratra         = opto.tra(chanidx, :)';
transmitteridx = tratra>0;
receiveridx    = tratra<0;
optodeidx      = (transmitteridx | receiveridx)'; % transpose to get back to tra order

% extract the wavelengths
wavelengths = opto.wavelength(tratra(transmitteridx));
wlidx = bsxfun(@minus, coefs(:, 1), wavelengths);

% find the relevant channel combinations
chancmb = logical([]); % we need to loop here to get size of chancmb
chanUsed = zeros(numel(chanidx), 1);
for c=1:numel(chanidx)
  % skip channel if it already was in another channel combination
  if chanUsed(c)
    continue;
  end
  % compute the channel combinations
  tupletidx = sum(bsxfun(@minus, optodeidx, optodeidx(c, :))~=0, 2)==0;
  chanUsed = chanUsed|tupletidx;
  chancmb(:, end+1) = tupletidx;
end

% do the transformation
tra      = zeros(size(chancmb, 2)*numel(chromophoreIdx), size(cfg.channel, 1));
labelnew = cell(size(chancmb, 2)*numel(chromophoreIdx), 1);

% transformation has to be done per channel combination
chanidx = []; % we will use this variable here for indexing channels (Rx-Tx cmbs)
optoidx = []; % we will use this variable here for indexing optodes (individual Rx and Tx)
for c=1:size(chancmb, 2)

  % find all other channels of the exact same transmitter / receiver pair
  chanidx = chancmb(:, c);

  % extract coefficient idx of these channels
  [coefidx, colidx] = find(wlidx(:, chanidx)==0);

  % compute the transmitter/receiver distance in cm
  optoidx = find(chanidx, 1, 'first'); % we can take 'first' because the transmitter-optodes have to be physically in the exact same spot to form "one" channel
  dist = sqrt(sum(diff(opto.optopos(optodeidx(optoidx, :), :)).^2));

  % select dpf
  dpf = mean(dpfs(chanidx));

  % compute distance factor for correct scaling
  distfactor = dist/10*dpf/10; % in mm

  % select the respective coefficients
  eps=coefs(coefidx,chromophoreIdx)/10; % in mm

  % add to transformation matrix
  traIdx = 1+(c-1)*numel(chromophoreIdx):(c)*numel(chromophoreIdx);
  tra(traIdx, chanidx) = pinv(eps)/(distfactor);

  % create the new labels for the channels
  chanorg = cfg.channel(find(chanidx, 1, 'first'));
  for n=1:numel(chromophoreIdx)
    labelnew(traIdx(n), 1) = regexprep(chanorg, sprintf('%s%s', regexptranslate('wildcard', '[*]'), '$'), sprintf('[%s]', chromophoreName{n}));
  end

end

%% create output
montage = [];
montage.labelorg = cfg.channel(:)';
montage.labelnew = labelnew;
montage.tra      = tra;


text:
% Forward and Inverse Fast Wavelet Transforms
%
% Usage:
%   output = wlet( signal, [output_type [, num_coeff [, direction]]])
%
%
%   Parameters:
%
%      output_type is one of:
%
%         1  Output is one dimensional, and represents the raw result
%            from the fast wavelet transform algorithm (wavelet coefficients).
%            This is the default.
%
%
%         2  Output is two dimensional. This output format is convenient for
%  	     displaying transformed data as a time-scale matrix (scalogram).
%
%      num_coeff
%
%         gives the number of filter coeffients used in approximating
%         the wavelet.  Possible values are 4, 12 or 20, and the default is 20.
%
%      direction
%
%         1  Transform is forward, from time domain to time-scale
%	     domain. (default)
%
%        -1  Transform is reverse, from time-scale domain to time domain.
%
%
% Two dimensional output can only be used with the forward transform, and
% the input to both forward and reverse transforms must be one-dimensional.
%
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
function [ output_struct ] = parsepdvinputs( varargin )
%PARSEPDVINPUTS Parse input arguments for pdv
%
%   Used in both in-memory and in-file version of PDV
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

p = inputParser;
p.KeepUnmatched=true;
p.addParamValue('deltax', 0.25, @(x) isscalar(x)&&(x>0));
p.addParamValue('filtersize', 5, @(x) isvector(x)&&(length(x)<=2)&&all(x>0));
p.addParamValue('filtertype', 'mean', @(x)any(strcmpi(x,{'mean','median'})));
p.addParamValue('dim', 1, @(x) isequal(x,1)||isequal(x,2));
p.addParamValue('framenumber',1);
p.FunctionName = 'PDV';
p.parse(varargin{:});

output_struct=p.Results;

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [trackA, trackB] = TrimTracks(trackA, trackB)
nFrame = min(trackA.nFrame, trackB.nFrame);

goodA = trackA.head <= nFrame;
trackA.points = trackA.points(goodA, 1:nFrame, :);
trackA.head = trackA.head(goodA);
trackA.tail = trackA.tail(goodA);
trackA.tail = min(trackA.tail, nFrame);
trackA.len = trackA.tail - trackA.head + 1;
trackA.nFrame = nFrame;
trackA.nTrack = length(trackA.head);

goodB = trackB.head <= nFrame;
trackB.points = trackB.points(goodB, 1:nFrame, :);
trackB.head = trackB.head(goodB);
trackB.tail = trackB.tail(goodB);
trackB.tail = min(trackB.tail, nFrame);
trackB.len = trackB.tail - trackB.head + 1;
trackB.nFrame = nFrame;
trackB.nTrack = length(trackB.head);

end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function [picks,P] = selectPF(tol,n)

halfn = round(n/2);
[I,J] = meshgrid(1:n,1:n);
P = abs(I - halfn*(rand(n) + .5)) < tol & ...
    abs(J - halfn*(rand(n) + .5)) < tol;

P(halfn+1:n,:) = 0;
P(halfn:halfn+1,:) = 1;
P(:,halfn:halfn+1) = 1;
P = ifftshift(P);
P(1,1) = 1;

picks = find(P);

repo: REST


text:
function a = c8mat_identity ( n )

%*****************************************************************************80
%
%% C8MAT_IDENTITY sets the square matrix A to the identity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, complex A(N,N), the matrix.
%
  a(1:n,1:n) = 0.0;

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = year_is_embolismic_greek ( y )

%*****************************************************************************80
%
%% YEAR_IS_EMBOLISMIC_GREEK returns TRUE if the Greek year was embolismic.
%
%  Discussion:
%
%    Apparently, the Greek calendar was emended haphazardly.  This
%    routine does not attempt to follow that historical pattern, and
%    just uses the Hebrew calendar pattern for now.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, logical VALUE, TRUE if the year was embolismic.
%
  if ( 12 <= i4_modp ( 7 * y + 13, 19 ) )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function output = SHA1 ()

fname=input('Input File (in ASCII format)? ','s');
hash_foutname=input('Output File for SHA Hash? ','s');
time1=clock;
%Open the input file and get the first line of data
fid=fopen(fname);
M = fread(fid);
fclose(fid);
ini=dec2bin(M(1),8);

for ii = 2:length(M)
    ini=cat(2,ini,dec2bin(M(ii),8));
end

s2=length(ini)/8;

    block_temp = [ ini, ...
                 '1', ...
                 num2str(zeros(mod(448-1-s2*8,512),1))' ...
                dec2bin(s2*8,64) ];
           nb=length(block_temp)/512;

block = reshape(block_temp,512,nb)';

% H(0) initialization
H = [ '67','45','23','01';...
      'ef','cd','ab','89';...
      '98','ba','dc','fe';...
      '10','32','54','76';...
      'c3','d2','e1','f0'] ;

for N = 1:length(block_temp)/512

	% Generate Wt
	Wt = WtgenSHA1( block(N,:) );

	%use previous last H's
	a_hex = H(1,:); a = hex2bin(a_hex);
	b_hex = H(2,:); b = hex2bin(b_hex);
	c_hex = H(3,:); c = hex2bin(c_hex);
	d_hex = H(4,:); d = hex2bin(d_hex);
	e_hex = H(5,:); e = hex2bin(e_hex);

	% Kt initialization

	for t = 1:80
        T1 = bin2dec2(cls(a,5) );
        %Number to Char
        if t <= 20
   f = xor(b&c, ~b&d);
elseif t <= 40
   f = xor(xor(b,c),d);
elseif t <= 60
    f = xor(xor(b&c,b&d),c&d);
elseif t <= 80
    f = xor(xor(b,c),d);
else
    error ('t must be in the range 1 to 80')
end

        T2 = bin2dec2(f);
        T3 = bin2dec2(e);
        if t>=60
            T4=floor(2^32*sqrt(10));
        elseif t>=40
            T4=floor(2^32*sqrt(5));
        elseif t>=20
            T4=floor(2^32*sqrt(3));
        elseif t>0
            T4=floor(2^32*sqrt(2));
        end
%        T4=bin2dec2(Kt(t,:));
        T5 = bin2dec2(Wt(t,:) );
        T = mod(T1+T2+T3+T4+T5,2^32);

        e = d;
        d = c;
        c = cls(b,30);
        b = a;
        a = dec2bin(T,32);
        %display_resultsSHA1(t,a,b,c,d,e)         %receives column bits
	end

	H = calculate_new_HSHA1(H,a,b,c,d,e);

end
output=H;

time2=clock;
disp(etime(time2,time1));

author: Sable

repo: mcbench-benchmarks


text:
function FrontNo = UpdateFront(PopObj,FrontNo,x)
% Update the front number of each solution when a solution is added or
% deleted

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

    [N,M] = size(PopObj);
    if nargin < 3
        %% Add a new solution (has been stored in the last of PopObj)
        FrontNo  = [FrontNo,0];
        Move     = false(1,N);
        Move(N)  = true;
        CurrentF = 1;
        % Locate the front No. of the new solution
        while true
            Dominated = false;
            for i = 1 : N-1
                if FrontNo(i) == CurrentF
                    m = 1;
                    while m <= M && PopObj(i,m) <= PopObj(end,m)
                        m = m + 1;
                    end
                    Dominated = m > M;
                    if Dominated
                        break;
                    end
                end
            end
            if ~Dominated
                break;
            else
                CurrentF = CurrentF + 1;
            end
        end
        % Move down the dominated solutions front by front
        while any(Move)
            NextMove = false(1,N);
            for i = 1 : N
                if FrontNo(i) == CurrentF
                    Dominated = false;
                    for j = 1 : N
                        if Move(j)
                            m = 1;
                            while m <= M && PopObj(j,m) <= PopObj(i,m)
                                m = m + 1;
                            end
                            Dominated = m > M;
                            if Dominated
                                break;
                            end
                        end
                    end
                    NextMove(i) = Dominated;
                end
            end
            FrontNo(Move) = CurrentF;
            CurrentF      = CurrentF + 1;
            Move          = NextMove;
        end
    else
        %% Delete the x-th solution
        Move     = false(1,N);
        Move(x)  = true;
        CurrentF = FrontNo(x) + 1;
        while any(Move)
            NextMove = false(1,N);
            for i = 1 : N
                if FrontNo(i) == CurrentF
                    Dominated = false;
                    for j = 1 : N
                        if Move(j)
                            m = 1;
                            while m <= M && PopObj(j,m) <= PopObj(i,m)
                                m = m + 1;
                            end
                            Dominated = m > M;
                            if Dominated
                                break;
                            end
                        end
                    end
                    NextMove(i) = Dominated;
                end
            end
            for i = 1 : N
                if NextMove(i)
                    Dominated = false;
                    for j = 1 : N
                        if FrontNo(j) == CurrentF-1 && ~Move(j)
                            m = 1;
                            while m <= M && PopObj(j,m) <= PopObj(i,m)
                                m = m + 1;
                            end
                            Dominated = m > M;
                            if Dominated
                                break;
                            end
                        end
                    end
                    NextMove(i) = ~Dominated;
                end
            end
            FrontNo(Move) = CurrentF - 2;
            CurrentF      = CurrentF + 1;
            Move          = NextMove;
        end
        FrontNo(x) = [];
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function resized_patch = sample_patch(im, pos, sample_sz, output_sz, gparams)

if nargin < 4
    output_sz = [];
end
if nargin < 5 || ~isfield(gparams, 'use_mexResize')
    gparams.use_mexResize = false;
end

% Pos should be integer when input, but floor in just in case.
pos = floor(pos);

% Downsample factor
resize_factor = min(sample_sz ./ output_sz);
df = max(floor(resize_factor - 0.1), 1);
if df > 1
    % pos = 1 + of + df * (npos - 1)

    % compute offset and new center position
    os = mod(pos - 1, df);
    pos = (pos - 1 - os) / df + 1;

    % new sample size
    sample_sz = sample_sz / df;

    % donwsample image
    im = im(1+os(1):df:end, 1+os(2):df:end, :);
end

% make sure the size is not too small and round it
sample_sz = max(round(sample_sz), 2);

xs = pos(2) + (1:sample_sz(2)) - floor((sample_sz(2)+1)/2);
ys = pos(1) + (1:sample_sz(1)) - floor((sample_sz(1)+1)/2);

%check for out-of-bounds coordinates, and set them to the values at
%the borders
xs(xs < 1) = 1;
ys(ys < 1) = 1;
xs(xs > size(im,2)) = size(im,2);
ys(ys > size(im,1)) = size(im,1);

%extract image
im_patch = im(ys, xs, :);

if isempty(output_sz) || isequal(sample_sz(:), output_sz(:))
    resized_patch = im_patch;
else
    if gparams.use_mexResize
        resized_patch = mexResize(im_patch, output_sz, 'linear');
    else
        resized_patch = imresize(im_patch, output_sz, 'bilinear', 'Antialiasing',false);
    end
end

end

author: martin-danelljan

repo: ECO


text:
function varargout = sample(f, varargin)
%SAMPLE   Values of a CHEBFUN3 object on a tensor product grid.
%   X = SAMPLE(F) returns the tensor of values of F on a tensor product
%   grid.
%
%   [CORE, C, R, T] = SAMPLE(F) returns the low rank representation of the
%   values of F on a tensor product grid so that X = CORE x_1 C x_2 R x_3 T.
%
%   [CORE, C, R, T] = SAMPLE(F, M, N, P) returns the values of F on an
%   M x N x P tensor product grid.
%
% See also CHEBFUN3/FEVAL.

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

% Empty check.
if ( isempty(f) )
    varargout = {[]};
    return
end

if ( nargin == 4 )
     m = varargin{1};
     n = varargin{2};
     p = varargin{3};
else
    [m,n,p] = length(f);
    m = max(m, 51);
    n = max(n, 51);
    p = max(p, 51);
end

% Use Slice-Tucker decomposition so we can keep it in low rank form:
[fCore, fCols, fRows, fTubes] = tucker(f);
Cvals = sample(fCols, m);
Rvals = sample(fRows, n);
Tvals = sample(fTubes, p);

% Evaluate:
if ( nargout <= 1 )
    varargout = {chebfun3.txm(chebfun3.txm(chebfun3.txm(fCore, Cvals, ...
        1), Rvals, 2), Tvals, 3)};
else
    varargout = {fCore, Cvals, Rvals, Tvals};
end

end

author: chebfun

repo: chebfun


text:
function pass = test_pcg( pref )
% Test operator PCG method.

% Tolerance:
if ( nargin == 0 )
    pref = cheboppref;
end
tol = 1e2*pref.bvpTol;

%% Example 1: (-u_xx=f, bc=0, mean(f)==0)
a = @(x) 1;
c = @(x) 1;
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
N = chebop( L );
N.bc = 0;
f = chebfun(@(x) 1-3*x.^2);

u = N\f;
% Chebfun-ready pcg:

v = pcg(N, f);

% Errors:
pass(1) = ( norm( u - v ) < tol );

%% Example 2: (-(a(x)*u_x)_x=f, bc=0, mean(f)==0)
a = @(x) 2+cos(pi*x);
c = @(x) 0*x;
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
N = chebop( L );
N.bc = 0;
f = chebfun(@(x) 1-3*x.^2);
u = N \ f;

% Chebfun-ready pcg:
v = pcg(N, f);

% Errors:
pass(2) = ( norm( u - v ) < tol );

%% Example 3: (-(a(x)*u_x)_x + c(x)*u=f, bc=0, mean(f)==0)

a = @(x) 2+cos(pi*x);
c = @(x) 1+10*x.^2;
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
N = chebop( L );
N.bc = 0;
f = chebfun(@(x) 1-3*x.^2);
u = N \ f;
% Chebfun-ready pcg:
v = pcg(N, f);

% Errors:
pass(3) = ( norm( u - v ) < tol );

%% Example 4: (-(a(x)*u_x)_x + c(x)*u=f, bc=0, sum(f)==0, piecewise)
a = @(x) 2+cos(pi*abs((x-0.5)));
c = @(x) 1+10*x.^2;
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
N = chebop( L );
N.bc = 0;
%f = chebfun(@(x) 1-2*x.^2);
f = chebfun(@(x) abs(x) - 0.5 ,'splitting','on');
u = N \ f;
% Chebfun-ready pcg:
v = pcg(N, f);
pass(4) = ( norm(u - v) < tol );

%% Example 5: (-(a(x)*u_x)_x + c(x)*u=f, bc~=0, mean(f)~=0)
a = @(x) 2+cos(pi*x);
c = @(x) 1+10*x.^2;
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
N = chebop( L );
N.lbc = 1;
N.rbc = -1;
f = chebfun(@(x) 1-2*x.^2);
u = N \ f;

% Chebfun-ready pcg:
v = pcg(N, f);

pass(5) = ( norm(u - v) < tol );

%% Example 6: (-(a(x)*u_x)_x + c(x)*u=f, bc~=0, mean(f)~=0, piecewise)
a = @(x) 2+cos(pi*x);
c = @(x) abs(cos(2*pi*x));
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
N = chebop( L );
N.lbc = 1;
N.rbc = -1;
f = chebfun(@(x) 1-2*x.^2);
u = N \ f;

% Chebfun-ready pcg:
v = pcg(N, f);

pass(6) = ( norm(u - v) < tol );

%% Example 7: (-(a(x)*u_x)_x + c(x)*u=f, bc~=0, mean(f)~=0, nonstandard domain)
a = @(x) 2+cos(pi*x);
c = @(x) abs(cos(2*pi*x));
L = @(x,u) -diff(a(x).*diff(u)) + c(x).*u;

% Chebop solve:
dom = [-3/2,2];
N = chebop( L, dom );
N.lbc = 3;
N.rbc = -1;
f = chebfun(@(x) 1-2*x.^2, dom);
u = N \ f;

% Chebfun-ready pcg:
v = pcg(N, f);

pass(7) = ( norm(u - v) < tol );

if ( all(pass) == 1 )
    pass = 1;
end

end


text:
clear all; close all; clc;

[filename1,pathname]=uigetfile('*.*','select the image');
image1=imread(num2str(filename1));
figure(1);
imshow(image1);	title('original image');     % orginal image for watermarking
image1=double(image1);

[row,col]=size(image1);

imagew=imread('dmg2.tif');

[marked]=blockdwt2(image1,imagew);  % generates the watermarked image

markedmax = max(marked(:));
markscale = marked/markedmax*255;

figure(2);
colormap(gray(256));
image(marked);							% shows the watermarked image
title('Watermarked image');
imwrite(marked,gray(256),'marked_image.bmp');	% saves the watermarked image as a bmp file
figure(3);
watermark=image1-marked;			% image adaptive watermark
watermark=watermark*255/max(watermark(:));

for i = 1:row                       % thresholding
     for j = 1:col
        if watermark(i,j) > 70
            watermark(i,j) = 255;
        end
        if watermark(i,j) < 70
            watermark (i,j) = 0;
        end
    end
end

colormap(gray(256));
image(watermark);						% shows the image adaptive watermark
title('watermark');
imwrite(marked,gray(256),'watermark.bmp');		% saves the image adaptive watermark as a bmp file


text:
function samples = mcmc_all(Q,lib,nsamp,varargin)
    % MCMC for computing local parameter variance.
    %
    % Input
    %  nsamp: the number of samples
    %
    % Include these string parameters to resample
    %   type
    %   token
    %   relations
    %
    %  By default these variables are not resampled
    %
    % Optional input strings
    %  debug: use full score for debuggin
    %  verbose: track sampling as we go
    %

    narg = numel(varargin);
    verbose = false;
    sample_type = false;
    sample_relations = false;
    sample_token = false;
    debug = false;
    for i=1:narg
        switch varargin{i}
            case 'verbose'
                verbose = true;
            case 'debug'
                debug = true;
            case {'relation','R','relations'}
                sample_relations = true;
            case 'token'
                sample_token = true;
            case 'type'
                sample_type = true;
            otherwise
                error('invalid input parameter');
        end
    end

    if ~(sample_type || sample_token || sample_relations)
       error('empty sampler');
    end

    MH = MCMC(debug);

    M = Q.copy();
    samples = cell(nsamp,1);
    samples_score = nan(nsamp,1);

    % for each sample
    for is = 1:nsamp

        if sample_type
           mcmc_iter_type(MH,M,lib);
        end

        if sample_token
           mcmc_iter_token(MH,M,lib);
        end

        if sample_relations
           mcmc_iter_relations(MH,M,lib);
        end

        samples{is} = M.copy();
        samples_score(is) = scoreMP(M,lib);

        assert(~isinf(samples_score(is)));

        % visualization
        if verbose
            fprintf(1,'Sample %d\n',is);
            figure(1);
            clf
            xlabel('iterations');
            ylabel('log-score');
            plot(1:is,samples_score(1:is));
            drawnow
            pause(0.1);
        end

    end

    if verbose
        MH.acceptance_ratios();
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function x = parse(inStr)

% parse 是一个函数，它接收一个由空格分隔的文本组成的字符串向量，
% 并将各个字符串项解析为一个 n 项矩阵，每个字符串一行。

% x     - the return matrix of strings
% inStr - the blank separated string vector

%%  切割字符串
strLen = size(inStr, 2);
x = blanks(strLen);
wordCount = 1;
last = 0;
for i = 1 : strLen
  if inStr(i) == ' '
    wordCount = wordCount + 1;
    x(wordCount, :) = blanks(strLen);
    last = i;
  else
    x(wordCount, i - last) = inStr(i);
  end
end

repo: Matlab-Machine

path: github-repos/MATLAB/Time9Y-Matlab-Machine/Matlab-Machine-6590d27c49bb2c9415613f132fc2b50fea982c60/009_Optimize time series prediction of BP neural network based on genetic algorithm/goat/parse.m


text:
function test_example_irasa()
% MEM 3gb
% WALLTIME 00:05:00
% DEPENDENCY ft_freqanalysis ft_specest_irasa external/signal/resample.m
% tested on MATLAB_R2022a, macOS_Monterey_12.4, FT_2f387ff (@JM fixed resample scaling issue)

clear all;

% set simulation parameters
A = 1; % scale of 1/f amplitude
C = 1; % 1/f slope
O = 1; % weight of oscillatory components of the simulated data

lf = 1; % lower bound of freq
hf = 300; % higher bound of freq
sl = 600; % spectral lines

fs = 1000; % sampling rate
n = 60000; % time pnts
t = (1:n)/fs; % time axis

% simulate data
for rpt = 1:1

    %     try
    %       % generate pink noise
    %         dspobj = dsp.ColoredNoise('Color', 'pink', 'SamplesPerFrame', length(t));
    %         fn = dspobj()';
    %     catch
    %         % use another method to make pink noise when dsp.ColoredNoise returns licence error
    %         fn = cumsum(randn(1,length(t)));
    %         fn = fn./max(abs(fn)); %%% @JM This scale doesn't seem right?
    %     end

    % another way to simulate pink noise which give users more control over the features of the noise
    freq = linspace(lf, hf, sl); % sampling frequesies
    fn = zeros(size(t));
    for i=1:length(freq) % cummulative sum over freq
        fn = fn + (A * 1/freq(i)^C) * cos(2*pi*freq(i)*t + rand*2*pi); % 1/f amplitude = a*(1/f^c)
    end

    % add a 10Hz and 60 Hz oscillation
    data.trial{1,rpt} = fn + O * cos(2*pi*10*t) + O * cos(2*pi*60*t);
    data.time{1,rpt}  = t;
    data.label{1}     = 'chan';
    data.trialinfo(rpt,1) = rpt;
end

% chunk 2-second segments (gives 1Hz frequency resolution) for long/continous trials
cfg           = [];
cfg.length    = 2; % freqency resolution = 1/2^floor(log2(cfg.length*0.9))
cfg.overlap   = 0.5;
data          = ft_redefinetrial(cfg, data);

% compute the fractal and original spectra
tic
cfg               = [];
cfg.foilim        = [1 200];
cfg.pad           = 'nextpow2';
cfg.method        = 'irasa';
cfg.output        = 'fractal';
fractal = ft_freqanalysis(cfg, data);
cfg.output        = 'original';
original = ft_freqanalysis(cfg, data);
toc % ~28s

% subtract the fractal component from the power spectrum
cfg               = [];
cfg.parameter     = 'powspctrm';
cfg.operation     = 'x2-x1';
oscillatory = ft_math(cfg, fractal, original);

% display the spectra in log-log scale
figure();
hold on;
plot(log(original.freq), log(original.powspctrm),'k');
plot(log(fractal.freq), log(fractal.powspctrm));
plot(log(fractal.freq), log(oscillatory.powspctrm));
xlabel('log-freq'); ylabel('log-power');
legend({'original','fractal','oscillatory'},'location','southwest');

if A~=0 && O==0
    title('pure fractal signal');
elseif A==0 && O~=0
    title('pure oscillatory signal');
elseif A~=0 && O~=0
    title('mixed signal');
end


text:
function [n,ind]=my_histcounts(X,edges,normalization)
persistent old_matlab
if isempty(old_matlab)
    old_matlab=verLessThan('matlab','8.4');
end
if old_matlab
    %We make histc behave like histcounts
    [n,ind]=histc(X,edges);
    ind(ind==length(n))=length(n)-1;
    n(end-1)=n(end-1)+n(end);
    n(end)=[];

    switch normalization
        case 'probability'
            n=n./sum(n);
        case 'count'
        otherwise
            warning('Other types of normalization are not supported on older Matlab versions')
    end
else
    [n, ~, ind]=histcounts(X,edges,'Normalization',normalization);
end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function a = r8row_sort_quick_a ( m, n, a )

%*****************************************************************************80
%
%% R8ROW_SORT_QUICK_A ascending quick sorts an R8ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of A.
%
%    Input, integer N, the number of columns of A.
%
%    Input, real A(M,N), the array to be sorted.
%
%    Output, real A(M,N), the sorted array.
%
  level_max = 30;

  if ( n <= 0 )
    return
  end

  if ( m < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8ROW_SORT_QUICK_A - Fatal error!\n' );
    fprintf ( 1, '  M < 1.\n' );
    error ( 'R8ROW_SORT_QUICK_A - Fatal error!' );
  end

  if ( m == 1 )
    return
  end

  level = 1;
  rsave(level) = m + 1;
  base = 1;
  m_segment = m;

  while ( 1 )
%
%  Partition the segment.
%
    [ a(base:base+m_segment-1,1:n), l_segment, r_segment ] = r8row_part_quick_a ( ...
      m_segment, n, a(base:base+m_segment-1,1:n) );
%
%  If the left segment has more than one element, we need to partition it.
%
    if ( 1 < l_segment )

      if ( level_max < level )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8ROW_SORT_QUICK_A - Fatal error!\n' );
        fprintf ( 1, '  Exceeding recursion maximum of %d\n', level_max );
        error ( 'R8ROW_SORT_QUICK_A - Fatal error!' );
      end

      level = level + 1;
      m_segment = l_segment;
      rsave(level) = r_segment + base - 1;
%
%  The left segment and the middle segment are sorted.
%  Must the right segment be partitioned?
%
    elseif ( r_segment < m_segment )

      m_segment = m_segment + 1 - r_segment;
      base = base + r_segment - 1;
%
%  Otherwise, we back up a level if there is an earlier one.
%
    else

      while ( 1 )

        if ( level <= 1 )
          return
        end

        base = rsave(level);
        m_segment = rsave(level-1) - rsave(level);
        level = level - 1;

        if ( 0 < m_segment )
          break
        end

      end

    end

  end

  return
end


text:
function [ftH2O] = Torr2ftH2O(Torr)
% Convert pressure from torr (same as mmHg) to feet of water at 4 degrees C
% Chad Greene 2012
ftH2O = Torr*0.0446033;