text:
function spm_defaults
% Set the defaults which are used by SPM
%__________________________________________________________________________
%
% If you want to customise some defaults for your installation, do not
% modify this file directly, but create a file named spm_my_defaults.m
% instead, accessible from MATLAB search path; e.g., it can be saved in
% MATLAB Startup Folder: userhome/Documents/MATLAB.
%
% Example: create the following file to change the image file extension:
% ----------- file /home/karl/Documents/MATLAB/spm_my_defaults.m -----------
% global defaults
% defaults.images.format = 'img';
%--------------------------------------------------------------------------
%
% spm_defaults should not be called directly in any script or function
% (apart from SPM internals).
% To load the defaults, use spm('Defaults',modality).
% To get/set the defaults, use spm_get_defaults.
%
%                 ** This file should not be edited **
%__________________________________________________________________________
% Copyright (C) 1994-2018 Wellcome Trust Centre for Neuroimaging

% SPM
% $Id: spm_defaults.m 7503 2018-12-03 12:20:12Z guillaume $

global defaults

% Command Line Usage default
%==========================================================================
defaults.cmdline  = 0;

% User Interface defaults
%==========================================================================
defaults.ui.monitor = NaN;
defaults.ui.colour  = [0.58 0.77 0.57];
defaults.ui.fs      = 14;  % unused
defaults.ui.print   = 'ps';
defaults.renderer   = 'opengl';

% File format specific
%==========================================================================
% Note that defaults.analyze.flip is no longer used.  Specifying the
% left/right handedness of the voxel indices is now done entirely by
% spm_flip_analyze_images.m
defaults.images.format  = 'nii'; % options: 'img', 'nii'
defaults.mat.format     = '-v6'; % options: '-mat', '-v6', '-v7.0', '-v7.3'

% Toolboxes defaults
%==========================================================================
defaults.tbx.dir = { fullfile(spm('Dir'),'toolbox') };

% DICOM Import defaults
%==========================================================================
defaults.dicom.root = 'flat'; % Folder hierarchy

% fMRI design defaults
%==========================================================================
defaults.stats.fmri.t   = 16;
defaults.stats.fmri.t0  = 8;
defaults.stats.fmri.hpf = 128;
defaults.stats.fmri.cvi = 'AR(1)';
defaults.stats.fmri.hrf = [6 16 1 1 6 0 32];

% Mask defaults
%==========================================================================
defaults.mask.thresh    = 0.8;

% Stats defaults
%==========================================================================
defaults.stats.maxmem      = 2^30;
defaults.stats.maxres      = 64;
defaults.stats.resmem      = false;
defaults.stats.fmri.ufp    = 0.001;  % Upper tail F-probability
defaults.stats.pet.ufp     = 0.05;
defaults.stats.eeg.ufp     = 0.05;
defaults.stats.topoFDR     = 1;
defaults.stats.rft.nonstat = 0;
defaults.stats.results.volume.distmin =  8;
defaults.stats.results.volume.nbmax   =  3;
defaults.stats.results.svc.distmin    =  4;
defaults.stats.results.svc.nbmax      = 16;
defaults.stats.results.mipmat         = {fullfile(spm('dir'),'MIP.mat')};

% Filename prefix defaults
%==========================================================================
defaults.slicetiming.prefix           = 'a';
defaults.realign.write.prefix         = 'r';
defaults.coreg.write.prefix           = 'r';
defaults.unwarp.write.prefix          = 'u';
defaults.normalise.write.prefix       = 'w';
defaults.deformations.modulate.prefix = 'm';
defaults.smooth.prefix                = 's';
defaults.imcalc.prefix                = 'i';

% Realignment defaults
%==========================================================================
defaults.realign.estimate.quality = 0.9;
defaults.realign.estimate.interp  = 2;
defaults.realign.estimate.wrap    = [0 0 0];
defaults.realign.estimate.sep     = 4;
defaults.realign.estimate.fwhm    = 5;
defaults.realign.estimate.rtm     = 1;
defaults.realign.write.mask       = 1;
defaults.realign.write.interp     = 4;
defaults.realign.write.wrap       = [0 0 0];
defaults.realign.write.which      = [2 1];

% Unwarp defaults
%==========================================================================
defaults.unwarp.estimate.rtm      = 0;
defaults.unwarp.estimate.fwhm     = 4;
defaults.unwarp.estimate.basfcn   = [12 12];
defaults.unwarp.estimate.regorder = 1;
defaults.unwarp.estimate.regwgt   = 1e5;
defaults.unwarp.estimate.foe      = [4 5];
defaults.unwarp.estimate.soe      = [];
defaults.unwarp.estimate.rem      = 1;
defaults.unwarp.estimate.jm       = 0;
defaults.unwarp.estimate.noi      = 5;
defaults.unwarp.estimate.expround = 'Average';
%
% Unwarp uses defaults.realign.write defaults for writing.
%

% Coregistration defaults
%==========================================================================
defaults.coreg.estimate.cost_fun = 'nmi';
defaults.coreg.estimate.sep      = [4 2];
defaults.coreg.estimate.tol      = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
defaults.coreg.estimate.fwhm     = [7 7];
defaults.coreg.write.interp      = 4;
defaults.coreg.write.wrap        = [0 0 0];
defaults.coreg.write.mask        = 0;

% Spatial Normalisation defaults
%==========================================================================
defaults.normalise.write.preserve   = 0;
defaults.normalise.write.bb         = [[-78 -112 -70];[78 76 85]];
defaults.normalise.write.vox        = [2 2 2];
defaults.normalise.write.interp     = 4;

% Old Spatial Normalisation defaults
%==========================================================================
defaults.old.normalise.estimate.smosrc  = 8;
defaults.old.normalise.estimate.smoref  = 0;
defaults.old.normalise.estimate.regtype = 'mni';
defaults.old.normalise.estimate.weight  = '';
defaults.old.normalise.estimate.cutoff  = 25;
defaults.old.normalise.estimate.nits    = 16;
defaults.old.normalise.estimate.reg     = 1;
defaults.old.normalise.write.preserve   = 0;
defaults.old.normalise.write.bb         = [[-78 -112 -70];[78 76 85]];
defaults.old.normalise.write.vox        = [2 2 2];
defaults.old.normalise.write.interp     = 1;
defaults.old.normalise.write.wrap       = [0 0 0];
defaults.old.normalise.write.prefix     = 'w';

% VBM Preprocessing defaults
%==========================================================================
defaults.old.preproc.tpm     = cellstr(char(...
    fullfile(spm('Dir'),'toolbox','OldSeg','grey.nii'),...
    fullfile(spm('Dir'),'toolbox','OldSeg','white.nii'),...
    fullfile(spm('Dir'),'toolbox','OldSeg','csf.nii'))); % Prior probability maps
defaults.old.preproc.ngaus          = [2 2 2 4]'; % Gaussians per class
defaults.old.preproc.warpreg        = 1;          % Warping Regularisation
defaults.old.preproc.warpco         = 25;         % Warp Frequency Cutoff
defaults.old.preproc.biasreg        = 0.0001;     % Bias regularisation
defaults.old.preproc.biasfwhm       = 60;         % Bias FWHM
defaults.old.preproc.regtype        = 'mni';      % Affine Regularisation
defaults.old.preproc.fudge          = 5;          % Fudge factor, can not be set in batch GUI
defaults.old.preproc.samp           = 3;          % Sampling distance
defaults.old.preproc.output.GM      = [0 0 1];
defaults.old.preproc.output.WM      = [0 0 1];
defaults.old.preproc.output.CSF     = [0 0 0];
defaults.old.preproc.output.biascor = 1;
defaults.old.preproc.output.cleanup = 0;

% Smooth defaults
%==========================================================================
defaults.smooth.fwhm = [8 8 8];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Prevent users from making direct calls to spm_defaults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent runOnce
try
    if ~isdeployed && isempty(runOnce)
        d = dbstack;
        if isempty(intersect({'spm','spm_get_defaults'},{d.name}))
            fprintf(['Direct calls to spm_defauts are deprecated.\n' ...
                'Please use spm(''Defaults'',modality) ' ...
                'or spm_get_defaults instead.\n']);
            runOnce = 1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Execute user-specified defaults files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def = defaults;
user_defaults = 'spm_my_defaults.m';
if exist(user_defaults,'file')
    if isdeployed && exist(fullfile(spm('Dir'),user_defaults),'file')
        user_defaults_file = cellstr(fullfile(spm('Dir'),user_defaults));
    else
        user_defaults_file = cellstr(which(user_defaults,'-ALL'));
    end
    for i=1:numel(user_defaults_file)
        try
            spm('run', user_defaults_file{i});
        catch
            lr = lasterror;
            warning(lr.message);
        end
    end
    if spm_check_version('matlab','8.0') >= 0, my_isequaln = @isequaln;
    else my_isequaln = @isequalwithequalnans; end
    if ~my_isequaln(def,defaults)
        fprintf('Defaults settings have been modified by file(s):\n');
        for i=1:numel(user_defaults_file)
            fprintf('  %s\n',user_defaults_file{i});
        end
        fn0 = fieldnames(def);
        mf = fn0(~cellfun(@(x) my_isequaln(def.(x),defaults.(x)),fn0));
        if ~isempty(mf)
            fprintf('Modified fields: ');
            for i=1:numel(mf)
                fprintf('%s ',mf{i});
            end
            fprintf('\n');
        end
    end
end


text:
function forest = forestTrain( data, hs, varargin )
% Train random forest classifier.
%
% Dimensions:
%  M - number trees
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  forest = forestTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .M        - [1] number of trees to train
%   .N1       - [5*N/M] number of data points for training each tree
%   .F1       - [sqrt(F)] number features to sample for each node split
%   .minCount - [1] minimum number of data points to allow split
%   .minChild - [1] minimum number of data points allowed at child nodes
%   .maxDepth - [64] maximum depth of tree
%   .dWts     - [] weights used for sampling and weighing each data point
%   .fWts     - [] weights used for sampling features
%
% OUTPUTS
%  forest   - learned forest model struct array w the following fields
%   .fids     - [Kx1] feature ids for each node
%   .thrs     - [Kx1] threshold corresponding to each fid
%   .child    - [Kx1] index of child for each node
%   .distr    - [KxH] prob distribution at each node
%   .count    - [Kx1] number of data points at each node
%   .depth    - [Kx1] depth of each node
%
% EXAMPLE
%  N=10000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  xs0=single(xs0); xs1=single(xs1);
%  pTrain={'maxDepth',50,'F1',2,'M',150,'minChild',5};
%  tic, forest=forestTrain(xs0,hs0,pTrain{:}); toc
%  hsPr0 = forestApply(xs0,forest);
%  hsPr1 = forestApply(xs1,forest);
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also forestApply, fernsClfTrain
%
% Piotr's Image&Video Toolbox      Version 3.10
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get additional parameters and fill in remaining parameters
dfs={ 'M',1, 'N1',[], 'F1',[], 'minCount',1, 'minChild',1, ...
  'maxDepth',64, 'dWts',[], 'fWts',[] };
[M,N1,F1,minCount,minChild,maxDepth,dWts,fWts]=getPrmDflt(varargin,dfs,1);
[N,F]=size(data); assert(length(hs)==N); assert(all(hs>0)); H=max(hs);
minChild=max(1,minChild); minCount=max([1 minCount minChild]);
if(isempty(N1)), N1=round(5*N/M); end; N1=min(N,N1);
if(isempty(F1)), F1=round(sqrt(F)); end; F1=min(F,F1);
if(isempty(dWts)), dWts=ones(1,N,'single'); end; dWts=dWts/sum(dWts);
if(isempty(fWts)), fWts=ones(1,F,'single'); end; fWts=fWts/sum(fWts);

% make sure data has correct types
if(~isa(data,'single')), data=single(data); end
if(~isa(hs,'uint32')), hs=uint32(hs); end
if(~isa(fWts,'single')), fWts=single(fWts); end
if(~isa(dWts,'single')), dWts=single(dWts); end

% train M random trees on different subsets of data
prmTree = {H,F1,minCount,minChild,maxDepth,fWts};
for i=1:M
  if(N==N1), data1=data; hs1=hs; dWts1=dWts; else
    d=wswor(dWts,N1,4); data1=data(d,:); hs1=hs(d);
    dWts1=dWts(d); dWts1=dWts1/sum(dWts1);
  end
  tree = treeTrain(data1,hs1,dWts1,prmTree);
  if(i==1), forest=tree(ones(M,1)); else forest(i)=tree; end
end

end

function tree = treeTrain( data, hs, dWts, prmTree )
% Train single random tree.
N=size(data,1); K=2*N-1;
[H,F1,minCount,minChild,maxDepth,fWts] = deal(prmTree{:});
thrs=zeros(K,1,'single'); distr=zeros(K,H,'single');
fids=zeros(K,1,'uint32'); child=fids; count=fids; depth=fids;
dids=cell(K,1); dids{1}=uint32(1:N); k=1; K=2;
while( k < K )
  % get node data and store distribution
  dids1=dids{k}; dids{k}=[]; hs1=hs(dids1); count(k)=length(dids1);
  pure=all(hs1(1)==hs1); if(pure), distr(k,hs1(1))=1; else
    distr(k,:)=histc(hs1,1:H)/single(count(k)); end
  % if pure node or insufficient data don't train split
  if( pure||count(k)<=minCount||depth(k)>maxDepth ), k=k+1; continue; end
  % train split and continue
  fids1=wswor(fWts,F1,4); data1=data(dids1,fids1);
  [~,order1]=sort(data1); order1=uint32(order1-1);
  [fid,thr,gini]=forestFindThr(data1,hs1,dWts(dids1),order1,H);
  fid=fids1(fid); left=data(dids1,fid)<thr; count0=nnz(left);
  if( gini<100 && count0>=minChild && (count(k)-count0)>=minChild )
    child(k)=K; fids(k)=fid-1; thrs(k)=thr;
    dids{K}=dids1(left); dids{K+1}=dids1(~left);
    depth(K:K+1)=depth(k)+1; K=K+2;
  end; k=k+1;
end; K=K-1;
% create output model struct
tree=struct('fids',fids(1:K),'thrs',thrs(1:K),'child',child(1:K),...
  'distr',distr(1:K,:),'count',count(1:K),'depth',depth(1:K));
end

function ids = wswor( prob, N, trials )
% Fast weighted sample without replacement. Alternative to:
%  ids=datasample(1:length(prob),N,'weights',prob,'replace',false);
M=length(prob); assert(N<=M); if(N==M), ids=1:N; return; end
if(all(prob(1)==prob)), ids=randperm(M,N); return; end
cumprob=min([0 cumsum(prob)],1); assert(abs(cumprob(end)-1)<.01);
cumprob(end)=1; [~,ids]=histc(rand(N*trials,1),cumprob);
[s,ord]=sort(ids); K(ord)=[1; diff(s)]~=0; ids=ids(K);
if(length(ids)<N), ids=wswor(cumprob,N,trials*2); end
ids=ids(1:N)';
end

author: zouchuhang

repo: LayoutNet


text:
function betweenScanMotEst(view, baseScan, targetScans)
%
% betweenScanMotEst(view, baseScan, targetScans)
%
% Robust 3D rigid body motion estimation between MEAN MAPS of different scans
%
% 04/05 Ress, rewritten from betweenScanMotComp.m

if ~exist('baseScan', 'var'), baseScan = selectScans(view, 'Select base scan'); end
if isempty(baseScan), return, end
baseScan = baseScan(1);
if ~exist('targetScans', 'var'), targetScans = selectScans(view, 'Select target scans'); end
% removes the base scan, if present
targetScans = targetScans(find(targetScans~=baseScan));
if isempty(targetScans), return, end

% Get or compute Mean Maps.
view = loadMeanMap(view);
meanMap = view.map;

% if the number of slices is too small, repeat the first and last slice
% to avoid running out of data (the derivative computation discards the
% borders in z, tipically 2 slices at the begining and 2 more at the end)
if size(meanMap,3)<=8
  meanMap = cat(3, meanMap(:,:,1,:), meanMap(:,:,1,:), meanMap,...
    meanMap(:,:,end,:), meanMap(:,:,end,:));
end

% get base mean map
baseMeanMap = meanMap{baseScan};

% Do motion estimation for each scan.
corrected = 0 * targetScans;
for iScan=1:length(targetScans)
  scan = targetScans(iScan);
  nSlices = length(sliceList(view,scan));
  dims = sliceDims(view,scan);

  % estimate motion between mean maps
  M = estMotionIter3(baseMeanMap,meanMap{scan},3,eye(4),1,1); % rigid body, ROBUST
  midX = [dims/2 nSlices/2]';
  midXp = M(1:3, 1:3) * midX; % Rotational motion
  rotMot = sqrt(sum((midXp - midX).^2));
  transMot = sqrt(sum(M(1:3, 4).^2)); % Translational motion
  totalMot = sqrt(rotMot^2 + transMot^2);
  disp(['Scan ', int2str(scan), ' - motion (voxels): rot = ', num2str(rotMot), '; trans = ', num2str(transMot), ...
      ' total = ', num2str(totalMot)])

end % scan LOOP

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Multichannel version of im2col. Rearranges image blocks into columns.
% The vectorized blocks of each channel are staked vertically in the output.
%
% USAGE: patches = im2col_ch(im, [ph pw], mode)
%
%  -> im      : input image
%  -> ph,pw   : patch size (ph x pw)
%  -> mode    : either 'distinct' of 'sliding'
%
%  <- patches : output patches (ph pw ch x n)
function patches = im2col_ch(im, psz, mode)

	if nargin < 3,
		mode = 'sliding';
	end

	pdim = prod(psz);
	ch = size(im,3);

	% extract patches from first channel
	patches = im2col(im(:,:,1),psz,mode);

	% number of patches
	n = size(patches,2);

	% allocate room for the other channels
	patches = [patches ; zeros((ch-1)*size(patches,1), size(patches,2))];
	for c = 2:ch,
		patches((c-1)*pdim + [1:pdim],:) = im2col(im(:,:,c), psz, mode);
	end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function sim = ImageSimilarity (im1, im2)
% This function computes the "similarity score" between two images. You
% should use the value for the similarity factor value when the two images
% are assigned the same character.
%
% Input:
%   im1, im2: Two images from the provided dataset (they should be 16x8
%     matrices of 0s and 1s).
%
% Output:
%   sim: The similarity score of those images.
%
% Copyright (C) Daphne Koller, Stanford University, 2012

a = im1(:);
b = im2(:);

meanSim = 0.283; % Avg sim score computed over held-out data.

cosDist = (a' * b) / (norm(a) * norm(b));

diff = (cosDist - meanSim) ^ 2;

if (cosDist > meanSim)
    sim = 1 + 5*diff;
else
    sim = 1 / (1 + 5*diff);
end

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function GUI = sessionGUI_removeStudy(study, forceRemove);
% Remove a study from the mrVista GUI.
%
% <GUI> = sessionGUI_removeStudy(<study=current>, <forceRemove=0>);
%
% This function will update or create the saved studies file
% (saved in the mrVista 2 repository: mrVista2/study/studies.mat).
%
% If the forceRemove flag is set to 1, will remove without prompting.
% Otherwise, asks the user to confirm first.
%
% ras, 07/06.
mrGlobals2;

if notDefined('study'),    study = GUI.settings.study;      end
if notDefined('forceRemove'), forceRemove = 0;              end

% N will index the index for this study
if isnumeric(study)
    N = study;
    study = GUI.studies(N).name;

elseif ischar(study)
    N = cellfind({GUI.studies.name}, study);

else
    error('Invalid study specification.')

end

% can't remove the first study (recent sessions); check that we're not
% trying:
if N==1
    myWarnDlg('Sorry, you can not remove the Recent Sessions study.');
    return
end

% Ask the user if he/she is sure, if not forcing the removal:
if ~forceRemove
    q = sprintf('Permanently delete study %s? ', study);
    resp = questdlg(q, mfilename);
    if ~isequal(resp, 'Yes'), return; end
end

% remove the study
ok = setdiff(1:length(GUI.studies), N);
GUI.studies = GUI.studies(ok);

% save the updated studies file
studyPath = studySave(GUI.studies);
fprintf('Studies file %s updated.\n', studyPath);

% select the previous study in the study list
GUI = sessionGUI_selectStudy(N-1);

return


text:
function s=get_name(a)
s=[get_name(a.algorithm)];
eval_name

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%CODEGENERATOR.GENMFUNFRICTION Generate M-function for joint friction
%
% cGen.genmfunfriction() generates a robot-specific M-function to compute
% joint friction.
%
% Notes::
% - Is called only if cGen has active flag genmfun
% - Access to generated function is provided via subclass of SerialLink
%   whose class definition is stored in cGen.robjpath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.gengravload.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [ F ] = genmfunfriction( CGen )

%% Does robot class exist?
if ~exist(fullfile(CGen.robjpath,[CGen.getrobfname,'.m']),'file')
    CGen.logmsg([datestr(now),'\tCreating ',CGen.getrobfname,' m-constructor ']);
    CGen.createmconstructor;
    CGen.logmsg('\t%s\n',' done!');
end

%% Generate m-function
CGen.logmsg([datestr(now),'\tGenerating joint friction m-function']);
symname = 'friction';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if exist(fname,'file')
    tmpStruct = load(fname);
else
    error ('genmfunfriction:SymbolicsNotFound','Save symbolic expressions to disk first!')
end

funfilename = fullfile(CGen.robjpath,[symname,'.m']);
[~,qd] = CGen.rob.gencoords;

matlabFunction(tmpStruct.(symname),'file',funfilename,...              % generate function m-file
    'outputs', {'F'},...
    'vars', {'rob',[qd]});
hStruct = createHeaderStruct(CGen.rob,symname);                 % replace autogenerated function header
replaceheader(CGen,hStruct,funfilename);
CGen.logmsg('\t%s\n',' done!');

end

%% Definition of the header contents for each generated file
function hStruct = createHeaderStruct(rob,fname)
[~,hStruct.funName] = fileparts(fname);
hStruct.shortDescription = ['Joint friction for the ',rob.name,' arm.'];
hStruct.calls = {['F = ',hStruct.funName,'(rob,qd)'],...
    ['F = rob.',hStruct.funName,'(qd)']};
hStruct.detailedDescription = {['Given a full set of generalized joint velocities the function'],...
    'computes the friction forces/torques.'};
hStruct.inputs = { ['rob: robot object of ', rob.name, ' specific class'],...
                   ['qd:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     velocities', ...
                   'Angles have to be given in radians!'};
hStruct.outputs = {['F:  [',int2str(rob.n),'x1] vector of joint friction forces/torques']};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
     'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {'gravload'};
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function setGuidataField(fig,fieldName,x)
% Sets a single field from fig's guidata(), assumed to be a scalar struct,
% to the value x.

g=guidata(fig);
g.(fieldName)=x;
guidata(fig,g);

end

author: kristinbranson

repo: JAABA


text:
function [loc_assort_pos,loc_assort_neg] = local_assortativity_wu_sign(W)
%LOCAL_ASSORTATIVITY_WU_SIGN     Local Assortativity
%
%   [loc_assort_pos,loc_assort_neg] = local_assortativity_wu_sign(W);
%
%   Local Assortativity measures the extent to which nodes are connected to
%   nodes of similar strength (vs. higher or lower strength). Adapted from
%   Thedchanamoorthy et al. (2014)'s formula to allow weighted/signed
%   networks (node degree replaced with node strength). Note, output values
%   sum to total assortativity.
%
%   Inputs:     W,        undirected connection matrix with positive and
%                         negative weights
%
%   Output:     loc_assort_pos, local assortativity from positive weights
%
%               loc_assort_neg, local assortativity from negative weights
%
%   Reference: Thedchanamoorthy G, Piraveenan M, Kasthuriratna D,
%              Senanayake U. Proc Comp Sci (2014) 29:2449-2461.
%
%
%   Jeff Spielberg, Boston University

%   Modification History:
%   May 2015: Original

W(1:(size(W,1)+1):end) = 0;
r_pos = assortativity_wei(W.*(W>0),0);
r_neg = assortativity_wei(-W.*(W<0),0);
[str_pos,str_neg] = strengths_und_sign(W);
loc_assort_pos = nan(size(W,1),1);
loc_assort_neg = nan(size(W,1),1);

for curr_node = 1:size(W,1)
    [~,j_pos] = find(W(curr_node,:)>0);
    loc_assort_pos(curr_node,1) = sum(abs(str_pos(j_pos)-str_pos(curr_node)))/str_pos(curr_node);

    [~,j_neg] = find(W(curr_node,:)<0);
    loc_assort_neg(curr_node,1) = sum(abs(str_neg(j_neg)-str_neg(curr_node)))/str_neg(curr_node);
end

loc_assort_pos = ((r_pos+1)/size(W,1))-(loc_assort_pos/sum(loc_assort_pos));
loc_assort_neg = ((r_neg+1)/size(W,1))-(loc_assort_neg/sum(loc_assort_neg));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ fea, out ] = ex_diffusion2( varargin )
%EX_DIFFUSION2 1D Time dependent diffusion equation examples.
%
%   [ FEA, OUT ] = EX_DIFFUSION2( VARARGIN ) 1D time dependent diffusion equation on
%   a line with exact solutions. Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       nu          scalar {1e-1}          Diffusion coefficient
%       icase       scalar {2}             Test case
%       hmax        scalar {1/25}          Max grid cell size
%       dt          scalar {0.1}           Time step size
%       ischeme     scalar {3}             Time stepping scheme
%       sfun        string {sflag1}        Shape function
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'nu',       1e-1; ...
  'icase',    2; ...
  'hmax',     1/25; ...
  'dt'        0.1; ...
  'ischeme'   3; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'tol',      1e-2; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

tmax = 1;
switch( opt.icase )
  case 1
    refsol = '(1-x)*2.2+x*3.3';
  case 2
    n = 1;
    refsol = ['exp(-',num2str(opt.nu*n^2*pi^2),'*t)*sin(',num2str(n*pi),'*x)'];
end

% Grid generation.
fea.grid = linegrid( 1/opt.hmax, 0, 1 );

% Problem definition.
fea.sdim  = { 'x' };
fea = addphys( fea, @convectiondiffusion );
fea.phys.cd.sfun = { opt.sfun };
fea.phys.cd.eqn.coef{2,4} = { opt.nu };
fea = parsephys(fea);

% Parse and solve problem.
x = fea.grid.p';
n = length(x);
if( strcmp( opt.sfun,'sflag2' ) )
  x = [ x; (x(2:end)+x(1:end-1))/2 ];
end
t   = 0;
u0  = eval( refsol );
fea = parseprob( fea );
if( opt.ischeme>0 )

  fea.bdr.d{1} = refsol;
  fea.bdr.d{2} = refsol;
  fea.bdr.n    = cell(1,2);
  [fea.sol.u,tlist] = solvetime( fea, 'fid', fid, 'init', u0, 'ischeme', opt.ischeme, 'tstep', opt.dt, 'tmax', tmax, 'nstbwe', 0, 'tstop', 0, 'imass', 2 );

else

  [M,A,f] = assembleprob( fea, 'f_m', 1, 'imass', 1, 'f_a', 1, 'f_f', 1, 'f_sparse', 1 );
  M = spdiags( full(sum(M')'), 0, size(M,1), size(M,1) );
  fea.sol.u = u0;
  dt = opt.dt;
  it = 0;
  tlist = 0;
  if( opt.ischeme==-1 )       % Crank-Nicolson.
    C = [ M + dt/2*A ];
  elseif( opt.ischeme==-2 )   % Backward Euler.
    C = [ M + dt*A ];
  end
  C(1,:) = 0;
  C(1,1) = 1;
  C(n,:) = 0;
  C(n,n) = 1;
  while 1
    t  = t  + dt;
    it = it + 1;

    u_r = eval( refsol );
    if( opt.ischeme==-1 )       % Crank-Nicolson.
      b = [ [ M - dt/2*A ]*u0 + dt*f ];
    elseif( opt.ischeme==-2 )   % Backward Euler.
      b = [ M*u0 + dt*f ];
    end
    b(1) = u_r(1);
    b(n) = u_r(n);

    u1 = C\b;

    if( t>=tmax )
      break
    end

    err = norm( u1 - u_r )/norm( u_r );
    errnm(it) = err;
    if( ~isempty(fid) )
      fprintf( fid, 'Time = %f, error norm = %d\n', t, err );
    end

    fea.sol.u = [ fea.sol.u u1 ];
    tlist = [ tlist t ];
    u0 = u1;
  end

end

% Postprocessing.
if( opt.iplot>0 )
  figure;
  if( opt.iplot>1 )
    i_sol_list = 1:numel(tlist);
  else
    i_sol_list = numel(tlist);
  end
  [~,ix] = sort( x );
  for i_sol=i_sol_list
    t = tlist(i_sol);
    clf
    postplot( fea, 'surfexpr', 'c', 'solnum', i_sol );
    hold on
    u_r = eval( refsol );
    plot( sort(x), u_r(ix), 'r--' );
    title( ['Solution at time ',num2str(t)])
    xlabel( 'x' )
    drawnow
  end
end

% Error checking.
for i_sol=1:numel(tlist)
  u_i = fea.sol.u(:,i_sol);
  t   = tlist(i_sol);
  u_r = eval( refsol );
  errnm(i_sol) = norm( u_i - u_r )/norm( u_r );
end
out.err  = errnm;
out.pass = all( errnm<opt.tol );

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function preproc = spm_cfg_preproc
% SPM Configuration file for toolbox 'Old Segment'
%______________________________________________________________________
% Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_preproc.m 6942 2016-11-21 13:17:44Z guillaume $

if ~isdeployed, addpath(fullfile(spm('dir'),'toolbox','OldSeg')); end

% ---------------------------------------------------------------------
% data Data
% ---------------------------------------------------------------------
data         = cfg_files;
data.tag     = 'data';
data.name    = 'Data';
data.help    = {'Select scans for processing. This assumes that there is one scan for each subject. Note that multi-spectral (when there are two or more registered images of different contrasts) processing is not yet implemented for this method.'};
data.filter = 'image';
data.ufilter = '.*';
data.num     = [1 Inf];
% ---------------------------------------------------------------------
% GM Grey Matter
% ---------------------------------------------------------------------
GM         = cfg_menu;
GM.tag     = 'GM';
GM.name    = 'Grey Matter';
GM.help    = {'Options to produce grey matter images: c1*, wc1* and mwc1*.'};
GM.labels = {
             'None'
             'Native Space'
             'Unmodulated Normalised'
             'Modulated Normalised'
             'Native + Unmodulated Normalised'
             'Native + Modulated Normalised'
             'Native + Modulated + Unmodulated'
             'Modulated + Unmodulated Normalised'
}';
GM.values = {[0 0 0] [0 0 1] [0 1 0] [1 0 0] [0 1 1] [1 0 1] [1 1 1]...
             [1 1 0]};
GM.def    = @(val)spm_get_defaults('old.preproc.output.GM', val{:});
% ---------------------------------------------------------------------
% WM White Matter
% ---------------------------------------------------------------------
WM         = cfg_menu;
WM.tag     = 'WM';
WM.name    = 'White Matter';
WM.help    = {'Options to produce white matter images: c2*, wc2* and mwc2*.'};
WM.labels = {
             'None'
             'Native Space'
             'Unmodulated Normalised'
             'Modulated Normalised'
             'Native + Unmodulated Normalised'
             'Native + Modulated Normalised'
             'Native + Modulated + Unmodulated'
             'Modulated + Unmodulated Normalised'
}';
WM.values = {[0 0 0] [0 0 1] [0 1 0] [1 0 0] [0 1 1] [1 0 1] [1 1 1]...
             [1 1 0]};
WM.def    = @(val)spm_get_defaults('old.preproc.output.WM', val{:});
% ---------------------------------------------------------------------
% CSF Cerebro-Spinal Fluid
% ---------------------------------------------------------------------
CSF         = cfg_menu;
CSF.tag     = 'CSF';
CSF.name    = 'Cerebro-Spinal Fluid';
CSF.help    = {'Options to produce CSF images: c3*, wc3* and mwc3*.'};
CSF.labels = {
              'None'
              'Native Space'
              'Unmodulated Normalised'
              'Modulated Normalised'
              'Native + Unmodulated Normalised'
              'Native + Modulated Normalised'
              'Native + Modulated + Unmodulated'
              'Modulated + Unmodulated Normalised'
}';
CSF.values = {[0 0 0] [0 0 1] [0 1 0] [1 0 0] [0 1 1] [1 0 1] [1 1 1]...
             [1 1 0]};
CSF.def    = @(val)spm_get_defaults('old.preproc.output.CSF', val{:});
% ---------------------------------------------------------------------
% biascor Bias Corrected
% ---------------------------------------------------------------------
biascor         = cfg_menu;
biascor.tag     = 'biascor';
biascor.name    = 'Bias Corrected';
biascor.help    = {'This is the option to produce a bias corrected version of your image. MR images are usually corrupted by a smooth, spatially varying artifact that modulates the intensity of the image (bias). These artifacts, although not usually a problem for visual inspection, can impede automated processing of the images.  The bias corrected version should have more uniform intensities within the different types of tissues.'};
biascor.labels = {
                  'Save Bias Corrected'
                  'Don''t Save Corrected'
}';
biascor.values = {1 0};
biascor.def    = @(val)spm_get_defaults('old.preproc.output.biascor', val{:});
% ---------------------------------------------------------------------
% cleanup Clean up any partitions
% ---------------------------------------------------------------------
cleanup         = cfg_menu;
cleanup.tag     = 'cleanup';
cleanup.name    = 'Clean up any partitions';
cleanup.help    = {
                   'This uses a crude routine for extracting the brain from segmented images.  It begins by taking the white matter, and eroding it acouple of times to get rid of any odd voxels.  The algorithm continues on to do conditional dilations for several iterations, where the condition is based upon gray or white matter being present. This identified region is then used to clean up the grey and white matter partitions, and has a slight influences on the CSF partition.'
                   ''
                   'If you find pieces of brain being chopped out in your data, then you may wish to disable or tone down the cleanup procedure.'
}';
cleanup.labels = {
                  'Dont do cleanup'
                  'Light Clean'
                  'Thorough Clean'
}';
cleanup.values = {0 1 2};
cleanup.def    = @(val)spm_get_defaults('old.preproc.output.cleanup', val{:});
% ---------------------------------------------------------------------
% output Output Files
% ---------------------------------------------------------------------
output         = cfg_branch;
output.tag     = 'output';
output.name    = 'Output Files';
output.val     = {GM WM CSF biascor cleanup };
output.help    = {
                  'This routine produces spatial normalisation parameters (*_seg_sn.mat files) by default. These can be used for writing spatially normalised versions of your data, via the "Normalise: Write" option. This mechanism may produce superior results than the "Normalise: Estimate" option (but probably not as good as those produced using DARTEL).'
                  ''
                  'In addition, it also produces files that can be used for doing inverse normalisation. If you have an image of regions defined in the standard space, then the inverse deformations can be used to warp these regions so that it approximately overlay your image. To use this facility, the bounding-box and voxel sizes should be set to non-finite values (e.g. [NaN NaN NaN] for the voxel sizes, and ones(2,3)*NaN for the bounding box. This would be done by the spatial normalisation module, which allows you to select a set of parameters that describe the nonlinear warps, and the images that they should be applied to.'
                  ''
                  'There are a number of options about what data you would like the routine to produce. The routine can be used for producing images of tissue classes, as well as bias corrected images. The native space option will produce a tissue class image (c*) that is in alignment with the original/* (see Figure \ref{seg1})*/.  You can also produce spatially normalised versions - both with (mwc*) and without (wc*) modulation/* (see Figure \ref{seg2})*/. The bounding box and voxel sizes of the spatially normalised versions are the same as that of the tissue probability maps with which they are registered. These can be used for doing voxel-based morphometry with (also see the ``Using DARTEL'' chapter of the manual). All you need to do is smooth them and do the stats (which means no more questions on the mailing list about how to do "optimized VBM").'
                  ''
                  'Modulation is to compensate for the effect of spatial normalisation.  When warping a series of images to match a template, it is inevitable that volumetric differences will be introduced into the warped images.  For example, if one subject''s temporal lobe has half the volume of that of the template, then its volume will be doubled during spatial normalisation. This will also result in a doubling of the voxels labelled grey matter.  In order to remove this confound, the spatially normalised grey matter (or other tissue class) is adjusted by multiplying by its relative volume before and after warping.  If warping results in a region doubling its volume, then the correction will halve the intensity of the tissue label. This whole procedure has the effect of preserving the total amount of grey matter signal in the normalised partitions.'
                  '/*\begin{figure} \begin{center} \includegraphics[width=140mm]{images/seg1} \end{center} \caption{Segmentation results. These are the results that can be obtained in the original space of the image (i.e. the results that are not spatially normalised). Top left: original image (X.img). Top right: bias corrected image (mX.img). Middle and bottom rows: segmented grey matter (c1X.img), white matter (c2X.img) and CSF (c3X.img). \label{seg1}} \end{figure} */'
                  '/*\begin{figure} \begin{center} \includegraphics[width=140mm]{images/seg2} \end{center} \caption{Segmentation results. These are the spatially normalised results that can be obtained (note that CSF data is not shown). Top row: The tissue probability maps used to guide the segmentation. Middle row: Spatially normalised tissue maps of grey and white matter (wc1X.img and wc2X.img). Bottom row: Modulated spatially normalised tissue maps of grey and white matter (mwc1X.img and mwc2X.img). \label{seg2}} \end{figure} */'
                  'A deformation field is a vector field, where three values are associated with each location in the field.  The field maps from co-ordinates in the normalised image back to co-ordinates in the original image.  The value of the field at co-ordinate [x y z] in the normalised space will be the co-ordinate [x'' y'' z''] in the original volume. The gradient of the deformation field at a co-ordinate is its Jacobian matrix, and it consists of a 3x3 matrix:'
                  ''
                  '%   /                      \'
                  '%   | dx''/dx  dx''/dy dx''/dz |'
                  '%   |                       |'
                  '%   | dy''/dx  dy''/dy dy''/dz |'
                  '%   |                       |'
                  '%   | dz''/dx  dz''/dy dz''/dz |'
                  '%   \                      /'
                  '/* \begin{eqnarray*}\begin{pmatrix}\frac{dx''}{dx} & \frac{dx''}{dy} & \frac{dx''}{dz}\cr\frac{dy''}{dx} & \frac{dy''}{dy} & \frac{dy''}{dz}\cr\frac{dz''}{dx} & \frac{dz''}{dy} & \frac{dz''}{dz}\cr\end{pmatrix}\end{eqnarray*}*/'
                  'The value of dx''/dy is a measure of how much x'' changes if y is changed by a tiny amount. The determinant of the Jacobian is the measure of relative volumes of warped and unwarped structures.  The modulation step simply involves multiplying by the relative volumes /*(see Figure \ref{seg2})*/.'
}';
% ---------------------------------------------------------------------
% tpm Tissue probability maps
% ---------------------------------------------------------------------
tpm         = cfg_files;
tpm.tag     = 'tpm';
tpm.name    = 'Tissue probability maps';
tpm.help    = {
               'Select the tissue probability images. These should be maps of grey matter, white matter and cerebro-spinal fluid probability. A nonlinear deformation field is estimated that best overlays the tissue probability maps on the individual subjects'' image. The default tissue probability maps are modified versions of the ICBM Tissue Probabilistic Atlases.These tissue probability maps are kindly provided by the International Consortium for Brain Mapping, John C. Mazziotta and Arthur W. Toga. http://www.loni.ucla.edu/ICBM/ICBM_TissueProb.html. The original data are derived from 452 T1-weighted scans, which were aligned with an atlas space, corrected for scan inhomogeneities, and classified into grey matter, white matter and cerebrospinal fluid. These data were then affine registered to the MNI space and downsampled to 2mm resolution.'
               ''
               'Rather than assuming stationary prior probabilities based upon mixing proportions, additional information is used, based on other subjects'' brain images.  Priors are usually generated by registering a large number of subjects together, assigning voxels to different tissue types and averaging tissue classes over subjects. Three tissue classes are used: grey matter, white matter and cerebro-spinal fluid. A fourth class is also used, which is simply one minus the sum of the first three. These maps give the prior probability of any voxel in a registered image being of any of the tissue classes - irrespective of its intensity.'
               ''
               'The model is refined further by allowing the tissue probability maps to be deformed according to a set of estimated parameters. This allows spatial normalisation and segmentation to be combined into the same model. This implementation uses a low-dimensional approach, which parameterises the deformations by a linear combination of about a thousand cosine transform bases. This is not an especially precise way of encoding deformations, but it can model the variability of overall brain shape. Evaluations by Hellier et al have shown that this simple model can achieve a registration accuracy comparable to other fully automated methods with many more parameters.'
}';
tpm.filter  = 'image';
tpm.dir     = fileparts(mfilename('fullpath'));
tpm.ufilter = '.*';
tpm.num     = [3 3];
tpm.def     = @(val)spm_get_defaults('old.preproc.tpm', val{:});
% ---------------------------------------------------------------------
% ngaus Gaussians per class
% ---------------------------------------------------------------------
ngaus         = cfg_entry;
ngaus.tag     = 'ngaus';
ngaus.name    = 'Gaussians per class';
ngaus.help    = {'The number of Gaussians used to represent the intensity distribution for each tissue class can be greater than one. In other words, a tissue probability map may be shared by several clusters. The assumption of a single Gaussian distribution for each class does not hold for a number of reasons. In particular, a voxel may not be purely of one tissue type, and instead contain signal from a number of different tissues (partial volume effects). Some partial volume voxels could fall at the interface between different classes, or they may fall in the middle of structures such as the thalamus, which may be considered as being either grey or white matter. Various other image segmentation approaches use additional clusters to model such partial volume effects. These generally assume that a pure tissue class has a Gaussian intensity distribution, whereas intensity distributions for partial volume voxels are broader, falling between the intensities of the pure classes. Unlike these partial volume segmentation approaches, the model adopted here simply assumes that the intensity distribution of each class may not be Gaussian, and assigns belonging probabilities according to these non-Gaussian distributions. Typical numbers of Gaussians could be two for grey matter, two for white matter, two for CSF, and four for everything else.'};
ngaus.strtype = 'n';
ngaus.num     = [4 1];
ngaus.def     = @(val)spm_get_defaults('old.preproc.ngaus', val{:});
% ---------------------------------------------------------------------
% regtype Affine Regularisation
% ---------------------------------------------------------------------
regtype         = cfg_menu;
regtype.tag     = 'regtype';
regtype.name    = 'Affine Regularisation';
regtype.help    = {
                   'The procedure is a local optimisation, so it needs reasonable initial starting estimates. Images should be placed in approximate alignment using the Display function of SPM before beginning. A Mutual Information affine registration with the tissue probability maps (D''Agostino et al, 2004) is used to achieve approximate alignment. Note that this step does not include any model for intensity non-uniformity. This means that if the procedure is to be initialised with the affine registration, then the data should not be too corrupted with this artifact.If there is a lot of intensity non-uniformity, then manually position your image in order to achieve closer starting estimates, and turn off the affine registration.'
                   ''
                   'Affine registration into a standard space can be made more robust by regularisation (penalising excessive stretching or shrinking).  The best solutions can be obtained by knowing the approximate amount of stretching that is needed (e.g. ICBM templates are slightly bigger than typical brains, so greater zooms are likely to be needed). For example, if registering to an image in ICBM/MNI space, then choose this option.  If registering to a template that is close in size, then select the appropriate option for this.'
}';
regtype.labels = {
                  'No Affine Registration'
                  'ICBM space template - European brains'
                  'ICBM space template - East Asian brains'
                  'Average sized template'
                  'No regularisation'
}';
regtype.values = {
                  ''
                  'mni'
                  'eastern'
                  'subj'
                  'none'
}';
regtype.def     = @(val)spm_get_defaults('old.preproc.regtype', val{:});
% ---------------------------------------------------------------------
% warpreg Warping Regularisation
% ---------------------------------------------------------------------
warpreg         = cfg_entry;
warpreg.tag     = 'warpreg';
warpreg.name    = 'Warping Regularisation';
warpreg.help    = {'The objective function for registering the tissue probability maps to the image to process, involves minimising the sum of two terms. One term gives a function of how probable the data is given the warping parameters. The other is a function of how probable the parameters are, and provides a penalty for unlikely deformations. Smoother deformations are deemed to be more probable. The amount of regularisation determines the tradeoff between the terms. Pick a value around one.  However, if your normalised images appear distorted, then it may be an idea to increase the amount of regularisation (by an order of magnitude). More regularisation gives smoother deformations, where the smoothness measure is determined by the bending energy of the deformations. '};
warpreg.strtype = 'e';
warpreg.num     = [1 1];
warpreg.def     = @(val)spm_get_defaults('old.preproc.warpreg', val{:});
% ---------------------------------------------------------------------
% warpco Warp Frequency Cutoff
% ---------------------------------------------------------------------
warpco         = cfg_entry;
warpco.tag     = 'warpco';
warpco.name    = 'Warp Frequency Cutoff';
warpco.help    = {'Cutoff of DCT bases.  Only DCT bases of periods longer than the cutoff are used to describe the warps. The number actually used will depend on the cutoff and the field of view of your image. A smaller cutoff frequency will allow more detailed deformations to be modelled, but unfortunately comes at a cost of greatly increasing the amount of memory needed, and the time taken.'};
warpco.strtype = 'e';
warpco.num     = [1 1];
warpco.def     = @(val)spm_get_defaults('old.preproc.warpco', val{:});
% ---------------------------------------------------------------------
% biasreg Bias regularisation
% ---------------------------------------------------------------------
biasreg         = cfg_menu;
biasreg.tag     = 'biasreg';
biasreg.name    = 'Bias regularisation';
biasreg.help    = {
                   'MR images are usually corrupted by a smooth, spatially varying artifact that modulates the intensity of the image (bias). These artifacts, although not usually a problem for visual inspection, can impede automated processing of the images.'
                   ''
                   'An important issue relates to the distinction between intensity variations that arise because of bias artifact due to the physics of MR scanning, and those that arise due to different tissue properties.  The objective is to model the latter by different tissue classes, while modelling the former with a bias field. We know a priori that intensity variations due to MR physics tend to be spatially smooth, whereas those due to different tissue types tend to contain more high frequency information. A more accurate estimate of a bias field can be obtained by including prior knowledge about the distribution of the fields likely to be encountered by the correction algorithm. For example, if it is known that there is little or no intensity non-uniformity, then it would be wise to penalise large values for the intensity non-uniformity parameters. This regularisation can be placed within a Bayesian context, whereby the penalty incurred is the negative logarithm of a prior probability for any particular pattern of non-uniformity.'
}';
biasreg.labels = {
                  'no regularisation (0)'
                  'extremely light regularisation (0.00001)'
                  'very light regularisation (0.0001)'
                  'light regularisation (0.001)'
                  'medium regularisation (0.01)'
                  'heavy regularisation (0.1)'
                  'very heavy regularisation (1)'
                  'extremely heavy regularisation (10)'
}';
biasreg.values = {0 1e-5 .0001 .001 .01 .1 1 10};
biasreg.def     = @(val)spm_get_defaults('old.preproc.biasreg', val{:});
% ---------------------------------------------------------------------
% biasfwhm Bias FWHM
% ---------------------------------------------------------------------
biasfwhm         = cfg_menu;
biasfwhm.tag     = 'biasfwhm';
biasfwhm.name    = 'Bias FWHM';
biasfwhm.help    = {'FWHM of Gaussian smoothness of bias. If your intensity non-uniformity is very smooth, then choose a large FWHM. This will prevent the algorithm from trying to model out intensity variation due to different tissue types. The model for intensity non-uniformity is one of i.i.d. Gaussian noise that has been smoothed by some amount, before taking the exponential. Note also that smoother bias fields need fewer parameters to describe them. This means that the algorithm is faster for smoother intensity non-uniformities.'};
biasfwhm.labels = {
                   '30mm cutoff'
                   '40mm cutoff'
                   '50mm cutoff'
                   '60mm cutoff'
                   '70mm cutoff'
                   '80mm cutoff'
                   '90mm cutoff'
                   '100mm cutoff'
                   '110mm cutoff'
                   '120mm cutoff'
                   '130mm cutoff'
                   '140mm cutoff'
                   '150mm cutoff'
                   'No correction'
}';
biasfwhm.values = {30 40 50 60 70 80 90 100 110 120 130 140 150 Inf};
biasfwhm.def     = @(val)spm_get_defaults('old.preproc.biasfwhm', val{:});
% ---------------------------------------------------------------------
% samp Sampling distance
% ---------------------------------------------------------------------
samp         = cfg_entry;
samp.tag     = 'samp';
samp.name    = 'Sampling distance';
samp.help    = {'The approximate distance between sampled points when estimating the model parameters. Smaller values use more of the data, but the procedure is slower.'};
samp.strtype = 'e';
samp.num     = [1 1];
samp.def     = @(val)spm_get_defaults('old.preproc.samp', val{:});
% ---------------------------------------------------------------------
% msk Masking image
% ---------------------------------------------------------------------
msk         = cfg_files;
msk.tag     = 'msk';
msk.name    = 'Masking image';
msk.val     = {{''}};
msk.help    = {'The segmentation can be masked by an image that conforms to the same space as the images to be segmented.  If an image is selected, then it must match the image(s) voxel-for voxel, and have the same voxel-to-world mapping.  Regions containing a value of zero in this image do not contribute when estimating the various parameters. '};
msk.filter  = 'image';
msk.ufilter = '.*';
msk.num     = [0 1];
% ---------------------------------------------------------------------
% opts Custom
% ---------------------------------------------------------------------
opts         = cfg_branch;
opts.tag     = 'opts';
opts.name    = 'Custom';
opts.val     = {tpm ngaus regtype warpreg warpco biasreg biasfwhm samp msk };
opts.help    = {'Various options can be adjusted in order to improve the performance of the algorithm with your data.  Knowing what works best should be a matter of empirical exploration.  For example, if your data has very little intensity non-uniformity artifact, then the bias regularisation should be increased.  This effectively tells the algorithm that there is very little bias in your data, so it does not try to model it.'};
% ---------------------------------------------------------------------
% oldseg Old Segment
% ---------------------------------------------------------------------
preproc         = cfg_exbranch;
preproc.tag     = 'oldseg';
preproc.name    = 'Old Segment';
preproc.val     = {data output opts};
preproc.help    = {
                   'Segment, bias correct and spatially normalise - all in the same model/* \cite{ashburner05}*/. This function can be used for bias correcting, spatially normalising or segmenting your data. Note that this module needs the images to be roughly aligned with the tissue probability maps before you begin.  If strange results are obtained, then this is usually because the images were poorly aligned beforehand.  The Display option can be used to manually reposition the images so that the AC is close to coordinate 0,0,0 (within a couple of cm) and the orientation is within a few degrees of the tissue probability map data.'
                   ''
                   'Many investigators use tools within older versions of SPM for a technique that has become known as "optimised" voxel-based morphometry (VBM). VBM performs region-wise volumetric comparisons among populations of subjects. It requires the images to be spatially normalised, segmented into different tissue classes, and smoothed, prior to performing statistical tests/* \cite{wright_vbm,am_vbmreview,ashburner00b,john_should}*/. The "optimised" pre-processing strategy involved spatially normalising subjects'' brain images to a standard space, by matching grey matter in these images, to a grey matter reference.  The historical motivation behind this approach was to reduce the confounding effects of non-brain (e.g. scalp) structural variability on the registration. Tissue classification in older versions of SPM required the images to be registered with tissue probability maps. After registration, these maps represented the prior probability of different tissue classes being found at each location in an image.  Bayes rule can then be used to combine these priors with tissue type probabilities derived from voxel intensities, to provide the posterior probability.'
                   ''
                   'This procedure was inherently circular, because the registration required an initial tissue classification, and the tissue classification requires an initial registration.  This circularity is resolved here by combining both components into a single generative model. This model also includes parameters that account for image intensity non-uniformity. Estimating the model parameters (for a maximum a posteriori solution) involves alternating among classification, bias correction and registration steps. This approach provides better results than simple serial applications of each component.'
                   ''
                   'Note that multi-spectral segmentation (e.g. from a registered T1 and T2 image) is not yet implemented, but is planned for a future SPM version.'
}';
preproc.prog = @spm_run_preproc;
preproc.vout = @vout;
%------------------------------------------------------------------------

%------------------------------------------------------------------------
function dep = vout(job)
opts  = job.output;

cdep(1)            = cfg_dep;
cdep(1).sname      = 'Norm Params Subj->MNI';
cdep(1).src_output = substruct('()',{1}, '.','snfile','()',{':'});
cdep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
cdep(2)            = cfg_dep;
cdep(2).sname      = 'Norm Params MNI->Subj';
cdep(2).src_output = substruct('()',{1}, '.','isnfile','()',{':'});
cdep(2).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
if opts.biascor,
    cdep(end+1)          = cfg_dep;
    cdep(end).sname      = 'Bias Corr Images';
    cdep(end).src_output = substruct('()',{1}, '.','biascorr','()',{':'});
    cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end;
sonames = {'GM','WM','CSF'};
for k1=1:3,
    if ~strcmp(opts.(sonames{k1}),'<UNDEFINED>')
        if opts.(sonames{k1})(3),
            cdep(end+1)          = cfg_dep;
            cdep(end).sname      = sprintf('c%d Images',k1);
            cdep(end).src_output = substruct('()',{1}, '.',sprintf('c%d',k1),'()',{':'});
            cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
        end;
        if opts.(sonames{k1})(2),
            cdep(end+1)          = cfg_dep;
            cdep(end).sname      = sprintf('wc%d Images',k1);
            cdep(end).src_output = substruct('()',{1}, '.',sprintf('wc%d',k1),'()',{':'});
            cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
        end;
        if opts.(sonames{k1})(1),
            cdep(end+1)          = cfg_dep;
            cdep(end).sname      = sprintf('mwc%d Images',k1);
            cdep(end).src_output = substruct('()',{1}, '.',sprintf('mwc%d',k1),'()',{':'});
            cdep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
        end;
    end
end;
dep = cdep;

author: spm


function [ compressIndex_new ] = anal_reindex_last(compressDir,  reIndexFile )
%ANAL_REINDEX_LAST Summary of this function goes here
%   Detailed explanation goes here

    fileExt = '*.txt';
    compress_files = dir(fullfile(compressDir,fileExt));
    compressIndex = [];
    for i = 0:length(compress_files)-1

        fileCnt = num2str(i);
        fileName = [compressDir, fileCnt, '.txt'];

        file_t = fopen(fileName);
        point_ID = fscanf(file_t, '%d');
        fclose(file_t);

        compressIndex = union(compressIndex, point_ID);

    end

    reIndex = dlmread(reIndexFile);

    isInIndex = ismember(reIndex(:,2), compressIndex);
    isSaved = find(isInIndex==1);
    old_index = reIndex(:,1);
    compressIndex_new = old_index(isSaved);

end


text:
function [ n_data, n, x, fx ] = l_polynomial_values ( n_data )

%*****************************************************************************80
%
%% L_POLYNOMIAL_VALUES returns some values of the Laguerre polynomial.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      LaguerreL[n,x]
%
%  First terms:
%
%      1
%     -X    +  1
%   (  X^2 -  4 X     +  2 ) / 2
%   ( -X^3 +  9 X^2 -  18 X    +    6 ) / 6
%   (  X^4 - 16 X^3 +  72 X^2 -   96 X +      24 ) / 24
%   ( -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -  600 x    +  120 ) / 120
%   (  X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720 ) / 720
%   ( -X^7 + 49 X^6 - 882 X^5 + 7350 X^4 - 29400 X^3 + 52920 X^2 - 35280 X
%     + 5040 ) / 5040
%
%  Recursion:
%
%    L(0,X) = 1,
%    L(1,X) = 1-X,
%    N * L(N,X) = (2*N-1-X) * L(N-1,X) - (N-1) * L(N-2,X)
%
%  Orthogonality:
%
%    Integral ( 0 <= X < oo ) exp ( - X ) * L(N,X) * L(M,X) dX
%    = 0 if N /= M
%    = 1 if N == M
%
%  Special values:
%
%    L(N,0) = 1.
%
%  Relations:
%
%    L(N,X) = (-1)^N / N% * exp ( x ) * (d/dx)^n ( exp ( - x ) * x^n )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the polynomial.
%
%    Output, real X, the point where the polynomial is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 17;

  fx_vec = [ ...
      0.1000000000000000E+01, ...
      0.0000000000000000E+00, ...
     -0.5000000000000000E+00, ...
     -0.6666666666666667E+00, ...
     -0.6250000000000000E+00, ...
     -0.4666666666666667E+00, ...
     -0.2569444444444444E+00, ...
     -0.4047619047619048E-01, ...
      0.1539930555555556E+00, ...
      0.3097442680776014E+00, ...
      0.4189459325396825E+00, ...
      0.4801341790925124E+00, ...
      0.4962122235082305E+00, ...
     -0.4455729166666667E+00, ...
      0.8500000000000000E+00, ...
     -0.3166666666666667E+01, ...
      0.3433333333333333E+02  ];

  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10, 11, ...
    12,  5,  5, ...
     5,  5 ];

  x_vec = [ ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     0.5E+00, ...
     3.0E+00, ...
     5.0E+00, ...
     1.0E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function t_opf_dc_ipopt(quiet)
%T_OPF_DC_MIPS  Tests for DC optimal power flow using MIPS solver.

%   MATPOWER
%   Copyright (c) 2004-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

num_tests = 43;

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end

t0 = 'DC OPF (IPOPT): ';
mpopt = mpoption('out.all', 0, 'verbose', verbose);
mpopt = mpoption(mpopt, 'opf.dc.solver', 'IPOPT');

if have_feature('ipopt')
    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved DC power flow case from MAT-file
    load soln9_dcopf;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    t = [t0 'w/angle diff lims : '];
    mpc = loadcase(casefile);
    mpc.branch(4, ANGMAX) = 3;
    mpc.branch(7, ANGMIN) = -4.5;
    r = rundcopf(mpc, mpopt);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(   f, 6456.7213, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,PG        ),    [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    e = zeros(size(branch, 1), 1);
    e(4) = 297.83776;
    e(7) = -26.94788;
    t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']);

    t = [t0 'w/ignored angle diff lims : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    r = rundcopf(mpc, mpopt1);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% two new z variables
    %%      0 <= z1, P3 - P1 <= z1
    %%      0 <= z2, P3 - P2 <= z2
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 4, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 4, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 3, [t 'user costs']);

    %% with A and N sized for AC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 2 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 4, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 4, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 3, [t 'user costs']);

    t = [t0 'infeasible : '];
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);   %% Pg1 + Pg2
    mpc.u = Inf;
    mpc.l = 600;
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(~success, [t 'no success']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc = loadcase(casefile);
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = rundcopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end
else
    t_skip(num_tests, 'IPOPT not available');
end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [Y,W,SetupStruc] = Process_ILRMA_PF(s,Transfer,SetupStruc)
K = SetupStruc.ILRMA_PF.K;
hop = SetupStruc.ILRMA_PF.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.ILRMA_PF.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(Num,frame_N,K);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
X_sp = zeros(Num,frame_N,K_m);
W_ILRMA = zeros(Num,Num,K_m);
V_sp = zeros(Num,N,K_m);
theta = 10^-6;
for i = 1:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Initialize W by PCA
    [E,D] = PCA(X_f,1,Num);
    V = sqrt(D)\E';
    V_sp(:,:,i) = V;
    X_sp(:,:,i) = V*X_f;
    %%%%%%%%%%%% Adjust amplitude of 'w'
    W_o = eye(Num);
    y_f = W_o*V*X_f;
    W_ILRMA(:,:,i) = W_o;
    Y_f(:,:,i) = y_f;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% ILRMA iterations
epsi = 1e-6;
max_iteration = 300;
L = 2*Num;  %%%%% the number of NMF basis
Z = max(rand(Num,L),epsi);
Z = Z./(ones(Num)*Z);
T = max(rand(K_m,L),epsi);
V = max(rand(L,frame_N),epsi);
R = zeros(K_m,frame_N,Num);
for i = 1:Num
    R(:,:,i) = (ones(K_m,1)*Z(i,:)).*T*V;
end
P = permute(Y_f(:,:,1:K_m),[3,2,1]);
P = abs(P).^2;
% K_sp = sqrt(K_m);
pObj = inf;
bZ = zeros(size(Z));
bT = zeros(size(T));
bV = zeros(size(V));
lamd = zeros(Num,1);
A = zeros(1001,2)-1; %%%% Show the decrease of the value of cost funtion, ILRMA max iterations 1000
for iteration = 1:max_iteration
    %%%%%% NMF with partioning function Z
    for i = 1:Num
        bZ(i,:) = sqrt((T'*(P(:,:,i).*(R(:,:,i).^(-2)))).*V*ones(frame_N,1)./((T'*R(:,:,i).^(-1)).*V*ones(frame_N,1)))';
    end
    Z = max(Z./(ones(Num)*Z),epsi);
    for i = 1:Num
        R(:,:,i) = (ones(K_m,1)*Z(i,:)).*T*V;
    end
    for i = 1:K_m
        P_temp = permute(P(i,:,:),[2 3 1]);
        R_temp = permute(R(i,:,:),[2 3 1]);
        bT(i,:) = sqrt((V*(P_temp.*(R_temp.^(-2)))).*Z'*ones(Num,1)./((V*R_temp.^(-1)).*Z'*ones(Num,1)))';
    end
    T = max(T.*bT,epsi);
    for i = 1:Num
        R(:,:,i) = (ones(K_m,1)*Z(i,:)).*T*V;
    end
    for i = 1:frame_N
        P_temp = permute(P(:,i,:),[1 3 2]);
        R_temp = permute(R(:,i,:),[1 3 2]);
        bV(:,i) = sqrt((T'*(P_temp.*(R_temp.^(-2)))).*Z'*ones(Num,1)./((T'*R_temp.^(-1)).*Z'*ones(Num,1)));
    end
    V = max(V.*bV,epsi);
    for i = 1:Num
        R(:,:,i) = (ones(K_m,1)*Z(i,:)).*T*V;
    end
    %%%%% AuxIVA
    dlw = 0;
    for i = 1:K_m
        W = W_ILRMA(:,:,i);
        X_f = X_sp(:,:,i);
       dlw = dlw +log(abs(det(W))+epsi);
        for i_n = 1:Num
            G_ = permute(R(i,:,i_n),[3 2 1]);
            G_ = repmat(G_+epsi,Num,1);
            Vk = (X_f./G_)*X_f'/frame_N;
            if rcond(Vk)<theta
                Vk = Vk+eye(Num)*min(eig(Vk))*theta;
            end
            wk = inv(W*Vk);
            wk = wk(:,i_n);
            wk = wk/(sqrt(wk'*Vk*wk)+epsi);
            W(i_n,:) = wk';
        end
        W_ILRMA(:,:,i) = W;
        Y_f(:,:,i) = W*X_f;
    end
    P = permute(Y_f(:,:,1:K_m),[3,2,1]);
    P = abs(P).^2;
    Obj = ((sum(sum(sum(log(R+epsi))))+sum(sum(sum(P./R))))/frame_N-2*dlw)/(Num*K_m);
    dObj = pObj-Obj;
    pObj = Obj;
    A(iteration,:) = [Obj,abs(dObj)/abs(Obj)];
    if(abs(dObj)/abs(Obj)<theta)
       break;
    end
    for i = 1:Num
        lamda = sqrt(sum(sum(P(:,:,i))/(frame_N*K_m)));
        lamd(i) = lamda;
        for i_k = 1:K_m
            W_ILRMA(i,:,i_k) = W_ILRMA(i,:,i_k)/lamda;
            for l = 1:L
                T(i_k,l) = T(i_k,l)*sum(Z(:,l))/lamda^2;
            end
        end
        P(:,:,i) = P(:,:,i)/lamda^2;
        R(:,:,i) = R(:,:,i)/lamda^2;
    end
    for l = 1:L
        Z_temp = sum(Z(:,l)./lamd.^(2));
        for i = 1:Num
            Z(i,l) = Z(i,l)/lamd(i)^2/Z_temp;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Post processing
W = zeros(Num,N,K_m);
Y_f(:,:,1) = zeros(Num,frame_N);
for i = 2:K_m
    W_inv = pinv(W_ILRMA(:,:,i)*V_sp(:,:,i));
    for ii = 1:Num
        Y_f(ii,:,i) = Y_f(ii,:,i)*W_inv(1,ii);
        W_ILRMA(ii,:,i) = W_ILRMA(ii,:,i)*W_inv(1,ii);
    end
    W(:,:,i) = W_ILRMA(:,:,i)*V_sp(:,:,i);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if(i~=K_m)
        Y_f(:,:,K+2-i) = conj(Y_f(:,:,i));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    y_temp = permute(Y_f(i,:,:),[3 2 1]);
    Y(:,i) = overlapadd(real(ifft(y_temp))',win,hop);
end
return;

repo: Voice-Separation-and-Enhancement


text:
function [p ener a enerFF aFF max_doseV doseV_wt] = optSourceModel_depth6cm_blurring(planC, energy, numBin, extraBin, depth, dose, ...
    depth1, x1, profile1,  depth2, x2, profile2, depth3, x3, profile3, depth4, x4, profile4, LB, UB);
% JC Dec 09,2005
% Optimize photon spectrum, using (Fatigue-life)*Fermi distribution
% Usage:
% energy = energy of the photon beam
% numBin = number of bins to divide the energy range [0 energy]
% depth is a row vector as the depth
% dose is the corresponding dose at each depth
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% Underline assumptions:
% The calculated dose sets are the same size as the uniformizedCTscan.
% 1. Interpolation the calculated dose on to the respectively dose profiles
% coordinates may cause problem: Since the individual dose profile
% measurements do not necessarily have the same coordinates/# of points,
%

% 2. The input depth/x coordinates are in cm
%    The depth dose (PDD) is normalized. dmax = 1
%    The lateral dose profiles at various depth are also normalized. Each
%       has a maximal value as 1. It will be scaled according to PDD.

% 3. The phantom is at least 25cm in the depth direction. Otherwise, the
% interpolation of PDD will have out-flow, meaning the non-meaningful
% values, N/A. An error will be thrown out.
%   Also, the input PDD's depth may not start from zero, this may generate
%   N/A if the yVals starts from zero. So chop yVals according to PDD's
%   depth.

% 4. All input measured data are in a nx1 vector, i.e. a column vector

% 5 Add more input parameters: adjustable to different machines
% LB=[0.6  -0.1 1.2   100  1  0.1  0.005    0.005];
% UB=[1.5  0.1   10  1000  4   0.25  0.05  0.05];

% Input parameters descriptions:
% LB, UB. Upper and lower bound for the parameters.
% p, the optimazied parameters.
% p(1:4) The parameters for "Fagitue-Fermi" function, to model the primary photon
% spectrum.
% p(5) The parameter to "shrink" the primary photon spectrum
% to get the spectrum for the extra-focal (FF) components.
% p(6) The parameter to "scale" the extra-focal (FF) components,
% to determine the relative fluence contribution of the extra-focal (FF)
% to the primary photon.
% p(7) used to determine "the slope of the horn", not used anymor i.e. It is how many "pixels"/"grid".
% The dimension of the ICsigma in (cm) should be p(9)*e.
% p(8) the electron contamination (fluence relative to the primary photon.
% p(9) the "ICsigma", to convolve the MC dose, to model the spearing of the dose by IC chamber.
% If the resolution of the dose calc. grid changes, i.e.
% The dimension of the sigma is p(9)*dx, where "dx" is the grid spacing for
% dose calc.

% 6. Add another parameter for the optimization, sigma. A gaussian filter
% to blur the penumbra.

% Sept 21, 2006
% Add electron contamination components.
% For PDD, includes depth < 5cm

% Nov 27, 2006
% correct
% doseP2_arrayFF(:,i) = dose3D(find(yVals == yDepth3),[xIndexStart:xIndexEnd],int16(s));
% doseP2_arrayFF(:,i) = dose3D(find(yVals == yDepth2),[xIndexStart:xIndexEnd],int16(s));

% Get the isocenter coordinate
bf = planC{7}.FractionGroupSequence.Item_1.ReferencedBeamSequence.Item_1;
bs = planC{7}.BeamSequence.Item_1;
iC = bs.ControlPointSequence.Item_1.IsocenterPosition;

position = {planC{7}.PatientSetupSequence.(['Item_' num2str(1)]).PatientPosition};

if strcmpi(position, 'HFP')
    isocenter.x = iC(1)/10;
    isocenter.y = iC(2)/10;
    isocenter.z = -iC(3)/10;
else
    isocenter.x = iC(1)/10;
    isocenter.y = -iC(2)/10;
    isocenter.z = -iC(3)/10;
end
% So iC is the isocenter. Need to find indices from the (x, y, z) values.

[xVals, yVals, zVals] = getUniformScanXYZVals(planC{3});
% Need to confirm that iC(2), i.e. the y value is right on the skin
% surface. Thus, the depth for dose along the center ray can be determined.
j = find(strcmpi({planC{4}.structureName}, 'Skin'));

if isempty(j)
    j = find(strcmpi({planC{4}.structureName}, 'Body'));
end

if isempty(j)
    error('Skin structure must be defined to input to DPM.');
else
    j = j(1);
end

skinMask=getUniformStr(j,planC);

%confirm that the isocenter.y is right on the skin surface
bbox=boundingBox(skinMask)
%[r,c,s] = xyztom(isocenter.x,isocenter.y,isocenter.z, planC)
[r,c,s] = xyztom(isocenter.x,isocenter.y,isocenter.z, 1, planC)

%Get the indices of where yVals==depth(end)/10 ; Failed???
%[r,c,s] = xyztom(isocenter.x,yVals(find(yVals==isocenter.y-min(isocenter.y-yVals))),isocenter.z, planC)

% Get the objective dose from input(Golden Beam Data), not from planC.

% get the yVals for different depth, to pick out the appropriate lateral dose profiles in the calcualted dose sets.
dy = yVals(2)-yVals(1);  % dy is < 0, because of y increases as row number decreases
yDepth0 = yVals(bbox(1));
yDepth1 = yVals(bbox(1)+int16(-depth1/dy));
yDepth2 = yVals(bbox(1)+int16(-depth2/dy));
yDepth3 = yVals(bbox(1)+int16(-depth3/dy));
yDepth4 = yVals(bbox(1)+int16(-depth4/dy));
yDepth25cm = yVals(bbox(1)+int16(-25/dy));

%depthTPS can't start from 0, which will cause intepolated data N/A.
depthTPS = [ceil(min(depth)/(-dy))*(-dy) : -dy: yVals(bbox(1))-yDepth25cm];
PDD = interp1(depth, dose, depthTPS);
doseV_obj = PDD';
doseV_array = zeros(length(PDD),numBin+extraBin);
max_doseV = zeros(numBin+extraBin,1);

% GET the start and end indices objective for dose profile at depth 5cm
xStart = max(min([x1 x2 x3 x4]));
xEnd = min(max([x1 x2 x3 x4]));
% use margin 1 to make sure the results xTPS is within the measured region.
xIndexStart = find(abs(xStart - xVals) == min(abs(xStart - xVals))) + 1;
xIndexEnd = find(abs(xEnd - xVals) == min(abs(xEnd - xVals))) - 1 ;

% Should correct xTPS for isocenter.x? Does, then xTPS maybe outside of the
% valid region.
 xTPS = [xVals(xIndexStart)-isocenter.x : (xVals(2)-xVals(1)): xVals(xIndexEnd)-isocenter.x];
% Currently, do not correct.
% xTPS = [xVals(xIndexStart): (xVals(2)-xVals(1)): xVals(xIndexEnd)];

doseProfile1 = interp1(x1, profile1, xTPS);
% Since the input measured data is the asdfScale it to be consistent with PDD
% Needed for Golden Beam data. Not for Murty's measurement.
doseProfile1_obj = doseProfile1' * PDD(find(abs(depth1-depthTPS) == min(abs(depth1-depthTPS))));
% Creat array to hold the dose profiles for the different energy bins.
doseP1_array = zeros(length(doseProfile1_obj), numBin+extraBin);

doseProfile2 = interp1(x2, profile2, xTPS);
% Scale it to be consistent with PDD
doseProfile2_obj = doseProfile2' * PDD(find(abs(depth2-depthTPS) == min(abs(depth2-depthTPS))));
doseP2_array = zeros(length(doseProfile2_obj), numBin+extraBin);

% GET the objective for dose profile at depth 10cm
doseProfile3 = interp1(x3, profile3, xTPS);
% Scale it to be consistent with PDD
doseProfile3_obj = doseProfile3' * PDD(find(abs(depth3-depthTPS) == min(abs(depth3-depthTPS))));
doseP3_array = zeros(length(doseProfile3_obj), numBin+extraBin);

% GET the objective for dose profile at depth 20cm
doseProfile4 = interp1(x4, profile4, xTPS);
% Scale it to be consistent with PDD
doseProfile4_obj = doseProfile4' * PDD(find(abs(depth4-depthTPS) == min(abs(depth4-depthTPS))));
doseP4_array = zeros(length(doseProfile4_obj), numBin+extraBin);

figure; plot(depthTPS, PDD);
% There's shift of doseProfile2_obj, x as half voxel size. 0.1875/2 =
% 0.094cm
hold on; plot(xTPS, doseProfile1, 'c', xTPS, doseProfile2, 'b', xTPS, doseProfile3, 'r', xTPS, doseProfile4, 'g');
legend('PDD', num2str(depth1), num2str(depth2), num2str(depth3), num2str(depth4));

% Read in the dose from the DPM calculation.
for i = 1:(numBin+extraBin),
    filename = ['dose3D_', num2str(i*energy/numBin),'MV.mat'];
    load (filename, 'dose3D');
    doseV_array(:,i) = dose3D([(bbox(1)+ceil(min(depth)/(-dy))):find(yVals == yDepth25cm)],int16(c),int16(s));
    doseP1_array(:,i) = dose3D(find(yVals == yDepth1),[xIndexStart:xIndexEnd],int16(s));
    doseP2_array(:,i) = dose3D(find(yVals == yDepth2),[xIndexStart:xIndexEnd],int16(s));
    doseP3_array(:,i) = dose3D(find(yVals == yDepth3),(xIndexStart:xIndexEnd),int16(s));
    doseP4_array(:,i) = dose3D(find(yVals == yDepth4),(xIndexStart:xIndexEnd),int16(s));

  end

% currentDir = pwd;
% cd ./FlatFilter
%NOW, for the flattenning filter part, need to do the exactly same thing,
%i.e. get PDD and doseP5cm, doseP10cm, doseP20cm....
numBinFF = ceil((numBin+extraBin)/2); %15
% JC. Use all energy bins.
numBinFF = numBin+extraBin; %15
doseV_arrayFF = zeros(length(PDD),numBinFF);
doseP1_arrayFF = zeros(length(doseProfile1), numBinFF);
doseP2_arrayFF = zeros(length(doseProfile2), numBinFF);
doseP3_arrayFF = zeros(length(doseProfile3), numBinFF);
doseP4_arrayFF = zeros(length(doseProfile4), numBinFF);
% Read in the dose from the DPM calculation, for the secondary source, i.e. Flattening Filter,FF.
for i = 1 : numBinFF,
    filename = ['dose3D_FF_', num2str((double(i*energy)/double(numBin))),'MV.mat'];
    load (filename, 'dose3D');
    doseV_arrayFF(:,i) = dose3D([(bbox(1)+ceil(min(depth)/(-dy))):find(yVals == yDepth25cm)],int16(c),int16(s));
    doseP1_arrayFF(:,i) = dose3D(find(yVals == yDepth1),[xIndexStart:xIndexEnd],int16(s));
    doseP2_arrayFF(:,i) = dose3D(find(yVals == yDepth2),[xIndexStart:xIndexEnd],int16(s));
    doseP3_arrayFF(:,i) = dose3D(find(yVals == yDepth3),(xIndexStart:xIndexEnd),int16(s));
    doseP4_arrayFF(:,i) = dose3D(find(yVals == yDepth4),(xIndexStart:xIndexEnd),int16(s));
end
% cd(currentDir);

% Read the dose of the electron contamination.
load dose3D_elec.mat dose3D
doseV_e = dose3D([(bbox(1)+ceil(min(depth)/(-dy))):find(yVals == yDepth25cm)],int16(c),int16(s));
% Need doseP5cm_e etc. is a collumn vector.
doseP1_e = dose3D(find(yVals == yDepth1),[xIndexStart:xIndexEnd],int16(s))';
doseP2_e = dose3D(find(yVals == yDepth2),[xIndexStart:xIndexEnd],int16(s))';
doseP3_e = dose3D(find(yVals == yDepth3),(xIndexStart:xIndexEnd),int16(s))';
doseP4_e = dose3D(find(yVals == yDepth4),(xIndexStart:xIndexEnd),int16(s))';

%%Need to pick the indices for dose profile at 5cm, 10cm, 20cm. And shift
%%the middle as x_coordinate = 0, as in the golden beam data.
%% Also need to scale the maximum dose at this depth according to the PDD
%% at this depth, because the profile dose is relative (max = 100).
%% Now only account the dose profile at depth 5cm.

%options = optimset('DISPLAY', 'notify', 'TolX', 1.e-12);
%%%% Using Ant Colony Opt. (api.m)
%%% O.K. parameters. JC

% Parameters: p(1) till p(4), the four parameters to determine the primary
% photon spectrum, sum(a) is the weight.
% p(5) "shrink" for the Flattening filter.
% p(6) is the weight of FF photon, aFF = aFF * p(6)
% p(7) is the Horn effect weight, respective to sum(a)
% p(8) is the electron contamination weight, respective to sum(a)

FUN='PDD_Profile_res';
% "working" LB, UB settings
% Add one more bounds for the penumbra
%LB=[0.6  -0.1 1.2   100  2  0.1  0.01    0.005 0.01];
%UB=[1.5  0.1   10  1000  4   0.2  0.05  0.04 5];

% expand the LB and UB by the suggesstions from JOD. On Oct 09, 2006
% LB=[0.6  -0.1 1.2   100  1  0.1  0.005    0.001];
% UB=[1.5  0.1   10  1000  4   0.25  0.05  0.05];

NumAnts=25;Nmoves=20;LocalMoves=30;
NONLNCON=[];rpmax=[]; RES=1e-8;
p=api(FUN,LB,UB,NumAnts,Nmoves,LocalMoves,NONLNCON,rpmax,RES, ...
    doseV_obj, doseV_array, doseV_arrayFF, doseV_e, ...
    doseProfile1_obj, doseP1_array, doseP1_arrayFF, doseP1_e, ...
    doseProfile2_obj, doseP2_array, doseP2_arrayFF, doseP2_e, ...
    doseProfile3_obj, doseP3_array, doseP3_arrayFF, doseP3_e, ...
    doseProfile4_obj, doseP4_array, doseP4_arrayFF, doseP4_e, ...
    energy, numBin, extraBin, numBinFF)

% Output results
ener = energy/numBin: energy/numBin:energy*(1+extraBin/numBin);
a = Fatigue(p(1:4), ener);
a(find(isnan(a))) = 0;

cutoff = 0.85;
Ef = energy*cutoff;
kt = energy*0.15;  % change kt
f = 1./(1+exp((ener-Ef)/kt));
a = a.*f;

doseV_wt = (doseV_array)  * real(a)';
doseProfile1_wt = (doseP1_array) * real(a)';
doseProfile2_wt = (doseP2_array) * real(a)';
doseProfile3_wt = (doseP3_array)  * real(a)';
doseProfile4_wt = (doseP4_array) * real(a)';

% Get the weights/spectrum for the flattening filter
enerFF = energy/numBin : energy/numBin : numBinFF*energy/numBin;
aFF = Fatigue(p(1:4), p(5)*enerFF);
aFF(find(isnan(aFF))) = 0;

% Get modified Flattening filter spectrum/weights
f = 1./(1+exp((p(5)*enerFF-Ef)/kt));
aFF = aFF.* f;
aFF = aFF * (sum(a)/sum(aFF)*p(6)) ;  % scale FF weight respect to sum(a)

doseV_wtFF = doseV_arrayFF * real(aFF)';
doseProfile1_wtFF = doseP1_arrayFF * real(aFF)';
doseProfile2_wtFF = doseP2_arrayFF * real(aFF)';
doseProfile3_wtFF = doseP3_arrayFF * real(aFF)';
doseProfile4_wtFF = doseP4_arrayFF * real(aFF)';

filterSize = 7;
filterWindow = [-3 -2 -1 0 1 2 3];
G1 = gauss(filterWindow, p(9));  %p(9) is sigma
G1 = G1/sum(G1);   % normalize, to make sum(G1)==1;

doseV_DPM = doseV_wt + doseV_wtFF + p(8) * sum(a)* doseV_e;
doseProfile1_DPM = doseProfile1_wt + doseProfile1_wtFF + p(8) * sum(a)* doseP1_e;
doseProfile1_DPM = conv(G1, doseProfile1_DPM);
doseProfile2_DPM = doseProfile2_wt + doseProfile2_wtFF + p(8) * sum(a)* doseP2_e;
doseProfile2_DPM = conv(G1, doseProfile2_DPM);
doseProfile3_DPM = doseProfile3_wt + doseProfile3_wtFF + p(8) * sum(a)* doseP3_e;
doseProfile3_DPM = conv(G1, doseProfile3_DPM);
doseProfile4_DPM = doseProfile4_wt + doseProfile4_wtFF + p(8) * sum(a)* doseP4_e;
doseProfile4_DPM = conv(G1, doseProfile4_DPM);

figure;
subplot(2,1,1); plot(depth,dose, '.');
hold on;
plot(depthTPS,doseV_DPM, 'r');
legend('measured', 'DPM wt');
axis([0 25 0 1.1])
subplot(2,1,2); plot([0 ener], [0 a], '+-b')
hold on;  plot([0 enerFF], [0 aFF], '+-r')

% Note: There's shift in the measured dose. Now fix the display
figure;
plot(xTPS, doseProfile1_obj, 'b', xTPS, doseProfile1_DPM(4:end-3), '-.r')
hold on; plot(xTPS, doseProfile3_obj, 'c', xTPS, doseProfile3_DPM(4:end-3), '-.m')
legend('measured', ['DPM ', num2str(depth1)], 'measured', ['DPM ', num2str(depth3)]);

figure; plot(xTPS, doseProfile2_obj, 'b', xTPS, doseProfile2_DPM(4:end-3), '-.r')
hold on; plot(xTPS, doseProfile4_obj, 'c', xTPS, doseProfile4_DPM(4:end-3), '-.m');
legend('measured', ['DPM ', num2str(depth2)], 'measured', ['DPM ', num2str(depth4)]);

filename = 'modelParameters';
save(filename, 'p', 'ener', 'a', 'enerFF', 'aFF', 'max_doseV', 'doseV_wt');
disp('ok');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function result=isControlLike(modifiers)

result= ( ( ~ismac() && any(strcmp('control',modifiers)) ) ...
          || ...
          (  ismac() && any(strcmp('command',modifiers)) ) ) ;

end


text:
function newFormat = openBabelConverter(origFormat, outputFormat, saveFileDir)
% This function converts chemoformatic formats using OpenBabel. It requires
% to have openbabel installed. The formats that can be converted are used
% MDL MOL, SMILES, InChI, InChIKey, MDL RXN, reaction SMILES and rInChI.
%
% USAGE:
%
%    newStructure = openBabelConverter(origFormat, outputFormat, saveFile)
%
% INPUT:
%    origFormat:    Original chemoinformatic format. Chemical tables such
%                   as MDL MOL or MDL RXN must be provided as files
%    outputFormat:  The format to be converted. Formats supported: smiles,
%                   mol, inchi, inchikey, rxn and rinchi.
%
% OPTIONAL INPUTS:
%    saveFileDir:	String with the directory where the new format will be
%                   saved. If is empty, the format is not saved.
%
% EXAMPLE:
%
%    Example 1 (MDL MOL to InChI):
%    origFormat = [pwd filesep 'alanine.mol'];
%    outputFormat = 'inchi';
%    newFormat = openBabelConverter(origFormat, outputFormat);
%
%    Example 2 (InChI to SMILES):
%    origFormat = 'InChI=1S/C3H7NO2/c1-2(4)3(5)6/h2H,4H2,1H3,(H,5,6)/t2-/m0/s1';
%    outputFormat = 'smiles';
%    newFormat = openBabelConverter(origFormat, outputFormat);
%
%    Example 3 (SMILES to mol):
%    origFormat = 'C[C@@H](C(=O)O)N';
%    outputFormat = 'mol';
%    newFormat = openBabelConverter(origFormat, outputFormat);

if nargin < 3 || isempty(saveFileDir)
    toSave = false;
else
    toSave = true;
end

% Check openbabel installation
if isunix || ispc
    obabelCommand = 'obabel';
else
    obabelCommand = 'openbabel.obabel';
end
[oBabelInstalled, cmdout] = system(obabelCommand);
if oBabelInstalled == 127
    error('To use this function, Open Babel must be installed; follow the installation instructions in https://openbabel.org/wiki/Category:Installation')
end

% Identify the chemoinformatic format of the input
if contains(origFormat, '.mol')
    if isfile(origFormat)
        inputType = 'mol';
    else
        error(['The file ' origFormat ' is missing'])
    end
elseif contains(origFormat, '.rxn')
    if isfile(origFormat)
        inputType = 'rxn';
    else
        error(['The file ' origFormat ' is missing'])
    end
elseif contains(origFormat, 'InChI=')
    inputType = 'inchi';
else
    inputType = 'smiles';
end

% Convert
if ismember(inputType, {'inchi'; 'smiles'})

    fid2 = fopen('tmp', 'w');
    fprintf(fid2, '%s\n', origFormat);
    fclose(fid2);
    [~, cmdout] = system([obabelCommand ' -i' inputType ' tmp -o' outputFormat]);

else
    switch inputType
        case 'mol'
            [~, cmdout] = system([obabelCommand ' -imol ' origFormat ' -o' outputFormat]);
        case 'rxn'
            [~, cmdout] = system([obabelCommand ' -irxn ' origFormat ' -o' outputFormat]);
    end
end

% Prepare the output
cmdout = splitlines(cmdout);
newFormat = [];
switch outputFormat
    case 'mol'
        startIdx = find(cellfun(@isempty, cmdout));
        endIdx = find(ismember(cmdout, 'M  END'));
        newFormat = cmdout(startIdx(1):endIdx);

    case 'inchikey'
        cmdout = split(cmdout{end - 2});
        newFormat = cmdout{1};

    case 'inchi'
        if any(contains(cmdout,'InChI=1S'))
            newFormat = cmdout{contains(cmdout,'InChI=1S')};
        else
            error('')
        end

    case 'smiles'
        cmdout = split(cmdout{end - 2});
        newFormat = cmdout{1};

    case 'rinchi'
        newFormat = cmdout{contains(cmdout,'RInChI=')};
end

% Save the file
if toSave && ~isempty(newFormat)
    fid2 = fopen(saveFileDir, 'w');
    fprintf(fid2, '%s\n', newFormat{:});
    fclose(fid2);
% elseif isempty(newFormat)
%     error(['The format ' origFormat ' couldn''t be converted'])
end

author: opencobra

repo: cobratoolbox


text:
function test_bug2597

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

% this function tests the functionality of ft_selectdata to preserve the
% ordering of the channels according to the first input argument (thus
% overruling alphabetization, as well as any order in the cfg).

% make some data
data        = [];
data.label  = {'b';'a';'c'};
data.avg    = repmat([2 1 3]',[1 2]);
data.dimord = 'chan_time';
data.time   = [1 2];

datanew = ft_selectdata([], data);
assert(isequal(datanew.label,{'b';'a';'c'}));
assert(isequal(datanew.avg(:,1),[2 1 3]'));

data2        = [];
data2.label  = {'a';'b';'c'};
data2.avg    = repmat([1 2 3]',[1 2]);
data2.dimord = 'chan_time';
data2.time   = [1 2];

[datanew, datanew2] = ft_selectdata([], data, data2);
assert(isequal(datanew.label,{'b';'a';'c'}));
assert(isequal(datanew2.label,datanew.label));

data2        = [];
data2.label  = {'b';'c';'d'};
data2.avg    = repmat([2 3 4]',[1 2]);
data2.dimord = 'chan_time';
data2.time   = [1 2];
[datanew, datanew2] = ft_selectdata([], data, data2);
assert(isequal(datanew.label,{'b';'c'}));
assert(isequal(datanew2.label,{'b';'c'}));
assert(isequal(datanew.avg(:,1),[2 3]'));
assert(isequal(datanew2.avg(:,1),[2 3]'));

cfg.select = 'union';
[datanew, datanew2] = ft_selectdata(cfg, data, data2);
assert(isequal(datanew.label,{'b';'a';'c';'d'}));
assert(isequal(datanew2.label,{'b';'a';'c';'d'}));
assert(isequaln(datanew.avg(:,1),[2 1 3 nan]'));
assert(isequaln(datanew2.avg(:,1),[2 nan 3 4]'));

data2        = [];
data2.label  = {'b';'c';'e'};
data2.avg    = repmat([2 3 5]',[1 2]);
data2.dimord = 'chan_time';
data2.time   = [1 2];

data3        = [];
data3.label  = {'d';'a';'c'};
data3.avg    = repmat([4 1 3]',[1 2]);
data3.dimord = 'chan_time';
data3.time   = [1 2];

cfg.select = 'intersect';
[datanew, datanew2, datanew3] = ft_selectdata(cfg, data, data2, data3);
assert(isequal(datanew.label, {'c'}));
assert(isequal(datanew2.label, {'c'}));
assert(isequal(datanew3.label, {'c'}));
assert(isequaln(datanew.avg(:,1), 3));
assert(isequaln(datanew2.avg(:,1), 3));
assert(isequaln(datanew3.avg(:,1), 3));

cfg.select = 'union';
[datanew, datanew2, datanew3] = ft_selectdata(cfg, data, data2, data3);
assert(isequal(datanew.label, {'b';'a';'c';'e';'d'}));
assert(isequal(datanew2.label, {'b';'a';'c';'e';'d'}));
assert(isequal(datanew3.label, {'b';'a';'c';'e';'d'}));
assert(isequaln(datanew.avg(:,1), [2 1 3 nan nan]'));
assert(isequaln(datanew2.avg(:,1), [2 nan 3 5 nan]'));
assert(isequaln(datanew3.avg(:,1), [nan 1 3 nan 4]'));