text:
function db_reload_subjects( iSubjects )
% DB_RELOAD_SUBJECTS: Reload a set of subjects of the current protocol.
%
% USAGE:  db_reload_subjects( iSubjects );
%
% INPUT:  iSubjects: indices of subjects to modify in protocol subjects array
%                    or 0, to specify the protocol's default subject
%
% DESCRIPTION:
%     For each subject : Parse the subject's directory again and
%     replace Subject structure in Brainstorm database.

% @=============================================================================
% 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, 2008-2010

% Get protocol information
ProtocolInfo     = bst_get('ProtocolInfo');
ProtocolSubjects = bst_get('ProtocolSubjects');
% Check input subject index
if isempty(ProtocolSubjects) || any(iSubjects < 0) || any(iSubjects > length(ProtocolSubjects.Subject))
    error('Invalid subject indices');
end

% If no progressbar is visible: create one
isProgressBar = ~bst_progress('isVisible');
if isProgressBar
    bst_progress('start', 'Reload subject', 'Reloading subjects...', 0, 50 * length(iSubjects));
end
% Loop on all subjects
for i = 1:length(iSubjects)
    iSubject = iSubjects(i);
    % Get subject directory (for default subject or regular subject)
    if (iSubject == 0)
        subjectSubDir = bst_fileparts(ProtocolSubjects.DefaultSubject.FileName);
    else
        subjectSubDir = bst_fileparts(ProtocolSubjects.Subject(iSubject).FileName);
    end
    % Check the existance of the subject's directory
    if ~file_exist(bst_fullfile(ProtocolInfo.SUBJECTS, subjectSubDir))
        db_fix_protocol();
        bst_progress('stop');
        return
    end
    % Parse subject directory
    sSubject = db_parse_subject(ProtocolInfo.SUBJECTS, subjectSubDir, 50);
    % If subject could not be loaded
    if isempty(sSubject)
        db_fix_protocol();
        bst_progress('stop');
        return
    end
    % Else subject was reloaded: update in ProtocolSubjects structure
    if (iSubject == 0)
        ProtocolSubjects.DefaultSubject = sSubject;
    else
        ProtocolSubjects.Subject(iSubject) = sSubject;
    end
end

% Update Brainstorm database
bst_set('ProtocolSubjects', ProtocolSubjects);
% Update display
panel_protocols('UpdateNode', 'Subject', iSubjects);
% Save database
db_save();
% Hide progress bar if it was started here
if isProgressBar
    bst_progress('stop');
end

author: brainstorm-tools


text:
% getSizeFromDimensions.m
% Brenden Kadota, Sunnybrook 2022
%
% Description: Helper function for MRSI processing. Takes an MRSI structure and
% cell array of dimension labels and returns their respective sizes.
% If the dimension does not exist, or is empty, set the dimension size to 1.
%
% Input:
%   MRSIStruct = MRSI structure used in FID-A
%   dimensionLabels = cell array of strings of dimension labels (ie. {'t', 'y', 'x'})
%
% Output:
%   dimensionSizes = dimension sizes from labels. Non-existent dimensions set to 1.

function dimensionSizes = getSizeFromDimensions(MRSIStruct, dimensionLabels)
    arguments
        MRSIStruct (1, 1) struct
        dimensionLabels (1, :) cell
    end
    % get dimension numbers from labels
    dimNumbers = getMultipleDimensions(MRSIStruct, dimensionLabels);
    dataSize = getSize(MRSIStruct);

    % get zero dimensions
    zeroIdx = dimNumbers == 0;
    biggerIndex = dimNumbers > length(dataSize);
    nonExistentIndexes = zeroIdx | biggerIndex;

    % get existing dimensions
    nonZeroIndexes = dimNumbers(~nonExistentIndexes);
    sizeOfNonZeroDimensions = dataSize(nonZeroIndexes);

    % initalize
    dimensionSizes = zeros(1, length(dimensionLabels));
    % set non existent dims to 1
    dimensionSizes(nonExistentIndexes) = 1;
    % set existing dimensions
    dimensionSizes(~nonExistentIndexes) = sizeOfNonZeroDimensions;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [logp, s] = kde(train, test, s)

%KDE Kernel Density Estimation.
%    This function computes a kernel density estimator from a set of examples,
%    by placing Gaussian kernels (with identical masses) on each training data
%    point and adjusting the "widths" to maximize the sum of leave-one-out log
%    densities. In multivariate data, either isotropic or diagonal covariance
%    matrix Gaussians are used.
%
%    usage: [logp, s] = kde(train, test, s)
%
%    inputs:  train (n by D) is a matrix of n training points in D dimensions
%             test  (N by D) is a matrix of test points
%             s     scalar or (1 by D) is a start guess for std devs (optional)
%
%    outputs: logp  (1 by N) is the vector of test log densities
%             s     final std devs
%
%    The size of the initial guess for s indicates whether isotropic or
%    diagonal covariance is desired. If no initial guess is supplied, diagonal
%    is assumed. The method uses a Newton scheme in the log of s, and usually
%    converges in very few iterations. The Newton steps are checked to ensure
%    that a reasonable fraction (half) of the expected improvement is achieved;
%    otherwise smaller steps are tried. The computational complexity is order
%    (nD)^2 + nND, memory requirement order Dn(N+n). The algorithm may
%    encounter numerical problems if initial guess is too small.
%
%    (C) Copyright Carl Edward Rasmussen, July 4th 2000.

[N, D] = size(test); [n, D] = size(train);  % get number of cases and dimension
if nargin == 2                                % if no start guess is given then
  s = std(train)/(n^(0.5/D));      % use scaled empirical axis aligned std devs
end
P = length(s);                                    % number of parameters to fit

t = repmat(train,[1,1,n]);
if P == 1                                    % if we are fitting a single width
  c = sum((permute(t,[1,3,2])-permute(t,[3,1,2])).^2,3);
else                                                 % else multiple parameters
  c = (permute(t,[1,3,2])-permute(t,[3,1,2])).^2;
end

G = 1; TINY = 1e-10;
ec = exp(-sum(c./repmat(permute(2*s.^2,[1,3,2]),[n,n,1]),3));
secd = repmat(sum(ec-eye(n),2),[1,P]);
f_old = sum(log(secd(:,1)))-n*sum(log(s));
while max(abs(G)) > TINY
  x = shiftdim(sum(repmat(ec,[1,1,P]).*c,1))./secd;
  DE = sum(x,1)./s.^2;
  xx = repmat(x,[1 1 P]);
  for i=1:P        % rewriting this loop as matrix expr would cost too much mem
    DDE(i,1:P) = sum(shiftdim(sum(c.*repmat(ec.*c(:,:,i),[1,1,P])))./secd);
  end
  DDE = (DDE - shiftdim(sum(xx.*permute(xx,[1,3,2]))))./(s'*s).^2;
  [v, l] = eig(DDE-2*diag(DE));                               % eigs of Hessian
  l = max(abs(diag(l)),min(l(:)/TINY));    % control sign and magnitude of eigs
  G = (n*D/P-DE)*v*diag(-1./l)*v';                        % compute Newton step
  G = G/max(1, sqrt(G*G'));                       % don't take too large a step
  eta = 1; s_old = s;
  while eta == 1 | f_new < f_old - eta*G*(n*D/P-DE')/2 - TINY    % improvement?
    s = s_old.*exp(G*eta);
    eta = eta/2;           % if we fail, then try smaller step next time around
    ec = exp(-sum(c./repmat(permute(2*s.^2,[1,3,2]),[n,n,1]),3));
    secd = repmat(sum(ec-eye(n),2),[1,P]);
    f_new = sum(log(secd(:,1)))-n*sum(log(s));
  end
  f_old = f_new;                   % remember function value for next iteration
end

x = repmat(2*s.^2,[n,D/P]);
logp = zeros(N, 1);
for i=1:N          % rewriting this loop as matrix expr would cost too much mem
  cc = sum((repmat(test(i,:),[n,1])-train).^2./x,2);
  hh = min(cc);
  logp(i) = log(mean(exp(hh-cc)))-hh;
end

logp = logp - D*log(2*pi)/2 - D*sum(log(s))/P;            % normalize densities

author: zhouyuanzxcv

repo: Hyperspectral


text:
function relaxMtFit(dataDir,outDir,voxRange)
%
% relaxMtFit(dataDir,outDir,voxRange)
%
% Calculates the k and f maps using the MT TOF sequence.
%
% See Yarnykh & Yuan (2004). Cross-relaxation imaging reveals
% detailed anatomy of white matter fiber tracts in the human
% brain. Neuroimage, 23(1):409-24. (PMID: 15325389)
%
% KEY DEPENDENCIES:
% lorentzian.m, relaxMtFitFunc.m
%
% SEE ALSO:
% relaxPreprocess.m to get raw DICOM data into the right format for
% this function.
%
% HISTORY:
% 2006.02.19 Nikola Stikov wrote it.
% 2007.01.?? NS recoded it to use lsqnonlin
% 2007.02.27 RFD rewrote the loop structure and data format.
% 2007.03.02 RFD slight performance enhancements. Also renamed
% several functions and added them all to the mrDiffusion
% repository.
%

if(~exist('dataDir','var')|isempty(dataDir))
  dataDir = pwd;
end

if(~exist('outDir','var')|isempty(outDir))
  outDir = dataDir;
end

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

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

disp(['Loading data from ' dataDir '...']);

ni = niftiRead(fullfile(dataDir,'T1.nii.gz'));
T1 = ni.data;
xform = ni.qto_xyz;
clear ni;
nz = T1>0;

if(exist(fullfile(dataDir,'brainMask.nii.gz'),'file'))
  ni = niftiRead(fullfile(dataDir,'brainMask.nii.gz'));
  brainMask = ni.data==1;
  clear ni;
else
  brainMask = nz;
end

%ni = niftiRead(fullfile(dataDir,'PD.nii.gz'));
%PD = ni.data;

ni = niftiRead(fullfile(dataDir,'S0.nii.gz'));
S0 = ni.data;
clear ni;

d = dir(fullfile(dataDir,'MT_*.nii.gz'));
delta = zeros(1,length(d));
for(ii=1:length(d))
  ni = niftiRead(fullfile(dataDir,d(ii).name));
  MT(:,:,:,ii) = ni.data;
  clear ni;
  % TODO: fix this crude hack. Maybe store offset freqs in nifti header?
  delta(ii) = sscanf(d(ii).name,'MT_%dkHz.nii.gz') ;
end
delta = delta'*1e3; % offset frequencies, in Hz

iterMax = 50; %maximum number of iterations allowed

sz = size(brainMask);

% Clip T1 and PD to reasonable values
% Tissue       T1 (s)          T2 (ms)        PD*
% CSF        0.8 - 20        110 - 2000    70 - 230
% White	    0.76 - 1.08       61 - 100 	   70 - 90
% Gray      1.09 - 2.15       61 - 109     85 - 125
% Meninges   0.5 - 2.2        50 - 165      5 - 44
% Muscle    0.95 - 1.82       20 - 67      45 - 90
% Adipose    0.2 - 0.75       53 - 94      50 - 100
% (PD values are based on PD=111 for 12mM aqueous NiCl2)
% (From http://www.cis.rit.edu/htbooks/mri/chap-8/chap-8.htm#8.7 )
% Our PDs seem to be scaled by 50 or so???
T1(T1<0.01) = 0.01;
T1(T1>10) = 10;
%PD(PD<0) = 0;
%PD(PD>12000) = 12000;

R1 = zeros(size(T1)); R1(nz) = 1./T1(nz);
%if(~exist(fullfile(dataDir,'R1.nii.gz'),'file'))
%  dtiWriteNiftiWrapper(single(R1),xform,fullfile(dataDir,'R1.nii.gz'));
%end

% Tidy-up the brain mask a bit
%brainMask(PD<1000|PD>10000|T1<0.2|T1>2|S0<20) = 0;
brainMask(T1<0.2|T1>2|any(MT==0,4)) = 0;

% RFD: empirically, [1 .08] is the median for brain tissue
x0 = [1 .08];
%x0 = [2.4, .10]'; %this is our initial guess, bese [3.4, .15]'

lb = [.1 .03]; % [k f]
ub = [5 .28];

t_m = 8e-3; %bese 8e-3
t_s = 5e-3; %bese 5e-3
t_r = 19e-3; %bese 19e-3
T2_B = 11e-6;
% Where does this come from? Should we estimate it from the data?
w1rms = 2400; % omega-1 RMS

for ii = 1:length(delta)
    W_B(ii) = pi*(w1rms^2)*lorentzian (delta(ii), T2_B);
end;

options = optimset('LevenbergMarquardt','on', 'Display', 'off');

% FOR DEBUGGING:
%showMontage(brainMask);

brainInds = find(brainMask);
numVoxelsPerUpdate = 10000;
nVoxAll = length(brainInds);
for(ii=1:size(MT,4))
  tmpVol = MT(:,:,:,ii);
  tmpMT(ii,:) = tmpVol(brainInds);
end
clear tmpVol;
MT = tmpMT;
clear tmpMT;
R1 = R1(brainInds);
S0 = S0(brainInds);

f = zeros(1,nVoxAll);
k = zeros(1,nVoxAll);
gof = zeros(1,nVoxAll);
totalSecs = 0;
if(~isempty(voxRange))
  if(any(voxRange<0))
	doVox = [floor(nVoxAll*voxRange(1))+1:ceil(nVoxAll*voxRange(2))];
  else
	doVox = [max(1,voxRange(1)):min(nVoxAll,voxRange(2))];
  end
  fprintf('Processing voxels %d - %d...\n', doVox(1), doVox(end));
else
   doVox = [1:nVoxAll];
   fprintf('Processing all %d voxels...\n',nVoxAll);
end

% What we compute here is actually W_F./R1_F. R1_F is computed in
% the fit function, but we precompute the rest out here to save a
% few cpu cycles in the loop below.
W_F = (w1rms./(2*pi*delta)).^2/.055;

warning off;
nVox = doVox(end)-doVox(1);
tmp = zeros(size(brainMask));
tmpName = tempname
tic;
for(ii=doVox)
  if(mod(ii,numVoxelsPerUpdate)==0)
    prevSecs = toc;
    totalSecs = totalSecs+prevSecs;
    secsPerVox = totalSecs./(ii-doVox(1));
    estTime = secsPerVox*(doVox(end)-ii);
    if(estTime>5400) estTime=estTime./3600; estTimeUnits='hours';
    elseif(estTime>90) estTime=estTime./60; estTimeUnits='minutes';
    else estTimeUnits='seconds'; end
    fprintf('Processed %d of %d voxels- %0.1f %s remaining (%0.3f secs per vox)...\n',ii-doVox(1),nVox,estTime,estTimeUnits,secsPerVox);
	tmp(brainInds) = f;
    m = makeMontage(tmp,[1:5:size(tmp,3)]);
	m = uint8(round(m./max(f).*255));
	imwrite(m,['/home/bob/public_html/f.png']);
	save(tmpName,'k','f','gof','voxRange','xform');
    tic;
  end

  % Some voxels produce a "Input to EIG must not contain NaN
  % or Inf" error in lsqnonl in. Tweaking the bounds or
  % starting estimate can fix it sometimes, but they are
  % probably junk voxels anyway, so we'll catch and skip them.
  try
    %[x, resnorm, residual, exitflag, output] = lsqnonlin(@(x) relaxMtFitFunc(x, MT(:,ii), W_B, W_F, T2_B, R1(ii), S0(ii), t_m, t_s, t_r), x0, lb, ub, options); %bese j-12;
    [x, resnorm, exitflag] = fminsearch(@(x) relaxMtFitFuncLs(x, MT(:,ii), W_B, W_F, T2_B, R1(ii), S0(ii), t_m, t_s, t_r), x0, options);
    if(exitflag>0)
      k(ii) = x(1);
      f(ii) = x(2);
      gof(ii) = resnorm;
    else
      gof(ii) = NaN;
    end
  catch
    % Leave the fit values at zero.
  end
end
warning on;

fprintf('Finished processing %d slices (%d voxels) in %0.1f seconds.\n\n',sz(3),nVox,totalSecs);

if(isempty(voxRange))
   outName = '';
else
   outName = sprintf('%0.2f-%0.2f',voxRange(1),voxRange(2));
end
im=zeros(size(brainMask)); im(brainInds) = f; f = im;
im=zeros(size(brainMask)); im(brainInds) = k; k = im;
im=zeros(size(brainMask)); im(brainInds) = gof; gof = im;
try
   dtiWriteNiftiWrapper(single(f),xform,fullfile(outDir,[outName 'f.nii.gz']));
   dtiWriteNiftiWrapper(single(k),xform,fullfile(outDir,[outName 'k.nii.gz']));
   dtiWriteNiftiWrapper(single(gof),xform,fullfile(outDir,[outName 'gof.nii.gz']));
catch
   outName = tempname
   save(outName,'k','f','gof','voxRange','xform');
end
% FOR DEBUGGING:
keyboard;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [y,zo]=v_filterbank(b,a,x,gd)
%V_FILTERBANK appply v_filterbank to a signal: [y,zo]=(b,a,x,gd)
%
% Inputs:
%    b    numerator coefficients, one row per filter
%    a    denominator coefficients, one row per filter
%    x    input signal
%    gd   group delay of each filter in samples [default=0]. The filter
%         outputs will be advanced to compensate for the group delays.
%         Alternatively, this input can be the zo output from a previous call.
%
% Outputa:
%    y    output signals, one column per filter
%    zo   output filter state

%      Copyright (C) Mike Brookes 2009-2010
%      Version: $Id: v_filterbank.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nf=size(b,1);           % number of filters
nz=max(size(b,2),size(a,2))-1;  % size of state  needed
zzo=zeros(nz,nf);
nx=length(x);           % number of input samples
xx=x(:);
if nargin<4 || ~numel(gd)
    gd=zeros(nf,1);
end
if isstruct(gd)
    zi=gd.zzo;          % get filter state
    qd=gd.qd;
    sd=qd(2)-qd(1);     % number of output samples we need to save
    qd(3)=qd(3)+nx;
    rd=gd.rd;
    yy=zeros(nx+sd,nf);
    yy(1:sd,:)=gd.sy;
    for i=1:nf
        [yy(sd+1:end,i),zzo(:,i)]=filter(b(i,:),a(i,:),xx,zi(:,i));
    end
else
    qd=zeros(1,3);      % [min-delay max-delay #samples]
    rd=round(gd);
    qd(1)=min(rd);
    qd(2)=max(rd);      % find the largest delay
    qd(3)=nx;           % number of filtered samples
    sd=qd(2)-qd(1);     % number of output samples we need to save
    yy=zeros(nx+sd,nf);
    for i=1:nf
        [yy(sd+1:end,i),zzo(:,i)]=filter(b(i,:),a(i,:),xx);
    end
end
ny=max(0,min(nx,qd(3)-qd(2)));    % numer of output samples
y=zeros(ny,nf);
if ny>0
    for i=1:nf
        off=rd(i)-qd(1);
        y(:,i)=yy(off+1:off+ny,i);
    end
end
if nargout>1
    zo.zzo=zzo;  % filter state
    zo.qd=qd;    % offsets
    zo.rd=rd; % rounded group delays
    zo.sy=yy(end-sd+1:end,:);  % save old outputs
end
if ~nargout   % plot pseudo spectrogram
    ng=300;         % target number of columns in image
    kd=max(1,floor(ny/ng));     % decimation factor
    jm=floor(ny/kd);  % % number f frames
yd=reshape(sum(reshape(y(1:kd*jm,:).^2,kd,nf*jm),1),jm,nf)/kd;
    ydm=max(yd(:));
    imagesc((1:jm)*kd+qd(3)-ny-(kd-1)/2,1:nf,10*log10(max(yd,ydm/1e4))');
    axis('xy');
    colorbar;
    v_cblabel('dB');
    xlabel('Sample Number');
    ylabel('Filter Channel');
    title('Filterbank output');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function outputImage = smooth(this, varargin)
% Smoothes image (or image time series) spatially with sinc-truncated Gaussian kernel,
% i.e., using spm_smooth and mimicking its functionality
%
%   Y = MrImageSpm4D()
%   sY = Y.smooth('fwhm', fwhm)
%
% outputImage is a method of class MrImageSpm4D.
%
% NOTE: For complex-valued data, use MrImage.smooth
%
% IN
%   fwhm    [1,1] or [1,3] Full width at half maximum of Gaussian kernel (in mm)
%           if single value is given, isotropic kernel is assumed
%           default: 8
%
% OUT
%
% EXAMPLE
%   smooth
%
%   See also MrImageSpm4D MrImage.smooth spm_smooth

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-02
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% outputImage file is part of the TAPAS UniQC 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/>.

defaults.fwhm = 8;

args = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

if length(fwhm) == 1
    fwhm = fwhm*[ 1 1 1];
end

outputImage = this.copyobj;

% save image file for processing as nii in SPM
outputImage.save('fileName', outputImage.get_filename('prefix', 'raw'));

matlabbatch = outputImage.get_matlabbatch('smooth', fwhm);
save(fullfile(outputImage.parameters.save.path, 'matlabbatch.mat'), ...
            'matlabbatch');
spm_jobman('run', matlabbatch);

% clean up: move/delete processed spm files, load new data into matrix
outputImage.finish_processing_step('smooth');

author: translationalneuromodeling

repo: tapas


text:
function experiments = stack_votlt2018()

set_global_variable('bundle', 'http://data.votchallenge.net/vot2018/longterm/description.json');

baseline.name = 'longterm';
baseline.converter = [];
baseline.type = 'unsupervised';
baseline.tags = {};

baseline.parameters.repetitions = 1;
baseline.parameters.burnin = 0;

baseline.analysis = {'precision_recall', 'speed'};

redetect.name = 'redetection';
redetect.converter = @(sequence) sequence_test_redetection(sequence, ...
    'Length', 200, 'Initialization', 5, 'Padding',2, ...
    'Scaling', 3);

redetect.type = 'unsupervised';
redetect.tags = {};

redetect.parameters.repetitions = 1;
redetect.parameters.burnin = 0;

redetect.analysis = {'redetection'};

experiments = {redetect, baseline};

end

author: votchallenge

repo: toolkit-legacy


text:
classdef MivGui < MimGuiBase
    % MivGui Main class for the MIV 3D Medical Image Viewer for Matlab
    %
    %     To run MIV, type
    %
    %         miv
    %
    %     at the Matlab command prompt
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    methods
        function obj = MivGui(splash_screen)
            obj = obj@MimGuiBase(MivAppDef, splash_screen);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function out = WL_dwtcoeff(y,wname,level)
% WL_dwtcoeff   Discrete wavelet transform coefficients.
%
% Decomposes the time series using a given wavelet and outputs statistics on the
% coefficients obtained up to a maximum level, level.
%
%---INPUTS:
%
% y, the input time series
%
% wname, the mother wavelet, e.g., 'db3', 'sym2' (see Wavelet Toolbox
%           Documentation)
%
% level, the level of wavelet decomposition (can be set to 'max' for the maximum
%               level determined by wmaxlev)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check that a Wavelet Toolbox license is available:
% ------------------------------------------------------------------------------
BF_CheckToolbox('wavelet_toolbox')

% ------------------------------------------------------------------------------
%% Check Inputs
% ------------------------------------------------------------------------------
doPlot = false; % Plot results to figures
N = length(y); % Length of the time series

if nargin < 2 || isempty(wname)
    wname = 'db3'; % Daubechies wavelet filter
end
if nargin < 3 || isempty(level)
    level = 3; % level of wavelet decomposition
end
if strcmp(level,'max')
    level = wmaxlev(N,wname);
end

maxLevelAllowed = wmaxlev(N,wname);
if maxLevelAllowed < level
    fprintf(1,'Chosen level is too large for this wavelet on this signal...\n');
end

% ------------------------------------------------------------------------------
%% Perform Wavelet Decomposition
% ------------------------------------------------------------------------------
% Computes the following:
%   (*) Wavelet decomposition vector, c
%   (*) Bookkeeping vector, l

if maxLevelAllowed < level
    [c, l] = wavedec(y, maxLevelAllowed, wname);
else
    [c, l] = wavedec(y, level, wname);
end

%% Expand DWT coefficients for visualization
% nbcol = 64; % color discretization steps
%
% cfd = zeros(level,N); % detail coefficients
% for k = 1:level
%     d = detcoef(c,l,k);
%     d = d(:)';
%     d = d(ones(1,2^k),:);
%     cfd(k,:) = wkeep1(d(:)',N);
% end
%
% cfd =  cfd(:);
% I = find(abs(cfd)<sqrt(eps));
% cfd(I) = zeros(size(I));
% cfd = reshape(cfd,level,N);
% cfd = wcodemat(cfd,nbcol,'row');

%-------------------------------------------------------------------------------
%% Plotting
%-------------------------------------------------------------------------------
if doPlot
    figure('color','w'); box('on');
    colormap(pink(nbcol));
    image(cfd);
    tics = 1:level;
    labs = int2str((1:level)');
    set(gca,'YTicklabelMode','manual','Ydir','normal', 'Box','On','Ytick',tics,'YTickLabel',labs);
    title('Discrete Wavelet Transform, Absolute Coefficients.');
    xlabel('Time (or Space)')
    ylabel('Level');
end

% ------------------------------------------------------------------------------
%% Get statistics on coefficients
% ------------------------------------------------------------------------------
for k = 1:level
    if k <= maxLevelAllowed
        d = detcoef(c,l,k); % detail coefficients at level k
        % maximum coefficient at this level:
        out.(sprintf('maxd_l%u',k)) = max(d);
        % minimum coefficient at this level:
        out.(sprintf('mind_l%u',k)) = min(d);
        % std coefficients at this level:
        out.(sprintf('stdd_l%u',k)) = std(d);
        % 1-D noise coefficient estimate (estimate of the noise std):
        out.(sprintf('noisestd_l%u',k)) = wnoisest(c,l,k);
    else
        out.(sprintf('maxd_l%u',k)) = NaN;
        out.(sprintf('mind_l%u',k)) = NaN;
        out.(sprintf('stdd_l%u',k)) = NaN;
        out.(sprintf('noisestd_l%u',k)) = NaN;
    end
end

% %% Compress Signal
% % Set approximation coefficients to zero
% % nc = wthcoef('a',c,l);
% NC = wthcoef('t',c,l,N,T,SORH);
%
%
% %% Single Level Reconstruction
% X = waverec(c,l,wname);
% plot(X);
% keyboard

%% Extract Approximation Coefficients from wavelet decomposition structure

% CA = appcoef(C,L,wname,level);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% PARSETXT - parse text input into cell array
%
% Usage: >> cellarray = parsetxt( txt, delims );
%
% Inputs:
%    txt    - input text
%    delims - optional char array of delimiters (default: [' ' ',' 9]);
%
% Note: commas, and simple quotes are ignored
%
% Author: Arnaud Delorme, CNL / Salk Institute, 18 April 2002

% Copyright (C) 18 April 2002 Arnaud Delorme, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function cellarray = parsetxt(txt, delims);

if nargin < 1
	help parsetxt;
	return;
end;

if nargin < 2
    delims = [' ' ',' 9 '"' '''' ];
end

cellarray = {};
tmptxt = '';
for index =1:length(txt)
    if ~isempty(findstr(txt(index), delims))
        if ~isempty(tmptxt), cellarray = { cellarray{:}, tmptxt }; end
        tmptxt = '';
    else
        tmptxt = [ tmptxt txt(index) ];
    end
end;
if ~isempty(tmptxt)
    cellarray = { cellarray{:}, tmptxt };
end
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


disp('Loading iris data');
X = load('data/iris.mtx');
y = load('data/iris.truth');

disp('Running ITML');
num_folds = 2;
knn_neighbor_size = 4;
acc = CrossValidateKNN(y, X, @(y,X) MetricLearningAutotuneKnn(@ItmlAlg, y, X), num_folds, knn_neighbor_size);

disp(sprintf('kNN cross-validated accuracy = %f', acc));


text:
function fr = comp_insdgfb(c,g,shift,Ls,dual)
%COMP_INSDGFB  Non-stationary Gabor filterbank synthesis
%   Usage: fr = comp_insdgfb(c,g,shift,Ls,dual)
%          fr = comp_insdgfb(c,g,shift,Ls)
%          fr = comp_insdgfb(c,g,shift)
%
%   Input parameters:
%         c         : Transform coefficients (matrix or cell array)
%         g         : Cell array of Fourier transforms of the analysis
%                     windows
%         shift     : Vector of frequency shifts
%         Ls        : Original signal length (in samples)
%         dual      : Synthesize with the dual frame
%   Output parameters:
%         fr        : Synthesized signal (Channels are stored in the
%                     columns)
%
%   Given the cell array *c* of non-stationary Gabor coefficients, and a
%   set of filters *g* and frequency shifts *shift* this function computes
%   the corresponding non-stationary Gabor filterbank synthesis.
%
%   If *dual* is set to 1 (default), an attempt is made to compute the
%   canonical dual frame for the system given by *g*, *shift* and the size
%   of the vectors in *c*. This provides perfect reconstruction in the
%   painless case, see the references for more information.
%
%   See also:  cqt, icqt, erblett, ierblett
%
%   References:  ltfatnote018 dogrhove12

% Author: Nicki Holighaus
% Date: 10.04.13

%% Check input arguments
if nargin < 5
    dual = 1;
    if nargin < 3
        error('Not enough input arguments');
    end
end

if iscell(c) == 0 % If matrix format coefficients were used, convert to
    % cell
    [M,N,CH] = size(c);
    c = reshape(c,N*M,CH);
    c = mat2cell(c,M*ones(N,1),CH);
else
    N = length(c);
    CH = size(c{1},2);
    M = cellfun(@(x) size(x,1),c);
end

timepos = cumsum(shift);        % Calculate positions from shift vector
NN = timepos(end);              % Reconstruction length before truncation
timepos = timepos-shift(1);     % Adjust positions

fr = zeros(NN,CH,assert_classname(c{1},g{1})); % Initialize output

if nargin < 4
    Ls = NN; % If original signal length is not given do not truncate
end

if dual == 1 % Attempt to compute canonical dual frame
    g = nsgabdual(g,shift,M,Ls);
end

%% The overlap-add procedure including multiplication with the synthesis
% windows

if numel(M) == 1
    M = M*ones(N,1);
end

for ii = 1:N
    Lg = length(g{ii});

    win_range = mod(timepos(ii)+(-floor(Lg/2):ceil(Lg/2)-1),NN)+1;

    temp = fft(c{ii})*M(ii);
    temp = temp(mod([end-floor(Lg/2)+1:end,1:ceil(Lg/2)]-1,M(ii))+1,:);

    fr(win_range,:) = fr(win_range,:) + ...
        bsxfun(@times,temp,g{ii}([Lg-floor(Lg/2)+1:Lg,1:ceil(Lg/2)]));
end

fr = ifft(fr);
fr = fr(1:Ls,:); % Truncate the signal to original length (if given)

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function intVal=monomialIntPrism(alpha)
%%MONOMIALINTPRISM Evaluate the integral of prod_{i=1}^3 x(i)^alpha(i)
%  over a standard prism in 3D. This is a triangle that has been extruded
%  upward. The vertices are (1,-1,-1), (-1,-1,-1), (-1,1,-1), (1,-1,1),
%  (-1,-1,1), and (-1,1,1).
%
%INPUTS: alpha A 3X1 or 1X3 vector of the integer exponents of the
%              monomial term. All elements must be >=0.
%
%OUTPUTS: intVal The value of the specified integral.
%
%The solution was obtained via analytic integration. These types of
%explicit moment formulae are useful for testing cubature integration
%points.
%
%EXAMPLE:
%An analytic solution for alpha=[10;18;6] is -4/3861. Here, we
%show that this function produces the same result within resonable finite
%precision limits.
% alpha=[19;12;8];
% intSol=monomialIntPrism(alpha);
% trueSol=-4/3861;
% RelErr=(intSol-trueSol)/trueSol
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(mod(alpha(3),2)==0)
    a=alpha(1);
    b=alpha(2);
    c=alpha(3);

    intVal=(2*(-1)^a*(1+a)+2*(-1)^b*(1+b+(-1)^a*(2+a+b)))/((1+a)*(1+b)*(2+a+b)*(1+c));
else
    intVal=0;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Y = spharm_fast(theta,phi,l,m,p0)
%
% Orthonormal spherical harmonic function (Arfken, p. 681) for quantum
% numbers l and m. Angles theta and phi should be given in radians.
% p0 contains the values of the Legendre polynomial
%
function Y = spharm_fast(theta,phi,l,m,p0)

% Copyright (c) 2016, Elekta Oy
% ---------------------------------------
%
% Redistribution and use of the Software in source and binary forms, with or without
% modification, are permitted for non-commercial use.
%
% The Software is provided "as is" without warranties of any kind, either express or
% implied including, without limitation, warranties that the Software is free of defects,
% merchantable, fit for a particular purpose. Developer/user agrees to bear the entire risk
% in connection with its use and distribution of any and all parts of the Software under this license.
%

%p0 = legendre(l,cos(theta));
p = p0(2:end);  % Values m > 1
%
% Arfken, p. 681. The Condon-Shortley phase is included in the
% Legendre polynomials
%
%scale = sqrt((2*l+1)*prod(1:(l-m))/(4*pi*prod(1:(l+m))));
if m < 0
  Y = ((-1)^m)*sqrt((2*l+1)*prod(1:(l+m))/(4*pi*prod(1:(l-m))))*p(-m);
  %Y = scale*((-1)^m)*(prod(1:(l-abs(m)))/prod(1:(l+abs(m))))*p(-m);
elseif m > 0
  Y = sqrt((2*l+1)*prod(1:(l-m))/(4*pi*prod(1:(l+m))))*p(m);
  %Y = scale*p(m);
else
  Y = sqrt((2*l+1)/(4*pi))*p0(1);
  %Y = scale*p0(1);
end
Y = Y*complex(cos(m*phi),sin(m*phi));

author: spm


text:
function [feaNew] = SRKDAtest(fea, model)
% SRKDAtest: Spectral Regression Kernel Discriminant Analysis Testing
%               SRKDAtest uses SRKDA as a dimensionality reduction tool.
%
%       [feaNew,elapse] = SRKDAtest(fea, model);
%
%             Input:
%
%               fea     - data matrix. Each row is a data point.
%             model     - model trained by SRKDAtrain.m
%
%             Output:
%
%             feaNew    - The data in the c-1 SRKDA subspace, where c is the
%                         number of classes.
%
%    Examples:
%
%
% See also SRKDAtrain, SRKDApredict, KSR, KSR_caller
%
%Reference:
%
%   [1] Deng Cai, Xiaofei He, and Jiawei Han. "Speed Up Kernel Discriminant
%   Analysis", The VLDB Journal, vol. 20, no. 1, pp. 21-33, January, 2011.
%
%   [2] Deng Cai, Xiaofei He and Jiawei Han, "SRDA: An Efficient Algorithm for
%   Large Scale Discriminant Analysis" IEEE Transactions on Knowledge and
%   Data Engineering, vol. 20, no. 1, pp. 1-12, January, 2008.
%
%   [3] Deng Cai, "Spectral Regression: A Regression Framework for
%   Efficient Regularized Subspace Learning", PhD Thesis, Department of
%   Computer Science, UIUC, 2009.
%
%   version 3.0 --Jan/2012
%   version 2.0 --December/2011
%   version 1.0 --May/2006
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

MAX_MATRIX_SIZE = 8000; % You can change this number based on your memory.

if ~strcmp(model.TYPE,'SRKDA')
    error('model does not match!');
end

nTrain = size(model.fea,1);
nTest = size(fea,1);
nBlock = ceil(MAX_MATRIX_SIZE*MAX_MATRIX_SIZE/nTrain);
if model.LARs
    feaNew = cell(length(model.LassoCardi),1);
    for i=1:length(model.LassoCardi)
        feaNew{i} = zeros(nTest,size(model.projection{i},2));
    end
    for j = 1:ceil(nTest/nBlock)
        if j == ceil(nTest/nBlock)
            smpIdx = (j-1)*nBlock+1:nTest;
        else
            smpIdx = (j-1)*nBlock+1:j*nBlock;
        end
        KTest= constructKernel(fea(smpIdx,:),model.fea,model.options);
        if model.bSemi
            KTest = KTest*model.KtestHat;
        end
        for i=1:length(model.LassoCardi)
            if model.bSemi
                feaNew{i}(smpIdx,:) = KTest(:,1:model.nLabel)*model.projection{i};
            else
                feaNew{i}(smpIdx,:) = KTest*model.projection{i};
            end
        end
        clear KTest;
    end
else
    feaNew = zeros(nTest,size(model.projection,2));
    for i = 1:ceil(nTest/nBlock)
        if i == ceil(nTest/nBlock)
            smpIdx = (i-1)*nBlock+1:nTest;
        else
            smpIdx = (i-1)*nBlock+1:i*nBlock;
        end
        KTest= constructKernel(fea(smpIdx,:),model.fea,model.options);
        if model.bSemi
            KTest = KTest*model.KtestHat;
            feaNew(smpIdx,:) = KTest(:,1:model.nLabel)*model.projection;
        else
            feaNew(smpIdx,:) = KTest*model.projection;
        end
        clear KTest;
    end
end

author: ZJULearning

repo: MatlabFunc


text:
function s = sudoku_box_invert ( s )

%*****************************************************************************80
%
%% SUDOKU_BOX_INVERT inverts the Sudoku boxes through the center box.
%
%  Discussion:
%
%    In terms of boxes, the transformation takes the Sudoku
%
%      A  B  C
%      D  E  F
%      G  H  I
%
%    to
%
%      I  H  G
%      F  E  D
%      C  B  A
%
%    which still has the "Sudoku" property.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer S(9,9), the Sudoku to be box-inverted.
%
%    Output, integer S(9,9), the box_inverted Sudoku.
%
  t(1:3,1:3) = s(1:3,1:3);
  s(1:3,1:3) = s(7:9,7:9);
  s(7:9,7:9) = t(1:3,1:3);

  t(1:3,1:3) = s(1:3,4:6);
  s(1:3,4:6) = s(7:9,4:6);
  s(7:9,4:6) = t(1:3,1:3);

  t(1:3,1:3) = s(1:3,7:9);
  s(1:3,7:9) = s(7:9,1:3);
  s(7:9,1:3) = t(1:3,1:3);

  t(1:3,1:3) = s(4:6,1:3);
  s(4:6,1:3) = s(4:6,7:9);
  s(4:6,7:9) = t(1:3,1:3);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dtiWarpToWhiteMatter(handles)

view = getSelectedVolume;
anatSize = size(view.anat);
wmVol = zeros(anatSize);

classFileName = [view.leftPath(1:end-5) '.class'];
c = readClassFile(classFileName);
voi = c.header.voi;
wmVol([voi(3):voi(4)], [voi(1):voi(2)], [voi(5):voi(6)]) = permute(c.data, [2,1,3])==c.type.white;

classFileName = [view.rightPath(1:end-5) '.class'];
c = readClassFile(classFileName);
voi = c.header.voi;
wmVolTmp = wmVol;
wmVolTmp([voi(3):voi(4)], [voi(1):voi(2)], [voi(5):voi(6)]) = permute(c.data, [2,1,3])==c.type.white;

% combine left and right
wmVol = wmVol | wmVolTmp;

wmVol = smooth3(wmVol, 'gaussian', 3);

[wmDti,mmDt,xformDt] = dtiGetNamedImage(handles.bg, 'white matter');
xform = handles.xformVAnatToAcpc;
bb = [-72 -110 -45; 72 85 85];
wmVolDti = dtiResliceVolume(wmVol, inv(xform), bb, mmDt);
wmVolDti = permute(wmVolDti, [2 1 3]);

figure;imagesc(wmVolDti(:,:,20)); colorbar;
figure;imagesc(wmDti(:,:,20)); colorbar;
figure;imagesc(wmVolDti(:,:,20)-wmDti(:,:,20)); colorbar;

% M = estMotionIter3(wmDti, wmVolDti, 3, eye(4), 0, 1);
%
% x = inv(dtiGetStandardXform(handles, xformDt)) * M;
% newWmDti = dtiResliceVolume(wmDti, x, bb, mmDt);
% newWmDti = permute(newWmDti, [2 1 3]);
% figure;imagesc(newWmDti(:,:,20)); colorbar;
% figure;imagesc(wmVolDti(:,:,20)-newWmDti(:,:,20)); colorbar;

[deformField, newImg, absDeform, MSE] = dtiDeformationScalar(wmVolDti, wmDti, 25, [7 5 3]);
figure;imagesc(newImg(:,:,20)); colorbar;
figure;imagesc(newImg(:,:,20)-wmVolDti(:,:,20)); colorbar;


text:
function ns3de_test ( )

%*****************************************************************************80
%
%% NS3DE_TEST tests the NS3DE library.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/navier_stokes_3d_exact/ns3de_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( )
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NS3DE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the NS3DE library.\n' );

  ns3de_test01 ( );
  ns3de_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NS3DE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( )

  return
end


text:
% run_hd_specialproc_auto.m
% Jamie Near, Sunnybrook Research Institute 2021.
%
% USAGE:
% [out1,out2,out1_w,out2_w]=run_hd_specialproc_auto(filestring,VoxIDs,aaDomain,tmaxin,iterin);
%
% DESCRIPTION:
% Processing script for Siemens Hadamard-encoded SPECIAL (HD-SPECIAL) MRS
% data in .dat format (twix raw data).  Includes combination of reciever
% channels, removal of bad averages, frequency drift correction,
% leftshifting, and Hadamard Reconstruction.
%
% INPUTS:
% filestring   = String variable for the name of the directory containing
%                   the water suppressed HD-SPECIAL .dat file.  Water unsuppressed
%                   .dat file should be contained in [filestring '_w/'];
% VoxIDs       = (Optional) Two-element cell array with elements that are string variables
%                   giving names to the voxel IDs. Default is {'Vox1','Vox2'}
% aaDomain     = (Optional) Perform the spectral registration (drift correction) using
%                   the full spectrum ('t'), or only a limited frequency
%                   range ('f').  Default is 'f'.
% tmaxin       = (Optional).  Duration (in sec.) of the time domain signal
%                   used in the spectral registration (drift correction).
%                   Default is 0.2 sec.
% iterin       = (Optional).  Maximum number of allowed iterations for the spectral
%                   registration to converge. Default is 20.
%
% OUTPUTS:
% out1         = Fully processed, water suppressed output spectrum from first voxel.
% out2         = Fully processed, water suppressed output spectrum from second voxel.
% out1_w       = Fully processed, water unsuppressed output spectrum from first voxel.
% out2_w       = Fully processed, water unsuppressed output spectrum from second voxel.

function [out1,out2,out1_w,out2_w]=run_hdspecialproc_auto(filestring,VoxIDs,aaDomain,tmaxin,iterin)

if nargin<5
    iterin=20;
    if nargin<4
        tmaxin=0.2;
        if nargin<3
            aaDomain='f';
            if nargin<2
                VoxIDs={'Vox1','Vox2'};
            end
        end
    end
end

%Ensure that filestring is an absolute path, otherwise you can get yourself
%in a tangle later inserting figures at the report stage.
if ~java.io.File(filestring).isAbsolute
    filestring = cd(cd(filestring));
end

%make a new directory for the output report and figures:
reportDir = fullfile(filestring,'report');
mkdir(reportDir);
reportFigsDir = fullfile(filestring,'report','figs');
mkdir(reportFigsDir);

%Find the filename of the SPECIAL dataset
close all
unixString1=fullfile(filestring,'*.dat');
[filename]=dir(unixString1);
filename=filename.name(1:end);

unixStringw=fullfile([filestring '_w'],'*.dat');
[filenamew]=dir(unixStringw);
if ~isempty(filenamew)
    filenamew=filenamew.name(1:end);
    water=true
else
    water=false
    out_w=struct();
    out_w_noproc=struct();
end

%read in the data:
out_raw=io_loadspec_twix(fullfile(filestring,filename));

%load water unsuppressed data and find the coil phases:
if water
    disp('***FOUND WATER UNSUPPRESSED DATA***');
    out_w_raw=io_loadspec_twix(fullfile([filestring '_w/'], filenamew));
    coilcombos1=op_getcoilcombos(op_combinesubspecs(op_fourStepCombine(out_w_raw,1),'diff'),2);
    coilcombos2=op_getcoilcombos(op_combinesubspecs(op_fourStepCombine(out_w_raw,3),'diff'),2);
else
    disp('***WATER UNSUPPRESSED DATA NOT FOUND***');
    coilcombos1=op_getcoilcombos_specReg(op_combinesubspecs(op_fourStepCombine(op_averaging(out_raw),1),'diff'),0,0.01,2);
    coilcombos2=op_getcoilcombos_specReg(op_combinesubspecs(op_fourStepCombine(op_averaging(out_raw),3),'diff'),0,0.01,2);
    out_w='';
    out_w_noproc='';
end

%now combine the coil channels:
[out1_cc,fid1_pre,spec1_pre,ph1,sig1]=op_addrcvrs(out_raw,2,'w',coilcombos1);
[out2_cc,fid2_pre,spec2_pre,ph2,sig2]=op_addrcvrs(out_raw,2,'w',coilcombos2);
if water
    [out1_w_cc,fid1_w_pre,spec1_w_pre,ph1_w,sig1_w]=op_addrcvrs(out_w_raw,2,'w',coilcombos1);
    [out2_w_cc,fid2_w_pre,spec2_w_pre,ph2_w,sig2_w]=op_addrcvrs(out_w_raw,2,'w',coilcombos2);
end

%Do the hadamard reconstruction:
out1_hd=op_fourStepCombine(out1_cc,1);
out2_hd=op_fourStepCombine(out2_cc,3);
if water
    out1_w_hd=op_fourStepCombine(out1_w_cc,1);
    out2_w_hd=op_fourStepCombine(out2_w_cc,3);
end

%make the un-processed spectra, which may be optionally output from this function:
out1_noproc=op_combinesubspecs(op_averaging(out1_hd),'diff');
out2_noproc=op_combinesubspecs(op_averaging(out2_hd),'diff');
if water
    out1_w_noproc=op_combinesubspecs(op_averaging(out1_w_hd),'diff');
    out2_w_noproc=op_combinesubspecs(op_averaging(out2_w_hd),'diff');
end

%plot the data before and after coil phasing:
%figure('position',[0 50 560 420]);
h=figure('visible','off');
subplot(1,3,1);
plot(out_raw.ppm,out_raw.specs(:,:,1,1));xlim([-1 7]);
set(gca,'FontSize',12);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude (a.u.)','FontSize',10);
title('Before correction','FontSize',12);
box off;
subplot(1,3,2);
plot(out_raw.ppm,spec1_pre(:,:,1,1));xlim([-1 7]);
set(gca,'FontSize',12);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude(a.u.)','FontSize',10);
title('After correction','FontSize',12);
box off;
subplot(1,3,3);
plot(out_raw.ppm,spec2_pre(:,:,1,1));xlim([-1 7]);
set(gca,'FontSize',12);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude(a.u.)','FontSize',10);
title('After correction','FontSize',12);
box off;
set(h,'PaperUnits','centimeters');
set(h,'PaperPosition',[0 0 20 10]);
saveas(h,fullfile(reportFigsDir,'coilReconFig'),'jpg');
saveas(h,fullfile(reportFigsDir,'coilReconFig'),'fig');
close(h);

if water
    %figure('position',[0 550 560 420]);
    h=figure('visible','off');
    subplot(1,3,1);
    plot(out_w_raw.ppm,out_w_raw.specs(:,:,1,1));xlim([4 5]);
    set(gca,'FontSize',12);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude (a.u.)','FontSize',10);
    title('Before correction','FontSize',12);
    box off;
    subplot(1,3,2);
    plot(out_w_raw.ppm,spec1_w_pre(:,:,1,1));xlim([4 5]);
    set(gca,'FontSize',12);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude(a.u.)','FontSize',10);
    title('After correction','FontSize',12);
    box off;
    subplot(1,3,3);
    plot(out_w_raw.ppm,spec2_w_pre(:,:,1,1));xlim([4 5]);
    set(gca,'FontSize',12);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude(a.u.)','FontSize',10);
    title('After correction','FontSize',12);
    box off;
    set(h,'PaperUnits','centimeters');
    set(h,'PaperPosition',[0 0 20 10]);
    saveas(h,fullfile(reportFigsDir,'coilReconFig_w'),'jpg');
    saveas(h,fullfile(reportFigsDir,'coilReconFig_w'),'fig');
    close(h);
end

%%%%%%%%%%%%%%%%%%%%%OPTIONAL ALIGNMENT OF SUBSPECTRA%%%%%%%%%%%%%%%%
close all;
fs1_ai=[];
phs1_ai=[];
fs2_ai=[];
phs2_ai=[];
alignISIS='y'
if strcmp(alignISIS,'y') || strcmp(alignISIS,'Y')
    %What we're actually doing is aligning the averages, then aligning the
    %subspectra, then aligning the averages again, and then aligning the
    %subspectra again.
    [out1_ai,fs1_temp,phs1_temp]=op_alignAverages(out1_hd,0.4,'y');
    [out2_ai,fs2_temp,phs2_temp]=op_alignAverages(out2_hd,0.4,'y');
    fs1_ai=fs1_temp;
    fs2_ai=fs2_temp;
    phs1_ai=phs1_temp;
    phs2_ai=phs2_temp;

    [out1_ai,fs1_temp,phs1_temp]=op_alignISIS(out1_ai,0.4);
    [out2_ai,fs2_temp,phs2_temp]=op_alignISIS(out2_ai,0.4);
    fs1_ai(:,2)=fs1_ai(:,2)+fs1_temp;
    fs2_ai(:,2)=fs2_ai(:,2)+fs2_temp;
    phs1_ai(:,2)=phs1_ai(:,2)+phs1_temp;
    phs2_ai(:,2)=phs2_ai(:,2)+phs2_temp;

    [out1_ai,fs1_temp,phs1_temp]=op_alignAverages(out1_ai,0.4,'y');
    [out2_ai,fs2_temp,phs2_temp]=op_alignAverages(out2_ai,0.4,'y');
    fs1_ai=fs1_ai+fs1_temp;
    fs2_ai=fs2_ai+fs2_temp;
    phs1_ai=phs1_ai+phs1_temp;
    phs2_ai=phs2_ai+phs2_temp;

    [out1_ai,fs1_temp,phs1_temp]=op_alignISIS(out1_ai,0.4);
    [out2_ai,fs2_temp,phs2_temp]=op_alignISIS(out2_ai,0.4);
    fs1_ai(:,2)=fs1_ai(:,2)+fs1_temp;
    fs2_ai(:,2)=fs2_ai(:,2)+fs2_temp;
    phs1_ai(:,2)=phs1_ai(:,2)+phs1_temp;
    phs2_ai(:,2)=phs2_ai(:,2)+phs2_temp;

    %for fs_ai and phs_ai, take the average across both subspecs:
    fs1_ai=mean(fs1_ai,2);
    fs2_ai=mean(fs2_ai,2);
    phs1_ai=mean(phs1_ai,2);
    phs2_ai=mean(phs2_ai,2);

    if water
        %Now repeat above for water unsuppressed data:
        [out1_w_ai,fs1_w_temp,phs1_w_temp]=op_alignAverages(out1_w_hd,0.4,'y');
        [out2_w_ai,fs2_w_temp,phs2_w_temp]=op_alignAverages(out2_w_hd,0.4,'y');
        fs1_w_ai=fs1_w_temp;
        fs2_w_ai=fs2_w_temp;
        phs1_w_ai=phs1_w_temp;
        phs2_w_ai=phs2_w_temp;

        [out1_w_ai,fs1_w_temp,phs1_w_temp]=op_alignISIS(out1_w_ai,0.4);
        [out2_w_ai,fs2_w_temp,phs2_w_temp]=op_alignISIS(out2_w_ai,0.4);
        fs1_w_ai(:,2)=fs1_w_ai(:,2)+fs1_w_temp;
        fs2_w_ai(:,2)=fs2_w_ai(:,2)+fs2_w_temp;
        phs1_w_ai(:,2)=phs1_w_ai(:,2)+phs1_w_temp;
        phs2_w_ai(:,2)=phs2_w_ai(:,2)+phs2_w_temp;

        [out1_w_ai,fs1_w_temp,phs1_w_temp]=op_alignAverages(out1_w_ai,0.4,'y');
        [out2_w_ai,fs2_w_temp,phs2_w_temp]=op_alignAverages(out2_w_ai,0.4,'y');
        fs1_w_ai=fs1_w_ai+fs1_w_temp;
        fs2_w_ai=fs2_w_ai+fs2_w_temp;
        phs1_w_ai=phs1_w_ai+phs1_w_temp;
        phs2_w_ai=phs2_w_ai+phs2_w_temp;

        [out1_w_ai,fs1_w_temp,phs1_w_temp]=op_alignISIS(out1_w_ai,0.4);
        [out2_w_ai,fs2_w_temp,phs2_w_temp]=op_alignISIS(out2_w_ai,0.4);
        fs1_w_ai(:,2)=fs1_w_ai(:,2)+fs1_w_temp;
        fs2_w_ai(:,2)=fs2_w_ai(:,2)+fs2_w_temp;
        phs1_w_ai(:,2)=phs1_w_ai(:,2)+phs1_w_temp;
        phs2_w_ai(:,2)=phs2_w_ai(:,2)+phs2_w_temp;

        %for fs_w_ai and phs_w_ai, take the average across both subspecs:
        fs1_w_ai=mean(fs1_w_ai,2);
        fs2_w_ai=mean(fs2_w_ai,2);
        phs1_w_ai=mean(phs1_w_ai,2);
        phs2_w_ai=mean(phs2_w_ai,2);
    end

    %Now check the plots to make sure that they look okay:
    %First make combined subspecs plots:
    out1_hd_temp=op_combinesubspecs(out1_hd,'diff');
    out2_hd_temp=op_combinesubspecs(out2_hd,'diff');
    out1_ai_cc_temp=op_combinesubspecs(out1_ai,'diff');
    out2_ai_cc_temp=op_combinesubspecs(out2_ai,'diff');
    if water
        out1_w_hd_temp=op_combinesubspecs(out1_w_hd,'diff');
        out2_w_hd_temp=op_combinesubspecs(out2_w_hd,'diff');
        out1_w_ai_cc_temp=op_combinesubspecs(out1_w_ai,'diff');
        out2_w_ai_cc_temp=op_combinesubspecs(out2_w_ai,'diff');
    end

    %Now plot them
    close all
    %figure('position',[0 0 560 420]);
    h=figure('visible','off');
    subplot(2,2,1);
    plot(out1_hd_temp.ppm,out1_hd_temp.specs);
    set(gca,'FontSize',12);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude(a.u.)','FontSize',10);
    title('Subspecs not aligned: (all averages)','FontSize',12);
    xlim([0.2 5.2]);
    box off;
    subplot(2,2,3);
    plot(out2_hd_temp.ppm,out2_hd_temp.specs);
    set(gca,'FontSize',12);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude(a.u.)','FontSize',10);
    title('Subspecs not aligned: (all averages)','FontSize',12);
    xlim([0.2 5.2]);
    box off;
    subplot(2,2,2);
    plot(out1_ai_cc_temp.ppm,out1_ai_cc_temp.specs);
    set(gca,'FontSize',8);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude(a.u.)','FontSize',10);
    title('Subspecs aligned: (all averages)','FontSize',12);
    xlim([0.2 5.2]);
    box off;
    subplot(2,2,4);
    plot(out2_ai_cc_temp.ppm,out2_ai_cc_temp.specs);
    set(gca,'FontSize',8);
    set(gca,'XDir','reverse');
    xlabel('Frequency (ppm)','FontSize',10);
    ylabel('Amplitude(a.u.)','FontSize',10);
    title('Subspecs aligned: (all averages)','FontSize',12);
    xlim([0.2 5.2]);
    box off;
    set(h,'PaperUnits','centimeters');
    set(h,'PaperPosition',[0 0 20 15]);
    saveas(h,[filestring '/report/figs/alignSubSpecs_prePostFig'],'jpg');
    saveas(h,[filestring '/report/figs/alignSubSpecs_prePostFig'],'fig');
    close(h);

    if water
        %figure('position',[0 550 560 420]);
        h=figure('visible','off');
        subplot(2,2,1);
        plot(out1_w_hd_temp.ppm,out1_w_hd_temp.specs);
        xlim([3.7 5.7]);
        set(gca,'FontSize',12);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Subspecs not aligned: (all averages)','FontSize',12);
        box off;
        subplot(2,2,3);
        plot(out2_w_hd_temp.ppm,out2_w_hd_temp.specs);
        xlim([3.7 5.7]);
        set(gca,'FontSize',12);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Subspecs not aligned: (all averages)','FontSize',12);
        box off;
        subplot(1,2,2);
        plot(out1_w_ai_cc_temp.ppm,out1_w_ai_cc_temp.specs);
        xlim([3.7 5.7]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Subspecs aligned: (all averages)','FontSize',12);
        box off;
        subplot(2,2,4);
        plot(out2_w_ai_cc_temp.ppm,out2_w_ai_cc_temp.specs);
        xlim([3.7 5.7]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Subspecs aligned: (all averages)','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 15]);
        saveas(h,fullfile(reportFigsDir,'alignSubSpecs_w_prePostFig'),'jpg');
        saveas(h,fullfile(reportFigsDir,'alignSubSpecs_w_prePostFig'),'fig');
        close(h);
    end

%     figure('position',[570 50 560 420]);
%     subplot(2,1,1);
%     plot([1:length(fs_ai)],fs_ai);
%     xlabel('Scan Number');
%     ylabel('Frequency Drift (Hz)');
%     title('Estimated Frequency Drift');
%     subplot(2,1,2);
%     plot([1:length(phs_ai)],phs_ai);
%     xlabel('Scan Number');
%     ylabel('Phase Drift (Degrees)');
%     title('Estimated Phase Drift');

    %sat=input('are you satisfied with alignment of subspecs? ','s');
    sat='y';
    if strcmp(sat,'n') || strcmp(sat,'N')
        out1_ai=out1_hd;
        out2_ai=out2_hd;
        if water
            out1_w_ai=out1_w_hd;
            out2_w_ai=out2_w_hd;
        end
    end

else
    out1_ai=out1_hd;
    out2_ai=out2_hd;
    out1_w_ai=out1_w_hd;
    out2_w_ai=out2_w_hd;
    fs1_ai=zeros(size(out1_hd.fids,out1_hd.dims.averages),1);
    fs2_ai=zeros(size(out2_hd.fids,out2_hd.dims.averages),1);
    phs1_ai=zeros(size(out1_hd.fids,out1_hd.dims.averages),1);
    phs2_ai=zeros(size(out2_hd.fids,out2_hd.dims.averages),1);
end

%Now combine the subspecs
out1_cs=op_combinesubspecs(out1_ai,'diff');
out2_cs=op_combinesubspecs(out2_ai,'diff');
if water
    out1_w_cs=op_combinesubspecs(out1_w_ai,'diff');
    out2_w_cs=op_combinesubspecs(out2_w_ai,'diff');
end

%%%%%%%%%%%%%%%%%%%%%END OF ALIGNMENT OF SUBSPECTRA%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%OPTIONAL REMOVAL OF BAD AVERAGES%%%%%%%%%%%%%%%%
close all;
out1_cs2=out1_cs;
out2_cs2=out2_cs;
nBadAvgTotal1=0;
nBadAvgTotal2=0;
nbadAverages1=1;
nbadAverages2=1;
allAveragesLeft1=[1:out1_cs.sz(out1_cs.dims.averages)]';
allAveragesLeft2=[1:out2_cs.sz(out2_cs.dims.averages)]';
allBadAverages1=[];
allBadAverages2=[];
rmbadav='y';
close all;
if rmbadav=='n' || rmbadav=='N'
    out1_rm=out1_cs;
    out2_rm=out2_cs;
    nsd='N/A';
else
    sat='n'
    while sat=='n'||sat=='N'
        nsd=3; %setting the number of standard deviations;
        iter1=1;
        iter2=1;
        nbadAverages1=1;
        nbadAverages2=1;
        nBadAvgTotal1=0;
        nBadAvgTotal2=0;
        allAveragesLeft1=[1:out1_cs.sz(out1_cs.dims.averages)]';
        allAveragesLeft2=[1:out2_cs.sz(out2_cs.dims.averages)]';
        allBadAverages1=[];
        allBadAverages2=[];
        out1_cs2=out1_cs;
        out2_cs2=out2_cs;
        while nbadAverages1>0;
            [out1_rm,metric1{iter1},badAverages1]=op_rmbadaverages(out1_cs2,nsd,'t');
            allBadAverages1=[allBadAverages1; allAveragesLeft1(badAverages1)];
            badavMask1_temp=zeros(length(allAveragesLeft1),1);
            badavMask1_temp(badAverages1)=1;
            allAveragesLeft1=allAveragesLeft1(~badavMask1_temp);
            nbadAverages1=numel(badAverages1);
            nBadAvgTotal1=nBadAvgTotal1+nbadAverages1;
            out1_cs2=out1_rm;
            iter1=iter1+1;
            disp([num2str(nbadAverages1) ' bad averages removed from spectrum #1 on this iteration.']);
            disp([num2str(nBadAvgTotal1) ' bad averages removed from spectrum #1 in total.']);
            close all;
        end
        while nbadAverages2>0;
            [out2_rm,metric2{iter2},badAverages2]=op_rmbadaverages(out2_cs2,nsd,'t');
            allBadAverages2=[allBadAverages2; allAveragesLeft2(badAverages2)];
            badavMask2_temp=zeros(length(allAveragesLeft2),1);
            badavMask2_temp(badAverages2)=1;
            allAveragesLeft2=allAveragesLeft2(~badavMask2_temp);
            nbadAverages2=numel(badAverages2);
            nBadAvgTotal2=nBadAvgTotal2+nbadAverages2;
            out2_cs2=out2_rm;
            iter2=iter2+1;
            disp([num2str(nbadAverages2) ' bad averages removed from spectrum #2 on this iteration.']);
            disp([num2str(nBadAvgTotal2) ' bad averages removed from spectrum #2 in total.']);
            close all;
        end

        %figure('position',[0 50 560 420]);
        h=figure('visible','off');
        subplot(2,2,1);
        plot(out1_cs.ppm,out1_cs.specs);xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Before removal of bad averages:','FontSize',12);
        box off;
        subplot(2,2,3);
        plot(out2_cs.ppm,out2_cs.specs);xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Before removal of bad averages:','FontSize',12);
        box off;
        subplot(2,2,2);
        plot(out1_rm.ppm,out1_rm.specs);xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('After removal of bad averages:','FontSize',12);
        box off;
        subplot(2,2,4);
        plot(out2_rm.ppm,out2_rm.specs);xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('After removal of bad averages:','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 15]);
        saveas(h,fullfile(reportFigsDir,'rmBadAvg_prePostFig'),'jpg');
        saveas(h,fullfile(reportFigsDir,'rmBadAvg_prePostFig'),'fig');
        close(h);

        %figure('position',[0 550 560 420]);
        h=figure('visible','off');
        plot(1:out1_cs.sz(out1_cs.dims.averages),metric1{1},'.r','MarkerSize',16)
        hold on
        removedAvgs1 = find(~ismember(1:out1_cs.sz(out1_cs.dims.averages),allAveragesLeft1));
        plot(removedAvgs1,metric1{1}(removedAvgs1,1),'ko','MarkerSize',20)
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Deviation Metric','FontSize',10);
        legend('Original averages','Removed Avg')
        legend boxoff;
        title('Deviation Metric','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 10]);
        saveas(h,fullfile(reportFigsDir,'rmBadAvg_vox1_scatterFig1'),'jpg');
        saveas(h,fullfile(reportFigsDir,'rmBadAvg_vox1_scatterFig1'),'fig');
        close(h);

        h=figure('visible','off');
        plot(1:out2_cs.sz(out2_cs.dims.averages),metric2{1},'.r','MarkerSize',16)
        hold on
        removedAvgs2 = find(~ismember(1:out2_cs.sz(out2_cs.dims.averages),allAveragesLeft2));
        plot(removedAvgs2,metric2{1}(removedAvgs2,1),'ko','MarkerSize',20)
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Deviation Metric','FontSize',10);
        legend('Original averages','Removed Avg')
        legend boxoff;
        title('Deviation Metric','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 10]);
        saveas(h,fullfile(reportFigsDir,'rmBadAvg_vox2_scatterFig'),'jpg');
        saveas(h,fullfile(reportFigsDir,'rmBadAvg_vox2_scatterFig'),'fig');
        close(h);

        %sat=input('are you satisfied with bad averages removal? ','s');
        sat='y';
    end
end

close all;

%Now remove the entries from fs_ai and phs_ai that correspond to
%the bad-averages that were removed.
BadAvgMask1=zeros(length(fs1_ai),1);
BadAvgMask2=zeros(length(fs2_ai),1);
BadAvgMask1(allBadAverages1)=1;
BadAvgMask2(allBadAverages2)=1;
size(BadAvgMask1)
size(BadAvgMask2)
size(fs1_ai)
size(fs2_ai)
fs1_ai=fs1_ai(~BadAvgMask1);
fs2_ai=fs2_ai(~BadAvgMask2);
phs1_ai=phs1_ai(~BadAvgMask1);
phs2_ai=phs2_ai(~BadAvgMask2);

%%%%%%%%%%%%%%%%%%%%END OF BAD AVERAGES REMOVAL%%%%%%%%%%%%%%%%%%%%

%now align averages;
%driftCorr=input('Would you like to perform the frequency drift correction?  ','s');
driftCorr='y';
if driftCorr=='n'|| driftCorr=='N'
    out1_aa=out1_rm;
    out2_aa=out2_rm;
    if water
        out1_w_aa=out1_w_cs;
        out2_w_aa=out2_w_cs;
    end
else
    sat='n'
    while sat=='n' || sat=='N'
        out1_rm2=out1_rm;
        out2_rm2=out2_rm;
        fsPoly1=100;
        fsPoly2=100;
        phsPoly1=1000;
        phsPoly2=1000;
        fsCum1=fs1_ai;
        fsCum2=fs2_ai;
        phsCum1=phs1_ai;
        phsCum2=phs2_ai;
        iter1=1;
        iter2=1;
        while (abs(fsPoly1(1))>0.001 || abs(phsPoly1(1))>0.01) && iter1<iterin
            close all
            if aaDomain=='t' || aaDomain=='T'
                tmax=input('input tmax for drift correction: ');
                [out1_aa,fs1,phs1]=op_alignAverages(out1_rm2,tmax,'n');
            elseif aaDomain=='f' || aaDomain=='F'
                tmax=tmaxin+0.04*randn(1);
                fmin=1.8+0.1*randn(1);
                fmaxarray=[2.4,2.85,3.35,4.2,4.4,5.2];
                fmax=fmaxarray(randi(6,1))
                [out1_aa,fs1,phs1]=op_alignAverages_fd(out1_rm2,fmin,fmax,tmax,'n');
            end
            if water
                [out1_w_aa,fs1_w,phs1_w]=op_alignAverages(out1_w_cs,5*tmax,'n');
            end

            fsCum1=fsCum1+fs1;
            phsCum1=phsCum1+phs1;
            fsPoly1=polyfit([1:out1_aa.sz(out1_aa.dims.averages)]',fs1,1)
            phsPoly1=polyfit([1:out1_aa.sz(out1_aa.dims.averages)]',phs1,1)
            iter1
            out1_rm2=out1_aa;
            if water;
                out1_w_cs=out1_w_aa;
            end
            iter1=iter1+1;
        end
        while (abs(fsPoly2(1))>0.001 || abs(phsPoly2(1))>0.01) && iter2<iterin
            close all
            if aaDomain=='t' || aaDomain=='T'
                tmax=input('input tmax for drift correction: ');
                [out2_aa,fs2,phs2]=op_alignAverages(out2_rm2,tmax,'n');
            elseif aaDomain=='f' || aaDomain=='F'
                tmax=tmaxin+0.04*randn(1);
                fmin=1.8+0.1*randn(1);
                fmaxarray=[2.4,2.85,3.35,4.2,4.4,5.2];
                fmax=fmaxarray(randi(6,1))
                [out2_aa,fs2,phs2]=op_alignAverages_fd(out2_rm2,fmin,fmax,tmax,'n');
            end
            if water
                [out2_w_aa,fs2_w,phs2_w]=op_alignAverages(out2_w_cs,5*tmax,'n');
            end

            fsCum2=fsCum2+fs2;
            phsCum2=phsCum2+phs2;
            fsPoly2=polyfit([1:out2_aa.sz(out2_aa.dims.averages)]',fs2,1)
            phsPoly2=polyfit([1:out2_aa.sz(out2_aa.dims.averages)]',phs2,1)
            iter2
            out2_rm2=out2_aa;
            if water;
                out2_w_cs=out2_w_aa;
            end
            iter2=iter2+1;
        end
        %Now plot the cumulative frequency drift correction:
        %figure('position',[0 550 1125 420]);
        h=figure('visible','off');
        subplot(2,2,1);
        plot(out1_rm.ppm,out1_rm.specs);xlim([0 6]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Before drift correction:');
        box off;
        subplot(2,2,2);
        plot(out2_rm.ppm,out2_rm.specs);xlim([0 6]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Before drift correction:');
        box off;
        subplot(2,2,3);
        plot(out1_aa.ppm,out1_aa.specs);xlim([0 6]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('After drift correction:');
        box off;
        subplot(2,2,4);
        plot(out2_aa.ppm,out2_aa.specs);xlim([0 6]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('After drift correction:');
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 15]);
        saveas(h,fullfile(reportFigsDir,'alignAvgs_prePostFig'),'jpg');
        saveas(h,fullfile(reportFigsDir,'alignAvgs_prePostFig'),'fig');
        close(h);

        %figure('position',[0 50 560 420]);
        h=figure('visible','off');
        subplot(2,1,1);
        plot([1:out1_aa.sz(out1_aa.dims.averages)],phsCum1);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Phase Drift (deg.)','FontSize',10);
        title('Estimated Phase Drift','FontSize',12);
        box off;
        subplot(2,1,2);
        plot([1:out2_aa.sz(out2_aa.dims.averages)],phsCum2);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Phase Drift (deg.)','FontSize',10);
        title('Estimated Phase Drift','FontSize',12);
        box off;
        legend('Phase drift','Location','SouthEast');
        legend boxoff;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 10 10]);
        saveas(h,fullfile(reportFigsDir,'phaseDriftFig'),'jpg');
        saveas(h,fullfile(reportFigsDir,'phaseDriftFig'),'fig');
        close(h);

        %figure('position',[570 50 560 420]);
        h=figure('visible','off');
        subplot(2,1,1);
        plot([1:out1_aa.sz(out1_aa.dims.averages)],fsCum1);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Frequency Drift [Hz]','FontSize',10);
        box off;
        subplot(2,1,2);
        plot([1:out2_aa.sz(out2_aa.dims.averages)],fsCum2);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Frequency Drift [Hz]','FontSize',10);
        box off;
        legend('Frequency Drift','Location','SouthEast');
        legend boxoff;
        title('Estimated Frequency Drift','FontSize',12);
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 10 10]);
        saveas(h,fullfile(reportFigsDir,'freqDriftFig'),'jpg');
        saveas(h,fullfile(reportFigsDir,'freqDriftFig'),'fig');
        close(h);

        %sat=input('Are you satisfied with the drift correction? ','s');
        sat='y';
        totalFreqDrift1=mean(max(fsCum1)-min(fsCum1));
        totalFreqDrift2=mean(max(fsCum2)-min(fsCum2));
        totalPhaseDrift1=mean(max(phsCum1)-min(phsCum1));
        totalPhaseDrift2=mean(max(phsCum2)-min(phsCum2));
    end

end

%now do the averaging and left shift to get rid of first order phase:
out1_av=op_leftshift(op_averaging(out1_aa),out1_aa.pointsToLeftshift);
out2_av=op_leftshift(op_averaging(out2_aa),out2_aa.pointsToLeftshift);
if water
    out1_w_av=op_leftshift(op_averaging(out1_w_aa),out1_w_aa.pointsToLeftshift);
    out2_w_av=op_leftshift(op_averaging(out2_w_aa),out2_w_aa.pointsToLeftshift);
end

%Now do an automatic phase correction (Use Creatine peak):
out1_av_zp=op_zeropad(out1_av,16);
out2_av_zp=op_zeropad(out2_av,16);
index1=find(abs(out1_av_zp.specs)==max(abs(out1_av_zp.specs(out1_av_zp.ppm>2.85 & out1_av_zp.ppm<3.15,1))));
index2=find(abs(out2_av_zp.specs)==max(abs(out2_av_zp.specs(out2_av_zp.ppm>2.85 & out2_av_zp.ppm<3.15,1))));
ph01=-phase(out1_av_zp.specs(index1,1))*180/pi;
ph02=-phase(out2_av_zp.specs(index2,1))*180/pi;
out1_ph=op_addphase(out1_av,ph01);
out2_ph=op_addphase(out2_av,ph02);
out1_noproc=op_addphase(op_leftshift(out1_noproc,out1_noproc.pointsToLeftshift),ph01);
out2_noproc=op_addphase(op_leftshift(out2_noproc,out2_noproc.pointsToLeftshift),ph02);
%And now for water unsuppressed data (use water peak):
if water
    out1_w_av_zp=op_zeropad(out1_w_av,16);
    out2_w_av_zp=op_zeropad(out2_w_av,16);
    index1w=find(abs(out1_w_av_zp.specs)==max(abs(out1_w_av_zp.specs(out1_w_av_zp.ppm>4 & out1_w_av_zp.ppm<5.5))));
    index2w=find(abs(out2_w_av_zp.specs)==max(abs(out2_w_av_zp.specs(out2_w_av_zp.ppm>4 & out2_w_av_zp.ppm<5.5))));
    ph01w=-phase(out1_w_av_zp.specs(index1w))*180/pi;
    ph02w=-phase(out2_w_av_zp.specs(index2w))*180/pi;
    out1_w_ph=op_addphase(out1_w_av,ph01w);
    out2_w_ph=op_addphase(out2_w_av,ph02w);
    out1_w_noproc=op_addphase(op_leftshift(out1_w_noproc,out1_w_noproc.pointsToLeftshift),ph01w);
    out2_w_noproc=op_addphase(op_leftshift(out2_w_noproc,out2_w_noproc.pointsToLeftshift),ph02w);
    out1_w_ph_zp=op_addphase(out1_w_av_zp,ph01w);
    out2_w_ph_zp=op_addphase(out2_w_av_zp,ph02w);
end

%Frequency shift all spectra so that Creatine appears at 3.027 ppm:
[~,frqShift1]=op_ppmref(out1_av_zp,2.9,3.1,3.027);
[~,frqShift2]=op_ppmref(out2_av_zp,2.9,3.1,3.027);
out1=op_freqshift(out1_ph,frqShift1);
out2=op_freqshift(out2_ph,frqShift2);
out1_noproc=op_freqshift(out1_noproc,frqShift1);
out2_noproc=op_freqshift(out2_noproc,frqShift2);
%And now for water unsuppressed data (user water peak and set to 4.65 ppm):
if water
    [~,frqShift1w]=op_ppmref(out1_w_ph_zp,4,5.5,4.65);
    [~,frqShift2w]=op_ppmref(out2_w_ph_zp,4,5.5,4.65);
    out1_w=op_freqshift(out1_w_ph,frqShift1w);
    out2_w=op_freqshift(out2_w_ph,frqShift2w);
    out1_w_noproc=op_freqshift(out1_w_noproc,frqShift1w);
    out2_w_noproc=op_freqshift(out2_w_noproc,frqShift2w);
end

%Make figure to show the final spectrum:
h=figure('visible','off');
subplot(1,2,1);
plot(out1.ppm,real(out1.specs),'linewidth',2);xlim([0.2 5.2]);
set(gca,'FontSize',8);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude(a.u.)','FontSize',10);
box off;
subplot(1,2,2);
plot(out2.ppm,real(out2.specs),'linewidth',2);xlim([0.2 5.2]);
set(gca,'FontSize',8);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude(a.u.)','FontSize',10);
box off;
title('Result: Final Spectra','FontSize',12);
set(h,'PaperUnits','centimeters');
set(h,'PaperPosition',[0 0 20 10]);
saveas(h,fullfile(reportFigsDir,'finalSpecFig'),'jpg');
saveas(h,fullfile(reportFigsDir,'finalSpecFig'),'fig');

wrt='y';
if wrt=='y' || wrt=='Y'
    RF1=io_writelcm(out1,fullfile(filestring,[VoxIDs{1} '_lcm']),out1.te);
    RF2=io_writelcm(out2,fullfile(filestring,[VoxIDs{2} '_lcm']),out2.te);
    RF1=io_writelcm(out1_noproc,fullfile(filestring,[VoxIDs{1} '_unprocessed_lcm']),out1_noproc.te);
    RF2=io_writelcm(out2_noproc,fullfile(filestring,[VoxIDs{2} '_unprocessed_lcm']),out2_noproc.te);
    if water
        RF1=io_writelcm(out1_w,fullfile(filestring,[VoxIDs{1} '_w_lcm']),out1_w.te);
        RF2=io_writelcm(out2_w,fullfile(filestring,[VoxIDs{2} '_w_lcm']),out2_w.te);
        RF1=io_writelcm(out1_w_noproc,fullfile(filestring,[VoxIDs{1} '_w_unprocessed_lcm']),out1_w_noproc.te);
        RF2=io_writelcm(out2_w_noproc,fullfile(filestring,[VoxIDs{2} '_w_unprocessed_lcm']),out2_w_noproc.te);
    end
end

close all

%write an html report:
fid=fopen(fullfile(reportDir,'report.html'),'w+');
fprintf(fid,'<!DOCTYPE html>');
fprintf(fid,'\n<html>');
logoPath=which('FID-A_LOGO.jpg');
fprintf(fid,'\n<img src= " %s " width="120" height="120"></body>',logoPath);
fprintf(fid,'\n<h1>FID-A Processing Report</h1>');
fprintf(fid,'\n<h2>Processing pipeline applied to HD-SPECIAL data using run_hdspecialproc_auto.m</h2>');
fprintf(fid,'\n<p>FILENAME: %s/%s/%s </p>',fullfile(filestring,filename));
fprintf(fid,'\n<p>DATE: %s </p>',date);
fprintf(fid,'\n\n<p> </p>');
fprintf(fid,'\n\n<h2>Results of multi-coil combination:</h2>');
fprintf(fid,'\n<img src= " %s " width="800" height="400"></body>',fullfile(reportFigsDir,'coilReconFig.jpg'));
fprintf(fid,'\n\n<p> </p>');
fprintf(fid,'\n\n<h2>Results of alignment of SPECIAL sub-spectra:</h2>');
fprintf(fid,'\n<img src= " %s " width="800" height="400"></body>',fullfile(reportFigsDir,'alignSubSpecs_prePostFig.jpg'));
fprintf(fid,'\n\n<p> </p>');
fprintf(fid,'\n\n<h2>Results of removal of bad averages:</h2>');
fprintf(fid,'\n<p>Original number of averages: \t%5.6f </p>',out_raw.sz(out_raw.dims.averages));
fprintf(fid,'\n<p>Number of bad Averages removed from vox1:  \t%5.6f </p>',nBadAvgTotal1);
fprintf(fid,'\n<p>Number of bad Averages removed from vox2:  \t%5.6f </p>',nBadAvgTotal2);
fprintf(fid,'\n<p>Number of remaining averages in processed dataset (vox1):  \t%5.6f </p>',out1_rm.sz(out1_rm.dims.averages));
fprintf(fid,'\n<p>Number of remaining averages in processed dataset (vox2):  \t%5.6f </p>',out2_rm.sz(out2_rm.dims.averages));
fprintf(fid,'\n<p>Bad Averages Removal Threshold was:  \t%2.2f </p>',nsd);
fprintf(fid,'\n<img src= " %s " width="800" height="600"><img src= " %s " width="800" height="400"><img src= " %s " width="800" height="400">',fullfile(reportFigsDir,'rmBadAvg_prePostFig.jpg'),fullfile(reportFigsDir,'rmBadAvg_vox1_scatterFig.jpg'),fullfile(reportFigsDir,'rmBadAvg_vox2_scatterFig.jpg'));
fprintf(fid,'\n\n<p> </p>');
fprintf(fid,'\n\n<h2>Results of spectral registration:</h2>');
fprintf(fid,'\n<p>Total frequency drift from vox1 was: \t%5.6f </p>',max(totalFreqDrift1));
fprintf(fid,'\n<p>Total frequency drift from vox2 was: \t%5.6f </p>',max(totalFreqDrift2));
fprintf(fid,'\n<p>Total phase drift from vox1 was: \t%5.6f </p>',max(totalPhaseDrift1));
fprintf(fid,'\n<p>Total phase drift from vox2 was: \t%5.6f </p>',max(totalPhaseDrift2));
fprintf(fid,'\n<img src= " %s " width="800" height="600">',fullfile(reportFigsDir,'alignAvgs_prePostFig.jpg'));
fprintf(fid,'\n\n<p> </p>');
fprintf(fid,'\n<img src= " %s " width="400" height="400"><img src=" %s " width="400" height="400">',fullfile(reportFigsDir,'freqDriftFig.jpg'),fullfile(reportFigsDir,'phaseDriftFig.jpg'));
fprintf(fid,'\n\n<p> </p>');
fprintf(fid,'\n\n<h2>Final Result:</h2>');
fprintf(fid,'\n<img src= " %s " width="800" height="400">',fullfile(reportFigsDir,'finalSpecFig.jpg'));
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function area = meshSurfaceArea(vertices, edges, faces)
%MESHSURFACEAREA Surface area of a polyhedral mesh.
%
%   S = meshSurfaceArea(V, F)
%   S = meshSurfaceArea(V, E, F)
%   Computes the surface area of the mesh specified by vertex array V and
%   face array F. Vertex array is a NV-by-3 array of coordinates.
%   Face array can be a NF-by-3 or NF-by-4 numeric array, or a Nf-by-1 cell
%   array, containing vertex indices of each face.
%
%   This functions iterates on faces, extract vertices of the current face,
%   and computes the sum of face areas.
%
%   This function assumes faces are coplanar and convex. If faces are all
%   triangular, the function "trimeshSurfaceArea" should be more efficient.
%
%
%   Example
%     % compute the surface of a unit cube (should be equal to 6)
%     [v f] = createCube;
%     meshSurfaceArea(v, f)
%     ans =
%         6
%
%   See also
%     meshes3d, trimeshSurfaceArea, meshVolume, meshFaceAreas,
%     meshFacePolygons, polygonArea3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-10-13,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% check input number
if nargin == 2
    faces = edges;
end

% pre-compute normals
normals = normalizeVector3d(meshFaceNormals(vertices, faces));

% init accumulator
area = 0;

if isnumeric(faces)
    % iterate on faces in a numeric array
    for i = 1:size(faces, 1)
        poly = vertices(faces(i, :), :);
        area = area + polyArea3d(poly, normals(i,:));
    end

else
    % iterate on faces in a cell array
    for i = 1:length(faces)
        poly = vertices(faces{i}, :);
        area = area + polyArea3d(poly, normals(i,:));
    end
end

function a = polyArea3d(v, normal)

nv = size(v, 1);
v0 = repmat(v(1,:), nv, 1);
products = sum(cross(v-v0, v([2:end 1], :)-v0, 2), 1);
a = abs(dot(products, normal, 2))/2;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                                                  %%%%%
%%%%    IEEE PES Power Grid Library - Optimal Power Flow - v21.07     %%%%%
%%%%          (https://github.com/power-grid-lib/pglib-opf)           %%%%%
%%%%             Benchmark Group - Active Power Increase              %%%%%
%%%%                         29 - July - 2021                         %%%%%
%%%%                                                                  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mpc = pglib_opf_case30_ieee__api
mpc.version = '2';
mpc.baseMVA = 100.0;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	 3	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	2	 2	 36.08	 12.70	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	3	 1	 3.99	 1.20	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	4	 1	 12.64	 1.60	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	5	 2	 156.64	 19.00	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	6	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	7	 1	 37.91	 10.90	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	8	 2	 49.89	 30.00	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	9	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 1.0	 1	    1.06000	    0.94000;
	10	 1	 9.64	 2.00	 0.0	 19.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	11	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 11.0	 1	    1.06000	    0.94000;
	12	 1	 18.62	 7.50	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	13	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 11.0	 1	    1.06000	    0.94000;
	14	 1	 10.31	 1.60	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	15	 1	 13.64	 2.50	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	16	 1	 5.82	 1.80	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	17	 1	 14.97	 5.80	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	18	 1	 5.32	 0.90	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	19	 1	 15.80	 3.40	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	20	 1	 3.66	 0.70	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	21	 1	 29.10	 11.20	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	22	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	23	 1	 5.32	 1.60	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	24	 1	 14.47	 6.70	 0.0	 4.3	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	25	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	26	 1	 5.82	 2.30	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	27	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	28	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 132.0	 1	    1.06000	    0.94000;
	29	 1	 3.99	 0.90	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
	30	 1	 17.63	 1.90	 0.0	 0.0	 1	    1.00000	    0.00000	 33.0	 1	    1.06000	    0.94000;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin
mpc.gen = [
	1	 175.5	 0.0	 176.0	 -176.0	 1.0	 100.0	 1	 351	 0.0; % NG
	2	 171.0	 0.0	 171.0	 -171.0	 1.0	 100.0	 1	 342	 0.0; % NG
	5	 0.0	 0.0	 40.0	 -40.0	 1.0	 100.0	 1	 0	 0.0; % SYNC
	8	 0.0	 0.0	 230.4	 -230.4	 1.0	 100.0	 1	 0	 0.0; % SYNC
	11	 0.0	 3.0	 24.0	 -18.0	 1.0	 100.0	 1	 0	 0.0; % SYNC
	13	 0.0	 9.0	 24.0	 -6.0	 1.0	 100.0	 1	 0	 0.0; % SYNC
];

%% generator cost data
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	 0.0	 0.0	 3	   0.000000	  18.421528	   0.000000; % NG
	2	 0.0	 0.0	 3	   0.000000	  52.182254	   0.000000; % NG
	2	 0.0	 0.0	 3	   0.000000	   0.000000	   0.000000; % SYNC
	2	 0.0	 0.0	 3	   0.000000	   0.000000	   0.000000; % SYNC
	2	 0.0	 0.0	 3	   0.000000	   0.000000	   0.000000; % SYNC
	2	 0.0	 0.0	 3	   0.000000	   0.000000	   0.000000; % SYNC
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	 2	 0.0192	 0.0575	 0.0528	 138.0	 138.0	 138.0	 0.0	 0.0	 1	 -30.0	 30.0;
	1	 3	 0.0452	 0.1652	 0.0408	 152.0	 152.0	 152.0	 0.0	 0.0	 1	 -30.0	 30.0;
	2	 4	 0.057	 0.1737	 0.0368	 139.0	 139.0	 139.0	 0.0	 0.0	 1	 -30.0	 30.0;
	3	 4	 0.0132	 0.0379	 0.0084	 135.0	 135.0	 135.0	 0.0	 0.0	 1	 -30.0	 30.0;
	2	 5	 0.0472	 0.1983	 0.0418	 144.0	 144.0	 144.0	 0.0	 0.0	 1	 -30.0	 30.0;
	2	 6	 0.0581	 0.1763	 0.0374	 139.0	 139.0	 139.0	 0.0	 0.0	 1	 -30.0	 30.0;
	4	 6	 0.0119	 0.0414	 0.009	 148.0	 148.0	 148.0	 0.0	 0.0	 1	 -30.0	 30.0;
	5	 7	 0.046	 0.116	 0.0204	 127.0	 127.0	 127.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 7	 0.0267	 0.082	 0.017	 140.0	 140.0	 140.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 8	 0.012	 0.042	 0.009	 148.0	 148.0	 148.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 9	 0.0	 0.208	 0.0	 142.0	 142.0	 142.0	 0.978	 0.0	 1	 -30.0	 30.0;
	6	 10	 0.0	 0.556	 0.0	 53.0	 53.0	 53.0	 0.969	 0.0	 1	 -30.0	 30.0;
	9	 11	 0.0	 0.208	 0.0	 142.0	 142.0	 142.0	 1.0	 0.0	 1	 -30.0	 30.0;
	9	 10	 0.0	 0.11	 0.0	 267.0	 267.0	 267.0	 1.0	 0.0	 1	 -30.0	 30.0;
	4	 12	 0.0	 0.256	 0.0	 115.0	 115.0	 115.0	 0.932	 0.0	 1	 -30.0	 30.0;
	12	 13	 0.0	 0.14	 0.0	 210.0	 210.0	 210.0	 1.0	 0.0	 1	 -30.0	 30.0;
	12	 14	 0.1231	 0.2559	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	12	 15	 0.0662	 0.1304	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	12	 16	 0.0945	 0.1987	 0.0	 30.0	 30.0	 30.0	 0.0	 0.0	 1	 -30.0	 30.0;
	14	 15	 0.221	 0.1997	 0.0	 20.0	 20.0	 20.0	 0.0	 0.0	 1	 -30.0	 30.0;
	16	 17	 0.0524	 0.1923	 0.0	 38.0	 38.0	 38.0	 0.0	 0.0	 1	 -30.0	 30.0;
	15	 18	 0.1073	 0.2185	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	18	 19	 0.0639	 0.1292	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	19	 20	 0.034	 0.068	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 20	 0.0936	 0.209	 0.0	 30.0	 30.0	 30.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 17	 0.0324	 0.0845	 0.0	 33.0	 33.0	 33.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 21	 0.0348	 0.0749	 0.0	 30.0	 30.0	 30.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 22	 0.0727	 0.1499	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	21	 22	 0.0116	 0.0236	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	15	 23	 0.1	 0.202	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	22	 24	 0.115	 0.179	 0.0	 26.0	 26.0	 26.0	 0.0	 0.0	 1	 -30.0	 30.0;
	23	 24	 0.132	 0.27	 0.0	 29.0	 29.0	 29.0	 0.0	 0.0	 1	 -30.0	 30.0;
	24	 25	 0.1885	 0.3292	 0.0	 27.0	 27.0	 27.0	 0.0	 0.0	 1	 -30.0	 30.0;
	25	 26	 0.2544	 0.38	 0.0	 25.0	 25.0	 25.0	 0.0	 0.0	 1	 -30.0	 30.0;
	25	 27	 0.1093	 0.2087	 0.0	 28.0	 28.0	 28.0	 0.0	 0.0	 1	 -30.0	 30.0;
	28	 27	 0.0	 0.396	 0.0	 75.0	 75.0	 75.0	 0.968	 0.0	 1	 -30.0	 30.0;
	27	 29	 0.2198	 0.4153	 0.0	 28.0	 28.0	 28.0	 0.0	 0.0	 1	 -30.0	 30.0;
	27	 30	 0.3202	 0.6027	 0.0	 28.0	 28.0	 28.0	 0.0	 0.0	 1	 -30.0	 30.0;
	29	 30	 0.2399	 0.4533	 0.0	 28.0	 28.0	 28.0	 0.0	 0.0	 1	 -30.0	 30.0;
	8	 28	 0.0636	 0.2	 0.0428	 140.0	 140.0	 140.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 28	 0.0169	 0.0599	 0.013	 149.0	 149.0	 149.0	 0.0	 0.0	 1	 -30.0	 30.0;
];

% INFO    : === Translation Options ===
% INFO    : Load Model:                  from file ./pglib_opf_case30_ieee.m.api.sol
% INFO    : Gen Active Capacity Model:   stat
% INFO    : Gen Reactive Capacity Model: al50ag
% INFO    : Gen Active Cost Model:       stat
% INFO    :
% INFO    : === Load Replacement Notes ===
% INFO    : Bus 2	: Pd=21.7, Qd=12.7 -> Pd=36.08, Qd=12.70
% INFO    : Bus 3	: Pd=2.4, Qd=1.2 -> Pd=3.99, Qd=1.20
% INFO    : Bus 4	: Pd=7.6, Qd=1.6 -> Pd=12.64, Qd=1.60
% INFO    : Bus 5	: Pd=94.2, Qd=19.0 -> Pd=156.64, Qd=19.00
% INFO    : Bus 7	: Pd=22.8, Qd=10.9 -> Pd=37.91, Qd=10.90
% INFO    : Bus 8	: Pd=30.0, Qd=30.0 -> Pd=49.89, Qd=30.00
% INFO    : Bus 10	: Pd=5.8, Qd=2.0 -> Pd=9.64, Qd=2.00
% INFO    : Bus 12	: Pd=11.2, Qd=7.5 -> Pd=18.62, Qd=7.50
% INFO    : Bus 14	: Pd=6.2, Qd=1.6 -> Pd=10.31, Qd=1.60
% INFO    : Bus 15	: Pd=8.2, Qd=2.5 -> Pd=13.64, Qd=2.50
% INFO    : Bus 16	: Pd=3.5, Qd=1.8 -> Pd=5.82, Qd=1.80
% INFO    : Bus 17	: Pd=9.0, Qd=5.8 -> Pd=14.97, Qd=5.80
% INFO    : Bus 18	: Pd=3.2, Qd=0.9 -> Pd=5.32, Qd=0.90
% INFO    : Bus 19	: Pd=9.5, Qd=3.4 -> Pd=15.80, Qd=3.40
% INFO    : Bus 20	: Pd=2.2, Qd=0.7 -> Pd=3.66, Qd=0.70
% INFO    : Bus 21	: Pd=17.5, Qd=11.2 -> Pd=29.10, Qd=11.20
% INFO    : Bus 23	: Pd=3.2, Qd=1.6 -> Pd=5.32, Qd=1.60
% INFO    : Bus 24	: Pd=8.7, Qd=6.7 -> Pd=14.47, Qd=6.70
% INFO    : Bus 26	: Pd=3.5, Qd=2.3 -> Pd=5.82, Qd=2.30
% INFO    : Bus 29	: Pd=2.4, Qd=0.9 -> Pd=3.99, Qd=0.90
% INFO    : Bus 30	: Pd=10.6, Qd=1.9 -> Pd=17.63, Qd=1.90
% INFO    :
% INFO    : === Generator Setpoint Replacement Notes ===
% INFO    : Gen at bus 1	: Pg=135.5, Qg=5.0 -> Pg=254.0, Qg=49.0
% INFO    : Gen at bus 2	: Pg=46.0, Qg=3.0 -> Pg=261.0, Qg=-46.0
% INFO    : Gen at bus 5	: Pg=0.0, Qg=0.0 -> Pg=0.0, Qg=32.0
% INFO    : Gen at bus 8	: Pg=0.0, Qg=15.0 -> Pg=0.0, Qg=192.0
% INFO    : Gen at bus 11	: Pg=0.0, Qg=9.0 -> Pg=0.0, Qg=15.0
% INFO    : Gen at bus 13	: Pg=0.0, Qg=9.0 -> Pg=0.0, Qg=-1.0
% INFO    :
% INFO    : === Generator Reactive Capacity Atleast Setpoint Value Notes ===
% INFO    : Gen at bus 1	: Qg 49.0, Qmin 0.0, Qmax 10.0 -> Qmin -58.8, Qmax 58.8
% INFO    : Gen at bus 2	: Qg -46.0, Qmin -40.0, Qmax 46.0 -> Qmin -55.2, Qmax 46.0
% INFO    : Gen at bus 8	: Qg 192.0, Qmin -10.0, Qmax 40.0 -> Qmin -230.4, Qmax 230.4
% INFO    : Gen at bus 11	: Qg 15.0, Qmin -6.0, Qmax 24.0 -> Qmin -18.0, Qmax 24.0
% INFO    :
% INFO    : === Generator Classification Notes ===
% INFO    : SYNC   4   -     0.00
% INFO    : NG     2   -   100.00
% INFO    :
% INFO    : === Generator Active Capacity Stat Model Notes ===
% INFO    : Gen at bus 1 - NG	: Pg=254.0, Pmax=271.0 -> Pmax=351   samples: 5
% INFO    : Gen at bus 2 - NG	: Pg=261.0, Pmax=92.0 -> Pmax=342   samples: 12
% INFO    : Gen at bus 5 - SYNC	: Pg=0.0, Pmax=0.0 -> Pmax=0   samples: 0
% INFO    : Gen at bus 8 - SYNC	: Pg=0.0, Pmax=0.0 -> Pmax=0   samples: 0
% INFO    : Gen at bus 11 - SYNC	: Pg=0.0, Pmax=0.0 -> Pmax=0   samples: 0
% INFO    : Gen at bus 13 - SYNC	: Pg=0.0, Pmax=0.0 -> Pmax=0   samples: 0
% INFO    :
% INFO    : === Generator Active Capacity LB Model Notes ===
% INFO    :
% INFO    : === Generator Reactive Capacity Atleast Max 50 Percent Active Model Notes ===
% INFO    : Gen at bus 1 - NG	: Pmax 351.0, Qmin -58.8, Qmax 58.8 -> Qmin -176.0, Qmax 176.0
% INFO    : Gen at bus 2 - NG	: Pmax 342.0, Qmin -55.2, Qmax 46.0 -> Qmin -171.0, Qmax 171.0
% INFO    :
% INFO    : === Generator Setpoint Replacement Notes ===
% INFO    : Gen at bus 1	: Pg=254.0, Qg=49.0 -> Pg=175.5, Qg=0.0
% INFO    : Gen at bus 1	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 2	: Pg=261.0, Qg=-46.0 -> Pg=171.0, Qg=0.0
% INFO    : Gen at bus 2	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 5	: Pg=0.0, Qg=32.0 -> Pg=0.0, Qg=0.0
% INFO    : Gen at bus 5	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 8	: Pg=0.0, Qg=192.0 -> Pg=0.0, Qg=0.0
% INFO    : Gen at bus 8	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 11	: Pg=0.0, Qg=15.0 -> Pg=0.0, Qg=3.0
% INFO    : Gen at bus 11	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 13	: Pg=0.0, Qg=-1.0 -> Pg=0.0, Qg=9.0
% INFO    : Gen at bus 13	: Vg=1.0 -> Vg=1.0
% INFO    :
% INFO    : === Writing Matpower Case File Notes ===

author: power-grid-lib

repo: pglib-opf