text:
% LINESEG - Form straight line segements from an edge list.
%
% Usage: [seglist, nedgelist] = lineseg(edgelist, tol, angtol, linkrad)
%
% Arguments:  edgelist - Cell array of edgelists (row col) coords.
%             tol      - Maximum deviation from straight line before a
%                        segment is broken in two (measured in pixels).
%             angtol   - Angle tolerance used when attempting to merge line
%                        segements (radians).
%             linkrad  - Maximum distance between end points of line
%                        segments for segments to be elegible for
%                        linking (pixels).
%  angtol and linkrad are optional.  If these parameters are omitted the
%  merging phase is omitted.
%
% Returns:  seglist - an Nx4 array storing line segments in the form
%                      [x1 y1 x2 y2
%                       x1 y1 x2 y2
%                          . . .    ] etc
%
%         nedgelist - a new cell array of edge lists where each edge list
%                     corresponds to each segment in seglist above.  The
%                     edgelist is in row,column coords in the form
%                     { [r1 c1   [r1 c1   etc }
%                        r2 c2    ...
%                        ...
%                        rN cN]   ....]
%
% * Note that a non-empty nedgelist is only returned if there is no segment
%   merging.
%
% This function takes each array of edgepoints in edgelist, finds the
% size and position of the maximum deviation from the line that joins the
% endpoints, if the maximum deviation exceeds the allowable tolerance the
% edge is shortened to the point of maximum deviation and the test is
% repeated.  In this manner each edge is broken down to line segments,
% each of which adhere to the original data with the specified tolerance.
%
% The optional final edge merging phase is provided because the initial
% edge linking phase may have separated edges at `Y' junctions in the
% image.  The merging phase can reconnect broken branches.
% Note however that the merging process can be slow.
%
% See also:  EDGELINK, MAXLINEDEV, MERGESEG, DRAWSEG
%

% Copyright (c) 2000-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% 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, 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.

% December 2000 - Original version
% February 2003 - Added the returning of nedgelist data.

% ** need to preallocate memory to improve speed ***

function [linelist, nedgelist] = lineseg(edgelist, tol, angtol, linkrad)

    linelist = [];
    Nline = 0;
    Nedge = length(edgelist);

    if nargin == 2
	merge = 0;
    else
	merge = 1;
    end

    for e = 1:Nedge
        y = edgelist{e}(:,1);
	x = edgelist{e}(:,2);

	fst = 1;                % Indecies of first and last points in edge
	lst = length(x);        % segment being considered.

	while  fst<lst
	    [m,i,d] = maxlinedev(x(fst:lst),y(fst:lst));  % Find size & posn of
                                                          % maximum deviation.

	    while m > tol       % While deviation is > tol  (m/d) ?
		lst = i+fst-1;  % Shorten line to point of max deviation by adjusting lst
		[m,i,d] = maxlinedev(x(fst:lst),y(fst:lst));
	    end

	    Nline = Nline+1;
            % Record line segment. Note that (c,r) corresponds to (x,y)
	    linelist(Nline,:) = [x(fst) y(fst) x(lst) y(lst)];
	    % Record edgelist that corresponds to the segment.
	    nedgelist{Nline} = [y(fst:lst) x(fst:lst)];

%	    fst = lst+1;        % Reset fst and lst.
	    fst = lst;        % make new fst match last lst so that segments
                              % share endpoints.
	    lst = length(x);
	end

    end
    % fprintf('No of segments = %d\n',length(linelist));
    % drawseg(linelist,1), title('Raw segments');

    if merge
%	fprintf('\nMerging Segments\n');
	linelist = mergeseg(linelist, angtol, linkrad, tol);
	nedgelist = {};
%	fprintf('No of merged segments = %d\n',length(linelist));
        % drawseg(linelist,2), title('Merged segments');
    end

author: zouchuhang

repo: LayoutNet


text:
function [a]=plus(b,c)
%A=B+C
%   [A]=PLUS(B,C) A is the sum of TT-matrices B and C
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
if ( isempty(b) )
   a=c;
   return
elseif (isempty(c) )
   a=b;
   return
end
    a=tt_matrix;
a.n=b.n;
a.m=b.m;
a.tt=b.tt+c.tt;
return
end

author: oseledets

repo: TT-Toolbox


text:
function pp = pchipd(x,y,d,xx)
%PCHIPD  Piecewise Cubic Hermite Interpolating Polynomial with Derivatives.
%   PP = PCHIPD(X,Y,D) provides the piecewise cubic polynomial which
%   interpolates values Y and derivatives D at the sites X.  This is meant
%   to augment the built-in Matlab function PCHIP, which does not allow the
%   user to specify derivatives.
%
%   X must be a vector.
%
%   If Y and D are vectors, then Y(i) and D(i) are the value and derivative
%   to be matched at X(i).
%
%   If Y and D are matrices, then size(Y,2) == size(D,2) == length(X).
%   Also, size(Y,1) == size(D,1).  Use this for interpolating vector valued
%   functions.
%
%   YY = PCHIPD(X,Y,D,XX) is the same as YY = PPVAL(PCHIPD(X,Y,D),XX), thus
%   providing, in YY, the values of the interpolant at XX.
%
%   Example comparing SPLINE, PCHIP, and PCHIPD
%     a = -10;
%     b = 10;
%     x = linspace(a,b,7);
%     f = @(x) 1./(1+exp(-x));  % logistic function
%     df = @(x) f(x).*(1-f(x)); % derivative of the logistic function
%     t = linspace(a,b,50);
%     r = f(t);
%     p = pchip(x,f(x),t);
%     s = spline(x,f(x),t);
%     q = pchipd(x,f(x),df(x),t);
%     plot(t,r,'k',x,f(x),'o',t,p,'-',t,s,'-.',t,q,'--')
%     legend('true','data','pchip','spline','pchipd',4)
%
%   See also INTERP1, SPLINE, PCHIP, PPVAL, MKPP, UNMKPP.
%

%
% 2010-10-04 (nwh) first version
%

% check inputs
% x must be a vector
if ~isvector(x)
  error('pchipd:input_error','x must be a vector of length > 2.')
end

% get size and orient
n = length(x);
x = x(:);

% make sure x is long enough, we can't construct an interpolating
% polynomial with just one point
if n < 2
  error('pchipd:input_error','x must be a vector of length > 2.')
end

% check y and d
if isvector(y) && isvector(d) && length(y) == n && length(d) == n
  % orient
  y = y(:);
  d = d(:);
  m = 1;
elseif size(y,2) == n && size(d,2) == n && size(y,1) == size(d,1)
  m = size(y,1);
  y = y';
  d = d';
else
  error('pchipd:input_error','y and d must be vectors or matrices of same size with length(x) columns.')
end

% sort breaks & data if needed
if ~issorted(x)
  [x x_ix] = sort(x);
  if m == 1
    y = y(x_ix);
    d = d(x_ix);
  else
    y = y(x_ix,:);
    d = d(x_ix,:);
  end
end

% compute coefficients
coef = zeros(m,n-1,4);
dx = diff(x);
for i = 1:m
  dy = diff(y(:,i));
  coef(i,:,4) = y(1:end-1,i)';
  coef(i,:,3) = d(1:end-1,i)';
  coef(i,:,2) = 3*dy./(dx.^2) - (2*d(1:end-1,i)+d(2:end,i))./dx;
  coef(i,:,1) = -2*dy./(dx.^3) + (d(1:end-1,i)+d(2:end,i))./(dx.^2);
end

% create the piecewise polynomial structure
pp = mkpp(x,coef,m);

% if user requests evaluations
if nargin > 3
  pp = ppval(pp,xx);
end

author: Sable

repo: mcbench-benchmarks


text:
%Disk logger test
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updaetd 2/3/2011
classdef DiskLogger_Test < TestCase
    properties
        simulation
    end

    methods
        function this = DiskLogger_Test(name)
            this = this@TestCase(name);
        end

        function setUp(this)
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], true);
            sim.applyOptions('lengthSec', 10);
            this.simulation = sim;
        end
    end

    methods
        function testLogging(this)
            import edu.stanford.covert.cell.sim.util.DiskLogger;

            sim = this.simulation;
            pm = sim.state('ProteinMonomer');

            metadata = edu.stanford.covert.cell.sim.util.ConditionSet.parseConditionSet(...
                sim, 'src_test/+edu/+stanford/+covert/+cell/+sim/+util/fixtures/runSimulations.xml').metadata;
            metadata.knowledgeBaseWID = 145;
            metadata.revision = edu.stanford.covert.util.revision;

            outputDirectory = 'output/runSmallTests/DiskLogger';
            delete([outputDirectory '/*']);
            segmentSizeSec = 2;
            logger = DiskLogger(outputDirectory, segmentSizeSec, metadata);
            logger.addMetadata(metadata);
            sim.run(logger);

            DiskLogger.reindexTimeCourses(outputDirectory);

            metadatafields = {'shortDescription'; 'longDescription';
                'email'; 'lastName'; 'firstName'; 'affiliation'; 'knowledgeBaseWID';
                'startTime'; 'endTime'; 'lengthSec'; 'revision'; 'differencesFromRevision';
                'userName'; 'hostName'; 'ipAddress';
                'outputDirectory'; 'segmentSizeStep'; 'downsampleStepSec';
                'stateNames'; 'dependentStateNames'};
            assertEqual(sort(metadatafields), sort(fieldnames(logger.metadata)));

            %% assert created files
            files = dir(outputDirectory);
            fileNames = {'.'; '..'; 'fittedConstants.mat'; 'metadata.mat'; 'options.mat'; 'parameters.mat';
                'state-0.mat'; 'state-1.mat'; 'state-2.mat'; 'state-3.mat'; 'state-4.mat'; 'state-5.mat';
                'randStreamStates.mat'};
            assertTrue(all(ismember(fileNames, {files.name}')));

            %% assert loading back from disk
            downsampleStepSec = 2;
            nTimepoints = sim.lengthSec / downsampleStepSec;

            %extraction downsampling
            [states, metadata, options, parameters, fittedConstants, randStreamStates] = ...
                DiskLogger.load(outputDirectory, '-all', [], [], downsampleStepSec, 'extract');
            metadata.downsampleStepSec = logger.metadata.downsampleStepSec;
            metadata.segmentSizeStep = logger.metadata.segmentSizeStep;

            assertEqual(sort(cellfun(@(x) x(7:end), sim.stateWholeCellModelIDs, 'UniformOutput', false)), sort(fieldnames(states)));
            assertEqual([1 1 nTimepoints], size(states.Time.values));
            assertEqual(permute([1:downsampleStepSec:sim.lengthSec-3 sim.lengthSec], [1 3 2]), states.Time.values);
            assertEqual([numel(sim.getRandStreamStates().simulation) sim.lengthSec+1], size(randStreamStates.simulation));

            assertEqual(logger.metadata, metadata);
            assertEqual(sim.getOptions, options);
            assertEqual(sim.getParameters, parameters);
            assertEqual(sim.getFittedConstants, fittedConstants);

            pmCounts = states.ProteinMonomer.counts;

            states = DiskLogger.load(outputDirectory, {'ProteinMonomer' 'counts' '-sum' ':'}, [], [], downsampleStepSec, 'extract');
            assertEqual(sum(pmCounts, 1), states.ProteinMonomer.counts);

            states = DiskLogger.load(outputDirectory, {'ProteinMonomer' 'counts' '-sum' '-sum'}, [], [], downsampleStepSec, 'extract');
            assertEqual(sum(sum(pmCounts, 1), 2), states.ProteinMonomer.counts);

            states = DiskLogger.load(outputDirectory, {'ProteinMonomer' 'counts' ':' [3 4]}, [], [], downsampleStepSec, 'extract');
            assertEqual(sum(pmCounts(:, [3 4], :), 2), states.ProteinMonomer.counts);

            states = DiskLogger.load(outputDirectory, {'ProteinMonomer' 'counts' [pm.matureIndexs pm.boundIndexs] [1 4]}, [], [], downsampleStepSec, 'extract');
            assertEqual(sum(pmCounts(pm.matureIndexs, [1 4], :) + pmCounts(pm.boundIndexs, [1 4], :), 2), states.ProteinMonomer.counts);

            %mean downsampling
            [states, metadata, options, parameters, fittedConstants] = ...
                DiskLogger.load(outputDirectory, '-all', [], [], downsampleStepSec, 'mean');
            metadata.downsampleStepSec = logger.metadata.downsampleStepSec;

            assertEqual(sort(cellfun(@(x) x(7:end), sim.stateWholeCellModelIDs, 'UniformOutput', false)), sort(fieldnames(states)));
            assertEqual([1 1 nTimepoints], size(states.Time.values));
            assertEqual(permute(0.5+(1:downsampleStepSec:sim.lengthSec), [1 3 2]), states.Time.values);

            assertEqual(logger.metadata, metadata);
            assertEqual(sim.getOptions, options);
            assertEqual(sim.getParameters, parameters);
            assertEqual(sim.getFittedConstants, fittedConstants);

            states = DiskLogger.load(outputDirectory, '-independent', [], [], downsampleStepSec, 'extract'); %#ok<NASGU>
            states = DiskLogger.load(outputDirectory, '-dependent', [], [], downsampleStepSec, 'extract'); %#ok<NASGU>

            %test loadSimulation
            [sim2, states, ~, ~, ~, metadata] = loadSimulation([pwd filesep outputDirectory]);
            metadata.downsampleStepSec = logger.metadata.downsampleStepSec;

            assertEqual(options, sim2.getOptions());
            assertEqual(parameters, sim2.getParameters());
            assertEqual(fittedConstants, sim2.getFittedConstants());
            assertEqual(metadata, logger.metadata);

            states2 = struct;
            for i = 1:numel(sim2.states)
                state = sim2.states{i};
                stateID = state.wholeCellModelID(7:end);
                states2.(stateID) = struct;
                for j = 1:numel(state.stateNames)
                    name = state.stateNames{j};
                    states2.(stateID).(name) = state.(name);
                end
            end
            assertEqual(states, states2);
        end

        function testEarlyTermination(this)
            import edu.stanford.covert.cell.sim.util.DiskLogger;

            %setup
            sim = this.simulation;

            metadata = edu.stanford.covert.cell.sim.util.ConditionSet.parseConditionSet(sim, 'src_test/+edu/+stanford/+covert/+cell/+sim/+util/fixtures/runSimulations.xml').metadata;
            metadata.knowledgeBaseWID = 145;

            outputDirectory = 'output/runSmallTests/DiskLogger';
            delete([outputDirectory '/*']);
            segmentSizeSec = 2;
            logger = DiskLogger(outputDirectory, segmentSizeSec, metadata);
            logger.addMetadata(metadata);

            %run simulation
            logger.initialize(sim);

            for i = 1:sim.getNumSteps / 2
                sim.evolveState();
                logger.append(sim);
            end

            logger.finalize(sim);

            %test
            files = dir(outputDirectory);
            fileNames = {'.'; '..'; 'fittedConstants.mat'; 'metadata.mat'; 'options.mat'; 'parameters.mat';
                'state-0.mat'; 'state-1.mat'; 'state-2.mat'; 'state-3.mat'};
            assertTrue(all(ismember(fileNames, {files.name}')));

            [states, metadata] = logger.load(outputDirectory, {'Time' 'values'}, [], [], 2, 'extract');
            assertEqual(sim.getNumSteps / 2, metadata.lengthSec);
            assertEqual(permute([1 3 5], [1 3 2]), states.Time.values);

            [states, metadata] = logger.load(outputDirectory, {'Time' 'values'}, [], [], 2, 'mean');
            assertEqual(sim.getNumSteps / 2, metadata.lengthSec);
            assertEqual(permute([1.5 3.5 5], [1 3 2]), states.Time.values);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%% platonic_solid
% Below is a demonstration of the features of the |platonic_solid| function

%%
clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceColor='b';
faceAlpha=1;
edgeColor='k';
edgeWidth=2;
markerSize=5;

%%

hf=cFigure; % Open figure for plotting
pColor=gjet(5);
for q=1:1:5
    %Defining the faces (F) and vertices (V) of a platonic solid
    [V,F]=platonic_solid(q,1); %q indicates solid type, r is the radius

    subplot(2,3,q); hold on;

    hp=gpatch(F,V,pColor(q,:),'k',faceAlpha,3);
    patchNormPlot(F,V);

    axisGeom(gca,fontSize);
    camlight('headlight'); lighting flat;
    axis off;
end
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% P3.11
colordef white; clear;  clc;
M = 250;    k = -M:M;   w = (pi/M)*k;   % [0, pi] axis divided into 501 points.
w0 = pi/2;  figure(1);
n = -50:50;

h = (0.9).^(abs(n));
H = dtft(h,n,w);
subplot(3,2,1); plot(w/pi,abs(H),'k');
xlabel('frequency in pi units');    ylabel('|H|');  title('Magnitude Response');
subplot(3,2,2); plot(w/pi,angle(H)/pi,'k');
xlabel('frequency in pi units');    ylabel('\Theta');   title('Phase Response');

h = sinc(0.2*n).*((stepseq(-20,-50,50)-stepseq(20,-50,50)));
H = dtft(h,n,w);
subplot(3,2,3); plot(w/pi,abs(H),'k');
xlabel('frequency in pi units');    ylabel('|H|');  title('Magnitude Response');
subplot(3,2,4); plot(w/pi,angle(H)/pi,'k');
xlabel('frequency in pi units');    ylabel('\Theta');   title('Phase Response');

h = sinc(0.2*n).*((stepseq(0,-50,50)-stepseq(40,-50,50)));
H = dtft(h,n,w);
subplot(3,2,5); plot(w/pi,abs(H),'k');
xlabel('frequency in pi units');    ylabel('|H|');  title('Magnitude Response');
subplot(3,2,6); plot(w/pi,angle(H)/pi,'k');
xlabel('frequency in pi units');    ylabel('\Theta');   title('Phase Response');

figure(2);
h = ((0.5).^n+(0.4).^n).*stepseq(0,-50,50);
H = dtft(h,n,w);
subplot(3,2,1); plot(w/pi,abs(H),'k');
xlabel('frequency in pi units');    ylabel('|H|');  title('Magnitude Response');
subplot(3,2,2); plot(w/pi,angle(H)/pi,'k');
xlabel('frequency in pi units');    ylabel('\Theta');   title('Phase Response');

h = (0.5).^(abs(n)).*cos(0.1*pi*n);
H = dtft(h,n,w);
subplot(3,2,3); plot(w/pi,abs(H),'k');
xlabel('frequency in pi units');    ylabel('|H|');  title('Magnitude Response');
subplot(3,2,4); plot(w/pi,angle(H)/pi,'k');
xlabel('frequency in pi units');    ylabel('\Theta');   title('Phase Response');


text:
function set9701a(m)
% sets parameters of DT3152 according to
% http://cmp.felk.cvut.cz/cmp/hardware/pulnix.html#TM-9701analog
% where m is a initialized handle to a framegrabber

% input parameters check
if nargin < 1,
	error('Not enough input arguments!');
end

if nargin > 1,
	error('Too many input arguments!');
end

for i=0:3
   setfg(m,'InputChannel',i)
   setfg(m,'FrameType','NonInterlaced')
   setfg(m,'VideoType','Composite')
   setfg(m,'ClockSource','internal')
   setfg(m,'ClockFreq',14318180)
   setfg(m,'TotalPixPerLine',909)

   setfg(m,'backporchstart',60)
   setfg(m,'clampstart',95)
   setfg(m,'clampend',100)

   setfg(m,'FirstActivePixel',126)
   setfg(m,'ActivePixelCount',760)

   setfg(m,'TotalLinesperFld',524)
   setfg(m,'FirstActiveLine',33)
   setfg(m,'ActiveLineCount',484)

   setfg(m,'FrameLeft',0)
   setfg(m,'HorFrameInc',1)
   setfg(m,'FrameWidth',760)
   setfg(m,'FrameTop',0)
   setfg(m,'VerFrameInc',1)
   setfg(m,'FrameHeight',484)
end

setfg(m,'InputChannel',0)


text:
function res = transpose(obj);

    res = ctranspose(obj, 0);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% Compose two vector fields
%  Changed: Dec 31st, 2011
%
function [vx,vy,vz] = compose(ax,ay,az, bx,by,bz)

    % Piggyback
    piggybackscale = 1.2; % just to get too many outside values
    [ax,lim] = piggyback(ax,piggybackscale);
    [ay,lim] = piggyback(ay,piggybackscale);
    [az,lim] = piggyback(az,piggybackscale);
    [bx,lim] = piggyback(bx,piggybackscale);
    [by,lim] = piggyback(by,piggybackscale);
    [bz,lim] = piggyback(bz,piggybackscale);

    % Coordinates
    nx = size(ax,1);
    ny = size(ax,2);
    nz = size(ax,3);

    [y,x,z] = ndgrid(1:nx, 1:ny, 1:nz); % coordinate image

    % Where points are going
    xp  = iminterpolate(x+ax, bx,by,bz);
    yp  = iminterpolate(y+ay, bx,by,bz);
    zp  = iminterpolate(z+az, bx,by,bz);

    % Update field
    vx = xp - x;
    vy = yp - y;
    vz = zp - z;

    % Zero vectors going outside the image
    zr  = (xp==0 & yp==0 & zp==0);
    vx(zr) = 0;
    vy(zr) = 0;
    vz(zr) = 0;

    % Unpiggyback
    vx = vx(lim(1):lim(2),lim(3):lim(4),lim(5):lim(6));
    vy = vy(lim(1):lim(2),lim(3):lim(4),lim(5):lim(6));
    vz = vz(lim(1):lim(2),lim(3):lim(4),lim(5):lim(6));

end


text:
function V = mttkrp(X,U,n)
%MTTKRP Matricized tensor times Khatri-Rao product for ktensor.
%
%   V = MTTKRP(X,U,n) efficiently calculates the matrix product of the
%   n-mode matricization of X with the Khatri-Rao product of all
%   entries in U, a cell array of matrices, except the nth.  How to
%   most efficiently do this computation depends on the type of tensor
%   involved.
%
%   See also KTENSOR, KTENSOR/TTV
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

N = ndims(X);

if (n==1)
    R = size(U{2},2);
else
    R = size(U{1},2);
end

% Compute matrix of weights
W = repmat(X.lambda,1,R);
for i = [1:n-1,n+1:N]
  W = W .* (X.u{i}' * U{i});
end

% Find each column of answer by multiplying columns of X.u{n} with weights
V = X.u{n} * W;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [model, partition_size] = fit_partitioned_model(...
    inputs, outputs, selectors, sel_sizes, min_size, partition_names, fn_name, varargin)
%function [models, partition_sizes] = fit_partitioned_model(...
%    inputs, outputs, selectors, sel_sizes, min_size, partition_names, fn_name, varargin)
%
% Fit models to different subsets (columns)  of the input/output data,
% as chosen by the selectors matrix. If there is only output data, set input=[].
% If there is less than min_size data in partition i,
% we set model{i} = []
%
% Example:
% selectors = [1 2 1 1 1
%              1 2 2 1 2]
% sel_sizes = [2 2] so there are 4 models: (1,1), (2,1), (1,2), (2,2)
% We fit model{1} to data from columns 1,4
% We fit model{2} to no data
% We fit model{3} to data from column 3,5
% We fit model{4} to data from column 2 (assuming min_size <= 1)
%
% For each partition, we call the specified function with the specified arguments
% as follows:
%  model{i}  = fn(input(:,cols{i}), output(:,cols{i}), args)
% (We omit input if [])
% partition_size(i) is the amount of data in the i'th partition.
%
% Example use: row 1 of selectors is whether an object is present/absent
% and row 2 is the location.
%
% Demo:
% inputs = 1:5; outputs = 6:10; selectors = as above
% fn = 'fit_partitioned_model_testfn';
% [model, partition_size] = fit_partitioned_model(inputs, outputs, selectors, [2 2], fn)
% should produce
% model{1}.input = [1 4], model{1}.output = [6 9]
% model{2} = []
% model{3}.input = [3 5], model{3}.output = [8 10],
% model{4}.input = [2], model{3}.output = [7],
% partition_size = [2 0 2 1]

sel_ndx = subv2ind(sel_sizes, selectors');
Nmodels = prod(sel_sizes);
model = cell(1, Nmodels);
partition_size = zeros(1, Nmodels);
for m=1:Nmodels
  ndx = find(sel_ndx==m);
  partition_size(m) = length(ndx);
  if ~isempty(partition_names) % & (partition_size(m) < min_size)
    fprintf('partition %s has size %d, min size = %d\n', ...
	    partition_names{m}, partition_size(m), min_size);
  end
  if partition_size(m) >= min_size
    if isempty(inputs)
      model{m} = feval(fn_name, outputs(:, ndx), varargin{:});
    else
      model{m} = feval(fn_name, inputs(:,ndx), outputs(:, ndx), varargin{:});
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%CODEGENERATOR.GENSLBLOCKINERTIA Generate Simulink block for inertia matrix
%
% cGen.genslbgenslblockinertia() generates a robot-specific Simulink block to compute
% robot inertia matrix.
%
% Notes::
% - Is called by CodeGenerator.geninertia if cGen has active flag genslblock
% - The Inertia matrix is stored row by row to avoid memory issues.
% - The Simulink block recombines the the individual blocks for each row.
% - The Simulink blocks are generated and stored in a robot specific block
%   library cGen.slib in the directory cGen.basepath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.geninertia.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function genslblockinertia(CGen)

%% Open or create block library
bdclose('all')                                                              % avoid problems with previously loaded libraries
load_system('simulink');
if ~(exist([CGen.slibpath,simulinkext]) == 2)                                  % Create new block library if none exists
 CGen.createnewblocklibrary;
end
open_system(CGen.slibpath);
set_param(CGen.slib,'lock','off');

q = CGen.rob.gencoords;

%% Generate Inertia Block
CGen.logmsg([datestr(now),'\tGenerating Simulink Block for the robot inertia matrix\n']);
nJoints = CGen.rob.n;

CGen.logmsg([datestr(now),'\t\t... enclosing subsystem ']);
symname = 'inertia';
InertiaBlock = [CGen.slib,'/',symname];
    if ~isempty(find_system(CGen.slib,'SearchDepth',1,'Name',symname))                    % Delete previously generated inertia matrix block
        delete_block(InertiaBlock);
        save_system;
    end
% Subsystem in which individual rows are concatenated
add_block('built-in/SubSystem',InertiaBlock);                               % Add new inertia matrix block
add_block('simulink/Math Operations/Matrix Concatenate'...
    , [InertiaBlock,'/inertia']...
    , 'NumInputs',num2str(nJoints)...
    , 'ConcatenateDimension','1');
add_block('simulink/Sinks/Out1',[InertiaBlock,'/out']);
add_block('simulink/Sources/In1',[InertiaBlock,'/q']);
add_line(InertiaBlock,'inertia/1','out/1');
CGen.logmsg('\t%s\n',' done!');

for kJoints = 1:nJoints
    CGen.logmsg([datestr(now),'\t\t... Embedded Matlab Function Block for joint ',num2str(kJoints),': ']);

    % Generate Embedded Matlab Function block for each row
    symname = ['inertia_row_',num2str(kJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    if exist(fname,'file')
        tmpStruct = load(fname);
    else
        error ('genslblockinertia:SymbolicsNotFound','Save symbolic expressions to disk first!')
    end

    blockaddress = [InertiaBlock,'/',symname];
    if doesblockexist(CGen.slib,symname)
        delete_block(blockaddress);
        save_system;
    end

    CGen.logmsg('%s',' block creation');
    symexpr2slblock(blockaddress,tmpStruct.(symname),'vars',{q});

    % connect output
    CGen.logmsg('%s',', output wiring');
    if ( verLessThan('matlab','7.11.0.584') ) && ( isequal(tmpStruct.(symname),zeros(1,nJoints)) )
        % There is a bug in earlier Matlab versions. If the symbolic
        % vector is a zero vector, then the Simulink Embedded Matlab
        % Function block outputs a scalar zero. We need to concatenate
        % a row vector of zeros here, which we have to construct on our
        % own.
        add_block('simulink/Math Operations/Matrix Concatenate'...              % Use a matrix concatenation block ...
            , [InertiaBlock,'/DimCorrection',num2str(kJoints)]...               % ... named with the current row number ...
            , 'NumInputs',num2str(nJoints),'ConcatenateDimension','2');         % ... intended to concatenate zero values for each joint ...
        % ... columnwise. This will circumvent the bug.

        for iJoints = 1:nJoints                                                 % Connect signal lines from the created block (which outputs
            add_line(InertiaBlock...                                            % a scalar zero in this case) with the bugfix block.
                , [symname,'/1']...
                , ['DimCorrection',num2str(kJoints),'/', num2str(iJoints)]);
        end

        add_line(InertiaBlock,['DimCorrection',num2str(kJoints)...              % Connect the fixed row with other rows.
            , '/1'],['inertia/', num2str(kJoints)]);

    else
        add_line(InertiaBlock,[symname,'/1']...          % In case that no bug occurs, we can just connect the rows.
            , ['inertia/', num2str(kJoints)]);
    end

    % Connect inputs
    CGen.logmsg('%s',', input wiring');
    add_line(InertiaBlock,'q/1',[symname,'/1']);
    CGen.logmsg('\t%s\n','row complete!');
end
addterms(InertiaBlock); % Add terminators where needed
distributeblocks(InertiaBlock);
CGen.logmsg([datestr(now),'\tInertia matrix block complete\n']);

%% Built inverse Inertia matrix block
CGen.logmsg([datestr(now),'\tGenerating Simulink Block for the inverse robot inertia matrix\n']);
CGen.logmsg([datestr(now),'\t\t... enclosing subsystem ']);
% block address
symname = 'invinertia';
invInertiaBlock = [CGen.slib,'/',symname];
% remove any existing blocks
    if ~isempty(find_system(CGen.slib,'SearchDepth',1,'Name',symname))                    % Delete previously generated block
        delete_block(invInertiaBlock);
        save_system;
    end
add_block('built-in/SubSystem',invInertiaBlock);
CGen.logmsg('\t%s\n',' done!');

% matrix inversion
CGen.logmsg([datestr(now),'\t\t... matrix inversion block ']);
add_block('simulink/Math Operations/Product',[invInertiaBlock,'/inverse']); % Use a product block...
set_param([invInertiaBlock,'/inverse'],'Inputs','/');                       % ... with single input '/'...
set_param([invInertiaBlock,'/inverse'],'Multiplication','Matrix(*)');       % ... and matrix multiplication
CGen.logmsg('\t%s\n',' done!');

% wire the input and output
CGen.logmsg([datestr(now),'\t\t... input and output ']);
add_block(InertiaBlock,[invInertiaBlock,'/inertiaMatrix']);
add_block('simulink/Sources/In1',[invInertiaBlock,'/q']);
add_block('simulink/Sinks/Out1',[invInertiaBlock,'/out']);

% wire the blocks among each other
CGen.logmsg('and internal wiring');
add_line(invInertiaBlock,'q/1','inertiaMatrix/1');
add_line(invInertiaBlock,'inertiaMatrix/1','inverse/1');
add_line(invInertiaBlock,'inverse/1','out/1');
CGen.logmsg('\t%s\n',' done!');

% add terminators where necessary
addterms(invInertiaBlock);
distributeblocks(invInertiaBlock);
CGen.logmsg([datestr(now),'\tInverse inertia matrix block complete.\n']);

%% Cleanup
% Arrange blocks
distributeblocks(CGen.slib);

% Lock, save and close library
set_param(CGen.slib,'lock','on');
save_system(CGen.slib,CGen.slibpath);
close_system(CGen.slib);

end

author: petercorke

repo: robotics-toolbox-matlab


text:
classdef EventTermCondDirectionEnum < matlab.mixin.SetGet
    %EventTermCondDirectionEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        NoDir('Either Increasing or Decreasing',0)
        Increasing('Increasing at Value',1)
        Decreasing('Decreasing at Value',-1)
    end

    properties
        name char = '';
        direction(1,1) double = 0;
    end

    methods
        function obj = EventTermCondDirectionEnum(name,direction)
            obj.name = name;
            obj.direction = direction;
        end
    end

    methods(Static)
        function listBoxStr = getListboxStr()
            m = enumeration('EventTermCondDirectionEnum');
            listBoxStr = {m.name};
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('EventTermCondDirectionEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function vw = percentTSeries(vw, scanNum, sliceNum, detrend,...
    inhomoCorrection, temporalNormalization, noMeanRemove)
%Convert raw time series into percentage change from mean.
% vw = percentTSeries(vw, [scanNum], [sliceNum], [detrend], ...
%           [inhomoCorrection], [temporalNormalization], [noMeanRemove])
%
% Checks the tSeriesScan and tSeriesSlice slots to see if the
% desired tSeries is already loaded. If so, don't do anything.
% Otherwise:
% 1) loads tSeries corresponding to scanNum/sliceNum.
% 2) removes the DC and baseline trend of a tSeries.
% 3) sets:
%    vw.tSeries = resulting percent tSeries
%    vw.tSeriesScan = scanNum
%    vw.tSeriesSlice = sliceNum
%
% INPUTS
% ------
%   vw:             mrVista view structure
%   scanNum:        integer scalar for scan number
%                       [default = viewGet(vw, 'current scan')]
%   sliceNum:       integer scalar or vector for slice number
%                       [default = viewGet(vw, 'current slice')]
%                       0: all slices: [1:viewGet(vw, 'num slices')]
%   detrend:        Options for how to remove the baseline:
%                       0: no trend removal
%                       1: highpass trend removal
%                       2: quadratic removal
%                      -1: linear trend removal
%                       [default: detrend = detrendFlag(vw,scanNum)]
%   inhomoCorrection: How to compensate for distance from the coil:
%                       0 do nothing
%                       1 divide by the mean, independently at each voxel
%                       2 divide by null condition
%                       3 divide by anything you like, e.g., robust
%                         estimate of intensity inhomogeneity
%                       For inhomoCorrection=3, you must compute the
%                           spatial gradient (from the Analysis menu) or
%                           load a previously computed spatial gradient
%                           (from the File/Parameter Map menu).
%   temporalNormalization: Boolean flag. If true, detrend each frame slice
%                           by slice so to that each frame has the same
%                           mean intensity. See doTemporalNormlalization
%                           [default = false]
%   noMeanRemove:       Hmm. Can someone explain this?
%
% OUTPUTS
% -------
%  vw:  Modified vw structure. The following fields get set in the vw:
%               tSeries, tSeriesScan, tSeriesSlice
%
% EXAMPLE
% -------
% dFolder = mrtInstallSampleData('functional', 'mrBOLD_01');
% cd(dFolder);
% vw = initHiddenInplane;
% vw = percentTSeries(vw, 1, 0);

% EDIT HISTORY:
% djh, 1/22/98
% arw, 12/05/99 Added option to remove quadratic function
% dbr, 8/1/00 Added high-pass baseline removal option
% dbr, 11/16/00 Made high-pass trend removal the default (detrendFlag = 1).
%               Linear trend removal is now detrendFlag = -1.
% djh, 11/00  Added option of dividing by spatialGradient
%             (estimate of intensity inhomogeneity) instead
%             of dividing by mean at each pixel.
% djh, 2/2001 Updated to mrLoadRet-3.0
%             Detrending is now done in detrendTSeries.m
%             This function now sets vw.tSeries (loadtSeries used to do this)
% djh, 9/28/2001 Subtract the mean (again) near the end to make sure it's zero
%             Otherwise, it messes up the correlation map.
% djh, 7/12/2002 Changed the options for inhomogeneity correction
%             Used to have only two options (0: divide by mean; 1: divide by robust est)
%             In the current code, current option 1 is the same as what used to be 0
%             and current option 2 is what used to be.
% dhb, 6/3/2003  Comment out code that checks slice and scan numbers and
%             assumes tSeries is cached if they match passed values.  This
%             check did not seem to be bulletproof.
% ras, 3/8/2007		reverted some changes suggested by Mark Schira about
%			  making singles. I agree to use single precision, but will
%			  make the change in loadtSeries and savetSeries, so that this
%			  code doesn't modify the input data type.

%% Argument check
if notDefined('scanNum'),               scanNum               = viewGet(vw,'current scan'); end
if notDefined('sliceNum'),              sliceNum              = viewGet(vw,'current slice'); end
if notDefined('detrend'),               detrend               = detrendFlag(vw,scanNum); end
if notDefined('inhomoCorrection'),      inhomoCorrection      = inhomoCorrectionFlag(vw,scanNum); end
if notDefined('temporalNormalization'), temporalNormalization = 0; end
if notDefined('noMeanRemove'),          noMeanRemove          = 0; end

if sliceNum == 0, sliceNum = 1:viewGet(vw, 'num slices', scanNum); end

%%
% load tSeries
tSeries = loadtSeries(vw, scanNum, sliceNum);

% also, if the tSeries is empty, return w/o erroring
if isempty(tSeries)
	vw = viewSet(vw, 'tSeries', []);
	vw = viewSet(vw, 'tSeriesScan', scanNum);
	vw = viewSet(vw, 'tSeriesSlice', sliceNum);
	return
end

nFrames = size(tSeries,1);

% Added by ARW
if (temporalNormalization)
	disp('Temporal normalization to first frame');
    tSeries=doTemporalNormalization(tSeries);
end

% Make the mean of all other frames the same as this.
% Divide by either the mean or the spatial gradient
%
switch inhomoCorrection
	case 0
		ptSeries = tSeries;
	case 1
		dc = nanmean(tSeries);
		dc(dc==0 | isnan(dc)) = Inf;  % prevents divide-by-zero warnings
		ptSeries = bsxfun(@rdivide, tSeries, dc);
	case 2
		myErrorDlg('Inhomogeneity correction by null condition not yet implement');
	case 3
		if ~isfield(vw, 'spatialGrad') || isempty(vw.spatialGrad)
			try
				vw = loadSpatialGradient(vw);
				updateGlobal(vw); % make sure it stays loaded
			catch                 %#ok<CTCH>
				myErrorDlg(['No spatial gradient map loaded. Either load '...
					'the spatial gradient map from File menu or edit ' ...
					'dataTypes to set inhomoCorrect = 0 or 1']);
			end
		end
		gradientImg = vw.spatialGrad{scanNum}(:,:,sliceNum);
		dc = gradientImg(:)';
		ptSeries = tSeries./(ones(nFrames,1)*dc);
	otherwise
		myErrorDlg(['Invalid option for inhomogeneity correction: ',num2str(inhomoCorrection)]);
end

% Remove trend
%
if detrend
	ptSeries = detrendTSeries(ptSeries,detrend,detrendFrames(vw,scanNum));
end

% Subtract the mean
% Used to just subtract 1 under the assumption that we had already divided by
% the mean, but now with the spatialGrad option the mean may not be exactly 1.
%
if noMeanRemove==0
	ptSeries = bsxfun(@minus, ptSeries, mean(ptSeries));
	% Multiply by 100 to get percent
	%
	ptSeries = 100*ptSeries;
end

% Set fields in view structure
vw = viewSet(vw, 'tSeries', ptSeries);
vw = viewSet(vw, 'tSeriesScan', scanNum);
vw = viewSet(vw, 'tSeriesSlice', sliceNum);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [sol,val] = coranaMin(sol,options)
% function [val,sol] = coranaMin(sol,options)
%
% Function to minimize the Corana function.
%
% val - the value of the Corana function at point sol
% sol - the location to evaluate the Corana function

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

numVar=size(sol,2)-1;
val = coranaEval(sol(1:numVar));
val = -val;

author: HuangCongQing


text:
% .*  Array multiply.
%    X.*Y denotes element-by-element multiplication. X and Y must have
%    compatible sizes. In the simplest cases, they can be the same size or
%    one can be a scalar. Two inputs have compatible sizes if, for every
%    dimension, the dimension sizes of the inputs are either the same or one
%    of them is 1.
%
%    C = TIMES(A,B) is called for the syntax 'A .* B' when A or B is an
%    object.
%
%    See also MTIMES.
%
%    Reference page in Doc Center
%       doc times
%
%    Other functions named times
%
%       calendarDuration/times    fints/times       tall/times
%       categorical/times         gpuArray/times    timeseries/times
%       codistributed/times       sym/times         ts/times
%       duration/times
%

author: jmaih


text:
function mu_pq = moment_central ( n, x, y, p, q )

%*****************************************************************************80
%
%% MOMENT_CENTRAL computes central moments of a polygon.
%
%  Discussion:
%
%    The central moment Mu(P,Q) is defined by
%
%      Mu(P,Q) = Integral ( polygon ) (x-Alpha(1,0))^p (y-Alpha(0,1))^q dx dy
%              / Area ( polygon )
%
%    where
%
%      Alpha(1,0) = Integral ( polygon ) x dx dy / Area ( polygon )
%      Alpha(0,1) = Integral ( polygon ) y dx dy / Area ( polygon )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carsten Steger,
%    On the calculation of arbitrary moments of polygons,
%    Technical Report FGBV-96-05,
%    Forschungsgruppe Bildverstehen, Informatik IX,
%    Technische Universitaet Muenchen, October 1996.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real X(N), Y(N), the vertex coordinates.
%
%    Input, integer P, Q, the indices of the moment.
%
%    Output, real MU_PQ, the unnormalized moment Mu(P,Q).
%
  alpha_10 = moment_normalized ( n, x, y, 1, 0 );
  alpha_01 = moment_normalized ( n, x, y, 0, 1 );

  mu_pq = 0.0;

  for i = 0 : p
    for j = 0 : q

      alpha_ij = moment_normalized ( n, x, y, i, j );

      mu_pq = mu_pq + r8_mop ( p + q - i - j ) ...
        * r8_choose ( p, i ) * r8_choose ( q, j ) ...
        * alpha_10 ^ ( p - i ) * alpha_01 ^ ( q - j ) * alpha_ij;

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mtPaths_ = dtiLoadMetrotracPathsFromStrVer3(str,i, numstats, mtPaths,xform)
% Loads version 3 pdb file
%
%   mtPaths_ = dtiLoadMetrotracPathsFromStrVer3(str,i, numstats,mtPaths,xform)
%
% str is the string buffer
% i is the current position in the string buffer
% numstats is the number of statistics in the pdb
% It appears that xform is never used.
%
% HISTORY:
% 2010.08.01
%
% Stanford VISTA
%
%%

% profile on

% Boolean.  Only Matlab transform what?  Appears to always be false and
% thus the code below is not used.

% How many bytes for int32 and double
int32size   = size(typecast(int32(0),'uint8'),2);
doublesize = size(typecast(double(0.0),'uint8'),2);

numpaths = typecast( uint8(str(i:i+int32size-1)),'uint32');
i = i + int32size;
points_per_fiber = typecast( uint8(str(i:i+int32size*numpaths-1)),'int32'); i = i + int32size*numpaths;

% compute the total number of points
total_pts = 0;
for k = 1:numpaths
    total_pts = total_pts + points_per_fiber(k);
end
total_pts = uint32(total_pts);

% read all the fiber points
fiber_points = typecast( uint8(str(i:i+doublesize*total_pts*3-1)),'double');
i = i + doublesize*total_pts*3;
points_read = 0;

%% Reading each pathway
for p = 1:numpaths
    numnodes = points_per_fiber(p);
    mtPaths.pathwayInfo(p).algo_type = 0; % algo type
    mtPaths.pathwayInfo(p).seed_point_index = 1; % seedpointindex ???

    % Reading path positions
    mtPaths.pathways{p,1} = fiber_points( 1+points_read*3 : (points_read+numnodes)*3);
    mtPaths.pathways{p,1} = reshape(mtPaths.pathways{p,1},3,numnodes);
    points_read = points_read + numnodes;

    if mod(p,10000)==0
        disp(['Loaded ' num2str(p) ' of ' num2str(numpaths)]);
    end
end

% Transform all fiber coords in one go (much faster!)
% I don't know why this offset is still necessary, isn't ACPC space all
% ACPC space??
bOnlyMatlabXForm=0;   % It appears that this is always false and could be deleted.
if( ~bOnlyMatlabXForm )
    xformToMatlab = eye(4);
    xformToAcpc = xformToMatlab*mtPaths.xform;
    if isfield(mtPaths,'pathways')
        mtPaths.pathways = dtiXformFiberCoords(mtPaths.pathways, xformToAcpc);
    else
        mtPaths.pathways = [];
        mtPaths.pathwayInfo = [];
    end
end

%% Read per fiber stats
for as = 1:numstats
    per_fiber_stat = typecast( uint8(str(i:i+doublesize*numpaths-1)),'double');
    i = i + doublesize*numpaths;
    for p = 1:numpaths
        mtPaths.pathwayInfo(p).pathStat(as) = per_fiber_stat(p);

    end
end
%% Read per point stats
for as = 1:numstats
    if(mtPaths.statHeader(as).is_computed_per_point)
        per_point_stat = typecast( uint8(str(i:i+doublesize*total_pts-1)),'double');
        i = i + doublesize*total_pts;
    end

    points_read = 0;
    % Assign per point stats value
    for p = 1:numpaths
        if(mtPaths.statHeader(as).is_computed_per_point)
            mtPaths.pathwayInfo(p).point_stat_array(as,:) = ...
                per_point_stat(points_read+1 : points_read+points_per_fiber(p));
            points_read = points_read + points_per_fiber(p);
        end
    end
end
mtPaths_ = mtPaths;

% profile off

end


text:
function out = colorMatchingFunctions
%colorMatchingFunctions Color-matching functions.
%   colorMatchingFunctions returns a table containing the color matching
%   functions for the CIE 1931 Standard Observer. The table contains the
%   variables lambda, x, y, and z. The color matching functions are
%   returned for wavelengths between 360nm and 830nm at intervals of 1nm.

%   Written by Steve Eddins to accompany Digital Image Processing Using
%   MATLAB, 3rd edition, Gatesmark Press, 2020,
%   http://imageprocessingplace.com.
%
%   Copyright 2019 The MathWorks, Inc.
%   License: https://github.com/mathworks/matlab-color-tools/blob/master/license.txt

persistent t

if isempty(t)
   t = readtable('CIE 1931 Standard Observer (1nm).xlsx',...
      'UseExcel',false);
end

out = t;

author: dipum

repo: dipum-toolbox

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


text:
function out = isSubset(A, B, tol)
%ISSUBSET   Test if A is a subset of B.
%   OUT = ISSUBSET(A, B, TOL) returns logical true if A is a subset of B up
%   to the tolerance TOL, where A and B are domains of a chebfun, chebfun2
%   or chebfun3 (i.e., vectors of size 1x2, 1x4 or 1x6 with A(1) <= A(2),
%   A(3) <= A(4), A(5) <= A(6)).

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

if ( isempty(A) )
    out = 1;
    return
end

if ( length(A) ~= length(B) )
    error('CHEBFUN:isSubset:size', ...
        'Domains must have the same number of entries.')
end

if ( ( A(1) < B(1) - tol ) || ( B(2) + tol < A(2) ) )
    out = 0;
    return
else
    % Recurse:
    out = isSubset(A(3:end), B(3:end), tol);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ x, w ] = hermite_ek_compute ( n )

%*****************************************************************************80
%
%% HERMITE_EK_COMPUTE computes a Gauss-Hermite quadrature rule.
%
%  Discussion:
%
%    The code uses an algorithm by Elhay and Kautsky.
%
%    The abscissas are the zeros of the N-th order Hermite polynomial.
%
%    The integral:
%
%      integral ( -oo < x < +oo ) exp ( - x * x ) * f(x) dx
%
%    The quadrature rule:
%
%      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2011
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the number of abscissas.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%

%
%  Define the zero-th moment.
%
  zemu = gamma ( 0.5 );
%
%  Define the Jacobi matrix.
%
  bj = zeros ( n, 1 );
  for i = 1 : n
    bj(i) = i / 2.0;
  end
  bj(1:n) = sqrt ( bj(1:n) );

  x = zeros ( n, 1 );

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  w(1:n) = w(1:n).^2;

  return
end


text:
% PURPOSE: Loads current ERP structure (if any) from de the workspace.
%          Otherwise, load ALLERP(CURRENTERP); Otherwise ERP = [];
%
%          To avoid clearing an already filled ERP structure after an interrupted
%          process (for instance, after an error)
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

function ALLERP = preloadALLERP
try
        ALLERP = evalin('base', 'ALLERP');
catch
        ALLERP = [];
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [lores hires] = collectSamplesScales(conf, ohires, numscales, scalefactor)

lores = [];
hires = [];

for scale = 1:numscales
    sfactor = scalefactor^(scale-1);
    chires = resize(ohires, sfactor, 'bicubic');

    chires = modcrop(chires, conf.scale); % crop a bit (to simplify scaling issues)
    % Scale down images
    clores = resize(chires, 1/conf.scale, conf.interpolate_kernel);
    midres = resize(clores, conf.upsample_factor, conf.interpolate_kernel);
    features = collect(conf, midres, conf.upsample_factor, conf.filters);
    clear midres

    interpolated = resize(clores, conf.scale, conf.interpolate_kernel);
    clear clores
    patches = cell(size(chires));
    for i = 1:numel(patches) % Remove low frequencies
        patches{i} = chires{i} - interpolated{i};
    end
    clear chires interpolated

    hires = [hires collect(conf, patches, conf.scale, {})];

    lores = [lores conf.V_pca' * features];
end

author: thomas-koehler

repo: SupER


text:
function updatelarvasamuel(hfly,hfly_extra,pos)

set(hfly,'XData',[pos.xcontour,nan,pos.xspine'],'YData',[pos.ycontour,nan,pos.yspine']);
set(hfly_extra,'XData',[pos.xhead,pos.xmid],'YData',[pos.yhead,pos.ymid]);

author: kristinbranson

repo: JAABA


text:
% OtQueryNew Connect to an Omnitrans query.
%
% Usage:
%    hquery = OtQueryNew()
% or:
%    hquery = OtQueryNew(<SQL statement>)
%    This form includes a call to OtQuerySetSQL()
%
% To release a query handle, call OtQueryFree().

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-02-28 Created (Johan Meijdam, VORtech)

function hquery = OtQueryNew(statement)
   if nargin == 1
      hquery = Ot('query', statement);
   else
      hquery = Ot('query');
   end
end

author: Sable


function G = get_payoff_G_matrix_from_ygrid_3d( y_1, y_2, y_3, S_0s, sigmas, R, contractParams)
%UNTITLED5 Summary of this function goes here
%   Detailed explanation goes here

payoff_type = contractParams.payoff_type;

rho12 = R(1,2);
rho23 = R(2,3);
rho13 = R(1,3);

gamma = (rho12*rho13 - rho23)/(1 - rho12^2);

if payoff_type == 1 || payoff_type == 2   % G = S_1, or G=S_2, or G=S_3
    dim = contractParams.dim;
    if dim == 1
        payoff = @(y1,y2,y3)S_0s(1)*exp(sigmas(1)*y1);
    elseif dim == 2
        payoff = @(y1,y2,y3)S_0s(2)*exp(sigmas(2)*(y2 + rho12*y1));
    else
        payoff = @(y1,y2,y3)S_0s(3)*exp(sigmas(3)*(rho13*y1 -gamma*y2 + y3));
    end

elseif payoff_type == 5  % Geometric Basket Call / Put:  G = (sqrt(S_1) * sqrt(S_2) - K)^+ and
    K = contractParams.K;
    if contractParams.call == 1
        payoff = @(y1,y2,y3) max(0, (S_0s(1)*exp(sigmas(1)*y1) * S_0s(2)*exp(sigmas(2)*(y2 + rho12*y1))* S_0s(3)*exp(sigmas(3)*(rho13*y1 -gamma*y2 + y3)))^(1/3) - K);
    else
        payoff = @(y1,y2,y3) max(0, K - (S_0s(1)*exp(sigmas(1)*y1) * S_0s(2)*exp(sigmas(2)*(y2 + rho12*y1))* S_0s(3)*exp(sigmas(3)*(rho13*y1 -gamma*y2 + y3)))^(1/3));
    end

elseif payoff_type == 6  % Arithmetic Basket Call / Put:  G = (sqrt(S_1) * sqrt(S_2) - K)^+
    K = contractParams.K;
    if contractParams.call == 1
        payoff = @(y1,y2,y3) max(0, (1/3)*(S_0s(1)*exp(sigmas(1)*y1) + S_0s(2)*exp(sigmas(2)*(y2 + rho12*y1)) + S_0s(3)*exp(sigmas(3)*(rho13*y1 -gamma*y2 + y3))) - K);
    else
        payoff = @(y1,y2,y3) max(0, K - (1/3)*(S_0s(1)*exp(sigmas(1)*y1) + S_0s(2)*exp(sigmas(2)*(y2 + rho12*y1)) + S_0s(3)*exp(sigmas(3)*(rho13*y1 -gamma*y2 + y3))));
    end

end

m_0 = length(y_1);
G = zeros(m_0, m_0, m_0);

for i=1:m_0
    for j=1:m_0
        for k=1:m_0
            G(i,j,k) = payoff(y_1(i), y_2(j), y_3(k));
        end
    end
end

end


text:
clear all
clc

iaf.designation='23012';
iaf.n=56;
iaf.HalfCosineSpacing=1;
iaf.wantFile=1;
iaf.datFilePath='./'; % Current folder
iaf.is_finiteTE=0;

af = naca5gen(iaf);

% set(af.hAF,'Marker','none')
hold on

plot(af.x,af.z,'bo-')

plot(af.xU,af.zU,'bo-')
plot(af.xL,af.zL,'ro-')

plot(af.xC,af.zC,'r--')

title(af.name)
axis equal


text:
function [ grSOrtho, H ] = msfmGradientLsml( anim, d )
% Compute the LSML gradient in MSFM from Rabaud CVPR08
%
% Just check Rabaud's CVPR08 paper
% the notations follow Rabaud's PhD thesis: Appendix D3
%
% USAGE
%  [ grSOrtho, H ] = msfmGradientLsml( anim, d )
%
% INPUTS
%  anim          - Animation object
%  d             - dimensionality of the shape manifold
%
% OUTPUTS
%  grSOrtho      - LSML gradient, orthogonal to the manifold
%  H             - [ 3nPoint x d x nFrame ] tangent basis at every point
%                  of the manifold
%
% EXAMPLE
%
% See also MSFMGRADIENTSRT
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

nFrame = anim.nFrame; nPoint = anim.nPoint;

X3 = reshape(anim.S,3*nPoint, nFrame);

% figure out neighbors in the data we have
neigh = computeNeighbor( X3,'k',3,'forceConn',2);
Nmat = neigh.Nmat;

% force temporally closeby frames to be considered together
row = zeros(1,nFrame); row(2) = 1;
Nmat(logical(toeplitz(row))) = 1;

% learn the manifold from here
pTh=struct('d',d,'rbfK',5,'nRestart',5,'nSamples',250,'show',0);
manifold=lsmlOptimizeTh(X3,Nmat,pTh);

% get the gradient for denoising
[Chi3, Chi] = lsmlDenoise( X3, Nmat, manifold, 'nItr', 1 );

grSOrtho = reshape(X3-Chi,3,nPoint,nFrame);

if nargout >=2; H = lsmlComputeH( X3, manifold, 1 ); end
end

author: vrabaud


text:
function [x, y, z, varargout] = stlread(filename)
% This function reads an STL file in binary format into matrixes X, Y and
% Z, and C.  C is optional and contains color rgb data in 5 bits.
%
% USAGE: [x, y, z, c] = stlread(filename);
%
% To plot use patch(x,y,z,c), or patch(x,y,z)
%
% Written by Doron Harlev

if nargout>4
    error('Too many output arguments')
end
use_color=(nargout==4);

fid=fopen(filename, 'r'); %Open the file, assumes STL Binary format.
if fid == -1
    error('File could not be opened, check name or path.')
end

ftitle=fread(fid,80,'uchar=>schar'); % Read file title
num_facet=fread(fid,1,'int32'); % Read number of Facets

%fprintf('\nTitle: %s\n', char(ftitle'));
%fprintf('Num Facets: %d\n', num_facet);

% Preallocate memory to save running time
x=zeros(3,num_facet); y=zeros(3,num_facet); z=zeros(3,num_facet);
if use_color
    c=uint8(zeros(3,num_facet));
end

%h = waitbar(0,'Please wait...');
for i=1:num_facet,
    norm=fread(fid,3,'float32'); % normal coordinates, ignored for now
    ver1=fread(fid,3,'float32'); % vertex 1
    ver2=fread(fid,3,'float32'); % vertex 2
    ver3=fread(fid,3,'float32'); % vertex 3
    col=fread(fid,1,'uint16'); % color bytes
    if (bitget(col,16)==1 & use_color)
        r=bitshift(bitand(2^16-1, col),-10);
        g=bitshift(bitand(2^11-1, col),-5);
        b=bitand(2^6-1, col);
        c(:,i)=[r; g; b];
    end
    x(:,i)=[ver1(1); ver2(1); ver3(1)]; % convert to matlab "patch" compatible format
    y(:,i)=[ver1(2); ver2(2); ver3(2)];
    z(:,i)=[ver1(3); ver2(3); ver3(3)];
    if mod(i,floor(num_facet/10))==0
        %waitbar(i/num_facet,h);
    end
end
if use_color
    varargout(1)={c};
end
fclose(fid);
%close(h);

% For more information http://rpdrc.ic.polyu.edu.hk/old_files/stl_binary_format.htm

author: alecjacobson

repo: gptoolbox


text:
function a = nnz(t)
%NNZ Number of nonzeros in sparse tensor.
%
%   NNZ(T) is the number of nonzero elements in T.
%
%   See also SPTENSOR, SPTENSOR/FIND.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isempty(t.subs)
    a = 0;
else
    a = size(t.subs,1);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%run this test case with the command
%results = runtests('modelSortingTests.m')
function tests = modelSortingTests
tests = functiontests(localfunctions);
end

function sortIdentifirs_and_permuteModelTest(testCase)

%Load the expected (i.e. sorted) model
sourceDir = fileparts(which(mfilename));
load([sourceDir,'/test_data/ecoli_textbook.mat'], 'model');
expModel = model;

%Create the actual model that will be permuted and sorted
actModel = expModel;

%Randomly permutate model, do not use RAVEN functions
rndIdx = randperm(numel(actModel.rxns));
fieldsToChange = {'rxns','lb','ub','rev','c','rxnNames','grRules','eccodes'};
for i=1:numel(fieldsToChange)
    actModel.(fieldsToChange{i}) = actModel.(fieldsToChange{i})(rndIdx);
end
actModel.S          = actModel.S(:,rndIdx);
actModel.rxnGeneMat = actModel.rxnGeneMat(rndIdx,:);

rndIdx = randperm(numel(actModel.mets));
fieldsToChange = {'mets','metNames','metComps','metFormulas','metMiriams'};
for i=1:numel(fieldsToChange)
    actModel.(fieldsToChange{i}) = actModel.(fieldsToChange{i})(rndIdx);
end
actModel.S     = actModel.S(rndIdx,:);

rndIdx = randperm(numel(actModel.genes));
fieldsToChange = {'genes','geneShortNames'};
for i=1:numel(fieldsToChange)
    actModel.(fieldsToChange{i}) = actModel.(fieldsToChange{i})(rndIdx);
end
actModel.rxnGeneMat = actModel.rxnGeneMat(:,rndIdx);

rndIdx = randperm(numel(actModel.comps));
fieldsToChange = {'comps','compNames'};
for i=1:numel(fieldsToChange)
    actModel.(fieldsToChange{i}) = actModel.(fieldsToChange{i})(rndIdx);
end
[~,J]=sort(rndIdx);
[toreplace, bywhat] = ismember(actModel.metComps,1:length(J));
actModel.metComps(toreplace) = J(bywhat(toreplace));

%Sort randomly permutated model
actModel = sortIdentifiers(actModel);

%Check that the actual model is the same as the expected model
verifyEqual(testCase,actModel,expModel)
end

function expandModel_and_contractModelTest(testCase)
%Load the expected model
sourceDir = fileparts(which(mfilename));
load([sourceDir,'/test_data/ecoli_textbook.mat'], 'model');

% Note that this does not work any model, as grRules might not be properly
% sorted (but should still be valid)
evalc('modelNew = expandModel(model);'); % Suppress warnings about complex grRules
modelNew = contractModel(modelNew);

verifyEqual(testCase,model,modelNew)
end

author: SysBioChalmers

repo: RAVEN


text:
function timer_etime_test02 ( )

%*****************************************************************************80
%
%% TEST02 times the vectorized EXP routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2005
%
%  Author:
%
%    John Burkardt
%
  n_log_min = 12;
  n_log_max = 22;
  n_min = 2^n_log_min;
  n_max = 2^n_log_max;
  n_rep = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Time vectorized operations:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    y(1:n) =        x(1:n)  \n' );
  fprintf ( 1, '    y(1:n) = PI *   x(1:n)  \n' );
  fprintf ( 1, '    y(1:n) = sqrt ( x(1:n) )\n' );
  fprintf ( 1, '    y(1:n) = exp  ( x(1:n) )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data vectors will be of minimum size %d\n', n_min );
  fprintf ( 1, '  Data vectors will be of maximum size %d\n', n_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', n_rep );

  for func = 1 : 4

    for i_rep = 1 : n_rep

      for n_log = n_log_min : n_log_max

        n = 2^n_log;

        x = rand(n,1);

        time1 = clock;

        if ( func == 1 )
          y(1:n) = x(1:n);
        elseif ( func == 2 )
          y(1:n) = pi * x(1:n);
        elseif ( func == 3 )
          y(1:n) = sqrt ( x(1:n) );
        elseif ( func == 4 )
          y(1:n) = exp ( x(1:n) );
        end

        time2 = clock;

        delta(n_log,i_rep) = etime ( time2, time1 );

      end

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, 'Test02 Results:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Vector Size  Rep #1        Rep #2        ' );
    fprintf ( 1, 'Rep #3        Rep #4        Rep #5\n' );
    fprintf ( 1, '\n' );

    for n_log = n_log_min : n_log_max
      n = 2^n_log;
      fprintf ( 1, '%10d', n );
      for j = 1 : n_rep
        fprintf ( 1, '%14f', delta(n_log,1:n_rep) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end


text:
% GET THE GLOBAL TRAJECTORY DATA FOR AN INDIVIDUAL OBJECT
function [currentFigure,figureSet] = GetFigure_objectTrajectory(SIM,DATA,currentFigure)

figureSet = [];
for n = 1:DATA.totalObjects
    [currentFigure,figureSet(n)] = BuildObjectTrajectoryFigure(SIM,DATA,currentFigure,n);
end
% ASSEMBLE TABBED FIGURE
windowHandle = GetTabbedFigure(figureSet,'OpenMAS Trajectory Overview');
set(windowHandle,'Position',DATA.figureProperties.windowSettings);          % Maximise the figure in the tab
savefig(windowHandle,[SIM.outputPath,'global-trajectory-overview']);        % Save the output figure
end

% Generate the trajectory figure for the object
function [currentFigure,figureHandle] = BuildObjectTrajectoryFigure(SIM,DATA,currentFigure,objectNum)

% Declare title string for figure
titlestr = sprintf('Global trajectory data for %s over a period of %ss',SIM.OBJECTS(objectNum).name,num2str(SIM.TIME.endTime));

objectLabel = sprintf('[ID-%s] %s',num2str(SIM.OBJECTS(objectNum).objectID),SIM.OBJECTS(objectNum).name);

% CONFIGURE THE PLOT ATTRIBUTES
figurePath  = strcat(SIM.outputPath,sprintf('global-trajectories-%s',sprintf('id-%s-%s',num2str(SIM.OBJECTS(objectNum).objectID),SIM.OBJECTS(objectNum).name)));
figureHandle = figure('Name',objectLabel);
set(figureHandle,'Position', DATA.figureProperties.windowSettings);        % [x y width height]
set(figureHandle,'Color',DATA.figureProperties.figureColor);               % Background colour
set(figureHandle,'Visible','off');
plotCellWidth = 4; plotCellA = 1;                                          % The width of each figure, the start of the plot
setappdata(figureHandle, 'SubplotDefaultAxesLocation', [0.1, 0.1, 0.85, 0.80]);

% EXTRACT TIME-STATE DATA FROM THE TRAJECTORY MATRIX
%[objectStates] = OMAS_getTrajectoryData_mex(DATA.globalTrajectories,SIM.globalIDvector,SIM.OBJECTS(objectNum).objectID,inf);
[objectStates] = OMAS_getTrajectoryData(DATA.globalTrajectories,SIM.globalIDvector,SIM.OBJECTS(objectNum).objectID,inf);

% STATE NAME VECTOR
stateTags = {{'$x$','$(m)$'},{'$y$','$(m)$'},{'$z$','$(m)$'},...
             {'$u$','$(m/s)$'},{'$v$','$(m/s)$'},{'$w$','$(m/s)$'},...
             {'$\phi$','$(rad)$'},{'$\theta$','$(rad)$'},{'$\psi$','$(rad)$'}};
% CONVERT QUATERNION STATE TO EULER ANGLE STATES
eulerStates = zeros(9,size(objectStates,2));
for i = 1:size(objectStates,2)
    eulerStates(1:6,i) = objectStates(1:6,i);
    eulerStates(7:9,i) = OMAS_geometry.quaternionToEulers(objectStates(7:10,i));
end

clear objectStates;

% FOR EACH OBJECTS PERSPECTIVE
stateVectorLength = size(eulerStates,1);
for n = 1:stateVectorLength
    % CALCULATE THE ABSOLUTE SEPERATION FROM ALL OTHER OBJECTS OVER TIME
    plotCellB = n*plotCellWidth;                                    % The end of the plot
    plotLocation = subplot(stateVectorLength,plotCellWidth,[plotCellA plotCellB]);

    % PLOT THE SEPERATION DATA ON CURRENT FIGURE
    plot(plotLocation,DATA.timeVector(1:SIM.TIME.endStep),...
         eulerStates(n,1:SIM.TIME.endStep),...
         'LineStyle','-',...
         'LineWidth',DATA.figureProperties.lineWidth,...
         'Color','b');
    % Y-axes
    ylabel(plotLocation,stateTags{n},...
        'Interpreter',DATA.figureProperties.interpreter,...
        'fontname',DATA.figureProperties.fontName,...
        'Fontweight',DATA.figureProperties.fontWeight,...
        'FontSize',DATA.figureProperties.axisFontSize,...
        'FontSmoothing','on');
    % Axes
    set(plotLocation,...
        'TickLabelInterpreter',DATA.figureProperties.interpreter,...
        'fontname',DATA.figureProperties.fontName,...
        'Fontweight',DATA.figureProperties.fontWeight,...
        'FontSize',DATA.figureProperties.axisFontSize,...
        'FontSmoothing','on',...
        'Color',DATA.figureProperties.axesColor);
    grid on; box on; grid minor;
    plotLocation.XAxis.MinorTickValues = plotLocation.XAxis.Limits(1):SIM.TIME.dt:plotLocation.XAxis.Limits(2);
    % ADD FIGURE TITLE TO FIRST PLOT HEADER
    if n == 1
        % Append title to first subplot
        title(plotLocation,titlestr,...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'Fontsize',DATA.figureProperties.titleFontSize,...
            'FontSmoothing','on');

    end
    % Prevent overlap of x-label
    if n == stateVectorLength
        % X-axis
        xlabel(plotLocation,'t (s)',...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'FontSize',DATA.figureProperties.axisFontSize,...
            'FontSmoothing','on');
    else
        set(plotLocation,'XTickLabel',[]);
    end
    % Define the x-limits
    xlim([SIM.TIME.startTime,SIM.TIME.endTime]);                           % Ensures plot alignment/sizing
    % Move to next subplot location
    plotCellA = plotCellA + plotCellWidth;
end
hold off;

% SAVE THE OUTPUT FIGURE
set(figureHandle,'Visible','on');                % Make it visable for saving
savefig(figureHandle,figurePath);                % As matlab figure
% Publish as .pdf if requested
if DATA.figureProperties.publish
	GetFigurePDF(figureHandle,figurePath);
end
set(figureHandle,'Visible','off');
currentFigure = currentFigure + 1;  % Iterate the plot count
end

repo: OpenMAS


text:
function fx = p06_fun ( n, x )

%*****************************************************************************80
%
%% P06_FUN evaluates the integrand for problem 6.
%
%  Interval:
%
%    -1 <= x <= 1
%
%  Integrand:
%
%    sqrt ( abs ( x + 0.5 ) )
%
%  Exact Integral:
%
%    ( sqrt ( 2 ) + 3 * sqrt ( 6 ) ) / 6 = 1.460447131787105
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Clenshaw, Alan Curtis,
%    A Method for Numerical Integration on an Automatic Computer,
%    Numerische Mathematik,
%    Volume 2, Number 1, December 1960, pages 197-205.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = sqrt ( abs ( x + 0.5 ) );

  return
end


text:
function write_nitf_security_tags(obj)
%WRITE_NITF_SECURITY_TAGS Writes the NITF security tags
%
% It's not expected that the typical user would call this function, rather
% it is called as part of a bigger NITF/SICD writer.
%
% Inputs:
%       obj:       This is the 'hidden' object-specific handle.  It is
%                  similar to 'this' in Java.
%
% References:
%       MIL-STD-2500C, Department of Defense Interface Standard
%                      National Imagery Transmission Format version 2.1
%                      01 May 2006
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Always writing version 2.1 NITF

% Security is limited in the SICD XML structure.  We fill in what we can.
if isfield(obj.sicdmeta,'CollectionInfo')&&isfield(obj.sicdmeta.CollectionInfo,'Classification')
    classification = obj.sicdmeta.CollectionInfo.Classification;
else
    classification = ' ';
end
code_index = regexp(classification,'/[^/]');
if isempty(code_index)
    code = '';
else
    code = classification((code_index(1)+1):end);
end
fwriten(obj.FID, classification(1), 1); % CLAS
fwriten(obj.FID, 'US',      2);  % CLSY
fwriten(obj.FID, code,      11); % CODE
fwriten(obj.FID, '',        2);  % CTLH
fwriten(obj.FID, '',        20); % REL
fwriten(obj.FID, '',        2);  % DCTP
fwriten(obj.FID, '',        8);  % DCDT
fwriten(obj.FID, '',        4);  % DCXM
fwriten(obj.FID, '',        1);  % DG
fwriten(obj.FID, '',        8);  % DGDT
fwriten(obj.FID, ''       , 43); % CLTX
fwriten(obj.FID, '',        1);  % CATP
fwriten(obj.FID, '',        40); % CAUT
fwriten(obj.FID, '',        1);  % CRSN
fwriten(obj.FID, '',        8);  % SRDT
fwriten(obj.FID, '',        15); % CTLN

end

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

author: ngageoint


text:
function []=fdisp(Y,n,T,endo,stringdates2,decimaldates2,Fstartlocation,Fendlocation,forecast_estimates,pref)

% function []=fdisp(Y,n,T,endo,stringdates2,decimaldates2,Fstartlocation,Fendlocation,forecast_estimates,datapath)
% plots the results for unconditional forecasts
% inputs:  - matrix 'Y': matrix of regressands for the VAR model (defined in 1.1.8)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - cell 'endo': list of endogenous variables of the model
%          - cell 'stringdates2': date strings for the sample+forecasts period
%          - vector 'decimaldates2': dates converted into decimal values, for the sample+forecasts period
%          - integer 'Fstartlocation': position of the forecast start date in stringdates2
%          - integer 'Fendlocation': position of the forecast end date in stringdates2
%          - cell 'forecast_estimates': lower bound, point estimates, and upper bound for the unconditional forecasts
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% first, create a cell that will store the data to be plot
plotdata=cell(n,1);

% each cell entry is a matrix of actual and forecast values
% this matrix comprises 4 rows: the first row is actual data, while the three other rows are the estimates (point estimates and confidence bands) for the forecasts
% also, this matrix has a number of rows equal to the dimension of decimaldates2, which comprises the total period sample+forecasts

% loop over variables

for ii=1:n
   plotdata{ii,1}=nan(4,size(decimaldates2,1));
   % record actual sample values
   plotdata{ii,1}(1,1:T)=Y(:,ii)';
   % copy the last point of the actual sample for the forecast part of the matrix (required to have a clean plot)
   plotdata{ii,1}(:,Fstartlocation-1)=repmat(Y(Fstartlocation-1,ii),4,1);
   % record forecast, lower bound
   plotdata{ii,1}(2,Fstartlocation:Fendlocation)=forecast_estimates{ii,1}(1,:);
   % record forecast, point estimate
   plotdata{ii,1}(3,Fstartlocation:Fendlocation)=forecast_estimates{ii,1}(2,:);
   % record forecast, upper bound
   plotdata{ii,1}(4,Fstartlocation:Fendlocation)=forecast_estimates{ii,1}(3,:);
end

if pref.plot

% create forecast figure
% then plot actual vs. fitted
forecast=figure('Tag','BEARresults');
set(forecast,'Color',[0.9 0.9 0.9]);
set(forecast,'name','unconditional forecasts');
ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
for ii=1:n
subplot(nrows,ncolumns,ii);
hold on
Xpatch=[decimaldates2(Fstartlocation-1:Fendlocation,1)' fliplr((decimaldates2(Fstartlocation-1:Fendlocation,1))')];
Ypatch=[plotdata{ii,1}(2,Fstartlocation-1:Fendlocation) fliplr(plotdata{ii,1}(4,Fstartlocation-1:Fendlocation))];
Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
set(Fpatch,'facealpha',0.6);
set(Fpatch,'edgecolor','none');
plot(decimaldates2,plotdata{ii,1}(3,:),'Color',[0.4 0.4 1],'LineWidth',2);
plot(decimaldates2,plotdata{ii,1}(1,:),'Color',[0 0 0],'LineWidth',2);
hold off
set(gca,'XLim',[decimaldates2(1,1) decimaldates2(end,1)],'FontName','Times New Roman');
set(gca,'XGrid','on');
set(gca,'YGrid','on');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal','Interpreter','none');
end

end % pref.plot
% finally, save on excel
bear.data.excelrecord5fcn(stringdates2, endo, plotdata, n, pref)

author: european-central-bank

repo: BEAR-toolbox


text:
function [s vw] = openVolumeWindow
%
% Calls openRawVolumeWindow to set up the VOLUME data structure,
% opens and initializes the volume window.  Then calls switch2Vol.
%
%   [s vw] = openVolumeWindow
%
% djh, 2/99
% Junjie 2004.02.14 let function return s

mrGlobals

s = openRawVolumeWindow;
VOLUME{s} = switch2Vol(VOLUME{s});

selectView(VOLUME{s});

% If user requested a view output, give it to them:
if nargout > 1,	vw = VOLUME{s}; end

return;


text:
%IHOUGH	Hough transform
%
%	params = IHOUGH
%	H = IHOUGH(IM)
%	H = IHOUGH(IM, params)
%
%	Compute the Hough transform of the image IM data.
%
%	The appropriate Hough accumulator cell is incremented by the
%	absolute value of the pixel value if it exceeds
%	params.edgeThresh times the maximum value found.
%
%	Pixels within params.border of the edge will not increment.
%
% 	The accumulator array has theta across the columns and offset down
%	the rows.  Theta spans the range -pi/2 to pi/2 in params.Nth increments.
%	Offset is in the range 1 to number of rows of IM with params.Nd steps.
%
%	Clipping is applied so that only those points lying within the Hough
%	accumulator bounds are updated.
%
%	The output argument H is a structure that contains the accumulator
%	and the theta and offset value vectors for the accumulator columns
%	and rows respectively.  With no output
%	arguments the Hough accumulator is displayed as a greyscale image.
%
%	H.h	the Hough accumulator
%	H.theta	vector of theta values corresponding to H.h columns
%	H.d	vector of offset values corresponding to H.h rows
%
%	For this version of the Hough transform lines are described by
%
%		d = y cos(theta) + x sin(theta)
%
%	where theta is the angle the line makes to horizontal axis, and d is
%	the perpendicular distance between (0,0) and the line.  A horizontal
%	line has theta = 0, a vertical line has theta = pi/2 or -pi/2
%
%	The parameter structure:
%
%	params.Nd number of offset steps (default 64)
%	params.Nth number of theta steps (default 64)
%	params.edgeThresh increment threshold (default 0.1)
%	params.border width of non-incrmenting border(default 8)
%
%
% SEE ALSO: xyhough testpattern isobel ilap

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

function H = ihough(IM, params)

	% process arguments
	%  the param structure is common to ihough()
	if nargin < 2,
		default.Nd = 64;
		default.Nth = 64;
		default.edgeThresh = 0.10;
		default.border = 8;

		% for houghpeaks
		default.houghThresh = 0.40;
		default.radius = 5;
		default.interpWidth = 5;
	end
	if nargin == 0,
		H = default;
		return;
	elseif nargin == 1,
		params = default;
	end

	[nr,nc] = size(IM);

	% find the significant edge pixels
	IM = abs(IM);
	globalMax = max(IM(:));
	i = find(IM > (globalMax*params.edgeThresh));
	[r,c]=ind2sub(size(IM), i);

	xyz = [c r IM(i)];

	% eliminate those near the edge
	k = (c < params.border) | (c>(nc-params.border)) | ...
		(r<params.border) | (r>(nr-params.border));
	xyz(k,:) = [];

	% now pass the x/y/strenth info to xyhough
	H = ihough_xy(xyz, [1 nr params.Nd], params.Nth);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [h,tf]=Jakes_Flat(fd,Ts,Ns,t0,E0,phi_N)
% Inputs:
%   fd      : Doppler frequency
%   Ts      : sampling period
%   Ns      : number of samples
%   t0      : initial time
%   E0      : channel power
%   phi_N  : inital phase of the maximum doppler frequency sinusoid
% Outputs:
%   h       : complex fading vector
%   t_state: current time
if nargin<6
    phi_N=0;
end
if nargin<5
    E0=1;
end
if nargin<4
    t0=0;
end
if nargin<3
    error('More arguments are needed for Jakes_Flat()');
end
N0=8;                 % As suggested by Jakes
N=4*N0+2;             % an accurate approximation
wd=2*pi*fd;           % Maximum doppler frequency[rad]
%t_state = t0;
%for i=1:Ns
%   ich=sqrt(2)*cos(phi_N)*cos(wd*t_state);
%   qch=sqrt(2)*sin(phi_N)*cos(wd*t_state);
%   for k=1:N0
%      phi_n=pi*k/(N0+1);
%      wn=wd*cos(2*pi*k/N);
%      ich=ich+2*cos(phi_n)*cos(wn*t_state);
%      qch=qch+2*sin(phi_n)*cos(wn*t_state);
%   end
%   h1(i) = E0/sqrt(2*N0+1)*complex(ich,qch);
%   t_state=t_state+Ts;             % save last time
%end
t = t0+[0:Ns-1]*Ts;
tf = t(end)+Ts;
coswt = [sqrt(2)*cos(wd*t); 2*cos(wd*cos(2*pi/N*[1:N0]')*t)]; % (2.32)
h = E0/sqrt(2*N0+1)*exp(j*[phi_N pi/(N0+1)*[1:N0]])*coswt;
% discrepancy = norm(h-h1)

author: LyricYang


text:
function [YM,Y]=multigauss(x,mi,sigm,c)
% [YM,Y]=multigauss(x,mu,sigm,c)
%
% computes multigaussian likelihood
%
% x   : data (columnwise vectors)
% sigm: variances vector  (diagonal of the covariance matrix)
% mu  : means
% c   :the weights
DEBUG=0;

[L,T]=size(x);

if DEBUG L,T,end

M=size(c,1);

if DEBUG M,end

% repeating, changing dimensions:
X=permute(repmat(x',[1,1,M]),[1,3,2]);      % (T,L) -> (T,M,L) one per mixture

Sigm=permute(repmat(sigm,[1,1,T]),[3,2,1]); % (L,M) -> (T,M,L)

Mu=permute(repmat(mi,[1,1,T]),[3,2,1]);     % (L,M) -> (T,M,L)

if DEBUG size(X),size(Mu),size(Sigm),pause;end

%Y=squeeze(exp( 0.5.*dot(X-Mu,(X-Mu)./Sigm))) % L dissapears: (L,T,M) -> (T,M)
lY=-0.5.*dot(X-Mu,(X-Mu)./Sigm,3);
% c,const -> (T,M) and then multiply by old Y
coef=(2.*pi).^(L./2).*sqrt(prod(sigm,1)); % c,const -> (T,M)
lcoef=repmat(log(c')-log(coef),[T,1]);

if DEBUG log(coef),lcoef,lY,pause;end

YM=exp(lcoef+lY);            % ( T,M ) one mixture per column
Y=sum(YM,2);                 % add mixtures

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example code for estimating the HRF using the Inverse-Logit Model, a
% Finte Impluse Response Model and the Canonical HRF with 2 derivatives.
% Also the code illustrates our code for detecting model misspecification.
%
% By Martin Lindquist and Tor Wager
% Created  10/02/09
% Last edited 05/20/13
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Load time course
%

mypath = which('ilogit');
if isempty(mypath), error('Cannot find directory with ilogit.m and other functions. Not on path?'); end
[mydir] = fileparts(mypath)

load(fullfile(mydir,'timecourse'))

tc = (tc- mean(tc))/std(tc);
len = length(tc);

%% Or: create your own
[xBF] = spm_get_bf(struct('dt', .5, 'name', 'hrf (with time and dispersion derivatives)', 'length', 32));
clear Xtrue
for i = 1:1, xx = conv(xBF.bf(:,i), [1 1 1 1 1 1 ]');
    Xtrue(:, i) = xx(1:66);
end
for i = 2:3, xx = conv(xBF.bf(:,i), [1 1]');
    Xtrue(:, i) = xx(1:66);
end
hrf = Xtrue * [1 .3 .2]';
xsecs = 0:.5:32;

hrf = [ 0; 0; hrf];
hrf = hrf(1:length(xsecs));
hrf = hrf ./ max(hrf);
figure; plot(xsecs, hrf, 'k')
%hrf = hrf(1:4:end); % downsample to TR, if TR is > 0.5

R = randperm(640); R = sort(R(1:36));
Run = zeros(640,1);
for i=1:length(R), Run(R(i)) = 1; end;
true_sig = conv(Run, hrf);
true_sig = true_sig(1:640);

tc_noise = noise_arp(640, [.3 0]);
tc = true_sig +  0.5 * tc_noise;
% tc = true_sig;
%figure; plot(tc);

Runc{1} = Run;

%%

create_figure; subplot(3,1,1); han = plot(tc);
title('Sample time course'); drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Settings
%

TR = 0.5;
%T = round(30/TR);
T = 30;
t = 1:TR:T;                        % samples at which to get Logit HRF Estimate
FWHM = 4;                       % FWHM for residual scan
pval = 0.01;
df = 600;
alpha = 0.001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create stick function (sample event onsets)
% Variable R contains onset times
% Variable Run contains stick (a.k.a. delta or indicator) function

% R = [3, 21, 56, 65, 109, 126, 163, 171, 216, 232, 269, 282, 323, 341, 376, 385, 429, 446, 483, 491, 536, 552, 589, 602];
% Run = zeros(640,1);
% for i=1:length(R), Run(R(i)) = 1; end;
%

try
    hold on;
    hh = plot_onsets(R,'k',-3,1, 1);
    drawnow
catch
    disp('Couldn''t find function to add onset sticks to plot. Skipping.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using IL-function

% Choose mode (deterministic/stochastic)

mode = 0;   % 0 - deterministic aproach
            % 1 - simulated annealing approach
            % Please note that when using simulated annealing approach you
            % may need to perform some tuning before use.

[h1, fit1, e1, param] = Fit_Logit2(tc,TR,Runc,T,mode);
[pv sres sres_ns1] = ResidScan(e1, FWHM);
[PowLoss1] = PowerLoss(e1, fit1, (len-7) , tc, TR, Runc, alpha);

hold on; han(2) = plot(fit1,'r');

disp('Summary: IL_function');

disp('Amplitude:'); disp(param(1));
disp('Time-to-peak:'); disp(param(2)*TR);
disp('Width:'); disp(param(3)*TR);

disp('MSE:'); disp((1/(len-1)*sum(e1.^2)));
disp('Mis-modeling:'); disp(pv);
disp('Power Loss:'); disp(PowLoss1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using FIR-model

% Choose mode (FIR/sFIR)

mode = 1;   % 0 - FIR
            % 1 - smooth FIR

[h2, fit2, e2, param] = Fit_sFIR(tc,TR,Runc,T,mode);
[pv sres sres_ns2] = ResidScan(e2, FWHM);
[PowLoss2] = PowerLoss(e2, fit2, (len-T) , tc, TR, Runc, alpha);

hold on; han(3) = plot(fit2,'g');

disp('Summary: FIR');

disp('Amplitude'); disp(param(1));
disp('Time-to-peak'); disp(param(2)*TR);
disp('Width'); disp(param(3)*TR);

disp('MSE:'); disp((1/(len-1)*sum(e2.^2)));
disp('Mis-modeling'); disp(pv);
disp('Power Loss:'); disp(PowLoss2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using Canonical HRF + 2 derivatives

p=1;

[h3, fit3, e3, param, info] = Fit_Canonical_HRF(tc,TR,Runc,30,p);
[pv sres sres_ns3] = ResidScan(e3, FWHM);
[PowLoss3] = PowerLoss(e3, fit3, (len-p) , tc, TR, Runc, alpha);

hold on; han(4) = plot(fit3,'m');

legend(han,{'Data' 'IL' 'sFIR' 'DD'})

disp('Summary: Canonical + 2 derivatives');

disp('Amplitude'); disp(param(1));
disp('Time-to-peak'); disp(param(2)*TR);
disp('Width'); disp(param(3)*TR);

disp('MSE:'); disp((1/(len-1)*sum(e3.^2)));
disp('Mis-modeling'); disp(pv);
disp('Power Loss:'); disp(PowLoss3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure;
%%

subplot(3,2,5); hold on;
plot(xsecs, hrf, 'k')
xsecs1 = xsecs(1:length(h1));
han2 = plot(xsecs1, h1,'r');
xsecs2 = xsecs(1:length(h2));
han2(2) = plot(xsecs2, h2,'g');
xsecs3 = xsecs(1:length(h3));
han2(3) = plot(xsecs3, h3,'m');
legend(han2,{'IL' 'sFIR' 'DD'})
title('Estimated HRF');

subplot(3,1,2); hold on;
hh = plot_onsets(R,'k',-3,1);
drawnow

han3 = plot(sres_ns1,'r');
hold on; han3(2) = plot(sres_ns2,'g');
hold on; han3(3) = plot(sres_ns3,'m');
hold on; plot((1:len),zeros(len,1),'--k');
legend(han3,{'IL' 'sFIR' 'DD'})
title('Mis-modeling (time course)');

subplot(3,2,6); hold on;

[s1] = Fit_sFIR(sres_ns1,TR,Runc,T,0);
[s2] = Fit_sFIR(sres_ns2,TR,Runc,T,0);
[s3] = Fit_sFIR(sres_ns3,TR,Runc,T,0);

han4 = plot(s1(1:T),'r');
hold on; han4(2) = plot(s2(1:T),'g');
hold on; han4(3) = plot(s3(1:T),'m');
hold on; plot((1:T),zeros(T,1),'--k');
legend(han4,{'IL' 'sFIR' 'DD'})
title('Mis-modeling (HRF)');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function atanh_test ( )

%*****************************************************************************80
%
%% ATANH_TEST tests R4_ATANH and R8_ATANH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ATANH_TEST:\n' );
  fprintf ( 1, '  Test ARCTANH_VALUES, R4_ATANH, R8_ATANH\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X     ARCTANH(X)\n' );
  fprintf ( 1, '                  R4_ATANH(X)         Diff\n' );
  fprintf ( 1, '                  R8_ATANH(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

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

    if ( n_data == 0 )
      break
    end

    fx2 = r4_atanh ( single ( x ) );
    fx3 = r8_atanh ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14f  %14g\n', x, fx1 );
    fprintf ( 1, '                  %14g  %14g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14g  %14g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer


text:
function [u,c]=getEllipsHarmAxes(pointHarmon,E)
%%GETELLIPSHARMAXES Find the Cartesian unit vectors corresponding to the
%                   normalized  gradients of the components of a point in
%                   ellipsoidal harmonic coordinates. Also, get the
%                   magnitudes of the gradients.
%
%INPUTS: pointHarmon A point given in terms of ellipsoidal harmonic reduced
%                    latitude and longitude in radians and a semi-major
%                    axis in meters.
%                  E The linear eccentricity defining the ellipsoidal
%                    harmonic coordinate system. If this parameter is
%                    omitted, then the linear eccentricity of the WGS84
%                    reference ellipsoid is used.
%
%OUTPUTS: u u(:,1), u(:,2) and u(:,3) are the unit vectors in the
%           respective directions of the respective gradients of reduced
%           latitude, longitude and the semi-major axis.
%         c c(1), c(2) and c(3) are the respective magnitudes of the
%           derivative of the Cartesian position with respect to reduced
%           latitude, longitude and the semi-major axis.
%
%The ellipsoidal harmonic coordinate system is described in Chapter 1.15
%of [1] and is an orthogonal coordinate system.
%
%REFERENCES:
%[1] B. Hofmann-Wellenhof and H. Moritz, Physical Geodesy, 2nd ed.
%    SpringerWienNewYork, 2006.
%
%January 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(E))
    a=Constants.WGS84SemiMajorAxis;
    f=Constants.WGS84Flattening;
    b=a*(1 - f);
    E=sqrt(a^2-b^2);
end

    beta=pointHarmon(1);
    lambda=pointHarmon(2);
    u=pointHarmon(3);

dBeta=[-sqrt(E^2+u^2)*sin(beta)*cos(lambda);
       -sqrt(E^2+u^2)*sin(beta)*sin(lambda);
        u*cos(beta)];
c1=norm(dBeta);
dBeta=dBeta/c1;

du=[u*cos(beta)*cos(lambda)/sqrt(E^2+u^2);
    u*cos(beta)*sin(lambda)/sqrt(E^2+u^2);
    sin(beta)];
c3=norm(du);
du=du/c3;

dLambda=[-sqrt(E^2+u^2)*cos(beta)*sin(lambda);
          sqrt(E^2+u^2)*cos(beta)*cos(lambda);
          0];
c2=norm(dLambda);

%If the point is too close to the poles, then it is possible that c2 is
%nearly equal to zero. However, a normalized dLambda can just be found by
%orthogonality: it is orthogonal to dBeta and du.
dLambda=cross(dBeta,du);

u=[dBeta,dLambda,du];
c=[c1;c2;c3];
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function sphere_triangle_quad_test05 ( )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_TEST05 tests SPHERE01_TRIANGLE_QUAD_ICOS1V.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    2, 2, 2; ...
    0, 2, 4; ...
    0, 0, 6; ...
    1, 2, 4; ...
    2, 4, 2; ...
    6, 2, 0; ...
    0, 0, 8; ...
    6, 0, 4; ...
    4, 6, 2; ...
    2, 4, 8; ...
   16, 0, 0 ]';
  n_mc1 = 1000;
  n_mc2 = 10000;
  n_mc3 = 100000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_TRIANGLE_QUAD_TEST05\n' );
  fprintf ( 1, '  SPHERE01_TRIANGLE_QUAD_ICOS1V approximates the\n' );
  fprintf ( 1, '  integral of a function over a spherical triangle on\n' );
  fprintf ( 1, '  the surface of the unit sphere using a vertex rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We do not have an exact result, so we compare each\n' );
  fprintf ( 1, '  estimate to the final one.\n' );
%
%  Choose three points at random to define a spherical triangle.
%
  [ v1, seed ] = sphere01_sample ( 1, seed );
  [ v2, seed ] = sphere01_sample ( 1, seed );
  [ v3, seed ] = sphere01_sample ( 1, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vertices of random spherical triangle:\n' );
  fprintf ( 1, '\n' );
  r8vec_transpose_print ( 3, v1, '  V1:' );
  r8vec_transpose_print ( 3, v2, '  V2:' );
  r8vec_transpose_print ( 3, v3, '  V3:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  FACTOR   N   RESULT\n' );

  for j = 1 : 17

    e(1:3) = e_test(1:3,j);

    e = polyterm_exponent ( 'SET', e );

    polyterm_exponent ( 'PRINT', e );

%   factor = 2 ^ 11;
    factor = 2 ^ 7;
    [ best, node_num ] = sphere01_triangle_quad_icos1v ( v1, v2, v3, factor, ...
      @polyterm_value_3d );

    factor = 1;
    for factor_log = 0 : 5

      [ result, node_num ] = sphere01_triangle_quad_icos1v ( v1, v2, v3, ...
        factor, @polyterm_value_3d );

      error = abs ( result - best );

      fprintf ( 1, '  %4d  %8d  %16.8g  %10.2e\n', ...
        factor, node_num, result, error );

      factor = factor * 2;

    end

  end

  return
end