text:
% Test file for trigtech/roots.m

function pass = test_roots(pref)

if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

%% Simple test:
f = testclass.make(@(x) cos(5*pi*x), [], pref);
r = roots(f);
exact = (-0.9:0.2:0.9).';
pass(1) = norm(r-exact,Inf) < 1e1*length(f)*eps;

%% More complicated:
k = 20;
f = testclass.make(@(x) sin(sin(pi*k*x)), [], pref);
r = roots(f);
pass(2) = norm(r-(-k:k)'/k, inf) < length(f)*eps;

%% No real roots.
f = testclass.make( @(x) 3./(5 - 4*cos(3*pi*x)), ...
    [], pref);
r = roots(f);
pass(3) = isempty(r);

%% Test that complex roots are now allowed.
f = testclass.make(@(x) 2 + cos(pi*x), [], pref);
r = roots(f, 'complex', 1);
pass(4) = norm( feval(f,r) ) < vscale(f)*eps;

f = testclass.make(@(x) sin(100*pi*x));
r1 = roots(f, 'complex', 1);
r2 = roots(f);
pass(5) = numel(r1) == 200 & numel(r2) >= 201;

%% Test an array-valued function:
f = testclass.make(@(x) [sin(pi*x), cos(pi*x)], [], pref);
r = roots(f);
r2 = [-1 0 1 -.5 .5 NaN].';
pass(6) = all( r(:) - r2 < 10*length(f)*eps | isnan(r2) );

f = testclass.make(@(x) [cos(2*pi*x), sin(pi*x)], [], pref);
r = roots(f, 'complex', 1);
r = r(:);
[temp, id] = sort(real(r));
r = r(id);
r2 = [-0.75 -0.25 0 0.25 0.75 1 NaN NaN].';
pass(7) = all( abs(r(:) - sort(r2)) < 1e1*length(f)*eps | isnan(r2) );

end


text:
function[varargout] = vmean(varargin)
%VMEAN  Mean over non-NaN elements along a specified dimension.
%
%   Y=VMEAN(X,DIM) takes the mean of all non-NaN elements of X along
%   dimension DIM.
%
%   [Y,NUM]=VMEAN(X,DIM) also outputs the number of non-NaN data points
%   NUM, which has the same dimension as X.
%
%   [Y1,Y2,...YN]=VMEAN(X1,X2,...XN,DIM) or also works, where all the XN
%   are the same size.
%
%   VMEAN(X1,X2,...XN,DIM);  with no output arguments overwrites the
%   original input variables.
%   __________________________________________________________________
%
%   Weighted means
%
%   VMEAN can also form a weighted mean.
%
%   Y=VMEAN(X,DIM,W) forms the weighted mean of X using weights W, an array
%   of the same size as X.
%
%   In this case [Y,NUM]=VMEAN(X,DIM,W) returns total weight used in
%   forming the mean at each point, rather than the number of data points.
%
%   Thus VMEAN(X,DIM,W) is the same as VMEAN(X.*W,DIM)./VMEAN(W,DIM).
%
%   [Y1,Y2,...YN]=VMEAN(X1,X2,...XN,DIM,W) also works.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2001--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  vmean_test;powermean_test;return
end

if ~isscalar(varargin{end})
    w=varargin{end};
    varargin=varargin(1:end-1);
else
    w=[];
end

dim=varargin{end};

for i=1:length(varargin)-1
  if isempty(w)
      [varargout{i},numi{i}]=vsum(varargin{i},dim);
      numi{i}(numi{i}==0)=nan;
      varargout{i}=varargout{i}./numi{i};
  else
      varargout{i}=vsum(varargin{i}.*w,dim);
      numi{i}=vsum(w,dim);
      varargout{i}=varargout{i}./numi{i};
  end
end

for i=length(varargin):nargout
  varargout{i}=numi{i-length(varargin)+1};
end

eval(to_overwrite(nargin-1))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[]=vmean_test
x1=[1 2 ; nan 4];
x2=[inf 6; nan 5];
ans1=[3/2 4]';
ans2=[6 5]';

vmean(x1,x2,2);
reporttest('VMEAN output overwrite', aresame(x1,ans1) && aresame(x2,ans2))

x1=[1 2 ; nan 4];
ans1=[3/2 4]';
ans2=[2 1]';

[y1,y2]=vmean(x1,2);
reporttest('VMEAN mean & num', aresame(y1,ans1) && aresame(y2,ans2))

function[varargout]=powermean(varargin)
%POWERMEAN  Power-weighted mean along a specified dimension.
%
%   FM=POWERMEAN(F,X,DIM) takes the mean of all finite elements of F along
%   dimension DIM, weighted by the squared magnitude of X.
%
%   The power-weighted mean of F is defined as
%
%        FM = SUM (ABS(X)^2.*F, DIM) / SUM(ABS(X)^2, DIM)
%
%   where X and F are arrays of the same size.
%
%   [FM1,FM2,...FMN]=POWERMEAN(F1,F2,...FN,X,DIM) also works.
%
%   POWERMEAN(F1,F2,...FN,X,DIM);  with no output arguments overwrites the
%   original input variables.
%
%   POWERMEAN with X a set of analytic signals is used to construct the
%   joint instantaneous moments.  For details on these quantities, see
%
%       Lilly and Olhede (2010),  Bivariate instantaneous frequency and
%           bandwidth.  IEEE Trans. Sig. Proc. 58 (2), 591--603.
%
%   See also VMEAN.
%
%   Usage: fm=powermean(f,x,dim);
%          [fm1,fm2]=powermean(f1,f2,x,dim);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2009--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    powermean_test,return
end

dim=varargin{end};
x=varargin{end-1};

power=vsum(abs(x).^2,dim);
vswap(power,0,nan);
for i=1:length(varargin)-2
  varargout{i}=vsum(abs(x).^2.*varargin{i},dim)./power;
end

eval(to_overwrite(nargin-2))

function[]=powermean_test

bool(1)=aresame(powermean([1 2],[2 3],2),frac(4+2*9,13));
bool(2)=aresame(powermean([1 2],[2 3],2),vmean([1 2],2,[2 3].^2));

reporttest('VMEAN weighted mean',allall(bool))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function partest(varargin)
%This function calculate the performance, based on Bayes theorem, of a
%clinical test.
%
% Syntax: 	PARTEST(X,ALPHA)
%
%     Input:
%           X is the following 2x2 matrix.
%           ALPHA - significance level for confidence intervals (default = 0.05).
%
%....................Affected(D+)..Healthy(D-)
%                    _______________________
%Positive Test(T+)  |   True    |  False    |
%                   | positives | positives |
%                   |___________|___________|
%                   |  False    |   True    |
%Negative Test(T-)  | negatives | negatives |
%                   |___________|___________|
%
%     Outputs:
%           - Prevalence
%           - Sensibility
%           - Specificity
%           - False positive and negative proportions
%           - False discovery and discovery rates
%           - Youden's Index and Number Needed to Diagnose (NDD)
%           - Positive predictivity
%           - Positive Likelihood Ratio
%           - Negative predictivity
%           - Negative Likelihood Ratio
%           - Predictive Summary Index (PSI) and Number Needed to Screen (NNS)
%           - Test Accuracy
%           - Mis-classification Rate
%           - F-Measure
%           - Test bias
%           - Error odds ratio
%           - Diagnostic odds ratio
%           - Discriminant Power
%
%      Example:
%
%
%           Calling on Matlab the function: partest
%           it will use a default matrix x=[731 270;78 1500] and alpha=0.05
%
%           Answer is:
%
%DIAGNOSTIC TEST PERFORMANCE PARAMETERS
% ----------------------------------------------------------------------------------------------------
% Prevalence: 31.4% (29.6% - 33.2%)
%
% Sensitivity (probability that test is positive on unhealthy subject): 90.4% (89.1% - 91.5%)
% False negative proportion: 9.6% (8.5% - 10.9%)
% False discovery rate: 27.0% (25.3% - 28.7%)
%
% Specificity (probability that test is negative on healthy subject): 84.7% (83.3% - 86.1%)
% False positive proportion: 15.3% (13.9% - 16.7%)
% False omission rate: 4.9% (4.2% - 5.9%)
%
% Youden's Index (a perfect test would have a Youden's index of +1): 0.7510
% Number Needed to Diagnose (NND): 1.33
% Around 14 persons need to be tested to return 10 positive tests for the presence of disease
%
% Precision or Predictivity of positive test
% (probability that a subject is unhealthy when test is positive): 73.0% (71.3% - 74.7%)
% Positive Likelihood Ratio: 5.9 (5.7 - 6.2)
% Moderate increase in possibility of disease presence
%
% Predictivity of negative test
% (probability that a subject is healthy when test is negative): 95.1% (94.1% - 95.8%)
% Negative Likelihood Ratio: 0.1138 (0.1094 - 0.1183)
% Moderate increase in possibility of disease absence
%
% Predictive Summary Index: 0.6808
% Number Needed to Screen (NNS): 1.47
% Around 15 persons need to be screened to avoid 10 events (i.e. death) for the presence of disease
%
% Accuracy or Potency: 86.5% (85.1% - 87.8%)
% Mis-classification Rate: 13.5% (12.2% - 14.9%)
% F-measure: 80.8% (79.2% - 82.3%)
%
% Test bias: 1.2373 (0.9474 - 1.6160)
% Test overestimates the phenomenon
% Error odds ratio: 1.6869 (1.2916 - 2.2032)
% Diagnostic odds ratio: 52.0655 (39.8649 - 68.0002)
% Discriminant Power: 2.2
%      A test with a discriminant value of 1 is not effective in discriminating between affected and unaffected individuals.
%      A test with a discriminant value of 3 is effective in discriminating between affected and unaffected individuals.
%
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2006). Clinical test performance: the performance of a
% clinical test based on the Bayes theorem.
% http://www.mathworks.com/matlabcentral/fileexchange/12705

%Input Error handling
args=cell(varargin);
nu=numel(args);
if nu>2
    error('Warning: Max two input data are required')
end
default.values = {[731 270;78 1500],0.05};
default.values(1:nu) = args;
[x alpha] = deal(default.values{:});

if nu>=1
    if isvector(x)
        error('Warning: PARTEST requires a 2x2 input matrix')
    end
    if ~all(isfinite(x(:))) || ~all(isnumeric(x(:))) || ~isequal(x,round(x))
        error('Warning: all matrix values must be numeric integer and finite')
    end
    [r,c] = size(x);
    if (r ~= 2 || c ~= 2)
        error('Warning: PARTEST requires a 2x2 input matrix')
    end
    clear r c
end
if nu>1
    if ~isscalar(alpha) || ~isnumeric(alpha) || ~isfinite(alpha) || isempty(alpha)
        error('Warning: it is required a numeric, finite and scalar ALPHA value.');
    end
    if alpha <= 0 || alpha >= 1 %check if alpha is between 0 and 1
        error('Warning: ALPHA must be comprised between 0 and 1.')
    end
end
clear args default nu

global z N

z=-realsqrt(2)*erfcinv(2-alpha);
cs=sum(x); %columns sum
rs=sum(x,2); %rows sums
N=sum(x(:)); %numbers of elements
d=diag(x); %true positives and true negatives

clc
disp('DIAGNOSTIC TEST PERFORMANCE PARAMETERS')
disp(repmat('-',1,100))
%Prevalence
%the prevalence of disease is estimated all D+/N
pr=(cs(1)/N);
%95% confidence interval critical value for Prevalence
prci=newcombe(pr);
fprintf('Prevalence: %0.1f%% (%0.1f%% - %0.1f%%)\n',pr*100,prci(1).*100,prci(2).*100)
disp(' ')

%Sensitivity and Specificity
%The Sensitivity is the probability that the test is positive on sick subjects: P(T+|D+)
%The Specificity is the probability that the test is negative on healthy subjects: P(T-|D-)
%In Matlab both parameters are obtained with only one instruction:
SS=d./cs'; %Sensitivity and Specificity
%Of course the false proportion is the complement to 1
fp=1-SS; %false proportions
%The false discovery rate is the probability that the disease is absent when the test is positive: P(D-|T+)
%The false omission rate is the probability that the disease is present when the test is negative: P(D+|T-)
fd=diag(rot90(x))./rs; %false discovery and omission rate

% 95% confidence intervals
Seci=newcombe(SS(1)); Spci=newcombe(SS(2));
fnci=newcombe(fp(1)); fpci=newcombe(fp(2));
fdci=newcombe(fd(1)); foci=newcombe(fd(2));
fprintf('Sensitivity (probability that test is positive on unhealthy subject): %0.1f%% (%0.1f%% - %0.1f%%)\n',SS(1)*100,Seci(1).*100,Seci(2).*100)
fprintf('False negative proportion: %0.1f%% (%0.1f%% - %0.1f%%)\n',fp(1)*100,fnci(1).*100,fnci(2).*100)
fprintf('False discovery rate: %0.1f%% (%0.1f%% - %0.1f%%)\n',fd(1)*100,fdci(1).*100,fdci(2).*100)
disp(' ')
fprintf('Specificity (probability that test is negative on healthy subject): %0.1f%% (%0.1f%% - %0.1f%%)\n',SS(2)*100,Spci(1).*100,Spci(2).*100)
fprintf('False positive proportion: %0.1f%% (%0.1f%% - %0.1f%%)\n',fp(2)*100,fpci(1).*100,fpci(2).*100)
fprintf('False omission rate: %0.1f%% (%0.1f%% - %0.1f%%)\n',fd(2)*100,foci(1).*100,foci(2).*100)
disp(' ')

%Youden's Index
%Youden's J statistics (also called Youden's index) is a single statistic that
%captures the performance of a diagnostic test. The use of such a single index
%is "not generally to be recommended". It is equal to the risk difference for a
%dichotomous test and it defined as: J = Sensitivity + Specificity - 1.
%A perfect test has J=1.
J=sum(SS)-1; %Youden's index
fprintf('Youden''s Index (a perfect test would have a Youden''s index of +1): %0.4f\n', J)

%The number needed to diagnose is defined as the number of patients that
%need to be tested to give one correct positive test.
NND=1/J; %Number needed to Diagnose (NND)
fprintf('Number Needed to Diagnose (NND): %0.2f\n',NND);
fprintf('Around %i persons need to be tested to return 10 positive tests for the presence of disease\n',ceil(NND*10))
disp(' ')

%Positive and Negative predictivity
%Positive predictivity is the probability that a subject is sick when test is positive: P(D+|T+)
%Negative predictivity is the probability that a subject is healthy when test is negative: P(D-|T-)
%Positive predictivity=Precision
%In Matlab both parameters are obtained with only one instruction:
PNp=d./rs; %Positive and Negative predictivity
% 95% confidence interval critical value for Positive and Negative predictivity
Ppci=newcombe(PNp(1));
Npci=newcombe(PNp(2));
%Positive and Negative Likelihood Ratio
%When we decide to order a diagnostic test, we want to know which test (or
%tests) will best help us rule-in or rule-out disease in our patient.  In the
%language of clinical epidemiology, we take our initial assessment of the
%likelihood of disease (“pre-test probability”), do a test to help us shift our
%suspicion one way or the other, and then determine a final assessment of the
%likelihood of disease (“post-test probability”).
%Likelihood ratios tell us how much we should shift our suspicion for a
%particular test result. Because tests can be positive or negative, there are at
%least two likelihood ratios for each test.  The “positive likelihood ratio”
%(LR+) tells us how much to increase the probability of disease if the test is
%positive, while the “negative likelihood ratio” (LR-) tells us how much to
%decrease it if the test is negative.
%You can also define the LR+ and LR- in terms of sensitivity and specificity:
%LR+ = sensitivity / (1-specificity)
%LR- = (1-sensitivity) / specificity
plr=SS(1)/fp(2); %Positive Likelihood Ratio
plrci=LRci(plr,x);
nlr=fp(1)/SS(2); %Negative Likelihood Ratio
nlrci=LRci(nlr,x);

fprintf('Precision or Predictivity of positive test\n')
fprintf('(probability that a subject is unhealthy when test is positive): %0.1f%% (%0.1f%% - %0.1f%%)\n', PNp(1)*100,Ppci(1).*100,Ppci(2).*100)
fprintf('Positive Likelihood Ratio: %0.1f (%0.1f - %0.1f)\n',plr,plrci(1),plrci(2))
dlr(plr)
disp(' ')
fprintf('Predictivity of negative test\n')
fprintf('(probability that a subject is healthy when test is negative): %0.1f%% (%0.1f%% - %0.1f%%)\n', PNp(2)*100,Npci(1).*100,Npci(2).*100)
if nlr<1e-4 || nlrci(1)<1e-4
    fprintf('Negative Likelihood Ratio: %0.4e (%0.4e - %0.4e)\n',nlr,nlrci(1),nlrci(2))
else
    fprintf('Negative Likelihood Ratio: %0.4f (%0.4f - %0.4f)\n',nlr,nlrci(1),nlrci(2))
end
dlr(nlr)
disp(' ')

%Predictive Summary Index (similar to Youden's Index)
PSI=sum(PNp)-1;
NNS=1/PSI; %Number needed to screen
fprintf('Predictive Summary Index: %0.4f\n', PSI)
fprintf('Number Needed to Screen (NNS): %0.2f\n',NNS);
fprintf('Around %i persons need to be screened to avoid 10 events (i.e. death) for the presence of disease\n',ceil(NNS*10))
disp(' ')

%Accuracy and Mis-classification rate
%Diagnostic accuracy (or Power) is defined as the proportion of all tests
%that give a correct result. The Mis-classification rate is its complement to 1.
%In statistics, the F1 score (also F-score or F-measure) is a measure of a
%test's accuracy. It considers both the Precision (positive predictivity)
%and the Sensitivity of the test to compute the score:
%P is the number of correct results divided by the number of all returned results
%S is the number of correct results divided by the number of results that should
%have been returned.
%The F1 score can be interpreted as a weighted average of the Precision and
%Sensitivity, where an F1 score reaches its best value at 1 and worst score at 0.
acc=trace(x)/N; %Accuracy
accci=newcombe(acc);
mcr=1-acc; %Mis-classification rate
mcrci=newcombe(mcr);
FMeasure=harmmean([SS(1) PNp(1)]); %F-measure
FMci=newcombe(FMeasure);
fprintf('Accuracy or Potency: %0.1f%% (%0.1f%% - %0.1f%%)\n',acc*100,accci(1).*100,accci(2).*100)
fprintf('Mis-classification Rate: %0.1f%% (%0.1f%% - %0.1f%%)\n',mcr*100,mcrci(1).*100,mcrci(2).*100)
fprintf('F-measure: %0.1f%% (%0.1f%% - %0.1f%%)\n',FMeasure*100,FMci(1).*100,FMci(2).*100)
disp(' ')

orse=realsqrt(sum(1./x(:))); %standard error of log(OR)
cv=([-1 1].*(z*orse));

%Test Bias (TB)
%A test which shows provable and systematic differences in the results of people
%based on group membership. For example, a test might be considered biased if
%members of one particular gender or race consistently and systematic have
%statistically different results from the rest of the testing population.
%It is defined as (T+)/(D+)=(TP+FP)/(TP+FN)
%A perfect test has a TB=1;
%If TB<1 the test underestimates the disease because there are more affected than positive test
%If TB>1 the test overestimates the disease because there are more positive test than affected
TB=rs(1)/cs(1); %Test Bias
orci=exp(reallog(TB)+cv); %OR confidence interval
fprintf('Test bias: %0.4f (%0.4f - %0.4f)\n',TB,orci(1),orci(2))
if TB>1
    disp('Test overestimates the phenomenon')
elseif TB<1
    disp('Test underestimates the phenomenon')
else
    disp('There is not test bias')
end

%Error Odds Ratio.
%Indicates if the probability of being wrongly classified is highest in the
%diseased or in the non-diseased group. If the error odds is higher than one the
%probability is highest in the diseased group (and the specificity of the test
%is better than the sensitivity), if the value is lower than one the probability
%of an incorrect classification is highest in the non-diseased group (and the
%sensitivity of the test is better than the specificity).
%It is defined as (Sensitivity/(1-Sensitivity))/(Specificity/(1-Specificity));
EOR=(SS(1)/fp(1))/(SS(2)/fp(2)); %Error odds ratio
orci=exp(reallog(EOR)+cv); %OR confidence interval
fprintf('Error odds ratio: %0.4f (%0.4f - %0.4f)\n',EOR,orci(1),orci(2))

%Diagnostic Odds Ratio.
%Diagnostic odds ratio is defined as how much more likely will the test
%make a correct diagnosis than an incorrect diagnosis in patients with the
%disease (Scott et al. 2008).
%Often used as a measure of the discriminative power of the test. Has the value
%one if the test does not discriminate between diseased and not diseased. Very
%high values above one means that a test discriminates well. Values lower than
%one mean that there is something wrong in the application of the test.
%It is defined as (Sensitivity/(1-Sensitivity))/((1-Specificity)/Specificity);
DOR=(SS(1)/fp(1))/(fp(2)/SS(2)); %Diagnostic odds ratio
orci=exp(reallog(DOR)+cv); %OR confidence interval
fprintf('Diagnostic odds ratio: %0.4f (%0.4f - %0.4f)\n',DOR,orci(1),orci(2))

%Discriminant power
%The discriminant power for a test, also termed the test effectiveness, is a
%measure of how well a test distinguishes between affected and unaffected
%persons. It is the sum of logs of Sensivity and Specificity over own false
%proportion, scaled by the standard deviation of the logistic normal
%distribution curve (square root of 3 divided by π). Test effectiveness is
%interpreted as the standardized distance between the means for both populations.
%A test with a discriminant value of 1 is not effective in discriminating between affected and unaffected individuals.
%A test with a discriminant value of 3 is effective in discriminating between affected and unaffected individuals.
dpwr=(realsqrt(3)/pi)*sum(log(SS./fp)); %Discriminant power
fprintf('Discriminant Power: %0.1f\n',dpwr)
disp([blanks(5) 'A test with a discriminant value of 1 is not effective in discriminating between affected and unaffected individuals.'])
disp([blanks(5) 'A test with a discriminant value of 3 is effective in discriminating between affected and unaffected individuals.'])

%Display graphs
xg=cs./N;
figure
hold on
H=zeros(1,4);
H(1)=fill(xg(1)+[0 xg(2) xg(2) 0],SS(2)+[0 0 fp(2) fp(2)],'r');
H(2)=fill([0 xg(1) xg(1) 0],fp(1)+[0 0 SS(1) SS(1)],'g');
H(3)=fill([0 xg(1) xg(1) 0],[0 0 fp(1) fp(1)],'y');
H(4)=fill(xg(1)+[0 xg(2) xg(2) 0],[0 0 SS(2) SS(2)],'b');
hold off
axis square
title('PARTEST GRAPH')
xlabel('Subjects proportion')
ylabel('Parameters proportion')
legend(H,'False Positive','True Positive (Sensibility)','False Negative','True Negative (Specificity)','Location','NorthEastOutside')

%The rose plot is a variation of the common pie chart. For both, we have k data
%points where each point denotes a frequency or a count. Pie charts and rose
%plots both use the area of segments of a circle to convey amounts.
%The pie chart uses a common radius and varies the central angle according to
%the data. That is, the angle is proportional to the frequency. So if the i-th
%point has count X and the total count is N, the i-th angle is 360*(X/N).
%For the rose plot, the angle is constant (i.e., divide 360 by the number of
%groups, k) and it is the square root of the radius that is proportional to the
%data.
%According to Wainer (Wainer (1997), Visual Revelations: Graphical Tales of Fate
%and Deception from Napolean Bonaporte to Ross Perot, Copernicus, Chapter 11.),
%the use of a common angle is the strength of the rose plot since it allows us
%to easily compare a sequence of rose plots (i.e., the corresponding segments in
%different rose plots are always in the same relative position).
%In particular, this makes rose plots an effective technique for displaying the
%data in contingency tables.
%As an interesting historical note, Wainer points out that rose plots were used
%by Florence Nightingale (she referred to them as coxcombs).
%color={'r','g','y','b'};
H=roseplot([fp(2) SS(1) fp(1) SS(2)]);
title('ROSEPLOT PARTEST GRAPH')
legend(H,'False Positive','True Positive (Sensibility)','False Negative','True Negative (Specificity)','Location','NorthEastOutside')
return
end

function dlr(lr) %Likelihood dialog
if lr==1
    disp('Test is not suggestive of the presence/absence of disease')
    return
end

if lr>10 || lr<0.1
    p1='Large (often conclusive)';
elseif (lr>5 && lr<=10) || (lr>0.1 && lr<=0.2)
    p1='Moderate';
elseif (lr>2 && lr<=5) || (lr>0.2 && lr<=0.5)
    p1='Low';
elseif (lr>1 && lr<=2) || (lr>0.5 && lr<=1)
    p1='Poor';
end

p2=' increase in possibility of disease ';

if lr>1
    p3='presence';
elseif lr<1
    p3='absence';
end
disp([p1 p2 p3])
return
end

function H=roseplot(x)
color={'r','g','y','b'};
k=length(x); H=zeros(1,k);
ang=(0:1:k)./k.*2.*pi;
figure
axis square
hold on
for I=1:k
    theta=[ang(I) linspace(ang(I),ang(I+1),500) ang(I+1)];
    rho=[0 repmat(realsqrt(x(I)),1,500) 0];
    [xg,yg]=pol2cart(theta,rho);
    H(I)=patch(xg,yg,color{I});
end
hold off
return
end

function ci=newcombe(p)
global z N
a=2*N*p+z^2;
b=z*realsqrt((z^2-2-1/N+4*p*(N*(1-p)+1)));
c=(2*(N+z^2));
%Of course, the critical interval lower bound cannot be less than 0 and the
%upper bound cant be greater than 1 and so:
ci(1)=max([0 (a-b-1)/c]);
ci(2)=min([1 (a+b+1)/c]);
return
end

function ci=LRci(lr,x)
global z
a=log(lr); b=realsqrt(sum(1./diag(x))-sum(1./sum(x,2)))*z;
ci=exp([a-b a+b]);
return
end

author: Sable

repo: mcbench-benchmarks


text:
function H= grid_addBars(fv, varargin)
%GRID_ADDBARS - Add a Colorbar to the Bottom of Each Subplot of a Grid Plot
%
%Synopsis:
% H= grid_addBars(FV, <opt>)
%
%Input:
% FV: data structure (like epo)
% OPT: struct or property/value list of optional properties:
%  .Height    - Height of the colorbar
%  .HScale    - Handle to the scale of the grid plot (returned by
%               grid_plot as .scale), a colorbar will be placed next to
%               the scale
%  .ShiftAxes - 'position' or 'ylim'
%  .CLim      - Color limit of colorbar, either a 1-by-2 vector giving
%               the limit or 'auto' (default)
%  .Colormap  - The colormap of colorbar, either a x-by-3 vector giving
%               the explicit map, or one of the predefined maps (see
%               help of 'colormap')
%  .Channels  - Channels for which bar should be added, either a cell
%               of strings giving the channels names, or en empty array
%               for all channels, or 'plus'
%  .Scale*    - with * in {'Height','Width','VPos','LeftShift',
%               'FontSize','Digits','Unit'} define the appearance of
%               the scale
%
%Output:
%   H - handle to graphic objects
%
%Caveat:
% You have to make sure that the FV data fits to the data that was displayed
% with grid_plot before (e.g. for the length of the vector in the time
% dimension). The channels are matched automatically. FV may have less
% channels than displayed by grid_plot.
%
%Example:
% H= grid_plot(epo, mnt, defopt_erps)
% epo_rsq= proc_r_square_signed(epo);
% grid_addBars(epo_rsq, 'HScale',H.scale)

props= {'VPos',                 0,                      'DOUBLE';
        'Height',               1/15,                   'DOUBLE';
        'ShiftAxes',            'ylim',                 'CHAR';
        'ShiftAlso',            {'scale'},              'CELL{CHAR}';
        'CLim',                 'auto',                 'CHAR|DOUBLE[2]';
        'Colormap',             flipud(gray(32)),       'DOUBLE[- 3]';
        'UseLocalColormap',     1,                      'BOOL';
        'MoveBack',             0,                      'BOOL';
        'Box',                  'on',                   'CHAR';
        'Visible',              'off',                  'CHAR';
        'AlphaStepsMode',       0,                      'BOOL';
        'HScale',               [],                     'STRUCT(ax)';
        'ScaleHeight',          .66,                    'DOUBLE';
        'ScaleWidth',           .075,                   'DOUBLE';
        'ScaleVPos',            .25,                    'DOUBLE';
        'ScaleLeftShift',       .05,                    'DOUBLE';
        'ScaleFontSize',        get(gca,'FontSize'),    'DOUBLE';
        'ScaleDigits',          4                       'DOUBLE';
        'ScaleUnit',            ''                      'CHAR';
        'Channels',             '*'                     'CHAR|DOUBLE'};

if nargin==0,
  H= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if ndims(fv.x)>2
  error('only one Class allowed');
end

if min(fv.x(:))<0,
  if isdefault.CLim,
    opt.CLim= 'sym';
  end
end
if isdefault.Colormap && ...
        (isequal(opt.CLim,'sym') || (isnumeric(opt.CLim) && opt.CLim(1)<0)),
  opt.Colormap= cmap_posneg(21);
end

if isdefault.AlphaStepsMode && ...
      isfield(fv, 'ClassName') && ...
      ~isempty(findstr(fv.ClassName{1},'alpha-steps')),
  opt.AlphaStepsMode= 1;
end
if opt.AlphaStepsMode,
  nValues= length(fv.alpha);
  [opt, isdefault]= ...
      opt_overrideIfDefault(opt, isdefault, ...
                            'CLim', [0 nValues], ...
                            'Colormap', flipud(gray(nValues+1)), ...
                            'ScaleUnit', '%');
end
if isdefault.ScaleUnit && isfield(fv, 'yUnit'),
  opt.ScaleUnit= fv.yUnit;
end
if isdefault.ScaleVPos && isempty(opt.ScaleUnit),
  opt.ScaleVPos= 0.5;
end
if isdefault.Visible && strcmpi(opt.Box, 'on'),
  opt.Visible= 'on';
end
if isdefault.UseLocalColormap,
  if visutil_isColormapUsed && ~isequal(opt.Colormap, get(gcf, 'colormap')),
    opt.UseLocalColormap= 1;
  end
end

[AxesStyle, lineStyle]= opt_extractPlotStyles(opt);

if opt.UseLocalColormap,
  iswhite= find(all(opt.Colormap==1,2));
  if ~isempty(iswhite),
    opt.Colormap(iswhite,:)= 0.9999*[1 1 1];
  end
  H.image= 'sorry: fakeimage - no handle';
else
  set(gcf, 'colormap',opt.Colormap);
end

if isnumeric(opt.Channels)
  ax= opt.Channels;
else
  ax= gridutil_getSubplots(opt.Channels);
end

% For image_localColormap we have to determine the CLim in advance.
% Therefore we need to determine the depicted channels.
% (If only image was used, this functions would be simpler.)
clab= cell(1, length(ax));
for ii= 1:length(ax),
  clab{ii}= getfield(get(ax(ii), 'userData'), 'chan');
  if iscell(clab{ii}),
    clab{ii}= clab{ii}{1};  % for multiple channels per ax, choose only
                            % the first one
  end
end
if strcmpi(opt.CLim, 'auto'),
  ci= util_chanind(fv, clab);
  mi= min(min(fv.x(:,ci)));
  if mi>=0 && isdefault.CLim,
    warning('know-it-all: switching to CLim mode ''0tomax''');
    opt.CLim= '0tomax';
  else
    opt.CLim= [mi max(max(fv.x(:,ci)))];
  end
elseif strcmpi(opt.CLim, 'sym'),
  ci= util_chanind(fv, clab);
  mi= min(min(fv.x(:,ci)));
  ma= max(max(fv.x(:,ci)));
  mm= max(abs(mi), ma);
  opt.CLim= [-mm mm];
end
if strcmpi(opt.CLim, '0tomax'),
  ci= util_chanind(fv, clab);
  opt.CLim= [0 max(max(fv.x(:,ci)))];
end

jj= 0;
for ii= 1:length(ax),
  set(ax(ii), 'YLimMode','manual');
  ud= get(ax(ii), 'userData');
  if iscell(ud.chan),
    ud.chan= ud.chan{1};  % for multiple channels per axis take only the first
  end
  ci= util_chanind(fv, ud.chan);
  if isempty(ci) && isempty(strmatch(ud.chan,opt.ShiftAlso,'exact')),
    continue;
  end
  pos= get(ax(ii), 'position');
  bar_pos= [pos(1:3) opt.Height*pos(4)];
  bar_pos(2)= pos(2) + opt.VPos*(1-opt.Height)*pos(4);
  switch(opt.ShiftAxes)
    case {1,'position'},
     if opt.VPos<0.5,
       new_pos= [pos(1) pos(2)+bar_pos(4) pos(3) pos(4)*(1-opt.Height)];
     else
       new_pos= [pos(1:3) pos(4)*(1-opt.Height)];
       if opt.VPos==1,
         axis_raiseTitle(ax(ii), opt.Height);
       end
     end
     set(ax(ii), 'position',new_pos);
    case {2,'ylim'},
     yLim= get(ax(ii), 'yLim');
     if opt.VPos<0.5,
       yLim(1)= yLim(1) - opt.Height*diff(yLim);
     else
       yLim(2)= yLim(2) + opt.Height*diff(yLim);
     end
     set(ax(ii), 'yLim',yLim);
   otherwise,
    error('ShiftAxes policy not known');
  end
  if isempty(ci),
    continue;
  end
  jj= jj+1;
  H.ax(jj)= axes('position', bar_pos);
  set(H.ax(jj), AxesStyle{:});
  hold on;      % otherwise axis properties like ColorOrder are lost
  if opt.UseLocalColormap,
    fv.x(fv.x > opt.CLim(2)) = opt.CLim(2);
    fv.x(fv.x < opt.CLim(1)) = opt.CLim(1);
    image_localColormap(fv.x(:,ci)', opt.Colormap, 'CLim',opt.CLim);
  else
    H.image(jj)= image(fv.x(:,ci)', 'cDataMapping','scaled');
  end
  set(H.ax(jj), AxesStyle{:});
  ud= struct('type','ERP plus: bar', 'chan', str_vec2str(fv.clab(ci)));
  set(H.ax(jj), 'userData',ud);
  hold off;
  if strcmp(get(H.ax(jj), 'box'), 'on'),
    set(H.ax(jj), 'LineWidth',0.3);
    axis_redrawFrame(H.ax(jj));
  end
end
if diff(opt.CLim)==0, opt.CLim(2)= opt.CLim(2)+eps; end
set(H.ax, 'xLim',[0.5 size(fv.x,1)+0.5], 'xTick',[], 'yTick',[], ...
          'CLim',opt.CLim);

if opt.MoveBack,
  obj_moveBack(H.ax);
end

if ~isempty(opt.HScale),
  pos= get(opt.HScale.ax, 'position');
  dh= pos(4)*(1-opt.ScaleHeight)*opt.ScaleVPos;
  pos_cb= [pos(1)+(1-opt.ScaleWidth-opt.ScaleLeftShift)*pos(3) pos(2)+dh ...
           opt.ScaleWidth*pos(3) opt.ScaleHeight*pos(4)];
  H.scale.ax= axes('position', pos_cb);
  CLim= opt.CLim;
  colbar= linspace(CLim(1), CLim(2), size(opt.Colormap,1));
  if opt.UseLocalColormap,
    image_localColormap(colbar', opt.Colormap, 'CLim',opt.CLim);
  else
    H.scale.im= imagesc(1, colbar, [1:size(opt.Colormap,1)]');
  end
  axis_redrawFrame(H.scale.ax);
  if opt.AlphaStepsMode,
    yTick= 0:nValues;
    if opt.UseLocalColormap, yTick= yTick+1; end %% Just a hack...
    set(H.scale.ax, 'yTick',yTick, 'yTickLabel',100*[1 fv.alpha]);
  else
    ticks= visutil_goodContourValues(CLim(1), CLim(2), -3);
    tickLabels= util_trunc(ticks, opt.ScaleDigits);
    if opt.UseLocalColormap,
      YLim= get(gca,'YLim');
      ticks= (ticks-CLim(1))*diff(YLim)/diff(CLim)+YLim(1);
    end
    set(H.scale.ax, 'yTick',ticks, 'yTickLabel',tickLabels);
  end
  set(H.scale.ax, 'xTick',[], 'tickLength',[0 0], 'YDir','normal', ...
                  'FontSize',opt.ScaleFontSize);
  if ~isempty(opt.ScaleUnit),
    yLim= get(H.scale.ax, 'yLim');
    H.scale.unit= text(1, yLim(2), opt.ScaleUnit);
    set(H.scale.unit, 'horizontalAli','center', 'verticalAli','bottom');
  end
end

if nargout==0,
  clear H;
end

author: bbci


text:
function [x,f,info]=convexQuadProgSOCP(G,a,C,b,numEqConst,params)
%%CONVEXQUADPROGSOCP Perform quadratic programming on a convex problem.
%                Specifically, this algorithm solves the optimization
%                problem
%                minimize_x a'*x+(1/2)*x'*G*x
%                such that C(:,1:numEqConst)'*x=b
%                      and C(:,(numEqConst+1):end)'*x>=b
%                Convexity means that G is positive (semi)definite. This
%                reformulates the problem as a second order cone problem
%                and then uses splittingConicSolver to solve the problem.
%                Unlike the active set method of [2], this formulation
%                tends to converge slowly on large problems (large being,
%                for example, 81 dimensions).
%
%INPUTS: G An nXn real, positive definite symmetric matrix (n>0).
%        a An nX1 real vector.
%        C An nXm real matrix (m>=0). The first numEqConst constraints (the
%          equality constraints) must be linearly independent or the
%          algorithm will identify the problem as infeasible.
%        b An mX1 vector.
% numEqConst The number of constraints in C that are equality constraints
%          (numEqConst<=m). If this parameter is omitted or an empty matrix
%          is passed, the default numEqConst=0 is used.
%   params Parameters that affect how the function splittingConicSolver,
%          which is used by this function, works. See the comments to
%          splittingConicSolver for more information.
%
%OUTPUTS: x The nX1 solution to the optimization problem or an empty matrix
%           if the problem is infeasible or unbounded.
%         f The value of the objective function at x, or an empty matrix
%           if the problem is infeasible or unbounded.
%      info The problem is formulated as a second-order cone optimization
%           problem. This is the information regarding the termination
%           state returned by the splittingConicSolver function.
%
%As shown in [1], the equivalent second order cone optimization problem is
%minimize t
%subject to norm(SG*x+SG\a)<=t
%       and C(:,1:numEqConst)'*x=b(1:numEqConst)
%       and C(:,(numEqConst+1):end)'*x>=b((numEqConst+1):end)
%       where SG=sqrtm(G)
%The function splittingConicSolver is used to solve the problem.
%
%EXAMPLE 1:
%This is the example worked out step-by-step in [2]:
% G=[4, -2;-2, 4];
% b=[0;0;2];
% a=[6;0];
% C=[1,0,1;
%    0,1,1];
% [x,f,info]=convexQuadProgSOCP(G,a,C,b)
%The optimal solution can be seen to be x=[0.5;1.5] and f=6.5;
%
%EXAMPLE 3:
%This is an example of a problem having both inequality and equality
%constraints:
% G=[4, -2;-2, 4];
% b=[2;1.6];
% a=[6;0];
% C=[1,0;
%    1,1];
% numEqConst=1;%One equality constraint.
% [x,f,exitCode]=convexQuadProgSOCP(G,a,C,b,numEqConst)
%The optimal solution is x=[0.4;1.6] and f=6.56.
%
%REFERENCES:
%[1] S. Lobo, Miguel, L. Vandenberg, S. Boyd, and H. Lebret, "Applications
%    of second-order cone programming," Linear Algebra and its
%    Applications, vol. 284, no. 1, pp. 193-228, 1998.
%[2] D.Goldfarb and A. Idnani, "A numerically stable dual method for
%    solving strictly convex quadratic programs," Mathematical Programming,
%    vol. 27, no. 1, pp. 1-33, Sep. 1983.
%
%February 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(numEqConst))
   numEqConst=0;
end

if(nargin<6)
   params=[];
end

if(numEqConst>size(C,1))
   error('numEqConst is larger than the number of constraints in C')
end

P0=G;
%The real command can help deal with finite precision issues.
P0Root=real(sqrtm(P0));
q0=a;
Aeq=C(:,1:numEqConst)';
beq=b(1:numEqConst);
ALeq=-C(:,(numEqConst+1):end)';
bLeq=-b((numEqConst+1):end);

numIneqConst=length(bLeq);

%The extra variable t is appended to the end.
numDim=length(a);
c=zeros(numDim+1,1);
c(end)=1;

%Allocate space.
numVar=numEqConst+numIneqConst+numDim+1;
F=zeros(numVar,numDim+1);
g=zeros(numVar,1);

%First, the equality constraints.
F(1:numEqConst,1:numDim)=Aeq;
g(1:numEqConst)=beq;
curRow=numEqConst+1;
cone=[];
cone.f=numEqConst;

%Next, the inequality constraints
sel=curRow:(curRow+numIneqConst-1);
F(sel,1:numDim)=ALeq;
g(sel)=bLeq;
cone.l=numIneqConst;
curRow=curRow+numIneqConst;

%Finally, the quadratic constraint.
F(curRow,:)=-c';
g(curRow)=0;
curRow=curRow+1;

sel=curRow:(curRow+numDim-1);
F(sel,1:numDim)=-P0Root;
%g(sel)=P0Root\q0;
g(sel)=pinv(P0Root)*q0;
cone.q=numDim+1;

[x,~,~,info]=splittingConicSolver(F,g,c,cone,params);

if(~isempty(x))
    x=x(1:numDim);%Get rid of the auxiliary variable.
    f=a'*x+(1/2)*x'*G*x;
else
    f=[];
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [xyz, XYZmm, Z, class] = cluster_local_maxima(cl, dthresh, verbose)
% Clusters are chosen so that they must be at least dthresh mm apart
% default is 10 mm
%
% :Usage:
% ::
%
%    [xyz, XYZmm, Z, class] = cluster_local_maxima(cl, [dthresh], [verbose])
%
% verbose output: 1/0, default is 0
%
% additional optional outputs (slower):
%
% class: vector of integers for which subcluster this cluster belongs to
%
% ..
%    tor wager
% ..

    if nargin < 2 || isempty(dthresh), dthresh = 10; end
    if nargin < 3, verbose = 0; end

    [N Z xyz] = spm_max(abs(cl.Z), cl.XYZ);

    XYZmm = voxel2mm(xyz, cl.M);

    ncoord = size(cl.XYZmm, 2);
    class = zeros(1, ncoord);

    if isempty(xyz), return, end   % no peak maximum

    if verbose, fprintf(1, 'Local maxima (initial): %3.0f\n', length(Z)); end

    % find maxima within 10 mm and collapse
    d = pdist(XYZmm');
    nd = length(d);
    tooclose = d < dthresh;

    while any(tooclose)

        [d, nd, tooclose, xyz, XYZmm, Z] = omit_lowz_of_closest_pair(d, nd, tooclose, xyz, XYZmm, Z, dthresh, verbose);

    end

    if verbose, fprintf(1, 'Local maxima (final): %3.0f\n', length(Z)); end

    if nargout > 3
        % assign each voxel a subcluster based on closest local max
        speaks = XYZmm';
        npeaks = length(Z);

        % this is slower even for only 1000 vox
        %d = squareform(pdist([speaks; cl.XYZmm']));
        %d = d(1:npeaks, npeaks+1:end);   % peaks x voxels

        for i = 1:ncoord
            d = distance(cl.XYZmm(:, i)', speaks);
            [mind, whclose] = min(d);
            class(i) = whclose(1);
        end
    end
end

function [d, nd, tooclose, xyz, XYZmm, Z] = omit_lowz_of_closest_pair(d, nd, tooclose, xyz, XYZmm, Z, dthresh, verbose)

    tooclose = tooclose .* d;

    % find minimum distance
    tooclose(tooclose == 0) = Inf;
    [mindist, wh] = min(tooclose);

    % get indices of two local maxima that are too close (in rows)
    closest = zeros(1, nd);
    closest(wh(1)) = mindist;
    closest = squareform(closest);
    [rows, cols] = find(closest);    % need to return cols to get correct behavior

    % find which has the lowest Z-score
    [lowZ, whz] = min(Z(rows));      % rows is sufficient b/c there is only one pair

    wh_to_omit = rows(whz); % which in list to omit

    if verbose
        fprintf(1, 'Omitting index: %3.0f with z-value of %3.2f, which is %3.2f mm from another max whose z is %3.2f\n', ...
            wh_to_omit, lowZ, mindist, max(Z(rows)));
    end

    % eliminate important variables
    xyz(:, wh_to_omit) = [];
    XYZmm(:, wh_to_omit) = [];
    Z(wh_to_omit) = [];

    % find maxima within 10 mm and collapse
    d = pdist(XYZmm');
    nd = length(d);
    tooclose = d < dthresh;

end

function d = distance(u, v)
% d = distance(u, v)
% Euclidean distance between u and v
% u and v should be column vectors or matrices in k-dim space, where k is
% the number of columns of u and v.
%
% if u is a single row, replicates u to size(v,1)

if size(u,1) < size(v,1)
    u = repmat(u,size(v,1),1);
    if size(u,1) < size(v,1), error('Inputs must have same number of rows, or one row for 1st input!'),end
end

delt = u - v;
d = (sum(delt .^ 2,2)) .^.5;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%
% This code calculates the interevent distances and the correlation integral
% of a subset selected with the sampling sphere. The code is called from circlefd.m.
% Francesco Pacchiani 1/2000
%
% Calculation of the 3D distances between all possible pairs
% (combination of n epicenters taken 2 at a time) of earthquakes of
% the given dataset.
%
%
% Variables
%
N = size(E,1);				% N= # of events in the catalogue; E= Earthquake catalogue
pairdist = []; 			% pairdist= Vector of interevent distances
j = nchoosek(N,2);			% j= # of interevent distances calculated
pairdist = zeros(j,1);
k = 0;
%E.Latitude= (max(Da(:,2))+min(Da(:,2)))/2;

Ho_Wb = waitbar(0,'Calculating the fractal dimension D');
Hf_Cfig = gcf;
Hf_child = get(groot,'children');
set(Hf_child,'pointer','watch','papertype','A4');
%
%
% Calculation of the interevent distances in 2D plus the depths differences.
%
%
for i = 1:(N-1)

    lon1 = repmat(E(i,1), [(N-i),1]);
    lat1 = repmat(E(i,2), [(N-i),1]);

    lon2 = E((i+1):end, 1);
    lat2 = E((i+1):end, 2);

    pairdist(k+1:k + size(lon1, 1)) = distance(lat1,lon1,lat2,lon2);

    k = k + size(lon1,1);

    waitbar((0.75/(N-1))*i, Ho_Wb);

end

clear i j k;
%
%
% Conversion of the interevent distances from degrees to kilometers and calculates
% the interevent distances in three dimensions.
%
%
if dtokm == 1
    pairdist = pairdist.*111;
end
%
%
% Calculation of the correlation integral using as input the
% pair distances computed above.
%
%
% Variables
%
d = 2;						%d = the dimension of the embedding volume.
rmax = max(pairdist);
rmin = min(pairdist);

if rmin == 0
    rmin = 0.01;
end

lrmin = log10(rmin);
lrmax = log10(max(pairdist));

%u = (log10(rmin):0.15:log10(rmax))';
%
% Defining the distance vector r in order that on the
% log-log graph all the points plot at equal distances from one another.
%
r = (logspace(lrmin, lrmax, 50))';
%r = zeros(size(u,1),1);
%r = 10.^u;
%
%
corint = [];						% corint= Vector of ?cumulative? correlation integral values for increasing interevent radius
corint = zeros(size(r,1),1);
k = 1;

for i = 1:size(r,1)

    j = [];
    j = pairdist < r(i);
    corint (k,1) = (2/(N*(N-1)))*sum(j);
    k = k + 1;
    waitbar(0.75 + (0.25/size(r,1))*i,Ho_Wb);

end

clear i j k;
close(Ho_Wb);
Hf_child = get(groot,'children');
set(Hf_child,'pointer','arrow');
%
%
dofd = 'fd';

author: CelsoReyes


text:
function polygon_properties_test11 ( )

%*****************************************************************************80
%
%% POLYGON_PROPERTIES_TEST11 tests POLYGON_LATTICE_AREA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_PROPERTIES_TEST11\n' );
  fprintf ( 1, '  POLYGON_LATTICE_AREA returns the "area"\n' );
  fprintf ( 1, '  of a polygon, measured in lattice points.\n' );

  i = 5;
  b = 6;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of interior lattice points = %d\n', i );
  fprintf ( 1, '  Number of boundary lattice points = %d\n', b );

  area = polygon_lattice_area ( i, b );

  fprintf ( 1, '  Area of polygon is %g\n', area );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [temp_planC] = loadPlanC_temp(fname,file)
%
%   Load and return a planC given a fullpath filename.  planC is returned as
%   temp_planC in the case where another planC is already open. The filename can
%   be a .mat or .mat.bz2 file. Bz2 files are extracted to a temporary .mat
%   file, loaded, and then the temporary .mat file is deleted.
%
%   18 Mar 06, KU
%
% 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/>.

planC = [];
temp_planC  = [];

h = waitbar(100,'Loading selected study...');

bzFile = 0;

% Test for compressed file (*.gz), GNUZIP (freeware).  Code by A. Blanco.
if ~isempty(strfind(file,'.bz2'))
    bzFile = 1;
    disp(['Uncompressing file ', fname]);
    outstr = gnuCERRCompression(file, 'uncompress',tempdir);
    file = file(1:end-3);
end

disp(['Loading file ', fname]);
planC = load(file,'planC');

temp_planC = planC.planC;   %Conversion from struct created by load

%Remove unzipped file after loading.
if bzFile
    delete(file);
end

close(h)

if ~exist('planC');
    error('.mat or .mat.bz2 file does not contain a planC variable.');
    return;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%  Figure 10.41      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_41.m is a script to create Figure 10.41, the step response of the
% attitude inner loop to a 2 degree step in pitch angle.
clf;
% the total inner loop system matrix is

ftq =[-0.0064    0.0263         0  -32.2000         0;
   -0.0941   -0.6240  761.1400 -196.2000         0;
   -0.0002   -0.0015   -4.4120  -12.4800         0;
         0         0    1.0000         0         0;
         0   -1.0000         0  830.0000         0];

g =[0;
  -32.7000;
   -2.0800;
         0;
         0];
h = [0     0    -1     0     0];
j =0;
% now remove the altitude to get the inner loop pitch angle
ft=ftq(1:4,1:4)
gt=g(1:4);
ht=[0 0 0 1];
% compute the dc gain
dcgain= -ht*inv(ft)*gt
% compute the reference input in radians from 2 degrees

r=2*pi/180
t=0:.1:8;
y=step(ft,r*gt/dcgain,ht,0,1,t);
plot(t,y);
xlabel('Time (sec)');
ylabel('\theta (t) (rad)');
title(' Fig. 10.41 step response of the pitch angle inner loop ');
%grid
nicegrid


text:
function kern = lfmaKernParamInit(kern)

% LFMAKERNPARAMINIT LFMA kernel parameter initialisation.
% The kernel is designed to interoperate with the multiple output block
% kernel so that f(t) can be inferred given several different
% instantiations of x(t).
%
% The parameters (m, c, delta and k) are constrained positive.
%
% FORMAT
% DESC initialises the latent force model kernel structure with some
% default parameters.
% RETURN kern : the kernel structure with the default parameters placed in.
% ARG kern : the kernel structure which requires initialisation.
%
% SEEALSO : kernCreate, kernParamInit, lfmKernCompute
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

% A wrapper that calls lfmKernParamInit to initialize the parameters of the
% model.

kern = lfmKernParamInit(kern);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function importfile(fileToRead1)
%IMPORTFILE(FILETOREAD1)
%  Imports data from the specified file
%  FILETOREAD1:  file to read

%  Auto-generated by MATLAB on 08-Mar-2011 12:24:52
%  Copyright 2011 The MathWorks, Inc.

% Import the file
sheetName='Sheet1';
[numbers, strings, raw] = xlsread(fileToRead1, sheetName);
if ~isempty(numbers)
    newData1.data =  numbers;
end

if ~isempty(strings) && ~isempty(numbers)
    [strRows, strCols] = size(strings);
    [numRows, numCols] = size(numbers);
    likelyRow = size(raw,1) - numRows;
    % Break the data up into a new structure with one field per column.
    if strCols == numCols && likelyRow > 0 && strRows >= likelyRow
        newData1.colheaders = strings(likelyRow, :);
    end
end

% Create new variables in the base workspace from those fields.
for i = 1:size(newData1.colheaders, 2)
    assignin('base', genvarname(newData1.colheaders{i}), newData1.data(:,i));
end


text:
function tone=note(keynum, dur)

fs=11025;
tt = 0:(1/fs):dur;

%  This generates white noise for whatever duration specified..  good snare
%  sound with the right envelope, but matlab developed timing problems with this.  I'll probably
%  try to get these running sometime after lab
 if keynum == 2
     tone=rand(1,length(tt));
     return;
 end

%generates rests
if keynum == 0
%for kk = 1:length(tt)
    tone([1:length(tt)]) = 0;
    %end
return;
end

%adding these octaves rounded out the sound.
freq=440*2^((keynum-49)/12);
freq3=freq*3;
freq5=freq*5;
freq9=freq*9;
freq7=freq*7;
tone1 = .75*sin(2*pi*freq*tt);
 tone3 = .65*sin(2*pi*freq3*tt);
 tone5=.5*sin(2*pi*freq5*tt);
 tone9 = .222*sin(2*pi*freq9*tt);
 tone7 = .12*sin(2*pi*freq7*tt);
 tone12= 1*sin(2*pi*freq*12*tt);

tone=tone1+tone3+tone5+tone7+tone9;%+tone12;


text:
function [cm3] = mi32cm3(mi3)
% Convert volume from cubic miles to cubic centimeters.
% Chad Greene 2012
cm3 = mi3*4168181825400000 ;


text:
function G = diffy( F, n )
%DIFFY   Differentiate a DISKFUNV with respect to its second argument
%   DIFFY(F) returns a DISKFUNV representing the derivative of F in its
%   second argument. This is the same as DIFF(F,2,1).
%
%   DIFFY(F,N) returns a DISKFUNV representing the Nth derivative of F in
%   its second argument. This is the same as DIFF(F, 2, N).
%
%   This command is for convenience as the syntax for DIFF, inherited from
%   the DIFF command for matrices, can be confusing.
%
% See also DIFFX, DIFF.

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

if ( nargin == 1 )
    % Default to first derivative.
    n = 1;
end

% Call DISKFUNV/DIFF:
G = diff( F, 2, n);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [x,y] = circlePoints(radSpacing)
%Make samples on a circle
%
%  [x,y] = circlePoints(radSpacing)
%
% Example
%    radSpacing = 0.11;
%    [x,y] = circlePoints(radSpacing);
%    plot(x,y,'o'); axis equal
%

theta = (0:radSpacing:2*pi);
x = cos(theta); y = sin(theta);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function iimg_write_images(dat, volInfo, outnames)
% Write a series of Analyze images given inputs.
%
% :Usage:
% ::
%
%     iimg_write_images(dat,volInfo,outnames)
%
% :Inputs:
%
%   **dat:**
%        voxels x images matrix of image data in index format
%
%   **volInfo:**
%        spm-style info structure; see iimg_read_img
%
%   **outnames:**
%        a string matrix (or cell array) of output filenames, or the name of a 4-D file to create

    if isempty(outnames)
        error('You must specify a string array of output names, or a single 4-D image file name');
    end

    if iscell(outnames)
        outnames = char(outnames{:});
    end

    n = size(dat,2);
    if size(outnames,1) ~= n
        error('Number of images (columns of dat) must == number of output img names.');
    end

    if volInfo.nvox ~= size(dat,1)
        error('Dims of dat (vox x images) and volInfo do not match.');
    end

    is4d = 0;
    if size(dat, 2) > 1 && size(outnames, 1) == 1, is4d = 1; end

    if is4d
        iimg_reconstruct_vols(dat,volInfo,'outname',deblank(outnames));

    else
        for i = 1:n
            %iimg_reconstruct_3dvol(dat(:,i),volInfo,'outname',deblank(outnames(i,:)));
            % This handles volume creation, etc.  It also works for a 4-D
            % image (but all volumes must be created at once)
            iimg_reconstruct_vols(dat(:,i),volInfo,'outname',deblank(outnames(i,:)));
        end
    end

end


text:
function ni_def = afni_ni_defs()
    persistent cached_ni_def

    % for efficency reasons generate the definitions only once;
    % upon each subsequent call the cached data is used
    if ~isstruct(cached_ni_def)
        cached_ni_def=generate_ni_def();
    end

    ni_def=cached_ni_def;

function ni_def=generate_ni_def()
    ni_def=struct();
    %work out types /* definitions from niml.h
    ni_def.NI_BYTE     =   0        ;       %/* == MRI_byte    */
    ni_def.NI_SHORT    =   1       ;        %/* == MRI_short   */
    ni_def.NI_INT      =   2       ;        %/* == MRI_int     */
    ni_def.NI_FLOAT32  =   3       ;        %/* == MRI_float   */
    ni_def.NI_FLOAT    =   ni_def.NI_FLOAT32;
    ni_def.NI_FLOAT64  =   4      ;         %/* == MRI_double  */
    ni_def.NI_DOUBLE   =   ni_def.NI_FLOAT64;
    ni_def.NI_COMPLEX64=   5   ;            %/* == MRI_complex */
    ni_def.NI_COMPLEX  =   ni_def.NI_COMPLEX64;
    ni_def.NI_RGB      =   6   ;            %/* == MRI_rgb     */
    ni_def.NI_RGBA     =   7  ;             %/* == MRI_rgba    */
    ni_def.NI_STRING   =   8  ;             %/* after "basic" types */

    %/*! Number of types of fixed size ("basic" types).
    %    Note that if this changes,
    %    the NI_rowtype stuff must be altered accordingly. */
    ni_def.NI_NUM_BASIC_TYPES = 8;

    %/*! One more than the last NI_ data type code defined above. */
    ni_def.NI_NUM_TYPES       = 9;

    ni_def.type_name = strvcat({
      'byte'  ; 'short'  ; 'int'     ;
      'float' ; 'double' ; 'complex' ;
      'rgb'   ; 'rgba'   ; 'String' });

    ni_def.type_alias = strvcat({   %/* aliases */
     'uint8'   ; 'int16'   ; 'int32'     ;
     'float32' ; 'float64' ; 'complex64' ;
     'rgb8'    ; 'rgba8'   ; 'CString' });

    % TODO: see how to deal with other types
    undef=@()error('Unsupported type');

    ni_def.type = {...
       @uint8  ; @int16    ; @int32      ; ...
       @single ; @double   ; undef       ; ...
       undef   ; undef     ; @char       };

    ni_def.size_bytes=[...
       1       ; 2          ; 4          ; ...
       4       ; 8          ; NaN        ; ...
       NaN     ; NaN        ; 1          ];

    ni_def.init = 1;

    ni_def.escape={'<'  ,'<' ;...
                   '>'  ,'>' ;...
                   '"','"' ;...
                   ''','''';...
                   '&' ,'&' };

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%define your parameters here
%{
%%mnist mlp
opts.dataset_name='mnist';
network_name='cnn'; %network name
use_gpu=1; %use gpu or not
PrepareDataFunc=@PrepareData_MNIST_MLP;
NetInit=@net_init_mlp_mnist;
use_selective_sgd=1;
ssgd_search_freq=10;
learning_method=@sgd; %training method: @sgd,@adagrad,@rmsprop,@adam
sgd_lr=1e-2;
opts.n_epoch=20;
opts.LoadResults=0;
%}

%{
n_epoch=20; %training epochs
dataset_name='mnist'; %dataset name
network_name='cnn'; %network name
use_gpu=1; %use gpu or not

%function handle to prepare your data
PrepareDataFunc=@PrepareData_MNIST_CNN;
%function handle to initialize the network
NetInit=@net_init_cnn_mnist;

%automatically select learning rates
use_selective_sgd=1;
%select a new learning rate every n epochs
ssgd_search_freq=10;
learning_method=@sgd; %training method: @sgd,@adagrad,@rmsprop,@adam

%sgd parameter
%(unnecessary if selective-sgd is used)
%sgd_lr=5e-2;
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end

%opts=[];
%addpath('./CoreModules');
opts.n_epoch=n_epoch; %training epochs
opts.dataset_name=dataset_name; %dataset name
opts.network_name=network_name; %network name
opts.use_gpu=use_gpu; %use gpu or not

if ~isfield(opts,'datatype')
    opts.datatype ='single';
end

if ~isfield(opts,'LoadNet')
    opts.LoadNet=0;
end

opts.dataDir=['./',opts.dataset_name,'/'];
opts=PrepareDataFunc(opts);

%opts.parameters=[];

opts.parameters.current_ep=1;

%%%parameters in the training

opts.parameters.learning_method=learning_method;
opts.parameters.selective_sgd=use_selective_sgd;%call selective-sgd

%selective-sgd parameters
if opts.parameters.selective_sgd==1
    if ~isfield(opts.parameters,'search_iterations')
        opts.parameters.search_iterations=30;%iterations used to determine the learning rate
    end
    opts.parameters.ssgd_search_freq=ssgd_search_freq;%%search every n epoch
    if ~isfield(opts.parameters,'lrs')

        opts.parameters.lrs =[1,.5];%initialize selection range
        if ~strcmp(func2str(opts.parameters.learning_method),'sgd')&&~strcmp(func2str(opts.parameters.learning_method),'sgd2')
            opts.parameters.lrs =opts.parameters.lrs.*1e-1;
        end
        opts.parameters.lrs=[opts.parameters.lrs,opts.parameters.lrs*1e-1,opts.parameters.lrs*1e-2,opts.parameters.lrs*1e-3];%initialize selection range
    end
    opts.parameters.selection_count=0;%initialize
    opts.parameters.selected_lr=[];%initialize
end

if opts.parameters.selective_sgd==0
    opts.parameters.lr =sgd_lr;
end

%%sgd parameters
if ~isfield(opts.parameters,'mom')
    opts.parameters.mom =0.9;
end

%adam parameters
if strcmp(func2str(opts.parameters.learning_method),'adam')
    if ~isfield(opts.parameters,'mom2')
        opts.parameters.mom2 =0.999;
    end
end

if ~isfield(opts.parameters,'batch_size')
    opts.parameters.batch_size=500;
end
if ~isfield(opts.parameters,'weightDecay')
    opts.parameters.weightDecay=1e-4;
end

opts=generate_output_filename(opts);

if ~isfield(opts,'plot')
    opts.plot =1;
end

if ~isfield(opts,'LoadResults')
    opts.LoadResults=0;
end

if ~opts.LoadResults
    TrainingScript();
end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
classdef MMAE < Metric
    %MAE static class to calculate the minimum mean absolute error (MAE) per
    %   class. Values range from 0 to J-1, where J is the number of classes.
    %
    %   MAE methods:
    %      CALCULATEMETRIC            - Computes the evaluation metric
    %      CALCULATECROSSVALMETRIC    - Computes the evaluation metric as an error
    %
    %   References:
    %     [1] M. Cruz-Ramírez, C. Hervás-Martínez, J. Sánchez-Monedero and
    %         P. A. Gutiérrez Metrics to guide a multi-objective evolutionary
    %         algorithm for ordinal classification, Neurocomputing, Vol. 135, July, 2014, pp. 21-31.
    %         https://doi.org/10.1016/j.neucom.2013.05.058
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.htmlml
    methods
        function obj = MMAE()
            obj.name = 'Max Mean Absolute Error';
        end
    end

    methods(Static = true)

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        % Function: calculateMetric (static)
        % Description: Computes the evaluation metric
        % Outputs: metric results
        % Arguments:
        %           argum1--> First argument (confusion matrix or predictions)
        %	    argum2--> Second argument (true labels)
        % 	    If there is only one argument, the results are computed
        %	    using the confusion matrix. In other case, with the
        %	    predictions and true labels.
        %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        function maxmae = calculateMetric(argum1,argum2)
            %CALCULATEMETRIC Computes the evaluation metric
            %   METRIC = CALCULATEMETRIC(CM) returns calculated metric from confussion
            %   matrix CM
            %   METRIC = CALCULATEMETRIC(actual, pred) returns calculated metric from
            %   real labels (ACTUAL) labels and predicted labels (PRE
            if nargin == 2
                argum1 = confusionmat(argum1,argum2);
            end
            n=size(argum1,1);
            cm = double(argum1);
            cost = abs(repmat(1:n,n,1) - repmat((1:n)',1,n));
            mae = zeros(n:1);
            cmt = cm';
            for i=0:n-1
                mae(i+1) = sum(cost(1+(i*n):(i*n)+n).*cmt(1+(i*n):(i*n)+n)) / sum(cmt(1+(i*n):(i*n)+n));
            end
            maxmae = max(mae);
        end

        function value = calculateCrossvalMetric(argum1,argum2)
            %CALCULATECROSSVALMETRIC Computes the evaluation metric and returns
            %it as an error.
            %   METRIC = CALCULATECROSSVALMETRIC(CM) returns calculated metric from confussion
            %   matrix CM
            %   METRIC = CALCULATECROSSVALMETRIC(actual, pred) returns calculated metric from
            %   real labels (ACTUAL) labels and predicted labels (PRED)
            if nargin == 2
                value = MMAE.calculateMetric(argum1,argum2);
            else
                value = MMAE.calculateMetric(argum1);
            end
        end

    end

end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
classdef RankAwareOMP < handle
    % Implements Rank Aware OMP algorithm for sparse approximation in MMV case

    properties
        % Maximum residual norm
        MaxResNorm = 1e-4
        % Indicates if we should stop on exceeding residual norm
        StopOnResidualNorm = true
        % Indicates if we should stop when residual norm stops improving
        StopOnResNormStable = true
        % Maximum number of iterations for approximation
        MaxIters
        % Indicates if log messages should be printed at each iteration
        Verbose = false
        % Minimum Sparsity
        MinK = 4
        % Ignored atom (which won't be considered in identification step)
        IgnoredAtom = -1
        % The norm to be chosen for rows
        P
    end

    properties(SetAccess=private)
        % The dictionary
        Dict
        % Ambient signal dimensions
        N
        % Number of atoms in dictionary
        D
        % Sparsity level of representations (may be negative)
        K
        % Result of a solver
        result
    end

    methods
        function self  = RankAwareOMP(Dict, K)
            % We assume that all the columns in dictionary are normalized.
            if isa(Dict, 'spx.dict.Operator')
                self.Dict = Dict;
            elseif ismatrix(Dict)
                self.Dict = spx.dict.MatrixOperator(Dict);
            else
                error('Unsupported operator.');
            end
            [self.N, self.D] = size(Dict);
            if ~exist('K', 'var')
                % No sparsity level has been pre-specified.
                K = -1;
            end
            self.K = K;
            % Maximum number of iterations
            maxIter = self.N;
            if K > 0
                % We have to consider pre-specified sparsity level
                maxIter = K;
            end
            self.MaxIters = maxIter;
            if K > 0 && self.MinK >= K
                self.MinK = 0;
            end
        end

        function result  = solve(self, Y)
            % Initialization
            % Solves approximation problem using OMP
            d = self.D;
            n = self.N;
            % The number of signals being approximated.
            s = size(Y, 2);
            % Initial residual
            R = Y;
            dict = self.Dict;
            % Active indices
            omega = [];
            % Estimate
            Z = zeros(d, s);
            if self.StopOnResNormStable
                oldResNorm = norm(R, 'fro');
            end
            maxIter = self.MaxIters;
            P = self.P;
            for iter=1:maxIter
                % Orthonormalize Residual
                U = orth(R); % R is n x s. U is n x r.
                % Compute inner products (DxN  * N x r = D x r )
                innerProducts = apply_ctranspose(dict, U);
                % Compute absolute values of inner products
                innerProducts = abs(innerProducts);
                % Compute l2 norm inner products over each row
                innerProducts = spx.norm.norms_l2_rw(innerProducts);
                % Mark the inner products of already selected columns as 0.
                innerProducts(omega) = 0;
                % Find the highest inner product
                [~, index] = max(innerProducts);
                % Add this index to support
                omega = [omega, index];
                % Solve least squares problem
                subdict = columns(dict, omega);
                tmp = linsolve(subdict, Y);
                % Updated solution
                Z(omega, :) = tmp;
                % Let us update the residual.
                R = Y - dict.apply(Z);
                resNorm = norm(R, 'fro');
                if self.StopOnResidualNorm || self.StopOnResNormStable
                    if resNorm < self.MaxResNorm
                        break;
                    end
                    if self.StopOnResNormStable
                        change = abs(oldResNorm  - resNorm);
                        if change/oldResNorm < .01
                            % No improvement
                            break;
                        end
                    end
                end
            end

            % Solution vector
            result.Z = Z;
            % Residual obtained
            result.R = R;
            % Number of iterations
            result.iterations = iter;
            % Solution support
            result.support = omega;
            % residual Frobenius norm
            result.residual_frobenius_norm = resNorm;
            self.result = result;
        end

    end
end

author: indigits

repo: sparse-plex


text:
function tr = multitrace(A)
% Computes the traces of the 2D slices in a 3D matrix.
%
% function tr = multitrace(A)
%
% For a 3-dimensional matrix A of size n-by-n-by-N, returns a column vector
% tr of length N such that tr(k) = trace(A(:, :, k));
%
% See also: multiprod multitransp multiscale

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    assert(ndims(A) <= 3, ...
           ['multitrace is only well defined for matrix arrays of 3 ' ...
            'or less dimensions.']);

	tr = diagsum(A, 1, 2);

end

author: andrewssobral

repo: lrslibrary


text:
function [scanFileNameC,doseFileNameC,maskFileNameC] = ...
    export_scan_dose_struct_from_planC(exportFormat,...
    scanNumV,doseNumV,structNameC,outdDir,planC)
% function [scanFileNameC,doseFileNameC,maskFileNameC] = ...
%     export_scan_dose_struct_from_planC(exportFormat,...
%     scanNumV,doseNumV,structNameC,outNrrdDir,planC)
%
% exportFormat: nii or nrrd
% scanNumV: indices of scans in planC to export
% doseNumV: indices of doses in planC to export
% structNameC: cell array of structure names to export
% outDir: directory to export to
% planC:
%
% APA, 9/13/2021

indexS = planC{end};
reorientFlag = 1;
dataType = [];

scanFileNameC = {};
doseFileNameC = {};
maskFileNameC = {};

% scan
for iScan = 1:length(scanNumV)
    scanNum = scanNumV(iScan);
    scanFileNameC = scan2imageOut(planC,scanNum,outdDir,reorientFlag,exportFormat,dataType);
end

% dose
for iDose = 1:length(doseNumV)
    doseNum = doseNumV(iDose);
    scanNum = getDoseAssociatedScan(doseNum,planC);
    doseFileNameC = dose2imageOut(planC, doseNum, scanNum, outdDir,reorientFlag,exportFormat);
end

% structure
maskFileNameC = mask2imageOut(planC,scanNum,structNameC,outdDir,reorientFlag,exportFormat);


text:
function elementalWeightMatrix = getElementalWeightMatrix(isotopeAbundance)
% Get the symbols for all chemical elements and their atomic weights
%
% USAGE:
%    elementalWeightMatrix = getElementalWeightMatrix(isotopeAbundance)
%
% OPTIONAL INPUTS:
%    isotopeAbundance:         * false to return standard atomic weights (default)
%                              * true to return the weights of naturally predominant isotopes for biological
%                                elements and standard weights for other elements.
%                              * `m` x 3 cell arrray with user defined isotope abundance:
%                                isotopeAbundance{i, 1} = 'Atomic_Symbol';
%                                isotopeAbundance{i, 2} = Mass_Number;
%                                isotopeAbundance{i, 3} = abundance;
%                                (where sum of abundances of all isotopes of an element must be one)
%    printWarning:             true to print warning if there are unknown elements in isotopeAbundance
%
% OUTPUT:
%    elementalWeightMatrix:    #elements x 2 cell array, the 1st coloumn being the symbols
%                              of chemial elements and the 2nd their standard atomic weights

if nargin == 0 || isempty(isotopeAbundance)
    isotopeAbundance = false;
end

atomicWeights = parse_Atomic_Weights_and_Isotopic_Compositions_for_All_Elements;
[symbols, ia] = unique([atomicWeights.data.AtomicSymbol]);
% remove XYZ (non-existing element) from the symbol list
elementRemove = strcmp(symbols, 'D') | strcmp(symbols, 'T') | strncmp(symbols, 'XYZ', 3);
[symbols, ia] = deal(columnVector(symbols(~elementRemove)), ia(~elementRemove));
% identify elements without standard weight
elementWoMW = arrayfun(@(x) isempty(x.StandardAtomicWeight), atomicWeights.data(ia));
% elements with standard weight come first
symbols = [symbols(~elementWoMW); symbols(elementWoMW)];
elementWeights = [[atomicWeights.data(ia(~elementWoMW)).StandardAtomicWeight], NaN(1, sum(elementWoMW))]';

if ~isscalar(isotopeAbundance) || isotopeAbundance

    % biological elements as defined in getMolecularMass.m
    %allBiologicalElements = {'C', 'O', 'P', 'N', 'S', 'H', 'Mg', 'Na', 'K', 'Cl', 'Ca', 'Zn', 'Fe', 'Cu', 'Mo', 'I'};
    if isscalar(isotopeAbundance)
        % use the defaulted most naturally predominant isotope of each biological element
        isotopeAbundance = {...
            'C', 12, 1; ...  % Carbon, all as C12
            'H', 1, 1; ...  % Hydrogen, all as H1
            'O', 16, 1; ...  % Oxygen, all as O16
            'N', 14, 1; ...  % Nitrogen, all as N14
            'S', 32, 1; ...  % Sulphur, all as S32 (94% naturally S32)
            'Mg', 24, 1; ...  % Magnesium, all as Mg24, (79 % naturally Mg24)
            'K', 39, 1; ...  % Potassium, all as K39, (93% naturally K39)
            'Cl', 35, 1; ...  % Chlorine, all as Cl35 (75.78% naturally Cl35)
            'Ca', 40, 1; ...  % Calcium, all as Ca40 (96% naturally Ca40)
            'Zn', 64, 1; ...  % Zinc, all as Zn64 (48% naturally Zn64)
            'Cu', 63, 1; ...  % Copper, all as Cu63 (69% naturally Cu63)
            'Mo', 98, 1; ...  % Molybdenum, all as Mo98 (24% naturally Mo98)
            'P', 31, 1; ...  % Phosphorus, all as P31
            'Na', 23, 1; ...  % Sodium, all as Na23
            'Fe', 56, 1; ...  % Iron, all as Fe56
            'I', 127, 1; ...  % Iodine, all as I127
            };
    end

    % get the list of all isotopes
    [isotopes, isoMassNums, isoWeights] = deal([atomicWeights.data.AtomicSymbol], ...
        [atomicWeights.data.MassNumber], [atomicWeights.data.RelativeAtomicMass]);
    % treat 'D' and 'T' as the same element as 'H'
    isotopes(strcmp(isotopes, 'D') | strcmp(isotopes, 'T')) = {'H'};
    isotopeAbundance(strcmp(isotopeAbundance(:, 1), 'D') | ...
        strcmp(isotopeAbundance(:, 1), 'T')) = {'H'};

    [eleWtIsoAb, ~, ib] = unique(isotopeAbundance(:, 1));
    [ynE, idE] = ismember(eleWtIsoAb, symbols);

    % sanity check
    if size(isotopeAbundance, 2) ~= 3
        error('Incorrect size. IsotopeAbundance must be a n-by-3 cell array.')
    elseif ~all(ynE)
        % unknown elements in isotope abundance profile
        error('%s in the supplied isotope abundance profile is/are not standard chemical elements.', strjoin(eleWtIsoAb(~ynE), ', '))
    end

    massNumbers = cell2mat(isotopeAbundance(:, 2));
    abundances = cell2mat(isotopeAbundance(:, 3));
    if any(massNumbers <= 0)  % positive masses
        error('Non-positive isotope mass in row(s) %s.', strjoin(cellstr(num2str(find(massNumbers <= 0))), ', '))
    elseif any(abundances < 0 | abundances > 1)  % 0 <= abundnace <= 1
        error('Invalid isotope abundance in row(s) %s. Must lie between 0 and 1.', ...
            strjoin(cellstr(num2str(find(abundances < 0 | abundances > 1))), ', '))
    end

    % get the order of the supplied isotopeAbundance in the atomic weight data
    eleWtIsoAbIDs = zeros(size(isotopeAbundance, 1), 1);
    for j = 1:size(isotopeAbundance, 1)
        f = strcmp(isotopeAbundance{j, 1}, isotopes(:)) & isoMassNums(:) == isotopeAbundance{j, 2};
        if sum(f) == 1
            eleWtIsoAbIDs(j) = find(f);
        end
    end
    % error if the isotopes cannot be matched
    if ~all(eleWtIsoAbIDs)
        notMatched = find(eleWtIsoAbIDs == 0);
        errMsg = '';
        for j = 1:numel(notMatched)
            errMsg = [errMsg, ['\n#' num2str(notMatched(j)) '\t' ...
                isotopeAbundance{notMatched(j), 1} '\t' num2str(isotopeAbundance{notMatched(j), 2})]];
        end
        error(['%s', errMsg], 'The following isotopes could not be identified:')
    end

    % get the corresponding relative atomic mass
    relativeMasses = columnVector(isoWeights(eleWtIsoAbIDs));

    % get the user-defined atomic masses
    sumIsoEq1 = false(numel(eleWtIsoAb), 1);
    averageWeight = zeros(numel(eleWtIsoAb), 1);
    for j = 1:numel(eleWtIsoAb)
        sumIsoEq1(j) = abs(sum(abundances(ib == j)) - 1) < 1e-7;
        averageWeight(j) = relativeMasses(ib == j)' * abundances(ib == j);
    end
    if any(~sumIsoEq1)  % sum(abundance) = 1
        error('Invalid isotope abundances for element(s) %s. The sum of isotope abundances for each element must be equal to 1.', ...
            strjoin(eleWtIsoAb(~sumIsoEq1), ', '))
    end

    % update the list of weights
    elementWeights(idE) = averageWeight;
end

elementalWeightMatrix = [symbols, num2cell(elementWeights)];

author: opencobra

repo: cobratoolbox


text:
function pred = FDDL_pred(Y, D, CoefM, opts) % GC
    vgamma = opts.gamma;
    opts.max_iter = 100;
    [X, ~] = lasso_fista(Y, D, zeros(size(D,2), size(Y,2)), vgamma, opts);
    C = size(CoefM,2);
    % w = 0.5;
    E = zeros(C, size(Y,2));
    for c = 1: C
        Dc = get_block_col(D, c, opts.D_range);
        Xc = get_block_row(X, c, opts.D_range);
        R1 = Y - Dc*Xc;
        E1 = sum(R1.^2);
        R2 = X - repmat(CoefM(:, c), 1, size(Y,2));
        E2 = sum(R2.^2);
        E(c,:) = E1 + opts.weight*E2;
    end
    [~, pred] = min(E);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function boundary = p05_boundary_nearest ( m, n, point )

%*****************************************************************************80
%
%% P05_BOUNDARY_NEAREST returns a nearest boundary point in problem 05.
%
%  Discussion:
%
%    The given input point need not be inside the region.
%
%    In some cases, more than one boundary point may be "nearest",
%    but only one will be returned.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real BOUNDARY(M,N), points on the boundary
%    that are nearest to each point.
%
  center1 = [  0.0, 0.0 ];
  center2 = [ -0.4, 0.0 ];
  r1 = 1.0;
  r2 = 0.55;

  for j = 1 : n

    dist_min = r8_huge ( );
%
%  Distance to the semicircle #1.
%
    if ( center1(2) <= point(2,j) )

      norm = sqrt ( ( point(1,j) - center1(1) ).^2 ...
                  + ( point(2,j) - center1(2) ).^2 );

      if ( 0.0 == norm )

        pn(1) = center1(1) + r1 * sqrt ( 0.5 );
        pn(2) = center1(2) + r1 * sqrt ( 0.5 );

      else

        pn(1) = center1(1) + ( point(1,j) - center1(1) ) / norm;
        pn(2) = center1(2) + ( point(2,j) - center1(2) ) / norm;

      end

    elseif ( point(1,j) <= center1(1) )

      pn(1) =  center1(1) - r1;
      pn(2) =  center1(2);

    elseif ( center1(1) <= point(1,j) )

      pn(1) =  center1(1) + r1;
      pn(2) =  center1(2);

    end

    dist = sqrt ( ( point(1,j) - pn(1) ).^2 + ( point(2,j) - pn(2) ).^2 );

    if ( dist < dist_min )
      dist_min = dist;
      boundary(1:2,j) = pn(1:2)';
    end
%
%  Distance to semicircle #2.
%
    if ( center2(2) <= point(2,j) )

      norm = sqrt ( ( point(1,j) - center2(1) ).^2 ...
                  + ( point(2,j) - center2(2) ).^2 );

      if ( 0.0 == norm )

        pn(1) = center2(1) + r2 * sqrt ( 0.5 );
        pn(2) = center2(2) + r2 * sqrt ( 0.5 );

      else

        pn(1) = center2(1) + ( point(1,j) - center2(1) ) / norm;
        pn(2) = center2(2) + ( point(2,j) - center2(2) ) / norm;

      end

    elseif ( point(1,j) <= center2(1) )

      pn(1) =  center2(1) - r2;
      pn(2) =  center2(2);

    elseif ( center2(1) <= point(1,j) )

      pn(1) =  center2(1) + r2;
      pn(2) =  center2(2);

    end

    dist = sqrt ( ( point(1,j) - pn(1) ).^2 + ( point(2,j) - pn(2) ).^2 );

    if ( dist < dist_min )
      dist_min = dist;
      boundary(1:2,j) = pn(1:2)';
    end
%
%  Distance to line segment #1: (P1,P2).
%
    p1(1:2) = [ center2(1) - r1, center2(1) ];
    p2(1:2) = [ center2(1) - r2, center2(2) ];

    [ pn, dist, t ] = segment_point_near_2d ( p1, p2, point(1:2,j)' );

    if ( dist < dist_min )
      dist_min = dist;
      boundary(1:2,j) = pn(1:2)';
    end
%
%  Distance to line segment #2: (Q1,Q2).
%
    q1(1:2) = [ center2(1) + r2, center2(2) ];
    q2(1:2) = [ center1(1) + r1, center1(2) ];

    [ pn, dist, t ] = segment_point_near_2d ( q1, q2, point(1:2,j)' );

    if ( dist < dist_min )
      dist_min = dist;
      boundary(1:2,j) = pn(1:2)';
    end

  end

  return
end


text:
function Population = ArchiveUpdate(Population,N)
% Update archive

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Select feasible solutions
    fIndex     = all(Population.cons <= 0,2);
    Population = Population(fIndex);
    if isempty(Population)
        return
    else
        if size(Population.objs,2)==2
            %% Non-dominated sorting
            [FrontNo,~] = NDSort(Population.objs,1);
            Next = (FrontNo == 1);
            Population = Population(Next);
            if sum(Next) > N
                %% Calculate the crowding distance of each solution
                CrowdDis   = CrowdingDistance(Population.objs);
                [~,Rank]   = sort(CrowdDis,'descend');
                Population = Population(Rank(1:N));
            end
        else
            Population = Population(NDSort(Population.objs,1)==1);
            Population = Population(randperm(length(Population)));
            PCObj = Population.objs;
            nND   = length(Population);
            %% Population maintenance
            if length(Population) > N
                % Normalization
                fmax  = max(PCObj,[],1);
                fmin  = min(PCObj,[],1);
                PCObj = (PCObj-repmat(fmin,nND,1))./repmat(fmax-fmin,nND,1);
                % Determine the radius of the niche
                d  = pdist2(PCObj,PCObj);
                d(logical(eye(length(d)))) = inf;
                sd = sort(d,2);
                r  = mean(sd(:,min(3,size(sd,2))));
                R  = min(d./r,1);
                % Delete solution one by one
                while length(Population) > N
                    [~,worst]  = max(1-prod(R,2));
                    Population(worst)  = [];
                    R(worst,:) = [];
                    R(:,worst) = [];
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function title = p01_title ( )

%*****************************************************************************80
%
%% P01_TITLE returns the name of problem p01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Oscillatory';

  return
end


text:
function fx = toValues(disc, f, flag)
%TOVALUES   Convert a CHEBFUN to its TRIGSPEC discretization.
%   C = TOVALUES(DISC, F) converts the smooth chebfun F to coefficients C for
%   use by a TRIGSPEC discretization DISC.
%
%   C = TOVALUES(DISC, F, 1) converts the (perhaps piecewise smooth) chebfun F
%   to coefficients C for use by a TRIGSPEC discretization DISC.

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

if ( isnumeric(f) )
    fx = f;
    return
end

dim = disc.dimension;
tech = disc.returnTech();
tech = tech();
xx = trigpts(dim, f.domain);
fx = tech.vals2coeffs(f(xx));

end


text:
function varargout = GradList(varargin)
% GRADLIST MATLAB code for GradList.fig
%      GRADLIST, by itself, creates a new GRADLIST or raises the existing
%      singleton*.
%
%      H = GRADLIST returns the handle to a new GRADLIST or the handle to
%      the existing singleton*.
%
%      GRADLIST('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GRADLIST.M with the given input arguments.
%
%      GRADLIST('Property','Value',...) creates a new GRADLIST or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before GradList_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to GradList_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help GradList

% Last Modified by GUIDE v2.5 24-May-2013 14:12:48

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GradList_OpeningFcn, ...
                   'gui_OutputFcn',  @GradList_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before GradList_figure is made visible.
function GradList_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to GradList_figure (see VARARGIN)

handles.Simuh=varargin{1};
Category=dir([handles.Simuh.MRiLabPath filesep 'Config' filesep 'Grad']);
ind=1;
for i=1:length(Category)
    if ~strcmp(Category(i).name,'.') & ~strcmp(Category(i).name,'..') & Category(i).isdir
        handles.CategoryDir{ind,1}=Category(i).name;
        ind=ind+1;
    end
end
set(handles.Category_listbox,'String',handles.CategoryDir);

handles.accept = 0;
% Choose default command line output for GradList_figure
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes GradList_figure wait for user response (see UIRESUME)
% uiwait(handles.GradList_figure);

% --- Outputs from this function are returned to the command line.
function varargout = GradList_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on selection change in Grad_listbox.
function Grad_listbox_Callback(hObject, eventdata, handles)
% hObject    handle to Grad_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns Grad_listbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Grad_listbox

if ~isempty(get(handles.GradPath_text,'String'))
    return;
end

if strcmp(handles.GradList,'Empty')
    return;
end

handles.GradXMLFile=[handles.Simuh.MRiLabPath filesep ...
                    'Config' filesep ...
                    'Grad' filesep ...
                    handles.CategoryDir{get(handles.Category_listbox,'Value')} filesep ...
                    handles.GradList{get(handles.Grad_listbox,'Value')} filesep...
                    handles.GradList{get(handles.Grad_listbox,'Value')} '.xml'];

handles.GradXMLDir=[handles.Simuh.MRiLabPath filesep ...
                    'Config' filesep ...
                    'Grad' filesep ...
                    handles.CategoryDir{get(handles.Category_listbox,'Value')} filesep ...
                    handles.GradList{get(handles.Grad_listbox,'Value')}];

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function Grad_listbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Grad_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in Category_listbox.
function Category_listbox_Callback(hObject, eventdata, handles)
% hObject    handle to Category_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns Category_listbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Category_listbox

if ~isempty(get(handles.GradPath_text,'String'))
    return;
end

if strcmp(handles.CategoryDir,'Empty')
    return;
end

set(handles.Grad_listbox,'Enable','on');

Grad=dir([handles.Simuh.MRiLabPath filesep ...
          'Config' filesep ...
          'Grad' filesep ...
          handles.CategoryDir{get(handles.Category_listbox,'Value')}]);
ind=1;
handles.GradList=[];
for i=1:length(Grad)
    if ~strcmp(Grad(i).name,'.') & ~strcmp(Grad(i).name,'..')
        handles.GradList{ind,1}=Grad(i).name;
        ind=ind+1;
    end
end
if ind==1
    handles.GradList='Empty';
end
set(handles.Grad_listbox,'String',handles.GradList);
set(handles.Grad_listbox,'Value',1);

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function Category_listbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Category_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in Accept_pushbutton.
function Accept_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Accept_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isempty(get(handles.GradPath_text,'String'))
    try
        if strcmp(handles.GradList{get(handles.Grad_listbox,'Value')},'Empty')
            DoUpdateInfo(handles.Simuh,'No gradient is selected!');
            return;
        end
    catch me
        DoUpdateInfo(handles.Simuh,'No gradient is selected!');
        return;
    end
end

%----update gradient selection
handles.Simuh.GradXMLFile=handles.GradXMLFile;
handles.Simuh.GradXMLDir=handles.GradXMLDir;
Grad=DoParseXML(handles.GradXMLFile);
handles.Simuh.GradStruct=Grad;
Grad_text=uicontrol(handles.Simuh.GradSel_uipanel,'Style','text');
set(Grad_text,'FontWeight','bold','FontSize',10,'Position',[0,0,140,18],...
    'String', Grad.Attributes(1).Value,'TooltipString',handles.GradXMLFile);
%----end
guidata(handles.Simuh.SimuPanel_figure,handles.Simuh);
DoUpdateInfo(handles.Simuh,'Grad file was successfully loaded!');
handles.accept = 1;
guidata(hObject, handles);
close(handles.GradList_figure);

% --- Executes on button press in Cancel_pushbutton.
function Cancel_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Cancel_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
GradList_figure_CloseRequestFcn(hObject, eventdata, handles);

% --- Executes when user attempts to close GradList_figure.
function GradList_figure_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to GradList_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

DoDisableButton([],[],handles.Simuh);
if handles.accept == 0
    DoUpdateInfo(handles.Simuh,'Grad file loading was cancelled!');
end
% Hint: delete(hObject) closes the figure
delete(handles.GradList_figure);

% --- Executes on button press in GradPath_pushbutton.
function GradPath_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to GradPath_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[filename,pathname,filterindex]=uigetfile({'Grad*.xml','XML-files (*.xml)'},'MultiSelect','off');
if filename~=0
    handles.GradXMLFile=[pathname filename];
    handles.GradXMLDir=pathname(1:end-1);
    set(handles.GradPath_text,'String',[pathname filename]);
    % Add searchig path
    path(path,handles.GradXMLDir);

else
    errordlg('No Grad was loaded!');
    return;
end
guidata(hObject, handles);

% --- Executes on button press in Default_pushbutton.
function Default_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Default_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles.Simuh,'GradXMLFile')
    handles.Simuh=rmfield(handles.Simuh,'GradXMLFile');
    delete(get(handles.Simuh.GradSel_uipanel,'Children'));
    guidata(handles.Simuh.SimuPanel_figure,handles.Simuh);
end
DoUpdateInfo(handles.Simuh,'Use default Gradient!');
handles.accept = 1;
guidata(hObject, handles);
close(handles.GradList_figure);

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab