text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcMovie(Temp,Event,handles)
handles=guidata(handles.MU_matrix_display);

if length(handles.V.DimSize)==2
    warndlg('At least 3D matrix is needed for generating movie.');
    return;
end

dimFlag = ['[:' repmat(',:', [1 length(handles.V.DimSize)-1]) ']'];
dim = inputdlg(['Please specify dimension flag for generating movie. Use '':'' to include all content in the dimension; Use ''start:end'' to include parts of the content in the dimension.'],'Specify Dimension Flag',1,{dimFlag});
if isempty(dim)
    warndlg('Making movie was cancelled.');
    return;
end

movie = inputdlg({'Please specify the frame rate (fps).', 'Please specify the movie quality (lowest 0 - 100 highest)'},'About Movie',1,{'30' ,'75'});
if isempty(movie)
    warndlg('Making movie was cancelled.');
    return;
end

MergeM=get(handles.Matrix_name_edit,'String');
name = inputdlg(['Please input an name for the movie file.'],'Specify File Name',1,{[MergeM '_mov']});
if isempty(dim)
    warndlg('Making movie was cancelled.');
    return;
end

pause(0.1);
try
    eval(['TMatrix=handles.TMatrix(' dim{1}(2:end-1) ');']);
    TSize=size(TMatrix);
    TMatrix=TMatrix(:,:,:);
    dimSize=size(TMatrix);
    TMatrix=reshape(TMatrix,[dimSize(1) dimSize(2) 1 dimSize(3)]);

    if isempty(handles.V2.Foreground_matrix) % no overlay

        Colormap = colormap(handles.Matrix_display_axes);
        TMatrix=double(TMatrix);
        TMatrix(TMatrix<=double(handles.V.C_lower)) = double(handles.V.C_lower);
        TMatrix(TMatrix>=double(handles.V.C_upper)) = double(handles.V.C_upper);
        TMatrix=TMatrix-double(handles.V.C_lower);
        TMatrix=(TMatrix./double(handles.V.C_upper - handles.V.C_lower))*(length(Colormap(:,1))-1)+1;

        M=immovie(TMatrix,Colormap);

    else % if overlay
        BImage=double(repmat(TMatrix,[1 1 3 1]));
        BImage=max(double(handles.V.C_lower),min(BImage,double(handles.V.C_upper))); % Make sure BImage is in the range of Color Bound
        BImage=(BImage-double(handles.V.C_lower))./double((handles.V.C_upper-handles.V.C_lower));

        eval(['mask=double(handles.Mask(' dim{1}(2:end-1) '));']); % multi-dimensional mask??????
        TSize2=size(mask);
        if length(TSize)>length(TSize2) & length(TSize2)==2
            mask=repmat(mask,[1,1,prod(TSize(3:end))]);
        elseif length(TSize)>length(TSize2) & length(TSize2)==3
            mask=repmat(mask,[1,1,1,prod(TSize(4:end))]);
        end
        mask=mask(:,:,:);
        FMatrix=((max(double(handles.V2.F_lower),min(mask,double(handles.V2.F_upper)))-handles.V2.F_lower)/(handles.V2.F_upper-handles.V2.F_lower))*(64-1)+1; % Make sure FMatrix is in the range of Color Bound
        FImage=ind2rgb(round(FMatrix(:)),handles.V2.Color_map);
        if handles.V2.Include0 ~=0
            FImage(repmat(mask(:),[1 1 3])==0)=0;
        end
        FImage=reshape(FImage,[size(FMatrix) 3]);
        FImage=permute(FImage,[1 2 4 3]);
        RGBImage=BImage* handles.V2.Backgroud_F+FImage* handles.V2.Foregroud_F;
        RGBImage=max(0,min(RGBImage,1)); % Make sure RGBImage is in the range of 0 and 1

        M=immovie(RGBImage);
    end

    writerObj=VideoWriter(name{1});
    writerObj.FrameRate=str2double(movie{1});
    writerObj.Quality=str2double(movie{2});
    open(writerObj);
    writeVideo(writerObj,M);
    msgbox(['Movie file ' char(39) writerObj.Filename char(39) ' has been generated at ' char(39) writerObj.Path char(39) '.'],'Movie File');
    close(writerObj);
catch me
    error_msg{1,1}='Making movie failed! Error occured during movie generation process.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    return;
end

% update current display matrix
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function offobj = offLoader(filename)

% function offobj = offLoader(filename)

offobj = struct();
fid = fopen(filename, 'r');
words = fscanf(fid, 'OFF %d %d %d');
nV = words(1); nF = words(2); nE = words(3);
offobj.vmat = reshape(fscanf(fid, '%f', nV*3), 3, nV);
fstr = textscan(fid, '%s', nF, 'delimiter', '\n', 'MultipleDelimsAsOne', 1);
tfmat = zeros(nF*2,4);
nf3 = 0;
for i=1:nF
    words = sscanf(fstr{1}{i}, '%d');
    if words(1) == 3
        nf3 = nf3 + 1;
        tfmat(nf3, :) = words([4,2,3,4]);
    elseif words(1) == 4
        nf3 = nf3 + 1;
        tfmat(nf3, :) = words([4,2,3,4]);
        nf3 = nf3 + 1;
        tfmat(nf3, :) = words([5,2,4,5]);
    else
        error('size of face is not 3 or 4');
    end
end
offobj.fmat = tfmat(1:nf3, :)';
fclose(fid);

end

author: jianxiongxiao

repo: ProfXkit


text:
function y = dec2base(varargin)
% Works for empty input
% Allows inputs to be char
if isempty(varargin{1})
    y = '';
else
    if ischar(varargin{1})
        varargin{1} = double(varargin{1});
    end
    if nargin>1 && ischar(varargin{2})
        varargin{2} = double(varargin{2});
    end
    if nargin>2 && ischar(varargin{3})
        varargin{3} = double(varargin{3});
    end
    y = builtin('dec2base', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function varargout = tapas_physio_report_contrasts(varargin)
% This function reports all relevant F-contrast-maps for physIO-created regressors
% in a given GLM
%
%  - Section plots on an anatomical overlay are created and saved as
%    post-script file with one page for each contrasts
%
%  - Input parameters are specified via name/value pairs, e.g.
%
%
%   args = tapas_physio_report_contrasts(...
%                   'fileReport', 'physio.ps', ...
%                   'fileSpm', 'analysisFolder/SPM.mat', ...
%                   'filePhysIO', 'analysisFolder/physio.mat', ...
%                   'fileStructural', 'anatomyFolder/warpedAnatomy.nii')
%
% IN
%
%   Required parameters:
%
%                  fileReport: sets the path and file name for the postscript
%                              file your contrast report is printed to.
%                              A reasonable (and default) location is the
%                              single subject analysis folder so the file
%                              resides with the GLM (fileSPM)
%                              default: 'pathSpm/physio_report_contrasts.ps';
%              fileStructural: structural underlay for results,
%                              e.g. 'mean.nii'
%                     fileSpm: SPM.mat holding physiological regressors,
%                              e.g.'SPM.mat'
%                  filePhysIO:  mat-file where PhysIO-object was saved
%                               e.g. physio.mat
%
%   Optional Parameters:
%
%          namesPhysContrasts: cell Array of contrast names in design matrix
%                              e.g. {'All Phys', 'Cardiac', 'Respiratory',
%                               'Card X Resp Interation',
%                               'HeartRateVariability',
%                               'RespiratoryVolumePerTime', 'Movement'}
%
%      indReportPhysContrasts: vector of contrast indices to be reported
%                               e.g. [1:7] for all physiological contrasts
%     reportContrastThreshold: 0.001
%    reportContrastCorrection: 'none' or 'FWE'
%    reportContrastMax:        maximum value of contrast colormap
%                              to scale different contrasts with equal
%                              F-values (default: Inf, scales to max F of
%                              map)
%      reportContrastPosition: 'max' or [1,3] vector of crosshair position
%                               (in mm)
%               fovMillimeter: field of view in mm; set to 0 for full FOV
%         doPlotSliceParallel: if true, slices are plotted parallel to
%                              their slice acquisition direction
%                       model: physio.model-substructure holding
%                              model orders, i.e. .c .r .cr
%                              See also tapas_physio_new
%         titleGraphicsWindow: additional title prepended to contrast name
%                              in each plot
%       saveTable              true of false (default)
%                              prints screenshot of results table to
%                              ps-file as well
%
% OUT
%   args    structure of default and updated arguments used in this
%           function; the fields of args hold all possible options of this
%           function
%
% EXAMPLE
%   tapas_physio_report_contrasts
%
%   See also tapas_physio_overlay_contrasts

% Author: Lars Kasper
% Created: 2014-10-16
% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% START #MOD =============================================================

% general paths study
defaults.titleGraphicsWindow = '';
% PhysIO Toolbox code should be in same folder as this file
defaults.filePhysIO      = 'physio.mat';
defaults.fileReport      = 'physio_report_contrasts.ps'; % where contrast maps are saved
defaults.fileStructural  = 'mean.nii';
defaults.fileSpm         = 'SPM.mat';
defaults.drawCrosshair   = true;

defaults.namesPhysContrasts = tapas_physio_get_contrast_names_default();

% selection of physiological contrasts to be reported, corresponding to
% namesPhysContrasts order
defaults.indReportPhysContrasts = 1:9;

defaults.reportContrastThreshold     = 0.001; % 0.05; 0.001;
defaults.reportContrastCorrection    = 'none'; % 'FWE'; 'none';
defaults.reportContrastMax           = Inf;
defaults.reportContrastPosition      = 'max'; % 'max' to jump to max; or [x,y,z] in mm
defaults.fovMillimeter               = 0; %mm; choose 0 to plot whole FOV (bounding box)

% if true, voxel space (parallel to slices), not world space (with interpolation) is used
defaults.doPlotSliceParallel          = true;
defaults.saveTable                    = false;
physio                                = tapas_physio_new('RETROICOR');
defaults.model                        = physio.model; % holding number of physiological regressors

% END #MOD
%% ========================================================================

args = tapas_physio_propval(varargin, defaults);
tapas_physio_strip_fields(args);

spm('defaults', 'FMRI');

% make sure to use absolute paths from now on...
fileSpm = tapas_physio_filename2path(fileSpm);
fileStructural = tapas_physio_filename2path(fileStructural);
fileReport = tapas_physio_filename2path(fileReport);
filePhysIO = tapas_physio_filename2path(filePhysIO);

load(fileSpm);
nContrasts = numel(indReportPhysContrasts);

if ~exist(fileStructural, 'file')
    % take standard structural from SPM, if overlay file does not exist
    pathSpm = fileparts(which('spm'));
    fileStructural = fullfile(pathSpm, 'canonical', 'avg152T1.nii');
end

% if input file given, load PhysIO-object
if exist(filePhysIO, 'file')
    load(filePhysIO, 'physio');
    model = physio.model;
end

% Temporarily set window style to undocked, so that SPM opens as usual
tmpWindowStyle = get(0, 'DefaultFigureWindowStyle');
set(0, 'DefaultFigureWindowStyle', 'normal');

% create physiological contrasts that don't exist so far and can be created
% due to the model components included
SPM = tapas_physio_create_missing_physio_contrasts(SPM, model, namesPhysContrasts);

%% report contrasts with overlays

for c = 1:nContrasts
    iC = indReportPhysContrasts(c);
    idxContrasts(c) = tapas_physio_check_get_xcon_index(SPM, ...
        namesPhysContrasts{iC});
end

idxContrasts = idxContrasts(idxContrasts > 0);

% TODO: harmonize option labels in both functions
tapas_physio_overlay_contrasts( ...
    'idxContrasts', idxContrasts, ...
    'titleGraphicsWindow', titleGraphicsWindow, ...
    'fileReport', fileReport, ...
    'fileSpm', fileSpm, ...
    'fileStructural', fileStructural, ...
    'threshold', reportContrastThreshold, ...
    'correction', reportContrastCorrection, ...
    'doPlotSliceParallel', doPlotSliceParallel, ...
    'fovMillimeter', fovMillimeter, ...
    'position', reportContrastPosition, ...
    'drawCrosshair', drawCrosshair, ...
    'colorbarMax', reportContrastMax, ...
    'saveTable', saveTable);

if nargout
    varargout{1} = args;
end

author: translationalneuromodeling

repo: tapas


text:
function model = changeRxnBounds(model, rxnNameList, value, boundType)
% Changes upper or lower bounds of a reaction or a set of reactions
%
% USAGE:
%
%    model = changeRxnBounds(model, rxnNameList, value, boundType)
%
% INPUTS:
%    model:          COBRA model structure
%    rxnNameList:    List of reactions (cell array or string)
%    value:          Bound values.
%                    Can either be a vector or a single scalar value if the same
%                    bound value is to be assinged to all reactions
%
% OPTIONAL INPUT:
%    boundType:      'u' - upper, 'l' - lower, 'b' - both (Default = 'b')
%                    Bound type can either be a cell array of strings or a
%                    string with as many letters as there are reactions in  `rxnNameList`
%
% OUTPUT:
%    model:          COBRA model structure with modified reaction bounds
%
% .. Author: -  Markus Herrgard 4/21/06

if (nargin < 4)
    boundType = 'b';
end

if ((length(value) ~= length(rxnNameList) & length(value) > 1) | (length(boundType) ~= length(rxnNameList) & length(boundType) > 1))
   error('Inconsistent lenghts of arguments: rxnNameList, value & boundType');
end

rxnID = findRxnIDs(model,rxnNameList);

% Remove reactions that are not in the model
if (iscell(rxnNameList))
    missingRxns = rxnNameList(rxnID == 0);
    for i = 1:length(missingRxns)
        warning('Reaction %s not in model',missingRxns{i});
    end
    if (length(boundType) > 1)
        boundType = boundType(rxnID ~= 0);
    end
    if (length(value) > 1)
        value = value(rxnID ~= 0);
    end
    rxnID = rxnID(rxnID ~= 0);
else
    if rxnID==0
        warning('Reaction %s not in model',rxnNameList)
    end
end

if (isempty(rxnID))
    %rxnNameList was a cell with a single reaction not found
    %Warning would have already been triggered in line 39 above
    %So do nothing!
elseif (sum(rxnID) == 0)
    warning('Reaction %s not in model',rxnNameList);
else
    nRxns = length(rxnID);
    if (length(boundType) > 1)
        if (length(value) == 1)
            value = repmat(value,nRxns,1);
        end
        for i = 1:nRxns
            switch lower(boundType{i})
                case 'u'
                    model.ub(rxnID(i)) = value(i);
                case 'l'
                    model.lb(rxnID(i)) = value(i);
                case 'b'
                    model.lb(rxnID(i)) = value(i);
                    model.ub(rxnID(i)) = value(i);
            end
        end
    else
        switch lower(boundType)
            case 'u'
                model.ub(rxnID) = value;
            case 'l'
                model.lb(rxnID) = value;
            case 'b'
                model.lb(rxnID) = value;
                model.ub(rxnID) = value;
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function a=d_logit(p)
%DREX subfunction
%Written by Issam El Naqa 2003-2005
%Extracted for generalized use 2005, AJH
%
% Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
% This file is part of the Dose Response Explorer System (DREES).
%
% DREES development has been led by:  Issam El Naqa, Aditya Apte, Gita Suneja, and Joseph O. Deasy.
%
% DREES has been financially supported by the US National Institutes of Health under multiple grants.
%
% DREES is distributed under the terms of the Lesser GNU Public License.
%
%     This version of DREES 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.
%
% DREES 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 DREES.  If not, see <http://www.gnu.org/licenses/>.

a = 1 ./ max(eps, (p .* (1-p)));

return

author: mvallieres

repo: radiomics


text:
function c = diff(c)

% DIFF differentiate each trace.
%
% C = DIFF(C) differentiate each trace through a call to WAVEFORM/DIFF.
% See WAVEFORM/DIFF for details

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

c.W = diff(c.W);

author: geoscience-community-codes

repo: GISMO


text:
function hdr = read_ctf_res4(fname)

% READ_CTF_RES4 reads the header in RES4 format from a CTF dataset
%
% hdr = read_ctf_res4(filename)
%
% Due to a non-dislosure agreement between CTF and the F.C. Donders Centre,
% the allwoed use of this function is limited. Do not distribute this function.

% This program is provided to users of CTF MEG systems as a courtesy only. Please
% do not redistribute it without permission from CTF Systems Inc.
% This program has no warranty whatsoever.

% Author(s): Jim McKay November 1999
% Last revision: Jim McKay
% Copyright (c) 1999-2000 CTF Systems Inc. All Rights Reserved.
%
% modifications Copyright (C) 2002, Ole Jensen
% modifications Copyright (C) 2003, Robert Oostenveld
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read header information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid = fopen(fname,'r','b');

% Check if header file exist
if fid == -1
    errMsg = strcat('Could not open header file:',fname);
    error(errMsg);
end

% first 8 bytes contain filetype

% Check is fileformat is correct
r_head=setstr(fread(fid,8,'char'))';
if (strcmp(r_head(1,1:7),'MEG41RS')==0),
  fclose(fid)
  errMsg = strcat('Resource file is not in CTF MEG4 format in file:',fname)
  error(errMsg);
end %if

% Read the initial parameters
appName       = setstr(fread(fid,256,'char'))' ;
dataOrigin    = setstr(fread(fid,256,'char'))' ;
dataDescrip   = setstr(fread(fid,256,'char'))' ;
no_trial_avgd = fread(fid,1,'int16')          ;
data_time     = setstr(fread(fid,255,'char'))';
data_date     = setstr(fread(fid,255,'char'))';

fseek(fid,1288,'bof');
% Read the general recording parameters
no_samples  = fread(fid,1,'int32');
no_channels = fread(fid,1,'int16');
fseek(fid,2,'cof');			% hole of 2 bytes due to improper alignment
sample_rate = fread(fid,1,'double');
epoch       = fread(fid,1,'double');
no_trials   = fread(fid,1,'int16');
fseek(fid,2,'cof');			% hole of 2 bytes due to improper alignment
preTrigpts=fread(fid,1,'int32');

fseek(fid,1360,'bof');
% read in the meg4Filesetup structure
run_name     = setstr(fread(fid,32,'char')');
run_title    = setstr(fread(fid,256,'char')');
instruments  = setstr(fread(fid,32,'char')');
coll_desc    = setstr(fread(fid,32,'char')');
subj_id      = setstr(fread(fid,32,'char')');
operator     = setstr(fread(fid,32,'char')') ;
sensFilename = setstr(fread(fid,60,'char')') ;

fseek(fid,1839,'bof');
% Read in the run description length
rd_len=fread(fid,1,'uint8');
% Go to the run description and read it in
fseek(fid,1844,'bof');
run_desc=setstr(fread(fid,rd_len,'char')');

% read in the filter information
temp=fread(fid,2,'uint8');
num_filt=temp(2);
for fi=0:(num_filt-1),
  filt_info=fread(fid,18,'uint8');
  num_fparm=filt_info(18);
  if num_fparm ~= 0,
    filt_parm=fread(fid,8*num_fparm,'uint8');
  end % if
end % for fi

% Read in the channel names
for i=1:no_channels,
    temp=fread(fid,32,'char')';
    temp(find(temp<32 )) = ' ';		% remove non-printable characters
    temp(find(temp>126)) = ' ';		% remove non-printable characters
    endstr = findstr(temp, '-'); temp(endstr:end) = ' ';	% cut off at '-'
    endstr = findstr(temp, ' '); temp(endstr:end) = ' ';	% cut off at ' '
    chan_name(i,:) = char(temp);		% as char array
    chan_label{i}  = deblank(char(temp));	% as cell array
end %for

% pre-allocate some memory space
sensGain = zeros([no_channels,1]);
qGain    = zeros([no_channels,1]);
ioGain   = zeros([no_channels,1]);
sensType = zeros([no_channels,1]);

% Read in the sensor information
fp = ftell(fid);
for chan=1:no_channels,
    fread(fid,1,'uint8');			% Read and ignore 1 byte from enum
    sensType(chan)=fread(fid,1,'uint8');	% Read sensor type
    fread(fid,2,'uint8');			% Read and ignore originalRunNum
    fread(fid,4,'uint8');			% Read and ignore coilShape
    sensGain(chan)=fread(fid,1,'double');	% Read sensor gain in Phi0/Tesla
    qGain(chan)=fread(fid,1,'double');		% Read qxx gain (usually 2^20 for Q20)
    ioGain(chan)=fread(fid,1,'double');		% Read i/o gain of special sensors (usually 1.0)
    ioOffset(chan)=fread(fid,1,'double');
    numCoils(chan)=fread(fid,1,'int16');
    grad_order_no(chan)=fread(fid,1,'int16');
    fread(fid,4,'char');

    % read the coil positions and orientations
    for i=1:8
      Chan(chan).coil(i).pos = fread(fid,3,'double')';
      fread(fid,1,'double');
      Chan(chan).coil(i).ori = fread(fid,3,'double')';
      fread(fid,3,'double');
    end

    % read the coil positions and orientations in head coordinates(?)
    for i=1:8
      Chan(chan).coilHC(i).pos = fread(fid,3,'double')';
      fread(fid,1,'double');
      Chan(chan).coilHC(i).ori = fread(fid,3,'double')';
      fread(fid,3,'double');
    end

    % jump to the next sensor info record
    fseek(fid, fp+chan*1328, 'bof');
end % for chan

% close the header file
fclose(fid);

% determine the different channel types
rowMEG  = [];
rowEEG  = [];
rowTRIG = [];
rowREF  = [];
rowALL  = 1:no_channels;
for k=1:no_channels
    if strcmp(char(chan_name(k,1)),'M')
        rowMEG = [rowMEG k];
    end
    if strcmp(char(chan_name(k,1:3)),'EEG')
        rowEEG = [rowEEG k];
    end
    if strcmp(char(chan_name(k,1:4)),'STIM')
        rowTRIG = [rowTRIG k];
    end
    if chan_name(k,1)=='B' | chan_name(k,1)=='G' | chan_name(k,1)=='P' | chan_name(k,1)=='Q' | chan_name(k,1)=='R'
        rowREF = [rowREF k];
    end
end

% assign all the variables that should be outputted as header information
hdr.Fs          = sample_rate;
hdr.nChans      = no_channels;
hdr.nSamples    = no_samples;
hdr.nSamplesPre = preTrigpts;
hdr.timeVec     = (1:no_samples)/sample_rate - preTrigpts/sample_rate;
hdr.nTrials     = no_trials;
hdr.gainV       = ioGain./(qGain.*sensGain);
hdr.label       = chan_label(:);
hdr.nameALL     = chan_name;
hdr.rowMEG      = rowMEG;
hdr.rowEEG      = rowEEG;
hdr.rowTRIG     = rowTRIG;
hdr.rowREF      = rowREF;
hdr.Chan        = Chan;
% hdr.nameEEG     = [];
% hdr.nameMEG     = [];
% hdr.nameEOG     = [];
% hdr.trigV       = [];
% hdr.SwapData    = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% construct the gradiometer system in DEWAR coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% combine the bottom and top coil of the MEG channels into hardware gradiometers
numMEG = length(hdr.rowMEG);
for i=1:numMEG
  coil1_pos(i,:) = hdr.Chan(hdr.rowMEG(i)).coil(1).pos;
  coil1_ori(i,:) = hdr.Chan(hdr.rowMEG(i)).coil(1).ori;
  coil2_pos(i,:) = hdr.Chan(hdr.rowMEG(i)).coil(2).pos;
  coil2_ori(i,:) = hdr.Chan(hdr.rowMEG(i)).coil(2).ori;
end

% apparently, some coils are oriented to the wrong side
tmp = coil2_pos - coil1_pos;
sel = find(dot(coil1_ori, tmp, 2)<1);
coil1_ori(sel,:) = -coil1_ori(sel,:);
coil2_ori(sel,:) = -coil2_ori(sel,:);

gradDEWAR.pnt = [coil1_pos; coil2_pos];
gradDEWAR.ori = [coil1_ori; coil2_ori];
gradDEWAR.tra = [eye(numMEG) -eye(numMEG)];
gradDEWAR.label = hdr.label(rowMEG);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% construct the gradiometer system in HEAD coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% combine the bottom and top coil of the MEG channels into hardware gradiometers
numMEG = length(hdr.rowMEG);
for i=1:numMEG
  coil1_pos(i,:) = hdr.Chan(hdr.rowMEG(i)).coilHC(1).pos;
  coil1_ori(i,:) = hdr.Chan(hdr.rowMEG(i)).coilHC(1).ori;
  coil2_pos(i,:) = hdr.Chan(hdr.rowMEG(i)).coilHC(2).pos;
  coil2_ori(i,:) = hdr.Chan(hdr.rowMEG(i)).coilHC(2).ori;
end

% apparently, some coils are oriented to the wrong side
tmp = coil2_pos - coil1_pos;
sel = find(dot(coil1_ori, tmp, 2)<1);
coil1_ori(sel,:) = -coil1_ori(sel,:);
coil2_ori(sel,:) = -coil2_ori(sel,:);

gradHC.pnt = [coil1_pos; coil2_pos];
gradHC.ori = [coil1_ori; coil2_ori];
gradHC.tra = [eye(numMEG) eye(numMEG)];
gradHC.label = hdr.label(rowMEG);

hdr.grad      = gradHC;		% default is in head-coordinates
hdr.grad.unit = 'cm';		% this is default in all CTF software

author: PatternRecognition

repo: OpenBMI


text:
classdef ENSMOEAD < ALGORITHM
% <multi/many> <real/integer>
% Ensemble of different neighborhood sizes based MOEA/D
% NS --- 25:25:100 --- Set of neighborhood sizes
% LP ---        50 --- Learning period

%------------------------------- Reference --------------------------------
% S. Zhao, P. N. Suganthan, and Q. Zhang, Decomposition-based multi-
% objective evolutionary algorithm with an ensemble of neighborhood sizes,
% IEEE Transactions on Evolutionary Computation, 2012, 16(3): 442-446.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [NS,LP] = Algorithm.ParameterSet(25:25:100,50);

            %% Generate the weight vectors
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            % Maximum number of solutions replaced by each offspring
            nr = ceil(Problem.N/100);

            %% Detect all the neighbours of each solution
            B = pdist2(W,W);
            [~,B] = sort(B,2);

            %% Generate random population
            Population = Problem.Initialization();
            Z          = min(Population.objs,[],1);
            % Utility for each subproblem
            Pi = ones(Problem.N,1);
            % Old Tchebycheff function value of each solution on its subproblem
            oldObj = max(abs((Population.objs-repmat(Z,Problem.N,1)).*W),[],2);

            %% Optimization
            p           = ones(1,length(NS))./length(NS);
            FEs         = zeros(1,length(NS));
            FEs_success = zeros(1,length(NS));
            while Algorithm.NotTerminated(Population)
                % Select neighborhood size for each subproblem
                ns = RouletteWheelSelection(Problem.N,1./p);

                % Apply MOEA/D-DRA for one generation
                for subgeneration = 1 : 5
                    % Choose I
                    Bounday = find(sum(W<1e-3,2)==Problem.M-1)';
                    I = [Bounday,TournamentSelection(10,floor(Problem.N/5)-length(Bounday),-Pi)];

                    % For each solution in I
                    for i = I
                        % Choose the parents
                        if rand < 0.9
                            P = B(i,randperm(NS(ns(i))));
                        else
                            P = randperm(Problem.N);
                        end

                        % Generate an offspring
                        Offspring = OperatorDE(Problem,Population(i),Population(P(1)),Population(P(2)));

                        % Update the ideal point
                        Z = min(Z,Offspring.obj);

                        % Update the solutions in P by Tchebycheff approach
                        g_old   = max(abs(Population(P).objs-repmat(Z,length(P),1)).*W(P,:),[],2);
                        g_new   = max(repmat(abs(Offspring.obj-Z),length(P),1).*W(P,:),[],2);
                        replace = find(g_old>=g_new,nr);
                        Population(P(replace)) = Offspring;
                        if ~isempty(replace)
                            FEs_success(ns(i)) = FEs_success(ns(i)) + 1;
                        end
                        FEs(ns(i)) = FEs(ns(i)) + 1;
                    end
                end
                if ~mod(ceil(Problem.FE/Problem.N),10)
                    % Update Pi for each solution
                    newObj    = max(abs((Population.objs-repmat(Z,Problem.N,1)).*W),[],2);
                    DELTA     = (oldObj-newObj)./oldObj;
                    Temp      = DELTA < 0.001;
                    Pi(~Temp) = 1;
                    Pi(Temp)  = (0.95+0.05*DELTA(Temp)/0.001).*Pi(Temp);
                    oldObj    = newObj;
                end
                if ~mod(ceil(Problem.FE/Problem.N),LP)
                    % Update the probability of choosing each neighborhood size
                    R           = FEs_success./FEs;
                    p           = R./sum(R);
                    FEs         = zeros(1,length(NS));
                    FEs_success = zeros(1,length(NS));
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function createHdf5File(filename, datasetName, inputData, varargin)
	% Creates an HDF5 file at filename under given datasetName hierarchy and saves inputData.
	% Biafra Ahanonu
	% started: 2014.01.07
	% inputs
		%
	% outputs
		%

	% changelog
		% 2019.03.25 [17:17:49] - Add support for custom user HDF5 chunking as opposed to previous automatic chunking
		% 2019.08.20 [11:38:54] - Added additional support for more data types.
        % 2021.02.02 [13:15:11] - Close space_id, dset_id, and fid with low-level HDF5 functions before appending data with hdf5write to avoid read/write issues.
        % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
        % 2021.12.08 [22:34:26] - Handle cases in which dataDimsChunkCopy that is input is too large for given matrix being saved to an HDF5 file.
	% TODO
		%

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

	%========================
	% Int: Defines gzip compression level (0-9). 0 = no compression, 9 = most compression.
	options.deflateLevel = 1;
	% Int: chunk size in [x y z] of the dataset, leave empty for auto chunking
	options.dataDimsChunkCopy = [];
	options.addInfo = [];
	% Char array: cell array of strings matching addInfo, e.g. '/movie/processingSettings'
	options.addInfoName = [];
	% get options
	options = getOptions(options,varargin);
	% % unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	%     eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	display(['creating: ' filename])
	% Create the HDF5 file
	fcpl_id = H5P.create('H5P_FILE_CREATE');
	fapl_id = H5P.create('H5P_FILE_ACCESS');

	fid = H5F.create(filename, 'H5F_ACC_TRUNC', fcpl_id, fapl_id);

	% Create the Space for the Dataset
	initDims = size(inputData);
	h5_initDims = fliplr(initDims);
	maxDims = [initDims(1) initDims(2) -1];
	h5_maxDims = fliplr(maxDims);
	space_id = H5S.create_simple(3, h5_initDims, h5_maxDims);

	% Create the Dataset
	% datasetName = '1';
	dcpl_id = H5P.create('H5P_DATASET_CREATE');

	if isempty(options.dataDimsChunkCopy)
		% set the last dimension to 1 for chunking
		chunkSize = [initDims(1) initDims(2) 1];
	else
		chunkSize = options.dataDimsChunkCopy;
	end

	% If chunk size is larger than the movie, reduce to size of the chunk.
	if any([chunkSize(1)>initDims(1), chunkSize(2)>initDims(2)])
		disp('Chunk size too large, using movie dimensions for H5 chunk size.')
		chunkSize = [initDims(1) initDims(2) 1];
	end

	h5_chunkSize = fliplr(chunkSize);

	H5P.set_chunk(dcpl_id, h5_chunkSize);
	fprintf('Set compression level to %d of 9\n',options.deflateLevel);
	H5P.set_deflate(dcpl_id,options.deflateLevel);

	inputClass = class(inputData);
	switch inputClass
		case 'single'
			dsetType_id = H5T.copy('H5T_IEEE_F32LE');
		case 'double'
			dsetType_id = H5T.copy('H5T_NATIVE_DOUBLE');
		case 'uint16'
			dsetType_id = H5T.copy('H5T_NATIVE_UINT16');
		case 'int16'
			dsetType_id = H5T.copy('H5T_NATIVE_INT16');
		case 'uint8'
			dsetType_id = H5T.copy('H5T_NATIVE_UINT8');
		case 'int8'
			dsetType_id = H5T.copy('H5T_NATIVE_INT8');
		otherwise
			disp(['Data type ' inputClass ' not supported.'])
			return;
			% body
	end

	dset_id = H5D.create(fid, datasetName, dsetType_id, space_id, dcpl_id);

	% Initial Data to Write
	% rowDim = initDims(1); colDim = initDims(2);
	% initDataToWrite = rand(initDims);

	% Write the initial data
	H5D.write(dset_id, 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', inputData);

    % Close the opened identifiers
	H5S.close(space_id);
	H5D.close(dset_id);
	H5F.close(fid);

	% Append movie relevant information to HDF5 file
	% if strcmp(options.writeMode,'new')
		% hdf5write(filename,'/movie/info/dimensions',dataDims,'WriteMode','append');
		currentDateTimeStr = datestr(now,'yyyymmdd_HHMM','local');
		hdf5write(filename,'/movie/info/date',currentDateTimeStr,'WriteMode','append');
		hdf5write(filename,'/movie/info/savePath',filename,'WriteMode','append');
		hdf5write(filename,'/movie/info/Deflate',options.deflateLevel,'WriteMode','append');
	% end

	% add information about data to HDF5 file
	if ~isempty(options.addInfo)
		if ~iscell(options.addInfo)
			options.addInfo = {options.addInfo};
			options.addInfoName = {options.addInfoName};
		end
		addInfoLen = length(options.addInfo);
		for addInfoStructNo = 1:addInfoLen
			thisAddInfo = options.addInfo{addInfoStructNo};
			infoList = fieldnames(thisAddInfo);
			nInfo = length(infoList);
			for fieldNameNo = 1:nInfo
				thisField = infoList{fieldNameNo};
				hdf5write(filename,[options.addInfoName{addInfoStructNo} '/' thisField],thisAddInfo.(thisField),'WriteMode','append');
			end
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function varargout = sech(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('bell-shape','callback');
        operator.derivative = @(x)(-tanh(x).*sech(x));
        operator.stationary = [0 1];
        operator.inflection = [-inf 1 -0.881493604 -1 0.881493604 1];
        operator.range = [0 1];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause

function [lbls, lblMap]=resupervise(umap, inData, newSubsetIdxs)
supervisors=umap.supervisors;
R=size(inData, 1);
subsetIds=unique(supervisors.labels);
subsetIds=subsetIds(subsetIds~=0);
mx=max(subsetIds)+1;
nSubsets=length(subsetIds);
a=zeros(nSubsets, R);
lblMap=java.util.TreeMap;
key=num2str(mx);
lblMap.put(java.lang.String(key), 'Previously unsupervised');
key=[key '.color'];
perc=.4;
clr=num2str([perc perc perc]);
lblMap.put(key, clr);
for i=1:nSubsets
    gid=subsetIds(i);
    updateMap(gid)
    r=umap.raw_data(supervisors.labels==gid,:);
    means_=mean(r);
    stds_=std(r);
    B=(abs(inData-means_(1,:)))./stds_(1,:);
    a(i,:)=sum(B,2);
end
lbls=zeros(R, 1);
[~, I]=min(a);
%lbls(unknownIdxs)=0;%mx;%0-subsetIds(I(unknownIdxs));
lbls(~newSubsetIdxs)=subsetIds(I(~newSubsetIdxs));
lbls(newSubsetIdxs)=mx;
disp([unique(lbls)'; UmapUtil.discreteCount(lbls, unique(lbls))]);
disp([unique(supervisors.labels)'; UmapUtil.discreteCount(supervisors.labels, unique(supervisors.labels))]);

    function updateMap(gid)
        key=num2str(gid);
        jKey=java.lang.String(key);
        name=supervisors.labelMap.get(jKey);
        lblMap.put(jKey, name);
        key=[key '.color'];
        lblMap.put(key, supervisors.labelMap.get(key));

    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Updated_Grid] = compute_transformation(i,Grid,H,E,Theta)
if (Grid.x(i) < E)
    Grid.y_s(i) = 0;
    Grid.h(i) = H;
else
    Grid.y_s(i) = (-Grid.x(i)*tan(Theta)) + (E*tan(Theta));
    Grid.h(i) = H + (Grid.x(i)*tan(Theta)) - (E*tan(Theta));
end
for j = 1:401,
    if (Grid.x(i) < E)
        Grid.m1(j,i) = 0;
    else
        Grid.m1(j,i) = (tan(Theta)/Grid.h(i)) - (Grid.y_t(j)*(tan(Theta)/Grid.h(i)));
    end
    Grid.y(j,i) = (Grid.y_t(j)*Grid.h(i)) + Grid.y_s(i);
    Grid.m2(j,i) = 1/Grid.h(i);
end
Updated_Grid = Grid;
end

author: Sable

repo: mcbench-benchmarks


text:
function out = getChebfunType(~)
%GETCHEBFUNTYPE   Returns SPHEREFUN for a SPINOPSPHERE.
%
% See also SPINOPSPHERE.

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

out = @spherefun;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function abram2_values_test ( )

%*****************************************************************************80
%
%% ABRAM2_VALUES_TEST demonstrates the use of ABRAM2_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ABRAM2_VALUES_TEST:\n' );
  fprintf ( 1, '  ABRAM2_VALUES stores values of \n' );
  fprintf ( 1, '  the Abramowitz function of order 2.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = abram2_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function TF = isax(h)
%ISAX Check if input is an axes handle.
% Usage:
%   TF = isax(h)
%
% Args:
%   h:
%
% See also: isfig, isgraphics

TF = isa(h,'matlab.graphics.axis.Axes');

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function lebesgue_test07 ( )

%*****************************************************************************80
%
%% LEBESGUE_TEST07 looks at Equidistant3 points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEBESGUE_TEST07:\n' );
  fprintf ( 1, '  Analyze Equidistant3 points.\n' );

  n_max = 11;
  l = zeros ( n_max, 1 );
  xfun = linspace ( -1.0, +1.0, 501 );
  xfun = xfun';

  for n = 1 : n_max
    x = equidistant3 ( n );
    l(n) = lebesgue_constant ( n, x, xfun );
  end

  r8vec_print ( n_max, l, '  Equidistant3 Lebesgue constants for N = 1 to 11:' )
%
%  Examine one case more closely.
%
  n = 11;
  x = equidistant3 ( n );
  r8vec_print ( n, x, '  Equidistant3 points for N = 11' );

  label = 'Equidistant3 points for N = 11';
  filename = 'equidistant3.png';
  lebesgue_plot ( n, x, xfun, label, filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved as "%s"\n', filename );

  return
end


text:
function [I, J] = ind2sub4up(IND)
%IND2SUB4UP Subscripts from linear index for upper triangular matrix (only
%elements above diagonal)
%   IND2SUB4UP determines the equivalent subscript values corresponding to
%   a given single index into a 2D upper triangular matrix, excluded all
%   elements over the diagonal.
%
%   [I, J] = IND2SUB4UP(IND) returns vectors I and J containing equivalent
%   row and column subscripts corresponding to the index vector IND.
%
%   Let ind be a vector of indexes for entries of some upper triangular
%   matrix. The entries are selected vertically so that:
%
%       ind = 1                is associated to entry      (1, 2)
%       ind = 2                is associated to entry      (1, 3)
%       ind = 3                is associated to entry      (2, 3)
%       ind = 4                is associated to entry      (1, 4)
%       ...
%       ind = N * (N - 1) / 2  is associated to entry      (N - 1, N)
%
% % ======================================================================
%
%    EXAMPLE
%
%       % Note that if
%             A = rand(10);
%       % and
%             b = A(find(triu(A, 1)));
%       % then, given indices
%             IND = [1:45];
%       % for vector b, these are equivalent to subscripts
%             [I, J] = ind2sub4up(IND);
%       % for matrix A. In fact:
%             all(A(sub2ind(size(A), I, J)) == b(IND))
%
%       %    ans =
%       %           1
%
%       % This is obtained without even knowing about size(A)
%
% % ======================================================================
%
%   See also SUB2IND, IND2SUB, FIND.
%
% % ======================================================================
%
%-*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-*%
%                                                                                %
%            Author: Liber Eleutherios                                             %
%            E-Mail: libereleutherios@gmail.com                             %
%            Date: 1 May 2010                                                       %
%                                                                                %
%-*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-*%
%
% % ======================================================================
%

% Check input
ctrl1 = isnumeric(IND) & isreal(IND);
if ctrl1
  IND = ceil(IND(:));
  ctrl2 = ~any(isnan(IND)) & ~any(isinf(IND)) & all(IND > 0);
  if ~ctrl2
  error('Check indexes: they need to be positive integers!')
  end
else
  error('Check indexes: they need to be positive integers!')
end

J = round(floor(-.5 + .5 * sqrt(1 + 8 * (IND - 1))) + 2);
I = round(J .* (3 - J) / 2 + IND - 1);

repo: gspbox


text:
function Y = slmulvec(X, v, d)
%SLMULVEC multiplies a vector to columns or rows of a matrix
%
% $ Syntax $
%   - Y = slmulvec(X, v, d)
%   - Y = slmulvec(X, v)
%
% $ Arguments $
%   - X:        The original matrix
%   - v:        The addend vector
%   - d:        The dimension along which the vector is to add
%   - Y:        The resultant matrix
%
% $ Description $
%   - Y = slmulvec(X, v, d) selects the most efficienct way to multiple a
%     vector v to every column/row of X. If d == 1, then v should be
%     a column vector, and is multiplied to each column of X, if d == 2,
%     then v should be a row vector, and is multiplied to each row of X.
%
%   - Y = slmulvec(X, v) will automatically determine d according to
%     the shape of v.
%
% $ Remarks $
%   - The implementation simply wraps the mex function vecop_core.
%
% $ History $
%   - Created by Dahua Lin, on Sep 10, 2006
%

if nargin < 3
    if size(v, 2) == 1
        d = 1;
    else
        d = 2;
    end
end

Y = vecop_core(X, v, d, 2);  % 2 is the opcode of multiplication in vecop_core

author: lmthang


text:
function c=ref_rdgtii(f,g,a,M)
%REF_RDGTII  Reference Real DGT type II
%   Usage:  c=ref_rdgt(f,g,a,M);
%
%   Linear algebra version of the algorithm. Create big matrix
%   containing all the basis functions and multiply with the transpose.

L=size(f,1);

b=L/M;
N=L/a;

Mhalf=ceil(M/2);

F=zeros(L,M*N);

l=(0:L-1).';

for n=0:N-1

  % Do the unmodulated coefficient.
  F(:,M*n+1)=circshift(g,n*a+floor(a/2));

  for m=1:Mhalf-1
    F(:,M*n+2*m)=sqrt(2)*cos(2*pi*m*(l+.5)/M).*circshift(g,n*a+floor(a/2));

    F(:,M*n+2*m+1)=sqrt(2)*sin(2*pi*m*(l+.5)/M).*circshift(g,n*a+floor(a/2));

  end;

  if mod(M,2)==0
    F(:,M*(n+1))=cos(pi*l).*circshift(g,n*a+floor(a/2));
  end;

end;

% dot-transpose will work because F is real.
c=F.'*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function epochJ=JulDate2JulEpoch(Jul1,Jul2)
%%JULDATE2JULEPOCH Convert two-part Julian dates given in a uniform
%                  timescale, such as terrestrial time (TT) date to a
%                  Julian epoch in the same scale. The timescale
%                  should be uniform (i.e. not UTC).
%
%INPUTS: Jul1, Jul2  Matrices of two parts of a Julian date given in the
%                    same timescale as the epoch The units of the date are
%                    days. The full date is the sum of both terms. The
%                    entries in the matrices correspond to different
%                    dates to be converted.
%
%OUTPUTS: epochJ A matrix of Julian epochs as fractional years in the same
%                timescale as the input.
%
%A Julian epoch is a factional year number denominated in terms of a
%year of exactly 365.25 days in TT. For example 2000.5 is a Julian epoch
%in the middle of the year 2000.
%
%This is a wrapper for the function iauEpj in the International
%Astronomical Union's (IAU) Standard's of Fundamental Astronomy library.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%epochJ=JulDate2JulEpoch(Jul1,Jul2);
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

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:
close all;
clearvars;
clc;
% The function
f=@(x,y) 2*x.^2 + 4*y.^2;;
% Construct a mesh of (x,y) values
[X,Y]=meshgrid(-3:0.1:3);
% compute f on all these (x, y) values
Z=f(X,Y);
% draw a contour plot
levels = [2:2:20];
[c, h] = contour3(X,Y,Z);
grid on;

author: indigits

repo: sparse-plex


text:
X=randn(64,200000);
A=sprand(200,200000,0.05);

tic
XAt=mexCalcXAt(X,A);
t=toc;
fprintf('mex-file time: %fs\n',t);

tic
XAt2=X*A';
t=toc;
fprintf('mex-file time: %fs\n',t);

sum((XAt(:)-XAt2(:)).^2)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function ulimit = SlopeLimitN(u);

% function ulimit = SlopeLimitN(u);
% Purpose: Apply slopelimiter (Pi^N) to u assuming u an N'th order polynomial

Globals1D;

% Compute cell averages
uh = invV*u; uh(2:Np,:)=0; uavg = V*uh; v = uavg(1,:);

% Apply slope limiter as needed.
ulimit = u; eps0=1.0e-8;

% find end values of each element
ue1 = u(1,:); ue2 = u(end,:);

% find cell averages
vk = v; vkm1 = [v(1),v(1:K-1)]; vkp1 = [v(2:K),v(K)];

% Apply reconstruction to find elements in need of limiting
ve1 = vk - minmod([(vk-ue1);vk-vkm1;vkp1-vk]);
ve2 = vk + minmod([(ue2-vk);vk-vkm1;vkp1-vk]);
ids = find(abs(ve1-ue1)>eps0 | abs(ve2-ue2)>eps0);

% Check to see if any elements require limiting
if(~isempty(ids))
  % create piecewise linear solution for limiting on specified elements
  uhl = invV*u(:,ids); uhl(3:Np,:)=0; ul = V*uhl;

  % apply slope limiter to selected elements
  ulimit(:,ids) = SlopeLimitLin(ul,x(:,ids),vkm1(ids),vk(ids),vkp1(ids));
end
return;

repo: Aero-matlab


text:
function eegplugin_miclust( fig, try_strings, catch_strings);
plotmenu = findobj(fig, 'tag', 'plot');
uimenu(plotmenu, 'label', 'Cluster dataset ICs', 'callback', ...
    [try_strings.check_ica '[EEG LASTCOM]= pop_miclust(EEG);' catch_strings.store_and_hist '[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);' ]);


text:
function illustrate_ppp(ResultsAllCellLines,mets,path,samples,label,fonts,tol)
% This function generates and saves heatmaps for the results of the function
% performPPP for all sample models.
%
% USAGE:
%
%     illustrate_ppp(ResultsAllCellLines, mets, path, samples, label, fonts, tol)
%
% INPUTS:
%    ResultsAllCellLines:     Result structure
%    mets:                    Metabolites that were tested in the phase plane analysis
%    step_size:               Step size of each metabololite tested
%    path:                    Path where output is saved
%    samples:                 Names of conditions
%    label:                   Defining label of X-axis, y-axis and z-axis, e.g., {`Glucose uptake (fmol/cell/hr)`; `Oxygen uptake
%                             (fmol/cell/hr)`; `Growth rate (hr-1)`} The z-axis is the growth rate which is color coded
%    fonts:                   Font size for labels on heatmap
%    tol:                     Fluxes are considered zero if they are below the tol
%
% Individual pdf files showing the result of the PPP are being saved automatically for each condition.
%
% .. Author: - Maike K. Aurich 19/02/15

k=1;

for i=1:length(samples)

    for p=1:size(mets,2)
              if contains(mets{p,1},'[')
         name  = ['phasePlane_' strtok(mets{1,p}, '[') '_' strtok(mets{2,p}, '[')];
        else
        name  = ['phasePlane_' strtok(mets{1,p}, '(') '_' strtok(mets{2,p}, '(')];
              end
        %load condition specific performPPP result
        growthRates = ResultsAllCellLines.(samples{i}).(name).growthRates;
        bounds = ResultsAllCellLines.(samples{k}).(name).bounds;

        %chop off the additional entries in bounds
        ivalues=bounds(1:size(growthRates,1),1);
        jvalues=bounds(1:size(growthRates,2),2);

        % to mark the reaction bounds in the pruned model onto the phase plane
        modelPruned = eval(['ResultsAllCellLines.' samples{i} '.modelPruned']);
        lb = modelPruned.lb(find(ismember(modelPruned.rxns, mets{1,p})));
        ub = modelPruned.ub(find(ismember(modelPruned.rxns, mets{1,p})));
        clear modelPruned

        %make title
        name = regexprep(samples{i},'_','-');
        title_name = ['Objective values under variation of ' regexprep(mets{1,p},'_','-') ' and ' regexprep(mets{2,p},'_','-') ' in ' name];
       % define max of z-axis, which is the growth rate and color coded
%        zmax = eval(['ResultsAllCellLines.' samples{i} '.maxBiomass.f'])+ 0.001;
        zmax=eval(['ResultsAllCellLines.' samples{i} '.maxBiomass.f'])+ 0.001;;
        B = growthRates;
        % set values below tol as zero
        B(find(abs(B)<abs(tol)))=0;
        B(find(B==0))=nan; % remove zero as value

        % 3D plot
        figure1=figure;
        axes1 = axes('Parent',figure1);
        surf(ivalues',jvalues',B','Parent',axes1,'SpecularStrength',0.4,...
            'AmbientStrength',0.55,...
            'LineStyle','none');

        colorbar('peer',axes1);

        lowToHigh = [2 15];
        caxis(lowToHigh + [-1 1]*diff(lowToHigh));

        % set viewpoint
        az = 0;
        el = 90;
        view(az, el);

        %hold off
        hold on
        plot3([lb,lb],[0,-10],[zmax,zmax],'Color',[0 0 0],'LineWidth',1.2);

        hold on
        plot3([ub,ub],[-10,0],[zmax,zmax],'Color',[0 0 0],'LineWidth',1.2);

        hold off
        xlabel(label(1,1),'FontSize',fonts);
        ylabel(label(2,1),'FontSize',fonts);
        zlabel(label(3,1),'FontSize',fonts);
        title(title_name,'FontSize',fonts);
        axis tight %% seems to limit the axis, excludes white edge
              set(gcf, 'PaperUnits', 'inches');
        zlim([0 32]) % z-axis now from 0 to 10
        x_width=8 ;y_width=6;
        set(gcf, 'PaperPosition', [0 0 x_width y_width]); %
        name  = [ 'phasePlane_' samples{i} '_' strtok(mets{1,p}, '(') '_' strtok(mets{2,p}, '(')];
        saveas(gcf,[path filesep name '.pdf']);
        % saveas(gcf,[path filesep name '.png']);
        clear ivalues
        clear jvalues
    end

end
end

author: opencobra

repo: cobratoolbox


text:
% function checking if a specific file is a script%

function bool = isscript(fileName);

fid = fopen(fileName, 'r');
cont = true;
while cont && ~feof(fid)
    l = strtok(fgetl(fid));

    if ~isempty(l) && l(1) ~= '%'
        if strcmpi(l, 'function')
            bool = false; return;
        else
            bool = true; return;
        end;
    end;
end;
bool = true; return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function handles = enableSimGraph(hObject, handles)

set(handles.btnSimGraphSave, 'Enable', 'on');
set(handles.btnClusterData, 'Enable', 'on');

handles = togglePlotSimGraph(hObject, handles, 'on');

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function s = imMin(img, varargin)
%IMMIN Minimum value of a grayscale image, or of each color component
%
%   S = imMin(IMG)
%   Computes the minimum value of pixels in image IMG. If image is grayscale
%   image, the result is a scalar. If image is a color image, the result is
%   1-by-3 row vector, each componenent corresponding to one color of the
%   image.
%
%   S = imMin(IMG, MASK)
%   Computes the minimum value only in the area specified by MASK.
%
%   S = imMin(..., 'color', COL)
%   Forces the function to consider the image as color (if COL is TRUE) or
%   as grascale (if COL is FALSE). This can be useful for vector image with
%   more than 3 color components.
%
%
%   Example
%   % apply to cameraman image
%   img = imread('cameraman.tif');
%   imMin(img)
%   ans =
%       253
%
%   % apply to a RGB image
%   img = imread('peppers.png');
%   imMin(img)
%   ans =
%       255   255   255
%
%   See also
%   imMax, imMedian, imMean
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-30,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Process input arguments

% detect if image is color
color = isColorImage(img);

% check if user specified 'color' option
if length(varargin)>1
    var = varargin{end-1};
    if ischar(var)
        if strcmpi(var, 'color')
            color = varargin{end};
            varargin(end-1:end) = [];
        end
    end
end

%% Process color image

if color
    % If image is color, process each band separately

    % compute image size and dimension (including color dimension)
    dim = size(img);
    nd = length(dim);

    % create idnexing structure
    inds = cell(1, nd);
    for i=1:nd
        inds{i} = 1:dim(i);
    end

    % iterate on colors
    nc = dim(3);
    s = zeros(1, nc);
    for i=1:nc
        % modify the indexing structure to work on the i-th component
        inds{3} = i;
        s(i) = imMin(img(inds{:}), varargin{:});
    end

    return;
end

%% process grayscale image

if isempty(varargin)
    % compute min over all image
    s = min(img(:));
else
    % use first argument as mask
    s = min(img(varargin{1}));
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [medianPrice,medianTime,totalVol,nObs] = realized_compute_median(price,time,volume)
% Computes the median price at each time stamp for a vector of prices which may have multiple
% observations with the same time stamp.
%
% USAGE:
%   [MEDIANPRICE,MEDIANTIME,TOTALVOL,NOBS] = realized_compute_median(PRICE,TIME,VOLUME)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - m by 1 vector of times where TIME(i) corresponds to PRICE(i), sorted in
%                        ascending order
%   VOLUME           - [OPTIONAL] m by 1 vector of transaction volumes or bid/ask size.  If omitted
%                        TOTALVOLUME will be a vector of 0's.
%
% OUTPUTS:
%   MEDIANPRICE      - n by 1 vector of median prices where n is the number of unique elements of TIME
%   MEDIANTIME       - n by 1 vector of time stamps corresponding to MEDIANPRICE
%   TOTALVOLUME      - n by 1 vector capturing the total volume at each time stamp
%   NOBS             - n by 1 vector indicating number of prices at each element of MEDIANTIME
%
% COMMENTS:
%   This is a helper function for most of the Realized toolkit.  Most function in the Realized
%   toolkit expect time stamps to be unique.  Median price is a reasonable and fairly robust (to
%   noise) method of computing a unique price for each time stamp.

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 2    Date: 6/12/2011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% InputChecking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2 || nargin>3
    error('Six or seven inputs required.')
end
if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    error('PRICE must be a m by 1 vector.')
end
m = size(price,1);

if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    error('TIME must be sorted and increasing')
end
if size(time,2)>1 || length(time)~=m
    error('TIME must be a m by 1 vector.')
end
time = double(time);

if nargin == 3
    if size(volume,2)>size(volume,1)
        volume=volume';
    end
    if size(volume,2)>1 || length(volume)~=m
        error('VOLUME must be a m by 1 vector.')
    end
else
    volume = zeros(m,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Make sure time and volume are double since I store it was uint32
time = double(time);
volume = double(volume);

% Quickly find the time change points
pl = find(diff(time));
pl = [1;pl+1;length(time)+1];

% Get the number of observations
N = length(pl);

% Pre-allocatte
medianPrice = zeros(N-1,1);
medianTime = time(pl(1:N-1));
totalVol = zeros(N-1,1);
nObs = zeros(N-1,1);

% Loop over unique sizes, with special cases for 1 and 2
[sizes,~,ind]=unique(diff(pl));
for i=1:length(sizes)
    j=find(ind==i);
    loc = bsxfun(@plus,pl(j),0:sizes(i)-1)';
    switch sizes(i)
        case 1
            medianPrice(j) = price(pl(j));
            totalVol(j) = volume(pl(j));
        case 2
            medianPrice(j) = (price(pl(j)) + price(pl(j)+1))/2;
        otherwise
            sortedPrice = sort(price(loc));
            if mod(sizes(i),2)==0
                medianPrice(j) = mean(sortedPrice([sizes(i)/2;sizes(i)/2+1],:))';
            else
                medianPrice(j) = sortedPrice(ceil(sizes(i)/2),:)';
            end
    end
    totalVol(j) = sum(volume(loc))';
    nObs(j) = sizes(i);
end

% Old, slower method
% Loop
% singleton = pldiff==1;
% medianPrice(singleton) = price(pl(singleton));
% totalVol(singleton) = volume(pl(singleton));
% nObs(singleton) = 1;
% notSingleton = find(~singleton);
%
% for i=1:length(notSingleton)
%     j = notSingleton(i);
%     medianPrice(j) = median(price(pl(j):pl(j+1)-1));
%     totalVol(j) = sum(volume(pl(j):pl(j+1)-1));
%     nObs(j) = pl(j+1)-pl(j)+1;
% end

author: bashtage

repo: mfe-toolbox


text:
function FRR = CreditAssignment(SW,D)
% Credit assignment

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

    K = 4;  % Number of operators
    Reward = zeros(1,K);
    for i = 1 : K
        Reward(i) = sum(SW(2,SW(1,:)==i));
    end
    [~,Rank] = sort(Reward,'descend');
    [~,Rank] = sort(Rank);
    Decay    = D.^Rank.*Reward;
    FRR      = Decay./sum(Decay);
end

author: BIMK

repo: PlatEMO


text:
function model2xls(model, fileName, compSymbols, compNames)
% Writes a model to and Excel spreadsheet.
%
% USAGE:
%
%    model2xls(model, fileName, compSymbols, compNames)
%
% INPUT:
%    model:          A COBRA model struct
%    fileName:       filename with an xsl extension.
%
% OPTIONAL INPUT:
%
%    compSymbols:    Symbols of compartments used in metabolite ids
%    compNames:      Names of the compartments identified by the symbols
%
% EXAMPLE:
%
%                   'Reaction List' tab headers (case sensitive):
%
%                     * Required:
%
%                       * 'Abbreviation':      HEX1
%                       * 'Reaction':          `1 atp[c] + 1 glc-D[c] --> 1 adp[c] + 1 g6p[c] + 1 h[c]`
%                       * 'GPR':               (3098.3) or (80201.1) or (2645.3) or ...
%                     * Optional:
%
%                       * 'Description':       Hexokinase
%                       * 'Subsystem':         Glycolysis
%                       * 'Reversible':        0 (false) or 1 (true)
%                       * 'Lower bound':       0
%                       * 'Upper bound':       1000
%                       * 'Objective':         0/1
%                       * 'Confidence Score':  0,1,2,3,4
%                       * 'EC Number':         2.7.1.1;2.7.1.2
%                       * 'KEGG ID':           R000001
%                       * 'Notes':             Reaction also associated with EC 2.7.1.2
%                       * 'References':        PMID:2043117;PMID:7150652,...
%
%                   'Metabolite List' tab: Required headers (case sensitive): (needs to be complete list of metabolites,
%                   i.e., if a metabolite appears in multiple compartments it has to be represented in multiple rows.
%                   Abbreviations need to overlap with use in Reaction List
%
%                     * Required
%
%                       * 'Abbreviation':      glc-D or glc-D[c]
%                     * Optional:
%
%                       * 'Charged formula' or formula:   C6H12O6
%                       * 'Charge':                       0
%                       * 'Compartment':                  cytosol
%                       * 'Description':                  D-glucose
%                       * 'KEGG ID':                      C00031
%                       * 'PubChem ID':                   5793
%                       * 'ChEBI ID':                     4167
%                       * 'InChI string':                 InChI=1/C6H12O6/c7-1-2-3(8)4(9)5(10)6(11)12-2/h2-11H,1H2/t2-,3-,4+,5-,6?/m1/s1
%                       * 'SMILES':                       OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O
%                       * 'HMDB ID':                      HMDB00122
%
% NOTE:
%
%    Optional inputs may be required for input on unix machines.
%
% .. Authors:
%        - extracted from writeCbModel
%        - Thomas Pfau June 2017

if ~exist('compSymbols','var')
    if isfield(model,'comps') && isfield(model,'compNames')
        compSymbols = model.comps;
        compNames = model.compNames;
    else
        [compSymbols,compNames] = getDefaultCompartments();
    end
end

%Set up xlwrite in case its not set up.
setupxlwrite();

ReactionXlsFields = {'Abbreviation','Description','Reaction','GPR','Lower bound','Upper bound',...
    'Objective','Confidence Score','Subsystem','Notes','EC Number','References','KEGG ID'};
ReactionModelFields = {'rxns','rxnNames','formulas','grRules','lb','ub',...
    'c','rxnConfidenceScore','subSystems','rxnNotes','rxnECNumbers','rxnReferences','rxnKEGGID'};
if isfield(model,'rules')
    model = creategrRulesField(model);
end
if isfield(model,'osenseStr')
    %Only do something if minimisation is requested. Assume that osenseStr
    %is valid.
    if strcmpi(model.osenseStr,'min')
        model.c = model.c*-1;
    end
end

%Explicit handling of subSystems:
if isfield(model,'subSystems')
    cellBool = cellfun(@(x) iscell(x), model.subSystems);
    if all(cellBool)
        %Merge the subSystems with separating ';'
        model.subSystems = cellfun(@(x) strjoin(x,';'),model.subSystems,'UniformOutput',0);
    end
end

model.formulas = printRxnFormula(model,'printFlag',0);
actualFieldNames = fieldnames(model);
usedFields = ismember(ReactionModelFields,actualFieldNames);
ExcelFields = ReactionXlsFields(usedFields);
modelFields = ReactionModelFields(usedFields);
tmpData = cell(numel(model.rxns) + 1,sum(usedFields));
tmpData(1,:) = ExcelFields;
for i = 1:numel(ExcelFields)
    modelField = model.(modelFields{i});
    if isnumeric(modelField)

        for j = 2:length(modelField)+1
            tmpData{j,i} = num2str(modelField(j-1));
        end
        %tmpData(2:end,i) = mat2str(modelField);
    else
        modelField = cellfun(@chopForExcel, modelField,'UniformOutput',0);
        tmpData(2:end,i) = modelField;
    end
end

xlwrite(fileName,tmpData,'Reaction List');

MetaboliteXlsFields = {'Abbreviation','Description','Charged formula','Charge','Compartment','KEGG ID',...
    'PubChem ID','ChEBI ID','InChi string','SMILES','HMDB ID'};
MetaboliteModelFields = {'mets','metNames','metFormulas','metCharges','metCompartments','metKEGGID',...
    'metPubChemID','metChEBIID','metInChIString','metSmiles','metHMDBID'};
%detemine compartments
[tokens tmp_met_struct] = regexp(model.mets,'(?<met>.+)\[(?<comp>.+)\]$','tokens','names'); % add the third type for parsing the string such as "M_10fthf5glu_c"
%if we have any compartment, we will use unknown as compartment ID for
%metabolites without compartment.
if any(cellfun(@isempty, tmp_met_struct))
    unknownComp = 'u';
    unknownCompName = 'Unknown';
else
    unknownComp = 'c';
    unknownCompName = 'Cytosol';
end
if ~any(ismember(compSymbols,unknownComp))
    compSymbols{end+1} = unknownComp;
    compNames{end+1} = unknownCompName;
end

model.metCompartments = cell(size(model.mets));

for i = 1:numel(tmp_met_struct)
    if isempty(tmp_met_struct{i})
        model.mets{i} = [model.mets{i} '[' unknownComp ']'];
        model.metCompartments{i} = compNames(ismember(compSymbols,unknownComp));
    else
        model.metCompartments{i} = compNames(ismember(compSymbols,tmp_met_struct{i}.comp));
    end
end
actualFieldNames = fieldnames(model);
usedFields = ismember(MetaboliteModelFields,actualFieldNames);
ExcelFields = MetaboliteXlsFields(usedFields);
modelFields = MetaboliteModelFields(usedFields);
tmpData = cell(numel(model.mets) + 1,sum(usedFields));
tmpData(1,:) = ExcelFields;
for i = 1:numel(ExcelFields)
    modelField = model.(modelFields{i});
    if isnumeric(modelField)
        tmpData(2:end,i) = num2cell(modelField);
    else
        modelField = cellfun(@chopForExcel, modelField,'UniformOutput',0);
        tmpData(2:end,i) = modelField;
    end
end

xlwrite(fileName,tmpData,'Metabolite List');

end

%% Chop strings for excel output
function strOut = chopForExcel(str)

if (length(str) > 5000)
    strOut = str(1:5000);
    fprintf('String longer than 5000 characters - truncated for Excel output\n%s\n',str);
else
    strOut = str;
end
end


text:
function [SeqR1s,SeqR5s]=MADAll(fea_Train,gnd_Train,fea_Test,gnd_Test,testDims,wgts)
%Nearest neighbor classifier with MAD measure and symmetric matching

nTrain = size(fea_Train,1);
nTest = size(fea_Test,1);
nDim=length(testDims);%number of feature dimensions to test
SeqR1s=zeros(nDim,1);%Rank 1 recognition rate based on matching gait sequences
SeqR5s=zeros(nDim,1);%Rank 5 recognition rate based on matching gait sequences

for iDim=1:nDim
    Dim=testDims(iDim);
    feaTrn=fea_Train(:,1:Dim);
    feaTst=fea_Test(:,1:Dim);
    DMat=MAD(feaTrn,feaTst, wgts);%This is the MAD distance matrix
    [SeqR1s(iDim,1),SeqR5s(iDim,1)]=SeqDist(DMat,gnd_Train,gnd_Test);
end

%Calculate matching scores between each training/test sample pair
function D = MAD(feaTrn, feaTst, wgts)
if size(wgts,2)==1, wgts=wgts';end
[numTst,p] = size(feaTst);
numTrn = size(feaTrn,1);
wgts=wgts(1:p);
D = zeros(numTst,numTrn);
for i=1:numTst
    for j=1:numTrn
        A=feaTst(i,:);
        B=feaTrn(j,:);
        % the following two lines calculate the MAD between A and B
        dist=sum(sum(sum(A.*B./wgts)));
        D(i,j)=-dist/(norm(A(:))*norm(B(:)));
    end
end

%The recognition rate between two sequences
function [SeqR1,SeqR5]=SeqDist(DMat,gnd_Train,gnd_Test)
trnSeqs=unique(gnd_Train);
numTrnSeq=length(trnSeqs);

tstSeqs=unique(gnd_Test);
numTstSeq=length(tstSeqs);

SeqDMat=zeros(numTstSeq,numTrnSeq);
for i=1:numTstSeq
    idxs_i=find(gnd_Test==tstSeqs(i));
    iMat=DMat(idxs_i,:);
    for j=1:numTrnSeq
        idxs_j=find(gnd_Train==trnSeqs(j));
        jMat=iMat(:,idxs_j);
        SeqDMat(i,j)=mean(min(jMat,[],1))+mean(min(jMat,[],2));
    end
end
[minDs,minIdxs]=min(SeqDMat, [], 2);%Minimum distance implies best match
IDs=trnSeqs(minIdxs);
SeqR1=sum(IDs==tstSeqs)/numTstSeq;%Rank 1

[stDs,stIdxs]=sort(SeqDMat,2);
IDs=trnSeqs(stIdxs(:,1:5));
tstSeqs5=repmat(tstSeqs,1,5);
SeqR5=sum(sum(IDs==tstSeqs5))/numTstSeq;%Rank 5


text:
function value = cast_value_to_type(type,value,cls) %#ok<INUSD>
% Internal function for GUI parameter casting.

if isstruct(type)
    type = PropertyType(type.PrimitiveType,type.Shape,type.Domain); end

validateattributes(type, {'PropertyType'}, {'scalar'});
assert(type.CanAccept(value), ...
    'PropertyGridField:InvalidOperation', ...
    'Setting type "%s" would invalidate current property value.', char(type));
value = type.ConvertFromMatLab(value);

repo: REST


text:
function [w3j,jmin, jmax] = Wigner3j_new(j2, j3, m1, m2, m3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	This subroutine will calculate the Wigner 3j symbols
%
%		j  j2 j3
%		m1 m2 m3
%
%	for all allowable values of j. The returned values in the array j are
%	calculated only for the limits
%
%		jmin = max(|j2-j3|, |m1|)
%		jmax = j2 + j3
%
%	To be non-zero, m1 + m2 + m3 = 0. In addition, it is assumed that all j and m are
%	integers. Returned values have a relative error less than ~1.d-8 when j2 and j3
%	are less than 103 (see below). In practice, this routine is probably usable up to 165.
%
%	This routine is based upon the stable non-linear recurence relations of Luscombe and
%	Luban (1998) for the "non classical" regions near jmin and jmax. For the classical
%	region, the standard three term recursion relationship is used (Schulten and Gordon 1975).
%	Note that this three term recursion can be unstable and can also lead to overflows. Thus
%	the values are rescaled by a factor "scalef" whenever the absolute value of the 3j coefficient
%	becomes greater than unity. Also, the direction of the iteration starts from low values of j
%	to high values, but when abs(w3j(j+2)/w3j(j)) is less than one, the iteration will restart
%	from high to low values. More efficient algorithms might be found for specific cases
%	(for instance, when all m's are zero).
%
%	Verification:
%
%	The results have been verified against this routine run in quadruple precision.
%	For 1.e7 acceptable random values of j2, j3, m2, and m3 between -200 and 200, the relative error
%	was calculated only for those 3j coefficients that had an absolute value greater than
%	1.d-17 (values smaller than this are for all practical purposed zero, and can be heavily
%	affected by machine roundoff errors or underflow). 853 combinations of parameters were found
%	to have relative errors greater than 1.d-8. Here I list the minimum value of max(j2,j3) for
%	different ranges of error, as well as the number of times this occured
%
%	1.d-7 < error  <=1.d-8 = 103	# = 483
%	1.d-6 < error <= 1.d-7 =  116	# = 240
%	1.d-5 < error <= 1.d-6 =  165	# = 93
%	1.d-4 < error <= 1.d-5 = 167	# = 36
%
%	Many times (maybe always), the large relative errors occur when the 3j coefficient
%	changes sign and is close to zero. (I.e., adjacent values are about 10.e7 times greater
%	in magnitude.) Thus, if one does not need to know highly accurate values of the 3j coefficients
%	when they are almost zero (i.e., ~1.d-10)  this routine is probably usable up to about 160.
%
%	These results have also been verified for parameter values less than 100 using a code
%	based on the algorith of de Blanc (1987), which was originally coded by Olav van Genabeek,
%	and modified by M. Fang (note that this code was run in quadruple precision, and
%	only calculates one coefficient for each call. I also have no idea if this code
%	was verified.) Maximum relative errors in this case were less than 1.d-8 for a large number
%	of values (again, only 3j coefficients greater than 1.d-17 were considered here).
%
%	The biggest improvement that could be made in this routine is to determine when one should
%	stop iterating in the forward direction, and start iterating from high to low values.
%
%	Calling parameters
%		IN
%			j2, j3, m1, m2, m3 	Integer values.
%		OUT
%			w3j			Array of length jmax - jmin + 1.
%			jmin, jmax		Minimum and maximum values
%						out output array.
%	Dependencies: None
%
%	Written by Mark Wieczorek August (2004)
%
%	August 2009: Based on the suggestions of Roelof Rietbroek, the calculation of RS has been slightly
%	modified so that division by zero will not cause a run time crash (this behavior depends on how the
%	compiler treats IEEE floating point exceptions). These values were never used in the original code
%	when this did occur.
%
%	Copyright (c) 2005-2009, Mark A. Wieczorek
%	All rights reserved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	implicit none
%	integer, intent(in) ::	j2, j3, m1, m2, m3
%	integer, intent(out) ::	jmin, jmax
%	real*8, intent(out) ::	w3j(:)
%	real*8 ::		wnmid, wpmid, scalef, denom, rs(j2+j3+1), &
%				wl(j2+j3+1), wu(j2+j3+1), xjmin, yjmin, yjmax, zjmax, xj, zj
%	integer :: 		j, jnum, jp, jn, k, flag1, flag2, jmid

%% some basic constants

flag1 = 0;
flag2 = 0;

scalef = 1;

jmin = max(abs(j2-j3), abs(m1));
jmax = j2 + j3;
jnum = jmax - jmin + 1;

w3j = zeros(1,jnum);

%% some simple functions
jindex = @(j) j-jmin+1;

a = @(j) sqrt((j.^2 - (j2-j3).^2) .* ((j2+j3+1).^2 - j.^2) * (j.^2-m1^2));

y= @(j) -(2*j+1) .* ( m1 .* (j2.*(j2+1) - j3*(j3+1)) - (m3-m2).*j.*(j+1));

x = @(j) j .* a(j+1);

z = @(j) (j+1) .* a(j);

%% exclude some basic cases	which give zero
if abs(m2) > j2 || abs(m3) > j3
  return
elseif m1 + m2 + m3 ~= 0
  return
elseif jmax < jmin
  return
end

%% Only one term is present

if jnum == 1

  w3j = 1 / sqrt(2 * jmin + 1);

  if (w3j < 0 && (-1)^(j2-j3+m2+m3) > 0) || ...
      (w3j > 0 && (-1)^(j2-j3+m2+m3) < 0)
    w3j = -w3j;
  end
  return

end

%% more then one term
%
% Calculate lower non-classical values for [jmin, jn]. If the second term
%	can not be calculated because the recursion relationsips give rise to a
%	1/0,  set flag1 to 1.  If all m's are zero,  this is not a problem
%	as all odd terms must be zero.
%

rs = 0;
wl = 0;

xjmin = x(jmin);
yjmin = y(jmin);

if (m1 == 0 && m2 == 0 && m3 == 0) % All m's are zero

  wl(jindex(jmin)) = 1;
  wl(jindex(jmin+1)) = 0;
  jn = jmin+1;

elseif yjmin == 0 % The second terms is either zero

  if xjmin == 0   % or undefined
    flag1 = 1;
    jn = jmin;
  else
    wl(jindex(jmin)) = 1;
    wl(jindex(jmin+1)) = 0;
    jn = jmin+1;
  end

elseif xjmin * yjmin >= 0 % The second term is outside of the non-classical region

  wl(jindex(jmin)) = 1;
  wl(jindex(jmin+1)) = -yjmin / xjmin;
  jn = jmin+1;

else							% Calculate terms in the non-classical region

  rs(jindex(jmin)) = -xjmin / yjmin;

  jn = jmax;
  for j = jmin + 1:jmax-1

    denom =  y(j) + z(j)*rs(jindex(j-1));
    xj = x(j);
    if abs(xj) > abs(denom) || xj * denom >= 0 || denom == 0
      jn = j-1;
      break
    else
      rs(jindex(j)) = -xj / denom;
    end

  end

  wl(jindex(jn)) = 1;

  for k = 1:jn - jmin

    wl(jindex(jn-k)) = wl(jindex(jn-k+1)) * rs(jindex(jn-k));

  end

  if jn == jmin	% Calculate at least two terms so that

    wl(jindex(jmin+1)) = -yjmin / xjmin;		% these can be used in three term
    jn = jmin+1;					% recursion

  end

end

if jn == jmax					% All terms are calculated

  w3j = wl;

  % normalize
  norm = sum((2*(jmin:jmax)+1) .* w3j.^2);
  w3j = w3j ./ sqrt(norm);

  % fix sign
  if (w3j(end) < 0 && (-1)^(j2-j3+m2+m3) > 0) || ...
      (w3j(end) > 0 && (-1)^(j2-j3+m2+m3) < 0)
    w3j = -w3j;
  end

  return
end

%%
%
% 	Calculate upper non-classical values for [jp, jmax].
%	If the second last term can not be calculated because the
%	recursion relations give a 1/0,  set flag2 to 1.
%	(Note, I don't think that this ever happens).
%

wu = 0;

yjmax = y(jmax);
zjmax = z(jmax);

if (m1 == 0 && m2 == 0 && m3 == 0)

  wu(jindex(jmax)) = 1;
  wu(jindex(jmax-1)) = 0;
  jp = jmax-1;

elseif yjmax == 0

  if zjmax == 0
    flag2 = 1;
    jp = jmax;
  else
    wu(jindex(jmax)) = 1;
    wu(jindex(jmax-1)) = - yjmax / zjmax;
    jp = jmax-1;
  end

elseif yjmax * zjmax >= 0

  wu(jindex(jmax)) = 1;
  wu(jindex(jmax-1)) = - yjmax / zjmax;
  jp = jmax-1;

else
  rs(jindex(jmax)) = -zjmax / yjmax;

  jp = jmin;
  for j = jmax-1:-1:jn

    denom = y(j) + x(j)*rs(jindex(j+1));
    zj = z(j);
    if abs(zj) > abs(denom) || zj * denom >= 0 || denom == 0
      jp = j+1;
      break
    else
      rs(jindex(j)) = -zj / denom;
    end

  end

  wu(jindex(jp)) = 1;

  for k=1:jmax - jp
    wu(jindex(jp+k)) = wu(jindex(jp+k-1))*rs(jindex(jp+k));
  end

  if jp == jmax
    wu(jindex(jmax-1)) = - yjmax / zjmax;
    jp = jmax-1;
  end

end

%%
%
% 	Calculate classical terms for [jn+1, jp-1] using standard three
% 	term rercusion relationship. Start from both jn and jp and stop at the
% 	midpoint. If flag1 is set,  perform the recursion solely from high to
% 	low values. If flag2 is set,  perform the recursion solely from low to high.

if flag1 == 0

  jmid = ceil((jn + jp)/2);

  for j = jn : jmid - 1

    wl(jindex(j+1)) = - (z(j)*wl(jindex(j-1)) +y(j)*wl(jindex(j))) / x(j);

    if abs(wl(jindex(j+1))) > 1 				% watch out for overflows.
      wl(jindex(jmin):jindex(j+1)) = wl(jindex(jmin):jindex(j+1)) / scalef;
    end

    % if values are decreasing
    if (abs(wl(jindex(j+1)) / wl(jindex(j-1))) < 1 && wl(jindex(j+1)) ~= 0)

      %  stop upward iteration
      jmid = j+1;	% and start with the downward
      break	% iteration.
    end
  end

  wnmid = wl(jindex(jmid));

  if (abs(wnmid/wl(jindex(jmid-1))) < 1e-6 && ...
      wl(jindex(jmid-1)) ~= 0) 				% Make sure that the stopping

    wnmid = wl(jindex(jmid-1));					% midpoint value is not a zero,
    jmid = jmid - 1;							% or close to it%
  end

  for j=jp:-1:jmid+1
    wu(jindex(j-1)) = - (x(j)*wu(jindex(j+1)) + y(j)*wu(jindex(j)) ) / z(j);
    if (abs(wu(jindex(j-1))) > 1)
      wu(jindex(j-1):jindex(jmax)) = wu(jindex(j-1):jindex(jmax)) / scalef;
    end

  end

  wpmid = wu(jindex(jmid));

  % rescale two sequences to common midpoint

  if jmid == jmax
    w3j(1:jnum) = wl(1:jnum);
  elseif jmid == jmin
    w3j(1:jnum) = wu(1:jnum);
  else
    w3j(1:jindex(jmid)) = wl(1:jindex(jmid)) * wpmid / wnmid;
    w3j(jindex(jmid+1):jindex(jmax)) = wu(jindex(jmid+1):jindex(jmax));
  end

elseif (flag1 == 1 && flag2 == 0) 	% iterature in downward direction only

  for j=jp:-1:jmin+1
    wu(jindex(j-1)) = - (x(j)*wu(jindex(j+1)) + y(j)*wu(jindex(j)) ) / z(j);
    if (abs(wu(jindex(j-1))) > 1)
      wu(jindex(j-1):jindex(jmax)) = wu(jindex(j-1):jindex(jmax)) / scalef;
    end
  end

  w3j(1:jnum) = wu(1:jnum);

elseif flag2 == 1 && flag1 == 0 % iterature in upward direction only

  for j = jn:jp-1
    wl(jindex(j+1)) = - (z(j)*wl(jindex(j-1)) +y(j)*wl(jindex(j))) / x(j);
    if abs(wl(jindex(j+1))) > 1
      wl(jindex(jmin):jindex(j+1)) = wl(jindex(jmin):jindex(j+1))/ scalef;
    end
  end

  w3j = wl;

elseif flag1 == 1 && flag2 == 1

  error('Can not calculate function for input values');

end

% normalize
norm = sum((2*(jmin:jmax)+1) .* w3j.^2);
w3j = w3j ./ sqrt(norm);

% fix sign
if (w3j(end) < 0 && (-1)^(j2-j3+m2+m3) > 0) || ...
    (w3j(end) > 0 && (-1)^(j2-j3+m2+m3) < 0)
  w3j = -w3j;
end

end

author: mtex-toolbox

repo: mtex

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


text:
report_this_filefun(mfilename('fullpath'));

figure
rect = [0.15,  0.60, 0.75, 0.35];
axes('position',rect)

pl =plot(bg(:,10),-bg(:,7),'gx');
set(pl,'Linewidth',1.5,'MarkerSize',7)
hold on

pl =plot(aft(:,12),-aft(:,7),'bx');
set(pl,'LineWidth',1.5,'MarkerSize',7)

pl =plot(maex,-maey,'hm');
set(pl,'LineWidth',1.5,'MarkerSize',15,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')

set(gca,'Ylim',[-14 0])
set(gca,'Xlim',[ 0 48 ])

set(gca,'Color',[color_bg([1,2]) 0.7])
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',16,'Linewidth',1.2)
ylabel('Depth in [km]')

set(gca,'XTickLabel',[]);

ax = axis;

rect = [0.15,  0.15, 0.75, 0.35];
axes('position',rect)

hold on
pco1 = pcolor(gx,gy,re4);
hold on

pl =plot(maex,-maey,'hm');
set(pl,'LineWidth',1.5,'MarkerSize',15,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')

colormap(hsv)

shading interp

if fre == 1
    caxis([fix1 fix2])
end

axis([ax]);
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',16,'Linewidth',1.8)

ylabel('Depth in  [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
xlabel('Distance along projection in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

vx =  (min(min(re3)):0.1:max(max(re3)));
if fre == 1; vx =  (fix1:0.1:fix2); end
v = [vx ; vx]; v = v';
rect = [0.94 0.15 0.01 0.35];
axes('position',rect)
pcolor((1:2),vx,v)
shading interp
set(gca,'XTickLabels',[])
set(gca,'FontSize',12,'FontWeight','bold',...
    'LineWidth',1.0,'YAxisLocation','right',...
    'Box','on','SortMethod','childorder','TickDir','out')
ax3 = gca;

whitebg(gcf,[0 0 0])
set(gcf,'Color','k','InvertHardcopy','off')

matdraw

author: CelsoReyes


text:
%GAUSSM Trainable mapping, mixture of Gaussians (MoG) density estimate
%
%  W = GAUSSM(A,K,R,S,M)
%  W = A*GAUSSM([],K,R,S,M)
%  W = A*GAUSSM(K,R,S,M)
%
% INPUT
%   A     Dataset
%   K     Number of Gaussians per class
%   R,S,M Regularization parameters, 0 <= R,S <= 1, see QDC
%
% OUTPUT
%   W     Mixture of Gaussians density estimate
%
% DESCRIPTION
% Estimation of a PDF by the dataset A by a mixture of Gaussians procedure.
% Use is made of EMCLUST(A,QDC,K). Unlabeled objects are neglected, unless
% A is entirely unlabeled or double. Then all objects are used. If A is a
% multi-class crisp labeled dataset the densities are estimated class by
% class and then weighted and combined according their prior probabilities.
% Use +A instead of A to obtain a single set of Gaussians. In all cases,
% just single density estimator W is returned.
%
% Note that it is necessary to set the label type of A to soft labels
% (A = LABTYPE(A,'soft') in order to use the traditional EM algorithm
% based on posterior probabilities instead of using crisp labels.
%
% The mapping W may be applied to a new dataset B using DENSITY = B*W.
%
%  W = A*GAUSSM
%
% uses a single Gaussian per class (K=1) and no regularisation. If
% regulariisation is desired, also K should be supplied.
%
% EXAMPLE
% a = gendatb;
% w = a*gaussm(2);
% scatterd(a)
% plotm(w)
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, QDC, MOGC, EMCLUST, PLOTM, TESTC

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: gaussm.m,v 1.8 2009/02/02 21:57:39 duin Exp $

function w = gaussm(varargin)

  mapname = 'Mixture of Gaussians';
	argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],1,0,0,[]);
  if mapping_task(argin,'definition')
    w = define_mapping(argin,'untrained');
    w = setname(w,mapname);
  else
    [a,n,r,s,dim] = deal(argin{:});

    if isa(a,'prdataset')
      labname = getname(a);
    else
      labname = '';
    end

    if ((~isdataset(a) && ~isdatafile(a)) || ...
        (getsize(a,3) ~= 1 && islabtype(a,'crisp')))
      % this handles the multiclass situation
      w = mclassm(a,prmapping(mfilename,n),'weight');
      w = setlabels(w,labname);
      w = setname(w,mapname);

    else
      % here we have a single class
      [m,k] = getsize(a);
      if n == 1
        % just one Gaussian, esitmate it
        [U,G] = meancov(a);
        res.mean = +U;
        res.cov  = G;
        res.prior= 1;
        w = normal_map(res,labname,k,1);
      else
        % multiple Gaussians, run EM
        [e,v] = emclust(a,qdc([],r,s,dim),n);
        ncomp0 = size(v.data.mean,1);
        iter = 0;
        while (ncomp0 ~= n & iter < 5)   % repeat until exactly n components are found
          [e,v1] = emclust(a,qdc([],r,s,m),n);
          ncomp1 = size(v1.data.mean,1);
          if ncomp1 > ncomp0
            v = v1;
            ncomp0 = ncomp1;
          end
          iter = iter + 1;
        end
        res = v.data;
        res.nlab = ones(n,1); % defines that all Gaussian components have to be
                              % combined into a single class.
        w = prmapping('normal_map','trained',res,labname,k,1);
      end
      w = setname(w,mapname);

    end

  end

return

author: marianux

repo: ecg-kit

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


text:
%% Misorientation Distribution Function
% Explains how to compute and analyze misorientation distribution
% functions.
%
%% TODO
% Please help to redo the section
%
%%
% When speaking about the misorientation distribution function (MDF) one
% has to differentiate to cases
%
% # the boundary (correlated) misorientation distribution function
% # the uncorelated misorientation distribution function
%
% While the first one considers only misorientations at grain boundaries
% the second one considers misorietation between arbitrary crystal
% orientations. To illustrate the difference lets consider the following
% EBSD data set

% Lets import some EBSD data and reconstruct the grains.

mtexdata forsterite
grains = calcGrains(ebsd)

%% The boundary misorientation distribution function
%
% In order to compute the boundary misorientation distribution function for
% the phase transition from Forsterite to Enstatite we first extract the
% misorientations along all Forsterite to Enstatite boundary segements

mori_boundary = grains.boundary('Fo','En').misorientation

%%
% and second compute the corresponding density function using the command
% <orientation.calcDensity.html calcDensity>

mdf_boundary = calcDensity(mori_boundary,'halfwidth',5*degree)

%%

adf_boundary = mdf_boundary.calcAxisDistribution

plot(adf_boundary)

%%

%%
% The misorientation distribution function can be processed as any other
% ODF. E.g. we can compute the prefered misorientation via

[v,mori] = max(mdf_boundary)

%%
% or plot the pole figure corresponding to the crystal axis (1,0,0)

plotPDF(mdf_boundary,Miller(1,0,0,ebsd('Fo').CS))

%% The uncorrelated misorientation distribution function
%
% Alternatively the uncorrelated misorientation distribution function can be
% computed by providing the option *uncorrelated*

mori = calcMisorientation(ebsd('En'),ebsd('Fo'))
mdf_uncor = calcDensity(mori)

%%
% Obviously it is different from the boundary misorientation distribution
% function.

plotPDF(mdf_uncor,Miller(1,0,0,ebsd('Fo').CS))

%% Computing the uncorrelated misorientation function from two ODFs
%
% Let given two odfs

odf_fo = calcDensity(ebsd('fo').orientations,'halfwidth',10*degree)
odf_en = calcDensity(ebsd('en').orientations,'halfwidth',10*degree)

%%
% Then the uncorrelated misorientation function between these two ODFs can
% be computed by

mdf = calcMDF(odf_en,odf_fo)

%%
% This misorientation distribution function should be similar to the
% uncorrelated misorientation function computed directly from the ebsd data

plotPDF(mdf,Miller(1,0,0,ebsd('Fo').CS))

%% Analyzing misorientation functions
%
%

%% SUB: Angle distribution
%
% Let us first compare the actual angle distribution of the boundary
% misorientations with the theoretical angle distribution of the
% uncorrelated MDF.

close all
plotAngleDistribution(grains.boundary('fo','en').misorientation)

hold on

plotAngleDistribution(mdf)

hold off

%%
% For computing the exact values see the commands
% <SO3Fun.calcAngleDistribution.html calcAngleDistribution(mdf)> and
% <orientation.calcAngleDistribution.html calcAngleDistribution(ori)>.

%% SUB: Axis distribution
%
% The same we can do with the axis distribution. First the actual angle distribution of the boundary
% misorientations

plotAxisDistribution(grains.boundary('fo','en').misorientation,'smooth')

%%
% Now the theoretical axis distribution of the
% uncorrelated MDF.

plotAxisDistribution(mdf)

%%
% For computing the exact values see the commands
% <SO3Fun.calcAxisDistribution.html calcAxisDistribution(mdf)> and
% <orientation.calcAxisDistribution.html calcAxisDistribution(grains)>.

aD = calcDensity(axis(grains.boundary('fo','en').misorientation))


text:
function ls=v_lpcar2ls(ar)
%V_LPCAR2LS convert ar polynomial to line spectrum pair frequencies LS=(AR)
% output vector elements will be in range 0 to 0.5
% the returned vector will be of length p

% This routine is nowhere near as efficient as it might be

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_lpcar2ls.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(ar);
p = p1-1;
p2 = fix(p/2);
d=0.5/pi;

if rem(p,2)		% odd order
  for k=1:nf
    aa=[ar(k,:) 0];
    r = aa + fliplr(aa);
    q = aa - fliplr(aa);
    fr = sort(angle(roots(r)));
    fq = [sort(angle(roots(deconv(q,[1 0 -1])))); 0];
    f = [fr(p2+2:p+1).' ; fq(p2+1:p).'];
    f(p+1) = [];
    ls(k,:) = d*f(:).';
  end
else
  for k=1:nf
    aa=[ar(k,:) 0];
    r = aa + fliplr(aa);
    q = aa - fliplr(aa);
    fr = sort(angle(roots(deconv(r,[1 1]))));
    fq = sort(angle(roots(deconv(q,[1 -1]))));
    f = [fr(p2+1:p).' ; fq(p2+1:p).'];
    ls(k,:) = d*f(:).';
  end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% //This software is licensed under the BSD 3 Clause license: http://opensource.org/licenses/BSD-3-Clause
%
%
% //Copyright (c) 2013, University of Oxford
% //All rights reserved.
%
% //Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%
% //Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% //Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% //Neither the name of the University of Oxford nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
% //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%   The method implemented in this file has been patented by their original
%   authors. Commercial use of this code is thus strongly not
%   recocomended.
%
% //Authors: 	Gari D Clifford -
% //            Roberta Colloca -
% //			Julien Oster	-

function [BC,PC,sZ] = BPcount( sZ )

% BPcount counts the number of non empty bins (BC) in the matrix Z and
% the number of {dRR(i),dRR(i-1)} in the same region (PC) and deletes
% the counted points from the matrix Z

%bdc is the BIN diagonal count: number of non empty bins contained in
%the i-th diagonal of Z
bdc=0;
BC=0;
%pdc is the POINTS diagonal count: number of {dRR(i),dRR(i-1)} contained in
%the i-th diagonal of Z
pdc=0;
PC=0;

for i=-2:2
     bdc=sum(diag(sZ,i)~=0);
     pdc=sum(diag(sZ,i));
     BC=BC+bdc;
     PC=PC+pdc;
     sZ=sZ-diag(diag(sZ,i),i);
 end

end

author: fernandoandreotti


text:
function test_bug1870

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_megplanar ft_datatype_sens ft_compute_leadfield

% this test is basically a small section of test_ft_megplanar
% the input data is consistent, but along the way the grad structure gets screwed up

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1870.mat'));
dataP = ft_megplanar(cfg, data);

% although it was not reported, also check on the combining of the planar channels
cfg = [];
dataC = ft_combineplanar(cfg, dataP);

% do some sanity checks
assert(isfield(dataP.grad, 'tra'));
ft_datatype_sens(dataP.grad);
ft_datatype_sens(dataC.grad);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip