text:
% MATLAB xUnit Test Framework
% Version 3.1 (R2010b) 19-Nov-2010
%
% Running Unit Tests
%   runtests                  - Run unit tests
%
% Writing Unit Tests
%   assertElementsAlmostEqual - Assert floating-point array elements almost equal
%   assertEqual               - Assert that inputs are equal
%   assertFilesEqual          - Assert that two files have the same content
%   assertExceptionThrown     - Assert that specified exception is thrown
%   assertFalse               - Assert that input condition is false
%   assertTrue                - Assert that input condition is true
%   assertVectorsAlmostEqual  - Assert floating-point vectors almost equal in norm sense
%   initTestSuite             - Utility script used for subfunction-based tests
%
% Framework Classes
%   CommandWindowTestRunDisplay - Print test suite results to command window
%   FunctionHandleTestCase    - Test case based on a function handle
%   TestCase                  - Class defining interface for test cases
%   TestCaseInDir             - Test case requiring temporary directory change
%   TestCaseWithAddPath       - Test case requiring temporary path modification
%   TestComponent             - Abstract base class for TestCase and TestSuite
%   TestComponentInDir        - Test component requiring temporary directory change
%   TestLogger                - Collect data (silently) from running test suite
%   TestRunDisplay            - Print test suite execution results
%   TestRunMonitor            - Abstract base class for monitoring test suite
%   TestSuite                 - Collection of TestComponent objects
%   TestSuiteInDir            - Test suite requiring temporary directory change
%   %VerboseTestRunDisplay    - Print test suite execution results

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

author: ricedsp


text:
function x = utl_prune_handles(x)
% Prune unreferenced workspace variables from anonymous functions in the given argument.
%
% In:
%   Data : an arbitrary data structure
%
% Out:
%   Data : the data structure with hidden references in anonymous functions removed
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-11-23

if iscell(x)
    % recurse into cell arrays...
    for k=1:numel(x)
        x{k} = utl_prune_handles(x{k}); end
elseif isstruct(x)
    % recurse into structs...
    for f=fieldnames(x)'
        if ~isempty(x)
            [x.(f{1})] = celldeal(utl_prune_handles({x.(f{1})})); end
    end
elseif isa(x,'function_handle')
    % inspect function handles
    f = char(x);
    if f(1) == '@'
        % anonymous function: take apart...
        parts = functions(x);
        % ... and put back together
        x = make_function(f,parts.workspace{1});
    end
end

% create a function handle
function f = make_function(decl__,workspace__)
% create workspace
for fn__=fieldnames(workspace__)'
    eval([fn__{1} ' = workspace__.(fn__{1}) ;']); end
clear workspace__ fn__;
% evaluate declaration
f = eval(decl__);

repo: REST


text:
function M = getmassmat(node,elem2dof,area,type,K)
%% GETMASSMAT Get mass matrix of the finite element space
%
% M = GETMASSMAT(elem2edge,area,Dlambda,elemType,K) get mass matrix of the finite element
% space specified by elemType.
%
% The type can be:
% - lump: lumped mass matrix for P1 element
% - HB: Hierarchical basis for P2 element
% - NB: Nodal basis for P2 element
% - NBB: Nodal basis with bubble functions for P2 element
%
% All cases for P2 element are added by Lin Zhong.

N = size(node,1);
NT = length(area);
Ndof = double(max(elem2dof(:)));

%% Coefficients
if ~exist('type','var'), type = 'P1'; end
if ~exist('K','var'), K = []; end

%% Assembling
n = size(elem2dof,2);
switch n
    case 3 %% P1 element
    %% Assemble the mass matrix by the mass lumping
    if strcmp(type,'lump')
        M = accumarray([elem2dof(:,1);elem2dof(:,2);elem2dof(:,3)],...
                       [area;area;area]/3,[Ndof,1]);
        if exist('K','var') && ~isempty(K) && ~isnumeric(K) % K is a function
            M = K(node).*M;
        elseif exist('K','var') && ~isempty(K) && isnumeric(K) && size(K,1) == N
            M = K.*M;
        end
        M = spdiags(M,0,N,N);
    else %% Assemble the full mass matrix
        if exist('K','var') && ~isempty(K) && ~isnumeric(K) % K is a function
            center = (node(elem2dof(:,1),:) + node(elem2dof(:,2),:) + ...
                      node(elem2dof(:,3),:))/3;
            area = K(center).*area;
        elseif exist('K','var') && ~isempty(K) && isnumeric(K) && size(K,1) == NT
            area = K.*area;
        end
        M = sparse(N,N);
        for i = 1:3
            for j = 1:3
                   Mij = area*((i==j)+1)/12;
                   M = M + sparse(elem2dof(:,i),elem2dof(:,j),Mij,Ndof,Ndof);
            end
        end
    end
    case 6  %% P2 element
    %% Mass matrices for NBB bases
    if strcmp(type,'NBB')
        elem2dof(:,7) = uint32(Ndof+(1:NT)'); % add bubble index
        Ndof = double(max(elem2dof(:)));
        elemM(1:NT,7)   = 9/20;
        elemM(1:NT,1:3) = 1/20;
        elemM(1:NT,4:6) = 2/15;
        elemM = elemM.*repmat(area,1,7);
        diagM = accumarray(elem2dof(:), elemM(:), [Ndof 1]);
        M = spdiags(double(diagM),0,Ndof,Ndof);
        return;
    end
    %% Mass matrices for HB or NB bases
    [lambda, w] = quadpts(4);
    nQuad = size(lambda,1);
    if strcmp(type,'HB')
        phi(:,1) = lambda(:,1);
        phi(:,2) = lambda(:,2);
        phi(:,3) = lambda(:,3);
    elseif strcmp(type,'NB')
        phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
        phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
        phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
    end
    phi(:,4) = 4*lambda(:,2).*lambda(:,3);
    phi(:,5) = 4*lambda(:,3).*lambda(:,1);
    phi(:,6) = 4*lambda(:,1).*lambda(:,2);

    Nbases = 6; NMhalf = 21;
    ii = zeros(NMhalf*NT,1);
    jj = zeros(NMhalf*NT,1);
    sM = zeros(NMhalf*NT,1);
    index = 0;
    for i = 1:Nbases
        for j = i:Nbases
            Mij = 0;
            for p = 1:nQuad
                Mij = Mij + w(p)*phi(p,i).*phi(p,j);
            end
            Mij = Mij.*area;
            ii(index+1:index+NT) = double(elem2dof(:,i));
            jj(index+1:index+NT) = double(elem2dof(:,j));
            sM(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end
    clear Mij
    diagIdx = (ii == jj);   upperIdx = ~diagIdx;
    M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),Ndof,Ndof);
    MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),Ndof,Ndof);
    M = M + MU + MU';
end

repo: Aero-matlab


text:
function EMfc
%This is an implementation of the EM algorithm for Clustering via Gaussian
%mixture models, using graphical on-line representation.

%Panagiotis Braimakis (s6990029)

%load simulated gaussian data.
clear;clc
%for p=1:1000
load data
%load M5    %uncomment only if you want to take results from k-means
            %algorithm

%if i load the M5 matrix then i get as starting values tha ones given via
%the kmeans algorithm, which as he saw detects only spherical clusters
%(which does not always hold)
%So lets give initial random id's to the observations.

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%OPTIONAL SUB-SAMPLING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5%%%%%%
 %sample from X
 t=10000;
%  Y=randperm(t);
%  X=X(Y,:);

cidx=unidrnd(5,t,1);
%in order to give a SEED for the Z matrix ,which classifies
%each row of data to it's cluster, we use the results of another clustering method
tic;
Z=zeros(t,5);   %Setting up the storing matrix.
for i=1:t
for j=1:5
if (cidx(i)==j)     %cidx was the (t*1) id-matrix of each record (from 1 to 5)
Z(i,j)=1;           %The new (t*5) Z matrix has only ones, where needed
end %if
end %for
end %for

%Since we have the seed we can now implement the M-step of the algorithm &
%get the loop started between the "Maximization" & "Expectation" steps.

    w=0;    %for the graphs.
    c=2;    %starting counter.
    ll(1)=-88888888;  %starting values just to set the first difference between the log-likelihood's
                %in the first step (remember l(-1) & l(0) does not exist!!!)
    ll(2)=-77777777;
while ll(c)-ll(c-1) > 10^(-10);

    c=c+1;
    g=c-2   %iteration number.
    %all these are the estimates for the parameters of mixture models with
    %normal components.(Geoffrey J. Mclachlan &Kaye E. Basford)

    %nk row matrix is an estimate of the number of "points" on each
    %cluster.
    nk=sum(Z,1);

    %just for the plots

    for i=1:t
max(i)=Z(i,1);
k=1;
for j=2:5
if  (Z(i,j)>max(i))
max(i)=Z(i,j);
P(i,k)=0;
k=j;
else
P(i,j)=0;
end
end
P(i,k)=1;
end

    n=sum(P,1);

        %Now the tk row-matrix estimates each time the probabilities of belonging
    %to the jth cluster (j=1:5)
    tk=nk./t;
    %%Now comes the mean's matrix if we suppose that i element of X belongs to the k-th
    %%cluster (that is done via the Z matrix each time).
    sumclusterX=X'*Z;
    diagnk=diag(nk);
    diag1nk=inv(diagnk);
    mean=sumclusterX*diag1nk;

    %Next one is the estimated intra-cluster VAriance-COvariance MAtrix of
    %the data.

 covma=cell(5,1);

    for j=1:5

    Q1=X-ones(t,1)*mean(:,j)';  %x's-means
    Q2=Q1';
A=cell(t,1);
    for i=1:t
        A(i)={Q1(i,:)};
    end
B=cell(t,1);
    for i=1:t
        B(i)={Q2(:,i)};
    end
        T=[Z(:,j) Z(:,j)];

