text:
% compute SVD of data and save to file
function [ops, U, Sv] = get_svdcomps(ops)

% load(sprintf('%s/%s/%s/regops_%s_%s_plane%d.mat', ops.ResultsSavePath, ops.mouse_name, ops.date, ...
%     ops.mouse_name, ops.date, ops.iplane))

iplane = ops.iplane;

[Ly, Lx] = size(ops.mimg);

ntotframes          = ceil(sum(ops.Nframes));
% number of frames used to compute SVD
ops.NavgFramesSVD   = min(ops.NavgFramesSVD, ntotframes);
% size of binning (in time)
nt0 = ceil(ntotframes / ops.NavgFramesSVD);

if isfield(ops, 'chunk_align') && ~isempty(ops.chunk_align); chunk_align   = ops.chunk_align(iplane);
else chunk_align = 1; end

if chunk_align>9
    nt0 =  ops.chunk_align;
end
ops.NavgFramesSVD = floor(ntotframes/nt0);
nimgbatch = nt0 * floor(2000/nt0);

%% load the data

ix = 0;
fid = fopen(ops.RegFile, 'r');
mov = zeros(numel(ops.yrange), numel(ops.xrange), ops.NavgFramesSVD, 'single');

while 1
    data = fread(fid,  Ly*Lx*nimgbatch, '*int16');
    if isempty(data)
        break;
    end
    data = reshape(data, Ly, Lx, []);

    % subtract off the mean of this batch
%     data = data - repmat(ops.mimg1, 1, 1, size(data,3));

    nSlices = nt0*floor(size(data,3)/nt0);
    if nSlices ~= size(data,3)
        data = data(:,:, 1:nSlices);
    end

    % bin data
    data = reshape(data, Ly, Lx, nt0, []);
    data = single(data);
    davg = squeeze(mean(data,3));

    mov(:,:,ix + (1:size(davg,3))) = davg(ops.yrange, ops.xrange, :);

    ix = ix + size(davg,3);
end
fclose(fid);
mov = mov(:, :, 1:ix);

%% SVD options

% number of SVD components kept
ops.nSVD = min(ops.nSVD, size(mov,3));
%
mov             = reshape(mov, [], size(mov,3));
% mov             = mov./repmat(mean(mov.^2,2).^.5, 1, size(mov,2));

% compute covariance matrix of frames
if ops.useGPU
    COV             = gpuBlockXtX(mov)/size(mov, 1);
else
    COV             = mov' * mov/size(mov,1);
end

ops.nSVD = min(size(COV,1)-2, ops.nSVD);

% take SVD of covariance matrix and keep ops.nSVD components
if ops.nSVD<1000 || size(COV,1)>1e4
    [V, Sv]          = eigs(double(COV), ops.nSVD);
else
    if ops.useGPU
        gpuCOV = gpuArray(double(COV));
        [V, Sv]         = svd(gpuCOV);
        clear gpuCOV;
        V = gather(single(V));
        Sv = gather(single(Sv));
    else
         [V, Sv]         = svd(COV);
    end
    V               = V(:, 1:ops.nSVD);
    Sv              = Sv(1:ops.nSVD, 1:ops.nSVD);
end
%%

% compute U (normalized spatial masks... pixels x components)
if ops.useGPU
    U               = normc(gpuBlockXY(mov, V));
else
    U               = normc(mov*V);
end
U               = single(U);
Sv              = single(diag(Sv));

if ~exist(ops.ResultsSavePath, 'dir')
    mkdir(ops.ResultsSavePath)
end

% project spatial masks onto raw data
fid = fopen(ops.RegFile, 'r');
ix = 0;
Fs = zeros(ops.nSVD, sum(ops.Nframes), 'single');
while 1
    data = fread(fid,  Ly*Lx*nimgbatch, '*int16');
    if isempty(data)
        break;
    end
    data = reshape(data, Ly, Lx, []);

    % subtract off the mean of this batch
    %         data = data - repmat(mean(data,3), 1, 1, size(data,3));
%     data = data - repmat(ops.mimg1, 1, 1, size(data,3));
    data = data(ops.yrange, ops.xrange, :);
    data = single(data);
    if ops.useGPU
        Fs(:, ix + (1:size(data,3))) = gpuBlockXtY(U, reshape(data, [], size(data,3)));
    else
        Fs(:, ix + (1:size(data,3))) = U' * reshape(data, [], size(data,3));
    end

    ix = ix + size(data,3);
end
fclose(fid);

if ~exist(ops.ResultsSavePath, 'dir')
    mkdir(ops.ResultsSavePath);
end

totF = [0 cumsum(ops.Nframes)];
for iexp = 1:length(ops.expts)
    Vcell{iexp} = Fs(:, (1+ totF(iexp)):totF(iexp+1));
end

%%% save SVDs
U = reshape(U, numel(ops.yrange), numel(ops.xrange), []);
try % this is faster, but is limited to 2GB files
    save(sprintf('%s/SVD_%s_%s_plane%d.mat', ops.ResultsSavePath, ...
        ops.mouse_name, ops.date, iplane), 'U', 'Sv', 'Vcell', 'ops', '-v6');
catch % this takes a bit less space, but is significantly slower
    save(sprintf('%s/SVD_%s_%s_plane%d.mat', ops.ResultsSavePath, ...
        ops.mouse_name, ops.date, iplane), 'U', 'Sv', 'Vcell', 'ops');
end

% keyboard;

author: cortex-lab


text:
function [ Q ] = UpdateSARSA( s, a, r, sp, ap, tab , alpha, gamma )
% UpdateQ update de Qtable and return it using Whatkins QLearing
% s1: previous state before taking action (a)
% s2: current state after action (a)
% r: reward received from the environment after taking action (a) in state
%                                             s1 and reaching the state s2
% a:  the last executed action
% tab: the current Qtable
% alpha: learning rate
% gamma: discount factor
% Q: the resulting Qtable

Q = tab;
Q(s,a) =  Q(s,a) + alpha * ( r + gamma*Q(sp,ap) - Q(s,a) );

author: epfl-lasa


text:
function [v, d, m] = scimat_pca(scimat)
% SCIMAT_PCA Principal Principal Component Analysis of the selected points
% in a SCIMAT segmentation mask.
%
% [V, D, M] = scimat_pca(SCIMAT)
%
%   This function computes Principal Component Analysis (PCA) on the
%   collection of points in a SCIMAT segmentation mask.
%
%   SCIMAT is a struct with the segmentation (see "help scimat" for
%   details).
%
%   V is a matrix with the column eigenvectors ordered in decreasing order
%   of the corresponding eigenvalues in D.
%
%   M is a 3-vector with the coordinates of the dataset centroid. The
%   centroid is the centre of the coordinate system defined by the
%   eigenvectors, and the centre of rotation if the data is rotated.
%
%   Assuming a 3D volume, V is a (3, 3)-matrix and D a 3-vector.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010,2014 University of Oxford
% Version: 0.2.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 1) ;
nargoutchk(0, 3) ;

% squeeze SCIMAT variables, if necessary
scimat = scimat_squeeze(scimat);

% extract linear indices of voxels in the segmentation
idx = find(scimat.data);

% get volume size
sz = size(scimat.data);

% convert linear index to multiple subscripts
[ir, ic, iz] = ind2sub(sz, idx);

% convert indices to real world coordinates
x = scimat_index2world([ir, ic, iz], scimat);

% compute centroid
m = mean(x);

% compute PCA
[v, d] = pts_pca(x');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function start_num = p00_start_num ( prob )

%*****************************************************************************80
%
%% P00_START_NUM returns the number of starting points for a problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the problem.
%
%    Output, integer START_NUM, the number of starting points.
%
  if ( prob == 1 )
    start_num = p01_start_num ( );
  elseif ( prob == 2 )
    start_num = p02_start_num ( );
  elseif ( prob == 3 )
    start_num = p03_start_num ( );
  elseif ( prob == 4 )
    start_num = p04_start_num ( );
  elseif ( prob == 5 )
    start_num = p05_start_num ( );
  elseif ( prob == 6 )
    start_num = p06_start_num ( );
  elseif ( prob == 7 )
    start_num = p07_start_num ( );
  elseif ( prob == 8 )
    start_num = p08_start_num ( );
  elseif ( prob == 9 )
    start_num = p09_start_num ( );
  elseif ( prob == 10 )
    start_num = p10_start_num ( );
  elseif ( prob == 11 )
    start_num = p11_start_num ( );
  elseif ( prob == 12 )
    start_num = p12_start_num ( );
  elseif ( prob == 13 )
    start_num = p13_start_num ( );
  elseif ( prob == 14 )
    start_num = p14_start_num ( );
  elseif ( prob == 15 )
    start_num = p15_start_num ( );
  elseif ( prob == 16 )
    start_num = p16_start_num ( );
  elseif ( prob == 17 )
    start_num = p17_start_num ( );
  elseif ( prob == 18 )
    start_num = p18_start_num ( );
  elseif ( prob == 19 )
    start_num = p19_start_num ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_START_NUM - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_START_NUM - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function O = renderCluster_ui4(O)
% :Usage:
% ::
%
%    O = renderCluster_ui4(O)
%
% :Inputs:
%
%   **O**
%
%   **O.head:**
%        name of head to use, no .img extension
%
%   **O.surf:**
%        'y' or 'n': image surface
%
%   **O.sets:**
%        number of sets
%
%   **O.getfrom:**
%        cell array of where to get clusters from: 'workspace' 'file' etc.
%
%   **O.varname:**
%        cell array of variable names
%
%   **O.color:**
%        cell array of colors - strings or vectors
%
%   **O.numbers:**
%        'y' or 'n': add numbers to plot
%
%   also need individual fields named contents of varname, which contain clusters

if ~isfield(O,'clusterAlpha'),O.clusterAlpha = .7;, end
if ~isfield(O,'coords'),O.coords = [30 30 -5];, end
if ~isfield(O,'whichcuts'),O.whichcuts = 'z';, end
if ~isfield(O,'dohead'),O.dohead = 1;, end

%------------------------------------------------------------------------------
% Get head, if not default
%------------------------------------------------------------------------------
cflag = 0;
while cflag == 0
    if ~isfield(O,'head')
        a = input('Use default head single_subj_T1...?  (type default or file name, no .img extension)','s');
    else a = O.head;
    end

    switch a
    case 'default', fileName = 'single_subj_T1';
    otherwise fileName = a;
    end
    if exist([fileName '.img']) == 2,
	cflag = 1;
    else
	disp(['Can''t locate file: ' fileName '.img'])
    O = rmfield(O,'head');
    end
end

%------------------------------------------------------------------------------
% Add surface?
%------------------------------------------------------------------------------
cflag = 0;
while cflag == 0
    if ~isfield(O,'surf')
        a = input('Add brain surface (y\n\filename,no .img extension)...?','s');
    else a = O.surf;
    end

    switch a
    case 'y', fileName2 = ['surf_' fileName];, dosurf = 1;
    case 'n', fileName2 = [];, cflag = 1;, dosurf = 0;
    otherwise fileName2 = a;, dosurf = 1;
    end
    if dosurf & exist([fileName2 '.mat']) == 2,
	    cflag = 1;
    elseif dosurf
	    disp(['Can''t locate file: ' fileName '.mat'])
        O = rmfield(O,'surf');
    end
end

if ~isfield(O,'sets')
    numSets = input('Enter number of cluster sets: ');
else numSets = O.sets;
end

%------------------------------------------------------------------------------
% Add numbers?
%------------------------------------------------------------------------------
cflag = 0;
while cflag == 0
    if ~isfield(O,'numbers')
        addNums = input('Add numbers to plot? ','s');
    else addNums = O.numbers;
    end
    cflag = 1;
end

for i = 1:numSets

    disp(['Cluster Set ' num2str(i)])
    disp('_______________________________________________')

    %------------------------------------------------------------------------------
    % Get clusters, if necessary
    %------------------------------------------------------------------------------
    cflag = 0;
    while cflag == 0
       if ~isfield(O,'getfrom')
            a = input('Get clusters from...?  (workspace,file,figure)','s');
        else
            a = O.getfrom{i};
        end

       switch a
       case 'workspace',
            cflag = 1;
            if ~isfield(O,'varname')
                varname = input('Enter name of workspace variable: ','s');
                eval(['clusters = ' varname])
            else
                varname = O.varname{i};
                eval(['clusters = O.' varname ';'])
            end

            if isempty(clusters),disp('No clusters or not a valid variable name.'),
                cflag=0;
                O = rmfield(O,'varname');
            end
       case 'file', clusters = imageCluster('getclusters');, cflag = 1;
       case 'figure', clusters = imageCluster('getfigclusters');, cflag = 1;
       otherwise disp(['Not a valid choice...type entry string with no quotes.'])
       end
    end

    if ~isfield(O,'color')
        mycolor = input('Enter color string in single quotes or color vector for clusters: ');
    else
        mycolor = O.color{i};
    end

    cH = [];
    for j = 1:length(clusters)
        cH(j) = imageCluster('cluster',clusters(j),'color',mycolor,'alpha',O.clusterAlpha); drawnow
        set(cH(j),'Tag',['set' num2str(i) 'cluster' num2str(j)])
    end

    clusterHandles{i} = cH;

    if strcmp(addNums,'y')
          numberHandles{i} = cl_line_plots(clusters);
    end

end

%------------------------------------------------------------------------------
% * lighting stuff
%------------------------------------------------------------------------------
lighting gouraud;
[az,el] = view;
myLight = lightangle(az,el);
set(myLight,'Tag','myLight')
rotate3d on

% set callback to light follow the camera
%------------------------------------------------------------------------------
if exist('lightFollowView') == 2
     set(gcf, 'WindowButtonUpFcn', 'lightFollowView');
else
     warning('Cannot find lightFollowView.m to set light position.')
end
axis vis3d
axis image
drawnow

%------------------------------------------------------------------------------
% * head surface
%------------------------------------------------------------------------------
if dosurf
    p = [];
    eval(['load ' fileName2])
    try
        p(end+1) = patch('Faces',faces,'Vertices',vertices,'FaceColor',[.5 .5 .5],'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',.5,'SpecularExponent',200);
    catch
        p(end+1) = patch('Faces',faces,'Vertices',vertices,'FaceColor',[.5 .5 .5],'EdgeColor','none','SpecularStrength',.2,'SpecularExponent',200);
    end
    axis vis3d
    axis image
    drawnow
    O.p = p;
end

O.clusterHandles = clusterHandles;
O.numberHandles = numberHandles;

if O.dohead
    [D,Ds,hdr,p,coords,X,Y,Z] = tor_3d('coords',O.coords,'whichcuts',O.whichcuts);
    set(gcf,'Color','k');axis off;colormap(gray)
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function out_spm_gii( GiiTessFile, GiiDataFile, Data )
% OUT_TESS_GII: Exports a source map to a set of SPM/GIfTI .gii/.dat file.
%
% USAGE:  out_spm_gii( GiiTessFile, GiiDataFile, Data )
%
% INPUT:
%    - GiiTessFile : Brainstorm tesselation matrix
%    - GiiDataFile : full path to output file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013

% Get file paths
[giiPath, giiBase, giiExt] = bst_fileparts(GiiTessFile);
[outPath, outBase, outExt] = bst_fileparts(GiiDataFile);
% If the two files are in the same folder: do not use the path for the .gii surface
if strcmp(giiPath, outPath)
    GiiTessFile = [giiBase, giiExt];
end
% .dat filename, to save the actual data
DatFileShort = [outBase, '.dat'];
DatFile = bst_fullfile(outPath, DatFileShort);

% ===== CREATE XML STRING =====
% Create XML file from a template
strXml = [...
    '<?xml version="1.0" encoding="UTF-8"?>' 13 10 ...
    '<!DOCTYPE GIFTI SYSTEM "http://www.nitrc.org/frs/download.php/115/gifti.dtd">' 13 10 ...
    '<GIFTI Version="1.0"  NumberOfDataArrays="1">' 13 10 ...
    '   <MetaData>' 13 10 ...
    '      <MD>' 13 10 ...
    '         <Name><![CDATA[SurfaceID]]></Name>' 13 10 ...
    '         <Value><![CDATA[' GiiTessFile ']]></Value>' 13 10 ...
    '      </MD>' 13 10 ...
    '   </MetaData>' 13 10 ...
    '   <LabelTable/>' 13 10 ...
    '   <DataArray  ArrayIndexingOrder="ColumnMajorOrder"' 13 10 ...
    '               DataType="NIFTI_TYPE_FLOAT32"' 13 10 ...
    '               Dim0="' num2str(size(Data,1)) '"' 13 10 ...
    '               Dim1="1"' 13 10 ...
    '               Dimensionality="2"' 13 10 ...
    '               Encoding="ExternalFileBinary"' 13 10 ...
    '               Endian="LittleEndian"' 13 10 ...
    '               ExternalFileName="' DatFileShort '"' 13 10 ...
    '               ExternalFileOffset="0"' 13 10 ...
    '               Intent="NIFTI_INTENT_NONE">' 13 10 ...
    '      <MetaData>' 13 10 ...
    '      </MetaData>' 13 10 ...
    '      <Data></Data>' 13 10 ...
    '   </DataArray>' 13 10 ...
    '</GIFTI>' 13 10];

% ===== SAVE XML FILE =====
% Open file for binary writing
[fid, message] = fopen(GiiDataFile, 'wb');
if (fid < 0)
    error(['Could not create file : ' message]);
end
% Write XML contents
fwrite(fid, strXml, 'char');
% Close file
fclose(fid);

% ===== SAVE DAT FILE =====
% Open file for binary writing
fid = fopen(DatFile,'w');
if (fid < 0)
    error(['Could not create file : ' message]);
end
% Write data vector
fwrite(fid, Data, 'float32');
% Close file
fclose(fid);

end

author: brainstorm-tools


text:
function fk_space = autoGen_fk_space(q1,q2,q3)
%AUTOGEN_FK_SPACE
%    FK_SPACE = AUTOGEN_FK_SPACE(Q1,Q2,Q3)

%    This function was generated by the Symbolic Math Toolbox version 8.4.
%    01-Jun-2020 11:59:03

t2 = cos(q1);
t3 = cos(q2);
t4 = sin(q1);
t5 = q2+q3;
t6 = cos(t5);
t7 = sin(t5);
t8 = t3.*(4.0./2.5e+1);
t9 = t6.*(5.7e+1./2.0e+2);
t10 = t8+t9;
fk_space = reshape([t2.*t6,t4.*t6,-t7,0.0,-t4,t2,0.0,0.0,t2.*t7,t4.*t7,t6,0.0,t2.*t10,t4.*t10,t7.*(-5.7e+1./2.0e+2)-sin(q2).*(4.0./2.5e+1),1.0],[4,4]);

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function calpak_test40 ( )

%*****************************************************************************80
%
%% CALPAK_TEST40 tests MONTH_TO_MONTH_NAME_HEBREW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST40\n' );
  fprintf ( 1, '  For the Hebrew calendar,\n' );
  fprintf ( 1, '  MONTH_TO_MONTH_NAME_HEBREW names the months.\n' );
  fprintf ( 1, '\n' );

  y = 1;
  months = year_length_months_hebrew ( y );

  for i = 1 : months
    m = i;
    month_name = month_to_month_name_hebrew ( y, m );
    fprintf ( 1, '  %2d  %s\n', m, month_name );
  end

  return
end


text:
function computeModelCoefficients_batchLGG(pathExperiments,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeModelCoefficients_batchLGG(pathExperiments,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the final logistic regression coefficients and
% bootstrap confidence intervals of the final models obtained for all
% outcomes analyzed in the LGG study.
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to the directory containing all experiments.
%                     --> Ex: '/myProject/WORKSPACE/LOGISTIC_REGRESSION'
% 2. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (formal reference to come).
%               --> Ex: 'IALR'
% 3. nBatch: Number of parallel batch.
%            --> Ex: 8
% 4. matlabPATH: Full path to the MATLAB executable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% 5. seed: Numerical number to use as seed for bootstrapping experiment
%          --> Ex: 54288
% -------------------------------------------------------------------------
% OUTPUTS: Final coefficients, model response and confidence intervals
%          saved in a folder named '/myProject/WORKSPACE/LOGISTIC_REGRESSION/FINAL_MODELS/nameOutcome/fSetName'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
time = 30; % Number of seconds to wait before checking if parallel computations are done

cd(pathExperiments), load('training')
cd('FINAL_MODELS'), pathFinalModels = pwd;
mkdir('batchLog_Coeff'), cd('batchLog_Coeff'), pathBatch = pwd;
nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
for o = 1:nOutcomes
    outcomes.(nameOutcomes{o}) = training.(nameOutcomes{o}).outcome;
end
setNames = fieldnames(training.(nameOutcomes{1}).text);
[param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

% PRODUCE BATCH COMPUTATIONS
save('workspace','pathFinalModels','training','param','pathBatch','imbalance','seed'), pause(3)
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'tic\n');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:numel(param{i})
        fprintf(fid,['cd(fullfile(pathFinalModels,param{',num2str(i),'}{',num2str(j),'}{2},param{',num2str(i),'}{',num2str(j),'}{1}))\n']);
        fprintf(fid,'load(''finalModel'')\n');
        fprintf(fid,['fprintf(''COMPUTING THE LOGISTIC REGRESSION COEFFICIENTS OF THE FINAL MODEL OF "',param{i}{j}{2},'" OUTCOME, "',param{i}{j}{1},'" FEATURE SET ... '')']);
        fprintf(fid,'\n');
        fprintf(fid,['[coeff,response,modelCI] = computeModelCoefficients(finalModel.Data,training.',param{i}{j}{2},'.outcome,imbalance,seed);\n']);
        fprintf(fid,['fprintf(''DONE!\\n'')']);
        fprintf(fid,'\n');
        fprintf(fid,'save(''coeff'',''coeff''), save(''response'',''response''), save(''modelCI'',''modelCI'')\n');
    end
    fprintf(fid,'cd(pathBatch)\n');
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all\n');
    fprintf(fid,'toc');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function c=ref_fftreal(f)
%REF_FFTREAL  Reference FFTREAL
%
%  FFTREAL is computed by doing an FFT, and keeping only the positive
%  frequencies.

L=size(f,1);
c=fft(f);

c=c(1:floor(L/2)+1,:);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% plot(x,NC(2012,:),'o');
plot(x,dC(2012,:),'o');

repo: Aero-matlab


text:
function [varargout] = rotate(axis, angle, varargin)
%align set alignment of the objects in the block
%
%   [...] = rotate(axis, angle, obj <, obj> ...);
%
%   Rotates the corresponding gradinet object(s) about the given axis by
%   the specified amount. Gradients parallel to the rotation axis and
%   non-gradient objects are not affected.
%   Possible rotation axes are 'x', 'y' or 'z'.
%
%   Returns either a cell-array of objects if one return parameter is
%   provided or an explicit list of objects if multiple parameters are
%   given. Can be used directly as a parameter of seq.addBlock().
%
%   See also  Sequence.addBlock
%
%   Maxim Zaitsev <maxim.zaitsev@uniklinik-freiburg.de>

axes={'x', 'y', 'z'};

% cycle through the objects and rotate gradients non-parallel to the
% given rotation axis. Rotated gradients assigned to the same axis are then
% added together.

% first create indexes of the objects to be bypassed or rotated
irotate1=[];
irotate2=[];
ibypass=[];
axes2rot=axes(~strcmp(axes,axis));
if length(axes2rot)~=2
    error('incorrect axis specification');
end

if 1==length(varargin) && iscell(varargin{1})
    va=varargin{1}; % we need this to allow for nested mr.rotate() calls
else
    va=varargin;
end

for i=1:length(va)
    event = va{i};
    if isempty(event)
        continue;
    end
    if  (~strcmp(event.type,'grad') && ...
            ~strcmp(event.type,'trap')) || ...
            strcmp(event.channel, axis)%['g' axis]
        ibypass=[ibypass i];
    else
        if strcmp(event.channel, axes2rot(1)) %['g' axes2rot(1)]
            irotate1=[irotate1 i];
        else
            if (strcmp(event.channel, axes2rot(2))) %['g' axes2rot(2)]
                irotate2=[irotate2 i];
            else
                ibypass=[ibypass i]; % should never happen
            end
        end
    end
end

% now every gradient to be rotated generates two new gradients, one on the
% original axis and one on the other from the axes2rot list

rotated1=cell(1,length(irotate1)+length(irotate2));
rotated2=cell(1,length(irotate1)+length(irotate2));
max_mag=0; % measure of the relevant amplitude
for i=1:length(irotate1)
    g=va{irotate1(i)};
    max_mag=max(max_mag, getGradAbsMag(g));
    rotated1{i}=mr.scaleGrad(g,cos(angle));
    g=mr.scaleGrad(g,sin(angle));
    g.channel=axes2rot{2};
    rotated2{i}=g;
end
o=length(irotate1);
for i=1:length(irotate2)
    g=va{irotate2(i)};
    max_mag=max(max_mag, getGradAbsMag(g));
    rotated2{i+o}=mr.scaleGrad(g,cos(angle));
    g=mr.scaleGrad(g,-sin(angle));
    g.channel=axes2rot{1};
    rotated1{i+o}=g;
end

% eliminate zero-amplitude gradients
thresh=1e-6*max_mag;
for i=length(rotated1):-1:1
    if getGradAbsMag(rotated1{i})<thresh
        rotated1(i)=[];
    end
end
for i=length(rotated2):-1:1
    if getGradAbsMag(rotated2{i})<thresh
        rotated2(i)=[];
    end
end

g=cell(1,2);
% now we add gradients on the corresponding axis together
if (length(rotated1)>1)
    g{1}=mr.addGradients(rotated1);
else
    if (~isempty(rotated1))
        g{1}=rotated1{1};
    end
end

if (length(rotated2)>1)
    g{2}=mr.addGradients(rotated2);
else
    if (~isempty(rotated2))
        g{2}=rotated2{1};
    end
end

% eliminate zero-amplitude gradients
for i=length(g):-1:1
    if isempty(g{i}) || getGradAbsMag(g{i})<thresh
        g(i)=[];
    end
end

% export
bypass=va(ibypass);
out={bypass{:},g{:}};

nout = nargout;
varargout = cell(1,nout);
if nout==1
    varargout{1}=out;
else
    nr=min(nout,length(out));
    if nout<length(out)
        warning('insufficient number of return parameters, some rotated gradient components might go lost');
    end
    for k=1:nr
        varargout{k}=out{k};
    end
end

end

function [out] = getGradAbsMag(grad)
    if strcmp(grad.type,'trap')
        out=abs(grad.amplitude);
    else
        out=max(abs(grad.waveform));
    end
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [regions, scores] = sampleRegions(pB, spLR, seg, niter, maxov, nregions)

stats = regionprops(seg, 'Area');
area = [stats.Area];
area = area / sum(area);

updateFactor = 0.9;

nsp = numel(area);

scores = cell(niter+1,1);
regions = cell(niter+1, 1);

pB2 = pB;

% get regions through hierarchical segmentation
for k = 1:niter

    %pB2 = pB2-min(pB2);  pB2 = pB2 / max(pB2); % normalize 0 to 1
    hier = boundaries2hierarchy(pB2, spLR, 'max', [], [], pB);
    if k==1
        [scores{k}, scoresSp] = getHierarchyRegionScores(hier);
        scores{niter+1} = scoresSp;
        regions{niter+1} = num2cell(1:nsp)'; % superpixel regions
    else
        % assign costs according to original boundary likelihoods
%         for k2 = 1:nsp-1
%             hier.thresh(k2) = max(pB(hier.edges_removed{k2}));
%         end
        scores{k} = getHierarchyRegionScores(hier);
    end
    regions{k} = hier.new_region;

%     if k>1
%         ts = [scores{k} ; scores{1}];
%         [sval, sind] = sort(ts, 'descend');
%         sum(sind(1:50)<=numel(scores{k}))
%     end

    % decrement likelihoods of boundaries that are removed late
    for k2 = 1:nsp-1
        pB2(hier.edges_removed{k2}) = pB2(hier.edges_removed{k2}) * updateFactor^(k2/(nsp-1));
        %pB2(hier.edges_removed{k2}) = pB2(hier.edges_removed{k2}) + rand(1)*0.05;
    end
end

iterind = cell(numel(scores), 1);
for k = 1:numel(scores)
    iterind{k} = k*ones(numel(scores{k}), 1);
end
iterind = cat(1, iterind{:});

scores = cat(1, scores{:});
regions = cat(1, regions{:});

for k = 1:numel(regions)
    scores(k) = scores(k) * sqrt(sum(area(regions{k})));
end

if isempty(nregions)
    return;
end

coverage = zeros(nsp, 1);
cov = zeros(nregions, 1);
isselected = false(numel(regions), 1);

[scores, sind] = sort(scores, 'descend');
regions = regions(sind);

% take top nregion regions that do not overlap too much
rind = zeros(nregions,1);
n = 0;
for k1 = 1:numel(regions)
    isov = false;
    for k2 = 1:k1-1
        if isselected(k2) && (regionOverlap(regions{k1}, regions{k2}, area) > maxov)
            isov = true;
            break;
        end
    end
    if ~isov
        n=n+1;
        coverage(regions{k1}) = coverage(regions{k1}) + 1;
        isselected(k1) = true;
        cov(n) = sum(area(coverage>0))/sum(area);
        rind(n) = k1;
        if n==nregions
            break;
        end
    end
end
%regions = regions(rind(1:n));

iterind = iterind(sind);
for k = 1:max(iterind)
    disp([num2str(k) ': ' num2str(sum(iterind(rind)==k))]);
end

rind2 = zeros(sum(coverage==0), 1);
n2 = 0;
if any(coverage==0)
    for k = 1:numel(regions)
        if any(coverage(regions{k})==0)
            n2 = n2 +1;
            rind2(n2) = k;
            coverage(regions{k}) = coverage(regions{k})+1;
            if ~any(coverage==0)
                break;
            end
        end
    end
end

regions = regions([rind(1:n) ; rind2(1:n2)]);
scores = scores([rind(1:n) ; rind2(1:n2)]);

disp(mean(scores))

author: Cloud-CV

repo: object-proposals


text:
function cfg = save_large_cfg_fields(cfg, highest_ft, reproducescript_dir, nowtime)

% SAVE_LARGE_CFG_FIELDS is a helper function for ft_postamble_savevar and ft_postamble_savefig, and
% is used for the cfg.reproducescript functionality.

fn = ignorefields('recursesize');
for i=1:numel(fn)
  if isfield(cfg.callinfo.usercfg, fn{i}) && isstruct((cfg.callinfo.usercfg.(fn{i})))
    outputfile = make_or_fetch_inputfile(reproducescript_dir,...
      sprintf('%s_%s_input_%s.mat', nowtime, highest_ft, ...
      fn{i}), fn{i}, cfg.callinfo.usercfg.(fn{i}));
    cfg.callinfo.usercfg.(fn{i})  = outputfile;
  end
end

end

author: spm


text:
function [ali, vocab] = readKaldiAlignmentText_struct(file_name)
lines = my_cat(file_name);

vocab = [];
for i=1:length(lines)
    PrintProgress(i,length(lines),1000);
    idx = regexp(lines{i}, ' ');
    uttName = lines{i}(1:idx(1)-1);
    ali.(['U_' uttName]) = str2num(lines{i}(idx(1)+1:end))';
    vocab = [vocab; ali.(['U_' uttName])];
    if mod(i,100)==0
        vocab = unique(vocab);
    end
end
vocab = unique(vocab);
end

author: singaxiong

repo: SignalGraph


text:
% SYNTAX:
% status = hmrE_tCCA(sessIdx, sessName, derivedFolder)
%
% UI NAME:
% hmrE_tCCA
%
% DESCRIPTION:
%
% INPUT:
% sessIdx       -
% sessName      - the name of the session in a subject
% derivedFolder -
%
% OUTPUTS:
% status - error status
%
% USAGE OPTIONS:
% tCCA_filter: status = hmrE_tCCA(iSess, name, outputDirname)
%
function status = hmrE_tCCA(sessIdx, sessName, derivedFolder)

status = 0;

filenamePrefixNum = sprintf('tCCAfilter_%d', sessIdx);
filenamePrefixName = sprintf('tCCAfilter_%s', sessName);

files = dir(['./', filenamePrefixNum, '*']);
try
    for ii = 1:length(files)
        if files(ii).isdir
            continue
        end
        k = findstr(filenamePrefixNum, files(ii).name);
        suffix = files(ii).name(k+length(filenamePrefixNum):end);
        if ~exist(sessName, 'dir')
            movefile(files(ii).name, [filenamePrefixName, suffix]);
            continue;
        end
        movefile(files(ii).name, [derivedFolder, '/', sessName, '/', filenamePrefixName, suffix]);
    end
catch
    status = -1;
end

author: BUNPC


text:
function plot_continuous_samples(SAMPLES,Y)

% plot results of MCMC sampler
% The mean calcium sample, spike sampler raster plot and samples for
% amplitude, number of spikes, discrete time constants, noise variance,
% baseline and initial concentration are generated, together with their
% autocorrelation functions. If the marginalized flag was used, then the
% posterior pdfs of baseline and initial concentration are plotted.

% Inputs:
% SAMPLES:  structure with SAMPLES obtained from cont_ca_sampler.m
% Y:        inpurt fluorescence trace

% Author: Eftychios A. Pnevmatikakis, 2016, Simons Foundation

T = length(Y);
N = length(SAMPLES.ns);
show_gamma = 1;
P = SAMPLES.params;
P.f = 1;
g = P.g(:);
p = min(length(g),2);
Dt = 1/P.f;                                     % length of time bin
if ~isfield(SAMPLES,'g');
    show_gamma = 0;
    SAMPLES.g = ones(N,1)*g';
end

if p == 1
    tau_1 = 0;
    tau_2 = -Dt/log(g);                         % continuous time constant
    G1 = speye(T);
    G2 = spdiags(ones(T,1)*[-g,1],[-1:0],T,T);
    ge = P.g.^((0:T-1)');
elseif p == 2
    gr = roots([1,-g']);
    p1_continuous = log(min(gr))/Dt;
    p2_continuous = log(max(gr))/Dt;
    tau_1 = -1/p1_continuous;                   %tau h - smaller (tau_d * tau_r)/(tau_d + tau_r)
    tau_2 = -1/p2_continuous;                   %tau decay - larger
    G1 = spdiags(ones(T,1)*[-min(gr),1],[-1:0],T,T);
    G2 = spdiags(ones(T,1)*[-max(gr),1],[-1:0],T,T);
    ge = G2\[1;zeros(T-1,1)];
else
    error('This order of the AR process is currently not supported');
end

if length(SAMPLES.Cb) == 2
    marg = 1;       % marginalized sampler
else
    marg = 0;       % full sampler
end

C_rec = zeros(N,T);
for rep = 1:N
    %trunc_spikes = ceil(SAMPLES.ss{rep}/Dt);
    tau = SAMPLES.g(rep,:);
    gr = exp(-1./tau);
    ge = max(gr).^(0:T-1)';
    s_1 =   sparse(ceil(SAMPLES.ss{rep}/Dt),1,exp((SAMPLES.ss{rep} - Dt*ceil(SAMPLES.ss{rep}/Dt))/tau(1)),T,1);
    s_2 =   sparse(ceil(SAMPLES.ss{rep}/Dt),1,exp((SAMPLES.ss{rep} - Dt*ceil(SAMPLES.ss{rep}/Dt))/tau(2)),T,1);
    if gr(1) == 0
        G1 = sparse(1:T,1:T,Inf*ones(T,1)); G1sp = zeros(T,1);
    else
        G1 = spdiags(ones(T,1)*[-min(gr),1],[-1:0],T,T); G1sp = G1\s_1(:);
    end
    G2 = spdiags(ones(T,1)*[-max(gr),1],[-1:0],T,T);
    Gs = (-G1sp+ G2\s_2(:))/diff(gr);
    if marg
        %C_rec(rep,:) = SAMPLES.Cb(1) + SAMPLES.Am(rep)*filter(1,[1,-SAMPLES.g(rep,:)],full(s_)+[SAMPLES.Cin(:,1)',zeros(1,T-p)]);
        C_rec(rep,:) = SAMPLES.Cb(1) + SAMPLES.Am(rep)*Gs + (ge*SAMPLES.Cin(:,1));
    else
        %C_rec(rep,:) = SAMPLES.Cb(rep) + SAMPLES.Am(rep)*filter(1,[1,-SAMPLES.g(rep,:)],full(s_)+[SAMPLES.Cin(rep,:),zeros(1,T-p)]);
        C_rec(rep,:) = SAMPLES.Cb(rep) + SAMPLES.Am(rep)*Gs + (ge*SAMPLES.Cin(rep,:)');
    end
end
Nc = 60;

if marg
    rows = 4;
else
    rows = 5;
end

figure;
    set(gcf, 'PaperUnits', 'inches','Units', 'inches')
    set(gcf, 'PaperPositionMode', 'manual')
    set(gcf, 'PaperPosition',[0,0, 14, 15])
    set(gcf, 'Position',[2,2, 14, 15])
    ha(1) = subplot(rows,4,[1:4]);plot(Dt*(1:T),Y); hold all; plot(Dt*(1:T),mean(C_rec,1),'linewidth',2);
        title('Calcium traces','fontweight','bold','fontsize',14)
        legend('Raw data','Mean sample');
    ha(2) = subplot(rows,4,[5:8]); imagesc((1:T)*Dt,1:N,samples_cell2mat(SAMPLES.ss,T));
        title('Spike raster plot','fontweight','bold','fontsize',14)
        linkaxes(ha,'x');
    subplot(rows,4,4+5); plot(1:N,SAMPLES.ns); title('# of spikes','fontweight','bold','fontsize',14)
    subplot(rows,4,4+6); plot(-Nc:Nc,xcov(SAMPLES.ns,Nc,'coeff')); set(gca,'XLim',[-Nc,Nc]);
        title('Autocorrelation','fontweight','bold','fontsize',14)

    if ~show_gamma
        subplot(rows,4,4+7); plot(1:N,SAMPLES.ld); title('Firing Rate','fontweight','bold','fontsize',14)
        subplot(rows,4,4+8); plot(-Nc:Nc,xcov(SAMPLES.ld,Nc,'coeff')); set(gca,'XLim',[-Nc,Nc])
        title('Autocorrelation','fontweight','bold','fontsize',14)
    else
        if gr(1) == 0
            subplot(rows,4,4+7);  plot(1:N,exp(-1./SAMPLES.g(:,2))); title('Decay Time Constant','fontweight','bold','fontsize',14);
            subplot(rows,4,4+8);  plot(-Nc:Nc,xcov(exp(-1./SAMPLES.g(:,2)),Nc,'coeff')); title('Autocorrelation','fontweight','bold','fontsize',14);
            set(gca,'XLim',[-Nc,Nc])
        else
            subplot(rows,4,4+7);  plot(1:N,exp(-1./SAMPLES.g)); title('Decay Time Constants','fontweight','bold','fontsize',14);
            g_cov = xcov(exp(-1./SAMPLES.g),Nc,'coeff');
            subplot(rows,4,4+8);  plot(-Nc:Nc,g_cov(:,[1,4])); title('Autocorrelation','fontweight','bold','fontsize',14);
            set(gca,'XLim',[-Nc,Nc])
        end
    end

    subplot(rows,4,4+9); plot(1:N,SAMPLES.Am); title('Spike Amplitude','fontweight','bold','fontsize',14)
    subplot(rows,4,4+10); plot(-Nc:Nc,xcov(SAMPLES.Am,Nc,'coeff')); set(gca,'XLim',[-Nc,Nc])
        title('Autocorrelation','fontweight','bold','fontsize',14)
    if marg
        xx = SAMPLES.Cb(1) + linspace(-4*SAMPLES.Cb(2),4*SAMPLES.Cb(2));
        subplot(4,4,15); plot(xx,normpdf(xx,SAMPLES.Cb(1),SAMPLES.Cb(2)));
            set(gca,'XLim',[xx(1),xx(end)])
            title('Marg. post. of baseline','fontweight','bold','fontsize',14)

        xx = SAMPLES.Cin(1) + linspace(-4*SAMPLES.Cin(2),4*SAMPLES.Cin(2));
        subplot(4,4,16); plot(xx,normpdf(xx,SAMPLES.Cin(1),SAMPLES.Cin(2)));
            set(gca,'XLim',[xx(1),xx(end)])
            title('Marg. post. of initial con','fontweight','bold','fontsize',14)
    else
        subplot(5,4,4+11); plot(1:N,SAMPLES.Cb); title('Baseline','fontweight','bold','fontsize',14)
        subplot(5,4,4+12); plot(-Nc:Nc,xcov(SAMPLES.Cb,Nc,'coeff')); set(gca,'XLim',[-Nc,Nc])
            title('Autocorrelation','fontweight','bold','fontsize',14)
        subplot(5,4,4+13); plot(1:N,SAMPLES.Cin); title('Initial Concentration','fontweight','bold','fontsize',14)
        xcov_Cin = xcov(SAMPLES.Cin,Nc,'coeff');
        subplot(5,4,4+14); plot(-Nc:Nc,xcov_Cin); set(gca,'XLim',[-Nc,Nc])
            title('Autocorrelation','fontweight','bold','fontsize',14)
        subplot(5,4,4+15); plot(1:N,SAMPLES.sn2); title('Noise variance','fontweight','bold','fontsize',14)
        subplot(5,4,4+16); plot(-Nc:Nc,xcov(SAMPLES.sn2,Nc,'coeff')); set(gca,'XLim',[-Nc,Nc])
            title('Autocorrelation','fontweight','bold','fontsize',14)
    end
    drawnow;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [b,stats,yadj] = repeated_ancova(X,Y,wicons,btwnnames,winames,ynames,varargin)
% Repeated measures ANCOVA with table and plot
% uses Robust IRLS
%
% :Usage:
% ::
%
%     [b,stats,yadj] = repeated_ancova(X,Y,wicons,btwnnames,winames,ynames,varargin)
%
% :Examples:
% ::
%
%    Y = rand(15,2);
%    X = Y + rand(15,2);
%    cons = [-1 1 -1 1; 1 1 -1 -1];  % placebo vs control, hot vs. warm
%
%    X = R.X(:,1:2);
%    Y = cl(2).CONTRAST.data;
%    cons = [-1 1 0 0; 0 0 -1 1];  % placebo vs control for heat then warm
%    repeated_ancova(X,Y,cons,{'Reported Placebo (C - P)' 'Order'},{'P-C Heat' 'P-C Warm'},{'CH' 'PH' 'CW' 'PW'});
%
% DOES NOT WORK WITH FIXED BTWN-SUBJECTS COVARIATES
%
% X must be a random variable that is observed multiple times for each
% subject, as does Y
%
% ..
%    tor wager, Aug. 06
% ..

dotable = 1;
doplot = 1;

for i = 1:length(varargin)
    if isstr(varargin{i})
        switch varargin{i}
            % reserved keywords
            case 'notable', dotable = 0;
            case 'noplot', doplot = 0;

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

[n,r] = size(Y);

y = Y(:);

obs = size(y,1); % total num observations

% get predictors for between-subjects effects
% x
% ---------------------------------------------------
if all(size(X) - [n r] == 0)
    % X is obs. x conditions for a single predictor
    nbtwn = 1;
    x = X(:);   %repmat(X,r,1);
    btwnnames = repmat(btwnnames,1,r);

elseif size(X,1) == n
    disp('Warning!!!!!! Results are not valid for fixed within-subjects contrasts.')
    % X is a between-Ss effect replicated for each measure
    nbtwn = size(X,2);
    for i = 1:nbtwn
        x(:,i) = repmat(X(:,i),r,1);
    end
end

x = scale(x,1);  % center x

% get predictors for repeated measures effects
% c
% ---------------------------------------------------
% effect of repeated measure on y (repeated condition difference in ANCOVA)
% at x = 0
c = expand_contrasts(wicons,n,r);
nwithin = size(c,2);

% get preds for interaction between continuous x and c
% xc
% ---------------------------------------------------
% difference in x-y slopes for different repeated measures in c
if nbtwn
    [xc,xcnames,ninteract] = make_interactions(nbtwn,nwithin,x,c,n,r,btwnnames,winames);
else
    xc = []; xcnames = {}; ninteract = 0;
end

% Subject effects, obs x n-1

%nused = 1 + nbtwn ;  % warning!  not sure if this is right; only needed if btwn effects are replicated for all conditions
% THE ABOVE IS BAD; CANNOT ESTIMATE MODEL WHEN X IS FIXED.

S = repmat([eye(n-1); zeros(1,n-1)], r, 1);
%S = S(:,nused+1:end);

% y is regressed on x, controlling for S, and testing c

% full model
M = [x c xc S]; M(:,end+1) = 1;
p = size(M,2);              % total num params estimated
xtxitx = pinv(M);

res = y - M * xtxitx * y;   % residuals
SSr = res' * res;           % sums of squared residuals
%dfe = obs - p;              % error df
%dfb = p - 1;                % full model df

% reduced model: Subject effects only
Mr = S; Mr(:,end+1) = 1;
res = y - Mr * pinv(Mr) * y;   % residuals
SSr2 = res' * res;              % sums of squared residuals (2)

pr = size(Mr,2);              % total num params estimated, including intercept = n
dfe = obs - pr;              % error df
MSerr = SSr2 ./ dfe;

% model comparison
SSrdiff = SSr2 - SSr;         % sums of squares explained by model params x c xc
dfdiff = p - pr;              % parameters used to explain SSrdiff
MSmodel = SSrdiff ./ dfdiff;    % mean square for non-subject, non-intercept effects in model

F = MSmodel ./ MSerr;
p = 1 - fcdf(F,dfdiff,dfe);   % all this is very approx now!!! it's late.

[b,stats] = robustfit(M,y,'bisquare',[],'off');

wh_reg = 1;

if nargout > 2
    % get which regs are of interest for fitted data
    z = get_which_of_interest(wh_reg,nbtwn,nwithin,ninteract,S);

    yadj = get_yadj(wh_reg,x,y,M,z,b,n,r);
end

% table
if dotable
    regtable(b,stats,nbtwn,nwithin,ninteract,btwnnames,winames,xcnames);
end

% plot
if doplot
    %regplot(wh_reg,b,x,y,X,Y,M,S,nbtwn,nwithin,ninteract,n,r,c,btwnnames,winames,ynames);
    regplot2(X,Y,r,stats.w,btwnnames,ynames);
    ylabel('Y');

    z = get_which_of_interest(wh_reg,nbtwn,nwithin,ninteract,S);
    yadj = get_yadj(wh_reg,x,y,M,z,b,n,r);
    regplot2(X,yadj,r,stats.w,btwnnames,ynames);

end

return

% setup: get within-subjects contrasts

function c = expand_contrasts(cons,n,r)

ncons = size(cons,1);
c = zeros(n*r,ncons);

for i = 1:ncons
    mycon = zeros(n*r,1);
    for j = 1:r
        mycon((j-1)*n+1 : j*n) = cons(i,j);
    end
    c(:,i) = mycon;
end

return

function [xc,xcnames,ninteract] = make_interactions(nbtwn,nwithin,x,c,n,r,btwnnames,winames)
ninteract = nbtwn * nwithin;
xc = zeros(n*r,ninteract);
xccounter = 1;
for i = 1:nbtwn
    for j = 1:nwithin
        xc(:,xccounter) = x(:,i) .* c(:,j);
        xcnames{xccounter} = [btwnnames{i} '*' winames{j}];
        xccounter = xccounter + 1;
    end
end
return

function regtable(b,stats,nbtwn,nwithin,ninteract,btwnnames,winames,xcnames)
fprintf(1,'Parameter\tb-hat\tt\tp\t    \n')

fprintf(1,'Between-subjects effects    \n')
wh = 1:nbtwn;       % which columns for these effects
for i=1:nbtwn
    fprintf(1,'%s\t%3.2f\t%3.2f\t%3.4f\t    \n', ...
        btwnnames{i},b(wh(i)),stats.t(wh(i)),stats.p(wh(i)));
end

fprintf(1,'Within-subjects condition effects    \n')
wh = nbtwn+1:nbtwn+nwithin;
for i=1:nwithin
    fprintf(1,'%s\t%3.2f\t%3.2f\t%3.4f\t    \n', ...
        winames{i},b(wh(i)),stats.t(wh(i)),stats.p(wh(i)));
end

fprintf(1,'Between * Within interactions    \n')
wh = nbtwn+nwithin+1 : nbtwn+nwithin+ninteract;
for i=1:ninteract
    fprintf(1,'%s\t%3.2f\t%3.2f\t%3.4f\t    \n', ...
        xcnames{i},b(wh(i)),stats.t(wh(i)),stats.p(wh(i)));
end

return

function z = get_which_of_interest(wh_reg,nbtwn,nwithin,ninteract,S)
z = zeros(1,nbtwn);
z(wh_reg) = 1;              % between of interest
z = [z ones(1,nwithin)];    % all within are of interest
zxc = zeros(1,ninteract);   % interaction effects
xccounter = 1;
for i = 1:nbtwn
    for j = 1:nwithin
        if i == wh_reg, zxc(xccounter) = 1; end
        xccounter = xccounter + 1;
    end
end
z = [z zxc];
z = [z zeros(1,size(S,2))]; % subject effects
z = [z 1];                  % intercept
z = logical(z);
return

function yadj = get_yadj(wh_reg,x,y,M,z,b,n,r);
% x (of interest) and adjusted y
xdata = reshape(x(:,wh_reg),n,r);
yadj = y - M(:,~z) * b(~z);
yadj = reshape(yadj,n,r);

function regplot2(X,Y,r,w,btwnnames,ynames)

tor_fig; set(gca,'FontSize',24)
colors = {'ro' 'gs' 'cv' 'mh'};
for i = 1:r
    h = plot_correlation(X(:,i),Y(:,i),'colors',colors(i),'robust','noprint','weights',w);
    han(i) = h{1}(1);
end

han2 = findobj(gcf,'Type','text');
delete(han2)

% range over which to get fits for line
xrange = max(X(:,1)) - min(X(:,1));
minx = min(X(:,1)) - .1*xrange;
maxx = max(X(:,1)) + .1*xrange;

set(gca,'XLim',[minx maxx]);

% Y range
yrange = max(Y(:)) - min(Y(:));
miny = min(Y(:)) - .1*yrange;
maxy = max(Y(:)) + .1*yrange;
set(gca,'YLim',[miny maxy]);

xlabel(btwnnames{1});
ylabel('Adjusted data');
%ynames = {'HC' 'HA' 'LC' 'LA'};
legend(han,ynames)

scn_export_papersetup(400);

return

% OLD FUNCTION--but a goodie.  this was for fixed regressors btwn.
function regplot(wh_reg,b,x,y,X,Y,M,S,nbtwn,nwithin,ninteract,n,r,c,btwnnames,winames,ynames)

npar = nbtwn+nwithin+ninteract;

% only significant betas
% sig = stats.p < .10;
% sig = sig(1:npar);
% bsig = b(1:npar) .* sig;
% bsig = [bsig; b(end)];

% all betas for btwn of interest and w/i conditions
% bsig = zeros(1,npar)';
% bsig(wh_reg) = b(wh_reg);  % regressor of interest
% wh = nbtwn+1:nbtwn+nwithin; % condition, within Ss effects
% bsig(wh) = b(wh);
% wh = nbtwn+nwithin+1:npar; % interactions with reg of interest
% wh = wh(1:2);           %%%%%kludgy fix for reg 1 only
% bsig(wh) = b(wh);
% bsig = [bsig; b(end)];

% range over which to get fits for line
xrange = max(X(:,1)) - min(X(:,1));
minx = min(X(:,1)) - .1*xrange;
maxx = max(X(:,1)) + .1*xrange;

% get predicted values for this btwn effect (xf)
x2 = linspace(minx,maxx,n)';
x2 = repmat(x2,1,r);
x2 = x2(:);
xf = x;
xf(:,wh_reg) = x2;

% get interactions
[xcf,xcnamesf,ninteract] = make_interactions(nbtwn,nwithin,xf,c,n,r,btwnnames,winames);

% Model for pred values
Mf = [xf c xcf S];
Mf(:,end+1) = 1;

% get which regs are of interest for fitted data
z = get_which_of_interest(wh_reg,nbtwn,nwithin,ninteract,S);

%fit lines
f = Mf(:,z) * b(z);          % fitted data (for lines)
f = reshape(f,n,r);         % fits
xfit = reshape(x2,n,r);   % x values for fits

% x (of interest) and adjusted y
xdata = X;
%xdata = reshape(x(:,wh_reg),n,r);
yadj = y - M(:,~z) * b(~z);
yadj = reshape(yadj,n,r);

% make figure
tor_fig; set(gca,'FontSize',24)
colors = {'ro' 'bs' 'kv' 'gp'};
for i = 1:r
    plot(xdata(:,i),yadj(:,i),colors{i},'MarkerFaceColor',colors{i}(1));  % adjusted data
    han(i) = plot(xfit(:,i),f(:,i),[colors{i}(1) '-']);    % fit line
end

xlabel(btwnnames{wh_reg});
ylabel('Adjusted data');
%ynames = {'HC' 'HA' 'LC' 'LA'};
legend(han,ynames)

set(gca,'XLim',[minx maxx]);

% Y range
yrange = max(yadj(:)) - min(yadj(:));
miny = min(yadj(:)) - .1*yrange;
maxy = max(yadj(:)) + .1*yrange;
set(gca,'YLim',[miny maxy]);

scn_export_papersetup(400);

return


text:
%% Template Matlab script to create an BIDS compatible participants.tsv file
% This example lists all required and optional fields.
% When adding additional metadata please use CamelCase
%
% DHermes, 2017

% Writing json files relies on the JSONio library
% https://github.com/bids-standard/bids-matlab
%
% Make sure it is in the matab/octave path
try
    bids.bids_matlab_version;
catch
    warning('%s\n%s\n%s\n%s', ...
            'Writing the JSON file seems to have failed.', ...
            'Make sure that the following library is in the matlab/octave path:', ...
            'https://github.com/bids-standard/bids-matlab');
end

%%
clear;

this_dir = fileparts(mfilename('fullpath'));
root_dir = fullfile(this_dir, '..', filesep, '..');

project_label = 'templates';

participants_tsv_name = fullfile(root_dir, project_label, 'participants.tsv');

%% make a participants table and save

tsv.participant_id = {'sub-01'; 'sub-epilepsy01'};
tsv.age = [0; 10];
tsv.sex = {'m'; 'f'};
tsv.handedness = {'l'; 'r'};

bids.util.tsvwrite(participants_tsv_name, tsv);

%% associated data dictionary

template = struct( ...
                  'LongName', '', ...
                  'Description', '', ...
                  'Levels', struct(), ...
                  'Units', '', ...
                  'TermURL', 'http://www.todo.com/fixme');

json.age = template;
json.age.Description = 'age of the participant';
json.age.Units = 'years';

json.sex = template;
json.sex.Description = 'sex of the participant as reported by the participant';
json.sex.Levels = struct( ...
                         'm', 'male', ...
                         'f', 'female');

json.handedness = template;
json.handedness.Description = 'handedness of the participant as reported by the participant';
json.handedness.Levels = struct( ...
                                'l', 'left', ...
                                'r', 'right');

json_name = fullfile(root_dir, project_label, 'participants.json');

% Write JSON
% Make sure the directory exists
bids.util.mkdir(fileparts(json_name));
bids.util.jsonencode(json_name, json);

author: bids-standard

repo: bids-starter-kit


text:
classdef imtool3DROI_ellipse < imtool3DROI_rect

    properties (SetAccess = protected, GetAccess = protected)
        nPoints = 20;   %number of points to use to define the polygon that makes the elliptical mask
    end

    methods
        %constructor
        function ROI = imtool3DROI_ellipse(varargin)

            switch nargin
                case 0  %use the current figure

                    %let the user draw the ROI
                    h = imellipse;

                    %get the parent axes
                    ha = get(h,'Parent');

                    %get the handle of the image
                    hi = imhandles(ha);
                    if length(hi)>1
                        for i=1:length(hi)
                            if ndims(get(hi(i),'CData'))<3
                                imageHandle = hi(i);
                            end
                        end
                    else
                        imageHandle = hi;
                    end

                    %get the position
                    pos = getPosition(h);
                    position = [pos(1)+pos(3)/2 pos(2)+pos(4)/2 pos(3) pos(4)];

                    %delete the imroi object
                    delete(h);
                case 1 %user inputs only the handle to the image
                    imageHandle = varargin{1};
                    parent = get(imageHandle,'Parent');
                    h = imellipse(parent);
                    pos = getPosition(h);
                    position = [pos(1)+pos(3)/2 pos(2)+pos(4)/2 pos(3) pos(4)];
                    delete(h);
                case 2 %user inputs both the parent handle and a position
                    imageHandle = varargin{1};
                    position = varargin{2};
                case 3
                    imageHandle = varargin{1};
                    position = varargin{2};
                    if isempty(position)
                        parent = get(imageHandle,'Parent');
                        h = imellipse(parent);
                        pos = getPosition(h);
                        position = [pos(1)+pos(3)/2 pos(2)+pos(4)/2 pos(3) pos(4)];
                        delete(h);
                    end
                    tool = varargin{3};

            end

            %contruct the rect ROI
            if ~exist('tool','var'), tool=[]; end
            ROI@imtool3DROI_rect(imageHandle,position, tool)

            %make the rectangle an ellipse
            set(ROI.graphicsHandles(1),'Curvature',[1 1]);

            %update the position
            newPosition(ROI,position)

            %Set the button down functions of the graphics
            for i=1:length(ROI.graphicsHandles)
                fun = @(hObject,evnt) ButtonDownFunction(hObject,evnt,ROI,i); set(ROI.graphicsHandles(i),'ButtonDownFcn',fun);
            end

        end

        function newPosition(ROI,position)

            %set the position property of the ROI
            ROI.position = position;

            %find the top left corner of the box
            pos = [position(1)-position(3)/2 position(2)-position(4)/2 position(3) position(4)];

            %get the graphics handles
            graphicsHandles = ROI.graphicsHandles;

            %get the corner positions
            t = pi/4:pi/2:2*pi-pi/4;
            [x,y] = getEllipsePoints(position,t,'');

            %set the new position of the rectangle and other graphics
            %objects
            set(graphicsHandles(1),'Position',pos);
            set(graphicsHandles(2),'Xdata',position(1),'Ydata',position(2));
            set(graphicsHandles(3),'Xdata',x(3),'Ydata',y(3));
            set(graphicsHandles(4),'Xdata',x(4),'Ydata',y(4));
            set(graphicsHandles(5),'Xdata',x(2),'Ydata',y(2));
            set(graphicsHandles(6),'Xdata',x(1),'Ydata',y(1));
            set(graphicsHandles(7),'Xdata',pos(1),'Ydata',position(2));
            set(graphicsHandles(8),'Xdata',pos(1)+pos(3),'Ydata',position(2));
            set(graphicsHandles(9),'Xdata',position(1),'Ydata',pos(2));
            set(graphicsHandles(10),'Xdata',position(1),'Ydata',pos(2)+pos(4));

            %get the ROI measurements
            stats = getMeasurements(ROI);

            %set the textbox
            V = get(gca,'View');
            if V(1)==-90
                x = pos(1) + pos(3) + ROI.tbuff;
            else
                x = pos(1);
            end
            y = pos(2)-ROI.tbuff;

            str = {['Mean: ' num2str(stats.mean,'%+.2f')], ['STD:     ' num2str(stats.STD,'%.2f')]};
            set(ROI.textHandle,'String',str,'Position',[x y]);

            %notify a new position
            notify(ROI,'newROIPosition');
        end

        function [x, y] = getPoly(ROI)
            %get the position
            position = ROI.position;

            %find the top left corner of the box
            pos = [position(1)-position(3)/2 position(2)-position(3)/2 position(3) position(4)];

            %make the polygon
            [x,y] = getEllipsePoints(position,ROI.nPoints,'nPoints');
        end

        function stats = getMeasurements(ROI)
            [x, y] = getPoly(ROI);
            im = double(get(ROI.imageHandle,'CData'));

            m = size(im,1);
            n = size(im,2);
            %Scale the polygon to match the size of the displayed image (in
            %the case that the displayed image is being upsampled to match
            %the screen resolution).
            x = x*n/ROI.imageHandle.XData(2);
            y = y*m/ROI.imageHandle.YData(2);
            mask = poly2mask(x,y,m,n);

            im = im(mask);

            stats.mean = mean(im);
            stats.STD = std(im);
            stats.min = min(im);
            stats.max = max(im);
            stats.mask = mask;
            stats.position = ROI.position;

        end

        function varargout = autoCenterROI(ROI,varargin)
            Tforeground = .2;
            Tbackground = .8;
            pos=ROI.position(1:2);
            switch nargin
                case 1
                    weighting = 'pixels';
                    mode = 'normal';
                case 2
                    weighting = varargin{1}; %Weighting can be 'pixels' or 'mask'
                    mode ='normal'; %can be 'normal' or 'force positive contrast';
                case 3
                    weighting = varargin{1};
                    mode = varargin{2};
            end

            %Get the image
            im = double(get(ROI.imageHandle,'CData'));
            %Get coordinate system
            [X, Y]=meshgrid((1:size(im,2))-pos(1),(1:size(im,1))-pos(2));
            R=sqrt((X).^2 +(Y.^2));
            %get the mask
            stats = getMeasurements(ROI);
            mask=stats.mask;

            amin=min(im(mask));
            amax=max(im(mask));
            sim=mat2gray(im,[amin amax]);
            t=graythresh(sim(mask)); %Otsu threshold
            bw=im2bw(sim,t);
            bw=bw & mask;
            %bw=bwareaopen(bw,100);

            switch mode
                case 'force positive contrast'
                    obj=mean(im(bw));
                    bkg=mean(im(~bw & mask));
                case 'normal'
                    bkg=mean(im(mask & R>quantile(R(mask),Tbackground)));
                    obj=mean(im(mask & R<quantile(R(mask),Tforeground)));
                otherwise
                    bkg=mean(im(mask & R>quantile(R(mask),Tbackground)));
                    obj=mean(im(mask & R<quantile(R(mask),Tforeground)));
            end
            contrast=obj-bkg;
            if contrast<0
                bw=~bw & mask;
            end

            if any(bw(:))
                success = true;
                %compute the centroid of the rod
                center=zeros(1,2);
                switch weighting
                    case 'pixels'
                        center(1) = sum(X(bw).*im(bw))/sum(im(bw));
                        center(2) = sum(Y(bw).*im(bw))/sum(im(bw));
                    case 'mask'
                        center(1) = sum(X(bw).*bw(bw))/sum(bw(bw));
                        center(2) = sum(Y(bw).*bw(bw))/sum(bw(bw));
                    otherwise
                        center(1) = sum(X(bw).*im(bw))/sum(im(bw));
                        center(2) = sum(Y(bw).*im(bw))/sum(im(bw));
                end

                position = [center 0 0];
                position=position+ROI.position;
                newPosition(ROI,position);
            else
                success = false;
            end

            switch nargout
                case 1
                    varargout{1} = success;
            end

        end

        function success = autoCenterROIFindCircleMethod(ROI,circleRange)
            rescale=4;
            %Get the image
            im = double(get(ROI.imageHandle,'CData'));
            stats = getMeasurements(ROI);
            mask=stats.mask;

            %Crop the image
            [Yind Xind]=ind2sub(size(im),find(mask(:)));
            im=im(min(Yind):max(Yind),min(Xind):max(Xind));

            %Upsample the image
            im=imresize(im,rescale);

            %Find the circles
            warning('off','images:imfindcircles:warnForLargeRadiusRange')
            center = imfindcircles(im, round(rescale*circleRange),'Sensitivity',.95);
            warning('on','images:imfindcircles:warnForLargeRadiusRange')

            %move the ROI
            if ~isempty(center)
                center=center(1,:);
                success=true;
                center=center/rescale;
                center=center + [min(Xind) min(Yind)];
                position = [center ROI.position(3:4)];
                newPosition(ROI,position);
            else
                success=false;
            end

        end

    end

end

function [x,y] = getEllipsePoints(position,t,mode)
%This function returns a list of vertices of an elliptical polygon with
%nPoints number of vertices;
if strcmp(mode,'nPoints')
    t=linspace(0,2*pi,t); %elliptical equation is parameterized by t
end
a = position(3)/2;
b = position(4)/2;
x = a*cos(t); x = x+position(1);
y = b*sin(t); y = y+position(2);
end

function ButtonDownFunction(hObject,evnt,ROI,n)

%get the parent figure handle
fig = ROI.figureHandle;

%get the type of click
click = get(fig,'SelectionType');

if strcmp(click,'normal')
    %get the current button motion and button up functions of the figure
    WBMF_old = get(fig,'WindowButtonMotionFcn');
    WBUF_old = get(fig,'WindowButtonUpFcn');

    %set the new window button motion function and button up function of the figure
    fun = @(src,evnt) ButtonMotionFunction(src,evnt,ROI,n);
    fun2=@(src,evnt)  ButtonUpFunction(src,evnt,ROI,WBMF_old,WBUF_old);
    set(fig,'WindowButtonMotionFcn',fun,'WindowButtonUpFcn',fun2);
end
end

function ButtonMotionFunction(src,evnt,ROI,n)
cp = get(ROI.axesHandle,'CurrentPoint'); cp=[cp(1,1) cp(1,2)];

position = getPosition(ROI);

switch n
    case 2                                          %middle cross
        position(1) = cp(1); position(2) = cp(2);

    case 3                                          %top left corner
        %find the x and y for the current position
        [x,y] = getEllipsePoints(position,5*pi/4,'');
        dx = x-cp(1); dy = y-cp(2);
        cp(1) = position(1)-(position(3)/2+dx); cp(2) = position(2)-(position(4)/2+dy);

         %find the bottom edge
        bottom = position(2)+position(4)/2;
        %get the new height
        height = bottom - cp(2);
        if height>1
            cy = cp(2)+height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)+position(3)/2-height*ROI.aspectRatio;
            end
        end
        %find the right edge
        right = position(1)+position(3)/2;
        %find the new width
        width = right-cp(1);
        if width>1
            %find the new center
            cx = cp(1)+width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 4                                          %top right corner
        %find the x and y for the current position
        [x,y] = getEllipsePoints(position,7*pi/4,'');
        dx = cp(1)-x; dy = y-cp(2);
        cp(1) = position(1)+(position(3)/2+dx); cp(2) = position(2)-(position(4)/2+dy);
         %find the bottom edge
        bottom = position(2)+position(4)/2;
        %get the new height
        height = bottom - cp(2);
        if height>1
            cy = cp(2)+height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)-position(3)/2+height*ROI.aspectRatio;
            end
        end

        %find the left edge
        left = position(1)-position(3)/2;
        %find the new width
        width = cp(1) - left;
        if width>1
            cx = cp(1)-width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 5                                          %bottom left corner
        %find the x and y for the current position
        [x,y] = getEllipsePoints(position,3*pi/4,'');
        dx = x-cp(1); dy = cp(2)-y;
        cp(1) = position(1)-(position(3)/2+dx); cp(2) = position(2)+(position(4)/2+dy);
        %find the top edge
        top = position(2)-position(4)/2;
        %get the new height
        height = cp(2) - top;
        if height>1
            cy = cp(2)-height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)+position(3)/2+-height*ROI.aspectRatio;
            end
        end

        %find the right edge
        right = position(1)+position(3)/2;
        %find the new width
        width = right-cp(1);
        if width>1
            %find the new center
            cx = cp(1)+width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 6                                          %bottom right corner
        %find the x and y for the current position
        [x,y] = getEllipsePoints(position,pi/4,'');
        dx = cp(1)-x; dy = cp(2)-y;
        cp(1) = position(1)+(position(3)/2+dx); cp(2) = position(2)+(position(4)/2+dy);
        %find the top edge
        top = position(2)-position(4)/2;
        %get the new height
        height = cp(2) - top;
        if height>1
            cy = cp(2)-height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)-position(3)/2+height*ROI.aspectRatio;
            end
        end

        %find the left edge
        left = position(1)-position(3)/2;
        %find the new width
        width = cp(1) - left;
        if width>1
            cx = cp(1)-width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 7                                          %left
        %find the right edge
        right = position(1)+position(3)/2;
        %find the new width
        width = right-cp(1);
        if width>1
            %find the new center
            cx = cp(1)+width/2;
            position(1) = cx;
            position(3) = width;
        end
        if ROI.fixedAspectRatio
            position(4) = width/ROI.aspectRatio;
        end

    case 8                                          %right
        %find the left edge
        left = position(1)-position(3)/2;
        %find the new width
        width = cp(1) - left;
        if width>1
            cx = cp(1)-width/2;
            position(1) = cx;
            position(3) = width;
        end

        if ROI.fixedAspectRatio
            position(4) = width/ROI.aspectRatio;
        end

    case 9                                          %top
        %find the bottom edge
        bottom = position(2)+position(4)/2;
        %get the new height
        height = bottom - cp(2);
        if height>1
            cy = cp(2)+height/2;
            position(2) = cy;
            position(4) = height;
        end

        if ROI.fixedAspectRatio
            position(3) = height*ROI.aspectRatio;
        end

    case 10                                         %bottom
        %find the top edge
        top = position(2)-position(4)/2;
        %get the new height
        height = cp(2) - top;
        if height>1
            cy = cp(2)-height/2;
            position(2) = cy;
            position(4) = height;
        end

        if ROI.fixedAspectRatio
            position(3) = height*ROI.aspectRatio;
        end

end

newPosition(ROI,position);

end

function ButtonUpFunction(src,evnt,ROI,WBMF_old,WBUF_old)
fig = ROI.figureHandle;

set(fig,'WindowButtonMotionFcn',WBMF_old,'WindowButtonUpFcn',WBUF_old);

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function tapas_ehgf_jget_plotTraj(r)
% Plots the estimated trajectories for the HGF perceptual model for
% the JGET project
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_ehgf_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Time axis
if size(r.u,2) > 1 && ~isempty(find(strcmp(fieldnames(r.c_prc),'irregular_intervals'))) && r.c_prc.irregular_intervals
    t = r.u(:,end)';
else
    t = ones(1,size(r.u,1));
end

ts = cumsum(t);
ts = [0, ts];

% Do we know the generative parameters?
if size(r.u,2) > 2
    genpar = true;
    mean   = r.u(:,2);
    sd     = r.u(:,3);
else
    genpar = false;
end

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = length(r.p_prc.p)/8;
end

% Upper levels
for j = 1:l-1

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Left subplot (x)                       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(l+1,2,2*j-1);

    if plotsd == true
        upperprior = r.p_prc.mux_0(l-j+1) +1.96*sqrt(r.p_prc.sax_0(l-j+1));
        lowerprior = r.p_prc.mux_0(l-j+1) -1.96*sqrt(r.p_prc.sax_0(l-j+1));
        upper = [upperprior; r.traj.mux(:,l-j+1)+1.96*sqrt(r.traj.sax(:,l-j+1))];
        lower = [lowerprior; r.traj.mux(:,l-j+1)-1.96*sqrt(r.traj.sax(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mux_0(l-j+1); r.traj.mux(:,l-j+1)], 'b', 'LineWidth', 1.5);
    hold all;
    plot(0, r.p_prc.mux_0(l-j+1), 'ob', 'LineWidth', 1.5); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu x_', num2str(l-j+1)]);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Right subplot (alpha)                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(l+1,2,2*j);

    if plotsd == true
        upperprior = r.p_prc.mua_0(l-j+1) +1.96*sqrt(r.p_prc.saa_0(l-j+1));
        lowerprior = r.p_prc.mua_0(l-j+1) -1.96*sqrt(r.p_prc.saa_0(l-j+1));
        upper = [upperprior; r.traj.mua(:,l-j+1)+1.96*sqrt(r.traj.saa(:,l-j+1))];
        lower = [lowerprior; r.traj.mua(:,l-j+1)-1.96*sqrt(r.traj.saa(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mua_0(l-j+1); r.traj.mua(:,l-j+1)], 'b', 'LineWidth', 1.5);
    hold all;
    plot(0, r.p_prc.mua_0(l-j+1), 'ob', 'LineWidth', 1.5); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of \alpha_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu \alpha_', num2str(l-j+1)]);
end

% Input level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Left subplot (x)                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(l+1,2,2*l-1);

if plotsd == true
    upperprior = r.p_prc.mux_0(1) +1.96*sqrt(r.p_prc.sax_0(1));
    lowerprior = r.p_prc.mux_0(1) -1.96*sqrt(r.p_prc.sax_0(1));
    upper = [upperprior; r.traj.mux(:,1)+1.96*sqrt(r.traj.sax(:,1))];
    lower = [lowerprior; r.traj.mux(:,1)-1.96*sqrt(r.traj.sax(:,1))];

    plot(0, upperprior, 'or', 'LineWidth', 1);
    hold all;
    plot(0, lowerprior, 'or', 'LineWidth', 1);
    fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
plot(ts, [r.p_prc.mux_0(1); r.traj.mux(:,1)], 'r', 'LineWidth', 1.5);
hold all;
plot(0, r.p_prc.mux_0(1), 'or', 'LineWidth', 1.5); % prior
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0.6 0]); % inputs
if genpar
    plot(ts(2:end), mean, '-', 'Color', 'k', 'LineWidth', 1); % mean of input distribution
    plot(ts(2:end), mean +1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
    plot(ts(2:end), mean -1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
end
xlim([0 ts(end)]);
title(['Input u (green) and posterior expectation of x_1 (red) for \kappa_x=', ...
       num2str(r.p_prc.kax), ', \omega_x=', num2str(r.p_prc.omx)], 'FontWeight', 'bold');
ylabel('u, \mu x_1');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Right subplot (alpha)                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(l+1,2,2*l);

if plotsd == true
    upperprior = r.p_prc.mua_0(1) +1.96*sqrt(r.p_prc.saa_0(1));
    lowerprior = r.p_prc.mua_0(1) -1.96*sqrt(r.p_prc.saa_0(1));
    upper = [upperprior; r.traj.mua(:,1)+1.96*sqrt(r.traj.saa(:,1))];
    lower = [lowerprior; r.traj.mua(:,1)-1.96*sqrt(r.traj.saa(:,1))];

    transupperprior = sqrt(exp(r.p_prc.kau *upperprior +r.p_prc.omu));
    translowerprior = sqrt(exp(r.p_prc.kau *lowerprior +r.p_prc.omu));
    transupper = sqrt(exp(r.p_prc.kau *upper +r.p_prc.omu));
    translower = sqrt(exp(r.p_prc.kau *lower +r.p_prc.omu));

    plot(0, transupperprior, 'or', 'LineWidth', 1);
    hold all;
    plot(0, translowerprior, 'or', 'LineWidth', 1);
    fill([ts, fliplr(ts)], [(transupper)', fliplr((translower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
transmuaprior = sqrt(exp(r.p_prc.kau *r.p_prc.mua_0(1) +r.p_prc.omu));
plot(ts, [transmuaprior; sqrt(exp(r.p_prc.kau *r.traj.mua(:,1) +r.p_prc.omu))], 'r', 'LineWidth', 1.5);
hold all;
plot(0, transmuaprior, 'or', 'LineWidth', 1.5); % prior
if genpar
    plot(ts(2:end), sd, '--', 'Color', 'k', 'LineWidth', 1);
end
xlim([0 ts(end)]);
title(['Belief on noise (red) for \kappa_\alpha=', ...
       num2str(r.p_prc.kaa), ', \omega_\alpha=', num2str(r.p_prc.oma)], 'FontWeight', 'bold');
ylabel('\mu \alpha_1');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Decision model                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(l+1,2,2*l+1);

if plotsd == true
    upper = r.traj.muxhat(:,1)+1.96*sqrt(r.p_obs.ze +r.traj.saxhat(:,1));
    lower = r.traj.muxhat(:,1)-1.96*sqrt(r.p_obs.ze +r.traj.saxhat(:,1));

    fill([ts(2:end), fliplr(ts(2:end))], [(upper)', fliplr((lower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    hold all;
end
plot(ts(2:end), r.traj.muxhat(:,1), 'Color', [153/256 17/256 153/256], 'LineWidth', 1.5);
hold all;
plot(ts(2:end), r.y(:,1), '.', 'Color', [1 0.65 0], 'MarkerSize', 15); % responses
if genpar
    plot(ts(2:end), mean, '-', 'Color', 'k', 'LineWidth', 1); % mean of input distribution
    plot(ts(2:end), mean +1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
    plot(ts(2:end), mean -1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
end
xlim([1 ts(end)]);
title('Decision model: prediction of decision (purple) and decision (orange)', 'FontWeight', 'bold');
ylabel('y, \^{\mu} x_1');
xlabel({'Trial number', ' '}); % A hack to get the relative subplot sizes right
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Learning rate                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if genpar
subplot(l+1,2,2*l+2);
[AX, H1, H2 ] = plotyy(ts(2:end), mean, ts(2:end), r.traj.lrx(:,1));
hold all;
ylim(AX(1), [min(mean -1.96.*sd-3), max(mean +1.96.*sd+3)]);
plot(AX(1), ts(2:end), mean +1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1.1); % 95% interval of input distribution
plot(AX(1), ts(2:end), mean -1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1.1); % 95% interval of input distribution
set(H1, 'Color', 'k', 'LineWidth', 1.1);
set(H2, 'Color', [178/256, 34/256, 34/256], 'LineWidth', 1.5);
set(AX(1), 'YColor', 'k');
set(AX(2), 'YColor', 'k');
xlim(AX(1), [1 ts(end)]);
xlim(AX(2), [1 ts(end)]);
title('Learning rate (bordeaux) and input sampling distribution (black)', 'FontWeight', 'bold');
ylabel(AX(1), 'Input');
ylabel(AX(2), '\^{\pi}_u/\pi_x');
xlabel('Trial number');
hold off;
end

author: translationalneuromodeling

repo: tapas


text:
function Y = acosh(X)
    %ACOSH  Symbolic inverse hyperbolic cosine.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % construct the operation string
    sstr = ['ArcCosh[' X.s ']'];

    % create a new object with the evaluated string
    Y = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
%% FILE_NAME_INC increments a partially numeric filename.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the string to be incremented.
%
%    Output, string FILE_NAME, the incremented string.
%
  lens = s_len_trim ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c && c <= '8' )

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end
%
%  No digits were found.  Return blank.
%
  if ( change == 0 )
    file_name = ' ';
  end

  return
end


text:
function h = mrViewer(mr, format, varargin)
% mrVista general-purpose MRI data viewer.
%
%   h = mrViewer([data file or mr struct], [format], [options]);
%
% This tool allows you to visualize various forms of 2D,  3D,  or
% 4D MRI data,  including structural and functional data.
%
% Returns a handle to the viewer UI. Info about the UI,  including
% the base MR object and any maps,  are contained in this object's
% UserData. To get the viewer struct alone,  without initializing a
% GUI,  use mrViewInit.
%
% Optional additional arguments:
%   'session', [val]:   session struct from which the data are taken,
%                       or handle to a session viewer. Enables saving
%                       ROIs drawn using the viewer to the session.
%
%   'tSeries', [val]:   paths to time series file/files (string or cell
%                       of strings),  to associate with the data.
%
%   'stim', [val]:      stimulus files (.par files) to load and associate
%                       with the tSeries. Enables plotting time courses of
%                       ROIs drawn with the viewer.
%
%   'dock':             always dock all uipanels to a single figure,
%                       rather than the default of putting some panels
%                       in their own figures.
%
%
% ras,  07/05.

% Code notes:
% this m-file opens the mrViewer figure. Other files refresh
% or perform operations on it.

if notDefined('format'),  format = ''; end

% if MRI data not specified,  prompt via mrLoad:
if notDefined('mr'),   mr = mrLoad([],  format);
elseif ischar(mr),     mr = mrLoad(mr,  format);
end

if ~ispref('VISTA', 'dockFlag'), setpref('VISTA', 'dockFlag', 1); end
dockFlag = getpref('VISTA', 'dockFlag');

% initialize the UI struct
ui = mrViewInit(mr);
clear mr;

% parse options here
for i = 1:length(varargin)
    switch lower(varargin{i})
        case 'dock',  dockFlag = 1;
        case 'undock',  dockFlag = 0;
        case 'session',  ui.session = varargin{i+1};
        case 'tseries',  ui = mrViewAttachTSeries(ui,  varargin{i+1});
        case 'stim',  ui = mrViewAttachStim(ui,  varargin{i+1});
    end
end

% open the figure
ui.fig = figure('Units', 'Normalized', ...
    'Position', [.12 .23 .4 .4], ...
    'UserData', ui, ...
    'Color', 'k', ...
    'NumberTitle', 'off', ...
    'MenuBar', 'none', ...
    'NextPlot', 'add', ...
    'Name', sprintf('%s [%s]', ui.mr.name, ui.tag), ...
    'CloseRequestFcn', sprintf('mrViewClose(''%s'');', ui.tag), ...
    'Tag', ui.tag);

% add main panel for axes
ui.panels.display = uipanel('Position', [0 0 1 1], ...
                           'BackgroundColor', [0 0 0], ...
                           'ShadowColor', [0 0 0], ...
                           'BorderType', 'none');       % [0 0 .4]

% add control panels
ui = mrViewNavPanel(ui, dockFlag);
ui = mrViewGrayscalePanel(ui, dockFlag);
ui = mrViewROIPanel(ui, 1); % dockFlag
ui = mrViewInfoPanel(ui, 1);
ui = mrViewMeshPanel(ui, 1);
ui = mrViewColorbarPanel(ui, 1);

% add menus
ui = mrViewFileMenu(ui);
ui = mrViewEditMenu(ui);
ui = mrViewViewMenu(ui);
ui = mrViewPlotMenu(ui);
% uimenu(ui.fig, 'Label', 'Window', 'Callback', winmenu('callback'));
ui = mrViewSpaceMenu(ui);
ui = mrViewMeshMenu(ui);
uimenu('Label', '     '); % spacer
ui = mrViewHelpMenu(ui);

% set the initial coordinate space to pixel space
ui = mrViewSet(ui, 'space', ui.settings.space);

% % guess a good initial clip value for the anatomies
% [binCnt binCenters] = hist(ui.mr.data(:), 100);
% clipvals = binCenters( mrvMinmax(find(binCnt>histThresh) ));
% ui = mrViewSetGrayscale(ui, 'clip', clipvals);

% refresh view
ui = mrViewRefresh(ui);

h = ui.fig;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Armstrong no

function n= is_armstrong(num)
if(num<0)
    n=0; %as negative no can't be armstrong
end
sum=0;
temp=num;
while(temp>0)
    sum=sum*10 + rem(temp,10);
    temp=fix(temp/10);
end
if(sum==num)
    n=1;    % the no is armstrong
else
    n=0;    % the no is not armstrong
end
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function h=convolve(f,g,varargin)
%CONVOLVE  Convolution
%   Usage:  h=convolve(f,g);
%
%   `convolve` has been deprecated. Please use |lconv| instead.
%
%   A call to `convolve(f,g)` can be replaced by ::
%
%     lconv(f,g);
%
%   See also: lconv

warning(['LTFAT: CONVOLVE has been deprecated, please use LCONV ' ...
         'instead. See the help on LCONV for more details.']);

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.L=[];
definput.keyvals.dim=[];

[flags,kv,L,dim]=ltfatarghelper({'L','dim'},definput,varargin);

[f,L1,Lf,Wf,dimout,permutedsize_f,order_f]=assert_sigreshape_pre(f,L,dim,'CONVOLVE');
[g,L2,Lg,Wg,dimout,permutedsize_g,order_g]=assert_sigreshape_pre(g,L,dim,'CONVOLVE');

Lh=Lf+Lg-1;

if (Wf>1) && (Wg>1)
  error('%s: Only one of the inputs can be multi-dimensional.',upper(mfilename));
end;

W=max(Wf,Wg);
if Wf<W
  f=repmat(f,1,W);
end;

if Wg<W
  g=repmat(g,1,W);
end;

if isreal(f) && isreal(g)
  h=comp_ifftreal(comp_fftreal(postpad(f,Lh)).*...
                  comp_fftreal(postpad(g,Lh)),Lh);
else
  h=ifft(fft(postpad(f,Lh)).*...
         fft(postpad(g,Lh)));

end;


text:
function krn = spm_smoothkern(fwhm,x,t)
% Generate a Gaussian smoothing kernel
% FORMAT krn = spm_smoothkern(fwhm,x,t)
% fwhm - full width at half maximum
% x    - position
% t    - either 0 (nearest neighbour) or 1 (linear).
%        [Default: 1]
%
% krn  - value of kernel at position x
%__________________________________________________________________________
%
% For smoothing images, one should really convolve a Gaussian with a sinc
% function. For smoothing histograms, the kernel should be a Gaussian
% convolved with the histogram basis function used. This function returns
% a Gaussian convolved with a triangular (1st degree B-spline) basis
% function (by default). A Gaussian convolved with a hat function (0th
% degree B-spline) can also be returned.
%
% Note that the convolution kernel returned by this function differ from
% the ones that other packages currently use for Gaussian smoothing -
% particularly when the FWHM is small compared with the voxel dimensions.
% The fact that SPM does it differently from other software does not mean
% that it is wrong.
%__________________________________________________________________________
% Copyright (C) 2005-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_smoothkern.m 7460 2018-10-29 15:55:12Z john $

if nargin<3, t = 1; end

% Variance from FWHM
s = (fwhm/sqrt(8*log(2)))^2+eps;

% The simple way to do it. Not good for small FWHM
% krn = (1/sqrt(2*pi*s))*exp(-(x.^2)/(2*s));

if t==0
    % Gaussian convolved with 0th degree B-spline
    % int(exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= -0.5..0.5)
    w1  = 1/sqrt(2*s);
    krn = 0.5*(erf(w1*(x+0.5))-erf(w1*(x-0.5)));
    krn(krn<0) = 0;

elseif t==1
    % Gaussian convolved with 1st degree B-spline
    %  int((1-t)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= 0..1)
    % +int((t+1)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t=-1..0)
    w1  =  0.5*sqrt(2/s);
    w2  = -0.5/s;
    w3  = sqrt(s/2/pi);
    krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) - 2*erf(w1*x   ).* x)...
          +w3*(exp(w2*(x+1).^2)     + exp(w2*(x-1).^2)     - 2*exp(w2*x.^2));
    krn(krn<0) = 0;

else
    error('Only defined for nearest neighbour and linear interpolation.');
    % This should probably be based on https://arxiv.org/pdf/1608.05854.pdf:
    %     Martin TB, Prunet S, Drissen L. Optimal fitting of Gaussian-apodized
    %     or under-resolved emission lines in Fourier transform spectra
    %     providing new insights on the velocity structure of NGC 6720. Monthly
    %     Notices of the Royal Astronomical Society. 2016 Sep 14;463(4):4223-38.
    % (thanks for the pointer Guillaume).
end


text:
function metGroupCont = createGroupContributionStruct(primaryFile, pH, secondaryFile)
% Generates a matlab structure out of the tab delimited group contribuion data
%
% The matlab structure with the group contibution data for each
% metabolite uses the `primaryFile` file in preference to the `secondaryFile` file but
% these can be any first and second preference files as long as they are
% in the correct format see `webCGMtoTabDelimitedFile.m`
%
% USAGE:
%
%    metGroupCont = createGroupContributionStruct(primaryFile, pH, secondaryFile)
%
% INPUT:
%    primaryFile:     tab delimited text file with group contribution data
%                     (`Janowski et al Biophysical Journal 95:1487-1499 (2008)`)
%                     i.e. output such as `webCGM.txt` from `webCGMtoTabDelimitedFile.m`
%
% OPTIONAL INPUTS:
%    pH:              ph at which group contribution data given for, `default = 7`
%
%    secondaryFile:    tab delimited text file with group contribution data
%                     (`Janowski et al Biophysical Journal 95:1487-1499 (2008)`)
%                     i.e. output such as `webCGM.txt` from `webCGMtoTabDelimitedFile.m`
%                     If the same metabolite abbreviation occurs in both files,
%                     then the data in the primary file takes precedence.
%
%
% Comment on input file format - the first two text columns in both files should correspond to:
% `abbreviation`, `formulaMarvin`,
% the next three columns in both files should correspond to:
% `delta_G_formation`, `delta_G_formation_Uncertainty`, `chargeMarvin`
%
% OUTPUTS:
%    metGroupCont:    structure with fields:
%
%                       * metGroupCont(m).abbreviation - metabolite abbreviation
%                       * metGroupCont(m).formulaMarvin - metabolite formula (Marvin)
%                       * metGroupCont(m).delta_G_formation
%                       * metGroupCont(m).delta_G_formation_uncertainty
%                       * metGroupCont(m).chargeMarvin - metabolite charge (Marvin)
%                       * metGroupCont(m).pH
%                       * metGroupCont(m).file - file data came from
%
% .. Author: - Ronan M. T. Fleming 9 July 2009

if ~exist('pH','var')
    pH=7;
end

fid=fopen(primaryFile,'r');
if fid==-1
    error(['Cannot open ' primaryFile]);
else
    fprintf('%s\n',['Reading primary metabolite group contribution data from: ' primaryFile]);
end
C = textscan(fid,'%s %s %f %f %f');
primaryAbbr=C{1};
C2=C{2};
C3=C{3};
C4=C{4};
C5=C{5};

%
for m=1:length(primaryAbbr)
    metGroupCont(m,1).abbreviation=primaryAbbr{m};
    metGroupCont(m,1).formulaMarvin=C2{m};
    metGroupCont(m,1).delta_G_formation=C3(m,1);
    metGroupCont(m,1).delta_G_formation_uncertainty=C4(m,1);
    metGroupCont(m,1).chargeMarvin=C5(m);
    metGroupCont(m,1).pH=pH;
    metGroupCont(m,1).file=primaryFile;
end

if exist('secondaryFile','var')
    fid=fopen(secondaryFile,'r');
    if fid==-1
        error(['Cannot open ' secondaryFile]);
    else
        fprintf('%s\n',['Reading secondary metabolite group contribution data from: ' secondaryFile]);
    end
    C = textscan(fid,'%s %s %f %f %f');
    secondaryAbbr=C{1};
    C2=C{2};
    C3=C{3};
    C4=C{4};
    C5=C{5};

    m=length(primaryAbbr)+1;
    for n=1:length(secondaryAbbr)
        ind=strcmp(secondaryAbbr{n},primaryAbbr);
        %add the data from the secondary group contribution file if missing from
        %the primary file
        if ~any(ind)
            metGroupCont(m,1).abbreviation=secondaryAbbr{n};
            metGroupCont(m,1).formulaMarvin=C2{n};
            metGroupCont(m,1).delta_G_formation=C3(n);
            metGroupCont(m,1).delta_G_formation_uncertainty=C4(n);
            metGroupCont(m,1).chargeMarvin=C5(n);
            metGroupCont(m,1).pH=pH;
            metGroupCont(m,1).file=secondaryFile;
            m=m+1;
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function str = encodevarname(str,varargin)
%
%    newname = encodevarname(name)
%
%    Encode an invalid variable name using a hex-format for bi-directional
%    conversions.

%    This function is sensitive to the default charset
%    settings in MATLAB, please call feature('DefaultCharacterSet','utf8')
%    to set the encoding to UTF-8 before calling this function.
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        name: a string, can be either a valid or invalid variable name
%
%    output:
%        newname: a valid variable name by converting the leading non-ascii
%              letter into "x0xHH_" and non-ascii letters into "_0xHH_"
%              format, where HH is the ascii (or Unicode) value of the
%              character.
%
%              if the encoded variable name CAN NOT be longer than 63, i.e.
%              the maximum variable name specified by namelengthmax, and
%              one uses the output of this function as a struct or variable
%              name, the name will be trucated at 63. Please consider using
%              the name as a containers.Map key, which does not have such
%              limit.
%
%    example:
%        encodevarname('_a')   % returns x0x5F_a
%        encodevarname('a_')   % returns a_ as it is a valid variable name
%        encodevarname('变量')  % returns 'x0xE58F98__0xE9878F_'
%
%    this file is part of EasyH5 Toolbox: https://github.com/NeuroJSON/easyh5
%
%    License: GPLv3 or 3-clause BSD license, see https://github.com/NeuroJSON/easyh5 for details
%

    if(~isvarname(str(1)))
        if(exist('unicode2native','builtin'))
            str=sprintf('x0x%s_%s',sprintf('%X',unicode2native(str(1))),str(2:end));
        else
            str=sprintf('x0x%X_%s',char(str(1))+0,str(2:end));
        end
    end
    if(isvarname(str))
        return;
    end
    if(exist('unicode2native','builtin'))
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        cpos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(cpos))
            return;
        end
        str0=str;
        pos0=[0 cpos(:)' length(str)];
        str='';
        for i=1:length(cpos)
            str=[str str0(pos0(i)+1:cpos(i)-1) sprintf('_0x%X_',str0(cpos(i))+0)];
        end
        if(cpos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end
end

author: fangq

repo: jsonlab

save_path: github-repos/MATLAB/fangq-jsonlab