text:
grayFileName='test.gray';

[gNodes, gEdges, gvSize] = readGrayGraph_progress(grayFileName,0);

% Do all the gray node stages in a loop so that we can have arbitrary
% numbers of layers.
numGrayLayers=max(gNodes(6,:));

% Get the indices for all the gnodes (all layers)
for t=1:numGrayLayers
    grayNodeIndices{t}=find(gNodes(6,:)==t);
end

% Extract the layer 1 nodes. These are special because they are used to map
% down to the unfolded boundary mesh.

l1gNodes=gNodes(:,grayNodeIndices{1});
l1mesh.vertices=l1gNodes(1:3,:);
l1mesh.indices=grayNodeIndices{1};

% How many gNodes are there?
nGnodes=length(gNodes);
% How many gEdges are there?
nGedges=length(gEdges);

% We want to make a grey node connection matrix - which grey nodes are connected to which other gnodes?
%statusStringAdd(statusHandle,'Finding grey connection matrix (slow)');

grayConMat=makeGrayConMat(gNodes,gEdges,0);

% We can assign layer 1 grey nodes to the white matter mesh using assignToNearest.dll (see assignToNearest.c)
% Can't do this for higher levels of grey matter 'cos they might get mis-assigned. (Also, potential problem near
% very crinkly edges. - Could we accidentally assign a l1 grey matter node to the wrong WM point?)
% (I think the answer is 'yes, rarely'. If a single layer or gray is sandwiched between two sides of a sulcus (say) : It's grown from
% one side but mrFlatMesh has no way of telling which one.
% Going deeper into mrGray's source code to determine the parentage of l1 gray nodes might be possible...)

% So for higher grey matter points, we have to restrict the possible sub node search space by assigning them >only< to
% points they are connected to. Note that a single layer2 grey node may be connected to several l1 nodes

% The gray may be defined over the entire mesh but we only want to deal with gray points over the
% unfolded part. The strategy should be....
% 1) do assignToNearest for each mesh point to find the nearest connected l1 node
% 2) Use the set of l1 nodes found in 1) to build up a list of other connected gray nodes
% 3) Repeat stage 2 for l3,l4

%statusStringAdd(statusHandle,'Mapping L1 to mesh.');

% Find 3D coords of all the l1 gnodes
l1GNodeCoords=l1gNodes(1:3,:)';

% Find conmat of l1 gNodes
l1conMat=grayConMat(grayNodeIndices{1},grayNodeIndices{1});

% Find the 3d neighbour distances
mesh.connectionMatrix=l1conMat;
mesh.uniqueVertices=l1GNodeCoords;
scaleFactor=[1 1 1];

D=find3DNeighbourDists(mesh,scaleFactor); % We're now in a voxel framework: everything has been scaled in the mrReadMrM function.
perimDist=30;
startNode=1000; % I don't know where this is
mesh.dist=dijkstra(D,startNode);

mesh.perimDist=perimDist;

[perimeterEdges,eulerCondition]=findLegalPerimeters(mesh,perimDist);
gplot3(mesh.connectionMatrix,mesh.uniqueVertices) ;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [sol, val] = gaDemo1Eval(sol,options)
% Demonstration evaluation function used in gademo1.
% f(x)=x+10sin(5x)+7cos(4x)
%
% function [val,sol] = gaDemo1Eval(sol,options)
%
% val - the fittness of this individual
% sol - the individual, returned to allow for Lamarckian evolution
% options - [current_generation]

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

x=sol(1);
val = x + 10*sin(5*x)+7*cos(4*x);

author: HuangCongQing


text:
function grains = minus(grains,xy)
% shift grains in x/y direction
%
% Syntax
%
%   % shift in x direction
%   grains = grains - [100,0]
%
% Input
%  grains- @grain2d
%  xy - x and y coordinates of the shift
%
% Output
%  grains - @grain2d

grains = grains + (-xy);

author: mtex-toolbox

repo: mtex

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


text:
%performs knn classification
%>
%> @param TestFeatureVector: features for test observation (length iNumFeatures)
%> @param TrainFeatureMatrix: features for all train observations (dimension iNumFeatures x iNumObservations)
%> @param TrainClassIndices: audio signal (length iNumObservations)
%> @param k: number of points taken into account (default = 3)
%>
%> @retval class index of the resulting class
% ======================================================================
function [class] = ToolSimpleKnn(TestFeatureVector, TrainFeatureMatrix, TrainClassIndices, k)

    % set order to 3 if not set
    if (nargin < 4)
        k = 3;
    end

    % compute distances to all training observations
    d = computeEucDist_I(TestFeatureVector, TrainFeatureMatrix);

    % sort the distances to find closest
    [dummy,idx] = sort(d);

    % pick the majority of the k closest training observations
    % note that for multi-class problems and even k, this needs to be
    % refined
    class = mode(TrainClassIndices(idx(1:k)));
end

function d = computeEucDist_I(A, B)
    d = sqrt(sum(A.^2, 2)*ones(1,size(B,1)) - ...
        2*A*B' + ...
        ones(size(A,1),1)*sum(B.^2, 2)');
end

author: alexanderlerch

repo: ACA-Code


text:
%PLOT_ELLIPSE Draw an ellipse or ellipsoid
%
% plot_ellipse(E, OPTIONS) draws an ellipse or ellipsoid defined by X'EX =
% 0 on the current plot, centred at the origin.  E (2x2) for an ellipse and
% E (2x3) for an ellipsoid.
%
% plot_ellipse(E, C, OPTIONS) as above but centred at C=[X,Y].  If
% C=[X,Y,Z] the ellipse is parallel to the XY plane but at height Z.
%
% H = plot_ellipse(...) as above but return graphic handle.
%
% Options::
% 'confidence',C   confidence interval, range 0 to 1
% 'alter',H        alter existing ellipses with handle H
% 'npoints',N      use N points to define the ellipse (default 40)
% 'edgecolor'      color of the ellipse boundary edge, MATLAB color spec
% 'fillcolor'      the color of the ellipses's interior, MATLAB color spec
% 'alpha'          transparency of the fillcolored ellipse: 0=transparent, 1=solid
% 'shadow'         show shadows on the 3 walls of the plot box
%
% - For an unfilled ellipse:
%   - any standard MATLAB LineStyle such as 'r' or 'b---'.
%   - any MATLAB LineProperty options can be given such as 'LineWidth', 2.
% - For a filled ellipse any MATLAB PatchProperty options can be given.
%
% Example::
%
%          H = plot_ellipse(diag([1 2]), [3 4]', 'r'); % draw red ellipse
%          plot_ellipse(diag([1 2]), [5 6]', 'alter', H); % move the ellipse
%          plot_ellipse(diag([1 2]), [5 6]', 'alter', H, 'LineColor', 'k'); % change color
%
%          plot_ellipse(COVAR, 'confidence', 0.95); % draw 95% confidence ellipse
%
% Notes::
% - The 'alter' option can be used to create a smooth animation.
% - If E (2x2) draw an ellipse, else if E (3x3) draw an ellipsoid.
% - The ellipse is added to the current plot irrespective of hold status.
% - Shadow option only valid for ellipsoids.
% - If a confidence interval is given then E is interpretted as a covariance
%   matrix and the ellipse size is computed using an inverse chi-squared function.
%   This requires CHI2INV in the Statistics and Machine Learning Toolbox or
%   CHI2INV_RTB from the Robotics Toolbox for MATLAB.
%
% See also PLOT_ELLIPSE_INV, PLOT_CIRCLE, PLOT_BOX, PLOT_POLY, CH2INV.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function handles = plot_ellipse(E, varargin)

    assert(size(E,1) == size(E,2), 'ellipse is defined by a square matrix');
    assert( size(E,1) == 2 || size(E,1) == 3, 'can only plot ellipsoid for 2 or 3 dimenions');

    opt.fillcolor = 'none';
    opt.alpha = 1;
    opt.edgecolor = 'k';
    opt.alter = [];
    opt.npoints = 40;
    opt.shadow = false;
    opt.confidence = [];

    [opt,arglist,ls] = tb_optparse(opt, varargin);

    % process some arguments

    if ~isempty(ls)
        opt.edgecolor = ls{1};
    end

    % process the probability
    if isempty(opt.confidence)
        s = 1;
    else
        if exist('chi2inv') == 2
            s = sqrt(chi2inv(opt.confidence, 2));
        elseif exist('chi2inv_rtb') == 2
            s = sqrt(chi2inv_rtb(opt.confidence, 2));
        else
            error('SMTB:missingfunc', 'Requires Stats Toolbox or RTB to be installed');
        end
    end

    if length(arglist) > 0 && isnumeric(arglist{1})
        % ellipse centre is provided
        centre = arglist{1};
        arglist = arglist(2:end);
    else
        % default to origin
        centre = zeros(1, size(E,1));
    end

    % check the ellipse to be altered
    if ~isempty(opt.alter) & ~ishandle(opt.alter)
        error('SMTB:plot_ellipse:badarg', 'argument to alter must be a valid graphic object handle');
    end

    holdon = ishold();
    hold on

    if size(E,1) == 3
        %% plot an ellipsoid

        % define mesh points on the surface of a unit sphere
        [Xs,Ys,Zs] = sphere();
        ps = [Xs(:) Ys(:) Zs(:)]';

        % warp it into the ellipsoid
        pe = sqrtm(E) * ps;

        % offset it to optional non-zero centre point
        if nargin > 1
            pe = bsxfun(@plus, centre(:), pe);
        end

        % put back to mesh format
        Xe = reshape(pe(1,:), size(Xs));
        Ye = reshape(pe(2,:), size(Ys));
        Ze = reshape(pe(3,:), size(Zs));

        if isempty(opt.alter)
              % plot it
%             Ce = ones(size(Xe));
%             Ce = cat(3, Ce*0.8, Ce*0.4, Ce*0.4);
            h = mesh(Xe, Ye, Ze, 'FaceColor', opt.fillcolor, ...
                        'FaceAlpha', opt.alpha, 'EdgeColor', opt.edgecolor, arglist{:});
        else
            % update an existing plot
            set(opt.alter, 'xdata', Xe, 'ydata', Ye, 'zdata', Ze,  ...
                        arglist{:});
        end

        % draw the shadow
        if opt.shadow
            I = ones(size(Xe));
            a = [xlim ylim zlim];
            mesh(a(1)*I, Ye, Ze, 'FaceColor', 0.7*[1 1 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
            mesh(Xe, a(3)*I, Ze, 'FaceColor', 0.7*[1 1 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
            mesh(Xe, Ye, a(5)*I, 'FaceColor', 0.7*[1 1 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
        end

    else
        %% plot an ellipse

        [V,D] = eig(E);

        % define points on a unit circle
        th = linspace(0, 2*pi, opt.npoints);
        pc = [cos(th);sin(th)];

        % warp it into the ellipse
        pe = sqrtm(E)*pc * s;

        % offset it to optional non-zero centre point
        centre = centre(:);
        if nargin > 1
            pe = bsxfun(@plus, centre(1:2), pe);
        end
        x = pe(1,:); y = pe(2,:);

%         if length(centre) > 2
%             % plot 3D data
%             z = ones(size(x))*centre(3);
%             if isempty(opt.alter)
%                 h = plot3(x', y', z', varargin{:});
%             else
%                 set(opt.alter, 'xdata', x, 'ydata', y, 'zdata', z, arglist{:});
%             end

            % plot 2D data

            if length(centre) > 2
                % plot 3D data
                z = ones(size(x))*centre(3);
            else
                z = zeros(size(x));
            end

            if strcmpi(opt.fillcolor, 'none')
                % outline only, draw a line

                if isempty(ls)
                    if ~isempty(opt.edgecolor)
                        arglist = ['Color', opt.edgecolor, arglist];
                    end
                else
                    arglist = [ls arglist];
                end

                if isempty(opt.alter)
                    h = plot3(x', y', z', arglist{:});
                else
                    set(opt.alter, 'xdata', x, 'ydata', y);
                end
            else
                % fillcolored, use a patch

                if ~isempty(opt.edgecolor)
                    arglist = ['EdgeColor', opt.edgecolor, arglist];
                end

                arglist = [ls, 'FaceAlpha', opt.alpha, arglist];

                if isempty(opt.alter)
                    h = patch(x', y', z', opt.fillcolor, arglist{:});
                else
                    set(opt.alter, 'xdata', x, 'ydata', y);
                end

            end
        end

  if ~holdon
      hold off
  end

    if nargout > 0
        handles = h;
    end
end

author: petercorke

repo: spatialmath-matlab


text:
function [ a, b ] = p04_ab ( m )

%*****************************************************************************80
%
%% P04_AB returns bounds for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real A(M,1), B(M,1), lower and upper bounds.
%
  a(1:m,1) = 0.0;
  b(1:m,1) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fm,am,iflaw]=doppler(N,Fs,f0,d,v,t0,c);
%DOPPLER Generate complex Doppler signal.
% 	 [FM,AM,IFLAW]=DOPPLER(N,FS,F0,D,V,T0,C)
%	 Returns the frequency modulation (FM), the amplitude
%	 modulation (AM) and the instantaneous frequency law (IFLAW)
%	 of the signal received by a fixed observer from a moving target
%	 emitting a pure frequency f0.
%
%	 N  : number of points.
%	 FS : sampling frequency (in Hertz).
%	 F0 : target   frequency (in Hertz).
%	 D  : distance from the line to the observer (in meters).
%	 V  : target velocity    (in m/s)
%	 T0 : time center                  (default : N/2).
%	 C  : wave velocity      (in m/s)  (default : 340).
%	 FM : Output frequency modulation.
%	 AM : Output amplitude modulation.
%	 IFLAW : Output instantaneous frequency law.
%
%	Example:
%	 N=512; [fm,am,iflaw]=doppler(N,200,65,10,50);
%	 subplot(211); plot(real(am.*fm));
%	 subplot(212); plot(iflaw);
%        [ifhat,t]=instfreq(sigmerge(am.*fm,noisecg(N),15),11:502,10);
%        hold on; plot(t,ifhat,'g'); hold off;
%
%	See also DOPNOISE

%	F. Auger, July 94, August 95 - O. Lemoine, October 95.
%	Copyright (c) 1996 by CNRS (France).
%
%  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 should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin <= 4),
 error ( 'At least 5 parameters are required' );
elseif (nargin == 5),
 t0=N/2; c=340.0;
elseif (nargin == 6),
 c=340.0;
end;

if (N <= 0),
 error ('The signal length N must be strictly positive' );
elseif (d <= 0.0),
 error ('The distance D must be positive' );
elseif (Fs < 0.0),
 error ('The sampling frequency FS must be positive' );
elseif (t0<1) | (t0>N),
 error ('T0 must be between 1 and N');
elseif (f0<0)|(f0>Fs/2),
 error ('F0 must be between 0 and FS/2');
elseif (v<0),
 error ('V must be positive');
else
 tmt0=((1:N)'-t0)/Fs;
 dist=sqrt(d^2+(v*tmt0).^2);
 fm = exp(j*2.0*pi*f0*(tmt0-dist/c));
 if (nargout>=2),
  if abs(f0)<eps,
   am=0;
  else
   am= 1.0 ./ sqrt(dist);
  end
 end;
 if (nargout==3), iflaw=(1-v^2*tmt0./dist/c)*f0/Fs; end;
end ;

author: HeLiangHIT


text:
clear

load('./results/results_clnf_wild.mat');

[clnf_error, ~,~,frontal_ids] = compute_error_small( experiment.labels,  experiment.shapes-1.0);
clnf_error_frontal = clnf_error(frontal_ids);
clnf_error_profile = clnf_error(~frontal_ids);

load('./results/results_ceclm_general.mat');

[ceclm_error,~,~,frontal_ids] = compute_error_small( experiment.labels,  experiment.shapes-1.0);
ceclm_error_frontal = ceclm_error(frontal_ids);
ceclm_error_profile = ceclm_error(~frontal_ids);
labels = experiment.labels;

load('results/CFAN_JANUS.mat');

[cfan_error,~,~,frontal_ids] = compute_error_small( labels_all,  shapes_all-1.0);
cfan_error_frontal = cfan_error(frontal_ids);
cfan_error_profile = cfan_error(~frontal_ids);

load('results/JANUS_3DDFA.mat');

[error_3ddfa,~,~,frontal_ids] = compute_error_small( labels_all,  shapes-1.0);
error_3ddfa_frontal = error_3ddfa(frontal_ids);
error_3ddfa_profile = error_3ddfa(~frontal_ids);

load('results/JANUS_pocr.mat');

[pocr_error,~,~,frontal_ids] = compute_error_small( labels,  experiments.shapes-0.5);
pocr_error_frontal = pocr_error(frontal_ids);
pocr_error_profile = pocr_error(~frontal_ids);

load('results/JANUS_chehra.mat');

[drmf_error,~,~,frontal_ids] = compute_error_small( labels,  shapes);
drmf_error_frontal = drmf_error(frontal_ids);
drmf_error_profile = drmf_error(~frontal_ids);

load('results/JANUS_sdm.mat');

[sdm_error,~,~,frontal_ids] = compute_error_small( labels,  experiments.shapes);
sdm_error_frontal = sdm_error(frontal_ids);
sdm_error_profile = sdm_error(~frontal_ids);

load('results/JANUS-CFSS.mat');
shapes = zeros(68,2, size(estimatedPose,1));

for i = 1:size(shapes, 3)
    shapes(:,1,i) = estimatedPose(i,1:68);
    shapes(:,2,i) = estimatedPose(i,69:end);
end

[cfss_error,~,~,frontal_ids] = compute_error_small( experiments.labels, shapes-1.0);
cfss_error_frontal = cfss_error(frontal_ids);
cfss_error_profile = cfss_error(~frontal_ids);

load('results/tcdcn_JANUS.mat');
shapes_c = shapes;
shapes = zeros(68,2, numel(shapes_c));

for i = 1:size(shapes, 3)
    shapes(:,1,i) = shapes_c{i}(:,1);
    shapes(:,2,i) = shapes_c{i}(:,2);
end

[tcdcn_error,~,~,frontal_ids] = compute_error_small( experiments.labels, shapes-1.0);
tcdcn_error_frontal = tcdcn_error(frontal_ids);
tcdcn_error_profile = tcdcn_error(~frontal_ids);

author: TadasBaltrusaitis

repo: OpenFace


text:
function [I_dx, I_dy] = sobel5x5(img)
%SOBEL5X5 Apply Sobel–Feldman filter over input image
% Sobel filters are discrete differentiation operator. They are used to
% compute an approximation of the gradient of an image intensity function
%
% INPUT:
%   - img(image): Given input image
%
% OUTPUT:
%   - I_dx: gradient of image along x-axis
%   - I_dy: gradient of image along y-axis

% sobel kernels
sobel_kernel1 = [1, 4, 6, 4, 1];
sobel_kernel2 = [1, 2, 0, -2, -1];
divisor = 48;

% gradient along x
I_dx = conv2(sobel_kernel2', sobel_kernel1, img, 'valid');
I_dx = uint8(I_dx / divisor + 128);

% gradient along y
I_dy = conv2(sobel_kernel1', sobel_kernel2, img, 'valid');
I_dy = uint8(I_dy / divisor + 128);

% imshowpair(I_dx, I_dy, 'montage');
end

repo: Stereo-Odometry-SOFT


text:
function y = deadZone(u, thr)

    % DEADZONE implements a dead zone for the input u. If an element of u is
    %          lower (absolute value) than a user-defined threshold,
    %          consider that element of u as zero. u can be either a vector
    %          or a scalar.
    %
    % FORMAT: y = deadZone(u, thr)
    %
    % INPUT:  - u   = [n * 1] input vector;
    %         - thr = user-defined threshold;
    %
    % OUTPUT: - y   = [n * 1] filtered output.
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2018
    %

    %% --- Initialization ---

    % initialize output
    y = u;

    for i = 1:length(u)

        if abs(u(i)) < thr

            y(i) = 0;
        end
    end
end

author: robotology

repo: whole-body-controllers


text:
function z=fun00(x)
z=0;
for i=1:length(x)
    z=z+x(i)*sin(abs(x(i)));
end

author: Sable

repo: mcbench-benchmarks


text:
  function [cost grad] = mri_r2_fit_costgrad_pd(pd, data)
%|function [cost grad] = mri_r2_fit_costgrad_pd(pd, data)
%| gradient wrt "proton density" (pd) image
%| for R2=1/T2 estimation from images

if nargin == 1 && streq(pd, 'test'), mri_r2_fit_costgrad_pd_test, return, end
if nargin < 2, ir_usage, end

yi = data.yi; % [(nd) nt] images
r2 = data.r2; % [(nd)] R2=1/T2 maps
te = data.te; % echo times

[yb yi] = mri_r2_fit_mean(pd, r2, yi, te);
cost = sum(abs(yi(:) - yb(:)).^2) / 2;

grad = mri_r2_fit_grad_pd(pd, r2, yi - yb, te);

if isfield(data, 'R')
	R = data.R;
	cost = cost + R.penal(R, pd);
	grad = grad + R.cgrad(R, pd);
end
grad = reshape(grad, size(pd));

% mri_r2_fit_mean()
function [yb yi] = mri_r2_fit_mean(pd, r2, yi, te)

nt = numel(te);
yi = reshapee(yi, [], nt); % [*nd nt]
yb = zeros(size(yi)); % [*nd nt]
for it=1:nt
	yb(:,it) = pd(:) .* exp(-te(it) * r2(:));
end

% mri_r2_fit_grad_pd()
function grad = mri_r2_fit_grad_pd(pd, r2, resid, te)

nt = numel(te);

grad = 0;
for it=1:nt
	tmp = exp(-te(it) * r2(:));
	grad = grad - tmp .* resid(:,it);
end

function mri_r2_fit_costgrad_pd_test

r2 = 0.01;
te = [0:3] * 20;
pd = 2;
yb = pd .* exp(-te * r2);
rng(0)
yi = yb + 0.1 * randn(size(yb));

%R = Reg1(true(1), 'beta', 3, 'order', 0, 'type_penal', 'mat');

data.yi = yi;
data.r2 = r2;
data.te = te;
%data.R = R;

pdlist = linspace(0, 4, 101);
cost = zeros(size(pdlist));
for ii=1:numel(pdlist)
	cost(ii) = mri_r2_fit_costgrad_pd(pdlist(ii), data);
end

p0 = 1.5;
[cost0 grad0] = mri_r2_fit_costgrad_pd(p0, data);

if im
	clf
	plot(pdlist, cost, '-', pdlist, cost0 + grad0*(pdlist-p0), '--')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function model=getMetaCycModelForOrganism(organismID,fastaFile,...
    keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond)
% getMetaCycModelForOrganism
%   Reconstructs a genome-scale metabolic model based on protein homology to the
%   MetaCyc pathway database
%
%   Input:
%   organismID          the query organism's abbreviation, which is defined
%                       by user
%   fastaFile           a FASTA file that contains the protein sequences of
%                       the organism for which to reconstruct a model
%   keepTransportRxns   include transportation reactions, which often have identical
%                       reactants and products that turn to be all-zero columns in
%                       the S matrix (opt, default false)
%   keepUnbalanced      include reactions cannot be unbalanced reactions, usually
%                       because they are polymeric reactions or because of a
%                       specific difficulty in balancing class structures
%                       (opt, default false)
%   keepUndetermined    include reactions that have substrates lack chemical
%                       structures or with non-numerical coefficients (e.g. n+1)
%                       (opt, default false)
%   minScore            minimum Bit scores of BLASTp search (opt, default 100)
%   minPositives        minimum Positives values of BLASTp search (opt, default 45 %)
%   useDiamond          use DIAMOND alignment tools to perform homology search
%                       if true, otherwise the BLASTP is used (opt, default true)
%
%   Output:
%   model               a model structure for the query organism
%
%   Usage: model=getMetaCycModelForOrganism(organismID,fastaFile,...
%    keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond)

organismID=char(organismID);
if nargin<2
    EM='No query protein fasta file is specified';
    dispEM(EM);
else
    fastaFile=char(fastaFile);
end
if nargin<3
    keepTransportRxns=false;
end
if nargin<4
    keepUnbalanced=false;
end
if nargin<5
    keepUndetermined=false;
end
if nargin<6
    minScore=100;
end
if nargin<7
    minPositives=45;
end
if nargin<8
    useDiamond=true;
end

if isempty(fastaFile)
    error('*** The query FASTA filename cannot be empty! ***');
else
    fprintf('\nChecking existence of query FASTA file... ');
    %Check if query fasta exists
    fastaFile=checkFileExistence(fastaFile,2); %Copy file to temp dir
    fprintf('done\n');
end

%First generate the full MetaCyc model
metaCycModel=getModelFromMetaCyc([],keepTransportRxns,keepUnbalanced,keepUndetermined);
fprintf('The full MetaCyc model loaded\n');

%Create the draft model from MetaCyc super model model=metaCycModel;
model.id=organismID;
model.name='Generated by homology with MetaCyc database';
model.rxns=metaCycModel.rxns;
model.rxnNames=metaCycModel.rxnNames;
model.eccodes=metaCycModel.eccodes;
model.subSystems=metaCycModel.subSystems;
model.rxnMiriams=metaCycModel.rxnMiriams;
model.rxnReferences=metaCycModel.rxnReferences;
model.lb=metaCycModel.lb;
model.ub=metaCycModel.ub;
model.rev=metaCycModel.rev;
model.c=metaCycModel.c;
model.equations=metaCycModel.equations;

%Get the root directory for RAVEN Toolbox.
ravenPath=findRAVENroot();

%Generate blast strcture by either DIAMOND or BLASTP
if useDiamond
    blastStruc=getDiamond(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa'));
else
    blastStruc=getBlast(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa'));
end

%Only look the query
blastStructure=blastStruc(2);

%Remove all blast hits that are below the cutoffs
indexes=blastStructure.bitscore>=minScore & blastStructure.ppos>=minPositives;
blastStructure.toGenes(~indexes)=[];
blastStructure.fromGenes(~indexes)=[];
blastStructure.evalue(~indexes)=[];
blastStructure.aligLen(~indexes)=[];
blastStructure.identity(~indexes)=[];
blastStructure.bitscore(~indexes)=[];
blastStructure.ppos(~indexes)=[];
fprintf('Completed searching against MetaCyc protein sequences.\n');

% Get the qualified genes of query organism from blast structure
model.genes=cell(10000,1);
model.proteins=cell(10000,1);
model.bitscore=zeros(10000,1);
model.ppos=zeros(10000,1);
num=1;

%Go through the strucutre and find out the hit with the best bit score
for i=1:numel(blastStructure.fromGenes)
    if num==1
        model.genes(num)=blastStructure.fromGenes(i);
        model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
        model.bitscore(num)=blastStructure.bitscore(i);
        model.ppos(num)=blastStructure.ppos(i);
        num=num+1;
        lastGene=blastStructure.fromGenes(i);
    else
        if ~isequal(lastGene,blastStructure.fromGenes(i))
            model.genes(num)=blastStructure.fromGenes(i);
            model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
            model.bitscore(num)=blastStructure.bitscore(i);
            model.ppos(num)=blastStructure.ppos(i);
            num=num+1;
            lastGene=blastStructure.fromGenes(i);
        else
            if model.bitscore(num)<blastStructure.bitscore(i)
                model.bitscore(num)=blastStructure.bitscore(i);
                model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
                model.ppos(num)=blastStructure.ppos(i);
            end
        end
    end
end
model.genes=model.genes(1:num-1);
model.proteins=model.proteins(1:num-1);
model.bitscore=model.bitscore(1:num-1);
model.ppos=model.ppos(1:num-1);

% Get the indexes of matched genes in the metaCycModel
% because some enzymes in the FASTA file cannot be found in the dump file
[hits, indexes]=ismember(model.proteins,metaCycModel.genes);
found = find(hits);
model.genes=model.genes(found);

% Restructure the rxnGeneMat matrix
model.rxnGeneMat=sparse(metaCycModel.rxnGeneMat(:,indexes(found)));

%Remove all reactions without genes
hasGenes=any(model.rxnGeneMat,2);
model=removeReactions(model,~hasGenes,true);

%Generate grRules, only consider the or relationship here Matched enzymes
%are stored in field rxnScores,
rxnNum=numel(model.rxns);
model.rxnConfidenceScores=NaN(rxnNum,1);
model.rxnConfidenceScores(:)=2;
model.grRules=cell(rxnNum,1);
%model.rxnScores=cell(rxnNum,1);
for j=1:rxnNum
    model.grRules{j}='';
    I=find(model.rxnGeneMat(j,:));
    for k=1:numel(I)
        if isempty(model.grRules{j})
            model.grRules(j)=model.genes(I(k));
            %model.rxnScores(j)=model.proteins(I(k));
        else
            model.grRules(j)=strcat(model.grRules(j),{' or '},model.genes(I(k)));
            %model.rxnScores(j)=strcat(model.rxnScores(j),{' or
            %'},model.proteins(I(k)));
        end
    end
end
%update genes field
model.genes=model.genes(any(model.rxnGeneMat));

%Construct the S matrix and list of metabolites
[S, mets, badRxns]=constructS(model.equations);
model.S=S;
model.mets=mets;

%model=removeReactions(model,badRxns,true,true);

%Then match up metabolites
metaCycMets=getMetsFromMetaCyc([]);

%Add information about all metabolites to the model
[a, b]=ismember(model.mets,metaCycMets.mets);
a=find(a);
b=b(a);

if ~isfield(model,'metNames')
    model.metNames=cell(numel(model.mets),1);
    model.metNames(:)={''};
end
model.metNames(a)=metaCycMets.metNames(b);

if ~isfield(model,'metFormulas')
    model.metFormulas=cell(numel(model.mets),1);
    model.metFormulas(:)={''};
end
model.metFormulas(a)=metaCycMets.metFormulas(b);

if ~isfield(model,'metCharges')
    model.metCharges=zeros(numel(model.mets),1);
end
model.metCharges(a)=metaCycMets.metCharges(b);

if ~isfield(model,'b')
    model.b=zeros(numel(model.mets),1);
end
%model.b(a)=metaCycMets.b(b);

if ~isfield(model,'inchis')
    model.inchis=cell(numel(model.mets),1);
    model.inchis(:)={''};
end
model.inchis(a)=metaCycMets.inchis(b);

if ~isfield(model,'metMiriams')
    model.metMiriams=cell(numel(model.mets),1);
end
model.metMiriams(a)=metaCycMets.metMiriams(b);

%Put all metabolites in one compartment called 's' (for system). This is
%done just to be more compatible with the rest of the code
model.comps={'s'};
model.compNames={'System'};
model.metComps=ones(numel(model.mets),1);

%It could also be that the metabolite names are empty for some reason In
%that case, use the ID instead
I=cellfun(@isempty,model.metNames);
model.metNames(I)=model.mets(I);

%Remove additional fields
model=rmfield(model,{'proteins','bitscore','ppos'});

%In the end fix grRules and rxnGeneMat
[grRules,rxnGeneMat] = standardizeGrRules(model,false); %Get detailed output
model.grRules = grRules;
model.rxnGeneMat = rxnGeneMat;
%Remove the temp fasta file
delete(fastaFile)
end

author: SysBioChalmers

repo: RAVEN


text:
function a = reshape(a,varargin)
%RESHAPE      Reshape for slope vectors/matrices
%
%   r = reshape(a,vector)  or  r = reshape(a,n1,n2,...)
%
% functionality as Matlab function reshape
%

% written  09/28/01     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/12     S.M. Rump  performance improvement
%

  index = reshape(1:prod(a.size),varargin{:});
  a.size = size(index);

repo: OpenMAS


text:
function [ cluass,cluNRG ] = bz_GradDescCluster( simMat, varargin )
%[cluass] = bz_GradDescCluster(simMat) clusters recording sites given a pairwise
%similarity matrix using gradient descent to find minimize within-cluster
%interaction energy. (i.e. maximize within-cluster coherence, see Berenyi
%et al 2014 for details).
%
%INPUT
%   simMat          undirected similarity matrix (coherence, correlation, or other)
%
%   (optional parameters)
%   'numsteps'      number of steps (default: 5e5)
%   'numinit'       number of initial clusters, should be larger than the
%                   number of clusters you expect (default: 20)
%   'stopthresh'    stability threshold to stop descending. probability of
%                   selected site changing cluster identity, P(switch)
%                   (default: 0.001)
%   'stopwin'       window of time for which the P(switch) must stay below
%                   stopthresh (default: 10000 steps)
%   'showplot'      true or false, to show the update plot as it goes
%   'brainmap'      cell array with two cells: {1} = linear indices, {2} = size map; 2 dim vector x and y dimensions of map

%OUTPUT
%   cluass      cluster assignments for each
%   cluNRG      negative energy (i.e. mean within-cluster coherence) for each cluster
%
%From Berenyi et al. (2014). Large-scale, high-density (up to 512
%   channels) recording of local circuits in behaving animals. Journal of
%   Neurophysiology, 111(5), 1132?1149. http://doi.org/10.1152/jn.00785.2013
%
%Code implementation by DLevenstein 2016. Updates DL and RS May 2017.
%
% TO DO
% - make brainmap general to any 2D ephys/imaging data
% - develop method for stopping num iterations
% - add consensus clustering option?

%% Parse the input parameters
parms = inputParser;
addParameter(parms,'numsteps',500000,@isnumeric);
addParameter(parms,'numinit',20,@isnumeric);
addParameter(parms,'showplot',true,@islogical);
addParameter(parms,'brainmap',{},@iscell);
addParameter(parms,'stopthresh',0.001,@(x) x>0 && x<=1)
addParameter(parms,'stopwin',10000,@isnumeric)

parse(parms,varargin{:})
numsteps = parms.Results.numsteps;
numinit = parms.Results.numinit;
SHOWPLOT = parms.Results.showplot;
LinearInds = parms.Results.brainmap{1};
SizeVid = parms.Results.brainmap{2};
stopthresh = parms.Results.stopthresh;
stopwin = parms.Results.stopwin;

%% Initialize and Run
numsites = size(simMat,1);
rng('shuffle') %Shuffle the random number generator seed... just in case
cluass = randi(numinit,numsites,1); %Start from random initial assignments

if SHOWPLOT; figure; switchtracker = zeros(numsteps,1); end
for ss = 1:numsteps
    %The time ticker...
    if mod(ss,1000)==1
        display(['Step: ',num2str(ss),' of ',num2str(numsteps)])
        [~,clusort] = sort(cluass);

        % Visuals
        if SHOWPLOT
            % plot changing clusters over time
            subplot(2,2,1)
                imagesc(simMat(clusort,clusort))
            %plot switching behavior
            subplot(2,2,2)
                plot(log10(smooth(switchtracker(1:ss),1000,'moving')))
                hold on
                plot(get(gca,'xlim'),log10(stopthresh).*[1 1],'r--')
                xlabel('Step #');ylabel('P(switch)')
                ylim([log10(stopthresh)-1 0])
                LogScale('y',10)
                hold off
            % plot map
            if ~isempty(LinearInds)
                map = ExpandVid(cluass, LinearInds, SizeVid);
                subplot(2,2,3)
                imagesc(map)
            end;
            %plot number of clusters
            subplot(2,2,4)
                plot(ss,length(unique(cluass)),'ko') %number clusters
                hold on
                xlabel('Step #');ylabel('# Clusters')
                drawnow
        end;
    end;

    % Fx meat
    [cluass,switchtracker(ss)] = DescOneStep(cluass,simMat);

    %Decide to exit the loop - if you've been consistent for a certain window
    if mean(switchtracker(max(ss-stopwin,1):ss)) <= stopthresh
        disp('Clustering has descended to the bottom of the (local) pit!')
        break
    end
end

%Calculate -Energy (mean similarity) for each final cluster
finalclus = unique(cluass);
numfinalclus = length(finalclus);
for ff = 1:numfinalclus
    %N = sum(cluass == finalclus(ff));
    clusmat = simMat(cluass == finalclus(ff),cluass == finalclus(ff));
    cluspairs = triu(clusmat,1);
    E(ff) = mean(cluspairs(:));
end
cluNRG = [finalclus,E'];

%% FUNCTION: OneStep
    function [newcluass,diditswitch] = DescOneStep(oldcluass,simMat)
        clus = unique(oldcluass);
        numclus = length(clus);

        %Pick a random site and calculate the energy for it's current cluster
        sitepick = randi(length(oldcluass),1);
        clu_A = oldcluass(sitepick); %Cluster the current site is in
        othersites = setdiff(1:length(oldcluass),sitepick);
        N_A = sum(oldcluass == clu_A)-1; %ALl OTHER elements of current cluster
        E_A = -(1./N_A).*sum(simMat(sitepick,oldcluass(othersites)==clu_A));

        %If only one site in cluster, energy is 0
        if N_A == 0
            E_A = 0;
        end

        %Calculate energy gap for switching site to other clusters
        energygap = zeros(size(clus));
        for cc = 1:numclus
            clu_B = clus(cc);
            if clu_B==clu_A
                continue
            end
            N_B = sum(oldcluass == clu_B);
            E_B = -(1./N_B).*(sum(simMat(sitepick,oldcluass==clu_B)));
            energygap(cc) = E_A - E_B;
        end

        %Switch site to cluster with largest energy gap
        newcluass = oldcluass;
        newcluass(sitepick) = clus(energygap==max(energygap));

        %Keep Track of switching - 1 if the channel stays in its cluster
        diditswitch = clus(energygap==max(energygap))~=clu_A;
    end
end

%% FUNCTION: ExpandVid

function map = ExpandVid(cluass, LinearInds, SizeVid)

    % Intialize
    map = single(nan(SizeVid(1)*SizeVid(2),1));

    % Put tmpR values int initialized 2D matrix
    map(LinearInds) = cluass;

    % Reshape into video
    map = reshape(map,SizeVid(1),SizeVid(2));

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ p4, flag ] = line_exp_perp_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% LINE_EXP_PERP_2D computes a line perpendicular to a line and through a point.
%
%  Discussion:
%
%    The explicit form of a line in 2D is:
%
%      ( P1, P2 ) = ( (X1,Y1), (X2,Y2) ).
%
%    The input point P3 should NOT lie on the line (P1,P2).  If it
%    does, then the output value P4 will equal P3.
%
%    P1-----P4-----------P2
%            |
%            |
%           P3
%
%    P4 is also the nearest point on the line (P1,P2) to the point P3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), two points on the line.
%
%    Input, real P3(2,1), a point (presumably not on the
%    line (P1,P2)), through which the perpendicular must pass.
%
%    Output, real P4(2,1), a point on the line (P1,P2),
%    such that the line (P3,P4) is perpendicular to the line (P1,P2).
%
%    Output, logical FLAG, is TRUE if the value could not be computed.
%
  dim_num = 2;

  bot = sum ( ( p2(1:dim_num,1) - p1(1:dim_num,1) ).^2 );

  if ( bot == 0.0 )
    p4(1:2,1) = Inf;
    flag = 1;
  end
%
%  (P3-P1) dot (P2-P1) = Norm(P3-P1) * Norm(P2-P1) * Cos(Theta).
%
%  (P3-P1) dot (P2-P1) / Norm(P3-P1)**2 = normalized coordinate T
%  of the projection of (P3-P1) onto (P2-P1).
%
  t = sum ( ( p1(1:dim_num,1) - p3(1:dim_num,1) ) ...
    .* ( p1(1:dim_num,1) - p2(1:dim_num,1) ) ) / bot;

  p4(1:dim_num,1) = p1(1:dim_num,1) + t * ( p2(1:dim_num,1) - p1(1:dim_num,1) );
  flag = 0;

  return
end


text:
function r_size  = nwspgr_size ( type, dim, k, sym, compress )

%*****************************************************************************80
%
%% NWSPGR_SIZE determines the size of a sparse grid rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Florian Heiss, Viktor Winschel,
%    Likelihood approximation by numerical integration on sparse grids,
%    Journal of Econometrics,
%    Volume 144, 2008, pages 62-80.
%
%  Parameters:
%
%    Input, string TYPE, selects the 1D integration rule:
%    "kpu": Nested rule for unweighted integral over [0,1];
%    "kpn": Nested rule for integral with Gaussian weight;
%    "gqu": Gaussian quadrature for unweighted integral over [0,1]
%    (Gauss-Legendre);
%    "gqn": Gaussian quadrature for integral with Gaussian weight
%    (Gauss-Hermite);
%    func:  any function name. Function must accept level l and return nodes N
%    and weights W for univariate quadrature rule with polynomial exactness
%    2*L-1 as [n, w] = feval(func,level)
%
%    Input, integer DIM, the dimension of the integration problem.
%
%    Input, integer K, the level.  When using the built in 1D rules, the
%    resulting sparse grid will be exact for polynomials up to total order
%    2*K-1.  When using the built-in quadrature rules, the maximum value of K
%    that is available is 25.
%
%    Input, integer SYM.  This is only useful when the user has
%    chosen to supply the 1D quadrature rule.  In that case, if the rule is
%    symmetric, specifying SYM to be 1 will allow the code to run faster.
%    But this also requires that the user quadrature rule function only
%    return the nonnegative abscissas and their weights.
%
%    Input, integer COMPRESS,
%    0, do not compress the rule ( = do not merge duplicate points.)
%    1, compress the rule.
%
%    Output, integer R_SIZE, the "size" of the rule.
%

%
%  Interpret the inputs.  Use STRCMPI so that case is ignored.
%
  builtinfct = ( strcmpi ( type, 'gqu' ) || ...
                 strcmpi ( type, 'gqn' ) || ...
                 strcmpi ( type, 'kpu' ) || ...
                 strcmpi ( type, 'kpn' ) );

  sym = logical ( sym );
%
%  Retrieve the 1D rules of levels 1 to K.
%  Create cell arrays X1D and W1D containing the node and weight information.
%  The array N1D stores the length or order of each rule.
%
%  The built in functions already return only the positive orthant.
%
%  If the rule is being supplied by the user, and the user indicates that
%  it's a symmetric rule, then similarly keep only the positive orthant.
%
%  The result will be that if SYM is true, whether for built in or user-supplied
%  functions, we have only the positive orthant data stored.
%
  n1d = zeros ( 1, k );
  x1d = cell ( k, 1 );
  w1d = cell ( k, 1 );

  for level = 1 : k

    [ x, w ] = feval ( type, level );

    if ( ( ~ builtinfct ) && sym )
      [ numnew, dummy ] = size ( x );
      [ x, sortvec ] = sortrows ( x );
      w = w ( sortvec );
      x = x((floor(numnew/2)+1):numnew,:);
      w = w((floor(numnew/2)+1):numnew,:);
    end

    n1d(level) = length ( w );
    x1d{level} = x;
    w1d{level} = w;

  end
%
%  Initialization.
%
  minq = max ( 0, k - dim );
  maxq = k - 1;
  nodes = [];
  weights = [];
%
%  Loop for max ( 0, K - DIM ) <= Q <= K - 1.
%
  for q = minq : maxq

    r = length ( weights );
%
%  BQ is the combinatorial coefficient applied to the component
%  product rules which have level Q.
%
    bq = ( -1 )^( maxq - q ) * nchoosek ( dim - 1, dim + q - k );
%
%  Compute the D-dimensional row vectors that sum to DIM+Q.
%
    is = get_seq ( dim, dim + q );
%
%  Preallocate new rows for nodes and weights.
%
    Rq = prod ( n1d(is), 2 );
    sRq = sum ( Rq );
    nodes   = [ nodes;   zeros(sRq,dim) ];
    weights = [ weights; zeros(sRq,1) ];
%
%  Generate each of the product rules indicated by IS, and
%  insert them into NODES and WEIGHTS.
%
    for j = 1 : size(is,1)
      midx = is(j,:);
      [ newn, neww ] = tensor_product ( x1d(midx), w1d(midx) );
      nodes((r+1):(r+Rq(j)),:) = newn;
      weights((r+1):(r+Rq(j))) = bq .* neww;
      r = r + Rq(j);
    end
%
%  Sort the nodes and merge repeated values.
%
    [ nodes, sortvec ] = sortrows ( nodes );
    weights = weights(sortvec);
    keep = 1;
    lastkeep = 1;

    for j = 2 : size(nodes,1)
      if ( compress )

        if ( nodes(j,:) == nodes(j-1,:) )
          weights(lastkeep) = weights(lastkeep) + weights(j);
        else
          lastkeep = j;
          keep = [ keep ; j ];
        end
      else
        lastkeep = j;
        keep = [ keep ; j ];
      end
    end

    nodes = nodes(keep,:);
    weights = weights(keep);

  end
%
%  The rule has been computed.
%  If we used symmetry, we now have to extend the rule.
%
  if ( sym )

    nr = length ( weights );
    m = x1d{1};

    for j = 1 : dim

      keep = zeros(nr,1);
      numnew = 0;

      for r = 1 : nr
        if ( nodes(r,j) ~= m )
          numnew = numnew + 1;
          keep(numnew) = r;
        end
      end

      if ( 0 < numnew )
        nodes = [nodes ; nodes(keep(1:numnew),:)];
        nodes(nr+1:nr+numnew,j) = 2*m - nodes(nr+1:nr+numnew,j);
        weights = [weights ; weights(keep(1:numnew))];
        nr = nr + numnew;
      end

    end

    [ nodes, sortvec ] = sortrows ( nodes );
    weights = weights(sortvec);

  end

  r_size = length ( weights );

  return
end


text:
function y = conFun(x)
y = [x(1)-1;x(2)];

repo: Optimization-Theory-and-Algorithm


text:
function [numDataSym] = offline(cfgNonHT)

mcsTable = wlan.internal.getRateTable(cfgNonHT);
Ntail = 6; Nservice = 16;
numDataSym = ceil((8*cfgNonHT.PSDULength + Nservice + Ntail)/mcsTable.NDBPS);

author: analogdevicesinc


text:
function writefile(DS, filename)
%WRITEFILE Writes a dataset to a DSDML file
%
% $ Syntax $
%   - writefile(DS, filename)
%
% $ Arguments $
%   - DS:           the dataset object
%   - filename:     the name of the file to be written to
%
% $ Description $
%   - writefile(DS, filename) writes a dataset object to the DSDML file
%     with the filename specified.
%
% $ History $
%   - Created by Dahua Lin on Jul 23rd, 2006
%

import com.mathworks.xml.*;

%% Initialize the document

xdoc = XMLUtils.createDocument('DataSet');
xdoc.setVersion('1.0');
xdoc.setEncoding('UTF-8');

S = struct(DS);

%% Set headers

docelem = xdoc.getDocumentElement;
docelem.setAttribute('version', char(S.version));
docelem.setAttribute('name', char(S.name));
docelem.setAttribute('unit', char(S.unittype));
docelem.setAttribute('format', char(S.format));
if ~isempty(S.author)
    docelem.setAttribute('author', char(S.author));
end
if ~isempty(S.description)
    docelem.setAttribute('description', char(S.description));
end
if ~isempty(S.attribs)
    attrnames = fieldnames(S.attribs);
    nattrs = length(attrnames);
    for i = 1 : nattrs
        docelem.setAttribute(attrnames{i}, S.attribs.(attrnames{i}));
    end
end

%% Create unit nodes

switch S.unittype
    case 'Sample'
        n = length(S.units);
        for i = 1 : n
            docelem.appendChild(CreateSampleNode(xdoc, S.units(i)));
        end

    case 'SampleGroup'
        n = length(S.units);
        for i = 1 : n
            docelem.appendChild(CreateSampleGroupNode(xdoc, S.units(i)));
        end

    otherwise
        error('dsdml:invalidutype', ...
            'Invalid unit type %s', DS.unittype);
end

%% Write document
xmlwrite(filename, xdoc);

%% Node Creating functions

function sampleNode = CreateSampleNode(xdoc, sample)

sampleNode = xdoc.createElement('Sample');
sampleNode.setAttribute('class_id', int2str(sample.class_id));
if ~isempty(sample.filename)
    sampleNode.setAttribute('filename', sample.filename);
end
if ~isempty(sample.attribs)
    attrnames = fieldnames(sample.attribs);
    nattrs = length(attrnames);
    for i = 1 : nattrs
        sampleNode.setAttribute(attrnames{i}, sample.attribs.(attrnames{i}));
    end
end

function groupNode = CreateSampleGroupNode(xdoc, grp)

groupNode = xdoc.createElement('SampleGroup');
groupNode.setAttribute('class_id', int2str(grp.class_id));
if ~isempty(grp.attribs)
    attrnames = fieldnames(grp.attribs);
    nattrs = length(attrnames);
    for i = 1 : nattrs
        sampleNode.setAttribute(attrnames{i}, grp.attribs.(attrnames{i}));
    end
end
nsamples = length(grp.samples);
for i = 1 : nsamples
    groupNode.appendChild(CreateSampleNode(xdoc, grp.samples(i)));
end

author: lmthang


text:
function [Gamma,vpath,error,Ypred] = tudadecode(X,Y,T,tuda,new_experiment,parallel_trials)
% Having estimated the TUDA model (i.e. the corresponding decoding models)
% in the same or a different data set, this function finds the model time
% courses (with no re-estimation of the decoding parameters)
%
% INPUT
% X: Brain data, (time by regions)
% Y: Stimulus, (time by q); q is no. of stimulus features
%               For binary classification problems, Y is (time by 1) and
%               has values -1 or 1
%               For multiclass classification problems, Y is (time by classes)
%               with indicators values taking 0 or 1.
%           If the stimulus is the same for all trials, Y can have as many
%           rows as trials, e.g. (trials by q)
% T: Length of series or trials
% tuda: Estimated TUDA model, using tudatrain
% new_experiment: Whether or not the estimated model is going to be applied
%   on data that follows the same paradigm used to train the model. If the
%   paradigm changes, then new_experiment should be 1. For example, that
%   would be the case if we train the model on perception and test it on
%   recalling.
% parallel_trials: if set to 1, then
%   all trials have the same experimental design and that the
%   time points correspond between trials; in this case, all trials
%   must have the same length. If set to 0, then there is not a fixed
%   experimental design for all trials.
%
% OUTPUT
% Gamma: Time courses of the states (decoding models) probabilities given data
% vpath: Most likely state path of hard assignments
% error: Error for each state on the new data;
%        if parallel_trials = 1, then error is
%          (trial tim by no. of stimuli by no. of states);
%        otherwise, this is
%          (total time by no. of stimuli by no. of states).
% Ypred: Predicted stimulus at each time point
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2018)

if nargin < 5, new_experiment = 0; end
if nargin < 6, parallel_trials = 0; end

max_num_classes = 5;
do_preproc = 1;

N = length(T); q = size(Y,2); ttrial = T(1); p = size(X,2); K = tuda.train.K;

% Check options and put data in the right format
tuda.train.parallel_trials = 0;
if isfield(tuda.train,'orders')
    orders = tuda.train.orders;
    tuda.train = rmfield(tuda.train,'orders');
else
    orders = [];
end
if isfield(tuda.train,'active')
    active = tuda.train.active;
    tuda.train = rmfield(tuda.train,'active');
else
    active = [];
end

classification = length(unique(Y(:))) < max_num_classes;
if classification
    Ycopy = Y;
    if size(Ycopy,1) == N
        Ycopy = repmat(reshape(Ycopy,[1 N q]),[ttrial 1 1]);
    end
    % no demeaning by default if this is a classification problem
    if ~isfield(tuda.train,'demeanstim'), tuda.train.demeanstim = 0; end
end

if do_preproc
    if isfield(tuda.train,'embeddedlags'), el = tuda.train.embeddedlags; end
    tuda.train.intercept = 0;
    [X,Y,T,options] = preproc4hmm(X,Y,T,tuda.train); % this demeans Y
    p = size(X,2);
    if classification && length(el) > 1
        Ycopy = reshape(Ycopy,[ttrial N q]);
        Ycopy = Ycopy(-el(1)+1:end-el(end),:,:);
        Ycopy = reshape(Ycopy,[T(1)*N q]);
    end
end
if ~isempty(active), tuda.train.active = active; end
if ~isempty(orders),  tuda.train.orders = orders;  end

if isfield(options,'A') % it is done in preproc4hmm
    options = rmfield(options,'A');
end
if isfield(options,'parallel_trials')
    options = rmfield(options,'parallel_trials');
end
if isfield(options,'add_noise')
    options = rmfield(options,'add_noise');
end

% Put X and Y together
Ttmp = T;
T = T + 1;
Z = zeros(sum(T),q+p,'single');
for j = 1:N
    t1 = (1:T(j)) + sum(T(1:j-1));
    t2 = (1:Ttmp(j)) + sum(Ttmp(1:j-1));
    if strcmp(options.classifier,'LDA') || (isfield(options,'encodemodel') && options.encodemodel)
        Z(t1(2:end),1:p) = X(t2,:);
        Z(t1(1:end-1),(p+1):end) = Y(t2,:);
    else
        Z(t1(1:end-1),1:p) = X(t2,:);
        Z(t1(2:end),(p+1):end) = Y(t2,:);
    end
end

if new_experiment
    off_diagonal = [tuda.P(triu(true(K),1)); tuda.P(tril(true(K),-1))];
    in_diagonal = tuda.P(eye(K)==1);
    tuda.P(triu(true(K),1)) = mean(off_diagonal);
    tuda.P(tril(true(K),-1)) = mean(off_diagonal);
    tuda.P(eye(K)==1) = mean(in_diagonal);
    tuda.Pi(:) = mean(tuda.Pi);
end

% Run TUDA inference
options.S = -ones(p+q);
if strcmp(options.classifier,'LDA') || (isfield(options,'encodemodel') && options.encodemodel)
    options.S(p+1:end,1:p) = 1;
else
    options.S(1:p,p+1:end) = 1;
end
options.updateObs = 0;
options.updateGamma = 1;
options.updateP = 1;
options.hmm = tuda;
options.repetitions = 0;
options.pca = 0;
options.cyc = 1;
options.tuda = 0;
[~,Gamma,~,vpath] = hmmmar(Z,T,options);
T = T - 1;

if parallel_trials
    if classification, error = zeros(max(T),K);
    else, error = zeros(max(T),q,K);
    end
else
    if classification, error = zeros(size(Gamma,1),K);
    else, error = zeros(size(Gamma,1),q,K);
    end
end
Betas = tudabeta(tuda);
if nargout>2
    for k = 1:K
      Ypred = X * Betas(:,:,k);
      if classification
          Ypred = continuous_prediction_2class(Ycopy,Ypred);
          Y = continuous_prediction_2class(Ycopy,Y);
          if q == 1
              e = abs(Y - Ypred) < 1e-4;
          else
              e = sum(abs(Y - Ypred),2) < 1e-4;
          end
      else
          e = (Y - Ypred).^2;
      end
      if parallel_trials
          maxT = max(T);
          me = zeros(maxT,1);
          ntrials = zeros(maxT,1);
          for j = 1:N
              t0 = sum(T(1:j-1));
              ind_1 = (1:T(j)) + t0;
              ind_2 = 1:length(ind_1);
              me(ind_2,:) = me(ind_2,:) + e(ind_1,:);
              ntrials(ind_2) = ntrials(ind_2) + 1;
          end
          e = me ./ ntrials;
      end
      if classification, error(:,k) = e;
      else, error(:,:,k) = e;
      end
    end

    error = squeeze(error);
end

end

author: OHBA-analysis

repo: HMM-MAR


function varargout = test_edit2(varargin)
% TEST_EDIT2 M-file for test_edit2.fig
%      TEST_EDIT2, by itself, creates a new TEST_EDIT2 or raises the existing
%      singleton*.
%
%      H = TEST_EDIT2 returns the handle to a new TEST_EDIT2 or the handle to
%      the existing singleton*.
%
%      TEST_EDIT2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in TEST_EDIT2.M with the given input arguments.
%
%      TEST_EDIT2('Property','Value',...) creates a new TEST_EDIT2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before test_edit2_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to test_edit2_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help test_edit2

% Last Modified by GUIDE v2.5 13-Mar-2004 14:59:53

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @test_edit2_OpeningFcn, ...
                   'gui_OutputFcn',  @test_edit2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before test_edit2 is made visible.
function test_edit2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to test_edit2 (see VARARGIN)

% Choose default command line output for test_edit2
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes test_edit2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = test_edit2_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function EditBox2_Callback(hObject, eventdata, handles)
% hObject    handle to EditBox2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of EditBox2 as text
%        str2double(get(hObject,'String')) returns contents of EditBox2 as a double

% Find the value typed into the edit box
str = get (handles.EditBox2,'String');

% Place the value into the text field
set (handles.TextBox2,'String',str);

% --- Executes during object creation, after setting all properties.
function EditBox2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to EditBox2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


text:
function [V,F] = triangle_strip(P,E)
  % TRIANGLE_STRIP build a strip of triangles given a curve
  %
  % [V,F] = triangle_strip(P,E)
  %
  % Inputs:
  %   P  #points by 2 list of 2d points
  %   E  #edges by 2 list of edge indices into P
  % Outputs:
  %   V  #points*2 by 3 list of 3d points
  %   F  #edges*2 by 3 list of triangle indices
  %
  n = size(P,1);
  V = [P -ones(n,1);P ones(n,1)];
  F = [E n+E(:,1);n+fliplr(E) E(:,2)];
end

author: alecjacobson

repo: gptoolbox


text:
function f = p24_f ( n, x )

%*****************************************************************************80
%
%% P24_F evaluates the objective function for problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = ( sum ( floor ( x(1:n) ) ) );

  return
end


text:
function [ xpnew ] = syncFlippedLandmarks( xp )

lookup = [6 5 4 3 2 1 12 11 10 9 8 7 17 16 15 14 13 20 19 18 21];
if size(xp,1)==2
    xpnew(1,:) = xp(1,lookup);
    xpnew(2,:) = xp(2,lookup);
else
    xpnew = xp(lookup);
end

end

author: anilbas


text:
function [ a, b ] = p40_lim ( )

%*****************************************************************************80
%
%% P40_LIM returns the integration limits for problem 40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
function [influenceM] = getInfluenceM(IM, structNum)
%"getInfluenceM"
%   Get the influence matrix (dose(:) = influenceM * weights(:)),
%   for a single structure.  The output influenceM is sparse, and has
%   a row for every calculated dose point in the structure.  If a
%   the sampleRate is 1, this is the number of points in the structure,
%   otherwise it is the number of points in the downsampled structure.
%
%   Uses the global influence matrix.  This is not the most efficent
%   system, should be changed later.
%
%   JRA 4/1/04
%
%Usage:
%   function [influenceM] = getInfluenceM(IM, structNum)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Get global influence for all points in plan.
gInfluenceM = getGlobalInfluenceM(IM, structNum);

%Get the structure mask, downsampled if necessary.
mask3M = getUniformStr(structNum);

%get indices of structures stored under beamlets
beamlets = IM.beams(1).beamlets;
structIndV = getAssociatedStr({beamlets(:,1).strUID});
structInd = find(ismember(structIndV,structNum));
sampleRate = beamlets(structInd(1),1).sampleRate;
%sampleRate = IM.beamlets(structNum,1).sampleRate;

mask3M = getDown3Mask(mask3M, sampleRate, 1) & mask3M;

%Find relevant indices into the gInfluenceM.
indV = find(mask3M);

%Extract the influence matrix for this structure.
influenceM = gInfluenceM(indV, :);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest, fBvalue, fAvalue] = calc_McCdflognormal(mCatalog, fBinning)
% function [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest, fBvalue, fAvalue] = calc_McCdflognormal(mCatalog, fBinning);
% -----------------------------------------------------------------------------------------------------------------------------
% Determine Mc using maximum likelihood score
% Fitting non-cumulative frequency magnitude distribution above and below Mc:
% below: Cumulative LOGNORMAL density function
% above: Gutenberg-Richter law
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Binning interval, usually 0.1
%
% Outgoing variables:
% mResult     : Solution matrix including
%               vProbability: maximum likelihood score
%               vMc         : Mc values
%               vX_res      : mu (of normal CDF), sigma (of normal CDF), residuum, exitflag
%               vNmaxBest   : Number of events in lowest magnitude bin considered complete
%               vABValue    : a and b-value
% fMls       : minimum maximum likelihood score --> best Mc
% fMc        : Best estimated magnitude of completeness
% mDatPredBest   : Matrix of non-cumulative FMD [Prediction, magnitudes, original distribution]
% vPredBest      : Matrix of non-cumulative FMD below Mc [magnitude, prediction, uncertainty of prediction]
% fBvalue        : b-value
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 03.11.03

% Initialize
vProbability = [];
vMc = [];
vABValue =[];
mFitRes = [];
vX_res = [];
vNCumTmp = [];
mDataPred = [];
vPredBest = [];
vDeltaBest = [];
vX_res = [];
vNmaxBest = [];
mResult=[];
mDatPredBest = [];

% Determine exact time period
fPeriod1 = max(mCatalog(:,3)) - min(mCatalog(:,3));

% Determine max. and min. magnitude
fMaxMag = ceil(10 * max(mCatalog(:,6))) / 10;

% Set starting value for Mc loop and LSQ fitting procedure
fMcTry= calc_Mc(mCatalog,1);
fSmu = abs(fMcTry/2);
fSSigma = abs(fMcTry/4);
if (fSmu > 1)
    fSmu = fMcTry/10;
    fSSigma = fMcTry/20;
end
fMcBound = fMcTry;

% Calculate FMD for original catalog
[vFMDorg, vNonCFMDorg] = calc_FMD(mCatalog);
fMinMag = min(vNonCFMDorg(1,:));

%% Shift to positive values
% if fMinMag ~= 0
%     fMcBound = fMcTry-fMinMag;
% end
% Loop over Mc-values
for fMc = fMcBound-0.6:0.1:fMcBound+0.6
    fMc = round(fMc*10)/10;
    vFMD = vFMDorg;
    vNonCFMD = vNonCFMDorg;
    vNonCFMD = fliplr(vNonCFMD);
    % Calculate a and b-value for GR-law and distribution vNCum
    [nIndexLo, fMagHi, vSel, vMagnitudes] = fMagToFitBValue(mCatalog, vFMD, fMc);
    if (length(mCatalog(vSel,1)) >= 20)
        [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog(vSel,:), fBinning);
        % Normalize to time period
        vFMD(2,:) = vFMD(2,:)./fPeriod1; % ceil taken out
        vNonCFMD(2,:) = vNonCFMD(2,:)./fPeriod1; % ceil removed
        % Compute quantity of earthquakes by power law
        fMaxMagFMD = max(vNonCFMD(1,:));
        fMinMagFMD = min(vNonCFMD(1,:));
        vMstep = [fMinMagFMD:0.1:fMaxMagFMD];
        vNCum = 10.^(fAValue-fBValue.*vMstep); % Cumulative number

        % Compute non-cumulative numbers vN
        fNCumTmp = 10^(fAValue-fBValue*(fMaxMagFMD+0.1));
        vNCumTmp  = [vNCum fNCumTmp ];
        vN = abs(diff(vNCumTmp));

        % Normalize vN
        vN = vN./fPeriod1;
        % Data selection
        % mData = Non-cumulative FMD values from GR-law and original data
        mData = [vN' vNonCFMD'];
        vSel = (mData(:,2) >= fMc);
        mDataTest = mData(~vSel,:);
        mDataTmp = mData.subset(vSel);
        % Choices of normalization
        fNmax = mDataTmp(1,3); % Frequency of events in Mc bin
        %fNmax = max(mDataTest(:,3));  % Use maximum frequency of events in bins below Mc
        %fNmax = mDataTest(length(mDataTest(:,1)),3); % Use frequency of events at bin Mc-0.1 -> best fit
        if (~isempty(fNmax) & ~isnan(fNmax) & fNmax ~= 0 & length(mDataTest(:,1)) > 4)
            mDataTest(:,3) = mDataTest(:,3)/fNmax; % Normalize datavalues for fitting with CDF
            % Move to M=0 to fit with lsq-algorithm
            fMinMagTmp = min(mDataTest(:,2));
            mDataTest(:,2) = mDataTest(:,2)-fMinMagTmp;
            % Curve fitting: Non cumulative part below Mc
            options = optimset;
            %options = optimset('Display','off','Tolfun',1e-7,'TolX',0.0001,'MaxFunEvals', 100000,'MaxIter',10000);
            options = optimset('Display','off','Tolfun',1e-5,'TolX',0.001,'MaxFunEvals', 1000,'MaxIter',1000);
            [vX, resnorm, resid, exitflag, output, lambda, jacobian]=lsqcurvefit(@calc_lognormal,[fSmu  fSSigma], mDataTest(:,2), mDataTest(:,3),[],[],options);
            mDataTest(:,1) = logncdf(mDataTest(:,2), vX(1), vX(2))*fNmax;
            if (length(mDataTest(:,2)) > length(vX(1,:)))
                %% Confidence interval determination
                % vPred : Predicted values of lognormal function
                % vPred+-delta : 95% confidence level of true values
                [vPred,delta] = nlpredci(@calc_lognormal,mDataTest(:,2),vX, resid, jacobian);
            else
                vPred = NaN;
                delta = NaN;
            end; % END: This section is due for errors produced with datasets less long than amount of parameters in vX
            % Results of fitting procedure
            mFitRes = [mFitRes; vX resnorm exitflag];
            % Move back to original magnitudes
            mDataTest(:,2) = mDataTest(:,2)+fMinMagTmp;
            % Set data together
            mDataTest(:,3) = mDataTest(:,3)*fNmax;
            mDataPred = [mDataTest; mDataTmp];
            % Denormalize to calculate probabilities
            mDataPred(:,1) = round(mDataPred(:,1).*fPeriod1);
            mDataPred(:,3) = mDataPred(:,3).*fPeriod1;
            vProb_ = calc_log10poisspdf2(mDataPred(:,3), mDataPred(:,1)); % Non-cumulative

            % Sum the probabilities
            fProbability = (-1) * sum(vProb_);
            vProbability = [vProbability; fProbability];
            % Move magnitude back
            mDataPred(:,2) = mDataPred(:,2)+fMinMag;
            vMc = [vMc; fMc];
            vABValue = [vABValue; fAValue fBValue];

            % Keep values
            vDeltaBest = [vDeltaBest; delta];
            vX_res = [vX_res; vX resnorm exitflag];
            vNmaxBest = [vNmaxBest; fNmax];

            % Keep best fitting model
            if (fProbability == min(vProbability))
                vDeltaBest = delta;
                vPredBest = [mDataTest(:,2) vPred*fNmax*fPeriod1 delta*fNmax*fPeriod1]; % Gives back uncertainty
                %fMc+fMinMag : Test procedure
                mDatPredBest = [mDataPred];
           end
        else
            %disp('Not enough data');
            % Setting values
            fProbability = NaN;
            fMc = NaN;
            vX(1) = NaN;
            vX(2) = NaN;
            resnorm = NaN;
            exitflag = NaN;
            delta = NaN;
            vPred = [NaN NaN NaN];
            fNmax = NaN;
            fAValue = NaN;
            fBValue = NaN;
            vProbability = [vProbability; fProbability];
            vMc = [vMc; fMc];
            vX_res = [vX_res; vX resnorm exitflag];
%             vDeltaBest = [vDeltaBest; NaN];
%             vPredBest = [vPredBest; NaN NaN NaN];
            vNmaxBest = [vNmaxBest; fNmax];
            vABValue = [vABValue; fAValue fBValue];
        end; % END of IF fNmax
    end; % END of IF length(mCatalog(vSel,1))

    % Clear variables
    vNCumTmp = [];
    mModelDat = [];
    vNCum = [];
    vSel = [];
    mDataTest = [];
    mDataPred = [];
end; % END of FOR fMc
% Result matrix
mResult = [mResult; vProbability vMc vX_res vNmaxBest vABValue];

% Find best estimate, excluding the case of mResult all NAN
if  ~isempty(nan(mResult))
    if ~isnan(nan(mResult(:,1)))
        vSel = find(nan(mResult(:,1)) == mResult(:,1));
        fMc = min(mResult(vSel,2));
        fMls = min(mResult(vSel,1));
        fMu = min(mResult(vSel,3));
        fSigma = min(mResult(vSel,4));
        fAvalue = min(mResult(vSel,8));
        fBvalue = min(mResult(vSel,9));
    else
        fMc = NaN;
        fMls = NaN;
        fMu = NaN;
        fSigma = NaN;
        fAvalue = NaN;
        fBvalue = NaN;
    end
else
    fMc = NaN;
    fMls = NaN;
    fMu = NaN;
    fSigma = NaN;
    fAvalue = NaN;
    fBvalue = NaN;
end

author: CelsoReyes


text:
function mtrEC2Tensor(subjectDirName, numRepeats, diffusionTime, numBootstrap)
% Converts eddy corrected DWI to tensor and related files.
%
%  mtrEC2Tensor(subjectDirName, numRepeats, [diffusionTime], [numBootstrap])
%
% subjectDirName: This conversion process follows a standardized file
% naming procedure. With files stored in the subject folder.
% numRepeats: Number of repeated volumes of DWIs recorded.
% numDWGrads: Number of diffusion weighted and non-diffusion weighted
%   measurements taken per repeat.
% diffusionTime: Time in seconds that diffusion was allowed to occur during
%   DWI experiment.  Default is 0.04s.
% numBootstrap: Number of synthetic bootstrap samples to generate during
%   tensor fitting and scanner noise estimation.  Default is 1000.
%
% With the following files and directories expected to already exist:
% INPUT
% rawAlignedFile - 'rawAligned.nii.gz': Eddy current corrected DWI images in compressed Nifti format
% bvalsFile - 'bvals': FSL format files specifying the b (s/mm^2) value of each DW measurement, where
%   b = q^2*delta and q = labeling strength (1/mm) and delta (s) is the diffusion time which we
%   assume as a default is 0.04s
% bvecsFile - 'rawAligned.bvecs': FSL format files specifying the vectors for DW measurements.
%   These vectors have a length of one.
%
%
% And the following files and directories produced:
% OUTPUT
% schemeFile - 'camino.scheme': Scheme file contains the same information
%   as the rawAligned.bvecs and bvals file combined with the additional diffusionTime
%   information.  This file is used during the bootstrap estimation of
%   scanner noise.
% rawAlignedVoxelFormatFile - 'rawDtiAligned.Bfloat': Big endian voxel order data file.
%   This file contains the same information as rawAlignedFile, except that the
%   data is ordered so that all voxel information is stored in the file in
%   order.
% bootstrapData - 'dtboot.Bdouble': Result of running the modelfit code
%   which uses bootstrap to generate 1000 sample
% brainMaskFile - 'brainMask.nii.gz': Output from FSL's BET program that
%   is a binary.
% tensorFile - 'bin\tensors.nii.gz': Tensor file in the format for DTIQuery.
%   Dxx,Dyy,Dzz,Dxy,Dxz,Dyz.
% meanB0File - 'bin\b0.nii.gz': Mean of the B0 data from the
%  raw eddy corrected files.
% faFile - 'bin\fa.nii.gz': Fractional anisotropy image
%   calculated from the tensorFile.
% wmMaskFile - 'bin\wmMask.nii.gz': Binary image indicating non-zero value
%   where white matter tissue is expected and zero otherwise.
% pdfFile - 'bin\pdf.nii.gz': Probability distribution function file that
%  is used during tractography in order to represent scanner noise per
%  voxel as well as shape information of the tensor fit per voxel. Data
%  stored: EV3,EV2,EV1,k1,k2,Cl,Cp.

%
% AUTOMATICALLY GENERATED OUTPUT
%
% SPECIFIED OUTPUT
%
% Examples:
%   Generate dt6 using dtiMakeDt6 <- (high level calls FromBin, FromTensorCalc, Muscle)
%
% AJS
%

% Go to subject directory
oldDirName = pwd;
cd(subjectDirName);

if( ieNotDefined('diffusionTime') )
    diffusionTime = 0.04;
end

if( ieNotDefined('numBootstrap') )
    numBootstrap = 1000;
end

schemeFile = fullfile('raw','camino.scheme');
rawAlignedFile = fullfile('raw','rawDti_aligned.nii.gz');
rawAlignedVoxelFormatFile = fullfile('raw','rawDti_aligned.Bfloat');
rawBvecsFile = fullfile('raw','rawDti_aligned.bvecs');
rawBvalsFile = fullfile('raw','rawDti_aligned.bvals');

% Check to see if which files need to be generated
%bGenBoot = queryOverwrite('dtboot.Bdouble');
bGenBoot = queryOverwrite(fullfile('bin','pdf.nii.gz'));
if bGenBoot
    bGenScheme = queryOverwrite(schemeFile);
else
    bGenScheme = 0;
end

% Convert bvecs and bvals into Camino scheme file
if bGenScheme
    if ispc
        error('Can not call fsl2scheme from matlab on Windows OS.')
    end
    delete(schemeFile);
    cmd = sprintf('fsl2scheme -bvalfile %s -bvecfile %s -diffusiontime %g -outputfile %s', rawBvalsFile, rawBvecsFile, diffusionTime, schemeFile);
    disp(cmd);
    [s, ret_info] = system(cmd,'-echo');
else
    disp('Using previously generated camino.scheme');
end

% Convert NIFTI raw image into Camino voxel format
if bGenBoot
    disp('Generating voxel format file, rawAligned.Bfloat');
    rawAligned = mtrNiftiToCamino(rawAlignedFile,'float', rawAlignedVoxelFormatFile);
    xformToAcPc = rawAligned.qto_xyz;
    mmPerVox = rawAligned.pixdim(1:3);
    numDWGrads = size(rawAligned.data,4) / numRepeats;
    if floor(numDWGrads) ~= numDWGrads
        msg = ['Num. volumes = ', num2str(size(rawAligned.data,4)), ' Num. repeats can''t be ', num2str(numRepeats)];
        error(msg);
    end
    % Calculate and save average B0 data
    img_b0 = mean(rawAligned.data(:,:,:,1:numDWGrads:end),4);
    dtiWriteNiftiWrapper(img_b0,xformToAcPc,fullfile('bin','b0.nii.gz'));
    clear rawAligned;
else
    nib0 = niftiRead(fullfile('bin','b0.nii.gz'));
    img_b0 = nib0.data;
    xformToAcPc = nib0.qto_xyz;
    mmPerVox = nib0.pixdim(1:3);
end

% Calculate and save brain mask
img_bm = mrAnatHistogramClip(img_b0,0.4,0.99);
img_bm = dtiCleanImageMask(img_bm > 0.25*max(img_bm(:)));
%img_bm = img_bm > 0.1;
%img_bm = dtiCleanImageMask(img_bm, 9);
dtiWriteNiftiWrapper(uint8(img_bm),xformToAcPc,fullfile('bin','brainMask.nii.gz'));

img_ten = zeros([size(img_bm) 6]);
img_fa = zeros([size(img_bm)]);
img_md = zeros([size(img_bm)]);

% Run bootstrap DTI fit
if bGenBoot

    % Store brain mask in a flat file for bootstrap fitting
    fid = fopen('brainMask.img','wb');
    fwrite(fid,img_bm,'short');fclose(fid);
    if ispc
        error('Can not call modelfitboot from matlab on Windows OS.')
    end
    delete(fullfile('raw','dtboot.Bdouble'));
    cmd = sprintf('modelfitboot -inputfile %s -inversion 1 -schemefile %s -outputfile raw/dtboot.Bdouble -bootstrap %g -bgmask brainMask.img -repeats %g -components %g', rawAlignedVoxelFormatFile, schemeFile, numBootstrap, numRepeats, numDWGrads);
    disp(cmd);
    [s, ret_info] = system(cmd,'-echo');

    % Remove flat file version of the brain mask
    delete('brainMask.img');

    % Load in the bootstrap results
    % XXX get format of data file, load data beforehand.
    fid = fopen(fullfile('raw','dtboot.Bdouble'),'rb','b');
    d = fread(fid,'double'); fclose(fid);

    % Write tensor file
    img_d = shiftdim(reshape(d,[21 81 106 76]),1);
    img_ten = zeros(size(img_d(:,:,:,3:8)));
    img_ten(:,:,:,1) = img_d(:,:,:,3);
    img_ten(:,:,:,2) = img_d(:,:,:,6);
    img_ten(:,:,:,3) = img_d(:,:,:,8);
    img_ten(:,:,:,4) = img_d(:,:,:,4);
    img_ten(:,:,:,5) = img_d(:,:,:,5);
    img_ten(:,:,:,6) = img_d(:,:,:,7);
    % Convert units from m^2/s to um^2 / ms
    img_ten(:) = img_ten(:)*1e9;
    dtiWriteNiftiWrapper(img_ten,xformToAcPc,fullfile('bin','tensors.nii.gz'));

    % Write fa file
    [eigVec, eigVal] = dtiSplitTensor(img_ten);
    [img_fa,img_md] = dtiComputeFA(eigVal);
    % Apply the brain mask
    img_fa(isnan(img_fa)) = 0;
    dtiWriteNiftiWrapper(img_fa,xformToAcPc,fullfile('bin','fa.nii.gz'));

    img_pdf = img_d(:,:,:,9:end-1);
    img_pdf(:,:,:,end+1) = eigVal(:,:,:,2);
    img_pdf(:,:,:,end+1) = eigVal(:,:,:,3);
    clear img_d;

    % Write pdf file for scanner fit uncertainty estimates
    % Make pdf 0 where it is outside of brain mask, this is necessary for
    % ConTrac program
    img_pdf( repmat(img_bm, [1 1 1 size(img_pdf,4)]) == 0 ) = 0;
    dtiWriteNiftiWrapper(img_pdf,xformToAcPc,fullfile('bin','pdf.nii.gz'));
    clear img_pdf;
    delete(fullfile('raw','dtboot.Bdouble'));

else
    disp('Using previously generated pdf.nii.gz');
    ni = niftiRead(fullfile('bin','tensors.nii.gz'));
    img_ten = ni.data;
    [eigVec, eigVal] = dtiSplitTensor(img_ten);
    [img_fa,img_md] = dtiComputeFA(eigVal);
    % Apply the brain mask
    img_fa(isnan(img_fa)) = 0;
    clear ni;
end

% Generate and write white matter mask file
bUseT1 = 0;
if bUseT1
    ni = niftiRead(fullfile('bin','t1.nii.gz'));
    disp('Please wait while calculating white matter from T1 image, this can take a few minutes ...');
    [wm, gm, csf] = mrAnatSpmSegment(ni.data,ni.qto_xyz,'MNIT1');
    bb = [-size(wm)/2; size(wm)/2-1];
    % now convert to mm space
    bb = ni.qto_xyz*[bb,[0;0]]';
    bb = bb(1:3,:)';
    img_wm = mrAnatResliceSpm(wm/max(wm(:)), inv(ni.qto_xyz), bb, [2 2 2]);
    dtiWriteNiftiWrapper(uint8(img_wm>0.5),ni.qto_xyz,fullfile('bin','wmMask.nii.gz'));
else
    %img_b0clip = mrAnatHistogramClip(img_b0,0.4,0.99);
    %img_brain = dtiCleanImageMask(img_b0clip > 0.4*max(img_b0clip(:)));
    % Do FA again after the clean image mask to make sure the ventricle holes
    % were not closed
    img_wm = img_bm & img_fa>.15 & (img_md<1.1 | img_fa>0.4);
    img_wm = dtiCleanImageMask(img_wm,0,0);
    dtiWriteNiftiWrapper(uint8(img_wm),xformToAcPc,fullfile('bin','wmMask.nii.gz'));
end

% Save out dt6File
l = license('inuse'); h = hostid;
created = ['Created at ' datestr(now,31) ' by ' l(1).user ' on ' h{1} '.'];
notes = '';
b0 = img_b0;
dtBrainMask = img_bm;
dt6 = img_ten;
disp('Writing dt6.mat');
save('dt6.mat', 'dt6', 'mmPerVox', 'notes', 'xformToAcPc', 'b0', 'dtBrainMask', 'created');

% Return to starting directory
cd(oldDirName);

function val = queryOverwrite(fileName)
val = 1;
if ~isempty(dir(fileName))
    [pathstr, strippedFileName, ext, versn] = fileparts(fileName);
    msg = sprintf('Do you want to overwrite %s ? Y/N [N]: ', [strippedFileName ext]);
    reply = input(msg, 's');
    if isempty(reply) || reply == 'n' || reply == 'N'
        val = 0;
    end
end
return;

author: vistalab

repo: vistasoft


text:
function t = sensitivityinfo(fid,seekset,channum)
% T = SENSITIVITYINFO(FID,SEEKSET,CHANNUM);
%   FID = valid filepointer to a sqdfile/ sqd-filename
%   SEEKSET = starting point of file read
%   CHANNUM = Channelnumber(s)
% Gets the info regarding senstivity of sensors from the file pointed
% to by fid for given channel_number and returns a structure

if nargin<1
    error('First argument must be valid file-pointer to a sqd-file');
elseif nargin<2
    seekset = -1;
end;
filename = '';
if ischar(fid)
    filename = fid;
    fid = fopen(filename,'rb','l');
end;

if fid~=-1
    % Get channelcount
    fseek( fid, 16, seekset);
    basic_offset = fread(fid,1,'long')+(3*32+8*128+8*128)/8;
    fseek( fid, basic_offset, seekset );
    channelcount	= fread(fid,1,'int');
    if nargin<3
        channum = 0:channelcount-1;
    end;
    if (~isempty(find(channum<0)))|(~isempty(find(channum>channelcount-1)))
        error([ 'Sorry - Invalid channel number, ' num2str(channum) ', given to sensorinfo']);
    end;
    % Get offset of sensitivity values
    fseek( fid, 80, seekset );
    sensitivity_offset = fread(fid,1,'long');
    % Read sensitivity data
    fseek( fid, sensitivity_offset, seekset ); %sizeof(double) = 8 bytes and read 2 doubles
    data = fread(fid,channelcount*2,'double');
    t.Offsets  = [data(channum*2+1)]';
%     t.Gains    = [data(channum*2+2)*10^12]';
    t.Gains    = [data(channum*2+2)*10^15]';
    if ~isempty(filename)
        fclose(fid);
    end;
else
    error('Invalid filename');
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%BCSDRAWDOWN BlueChipStock drawdown analysis

addpath ./source

load BlueChipBacktest

RefIndex = 45;

[NumMonths, NumAssets] = size(RetHistory);
NumStocks = NumAssets - 3;

NumPortfolios = 40;
Periodicity = 12;

NumPeriods = floor(NumMonths/Periodicity);
StartIndex = NumMonths - NumPeriods * Periodicity;
if StartIndex < 1
	NumPeriods = NumPeriods - 1;
	StartIndex = StartIndex + Periodicity;
end
EndIndex = NumMonths;

iend = StartIndex;

PortRet = ones(NumPortfolios, NumPeriods);
IndexRet = ones(1,NumPeriods);

for k = 1:NumPeriods
	istart = iend;
	iend = istart + Periodicity;

	% calculate asset returns at specified periodicity

	A = ones(NumAssets,1);
	for i = (istart+1):iend
		for j = 1:NumAssets;
			A(j) = A(j) * (1.0 + RetHistory(i,j));
		end
	end

	for j = 1:NumAssets
		if isnan(A(j))
			A(j) = 0.0;
		end
	end

	% calculate portfolio returns at specified periodicity

	H = PortHistory{istart};
	P = H * A(1:NumStocks);

	if (k > 1)
		for i = 1:NumPortfolios
			PortRet(i,k) = PortRet(i,k - 1) * P(i);
		end
		IndexRet(k) = IndexRet(k - 1) * A(RefIndex);
	end
end

MaxDD = maxdrawdown(PortRet');
IDD = maxdrawdown(IndexRet');
IDD = repmat(IDD,1,NumPortfolios);

figure(1);
plot(MaxDD);
hold all
plot(IDD)
set(gca,'ylim',[-1 0]);
ylabel('\bfMaximum Drawdown');
xlabel('\bfPortfolio Number on Efficient Frontier');
title('\bfDrawdown of Portfolio Sequence');
hold off


text:
function fem1d_heat_explicit_test01 ( )

%*****************************************************************************80
%
%% FEM1D_HEAT_EXPLICIT_TEST01 runs a simple test.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_HEAT_EXPLICIT_TEST01:\n' );
  fprintf ( 1, '  The time dependent 1D heat equation is\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Ut - k * Uxx = f(x,t)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  for space interval A <= X <= B with boundary conditions\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    U(A,t) = UA(t)\n' );
  fprintf ( 1, '    U(B,t) = UB(t)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  and time interval T0 <= T <= T1 with initial condition\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    U(X,T0) = U0(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To compute an approximate solution:\n' );
  fprintf ( 1, '    the interval [A,B] is replace by a discretized mesh Xi;\n' );
  fprintf ( 1, '    a set of finite element functions PSI(X) are determined,\n' );
  fprintf ( 1, '    the solution U is written as a weighted sum of the basis functions,\n' );
  fprintf ( 1, '    the weak form of the differential equation is written,\n' );
  fprintf ( 1, '    a time grid Tj is imposed, and time derivatives replaced by\n' );
  fprintf ( 1, '    an explicit forward Euler first difference,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    The continuous PDE has now been transformed into a set of algebraic\n' );
  fprintf ( 1, '    equations for the coefficients C(Xi,Tj).\n' );
%
%  Set the nodes.
%
  x_num = 21;
  x_min = 0.0;
  x_max = 1.0;
  dx = ( x_max - x_min ) / ( x_num - 1 );
  x = ( linspace ( x_min, x_max, x_num ) )';
%
%  Set the times.
%
  t_num = 401;
  t_min = 0.0;
  t_max = 80.0;
  dt = ( t_max - t_min ) / ( t_num - 1 );
  t = ( linspace ( t_min, t_max, t_num ) )';
%
%  Set finite element information.
%
  element_num = x_num - 1;
  element_node(1,1:element_num) = 1 : x_num - 1;
  element_node(2,1:element_num) = 2 : x_num;
  quad_num = 3;
  mass = assemble_mass ( x_num, x, element_num, element_node, quad_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of X nodes = %d\n', x_num );
  fprintf ( 1, '  X interval = [ %f, %f ]\n', x_min, x_max );
  fprintf ( 1, '  X step size = %f\n', dx );
  fprintf ( 1, '  Number of T steps = %d\n', t_num );
  fprintf ( 1, '  T interval = [ %f, %f ]\n', t_min, t_max );
  fprintf ( 1, '  T step size = %f\n', dt );
  fprintf ( 1, '  Number of elements = %d\n', element_num );
  fprintf ( 1, '  Number of quadrature points = %d\n', quad_num );

  u_mat = zeros ( x_num, t_num );

  for j = 1 : t_num

    if ( j == 1 )

      u(1:x_num,1) = ic_test01 ( x_num, x, t(j) );
      u(1:x_num,1) = bc_test01 ( x_num, x, t(j), u(1:x_num,1) );

    else

      u(1:x_num,1) = fem1d_heat_explicit ( x_num, x, t(j-1), dt, @k_test01, ...
        @rhs_test01, @bc_test01, element_num, element_node, quad_num, mass, u );

    end

    u_mat(1:x_num,j) = u(1:x_num,1);

  end
%
%  Make a product grid of T and X for plotting.
%
  [ t_mat, x_mat ] = meshgrid ( t, x );
%
%  Make a mesh plot of the solution.
%
  mesh ( t_mat, x_mat, u_mat )
  xlabel ( '<--Time-->' )
  ylabel ( '<--X-->' )
  zlabel ( '<--H(X,T)-->' )
  title ( 'H(X,T) for TEST01, by FEM1D\_HEAT\_EXPLICIT' )
%
%  Write the data to files.
%
  r8mat_write ( 'h_test01.txt', x_num, t_num, u_mat );
  r8vec_write ( 't_test01.txt', t_num, t );
  r8vec_write ( 'x_test01.txt', x_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  H(X,T) written to "h_test01.txt"\n' );
  fprintf ( 1, '  T values written to "t_test01.txt"\n' );
  fprintf ( 1, '  X values written to "x_test01.txt"\n' );

  return
end

author: johannesgerer


text:
 function sls = de_ftab_sls(varargin)
%function sls = de_ftab_sls(varargin)
%|
%| Determine structure that characterizes the "s" limits for polyenergetic CT,
%| where s_l is a line integral through the lth material type, l=1...L.
%| L is the number of material components.
%|
%| in
%| option
%|	'sl'	cell{LL} sample thicknesses for each of LL materials
%|	'n'	[L,1]	number of samples of lth material integrals
%|			default: [50 30] or [length(sl{ll})]
%|	'max'	[L,1]	maximum material density line integrals, units: g/cm^2
%|			default: [45 43] or [max(sl{ll})]
%|	'min'	[L,1]	minimum material density line integrals, units: g/cm^2
%|			default: [0] or [min(sl{ll})]
%|
%| out
%|	sls	strum
%|		data:
%|		sls.sl,n,max,min	see above
%|		methods:
%|		.sll			[n1,n2,...,nL,L] ndgrid of sl{*} values
%|
%| Copyright 2008-6-15, Jeff Fessler, University of Michigan

if nargin == 1 && streq(varargin{1}, 'test')
	de_ftab_sls_test
return
end
if ~nargout, ir_usage, end

sls.sl = {};
sls.n = [];
sls.min = [];
sls.max = [];

sls = vararg_pair(sls, varargin);

[sls.sl sls.n sls.min sls.max] = ...
	de_ftab_sls_do(sls.sl, sls.n, sls.min, sls.max);

sls = strum(sls, {'sll', @de_ftab_sls_sll, '()'});

% de_ftab_sls_sll()
function out = de_ftab_sls_sll(sls, varargin)
sll = ndgrid_jf('mat', sls.sl);
out = sll(varargin{:});

% de_ftab_sls_do()
function [sl, s_n, s_min, s_max] = de_ftab_sls_do(sl, s_n, s_min, s_max)

% number of samples of the material integrals "s"
if isempty(s_n)
	if isempty(sl)
		s_n = [45 43];
%		s_n = [51 31]; % makes plot_jac look nice
	else
		for ll=1:length(sl)
			s_n(1,ll) = length(sl{ll});
		end
	end
end

% minimum material "integrals"
if isempty(s_min)
	if isempty(sl)
		s_min = [0 0];
	else
		for ll=1:length(sl)
			s_min(1,ll) = min(sl{ll});
		end
	end
end

% maximum material "integrals"
if isempty(s_max)
	if isempty(sl)
		% soft max: 50cm * 1g/cc
		% bone max: 15cm * 2g/cc (for now)
		s_max = [50 30];
	else
		for ll=1:length(sl)
			s_max(1,ll) = max(sl{ll});
		end
	end
end

if isempty(sl)
	for ll=1:length(s_n)
		sl{ll} = linspace(s_min(ll), s_max(ll), s_n(ll))';
	end
end

% de_ftab_sls_test
function de_ftab_sls_test
sls = de_ftab_sls;
pr 'size(sls.sll)'
pr 'size(sls.sll(:,:,1))'


text:
function cl=cleanRequired(seqFolder)

cl = ~isempty(strfind(seqFolder,'MOT16')) || ~isempty(strfind(seqFolder,'MOT17'));

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
% LFCalRectifyLF - rectify a light field using a calibrated camera model, called as part of LFUtilDecodeLytroFolder
%
% Usage:
%     [LF, RectOptions] = LFCalRectifyLF( LF, CalInfo, RectOptions )
%     [LF, RectOptions] = LFCalRectifyLF( LF, CalInfo )
%
% This function is called by LFUtilDecodeLytroFolder to rectify a light field. It follows the
% rectification procedure described in:
%
% D. G. Dansereau, O. Pizarro, and S. B. Williams, "Decoding, calibration and rectification for
% lenslet-based plenoptic cameras," in Computer Vision and Pattern Recognition (CVPR), IEEE
% Conference on. IEEE, Jun 2013.
%
% Minor differences from the paper: light field indices [i,j,k,l] are 1-based in this
% implementation, and not 0-based as described in the paper.
%
% Note that a calibration should only be applied to a light field decoded using the same lenslet
% grid model. This is because the lenslet grid model forms an implicit part of the calibration,
% and changing it will invalidate the calibration.
%
% LFCalDispRectIntrinsics is useful for visualizing the sampling pattern associated with a requested
% intrinsic matrix.
%
% Inputs:
%
%         LF : The light field to rectify; should be floating point
%
%     CalInfo struct contains a calibrated camera model, see LFUtilCalLensletCam:
%            .EstCamIntrinsicsH : 5x5 homogeneous matrix describing the lenslet camera intrinsics
%            .EstCamDistortionV : Estimated distortion parameters
%
%    [optional] RectOptions struct (all fields are optional) :
%        .NInverse_Distortion_Iters : Number of iterations in inverse distortion estimation
%                        .Precision : 'single' or 'double'
%               .RectCamIntrinsicsH : Requests a specific set of intrinsics for the rectified light
%                                     field. By default the rectified intrinsic matrix is
%                                     automatically constructed from the calibrated intrinsic
%                                     matrix, but this process can in some instances yield poor
%                                     results: excessive black space at the edges of the light field
%                                     sample space, or excessive loss of scene content off the edges
%                                     of the space. This parameters allows you to fine-tune the
%                                     desired rectified intrinsic matrix.
%
% Outputs :
%
%     LF : The rectified light field
%     RectOptions : The rectification options as applied, including any default values employed.
%
% User guide: LFToolbox.pdf
% See also: LFCalDispRectIntrinsics, LFUtilCalLensletCam, LFUtilDecodeLytroFolder, LFUtilProcessWhiteImages

% Copyright (c) 2013-2020 Donald G. Dansereau

function [LF, RectOptions] = LFCalRectifyLF( LF, CalInfo, RectOptions )

%---Defaults---
RectOptions = LFDefaultField( 'RectOptions', 'Precision', 'single' );
RectOptions = LFDefaultField( 'RectOptions', 'NInverse_Distortion_Iters', 2 );
RectOptions = LFDefaultField( 'RectOptions', 'MaxUBlkSize', 32 );
LFSize = size(LF);
RectOptions = LFDefaultField( 'RectOptions', 'RectCamIntrinsicsH', LFDefaultIntrinsics( LFSize, CalInfo ) );

%---Build interpolation indices---
fprintf('Generating interpolation indices...\n');
NChans = LFSize(5);
LF = cast(LF, RectOptions.Precision);

%---chop up the LF along u---
fprintf('Interpolating...');
UBlkSize = RectOptions.MaxUBlkSize;
LFOut = LF;
for( UStart = 1:UBlkSize:LFSize(4) )
    UStop = UStart + UBlkSize - 1;
    UStop = min(UStop, LFSize(4));

    t_in=cast(1:LFSize(1), 'uint16'); % saving some mem by using uint16
    s_in=cast(1:LFSize(2), 'uint16');
    v_in=cast(1:LFSize(3), 'uint16');
    u_in=cast(UStart:UStop, 'uint16');
    [tt,ss,vv,uu] = ndgrid(t_in,s_in,v_in,u_in);

    % InterpIdx initially holds the index of the desired ray, and is evolved through the application
    % of the inverse distortion model to eventually hold the continuous-domain index of the undistorted
    % ray, and passed to the interpolation step.
    InterpIdx = [ss(:)'; tt(:)'; uu(:)'; vv(:)'; ones(size(ss(:)'))];
    DestSize = size(tt);
    clear tt ss vv uu

    InterpIdx = LFMapRectifiedToMeasured( InterpIdx, CalInfo, RectOptions );

    for( ColChan = 1:NChans )
        % todo[optimization]: use a weighted interpolation scheme to exploit the weight channel
        InterpSlice = interpn(squeeze(LF(:,:,:,:,ColChan)), InterpIdx(2,:),InterpIdx(1,:), InterpIdx(4,:),InterpIdx(3,:), 'linear');
        InterpSlice = reshape(InterpSlice, DestSize);
        LFOut(:,:,:,UStart:UStop,ColChan) = InterpSlice;
    end

    fprintf('.')
end
LF = LFOut;
clear LFOut;

%---Clip interpolation result, which sometimes rings slightly out of range---
LF(isnan(LF)) = 0;
LF = max(0, min(1, LF));

fprintf('\nDone\n');

end

repo: LFToolbox


text:
clc; clear;

% Name of the video
videoString = "01";

% Define path of background image. Can be empty.
backgroundImagePath = sprintf('data/%s_highway.jpg', videoString);
% backgroundImagePath = "";

% Read tracks by using tracks file and static tracks file.
tracksFilename = sprintf('data/%s_tracks.csv', videoString);
tracksStaticFilename = sprintf('data/%s_tracksMeta.csv', videoString);
tracks = readInTracksCsv(tracksFilename, tracksStaticFilename);
% Read video meta data from video meta file
videoMetaFilename = sprintf('data/%s_recordingMeta.csv', videoString);
videoMeta = readInVideoCsv(videoMetaFilename);

% The visualization needs the tracks, the videoMeta data and the
% backgroundImagePath (which can be empty).
argumentsInput.tracks = tracks;
argumentsInput.videoMeta = videoMeta;
argumentsInput.backgroundImagePath = backgroundImagePath;

% Plot variables initialization. These are the default values. Please use
% them or modify them for your own purpose.
argumentsInput.removingList = [];
argumentsInput.trackWidth = 408; % Width of the video frame
argumentsInput.laneThickness = 0.3; % thickness of the lane markings
argumentsInput.laneColor = [1 1 1]; % color of the lane markings (white)
argumentsInput.currentFrame = 1; % Starting frame
argumentsInput.minimumFrame = 1; % Minimum frame that exists (1)
argumentsInput.maximumFrame = tracks(end).finalFrame-1; % Maximum frame
argumentsInput.playRate = 1; % The play rate (going "playRate" steps forward)
argumentsInput.stop = false; % Boolean for
argumentsInput.streetColor = [0.6 0.6 0.6]; % The street color
argumentsInput.boundingBoxColor = [1 0 0]; % The color of the bounding box
argumentsInput.plotTextAnnotation = true; % Whether the text annotations shall be plotted
argumentsInput.plotTrackingLine = true; % Whether the tracking line shall be plotted

argumentsInput.geoParams.meterPerPixel = 0.10106;
argumentsInput.geoParams.kmhPerPixelFrame = 9.0954;

% Start the visualization
trackVisualization(argumentsInput);

author: RobertKrajewski

repo: highD-dataset


text:
function varargout = ft_info(varargin)

% FT_INFO prints an info message on screen, depending on the verbosity
% settings of the calling high-level FieldTrip function.
%
% Use as
%   ft_info(...)
% with arguments similar to fprintf, or
%   ft_info(msgId, ...)
% with arguments similar to warning.
%
% You can switch of all messages using
%   ft_info off
% or for specific ones using
%   ft_info off msgId
%
% To switch them back on, you would use
%   ft_info on
% or for specific ones using
%   ft_info on msgId
%
% Messages are only printed once per timeout period using
%   ft_info timeout 60
%   ft_info once
% or for specific ones using
%   ft_info once msgId
%
% You can see the most recent messages and identifier using
%   ft_info last
%
% You can query the current on/off/once state for all messages using
%   ft_info query
%
% See also FT_ERROR, FT_WARNING, FT_NOTICE, FT_INFO, FT_DEBUG, ERROR, WARNING

% Copyright (C) 2017, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargout
  varargout{:} = ft_notification(varargin{:});
elseif isequal(varargin, {'last'})
  % return an answer anyway
  varargout{1} = ft_notification(varargin{:});
else
  ft_notification(varargin{:});
end

author: spm