text:
function [ctHandle,cMap,window] = matRad_plotCtSlice(axesHandle,ctCube,cubeIdx,plane,slice,cMap,window)
% matRad function that generates the plot for the CT in the GUI
% The function can also be used in personal matlab figures by passing the
% corresponding axes handle
%
% call
%   [ctHandle,cMap,window] = matRad_plotCtSlice(axesHandle,ctCube,cubeIdx,plane,slice)
%   [ctHandle,cMap,window] = matRad_plotCtSlice(axesHandle,ctCube,cubeIdx,plane,slice,cMap)
%   [ctHandle,cMap,window] = matRad_plotCtSlice(axesHandle,ctCube,cubeIdx,plane,slice,window)
%   [ctHandle,cMap,window] = matRad_plotCtSlice(axesHandle,ctCube,cubeIdx,plane,slice,cMap,window)
%
% input
%   axesHandle  handle to axes the slice should be displayed in
%   ctCube      the cell of ct cubes
%   cubeIdx     Index of the desired cube in the ct struct
%   plane       plane view (coronal=1,sagittal=2,axial=3)
%   slice       slice in the selected plane of the 3D cube
%   cMap        optional argument defining the colormap, default is bone
%               if you want to use the default map with the window argument
%               you can use an empty array []
%   window      optional argument defining the displayed range. default is
%               [min(ctCube(:)) max(ctCube(:))]
%
% output
%   ctHandle    handle of the plotted CT axes
%   cMap        used colormap (same as input if set)
%   window      used window (same as input if set)
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

%Use default colormap?
if nargin < 6 || isempty(cMap)
    cMap = bone(64);
end

if nargin < 7 || isempty(window)
    window = [min(ctCube{cubeIdx}(:)) max(ctCube{cubeIdx}(:))];
end

cMapScale = size(cMap,1) - 1;

%Prepare the slice and convert it to uint8
if plane == 1 % Coronal plane
	ctIndexed = uint8(cMapScale*(squeeze((ctCube{cubeIdx}(slice,:,:)-window(1))/(window(2) - window(1)))));
elseif plane == 2 % sagittal plane
    ctIndexed = uint8(cMapScale*(squeeze((ctCube{cubeIdx}(:,slice,:)-window(1))/(window(2) - window(1)))));
elseif plane == 3 % Axial plane
    ctIndexed = uint8(cMapScale*(squeeze((ctCube{cubeIdx}(:,:,slice)-window(1))/(window(2) - window(1)))));
else
	matRad_cfg.dispError('Invalid plane ''%d'' selected for visualization!',plane);
end

%This circumenvents a bug in Octave when the index in the image hase the maximum value of uint8
if matRad_cfg.isOctave
	ctIndexed(ctIndexed == 255) = 254;
end

ct_rgb = ind2rgb(ctIndexed,cMap);

ctHandle = image('CData',ct_rgb,'Parent',axesHandle);

end

repo: matRad


text:
function [ y, j, ierror ] = yj_check_julian ( y, j )

%*****************************************************************************80
%
%% YJ_CHECK_JULIAN checks a Julian YJ date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, the YJ date.
%
%    Output, integer IERROR, is 0 if no serious error was found in
%    the date, and 1 otherwise.
%

%
%  Check the year.
%
  [ y, ierror ] = y_check_julian ( y );

  if ( ierror ~= 0 )
    return
  end
%
%  Make sure J is not too small or too big.
%
  [ y, j ] = j_borrow_julian ( y, j );

  [ y, j ] = j_carry_julian ( y, j );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function prob_test1486 ( )

%*****************************************************************************80
%
%% TEST1486 tests UNIFORM_01_MEAN, *_SAMPLE, *_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1486\n' );
  fprintf ( 1, '  For the Uniform 01 PDF:\n' );
  fprintf ( 1, '  UNIFORM_01_MEAN computes mean;\n' );
  fprintf ( 1, '  UNIFORM_01_SAMPLE samples;\n' );
  fprintf ( 1, '  UNIFORM_01_VARIANCE computes variance.\n' );

  mean = uniform_01_mean ( );
  variance = uniform_01_variance ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF mean =            %14f\n', mean );
  fprintf ( 1, '  PDF variance =        %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = uniform_01_sample ( seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function [tt]=tt_qreshape(tt,s,sz)

% Reshapes the cores of the input (Q)TT representation.
% s (scalar) is the number of mode indices in each core of the input representation;
% sz is a matrix of size d x s1, where d is the number of cores and s1 is the number of mode indices in each core of the output representation;
% sz(k,alpha) is the mode size of the alpha-th index in the k-th core of the output representation.
% k-th component of prod(sz,2) is the overall numbers of mode degrees of freedom of the k-th core of the output representation
%
% April 26, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%

d=size(tt,1);
for k=1:d
	szk=size(tt{k});
	szk=[szk,ones(1,s)];
	szkr=szk(s+1:numel(szk));
	tt{k}=reshape(tt{k},[sz(k,:),szkr]);
end

return
end

author: oseledets

repo: TT-Toolbox


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp( -1 * sum((x1 - x2).^2) / (2 * (sigma ^ 2)));

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

end

author: gopaczewski

repo: coursera-ml


text:
function demoRevolutionSurface
%DEMOREVOLUTIONSURFACE Demo of revolutionSurface
%
%   Example
%     demoRevolutionSurface
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-04-20
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Draw a torus with vertical axis as revolution axis

circle  = circleToPolygon([10 0 3], 50);
[x, y, t] = revolutionSurface(circle, linspace(0, 4*pi/3, 50));

figure;
surf(x, y, t);
axis equal;

%% Draw a torus with horizontal axis as revolution axis

circle  = circleToPolygon([0 10 3], 50);
[x, y, t] = revolutionSurface(circle, [0 0 1 0], linspace(0, pi/3, 50));

figure;
surf(x, y, t);
axis equal;
drawnow;

end

%% Inner function to avoid call to geom2d module
function varargout = circleToPolygon(circle, varargin)
%CIRCLETOPOLYGON Convert a circle into a series of points
%
%   P = circleToPolygon(CIRC, N);
%   Converts the circle CIRC into an array of  (N+1)-by-2 of double,
%   containing x and y positions of vertices.
%   CIRC is given as [x0 y0 r], where x0 and y0 are coordinate of center,
%   and r is the radius.
%   The resulting polygon is closed (first and last vertices are the same).
%
%   P = circleToPolygon(CIRCLE);
%   uses a default value of N=64 vertices.
%
%   Example
%   circle = circleToPolygon([10 0 5], 16);
%   figure;
%   drawPolygon(circle);
%
%   See also:
%   circles2d, polygons2d, circleArcToPolyline, ellipseToPolygon
%
%
% ---------
% author : David Legland
% created the 06/04/2005.
% Copyright 2010 INRA - Cepia Software Platform.
%

% HISTORY
% 2007-04-20 return a closed polygon with N+1 vertices, use default N=64
% 2011-12-09 rename to 'circleToPolygon'

% determines number of points
N = 64;
if ~isempty(varargin)
    N = varargin{1};
end

% create circle
t = linspace(0, 2*pi, N+1)';
x = circle(1) + circle(3) * cos(t);
y = circle(2) + circle(3) * sin(t);

if nargout == 1
    varargout{1} = [x y];
elseif nargout == 2
    varargout{1} = x;
    varargout{2} = y;
end

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function d = divisor(n)
%% divisor : provides a list of integer divisors of a number.
% divisor(n) : row vector of all distinct divisors of a positive integer N,
%               including 1 and N.
%
% Remark:
%   This function uses the default factor() routine in Matlab and hence is
% limited to input values upto 2^32. However if factor() routine does get
% updated for larger integers, this function will still work fine.
%   Using factor() provides a significant speed improvement over manually
% seaching for the each divisor of n.
%
% Example:
%   a = divisor(12);
%   returns -> a = [1, 2, 3, 4, 6, 12];
%
% See Also:
%   factor, primes

% Author: Yash Kochar ( yashkochar@yahoo.com )
% Last modified: 21st June 2009
%-------------------------------------------------------------------------------

%% Input error check :
%   Check whether input is positive integer and scalar.
if ~isscalar(n)
    error('divisor:NonScalarInput','Input must be a scalar.');
end
if (n < 1) || (floor(n) ~= n)
  error('divisor:PositiveIntegerOnly', 'Input must be a positive integer.');
end

%% Find prime factors of number :
pf = factor(n);         % Prime factors of n
upf = unique(pf);       % Unique

%% Calculate the divisors
d = upf(1).^(0:1:sum(pf == upf(1)))';
for f = upf(2:end)
    d = d*(f.^(0:1:sum(pf == f)));
    d = d(:);
end
d = sort(d)';   % To further improve the speed one may remove this sort command
                %   Just remember to take the transpose of "d" to get a result
                %   as a row vector instead of a column vector.

author: goGPS-Project


text:
function [index, new_cutoff] = find_cutoff_index(cutoff, t2_vals)
%-----------------------------------------------------------------------------------
%
%  [index, new_cutoff] = find_cutoff_index(cutoff, t2_vals)
%
%  Function to find index corresponding to closest T2 value to cutoff in range
%
%	Ives Levesque, Feb 2007
%
%-----------------------------------------------------------------------------------

if cutoff > max(t2_vals) | cutoff < min(t2_vals)
   error('Cutoff value not in range.')
end

index = find(t2_vals==cutoff);

if isempty(index)
    temp_index = max(find(t2_vals<cutoff));
    if abs(t2_vals(temp_index+1)-cutoff) < abs(t2_vals(temp_index)-cutoff)
        index = temp_index + 1;
    else
        index = temp_index;
    end
end

new_cutoff = t2_vals(index);

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

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [f,residuals] = corelucydenoiseNew(Y,H,DAMPAR22,wI,READOUT,SUBSMPL,idx,vec,num)
%CORELUCY Accelerated Damped Lucy-Richarson Operator.
%  Calculates function that when used with the scaled projected array
%  produces the next iteration array that maximizes the likelihood that
%  the entire suite satisfies the Poisson statistics.
%
% See also DECONVLUCY and DECONVBLIND.

%  Copyright 1993-2003 The MathWorks, Inc.
%  $Revision: 1.2.4.2 $

%   References
%   ----------
%   "Acceleration of iterative image restoration algorithms, by D.S.C. Biggs
%   and M. Andrews, Applied Optics, Vol. 36, No. 8, 1997.
%   "Deconvolutions of Hubble Space Telescope Images and Spectra",
%   R.J. Hanisch, R.L. White, and R.L. Gilliland. in "Deconvolution of Images
%   and Spectra", Ed. P.A. Jansson, 2nd ed., Academic Press, CA, 1997.

ReBlurred = real(ifftn(H.*fftn(Y)));

% 1. Resampling if needed
if SUBSMPL ~= 1,% Bin ReBlurred back to the sizeI for non-singleton dims

  %1.Reshape so that the-to-binned dimension separates into two
  %dimensions, with one of them consisting of elements of a single bin.
  ReBlurred = reshape(ReBlurred,vec);

  %2. Bin (==calculate mean) along the first of the-to-binned dimension,
  %that dimension consists of the bin elements. Reshape to get rid off
  for k = num,% new appeared singleton.
    vec(k) = [];
    ReBlurred = reshape(mean(ReBlurred,k),vec);
  end

end;

% Define Residual (This part we added, Andre Diamant)

Residual=wI-ReBlurred;
%fprintf('Residual calculated and mean equals %.5f \n',mean(Residual(:)))

% Apply the wavelet xform denoising method

Residual=denoise(Residual,3);
fprintf('DONE: Mean residual equals %.5f \n',mean(Residual(:)))

% 2. An Estimate for the next step
ReBlurred = ReBlurred + READOUT;
ReBlurred(ReBlurred == 0) = eps;
AnEstim = (ReBlurred+Residual)./ReBlurred + eps;

% 3. Damping if needed
if DAMPAR22 == 0,% No Damping
  ImRatio = AnEstim(idx{:});
else % Damping of the image relative to DAMPAR22 = (N*sigma)^2
  gm = 10;
  g = (wI.*log(AnEstim)+ ReBlurred - wI)./DAMPAR22;
  g = min(g,1);
  G = (g.^(gm-1)).*(gm-(gm-1)*g);
  ImRatio = 1 + G(idx{:}).*(AnEstim(idx{:}) - 1);
end;

f = fftn(ImRatio);
residuals = sum(Residual(:));

author: mvallieres

repo: radiomics


text:
function out = cfg_run_mkdir(job)

% Make a directory and return its path in out.dir{1}.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_mkdir.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

out.dir{1} = fullfile(job.parent{1}, job.name);
if ~exist(out.dir{1}, 'dir')
    mkdir(out.dir{1});
end;

author: spm


text:
function test_bug731

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

% this function checks problems encountered with the neuromag digital trigger reading

cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/test/bug731/test_bug731.fif');
event       = ft_read_event(cfg.dataset);

% issue confirmed -> updated bug but will wait for Alex' input before proceeding

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p = addMultipliedEqualityCuts(p)

newRows = [];
for i =  p.linears
    for j = 1:p.K.f
        row = p.F_struc(j,:);
        used = find(row(2:end));
        if ~any(p.variabletype(used))
            % Linear constraint
            ok = 1;
            for k = 1:length(used)
                find1 = findrows(p.bilinears(:,2:end),[i used(k)]);
                find2 = findrows(p.bilinears(:,2:end),[used(k) i]);
                if (isempty(find1) & isempty(find2))
                    ok = 0;
                    break
                else
                    monom(k) = unique([find1 find2]);
                end
            end
            if ok
                monom = p.bilinears(monom,1);
                row(1+monom) = row(1+used);
                row(1+used) = 0;
                row(1+i) = row(1);
                row(1)=0;
                newRows = [newRows;row];
            end
        end
    end
end
if ~isempty(newRows)
    p.F_struc = [newRows;p.F_struc];
    p.K.f = p.K.f + size(newRows,1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function k = WalravenValeton_k(Lwa, wv_sigma)
%
%       k = WalravenValeton_k(Lwa, wv_sigma)
%
%
%        Input:
%           -Lwa: world adaptation luminance in cd/m^2
%           -wv_sigma: cd/m^2
%
%        Output:
%           -k: WalravenValeton's k value
%
%     Copyright (C) 2016 Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "A Model of Visual Adaptation for Realistic Image Synthesis"
% 	  by James A. Ferwerda, Sumanta N. Pattanaik, Peter Shirley, Donald P. Greenberg
%     in Proceedings of SIGGRAPH 1996
%

if(~exist('wv_sigma', 'var'))
    wv_sigma = 100; %cd/m^2
end

k = (wv_sigma - Lwa / 4) ./ (wv_sigma + Lwa);

k(k < 0.0) = 0.0;

end

author: banterle


text:
% GPSIMTOOLBOXES Toolboxes for the GPSIM software.

importLatest('netlab');
importLatest('minimize');
importLatest('gpmat');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [center, U, obj_fcn] = fcm(data, cluster_n, options)
%FCM Data set clustering using fuzzy c-means clustering.
%
%   [CENTER, U, OBJ_FCN] = FCM(DATA, N_CLUSTER) finds N_CLUSTER number of
%   clusters in the data set DATA. DATA is size M-by-N, where M is the number of
%   data points and N is the number of coordinates for each data point. The
%   coordinates for each cluster center are returned in the rows of the matrix
%   CENTER. The membership function matrix U contains the grade of membership of
%   each DATA point in each cluster. The values 0 and 1 indicate no membership
%   and full membership respectively. Grades between 0 and 1 indicate that the
%   data point has partial membership in a cluster. At each iteration, an
%   objective function is minimized to find the best location for the clusters
%   and its values are returned in OBJ_FCN.
%
%   [CENTER, ...] = FCM(DATA,N_CLUSTER,OPTIONS) specifies a vector of options
%   for the clustering process:
%       OPTIONS(1): exponent for the matrix U             (default: 2.0)
%       OPTIONS(2): maximum number of iterations          (default: 100)
%       OPTIONS(3): minimum amount of improvement         (default: 1e-5)
%       OPTIONS(4): info display during iteration         (default: 1)
%   The clustering process stops when the maximum number of iterations
%   is reached, or when the objective function improvement between two
%   consecutive iterations is less than the minimum amount of improvement
%   specified. Use NaN to select the default value.
%
%   Example
%       data = rand(100,2);
%       [center,U,obj_fcn] = fcm(data,2);
%       plot(data(:,1), data(:,2),'o');
%       hold on;
%       maxU = max(U);
%       % Find the data points with highest grade of membership in cluster 1
%       index1 = find(U(1,:) == maxU);
%       % Find the data points with highest grade of membership in cluster 2
%       index2 = find(U(2,:) == maxU);
%       line(data(index1,1),data(index1,2),'marker','*','color','g');
%       line(data(index2,1),data(index2,2),'marker','*','color','r');
%       % Plot the cluster centers
%       plot([center([1 2],1)],[center([1 2],2)],'*','color','k')
%       hold off;
%
%   See also FCMDEMO, INITFCM, IRISFCM, DISTFCM, STEPFCM.

%   Roger Jang, 12-13-94, N. Hickey 04-16-01
%   Copyright 1994-2018 The MathWorks, Inc.

if nargin ~= 2 && nargin ~= 3
	error(message("fuzzy:general:errFLT_incorrectNumInputArguments"))
end

data_n = size(data, 1);

% Change the following to set default options
default_options = [2;	% exponent for the partition matrix U
		100;	% max. number of iteration
		1e-5;	% min. amount of improvement
		1];	% info display during iteration

if nargin == 2
	options = default_options;
else
	% If "options" is not fully specified, pad it with default values.
	if length(options) < 4
		tmp = default_options;
		tmp(1:length(options)) = options;
		options = tmp;
	end
	% If some entries of "options" are nan's, replace them with defaults.
	nan_index = find(isnan(options)==1);
	options(nan_index) = default_options(nan_index);
	if options(1) <= 1
		error(message("fuzzy:general:errFcm_expMustBeGtOne"))
	end
end

expo = options(1);		% Exponent for U
max_iter = options(2);		% Max. iteration
min_impro = options(3);		% Min. improvement
display = options(4);		% Display info or not

obj_fcn = zeros(max_iter, 1);	% Array for objective function

U = initfcm(cluster_n, data_n);			% Initial fuzzy partition
% Main loop
for i = 1:max_iter
	[U, center, obj_fcn(i)] = stepfcm(data, U, cluster_n, expo);
	if display
		fprintf('Iteration count = %d, obj. fcn = %f\n', i, obj_fcn(i));
	end
	% check termination condition
	if i > 1
		if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro, break; end
	end
end

iter_n = i;	% Actual number of iterations
obj_fcn(iter_n+1:max_iter) = [];

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%NEWFIG Create new figure on given position
%
%	NEWFIG(FIGURE_NUMBER,FIGURE_PER_ROW)
%
% INPUT
%   FIGURE_NUMBER   Number of the figure
%   FIGURE_PER_ROW  Figures per row (default: 4)
%
% OUTPUT
%
% DESCRIPTION
% Creates figure number FIGURE_NUMBER and places it on the screen,
% such that (when sufficient figures are created), the screen is
% covered by an array of figures, with FIGURE_PER_ROW figures per row.
%
% SEE ALSO (PRTools Guide)
% SCATTERD, PLOTC

% $Id: newfig.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function newfig(fign,n);

	if nargin < 2, n=4; end
	if nargin < 1, fign=gcf+1; end

	% if the figure already exist, kill it!
	if any(get(0,'children') == fign)
		delete(fign);
	end

	% Create the figure,
	figure(fign);
	set(gcf,'menubar','none');
	% and place it in the correct position:
	ny = 1-ceil(fign/n)/n;
	nx = (fign -1 - n*floor((fign-1)/n))/n;
	d = 1/n;
	set(gcf,'units','normalized','position',[nx ny d*0.95 d*0.85]);

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function model = cgpsimExpandParam(model, params)

% CGPSIMEXPANDPARAM Expand params into model structure.
% FORMAT
% DESC this is just a wrapper for learning three GP models
% simultaneously.
% ARG model : input model from which parameters will be extracted.
% ARG params : parameters to place in the model.
% RETURN model : model filled with the given parameters.
%
% SEEALSO : gpsimExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2006

% SHEFFIELDML

for i = 1:length(model.comp)
  model.comp{i} = gpsimExpandParam(model.comp{i}, params);
end


text:
function [Gc,Kp,Ti,Td,H]=optPID(key,typ,vars)
k=vars(1); L=vars(2); T=vars(3); N=vars(4);
Td=[];H=1;
if length(vars)==5, iC=vars(5); tt=0;
else,
   Kc=vars(5); Tc=vars(6); kappa=vars(7); tt=1;
end
if tt==0
   if key==2
PIDtab=[0.980, 0.712, 0.569, 1.072, 0.786, 0.628;
       -0.892,-0.921,-0.951,-0.560,-0.559,-0.583;
        0.690, 0.968, 1.023, 0.648, 0.883, 1.007;
       -0.155,-0.247,-0.179,-0.114,-0.158,-0.167];
   elseif key==3
PIDtab=[1.048, 1.042, 0.968, 1.154, 1.142, 1.061;
       -0.897,-0.897,-0.904,-0.567,-0.579,-0.583;
        1.195, 0.987, 0.977, 1.047, 0.919, 0.892;
       -0.368,-0.238,-0.253,-0.220,-0.172,-0.165;
        0.489, 0.385, 0.316, 0.490, 0.384, 0.315;
        0.888, 0.906, 0.892, 0.708, 0.839, 0.832];
   elseif key==4
PIDtab=[1.260, 1.053, 0.942, 1.295, 1.120, 1.001;
       -0.887,-0.930,-0.933,-0.619,-0.625,-0.624;
        0.701, 0.736, 0.770, 0.661, 0.720, 0.754;
       -0.147,-0.126,-0.130,-0.110,-0.114,-0.116;
        0.375, 0.349, 0.308, 0.378, 0.350, 0.308;
        0.886, 0.907, 0.897, 0.756, 0.811, 0.813];
   end
   ii=0; if (L/T>1) ii=3; end; tt=L/T;
   a1=PIDtab(1,ii+iC); b1=PIDtab(2,ii+iC);
   a2=PIDtab(3,ii+iC); b2=PIDtab(4,ii+iC);
   Kp=a1/k*tt^b1; Ti=T/(a2+b2*tt);
   if key==3| key==4
      a3=PIDtab(5,ii+iC); b3=PIDtab(6,ii+iC);
      Td=a3*T*tt^b3;
   end
else
   if key==2,
      Kp=0.361*Kc; Ti=0.083*(1.935*kappa+1)*Tc;
   elseif key==3,
      Kp=0.509*Kc; Td=0.125*Tc;
      Ti=0.051*(3.302*kappa+1)*Tc;
   elseif key==4,
      Kp=(4.437*kappa-1.587)...
         /(8.024*kappa-1.435)*Kc;
      Ti=0.037*(5.89*kappa+1)*Tc;
      Td=0.112*Tc;
   end
end
if key==2, Gc=tf(Kp*[Ti,1],[Ti,0]);
elseif key==3
   nn=[Kp*Ti*Td*(N+1)/N, Kp*(Ti+Td/N), Kp];
   dd=Ti*[Td/N,1,0]; Gc=tf(nn,dd);
elseif key==4
   Gc=tf(Kp*[Ti,1],[Ti,0]);
   nH=[(1+Kp/N)*Ti*Td, Kp*(Ti+Td/N), Kp];
   dH=Kp*conv([Ti,1],[Td/N,1]); H=tf(nH,dH);
end

author: Sable

repo: mcbench-benchmarks


text:
function hFig = view_spheres(HeadModelFile, ChannelFile, sSubject)
% VIEW_SPHERES: Show all the spheres from a "overlapping spheres" forward model.
%
% USAGE:  hFig = view_spheres(HeadModelFile, ChannelFile, sSubject)

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

%% ===== PARSE INPUTS =====
% Get reference surface
if ~isempty(sSubject.iScalp)
    sScalp = in_tess_bst(sSubject.Surface(sSubject.iScalp).FileName);
else
    sScalp = [];
end
if ~isempty(sSubject.iInnerSkull)
    sInnerSkull = in_tess_bst(sSubject.Surface(sSubject.iInnerSkull).FileName);
elseif ~isempty(sSubject.iCortex)
    CortexFile = sSubject.Surface(sSubject.iCortex).FileName;
    sInnerSkull = tess_envelope(CortexFile, 'convhull', 1082, .003);
else
    sInnerSkull = [];
end
hFig = [];

%% ===== LOAD DATA =====
% Load HeadModel ('Param' field only)
HeadModelMat = in_bst_headmodel(HeadModelFile, 0, 'Param', 'MEGMethod', 'EEGMethod');
Spheres = HeadModelMat.Param;
% Load Channels
ChannelMat = in_bst_channel(ChannelFile);
Channel = ChannelMat.Channel;

%% ===== SELECT SPHERES TO DISPLAY =====
iSel = [];
if ~isempty(HeadModelMat.MEGMethod)
    iMeg = good_channel(Channel, [], 'MEG');
    if ismember(HeadModelMat.MEGMethod, {'meg_sphere', 'singlesphere'})
        iSel = [iSel, iMeg(1)];
    elseif ismember(HeadModelMat.MEGMethod, {'os_meg', 'localspheres'})
        iSel = [iSel, iMeg];
    end
end
if ismember(HeadModelMat.EEGMethod, {'eeg_3sphereberg', 'singlesphere', 'concentricspheres'})
    iEeg = good_channel(Channel, [], 'EEG');
    iSel = [iSel, iEeg(1)];
    Channel(iEeg(1)).Name = 'EEG';
end
Spheres = Spheres(iSel);
Channel = Channel(iSel);

%% ===== DISPLAY SURFACES =====
hFig = [];
% Display scalp
if ~isempty(sScalp)
    hFig = view_surface_matrix(sScalp.Vertices, sScalp.Faces, .9, []);
end
% Display inner skull (or inflated cortex envelope)
if ~isempty(sInnerSkull)
    if isempty(hFig)
        hFig = view_surface_matrix(sInnerSkull.Vertices, sInnerSkull.Faces, .2, []);
    else
        view_surface_matrix(sInnerSkull.Vertices, sInnerSkull.Faces, .2, [], hFig);
    end
end
% No figure create: error
if isempty(hFig)
    error('No reference surface available');
end
% Set orientation: left
figure_3d('SetStandardView', hFig, 'left');
% Update figure name
set(hFig, 'Name', ['Check spheres: ' HeadModelFile]);
% Get axes handles
hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');

%% ===== DISPLAY SPHERES+SENSORS =====
% Hack keyboard callback
KeyPressFcn_bak = get(hFig, 'KeyPressFcn');
set(hFig, 'KeyPressFcn', @KeyPress_Callback);
% Create legend (Sphere)
hLabel = uicontrol('Style',               'text', ...
                   'String',              '...', ...
                   'Units',               'Pixels', ...
                   'Position',            [6 0 400 18], ...
                   'HorizontalAlignment', 'left', ...
                   'FontUnits',           'points', ...
                   'FontSize',            bst_get('FigFont'), ...
                   'ForegroundColor',     [.3 1 .3], ...
                   'BackgroundColor',     [0 0 0], ...
                   'Parent',              hFig);
% Plot sensors
if length(Channel) > 10
    markersLocs   = cell2mat(cellfun(@(c)c(:,1), {Channel.Loc}, 'UniformOutput', 0))';
    hSensors = figure_3d('PlotSensorsNet', hAxes, markersLocs, 0, 0);
    set(hSensors, 'LineWidth', 1, 'MarkerSize', 2);
end

% Base of spheres
[X,Y,Z] = sphere(20);
XYZ = tess_sphere(812);
[TH,PHI,R] = cart2sph(X,Y,Z);

% Current sphere
iSphere = 1;
% Draw first sphere
DrawSphere();

%% ===== KEYBOARD CALLBACK =====
    function KeyPress_Callback(h, keyEvent)
        switch (keyEvent.Key)
            % === LEFT, RIGHT, PAGEUP, PAGEDOWN : Processed by TimeWindow  ===
            case {'leftarrow', 'space', 'uparrow'}
                iSphere = iSphere - 1;
            case 'pagedown'
                iSphere = iSphere - 10;
            case {'rightarrow', 'downarrow'}
                iSphere = iSphere + 1;
            case 'pageup'
                iSphere = iSphere + 10;
            otherwise
                KeyPressFcn_bak(h, keyEvent);
                return;
        end
        % Redraw sphere
        if (iSphere <= 0)
            iSphere = length(Spheres);
        end
        if (iSphere > length(Spheres))
            iSphere = 1;
        end
        DrawSphere();
    end

%% ===== DRAW CURRENT SPHERE =====
    function DrawSphere()
        % Delete previous spheres and sensors
        delete(findobj(hAxes, '-depth', 1, 'Tag', 'OverlapSphere'));
        % Draw current sphere
        [X,Y,Z] = sph2cart(TH, PHI, R * Spheres(iSphere).Radii(end));
        hSphere = patch(surf2patch(X + Spheres(iSphere).Center(1),...
                                   Y + Spheres(iSphere).Center(2),...
                                   Z + Spheres(iSphere).Center(3),...
                                   Z + Spheres(iSphere).Center(3)));
        set(hSphere, 'Parent',    hAxes, ...
                     'Facecolor', 'none', ...
                     'EdgeColor', [.8 .8 .8], ...
                     'EdgeAlpha', 1, ...
                     'LineWidth', 1, ...
                     'Tag',       'OverlapSphere');

        % Remove previous selected sensor
        delete(findobj(hAxes, '-depth', 1, 'Tag', 'SelSphereChannel'));
        % Plot selected sensor
        if ~isempty(Channel(iSphere).Loc) && ~ismember(Channel(iSphere).Name, {'EEG','MEG','MEG MAG', 'MEG GRAD'})
            line(Channel(iSphere).Loc(1,1), Channel(iSphere).Loc(2,1), Channel(iSphere).Loc(3,1), ...
                 'Parent',          hAxes, ...
                 'LineWidth',       2, ...
                 'LineStyle',       'none', ...
                 'Marker',          'o', ...
                 'MarkerFaceColor', [1 0 0], ...
                 'MarkerEdgeColor', [.4 .4 .4], ...
                 'MarkerSize',      8, ...
                 'Tag',             'SelSphereChannel');
        end

         % Update legend
        newLegend = sprintf('Sphere #%d/%d  (%s)', iSphere, length(Spheres), Channel(iSphere).Name);
        if (iSphere == 1) && (length(Spheres) > 1)
            newLegend = [newLegend, '       [Press arrows for next/previous sphere...]'];
        end
        set(hLabel, 'String', newLegend);
    end

end

author: brainstorm-tools


text:
function SaveDNN(layer, para, LOG, modelfile)
layer = clean_network_layer(layer);
if para.useGPU
    for layer_i=1:length(layer)
        layer_fields = fields(layer{layer_i});
        for i=1:length(layer_fields)
            layer{layer_i}.(layer_fields{i}) = gather(layer{layer_i}.(layer_fields{i}));
        end
    end
    for i=1:length(para.preprocessing)
        para.preprocessing{i} = FeaturePipe2Memory(para.preprocessing{i});
    end
end
save(modelfile, 'layer', 'para', 'LOG');

end

author: singaxiong

repo: SignalGraph


text:
% compile mex file
mex mex/perform_moment_transform.cpp mex/perform_moment_transform_mex.cpp

author: gpeyre

repo: matlab-toolboxes


text:
classdef nme_gen3p_acc < mp.nme_gen3p & mp.form_acc

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

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% pdfbcreate_equiv_models
% written by:   Duncan Po
% Date:         December 5/2002
% Usage:        models = pdfbcreate_equiv_models(model)
% Input:    model   -   original model
% Output:   models  -   an array of equivalent models
%
% create equivalent models by flipping the states

function models = pdfbcreate_equiv_models(model)

l = 1;
temp_model{l}{1} = model;
for j = 1:length(model.stdv)
    for k= 1:length(model.stdv{j})
        l = l+1;
        for mm = 1:2.^(l-1)
            if mod(mm, 2) == 1
                temp_model{l}{mm} = pdfbflip_model(temp_model{l-1}{ceil(mm/2)}, j, k);
            else
                temp_model{l}{mm} = temp_model{l-1}{ceil(mm/2)};
            end;
        end;
    end;
end;

num = length(temp_model{end});
for ddd = 1:num
    models{ddd} = temp_model{end}{ddd};
end;


text:
function [displacements,rois_dic,seedinfo,outstate] = ncorr_alg_dicanalysis(imgs,radius,spacing,cutoff_diffnorm,cutoff_iteration,total_threads,enabled_stepanalysis,subsettrunc,num_img,total_imgs,pos_parent,params_init)
% This function performs RG-DIC.
%
% Inputs -----------------------------------------------------------------%
%   imgs - struct; Contains struct('imginfo',{},'roi',{}). Contains the
%   reference and current images packaged together. The imginfo field
%   contains an ncorr_class_img and the ROI field contains an ncorr_class_roi.
%   radius - integer; subset radius
%   spacing - integer; subset spacing
%   cutoff_diffnorm - double; cutoff for the norm of difference vector
%   cutoff_iteration - integer; cutoff for the number of IC-GN iterations
%   total_threads - integer; total number of threads
%   enabled_stepanalysis - logical; if true, then process as many seeds as
%   possible. If false, process all the seeds.
%   subsettrunc - logical; if true, then subset truncation is enabled
%   num_img - integer; reference image number
%   total_imgs - integer; total number of images
%   pos_parent - integer array; this is the position of the parent figure
%   which determines where to position this figure
%   params_init - struct; contains struct('paramvector',{},'num_region',{},
%   'num_thread',{},'computepoints',{}). If its not empty, it contains last
%   set of seeds used for the previous iteration
%
% Outputs ----------------------------------------------------------------%
%   displacements - struct; contains
%   struct('plot_u',{},'plot_v',{},'plot_corrcoef',{},'plot_validpoints',{})
%   rois_dic - ncorr_class_img; ROIs after being unioned with reference roi
%   and validpoints from DIC analysis.
%   seedinfo - struct; contains struct('paramvector',{},'num_region',{},
%   'num_thread',{},'computepoints',{})
%   outstate - integer; returns either out.cancelled, out.failed, or
%   out.success.
%
% outstate will only return failed if ncorr_alg_rgdic throws an exception.

    % Initialize outputs
    outstate = out.cancelled;
    displacements = struct('plot_u',{},'plot_v',{},'plot_corrcoef',{},'plot_validpoints',{});
    rois_dic = ncorr_class_roi.empty;
    seedinfo = struct('paramvector',{},'num_region',{},'num_thread',{},'computepoints',{});

    % --------------------------------------------------------------------%
    % Get Seeds ----------------------------------------------------------%
    % --------------------------------------------------------------------%

    if (isempty(params_init))
        % Get new seeds manually through GUI. outstate_seeds will either be
        % success or cancelled.
        [seedinfo_prelim,threaddiagram,outstate_seeds] = ncorr_gui_seedanalysis(imgs(1).imginfo, ...
                                                                                [imgs(2:end).imginfo], ...
                                                                                imgs(1).roi, ...
                                                                                radius, ...
                                                                                spacing, ...
                                                                                cutoff_diffnorm, ...
                                                                                cutoff_iteration, ...
                                                                                total_threads, ...
                                                                                enabled_stepanalysis, ...
                                                                                subsettrunc, ...
                                                                                num_img, ...
                                                                                total_imgs, ...
                                                                                pos_parent);
        % See if analysis was cancelled
        if (outstate_seeds ~= out.success)
            return;
        end
    else
        % Use params init to place seeds directly
        % Initialize seedinfo_prelim
        seedinfo_prelim = struct('paramvector',{},'num_region',{},'num_thread',{},'computepoints',{});

        % Initialize Thread diagram buffers - some buffers are unused
        ref_reduced = imgs(1).imginfo.reduce(spacing); % Unused
        roi_reduced = imgs(1).roi.reduce(spacing); % Unused
        preview_threaddiagram = zeros(size(roi_reduced.mask)); % Unused
        threaddiagram_buffer = -ones(size(roi_reduced.mask));
        threaddiagram = -ones(size(roi_reduced.mask)); % This must be negative 1

        % Calculate seeds - cycle over regions
        num_imgs_success = inf; % Initialize to impossibly high value
        manualseed = false; % Set this parameter true if there is trouble. This causes the user to reselect the seeds manually through a GUI.
        % Cycle over regions
        for i = 0:size(params_init,2)-1
            % Get updated seed positions based on seed locations of the
            % last successfully seeded image for this region
            % Initialize
            pos_seed = zeros(size(params_init,1),2);
            % Cycle over threads
            for j = 0:size(params_init,1)-1
                % Must round displacements so that the updated position
                % lies properly on "spaced grid" - THIS IS IMPORTANT!!!
                pos_seed(j+1,1) = params_init(j+1,i+1,end).paramvector(1) + round(params_init(j+1,i+1,end).paramvector(3)/(spacing+1))*(spacing+1);
                pos_seed(j+1,2) = params_init(j+1,i+1,end).paramvector(2) + round(params_init(j+1,i+1,end).paramvector(4)/(spacing+1))*(spacing+1);
            end

            % Get num_region - just use seed from the first thread
            num_region = params_init(1,i+1,end).num_region;

            % Make sure seed positions are within the ROI and are unique.
            % If they arent then prompt the user and have them replace the
            % seeds manually.
            regionmask = roi_reduced.get_regionmask(num_region);
            if (size(pos_seed,1) == size(unique(pos_seed,'rows'),1) && ...
                all(pos_seed(:,1)./(spacing+1) >= 0) && all(pos_seed(:,1)./(spacing+1) < size(regionmask,2)) && ...
                all(pos_seed(:,2)./(spacing+1) >= 0) && all(pos_seed(:,2)./(spacing+1) < size(regionmask,1)) && ...
                all(regionmask(sub2ind(size(regionmask),pos_seed(:,2)/(spacing+1)+1,pos_seed(:,1)/(spacing+1)+1))))
                % Get seeds - outstate will be successful if at least one
                % image is seeded
                [seedinfo_buffer,convergence_buffer,outstate_seeds] = ncorr_alg_seedanalysis(imgs(1).imginfo, ...
                                                                                [imgs(2:end).imginfo], ...
                                                                                imgs(1).roi, ...
                                                                                num_region, ...
                                                                                pos_seed, ...
                                                                                radius, ...
                                                                                cutoff_diffnorm, ...
                                                                                cutoff_iteration, ...
                                                                                enabled_stepanalysis, ...
                                                                                subsettrunc, ...
                                                                                num_img, ...
                                                                                total_imgs); %#ok<ASGLU>
                % See if analysis was cancelled
                if (outstate_seeds == out.cancelled)
                    return;
                end

                % See if analysis was successful or failed
                if (outstate_seeds == out.success)
                    % Form thread diagram for these seeds
                    ncorr_alg_formthreaddiagram(threaddiagram_buffer,preview_threaddiagram,int32(pos_seed/(spacing+1)),regionmask,ref_reduced.formatted());

                    % Get compute points
                    for j = 0:size(seedinfo_buffer,3)-1
                        for k = 0:size(seedinfo_buffer,1)-1
                            seedinfo_buffer(k+1,1,j+1).computepoints = length(find(threaddiagram_buffer == k));
                        end
                    end
                else
                    % Seed analysis failed. Alert user and then ask him/her
                    % to manually place seeds.
                    h_error = errordlg('Not a single image was seeded correctly; please replace seeds manually.','Error','modal');
                    uiwait(h_error);

                    manualseed = true;
                end
            else
                % Alert user
                h_error = errordlg('One or more seeds went outside of the ROI **OR** converged on top of each other, please replace manually. If this happens often then dont place seeds near each other or near the boundary.','Error','modal');
                uiwait(h_error);
                manualseed = true;
            end

            if (~manualseed)
                % Take minimum of num_imgs_success and
                % seedinfo. Buffer can be more or less than
                % num_imgs_success
                num_imgs_success = min(num_imgs_success, size(seedinfo_buffer,3));

                % Clear out other images in buffer and prelim
                if (~isempty(seedinfo_prelim))
                    seedinfo_prelim = seedinfo_prelim(:,:,1:num_imgs_success);
                end
                seedinfo_buffer = seedinfo_buffer(:,:,1:num_imgs_success);

                % Append seedinfo_buffer - append along 2nd dimension
                seedinfo_prelim = horzcat(seedinfo_prelim,seedinfo_buffer); %#ok<AGROW>

                % Merge threaddiagram from previous iteration
                threaddiagram(threaddiagram_buffer ~= -1) = threaddiagram_buffer(threaddiagram_buffer ~= -1);
            else
                % If there's a problem, just ask user to manually reset
                % seeds. Outstate will either be success or cancelled.
                [seedinfo_prelim,threaddiagram,outstate_seeds] = ncorr_gui_seedanalysis(imgs(1).imginfo, ...
                                                                                [imgs(2:end).imginfo], ...
                                                                                imgs(1).roi, ...
                                                                                radius, ...
                                                                                spacing, ...
                                                                                cutoff_diffnorm, ...
                                                                                cutoff_iteration, ...
                                                                                total_threads, ...
                                                                                enabled_stepanalysis, ...
                                                                                subsettrunc, ...
                                                                                num_img, ...
                                                                                total_imgs, ...
                                                                                pos_parent);

                % Check if analysis was cancelled
                if (outstate_seeds ~= out.success)
                    return;
                end

                % Break since seeds for all regions are placed in GUI.
                break;
            end
        end

        % Debug ----------------------------------------------------------%
        %{
        figure, imshow(threaddiagram,[]); hold on;
        for i = 0:size(seedinfo_prelim,2)-1
            for j = 0:size(seedinfo_prelim,1)-1
                plot(seedinfo_prelim(j+1,i+1,1).paramvector(1)/(spacing+1),seedinfo_prelim(j+1,i+1,1).paramvector(2)/(spacing+1),'ro');
            end
        end
        hold off;
        %}
        % ----------------------------------------------------------------%
    end

    % --------------------------------------------------------------------%
    % Perform DIC Analysis -----------------------------------------------%
    % --------------------------------------------------------------------%

    % Format Seeds ---------------------------------------------------%
    seedinfo_prelim_f = seedinfo_prelim;
    for i = 0:size(seedinfo_prelim,1)-1
        for j = 0:size(seedinfo_prelim,2)-1
            for k = 0:size(seedinfo_prelim,3)-1
                seedinfo_prelim_f(i+1,j+1,k+1).num_region = int32(seedinfo_prelim(i+1,j+1,k+1).num_region);
                seedinfo_prelim_f(i+1,j+1,k+1).num_thread = int32(seedinfo_prelim(i+1,j+1,k+1).num_thread);
                seedinfo_prelim_f(i+1,j+1,k+1).computepoints = int32(seedinfo_prelim(i+1,j+1,k+1).computepoints);
            end
        end
    end

    % Begin DIC analysis ---------------------------------------------%
    displacements_prelim = struct('plot_u',{},'plot_v',{},'plot_corrcoef',{},'plot_validpoints',{});
    rois_dic_prelim = ncorr_class_roi.empty;
    for i = 0:size(seedinfo_prelim,3)-1
        % Put in try block because DIC can run
        % out of memory
        try
            tic
            % No need to send in the number of regions or total number
            % of threads because this is encoded in the size of the
            % seeds.
            % FYI: seeinfo_prelim(i,j,k) where i refers to the thread
            % number, j refers to the region, and k refers to the image.
            % ncorr_alg_rgdic will either return success or cancelled or
            % return an exception.
            [displacements_prelim(i+1),outstate_dic] = ncorr_alg_rgdic(imgs(1).imginfo.formatted(), ...
                                                                       imgs(i+2).imginfo.formatted(), ...
                                                                       imgs(1).roi.formatted(), ...
                                                                       seedinfo_prelim_f(:,:,i+1), ...
                                                                       int32(threaddiagram), ...
                                                                       int32(radius), ...
                                                                       int32(spacing), ...
                                                                       cutoff_diffnorm, ...
                                                                       int32(cutoff_iteration), ...
                                                                       logical(subsettrunc), ...
                                                                       int32(num_img+i), ...
                                                                       int32(total_imgs));
            toc
        catch %#ok<CTCH>
            % Only time an exception is returned is either if rgdic has
            % incorrect inputs, ran out of memory, or there is a bug in
            % the code. The first option never happens if rgdic
            % is called through the program. The second option can
            % happen and is thus the only real handleable exception
            % which will get thrown by rgdic. The last case was not
            % intended, so it is not handled here.
            h_error = errordlg('Ncorr most likely ran out of memory while performing DIC. Please clear memory in workspace, restart Ncorr, or crop/use smaller images before doing analysis.','Error','modal');
            uiwait(h_error);
            return;
        end

        % See if analysis was cancelled
        if (outstate_dic ~= out.success)
            return;
        end

        % Take union of reference ROI with validpoints to get the new
        % ROI.
        rois_dic_prelim(i+1) = imgs(1).roi.get_union(displacements_prelim(i+1).plot_validpoints,spacing);
    end

    % Set outputs
    for i = 0:length(displacements_prelim)-1
        displacements(i+1) = displacements_prelim(i+1);
        rois_dic(i+1) = rois_dic_prelim(i+1);
    end
    for i = 0:size(seedinfo_prelim,1)-1
        for j = 0:size(seedinfo_prelim,2)-1
            for k = 0:size(seedinfo_prelim,3)-1
                seedinfo(i+1,j+1,k+1) = seedinfo_prelim(i+1,j+1,k+1);
            end
        end
    end
    outstate = out.success;
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [Archive,znad] = UpdateArchive(Archive,N)
% Update the archive in DMOEA-eC

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Detect the non-dominated solutions
    Archive = Archive(NDSort(Archive.objs,1)==1);
    AObj    = Archive.objs;

    %% Select the extreme solutions
    Choose           = false(1,length(Archive));
    [~,extreme1]     = max(AObj,[],1);
    [~,extreme2]     = min(AObj,[],1);
    Choose(extreme1) = true;
    Choose(extreme2) = true;

    %% Select other solutions by truncation
    if sum(Choose) > N
        selected = find(Choose);
        Choose   = selected(randperm(length(selected),N));
    else
        Distance = pdist2(AObj,AObj);
        Distance(logical(eye(length(Distance)))) = inf;
        while sum(Choose) < N && ~all(Choose)
            unSelected = find(~Choose);
            [~,x]      = max(min(Distance(~Choose,Choose),[],2));
            Choose(unSelected(x)) = true;
        end
    end
    Archive = Archive(Choose);

    %% Update the nadir point
    znad = max(Archive.objs,[],1);
end


text:
function zp=nearestPointOnPlane(z,n,p0)
%%NEARESTPOINTONPLANE Given a point z, determine the nearest point on the
%           plane given by the equation n'*(x-p0)=0, where n is a normal to
%           the plane, p0 is a point and x is the free variable. All
%           quantities are real.
%
%INPUTS: z A numDimXnumPoints set of points to project.
%        n A numDimX1 normal to the plane.
%       p0 A numDimX1 point on the plane.
%
%OUTPUTS: zp The numDimXnumPoints set of projections of the points in z
%            onto the plane.
%
%The expression for the nearest point to a plane is given in Chapter 12.1
%of [1].
%
%EXAMPLE:
% n=[1;2;3];
% p0=[0;-1;2];
% z=[4;4;4];
% zp=nearestPointOnPlane(z,n,p0)
% %One will get zp=[2.571428571428572;1.142857142857143;-0.285714285714286].
% %It can be verified that
% n'*(zp-p0)
% %is zero within finite precision limits and
% cross(n,(z-zp))
% %is also zero within finite precision limits.
%
%REFERENCES:
%[1] P. J. Schneider and D. H. Eberly, Geometric Tools for Computer
%    Graphics. Amsterdam: Morgan Kaufmann Publishers, 2003.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

d=-n'*p0;
%The equation for the plane is now:
%x'*n+d=0;

u=n/(n'*n);

%For a single point z, the following line is equivalent to
%zp=z-(z'*n+d)*u;
zp=bsxfun(@minus,z,bsxfun(@times,sum(bsxfun(@times,z,n),1)+d,u));

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 [ ntree, itree, seed ] = tree_rooted_random ( nnode, seed )

%*****************************************************************************80
%
%% TREE_ROOTED_RANDOM selects a random unlabeled rooted tree.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    FORTRAN90 version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer NNODE, the number of nodes.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, integer NTREE(NNODE).  NTREE(I) is the number of
%    rooted, unlabeled trees on I nodes, for I = 1, 2, ... NNODE.
%
%    Output, integer ITREE(NNODE).  (I,ITREE(I)) is the I-th edge
%    of the output tree for I = 2,NNODE.  ITREE(1)=0.
%
  if ( nnode <= 0  )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TREE_ROOTED_RANDOM - Fatal error!\n' );
    fprintf ( 1, '  NNODE = %d\n', nnode );
    fprintf ( 1, '  but NNODE must be at least 1.\n' );
    error ( 'TREE_ROOTED_RANDOM - Fatal error!' );
  end
%
%  Compute a table of the number of such trees for a given number of nodes.
%
  ntree = tree_rooted_enum ( nnode );
%
%  Now select one such tree at random.
%
  l = 0;

  nval = nnode;
  is1 = 0;
  is2 = 0;

  while ( 1 )

    while ( 2 < nval )

      [ r, seed ] = r8_uniform_01 ( seed );

      iz = floor ( ( nval - 1 ) * ntree(nval) * r );

      id = 0;

      id = id + 1;
      itd = id * ntree(id);
      m = nval;
      j = 0;

      while ( 1 )

        j = j + 1;
        m = m - id;

        if ( m < 1 )
          id = id + 1;
          itd = id * ntree(id);
          m = nval;
          j = 0;
          continue
        end

        iz = iz - ntree(m) * itd;

        if ( iz < 0 )
          break
        end

      end

      is1 = is1 + 1;
      stack(1,is1) = j;
      stack(2,is1) = id;
      nval = m;

    end

    itree(is2+1) = l;
    l = is2 + 1;
    is2 = is2 + nval;

    if ( 1 < nval )
      itree(is2) = is2 - 1;
    end

    while ( 1 )

      nval = stack(2,is1);

      if ( nval ~= 0 )
        stack(2,is1) = 0;
        break
      end

      j = stack(1,is1);
      is1 = is1 - 1;
      m = is2 - l + 1;
      ll = itree(l);
      ls = l + ( j - 1 ) * m - 1;

      if ( j ~= 1 )
        for i = l : ls
          itree(i+m) = itree(i) + m;
          if ( mod(i-l,m) == 0 )
            itree(i+m) = ll;
          end
        end
      end

      is2 = ls + m;

      if ( is2 == nnode )
        return
      end

      l = ll;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all; close all;

addpath('../matlab/');

model_root_dir = 'DSC_removal_SDR/';
definition_file = [model_root_dir 'deploy.prototxt'];
binary_file = [model_root_dir 'snapshot/' 'DSC_iter_160000_lab.caffemodel'];

assert(exist(binary_file, 'file') ~= 0);
assert(exist(definition_file, 'file') ~= 0);

caffe.reset_all();

caffe.set_mode_gpu();
caffe.set_device(0);

% Initialize a network
net = caffe.Net(definition_file, binary_file, 'test');

%root_dir = '/home/xwhu/dataset/SRD/train_argu/shadow/';
% root_dir = '/home/xwhu/dataset/ISTD/train_argu/train_A/';
%imgFiles=dir([root_dir '*.jpg']);

root_dir = '/home/xwhu/dataset/SRD/test_data/shadow/';
image_list=textread('../../data/SRD/test.txt', '%s');
%root_dir = '/home/xwhu/dataset/ISTD/test/test_A/';
%image_list=textread('../../data/ISTD/test.txt', '%s');

save_root = [model_root_dir 'result/'];

if exist(save_root, 'dir') == 0
    mkdir(save_root);
end

%nImg=length(imgFiles);
nImg=length(image_list);

imgW = 400; imgH = 400;
scale = 0.0039212686;

usedtime = 0;
show = 0;

colorTransform = makecform('srgb2lab');
colorTransform2 = makecform('lab2srgb');

for k = 1 : nImg

    test_image = imread([root_dir image_list{k}]);
    %test_image = imread([root_dir imgFiles(k).name]);

    test_image = applycform(test_image, colorTransform);

    if (show)
        imshow(test_image);
    end

    ori_size = [size(test_image,1), size(test_image,2)];
    test_image = imresize(test_image,[imgH imgW]);
    test_image = single(test_image(:,:,[3 2 1]));
    test_image = test_image.*scale;
    test_image = permute(test_image, [2 1 3]);

    % network forward
    tic; outputs = net.forward({test_image}); pertime=toc;
    usedtime=usedtime+pertime; avgtime=usedtime/k;

    res_fuse = net.blobs('upscore-fuse').get_data();
    res_global = net.blobs('res_g').get_data();
    final = (res_fuse + res_global)./2;

    final = permute(final, [2 1 3]);
    final = final./scale;
    final = uint8(final(:,:,[3 2 1]));
    final = imresize(final, ori_size);

    final = applycform(final, colorTransform2);

    %imwrite(final,[save_root imgFiles(k).name]);
    imwrite(final,[save_root image_list{k}]);

    if (mod(k,100)==0), fprintf('idx %i/%i, avgtime=%.4fs\n',k,nImg,avgtime); end

end

fprintf('idx %i/%i, avgtime=%.4fs\n',k,nImg,avgtime);

author: xw-hu

repo: DSC

save_path: github-repos/MATLAB/xw-hu-DSC


text:
function obj = modelDownsampleRawMovies(obj)
	% Downsamples calcium imaging files and moves them to appropriate folders. Best for Inscopix nVista v1.0 and v2.0 data but should work for any HDF5 file.
	% Biafra Ahanonu
	% inputs
		%
	% outputs
		%

	% changelog
		% 2015.11.25 - Added ability to downsample TIFFs, at the moment assumes that they will be done the inscopix way.
		% 2017.02.18 - slight change in how determination of already existing movies is performed
		% 2019.07.11 [19:53:47] - Added support for ISXD Inscopix files via Inscopix Data Processing Software. Also changed way program is run to loop over files and check for each extension in the case a folder as multiple types of files that need to be downsampled (e.g. HDF5, TIF, ISXD).
		% 2019.11.18 [17:28:48] - Input dialog for folder list or GUI to manually entering multiple folders to downsample. Easier for users than current purely comma separated single line list.
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2021.12.08 [22:11:00] - Additional updates to handle CIAtah v4.0 API switch.
		% 2022.02.25 [15:40:56] - Change TIF finding to split apart the information file before using for regular expression.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	try
		[newFolderListCell] = subfxnSelectFolders();
		if ~isempty(newFolderListCell)
			obj.downsampleRawOptions.folderListInfo = strjoin(newFolderListCell,',');
		end

		downsampleSettings = inputdlg(...
			{...
				'Folder(s) where raw HDF5s are located [from previous prompts]. Use comma to separate multiple source folders:',...
				'Folder to save downsampled HDF5s to:',...
				'Decompression source root folder(s). Use comma to separate multiple source folders:',...
				'Downsample factor (in x-y):',...
				'Regexp for HDF5 (e.g. "recording.*.hdf5"), TIFF (e.g. "recording.*.tif"), ISXD/Inscopix (e.g. ".*.isxd") files:',...
				'HDF5 hierarchy name where movie is stored:',...
				'Max chunk size (MB)',...
				'Regexp for folders:',...
				'Regexp for base filename (use txt log name):',...
				'Extension for base filename:',...
				'Location to store raw HDF5 (leave empty to skip):',...
				'Raw downsample factor (in x-y):',...
				'Output HDF5 dataset name:'
			},...
			'downsample settings',[1 100],...
			{...
				obj.downsampleRawOptions.folderListInfo,...
				obj.downsampleRawOptions.downsampleSaveFolder,...
				obj.downsampleRawOptions.downsampleSrcFolder,...
				obj.downsampleRawOptions.downsampleFactor,...
				obj.downsampleRawOptions.fileFilterRegexp,...
				obj.downsampleRawOptions.datasetName,...
				obj.downsampleRawOptions.maxChunkSize,...
				obj.downsampleRawOptions.srcFolderFilterRegexp,...
				obj.downsampleRawOptions.srcSubfolderFileFilterRegexp,...
				obj.downsampleRawOptions.srcSubfolderFileFilterRegexpExt,...
				obj.downsampleRawOptions.downsampleSaveFolderTwo,...
				obj.downsampleRawOptions.downsampleFactorTwo,...
				obj.downsampleRawOptions.outputDatasetName,...
			});
		obj.downsampleRawOptions.folderListInfo = downsampleSettings{1};
		obj.downsampleRawOptions.downsampleSaveFolder = downsampleSettings{2};
		obj.downsampleRawOptions.downsampleSrcFolder = downsampleSettings{3};
		obj.downsampleRawOptions.downsampleFactor = downsampleSettings{4};
		obj.downsampleRawOptions.fileFilterRegexp = downsampleSettings{5};
		obj.downsampleRawOptions.datasetName = downsampleSettings{6};
		obj.downsampleRawOptions.maxChunkSize = downsampleSettings{7};
		obj.downsampleRawOptions.srcFolderFilterRegexp = downsampleSettings{8};
		obj.downsampleRawOptions.srcSubfolderFileFilterRegexp = downsampleSettings{9};
		obj.downsampleRawOptions.srcSubfolderFileFilterRegexpExt = downsampleSettings{10};
		obj.downsampleRawOptions.downsampleSaveFolderTwo = downsampleSettings{11};
		obj.downsampleRawOptions.downsampleFactorTwo = downsampleSettings{12};
		obj.downsampleRawOptions.outputDatasetName = downsampleSettings{13};

		% dsOpt = obj.downsampleRawOptions;

		folderListInfo = strsplit(downsampleSettings{1},',');
		cellfun(@display,folderListInfo);
		nFolders = length(folderListInfo);
		% ~isempty(regexp(downsampleOptions.fileFilterRegexp,'(.hdf5|.h5)'))
		% ~isempty(regexp(downsampleOptions.fileFilterRegexp,'(.tiff|.tif)'))
		% ============================
		for folderNo = 1:nFolders
			display(repmat('=',1,42))
			display([num2str(folderNo) '\' num2str(nFolders) ': ' folderListInfo{folderNo}])
			% downsample if all decompressed files are in the same folder
			downsampleOptions.downsampleSaveFolder = downsampleSettings{2};
			downsampleOptions.datasetName = downsampleSettings{6};
			downsampleOptions.maxChunkSize = str2num(downsampleSettings{7});
			downsampleOptions.downsampleFactor = str2num(downsampleSettings{4});
			downsampleOptions.fileFilterRegexp = downsampleSettings{5};
			downsampleOptions.outputDatasetName = downsampleSettings{13};

			downsampleOptions.downsampleSaveFolderTwo = downsampleSettings{11};
			downsampleOptions.downsampleFactorTwo = str2num(downsampleSettings{12});

			movieList = getFileList([folderListInfo{folderNo} filesep], downsampleOptions.fileFilterRegexp);
			nMoviesH = length(movieList);

			destFolderAll = {};
			destFolderAll2 = {};
			for movieNo = 1:nMoviesH
				movieName = movieList{movieNo};
				[pathstr,name,ext] = fileparts(movieName);
				fprintf('Downsampling %d/%d: %s\n',movieNo,nMoviesH,movieName);
				if ~isempty(regexp(ext,'(.hdf5|.h5)'))
					display('Downsampling HDF5s...')
					downsampleHDFMovieFxnObj({movieName},downsampleOptions);
				elseif ~isempty(regexp(ext,'(.isxd)'))
					display('Downsampling ISXDs...')
					[downsampleSaveFolderMod, downsampleSaveFolderTwoMod] = downsampleIsxdMovieFxnObj({movieName},downsampleOptions);
					destFolderAll{end+1} = downsampleSaveFolderMod;
					destFolderAll2{end+1} = downsampleSaveFolderTwoMod;
				elseif ~isempty(regexp(ext,'(.tiff|.tif)'))
					if movieNo==1
						display('Downsampling TIFFs...')
						folderPath = folderListInfo{folderNo};
						downsampleOptions.srcSubfolderFileFilterRegexp = downsampleSettings{9};
						downsampleOptions.srcSubfolderFileFilterRegexpExt = downsampleSettings{10};
						downsampleTiffMovieFxnObj(folderPath,downsampleOptions);
					end
				end
			end

			% if ~isempty(regexp(downsampleOptions.fileFilterRegexp,'(.hdf5|.h5)'))
			% 	display('Downsampling HDF5s...')
			% 	downsampleHDFMovieFxnObj(movieList,downsampleOptions);
			% elseif ~isempty(regexp(downsampleOptions.fileFilterRegexp,'(.isxd)'))
			% 	display('Downsampling ISXDs...')
			% 	downsampleIsxdMovieFxnObj(movieList,downsampleOptions);
			% elseif ~isempty(regexp(downsampleOptions.fileFilterRegexp,'(.tiff|.tif)'))
			% 	display('Downsampling TIFFs...')
			% 	folderPath = folderListInfo{folderNo};
			% 	downsampleOptions.srcSubfolderFileFilterRegexp = downsampleSettings{9};
			% 	downsampleOptions.srcSubfolderFileFilterRegexpExt = downsampleSettings{10};
			% 	downsampleTiffMovieFxnObj(folderPath,downsampleOptions);
			% end
			% ioptions.folderListInfo = [folderListInfo{folderNo} filesep];
			% ioptions.downsampleSaveFolder = [downsampleSettings{2} filesep];
			% ioptions.downsampleFactor = str2num(downsampleSettings{4});
			% ioptions.fileFilterRegexp = downsampleSettings{5};
			% ioptions.datasetName = downsampleSettings{6};
			% ioptions.maxChunkSize = str2num(downsampleSettings{7});
			% ioptions.runArg = 'downsampleMovie';
			% ostruct = controllerAnalysis('options',ioptions);
		end
		% re-create folder structure
		clear ioptions
		% run HDF5 or TIFF
		%
		downsampleSrcFolder = strsplit(downsampleSettings{3},',');
		downsampleSaveFolder = [downsampleSettings{2} filesep];
		% used to determine which folders to copy from src to dest
		ioptions.srcFolderFilterRegexp = downsampleSettings{8};
		% this regexp is used to search the destination directory
		ioptions.srcSubfolderFileFilterRegexp = downsampleSettings{9};
		%
		ioptions.srcSubfolderFileFilterRegexpExt = downsampleSettings{10};
		[success destFolders] = moveFilesToFolders(downsampleSrcFolder,char(downsampleSaveFolder(:))','options',ioptions);
		destFolders = unique(destFolders);
		if isempty(destFolderAll)
		else
			destFolders = unique(destFolderAll);
		end
		% ============================
		if ~isempty(downsampleSettings{11})
			% movie second set of movie files
			downsampleSaveFolder = [downsampleSettings{11} filesep];
			% used to determine which folders to copy from src to dest
			ioptions.srcFolderFilterRegexp = downsampleSettings{8};
			% this regexp is used to search the destination directory
			ioptions.srcSubfolderFileFilterRegexp = downsampleSettings{9};
			%
			ioptions.srcSubfolderFileFilterRegexpExt = downsampleSettings{10};
			[success destFolders2] = moveFilesToFolders(downsampleSrcFolder,char(downsampleSaveFolder(:))','options',ioptions);

			if isempty(destFolderAll2)
			else
				destFolders2 = unique(destFolderAll2);
			end
		end
		% ============================
		% if ~isempty(downsampleSettings{11})
		%     % now same the semi-downsampled raw files
		%     for folderNo = 1:nFolders
		%         % downsample if all decompressed files are in the same folder
		%         downsampleOptions.downsampleSaveFolder = [downsampleSettings{11} filesep];
		%         downsampleOptions.datasetName = downsampleSettings{6};
		%         downsampleOptions.maxChunkSize = str2num(downsampleSettings{7});
		%         downsampleOptions.downsampleFactor = str2num(downsampleSettings{12});
		%         downsampleOptions.fileFilterRegexp = downsampleSettings{5};

		%         movieList = getFileList([folderListInfo{folderNo} filesep], downsampleOptions.fileFilterRegexp);

		%         downsampleHDFMovieFxnObj(movieList,downsampleOptions);
		%     end
		%     % re-create folder structure
		%     clear ioptions
		%     %
		%     downsampleSaveFolder = [downsampleSettings{11} filesep];
		%     % used to determine which folders to copy from src to dest
		%     ioptions.srcFolderFilterRegexp = downsampleSettings{8};
		%     % this regexp is used to search the destination directory
		%     ioptions.srcSubfolderFileFilterRegexp = downsampleSettings{9};
		%     %
		%     ioptions.srcSubfolderFileFilterRegexpExt = downsampleSettings{10};
		%     [success destFolders2] = moveFilesToFolders(strsplit(downsampleSettings{3},','),char(downsampleSaveFolder(:))','options',ioptions);
		%     % destFolders = unique(destFolders);
		% end
		% destFolders
		% class(destFolders)

		% add folder info to inputFolders class property
		if isempty(obj.inputFolders)
			if iscell(destFolders)
				nAddFolders = length(destFolders);
				for addFolderNo = 1:nAddFolders
					obj.inputFolders{addFolderNo,1} = destFolders{addFolderNo};
				end
			else
				nAddFolders = 1;
				obj.inputFolders{1,1} = destFolders;
			end
			obj.dataPath = obj.inputFolders;
		else
			if iscell(destFolders)
				nAddFolders = length(destFolders);
				for addFolderNo = 1:nAddFolders
					obj.inputFolders{end+addFolderNo,1} = destFolders{addFolderNo};
				end
			else
				nAddFolders = 1;
				obj.inputFolders{end+1,1} = destFolders;
			end
			% obj.inputFolders = [obj.inputFolders(:); destFolders(:)];
			% obj.inputFolders = cat(obj.inputFolders,destFolders);
			obj.dataPath = obj.inputFolders;
			% obj.dataPath = cat(obj.dataPath,destFolders);
		end
		obj.modelGetFileInfo();
	catch err
		obj.foldersToAnalyze = [];
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end
function downsampleTiffMovieFxnObj(folderPath,options)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	options.deflateLevel = 1;

	% get list of regular expressions
	fileRegexpListOriginal = getFileList([folderPath filesep],options.srcSubfolderFileFilterRegexp);
	fileRegexpList = regexp(fileRegexpListOriginal, options.srcSubfolderFileFilterRegexp,'match');
	options.srcSubfolderFileFilterRegexp
	cellfun(@display,fileRegexpList);
	% fileRegexpList = getFileList([folderPath filesep], options.srcFolderFilterRegexp);
	nFileRegexp = length(fileRegexpList);

	% go over each set of regexp
	for fileRegexpNo=1:nFileRegexp
		disp('==========================================')
		thisFileRegexp = fileRegexpList{fileRegexpNo};
		thisFileRegexp = regexprep(thisFileRegexp,options.srcSubfolderFileFilterRegexpExt,'');
		thisFileRegexp = thisFileRegexp{1};
		display(thisFileRegexp)
		[~,thisFileRegexp,~] = fileparts(thisFileRegexp);
		% get a list of tiff associated with regexp
		[folderPath filesep]
		[thisFileRegexp '.*.tif$']
		movieList = getFileList([folderPath filesep], [thisFileRegexp '.*.tif$']);
		% due to naming convention, shift last to 1st
		movieList = circshift(movieList,[0 1]);
		cellfun(@display,movieList);
		% get path string
		[pathstr,name,ext] = fileparts(fileRegexpListOriginal{fileRegexpNo});
		newFilename = [options.downsampleSaveFolder filesep 'concat_' name '.h5'];
		display(['save HDF5: ' newFilename])
		% load and downsample each tiff
		nMovies = length(movieList);
		nMovies
		for movieNo=1:nMovies
			display(repmat('+',1,21))
			display(['downsampling ' num2str(movieNo) '/' num2str(nMovies)])
			inputFilePath = movieList{movieNo};
			display(['input: ' inputFilePath]);
			% [pathstr,name,ext] = fileparts(inputFilePath);

			% downsampleFilename = [pathstr '\concat_' name '.h5']
			doptions.downsampleDimension = 'space';
			doptions.downsampleType = 'bilinear';
			doptions.downsampleFactor = options.downsampleFactor;
			[inputMovie] = downsampleMovie(inputFilePath, 'options',doptions);
			% save the movie
			% save or append each tiff to corresponding HDF5 file
			if movieNo==1
				createHdf5File(newFilename, options.outputDatasetName, inputMovie,'deflateLevel',options.deflateLevel);
			else
				appendDataToHdf5(newFilename, options.outputDatasetName, inputMovie);
			end
		end
		srcFilenameTxt = [pathstr filesep name '.txt'];
		srcFilenameXml = [pathstr filesep name '.xml'];
		destFilenameTxt = [options.downsampleSaveFolder filesep name '.txt']
		destFilenameXml = [options.downsampleSaveFolder filesep name '.xml']
		if exist(srcFilenameTxt,'file')
			display([srcFilenameTxt ' > ' destFilenameTxt])
			copyfile(srcFilenameTxt,destFilenameTxt)
		elseif exist(srcFilenameXml,'file')
			display([srcFilenameXml ' > ' destFilenameXml])
			copyfile(srcFilenameXml,destFilenameXml)
		end
	end
end
function downsampleHDFMovieFxnObj(movieList,options)
	% downsamples an HDF5 movie, normally the raw recording files

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% display(movieList)
	cellfun(@display,movieList);
	nMovies = length(movieList);
	for movieNo=1:nMovies
		display(repmat('+',1,21))
		display(['downsampling ' num2str(movieNo) '/' num2str(nMovies)])
		inputFilePath = movieList{movieNo};
		display(['input: ' inputFilePath]);
		[pathstr,name,ext] = fileparts(inputFilePath);
		downsampleFilename = [pathstr '\concat_' name '.h5'];
		srcFilenameTxt = [pathstr filesep name '.txt'];
		srcFilenameXml = [pathstr filesep name '.xml'];
		downsampleFilenameAlt = [options.downsampleSaveFolder '\concat_' name '.h5'];

		fprintf('downsampleFilename = %s\nsrcFilenameTxt = %s\nsrcFilenameXml = %s\n',downsampleFilename,srcFilenameTxt,srcFilenameXml)

		% currentDateTimeStr = datestr(now,'yyyymmdd_HHMM','local');
		% mkdir([thisDir filesep 'processing_info'])
		% diarySaveStr = [thisDir filesep 'processing_info' filesep currentDateTimeStr '_preprocess.log'];
		% display(['saving diary: ' diarySaveStr])
		% diary(diarySaveStr);

		try
			% see if file already exists at the destination
			if isempty(options.downsampleSaveFolder)&exist(downsampleFilename,'file')==0
				% check whether to downsample in different location
				if isempty(options.downsampleSaveFolder)
					downsampleHdf5Movie(inputFilePath, 'inputDatasetName', options.datasetName, 'maxChunkSize', options.maxChunkSize,'downsampleFactor',options.downsampleFactor,'outputDatasetName',options.outputDatasetName,'saveFolderTwo',options.downsampleSaveFolderTwo,'downsampleFactorTwo',options.downsampleFactorTwo);
				else
					downsampleHdf5Movie(inputFilePath, 'inputDatasetName', options.datasetName, 'maxChunkSize', options.maxChunkSize,'saveFolder',options.downsampleSaveFolder,'downsampleFactor',options.downsampleFactor,'outputDatasetName',options.outputDatasetName,'saveFolderTwo',options.downsampleSaveFolderTwo,'downsampleFactorTwo',options.downsampleFactorTwo);

					% copy information files to new downsample folder location
					destFilenameTxt = [options.downsampleSaveFolder filesep name '.txt'];
					destFilenameXml = [options.downsampleSaveFolder filesep name '.xml'];
					fprintf('destFilenameTxt = %s\destFilenameXml = %s\n',destFilenameTxt,destFilenameXml)
					if exist(srcFilenameTxt,'file')
						copyfile(srcFilenameTxt,destFilenameTxt)
					elseif exist(srcFilenameXml,'file')
						copyfile(srcFilenameXml,destFilenameXml)
					end
					if ~isempty(options.downsampleSaveFolderTwo)
						destFilenameTxt = [options.downsampleSaveFolderTwo filesep name '.txt']
						destFilenameXml = [options.downsampleSaveFolderTwo filesep name '.xml']
						if exist(srcFilenameTxt,'file')
							copyfile(srcFilenameTxt,destFilenameTxt)
						elseif exist(srcFilenameXml,'file')
							copyfile(srcFilenameXml,destFilenameXml)
						end
					end
				end
			elseif ~isempty(options.downsampleSaveFolder)&exist(downsampleFilenameAlt,'file')==0
				downsampleHdf5Movie(inputFilePath, 'inputDatasetName', options.datasetName, 'maxChunkSize', options.maxChunkSize,'saveFolder',options.downsampleSaveFolder,'downsampleFactor',options.downsampleFactor,'outputDatasetName',options.outputDatasetName,'saveFolderTwo',options.downsampleSaveFolderTwo,'downsampleFactorTwo',options.downsampleFactorTwo);
				destFilenameTxt = [options.downsampleSaveFolder filesep name '.txt']
				destFilenameXml = [options.downsampleSaveFolder filesep name '.xml']
				if exist(srcFilenameTxt,'file')
					copyfile(srcFilenameTxt,destFilenameTxt)
				elseif exist(srcFilenameXml,'file')
					copyfile(srcFilenameXml,destFilenameXml)
				end

				% copy information files to new downsample folder location
				if ~isempty(options.downsampleSaveFolderTwo)
					destFilenameTxt = [options.downsampleSaveFolderTwo filesep name '.txt']
					destFilenameXml = [options.downsampleSaveFolderTwo filesep name '.xml']
					if exist(srcFilenameTxt,'file')
						copyfile(srcFilenameTxt,destFilenameTxt)
					elseif exist(srcFilenameXml,'file')
						copyfile(srcFilenameXml,destFilenameXml)
					end
				end
			else
				display(['skipping: ' inputFilePath])
			end
		catch err
			display(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			display(repmat('@',1,7))
		end
	end
end
function [downsampleSaveFolderMod, downsampleSaveFolderTwoMod] = downsampleIsxdMovieFxnObj(movieList,options)
	% downsamples an HDF5 movie, normally the raw recording files

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% display(movieList)
	cellfun(@display,movieList);
	nMovies = length(movieList);
	for movieNo=1:nMovies
		display(repmat('+',1,21))
		display(['downsampling ' num2str(movieNo) '/' num2str(nMovies)])
		inputFilePath = movieList{movieNo};
		display(['input: ' inputFilePath]);
		[pathstr,name,ext] = fileparts(inputFilePath);
		downsampleFilename = [pathstr '\concat_' name '.h5'];
		srcFilenameTxt = [pathstr filesep name '.txt'];
		srcFilenameXml = [pathstr filesep name '.xml'];
		srcFilenameJson = [pathstr filesep 'session.json'];
		downsampleFilenameAlt = [options.downsampleSaveFolder '\concat_' name '.h5'];

		[~,folderName,~] = fileparts(strip(pathstr,'right',filesep));
		if isempty(options.downsampleSaveFolder)
			downsampleSaveFolderMod = options.downsampleSaveFolder;
		else
			downsampleSaveFolderMod = [options.downsampleSaveFolder filesep folderName];
		end
		if isempty(options.downsampleSaveFolderTwo)
			downsampleSaveFolderTwoMod = options.downsampleSaveFolderTwo;
		else
			downsampleSaveFolderTwoMod = [options.downsampleSaveFolderTwo filesep folderName];
		end
		if ~exist(downsampleSaveFolderMod,'dir');mkdir(downsampleSaveFolderMod);end
		if ~exist(downsampleSaveFolderTwoMod,'dir');mkdir(downsampleSaveFolderTwoMod);end

		if ~isempty(options.downsampleSaveFolder)
			% Copy JSON files
			jsonFileList = getFileList(pathstr,'.*.json');
			if ~isempty(jsonFileList)
				for jsonNo = 1:length(jsonFileList)
					jsonSrcFile = jsonFileList{jsonNo};
					[~,jsonName,ext] = fileparts(jsonSrcFile);
					jsonDestFile = [downsampleSaveFolderMod filesep jsonName ext];
					if exist(jsonSrcFile,'file')
						fprintf('jsonSrcFile = %s\jsonDestFile = %s\n',jsonSrcFile,jsonDestFile)
						copyfile(jsonSrcFile,jsonDestFile)
					end
					if ~isempty(options.downsampleSaveFolderTwo)
						jsonDestFile = [downsampleSaveFolderTwoMod filesep jsonName ext];
						if exist(jsonSrcFile,'file')
							fprintf('jsonSrcFile = %s\jsonDestFile = %s\n',jsonSrcFile,jsonDestFile)
							copyfile(jsonSrcFile,jsonDestFile)
						end
					end
				end
			end
		end
		% fprintf('downsampleFilename = %s\nsrcFilenameTxt = %s\nsrcFilenameXml = %s\n',downsampleFilename,srcFilenameTxt,srcFilenameXml)

		% currentDateTimeStr = datestr(now,'yyyymmdd_HHMM','local');
		% mkdir([thisDir filesep 'processing_info'])
		% diarySaveStr = [thisDir filesep 'processing_info' filesep currentDateTimeStr '_preprocess.log'];
		% display(['saving diary: ' diarySaveStr])
		% diary(diarySaveStr);

		try
			% see if file already exists at the destination
			if isempty(options.downsampleSaveFolder)&exist(downsampleFilename,'file')==0
				% check whether to downsample in different location
				if isempty(options.downsampleSaveFolder)
					convertInscopixIsxdToHdf5(inputFilePath, 'inputDatasetName', options.datasetName, 'maxChunkSize', options.maxChunkSize,'downsampleFactor',options.downsampleFactor,'outputDatasetName',options.outputDatasetName,'saveFolderTwo',downsampleSaveFolderTwoMod,'downsampleFactorTwo',options.downsampleFactorTwo);
				else
					convertInscopixIsxdToHdf5(inputFilePath, 'inputDatasetName', options.datasetName, 'maxChunkSize', options.maxChunkSize,'saveFolder',downsampleSaveFolderMod,'downsampleFactor',options.downsampleFactor,'outputDatasetName',options.outputDatasetName,'saveFolderTwo',downsampleSaveFolderTwoMod,'downsampleFactorTwo',options.downsampleFactorTwo);

					% copy information files to new downsample folder location
					% destFilenameTxt = [options.downsampleSaveFolder filesep name '.txt'];
					% destFilenameXml = [options.downsampleSaveFolder filesep name '.xml'];
					% fprintf('destFilenameTxt = %s\destFilenameXml = %s\n',destFilenameTxt,destFilenameXml)
					% if exist(srcFilenameTxt,'file')
					% 	copyfile(srcFilenameTxt,destFilenameTxt)
					% elseif exist(srcFilenameXml,'file')
					% 	copyfile(srcFilenameXml,destFilenameXml)
					% end
					% if ~isempty(options.downsampleSaveFolderTwo)
					% 	destFilenameTxt = [options.downsampleSaveFolderTwo filesep name '.txt']
					% 	destFilenameXml = [options.downsampleSaveFolderTwo filesep name '.xml']
					% 	if exist(srcFilenameTxt,'file')
					% 		copyfile(srcFilenameTxt,destFilenameTxt)
					% 	elseif exist(srcFilenameXml,'file')
					% 		copyfile(srcFilenameXml,destFilenameXml)
					% 	end
					% end
				end
			elseif ~isempty(options.downsampleSaveFolder)&exist(downsampleFilenameAlt,'file')==0
				convertInscopixIsxdToHdf5(inputFilePath, 'inputDatasetName', options.datasetName, 'maxChunkSize', options.maxChunkSize,'saveFolder',downsampleSaveFolderMod,'downsampleFactor',options.downsampleFactor,'outputDatasetName',options.outputDatasetName,'saveFolderTwo',downsampleSaveFolderTwoMod,'downsampleFactorTwo',options.downsampleFactorTwo);
				% destFilenameTxt = [options.downsampleSaveFolder filesep name '.txt']
				% destFilenameXml = [options.downsampleSaveFolder filesep name '.xml']
				% if exist(srcFilenameTxt,'file')
				% 	copyfile(srcFilenameTxt,destFilenameTxt)
				% elseif exist(srcFilenameXml,'file')
				% 	copyfile(srcFilenameXml,destFilenameXml)
				% end

				% copy information files to new downsample folder location
				% if ~isempty(options.downsampleSaveFolderTwo)
				% 	destFilenameTxt = [options.downsampleSaveFolderTwo filesep name '.txt']
				% 	destFilenameXml = [options.downsampleSaveFolderTwo filesep name '.xml']
				% 	if exist(srcFilenameTxt,'file')
				% 		copyfile(srcFilenameTxt,destFilenameTxt)
				% 	elseif exist(srcFilenameXml,'file')
				% 		copyfile(srcFilenameXml,destFilenameXml)
				% 	end
				% end
			else
				display(['skipping: ' inputFilePath])
			end
		catch err
			display(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			display(repmat('@',1,7))
		end
	end
end

function [newFolderListCell] = subfxnSelectFolders()
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	nExistingFolders = 0;
	usrIdxChoiceStr = {'manually enter folders to list','GUI select folders','Use existing folder list'};
	scnsize = get(0,'ScreenSize');
	[sel, ok] = listdlg('ListString',usrIdxChoiceStr,'ListSize',[scnsize(3)*0.3 scnsize(4)*0.3],'Name','How to add folders with movies to downsample?');
	inputMethod = usrIdxChoiceStr{sel};

	switch inputMethod
		case 'Use existing folder list'
			newFolderListCell = {};
		case 'manually enter folders to list'
			AddOpts.Resize='on';
			AddOpts.WindowStyle='normal';
			AddOpts.Interpreter='tex';
			% inputdlg
			newFolderList = inputdlgcol('One new line per folder path. Enter folder path WITHOUT any single/double quotation marks around the path.','Adding folders to calciumImagingAnalysis object for downsampling.',[21 150],{''},AddOpts);
			if isempty(newFolderList)
				warning('No folders given. Please re-run modelDownsampleRawMovies.')
				return
			end
			newFolderList = newFolderList{1,1};
			% size(newFolderList)
			% class(newFolderList)

			% if obj.guiEnabled==1
			% 	scnsize = get(0,'ScreenSize');
			% 	[fileIdxArray, ok] = listdlg('ListString',obj.fileIDNameArray,'ListSize',[scnsize(3)*0.2 scnsize(4)*0.25],'Name','which folders to analyze?');
			% else
			% 	fileIdxArray = 1:length(obj.fileIDNameArray);
			% end
			nNewFolders = size(newFolderList,1);
			fileIdxArray = (nExistingFolders+1):(nExistingFolders+nNewFolders);
			nFolders = length(fileIdxArray);
			newFolderListCell = {};
			for thisFileNumIdx = 1:nFolders
				% strtrim(newFolderList(thisFileNumIdx,:))
				% class(strtrim(newFolderList(thisFileNumIdx,:)))
				newFolderListCell{thisFileNumIdx} = strtrim(newFolderList(thisFileNumIdx,:));
			end
		case 'GUI select folders'
			thisFileNumIdx = 1;
			newFolderListCell = {};
			pathToAdd = '';
			while ~isempty(pathToAdd)|length(newFolderListCell)<1
				try
					if ischar(pathToAdd)
						disp('Select a folder to add for downsampling. Press cancel to stop adding folders.')
						pathToAdd = uigetdir(pathToAdd,'Select a folder to add for downsampling. Press cancel to stop adding folders.');
						% If user cancels, do not add folder.
						if pathToAdd~=0
							fprintf('Adding folder #%d: %s.\n',thisFileNumIdx,pathToAdd);
							newFolderListCell{thisFileNumIdx} = pathToAdd;
							thisFileNumIdx = thisFileNumIdx+1;
						end
					else
						% Force exit
						pathToAdd = [];
					end
				catch err
					disp(repmat('@',1,7))
					disp(getReport(err,'extended','hyperlinks','on'));
					disp(repmat('@',1,7))
				end
			end
			nNewFolders = length(newFolderListCell);
			fileIdxArray = (nExistingFolders+1):(nExistingFolders+nNewFolders);
		otherwise
			% body
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [w, infos] = sag(problem, in_options)
% Stochastic average descent (SAG) algorithm.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       N. L. Roux, M. Schmidt, and F. R. Bach,
%       "A stochastic gradient method with an exponential convergence rate for finite training sets,"
%       NIPS, 2012.
%
% This file is part of SGDLibrary.
%
% Created by H.Kasai on Feb. 15, 2016
% Modified by H.Kasai on Mar. 25, 2018

    % set dimensions and samples
    d = problem.dim();
    n = problem.samples();

    % set local options
    local_options.sub_mode = 'SAG';

    % merge options
    options = mergeOptions(get_default_options(d), local_options);
    options = mergeOptions(options, in_options);

    % initialize
    total_iter = 0;
    epoch = 0;
    grad_calc_count = 0;
    w = options.w_init;
    num_of_bachces = floor(n / options.batch_size);

    % prepare an array of gradients, and a valiable of average gradient
    grad_array = zeros(d, num_of_bachces);
    grad_ave = mean(grad_array, 2);

    % store first infos
    clear infos;
    [infos, f_val, optgap] = store_infos(problem, w, options, [], epoch, grad_calc_count, 0);

    % display infos
    if options.verbose > 0
        fprintf('%s: Epoch = %03d, cost = %.16e, optgap = %.4e\n', options.sub_mode, epoch, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % permute samples (ToDo)
    perm_idx = 1:n;

    % main loop
    while (optgap > options.tol_optgap) && (epoch < options.max_epoch)

        for j = 1 : num_of_bachces

            % update step-size
            step = options.stepsizefun(total_iter, options);

            % calculate gradient
            start_index = (j-1) * options.batch_size + 1;
            indice_j = perm_idx(start_index:start_index+options.batch_size-1);
            grad = problem.grad(w, indice_j);

            % update average gradient
            if strcmp(options.sub_mode, 'SAG')
                grad_ave = grad_ave + (grad - grad_array(:, j)) / num_of_bachces;
            else % SAGA
                grad_ave = grad_ave + (grad - grad_array(:, j));
            end
            % replace with new grad
            grad_array(:, j) = grad;

            % update w
            w = w - step * grad_ave;

            % proximal operator
            if ismethod(problem, 'prox')
                w = problem.prox(w, step);
            end

            total_iter = total_iter + 1;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + num_of_bachces * options.batch_size;
        epoch = epoch + 1;

        % store infos
        [infos, f_val, optgap] = store_infos(problem, w, options, infos, epoch, grad_calc_count, elapsed_time);

        % display infos
        if options.verbose > 0
            fprintf('%s: Epoch = %03d, cost = %.16e, optgap = %.4e\n', options.sub_mode, epoch, f_val, optgap);
        end
    end

    if optgap < options.tol_optgap
        fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', options.tol_optgap);
    elseif epoch == options.max_epoch
        fprintf('Max epoch reached: max_epochr = %g\n', options.max_epoch);
    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
% The COBRAToolbox: Lp9.m
%
% Purpose:
%     - test if bug fix in LP9 affects the results for test model (without
%     coupling constraints)
%
% Authors:
%     Agnieszka Wegrzyn 2019/04/09, fix compatibility for models with coupling constraints
%

global CBTDIR

% define the features required to run the test
requireOneSolverOf = {'gurobi','ibm_cplex'};

% require the specified toolboxes and solvers, along with a UNIX OS
solversPkgs = prepareTest('requireOneSolverOf', requireOneSolverOf, 'excludeSolvers', {'matlab', 'lp_solve','pdco'});

% save the current path and initialize the test
currentDir = cd(fileparts(which(mfilename)));

% determine the test path for references
testPath = pwd;

% load the model
model = getDistributedModel('ecoli_core_model.mat'); %For all models in the test/models folder and subfolders

% load reference data and input variables
load('testData_LP9.mat');

% load
for k = 1:length(solversPkgs.LP)
    fprintf(' -- Running testLP9.m using the solver interface: %s ... ', solversPkgs.LP{k});

    solverLPOK = changeCobraSolver(solversPkgs.LP{k}, 'LP', 0);

    if solverLPOK
        % created test data for new version
        V = LP9(options.K, options.P, model, options.LPproblem, options.epsilon);
        solTest = model.c'*V;
    end
    assert(isequal(solOri,solTest), '\nResults are not consistent between old and new version of LP9\n')
    % output a success message
    fprintf('\nDone.\n');
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
classdef RWMOP16 < PROBLEM
% <multi> <real> <constrained>
% Cantilever beam design problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- 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
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 2;
            obj.lower    = [0.01 0.20];
            obj.upper    = [0.05 1];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1);
            x2 = x(:,2);
            P  = 1;
            E  = 207000000;
            Sy = 300000;
            delta_max = 0.005;
            rho = 7800;
            % Objectives
            f(:,1) = 0.25 .* rho .* pi .* x2 .* x1.^2;
            f(:,2) = (64 .* P .* x2.^3)./(3 .* E .* pi .* x1.^4);
            % Constraints
            g(:,1) = -Sy + (32 .* P .* x2)./(pi .* x1.^3);
            g(:,2) = -delta_max + (64 .* P .* x2.^3)./(3 .* E .* pi .* x1.^4);
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [3.0630528e+00   2.0408763e-03];
        end
    end
end


text:
function out=convert(in,varargin)
%OPTS Set gradient limits of the MR system.
%   out=convert(in,from,to) Convert the numerical data in, given in
%   specificed units 'from' to units specified in 'to'.
%
%   Valid unit strings are:
%    'Hz/m' 'mT/m' 'rad/ms/mm' 'Hz/m/s' 'mT/m/ms' 'T/m/s' 'rad/ms/mm/ms'

persistent parser
validGradUnits={'Hz/m','mT/m','rad/ms/mm'};
validSlewUnits={'Hz/m/s','mT/m/ms','T/m/s','rad/ms/mm/ms'};
validUnits=cat(2,validGradUnits,validSlewUnits);
if isempty(parser)
    parser = inputParser;
    parser.FunctionName = 'convert';
    parser.addRequired('in',@isnumeric);
    parser.addRequired('fromUnit',...
        @(x) any(validatestring(x,validUnits)));
    parser.addOptional('toUnit',[],...
        @(x) any(validatestring(x,validUnits)));
    parser.addParamValue('gamma',42.576e6,@isnumeric); % Hz/T
end
parse(parser,in,varargin{:});
opt = parser.Results;

% Set default output unit if not given
if isempty(opt.toUnit)
    if ismember(opt.fromUnit,validGradUnits)
        opt.toUnit = validGradUnits{1};
    elseif ismember(opt.fromUnit,validSlewUnits)
        opt.toUnit = validSlewUnits{1};
    end
end

% Convert to standard units
switch opt.fromUnit
    % Grad units
    case 'Hz/m'
        standard = in;
    case 'mT/m'
        standard = in*1e-3*opt.gamma;
    case 'rad/ms/mm'
        standard = in*1e6/(2*pi);
    % Slew units
    case 'Hz/m/s'
        standard = in;
    case {'mT/m/ms','T/m/s'}
        standard = in*opt.gamma;
    case 'rad/ms/mm/ms'
        standard = in*1e9/(2*pi);
end

% Convert from standard units
switch opt.toUnit
    % Grad units
    case 'Hz/m'
        out = standard;
    case 'mT/m'
        out = 1e3*standard/opt.gamma;
    case 'rad/ms/mm'
        out = standard*2*pi*1e-6;
    % Slew units
    case 'Hz/m/s'
        out = standard;
    case {'mT/m/ms','T/m/s'}
        out = standard/opt.gamma;
    case 'rad/ms/mm/ms'
        out = standard*2*pi*1e-9;
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq