text:
function view=analysisFlatMenu(view)
%
% view=analysisFlatMenu(view)
%
% Set up the callbacks of the ANALYSIS menu for the FLAT view.
%
% wap, 5/18/99
% arp, 052199
%

analysismenu = uimenu('Label','Analysis','separator','on');

% Cortical distance measurement (prompt for user input)
callBackstr=['measureCorticalDistance(' view.name ');'];
uimenu(analysismenu,'Label','Measure cortical distance (select points)','Separator','off',...
    'CallBack',callBackstr);

% Cortical distance measurement (use current ROI)
callBackstr=['measureCorticalDistance(' view.name ', getCurROIcoords(' view.name ') );'];
uimenu(analysismenu,'Label','Measure cortical distance (current ROI)','Separator','off',...
    'CallBack',callBackstr);

% Cortical distance measurement- sub-sampling method (use current line ROI)
callBackstr=['measureCorticalDistanceLineROI(' view.name ', getCurROIcoords(' view.name ') );'];
uimenu(analysismenu,'Label','Measure cortical distance (subsamples ROI)','Separator','off',...
    'CallBack',callBackstr);

% Cortical distance measurement- Binning method (use current line ROI)
callBackstr=['measureCorticalDistanceBins(' view.name ', getCurROIcoords(' view.name ') , 4);'];
uimenu(analysismenu,'Label','Measure cortical distance (bins ROI, binDist=4mm)','Separator','off',...
    'CallBack',callBackstr);

% Cortical distance measurement- Binning method (use current line ROI)
callBackstr=['measureCorticalDistanceBins(' view.name ', getCurROIcoords(' view.name '));'];
uimenu(analysismenu,'Label','Measure cortical distance (bins ROI, specify binDist)','Separator','off',...
    'CallBack',callBackstr);

% Current ROI surface area measurement
callBackstr=['measureFlatROIArea(' view.name ');'];
uimenu(analysismenu,'Label','ROI surface area (old, slow, inaccurate!)','Separator','on',...
    'CallBack',callBackstr);
callBackstr=[ view.name ' = measureFlatROIAreaMesh(' view.name ');'];
uimenu(analysismenu,'Label','ROI surface area (new mesh method)','Separator','off',...
    'CallBack',callBackstr);

%---------Atlas Version 2----------------------------------------------------

% atlasV2 = uimenu(analysismenu,'Label','Atlas (v. 2)','Separator','on');

% view = atlasCreate(atlasName, ringWedgeScans,view);
atlasCreate = uimenu(analysismenu,'Label','Create atlas','Separator','on');

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''Hemifield'',[],', view.name ');'];
uimenu(atlasCreate,'Label','Hemifield','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''uqf'',[],', view.name ');'];
uimenu(atlasCreate,'Label','upperquarterfield','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''lqf'',[],', view.name ');'];
uimenu(atlasCreate,'Label','lowerquarterfield','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''HemiLowerVF'',[],', view.name ');'];
uimenu(atlasCreate,'Label','HemiLowerVF','Separator','off','CallBack',callBackstr);

callBackstr = [ '[' view.name ' corners] = atlasCreate(''HemiUpperVF'',[],', view.name ');'];
uimenu(atlasCreate,'Label','HemiUpperVF (V3v/hV4)','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''V1/V2'',[],', view.name ');'];
uimenu(atlasCreate,'Label','V1/V2','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''V1/V2/V3'',[],', view.name ');'];
uimenu(atlasCreate,'Label','V1/V2/V3','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''ventralV2V3V4'',[],', view.name ');'];
uimenu(atlasCreate,'Label','V2v/V3v/V4v','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''V4/V8'',[],', view.name ');'];
uimenu(atlasCreate,'Label','V4/V8','Separator','off','CallBack',callBackstr);

% view = atlasCreate(atlasName, ringWedgeScans,view);
callBackstr = [ '[' view.name ' corners] = atlasCreate(''2hemifields'',[],', view.name ');'];
uimenu(atlasCreate,'Label','2 Hemifields (VO1/VO2)','Separator','off','CallBack',callBackstr);

callBackstr = [ '[' view.name ' corners] = atlasCreate(''4 hemifields'',[],', view.name ');'];
uimenu(atlasCreate,'Label','4 Hemifields (LO-1/2, TO-1/2)','Separator','off','CallBack',callBackstr);

%---------------------------------------------------------
% [view images] = makeAtlasData(view,[],[],corners);
callBackstr=[ '[' view.name  ' images] = makeAtlasData(' view.name ');'];
uimenu(analysismenu,'Label','Create data for the fitting GUI','Separator','off','CallBack',callBackstr);

% gui(images);
callBackstr=[ 'atlasfitGUI(images)'];
uimenu(analysismenu,'Label','Open the fitting GUI','Separator','off','CallBack',callBackstr);

% [images2 images3 images]= RawToCommonAtlas(images, corners, '/biac3/wandell4/data/MT/YM070815/safirData_left0706', '/biac3/wandell4/data/MT/YM070815/safirResult_left0706');
callBackstr=['images = loadPRFsize(' view.name ', images); [images images2 images3] = RawToCommonAtlas(images, corners, ''' pwd '/safirData'',''' pwd '/safirResult''); atlasStore(' view.name ',images)'];
uimenu(analysismenu,'Label','Create data type for the fitted atlas','Separator','off','CallBack',callBackstr);

% view = warpFromViews(view);
callBackstr=[ view.name ' = warpFromViews(' view.name ');'];
uimenu(analysismenu,'Label','Warp Atlas (Atlas view)','Separator','on','CallBack',callBackstr);

% Visual area ROI from atlas
% view = makeAreaROI(view);
callBackstr=[ view.name ' = makeAreaROI(' view.name ');'];
uimenu(analysismenu,'Label','Create Area ROI (Data view)','Separator','off','CallBack',callBackstr);

% view= atlasIsoAngleROIs(view);
callBackstr=[ view.name ' = atlasIsoAngleROIs(' view.name ');'];
uimenu(analysismenu,'Label','Iso-angle ROI (Data view)','Separator','off','CallBack',callBackstr);

% view= atlasIsoEccROIs(view);
callBackstr=[ view.name ' = atlasIsoEccROIs(' view.name ');'];
uimenu(analysismenu,'Label','Iso-eccentricity ROI (Data view)','Separator','off','CallBack',callBackstr);

% Atlas Error
callBackstr=[ view.name ' = atlasError(' view.name ');'];
uimenu(analysismenu,'Label','Compute error (data view)','Separator','off',...
    'CallBack',callBackstr);

%--------------------------------------------------------------------------
% Cortical Magnification computation
callBackstr=[ view.name ' = measureCortMag(' view.name ');'];
uimenu(analysismenu, 'Label', 'Compute cortical magnification', 'Separator', 'on',...
    'CallBack', callBackstr);
callBackstr=[ view.name ' = measureCortMag(' view.name ', [], uigetfile);'];
uimenu(analysismenu, 'Label', 'Compute cortMag from file', 'Separator', 'off',...
    'CallBack', callBackstr);

%------------Atlas V. 1 (Old)---------------------------------------------------------
atlasOldmenu = uimenu(analysismenu,'Label','Atlas (v. 1)','Separator','on');
% Create Rectangle ROI with default name callback
%   view=newROI(view);
%   view=addROIrect(view,1);
%   view=refreshScreen(view,1);

% 2-D Atlas fitting
callBackstr=['createAtlas(' view.name ');'];
uimenu(atlasOldmenu,'Label','Create Retinotopy Atlas (data view)','Separator','on',...
    'CallBack',callBackstr);

% 2-D Atlas phase adjustment
callBackstr=[ view.name ' = fixPhase(' view.name ');'];
uimenu(atlasOldmenu,'Label','Fix atlas phase (atlas view)','Separator','off',...
    'CallBack',callBackstr);

% 2-D Atlas affine transform
callBackstr=[ view.name ' = atlasAffineTransform(' view.name ');'];
uimenu(atlasOldmenu,'Label','Affine transform Atlas (atlas view)','Separator','off',...
    'CallBack',callBackstr);

% 2-D Atlas warping
callBackstr=[ view.name ' = warpFromViews(' view.name ');'];
uimenu(atlasOldmenu,'Label','Warp atlas (atlas view)','Separator','off',...
    'CallBack',callBackstr);
% Visual area ROI from atlas
callBackstr=[ view.name ' = makeAreaROI(' view.name ');'];
uimenu(atlasOldmenu,'Label','Create area ROI (data view)','Separator','off',...
    'CallBack',callBackstr);

% Iso-angle ROI from atlas
callBackstr=[ view.name ' = makeIsoAngleROIs(' view.name ');'];
uimenu(atlasOldmenu,'Label','Create iso-angle ROIs (data view)','Separator','off',...
    'CallBack',callBackstr);

% Iso-angle ROI from atlas
callBackstr=[ view.name ' = makeIsoEccenROIs(' view.name ');'];
uimenu(atlasOldmenu,'Label','Create iso-eccentricity ROIs (data view)','Separator','off',...
    'CallBack',callBackstr);

% Atlas Error Map
callBackstr=[ view.name ' = atlasErrorMap(' view.name ');'];
uimenu(atlasOldmenu,'Label','Compute error map (data view)','Separator','off',...
    'CallBack',callBackstr);

return;


text:
function [success] = omparameter (modelname,name,wert);
% Change parameter in *_init.txt file; Returns linenumber replaced if
% succeeded, -1 if not.
% SYNTAX: [linenumber] = omparameter(modelname,parameter,value)
% Beispiel: omparameter('package.model','stop value',10)
%
% Feedback/problems: Christian Schaad, ingenieurbuero@christian-schaad.de

success=-1;
inputfile=[modelname,'_init_original.xml'];
outputfile=[modelname,'_init.xml'];
%def_output=['outputFormat   = "mat"'];
%plt_output=['outputFormat   = "plt"'];

fid=fopen(inputfile);
i=0;
while 1
    i=i+1;
    tline0 = fgetl(fid);
    if ~ischar(tline0), break, end
    tline= [tline0,char(10)];
    %if ~isempty(strfind(tline,def_output));
    %tline=[plt_output,char(10)];
    %disp('MAT-Output changed to PLT-Output')
    %end
    if ~isempty(strfind(tline,['name = "',name,'"',char(10)]))
    for k=0:7
    dataset(i+k).string=tline;
    tline0 = fgetl(fid);
    tline= [tline0,char(10)];
    end
    i=i+k+1;

    tline=['<Real start="',wert,'" fixed="true" />',char(10)	];
    %disp(tline)
    success=i;
    end

    dataset(i).string=tline;
    clear tline tline0;

end
fclose(fid);

fid=fopen(outputfile,'w');
for j=1:i-1
fprintf(fid,[dataset(j).string]);
end
fclose(fid);
switch success
    case -1
        error(['Parameter ',name,' not found!']);
end;
disp(['Replaced line number ',num2str(success),', parameter ',num2str(name),' with ',num2str(wert) ]);

author: Sable

repo: mcbench-benchmarks


%  Script file: calc_decay.m
%
%  Purpose:
%    This program calculates the amount of Thorium 227 and
%    Radium 223 left as a function of time, given an inital
%    concentration of 1 gram of Thorium 227 and no grams of
%    Radium 223.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/15/07    S. J. Chapman        Original code
%
% Define variables:
%   odefun_handle -- Handle to function that defines the derivative
%   tspan         -- Duration to solve equation for
%   yo            -- Initial condition for equation
%   t             -- Array of solution times
%   y             -- Array of solution values

% Get a handle to the function that defines the derivative.
odefun_handle = @decay1;

% Solve the equation over the period 0 to 100 days
tspan = [0 100];

% Set the initial conditions
y0(1) = 1000000;     % Atoms of Thorium 227
y0(2) = 0;           % Atoms of Radium 223

% Call the differential equation solver.
[t,y] = ode45(odefun_handle,tspan,y0);

% Plot the result
figure(1);
plot(t,y(:,1),'b-','LineWidth',2);
hold on;
plot(t,y(:,2),'k--','LineWidth',2);
title('\bfAmount of Thorium 227 and Radium 223 vs Time');
xlabel('\bfTime (days)');
ylabel('\bfNumber of Atoms');
legend('Thorium 227','Radium 223');
grid on;
hold off;


text:
function [ function_value] = SimpleQPObjective( values, Problem)
%SimpleQPObjective Calculates obj.*sum(x.^2)
%   Where obj is provided as objArguments{1} of the Problem struct.

obj = Problem.objArguments{1};

function_value = sum(obj.*values.^2);
end

author: opencobra

repo: cobratoolbox


text:
function tSeries = meanTSeriesForPlotting(vw, selectedROIs, getRawData)
% Return the mean t-series for one or more ROIs. Start by checking for data
% in Inplane. If this fails, try the gray view.
%
%   tSeries = meanTSeriesForPlotting(vw, selectedROIs)
%
%   Example:
%
%       tSeries = meanTSeriesForPlotting(vw);
%
%   See: plotMultipleSingleCycleErr.m
%
%   JW, 8/2009: Split off from various plotting functions that each
%   duplicated this code
%   JW, 9/2010: Added flag 'getRawData' in case you want to plot without
%               converting to % signal change, replaced GLOBAL INPLANE with
%               ip
mrGlobals;

viewType = viewGet(vw, 'viewType');

if ~exist('selectedROIs', 'var') || isempty(selectedROIs),
    selectedROIs = roiGetList(vw);
end

if ~exist('getRawData', 'var') || isempty(getRawData),
    getRawData  = false;
end

nROIs = length(selectedROIs);
scan  = viewGet(vw,'currentScan');

switch viewType
    case {'Volume' 'Gray'}   %%%For ROIs in Gray view - xform to inplane

        selectedVOLUME = viewSelected('volume');

        %initiate hiddenINPLANE
        ip = initHiddenInplane;
        ip = viewSet(ip,'name','hidden');

        % Set the Inplane scan number and datatype to match the Volume view.
        curDataType = viewGet(vw,'datatypenumber');

        ip = viewSet(ip,'datatypenumber',curDataType);
        ip = viewSet(ip,'currentscan',scan);

        %Transfer current VOLUME ROI to INPLANE
        for ii = 1:nROIs
            vw = viewSet(vw,'selectedROI', selectedROIs(ii));
            ip = vol2ipCurROI(vw,ip);
            selectedInplaneROIs(ii) = viewGet(ip, 'curroi');
        end

    case {'Flat'} %%%For ROIs in Flat view - xform to inplane

        %initiate and / or select VOLUME and INPLANE windows
        if isempty(VOLUME),
            VOLUME{1} = initHiddenGray;
            VOLUME{1} = viewSet(VOLUME{1},'name','hidden');
        else
            selectedVOLUME = viewSelected('volume');
        end

        ip = initHiddenInplane;
        ip = viewSet(ip,'name','hidden');

        % Set the Inplane scan number and datatype to match the Flat view.
        %         curScan =     viewGet(FLAT{selectedFLAT},'currentscan');
        curDataType = viewGet(vw,'datatypenumber');
        ip = viewSet(ip,'datatypenumber',curDataType);
        ip = viewSet(ip,'currentscan',scan);

        %Transfer current FLAT ROI to INPLANE
        for ii = 1: nROIs
            vw = viewSet(vw,'selectedROI', selectedROIs(ii));
            ip = flat2ipCurROI(vw,ip,VOLUME{selectedVOLUME});
            selectedInplaneROIs(ii) = viewGet(ip, 'curroi');
        end

    case {'Inplane'}   %%%For ROIs in INPLANE view - select inplane
        selectedInplaneROIs = selectedROIs;
        ip = vw;
end

% Compute meanTSeries for each ROI
ROIcoords = cell(1,nROIs);
try
    for r=1:nROIs
        ROIcoords{r}=ip.ROIs(selectedInplaneROIs(r)).coords;
    end
    tSeries = meanTSeries(ip,scan,ROIcoords, getRawData);
catch ME
    warning(ME.identifier, ME.message)
    for r=1:nROIs
        ROIcoords{r}=VOLUME{selectedVOLUME}.ROIs(selectedROIs(r)).coords;
    end
    tSeries = meanTSeries(VOLUME{selectedVOLUME},scan,ROIcoords, getRawData);
end

if ~iscell(tSeries)
    tmp{1}=tSeries;
    tSeries=tmp;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% MULTIGRID OF FOR THE STOKES EQNS IN 2D
%
% This example is to show the convergence of multigrid methods for various
% finite element approximation of the Stokes equation on the unit square:
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%                              u = g_D  on \Gamma.
%
% with the pure Dirichlet boundary condition.
%

%% Setting
clear all; close all;
% [node,elem] = squaremesh([0,1,0,1],0.5);
[node,elem] = circlemesh(0,0,1,0.25);
% [node,elem] = uniformrefine(node,elem);
showmesh(node,elem);
pde = Stokesdata1;
%  pde = StokesZulehnerdata;
bdFlag = setboundary(node,elem,'Dirichlet');
option.L0 = 1;
option.maxIt = 3;
option.printlevel = 1;
option.plotflag = 0;
option.rateflag = 0;

%% MG options
option.solver = 'asmg';
option.smoothingstep = 2;
option.smootherbarSp = 'SGS';

%% CR-P0 element
display('CR-P0')
option.elemType = 'CRP0';
femStokes(node,elem,pde,bdFlag,option);

%% P2-P0 element
display('P2-P0')
option.elemType = 'P2P0';
femStokes(node,elem,pde,bdFlag,option);

%% P2-P1 element
display('P2-P1')
option.elemType = 'P2P1';
option.solver = 'asmg';
femStokes(node,elem,pde,bdFlag,option);

%% isoP2-P0 element
display('isoP2-P0')
option.elemType = 'isoP2P0';
femStokes(node,elem,pde,bdFlag,option);

%% isoP2-P1 element
display('isoP2-P1')
option.elemType = 'isoP2P1';
femStokes(node,elem,pde,bdFlag,option);

%% P1b-P1 element
display('P1b-P0')
option.elemType = 'P1bP1';
femStokes(node,elem,pde,bdFlag,option);

%% Conclusion
%

repo: Aero-matlab


text:
function value = cos_deg ( angle )

%*****************************************************************************80
%
%% COS_DEG returns the cosine of an angle given in degrees.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ANGLE, the angle, in degrees.
%
%    Output, real COS_DEG, the cosine of the angle.
%
  angle_rad = pi * angle / 180.0;

  value = cos ( angle_rad );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [U,G,I,J] = remove_small_components(V,F,varargin)
  % REMOVE_SMALL_COMPONENTS
  %
  % [U,G,I,J] = remove_small_components(V,F)
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of face indices into rows of V
  %   Optional:
  %     'Tol' followed by minimum volume (or surface area) to use OR the word
  %       'max' to keep only biggest component {0.0001*total_vol}
  %     'Volumetric' followed by whether to use volume or surface area
  % Outputs:
  %   U  #U by 3 list of vertex positions
  %   G  #G by 3 list of face indices into rows of U
  %   I  #V by 1 list of indices such that: G = I(F)
  %   J  #G by 1 list of indices into F
  %

  tol = [];
  volumetric = true;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Tol','Volumetric'}, {'tol','volumetric'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if isempty(tol)
    if volumetric
      [~,total] = centroid(V,F);
    else
      total = sum(doublearea(V,F))*0.5;
    end
    tol = 0.0001*total;
  end

  [~,C] = connected_components(F);
  nc = max(C);
  val = zeros(nc,1);
  for i = 1:nc
    Fi = F(C==i,:);
    if volumetric
      [~,val(i)] = centroid(V,Fi);
    else
      val(i) = 0.5*sum(doublearea(V,Fi));
    end
  end

  if isnumeric(tol)
    J = find(ismember(C,find(val>tol)));
  else
    assert(strcmp(tol,'max'));
    [~,max_c] = max(val);
    J = find(C==max_c);
  end
  F = F(J,:);
  [U,I] = remove_unreferenced(V,F);
  G = I(F);
end

author: alecjacobson

repo: gptoolbox


text:
function x = p39_start ( n )

%*****************************************************************************80
%
%% P39_START returns a starting point for optimization for problem 39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ 0.6, 1.3 ]';

  return
end


text:
function [X_den,P]=denoise_TV_One(Xobs,lambda,l,u,P_init,pars)
%This function implements the FISTA method for TV denoising problems.
%
% INPUT
% Xobs ..............................an observed noisy image.
% lambda ........................ parameter
% pars.................................parameters structure
% pars.MAXITER ..................... maximum number of iterations
%                                                      (Default=100)
% pars.epsilon ..................... tolerance for relative error used in
%                                                       the stopping criteria (Default=1e-4)
% pars.print ..........................  1 if a report on the iterations is
%                                                       given, 0 if the  report is silenced
% pars.tv .................................. type of total variation
%                                                      penatly.  'iso' for isotropic (default)
%                                                      and 'l1' for nonisotropic
%
% OUTPUT
% X_den ........................... The solution of the problem
%                                            min{||X-Xobs||^2+2*lambda*TV(X)}

%Define the Projection onto the box
if((l==-Inf)&&(u==Inf))
    project=@(x)x;
elseif (isfinite(l)&&(u==Inf))
    project=@(x)(((l<x).*x)+(l*(x<=l)));
elseif (isfinite(u)&&(l==-Inf))
     project=@(x)(((x<u).*x)+((x>=u)*u));
elseif ((isfinite(u)&&isfinite(l))&&(l<u))
    project=@(x)(((l<x)&(x<u)).*x)+((x>=u)*u)+(l*(x<=l));
else
    error('lower and upper bound l,u should satisfy l<u');
end

% Assigning parameres according to pars and/or default values
flag=exist('pars', 'var');
if (flag&&isfield(pars,'MAXITER'))
    MAXITER=pars.MAXITER;
else
    MAXITER=100;
end
% if (flag&&isfield(pars,'epsilon'))
%     epsilon=pars.epsilon;
% else
%     epsilon=1e-4;
% end
% if(flag&&isfield(pars,'print'))
%     prnt=pars.print;
% else
%     prnt=1;
% end
if(flag&&isfield(pars,'tv'))
    tv=pars.tv;
else
    tv='iso';
end

[m,n]=size(Xobs);
% clear P; clear R;
if(isempty(P_init))
    P{1}=zeros(m-1,n);    P{2}=zeros(m,n-1);
    R{1}=zeros(m-1,n);    R{2}=zeros(m,n-1);
else
    P{1}=P_init{1};    P{2}=P_init{2};
    R{1}=P_init{1};    R{2}=P_init{2};
end
tk=1;tkp1=1;count=0;i=0;

D=zeros(m,n);%fval=inf;fun_all=[];
while((i<MAXITER)&&(count<5))
%    fold=fval;
    i=i+1;
%     Dold=D;
    Pold=P;
    tk=tkp1;
    D=project(Xobs-lambda*Lforward(R, m, n));
    Q=Ltrans(D, m, n);
    %%%%%%%%%%
    % Taking a step towards minus of the gradient
    P{1}=R{1}+1/(8*lambda)*Q{1};
    P{2}=R{2}+1/(8*lambda)*Q{2};

    %%%%%%%%%%
    % Peforming the projection step
    switch tv
        case 'iso'
            A=[P{1}.^2;zeros(1,n)]+[P{2}.^2,zeros(m,1)];
            A=sqrt(max(A,1));
            P{1}=P{1}./A(1:m-1,:); P{2}=P{2}./A(:,1:n-1);
        case 'l1'
            P{1}=P{1}./(max(abs(P{1}),1));
            P{2}=P{2}./(max(abs(P{2}),1));
        otherwise
            error('unknown type of total variation. should be iso or l1');
    end

    %%%%%%%%%%
    %Updating R and t
    tkp1=(1+sqrt(1+4*tk^2))/2;

    R{1}=P{1}+(tk-1)/(tkp1)*(P{1}-Pold{1});
    R{2}=P{2}+(tk-1)/tkp1*(P{2}-Pold{2});

%     re=norm(D-Dold,'fro')/norm(D,'fro');
%     if (re<epsilon)
%         count=count+1;
%     else
%         count=0;
%     end
    C=Xobs-lambda*Lforward(P, m, n);
    PC=project(C);
%     fval=-norm(C-PC,'fro')^2+norm(C,'fro')^2;
%     fun_all=[fun_all;fval];
end
X_den=D;iter=i;

function X=Lforward(P, m, n)

%       [m2,n2]=size(P{1});
%       [m1,n1]=size(P{2});
%
%       if (n2~=n1+1)
%           error('dimensions are not consistent')
%       end
%       if(m1~=m2+1)
%           error('dimensions are not consistent')
%       end
%
%       m=m2+1;
%       n=n2;

      X=zeros(m,n);
      X(1:m-1,:)=P{1};
      X(:,1:n-1)=X(:,1:n-1)+P{2};
      X(2:m,:)=X(2:m,:)-P{1};
      X(:,2:n)=X(:,2:n)-P{2};
   end

   function P=Ltrans(X, m, n)
%       [m,n]=size(X);
      P{1}=X(1:m-1,:)-X(2:m,:);
      P{2}=X(:,1:n-1)-X(:,2:n);
   end
end

author: thomaskuestner


text:
function indx = r8vec_indexed_heap_d ( n, a, indx )

%*****************************************************************************80
%
%% R8VEC_INDEXED_HEAP_D creates a descending heap from an indexed R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    An indexed R8VEC is an R8VEC of data values, and an R8VEC of N indices,
%    each referencing an entry of the data vector.
%
%    The function adjusts the index vector INDX so that, for 1 <= J <= N/2,
%    we have:
%      A(INDX(2*J))   <= A(INDX(J))
%    and
%      A(INDX(2*J+1)) <= A(INDX(J))
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms for Computers and Calculators,
%    Academic Press, 1978,
%    ISBN: 0-12-519260-6,
%    LC: QA164.N54.
%
%  Parameters:
%
%    Input, integer N, the size of the index array.
%
%    Input, real A(*), the data vector.
%
%    Input, integer INDX(N), the index array.
%    Each entry of INDX must be a valid index for the array A.
%
%    Output, integer INDX(N), the indices have been reordered into a
%    descending heap.
%

%
%  Only nodes N/2 down to 1 can be "parent" nodes.
%
  for i = floor ( n / 2 ) : -1 : 1
%
%  Copy the value out of the parent node.
%  Position IFREE is now "open".
%
    key = indx(i);
    ifree = i;

    while ( 1 )
%
%  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
%  IFREE.  (One or both may not exist because they exceed N.)
%
      m = 2 * ifree;
%
%  Does the first position exist?
%
      if ( n < m )
        break
      end
%
%  Does the second position exist?
%
      if ( m + 1 <= n )
%
%  If both positions exist, take the larger of the two values,
%  and update M if necessary.
%
        if ( a(indx(m)) < a(indx(m+1)) )
          m = m + 1;
        end

      end
%
%  If the large descendant is larger than KEY, move it up,
%  and update IFREE, the location of the free position, and
%  consider the descendants of THIS position.
%
      if ( a(indx(m)) <= a(key) )
        break
      end

      indx(ifree) = indx(m);
      ifree = m;

    end
%
%  Once there is no more shifting to do, KEY moves into the free spot IFREE.
%
    indx(ifree) = key;

  end

  return
end


text:
% relperm and capillary pressure curves
% to be used with FVTool package
% Written by Ali A. Eftekhari
% This is what this function does:
% res=zeros(size(sw));
% for i=1:numel(sw)
%   if swc(i)<=sw(i) && sw(i)<=1-sor(i)
%     res(i)=kro0(i)*((1-sw(i)-sor(i))/(1-sor(i)-swc(i)))^no(i);
%   elseif 0.0<sw(i) && sw(i)<swc(i)
%     res(i)=1+(kro0(i)-1)/swc(i)*sw(i);
%   elseif sw(i)>1-sor(i)
%     res(i)=0.0;
%   elseif sw(i)<=0.0
%     res(i)=1.0;
%   end
% end
%
% However, it is slow in this form. Therefore, I have rewritten it in
% vectorized from. This one is kept here for reference. You can accelerate
% it by @parfor loop, but the vectorized form is substantially faster
function res=kro_loop(sw, kro0, sor, swc, no)
res=zeros(size(sw));
for i=1:numel(sw)
  if swc(i)<=sw(i) && sw(i)<=1-sor(i)
    res(i)=kro0(i)*((1-sw(i)-sor(i))/(1-sor(i)-swc(i)))^no(i);
  elseif 0.0<sw(i) && sw(i)<swc(i)
    res(i)=1+(kro0(i)-1)/swc(i)*sw(i);
  elseif sw(i)>1-sor(i)
    res(i)=0.0;
  elseif sw(i)<=0.0
    res(i)=1.0;
  end
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function p = addComplementarityCuts(p)

if size(p.complementary,1)>0
    p_cut = emptyNumericalModel;
    for i = 1:size(p.complementary,1)
        variables = p.complementary(i,:);
        if p.lb(variables(1))>0 & p.lb(variables(2))>0
            p_cut.F_struc(end+1,1)=p.ub(variables(2));
            p_cut.F_struc(end,1+variables(1))=-p.ub(variables(2))/p.ub(variables(1));
            p_cut.F_struc(end,1+variables(2))=-1;
            p_cut.K.l = p_cut.K.l + 1;
        end
    end
    p = mergeNumericalModels(p,p_cut);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
%
% This file is part of HRVAS
%
% HRVAS is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% HRVAS 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 HRVAS.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function batchHRV(opt)
%batchHRV: function that calculates hrv on all ibi files in selected
%           directory
%
%Inputs:
%
%Outputs:
%

%check inputs
% if (nargin < 1) || ~isstruct(opt)
%
% end

    %choose file containing HRV analysis options
    reply = input('Enter full file path of HRV analysis option file: ', 's');
    if isempty(reply) || exist(reply,'file')~=2
        error('Please choose valid file path!')
        return;
    end
    opt=load(reply);
    opt=opt.settings;

    %choose dir containing input files
    reply = input('Enter dir path containing IBI files: ', 's');
    if isempty(reply) || exist(reply,'dir')~=7
        error('Please choose valid input dir path!')
        return;
    end

    %check input directory
    fileList = dir(fullfile(reply, '*.ibi')); %get list of files
    fileList(any([fileList.isdir],1))=[]; %remove folders/dir from the list
    fnames = {fileList.name}; %get file names only from fileList structure
    %build array of full file paths
    fpaths=cell(length(fnames),1);
    for ff=1:length(fnames)
        fpaths{ff}=fullfile(reply,fnames{ff});
    end

    %display results of input dir
    disp([reply ' contains ' num2str(length(fpaths)) ' .ibi files.'])

    %choose output file name to save hrv reslts
    outPath = input('Enter file name for exported HRV data: ', 's');
    if isempty(outPath)
        error('Please choose valid input dir path!')
        return;
    end

    %Batch process HRV
    hrv=batchGetHRV(fpaths,opt);

    %Export/save HRV
    saveHRV(hrv,opt,fpaths,outPath)
    disp('batchHRV Done!')

end

function saveHRV(hrv,opt,fList,outPath)
       exportHRV(outPath,fList,hrv,opt)
end

function hrv=batchGetHRV(fList,opt)
    %Preallocate hrv array
    %hrv=repmat(struct('b',0,'a',0),1,length(fList));

    %display time remaining
    disp('Processing: 0% complete.')
    tic; avgElapsed=0;

    %get hrv for each file
    nFiles=length(fList);
    for f=1:nFiles
        hrv(f)=getHRV(fList{f},opt);

        %calculate remaining time for waitbar
        elapsedTime=toc; %total time elapsed since start
        %average time elapsed to process a single file
        avgElapsed=elapsedTime/f;
        %time remaining in processing (min)
        timeRem=ceil(avgElapsed*(nFiles-f)/60);
        bar=round(f/nFiles*100);
        disp(['Processing: ' num2str(bar) ...
            '% complete. (~' num2str(timeRem) ' min)'])
    end

end

function output=getHRV(f,settings)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % LOAD IBI
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nIBI=[]; dIBI=[];
    IBI=loadIBI(f,settings);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Preprocess Data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %build cell array of locate artifacts methods
    methods={}; methInput=[];
    if settings.ArtLocatePer
        methods=[methods,'percent'];
        methInput=[methInput,settings.ArtLocatePerVal];
    end
    if settings.ArtLocateSD
        methods=[methods,'sd'];
        methInput=[methInput,settings.ArtLocateSDVal];
    end
    if settings.ArtLocateMed
        methods=[methods,'median'];
        methInput=[methInput,settings.ArtLocateMedVal];
    end
    %determine which window/span to use
    if strcmpi(settings.ArtReplace,'mean')
        replaceWin=settings.ArtReplaceMeanVal;
    elseif strcmpi(settings.ArtReplace,'median')
        replaceWin=settings.ArtReplaceMedVal;
    else
        replaceWin=0;
    end

    %Note: We don't need to use all the input arguments,but we will let the
    %function handle all inputs
    [dIBI,nIBI,trend,art] = preProcessIBI(IBI, ...
        'locateMethod', methods, 'locateInput', methInput, ...
        'replaceMethod', settings.ArtReplace, ...
        'replaceInput',replaceWin, ...
        'detrendMethod', settings.Detrend, ...
        'smoothMethod', settings.SmoothMethod, ...
        'smoothSpan', settings.SmoothSpan, ...
        'smoothDegree', settings.SmoothDegree, ...
        'polyOrder', settings.PolyOrder, ...
        'waveletType', ...
        [settings.WaveletType num2str(settings.WaveletType2)], ...
        'waveletLevels', settings.WaveletLevels, ...
        'lambda', settings.PriorsLambda,...
        'resampleRate',settings.Interp,...
        'meanCorrection',true);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate HRV
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    output.ibiinfo.count=size(IBI,1); % total # of ibi
    output.ibiinfo.outliers=sum(art); % number of outliers

    %Time-Domain (using non-detrented ibi)
    output.time=timeDomainHRV(nIBI,settings.SDNNi*60,settings.pNNx);
    %output.time.mean=round(mean(nIBI(:,2).*1000)*10)/10;
    %output.time.meanHR=round(mean(60./nIBI(:,2))*10)/10;

    %Freq-Domain
    output.freq=freqDomainHRV(dIBI,settings.VLF,settings.LF, ...
        settings.HF,settings.AROrder,settings.WinWidth, ...
        settings.WinOverlap,settings.Points,settings.Interp);

    %Nonlinear (using non-detrented ibi)
    output.nl=nonlinearHRV(nIBI,settings.m,settings.r, ...
        settings.n1,settings.n2,settings.breakpoint);
    %Poincare
    output.poincare=poincareHRV(nIBI);

    %Time-Freq
    output.tf=timeFreqHRV(dIBI,nIBI,settings.VLF,settings.LF, ...
        settings.HF,settings.AROrder, settings.tfWinSize, ...
        settings.tfOverlap,settings.Points,settings.Interp, ...
        {'ar','lomb','wavelet'});
    %%%%%%%%%%%%%%%%%%%%%%%

    clear  nIBI dIBI IBI

end

function ibi=loadIBI(f,opt)
    if ~exist(f,'file')
        error(['Error opening file: ' f])
        return
    end

    ibi=[];
    DELIMITER = ',';
    HEADERLINES = opt.headerSize;

    %read ibi
    tmpData = importdata(f, DELIMITER, HEADERLINES);
    if HEADERLINES>0
        tmpData=tmpData.data;
    end

    %check ibi dimentions
    [rows cols] = size(tmpData);
    if rows==1
        tmpData=tmpData';
        ibi=zeros(cols,2);
        tmp=cumsum(tmpData);
        ibi(2:end,1)=tmp(1:end-1);
        ibi(:,2)=tmpData;
    elseif cols==1
        ibi=zeros(rows,2);
        tmp=cumsum(tmpData);
        ibi(2:end,1)=tmp(1:end-1);
        ibi(:,2)=tmpData;
    elseif rows<cols
        ibi=tmpData';
    else
        ibi=tmpData;
    end
    clear tmpData
end

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


text:
function pass = test_scl( pref )
% Check correct vertical scaling.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 10*pref.cheb2Prefs.chebfun2eps;
j = 1;

% Scale invariant
f = chebfun2( @(x,y) cos(x.*y) );
g = chebfun2( @(x,y) eps*cos(x.*y) );
err = norm( eps*f - g);
pass(j) = ( err  <  tol ); j = j + 1;

% hscale invariant
f = chebfun2( @(x,y) cos(x.*y) );
g = chebfun2( @(x,y) cos(x/eps.*y/eps), eps*[-1 1 -1 1] );
err = abs( f(1,1) - g(eps,eps) );
pass(j) = ( err  < tol ); j = j + 1;
err = abs( f(pi/6,1) - g(eps*pi/6,eps) );
pass(j) = ( err  < tol ); j = j + 1;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [montage, cfg] = ft_prepare_montage(cfg, data)

% FT_PREPARE_MONTAGE creates a referencing scheme based on the input configuration
% options and the channels in the data structure. The resulting montage can be
% given as input to FT_APPLY_MONTAGE, or as cfg.montage to FT_PREPROCESSING.
%
% Use as
%   montage = ft_prepare_montage(cfg, data)
%
% The configuration can contain the following fields:
%   cfg.refmethod     = 'avg', 'comp', 'bipolar', 'laplace', 'doublebanana', 'longitudinal', 'circumferential', 'transverse' (default = 'avg')
%   cfg.implicitref   = string with the label of the implicit reference, or empty (default = [])
%   cfg.refchannel    = cell-array with new EEG reference channel(s), this can be 'all' for a common average reference
%   cfg.groupchans    = 'yes' or 'no', should channels be rereferenced in separate groups
%                       for bipolar and laplace methods, this requires channnels to be
%                       named using an alphanumeric code, where letters represent the
%                       group and numbers represent the order of the channel whithin
%                       its group (default = 'no')
%
% The implicitref option allows adding the implicit reference channel to the data as
% a channel with zeros.
%
% The resulting montage is a structure with the fields
%   montage.tra      = MxN matrix
%   montage.labelold = Nx1 cell-array
%   montage.labelnew = Mx1 cell-array
%
% As an example, an output bipolar montage could look like this
%   bipolar.labelold  = {'1',   '2',   '3',   '4'}
%   bipolar.labelnew  = {'1-2', '2-3', '3-4'}
%   bipolar.tra       = [
%     +1 -1  0  0
%      0 +1 -1  0
%      0  0 +1 -1
%   ];
%
% See also FT_PREPROCESSING, FT_APPLY_MONTAGE

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble loadvar    data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the data can be passed as input argument or can be read from disk
hasdata = exist('data', 'var') && ~isempty(data);

