text:
% OP_CURLV_P: assemble the matrix B = [b(i,j)], b(i,j) = (coeff p_i, curl v_j), exploiting the tensor product structure.
%
%   mat = op_curlv_p_tp (spv, spp, msh, [coeff]);
%   [rows, cols, values] = op_curlv_p (spv, spp, msh, [coeff]);

% INPUT:
%
%   spv:     object that defines the vector-valued space of trial functions (see sp_vector)
%   spp:     object that defines the scalar-valued space of the multiplier (see sp_scalar)
%   msh:     object that defines the domain partition and the quadrature rule (see msh_cartesian)
%   epsilon: function handle to compute some physical coefficient (optional)
%
% OUTPUT:
%
%   mat:    assembled stiffness matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2011, 2016, 2017 Rafael Vazquez
%
%    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/>.

function varargout = op_curlv_p_tp (spv, spp, msh, coeff)

  for idim = 1:msh.ndim
    size2 = size (spp.sp_univ(idim).connectivity);
    for icomp = 1:spv.ncomp_param
      size1 = size (spv.scalar_spaces{icomp}.sp_univ(idim).connectivity);
      if (size1(2) ~= size2(2) || size1(2) ~= msh.nel_dir(idim))
        error ('One of the discrete spaces is not associated to the mesh')
      end
    end
  end

  A = spalloc (spp.ndof, spv.ndof, 3*spv.ndof);

  ndim = numel (msh.qn);

  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    spv_col = sp_evaluate_col (spv, msh_col, 'value', false, 'curl', true);
    spp_col = sp_evaluate_col (spp, msh_col);

    if (nargin == 4)
      for idim = 1:msh.rdim
        x{idim} = reshape (msh_col.geo_map(idim,:,:), msh_col.nqn, msh_col.nel);
      end
      coeffs = coeff (x{:});
    else
      coeffs = ones (msh_col.nqn, msh_col.nel);
    end

    A = A + op_curlv_p (spv_col, spp_col, msh_col, coeffs);
  end

  if (nargout == 1)
    varargout{1} = A;
  elseif (nargout == 3)
    [rows, cols, vals] = find (A);
    varargout{1} = rows;
    varargout{2} = cols;
    varargout{3} = vals;
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
 function centers = lloyd_max_hist(data, centers, MM, tol, max_iter, chat)
%function centers = lloyd_max_hist(data, centers, MM, tol, max_iter, chat)
%|
%| "improved" version of the lloyd-max algorithm for scalar quantizer design
%| that saves computation by histogramming the data first.
%| Before using this, try using highrate_centers() first!
%|
%| in
%|	data	[N 1]	training data
%|	centers	[K 1]	initial guess of centroids (codebook)
%|	MM		# of histogram bins: K < M < N
%| out
%|	centers [K 1]	final centroids
%|
%| Copyright 2004-7-1, Jeff Fessler, University of Michigan

if nargin == 1 && streq(data, 'test'), lloyd_max_hist_test, return, end
if nargin < 2, ir_usage, end
if nargin < 3, MM = 0; end
if nargin < 4, tol = 1e-3; end
if nargin < 5, max_iter = 40; end
if nargin < 6, chat = 0; end

data = data(:);
Norig = length(data);
if MM
	if ~isreal(data)
		if length(MM) == 1
			MM = [MM MM];
		end
		[wt data] = hist_equal([real(data) imag(data)], MM);
		[dr di] = ndgrid(data(:,1), data(:,2));
		data = dr + 1i * di;
	else
		[wt data] = hist(data, MM);
	end
	data = data(:);
	wt = wt(:);
	% eliminate data bins with 0 wt
	data = data(wt ~= 0);
	wt = wt(wt ~= 0);
else
	wt = ones(size(data));
end

K = length(centers);
if length(unique(centers)) ~= K
	error 'initial centers are not unique'
end
redundant = 0;
if K > length(data)
	if K > Norig
		warning(sprintf('#centroids %d > #data %d!?', K, Norig))
	else
		printf('Warn: #centroids %d > #unique(data) %d.', K, length(data))
	end
	redundant = 1;
end

tol = tol * max(abs(data));
iter = 1;
change = inf;
while iter <= max_iter && change > tol
	if ~redundant && (length(unique(centers)) ~= K)
		warning 'centers became not unique'
	end
	index = quant1_index(data, centers);
	old = centers;
	for kk=1:K
		ik = index == kk;
		if sum(ik)
			if ~sum(wt(ik))
				warning 'bug?'
				keyboard
			end
			centers(kk) = sum(data(ik) .* wt(ik)) ./ sum(wt(ik));
		else
			centers(kk) = NaN;
		end
	end

	% assign any unused centers to the data point(s) furthest from centroids
	if redundant
		centers(isnan(centers)) = 0;
	elseif any(isnan(centers))
		printm('fixing unused centers %d', sum(isnan(centers)))
		while any(isnan(centers))
			cgood = col(centers(~isnan(centers)));
			index = quant1_index(data, cgood);
			dhat = cgood(index);
			iworst = imax(abs(data - dhat));
			knan = find(isnan(centers));
			centers(knan(1)) = data(iworst);
		end
	end
%	disp([iter centers])
	change = max(abs(centers - old));
	iter = iter + 1;
end
if iter == max_iter + 1
	warning 'max %d iterations reached'
end
if chat
	printf('%s: %d iterations', mfilename, iter)
end

% quant1_index()
% find index of nearest centroid.
% this version works even for complex data / centroids
%
function index = quant1_index(x, centers)
[dummy index] = min(abs(outer_sum(x, -centers)), [], 2);
%breaks = (centers(2:end) + centers(1:end-1)) / 2;
%index0 = 1 + sum(outer_sum(x, -breaks) > 0, 2);
%minmax(index-index0)

% quant1_rms()
% rms error between data and its quantized version (complex ok)
%
function rms = quant1_rms(x, centers)
index = quant1_index(x, centers);
rms = sqrt(mean(abs(x(:) - col(centers(index))).^2));

% lloyd_max_hist_test
% self test: compare this approach to matlab's lloyds routine
function lloyd_max_hist_test
rng(0)
x = [10*randn(10^4,1); 50 + 15*randn(10^4,1)];
%x = [10*rand(10^4,1); 50 + 15*rand(10^4,1)];

L = 5;
%c0 = 10*linspace(-1,1,L);
pn = jf_protected_names;
c0 = pn.prctile(x, 100*([1:L]-0.5)/L);

[nx cx] = hist(x, 100);

ch = highrate_centers(x, L);

if exist('lloyds') == 2
	tic
	[p1 c1] = lloyds(x, c0);
	t1 = toc;
else
	p1 = nan(L,1);
	c1 = nan(L,1);
	t1 = nan;
end

tic
c2 = lloyd_max_hist(x, c0, 50);
t2 = toc;

o = ones(L,1);
if im
	plot(cx, nx, '-', c0, 80*o, 'yo', ch, 60*o, 'yx', ...
		c1, 40*o, 'ys', c2, 20*o, 'y^')
	legend('hist', 'c0', 'highrate', 'lloyds', 'hist')
end

tic
c3 = lloyd_max_hist(x, c0);
t3 = toc;

if exist('quantiz') == 2
	[dum1 dum2 distor] = quantiz(x, p1, c1);
else
	distor = nan;
end

rms0 = quant1_rms(x, c0);
rms1 = quant1_rms(x, c1);
rms2 = quant1_rms(x, c2);
rms3 = quant1_rms(x, c3);
rmsh = quant1_rms(x, ch);
printf('lloyds matlab time=%g rms=%g distor=%g', t1, rms1, sqrt(distor))
printf('lloyd_max_hist time=%g rms=%g', t2, rms2)
printf('lloyd_max (no hist) time=%g rms=%g', t3, rms3)
printf('rms0=%g rmsh=%g', rms0, rmsh)
if im
	disp([c0, c1, c2, c3])
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function rxnsOut = plotPredictedExchange(predictedFluxes, otherFluxes, extraFluxes, param)
% plot a comparison of experimental and predicted exchange reaction rates
%
% INPUT
% predictedFluxes:   table with the following variables
%           * .rxns:    k x 1 cell array of reaction identifiers
%           * .v:       k x 1 double of predicted mean flux
%           * .lb:      k x 1 double lower bound on reaction flux
%           * .ub:      k x 1 double upper bound on reaction flux
%           * .labels:  k x 1 cell array of labels to display on y axis
%           * .Properties.Description: string describing the data, used in plot legend
%
% OPTIONAL INPUT
% otherFluxes:  table with the following variables
%           * .rxns:    k x 1 cell array of reaction identifiers
%           * .Properties.Description: sring describing the data, used in plot legend
%           * .labels:  k x 1 cell array of labels to display on y axis
%   and
%           * .mean:    k x 1 double of measured mean flux
%           * .SD:      k x 1 double standard deviation of the measured flux
%
%   or
%           * .v:       k x 1 double of predicted mean flux
%           * .lb:      k x 1 double lower bound on reaction flux
%           * .ub:      k x 1 double upper bound on reaction flux
%
%
% param: paramters structure with the following fields
%           * .rxns: subset of reactions to plot
%           * .measuredFluxes  {(1),0} set to zero to not display experimental fluxes
%
%           * .labelType determines which y axes lables to display
%                       'metabolitePlatform' = platform left and metabolite name right
%                       'metabolite'         = metabolite only right  (default)
%                       'platform'           = platform only left
%           * .expOrder {(1),0}
%                       1  = order reactions by value of experimental flux
%                       0  = order reactions by value of predicted flux (predictedFluxes) for each reaction corresponding to an experimental flux
%           * .saveFigures = 1 saves figure as .fig and .png to current directory
%           * .plotTitle: title of the plot

if exist('predictedFluxes','var')
    predDescription = predictedFluxes.Properties.Description;
else
    error('predictedFluxes missing: must provide a predicted solution');
end
if exist('otherFluxes','var') && ~isempty(otherFluxes)
    otherDescription = otherFluxes.Properties.Description;
else
    otherFluxes = [];
    otherDescription = '';
end
if exist('extraFluxes','var') && ~isempty(extraFluxes)
    extraDescription = extraFluxes.Properties.Description;
else
    extraFluxes = [];
    extraDescription = '';
end
if ~exist('param','var')
    param = struct();
end
if ~isfield(param,'labelType')
    param.labelType = 'metabolite';
end
if ~isfield(param,'saveFigures')
    param.saveFigures = 0;
end
if ~isfield(param,'expOrder')
    param.expOrder = 1;
end
if ~isfield(param,'plotTitle')
    param.plotTitle = [];
end
if ~isfield(param,'measuredBounds')
    param.measuredBounds = 1;
end
if ~isfield(param,'predictedBounds')
    param.predictedBounds = 1;
end
if ~isfield(param,'otherBounds')
    param.otherBounds = 1;
end
if ~isfield(param,'extraBounds')
    param.extraBounds = 1;
end
if ~isfield(param,'exchangesOnly')
    param.exchangesOnly = 1;
end
if ~isvar(predictedFluxes,'labels')
    if isvar(predictedFluxes,'metNames')
        predictedFluxes.labels = predictedFluxes.metNames;
    end
end
if isvar(predictedFluxes,'rxns')
    predictedFluxes.rxns = predictedFluxes.rxns;
end
if ~isempty(otherFluxes)
    if isvar(otherFluxes,'rxns')
        otherFluxes.rxns = otherFluxes.rxns;
    end
    if ~isvar(otherFluxes,'rxns')
        error('otherFluxes.rxns missing - needed to map to other variables')
    end
end
if ~isempty(extraFluxes)
    if isvar(extraFluxes,'rxns')
        extraFluxes.rxns = extraFluxes.rxns;
    end
    if ~isvar(extraFluxes,'rxns')
        error('extraFluxes.rxns missing - needed to map to other variables')
    end
end

%restrict to a subset of fluxes
if isfield(param,'rxns')
    if isempty(param.rxns)
        fprintf('%s\n','plotPredictedExchange: param.rxns is empty, nothing to plot')
        return
    else
        %hack
        predictedFluxes = predictedFluxes(ismember(predictedFluxes.rxns,param.rxns) & ~ismember(predictedFluxes.labels,'Facilitated Diffusion'),:);
    end
end

% for i=1:length(predictedFluxes.labels)
%     labels2 = strrep(predictedFluxes.labels{i},' for ', ' of ');
%     predictedFluxes.labels{i} = [predictedFluxes.labels{i} ', ' labels2];
% end

%order the predicted exchanges
predictedFluxes.zerobool = (predictedFluxes.v-1e-5)<=0 & (predictedFluxes.v+1e-5)>=0;
predictedFluxes.posbool = sign(predictedFluxes.v)==1 & ~predictedFluxes.zerobool;
predictedFluxes.negbool = sign(predictedFluxes.v)==-1 & ~predictedFluxes.zerobool;
predictedFluxes.uptakeProbability = normcdf(0,predictedFluxes.v,abs(predictedFluxes.v)/10);

toRankOrder = zeros(length(predictedFluxes.v),1);
toRankOrder(predictedFluxes.zerobool) = predictedFluxes.uptakeProbability(predictedFluxes.zerobool);
toRankOrder(predictedFluxes.posbool) =  predictedFluxes.v(predictedFluxes.posbool)+1;
toRankOrder(predictedFluxes.negbool)=predictedFluxes.v(predictedFluxes.negbool);

if 1
    [sortedUptakeProbability,xi] = sort(toRankOrder);
else
    [sortedUptakeProbability,xi] = sort(predictedFluxes.uptakeProbability,'descend');
end
predictedFluxes = predictedFluxes(xi,:);
zerobool = predictedFluxes.zerobool;
posbool = predictedFluxes.posbool;
negbool = predictedFluxes.negbool;

%use log modulus transformation
v = logmod(predictedFluxes.v,10);
if isvar(predictedFluxes,'lb') && isvar(predictedFluxes,'ub')
    vl = logmod(predictedFluxes.lb,10);
    vl = v - vl;
    vu = logmod(predictedFluxes.ub,10);
    vu = vu - v;
else
    xl = logmod(measuredFluxes.mean - measuredFluxes.SD,10);
    xl = x - xl;
    xu = logmod(measuredFluxes.mean + measuredFluxes.SD,10);
    xu = xu - x;
end

if isempty(otherFluxes)
    min_v = min(v);
    max_v = max(v);
    bool=~isnan(v);
else
    %place the other fluxes in the same order as the predictions
    otherFluxes = mapAontoB(otherFluxes.rxns,predictedFluxes.rxns, otherFluxes);

    if isvar(otherFluxes,'v') && isvar(otherFluxes,'lb') && isvar(otherFluxes,'ub')
        %use log modulus transformation
        v2 = logmod(otherFluxes.v,10);
        v2l = logmod(otherFluxes.lb,10);
        v2l = v2 - v2l;
        v2u = logmod(otherFluxes.ub,10);
        v2u = v2u - v2;
    elseif isvar(otherFluxes,'mean') && isvar(otherFluxes,'SD')
        %use log modulus transformation
        v2 = logmod(otherFluxes.mean,10);
        v2l = logmod(otherFluxes.mean - otherFluxes.SD,10);
        v2l = v2 - v2l;
        v2u = logmod(otherFluxes.mean + otherFluxes.SD,10);
        v2u = v2u - v2;
    else
        error('Unrecognised input: otherFluxes')
    end

    %limits
    min_v = min([v;v2]);
    max_v = max([v;v2]);

    bool=~isnan(v) & ~isnan(v2);
end

if isempty(extraFluxes)
    if isempty(otherFluxes)
        min_v = min(v);
        max_v = max(v);
        bool=~isnan(v);
    else
        min_v = min([v;v2]);
        max_v = max([v;v2]);
        bool=~isnan(v) & ~isnan(v2);
    end
else
    %place the other fluxes in the same order as the predictions
    extraFluxes = mapAontoB(extraFluxes.rxns,predictedFluxes.rxns, extraFluxes);

    if isvar(extraFluxes,'v') && isvar(extraFluxes,'lb') && isvar(extraFluxes,'ub')
        %use log modulus transformation
        v3 = logmod(extraFluxes.v,10);
        v3l = logmod(extraFluxes.lb,10);
        v3l = v3 - v3l;
        v3u = logmod(extraFluxes.ub,10);
        v3u = v3u - v3;
    elseif isvar(extraFluxes,'mean') && isvar(extraFluxes,'SD')
        %use log modulus transformation
        v3 = logmod(extraFluxes.mean,10);
        v3l = logmod(extraFluxes.mean - extraFluxes.SD,10);
        v3l = v3 - v3l;
        v3u = logmod(extraFluxes.mean + extraFluxes.SD,10);
        v3u = v3u - v3;
    else
        error('Unrecognised input: extraFluxes')
    end

    min_v = min([v;v2;v3]);
    max_v = max([v;v2;v3]);

    bool=~isnan(v) & ~isnan(v2) & ~isnan(v3);
end

%figure('Renderer', 'painters', 'Position', [10 10 900 1500])
%figure('units','normalized','outerposition',[0 0 1 1])
figure('units','normalized','Position',[0 0 0.5 1.5])
hold on

%patch in the background
X = [min_v max_v max_v min_v];

Y = [nnz(negbool & bool)+nnz(zerobool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5 length(predictedFluxes.v(bool))+0.5 length(predictedFluxes.v(bool))+0.5];

patch(X,Y,[221/255, 192/255, 237/255],'FaceAlpha',0.5)

if any(zerobool)
    Y = [nnz(negbool & bool)+0.5 nnz(negbool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5];
    patch(X,Y,[192/255, 237/255, 207/255],'FaceAlpha',0.5)
end

Y = [0 0 nnz(negbool & bool)+0.5 nnz(negbool & bool)+0.5];
patch(X,Y,[192/255, 233/255, 237/255],'FaceAlpha',0.9)

if param.predictedBounds
    %predicted fluxes
    errorbar(v(bool),(1:length(v(bool))),vl(bool),vu(bool),'horizontal','dk','LineWidth',1.5,'MarkerSize',10)%,'MarkerFaceColor',[1 1 1])
else
    %predicted fluxes
    plot(v(bool),1:length(v(bool)),'dk','MarkerSize',12)%,'MarkerFaceColor',[1 1 1])
end

%x & y limits
ylim([0,length(predictedFluxes.v(bool))+0.5])

%other fluxes
if ~isempty(otherFluxes)
    if param.otherBounds
        errorbar(v2(bool),(1:length(v2(bool)))-0.2,v2l(bool),v2u(bool),'horizontal','o','LineWidth',1.5,'MarkerSize',10, 'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])%,'MarkerFaceColor',[1 1 1])
    else
        %other fluxes
        plot(v2(bool),(1:length(v2(bool)))-0.2,'o','MarkerSize',10,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
    end

    xlim([min_v max_v])
    %xlim([-0.5 0.5])
else
    xlim([min_v max_v])
    %xlim([-0.5 0.5])
    %title(['Comparison of \color{blue}measured \color{black}and \color{red}predicted \color{black}net exchange reaction rates. (Objective ' solution1 ')'])
end

%extra fluxes
if ~isempty(extraFluxes)
    if param.extraBounds
        errorbar(v3(bool),(1:length(v3(bool)))-0.4,v3l(bool),v3u(bool),'horizontal','x','LineWidth',1.5,'MarkerSize',10, 'MarkerEdgeColor','magenta','MarkerFaceColor',[1 .6 .6])%,'MarkerFaceColor',[1 1 1])
    else
        %extra fluxes
        plot(v3(bool),(1:length(v3(bool)))-0.4,'x','MarkerSize',10,'MarkerEdgeColor','magenta','MarkerFaceColor',[1 .6 .6])
    end

    xlim([min_v max_v])
    %xlim([-0.5 0.5])
else
    xlim([min_v max_v])
    %xlim([-0.5 0.5])
    %title(['Comparison of \color{blue}measured \color{black}and \color{red}predicted \color{black}net exchange reaction rates. (Objective ' solution1 ')'])
end

if isempty(param.plotTitle)
    %title(['Predicted net reaction rates: ' predDescription ' vs ' otherDescription])
else
    title(param.plotTitle);
end

if param.exchangesOnly
    if any(zerobool)
        firstLegendNames = {'secretion','unchanged','uptake'};
    else
        firstLegendNames = {'secretion','uptake'};
    end
else
    if any(zerobool)
        firstLegendNames = {'secretion/forward flux','unchanged/zero','uptake/negative flux'};
    else
        firstLegendNames = {'secretion/forward flux','uptake/negative flux'};
    end
end

if isempty(otherFluxes)
    legend([firstLegendNames, predDescription],'location','northwest');
else
    if isempty(extraFluxes)
        legend([firstLegendNames,predDescription,otherDescription],'location','northwest');
    else
         legend([firstLegendNames,predDescription,otherDescription,extraDescription],'location','northwest');
    end
end
xlabel('Net flux (uMol/gDW/hr) on $\log_{10}$ scale: $\textrm{sign}(x) \times \log_{10}(1+x)$','interpreter','latex');

%labels
ax = gca;

%xtick=[-100:10:100];
%xticklab = cellstr(num2str(sign(xtick).*log10(1+abs(xtick)), '%d'));
%set(gca,'XTick',xtick,'XTickLabel',xticklab,'TickLabelInterpreter','tex')

ax.YTick = 1:length(predictedFluxes.v(bool));
ax.YTickLabel = predictedFluxes.labels(bool);
rxnsOut = predictedFluxes.rxns(bool);
%hack to fix some metabolite lables

ax.YTickLabel = strrep(ax.YTickLabel,'Acetyl Isoleucine (Chemspider Id: 9964364)','Acetyl Isoleucine');
ax.YTickLabel = strrep(ax.YTickLabel,'3, 4-dihydroxy','3,4-dihydroxy');
ax.YTickLabel = strrep(ax.YTickLabel,'N-Acetyl-Tyrosine','N-Acetyl-L-tyrosine');
ax.YTickLabel = strrep(ax.YTickLabel,'Oxugen','Oxygen');
ax.YTickLabel = strrep(ax.YTickLabel,'Excahnge of ','');
ax.YTickLabel = strrep(ax.YTickLabel,'CYOR_u10mi','Oxidative phosphorylation (CYOR_u10mi)');
ax.YTickLabel = strrep(ax.YTickLabel,'NADH2_u10mi','NADH Dehydrogenase, Mitochondrial (NADH2_u10m)');

if param.exchangesOnly
    ax.YTickLabel = strrep(ax.YTickLabel,'Exchange of ','');
    ax.YTickLabel = strrep(ax.YTickLabel,' demand','');
end

ax.TickLabelInterpreter='none';
ax.FontSize = 14;
switch param.labelType
    case 'metabolitePlatform'
        ax.YAxisLocation = 'right';
    case 'metabolite'
        ax.YAxisLocation = 'right';
    case 'platform'
        ax.YAxisLocation = 'left';
end

%         %right labels
%         plot(v,1:nnz(predictedExLOCB),'*r')
%         axL = gca;
%         axL.YAxisLocation = 'right';
%         axL.Color = 'none';
%         axL.Box = 'off';
%         yticks(1:length(measuredFluxes.mean))
%         axL.TickLabelInterpreter='latex';
%         axL.YColor ='black';
%         axL.FontSize = 10;
%         yt = yticks;
%         yticklabels(measuredFluxes.labelsR)

grid on
ax = gca;
ax.GridColor = [0 .5 .5];
ax.GridLineStyle = '--';
ax.GridAlpha = 0.5;
ax.Layer = 'top';

if 0
    xlim([-0.25 0.15])
    ylim([15.5 69.5])
    zoomed = 'zoom';
else
    zoomed = '';
end
if isempty(otherFluxes)
    descriptions = predDescription;
else
    descriptions = [predDescription '_vs_' otherDescription];
end

%%%% TODO %%% take the model name from the model structure
%     savefig(['iDN1_measured_vs_predicted_exchanges_' names '_' labelType zoomed])
%     saveas(gcf,['iDN1_measured_vs_predicted_exchanges_' names '_' labelType zoomed],'png')
if param.saveFigures
    savefig(['predicted_exchanges_' names '_' param.labelType zoomed])
    saveas(gcf,['predicted_exchanges_' names '_' param.labelType zoomed],'png')
end

author: opencobra

repo: cobratoolbox


text:
function [x, funVal, ValueL]=LeastC(A, y, z, opts)
%
%%
% Function LeastC
%      Least Squares Loss with the L1 ball constraint (also known as Lasso)
%
%% Problem
%
%  min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2
%  s.t. ||x||_1 <= z
%
%  By default, rsL2=0.
%  When rsL2 is nonzero, this correspons the well-know elastic net.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        Radius of the L1 ball (z >0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         Solution
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 18, 2010.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
%% Related functions
%
%  sll_opts, initFactor, pathSolutionLeast,
%  LeastR, nnLeastR, nnLeastC,
%  eplb
%
%%

%% Verify and initialize the parameters
%%

% Verify the number of input parameters
if (nargin <3)
    error('\n Inputs: A, y and z should be specified!\n');
elseif (nargin==3)
    opts=[];
end

% Get the size of the matrix A
[m,n]=size(A);

% Verify the length of y
if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

% Verify the value of z
if (z<=0)
    error('\n z should be positive!\n');
end

% run sll_opts to set default values (flags)
opts=sll_opts(opts);

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATy=A'*y;
elseif (opts.nFlag==1)
    ATy=A'*y - sum(y) * mu';  ATy=ATy./nu;
else
    invNu=y./nu;              ATy=A'*invNu-sum(invNu)*mu';
end

% L2 norm regularization
if isfield(opts,'rsL2')
    rsL2=opts.rsL2;
    if (rsL2<0)
        error('\n opts.rsL2 should be nonnegative!');
    end
else
    rsL2=0;
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
end

if (opts.init==0) % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail

    x_norm=sum(abs(x)); % L1 norm of x
    x_2norm=x'*x;       % squared two norm of x
    if x_norm>=1e-6
        ratio=initFactor(x_norm, Ax, y, z,'LeastC', rsL2, x_2norm);
        x=ratio*x;    Ax=ratio*Ax;
    end
end

%% The main program

bFlag=0; % this flag tests whether the gradient step only changes a little

%% The Armijo Goldstein line search schemes
if (opts.lFlag==0)

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax.
    % xxp=x - xp
    xp=x; Axp=Ax; xxp=zeros(n,1);

    % alphap and alpha are used for computing the weight in forming search point
    alphap=0; alpha=1;

    % lambda0 is a guess of the root in the Euclidean projection
    lambda0=0;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
        else
            invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end

        % obtain the gradient g
        g=ATAs-ATy + rsL2 * s;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and project v onto the L1 ball
            v=s-g/L;

            % projection
            [x, lambda, zf_step]=eplb(v, n, z, lambda0);
            lambda0=lambda;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            Av=Ax -As;
            r_sum=v'*v; l_sum=Av'*Av;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= (L - rsL2) * ||v||_2^2
            if(l_sum <= r_sum * (L-rsL2))
                break;
            else
                L=max(2*L, l_sum/r_sum + rsL2);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=Axy' * Axy/2 + rsL2/2 * x'*x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end % end of .lFlag=0
end

%% Adaptive Line Search

% we set gamma_0 to the L that is appropriate for the starting point
% opts.x0

if (opts.lFlag==1)

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

    lambda0=0;
    % lambda0 is a guess of the root in the Euclidean projection

    gamma=1;
    % we shall set the value of gamma = L,
    % and L is appropriate for the starting point

    xp=x; Axp=Ax;
    % store x and Ax
    xxp=zeros(n,1);
    % the difference of x and xp

    % compute AT Ax
    if (opts.nFlag==0)
        ATAx=A'*Ax;
    elseif (opts.nFlag==1)
        ATAx=A'*Ax - sum(Ax) * mu';  ATAx=ATAx./nu;
    else
        invNu=Ax./nu;                ATAx=A'*invNu-sum(invNu)*mu';
    end

    % We begin the adaptive line search in the following
    %
    % Note that, in the line search, L and beta are changing

    for iterStep=1:opts.maxIter

        ATAxp=ATAx;
        % store ATAx to ATAxp

        if (iterStep~=1)
            % compute AT Ax
            if (opts.nFlag==0)
                ATAx=A'*Ax;
            elseif (opts.nFlag==1)
                ATAx=A'*Ax - sum(Ax) * mu';  ATAx=ATAx./nu;
            else
                invNu=Ax./nu;                ATAx=A'*invNu-sum(invNu)*mu';
            end
        end

        %--------- Line Search for L begins
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp;
                As=Ax + beta* (Ax-Axp);
                ATAs=ATAx + beta* (ATAx-ATAxp);
                % compute the search point s, A * s, and A' * A * s
            else
                alpha= (-1+ sqrt(5)) / 2;
                beta=0; s=x;  As=Ax; ATAs=ATAx;
            end

            g=ATAs-ATy + rsL2 * s;
            % compute the gradient g

            v=s-g/L;
            % a gradient step based on the search point s

            % projection
            [xnew, lambda, zf_step]=eplb(v, n, z, lambda0);
            lambda0=lambda;

            v=xnew-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A xnew
            if (opts.nFlag==0)
                Axnew=A* xnew;
            elseif (opts.nFlag==1)
                invNu=xnew./nu; mu_invNu=mu * invNu;
                Axnew=A*invNu -repmat(mu_invNu, m, 1);
            else
                Axnew=A*xnew-repmat(mu*xnew, m, 1);     Axnew=Axnew./nu;
            end

            Av=Axnew -As;
            r_sum=v'*v; l_sum=Av'*Av + v'*v * rsL2;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2  + rsL2 * ||v||_2^2<= L* ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end
        %--------- Line Search for L ends

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;
        % store values for L

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;   x=xnew;  xxp=x-xp;
        Axp=Ax; Ax=Axnew;
        % update x and Ax with xnew and Axnew

        Axy=Ax-y;
        funVal(iterStep)=Axy' * Axy/2 + rsL2/2 * x'*x;
        % compute function value

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

author: thomaskuestner


text:
% this function is dedicated for crossover operation used by MATLAB ga
% function.
function mutationChildren = SP_mutation(parents ,options,NVARS, ...
    FitnessFcn, state, thisScore,thisPopulation,mutationRate)

mutationChildren = cell(length(parents),1);% Normally zeros(length(parents),NVARS);
for i=1:length(parents)
    parent = thisPopulation{parents(i)}; % Normally thisPopulation(parents(i),:)
    p = ceil(length(parent) * rand(1,2));
    child = parent;
    child(p(1)) = parent(p(2));
    child(p(2)) = parent(p(1));
    mutationChildren{i} = child; % Normally mutationChildren(i,:)
end

author: Sable

repo: mcbench-benchmarks


text:
function [depFiles] = lexdepfun(rootFunction,varargin)
%
% lexdepfun.m--Recursively searches for file dependencies for the specified
% mfile. If the mfile name specified by the user does not include any path
% information, which.m will be used to acquire it.
%
% Only user-defined mfiles are returned as dependencies.
%
% If optional, second argument is '-quiet', then progress is not displayed
% as search happens.
%
% Syntax: depFiles = lexdepfun(rootFunction,<quietFlag>)
%
% e.g.,   depFiles = lexdepfun('myplot')

% Developed in Matlab 7.6.0.324 (R2008a) on GLNX86.
% Kevin Bartlett (kpb@uvic.ca), 2008-06-25 13:45
%-------------------------------------------------------------------------

% Development note. Rewritten largely from scratch in June 2008 to work
% with class-definition directories ("@-directories").

isQuiet = false;

if nargin > 1
    if strcmpi(varargin{1},'-quiet')
        isQuiet = true;
    end % if
end % if

recursionLevel = 0;
alreadySearched = {};

% Get list of paths to search and make first recursive call.
[pathName,rootBaseName,ext] = fileparts(rootFunction);

if isempty(pathName)
    % If no path information, use which.m to get path to root function.
    rootFunction = whichfile(rootFunction);

    if isempty(rootFunction)
        error([mfilename '.m--rootFunction must be on Matlab''s path.']);
    end % if

end % if

% Get list of directories to search.
pathStr = path;

if isunix
    pathCell = textscan(pathStr,'%s','delimiter',':');
else
    pathCell = textscan(pathStr,'%s','delimiter',';');
end % if

pathCell = pathCell{1};

% Exlude the current directory and parent directory indicators ('.' and
% '..', respectively) from the list of directories to search.
matchIndex = strmatch('.',pathCell);
pathCell(matchIndex) = [];
matchIndex = strmatch('..',pathCell);
pathCell(matchIndex) = [];

% Only interested in user-defined functions, so exclude Mathworks
% mfiles.
matPath = matlabroot;
matchIndex = strmatch(matPath,pathCell);
pathCell(matchIndex) = [];

% Get list of base names of all mfiles on path. Only the base names are
% collected because these are what are searched for lexically in the
% file whose dependencies are being determined.
onPathBaseNames = {};

for iPath = 1:length(pathCell)

    thisDir = pathCell{iPath};
    dirOut = dir(fullfile(thisDir,'*.m'));
    numFiles = length(dirOut);
    thisPathBaseNames = cell(1,numFiles);

    for iFile = 1:numFiles
        [dummy,thisBaseName] = fileparts(dirOut(iFile).name);
        thisPathBaseNames{iFile} = thisBaseName;
    end % for

    onPathBaseNames = union(onPathBaseNames,thisPathBaseNames);

    % If there are any "@-directories" (class-definition directories) in
    % this directory, then the class-definition files they contain are also
    % on Matlab's search path, although they do not appear in the output
    % from path.m. For example, if a sub-directory named "@dial" exists in
    % this directory, then the file ".../@dial/dial.m" should be included
    % in onPathBaseNames.
    dirOut = dir(fullfile(thisDir,'@*'));
    numFiles = length(dirOut);
    atDirNames = cell(1,numFiles);

    for iFile = 1:numFiles
        [dummy,thisBaseName] = fileparts(dirOut(iFile).name);
        thisBaseName = strrep(thisBaseName,'@','');
        atDirNames{iFile} = thisBaseName;
    end % for

    onPathBaseNames = union(onPathBaseNames,atDirNames);

end % for

% Make recursive call to search routine and return results.
if isQuiet == false
    progressStr = ['RecursionLevel ' num2str(recursionLevel) ' (' rootBaseName ') '];
    fprintf(1,'%s',progressStr);
end % if

[depFiles,recursionLevel] = recursive_search(rootFunction,recursionLevel,onPathBaseNames,alreadySearched,isQuiet);

if isQuiet == false
    fprintf(1,'\n');
end % if

%-------------------------------------------------------------------------
function [depFiles,recursionLevel,alreadySearched] = recursive_search(rootFunction,recursionLevel,onPathBaseNames,alreadySearched,isQuiet)
%
% recursive_search.m--Recursively searches the specified function
% (fully-qualified filename) for dependencies.
%
% Syntax: [depFiles,recursionLevel,alreadySearched] = recursive_search(rootFunction,recursionLevel,onPathBaseNames,alreadySearched,isQuiet)

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

depFiles = {rootFunction};
candidates = onPathBaseNames;

% If there is a "private" directory under the root function's directory,
% then the mfiles it contains will need to be considered as possible
% dependencies (private directories don't appear in Matlab's path).
[pathName,rootBaseName,ext] = fileparts(rootFunction);
privDirName = fullfile(pathName,'private');
privDirBaseNames = {};

if exist(privDirName,'dir')

    dirOut = dir(fullfile(privDirName,'*.m'));

    numPrivFiles = length(dirOut);
    privDirBaseNames = cell(1,numPrivFiles);

    for iFile = 1:numPrivFiles
        [dummy,thisBaseName] = fileparts(dirOut(iFile).name);
        privDirBaseNames{iFile} = thisBaseName;
    end % for

    candidates = union(candidates,privDirBaseNames);

end % if

% If the root function lives in an "@-directory" (a class-definition
% directory), then any other files in that directory should be considered
% as possible dependencies, just like the contents of a "private" directory
% (parent directories of @-directories appear in Matlab's path, but the
% @-directories themselves do not).
[dummy,parentDirName] = fileparts(pathName);
classDefBaseNames = {};

if strncmp(parentDirName,'@',1)

    dirOut = dir(fullfile(pathName,'*.m'));

    numClassDefFiles = length(dirOut);
    classDefBaseNames = cell(1,numClassDefFiles);

    for iFile = 1:numClassDefFiles
        [dummy,thisBaseName] = fileparts(dirOut(iFile).name);
        classDefBaseNames{iFile} = thisBaseName;
    end % for

    classDefBaseNames = setdiff(classDefBaseNames,rootBaseName);
    candidates = union(candidates,classDefBaseNames);

end % if

% Look in the mfile being searched for strings that match the candidate
% dependencies.
lexicalMatches = lexical_search(rootFunction,candidates);

% Each of the lexical matches is a possible dependency that itself must be
% searched for dependencies.
numMatches = length(lexicalMatches);

if isQuiet == false
    if numMatches == 0
        progressStr = 'no lexical matches found.';
        fprintf(1,'%s',progressStr);
    end % if
end % if

for iMatch = 1:numMatches

    % Show progress.
    if isQuiet == false
        progressStr = ['lexical match ' num2str(iMatch) ' of ' num2str(numMatches) '.'];
        fprintf(1,'%s',progressStr);
    end % if

    thisMatch = lexicalMatches{iMatch};
    recursiveRootFunction = '';

    if ismember(thisMatch,privDirBaseNames)
        recursiveRootFunction = fullfile(privDirName,[thisMatch '.m']);
    elseif ismember(thisMatch,classDefBaseNames)
        recursiveRootFunction = fullfile(pathName,[thisMatch '.m']);
    elseif ismember(thisMatch,onPathBaseNames)
        % Dependency is on Matlab's path, so can use which.m to find it.
        recursiveRootFunction = whichfile(thisMatch);
    end % if

    if ~isempty(recursiveRootFunction)

        if ~ismember(recursiveRootFunction,alreadySearched)

            recursionLevel = recursionLevel + 1;

            if recursionLevel > 0
                progressStr = ['RecursionLevel ' num2str(recursionLevel) ' (' thisMatch ') '];
                leadingBlank = repmat('.',1,recursionLevel*6);
                progressStr = sprintf('%s%s',leadingBlank,progressStr);
            else
                progressStr = ['RecursionLevel ' num2str(recursionLevel) ' (' thisMatch ') '];
                progressStr = sprintf('\n%s',progressStr);
            end % if

            if isQuiet == false
                fprintf(1,'\n%s',progressStr);
            end % if

            alreadySearched = union(alreadySearched,recursiveRootFunction);

            [newDepFiles,recursionLevel,alreadySearched] = ...
                recursive_search(recursiveRootFunction,recursionLevel,candidates,alreadySearched,isQuiet);

            recursionLevel = recursionLevel - 1;

            if recursionLevel > 0
                progressStr = ['RecursionLevel ' num2str(recursionLevel) ' (' thisMatch ') '];
                leadingBlank = repmat('.',1,recursionLevel*6);
                progressStr = sprintf('%s%s',leadingBlank,progressStr);
            else
                progressStr = ['RecursionLevel ' num2str(recursionLevel) ' (' thisMatch ') '];
                progressStr = sprintf('\n%s',progressStr);
            end % if

            if isQuiet == false
                fprintf(1,'\n%s',progressStr);
            end % if

            depFiles = union(depFiles,newDepFiles);

        end % if

    end % if

end % for

%-------------------------------------------------------------------------
function [lexicalMatches] = lexical_search(filename,candidates)
%
% lexical_search.m--Searches the specified function (fully-qualified
% filename) for strings that match candidates for dependencies. The
% basename of the function itself is not returned in the list of lexical
% matches.
%
% Syntax: lexicalMatches = lexical_search(filename,candidates)

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

% Get list of characters that are not permitted in an mfilename.

% ...The list of permitted characters includes a whitespace. This is so
% that the whitespace separating words in a line of the mfile are not
% removed along with the '(', '}', etc.
goodChars = [ ['a':'z'] ['A':'Z'] ['0':'9'] '_' ' ' ];
badChars = setdiff(char(1:127),goodChars);
tokens = {};

% Read in the mfile and extract string tokens from it. "Tokens" are
% space-delimited strings that conform to our expectations of what an mfile
% basename should look like (no special characters or leading numeric
% characters, etc.)
fullRootMfile = whichfile(filename);

if ~isempty(fullRootMfile)

    fid = fopen(fullRootMfile);

    if fid < 1
        error([mfilename '.m--Could not open mfile ' fullRootMfile '.']);
    end % if

    textLines = textscan(fid,'%s','commentstyle','%','delimiter','\n','whitespace','');
    fclose(fid);
    textLines = strtrim(textLines{1});
    numLines = length(textLines);

    % The 'commentstyle' parameter for textscan isn't smart enough to
    % detect comment lines when preceded by spaces. Now that textLines has
    % been deblanked, seek out comment lines and delete them.
    for lineCount = 1:numLines

        thisLine = textLines{lineCount};

        if ~isempty(thisLine)
            if strcmp(thisLine(1),'%')
                textLines{lineCount} = '';
            end % if
        end % if

    end % for

    % Cycle through the lines and extract the strings from them.
    for lineCount = 1:numLines

        thisLine = textLines{lineCount};

        % Replace any characters that aren't allowed in mfilenames with
        % whitespaces.
        thisLine(ismember(thisLine,badChars)) = ' ';

        if isempty(thisLine)
            thisLineStrs = '';
        else
            thisLineStrs = textscan(thisLine,'%s');
            thisLineStrs = thisLineStrs{1};
        end % if

        % Ignore line if it is empty.
        if ~isempty(thisLineStrs)

            % Ignore line if it starts with 'function'.
            if strcmp(thisLineStrs{1},'function')==0
                [tokens{length(tokens)+1: ...
                    length(tokens)+length(thisLineStrs)}] = ...
                    deal(thisLineStrs{:});
            end % if

        end % if

    end % for

end % if

% Remove duplicate tokens.
tokens = unique(tokens);
numTokens = length(tokens);

% Remove tokens that have leading numeric characters (illegal for
% variable names or mfile names in Matlab). Perform the removal by
% replacing individual elements with blanks (simplifies indexing).
for iToken = 1:numTokens

    thisToken = tokens{iToken};

    if ismember(thisToken(1),'0123456789')
        tokens{iToken} = '';
    end % if

end % for

% Remove known Matlab keywords.
for iToken = 1:numTokens

    thisToken = tokens{iToken};

    if iskeyword(thisToken)
        tokens{iToken} = '';
    end % if

end % for

% Get rid of the blanks introduced by the removal of unwanted elements.
isGoodToken = true(1,numTokens);

for iToken = 1:numTokens

    thisToken = tokens{iToken};

    if isempty(thisToken)
        isGoodToken(iToken) = false;
    end % if

end % for

numTokens = length(isGoodToken(isGoodToken==true));
[tmp{1:numTokens}] = deal(tokens{isGoodToken});
tokens = tmp;

% Match the tokens with the candidate mfile basenames.
lexicalMatches = intersect(tokens,candidates);

% Remove the filename itself from its list of candidate dependencies.
[dummy,baseName] = fileparts(filename);
lexicalMatches = setdiff(lexicalMatches,baseName);

%-------------------------------------------------------------------------
function [whichFile] = whichfile(UNLIKELYFILENAME)
%
% whichfile.m--Like Matlab's which.m, but ignores variable names matching
% the input string.
%
% Matlab's which.m returns the string 'variable' when given an input string
% that matches the name of a variable on the workspace, even if the string
% also matches an m-file on the path. For example,
%   name = which('jump')
% where an m-file named jump.m exists on the Matlab path AND 'jump' is the
% name of a variable results in the variable "name" being assigned the
% value 'variable'.
%
% whichfile.m avoids this counter-intuitive behaviour by ignoring variables
% that match the input string.
%
% Syntax: whichFile = whichfile(fileName)
%
% e.g.,  ls = 99;
%        badOutput = which('ls') % (returns the string 'variable')
%        whichFile = whichfile('ls')

% Developed in Matlab 7.2.0.232 (R2006a) on PCWIN.
% Kevin Bartlett (kpb@uvic.ca), 2008-02-12 23:38
%-------------------------------------------------------------------------

% This program makes use of a variable named "UNLIKELYFILENAME", under the
% assumption that it is extremely unlikely for a user to want to trace a
% file named "UNLIKELYFILENAME.m". Just in case, check to make sure that
% the variable UNLIKELYFILENAME does not contain the string
% 'UNLIKELYFILENAME'.
if strcmp(UNLIKELYFILENAME,'UNLIKELYFILENAME')
    error([mfilename '.m--The string ''UNLIKELYFILENAME'' is the ONLY string that is not allowed! Buy a lottery ticket!']);
end % if

% Because UNLIKELYFILENAME is the only variable in this workspace, can now
% use Matlab's which.m without inadvertently returning the string
% 'variable'.
whichFile = which(UNLIKELYFILENAME);


text:
function genotypeFactor = genotypeGivenAlleleFreqsFactor(alleleFreqs, genotypeVar)
% This function computes the probability of each genotype given the allele
% frequencies in the population.

% Note that we assume that the copies of the gene are independent.  Thus,
% knowing the allele for one copy of the gene does not affect the
% probability of having each allele for the other copy of the gene.  As a
% result, the probability of a genotype is the product of the frequencies
% of its constituent alleles (or twice that product for heterozygous
% genotypes).

% Input:
%   alleleFreqs: An n x 1 vector of the frequencies of the alleles in the
%   population, where n is the number of alleles
%   genotypeVar: The variable number for the genotype (goes in the .var
%   part of the factor)
%
% Output:
%   genotypeFactor: Factor in which the val has the probability of having
%   each genotype (note that this is the FULL CPD with no evidence
%   observed)

% The number of genotypes is (number of alleles choose 2) + number of
% alleles -- need to add number of alleles at the end to account for
% homozygotes

genotypeFactor = struct('var', [], 'card', [], 'val', []);
numAlleles = length(alleleFreqs);

% Each allele has an ID that is the index of its allele frequency in the
% allele frequency list.  Each genotype also has an ID.  We need allele and
% genotype IDs so that we know what genotype and alleles correspond to each
% probability in the .val part of the factor.  For example, the first entry
% in .val corresponds to the probability of having the genotype with
% genotype ID 1, which consists of having two copies of the allele with
% allele ID 1.  There is a mapping from a pair of allele IDs to genotype
% IDs and from genotype IDs to a pair of allele IDs below; we compute this
% mapping using generateAlleleGenotypeMappers(numAlleles). (A genotype
% consists of 2 alleles.)

[allelesToGenotypes, genotypesToAlleles] = generateAlleleGenotypeMappers(numAlleles);

% One or both of these matrices might be useful.
%
%   1.  allelesToGenotypes: n x n matrix that maps pairs of allele IDs to
%   genotype IDs, where n is the number of alleles -- if
%   allelesToGenotypes(i, j) = k, then the genotype with ID k comprises of
%   the alleles with IDs i and j
%
%   2.  genotypesToAlleles: m x 2 matrix of allele IDs, where m is the
%   number of genotypes -- if genotypesToAlleles(k, :) = [i, j], then the
%   genotype with ID k is comprised of the allele with ID i and the allele
%   with ID j

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%INSERT YOUR CODE HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fill in genotypeFactor.var.  This should be a 1-D row vector.
% Fill in genotypeFactor.card.  This should be a 1-D row vector.
m=size(genotypesToAlleles,1);
genotypeFactor.var = [genotypeVar] ;
genotypeFactor.card = [m];
genotypeFactor.val = zeros(1, prod(genotypeFactor.card));
% Replace the zeros in genotypeFactor.val with the correct values.

for i=1:m
	k = genotypesToAlleles(i,:);
	genotypeFactor.val(i) = alleleFreqs(k(1))*alleleFreqs(k(2));
	if k(1)!=k(2)
		genotypeFactor.val(i)*=2;
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function s=phon2sone(p)
%PHON2SONE convert PHON loudness values to SONEs s=(p)
%Inputs:    p is a matrix of phon values
%
%Outputs:   s is a matrix, the same size as p, of sone values
%
% The phon scale measures perceived loudness in dB; at 1 kHz it is identical to dB SPL
% relative to 20e-6 Pa sound pressure. The sone scale is proportional to apparent loudness
% and, by definition, equals 1 at 40 phon. The form of the loudness curve is taken from [1].
% The hearing threshold at 1 kHz for 18 to 25 year olds with normal hearing is taken from [2].
%
% Refs: [1]	J. Lochner and J. Burger. Form of the loudness function in the presence of masking noise.
%           The Journal of the Acoustical Society of America, 33: 1705, 1961.
%       [2]	ISO/TC43. Acoustics  normal equal-loudness-level contours.
%           Standard ISO 226:2003, Aug. 2003.

%      Copyright (C) Mike Brookes 2012-2013
%      Version: $Id: phon2sone.m 3295 2013-08-02 14:03:11Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent a b d
if isempty(a)
    b=log(10)*0.1*0.27; % 0.27 is the exponent from [1] and [2]
    d=exp(b*2.4); % 2.4 dB is teh hearing threshold from [2]
    a=1/(exp(b*40)-d); % scale factor to make p=40 give s=1
end
if nargout>0
    s=a*(exp(b*p)-d);
else
    if nargin<1 || isempty(p)
        pp=linspace(5,90,100)'; % phon values
    else
        pp=p;
    end
    semilogy(pp,phon2sone(pp));
    axisenlarge(-1);
    yticksi;
    xlabel('phon = dB SPL @ 1 kHz');
    ylabel('sone');
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function asksys(bin,f)
disp('========================================');
disp(' HAM DIEU CHE DICH BIEN: ASK');
disp(' VI DU: asksys([0 1 0 1 1 0 1 1 0],3)');
disp('Written by Nguyen Hoang Minh DHCNTPHCM. he..he..');
disp('========================================');

bin=[0 1 0 1 1 0 1 1 1 0];f=4;k=1000;

t=0:2*pi/(k-1):2*pi;

L = length(bin);sig=cos(f*t);
bit1=ones(1,k);bit0=zeros(1,k);
mbit=[];mcw=[];
%==============================================
for n=1:L;
    if bin(n)==0;
       cw=sig;bit=bit0;
    else
       cw=2*sig;bit=bit1;
    end
   mbit=[mbit bit];
   mcw=[mcw cw];
end
ask=mcw;

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

s=length(mcw);mrec=[];

for m=1:k:s
    vm=abs(mcw(m));
   if vm > 1.5
      rec=bit1;
   elseif vm < 1.2
      rec=bit0;
    end
    mrec=[mrec rec];
end
deask=mrec;
subplot(3,1,1);plot(mbit,'r','linewidth',2);axis([0  k*L -0.5 1.5]);grid on;title('Data in');
subplot(3,1,2);plot(ask,'m','linewidth',1.5);axis([0  k*L -2.5 2.5]);grid on;title('ASK modulation');
subplot(3,1,3);plot(deask,'g','linewidth',1.5);axis([0  k*L -0.5 1.5]);grid on;title('ASK demodulation,Data out');


text:
function perimeterEdges=findPerimeter(mesh)
% function perimeterPoints=findPerimeter(mesh,faceIndexList)
% Finds perimeter points.
% Works by looking at the connection matrix (which is a list of edges)
% The perimeter points are those points that are on edges that are used by only one face

%nVerts=length(mesh.uniqueVertices);

% Create a list of edges...
[edgeList1,edgeList2]=find(triu(mesh.connectionMatrix));

edgeList=[edgeList1,edgeList2];

disp ('Number of edges found:');
numUniqueEdges=length(edgeList);

disp (length(edgeList));

sortedEdgeList=sort(edgeList')';
sortedEdgeList=sub2ind([numUniqueEdges,numUniqueEdges],sortedEdgeList(:,1),sortedEdgeList(:,2));

% Now convert the list of faces into three sets of (sorted!) edges
% And convert the x,y coordinates into unique indices
FaceEdges1=[mesh.uniqueFaceIndexList(:,1),mesh.uniqueFaceIndexList(:,2)];
FaceEdges2=[mesh.uniqueFaceIndexList(:,1),mesh.uniqueFaceIndexList(:,3)];
FaceEdges3=[mesh.uniqueFaceIndexList(:,2),mesh.uniqueFaceIndexList(:,3)];

% Sort them so that [a b] and [b a] are seen as the same edge.
sortedFaceEdges1=sort(FaceEdges1')';
sortedFaceEdges2=sort(FaceEdges2')';
sortedFaceEdges3=sort(FaceEdges3')';

% Make them into unique index numbers to help searching and sorting
[FaceEdges1]=sub2ind([numUniqueEdges,numUniqueEdges],sortedFaceEdges1(:,1),sortedFaceEdges1(:,2));
[FaceEdges2]=sub2ind([numUniqueEdges,numUniqueEdges],sortedFaceEdges2(:,1),sortedFaceEdges2(:,2));
[FaceEdges3]=sub2ind([numUniqueEdges,numUniqueEdges],sortedFaceEdges3(:,1),sortedFaceEdges3(:,2));

% concatenate the list of Face edges into one long list of indices
FaceEdges=[FaceEdges1,FaceEdges2,FaceEdges3]';

% Sort them - should produce doublets of most edges 'cos most edges are part of two faces
sortedFaceEdges=sort(FaceEdges(:));

% SortedFaceEdges is a sorted list of all the edges in the face list. If all the edges are
% members of two faces, it'll just be pairs of numbers.
% Lone edges (part of only one face) will be on their own here as well.
% So the list might look like
% 1,1,2,2,3,3,4,4,5,5,6,6,7,8,9,9,10,10....
% Where edge 7 and edge 8 are lone...
% if n(i)==n(i+1) or n(i)==n(i-1) then n(i) is an internal edge
% So in fact we're looking for cases where the above fails...
diffFromUpper=sortedFaceEdges-shift(sortedFaceEdges,[1,0]);
diffFromLower=sortedFaceEdges-shift(sortedFaceEdges,[-1,0]);
loneEdges=sortedFaceEdges(find(diffFromUpper.*diffFromLower));

[a b c]=unique(sortedFaceEdges(:));
% Just a quick check:
% 'a' tells you what unique edges there are in the faces. This should be the same as edgeList

nFound=length(loneEdges);
if (nFound>0)
   disp(nFound);
   disp ('Perimeter edges found');

% In the test cases, the mesh is closed so there's no perimeter, loneEdges should be empty
% But let's pretend we've found some...
[perimeterEdges(:,1),perimeterEdges(:,2)]=ind2sub([numUniqueEdges,numUniqueEdges],loneEdges);
else
   disp ('No Perimeter edges found - mesh is closed');
   perimeterEdges=[];
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%
% Comparison of various dual norms (aka integral probability metrics).

SetAR = @(ar)set(gca, 'PlotBoxAspectRatio', [1 ar 1]);
rep = ['results/dual-norms/'];
[~,~] = mkdir(rep);

% random test
normalize = @(x)x/sum(x(:));
k = @(x,y)-abs(x-y);
nx = 3; ny = 5;
x = randn(nx,1);
y = randn(ny,1);
p = normalize(rand(nx,1));
q = normalize(rand(ny,1));
d = rkhs_norm(k,x,p,x,p);
d = rkhs_norm(k,x,p,y,q);

name = '1dirac-1dirac';
name = '1dirac-2dirac';

% delta_0 <-> delta_t
xmax = 3;
tlist = linspace(-xmax,xmax,257)';

switch name
    case '1dirac-1dirac'
        p = 1; q = 1;
        x = @(t)0;
        y = @(t)t;
    case '1dirac-2dirac'
        p = 1; q = [1;1]/2;
        x = @(t)0;
        y = @(t)[-t/2;t/2];
    otherwise
        error('Unknown');
end

D = []; lgd = {};
% energy distance
k = @(x,y)-abs(x-y);
f = @(t)sqrt( rkhs_norm(k, x(t),p, y(t),q) );
D(:,end+1) = arrayfun(f,tlist);
lgd{end+1} = 'Energy';
% Gaussian RKHS
sigma = .3;
k = @(x,y)exp( -(x-y).^2 / (2*sigma^2) );
f = @(t)sqrt( rkhs_norm(k, x(t),p, y(t),q) );
D(:,end+1) = arrayfun(f,tlist);
lgd{end+1} = 'Gauss';
% W1 distance
D(:,end+1) = abs(tlist);
lgd{end+1} = 'W_1';
% Flat distance
D(:,end+1) = min(abs(tlist),1);
lgd{end+1} = 'Flat';

%% display
clf;
plot(tlist, D, 'LineWidth', 2);
legend(lgd, 'Location', 'NorthWest');
axis tight;
SetAR(1/2);
set(gca, 'FontSize', 20);
saveas(gcf, [rep name '.eps'], 'epsc');

author: optimaltransport


text:
function [x,fval,exitflag,info] = opti_lpsolve(f,A,b,Aeq,beq,lb,ub,int,sos,opts)
%OPTI_LPSOLVE Solve a LP or MILP using LP_SOLVE
%
%   [x,fval,exitflag,info] = opti_lpsolve(f,A,b,Aeq,beq,lb,ub,int,sos) solves
%   the linear program min f'x where A,b are the inequality constraints, Aeq,
%   beq are the equality constraints, lb,ub are the bounds and int is a
%   string of the integer variables in the form 'CIB' (Continuous / Integer
%   Binary). sos is a structure with fields sostype, sosind, soswt for SOS.
%
%   THIS IS A WRAPPER FOR LP_SOLVE USING THE MEX INTERFACE
%   See supplied LGPL License

%   Copyright (C) 2011 Jonathan Currie (IPL)

% Modified Version of lp_solve.m supplied with LP_SOLVE distribution

t = tic;

% Handle missing arguments
if nargin < 10, opts = optiset; end
if nargin < 9, sos = []; end
if nargin < 8, int = []; end
if nargin < 7, ub = []; end
if nargin < 6, lb = []; end
if nargin < 5, beq = []; end
if nargin < 4, Aeq = []; end
if nargin < 3, error('You must supply at least 3 arguments to opti_lpsolve'); end

%Check for linear constraints & augment as required
if(~isempty(A) || ~isempty(Aeq))
    %All <= for A,b
    con_type = ones(1,size(A,1)); %1 <=
    %Augment A, Aeq
    if(~isempty(Aeq))
        A = [A;Aeq]; b = [b;beq];
        con_type = [con_type 3*ones(1,size(Aeq,1))]; %3 ==
    end
    %Warn re sparse
    if(~issparse(A))
        if(strcmpi(opts.warnings,'all'))
            warning('opti:sparse','The A matrix should be sparse, correcting: [sparse(A)]');
        end
        A = sparse(A);
    end
end

%Load LP_SOLVE Problem
n = length(f);
m = size(A,1);
lp = lp_solve('make_lp', m, n);
%Load Objective
lp_solve('set_obj_fn', lp, f);
lp_solve('set_minim', lp); % default is solving minimum lp.
%Load Linear Constraints
if(~isempty(A))
    lp_solve('set_mat', lp, A);
    lp_solve('set_rh_vec', lp, b);
    lp_solve('set_constr_type', lp, con_type);
end
%Load Bounds
if(~isempty(lb)), lp_solve('set_lowbo', lp, lb); end
if(~isempty(ub)), lp_solve('set_upbo', lp, ub); end
%Set Integer Constraints
ind = find(int ~= 'C');
ismip = 0;
if(any(ind))
    for i = 1:length(ind)
      if(int(ind(i)) == 'B')
          lp_solve('set_binary', lp, ind(i), 1);
      elseif(int(ind(i)) == 'I')
          lp_solve('set_int', lp, ind(i), 1);
      end
    end
    ismip = 1;
end
%Set SOS
if(~isempty(sos) && isfield(sos,'type') && ~isempty(sos.type))
    if(length(sos.type) > 1)
        [r,c] = size(sos.type);
        if(~iscell(sos.index) || ~iscell(sos.weight))
            error('When adding multiple SOS they must be supplied as cell arrays!');
        end
        if(r > c)
            tp = str2num(sos.type); %#ok<*ST2NM>
        else
            tp = str2num(sos.type');
        end
        for i = 1:length(t) %not sure if this is actually working?
            lp_solve('add_SOS', lp, ['set' num2str(i)], tp(i), i, sos.index{i}, sos.weight{i});
        end
    else
        lp_solve('add_SOS', lp, 'set1', str2num(sos.type), 1, sos.index, sos.weight);
    end
end
%Set Options
lp_solve('set_verbose', lp, dispLevel(opts.display));
lp_solve('set_timeout', lp, opts.maxtime); %converts internally to long?
lp_solve('set_bb_depthlimit', lp, opts.maxnodes); %not exactly number of nodes..
lp_solve('set_epsint', lp, opts.tolint);
%Hard Coded Options
lp_solve('set_scaling', lp, 7); %curtis-reid scaling
lp_solve('set_presolve', lp, 1+2+4); %presolve rows + cols + elim linear dependent rows
%Solve
exitflag = lp_solve('solve', lp);
%Get Outputs
[fval, x] = lp_solve('get_solution', lp);
iter = lp_solve('get_total_iter', lp);
msg = lp_solve('get_statustext', lp, exitflag);
if(~ismip)
    try
        lam = lp_solve('get_dual_solution',lp);
    catch %#ok<CTCH>
        lam = struct('ineqlin',[],'eqlin',[],'bounds',[]);
    end
end
%Delete Problem
lp_solve('delete_lp', lp);
%Check for deleted vars
if(length(x) ~= length(f))
    optiwarn('LP_SOLVE presolve has removed %d variable(s) from the optimization problem, thus the returned decision variable vector will be shorter than expected. OPTI does not currently insert the missing variables.',abs(length(x) ~= length(f)));
end

%Assign Outputs
info.Iterations = iter;
info.Time = toc(t);
info.Algorithm = 'LP_SOLVE: Simplex';
info.Status = msg;

switch(exitflag)
    case {0,1} %suboptimal as well...?
        exitflag = 1;
    case {2,3} %infeasible
        exitflag = -1;
    case 7 %timeout
        exitflag = 0;
    otherwise %unknown
        exitflag = -2;
end

%Add constant objective term
if(~isempty(fval) && isfield(opts,'objbias') && ~isempty(opts.objbias))
    fval = fval + opts.objbias;
end

%Assign Lambda
if(~ismip)
    eq = con_type == 3;
    info.Lambda = struct('ineqlin',lam(~eq),'eqlin',lam(eq),'bounds',lam(length(b)+1:end));
else
    info.Lambda = [];
end

function  print_level = dispLevel(lev)
%Return LPSOLVE compatible display level
switch(lower(lev))
    case'off'
        print_level = 3;
    case 'iter'
        print_level = 6;
    case 'final'
        print_level = 4;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function plotOptTorques(model, nlp, gait, indices)

    act_joint_idx = find(arrayfun(@(x)~isempty(x.Actuator),model.Joints));

    if nargin < 4
        indices = act_joint_idx;
    else
        if isempty(indices), return; end
    end

    cont_domain_idx = find(cellfun(@(x)isa(x,'ContinuousDynamics'),{nlp.Phase.Plant}));
    joint_names = {model.Joints.Name};

    t = [];
    u_opt = []; u_lb = [];  u_ub = [];

    for j=cont_domain_idx
        t = [t,gait(j).tspan];

        u_opt = [u_opt,gait(j).inputs.u];
        u_lb  = [u_lb,[nlp.Phase(j).OptVarTable.u.LowerBound]];
        u_ub  = [u_ub,[nlp.Phase(j).OptVarTable.u.UpperBound]];
    end

    ax = [];
    for i=1:length(indices)
        idx = indices(i);
        if ~ismember(idx,act_joint_idx)
            continue;
        end
        f = figure;clf;
        set(f, 'WindowStyle', 'docked');
        %         f.Position = [680 558 560 420];
        ax = [ax, axes(f)]; %#ok<LAXES,*AGROW>
        hold on;
        plot(t, u_opt(i,:), 'b');
        plot(t, u_lb(i,:), 'r--');
        plot(t, u_ub(i,:), 'g--');

        title('Torque');
        legend('u', 'lb', 'ub');

        f.Name = [joint_names{idx},'_torque'];
    end

    linkaxes(ax, 'x');

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function pdepetrans
% 1D transport - modelling with extensions for decay and fast sorption
%    using MATLAB pdepe
%
%   $Ekkehard Holzbecher  $Date: 2006/03/16 $
%--------------------------------------------------------------------------
T = 1;                     % maximum time [s]
L = 1;                     % length [m]
D = 1;                     % diffusivity [m*m/s]
v = 1;                     % velocity [m/s]
lambda = 0.0;              % decay constant [1/s]
sorption = 2;              % sorption-model: no sorption (0), linear (1),
                           %                 Freundlich (2), Langmuir (3)
k1 = 0.0004;               % sorption parameter 1 (R=0 for linear isotherm with Kd, else k1=R)
k2 = 0.5;                  % sorption parameter 2 (Kd for linear isotherm with Kd)
rhob = 1300;               % porous medium bulk density [kg/m*m*m]
theta = 0.2;               % porosity [-]
c0 = 0.0;                  % initial concentration [kg/m*m*m]
cin = 1;                   % boundary concentration [kg/m*m*m]

M = 40;                    % number of timesteps
N = 40;                    % number of nodes
%-------------------------- output parameters
gplot = 0;                 % =1: breakthrough curves; =2: profiles
gsurf = 0;                 % surface
gcont = 0;                 % =1: contours; =2: filled contours
ganim = 2;                 % animation of profiles; =1: single line; =2: all lines

t = linspace (T/M,T,M);    % time discretization
x = linspace (0,L,N);      % space discretization

%----------------------execution-------------------------------------------
if sorption == 1 & k1 <=0
    k1 = 1+k2*rhob/theta;
else
    if sorption > 1 k1 = rhob*k1/theta; end
end
options = odeset; if (c0 == 0) c0 = 1.e-20; end
c = pdepe(0,@transfun,@ictransfun,@bctransfun,x,[0 t],options,D,v,lambda,sorption,k1,k2,c0,cin);

%---------------------- graphical output ----------------------------------
switch gplot
    case 1
        plot ([0 t],c)        % breakthrough curves
        xlabel ('time'); ylabel ('concentration');
    case 2
        plot (x,c','--')      % profiles
        xlabel ('space'); ylabel ('concentration');
end
if gsurf                      % surface plot
    figure; surf (x,[0 t],c);
    xlabel ('space'); ylabel ('time'); zlabel('concentration');
end
if gcont figure; end
switch gcont
    case 1
        contour (x,[0 t],c)   % contours
        grid on; xlabel ('space'); ylabel ('time');
    case 2
        contourf(x,[0 t],c)   % filled contours
        colorbar; xlabel ('space'); ylabel ('time');
end
if (ganim)
    [FileName,PathName] = uiputfile('*.mpg');
    figure; if (ganim > 1) hold on; end
    for j = 1:size(c,1)
        axis manual;  plot (x,c(j,:),'r','LineWidth',2);
        YLim = [min(c0,cin) max(c0,cin)];
        legend (['t=' num2str(T*(j-1)/M)]);
        Anim(j) = getframe;
        plot (x,c(j,:),'b','LineWidth',2);
    end
    mpgwrite (Anim,colormap,[PathName '/' FileName]);     % mgwrite not standard MATLAB
    movie (Anim,0);   % play animation
end

%----------------------functions------------------------------
function [c,f,s] = transfun(x,t,u,DuDx,D,v,lambda,sorption,k1,k2,c0,cin)
switch sorption
    case 0
        R = 1;
    case 1
        R = k1;
    case 2
        R = 1+k1*k2*u^(k2-1);
    case 3
        R = 1+k1*k2*u/(k2+u)/(k2+u);
end
c = R;
f = D*DuDx;
s = -v*DuDx -lambda*R*u;
% --------------------------------------------------------------
function u0 = ictransfun(x,D,v,lambda,sorption,k1,k2,c0,cin)
u0 = c0;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = bctransfun(xl,ul,xr,ur,t,D,v,lambda,sorption,k1,k2,c0,cin)
pl = ul-cin;
ql = 0;
pr = 0;
qr = 1;


text:
%CODEGENERATOR.GENFDYN Generate code for forward dynamics
%
% Iqdd = cGen.genfdyn() is a symbolic vector (1xN) of joint inertial
% reaction forces/torques.
%
% Notes::
% - Side effects of execution depends on the cGen flags:
%   - saveresult: the symbolic expressions are saved to
%     disk in the directory specified by cGen.sympath
%   - genmfun: ready to use m-functions are generated and
%     provided via a subclass of SerialLink stored in cGen.robjpath
%   - genslblock: a Simulink block is generated and stored in a
%     robot specific block library cGen.slib in the directory
%     cGen.basepath
%   - genccode: generates C-functions and -headers in the directory
%     specified by the ccodepath property of the CodeGenerator object.
%   - mex: generates robot specific MEX-functions as replacement for the
%     m-functions mentioned above. Access is provided by the SerialLink
%     subclass. The MEX files rely on the C code generated before.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.geninertia, CodeGenerator.genfkine.

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

function [Iqdd] = genfdyn(CGen)
[q,qd] = CGen.rob.gencoords;
tau = CGen.rob.genforces;
nJoints = CGen.rob.n;

CGen.logmsg([datestr(now),'\tLoading required symbolic expressions\n']);

%% Inertia matrix
CGen.logmsg([datestr(now),'\tLoading inertia matrix row by row']);

I = sym(zeros(nJoints));
for kJoints = 1:nJoints
    CGen.logmsg(' %s ',num2str(kJoints));
    symname = ['inertia_row_',num2str(kJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    if ~exist(fname,'file')
        CGen.logmsg(['\n',datestr(now),'\t Symbolics not found, generating...\n']);
        CGen.geninertia;
    end
    tmpstruct = load(fname);
    I(kJoints,:)=tmpstruct.(symname);

end
CGen.logmsg('\t%s\n',' done!');

%% Matrix of centrifugal and Coriolis forces/torques matrix
CGen.logmsg([datestr(now),'\t\tCoriolis matrix by row']);

C = sym(zeros(nJoints));
for kJoints = 1:nJoints
    CGen.logmsg(' %s ',num2str(kJoints));
    symname = ['coriolis_row_',num2str(kJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    if ~exist(fname,'file')
        CGen.logmsg(['\n',datestr(now),'\t Symbolics not found, generating...\n']);
        CGen.gencoriolis;
    end
    tmpstruct = load(fname);
    C(kJoints,:)=tmpstruct.(symname);

end
CGen.logmsg('\t%s\n',' done!');

%% Vector of gravitational load
CGen.logmsg([datestr(now),'\t\tvector of gravitational forces/torques']);
symname = 'gravload';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if ~exist(fname,'file')
    CGen.logmsg(['\n',datestr(now),'\t Symbolics not found, generating...\n']);
    CGen.gengravload;
end
tmpstruct = load(fname);
G = tmpstruct.(symname);

CGen.logmsg('\t%s\n',' done!');

%% Joint friction
CGen.logmsg([datestr(now),'\t\tjoint friction vector']);
symname = 'friction';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if ~exist(fname,'file')
    CGen.logmsg(['\n',datestr(now),'\t Symbolics not found, generating...\n']);
    CGen.genfriction;
end
tmpstruct = load(fname);
F = tmpstruct.(symname);

CGen.logmsg('\t%s\n',' done!');

% Full inverse dynamics
CGen.logmsg([datestr(now),'\tGenerating symbolic inertial reaction forces/torques expression\n']);
Iqdd = tau.'-C*qd.' -G.' +F.';
Iqdd = Iqdd.';

%% Save symbolic expressions
if CGen.saveresult
    CGen.logmsg([datestr(now),'\tSaving symbolic inertial reaction forces/torques expression']);

    CGen.savesym(Iqdd,'Iqdd','Iqdd.mat')

    CGen.logmsg('\t%s\n',' done!');
end

%% M-Functions
if CGen.genmfun
    CGen.genmfunfdyn;
end

%% Embedded Matlab Function Simulink blocks
if CGen.genslblock
    CGen.genslblockfdyn;
end

%% C-Code
if CGen.genccode
    CGen.genccodefdyn;
end

%% MEX
if CGen.genmex
    CGen.genmexfdyn;
end

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function iNewFiles = import_video(iStudy, VideoFiles)
% IMPORT_VIDEO Link video files to the database.
%
% USAGE:  iNewFiles = import_dipoles(iStudy, VideoFiles=[ask])
%
% INPUT:
%    - iStudy       : Index of the study where to import the DipolesFiles
%    - DipolesFiles : Full filename, or cell array of filenames, of video files to import

% @=============================================================================
% 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, 2015-2019

%% ===== PARSE INPUTS =====
if (nargin < 2)
    VideoFiles = [];
elseif ~iscell(VideoFiles)
    VideoFiles = {VideoFiles};
end
% Returned variables
iNewFiles = [];

%% ===== SELECT FILES =====
% If file to load was not defined : open a dialog box to select it
if isempty(VideoFiles)
    % Get default import directory and formats
    LastUsedDirs   = bst_get('LastUsedDirs');
    % Get Dipoles file
    [VideoFiles, FileFormat] = java_getfile('open', ...
            'Import dipoles...', ...   % Window title
            LastUsedDirs.ImportData, ...    % Last used directory
            'multiple', 'files', ...        % Selection mode
            {{'.avi','.mpg','.mpeg','.mp4','.mp2','.mkv','.wmv','.divx','.mov','.asf'}, 'Video files (*.avi;*.mpg;*.mpeg;*.mp4;*.mp2;*.mkv;*.wmv;*.divx;*.mov;.asf)', 'VIDEO'}, []);
    % If no file was selected: exit
    if isempty(VideoFiles)
        return
    % Make sure it's an array of files
    elseif ~iscell(VideoFiles)
        VideoFiles = {VideoFiles};
    end
    % Save default import directory
    LastUsedDirs.ImportData = bst_fileparts(VideoFiles{1});
    bst_set('LastUsedDirs', LastUsedDirs);
end

%% ===== GET RECORDINGS SYNCHRO =====
RecStart = 0;
% Get study structure
sStudy = bst_get('Study', iStudy);
% Look for the first RAW file in the folder
if ~isempty(sStudy.Data)
    iDataRaw = find(strcmpi({sStudy.Data.DataType}, 'raw'),1);
else
    iDataRaw = [];
end
% Read the syncro information from the RAW link
if ~isempty(iDataRaw)
    % Read raw link
    DataMat = in_bst_data(sStudy.Data(iDataRaw).FileName, 'F');
    sFile = DataMat.F;
    % Get the start time of the recordings
    if isfield(sFile, 'header') && isfield(sFile.header, 'FirstTimeStamp') && ~isempty(sFile.header.FirstTimeStamp)
        RecStart = double(sFile.header.FirstTimeStamp) * 1e-6;
    end
end

%% ===== CREATE LINKS =====
% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');
% Get output study
sStudy = bst_get('Study', iStudy);
% Loop on each file
for iFile = 1:length(VideoFiles)
    % === SYNCHRO WITH SUBTITLES ===
    VideoStart = [];
    % Get file name
    [fPath, fBase, fExt] = bst_fileparts(VideoFiles{iFile});
    % Try to find an .smi file corresponding to the video and use it for the synchronization
    SmiFile = bst_fullfile(fPath, [fBase, '.smi']);
    if file_exist(SmiFile)
        % Open the file
        fid = fopen(SmiFile, 'r');
        if (fid < 0)
            disp(['BST> Error: Cannot open file "' SmiFile '"']);
            continue;
        end
        % Read the file line by line to get the first syncro tag
        while ~feof(fid)
            % Read one line
            newLine = fgetl(fid);
            if ~ischar(newLine)
                break;
            end
            % Look for line: "<SYNC Start=0><P Class=ENUSCC>XXXXX</SYNC>"
            if (length(newLine) > 21) && strcmpi(newLine(1:14), '<SYNC Start=0>')
                VideoStart = str2num(str_striptag(newLine)) .* 1e-6;
                break;
            end
        end
        % Close the file
        fclose(fid);
    end

    % === CREATE STRUCTURE ===
    % Create new empty videolink structure
    VideoLinkMat = struct();
    VideoLinkMat.LinkTo     = VideoFiles{iFile};
    VideoLinkMat.Comment    = fBase;
    VideoLinkMat.VideoStart = VideoStart - RecStart;
    % Add History field
    VideoLinkMat = bst_history('add', VideoLinkMat, 'import', ['Import from: ' VideoFiles{iFile}]);

    % === SAVE NEW FILE ===
    % Create output filename
    OutputFile = bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(sStudy.FileName), ['videolink_', file_standardize(fBase), '.mat']);
    OutputFile = file_unique(OutputFile);
    % Save new file in Brainstorm format
    bst_save(OutputFile, VideoLinkMat, 'v7');

    % === UPDATE DATABASE ===
    % Create structure
    sImage = db_template('image');
    sImage.FileName = file_short(OutputFile);
    sImage.Comment  = VideoLinkMat.Comment;
    % Add to study
    iImage = length(sStudy.Image) + 1;
    sStudy.Image(iImage) = sImage;
    iNewFiles = [iNewFiles, iImage];
end

% Save study
bst_set('Study', iStudy, sStudy);
% Update tree
panel_protocols('UpdateNode', 'Study', iStudy);
% Save database
db_save();

author: brainstorm-tools


text:
function lorenz_cluster_test ( )

%*****************************************************************************80
%
%% LORENZ_CLUSTER_TEST tests the LORENZ_CLUSTER library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LORENZ_CLUSTER_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LORENZ_CLUSTER library.\n' );

  lorenz_cluster_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LORENZ_CLUSTER_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bandwidth pm] = pll_simulation(cap, res, ipump, vco_sensitivity, fout, fcomp )

% Simulates 2nd, 3rd and 4th order PLL loops for the topologies shown
% below using basic control systems theory.  This is useful for design
% verification.
%
% Input
%       - cap is the capacitors of the loop in Farads [C1, C2, C3, C4].  If
%           these components are not present in the loop, set their value
%           to zero.
%       - res is the resistors of the  loop in Ohms [ R2, R3, R4].  If
%           these components are not present in the loop, set their value to zero.
%       - order of the loop, either 2,3, or 4.
%       - ipump is the charge pump current in Amperes
%       - vco_sensitivity is the VCO sensitivity in Hertz/Volt
%       - fout is the output frequency in Hertz
%       - fcomp is the comparison frequency in Hertz
%
% Output
%       - bandwidth is the open loop bandwidth in Hertz
%       - pm is the phase margin in degrees
%
% The methods used here are derived from those presented in Dean Banerjee's
% Book "PLL Performance, Simulation, and Design" 4th Ed available at
% National Semiconductors site www.national.com.  Most of the work is
% derived from Chapter 8 pp.43-47 and Chapter 9 pp. 48-53.
%
%
% Loop Topologies
% 2nd Order
%        +  _____                                   ______
% fcomp -->|Phase|----------------------------------| VCO |---->fout
%          |Det. |     |      |                     |     |  |
%           -----      |      |                      -----   |
%            ^ -      C1     R2                              |
%            |         |      |                              |
%            |        GND    C2                              |
%            |                |                              |
%            |               GND         _______             |
%            ----------------------------| 1/N |--------------
%                                        -------
%
% 3rd Order
%        +  _____                                   ______
% fcomp -->|Phase|-----------------R3---------------| VCO |---->fout
%          |Det. |     |      |         |           |     |  |
%           -----      |      |         |            -----   |
%            ^ -      C1     R2        C3                    |
%            |         |      |         |                    |
%            |        GND    C2        GND                   |
%            |                |                              |
%            |               GND         _______             |
%            ----------------------------| 1/N |--------------
%                                        -------
%
%
% 4th Order
%        +  _____                                   ______
% fcomp -->|Phase|-----------------R3------R4-------| VCO |---->fout
%          |Det. |     |      |         |     |     |     |  |
%           -----      |      |         |     |      -----   |
%            ^ -      C1     R2        C3     C4             |
%            |         |      |         |     |              |
%            |        GND    C2        GND   GND             |
%            |                |                              |
%            |               GND         _______             |
%            ----------------------------| 1/N |--------------
%                                        -------
%
%
%   Author: Ben Gilbert
%   Homepage: http://nicta.com.au/people/gilbertb
%   Email: ben.gilbert (wibble) nicta.com.au
%   (c) 2009 by National ICT Australia (NICTA)
%

%% Setup Parameters
C1 = cap(1);
C2 = cap(2);
C3 = cap(3);
C4 = cap(4);

R2 = res(1);
R3 = res(2);
R4 = res(3);

% Conversion of parameters to more convenient units
Kpd = ipump/2/pi; % phase detector gain
Kvco = vco_sensitivity*2*pi; % vco gain

%% Plot Setup
% Generates logarithmic spaced points for the calculations
fplotstart = 10; % Hz
fplotstop = 10E6; % Hz
plotpoints = 100;
pltfreqs = [];
for ppts=0:plotpoints
    pltfreqs = [ pltfreqs fplotstart * 10 ^ (ppts/plotpoints* log10(fplotstop/fplotstart)) ];
end

%% Loop Poll's Polynomial Coefficients

A0 = C1 + C2 + C3 + C4;
A1 = C2*R2*(C1+C3+C4) + R3*(C1+C2)*(C3+C4) + C4*R4*(C1+C2+C3);
A2 = C1*C2*R2*R3*(C3+C4) + C4*R4*(C2*C3*R3+C1*C3*R3+C1*C2*R2+C2*C3*R2);
A3 = C1*C2*C3*C4*R2*R3*R4;

T2 = R2*C2;

%% Filter Transfer Function
Zfilt = @(s)  (1 + s.*T2)./s./(A3.*s.^3 + A2.*s.^2 + A1.*s + A0); % filter transfer function
zfilt = @(f) Zfilt(2*pi*1i*f); % filter transfer function as a function of frequency

%% VCO Transfer Function
Gvco = @(s) Kvco./s; % vco transfer function
gvco = @(f) Gvco(2*pi*1i*f); % vco transfer function expressed as a function of frequency

%% Forward Path Transfer Function
G = @(s)  Kpd * Gvco(s) .* Zfilt(s); % forward path transfer function
%% Reverse (Feedback) Transfer Function
H = fcomp/fout; % feedback path transfer function
%% Open Loop Transfer Function
GH = @(s) G(s)*H; % open loop transfer function
gh = @(f) GH(2*pi*1i*f); % open loop transfer function expressed as a function of frequency

%% Gain Plots
% Plots of transfer function magnitudes
% figure;
% loglog( ...
%     pltfreqs, (abs(gh(pltfreqs))), ...
%     pltfreqs, (abs(zfilt(pltfreqs))), ...
%     pltfreqs, (abs(gvco(pltfreqs))) ...
% );
% grid on; title('Open Loop Gain and Contributing Factors');
% xlabel('Frequency [Hz]');
% legend('Open Loop Gain','Loop Filter','VCO', 'Location', 'SouthWest');

%% Bode Plot
figure;
    subplot(2,1,1);
        semilogx(pltfreqs, 20*log10(abs(gh(pltfreqs))));
        grid on;
        title('Open Loop Magnitude');
    subplot(2,1,2);
        semilogx(pltfreqs, angle(gh(pltfreqs)).*180/pi);
        grid on;
        title('Open Loop Phase');
        ylim([-180 180]);

%% Find open loop bandwidth and phase margin
ghdB = @(f) 20*log10(abs(gh(f)));
bandwidth = fzero(ghdB, [fplotstart fplotstop]); % find bandwidth numerically
pm = 180 + angle(gh(bandwidth)).*180/pi;

%% Closed Loop Transfer Function
Gcl = @(s) GH(s)./(1 + GH(s)); % closed loop transfer function
gcl = @(f) Gcl(2*pi*1i*f);
%% Closed Loop Plot (Magnitude)
% figure;
%     subplot(2,1,1);
%         semilogx(pltfreqs, 20*log10(abs(gcl(pltfreqs))));
%         grid on;
%         title('Closed Loop Magnitude');
%     subplot(2,1,2);
%         semilogx(pltfreqs, angle(gcl(pltfreqs)).*180/pi);
%         grid on;
%         title('Closed Loop Phase');
%         ylim([-180 180]);


text:
function noise = heavisideNoiseParamInit(noise, y)

% HEAVISIDENOISEPARAMINIT Heaviside classification model's parameter initialisation.

% IVM

if nargin > 1
  nClass1 = sum(y==1);
  nClass2 = sum(y==-1);
  noise.bias = invCumGaussian(nClass1./(nClass2+nClass1));
  noise.numProcess = size(y, 2);
else
  noise.bias = zeros(1, noise.numProcess);
end

noise.eta = 0.01;
noise.nParams = length(noise.bias) + length(noise.eta);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [w, infos] = slbfgs(problem, in_options)
% Stochastic limited-memory quasi-newton methods (Stochastic L-BFGS) algorithms.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       sub_mode:   SQN:
%                   Byrd, R. H., Hansen, S. L., Nocedal, J., & Singer, Y.
%                   "A stochastic quasi-Newton method for large-scale optimization,"
%                   SIAM Journal on Optimization, 26(2), 1008-1031, 2016.
%
%       sub_mode:   SVRG-SQN:
%                   Philipp Moritz, Robert Nishihara, Michael I. Jordan,
%                   "A Linearly-Convergent Stochastic L-BFGS Algorithm,"
%                   Artificial Intelligence and Statistics (AISTATS), 2016.
%
%       sub_mode:   SVRG LBFGS:
%                   R. Kolte, M. Erdogdu and A. Ozgur,
%                   "Accelerating SVRG via second-order information,"
%                   OPT2015, 2015.
%
%
% Created by H.Kasai on Oct. 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 = 'SQN';  % SQN or SVRG-SQN or SVRG-LBFGS
    local_options.mem_size = 20;

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

    % set paramters
    if options.batch_size > n
        options.batch_size = n;
    end

    if ~isfield(in_options, 'batch_hess_size')
        options.batch_hess_size = 20 * options.batch_size;
    end

    if options.batch_hess_size > n
        options.batch_hess_size = n;
    end

    if strcmp(options.sub_mode, 'SQN') || strcmp(options.sub_mode, 'SVRG-SQN')
        if ~isfield(options, 'L')
            options.L = 20;
        else
            options.L = in_options.L;
        end
    elseif strcmp(options.sub_mode, 'SVRG-LBFGS')
        options.L = Inf;
    end

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

    s_array = [];
    y_array = [];
    u_old = w;
    u_new = zeros(d,1);

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

    % set start time
    start_time = tic();

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

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

        % permute samples
        if options.permute_on
            perm_idx = randperm(n);
        else
            perm_idx = 1:n;
        end

        if strcmp(options.sub_mode, 'SVRG-SQN') || strcmp(options.sub_mode, 'SVRG-LBFGS')
            % compute full gradient
            %full_grad_new = problem.grad(w,1:n);
            full_grad_new = problem.full_grad(w);
            % count gradient evaluations
            grad_calc_count = grad_calc_count + n;
        end

        if strcmp(options.sub_mode, 'SVRG-LBFGS')
            if epoch > 0
                % store cavature pair
                s_array = [s_array w - w0];
                y_array = [y_array full_grad_new - full_grad];

                % remove overflowed pair
                if(size(s_array,2)>options.mem_size)
                    s_array(:,1) = [];
                    y_array(:,1) = [];
                end
            end
        end

        if strcmp(options.sub_mode, 'SVRG-SQN') || strcmp(options.sub_mode, 'SVRG-LBFGS')
            % store w for SVRG
            w0 = w;
            full_grad = full_grad_new;
        end

        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);

            % calculate variance reduced gradient
            if strcmp(options.sub_mode, 'SVRG-SQN') || strcmp(options.sub_mode, 'SVRG-LBFGS')
                grad_w0 = problem.grad(w0,indice_j);
                grad = full_grad + grad - grad_w0;
            end

            if epoch > 0
                % perform LBFGS two loop recursion
                Hg = lbfgs_two_loop_recursion(grad, s_array, y_array);
                % update w
                w = w + (step*Hg);
            else
                w = w - (step*grad);
            end

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

            % calculate averaged w
            u_new = u_new + w/options.L;

            % update LBFGS vectors Hessian at every L iteration for 'SQN' or 'SVRG-SQN'
            % 'SVRG-LBFGS' does nothing because of L = Inf
            if(mod(total_iter,options.L)==0 && total_iter)

                % calcluate Hessian-vector product using subsamples
                %sub_indices = datasample((1:n), options.batch_hess_size);
                % "datasample" is supported only in statistics package in Octave.
                % To avoid the packege, the following is an alternative. Modified by H.K. on Mar. 27, 2018.
                perm_sub_idx_hessian = randperm(n);
                sub_indices = perm_sub_idx_hessian(1:options.batch_hess_size);

                % calculate hessian
                %H = problem.hess(w, sub_indices);
                %Hv = H*(u_new - u_old);
                % calculate hessian-vector product
                Hv = problem.hess_vec(w, u_new-u_old, sub_indices);

                % store cavature pair
                % 'y' curvature pair is calculated from a Hessian-vector product.
                s_array = [s_array u_new - u_old];
                y_array = [y_array Hv];

                % remove overflowed pair
                if(size(s_array,2)>options.mem_size)
                    s_array(:,1) = [];
                    y_array(:,1) = [];
                end

                u_old = u_new;
                u_new = zeros(d,1);

                % count gradient evaluations
                grad_calc_count = grad_calc_count + options.batch_hess_size;
            end

            total_iter = total_iter + 1;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + j * 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:
function s = collatz ( seed )

%*****************************************************************************80
%
%% COLLATZ computes the Collatz sequence for a given starting point.
%
%  Discussion:
%
%    The Collatz sequence is defined recursively as follows:
%
%      Let T be the current entry of the sequence, and U the next:
%
%      If T = 1, the sequence terminates (or U = 1, your choice);
%      Else if T is even, U = T / 2;
%      Else (if T is odd, and greater than 1) U = 3 * T + 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, the starting point for the sequence.
%    SEED must be positive.
%
%    Output, integer S[*], the Collatz sequence, starting with SEED.
%    Once the sequence reaches the value 1, no more entries are computed.
%
  if ( seed <= 0 )
    s = [];
    return
  end

  n = 1;
  t = seed;
  s(n) = t;

  while ( t ~= 1 )

    if ( 2 * round ( t / 2 ) == t )
      t = t / 2;
    else
      t = 3 * t + 1;
    end

    n = n + 1;
    s(n) = t;

  end

  return
end


text:
%
sz = [42,31];
res = 3;
j = 2;
min_margin = 2^j*[1,1];

% the smallest multiple of 2^res that is larger than sz by at least 2*min_margin
sz_paded = 2^res*ceil((sz + 2*min_margin)/2^res)

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet