text:
function [mu,dmu,k,gamma] = sparse_learning(Phi,T,lambda,iters,flag1,flag2,flag3)
% *************************************************************************
%
% *** PURPOSE ***
% Implements generalized versions of SBL and FOCUSS for learning sparse
% representations from possibly overcomplete dictionaries.
%
%
% *** USAGE ***
% [mu,dmu,k,gamma] = sparse_learning(Phi,T,lambda,iters,flag1,flag2,flag3);
%
%
% *** INPUTS ***
% Phi       = N X M dictionary
% T         = N X L data matrix
% lambda    = scalar trade-off parameter (balances sparsity and data fit)
% iters     = maximum number of iterations
%
% flag1     = 0: fast Mackay-based SBL update rules
% flag1     = 1: fast EM-based SBL update rule
% flag1     = 2: traditional (slow but sometimes better) EM-based SBL update rule
% flag1     = [3 p]: FOCUSS algorithm using the p-valued quasi-norm
%
% flag2     = 0: regular initialization (equivalent to min. norm solution)
% flag2     = gamma0: initialize with gamma = gamma0, (M X 1) vector
%
% flag3     = display flag; 1 = show output, 0 = supress output
%
% *** OUTPUTS ***
% mu        = M X L matrix of weight estimates
% dmu       = delta-mu at convergence
% k         = number of iterations used
% gamma     = M X 1 vector of hyperparameter values
%
%
% *************************************************************************
% Written by:  David Wipf, david.wipf@mrsc.ucsf.edu
% *************************************************************************

% *** Control parameters ***
MIN_GAMMA       = 1e-16;
MIN_DMU         = 1e-12;
MAX_ITERS       = iters;
DISPLAY_FLAG    = flag3;     % Set to zero for no runtime screen printouts

% *** Initializations ***
[N M] = size(Phi);
[N L] = size(T);

if (~flag2)         gamma = ones(M,1);
else                gamma = flag2;  end;

keep_list = [1:M]';
m = length(keep_list);
mu = zeros(M,L);
dmu = -1;
k = 0;

% *** Learning loop ***
for k=0:MAX_ITERS-1

    % *** Prune things as hyperparameters go to zero ***
    if (min(gamma) < MIN_GAMMA )
		index = find(gamma > MIN_GAMMA);
		gamma = gamma(index);
		Phi = Phi(:,index);
		keep_list = keep_list(index);
        m = length(gamma);

        if (m == 0)   break;  end;
    end;

    % *** Compute new weights ***
    G = repmat(sqrt(gamma)',N,1);
    PhiG = Phi.*G;
    [U,S,V] = svd(PhiG,'econ');

    [d1,d2] = size(S);
    if (d1 > 1)     diag_S = diag(S);
    else            diag_S = S(1);      end;

    U_scaled = U(:,1:min(N,m)).*repmat((diag_S./(diag_S.^2 + lambda + 1e-16))',N,1);
    Xi = G'.*(V*U_scaled');

    mu = Xi*T;

    % *** Update hyperparameters ***
    gamma_old = gamma;
    mu2_bar = sum(mu.^2,2);

    if (flag1(1) == 0)
        % MacKay fixed-point SBL
        R_diag = real( (sum(Xi.'.*Phi)).' );
        gamma = mu2_bar./(L*R_diag);

    elseif (flag1(1) == 1)
        % Fast EM SBL
        R_diag = real( (sum(Xi.'.*Phi)).' );
        gamma = sqrt( gamma.*real(mu2_bar./(L*R_diag)) );

    elseif (flag1(1) == 2)
        % Traditional EM SBL
        PhiGsqr = PhiG.*G;
        Sigma_w_diag = real( gamma - ( sum(Xi.'.*PhiGsqr) ).' );
        gamma = mu2_bar/L + Sigma_w_diag;

    else
        % FOCUSS
        p = flag1(2);
        gamma = (mu2_bar/L).^(1-p/2);
    end;

    % *** Check stopping conditions, etc. ***
    if (DISPLAY_FLAG) disp(['iters: ',num2str(k),'   num coeffs: ',num2str(m), ...
            '   gamma change: ',num2str(max(abs(gamma - gamma_old)))]); end;
end;

% *** Expand weights, hyperparameters ***
temp = zeros(M,1);
if (m > 0) temp(keep_list,1) = gamma;  end;
gamma = temp;

temp = zeros(M,L);
if (m > 0) temp(keep_list,:) = mu;  end;
mu = temp;

return;

repo: REST


text:
%--- help for dsge/check_derivatives ---
%
%  Compares the derivatives and the solutions from various differentiation techniques
%
%  ::
%
%    check_derivatives(obj)
%    retcode=check_derivatives(obj)
%
%  Args:
%     obj (rise | dsge): model object or vectors of model objects
%
%  Returns:
%     :
%
%     - **retcode** [numeric]: 0 if no problem is encountered during the
%       comparisons. Else the meaning of recode can be found by running
%       decipher(retcode)
%
%  Note:
%
%     - The derivatives computed are 'automatic', 'symbolic' or 'numerical'
%     - The comparisons are done relative to automatic derivatives, which are
%       assumed to be the most accurate.
%
%

author: jmaih


text:
function obj = updateStem3(obj,dataIndex)

  %-------------------------------------------------------------------------%

  %-AXIS INDEX-%
  axIndex = obj.getAxisIndex(obj.State.Plot(dataIndex).AssociatedAxis);

  %-PLOT DATA STRUCTURE- %
  data = get(obj.State.Plot(dataIndex).Handle);

  %-CHECK FOR MULTIPLE AXES-%
  [xsource, ysource] = findSourceAxis(obj,axIndex);

  %------------------------------------------------------------------------%

  %-get coordenate x,y,z data-%
  xdata = data.XData;
  ydata = data.YData;
  zdata = data.ZData;
  npoints = length(xdata);

  %------------------------------------------------------------------------%

  %-check if stem3-%
  isstem3 = ~isempty(zdata);

  %------------------------------------------------------------------------%

  %-SCENE-%
  if isstem3
    eval(['scene = obj.layout.scene' num2str(xsource) ';']);
  else
    eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
    eval(['yaxis = obj.layout.yaxis' num2str(xsource) ';']);
  end

  %-------------------------------------------------------------------------%

  %-scatter3d scene-%
  if isstem3
    stem_data.scene = sprintf('scene%d', xsource);
  else
    stem_data.xaxis = sprintf('x%d', xsource);
    stem_data.yaxis = sprintf('y%d', xsource);
  end

  %------------------------------------------------------------------------%

  %-scatter3d type-%
  if isstem3
    stem_data.type = 'scatter3d';
  else
    stem_data.type = 'scatter';
  end

  %-------------------------------------------------------------------------%

  %-scatter3d visible-%
  stem_data.visible = strcmp(data.Visible,'on');

  %-------------------------------------------------------------------------%

  %-scatter3d name-%
  stem_data.name = data.DisplayName;

  %------------------------------------------------------------------------%

  %-scatter mode-%
  stem_data.mode = 'lines+markers';

  %------------------------------------------------------------------------%

  %-allocated space for extended data-%
  xdata_extended = zeros(3*npoints, 1);
  ydata_extended = zeros(3*npoints, 1);

  if isstem3
    zdata_extended = zeros(3*npoints, 1);
  end

  %-format data-%
  m = 1;
  for n = 1:npoints
    %-x data-%
    xdata_extended(m) = xdata(n);
    xdata_extended(m+1) = xdata(n);
    xdata_extended(m+2) = nan;

    %-y data-%
    ydata_extended(m) = 0;
    if isstem3
      ydata_extended(m) = ydata(n);
    end
    ydata_extended(m+1) = ydata(n);
    ydata_extended(m+2) = nan;

    %-z data-%
    if isstem3
      zdata_extended(m) = 0;
      zdata_extended(m+1) = zdata(n);
      zdata_extended(m+2) = nan;
    end

    m = m + 3;
  end

  %------------------------------------------------------------------------%

  %-scatter3d line-%
  stem_data.line = extractLineLine(data);

  %-------------------------------------------------------------------------%

  %-scatter3d marker-%
  stem_data.marker = extractLineMarker(data);

  %-------------------------------------------------------------------------%

  if isstem3

    %-fix marker symbol-%
    symbol = stem_data.marker.symbol;

    if strcmpi(symbol, 'asterisk-open') || strcmpi(symbol, 'cross-thin-open')
      stem_data.marker.symbol = 'cross';
    end

    stem_data.marker.size = stem_data.marker.size * 0.6;

    %-fix dash line-%
    dash = stem_data.line.dash;

    if strcmpi(dash, 'dash')
      stem_data.line.dash = 'dot';
    end
  end

  %-------------------------------------------------------------------------%

  %-hide every other marker-%
  markercolor = cell(3*npoints,1);
  linecolor = cell(3*npoints,1);
  hidecolor = 'rgba(0,0,0,0)';

  linecolor(1:3:3*npoints) = {hidecolor};
  markercolor(1:3:3*npoints) = {hidecolor};

  try
    linecolor(2:3:3*npoints) = {stem_data.marker.line.color};
    markercolor(2:3:3*npoints) = {stem_data.marker.color};
  catch
    linecolor(2:3:3*npoints) = {stem_data.marker.color};
    markercolor(2:3:3*npoints) = {hidecolor};
  end

  linecolor(3:3:3*npoints) = {hidecolor};
  markercolor(3:3:3*npoints) = {hidecolor};

  %-add new marker/line colors-%
  stem_data.marker.color = markercolor;
  stem_data.marker.line.color = linecolor;

  stem_data.marker.line.width = stem_data.marker.line.width * 2;
  stem_data.line.width = stem_data.line.width * 2;

  %------------------------------------------------------------------------%

  %-set x y z data-%
  stem_data.x = xdata_extended;
  stem_data.y = ydata_extended;

  if isstem3
    stem_data.z = zdata_extended;
  end

  %------------------------------------------------------------------------%

  %-set plotly data-%
  obj.data{dataIndex} = stem_data;

  %------------------------------------------------------------------------%

  %-SETTING SCENE-%

  if isstem3
    %-aspect ratio-%
    asr = obj.PlotOptions.AspectRatio;

    if ~isempty(asr)
      if ischar(asr)
        scene.aspectmode = asr;
      elseif isvector(asr) && length(asr) == 3
        xar = asr(1);
        yar = asr(2);
        zar = asr(3);
      end
    else

        %-define as default-%
        xar = max(xdata(:));
        yar = max(ydata(:));
        xyar = max([xar, yar]);
        xar = xyar; yar = xyar;
        zar = 0.7*max([xar, yar]);
    end

    scene.aspectratio.x = xar;
    scene.aspectratio.y = yar;
    scene.aspectratio.z = zar;

    %---------------------------------------------------------------------%

    %-camera eye-%
    ey = obj.PlotOptions.CameraEye;

    if ~isempty(ey)
        if isvector(ey) && length(ey) == 3
            scene.camera.eye.x = ey(1);
            scene.camera.eye.y = ey(2);
            scene.camera.eye.z = ey(3);
        end
    else

        %-define as default-%
        xey = - xar; if xey>0 xfac = 0.2; else xfac = -0.2; end
        yey = - yar; if yey>0 yfac = -0.2; else yfac = 0.2; end
        if zar>0 zfac = 0.2; else zfac = -0.2; end

        scene.camera.eye.x = xey + xfac*xey;
        scene.camera.eye.y = yey + yfac*yey;
        scene.camera.eye.z = zar + zfac*zar;
    end

    %---------------------------------------------------------------------%

    %-zerolines hidded-%
    scene.xaxis.zeroline = false;
    scene.yaxis.zeroline = false;
    scene.zaxis.zeroline = false;

    scene.xaxis.linecolor = 'rgba(0,0,0,0.8)';
    scene.yaxis.linecolor = 'rgba(0,0,0,0.8)';
    scene.zaxis.linecolor = 'rgba(0,0,0,0.5)';

    scene.xaxis.ticklen = 5;
    scene.yaxis.ticklen = 5;
    scene.zaxis.ticklen = 5;

    scene.xaxis.tickcolor = 'rgba(0,0,0,0.8)';
    scene.yaxis.tickcolor = 'rgba(0,0,0,0.8)';
    scene.zaxis.tickcolor = 'rgba(0,0,0,0.8)';

    scene.xaxis.range = data.Parent.XLim;
    scene.yaxis.range = data.Parent.YLim;
    scene.zaxis.range = data.Parent.ZLim;

    scene.xaxis.tickvals = data.Parent.XTick;
    scene.yaxis.tickvals = data.Parent.YTick;
    scene.zaxis.tickvals = data.Parent.ZTick;

    scene.xaxis.title = data.Parent.XLabel.String;
    scene.yaxis.title = data.Parent.YLabel.String;
    scene.zaxis.title = data.Parent.ZLabel.String;

    %-update scene-%
    obj.layout = setfield(obj.layout, sprintf('scene%d', xsource), scene);

  else
    yaxis.zeroline = true;

    xaxis.linecolor = 'rgba(0,0,0,0.4)';
    yaxis.linecolor = 'rgba(0,0,0,0.4)';

    xaxis.tickcolor = 'rgba(0,0,0,0.4)';
    yaxis.tickcolor = 'rgba(0,0,0,0.4)';

    xaxis.tickvals = data.Parent.XTick;
    yaxis.tickvals = data.Parent.YTick;

    xaxis.title = data.Parent.XLabel.String;
    yaxis.title = data.Parent.YLabel.String;

    %-update axis-%
    obj.layout = setfield(obj.layout, sprintf('xaxis%d', xsource), xaxis);
    obj.layout = setfield(obj.layout, sprintf('yaxis%d', ysource), yaxis);
  end

  %------------------------------------------------------------------------%

end

author: plotly


text:
  function img = easyhelix(cg, ig, proj, varargin)
%|
%| FBP reconstruction of cone-beam tomography data collected with
%| a circular source trajectory.
%| See feldkamp_example.m for example.
%|
%| in
%|	cg			ct_geom()
%|	ig			image_geom()
%|	proj	[ns nt na]	cone-beam projection views (line integrals)
%|
%|
%| out
%|	img	[nx ny nz]	reconstructed image
%|

% defaults
img = helix_do(proj, cg, ig, cg.na, cg.dt, cg.dsd, cg.dso, cg.orbit, cg.orbit_start);
end % feldkamp()

%
% feldkamp_do()
%
function img = helix_do(proj, cg, ig, na, dt, dsd, dso, orbit, orbit_start)

% step 1: fix z-sampling; for each CB source point, determine multifan
% betas = deg2rad(orbit_start + orbit * [0:na-1] / na); % [na] source angles

% assumes that na is even
num_turns = orbit/360;
phis = mod(deg2rad(orbit_start + orbit * [0:(na/num_turns-1)]/na), 2*pi);
numPhis = size(phis',1);

if ig.dz < .9*(dso/dsd * dt)
    min_spacing = ig.dz;
else
    min_spacing = .9*(dso/dsd * dt);
end

num_zsamp = ig.nz;
% num_zsamp = ceil((ig.z(size(ig.z,1))-ig.z(1))/min_spacing);
% zSamples = ig.z(1) + min_spacing * [0:num_zsamp-1];
zSamples = ig.z;
% -cg.source_zs(1)+cg.source_zs(2)
% myPitch = (-cg.source_zs(1)+cg.source_zs(2))*(na/2-1)
myPitch = cg.pitch * cg.nt * dso / dsd * cg.dt;

delta = asin(cg.rmax/dso);
dist = .5 * myPitch * (pi + 2 * delta)/(2*pi);

%step 2: rebin the cone beam data to fit a fan beam projection
[ns, nt, na] = size(proj);
fanBeam = zeros(ns, numPhis, num_zsamp);
for i=0:na-1
    currentZ = cg.source_zs(i+1);
    upper = currentZ + dist;
    lower = currentZ - dist;

    lcount = 1;

    while lcount < num_zsamp+1 && zSamples(lcount) < lower
        lcount = lcount + 1;
    end

    while lcount < num_zsamp+1 && zSamples(lcount) < upper
        deltaZ = zSamples(lcount)-currentZ;
        shortscan = (pi/2+delta)*(1-deltaZ/dist);

        is = ndgrid(1:cg.ns,1);
        spoints = cg.s;

        % tpoints = ((spoints).^2+dsd*dsd)./ (dso*dsd).* deltaZ;
        tpoints = ones(cg.ns,1).* deltaZ;

        scale = sqrt(spoints.^2+dsd^2)./sqrt(spoints.^2+tpoints.^2+dsd^2);

        t_in = (tpoints./cg.dt + cg.wt);

        x0 = floor(t_in);
        x1 = 1 + x0;

        alpha = t_in - x0;

        for j = 1:size(x0,1)
            weight = myParker(spoints(j),shortscan,delta,dsd);
            fanBeam(is(j), mod(i,numPhis)+1, lcount) = weight*proj(is(j), nt/2, i+1);
        end

        lcount = lcount + 1;
    end

end

% Everything below here goes really fast
disp('we are here')
    orbitStart = zeros(num_zsamp,1);
    % decide on beginning orbit
    orbit = deg2rad(cg.orbit_start);;
    zloc = cg.source_z0;
    change_phis = phis(2)-phis(1);
    for z = 1: num_zsamp
        orbitFound = 0;
        while orbitFound == 0
            if zSamples(z) < zloc + dist
                orbitStart(z) = orbit;
                orbitFound = 1;
            else
                orbit = orbit + change_phis;
                zloc = zloc + (cg.source_zs(2) - cg.source_zs(1));
            end
        end
    end

    newNA = floor(2*dist/myPitch*numPhis);
    newFanBeam = zeros(ns, newNA, num_zsamp);

    for z = 1 : num_zsamp
        angle = round((orbitStart(z)-deg2rad(cg.orbit_start)) / change_phis);
        for o = 0: newNA-1
            nextAngle = angle + o;
            newFanBeam(:, o+1, z) = fanBeam(:, mod(nextAngle, numPhis)+1, z);
        end
    end

    down = cg.down;
    wantedZ=ceil(ig.nz/2);
     im(newFanBeam(:,:,wantedZ), 'test'), cbar

    for z = 1:num_zsamp
        sinost = sino_geom('fan', 'ns', cg.ns*down, 'na', newNA*down, ...
            'ds', cg.ds/down, 'down', down, 'orbit', 'short', ...
            'orbit_start', rad2deg(orbitStart(z)), ...
            'offset_s', cg.offset_s, ...
            'dsd', cg.dsd, 'dod', cg.dod, 'dfs', cg.dfs);
        wt = fbp_fan_short_wt(sinost);

        for i = 1:ns
            for j = 1:newNA
                newFanBeam(i, j, z) = wt(i,j)*newFanBeam(i, j, z);
            end
        end
    end

    im(newFanBeam(:,:,:), 'test'), cbar

   % same image geometry as before, except for the z direction
    ig = image_geom('nx', ig.nx*ig.down, 'ny', ig.ny*ig.down, 'dx', ...
        ig.dx/ig.down, 'down', ig.down);
    mask2 = true([ig.nx ig.ny]);
    mask2(end) = 0; % trick: test it
    ig.mask = repmat(mask2, [1 1]);
    clear mask2

    img = zeros(ig.nx,ig.ny,num_zsamp);

    for i = 1:num_zsamp
        sinost = sino_geom('fan', 'ns', cg.ns*down, 'na', newNA*down, ...
            'ds', cg.ds/(down), 'down', down, 'orbit', 'short', ...
            'orbit_start', rad2deg(orbitStart(i)), ...
            'offset_s', cg.offset_s, ...
            'dsd', cg.dsd, 'dod', cg.dod, 'dfs', cg.dfs);
        geomsino = fbp2(sinost, ig);
        img(:,:,i) = fbp2(newFanBeam(:,:,i), geomsino);
    end

    im(img(:,:,round(num_zsamp/2)+1)), cbar;

end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function cae = caenumgradcheck(cae, x, y)
    epsilon = 1e-4;
    er = 1e-6;
    disp('performing numerical gradient checking...')
    for i = 1 : numel(cae.o)
        p_cae = cae; p_cae.c{i} = p_cae.c{i} + epsilon;
        m_cae = cae; m_cae.c{i} = m_cae.c{i} - epsilon;

        [m_cae, p_cae] = caerun(m_cae, p_cae, x, y);
        d = (p_cae.L - m_cae.L) / (2 * epsilon);

        e = abs(d - cae.dc{i});
        if e > er
            disp('OUTPUT BIAS numerical gradient checking failed');
            disp(e);
            disp(d / cae.dc{i});
            keyboard
        end
    end

    for a = 1 : numel(cae.a)

        p_cae = cae; p_cae.b{a} = p_cae.b{a} + epsilon;
        m_cae = cae; m_cae.b{a} = m_cae.b{a} - epsilon;

        [m_cae, p_cae] = caerun(m_cae, p_cae, x, y);
        d = (p_cae.L - m_cae.L) / (2 * epsilon);
%        cae.dok{i}{a}(u) = d;
        e = abs(d - cae.db{a});
        if e > er
            disp('BIAS numerical gradient checking failed');
            disp(e);
            disp(d / cae.db{a});
            keyboard
        end

        for i = 1 : numel(cae.o)
            for u = 1 : numel(cae.ok{i}{a})
                p_cae = cae; p_cae.ok{i}{a}(u) = p_cae.ok{i}{a}(u) + epsilon;
                m_cae = cae; m_cae.ok{i}{a}(u) = m_cae.ok{i}{a}(u) - epsilon;

                [m_cae, p_cae] = caerun(m_cae, p_cae, x, y);
                d = (p_cae.L - m_cae.L) / (2 * epsilon);
%                cae.dok{i}{a}(u) = d;
                e = abs(d - cae.dok{i}{a}(u));
                if e > er
                    disp('OUTPUT KERNEL numerical gradient checking failed');
                    disp(e);
                    disp(d / cae.dok{i}{a}(u));
%                    keyboard
                end
            end
        end

        for i = 1 : numel(cae.i)
            for u = 1 : numel(cae.ik{i}{a})
                p_cae = cae;
                m_cae = cae;
                p_cae.ik{i}{a}(u) = p_cae.ik{i}{a}(u) + epsilon;
                m_cae.ik{i}{a}(u) = m_cae.ik{i}{a}(u) - epsilon;
                [m_cae, p_cae] = caerun(m_cae, p_cae, x, y);
                d = (p_cae.L - m_cae.L) / (2 * epsilon);
%                cae.dik{i}{a}(u) = d;
                e = abs(d - cae.dik{i}{a}(u));
                if e > er
                    disp('INPUT KERNEL numerical gradient checking failed');
                    disp(e);
                    disp(d / cae.dik{i}{a}(u));
                end
            end
        end
    end

    disp('done')

end

function [m_cae, p_cae] = caerun(m_cae, p_cae, x, y)
    m_cae = caeup(m_cae, x); m_cae = caedown(m_cae); m_cae = caebp(m_cae, y);
    p_cae = caeup(p_cae, x); p_cae = caedown(p_cae); p_cae = caebp(p_cae, y);
end

%function checknumgrad(cae,what,x,y)
%    epsilon = 1e-4;
%    er = 1e-9;
%
%    for i = 1 : numel(eval(what))
%        if iscell(eval(['cae.' what]))
%            checknumgrad(cae,[what '{' num2str(i) '}'], x, y)
%        else
%            p_cae = cae;
%            m_cae = cae;
%            eval(['p_cae.' what '(' num2str(i) ')']) = eval([what '(' num2str(i) ')']) + epsilon;
%            eval(['m_cae.' what '(' num2str(i) ')']) = eval([what '(' num2str(i) ')']) - epsilon;
%
%            m_cae = caeff(m_cae, x); m_cae = caedown(m_cae); m_cae = caebp(m_cae, y);
%            p_cae = caeff(p_cae, x); p_cae = caedown(p_cae); p_cae = caebp(p_cae, y);
%
%            d = (p_cae.L - m_cae.L) / (2 * epsilon);
%            e = abs(d - eval(['cae.d' what '(' num2str(i) ')']));
%            if e > er
%                error('numerical gradient checking failed');
%            end
%         end
%     end
%
% end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function [V,F]=platonic_solid(varargin)

% function [V,F]=platonic_solid(n,r)
% ------------------------------------------------------------------------
% PLATONIC_SOLID Creates the PATCH data, the vertices (V) and faces (F) for
% a given platonic solid (according to "n" see below) with radius (r)
%
% n=1 -> Tetrahedron
% n=2 -> Cube
% n=3 -> Octahedron
% n=4 -> Icosahedron
% n=5 -> Dodecahedron
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 12/08/2014 Updated for GIBBON
%------------------------------------------------------------------------
%% Parse input

switch nargin
    case 1
        n=varargin{1};
        r=1;
    case 2
        n=varargin{1};
        r=varargin{2};
end

if isempty(r)
    r=1;
end

%%
switch n
    case 1 % Tetrahedron
        X=[-0.5;0.5;0;0;];
        Y=[-sqrt(3)/6;  -sqrt(3)/6; sqrt(3)/3; 0];
        Z=[-0.25.*sqrt(2/3); -0.25.*sqrt(2/3); -0.25.*sqrt(2/3);  0.75.*sqrt(2/3)];
        X=X([1 2 4 3]);
        Y=Y([1 2 4 3]);
        Z=Z([1 2 4 3]);
        F=[3,2,1;1,2,4;2,3,4;4,3,1;];
        F=fliplr(F);
    case 2 % Cube
        X=[-1;  1; 1; -1; -1;  1; 1; -1;];
        Y=[-1; -1; 1;  1; -1; -1; 1;  1;];
        Z=[-1; -1;-1; -1;  1;  1; 1;  1;];
        F=[4,3,2,1;
            1,2,6,5;
            2,3,7,6;
            3,4,8,7;
            4,1,5,8;
            5,6,7,8;];
    case 3 % Octahedron
        X=[-1;  1; 1; -1;  0;   0;];
        Y=[-1; -1; 1;  1;  0;   0;];
        Z=[ 0;   0; 0;  0; -1;  1;];
        F=[5,2,1;5,3,2;5,4,3;5,1,4;1,2,6;2,3,6;3,4,6;4,1,6;];
    case 4 % Icosahedron
        phi=(1+sqrt(5))/2;
        X=[0;0;0;0;-1;-1;1;1;-phi;phi;phi;-phi;];
        Y=[-1;-1;1;1;-phi;phi;phi;-phi;0;0;0;0;];
        Z=[-phi;phi;phi;-phi;0;0;0;0;-1;-1;1;1;];
        F=[9,4,1;1,5,9;1,8,5;10,8,1;4,10,1;5,2,12;12,2,3;12,3,6;12,6,9;12,9,5;10,7,11;8,10,11;2,8,11;3,2,11;7,3,11;2,5,8;10,4,7;7,6,3;6,7,4;6,4,9;];
    case 5 % Dodecahedron
        phi=(1+sqrt(5))/2;
        X=[1;(1/phi);-phi;phi;-1;0;-phi;1;-1;-1;1;(1/phi);-1;0;0;-(1/phi);phi;-(1/phi);1;0;];
        Y=[1;0;-(1/phi);(1/phi);1;-phi;(1/phi);-1;1;-1;-1;0;-1;-phi;phi;0;-(1/phi);0;1;phi;];
        Z=[1;phi;0;0;-1;-(1/phi);0;1;1;1;-1;-phi;-1;(1/phi);-(1/phi);phi;0;-phi;-1;(1/phi);];
        F=[20,9,16,2,1;2,16,10,14,8;16,9,7,3,10;7,9,20,15,5;18,13,3,7,5;3,13,6,14,10;6,13,18,12,11;6,11,17,8,14;11,12,19,4,17;1,2,8,17,4;1,4,19,15,20;12,18,5,15,19;];
    otherwise
        warning('False input for n')
end

%Altering radius
[THETA,PHI,~]=cart2sph(X,Y,Z);
[X,Y,Z]=sph2cart(THETA,PHI,r.*ones(size(X)));
V=[X(:) Y(:) Z(:)];

%%
% _*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:
%                                      [zhi,zlo] = quadadd(xhi,xlo,y)
% QUADADD   Compute (zhi+zlo) = (xhi+xlo) + y.
%   x an z are in double-double format (quad precision)
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

function [zhi,zlo] = quadadd(xhi,xlo,y) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')


text:
function [outputTable] = readExternalTable(inputTableFilename,varargin)
	% Read in table, decides whether to do a single or multiple tables, if multiple tables, should have the same column names.
	% Biafra Ahanonu
	% branched from controllerAnalysis: 2014.08.01 [16:09:16]
	% inputs
		% inputTableFilename - string pointing toward file
	% outputs
		%

	[outputTable] = ciapkg.io.readExternalTable(inputTableFilename,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function test_bug1450

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY mergestruct

% ft_checkconfig used to be responsible for mergestruct, but on 19-10-2015
% that moved to ft_preable_init which is called prior to ft_checkconfig

% this test script should be explicit about ft_default
global ft_default
ft_default = [];
ft_default.field1 = 1;
ft_default.field2 = 2;
ft_default.field3 = 3;
ft_default.sub.field1 = 1;
ft_default.sub.field2 = 2;
ft_default.sub.field3 = 3;
ft_default.sub.sub.field1 = 1;
ft_default.sub.sub.field2 = 2;
ft_default.sub.sub.field3 = 3;
ft_default.sub.sub.sub.field1 = 1;
ft_default.sub.sub.sub.field2 = 2;
ft_default.sub.sub.sub.field3 = 3;
ft_default.sub.sub.sub.sub(1).field1 = 1.1;
ft_default.sub.sub.sub.sub(2).field1 = 1.2;
ft_default.sub.sub.sub.sub(3).field1 = 1.3;

cfg = [];
cfg = mergestruct(cfg, ft_default);
assert(isequal(cfg, ft_default));

cfg = [];
cfg.field1 = 1;
cfg = mergestruct(cfg, ft_default);
% field1 should remain as it is, field 2 and 3 should have been added
assert(isequal(cfg, ft_default));

cfg = [];
cfg.sub = [];
cfg = mergestruct(cfg, ft_default);
% the subfields should have been added
assert(isequal(cfg, ft_default));

cfg = [];
cfg.sub.field1 = 1;
cfg = mergestruct(cfg, ft_default);
% the sub-fields 2 and 3 should have been added, sub-field 1 should remain as it is
assert(isequal(cfg, ft_default));

cfg = [];
cfg.sub.sub.field1 = 1;
cfg = mergestruct(cfg, ft_default);
% do the same test, but now two levels deep
assert(isequal(cfg, ft_default));

cfg = [];
cfg.sub.sub.sub.field1 = 1;
cfg = mergestruct(cfg, ft_default);
% do the same test, but now three levels deep
assert(isequal(cfg, ft_default));

cfg = [];
cfg.sub.sub.sub.sub.field2 = 2;
cfg = mergestruct(cfg, ft_default);
% should be concatenated
assert(numel(cfg.sub.sub.sub.sub)==4);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef prtPreProc < prtAction
    % prtPreProc  Base class for prt pre processing objects
    %
    % PreProcessors are currently have no additional properties or methods
    %
    % This is an abstract class and cannot be instantiated.
    %
    % A prtPreProc object inherits all methods and properties from the
    % prtAction class
    %
    %   See Also: prtPreProc, prtPreProcPca, prtPreProcPls,
    %   prtPreProcHistEq, prtPreProcZeroMeanColumns, prtPreProcLda,
    %   prtPreProcZeroMeanRows, prtPreProcLogDisc, prtPreProcZmuv,
    %   prtPreProcMinMaxRows

    properties (SetAccess = protected)
        isSupervised = false;  % False
        isCrossValidateValid = true; % True
    end

    methods
        function obj = prtPreProc()
            % As an action subclass we must set the properties to reflect
            % our dataset requirements
            obj.classTrain = 'prtDataSetStandard';
            obj.classRun = 'prtDataSetStandard';
            obj.classRunRetained = true;
        end
	end
	methods (Hidden = true)
        function featureNameModificationFunction = getFeatureNameModificationFunction(obj)
			featureNameModificationFunction = prtUtilFeatureNameModificationFunctionHandleCreator('#strIn#_{%s}', obj.nameAbbreviation);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function T_hat = run_tc(params)
  T = params.T;
  Idx = params.Idx;

  Omega = find(Idx);
  Ak = T(Omega);

  N = ndims(T);
  Nway = size(T); % dimension of tensor
  coreNway = [Nway(1),Nway(2),1];

  % rank_dec strategy
  opts = [];
  opts.alpha_adj = 0;
  opts.rank_adj = -1*ones(1,3);
  opts.rank_min = 5*ones(1,3);
  opts.rank_max = 20*ones(1,3);
  EstCoreNway = round(1.25*coreNway);
  coNway = zeros(1,N);
  for n = 1:N
    coNway(n) = prod(Nway)/Nway(n);
  end
  % use random generated starting point
  for i = 1:3
    X0{i} = randn(Nway(i),EstCoreNway(i));
    Y0{i} = randn(EstCoreNway(i),coNway(i));
  end
  opts.X0 = X0; opts.Y0 = Y0;
  [X_dec,Y_dec,Out_dec] = TMac(Ak,Omega,Nway,EstCoreNway,opts);
  % use the weighted sum of all mode matrix factorizations as estimated tensor
  T_hat = zeros(Nway);
  for i = 1:N
    T_hat = T_hat+Out_dec.alpha(i)*Fold(X_dec{i}*Y_dec{i},Nway,i);
  end
end

author: andrewssobral


text:
function [t,u,v,idx,xnode]=raysurf(p0,v0,node,face)
%
% [t,u,v,idx,xnode]=raysurf(p,v,node,face)
%
% perform a Havel-styled ray tracing for a triangular surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   p0: list of starting points of the rays
%   v0: directional vector of the rays,
%   node: a list of node coordinates (nn x 3)
%   face: a surface mesh triangle list (ne x 3)
%
% output:
%   t: distance from p0 to the intersection point for each surface
%      triangle, if t(i)=NaN, no intersection was found for that ray
%   u: bary-centric coordinate 1 of all intersection points
%   v: bary-centric coordinate 2 of all intersection points
%      the final bary-centric triplet is [u,v,1-u-v]
%   idx: idx lists the IDs of the face elements that intersects
%      each ray
%   xnode: optional output, if requested, xnode gives the intersection
%      point coordinates; to compute manually, xnode=p0+repmat(t,1,3).*v0
%
% Reference:
%  [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using
%          SSE4)," IEEE Trans. on Visualization and Computer Graphics,
%          16(3):434-438 (2010)
%  [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous
%          Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
%          press)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=size(p0,1);
if(len==0)
   error('p0 can not be empty');
end
if(size(node,2)<3)
   error('node must contain at least 3 columns');
end
if(size(face,2)<3)
   error('face must contain at least 3 columns');
end

if(size(v0,1)==1 || size(v0,2)==1 && len>1)
   v0=repmat(v0(:)',len,1);
end

t=zeros(len,1)*nan;
u=t;
v=t;
idx=t;

for i=1:len
   [ti,ui,vi,id]=raytrace(p0(i,:),v0(i,:),node,face);
   if(isempty(id)) continue; end
   ti=ti(id);
   tpid=find(ti>=0);
   if(isempty(tpid)) continue; end
   [tmin,tloc]=min(ti(find(ti>=0)));
   t(i)=tmin;
   u(i)=ui(id(tpid(tloc)));
   v(i)=vi(id(tpid(tloc)));
   idx(i)=id(tpid(tloc));
end

if(nargout>=5)
   xnode=p0+repmat(t,1,3).*v0;
end

author: fangq


text:
function outNii = niftiSelect(nii, varargin)
% Preserve data within a spatial cube and zero everything outside it
%
%    boxNii = niftiSelect(nii, varargin)
%
% Typical uses
%    Select a portion of a volume in a white matter mask
%    Extract certain volumes from a 4D data set (e.g., bvals = 2000)
%
% This function writes and saves a new nifti where specified dimensions of
% the data are overwritten with zeros. Useful for unit testing.
%
%  Parameters
%   'nii'          filename or nii struct
%   'saveFlag'   logical
%   'niiNewName' string
%   'keepLR'     vector
%   'keepPA'     vector
%   'keepIS'     vector
%
% Outputs
%   boxNii: the new nifti file (in the form of readFileNifti)
%
% Inputs
%   nii: the nifti file to be copied, and an output of readFileNifti
%   keepLR: the dimensions to keep in the 1st dimension
%   keepPA: the dimensions to keep in the 2nd dimension
%   keepIS: the dimensions to keep in the 3rd dimension
%   niiNewName: Full path to the new nifti file
%   saveFlag:  Save the nii as niiNewName (default: boxNifti)
%
% Examples
%  Extract a small volume from a T1
%    nii = niftiRead(fullfile(vistaRootPath,'rdt','inplane.nii.gz'));
%    niftiView(nii);
%    clear p;  p.keepLR = 16:64; p.keepPA = 96:159; p.keepIS = 5:15;  p.saveFlag = false;
%    boxNii = niftiSelect(nii,p);
%    niftiView(boxNii);
%
%    p.saveFlag = true; p.niiNewName = 'deleteMe';
%    niftiSelect(nii,p);  boxNii = niftiRead(p.niiNewName); niftiView(boxNii);
%
% RL/BW Vistasoft Team, 2016

%% Parameters
p = inputParser;

% Could be a file name or a nifti struct
vFunc = @(x)((ischar(x) && exist(x,'file')) || (isstruct(x) && isfield(x,'nifti_type')));
p.addRequired('nii',vFunc);
if ischar(nii), nii = niftiRead(nii); end

p.addParameter('saveFlag',true,@islogical);
p.addParameter('niiNewName','boxNifti',@ischar);

sz = size(nii.data);
p.addParameter('keepLR',1:sz(1),@isvector);
p.addParameter('keepPA',1:sz(2),@isvector);
p.addParameter('keepIS',1:sz(3),@isvector);

p.parse(nii,varargin{:});
niiNewName = p.Results.niiNewName;
saveFlag   = p.Results.saveFlag;
keepLR = p.Results.keepLR;
keepPA = p.Results.keepPA;
keepIS = p.Results.keepIS;

%% the original nifti
niiDims = nii.dim;

%% the dimension we want to ZERO out
zeroLR = setdiff(1:niiDims(1), keepLR);
zeroPA = setdiff(1:niiDims(2), keepPA);
zeroIS = setdiff(1:niiDims(3), keepIS);

% volume of the original and the cropped
% vol = prod(niiDims);
% volCropped = length(keepLR) * length(keepPA) * length(keepIS);
% fracCrop = volCropped/vol;
% display(['The cropped volume is ' num2str(fracCrop) ' of the original. '])

%% the new data field
newData = nii.data;
newData(zeroLR,:, :) = 0;
newData(:,zeroPA, :) = 0;
newData(:,:, zeroIS) = 0;

%% Make the new nifti struct
outNii = nii;
outNii.fname = niiNewName;
outNii.data  = newData;

if saveFlag, niftiWrite(outNii); end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = sabr_complex(f,k,T,alpha,beta,rho,nu)
% sabr price based on the complex Black-Scholes formula
% from a paper by Benhamou et al.

    f0 = k; % projected fwd as suggested by Benhamou et al. makes z0=0
    z = (f^(1-beta) - k.^(1-beta))./(alpha * (1-beta));
    z0 = (f0.^(1-beta) - k.^(1-beta))/(alpha * (1-beta));

    x = log((-rho + nu*z+sqrt(1-2*nu*rho*z+nu^2*z.^2))./(1-rho))/nu;
    b1 = beta ./(alpha * (1-beta) * z0 + k.^(1-beta));

    theta = 0.25 * rho*nu*alpha*b1.*z.^2 ...
        + log(alpha*(f*k).^(0.5*beta).*z./(f-k)) ...
        + log(x./z.*(1-2*nu*rho*z+nu^2*z.^2).^0.25);

    kappa = 0.125 .* (alpha^2*(beta-2)*beta*k.^(2*beta) ...
        ./(k-alpha*(beta-1)*k.^beta.*z0).^2 ...
        + 6*alpha*beta*k.^beta*nu*rho./(k-alpha*(beta-1)*k.^beta.*z0) ...
        + nu^2*(2-3*rho^2+2*nu*rho*z0-nu^2*z0.^2) ...
        ./(1-2*nu*rho*z0+nu^2*z0.^2));

    integral = sqrt(pi)./(2i*sqrt(kappa))...
        .*(exp( 1i * sqrt(2)* x.* sqrt(kappa)) ...
        .*(erfcomplex(x/sqrt(2*T)+1i*sqrt(kappa*T))-1) ...
        +  exp(-1i * sqrt(2)* x.* sqrt(kappa)) ...
        .*(erfcomplex(-x/sqrt(2*T)+1i*sqrt(kappa*T))+1));

    y = f-k+(f-k)./(2*x*sqrt(2*pi)).* real(exp(theta) .* integral);
end

author: Sable

repo: mcbench-benchmarks


text:
function vargout = subplot_tight(m, n, p, margins, varargin)
    %% subplot_tight
    % A subplot function substitude with margins user tunabble parameter.
    %
    %% Syntax
    %  h=subplot_tight(m, n, p);
    %  h=subplot_tight(m, n, p, margins);
    %  h=subplot_tight(m, n, p, margins, subplotArgs...);
    %
    %% Description
    % Our goal is to grant the user the ability to define the margins between neighbouring
    %  subplots. Unfotrtunately Matlab subplot function lacks this functionality, and the
    %  margins between subplots can reach 40% of figure area, which is pretty lavish. While at
    %  the begining the function was implememnted as wrapper function for Matlab function
    %  subplot, it was modified due to axes del;etion resulting from what Matlab subplot
    %  detected as overlapping. Therefore, the current implmenetation makes no use of Matlab
    %  subplot function, using axes instead. This can be problematic, as axis and subplot
    %  parameters are quie different. Set isWrapper to "True" to return to wrapper mode, which
    %  fully supports subplot format.
    %
    %% Input arguments (defaults exist):
    %   margins- two elements vector [vertical,horizontal] defining the margins between
    %        neighbouring axes. Default value is 0.04
    %
    %% Output arguments
    %   same as subplot- none, or axes handle according to function call.
    %
    %% Issues & Comments
    %  - Note that if additional elements are used in order to be passed to subplot, margins
    %     parameter must be defined. For default margins value use empty element- [].
    %  -
    %
    %% Example
    % close all;
    % img=imread('peppers.png');
    % figSubplotH=figure('Name', 'subplot');
    % figSubplotTightH=figure('Name', 'subplot_tight');
    % nElems=17;
    % subplotRows=ceil(sqrt(nElems)-1);
    % subplotRows=max(1, subplotRows);
    % subplotCols=ceil(nElems/subplotRows);
    % for iElem=1:nElems
    %    figure(figSubplotH);
    %    subplot(subplotRows, subplotCols, iElem);
    %    imshow(img);
    %    figure(figSubplotTightH);
    %    subplot_tight(subplotRows, subplotCols, iElem, [0.0001]);
    %    imshow(img);
    % end
    %
    %% See also
    %  - subplot
    %
    %% Revision history
    % First version: Nikolay S. 2011-03-29.
    % Last update:   Nikolay S. 2012-05-24.
    %
    % *List of Changes:*
    % 2012-05-24
    %  Non wrapping mode (based on axes command) added, to deal with an issue of disappearing
    %     subplots occuring with massive axes.

    %% Default params
    isWrapper=false;
    if (nargin<4) || isempty(margins)
        margins=[0.04,0.04]; % default margins value- 4% of figure
    end
    if length(margins)==1
        margins(2)=margins;
    end

    %note n and m are switched as Matlab indexing is column-wise, while subplot indexing is row-wise :(
    [subplot_col,subplot_row]=ind2sub([n,m],p);

    height=(1-(m+1)*margins(1))/m; % single subplot height
    width=(1-(n+1)*margins(2))/n;  % single subplot width

    % note subplot suppors vector p inputs- so a merged subplot of higher dimentions will be created
    subplot_cols=1+max(subplot_col)-min(subplot_col); % number of column elements in merged subplot
    subplot_rows=1+max(subplot_row)-min(subplot_row); % number of row elements in merged subplot

    merged_height=subplot_rows*( height+margins(1) )- margins(1);   % merged subplot height
    merged_width= subplot_cols*( width +margins(2) )- margins(2);   % merged subplot width

    merged_bottom=(m-max(subplot_row))*(height+margins(1)) +margins(1); % merged subplot bottom position
    merged_left=min(subplot_col)*(width+margins(2))-width;              % merged subplot left position
    pos=[merged_left, merged_bottom, merged_width, merged_height];

    if isWrapper
        h=subplot(m, n, p, varargin{:}, 'Units', 'Normalized', 'Position', pos);
    else
        h=axes('Position', pos, varargin{:});
    end

    if nargout==1
        vargout=h;
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function logAxis(ax,xLimits,yLimits)

if nargin == 0, ax = gca; end
set(ax,'XScale','log','YScale','log');

if (nargin < 2) || isempty(xLimits), xLimits = get(ax,'XLim'); end

if (nargin < 3) || isempty(yLimits), yLimits = get(ax,'YLim'); end

logScale = diff(yLimits)/diff(xLimits);
powerScale = diff(log10(yLimits))/diff(log10(xLimits));

set(ax,'Xlim',xLimits,...
  'YLim',yLimits,...
  'DataAspectRatio',[1 logScale/powerScale 1]);

end

author: mtex-toolbox

repo: mtex

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


text:
%% Copyright (C) 2014, 2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @documentencoding UTF-8
%% @deftypefun  {@var{F} =} fibonacci (@var{n})
%% @deftypefunx {@var{p} =} fibonacci (@var{n}, @var{x})
%% Return symbolic Fibonacci numbers or Fibonacci polynomials.
%%
%% Examples:
%% @example
%% @group
%% fibonacci(15)
%%   @result{} (sym) 610
%% syms n
%% fibonacci(n)
%%   @result{} (sym)
%%       F
%%        n
%% @end group
%% @end example
%%
%% Polynomial example:
%% @example
%% @group
%% syms x
%% fibonacci(10, x)
%%   @result{} (sym)
%%        9      7       5       3
%%       x  + 8⋅x  + 21⋅x  + 20⋅x  + 5⋅x
%% @end group
%% @end example
%%
%% @seealso{euler, bernoulli}
%% @end deftypefun

function r = fibonacci(n, x)

  if (nargin == 1)
    r = pycall_sympy__ ('return sp.fibonacci(*_ins),', sym(n));
  elseif (nargin == 2)
    r = pycall_sympy__ ('return sp.fibonacci(*_ins),', sym(n), sym(x));
  else
    print_usage ();
  end

end

%!assert (isequal ( fibonacci (sym(0)), 0))
%!assert (isequal ( fibonacci (sym(14)), sym(377)))
%!assert (isequal ( fibonacci (14), 377))
%!test syms x
%! assert (isequal (fibonacci (5,x), x^4 + 3*x^2 + 1))

repo: octsympy


text:
clear;

addpath('../PDM_helpers/');

Reconstruct_Torresani;

%% Create the PDM using PCA, from the recovered 3D data
clear
load('Torr_menpo_wild.mat');

% need to still perform procrustes though

x = P3(1:end/3,:);
y = P3(end/3+1:2*end/3,:);

% To make sure that PDM faces the right way (positive Z towards the screen)
z = P3(2*end/3+1:end,:);

[ normX, normY, normZ, meanShape, Transform ] = ProcrustesAnalysis3D(x,y,z, true);
observations = [normX normY normZ];

[princComp, score, eigenVals] = princomp(observations,'econ');
% Keeping most variation
totalSum = sum(eigenVals);
count = numel(eigenVals);
for i=1:numel(eigenVals)
   if ((sum(eigenVals(1:i)) / totalSum) >= 0.999)
      count = i;
      break;
   end
end

V = princComp(:,1:count);
E = eigenVals(1:count);
M = meanShape(:);

% Now normalise it to have actual world scale
lPupil = [(M(37) + M(40))/2; (M(37 + 68) + M(40 + 68))/2];
rPupil = [(M(43) + M(46))/2; (M(43 + 68) + M(46 + 68))/2];

dist = norm(lPupil - rPupil,2);

% average human interocular distance is 65mm
scaling = 65 / dist;

% normalise the mean values as well
M(1:end/3) = M(1:end/3) - mean(M(1:end/3));
M(end/3+1:2*end/3) = M(end/3+1:2*end/3) - mean(M(end/3+1:2*end/3));
M(2*end/3+1:end) = M(2*end/3+1:end) - mean(M(2*end/3+1:end));

M = M * scaling;
E = E * scaling .^ 2;

% orthonormalise V
scalingV = sqrt(sum(V.^2));
V = V ./ repmat(scalingV, numel(M),1);

E = E .* (scalingV') .^ 2;

%% Also align the model to a Multi-PIE one for accurate head orientation
% also align the mean shape model (an aligned 68 point PDM from Multi-PIE dataset)
M_wild = M;

load('pdm_68_multi_pie.mat', 'M');
M_align = M;

n_points = numel(M)/3;

M_wild = reshape(M_wild, n_points, 3);

% work out the translation and rotation needed
[ R, T ] = AlignShapesKabsch(M_wild, reshape(M_align, n_points,3));

% Transform the wild one to be in the same reference as the Multi-PIE one
M_aligned = (R * M_wild')';

M_aligned(:,1) = M_aligned(:,1) + T(1);
M_aligned(:,2) = M_aligned(:,2) + T(2);
M_aligned(:,3) = M_aligned(:,3) + T(3);

M_aligned = M_aligned(:);
V_aligned = V;

% need to align the principal components as well
for i=1:size(V,2)

    V_aligned_pie_curr = (R * reshape(V(:,i), n_points, 3)')';
    V_aligned(:,i) = V_aligned_pie_curr(:);

end

V = V_aligned;
M = M_aligned;

save('pdm_68_aligned_menpo.mat', 'E', 'M', 'V');
writePDM(V, E, M, 'pdm_68_aligned_menpo.txt');

% also save this to model location
if(exist('../../models/pdm/', 'file'))
    save('../../models/pdm/pdm_68_aligned_menpo.mat', 'E', 'M', 'V');
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function check_shape(shape,value,argname,funcname)
% Check whether a given value can be assigned to a field of a given shape.
% check_type(ShapeString,Value,ArgumentName,FunctionName)
%
% This function throws an error message if the check fails.
%
% In:
%   ShapeString : a string that identifies the type of the argument (as in deduce_shape).
%
%   Value : the value to check
%
%   ArgumentName : name of the affected argument, for diagnostic messages.
%
%   FunctionName : name of the function that defines the argument.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2014-02-26

% Copyright (C) Christian Kothe, SCCN, 2014, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

switch shape
    case 'scalar'
        if any(size(value) > 1)
            error('BCILAB:arg:shapecheck','The value assigned to %s in %s must be scalar, but was: %s.',argname,funcname,hlp_tostring(value,1000)); end
    case 'empty'
        if ~isempty(value)
            error('BCILAB:arg:shapecheck','The value assigned to %s in %s must be empty, but was: %s.',argname,funcname,hlp_tostring(value,1000)); end
    case 'row'
        if size(value,1) > 1
            error('BCILAB:arg:shapecheck','The value assigned to %s in %s must be a row vector, but was: %s.',argname,funcname,hlp_tostring(value,1000)); end
    case 'column'
        if size(value,2) > 1
            error('BCILAB:arg:shapecheck','The value assigned to %s in %s must be a column vector, but was: %s.',argname,funcname,hlp_tostring(value,1000)); end
    case 'matrix'
        if ndims(value) > 2 %#ok<ISMAT>
            error('BCILAB:arg:shapecheck','The value assigned to %s in %s must be a matrix, but was: %s.',argname,funcname,hlp_tostring(value,1000)); end
end


text:
function value = lens_half_area_2d ( r, theta1, theta2 )

%*****************************************************************************80
%
%% LENS_HALF_AREA_2D returns the area of a circular half lens in 2D.
%
%  Discussion:
%
%    A circular half lens is formed by drawing a circular arc, and joining its endpoints.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real THETA1, THETA2, the angles of the rays
%    that begin and end the arc.
%
%    Output, real HALF_LENS_AREA_2D, the area of the half lens.
%
  sector = circle_sector_area_2d ( r, theta1, theta2 );
  triangle = circle_triangle_area_2d ( r, theta1, theta2 );
  value = sector - triangle;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function handles = distributionPlot(varargin)
%DISTRIBUTIONPLOT creates violin plots for convenient visualization of multiple distributions
% 2019.03.28 [13:29:26] - Biafra fix % # issue
% SYNOPSIS: handles = distributionPlot(data,propertyName,propertyValue,...)
%           handles = distributionPlot(ah,...)
%
% INPUT data : m-by-nData array of values, or vector of grouped data (use
%           the 'groups' property to specify the grouping variable), or
%           cell array of length nData.
%           The cell array can either contain vectors with values, or
%           m-by-2 arrays with [bins,counts] if you want to determine the
%           histograms by yourself (m can be different between cell
%           elements). Note that arrays inside cells with any
%           other shape than m-by-2 are reshaped to vector an a warning is
%           thrown (DISTRIBUTIONPLOT:AUTORESHAPE).
%
%       DISTRIBUTIONPLOT accepts the following propertyName/propertyValue
%           pairs (all are optional):
%
%       distWidth :  width of distributions; ideally between 0 and 1.
%           1 means that adjacent distributions might touch. Default: 0.9
%       variableWidth : If true, the width of the distribution changes,
%           reflecting the shape of the histogram of the data. If false,
%           the distribution is only encoded by color levels. Default: true
%       color : uniform coloring of histograms. Supply either a color
%           string ('r'), or a truecolor vector ([1 0 0]). Use a
%           cell array of length nData to specify one color per
%           distribution. Default: 'k'
%           If variableWidth is set to false, a colormap is generated that
%           goes from white to the chose color (or from black, if
%           invert==true).
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       colormap : colormap used to describe the distribution (first row
%           corresponds to bins with least data, last row corresponds to
%           bins with most data (invert the grayscale colormap to have
%           black indicate the most data).
%           Supply a cell array of length nData to color distributions
%           individually. Note that using multiple colormaps means that
%           the colorbar doesn't contain much useful information.
%           Default: []
%           Colormap will index into the figure colormap, which will be
%           modified by distributionPlot. This is done to allow editing the
%           distributions in e.g. Adobe Illustrator.
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       globalNorm : normalization for bin width (x-direction)
%           0 : every histogram is normalized individually so that the
%               maximum bin width is equal to distWidth. This is best
%               suited to comparing distribution shapes. Default.
%           1 : histograms are normalized such that equal bin width
%               reports equal numbers of counts per bin.
%           2 : histograms are normalized so that the relative areas
%               covered by the histograms reflect the relative total number
%               of data points.
%           3 : histograms areas are normalized so that relative densities
%               are the same across histograms. Thus, if
%               data = {rand(100,1),rand(500,1)},
%               then
%               distributionPlot(data,'globalNorm',2,'histOpt',0,'divFactor',10)
%               shows the left histogram 5x as wide as the right, while
%               distributionPlot(data,'globalNorm',3,'histOpt',0,'divFactor',10)
%               displays both histograms equally wide, since each bin
%               contains ~10% of the data.
%           Options 1 and 2 produce similar results if the bins are spaced
%           equally for the distributions. Options 0 and 3 produce similar
%           results if the data are drawn from the same distributions.
%           Note that colormaps currently always report the number of data
%           points per bin; 'globalNorm' only applies to the distribution
%           shape.
%
%       groups : grouping variable for grouped data. Grouping will be
%                   resolved by calling grp2idx, and unless xNames have
%                   been supplied, group names determine the x-labels.
%                   If the grouping variable is numeric, group labels also
%                   determine x-values, unless the parameter xValues has
%                   been specified.
%       histOpt : histogram type to plot
%                   0 : use hist command (no smoothing, fixed number of
%                       bins)
%                   1 : smoothened histogram using ksdensity with
%                       Normal kernel. Default.
%                   1.1: smoothened histogram using ksdensity where the
%                       kernel is robustly estimated via histogram.m.
%                       Normal kernel.
%                   2 : histogram command (no smoothing, automatic
%                       determination of thickness (y-direction) of bins)
%       divFactor : Parameter dependent on histOpt. If...
%                   histOpt == 0: divFactor = # of bins. Default: 25.
%                       Alternatively, pass a vector which will be
%                       interpreted as bin centers.
%                   histOpt == 1: divFactor decides by how much the default
%                       kernel-width is multiplied in order to avoid an
%                       overly smooth histogram. Default: 1/2
%                   histOpt == 2: divFactor decides by how much the
%                       automatic bin width is multiplied in order to have
%                       more (<1) or less (>1) detail. Default: 1
%       addSpread : if 1, data points are plotted with plotSpread.
%                   distWidth is ideally set to 0.95
%                   This option is not available if the data is supplied as
%                   histograms.
%                   Please download plotSpread.m separately from the File
%                   Exchange using the link in the remarks
%       showMM : if 1, mean and median are shown as red crosses and
%                green squares, respectively. This is the default
%                2: only mean
%                3: only median
%                4: mean +/- standard error of the mean (no median)
%                5: mean +/- standard deviation (no median)
%                6: draw lines at the 25,50,75 percentiles (no mean)
%                0: plot neither mean nor median
%       xValues: x-coordinate where the data should be plotted.
%                If xValues are given, "distWidth" is scaled by the median
%                difference between adjacent (sorted) x-values. Note that
%                this may lead to overlapping distributions. Default:
%                1:nData
%       xNames : cell array of length nData containing x-tick names
%               (instead of the default '1,2,3')
%       xMode  : if 'auto', x-ticks are spaced automatically. If 'manual',
%                there is a tick for each distribution. If xNames is
%                provided as input, xMode is forced to 'manual'. Default:
%                'manual'.
%          NOTE: SPECIFYING XNAMES OR XVALUES OR XMODE WILL ERASE PREVIOUS
%                LABELS IF PLOTTING INTO EXISTING AXES
%       yLabel : string with label for y-axis. Default : ''
%                If empty and data is histograms, ylabel is set to 'counts'
%       invert : if 1, axes color is changed to black, and colormap is
%                   inverted.
%       histOri: Orientation of histogram. Either 'center', 'left', or
%                'right'. With 'left' or 'right', the left or right half of
%                the standard violin plot is shown. Has no effect if
%                variableWidth is false. Default: center
%       xyOri  : orientation of axes. Either 'normal' (=default), or
%                'flipped'. If 'flipped', the x-and y-axes are switched, so
%                that violin plots are horizontal. Consequently,
%                axes-specific properties, such as 'yLabel' are applied to
%                the other axis.
%       widthDiv : 1-by-2 array with [numberOfDivisions,currentDivision]
%                widthDiv allows cutting the stripe dedicated to a single
%                distribution into multible bands, which can be filled with
%                sequential calls to distributionPlot. This is one way
%                to compare two (or more) sequences of distributions. See
%                example below.
%       ah : axes handle to plot the distributions. Default: gca
%
% OUTPUT handles : 1-by-4 cell array with patch-handles for the
%                  distributions, plot handles for mean/median, the
%                  axes handle, and the plotSpread-points handle
%
%
% EXAMPLES
%         %--Distributions contain more information than boxplot can capture
%         r = rand(1000,1);
%         rn = randn(1000,1)*0.38+0.5;
%         rn2 = [randn(500,1)*0.1+0.27;randn(500,1)*0.1+0.73];
%         rn2=min(rn2,1);rn2=max(rn2,0);
%         figure
%         ah(1)=subplot(3,4,1:2);
%         boxplot([r,rn,rn2])
%         ah(2)=subplot(3,4,3:4);
%         distributionPlot([r,rn,rn2],'histOpt',2); % histOpt=2 works better for uniform distributions than the default
%         set(ah,'ylim',[-1 2])
%
%         %--- additional options
%
%         data = [randn(100,1);randn(50,1)+4;randn(25,1)+8];
%         subplot(3,4,5)
%
%         %--- defaults
%         distributionPlot(data);
%         subplot(3,4,6)
%
%         %--- show density via custom colormap only, show mean/std,
%         distributionPlot(data,'colormap',copper,'showMM',5,'variableWidth',false)
%         subplot(3,4,7:8)
%
%         %--- auto-binwidth depends on # of datapoints; for small n, plotting the data is useful
%         % note that this option requires the additional installation
%         % of plotSpread from the File Exchange (link below)
%         distributionPlot({data(1:5:end),repmat(data,2,1)},'addSpread',true,'showMM',false,'histOpt',2)
%
%         %--- show quantiles
%         subplot(3,4,9),distributionPlot(randn(100,1),'showMM',6)
%
%         %--- horizontal orientation
%         subplot(3,4,10:11),
%         distributionPlot({chi2rnd(3,1000,1),chi2rnd(5,1000,1)},'xyOri','flipped','histOri','right','showMM',0),
%         xlim([-3 13])
%
%         %--- compare distributions side-by-side (see also example below)
%         % plotting into specified axes will throw a warning that you can
%         % turn off using " warning off DISTRIBUTIONPLOT:ERASINGLABELS "
%         ah = subplot(3,4,12);
%         subplot(3,4,12),distributionPlot(chi2rnd(3,1000,1),'histOri','right','color','r','widthDiv',[2 2],'showMM',0)
%         subplot(3,4,12),distributionPlot(chi2rnd(5,1000,1),'histOri','left','color','b','widthDiv',[2 1],'showMM',0)
%
%         %--Use globalNorm to generate meaningful colorbar
%         data = {randn(100,1),randn(500,1)};
%         figure
%         distributionPlot(data,'globalNorm',true,'colormap',1-gray(64),'histOpt',0,'divFactor',[-5:0.5:5])
%         colorbar
%
%         %--Use widthDiv to compare two series of distributions
%         data1 = randn(500,5);
%         data2 = bsxfun(@plus,randn(500,5),0:0.1:0.4);
%         figure
%         distributionPlot(data1,'widthDiv',[2 1],'histOri','left','color','b','showMM',4)
%         distributionPlot(gca,data2,'widthDiv',[2 2],'histOri','right','color','k','showMM',4)
%
%         %--Christmas trees!
%           x=meshgrid(1:10,1:10);
%           xx = tril(x);
%           xx = xx(xx>0);
%           figure
%           hh=distributionPlot({xx,xx,xx},'color','g','addSpread',1,'histOpt',2,'showMM',0);
%           set(hh{4}{1},'color','r','marker','o')
% END
%
% REMARKS To show distributions as clouds of points (~beeswarm plot),
%         and/or to use the option "addSpread", please download the
%         additional function plotSpread.m from the File Exchange
%         http://www.mathworks.com/matlabcentral/fileexchange/37105-plot-spread-points-beeswarm-plot
%
%         I used to run ksdensity with the Epanechnikov kernel. However,
%         for integer data, the shape of the kernel can produce peaks
%         between the integers, which is not ideal (use histOpt=2 for
%         integer valued data).
%
%         A previous iteration of distributionPlot used the input
%         specifications below. They still work to ensure backward
%         compatibility, but are no longer supported or updated.
%           handles = distributionPlot(data,distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm)
%           where distWidth of 1 means that the maxima
%           of  two adjacent distributions might touch. Negative numbers
%           indicate that the distributions should have constant width, i.e
%           the density is only expressed through greylevels.
%           Values between 1 and 2 are like values between 0 and 1, except
%           that densities are not expressed via graylevels. Default: 1.9
%
%
% SEE ALSO histogram, ksdensity, plotSpread, boxplot, grp2idx
%

% created with MATLAB ver.: 7.6.0.324 (R2008a) on Windows_NT
%
% created by: Jonas Dorn; jonas.dorn@gmail.com
% DATE: 08-Jul-2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%====================================
%% TEST INPUT
%====================================

% set defaults
def.xNames = [];
def.showMM = 1;
def.distWidth = 0.9;
def.histOpt = 1;
def.divFactor = [25,2,1];
def.invert = false;
def.colormap = [];
def.color = 'k';
def.addSpread = false;
def.globalNorm = false;
def.variableWidth = true;
def.groups = [];
def.yLabel = '';
def.xValues = '';
def.xMode = 'manual';
def.histOri = 'center';
def.xyOri = 'normal';
def.widthDiv = [1 1];
isHistogram = false; % # this parameter is not set by input

if nargin == 0 || isempty(varargin{1})
    error('not enough input arguments')
end

% check for axes handle
if ~iscell(varargin{1}) && isscalar(varargin{1}) == 1 && ...
        ishandle(varargin{1}) && strcmp(get(varargin{1},'Type'),'axes')
    ah = varargin{1};
    data = varargin{2};
    varargin(1:2) = [];
    newAx = false;

else
    ah = gca;
    data = varargin{1};
    varargin(1) = [];
    newAx = true;
end

% check for current axes limits. Set NaN if the axes have no children
% yet - we need that in case we're building a complicated set of
% distributions
if ~isempty(get(ah,'children'))
    xAxLim = xlim;
    yAxLim = ylim;
else
    [xAxLim,yAxLim] = deal([NaN NaN]);
end

fh = get(ah,'Parent');

% check data. If not cell, convert
if ~iscell(data)
    [nPoints,nData] = size(data);
    data = mat2cell(data,nPoints,ones(nData,1));
else
    % get nData
    data = data(:);
    nData = length(data);
    % make sure all are vectors
    badCol = ~cellfun(@isvector,data) & ~cellfun(@isempty,data);
    if any(badCol)
        nCols = cellfun(@(x)(size(x,2)),data(badCol));
        if all(nCols==2)
            % bins,counts
            isHistogram = true;
        else
            warning('DISTRIBUTIONPLOT:AUTORESHAPE',...
                'Elements %s of the cell array are not vectors. They will be reshaped automatically',...
                num2str(find(badCol)'));
            data(badCol) = cellfun(@(x)(x(:)),data(badCol),'UniformOutput',false);
        end
    end
end

parserObj = inputParser;
parserObj.FunctionName = 'distributionPlot';
stdWidth = 1; % scaling parameter for variableWidth with uneven x-values
% check whether we're dealing with pN/pV or straight arguments
if ~isempty(varargin) && ~ischar(varargin{1}) && ~isstruct(varargin{1})
    % use old format
    % distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm
    def.distWidth = 1.9;
    parserObj.addOptional('distWidth',def.distWidth);
    parserObj.addOptional('showMM',def.showMM);
    parserObj.addOptional('xNames',def.xNames);
    parserObj.addOptional('histOpt',def.histOpt);
    parserObj.addOptional('divFactor',def.divFactor);
    parserObj.addOptional('invert',def.invert);
    parserObj.addOptional('addSpread',def.addSpread);
    parserObj.addOptional('globalNorm',def.globalNorm);
    parserObj.addOptional('groups',def.groups);
    parserObj.addOptional('yLabel',def.yLabel);
    parserObj.addOptional('color',def.color);

    parserObj.parse(varargin{:});
    opt = parserObj.Results;
    % fill in defaults that are not supported in the old version of the
    % code
    opt.colormap = [];
    opt.variableWidth = true;
    opt.histOri = 'center';
    opt.xValues = [];
    opt.xMode = 'auto';
    opt.xyOri = 'normal';
    opt.widthDiv = [1 1];

    % overwrite empties with defaults - inputParser considers empty to be a
    % valid input.
    fnList = fieldnames(opt);
    for fn = fnList'
        if isempty(opt.(fn{1}))
            opt.(fn{1}) = def.(fn{1});
        end
    end

    % fix a few parameters
    if opt.distWidth > 1
        opt.distWidth = opt.distWidth - 1;
    else
        opt.colormap = 1-gray(128);
    end
    if opt.distWidth < 0
        opt.variableWidth = false;
        opt.distWidth = abs(opt.distWidth);
    end

    if ~isempty(opt.xNames)
        opt.xMode = 'manual';
    end

else
    defNames = fieldnames(def);
    for dn = defNames(:)'
        parserObj.addParamValue(dn{1},def.(dn{1}));
    end

    parserObj.parse(varargin{:});
    opt = parserObj.Results;

    % if groups: deal with data
    if ~isempty(opt.groups)
        [idx,labels,vals] = grp2idx(opt.groups);
        % convert data to cell array
        data = accumarray(idx,data{1},[],@(x){x});
        nData = length(data);
        % if not otherwise provided, use group labels for xnames
        if isempty(opt.xNames)
            opt.xNames = labels;
            if ~iscell(opt.xNames)
                opt.xNames = num2cell(opt.xNames);
            end
        end
        if isnumeric(vals) && isempty(opt.xValues)
            opt.xValues = vals;
        end

    end

    if ~ischar(opt.xyOri) || ~any(ismember(opt.xyOri,{'normal','flipped'}))
        error('option xyOri must be either ''normal'' or ''flipped'' (is ''%s'')',opt.xyOri);
    end

end
% common checks

% default x-values: 1:n
if isempty(opt.xValues)
    opt.xValues = 1:nData;
elseif length(opt.xValues) ~= nData
    error('please supply as many x-data values as there are data entries')
elseif length(opt.xValues) > 1 % only check for scale if more than 1 value
    % scale width
    stdWidth = median(diff(sort(opt.xValues)));
    opt.distWidth = opt.distWidth * stdWidth;
end

if ~isscalar(opt.divFactor) && length(opt.divFactor) == 3 && all(opt.divFactor==def.divFactor)
    opt.divFactor = opt.divFactor(floor(opt.histOpt)+1);
end
if isHistogram
    opt.histOpt = 99;
    if isempty(opt.yLabel)
        opt.yLabel = 'counts';
    end
end

% check colors/colormaps: do we need to expand colormap?
if ~iscell(opt.colormap)
    opt.colormap = {opt.colormap};
end
if ~iscell(opt.color)
    opt.color = {opt.color};
end
for iColor = 1:length(opt.color)
    if ischar(opt.color{iColor})
        opt.color{iColor} = colorCode2rgb(opt.color{iColor});
    end
end

% expand - if only single colormap specified, we expand only once
if ~opt.variableWidth
    missingColormaps = find(cellfun(@isempty,opt.colormap));
    for iMissing = missingColormaps(:)'

        endColor = opt.color{max(iMissing,length(opt.color))};
        % normally, we go from white to color
        cmap = zeros(128,3);
        for rgb = 1:3
            cmap(:,rgb) = linspace(1,endColor(rgb),128);
        end
        opt.colormap{iMissing} = cmap;

    end
end

% if we have colormaps, we need to create a master which we add to the
% figure. Invert if necessary, and expand the cell array to nData
colormapLength = cellfun(@(x)size(x,1),opt.colormap);
if any(colormapLength>0)

    colormap = cat(1,opt.colormap{:});
    if opt.invert
        colormap = 1-colormap;
    end
    set(fh,'Colormap',colormap)
    if length(opt.colormap) == 1
        opt.colormap = repmat(opt.colormap,nData,1);
        colormapLength = repmat(colormapLength,nData,1);
        colormapOffset = zeros(nData,1);
        singleMap = true;
    else
        colormapOffset = [0;cumsum(colormapLength(1:end-1))];
        singleMap = false;
    end

else

    colormapLength = zeros(nData,1);
    if length(opt.color) == 1
        opt.color = repmat(opt.color,nData,1);
    end
    if opt.invert
        opt.color = cellfun(@(x)1-x,opt.color,'uniformOutput',false);
    end
end

% set hold on
holdState = get(ah,'NextPlot');
set(ah,'NextPlot','add');

% if new axes: invert
if newAx && opt.invert
    set(ah,'Color','k')
end

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

%===================================
%% PLOT DISTRIBUTIONS
%===================================

% assign output
hh = NaN(nData,1);
[m,md,sem,sd] = deal(nan(nData,1));
if opt.showMM == 6
    md = nan(nData,3,3); % md/q1/q3, third dim is y/xmin/xmax
end

% get base x-array
% widthDiv is a 1-by-2 array with
% #ofDivs, whichDiv
% The full width (distWidth) is split into
% #ofDivs; whichDiv says which "stripe" is active
xWidth = opt.distWidth/opt.widthDiv(1);
xMin = -opt.distWidth/2;
xLow = xMin + xWidth * (opt.widthDiv(2)-1);
xBase = [-xWidth;xWidth;xWidth;-xWidth]/2;
xOffset = xLow + xWidth/2;

% b/c of global norm: loop twice
plotData = cell(nData,2);

% loop through data. Prepare patch input, then draw patch into gca
for iData = 1:nData
    currentData = data{iData};
    % only plot if there is some finite data
    if ~isempty(currentData(:)) && any(isfinite(currentData(:)))

        switch floor(opt.histOpt)
            case 0
                % use hist
                [xHist,yHist] = hist(currentData,opt.divFactor);

            case 1
                % use ksdensity

                if opt.histOpt == 1.1
                    % use histogram to estimate kernel
                    [dummy,x] = histogramDistributionPlot(currentData); %#ok<ASGLU>
                    if length(x) == 1
                        % only one value. Make fixed distribution
                        dx = 0.1;
                        yHist = x;
                        xHist = sum(isfinite(currentData));
                    else
                        dx = x(2) - x(1);

                    % make sure we sample frequently enough
                    x = min(x)-dx:dx/3:max(x)+dx;
                    [xHist,yHist] = ksdensity(currentData,x,'kernel','normal','width',dx/(1.5*opt.divFactor));
                    end
                else

                    % x,y are switched relative to normal histogram
                    [xHist,yHist,u] = ksdensity(currentData,'kernel','normal');
                    % take smaller kernel to avoid over-smoothing
                    if opt.divFactor ~= 1
                        [xHist,yHist] = ksdensity(currentData,'kernel','normal','width',u/opt.divFactor);
                    end
                end

                % modify histogram such that the sum of bins (not the
                % integral under the curve!) equals the total number of
                % observations, in order to be comparable to hist
                xHist = xHist/sum(xHist)*sum(isfinite(currentData));

            case 2
                % use histogram - bar heights are counts as in hist
                [xHist,yHist] = histogramDistributionPlot(currentData,opt.divFactor,0);
            case 99
                % bins,counts already supplied
                xHist = currentData(:,2)';
                yHist = currentData(:,1)';
        end
        plotData{iData,1} = xHist;
        plotData{iData,2} = yHist;
    end
end

goodData = find(~cellfun(@isempty,plotData(:,1)));
% get norm
switch opt.globalNorm
    case 3
        % #3 normalizes relative densities
        xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
        xNorm(goodData) = xNorm(goodData) .* cellfun(@sum,plotData(goodData,1))';
        maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
        xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));

    case 2
        % #2 should normalize so that the integral of the
        % different histograms (i.e. area covered) scale with the
        % respective sum of counts across all bins. Requires evenly spaced
        % histograms at the moment
        xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
        maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
        xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));
    case 1
        xNorm(goodData) = max(cat(2,plotData{:,1}));
    case 0
        xNorm(goodData) = cellfun(@max,plotData(goodData,1));
end

for iData = goodData'

    % find current data again
    currentData = data{iData};

    xHist = plotData{iData,1};
    yHist = plotData{iData,2};

    % find y-step
    dy = min(diff(yHist));
    if isempty(dy)
        dy = 0.1;
    end

    % create x,y arrays
    nPoints = length(xHist);
    xArray = repmat(xBase,1,nPoints);
    yArray = repmat([-0.5;-0.5;0.5;0.5],1,nPoints);

    % x is iData +/- almost 0.5, multiplied with the height of the
    % histogram
    if opt.variableWidth

        tmp = xArray.*repmat(xHist,4,1)./xNorm(iData);

        switch opt.histOri
            case 'center'
                % we can simply use xArray
                xArray = tmp;
            case 'right'
                % shift everything to the left
                delta = tmp(1,:) - xArray(1,:);
                xArray = bsxfun(@minus,tmp,delta);
            case 'left'
                % shift everything to the right
                delta = tmp(1,:) - xArray(1,:);
                xArray = bsxfun(@plus,tmp,delta);
        end

        xArray = xArray + opt.xValues(iData);

    else
        xArray = xArray + iData;
    end

    % add offset (in case we have multiple widthDiv)
    xArray = xArray + xOffset;

    % yData is simply the bin locations
    yArray = repmat(yHist,4,1) + dy*yArray;

    % add patch
    vertices = [xArray(:),yArray(:)];
    faces = reshape(1:numel(yArray),4,[])';

    if colormapLength(iData) == 0
        colorOpt = {'FaceColor',opt.color{iData}};
    else
        % calculate index into colormap
        if singleMap
            % use scaled mapping so that colorbar is meaningful
            if opt.globalNorm > 0
                colorOpt = {'FaceVertexCData',xHist','CDataMapping','scaled','FaceColor','flat'};
            else
                colorOpt = {'FaceVertexCData',xHist'/xNorm(iData),'CDataMapping','scaled','FaceColor','flat'};
            end

        else
            idx = round((xHist/xNorm(iData))*(colormapLength(iData)-1))+1;
            colorOpt = {'FaceVertexCData',idx'+colormapOffset(iData),'CDataMapping','direct','FaceColor','flat'};
        end
    end

    switch opt.xyOri
        case 'normal'
            hh(iData)= patch('Vertices',vertices,'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
        case 'flipped'
            hh(iData)= patch('Vertices',vertices(:,[2,1]),'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
    end

    if opt.showMM > 0
        if isHistogram
            [m(iData),sem(iData)] = weightedStats(currentData(:,1),currentData(:,2),'w');
            sd(iData) = sem(iData) * sqrt(sum(currentData(:,2)));
            % weighted median: where we're at middle weight
            % may need some tweaking
            goodCurrentData = sortrows(currentData(all(isfinite(currentData),2),:),1);
            weightList = cumsum(goodCurrentData(:,2));
            weightList = weightList / weightList(end);
            md(iData) = goodCurrentData(find(weightList>0.5,1,'first'),1);
        else
            m(iData) = nanmean(currentData);
            md(iData) = nanmedian(currentData);
            sd(iData) = nanstd(currentData);
            sem(iData) = sd(iData)/sqrt(sum(isfinite(currentData)));
        end

        if opt.showMM == 6
            % read quantiles - "y"-value, plus x-start-stop
            % re-use md array which allows using a loop below instead of
            % lots of copy-paste
            % md array is md/q1/q3, with third dimension y/xmin/xmax

            md(iData,2,1) = prctile(currentData,25);
            md(iData,3,1) = prctile(currentData,75);

            for qq = 1:3
                % find corresponding y-bin
                yLoc =  repmat(...
                    any(yArray>md(iData,qq,1),1) & any(yArray<=md(iData,qq,1),1),...
                    [4 1]);
                % look up corresponding x-values. Note that there is a bit
                % of a risk that the line will be exactly between two very
                % different bins - but if we make the line longer, it will
                % be ugly almost all the time
                md(iData,qq,2) = min( xArray( yLoc ) );
                md(iData,qq,3) = max( xArray( yLoc ) );
            end

        end
    end
end % loop

sh = [];
if opt.addSpread
    if isHistogram
        disp('Option addSpread is unavailable if data is supplied as histograms. Call plotSpread separately')
    else
        % add spread
        try
            sh = plotSpread(ah,data,'xValues',opt.xValues,'xyOri',opt.xyOri);
            set(sh{1},'color',[0,128,255]/255);
        catch me
            if strcmp(me.identifier,'MATLAB:UndefinedFunction')
                error('plotSpread not found. Please download it from the Matlab File Exchange')
            else
                rethrow(me)
            end
        end
    end
end

mh = [];mdh=[];
if opt.showMM
    % plot mean, median. Mean is filled red circle, median is green square
    % I don't know of a very clever way to flip xy and keep everything
    % readable, thus it'll be copy-paste
    switch opt.xyOri
        case 'normal'
            if any(opt.showMM==[1,2])
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
            end
            if any(opt.showMM==[1,3])
                mdh = plot(ah,opt.xValues+xOffset,md,'sg','MarkerSize',12);
            end
            if opt.showMM == 4
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,opt.xValues+xOffset,m,sem);
            end
            if opt.showMM == 5
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,opt.xValues+xOffset,m,sd);
            end
            if opt.showMM == 6
                mdh(1,:) = plot(ah,squeeze(md(:,1,2:3))',repmat(md(:,1,1)',2,1),'color','r','lineWidth',2);%,'lineStyle','--');
                mdh(2,:) = plot(ah,squeeze(md(:,2,2:3))',repmat(md(:,2,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
                mdh(3,:) = plot(ah,squeeze(md(:,3,2:3))',repmat(md(:,3,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
            end
        case 'flipped'
            if any(opt.showMM==[1,2])
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
            end
            if any(opt.showMM==[1,3])
                mdh = plot(ah,md,opt.xValues+xOffset,'sg','MarkerSize',12);
            end
            if opt.showMM == 4
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sem,NaN(size(sem))]);
            end
            if opt.showMM == 5
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sd,NaN(size(sd))]);
            end
            if opt.showMM == 6
                mdh(1,:) = plot(ah,repmat(md(:,1,1)',2,1),squeeze(md(:,1,2:3))','color','r','lineWidth',2);%,'lineStyle','--');
                mdh(2,:) = plot(ah,repmat(md(:,2,1)',2,1),squeeze(md(:,2,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
                mdh(3,:) = plot(ah,repmat(md(:,3,1)',2,1),squeeze(md(:,3,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
            end
    end
end

% find extents of x-axis (or y-axis, if flipped)
minX = min(opt.xValues)-stdWidth;
maxX = max(opt.xValues)+stdWidth;

if ~isnan(xAxLim(1))
    % we have previous limits
    switch opt.xyOri
        case 'normal'
            minX = min(minX,xAxLim(1));
            maxX = max(maxX,xAxLim(2));
        case 'flipped'
            minX = min(minX,yAxLim(1));
            maxX = max(maxX,yAxLim(2));
    end
end

% if ~empty, use xNames
switch opt.xyOri
    case 'normal'
        switch opt.xMode
            case 'manual'
                if newAx == false
                    warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
                end
                set(ah,'XTick',opt.xValues);
                if ~isempty(opt.xNames)
                    set(ah,'XTickLabel',opt.xNames)
                end
            case 'auto'
                % no need to do anything
        end
        if ~isempty(opt.yLabel)
            ylabel(ah,opt.yLabel);
        end
        % have plot start/end properly
        xlim([minX,maxX])
    case 'flipped'
        switch opt.xMode
            case 'manual'
                if newAx == false
                    warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
                end
                set(ah,'YTick',opt.xValues);
                if ~isempty(opt.xNames)
                    set(ah,'YTickLabel',opt.xNames)
                end
            case 'auto'
                % no need to do anything
        end
        if ~isempty(opt.yLabel)
            xlabel(ah,opt.yLabel);
        end
        % have plot start/end properly
        ylim([minX,maxX])
end

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

%==========================
%% CLEANUP & ASSIGN OUTPUT
%==========================

if nargout > 0
    handles{1} = hh;
    handles{2} = [mh;mdh];
    handles{3} = ah;
    handles{4} = sh;
end

set(ah,'NextPlot',holdState);


text:
function [priormeans,posteriormeans,covmatrix] = gpnddisimPredictRNAOnly(model,predtimes);

% GPASIMPREDICT Compute predictions (means and a covariance matrix)
% of RNA values for the GPASIM model, conditional on the existing
% RNA values in the model but not on any existing POL2 values.
%
% FORMAT
%---------------------------------
% DESC computes predictions for the asynchronous Gaussian
% process single input motif model.
%
% ARG model : the model for which the gradient is computed.
%
% ARG pol2times : the time points where predictions for POL2 are needed
%
% ARG rnatime : the time points where predictions for RNA are needed
%
% RETURN means : the predicted mean values, first the POL2
% predictions and then the RNA predictions.
%
% RETURN covmatrix : the covariance matrix between the
% predictions; for example, the diagonal values are the variances
% of each prediction.
%---------------------------------
%
% SEEALSO : gpasimCreate
%
% COPYRIGHT : Jaakko Peltonen, 2011

% GPASIMPREDICT

numGenes=model.numGenes;
if numGenes==0,
  error('gpnddisimPredictRNAOnly requires a model that has parameters for RNA modeling\n');
end;

% compute prior means
%pol2priormeans=[1:size(predtimes,1)]'*0;
pol2priormeans=ones(size(predtimes,1),1)*model.simMean;

% Mean for the mRNA is nonconstant over time and depends on the
% B,D,S parameters and on the POL2 mean
Bj=model.B(1);
Dj=model.D(1);
Sj=model.S(1);
if model.use_disimstartmean==1,
  disimStartMean=model.disimStartMean(1);
end;

% compute the RNA mean curve

rnapriormeans=[];
tempind1=1;
for k=1:numGenes,
  nt=length(predtimes);
  rnapriormeans=[rnapriormeans;nan*ones(nt,1)];
  if (model.use_disimstartmean==1),
    tempt=predtimes-model.delay(k);
    I=find(tempt<0);
    tempt(I)=0;
    rnapriormeans(tempind1:tempind1+nt-1)=...
        model.disimStartMean(k)*exp(model.D(k)*(-predtimes)) ...
        +(model.B(k)/model.D(k))*(1-exp(-model.D(k)*predtimes)) ...
        +(model.simMean*model.S(k)/model.D(k))*(1-exp(-model.D(k)*tempt));
  else
    tempt=predtimes-model.delay(k);
    I=find(tempt<0);
    tempt(I)=0;
    rnapriormeans(tempind1:tempind1+nt-1)=...
        ((model.B(k)+model.simMean*model.S(k))/model.D(k))*exp(model.D(k)*(-predtimes))...
        +((model.B(k)+model.simMean*model.S(k))/model.D(k))*(1-exp(-model.D(k)*tempt));
  end;
  tempind1=tempind1+nt;
end;

%if model.use_disimstartmean==1,
%  rnapriormeans=(Bj+model.simMean*Sj)/Dj+(disimStartMean-(Bj+model.simMean*Sj)/Dj)*exp(Dj*(-predtimes));
%else
%  rnapriormeans=(Bj+model.simMean*Sj)/Dj*ones(size(predtimes));
%end;

if 1,
K_new=kernCompute(model.kern, predtimes);
K_new=K_new(length(predtimes)+1:2*length(predtimes),length(predtimes)+1:2*length(predtimes));

K_new_old=kernCompute(model.kern, predtimes, model.t);
K_new_old=K_new_old(length(predtimes)+1:2*length(predtimes),length(model.t)+1:2*length(model.t));

K_old=model.K(length(model.t)+1:2*length(model.t),length(model.t)+1:2*length(model.t));
K_old_new=K_new_old';
end;

tempm=model.m(length(model.t)+1:2*length(model.t));

priormeans=rnapriormeans;
size(tempm)
size(K_new_old)
size(K_old)
size(priormeans)
posteriormeans=priormeans+K_new_old*(K_old\tempm);
%covmatrix=K_new-K_new_old*inv(K_old)*K_old_new;

covmatrix=K_new-K_new_old*(K_old\K_old_new);

%posteriormeans=real(posteriormeans);
%covmatrix=real(covmatrix);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%KDGAUSS Derivative of Gaussian kernel
%
% K = KDGAUSS(SIGMA) is a 2-dimensional derivative of Gaussian kernel (WxW)
% of width (standard deviation) SIGMA and centred within the matrix K whose
% half-width H = 3xSIGMA and W=2xH+1.
%
% K = KDGAUSS(SIGMA, H) as above but the half-width is explictly specified.
%
% Notes::
% - This kernel is the horizontal derivative of the Gaussian, dG/dx.
% - The vertical derivative, dG/dy, is K'.
% - This kernel is an effective edge detector.
%
% See also KGAUSS, KDOG, KLOG, ISOBEL, ICONV.

% Copyright (C) 1993-2011, 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 m = kdgauss(sigma, w)

    if nargin == 1,
        w = ceil(3*sigma);
    end
    ww = 2*w + 1;

    [x,y] = meshgrid(-w:w, -w:w);

    % This should properly be
    %   m = -x/sigma^4 /(2*pi) .*  exp( -(x.^2 + y.^2)/2/sigma^2);
    % but the effect of the error is simply to scale the result by sigma^2.
    %
    % Too many results in the book depend on this...

    m = -x/sigma^2 /(2*pi) .*  exp( -(x.^2 + y.^2)/2/sigma^2);

author: petercorke

repo: machinevision-toolbox-matlab


text:
close all;
seqNames = {'14_14', '14_06', '14_20','13_29', '13_30', '13_31'};
jogTsteps = { [42:82], '', '', '', [165:185], [145:165]};

jobNames = {'SM+zDD+AnnealLin', 'SM+zDD', 'SM+cDD+AnnealLin', 'SM+cDD', 'Prior'};
jobIDs = [2134873 2134874 2134875 2134876 1515159];

EXPORT_DIR = '/home/mhughes/git/BPARHMM-NEW/figs/Mocap6/';
taskID = 1;

jj =0;
for jobID = jobIDs( [1 end] );
    jj = jj+1;
    ss = 0;
    for seqID = [1 5 6]
        ss = ss + 1;
        plotAlignedStateSeqSegment( jobID, taskID, seqNames{seqID}, jogTsteps{ seqID }  );

        fname = fullfile(EXPORT_DIR, sprintf('JogSegmentation_seq%d_%s', seqID, jobNames{jj}) );
        export_fig( fname, '-eps');
        fprintf( 'Exporting %s\n', fname );
    end
end

author: michaelchughes

repo: NPBayesHMM


text:
function ft_realtime_fileproxy(cfg)

% FT_REALTIME_FILEPROXY reads continuous data from an EEG/MEG file and writes it to a
% FieldTrip buffer. This works for any file format that is supported by FieldTrip.
%
% The FieldTrip buffer is a network transparent server that allows the acquisition
% client to stream data to it. An analysis client can connect to read the data upon
% request. Multiple clients can connect simultaneously, each analyzing a specific
% aspect of the data concurrently.
%
% Use as
%   ft_realtime_fileproxy(cfg)
% with the following configuration options
%   cfg.minblocksize         = number, in seconds (default = 0)
%   cfg.maxblocksize         = number, in seconds (default = 1)
%   cfg.channel              = cell-array, see FT_CHANNELSELECTION (default = 'all')
%   cfg.jumptoeof            = jump to end of file at initialization (default = 'no')
%   cfg.readevent            = whether or not to copy events (default = 'no'; event type can also be specified; e.g., 'UPPT002')
%   cfg.speed                = relative speed at which data is written (default = inf)
%
% The source of the data is configured as
%   cfg.source.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.source.datafile      = string
%   cfg.source.headerfile    = string
%   cfg.source.eventfile     = string
%   cfg.source.dataformat    = string, default is determined automatic
%   cfg.source.headerformat  = string, default is determined automatic
%   cfg.source.eventformat   = string, default is determined automatic
%
% The target to write the data to is configured as
%   cfg.target.datafile      = string, target destination for the data (default = 'buffer://localhost:1972')
%   cfg.target.dataformat    = string, default is determined automatic
%
% To stop this realtime function, you have to press Ctrl-C
%
% See also FT_REALTIME_SIGNALPROXY, FT_REALTIME_SIGNALVIEWER

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

% set the defaults
if ~isfield(cfg, 'source'),               cfg.source = [];                                  end
if ~isfield(cfg, 'target'),               cfg.target = [];                                  end
if ~isfield(cfg.source, 'headerformat'),  cfg.source.headerformat = [];                     end % default is detected automatically
if ~isfield(cfg.source, 'dataformat'),    cfg.source.dataformat = [];                       end % default is detected automatically
if ~isfield(cfg.target, 'headerformat'),  cfg.target.headerformat = [];                     end % default is detected automatically
if ~isfield(cfg.target, 'dataformat'),    cfg.target.dataformat = [];                       end % default is detected automatically
if ~isfield(cfg.target, 'datafile'),      cfg.target.datafile = 'buffer://localhost:1972';  end
if ~isfield(cfg, 'minblocksize'),         cfg.minblocksize = 0;                             end % in seconds
if ~isfield(cfg, 'maxblocksize'),         cfg.maxblocksize = 1;                             end % in seconds
if ~isfield(cfg, 'channel'),              cfg.channel = 'all';                              end
if ~isfield(cfg, 'jumptoeof'),            cfg.jumptoeof = 'no';                             end % jump to end of file at initialization
if ~isfield(cfg, 'readevent'),            cfg.readevent = 'no';                             end % capture events?
if ~isfield(cfg, 'speed'),                cfg.speed = inf ;                                 end % inf -> run as fast as possible

% translate dataset into datafile+headerfile
cfg.source = ft_checkconfig(cfg.source, 'dataset2files', 'yes');
cfg.target = ft_checkconfig(cfg.target, 'dataset2files', 'yes');
ft_checkconfig(cfg.source, 'required', {'datafile' 'headerfile'});
ft_checkconfig(cfg.target, 'required', {'datafile' 'headerfile'});

if ~isfield(cfg.source,'eventfile') || isempty(cfg.source.eventfile)
    cfg.source.eventfile = cfg.source.datafile;
end

if ~isfield(cfg.target,'eventfile') || isempty(cfg.target.eventfile)
    cfg.target.eventfile = cfg.target.datafile;
end

% ensure that the persistent variables related to caching are cleared
clear ft_read_header
% read the header for the first time
hdr = ft_read_header(cfg.source.headerfile);
fprintf('updating the header information, %d samples available\n', hdr.nSamples*hdr.nTrials);

% define a subset of channels for reading
cfg.channel = ft_channelselection(cfg.channel, hdr.label);
chanindx    = match_str(hdr.label, cfg.channel);
nchan       = length(chanindx);
if nchan==0
    ft_error('no channels were selected');
end

minblocksmp = round(cfg.minblocksize*hdr.Fs);
minblocksmp = max(minblocksmp, 1);
maxblocksmp = round(cfg.maxblocksize*hdr.Fs);
count       = 0;

if strcmp(cfg.jumptoeof, 'yes')
    prevSample = hdr.nSamples * hdr.nTrials;
else
    prevSample  = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

evt = [];
while true

    % determine number of samples available in buffer
    hdr = ft_read_header(cfg.source.headerfile, 'cache', true);

    % see whether new samples are available
    newsamples = (hdr.nSamples*hdr.nTrials-prevSample);

    if newsamples>=minblocksmp

        begsample  = prevSample+1;
        endsample  = prevSample+min(newsamples,maxblocksmp);

        % remember up to where the data was read
        count       = count + 1;
        fprintf('processing segment %d from sample %d to %d\n', count, begsample, endsample);

        % read data segment
        dat = ft_read_data(cfg.source.datafile,'header', hdr, 'dataformat', cfg.source.dataformat, 'begsample', begsample, 'endsample', endsample, 'chanindx', chanindx, 'checkboundary', false);

        % it only makes sense to read those events associated with the currently processed data
        if ~strcmp(cfg.readevent,'no')
          evt = ft_read_event(cfg.source.eventfile, 'header', hdr, 'minsample', begsample, 'maxsample', endsample);

          if ~strcmp(cfg.readevent,'yes')
            evt = ft_filter_event(evt, 'type', cfg.readevent);
          end

        end

        prevSample  = endsample;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % from here onward it is specific to writing the data to another stream
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if count==1
            % the input file may have a different offset than the output file
            offset = begsample - 1;
            % flush the file, write the header and subsequently write the data segment
            ft_write_data(cfg.target.datafile, dat, 'header', hdr, 'dataformat', cfg.target.dataformat, 'chanindx', chanindx, 'append', false);
            if ~strcmp(cfg.readevent,'no')
                for i=1:numel(evt)
                    evt(i).sample = evt(i).sample - offset;
                end
                ft_write_event(cfg.target.eventfile,evt,'append',false);
            end
        else
            % write the data segment
            ft_write_data(cfg.target.datafile, dat, 'header', hdr, 'dataformat', cfg.target.dataformat, 'chanindx', chanindx, 'append', true);
            if ~strcmp(cfg.readevent,'no')
                for i=1:numel(evt)
                    evt(i).sample = evt(i).sample - offset;
                end
                ft_write_event(cfg.target.eventfile,evt,'append',true);
            end
        end % if count==1

        % wait for a realistic amount of time
        pause(((endsample-begsample+1)/hdr.Fs)/cfg.speed);

    end % if enough new samples
end % while true


text:
function gauss = get_gauss_kernel_by_sigma(sigma, d)
% generate gaussian kernel with given pixel resolution,
% parameter d controls nonsingular dimensions.
    wd = round(8 * max(sigma));
    N = round((wd-1)/2);
    x = -N:N;
    size(x)

    ndims = numel(sigma);

    y = cell(ndims, 1);
    for i = 1 : ndims
        y{i} = 1/(sqrt(2*pi)*sigma(i)) * exp( - (x).^2/(2*sigma(i)^2) );
        if d(i) == 0
            y{i} = y{i}*0;
            y{i}((numel(y{i}) - 1)/2 + 1) = 1;
        end
    end

    if ndims == 2
%         K = y{2}' * y{1};
        K = y{1}' * y{2};
    elseif ndims == 3
        t = y{1}' * y{2};
        K = zeros((2*N+1)*[1,1,1]);
        for i = 1 : 2*N+1
            K(:,:,i) = t * y{3}(i);
        end
    end
    size(K)
    N

    gauss = K / sum(K(:));
end

repo: pTVreg