text:
function h=subaxis(varargin)
%SUBAXIS Create axes in tiled positions. (just like subplot)
%   Usage:
%      h=subaxis(rows,cols,cellno[,settings])
%      h=subaxis(rows,cols,cellx,celly[,settings])
%      h=subaxis(rows,cols,cellx,celly,spanx,spany[,settings])
%
% SETTINGS: Spacing,SpacingHoriz,SpacingVert
%           Padding,PaddingRight,PaddingLeft,PaddingTop,PaddingBottom
%           Margin,MarginRight,MarginLeft,MarginTop,MarginBottom
%           Holdaxis
%
%           all units are relative (e.g from 0 to 1)
%
%           Abbreviations of parameters can be used.. (Eg MR instead of MarginRight)
%           (holdaxis means that it wont delete any axes below.)
%
%
% Example:
%
%   >> subaxis(2,1,1,'SpacingVert',0,'MR',0);
%   >> imagesc(magic(3))
%   >> subaxis(2,'p',.02);
%   >> imagesc(magic(4))
%
% 2001 / Aslak Grinsted  (Feel free to modify this code.)
f=gcf;

Args=[];
UserDataArgsOK=0;
Args=get(f,'UserData');
if isstruct(Args)
    UserDataArgsOK=isfield(Args,'SpacingHorizontal')&isfield(Args,'Holdaxis')&isfield(Args,'rows')&isfield(Args,'cols');
end
OKToStoreArgs=isempty(Args)|UserDataArgsOK;

if isempty(Args)&(~UserDataArgsOK)
    Args=struct('Holdaxis',0, ...
        'SpacingVertical',0.05,'SpacingHorizontal',0.05, ...
        'PaddingLeft',0,'PaddingRight',0,'PaddingTop',0,'PaddingBottom',0, ...
        'MarginLeft',.1,'MarginRight',.1,'MarginTop',.1,'MarginBottom',.1, ...
        'rows',[],'cols',[]);
end
Args=parseArgs(varargin,Args,{'Holdaxis'},{'Spacing' {'sh','sv'}; 'Padding' {'pl','pr','pt','pb'}; 'Margin' {'ml','mr','mt','mb'}});

if (length(Args.NumericArguments)>1)
    Args.rows=Args.NumericArguments{1};
    Args.cols=Args.NumericArguments{2};
%remove these 2 numerical arguments
    Args.NumericArguments={Args.NumericArguments{3:end}};
end

if OKToStoreArgs
    set(f,'UserData',Args);
end

switch length(Args.NumericArguments)
   case 0
       return % no arguments but rows/cols....
   case 1
      x1=mod((Args.NumericArguments{1}-1),Args.cols)+1; x2=x1;
      y1=floor((Args.NumericArguments{1}-1)/Args.cols)+1; y2=y1;
   case 2
      x1=Args.NumericArguments{1};x2=x1;
      y1=Args.NumericArguments{2};y2=y1;
   case 4
      x1=Args.NumericArguments{1};x2=x1+Args.NumericArguments{3}-1;
      y1=Args.NumericArguments{2};y2=y1+Args.NumericArguments{4}-1;
   otherwise
      error('subaxis argument error')
end

cellwidth=((1-Args.MarginLeft-Args.MarginRight)-(Args.cols-1)*Args.SpacingHorizontal)/Args.cols;
cellheight=((1-Args.MarginTop-Args.MarginBottom)-(Args.rows-1)*Args.SpacingVertical)/Args.rows;
xpos1=Args.MarginLeft+Args.PaddingLeft+cellwidth*(x1-1)+Args.SpacingHorizontal*(x1-1);
xpos2=Args.MarginLeft-Args.PaddingRight+cellwidth*x2+Args.SpacingHorizontal*(x2-1);
ypos1=Args.MarginTop+Args.PaddingTop+cellheight*(y1-1)+Args.SpacingVertical*(y1-1);
ypos2=Args.MarginTop-Args.PaddingBottom+cellheight*y2+Args.SpacingVertical*(y2-1);

if Args.Holdaxis
    h=axes('position',[xpos1 1-ypos2 xpos2-xpos1 ypos2-ypos1]);
else
    h=subplot('position',[xpos1 1-ypos2 xpos2-xpos1 ypos2-ypos1]);
end

set(h,'box','on');
%h=axes('position',[x1 1-y2 x2-x1 y2-y1]);
set(h,'units',get(gcf,'defaultaxesunits'));
set(h,'tag','subaxis');

if (nargout==0) clear h; end;

author: konstantinberlin

repo: malware-windows-audit-log-detection


text:
function rot = affineExtractRotation(A)
%
% rot = affineExtractRotation(A)
%
% Extracts the rotation component from an affine transform matrix.
% Returns the 3x3 rotation matrix. To get the three Euler angles, use
% affineDecompose.
%
% (See also affineBuild, affineDecompose.)
%
% HISTORY:
%   2007.05.11 RFD (bob@white.stanford.edu) wrote it.

% Extract 3x3 rotation/scale/skew component
rot = A(1:3,1:3);
C = chol(rot'*rot);
s = diag(C)';
if det(rot)<0, s(1) = -s(1); end % Fix negative determinants

% Remove skews and scales, leaving just the rotations
C = diag(diag(C))\C;
k = C([4,7,8]);
sk = [ s(1)  0    0;
         0   s(2)  0;
         0    0   s(3)] ...
     *[  1   k(1) k(2);
         0   1    k(3);
         0   0     1  ];
rot = rot/sk;

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

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

% Last Modified by GUIDE v2.5 12-Nov-2008 20:42:52

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

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

close('credits')

author: Sable

repo: mcbench-benchmarks


text:
%% S2AxisFieldHarmonic
%
% S2AxisFieldharmonic handles axis fields on the sphere.
% Axis can be understood as three-dimensional vectors without direction or length.
%%
% S2AxisFieldHarmonic handles functions of the form
%
% $$ f\colon \bf S^2\to\bf R^3_{/<\pm \mathrm{Id}>}. $$
%
%% Defining a S2AxisFieldHarmonic
%
%%
% *Definition via function values*
%
% At first you need some vertices
nodes = equispacedS2Grid('points', 1e5);
nodes = nodes(:);
%%
% Next you define function values for the vertices
y = vector3d(sin(5*nodes.x), 1, nodes.y, 'antipodal');
%%
% Now the actual command to get |sAF1| of type <S2AxisFieldHarmonic.S2AxisFieldHarmonic |S2AxisFieldHarmonic|>
sAF1 = S2AxisFieldHarmonic.approximation(nodes, y)

%%
% *Definition via function handle*
%
% If you have a function handle for the function you could create a
% |S2AxisFieldHarmonic| via quadrature. At first lets define a function
% handle which takes <vector3d.vector3d.html |vector3d|> as an argument and
% returns antipodal <vector3d.vector3d.html |vector3d|>:

f = @(v) vector3d(v.x, v.y, 0*v.x, 'antipodal');
%%
% Now you can call the quadrature command to get |sAF2| of type <S2AxisFieldHarmonic.S2AxisFieldHarmonic |S2AxisFieldHarmonic|>
sAF2 = S2AxisFieldHarmonic.quadrature(@(v) f(v))

%% Visualization
%
% One can use the default |plot|-command

plot(sAF1);
%%
% * same as quiver(sAF1)

%%
% or use the 3D plot of a sphere with the axis on itself
clf;
quiver3(sAF2);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [img,nrg] = rayCubeAnalytic(L,mat,air,Xsrc,Xmic,Rmax)
%+========================================================================+
%|                                                                        |
%|           OPENRAY - LIBRARY FOR TRI-DIMENSIONAL RAY TRACING            |
%|           openRay is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : rayCubeAnalytic.m                             |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Analytical source image method for a cube     |
%|  `---'  |                                                              |
%+========================================================================+

% Indices
i = -30:30;
a = i + 0.5 - 0.5*(-1).^i;
b = (-1).^i;

% Position relatives des sources images
x       = b*Xsrc(1) + a*L(1) - Xmic(1);
y       = b*Xsrc(2) + a*L(2) - Xmic(2);
z       = b*Xsrc(3) + a*L(3) - Xmic(3);
[x,y,z] = meshgrid(x,y,z);
img     = [x(:),y(:),z(:)];

% Dissipation de l'energie par propagation spherique
dst = sqrt(sum(img.^2,2));
nrg = 1./(dst.^2);

% Dissipation de l'energie par les parois
nrg   = nrg * ones(1,size(mat,2));
rfl   = (1-mat);
i0    = abs(0.5*i - 0.25 + 0.25*(-1).^i);
i1    = abs(0.5*i + 0.25 - 0.25*(-1).^i);
for j = 1:size(nrg,2)
    [rx,ry,rz] = meshgrid(...
        (rfl(1,j).^i0).*(rfl(2,j).^i1),...
        (rfl(3,j).^i0).*(rfl(4,j).^i1),...
        (rfl(5,j).^i0).*(rfl(6,j).^i1));
    nrg(:,j) = (rx(:).*ry(:).*rz(:)) .* nrg(:,j);
end

% Energie dissipee par l'air en fonction de la distance
dst = sqrt(sum(img.^2,2));
nrg = nrg .* exp(-dst*air);

% Selections des images selon l'ordre initial
ind = find(dst<Rmax);
dst = dst(ind);
img = img(ind,:);
nrg = nrg(ind,:);

% Sort in phase
[~,ind] = sort(dst);
img       = img(ind,:);
nrg       = nrg(ind,:);

% Normalisation de l'energie
nrg = nrg./max(max(nrg));

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = isempty(varargin)
%ISEMPTY   True for empty SPHEREFUN.
%   ISEMPTY(F) returns 1 if F is an empty SPHEREFUN object and 0 otherwise.

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

[varargout{1:nargout}] = isempty@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = r8pp_indicator ( n )

%*****************************************************************************80
%
%% R8PP_INDICATOR sets up a R8PP indicator matrix.
%
%  Discussion:
%
%    The R8PP storage format is appropriate for a symmetric positive
%    definite matrix.  Only the upper triangle of the matrix is stored,
%    by successive partial columns, in an array of length (N*(N+1))/2,
%    which contains (A11,A12,A22,A13,A23,A33,A14,...,ANN)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Output, real A((N*(N+1))/2), the R8PP matrix.
%
  fac = 10^( i4_log_10 ( n ) + 1 );

  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      a(k) = fac * i + j;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [img, v2smap]=surf2vol(node,face,xi,yi,zi,varargin)
%
% [img, v2smap]=surf2vol(node,face,xi,yi,zi,'options',values,...)
%
% convert a triangular surface to a shell of voxels in a 3D image
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%	 node: node list of the triangular surface, 3 columns for x/y/z
%	 face: triangle node indices, each row is a triangle
%              if face contains the 4th column, it indicates the label of
%              the face triangles (each face componment must be closed); if
%              face contains 5 columns, it stores a tetrahedral mesh with
%              labels, where the first 4 columns are the element list and
%              the last column is the element label;
%	 xi,yi,zi: x/y/z grid for the resulting volume
%        options: 'fill', if set to 1, the enclosed voxels are labeled by 1
%                 'label', if set to 1, the enclosed voxels are labeled by
%                          the corresponding label of the face or element;
%                          setting 'label' to 1 also implies 'fill'.
%
% output:
%	 img: a volumetric binary image at position of ndgrid(xi,yi,zi)
%        v2smap (optional): a 4x4 matrix denoting the Affine transformation to map
%             the voxel coordinates back to the mesh space.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fprintf(1,'converting a closed surface to a volumetric binary image ...\n');

opt=varargin2struct(varargin{:});
label=jsonopt('label',0,opt);

elabel=[1];
if(size(face,2)>=4)
        elabel=unique(face(:,end));
        if(size(face,2)==5)
            label=1;
            el=face;
            face=[];
            for i=1:length(elabel)
                    fc=volface(el(el(:,5)==elabel(i),1:4));
                    fc(:,4)=elabel(i);
                    face=[face ; fc];
            end
        end
else
        fc=face;
end

img=zeros(length(xi),length(yi),length(zi));

for i=1:length(elabel)
        if(size(face,2)==4)
            fc=face(face(:,4)==elabel(i),1:3);
        end
        im=surf2volz(node(:,1:3),fc(:,1:3),xi,yi,zi);
        im=im | shiftdim(surf2volz(node(:,[3 1 2]),fc(:,1:3),zi,xi,yi),1);
        im=im | shiftdim(surf2volz(node(:,[2 3 1]),fc(:,1:3),yi,zi,xi),2);

        v2smap=[];

        % here we assume the grid is uniform; surf2vol can handle non-uniform grid,
        % but the affine output is not correct in this case

        if(jsonopt('fill',0,opt) || label)
                im=imfill(im,'holes');
                if(label)
                    im=im*elabel(i);
                end
        end
        img=max(im,img);
end

if(nargout>1)
        dlen=abs([xi(2)-xi(1) yi(2)-yi(1) zi(2)-zi(1)]);
        p0=min(node);
        offset=p0;
        v2smap=diag(abs(dlen));
        v2smap(4,4)=1;
        v2smap(1:3,4)=offset';
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% mtimesx_build compiles mtimesx.c with BLAS libraries
%******************************************************************************
%
%  MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
%  Function:    mtimesx_build
%  Filename:    mtimesx_build.m
%  Programmer:  James Tursa
%  Version:     1.40
%  Date:        October 4, 2010
%  Copyright:   (c) 2009, 2010 by James Tursa, All Rights Reserved
%
%  This code uses the BSD License:
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  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 OWNER 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.
%
%--
%
%  mtimesx_build compiles mtimesx.c and mtimesx_RealTimesReal.c with the BLAS
%  libraries libmwblas.lib (if present) or libmwlapack.lib (if libmwblas.lib
%  is not present). This function basically works as follows:
%
%  - Opens the current mexopts.bat file in the directory [prefdir], and
%    checks to make sure that the compiler selected is cl or lcc. If it
%    is not, then a warning is issued and the compilation continues with
%    the assumption that the microsoft BLAS libraries will work.
%
%  - Looks for the file libmwblas.lib or libmwlapack.lib files in the
%    appropriate directory: [matlabroot '\extern\lib\win32\microsoft']
%                       or  [matlabroot '\extern\lib\win32\lcc']
%                       or  [matlabroot '\extern\lib\win64\microsoft']
%                       or  [matlabroot '\extern\lib\win64\lcc']
%
%  - Changes directory to the directory of the file mtimesx.m.
%
%  - Compiles mtimesx.c (which includes mtimesx_RealTimesReal.c) along with
%    either libmwblas.lib or libmwlapack.lib depending on the version of
%    MATLAB. The resulting exedcutable mex file is placed in the same
%    directory as the source code. The files mtimesx.m, mtimesx.c, and
%    mtimesx_RealTimesReal.c must all be in the same directory.
%
%  - Changes the directory back to the original directory.
%
%  Change Log:
%  2009/Sep/27 --> 1.00, Initial Release
%  2010/Feb/15 --> 1.10, Fixed largearrardims typo to largeArrayDims
%  2010/Oct/04 --> 1.40, Updated support for OpenMP compiling
%
%**************************************************************************

function mtimesx_build(x)
disp(' ');
disp('... Build routine for mtimesx');

TRUE = 1;
FALSE = 0;

%\
% Check for number of inputs & outputs
%/

noopenmp = FALSE;
if( nargin == 1 )
    if( isequal(upper(x),'NOOPENMP') )
        noopenmp = TRUE;
    else
        error('Invalid input.');
    end
elseif( nargin ~= 0 )
    error('Too many inputs. Expected none.');
end
if( nargout ~= 0 )
    error('Too many outputs. Expected none.');
end

%\
% Check for non-PC
%/

disp('... Checking for PC');
try
    % ispc does not appear in MATLAB 5.3
    pc = ispc ;
catch
    % if ispc fails, assume we are on a Windows PC if it's not unix
    pc = ~isunix ;
end

if( ~pc )
    disp('building linux version');
    mex('mtimesx.c','-DDEFINEUNIX','-largeArrayDims','-lmwblas');
    return;
end

%\
% Check to see that mtimesx.c source code is present
%/

disp('... Finding path of mtimesx C source code files');
try
    mname = mfilename('fullpath');
catch
    mname = mfilename;
end
cname = [mname(1:end-6) '.c'];
if( isempty(dir(cname)) )
    disp('Cannot find the file mtimesx.c in the same directory as the');
    disp('file mtimesx_build.m. Please ensure that they are in the same');
    disp('directory and try again. The following file was not found:');
    disp(' ');
    disp(cname);
    disp(' ');
    error('Unable to compile mtimesx.c');
end
disp(['... Found file mtimesx.c in ' cname]);

%\
% Check to see that mtimesx_RealTimesReal.c source code is present
%/

rname = [mname(1:end-13) 'mtimesx_RealTimesReal.c'];
if( isempty(dir(rname)) )
    disp('Cannot find the file mtimesx_RealTimesReal.c in the same');
    disp('directory as the file mtimesx_build.m. Please ensure that');
    disp('they are in the same directory and try again. The');
    disp('following file was not found:');
    disp(' ');
    disp(rname);
    disp(' ');
    error('Unable to compile mtimesx.c');
end
disp(['... Found file mtimesx_RealTimesReal.c in ' rname]);

%\
% Open the current mexopts.bat file
%/

matlab_ver = version('-release');
if str2num(matlab_ver(1:4)) >= 2014
    mexopts = [prefdir '\mex_C++_win64.xml'];
else
    mexopts = [prefdir '\mexopts.bat'];
end
fid = fopen(mexopts);
if( fid == -1 )
    error('A C/C++ compiler has not been selected with mex -setup');
end
disp(['... Opened the mexopts.bat file in ' mexopts]);
disp('... Reading the mexopts.bat file to find the compiler and options used.');

%\
% Check for the correct compiler selected.
%/

ok_cl = FALSE;
ok_lcc = FALSE;
omp_option = '';
compiler = '(unknown)';
compilername = '';
while( TRUE )
    tline = fgets(fid);
    if( isequal(tline,-1) )
        break;
    else
        if( isempty(compilername) )
            y = findstr(tline,'OPTS.BAT');
            if( ~isempty(y) )
                x = findstr(tline,'rem ');
                if( ~isempty(x) )
                    compilername = tline(x+4:y-1);
                end
            end
        end
        x = findstr(tline,'COMPILER=lcc');
        if( ~isempty(x) )
            ok_lcc = TRUE;
            libdir = 'lcc';
            compiler = 'LCC';
            disp(['... ' compiler ' is the selected compiler']);
            break;
        end
        x = findstr(tline,'COMPILER=cl');
        if( ~isempty(x) )
            ok_cl = TRUE;
            libdir = 'microsoft';
            compiler = ['Microsoft_' compilername '_cl'];
            omp_option = ' /openmp';
            disp(['... ' compiler ' is the selected compiler']);
            break;
        end
        x = findstr(tline,'COMPILER=bcc32');
        if( ~isempty(x) )
            ok_cl = TRUE;
            libdir = 'microsoft';
            compiler = ['Borland_' compilername '_bcc32'];
            disp(['... ' compiler ' is the selected compiler']);
            disp('... Assuming that Borland will link with Microsoft libraries');
            break;
        end
        x = findstr(tline,'COMPILER=icl');
        if( ~isempty(x) )
            ok_cl = TRUE;
            if( pc )
                omp_option = ' -Qopenmp';
            else
                omp_option = ' -openmp';
            end
            libdir = 'microsoft';
            compiler = ['Intel_' compilername '_icl'];
            disp(['... ' compiler ' is the selected compiler']);
            disp('... Assuming that Intel will link with Microsoft libraries');
            break;
        end
        x = findstr(tline,'COMPILER=wc1386');
        if( ~isempty(x) )
            ok_cl = TRUE;
            libdir = 'microsoft';
            compiler = ['Watcom_' compilername '_wc1386'];
            disp(['... ' compiler ' is the selected compiler']);
            disp('... Assuming that Watcom will link with Microsoft libraries');
            break;
        end
        x = findstr(tline,'COMPILER=gcc');
        if( ~isempty(x) )
            ok_cl = TRUE;
            libdir = 'microsoft';
            omp_option = ' -fopenmp';
            compiler = 'GCC';
            disp(['... ' compiler ' is the selected compiler']);
            disp('... Assuming that GCC will link with Microsoft libraries');
            break;
        end
    end
end
fclose(fid);

%\
% MS Visual C/C++ or lcc compiler has not been selected
%/

if( ~(ok_cl | ok_lcc) )
    warning('... Supported C/C++ compiler has not been selected with mex -setup');
    warning('... Assuming that Selected Compiler will link with Microsoft libraries');
    warning('... Continuing at risk ...');
    libdir = 'microsoft';
end

%\
% If an OpenMP supported compiler is potentially present, make sure that the
% necessary compile option is present in the mexopts.bat file on the COMPFLAGS
% line.  If necessary, build a new mexopts.bat file with the correct option
% added to the COMPFLAGS line.
%/

while( TRUE )
ok_openmp = FALSE;
ok_compflags = FALSE;
xname = '';
if( isempty(omp_option) )
    disp('... OpenMP compiler not detected ... you may want to check this website:');
    disp('    http://openmp.org/wp/openmp-compilers/');
elseif( noopenmp )
    disp(['... OpenMP compiler potentially detected, but not checking for ''' omp_option ''' compile option']);
else
    disp('... OpenMP compiler potentially detected');
    disp(['... Checking to see that the ''' omp_option ''' compile option is present']);
    fid = fopen(mexopts);
    while( TRUE )
        tline = fgets(fid);
        if( isequal(tline,-1) )
            break;
        else
            x = findstr(tline,'set COMPFLAGS');
            if( ~isempty(x) )
                ok_compflags = TRUE;
                x = findstr(tline,omp_option);
                if( ~isempty(x) )
                    ok_openmp = TRUE;
                end
                break;
            end
        end
    end
    fclose(fid);
    if( ~ok_compflags )
        warning(['... COMPFLAGS line not found ... ''' omp_option ''' will not be added.']);
    elseif( ~ok_openmp )
        disp(['... The ''' omp_option ''' compile option is not present ... adding it']);
        xname = [mname(1:end-6) '_mexopts.bat'];
        disp(['... Creating custom options file ' xname ' with the ''' omp_option ''' option added.']);
        fid = fopen(mexopts);
        fidx = fopen(xname,'w');
        if( fidx == -1 )
            xname = '';
            warning(['... Unable to create custom mexopts.bat file ... ''' omp_option ''' will not be added']);
        else
            while( TRUE )
                tline = fgets(fid);
                if( isequal(tline,-1) )
                    break;
                else
                    x = findstr(tline,'set COMPFLAGS');
                    if( ~isempty(x) )
                        n = numel(tline);
                        e = n;
                        while( tline(e) < 32 )
                            e = e - 1;
                        end
                        tline = [tline(1:e) omp_option tline(e+1:n)];
                    end
                    fwrite(fidx,tline);
                end
            end
            fclose(fidx);
        end
        fclose(fid);
    end
end

%\
% Construct full file name of libmwblas.lib and libmwlapack.lib. Note that
% not all versions have both files. Earlier versions only had the lapack
% file, which contained both blas and lapack routines.
%/

comp = computer;
mext = mexext;
lc = length(comp);
lm = length(mext);
cbits = comp(max(1:lc-1):lc);
mbits = mext(max(1:lm-1):lm);
if( isequal(cbits,'64') | isequal(mbits,'64') )
    compdir = 'win64';
    largearraydims = '-largeArrayDims';
else
    compdir = 'win32';
    largearraydims = '';
end

lib_blas = [matlabroot '\extern\lib\' compdir '\' libdir '\libmwblas.lib'];
d = dir(lib_blas);
if( isempty(d) )
    disp('... BLAS library file not found, so linking with the LAPACK library');
    lib_blas = [matlabroot '\extern\lib\' compdir '\' libdir '\libmwlapack.lib'];
end
disp(['... Using BLAS library lib_blas = ''' lib_blas '''']);

%\
% Save old directory and change to source code directory
%/

cdold = cd;
if( length(mname) > 13 )
    cd(mname(1:end-13));
end

%\
% Do the compile
%/

disp('... Now attempting to compile ...');
disp(' ');
try
    if( isunix )
        if( isempty(largearraydims) )
            if( isempty(xname) )
                disp(['mex(''-DDEFINEUNIX'',''' cname ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex('-DDEFINEUNIX',cname,lib_blas,['-DCOMPILER=' compiler]);
            else
                disp(['mex(''-f'',''' xname ''',''-DDEFINEUNIX'',''' cname ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex('-f',xname,'-DDEFINEUNIX',cname,lib_blas,['-DCOMPILER=' compiler]);
            end
        else
            if( isempty(xname) )
                disp(['mex(''-DDEFINEUNIX'',''' cname ''',''' largearraydims ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex('-DDEFINEUNIX',largearraydims,cname,lib_blas,['-DCOMPILER=' compiler]);
            else
                disp(['mex(''-f'',''' xname ''',''-DDEFINEUNIX'',''' cname ''',''' largearraydims ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex('-f',xname,'-DDEFINEUNIX',largearraydims,cname,lib_blas,['-DCOMPILER=' compiler]);
            end
        end
    else
        if( isempty(largearraydims) )
            if( isempty(xname) )
                disp(['mex(''' cname ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex(cname,lib_blas,['-DCOMPILER=' compiler]);
            else
                disp(['mex(''-f'',''' xname ''',''' cname ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex('-f',xname,cname,lib_blas,['-DCOMPILER=' compiler]);
            end
        else
            if( isempty(xname) )
                disp(['mex(''' cname ''',''' largearraydims ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex(cname,largearraydims,lib_blas,['-DCOMPILER=' compiler]);
            else
                disp(['mex(''-f'',''' xname ''',''' cname ''',''' largearraydims ''',lib_blas,''-DCOMPILER=' compiler ''')']);
                disp(' ');
                mex('-f',xname,cname,largearraydims,lib_blas,['-DCOMPILER=' compiler]);
            end
        end
    end
    disp('... mex mtimesx.c build completed ... you may now use mtimesx.');
    disp(' ');
    mtimesx;
    break;
catch
    if( noopenmp )
        cd(cdold);
        disp(' ');
        disp('... Well, *that* didn''t work either!');
        disp(' ');
        disp('The mex command failed. This may be because you have already run');
        disp('mex -setup and selected a non-C compiler, such as Fortran. If this');
        disp('is the case, then rerun mex -setup and select a C/C++ compiler.');
        disp(' ');
        error('Unable to compile mtimesx.c');
    else
        disp(' ');
        disp('... Well, *that* didn''t work ...');
        disp(' ');
        if( isequal(omp_option,' /openmp') )
            disp('This may be because an OpenMP compile option was added that the');
            disp('compiler did not like. For example, the Standard versions of the');
            disp('Microsoft C/C++ compilers do not support OpenMP, only the');
            disp('Professional versions do. Attempting to compile again but this');
            disp(['time will not add the ''' omp_option ''' option.'])
        else
            disp('This may be because an OpenMP compile option was added that the');
            disp('compiler did not like. Attempting to compile again, but this time');
            disp(['will not add the ''' omp_option ''' option.'])
        end
        disp(' ');
        noopenmp = TRUE;
    end
end
end

%\
% Restore old directory
%/

cd(cdold);

return
end

author: martin-danelljan

repo: ECO


text:
function [Rlatt,Rrp,ind_rp,eff] = latmio_und(R,ITER,D)
%LATMIO_UND     Lattice with preserved degree distribution
%
%   [Rlatt,Rrp,ind_rp,eff] = latmio_und(R,ITER,D);
%
%   This function "latticizes" an undirected network, while preserving the
%   degree distribution. The function does not preserve the strength
%   distribution in weighted networks.
%
%   Input:      R,      undirected (binary/weighted) connection matrix
%               ITER,   rewiring parameter
%                       (each edge is rewired approximately ITER times)
%               D,      distance-to-diagonal matrix
%
%   Output:     Rlatt,  latticized network in original node ordering
%               Rrp,    latticized network in node ordering used for
%                       latticization
%               ind_rp, node ordering used for latticization
%               eff,    number of actual rewirings carried out
%
%   References: Maslov and Sneppen (2002) Science 296:910
%               Sporns and Zwi (2004) Neuroinformatics 2:145
%
%   2007-2012
%   Mika Rubinov, UNSW
%   Jonathan Power, WUSTL
%   Olaf Sporns, IU

%   Modification History:
%   Jun 2007: Original (Mika Rubinov)
%   Apr 2008: Edge c-d is flipped with 50% probability, allowing to explore
%             all potential rewirings (Jonathan Power)
%   Feb 2012: limit on number of attempts, distance-to-diagonal as input,
%             count number of successful rewirings (Olaf Sporns)
%   Feb 2012: permute node ordering on each run, to ensure lattices are
%             shuffled across mutliple runs (Olaf Sporns)

n=size(R,1);

% randomly reorder matrix
ind_rp = randperm(n);
R = R(ind_rp,ind_rp);

% create 'distance to diagonal' matrix
if nargin<3 %if D is not specified by user
    D=zeros(n);
    u=[0 min([mod(1:n-1,n);mod(n-1:-1:1,n)])];
    for v=1:ceil(n/2)
        D(n-v+1,:)=u([v+1:n 1:v]);
        D(v,:)=D(n-v+1,n:-1:1);
    end
end
%end create

[i,j]=find(tril(R));
K=length(i);
ITER=K*ITER;

% maximal number of rewiring attempts per 'iter'
maxAttempts= round(n*K/(n*(n-1)/2));
% actual number of successful rewirings
eff = 0;

for iter=1:ITER
    att=0;
    while (att<=maxAttempts)                                     %while not rewired
        while 1
            e1=ceil(K*rand);
            e2=ceil(K*rand);
            while (e2==e1),
                e2=ceil(K*rand);
            end
            a=i(e1); b=j(e1);
            c=i(e2); d=j(e2);

            if all(a~=[c d]) && all(b~=[c d]);
                break           %all four vertices must be different
            end
        end

        if rand>0.5
            i(e2)=d; j(e2)=c; 	%flip edge c-d with 50% probability
            c=i(e2); d=j(e2); 	%to explore all potential rewirings
        end

        %rewiring condition
        if ~(R(a,d) || R(c,b))
            %lattice condition
            if (D(a,b)*R(a,b)+D(c,d)*R(c,d))>=(D(a,d)*R(a,b)+D(c,b)*R(c,d))
                R(a,d)=R(a,b); R(a,b)=0;
                R(d,a)=R(b,a); R(b,a)=0;
                R(c,b)=R(c,d); R(c,d)=0;
                R(b,c)=R(d,c); R(d,c)=0;

                j(e1) = d;          %reassign edge indices
                j(e2) = b;
                eff = eff+1;
                break;
            end %lattice condition
        end %rewiring condition
        att=att+1;
    end %while not rewired
end %iterations

% lattice in node order used for latticization
Rrp = R;
% reverse random permutation of nodes
[~,ind_rp_reverse] = sort(ind_rp);
Rlatt = Rrp(ind_rp_reverse,ind_rp_reverse);


text:
function ndsimKernDisplay(kern, spacing)

% NDSIMKERNDISPLAY Display parameters of the NDSIM kernel.
% FORMAT
% DESC displays the parameters of the single input motif
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : ndsimKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFICATIONS : David Luengo, 2009
%
% COPYRIGHT : Antti Honkela, 2011

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
if kern.isStationary
    fprintf('Stationary version of the kernel\n');
else
    fprintf('Non-stationary version of the kernel\n');
end
fprintf(spacing);
if isfield(kern, 'isNormalised') && kern.isNormalised
    fprintf('Normalised version of the kernel\n');
else
    fprintf('Unnormalised version of the kernel\n');
end
fprintf(spacing);
if isfield(kern, 'isNegativeS') && kern.isNegativeS
    fprintf('Sensitivities allowed to be negative.\n');
else
    fprintf('Sensitivities constrained positive.\n');
end
fprintf(spacing);
fprintf('NDSIM inverse width: %2.4f (length scale %2.4f)\n', ...
        kern.inverseWidth, 1/sqrt(kern.inverseWidth));
fprintf(spacing);
fprintf('NDSIM Variance: %2.4f\n', kern.variance)
if isfield(kern, 'gaussianInitial') && kern.gaussianInitial,
  fprintf(spacing);
  fprintf('NDSIM Initial Variance: %2.4f\n', kern.initialVariance)
end
if isfield(kern, 'isNegativeS') && kern.isNegativeS
  fprintf(spacing);
  fprintf('NDSIM Sensitivity: %2.4f\n', kern.sensitivity)
end
%fprintf(spacing);
%fprintf('NDSIM delay: %2.4f\n', kern.delay)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


%   getBox:  Given the current state, returns a number from 1 to 162
%             designating the region of the state space encompassing the current state.
%             Returns a value of -1 if a failure state is encountered.

function box = getBox(theta,thetaDot,x,xDot)
theta = rad2deg(theta);
thetaDot = rad2deg(thetaDot);
if (x < -2.4 || x > 2.4  || theta < -12 || theta > 12)
    box = -1;
else

if (theta<-6&&theta>=-12)
	thetaBucket = 1;
elseif (theta<-1&&theta>=-6)
	thetaBucket = 2;
elseif (theta<0&&theta>=-1)
	thetaBucket = 3;
elseif (theta<1&&theta>=0)	% zero included
	thetaBucket = 4;
elseif (theta<6&&theta>=1)
	thetaBucket = 5;
elseif (theta<=12&&theta>=6)
	thetaBucket = 6;
end

if (x<-0.8&&x>=-2.4)
	xBucket = 1;
elseif (x<=0.8&&x>=-0.8)
	xBucket = 2;
elseif (x<=2.4&&x>0.8)
	xBucket = 3;
end

if (xDot<-0.5)
	xDotBucket = 1;
elseif (xDot>=-0.5&&xDot<=0.5)
	xDotBucket = 2;
else
	xDotBucket = 3;
end

if (thetaDot<-50)
	thetaDotBucket = 1;
elseif (thetaDot>=-50&&thetaDot<=50)
	thetaDotBucket = 2;
else
	thetaDotBucket = 3;
end

box = sub2ind([6,3,3,3],thetaBucket, thetaDotBucket,xBucket,xDotBucket);
end
return;


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nonRegressionTest.m                           |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : run all non regression test                   |
%|  `---'  |                                                              |
%+========================================================================+

clear all
close all
clc

% Gypsilab path
run('../addpathGypsilab.m')

% Mesh management
cd('meshManagement/')
run('nrtMshBitree.m')
run('nrtMshClean.m');
run('nrtMshCube.m');
run('nrtMshOctree.m')
run('nrtMshRead.m');
run('nrtMshRefine.m');
run('nrtMshSegment.m');
run('nrtMshSquare.m');
run('nrtMshTransfo.m');
run('nrtMshWrite.m');

% Domain quadrature
cd('../domainQuadrature/')
run('nrtDom1D.m');
run('nrtDomEdge.m');
run('nrtDomHalfSquare.m');
run('nrtDomIntSing2D.m');
run('nrtDomIntSing3D.m');
run('nrtDomND.m');
run('nrtDomTetra.m');
run('nrtDomTrace.m');
run('nrtDomTriangle.m');

% Finite element
cd('../finiteElement/')
run('nrtFemAllWithShuffle.m');
run('nrtFemContinuity1.m');
run('nrtFemContinuity2.m');
run('nrtFemConvergence2D.m');
run('nrtFemDirichletCube.m');
run('nrtFemDirichletDisk.m');
run('nrtFemDirichletSquare.m');
run('nrtFemDirichletString.m');
run('nrtFemElasticite2D')
run('nrtFemJunction.m');
run('nrtFemLaplace.m');
run('nrtFemOperators.m');
run('nrtFemRwgNed.m');
run('nrtFemWave1D.m');
run('nrtFemWave2D.m');

% Hierarchical matrices
cd('../hierarchicalMatrix')
run('nrtHmxAlgebra.m');
run('nrtHmxBEMConvergence.m');
run('nrtHmxBuilder.m');
run('nrtHmxBuilderFem.m');
run('nrtHmxCompareLU.m');
run('nrtHmxCompressorPartial.m');
run('nrtHmxCompressorsBox.m');
run('nrtHmxCompressorSingular');
run('nrtHmxCompressorTotal.m');
run('nrtHmxCriticalDimension.m');
run('nrtHmxCriticalForm.m');
run('nrtHmxLowrank.m');
run('nrtHmxNonSquare.m');

% Fast & Free memory Method
cd('../fastFreeMemory')
run('nrtFfmAlgebra.m');
run('nrtFfmBuilder.m');
run('nrtFfmBuilderFem.m');
run('nrtFfmHelmholtzBWdir.m');
run('nrtFfmHelmholtzBWneu.m');
run('nrtFfmMaxwellCFIEpec.m');

% Block matrix
cd('../blockMatrix')
run('nrtBmmAlgebra.m');
run('nrtBmmStokesConv.m');

% Scattering 2D
cd('../scattering2d')
run('nrtLaplace2dSDrad.m')
run('nrtHelmholtz2dSDrad.m')
run('nrtHmxHelmholtz2dS.m');
run('nrtHmxHelmholtz2dD.m');
run('nrtHmxHelmholtz2dDt.m');
run('nrtHmxHelmholtz2dH.m');
run('nrtHmxHelmholtz2dBWdir.m');
run('nrtHmxHelmholtz2dBWneu.m');

% Scattering 3D
cd('../scattering3d')
run('nrtHelmholtzCalderon.m')
run('nrtHelmholtzSDrad.m');
run('nrtHmxHelmholtzS.m');
run('nrtHmxHelmholtzD.m');
run('nrtHmxHelmholtzDxy.m');
run('nrtHmxHelmholtzDt.m');
run('nrtHmxHelmholtzH.m');
run('nrtHmxHelmholtzBWdir.m');
run('nrtHmxHelmholtzBWneu.m');
run('nrtHmxMaxwellT.m');
run('nrtHmxMaxwellNxK.m');
run('nrtHmxMaxwellCFIE.m');

% Fem-Bem dielectrique
cd('../femBemDielectrique')
run('nrtHmxFemBemEFIE.m');
run('nrtHmxFemBemEFIEhalf.m');
run('nrtHmxFemBemCFIE.m');
run('nrtHmxFemBemCFIEhalf.m');

% Inverse problem
cd('../inverseProblem')
run('nrtIpbHelmholtz.m');
run('nrtIpbHelmholtz0.m');

% Ray-tracing
cd('../rayTracing');
run('nrtRayCube.m');
run('nrtRayFabryPerot.m');
run('nrtRayLabyrinthe.m');
run('nrtRaySource.m');
run('nrtRaySphere.m');
run('nrtRayTheatre');

% Stokes
cd('../stokes');
run('nrtHmxStkConvergence.m');
run('nrtStkConvergence.m');
run('nrtStkRadiation.m');
run('translatingSphere/nrtStkTranslatingSphere.m');

% Vibro-acoustic
cd('../vibroAcoustic');
run('nrtVibroSlab2d.m');
run('nrtHmxVibroSlab2d.m');

% Operators
cd('../operators');
run('nrtOprValidation');

% End
disp('~~> Non Regresion Test are done. Michto gypsilab !');


text:
function [rd, boardName] = hdlcoder_ref_design_customization
% Reference design plugin registration file
% 1. The registration file with this name inside of a board plugin folder
%    will be picked up
% 2. Any registration file with this name on MATLAB path will also be picked up
% 3. The registration file returns a cell array pointing to the location of
%    the reference design plugins
% 4. The registration file also returns its associated board name
% 5. Reference design plugin must be a package folder accessible from
%    MATLAB path, and contains a reference design definition file

%   Copyright 2013-2014 The MathWorks, Inc.

rd = {...
    'AnalogDevices.fmcomms5.zc706.plugin_rd_rx', ...
    'AnalogDevices.fmcomms5.zc706.plugin_rd_tx', ...
    'AnalogDevices.fmcomms5.zc706.plugin_rd_rxtx', ...
    };

boardName = 'AnalogDevices FMCOMMS5 ZC706';

end

author: analogdevicesinc


text:
function test_ft_denoise_hfc(testfile)

% WALLTIME 00:10:00
% MEM 8gb
% DEPENDENCY test_ft_denoise_hfc

% switch nargin
%   case 0
%     locate_function = which('test_ft_denoise_hfc');
%     % load data
%     load(fullfile(locate_function,'..','..','..','project_zeus',...
%       'data_UCL_OPM.mat'));
%   case 1
%     % If user inputs a filename, use path specific for DCCN
%     % Find this data: https://doi.org/10.17605/OSF.IO/CJNXH
%     load(dccnpath('/home/common/matlab/fieldtrip/data/test/pull2123/data_UCL_OPM.mat'));
% end
if nargin == 0
  testfile = dccnpath('/home/common/matlab/fieldtrip/data/test/pull2123/data_UCL_OPM.mat');
else
  % testfile is required in the input, as a full path to a mat-file
end
load(testfile);

% Test of Harmonic Field Correction (HFC) where L = 1
cfg               = [];
cfg.order         = 1;
cfg.residualcheck = 'yes';
[data_out]        = ft_denoise_hfc(cfg,data);

% Check L = 2
cfg               = [];
cfg.order         = 2;
cfg.residualcheck = 'yes';
[data_out]        = ft_denoise_hfc(cfg,data);

% Check L = 3
cfg               = [];
cfg.order         = 3;
cfg.residualcheck = 'yes';
[data_out]        = ft_denoise_hfc(cfg,data);

% Check L = 1 without residual check
cfg               = [];
cfg.order         = 1;
cfg.residualcheck = 'no';
[data_out]        = ft_denoise_hfc(cfg,data);


text:
function DoManualShim

global VCtl
global VObj
global VMag

if ~isfield(VCtl,'Sh_X2_Y2')
    warndlg('Please load Shim tab. Manual Shimming was not performed!');
    return;
end

Mxdims=size(VObj.Rho);
[xgrid,ygrid,zgrid]=meshgrid((-(Mxdims(2)-1)/2)*VObj.XDimRes:VObj.XDimRes:((Mxdims(2)-1)/2)*VObj.XDimRes,...
                            (-(Mxdims(1)-1)/2)*VObj.YDimRes:VObj.YDimRes:((Mxdims(1)-1)/2)*VObj.YDimRes,...
                            (-(Mxdims(3)-1)/2)*VObj.ZDimRes:VObj.ZDimRes:((Mxdims(3)-1)/2)*VObj.ZDimRes);

B0ShimField = VCtl.Sh_X .* xgrid + ...
              VCtl.Sh_Y .* ygrid + ...
              VCtl.Sh_Z .* zgrid + ...
              VCtl.Sh_ZX.* zgrid .* xgrid + ...
              VCtl.Sh_ZY.* zgrid .* ygrid + ...
              VCtl.Sh_Z2.* zgrid .^ 2 + ...
              VCtl.Sh_XYZ.* xgrid .* ygrid .* zgrid+ ...
              VCtl.Sh_X2_Y2.* (xgrid .* ygrid) .^ 2;

VMag.dB0 = VMag.dB0 + B0ShimField;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [kPa] = hPa2kPa(hPa)
% Convert pressure from hectopascals to kilopascals.
% Chad Greene 2012
kPa = hPa/10;


text:
% Prabhakar. S
% AE 110 Lab Assignment
% Satellite Tracking Program
% Satellite:  NOAA 12
% Program Name elaz : Elevation-Azimuth being called by the program sts

% Orbital Elements ready for degree to radians conversion!
% Two Line Elements and the Epoch Data are obtained from the GUI Interface

load variables.mat; % Loading the variables from the GUI Interface sts
I2 = yoe;		% Year of epoch
J2 = moe;			% Month of epoch
%h2 = hoe;       % Hour of epoch
%m2 = mioe;       % Minute of epoch
%K2 = doe + hoe/24 + mioe/(24*60);		% Day of epoch
K2=doe;
zero=0; % For printing a single digit month or day preceeded by a zero.
%ut=.08235372;
%frac=.08235372; % fraction of the day from the TLE

% This calculation is to find the Julian Date,
% JD Fortran code written by Fliegel and Van Flandern [1968]
% Handout by Dr. P

JD = (367*I2 - (fix(7*(I2 + fix((J2 + 9)/12))/4)) + fix((275*J2)/9) + K2 + 1721013.5 + frac);

%JD = (K2 - 32075 + round(1461 * (I2 + 4800 + (J2 -14)/12)/4 ...
%+ 367* (J2 - 2 - (J2 - 14)/12*12)/12 - 3*((I2 + 4900 + (J2 - 14)/12)/100)/4)) ;
% fprintf('\n Two Line Elements Data - 28 Dec 2000, Time 11:30:00\n');
% fprintf(' The Julian Date is %f',JD);

% Propagation begins
I = yoe;								% Year
J = moe;								% Month
%h = hoe;   				   		% Hour
%m = mioe;     				 		% Minute
%K = doe + hoe/24 + mioe/(24*60);	% Day
K=doe;

JD0 = (367*I - (fix(7*(I + fix((J + 9)/12))/4)) + fix((275*J)/9) + K + 1721013.5 + frac);

%JD0 = (K - 32075 + round(1461 * (I + 4800 + (J -14)/12)/4 ...
%+ 367 * (J - 2 - (J - 14)/12*12)/12 - 3*((I + 4900 + (J - 14)/12)/100)/4)) ;

I3 = yoe;							% Year
J3 = moe;							% Month
%h3 = hoe;       					% Hour
%m3 = mioe;       				% Minute
%K3 = (doe+pro) + h3/24 + m3/(24*60);	% Day
K3 = (doe+pro);

JD1 = (367*I3 - (fix(7*(I3 + fix((J3 + 9)/12))/4)) + fix((275*J3)/9) + K3 + 1721013.5 + frac);

%JD1 = K3 - 32075 + round(1461 * (I3 + 4800 + (J3 -14)/12)/4 ...
%  + 367 * (J3 - 2 - (J3 - 14)/12*12)/12 - 3*((I3 + 4900 + (J3 - 14)/12)/100)/4);
% Propagation ends here.

% Propagation step Size
step = 5/(24*60); % right now for every 5 minutes
j=1;

fprintf('\n\t\t\t SATELLITE    TRACKING    SYSTEM  \n');
fprintf(' ---------------------------------------------------------------------------');
fprintf('\n				             G r e g o r i a n   D a t e \n');
fprintf(' ---------------------------------------------------------------------------');
fprintf('\n  Julian Date      Azimuth     Elevation  Month/Day/Year Hour:Minute:Seconds\n');
fprintf(' ---------------------------------------------------------------------------');

while JD0<JD1;

% Declaring variables and initializing values

degtorad = pi/180;		% Radians conversion factor
i = incl*degtorad;		% Inclination
CapOmega = raan*degtorad;		% Right Ascension of the Ascending Node - denoted by symbol Capital Omega
ecc = e;					% Eccentricity
SmallOmega = aop*degtorad;						% Argument of Periapsis - denoted by  small Omega
Mo = ma * degtorad;		% Initial Mean anomaly
n = revs;	   		% Mean Motion (revs/day)
GP = 3.986e5;				% Gravitational parameter - denoted by nu
a = 24000;					% Semi-major Axis
% a = (GP/n^2)^(1/3);
% fprintf('Semi-major axis is computed to be :\n',a);

Latitude = 0 * degtorad;
Longitude = 280 * degtorad;
L = ((90*degtorad) - Latitude);

We = 6.300388097;		% Angular velocity of planet Earth
perturb = 98.9246096453622 * degtorad; % Perturbation due to the Oblateness
tq = 2448621.5; % JD for
Qg = perturb + We * (JD0 - tq);
Q = Qg + Longitude;

tp = JD - Mo/n;
M = n * (JD0 - tp);

E = M;				          		% Use M for the first value of E
	f=M-E+ecc*sin(E);			 		% Kepler's Equation
	while abs(f)>0.00000001; 		% Loop for Kepler's equation begins
		f=M-E+ecc*sin(E);			 	%
		fd=-1+ecc*cos(E);			 	%
		E=E-f/fd;						%
	end;							 		% loop for Kepler's eqn ends here!

% Computing true anomaly and radius

TA=2*atan((((1+ecc)/(1-ecc))^0.5)*tan(E/2));
    if TA<0;
    	  TA=2*pi+TA;
    end;

r=a*(1-ecc*cos(E)); % radius

% Coordinate transformations

adot = (-(3/2) * n * 1.0827e-3 * (6378.14/a)^2 * cos(i))/(1-ecc^2)^2;
A = CapOmega + adot * (JD0 - JD);
SmallOmega = SmallOmega + adot * (JD0 - JD);

Rseu = [0; 0; 6378.14];
Txyz = [ cos(A) -sin(A) 0; sin(A) cos(A) 0; 0 0 1 ] * [ 1 0 0; 0 cos(i) -sin(i); 0 sin(i) cos(i)] * ...
   [ cos(SmallOmega) -sin(SmallOmega) 0; sin(SmallOmega) cos(SmallOmega) 0; 0 0 1 ];
Tseu = [ cos(L) 0 -sin(L); 0 1 0; sin(L) 0 cos(L)] * [ cos(Q) sin(Q) 0; -sin(Q) cos(Q) 0; 0 0 1];
ruvw = [ (r * cos(TA)); (r * sin(TA)); 0 ];
Pseu = Tseu*(Txyz*ruvw) - Rseu;
P = (Pseu(1)^2 + Pseu(2)^2 + Pseu(3)^2)^0.5;

El = (asin ( Pseu(3)/P ))/degtorad;
Az = (atan2 ( Pseu(2), -Pseu(1)))/degtorad;
	if Az < 0;
   	Az = 360 + Az;
   end;

%Loop begins here  %if El>0;

array(j,1) = JD0;
array(j,3) = Az;
array(j,4) = El;
%array(j,5) = cat(year,month,day);
j = j+1;
% loop ends here

JD0 = JD0 + step;
			% Gregorian Date from Julian Date
			% This calculation is valid for any Julian Day Number including negative JDN and
			% produces a Gregorian date (or possibly a proleptic Gregorian date)

			Z = floor(JD0 - 1721118.5); % JD0 is the Julian Date in propagation
			R = (JD0 - 1721118.5 - Z);  % R is the fractional part of JD0
			G = (Z - .25);
			A = (floor(G / 36524.25));  % Calculate the value of A which is the number of full centuries
			B = (A - (A / 4)); % The value of B is this number of days minus a constant
			year = (floor((B+G) / 365.25)); % Calculate the value of Y, the year in a calendar whose years start on March 1
			C = (B + Z - floor(365.25 * year));  % Day count
         month = (fix((5 * C + 456) / 153));  % Month
         UT = (C - fix((153 * month - 457) / 5) + R); % Calculation for UTC
         day = floor(UT);  % Gregorian Day

						if month > 12;
		      			year = year + 1 ;
		      			month = month - 12;
                  end;
                  UT = UT - floor(UT);
         UT = UT*24;
         hr = floor(UT);   %  hour
         UT = UT-floor(UT);
	      UT =UT* 60;
	      min = floor(UT);  % minute
   	   UT =UT- floor(UT);
      	UT =UT* 60;
         secs = round(UT); % seconds
         if secs==60;
            min=min+1;
            secs=0;
         end;
         % Gregorian Date conversion algorithm ends here

         % computations for printintg two digits in case year,month,day,min and secs is < 10
         % (i.e) to print third month, 3 as 03 and day 5 as 05 etc..
         ZERO=int2str(zero);
			MIN=int2str(min); % converting minute into a string
			if min <10;
		  			 MinNew=strcat(ZERO,MIN);
                  min=MinNew;
               else;
                  min=MIN;
               end;

         	DAY=int2str(day); % converting DAY into a string
				 if day <10;
		  				 DayNew=strcat(ZERO,DAY);
                   day=DayNew;
                else;
                  day=DAY;
         	 end;
          HR=int2str(hr); % converting Hour into a string
					if hr <10;
		  				 HrNew=strcat(ZERO,HR);
                   hr=HrNew;
               else;
                 	 hr=HR;
		          end;

                MONTH=int2str(month); % converting month into a string
						if month <10;
			  				 MonthNew=strcat(ZERO,MONTH);
   	                month=MonthNew;
            	   else
         	         month=MONTH;
                  end;

                 SECS=int2str(secs); % converting seconds into a string
						if secs <10;
			  				 SecsNew=strcat(ZERO,SECS);
   	                secs=SecsNew;
            	   else
         	         secs=SECS;
			         end;

         % computations for printing two digits ends here!

         if El>0;
            fprintf('\n %f | %f | %f |  %s/%s/%4d  |  %s:%s:%s  |',JD0, Az, El, month,day,year, hr,min,secs);
            %plot(El,Az);
         end;
end;
fprintf('\n ---------------------------------------------------------------------------\n');
fprintf('\n\t\t\t       End of Output\!!! \n\n\n');


text:
function Y = as_m(tr,X)
%
%  Evaluates matrix-matrix products with the real, symmetric matrix A,
%  i.e., Y = A*X.
%
%  A is provided as global data. This data must be generated by calling
%  'as_m_i' before calling this routine!
%
%  If called without input parameters, this routine returns the order of
%  the matrix A.
%
%  Calling sequence:
%
%    Y = as_m(tr,X)
%    n = as_m
%
%  Input:
%
%    tr        is not referenced;
%    X         a matrix of proper size.
%
%  Output:
%
%    Y         the resulting product;
%    n         the order of the matrix A.
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

ni = nargin;

if ni~=2 & ni~=0
  error('Wrong number of input arguments.');
end

global LP_A

if ~length(LP_A)
  error('This routine needs global data which must be generated by calling ''as_m_i'' first.');
end

if ni==0
  Y = size(LP_A,1);
else
  Y = LP_A*X;
end


text:
function matlab_vectors = hypreParVectors2matlab(filename,num_of_vectors)
%HYPREPARVECTORS2MATLAB Converts Hypre IJ Vectors to Matlab.
%
%   This program was tested on Matlab 7.1.
%
%   Description:
%     [matlab_vectors] = hypreParVectors2matlab(filename, num_of_vectors)
%     converts hypre formatted vectors, specified by the argument 'filename'
%     to matlab formatted vectors. Number of vectors to be converted is
%     given by the argument num_of_vectors
%
%   Inputs:
%     filename = name of files without suffix
%     num_of_vectors = number of files containing vector description
%
%    Example:
%      [matlab_vectors] = hypreParVectors2matlab( 'vectors', 5 ) converts
%      HYPRE vector represented in 5 input HYPRE files to MATLAB format.
%      HYPRE formatted vectors (filename is 'vectors', for 2 CPUs) are
%        represented as follows:
%          vectors.0.0  vectors.0.1  vectors.INFO.0.0  vectors.INFO.0.1
%          vectors.1.0  vectors.1.1  vectors.INFO.1.0  vectors.INFO.1.1
%          vectors.2.0  vectors.2.1  vectors.INFO.2.0  vectors.INFO.2.1
%          ...
%      The first line in vectors.0.0 and vectors.0.1 represents the
%        number of elements in that particular file
%        vectors.INFO.0.0 and vectors.INFO.0.1 contain information
%        about the size of the vector and the size of the particular file
%      For these two HYPRE formatted vectors
%        vectors.0.0  vectors.0.1  vectors.INFO.0.0  vectors.INFO.0.1
%        vectors.1.0  vectors.1.1  vectors.INFO.1.0  vectors.INFO.1.1
%        vectors.2.0  vectors.2.1  vectors.INFO.2.0  vectors.INFO.2.1
%        ...
%        the program will create a (size_of_vector)-by-2 matrices in matlab.
%
%   See also matlabIJ2hypre, matlab2hypreParVectors.m,
%      hypreIJ2matlab.m, testIJmatlabhypre.m
%
%   Author: Diana Zakaryan, Dept of Mathematics, University of Colorado,
%      Denver, 15-Mar-2005.
%

% empty vector
matlab_vectors = [];

% calls the actual program num_of_vectors time
% in order to convert all num_of_vectors hypre vectors to matlab format
for j=1:num_of_vectors
    s=int2str(j-1);
    filename2=strcat(filename,'.',s);
    % calls single vector convertor
    matlab_vector = hypreParVector2matlab(filename2);
    matlab_vectors = cat(2,matlab_vectors,matlab_vector);
end

function matlab_vector = hypreParVector2matlab(filename2);

%read all the attributes of the specified filename
%in the specified directory
[stat,mess]=fileattrib(strcat(filename2,'.*'));
%[pathstr,name,ext,versn]= fileparts(mess(1).Name);
%empty vector
matlab_vector=[];
for i=1:size(mess,2)
    [pathstr,name,extension,versn]= fileparts(mess(i).Name);
    k = findstr('.INFO',name);
    if size(k) == 1
        continue;
    end
    % fills in with data
    filename_temp=strcat(name,extension);
    hypre_data = dlmread(filename_temp,'',1,0);
    if i==1
        matlab_vector = hypre_data;
    else
        matlab_vector=cat(1,matlab_vector,hypre_data);
    end
    clear ('hypre_data');
end


text:
function i4_is_triangular_test ( )

%*****************************************************************************80
%
%% I4_IS_TRIANGULAR_TEST tests I4_IS_TRIANGULAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_IS_TRIANGULAR_TEST\n' );
  fprintf ( 1, '  I4_IS_TRIANGULAR returns T or F depending on\n' );
  fprintf ( 1, '  whether I is triangular.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I   T/F\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 20

    l = i4_is_triangular ( i );

    fprintf ( 1, '  %2d    %1d\n',  i, l );

  end

  return
end


text:
function X = sampleimages( winsize, samples )
% sampleimages - gathers image patches from Olshausens images.
%

% Start by loading the images
fprintf('Loading images...\n');
load ../data/IMAGES;
num_images=size(IMAGES,2);
image_size=sqrt(size(IMAGES,1));

% This will hold the patches
X=zeros(winsize^2,samples);
totalsamples = 0;

for i=1:num_images

  % Choose an image for this batch
  this_image=reshape(IMAGES(:,i),image_size,image_size)';
  BUFF=4;

  % Determine how many patches to take
  getsample = floor(samples/num_images);
  if i==num_images, getsample = samples-totalsamples; end

  % Extract patches at random from this image to make data vector X
  for j=1:getsample
    r=BUFF+ceil((image_size-winsize-2*BUFF)*rand);
    c=BUFF+ceil((image_size-winsize-2*BUFF)*rand);
    totalsamples = totalsamples + 1;
    X(:,totalsamples) = ...
	reshape( this_image(r:r+winsize-1,c:c+winsize-1),winsize^2,1);
  end

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function hermite_cubic_test13 ( )

%*****************************************************************************80
%
%% HERMITE_CUBIC_TEST13 tests HERMITE_CUBIC_LAGRANGE_INTEGRATE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_CUBIC_TEST13:\n' );
  fprintf ( 1, '  HERMITE_CUBIC_LAGRANGE_INTEGRATE integrates a Hermite cubic\n' );
  fprintf ( 1, '  Lagrange polynomial from A to B.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute each result TWICE:\n' );
  fprintf ( 1, '  First row computed using HERMITE_CUBIC_INTEGRATE.\n' );
  fprintf ( 1, '  Second row computed using HERMITE_CUBIC_LAGRANGE_INTEGRATE.\n' );

  x1 = 0.0;
  x2 = 10.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '        A           B           LF1         LD1         LF2         LD2\n' );
  fprintf ( 1, '\n' );

  a = x1 - 1.0;

  for j = - 3 : 12

    b = ( ( 10 - j ) * x1   ...
        +        j   * x2 ) ...
        /   10.0;

    p(1) = hermite_cubic_integrate ( x1, 1.0, 0.0, x2, 0.0, 0.0, a, b );
    p(2) = hermite_cubic_integrate ( x1, 0.0, 1.0, x2, 0.0, 0.0, a, b );
    p(3) = hermite_cubic_integrate ( x1, 0.0, 0.0, x2, 1.0, 0.0, a, b );
    p(4) = hermite_cubic_integrate ( x1, 0.0, 0.0, x2, 0.0, 1.0, a, b );

    q = hermite_cubic_lagrange_integrate ( x1, x2, a, b );

    fprintf ( 1, '  %10.4f  %10.4f  %10.4f  %10.4f  %10.4f  %10.4f\n', ...
      a, b, p(1:4) );
    fprintf ( 1, '                          %10.4f  %10.4f  %10.4f  %10.4f\n', ...
      q(1:4) );

  end

  return
end


text:
% Calculate head direction.
%
% Calculates the head direction for each position sample pair. Direction
% is defined as east = 0 degrees, north = 90 degrees, west = 180 degrees,
% south = 270 degrees. Direction is set to NaN for missing samples.
% Position matrix contains information about front and back LED. Head
% direction is the counter-clockwise direction from back LED to the front.
%
%  USAGE
%    hd = calcHeadDirection(positions)
%
%    positions          Animal's position data, Nx5. Position data should
%                       contain timestamps (1 column), X/Y coordinates of
%                       first LED (2 and 3 columns correspondingly), X/Y
%                       coordinates of the second LED (4 and 5 columns
%                       correspondingly).
%                       it is assumed that positions(:, 2:3) correspond to
%                       front LED, and positions(:, 4:5) to the back LED.
%                       The resulting hd is the directon from back LED to
%                       the front LED.
%    hd                 Vector of head directions in degrees.
%
function hd = calcHeadDirection(positions)

    if size(positions, 2) < 5
        error('Position data should be 2D (type ''help analyses.calcHeadDirection'' for details).');
    end
    % 1 column contains timestamps
    x1 = positions(:, 2);
    y1 = positions(:, 3);
    x2 = positions(:, 4);
    y2 = positions(:, 5);

    hd = rem(atan2d(y2-y1, x2-x1) + 180, 360);
end

author: kavli-ntnu


text:
function [X Z G U histo histo_R] = SPC(T,Q,TV_SV,rho,K,SNR,rate,maxiter,tol,out)

  warning off;

  histo=[];
  histo_R=[];
  N  = ndims(T);
  II = size(T);
  NN = sum(Q(:));
  epsiron = 10^(-SNR/10) * sum_all(T(Q).^2);

  %% initialization
  E  = zeros(II);
  X  = zeros(II);
  Z  = zeros(II);
  ONE= zeros(II);

  R  = 1;

  X(Q) = T(Q);
  X(~Q) = sum(X(:))/NN;
  for n = 1:N
    U{n} = randn(II(n),R);
    U{n} = U{n}/norm(U{n});
    Pp  = eye(II(n)-1,II(n));
    Pm  = eye(II(n)); Pm(1,:) = [];
    P{n} = Pp - Pm;
    PtP{n} = P{n}'*P{n};
    Pr{n} = inv(rho(n)*PtP{n} + eye(II(n)));
  end
  G = tensor_allprod(X,U,1);

  for r = 1:R
    Z = Z + G(r)*outerprod(U,r);
  end
  obj = sum_all((T(Q) - Z(Q)).^2);
  X(~Q) = Z(~Q);
  E = X - Z;

  %% output some figures

  h1 = figure();clf;
  subplot(1,2,1);cla;hold on;
  subplot(1,2,2);
  drawnow;
  if out == 1
    h2 = figure();clf;
    imagesc(uint8(X));drawnow;
    imwrite(uint8(X),['saved/' TV_SV '_iter_0.png']);
  end

  %% start main algorithm
  for iter = 1:maxiter

    [val ID] = sort(G);
    for k = ID(1:min(end,K))

      ONE = G(k)*outerprod(U,k);
      Z = Z - ONE;
      E = E + ONE;

      div = 1;
      for n = 1:N

        if strcmp(TV_SV,'tv')

          u = innerprod_one_exc(E,U,k,n);
          u = u(:);
          a = u(:);
          object = 0.5*G(k)*sum(abs(P{n}*a)) + a'*u + 0.5*G(k);
          for nn = 1:1000
            df = P{n}'*sign(P{n}*a);
            dL = (0.5*G(k)*rho(n)*df - u + G(k)*a)/G(k);
            al = [0 0.1 0.01 0.001 0.0001 0.00001];
            for ai = 1:length(al)
              a2 = a - al(ai)*dL;
              score(ai) = 0.5*G(k)*sum(abs(P{n}*a2)) + a2'*u + 0.5*G(k);
            end
            [object2 ai] = min(score);
            a2 = a - al(ai)*dL;
            if abs(object2 - object)/II(n) < 1e-3
              break;
            else
              a = a2;
              object = object2;
            end
          end
          lam = norm(a);
          u = a/lam;
          U{n}(:,k) = u;
          v{n} = u;
          div  = div + rho(n)*sum(abs(P{n}*u));

        elseif strcmp(TV_SV,'sv')

          u = innerprod_one_exc(E,U,k,n);
          u = Pr{n} * u(:);
          lam = norm(u);
          u = u/lam;
          U{n}(:,k) = u;
          v{n} = u;
          div  = div + rho(n)*u'*PtP{n}*u;

        else
          error('2rd input is ''tv'' or ''qc'' ');
        end

      end
      G(k) = tensor_allprod(E,v,1)/div;
      if G(k) < 0
        G(k) = -G(k);
        U{1}(:,k) = - U{1}(:,k);
      end

      ONE = G(k)*outerprod(U,k);
      E = E - ONE;
      E(~Q) = 0;
      Z = Z + ONE;
      X(~Q) = Z(~Q);

    end

    %% calculate MSE
    obj2 = sum_all(E.^2);

    %% convergence speed
    speed = abs(obj2 - obj)/abs(epsiron - obj2);

    %% step for R <= R + 1
    if speed < rate || abs(obj2-obj)/NN < tol

      R = R + 1;
      for n = 1:N
        u = randn(II(n),1);
        u = u/norm(u);
        U{n}(:,R) = u;
        v{n} = u;
      end
      G(R) = tensor_allprod(E,v,1);
      E = E - G(R)*outerprod(U,R);
      E(~Q) = 0;

    end

    %% checking convergence
    if obj2 < epsiron
      break;
    else
      obj = obj2;
      if mod(iter,5)==0
        fprintf('%d:  %f :: %f :: %f :: Pid %d \n',iter,obj2/NN,epsiron/NN,speed,R);
      end
      histo(iter) = obj;
      histo_R(iter,:) = size(G);
    end

    %% output figures
    set(0,'CurrentFigure',h1);
    subplot(1,2,1);cla;hold on;
    plot((histo));
    plot((epsiron)*ones(1,length(histo)))
    grid on;
    set(gca,'YScale','log')
    title('MSE')
    subplot(1,2,2);
    plot(histo_R(:,2));
    title('number of components R')

    drawnow;
    if out == 1
      set(0,'CurrentFigure',h2);
      imagesc(uint8(X));
      title(['Number of R = ' num2str(R)]);
      drawnow;
    end
    if mod(iter,10) == 0 & out == 1
      imwrite(uint8(Z),['saved/' TV_SV '_iter_' num2str(iter) '.png']);
    end

    if mod(iter,100) == 0
      pack;
    end

  end

author: andrewssobral


text:
function [points3d, covariance] = myEsimatePosAmers(pointTracks, ...
    camPoses, cameraParams,ParamFilter)
chiC = ParamFilter.chiC;
RotC = chiC(1:3,1:3);
xC = chiC(1:3,4);
R = eye(2);
k = 2*length(camPoses);
q = 6*length(camPoses);
N_aug = q+k;
Rc = chol(kron(eye(k/2),R));

S = zeros(q);
chi = cell(k/2,1);
for i = 1:k/2
    S(6*i-5:6*i,6*i-5:6*i) = camPoses(i).S([1:3 7:9],[1:3 7:9]);
    chi{i} = [camPoses(i).Orientation camPoses(i).Location;0 0 0 1];
end

S_aug = blkdiag(S,Rc);

% scaled unsented transform
W0 = 1-N_aug/3;
Wj = (1-W0)/(2*N_aug);
gamma = sqrt(N_aug/(1-W0));
alpha = 1;
beta = 2;

% Compute transformed measurement
X = zeros(N_aug,2*N_aug+1); %sigma points
Y = zeros(3,2*N_aug+1);
ybar = zeros(3,1); % Measurement mean
for j = 1:2*N_aug+1
    if j == 1
    elseif j < N_aug+2
        X(:,j) = gamma*S_aug(j-1,:)';
    else
        X(:,j) = -gamma*S_aug(j-N_aug-1,:)';
    end
    camPoses_j = camPoses;
    pointTracks_j = pointTracks;
    xi_j = X(1:q,j);
    v_j = X(q+1:N_aug,j);
    for i = 1:k/2
        chi_j = expSE3(xi_j(6*i-5:6*i))*chi{i};
        camPoses_j(i).Orientation = RotC'*chi_j(1:3,1:3)';
        camPoses_j(i).Location = chi_j(1:3,4) +chi_j(1:3,1:3)*RotC*xC;
        pointTracks_j.Points(i,:) = pointTracks.Points(i,:) + v_j(2*i-1:2*i)';
    end

    Y(:,j) = EsimatePosAmers(pointTracks_j,camPoses_j,cameraParams)';
    if j == 1
        ybar = ybar + W0*Y(:,j);
    else
        ybar = ybar + Wj*Y(:,j);
    end
end

Y(:,1) = sqrt(abs(W0+(1-alpha^2+beta)))*(Y(:,1)-ybar);
YY = sqrt(Wj)*(Y(:,2:2*N_aug+1)-ybar*ones(1,2*N_aug));
[~,Rs] = qr(YY');
Ss = Rs(1:3,1:3);
[Sy,~] = cholupdate(Ss,Y(:,1),'-'); % Sy'*Sy = Pyy

Pxy = zeros(q,3);
for j = 2:2*N_aug+1
    Pxy = Pxy + Wj*X(1:q,j)*(Y(:,j)-ybar)';
end

points3d = ybar;
S = S(end-5:end,end-5:end);
covariance = [S'*S Pxy(end-5:end,:); Pxy(end-5:end,:)' Sy'*Sy];
end

author: mbrossar


text:
function [totSec hh mm ss fract] = dcm_hhmmss(dateStr)
if ~ischar(dateStr)
    dateStr = num2str(dateStr);
end
hh = str2double(dateStr(1,1:2));
mm = str2double(dateStr(1,3:4));
ss = str2double(dateStr(1,5:6));
fract_start = find(dateStr,'.');
fract = str2double(dateStr(1,fract_start+1:end));
totSec = hh*60*60 + mm*60 + ss;
end

author: mvallieres

repo: radiomics


text:
%Transcriptional Regulation
%
% @wholeCellModelID Process_TranscriptionalRegulation
% @name             Transcriptional Regulation
% @description
%   Biology
%   ==============
%   The the rate of transcription of each transcription unit is known to be
%   regulated by proteins referred to as transcription factors which modulate
%   the affinity of the RNA polymerase for each promoter. Transcription factors
%   can both positively and negatively modulate RNA polymerase - promoter
%   affinity. Transcription factors can stabilize RNA polymerase - promoter
%   complexes by contributing additional negative free energy to the complex by
%   providing an additional surface for the RNA polymerase binding.
%   Transcription factors can destabilize RNA polymerase - promoter complexes
%   for example by sterically blocking promoters and preventing the RNA
%   polymerase from binding promoters.
%
%   This process simulates the binding of transcription factors to the promoters
%   of each transcription unit (affinity), as well as the affect of the
%   transcription factors on the recruiment of the RNA polymerase to each
%   promoter (activity). The process was built using experimentally observed fold
%   change expression effect of each transcription factor on each promoter.
%
%   Affinity
%   ++++++++++++++
%   Binds enzymes to promoters assuming:
%   1) Transcription factors have high affinity for promoters
%   2) Transcription factors bind promoters rapidly
%   3) Transcription factors bind promoters stably over time
%   4) Only 1 copy of a transcription factor can bind each promoter
%   5) Transcription factors only compete within their species for
%      promoters. That is for each transcription unit we assume each
%      transcription factor species binds a distinct promoters region.
%
%   Consequently, at each time step we simulate that each free
%   transcription factor binds randomly binds unoccupied promoters (no copy
%   of that transcription factor is already bound to the promoter). Random
%   transcription factor-promoter binding is weighted by the affinity of
%   each transcription factor for each promoter.
%
%   Because transcription factor-promoter affinities are generally not
%   experimentally observed, we base them on the transcription factor fold
%   change activities.
%
%   Activity
%   ++++++++++++++
%   The effect of bound transcription factors on the recruitment of RNA
%   polymerase and the expression of transcription units is simulated here
%   and incorporated into the calculation of the RNA polymerase
%   transcription unit promoter binding probabilities in the transcription
%   process. Specifically, the wild-type average RNA polymerase
%   transcription unit promoter binding probabilities are multiplied by the
%   binding probability fold change effects simulated in this process.
%
%   The RNA polymerase binding probability fold change is simulated for
%   each promoter as the product of the observed expression fold change
%   effects of each bound transcription factor. When a promoter is bound by
%   a single transcription factor, the net RNA polymerase binding
%   probability fold change is the observed expression fold change of that
%   transcription factor. When a promoter is bound by multiple
%   transcription units, the net RNA polymerase binding probability fold
%   change is given by the product of the individual fold change effects of
%   the bound transcription factors.
%
%   Knowledge Base
%   ==============
%   The list of transcriptional regulatory relationships is maintained in the
%   the knowledge base. As of 8/10/2010, it contained 31 such relationships
%   between 5 transcription factors and 29 transcription units containing 37
%   genes. The knowledge base was built from a variety of literature sources
%   and databases including:
%   - PUB_0096
%   - PUB_0110
%   - PUB_0112
%   - PUB_0196
%   - PUB_0418-20
%   - PUB_0433-8
%   - PUB_0505
%
%   Initialization
%   =================
%   Because we assume that transcription factors have high affinity for DNA and
%   bind DNA stably, we initialize as many transcription factors as possible to
%   the promoter-bound state. We randomly assign transcription factors to
%   promoters using their relative affinities.
%
%   Simulation
%   ==============
%   For each kind of transcription factor, bind any free transcription factors
%   to promoters that aren't already occupied by this kind of transcription
%   factor. Choose the binding sites randomly, weighted by the transcription
%   factor's affinity to them.
%
%   References
%   ==============
%   1) Lacramioara Bintu and Nicolas E Buchler and Hernan G Garcia and
%      Ulrich Gerland and Terence Hwa, Jane Kondev and Rob Phillips (2005).
%      Transcriptional regulation by the numbers: models. Curr Opin Genet
%      Dev. 15:116-24.
%   2) Lacramioara Bintu and Nicolas E Buchler and Hernan G Garcia and
%      Ulrich Gerland and Terence Hwa, Jane Kondev and Rob Phillips (2005).
%      Transcriptional regulation by the numbers: applications. Curr Opin
%      Genet Dev. 15:125-35.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 10/19/2010
classdef TranscriptionalRegulation < edu.stanford.covert.cell.sim.Process & edu.stanford.covert.cell.sim.ChromosomeProcessAspect
    %property annotations
    properties (Constant)
        optionNames__              = {}; %names of option properties
        fixedConstantNames__       = {   %names of fixed constant properties
		    'tfIndexs';
            'tuIndexs';
            'tfAffinities';
            'tfActivities';
            'tfPositionStrands';
            'otherActivities';
		    };
        fittedConstantNames__      = {}; %names of fitted constant properties
        localStateNames__          = {}; %names of simulation state properties redundant with timecourses in this or other processes or the simulation
    end

    %IDs, names, and local indices
    properties
        stimuliWholeCellModelIDs = {};
        substrateWholeCellModelIDs = {};

        enzymeWholeCellModelIDs
        enzymeIndexs_fur         %MG_236_MONOMER  ferric uptake repressor
        enzymeIndexs_gntR        %MG_101_MONOMER  Uncharacterized HTH-type transcriptional regulator
        enzymeIndexs_hrcA        %MG_205_DIMER    heat-inducible transcription repressor HrcA, putative
        enzymeIndexs_luxR        %MG_428_DIMER    LuxR bacterial regulatory protein, putative
        enzymeIndexs_spx         %MG_127_MONOMER  Spx subfamily protein

        transcriptionUnitWholeCellModelIDs
    end

    %fixed biological constants
    properties
        tfIndexs           %identifying index of transcription factor that binds each site [sites X 2]
        tuIndexs           %identifying index of transcription unit associated with each TF binding site [sites X 2]
        tfAffinities       %affinity of each transcription factor for each promoter it binds [sites X 2]
        tfActivities       %fold change expression of each transcription unit when its promoter is bound by a particular transcription factor [sites X 2]
        tfPositionStrands  %start position and strand of each transcription factor DNA binding site [2*sites X 2]
        otherActivities    %fold change each non-DNA-binding transcription factor has on each transcription unit's expression [nTFs X nTUs]
    end

    %dependent local state (implemented as dependent property for convenience)
    properties (Dependent, SetAccess = protected)
        tfBoundPromoters       %number of each transcription factor bound to promoter of each transcription unit [transcription factors X promoters X 2]
        boundTFs               %count of how many of each transcription factor are bound [transcription factors X 1]
    end

    %references to cell state
    properties
        rnaPolymerase             %reference to RNA polymerase state
    end

    %constructor
    methods
        function this = TranscriptionalRegulation(wholeCellModelID, name)
            this = this@edu.stanford.covert.cell.sim.Process(wholeCellModelID, name);
        end
    end

    %communication between process/simulation
    methods
        %set references to state objects
        function storeObjectReferences(this, simulation)
            this.storeObjectReferences@edu.stanford.covert.cell.sim.Process(simulation);
            this.storeObjectReferences@edu.stanford.covert.cell.sim.ChromosomeProcessAspect(simulation);

            this.rnaPolymerase = simulation.state('RNAPolymerase');
            this.states = [this.states; {this.rnaPolymerase}];
        end

        %initialize constants
        function initializeConstants(this, knowledgeBase, varargin)
            %build object arrays of transcription factor monomers and complexes
            monomers = [];
            complexs = [];
            for i = 1:knowledgeBase.numTranscriptionUnits
                m = knowledgeBase.transcriptionUnits(i).transcriptionFactorProteinMonomers;
                c = knowledgeBase.transcriptionUnits(i).transcriptionFactorProteinComplexs;
                if ~isempty(m)
                    monomers = [monomers; m]; %#ok<AGROW>
                end
                if ~isempty(c)
                    complexs = [complexs; c]; %#ok<AGROW>
                end
            end

            %whole cell model ids of all transcription factors
            this.enzymeWholeCellModelIDs = unique({
                monomers.wholeCellModelID...
                complexs.wholeCellModelID}');

            %transcription factors indices
            this.enzymeIndexs_fur  = this.enzymeIndexs({'MG_236_MONOMER'});
            this.enzymeIndexs_gntR = this.enzymeIndexs({'MG_101_MONOMER'});
            this.enzymeIndexs_hrcA = this.enzymeIndexs({'MG_205_DIMER'});
            this.enzymeIndexs_luxR = this.enzymeIndexs({'MG_428_DIMER'});
            this.enzymeIndexs_spx  = this.enzymeIndexs({'MG_127_MONOMER'});

            %call super class method to compute mapping between enzymes of
            %this process, and of protein monomers and complexes of the simulation
            this.initializeConstants@edu.stanford.covert.cell.sim.Process(...
                knowledgeBase, varargin{:});
            this.initializeConstants@edu.stanford.covert.cell.sim.ChromosomeProcessAspect(...
                knowledgeBase, varargin{:});

            %whole cell model ids of transcription units
            this.transcriptionUnitWholeCellModelIDs = ...
                {knowledgeBase.transcriptionUnits.wholeCellModelID}';

            %build matrices of transcription factor positions, affinities, activities
            n = length(monomers) + length(complexs);
            this.tuIndexs = zeros(n, 2);
            this.tfIndexs = zeros(n, 2);
            this.tfActivities = zeros(n, 2);
            this.tfPositionStrands = zeros(n*2, 2);
            c = this.chromosome;
            j = 1;
            for i = 1:knowledgeBase.numTranscriptionUnits
                tu = knowledgeBase.transcriptionUnits(i);

                if ~isempty(tu.transcriptionFactorProteinMonomers)
                    [tf, proteinIndexs] = ismember(...
                        [tu.transcriptionFactorProteinMonomers.idx]',...
                        this.enzymeMonomerGlobalIndexs);
                    k = nnz(tf) - 1;
                    this.tuIndexs(j:j+k) = i;

                    proteinIndexs = proteinIndexs(tf);
                    this.tfIndexs(j:j+k) = this.enzymeMonomerLocalIndexs(proteinIndexs);

                    regulationIndexs = ...
                        this.enzymeMonomerCompartmentIndexs(proteinIndexs) == ...
                        [tu.transcriptionFactorProteinMonomerCompartments.idx]';
                    this.tfActivities(j:j+k) = ...
                        tu.transcriptionFactorProteinMonomerActivitys(regulationIndexs);

                    footprints = c.monomerDNAFootprints([tu.transcriptionFactorProteinMonomers.idx]');
                    startCoordinates = tu.transcriptionFactorProteinMonomerBindingSiteStartCoordinates;
                    lengths = tu.transcriptionFactorProteinMonomerBindingSiteLengths;
                    this.tfPositionStrands(j:j+k) = ...
                        mod(ceil(startCoordinates + lengths/2 - footprints/2) - 1, c.sequenceLen) + 1;
                    j = j + 1 + k;
                end

                if ~isempty(tu.transcriptionFactorProteinComplexs)
                    [tf, proteinIndexs] = ismember(...
                        [tu.transcriptionFactorProteinComplexs.idx]',...
                        this.enzymeComplexGlobalIndexs);
                    k = nnz(tf) - 1;
                    this.tuIndexs(j:j+k) = i;

                    proteinIndexs = proteinIndexs(tf);
                    this.tfIndexs(j:j+k) = this.enzymeComplexLocalIndexs(proteinIndexs);

                    regulationIndexs = ...
                        this.enzymeComplexCompartmentIndexs(proteinIndexs) == ...
                        [tu.transcriptionFactorProteinComplexCompartments.idx]';
                    this.tfActivities(j:j+k) = ...
                        tu.transcriptionFactorProteinComplexActivitys(regulationIndexs);

                    footprints = c.complexDNAFootprints([tu.transcriptionFactorProteinComplexs.idx]');
                    startCoordinates = tu.transcriptionFactorProteinComplexBindingSiteStartCoordinates;
                    lengths = tu.transcriptionFactorProteinComplexBindingSiteLengths;
                    this.tfPositionStrands(j:j+k) = ...
                        mod(ceil(startCoordinates + lengths/2 - footprints/2) - 1, c.sequenceLen) + 1;
                    j = j + 1 + k;
                end
            end

            validateattributes(this.tfActivities, {'numeric'}, {'nonnegative'});

            this.tuIndexs(:,2) = this.tuIndexs(:,1);
            this.tfIndexs(:,2) = this.tfIndexs(:,1);
            this.tfActivities(:,2) = this.tfActivities(:,1);

            this.tfPositionStrands(1:n,2) = 1;
            this.tfPositionStrands(n+1:2*n,1) = this.tfPositionStrands(1:n,1);
            this.tfPositionStrands(n+1:2*n,2) = 3;

            this.tfAffinities = this.tfActivities;
            this.tfActivities(this.tfActivities<1 & this.tfAffinities>0) = ...
                this.tfActivities(this.tfActivities<1 & this.tfAffinities>0) .^ -1;

            bind = ~isnan(this.tfPositionStrands(1:n,1));
            this.otherActivities = ones(...
                length(this.transcriptionUnitWholeCellModelIDs), ...
                length(this.enzymeWholeCellModelIDs));
            this.otherActivities(sub2ind(...
                size(this.otherActivities), ...
                this.tuIndexs(~bind), ...
                this.tfIndexs(~bind))) = ...
                this.tfActivities(~bind);
            this.tuIndexs = this.tuIndexs(bind,:);
            this.tfIndexs = this.tfIndexs(bind,:);
            this.tfAffinities = this.tfAffinities(bind,:);
            this.tfActivities = this.tfActivities(bind,:);
            this.tfPositionStrands = this.tfPositionStrands(repmat(bind, [2 1]), :);
        end
    end

    %model
    methods
        %Calculate
        %- contribution to FBA objective
        %- minimum expression consistent with cell cycle length
        function [bmProd, byProd, minEnzExp, maxEnzExp] = calcResourceRequirements_LifeCycle(this, ~, ~)
            %no substrate and byproducts required
            bmProd = zeros(size(this.substrateWholeCellModelIDs));
            byProd = zeros(size(this.substrateWholeCellModelIDs));

            %no enzymes required
            minEnzExp = zeros(size(this.enzymeWholeCellModelIDs));
            maxEnzExp = Inf(size(this.enzymeWholeCellModelIDs));
        end

        %initialization:
        %- Simulation initializeState method initializes 1 undamaged chromosome
        %- Various processes may bind proteins to chromosome
        %- Here we bind transcription factors to the initial chromosome
        function initializeState(this)
            this.evolveState();
        end

        %resource requirements
        function result = calcResourceRequirements_Current(this)
            result = zeros(size(this.substrates));
        end

        %simulation
        function evolveState(this)
            %Bind transcription factors to accessible sites
            tfBoundPromoters = this.bindTranscriptionFactors();

            %Compute new fold changes
            this.rnaPolymerase.transcriptionFactorBindingProbFoldChange = ...
                this.calcBindingProbabilityFoldChange(tfBoundPromoters);
        end

        %Bind transcription factors to accessible sites.
        function tfBoundPromoters = bindTranscriptionFactors(this)
            %Estimate which binding sites exist and are accessible to the
            %transcription unit. (Note chromosome class which ensure that sites
            %are truly accessible, and only allow transcription factors to bind
            %to these sites).
            tfBoundPromoters = this.tfBoundPromoters;
            accessible = ~tfBoundPromoters & reshape(this.chromosome.isRegionPolymerized(this.tfPositionStrands, 1, false), [], 2);

            for i = 1:size(this.enzymes, 1)
                %For each transcription factor:
                %  Stochastically pick transcription unit promoters elements to
                %  bind from among free promoter elements weighted by the
                %  transcription factors' affinities for the promoters. (Note:
                %  each transcription factor is assumed to bind a distinct
                %  promoter region near each transcription unit such that
                %  transcription factors only compete within their species
                %  for binding to promoters.)

                if this.enzymes(i) <= 0
                    continue;
                end

                sites = find((this.tfIndexs == i) & accessible);
                if isempty(sites)
                    continue;
                end

                tf = this.bindProteinToChromosome(...
                    this.tfPositionStrands(sites, :), i, this.enzymes(i), this.tfAffinities(sites), [], false, 1, false, []);

                tfBoundPromoters(sites(tf)) = true;
            end
        end

        function foldChange = calcBindingProbabilityFoldChange(this, tfBoundPromoters)
            %reset fold change
            foldChange = ones(size(this.chromosome.transcriptionUnitStartCoordinates, 1), 2);

            %fold change effect of bound TFs
            for i = 1:size(this.tfIndexs, 1)
                for j = 1:size(this.tfIndexs, 2)
                    if ~tfBoundPromoters(i, j)
                        continue;
                    end

                    foldChange(this.tuIndexs(i, j), j) = foldChange(this.tuIndexs(i, j), j) * this.tfActivities(i, j);
                end
            end

            otherFoldChanges = prod(this.otherActivities .^ (this.enzymes(:, ones(1, size(foldChange, 1)))' > 0), 2);
            foldChange = foldChange .* [otherFoldChanges otherFoldChanges];
        end
    end

    %get methods of dependent local state
    methods
        function result = get.tfBoundPromoters(this)
            result = reshape(this.isDnaBound(this.tfPositionStrands, this.tfIndexs(:)), [], 2);
        end

        function result = get.boundTFs(this)
            result = histc(this.tfIndexs(logical(this.tfBoundPromoters(:,1))), 1:numel(this.enzymes));
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
classdef MaF8 < PROBLEM
% <multi/many> <real>
% MP-DMP

%------------------------------- Reference --------------------------------
% R. Cheng, M. Li, Y. Tian, X. Zhang, S. Yang, Y. Jin, and X. Yao, A
% benchmark test suite for evolutionary many-objective optimization,
% Complex & Intelligent Systems, 2017, 3(1): 67-81.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties(Access = private)
        Points; % Vertexes
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Parameter setting
            if isempty(obj.M); obj.M = 10; end
            obj.M        = max(obj.M,3);
            obj.D        = 2;
            obj.lower    = [-10000,-10000];
            obj.upper    = [10000,10000];
            obj.encoding = ones(1,obj.D);
            % Generate vertexes
            obj.Points  = [];
            [thera,rho] = cart2pol(0,1);
            [obj.Points(:,1),obj.Points(:,2)] = pol2cart(thera-(1:obj.M)*2*pi/obj.M,rho);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = pdist2(PopDec,obj.Points);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            [X,Y] = ndgrid(linspace(-1,1,ceil(sqrt(N))));
            ND    = inpolygon(X(:),Y(:),obj.Points(:,1),obj.Points(:,2));
            R     = pdist2([X(ND),Y(ND)],obj.Points);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 3
                [X,Y]    = ndgrid(linspace(-1,1,40));
                R        = pdist2([X(:),Y(:)],obj.Points);
                ND       = inpolygon(X(:),Y(:),obj.Points(:,1),obj.Points(:,2));
                R(~ND,:) = nan;
                R = {reshape(R(:,1),size(X)),reshape(R(:,2),size(X)),reshape(R(:,3),size(X))};
            else
                R = [];
            end
        end
        %% Display a population in the decision space
        function DrawDec(obj,Population)
            Draw(obj.Points([1:end,1],:),'-k','LineWidth',1.5,{'\it x\rm_1','\it x\rm_2',[]});
            Draw(obj.Points,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 1 1],'Markeredgecolor',[.4 .4 .4]);
            Draw(Population.decs);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [faces,vertices] = struct2mesh(mask3M,stlfile,qOffset,voxel_size)
%
% nii = load_untouch_nii(niifile);
% logicalmask = zeros(size(nii.img));
% logicalmask(find(nii.img)) = 1;

if nargin < 4 || ~exist('voxel_size','var') || isempty(voxel_size)
    voxel_size = [1 1 1];
end

if nargin < 3 || ~exist('qOffset','var') || isempty(qOffset)
    sz3M = size(mask3M);
    qOffset = -[voxel_size(1)*sz3M(1) voxel_size(2)*sz3M(2) voxel_size(3)*sz3M(3)]/2;
end

logicalmask = flip(flip(mask3M,1),2);

gridX =voxel_size(1)* [1:size(logicalmask,1)]  + qOffset(1);
gridY = voxel_size(2)*[1:size(logicalmask,2)] + qOffset(2);
gridZ = voxel_size(3)*[1:size(logicalmask,3)] + qOffset(3);

[faces,vertices] = CONVERT_voxels_to_stl(stlfile,logicalmask,gridX,gridY,gridZ,'binary');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [TorF, vstr, rdate] = have_feature_syngrid()
%HAVE_FEATURE_SYNGRID  Detect availability/version info for SynGrid
%
%   Feature detection function implementing 'syngrid' tag for HAVE_FEATURE
%   to detect availability/version of SynGrid, Synthetic Grid Creation for
%   MATPOWER (https://github.com/MATPOWER/mx-syngrid).
%
%   See also HAVE_FEATURE, SYNGRID.

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

TorF = exist('syngrid', 'file') == 2;
if TorF
    v = sgver('all');
    vstr = v.Version;
    rdate = v.Date;
else
    vstr = '';
    rdate = '';
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [scaledSel, sel, prefCond] = mv_selectivity(mv, conditions, shiftBaseline);
% [scaledSel sel, prefCond] = mv_selectivity(mv, <conditions=all>, <shiftBaseline=1> );
%
%
% Compute the selectivity index for a multi-voxel UI, given the selected
% conditions and threshold.
%
% Selectivity is defined as:
%   sel = (max - nonmax) / (max + abs(nonmax))
% where [max] is the amplitude of response to the "preferred" condition --
% by definition, the condition which produced the maximal response;
% [nonmax] is the set of response amplitudes to all other selected
% conditions. Response amplitudes are computed according the to
% event-related paramter 'ampType': see er_setParams, er_defaultParams.
%
% In addition to coding the degree of response, this also encodes the
% preferred condition in the following manner: if the first selected
% condition is preferred, the map ranges from 0-1; if the second, from 1-2;
% and so on. In general, the main value of the map is
%   (preferred condition-1) + sel.
%
% 'shiftBaseline': use the correction suggested by Alex Martin to enforce
% non-negative amplitudes, by shifting the amplitudes for such voxels to
% ensure they're positive. Default is 1.
%
% ras 05/06: broken off of mv_exportSelectivity into a separate function.

if notDefined('conditions'),
    conditions = mv.trials.condNums(mv.trials.condNums>0);
end
if notDefined('shiftBaseline'), shiftBaseline = 1; end

mv.params.selConds = conditions;
amps = mv_amps(mv);
nVoxels = size(amps,1);
nConds = size(amps,2);

%%%%% perform Alex Martin's correction if requested
if shiftBaseline
	mn = min(amps,[],2);
	hasNegResponse = find(mn < 0);
	offset = zeros(size(mn));
	offset(hasNegResponse) = -mn(hasNegResponse);
	offset = repmat(offset, [1 nConds]);
	amps = amps + offset;
end

%%%%% core part: compute selectivity
mx = max(amps,[],2); % max values
for i = 1:nVoxels
    if isnan(mx(i))
        prefCond(i) = 0;
        other(i,:) = 0;
        disp('NaN')
        continue;
    end

    % preferred condition
    prefCond(i) = find(amps(i,:)==mx(i));

    % amplitudes of other conditions
    other(i,:) = amps(i,setdiff(1:nConds,prefCond(i)));
end

nonmx = mean(other,2);
sel = (mx-nonmx) ./ (mx+abs(nonmx));

scaledSel = sel + prefCond' - 1;
for i=1:nConds
    ii=find(scaledSel>i-1);
    jj=find(scaledSel(ii) <i);
    fprintf(1,'cond %i numvoxels %i\n', i, length(jj));

end

return

author: vistalab