% basic check/initialization of input arguments
if ~hasdata
  data = struct([]);
else
  data = ft_checkdata(data, 'datatype', {'raw', 'raw+comp', 'timelock', 'timelock+comp'});
end

% do a sanity check for incompatible options which are used in ft_preprocessing and elsewhere
cfg = ft_checkconfig(cfg, 'forbidden', {'montage', 'method'});

% set default configuration options
cfg.refmethod    = ft_getopt(cfg, 'refmethod', 'avg');
cfg.channel      = ft_getopt(cfg, 'channel', 'all');
cfg.implicitref  = ft_getopt(cfg, 'implicitref');
cfg.refchannel   = ft_getopt(cfg, 'refchannel');
cfg.groupchans   = ft_getopt(cfg, 'groupchans', 'no');

if isfield(cfg, 'reref')
  assert(istrue(cfg.reref), 'cannot create a montage without cfg.reref=''yes''');
end

% here the actual work starts
if hasdata
  cfg.channel    = ft_channelselection(cfg.channel, data.label);
  cfg.refchannel = ft_channelselection(cfg.refchannel, cat(1, data.label(:), cfg.implicitref));
else
  if ischar(cfg.channel)
    cfg.channel = {cfg.channel};
  end
  if ischar(cfg.refchannel)
    cfg.refchannel = {cfg.refchannel};
  end
end

% the first montage serves to add the implicit reference channel
% the last row for the implicitref will be all zero
montage1 = [];
montage1.labelold = cfg.channel(:);
montage1.labelnew = cfg.channel(:);
if ~isempty(cfg.implicitref)
  montage1.labelnew{end+1} = cfg.implicitref;
end
montage1.tra = eye(numel(montage1.labelnew), numel(montage1.labelold));

% the second montage serves to subtract the selected reference channels
switch cfg.refmethod
  case 'avg'
    % make a montage that subtracts the mean of all channels specified in cfg.refchannel
    montage2 = [];
    montage2.labelold = montage1.labelnew;
    montage2.labelnew = montage1.labelnew;
    montage2.tra      = eye(numel(montage2.labelnew));
    refsel = match_str(montage2.labelold, cfg.refchannel);
    montage2.tra(:,refsel) = montage2.tra(:,refsel) - 1/numel(refsel);

    % apply montage2 to montage1, the result is the combination of both
    montage = ft_apply_montage(montage1, montage2);

  case 'comp'
    assert(isempty(cfg.implicitref), 'implicitref not supported with refmethod=''%s''', cfg.refmethod);
    % construct an linear projection from channels to components
    montage = [];
    montage.labelold = data.topolabel;
    montage.labelnew = data.label;
    montage.tra      = data.unmixing;

  case 'invcomp'
    assert(isempty(cfg.implicitref), 'implicitref not supported with refmethod=''%s''', cfg.refmethod);
    % construct an linear projection from components to channels
    montage = [];
    montage.labelold = data.label;
    montage.labelnew = data.topolabel;
    montage.tra      = data.topo;

  case 'bipolar'
    if istrue(cfg.groupchans)
      % this uses recursion to make the montage for each group
      label = montage1.labelnew;
      group = cell(size(label));
      for i=1:numel(label)
        group(i) = regexp(label{i}, '^[a-zA-Z]+', 'match');
      end
      ugroup = unique(group, 'stable');

      tmpcfg = keepfields(cfg, {'refmethod', 'implicitref', 'refchannel', 'channel'});
      tmpcfg.groupchans = 'no';
      tmpcfg.showcallinfo = 'no';

      groupmontage = cell(size(ugroup));
      for g = 1:length(ugroup) % for each group of channels
        % select the channels within the group
        tmpcfg.channel = label(strcmp(group, ugroup{g}));
        groupmontage{g} = ft_prepare_montage(tmpcfg);
      end

      % start with a montage that keeps all channels the same
      montage2 = [];
      montage2.labelold = label(:)';
      montage2.labelnew = label(:);
      montage2.tra = eye(length(label));

      for g = length(ugroup):-1:1
        % apply the montage for each groep (i.e., combine them), keep all other channels as they are
        % unused channels end up at the end, hence doing it in reverse order keeps the channel order consistent
        montage2 = ft_apply_montage(montage2, groupmontage{g}, 'keepunused', true);
      end

    else
      % create a montage for the bipolar derivation of sequential channels
      montage2          = [];
      montage2.labelold = montage1.labelnew;
      montage2.labelnew = strcat(montage1.labelnew(1:end-1),'-',montage1.labelnew(2:end));
      tra_neg           = diag(-ones(numel(montage1.labelnew)-1,1), 1);
      tra_plus          = diag( ones(numel(montage1.labelnew)-1,1),-1);
      montage2.tra      = tra_neg(1:end-1,:)+tra_plus(2:end,:);

    end % if groupchans

    % apply montage2 to montage1, the result is the combination of both
    montage = ft_apply_montage(montage1, montage2);

  case 'laplace'
    if istrue(cfg.groupchans)
      % this uses recursion to make the montage for each group
      label = montage1.labelnew;
      group = cell(size(label));
      for i=1:numel(label)
        group(i) = regexp(label{i}, '^[a-zA-Z]+', 'match');
      end
      group = unique(group, 'stable');

      tmpcfg = keepfields(cfg, {'refmethod', 'implicitref', 'refchannel', 'channel'});
      tmpcfg.groupchans = 'no';
      tmpcfg.showcallinfo = 'no';

      groupmontage = cell(size(group));
      for g = 1:length(group) % for each group of channels
        % select the channels within the group
        tmpcfg.channel = label(startsWith(label, group{g}));
        groupmontage{g} = ft_prepare_montage(tmpcfg);
      end

      % start with a montage that keeps all channels the same
      montage2 = [];
      montage2.labelold = label(:)';
      montage2.labelnew = label(:);
      montage2.tra = eye(length(label));

      for g = length(group):-1:1
        % apply the montage for each groep (i.e., combine them), keep all other channels as they are
        % unused channels end up at the end, hence doing it in reverse order keeps the channel order consistent
        montage2 = ft_apply_montage(montage2, groupmontage{g}, 'keepunused', true);
      end

    else
      % create a montage for the laplacian derivation of neighboring channels where
      % each channel is re-referenced against the mean of the two nearest channels,
      % while channels at the extremities will be re-referenced against their closest
      % neighbor; this is particularly useful for sEEG shafts

      montage2          = [];
      montage2.labelold = montage1.labelnew;
      tra_neg1          = diag(-0.5*ones(numel(montage1.labelnew)-1,1), 1);
      tra_neg2          = diag(-0.5*ones(numel(montage1.labelnew)-1,1), -1);

      if numel(montage2.labelold) > 2
        montage2.labelnew = montage1.labelnew;
        tra_plus          = diag(ones(numel(montage1.labelnew),1),0);
        montage2.tra      = tra_neg1+tra_neg2+tra_plus;
        montage2.tra(1,2) = -1;
        montage2.tra(end,end-1) = -1;
      else
        % The following commands are just meant to avoid errors, in case
        % only one channel is given as input: the algorithm will behave
        % exactly as it would in the 'bipolar' case
        montage2.labelnew = cell(numel(montage1.labelnew),0);
        tra_plus          = diag(ones(numel(montage1.labelnew)-1,1),-1);
        montage2.tra      = tra_neg1+tra_neg2+tra_plus(2:end,:);
      end

    end % if groupchans

    % apply montage2 to montage1, the result is the combination of both
    montage = ft_apply_montage(montage1, montage2);

  case {'longitudinal', 'doublebanana'}
    % see https://www.learningeeg.com/montages-and-technical-components
    montage = [];
    montage.labelold = {'Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'Oz', 'O2'};
    montage.labelnew = {
      % left temporal chain
      'Fp1-F7'
      'F7-T3'
      'T3-T5'
      'T5-O1'
      % left parasagittal chain
      'Fp1-F3'
      'F3-C3'
      'C3-P3'
      'P3-O1'
      % central chain
      'Fz-Cz'
      'Cz-Pz'
      % right parasagittal chain
      'Fp2-F4'
      'F4-C4'
      'C4-P4'
      'P4-O2'
      % right temporal chain
      'Fp2-F8'
      'F8-T4'
      'T4-T6'
      'T6-O2'
      };
    % construct the montage from new channel labels
    montage.tra = zeros(length(montage.labelnew), length(montage.labelold));
    for i=1:length(montage.labelnew)
      lab = split(montage.labelnew{i}, '-');
      montage.tra(i, strcmp(montage.labelold, lab{1})) = +1;
      montage.tra(i, strcmp(montage.labelold, lab{2})) = -1;
    end
    % do a sanity check on the montage
    assert(all(sum(montage.tra, 2)==0));
    assert(all(max(montage.tra, [], 2)==+1));
    assert(all(min(montage.tra, [], 2)==-1));

  case 'circumferential'
    % see https://www.learningeeg.com/montages-and-technical-components
    montage = [];
    montage.labelold = {'Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'Oz', 'O2'};
    montage.labelnew = {
      'Fp1-F7'
      'F7-T3'
      'T3-T5'
      'T5-O1'
      'O1-O2'
      'O2-T6'
      'T6-T4'
      'T4-F8'
      'F8-Fp2'
      'Fp2-Fp1'
      };
    % construct the montage from new channel labels
    montage.tra = zeros(length(montage.labelnew), length(montage.labelold));
    for i=1:length(montage.labelnew)
      lab = split(montage.labelnew{i}, '-');
      montage.tra(i, strcmp(montage.labelold, lab{1})) = +1;
      montage.tra(i, strcmp(montage.labelold, lab{2})) = -1;
    end
    % do a sanity check on the montage
    assert(all(sum(montage.tra, 2)==0));
    assert(all(max(montage.tra, [], 2)==+1));
    assert(all(min(montage.tra, [], 2)==-1));

  case 'transverse'
    % this is inspired by https://doi.org/10.1016/j.earlhumdev.2011.08.008 but probably better documented elsewhere
    % see https://www.learningeeg.com/montages-and-technical-components
    % this particular implementation only considers the 10% distances
    montage = [];
    montage.labelold = {'Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'Oz', 'O2'};
    montage.labelnew = {
      'Fp1-Fp2'
      'F7-F3'
      'F3-Fz'
      'Fz-F4'
      'F4-F8'
      'T3-C3'
      'C3-Cz'
      'Cz-C4'
      'C4-T4'
      'T5-P3'
      'P3-Pz'
      'Pz-P4'
      'P4-T6'
      'O1-O2'
      };
    % construct the montage from new channel labels
    montage.tra = zeros(length(montage.labelnew), length(montage.labelold));
    for i=1:length(montage.labelnew)
      lab = split(montage.labelnew{i}, '-');
      montage.tra(i, strcmp(montage.labelold, lab{1})) = +1;
      montage.tra(i, strcmp(montage.labelold, lab{2})) = -1;
    end
    % do a sanity check on the montage
    assert(all(sum(montage.tra, 2)==0));
    assert(all(max(montage.tra, [], 2)==+1));
    assert(all(min(montage.tra, [], 2)==-1));

  otherwise
    error('unsupported refmethod=''%s''', cfg.refmethod);
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble provenance
ft_postamble previous data
ft_postamble history montage

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function returnPathway = colorPathway(pathway, reactionIDs, fluxes, referenceFluxes, ...
    cutOff, defaultColor, upColor, downColor)
% colorPathway
%	Calculates the color for those reactions that should be marked in the
%	map and returns an updated pathway structure.
%
%   pathway         pathway structure of the metabolic network
%   reactionsIDs    cell array with the names of the reactions in the model
%   fluxes          vector with flux values from the simulation
%   referenceFluxes vector with fluxes for the reference simulation
%   cutOff          a reaction is only colored if the absolute value of
%                   at least one of the fluxes is above the cutoff value
%                   (opt, default 0)
%   defaultColor    a color in Matlab format to be used if there are no
%                   changes between the fluxes. This color is also used to
%                   calculate the transition between the colors for up and
%                   down regulated fluxes (opt, default [1 1 1])
%   upColor         a color in Matlab format to be used if the flux is
%                   larger than the reference flux (opt, default [0 1 0])
%   downColor       a color in Matlab format to be used if the flux is
%                   smaller than the reference flux (opt, default [1 0 0])
%
%   returnPathway   updated pathway structure which contains coloring data
%       color       array indicating the coloring of the enzyme
%       signChange  true if the reaction has changed direction
%
%	Usage: returnPathway = colorPathway(pathway, reactionIDs, fluxes, referenceFluxes, ...
%           cutOff, defaultColor, upColor, downColor)

if nargin<8
    downColor=[1 0 0];
end
if nargin<7
    upColor=[0 1 0];
end
if nargin<6
    defaultColor=[1 1 1];
end
if nargin<5
    cutOff=0;
end

%Loop through the components in the pathway. Check if any component has a
%note and then check for the corresponding reaction id
returnPathway=pathway;

for i=1:length(pathway.listOfSpecies)
    if strcmpi(pathway.listOfSpecies(i).type,'PROTEIN')
        if isfield(pathway.listOfSpecies(i),'note')
            if ~isempty(pathway.listOfSpecies(i).note)
                %If there is a note check if there is a corresponding
                %reaction id
                correspondingReaction=find(strcmpi(reactionIDs,pathway.listOfSpecies(i).note));
                if correspondingReaction>0
                    %Check if either flux is above the cutoff value
                    if abs(referenceFluxes(correspondingReaction))>=cutOff ||...
                            abs(fluxes(correspondingReaction))>=cutOff
                        %Calculate the corresponding color
                        [color, signChange]=getColorCodes(referenceFluxes(correspondingReaction)...
                            ,fluxes(correspondingReaction),1 , defaultColor, upColor, downColor);
                        if ~isempty(color)
                            returnPathway.listOfSpecies(i).color=color{1,1};
                            returnPathway.listOfSpecies(i).signChange=signChange{1,1};
                        end
                    end
                end
            end
        end
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
function tagS = RT_DVH_module_tags
%"RT_DVH_module_tags"
%   Return the tags used to represent a DVHs as specified by C.8.8.4 in
%   PS3.3 of 2006 DICOM specification.
%
%   Tags are returned in a struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%JRA 06/06/06
%
%Usage:
%   tagS = RT_DVH_module_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Initialize the tagS structure.
tagS = struct('tag', {}, 'tagdec', {}, 'type', {}, 'children', {});

%Create an empty tagS template for sequence creation.
template = tagS;

%Add tags based on PS3.3 attribute lists.

%Referenced Structure Set Sequence
tagS(end+1) = struct('tag', '300C0060', 'tagdec', 806092896, 'type', '1', 'children', []);
child_1 = template;

%Referenced SOP Class UID
child_1(end+1) = struct('tag', '00081150', 'tagdec', 528720, 'type', '1', 'children', []);

%Referenced SOP Instance UID
child_1(end+1) = struct('tag', '00081155', 'tagdec', 528725, 'type', '1', 'children', []);
tagS(end).children = child_1;

%DVH Normalization Point
tagS(end+1) = struct('tag', '30040040', 'tagdec', 805568576, 'type', '3', 'children', []);

%DVH Normalization Dose Value
tagS(end+1) = struct('tag', '30040042', 'tagdec', 805568578, 'type', '3', 'children', []);

%DVH Sequence
tagS(end+1) = struct('tag', '30040050', 'tagdec', 805568592, 'type', '1', 'children', []);
child_1 = template;

%DVH Referenced ROI Sequence
child_1(end+1) = struct('tag', '30040060', 'tagdec', 805568608, 'type', '1', 'children', []);
child_2 = template;

%Referenced ROI Number
child_2(end+1) = struct('tag', '30060084', 'tagdec', 805699716, 'type', '1', 'children', []);

%DVH ROI Contribution Type
child_2(end+1) = struct('tag', '30040062', 'tagdec', 805568610, 'type', '1', 'children', []);
child_1(end).children = child_2;

%DVH Type
child_1(end+1) = struct('tag', '30040001', 'tagdec', 805568513, 'type', '1', 'children', []);

%Dose Units
child_1(end+1) = struct('tag', '30040002', 'tagdec', 805568514, 'type', '1', 'children', []);

%Dose Type
child_1(end+1) = struct('tag', '30040004', 'tagdec', 805568516, 'type', '1', 'children', []);

%DVH Dose Scaling
child_1(end+1) = struct('tag', '30040052', 'tagdec', 805568594, 'type', '1', 'children', []);

%DVH Volume Units
child_1(end+1) = struct('tag', '30040054', 'tagdec', 805568596, 'type', '1', 'children', []);

%DVH Number of Bins
child_1(end+1) = struct('tag', '30040056', 'tagdec', 805568598, 'type', '1', 'children', []);

%DVH Data
child_1(end+1) = struct('tag', '30040058', 'tagdec', 805568600, 'type', '1', 'children', []);

%DVH Minimum Dose
child_1(end+1) = struct('tag', '30040070', 'tagdec', 805568624, 'type', '1', 'children', []);

%DVH Maximum Dose
child_1(end+1) = struct('tag', '30040072', 'tagdec', 805568626, 'type', '1', 'children', []);

%DVH Mean Dose
child_1(end+1) = struct('tag', '30040074', 'tagdec', 805568628, 'type', '1', 'children', []);

tagS(end).children = child_1;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function A = A_Matrix(phi,th,psi,v,Lt)

%This is the jacobian of the system dynamics. These equations are derived
%in the function Derive_EoM.m

A = [...

 0, 0, -v*cos(phi)*cos(psi)*cos(th),  v*cos(psi)*sin(phi)*sin(th);
 0, 0, -v*cos(phi)*cos(psi)*sin(th), -v*cos(psi)*cos(th)*sin(phi);
 0, 0,                            0,     (v*cos(phi)*cos(psi))/Lt;
 0, 0,                            0,    -(v*cos(phi)*cos(psi))/Lt];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
tic;
doChild = true;
if(ispc)
    tdir = '\\white.stanford.edu\biac2-wandell2\data\reading_longitude\templates\';
    baseDir = '\\white.stanford.edu\biac2-wandell2\data\reading_longitude\';
else
    %tdir = '/biac2/wandell2/data/reading_longitude/templates/';
    tdir = '/silver/scr1/data/templates/';
    baseDir = '/biac2/wandell2/data/reading_longitude/';
end
if(doChild)
    baseDir = fullfile(baseDir,'dti','*0*');
    tdir = fullfile(tdir,'child_new');
    % es's data are OK now.
    %excludeList = {'es041113','tk040817'};
    excludeList = {'tk040817'};
    tBaseName = 'SIRL%02d';
else
    baseDir = fullfile(baseDir,'dti_adults','*0*');
    tdir = fullfile(tdir,'adult_new');
    excludeList = {'ah051003','gf050826','da050311','bw040806'};
    tBaseName = 'SIRL%02dadult';
end

if(isempty(mfilename)) thisfilename = 'dtiBuildAdultTemplate.m';
else  thisfilename = mfilename; end
%[files,subCodes] = findSubjects(baseDir,
%'*_dt6',{'ah051003','gf050826'});
% 2006.03.28 RFD: we now excude da050311 because the top of the
% brain is cut-off in the dti and 'bw040806' as an age outlier.
[files,subCodes] = findSubjects(baseDir, '*_dt6_noMask',excludeList);
N = length(files);

forceBuildAll = false;
useBrainMaskScales = true;
tensorInterpMethod = 1;
tBaseName = sprintf(tBaseName,N);
useBrainMask = true;
if(useBrainMaskScales) suffix{1} = '_bmms';
else suffix{1} = '_ms'; end
suffix{2} = 'warp1'; suffix{3} = 'warp2'; suffix{4} = 'warp3';
mmPerVox = [1 1 1];

if(~exist(tdir,'dir'))
    mkdir(tdir);
end

spm_defaults; global defaults; defaults.analyze.flip = 0;
params = defaults.normalise.estimate;
params.smosrc = 4;

% 1. BOOTSTRAP THE INITIAL TEMPLATE
curName = [tBaseName suffix{1}];
curTemplate = fullfile(tdir,curName);
if(forceBuildAll | ~exist([curTemplate '.img'],'file'))
    % Build the zeroth-iteration template from a simple average
    % of all the ac-pc aligned brains

    % measure brain masks to get the minimumand maximum brain extent in
    % each dimension.
    clear minAcPc maxAcPc;
    for(ii=1:N)
        fprintf(['Measuring ' subCodes{ii} ':\t']);
        d = load(files{ii},'anat');
        % find x,y limits
        tmp = sum(d.anat.brainMask,3);
        x = find(sum(tmp,1)); x = [x(1) x(end)];
        y = find(sum(tmp,2)); y = [y(1) y(end)];
        % find z limits
        tmp = squeeze(sum(d.anat.brainMask,1));
        z = find(sum(tmp,1)); z = [z(1) z(end)];
        curSizeAcPc = mrAnatXformCoords(d.anat.xformToAcPc,[y;x;z]');
        minAcPc(ii,:) = curSizeAcPc(1,:);
        maxAcPc(ii,:) = curSizeAcPc(2,:);
        fprintf('RAS = [% 3.1f % 3.1f % 3.1f] mm, LPI: [% 3.1f % 3.1f % 3.1f] mm.\n',curSizeAcPc(2,:),abs(curSizeAcPc(1,:)));
    end
    meanBrainSizeFromMask = [mean(minAcPc); mean(maxAcPc)];

    % Load all the talscales
    for(ii=[1:N])
        fprintf('Loading scale factors from %s (%d of %d)...\n',subCodes{ii}, ii, N);
        d = load(files{ii},'anat');
        sf(ii) = d.anat.talScale;
    end
    % To get the mean extent, we average the Talairach scale factors
    % and then divide Talairach's actual dimensions (in mm) by these
    % scales.
    tal = mrAnatGetTalairachDists;
    if(useBrainMaskScales)
        meanExtent.sac = abs(meanBrainSizeFromMask(2,3));
        meanExtent.iac = abs(meanBrainSizeFromMask(1,3));
        meanExtent.lac = abs(meanBrainSizeFromMask(1,1));
        meanExtent.rac = abs(meanBrainSizeFromMask(2,1));
        meanExtent.aac = abs(meanBrainSizeFromMask(2,2));
        % We still need the AC-to-PC talScale, since the brain mask doesn't
        % give us that one. But, that one is also immune to the ambiguity
        % that plagues marking the brain-edge.
        meanExtent.acpc = abs(tal.acpc)/mean([sf(:).acpc]);
        meanExtent.ppc = abs(meanBrainSizeFromMask(1,2)) - meanExtent.acpc;
    else
        meanExtent.sac = abs(tal.sac)/mean([sf(:).sac]);
        meanExtent.iac = abs(tal.iac)/mean([sf(:).iac]);
        meanExtent.lac = abs(tal.lac)/mean([sf(:).lac]);
        meanExtent.rac = abs(tal.rac)/mean([sf(:).rac]);
        meanExtent.aac = abs(tal.aac)/mean([sf(:).aac]);
        meanExtent.acpc = abs(tal.acpc)/mean([sf(:).acpc]);
        meanExtent.ppc = abs(tal.ppc)/mean([sf(:).ppc]);
    end
    for(ii=[1:N])
        if(useBrainMaskScales)
            extent(ii).sac = abs(maxAcPc(ii,3));
            extent(ii).iac = abs(minAcPc(ii,3));
            extent(ii).lac = abs(minAcPc(ii,1));
            extent(ii).rac = abs(maxAcPc(ii,1));
            extent(ii).aac = abs(maxAcPc(ii,2));
            extent(ii).acpc = abs(tal.acpc)/sf(ii).acpc;
            extent(ii).ppc = abs(minAcPc(ii,2)) - meanExtent.acpc;
        else
            extent(ii).sac = abs(tal.sac)/sf(ii).sac;
            extent(ii).iac = abs(tal.iac)/sf(ii).iac;
            extent(ii).lac = abs(tal.lac)/sf(ii).lac;
            extent(ii).rac = abs(tal.rac)/sf(ii).rac;
            extent(ii).aac = abs(tal.aac)/sf(ii).aac;
            extent(ii).acpc = abs(tal.acpc)/sf(ii).acpc;
            extent(ii).ppc = abs(tal.ppc)/sf(ii).ppc;
        end
        newScale(ii).sac = meanExtent.sac/extent(ii).sac;
        newScale(ii).iac = meanExtent.iac/extent(ii).iac;
        newScale(ii).lac = meanExtent.lac/extent(ii).lac;
        newScale(ii).rac = meanExtent.rac/extent(ii).rac;
        newScale(ii).aac = meanExtent.aac/extent(ii).aac;
        newScale(ii).acpc = meanExtent.acpc/extent(ii).acpc;
        newScale(ii).ppc = meanExtent.ppc/extent(ii).ppc;
        newScale(ii).pcReference = -meanExtent.acpc;
    end

    % Reslice using the talscale measurements.
    d = load(files{1},'anat');
    ac = round(mrAnatXformCoords(inv(d.anat.xformToAcPc),[0 0 0]));
    im = double(d.anat.img);
    im = im./max(im(:));
    im(~d.anat.brainMask) = 0;
    sz = size(im);
    % t1's from dt6 files shouldn't need histogram clipping. But we will
    % scale every brain to 0-1.
    %imSkull = mrAnatHistogramClip(double(imSkull), 0.5, 0.99);
    %proportionNonBrain = sum(~d.anat.brainMask(:))./prod(size(d.anat.brainMask));
    %im = mrAnatHistogramClip(double(im), proportionNonBrain, 0.99);
    bb = [-d.anat.mmPerVox.*(ac-1); d.anat.mmPerVox.*(sz-ac)];
    ts = newScale(1); ts.talScaleDir = 'tal2acpc';
    % warp T1 to tal space
    ts.outMat = inv(d.anat.xformToAcPc);
    [im,newXform] = mrAnatResliceSpm(im,ts,bb,mmPerVox,[7 7 7 0 0 0],0);
    im(isnan(im)) = 0;
    im(im<0) = 0; im(im>1) = 1;
    allIm = zeros([sz N]);
    allIm(:,:,:,1) = im;
    for(ii=[2:N])
        fprintf('Processing %s (%d of %d)...\n',subCodes{ii}, ii, N);
        d = load(files{ii},'anat');
        % All coords should be the same.
        ac = round(mrAnatXformCoords(inv(d.anat.xformToAcPc),[0 0 0]));
        im = double(d.anat.img);
        im = im./max(im(:));
        im(~d.anat.brainMask) = 0;
        ts = newScale(ii); ts.talScaleDir = 'tal2acpc';
        ts.outMat = inv(d.anat.xformToAcPc);
        [im,newXform] = mrAnatResliceSpm(im,ts,bb,mmPerVox,[7 7 7 0 0 0],0);
        im(isnan(im)) = 0;
        im(im<0) = 0; im(im>1) = 1;
        allIm(:,:,:,ii) = im;
    end
    % Save all the resliced brains in a temp file.
    save(fullfile(tempdir, [curName '_all.mat']), 'allIm');
    for(ii=1:N)
        allIm(:,:,:,ii) = mrAnatHistogramClipOptimal(allIm(:,:,:,ii));
    end
    meanIm = mean(allIm,4);
    figure; image(makeMontage(uint8(meanIm*255))); colormap(gray(256)); axis equal tight off;

    ax = flipdim(permute(squeeze(allIm(:,:,ac(3),:)),[2,1,3]),1);
    cr = flipdim(permute(squeeze(allIm(:,ac(2),:,:)),[2,1,3]),1);
    sg = flipdim(permute(squeeze(allIm(ac(1),:,:,:)),[2,1,3]),1);
    ax(:,:,N+1) = mean(ax(:,:,[1:N]), 3);
    cr(:,:,N+1) = mean(cr(:,:,[1:N]), 3);
    sg(:,:,N+1) = mean(sg(:,:,[1:N]), 3);
    ax = uint8(ax.*255+0.5);
    cr = uint8(cr.*255+0.5);
    sg = uint8(sg.*255+0.5);

    figure;image(makeMontage(ax));colormap(gray(256));axis equal tight off;
    figure;image(makeMontage(cr));colormap(gray(256));axis equal tight off;
    figure;image(makeMontage(sg));colormap(gray(256));axis equal tight off;

    imwrite(makeMontage(ax), gray(256), fullfile(tdir, [tBaseName '_ax.png']));
    imwrite(makeMontage(cr), gray(256), fullfile(tdir, [tBaseName '_cr.png']));
    imwrite(makeMontage(sg), gray(256), fullfile(tdir, [tBaseName '_sg.png']));

    % Now build an atlas from the mean image
    V.dat = int16(meanIm.*(2^15-1)+0.5);
    notes = [curName ': average of ' num2str(N) ' mean-Tal-scaled brains. Created at ' datestr(now,31)];
    hdr = saveAnalyze(V.dat, curTemplate, mmPerVox, notes, ac);
    save(fullfile(tdir,[curName '_details.mat']), 'meanIm','notes','extent','meanExtent','ac','mmPerVox');
else
   disp([curTemplate ' exists- skipping. (Set forceBuildAll to rebuild.)']);
end

% 2. BUILD THE INTERMEDIATE TEMPLATE BY WARPING TO THE INITIAL TEMPLATE
prevTemplate = curTemplate;
curName = [tBaseName suffix{2}];
curTemplate = fullfile(tdir,curName);
if(forceBuildAll | ~exist([curTemplate '.img'],'file'))
    % Now use the mean-scaled template to generate a refined template using
    % non-linear warping.

    [im,mm,hdr] = loadAnalyze(prevTemplate);
    ac = round(mrAnatXformCoords(inv(hdr.mat),[0 0 0]));

    d = load(files{1},'anat');
    im = double(d.anat.img);
    im = im./max(im(:));
    sz = size(im);
    bb = [-d.anat.mmPerVox.*(ac-1); d.anat.mmPerVox.*(sz-ac)];
    % warp T1 to template space
    xform = d.anat.xformToAcPc;
    if(useBrainMask) im(~d.anat.brainMask) = 0; end
    sn = mrAnatComputeSpmSpatialNorm(im, xform, [prevTemplate '.img'], params);
    sn.outMat = inv(sn.VF.mat);
    [im,newXform] = mrAnatResliceSpm(im,sn,bb,mmPerVox,[7 7 7 0 0 0],0);
    im(isnan(im)) = 0;
    im(im<0) = 0; im(im>1) = 1;
    allIm = zeros([sz N]);
    allIm(:,:,:,1) = im;
    for(ii=[2:N])
        fprintf('Processing %s (%d of %d)...\n',subCodes{ii}, ii, N);
        d = load(files{ii},'anat');
        [p,f,e] = fileparts(files{ii});
        im = double(d.anat.img);
        im = im./max(im(:));
        xform = d.anat.xformToAcPc;
        if(useBrainMask) im(~d.anat.brainMask) = 0; end
        sn = mrAnatComputeSpmSpatialNorm(im, xform, [prevTemplate '.img'], params);
        sn.outMat = inv(sn.VF.mat);
        [im,newXform] = mrAnatResliceSpm(im,sn,bb,mmPerVox,[7 7 7 0 0 0],0);
        im(isnan(im)) = 0;
        im(im<0) = 0; im(im>1) = 1;
        allIm(:,:,:,ii) = im;
    end
    % Save all the resliced brains in a temp file.
    save(fullfile(tempdir, [curName '_all.mat']), 'allIm');
    for(ii=1:N)
        allIm(:,:,:,ii) = mrAnatHistogramClipOptimal(allIm(:,:,:,ii));
    end

    meanIm = mean(allIm,4);
    figure; image(makeMontage(uint8(meanIm*255)));
    colormap(gray(256)); axis equal tight off;

    ax = flipdim(permute(squeeze(allIm(:,:,ac(3),:)),[2,1,3]),1);
    cr = flipdim(permute(squeeze(allIm(:,ac(2),:,:)),[2,1,3]),1);
    sg = flipdim(permute(squeeze(allIm(ac(1),:,:,:)),[2,1,3]),1);
    ax(:,:,N+1) = mean(ax(:,:,[1:N]), 3);
    cr(:,:,N+1) = mean(cr(:,:,[1:N]), 3);
    sg(:,:,N+1) = mean(sg(:,:,[1:N]), 3);
    ax = uint8(ax.*255+0.5);
    cr = uint8(cr.*255+0.5);
    sg = uint8(sg.*255+0.5);
    figure;image(makeMontage(ax));colormap(gray(256));axis equal tight off;
    figure;image(makeMontage(cr));colormap(gray(256));axis equal tight off;
    figure;image(makeMontage(sg));colormap(gray(256));axis equal tight off;
    imwrite(makeMontage(ax), gray(256), fullfile(tdir, [tBaseName '_ax.png']));
    imwrite(makeMontage(cr), gray(256), fullfile(tdir, [tBaseName '_cr.png']));
    imwrite(makeMontage(sg), gray(256), fullfile(tdir, [tBaseName '_sg.png']));

    % Now build an atlas from the mean image
    V.dat = int16(meanIm.*(2^15-1)+0.5);
    [p,f,e]=fileparts(prevTemplate); tname = [f e];
    notes = [curName ': average of ' num2str(N) ' brains warped to ' tname '. Created at ' datestr(now,31)];
    hdr = saveAnalyze(V.dat, curTemplate, mmPerVox, notes, ac);
else
    disp([curTemplate ' exists- skipping. (Set forceBuildAll to rebuild.)']);
end

%
% 3. BUILD THE FINAL TEMPLATE BY WARPING TO THE INTERMEDIATE TEMPLATE
%
prevTemplate = curTemplate;
curName = [tBaseName];
curTemplate = fullfile(tdir,curName);
if(forceBuildAll | ~exist([curTemplate '.img'],'file'))
    outDir = fullfile(tdir, [tBaseName suffix{3}]);
    if(~exist(outDir)) mkdir(outDir); end
    for(ii=[1:N])
        fprintf('Processing %s (%d of %d)...\n',subCodes{ii}, ii, N);
        d = load(files{ii});
        im = double(d.anat.img);
        im = im./max(im(:));
        xform = d.anat.xformToAcPc;
        if(useBrainMask) im(~d.anat.brainMask) = 0; end
        sn = mrAnatComputeSpmSpatialNorm(im, xform, [prevTemplate '.img'], params);
        t1Rng = [min(d.anat.img(:)) max(d.anat.img(:))];
        b0Rng = [min(d.b0(:)) max(d.b0(:))];
        d = dtiSpmDeformer(d, sn, tensorInterpMethod, mmPerVox);
        d.anat.img(d.anat.img<t1Rng(1)) = t1Rng(1);
        d.anat.img(d.anat.img>t1Rng(2)) = t1Rng(2);
        d.anat.img = int16(d.anat.img+0.5);
        d.b0(d.b0<b0Rng(1)) = b0Rng(1);
        d.b0(d.b0>b0Rng(2)) = b0Rng(2);
        d.b0 = int16(d.b0+0.5);
        newFileName = fullfile(outDir, [subCodes{ii} '_sn']);
        disp(['Saving to ' newFileName '...']);
        dtiSaveStruct(d, newFileName);
    end

    % Now load all the warped dt6 files and build a template from them.
    snFiles = findSubjects(outDir, '*_sn*',{});
    disp(['Loading ' snFiles{1} '...']);
    dt = load(snFiles{1});
    avg = dt;
    avg.b0 = mrAnatHistogramClip(double(avg.b0),0.5,0.99);
    avg.anat.img = mrAnatHistogramClip(double(avg.anat.img),0.5,0.99);
    if(isfield(avg.anat,'brainMask'))
        avg.anat.brainMask = double(avg.anat.brainMask);
    end
    Nb0 = double(avg.b0>0);
    Nt1 = double(avg.anat.img>0);
    for(ii=2:length(snFiles))
        disp(['Loading ' snFiles{ii} '...']);
        dt = load(snFiles{ii});
        dt.dt6(isnan(dt.dt6)) = 0;
        b0 = mrAnatHistogramClip(double(dt.b0),0.4,0.99);
        b0(isnan(b0)) = 0;
        avg.b0 = avg.b0+b0;
        avg.dt6 = avg.dt6+dt.dt6;
        t1 = mrAnatHistogramClip(double(dt.anat.img),0.4,0.99);
        t1(isnan(t1)) = 0;
        avg.anat.img = avg.anat.img+t1;
        if(isfield(avg.anat,'brainMask') && isfield(dt.anat,'brainMask'))
            avg.anat.brainMask = avg.anat.brainMask+double(dt.anat.brainMask);
        end
        Nb0 = Nb0+double(b0>0);
        Nt1 = Nt1+double(t1>0);
    end
    avg.b0(Nb0>0) = avg.b0(Nb0>0)./Nb0(Nb0>0);
    for(ii=1:6)
        tmp = avg.dt6(:,:,:,ii);
        tmp(Nb0>0) = tmp(Nb0>0)./Nb0(Nb0>0);
        avg.dt6(:,:,:,ii) = tmp;
    end
    avg.anat.img(Nt1>0) = avg.anat.img(Nt1>0)./Nt1(Nt1>0);
    if(isfield(avg.anat,'brainMask'))
        avg.anat.brainMask = avg.anat.brainMask./N;
    end
    % Apply a new dt mask based on the average B0
    mask = dtiCleanImageMask(avg.b0>0.3, 5);
    avg.dt6(repmat(~mask, [1,1,1,6])) = 0;
    avg.notes.dataFiles = snFiles;
    avg.notes.normalizationTemplate = prevTemplate;
    avg.notes.normalizationParams = params;
    avg.notes.createdOn = datestr(now,31);
    fid=fopen(which(mfilename),'rt');thisScript=fread(fid);fclose(fid);
    avg.notes.buildScript = char(thisScript');
    newFileName = fullfile(outDir, ['average_dt6']);
    disp(['Saving to ' newFileName '...']);
    dtiSaveStruct(avg, newFileName);

    % Save mean T1 in analyze format
    V.dat = int16(avg.anat.img.*(2^15-1)+0.5);
    [p,f,e]=fileparts(prevTemplate); tname = [f e];
    notes = ['average of ' num2str(N) ' brains warped to ' tname '. Created at ' datestr(now,31)];
    origin = round(mrAnatXformCoords(inv(avg.anat.xformToAcPc),[0 0 0]));
    saveAnalyze(V.dat, fullfile(tdir,curName), mmPerVox, notes, origin);

    % Save mean T1 (skull-stripped) in analyze format
    % V.dat(avg.anat.brainMask<0.5) = 0;
    % saveAnalyze(V.dat, fullfile(tdir,'SIRL54ms_warp2_brain_thresh'), mmPerVox, notes, origin);
    % An alternative method of skull-stipping. This method essentially
    % treats the brainMask as a probability map rather than a hard
    % threshold. The results are very similar to the hard-thresh
    % version, but it leaves a softer edge in regions of uncertainty
    % (eg. orbital frontal and cerebellium).
    V.dat = int16(avg.anat.img.*avg.anat.brainMask.*(2^15-1)+0.5);
    saveAnalyze(V.dat, fullfile(tdir,[curName '_brain']), mmPerVox, notes, origin);

    % Save mean b=0 in analyze format
    V.dat = int16(avg.b0.*(2^15-1)+0.5);
    origin = round(mrAnatXformCoords(inv(avg.xformToAcPc),[0 0 0]));
    saveAnalyze(V.dat, fullfile(tdir,[curName '_EPI']),  avg.mmPerVox, notes, origin);
else
    disp([curTemplate ' exists- skipping. (Set forceBuildAll to rebuild.)']);
end

%
% 4. BUILD AN AVERAGE DATASET BY WARPING EACH TO THE WARP2 TEMPLATE
%
prevTemplate = fullfile(tdir,[curName '_brain']);
curName = [tBaseName suffix{4}];
outDir = fullfile(tdir, curName);
if(~exist(outDir)) mkdir(outDir); end
for(ii=[1:N])
  fprintf('Processing %s (%d of %d)...\n',subCodes{ii}, ii, N);
  d = load(files{ii});
  im = double(d.anat.img);
  im = im./max(im(:));
  xform = d.anat.xformToAcPc;
  if(useBrainMask) im(~d.anat.brainMask) = 0; end
  sn = mrAnatComputeSpmSpatialNorm(im, xform, [prevTemplate '.img'], params);
  t1Rng = [min(d.anat.img(:)) max(d.anat.img(:))];
  b0Rng = [min(d.b0(:)) max(d.b0(:))];
  d = dtiSpmDeformer(d, sn, tensorInterpMethod, mmPerVox);
  d.anat.img(d.anat.img<t1Rng(1)) = t1Rng(1);
  d.anat.img(d.anat.img>t1Rng(2)) = t1Rng(2);
  d.anat.img = int16(d.anat.img+0.5);
  d.b0(d.b0<b0Rng(1)) = b0Rng(1);
  d.b0(d.b0>b0Rng(2)) = b0Rng(2);
  d.b0 = int16(d.b0+0.5);

  newFileName = fullfile(outDir, [subCodes{ii} '_sn']);
  disp(['Saving to ' newFileName '...']);
  dtiSaveStruct(d, newFileName);
end
% Now load all the warped dt6 files and build an average brain
snFiles = findSubjects(outDir, '*_sn*',{});
N = length(snFiles);
disp(['Loading ' snFiles{1} '...']);
dt = load(snFiles{1});
avg = dt;
avg.b0 = mrAnatHistogramClip(double(avg.b0),0.5,0.99);
avg.anat.img = mrAnatHistogramClip(double(avg.anat.img),0.5,0.99);
if(isfield(avg.anat,'brainMask'))
  avg.anat.brainMask = double(avg.anat.brainMask);
end
Nb0 = double(avg.b0>0);
Nt1 = double(avg.anat.img>0);
for(ii=2:N)
  disp(['Loading ' snFiles{ii} '...']);
  dt = load(snFiles{ii});
  dt.dt6(isnan(dt.dt6)) = 0;
  b0 = mrAnatHistogramClip(double(dt.b0),0.4,0.99);
  b0(isnan(b0)) = 0;
  avg.b0 = avg.b0+b0;
  avg.dt6 = avg.dt6+dt.dt6;
  t1 = mrAnatHistogramClip(double(dt.anat.img),0.4,0.99);
  t1(isnan(t1)) = 0;
  avg.anat.img = avg.anat.img+t1;
  if(isfield(avg.anat,'brainMask') && isfield(dt.anat,'brainMask'))
    avg.anat.brainMask = avg.anat.brainMask+double(dt.anat.brainMask);
  end
  Nb0 = Nb0+double(b0>0);
  Nt1 = Nt1+double(t1>0);
end
avg.b0(Nb0>0) = avg.b0(Nb0>0)./Nb0(Nb0>0);
for(ii=1:6)
  tmp = avg.dt6(:,:,:,ii);
  tmp(Nb0>0) = tmp(Nb0>0)./Nb0(Nb0>0);
  avg.dt6(:,:,:,ii) = tmp;
end
avg.anat.img(Nt1>0) = avg.anat.img(Nt1>0)./Nt1(Nt1>0);
if(isfield(avg.anat,'brainMask'))
  avg.anat.brainMask = avg.anat.brainMask./N;
end
% Apply a new dt mask based on the average B0
mask = dtiCleanImageMask(avg.b0>0.3, 5);
avg.dt6(repmat(~mask, [1,1,1,6])) = 0;
avg.notes.dataFiles = snFiles;
avg.notes.normalizationTemplate = prevTemplate;
avg.notes.normalizationParams = params;
avg.notes.createdOn = datestr(now,31);
fid=fopen(which(mfilename),'rt');thisScript=fread(fid);fclose(fid);
avg.notes.buildScript = char(thisScript');
newFileName = fullfile(outDir, ['average_dt6']);
disp(['Saving to ' newFileName '...']);
dtiSaveStruct(avg, newFileName);
disp(['Finished (' num2str(toc) ' secs)']);

error('stop here');

%% Get demographics for subjects
%
%
baseDir = '/biac2/wandell2/data/reading_longitude/';
tdir = fullfile(baseDir,'templates','adult');
[files,subCodes] = findSubjects(fullfile(tdir,'SIRL20adultwarp3_averageDataset'), '*_sn', {});
fid = fopen(fullfile(tdir,'subjectInfo.csv'),'w');
fprintf(fid,'Name,Sex,Age,ExactAge,birthDate,scanDate\n');
for(ii=1:length(subCodes))
    if(strcmp(subCodes{ii}(1:2),'ab')) subCodes{ii} = ['aab' subCodes{ii}(3:end)]; end
    if(strcmp(subCodes{ii}(1:2),'rd')) subCodes{ii} = ['rfd' subCodes{ii}(3:end)]; end
    dtDir = fullfile(baseDir,'dti_adults',subCodes{ii});
    if(~exist(dtDir,'dir'))
        dtDir = fullfile(baseDir,'dti_adults','additional_runs',subCodes{ii});
    end
    dataDir = fullfile(dtDir,'dti');
    if(~exist(dataDir,'dir'))
        dataDir = fullfile(dtDir,'dti_6dir');
    end
    f = fullfile(dataDir,'B0_001.dcm');
    if(~exist(f,'file'))
        f = [f '.gz'];
        gunzip(f,'/tmp/');
        f = '/tmp/B0_001.dcm';
    end
    di = dicominfo(f);
    ad = datenum(di.AcquisitionDate,'yyyymmdd');
	bd = datenum(di.PatientBirthDate,'yyyymmdd');
    exactAge(ii) = (ad-bd)/364.25;
    sex(ii) = di.PatientSex=='M';
    fprintf(fid,'%s,%s,%s,%s,%0.2f,%s,%s\n',...
        subCodes{ii}, di.PatientName.FamilyName, di.PatientSex, di.PatientAge, exactAge(ii), di.PatientBirthDate, di.AcquisitionDate);
end
fclose(fid);
fprintf('Age: %0.0f - %0.0f years (mean %0.1f), %d females out of %d total.\n',...
    min(exactAge),max(exactAge),mean(exactAge),sum(~sex),length(sex));

%
% 5. BUILD AN MNI-EPI SET FOR COMPARISON
%
tic;
if(ispc)
    tdir = '\\white.stanford.edu\biac2-wandell2\data\reading_longitude\templates\adult';
    baseDir = '\\white.stanford.edu\biac2-wandell2\data\reading_longitude\dti_adults\*0*';
else
    tdir = '/biac2/wandell2/data/reading_longitude/templates/adult';
    baseDir = '/biac2/wandell2/data/reading_longitude/dti_adults/*0*';
end
tensorInterpMethod = 1;
tBaseName = 'SIRL20_adult';
mmPerVox = [1 1 1];
spm_defaults; global defaults; defaults.analyze.flip = 0;
params = defaults.normalise.estimate;
%params.smosrc = 4;
[files,subCodes] = findSubjects(baseDir, '*_dt6',{'ah051003','gf050826'});
N = length(files);

spmDir = fileparts(which('spm_normalise'));
prevTemplate = fullfile(spmDir,'templates','EPI.mnc');
curName = [tBaseName '_MNI_EPI'];
outDir = fullfile(tdir, curName);
if(~exist(outDir)) mkdir(outDir); end
for(ii=[1:N])
  fprintf('Processing %s (%d of %d)...\n',subCodes{ii}, ii, N);
  d = load(files{ii});
  im = double(d.b0);
  im = mrAnatHistogramClip(im);
  xform = d.anat.xformToAcPc;
  if(useBrainMask) im(~d.anat.brainMask) = 0; end
  sn = mrAnatComputeSpmSpatialNorm(im, xform, [prevTemplate '.img'], params);
  t1Rng = [min(d.anat.img(:)) max(d.anat.img(:))];
  b0Rng = [min(d.b0(:)) max(d.b0(:))];
  d = dtiSpmDeformer(d, sn, tensorInterpMethod, mmPerVox);
  d.anat.img(d.anat.img<t1Rng(1)) = t1Rng(1);
  d.anat.img(d.anat.img>t1Rng(2)) = t1Rng(2);
  d.anat.img = int16(d.anat.img+0.5);
  d.b0(d.b0<b0Rng(1)) = b0Rng(1);
  d.b0(d.b0>b0Rng(2)) = b0Rng(2);
  d.b0 = int16(d.b0+0.5);

  newFileName = fullfile(outDir, [subCodes{ii} '_sn']);
  disp(['Saving to ' newFileName '...']);
  dtiSaveStruct(d, newFileName);
end
% Now load all the warped dt6 files and build an average brain
files = findSubjects(outDir, '*_sn*',{});
N = length(files);
disp(['Loading ' files{1} '...']);
dt = load(files{1});
avg = dt;
avg.b0 = mrAnatHistogramClip(double(avg.b0),0.5,0.99);
avg.anat.img = mrAnatHistogramClip(double(avg.anat.img),0.5,0.99);
if(isfield(avg.anat,'brainMask'))
  avg.anat.brainMask = double(avg.anat.brainMask);
end
avg.dtBrainMask = double(avg.dtBrainMask);
Nb0 = double(avg.b0>0);
Nt1 = double(avg.anat.img>0);
for(ii=2:N)
  disp(['Loading ' files{ii} '...']);
  dt = load(files{ii});
  dt.dt6(isnan(dt.dt6)) = 0;
  b0 = mrAnatHistogramClip(double(dt.b0),0.4,0.99);
  b0(isnan(b0)) = 0;
  avg.b0 = avg.b0+b0;
  avg.dt6 = avg.dt6+dt.dt6;
  avg.dtBrainMask = avg.dtBrainMask+dt.dtBrainMask;
  t1 = mrAnatHistogramClip(double(dt.anat.img),0.4,0.99);
  t1(isnan(t1)) = 0;
  avg.anat.img = avg.anat.img+t1;
  if(isfield(avg.anat,'brainMask') && isfield(dt.anat,'brainMask'))
    avg.anat.brainMask = avg.anat.brainMask+double(dt.anat.brainMask);
  end
  Nb0 = Nb0+double(b0>0);
  Nt1 = Nt1+double(t1>0);
end
avg.b0(Nb0>0) = avg.b0(Nb0>0)./Nb0(Nb0>0);
for(ii=1:6)
  tmp = avg.dt6(:,:,:,ii);
  tmp(Nb0>0) = tmp(Nb0>0)./Nb0(Nb0>0);
  avg.dt6(:,:,:,ii) = tmp;
end
avg.anat.img(Nt1>0) = avg.anat.img(Nt1>0)./Nt1(Nt1>0);
if(isfield(avg.anat,'brainMask'))
  avg.anat.brainMask = avg.anat.brainMask./N;
end
avg.dtBrainMask = avg.dtBrainMask./N;
% Apply a new dt mask based on the average B0
mask = dtiCleanImageMask(avg.b0>0.3, 5);
avg.dt6(repmat(~mask, [1,1,1,6])) = 0;
avg.notes.dataFiles = files;
avg.notes.normalizationTemplate = template;
avg.notes.normalizationParams = params;
avg.notes.createdOn = datestr(now,31);
fid=fopen(which(thisfilename),'rt');thisScript=fread(fid);fclose(fid);
avg.notes.buildScript = char(thisScript');
newFileName = fullfile(outDir, ['average_dt6']);
disp(['Saving to ' newFileName '...']);
dtiSaveStruct(avg, newFileName);
disp(['Finished (' num2str(toc) ' secs)']);


text:
% Batch run autoclusting, Sweep master threshold with 2-fold CV
% ie all clustering thresholds are the same
% Need to run Batch_k20_CV12_AK first

%% Initialize Parameters

range_fish = 8:15; % select which fish to analyze (can be multiple)
M_stim = 2; % select which stimuli to use when clustering
masterThresh_sweep = 0.5:.1:0.9; % values of thresh to sweep

isFullData = 1; % If using all cells

masterDef = 0.7; % Default values of clusParams
mergeDef = masterDef;
capDef = masterDef;
reg1Def = masterDef;
reg2Def = masterDef;
minSizeDef = 10;
k1Def = 20;
k2Def = 500;

clusParams = struct('merge',mergeDef,'cap',capDef,'reg1',reg1Def,...
    'reg2',reg2Def,'minSize',minSizeDef,'k1',k1Def,'k2',k2Def);

masterThresh_nClus = zeros(length(masterThresh_sweep),length(range_fish));
masterThresh_CVscore = zeros(length(masterThresh_sweep),length(range_fish));
masterThresh_CVscore_raw = cell(length(masterThresh_sweep),length(range_fish));
masterThresh_nCells = zeros(length(masterThresh_sweep),length(range_fish));
masterThresh_compTime = zeros(length(masterThresh_sweep),length(range_fish));

masterThresh_data = struct('clusParams',clusParams,'nClus',masterThresh_nClus,'CVscore',masterThresh_CVscore,...
    'nCells',masterThresh_nCells,'compTime',masterThresh_compTime,'clusA',[],'clusB',[]);

data_masterdir = GetCurrentDataDir();

%%
disp(['Batch Script: Sweep masterThresh, iFish']);
scriptStart = tic;
clear leg;
for k_fish = 1:length(range_fish),

    %% Load fish data
    i_fish = range_fish(k_fish);
    leg{k_fish} = ['Fish ' num2str(range_fish(k_fish))];
    LoadFullFish(hfig,i_fish,isFullData); % Load the data
    absIX = getappdata(hfig,'absIX');

    %% Partitions for CV
    timelists = getappdata(hfig,'timelists');
    timelists_names = getappdata(hfig,'timelists_names');
    periods = getappdata(hfig,'periods');

    timelistsCV = cell(length(M_stim),2);

    % Divides each stimulous in half
    for k_stim = 1:length(M_stim),
        i_stim = M_stim(k_stim);
        TL = timelists{i_stim};
        period = periods(i_stim);
        nrep = size(TL,2)/periods(i_stim); % integer
        n = floor(nrep/2);
        timelistsCV{k_stim,1} = TL(1):TL(n*period);
        timelistsCV{k_stim,2} = TL(1+n*period):TL(2*n*period);
    end

    %% Sweep through k2 and do clustering
    for k_masterThresh = 1:length(masterThresh_sweep),

        sweepStart = tic;
        clusParams.merge = masterThresh_sweep(k_masterThresh);
        clusParams.cap = masterThresh_sweep(k_masterThresh);
        clusParams.reg1 = masterThresh_sweep(k_masterThresh);
        clusParams.reg2 = masterThresh_sweep(k_masterThresh);
        Score = zeros(1,2);%(length(M_stim),2);
        %%
        nClus = zeros(1,2);
        nCells = zeros(1,2);
        CIX = cell(1,2);
        GIX = cell(1,2);
        for half = 1:2, % CV halves, load them from saved
            %%
            i_ClusGroup = 2;
            i_Cluster = half+2;
            % this matches the assignments
            % in Batch_k20_CV12_AK
            [cIX,gIX] = LoadCluster_Direct(i_fish,i_ClusGroup,i_Cluster,absIX);

            tIX = timelistsCV{k_stim,half};
            M_0 = GetTimeIndexedData_Default_Direct(hfig,[],tIX,'isAllCells');

            isWkmeans = 0;
            [cIX,gIX] = AutoClustering(cIX,gIX,M_0,isWkmeans,clusParams);
%%
            nClus(half) = length(unique(gIX));
            nCells(half) = length(unique(cIX));
            CIX{half} = cIX;
            GIX{half} = gIX;
        end

        % Perform 2x Cross Validation
        Score(1) = HungarianCV(nClus(1),nClus(2),CIX{1},CIX{2},GIX{1},GIX{2});% true,timelists_names{i_stim});
        Score(2) = HungarianCV(nClus(2),nClus(1),CIX{2},CIX{1},GIX{2},GIX{1});% true,timelists_names{i_stim});
        masterThresh_CVscore_raw{k_masterThresh,k_fish} = Score;

        % Save clusters and scores

        masterThresh_CVscore(k_masterThresh,k_fish) = mean(Score);
        masterThresh_nClus(k_masterThresh,k_fish) = mean([nClus(1),nClus(2)]);
        masterThresh_nCells(k_masterThresh,k_fish) = mean([nCells(1),nCells(2)]);
        masterThresh_compTime(k_masterThresh,k_fish) = toc(sweepStart);

        masterThresh_data(k_masterThresh,k_fish).CVscore = mean(Score);
        masterThresh_data(k_masterThresh,k_fish).nClus = mean([nClus(1),nClus(2)]);
        masterThresh_data(k_masterThresh,k_fish).nCells = mean([nCells(1),nCells(2)]);
        masterThresh_data(k_masterThresh,k_fish).clusParams = clusParams;
        masterThresh_data(k_masterThresh,k_fish).compTime = toc(sweepStart);

        masterThresh_data(k_masterThresh,k_fish).clusA = struct('cIX',CIX{1},'gIX',GIX{1});
        masterThresh_data(k_masterThresh,k_fish).clusB = struct('cIX',CIX{2},'gIX',GIX{2});
        end
end
scriptTime = toc(scriptStart);
disp(['Batch Script Finished, Took ' num2str(scriptTime) ' sec']);

%%
figure;
subplot(2,2,1)
plot(masterThresh_sweep,masterThresh_nClus,'-o')
% legend('Fish8','Fish9');
%xlabel('total k for kmeans')
ylabel('# of auto-clusters')

subplot(2,2,2)
plot(masterThresh_sweep,masterThresh_CVscore,'-o')
ylim([0,1])
ylabel('CV (overlapping cell %)')

subplot(2,2,3)
plot(masterThresh_sweep,masterThresh_nCells,'-o');
%xlabel('total k for kameans')
ylabel('total # of cells included in clusters')
xlabel('master threshold')

subplot(2,2,4)
plot(masterThresh_sweep,masterThresh_nCells.*masterThresh_CVscore,'-o');
%xlabel('total k for kameans')
ylabel('total # of cells passing CV')
xlabel('master threshold')
legend(leg,'location','best');

author: xiuyechen

repo: FishExplorer


text:
function [cfg] = ft_multiplotCC(cfg, data)

% FT_MULTIPLOTCC visualises the coherence between channels by using
% multiple topoplots. The topoplot at a given channel location shows the
% coherence of that channel with all other channels.
%
% Use as
%   ft_multiplotCC(cfg, data)
%
% See also FT_PREPARE_LAYOUT, FT_TOPOPLOTCC, FT_CONNECTIVITYPLOT

% Undocumented local options:
% cfg.layout  = layout filename or a structure produced by prepare_layout
% cfg.xlim
% cfg.parameter
% This function requires input from FT_FREQSTATISTICS_SHIFTPREDICT
% This function should be rewritten, using the clean topoplot implementation

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
data = ft_checkdata(data);

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'zparam', 'parameter'});
cfg = ft_checkconfig(cfg, 'deprecated',  {'xparam'});

% set the defaults
cfg.xlim      = ft_getopt(cfg, 'xlim', 'all');
cfg.parameter = ft_getopt(cfg, 'parameter', 'avg.icohspctrm');
cfg.renderer  = ft_getopt(cfg, 'renderer'); % let MATLAB decide on the default

if strcmp(cfg.parameter, 'avg.icohspctrm') && ~issubfield(data, 'avg.icohspctrm')
  data.avg.icohspctrm = abs(imag(data.avg.cohspctrm));
end

if strcmp(data.dimord, 'chan_chan_freq')
  % reshape input-data, such that ft_topoplotTFR will take it
  cnt = 1;
  siz = size(data.prob);
  data.labelcmb = cell(siz(1)*siz(2),2);
  data.prob = reshape(data.prob, [siz(1)*siz(2) siz(3)]);
  data.stat = reshape(data.stat, [siz(1)*siz(2) siz(3)]);
  for j = 1:length(data.label)
    for k = 1:length(data.label)
      data.labelcmb(cnt,:) = [data.label(k) data.label(j)];
      cnt = cnt + 1;
    end
  end
  tmpdata = data;
else
  dat   = getsubfield(data, cfg.parameter);
  scale = [0 max(dat(:))-0.2];
end

if isfield(cfg, 'xparam')
  xparam = getsubfield(data, cfg.xparam);
  if ~strcmp(cfg.xlim, 'all')
    fbin = [nearest(xparam, cfg.xlim(1)) nearest(xparam, cfg.xlim(2))];
  else
    fbin = [xparam(1) xparam(end)];
  end
end

% read or create the layout that will be used for plotting
tmpcfg = keepfields(cfg, {'layout', 'channel', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
layout = ft_prepare_layout(tmpcfg);
ft_plot_layout(layout, 'box', false, 'label', 'no', 'point', 'no');

X = layout.pos(:,1);
Y = layout.pos(:,2);
Width = layout.width;
Height = layout.height;
Lbl = layout.label;
chNum = numel(layout.label);

xScaleFac = 1/(max(Width)+ max(X) - min(X));
yScaleFac = 1/(max(Height)+ max(Y) - min(Y));

Xpos = xScaleFac*(X-min(X));
Ypos = 0.9*yScaleFac*(Y-min(Y));

for k=1:length(chNum) - 2
  subplotOL('position',[Xpos(k) Ypos(k)+(Height(k)*yScaleFac) Width(k)*xScaleFac*2 Height(k)*yScaleFac*2])
  tmpcfg = keepfields(cfg, {'parameter', 'layout', 'renderer', 'visible'});
  if exist('tmpdata', 'var')
    tmpcfg.style      = 'straight';
    tmpcfg.marker     = 'off';
    try
      tmpcfg.refmarker = strmatch(Lbl(k), data.reflabel);
    catch
      tmpcfg.refmarker = strmatch(Lbl(k), data.label);
    end
    tmpcfg.interplimits = 'electrodes';
    if isfield(cfg, 'xparam')
      tmpcfg.xparam = cfg.xparam;
      tmpcfg.xlim   = xparam;
    else
      tmpcfg.xparam = 'freq';
      tmpcfg.xlim   = [k-0.5 k+0.5];
    end
    tmpcfg.refchannel = Lbl(k);
    tmpcfg.colorbar   = 'no';
    tmpcfg.zlim       = scale;
    tmpcfg.grid_scale = 30;
    ft_topoplotTFR(tmpcfg, data);
    drawnow
  end
end

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% set the figure window title
if ~isempty(dataname)
  set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), mfilename, join_str(', ', dataname)));
else
  set(gcf, 'Name', sprintf('%d: %s', double(gcf), mfilename));
end
set(gcf, 'NumberTitle', 'off');

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end