text:
report_this_filefun(mfilename('fullpath'));

bz2 = [];
p1 = []; p1b=[];p99b=[];p50b=[];pmib=[];pmab=[];
p99 = [];
p50 = [];pma=[]; pmi=[];
niv = 1:0.2:4;

tdiff = round((teb - t0b)*365/par1);
ni = str2double(prmptdlg('Number of events in each window?','100'));
na = str2double(prmptdlg('Number of random samples drawn ?','30'));
nr = str2double(prmptdlg('Number of repeats ?','30'));
iwl = str2double(prmptdlg('windowlength ?','2'));
iwl = iwl*365/par1;
zr = (-15:0.1:15)*0;
wai = waitbar(0,' Please Wait ...  ');
set(wai,'NumberTitle','off','Name','Makegrid  -Percent done');;
zr = [];
for k=1:nr
    l = ceil(rand([ni na])*a.Count);
    [cumu, xt] = hist(reshape(a(l,3),[ni,na]),(t0b:par1/365:teb));
    for ti = 1:ni
        mean1 = mean([cumu(1:ti-1,:) ; cumu(ti+iwl+1:ni,:)]);
        mean2 = mean(cumu(ti:ti+iwl,:));
        var1 = cov([cumu(1:ti-1,:) ; cumu(ti+iwl+1:ni,:)]);
        var2 = cov(cumu(ti:ti+iwl,:));
    end     % for i
    as = (mean1 - mean2)./(sqrt(var1/(len-iwl)+var2/iwl));

    pmab = [pmab max(pma) ];
    pma=[];
    waitbar(k/nr);
end

close(wai)
figure
histogram(pmab)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
grid
xlabel('Windowlength in [years]')
ylabel('Range of z')
title(['ni  =  ' num2str(ni) 'events, ' num2str(na) ' random samples'])

matdraw

author: CelsoReyes


text:
classdef StiffnessVoigt2TensorConverter < FourthOrderVoigt2TensorConverter

    properties
    end

    methods (Access = public)

        function obj = StiffnessVoigt2TensorConverter(tensor)
           obj.computeConversion(tensor)
        end

    end

    methods (Access = protected)

        function selectTensorClass(obj)
            obj.tensor = Stiffness3DTensor();
        end

    end

    methods (Access = protected,Static)

        function f = getVoigtFactor(iv,jv)
           f = 1;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function fI = faceinterpolate3(f,node,face,quadOrder)
%% FACEINTERPOLATE3 interpolate to face elements RT0.
%
% uI = faceinterpolate3(u,node,face) interpolates a given function u
% into the lowesr order RT0. The coefficient is given by the face integral
% int_f u*n ds.
%
% Example
%
%
%    [node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
%    maxIt = 3;
%    pde = polynomialdata0;
%    err = zeros(maxIt,1);
%    N = zeros(maxIt,1);
%    for i =1:maxIt
%      [node,elem] = uniformrefine3(node,elem);
%      uI = faceinterpolate3(pde.u,node,elem);
%      err(i) = getL2error3RT0(node,elem,pde.u,uI);
%      N(i) = size(u,1);
%    end
%     figure;
%     showrate(N,err,2,'r-+','||u - u_I||');
%
% See also edgeinterpolate, edgeinterpolate1, edgeinterpolate2
%
% Created by Lin Zhong and Ming Wang at July, 2012.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%%
NF = size(face,1);
v12 = node(face(:,2),:) - node(face(:,1),:);
v13 = node(face(:,3),:) - node(face(:,1),:);
nVec = mycross(v12,v13);

%% assemble the right hand side
if ~exist('quadOrder','var'), quadOrder = 2; end
[lambda,weight] = quadpts(quadOrder);
nQuad = size(lambda,1);
bt = zeros(NF,3);
for p = 1:nQuad
    pxy = lambda(p,1)*node(face(:,1),:) ...
		+ lambda(p,2)*node(face(:,2),:) ...
		+ lambda(p,3)*node(face(:,3),:);
    bt = bt + weight(p)*f(pxy);
end
fI = dot(nVec,bt,2)/2;

repo: Aero-matlab


text:
 function [fun, d1, d2] = ir_poly2_fun(order, varargin)
%function [fun, d1, d2] = ir_poly2_fun(order, [options])
%|
%| create anonymous functions @(x,y) with all the terms in a 2d polynomial
%| up to given order, e.g. [1+0*x x y ...]
%| also return partial derivatives w.r.t. x and y
%|
%| in
%|	order
%|
%| options
%|	maxdegree	maximum of sum of degrees (default: 2*order)
%|	dc		set to 1 to include dc (constant) term (default: 1)
%|
%| out
%|	fun		use fun(x, y) to evaluate the polynomial
%|	d1,d2		likewise, for 1st partial derivatives thereof
%|
%| Copyright 2005-6-18, Jeff Fessler, The University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(order, 'test'), ir_poly2_fun_test, return, end

arg.maxdegree = []; % se below
arg.dc = true;
arg = vararg_pair(arg, varargin);

if isempty(arg.maxdegree), arg.maxdegree = 2 * order; end

str = '';
d1 = '';
d2 = '';
for i2=0:order
	for i1=0:order
		if i1+i2 == 0 && ~arg.dc
			continue
		end
		if i1 + i2 > arg.maxdegree
			continue
		end

		str = [str sprintf(' (x.^%d).*(y.^%d)', i1, i2)];
		if i1 == 0
			s1 = '0*x';
		else
			s1 = sprintf('%d*x.^%d', i1, i1-1);
		end
		if i2 == 0
			s2 = '0*y';
		else
			s2 = sprintf('%d*y.^%d', i2, i2-1);
		end
		d1 = [d1 sprintf(' (%s).*(y.^%d)', s1, i2)];
		d2 = [d2 sprintf(' (x.^%d).*(%s)', i1, s2)];
	end
end
str = ['[ ' str ' ]'];
d1 = ['[ ' d1 ' ]'];
d2 = ['[ ' d2 ' ]'];

% create appropriate anonymous functions
fun = ir_str2func(['@(x,y) ' str]);
d1 = ir_str2func(['@(x,y) ' d1]);
d2 = ir_str2func(['@(x,y) ' d2]);

% ir_poly2_fun_test()
function ir_poly2_fun_test
[fun d1 d2] = ir_poly2_fun(1, 'dc', 0);
x = 1:5; y = 2:6; fun(x,y); d1(x,y); d2(x,y);
[fun d1 d2] = ir_poly2_fun(2, 'maxdegree', 3);
fun(x,y); d1(x,y); d2(x,y);
if im
	pr fun
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function dX = apply_recursive_differentiation(model,x,requested,recursivederivativeprecompute)

global newmodel
persistent dX0
persistent index
persistent mtT
persistent monomTablePattern
persistent ss
persistent lastrequested

% Compute all evaluation-based derivatives df(x)
dxi = [];
dxj = [];
dxs = [];
for i = 1:length(model.evaluation_scheme)
    if strcmp(model.evaluation_scheme{i}.group,'eval')
        for j = model.evaluation_scheme{i}.variables
            if any(requested(model.evalMap{j}.computes))
                k = model.evalMap{j}.variableIndex;
                derivative = model.evalMap{j}.properties.derivative(x(k));
                z{i,j} = derivative;
                if i == 1
                    dxj = [dxj model.evalMap{j}.variableIndex];
                    if length(derivative)>1 && length(model.evalMap{j}.computes)==1
                        dxi = [dxi repmat(model.evalMap{j}.computes,1,length(derivative))];
                    else
                        dxi = [dxi model.evalMap{j}.computes];
                    end
                    dxs = [dxs;derivative(:)];
                end
            end
        end
    end
end

% Apply chain-rule. This code is horrible

if newmodel || ~isequal(requested,lastrequested);
    % Some precalc save over iterations
    ss = any(model.deppattern(requested,model.linearindicies),1);
    monomTablePattern = logical(model.monomtable);
    mtT = model.monomtable';
    [~,dxj] = ismember(dxj,model.linearindicies);
    dxi1 = [dxi model.linearindicies];
    dxj1 = [dxj 1:length(model.linearindicies)];
    dxs1 = [dxs(:)' ones(length(model.linearindicies),1)'];
    dX = sparse(dxi1,dxj1,dxs1,length(model.c),length(model.linearindicies));
    dX0 = dX;
    index = sub2ind([length(model.c),length(model.linearindicies)],dxi,dxj);
    newmodel = 0;
    lastrequested = requested;
else
    dX = dX0;
    dX(index) = dxs;
end

for i = 1:length(model.evaluation_scheme)
	group = model.evaluation_scheme{i}.group;
	if strcmp(group, 'eval')
		if i>1
			for variable = find(ss)
				for j = recursivederivativeprecompute{variable,i}
					 k = model.evalMap{j}.variableIndex;
					 if length(model.evalMap{j}.computes) == 1
						  dX(model.evalMap{j}.computes,variable) = dX(k,variable)'*z{i,j};
					 else
					  val = dX(k,variable).*z{i,j};
					  [idx,~,vdx] = find(val);
					  if ~isempty(vdx)
							dX(model.evalMap{j}.computes(idx),variable) = vdx;
					  end
					 end
				end
			end
		end
	elseif strcmp(group, 'monom')
		computed = model.monomials(model.evaluation_scheme{i}.variables);

		% quadratic and bilinear expressions in the bottom layer of
		% the computational tree can be differentiated very easily
		if i == 1
			 BilinearIndex = find(model.variabletype(computed)==1);
			 Bilinears = computed(BilinearIndex);
			 if ~isempty(Bilinears)
				  x1 = model.BilinearsList(Bilinears,1);
				  x2 = model.BilinearsList(Bilinears,2);
				  [~,x1loc] = ismember(x1,model.linearindicies);
				  [~,x2loc] = ismember(x2,model.linearindicies);
				  dX(sub2ind(size(dX),[Bilinears(:);Bilinears(:)],[x1loc;x2loc]))=x([x2;x1]);
			 end
			 QuadraticIndex = find(model.variabletype(computed)==2);
			 Quadratics = computed(QuadraticIndex);
			 if ~isempty(Quadratics)
				  x1 =  model.QuadraticsList(Quadratics,1);
				  [~,x1loc] = ismember(x1,model.linearindicies);
				  dX(sub2ind(size(dX),Quadratics',x1loc))=2*x(x1);
			 end
			 computed([BilinearIndex QuadraticIndex])=[];
		end

		if i > 1
			 % We might have inner derivatives from earlier
			 isBilinear = model.variabletype==1;
			 if all(isBilinear(computed))
				  % Special case quicker. Only bilinear monomials
				  b1 = model.BilinearsList(computed,1);
				  b2 = model.BilinearsList(computed,2);
				  dp = repmat(x(b2),1,length(model.linearindicies)).*dX(b1,:)+repmat(x(b1),1,length(model.linearindicies)).*dX(b2,:);
				  dX(computed,:) = dp;
			 else
				  compMTP = monomTablePattern(computed,:);
				  for variable = find(ss)
						[dXRows,~,dXVals] = find(dX(:,variable));
						hh = any(compMTP(:,dXRows),2);
						compInd = computed(hh);
						newIndices = compInd(requested(compInd)); % all values for loop
						% Create a working copy of dX(:,variable). Real size isn't
						% length(dXRows)+length(newIndices) but length(unique([dXRows;newIndices']).
						% But it's more efficient to allocate some more memory than calculate the true length.
						dx = sparse(dXRows,1,dXVals, size(dX,1),1, length(dXRows)+length(newIndices));
						for j = newIndices(:)'
							if isBilinear(j)
								 b1 = model.BilinearsList(j,1);
								 b2 = model.BilinearsList(j,2);
								 dp = x(b2)*full(dx(b1))+x(b1)*full(dx(b2));
							else
								 dp = 0;
								 monomsj = mtT(:,j);
								 active = find(monomsj & dx);
								 for k = active(:)'
									monoms = monomsj;
									[s,~,v] = find(monoms);
									k2s = s==k;
									r = v(k2s);
									v(k2s) = r-1;
									ztemp = prod(x(s).^v);
									aux = full(dx(k));
									dp = dp + r*aux*ztemp;
								 end
							end
							dx(j) = real(dp);
						end
						dX(:,variable) = dx(:);
				  end
			 end
		else
			% We are at the bottom layer, so no inner derivatives
			compMTP = monomTablePattern(computed,:);
			for variable = find(ss)
				dx = dX(:,variable);
				fdX = find(dx); % logical(dx); for sparse-matrices, FIND is about as fast.
				hh = any(compMTP(:,fdX),2); %#ok<FNDSB>
				compInd = computed(hh);
				for j = compInd(requested(compInd))
					k = model.linearindicies(variable);
					monoms = mtT(:,j);

					[s,~,v] = find(monoms);
					k2s = s==k;
					r = v(k2s);
					v(k2s) = r-1;
					dp = r*prod(x(s).^v);
					dX(j,variable) = real(dp);
				end
			end
		end
	end
end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Functions enabling use of multiple colormaps per figure.
%
% version 2.3, 3/2007
%
%   freezeColors    - Lock colors of plot, enabling multiple colormaps per figure.
%   unfreezeColors  - Restore colors of a plot to original indexed color.
%
% demo/html
%   freezeColors_pub.html   - Published demonstration.
%
% test
%   test_main       - Test function.

% AUTHOR
% John Iversen, 2005-10
% john_iversen@post.harvard.edu
%
% Free for any use, so long as AUTHOR information remains with code.
%
% HISTORY
%
% JRI 6/2005  (v1)
% JRI 9/2006  (v2.1) now operates on all objects with CData (not just images as before)
% JRI 11/2006 handle NaNs in images/surfaces (was not posted on file exchange, included in v2.3)
% JRI 3/2007  (v2.3) Fix bug in handling of caxis that made colorbar scales on frozen plots incorrect.
% JRI 4/2007  Add 'published' html documentation.
% JRI 9/2010  Changed documentation for colorbars

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
classdef RigidJoint < CoordinateFrame
    % A coordinate frame that is rigidly attached to the origin a rigid
    % joint and the child link of the joint. In addition, this class stores
    % the type, rotation axis, child link, and limit of the rigid joints.
    %
    % Copyright (c) 2016-2017, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    properties (SetAccess=protected, GetAccess=public)
        % The type of the rigid joint
        %
        % @note valid joint types include:
        % revolute, prismatic, continuous, fixed.
        %
        % @type char
        Type

        % The rotation axis vector of the rigid joint
        %
        % @note it must be a 3x1 vector.
        %
        % @type rowvec
        Axis

        % The name of the child link
        %
        % @type char
        Child

        % The name of the parent link
        %
        % @type char
        Parent

        % A structure data stores the physical limits of the rigid joints
        %
        % Required fields of limit:
        % effort: the maximum actuator effort @type double
        % lower: the minimum allowable displacement @type double
        % upper: the maximum allowable displacement @type double
        % velocity: the maximum allowable velocity @type double
        %
        % @type struct
        Limit

        % A structure contains information about the actuator for the joint
        %
        % An empty property indicates the joint is not actuated.
        %
        % @type struct
        Actuator

    end

    properties (Hidden, SetAccess=protected, GetAccess=public)

        % The indices of kinematic chains from the base to the current
        % joint
        %
        % @type rowvec
        ChainIndices

        % The twist from the base frame
        %
        % @type rowvec
        Twist

        % The twist pairs of all precedent joints (coordinate frame)
        %
        % @type SymExpression
        TwistPairs
    end
    methods

        function obj = RigidJoint(varargin)
            % The class constructor function
            %
            % Parameters:
            %  varargin: variable nama-value pair input arguments, in detail:
            %    Name: the name of the frame @type char
            %    Reference: the reference frame @type CoordinateFrame
            %    Offset: the offset of the origin @type rowvec
            %    R: the rotation matrix or the Euler angles @type rowvec
            %    Axis: the rotation axis vector @type rowvec
            %    Type: the type of the joints @type char
            %    Child: the child link frame @type RigidLink
            %    Limit: the joint physical limits @type struct

            % consruct the superclass object first
            obj = obj@CoordinateFrame(varargin{:});
            if nargin == 0
                return;
            end
            argin = struct(varargin{:});

            % validate and assign the joint type
            if isfield(argin, 'Type') && ~isempty(argin.Type)
                obj = obj.setType(argin.Type);
            else
                error('The joint type is not defined.');
            end

            % validate and assign the joint axis
            if isfield(argin, 'Axis') && ~isempty(argin.Axis)
                obj = obj.setAxis(argin.Axis);
            elseif strcmp(obj.Type, 'fixed')
                obj = obj.setAxis([0,0,1]);
            else
                error('The joint rotation axis is not defined.');
            end

            % validate and assign the child link
            if isfield(argin, 'Child') && ~isempty(argin.Child)
                obj = obj.setChild(argin.Child);
            else
                error('The child link is not defined.');
            end

            % validate and assign the parent link
            if isfield(argin, 'Parent') && ~isempty(argin.Parent)
                obj = obj.setParent(argin.Parent);
            else
                error('The parent link is not defined.');
            end

            % validate and assign the physical limits
            if isfield(argin, 'Limit') && ~isempty(argin.Limit)
                obj = obj.setLimit(argin.Limit);
            else
                warning('The joint limits are not defined. Using default values.');
                default_limit = struct(...
                    'effort',0,...
                    'lower',-inf,...
                    'upper',inf,...
                    'velocity',inf);
                obj = obj.setLimit(default_limit);
            end

            % validate and assign the actuator info
            if isfield(argin, 'Actuator') && ~isempty(argin.Actuator)
                obj = obj.setActuator(argin.Actuator);
            end
        end

    end

    %% methods defined in external files
    methods
        obj = setAxis(obj, axis);

        obj = setType(obj, type);

        obj = setParent(obj, parent);

        obj = setChild(obj, child);

        obj = setLimit(obj, varargin);

        obj = setActuator(obj, varargin);

        obj = setChainIndices(obj, indices);

        xi = getTwist(obj);

        obj = computeTwist(obj);

        obj = setTwistPairs(obj, dofs, q);
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
%HCLUST hierarchical clustering
%
%  [LABELS, DENDROGRAM] = HCLUST(D,TYPE,K)
%   DENDROGRAM = HCLUST(D,TYPE)
%
% INPUT
%  D     dissimilarity matrix
%  TYPE  string name of clustering criterion (optional)
%        - 's' or 'single'   : single linkage (default)
%        - 'c' or 'complete' : complete linkage
%        - 'a' or 'average'  : average linkage (weighted over cluster sizes)
%  K     number of clusters (optional)
%
% OUTPUT
%  LABELS       vector with labels
%  DENDROGRAM   matrix with dendrogram
%
% DESCRIPTION
% Computation of cluster labels and a dendrogram between the clusters for
% objects with a given distance matrix D. K is the desired number of
% clusters. The dendrogram is a 2*K matrix. The first row yields all
% cluster sizes. The second row is the cluster level on which the set of
% clusters starting at that position is merged with the set of clusters
% just above it in the dendrogram. A dendrogram may be plotted by PLOTDG.
%
%   DENDROGRAM = HCLUST(D,TYPE)
%
% As in this case no clustering level is supplied, just the entire
% dendrogram is returned. The first row now contains the object indices.
%
% EXAMPLE
% a = gendats([25 25],20,5);    % 50 points in 20 dimensional feature space
% d = sqrt(distm(a));           % Euclidean distances
% dendg = hclust(d,'complete'); % dendrogram
% plotdg(dendg)
% lab = hclust(d,'complete',2); % labels
% confmat(lab,getlabels(a));    % confusion matrix
%
% SEE ALSO (PRTools Guide)
% PLOTDG, PRKMEANS, KCENTRES, MODESEEK, EMCLUST

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

% $Id: hclust.m,v 1.3 2008/02/14 11:54:43 duin Exp $

function [labels, dendrogram] = hclust(D,type,k)

	if nargin < 3, k = []; end
	if nargin < 2, type = 's'; end
	[m,m1] = size(D);
	if m ~= m1
		error('Input matrix should be square')
	end
	D = D + diag(inf*ones(1,m));     % set diagonal at infinity.
	W = [1:m+1];               % starting points of clusters in linear object set.
	V = [1:m+2];               % positions of objects in final linear object set.
	F = inf * ones(1,m+1);     % distance of next cluster to previous cluster to be stored at first point of second cluster
	Z = ones(1,m); % number of samples in a cluster (only for average linkage)

  t = sprintf('Analysing %i cluster levels: ',m);
  prwaitbar(m,t);
	for n = 1:m-1
    prwaitbar(m,n,[t int2str(n)]);
		% find minimum distance D(i,j) i<j
		[di,I] = min(D); [dj,j] = min(di); i = I(j);
		if i > j, j1 = j; j = i; i = j1; end
		% combine clusters i,j
		switch type
		case {'s','single'}
			D(i,:) = min(D(i,:),D(j,:));
		case {'c','complete'}
			D(i,:) = max(D(i,:),D(j,:));
		case {'a','average'}
			D(i,:) = (Z(i)*D(i,:) + Z(j)*D(j,:))/(Z(i)+Z(j));
			Z(i:j-1) = [Z(i)+Z(j),Z(i+1:j-1)]; Z(j) = [];
		otherwise
			error('Unknown clustertype desired')
		end
		D(:,i) = D(i,:)';
		D(i,i) = inf;
		D(j,:) = []; D(:,j) = [];
		% store cluster distance
		F(V(j)) = dj;
		% move second cluster in linear ordering right after first cluster
		IV = [1:V(i+1)-1,V(j):V(j+1)-1,V(i+1):V(j)-1,V(j+1):m+1];
		W = W(IV); F = F(IV);
		% keep track of object positions and cluster distances
		V = [V(1:i),V(i+1:j) + V(j+1) - V(j),V(j+2:m-n+3)];
  end
  prwaitbar(0);
	if ~isempty(k) | nargout == 2
		if isempty(k), k = m; end
		labels = zeros(1,m);
		[S,J] = sort(-F);      % find cluster level
		I = sort(J(1:k+1));    % find all indices where cluster starts
		for i = 1:k            % find for all objects cluster labels
			labels(W(I(i):I(i+1)-1)) = i * ones(1,I(i+1)-I(i));
		end                    % compute dendrogram
		dendrogram = [I(2:k+1) - I(1:k); F(I(1:k))];
		labels = labels';
	else
		labels = [W(1:m);F(1:m)]; % full dendrogram
	end
return

author: marianux

repo: ecg-kit

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


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%

randidx = randperm(size(X , 1));
centroids = X(randidx(1 : K) ,:);

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function tests = test_ft_preproc_standardize

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_standardize

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

state = [];

result = [];
result{end+1} = ft_preproc_standardize(dat, 1, 1000, state);
result{end+1} = ft_preproc_standardize(dat, 1, 900, state);
result{end+1} = ft_preproc_standardize(dat, 101, 1000, state);
result{end+1} = ft_preproc_standardize(dat, 101, 900, state);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mining_queue = esvm_initialize_mining_queue(imageset, ordering)
%Initialize the mining queue with ordering (random by default)
%function mining_queue = esvm_initialize_mining_queue(imageset, ordering)
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if ~exist('ordering','var')
  fprintf(1,'Randomizing mining queue\n');
  myRandomize;
  ordering = randperm(length(imageset));
end

mining_queue = cell(0,1);
for zzz = 1:length(ordering)
  mining_queue{zzz}.index = ordering(zzz);
  mining_queue{zzz}.num_visited = 0;
end

author: quantombone

repo: exemplarsvm


text:
function [result] = sv_bval(params)
% function [params] = pf_bvalues(params)
% --------------------------------------
% Calculation of b-values for the probabilistic forecast test.
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bMap               Calculate a map (=1) or a cross-section (=0)
%   params.bNumber            Use constant number (=1) or constant radius (=0)
%   params.nNumberEvents      Number of earthquakes if bNumber == 1
%   params.fRadius            Radius of gridnode if bNumber == 0
%   params.nMinimumNumber     Minimum number of earthquakes per node
%   params.fSplitTime         Time at which the catalog will be split
%   params.bLearningPeriod    Fix duration of the learning period (1 = fix, 0 = use catalog from start)
%   params.fLearningPeriod    Duration of the learning period
%   params.bForecastPeriod    Fix duration of the forecasting period (1 = fix, 0 = use catalog till end)
%   params.fForecastPeriod    Duration of the forecasting period
%
% Output parameters:
%   Same as input parameters including
%   result.mValueGrid         Matrix of calculated b-values
%   result.vcsGridNames       Strings describing the values in mValueGrid
%
% Danijel Schorlemmer
% June 4, 2001

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Init result matrix
result = params;
result.mValueGrid = [];

for nNode_ = 1:length(params.mPolygon(:,1))
  % Create node catalog
  mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
  % Calculate the b-values
  [vDummy, fBValueObserved, fStdDevObserved, vDummy] =  bmemag(mNodeCatalog_);
  result.mValueGrid = [result.mValueGrid; fBValueObserved fStdDevObserved];
end; % for nNode
% Create the description-strings for the output-window
result.vcsGridNames = cellstr(char('b-value', 'Std. dev. of b-value'));


text:
%% Script to test the optical flow estimation
clear all
close all

img1f = '/datos/research/deblurring/data/jue/videoDeblur/snow/input/0001.png';
img2f = '/datos/research/deblurring/data/jue/videoDeblur/snow/input/0002.png';

%img1 = single(rgb2gray(imread(img1f)));
%img2 = single(rgb2gray(imread(img2f)));

img1 = double((imread(img1f)));
img2 = double((imread(img2f)));

img1 = .299*img1(:,:,1) + .587*img1(:,:,2) + .114*img1(:,:,3);
img2 = .299*img2(:,:,1) + .587*img2(:,:,2) + .114*img2(:,:,3);

img1 = single(floor(img1));
img2 = single(floor(img2));

figure(1)
imagesc(img2)

%%
tau=0.25;
lambda = 0.1;
theta=0.3;
nscales=4;
zfactor=0.5;
nwarps=5;
epsilon=0.01;
verbose=1;

[m,n,~] = size(img1);
tagstruct.ImageLength = m;
tagstruct.ImageWidth = n;
tagstruct.Compression = Tiff.Compression.None;
tagstruct.SampleFormat = Tiff.SampleFormat.IEEEFP;
tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
tagstruct.BitsPerSample =  32;                        % float data
tagstruct.SamplesPerPixel = 2; %u and v (2 channels)
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;

flowF = tvl1flow(img1,img2,tau,lambda,theta,nscales,zfactor,nwarps, epsilon,verbose);

%Save into a TIFF file
t = Tiff('aux1.tiff', 'w');
t.setTag(tagstruct);
t.write(single(flowF));
t.close();

%% OLD

setenv('DYLD_LIBRARY_PATH', '/usr/local/lib'),
tvl1flowS='/datos/research/deblurring/code/heat/videodeblur/tvl1flow_3_mex/tvl1flow';

flowB = 'aux.tiff';
system(sprintf('%s %s %s %s 0 0.25 0.1 0.3 4 0.5 5 0.01 1',tvl1flowS,img1f,img2f,flowB)); % flowB.tiff

%%
aux1 = imread('aux1.tiff');
aux = imread('aux.tiff');

figure(1)
imagesc(flowF(:,:,1)-aux(:,:,1))

figure(2)
imagesc(flowF(:,:,2)-aux(:,:,2))

author: shuochsu

repo: DeepVideoDeblurring


text:
function db_open(user, password, server, port, database)

% DB_OPEN opens the connection to the database
%
% Use as
%    db_open
%    db_open(user, password, server, port, database)
%    db_open('mysql://<user>:<password>@<host>:<port>')

% Copyright (C) 2007, Robert Oostenveld
%
% $Log: db_open.m,v $
% Revision 1.2  2007/11/27 14:38:49  roboos
% added query "use database"
%
% Revision 1.1  2007/11/07 10:50:46  roboos
% created helper functions for easy access to a MySQL database table using a structure for representing the data
%

% persistent variables should be defined at the beginning of the function
persistent prev_filename

switch nargin
  case 0
    % use the default settings and combine them in a string of the form mysql://<user>:<password>@<host>:<port>
    user     = 'fieldtrip';
    password = 'fieldtrip';
    server   = 'odin';
    port     = 3306;
    filename = sprintf('mysql://<%s>:<%s>@<%s>:<%d>', user, password, server, port);
  case 1
    % the database settings are specified as mysql://<user>:<password>@<host>:<port>
    filename = user;
    [user, password, server, port] = filetype_check_uri(filename);
  case 3
    % the port was not specified, combine the settings in a string of the form mysql://<user>:<password>@<host>
    filename = sprintf('mysql://<%s>:<%s>@<%s>', user, password, server);
  case 4
    % combine the settings in a string of the form mysql://<user>:<password>@<host>:<port>
    filename = sprintf('mysql://<%s>:<%s>@<%s>:<%d>', user, password, server, port);
  case 5
    % combine the settings in a string of the form mysql://<user>:<password>@<host>:<port>
    filename = sprintf('mysql://<%s>:<%s>@<%s>:<%d>/<%s>', user, password, server, port, database);
  otherwise
    ft_error('incorrect input arguments');
end

if ~strcmp(filename, prev_filename)
  % close the database, it will be reopened furter down in the code
  mysql('close');
  prev_filename = [];
end

if ~strcmp(filename, prev_filename)
  % open the database
  if ~isempty(port)
    server = sprintf('%s:%d', server, port);
  end
  try
    mysql('open', server, user, password);
    if nargin==5
      mysql('use', database);
    end
    prev_filename = filename;
  catch
    prev_filename = [];
    ft_warning(lasterr);
  end
end


text:
%%test network
cam1_size=size(test_data_cam1);
cam2_size=size(test_data_cam2);
cam1_feature1=[];
cam1_feature2=[];
cam1_feature3=[];
cam1_score=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam1_size(1)/param.test_batch_size)-1
    for n=1:param.test_batch_size
        im_data=imresize((reshape(single(test_data_cam1(m*param.test_batch_size+n,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam1_feature1=[cam1_feature1;(squeeze(net.blobs('data1_p_1').get_data))'];
    cam1_feature2=[cam1_feature2;(squeeze(net.blobs('data1_p_2').get_data))'];
    cam1_feature3=[cam1_feature3;(squeeze(net.blobs('data1_p_3').get_data))'];
    cam1_score=[cam1_score;(squeeze(net.blobs('sig').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size+1:cam1_size(1)
        count=m-(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imresize((reshape(single(test_data_cam1(m,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam1_size(1)-(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('data1_p_1').get_data))';
result2=(squeeze(net.blobs('data1_p_2').get_data))';
result3=(squeeze(net.blobs('data1_p_3').get_data))';
score=(squeeze(net.blobs('sig').get_data))';
cam1_feature1=[cam1_feature1;result1(1:index,:)];
cam1_feature2=[cam1_feature2;result2(1:index,:)];
cam1_feature3=[cam1_feature3;result3(1:index,:)];
cam1_score=[cam1_score;score(1:index,:)];
prob_feature=[];
name_cam1={};
for m=1:param.train_person_num
    test_label_index=find(label_test_cam1==(m+param.train_person_num));
    a=bsxfun(@times,cam1_feature1(test_label_index,:),cam1_score(test_label_index,1));
    b=bsxfun(@times,cam1_feature2(test_label_index,:),cam1_score(test_label_index,2));
    c=bsxfun(@times,cam1_feature3(test_label_index,:),cam1_score(test_label_index,3));
    prob_feature=[prob_feature;sum(a/sum(cam1_score(test_label_index,1))),sum(b/sum(cam1_score(test_label_index,2))),sum(c/sum(cam1_score(test_label_index,3)))];
%     prob_feature=[prob_feature;sum(a/sum(cam1_score(test_label_index,1)))];
%     prob_feature=[prob_feature;sum(b/sum(cam1_score(test_label_index,2)))];
%     prob_feature=[prob_feature;sum(c/sum(cam1_score(test_label_index,3)))];
name_cam1=[name_cam1,test_image_name_cam1{test_label_index(1)}];
end
%extract feature cam2
cam2_feature1=[];
cam2_feature2=[];
cam2_feature3=[];
cam2_score=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam2_size(1)/param.test_batch_size)-1
    for n=1:param.test_batch_size
        im_data=imresize((reshape(single(test_data_cam2(m*param.test_batch_size+n,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam2_feature1=[cam2_feature1;(squeeze(net.blobs('data1_p_1').get_data))'];
    cam2_feature2=[cam2_feature2;(squeeze(net.blobs('data1_p_2').get_data))'];
    cam2_feature3=[cam2_feature3;(squeeze(net.blobs('data1_p_3').get_data))'];
    cam2_score=[cam2_score;(squeeze(net.blobs('sig').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size+1:cam2_size(1)
        count=m-(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imresize((reshape(single(test_data_cam2(m,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam2_size(1)-(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('data1_p_1').get_data))';
result2=(squeeze(net.blobs('data1_p_2').get_data))';
result3=(squeeze(net.blobs('data1_p_3').get_data))';
score=(squeeze(net.blobs('sig').get_data))';
cam2_feature1=[cam2_feature1;result1(1:index,:)];
cam2_feature2=[cam2_feature2;result2(1:index,:)];
cam2_feature3=[cam2_feature3;result3(1:index,:)];
cam2_score=[cam2_score;score(1:index,:)];
gallery_feature=[];
name_cam2={};
for m=1:param.train_person_num
    test_label_index=find(label_test_cam2==(m+param.train_person_num));
    a=bsxfun(@times,cam2_feature1(test_label_index,:),cam2_score(test_label_index,1));
    b=bsxfun(@times,cam2_feature2(test_label_index,:),cam2_score(test_label_index,2));
    c=bsxfun(@times,cam2_feature3(test_label_index,:),cam2_score(test_label_index,3));
    gallery_feature=[gallery_feature;sum(a/sum(cam2_score(test_label_index,1))),sum(b/sum(cam2_score(test_label_index,2))),sum(c/sum(cam2_score(test_label_index,3)))];
%     gallery_feature=[gallery_feature;sum(a/sum(cam2_score(test_label_index,1)))];
%     gallery_feature=[gallery_feature;sum(b/sum(cam2_score(test_label_index,2)))];
%     gallery_feature=[gallery_feature;sum(c/sum(cam2_score(test_label_index,3)))];
name_cam2=[name_cam2,test_image_name_cam2(test_label_index(1))];
end
% cal cmc
prob_norm=bsxfun(@rdivide,prob_feature,sum(abs(prob_feature).^2,2).^(1/2));
gallery_norm=bsxfun(@rdivide,gallery_feature,sum(abs(gallery_feature).^2,2).^(1/2));
% [~,it]=sort(prob_norm*gallery_norm',2,'descend');
score_matrix=prob_norm*gallery_norm';
rank1_hit=0;
rank5_hit=0;
rank10_hit=0;
rank20_hit=0;
problen=size(prob_feature,1);
for m=1:problen
    [~,location]=sort(score_matrix(m,:),2,'descend');
    if find(location(1)==m)
        rank1_hit=rank1_hit+1;
    end
    if find(location(1:5)==m)
        rank5_hit=rank5_hit+1;
    end
    if find(location(1:10)==m)
        rank10_hit=rank10_hit+1;
    end
    if find(location(1:20)==m)
        rank20_hit=rank20_hit+1;
    end
end
rank_acc=[rank1_hit/problen,rank5_hit/problen,rank10_hit/problen,rank20_hit/problen];
fin=fopen(param.result_save_file,'a');
fprintf(fin,'2below split_index:%d,iter:%d, rank1:%f,rank5:%f,rank10:%f,rank20:%f\n',split_index,iter,rank_acc(1),rank_acc(2),rank_acc(3),rank_acc(4));
fprintf('split_index:%d,iter:%d, rank1:%f,rank5:%f,rank10:%f,rank20:%f\n',split_index,iter,rank_acc(1),rank_acc(2),rank_acc(3),rank_acc(4));
fclose(fin);

author: liuyuisanai

repo: Quality-Aware-Network


text:
function arrivals_per_event(cobj)
    N = numel(cobj.arrivals);
    na = [];
    if N>1
        for c=1:N
            na = [na numel(cobj.arrivals{c}.time)];
        end
        m = max(na);
        for v=1:m
            disp(sprintf('arrivals = %d, count = %d', v, sum(na==v) ));
        end
    end

end

author: geoscience-community-codes

repo: GISMO


text:
% :Usage:
% ::
%
%     OUT = mean_image(input file names, output file name, [weights], [string command args]);
%
% Creates a weighted mean of images based on weights for each image
%
% :Inputs:
%
%   **'normlike':**
%        do in style of mean_warped_image, for norm toolbox;
%        ignore other string arguments
%
%
%   **'reweight':**
%        weight by distance from median in 2nd round of averaging
%
%   **'plot':**
%        plot output
%
%   **'sharpen':**
%        segments and smooths within tissue classes
%
% :Examples:
% ::
%
%    w = filenames('w*T1.img', 'char')
%    w = mean_image(w, 'final_mean.img', [], 'normlike'); % trimmed, weighted best 50%
%
%    worig = filenames('hr*/wT1.img', 'char')           % no weights
%    w = mean_image(worig, 'orig_mean_noweights.img', []);
%
%    mean_image(DB.PP, 'Activation.img', DB.studyweight);
%
%    mean_image(P, 'mean_wmeanw1007.img', ones(size(P, 1), 1), 'sharpen', 'plot', 'reweight');
%
%    P = filenames('w*mr.img', 1);
%
%    % Enter empty inputs for default settings.
%    OUT.P = filenames
%    OUT.dist

function OUT = mean_image(VP, Pout, w, varargin)
    % --------------------------------------
    % inputs
    % --------------------------------------
    dosharp = 0;
    dodist = 0;
    doplot = 0;
    doreweight = 0;
    donormlike = 0;

    dozscore = 0;

    for i = 1:length(varargin)
        if ischar(varargin{i})
            if strcmp(varargin{i}, 'sharpen'), dosharp = 1; end
            %if strcmp(varargin{i}, 'dist'), dodist = 1; end
            if strcmp(varargin{i}, 'plot'), doplot = 1; end
            if strcmp(varargin{i}, 'reweight'), doreweight = 1; end
            if strcmp(varargin{i}, 'normlike'), donormlike = 1; end

            if strcmp(varargin{i}, 'zscore'), dozscore = 1; end
        end
    end

    if isempty(VP)
        VP = spm_vol(spm_get(1, '*', 'Select images'));
    else
        if(iscellstr(VP)),
            VP = char(VP);
        end
        P = VP;
        VP = spm_vol(VP);
    end
    n = length(VP);

    if ~exist('Pout', 'var') || isempty(Pout), Pout = 'mean.nii'; end
    if ~exist('w', 'var') || isempty(w), w = ones(n, 1); end

    OUT.P = VP;

    % --------------------------------------
    % do mean in style of mean_warped_image
    % --------------------------------------
    % output is weights
    if donormlike
        OUT = do_mean_like_mean_warped_image(P, Pout);
        return
    end

    % --------------------------------------
    % make mean image
    % --------------------------------------

    status_string = sprintf('%3d', 0);
    fprintf(['Percent done: ' status_string]);

    switch spm('Ver')
        case 'SPM2'
            % spm_defaults is a script
            disp('WARNING: spm defaults not set for spm2. Make sure your defaults are set correctly');

        case 'SPM5'
            % spm_defaults is a function
            spm_defaults()
        case 'SPM8'
            spm_defaults()
    end

    % initialize output
    V = VP(1);
    V.fname = Pout;
    V.n(1) = 1;    % write to first volume

    dat = zeros(V.dim(1:3));

    % read and weight

    for i = 1:n
        v = iimg_read_vols(VP(i));

        if dozscore
            v = center_scale(v);
        end

        dat = dat + v .* w(i);

        erase_string(status_string);
        status_string = sprintf('%3d', round(i/n*100));
        fprintf(status_string);
    end

    % divide by N
    dat = dat ./ n;

    fprintf('\n');

    % write output
    %strrep required to make file extension correct
    V.fname = strrep(V.fname, '.nii,1', '.nii');
    V.fname = strrep(V.fname, '.img,1', '.img');
    spm_write_vol(V, dat);

    if doplot
        z = round(V.dim(3)./2);
        f1 = figure;
        imagesc(dat(:,:,z));
        axis image;
        axis off;
        title('Mean');  drawnow
    end

    % --------------------------------------
    % Reweight
    % --------------------------------------

    if doreweight
        fprintf(1, ' Reweighting.000');

        for i = 1:n
            fitness(i) = norm_eval_fitness(tempdat, normdat, 0);
        end

        imgMedian = zeros(n, 1);
        imgMean = imgMedian;
        imgStd = imgMedian;

        medianDev = zeros(n, 5);
        newdat = zeros(V.dim(1:3));

        % subtract median from mean image
        [dat, datMedian] = subtract_median(dat);

        % get medians and subtract
        for i = 1:n

            fprintf(1, '\b\b\b%03d', i);

            v = spm_read_vols(spm_vol(deblank(P(i,:))));

            [v, imgMedian(i)] = subtract_median(v);

            %imgMean(i) = mean(vi);
            %imgStd(i) = std(vi);

            % distance from median-centered old mean image
            vi = v - dat;  vi = abs(vi(:));
            vi(isnan(vi) | vi==0) = [];
            if isempty(vi), warning('Image is EXACTLY the median image.'); end

            % get median deviation, and 5th and 95th prctiles
            medianDev(i,:) = prctile(vi, [50 5 20 80 95]);

            % get new weights based on distance from mean image
            w(i) = datMedian./medianDev(i, 1);

            newdat = newdat + v .* w(i);

        end

        newdat = newdat .* sum(w);          % normalize to keep in same scale

        whnull = isnan(newdat) | newdat==0;
        %newdat = newdat + mean(imgMedian);  % add median back in to average
        newdat(whnull) = 0;

        dat = newdat;
        OUT.medianDev = medianDev;
        OUT.prctiles = [50 5 20 80 95];

        if doplot
            tor_fig;
            plot(medianDev(:,1), 'o-');
            ylabel('Med. Abs. Dev from median image');
            xlabel('Image number'); grid on;

            figure(f1); subplot(4, 1, 2);
            imagesc(dat(:,:,z));
            axis image;
            axis off;
            title('Reweighted');drawnow
        end

        % write output
        spm_write_vol(V, dat);
    end

    % --------------------------------------
    % Sharpen
    % --------------------------------------

    if dosharp
        fprintf(1, ' Sharpening: segmenting.');

        defaults.segment.estimate.reg = .1;% was .01;
        VO = spm_segment(Pout, which('T1.mnc'), defaults.segment);

        spm_write_vol(VO(1), VO(1).dat);
        spm_write_vol(VO(2), VO(2).dat);
        spm_write_vol(VO(3), VO(3).dat);

        fprintf(1, ' Lightening white matter.');

        thr = prctile(VO(2).dat(:), 90);
        whiteMask = double(VO(2).dat > thr);

        mystd = .2 .* nanstd(dat(:));
        dat = dat + whiteMask .* mystd;

        % write output
        name = [Pout(1:end-4) '_sharp.img'];
        V.fname = name;
        spm_write_vol(V, dat);

        % make brain mask
        brainMask = double(VO(1).dat>0 | VO(2).dat>0);

        name = [Pout(1:end-4) '_brain.img'];
        V.fname = name;
        spm_write_vol(V, brainMask);

        ovl = dat .* brainMask;
        name = [Pout(1:end-4) '_overlay.img'];
        V.fname = name;
        spm_write_vol(V, ovl);

        if doplot
            figure(f1); subplot(4, 1, 3);
            imagesc(dat(:,:,z));
            axis off;
            axis image;
            title('Sharpened');drawnow

            subplot(4, 1, 4);
            imagesc(ovl(:,:,z));
            axis off;
            axis image;
            title('Overlay');drawnow
        end
    end
end

function [dat, datMedian] = subtract_median(dat)
    vi = dat(:);
    vi(isnan(vi) | vi==0) = [];
    datMedian = median(vi);

    whnull = isnan(dat) | dat==0;

    dat = dat - datMedian;
    dat(whnull) = 0;
end

function v = center_scale(v)
    vi = v(:);
    vi(isnan(vi) | vi==0) = [];
    v = (v - mean(vi)) ./ std(vi);
end

function weights = do_mean_like_mean_warped_image(P, outname)
    spm_defaults;
    defaults.analyze.flip = 0;
    [volTemplate, normdat] = iimg_read_img(P);

    t1 = clock;
    n = size(P, 1);

    % Provisional mean
    % ------------------------------------
    fprintf(1, 'Initial. ');

    normdat(isnan(normdat)) = 0;

    % convert in-mask voxels to mean = gm, so averaging is not distorted by
    % overall intensity differences among images
    % and fitness weights determine relative contribution of images
    means = mean(normdat);      % mean of each image
    gm = mean(means);           % global mean
    scalef = gm ./ means;       % scaling factor to norm each image to grand mean

    for i = 1:n
        normdat(:,i) = normdat(:,i) .* scalef(i);
    end

    meandat = mean(normdat, 2);
    vox = size(meandat, 1);

    fprintf(1, 'Weighted mean of best 50%%. ');
    % Weights based on closeness to mean
    % ------------------------------------
    for i = 1:n
        fitness(i) = norm_eval_fitness(meandat, normdat(:,i), 0);
    end

    mfit = nanmedian(fitness);
    weights = fitness;
    weights(weights < mfit) = 0;

    weights = weights ./ sum(weights);

    if all(weights == 0)
        warning('All weights are zero!')
        weights = ones(1, n);
    end

    % Weighted mean
    % ------------------------------------
    meandat = zeros(vox, 1);
    for i = 1:n

        if weights(i) > 0   % just to save time
            meandat = meandat + weights(i) .* normdat(:,i);
        end

    end

    fprintf(1, 'Done: Writing\n');
    % Write output image
    % This will be the new normalization template
    % ------------------------------------
    meandata = iimg_reconstruct_3dvol(meandat, volTemplate, 'outname', outname);

    fprintf(1, 'Mean warped image completed: %3.0f s\n', etime(clock, t1));

    fprintf(1, '-----------------------------\n');
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function qout=qnorm(qin)
% QNORM(Q) normalizes quaternions.
%     Works on vectors of quaternions too.  If input is a vector of four
%     quaternions, QNORM will determine whether the quaternions are row or
%     column vectors according to ISQ.
%
% See also ISQ.

% Release: $Name: quaternions-1_2_2 $
% $Revision: 1.9 $
% $Date: 2001/05/01 20:20:31 $

% Copyright (C) 2001, Jay A. St. Pierre.  All rights reserved.

if nargin~=1
  error('qnorm() requires one input argument');
else
  qtype = isq(qin);
  if ( qtype == 0 )
    error(['Invalid input: must be a quaternion or a vector of' ...
          ' quaternions'])
  elseif ( qtype==3 )
    warning(['Component quaternion shape indeterminate... assuming row' ...
             ' vectors'])
  end
end

% Make sure qin is a column of quaternions
if( qtype == 1 )
  qin=qin.';
end

% Find the magnitude of each quaternion
qmag=sqrt(sum(qin.^2,2));

% Make qmag the same size a q
qmag=[qmag qmag qmag qmag];

% Divide each element of q by appropriate qmag
qout=qin./qmag;

% Make sure output is same shape as input
if( qtype == 1 )
  qout=qout.';
end

author: christianwengert


text:
%% (Internal) Calculate an estimate of the SNR of ECG signal, considering signal the ensemble average, and noise each realization minus the ensemble average.
%
%   corr_gain = calc_correlation_gain(signal, heasig, references, limits, bRobust)
%
% Arguments:
%
%      + signal: signal to use
%
%      + heasig: header information about the signal.
%
%      + references: temporal samples where the synch occurs (QRS complex locations)
%
%      + limits: the amount of samples w.r.t. references(i) to consider a a
%                time window around references(i), from references(i) -
%                limits(1) to references(i) + limits(2)
%
%      + bRobust: Calculate the esemble average using mean or median.
%                 Default: use median.
%
% Output:
%
%      + corr_gain: Estimate of the SNR considering signal the ensemble
%                   average, and noise each realization minus the ensemble
%                   average.
%
% Example:
%
% See also CalcRRserieQuality
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function corr_gain = calc_correlation_gain(signal, heasig, references, limits, bRobust)

corr_gain = nan;
[nsamp, nsig] = size(signal);
pattern_size = sum(limits);
lreferences = length(references);

sig_pack = pack_signal(signal, references(references > limits(1) & references < (nsamp-limits(2) )), limits, true);

if( isempty(sig_pack) )
    noise_power = nan;
    return
end

if( nargin < 5 || isempty(bRobust) )
    pattern_avg = squeeze(median(sig_pack,3));
else
    pattern_avg = squeeze(mean(sig_pack,3));
end

sig_pack_sub_avg = bsxfun(@minus, sig_pack, pattern_avg);

corr_gain = 10*log10( mean(pattern_avg.^2) ./ 10.^squeeze(mean(log10(mean(sig_pack_sub_avg.^2,1)),3)) );

% % Obsolete: estimate noise from a random sampling simulating guessing the QRS locations.
%
% if( nargin < 4 || isempty(noise_power) )
%     noise_power = nan(30,nsig);
%     noise_disp = nan(30,nsig);
%     for ii = 1:30
%         noise_pack = pack_signal(signal, sort(randsample(limits(1):(nsamp-limits(2)), lreferences)), limits);
% %         noise_avg = flipud(reshape(squeeze(median(noise_pack,3)), pattern_size, nsig));
%         noise_limits = prctile(noise_pack, [5 50 95], 3);
%
%         noise_down = squeeze(noise_limits(:,:,1));
%         noise_avg = squeeze(noise_limits(:,:,2));
%         noise_up = squeeze(noise_limits(:,:,3));
%
%         noise_power(ii,:) = mean(noise_avg.^2);
%         noise_disp(ii,:) = mean(noise_up - noise_down);
%     end
%     noise_power = nanmedian(noise_power);
%     noise_disp = nanmedian(noise_disp);
% end
%
% corr_gain = 10*log10(mean(pattern_avg.^2) ./ noise_power );

% pattern_limits = prctile(sig_pack, [5 50 95], 3);
%
% pattern_down = squeeze(pattern_limits(:,:,1));
% pattern_up = squeeze(pattern_limits(:,:,3));
% mean_pattern_disp = mean(pattern_up - pattern_down);

% corr_gain_2 = 20*log10( noise_disp ./ mean_pattern_disp );

% debug
% aux_val = [];
% for ii  = 1:100
%     noise_pack = pack_signal(signal, sort(randsample(limits(1):(nsamp-limits(2)), length(references))), limits);
%     noise_avg = flipud(reshape(squeeze(median(noise_pack)), pattern_size, nsig));
%     aux_val = [aux_val; mean(noise_avg.^2)];
% end
%
% figure
% hist( aux_val, 100)
% legend()


text:
%
% Usage:  [Z [A,B]]=mexArchetypalAnalysis(X,param);
% Usage:  [Z [A,B]]=mexArchetypalAnalysis(X,Z,param);
%
% Name: mexArchetypalAnalysis
%
% Description: documentation to appear soon
%
% Inputs: X:  double m x n matrix   (input signals)
%               m is the signal size
%               n is the number of signals to decompose
% Output: Z: double %
%
% Author: Yuansi Chen and Julien Mairal, 2014

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function out = numeric(fa)
% Convert to numeric form
% FORMAT numeric(fa)
% fa - a file_array
%__________________________________________________________________________
% Copyright (C) 2005-2018 Wellcome Trust Centre for Neuroimaging

%
% $Id: numeric.m 7501 2018-11-30 12:16:58Z guillaume $

out = full(fa);

author: spm


text:
function [u,v]=pow2cep(m,c,mode)
%CEP2POW convert cepstral means and variances to the power domain
% Inputs:
%    m: vector giving means in the power domain
%    c: covariance matrix in the power domain
% mode: 'c'  pow=exp(irdct(cep))   [default]
%       'f'  pow=exp(rsfft(cep)/n)  [fft length even]
%       'fo' pow=exp(rsfft(cep)/n)  [fft length odd]
%       'i'  pow=exp(cep)           [ no transformation ]
%
% Outputs:
%    u: row vector giving the cepstral means with u(1) the 0'th cepstral coefficient
%    v: cepstral covariance matrix

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: pow2cep.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3 mode='c'; end
if min(size(c))==1
   v=diag(c);
end
m=m(:)';        % force to be a row vector
q=log(1+c./(m'*m));
p=log(m)-0.5*diag(q)';
if any(mode=='f')
   n=2*length(m)-2;
   if any(mode=='o')
      n=n+1;
   end
   u=rsfft(p,n);
   v=rsfft(rsfft(q,n)',n);
elseif any(mode=='i')
    u=p;
    v=q;
else
   u=rdct(p);
   v=rdct(rdct(q)');
end

author: decouples


text:
function [ heat_long, heat_lat, mse_doa ] = create_heatmap( doa_meters12, doa_meters13, doa_meters23, rx1_lat, rx1_long, rx2_lat, rx2_long, rx3_lat, rx3_long, resolution, geo_ref_lat, geo_ref_long )
%create_heatmap_kl Creates a heatmap for based on mean squared error

%    returns:
%    heat_long: longitudes of heatmap points
%    heat_lat: latitudes of heatmap points
%    mse_doa: heatmap magnitudes

    disp('creating heatmap... ');

    num_points = resolution; % points in one dimension (creates squared area)

    % defines the area, where the heatmap is displayed (around the geodetic
    % reference point)
    lat_span = 0.05;
    start_lat = geo_ref_lat - lat_span;
    stop_lat  = geo_ref_lat + lat_span;

    long_span = 0.2;
    start_long = geo_ref_long - long_span;
    stop_long  = geo_ref_long + long_span;

    % create heatmap
    heat_lat  = linspace(start_lat,  stop_lat,  num_points);
    heat_long = linspace(start_long, stop_long, num_points);
    mse_doa = zeros(num_points, num_points);

    for lat_idx = 1:num_points
        for long_idx = 1:num_points
            % calculate mean squared error of current point in terms of tdoa

            % distance current point to receivers
            dist_to_rx1 = dist_latlong( heat_lat(lat_idx), heat_long(long_idx), rx1_lat, rx1_long, geo_ref_lat, geo_ref_long );
            dist_to_rx2 = dist_latlong( heat_lat(lat_idx), heat_long(long_idx), rx2_lat, rx2_long, geo_ref_lat, geo_ref_long );
            dist_to_rx3 = dist_latlong( heat_lat(lat_idx), heat_long(long_idx), rx3_lat, rx3_long, geo_ref_lat, geo_ref_long );

            % current doa in meters
            current_doa12 = dist_to_rx1 - dist_to_rx2;
            current_doa13 = dist_to_rx1 - dist_to_rx3;
            current_doa23 = dist_to_rx2 - dist_to_rx3;

            % error doa
            doa_error = (current_doa12 - doa_meters12)^2 + (current_doa13 - doa_meters13)^2 + (current_doa23 - doa_meters23)^2;
            mse_doa(long_idx, lat_idx) = doa_error;
        end
    end

    mse_doa = 1./mse_doa;
    mse_doa = mse_doa .* (1/max(max(mse_doa)));

    disp('creating heatmap done! ');
end

repo: tdoa-evaluation-rtlsdr


text:
function mat = createRotation90(varargin)
%CREATEROTATION90  Matrix of a rotation for 90 degrees multiples
%
%   MAT = createRotation90
%   Returns the 3-by-3 matrix corresponding to a rotation by 90 degrees.
%   As trigonometric functions are explicitley converted to +1 or -1, the
%   resulting matrix obtained with this function is more precise than
%   the one obtained with createRotation.
%
%   MAT = createRotation90(NUM)
%   Specifies the number of rotations to performs. NUM should be an integer
%   (possibly negative).
%
%   Example
%     poly = [10 0;20 0;10 10];
%     rot = createRotation90;
%     poly2 = transformPoint(poly, rot);
%     figure; hold on; axis equal;
%     drawPolygon(poly);
%     drawPolygon(poly2, 'm');
%     legend('original', 'rotated');
%
%     % specify number of rotations, and center
%     rot = createRotation90(2, [10 10]);
%     poly3 = transformPoint(poly, rot);
%     drawPolygon(poly3, 'g');
%
%   See also
%   transforms2d, createRotation
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-06-20,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

% default values
num = 1;
center = [0 0];

% process input arguments
while ~isempty(varargin)
    var = varargin{1};
    if isnumeric(var) && isscalar(var)
        % extract number of rotations
        num = mod(mod(var, 4) + 4, 4);

    elseif isnumeric(var) && length(var) == 2
        % extract rotation center
        center = var;

    else
        % unknown argument
        error('MatGeom:createRotation90', ...
            'Unable to parse input arguments');
    end
    varargin(1) = [];
end

% determine rotation parameters
switch num
    case 0
        ct = 1;
        st = 0;
    case 1
        ct = 0;
        st = 1;
    case 2
        ct = -1;
        st = 0;
    case 3
        ct = 0;
        st = -1;
end

% compute transform matrix
mat = [ ...
    ct -st 0; ...
    st ct 0; ...
    0 0 1];

% change center if needed
if sum(center ~= [0 0]) > 0
    tra = createTranslation(center);
    mat = tra * mat / tra;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function x = r83_cr_sl ( n, a_cr, b )

%*****************************************************************************80
%
%% R83_CR_SL solves a real linear system factored by R83_CR_FA.
%
%  Discussion:
%
%    The matrix A must be tridiagonal.  R83_CR_FA is called to compute the
%    LU factors of A.  It does so using a form of cyclic reduction.  If
%    the factors computed by R83_CR_FA are passed to R83_CR_SL, then one or many
%    linear systems involving the matrix A may be solved.
%
%    Note that R83_CR_FA does not perform pivoting, and so the solution
%    produced by R83_CR_SL may be less accurate than a solution produced
%    by a standard Gauss algorithm.  However, such problems can be
%    guaranteed not to occur if the matrix A is strictly diagonally
%    dominant, that is, if the absolute value of the diagonal coefficient
%    is greater than the sum of the absolute values of the two off diagonal
%    coefficients, for each row of the matrix.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Hockney,
%    A fast direct solution of Poisson's equation using Fourier Analysis,
%    Journal of the ACM,
%    Volume 12, Number 1, pages 95-113, January 1965.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A_CR(3,2*N+1), factorization information computed by R83_CR_FA.
%
%    Input, real B(N), the right hand side vector.
%
%    Output, real X(N), the solution of the linear system.
%
  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R83_CR_SL - Fatal error!\n' );
    fprintf ( 1, '  Nonpositive N = %d\n', n );
    return
  end

  x = zeros ( n, 1 );

  if ( n == 1 )
    x(1) = a_cr(2,2) * b(1);
    return
  end
%
%  Set up RHS.
%
  rhs = zeros ( 2 * n + 1, 1 );

  rhs(1) = 0.0;
  rhs(2:n+1) = b(1:n);
  rhs(n+2:2*n+1) = 0.0;

  il = n;
  ndiv = 1;
  ipntp = 0;

  while ( 1 < il )

    ipnt = ipntp;
    ipntp = ipntp + il;
    il = floor ( il / 2 );
    ndiv = ndiv * 2;
    ihaf = ipntp;

    for iful = ipnt + 2 : 2 : ipntp
      ihaf = ihaf + 1;
      rhs(ihaf+1) = rhs(iful+1) - a_cr(3,iful) * rhs(iful) ...
        - a_cr(1,iful+1) * rhs(iful+2);
    end

  end

  rhs(ihaf+1) = rhs(ihaf+1) * a_cr(2,ihaf+1);
  ipnt = ipntp;

  while ( 0 < ipnt )

    ipntp = ipnt;
    ndiv = floor ( ndiv / 2 );
    il = floor ( n / ndiv );
    ipnt = ipnt - il;
    ihaf = ipntp;

    for ifulm = ipnt + 1 : 2 : ipntp
      iful = ifulm + 1;
      ihaf = ihaf + 1;
      rhs(iful+1) = rhs(ihaf+1);
      rhs(ifulm+1) = a_cr(2,ifulm+1) * ( rhs(ifulm+1) ...
        - a_cr(3,ifulm) * rhs(ifulm) ...
        - a_cr(1,ifulm+1) * rhs(iful+1) );
    end

  end

  x(1:n) = rhs(2:n+1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function frobenius_number_order2_values_test ( )

%*****************************************************************************80
%
%% FROBENIUS_NUMBER_ORDER2_VALUES_TEST tests FROBENIUS_NUMBER_ORDER2_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FROBENIUS_NUMBER_ORDER2_VALUES_TEST:\n' );
  fprintf ( 1, '  FROBENIUS_NUMBER_ORDER2_VALUES returns values of\n' );
  fprintf ( 1, '  the Frobenius number of order 2.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         C1        C2          F(C1,C2)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, c1, c2, f ] = frobenius_number_order2_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %8d  %8d  %8d\n', c1, c2, f );

  end

  return
end


text:
function fg = dtiClipFiberGroup(fg, rlClip, apClip, siClip,saveFlag)
%
% [fg fgNot] = dtiClipFiberGroup(fg, rlClip, apClip, siClip,[saveFlag])
%
% Clips a MrD fiber group. Adapted from dtiFiberUI (RFD)
% If saveFlag == 1 then the fiber group will be written to the working
% directory with the fg.name'_clip'.
%
% HISTORY:
% 03.26.2009 - LMP wrote the thing.
% 08.2009 ER added proper treatment of "subgroup" field.

if notDefined('saveFlag'), saveFlag = 0; end
if(nargin<4) siClip = []; end
if(nargin<3) apClip = []; end
newName = [fg.name '_clip'];
if(nargin==1)
    prompt = {'Left (-80) Right (+80) clip (blank for none):',...
          'Posterior (-120) Anterior (+80) clip (blank for none):',...
          'Inferior (-50) Superior (+90) clip (blank for none):',...
          'New Fiber Group name:'};
    defAns = {'','','',newName};
    resp = inputdlg(prompt,'Clip Current FG',1,defAns);
    if(isempty(resp))
        disp('User cancelled clip.');
        return;
    end
    rlClip = str2num(resp{1});
    apClip = str2num(resp{2});
    siClip = str2num(resp{3});
    fg.name = resp{4};
end
empty = zeros(size(fg.fibers));
for ii=1:length(fg.fibers)
    keep = ones(size(fg.fibers{ii}(1,:)));
    if(~isempty(rlClip))
        keep = keep & (fg.fibers{ii}(1,:)<rlClip(1) | fg.fibers{ii}(1,:)>rlClip(2));
    end
    if(~isempty(apClip))
        keep = keep & (fg.fibers{ii}(2,:)<apClip(1) | fg.fibers{ii}(2,:)>apClip(2));
    end
    if(~isempty(siClip))
        keep = keep & (fg.fibers{ii}(3,:)<siClip(1) | fg.fibers{ii}(3,:)>siClip(2));
    end
    fg.fibers{ii} = fg.fibers{ii}(:,keep);
    empty(ii) = isempty(fg.fibers{ii})|size(fg.fibers{ii},2)==1;
end
fg.fibers = fg.fibers(~empty);
if isfield(fg, 'subgroup')
fg.subgroup = fg.subgroup(~empty);
end
fg.name = newName;

if saveFlag == 1
   dtiWriteFiberGroup(fg,fg.name);
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [GWpos,GWneg] = gateway_coef_sign(W,Ci,centtype)

%   Gateway coefficient
%
%   [Gpos,Gneg] = gateway_coef_sign(W,Ci,centtype);
%
%   Gateway coefficient is a variant of participation coefficient. Similar
%   to participation coefficient, gateway coefficient measures the
%   diversity of intermodular connections of individual nodes, but this is
%   weighted by how critical these connections are to intermodular
%   connectivity (e.g., if a node is the only connection between it's
%   module and another module, it will have a higher gateway coefficient).
%
%   Inputs:     W,        undirected connection matrix with positive and
%                         negative weights
%
%               Ci,       community affiliation vector
%
%               centtype, centrality measure to use
%                         1 = Node Strength
%                         2 = Betweenness Centrality
%
%   Output:     Gpos,     gateway coefficient for positive weights
%               Gneg,     gateway coefficient for negative weights
%
%   Reference: Vargas ER, Wahl LM. Eur Phys J B (2014) 87:1-10.
%
%   Jeff Spielberg, University of Delaware

%   Modification History:
%   May 2015:  Original (originally adapted from participation_coef_sign.m)
%   July 2018: Bugfix, change in how weighted matrices are handled,
%              improvements for efficiency, additional line documentation

[~,~,Ci]       = unique(Ci);                                        % Remap module indices to consecutive numbers
n              = length(W);                                         % Number of nodes
W(1:(n+1):end) = 0;                                                 % Ensure diagonal is zero
GWpos          = gcoef(W.*(W>0));                                   % Compute gateway coefficient for positive weights
GWneg          = gcoef(-W.*(W<0));                                  % Compute gateway coefficient for negative weights

    function GW = gcoef(W_)
        k    = sum(W_,2);                                           % Compute node strength
        Gc   = (W_~=0)*diag(Ci);                                    % Create neighbor community affiliation matrix
        nmod = max(Ci);                                             % Find # of modules
        ks   = zeros(n,nmod);                                       % Preallocate space
        kjs  = zeros(n,nmod);                                       % Preallocate space
        cs   = zeros(n,nmod);                                       % Preallocate space
        switch centtype                                             % Which centrality measure to use?
            case 1                                                  % Node Strength
                cent = sum(W_,2);
            case 2                                                  % Betweenness Centrality
                L    = weight_conversion(W_,'lengths');
                cent = betweenness_wei(L);
        end
        mcn = 0;                                                    % Set max summed centrality per module to 0
        for i = 1:nmod                                              % For each module
            if sum(cent(Ci==i))>mcn                                 % If current module has a higher sum
                mcn = sum(cent(Ci==i));                             % Reassign value
            end
            ks(:,i) = sum(W_.*(Gc==i),2);                           % Compute the total weight of the connections per node to each module
        end
        for i = 1:nmod                                              % For each module
            if sum(Ci==i)>1                                         % If there is more than 1 node in a module
                kjs(Ci==i,:) = ones(sum(Ci==i),1)*sum(ks(Ci==i,:)); % Compute total module-module connections
                kjs(Ci==i,i) = kjs(Ci==i,i)/2;                      % Account for redundancy due to double counting within-network work weights
            end
        end
        for i = 1:n                                                 % For each node
            if k(i)>0                                               % If node is connected
                for ii = 1:nmod                                     % For each module
                    cs(i,ii) = sum(cent((Ci.*(W_(:,i)>0))==ii));    % Sum of centralities of neighbors of a node within a module
                end
            end
        end
        ksm           = ks./kjs;                                    % Normalize by total connections
        ksm(kjs==0)   = 0;                                          % Account for division by 0
        csm           = cs./mcn;                                    % Normalize by max summed centrality
        gs            = (1-(ksm.*csm)).^2;                          % Calculate total weighting
        GW            = 1-sum((ks.^2)./(k.^2).*gs,2);               % Compute gateway coefficient
        GW(isnan(GW)) = 0;                                          % Account for division by 0
        GW(~GW) = 0;                                                % Set to 0 if no neighbors
    end
end


text:
function [ varargout ] = TrimmedICP(Md, MovData, RefData, Tf0, MaxIter,TrMin,TrMax,lamda)
if nargin == 0
    clc; close all;
    %load ../../StanfordData/bunny;
    RefData = rand(3, 100); %bunny{1}';
    MovData = rand(3, 100);% bunny{2}';
    Md = createns(RefData');
    Tf0 = [eye(3) zeros(3, 1)];
    MaxIter = 50;
    TrMin = 0.35;
    TrMax = 1.0;
    lamda = 2.0;
end
Dim = size(MovData, 1);
PreMSE= 10^5;   CurMSE= 10^6;  Iter = 1;
RelErr = 1.0;
R0 = Tf0(1:Dim, 1:Dim);
T0 = Tf0(1:Dim, end);
Tf = [R0 T0; zeros(1, Dim) 1 ];
for Iter = 1 : 1 : MaxIter       % (abs(CurMSE-PreMSE)>10^(-9))   % this threshold is sensitive to resolutions....for bunny, this value should be 1e-12.
    TData = Loc2Glo(MovData, Tf(1:Dim, 1:Dim)', Tf(1:Dim, end) );
    [corr,TD] = knnsearch( Md,TData(1:Dim, :)');
    SortTD2 = sortrows(TD.^2); % Sort the correspongding points
    minTDIndex = floor(TrMin*length(TD)); % Get minimum index of TD
    maxTDIndex = ceil(TrMax*length(TD)); % Get maxmum index of TD
    TDIndex = [minTDIndex : maxTDIndex]';
    mTr = TDIndex./length(TD);
    mCumTD2 = cumsum(SortTD2);
    mMSE = mCumTD2(minTDIndex : maxTDIndex)./TDIndex;
    mPhi = ObjectiveFunction(mMSE, mTr);
    PreMSE=CurMSE;
    [CurMSE, nIndex] = min(mPhi);
    Trim = mTr(nIndex); % Update Tr for next step
    corr(:,2) = [1 : length(corr)]';
    % Sort the corresponding points
    corrTD = [corr, TD];
    SortCorrTD = sortrows(corrTD, 3);

    TrLength = floor(Trim*size(SortCorrTD,1)); % The number of corresponding points after trimming
    TCorr = SortCorrTD(1:TrLength, 1:2);     % Trim the corresponding points according to overlap parameter Tr
    % Register MData with TData
    dM = RegFun(RefData(:, TCorr(:, 1)), TData(:, TCorr(:, 2)) );
    % dM = reg(RefData, TData, TCorr);
    Tf = dM * Tf;
    % [M, TCorr, scan] = CalRtPhi(model, data, SortCorrTD, Trim);
    Err = [norm(dM(1:Dim, 1:Dim)-eye(Dim)) norm(dM(1:Dim, end))];
    if max(Err) <= 1e-5
        break;
    end
    RelErr = (PreMSE - CurMSE) / abs(PreMSE);
    if abs(RelErr) <= 1e-3
        break;
    end
end
IS_SHOW = 0;
if IS_SHOW
    Res = CalRes(MovData);
    h = figure;
    ICP_PlotFun(MovData, RefData, Tf, Res, h);
    title('Trimmed ICP');
end;
if nargout == 1
    varargout{1} = Tf;
end
if nargout == 2
    varargout{1} = Tf(1:Dim, 1:Dim);
    varargout{2} = Tf(1:Dim, end);
end
bTest = 1;

end

%%%%%%%%%%%%%%%%%%%%Integrated Function%%%%%%%%%%%%%%%%%%%%
% %% Calculate R,t,Phi based on current overlap parameter
% function [M,TCorr,TData] = CalRtPhi(Model, scan, SortCorrTD,Tr)
%
% TrLength = floor(Tr*size(SortCorrTD,1)); % The number of corresponding points after trimming
% TCorr = SortCorrTD(1:TrLength, 1:2);     % Trim the corresponding points according to overlap parameter Tr
% % Register MData with TData
% [M] = reg(Model(1:3,:), scan(1:3,:), TCorr);
% % To obtain the transformation data
% TData = M*scan;
%
% end
%
% function [M,TCorr,TData] = CalRtPhi(Model, scan, SortCorrTD,Tr)
%
% TrLength = floor(Tr*size(SortCorrTD,1)); % The number of corresponding points after trimming
% TCorr = SortCorrTD(1:TrLength, 1:2);     % Trim the corresponding points according to overlap parameter Tr
% % Register MData with TData
% [M] = reg(Model(1:3,:), scan(1:3,:), TCorr);
% % To obtain the transformation data
% TData = M*scan;
%
% end
%%%%%%%%%%%%%%% Calculate the registration matrix %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% T(TData)->MData %%%%%%%%%%%%%%%%%%%%%%%%%
% SVD solution
function [M] = reg(Model, Data, corr)

n = length(corr);
M = Model(:,corr(:,1));
mm = mean(M,2);
S = Data(:,corr(:,2));
ms = mean(S,2);
Sshifted = [S(1,:)-ms(1); S(2,:)-ms(2); S(3,:)-ms(3)];
Mshifted = [M(1,:)-mm(1); M(2,:)-mm(2); M(3,:)-mm(3)];
K = Sshifted*Mshifted';
K = K/n;
[U A V] = svd(K);
R1 = V*U';
if det(R1)<0
    B = eye(3);
    B(3,3) = det(V*U');
    R1 = V*B*U';
end
t1 = mm - R1*ms;
M=[];
M(1:3,1:3)=R1;
M(1:3,4)=t1;
M(4,:)=[0,0,0,1];
end

function [Phi] = ObjectiveFunction(MSE, TrB)
lamga= 2;
Phi = MSE./((TrB).^((1+lamga)));
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
% Approximates X = INV(A)*Y by a truncated series.
function x = linsolve_neumann(A,y,D,n)

if nargin < 3 && isnumeric(A)
  tol = 1e-6;
  D = (1+tol)*normest(A,tol); % this is probably veeeery bad...
end
if nargin < 4
  n = 100;
end

if isnumeric(D)
  if numel(D) == 1
    D = D*speye(length(y));
  end
  D = @(x) D\x;
end

if isnumeric(A)
  A = @(x) A*x;
end

z = D(y);
x = z;
for i=1:n
  z = z - D(A(z));
  x = x + z;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function gmt_map(sOutput, sPrefix, hAxes)
% function gmt_map(sOutput sPrefix, hAxes)
% ----------------------------------------
% Wrapper function for conversion of Matlab map-figures to GMT
%
% Input parameters:
%   sOutput         Filename without path information of the GMT script
%   sPrefix         Prefix of names of files generated by gmt_map (can contain path information)
%                   Default: empty string
%   hAxes           Handle of the axes for conversion (default: gca)
%
% Danijel Schorlemmer
% March 11, 2004

if ~exist('sPrefix', 'var')
  sPrefix = '';
end

if ~exist('hAxes', 'var')
  hAxes = gca;
end

gmt_init(hAxes, sOutput);

% Iterate thru the layers of the plot
vKids = get(hAxes, 'Children');
for nCnt = length(vKids):-1:1
  sType = lower(get(vKids(nCnt), 'Type'));
  switch sType
  case 'surface'
    [fXDiff, fYDiff] = gmt_pcolor(vKids(nCnt), sOutput, num2str(nCnt));
  case 'line'
    gmt_xy(vKids(nCnt), sOutput, num2str(nCnt));
  case 'text'
    gmt_text(vKids(nCnt), sOutput, num2str(nCnt));
  end
end

% Search for a colorbar and plot it
gmt_scale(hAxes, sOutput, '')

gmt_done(sOutput);

unix(['chmod u+x ' sOutput]);

author: CelsoReyes


text:
function [metabolite_structure,errorFlag] = sanityCheckMetIds(metabolite_structure,removeErrors)
% perform basic sanity checks

if ~exist('removeErrors', 'var')
    removeErrors = 0;
end

errorFlag = '';
Mets = fieldnames(metabolite_structure);
cnt = 1;
for i = 1 : length(Mets)
    % charge must be numeric
    if ~isnumeric(metabolite_structure.(Mets{i}).charge)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).charge;
        errorFlag{cnt,3} = 'charge is not numeric';
        cnt = 1 + cnt;
    end
    % avgmolweight must be numeric
    if  isfield(metabolite_structure.(Mets{i}),'avgmolweight')&& ~isnumeric(metabolite_structure.(Mets{i}).avgmolweight)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).avgmolweight;
        errorFlag{cnt,3} = 'avgmolweight is not numeric';
        cnt = 1 + cnt;

        if removeErrors
            % remove error
            metabolite_structure.(Mets{i}).avgmolweight = NaN;
            metabolite_structure.(Mets{i}).avgmolweight_source = [];
        end
    end
    % monoisotopicweight must be numeric
    if isfield(metabolite_structure.(Mets{i}),'monoisotopicweight')&& ~isnumeric(metabolite_structure.(Mets{i}).monoisotopicweight)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).monoisotopicweight;
        errorFlag{cnt,3} = 'monoisotopicweight is not numeric';
        cnt = 1 + cnt;
        if removeErrors
            % remove error
            metabolite_structure.(Mets{i}).monoisotopicweight = NaN;
            metabolite_structure.(Mets{i}).monoisotopicweight_source = [];
        end
    end
    % pubChemId must be numeric
    if isfield(metabolite_structure.(Mets{i}),'pubChemId')&&  ~isnumeric(metabolite_structure.(Mets{i}).pubChemId)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).pubChemId;
        errorFlag{cnt,3} = 'pubChemId is not numeric';
        cnt = 1 + cnt;

        if removeErrors
            % remove error
            metabolite_structure.(Mets{i}).pubChemId = NaN;
            metabolite_structure.(Mets{i}).pubChemId_source = [];
        end
    end
    % cheBIId must be numeric
    if isfield(metabolite_structure.(Mets{i}),'cheBIId')&&   ~isnumeric(metabolite_structure.(Mets{i}).cheBIId)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).cheBIId;
        errorFlag{cnt,3} = 'cheBIId is not numeric';
        cnt = 1 + cnt;

        if removeErrors
            % remove error
            metabolite_structure.(Mets{i}).cheBIId = NaN;
            metabolite_structure.(Mets{i}).cheBIId_source = [];
        end
    end
    % chemspider must be numeric
    if isfield(metabolite_structure.(Mets{i}),'chemspider')&& ~isnumeric(metabolite_structure.(Mets{i}).chemspider)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).chemspider;
        errorFlag{cnt,3} = 'chemspider is not numeric';
        cnt = 1 + cnt;

        if removeErrors
            % remove error
            metabolite_structure.(Mets{i}).chemspider = NaN;
            metabolite_structure.(Mets{i}).chemspider_source = [];
        end
    end
    % metlin must be numeric
    if isfield(metabolite_structure.(Mets{i}),'metlin')&& ~isnumeric(metabolite_structure.(Mets{i}).metlin)
        errorFlag{cnt,1} = Mets{i};
        errorFlag{cnt,2} = metabolite_structure.(Mets{i}).metlin;
        errorFlag{cnt,3} = 'metlin is not numeric';
        cnt = 1 + cnt;

        if removeErrors
            % remove error
            metabolite_structure.(Mets{i}).metlin = NaN;
            metabolite_structure.(Mets{i}).metlin_source = [];
        end
    end
    % keggId must start with C,D, or G
    if  isfield(metabolite_structure.(Mets{i}),'keggId') &&  isempty(find(isnan(metabolite_structure.(Mets{i}).keggId),1))
        if isnumeric(metabolite_structure.(Mets{i}).keggId) || isempty(regexp(metabolite_structure.(Mets{i}).keggId,'^C')) && isempty(regexp(metabolite_structure.(Mets{i}).keggId,'^D')) ...
                && isempty(regexp(metabolite_structure.(Mets{i}).keggId,'^G'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).keggId;
            errorFlag{cnt,3} = 'keggId does not start with C,D, or G';
            cnt = 1 + cnt;

            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).keggId = NaN;
                metabolite_structure.(Mets{i}).keggId_source = [];
            end
        end
    end
    % hmdb must start with HMDB
    if  isfield(metabolite_structure.(Mets{i}),'hmdb')&&  isempty(find(isnan(metabolite_structure.(Mets{i}).hmdb),1))
        if isnumeric(metabolite_structure.(Mets{i}).hmdb) || isempty(regexp(metabolite_structure.(Mets{i}).hmdb,'^HMDB'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).hmdb;
            errorFlag{cnt,3} = 'hmdb does not start with HMDB';
            cnt = 1 + cnt;

            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).hmdb = NaN;
                metabolite_structure.(Mets{i}).hmdb_source = [];
            end
        end
    end
    % reconMap3 must be equal to VMHId
    %     if isnumeric(metabolite_structure.(Mets{i}).reconMap3) || isempty(strmatch(metabolite_structure.(Mets{i}).VMHId,metabolite_structure.(Mets{i}).reconMap3,'exact')) && ...
    %             ~isempty(find(isnan(metabolite_structure.(Mets{i}).reconMap3)))
    %         errorFlag{cnt,1} = Mets{i};
    %         errorFlag{cnt,2} = metabolite_structure.(Mets{i}).reconMap3;
    %         errorFlag{cnt,3} = 'reconMap3 is not equal to VMHId';
    %         cnt = 1 + cnt;
    %     end

    % food_db must start with FDB
    if isfield(metabolite_structure.(Mets{i}),'food_db')&&   isempty(find(isnan(metabolite_structure.(Mets{i}).food_db),1))
        if isnumeric(metabolite_structure.(Mets{i}).food_db) ||isempty(regexp(metabolite_structure.(Mets{i}).food_db,'^FDB'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).food_db;
            errorFlag{cnt,3} = 'food_db does not start with FDB';
            cnt = 1 + cnt;

            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).food_db = NaN;
                metabolite_structure.(Mets{i}).food_db_source = [];
            end
        end
    end
    % metanetx must start with MNXM
    if isfield(metabolite_structure.(Mets{i}),'metanetx')&&  isempty(find(isnan(metabolite_structure.(Mets{i}).metanetx),1))
        if isnumeric(metabolite_structure.(Mets{i}).metanetx) ||isempty(regexp(metabolite_structure.(Mets{i}).metanetx,'^MNXM'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).metanetx;
            errorFlag{cnt,3} = 'metanetx does not start with MNXM';
            cnt = 1 + cnt;

            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).metanetx = NaN;
                metabolite_structure.(Mets{i}).metanetx_source = [];
            end
        end
    end
    % epa_id must start with DTX
    if  isfield(metabolite_structure.(Mets{i}),'epa_id')&&  isempty(find(isnan(metabolite_structure.(Mets{i}).epa_id),1))
        if isnumeric(metabolite_structure.(Mets{i}).epa_id) ||isempty(regexp(metabolite_structure.(Mets{i}).epa_id,'^DTX'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).epa_id;
            errorFlag{cnt,3} = 'epa_id does not start with DTX';
            cnt = 1 + cnt;

            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).epa_id = NaN;
                metabolite_structure.(Mets{i}).epa_id_source = [];
            end
        end
    end
    % seed must start with cpd
    if  isfield(metabolite_structure.(Mets{i}),'seed')&&  isempty(find(isnan(metabolite_structure.(Mets{i}).seed),1))
        if isnumeric(metabolite_structure.(Mets{i}).seed) ||isempty(regexp(metabolite_structure.(Mets{i}).seed,'^cpd'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).seed;
            errorFlag{cnt,3} = 'seed does not start with cpd';
            cnt = 1 + cnt;
            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).seed = NaN;
                metabolite_structure.(Mets{i}).seed_source = [];
            end
        end
    end
    % casRegistry must contain 2 -
    if isfield(metabolite_structure.(Mets{i}),'casRegistry')&& isempty(find(isnan(metabolite_structure.(Mets{i}).casRegistry),1))
        if isnumeric(metabolite_structure.(Mets{i}).casRegistry) ||isempty(regexp(metabolite_structure.(Mets{i}).casRegistry,'\d*-\d*-\d*'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).casRegistry;
            errorFlag{cnt,3} = 'casRegistry does not follow \d*-\d*-\d*';
            cnt = 1 + cnt;
            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).casRegistry = NaN;
                metabolite_structure.(Mets{i}).casRegistry_source = [];
            end
        end
    end
    % inchiKey must contain 2 -
    if  isfield(metabolite_structure.(Mets{i}),'inchiKey')&&  isempty(find(isnan(metabolite_structure.(Mets{i}).inchiKey),1))
        if isnumeric(metabolite_structure.(Mets{i}).inchiKey) ||isempty(regexp(metabolite_structure.(Mets{i}).inchiKey,'[A-Z]+-[A-Z]+-[A-Z]'))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).inchiKey;
            errorFlag{cnt,3} = 'inchiKey does not follow [A-Z]+-[A-Z]+-[A-Z]+';
            cnt = 1 + cnt;
            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).inchiKey = NaN;
                metabolite_structure.(Mets{i}).inchiKey_source = [];
            end
        end
    end
    % inchiString must start with InChI=
    if isfield(metabolite_structure.(Mets{i}),'inchiString')&& isempty(find(isnan(metabolite_structure.(Mets{i}).inchiString),1))
        if isnumeric(metabolite_structure.(Mets{i}).inchiString) ||isempty(regexp(metabolite_structure.(Mets{i}).inchiString,'^InChI='))
            errorFlag{cnt,1} = Mets{i};
            errorFlag{cnt,2} = metabolite_structure.(Mets{i}).inchiString;
            errorFlag{cnt,3} = 'inchiString does not follow ^InChI=';
            cnt = 1 + cnt;
            if removeErrors
                % remove error
                metabolite_structure.(Mets{i}).inchiString = NaN;
                metabolite_structure.(Mets{i}).inchiString_source = NaN;
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function[nui,om2i]=morlfreq(omi)
%MORLFREQ  Compute Morlet wavelet carrier frequency given peak frequency.
%
%   FNU=MORLFREQ(FMAX), where FMAX is the desired peak frequency of a
%   Morlet wavelet, returns the appropriate carrier wave frequency FNU.
%
%   The Morlet wavelet carrier wave and peak frequency differ on account
%   of the correction term which ensures zero mean.
%
%   Note that FMAX and FNU are both *radian* as in cos(omega t) and not
%   cyclic as in in cos(2 pi f t).
%
%   [FNU,FMIN]=MORLFREQ(FMAX) also returns the (radian) frequency at which
%   the Morlet wavelet obtains a minimum, which occurs for FMIN<0.
%
%   For details see Appendix A of
%
%       Lilly and Olhede (2009).  Higher-order properties of analytic
%          wavelets.  IEEE Trans. Sig. Proc., 57 (1), 146--160.
%
%   'morlfreq --t' runs a test.
%   'morlfreq --f' generates a sample figure.
%
%   Usage: fnu=morlfreq(fmax);
%          [fnu,fmin]=morlfreq(fmax);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008--2013 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(omi, '--t')
    morlfreq_test,return
end
if strcmpi(omi, '--f')
    type makefigs_morlfreq
    makefigs_morlfreq;
    return
end

sizeom=size(omi);
omi=omi(:);

N=10000;

%This gives a tighter spacing near one than near zero
nutilde=log(log(linspace(exp(exp(1./1e9)),exp(exp(1-1./1e9)),N)));
%nutilde=linspace(1./1e9,1-1./1e9,N);

om=(sqrt(-frac((log(1-nutilde)),nutilde)));
nu=nutilde.*om;

bool=(omi>=maxmax(om));
index1=find(bool);
if ~isempty(index1)
    nui(index1)=omi(index1);
end
index2=find(~bool);
if ~isempty(index2)
    nui(index2)=interp1(om,nu,omi(index2),'linear');
end

%I chose this grid via trial and error
nutilde=1./[linspace(1./10,10,N) logspace(log10(10+0.0000001),log10(1000),N)];
om2=-(sqrt(frac((log(1+nutilde)),nutilde)));
nu2=(nutilde).*abs(om2);

om2i=zeros(size(nui));%figure,plot(nui),maxmax(nu2)
bool=(nui>=maxmax(nu2));
index1=find(bool);
if ~isempty(index1)
    om2i(index1)=nan;
end
index2=find(~bool);%length(index2)
if ~isempty(index2)
    om2i(index2)=interp1(nu2,om2,nui(index2),'linear');
end
%x=om2i-nui-om2i.*exp(-om2i.*nui);
%figure,plot(x)

nui=reshape(nui,sizeom);
om2i=reshape(om2i,sizeom);

function[]=morlfreq_test
om=linspace(0,4,1000);
[nu,fneg]=morlfreq(om);

psi1a=exp(-frac(1,2).*(nu-om).^2).*(-(om-nu))+exp(-frac(1,2).*(nu.^2+om.^2)).*om;
reporttest('MORLFREQ first derivative vanishes at peak',aresame(psi1a,0*psi1a,1e-4))

omneg=fneg*2*pi;
psi1b=exp(-frac(1,2).*(nu-omneg).^2).*(-(omneg-nu))+exp(-frac(1,2).*(nu.^2+omneg.^2)).*omneg;

reporttest('MORLFREQ first derivative vanishes at negative peak',aresame(psi1b,0*psi1b,1e-4))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab