text:
function [errorCode, ima] = pco_sdk_example_single_EDITED(exposure_time,triggermode)
% set variables grab and display a single images
%
%   [errorCode] = pco_sdk_example_single(exposure_time,triggermode)
%
% * Input parameters :
%    exposure_time           camera exposure time (default=10ms)
%    triggermode             camera trigger mode (default=AUTO)
%
% * Output parameters :
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%
%grab images from a recording pco camera
%using script function pco_camera_stack
%display the grabbed images
%
%function workflow
%open camera
%set variables
%start camera
%grab image
%show image
%stop camera
%close camera
%

glvar=struct('do_libunload',0,'do_close',0,'camera_open',0,'out_ptr',[]);

if(~exist('exposure_time','var'))
 exposure_time = 10;
end

if(~exist('triggermode','var'))
 triggermode = 0;
end

%reduce_display_size=1: display only top-left corner 800x600Pixel
reduce_display_size=0;

pco_camera_load_defines();
subfunc=pco_camera_subfunction();
subfunc.fh_lasterr(0);

[errorCode,glvar]=pco_camera_open_close(glvar);
pco_errdisp('pco_camera_setup',errorCode);
disp(['camera_open should be 1 is ',int2str(glvar.camera_open)]);
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

out_ptr=glvar.out_ptr;

try

subfunc.fh_stop_camera(out_ptr);

cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,~,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
pco_errdisp('PCO_GetCameraDescription',errorCode);

if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
 subfunc.fh_enable_timestamp(out_ptr,TIMESTAMP_MODE_OFF); %TIMESTAMP_MODE_BINARYANDASCII
end

%enable MetaData if available
%if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_METADATA))
% subfunc.fh_set_metadata_mode(out_ptr,1);
%end

subfunc.fh_set_exposure_times(out_ptr,exposure_time,2,0,2)
subfunc.fh_set_triggermode(out_ptr,triggermode);

%if PCO_ArmCamera does fail no images can be grabbed

errorCode = calllib('PCO_CAM_SDK', 'PCO_ArmCamera', out_ptr);
if(errorCode~=PCO_NOERROR)
 pco_errdisp('PCO_ArmCamera',errorCode);
 ME = MException('PCO_ERROR:ArmCamera','Cannot continue script with not armed camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

%adjust transfer parameter if necessary
subfunc.fh_set_transferparameter(out_ptr);
subfunc.fh_get_triggermode(out_ptr);
subfunc.fh_show_frametime(out_ptr);

disp('get single image');
%subfunc.fh_start_camera(out_ptr);
[errorCode,ima,metadata,glvar]=pco_camera_stack(1,glvar);
if(errorCode==PCO_NOERROR)

% m=max(max(imar(10:end-10,10:end-10)));
% imshow(ima',[0,m+100]);
figure
imagesc(ima)
%draw_image(imar,[0 m+100]);
% disp(['found max ',int2str(m)]);
 %disp('Press "Enter" to proceed')
 %pause();
 %close()
end

subfunc.fh_stop_camera(out_ptr);

if(~isempty(metadata))
 [metadatasize,~]=size(metadata);
 [errorCode_m,metastruct]=subfunc.fh_get_struct_metadata(metadata,metadatasize);
 if(errorCode_m)
  pco_errdisp('get_struct_metadata',errorCode_m);
 else
  txt=subfunc.fh_print_meta_timestamp(metastruct);
  disp(['Timestamp of meta (',num2str(1,'%04d'),'):  ',txt]);
  subfunc.fh_print_meta_struct(metastruct);
 end
end

catch ME
 errorCode=subfunc.fh_lasterr();
 txt=blanks(101);
 txt=calllib('PCO_CAM_SDK','PCO_GetErrorTextSDK',pco_uint32err(errorCode),txt,100);

 clearvars -except ME glvar;

 if(glvar.camera_open==1)
  glvar.do_close=1;
  glvar.do_libunload=1;
  pco_camera_open_close(glvar);
 end

 if strfind(ME.identifier,'PCO_ERROR:')
  msg=[ME.identifier,' ',ME.message];
  disp(txt);
  warning('off','backtrace')
  warning(msg)
  for k=1:length(ME.stack)
    disp(['from file ',ME.stack(k).file,' at line ',num2str(ME.stack(k).line)]);
  end

  clearvars;
  commandwindow;
  return;
 else
  clearvars -except ME;
  rethrow(ME)
 end
end

clearvars -except glvar ima errorCode;

if(glvar.camera_open==1)
 glvar.do_close=1;
 glvar.do_libunload=1;
 pco_camera_open_close(glvar);
end

clearvars glvar;
commandwindow;
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function varargout = loadBiGGModel( model_ids, format, multichoice)
% Reads from url http://bigg.ucsd.edu/api/v2/models
%
% USAGE:
%
%    model = listBiGGModel(model_ids)
%
% OPTIONAL INPUTS:
%
%    model_ids:     The BiGG ID(s) of the model(s). If no ID is provided (either empty or no arguments) a gui
%                   will ask for the model to load.
%    format:        The format from Bigg to load. Either 'sbml' or ('mat').
%                   Default('mat')
%    multichoice:   Whether multiple models can be loaded if no ids are
%                   given. (Default: true)
% OUTPUT:
%    varargout:     The models in the order of the selected ids. If
%                   multiple models are selected in the dialog, a struct
%                   is provided with the models being fields of that
%                   struct.
%
% EXAMPLES:
%
%    1) Load a specific model from the mat file of the BiGG database
%           model = loadBiGGModel('e_coli_core');
%
%    2) Load several models from the BiGG database at once
%           [model1, model2, ... = loadBiGGModel({'bigg_id1','bigg_id2',...});
%
%    3) Load model from the BiGG database sbml files
%           model = loadBiGGModel('iJR904','sbml');
%
%
% .. Authors:
%       - Thomas Pfau Sep 2017

BIGG_FolderName = tempname;
mkdir(BIGG_FolderName);
finishup = onCleanup(@() rmdir(BIGG_FolderName,'s'));

if ~exist('format','var')
    format = 'mat';
end

if ~exist('multichoice','var')
    multichoice = true;
end
if multichoice
    selectionMode = 'multichoice';
else
    selectionMode = 'single';
end

if strcmp(format,'sbml')
    format = 'xml';
end

%Do not create a struct for multiple ids
structout = false;

modelfiles = {};

modelOptions = webread('http://bigg.ucsd.edu/api/v2/models');
modelOptions = modelOptions.results;
modelIDs = {modelOptions.bigg_id}';

if nargin == 0 || isempty(model_ids)
    orgNames = cellfun(@(x,y) strcat(x,' [',y,']'), {modelOptions.organism}, {modelOptions.bigg_id}, 'UniformOutput',0)' ;
    [sortedNames, order] = sort(orgNames);
    maxNameSize = max(cellfun(@numel,orgNames));

    if usejava('desktop')
        [s,v] = listdlg('Name','Model Selection','PromptString','Select Model(s) to download','SelectionMode',selectionMode, ...
                    'OKString','Load', 'ListString', sortedNames, 'ListSize', [maxNameSize*7,160] );
    else
        for i = 1:numel(sortedNames)
            fprintf('%i: %s\n',i,sortedNames{i});
        end
        s = input(['Please select a model (e.g. type 3 for ' sortedNames{3} '):'],'s');
        s = str2num(s);
    end

    model_ids = modelIDs(order(s));
    if numel(model_ids) > 1
        structout = true;
    end
else
    if ~iscell(model_ids)
        model_ids = {model_ids};
    end

    present = ismember(model_ids,modelIDs);
    if any(~present)
        nonpresent = model_ids(~present);
        model_ids = model_ids(present);
        fprintf('Could not find the following model ids in BiGG:\n');
        for i = 1:numel(nonpresent)
            fprintf('%s\n',nonpresent{i});
        end
    end
end

%If nothing is selected return nothing.
if numel(model_ids) == 0
       warning('No model selected.')
       varargout = {};
       return
end

varargout = cell(numel(model_ids),1);
for i = 1:numel(model_ids)
    filename = [BIGG_FolderName filesep model_ids{i} '.' format];
    url = ['http://bigg.ucsd.edu/static/models/' model_ids{i} '.' format];
    websave(filename,url);
    varargout{i} = readCbModel(filename);
end

if structout
    outstruct = struct();
    for i = 1:numel(model_ids)
        fieldID = regexprep(model_ids{i},'[^a-zA-Z0-9_]','');
        if isempty(regexp(fieldID,'^[a-zA-Z]'))
            fieldID = ['M_' fieldID];
        end
        outstruct.(fieldID) = varargout{i};
    end
    varargout = {outstruct};
end

author: opencobra

repo: cobratoolbox


text:
function [param]=sine_fit(x,y,fixed_params,initial_params,plot_flag)
% Optimization of parameters of the sine function to time series
%
% Syntax:
%       [param]=sine_fit(x,y)
%
%       that is the same that
%       [param]=sine_fit(x,y,[],[],[])     % no fixed_params, automatic initial_params
%
%       [param]=sine_fit(x,y,fixed_params)        % automatic initial_params
%       [param]=sine_fit(x,y,[],initial_params)   % use it when the estimation is poor
%       [param]=sine_fit(x,y,fixed_params,initial_params,plot_flag)
%
% param = [offset, amplitude, phaseshift, frequency]
%
% if fixed_params=[NaN, NaN , NaN , NaN]        % or fixed_params=[]
% optimization of offset, amplitude, phase shift and frequency (default)
%
% if fixed_params=[NaN, 1 , NaN , 1/(2*pi)]
% optimization of offset and phase shift of a sine of amplitude=1 and frequency=1/(2*pi)
%
% if fixed_params=[0, NaN , NaN , NaN]
% optimization of amplitude, phase shift and frequency of a sine of offset=0
%
%
% Example:
% %% generate data vectors (x and y)
% fsine = @(param,timeval) param(1) + param(2) * sin( param(3) + 2*pi*param(4)*timeval );
% param=[0 1 0 1/(2*pi)];  % offset, amplitude, phaseshift, frequency
% timevec=0:0.1:10*pi;
% x=timevec;
% y=fsine(param,timevec) + 0.1*randn(size(x));
%
% %% standard parameter estimation
% [estimated_params]=sine_fit(x,y)
%
% %% parameter estimation with forced 1.5 fixed amplitude
% [estimated_params]=sine_fit(x,y,[NaN 1.5 NaN NaN])
%
% %% parameter estimation without plotting
% [estimated_params]=sine_fit(x,y,[],[],0)
%
%
% Doubts, bugs: rpavao@gmail.com

% warning off

if nargin<=1 %fail
    fprintf('');
    help sine_fit
    return
end

automatic_initial_params=NaN(1,4);
automatic_initial_params(1)=median(y);
automatic_initial_params(2)=mean(minmax(y));
automatic_initial_params(3)=0;
[~,pos]=findpeaks(smooth(y,10));
automatic_initial_params(4)=1/max(diff(x(pos)));

% phas_shift estimator, not working...
% [~,pos]=max(y);
% automatic_initial_params(3)=mod(x(pos),automatic_initial_params(4));

if nargin==2 %simplest valid input
    fixed_params=NaN(1,4);
    initial_params=automatic_initial_params;
    plot_flag=1;
end
if nargin==3
    initial_params=automatic_initial_params;
    plot_flag=1;
end
if nargin==4
    plot_flag=1;
end

if exist('fixed_params','var')
    if isempty(fixed_params)
        fixed_params=NaN(1,4);
    end
end
if exist('initial_params','var')
    if isempty(initial_params)
        initial_params=automatic_initial_params;
    end
end
if exist('plot_flag','var')
    if isempty(plot_flag)
        plot_flag=1;
    end
end

f_str='f = @(param,timeval)';
free_param_count=0;
bool_vec=NaN(1,4);
for i=1:4;
    if isnan(fixed_params(i))
        free_param_count=free_param_count+1;
        f_str=[f_str ' param(' num2str(free_param_count) ')'];
        bool_vec(i)=1;
    else
        f_str=[f_str ' ' num2str(fixed_params(i))];
        bool_vec(i)=0;
    end
    if i==1; f_str=[f_str ' +']; end
    if i==2; f_str=[f_str ' * sin(']; end
    if i==3; f_str=[f_str ' + 2*pi*']; end
    if i==4; f_str=[f_str '*timeval );']; end
end

eval(f_str)

[BETA,RESID,J,COVB,MSE] = nlinfit(x,y,f,initial_params(bool_vec==1));

free_param_count=0;
for i=1:4;
    if isnan(fixed_params(i))
        free_param_count=free_param_count+1;
        param(i)=BETA(free_param_count);
    else
        param(i)=fixed_params(i);
    end
end

if plot_flag==1
    x_vector=min(x):mean(diff(x))/2:max(x);
    plot(x,y,'k.',x_vector,f(BETA,x_vector),'r-')
    xlim(minmax(x_vector))
end

author: Sable

repo: mcbench-benchmarks


text:
function [varargout] = isfile(varargin)

%ISFILE Determine if the input points to a file
%   TF = ISFILE(INPUT) returns true if INPUT points to a file and false otherwise.
%
% This is a compatibility function that should only be added to the path on
% MATLAB versions prior to R2017b. Note that currently this function only allows
% for a single string in the input, as opposed to the isfile function from MATLAB,
% which allows for multiple inputs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% see https://github.com/fieldtrip/fieldtrip/issues/899

alternatives = which(mfilename, '-all');
if ~iscell(alternatives)
  % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
  alternatives = {alternatives};
end

keep = true(size(alternatives));
for i=1:numel(alternatives)
  keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
  keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
end
alternatives = alternatives(keep);

if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
  % remove this directory from the path
  p = fileparts(mfilename('fullpath'));
  sw = warning('off','backtrace');
  warning('Removing "%s" from your path.\nSee http://www.fieldtriptoolbox.org/faq/should_i_add_fieldtrip_with_all_subdirectories_to_my_matlab_path/', p);
  warning(sw);
  rmpath(p);
  % call the original MATLAB function
  if exist(mfilename, 'builtin')
    [varargout{1:nargout}] = builtin(mfilename, varargin{:});
  else
    [varargout{1:nargout}] = feval(mfilename, varargin{:});
  end
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = isfile(input)

% deal with the input arguments
if nargin==1
  [input] = deal(varargin{1:1});
else
  error('incorrect number of input arguments')
end

tf = ismember(exist(input),[2 3 4 5 6]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};

author: spm


text:
function randOrthoMat=randOrthoMat(numDim)
%RANDORTHOMAT Generate a random orthogonal matrix of the desired
%             dimensionality (with respect to the Haar distribution over
%             the orthogonal group --kind of a "uniform" distribution).
%             This means that the dot product of any two columns is zero.
%             The orthogonal matrices generated have determinants of +1 or
%             -1.
%
%INPUTS: numDim The scalar number of dimensions that the random orthogonal
%               matrix should have.
%
%OUTPUTS: randOrthoMat A numDimXnumDim random orthogonal matrix with a
%                     determinant of either +1 or -1.
%
%As noted in [1], random orthogonal matrices can be obtained by first
%gennerating a matrix of independent standard normal random variables, then
%takign the QR decomposition of the matrix such that R has all positive
%elements along the diagonal. However, since Matlab's QR command can
%produce R matrices with negative diagonal elements, we have to effectively
%transfer the sign of R onto Q to make the elements of the diagonal of R
%all positive. This means multiplying Q by diag(sign(diag(R))).
%
%EXAMPLE:
%We will plot the random directions in 2D on the unit circle to show that
%this produces a relatively uniform set of rotation matrices.
% numPoints=500;
% figure(1)
% clf
% hold on
% for curPoint=1:numPoints
%     thePoint=randOrthoMat(2)*[1;0];
%     scatter(thePoint(1),thePoint(2),'.k','linewidth',2)
% end
%
%REFERENCE:
%[1] G. W. Stewart, "The efficient generation of random orthogonal matrices
%    with an application to condition estimators," SIAM Journal on
%    Numerical Analysis, vol. 17, no. 3, pp. 403-409, Jun. 1980.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

X=randn(numDim,numDim);
[Q,R]=qr(X);

randOrthoMat=Q*diag(sign(diag(R)));

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [infoPackets, version] = openCCF(filename)

% openCCF
%
% This script reads a .ccf file and outputs it in a structure.
%
%   filename:  Name of the file to be opened. If the fname is omitted
%              the user will be prompted to select a file.
%              DEFAULT: Will open Open File UI.
%
%   Channels   1-128 are neural channels
%   Channels 129-144 are external analog inputs
%   Channels 145-148 are external analog outpus
%   Channels
%
%   Kian Torab
%   kian@blackrockmicro.com
%   Blackrock Microsystems
%   Version 2.1.0.0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0:
%   - Initial release.
%
% 1.2.3.0:
%   - Minor bug fix that led to a crash in certain cases.
%
% 1.2.4.0:
%   - Minor bug fix regarding passing a filename variable to the function.
%
% 2.0.0.0:
%   - Implemented XML CCF file format.
%
% 2.1.0.0:
%   - Fixed a bug in loading nTrode groups with a base of 0.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

infoPackets.openCCFVersion = '2.1.0.0';

%% Parameters

% Which versions of .ccf's are supported
supportedVersions = {'3.6', '3.7', '3.8', '3.9'};
% Number of cbPKT_CHANINFO packets to read in the file
nValidPKT_CHANINFO = 160;

if( nargin ~= 1 )
    [filename, pathname] = getFile('*.ccf');
    if ~filename
        disp('No file was selected.');
        return;
    end
    fullfilename = [pathname filename];
elseif (nargin == 1)
    fullfilename = filename;
else
    if strcmpi(filename, 'ver')
    disp(['Version is ' infoPackets.openCCFVersion]);
    infoPackets = [];
    return;
    end
end

%% Open file
[fid, message] = fopen( fullfilename, 'rb', 'l');
if( fid == -1 ),
   warning(['Unable to open file: ' fullfilename '. Error message ' message]);
   return;
end;

%% Read header incl. version
head = fread(fid, 16, '*char*1' )';

if strcmp(head(1:5), 'cbCCF')
    version = deblank(strtrim(head(6:end))); % Strip possible leading space and trailing nulls
elseif strcmp(head(1:5), '<?xml')
    version = '3.9';
else
  warning('Not a well-formed .ccf file');
  fclose(fid);
  return;
end

% Save version info in the structure.
infoPackets.Version = version;

% Verify to see if the file version is supported
if ~ismember(version, supportedVersions)
  warning(sprintf('Unsupported file version: %s', version));
  fclose(fid);
  return;
end

if strcmpi(version, '3.9')
    infoPackets = parseCCF(fullfilename);
    for nTrodeIDX = 1:length(infoPackets.Children(7).Children)
        for trodeIDX = 1:length(infoPackets.Children(7).Children(nTrodeIDX).Children(9).Children)
            tempTrodes(trodeIDX) =  ...
            str2double(infoPackets.Children(7).Children(nTrodeIDX).Children(9).Children(trodeIDX).Children.Data) + 1;
        end
        if ~all(tempTrodes == 1)
            infoPackets.NTrodeInfo.NTrodeID(nTrodeIDX) = nTrodeIDX;
            infoPackets.NTrodeInfo.NTrodeMembers{nTrodeIDX} = tempTrodes;
        end
    end
    return;
end

%% Read packets
for p = 1:nValidPKT_CHANINFO
    if p == 128
%    disp('w');
end

  infoPackets.ChanInfo(p) = readInfoPacket(fid, version);
end
infoPackets.AdaptiveFilterInfo = readAdapFilterPacket(fid);
infoPackets.SortingInfo = readSortingPacket(fid);
infoPackets.SystemInfo = []; %readSystemInfoPacket(fid);
infoPackets.NTrodeInfo = readNTrodeInfoPacket(infoPackets);

%% Close file
fclose(fid);

function p = readInfoPacket(fid, version)
% Read a cbPKT_CHANINFO packet
    p = [];
    cbLEN_STR_LABEL = 16;
    cbMAXUNITS = 5;
    cbMAXHOOPS = 4;
    p.SystemClockTimestamp = fread(fid, 1, '*uint32');
    p.Child = fread(fid, 1, '*uint16');
    p.AnaInPacketType = fread(fid, 1, '*uint8');
    p.SystemDataLength = fread(fid, 1, '*uint8');
    p.AnaInChannelID = fread(fid, 1, '*uint32');
    p.AnaInChannelProcessor = fread(fid, 1, '*uint32');
    p.AnaInBankID = fread(fid, 1, '*uint32');
    p.AnaInPinID = fread(fid, 1, '*uint32');
    h = fread(fid, 1, '*uint32');
    p.CapabilitiesGeneralChannel = setFlagsCBCHAN(h);
    h = fread(fid, 1, '*uint32');
    p.CapabilitiesDigOutChannel = setFlagsCBDOUT(h);
    h = fread(fid, 1, '*uint32');
    p.CapabilitiesDigInChannel = setFlagsCBDINP(h);
    h = fread(fid, 1, '*uint32');
    p.CapabilitiesAnaOutChannel = setFlagsCBAOUT(h);
    h = fread(fid, 1, '*uint32');
    p.CapabilitiesAnaInChannel = setFlagsCBAINP(h);
    p.CapabilitiesSpikeProcessing = fread(fid, 1, '*uint32');

    p.InputPhysicalChannelScalingInfo       = readCbScaling(fid);
    p.InputPhysicalChannelFilterDefinition  = readFiltDesc(fid);
    p.OutputPhysicalChannelScalingInfo      = readCbScaling(fid);
    p.OutputPhysicalChannelFilterDefinition = readFiltDesc(fid);

    p.ChannelLabel     = fread(fid, cbLEN_STR_LABEL, '*char*1')';
    p.ChannelUserFlags = fread(fid, 1, '*uint32');
    p.ChannelPosition  = fread(fid, 4, '*int32')';

    p.InputPhysicalChannelScalingInfoCustom  = readCbScaling(fid);
    p.OuputPhysicalChannelScalingInfoCustom = readCbScaling(fid);

    p.doutopts = fread(fid, 1, '*uint32');
    p.dinpopts = fread(fid, 1, '*uint32');
    p.aoutopts = fread(fid, 1, '*uint32');
    p.eopchar  = fread(fid, 1, '*uint32');

    % Below two lines are actually part of a union, not sure how to figure out
    % which option is actually being used
    p.AnaInMonitorChannelAddress             = fread(fid, 1, '*uint32');
    p.AnaInOutputValue                       = fread(fid, 1, '*int32');
    p.AnaInOptions                           = fread(fid, 1, '*uint32');
    p.AnaInLNCAdaptationRate                 = fread(fid, 1, '*uint32');
    p.AnaInContinuousStreamFilterID          = fread(fid, 1, '*uint32');
    p.AnalogInputContinuousStreamSampleGroup = fread(fid, 1, '*uint32');
    p.AnaInContinuousStreamDisplayFactorMin  = fread(fid, 1, '*int32');
    p.AnaInContinuousStreamDisplayFactorMax  = fread(fid, 1, '*int32');
    p.AnaInSpikeStreamFilterID               = fread(fid, 1, '*uint32');
    p.AnaInSpikeStreamDisplayFactorMax       = fread(fid, 1, '*int32');
    p.AnaInLNCDisplayFactorMax               = fread(fid, 1, '*int32');
    p.AnaInSpikeProcessingOptions            = fread(fid, 1, '*uint32');
    p.AnaInSpikeStreamThresholdLevel         = fread(fid, 1, '*int32');
    p.AnaInSpikeStreamThresholdLimit         = fread(fid, 1, '*int32');
    p.AnaInNTrodeGroupID                     = fread(fid, 1, '*uint32');
    p.AnaInAmplitudeRejectionValuePositive   = fread(fid, 1, '*int16');
    p.AnaInAmplitudeRejectionValueNegative   = fread(fid, 1, '*int16');
    p.AnaInDigitalReferencingChannel         = fread(fid, 1, '*uint32');
    p.AnaInManualMappingUnit                 = readCbManualUnitMappings(fid, cbMAXUNITS, version);
    p.AnaInSpikeHoopsSorting                 = readCbHoops(fid, cbMAXUNITS, cbMAXHOOPS);

function p = readAdapFilterPacket(fid)
    p = [];
    p.Time = fread(fid, 1, '*uint32');
    p.Child = fread(fid, 1, '*uint16');
    p.Type = fread(fid, 1, '*uint8');
    p.DataLength = fread(fid, 1, '*uint8');
    p.Chan = fread(fid, 1, '*uint32');
    p.Mode = fread(fid, 1, '*uint32');
    p.LearningRate = fread(fid, 1, '*float32');
    p.RefChan1 = fread(fid, 1, '*uint32');
    p.RefChan2 = fread(fid, 1, '*uint32');

function p = readSortingPacket(fid)
    p = [];
    p.Time = fread(fid, 1, '*uint32');
    p.Child = fread(fid, 1, '*uint16');
    p.Type = fread(fid, 1, '*uint8');
    p.DataLength = fread(fid, 1, '*uint8');
    p.NumMaxSimultChans = fread(fid, 1, '*uint32');
    p.RefractoryPeriodInSamples = fread(fid, 1, '*uint32');

function p = readSystemInfoPacket(fid)
    p = [];

function p = readNTrodeInfoPacket(infoPackets)
    p = [];
    NTrodeGroups = unique([infoPackets.ChanInfo.AnaInNTrodeGroupID]); NTrodeGroups(1) = [];
    NTrodeGroupsCount = length(NTrodeGroups);

for ntrodeIDX = 1:NTrodeGroupsCount
    p.NTrodeID(ntrodeIDX) = NTrodeGroups(ntrodeIDX);
    p.NTrodeMembers{ntrodeIDX} = find([infoPackets.ChanInfo.AnaInNTrodeGroupID] == ntrodeIDX);
end

%% Sub-struct readers
function sc = readCbScaling(fid)
% Read a cbSCALING struct

cbLEN_STR_UNIT = 8;

sc.AnalogChannelScaleMin = fread(fid, 1, '*int16');
sc.AnalogChannelScaleMax = fread(fid, 1, '*int16');
sc.AnalogChannelScaleMin = fread(fid, 1, '*int32');
sc.AnalogChannelScaleMax = fread(fid, 1, '*int32');
sc.AnalogChannelGain = fread(fid, 1, '*int32');
sc.AnalogChannelUnit = fread(fid, cbLEN_STR_UNIT, '*char*1')';

function fd = readFiltDesc(fid)
% Read a cbFILTDESC struct

cbLEN_STR_FILT_LABEL = 16;

fd.FilterLabel = fread(fid, cbLEN_STR_FILT_LABEL, '*char*1')';
fd.FilterHighPassFrequency = fread(fid, 1, '*uint32');
fd.FilterHighPassOrder = fread(fid, 1, '*uint32');
fd.FilterHighPassType = fread(fid, 1, '*uint32');
fd.FilterLowPassFrequency = fread(fid, 1, '*uint32');
fd.FilterLowPassOrder = fread(fid, 1, '*uint32');
fd.FilterLowPassType = fread(fid, 1, '*uint32');

function ums = readCbManualUnitMappings(fid, maxUnits, version)
% Loop to read all the cbMANUALUNITMAPPING packets. Note that this requires
% knowing the file version: v3.6 is different from 3.7 and 3.8

% Loop growing is inefficient, but it's just not that much data.
if strcmp(version, '3.6')
    for i = 1:maxUnits
        ums(i) = readCbManualUnitMapping3_6(fid);
    end
else
    for i = 1:maxUnits
        ums(i) = readCbManualUnitMapping(fid);
    end
end

function um = readCbManualUnitMapping(fid)
% Read a cbMANUALUNITMAPPING packet for a v3.7 or 3.8 file

um.nOverride = fread(fid, 1,     '*int16');
um.afOrigin  = fread(fid, 3,     '*int16')';
um.afShape   = fread(fid, [3 3], '*int16')';
um.aPhi      = fread(fid, 1,     '*int16')';
um.bValid    = fread(fid, 1,     '*uint32')';

function um = readCbManualUnitMapping3_6(fid)
% Read a cbMANUALUNITMAPPING packet for a v3.6 file

um.nOverride = fread(fid, 1,     '*uint32');
um.afOrigin  = fread(fid, 3,     '*float32')';
um.afShape   = fread(fid, [3 3], '*float32')';
um.aPhi      = fread(fid, 1,     '*float32')';
um.bValid    = fread(fid, 1,     '*uint32')';

function hs = readCbHoops(fid, maxUnits, maxHoops)
% Loop to read the cbHOOP packets

% Loop growing is inefficient, but it's just not that much data.
for u = 1:maxUnits
  for h = 1:maxHoops
    hs(u, h) = readCbHoop(fid);
  end
end

function h = readCbHoop(fid)
% Read a cbHOOP packet
h.valid = fread(fid, 1, '*uint16')';
h.time = fread(fid, 1, '*int16')';
h.min = fread(fid, 1, '*int16')';
h.max = fread(fid, 1, '*int16')';

function h = setFlagsCBCHAN(flags)

binaryValue = dec2bin(flags,12);

% define constants
h.ChanExists     = 0;  % Channel id is allocated
h.ChanConnected  = 0;  % Channel is connected and mapped and ready to use
h.ChanIsolated   = 0;  % Channel is electrically isolated
h.ChanType       = '';

if strcmpi(binaryValue(1), '1')
    h.ChanType = 'DigitalOut';
end
if strcmpi(binaryValue(2), '1')
    h.ChanType = 'DigitalInput';
end
if strcmpi(binaryValue(3), '1')
    h.ChanType = 'AnalogOutput';
end
if strcmpi(binaryValue(4), '1')
    h.ChanType = 'AnalogInput';
end
if strcmpi(binaryValue(10), '1')
    h.ChanIsolated = 1;
end
if strcmpi(binaryValue(11), '1')
    h.ChanConnected = 1;
end
if strcmpi(binaryValue(12), '1')
    h.ChanExists = 1;
end

function h = setFlagsCBDOUT(flags)

binaryValue = dec2bin(flags,30);

h.BaudRate = 0;
h.NumberOfBits = 0;
h.CanBeManuallyConfigured = 0;
h.TrackMostRecentlySelected = 0;
h.CanOutputFrequency = 0;
h.ChanMonitorUnit0 = 0;
h.ChanMonitorUnit1 = 0;
h.ChanMonitorUnit2 = 0;
h.ChanMonitorUnit3 = 0;
h.ChanMonitorUnit4 = 0;
h.ChanMonitorUnit5 = 0;
h.ChanMonitorUnitAll = 0;

if strcmpi(binaryValue(1:6), '111111')
    h.ChanMonitorUnitAll = 1;
else
    if strcmpi(binaryValue(1), '1')
        h.ChanMonitorUnit5 = 1;
    end
    if strcmpi(binaryValue(2), '1')
        h.ChanMonitorUnit4 = 1;
    end
    if strcmpi(binaryValue(3), '1')
        h.ChanMonitorUnit3 = 1;
    end
    if strcmpi(binaryValue(4), '1')
        h.ChanMonitorUnit2 = 1;
    end
    if strcmpi(binaryValue(5), '1')
        h.ChanMonitorUnit1 = 1;
    end
    if strcmpi(binaryValue(6), '1')
        h.ChanMonitorUnit0 = 1;
    end
end
if strcmpi(binaryValue(12), '1')
    h.CanOutputFrequency = 1;
end
if strcmpi(binaryValue(13), '1')
    h.TrackMostRecentlySelected = 1;
end
if strcmpi(binaryValue(18), '1')
    h.CanBeManuallyConfigured = 1;
end
if strcmpi(binaryValue(19), '1')
    h.NumberOfBits = 32;
end
if strcmpi(binaryValue(20), '1')
    h.NumberOfBits = 16;
end
if strcmpi(binaryValue(21), '1')
    h.NumberOfBits = 8;
end
if strcmpi(binaryValue(22), '1')
    h.NumberOfBits = 1;
end
if strcmpi(binaryValue(25), '1')
    h.BaudRate = 115200;
end
if strcmpi(binaryValue(26), '1')
    h.BaudRate = 57600;
end
if strcmpi(binaryValue(27), '1')
    h.BaudRate = 38400;
end
if strcmpi(binaryValue(28), '1')
    h.BaudRate = 19200;
end
if strcmpi(binaryValue(29), '1')
    h.BaudRate = 9600;
end
if strcmpi(binaryValue(30), '1')
    h.BaudRate = 2400;
end

function h = setFlagsCBDINP(flags)

binaryValue = dec2bin(flags,22);

h.BaudRate = 0;
h.NumberOfBits = 0;
h.PortReadType = 0;
h.CapturePacketMethod = 0;
h.PortControlsEvent = 0;

if strcmpi(binaryValue(1), '1')
    h.PortReadType = '8-bit strobe/8-bit falling edge';
end
if strcmpi(binaryValue(2), '1')
    h.PortReadType = '8-bit strobe/8-bit rising edge';
end
if strcmpi(binaryValue(3), '1')
    h.PortReadType = '8-bit strobe/8-bit any edge';
end
if strcmpi(binaryValue(4), '1')
    h.PortReadType = '8-bit any bit falling edge';
end
if strcmpi(binaryValue(5), '1')
    h.PortReadType = '8-bit any bit rising edge';
end
if strcmpi(binaryValue(6), '1')
    h.PortControlsEvent = 1;
end
if strcmpi(binaryValue(7), '1')
    h.CapturePacketMethod = 'Logic input';
end
if strcmpi(binaryValue(8), '1')
    h.CapturePacketMethod = 'Character';
end
if strcmpi(binaryValue(9), '1')
    h.CapturePacketMethod = 'Word strob';
end
if strcmpi(binaryValue(10), '1')
    h.CapturePacketMethod = 'Any bit change';
end
if strcmpi(binaryValue(11), '1')
    h.NumberOfBits = 32;
end
if strcmpi(binaryValue(12), '1')
    h.NumberOfBits = 16;
end
if strcmpi(binaryValue(13), '1')
    h.NumberOfBits = 8;
end
if strcmpi(binaryValue(14), '1')
    h.NumberOfBits = 1;
end
if strcmpi(binaryValue(17), '1')
    h.BaudRate = 115200;
end
if strcmpi(binaryValue(18), '1')
    h.BaudRate = 57600;
end
if strcmpi(binaryValue(19), '1')
    h.BaudRate = 38400;
end
if strcmpi(binaryValue(20), '1')
    h.BaudRate = 19200;
end
if strcmpi(binaryValue(21), '1')
    h.BaudRate = 9600;
end
if strcmpi(binaryValue(22), '1')
    h.BaudRate = 2400;
end

function h = setFlagsCBAOUT(flags)

binaryValue = dec2bin(flags,22);

h.Audio = 0;
h.AnalogScale = 0;
h.ChanFunction = '';

if strcmpi(binaryValue(1), '1')
    h.ChanFunction = 'Waveform generator';
end
if strcmpi(binaryValue(2), '1')
    h.ChanFunction = 'Monitors spikes';
end
if strcmpi(binaryValue(3), '1')
    h.ChanFunction = 'Monitors SMP';
end
if strcmpi(binaryValue(4), '1')
    h.ChanFunction = 'Monitors LNC';
end
if strcmpi(binaryValue(5), '1')
    h.ChanFunction = 'Monitors raw';
end
if strcmpi(binaryValue(6), '1')
    h.AnalogScale = 'Waveform generator';
end
if strcmpi(binaryValue(7), '1')
    h.AnalogScale = 'Waveform generator';
end
if strcmpi(binaryValue(8), '1')
    h.AnalogScale = 'Waveform generator';
end
if strcmpi(binaryValue(9), '1')
    h.Audio = 1;
end

function h = setFlagsCBAINP(flags)

binaryValue = dec2bin(flags,22);

h.Audio = 0;
h.AnalogScale = 0;
h.ChanFunction = '';

if strcmpi(binaryValue(1), '1')
    h.ChanFunction = 'Waveform generator';
end
if strcmpi(binaryValue(2), '1')
    h.ChanFunction = 'Monitors spikes';
end
if strcmpi(binaryValue(3), '1')
    h.ChanFunction = 'Monitors SMP';
end
if strcmpi(binaryValue(4), '1')
    h.ChanFunction = 'Monitors LNC';
end
if strcmpi(binaryValue(5), '1')
    h.ChanFunction = 'Monitors raw';
end
if strcmpi(binaryValue(6), '1')
    h.AnalogScale = 'Waveform generator';
end
if strcmpi(binaryValue(7), '1')
    h.AnalogScale = 'Waveform generator';
end
if strcmpi(binaryValue(8), '1')
    h.AnalogScale = 'Waveform generator';
end
if strcmpi(binaryValue(9), '1')
    h.Audio = 1;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [c, acc] = compute_confusion(numClasses, gts, preds, areas, doNotNormalizePerClass)
if ~exist('doNotNormalizePerClass')
  doNotNormalizePerClass = false ;
end
if nargin <= 3, areas = ones(size(gts)) ; end
c = accumarray([gts(:), preds(:)], areas(:), numClasses*[1,1]) ;
if ~doNotNormalizePerClass
  c = bsxfun(@times, 1./sum(c,2), c) ;
  acc = mean(diag(c)) ;
else
  c = c / sum(c(:)) ;
  acc = sum(diag(c)) ;
end

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
% Convert an adjacency matrix of a general graph to the adjacency matrix of
%          a simple graph (symmetric, no loops, no double edges, no weights)
%
% INPUTS: adjacency matrix, nxn
% OUTPUTs: adjacency matrix (nxn) of the corresponding simple graph
%
% Other routines used: symmetrize.m
% GB: last updated, Sep 6 2014

function adj=adj2simple(adj)

adj=adj>0; % make all edges weight 1
adj = symmetrize(adj);
adj = adj - diag(diag(adj)); % clear the diagonal (selfloops)

author: aeolianine

repo: octave-networks-toolbox


text:
% Converts a frequency to some scale and back again
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef ScalingFunction < handle
    methods (Abstract)
        scale_to_hertz(obj, scale)
        hertz_to_scale(obj, hertz)
    end

end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%ReLU function

x = -10:0.01:10;
y = x;
y(x<0) = 0;
plot(x,y);
xlabel('x');
ylabel('y');
grid on

author: TheAlgorithms

repo: MATLAB-Octave


function [ywc1, ywt1, ywc2, ywt2] = wavefilter(X, fig)

% wavelets (Haar) filtering: see Lubik, Matthes, Verona (2019)

%  ywc1 has the  cycles at frequencies 8-32 quarters (loose 16 datapoints)
%  ywc2 has the  cycles at frequencies 8-64 quarters (loose 32 datapoints)
%  ywt1 and  ywt2 have  the trends  and are defined  as  residuals

enddT=size(X,1);
ywc1=zeros(enddT,size(X,2));   ywt1=zeros(enddT,size(X,2));
ywc2=zeros(enddT,size(X,2));   ywt2=zeros(enddT,size(X,2));

for qq=1:size(X,2)
%    y=zeros(enddT,1);
    xx1=zeros(enddT,1); xx2=zeros(enddT,1); xx3=zeros(enddT,1);
    xx4=zeros(enddT,1); xx5=zeros(enddT,1);

    y=squeeze(X(:,qq));

    % 2-4  quarters cycles
    for  tt=2:length(y)
    %xx1(tt,1)=(1/2)*(y(tt)-(y(tt)-y(tt-1)));
    xx1(tt,1)=(1/2)*(y(tt)-y(tt-1));
    end
    % 4-8  quarters cycles
    for  tt=4:length(y)
    xx2(tt,1)=(1/4)*(y(tt)+y(tt-1)-(y(tt-2)+y(tt-3)) );
    end
    % 8-16  quarters cycles
    for  tt=8:length(y)
    xx3(tt,1)=(1/8)*(y(tt)+y(tt-1)+y(tt-2)+y(tt-3)- ...
                  (y(tt-4)+y(tt-5)+y(tt-6)+y(tt-7)) );
    end
    % 16-32 quarters cycles
    for  tt=16:length(y)
    xx4(tt,1)=(1/16)*(y(tt)+y(tt-1)+y(tt-2)+y(tt-3)+y(tt-4)+ ...
                    y(tt-5)+y(tt-6)+y(tt-7)-...
                   (y(tt-8)+y(tt-9)+y(tt-10)+y(tt-11)+y(tt-12)+ ...
                    y(tt-13)+y(tt-14)+y(tt-15)) );
    end
    %  32-64 quarters cycles
    for  tt=32:length(y)
     xx5(tt,1)=(1/32)*(y(tt)+y(tt-1)+y(tt-2)+y(tt-3)+y(tt-4)+y(tt-5)+ ...
        y(tt-6)+y(tt-7)+y(tt-8)+y(tt-9)+y(tt-10)+ y(tt-11)+y(tt-12)+ ...
         y(tt-13)+ y(tt-14)+y(tt-15) - ...
        (y(tt-16)+y(tt-17)+y(tt-18)+y(tt-19)+y(tt-20)+y(tt-21)+ ...
         y(tt-22)+y(tt-23)+y(tt-24)+y(tt-25)+y(tt-26)+y(tt-27)+ ...
         y(tt-28)+y(tt-29)+y(tt-30)+y(tt-31) ) );
    end
   % cyclical components
   ywc1(16:length(y),qq)=xx3(16:length(y))+xx4(16:length(y));
   ywc2(32:length(y),qq)=xx3(32:length(y))+...
                       xx4(32:length(y))+xx5(32:length(y));
   % trend components (actual-BC cycles-high frequnecy cycles)
   ywt1(16:length(y),qq)=y(16:length(y))-ywc1(16:length(y),qq) ...
                         -xx1(16:length(y))-xx2(16:length(y));
   ywt2(32:length(y),qq)=y(32:length(y))-ywc2(32:length(y),qq) ...
                         -xx1(32:length(y))-xx2(32:length(y));

% xxx=[xx1 xx2 xx3 xx4  xx5]
% figure(100)
%  plot(xxx)

  if fig==1
      % plot  cyclical components
   figure(1)
   plot(ywc1(1:length(y),qq),'r', 'linewidth',2); hold  on;
   plot(ywc2(1:length(y), qq),'b','linewidth',2); hold  off; axis  tight;
   legend('BC(8-32)','BC+LOW(8-64)')
   pause
   end
 end

end


text:
function write_mrtrix (image, filename)

% function: write_mrtrix (image, filename)
%
% write the data contained in the structure 'image' in the MRtrix
% format image 'filename' (i.e. files with the extension '.mif' or '.mih').
%
% 'image' is either a N-dimensional array (N <= 16), or a structure containing
% the following fields:
%    image.data:            a N-dimensional array (N <= 16)
%    image.vox:             N-vector of voxel sizes (in mm) (default: { 2 }) [optional]
%    image.comments:        a cell array of strings [optional]
%    image.datatype:        the datatype specifier (default: float32) [optional]
%    image.mrtrix_version:  a character array [optional]
%    image.transform:       a 4x4 matrix [optional]
%    image.dw_scheme:       a NDWx4 matrix of gradient directions [optional]

fid = fopen (filename, 'w');
assert(fid ~= -1, 'error opening %s', filename);
fprintf (fid, 'mrtrix image\ndim: ');

if isstruct(image)
  dim = size(image.data);
else
  dim = size(image);
end

while prod(size(dim)) < 3
  dim(end+1) = 1;
end

fprintf (fid, '%d', dim(1));
fprintf (fid, ',%d', dim(2:end));

fprintf (fid, '\nvox: ');
if isstruct (image) && isfield (image, 'vox')
  fprintf (fid, '%.3f', image.vox(1));
  fprintf (fid, ',%.3f', image.vox(2:end));
else
  fprintf(fid, '2');
  fprintf(fid, ',%d', 2*ones(1,size(dim,2)-1));
end

fprintf (fid, '\nlayout: +0');
fprintf (fid, ',+%d', 1:(size(dim,2)-1));

[computerType, maxSize, endian] = computer;
if isstruct (image) && isfield (image, 'datatype')
  datatype = lower(image.datatype);
  byteorder = datatype(end-1:end);

  if strcmp (byteorder, 'le')
    precision = datatype(1:end-2);
    byteorder = 'l';
  elseif strcmp(byteorder, 'be')
    precision = datatype(1:end-2);
    byteorder = 'b';
  else
    if strcmp(datatype, 'bit')
      precision = 'bit1';
      byteorder = 'n';
    elseif strcmp (datatype, 'int8') || strcmp (datatype, 'uint8')
      precision = datatype;
      byteorder = 'n';
      if endian == 'L'
        datatype(end+1:end+3) = 'le';
      else
        datatype(end+1:end+3) = 'be';
      end
    end
  end
else
  if endian == 'L'
    datatype = 'float32le';
  else
    datatype = 'float32be';
  end
  precision = 'float32';
  byteorder = 'n';
end
fprintf (fid, [ '\ndatatype: ' datatype ]);

fprintf (fid, '\nmrtrix_version: %s', 'matlab');

if isstruct (image) && isfield (image, 'transform')
  fprintf (fid, '\ntransform: %.6f', image.transform(1,1));
  fprintf (fid, ',%.6f', image.transform(1,2:4));
  fprintf (fid, '\ntransform: %6f', image.transform(2,1));
  fprintf (fid, ',%.6f', image.transform(2,2:4));
  fprintf (fid, '\ntransform: %.6f', image.transform(3,1));
  fprintf (fid, ',%.6f', image.transform(3,2:4));
end

if isstruct (image) && isfield (image, 'dw_scheme')
  for i=1:size(image.dw_scheme,1)
    fprintf (fid, '\ndw_scheme: %.6f', image.dw_scheme(i,1));
    fprintf (fid, ',%.6f', image.dw_scheme(i,2:4));
   end
end

% write out any other fields:
if isstruct (image)
  f = fieldnames(image);
  % don't worry about fields that have already been dealt with:
  f(strcmp(f,'dim')) = [];
  f(strcmp(f,'vox')) = [];
  f(strcmp(f,'layout')) = [];
  f(strcmp(f,'datatype')) = [];
  f(strcmp(f,'transform')) = [];
  f(strcmp(f,'dw_scheme')) = [];
  f(strcmp(f,'data')) = [];

  % write out contents of the remainging fields:
  for n = 1:numel(f)
    val = getfield (image, f{n});
    if iscell (val)
      for i=1:numel(val)
        fprintf (fid, '\n%s: %s', f{n}, val{i});
      end
    else
      fprintf (fid, '\n%s: %s', f{n}, val);
    end
  end
end

if strcmp(filename(end-3:end), '.mif')
  datafile = filename;
  dataoffset = ftell (fid) + 24;
  fprintf (fid, '\nfile: . %d\nEND\n                         ', dataoffset);
elseif strcmp(filename(end-3:end), '.mih')
  datafile = [ filename(end-3:end) '.dat' ];
  dataoffset = 0;
  fprintf (fid, '\nfile: %s %d\nEND\n', datafile, dataoffset);
else
  fclose(fid);
  error('unknown file suffix - aborting');
end

fclose(fid);

fid = fopen (datafile, 'r+', byteorder);
fseek (fid, dataoffset, -1);

if isstruct(image)
  fwrite (fid, image.data, precision);
else
  fwrite (fid, image, precision);
end
fclose (fid);


text:
function tests = test_kmeans_pp
  tests = functiontests(localfunctions);
end

function test_1(testCase)
    X = [1 2 7 8 20 21
         1 2 7 8 20 21];
    % Capture current random number generator state
    st = rng;
    % Go to default state
    rng('default');
    % Perform testing
    [seeds, labels] = spx.cluster.kmeans.pp_initialize(X, 3);
    verifyEqual(testCase, seeds, [
        20     8     1
        20     8     1]);
    verifyEqual(testCase, labels, [3     3     2     2     1     1]);
    % restore previous state
    rng(st);
end

author: indigits

repo: sparse-plex


text:
% rbm_get_visible
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [x] = rbm_get_visible(h0, R)

n_visible = size(R.W, 1);
n_hidden = size(R.W, 2);

if R.data.binary == 1
    x = sigmoid(bsxfun(@plus,binornd(1, h0) * R.W', R.vbias'));
else
    x = bsxfun(@plus, h0 * R.W', R.vbias');
end
end

author: kyunghyuncho

repo: deepmat


text:
% Calculate arena coverage (the amount of space an animal has covered during experiment)
%
% Calculates arena coverage.
%
%  USAGE
%   coverage = analyses.arenaCoverage(pos, binWidth, shape, dimensions)
%   pos         Position samples, matrix of size at least Nx2. Format is either [t x] or [t x y].
%   binWidth    Width of horizontal and vertical bins in cm. If only one value is provided, then
%               the same bin width is used in both directions.
%   shape       Arena shape, integer. One of the values of bntConstants.ArenaShape.
%   dimensions  Vector of arena dimensions, the actual number of elements depends on arena shape.
%
%   coverage    Arena coverage, float in range [0..100].
%
function coverage =  arenaCoverage(pos, binWidth, shape, dimensions)
    inp = inputParser;

    checkPosDimensions = @(x) size(x, 2) >= 2;
    checkBinWidth = @(x) helpers.isdvector(x, '>=0') && length(x) <= 2;
    checkShape = @(x) ismember(x, helpers.ArenaShape.allShapes());
    checkDims = @(x) helpers.isdvector(x, '>=0') && (x(1) > 0);

    addRequired(inp, 'pos', checkPosDimensions);
    addRequired(inp, 'binWidth', checkBinWidth);
    addRequired(inp, 'shape', checkShape);
    addRequired(inp, 'dimensions', checkDims);

    parse(inp, pos, binWidth, shape, dimensions);

    t = pos(:, 1);
    x = pos(:, 2);
    if size(pos, 2) > 2
        y = pos(:, 3);
    else
        y = [];
    end

    binWidthX = binWidth(1);
    if length(binWidth) == 1
        binWidthY = binWidthX;
    else
        binWidthY = binWidth(2);
    end

    % filter out points that lie outside arena. This is important for circles
    % because they do not fully cover all the bins.
    if shape == bntConstants.ArenaShape.Circle
        radius = dimensions(1)/2;
        ind = sqrt(x.^2 + y.^2) > radius;
        x(ind) = nan;
        y(ind) = nan;
    end

    limitsX = [-dimensions(1)/2 dimensions(1)/2];
    [xBinned, nBinsX, edgesX] = helpers.bin(x, limitsX, binWidthX);
    nBins = nBinsX;

    if ~isempty(y)
        if length(dimensions) == 1
            if shape == bntConstants.ArenaShape.Track
                yLength = 1;
            else
                yLength = dimensions(1);
            end
        else
            if shape == bntConstants.ArenaShape.Track && dimensions(2) == 0
                yLength = 1;
            else
                yLength = dimensions(2);
            end
        end

        limitsY = [-yLength/2 yLength/2];
        [yBinned, nBinsY, edgesY] = helpers.bin(y, limitsY, binWidthY);
        nBins = [nBinsX nBinsY];
    end

    dt = diff(t);
    dt(end+1) = dt(end);

    if isempty(y)
        occupancy = general.accumulate(xBinned, dt, nBinsX)';
    else
        occupancy = general.accumulate([xBinned yBinned], dt, nBins)';
    end

    switch shape
        case bntConstants.ArenaShape.Circle
            radius = dimensions(1)/2;
            radiusBin = ceil(radius / binWidthX) + 1;

            halfSize = ceil(size(occupancy)/2);
            [rr, cc] = meshgrid(1:nBinsX, 1:nBinsY);

            distMap = sqrt((cc - halfSize(2)).^2 + (rr - halfSize(1)).^2); % each element is the distance
                                                % from the middle of the map to current point
            outerCircle = distMap >= radiusBin;
            distMap(outerCircle) = nan;

            numBins = sum(sum(isfinite(distMap)));
            coverage = (length(find(occupancy > 0)) / numBins) * 100;
        otherwise
            coverage = length(find(occupancy > 0)) / prod(nBins) * 100;
    end
end

author: kavli-ntnu


text:
function opt_vect = foptions()
% FOPTIONS Sets default parameters for optimisation routines
% For compatibility with MATLAB's foptions()
%
% Copyright (c) Dharmesh Maniyar, Ian T. Nabney (2004)

opt_vect      = zeros(1, 18);
opt_vect(2:3) = 1e-4;
opt_vect(4)   = 1e-6;
opt_vect(16)  = 1e-8;
opt_vect(17)  = 0.1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% dslam_open is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with dslam_open. If not, see <http://www.gnu.org/licenses/>.

function opt_handle = initOptHandle(num_robots)

for i = 1:num_robots
    opt_handle.optimized_groups(i).members = [i];
end

opt_handle.running = false;

end

author: uzh-rpg


text:
function out = scatplot(x,y,method,radius,N,n,po,ms,toplot)
% Scatter plot with color indicating data density
%
% USAGE:
%   out = scatplot(x,y,method,radius,N,n,po,ms)
%   out = scatplot(x,y,dd)
%
% DESCRIPTION:
%   Draws a scatter plot with a colorscale
%   representing the data density computed
%   using three methods
%
% INPUT VARIABLES:
%   x,y - are the data points
%   method - is the method used to calculate data densities:
%       'circles' - uses circles with a determined area
%               centered at each data point
%       'squares' - uses squares with a determined area
%               centered at each data point
%       'voronoi' - uses voronoi cells to determin data densities
%               default method is 'voronoi'
%   radius - is the radius used for the circles or squares
%       used to calculate the data densities if
%       (Note: only used in methods 'circles' and 'squares'
%           default radius is sqrt((range(x)/30)^2 + (range(y)/30)^2)
%   N - is the size of the square mesh (N x N) used to
%       filter and calculate contours
%       default is 100
%   n - is the number of coeficients used in the 2-D
%       running mean filter
%       default is 5
%       (Note: if n is length(2), n(2) is tjhe number of
%       of times the filter is applied)
%   po - plot options:
%       0 - No plot
%       1 - plots only colored data points (filtered)
%       2 - plots colored data points and contours (filtered)
%       3 - plots only colored data points (unfiltered)
%       4 - plots colored data points and contours (unfiltered)
%           default is 1
%   ms - uses this marker size for filled circles
%       default is 4
%
% OUTPUT VARIABLE:
%   out - structure array that contains the following fields:
%       dd - unfiltered data densities at (x,y)
%       ddf - filtered data densities at (x,y)
%       radius - area used in 'circles' and 'squares'
%               methods to calculate densities
%       xi - x coordenates for zi matrix
%       yi - y coordenates for zi matrix
%       zi - unfiltered data densities at (xi,yi)
%       zif - filtered data densities at (xi,yi)
%       [c,h] = contour matrix C as described in
%           CONTOURC and a handle H to a contourgroup object
%       hs = scatter points handles
%
%Copy-Left, Alejandro Sanchez-Barba, 2005

if nargin==0
    scatplotdemo
    return
end
if nargin<3 | isempty(method)
    method = 'vo';
end
if nargin<9
    toplot = 1;
end
if isnumeric(method)
   gsp(x,y,method,2)
   return
else
    method = method(1:2);
end
if nargin<4 | isempty(n)
    n = 5; %number of filter coefficients
end
if nargin<5 | isempty(radius)
    radius = sqrt((range(x)/30)^2 + (range(y)/30)^2);
end
if nargin<6 | isempty(po)
    po = 1; %plot option
end
if nargin<7 | isempty(ms)
    ms = 4; %markersize
end
if nargin<8 | isempty(N)
    N = 100; %length of grid
end
%Correct data if necessary
x = x(:);
y = y(:);
%Asuming x and y match
idat = isfinite(x);
x = x(idat);
y = y(idat);
holdstate = ishold;
if holdstate==0
    cla
end
hold on
%--------- Caclulate data density ---------
dd = datadensity(x,y,method,radius);
%------------- Gridding -------------------
xi = repmat(linspace(min(x),max(x),N),N,1);
yi = repmat(linspace(min(y),max(y),N)',1,N);
zi = griddata(x,y,dd,xi,yi);
%----- Bidimensional running mean filter -----
zi(isnan(zi)) = 0;
coef = ones(n(1),1)/n(1);
zif = conv2(coef,coef,zi,'same');
if length(n)>1
    for k=1:n(2)
        zif = conv2(coef,coef,zif,'same');
    end
end
%-------- New Filtered data densities --------
ddf = griddata(xi,yi,zif,x,y);
%----------- Plotting --------------------
if toplot
    switch po
        case {1,2}
            if po==2
                [c,h] = contour(xi,yi,zif);
                out.c = c;
                out.h = h;
            end %if
            hs = gsp(x,y,ddf,ms);
            out.hs = hs;
            colorbar
        case {3,4}
            if po>3
                [c,h] = contour(xi,yi,zi);
                out.c = c;
            end %if
            hs = gsp(x,y,dd,ms);
            out.hs = hs;
            colorbar
    end %switch
end
%------Relocate variables and place NaN's ----------
dd(idat) = dd;
dd(~idat) = NaN;
ddf(idat) = ddf;
ddf(~idat) = NaN;
%--------- Collect variables ----------------
out.dd = dd;
out.ddf = ddf;
out.radius = radius;
out.xi = xi;
out.yi = yi;
out.zi = zi;
out.zif = zif;
if ~holdstate
    hold off
end
return
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function scatplotdemo
po = 2;
method = 'squares';
radius = [];
N = [];
n = [];
ms = 5;
x = randn(1000,1);
y = randn(1000,1);

out = scatplot(x,y,method,radius,N,n,po,ms)

return
%~~~~~~~~~~ Data Density ~~~~~~~~~~~~~~
function dd = datadensity(x,y,method,r)
%Computes the data density (points/area) of scattered points
%Striped Down version
%
% USAGE:
%   dd = datadensity(x,y,method,radius)
%
% INPUT:
%   (x,y) -  coordinates of points
%   method - either 'squares','circles', or 'voronoi'
%       default = 'voronoi'
%   radius - Equal to the circle radius or half the square width
Ld = length(x);
dd = zeros(Ld,1);
switch method %Calculate Data Density
    case 'sq'  %---- Using squares ----
        for k=1:Ld
            dd(k) = sum( x>(x(k)-r) & x<(x(k)+r) & y>(y(k)-r) & y<(y(k)+r) );
        end %for
        area = (2*r)^2;
        dd = dd/area;
    case 'ci'
        for k=1:Ld
            dd(k) = sum( sqrt((x-x(k)).^2 + (y-y(k)).^2) < r );
        end
        area = pi*r^2;
        dd = dd/area;
    case 'vo'  %----- Using voronoi cells ------
        [v,c] = voronoin([x,y]);
        for k=1:length(c)
            %If at least one of the indices is 1,
            %then it is an open region, its area
            %is infinity and the data density is 0
            if all(c{k}>1)
                a = polyarea(v(c{k},1),v(c{k},2));
                dd(k) = 1/a;
            end %if
        end %for
end %switch
return
%~~~~~~~~~~ Graf Scatter Plot ~~~~~~~~~~~
function varargout = gsp(x,y,c,ms)
%Graphs scattered poits
map = colormap;
ind = fix((c-min(c))/(max(c)-min(c))*(size(map,1)-1))+1;
h = [];
%much more efficient than matlab's scatter plot
for k=1:size(map,1)
    if any(ind==k)
        h(end+1) = line('Xdata',x(ind==k),'Ydata',y(ind==k), ...
            'LineStyle','none','Color',map(k,:), ...
            'Marker','.','MarkerSize',ms);
    end
end
if nargout==1
    varargout{1} = h;
end
return

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function [hpatch, cl] = cluster_image_sphere(cl, varargin)
% :Usage:
% ::
%
%    [hpatch, cl] = cluster_image_sphere(cl or [k x 3 list of mm coords], varargin)
%
% Images spheres at cluster centers
% Combine with addbrain and cluster_tools('connect3d')
% or cluster_nmdsfig_glassbrain to get 3-D plots of connected blobs
% or surface_cutaway
%
% :Optional Inputs:
%
%   {'color', 'colors'}, mycolor = varargin{i+1};
%
%   'radius', myradius = varargin{i+1};
%
% :Outputs:
%
%   patch handles, and new cl with sphere coordiates in XYZmm and XYZ
%
% :Example: Usage
% ::
%
%    function [hpatch, cl] = cluster_image_sphere(cl)
%
%    % With optional arguments:
%    [hpatch, cl] = cluster_image_sphere(cl, 'color', 'b', 'radius', 10)
%    [hpatch, cl] = cluster_image_sphere(cl, 'color', {'r' 'g' 'b' etc}, 'radius', 10)
%
% :Examples: Given an MNI coordinate, plot a sphere on a brain surface
% ::
%
%    my_mm_coord = [40, 46, 22]';
%    create_figure('surface')
%    cl = [];
%    cl.XYZmm = my_mm_coord;
%    cl.mm_center = my_mm_coord';
%    V = spm_vol(which('brainmask.nii'));
%    cl.M = V.mat;
%    [hpatch, cl] = cluster_image_sphere(cl, 'color', 'g', 'radius', 10)
%    p = addbrain;
%    set(p, 'FaceAlpha', 1);
%    axis image
%    view(135, 30); lighting gouraud; lightRestoreSingle; material dull;
%
% Example: Turn xyz mm coordinates into clusters and image them
% ::
%
%    my_mm_coord = [40 46 22; 50 26 40; 45 36 50; 60 12 0]
%    [hpatch, cl] = cluster_image_sphere(my_mm_coord, 'color', 'b', 'radius', 4);
%
% ..
%    Tor Wager, July 2007 (original version)
%
%    Programmers' notes:
%    Tor Wager, July 2007
%    updated April 2011 for flexible radius
%    updated 12/2012 for xyz to spheres
% ..

    mycolor = 'r';
    myradius = 8;

    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % reserved keywords
                case {'color', 'colors'}, mycolor = varargin{i+1}; varargin{i+1} = [];

                case 'radius', myradius = varargin{i+1};

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
    end

    if isempty(cl)
        hpatch = [];
        return
    end

    if isa(cl, 'region')
       cl = region2struct(cl);
    end

    % if cl is actually a series of xyz coordinates...
    if ~isstruct(cl)
        cl = create_cl_struct(cl);
    end

    if length(myradius) == 1
        myradius = repmat(myradius, 1, length(cl));
    end

    % get colors for each sphere
    if length(mycolor) == 1 || ~iscell(mycolor)
        mycolor = repmat({mycolor}, 1, length(cl));
    end

   %newcl = convert_cl_xyz_to_sphere(cl(1), myradius);

    for i = 1:length(cl)

        newcl(i) = convert_cl_xyz_to_sphere(cl(i), myradius(i));

        hpatch(i) = image_isosurface(newcl(i), mycolor{i});

    end

    cl = newcl;

end

function cl = create_cl_struct(xyz)

V = spm_vol(which('brainmask.nii'));
if size(xyz, 2) ~= 3, xyz = xyz'; end
if size(xyz, 2) ~= 3, error('Bad input xyz coordinates?'); end

cl = [];
for i = 1:size(xyz, 1)
    my_mm_coord = xyz(i, :);
    cl(i).XYZmm = my_mm_coord';
    cl(i).mm_center = my_mm_coord;
    cl(i).M = V.mat;

    cl(i).voxSize = diag(cl(i).M(1:3, 1:3))';
end

end

function cl = convert_cl_xyz_to_sphere(cl, myradius)

    newXYZmm = round(points_in_sphere(cl(1).mm_center, myradius)');

    cl(1).Z = ones(1, size(newXYZmm,2));
    cl(1).XYZmm = newXYZmm;
    cl(1).XYZ = mm2voxel(cl(1).XYZmm, cl(1).M, 1)';
    cl(1).numVox = size(cl(1).XYZ, 2);

end

function hpatch = image_isosurface(cl, mycolor)

    % controls padding to make sure we cover whole area
    padval = 1;

    % controls smoothness, etc.
    mythresh = .5;
    mysmoothbox = 3;
    mygaussstd = 1;

    mm_coords = cl(1).XYZmm;

    xyzmin = min(mm_coords') - padval;     % minus/plus for padding
    xyzmax = max(mm_coords') + padval;

    [X, Y, Z] = meshgrid(xyzmin(1):xyzmax(1), xyzmin(2):xyzmax(2), xyzmin(3):xyzmax(3));

    % construct volume data for area to image

    xvox = mm_coords(1,:) - xyzmin(1) + 1;
    yvox = mm_coords(2,:) - xyzmin(2) + 1;
    zvox = mm_coords(3,:) - xyzmin(3) + 1;

    V = zeros(size(X));

    for i = 1:size(xvox, 2)
        V(xvox(i), yvox(i), zvox(i)) = cl(1).Z(i);
    end

    % not needed if we have all mm points in sphere
    %VI = interp3(cl(1).XYZmm(1,:), cl(1).XYZmm(2,:), cl(1).XYZmm(3,:), cl(1).Z, X, Y, Z);

    V = smooth3(V, 'gaussian', mysmoothbox, mygaussstd);
    FV = isosurface(X,Y,Z,V, mythresh);

    hpatch = patch(FV);
    set(hpatch, 'EdgeColor', 'none', 'FaceColor', mycolor);

    drawnow

    %lighting gouraud

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function t = structn(obj)
% Convert a NIFTI-1 object into a form of struct
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

if numel(obj)~=1,
    error('Too many elements to convert');
end;
fn = fieldnames(obj);
for i=1:length(fn)
    tmp = subsref(obj,struct('type','.','subs',fn{i}));
    if ~isempty(tmp)
        t.(fn{i}) = tmp;
    end;
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       : nrtMshSegment.m                               |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 25.11.2018                                    |
%| ( === ) |   SYNOPSIS   : Edge mesh of a segment                        |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

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

% Create mesh
Nvtx = 21;
L    = 4;
mesh = mshSegment(Nvtx,L);

% Colours
ctr = 1/2 .* (...
        mesh.vtx(mesh.elt(:,1),:) + ...
        mesh.vtx(mesh.elt(:,2),:) ) ;
mesh.col(ctr(:,1)<0)               = 1;
mesh.col(ctr(:,1)>=0 & ctr(:,1)<1) = 2;
mesh.col(ctr(:,1)>=1)              = 3;

% Graphical representation
figure
plot(mesh)
colorbar

% Sub-meshing
mesh1 = mesh.sub(mesh.col==1);
mesh2 = mesh.sub(mesh.col==2);
mesh3 = mesh.sub(mesh.col==3);
norm(mesh1.col-1,'inf')
norm(mesh2.col-2,'inf')
norm(mesh3.col-3,'inf')

% Graphical representation
figure
plot(mesh1)
hold on
plot(mesh2)
plot(mesh3)
colorbar

% Center
Xctr = mesh.ctr;
hold on
plot3(Xctr(:,1),Xctr(:,2),Xctr(:,3),'*y')

% Normals
plotNrm(mesh,'r')

% Volume
V = mesh.ndv;
norm(sum(V)-L)/L

% Length
l = mesh.stp;

% Edges
figure
plot(mesh.edg)
colorbar

% Particles
figure
plot(mesh.prt)
colorbar

% Unicity
tmp     = mesh;
tmp.vtx = [tmp.vtx ; tmp.vtx];
tmp.elt = [tmp.elt ; tmp.elt + size(mesh.vtx,1)];
tmp.col = [tmp.col ; tmp.col];
[~,I]   = unique(tmp);
norm(I-(1:size(mesh.elt,1))','inf')

% Intersection
[~,I] = intersect(mesh1,mesh);
norm(I-(1:size(mesh1.elt,1))','inf')

% Union
[tmp,I] = union(mesh,mesh1);
norm(I-(1:size(mesh.elt,1))','inf')

% Difference
[mesh4,~] = setdiff(mesh,mesh1);
figure
plot(mesh4)
colorbar

% Reconstruction
mesh5 = union(mesh1,mesh2);
mesh5 = union(mesh5,mesh3);
setdiff(mesh,mesh5)

% Boundary
bound = mesh5.bnd;
figure
plot(bound)
colorbar

% Clean degenerated mesh
mesh6 = mesh;
ind   = (mesh6.vtx(:,1)>0);
mesh6.vtx(ind,1) = 0;
mesh6 = clean(mesh6,1e-6);
setdiff(mesh6,mesh1)

% Clean degenerated particle mesh
vtx = mesh.vtx;
vtx(vtx(:,1)>0,1) = 0;
mesh7 = msh(vtx);
mesh7 = clean(mesh7,1e-6);
setdiff(mesh7,msh(mesh1.vtx))

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [A, B] = zeroKronApprox( K, m, n )
%
%       [A, B] = zeroKronApprox( K, m, n );
%
% computes a kronecker sum approximation to the blurring matrix K
% that arises from the input PSF under zero boundary conditions.
% This approximation is done a-la Kamm-Nagy (see paper for details).
%
%  Input:
%         K - psfMatrix object
%         m - size of matrix A (assumed square)
%         n - size of matrix B (assumed square)
%
%  Output:
%      Matrices A and B such that K \aprrox A \otimes B.
%

%  J. Nagy  2/11/02

%  Modifications:  This used to be zeroSepMatrix, which did not compute
%                  A and B.
%  5/25/02, J. Nagy
%           This now computes A and B.
%
%  11/17/02, J. Nagy
%            This was designed for image processing problems, where
%            it is common to use lexicographical (row) ordering.
%            But @kronMatrix functions where designed using vec (column)
%            ordering.  This inconsistency has been fixed.

P1 = K.psf;
P2 = P1.image;
PSF = P2{1};
c1 = P1.center;
center = c1{1};

[mp, np] = size(PSF);

if ( mp ~= np )
  error('For now, we expect PSF to be square')
end

%
% Compute weighted PSF.
%
for i = 1:center(1)
  Aweights(i,1) = sqrt( i+mp-center(1) );
end;
for i = center(1)+1:mp
  Aweights(i,1) = sqrt( mp+center(1)-i);
end;
for i = 1:center(2)
  Bweights(i,1) = sqrt( i+np-center(2) );
end;
for i = center(2)+1:np
  Bweights(i,1) = sqrt( np+center(2)-i );
end;

Phat = (Aweights*Bweights').*PSF;
%
% Compute SVD of weighted PSF, which is then used to construct
% the separable approximation.
%
[U,S,V] = svd( Phat );

%
% check to make sure first column looks like
% a Gaussian, and is not inverted.
%
minU = abs(min(min(U(:,1))));
maxU = max(max(abs(U(:,1))));
if minU == maxU
  U = -U;
  V = -V;
end

%
% Construct approximation.
%
a = ( U(:,1) * sqrt(S(1,1)) ) ./ Aweights;
b = ( V(:,1) * sqrt(S(1,1)) ) ./ Bweights;

%
%  This construction corresponds to lexicographical ordering,
%  but kronMatrix does everything corresponding to vec ordering.
%  Thus, these two do not work together.
%%
%%A = build_toep(a, center(1), n);
%%B = build_toep(b, center(2), m);

%  To fix this, we just need to switch A and B.  Then kronMatrix
%  can continue to be used corresponding to a vec ordering.
%
A = build_toep(b, center(2), m);
B = build_toep(a, center(1), n);

repo: IRtools


text:
% this function convert the input to a correspounding kVA.

function OutVal = SP_2kva(InVal)

global KWTransformersTypes TransformersTypes
Tag = find(KWTransformersTypes == InVal);
OutVal = TransformersTypes(1,Tag);


text:
function l = gammaPriorLogProb(prior, x)

% GAMMAPRIORLOGPROB Log probability of Gamma prior.

% PRIOR

% Compute log prior
D = length(x);
l = D*prior.a*log(prior.b)-D*gammaln(prior.a)+(prior.a-1)*sum(log(x))-prior.b*sum(x);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Ml = hmxPlus(Ml,Mr)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx 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       : hmxPlus.m                                     |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Sum of H-Matrix with the rule                 |
%|  `---'  |                Full > H-Matrix > Compr                       |
%+========================================================================+

% Check dimensions
if sum(size(Ml) ~= size(Mr))
    error('hmxPlus.m : matrix dimensions must agree.')
end

%%% H-Matrix + H-Matrix -> H-Matrix
if isa(Ml,'hmx') && isa(Mr,'hmx')
    % H-Matrix + H-Matrix --> H-Matrix (recursion)
    if (Ml.typ == 0) && (Mr.typ == 0)
        % Construction
        for i = 1:4
            Ml.chd{i} = hmxPlus(Ml.chd{i},Mr.chd{i});
        end

        % Fusion
        Ml = hmxFusion(Ml);

    % H-Matrix + Compr -> H-Matrix
    elseif (Ml.typ == 0) && (Mr.typ == 1)
        Ml = hmxPlusAB(Ml,Mr.dat{1},Mr.dat{2});

    % H-Matrix + Full -> Unavailable
    elseif (Ml.typ == 0) && (Mr.typ == 2)
        error('hmxPlus : unvailable case')

    % Compr + --- -> ---
    elseif (Ml.typ == 1)
        Ml = hmxPlusAB(Mr,Ml.dat{1},Ml.dat{2});

    % Full + H-Matrix -> Unavailable
    elseif (Ml.typ == 2) && (Mr.typ == 0)
        error('hmxPlus : unvailable case')

    % Full + Compr -> Full
    elseif (Ml.typ == 2) && (Mr.typ == 1)
        Ml = hmxPlusAB(Ml,Mr.dat{1},Mr.dat{2});

    % Full + Full -> Full
    elseif (Ml.typ == 2) && (Mr.typ == 2)
        Ml.dat = Ml.dat + Mr.dat;

    else
        error('hmxPlus : unvailable case')
    end

%%% H-Matrix + Matrix -> H-Matrix
elseif isa(Ml,'hmx')
    Ml = Ml + hmx(Ml.pos{1},Ml.pos{2},Mr,Ml.tol);

%%% Matrix + H-Matrix -> Matrix
elseif isa(Mr,'hmx')
    Ml = hmx(Mr.pos{1},Mr.pos{2},Ml,Mr.tol) + Mr;

%%% Unavailable
else
    error('hmxPlus.m : unavailable case')
end
end


text:
clear all; close all; clear classes; clc;

%% Ag in Table I on p.4374 of Johnson and Christy, 1972 Phys. Rev. B
eV = [0.64 0.77 0.89 1.02 1.14 1.26 1.39 1.51 1.64 1.76 1.88 2.01 2.13 2.26 2.38 2.50 2.63 2.75 2.88 3.00 3.12 3.25 3.37 3.50 3.62 3.74 3.87 3.99 4.12 4.24 4.36 4.49 4.61 4.74 4.86 4.98 5.11 5.23 5.36 5.48 5.60 5.73 5.85 5.98 6.10 6.22 6.35 6.47 6.60];
n = [0.24 0.15 0.13 0.09 0.04 0.04 0.04 0.04 0.03 0.04 0.05 0.06 0.05 0.06 0.05 0.05 0.05 0.04 0.04 0.05 0.05 0.05 0.07 0.10 0.14 0.17 0.81 1.13 1.34 1.39 1.41 1.41 1.38 1.35 1.33 1.31 1.30 1.28 1.28 1.26 1.25 1.22 1.20 1.18 1.15 1.14 1.12 1.10 1.07];
k = [14.08 11.85 10.10 8.828 7.795 6.992 6.312 5.727 5.242 4.838 4.483 4.152 3.858 3.586 3.324 3.093 2.869 2.657 2.462 2.275 2.070 1.864 1.657 1.419 1.142 0.829 0.392 0.616 0.964 1.161 1.264 1.331 1.372 1.387 1.393 1.389 1.378 1.367 1.357 1.344 1.342 1.336 1.325 1.312 1.296 1.277 1.255 1.232 1.212];

%% Convert the photon energies to the wavelengths.
wvlen = PhysC.h * PhysC.c0 * 1e9 ./ eV;

%% Make the row vectors into column vectors to be used in interp1q().
eV = eV.';
n = n.';
k = k.';
wvlen = wvlen.';

%% Calculate the permittivity from n and k following the exp(+i w t) time dependence.
eps = (n - 1i*k).^2;

%% Plot real(eps) and imag(eps).  Compare with Fig.3 on p.4375 of Johnson.

nk_wvlen = 1;
eps_eV = 2;
eps_wvlen = 3;
eps_omega = 4;
q_omega = 5;
plotstyle = eps_eV;
switch plotstyle
    case nk_wvlen  % plot n and k
        loglog(wvlen, n, 'o-', wvlen, k, 'o-')
        %plot(wvlen, n, wvlen, k)
        legend('n', 'k', 'Location', 'SouthEast');
        xlabel 'wavelength (nm)'
        %axis([1e2 1e4 1e-2 1e2])
    case eps_eV  % plot real(eps) and -imag(eps)
        plot(eV, real(eps), 'o-', eV, -imag(eps), 'o-')
        legend('\epsilon_1', '\epsilon_2', 'Location', 'SouthEast');
        xlabel 'Photon Energy (eV)'
%         axis([0.5 6.5 -7 7]);
    case eps_wvlen
        plot(wvlen, real(eps), 'o-', wvlen, -imag(eps), 'o-')
        legend('\epsilon_1', '\epsilon_2', 'Location', 'SouthEast');
        xlabel 'wavelength (nm)'
        %axis([1e2 1e4 1e-2 1e2])
    case eps_omega
        plot(2*pi./wvlen, real(eps), 'o-', 2*pi./wvlen, -imag(eps), 'o-')
        legend('\epsilon_1', '\epsilon_2', 'Location', 'SouthEast');
        xlabel '\omega (c/nm)'
    case q_omega
        omega = 2*pi./wvlen;
        deps1 = real(eps(2:end)) - real(eps(1:end-1));
        domega = omega(2:end) - omega(1:end-1);
        eps1_inter = (real(eps(2:end)) + real(eps(1:end-1)))/2;
        omega_inter = (omega(2:end) + omega(1:end-1))/2;
        numer = eps1_inter + omega_inter .* (deps1./domega);
        denom = -(imag(eps(2:end)) + imag(eps(1:end-1)));  % extra factor 2
        plot(2*pi./omega_inter, numer./denom, 'o-');
        xlabel 'wavelength (nm)'
        ylabel 'electric Q';
end
title(mfilename)

%% Save data.
save(mfilename, 'eV', 'n', 'k');

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
% Two examples to show how to solve an MILP and an LP with interior method

disp('-- Integer problem --');
s=1;
c=[-1,-1]';
a=[-2,5;2,-2];
b=[5;1];
ctype=['U','U']';
lb=[0;0]; ub=[];
vartype=['B';'B'];
param.msglev=3;
[xmin,fmin,status,extra]=glpk(c,a,b,lb,ub,ctype,vartype,s,param)
% --- OBSOLETE ---
% [xmin,fmin,status,extra]=glpkmex(s,c,a,b,ctype,lb,ub,vartype,param)
pause;

disp('3rd problem');
s=1;
c=[0 0 0 -1 -1]';
a=[-2 0 0 1 0;...
    0 1 0 0 2;...
    0 0 1 3 2];
b=[4 12 18]';
ctype=['S','S','S']';
lb=[0,0,0,0,0]'; ub=[];
vartype=['C','C','C','C','C']';
param.lpsolver=2;
[xmin,fmin,status,extra]=glpk(c,a,b,lb,ub,ctype,vartype,s,param)
% --- OBSOLETE ---
% [xmin,fmin,status,extra]=glpkmex(s,c,a,b,ctype,lb,ub,vartype)

author: opencobra

repo: cobratoolbox


text:
impts=zeros(yw*exf,xw*exf);
n_rendered=0;
weight=str2double(get(handles.weight,'String'));
size_fac=str2double(get(handles.size_fac_edit,'String'));

if get(handles.adv_cb,'Value') %if using advanced thresholding don't need r0 tolerances
    if get(handles.lp_N_cb,'Value') %if specifying tolerances for loc. prec.
        for i=nstart:nend
            if xc(i)>=1 && yc(i)>=1 && xc(i)<xw*exf && yc(i)<yw*exf && lp(i)*1000 >= lp_tol_min && lp(i)*1000 <= lp_tol_max && a0_err_all(i)/a0_all(i) <= fr_un && xf_err_all(i) <= max_unc && yf_err_all(i) <= max_unc
              wide=ceil(size_fac*lppix(i)*1.5+1);
              if xc(i)-wide>=1 && xc(i)+wide<xw*exf && yc(i)-wide>=1 && yc(i)+wide<yw*exf
                n_rendered=n_rendered+1;
                for j=xc(i)-wide:xc(i)+wide
                  for k=yc(i)-wide:yc(i)+wide
                    dx=double(j)-xf(i);
                    dy=double(k)-yf(i);
                    int=pi*lp2pix(i)*size_fac;
                    a=exp(-2*(dx*dx+dy*dy)/(size_fac*size_fac*lp2pix(i)))*N(i)*weight/int;
                    impts(k,j)=impts(k,j)+a;
                  end
                end
              end
            end
            waitbarxmod(i/nend,w); %update
        end
    else %if specifying tolerances for # of photons/molecule (N)
        for i=nstart:nend
            if xc(i)>=1 && yc(i)>=1 && xc(i)<xw*exf && yc(i)<yw*exf && N(i) >= N_tol_min && N(i) <= N_tol_max && a0_err_all(i)/a0_all(i) <= fr_un && xf_err_all(i) <= max_unc && yf_err_all(i) <= max_unc
              wide=ceil(size_fac*lppix(i)*1.5+1);
              if xc(i)-wide>=1 && xc(i)+wide<xw*exf && yc(i)-wide>=1 && yc(i)+wide<yw*exf
                n_rendered=n_rendered+1;
                for j=xc(i)-wide:xc(i)+wide
                  for k=yc(i)-wide:yc(i)+wide
                    dx=double(j)-xf(i);
                    dy=double(k)-yf(i);
                    int=pi*lp2pix(i)*size_fac;
                    a=exp(-2*(dx*dx+dy*dy)/(size_fac*size_fac*lp2pix(i)))*N(i)*weight/int;
                    impts(k,j)=impts(k,j)+a;
                  end
                end
              end
            end
            waitbarxmod(i/nend,w); %update
        end
    end
else %if not using advanced thresholding include r0 in tolerances
    if get(handles.lp_N_cb,'Value') %if specifying tolerances for loc. prec.
        for i=nstart:nend
            if xc(i)>=1 && yc(i)>=1 && xc(i)<xw*exf && yc(i)<yw*exf && lp(i)*1000 >= lp_tol_min && lp(i)*1000 <= lp_tol_max && r0_all(i) >= r0_tol_min && r0_all(i) <= r0_tol_max && r0_err_all(i)/r0_all(i) <= fr_un && a0_err_all(i)/a0_all(i) <= fr_un && xf_err_all(i) <= max_unc && yf_err_all(i) <= max_unc
              wide=ceil(size_fac*lppix(i)*1.5+1);
              if xc(i)-wide>=1 && xc(i)+wide<xw*exf && yc(i)-wide>=1 && yc(i)+wide<yw*exf
                n_rendered=n_rendered+1;
                for j=xc(i)-wide:xc(i)+wide
                  for k=yc(i)-wide:yc(i)+wide
                    dx=double(j)-xf(i);
                    dy=double(k)-yf(i);
                    int=pi*lp2pix(i)*size_fac;
                    a=exp(-2*(dx*dx+dy*dy)/(size_fac*size_fac*lp2pix(i)))*N(i)*weight/int;
                    impts(k,j)=impts(k,j)+a;
                  end
                end
              end
            end
            waitbarxmod(i/nend,w); %update
        end
    else %if specifying tolerances for # of photons/molecule (N)
        for i=nstart:nend
            if xc(i)>=1 && yc(i)>=1 && xc(i)<xw*exf && yc(i)<yw*exf && N(i) >= N_tol_min && N(i) <= N_tol_max && r0_all(i) >= r0_tol_min && r0_all(i) <= r0_tol_max && r0_err_all(i)/r0_all(i) <= fr_un && a0_err_all(i)/a0_all(i) <= fr_un && xf_err_all(i) <= max_unc && yf_err_all(i) <= max_unc
              wide=ceil(size_fac*lppix(i)*1.5+1);
              if xc(i)-wide>=1 && xc(i)+wide<xw*exf && yc(i)-wide>=1 && yc(i)+wide<yw*exf
                n_rendered=n_rendered+1;
                for j=xc(i)-wide:xc(i)+wide
                  for k=yc(i)-wide:yc(i)+wide
                    dx=double(j)-xf(i);
                    dy=double(k)-yf(i);
                    int=pi*lp2pix(i)*size_fac;
                    a=exp(-2*(dx*dx+dy*dy)/(size_fac*size_fac*lp2pix(i)))*N(i)*weight/int;
                    impts(k,j)=impts(k,j)+a;
                  end
                end
              end
            end
            waitbarxmod(i/nend,w); %update
        end
    end
end

thresh=impts*0+1;
impts=thresh.*(impts>thresh)+impts.*(impts<=thresh);
clear thresh;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function rootdirexamples = findexamplesdir()

rootdir = fileparts(which('SnirfClass.m'));
rootdirexamples = filesepStandard([rootdir, '/Examples/']);

author: BUNPC


text:
function geometry_test0236 ( )

%*****************************************************************************80
%
%% TEST0236 tests DODEC_SIZE_3D, DODEC_SHAPE_3D, SHAPE_PRINT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0236\n' );
  fprintf ( 1, '  For the dodecahedron,\n' );
  fprintf ( 1, '  DODEC_SIZE_3D returns dimension information;\n' );
  fprintf ( 1, '  DODEC_SHAPE_3D returns face and order information.\n' );
  fprintf ( 1, '  SHAPE_PRINT_3D prints this information.\n' );

  [ point_num, edge_num, face_num, face_order_max ] = dodec_size_3d ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num );
  fprintf ( 1, '  Number of faces =    %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max );

  [ point_coord, face_order, face_point ] = dodec_shape_3d ( point_num, ...
    face_num, face_order_max );

  shape_print_3d ( point_num, face_num, face_order_max, ...
    point_coord, face_order, face_point );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% LFSelectFromDatabase - support function for selecting white image/calibration by matching serial/zoom/focus
%
% Usage:
%
%     SelectedCamInfo = LFSelectFromDatabase( DesiredCamInfo, DatabaseFname )
%
% This helper function is used when decoding a light field to select an appropriate white image,
% and when rectifying a light field to select the appropriate calibration. It works by parsing a
% database (white image or calibration), and searching its CamInfo structure for the best match to
% the requested DesiredCamInfo. DesiredCamInfo is set based on the camera settings used in
% measuring the light field. See LFUtilDecodeLytroFolder / LFLytroDecodeImage for example usage.
%
% Selection prioritizes the camera serial number, then zoom, then focus. It is unclear whether this
% is the optimal approach.
%
% The output SelectedCamInfo includes all fields in the database's CamInfo struct, including the
% filename of the selected calibration or white image. This facilitates decoding / rectification.
%
% User guide: LFToolbox.pdf
% See also: LFUtilProcessWhiteImages, LFUtilProcessCalibrations, LFUtilDecodeLytroFolder, LFLytroDecodeImage

% Copyright (c) 2013-2020 Donald G. Dansereau

function SelectedCamInfo = LFSelectFromDatabase( DesiredCamInfo, DatabaseFname )

%---Load the database---
load(DatabaseFname, 'CamInfo');

%---Find the closest to the desired settings, prioritizing serial, then zoom, then focus---
ValidSerial = find( ismember({CamInfo.CamSerial}, {DesiredCamInfo.CamSerial}) );

% Discard non-matching serials
CamInfo = CamInfo(ValidSerial);
OrigIdx = ValidSerial;

% Find closest zoom
ZoomDiff = abs([CamInfo.ZoomStep] - DesiredCamInfo.ZoomStep);
BestZoomDiff = min(ZoomDiff);
BestZoomIdx = find( ZoomDiff == BestZoomDiff ); % generally multiple hits

% Retain the (possibly multiple) matches for the best zoom setting
% CamInfo.ZoomStep = CamInfo(BestZoomIdx).ZoomStep;
CamInfo = CamInfo(BestZoomIdx);
OrigIdx = OrigIdx(BestZoomIdx);

% Of those that are closest in zoom, find the one that's closest in focus
FocusDiff = abs([CamInfo.FocusStep] - DesiredCamInfo.FocusStep);
[~,BestFocusIdx] = min(FocusDiff);

% Retrieve the index into the original
BestOriglIdx = OrigIdx(BestFocusIdx);
SelectedCamInfo = CamInfo(BestFocusIdx);

repo: LFToolbox


text:
function Tr = loadCalibrationRigid(filename)

% open file
fid = fopen(filename,'r');

if fid<0
  error(['ERROR: Could not load: ' filename]);
end

% read calibration
R  = readVariable(fid,'R',3,3);
T  = readVariable(fid,'T',3,1);
Tr = [R T;0 0 0 1];

% close file
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = readVariable(fid,name,M,N)

% rewind
fseek(fid,0,'bof');

% search for variable identifier
success = 1;
while success>0
  [str,success] = fscanf(fid,'%s',1);
  if strcmp(str,[name ':'])
    break;
  end
end

% return if variable identifier not found
if ~success
  A = [];
  return;
end

% fill matrix
A = zeros(M,N);
for m=1:M
  for n=1:N
    [val,success] = fscanf(fid,'%f',1);
    if success
      A(m,n) = val;
    else
      A = [];
      return;
    end
  end
end

repo: MSCKF


text:
classdef OptimalSuperEllipsePrinter < handle

    properties (Access = private)
        mesh
        meshBackground
        levelSet
        mx
        my
        iter
    end

    properties (Access = private)
        mxV
        myV
        fileName
    end

    methods (Access = public)

        function obj = OptimalSuperEllipsePrinter()
            obj.init();
            obj.compute();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.mxV = [0.8462 0.99 0.01 0.99 0.30108 0.1 0.2 0.2  0.85];
            obj.myV = [0.8462 0.2  0.01 0.99 0.72687 0.1 0.6 0.95 0.85];
            obj.fileName = 'OptimaSuperEllipseMicroEllipse';
          %  iter = 0;
        end

        function compute(obj)
            for iTest = 1:length(obj.mxV)
                obj.iter = iTest;

                rho = SuperEllipseParamsRelator.rho(obj.mxV(obj.iter),obj.myV(obj.iter),32)

                obj.createMeshBackground();
                obj.createLevelSet();
                obj.createMesh();
                obj.print();
            end
        end

        function createMeshBackground(obj)
            s.testName = 'RVE_Square_Triangle_FineFine';
            obj.meshBackground = Mesh().createFromFile(s);
        end

        function createLevelSet(obj)
            sM.coord  = obj.meshBackground.coord;
            sM.connec = obj.meshBackground.connec;
            s.mesh = Mesh_Total(sM);
            s.widthH = obj.mxV(obj.iter);
            s.widthV = obj.myV(obj.iter);
            %s.pnorm  = obj.computeSmoothingExponent();
            s.pnorm  = 2;
            s.type = 'smoothRectangle';
            s.levelSetCreatorSettings = s;
            s.type = 'LevelSet';
            s.scalarProductSettings.epsilon = 1;
            obj.levelSet = LevelSet(s);
        end

        function q = computeSmoothingExponent(obj)
            s.m1 = obj.mxV(obj.iter);
            s.m2 = obj.myV(obj.iter);
            s.type = 'Optimal';
            qComputer = SmoothingExponentComputer.create(s);
            q = qComputer.compute();
        end

        function createMesh(obj)
            s.fileName = obj.fileName;
            s.levelSet = obj.levelSet.value;
            s.meshBackground = obj.meshBackground;
            mCreator = MeshCreatorFromLevelSetWithMMG(s);
            obj.mesh = mCreator.create();
        end

        function print(obj)
            outputName = [obj.fileName,'Print',num2str(obj.iter)];
            s.mesh       = obj.mesh;
            s.outPutName = outputName;
            printer = SuperEllipsePrinter(s);
            printer.print();
            printer.captureImage();
        end

    end

end


text:
function [pcaDims, lambda] = GridSearch(trainX, trainY, testX, testY)
%% [pcaDims, lambda] = GridSearch(trainX, trainY, testX, testY)
%
% This function demostrates how to use the development set to learn the
% optimal parameters.
%
% Inputs:
%   trainX, trainY: training data and class labels.
%   testX, testY: test data and class labels.
% Outputs:
%   pcaDims: the optimal PCA dimensions.
%   lambda: the optimal ridge regression parameter.
%
% Version: 1.0
% Date: 2014-07-22
%
% Author: Shengcai Liao
% Institute: National Laboratory of Pattern Recognition,
%   Institute of Automation, Chinese Academy of Sciences
% Email: scliao@nlpr.ia.ac.cn

candiPcaDims = 100:100:500; % candidate parameter values for the pca dimensions.
candiLambda = 10.^(-4:0); % candidate parameter values for the ridge regression.
veriFarPoints = [0, kron(10.^(-8:-1), 1:9), 1]; % FAR points for verification performance evaluation

fprintf('Learn parameters on the development set...\n');

%% Learn a PCA subspace.
W = PCA(trainX);

%% Transform both the training and test data into the learned PCA subspace.
trainX = trainX * W;
testX = testX * W;

%% Select classes which have at least two images.
hst = hist(trainY, 1 : max(trainY));
classIndex = find(hst >= 2);
[sampleIndex, trainY] = ismember(trainY, classIndex); % class labels are continuously relabeled starting from 1
trainY = trainY(sampleIndex);
trainX = trainX(sampleIndex,:);
numTrainSamples = length(trainY);

%% Construct the discriminant response matrix
Y = zeros(numTrainSamples, length(classIndex));
Y( sub2ind(size(Y), 1 : numTrainSamples, trainY') ) = 1;

%% Pre-compute the working variables for the ridge regression.
R = trainX' * trainX;
Z = trainX' * Y;

numPara1 = length(candiPcaDims);
numPara2 = length(candiLambda);
auc = zeros(numPara1, numPara2);

%% Learn the optimal parameters by grid search
for i = 1 : numPara1
    d = candiPcaDims(i);

    for j = 1 : numPara2
        % Learn a linear subspace by ridge regression. You can replace this with
        % your own supervised learning algorithm here.
        W = (R(1:d, 1:d) + candiLambda(j) * numTrainSamples * eye(d)) \ Z(1:d,:);

        % Transform the test data into the learned subspace.
        X = testX(:, 1:d) * W;

        % Normlize each row to unit length. If you do not have this function,
        % do it manually.
        X = normr(X);

        % Compute the cosine similarity score between the test samples.
        score = X * X';

        % Evaluate the verification performance.
        VR = EvalROC(score, testY, [], veriFarPoints);

        % Average the verification rates as the AUC performance
        auc(i,j) = mean(VR);
    end
end

%% Plot the AUC performance w.r.t. different PCA dimensions.
figure; semilogx(candiPcaDims, auc, 'LineWidth', 2);
grid on;
xlabel('PCA Dimension');
ylabel('AUC');
title('AUC performance w.r.t. different PCA dimensions');
labels = cellstr( [repmat('lambda = ', [length(candiLambda),1]), num2str(candiLambda')] );
legend(labels, 'Location', 'NorthWest');
set(gca, 'XTick', candiPcaDims);
set(gca, 'XTickLabel', cellstr(num2str(candiPcaDims')));
drawnow;

%% Get the best parameters.
[bestAuc, index] = max(auc(:));
[r,c] = ind2sub(size(auc), index);
pcaDims = candiPcaDims(r);
lambda = candiLambda(c);
fprintf('The best AUC: %g. The best PCA dimensions: %d. The best lambda: %g.\n\n', bestAuc, pcaDims, lambda);

author: happynear

repo: FaceVerification


text:
function [fMc, fBvalue, fAvalue, mu, fSigma] = calc_McEMR(catalog, binInterval)
    % Determine Mc using Entire Magnitude Range (EMR)-method. Calculates also a- and b-value.
    % [fMc, fBvalue, fAvalue, mu, fSigma] = calc_McEMR(catalog, binInterval);
    % -----------------------------------------------------------------------------------------------------
    % Determine Mc using EMR-method. Calculates also a- and b-value.
    % Fitting non-cumulative frequency magnitude distribution above and below Mc:
    % below: Cumulative NORMAL distribution function
    % above: Gutenberg-Richter law
    %
    % Incoming variables:
    % catalog   : EQ catalog
    % binInterval   : Binning interval, usually 0.1
    %
    % Outgoing variables:
    % fMc        : Best estimated magnitude of completeness
    % fBvalue    : b-value
    % fAvalue    : a-value
    % mu        : mu-value of the normal CDF
    % fSigma     : sigma-values of the normal CDF
    %
    % J. Woessner: woessner@seismo.ifg.ethz.ch

    narginchk(2,2);

    % Initialize
    vProbability = [];
    vMc = [];
    vABValue =[];
    mFitRes = [];
    vDeltaBest = [];
    vX_res = [];
    vNmaxBest = [];
    mResult=[];

    % Determine exact time period
    timespan = years(max(catalog.Date) - min(catalog.Date)); % guessing it should be years

    % Set starting value for Mc loop and LSQ fitting procedure
    fMcTry= calc_Mc(catalog, McMethods.MaxCurvature);
    fSmu = abs(fMcTry / 2);
    fSSigma = abs(fMcTry / 4);
    if (fSmu > 1)
        fSmu = fMcTry / 10;
        fSSigma = fMcTry / 20;
    end
    fMcBound = fMcTry;

    % Calculate FMD for original catalog
    [vFMDorga, vNonCFMDorg, fmdbins] = calc_FMD(catalog.Magnitude);
    % convert answer back to this file's expectations...
    vFMDorg = [fmdbins'; vFMDorga']; % as rows
    vNonCFMDorg = [fmdbins'; vNonCFMDorg'];

    fMinMag = min(vNonCFMDorg(1,:));

    % %% Shift to positive values
    % if fMinMag ~= 0
    %     fMcBound = fMcTry-fMinMag;
    % end
    % Loop over Mc-values
    for fMc = round(fMcBound-0.4 : 0.1 : fMcBound+0.4, -1)
        vFMD = vFMDorg;
        vNonCFMD = vNonCFMDorg;
        vNonCFMD = fliplr(vNonCFMD);
        % Calculate a and b-value for GR-law and distribution vNCum
        [~, ~, vSel, ~] = fMagToFitBValue(catalog, vFMD, fMc);
        if (length(catalog.Longitude(vSel)) >= 20)
            %[ fBValue, fStdDev, fAValue] =  calc_bmemag(catalog.Magnitude(vSel), binInterval);
            [fBValue, ~, fAValue] =  calc_bmemag(catalog.Magnitude(vSel), binInterval);
            % Normalize to time period
            vFMD(2,:)       = vFMD(2,:)./timespan; % ceil taken out
            vNonCFMD(2,:)   = vNonCFMD(2,:)./timespan; % ceil removed
            % Compute quantity of earthquakes by power law
            fMaxMagFMD  = max(vNonCFMD(1,:));
            fMinMagFMD  = min(vNonCFMD(1,:));
            vMstep      = fMinMagFMD:0.1:fMaxMagFMD;
            vNCum       = 10.^(fAValue-fBValue.*vMstep); % Cumulative number

            % Compute non-cumulative numbers vN
            fNCumTmp    = 10^(fAValue - fBValue * (fMaxMagFMD + 0.1));
            vNCumTmp    = [vNCum, fNCumTmp ];
            vN          = abs(diff(vNCumTmp));

            % Normalize vN
            vN = vN./timespan;
            % Data selection
            % mData = Non-cumulative FMD values from GR-law and original data
            mData = [vN' vNonCFMD'];
            vSel = (mData(:,2) >= fMc);
            mDataTest = mData(~vSel,:);
            mDataTmp = mData(vSel,:);
            % Check for zeros in observed data
            vSelCheck = (mDataTest(:,3) == 0);
            mDataTest = mDataTest(~vSelCheck,:);
            % Choices of normalization
            fNmax = mDataTmp(1,3); % Frequency of events in Mc bin

            if (~isempty(isempty(fNmax)) &&  ~isnan(fNmax) & fNmax ~= 0 & length(mDataTest(:,1)) > 4)
                mDataTest(:,3) = mDataTest(:,3)/fNmax; % Normalize datavalues for fitting with CDF
                % Move to M=0 to fit with lsq-algorithm
                fMinMagTmp = min(mDataTest(:,2));
                mDataTest(:,2) = mDataTest(:,2)-fMinMagTmp;
                % Curve fitting: Non cumulative part below Mc
                options = optimset('Display','off',...
                    'Tolfun'        , 1e-5,...
                    'TolX'          , 0.001,...
                    'MaxFunEvals'   , 1000,...
                    'MaxIter'       , 1000);
                [vX, resnorm, resid, exitflag, output, lambda, jacobian] = lsqcurvefit(...
                    @calc_normalCDF,[fSmu  fSSigma], mDataTest(:,2), mDataTest(:,3),[],[], options);
                mDataTest(:,1) = normcdf(mDataTest(:,2), vX(1), vX(2))*fNmax;
                if (length(mDataTest(:,2)) > length(vX(1,:)))
                    %% Confidence interval determination
                    % vPred : Predicted values of lognormal function
                    % vPred+-delta : 95% confidence level of true values
                    [vPred,delta] = nlpredci(@calc_normalCDF, mDataTest(:,2), vX, resid, jacobian);
                else
                    vPred = NaN;
                    delta = NaN;
                end % END: This section is due for errors produced with datasets less long than amount of parameters in vX
                % Results of fitting procedure
                mFitRes = [mFitRes; vX resnorm exitflag];
                % Move back to original magnitudes
                mDataTest(:,2) = mDataTest(:,2)+fMinMagTmp;
                % Set data together
                mDataTest(:,3) = mDataTest(:,3)*fNmax;
                mDataPred = [mDataTest; mDataTmp];
                % Denormalize to calculate probabilities
                mDataPred(:,1) = round(mDataPred(:,1).*timespan);
                mDataPred(:,3) = mDataPred(:,3).*timespan;
                vProb_ = calc_log10poisspdf2(mDataPred(:,3), mDataPred(:,1)); % Non-cumulative

                % Sum the probabilities
                fProbability    = (-1) * sum(vProb_,'omitnans');
                vProbability    = [vProbability; fProbability];
                % Move magnitude back
                mDataPred(:,2)  = mDataPred(:,2)+fMinMag;
                vMc             = [vMc; fMc];
                vABValue        = [vABValue; fAValue fBValue];

                % Keep values
                vDeltaBest  = [vDeltaBest; delta];
                vX_res      = [vX_res; vX resnorm exitflag];
                vNmaxBest   = [vNmaxBest; fNmax];

                % Keep best fitting model
                if (fProbability == min(vProbability))
                    vDeltaBest      = delta;
                    vPredBest       = [mDataTest(:,2) vPred*fNmax*timespan delta*fNmax*timespan]; % Gives back uncertainty
                    mDatPredBest    = [mDataPred];
                end
            else
                %disp('Not enough data');
                % Setting values
                fProbability = NaN;
                fMc = NaN;
                vX(1) = NaN;
                vX(2) = NaN;
                resnorm = NaN;
                exitflag = NaN;
                delta = NaN;
                vPred = [NaN NaN NaN];
                fNmax = NaN;
                fAValue = NaN;
                fBValue = NaN;
                vProbability = [vProbability; fProbability];
                vMc = [vMc; fMc];
                vX_res = [vX_res; vX resnorm exitflag];
                %             vDeltaBest = [vDeltaBest; NaN];
                %             vPredBest = [vPredBest; NaN NaN NaN];
                vNmaxBest = [vNmaxBest; fNmax];
                vABValue = [vABValue; fAValue fBValue];
            end % END of IF fNmax
        end % END of IF length(catalog.Longitude(vSel))

    end % END of FOR fMc
    % Result matrix
    mResult = [mResult; vProbability vMc vX_res vNmaxBest vABValue];

    % Find best estimate, excluding the case of mResult all NAN
    if  ~isempty(min(mResult)) && ~isnan(min(mResult(:,1)))
        vSel = find(min(mResult(:,1)) == mResult(:,1));
        fMc = min(mResult(vSel,2));
        %fMls = min(mResult(vSel,1));
        mu = min(mResult(vSel,3));
        fSigma = min(mResult(vSel,4));
        fAvalue = min(mResult(vSel,8));
        fBvalue = min(mResult(vSel,9));
    else
        fMc = NaN;
        %fMls = NaN;
        mu = NaN;
        fSigma = NaN;
        fAvalue = NaN;
        fBvalue = NaN;
    end
end

author: CelsoReyes