text:
function r8vec_sorted_unique_count_test ( )

%*****************************************************************************80
%
%% R8VEC_SORTED_UNIQUE_COUNT_TEST tests R8VEC_SORTED_UNIQUE_COUNT;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2005
%
%  Author:
%
%    John Burkardt
%
  n = 30;
  tol = 0.25;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_SORTED_UNIQUE_COUNT_TEST\n' );
  fprintf ( 1, '  R8VEC_SORTED_UNIQUE_COUNT counts the unique entries\n' );
  fprintf ( 1, '  of a sorted R8VEC;\n' );

  b = 0.0;
  c = n;
  seed = 123456789;

  [ a, seed ] = r8vec_uniform_ab ( n, b, c, seed );

  a(1:n) = floor ( a(1:n) );

  unique_num = r8vec_sorted_unique_count ( n, a, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using a tolerance of %f\n', tol );
  fprintf ( 1, ...
    '  R8VEC_SORTED_UNIQUE_COUNT counts %d unique entries in A.\n', ...
    unique_num );

  return
end


text:
function [ p, l, u ] = gfpp_plu ( n, alpha )

%*****************************************************************************80
%
%% GFPP_PLU returns the PLU factorization of the GFPP matrix.
%
%  Discussion
%
%    This factorization assumes that Gaussian elimination is performed
%    without pivoting.  If ALPHA is not between 0 and 1, then the
%    PLU factors returned here will not be the PLU factors derived
%    from Gaussian elimination with pivoting.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA, determines subdiagonal elements.
%
%    Output, real P(N,N), L(N,N), U(N,N), the P, L, U factors
%    of the matrix.
%
  p = zeros ( n, n );
  for i = 1 : n
    p(i,i) = 1.0;
  end

  l = zeros ( n, n );
  for i = 1 : n
    l(i,1:i-1) = - abs ( alpha );
    l(i,i) = 1.0;
    l(i,i+1:n) = 0.0;
  end

  u = zeros ( n, n );
  u(1:n,1:n) = 0.0;
  for i = 1 : n
    u(i,i) = 1.0;
  end

  u(1,n) = 1.0;
  for i = 2 : n
    u(i,n) = 1.0 + abs ( alpha ) * sum ( u(1:i-1,n) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y=scInv(z, LinCol)
%scInv : Inverts a given impedance (admitance) by mirroring it about the origin
%
%  SYNOPSIS:
%     This function inverts the given impedance (admittance) by mirroring it about
%     the origin of the smith chart.
%
%     See also scDraw, scMove, scConCirc, scMatchCirc
%
%  SYNTAX:
%     [y] = scInv(z)
%
%  INPUT ARGUMENTS:
%     z     : Impedance (or admittance)
%
%  OUTPUT ARGUMENT:
%     y : Inverted admittance (or impedance)
%
%  EXAMPLE:
%        y = scInv([2 3])
%        y =
%        0.1538   -0.2308
%
%
%     Mohammad Ashfaq - (31-05-2000)
%     Mohammad Ashfaq - (13-04-2006) Modified (example included)
%

 if nargin < 1
    error('scInv.m: One input argument required...You may give LinCol as second argument as well');
 end

 if nargin == 1
    LinCol = 'm';
 end

 if (size(z)==[1, 1])
    if conj(z) == z
       r = z;
       x = 0;
    else
       r = real (z);
       x = imag (z);
    end
 elseif max(size(z)) == 2 && min(size(z)) == 1
    r = z(1);
    x = z(2);
 else
    error('scInv.m: Input must either be a complex number or a 1x2 vector  [r x]');
 end

 if (r+j*x)~=0
    y1   = 1/(r+j*x);
    y(1) = real(y1);
    y(2) = imag(y1);

 else
    y = inf;
 end

author: Sable

repo: mcbench-benchmarks


text:
function S = mixgetcovariance(mix,inds,centers)

if ~exist('inds','var'),
  inds = 1:mix.nin;
end;
if ~exist('centers','var'),
  centers = 1:mix.ncentres;
end;

inds = inds(:)';
ninds = length(inds);
centers = centers(:)';
ncenters = length(centers);
S = zeros(ninds,ninds,ncenters);

switch mix.covar_type,

 case 'spherical',
  for i = 1:ncenters,
    j = centers(i);
    S(:,:,i) = eye(ninds)*mix.covars(j);
  end;

 case 'diag',
  for i = 1:ncenters,
    j = centers(i);
    S(:,:,i) = diag(mix.covars(j,:));
  end;

 case {'full','sparse'},
  S = mix.covars(inds,inds,centers);

 case 'ppca',
  for i = 1:ncenters,
    j = centers(i);
    covars = mix.covars(j) * eye(mix.nin) + ...
	     mix.U(:, :, j)*(diag(mix.lambda(j, :))-diag(mix.covars(j)))* ...
	     (mix.U(:, :, j)');
    S(:,:,i) = covars(inds,inds);
  end;

 case 'zero',
  S(:) = 0;

 otherwise,
  error(['Unknown covar_type ',mix.covar_type]);
end;

author: kristinbranson

repo: JAABA


text:
function ksvddenoisedemo
%KSVDDENOISEDEMO K-SVD denoising demonstration.
%  KSVDDENISEDEMO reads an image, adds random white noise and denoises it
%  using K-SVD denoising. The input and output PSNR are compared, and the
%  trained dictionary is displayed.
%
%  To run the demo, type KSVDDENISEDEMO from the Matlab prompt.
%
%  See also KSVDDEMO.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  August 2009

disp(' ');
disp('  **********  K-SVD Denoising Demo  **********');
disp(' ');
disp('  This demo reads an image, adds random Gaussian noise, and denoises the image');
disp('  using K-SVD denoising. The function displays the original, noisy, and denoised');
disp('  images, and shows the resulting trained dictionary.');
disp(' ');

%% prompt user for image %%

pathstr = fileparts(which('ksvddenoisedemo'));
dirname = fullfile(pathstr, 'images', '*.png');
imglist = dir(dirname);

disp('  Available test images:');
disp(' ');
for k = 1:length(imglist)
  printf('  %d. %s', k, imglist(k).name);
end
disp(' ');

imnum = 0;
while (~isnumeric(imnum) || ~iswhole(imnum) || imnum<1 || imnum>length(imglist))
  imnum = input(sprintf('  Image to denoise (%d-%d): ', 1, length(imglist)), 's');
  imnum = sscanf(imnum, '%d');
end

imgname = fullfile(pathstr, 'images', imglist(imnum).name);

%% generate noisy image %%

sigma = 20;

disp(' ');
disp('Generating noisy image...');

im = imread(imgname);
im = double(im);

n = randn(size(im)) * sigma;
imnoise = im + n;

%% set parameters %%

params.x = imnoise;
params.blocksize = 8;
params.dictsize = 256;
params.sigma = sigma;
params.maxval = 255;
params.trainnum = 40000;
params.iternum = 20;
params.memusage = 'high';

% denoise!
disp('Performing K-SVD denoising...');
[imout, dict] = ksvddenoise(params);

% show results %

dictimg = showdict(dict,[1 1]*params.blocksize,round(sqrt(params.dictsize)),round(sqrt(params.dictsize)),'lines','highcontrast');
figure; imshow(imresize(dictimg,2,'nearest'));
title('Trained dictionary');

figure; imshow(im/params.maxval);
title('Original image');

figure; imshow(imnoise/params.maxval);
title(sprintf('Noisy image, PSNR = %.2fdB', 20*log10(params.maxval * sqrt(numel(im)) / norm(im(:)-imnoise(:))) ));

figure; imshow(imout/params.maxval);
title(sprintf('Denoised image, PSNR: %.2fdB', 20*log10(params.maxval * sqrt(numel(im)) / norm(im(:)-imout(:))) ));

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [faces,idx,facemap]=uniqfaces(elem)
%
% [faces,idx,facemap]=uniqfaces(elem)
%
% return the unique face list from a or tetrahedral mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%     elem: a list of elements, each row is a list of nodes for an element.
%           elem can have 2, 3 or 4 columns
%
% output:
%     face: unique faces in the mesh, denoted by a triplet of node indices
%     idx:  index of the output in the raw face list (returned by meshface)
%     facemap: index of the raw faces in the output list (for triangular mesh)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(elem)==3)
   faces=elem;
elseif(size(elem)>=4)
   faces=meshface(elem);
else
   error('invalid input');
end

[ufaces,idx,jdx]=unique(sort(faces,2),'rows');
faces=faces(idx,:);
if(nargout>2)
    facemap=reshape(jdx,[size(elem,1) nchoosek(size(elem,2),3)]);
end

author: fangq


text:
function [ patches, pdm, clmParams ] = Load_CLNF_wild()
%LOAD_CLNF_GENERAL Summary of this function goes here
%   Detailed explanation goes here
    clmParams = struct;

    clmParams.window_size = [25,25; 23,23; 21,21;21,21];
    clmParams.numPatchIters = size(clmParams.window_size,1);

    [patches] = Load_Patch_Experts( '../models/wild/', 'ccnf_patches_*_wild.mat', [], [], clmParams);

    % the default PDM to use
    pdmLoc = ['../models/pdm/pdm_68_aligned_wild.mat'];
    load(pdmLoc);

    pdm = struct;
    pdm.M = double(M);
    pdm.E = double(E);
    pdm.V = double(V);

    % the default model parameters to use
    clmParams.regFactor = [35, 27, 20, 20];
    clmParams.sigmaMeanShift = [1.25, 1.375, 1.5, 1.5];
    clmParams.tikhonov_factor = [2.5, 5, 7.5, 7.5];

    clmParams.startScale = 1;
    clmParams.num_RLMS_iter = 10;
    clmParams.fTol = 0.01;
    clmParams.useMultiScale = true;
    clmParams.use_multi_modal = 1;
    clmParams.multi_modal_types  = patches(1).multi_modal_types;

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function pass = test_constructor( pref )
% Test the Chebfun2v constructor when performing simple arithmetic
% operations.

if ( nargin < 1 )
    pref = chebfunpref;
end

tol = 100*pref.cheb2Prefs.chebfun2eps;

D = [-1 1 -1 1; -2 3 0 1];

for kk = 1:size(D,1)
    d = D(kk,:);
    % Check the constructor works with lots of different syntax:
    f1 = @(x,y) cos(x.*y);
    f2 = @(x,y) sin(x) + cos(y);
    f3 = @(x,y) exp(x) - exp(-y) + x;
    g1 = chebfun2(f1, d);
    g2 = chebfun2(f2, d);
    g3 = chebfun2(f3, d);

    H1 = chebfun2v(f1, f2, d);
    H2 = chebfun2v(f1, g2, d);
    H3 = chebfun2v(g1, f2, d);
    H4 = chebfun2v(g1, g2, d);

    pass(1) = norm( H1 - H2 ) < tol;
    pass(2) = norm( H1 - H3 ) < tol;
    pass(3) = norm( H1 - H4 ) < tol;

    % Now try the same thing for three components:
    H1 = chebfun2v(f1, f2, f3, d);
    H2 = chebfun2v(f1, f2, g3, d);
    H3 = chebfun2v(f1, g2, f3, d);
    H4 = chebfun2v(g1, f2, f3, d);
    H5 = chebfun2v(g1, g2, f3, d);
    H6 = chebfun2v(g1, f2, g3, d);
    H7 = chebfun2v(f1, g2, g3, d);
    H8 = chebfun2v(g1, g2, g3, d);

    pass(4) = norm( H1 - H2 ) < tol;
    pass(5) = norm( H1 - H3 ) < tol;
    pass(6) = norm( H1 - H4 ) < tol;
    pass(7) = norm( H1 - H5 ) < tol;
    pass(8) = norm( H1 - H6 ) < tol;
    pass(9) = norm( H1 - H7 ) < tol;
    pass(10) = norm( H1 - H8 ) < tol;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function d = dpoly(obj,feat,varargin)
% INPUTS:
% obj contains the mesh points and the outer, mainland, inner polygons,
% bbox (the bounding box)
% feat contains inpoly_flip to check whether to flip the inpoly result or not
%
% OUTPUTS:
% d is the distance from point, p to closest point on polygon
% (d is  negative if inside the bounded polygon, pv and positive if outside)
% by Keith Roberts and William Pringle 2017-2018.
if nargin == 3
    p  = varargin{1};
else
    [xg,yg] = CreateStructGrid(obj);
    p = [xg(:),yg(:)];
    clearvars xg yg
end
%% Doing the distance + inpoly calc
if ~isempty(obj.lmsl)
    pv = [obj.lmsl.mainland; obj.lmsl.inner];
    pg = [obj.lmsl.outer; obj.lmsl.inner];
else
    pv = [feat.mainland; feat.inner];
    pg = [feat.outer; feat.inner];
end
edges = Get_poly_edges( pg );
pv1 = pv; % <-dup for inpoly to work
pv1(isnan(pv(:,1)),:) = []; clear pv;
d = 0*p(:,1); in = false(size(d));
noblks = ceil(length(p)*2*8*1e-9);
blklen = floor(length(p)/noblks);
ns = 1;
disp(['Memory management: number of loops in dpoly = ' num2str(noblks)])
tic
for blks = 1:noblks
    if blks == noblks
        ne = length(p);
    else
        ne = ns + blklen - 1;
    end
    %% Doing the distance check
    [~,d(ns:ne)] = WrapperForKsearch(pv1, p(ns:ne,:), 1);

    %% Doing the inpoly check
    in(ns:ne) = inpoly(p(ns:ne,:), pg, edges);

    ns = ne + 1;
end
toc

% d is negative if inside polygon and vice versa.
if feat.inpoly_flip
    d = (-1).^(~in).*d;
else
    d = (-1).^( in).*d;
end
end

author: CHLNDDEV


text:
function hf = r8_hyper_2f1 ( a, b, c, x )

%*****************************************************************************80
%
%% R8_HYPER_2F1 evaluates the hypergeometric function F(A,B,C,X).
%
%  Discussion:
%
%    A minor bug was corrected.  The HW variable, used in several places as
%    the "old" value of a quantity being iteratively improved, was not
%    being initialized.  JVB, 11 February 2008.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
%    MATLAB version by John Burkardt.
%
%    The F77 original version of this routine is copyrighted by
%    Shanjie Zhang and Jianming Jin.  However, they give permission to
%    incorporate this routine into a user program provided that the copyright
%    is acknowledged.
%
%  Reference:
%
%    Shanjie Zhang, Jianming Jin,
%    Computation of Special Functions,
%    Wiley, 1996,
%    ISBN: 0-471-11963-6,
%    LC: QA351.C45
%
%  Parameters:
%
%    Input, real A, B, C, X, the arguments of the function.
%    C must not be equal to a nonpositive integer.
%    X < 1.
%
%    Output, real HF, the value of the function.
%
  el = 0.5772156649015329;

  l0 = ( c == floor ( c ) ) && ( c < 0.0 );
  l1 = ( 1.0 - x < 1.0E-15 ) && ( c - a - b <= 0.0 );
  l2 = ( a == floor ( a ) ) && ( a < 0.0 );
  l3 = ( b == floor ( b ) ) && ( b < 0.0 );
  l4 = ( c - a == floor ( c - a ) ) && ( c - a <= 0.0 );
  l5 = ( c - b == floor ( c - b ) ) && ( c - b <= 0.0 );

  if ( l0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_HYPER_2F1 - Fatal error!\n' );
    fprintf ( 1, '  The hypergeometric series is divergent.\n' );
    fprintf ( 1, '  C is integral and negative.\n' );
    fprintf ( 1, '  C = %f\n', c );
    error ( 'R8_HYPER_F1 - Fatal error!' );
  end

  if ( l1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_HYPER_2F1 - Fatal error!\n' );
    fprintf ( 1, '  The hypergeometric series is divergent.\n' );
    fprintf ( 1, '  1 = X < 0, C - A - B <= 0.\n' );
    fprintf ( 1, '  A = %f\n', a );
    fprintf ( 1, '  B = %f\n', b );
    fprintf ( 1, '  C = %f\n', c );
    fprintf ( 1, '  X = %f\n', x );
    error ( 'R8_HYPER_F1 - Fatal error!' );
  end

  if ( 0.95 < x )
    eps = 1.0E-08;
  else
    eps = 1.0E-15;
  end

  if ( x == 0.0 || a == 0.0 || b == 0.0 )

    hf = 1.0;
    return

  elseif ( 1.0 - x == eps && 0.0 < c - a - b )

    gc = gamma ( c );
    gcab = gamma ( c - a - b );
    gca = gamma ( c - a );
    gcb = gamma ( c - b );
    hf = gc * gcab / ( gca * gcb );
    return

  elseif ( 1.0 + x <= eps && abs ( c - a + b - 1.0 ) <= eps )

    g0 = sqrt ( pi ) * 2.0^( - a );
    g1 = gamma ( c );
    g2 = gamma ( 1.0 + a / 2.0 - b );
    g3 = gamma ( 0.5 + 0.5 * a );
    hf = g0 * g1 / ( g2 * g3 );
    return

  elseif ( l2 || l3 )

    if ( l2 )
      nm = floor ( abs ( a ) );
    end

    if ( l3 )
      nm = floor ( abs ( b ) );
    end

    hf = 1.0;
    r = 1.0;

    for k = 1 : nm
      r = r * ( a + k - 1.0 ) * ( b + k - 1.0 ) ...
        / ( k * ( c + k - 1.0 ) ) * x;
      hf = hf + r;
    end

    return

  elseif ( l4 || l5 )

    if ( l4 )
      nm = floor ( abs ( c - a ) );
    end

    if ( l5 )
      nm = floor ( abs ( c - b ) );
    end

    hf = 1.0;
    r  = 1.0;
    for k = 1 : nm
      r = r * ( c - a + k - 1.0 ) * ( c - b + k - 1.0 ) ...
        / ( k * ( c + k - 1.0 ) ) * x;
      hf = hf + r;
    end
    hf = ( 1.0 - x )^( c - a - b ) * hf;
    return

  end

  aa = a;
  bb = b;
  x1 = x;

  if ( x < 0.0 )
    x = x / ( x - 1.0 );
    if ( a < c && b < a && 0.0 < b )
      a = bb;
      b = aa;
    end
    b = c - b;
  end

  if ( 0.75 <= x )

    gm = 0.0;

    if ( abs ( c - a - b - floor ( c - a - b ) ) < 1.0E-15 )

      m = floor ( c - a - b );
      ga = gamma ( a );
      gb = gamma ( b );
      gc = gamma ( c );
      gam = gamma ( a + m );
      gbm = gamma ( b + m );

      pa = r8_psi ( a );
      pb = r8_psi ( b );

      if ( m ~= 0 )
        gm = 1.0;
      end

      for j = 1 : abs ( m ) - 1
        gm = gm * j;
      end

      rm = 1.0;
      for j = 1 : abs ( m )
        rm = rm * j;
      end

      f0 = 1.0;
      r0 = 1.0;
      r1 = 1.0;
      sp0 = 0.0;
      sp = 0.0;

      if ( 0 <= m )

        c0 = gm * gc / ( gam * gbm );
        c1 = - gc * ( x - 1.0 )^m / ( ga * gb * rm );

        for k = 1 : m - 1
          r0 = r0 * ( a + k - 1.0 ) * ( b + k - 1.0 ) ...
            / ( k * ( k - m ) ) * ( 1.0 - x );
          f0 = f0 + r0;
        end

        for k = 1 : m
          sp0 = sp0 + 1.0 / ( a + k - 1.0 ) ...
            + 1.0 / ( b + k - 1.0 ) - 1.0 / k;
        end

        f1 = pa + pb + sp0 + 2.0 * el + log ( 1.0 - x );
        hw = f1;

        for k = 1 : 250

          sp = sp + ( 1.0 - a ) / ( k * ( a + k - 1.0 ) ) ...
            + ( 1.0 - b ) / ( k * ( b + k - 1.0 ) );

          sm = 0.0;
          for j = 1 : m
            sm = sm + ( 1.0 - a ) ...
              / ( ( j + k ) * ( a + j + k - 1.0 ) ) ...
              + 1.0 / ( b + j + k - 1.0 );
          end

          rp = pa + pb + 2.0 * el + sp + sm + log ( 1.0 - x );

          r1 = r1 * ( a + m + k - 1.0 ) * ( b + m + k - 1.0 ) ...
            / ( k * ( m + k ) ) * ( 1.0 - x );

          f1 = f1 + r1 * rp;

          if ( abs ( f1 - hw ) < abs ( f1 ) * eps )
            break
          end

          hw = f1;

        end

        hf = f0 * c0 + f1 * c1;

      elseif ( m < 0 )

        m = - m;
        c0 = gm * gc / ( ga * gb * ( 1.0 - x )^m );
        c1 = - ( - 1 )^m * gc / ( gam * gbm * rm );

        for k = 1 : m - 1
          r0 = r0 * ( a - m + k - 1.0 ) * ( b - m + k - 1.0 ) ...
            / ( k * ( k - m ) ) * ( 1.0 - x );
          f0 = f0 + r0;
        end

        for k = 1 : m
          sp0 = sp0 + 1.0 / k;
        end

        f1 = pa + pb - sp0 + 2.0 * el + log ( 1.0 - x );
        hw = f1;

        for k = 1 : 250

          sp = sp + ( 1.0 - a ) ...
            / ( k * ( a + k - 1.0 ) ) ...
            + ( 1.0 - b ) / ( k * ( b + k - 1.0 ) );

          sm = 0.0;
          for j = 1 : m
            sm = sm + 1.0 / ( j + k );
          end

          rp = pa + pb + 2.0 * el + sp - sm + log ( 1.0 - x );

          r1 = r1 * ( a + k - 1.0 ) * ( b + k - 1.0 ) ...
            / ( k * ( m + k ) ) * ( 1.0 - x );

          f1 = f1 + r1 * rp;

          if ( abs ( f1 - hw ) < abs ( f1 ) * eps )
            break
          end

          hw = f1;

        end

        hf = f0 * c0 + f1 * c1;

      end

    else

      ga = gamma ( a );
      gb = gamma ( b );
      gc = gamma ( c );
      gca = gamma ( c - a );
      gcb = gamma ( c - b );
      gcab = gamma ( c - a - b );
      gabc = gamma ( a + b - c );
      c0 = gc * gcab / ( gca * gcb );
      c1 = gc * gabc / ( ga * gb ) * ( 1.0 - x )^( c - a - b );
      hf = 0.0;
      hw = hf;
      r0 = c0;
      r1 = c1;

      for k = 1 : 250

        r0 = r0 * ( a + k - 1.0 ) * ( b + k - 1.0 ) ...
          / ( k * ( a + b - c + k ) ) * ( 1.0 - x );

        r1 = r1 * ( c - a + k - 1.0 ) * ( c - b + k - 1.0 ) ...
          / ( k * ( c - a - b + k ) ) * ( 1.0 - x );

        hf = hf + r0 + r1;

        if ( abs ( hf - hw ) < abs ( hf ) * eps )
          break
        end

        hw = hf;

      end

      hf = hf + c0 + c1;

    end

  else

    a0 = 1.0;

    if ( a < c && c < 2.0 * a && b < c && c < 2.0 * b )

      a0 = ( 1.0 - x )^( c - a - b );
      a = c - a;
      b = c - b;

    end

    hf = 1.0;
    hw = hf;
    r = 1.0;

    for k = 1 : 250

      r = r * ( a + k - 1.0 ) * ( b + k - 1.0 ) ...
        / ( k * ( c + k - 1.0 ) ) * x;

      hf = hf + r;

      if ( abs ( hf - hw ) <= abs ( hf ) * eps )
        break
      end

      hw = hf;

    end

    hf = a0 * hf;

  end

  if ( x1 < 0.0 )
    x = x1;
    c0 = 1.0 / ( 1.0 - x )^aa;
    hf = c0 * hf;
  end

  if ( 120 < k )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_HYPER_2F1 - Warning!\n' );
    fprintf ( 1, '  A large number of iterations were needed.\n' );
    fprintf ( 1, '  The accuracy of the results should be checked.\n' );
  end

  return
end


text:
function distance = cylinder_point_dist_3d ( p1, p2, r, p )

%*****************************************************************************80
%
%% CYLINDER_POINT_DIST_3D determines the distance from a cylinder to a point in 3D.
%
%  Discussion:
%
%    We are computing the distance to the SURFACE of the cylinder.
%
%    The surface of a (right) (finite) cylinder in 3D is defined by an axis,
%    which is the line segment from point P1 to P2, and a radius R.  The points
%    on the surface of the cylinder are:
%    * points at a distance R from the line through P1 and P2, and whose nearest
%      point on the line through P1 and P2 is strictly between P1 and P2,
%    PLUS
%    * points at a distance less than or equal to R from the line through P1
%      and P2, whose nearest point on the line through P1 and P2 is either
%      P1 or P2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3), P2(3), the first and last points
%    on the axis line of the cylinder.
%
%    Input, real R, the radius of the cylinder.
%
%    Input, real P(3), the point.
%
%    Output, real DISTANCE, the distance from the point
%    to the cylinder.
%
  dim_num = 3;

  axis(1:dim_num) = p2(1:dim_num) - p1(1:dim_num);

  axis_length = r8vec_norm ( dim_num, axis );

  if ( axis_length == 0.0 )
    distance = - r8_huge ( );
    return
  end

  axis(1:dim_num) = axis(1:dim_num) / axis_length;

  p_dot_axis = ( p(1:dim_num) - p1(1:dim_num) ) * axis(1:dim_num)';
%
%  Case 1: Below bottom cap.
%
  if ( p_dot_axis <= 0.0 )

    distance = disk_point_dist_3d ( p1, r, axis, p );
%
%  Case 2: between cylinder planes.
%
  elseif ( p_dot_axis <= axis_length )

    p_length = r8vec_norm ( dim_num, p(1:dim_num) - p1(1:dim_num) );
    off_axis_component = sqrt ( p_length.^2 - p_dot_axis.^2 );

    distance = abs ( off_axis_component - r );

    if ( off_axis_component < r )
      distance = min ( distance, axis_length - p_dot_axis );
      distance = min ( distance, p_dot_axis );
    end
%
%  Case 3: Above the top cap.
%
  elseif ( axis_length < p_dot_axis )

    distance = disk_point_dist_3d ( p2, r, axis, p );

  end

  return
end


text:
function [event] = PL2EventTs(filename, channel)
% PL2EventTs(filename, channel): read event timestamps from a .pl2 file
%
% [event] = PL2EventTs(filename, channel)
% [event] = PL2EventTs(filename, 1)
% [event] = PL2EventTs(filename, 'KBD2')
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   channel - event channel name or 1-based channel number
%           1-based channel number is the index of the event channel in pl2.EventChannels cell array
%           shown in the Name column of the output produced by PL2Print(pl2.EventChannels)
%           To print the list of event channels, use these commands:
%
%           pl2 = PL2GetFileIndex(filename);
%           PL2Print(pl2.EventChannels);
%
% OUTPUT:
%   event.Ts - array of timestamps (in seconds)
%   event.Strobed - array of strobed event values

event = internalPL2EventTs(filename, channel);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% This demo shows examples of data-driven reversible jump (RJ) moves
%  performed on a toy dataset with Gaussian emissions (MANY dims per step)
% We compare and constrast  state-sequence based moves
%   and continuous parameter based moves

clear all;
close all;

doAR = 1;
doDataDriven = 1;
obsDim = 25;

% Generate Gaussian toy data for the BP-HMM
Ktrue = 8;
T = 500;
if doAR
dataP = {'SynthAR', 'nStates', Ktrue, 'obsDim', obsDim, 'T', T};
else
dataP = {'SynthGaussian', 'nStates', Ktrue, 'obsDim', obsDim, 'T', T};
end
mP = {};
initP = {'F.nTotal', 1};
aalgP = {'Niter', 30, ...
         'RJ.birthPropDistr', 'DataDriven', 'BP.doSampleMass', 0, ...
         'BP.doSampleConc', 0, 'HMM.doSampleHypers', 0, ...
          'RJ.minW', obsDim};

for doZMove = [0 1]
    algP = aalgP;
    algP(end+1:end+2) = {'doSampleUniqueZ', doZMove};
    algP(end+1:end+2) = {'doSampleFUnique', ~doZMove};

    outP  = {obsDim+(doZMove+1)*100,1};

    runBPHMM( dataP, mP, outP, algP, initP );
end

plotLogPr( obsDim+(1:2)*100, 1, {'Unique C', 'Unique Z'} );

author: michaelchughes

repo: NPBayesHMM


text:
function p = real(p);
%REAL         Real part of (interval) polynomial
%
%   r = real(p)
%

% written  08/28/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  p.c = real(p.c);
  p = normalize(p);

repo: OpenMAS


text:
% std_plotcurve() - plot ERP or spectral traces for a STUDY component
%                   or channel cluster
% Usage:
%          >> std_plotcurve( axvals, data, 'key', 'val', ...)
% Inputs:
%  axvals - [vector or cell array] axis values for the data.
%  data  -  [cell array] mean data for each subject group and/or data
%           condition. For example, to plot mean ERPs from a STUDY
%           for epochs of 800 frames in two conditions from three groups
%           of 12 subjects:
%
%           >> data = { [800x12] [800x12] [800x12];... % 3 groups, cond 1
%                       [800x12] [800x12] [800x12] };  % 3 groups, cond 2
%           >> std_plotcurve(erp_ms,data);
%
%           By default, parametric statistics are computed across subjects
%           in the three groups. (group,condition) ERP averages are plotted.
%           See below and >> help statcond
%           for more information about the statistical computations.
%
% Optional display parameters:
%  'datatype'    - ['erp'|'spec'] data type {default: 'erp'}
%  'titles'      - [cell array of string] titles for each of the subplots.
%                  { default: none}
%
% Statistics options:
%  'groupstats'  - [cell] One p-value array per group {default: {}}
%  'condstats'   - [cell] One p-value array per condition {default: {}}
%  'interstats'  - [cell] Interaction p-value arrays {default: {}}
%  'threshold'   - [NaN|real<<1] Significance threshold. NaN -> plot the
%                  p-values themselves on a different figure. When possible,
%                  significance regions are indicated below the data.
%                  {default: NaN}
%
% Curve plotting options (ERP and spectrum):
%  'plotgroups'  - ['together'|'apart'] 'together' -> plot mean results
%                  for subject groups in the same figure panel in different
%                  colors. 'apart' -> plot group results on different figure
%                  panels {default: 'apart'}
%  'plotconditions' - ['together'|'apart'] 'together' -> plot mean results
%                  for data conditions on the same figure panel in
%                  different
%                  colors. 'apart' -> plot conditions on different figure
%                  panel. Note: 'plotgroups' and 'plotconditions' arguments
%                  cannot both be 'together' {default: 'apart'}
%  'legend'      - ['on'|'off'] turn plot legend on/off {default: 'off'}
%  'plotdiff'    - ['on'|'off'] plot difference between two groups
%                  or conditions plotted together.
%  'plotstderr'  - ['on'|'off'|'diff'|'nocurve'|'diffnocurve'] plots in
%                  a surface indicating the standard error. 'diff' only
%                  does it for the difference (requires 'plotdiff' 'on'
%                  above). 'nocurve' does not plot the mean. This functionality
%                  does not work for all data configuration {default: 'off'}
%  'figure'      - ['on'|'off'] creates a new figure ('on'). The 'off' mode
%                  plots all of the groups and conditions on the same pannel.
% 'plotsubjects' - ['on'|'off'] overplot traces for individual components
%                  or channels {default: 'off'}
% 'singlesubject' - ['on'|'off'] set to 'on' to plot single subject.
%                  {default: 'off'}
% 'ylim'         - [min max] ordinate limits for ERP and spectrum plots
%                  {default: all available data}
%
% Scalp map plotting options:
%  'chanlocs'    - [struct] channel locations structure
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: pop_erspparams(), pop_erpparams(), pop_specparams(), statcond()

function std_plotcurve(allx, data, varargin)

pgroup = [];
pcond  = [];
pinter = [];
if nargin < 2
    help std_plotcurve;
    return;
end;

opt = finputcheck( varargin, { 'ylim'          'real'   []              [];
                               'filter'        'real'   []              [];
                               'threshold'     'real'   []              NaN;
                               'unitx'         'string' []              'ms';
                               'chanlocs'      'struct' []              struct('labels', {});
                               'plotsubjects'  'string' { 'on' 'off' }  'off';
                               'condnames'     'cell'   []              {}; % just for legends
                               'groupnames'    'cell'   []              {}; % just for legends
                               'groupstats'    'cell'   []              {};
                               'condstats'     'cell'   []              {};
                               'interstats'    'cell'   []              {};
                               'titles'        'cell'   []              {};
                               'figure'        'string' { 'on' 'off' }   'on';
                               'plottopo'      'string' { 'on' 'off' }   'off';
                               'plotstderr'    'string' { 'on' 'off' 'diff' 'nocurve' }   'off';
                               'plotdiff'      'string' { 'on' 'off' }   'off';
                               'legend'        { 'string' 'cell' } { { 'on' 'off' } {} }  'off';
                               'datatype'      'string' { 'ersp' 'itc' 'erp' 'spec' }    'erp';
                               'plotgroups'    'string' { 'together' 'apart' }  'apart';
                               'plotmode'      'string' { 'test' 'condensed' }  'test'; % deprecated
                               'plotconditions'    'string' { 'together' 'apart' }  'apart' }, 'std_plotcurve');

% opt.figure =  'off'; % test by nima
if isstr(opt), error(opt); end;
opt.singlesubject = 'off';
if strcmpi(opt.plottopo, 'on') && size(data{1},3) == 1, opt.singlesubject = 'on'; end;
%if size(data{1},2) == 1,                              opt.singlesubject = 'on'; end;
if all(all(cellfun('size', data, 2)==1))               opt.singlesubject = 'on'; end;
if any(any(cellfun('size', data, 2)==1)), opt.groupstats = {}; opt.condstats = {}; end;
if strcmpi(opt.datatype, 'spec'), opt.unit = 'Hz'; end;
if strcmpi(opt.plotsubjects, 'on')
    opt.plotgroups      = 'apart';
    opt.plotconditions  = 'apart';
end;
if strcmpi(opt.plotconditions, 'together') &&  ~isempty(opt.groupstats), opt.plotconditions = 'apart'; end;
if strcmpi(opt.plotgroups,     'together') &&  ~isempty(opt.condstats) , opt.plotgroups     = 'apart'; end;
if isstr(opt.legend), opt.legend = {}; end;
if isempty(opt.titles), opt.titles = cell(10,10); opt.titles(:) = { '' }; end;
if length(data(:)) == length(opt.legend(:)),
    opt.legend = reshape(opt.legend, size(data))';
    opt.legend(cellfun(@isempty, data)) = [];
    opt.legend = (opt.legend)';
end;

% plot
% ----
% if strcmpi(opt.figure, 'on'),
% else
%     % all groups and conditions in the same figureopt.chanlocs
%     ncol = size(tmpdata,3);
%     if ncol == 1, ncol = size(tmpdata,1); end;
%     tmpcol = col(colcount:colcount+ncol-1);
%     colcount = mod(colcount+ncol-1, length(col))+1;
%     %if strcmpi(opt.plotgroups, 'together') && isempty(tmpdata{1})
%
% end;
%

% color matrix
% -----------------------
onecol  = { 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' };
manycol = { 'b' 'g' 'm' 'c' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' ...
                   'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' };
if strcmpi(opt.plotgroups, 'together') || strcmpi(opt.plotconditions, 'together') || strcmpi(opt.figure, 'off')
     col = manycol;
else col = onecol;
end;
nonemptycell = find(~cellfun(@isempty, data));
maxdim = max(length(data(:)), size(data{nonemptycell(1)}, ndims(data{nonemptycell(1)})));
tmpcol = col;
if strcmpi(opt.plotsubjects, 'off')
    coldata = col(mod([0:maxdim-1], length(col))+1);
    coldata = reshape(coldata(1:length(data(:))), size(data));
else
    coldata = cell(size(data));
end;

% remove empty entries
% --------------------
datapresent = ~cellfun(@isempty, data);
if size(data,1) > 1, for c = size(data,1):-1:1, if sum(datapresent(c,:)) == 0, data(c,:) = []; coldata(c,:) = []; if ~strcmpi(opt.plotconditions, 'together') opt.titles(c,:) = []; end; if ~isempty(opt.groupstats), opt.groupstats(c) = []; end; end; end; end;
if size(data,2) > 1, for g = size(data,2):-1:1, if sum(datapresent(:,g)) == 0, data(:,g) = []; coldata(:,g) = []; if ~strcmpi(opt.plotgroups    , 'together') opt.titles(:,g) = []; end; if ~isempty(opt.condstats ), opt.condstats( g) = []; end; end; end; end;
if strcmpi(opt.plotsubjects, 'off'), tmpcol = coldata'; tmpcol = tmpcol(:)'; end;

% number of columns and rows to plot
% ----------------------------------
nc = size(data,1);
ng = size(data,2);
if strcmpi(opt.plotgroups, 'together'),      ngplot = 1; else ngplot = ng; end;
if strcmpi(opt.plotconditions,  'together'), ncplot = 1; else ncplot = nc; end;
if nc >= ng, opt.subplot = 'transpose';
else         opt.subplot = 'normal';
end;
if isempty(opt.condnames)
    for c=1:nc, opt.condnames{c} = sprintf('Cond. %d', c); end;
    if nc == 1, opt.condnames = { '' }; end;
end;
if isempty(opt.groupnames)
    for g=1:ng, opt.groupnames{g} = sprintf('Group. %d', g); end;
    if ng == 1, opt.groupnames = { '' }; end;
end;

% plotting paramters
% ------------------
if ng > 1 && ~isempty(opt.groupstats), addc = 1; else addc = 0; end;
if nc > 1 && ~isempty(opt.condstats ), addr = 1; else addr = 0; end;
if strcmpi(opt.singlesubject, 'off') ...
        && ( ~isempty(opt.condstats) || ~isempty(opt.groupstats) ) % only for curves
    plottag = 0;
    if strcmpi(opt.plotgroups, 'together') && isempty(opt.condstats) && ~isempty(opt.groupstats) && ~isnan(opt.threshold), addc = 0; plottag = 1; end;
    if strcmpi(opt.plotconditions , 'together') && ~isempty(opt.condstats) && isempty(opt.groupstats) && ~isnan(opt.threshold), addr = 0; plottag = 1; end;
    if ~isnan(opt.threshold) && plottag == 0 && strcmpi(opt.figure, 'on')
        disp('Warning: cannot plot condition/group on the same panel while using a fixed');
        disp('         threshold, unless you only compute statistics for ether groups or conditions');
        opt.plotgroups = 'apart';
        opt.plotconditions  = 'apart';
    end;
end;

% compute significance mask
% --------------------------
if ~isempty(opt.interstats), pinter = opt.interstats{3}; end;

if ~isnan(opt.threshold) && ( ~isempty(opt.groupstats) || ~isempty(opt.condstats) )
    pcondplot  = opt.condstats;
    pgroupplot = opt.groupstats;
    pinterplot = pinter;
    maxplot = 1;
else
    warning off;
    for ind = 1:length(opt.condstats),  pcondplot{ind}  = -log10(opt.condstats{ind}); end;
    for ind = 1:length(opt.groupstats), pgroupplot{ind} = -log10(opt.groupstats{ind}); end;
    if ~isempty(pinter), pinterplot = -log10(pinter); end;
    maxplot = 3;
    warning on;
end;

% labels
% ------
if strcmpi(opt.unitx, 'ms'), xlab = 'Time (ms)';      ylab = 'Potential (\muV)';
else                         xlab = 'Frequency (Hz)'; ylab = 'Power (10*log_{10}(\muV^{2}/Hz))';
end;
if ~isnan(opt.threshold), statopt = {  'xlabel' xlab };
else                      statopt = { 'logpval' 'on' 'xlabel' xlab 'ylabel' '-log10(p)' 'ylim' [0 maxplot] };
end;

% adjust figure size
% ------------------
if strcmpi(opt.figure, 'on')
    figure('color', 'w');
    pos = get(gcf, 'position');
    basewinsize = 200/max(nc,ng)*3;
    if strcmpi(opt.plotgroups, 'together') pos(3) = 200*(1+addc);
    else                                   pos(3) = 200*(ng+addc);
    end;
    if strcmpi(opt.plotconditions , 'together') pos(4) = 200*(1+addr);
    else                                        pos(4) = 200*(nc+addr);
    end;
    if strcmpi(opt.subplot, 'transpose'), set(gcf, 'position', [ pos(1) pos(2) pos(4) pos(3)]);
    else                                  set(gcf, 'position', pos);
    end;
else
    opt.subplot = 'noplot';
end;

tmplim = [Inf -Inf];
colcount = 1; % only when plotting all conditions on the same figure
for c = 1:ncplot
    for g = 1:ngplot
        if strcmpi(opt.plotgroups, 'together'),         hdl(c,g)=mysubplot(ncplot+addr, ngplot+addc, 1 + (c-1)*(ngplot+addc), opt.subplot); ci = g;
        elseif strcmpi(opt.plotconditions, 'together'), hdl(c,g)=mysubplot(ncplot+addr, ngplot+addc, g, opt.subplot); ci = c;
        else                                            hdl(c,g)=mysubplot(ncplot+addr, ngplot+addc, g + (c-1)*(ngplot+addc), opt.subplot); ci = 1;
        end;

        if ~isempty(data{c,g})

            % read all data from one condition or group
            % -----------------------------------------
            dimreduced_sizediffers = 0;
            if ncplot ~= nc && ngplot ~= ng
                maxdim = max(max(cellfun(@(x)(size(x, ndims(x))), data)));
                for cc = 1:size(data,1)
                    for gg = 1:size(data,2)
                        tmptmpdata = real(data{cc,gg});
                        if ndims(tmptmpdata) == 3,
                            if cc == 1 && gg == 1, tmpdata = NaN*zeros([size(tmptmpdata,1) size(tmptmpdata,2) maxdim length(data(:))]); end;
                            tmpdata(:,:,1:size(tmptmpdata,3),gg+((cc-1)*ng)) = tmptmpdata;
                        else
                            if cc == 1 && gg == 1, tmpdata = NaN*zeros([size(tmptmpdata,1) maxdim length(data(:))]); end;
                            tmpdata(:,1:size(tmptmpdata,2),gg+((cc-1)*ng))   = tmptmpdata;
                        end;
                    end;
                end;
            elseif ncplot ~= nc % plot conditions together
                for ind = 2:size(data,1), if any(size(data{ind,1}) ~= size(data{1})), dimreduced_sizediffers = 1; end; end;
                for cc = 1:nc
                    tmptmpdata = real(data{cc,g});
                    if dimreduced_sizediffers
                        tmptmpdata = nan_mean(tmptmpdata,ndims(tmptmpdata));
                    end;
                    if cc == 1, tmpdata = zeros([size(tmptmpdata) nc]); end;
                    if ndims(tmptmpdata) == 3, tmpdata(:,:,:,cc) = tmptmpdata;
                    else                       tmpdata(:,:,cc)   = tmptmpdata;
                    end;
                end;
            elseif ngplot ~= ng % plot groups together
                for ind = 2:size(data,2), if any(size(data{1,ind}) ~= size(data{1})), dimreduced_sizediffers = 1; end; end;
                for gg = 1:ng
                    tmptmpdata = real(data{c,gg});
                    if dimreduced_sizediffers
                        tmptmpdata = nan_mean(tmptmpdata,ndims(tmptmpdata));
                    end;
                    if gg == 1, tmpdata = zeros([size(tmptmpdata) nc]); end;
                    if ndims(tmptmpdata) == 3, tmpdata(:,:,:,gg) = tmptmpdata;
                    else                       tmpdata(:,:,gg)   = tmptmpdata;
                    end;
                end;
            else tmpdata = real(data{c,g});
                % nothing
            end;

            % plot difference
            % ---------------
            if ~strcmpi(opt.plotdiff, 'off')
                if ngplot ~= ng || ncplot ~= nc
                    if size(tmpdata,3) == 2
                        tmpdata(:,:,end+1) = tmpdata(:,:,2)-tmpdata(:,:,1);
                        opt.legend{end+1} = [ opt.legend{2} '-' opt.legend{1} ];
                    elseif size(tmpdata,4) == 2
                        tmpdata(:,:,:,end+1) = tmpdata(:,:,:,2)-tmpdata(:,:,:,1);
                        opt.legend{end+1} = [ opt.legend{2} '-' opt.legend{1} ];
                    else
                        disp('Cannot plot difference, more than 2 indep. variable values');
                    end;
                else
                    disp('Cannot plot difference, indep. variable value must be plotted together');
                end;
            end;

            if ~isempty(opt.filter), tmpdata = myfilt(tmpdata, 1000/(allx(2)-allx(1)), 0, opt.filter); end;

            % plotting options
            % ----------------
            plotopt = { allx };
            if ~dimreduced_sizediffers
                if strcmpi(opt.plottopo, 'on'),
                    if strcmpi(opt.plotsubjects, 'off') tmpstd = squeeze(real(std(tmpdata,[],3)))/sqrt(size(tmpdata,3)); tmpstd = squeeze(permute(tmpstd, [2 1 3])); tmpdata = squeeze(real(nan_mean(tmpdata,3))); end;
                elseif strcmpi(opt.plotsubjects, 'off') tmpstd = squeeze(real(std(tmpdata,[],2)))/sqrt(size(tmpdata,2)); tmpstd = squeeze(permute(tmpstd, [2 1 3])); tmpdata = squeeze(real(nan_mean(tmpdata,2)));
                end;
            end;
            tmpdata = squeeze(permute(tmpdata, [2 1 3]));
            if strcmpi(opt.plottopo, 'on'), highlight = 'background'; else highlight = 'bottom'; end;
            if strcmpi(opt.plotgroups, 'together') &&  isempty(opt.condstats) && ...
                             ~isnan(opt.threshold) && ~isempty(opt.groupstats)
                plotopt = { plotopt{:} 'maskarray' };
                tmpdata = { tmpdata pgroupplot{c}' };
            elseif strcmpi(opt.plotconditions, 'together') &&  isempty(opt.groupstats) && ...
                                     ~isnan(opt.threshold) && ~isempty(opt.condstats)
                plotopt = { plotopt{:} 'maskarray' };
                tmpdata = { tmpdata pcondplot{g}' };
            end;
            plotopt = { plotopt{:} 'highlightmode', highlight };
            if strcmpi(opt.plotsubjects, 'on')
                plotopt = { plotopt{:} 'plotmean' 'on' 'plotindiv' 'on' };
            else
                plotopt = { plotopt{:} 'plotmean' 'off' };
            end;
            plotopt = { plotopt{:} 'ylim' opt.ylim 'xlabel' xlab 'ylabel' ylab };
            if ncplot ~= nc || ngplot ~= ng
                plotopt = { plotopt{:} 'legend' opt.legend };
            end;

%             % plot
%             % ----
%             if strcmpi(opt.figure, 'on'),
%                 tmpcol = col;
%             else
%                 % all groups and conditions in the same figureopt.chanlocs
%                 ncol = size(tmpdata,3);
%                 if ncol == 1, ncol = size(tmpdata,1); end;
%                 tmpcol = col(colcount:colcount+ncol-1);
%                 colcount = mod(colcount+ncol-1, length(col))+1;
%                 %if strcmpi(opt.plotgroups, 'together') && isempty(tmpdata{1})
%
%             end;
%
            if strcmpi(opt.plottopo, 'on') && length(opt.chanlocs) > 1
                metaplottopo(tmpdata, 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                    'plotargs', { plotopt{:} }, 'datapos', [2 3], 'title', opt.titles{c,g});
            elseif iscell(tmpdata)
                plotcurve( allx, tmpdata{1}, 'colors', tmpcol, 'maskarray', tmpdata{2}, plotopt{3:end}, 'title', opt.titles{c,g});
            else
                if isempty(findstr(opt.plotstderr, 'nocurve'))
                    plotcurve( allx, tmpdata, 'colors', tmpcol, plotopt{2:end}, 'title', opt.titles{c,g});
                end;
                if ~strcmpi(opt.plotstderr, 'off')
                    if ~dimreduced_sizediffers
                        if ~isempty(findstr(opt.plotstderr, 'diff')), begind = 3; else begind = 1; end;
                        set(gcf, 'renderer', 'OpenGL')
                        for tmpi = begind:size(tmpdata,1)
                            hold on; chandle = fillcurves( allx, tmpdata(tmpi,:)-tmpstd(tmpi,:), tmpdata(tmpi,:)+tmpstd(tmpi,:), tmpcol{tmpi}); hold on;
                            numfaces = size(get(chandle(1), 'Vertices'),1);
                            set(chandle(1), 'FaceVertexCData', repmat([1 1 1], [numfaces 1]), 'Cdatamapping', 'direct', 'facealpha', 0.3, 'edgecolor', 'none');
                        end;
                    else
                        disp('Some conditions have more subjects than others, cannot plot standard error');
                    end;
                end;
            end;
        end;

        if strcmpi(opt.plottopo, 'off'), % only non-topographic
            xlim([allx(1) allx(end)]); hold on;
            if isempty(opt.ylim)
                tmp = ylim;
                tmplim = [ min(tmplim(1), tmp(1)) max(tmplim(2), tmp(2)) ];
            else
                ylim(opt.ylim);
            end;
        end;

        % statistics accross groups
        % -------------------------
        if g == ngplot && ng > 1 && ~isempty(opt.groupstats)
            if ~strcmpi(opt.plotgroups, 'together') || ~isempty(opt.condstats) || isnan(opt.threshold)
                if strcmpi(opt.plotgroups, 'together'),         mysubplot(ncplot+addr, ngplot+addc, 2 + (c-1)*(ngplot+addc), opt.subplot); ci = g;
                elseif strcmpi(opt.plotconditions, 'together'), mysubplot(ncplot+addr, ngplot+addc, ngplot + 1, opt.subplot); ci = c;
                else                                            mysubplot(ncplot+addr, ngplot+addc, ngplot + 1 + (c-1)*(ngplot+addc), opt.subplot); ci = 1;
                end;
                if strcmpi(opt.plotconditions, 'together'), condnames = 'Conditions'; else condnames = opt.condnames{c}; end;
                if ~isnan(opt.threshold)
                     if strcmpi(opt.plottopo, 'on'),
                          metaplottopo({zeros(size(pgroupplot{c}')) pgroupplot{c}'}, 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                              'plotargs', { allx 'maskarray' statopt{:} }, 'datapos', [2 3], 'title', opt.titles{c, g+1});
                     else plotcurve(allx, zeros(size(allx)), 'maskarray', mean(pgroupplot{c},2), 'ylim', [0.1 1], 'title', opt.titles{c, g+1}, statopt{:});
                     end;
                else
                     if strcmpi(opt.plottopo, 'on'),
                          metaplottopo(pgroupplot{c}', 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                              'plotargs', { allx statopt{:} }, 'datapos', [2 3], 'title', opt.titles{c, g+1});
                     else plotcurve(allx, mean(pgroupplot{c},2), 'title', opt.titles{c, g+1}, statopt{:});
                     end;
                end;
            end;
        end;
    end;
end;

for g = 1:ng
    % statistics accross conditions
    % -----------------------------
    if ~isempty(opt.condstats) && nc > 1
        if ~strcmpi(opt.plotconditions, 'together') || ~isempty(opt.groupstats) || isnan(opt.threshold)
            if strcmpi(opt.plotgroups, 'together'),         mysubplot(ncplot+addr, ngplot+addc, 1 + c*(ngplot+addc), opt.subplot); ci = g;
            elseif strcmpi(opt.plotconditions, 'together'), mysubplot(ncplot+addr, ngplot+addc, g + ngplot+addc, opt.subplot); ci = c;
            else                                            mysubplot(ncplot+addr, ngplot+addc, g + c*(ngplot+addc), opt.subplot); ci = 1;
            end;
            if strcmpi(opt.plotgroups, 'together'), groupnames = 'Groups'; else groupnames = opt.groupnames{g}; end;
            if ~isnan(opt.threshold)
                 if strcmpi(opt.plottopo, 'on'),
                      metaplottopo({zeros(size(pcondplot{g}')) pcondplot{g}'}, 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                          'plotargs', { allx 'maskarray' statopt{:} }, 'datapos', [2 3], 'title', opt.titles{end, g});
                 else plotcurve(allx, zeros(size(allx)), 'maskarray', mean(pcondplot{g},2), 'ylim', [0.1 1], 'title', opt.titles{end, g}, statopt{:});
                 end;
            else
                 if strcmpi(opt.plottopo, 'on'),
                      metaplottopo(pcondplot{g}', 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                          'plotargs', { allx statopt{:} }, 'datapos', [2 3], 'title', opt.titles{end, g});
                 else plotcurve(allx, mean(pcondplot{g},2), 'title',  opt.titles{end, g}, statopt{:});
                 end;
            end;
        end;
    end;
end;

% statistics accross group and conditions
% ---------------------------------------
if ~isempty(opt.groupstats) && ~isempty(opt.condstats) && ng > 1 && nc > 1
    mysubplot(ncplot+addr, ngplot+addc, ngplot + 1 + ncplot*(ngplot+addr), opt.subplot);
    if ~isnan(opt.threshold)
         if strcmpi(opt.plottopo, 'on'),
              metaplottopo({zeros(size(pinterplot')) pinterplot'}, 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                  'plotargs', { allx 'maskarray' statopt{:} }, 'datapos', [2 3], 'title', opt.titles{end, end});
         else plotcurve(allx, zeros(size(allx)), 'maskarray', mean(pinterplot,2), 'ylim', [0.1 1], 'title', opt.titles{end, end}, statopt{:});
              xlabel(xlab); ylabel('-log10(p)');
        end;
    else
         if strcmpi(opt.plottopo, 'on'),
              metaplottopo(pinterplot', 'chanlocs', opt.chanlocs, 'plotfunc', 'plotcurve', ...
                  'plotargs', { allx statopt{:} }, 'datapos', [2 3], 'title', opt.titles{end, end});
         else plotcurve(allx, mean(pinterplot,2), 'title', opt.titles{end, end}, statopt{:});
         end;
    end;
end;

% axis limit
% ----------
for c = 1:ncplot
    for g = 1:ngplot
        if isempty(opt.ylim) && strcmpi(opt.plottopo, 'off')
            set(hdl(c,g), 'ylim', tmplim);
        end;
    end;
end;

if strcmpi(opt.plottopo, 'off') && length(hdl(:)) > 1
    axcopy;
    % remove axis labels (for most but not all)
    % ------------------
    if strcmpi(opt.subplot, 'transpose')
        for c = 1:size(hdl,2)
            for g = 1:size(hdl,1)
                axes(hdl(g,c));
                if c ~= 1 && size(hdl,2) ~=1, xlabel(''); legend off; end;
                if g ~= 1 && size(hdl,1) ~= 1, ylabel(''); legend off; end;
            end;
        end;
    else
        for c = 1:size(hdl,1)
            for g = 1:size(hdl,2)
                axes(hdl(c,g));
                if g ~= 1 && size(hdl,2) ~=1, ylabel(''); legend off; end;
                if c ~= size(hdl,1) && size(hdl,1) ~= 1, xlabel(''); legend off; end;
            end;
        end;
    end;
end;

% mysubplot (allow to transpose if necessary)
% -------------------------------------------
function hdl = mysubplot(nr,nc,ind,subplottype);

    r = ceil(ind/nc);
    c = ind -(r-1)*nc;
    if strcmpi(subplottype, 'transpose'),  hdl = subplot(nc,nr,(c-1)*nr+r);
    elseif strcmpi(subplottype, 'normal'), hdl = subplot(nr,nc,(r-1)*nc+c);
    elseif strcmpi(subplottype, 'noplot'), hdl = gca;
    else error('Unknown subplot type');
    end;

% rapid filtering for ERP
% -----------------------
function tmpdata2 = myfilt(tmpdata, srate, lowpass, highpass)

    tmpdata2 = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3)*size(tmpdata,4));
    tmpdata2 = eegfiltfft(tmpdata2',srate, lowpass, highpass)';
    tmpdata2 = reshape(tmpdata2, size(tmpdata,1), size(tmpdata,2), size(tmpdata,3), size(tmpdata,4));

author: PatternRecognition

repo: OpenBMI


text:
function varargout = real(varargin)
%REAL      Real part of a SPHEREFUN.
%
% See also SPHEREFUN/IMAG.

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

    [varargout{1:nargout}] = real@separableApprox(varargin{:});
end


text:
function obj = resetMijiClass(obj)
	% This clears Miji from Java's dynamic path and then re-initializes. Use if Miji is not loading normally.
	% Biafra Ahanonu
	% started: 2014.07.31 - branch from calciumImagingAnalysis 2020.05.07 [15:47:29]
	% inputs
		%
	% outputs
		%

	% changelog
		%
	% TODO
		%

	ciapkg.api.resetMiji;
	% success = 0;

	% for i = 1:2
	% 	try
	% 		% clear MIJ miji Miji mij;
	% 		javaDyna = javaclasspath('-dynamic');
	% 		matchIdx = ~cellfun(@isempty,regexpi(javaDyna,'Fiji'));
	% 		% cellfun(@(x) javarmpath(x),javaDyna(matchIdx));
	% 		javaDynaPathStr = join(javaDyna(matchIdx),''',''');
	% 		if ~isempty(javaDynaPathStr)
	% 			eval(sprintf('javarmpath(''%s'');',javaDynaPathStr{1}))
	% 		end
	% 		clear MIJ miji Miji mij;
	% 		% pause(1);
	% 		% java.lang.Runtime.getRuntime().gc;
	% 		% Miji;
	% 		% MIJ.exit;
	% 	catch err
	% 		disp(repmat('@',1,7))
	% 		disp(getReport(err,'extended','hyperlinks','on'));
	% 		disp(repmat('@',1,7))
	% 	end
	% end

	% success = 1;
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [tt]=tt_qlaplacez_dn(d)

% returns a rank-8 QTT decomposition of
% Delta_{1}^{-1} \otimes \Id_{2} \ otimes \ldots \otimes \Id_{D} + \ldots +
%  + \Id_{1} \ otimes \ldots \otimes \Id_{D-1} \otimes Delta_{D}^{-1},
% Delta_{k} being a discretization of Laplace operator on 2^{d(k)} points
% uniform grid,
% Dirichlet boundary conditions being imposed
%
% D=size(d,2) must be >= 1
%
% September 3, 2010
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
% Look for details in the Preprint No. 75, 2010 of
% Max-Planck Institute for Mathematics in the Sciences
% Vladimir A. Kazeev and Boris N. Khoromskij
% On explicit QTT representation of Laplace operator and its inverse
% http://www.mis.mpg.de/publications/preprints/2010/prepr2010-75.html

d=fliplr(d);
D=size(d,2);
tt=cell(sum(d),1);
I=eye(2);
J=zeros(2);
J(1,2)=1;
I2=zeros(2);
I2(2,2)=1;
E=ones(2);

if (D == 1)
	for key=1 : d
		tt{key}=eye(2);
	end
else

	key=0;
	for k=1 : D
		for kappa=1 : d(k)
			key=key+1;
			if (kappa == 1)
				if (k == 1)
					tt{key}=zeros(2,2,5);
					tt{key}(:,:,1)=I+I2+J+J';
					tt{key}(:,:,2)=2*E;
					tt{key}(:,:,3)=I2+J'+E;
					tt{key}(:,:,4)=I2+J+E;
					tt{key}(:,:,5)=I;
				elseif (k == D)
					tt{key}=zeros(2,2,2,4);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,1)=I+I2+J+J';
					tt{key}(:,:,2,2)=2*E;
					tt{key}(:,:,2,3)=I2+J'+E;
					tt{key}(:,:,2,4)=I2+J+E;
				else
					tt{key}=zeros(2,2,2,8);
					tt{key}(:,:,1,1)=I+I2+J+J';
					tt{key}(:,:,1,2)=2*E;
					tt{key}(:,:,1,3)=I2+J'+E;
					tt{key}(:,:,1,4)=I2+J+E;
					tt{key}(:,:,1,5)=I;
					tt{key}(:,:,2,5)=I+I2+J+J';
					tt{key}(:,:,2,6)=2*E;
					tt{key}(:,:,2,7)=I2+J'+E;
					tt{key}(:,:,2,8)=I2+J+E;
				end
			elseif (kappa == d(k))
				if (k == D)
					tt{key}=zeros(2,2,4);
					tt{key}(:,:,1)=I;
					tt{key}(:,:,2)=I2;
					tt{key}(:,:,3)=J;
					tt{key}(:,:,4)=J';
				elseif (k == 1)
					tt{key}=zeros(2,2,5,2);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,1)=I2;
					tt{key}(:,:,3,1)=J;
					tt{key}(:,:,4,1)=J';
					tt{key}(:,:,5,2)=I;
				else
					tt{key}=zeros(2,2,8,2);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,1)=I2;
					tt{key}(:,:,3,1)=J;
					tt{key}(:,:,4,1)=J';
					tt{key}(:,:,5,2)=I;
					tt{key}(:,:,6,2)=I2;
					tt{key}(:,:,7,2)=J;
					tt{key}(:,:,8,2)=J';
				end
			else
				if (k == D)
					tt{key}=zeros(2,2,4,4);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,2)=2*E;
					tt{key}(:,:,3,3)=E;
					tt{key}(:,:,4,4)=E;
					tt{key}(:,:,2,1)=I2;
					tt{key}(:,:,3,1)=J;
					tt{key}(:,:,4,1)=J';
					tt{key}(:,:,2,3)=I2+J';
					tt{key}(:,:,2,4)=I2+J;
				elseif (k == 1)
					tt{key}=zeros(2,2,5,5);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,2)=2*E;
					tt{key}(:,:,3,3)=E;
					tt{key}(:,:,4,4)=E;
					tt{key}(:,:,2,1)=I2;
					tt{key}(:,:,3,1)=J;
					tt{key}(:,:,4,1)=J';
					tt{key}(:,:,2,3)=I2+J';
					tt{key}(:,:,2,4)=I2+J;
					tt{key}(:,:,5,5)=I;
				else
					tt{key}=zeros(2,2,8,8);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,2)=2*E;
					tt{key}(:,:,3,3)=E;
					tt{key}(:,:,4,4)=E;
					tt{key}(:,:,2,1)=I2;
					tt{key}(:,:,3,1)=J;
					tt{key}(:,:,4,1)=J';
					tt{key}(:,:,2,3)=I2+J';
					tt{key}(:,:,2,4)=I2+J;
					tt{key}(:,:,5,5)=I;
					tt{key}(:,:,6,6)=2*E;
					tt{key}(:,:,7,7)=E;
					tt{key}(:,:,8,8)=E;
					tt{key}(:,:,6,5)=I2;
					tt{key}(:,:,7,5)=J;
					tt{key}(:,:,8,5)=J';
					tt{key}(:,:,6,7)=I2+J';
					tt{key}(:,:,6,8)=I2+J;
				end
			end
		end
	end
end

tt=tt_matrix(tt); % @Bydlocode
return
end

author: oseledets

repo: TT-Toolbox


text:
function results = poly_segmentation_3d(record)
if nargin < 1
  record = false;
end
import iris.drawing.*;

lb = [0;0;0];
ub = [10;10;10];
dim = 3;

n_obs = 20;
% obstacles = iris.test.random_obstacles(dim, n_obs, lb, ub,3);
obstacles = iris.test.random_cubic_obstacles(dim, n_obs, lb, ub);

A_bounds = [-1,0,0;
            0,-1,0;
            0,0,-1;
            1,0,0;
            0,1,0;
            0,0,1];
b_bounds = [-lb;ub];

start = 0.5 * (ub + lb);

% profile on
[A,b,C,d,results] = iris.inflate_region(obstacles, A_bounds, b_bounds, start);
% profile viewer
if n_obs < 50
  animate_results(results, record);
end

end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

X = [ones(m, 1) X];

for i=1:m
    a1 = X(i,:)';
    a2 = [1; sigmoid(Theta1 * a1)];
    a3 = sigmoid(Theta2 * a2);
    [m, idx] = max(a3');
    p(i) = idx;
end;

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

end

author: zsiciarz

repo: ml-coursera


text:
function symbols = computeSymbolsFromParams(par, sampleRate, iq)

s0 = par(1);
s = par(2);
skew = par(3);
rot = par(4);

%% Adjust the scale

iq = iq * s;

%% Isolate the Symbols
% Because of differences between the rate at which symbols are broadcast,
% the rate at which samples are taken, and the blurring effect of
% transmission, the symbols are uniformly spaced throughout the acquired IQ
% data.  The remainder of the samples are transitions between symbols.
% From the capture time and the symbol rate, we know how many symbols.  By
% combining that with the sample rate, we can calculate the symbol spacing
% in our IQ cloud.
% By using the value manipulation tools in the editor and Cell mode, we can
% immediately visualize the results of different symbol choices.
% I know the symbol spacing - knows how sig broadcasted and how signature
% is sampling

manualSymbols= iq(round(s0:sampleRate:length(iq)));

%% What's Wrong Now!
% At this point, we have located our symbols in the IQ data.  However, the
% symbols still don't cluster into eight nice neat groupings.  The problem
% is that we need to account for skew in the signal.  If everything is not
% synchronized correctly, the symbols appear to "walk" around between their
% intended locations.  We can correct for that with a linear phase
% correction.

skewCorrection = exp(2*pi*i*linspace(0,skew,length(manualSymbols)));
manualSymbolsDeskewed = manualSymbols.*skewCorrection; %*exp(2*pi*i*0.04) adds global rotation correction

%% Global Rotation
% Now that we have located our symbols and corrected for skew between the
% clocks, if our clusters aren't aligned correctly globally (all shited
% clockwise or counterclockwise) we can apply a global rotation.

symbols = manualSymbolsDeskewed.*exp(2*pi*i*rot);


text:
% std_readdata() - LEGACY FUNCTION, SHOULD NOT BE USED ANYMORE. INSTEAD
%                  USE std_readerp, std_readspec, ...
%                  load one or more requested measures
%                  ['erp'|'spec'|'ersp'|'itc'|'dipole'|'map']
%                  for all components of a specified cluster.
%                  Called by cluster plotting functions
%                  std_envtopo(), std_erpplot(), std_erspplot(), ...
% Usage:
%         >> [STUDY, clustinfo, finalinds] = std_readdata(STUDY, ALLEEG);
%                                              % return all measures
%         >> [STUDY, clustinfo, finalinds] = std_readdata(STUDY, ALLEEG, ...
%                                              cluster, infotype,varargin);
% Inputs:
%       STUDY - studyset structure containing some or all files in ALLEEG
%      ALLEEG - vector of loaded EEG datasets
%
% Optional inputs:
%  'infotype'  - ['erp'|'spec'|'ersp'|'itc'|'dipole'|'map'|'data'|'event']
%                type of stored data or cluster information to read. May
%                also be a cell array of these types, for example: { 'erp'
%                'map' 'dipole' }. {default: 'erp'}
%  'channels'  - [cell] list of channels to import {default: all}
%  'clusters'  - [integer] list of clusters to import {[]|default: all but
%                the parent cluster (1) and any 'NotClust' clusters}
%  'freqrange' - [min max] frequency range {default: whole measure range}
%  'timerange' - [min max] time range {default: whole measure epoch}
%  'statmode'  - ['individual'|'trials'] statistical mode for ERSP (also
%                specify what type of data to import -- mean (of individual
%                subjects), or trials. This functionality is still unstable
%                for 'trials' { default: 'individual'}
%  'rmsubjmean' - ['on'|'off'] remove mean subject spectrum from every component
%                 or channel spectrum, making them easier to compare
%                 { default: 'off' }
%  'subbaseline' - ['on'|'off'] remove all condition and spectrum baselines
%                  for ERSP data { default: 'on' }
%
% Optional inputs for events:
%  'type','timewin','fieldname' - optional parameters for the function
%                  eeg_geteventepoch may also be used to select specific
%                  events.
%
% Optional inputs for raw data:
%  'rmclust'  - [integer] list of artifact cluster indices to subtract
%               when reading raw data.
%
% Output:
%    STUDY     - (possibly) updated STUDY structure
%    clustinfo - structure of specified cluster information.
%                This is the same as STUDY.cluster(cluster_number)
%    Fields:
%     clustinfo.erpdata    % (ncomps, ntimes) array of component ERPs
%     clustinfo.erptimes   % vector of ERP epoch latencies (ms)
%
%     clustinfo.specdata   % (ncomps, nfreqs) array of component spectra
%     clustinfo.specfreqs  % vector of spectral frequencies (Hz)
%
%     clustinfo.erspdata   % (ncomps,ntimes,nfreqs) array of component ERSPs
%      clustinfo.ersptimes % vector of ERSP latencies (ms)
%      clustinfo.erspfreqs % vector of ERSP frequencie (Hz)
%
%     clustinfo.itcdata    % (ncomps,ntimes,nfreqs) array of component ITCs
%      clustinfo.itctimes  % vector of ITC latencies (ms)
%      clustinfo.itc_freqs % vector of ITC frequencies (Hz)
%
%     clustinfo.topo       % (ncomps,65,65) array of component scalp map grids
%       clustinfo.topox    % abscissa values for columns of the scalp maps
%       clustinfo.topoy    % ordinate values for rows of the scalp maps
%
%     clustinfo.data           % raw data arrays
%       clustinfo.datatimes    % time point for raw data
%       clustinfo.datasortvals % sorting values (one per trial)
%       clustinfo.datacontinds % dataset indices (contacted for all trials)
%
%     clustinfo.dipole     % array of component dipole information structs
%                          % with the same format as EEG.dipfit.model
%
%   finalinds - either the cluster(s) or channel(s) indices selected.
%
% Example:
%         % To obtain the ERPs for all Cluster-3 components from a STUDY
%         %
%         [STUDY clustinfo] = std_readdata(STUDY, ALLEEG, 'clusters',3, 'infotype','erp');
%         figure; plot(clustinfo.erptimes, mean(clustinfo.erpdata,2));
%
% Author: Arnaud Delorme, CERCO, 2006-

% Copyright (C) Arnaud Delorme, arno@salk.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

function [STUDY, clustinfo, finalinds] = std_readdata(STUDY, ALLEEG, varargin);

if nargin < 2
    help std_readdata;
    return;
end
if nargin < 3
    channel = {};
end

[opt moreopts] = finputcheck( varargin, { ...
    'type'       { 'string' 'cell' } { [] [] } '';
    'timewin'    'real'    []        [-Inf Inf];
    'fieldname'  'string'  []        'latency';
    'condition'  'cell'    []       {};
    'channels'   'cell'    []       {};
    'clusters'   'integer' []       [];
    'rmclust'    'integer' []       [];
    'freqrange'  'real'    []       [];
    'timerange'  'real'    []       [];
    'statmode'   'string'  { 'subjects' 'individual' 'common' 'trials' }       'individual';
    'rmicacomps'   'string'  { 'on' 'off' }       'off';
    'subbaseline'  'string'  { 'on' 'off' }       'off';
    'rmsubjmean'   'string'  { 'on' 'off' }       'off';
    'singletrials' 'string'  { 'on' 'off' }       'on';
    'infotype'   'string'  { 'erp' 'spec' 'ersp' 'itc' 'map' 'topo' 'dipole' 'scalp' 'data' 'event' '' } '' }, ...
    'std_readdata', 'ignore');
if isstr(opt), error(opt); end;
if isempty(opt.infotype), disp('No measure selected, returning ERPs'); opt.infotype = 'erp'; end;

if strcmpi(opt.infotype, 'erp'),
    STUDY = pop_erpparams(STUDY, 'default');
    if isempty(opt.timerange), opt.timerange = STUDY.etc.erpparams.timerange; end;
elseif strcmpi(opt.infotype, 'spec'),
    STUDY = pop_specparams(STUDY, 'default');
    if isempty(opt.freqrange), opt.freqrange = STUDY.etc.specparams.freqrange; end;
elseif strcmpi(opt.infotype, 'ersp') | strcmpi(opt.infotype, 'itc')
    STUDY = pop_erspparams(STUDY, 'default');
    if isempty(opt.freqrange), opt.freqrange   = STUDY.etc.erspparams.freqrange; end;
    if isempty(opt.timerange), opt.timerange   = STUDY.etc.erspparams.timerange; end;
    if strcmpi(opt.statmode, 'individual') |  strcmpi(opt.statmode, 'subjects'), opt.statmode    = STUDY.etc.erspparams.statmode;end;
    if strcmpi(opt.subbaseline, 'on'),      opt.subbaseline = STUDY.etc.erspparams.subbaseline; end;
end;

nc = max(length(STUDY.condition),1);
ng = max(length(STUDY.group),1);
zero = single(0);
tmpver = version;
if tmpver(1) == '6' | tmpver(1) == '5'
    zero = double(zero);
end;

% find channel indices
% --------------------
if ~isempty(opt.channels)
    finalinds = std_chaninds(STUDY, opt.channels);
else
    finalinds = opt.clusters;
end;

% read topography with another function
% -------------------------------------
if strcmpi(opt.infotype, 'map') | strcmpi(opt.infotype, 'scalp') | strcmpi(opt.infotype, 'topo')
    [STUDY tmpclust] = std_readtopoclust(STUDY, ALLEEG, opt.clusters);
    clustinfo = [];
    for index = 1:length(tmpclust)
        if index == 1, clustinfo        = tmpclust{index};
        else           clustinfo(index) = tmpclust{index};
        end;
    end;
    return;
end;

% read indices for removing component clusters from data
% ------------------------------------------------------
% NOT USED ANY MORE
if ~isempty(opt.rmclust)
    for ind = 1:length(opt.rmclust)
        [tmp setindsrm{ind} allindsrm{ind}] = std_setinds2cell(STUDY, opt.rmclust(ind));
    end;
end;

for ind = 1:length(finalinds)

    % find indices
    % ------------
    if ~isempty(opt.channels)
        tmpstruct = STUDY.changrp(finalinds(ind));
        alldatasets = 1:length(STUDY.datasetinfo);
        %allchanorcomp = -tmpstruct.chaninds;
        allinds       = tmpstruct.allinds;
        for i=1:length(allinds(:)), allinds{i} = -allinds{i}; end; % invert sign for reading
        setinds       = tmpstruct.setinds;
    else
        % Use the 3 lines below and comment the last one when ready
        % to switch and remove all .comps and .inds
        %tmpstruct = STUDY.cluster(finalinds(ind));
        %allinds = tmpstruct.allinds;
        %setinds = tmpstruct.setinds;
        [ tmpstruct setinds allinds ] = std_setcomps2cell(STUDY, finalinds(ind));
        %[ STUDY.cluster(finalinds(ind)) ] = std_setcomps2cell(STUDY, finalinds(ind));
        %STUDY.cluster(finalinds(ind))     = std_cell2setcomps( STUDY, ALLEEG, finalinds(ind));
        %[ STUDY.cluster(finalinds(ind)) setinds allinds ] = std_setcomps2cell(STUDY, finalinds(ind));
    end;

    dataread = 0;
    switch opt.infotype
        case 'event',
            % check if data is already here
            % -----------------------------
            if isfield(tmpstruct, 'datasortvals')
                if ~isempty(tmpstruct.datasortvals)
                    dataread = 1;
                end;
            end;

            if ~dataread
                % reserve arrays
                % --------------
                datasortvals   = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                datacontinds    = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                for c = 1:nc
                    for g = 1:ng
                        datasortvals{c, g} = repmat(zero, [1, sum([ ALLEEG(setinds{c,g}).trials ] )]);
                        datacontinds{c, g} = repmat(zero, [1, sum([ ALLEEG(setinds{c,g}).trials ] )]);
                    end;
                end;

                % read the data and select channels
                % ---------------------------------
                if ~isempty(opt.type)
                    fprintf('Reading events:');
                    for c = 1:nc
                        for g = 1:ng
                            counttrial = 1;
                            for indtmp = 1:length(allinds{c,g})
                                tmpdata  = eeg_getepochevent(ALLEEG(setinds{c,g}(indtmp)), ...
                                    'type', opt.type, 'timewin', opt.timewin, 'fieldname', opt.fieldname);
                                datasortvals{c, g}(:,counttrial:counttrial+ALLEEG(setinds{c,g}(indtmp)).trials-1) = squeeze(tmpdata);
                                datacontinds{c, g}(:,counttrial:counttrial+ALLEEG(setinds{c,g}(indtmp)).trials-1) = setinds{c,g}(indtmp);
                                counttrial = counttrial+ALLEEG(setinds{c,g}(indtmp)).trials;
                                fprintf('.');
                            end;
                        end;
                    end;
                end;
                fprintf('\n');
                tmpstruct.datasortvals = datasortvals;
                tmpstruct.datacontinds = datacontinds;
            end;

        case 'data',
            % check if data is already here
            % -----------------------------
            if isfield(tmpstruct, 'data')
                if ~isempty(tmpstruct.data)
                    dataread = 1;
                end;
            end;

            if ~dataread
                % reserve arrays
                % --------------
                alldata   = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                for c = 1:nc
                    for g = 1:ng
                        alldata{c, g} = repmat(zero, [ALLEEG(1).pnts, sum([ ALLEEG(setinds{c,g}).trials ] )]);
                    end;
                end;

                % read the data and select channels
                % ---------------------------------
                fprintf('Reading data/activity:');
                for c = 1:nc
                    for g = 1:ng
                        counttrial = 1;
                        for indtmp = 1:length(allinds{c,g})
                            settmpind = STUDY.datasetinfo(setinds{c,g}(indtmp)).index; % current dataset

                            if isempty(opt.channels)
                                tmpdata = eeg_getdatact(ALLEEG(settmpind), 'component', allinds{c,g}(indtmp), 'verbose', 'off');
                            elseif isempty(opt.rmclust) & strcmpi(opt.rmicacomps, 'off')
                                tmpdata = eeg_getdatact(ALLEEG(settmpind),  'channel', -allinds{c,g}(indtmp), 'verbose', 'off');
                            elseif strcmpi(opt.rmicacomps, 'on')
                                rmcomps = find(ALLEEG(settmpind).reject.gcompreject);
                                tmpdata = eeg_getdatact(ALLEEG(settmpind),  'channel', -allinds{c,g}(indtmp), 'rmcomps', rmcomps, 'verbose', 'off');
                            else
                                rmcomps = getclustcomps(STUDY, opt.rmclust, settmpind);
                                tmpdata = eeg_getdatact(ALLEEG(settmpind),  'channel', -allinds{c,g}(indtmp), 'rmcomps', rmcomps, 'verbose', 'off');
                            end;

                            alldata{c, g}(:,counttrial:counttrial+ALLEEG(setinds{c,g}(indtmp)).trials-1) = squeeze(tmpdata);
                            counttrial = counttrial+ALLEEG(setinds{c,g}(indtmp)).trials;
                            fprintf('.');
                        end;
                    end;
                end;
                fprintf('\n');
                tmpstruct.datatimes = ALLEEG(1).times;
                tmpstruct.data      = alldata;
            end;

        case 'erp',
            disp('std_readdata is a legacy function, it should not be used to read ERP data, use std_readerp instead');

            % check if data is already here
            % -----------------------------
            if isfield(tmpstruct, 'erpdata')
                if isequal( STUDY.etc.erpparams.timerange, opt.timerange) & ~isempty(tmpstruct.erpdata)
                    dataread = 1;
                end;
            end;

            if ~dataread
                % reserve arrays
                % --------------
                allerp   = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                tmpind = 1; while(isempty(setinds{tmpind})), tmpind = tmpind+1; end;
                [ tmp alltimes ] = std_readerp( ALLEEG, setinds{tmpind}(1), allinds{tmpind}(1), opt.timerange);
                for c = 1:nc
                    for g = 1:ng
                        allerp{c, g} = repmat(zero, [length(alltimes), length(allinds{c,g})]);
                    end;
                end;

                % read the data and select channels
                % ---------------------------------
                fprintf('Reading ERP data:');
                for c = 1:nc
                    for g = 1:ng
                        for indtmp = 1:length(allinds{c,g})
                            [ tmperp alltimes ] = std_readerp( ALLEEG, setinds{c,g}(indtmp), allinds{c,g}(indtmp), opt.timerange);
                            allerp{c, g}(:,indtmp) = tmperp(:);
                            fprintf('.');
                        end;
                    end;
                end;
                fprintf('\n');
                tmpstruct.erptimes = alltimes;
                tmpstruct.erpdata  = allerp;
            end;

        case 'spec',
            disp('std_readdata is a legacy function, it should not be used to read SPECTRUM data, use std_readspec instead');
            % check if data is already here
            % -----------------------------
            if isfield(tmpstruct, 'specdata')
                if isequal( STUDY.etc.specparams.freqrange, opt.freqrange) & ~isempty(tmpstruct.specdata)
                    dataread = 1;
                end;
            end;

            if ~dataread
                % reserve arrays
                % --------------
                allspec  = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                filetype = 'spec';
                tmpind = 1; while(isempty(setinds{tmpind})), tmpind = tmpind+1; end;
                try,
                    [ tmp allfreqs ] = std_readspec( ALLEEG, setinds{tmpind}(1), allinds{tmpind}(1), opt.freqrange);
                catch
                    filetype = 'ersp';
                    disp('Cannot find spectral file, trying ERSP baseline file instead');
                    [ tmpersp allfreqs alltimes tmpparams tmpspec] = std_readersp( ALLEEG, setinds{tmpind}(1), allinds{tmpind}(1), [], opt.freqrange);
                end;
                for c = 1:nc
                    for g = 1:ng
                        allspec{c, g} = repmat(zero, [length(allfreqs), length(allinds{c,g}) ]);
                    end;
                end;

                % read the data and select channels
                % ---------------------------------
                if strcmpi(filetype, 'spec')
                    fprintf('Reading Spectrum data...');
                    for c = 1:nc
                        for g = 1:ng
                            allspec{c, g} = std_readspec( ALLEEG, setinds{c,g}(:), allinds{c,g}(:), opt.freqrange)';
                        end;
                    end;
                else % std_readersp cannot be converted to read multiple datasets since it subtracts data between conditions
                    fprintf('Reading Spectrum data:');
                    for c = 1:nc
                        for g = 1:ng
                            for indtmp = 1:length(allinds{c,g})
                                [ tmpersp allfreqs alltimes tmpparams tmpspec] = std_readersp( ALLEEG, setinds{c,g}(indtmp), allinds{c,g}(indtmp), [], opt.freqrange);
                                allspec{c, g}(:,indtmp) = 10*log(tmpspec(:));
                                fprintf('.');
                            end;
                        end;
                    end;
                end;
                fprintf('\n');

                % remove mean of each subject across groups and conditions
                if strcmpi(opt.rmsubjmean, 'on') & ~isempty(opt.channels)
                    disp('Removing mean spectrum accross subjects');
                    for indtmp = 1:length(allinds{c,g}) % scan subjects
                       meanspec =zeros(size( allspec{1, 1}(:,indtmp) ));
                       for c = 1:nc
                            for g = 1:ng
                                meanspec = meanspec + allspec{c, g}(:,indtmp)/(nc*ng);
                            end;
                       end;
                       for c = 1:nc
                            for g = 1:ng
                                allspec{c, g}(:,indtmp) = allspec{c, g}(:,indtmp) - meanspec; % subtractive model
                                % allspec{c, g}(:,indtmp) = allspec{c, g}(:,indtmp)./meanspec; % divisive model
                            end;
                       end;
                    end;
                end;

                tmpstruct.specfreqs = allfreqs;
                tmpstruct.specdata  = allspec;
            end;

        case { 'ersp' 'itc' 'pac' },
            disp('std_readdata is a legacy function, it should not be used to read ERSP/ITC data, use std_readersp instead');
            % check if data is already here
            % -----------------------------
            if strcmpi(opt.infotype, 'ersp')
                if isfield(tmpstruct, 'erspdata')
                    if isequal( STUDY.etc.erspparams.timerange, opt.timerange) & ...
                            isequal( STUDY.etc.erspparams.freqrange, opt.freqrange) & ~isempty(tmpstruct.erspdata)
                        dataread = 1;
                    end;
                end;
            elseif strcmpi(opt.infotype, 'itc')
                if isfield(tmpstruct, 'itcdata')
                    if isequal( STUDY.etc.erspparams.timerange, opt.timerange) & ...
                            isequal( STUDY.etc.erspparams.freqrange, opt.freqrange) & ~isempty(tmpstruct.itcdata)
                        dataread = 1;
                    end;
                end;
            else
                if isfield(tmpstruct, 'pacdata')
                    dataread = 1;
                end;
            end;

            if ~dataread
                % find total nb of trials
                % -----------------------
                if strcmpi(opt.statmode, 'trials')
                    tottrials = cell( length(STUDY.condition), length(STUDY.group) );
                    for index = 1:length( STUDY.datasetinfo )
                        condind = strmatch( STUDY.datasetinfo(index).condition, STUDY.condition );
                        grpind  = strmatch( STUDY.datasetinfo(index).group    , STUDY.group     );
                        if isempty(tottrials{condind, grpind}), tottrials{condind, grpind} = ALLEEG(index).trials;
                        else       tottrials{condind, grpind} = tottrials{condind, grpind} + ALLEEG(index).trials;
                        end;
                    end;
                end;

                % reserve arrays
                % --------------
                ersp     = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                erspbase = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                erspinds = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
                tmpind = 1; while(isempty(setinds{tmpind})), tmpind = tmpind+1; end;
                [ tmp allfreqs alltimes ] = std_readersp( ALLEEG, setinds{tmpind}(1), Inf*allinds{tmpind}(1), opt.timerange, opt.freqrange);
                for c = 1:nc
                    for g = 1:ng
                        ersp{c, g}     = repmat(zero, [length(alltimes), length(allfreqs), length(allinds{c,g}) ]);
                        erspbase{c, g} = repmat(zero, [               1, length(allfreqs), length(allinds{c,g}) ]);
                        if strcmpi(opt.statmode, 'trials')
                            ersp{ c, g} = repmat(zero, [length(alltimes), length(allfreqs), tottrials{c, g} ]);
                            count{c, g} = 1;
                        end;
                    end;
                end;

                % read the data and select channels
                % ---------------------------------
                fprintf('Reading all %s data:', upper(opt.infotype));
                for c = 1:nc
                    for g = 1:ng
                        for indtmp = 1:length(allinds{c,g})
                            try     % this 'try' is a poor solution the problem of attempting
                                % to read specific channel/component data that doesn't exist
                                % called below by: allinds{c,g}(indtmp)
                                if strcmpi(opt.statmode, 'trials')
                                    [ tmpersp allfreqs alltimes tmpparams] = std_readtimef( ALLEEG, setinds{c,g}(indtmp), allinds{c,g}(indtmp), ...
                                        opt.timerange, opt.freqrange);
                                    indices = [count{c, g}:count{c, g}+size(tmpersp,3)-1];
                                    ersp{    c, g}(:,:,indices) = permute(tmpersp, [2 1 3]);
                                    erspinds{c, g}(1:2,indtmp) = [ count{c, g} count{c, g}+size(tmpersp,3)-1 ];
                                    count{c, g} = count{c, g}+size(tmpersp,3);
                                    if size(tmpersp,3) ~= ALLEEG(STUDY.datasetinfo(index).index).trials
                                        error( sprintf('Wrong number of trials in datafile for dataset %d\n', STUDY.datasetinfo(index).index));
                                    end;
                                elseif strcmpi(opt.infotype, 'itc')
                                    [ tmpersp allfreqs alltimes tmpparams] = std_readitc( ALLEEG, setinds{c,g}(indtmp), allinds{c,g}(indtmp), ...
                                        opt.timerange, opt.freqrange);
                                    ersp{c, g}(:,:,indtmp)     = abs(permute(tmpersp    , [2 1]));
                                elseif strcmpi(opt.infotype, 'pac')
                                    [ tmpersp allfreqs alltimes tmpparams] = std_readpac( ALLEEG, setinds{c,g}(indtmp), allinds{c,g}(indtmp), ...
                                        opt.timerange, opt.freqrange);
                                    ersp{c, g}(:,:,indtmp)     = abs(permute(tmpersp    , [2 1]));
                                else
                                    [ tmpersp allfreqs alltimes tmpparams tmperspbase] = std_readersp( ALLEEG, setinds{c,g}(indtmp), allinds{c,g}(indtmp), ...
                                        opt.timerange, opt.freqrange);
                                    ersp{c, g}(    :,:,indtmp) = permute(tmpersp    , [2 1]);
                                    erspbase{c, g}(:,:,indtmp) = 10*log(permute(tmperspbase, [2 1]));
                                end
                            catch
                            end
                            fprintf('.');
                        end;
                    end;
                end;
                fprintf('\n');

                % compute ERSP or ITC if trial mode
                % (since only the timef have been loaded)
                % ---------------------------------------
                if strcmpi(opt.statmode, 'trials')
                    for c = 1:nc
                        for g = 1:ng
                            if strcmpi(opt.infotype, 'itc')
                                ersp{c,g} = ersp{c,g}./abs(ersp{c,g});
                            else
                                ersp{c,g} = 20*log10(abs(ersp{c,g}));
                                % remove baseline (each trial baseline is removed => could
                                % also be the average across all data trials)
                                [tmp indl] = min( abs(alltimes-0) );
                                erspbase{c,g} = mean(ersp{c,g}(1:indl,:,:,:));
                                ersp{c,g} = ersp{c,g} - repmat(erspbase{c,g}, [size(ersp{c,g},1) 1 1 1]);
                            end;
                        end;
                    end;
                end;

                % compute average baseline across groups and conditions
                % -----------------------------------------------------
                if strcmpi(opt.subbaseline, 'on') & strcmpi(opt.infotype, 'ersp')
                    disp('Recomputing baseline...');
                    for g = 1:ng        % ng = number of groups
                        for c = 1:nc    % nc = number of components
                            if strcmpi(opt.statmode, 'trials')
                                if c == 1, meanpowbase = abs(mean(erspbase{c,g}/nc,3));
                                else       meanpowbase = meanpowbase + abs(mean(erspbase{c,g}/nc,3));
                                end;
                            else
                                if c == 1, meanpowbase = abs(erspbase{c,g}/nc);
                                else       meanpowbase = meanpowbase + abs(erspbase{c,g}/nc);
                                end;
                            end;
                        end;

                        % subtract average baseline
                        % -------------------------
                        for c = 1:nc
                            if strcmpi(opt.statmode, 'trials'), tmpmeanpowbase = repmat(meanpowbase, [length(alltimes) 1 tottrials{c,g}]);
                            else                                tmpmeanpowbase = repmat(meanpowbase, [length(alltimes) 1 1]);
                            end;
                            ersp{c,g} = ersp{c,g} - repmat(abs(erspbase{c,g}), [length(alltimes) 1 1 1]) + tmpmeanpowbase;
                        end;
                        clear meanpowbase;
                    end;
                end;

                if strcmpi(opt.statmode, 'common')
                    % collapse the two last dimensions before computing significance
                    % i.e. 18 subject with 4 channels -> same as 4*18 subjects
                    % --------------------------------------------------------------
                    disp('Using all channels for statistics...');
                    for c = 1:nc
                        for g = 1:ng
                            ersp{c,g} = reshape( ersp{c,g}, size(ersp{c,g},1), size(ersp{c,g},2), size(ersp{c,g},3)*size(ersp{c,g},4));
                        end;
                    end;
                end;

                % copy data to structure
                % ----------------------
                if strcmpi(opt.infotype, 'ersp')
                    tmpstruct.erspfreqs = allfreqs;
                    tmpstruct.ersptimes = alltimes;
                    tmpstruct.erspdata  = ersp;
                    tmpstruct.erspbase  = erspbase;
                    if strcmpi(opt.statmode, 'trials')
                        tmpstruct.erspsubjinds  = erspinds;
                    end;
                elseif strcmpi(opt.infotype, 'itc')
                    tmpstruct.itcfreqs = allfreqs;
                    tmpstruct.itctimes = alltimes;
                    tmpstruct.itcdata  = ersp;
                    if strcmpi(opt.statmode, 'trials')
                        tmpstruct.itcsubjinds  = erspinds;
                    end;
                else
                    tmpstruct.pacfreqs = allfreqs;
                    tmpstruct.pactimes = alltimes;
                    tmpstruct.pacdata  = ersp;
                end;
            end;
        case 'dipole',
            fprintf('Reading dipole data...\n');
%           old format EEGLAB 8.3
%             alldips = {};
%             for c = 1:nc
%                 for g = 1:ng
%                     for indtmp = 1:size(tmpstruct.sets,2)
%                         alldips{c, g}(indtmp).posxyz = ALLEEG(tmpstruct.sets(1,indtmp)).dipfit.model(tmpstruct.comps(1,indtmp)).posxyz;
%                         alldips{c, g}(indtmp).momxyz = ALLEEG(tmpstruct.sets(1,indtmp)).dipfit.model(tmpstruct.comps(1,indtmp)).momxyz;
%                         alldips{c, g}(indtmp).rv     = ALLEEG(tmpstruct.sets(1,indtmp)).dipfit.model(tmpstruct.comps(1,indtmp)).rv;
%                     end;
%                 end;
%             end;
%            tmpstruct.dipoles = alldips;
            alldips = [];
            for indtmp = 1:size(tmpstruct.sets,2)
                alldips(indtmp).posxyz = ALLEEG(tmpstruct.sets(1,indtmp)).dipfit.model(tmpstruct.comps(1,indtmp)).posxyz;
                alldips(indtmp).momxyz = ALLEEG(tmpstruct.sets(1,indtmp)).dipfit.model(tmpstruct.comps(1,indtmp)).momxyz;
                alldips(indtmp).rv     = ALLEEG(tmpstruct.sets(1,indtmp)).dipfit.model(tmpstruct.comps(1,indtmp)).rv;
            end;
            tmpstruct.alldipoles = alldips;

        case { 'map' 'scalp' 'topo' }
            % this is currenlty being done by the function std_readtopoclust
            % at the beginning of this function
        otherwise, error('Unrecognized ''infotype'' entry');
    end; % end switch

    % copy results to structure
    % -------------------------
    fieldnames = { 'erpdata' 'erptimes' 'specdata' 'specfreqs' 'erspdata' 'erspbase' 'erspfreqs' 'ersptimes' ...
        'itcfreqs' 'itctimes' 'itcdata' 'erspsubjinds' 'itcsubjinds' 'allinds' 'setinds' 'dipoles' 'alldipoles' ...
        'data' 'datatimes' 'datasortvals' 'datacontinds' };
    for f = 1:length(fieldnames)
        if isfield(tmpstruct, fieldnames{f}),
            tmpdata = getfield(tmpstruct, fieldnames{f});
            if ~isempty(opt.channels)
                STUDY.changrp = setfield(STUDY.changrp, { finalinds(ind) }, fieldnames{f}, tmpdata);
            else STUDY.cluster = setfield(STUDY.cluster, { finalinds(ind) }, fieldnames{f}, tmpdata);
            end;
        end;
    end;
end;

% return structure
% ----------------
if ~isempty(opt.channels)
    clustinfo = STUDY.changrp(finalinds);
else clustinfo = STUDY.cluster(finalinds);
end;

% find components in cluster for specific dataset
% -----------------------------------------------
function rmcomps = getclustcomps(STUDY, rmclust, settmpind);

rmcomps   = [ ];
for rmi = 1:length(rmclust)
    findind   = find(settmpind == STUDY.cluster(rmclust(rmi)).setinds{c,g});
    rmcomps   = [ rmcomps STUDY.cluster(rmclust(rmi)).allinds{c,g}(findind) ];
end;