C=cell(1,1);
    for e=1:t
       C{e}=(T(e,:)'.*B{e})*A{e};
    end

    nom=[0 0;0 0];
    for k = 1:length(C)
    nom=C{k}+nom;
    end %for

    covma{j}=nom./nk(j);
    end %for

  %PLOTS the 95% confidence ellipses & the (X,Y) pairs of each cluster in each iteration.

%   w=w+1;
%   subplot(2,2,w);
%   hold on,
%
%
  plot(X(:,1),X(:,2),'m.');
  confreg(mean(:,1),covma{1},n(1),0.05);
  confreg(mean(:,2),covma{2},n(2),0.05);
  confreg(mean(:,3),covma{3},n(3),0.05);
  confreg(mean(:,4),covma{4},n(4),0.05);
  confreg(mean(:,5),covma{5},n(5),0.05);
  axis equal;
  zoom out;
  xlabel(g);
  title('Progress Graph');
  ylabel('Data vs Clusters');
  grid on;
  hold off
pause(0.1);

  if w==4    %tests' whether your initial graph counter is 4 in order to reset him and subplot to
    w=0;    %the correct position.
  else
    w=w;
  end %if

   %f(xi|theta-hat) matrix (n*k) -which is for every xi and every
    %mean,var-cov matrix.

 for k=1:5
     for i=1:t
 Q11=X(i,:)-mean(:,k)';  %xi-mean
 Q22=Q11';
   f(i,k)=(((2*pi)^-1)*(det(covma{k}))^(-1/2))*exp(-(1/2)*(Q11*(covma{k}^-1))*Q22);
 %The l(i,k) matrix is only computed in these two for's in order to avoid extra aggravation
 %of our program.

 l(i,k)=tk(k)*f(i,k);

    end %for
 end %for

%the log-likelihood of our data given the iteration c is equal to:
ll(c)=sum(log(sum(l,2)));

%  E-step (posterior expected value)

for i=1:t
    for k=1:5

        Z(i,k)=(tk(k)*f(i,k))/sum(tk'.*f(i,:)');
    end
end
save sampleresults1 Z mean covma n P ll ;

 end %while

% disp('Time to Find Clusters: ');
% Q(p)=toc
% LL{p}=ll;
% MEANS{p}=mean;
% COVS{p}=covma;
% save MEANS&COVS100 MEANS COVS Q LL;
% p
%end %for p
disp('Have a nice morning')

author: Sable

repo: mcbench-benchmarks


text:
function x = perform_histogram_equalization(x,y,options)

% perform_histogram_equalization - perform histogram equalization
%
%   x = perform_histogram_equalization(x,y,options);
%
%   Change the values of x so that its ordered values match
%   the ordered values of y.
%
%   You can set
%       options.cols=1 to operate along columns
%       options.rows=1 to operate along rows
%       options.absval to operate only on absolute values.
%
%   Copyright (c) 2006 Gabriel Peyr?

options.null = 0;

if isfield(options ,'absval')
    absval = options.absval;
else
    absval = 0;
end
if isfield(options ,'cols')
    cols = options.cols;
else
    cols = 0;
end
if isfield(options ,'rows')
    rows = options.rows;
else
    rows = 0;
end
if isfield(options, 'match_ycbcr')
    match_ycbcr = options.match_ycbcr;
else
    match_ycbcr = 1;
end

% for color images
if size(x,3)>1
    if size(x,3)~=size(y,3)
        error('x and y images must have the same number of components.');
    end

    if size(x,3)==3 && match_ycbcr
        x = rgb2ycbcr(x);
        y = rgb2ycbcr(y);
    end
    % match each color
    for i=1:3
        x(:,:,i) = perform_histogram_equalization(x(:,:,i), y(:,:,i), options);
    end
    if size(x,3)==3 && match_ycbcr
        x = ycbcr2rgb(x);
    end
    return;
end

if cols && rows
    error('You cannote specify both cols and rows');
end
if cols && size(x,2)>1
    if not(size(x,2)==size(y,2))
        error('x and y must have same number of columns');
    end
    for i=1:size(x,2)
        x(:,i) = perform_histogram_equalization(x(:,i),y(:,i),options);
    end
    return;
end
if rows && size(x,1)>1
    if not(size(x,1)==size(y,1))
        error('x and y must have same number of rows');
    end
    for i=1:size(x,1)
        x(i,:) = perform_histogram_equalization(x(i,:),y(i,:),options);
    end
    return;
end

sx = size(x);
x = x(:);
y = y(:);

if absval
    s = sign(x);
    x = abs(x);
    y = abs(y);
end

[vx,Ix] = sort(x);
[vy,Iy] = sort(y);

nx = length(x);
ny = length(y);

ax = linspace(1,ny,nx);
ay = 1:ny;
vx = interp1(ay,vy,ax);
x(Ix) = vx;

if absval
    x = x .* s;
end

x = reshape(x,sx);

function rgb = ycbcr2rgb(in)
%YCBCR2RGB Convert YCbCr values to RGB color space.
%   RGBMAP = YCBCR2RGB(YCBCRMAP) converts the YCbCr values in the
%   colormap YCBCRMAP to the RGB color space. If YCBCRMAP is M-by-3 and
%   contains the YCbCr luminance (Y) and chrominance (Cb and Cr) color
%   values as columns, then RGBMAP is an M-by-3 matrix that contains
%   the red, green, and blue values equivalent to those colors.
%
%   RGB = YCBCR2RGB(YCBCR) converts the YCbCr image to the equivalent
%   truecolor image RGB.
%
%   Class Support
%   -------------
%   If the input is a YCbCr image, it can be of class uint8, uint16,
%   or double; the output image is of the same class as the input
%   image.  If the input is a colormap, the input and output
%   colormaps are both of class double.
%
%   Example
%   -------
%   Convert image from RGB space to YCbCr space and back.
%
%       rgb = imread('board.tif');
%       ycbcr = rgb2ycbcr(rgb);
%       rgb2 = ycbcr2rgb(ycbcr);
%
%   See also NTSC2RGB, RGB2NTSC, RGB2YCBCR.

%   Copyright 1993-2003 The MathWorks, Inc.
%   $Revision: 1.15.4.2 $  $Date: 2003/08/23 05:54:51 $

%   Reference:
%     Charles A. Poynton, "A Technical Introduction to Digital Video",
%     John Wiley & Sons, Inc., 1996

%initialize variables
isColormap = false;
classin = class(in);

%must reshape colormap to be m x n x 3 for transformation
if (ndims(in)==2)
  %colormap
  isColormap=true;
  colors = size(in,1);
  in = reshape(in, [colors 1 3]);
end

%initialize output
rgb = in;

% set up constants for transformation. T alone will transform YCBCR in [0,255]
% to RGB in [0,1]. We must scale T and the offsets to get RGB in the appropriate
% range for uint8 and for uint16 arrays.
T = [65.481 128.553 24.966;...
     -37.797 -74.203 112; ...
     112 -93.786 -18.214];
Tinv = T^-1;
offset = [16;128;128];
Td = 255 * Tinv;
offsetd = Tinv * offset;
T8 = Td;
offset8 = T8 * offset;
T16 = (65535/257) * Tinv;
offset16 = 65535 * Tinv * offset;

switch classin
 case 'double'
  for p = 1:3
    rgb(:,:,p) = Td(p,1) * in(:,:,1) + Td(p,2) * in(:,:,2) + ...
        Td(p,3) * in(:,:,3) - offsetd(p);
  end

 case 'uint8'
  for p = 1:3
    rgb(:,:,p) = imlincomb(T8(p,1),in(:,:,1),T8(p,2),in(:,:,2), ...
                           T8(p,3),in(:,:,3),-offset8(p));
  end

 case 'uint16'
  for p = 1:3
    rgb(:,:,p) = imlincomb(T16(p,1),in(:,:,1),T16(p,2),in(:,:,2), ...
                           T16(p,3),in(:,:,3),-offset16(p));
  end
end

if isColormap
   rgb = reshape(rgb, [colors 3 1]);
end

if isa(rgb,'double')
  rgb = min(max(rgb,0.0),1.0);
end

%%%
%Parse Inputs
%%%
function X = parse_inputs(varargin)

checknargin(1,1,nargin,mfilename);
X = varargin{1};

if ndims(X)==2
  checkinput(X,{'uint8','uint16','double'},'real nonempty',mfilename,'MAP',1);
  if (size(X,2) ~=3 || size(X,1) < 1)
    eid = sprintf('Images:%s:invalidSizeForColormap',mfilename);
    msg = 'MAP must be a m x 3 array.';
    error(eid,'%s',msg);
  end
elseif ndims(X)==3
  checkinput(X,{'uint8','uint16','double'},'real',mfilename,'RGB',1);
  if (size(X,3) ~=3)
    eid = sprintf('Images:%s:invalidTruecolorImage',mfilename);
    msg = 'RGB must a m x n x 3 array.';
    error(eid,'%s',msg);
  end
else
  eid = sprintf('Images:%s:invalidInputSize',mfilename);
  msg = ['RGB2GRAY only accepts two-dimensional or three-dimensional ' ...
         'inputs.'];
  error(eid,'%s',msg);
end

function out = rgb2ycbcr(in)
%RGB2YCBCR Convert RGB values to YCBCR color space.
%   YCBCRMAP = RGB2YCBCR(MAP) converts the RGB values in MAP to
%   the YCBCR color space.  YCBCRMAP is a M-by-3 matrix that contains
%   the YCBCR luminance (Y) and chrominance (Cb and Cr) color values as
%   columns.  Each row represents the equivalent color to the
%   corresponding row in the RGB colormap.
%
%   YCBCR = RGB2YCBCR(RGB) converts the truecolor image RGB to the
%   equivalent image in the YCBCR color space.
%
%   Class Support
%   -------------
%   If the input is an RGB image, it can be of class uint8, uint16,
%   or double; the output image is of the same class as the input
%   image.  If the input is a colormap, the input and output colormaps
%   are both of class double.
%
%   Examples
%   --------
%   Convert RGB image to YCbCr.
%
%      RGB = imread('board.tif');
%      YCBCR = rgb2ycbcr(RGB);
%
%   Convert RGB color space to YCbCr.
%
%      map = jet(256);
%      newmap = rgb2ycbcr(map);

%   See also NTSC2RGB, RGB2NTSC, YCBCR2RGB.

%   Copyright 1993-2003 The MathWorks, Inc.
%   $Revision: 1.13.4.2 $  $Date: 2003/08/23 05:54:37 $

%   Reference:
%     C.A. Poynton, "A Technical Introduction to Digital Video", John Wiley
%     & Sons, Inc., 1996, p. 175

%initialize variables
isColormap = false;
classin = class(in);

%must reshape colormap to be m x n x 3 for transformation
if (ndims(in)==2)
  %colormap
  isColormap=true;
  colors = size(in,1);
  in = reshape(in, [colors 1 3]);
end

% set up constants for transformation
T = [65.481 128.553 24.966;...
     -37.797 -74.203 112; ...
     112 -93.786 -18.214];
offset = [16;128;128];
offset16 = 257 * offset;
fac8 = 1/255;
fac16 = 257/65535;

%initialize output
out = in;

% do transformation
switch classin
 case 'uint8'
  for p=1:3
    out(:,:,p) = imlincomb(T(p,1)*fac8,in(:,:,1),T(p,2)*fac8,in(:,:,2),...
                         T(p,3)*fac8,in(:,:,3),offset(p));
  end
 case 'uint16'
  for p=1:3
    out(:,:,p) = imlincomb(T(p,1)*fac16,in(:,:,1),T(p,2)*fac16,in(:,:,2),...
                         T(p,3)*fac16,in(:,:,3),offset16(p));
  end
 case 'double'
  % These equations transform RGB in [0,1] to YCBCR in [0, 255]
  for p=1:3
    out(:,:,p) =  T(p,1) * in(:,:,1) + T(p,2) * in(:,:,2) + T(p,3) * ...
        in(:,:,3) + offset(p);
  end
  out = out / 255;
end

if isColormap
   out = reshape(out, [colors 3 1]);
end

author: gpeyre

repo: matlab-toolboxes


text:
function [mmHg] = mbar2mmHg(mbar)
% Convert pressure from millibars to millimeters of mercury.
% Chad Greene 2012
mmHg = mbar*0.750062;


text:
function model = glmOmnibus(model, units);
%
% model = glmOmnibus(model, [units]);
%
% Compute the 'omnibus' contrast for a GLM -- an F test of the significance
% of all experimental conditions against baseline -- appending the results
% to a GLM model struct.
%
% model: struct produced by the glm function.
%
% units: 'F' or 'p', specifies the units to return for each voxel. [Default
% 'p', p-value associated w/ the F distribution.]
%
% Appends a model.omnibus field, containing the results of the contrast
% for each voxel, as well as an omnibus_units field, specifying the units
% used.
%
% For details see Burock and Dale, "Estimation and Detection of
% Event-Related fMRI Signals with Temporally Correlated Noise: A
% Statistically Efficient and Unbiased Approach", HBM, 2001.
%
% ras, 01/04/2006.
if ieNotDefined('units'), units = 'p'; end

% Omnibus Significance Test %  (from er_selxavg)
R = eye(model.nh, Navgs_tot);
q = R*hhat;
if model.nh > 1
    qsum = sum(q); % To get a sign %
else
    qsum = q;
end

if (model.nh == 1)
    Fnum = inv(R * model.C * R') * (q.^2) ;  %'
else
    Fnum = sum(q .* (inv(R*model.C*R') * q));  %'
end

Fden = model.nh * (eres_std.^2);
ind = find(Fden == 0);
Fden(ind) = 10^10;
F = sign(qsum) .* Fnum ./ Fden;

if (~isempty(fomnibusstem))
    fname = sprintf('%s_%03d.mat',fomnibusstem,slice);
    tmp = reshape(F,[nrows ncols]);
    er_svtfile(tmp,fname,override);
end

if isequal(lower(units), 'p')
    % convert to p-value, looking it up for the associated F distribution:
    % This will have [J, n-K] degrees of freedom, where J is the number
    % of rows in our restriction matrix R, and n and K are the # of rows
    % and columns in the design matrix X, respectively.
    p = sign(F) .* er_ftest(model.nh, model.dof, abs(F));

    % this prevents values from becoming infinite -- but I'm not
    % sure why p(indz) is set to 1 rather than s
    indz = find(p==0);
    p(indz) = 1;
    p = sign(p).*(-log10(abs(p)));

    omnibus(1:nrows, 1:ncols, slice) = reshape(p, [nrows ncols]);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function n = p13_n ( )

%*****************************************************************************80
%
%% P13_N returns the number of variables for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = - 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [FrontNo,MaxFNo] = NDSort_SDR(PopObj,nSort)
% Do non-dominated sorting by strengthened dominance relation (SDR)

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

    N      = size(PopObj,1);
    NormP  = sum(PopObj,2);
    cosine = 1 - pdist2(PopObj,PopObj,'cosine');
    cosine(logical(eye(length(cosine)))) = 0;
    Angle  = acos(cosine);

    temp  = sort(unique(min(Angle,[],2)));
    minA  = temp(min(ceil(N/2),end));
    Theta = max(1,(Angle./minA).^1);

    dominate = false(N);
    for i = 1 : N-1
        for j = i+1 : N
            if NormP(i)*Theta(i,j) < NormP(j)
                dominate(i,j) = true;
            elseif NormP(j)*Theta(j,i) < NormP(i)
                dominate(j,i) = true;
            end
        end
    end

    FrontNo = inf(1,N);
    MaxFNo  = 0;
    while sum(FrontNo~=inf) < min(nSort,N)
        MaxFNo  = MaxFNo + 1;
        current = ~any(dominate,1) & FrontNo==inf;
        FrontNo(current)    = MaxFNo;
        dominate(current,:) = false;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
target_folder = 'aligned_images';

pdollar_toolbox_path='D:/face project/pdollar-toolbox';
addpath(genpath(pdollar_toolbox_path));

MTCNN_path = 'D:\face project\MTCNN_face_detection_alignment\code\codes\MTCNNv2';
caffe_model_path=[MTCNN_path , '/model'];
addpath(genpath(MTCNN_path));

coord5points = [30.2946, 65.5318, 48.0252, 33.5493, 62.7299; ...
                51.6963, 51.5014, 71.7366, 92.3655, 92.2041];
imgSize = [112, 96];

%caffe.set_mode_cpu();
gpu_id=0;
caffe.set_mode_gpu();
caffe.set_device(gpu_id);
caffe.reset_all();

%three steps's threshold
threshold=[0.6 0.7 0.7]

minsize = 100;

%scale factor
factor=0.709;

%load caffe models
prototxt_dir =strcat(caffe_model_path,'/det1.prototxt');
model_dir = strcat(caffe_model_path,'/det1.caffemodel');
PNet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir = strcat(caffe_model_path,'/det2.prototxt');
model_dir = strcat(caffe_model_path,'/det2.caffemodel');
RNet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir = strcat(caffe_model_path,'/det3.prototxt');
model_dir = strcat(caffe_model_path,'/det3.caffemodel');
ONet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir =  strcat(caffe_model_path,'/det4.prototxt');
model_dir =  strcat(caffe_model_path,'/det4.caffemodel');
LNet=caffe.Net(prototxt_dir,model_dir,'test');
faces=cell(0);

identity_list = dir('./Images');
identity_list = identity_list(3:end);

for i=1:length(identity_list)
    subset_list = dir(fullfile('./Images',identity_list(i).name));
    subset_list = subset_list(3:end);
    for j=1:length(subset_list)
        image_list = dir(fullfile('./Images',identity_list(i).name, subset_list(j).name));
        image_list = image_list(3:end);
        disp([fullfile('./Images',identity_list(i).name, subset_list(j).name) ' ' num2str(length(image_list)) ' files']);
        for k=1:length(image_list)
            filename = fullfile('./Images',identity_list(i).name, subset_list(j).name, image_list(k).name);
            [~, ~, c] = fileparts(filename);
            if strcmp(c,'.png')==0
                continue;
            end;
            img = imread(filename);
            if size(img,1) == 112 && size(img,2) == 96
                continue;
            end;
            [boundingboxes points]=detect_face(img,minsize,PNet,RNet,ONet,LNet,threshold,false,factor);
            if isempty(boundingboxes)
                continue;
            end;
            default_face = 1;
            if size(boundingboxes,1) > 1
                for bb=2:size(boundingboxes,1)
                    if abs((boundingboxes(bb,1) + boundingboxes(bb,3))/2 - size(img,2) / 2) + abs((boundingboxes(bb,2) + boundingboxes(bb,4))/2 - size(img,1) / 2) < ...
                            abs((boundingboxes(default_face,1) + boundingboxes(default_face,3))/2 - size(img,2) / 2) + abs((boundingboxes(default_face,2) + boundingboxes(default_face,4))/2 - size(img,1) / 2)
                        default_face = bb;
                    end;
                end;
            end;
            facial5points = double(reshape(points(:,default_face),[5 2])');
            Tfm =  cp2tform(facial5points', coord5points', 'similarity');
            cropImg = imtransform(img, Tfm, 'XData', [1 imgSize(2)],...
                                          'YData', [1 imgSize(1)], 'Size', imgSize);
            imwrite(cropImg, filename);
            %show detection result
        % 	numbox=size(boundingboxes,1);
        %     figure(1);
        % 	imshow(img)
        % 	hold on;
        % 	for j=1:numbox
        % 		plot(points(1:5,j),points(6:10,j),'g.','MarkerSize',10);
        % 		r=rectangle('Position',[boundingboxes(j,1:2) boundingboxes(j,3:4)-boundingboxes(j,1:2)],'Edgecolor','g','LineWidth',3);
        %     end
        %     hold off;
        %     figure(2);
        %     imshow(cropImg);
        % 	pause

        end;
    end;
end;

author: happynear

repo: FaceVerification


text:
function val = viewGetUI(vw,param,varargin)
% Get data from various view structures
%
% This function is wrapped by viewGet. It should not be called by anything
% else other than viewGet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'), vw = getCurView; end
if notDefined('param'), error('No parameter defined'); end

mrGlobals;
val = [];

switch param

    case 'ishidden'
        val = strcmp(viewGet(vw, 'name'), 'hidden');
    case 'ui'
        if (checkfields(vw, 'ui'))
            val = vw.ui;
        else
            warning('vista:viewError','No user interface found. Returning empty...');
        end
    case 'fignum'
        if (checkfields(vw, 'ui', 'figNum'))
            if isnumeric(vw.ui.figNum), val = vw.ui.figNum;
            else                           val = vw.ui.figNum.Number; end
        else
            warning('vista:viewError','No figure number found. Returning empty...');
        end
    case 'windowhandle'
        if (checkfields(vw, 'ui', 'windowHandle'))
            val = vw.ui.windowHandle;
        else
            warning('vista:viewError','No window handle found. Returning empty...');
        end
    case 'displaymode'
        if (checkfields(vw, 'ui', 'displayMode'))
            val = vw.ui.displayMode;
        else
            warning('vista:viewError', 'No display mode found. Returning empty...');
        end
    case 'anatomymode'
        if (checkfields(vw, 'ui', 'anatMode'))
            val = vw.ui.anatMode;
        else
            warning('vista:viewError', 'No anatomy mode found. Returning empty...');
        end
    case 'coherencemode'
        if (checkfields(vw, 'ui', 'coMode'))
            val = vw.ui.coMode;
        else
            warning('vista:viewError', 'No coherence mode found. Returning empty...');
        end
    case 'correlationmode'
        if (checkfields(vw, 'ui', 'corMode'))
            val = vw.ui.corMode;
        else
            warning('vista:viewError', 'No correlation mode found. Returning empty...');
        end
    case 'phasemode'
        if (checkfields(vw, 'ui', 'phMode'))
            val = vw.ui.phMode;
        else
            warning('vista:viewError', 'No phase mode found. Returning empty...');
        end
    case 'amplitudemode'
        if (checkfields(vw, 'ui', 'ampMode'))
            val = vw.ui.ampMode;
        else
            warning('vista:viewError', 'No amplitude mode found. Returning empty...');
        end
    case 'projectedamplitudemode'
        if (checkfields(vw, 'ui', 'projampMode'))
            val = vw.ui.projampMode;
        else
            warning('vista:viewError', 'No projected amplitude mode found. Returning empty...');
        end
    case 'mapmode'
        if (checkfields(vw, 'ui', 'mapMode'))
            val = vw.ui.mapMode;
        else
            warning('vista:viewError', 'No map mode found. Returning empty...');
        end
    case 'zoom'
        if checkfields(vw, 'ui', 'zoom')
            val = vw.ui.zoom;
        else
            warning('vista:viewError', 'No UI zoom setting found. Returning empty...');
        end
    case 'crosshairs'
        if checkfields(vw, 'ui', 'crosshairs')
            val = vw.ui.crosshairs;
        else
            warning('vista:viewError', 'No crosshairs found. Returning empty...');
        end
    case 'locs'
        if ~isfield(vw, 'loc')
            if checkfields(vw, 'ui', 'sliceNumFields')
                % single-orientation vw: get from slice # fields
                str = (get(vw.ui.sliceNumFields, 'String'))';
                for n=1:3, val(n) = str2num(str{n}); end %#ok<ST2NM>
            else
                warning('vista:viewError', 'No cursor location found. Returning empty...');
            end
        else
            val = vw.loc;
        end

    case 'phasecma'
        % This returns only the color part of the map
        nGrays = viewGet(vw, 'curnumgrays');
        if (isempty(nGrays)),
            warning('vista:viewError', 'Number of grays necessary to retrieve phase color map. Returning empty...');
            return;
        end
        if (checkfields(vw, 'ui', 'phMode', 'cmap'))
            val = round(vw.ui.phMode.cmap(nGrays+1:end,:) * 255)';
        else
            warning('vista:viewError', 'No phase color map found. Returning empty...');
        end
    case 'cmapcurrent'
        nGrays = viewGet(vw, 'curnumgrays');
        displayMode = viewGet(vw, 'displayMode');
        if (isempty(displayMode) || isempty(nGrays)),
            warning('vista:viewError', 'Display mode/number of grays necessary to retrieve color map. Returning empty...');
            return;
        end

        displayMode = [displayMode 'Mode'];
        if (checkfields(vw, 'ui', displayMode, 'cmap'))
            val = round(vw.ui.(displayMode).cmap(nGrays+1:end,:) * 255)';
        else
            warning('vista:viewError', 'No color map found for display mode ''%s''. Returning empty...', displayMode);
        end

    case 'cmapcurmodeclip'
        displayMode = viewGet(vw, 'displayMode');
        if (isempty(displayMode)),
            warning('vista:viewError', 'Display mode necessary to retrieve clip mode. Returning empty...');
            return;
        end

        displayMode = [displayMode 'Mode'];
        if (checkfields(vw, 'ui', displayMode, 'clipMode'))
            val = vw.ui.(displayMode).clipMode;
        else
            warning('vista:viewError', 'No clip mode found for display mode ''%s''. Returning empty...', displayMode);
        end
    case 'cmapcurnumgrays'
        displayMode = viewGet(vw, 'displayMode');
        if (isempty(displayMode)),
            warning('vista:viewError', 'Display mode necessary to retrieve number of grays. Returning empty...');
            return;
        end

        displayMode = [displayMode 'Mode'];
        if (checkfields(vw, 'ui', displayMode, 'numGrays'))
            val = vw.ui.(displayMode).numGrays;
        else
            warning('vista:viewError', 'No number of grays found for display mode ''%s''. Returning empty...', displayMode);
        end
    case 'cmapcurnumcolors'
        displayMode = viewGet(vw, 'displayMode');
        if (isempty(displayMode)),
            warning('vista:viewError', 'Display mode necessary to retrieve number of colors. Returning empty...');
            return;
        end

        displayMode = [displayMode 'Mode'];
        if (checkfields(vw, 'ui', displayMode, 'numColors'))
            val = vw.ui.(displayMode).numColors;
        else
            warning('vista:viewError', 'No number of grays found for display mode ''%s''. Returning empty...', displayMode);
        end
    case 'flipud'
        % Return the boolean indicating whether to invert the image u/d in
        % the graphical user interface.  It is sometimes convenient to do
        % this in the Inplane view if the top of the slice corresponds to
        % the bottom of the brain.
        % Example:
        %   flipud = viewGet(vw, 'flip updown');
        if checkfields(vw, 'ui', 'flipUD'), val = vw.ui.flipUD;
        else                                val = 0; end

    otherwise
        error('Unknown viewGet parameter');

end

return


text:
function mb_iGMM1D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2011 by Mark Bangert markbangert@gmail.com
% inference with an infinite Gaussian mixture model
% http://www.gatsby.ucl.ac.uk/~edward/pub/inf.mix.nips.99.pdf

% this routine may be called without input data.
% a test data set is produced from a mixture of Gaussians

numOfIterations = 1000;
numOfBurnInIter =  300;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% produce data
numOfDataPoints = 1000; % number of data points

% true values of the generating parameters
true_numOfClasses = 10;
true_alpha        = 0.5;
true_lambda       = 180;
true_r            = 0.008;
true_beta         = 3;
true_w            = 100;

% generate weights of the clusters, their means and their variances
true_mus          = randn(1,true_numOfClasses)/true_r + true_lambda;
true_precs        = gamrnd(true_beta,1/true_w,1,true_numOfClasses);
tmp               = gamrnd(true_alpha,1,1,true_numOfClasses);
true_weights      = tmp/sum(tmp); clear tmp;

% generate samples
missingData       = mb_mvrnd(true_weights,numOfDataPoints);
data              = randn(1,numOfDataPoints)./(true_precs(missingData)).^.5 + true_mus(missingData);

        % visualize data and ground truth
        scrsz = get(0,'ScreenSize');

        clusterFigHandle = figure('Position',[1 scrsz(2) scrsz(3)/2 scrsz(4)],'Name','Clustering');

        colorMx = colormap(jet);
        colorMx = colorMx(randperm(64),:);

        hold on
        [counts,bins] = hist(data,floor(numOfDataPoints/10));
        counts = counts/(sum(counts)*(bins(2)-bins(1))); % normalization of histograms
        bar(bins,counts,'EdgeColor','none');
        for i = 1:true_numOfClasses
            y = (true_precs(i)/2/pi)^.5 * exp(-(bins-true_mus(i)).^2*true_precs(i)/2); % normalized gaussian
            plot(bins,y*true_weights(i),'k--','LineWidth',2); % normalize according to weights
        end
        currIgmmVisHandle = [];
        currSolution = 0*y;
        currSolutionVisHandle = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% start inference

tic

dataMean           = mean(data);
dataVar            = var(data);
numOfPointsInClass = numOfDataPoints;

% start the markov chain with one class
alphaProbDensity   = @(x)x^(-3/2)*exp(-1/2/x); % eq 14
alpha              = mb_sliceSampling(alphaProbDensity,rand,1);
lam                = randn*dataVar^.5 + dataMean; % eq 3
r                  = gamrnd(1/2,2/dataVar,1,1); % eq 3
betaProbDensity    = @(x)x^(-3/2)*exp(-1/2/x); % eq 7
beta               = mb_sliceSampling(betaProbDensity,rand,1);
w                  = gamrnd(1/2,2*dataVar,1,1); % eq 7
mus                = randn/r^.5 +lam; % eq 2
precs              = gamrnd(beta/2,2/w/beta,1,1); % eq 6
cs                 = ones(1,numOfDataPoints);
numOfClasses       = 1;

% perform a large number of MCMC sampling steps
iter = 0;
while iter < numOfIterations

    iter = iter + 1;

            % visualization
            figure(clusterFigHandle);
            delete(currIgmmVisHandle);
            currIgmmVisHandle = [];

            for i = 1:numOfClasses

                currIgmmVisHandle = [currIgmmVisHandle plot(bins,(precs(i)/2/pi)^.5*numOfPointsInClass(i)/numOfDataPoints* ...
                    exp(-(bins-mus(i)).^2*precs(i)/2),'Color',colorMx(mod(i-1,size(colorMx,1))+1,:),'LineWidth',2)];

            end

    % collapsed gibbs sampling
    for i = 1:numOfClasses

        numOfDataPointsInCurrClass = numOfPointsInClass(i);
        dataOfCurrClass            = data(cs==i);
        meanOfCurrClass            = mean(dataOfCurrClass);
        precOfCurrClass            = precs(i);

        % for mu - eq 4
        newLoc   = (meanOfCurrClass*numOfDataPointsInCurrClass*precOfCurrClass + lam*r) / (numOfDataPointsInCurrClass*precOfCurrClass + r);
        newScale = 1/(numOfDataPointsInCurrClass*precOfCurrClass + r);
        mus(i)   = randn*newScale^.5 + newLoc;

        % for s - eq 8
        newLoc   = (beta + numOfDataPointsInCurrClass)/2;
        newScale = 2/(w*beta + sum((dataOfCurrClass-mus(i)).^2));
        precs(i) = gamrnd(newLoc,newScale,1,1);

    end

    % for lambda - eq 5
    newLoc   = (dataMean/dataVar + r*sum(mus)) / (1/dataVar + numOfClasses*r);
    newScale = 1 / (1/dataVar + numOfClasses*r);
    lam      = randn*newScale^.5 + newLoc;

    % for r - eq 5
    newLoc   = (1 + numOfClasses)/2;
    newScale = 2/(dataVar + sum((mus-lam).^2));
    r        = gamrnd(newLoc,newScale,1,1);

    % for w - eq 9
    newLoc   = (numOfClasses*beta + 1)/2;
    newScale = 2/(1/dataVar + beta*sum(precs));
    w        = gamrnd(newLoc,newScale,1,1);

    % for beta - eq 9 - stay in log space for a long time
    betaProbDensity = @(x)exp(-numOfClasses*gammaln(x/2) - 1/2/x - .5*(numOfClasses*x-3)*log(x/2) + sum( (x/2)*log(precs*w) - x*w*precs/2 ));
    beta = mb_sliceSampling(betaProbDensity,beta,1);

    % for alpha - eq 15 - stay in logspace for a long time (gammaln(numOfDataPoints) for numerical stability)
    alphaProbDensity = @(x)exp( (numOfClasses-3/2)*log(x) - .5/x + gammaln(x) - gammaln(x+numOfDataPoints) + gammaln(numOfDataPoints));
    alpha = mb_sliceSampling(alphaProbDensity,alpha,1);

    % assign classes to data with precomputed randomnumbers (faster)
    normRandVar = randn(1,numOfDataPoints)/r^.5 + lam;
    gamRandVar  = gamrnd(beta/2,2/w/beta,1,numOfDataPoints);
    randVar     = rand(1,numOfDataPoints);

    for i = 1:numOfDataPoints

        prob = NaN*zeros(1,numOfClasses+1);

        % calculate probability for membership in existing classes
        for j = 1:numOfClasses

            if cs(i) == j; % if data point member of this class
                nij               = numOfPointsInClass(j) - 1;
            else
                nij               = numOfPointsInClass(j);
            end

            if nij > 0
                prob(j) = nij/(numOfDataPoints-1+alpha)*(precs(j))^.5 * exp(-precs(j)/2 * (data(i)-mus(j))^2);
            else % if data point the only member in this class
                prob(j) = alpha/(numOfDataPoints-1+alpha)*(precs(j))^.5 * exp(-precs(j)/2 * (data(i)-mus(j))^2);
            end

        end

        % probability for membership in new class - eq 17
        % sample paramters for this class from priors
        mu_new    = normRandVar(i);
        s_new     = gamRandVar(i);
        prob(end) = alpha/(numOfDataPoints-1+alpha)*s_new^.5*exp(-s_new/2 * (data(i)-mu_new)^2);

        cdf       = cumsum(prob);
        rndNum    = randVar(i)*cdf(end);

        % assign data point to class according to probability
        for j = 1:numOfClasses + 1
           if cdf(j) >= rndNum % bingo this data point will be assigned to class j
               numOfPointsInClass(cs(i)) = numOfPointsInClass(cs(i)) - 1;
               if numOfPointsInClass(cs(i)) < 1
                    fprintf(['Removing class # ' num2str(cs(i)) '\n']);
                    % remove parameters
                    mus(cs(i))                = [];
                    precs(cs(i))              = [];
                    numOfPointsInClass(cs(i)) = [];
                    numOfClasses              = numOfClasses - 1;
                    j(j>cs(i))                = j - 1;
                    % rename all higher classes
                    cs(cs>cs(i))              = cs(cs>cs(i)) - 1;
               end
               cs(i) = j;
               break;
           end
        end

        if j == numOfClasses + 1
            fprintf(['Adding new class # ' num2str(j) '\n']);
            % add parameters
            mus                   = [mus mu_new];
            precs                 = [precs s_new];
            numOfPointsInClass(j) = 1;
            numOfClasses          = numOfClasses + 1;
        else
            numOfPointsInClass(j) = numOfPointsInClass(j) + 1;
        end

    end

            % visualization of final solution - start after burn in
            if iter >= numOfBurnInIter
                delete(currSolutionVisHandle);
                contributionOfCurrentSample2Solution = 0*bins;
                for i = 1:numel(mus)

                    contributionOfCurrentSample2Solution = contributionOfCurrentSample2Solution + ...
                        (precs(i)/2/pi)^.5*numOfPointsInClass(i)/numOfDataPoints*exp(-(bins-mus(i)).^2*precs(i)/2);

                end

                currSolution = (iter-numOfBurnInIter)/(iter-numOfBurnInIter+1)*currSolution + ...
                    1/(iter-numOfBurnInIter+1)*contributionOfCurrentSample2Solution;

                currSolutionVisHandle = plot(bins,currSolution,'c','LineWidth',4);
            end
            drawnow;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% helper functions

function rndVar = mb_mvrnd(p,n)

r      = rand(1,n);
p      = cumsum(p);
if p(end)~=1
    p = p/p(end);
end

rndVar = ones(1,n);
for i = 1:n
    rndVar(i) = find(r(i)<=p,1,'first');
end

end

function randNumbers = mb_sliceSampling(customFunc,xStart,numOfRandNumbers)
% slice sampling implementation according to section 29.7 from
% http://www.inference.phy.cam.ac.uk/itprnn/book.pdf

tic

if nargin < 2
    xStart = rand;
end
if nargin < 3
    numOfRandNumbers = 1;
end

randNumbers = NaN*ones(1,numOfRandNumbers);

for i = 1:numOfRandNumbers

    yStart     = customFunc(xStart);

    randStart  = rand*yStart;

    randStep   = rand;
    stepLength = .1;

    xLeft      = xStart - randStep*stepLength;
    xRight     = xStart + (1-randStep)*stepLength;

    counter = 0;
    while customFunc(xLeft) > randStart
        xLeft = xLeft - stepLength*2^counter;
        counter = counter + 1;
    end
    counter = 0;
    while customFunc(xRight) > randStart
        xRight = xRight + stepLength*2^counter;
        counter = counter + 1;
    end

    while 1
        randNumber = rand*(xRight-xLeft) + xLeft;
        randNumberFuncVal = customFunc(randNumber);
        if randNumberFuncVal > randStart && imag(randNumberFuncVal) == 0
            break;
        else
            if randNumber > xStart
                xRight = randNumber;
            else
                xLeft = randNumber;
            end
        end

        elTime = toc;
        if elTime > 10
            randNumbers = NaN;
            return;
        end

    end

    randNumbers(i) = randNumber;
    xStart = randNumber;

end

end


text:
function spline_test125 ( )

%*****************************************************************************80
%
%% TEST125 tests LEAST_SET_OLD, LEAST_VAL_OLD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2007
%
%  Author:
%
%    John Burkardt
%
  maxdeg = 6;
  ntab = 21;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST125\n' );
  fprintf ( 1, '  LEAST_SET_OLD sets a least squares polynomial,\n' );
  fprintf ( 1, '  LEAST_VAL_OLD evaluates it.\n' );

  for i = 1 : ntab
    xtab(i) = ( ( ntab - i     ) * ( -1.0 )   ...
              + (        i - 1 ) * ( +1.0 ) ) ...
              / ( ntab     - 1 );
    ytab(i) = ( floor ( exp ( xtab(i) ) * 100.0 + 0.5 ) ) / 100.0;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data to be interpolated:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of data values = %d\n', ntab );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       X             Y\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntab
    fprintf ( 1, '%14f  %14f\n', xtab(i), ytab(i) );
  end

  for ndeg = 1 : maxdeg

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Use a polynomial of degree: %d\n', ndeg );
    fprintf ( 1, '\n' );

    [ ptab, b, c, d, eps, ierror ] = least_set_old ( ntab, xtab, ytab, ndeg );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Total approximation error = %f\n', eps );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  X, F(X), P(X), Error\n' );
    fprintf ( 1, '\n' );

    for i = 1 : ntab

      if ( i < ntab )
        jhi = 2;
      else
        jhi = 0;
      end

      for j = 0 : jhi

        if ( i < ntab )

          xval = ( ( 3 - j ) * xtab(i)     ...
                 + (     j ) * xtab(i+1) ) ...
                 / ( 3     );

        else

          xval = xtab(i);

        end

        yval = least_val_old ( xval, ndeg, b, c, d );

        ytrue = ( floor ( exp ( xval ) * 100.0 + 0.5 ) ) / 100.0;

        fprintf ( 1, '%14f  %14f  %14f  %14f\n', xval, yval, ytrue, ...
          yval - ytrue );

      end

    end

  end

  return
end


text:
classdef spinoperator
%SPINOPERATOR   Abstract class for representing the spatial part of differential
%operators for time-dependent PDEs.
%   SPINOPERATOR is a class for representing the spartial part S of a
%   time-dependent PDE of the form u_t = S(u) = Lu + N(u), where L is a linear
%   operator and N is a nonlinear operator. SPINOP (in 1D), SPINOP2 (in 2D),
%   SPINOP3 (in 3D) and SPINOPSPHERE (on the sphere) are full implementations.
%
% See also SPINOP, SPINOP2, SPINOP3, SPINOPSPHERE.

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        domain              % Spatial domain (1x2 DOUBLE in 1D, 1x4 in 2D and on
                            % the sphere, 1x6 in 3D)
        init                % Initial condition (CHEBFUN in 1D, CHEBFUN2 in 2D,
                            % CHEBFUN3 in 3D, SPHEREFUN on the SPHERE,
                            % CHEBMATRIX for systems)
        lin                 % Linear part of the operator (FUNCTION HANDLE)
        nonlin              % Nonlinear part of the operator (FUNCTION HANDLE)
        tspan               % Vector of strictly increasing time
                            % samples (DOUBLE)
    end

    % DEPENDENT PROPERTIES:
    properties ( Access = public, Dependent = true )
        numVars             % Number of unknown functions (1x1 INT)
    end

    % DEPENDENT AND HIDDEN PROPERTIES:
    properties ( Access = public, Hidden = true, Dependent = true )
        nonlinearPartCoeffs % Differential part of the operator in coefficient
                            % space (FUNCTION HANDLE)
        nonlinearPartVals   % Nondifferential part of the operator in value
                            % space (FUNCTION HANDLE)
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% METHODS FOR DEPENDENT PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % METHOD for numVars:
        function nVars = get.numVars(S)

            nVars = nargin(S.lin);

        end

        % METHOD for nonlinearPartCoeffs:
        function Nc = get.nonlinearPartCoeffs(S)

            % Extract the nonlinear part:
            N = S.nonlin;

            % If N is empty, we are done:
            if ( isempty(N) == 1 )
                Nc = @(u) 0*u;
                return
            end

            % Else, get the variables of the workspace:
            func = functions(N);
            wrk = func.workspace{1};
            names = fieldnames(wrk);
            if ( isempty(names) == 0 )
                lengthNames = size(names, 1);
                for iter = 1:lengthNames
                    eval(sprintf('%s = wrk.(names{iter});', names{iter}));
                end
            end

            % Convert FUNCN to a string STRN:
            strN = func2str(N);

            % Get the number of variables NVARS:
            nVars = nargin(N);

            % For scalar equations in 1D, we support nonlinearities of the form
            % diff(f(u),m) with m>=0:
            if ( isa(S, 'spinop') == 1 && nVars == 1 )

                % Compute the differentiation order to get NC=diff(u,m):
                diffOrderTwoOrGreater = regexp(strN, ',\d*', 'match');
                diffOrderOne = regexp(strN, 'diff(', 'match');
                if ( isempty(diffOrderTwoOrGreater) == 0 )
                    diffOrder = diffOrderTwoOrGreater{1}(2:end);
                elseif ( isempty(diffOrderOne) == 0 )
                    diffOrder = '1';
                else
                    diffOrder = '0';
                end
                Nc = ['@(u) diff(u,', diffOrder, ')'];
                Nc = eval(Nc);

            % For scalar equations in 2D/3D and on the sphere, and for systems
            % of equations in 1D/2D/3D and on the sphere, we only support
            % nonlinearities of the form f_i(u_1,...,u_n), i.e., with no
            % differentiation, so Nc=1:
            else

                % Nc=1:
                Nc = 1;

            end

        end

        % METHOD for nonlinearPartVals:
        function Nv = get.nonlinearPartVals(S)

            % Extract the nonlinear part:
            N = S.nonlin;

            % If N is empty, we are done:
            if ( isempty(N) == 1 )
                Nv = @(u) 0*u;
                return
            end

            % Else, get the variables of the workspace:
            func = functions(N);
            wrk = func.workspace{1};
            names = fieldnames(wrk);
            if ( isempty(names) == 0 )
                lengthNames = size(names, 1);
                for iter = 1:lengthNames
                    eval(sprintf('%s = wrk.(names{iter});', names{iter}));
                end
            end

            % Convert FUNCN to a string STRN:
            strN = func2str(N);

            % Get the number of variables NVARS:
            nVars = nargin(N);

            % For scalar equations in 1D, we support nonlinearities of the form
            % diff(f(u),m) with m>=0:
            if ( isa(S, 'spinop') == 1 && nVars == 1 )

                % Get rid of the differentiation part in STRN to get NV=f(u):
                oldString = {'diff', ',\d*)'};
                newString = {'', ')'};
                Nv = regexprep(strN, oldString, newString);
                Nv = eval(Nv);

            % For scalar equations in 2D/3D and on the sphere, and for systems
            % of equations in 1D/2D/3D and on the sphere, we only support
            % nonlinearities of the form f_i(u_1,...,u_n), i.e., with no
            % differentiation, so Nv=N:
            else

                % Nv=N:
                strNv = func2str(N);

                % We're going to relabel the variables, e.g., @(u,v) u + v is
                % going to be relabelled @(u) u(1:length(u)/2) + ...
                % u(length(u)/2+1:end). First, get the names of the variables:
                openParenthesis = strfind(strNv, '(');
                openParenthesis = openParenthesis(1);
                closeParenthesis = strfind(strNv, ')');
                closeParenthesis = closeParenthesis(1);
                variablesNames = strNv(openParenthesis+1:closeParenthesis-1);
                variablesNames = regexp(variablesNames,  ',', 'split');

                % Second, relabel the variables:
                strNv = strNv(closeParenthesis+1:end);
                for iter = 1:nVars
                    idx1 = [num2str(rats((iter-1)/nVars)), '*', 'length(u)', '+1'];
                    idx2 = [num2str(rats(iter/nVars)), '*', 'length(u)'];
                    strNvNew = strrep(strNv, variablesNames{iter}, ...
                        ['u(', idx1, ':', idx2, ',:,:)']);
                    strNv = strNvNew;
                end
                Nv = ['@(u)', strNvNew];
                Nv = eval(Nv);

            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% ABSTRACT AND NON-STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Abstract = true, Static = false )

        % Discretize a SPINOPERATOR:
        [L, Nc] = discretize(S, N)

        % Returns indexes for dealiasing procedure (2/3-rule):
        idx = getDealiasingIndexes(S, N, nVars)

        % Returns the type of CHEBFUN on which the SPINOPERATOR acts:
        out = getChebfunType(S)

        % Returns the transform coeffs -> values:
        F = getCoeffs2ValsTransform(S)

        % Returns the spatial dimension:
        dim = getDimension(S)

        % Returns a grid correspoding to a SPINOPRERATOR object:
        grid = getGrid(S, N, dom)

        % Returns the adequate SPINPREFERENCE object:
        pref = getPreference(S)

        % Returns the transform values -> coeffs:
        F = getVals2CoeffsTransform(S)

        % Returns 1 if the linear part of the SPINOPERATOR is diagonal,
        % 0 otherwise:
        out = isDiag(S)

        % Initialize a movie when solving a PDE specified by a SPINOPERATOR:
        [p, opts] = initializeMovie(S, dt, pref, v, compGrid, plotGrid)

        % Update the movie when solving a PDE specified by a SPINOPERATOR:
        opts = updateMovie(S, dt, p, options, t, v, compGrid, plotGrid)

        % Reshape the data that will be used for constructing the solution at
        % the end of the time-stepping: (For example, for SPINOPSPEHRE, it
        % extracts half of the data, since the data has been doubled-up with
        % the DFS method.)
        data = reshapeData(S, data, nVars)

        % Add the (repeated) endpoints to a periodic grid:
        grid = reshapeGrid(S, grid)

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CONCRETE AND STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Abstract = false, Static = true )

        % Solve a PDE defined by a SPINOPERATOR:
        [uout, tout, computingTime] = solvepde(varargin)

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CONCRETE AND NON-STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Abstract = false, Static = false )

        % Create a contour around each eigenvalue of the linear part of a
        % SPINOPERATOR:
        LR = computeLR(S, dt, L, M)

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% [PYR, INDICES] = buildSCFpyrLevs(LODFT, LOGRAD, XRCOS, YRCOS, ANGLE, HEIGHT, NBANDS)
%
% Recursive function for constructing levels of a steerable pyramid.  This
% is called by buildSCFpyr, and is not usually called directly.

% Original code: Eero Simoncelli, 5/97.
% Modified by Javier Portilla to generate complex bands in 9/97.

function [pyr,pind] = buildSCFpyrLevs(lodft,log_rad,Xrcos,Yrcos,angle,ht,nbands);

if (ht <= 0)

  lo0 = ifft2(ifftshift(lodft));
  pyr = real(lo0(:));
  pind = size(lo0);

else

  bands = zeros(prod(size(lodft)), nbands);
  bind = zeros(nbands,2);

%  log_rad = log_rad + 1;
  Xrcos = Xrcos - log2(2);  % shift origin of lut by 1 octave.

  lutsize = 1024;
  Xcosn = pi*[-(2*lutsize+1):(lutsize+1)]/lutsize;  % [-2*pi:pi]
  order = nbands-1;
  %% divide by sqrt(sum_(n=0)^(N-1)  cos(pi*n/N)^(2(N-1)) )
  %% Thanks to Patrick Teo for writing this out :)
  const = (2^(2*order))*(factorial(order)^2)/(nbands*factorial(2*order));

%
%  Ycosn = sqrt(const) * (cos(Xcosn)).^order;
%
  % analityc version: only take one lobe
  alfa=	mod(pi+Xcosn,2*pi)-pi;
  Ycosn = 2*sqrt(const) * (cos(Xcosn).^order) .* (abs(alfa)<pi/2);

  himask = pointOp(log_rad, Yrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);

  for b = 1:nbands
    anglemask = pointOp(angle, Ycosn, Xcosn(1)+pi*(b-1)/nbands, Xcosn(2)-Xcosn(1));
    banddft = ((-i)^(nbands-1)) .* lodft .* anglemask .* himask;
    band = ifft2(ifftshift(banddft));

%    bands(:,b) = real(band(:));
    % analytic version: full complex value
    bands(:,b)=band(:);
    bind(b,:)  = size(band);
  end

  dims = size(lodft);
  ctr = ceil((dims+0.5)/2);
  lodims = ceil((dims-0.5)/2);
  loctr = ceil((lodims+0.5)/2);
  lostart = ctr-loctr+1;
  loend = lostart+lodims-1;

  log_rad = log_rad(lostart(1):loend(1),lostart(2):loend(2));
  angle = angle(lostart(1):loend(1),lostart(2):loend(2));
  lodft = lodft(lostart(1):loend(1),lostart(2):loend(2));
  YIrcos = abs(sqrt(1.0 - Yrcos.^2));
  lomask = pointOp(log_rad, YIrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);

  lodft = lomask .* lodft;

  [npyr,nind] = buildSCFpyrLevs(lodft, log_rad, Xrcos, Yrcos, angle, ht-1, nbands);

  pyr = [bands(:); npyr];
  pind = [bind; nind];

end

repo: SelfExSR


text:
function [OMNISol,bilevelMILPProblem] = OMNI(model, selectedRxnList, options, constrOpt, measOpt, prevSolutions, verbFlag)
%
%% ***********************NOT WORKING**************************************
%

%function [OMNISol,bilevelMILPProblem] = OMNI(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)

%OMNI Run OMNI in the most general form
%
% OMNI(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileName)
%
%INPUTS
% model                 Structure containing all necessary variables to
%                       describe a stoichiometric model
%   rxns                  Rxns in the model
%   mets                  Metabolites in the model
%   S                     Stoichiometric matrix (sparse)
%   b                     RHS of Sv = b (usually zeros)
%   c                     Objective coefficients
%   lb                    Lower bounds for fluxes
%   ub                    Upper bounds for fluxes
%   rev                    Reversibility of fluxes
% selectedRxnList       List of reactions that can be knocked-out in OMNI
% options               OMNI options
%   numDel                # of bottlenecks
%   numDelSense           Direction of # of bottleneck constraint (G/E/L)
%   vMax                  Max flux
%   solveOMNI             Solve problem within Matlab
%   createGams            Create GAMS input file
%   gamsFile              GAMS input file name
% constrOpt             Explicitly constrained reaction options
%   rxnList               Reaction list
%   values                Values for constrained reactions
%   sense                 Constraint senses for constrained reactions
%                         (G/E/L)
% measOpt               Measured flux options
%   rxnSel                Names of measured reactions
%   values                Flux values of measured reactions
%   weights               Weights for measured fluxes
%
%OPTIONAL INPUTS
% prevSolutions         Previous solutions
% verbFlag              Verbose flag
% solutionFileName      File name for storing temporary solutions
%
%OUTPUTS
% OMNISol               OMNI solution structure
% bilevelMILPProblem    bi-level MILP problem structure used
%
% Markus Herrgard 3/28/05

% Set these for MILP callbacks
global MILPproblemType;
global selectedRxnIndIrrev;
%global rxnList;
global irrev2rev;
%global solutionFileName;
%global biomassRxnID;
%global OMNIKOrxnList;
%global OMNIObjective;
%global OMNIGrowth;
%global solID;

if (nargin < 5)
    prevSolutions = [];
end
if (nargin < 6)
    verbFlag = false;
end
% if (nargin < 7)
%     solutionFileName = 'OMNISolutions.mat';
% else
%     solutionFileName = solutionFileNameTmp;
% end

% Convert to irreversible rxns
[modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);

% Create the index of the previous KO's suggested by OMNI to avoid obtaining the same
% solution again
selPrevSolIrrev = [];
for i = 1:size(prevSolutions,2)
    prevSolRxnList = model.rxns(prevSolutions(:,i)==1);
    selPrevSol = ismember(model.rxns,prevSolRxnList);
    selPrevSolIrrev(:,i) = selPrevSol(irrev2rev);
end

[nMets,nRxns] = size(modelIrrev.S);

% Create matchings for reversible reactions in the set selected for KOs
% This is to ensure that both directions of the reaction are knocked out
selSelectedRxn = ismember(model.rxns,selectedRxnList);
selSelectedRxnIrrev = selSelectedRxn(irrev2rev);
selectedRxnIndIrrev = find(selSelectedRxnIrrev);
cnt = 0;
%prevRxnID = -10;
nSelected = length(selectedRxnIndIrrev);
selRxnCnt = 1;
while selRxnCnt <= nSelected
    rxnID = selectedRxnIndIrrev(selRxnCnt);
    if (matchRev(rxnID)>0)
        cnt = cnt + 1;
        selectedRxnMatch(cnt,1) = selRxnCnt;
        selectedRxnMatch(cnt,2) = selRxnCnt+1;
        selRxnCnt = selRxnCnt + 1;
    end
    selRxnCnt = selRxnCnt + 1;
end

% Set inner constraints for the LP
constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev);
% constrOptIrrev = model;
% constrOptIrrev = [];

% Set objectives for linear and integer parts
cLinear = zeros(nRxns,1);
cInteger = zeros(sum(selSelectedRxnIrrev),1);

% Set the correct objective coefficient (not necessary for OMNI)
% targetRxnID = find(ismember(model.rxns,options.targetRxn));
% targetRxnIDirrev = rev2irrev{targetRxnID}(1);
% cLinear(targetRxnIDirrev) = 1;

% Set measured reaction in objective
sel_meas_rxn = measOpt.rxnSel';
b_meas_rxn = measOpt.values';
wt_meas_rxn = measOpt.weights';
n_m = length(sel_meas_rxn);

% Create selection vector in the decoupled representation
% This is to ensure that the objective function for measured reversible
% reactions is constructed correctly
sel_m = zeros(nRxns,1);
ord_ir = [];
b_meas_tmp = [];
wt_meas_tmp = [];
for i = 1:n_m
    rxn_name = sel_meas_rxn{i};
    rxn_id = find(strcmp(model.rxns,rxn_name));
    if (~isempty(rxn_id)) % Protect against measured fluxes that are not part of the model
        b_meas_tmp = [b_meas_tmp;b_meas_rxn(i)];
        wt_meas_tmp = [wt_meas_tmp;wt_meas_rxn(i)];
        % Reversible rxns
        if (model.rev(rxn_id))
            rxn_id_ir = rev2irrev{rxn_id}(1);
            sel_m(rxn_id_ir) = 1;
            sel_m(rxn_id_ir+1) = -1;
        else
            % Irrev rxns
            rxn_id_ir = rev2irrev{rxn_id};
            sel_m(rxn_id_ir) = 1;
        end
        % Figure out ordering in decoupled representation
        ord_ir = [ord_ir rxn_id_ir];
    end
end
% Get ordering indices
[tmp,ord_ind] = sort(ord_ir);
% Reorder or create weights
if (sum(wt_meas_rxn) == 0)
    measOpts.weights = ones(n_m,1);
else
    measOpts.weights = wt_meas_tmp(ord_ind);
end
% Reorder measured flux values
measOpts.values = b_meas_tmp(ord_ind);

measOpts.rxnSel = sel_m;

% Create the constraint matrices for the bilevel MILP
bilevelMILPProblem = createBilevelMILPproblem(modelIrrev,cLinear,cInteger,selSelectedRxnIrrev,...
    selectedRxnMatch,constrOptIrrev,measOpts,options,selPrevSolIrrev);

% Initial guess (random)
%bilevelMILPProblem.x0 = round(rand(length(bilevelMILPProblem.c),1));
if isfield(options,'initSolution')
    if (length(options.initSolution) > options.numDel | ~all(ismember(options.initSolution,selectedRxnList)))
        warning('Initial solution not valid - starting from a random initial solution')
        bilevelMILPProblem.x0 = [];
    else
        % Set initial integer solution
        selInitRxn = ismember(model.rxns,options.initSolution);
        selInitRxnIrrev = selInitRxn(irrev2rev);
        initRxnIndIrrev = find(selInitRxnIrrev);
        initIntegerSol = ~ismember(selectedRxnIndIrrev,initRxnIndIrrev);
        selInteger = bilevelMILPProblem.vartype == 'B';
        [nConstr,nVar] = size(bilevelMILPProblem.A);
        bilevelMILPProblem.x0 = nan(nVar,1);
        bilevelMILPProblem.x0(selInteger) = initIntegerSol;

%         LPproblem.b = bilevelMILPProblem.b - bilevelMILPProblem.A(:,selInteger)*initIntegerSol;
%         LPproblem.A = bilevelMILPProblem.A(:,bilevelMILPProblem.vartype == 'C');
%         LPproblem.c = bilevelMILPProblem.c(bilevelMILPProblem.vartype == 'C');
%         LPproblem.lb = bilevelMILPProblem.lb(bilevelMILPProblem.vartype == 'C');
%         LPproblem.ub = bilevelMILPProblem.ub(bilevelMILPProblem.vartype == 'C');
%         LPproblem.osense = -1;
%         LPproblem.csense = bilevelMILPProblem.csense;
%         LPsol = solveCobraLP(LPproblem);
%
%         bilevelMILPProblem.x0(~selInteger) = LPsol.full;
    end
else
    bilevelMILPProblem.x0 = [];
end

% Minimize
bilevelMILPProblem.osense = 1;

if (verbFlag)
    [nConstr,nVar] = size(bilevelMILPProblem.A);
    nInt = length(bilevelMILPProblem.intSolInd);
    fprintf('MILP problem with %d constraints %d integer variables and %d continuous variables\n',...
        nConstr,nInt,nVar);
end

bilevelMILPProblem.model = modelIrrev;

% Set these for CPLEX callbacks
MILPproblemType = 'OMNI';
% rxnList = model.rxns;
% biomassRxnID = find(modelIrrev.c==1);
% solID = 0;
% OMNIObjective = [];
% OMNIGrowth = [];
% OMNIKOrxnList = {};

% Solve problem
if (options.solveOMNI)
    OMNISol = solveCobraMILP(bilevelMILPProblem,'printLevel',0);
    if OMNISol.stat~=0
        if (~isempty(OMNISol.cont))
            OMNISol.fluxes = convertIrrevFluxDistribution(OMNISol.cont(1:length(matchRev)),matchRev);
        end
        if (~isempty(OMNISol.int))
            % Figure out the KO reactions
            OMNIRxnInd = selectedRxnIndIrrev(OMNISol.int < 1e-4);
            OMNISol.kos = model.rxns(unique(irrev2rev(OMNIRxnInd)));

%             %sanity check
%             modelTemp = changeRxnBounds(model,OMNISol.kos,0,'b');
%             solTemp = optimizeCbModel(modelTemp);
%             if abs(solTemp.f - OMNISol.obj) > 1e-4
%                 [OMNISol,bilevelMILPProblem] = OMNI(model, selectedRxnList, options, constrOpt, measOpt, prevSolutions, verbFlag);
%                 previous_solutions(:,end+1) = zeros(length(model.rxns),1);
%             end
        end
    else
        OMNISol.fluxes=[];
        OMNISol.kos={};
    end
else
    OMNISol.rxnList = {};
    OMNISol.fluxes = [];
end

author: opencobra

repo: cobratoolbox


text:
function [ xmin, ixmin ] = i4row_min ( m, n, x )

%*****************************************************************************80
%
%% I4ROW_MIN returns the minimums of rows of an I4ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer X(M,N), the array to be examined.
%
%    Output, integer XMIN(M), the minimums of the rows of X.
%
%    Output, integer IXMIN(M); IXMIN(I) is the column of X in which
%    the minimum for row I occurs.
%
  for i = 1 : m

    ixmin(i) = 1;
    xmin(i) = x(i,1);
    for j = 2 : n
      if ( x(i,j) < xmin(i) )
        ixmin(i) = j;
        xmin(i) = x(i,j);
      end
    end

  end

  return
end


text:
function[res] = matFileHasField(path, fieldName)
% [res] = matFileHasField(path, fieldName)
%
% Check whether a mat file includes a variable, without loading it.
% Returns true if it is included or else false.
%
% Copyright by Holger Caesar, 2015

matObj = matfile(path, 'Writable', false);
info = whos(matObj, fieldName);
res = numel(info) == 1;

author: nightrome

repo: matconvnet-calvin


text:
function s32 = i4_to_s32 ( i4 )

%*****************************************************************************80
%
%% I4_TO_S32 converts an I4 to an S32.
%
%  Discussion:
%
%    An I4 is a 32 bit integer.
%
%    An S32 is a 32 character string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I4, the integer to be coded.
%
%    Output, string S32, the string.
%
  s32 = [];

  i4_copy = abs ( i4 );
%
%  Binary digits:
%
  for i = 32 : -1 : 2

    if ( mod ( i4_copy, i4_two ) == 1 )
      s32 = strcat ( '1', s32 );
    else
      s32 = strcat ( '0', s32 );
    end

    i4_copy = floor ( i4_copy / 2 );

  end
%
%  Sign bit
%
  s32 = strcat ( '0', s32 );
%
%  If original number was negative, then reverse all bits.
%
  if ( i4 < 0 )
    for i = 1 : 32
      if ( s32(i) == '0' )
        s32(i) = '1';
      else
        s32(i) = '0';
      end
    end
  end

  return
end


text:
function [ux,u,m] = applypca2(X)
% [ux,u,m] = applypca(X)
%
% finds principal components of
%
% input:
%  x  input data (each column is an input vector)
%
% output:
%  ux  data projected onto pca (1st row=leading principal component)
%  u   the orthogonal matrix with all principal components
%  m   mean of the data
%
% You can recover the original data with
%
% xoriginal=u*ux+repmat(m,1,size(ux,2));
%
% See also pca
%
% copyright by Kilian Q. Weinberger, 2006

% [u,v]=pca(X);
[d,N]  = size(X);

m = mean(X,2);
X = X - m*ones(1,N); % remove mean from data

cc = cov(X',1); % compute covariance
[cvv,cdd] = eig(cc); % compute eignvectors
[zz,ii] = sort(diag(-cdd)); % sort according to eigenvalues
u = cvv(:,ii); % pick leading eigenvectors

ux=u'*X;

repo: MARS-evaluation


text:
function coeffS = setGaussOrder(coeffS,derivativeOrder)
% function coeffS = setOrder(coeffS,derivativeOrder)
%
% APA, 6/13/2018

A1(1) = 1.3530;
B1(1) = 1.8151;
coeffS.W1    = 0.6681;
coeffS.L1    = -1.3932;
A2(1) = -0.3531;
B2(1) = 0.0902;
coeffS.W2    = 2.0787;
coeffS.L2    = -1.3732;

A1(2) = -0.6724;
B1(2) = -3.4327;
A2(2) = 0.6724;
B2(2) = 0.6100;

A1(3) = -1.3563;
B1(3) = 5.2318;
A2(3) = 0.3446;
B2(3) = -2.2355;

switch derivativeOrder

    case 'zero'
        % Zero'th order (Gaussian)
        coeffS.across_scale_normalization = 1.0;

        coeffS = computeDcoeffs(coeffS);

        coeffS.A1 = A1(1);
        coeffS.B1 = B1(1);
        coeffS.A2 = A2(1);
        coeffS.B2 = B2(1);
        coeffS = computeNcoeffs(coeffS);
        coeffS.alpha0 = 2 * coeffS.SN / coeffS.SD - coeffS.N0;
        coeffS.N0 = coeffS.N0 * coeffS.across_scale_normalization / coeffS.alpha0;
        coeffS.N1 = coeffS.N1 * coeffS.across_scale_normalization / coeffS.alpha0;
        coeffS.N2 = coeffS.N2 * coeffS.across_scale_normalization / coeffS.alpha0;
        coeffS.N3 = coeffS.N3 * coeffS.across_scale_normalization / coeffS.alpha0;
        symmetric = true;
        coeffS = computeMcoeffs(coeffS,symmetric);

    case 'second'

        % Second order (Laplacian Of Gaussian)
        %coeffS.across_scale_normalization = coeffS.sigma^2;
        coeffS.across_scale_normalization = 1.0;

        coeffS = computeDcoeffs(coeffS);

        coeffS.A1 = A1(1);
        coeffS.B1 = B1(1);
        coeffS.A2 = A2(1);
        coeffS.B2 = B2(1);
        coeffS.N0 = 0;
        coeffS.N1 = 0;
        coeffS.N2 = 0;
        coeffS.N3 = 0;
        coeffS.SN = 0;
        coeffS.DN = 0;
        coeffS.EN = 0;
        coeffS = computeNcoeffs(coeffS);
        coeffS.N0_0 = coeffS.N0;
        coeffS.N1_0 = coeffS.N1;
        coeffS.N2_0 = coeffS.N2;
        coeffS.N3_0 = coeffS.N3;
        coeffS.SN0 = coeffS.SN;
        coeffS.DN0 = coeffS.DN;
        coeffS.EN0 = coeffS.EN;

        coeffS.A1 = A1(3);
        coeffS.B1 = B1(3);
        coeffS.A2 = A2(3);
        coeffS.B2 = B2(3);
        coeffS.N0 = 0;
        coeffS.N1 = 0;
        coeffS.N2 = 0;
        coeffS.N3 = 0;
        coeffS.SN = 0;
        coeffS.DN = 0;
        coeffS.EN = 0;
        coeffS = computeNcoeffs(coeffS);
        coeffS.N0_2 = coeffS.N0;
        coeffS.N1_2 = coeffS.N1;
        coeffS.N2_2 = coeffS.N2;
        coeffS.N3_2 = coeffS.N3;
        coeffS.SN2 = coeffS.SN;
        coeffS.DN2 = coeffS.DN;
        coeffS.EN2 = coeffS.EN;

        coeffS.beta = -( 2 * coeffS.SN2 - coeffS.SD * coeffS.N0_2 ) / ( 2 * coeffS.SN0 - coeffS.SD * coeffS.N0_0 );
        coeffS.N0 = coeffS.N0_2 + coeffS.beta * coeffS.N0_0;
        coeffS.N1 = coeffS.N1_2 + coeffS.beta * coeffS.N1_0;
        coeffS.N2 = coeffS.N2_2 + coeffS.beta * coeffS.N2_0;
        coeffS.N3 = coeffS.N3_2 + coeffS.beta * coeffS.N3_0;
        coeffS.SN = coeffS.SN2 + coeffS.beta * coeffS.SN0;
        coeffS.DN = coeffS.DN2 + coeffS.beta * coeffS.DN0;
        coeffS.EN = coeffS.EN2 + coeffS.beta * coeffS.EN0;

        coeffS.alpha2  = coeffS.EN * coeffS.SD * coeffS.SD - ...
            coeffS.ED * coeffS.SN * coeffS.SD - ...
            2 * coeffS.DN * coeffS.DD * coeffS.SD + ...
            2 * coeffS.DD * coeffS.DD * coeffS.SN;

        coeffS.alpha2 = coeffS.alpha2 / (coeffS.SD * coeffS.SD * coeffS.SD);

        coeffS.N0 = coeffS.N0 * coeffS.across_scale_normalization / coeffS.alpha2;
        coeffS.N1 = coeffS.N1 * coeffS.across_scale_normalization / coeffS.alpha2;
        coeffS.N2 = coeffS.N2 * coeffS.across_scale_normalization / coeffS.alpha2;
        coeffS.N3 = coeffS.N3 * coeffS.across_scale_normalization / coeffS.alpha2;

        symmetric = true;
        coeffS = computeMcoeffs(coeffS,symmetric);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ xval, yval ] = bc_val ( n, t, xcon, ycon )

%*****************************************************************************80
%
%% BC_VAL evaluates a parameterized Bezier curve.
%
%  Discussion:
%
%    BC_VAL(T) is the value of a vector function of the form
%
%      BC_VAL(T) = ( X(T), Y(T) )
%
%    where
%
%      X(T) = Sum ( 0 <= I <= N ) XCON(I) * BERN(I,N)(T)
%      Y(T) = Sum ( 0 <= I <= N ) YCON(I) * BERN(I,N)(T)
%
%    BERN(I,N)(T) is the I-th Bernstein polynomial of order N
%    defined on the interval [0,1],
%
%    XCON(0:N) and YCON(0:N) are the coordinates of N+1 "control points".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kahaner, Moler, and Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1989.
%
%  Parameters:
%
%    Input, integer N, the order of the Bezier curve, which
%    must be at least 0.
%
%    Input, real T, the point at which the Bezier curve should
%    be evaluated.  The best results are obtained within the interval
%    [0,1] but T may be anywhere.
%
%    Input, real XCON(0:N), YCON(0:N), the X and Y coordinates
%    of the control points.  The Bezier curve will pass through
%    the points ( XCON(0), YCON(0) ) and ( XCON(N), YCON(N) ), but
%    generally NOT through the other control points.
%
%    Output, real XVAL, YVAL, the X and Y coordinates of the point
%    on the Bezier curve corresponding to the given T value.
%
  bval = bp01 ( n, t );

  xval = xcon(1:n+1) * bval(1:n+1)';
  yval = ycon(1:n+1) * bval(1:n+1)';

  return
end


text:
classdef StarDetector < handle
    %STARDETECTOR  The class implements the Star keypoint detector
    %
    % ## References
    % [Agrawal08]:
    % > Motilal Agrawal, Kurt Konolige, and Morten Rufus Blas.
    % > "Censure: Center surround extremas for realtime feature detection and matching".
    % > In Computer Vision-ECCV 2008, pages 102-115. Springer, 2008.
    %
    % See also: cv.StarDetector.StarDetector, cv.FeatureDetector
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = StarDetector(varargin)
            %STARDETECTOR  The full constructor
            %
            %     obj = cv.StarDetector()
            %     obj = cv.StarDetector(..., 'OptionName',optionValue, ...)
            %
            % ## Options
            % * __MaxSize__ maximum size of the features. The following values
            %   are supported: 4, 6, 8, 11, 12, 16, 22, 23, 32, 45, 46, 64,
            %   90, 128. In the case of a different value the result is
            %   undefined. default 45
            % * __ResponseThreshold__ threshold for the approximated laplacian,
            %   used to eliminate weak features. The larger it is, the less
            %   features will be retrieved. default 30
            % * __LineThresholdProjected__ another threshold for the laplacian
            %   to eliminate edges. default 10
            % * __LineThresholdBinarized__ yet another threshold for the
            %   feature size to eliminate edges. The larger the 2nd threshold,
            %   the more points you get. default 8
            % * __SuppressNonmaxSize__ default 5
            %
            % See also: cv.StarDetector.detect
            %
            this.id = StarDetector_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.StarDetector
            %
            if isempty(this.id), return; end
            StarDetector_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = StarDetector_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.StarDetector.empty, cv.StarDetector.load
            %
            StarDetector_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.StarDetector.clear
            %
            b = StarDetector_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.StarDetector.load
            %
            StarDetector_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.StarDetector.save
            %
            StarDetector_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.StarDetector.save, cv.StarDetector.load
            %
            name = StarDetector_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: FeatureDetector
    methods
        function keypoints = detect(this, img, varargin)
            %DETECT  Detects keypoints in an image or image set
            %
            %     keypoints = obj.detect(img)
            %     keypoints = obj.detect(imgs)
            %     [...] = obj.detect(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit/16-bit grayscale image.
            % * __imgs__ Image set (second variant), cell array of images.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. In the first variant, a
            %   1-by-N structure array. In the second variant of the method,
            %   `keypoints{i}` is a set of keypoints detected in `imgs{i}`.
            %
            % ## Options
            % * __Mask__ A mask specifying where to look for keypoints
            %   (optional). It must be a logical or 8-bit integer matrix with
            %   non-zero values in the region of interest. In the second
            %   variant, it is a cell-array of masks for each input image,
            %   `masks{i}` is a mask for `imgs{i}`. Not set by default.
            %
            % See also: cv.StarDetector.StarDetector
            %
            keypoints = StarDetector_(this.id, 'detect', img, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% EEG_COUNTEPOCHS Count how many epochs there are of each type
%
% Usage:
%  >> eeg_countepochs(EEG);
%
% Inputs:
%   EEG            - input dataset
%   epochmarker    - ['type'|'eventtype'] indicates which part of the
%   EEG.epoch structure the different trial types are stored in. Depending
%   on what system the data are from and how they were preprocessed, this
%   may be either EEG.epoch.type or EEG.epoch.eventtype. Defaults to
%   'type'.
%
%  Outputs:
%    sweeps          - Scalar structure showing, for each epoch type
%    (sweeps.types) the number of sweeps in the dataset with that type
%    (sweeps.counts)
%
% Example:
% eeg_countepochs( EEG, 'type' )
% eeg_countepochs( EEG, 'eventtype' )
%
% Author: Stephen Politzer-Ahles, University of Kansas, 2013

% Copyright: Stephen Politzer-Ahles, University of Kansas, 2013
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [sweeps] = eeg_countepochs(EEG, epochmarker);

if nargin < 1
	help eeg_countepochs;
	return;
end

% Initialize an array which will keep counts
clearvars types counts
types{1} = '';
counts = [0];

% Iterate through all trials
for trial=1:length(EEG.epoch)

    % Default 'epochmarker' to 'type' if no input was provided.
    if nargin < 2
        epochmarker = 'type';
    end

    % Look for epochs that have >1 event and find which event is epoched around
    eventindex = 1; % if only 1 event, eventindex=1
    if length(EEG.epoch(trial).eventlatency) > 1
        for eventnum = 1:length(EEG.epoch(trial).eventlatency)
            if EEG.epoch(trial).eventlatency{eventnum}==0
                eventindex = eventnum;
            end % end if
        end % end for
    end % end if

	% Get the trigger number for this trial
    if strcmp(epochmarker,'eventtype')

        % change event type from cell/number event to string
        if iscell(EEG.epoch(trial).eventtype(eventindex))
            type = EEG.epoch(trial).eventtype{eventindex};
        elseif isnumeric(EEG.epoch(trial).eventtype(eventindex))
            type = num2str( EEG.epoch(trial).eventtype(eventindex) );
        else % is char
            type = EEG.epoch(trial).eventtype;
        end
    elseif strcmp(epochmarker,'type')

        % change cell/number event to string
        if iscell(EEG.epoch(trial).type(eventindex))
            type = EEG.epoch(trial).type{eventindex};
        elseif isnumeric(EEG.epoch(trial).type(eventindex))
            type = num2str( EEG.epoch(trial).type(eventindex) );
        else % is char
            type = EEG.epoch(trial).type;
        end
    end; % end if

	% Find the row of the array that corresponds to this trigger (or if this trigger has not been counted yet)
	index_of_type = find( cellfun(@(x) strcmp(x,type), types) );

	% If this trial has a trigger that has already been counted before,
	% 'row_of_type' will be a number. If not, it will be an empty
	% array.
	if index_of_type

		% If we already have a count for this trigger, increment that count by 1
		counts(index_of_type) = counts(index_of_type) + 1;
	else

		% If not, create a new count for this type
		types = [types type];
        counts(end+1) = 1;
	end; % end if

end; % end for

% Remove the unnecessary first item
types = types(2:end);
counts = counts(2:end);

% Sorts
[Y,I] = sort( types );

% Add types and sweeps to a scalar structure
sweeps.types = types(I);
sweeps.counts = counts(I);

end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%% parameters for simulation
gam = [1.58, -0.6];
noise = .5;
T = 455;
framerate = 30;
firerate = 2;
seed = 3;

%% simulate calcium fluorescience
[Y, trueC, trueSpikes] = gen_data(gam, noise, T, framerate, ...
    firerate, [], [], seed);

%% plot results
figure('papersize', [15, 2.5]);
init_fig;
hold on;
col = {[0, 114, 176]; ...
    [0, 158, 115]; ...
    [213, 94, 0]};

% fluorescence trace
plot(1:T, Y(1,:)/3, 'o', 'color',  uint8(col{2}));

% calcium trace
plot(1:T, trueC(1,:)/3, 'color', 'k'); % uint8(col{1}));
% spike train
tsp = find(trueSpikes(1, :));
for m=1:length(tsp)
    plot([1,1]*tsp(m), [0, 1], 'color', uint8(col{3}));
end
axis tight;
xlabel('Time');
ylabel('Fluorescence');
legend('y', 'c', 's');
set(gco, 'fontweigth', 'bold');
saveas(gcf, 'fig/model.pdf');

author: zhoupc


text:
function dataOut=plotMultipleProjectedAmps_PerCondition(view,projectionPhase)
%
% plotMultipleProjectedAmps_PerCondition(view,projectionPhase)
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in a selection of ROIs. All y-axes are made the same. The bar heights
% and a coarse SEM can be obtained from get(gca,'UserData').
% Amplitudes are projected against a single phase (supplied).

% gmb  5/25/98
% bw   2/19/99  Added seY field to the UserData field.
%	    seY is an estimate of the variability in the
%      amplitudes.  It is the SEM of the in the complex
%      (amp*exp(-i*ph)) representation.  The values are
%      computed in vectorMean.m
% fwc   11/07/02 plots data relative to current view
%       added plotting of multiple ROIs
%       ROI selection copied from plotMultipleTSeries.m
% arw   071204 Correctly extracts projections phase from current scan
% Compute means across scans, for all pixels in the
% currently selected ROI.  The seZ value is the mean
% distance from the mean.
% projectionPhase (if supplied) is in radians.

mrGlobals;

refScan = getCurScan(view);
if (~exist('projectionPhase','var'))
    computeProjPhase=1;
else
    if (isnan(projectionPhase))
        % User input requested
        projectionPhase=input('Enter projection phase in radians:');
    end

    computeProjPhase=0;
end

% Select ROIs
nROIs=size(view.ROIs,2);
roiList=cell(1,nROIs);
for r=1:nROIs
    roiList{r}=view.ROIs(r).name;
end
selectedROIs = find(buttondlg('ROIs to Plot',roiList));
nROIs=length(selectedROIs);
if (nROIs==0)
    error('No ROIs selected');
end

% Plot it
selectGraphWin
clf
fontSize = 8;
headerStr = ['Mean Amplitudes'];
set(gcf,'Name',headerStr);

minylim=0;
maxylim=0;
nrows=0;
ncols=0;

nscans = numScans(view);
ROIamps=zeros(nscans,nROIs);
ROIseZ=zeros(nscans,nROIs);
ROImeanPhs=zeros(nscans,nROIs);

for r=1:nROIs

    n=selectedROIs(r);
    view = selectROI(view,n); % is there another way?
    [meanAmps,meanPhs,seZ] = vectorMeans(view);

    if (computeProjPhase)
        projectionPhase=meanPhs(refScan);
    end

    % Compute the amplitude projected onto the reference phase
    meanAmps = meanAmps.*cos(meanPhs-projectionPhase);

    ROIamps(:,r)=meanAmps(:);
    ROIseZ(:,r)=seZ(:);

    ROImeanPhs(:,r)=meanPhs(:);
    %xstr{r}=[view.ROIs(selectedROIs(r)).name];
    xstr{r}=int2str(r);
    roiName{r}=view.ROIs(selectedROIs(r)).name;
    fprintf(['\n#%d :',roiName{r}],r);
end

dataOut.ROIamps=ROIamps;
dataOut.ROIseZ=ROIseZ;
dataOut.ROIname=roiName;

% Now do the plotting

if nscans<=3
    nrows=1;
    ncols=nscans;
    fontSize = 9;
elseif nscans<=8
    nrows=2;
    ncols=ceil(nscans/nrows);
    fontSize = 8;
else
    nrows=ceil(sqrt(nscans));
    ncols=ceil(nscans/nrows);
    fontSize = 6;
end

    scanList = [1:numScans(view)];

for r=1:nscans

    subplot(nrows,ncols,r);

    %plot the bar graph
    if(r==refScan)
        h=mybar(ROIamps(r,:)',ROIseZ(r,:)',xstr,[],[1 0 0]);
    else
        h=mybar(ROIamps(r,:)',ROIseZ(r,:)',xstr,[],[0 0 1]);
    end

    xlabel('ROI','FontSize',fontSize);
    ylabel('Mean Amplitude','FontSize',fontSize);
    set(gca,'FontSize',ceil(fontSize*1.2));
    conditionName{r}=dataTYPES(view.curDataType).scanParams(r).annotation;
    fprintf(['\nCondition #%d :',conditionName{r}],r);

    title(['Condition #: ' int2str(r)]);
    yl=ylim;

    if yl(1)< minylim
        minylim=yl(1);
    end
    if yl(2)> maxylim
        maxylim=yl(2);
    end
end

% 	xlabel('Scan','FontSize',fontSize);
% 	ylabel('Mean Amplitude','FontSize',fontSize);
% 	ylim =get(gca,'YLim');
% 	set(gca,'YLim',ylim*1.1);
%     % slightly bigger title
% 	set(gca,'FontSize',ceil(fontSize*1.2));
% 	title(['ROI: ' view.ROIs(selectedROIs(r)).name]);

%  foo=cell2struct(h,'bar');
% 	hbar=foo.bar(refScan);
% 	set(hbar,'FaceColor','r')

%Save the data in gca('UserData')
data.y =ROIamps(r,:);
data.refScan = refScan;
data.seY = ROIseZ(r,:); % this should probably be adapted

set(gca,'UserData',data);

% give all plots same y-axis

for r=1:nscans
    subplot(nrows,ncols,r);
    ylim([minylim maxylim]);
end

return;


text:
function [logbf] = spm_mlm_posthoc (mlm,c,a)
% Post-hoc model comparison of multivariate linear models
% FORMAT [logbf] = spm_mlm_posthoc (mlm,c,a)
%
% mlm          MLM data structure - see spm_mlm_bayes.m
%              This contains eg. the [p x d] posterior mean regression
%              coefficient matrix mlm.wmean.
%
% c            [k x p*d] contrast matrix defining k-dimensional subspace
% a            hypothesized value (zeros(k,1) by default)
%
% The contrast matrix and hypothesized value define the reduced model.
% The contrast is applied to the vectorised parameters w = vec(mlm.wmean)
%
% The Bayes Factor in favour of the alternative hypothesis over the null
% is computed using a Savage-Dickey ratio (the probability of the
% hypothesized value under the prior versus its probability under the
% posterior)
%
% bf = p(c*w=a|mlm)/p(c*w=a|Y,mlm)
%
% logbf        Log Bayes Factor
%___________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mlm_posthoc.m 4651 2012-02-09 16:03:39Z will $

wpost=mlm.wmean(:);

post_mean=c*wpost;
k=size(post_mean,1);
prior_mean=zeros(k,1);

post_cov=c*mlm.wcov*c';

% Get prior covariance matrix
prec_prior=size(mlm.wcov,1);
for g=1:mlm.prior.groups,
    prec_prior=prec_prior+mlm.prior.group(g).mean_alpha*mlm.prior.group(g).index;
end
w_prior_cov=diag(1./prec_prior);
prior_cov=c*w_prior_cov*c';

if nargin < 3 | isempty (a)
    a=zeros(k,1);
end

% Enforce symmetry - to avoid numerical issues
prior_cov=(prior_cov+prior_cov')/2;
post_cov=(post_cov+post_cov')/2;

% For now just use spm_mvNpdf. This can be made more efficient !
prior_prob=spm_mvNpdf(a,prior_mean,prior_cov);
post_prob=spm_mvNpdf(a,post_mean,post_cov);

logbf=log(prior_prob)-log(post_prob);

author: spm


text:
%%%
%
% Gather src hidden vectors for attention-based models.
% For each sentence ii, we extract a set of vectors to pay attention to
%   srcVecsAll(:, ii, srcPositions(ii)-posWin:srcPositions(ii)+posWin)
%   and put it into srcVecs(:, ii, :). Boundary cases are handled as well.
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
%%%

% TODO move srcMaxLen out
% IMPORTANT: we assume the sentences are reversed here
% srcPositions = srcMaxLen - srcPositions
function [srcVecsSub, h2sInfo] = buildSrcVecs(srcVecsAll, srcPositions, curMask, srcLens, srcMaxLen, params, h2sInfo)
  posWin = params.posWin;
  numPositions = 2*posWin+1;
  [lstmSize, batchSize, numSrcHidVecs] = size(srcVecsAll);

  % masking
  srcPositions(curMask.maskedIds) = [];
  srcLens(curMask.maskedIds) = [];
  unmaskedIds = curMask.unmaskedIds;

  % flatten matrices of size batchSize*numPositions (not exactly batch size but close)
  % init. IMPORTANT: don't swap these two lines
  % assume unmaskedIds = [1, 2, 3, 4, 5], numPositions=3
  indicesSub = reshape(repmat(1:numPositions, length(unmaskedIds), 1), 1, []); % 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3
  unmaskedIds = repmat(unmaskedIds, 1, numPositions); % 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

  % Note: generate multiple sequences of the same lengths without using for loop, see this post for many elegant solutions
  % http://www.mathworks.com/matlabcentral/answers/217205-fast-ways-to-generate-multiple-sequences-without-using-for-loop
  % The below version is the only solution that is faster than for loop (3 times).
  % If startAttnIds = [ 2 4 6 8 10 ] and numPositions=3
  % then indicesAll = [ 2 4 6 8 10 3 5 7 9 11 4 6 8 10 12 ].
  srcLens = repmat(srcLens, [1, numPositions]);
  startAttnIds = srcPositions-posWin;
  indicesAll = reshape(bsxfun(@plus, startAttnIds(:), 0:(numPositions-1)), 1, []);

  % check those that are out of boundaries
  excludeIds = find(indicesAll>numSrcHidVecs | indicesAll<(srcMaxLen-srcLens+1));
  if ~isempty(excludeIds)
    indicesAll(excludeIds) = []; unmaskedIds(excludeIds) = []; indicesSub(excludeIds) = [];
  end

  h2sInfo.indicesAll = indicesAll;
  h2sInfo.unmaskedIds = unmaskedIds;

  if ~isempty(unmaskedIds)
    srcVecsSub = zeroMatrix([lstmSize, batchSize*numPositions], params.isGPU, params.dataType);

    % create linear indices
    h2sInfo.linearIdSub = sub2ind([batchSize, numPositions], unmaskedIds, indicesSub);
    h2sInfo.linearIdAll = sub2ind([batchSize, numSrcHidVecs], unmaskedIds, indicesAll);

    % create srcVecs
    srcVecsAll = reshape(srcVecsAll, lstmSize, []);
    srcVecsSub(:, h2sInfo.linearIdSub) = srcVecsAll(:, h2sInfo.linearIdAll);
    srcVecsSub = reshape(srcVecsSub, [lstmSize, batchSize, numPositions]);
  else
    srcVecsSub = zeroMatrix([lstmSize, batchSize, numPositions], params.isGPU, params.dataType);
    h2sInfo.linearIdSub = [];
    h2sInfo.linearIdAll = [];
  end

  h2sInfo.alignMask = zeroMatrix([batchSize, numPositions], params.isGPU, params.dataType);
  h2sInfo.alignMask(h2sInfo.linearIdSub) = 1;
  h2sInfo.alignMask = h2sInfo.alignMask'; % numPositions * batchSize
  h2sInfo.srcPositions = srcPositions;
end

author: lmthang