text:
function board = designboard(file)

%
% This program allows you to draw a board, based on some image (file). Each
% territory is drawn by clicking the coordinates of its vertices. On
% clicking the first vertex again, the patch is closed. To finish the
% territory, click again somewhere in the middle. This determines the
% position where the number of troops will be shown. You can ofcourse draw
% any kind of board, having arbitrary names and continents. However, to be
% fully compatible with the playrisk-program, the board should have the
% same continents as the original game: Asia, North-America, Europe,
% Africa, Australia and South-America. Moreover, the map should be fully
% interconnected. This means that it should be possible to travel from any
% territory to any other territory: no isolated patches. This program keeps
% track of the neighbouring patches, but oversea connections should be
% specified manually using the 'connect()' function. If you downloaded this
% file from the Mathworks Central File Exchange, there should be a
% predesigned board included meeting these requirements. An optional
% feature is the display of lines between two patches. This should also
% be done outside this program using the lines() function.
%

sz = [800, 600];
close(findobj('Name', 'Board'))
board = struct('Name',[],'Continent',[],'Neighbors',[],'xy',[],'TextPos',[], 'lines', []);
res = get(0, 'ScreenSize');
fig = figure('Name', 'Board', 'Position', res);
ax  = axes('Parent', fig, 'XLim', [0, sz(1)], 'YLim', [0, sz(2)]);

if nargin == 1 && exist(file, 'file')
    pic = imread(file);
    imshow(pic, 'Parent', ax, 'InitialMagnification', 'fit');
end

eps = 5;    % Welding treshold

%% MAKE BOARD
count = 1;
while true
    prompt = {'Name','Continent'};
    data = inputdlg(prompt, 'Name', 1, {'',''});
    if isempty(data)
        break
    end
    [x, y, nblist, tp] = getxy(board, eps);
    poly = patch(x, y, [0 0 0]);
    set(poly, 'EdgeColor', [1 0 0], 'Marker', '.', 'MarkerEdgeColor', [0 0 1]);
    ok = questdlg('OK?', 'Verify', 'Yes', 'No', 'Yes');
    if strcmp(ok, 'Yes')
        board(count) = struct('Name', data{1},          ...
                              'Continent', data{2},     ...
                              'Neighbors', nblist,      ...
                              'xy', [x,y],              ...
                              'TextPos', tp, ...
                              'lines', []);
        count = count + 1;
        set(poly, 'FaceColor', [.5 .5 .5])
    else
        delete(poly);
    end
end

%% SYNCHRONIZE NEIGHBORS
for i = 1:numel(board)
    for j = board(i).Neighbors
        if isempty(find(board(j).Neighbors == i, 1))
            board(j).Neighbors = [board(j).Neighbors, i];
        end
    end
end

%% FINALIZE
if isempty(board(1).Name)
    close(fig);
    return
end

for i = 1:numel(board)
    board(i).xy(:,2) =  sz(2) - board(i).xy(:,2);
    board(i).TextPos(2) = sz(2) - board(i).TextPos(2);
end
close(fig);

author: Sable

repo: mcbench-benchmarks


text:
function computeAllPrediction_batchVASARI_LGG(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeAllPrediction_batchVASARI_LGG(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes prediction performance estimation for a given
% feature set type, and for all model orders of all experiments with
% different degrees of freedom. See ref. [1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to where all experiments need to be
%                     performed.
%                     --> Ex: /myProject/WORKSPACE/VASARI
% 2. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 3. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 4. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (formal reference to come).
%               --> Ex: 'IALR'
% 5. nBatch: Number of parallel batch.
%            --> Ex: 8
% 6. matlabPATH: Full path to the MATLAB executable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% 7. seed: Random generator seed for reproducibility of experiments.
% -------------------------------------------------------------------------
% OUTPUTS: Prediction performance results are saved in a folder named 'RESULTS' in the
% corresponding folder of 'pathExperiments'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

% INITIALIZATON
time = 60; % Number of seconds to wait before checking if parallel computations are done
cd(pathExperiments), load('outcomes')
pathModels = fullfile(pwd,'MODELS'); mkdir('RESULTS'), cd('RESULTS'), pathResults = pwd;
mkdir('batchLog_Results'), cd('batchLog_Results'), pathBatch = pwd;
setNames = {'VASARI'};
[param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

% PRODUCE BATCH COMPUTATIONS
save('workspace','pathModels','pathResults','outcomes','param','maxOrder','nBoot','imbalance','seed'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:numel(param{i})
        fprintf(fid,['computeAllPrediction_VASARI_LGG(pathModels,pathResults,outcomes,param{',num2str(i),'}{',num2str(j),'},maxOrder,nBoot,imbalance,seed)\n']);
    end
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [exchangedMetsTable,debugMets,fixedModel,exchangeMetsL] = getExchangedMetaboliteTable(model,exRxnsOrdered,L)

[m,n]=size(model.S);

fixedModel = model;

debugMets=[];

exRxns = cell(m,1);
for i=1:m
    j = find(strcmp(model.rxns,['EX_' model.mets{i}]));
    if isempty(j)
        exRxns{i}='';
    else
        if model.S(i,j)~=0
            exRxns{i} = model.rxns{j};
        else
            debugMets = [debugMets;model.mets{i}];

            warning(['missing exchange reaction coefficient for ' model.mets{i}])
            disp([model.mets{i},model.metNames{i},model.rxns(contains(model.rxns,model.mets{i})), model.mets{model.S(:,j)~=0}])
            fixedModel.S(i,:)=0;
            fixedModel.S(i,j)=-1;
        end
    end
end
exchangedMetsTable = table(model.mets,exRxns,(1:m)','VariableNames',{'mets','rns','ind'});

exMetsOrdered=cell(length(exRxnsOrdered),1);
if ~isempty(exRxnsOrdered)
    for k=1:length(exRxnsOrdered)
        exMetsOrdered{k} = strrep(exRxnsOrdered{k},'EX_','');
    end
    exchangedMetsTable = mapAontoB(exchangedMetsTable.mets,exMetsOrdered,exchangedMetsTable);
    bool = exchangedMetsTable.ind==0;
    if any(bool)
        L = [L,zeros(size(L,1),max(exchangedMetsTable.ind)-size(L,2))];
    end
    exchangedMetsTable.ind(bool)=size(L,2);
    %hack to deal with L from different model

    exchangeMetsL = L(:,exchangedMetsTable.ind);

    exchangeMetsL = exchangeMetsL(any(exchangeMetsL,2),:);

end

figure
spy(exchangeMetsL)
ax = gca;
outerpos = ax.OuterPosition;
ti = ax.TightInset;
left = outerpos(1) + ti(1);
bottom = outerpos(2) + ti(2);
ax_width = outerpos(3) - ti(1) - ti(3);
ax_height = outerpos(4) - ti(2) - ti(4);
ax.Position = [left bottom ax_width ax_height];

% did not work - ordering not diagonal.
% table(model.mets((getCorrespondingRows(model.S, true(size(model.S,1),1), ~model.SConsistentRxnBool, 'inclusive'))~=0),model.rxns(~model.SConsistentRxnBool),'VariableNames',{'mets','rns'});

author: opencobra

repo: cobratoolbox


text:
function [U,c] = MgnCalibration(X)
% performs magnetometer calibration from a set of data
% using Merayo technique with a non iterative algoritm
% J.Merayo et al. "Scalar calibration of vector magnemoters"
% Meas. Sci. Technol. 11 (2000) 120-132.
%
%   X      : a Nx3 (or 3xN) data matrix
%              each row (columns) contains x, y, z measurements
%              N must be such that the data set describes
%              as completely as possible the 3D space
%              In any case N > 10
%
%    The calibration tries to find the best 3D ellipsoid that fits the data set
%    and returns the parameters of this ellipsoid
%
%    U     :  shape ellipsoid parameter, (3x3) upper triangular matrix
%    c      : ellipsoid center, (3x1) vector
%
%    Ellipsoid equation : (v-c)'*(U'*U)(v-c) = 1
%    with v a rough triaxes magnetometer  measurement
%
%    calibrated measurement w = U*(v-c)
%
%   author : Alain Barraud, Suzanne Lesecq 2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[N,m] = size(X);
if m>3&&N==3,X = X';N = m;m = 3;end;%check that X is not transposed
if N<=10,U = [];c = [];return;end;%not enough data no calibration !!
% write  the ellipsoid equation as D*p=0
% the best parameter is the solution of min||D*p|| with ||p||=1;
% form D matrix from X measurements
x = X(:,1); y = X(:,2); z = X(:,3);
D = [x.^2, y.^2, z.^2, x.*y, x.*z, y.*z, x, y, z, ones(N,1)];
D=triu(qr(D));%avoids to compute the svd of a large matrix
[U,S,V] = svd(D);%because usually N may be very large
p = V(:,end);if p(1)<0,p =-p;end;
% the following matrix A(p) must be positive definite
% The optimization done by svd does not include such a constraint
% With "good" data the constraint is allways satisfied
% With too poor data A may fail to be positive definite
% In this case the calibration fails
%
A = [p(1) p(4)/2 p(5)/2;
       p(4)/2 p(2) p(6)/2;
       p(5)/2 p(6)/2 p(3)];
[U,ok] = fchol(m,A);
if ~ok,U = [];c = [];return;end%calibration fails too poor data!!
b = [p(7);p(8);p(9)];
v = Utsolve(U,b/2,m);
d = p(10);
s = 1/sqrt(v*v'-d);
c =-Usolve(U,v,m)';%ellipsoid center
U = s*U;%shape ellipsoid parameter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,ok] = fchol(n,A)
% performs Cholesky factoristation
A(1,1:n) = A(1,1:n)/sqrt(A(1,1));
A(2:n,1) = 0;
for j=2:n
  A(j,j:n) = A(j,j:n) - A(1:j-1,j)'*A(1:j-1,j:n);
  if A(j,j)<=0,ok=0;break;end%A is not positive definite
  A(j,j:n) = A(j,j:n)/sqrt(A(j,j));
  A(j+1:n,j) = 0;
end
ok=1;
function x=Utsolve(U,b,n)
% solves U'*x=b
x(1) = b(1)/U(1,1);
for k=2:n
    x(k) = (b(k)-x(1:k-1)*U(1:k-1,k))/U(k,k);
end
function x=Usolve(U,b,n)
% solves U*x=b
x(n) = b(n)/U(n,n);
for k=n-1:-1:1
    x(k) = (b(k)-U(k,k+1:n)*x(k+1:n)')/U(k,k);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function varargout = showXSlice(img, sliceIndex)
%SHOWXSLICE Show YZ slice of a 3D image
%
%   showXSlice(IMG, INDEX)
%   Display the given slice as a 3D planar image. INDEX is the slice index,
%   between 1 and stackSize(img, 1).
%
%   Example
%   % Display 3D orthoslices of a humain head
%   img = analyze75read(analyze75info('brainMRI.hdr'));
%   figure(1); clf; hold on;
%   showZSlice(img, 13);
%   showXSlice(img, 60);
%   showYSlice(img, 80);
%   axis equal
%   xlabel('x'); ylabel('y'); zlabel('z');
%
%   See also
%   showYSlice, showZSlice, getSlice

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-06-30,    using Matlab 7.9.0.529 (R2009b)
% http://www.pfl-cepia.inra.fr/index.php?page=slicer
% Copyright 2010 INRA - Cepia Software Platform.

%% Extract image info

dim = stackSize(img);

% compute voxel positions
lx = 1:dim(1);

% position vectors of voxel corners
vy = ((0:dim(2)) + .5);
vz = ((0:dim(3)) + .5);

% global parameters for surface display
params = {'facecolor', 'texturemap', 'edgecolor', 'none'};

% compute position of voxel vertices in 3D space
[yz_y, yz_z] = meshgrid(vy, vz);
yz_x = ones(size(yz_y)) * lx(sliceIndex);

% extract slice in x direction
slice = stackSlice(img, 'x', sliceIndex);

% eventually converts to uint8, rescaling data between 0 and max value
if ~isa(slice, 'uint8')
    slice = double(slice);
    slice = uint8(slice * 255 / max(slice(:)));
end

% convert grayscale to rgb (needed by 'surface' function)
if length(size(slice)) == 2
    slice = repmat(slice, [1 1 3]);
end

% display voxel values in appropriate reference space
hs = surface(yz_x, yz_y, yz_z, slice, params{:});

%% process output arguments

if nargout > 0
    varargout = {hs};
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function phy = simplex_unit_to_general ( dim_num, point_num, t, ref )

%*****************************************************************************80
%
%% SIMPLEX_UNIT_TO_GENERAL maps the unit simplex to a general simplex.
%
%  Discussion:
%
%    Given that the unit simplex has been mapped to a general simplex
%    with vertices T, compute the images in T, under the same linear
%    mapping, of points whose coordinates in the unit simplex are REF.
%
%    The vertices of the unit simplex are listed as suggested in the
%    following:
%
%      (0,0,0,...,0)
%      (1,0,0,...,0)
%      (0,1,0,...,0)
%      (0,0,1,...,0)
%      (...........)
%      (0,0,0,...,1)
%
%    Thanks to Andrei ("spiritualworlds") for pointing out a mistake in the
%    previous implementation of this routine, 02 March 2008.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points to transform.
%
%    Input, real T(DIM_NUM,DIM_NUM+1), the vertices of the
%    general simplex.
%
%    Input, real REF(DIM_NUM,POINT_NUM), points in the
%    reference triangle.
%
%    Output, real PHY(DIM_NUM,POINT_NUM), corresponding points
%    in the physical triangle.
%

%
%  The image of each point is initially the image of the origin.
%
%  Insofar as the pre-image differs from the origin in a given vertex
%  direction, add that proportion of the difference between the images
%  of the origin and the vertex.
%
  for dim = 1 : dim_num

    phy(dim,1:point_num) = t(dim,1);

    for vertex = 2 : dim_num + 1

      phy(dim,1:point_num) = phy(dim,1:point_num) ...
        + ( t(dim,vertex) - t(dim,1) ) * ref(vertex-1,1:point_num);

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g=BaylissLinearTapering(sidelobedB,N,xPoints,a)
%%BAYLISSLINEARTAPERING The Bayliss tapering is a set of complex amplitude
%          weights for a continuous linear (narrowband) aperture antenna
%          that will form a difference beam (odd symmetry about an axis)
%          and hold the four closest sidelobes to a desired level. Such a
%          tapering can be discretized and applied to the elements in a
%          circular phased array (An array of antenna elements can be
%          viewed as a discrete approximation to a continuous aperture).
%          This function will provide the tapering values at a set of
%          discrete points given by xPoints (the origin is taken to be the
%          center of the aperture). The radius of the aperture can either
%          be provided or is taken as value of the farthest point provided.
%
%INPUTS: sidelobedB The number of decibels of the ratio of the close-in
%             sidelobe voltages to the main lobe voltage. This must be a
%             negative number. A typical value is -30.
%           N The Bayliss tapering is computed using a certain number of
%             terms. Using too many terms can be undesirable as edge
%             illumination increases, as noted in [1]. If this parameter is
%             omitted or an empty matrix is passed, then the default of 17
%             is used. In [1], it is suggested that N be chosen to be
%             <2*a/lambda, where a is the radius of the aperture and
%             lambda the wavelength.
%    xyPoints An NX1 or 1XN set of N points at which the tapering values
%             should be evaluated. The center of the aperture is taken to
%             be the origin.
%           a The radius of the aperture. Tapering weights for points in
%             xPoints outside of the aperture are taken to be zero. If
%             this parameter is omitted or an empty matrix is passed, then
%             the radius is taken to be the distance of the farthest point
%             from the origin in xPoints.
%
%OUTPUTS: g The NX1 set of discretized Bayliss tapering values evaluated at
%           the points given in xPoints. All Bayliss tapering values are
%           imaginary. The coefficients are not normalized.
%
%This function implements the algorithm given in the appendix of in [1]
%using the polynomial interpolation values in the table below Figure 4.
%This approximation means that low sidelobe patterns (-45 dB and below)
%will not have good fidelity sidelobes.
%
%EXAMPLE 1:
%Here, we evaluate the tapering values for 30dB down on a fine grid of
%points to plot what the amplitude and phase of the tapering weights looks
%like.
% numPoints=300;
% xPoints=linspace(-1,1,numPoints);
% %The Bayliss tapering weights, evaluated across the aperture.
% a=1;%Aperture radius=1.
% %gTest=BaylissLinear();
% gBayliss=BaylissLinearTapering(-30,17,xPoints,a);
%
% figure(1)
% clf
% plot(xPoints,abs(gBayliss),'-b','linewidth',2)
% h1=xlabel('x');
% h2=ylabel('Imaginary Weight');
% title('Bayliss Tapering Amplitude')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
% figure(2)
% clf
% plot(xPoints,angle(gBayliss),'-b','linewidth',2)
% h1=xlabel('x');
% h2=ylabel('Imaginary Weight');
% title('Bayliss Tapering Phase')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%Here, we consider the array response when using tapering values for 30dB
%sidelobes on a linear array with lambda/2 spacing between elements.
% N=17;
% sidelobedB=-30;
% Nx=41;%There are 2*Nx+1 points total.
% %Generate points symmetric about the origin.
% xPoints=(-(Nx-1)/2):1/2:((Nx-1)/2);
% g=BaylissLinearTapering(sidelobedB,N,xPoints);
%
% %Now, display the response with the tapering
% T=diag(g);
% [Rsp,U]=standardUVBeamPattern(T,xPoints,'NormRealVal');
%
% figure(2)
% clf
% plot(U,Rsp,'-b','LineWidth',2);
% h1=xlabel('u');
% h2=ylabel('Array Response');
% title('Bayliss Weighted Array Response')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] E. T. Bayliss, "Design of monopulse antenna difference patterns with
%    low sidelobes," The Bell System Technical Journal, vol. 47, no. 5, pp.
%    623-650, May-Jun. 1968.
%
%August 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(N))
   N=17;
end

%Defined below Equation 41 in [1].
mu=((0:(N-1))+1/2)';

%This holds the coefficients for the interpolating polynomials given below
%Figure 4 in [1]. The polynomials all take the desired sidelobe level in
%decibels as an input parameter. The first row is for the term A, which
%is a translation of the SNR parameter to a parameter in the paper. The
%next four rows are xi_1 to xi_4, which are the locations of the first four
%zeroes in the modified pattern. The final row is for p_0, which is related
%to the point at which the peak of the asymptotic difference pattern=1.
polyCoeffTable=[0.30387530,-0.05042922,-0.00027989,-0.00000343,-0.00000002;
                0.98583020,-0.03338850, 0.00014064, 0.00000190, 0.00000001;
                2.00337487,-0.01141548, 0.00041590, 0.00000373, 0.00000001;
                3.00636321,-0.00683394, 0.00029281, 0.00000161, 0;
                4.00518423,-0.00501795, 0.00021735, 0.00000088, 0;
                0.47972120,-0.01456692,-0.00018739,-0.00000218,-0.00000001];

A  =polyCoeffTable(1,1)+sidelobedB*(polyCoeffTable(1,2)+sidelobedB*(polyCoeffTable(1,3)+sidelobedB*(polyCoeffTable(1,4)+sidelobedB*polyCoeffTable(1,5))));
xi1=polyCoeffTable(2,1)+sidelobedB*(polyCoeffTable(2,2)+sidelobedB*(polyCoeffTable(2,3)+sidelobedB*(polyCoeffTable(2,4)+sidelobedB*polyCoeffTable(2,5))));
xi2=polyCoeffTable(3,1)+sidelobedB*(polyCoeffTable(3,2)+sidelobedB*(polyCoeffTable(3,3)+sidelobedB*(polyCoeffTable(3,4)+sidelobedB*polyCoeffTable(3,5))));
xi3=polyCoeffTable(4,1)+sidelobedB*(polyCoeffTable(4,2)+sidelobedB*(polyCoeffTable(4,3)+sidelobedB*(polyCoeffTable(4,4)+sidelobedB*polyCoeffTable(4,5))));
xi4=polyCoeffTable(5,1)+sidelobedB*(polyCoeffTable(5,2)+sidelobedB*(polyCoeffTable(5,3)+sidelobedB*(polyCoeffTable(5,4)+sidelobedB*polyCoeffTable(5,5))));
%p0 =polyCoeffTable(6,1)+sldelobedB*(polyCoeffTable(6,2)+sldelobedB*(polyCoeffTable(6,3)+sldelobedB*(polyCoeffTable(6,4)+sldelobedB*polyCoeffTable(6,5))));

Z=zeros(N+1,1);
%Equation 15
Z(1)=0;%The Z(0) term
%Now, the moved zeros
Z(2)=xi1;
Z(3)=xi2;
Z(4)=xi3;
Z(5)=xi4;
for k=5:N
    %The location of the non-moved zeros as given by Equation 13.
    Z(k+1)=sqrt(A^2+k^2);
end

sigma=(N+1/2)/Z(N+1);

B=zeros(N,1);
for m=0:(N-1)
    %This loop implements Equation 47.
    n=1:(N-1);
    num=prod(1-((m+1/2)./(sigma*Z(1+n))).^2);
    l=[0:(m-1),(m+1):(N-1)];
    denom=prod(1-((m+1/2)./(l+1/2)).^2);
    %We just use C=1.
    B(m+1)=1/(2*1j)*(-1)^m*(m-1/2)^2*num/denom;
end

numPoints=length(xPoints);

if(nargin<4||isempty(a))
    %The maximum distance from the origin to a point is taken to be the
    %radius of the aperture.
    a=sqrt(max(sum(xPoints.^2,1)));
end

%The loop below implements Equation 41 in [1].
g=zeros(numPoints,1);
for curPoint=1:numPoints
    if(abs(xPoints(curPoint))<=a)
        %The normalized radius at this point.
        p=pi*xPoints(curPoint)/a;
        g(curPoint)=g(curPoint)+sum(B.*sin(mu*p));
    end
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 [parameters, fitnessScore, exitFlag, newModel]=fitParameters(model,xRxns,xValues,rxnsToFit,valuesToFit,parameterPositions,fitToRatio,initialGuess,plotFitting)
% fitParameters
%   Fits parameters such as maintenance ATP by quadratic programming
%
%   model                 a model structure
%   xRxns                 cell array with the IDs of the reactions that will be
%                         fixed for each data point
%   xValues               matrix with the corresponding values for each
%                         xRxns (columns are reactions)
%   rxnsToFit             cell array with the IDs of reactions that will be fitted to
%   valuesToFit           matrix with the corresponding values for each
%                         rxnsToFit (columns are reactions)
%   parameterPositions    stucture that determines where the parameters are in the
%                         stoichiometric matrix. Contains the fields:
%   	position          cell array of vectors where each element contains
%                         the positions in the S-matrix for that parameter
%   	isNegative        cell array of vectors where the elements are true
%                         if that position should be the negative of the
%                         fitted value (to differentiate between
%                         production/consumption)
%	fitToRatio            if the ratio of simulated to measured values should
%                         be fitted instead of the absolute value. Used to prevent
%                         large fluxes from having too large impact (opt,
%                         default true)
%   initialGuess          initial guess of the parameters (opt)
%   plotFitting           true if the resulting fitting should be plotted
%                         (opt, default false)
%
%   parameters            fitted parameters in the same order as in
%                         parameterPositions
%   fitnessScore          the correponding residual sum of squares
%   newModel              updated model structure with the fitted parameters
%
%   Usage: [parameters, fitnessScore, exitFlag, newModel]=fitParameters(model,...
%           xRxns,xValues,rxnsToFit,valuesToFit,parameterPositions,fitToRatio,...
%           initialGuess,plotFitting)

if nargin<7
    fitToRatio=true;
end
if nargin<8
    initialGuess=ones(numel(parameterPositions.position),1);
end
if isempty(initialGuess)
    initialGuess=ones(numel(parameterPositions.position),1);
end
if nargin<9
    plotFitting=false;
end

xRxns=convertCharArray(xRxns);
rxnsToFit=convertCharArray(rxnsToFit);

%Find the indexes of reactions that will be fitted
[I, rxnsToFitIndexes]=ismember(rxnsToFit,model.rxns);

if ~all(I)
    EM='Could not find all reactions in rxnsToFit';
    dispEM(EM);
end

%Find the indexes of reactions that will be used for constraints.
[I, xRxnsIndexes]=ismember(xRxns,model.rxns);

if ~all(I)
    EM='Could not find all reactions in xRxns';
    dispEM(EM);
end

[parameters, fitnessScore, exitFlag]=fminsearch(@(parameters) getRSS(parameters,model,xRxnsIndexes,xValues,rxnsToFitIndexes,valuesToFit,parameterPositions,fitToRatio),initialGuess);

parameters=abs(parameters);

if plotFitting==true
    %Set the resulting parameters
    [~, resultingFluxes, newModel]=getRSS(parameters,model,xRxnsIndexes,xValues,rxnsToFitIndexes,valuesToFit,parameterPositions,true);
    plot(xValues,valuesToFit,'o',xValues,resultingFluxes,'-*');
end
end

function [rss, resultingFluxes, newModel]=getRSS(parameters,model,xRxnsIndexes,xValues,rxnsToFitIndexes,valuesToFit,parameterPositions,fitToRatio)
parameters=abs(parameters);

%Set the parameters at the positions specified in parameterPositions
for i=1:numel(parameterPositions.position)
    %Set positive
    model.S(parameterPositions.position{i}(parameterPositions.isNegative{i}==false))=parameters(i);

    %Set negative
    model.S(parameterPositions.position{i}(parameterPositions.isNegative{i}==true))=parameters(i)*-1;
end

%Also return an updated model
newModel=model;

%Loop through each data point, set xRxns to xValues and calculate the sum
%of squares for the rxnsToFit
rss=0;
resultingFluxes=[];
for i=1:size(xValues,1)
    %Fix for more xRxns!
    model.lb(xRxnsIndexes)=xValues(i,:);
    model.ub(xRxnsIndexes)=xValues(i);

    sol=solveLP(model);

    %Calculate the rss
    if fitToRatio==false
        rs=sol.x(rxnsToFitIndexes)'-valuesToFit(i,:);
    else
        rs=sol.x(rxnsToFitIndexes)'./valuesToFit(i,:)-ones(1,size(valuesToFit,2));
    end
    rss=rss+rs*rs';
    resultingFluxes=[resultingFluxes sol.x(rxnsToFitIndexes)];
end
end

author: SysBioChalmers

repo: RAVEN


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%
v=sigmoid(X*theta);
for i=1:m
	if v(i,1)>=0.5
		p(i,1)=1;
	else
		p(i,1)=0;
end

% =========================================================================

end

author: loserChen

repo: Coursera-MachineLearning


text:
function t = alldifferent(X)
% alldifferent (overloaded)

t = alldifferent(sdpvar(X));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function region = p06_region ( )

%*****************************************************************************80
%
%% P06_REGION returns the name of the integration region for problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end


text:
function []= mdp_verbose()

% mdp_verbose   Ask for running resolution functions of the MDP Toolbox
%               in verbose mode

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

global mdp_VERBOSE;

mdp_VERBOSE = 1;


text:
%CCG - Compute multiple cross- and auto-correlograms
%
%  USAGE
%
%    [ccg,t] = CCG(times,groups,<options>)
%
%    times          times of all events
%                   NOTE: spiketimes in SECONDS.
%    groups         group IDs for each event in time list (should be
%                   integers 1:nGroups)
%
%   alternate buzcode usage:
%   times           {Ncells} array of [Nspikes]  spiketimes for each cell.
%   groups          []
%
%                   spikes = bz_GetSpikes
%                   [ccg,t] = CCG(spikes.times,[])
%
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'binSize'     bin size in s (default = 0.01)
%     'duration'    duration in s of each xcorrelogram (default = 2)
%     'norm'        normalization of the CCG, 'counts' or 'rate' (DL added 8/1/17)
%                   'counts' gives raw event/spike count,
%                   'rate' returns CCG in units of spks/second (default: counts)
%    =========================================================================
%
%
%  OUTPUT
%   ccg     [t x ngroups x ngroups] matrix where ccg(t,i,j) is the
%           number (or rate) of events of group j at time lag t with
%           respect to reference events from group i
%   t       time lag vector (units: seonds)
%
%  SEE
%
%    See also ShortTimeCCG.

% Copyright (C) 2012 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function [ccg,t] = CCG(times,groups,varargin)

% Default values
duration = 2;
binSize = 0.01;
Fs = 1/20000;
normtype = 'counts';

% Option for spike times to be in {Ncells} array of spiketimes DL2017
if iscell(times) && isempty(groups)
    numcells = length(times);
    for cc = 1:numcells
        groups{cc}=cc.*ones(size(times{cc}));
    end
    times = cat(1,times{:}); groups = cat(1,groups{:});
end

%Sort
[times,sortidx] = sort(times);
groups = groups(sortidx);

% Check parameters
if nargin < 2,
  error('Incorrect number of parameters (type ''help CCG'' for details).');
end
%if ~isdvector(times),
%	error('Parameter ''times'' is not a real-valued vector (type ''help CCG'' for details).');
%end
if ~isdscalar(groups) && ~isdvector(groups),
	error('Parameter ''groups'' is not a real-valued scalar or vector (type ''help CCG'' for details).');
end
if ~isdscalar(groups) && length(times) ~= length(groups),
	error('Parameters ''times'' and ''groups'' have different lengths (type ''help CCG'' for details).');
end
groups = groups(:);
times = times(:);

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help CCG'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'binsize',
			binSize = varargin{i+1};
			%if ~isdscalar(binSize,'>0'),
		%		error('Incorrect value for property ''binSize'' (type ''help CCG'' for details).');
	%		end
		case 'duration',
			duration = varargin{i+1};
			if ~isdscalar(duration,'>0'),
				error('Incorrect value for property ''duration'' (type ''help CCG'' for details).');
            end

            case 'Fs',
			Fs = varargin{i+1};
			if ~isdscalar(Fs,'>0'),
				error('Incorrect value for property ''Fs'' (type ''help CCG'' for details).');
            end
        case 'norm'
            normtype = varargin{i+1};

  end
end

% Number of groups, number of bins, etc.
if length(groups) == 1,
	groups = ones(length(times),1);
	nGroups = 1;
else
	nGroups = max(unique(groups));
end

halfBins = round(duration/binSize/2);
nBins = 2*halfBins+1;
t = (-halfBins:halfBins)'*binSize;
times = round(times/Fs);
binSize_Fs = round(binSize/Fs);
if length(times) <= 1,
	% ---- MODIFIED BY EWS, 1/2/2014 ----
    % *** Use unsigned integer format to save memory ***
    ccg = uint16(zeros(nBins,nGroups,nGroups));
    % -----------------------------------
    return
end

% Compute CCGs
nEvents = length(times);
%

counts = double(CCGHeart(times,uint32(groups),binSize_Fs,uint32(halfBins)));
% -----------------------------------
%
% Reshape the results
n = max(groups);
counts = reshape(counts,[nBins n n]);

if n < nGroups,
	counts(nBins,nGroups,nGroups) = 0;
end

%Rate normalization: counts/numREFspikes/dt to put in units of spikes/s. DL
switch normtype
    case 'rate'
        for gg = 1:nGroups
            numREFspikes = sum(groups==gg);%number of reference events for group
            counts(:,gg,:) = counts(:,gg,:)./numREFspikes./binSize;
        end
end

ccg = flipud(counts);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ res, cmd_out ] = Dropbox( action )
% Function to start and kill dropbox from MATLAB
%
% Syntax:   [ res, cmd_out ] = Dropbox( action )
%
% Inputs:
% 	action - Action to perform. Either 'start' or 'kill.
%
% Outputs:
% 	res - Command exit status
% 	cmd_out - Output of the operating system command
%
% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2017
% Date: 3 October 2015
% Revision: 0.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch lower(action)
    case 'start'
        [res, cmd_out] = dos('start "" "%APPDATA%\Dropbox\bin\Dropbox.exe" & exit');
        disp('... Started Dropbox');

    case 'kill'
        [res, cmd_out] = dos('taskkill /F /IM Dropbox.exe');
        if res==0, disp('Terminated Dropbox ...'); end

    otherwise
        error('Dropbox action not supported...');
end

end

author: jtkim-kaist

repo: Speech-enhancement


text:
function y = sigmoid(x,c,a)
% sigmoid evaluates a simple sigmoid function along x:
%
%         ______1________
%     y =        -a(x-c)
%          1 + e^
%
%% Syntax
%
% y = sigmoid(x)
% y = sigmoid(x,c)
% y = sigmoid(x,c,a)
%
%% Description
%
% y = sigmoid(x) generates a sigmoid function along x.
%
% y = sigmoid(x,c) makes a sigmoid that scaled from zero to one, where
% c corresponds to the x value where y = 0.5. If c is not specified, a
% default value of c = 0 is assumed.
%
% y = sigmoid(x,c,a) specifies a, the rate of change. If a is close to
% zero, the sigmoid function will be gradual. If a is large, the sigmoid
% function will have a steep or sharp transition. If a is negative, the
% sigmoid will go from 1 to zero. A default value of a = 1 is assumed if
% a is not declared.
%
%% Example 1
% A simple sigmoid:
%
% x = -10:.01:10;
% plot(x,sigmoid(x))
%
%% Example 2:
% Make a sigmoid function along x = 1 to 100, such that y(x=60) = 0.5:
%
% x = 1:100;
% y = sigmoid(x,60);
% plot(x,y,'b','linewidth',2)
%
% Now do the same thing as above, but make the transition more gradual:
%
% y2 = sigmoid(x,60,0.1);
% hold on
% plot(x,y2,'r','linewidth',2)
% legend('default a = 1','a = 1/10')
% legend boxoff
%
%% Author Info:
% Chad Greene, May 28, 2015.
% http://www.chadagreene.com

%% Parse Inputs:

narginchk(1,3)
if nargin<3
    a = 1;
else
    assert(isscalar(a)==1,'a must be a scalar.')
end

if nargin<2
    c = 0;
else
    assert(isscalar(c)==1,'c must be a scalar.')
end

%% Perform mathematics:

y = 1./(1 + exp(-a.*(x-c)));


text:
function [ edge_data, node_num2 ] = tet_mesh_order4_to_order10_size ( ...
  tetra_num, tetra_node1, node_num1 )

%*****************************************************************************80
%
%% TET_MESH_ORDER4_TO_ORDER10_SIZE sizes a quadratic tet mesh from a linear one.
%
%  Discussion:
%
%    A quadratic (10 node) tet mesh can be derived from a linear
%    (4 node) tet mesh by interpolating nodes at the midpoint of
%    every edge of the mesh.
%
%    The mesh is described indirectly, as the sum of individual
%    tetrahedrons.  A single physical edge may be a logical edge of
%    any number of tetrahedrons.  It is important, however, that a
%    new node be created exactly once for each edge, assigned an index,
%    and associated with every tetrahedron that shares this edge.
%
%    This routine handles that problem.
%
%    The primary amount of work occurs in sorting a list of 6 * TETRA_NUM
%    data items, one item for every edge of every tetrahedron.  Each
%    data item records, for a given tetrahedron edge, the global indices
%    of the two endpoints, the local indices of the two endpoints,
%    and the index of the tetrahedron.
%
%    Through careful sorting, it is possible to arrange this data in
%    a way that allows the proper generation of the interpolated nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TETRA_NUM, the number of tetrahedrons in the
%    linear mesh.
%
%    Input, integer TETRA_NODE1(4,TETRA_NUM), the nodes that make up
%    each tetrahedron in the linear mesh.
%
%    Input, integer NODE_NUM1, the number of nodes for the linear mesh.
%
%    Output, integer EDGE_DATA(5,6*TETRA_NUM), edge data.
%
%    Output, integer NODE_NUM2, the number of nodes for the quadratic mesh.
%

%
%  Step 1.
%  From the list of nodes for tetrahedron T, of the form: (I,J,K,L)
%  construct the six edge relations:
%
%    (I,J,1,2,T)
%    (I,K,1,3,T)
%    (I,L,1,4,T)
%    (J,K,2,3,T)
%    (J,L,2,4,T)
%    (K,L,3,4,T)
%
%  In order to make matching easier, we reorder each pair of nodes
%  into ascending order.
%
  for tetra = 1 : tetra_num

    i = tetra_node1(1,tetra);
    j = tetra_node1(2,tetra);
    k = tetra_node1(3,tetra);
    l = tetra_node1(4,tetra);

    [ a, b ] = i4i4_sort_a ( i, j );

    edge_data(1:5,6*(tetra-1)+1) = [ a, b, 1, 2, tetra ]';

    [ a, b ] = i4i4_sort_a ( i, k );

    edge_data(1:5,6*(tetra-1)+2) = [ a, b, 1, 3, tetra ]';

    [ a, b ] = i4i4_sort_a ( i, l );

    edge_data(1:5,6*(tetra-1)+3) = [ a, b, 1, 4, tetra ]';

    [ a, b ] = i4i4_sort_a ( j, k );

    edge_data(1:5,6*(tetra-1)+4) = [ a, b, 2, 3, tetra ]';

    [ a, b ] = i4i4_sort_a ( j, l );

    edge_data(1:5,6*(tetra-1)+5) = [ a, b, 2, 4, tetra ]';

    [ a, b ] = i4i4_sort_a ( k, l );

    edge_data(1:5,6*(tetra-1)+6) = [ a, b, 3, 4, tetra ]';

  end
%
%  Step 2. Perform an ascending dictionary sort on the neighbor relations.
%  We only intend to sort on rows 1:2; the routine we call here
%  sorts on the full column but that won't hurt us.
%
%  What we need is to find all cases where tetrahedrons share an edge.
%  By sorting the columns of the EDGE_DATA array, we will put shared edges
%  next to each other.
%
  edge_data = i4col_sort_a ( 5, 6*tetra_num, edge_data );
%
%  Step 3. All the tetrahedrons which share an edge show up as consecutive
%  columns with identical first two entries.  Figure out how many new
%  nodes there are, and allocate space for their coordinates.
%
  node_num2 = node_num1;

  n1_old = -1;
  n2_old = -1;

  for edge = 1 : 6 * tetra_num
    n1 = edge_data(1,edge);
    n2 = edge_data(2,edge);
    if ( n1 ~= n1_old | n2 ~= n2_old )
      node_num2 = node_num2 + 1;
      n1_old = n1;
      n2_old = n2;
    end
  end

  return
end


text:
classdef prtClusterKmodes < prtCluster %prtClass %prtAction %should extent prtCluster
    %
    %   Example:
    %
    %   ds = prtDataGenMoon;                 % Load a prtDataSet
    %   clusterAlgo = prtClusterKmodes;      % Create a prtClusterKmeans object
    %   clusterAlgo.nClusters = 3;           % Set the number of desired clusters
    %
    %   % Set the internal decision rule to be MAP. Not required for
    %   % clustering, but necessary to plot the results.
    %   clusterAlgo.internalDecider = prtDecisionMap;
    %   clusterAlgo = clusterAlgo.train(ds); % Train the cluster algorithm
    %   plot(clusterAlgo);                   % Plot the results
    %
    %   See also prtCluster, prtClusterKMeans, prtClusterGmm

% Copyright (c) 2014 CoVar Applied Technologieshnb
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.

    properties (SetAccess=private)
        name = 'K-Modes Clustering' % K-Means Clustering
        nameAbbreviation = 'K-ModesCluster' % K-MeansCluster
    end

    properties
        nClusters = 3;  % The number of clusters to find

        kdeRv = prtRvKde;

        handleEmptyClusters = 'remove';  % Action to take when an empty cluster occurs
        distanceMetricFn = @prtDistanceEuclidean;  % The distance metric; should be a function like D = prtDistanceEuclidean(dataSet1,dataSet2)
        initialClusterCenters = 'plusplus';
        trainingPlotVisualization = false;
    end

    properties (SetAccess = protected)
        clusterCenters = [];   % The cluster centers
    end
    properties (SetAccess = private, Hidden = true)
        uY = [];
    end

    methods

        function self = set.nClusters(self,value)
            if ~prtUtilIsPositiveScalarInteger(value)
                error('prt:prtClusterKmodes:nClusters','value (%s) must be a positive scalar integer',mat2str(value));
            end
            self.nClusters = value;
        end

        function self = set.handleEmptyClusters(self,value)
            if ~isa(value,'char') || ~any(strcmpi(value,{'exit','regularize'}))
                error('prt:prtClusterKmodes:handleEmptyClusters','value (%s) must be one of ''remove'', or ''random''',mat2str(value));
            end
            self.kmeansHandleEmptyClusters = value;
        end
        function self = set.initialClusterCenters(self,value)
            if isa(value,'char')
                if ~any(strcmpi(value,{'random','plusplus'}))
                    error('prt:prtClusterKmodes:initialClusterCenters','If specified as a string initialClusterCenters must be one of ''random'', or ''plusplus''');
                end
            else
                % initialMeans were specified as a matrix. Assume
                % everything will be ok and let prtUtilKMeans catch errors
            end
            self.initialClusterCenters = value;
        end

        function self = set.distanceMetricFn(self,value)
            if ~isa(value,'function_handle')
                error('prt:prtClusterKmeans:distanceMetricFn','distanceMetricFn must be a function handle');
            elseif nargin(value) ~= 2
                error('prt:prtClusterKmeans:distanceMetricFn','distanceMetricFn must be a function handle that takes two input arguments');
            end
            self.distanceMetricFn = value;
        end
        function self = prtClusterKmodes(varargin)
            % Allow for string, value pairs
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,ds)
            %self.clusterCenters = prtUtilKmeans(ds.getObservations,slf.nClusters,'distanceMetricFn',Obj.distanceMetricFn,'handleEmptyClusters',Obj.kmeansHandleEmptyClusters,'initialMeans',Obj.initialMeans);

            if ischar(self.initialClusterCenters)
                nInitializeKMeansIterations = 1;
                self.clusterCenters = prtUtilKmeans(ds.getObservations,self.nClusters,'distanceMetricFn',self.distanceMetricFn,'handleEmptyClusters',self.handleEmptyClusters,'initialMeans',self.initialClusterCenters,'maxIterations',nInitializeKMeansIterations);
            else
                self.clusterCenters = self.initialClusterCenters;e
            end

            data = ds.data;

            [nSamples,nDimensions] = size(data);
            clusterIndexOld = nan(nSamples,1);

            nMaxIterations = 100;
            for iteration = 1:nMaxIterations

                distanceMat = self.distanceMetricFn(data,self.clusterCenters);
                [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>

                if self.trainingPlotVisualization
                    self.trainingIterationPlotVisualization(data,self.clusterCenters,clusterIndex,iteration);
                end
                %Handle empty clusters:
                nMaxFixSteps = 10;
                for iFix = 1:nMaxFixSteps
                    if length(unique(clusterIndex)) ~= self.nClusters
                        invalidClusters = setdiff(1:self.nClusters,clusterIndex);
                        validClusters = intersect(1:self.nClusters,clusterIndex);
                        switch lower(inputStructure.handleEmptyClusters)
                            case 'remove'
                                classMeans = classMeans(validClusters,:);
                                self.nClusters = self.nClusters - length(invalidClusters);
                                distanceMat = distanceMat(:,validClusters);
                                [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>
                            case 'random'
                                randInds = max(1,ceil(rand(1,length(invalidClusters))*nSamples));
                                classMeans(invalidClusters,:) = data(randInds,:);
                                distanceMat =  inputStructure.distanceMetricFn(data,classMeans);
                                [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>
                            otherwise
                                error('invalid');
                        end
                    else
                        break
                    end
                end
                nSamplesForCdf = 1000;
                newClusters = zeros(self.nClusters,nDimensions);
                for clusterInd = 1:self.nClusters

                    cData = data(clusterIndex == clusterInd,:);
                    for iDim = 1:nDimensions
                        % newClusters(clusterInd,iDim) = % mode(cData(:,iDim)); % would work for discrete data only

                        cRv = mle(self.kdeRv,cData(:,iDim));

                        cXSamples = sort(cat(1,linspace(min(cData(:,iDim)), max(cData(:,iDim)),nSamplesForCdf)',cData(:,iDim)),'ascend');
                        cCdf = cRv.cdf(cXSamples);

                        %cModeInd = find(cCdf >= 0.5,1,'first');
                        [~,cModeInd] = max(cCdf);
                        if isempty(cModeInd)
                            % Something went horribly wrong
                            cMode = mean(cData(:,iDim));
                        else
                            cMode = cXSamples(cModeInd);
                        end

                        newClusters(clusterInd,iDim) = cMode;
                    end
                end
                self.clusterCenters = newClusters;

                if iFix == nMaxFixSteps
                    %maxIterReached = true;
                    return
                end

                %Check convergence:
                if all(clusterIndexOld == clusterIndex)
                    %maxIterReached = false;
                    return;
                else
                    clusterIndexOld = clusterIndex;
                end
            end

            % maxIterReached = iteration == nMaxIterations;

        end

        function ds = runAction(self,ds)

            fn = self.distanceMetricFn;
            distance = fn(ds.getObservations,self.clusterCenters);

            if size(distance,1) ~= ds.nObservations || size(distance,2) ~= size(self.clusterCenters,1)
                error('prt:prtClusterKmodes:badDistanceMetric','Expected a matrix of size %s from the distance metric, but distance metric function output a matrix of size %s',mat2str([DataSet.nObservations, size(Obj.clusterCenters,1)]),mat2str(size(distance)));
            end

            %The smallest distance is the expected class:
            [dontNeed,clusters] = min(distance,[],2);  %#ok<ASGLU>

            binaryMatrix = zeros(size(clusters,1),self.nClusters);
            for i = 1:self.nClusters
                binaryMatrix(i == clusters,i) = 1;
            end
            ds = ds.setObservations(binaryMatrix);
        end
        end
    methods (Hidden)
        function trainingIterationPlotVisualization(~, data,classMeans,clusterIndex,iter)

            [nSamples,nDimensions] = size(data); %#ok<ASGLU>

            if nDimensions > 3
                warning('prt:prtClusterKmods:iterationPlot','plotVisualization is true, but dimensionality of data (%d) is > 3',nDimensions);
                return;
            end
            ds = prtDataSetClass(data,clusterIndex);

            plot(ds);
            hold on;
            switch nDimensions
                case 1
                    plot(classMeans,'b.');
                case 2
                    plot(classMeans(:,1),classMeans(:,2),'b.');
                case 3
                    plot3(classMeans(:,1),classMeans(:,2),classMeans(:,3),'b.');
            end
            hold off;
            title(sprintf('Iteration %d',iter));
            drawnow; %pause;

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef RMMEDA_F2 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing RM-MEDA

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, and Y. Jin, RM-MEDA: A regularity model-based
% multiobjective estimation of distribution algorithm, IEEE Transactions on
% Evolutionary Computation, 2008, 12(1): 41-63.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            g = 1 + 9*mean((X(:,2:end)-repmat(X(:,1),1,size(X,2)-1)).^2,2);
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-(PopObj(:,1)./g).^2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% imbibition capillary pressure curve
% Written by Ali A. Eftekhari
function res=dpc_imb(sw, pce_w, pce_o, swc, sor, labda_w, labda_o, pc_max_w, pc_max_o)
  dpc1=dpc_drain(sw, pce_w, swc, labda_w, pc_max_w);
  dpc2=dpc_drain(1-sw, pce_o, sor, labda_o, pc_max_o);
  res=dpc1-dpc2;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function tf = isTestString(str)
%isTestString True if string looks like the name of a test
%   tf = isTestString(str) returns true if the string str looks like the name of
%   a test.  If str is a cell array of strings, then isTestString tests each
%   string in the cell array, returning the results in a logical array with the
%   same size as str.

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

test_exp = '^[tT]est';
tf = mtest.utils.containsRegexp(str, test_exp);

author: ricedsp


text:
function varargout = Data_Analyzer(varargin)
% DATA_ANALYZER MATLAB code for Data_Analyzer.fig
%      DATA_ANALYZER, by itself, creates a new DATA_ANALYZER or raises the existing
%      singleton*.
%
%      H = DATA_ANALYZER returns the handle to a new DATA_ANALYZER or the handle to
%      the existing singleton*.
%
%      DATA_ANALYZER('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in DATA_ANALYZER.M with the given input arguments.
%
%      DATA_ANALYZER('Property','Value',...) creates a new DATA_ANALYZER or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before Data_Analyzer_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to Data_Analyzer_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 Data_Analyzer

% Last Modified by GUIDE v2.5 25-Apr-2014 11:25:51

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Data_Analyzer_OpeningFcn, ...
                   'gui_OutputFcn',  @Data_Analyzer_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 Data_Analyzer is made visible.
function Data_Analyzer_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 Data_Analyzer (see VARARGIN)

handles.fSelect = 1;
movegui('center'); % Moves GUI to center of screen

% Choose default command line output for Data_Analyzer
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = Data_Analyzer_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 button press in loadButton.
function loadButton_Callback(hObject, eventdata, handles)
% hObject    handle to loadButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

evalin('base','clear all');
uiopen('load');

if exist('Seconds')
assignin('base', 'Seconds', Seconds);
assignin('base', 'Throttle', Throttle);
assignin('base', 'RPM', RPM);

if exist('gramsMeas')
assignin('base', 'gramsMeas', gramsMeas);
end

end

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

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

handles.fSelect = get(hObject,'Value');

guidata(hObject,handles);

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

% Hint: popupmenu 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

function bedit_Callback(hObject, eventdata, handles)
% hObject    handle to bedit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of bedit as text
%        str2double(get(hObject,'String')) returns contents of bedit as a double

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

% Hint: edit 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

function hedit_Callback(hObject, eventdata, handles)
% hObject    handle to hedit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of hedit as text
%        str2double(get(hObject,'String')) returns contents of hedit as a double

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

% Hint: edit 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

function lowCutoff_Callback(hObject, eventdata, handles)
% hObject    handle to lowCutoff (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of lowCutoff as text
%        str2double(get(hObject,'String')) returns contents of lowCutoff as a double

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

% Hint: edit 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

function highCutoff_Callback(hObject, eventdata, handles)
% hObject    handle to highCutoff (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of highCutoff as text
%        str2double(get(hObject,'String')) returns contents of highCutoff as a double

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

% Hint: edit 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

function inputValues_Callback(hObject, eventdata, handles)
% hObject    handle to inputValues (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of inputValues as text
%        str2double(get(hObject,'String')) returns contents of inputValues as a double

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

% Hint: edit 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 runButton.
function runButton_Callback(hObject, eventdata, handles)
% hObject    handle to runButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

Seconds = evalin('base', 'Seconds');
Throttle = evalin('base', 'Throttle');
RPM = evalin('base', 'RPM');
existCheck = 'exist(''gramsMeas'',''var'')'; % String to evaluate from context of base workspace
if evalin('base',existCheck)
gramsMeas = evalin('base', 'gramsMeas');
end

% Check to see if the user inputs low/high cutoffs:
inputCheck = get(handles.lowCutoff, 'String');

% If the user doesn't input anything:
if isempty(inputCheck)

    LowCutoff = min(Throttle);
    HighCutoff = max(Throttle);

% Or, if they manually input the low/high cutoffs:
else
    LowCutoff = str2num(get(handles.lowCutoff, 'String'));
    HighCutoff = str2num(get(handles.highCutoff, 'String'));
end

b = str2num(get(handles.bedit, 'String')); % assumed inches
h = str2num(get(handles.hedit, 'String')); % assumed inches
L = str2num(get(handles.Ledit, 'String')); % assumed inches, converted within function
switch handles.fSelect % Get Tag of selected object
    case 1

        calculate_CR_B(Seconds, Throttle, RPM, LowCutoff, HighCutoff);

    case 2

        calculate_CR_B(Seconds, Throttle, RPM, LowCutoff, HighCutoff);
        calculate_CT(Seconds, Throttle, RPM, b, h, gramsMeas, LowCutoff, HighCutoff);

    case 3

        calculate_CR_B(Seconds, Throttle, RPM, LowCutoff, HighCutoff);
        calculate_CQ(Seconds, Throttle, RPM, L, gramsMeas, LowCutoff, HighCutoff);

    case 4

        calculate_TC(Seconds, Throttle, RPM);
end

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

% Hint: edit 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

function Ledit_Callback(hObject, eventdata, handles)
% hObject    handle to Ledit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Ledit as text
%        str2double(get(hObject,'String')) returns contents of Ledit as a double

repo: Quad-Sim

path: github-repos/MATLAB/dch33-Quad-Sim/Quad-Sim-961bc69d4939c8d0661eeb9820de668994262f65/Quadcopter Dynamic Modeling and Simulation/Data Acquisition & Analysis/MATLAB Data Analysis/Data_Analyzer.m


text:
function out = reorderLegend(order,h)
% reorderLegend
%   Re-orders the legend entries
%   (without modifying the actual order the curves are stacked on the axes)
%
%   order  - an array specifying the order of legend entries
%   h      - axis handle (optional, defaults to current axis)
%
%   Invoking the function without an argument or with order==[]
%   causes the labels of the curves to be returned (or displayed, if
%   nargout==0).
%       e.g.
%           >> x = 0:.1:10;
%           >> plot(x,x,x,2*x,x,3*x,x,4*x); legend('x','2x','3x','4x');
%           >> reorderLegend
%           [1] 4x
%           [2] 3x
%           [3] 2x
%           [4] x
%
%   The indices of the labels (i.e., 1 to 4 in the above e.g.)
%   can be used to specify the desired order of legend entries.
%       e.g.
%           >> reorderLegend([4,2,1,3]);
%               % legend: 'x','3x','4x','2x' (from top to bottom)
%
%   Notes:
%   - By default, the legend entries are in the reverse order of the
%   indices (because the curves are in a stack).
%   - `reorderLegend' doesn't change the order of the curves
%   (i.e., the indices are not changed when it's called).
%
%   Author        : Damith Senaratne, (http://www.damiths.info)
%   Released date : 21th June 2012

    if nargin < 2
        h = gca;
    end

    hC = findobj(h,'Type','line');

    if nargin == 0 || isempty(order)
        % output/display the internal order
        lbl = get(hC,{'DisplayName'});
        if nargout == 1
            out = lbl;
        else
            for k = 1:length(lbl)
               display([sprintf('[%*d] ',fix(log10(length(lbl)))+1,k) lbl{k}]);
            end
        end

    else
       % reorder the legend entries
       legend(hC(order));
    end
end


text:
 function g = gibbsperiodicKernDiagGradient(kern, x, covDiag)

% GIBBSPERIODICKERNDIAGGRADIENT Compute the gradient of the GIBBSPERIODIC kernel's diagonal wrt parameters.
% FORMAT
% DESC computes the gradient of functions of the diagonal of the
% Gibbs-kernel derived periodic kernel matrix with respect to the parameters of the kernel. The
% parameters' gradients are returned in the order given by the
% gibbsperiodicKernExtractParam command.
% ARG kern : the kernel structure for which the gradients are
% computed.
% ARG x : the input data for which the gradient is being computed.
% ARG factors : partial derivatives of the function of interest with
% respect to the diagonal elements of the kernel.
% RETURN g : gradients of the relevant function with respect to each
% of the parameters. Ordering should match the ordering given in
% gibbsperiodicKernExtractParam.
%
% SEEALSO : gibbsperiodicKernParamInit, kernDiagGradient, gibbsperiodicKernExtractParam, gibbsperiodicKernGradient
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

g = zeros(1, kern.nParams);
g(end) = sum(covDiag);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [maps, mapeq, maps_b, mapeq_b] = genSDPConstraintMap(n,k, W)

    maps =  {}; %cell(2*n,1);
    mapeq = {};
    eqCount = 1; ieqCount = 1;
    for i=1:n
        A = zeros(n,n);
        A(i,:) = 1;
        Av = [A(:); 0];

        M = sparse(n*n+1, n*n+1);
        M(n*n+1, :) = Av;
        maps{ieqCount} = M; ieqCount = ieqCount + 1;
    end

    for i=1:n
        A = zeros(n,n);
        A(:,i) = 1;
        Av = [A(:); 0];

        M = sparse(n*n+1, n*n+1);
        M(n*n+1, :) = Av;
        maps{ieqCount} = M; ieqCount = ieqCount + 1;
    end
    maps_b = ones (2*n, 1);

    %----------------------------------------------------

    % 1^TX1 = k
    M = sparse(n*n+1, n*n+1);
    M(n*n+1,1:n*n) = 1;
    mapeq{eqCount} = M; eqCount = eqCount+1;
    mapeq_b = k;

    % Sum Yqrst = k
%     M = ones(n*n, n*n);
%     M = [M; zeros(1,n*n)];
%     M = [M zeros(n*n+1,1)];
%     mapeq{eqCount} = M; eqCount = eqCount + 1;
%     mapeq_b = [mapeq_b; k*k];

    % Wherever W(i,j) < 0, force Y(i,j) to 0
%     for i=1:size(W,1)
%         for j=1:size(W,2)
%            if (W(i,j) < 0)
%                M = sparse(n*n+1, n*n+1); M(i,j) = 1;
%                mapeq{eqCount} = M; eqCount = eqCount + 1;
%                mapeq_b = [mapeq_b; 0];
%             end
%         end
%     end
% %
%
    % Force Yqrst = 0;
%      for p=1:n
%        for q=1:n
%            for s=1:n
%                for t=1:n
%                    i = p + (q-1)*n;
%                    j = s + (t-1)*n;
%                    if (p==s && q~=t)||(p~=s && q==t)
%                        M = sparse(n*n+1, n*n+1); M(i,j) = 1;
%                        mapeq{eqCount} = M; eqCount = eqCount + 1;
%                        mapeq_b = [mapeq_b; 0];
%                    else
% %                        M = sparse(n*n+1, n*n+1);
% %                        M(i,j) = 1; M(i,n*n+1) = -1;
% %                        maps_b = [maps_b;0];
% %                        maps{ieqCount} = M; ieqCount = ieqCount + 1;
% %
% %                        if (i~=j)
% %                            M = sparse(n*n+1, n*n+1);
% %                            M(i,j) = 1; M(j,n*n+1) = -1;
% %                            maps_b = [maps_b;0];
% %                            maps{ieqCount} = M; ieqCount = ieqCount + 1;
% %                        end
%                    end
%
%
%
%                 end
%             end
%         end
%     end

end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function M = setdiag(M, v)
% SETDIAG Set the diagonal of a matrix to a specified scalar/vector.
% M = set_diag(M, v)

n = length(M);
if length(v)==1
  v = repmat(v, 1, n);
end
J = 1:n+1:n^2;
M(J) = v;

%M = triu(M,1) + tril(M,-1) + diag(v);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [S,C]=allsums(varargin)
% ALLSUMS Distribution of unique sums among all combinations of vectors.
%   [S,C]=ALLSUMS(IN) depending on input IN, will do one of two different
%   procedures for calculating all unique combinations of sums.
%
%   The output vectors, S and C, are respectively the list of unique
%   sums, and the count distribution for those sums among all combinations.
%
%
%   SINGLE INPUT:
%
%   [S,C]=ALLSUMS(VEC) for a single input vector of length N, returns the
%   distribution of unique sums among the (2^N)-1 combinations of
%   non-empty subsets of elements within VEC.
%
%    EXAMPLE:
%      Number of ways to have a total of 250 using only the odd numbers
%      from 1 to 100.
%      [S,C]=allsums(1:2:100);
%      num=C(S==250)
%
%
%   MULTIPLE INPUTS:
%
%   [S,C]=ALLSUMS(VEC,K) with vector VEC of length N, assumes the
%   equivalent scenario of rolling a N-sided die, with values VEC,
%   K times and returns the distributions of unique sums among all
%   N^K combinations of values.
%
%   [S,C]=ALLSUMS(VEC1,K1,VEC2,K2,...) with vectors of length N1, N2, ...
%   assumes the equivalent scenario of rolling a N1-sided die, with
%   values VEC1, K1 times; then rolling a N2-sided dice, with values VEC2,
%   K2 times; etc. and returns the distribution of unique sums among
%   all (N1^K1)*(N2^K2)*... combinations of values
%
%    SINGLE VECTOR, REPEATED
%      Probability distribution of rolling a regular 6-sided die 100 times
%      [S,C]=allsums(1:6,100);
%      plot(S,C./sum(C)), grid on
%
%    MULTIPLE VECTORS, REPEATED
%      Probability of the total yielding exactly zero when rolling a
%      A) six-sided die with values [-5:0] ten times, and then a
%      B) ten-sided die with values [-4:5] twenty times.
%      [S,C]=allsums([-5:0],10,[-4:5],20);
%      prob=C(S==0)/sum(C)
%

%   Mike Sheppard
%   Last Modified: 27-Oct-2011

%ERROR CHECKING
ninputs=numel(varargin);
switch ninputs
    case 0
        error('allsums:Inputs','Requires at least one input');
    case 1
        VEC=varargin{1}; %Input vector
        if any(~isfinite(VEC))
            error('allsums:Inputs','Requires input vector to be numeric');
        end
    otherwise
        num=ninputs/2;
        if mod(ninputs,2)~=0
            error('allsums:Inputs','Requires inputs to be in [Vector],[Number] order');
        end
        for k=1:num
            ALLVEC{k}=varargin{2*k-1}; ALLNUMR{k}=varargin{2*k}; numr=ALLNUMR{k};
            if any(~isfinite(ALLVEC{k}))
                error('allsums:Inputs','Requires all input vectors to be numeric');
            end
            if any([~isfinite(numr) ~isscalar(numr) numr~=round(numr) numr<=0])
                error('allsums:Inputs','Requires number of repetitions to be a positive scalar integer');
            end
        end
end

%ALGORITHM
switch ninputs
    case 1
       % Distribution of unique sums among the (2^N)-1 combinations of
       % non-empty subsets of elements within VEC.
        S=[]; C=[];
        for k=1:numel(VEC)
            [S,ign,indx]=unique([S VEC(k) S+VEC(k)]);
            C=accumarray(indx',[C 1 C])';
        end
        S=S(:); C=C(:);
    otherwise
        % Distribution of unique sums among all (N1^K1)*(N2^K2)*...
        % combinations of vectors of length N_k repeated K_k times,
        % with values VEC_k each.
        S=0; C=1;
        for k=1:num
            VEC=ALLVEC{k}; VEC=VEC(:); numr=ALLNUMR{k};
            for n=1:numr
                tempv=bsxfun(@plus,S',VEC);
                tempc=repmat(C',length(VEC),1);
                [S,ign,indx]=unique(tempv(:));
                C=accumarray(indx,tempc(:));
            end
        end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function [mergedImage, newDimLabel] = ...
    merge(this, mergeDims, varargin)
% Merges multiple dimensions into one dimensions
%
%   Y = MrImage()
%   mergedImage = Y.merge({'echo', 't'})
%
% This is a method of class MrImage.
%
% IN
%   mergeDims       which dims should be merged into one dim
%                   can be numeric [1,4] or character {'coil', 't'}
%                   Note that the order of mergeDims in the input
%                   determines their order in the mergedImage, i.e.
%                   I.merge({'coil', 't'}) is different from I.merge({'t',
%                   'coil'}). Also, the mergedDim is always the last in the
%                   new image; use permute to change that.
%
%   varargin        prop/val pairs to describe the new dimension
%                   including resolutions, ranges, dimLabels, units,
%                   samplingPoints, samplingWidths
%                   defaults:   samplingPoints = 1:nSamples
%                               dimLabels = [oldDimLabels(1), '_', ...
%                               '_', olDimLabels(end)]
%                               units = ''
% OUT
%
% EXAMPLE
%   merge
%
%   See also MrImage

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-12-23
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

mergedImage = this.copyobj();
newDimLabel = [];
if ~isempty(mergeDims)
    % create merged dimInfo first - this also give the new labels and sampling
    % points
    [mergedDimInfo, ~, newDimLabel] = this.dimInfo.merge(mergeDims, varargin);

    % only split if more than one mergeDim is actually given
    if numel(mergeDims) > 1
        % split into individual images
        split_array = this.split('splitDims', mergeDims);
        split_array = reshape(split_array, 1, []);

        % overwrite the sampling points to the newly created ones
        for n = 1:numel(split_array)
            split_array{n}.dimInfo = mergedDimInfo.select(newDimLabel, n);
        end

        % use combine to merge image along non-singleton dimension
        mergedImage = split_array{1}.combine(split_array);
    else % otherwise just change the dimLabel
        mergedImage.dimInfo = mergedDimInfo;
    end
end
end

author: translationalneuromodeling

repo: tapas


text:
function [RMAOV1MS] = RMAOV1MS(X,alpha)
%RMAOV1MS Repeated measures single-factor analysis of variance multiple samples
% test.
% This m-file deals with repeated measurements at t time points obtained from
% s groups of subjects or multiple samples. With or without same size. The
% simplest model is:
%
%           x_hij = m + g_h +t_j + (gt)_hj + p_i(h) + e_hij
%
% where: x_hij = response at time j from the ith subject in group h; m =
%        overall mean; g_h = fixed effect of group h; t_j = fixed effect of
%        time j; (gt)_hj = fixed effect for the interaction of the hth group
%        with the jth time; p_i(h) = random effects for the ith subject in
%        the hth group, and e_hij = independent random error terms.
%
% For we are interested test for differences among groups, test differences
% among time points, and test the significance of the group x time interaction.
% The first test requires the assumption that the within-group covariance
% matrices are equal. The two last-ones tests require this same assumption
% and that the sphericity condition is satisfied. In general, these assumtions
% are required for all tests of within-subjects effects (repeated measures).
%
% Syntax: RMAOV1MS(X,alpha)
%
% Inputs:
%     X - data matrix (Size of matrix must be n-by-4;dependent variable=
%         column 1;group=column 2;independent variable=column 3;subject=
%         column 4).
% alpha - significance level (default = 0.05).
%
% Outputs:
%       - Complete Analysis of Variance Table.
%
% Example: From the example 5.4.2 of Charles S. Davis (2002). from an experiment
%          comparing dissolution times of pills. We use the logarithms of the
%          times, in seconds, for each pill to reach fractions remaining of
%          0.9,0.7,0.5,0.3,0.25, and 0.10. The four groups correspond to
%          different storage conditions. Under the assumption thet the
%          logarithms of the times are normally distributed, repeated measures
%          ANOVA can be used to test for effects due to group, fraction
%          remaining, and interaction effect between group and fraction
%          remaining. Tha data are shown on the below table.
%
%          -----------------------------------------------------------
%                                        Fraction remainig
%                           ------------------------------------------
%          Group    Tablet    0.90   0.70   0.50   0.30   0.25   0.10
%          -----------------------------------------------------------
%            1         1       2.56  2.77   2.94   3.14   3.18   3.33
%                      2       2.64  2.89   3.09   3.26   3.33   3.47
%                      3       2.94  3.18   3.33   3.50   3.50   3.66
%                      4       2.56  2.83   3.04   3.22   3.26   3.37
%                      5       2.64  2.77   2.94   3.14   3.22   3.30
%                      6       2.56  2.77   2.94   3.14   3.18   3.26
%            2         1       2.56  2.83   3.07   3.26   3.33   3.51
%                      2       2.44  2.74   3.02   3.20   3.28   3.44
%                      3       2.34  2.67   2.91   3.16   3.22   3.39
%                      4       2.41  2.71   2.94   3.16   3.21   3.36
%            3         1       2.46  2.83   3.09   3.32   3.37   3.54
%                      2       2.60  2.93   3.21   3.40   3.46   3.62
%                      3       2.48  2.84   3.12   3.35   3.41   3.58
%                      4       2.49  2.82   3.05   3.29   3.37   3.52
%            4         1       2.40  2.67   2.94   3.20   3.26   3.47
%                      2       2.64  2.94   3.18   3.40   3.45   3.66
%                      3       2.40  2.64   2.86   3.09   3.16   3.38
%          -----------------------------------------------------------
%
%     Data matrix must be:
% X=[2.56 1 1 1;2.77 1 1 2;2.94 1 1 3;3.14 1 1 4;3.18 1 1 5;3.33 1 1 6;2.64 1 2 1;2.89 1 2 2;3.09 1 2 3;3.26 1 2 4;3.33 1 2 5;3.47 1 2 6;
% 2.94 1 3 1;3.18 1 3 2;3.33 1 3 3;3.50 1 3 4;3.50 1 3 5;3.66 1 3 6;2.56 1 4 1;2.83 1 4 2;3.04 1 4 3;3.22 1 4 4;3.26 1 4 5;3.37 1 4 6;
% 2.64 1 5 1;2.77 1 5 2;2.94 1 5 3;3.14 1 5 4;3.22 1 5 5;3.30 1 5 6;2.56 1 6 1;2.77 1 6 2;2.94 1 6 3;3.14 1 6 4;3.18 1 6 5;3.26 1 6 6;
% 2.56 2 1 1;2.83 2 1 2;3.07 2 1 3;3.26 2 1 4;3.33 2 1 5;3.51 2 1 6;2.44 2 2 1;2.74 2 2 2;3.02 2 2 3;3.20 2 2 4;3.28 2 2 5;3.44 2 2 6;
% 2.34 2 3 1;2.67 2 3 2;2.91 2 3 3;3.16 2 3 4;3.22 2 3 5;3.39 2 3 6;2.41 2 4 1;2.71 2 4 2;2.94 2 4 3;3.16 2 4 4;3.21 2 4 5;3.36 2 4 6;
% 2.46 3 1 1;2.83 3 1 2;3.09 3 1 3;3.32 3 1 4;3.37 3 1 5;3.54 3 1 6;2.60 3 2 1;2.93 3 2 2;3.21 3 2 3;3.40 3 2 4;3.46 3 2 5;3.62 3 2 6;
% 2.48 3 3 1;2.84 3 3 2;3.12 3 3 3;3.35 3 3 4;3.41 3 3 5;3.58 3 3 6;2.49 3 4 1;2.82 3 4 2;3.05 3 4 3;3.29 3 4 4;3.37 3 4 5;3.52 3 4 6;
% 2.40 4 1 1;2.67 4 1 2;2.94 4 1 3;3.20 4 1 4;3.26 4 1 5;3.47 4 1 6;2.64 4 2 1;2.94 4 2 2;3.18 4 2 3;3.40 4 2 4;3.45 4 2 5;3.66 4 2 6;
% 2.40 4 3 1;2.64 4 3 2;2.86 4 3 3;3.09 4 3 4;3.16 4 3 5;3.38 4 3 6];
%
% Calling on Matlab the function:
%             RMAOV1MS(X)
%
% Answer is:
%
% The number of groups are: 4
% The maximum number of subjects in a group are: 6
% The minimum number of subjects in a group are: 3
% The number of IV levels are: 6
%
% ANOVA table for repeated measures single-factor multiple samples.
% ------------------------------------------------------------------------------
% SOV                        SS          df          MS            F        P
% ------------------------------------------------------------------------------
% Groups                   0.2038         3        0.0679         0.87   0.4798
% Subjts(Gps.)             1.0108        13        0.0778
% IV                      10.0745         5        2.0149      3198.78   0.0000
% Group x IV               0.2045        15        0.0136        21.64   0.0000
% Residual                 0.0409        65        0.0006
% ------------------------------------------------------------------------------
% Total                   11.5344       101
% ------------------------------------------------------------------------------
% If the P results are smaller than 0.05
% the corresponding Ho's tested result statistically significant. Otherwise,
% are not significative.
%
% Created by A. Trujillo-Ortiz, R. Hernandez-Walls, K. Barba-Rojo and
%            A. Castro-Perez
%            Facultad de Ciencias Marinas
%            Universidad Autonoma de Baja California
%            Apdo. Postal 453
%            Ensenada, Baja California
%            Mexico.
%            atrujo@uabc.mx
%
% Copyright.December 13, 2006.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls, K. Barba-Rojo and A. Castro-Perez. (2006).
%   RMAOV1MS:Repeated measures single-factor analysis of variance multiple samples test.
%   A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/
%   fileexchange/loadFile.do?objectId=13380
%
% Reference:
%    Davis, S. Ch. (2002), Statistical Methods for the Analysis of Repeated
%             Measurements. Springer-Verlag:New York. Chapter 5.
%

if nargin < 2,
   alpha = 0.05; %(default)
end;

if (alpha <= 0 | alpha >= 1)
   fprintf('Warning: significance level must be between 0 and 1\n');
   return;
end;

if nargin < 1,
   error('Requires at least one input argument.');
   return;
end;

a = max(X(:,2)); %Number of groups.
b = max(X(:,3)); %Maximum number of subjects in a group.
t = max(X(:,4)); %Number of levels of independent variable.

indice = X(:,2);
for i = 1:a
   Xe = find(indice==i);
   eval(['G' num2str(i) ' = X(Xe,1);']);
end
indice = X(:,4);
for k = 1:t
    Xe = find(indice==k);
    eval(['IV' num2str(k) '=X(Xe,1);']);
end

C = (sum(X(:,1)))^2/length(X(:,1));  %Correction term.
SSTO = sum(X(:,1).^2)-C;  %Total sum of squares.
dfTO = length(X(:,1))-1;  %Total degrees of freedom.

for i = 1:a
   for j = 1:b
      Xe = find((X(:,2)==i) & (X(:,3)==j));
      eval(['T' num2str(i) num2str(j) ' = X(Xe,1);']);
   end
end

for i = 1:a
    for k = 1:t
        Xe = find((X(:,2)==i) & (X(:,4)==k));
        eval(['GI' num2str(i) num2str(k) '=X(Xe,1);']);
     end
end

%Procedure related to the subjects in groups.
NSG = [];
for i = 1:a
   for j = 1:b
   eval(['n' num2str(i) num2str(j) ' = length(T' num2str(i) num2str(j) ');'])
   eval(['x =  n' num2str(i) num2str(j) ';'])
   NSG = [NSG,x];
   end
end

%Due it is an unequal sample sizes design, some subgroups arithmetic operations
%could been not-a-number (like 0.0/0.0). So it is necessary to make the estimations
%ignoring NaNs.
warning off;
SG = [];
for i = 1:a
   for j = 1:b
   eval(['x = ((sum(T' num2str(i) num2str(j) ').^2)/length(T' num2str(i) num2str(j) '));'])
   SG = [SG,x];
   nans = isnan(SG);
   r = find(nans);
   SG(r) = zeros(size(r));
   end
end
warning on;
SSSG = sum(SG)-C;
dfSG = (a*b)-1;

%Procedure related to the groups.
G = [];
for i = 1:a
   eval(['x = ((sum(G' num2str(i) ').^2)/length(G' num2str(i) '));'])
   G = [G,x];
end
SSG = sum(G)-C;
dfG = a-1;
MSG = SSG/dfG;

%Procedure related to the IV (independent variable).
IV = [];
for k = 1:t
    eval(['x =((sum(IV' num2str(k) ').^2)/length(IV' num2str(k) '));']);
    IV = [IV,x];
 end
SSIV = sum(IV)-C;  %sum of squares for the IV
dfIV = t-1;  %degrees of freedom for the IV
MSIV = SSIV/dfIV;  %mean square for the IV

%Procedure related to the interaction group-independent variable.
GI = [];
for i = 1:a
    for k = 1:t
        eval(['x =((sum(GI' num2str(i) num2str(k) ').^2)/length(GI' num2str(i) num2str(k) '));']);
        GI = [GI,x];
     end
 end
SSGI = sum(GI)-sum(G)-sum(IV)+C;  %sum of squares of the Group x IV
dfGI = dfG*dfIV;  %degrees of freedom of the Group x IV
MSGI = SSGI/dfGI;  %mean square for the Group x IV

%Procedure for estimation of sample size for unequal sample sizes.
no1 = []; no2 = [];
for i = 1:a
    for j = 1:b
        eval(['x = length(T' num2str(i) num2str(j) ');']);
        no1 = [no1,x];
        eval(['x = length(T' num2str(i) num2str(j) ').^2;']);
        no2 = [no2,x];
    end
end

%Procedure related to the within subgroups (error).
SSE = SSTO-SSSG;
nneg = length(find((no1-1)<0));
dfE = sum(no1-1)+nneg;
MSE = SSE/dfE;

%Procedure related to the among subgroups within groups.
SSSGG = SSTO-SSG-SSE;
dfSGG = dfTO-dfG-dfE;
MSSGG = SSSGG/dfSGG;

%Procedure related to the residual.
SSR = SSTO-SSG-SSSGG-SSIV-SSGI;
dfR = dfTO-dfG-dfSGG-dfIV-dfGI;
MSR = SSR/dfR;

fprintf('The number of groups are:%2i\n', a);
fprintf('The maximum number of subjects in a group are:%2i\n', b);
c = min(full(sum(spones(reshape(no2,b,a))))); %Minimum number of subjects in a group.
fprintf('The minimum number of subjects in a group are:%2i\n', c);
fprintf('The number of IV levels are:%2i\n\n', t);
FG = MSG/MSSGG;
FIV = MSIV/MSR;
FGI = MSGI/MSR;
PrG = 1 - fcdf(FG,dfG,dfSGG);  %Probability associated to the groups F statistic.
PrIV = 1 - fcdf(FIV,dfIV,dfR);  %Probability associated to the IV F statistic.
PrGI = 1 - fcdf(FGI,dfGI,dfR); %Probability associated to the interaction group-IV F statistic.
disp('ANOVA table for repeated measures single-factor multiple samples.')
fprintf('------------------------------------------------------------------------------\n');
disp('SOV                        SS          df          MS            F        P       ')
fprintf('------------------------------------------------------------------------------\n');
fprintf('Groups       %18.4f%10i%14.4f%13.2f%9.4f\n',SSG,dfG,MSG,FG,PrG);
fprintf('Subjts(Gps.)%19.4f%10i%14.4f\n',SSSGG,dfSGG,MSSGG);
fprintf('IV    %25.4f%10i%14.4f%13.2f%9.4f\n',SSIV,dfIV,MSIV,FIV,PrIV);
fprintf('Group x IV %20.4f%10i%14.4f%13.2f%9.4f\n',SSGI,dfGI,MSGI,FGI,PrGI);
fprintf('Residual %22.4f%10i%14.4f\n',SSR,dfR,MSR);
fprintf('------------------------------------------------------------------------------\n');
fprintf('Total %25.4f%10i\n',SSTO,dfTO);
fprintf('------------------------------------------------------------------------------\n');
fprintf('If the P results are smaller than% 3.2f\n', alpha );
disp('the corresponding Ho''s tested result statistically significant. Otherwise,');
disp('are not significative.');
disp('     ')

return,


text:
classdef ShFunWithElasticPdes < ShapeFunctional

    properties (Access = protected)
        interpolation
        physicalProblem
    end

    properties (Access = private)
        orientationUpdater
        alpha
    end

    properties (Access = public)
        regDesignVariable
    end

    methods (Access = public)

        function computeFunctionAndGradient(obj)
            obj.updateAlpha();
            obj.computeFunction();
            obj.updateAlpha();
            obj.computeGradient();
            obj.updateAlpha();
        end

        function computeFunction(obj)
            obj.updateHomogenizedMaterialProperties();
            obj.solveState();
            obj.computeFunctionValue();
            obj.normalizeFunction();
        end

        function computeGradient(obj)
            obj.solveAdjoint();
            obj.computeGradientValue();
            obj.filterGradient();
            obj.normalizeGradient();
        end

        function f = getPhysicalProblems(obj)
            f{1} = obj.physicalProblem;
        end

        function f = getRegularizedDesignVariable(obj)
            f = obj.regDesignVariable{1:end-1};
        end

        function q = getQuad(obj)
            q = obj.physicalProblem.getQuadrature();
        end

    end

    methods (Access = protected)

        function createOrientationUpdater(obj)
            cParams.type = 'MinimumEigenValue';
            obj.orientationUpdater = OrientationUpdater.create(cParams);
        end

        function createEquilibriumProblem(obj,fileName)
            a.fileName = fileName;
            s = FemDataContainer(a);
            obj.physicalProblem = FEM.create(s);
            obj.initPrincipalDirections();
        end

        function updateHomogenizedMaterialProperties(obj)
            obj.filterDesignVariable();
            obj.homogenizedVariablesComputer.computeCtensor(obj.regDesignVariable);
        end

        function filterDesignVariable(obj)
            nx = length(obj.designVariable.value)/obj.designVariable.nVariables;
            x  = obj.designVariable.value;
            xf = cell(obj.designVariable.nVariables,1);
            for ivar = 1:obj.designVariable.nVariables
                i0 = nx*(ivar-1) + 1;
                iF = nx*ivar;
                xs = x(i0:iF);
                xf{ivar} = obj.filter.getP0fromP1(xs);
            end
            xf{ivar+1} = obj.designVariable.alpha;
            obj.regDesignVariable = xf;
        end

        function filterGradient(obj)
            g = obj.gradient;
            gf = zeros(size(obj.Msmooth,1),obj.nVariables);
            for ivar = 1:obj.nVariables
                gs = g(:,:,ivar);
                gf(:,ivar) = obj.filter.getP1fromP0(gs);
            end
            gf = obj.Msmooth*gf;
            g = gf(:);
            obj.gradient = g;
        end

        function v = getPdeVariableToPrint(obj,p)
            cParams.physicalProblem = p;
            g = PdeVariableToPrintGetter(cParams);
            v = g.compute();
        end

        function fP = addHomogPrintVariablesNames(obj,fP)
            fH = obj.homogenizedVariablesComputer.createPrintVariables();
            nP = numel(fP);
            for i = 1:numel(fH)
                fP{nP+i} = fH{i};
            end
        end

      function updateAlpha(obj)
          if isequal(obj.designVariable.type,'MicroParams')
            if isfield(obj.physicalProblem.variables,'principalStress')
            cParams.pD = obj.physicalProblem.variables.principalDirections;
            cParams.pS = obj.physicalProblem.variables.principalStress;
            obj.orientationUpdater.compute(cParams);
            alpha = obj.orientationUpdater.alpha;
            obj.designVariable.alpha = alpha;
            end
          end
      end

      function fP = addHomogVariables(obj,fP)
          fH = obj.homogenizedVariablesComputer.addPrintableVariables(obj.designVariable);
          for i = 1:numel(fH)
              fP{end+1} = fH{i};
          end
      end

    end

    methods (Access = private)

        function s = createFEMparameters(obj, file)
            gidParams = obj.createGiDparameters(file);
            s.dim       = gidParams.pdim;
            s.type      = gidParams.ptype;
            s.scale     = gidParams.scale;
            s.mesh      = gidParams.mesh;
            s.dirichlet = gidParams.dirichlet;
            s.pointload = gidParams.pointload;
        end

        function gidParams = createGiDparameters(obj,file)
            gidReader = FemInputReader_GiD();
            gidParams = gidReader.read(file);
        end

        function initPrincipalDirections(obj)
            if isempty(obj.designVariable.alpha)
                dim = obj.physicalProblem.getDimensions();
                nelem = size(obj.dvolu,1);
                ndim = dim.ndimf; %dim.ndim
                alpha0 = zeros(ndim,nelem);
                alpha0(1,:) = 1;
                obj.designVariable.alpha = alpha0;
            end
            obj.physicalProblem.variables.principalDirections = obj.designVariable.alpha;
        end

    end

    methods (Access = protected, Abstract)
        computeFunctionValue(obj)
        computeGradientValue(obj)
        solveState(obj)
        solveAdjoint(obj)
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan